EI论文复现：考虑储能电站分布特性的微电网多时间尺度协同优化程序代码！

本程序参考EI论文《考虑特性分布的储能电站接入的电网多时间尺度源储荷协调调度策略》，文章综合储能电站、负荷侧各类需求响应资源的多时间尺度特性制定日前调度计划，并通过日内滚动与实时修正实现对预测数据更高精度的保证。程序中算例丰富、注释清晰、干货满满，创新性很高！下面对文章和程序做简要介绍！

创新点：

1. 源-荷-储协同调度；

2. 多时间尺度：

3. 需求响应；

为提高区域电网新能源消纳率，使电力系统运行成本最小，提高电力系统供电可靠性。提出考虑特性分布的储能电站接入的电网多时间尺度“源–储–荷”协调调度策略。综合储能电站、 负荷侧各类需求响应资源的多时间尺度特性制定日前调度计划，并通过日内滚动与实时修正实现对预测数据更高精度的保证。以系统运行成本、弃风惩罚成本、失负荷惩罚成本为目标函数，建立多时间尺度调度模型。最后通过MATLAB平台调用商用软件包CPLEX完成混合整数规划优化，结果表明该调度策略可以达到系统运行成本最优，降低弃风率的效果。有效提高了“源–储–荷”系统的新能源消纳率。

程序结果：

上图为系统成本运行结果图，日内小时调度结果​！

部分程序：

%目标函数
%常规机组运行成本a2*p+bp+c分段线性化
tn=T;
Pmax=[mp(1,2),mp(2,2),mp(3,2),mp(4,2),mp(5,2),mp(6,2)];
Pmin=[mp(1,3),mp(2,3),mp(3,3),mp(4,3),mp(5,3),mp(6,3)];
gn=5;%分段函数线性化，下同
x_pf=sdpvar(6, tn,'full');
gw1=sdpvar(gn+1,tn,'full');
gw2=sdpvar(gn+1,tn,'full');
gw3=sdpvar(gn+1,tn,'full');
gw4=sdpvar(gn+1,tn,'full');
gw5=sdpvar(gn+1,tn,'full');
gw6=sdpvar(gn+1,tn,'full');
gz1=binvar(gn, tn,'full');
gz2=binvar(gn, tn,'full');
gz3=binvar(gn, tn,'full');
gz4=binvar(gn, tn,'full');
gz5=binvar(gn, tn,'full');
gz6=binvar(gn, tn,'full');
gl1=(Pmax-Pmin)./gn;
gl2=zeros(6,gn+1);
for i=1:6
gl2(i,:)=Pmin(i):gl1(i):Pmax(i);
end
con = [con, x_pf(1,:)==gl2(1,:).^2*gw1];
con = [con, x_pf(2,:)==gl2(2,:).^2*gw2];
con = [con, x_pf(3,:)==gl2(3,:).^2*gw3];
con = [con, x_pf(4,:)==gl2(4,:).^2*gw4];
con = [con, x_pf(5,:)==gl2(5,:).^2*gw5];
con = [con, x_pf(6,:)==gl2(6,:).^2*gw6];
con = [con, gw1(1,:)<=gz1(1,:)];
for i=2:gn
    con = [con, gw1(i,:)<=gz1(i-1,:)+gz1(i,:)];
end
con = [con, gw1(gn+1,:)<=gz1(gn,:)];
con = [con, sum(gz1)==ones(1,tn)];
con = [con, gw2(1,:)<=gz2(1,:)];
for i=2:gn
    con = [con, gw2(i,:)<=gz2(i-1,:)+gz2(i,:)];
end
con = [con, gw2(gn+1,:)<=gz2(gn,:)];
con = [con, sum(gz2)==ones(1,tn)];
con = [con, gw3(1,:)<=gz3(1,:)];
for i=2:gn
    con = [con, gw3(i,:)<=gz3(i-1,:)+gz3(i,:)];
end
con = [con, gw3(gn+1,:)<=gz3(gn,:)];
con = [con, sum(gz3)==ones(1,tn)];
con = [con, gw4(1,:)<=gz4(1,:)];
for i=2:gn
    con = [con, gw4(i,:)<=gz4(i-1,:)+gz4(i,:)];
