考虑区域多能源系统集群协同优化的联合需求侧响应模型(matlab代码)-阿里云开发者社区

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);