基于储能电站服务的冷热电多微网系统双层优化配置

1. 文章介绍
        摘要：随着储能技术的进步和共享经济的发展，共享储能电站服务模式将成为未来用户侧储能应用的新形态。提出基于储能电站服务的冷热电联供型多微网系统双层优化配置方法。首先，提出储能电站服务这种新型的共享储能模式，分析共享储能电站的运行方式和盈利机制。其次，将储能电站服务应用到冷热电联供型多微网系统中，建立考虑两个不同时间尺度问题的双层规划模型，上层模型负责求解长时间尺度的储能电站配置问题，下层模型负责求解短时间尺度的多微网系统优化运行问题。再次，根据下层优化模型的 Karush-Kuhn-Tucher(KKT)条件将下层模型转换为上层模型的约束条件，采用 Big-M 法对非线性问题线性化。最后，通过3个场景的算例分析验证所提双层规划模型的合理性和有效性，并证明所提出的共享储能服务能够有效降低用户成本，节约储能资源，实现用户与储能电站运营商的互利共赢。
1.1 储能电站服务模式

1.2 上层储能电站规划

1.2.1 目标函数
  

1.2.2. 约束条件

•     储能电站倍率约束 ​
•     · 储能电站荷电状态和充放电功率约束​ 

1.3 下层冷热电多微网系统优化运行模型
    1.3.1 目标函数 ​
    1.3.2 约束条件

•    电功率平衡约束 ​
•     冷功率平衡约束​ 
•    热功率平衡约束​
•   余热锅炉余热平衡约束​
•   储能电站充放电功率平衡约束​
•   微网设备出力上下限约束​
•   从电网购电功率约束​
• 微网与储能电站间购售电功率约束​

  2. 算例分析：
        1）场景 1，冷热电联供型多微网系统不配置储能，独立运行，多余的电能直接弃电，不足的电能从电网购买。
        2）场景 2，冷热电联供型多微网系统选择自行投资建设电储能装置，电储能充放电效率等参数同共享储能电站。
        3）场景 3，冷热电联供型多微网系统参与共享储能电站服务，使用储能电站的储能充放电服务。
  3. 代码

clc;
clear;
close all;

%% 选择用哪个求解器求解
ANSWER=questdlg('选择Yalmip使用的求解器:','电网技术论文复现',...
    'Cplex求解器','Gurobi求解器','Cplex求解器');

UseCplex=strcmp(ANSWER,'Cplex求解器');
UseGurobi=strcmp(ANSWER,'Gurobi求解器');

%% 模型参数设定
W=1;              %典型日个数为1
Tw=91;            %典型日天数为91
N=3;              %微网数为3
M=1E8;            %Big-M法中的M
NT=24;            %调度时段数为24；
%微网向储能电站的售电电价
sddj=[0.20,0.20,0.20,0.20,0.20,0.20,0.20,0.20,0.95,0.95,0.95,0.95,0.55,0.55,0.55,0.55,0.95,0.95,0.95,0.95,0.95,0.55,0.55,0.55];
%微网从储能电站购电的电价
gddj=[0.40,0.40,0.40,0.40,0.40,0.40,0.40,0.40,1.15,1.15,1.15,1.15,0.75,0.75,0.75,0.75,1.15,1.15,1.15,1.15,1.15,0.75,0.75,0.75];
%微网向储能电站缴纳的服务费单价
fwf=0.05*ones(1,24);   %0.05
%微网从电网的购电电价
dwdj=[0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37,1.36,1.36,1.36,1.36,0.82,0.82,0.82,0.82,1.36,1.36,1.36,1.36,1.36,0.82,0.82,0.82];

%% 决策变量初始化
Pessmax=sdpvar(1);     %储能电站最大充放电功率
Eessmax=sdpvar(1);     %储能电站的最大容量
Pessswi=sdpvar(3,24);  %第i个微网向储能电站的售电功率
Pessbwi=sdpvar(3,24);  %第i个微网向储能电站的购电功率
PGTwi=sdpvar(3,24);    %第i个微网的燃气轮机输出功率
Pgridwi=sdpvar(3,24);  %第i个微网从电网的购电功率
PECwi=sdpvar(3,24);    %第i个微网电制冷机消耗的电功率
QACwi=sdpvar(3,24);    %第i个微网制冷机的输出制冷功率
QGBwi=sdpvar(3,24);    %第i个电网的燃气轮机输出热功率
PHXwi=sdpvar(3,24);    %第i个电网的换热装置输出制热功率
Eess=sdpvar(1,24);     %储能电站存储的能量
Pessabs=sdpvar(1,24);  %储能电站充电功率
Pessrelea=sdpvar(1,24);%储能电站放电功率
Uabs=binvar(1,24);     %储能电站充放电状态，0-1变量
Urelea=binvar(1,24);   %储能电站充放电状态，0-1变量

