代码链接: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,可以停止迭代求解。关于对偶理论的理解可以参考如下的文章和视频:
原问题:
对原问题进行分解,得到的主问题形式为:
进一步将原问题表达为如下形式:
引入各约束的对偶变量:
转化为标准形式:
写出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