程序员经验分享-hwhale

鲁棒优化（4）：通过yalmip中的kkt命令实现CCG两阶段鲁棒优化

2023-11-15

两阶段鲁棒优化的原理推导部分，已经较多的文章进行分析。目前大部分同学面临的问题是，子问题模型中存在的双线性项该如何处理？

目前，主流方式是，采用对偶定理或KKT条件，将第二阶段的双层问题变成单层问题。
简略的思想如下：
首先是原始的两阶段模型：
在这里插入图片描述
对上述的两阶段模型，展开分成主问题与子问题：
在这里插入图片描述
主问题与子问题相互迭代，当两个问题的最优解不断收敛并相等时，两阶段鲁棒CCG问题求解完成。

更具体原理推导过程详见：
鲁棒优化| C&CG算法求解两阶段鲁棒优化：全网最完整、最详细的【入门-完整推导-代码实现】笔记
微电网两阶段鲁棒优化经济调度方法
列与约束生成(Column and Constraint Generation, C&CG/CCG)算法

以微电网经济调度为例，编写如下案例。程序中yalmip KKT命令的使用方法详见yalmip官网https://yalmip.github.io/command/kkt/

主体程序如下：

LB_recoder = -3000;
UB_recoder = +3000;

for i=1:6
    if i==1
        Load_u1 = [88.24 	83.01 	80.15 	79.01 	76.07 	78.39 	89.95 	128.85 	155.45 	176.35 	193.71 	182.57 	179.64 	166.31 	164.61 	164.61 	174.48 	203.93 	218.99 	238.11 	216.14 	173.87 	131.07 	94.04];
        [LB, Temp_net, Temp_cha, Temp_dis,X_1] = MP1(Load_u1); %给定初始场景Load，获得下界，以及01变量。
    else
        [LB, Temp_net, Temp_cha, Temp_dis] = MP(X); %给定初始场景Load，获得下界，以及01变量
    end
    LB_recoder = [LB_recoder, LB];
    [UB, X] = SP(Temp_net, Temp_cha, Temp_dis); %传入01变量，获得最坏的场景Load
    UB_recoder = [UB_recoder, UB];
    if UB-LB<=3
        break;
    end
end

plot(LB_recoder); %画图
hold on
plot(UB_recoder);

主问题如下：

function [LB, Temp_net, Temp_cha, Temp_dis] = MP(X) %传入子问题产生的割集
Load_u = X(1, :);
Pbuy_1 = X(2, :);
Psell_1 = X(3, :);
Pcha_1 = X(4, :);
Pdis_1 = X(5, :);

%-------------------------常量定义-----------------------%
%风机预测出力
Pw = [66.9	68.2	71.9	72	78.8	94.8	114.3	145.1	155.5	142.1	115.9	127.1	141.8	145.6...
    145.3	150	206.9	225.5	236.1	210.8	198.6	177.9	147.2	58.7];
%光伏预测出力
Ppv = [0	0	0	0	0.06	6.54	20.19	39.61	49.64	88.62	101.59	66.78	110.46	67.41	31.53...
    50.76	20.6	22.08	2.07	0	0	0	0	0];
%分时电价
C_buy = [0.25	0.25	0.25	0.25	0.25	0.25	0.25	0.53	0.53	0.53	0.82	0.82...
    0.82	0.82	0.82	0.53	0.53	0.53	0.82	0.82	0.82	0.53	0.53	0.53];
C_sell = [0.22	0.22	0.22	0.22	0.22	0.22	0.22	0.42	0.42	0.42	0.65	0.65...
    0.65	0.65	0.65	0.42	0.42	0.42	0.65	0.65	0.65	0.42	0.42	0.42];
