考虑阶梯式碳交易机制与电制氢的综合能源系统热电优化MATLAB复现程序

%% [复现]考虑阶梯式碳交易机制与电制氢的综合能源系统热电优化 %电力自动化设备,41(09) clc clear close all %% 决策变量初始化 P_CHP_e=sdpvar(1,24); %CHP的输出电功率 P_CHP_h=sdpvar(1,24); %CHP的输出热功率 P_g_CHP=sdpvar(1,24); %CHP消耗天然气功率 P_e_EL=sdpvar(1,24); %EL设备的耗电量 P_EL_H=sdpvar(1,24); %EL电解槽的产氢功率 P_H_MR=sdpvar(1,24); %输入MR设备的氢能功率 P_MR_g=sdpvar(1,24); %MR设备输出的天然气功率 P_H_HFC=sdpvar(1,24); %输入HFC设备的氢能功率 P_HFC_e=sdpvar(1,24); %HFC设备输出的电功率 P_HFC_h=sdpvar(1,24); %HFC设备输出的热功率 P_DG=sdpvar(1,24); %风电消纳功率 P_g_GB=sdpvar(1,24); %输入GB设备的天然气功率 P_GB_h=sdpvar(1,24); %GB设备输出的热功率 %储能部分(电ES1、热ES2、气ES3、氢ES4) P_ES1_cha=sdpvar(1,24);P_ES2_cha=sdpvar(1,24);P_ES3_cha=sdpvar(1,24);P_ES4_cha=sdpvar(1,24); %充放功率 P_ES1_dis=sdpvar(1,24);P_ES2_dis=sdpvar(1,24);P_ES3_dis=sdpvar(1,24);P_ES4_dis=sdpvar(1,24); S_1=sdpvar(1,24);S_2=sdpvar(1,24);S_3=sdpvar(1,24);S_4=sdpvar(1,24); %各储能的实时容量状态 %引入充放标志二进制变量 B_ES1_cha=binvar(1,24);B_ES2_cha=binvar(1,24);B_ES3_cha=binvar(1,24);B_ES4_cha=binvar(1,24); %充标志 B_ES1_dis=binvar(1,24);B_ES2_dis=binvar(1,24);B_ES3_dis=binvar(1,24);B_ES4_dis=binvar(1,24); %放标志 P_e_buy=sdpvar(1,24); %购电功率 P_g_buy=sdpvar(1,24); %购气功率 %% 导入常数参数 %电负荷 P_e_load=[717.451523545706;695.290858725762;689.750692520776;698.060941828255;745.152354570637;808.864265927978;836.565096952909;872.576177285319;886.426592797784;900.277008310249;894.736842105263;883.656509695291;875.346260387812;864.265927977839;864.265927977839;868.421052631579;876.731301939058;889.196675900277;880.886426592798;864.265927977839;836.565096952909;817.174515235457;772.853185595568;745.152354570637]'; %热负荷 P_h_load=[864.265927977839;941.828254847645;958.448753462604;955.678670360111;988.919667590028;997.229916897507;903.047091412742;833.795013850416;786.703601108033;703.601108033241;664.819944598338;626.038781163435;595.567867036011;590.027700831025;565.096952908587;639.889196675900;714.681440443213;806.094182825485;811.634349030471;831.024930747922;811.634349030471;808.864265927978;800.554016620499;808.864265927978]'; %天然气负荷 P_g_load=[229.916897506925;224.376731301939;216.066481994460;221.606648199446;224.376731301939;252.077562326870;268.698060941828;288.088642659280;299.168975069252;288.088642659280;293.628808864266;282.548476454294;279.778393351801;271.468144044321;271.468144044321;268.698060941828;277.008310249307;293.628808864266;307.479224376731;304.709141274238;293.628808864266;285.318559556787;277.008310249307;265.927977839335]'; %风电预测出力 P_DG_max=[850.415512465374;864.265927977839;886.426592797784;891.966759002770;894.736842105263;849.030470914127;833.795013850416;653.739612188366;556.786703601108;501.385041551247;432.132963988920;310.249307479224;240.997229916897;252.077562326870;265.927977839335;296.398891966759;343.490304709141;354.570637119114;426.592797783934;526.315789473684;675.900277008310;742.382271468144;854.570637119114;878.116343490305]'; %购电价格 c_e_buy=[0.38*ones(1,7),0.68*ones(1,4),1.2*ones(1,3),0.68*ones(1,4),1.2*ones(1,4),0.38*ones(1,2)]; %购气价格 c_g_buy=0.35*ones(1,24); %% 导入约束条件 C=[]; C=[C, P_CHP_e==0.92*P_g_CHP, %CHP的电-气能量转换约束 0<=P_g_CHP<=600, %CHP消耗的气功率上下限约束 0.5*P_CHP_e<=P_CHP_h, %热电比上下限 P_CHP_h<=2.1*P_CHP_e, %热电比上下限 -0.2*600<=P_g_CHP(2:24)-P_g_CHP(1:23)<=0.2*600, %CHP的爬坡约束(1-24时段) ]; C=[C, P_EL_H==0.87*P_e_EL, %EL(电解槽)的氢-电能量转换约束 0<=P_e_EL<=500, %EL的消耗电功率的上下限约束 -0.2*500<=P_e_EL(2:24)-P_e_EL(1:23)<=0.2*500, %EL的爬坡约束(1-24时段) ]; C=[C, P_MR_g==0.6*P_H_MR, %MR(甲烷反应器)的气-氢能量转换约束 0<=P_H_MR<=250, %MR消耗的氢功率的上下限约束 -0.2*250<=P_H_MR(2:24)-P_H_MR(1:23)<=0.2*250, %MR的爬坡约束(1-24时段) ]; C=[C, P_HFC_e==0.95*P_H_HFC, %HFC(氢燃料电池)的电-氢能量转换约束 0<=P_H_HFC<=250, %HFC消耗的氢功率上下限约束 0.5*P_HFC_e<=P_HFC_h, %HFC的热电比上下限 P_HFC_h<=2.1*P_HFC_e, %HFC的热电比上下限 -0.2*250<=P_H_HFC(2:24)-P_H_HFC(1:23)<=0.2*250, %HFC的爬坡约束(1-24时段) ]; C=[C, 0<=P_DG<=P_DG_max, %风电出力约束 P_GB_h==0.95*P_g_GB, %GB的热-气能量转换约束 0<=P_g_GB<=800, %GB的出力上下限约束 -0.2*800<=P_g_GB(2:24)-P_g_GB(1:23)<=0.2*800, %GB的爬坡约束(1-24时段) ]; C=[C, 0<=P_ES1_cha<=B_ES1_cha*0.5*450, %储电设备的最大充电功率约束 0<=P_ES2_cha<=B_ES2_cha*0.5*500, %储热设备的最大充热功率约束 0<=P_ES3_cha<=B_ES3_cha*0.5*150, %储气设备的最大充气功率约束 0<=P_ES4_cha<=B_ES4_cha*0.5*200, %储氢设备的最大充氢功率约束 0<=P_ES1_dis<=B_ES1_dis*0.5*450, %储电设备的最大放电功率约束 0<=P_ES2_dis<=B_ES2_dis*0.5*500, %储热设备的最大放热功率约束 0<=P_ES3_dis<=B_ES3_dis*0.5*150, %储气设备的最大放气功率约束 0<=P_ES4_dis<=B_ES4_dis*0.5*200, %储氢设备的最大放氢功率约束 S_1(1)==0.3*450, %储电设备的初始容量 S_2(1)==0.3*500, %储热设备的初始容量 S_3(1)==0.3*150, %储气设备的初始容量 S_4(1)==0.3*200, %储氢设备的初始容量 %始末状态守恒约束 S_1(24)==S_1(1), S_2(24)==S_2(1), S_3(24)==S_3(1), S_4(24)==S_4(1), %充放状态唯一 B_ES1_cha+B_ES1_dis<=1, B_ES2_cha+B_ES2_dis<=1, B_ES3_cha+B_ES3_dis<=1, B_ES4_cha+B_ES4_dis<=1, %储能容量上下限约束 0.2*450<=S_1<=0.9*450, 0.2*500<=S_2<=0.9*500, 0.2*150<=S_3<=0.9*150, 0.2*200<=S_4<=0.9*200, %储能容量变化约束 S_1(2:24)==S_1(1:23)+0.95*P_ES1_cha(2:24)-P_ES1_dis(2:24)/0.95, S_2(2:24)==S_2(1:23)+0.95*P_ES2_cha(2:24)-P_ES2_dis(2:24)/0.95, S_3(2:24)==S_3(1:23)+0.95*P_ES3_cha(2:24)-P_ES3_dis(2:24)/0.95, S_4(2:24)==S_4(1:23)+0.95*P_ES4_cha(2:24)-P_ES4_dis(2:24)/0.95, ]; C=[C, P_e_buy==P_e_load+P_e_EL+P_ES1_cha-P_ES1_dis-P_DG-P_CHP_e-P_HFC_e, %电功率平衡约束 P_HFC_h+P_CHP_h+P_GB_h==P_h_load+P_ES2_cha-P_ES2_dis, %热功率平衡约束 P_g_buy==P_g_load+P_ES3_cha-P_ES3_dis+P_g_CHP+P_g_GB-P_MR_g, %气功率平衡约束 P_EL_H==P_H_MR+P_H_HFC+P_ES4_cha-P_ES4_dis, %氢功率平衡约束 0<=P_e_buy<=5000, %购电功率约束 0<=P_g_buy<=5000, %购气功率约束 ]; %% 导入目标函数 %公式5 碳排放权配额模型 E_e_buy=0.728*sum(P_e_buy); %购电配额 E_CHP=0.102*3.6*sum(P_CHP_h+6/3.6*P_CHP_e); %CHP配额 E_GB=0.102*3.6*sum(P_GB_h); %GB配额 E_IES=E_e_buy+E_CHP+E_GB; %IES总碳排放配额 %公式6 实际碳排放模型 E_e_buy_a=1.08*sum(P_e_buy); E_CHP_a=0.065*3.6*sum(P_CHP_h+6/3.6*P_CHP_e); E_GB_a=0.065*3.6*sum(P_GB_h); E_MR_a=1*sum(P_MR_g); %实际MR减少的碳排放 E_IES_a=E_e_buy_a+E_CHP_a+E_GB_a; E=E_IES_a-E_MR_a; %实际IES总碳排放 %阶梯碳交易成本(分段线性化) E_v=sdpvar(1,5)%每段区间内的长度,分为5段,每段长度是2000 lamda=0.250;%碳交易基价 C=[C, E==sum(E_v),%总长度等于E 0<=E_v(1:4)<=2000,%除了最后一段，每段区间长度小于等于2000 0<=E_v(5), ]; %碳交易成本 C_CO2=0; for v=1:5 C_CO2=C_CO2+(lamda+(v-1)*0.25*lamda)*E_v(v); end %总目标函数 obj=c_e_buy*P_e_buy'+c_g_buy*P_g_buy'+C_CO2+0.2*sum(P_DG_max-P_DG); %购能成本+碳交易成本+弃风成本 %% 求解器相关配置 ops=sdpsettings('solver','cplex','verbose',2,'usex0',0); ops.cplex.mip.tolerances.mipgap=1e-6; %% 进行求解计算 result=optimize(C,obj,ops); if result.problem==0 % problem =0 代表求解成功 % do nothing! else error('求解出错'); end %% 画图输出运行结果 k_CHP=double(P_CHP_h./P_CHP_e); k_HFC=double(P_HFC_h./P_HFC_e); figure plot(k_CHP,'-pg','LineWidth',1.5) xlabel('时段') ylabel('CHP热电比') %电功率情况 Plot_e=[]; for t=1:24 Plot_e(1,t)=P_e_buy(t); Plot_e(2,t)=-P_e_EL(t); Plot_e(3,t)=-P_ES1_cha(t)+P_ES1_dis(t); Plot_e(4,t)=P_DG(t); Plot_e(5,t)=P_CHP_e(t); Plot_e(6,t)=P_HFC_e(t); end Plot_e=Plot_e'; figure bar(Plot_e,'stacked'); hold on plot(P_e_load,'r-o','LineWidth',1.5); xlabel('时间/h'); ylabel('功率/kW'); title('最优电负荷情况'); legend('购电','电解槽','储电','风电','CHP','氢燃料电池','电负荷'); %热功率情况 Plot_h=[]; for t=1:24 Plot_h(1,t)=P_H