微电网优化(Matlab复现)— 微电网两阶段鲁棒优化经济调度方法_刘一欣

论文链接：微电网两阶段鲁棒优化经济调度方法 - 中国知网

代码链接：https://m.tb.cn/h.5Mg7fCo?tk=hnpmWgZiv2R 

复现效果：

运行环境：Matlab 2020b+Cplex+yalmip

1 微电网结构

图 1 所示为典型的微电网结构，由可控分布式电源、可再生分布式电源、储能及本地负荷集成而成。 此外， 考虑微电网内包含需求响应负荷的情况，微电网可通过灵活调整需求响应负荷的用电计划，降低运行成本。

图 1 微电网结构

2 系统模型

2.1 微型燃气轮机

微型燃气轮机的成本函数，

微型燃气轮机的约束条件，由于微型燃气轮机的功率响应速度相对于小时级调度而言较快，
因此不考虑其爬坡率约束，仅考虑输出功率约束：

输出功率最大/最小约束：

2.2 储能

储能的成本函数，

储能的约束条件，

储能的充电功率最大/最小约束：
储能的放电功率最大/最小约束：
储能在调度的始末时刻容量相等：
储能各时段的剩余容量约束：

2.3 需求响应负荷

需求响应负荷的成本函数，式(10)中的绝对值项用于表示实际调度功率和期望用电功率之间的偏差，通过引入辅助变量及约束可将其化为式(11)所示的线性形式。

需求响应负荷的约束条件，

需求响应负荷最大/最小约束：
需求响应功率平衡约束：
需求响应功率最小约束：

2.4 配电网交互功率

配电网之间的交互功率的成本函数，

配电网之间的交互功率的约束条件，

功率平衡约束：
微电网向配电网购买的功率约束：
微电网向配电网出售的功率约束：

3 两阶段鲁棒优化模型

3.1 目标函数

微电网的运行目标为日运行成本最小化，如式(18)所示，包括上述的约束条件。

3.2 确定性优化模型

当不考虑光伏出力和负荷功率的不确定性时，可得到上述微电网经济调度问题的确定性优化模
型，其紧凑形式可表述为：

式中 x、 y 为优化变量，具体表达式为：

在确定性模型中，光伏出力和负荷功率的为预测值，并没有考虑预测误差的影响。

参数	参数说明
	目标函数(18)对应的系数列向量
	对应约束下变量的系数矩阵
	常数列向量
	不等式约束，包括式(2)、(7)、(9)、(13)
	等式约束，等式不一定为0，代码中为Ky=g，包括式(6)、(8)、(12)、(14)
	双向约束，包括式(4)、(5)、(15)、(16)
	确定性优化模型约束，光伏出力和负荷功率的取值为各时段相应的预测值
	光伏出力、负荷功率的预测值

3.3 不确定性优化模型

确定性优化模型得到的方案往往显得过于“冒险”，需要在模型中计及不确定性的影响。考虑光伏出力和负荷功率的波动范围位于式(22)所构建的箱型不确定集 U 内：

文章搭建的两阶段鲁棒优化模型的目的在于找到不确定变量u在不确定集U内朝着最恶劣场景变化时经济性最优的调度方案，具有如下形式：

公式理解：外层的最小化为第一阶段问题，优化变量为x；内层的最大最小化为第二阶段问题， 优化变量为 u 和 y，其中的最小化问题等同于式(19)的目标函数，表示最小化运行成本；max 结构的目的就在于找到导致运行成本最大的最恶劣场景。

整体的思路就是给定一组(x,u)，但此时的u是不确定的，要考虑光伏和负荷的预测误差，考虑到此时最恶劣的场景也就确定了u，确定了(x,u)后此时的问题也就转变为确定性的场景，最后求解y，使运行成本最小化，得到各个设备的出力结果。

4 求解推导

 对偶理论理解：我们在求解最优化问题时，常常不知道什么时候优化结束，也不知道该问题的最优值究竟是多少，求解的最优值可能与真实的最优值差距很大。对偶理论在一定程度上为我们获得的这个解提供一个下界。例如，我们求解一个最优化问题时，发现目标函数值迭代到100附近就迭代不动了。这时，你去求解它的对偶问题，发现对偶问题的最优函数值是95。可以确定原问题的最优值一定大于95，目前已经求解到100，可以停止迭代求解。关于对偶理论的理解可以参考如下的文章和视频：

文章：最优化3: 对偶 - 知乎

视频：最优化理论与方法-第十讲-约束优化（三）：对偶理论_哔哩哔哩_bilibili

原问题：

对原问题进行分解，得到的主问题形式为：                 

进一步将原问题表达为如下形式：

引入各约束的对偶变量：

