使用Yalmip+matlab求解主从博弈（双层规划）问题 | 教程（二）

yalmip官方学习网址：
https://yalmip.github.io/command

作为yalmip小白，本文主要根据实际代码运用，进行yalmip学习。
学习以yalmip官方网站为主，其他网站和教程为辅。

yalmip下的KKT使用

以某电动汽车充电背景下，代理商与电动汽车车主博弈的代码为例：

首先定义参数和变量

%参数
price_day_ahead=[0.35;0.33;0.3;0.33;0.36;0.4;0.44;0.46;0.52;0.58;0.66;0.75;0.81;0.76;0.8;0.83;0.81;0.75;0.64;0.55;0.53;0.47;0.40;0.37];
price_b=1.2*price_day_ahead;
price_s=1.2*price_day_ahead;
lb=0.8*price_day_ahead;
ub=1.2*price_day_ahead;
T_1=[1;1;1;1;1;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1];
T_2=[1;1;1;1;1;1;1;1;0;0;0;0;1;1;1;0;0;0;0;1;1;1;1;1];
T_3=[0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0];
index1=find(T_1==0);index2=find(T_2==0);index3=find(T_3==0);
%定义变量
Ce=sdpvar(24,1);%电价
z=binvar(24,1);%购售电状态
u=binvar(24,1);%储能状态
Pb=sdpvar(24,1);%日前购电
Pb_day=sdpvar(24,1);%实时购电
Ps_day=sdpvar(24,1);%实时售电
Pdis=sdpvar(24,1);%储能放电
Pch=sdpvar(24,1);%储能充电
Pc1=sdpvar(24,1);%一类车充电功率
Pc2=sdpvar(24,1);%二类车充电功率
Pc3=sdpvar(24,1);%三类车充电功率
S=sdpvar(24,1);%储荷容量

分为内外两层，对主从博弈下的双层规划问题进行求解：

%内层
CI=[sum(Pc1)==50*(0.9*24-9.6),sum(Pc2)==20*(0.9*24-9.6),sum(Pc3)==10*(0.9*24-9.6),Pc1>=0,Pc2>=0,Pc3>=0,Pc1<=50*3,Pc2<=20*3,Pc3<=10*3,Pc1(index1)==0,Pc2(index2)==0,Pc3(index3)==0];%电量需求约束
OI=sum(Ce.*(Pc1+Pc2+Pc3));
ops=sdpsettings('solver','gurobi','kkt.dualbounds',0);
[K,details] = kkt(CI,OI,Ce,ops);%建立KKT系统，Ce为参量

标准句式：

[KKTsystem, details] = kkt(Constraint,Objective,z) %z为变量，需要有上下限

对应即得，CI为约束（不同约束之间可直接用 ， 分隔），OI为目标函数，Ce为变量。

%外层
CO=[lb<=Ce<=ub,mean(Ce)==0.5,Pb>=0,Ps_day<=Pdis,Pb_day>=0,Pb_day<=1000*z,Ps_day>=0,Ps_day<=1000*(1-z),Pch>=0,Pch<=1000*u,Pdis>=0,Pdis<=1000*(1-u)];%边界约束
CO=[CO,Pc1+Pc2+Pc3+Pch-Pdis==Pb+Pb_day-Ps_day];%能量平衡
CO=[CO,sum(0.9*Pch-Pdis/0.9)==0,S(24)==2500,S>=0,S<=5000];%SOC约束
OO=-(details.b'*details.dual+details.f'*details.dualeq)+sum(price_s.*Ps_day-price_day_ahead.*Pb-price_b.*Pb_day);%目标函数
optimize([K,CI,CO,boundingbox([CI,CO]),details.dual<=1],-OO)

因此，整体代码分为内外两层，其针对的问题为：推导内部问题的KKT条件来求解双层二次规划。此时KKT无法导出对偶的边界，因此需要进行手动添加。原始边界需要通过KKT中的内部问题进行推导。

官网上的范例如下：

sdpvar x1 x2 y1 y2 y3
OO = -8x1-4x2+4y1-40y2-4*y3;
OO = OO+OO^2;
CO = [x1>=0, x2>=0];

OI = (x1+2x2+y1+y2+2y3)^2; CI = [[y1 y2 y3] >= 0,
-y1+y2+y3 <= 10,
2x1-y1+2y2-0.5y3 <= 10,
2x2+2y1-y2-0.5y3 <= 9.7];

[K,details] = kkt(CI,OI,[x1 x2])
optimize([K,CO,boundingbox([CI,CO]),details.dual<=100],OO)

同时boundingbox的标准句式如下：

[B,L,U] = boundingbox(Constraint,options,x)

optimize([K,CI,CO,boundingbox([CI,CO]),details.dual<=1],-OO)

其中boundingbox([CI,CO])表示两个约束
‘-oo’为目标函数