两阶段鲁棒优化的原理推导部分,已经较多的文章进行分析。目前大部分同学面临的问题是,子问题模型中存在的双线性项该如何处理?
目前,主流方式是,采用对偶定理或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
负荷预测场景(蓝)与负荷最劣场景(红)
迭代收敛图