转化为标准形式：

写出Lagrange函数：

   进一步求Lagrange函数的最值：

如果可行，这里如果为等式约束，对偶变量没有大小限制，作者对该问题做了修正，文章中的可以去除，推出：

最后将其转化为 max 形式，并与外层的 max 问题合并，最后得到如下对偶问题：

该对偶问题最优解所对应的 u*为不确定集 U 的一个极点，也就是说，式(27)取到最大值时，不确定变量 u 的取值应为上式所描述的波动区间的边界。此外，在本文所研究的微电网中，光伏出力取到区间的最小值和负荷功率取到区间的最大值时，微电网的运行成本更高，更符合“最恶劣”场景的定义，因此可将式(22)改写成如下形式：

将式(28)中的不确定变量表达式代入式(27)后，将出现二进制变量和连续变量乘积的形式，通过引入辅助变量和相关约束对其进行线性化，最后得到：

5 代码详解

5.1 设置参数

%% 设参
%微型燃气轮机
pm_max=1500;    %联络线功率上限
eta=0.95;       %储能充放电效率
p_g_max=800;    %出力上限
p_g_min=80;     %出力下限
a=0.67;         %成本系数
b=0;            %成本系数
%储能
ps_max=500;     %充放电功率上限
ES_max=1800;    %荷电状态上限
ES_min=400;     %荷电状态下限
ES0=1000;       %初始荷电状态
KS=0.38;        %单位充放电成本
%需求响应
DDR=2940;       %总用电需求
DR_max=200;     %最大用电需求
DR_min=50;      %最小用电需求
KDR=0.32;       %单位调度成本
%分时电价
price = [0.48;0.48;0.48;0.48;0.48;0.48;0.48;0.9;1.35;1.35;1.35;0.9;0.9;0.9;0.9;0.9;0.9;0.9;1.35;1.35;1.35;1.35;1.35;0.48];
%风电、光伏、负荷预测值
PW_=[0.6564    0.6399    0.6079    0.5594    0.5869    0.5794    0.6138    0.6192   0.6811    0.6400    0.7855    0.7615    0.6861    0.8780    0.6715    0.7023    0.6464    0.6321    0.6819    0.6943    0.7405    0.6727    0.6822    0.6878];
p_pv=1500*[     0         0         0         0         0    0.0465    0.1466    0.3135     0.4756    0.5213    0.6563    1.0000    0.7422    0.6817    0.4972    0.4629    0.2808    0.0948    0.0109         0         0         0         0         0];
PL=1500*[ 0.4658    0.4601    0.5574    0.5325    0.5744    0.6061    0.6106    0.6636    0.7410    0.7080    0.7598    0.8766    0.7646    0.7511    0.6721    0.5869    0.6159    0.6378    0.6142    0.6752    0.6397    0.5974    0.5432    0.4803];
%需求响应负荷的期望用电功率  
P_DR=1*[110 100 90 80 100 100 130 100 120 160 175 200 140 100 100 120 140 150 190 200 200 190 80 60];

决策变量也就是待求的变量

%% 设决策变量
%储能
p_ch=sdpvar(1,24,'full');   %储能充电
p_dis=sdpvar(1,24,'full');  %储能放电
us=binvar(1,24,'full');     %充放电标识
%配电网
p_buy=sdpvar(1,24,'full');  %配网购电
p_sell=sdpvar(1,24,'full'); %配网售电
um=binvar(1,24,'full');     %购售电标识
%分布式电源
%p_pv=sdpvar(1,24,'full');  %光伏发电
%pL=sdpvar(1,24,'full');    %固定日负荷
p_g=sdpvar(1,24,'full');    %分布式电源
%负荷
PDR=sdpvar(1,24,'full');    %可转移负荷
PDR1=sdpvar(1,24,'full');   %可转移负荷辅助变量
PDR2=sdpvar(1,24,'full');   %可转移负荷辅助变量
afa=sdpvar(1,1,'full');

5.2 约束条件

%% 约束条件
%分布式电源约束
C=[-Q1*y>=-p_g_max];                    %分布式电源的功率小于最大功率
C=C+[Q1*y>=p_g_min];                    %分布式电源的功率大于最小功率

%储能约束
C=C+[-Q31*y-ps_max*Q01*x>=-ps_max];     %储能充电功率不能超过最大充电功率
C=C+[-Q32*y>=-Q01*x*ps_max];            %储能放电功率不能超过最大放电功率
C=C+[Q31*y>=0];                         %储能充电功率必须非负
C=C+[Q32*y>=0];                         %储能放电功率必须非负
C=C+[Q2*y==0];                          %保证储能在调度前后能量相同
C=C+[-Q4*y>=-(ES_max-ES0)];             %储能量不能超过最大储能量
C=C+[Q4*y>=(ES_min-ES0)];               %储能量不能低于最小储能量