end
con = [con, gw4(gn+1,:)<=gz4(gn,:)];
con = [con, sum(gz4)==ones(1,tn)];
con = [con, gw5(1,:)<=gz5(1,:)];
for i=2:gn
    con = [con, gw5(i,:)<=gz5(i-1,:)+gz5(i,:)];
end
con = [con, gw5(gn+1,:)<=gz5(gn,:)];
con = [con, sum(gz5)==ones(1,tn)];
con = [con, gw6(1,:)<=gz6(1,:)];
for i=2:gn
    con = [con, gw6(i,:)<=gz6(i-1,:)+gz6(i,:)];
end
con = [con, gw6(gn+1,:)<=gz6(gn,:)];
con = [con, sum(gz6)==ones(1,tn)];
con = [con, Pgone(1,:)==gl2(1,:)*gw1];
con = [con, Pgone(2,:)==gl2(2,:)*gw2];
con = [con, Pgone(3,:)==gl2(3,:)*gw3];
con = [con, Pgone(4,:)==gl2(4,:)*gw4];
con = [con, Pgone(5,:)==gl2(5,:)*gw5];
con = [con, Pgone(6,:)==gl2(6,:)*gw6];
con = [con, gw1>=0];
con = [con, gw2>=0];
con = [con, gw3>=0];
con = [con, gw4>=0];
con = [con, gw5>=0];
con = [con, gw6>=0];
%场景2线性化
x_pfz=sdpvar(6, tn,'full');
gw1z=sdpvar(gn+1,tn,'full');
gw2z=sdpvar(gn+1,tn,'full');
gw3z=sdpvar(gn+1,tn,'full');
gw4z=sdpvar(gn+1,tn,'full');
gw5z=sdpvar(gn+1,tn,'full');
gw6z=sdpvar(gn+1,tn,'full');
gz1z=binvar(gn, tn,'full');
gz2z=binvar(gn, tn,'full');
gz3z=binvar(gn, tn,'full');
gz4z=binvar(gn, tn,'full');
gz5z=binvar(gn, tn,'full');
gz6z=binvar(gn, tn,'full');
gl1z=(Pmax-Pmin)./gn;
gl2z=zeros(6,gn+1);
for i=1:6
gl2z(i,:)=Pmin(i):gl1z(i):Pmax(i);
end
con = [con, x_pfz(1,:)==gl2z(1,:).^2*gw1z];
con = [con, x_pfz(2,:)==gl2z(2,:).^2*gw2z];
con = [con, x_pfz(3,:)==gl2z(3,:).^2*gw3z];
con = [con, x_pfz(4,:)==gl2z(4,:).^2*gw4z];
con = [con, x_pfz(5,:)==gl2z(5,:).^2*gw5z];
con = [con, x_pfz(6,:)==gl2z(6,:).^2*gw6z];
con = [con, gw1z(1,:)<=gz1z(1,:)];
for i=2:gn
    con = [con, gw1z(i,:)<=gz1z(i-1,:)+gz1z(i,:)];
end
con = [con, gw1z(gn+1,:)<=gz1z(gn,:)];
con = [con, sum(gz1z)==ones(1,tn)];
con = [con, gw2z(1,:)<=gz2z(1,:)];
for i=2:gn
    con = [con, gw2z(i,:)<=gz2z(i-1,:)+gz2z(i,:)];
end
con = [con, gw2z(gn+1,:)<=gz2z(gn,:)];
con = [con, sum(gz2z)==ones(1,tn)];
con = [con, gw3z(1,:)<=gz3z(1,:)];
for i=2:gn
    con = [con, gw3z(i,:)<=gz3z(i-1,:)+gz3z(i,:)];
end
con = [con, gw3z(gn+1,:)<=gz3z(gn,:)];
con = [con, sum(gz3z)==ones(1,tn)];
con = [con, gw4z(1,:)<=gz4z(1,:)];
for i=2:gn
    con = [con, gw4z(i,:)<=gz4z(i-1,:)+gz4z(i,:)];
end
con = [con, gw4z(gn+1,:)<=gz4z(gn,:)];
con = [con, sum(gz4z)==ones(1,tn)];
con = [con, gw5z(1,:)<=gz5z(1,:)];
for i=2:gn
    con = [con, gw5z(i,:)<=gz5z(i-1,:)+gz5z(i,:)];
