目录
1 主要内容
该程序复现《考虑区域多能源系统集群协同优化的联合需求侧响应模型》文献模型,程序的核心是对多个区域级多能源系统互联系统进行多目标优化,并且考虑联合需求侧响应,以多个区域多能源系统运行总成本最小、碳排放最小为目标,建立多区域电气热(冷)互联系统多目标优化模型,和原文的区别是:多目标求解原文献用的是NSGA_Ⅱ算法,但是程序采用的是混合整数规划算法,即直接采用yalmip求解器进行求解,程序步骤清晰,不仅给出了多目标求解代码,同时给出了单目标求解代码,采用matlab+yalmip(cplex或gurobi)进行求解,必要注释清晰,方便学习!
- 区域多能源系统用户激励模型
该部分用户激励模型和常规的需求侧响应区别不大,采用可削减、可转移和可替代负荷作为需求响应变量,并且该三类负荷满足调节比例限制,该部分成本在程序中以约束的形式表达,但是实际上尽量选取变量直接计算的方式更为妥当。对应的程序代码如下:
C=[C, 0<=Ekexuejian1,Ekexuejian1<=0.1*E_load1D ];%需求响应 公式(2) C=[C, 0<=Ekexuejian2,Ekexuejian2<=0.1*E_load2D ]; C=[C, 0<=Ekexuejian3,Ekexuejian3<=0.1*E_load3D ]; C=[C, -0.1*E_load1D<=Ekepingyi1,Ekepingyi1<=0.1*E_load1D ]; C=[C, -0.1*E_load2D<=Ekepingyi2,Ekepingyi2<=0.1*E_load2D ]; C=[C, -0.1*E_load3D<=Ekepingyi3,Ekepingyi3<=0.1*E_load3D ]; C=[C, sum(Ekepingyi1)==0 ]; C=[C, sum(Ekepingyi2)==0 ]; C=[C, sum(Ekepingyi3)==0 ]; C=[C, 0<=Eketidai1,Eketidai1<=0.1*E_load1D ]; C=[C, 0<=Eketidai2,Eketidai2<=0.1*E_load2D ]; C=[C, 0<=Eketidai3,Eketidai3<=0.1*E_load3D ]; C=[C, Eketidai1==Gketidai1*10 ]; %10是单位m3体积的天然气热值36MJ C=[C, Eketidai2==Gketidai2*10 ]; %10是单位m3体积的天然气热值36MJ C=[C, Eketidai3==Gketidai3*10 ]; %10是单位m3体积的天然气热值36MJ %% 需求侧响应补偿成本 公式(1) F_IDSR = sdpvar(1,1); C=[C,F_IDSR==0.432*sum(abs(Ekexuejian1)+abs(Ekexuejian2)+abs(Ekexuejian3))+0.060*sum(abs(Ekepingyi1)+abs(Ekepingyi2)+abs(Ekepingyi3))/2+0.120*sum(abs(Eketidai1)+abs(Eketidai2)+abs(Eketidai3)) ];
- 目标函数
该程序目标函数有两个,分别是总运行成本和总排放量,程序采用的是yalmip直接求解帕累托解集的方式,之前求解帕累托解集均采用智能算法,这种方式给我们提供了一种思路:采用两个目标变动加权的形式寻找模型的前沿面。
2 部分代码
%空调 EKongTiao1=sdpvar(1,24); EKongTiao1min= 0; EKongTiao1max= 1000*200; HKongTiao1=sdpvar(1,24); LKongTiao1=sdpvar(1,24); HBKongTiao1=binvar(1,24); LBKongTiao1=binvar(1,24); %电制冷 EDianZhileng1=sdpvar(1,24); EDianZhileng1min= 0; EDianZhileng1max= 1000*500; LDianZhileng1=sdpvar(1,24); %吸热制冷 HReZhileng1=sdpvar(1,24); LReZhileng1min= 0; LReZhileng1max= 1000*200; LReZhileng1=sdpvar(1,24); %P2G EP2G1=sdpvar(1,24); EP2G1min= 0; EP2G1max= 1000*100; GP2G1=sdpvar(1,24); %需求响应 Ekexuejian1= sdpvar(1,24); Ekepingyi1= sdpvar(1,24); Eketidai1= sdpvar(1,24); Gketidai1= sdpvar(1,24); %% 区域②的负荷 冬季 E_load2D = kL*1000*[375,300,250,250,250,250,375,500,300,300,300,300,400,300,300,300,300,500,700,750,750,700,660,660]; G_load2D = kL*1000*[50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,100,100,50,50,50,50 ]; H_load2D = kL*1000*[450,500,500,525,550,575,600,625,375,325,250,200,150,125,125,125,125,625,500,440,375,1000,750,650 ]; L_load2D = kL*1000*[65,125,200,250,300,375,250,125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ]; % 区域②的负荷 夏季 E_load2X = kL*1000*[250,250,250,250,250,250,375,375,250,250,250,250,325,250,250,250,250,325,500,875,875,825,750,625]; G_load2X = kL*1000*[50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,50,100,100,50,50,50,50 ]; H_load2X = kL*1000*[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 ]; L_load2X = kL*1000*[60,60,60,60,60,60,60,125,150,175,200,225,250,375,425,300,250,250,250,750,650,425,375,325]; % 区域②变量初始化 %电网交互 Egrid2 = sdpvar(1,24); Egrid2min = 0; Egrid2max = 0; %热网交互 Hgrid2 = sdpvar(1,24); Hgrid2min= 0; Hgrid2max= 0; %天然气网交互 Ggrid2 = sdpvar(1,24); Ggrid2min= 0; Ggrid2max= 1000*2000; %燃气轮机 GMT2=sdpvar(1,24); EMT2=sdpvar(1,24); EMT2max= 1000*1200; EMT2min= 1000*120; %燃气锅炉 GGB2=sdpvar(1,24); HGB2=sdpvar(1,24); HGB2max= 1000*200; HGB2min= 0; %空调 EKongTiao2=sdpvar(1,24); EKongTiao2min= 0; EKongTiao2max= 1000*1000; HKongTiao2=sdpvar(1,24); LKongTiao2=sdpvar(1,24); HBKongTiao2=binvar(1,24); LBKongTiao2=binvar(1,24); %电制冷 EDianZhileng2=sdpvar(1,24); EDianZhileng2min= 0; EDianZhileng2max= 1000*500; LDianZhileng2=sdpvar(1,24); %吸热制冷 HReZhileng2=sdpvar(1,24); LReZhileng2min= 0; LReZhileng2max= 1000*200; LReZhileng2=sdpvar(1,24); %P2G EP2G2=sdpvar(1,24); EP2G2min= 0; EP2G2max= 1000*100; GP2G2=sdpvar(1,24); %需求响应 Ekexuejian2= sdpvar(1,24); Ekepingyi2= sdpvar(1,24); Eketidai2= sdpvar(1,24); Gketidai2= sdpvar(1,24); %% 区域③的负荷 冬季 E_load3D = kL*1000*[300,300,300,300,300,300,300,300,300,875,875,875,875,875,875,875,875,875,300,300,300,300,300,300]; G_load3D = kL*1000*[150,150,150,150,150,150,200,250,375,365,355,345,355,365,375,375,345,200,200,200,150,150,150,150]; H_load3D = kL*1000*[375,375,375,375,375,375,375,375,375,625,625,625,625,625,625,625,625,625,375,375,375,375,375,375]; L_load3D = kL*1000*[125,125,125,125,125,125,125,125,375,375,375,375,375,375,375,375,375,375,125,125,125,125,125,125 ]; % 区域③的负荷 夏季 E_load3X = kL*1000*[375,375,375,375,375,375,375,375,375,875,875,875,875,875,875,875,875,875,375,375,375,375,375,375]; G_load3X = kL*1000*[150,150,150,150,150,150,200,350,375,365,355,345,355,365,375,375,345,200,200,200,150,150,150,150 ]; H_load3X = kL*1000*[375,375,375,375,375,375,375,375,375,500,500,500,500,500,500,500,500,500,375,375,375,375,375,375 ]; L_load3X = kL*1000*[200,200,200,200,200,200,200,200,200,625,625,625,625,625,625,625,625,625,200,200,200,200,200,200]; % 区域③变量初始化 %电网交互 Egrid3 = sdpvar(1,24); Egrid3min = 0; Egrid3max = 1000*2000; %热网交互 Hgrid3 = sdpvar(1,24); Hgrid3min= 0; Hgrid3max= 0; %天然气网交互 Ggrid3 = sdpvar(1,24); Ggrid3min= 0; Ggrid3max= 0; %燃气轮机 GMT3=sdpvar(1,24); EMT3=sdpvar(1,24); EMT3max= 1000*1100; EMT3min= 1000*110; %燃气锅炉 GGB3=sdpvar(1,24); HGB3=sdpvar(1,24); HGB3max= 1000*500; HGB3min= 0; %空调 EKongTiao3=sdpvar(1,24); EKongTiao3min= 0; EKongTiao3max= 1000*1000; HKongTiao3=sdpvar(1,24); LKongTiao3=sdpvar(1,24); HBKongTiao3=binvar(1,24); LBKongTiao3=binvar(1,24); %电制冷 EDianZhileng3=sdpvar(1,24); EDianZhileng3min= 0; EDianZhileng3max= 1000*500; LDianZhileng3=sdpvar(1,24); %吸热制冷 HReZhileng3=sdpvar(1,24); LReZhileng3min= 0; LReZhileng3max= 1000*200; LReZhileng3=sdpvar(1,24); %P2G EP2G3=sdpvar(1,24); EP2G3min= 0; EP2G3max= 1000*100; GP2G3=sdpvar(1,24); %需求响应 Ekexuejian3= sdpvar(1,24); Ekepingyi3= sdpvar(1,24); Eketidai3= sdpvar(1,24); Gketidai3= sdpvar(1,24); %% 区域①②③的设备运行约束建模 C=[]; C=[C,Egrid1min<=Egrid1,Egrid1<=Egrid1max];%电网交互约束 C=[C,Egrid2min<=Egrid2,Egrid2<=Egrid2max]; C=[C,Egrid3min<=Egrid3,Egrid3<=Egrid3max]; C=[C,Hgrid1min<=Hgrid1,Hgrid1<=Hgrid1max];%热网交互约束 C=[C,Hgrid2min<=Hgrid2,Hgrid2<=Hgrid2max]; C=[C,Hgrid3min<=Hgrid3,Hgrid3<=Hgrid3max]; C=[C,Ggrid1min<=Ggrid1,Ggrid1<=Ggrid1max];%气网交互约束 C=[C,Ggrid2min<=Ggrid2,Ggrid2<=Ggrid2max]; C=[C,Ggrid3min<=Ggrid3,Ggrid3<=Ggrid3max]; C=[C,EMT1min<=EMT1,EMT1<=EMT1max];%燃气轮机约束 C=[C,EMT2min<=EMT2,EMT2<=EMT2max]; C=[C,EMT3min<=EMT3,EMT3<=EMT3max]; C=[C,EMT1==GMT1*10*0.35 ]; %nMT=0.35 C=[C,EMT2==GMT2*10*0.35 ]; C=[C,EMT3==GMT3*10*0.35 ]; C=[C,HGB1min<=HGB1,HGB1<=HGB1max];%燃气锅炉约束 C=[C,HGB2min<=HGB2,HGB2<=HGB2max]; C=[C,HGB3min<=HGB3,HGB3<=HGB3max]; C=[C,HGB1==GGB1*10*0.8 ]; %nGB=0.8 C=[C,HGB2==GGB2*10*0.8 ]; C=[C,HGB3==GGB3*10*0.8 ]; C=[C,EKongTiao1min<=EKongTiao1,EKongTiao1<=EKongTiao1max];%空调约束 C=[C,EKongTiao2min<=EKongTiao2,EKongTiao2<=EKongTiao2max]; C=[C,EKongTiao3min<=EKongTiao3,EKongTiao3<=EKongTiao3max]; C=[C,HBKongTiao1+LBKongTiao1<=1 ]; C=[C,HBKongTiao2+LBKongTiao2<=1 ]; C=[C,HBKongTiao3+LBKongTiao3<=1 ]; M=1e7; C=[C,0<=HKongTiao1,HKongTiao1<=HBKongTiao1*M ]; C=[C,0<=HKongTiao2,HKongTiao2<=HBKongTiao2*M ]; C=[C,0<=HKongTiao3,HKongTiao3<=HBKongTiao3*M ]; C=[C,0<=LKongTiao1,LKongTiao1<=LBKongTiao1*M ]; C=[C,0<=LKongTiao2,LKongTiao2<=LBKongTiao2*M ]; C=[C,0<=LKongTiao3,LKongTiao3<=LBKongTiao3*M ]; C=[C,EKongTiao1==(HKongTiao1+LKongTiao1)*2.6 ]; %电空调转热转冷效率1/2.6 C=[C,EKongTiao2==(HKongTiao2+LKongTiao2)*2.6 ]; C=[C,EKongTiao3==(HKongTiao3+LKongTiao3)*2.6 ]; C=[C,EDianZhileng1min<=EDianZhileng1,EDianZhileng1<=EDianZhileng1max];%电制冷约束 C=[C,EDianZhileng2min<=EDianZhileng2,EDianZhileng2<=EDianZhileng2max]; C=[C,EDianZhileng3min<=EDianZhileng3,EDianZhileng3<=EDianZhileng3max]; C=[C,EDianZhileng1==LDianZhileng1*4 ]; %电制冷1/4 C=[C,EDianZhileng2==LDianZhileng2*4 ]; C=[C,EDianZhileng3==LDianZhileng3*4 ]; C=[C,LReZhileng1min<=LReZhileng1,LReZhileng1<=LReZhileng1max];%吸热制冷约束 C=[C,LReZhileng2min<=LReZhileng2,LReZhileng2<=LReZhileng2max]; C=[C,LReZhileng3min<=LReZhileng3,LReZhileng3<=LReZhileng3max]; C=[C,HReZhileng1==LReZhileng1/0.7 ]; %吸热制冷0.7 C=[C,HReZhileng2==LReZhileng2/0.7 ]; C=[C,HReZhileng3==LReZhileng3/0.7 ]; C=[C,EP2G1min<=EP2G1,EP2G1<=EP2G1max];%P2G约束 C=[C,EP2G2min<=EP2G2,EP2G2<=EP2G2max]; C=[C,EP2G3min<=EP2G3,EP2G3<=EP2G3max]; C=[C,EP2G1==GP2G1*10/0.5 ]; %nP2G=0.5 C=[C,EP2G2==GP2G2*10/0.5 ]; C=[C,EP2G3==GP2G3*10/0.5 ]; C=[C, 0<=Ekexuejian1,Ekexuejian1<=0.1*E_load1D ];%需求响应 公式(2) C=[C, 0<=Ekexuejian2,Ekexuejian2<=0.1*E_load2D ]; C=[C, 0<=Ekexuejian3,Ekexuejian3<=0.1*E_load3D ]; C=[C, -0.1*E_load1D<=Ekepingyi1,Ekepingyi1<=0.1*E_load1D ]; C=[C, -0.1*E_load2D<=Ekepingyi2,Ekepingyi2<=0.1*E_load2D ]; C=[C, -0.1*E_load3D<=Ekepingyi3,Ekepingyi3<=0.1*E_load3D ]; C=[C, sum(Ekepingyi1)==0 ]; C=[C, sum(Ekepingyi2)==0 ]; C=[C, sum(Ekepingyi3)==0 ]; C=[C, 0<=Eketidai1,Eketidai1<=0.1*E_load1D ]; C=[C, 0<=Eketidai2,Eketidai2<=0.1*E_load2D ]; C=[C, 0<=Eketidai3,Eketidai3<=0.1*E_load3D ]; C=[C, Eketidai1==Gketidai1*10 ]; %10是单位m3体积的天然气热值36MJ C=[C, Eketidai2==Gketidai2*10 ]; %10是单位m3体积的天然气热值36MJ C=[C, Eketidai3==Gketidai3*10 ]; %10是单位m3体积的天然气热值36MJ %% 区域①②③平衡约束 %电平衡 C=[C, Egrid1+Egrid2+Egrid3+EMT1+EMT2+EMT3-EKongTiao1-EKongTiao2-EKongTiao3-EP2G1-EP2G2-EP2G3-EDianZhileng1-EDianZhileng2-EDianZhileng3==E_load1D+E_load2D+E_load3D+Ekepingyi1+Ekepingyi2+Ekepingyi3-Ekexuejian1-Ekexuejian2-Ekexuejian3-Eketidai1-Eketidai2-Eketidai3 ]; %热平衡 C=[C,Hgrid1+Hgrid2+Hgrid3+HKongTiao1+HKongTiao2+HKongTiao3+HGB1+HGB2+HGB3-HReZhileng1-HReZhileng2-HReZhileng3==H_load1D+H_load2D+H_load3D ]; %冷平衡 C=[C,LDianZhileng1+LDianZhileng2+LDianZhileng3+LReZhileng1+LReZhileng2+LReZhileng3+LKongTiao1+LKongTiao2+LKongTiao3==L_load1D+L_load2D+L_load3D ]; %气平衡 C=[C,Ggrid1+Ggrid2+Ggrid3+GP2G1+GP2G2+GP2G3 == G_load1D+G_load2D+G_load3D+Gketidai1+Gketidai2+Gketidai3]; %% 需求侧响应补偿成本 公式(1) F_IDSR = sdpvar(1,1); C=[C,F_IDSR==0.432*sum(abs(Ekexuejian1)+abs(Ekexuejian2)+abs(Ekexuejian3))+0.060*sum(abs(Ekepingyi1)+abs(Ekepingyi2)+abs(Ekepingyi3))/2+0.120*sum(abs(Eketidai1)+abs(Eketidai2)+abs(Eketidai3)) ]; % %% 电热冷气能源互传成本 E1yuan = sdpvar(1,24);
3 程序结果
4 下载链接
见下方二维码!