复现文章：

基于共享储能电站的工业用户日前优化经济调度——李淋（电力建设2020）

摘要：

文章提出一种基于共享储能电站的工业用户日前优化经济调度方法。首先提出共享储能电站的概念，分析其商业运营模式。然后将共享储能电站应用到工业用户经济优化调度中，通过协调各用户使用共享储能电站进行充电和放电的功率，实现用户群日运行成本最优。最后以江苏省3 个工业用户进行算例仿真，与用户不配置储能和用户独立配置储能场景对比，得出引入共享储能电站可以显著降低用户群日运行成本，并对储能电站年服务费收益、静态投资回收年限和投资回报率与共享储能电站服务费定价间的关系做进一步的研究。

部分程序：

clc
clear
close all

%% 参数设置
P_grid=sdpvar(3,24); %用户从电网购买的电功率
P_ess_b=sdpvar(3,24); %用户使用自身内部独立储能的放电的功率
P_ess_s=sdpvar(3,24); %用户使用自身内部独立储能的充电的功率
U_ess_b=binvar(3,24); %用户使用自身内部独立储能的放电状态位,取1时为放电,0为未放电
U_ess_s=binvar(3,24); %用户使用共自身内部独立储能的充电状态位,取1时为充电,0为未充电
E=sdpvar(3,24); %用户自身内部独立储能的荷电状态
E_init=sdpvar(3,1); %用户自身内部独立储能的初始容量
E_max=sdpvar(3,1); %用户自身内部独立储能的预配置容量
P_max=sdpvar(3,1); %用户自身内部独立储能的预配置最大充放电功率
%% 用户参数导入
gamma=[0.37*ones(1,8),1.36*ones(1,4),0.82*ones(1,5),1.36*ones(1,4),0.82*ones(1,3)]; %电网的"峰-平-谷"电价
P_load(1,:)=[80,75,75,75,75,75,75,75,95,155,160,180,160,165,150,165,160,175,160,130,140,100,90,85];
P_load(2,:)=[200,50,150,190,25,50,30,20,40,170,220,150,65,160,170,150,230,125,140,200,125,70,110,65];
P_load(3,:)=[105,100,95,100,100,105,110,110,120,110,115,120,135,110,120,125,130,125,115,130,120,120,115,120];
P_pv(1,:)=[0,0,0,0,0,0,0,0,75,275,540,675,750,750,330,170,75,10,0,0,0,0,0,0];
P_pv(2,:)=zeros(1,24);
P_pv(3,:)=zeros(1,24);
P_wind(1,:)=zeros(1,24);
P_wind(2,:)=[150,150,175,165,140,130,120,90,50,55,80,100,135,135,130,110,60,50,50,70,90,120,140,120];
P_wind(3,:)=[130,145,155,135,120,125,115,95,40,30,50,85,110,115,110,70,35,40,40,50,80,80,115,95];
%% 约束条件
C=[]; %约束条件矩阵初始化
for i=1:3
for t=1:24
C=[C,
P_pv(i,t)+P_wind(i,t)+P_grid(i,t)+P_ess_b(i,t)-P_ess_s(i,t)-P_load(i,t)==0, %电功率平衡约束
U_ess_b(i,t)+U_ess_s(i,t)<=1, %避免同时出现充放电的情况,故增加此约束
P_grid(i,t)>=0,
]; %电功率平衡约束
end
end
%储能电站荷电状态连续性约束
for i=1:3
C=[C,E(i,1)==E_init(i)+0.95*P_ess_s(i,1)-P_ess_b(i,1)/0.95,]; %1时段和其他时段不同,涉及E(0)
for t=2:24
C=[C,E(i,t)==E(i,t-1)+0.95*P_ess_s(i,t)-P_ess_b(i,t)/0.95,];
end
for t=1:24
C=[C,0.1*E_max(i)<=E(i,t)<=0.9*E_max(i),];
end
C=[C,E_init(i)==0.2*E_max(i),]; %用户独立储能的初始电量约束
end
%用户自身内部独立储能的充放电功率约束,用类似于文章2.3中的Big-M法进行线性化处理
M=1E8; %这里的M是个很大的数
for i=1:3
for t=1:24
C=[C,
0<=P_ess_b(i,t)<=P_max(i),
0<=P_ess_b(i,t)<=U_ess_b(i,t)*M,
0<=P_ess_s(i,t)<=P_max(i),
0<=P_ess_s(i,t)<=U_ess_s(i,t)*M,
];
end
end