考虑风光燃储的综合能源系统鲁棒优化调度

clc; clear; yalmip; %常数 T=24;%时间24小时 a=0.57;%燃气轮机成本系数 PGmax=30;%微型燃气轮机的最大输出功率 PGmin=5;%微型燃气轮机的最小输出功率 Ks=0.08;%折算后的单位充放电成本 y=0.95;%储能单元的充放电效率 PSmax=10;%储能允许的最大充放电功率 ES0=25;% ESmax=40;%储能为50 ESmin=10;% DDR=335.53;% KDR=0.36;% nlo=0.58;% nwt=0.4;% npv=0.4;% %决策变量 PG=sdpvar(1,T,'full');%燃气轮机出力变量 Pdis=sdpvar(1,T,'full');%放电 Pcha=sdpvar(1,T,'full');%充电 PDR=sdpvar(1,T,'full');%需求响应调度负荷 PDR1=sdpvar(1,T); PDR2=sdpvar(1,T); rlo=sdpvar(1,T,'full');%弃负荷量 rwt=sdpvar(1,T,'full');%弃风量 rpv=sdpvar(1,T,'full');%弃光量 PE=sdpvar(1,24,'full');%蓄电池出力 Temp_cha=binvar(1,24,'full'); %充电标志 Temp_dis=binvar(1,24,'full'); %放电标志 Temp_static=binvar(1,24,'full'); %电池静置标志 %风光负荷预测 xx=1:T; PV=[0 0 0 0 0 0 2.86 14.21 18.53 29.5 38.59 44.58 46.86 44.55 39.51 31.33 20.84 15.26 4.53 0.81 0 0 0 0 ]; WT=[61.92 59.45 56.65 52.53 48.82 43.11 39.35 37.77 36.44 33.75 32.26 33.11 34.02 32.62 33.51 37.38 40.36 43.01 45.48 45.86 49.8 54.29 55.87 59.87 ]; PL=[27.73 24.37 19.82 18.27 15.12 14.61 25.81 41.74 57.32 69.85 75.01 76.71 64 67.8 66.48 69.16 70.27 70.52 71.17 68.61 71.74 57.47 46.04 29.28 ]; PDRmax=[34.19 35.08 36.83 34.26 33.7 28.5 16.4 10.24 -2.35 -6.6 -4.16 0.98 16.88 9.37 6.54 -0.45 -9.07 -12.25 -21.16 -21.94 -21.94 -3.18 9.83 30.59 ]; PDRmin=[1.156 0.442 0.367 0.286 0.467 0.721 0.942 1.154 1.47 1.912 2.141 2.381 1.752 1.143 1.103 1.575 1.614 1.796 2.225 2.738 2.772 1.994 1.364 0.616 ]; PDRx=[5.78 4.42 3.67 2.86 4.67 7.21 9.42 11.54 14.7 19.12 21.41 23.81 17.52 11.43 11.03 15.75 16.14 17.96 22.25 27.38 27.72 19.94 13.64 6.16 ];%PD* bij=[12.061 11.3545 10.4795 9.7065 8.835 7.9275 8.9125 11.971 13.9775 16.4725 18.1285 19.3245 18.532 18.3555 17.601 17.2225 16.207 15.7925 14.6185 13.8615 14.644 13.8905 12.9845 11.9085 ];%最大预测误差幅值和数据 rlomax=[2.773 2.437 1.982 1.827 1.512 1.461 2.581 4.174 5.732 6.985 7.501 7.671 6.4 6.78 6.648 6.916 7.027 7.052 7.117 6.861 7.174 5.747 4.604 2.928 ];%弃负荷最大 rwtmax=[6.192 5.945 5.665 5.253 4.882 4.311 3.935 3.777 3.644 3.375 3.226 3.311 3.402 3.262 3.351 3.738 4.036 4.301 4.548 4.586 4.98 5.429 5.587 5.987 ];%弃风最大 rpvmax=[0 0 0 0 0 0 0.286 1.421 1.853 2.95 3.859 4.458 4.686 4.455 3.951 3.133 2.084 1.526 0.453 0.081 0 0 0 0 ];%弃光最大 st =[];%约束条件 %燃气轮机出力约束 for t=1:T st = [st,PGmin<=PG(1,t)<=PGmax]; end %储能出力约束 for t=1:T st = [st, -10<=PE(1,t)<=10,0<=Pcha(1,t)<=PSmax,-PSmax<=Pdis(1,t)<=0];%电池充放电约束 st = [st, implies(Temp_cha(1,t),[PE(1,t)>=0,Pcha(1,t)==PE(1,t),Pdis(1,t)==0])];%充电情况约束 st = [st, implies(Temp_dis(1,t),[PE(1,t)<=0,Pdis(1,t)==PE(1,t),Pcha(1,t)==0])];%放电情况约束 st = [st, implies(Temp_static(1,t),[PE(1,t)==0,Pdis(1,t)==0,Pcha(1,t)==0])];%静置情况约束 st = [st,Temp_cha(1,t)+Temp_dis(1,t)+Temp_static(1,t)==1]; st=[st,-15<=(1/y)*sum(Pdis(1,1:t))+y*sum(Pcha(1,1:t))<=25] ;%SOC约束，电池容量300kwh，初始S0C为0.4，0.3<=SOC<=0.95 end st=[st,(1/y)*sum(Pdis)+y*sum(Pcha)==0] ;%ST=S0，始末E(t)相等约束 %需求响应约束 for t=1:T st=[st,sum(PDR)==DDR]; end for t=1:T st=[st,PDRmin(1,t)<=PDR(1,t)<=PDRmax(1,t)+PG(1,t)]; end for t=1:T st=[st,PDR(1,t)-PDRx(1,t)+PDR1(1,t)-PDR2(1,t)==0]; end for t=1:T st=[st,PDR1(1,t)>=0]; end for t=1:T st=[st,PDR2(1,t)>=0]; end %需求响应约束 for t=1:T st=[st,sum(PDR)==DDR]; end for t=1:T st=[st,PDRmin(1,t)<=PDR(1,t)<=PDRmax(1,t)+PG(1,t)]; end for t=1:T st=[st,PDR(1,t)-PDRx(1,t)+PDR1(1,t)-PDR2(1,t)==0]; end for t=1:T st=[st,PDR1(1,t)>=0]; end for t=1:T st=[st,PDR2(1,t)>=0]; end %弃风弃光弃负荷约束 for t=1:T st=[st,0<=rlo(1,t)<=rlomax(1,t)]; end for t=1:T st=[st,0<=rwt(1,t)<=rwtmax(1,t)]; end for t=1:T st=[st,0<=rpv(1,t)<=rpvmax(1,t)]; end %功率平衡约束 for t=1:T st=[st,PL(1,t)-rlo(1,t)==WT(1,t)+PV(1,t)+PG(1,t)-PDR(1,t)-rwt(1,t)-rpv(1,t)-PE(1,t)]; end %目标函数 z=0; for t=1:T z=z+(a*PG(1,t)+nlo*rlo(1,t)+KDR*(PDR1(1,t)+PDR2(1,T))+nwt*rwt(1,t)+npv*rpv(1,t))-Ks*Pdis(1,t)/y+Ks*Pcha(1,t)*y; end %配置求解器 ops = sdpsettings('solver','cplex');%求解 optimize(st,z,ops) z=value(z) PG=value(PG); rpv=value(rpv); rlo=value(rlo); rwt=value(rwt); PDR=value(PDR); PDR1=value(PDR1); PDR2=value(PDR2); Pdis=value(Pdis); Pcha=value(Pcha); %% 电功率平衡 subplot(1,1,1) bar([PG',-Pdis',-Pcha'],'stack') %阶梯图 legend('PG','Pdis','Pcha'); %在坐标轴上添加图例 set(0,'ShowHiddenHandles','On') set(gcf,'menubar','figure') xlabel('时刻（t）'); ylabel('功率（kW）');