end
con = [con, gw5z(gn+1,:)<=gz5z(gn,:)];
con = [con, sum(gz5z)==ones(1,tn)];
con = [con, gw6z(1,:)<=gz6z(1,:)];
for i=2:gn
    con = [con, gw6z(i,:)<=gz6z(i-1,:)+gz6z(i,:)];
end
con = [con, gw6z(gn+1,:)<=gz6z(gn,:)];
con = [con, sum(gz6z)==ones(1,tn)];
con = [con, Pgtwo(1,:)==gl2z(1,:)*gw1z];
con = [con, Pgtwo(2,:)==gl2z(2,:)*gw2z];
con = [con, Pgtwo(3,:)==gl2z(3,:)*gw3z];
con = [con, Pgtwo(4,:)==gl2z(4,:)*gw4z];
con = [con, Pgtwo(5,:)==gl2z(5,:)*gw5z];
con = [con, Pgtwo(6,:)==gl2z(6,:)*gw6z];
con = [con, gw1z>=0];
con = [con, gw2z>=0];
con = [con, gw3z>=0];
con = [con, gw4z>=0];
con = [con, gw5z>=0];
con = [con, gw6z>=0];    
for gt=1:24
    for gi=1:6
        fpg1(gi,gt)=mp(gi,4).*x_pf(gi,gt)+mp(gi,5).*Pgone(gi,gt)+mp(gi,6);
        fpg2(gi,gt)=mp(gi,4).*x_pfz(gi,gt)+mp(gi,5).*Pgtwo(gi,gt)+mp(gi,6);
    end
end
 yyone=binvar(6,23,'full');
 yytwo=binvar(6,23,'full');
 for i=1:23
     for j=1:6
     con=[con,yyone(j,i)<=ugone(j,i+1)];
     con=[con,yyone(j,i)<=1-ugone(j,i)];
     con=[con,yyone(j,i)>=ugone(j,i+1)-ugone(j,i)];
     con=[con,yytwo(j,i)<=ugtwo(j,i+1)];
     con=[con,yytwo(j,i)<=1-ugtwo(j,i)];
     con=[con,yytwo(j,i)>=ugtwo(j,i+1)-ugtwo(j,i)];
     end
 end
  yygone=binvar(1,23,'full');
  yygtwo=binvar(1,23,'full');
 for i=1:23
     con=[con,yygone(1,i)<=upwone(1,i+1)];
     con=[con,yygone(1,i)<=1-upwone(1,i)];
     con=[con,yygone(1,i)>=upwone(1,i+1)-upwone(1,i)];
     con=[con,yygtwo(1,i)<=upwtwo(1,i+1)];
     con=[con,yygtwo(1,i)<=1-upwtwo(1,i)];
     con=[con,yygtwo(1,i)>=upwtwo(1,i+1)-upwtwo(1,i)];
 end
 con=[con,0<=plossone<=pload1,0<=plosstwo<=pload1];
fg=sum(sum(ps(1).*fpg1+ps(2).*fpg2))+0.9*sum(sum(ps(1).*yyone+ps(2).*yytwo));
fess=sum(ps(1)*(0.53*abs(Perssone)+0.12*abs(Perssone))+ps(2)*(0.53*abs(Persstwo)+0.12*abs(Persstwo)))+sum(0.012.*(uersstwo+uersstwo));
fdg=ps(1)*sum(0.02.*Pwone)+ps(2)*(sum(0.02.*Pwtwo)+0.9*sum(ps(1).*yygone+ps(2).*yygtwo))+ps(1)*5.2*sum(pw1-Pwone)+ps(2)*5.2*sum(pw2-Pwtwo);
fload=ps(1)*(sum(100.*abs(Pidraone)+125.*abs(Pidrbone))+15.9.*sum(sum(plossone)))+ps(2)*(sum(100.*abs(Pidratwo)+125.*abs(Pidrbtwo))+15.9.*sum(sum(plosstwo)));
f1=fg+fess+fdg+fload;
ops=sdpsettings('solver','cplex');
result=optimize(con,f1,ops);

部分图片来源于网络，侵权联系删除！

欢迎感兴趣的小伙伴关注，小编会继续推送更有质量的学习资料、文章和程序代码！