%% 各变量及常量定义
%------------------------变量定义-----------------------%
Pbuy = sdpvar(1, 24, 'full'); %从电网购电电量
Psell = sdpvar(1, 24, 'full'); %向电网售电电量
Pcha = sdpvar(1, 24);
Pdis = sdpvar(1, 24);
Temp_net = binvar(1, 24, 'full'); % 购|售电标志
Temp_cha = binvar(1, 24, 'full'); %充电标志
Temp_dis = binvar(1, 24, 'full'); %放电标志
a = sdpvar(1, 1);
st = [Pbuy - Psell + Pw + Ppv == Load_u + Pcha - Pdis, ... %主网功率交换约束
    0 <= Pbuy <= Temp_net .* 160, ...
    0 <= Psell <= (1 - Temp_net) .* 160, ...
    0 <= Pcha <= Temp_cha .* 40, ...
    0 <= Pdis <= Temp_dis .* 40, ...
    Temp_cha + Temp_dis <= 1, ...
    sum(Pcha - Pdis) == 0, ...
    C_buy * Pbuy' - C_sell * Psell' - 0.2 * ones(1, 24) * Pdis' <= a];
for t = 1 : 24
    st = [st, -30 <= sum(Pcha(1, 1 : t) - Pdis(1, 1 : t)) <= 165]; %SOC约束，电池容量300kwh，初始S0C为0.4，0.3<=SOC<=0.95
end

st = [st, C_buy * Pbuy_1' - C_sell * Psell_1' - 0.2 * ones(1, 24) * Pdis_1' <= a]; %需要更新的约束

%% 目标函数
obj = a;
ops = sdpsettings('solver', 'gurobi'); %参数指定程序用gurobi求解器
optimize(st, obj, ops)
Temp_net = value(Temp_net);
Temp_cha = value(Temp_cha);
Temp_dis = value(Temp_dis);
LB = value(obj);
end


function [LB, Temp_net, Temp_cha, Temp_dis, X_1] = MP1(Load_u)
%-------------------------常量定义-----------------------%
% Load_u=[88.24 	83.01 	80.15 	79.01 	76.07 	78.39 	89.95 	128.85 	155.45 	176.35 	193.71 	182.57 	179.64 	166.31 	164.61 	164.61 	174.48 	203.93 	218.99 	238.11 	216.14 	173.87 	131.07 	94.04];
%风机预测出力
Pw = [66.9	68.2	71.9	72	78.8	94.8	114.3	145.1	155.5	142.1	115.9	127.1	141.8	145.6...
    145.3	150	206.9	225.5	236.1	210.8	198.6	177.9	147.2	58.7];
%光伏预测出力
Ppv = [0	0	0	0	0.06	6.54	20.19	39.61	49.64	88.62	101.59	66.78	110.46	67.41	31.53...
    50.76	20.6	22.08	2.07	0	0	0	0	0];
%分时电价
C_buy = [0.25	0.25	0.25	0.25	0.25	0.25	0.25	0.53	0.53	0.53	0.82	0.82...
    0.82	0.82	0.82	0.53	0.53	0.53	0.82	0.82	0.82	0.53	0.53	0.53];
C_sell = [0.22	0.22	0.22	0.22	0.22	0.22	0.22	0.42	0.42	0.42	0.65	0.65...
    0.65	0.65	0.65	0.42	0.42	0.42	0.65	0.65	0.65	0.42	0.42	0.42];
%% 各变量及常量定义
%------------------------变量定义-----------------------%
Pbuy = sdpvar(1, 24, 'full'); %从电网购电电量
Psell = sdpvar(1, 24, 'full'); %向电网售电电量
Temp_net = binvar(1, 24, 'full'); % 购|售电标志
Temp_cha = binvar(1, 24, 'full'); %充电标志
Temp_dis = binvar(1, 24, 'full'); %放电标志
Pcha = sdpvar(1, 24);
Pdis = sdpvar(1, 24);
a = sdpvar(1, 1);
st = [Pbuy - Psell + Pw + Ppv == Load_u + Pcha - Pdis, ...
    0 <= Pbuy <= Temp_net .* 160, ...
    0 <= Psell <= (1 - Temp_net) .* 160, ...
    0 <= Pcha <= Temp_cha .* 40, ...
    0 <= Pdis <= Temp_dis .* 40, ...
    Temp_cha + Temp_dis <= 1, ...
    sum(Pcha - Pdis) == 0, ...
    C_buy * Pbuy' - C_sell * Psell' - 0.2 * ones(1, 24) * Pdis' <= a]; %主网功率交换约束,...
for t = 1 : 24
    st = [st, -30 <= sum(Pcha(1, 1 : t) - Pdis(1, 1 : t)) <= 165]; %SOC约束，电池容量300kwh，初始S0C为0.4，0.3<=SOC<=0.95
end
%% 目标函数
obj = a;
ops = sdpsettings('solver', 'gurobi'); %参数指定程序用cplex求解器
optimize(st, obj, ops)
Temp_net = value(Temp_net);
Temp_cha = value(Temp_cha);
Temp_dis = value(Temp_dis);
LB = value(obj);

Pbuy = value(Pbuy); %从电网购电电量
Psell = value(Psell); %向电网售电电量
Pcha = value(Pcha);
Pdis = value(Pdis);
X_1 = [Load_u; Pbuy; Psell; Pcha; Pdis];
end

子问题如下：

function [UB, X] = SP(Temp_net, Temp_cha, Temp_dis)
%%目标函数值中不能包含风光出力收益，则原问题和对偶问题都等价
%% 各变量及常量定义
Load = [88.24 	83.01 	80.15 	79.01 	76.07 	78.39 	89.95 	128.85 	155.45 	176.35 	193.71 	182.57 	179.64 	166.31 	164.61 	164.61 	174.48 	203.93 	218.99 	238.11 	216.14 	173.87 	131.07 	94.04];
%风机预测出力
Pw = [66.9	68.2	71.9	72	78.8	94.8	114.3	145.1	155.5	142.1	115.9	127.1	141.8	145.6...
    145.3	150	206.9	225.5	236.1	210.8	198.6	177.9	147.2	58.7];
%光伏预测出力
Ppv = [0	0	0	0	0.06	6.54	20.19	39.61	49.64	88.62	101.59	66.78	110.46	67.41	31.53...
    50.76	20.6	22.08	2.07	0	0	0	0	0];
%分时电价
C_buy = [0.25	0.25	0.25	0.25	0.25	0.25	0.25	0.53	0.53	0.53	0.82	0.82...
    0.82	0.82	0.82	0.53	0.53	0.53	0.82	0.82	0.82	0.53	0.53	0.53];
C_sell = [0.22	0.22	0.22	0.22	0.22	0.22	0.22	0.42	0.42	0.42	0.65	0.65...
    0.65	0.65	0.65	0.42	0.42	0.42	0.65	0.65	0.65	0.42	0.42	0.42];
%------------------------变量定义-----------------------%
Pbuy = sdpvar(1, 24); %从电网购电电量
Psell = sdpvar(1, 24); %向电网售电电量
Pcha = sdpvar(1, 24);
Pdis = sdpvar(1, 24);
g = binvar(1, 24);
u = sdpvar(1, 24);
outerst = [u == Load + 20 .* g, sum(g) <= 8];%sum（g）<=8，表示保守度控制
innerst = [Pbuy - Psell + Pw + Ppv == u + Pcha - Pdis]; %内层约束
innerst = [innerst, ...
    0 <= Pbuy <= Temp_net .* 160, ...
    0 <= Psell <= (1 - Temp_net) .* 160, ...
    0 <= Pcha <= Temp_cha .* 40, ...
    0 <= Pdis <= Temp_dis .* 40, ...
    sum(Pcha - Pdis) == 0]; %主网功率交换约束
for t = 1 : 24
    innerst = [innerst, -30 <= sum(Pcha(1, 1 : t) - Pdis(1, 1 : t)) <= 165]; %SOC约束，电池容量300kwh，初始S0C为0.4，0.3<=SOC<=0.95
end



%% 目标函数
%------------------总费用--------------------%
obj = C_buy * Pbuy' - C_sell * Psell' - 0.2 * ones(1, 24) * Pdis';
[KKTsystem, details] = kkt(innerst, obj, u); %kkt命令，获取kkt条件
optimize([KKTsystem, outerst], -obj);
UB = value(obj);
Load_u = value(u);
Pbuy = value(Pbuy); %从电网购电电量
Psell = value(Psell); %向电网售电电量
Pcha = value(Pcha);
Pdis = value(Pdis);
X = [Load_u; Pbuy; Psell; Pcha; Pdis];
end

负荷预测场景（蓝）与负荷最劣场景（红）
负荷预测场景（蓝）与负荷最劣场景（红）

迭代收敛图
目标函数收敛图

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

yalmip

gurobi

算法

MATLAB

鲁棒优化（4）：通过yalmip中的kkt命令实现CCG两阶段鲁棒优化 的相关文章

  • 通过 h5py 将 matlab v7.3 文件读入 python numpy 数组列表

    我知道以前已经有人问过这个问题 但在我看来 仍然没有答案可以解释正在发生的事情 并且不适用于我的情况 我有一个 matlab v7 3 文件 其结构如下 gt rank lt 1x454 cell gt gt each element is

  • 有没有办法在 MATLAB 中查看 pcode 文件 (.p) 的源代码？

    有没有办法在 MATLAB 中打开 pcode 文件 p 如果 开放 是指edit 那么当然不是 pcode 中的 p 代表 受保护 其主要设计目标是在保护其源代码的同时部署功能组件 如果 开放 是指run 那么当然是的 引用手册 http

  • 图像分析-光纤识别

    我是图像分析新手 您知道如何以仅获取纤维的方式对该图像进行二值化吗 我尝试过不同的阈值技术等 但没有成功 我不介意应该使用什么工具 但我更喜欢 NET or Matlab PS 我不知道该把答案放在哪里 所以我把它放在StackOverfl

  • 继续在 Matlab 中一遍又一遍地播放声音？

    我正在尝试创建一个 MATLAB 程序来每隔几分钟一遍又一遍地播放声音 现在我将其设置为每隔几秒播放一次 只是为了消除系统中的一些错误 但是 当我的程序尝试重播声音时 我收到此错误 Error using gt audioplayer au

  • 如何在Matlab中自定义轮廓线？

    我正在准备一个等高线图 我应该在其中突出显示特定级别的等高线 例如 我的轮廓线值位于 1 和 1 之间 我想突出显示与值 0 相对应的线 我尝试使用以下过程来执行此操作 M c contourf longitude latitude del

  • Matlab：如何显示数组的“真实”值？

    我有一个在脚本中计算的向量 计算后 我将值显示到命令窗口 显示如下 finalResults 1 0e 05 0 0001 0 0 0005 0 0002 0 0001 0 0027 0 0033 0 0001 0 0000 0 0000

  • 按元素出现的频率对数组元素进行排序

    是否可以在 matlab octave 中使用sort函数根据元素的相对频率对数组进行排序 例如数组 m 4 4 4 10 10 10 4 4 5 应该产生这个数组 5 10 10 10 4 4 4 4 4 5是出现频率较低的元素 位于顶部

  • Matlab Mex文件编译

    我正在尝试编译一个 mex 文件以在 matlab 中使用套接字连接 问题是它总是说我没有安装sdk或编译器 但我已经安装了 Visual Studio 2010 Express Visual Studio 2012 Express Vis

  • MATLAB 滚动图

    我有一个脑电图数据库 我想绘制它 数据库是一个19 1000 134的矩阵 其中 19 是通道数 在第一种方法中 我只使用一个渠道 1000 个样本大小 采样率为 500 Hz 时为 1000 个点 即 2 秒数据 134 epochs的数

  • 在 matlab 中求 3d 峰的体积

    现在我有一个带有峰值的 3D 散点图 我需要找到其体积 我的数据来自图像 因此 x 和 y 值表示 xy 平面上的像素位置 z 值是每个像素的像素值 这是我的散点图 scatter3 x y z 20 z filled 我试图找到数据峰值的

  • 使用mat2cell将MxN的矩阵划分为1xN大小的M矩阵

    我有一个大小为 MxN 的矩阵 比方说 1867x3 1867 行和 3 列 我想将其分成 1867 个大小为 1x3 的单元格 我使用了mat2cell X 1 1866 这里X是矩阵 1867x3 结果给出了两个单元格 一个单元格的大小

  • 句柄类和值类的区别

    我有一些 C 背景 想使用 Matlab 中的类 句柄和值类有什么区别 我知道如果我想定义一个带有重载运算符 例如 和 的矩阵类 我会使用值类 然而 有时 当我选择一个手柄类时 事情似乎只对我有用 MathWorks 提供了一些有关其用途的

  • Simulink 仿真引擎如何工作？

    我想了解 Simulink 仿真引擎的工作原理 它是否使用离散事件模拟机制 那么如何处理连续时间 它是否依赖于基于静态循环的代码生成 或者 在第一个周期之前 它会计算出块的执行顺序 从不需要任何其他块输入的块开始 每个周期 它都会根据输入和

  • 在 MATLAB 中重命名文件

    我正在尝试以编程方式重命名工作目录中的文件a temp txt to b hello txt 您建议如何这样做 MATLAB中有一个简单的文件重命名函数吗 我认为您正在寻找 MOVEFILE

  • MATLAB 特征函数

    我很好奇哪里可以找到完整的描述FEATURE功能 它接受哪些论点 没有找到文档 我只听说过memstats and getpid 还要别的吗 gt gt which feature built in undocumented 注意 更完整的

  • Matlab 和 Python 中的优化算法（dog-leg trust-region）

    我正在尝试使用 Matlab 和 Python 中的狗腿信赖域算法求解一组非线性方程 在Matlab中有fsolve https www mathworks com help optim ug fsolve html其中此算法是默认算法 而

  • Matlab：如何更改矩阵的存储方式？从 1x1x3 到 1x3？

    我目前有 val 1 0 7216 val 2 0 7216 val 3 0 7216 但我想要 0 7216 0 716 0 721 我可以做什么样的操作来做到这一点 The reshape函数将在这里解决问题 Arrange the e

  • 在 Matlab 中高效获取像素坐标

    我想在 Matlab 中创建一个函数 给定一个图像 该函数将允许人们通过单击图像中的像素来选择该像素并返回该像素的坐标 理想情况下 人们能够连续单击图像中的多个像素 并且该函数会将所有相应的坐标存储在一个矩阵中 有没有办法在Matlab中做

  • 如何使用 MATLAB 的 substruct 函数创建表示使用“end”的引用的结构？

    我想使用substruct http www mathworks com help matlab ref substruct html函数创建一个结构体以供使用subsref 目的是使用索引字符串subsref而不是通常的 符号 因为我正在

  • 如何在 matlab 中创建由多个 3d 图像数据数组组成的数组

    我正在阅读 15 张图片imagedata imread imagename jpg 它的大小总是320 by 320 by 3 如何将数据放入数组中 使用 for for 循环 以便在访问新数组的第一个元素时获得输入的第一个图像的 RGB

随机推荐

  • Linux终端与SSH

    1 终端 当在我们需要操作服务器时要接上显示器和键盘 在键盘上输入 可以在显示器上看到 我们把 显示器 键盘 鼠标等这些统称为 输入 输出的终端设备 或者把显示器和键盘这些统称作 终端 当我们通过终端设备连接到服务器上并打开这些设备 如显示

  • ES6 新增的循环方法

    在 ES6 ECMAScript 2015 中 新增了一些循环方法 这些方法可以帮助我们更方便地遍历数组 字符串 Set Map 等数据结构 本文将介绍一些常用的 ES6 循环方法 for of 循环 for of 循环是一种遍历可迭代对象

  • 防火墙安全策略①

    防火墙 基本定义 防火墙是一种隔离 非授权用户在区域间 并过滤 对受保护网络有害流量或数据包 的设备 防火墙主要是借助硬件和软件的作用于内部和外部网络的环境间产生一种保护的屏障 从而实现对计算机不安全网络因素的阻断 只有在防火墙同意情况下

  • 如何将类数组转换为真正的数组

    开发过程中 有很多时候获取到的是类似数组的对象 比如元素的集合 elementcollection nodeList 以及今天开发时发现classList也是类数组 有时我们需要类数组去调用数组的方法 怎么办 一 遍历类数组 依次将元素放入

  • vue使用高德地图获取当前经纬度

    vue使用高德地图Api获取当前经纬度信息 在utils里面新建getLocation js 封装获取经纬度的公用方法 优化加载速度动态cdn引入高德地图 由于高德Api方法获取当前经纬度比较慢 如果需求是在获取到当前经纬度数据之后请求一些

  • JDBC数据库连接MySQL5.7

    1 首先准备mysql 和eclipse环境 在环境搭建好之后 从eclipse官网下载jdbc的驱动包 下载地址http dev mysql com downloads connector j 2 从下载的文件中取出mysql conne

  • 做电商的您还在为淘宝问题订单备注分类烦恼吗？

    做电商的您还在为淘宝问题订单备注分类烦恼吗 能够自动给淘宝订单添加旗子类型与备注的电商运营机器人 实现轻松管理淘宝问题订单 支持自动给订单添加旗子类型 自动给订单添加备注 解决问题订单处理缓慢 备注修改不及时等难题 使用UiBot之前 问题

  • Python javascript操作DOM

    文档对象模型 Document Object Model DOM 是一种用于HTML和XML文档的编程接口 它给文档提供了一种结构化的表示方法 可以改变文档的内容和呈现方式 我们最为关心的是 DOM把网页和脚本以及其他的编程语言联系了起来

  • 【元宇宙】智能手机万岁

    凭借出色的新设备 我们很快就能进人元字宙 想象这样的情景是很趣的 但是 至少到21世纪20年代 元宇宙时代的大多数设备很可能是我们已经在使用的设备 AR 和 VR 设备不仅面临重大的技术 财务和体验障碍 而且它们在上市后同样会面临币场反响冷

  • Esp8266 Node Mcu 一直乱码的问题详解

    最近一直在做项目 遇到的这个问题花了我很长时间 因此在这里写出自己的经历供大家参考 喜欢的可以点个赞 比较简单的方案 在Arduino上设置Node Mcu 1 打开文件 gt 首选项 复制这样一个网址 http arduino esp82

  • 关于Cache-Control: no-cache和no-store

    在公司上班的正真上班的第一天 发现的jsp页面上 设置了response HTTP头 是设置了这三个属性 Cache Control no cache Cache Control no store Expires 这三个属性都是和网页的缓存

  • OpenMPI:介绍与an'zhuang

    前言 在跑GLowdemo时发现没装OpenMPI 因此 只有装了吧 介绍 OpenMPI 1 是一种高性能消息传递库 最初是作为融合的技术和资源从其他几个项目 FT MPI LA MPI LAM MPI 以及 PACX MPI 它是MPI

  • pytest.fixture详解

    scope分为session package module class function级别 其中function也是包含测试类的测试方法的 package指在当前contest py的package下只执行一次 如果想在某个package

  • ＜深度学习基础＞ Batch Normalization

    Batch Normalization批归一化 BN优点 减少了人为选择参数 在某些情况下可以取消dropout和L2正则项参数 或者采取更小的L2正则项约束参数 减少了对学习率的要求 现在我们可以使用初始很大的学习率或者选择了较小的学习率

  • Unity3D中Update和FixedUpdate、LateUpdate的区别

    1 MonoBehaviour Update 更新 当MonoBehaviour启用时 其Update在每一帧被调用 2 MonoBehaviour FixedUpdate 固定更新 当MonoBehaviour启用时 其 FixedUpd

  • NOI2017搞基记

    一篇很长的流水账 写的长也就是因为考得好吧 去年写NOI游记的时候 就想着快点写完就好了 以后大概会再写一篇比较矫情的回顾一下竞赛的历程的吧 虽然我这种狗逼划水语文课代表的水平 大概激励人心的效果肯定赶不上hzwer的吧 DAY 2 在家排

  • linux下sublimetext的中文输入问题解决方法

    InputHelper插件 好奇怪哦 竟然一个文本编辑器在linux平台下竟然原生不能切换并使用系统自带的输入法 所以就有了一系列插件 看到网上各种方法 我觉得还是使用inputhelper这个插件最简单 使用这个插件可以通过Package

  • Java实现标题相似度计算，文本内容相似度匹配，Java通过SimHash计算标题文本内容相似度

    目录 一 前言 二 关于SimHash 补充知识 一 什么是海明距离 二 海明距离的应用 三 什么是编辑距离 三 SimHash算法的几何意义和原理 一 SimHash算法的几何意义 二 SimHash的计算原理 三 文本的相似度计算 四

  • Linux添加组播

    sudo route add net 224 1 1 0 netmask 255 255 255 0 dev ens33 转载于 https www cnblogs com tiandsp p 10985838 html

  • 鲁棒优化（4）：通过yalmip中的kkt命令实现CCG两阶段鲁棒优化

    两阶段鲁棒优化的原理推导部分 已经较多的文章进行分析 目前大部分同学面临的问题是 子问题模型中存在的双线性项该如何处理 目前 主流方式是 采用对偶定理或KKT条件 将第二阶段的双层问题变成单层问题 简略的思想如下 首先是原始的两阶段模型 对

热门标签