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

2023-12-20 21:36:49

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

代码链接: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 确定性优化模型

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

式中 xy 为优化变量,具体表达式为:

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

参数

参数说明

目标函数(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 运行可视化

文章来源:https://blog.csdn.net/qq_41921826/article/details/135036868
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。