1 主要内容
该程序参考文献《多能互补微网两阶段鲁棒优化调度研究》，在考虑风光不确定集的基础上提出采用计及DR响应的多能互补微网两阶段鲁棒备用调度模型，深入研究电-气-冷-热联供型微网在充分运用新能源技术的效果，以期通过对系统鲁棒性分析实现多种能源设备的协同优化调度，有效应对多能互补微网系统的不确定性影响，同时提升系统的经济性和灵活性。
建立含电转气的多能互补微网经济调度模型，以微网系统运行费用最低作为目标，考虑多能微网系统存在的冷热电能量平衡约束以及相关设备运行约束。
通过建立多能互补微网日前-日内两阶段鲁棒优化模型，分析在系统风光不确定下的最佳经济调度方案，并考虑电网备用容量和IDR容量。日前阶段模型采用确定性调度，根据风电、光伏出力预测值，建立目标函数包括购售电成本、各类机组运行成本、备用成本以及IDR成本的调度模型，在该阶段充分考虑随机事件可能造成的影响而留有电网备用容量和IDR容量。日内阶段根据不确定集范围，调用备用资源确保多能微网系统安全稳定运行，该阶段通过对偶模型求解多能微网最恶劣运行场景，并将该结果反馈至日前阶段再次优化调度。经过两阶段的协同调度求解得到风光不确定性情况下的最佳优化结果，既保证了系统的安全性，又能得到最佳经济方案。
- 求解流程图
2 部分程序
clc;clear all; parameter; Z=[ones(8,1);zeros(40,1);ones(8,1);zeros(40,1)];%不确定变量初始值 UB=inf;%初始上界 kloop=10;%循环次数 Z0=zeros(4*T,kloop); for kl=1:kloop yalmip('clear'); Z0(:,kl)=Z;%更新不确定变量 MP;%运行第一阶段 LB=F;%更新下界 yalmip('clear'); dualtest;%运行第二阶段 UB=min(UB,LB-yita+F2);%更新上界 cul(kl)=UB-LB;%判断是否收敛 end %结果 figure; plot(cul,'LineWidth',1) xlabel('迭代次数'); ylabel('UB-LB'); grid on %典型日电转气 figure; yg2q=[x0_Pegin;x0_Pmrin;x0_Phin]'; bar(yg2q); xlabel('调度时段/h'); ylabel('功率/kW'); legend('电解槽出力','甲烷反应器出力','储氢罐输入功率'); %冷热电耦合 x0_Phfch=eta_hfch.*x0_Phout; x0_Phfce=eta_hfce.*x0_Phout; x0_Pgth=eta_gth.*x0_Pgtin; x0_Pgte=eta_gte.*x0_Pgtin; x0_Pgb=eta_gb.*x0_Pgbin; figure; yh=[x0_Phfch;x0_Pgth;x0_Pgb]'; bar(yh); hold on plot(x0_Phfce,'r-','LineWidth',1); plot(x0_Pgte,'m-','LineWidth',1); plot(eta_ac.*x0_Pacin,'b-p','LineWidth',1); xlabel('调度时段/h'); ylabel('功率/kW'); h4=legend('氢燃料电池热输出','燃气轮机热输出','燃气锅炉输出','氢燃料电池电输出','燃气轮机电输出','吸收式制冷机输出功率'); set(h4, 'Orientation', 'horizon') % set(h1, 'Box', 'off') set(h4, 'NumColumns',2); ylim([0 5500]); %电功率平衡 figure; yef=[-x0_Pbtc;-x0_Pegin;-x0_Pecin;min(x0_Pex,0)]'; b1=bar(yef,'stack'); b1(1).FaceColor = [0.1 0.5 0.9]; b1(2).FaceColor = [0.9 0.1 0.5]; b1(3).FaceColor = [0.5 0.9 0.1]; hold on yez=[x0_Pbtd;max(x0_Pex,0);x0_Pw;x0_Pv;eta_gte.*x0_Pgtin;eta_hfce.*x0_Phout]'; b2=bar(yez,'stack'); b2(1).FaceColor = [0.1 0.5 0.9]; b2(2).FaceColor = [0.5,0.3,0.5]; b2(3).FaceColor = [1,0.2,1]; b2(4).FaceColor = [1,0.6,0.1]; b2(5).FaceColor = [0.2,0.5,0.2]; b2(6).FaceColor = [0.3,0.8,0.8]; b3=plot(Pel,'r','LineWidth',1.5); h2=legend([b1(2:4),b2(1:6),b3],'电解槽','电制冷机','向电网供应功率','蓄电池','从电网吸收功率','风机','光伏','燃气轮机电功率','氢燃料电池电功率','电负荷'); set(h2, 'Orientation', 'horizon') % set(h1, 'Box', 'off') set(h2, 'NumColumns',2); xlabel('调度时段/h'); ylabel('功率/kW'); ylim([-10000 25000]); %热功率平衡 figure; yhf=-x0_Pbthc; b3=bar(yhf,'stack'); b3(1).FaceColor = [0.1 0.5 0.9]; hold on yhz=[x0_Pbthd;eta_eh.*x0_Pehin]'; b4=bar(yhz,'stack'); b4(1).FaceColor = [0.1 0.5 0.9]; b5=plot(Phl,'r','LineWidth',1.5); legend([b4(1:2),b5],'储热','余热锅炉','热负荷'); xlabel('调度时段/h'); ylabel('功率/kW'); %冷功率平衡 figure; ycf=-x0_Pbtcc; b3=bar(ycf,'stack'); b3(1).FaceColor = [0.1 0.5 0.9]; hold on ycz=[x0_Pbtcd;eta_ec.*x0_Pecin;eta_ac.*x0_Pacin]'; b4=bar(ycz,'stack'); b4(1).FaceColor = [0.1 0.5 0.9]; b5=plot(Pcl,'r','LineWidth',1.5); legend([b4(1:3),b5],'储冷','电制冷机','吸收式制冷机','冷负荷'); % legend('冷负荷'); xlabel('调度时段/h'); ylabel('功率/kW'); ylim([-1000 8000]); %负荷及风光预测值 figure; plot(Pel,'b-o','LineWidth',1.5); hold on plot(Phl,'r-*','LineWidth',1.5); plot(Pcl,'g-.','LineWidth',1.5); plot(Pwt,'c->','LineWidth',1.5); plot(Ppv,'m--','LineWidth',1.5); legend('电负荷','热负荷','冷负荷','风电预测值','光伏预测值'); xlabel('调度时段/h'); ylabel('功率/kW'); %鲁棒优化后风电出力 figure; plot(Pwt,'m-','LineWidth',1.5) hold on plot(Pwt+Pwp,'m--','LineWidth',1.5) plot(Pwt-Pwm,'m--','LineWidth',1.5) plot(x0s(73:96,1),'b-','LineWidth',1.5) legend('风电预测值','风电上限','风电下限','最恶劣场景'); xlabel('调度时段/h'); ylabel('功率/kW'); %光伏 figure; plot(Ppv,'m-','LineWidth',1.5) hold on plot(Ppv+Pvp,'m--','LineWidth',1.5) plot(Ppv-Pvm,'m--','LineWidth',1.5) plot(x0s(97:120,1),'b-','LineWidth',1.5) legend('光伏预测值','光伏上限','光伏下限','最恶劣场景'); xlabel('调度时段/h'); ylabel('功率/kW'); %电价 figure; stairs(Cph,'LineWidth',1.5) xlabel('调度时段/h'); ylabel('电价(元/(kWh))'); grid on %储氢罐容量 figure; bar(x0_Phin,'stack'); hold on bar(-x0_Phout,'stack'); ylabel('功率/kW'); yyaxis right plot(x0_Wh,'r--','LineWidth',1.5) ylabel('储氢罐容量/kWh'); legend('储存氢气','释放氢气','储氢容量') xlabel('调度时段/h'); grid on %储冷容量 figure; bar(x0_Pbtcc,'stack'); hold on bar(-x0_Pbtcd,'stack'); ylabel('功率/kW'); yyaxis right plot(x0_Wbtc,'r--','LineWidth',1.5) ylabel('储冷罐容量/kWh'); legend('储存冷能','释放冷能','储冷容量') xlabel('调度时段/h'); grid on %储热容量 figure; bar(x0_Pbthc,'stack'); hold on bar(-x0_Pbthd,'stack'); ylabel('功率/kW'); yyaxis right plot(x0_Wbth,'r--','LineWidth',1.5) ylabel('储热罐容量/kWh'); legend('储存热能','释放热能','储热容量') xlabel('调度时段/h'); grid on %蓄电池容量 figure; bar(x0_Pbtc,'stack'); hold on bar(-x0_Pbtd,'stack'); ylabel('功率/kW'); yyaxis right plot(x0_Wbt,'r--','LineWidth',1.5) ylabel('蓄电池容量/kWh'); legend('储存电能','释放电能','蓄电池容量') xlabel('调度时段/h'); grid on
3 程序结果
含保守度的迭代末次风电最恶劣场景图，在两阶段鲁棒优化过程中，随着主子问题迭代，二阶段可再生能源的最恶劣场景也会随之变化，将这些最恶劣场景叠加到一阶段进行综合优化求解，因此遇到有最恶劣场景图的可以注意一下，程序内部是不是实现了叠加循环。