%% 上层储能电站的目标函数
Cinvw=0.95*(Pessmax+Eessmax)/(8*365)+72*Pessmax/365;                       %储能电站日平均投资和维护成本
Cesssw=Pessswi(1,:)*sddj'+ Pessswi(2,:)*sddj'+ Pessswi(3,:)*sddj';         %每个典型日储能电站从微网购电成本
Cessbw=Pessbwi(1,:)*gddj'+ Pessbwi(2,:)*gddj'+ Pessbwi(3,:)*gddj';         %每个典型日储能电站向微网的售电效益
Cservew=Pessswi(1,:)*fwf'+ Pessswi(2,:)*fwf'+ Pessbwi(1,:)*fwf'+ Pessbwi(2,:)*fwf'...
    + Pessswi(3,:)*fwf'+ Pessbwi(3,:)*fwf';                                %每个典型日储能电站从微网收取服务费
objective=Tw*(Cinvw + Cesssw - Cessbw - Cservew);                          %上层优化目标为储能电站年运行成本最小

%% 上层储能电站的约束条件
Con=[Eessmax==2.6658*Pessmax,     %式（6）中储能电站容量和额定功率约束
    Eess(1)==0.2*Eessmax,];       %式（7）中储能电站初始存储的能量约束
for i=2:24
    Con=[Con,
        0.1*Eessmax<=Eess(i)<=0.9*Eessmax, %储能电站容量上下限，式7-3
        Eess(i)==(Eess(i-1)+(0.95*Pessabs(i)-(1/0.95)*Pessrelea(i))), %储能电站与上一时段的能量变化约束，式7-3
        ];
end
for i=1:24
    %这里做一点说明，和论文中式(7)不太一样，因为不能出现两个决策变量相乘，故这里用类似于大M法的思路代替
    Con=[Con,
        0<=Pessabs(i)<=Pessmax,     %式7-4
        0<=Pessabs(i)<=Uabs(i)*M,
        0<=Pessrelea(i)<=Pessmax,   %式7-5
        0<=Pessrelea(i)<=Urelea(i)*M,
        Uabs(i) + Urelea(i)<=1,       %确定充放电状态不可能同时发生，式7-6，7
        ];
end

%% 进行下层的KKT条件处理
%下层的三个微网的运行参数导入
%这里以春季作为典型日进行分析
Pcoolw=zeros(24,3);    %下层的3个微网的冷负荷功率
Pheatw=zeros(24,3);    %下层的3个微网的热负荷功率
Ploadw=zeros(24,3);    %下层的3个微网的电负荷功率
PPVw=zeros(24,3);      %下层的3个微网的光伏发电功率
PWTw=zeros(24,3);      %下层的3个微网的风力发电功率
%导入微网1的参数
Pcoolw(:,1)=xlsread('data','微网1','I2:I25');     
Pheatw(:,1)=xlsread('data','微网1','J2:J25');
Ploadw(:,1)=xlsread('data','微网1','K2:K25');
PPVw(:,1)=xlsread('data','微网1','L2:L25');
PWTw(:,1)=xlsread('data','微网1','M2:M25');
%导入微网2的参数(MG2中没有风电)
Pcoolw(:,2)=xlsread('data','微网2','J2:J25');     
Pheatw(:,2)=xlsread('data','微网2','K2:K25');
Ploadw(:,2)=xlsread('data','微网2','L2:L25');
PPVw(:,2)=xlsread('data','微网2','M2:M25');
%导入微网3的参数
Pcoolw(:,3)=xlsread('data','微网3','I2:I25'); 
Pheatw(:,3)=xlsread('data','微网3','J2:J25');
Ploadw(:,3)=xlsread('data','微网3','K2:K25');
PPVw(:,3)=xlsread('data','微网3','L2:L25');
PWTw(:,3)=xlsread('data','微网3','M2:M25');

