程序员经验分享-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两阶段鲁棒优化 的相关文章

  • 单元格的 Fieldnames 函数的等效项

    正如标题所说 只是想知道是否有一个函数可以用作字段名 http www mathworks co uk help matlab ref fieldnames html 但适用于单元格 所以如果我有类似的东西 a imread redsqua

  • 为什么 mex 文件中的 OpenMP 仅产生 1 个线程？

    我是 OpenMP 新手 我有以下代码 使用配置了 MSVS2010 的 Matlab mex 可以正常编译 计算机有 8 个可用处理器 我也使用 matlabpool 检查过 include mex h include

  • 归一化互相关的基础知识

    我正在尝试使用范数校正2 归一化互相关 http en wikipedia org wiki Cross correlation Normalized cross correlation 来自 MATLAB 用于计算发育中胚胎中移动形状的速

  • MATLAB 可执行文件太慢

    我使用以下命令将 MATLAB 程序转换为基于控制台的应用程序deploytool在 MATLAB 中 MATLAB m文件执行大约需要 2 秒 但在我将其转换为可执行文件并调用 exe 执行需要45秒 太长了 我想将 MATLAB 程序与

  • 动态调整自定义刻度数

    Taking SO 的一个例子 https stackoverflow com a 7139485 97160 我想根据当前视图调整轴刻度 这是默认行为 除非设置自定义的刻度数 下图展示了由此产生的行为 左侧是默认行为 右侧是带有自定义刻度

  • 如何为已编译的 MATLAB 创建安装程序并要求用户接受我们的许可条款？

    我正在 MATLAB 中编写程序分发给 Windows 用户 我使用 MATLAB 编译器和 MATLAB r2014a 版本来创建程序 我可以使用 MATLAB 应用程序编译器创建 Windows 安装程序 并且它的工作效果可以接受 但是

  • 直方图均衡结果

    I am trying to code histogram equalization by my self but the results are different from the built in function in matlab

  • 如何在向量中的所有点之间绘制线？

    我有一个包含二维空间中一些点的向量 我希望 MATLAB 用从每个点到每个其他点绘制的线来绘制这些点 基本上 我想要一个所有顶点都连接的图 你能用情节来做到这一点吗 如果可以 怎么做 一种解决方案是使用该函数为每个点组合创建一组索引MESH

  • Numpy 相当于 MATLAB 的 hist [重复]

    这个问题在这里已经有答案了 由于某种原因 Numpy 的 hist 总是返回比 MATLAB 的 hist 少 1 个 bin 例如在 MATLAB 中 x 1 2 2 2 1 4 4 2 3 3 3 3 Rep Val hist x un

  • FMINCON 的替代方案

    除了 fmincon 之外还有其他更快 更高效的求解器吗 我正在使用 fmincon 来解决特定问题 但对于中等大小的向量变量来说 我的内存不足 我也没有任何超级计算机或云计算选项可供使用 我知道任何替代解决方案仍然会耗尽内存 但我只是想看

  • 如何将数据传递给 MATLAB oncleanup 函数？

    我有一个编译好的 matlab 程序 可以自动调整机器参数 在调整周期结束时 我需要恢复一些原始设置 有时会发生意外错误 有时用户会发现调整算法未正常工作 因此应终止 使用 control C 如果发生可预测的错误 我可以使用 try ca

  • 如何选择面积最大的对象？

    我用过bwconvhull检测图像的某个部分 正如您在图像中看到的那样 有许多具有特定质心的对象 我想做的是检测面积最大的物体 左起第一个大物体 并忽略其他物体 我应该遵循哪种方法 我将非常感谢您的帮助 以下是代码 由于我仍在努力 所以写得

  • 使用简单矩阵乘法时出错

    我在一次简单的乘法运算中偶然发现了一个错误 这让我感到非常惊讶 我一直以为这里发生了什么 只为矩阵乘法 http www mathworks nl help matlab matlab prog operators html x 2 y z

  • MATLAB 变量传递和惰性赋值

    我知道在 Matlab 中 当将新变量分配给现有变量时 会进行 惰性 评估 例如 array1 ones 1 1e8 array2 array1 的价值array1不会被复制到array2除非元素array2被修改 由此我推测Matlab中

  • 帮助我理解FFT函数（Matlab）

    1 除了负频率之外 FFT 函数提供的最小频率是多少 是零吗 2 如果它为零 我们如何在对数刻度上绘制零 3 结果总是对称的 或者只是看起来是对称的 4 如果我使用abs fft y 来比较2个信号 我是否会失去一些准确性 1 除了负频率之

  • 优先连接，Matlab 中的复杂网络

    大家好 我现在正在 MATLAB 中研究优先附件模型 在理解以下内容时遇到一些困难 假设我一开始有 4 个节点 连接如下 time 0 1 lt gt 2 3 lt gt 4 在下一个时间步骤中 我添加一个节点和 4 个连接 然后添加另一个

  • 2D 网格的纹理贴图

    我有一组点 x y meshgrid 1 N 1 M 在常规二维上定义 N x M网格 我还有另一组要点 u v 这是原始网格的一些变形 即 u v f x y 但是我没有实际的f导致变形 如何将纹理映射到由定义的 变形 网格u v 即 给

  • 像matlab一样在python中连接数组而不知道输出数组的大小

    我正在尝试在 python 中连接数组 类似于 matlab array1 zeros 3 500 array2 ones 3 700 array array1 array2 我在 python 中做了以下操作 array1 np zero

  • MATLAB parfor 和 C++ 类 mex 包装器（需要复制构造函数？）

    我正在尝试使用概述的方法将 C 类包装在 matlab mex 包装器中here http www mathworks com matlabcentral newsreader view thread 278243 基本上 我有一个初始化

  • 绘制布朗运动 matlab

    首先 我只想说我不太习惯使用matlab 但我需要一个作业 我应该创建一个 布朗运动 我的代码目前如下所示 clf hold on prompt Ge ett input size input prompt numParticles inp

随机推荐

  • 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条件 将第二阶段的双层问题变成单层问题 简略的思想如下 首先是原始的两阶段模型 对

热门标签