%配电网交互约束
C=C+[-Q52*y-pm_max*Q02*x>=-pm_max];     %出售电力的功率减去电力市场交易的功率不小于最大功率
C=C+[-Q51*y>=-Q02*x*pm_max];            %出售电力的功率不小于电力市场交易的功率乘以最大功率
C=C+[Q51*y>=0];                         %出售电力的功率大于等于零
C=C+[Q52*y>=0];                         %购买电力的功率大于等于零
C=C+[Q6*y+G1*u==0];                     %功率平衡的约束条件

%可转移负荷约束
C=C+[Q8*y==DDR];                        %表示可转移负荷的功率需求的上下限约束，确保可转移负荷的功率在一定范围内
C=C+[-Q9*y>=-DR_max];                   %可转移负荷的功率小于最大用电需求
C=C+[Q9*y>=DR_min];                     %可转移负荷的功率大于最小用电需求
C=C+[Q101*y>=0];                        %大于等于零
C=C+[Q102*y>=0];                        %大于等于零
%C=C+[Q9*y+Q101*y-Q102*y=P_DR];
C=C+[Q103*y==P_DR'];

5.3 目标函数

求解问题可以表述为如下的紧凑形式，然后构建紧凑形式的系数矩阵，进一步通过矩阵来表达求解的问题。

%% 紧凑形式
%cy
%Dy>=d
%Ky=g
%Fx+Gy>=h
%Ly+Yu=0

%紧凑后的目标函数
c=QC;           
%紧凑后的约束条件
C=C+[D*y>=d];  
C=C+[K*y==g];
C=C+[F*x+G*y>=h];
C=C+[L*y+Y*u==0];
C=C+[afa>=c*y];

%系数矩阵定义
x=[us,um]';                                     %第一阶段变量，使用转置操作将其转换为列向量
y=[p_g,p_ch,p_dis,PDR,PDR1,PDR2,p_buy,p_sell]';    %第二阶段变量，同样使用转置操作

D=[-Q1;Q1;Q31;Q32;-Q4;Q4;Q51;Q52;-Q9;Q9;Q101;Q102];
d=[-p_g_max*ones(24,1);p_g_min*ones(24,1);0*ones(24,1);0*ones(24,1);-(ES_max-ES0)*ones(24,1);(ES_min-ES0)*ones(24,1);0*ones(24,1);0*ones(24,1);-DR_max*ones(24,1);DR_min*ones(24,1);0*ones(24,1);0*ones(24,1)];

K=[Q2;Q8;Q103];
g=[0;DDR;P_DR'.*ones(24,1)];

F=[-ps_max*Q01;ps_max*Q01;-pm_max*Q02;pm_max*Q02];
G=[-Q31;-Q32;-Q52;-Q51];
h=[-ps_max*ones(24,1);0*ones(24,1);-pm_max*ones(24,1);0*ones(24,1)];

L=[Q6];
Y=[G1];

决策变量x是一个48维的列向量：

x=[US,UM]T

决策变量y是一个240维的列向量： 

y=[PG,PSch,PSdis,PDR,PDR1,PDR2,PMbuy,PMsell,PPV,PL]T

D是一个144行240列的矩阵，d是一个144维的列向量，E24表示24阶的单位矩阵，L24表示下三角和对角线全为1的矩阵，上三角全为0的矩阵。Dy≥d写成矩阵形式： 

 K是一个50行240列的矩阵，k是一个50维的列向量，Ky=g写成矩阵形式：

F是一个96行48列的矩阵，G是一个96行240列的矩阵，h为96维的列向量，Fx+Gy≥h写成矩阵形式：

是一个48行240列的矩阵，是48维的列向量，写成矩阵形式：

5.4 问题求解

%% 迭代求解两阶段鲁棒优化问题
% 先运行一次
[x,LB,y] = MP2();   %调用函数MP2
[u,UB] = SP(x);     %调用函数SP()，传入MP2()中返回的参数x
UB1 = UB;         %将函数SP()的返回值UB的值赋给UB1
p(1)= UB1 - LB;     %计算UB1-LB，并将结果赋给数组p的第一个元素。
%开始迭代
for k=1:4
    [x,LB,y] = MP(u);    %调用函数MP()
    [u,UB] = SP(x);      %调用函数SP()，传入MP2()中返回的参数x
    UB = min(UB1,UB);  %取UB1和UB中较小的值
    p(k+1) = UB-LB;     %计算UB-LB，并将结果赋给数组p的其他元素
end

5.5 运行可视化