% 下层等号约束导入
Con=[Con,              
    PGTwi+PWTw'+PPVw'+Pgridwi+Pessbwi-Pessswi-PECwi-Ploadw'==0, %式(11)电功率平衡约束
    4*PECwi+QACwi-Pcoolw'==0, %式(12)冷功率平衡约束
    QGBwi+PHXwi-Pheatw'==0,   %式(13)热功率平衡约束
    PHXwi/0.9+QACwi/1.2-1.47*0.8*PGTwi==0, %式(14)余热锅炉余热平衡约束
    Pessbwi(1,:)+Pessbwi(2,:)+Pessbwi(3,:)-Pessswi(1,:)-Pessswi(2,:)-Pessswi(3,:)==Pessrelea-Pessabs, %式(15)储能电站充放电功率平衡约束
    ];

%% 处理附录式（5）-（14）KKT条件
%定义等式约束的拉格朗日乘子
l1itw=sdpvar(3,24);   
l2itw=sdpvar(3,24);
l3itw=sdpvar(3,24);
l4itw=sdpvar(3,24);
l5itw=sdpvar(1,24);
%定义不等式约束的拉格朗日乘子
u1itmin=sdpvar(3,24);    %定义最小值不等式约束的拉格朗日乘子矩阵
u2itmin=sdpvar(3,24);
u3itmin=sdpvar(3,24);
u4itmin=sdpvar(3,24);
u5itmin=sdpvar(3,24);
u6itmin=sdpvar(3,24);
u7itmin=sdpvar(3,24);
u8itmin=sdpvar(3,24);
u1itmax=sdpvar(3,24);    %定义最大值不等式约束的拉格朗日乘子矩阵
u2itmax=sdpvar(3,24);
u3itmax=sdpvar(3,24);
u4itmax=sdpvar(3,24);
u5itmax=sdpvar(3,24);
u6itmax=sdpvar(3,24);
u7itmax=sdpvar(3,24);
u8itmax=sdpvar(3,24);
u9itmax=sdpvar(3,24);
Usalewi=binvar(3,24);      %定义0-1矩阵变量/布尔矩阵
Ubuywi=binvar(3,24);
%依次添加附录（5）-（14）KKT平衡条件约束
for i=1:3           
    for t=1:24
        Con=[Con,
            Tw*2.2/(0.3*9.7)+l1itw(i,t)-l4itw(i,t)*1.47*0.8+u1itmax(i,t)-u1itmin(i,t)==0,
            Tw*dwdj(t)+l1itw(i,t)+u6itmax(i,t)-u6itmin(i,t)==0,
            Tw*(gddj(t)+fwf(t))+l1itw(i,t)+u8itmax(i,t)-u8itmin(i,t)+l5itw(1,t)==0,
            Tw*(-sddj(t)+fwf(t))-l1itw(i,t)+u7itmax(i,t)-u7itmin(i,t)-l5itw(1,t)==0,
            -l1itw(i,t)+4*l2itw(i,t)+u3itmax(i,t)-u3itmin(i,t)==0,
            l3itw(i,t)+l4itw(i,t)/0.9+u5itmax(i,t)-u5itmin(i,t)==0,
            Tw*2.2/(0.3*9.7)+l3itw(i,t)+u4itmax(i,t)-u4itmin(i,t)==0,
            l2itw(i,t)+l4itw(i,t)/0.9+u2itmax(i,t)-u2itmin(i,t)==0,
            -u8itmax(i,t)*4000+u9itmax(i,t)==0,
            -u7itmax(i,t)*4000+u9itmax(i,t)==0,
            ];
    end
