鲁棒优化（4）：通过yalmip中的kkt命令实现CCG两阶段鲁棒优化
2024-04-01 01:02:01
两阶段鲁棒优化的原理推导部分，已经较多的文章进行分析。目前大部分同学面临的问题是，子问题模型中存在的双线性项该如何处理？
目前，主流方式是，采用对偶定理或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:6if i==1Load_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变量endLB_recoder = [LB_recoder, LB];[UB, X] = SP(Temp_net, Temp_cha, Temp_dis); %传入01变量，获得最坏的场景LoadUB_recoder = [UB_recoder, UB];if UB-LB<=3break;end
endplot(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 : 24st = [st, -30 <= sum(Pcha(1, 1 : t) - Pdis(1, 1 : t)) <= 165]; %SOC约束，电池容量300kwh，初始S0C为0.4，0.3<=SOC<=0.95
endst = [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);
endfunction [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 : 24st = [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 : 24innerst = [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负荷预测场景（蓝）与负荷最劣场景（红）
迭代收敛图
鲁棒优化（4）：通过yalmip中的kkt命令实现CCG两阶段鲁棒优化相关推荐
- 鲁棒优化入门(4)-两阶段鲁棒优化及行列生成算法(CCG)超详细讲解(附matlab代码)
本文的主要参考文献: Zeng B , Zhao L . Solving Two-stage Robust Optimization Problems by A Constraint-and-Colu ...
- 基于两阶段鲁棒优化的微网电源储能容量优化配置
目录 1 主要内容 1.1 目标函数 1.2 约束条件 1.3 不确定变量 2 部分代码 3 程序效果 4 程序链接 1 主要内容 程序主要复现的是<考虑寿命损耗的微网电池储能容量优化配置> ...
- 电网两阶段鲁棒优化调度模型（含matlab程序）
目录 一 两阶段鲁棒优化理论 二 两阶段鲁棒优化程序实现 2.1 主/子问题变量要分清 2.2 对偶问题 2.3 线性化处理 2.4 编程小技巧 2.5 迭代问题 三 程序运行效果 视频讲解 两阶段鲁 ...
- 【重磅】IEEE33配电网两阶段鲁棒优化调度CCG
目录 1 前言 2基本内容 2.1 配网两阶段鲁棒模型 2.2 求解步骤 3部分程序 4程序结果 5程序链接 1 前言 鲁棒优化是电力系统研究的热点,而两阶段鲁棒和分布鲁棒研究就成为各类期刊(sci/ ...
- 基于CCG算法的IEEE33配电网两阶段鲁棒优化调度matlab
目录 1 前言 2基本内容 2.1 配网两阶段鲁棒模型 2.2 求解步骤 3部分程序 4程序结果 5程序链接 1 前言 鲁棒优化是电力系统研究的热点,而两阶段鲁棒和分布鲁棒研究就成为各类期刊(sci ...
- 微电网两阶段鲁棒优化经济调度方法(完美复现)
针对微电网内可再生能源和负荷的不确定性,建立了 min-max-min 结构的两阶段鲁棒优化模型,可得到最恶劣场 景下运行成本最低的调度方案.模型中考虑了储能.需求侧 负荷及可控分布式电源等的运行约束 ...
- 9微电网两阶段鲁棒优化经济调度方法（MATLAB程序）
15R 复现文章: 微电网两阶段鲁棒优化经济调度方法--刘一欣(中国电机工程学报) 主要内容: 针对微电网内可再生能源和负荷的不确定性,建立了min-max-min 结构的两阶段鲁棒优化模型,可得到最 ...
- 两阶段鲁棒优化程序 采用微网为模型，主要将安装成本、运营成本以及综合效益三个方面纳入考虑范围
两阶段鲁棒优化程序 采用微网为模型,主要将安装成本.运营成本以及综合效益三个方面纳入考虑范围,建立两阶段鲁棒优化模型,采用的是CCG方法,本程序为matlab编制,有售后,可以进行!另外本程序考虑发电 ...
- 采用CCG和kkt条件编制两阶段鲁棒优化程序，以储能、发电、风电和光伏容量作为第一阶段变量
#电力系统##微电网#微网两阶段鲁棒优化matlab版 采用CCG和kkt条件编制两阶段鲁棒优化程序,以储能.发电.风电和光伏容量作为第一阶段变量,以主体出力作为第二阶段变量,以负荷.风电和光伏出力作 ...
最新文章
- 2022-2028年中国刨花板市场投资分析及前景预测报告(全卷)
- C++运算符重载（成员函数方式）
- 【Android】进入Material Design时代
- Sumdiv POJ - 1845
- jquery-文档宽高-窗口宽高-scrollTop与scrollLeft
- JQuery + XML作为前后台数据交换格式实践
- ubuntu系统debootstrap的再三实验
- 已解决：mysql: [Warning] Using a password on the command line interface can be insecure.
- 动态迁移到另一点出现电源关闭的原因及解决方法
- 【路径规划】基于matlab麻雀算法求解机器人栅格地图最短路径规划问题【含Matlab源码 1582期】
- 2019 CCF 推荐 会议 列表
- 告别手机自带浏览器，分享2022年好用的手机浏览器
- 利用偏最小二乘法选出最重要的特征波段Matlab
- 系统论重要模型_笔记
- 女人不适合做产品经理？
- python回归分析结果怎么看_Python数据分析｜线性回归
- OpenHarmony恢复启动子系统init进程之服务管理与发布
- 基本面分析中必须了解的88条避雷常识
- MG本地服务器如何改ip
- 在IDEA中给项目同时配置git和svn