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

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

目前，主流方式是，采用对偶定理或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);

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

主问题如下：

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

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



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

子问题如下：

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

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

迭代收敛图