end
%处理附录式（15）-（31）KKT条件
%引入Big-M法所需的布尔变量
v15=binvar(3,24);
v16=binvar(3,24);
v17=binvar(3,24);
v18=binvar(3,24);
v19=binvar(3,24);
v20=binvar(3,24);
v21=binvar(3,24);
v22=binvar(3,24);
v23=binvar(3,24);
v24=binvar(3,24);
v25=binvar(3,24);
v26=binvar(3,24);
v27=binvar(3,24);
v28=binvar(3,24);
v29=binvar(3,24);
v30=binvar(3,24);
v31=binvar(3,24);
%Big-M法处理割平面约束
for i=1:3
    for t=1:24
        Con=[Con,
            0<=u1itmin(i,t)<=M*v15(i,t),
            0<=PGTwi(i,t)-5<=M*(1-v15(i,t)),
            0<=u1itmax(i,t)<=M*v16(i,t),
            0<=3000-PGTwi(i,t)<=M*(1-v16(i,t)),
            0<=u2itmin(i,t)<=M*v17(i,t),
            0<=QACwi(i,t)-5<=M*(1-v17(i,t)),
            0<=u2itmax(i,t)<=M*v18(i,t),
            0<=4000-QACwi(i,t)<=M*(1-v18(i,t)),
            0<=u3itmin(i,t)<=M*v19(i,t),
            0<=PECwi(i,t)-5<=M*(1-v19(i,t)),
            0<=u3itmax(i,t)<=M*v20(i,t),
            0<=4000-PECwi(i,t)<=M*(1-v20(i,t)),
            0<=u4itmin(i,t)<=M*v21(i,t),
            0<=QGBwi(i,t)-5<=M*(1-v21(i,t)),
            0<=u4itmax(i,t)<=M*v22(i,t),
            0<=4000-QGBwi(i,t)<=M*(1-v22(i,t)),
            0<=u5itmin(i,t)<=M*v23(i,t),
            0<=PHXwi(i,t)-5<=M*(1-v23(i,t)),
            0<=u5itmax(i,t)<=M*v24(i,t),
            0<=4000-PHXwi(i,t)<=M*(1-v24(i,t)),
            0<=u6itmin(i,t)<=M*v25(i,t),
            0<=Pgridwi(i,t)<=M*(1-v25(i,t)),
            0<=u6itmax(i,t)<=M*v26(i,t),
            0<=4000-Pgridwi(i,t)<=M*(1-v26(i,t)),
            0<=u7itmin(i,t)<=M*v27(i,t),
            0<=Pessswi(i,t)<=M*(1-v27(i,t)),
            0<=u7itmin(i,t)<=M*v28(i,t),
            0<=(4000*Usalewi(i,t)-Pessswi(i,t))<=M*(1-v28(i,t)),
            0<=u8itmin(i,t)<=M*v29(i,t),
            0<=Pessbwi(i,t)<=M*(1-v29(i,t)),
            0<=u8itmin(i,t)<=M*v30(i,t),
            0<=(4000*Ubuywi(i,t)-Pessbwi(i,t))<=M*(1-v30(i,t)),
            0<=u9itmax(i,t)<=M*v31(i,t),
            0<=(1-Ubuywi(i,t)-Usalewi(i,t))<=M*v31(i,t),
            ];
    end
end

%% 求解器相关配置
if UseCplex
    ops = sdpsettings('solver','cplex','verbose',2,'usex0',0);
    ops.cplex.mip.tolerances.mipgap = 1e-6;
end
if UseGurobi
    ops = sdpsettings('solver','gurobi','verbose',2,'usex0',0);
    ops.gurobi.MIPGap = 1e-6;
end

%% 进行求解计算         
result = optimize(Con,objective,ops);

if result.problem == 0 
    % problem =0 代表求解成功 
    % do nothing!
else
    error('求解出错');
end  
%文字输出运行结果
Pessmax=value(Pessmax);
Eessmax=value(Eessmax);
objective=value(objective);
display(['通过Yalmip求得的储能电站在该季度收益为 : ', num2str(-objective)]);
display(['储能电站在该季度的最优最大充放电功率为 : ', num2str(Pessmax)]);
display(['储能电站在该季度的最优容量为 : ', num2str(Eessmax)]);
%% 数据分析与画图
Eess=value(Eess);
Pessabs=value(Pessabs);
Pessrelea=value(Pessrelea);
%
for t=1:24
    Plot_SESS(1,t)=Pessabs(t);
    Plot_SESS(2,t)=-1*Pessrelea(t);
end
Plot_SESS=Plot_SESS';
figure
bar(Plot_SESS,'stacked');
hold on
plot(Eess,'r-o','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/kW');
title('共享储能电站充放电行为优化结果图');
hold on
legend('充电量','放电量','容量');
box off

%%结果展示-更新
figure
bar(Pcoolw,'stacked');
legend('MG1冷负荷','MG2冷负荷','MG3冷负荷');
xlabel('时间/h');
ylabel('功率/kW');
title('微网冷负荷参数');
box off

figure
bar(Pheatw,'stacked');
legend('MG1热负荷','MG2热负荷','MG3热负荷');
xlabel('时间/h');
ylabel('功率/kW');
title('微网热负荷参数');
box off

figure
bar(Ploadw,'stacked');
legend('MG1电负荷','MG2电负荷','MG3电负荷');
xlabel('时间/h');
ylabel('功率/kW');
title('微网电负荷参数');
box off

%%优化结果绘图
Pessswi=value(Pessswi);%第i个微网向储能电站的售电功率
Pessbwi=value(Pessbwi);%第i个微网向储能电站的购电功率
PGTwi=value(PGTwi); %第i个微网的燃气轮机输出功率
Pgridwi=value(Pgridwi);%第i个微网从电网的购电功率
PECwi=value(PECwi);%第i个微网电制冷机消耗的电功率
QACwi=value(QACwi);%第i个微网制冷机的输出制冷功率
QGBwi=value(QGBwi);%第i个电网的燃气轮机输出热功率
PHXwi=value(PHXwi);%第i个电网的换热装置输出制热功率
Eess=value(Eess);%储能电站存储的能量
Pessabs=value(Pessabs);%储能电站充电功率
Pessrelea=value(Pessrelea);%储能电站放电功率

figure
bar(-Pessswi(1,:));
hold on
bar(Pessbwi(1,:));
hold on 
legend('MG1向储能电站售电功率','MG1向储能电站购电功率');
xlabel('时间/h');
ylabel('功率/kW');
title('微网向储能电站购售电功率结果图');
box off

figure
bar(-Pessswi(3,:));
hold on
bar(Pessbwi(3,:));
hold on 
legend('MG3向储能电站售电功率','MG3向储能电站购电功率');
xlabel('时间/h');
ylabel('功率/kW');
title('微网向储能电站购售电功率结果图');
box off

figure
bar(-Pessswi(2,:));
hold on
bar(Pessbwi(2,:));
hold on 
legend('MG2向储能电站售电功率','MG2向储能电站购电功率');
xlabel('时间/h');
ylabel('功率/kW');
title('微网向储能电站购售电功率结果图');
box off

figure
bar(PGTwi(1,:),'stacked');%第i个微网的燃气轮机输出功率
% hold on
% bar(-Pgridwi(1,:),'stacked');%第i个微网从电网的购电功率
hold on 
plot(PECwi(1,:),'r-o','LineWidth',1.5)%第i个微网电制冷机消耗的电功率
hold on
plot(QACwi(1,:),'b-*','LineWidth',1.5)%第i个微网制冷机的输出制冷功率
hold on
plot(QGBwi(1,:),'g-*','LineWidth',1.5)%第i个电网的燃气轮机输出热功率
hold on
plot(PHXwi(1,:),'b-d','LineWidth',1.5)%第i个电网的换热装置输出制热功率
legend('燃气轮机输出功率','制冷机消耗的电功率','制冷机的输出制冷功率','燃气轮机输出热功率','换热装置输出制热功率');
xlabel('时间/h');
ylabel('功率/kW');
title('MG1调度结果');
box off

figure
bar(PGTwi(2,:),'stacked');%第i个微网的燃气轮机输出功率
% hold on
% bar(-Pgridwi(1,:),'stacked');%第i个微网从电网的购电功率
hold on 
plot(PECwi(2,:),'r-o','LineWidth',1.5)%第i个微网电制冷机消耗的电功率
hold on
plot(QACwi(2,:),'b-*','LineWidth',1.5)%第i个微网制冷机的输出制冷功率
hold on
plot(QGBwi(2,:),'g-*','LineWidth',1.5)%第i个电网的燃气轮机输出热功率
hold on
plot(PHXwi(2,:),'b-d','LineWidth',1.5)%第i个电网的换热装置输出制热功率
legend('燃气轮机输出功率','制冷机消耗的电功率','制冷机的输出制冷功率','燃气轮机输出热功率','换热装置输出制热功率');
xlabel('时间/h');
ylabel('功率/kW');
title('MG2调度结果');
box off

figure
bar(PGTwi(3,:),'stacked');%第i个微网的燃气轮机输出功率
% hold on
% bar(-Pgridwi(1,:),'stacked');%第i个微网从电网的购电功率
hold on 
plot(PECwi(3,:),'r-o','LineWidth',1.5)%第i个微网电制冷机消耗的电功率
hold on
plot(QACwi(3,:),'b-*','LineWidth',1.5)%第i个微网制冷机的输出制冷功率
hold on
plot(QGBwi(3,:),'g-*','LineWidth',1.5)%第i个电网的燃气轮机输出热功率
hold on
plot(PHXwi(3,:),'b-d','LineWidth',1.5)%第i个电网的换热装置输出制热功率
legend('燃气轮机输出功率','制冷机消耗的电功率','制冷机的输出制冷功率','燃气轮机输出热功率','换热装置输出制热功率');
xlabel('时间/h');
ylabel('功率/kW');
title('MG3调度结果');
box off

4. 部分程序运行结果