整理以下代码:%% 含电转气-碳捕集(P2G-CCS)耦合的综合能源系统低碳经济调度
clc
clear
%% 参数设定
P_C2P=sdpvar(1,24); %CCPP-P2G系统的总能耗
P_P2G=sdpvar(1,24); %P2G设备的能耗
P_CC=sdpvar(1,24); %碳捕集能耗
P_OP=sdpvar(1,24); %碳捕集运行能耗
P_GN=sdpvar(1,24); %碳捕集电厂净出力
P_G=sdpvar(1,24); %碳捕集电厂等效出力
P_GC=sdpvar(1,24); %碳捕集电厂提供的碳捕集能耗
P_Galpha=sdpvar(1,24); %碳捕集电厂的烟气处理能耗
Q_CC=sdpvar(1,24); %CCPP-P2G系统捕集的CO2的总量,单位为吨
Q_P2Gsum=sdpvar(1,24); %P2G消耗的CO2的总量,单位为吨
V_P2G=sdpvar(1,24); %P2G生成的天然气
P_alpha=sdpvar(1,24); %总烟气处理能耗
alpha_1=sdpvar(1,24); %t时刻进行烟气处理的烟气量中由 碳捕集电厂运行产生的烟气提供的部分
alpha_2=sdpvar(1,24); %流入烟气存储罐的气量
alpha_3=sdpvar(1,24); %t时刻进行烟气处理的烟气量中由储气装置提供的烟气量
P_WC=sdpvar(1,24); %风电机组提供的碳捕集能耗
P_VC=sdpvar(1,24); %光伏机组提供的碳捕集能耗
P_WIC=sdpvar(1,24); % 碳捕集电厂中火电部分提供的碳捕集能耗
P_Valpha=sdpvar(1,24); %光伏机组提供的烟气处理能耗
P_Walpha=sdpvar(1,24); %风力机组提供的烟气处理能耗
P_WIalpha=sdpvar(1,24); % 碳捕集电厂中火电部分提供的烟气处理能耗
P_WN=sdpvar(1,24); %风力机组的上网功率
P_VN=sdpvar(1,24); %光伏机组的上网功率
P_WIN=sdpvar(1,24); % 碳捕集电厂的上网功率
P_WI=sdpvar(1,24); % 碳捕集发电出力
Q_N=sdpvar(1,24); %碳捕集电厂的CO2的净排放量,单位为t/h
P_PH=sdpvar(1,24); %CHP机组的输出总功率
P_CHP=sdpvar(1,24); %CHP机组的输出电功率
H_CHP=sdpvar(1,24); %CHP机组的输出热功率
V_CHP=sdpvar(1,24); %CHP机组消耗的天然气量
H_GB=sdpvar(1,24); %燃气锅炉的输出热功率
V_GB=sdpvar(1,24); %燃气锅炉消耗的天然气量
S_ES=sdpvar(1,24); %电储能在t时段末的蓄电量,单位为MW
P_ESC=sdpvar(1,24); %电储能的充电功率
P_ESD=sdpvar(1,24); %电储能的放电功率
S_TS=sdpvar(1,24); %热储能在t时段末的蓄热量,单位为MW
H_TSC=sdpvar(1,24); %热储能的充热功率
H_TSD=sdpvar(1,24); %热储能的放热功率
P_EM=sdpvar(1,24); %系统在电网的购电量
P_CUT=sdpvar(1,24); %各级中断负荷功率之和
P_Cmax=sdpvar(1,24); %碳捕集系统的运行能耗上限
miu_ESC=binvar(1,24); %充电的布尔变量
miu_ESD=binvar(1,24); %放电的布尔变量
miu_TSC=binvar(1,24); %充热的布尔变量
miu_TSD=binvar(1,24); %放热的布尔变量
V_WIalpha=sdpvar(1,24); %烟气存储罐储气量
Q_CS=sdpvar(1,24); %流入烟气存储罐的气量
Q_P2G=sdpvar(1,24); %流入烟气存储罐的气量
P_A=15*ones(1,24); %CCPP-P2G系统固定能耗(因占比较少设为定值),单位为MW
%风电机组的预测出力
P_W=[232.75,247.44,219.09,188.78,239.58,232.84,188.52,159.84,111.45,51.23,119.88,137.29,141.39,115.78,135.24,143.44,151.64,195.69,159.70,180.94,203.38,193.64,155.32,247.43];
%光伏机组的预测出力
P_V=[0,0,0,0,0,22,63,97,110,118,128,132,133,136,131,133,120,85,37,0,0,0,0,0];
%电负荷
P_EL=[457,319,296,228,184,297,406,509,607,687,803,857,845,793,832,801,795,731,640,593,554,518,525,409];
%热功率
H_HL=[109,131,158,153,139,121,111,98,82,57,22,12,42,62,89,99,122,131,148,160,139,131,119,74];
%购电价
k_EM=[38.85,39.18,36.89,35.57,39.84,43.77,51.31,64.10,74.59,77.21,85.41,89.02,82.46,80.49,83.11,81.80,78.52,73.93,69.67,76.89,74.26,66.39,55.57,46.72];
S_ES_init=60;
S_TS_init=30;
P_WA=0.2*P_W; %t时段弃风功率
P_VA=0.2*P_V; %t时段弃光功率
%% 约束条件
C=[];
%CCPP-P2G系统能耗以及CCPP出力
for t=1:24
C=[C,
P_C2P(t)==P_P2G(t)+P_CC(t), %CCPP-P2G系统总能耗约束
P_P2G(t)==P_WA(t)+P_VA(t), %P2G消纳的弃风、弃光量约束
P_CC(t)==P_A(t)+P_OP(t), %碳捕集能耗约束
P_GN(t)==P_G(t)-P_GC(t)-P_Galpha(t), %碳捕集电厂功率约束
];
end
%CCPP-P2G系统碳利用量以及天然气生成量
for t=1:24
C=[C,
Q_CC(t)==P_OP(t)/0.269, %CCPP-P2G系统捕集的CO2的总量与耗能约束 P_OP是碳捕集运行能耗
Q_P2Gsum(t)==0.2*0.6*P_P2G(t), %P2G设备的消耗的CO2量和电功率约束
V_P2G(t)==3.6*0.6*P_P2G(t)/39, %P2G设备的生成天然气
];
end
% 碳捕集电厂烟气处理模型
for t=1:24
C=[C,P_alpha(t)==0.513*(alpha_1(t)+alpha_3(t)),]; %烟气处理系统能耗
end
%碳捕集-风电-光伏联合运行策略
for t=1:24
C=[C,
P_GC(t)+P_WC(t)+P_VC(t)+P_WIC(t)==P_CC(t), %碳捕集能耗等式约束
P_OP(t)==0.269*Q_CC(t), %CCPP-P2G系统捕集的CO2的总量与耗能约束
P_Valpha(t)+P_Walpha(t)+P_Galpha(t)+P_WIalpha(t)==P_alpha(t), %烟气处理能耗等式约束
P_WN(t)+P_WC(t)+P_Walpha(t)==P_W(t), %风力机组的出力约束
P_VN(t)+P_VC(t)+P_Valpha(t)==P_V(t), %光伏机组的出力约束
P_WIN(t)+P_WIC(t)+P_WIalpha(t)==P_WI(t), % 碳捕集电厂的出力约束
Q_N(t)==0.96*P_G(t)-Q_CC(t), %碳捕集电厂的碳排放约束
];
end
%CHP机组和燃气锅炉模型
for t=1:24
C=[C,
P_PH(t)==P_CHP(t)+H_CHP(t), %CHP机组的输出功率约束
P_CHP(t)==V_CHP(t)*39*0.35, %CHP机组的输出电功率约束
H_CHP(t)==V_CHP(t)*39*0.40, %CHP机组的输出热功率约束
H_GB(t)==V_GB(t)*39*0.40, %CHP机组的输出热功率约束
];
end
%储能装置模型
for t=2:24
C=[C,
S_ES(t)==S_ES(t-1)*(1-0.001)+0.95*P_ESC(t)-P_ESD(t)/0.95, %电储能运行约束
S_TS(t)==S_TS(t-1)*(1-0.01)+0.88*H_TSC(t)-H_TSD(t)/0.88, %热储能运行约束
];
end
%电功率和热功率平衡约束
for t=1:24
C=[C,
P_GN(t)+P_WIN(t)+P_CHP(t)+P_WN(t)+P_VN(t)+P_ESD(t)+P_EM(t)==P_EL(t)+P_ESC(t), %电功率平衡约束
H_CHP(t)+H_GB(t)+H_TSD(t)==H_HL(t)+H_TSC(t), %热功率平衡约束
];
end
%碳捕集电厂约束
for t=1:24
C=[C,
100<=P_G(t)<=400, %碳捕集电厂出力上下限约束
0<=Q_CC(t)<=0.96*400,
15<=P_GC(t)+P_WC(t)+P_VC(t)+P_WIC(t)<=P_Cmax(t), %碳捕集系统的运行能耗上下限
P_Cmax(t)==0.269*0.96*P_G(t), %碳捕集系统的运行能耗上限赋值
];
end
for t=2:24
C=[C,
-60<=P_G(t)-P_G(t-1)<=60, %碳捕集电厂出力爬坡速率约束
-65<=P_GC(t)+P_WC(t)+P_VC(t)+P_WIC(t)-P_GC(t-1)-P_WC(t-1)-P_VC(t-1)-P_WIC(t-1)<=65, %碳捕集电厂碳捕集能耗爬坡速率约束
];
end
%CHP机组电热出力以及爬坡约束
for t=1:24
C=[C,
0<=P_CHP(t)<=140, %CHP的电功率出力约束
0<=H_CHP(t)<=160, %CHP的热功率出力约束
];
end
for t=2:24
C=[C,-120<=P_PH(t)-P_PH(t-1)<=120,]; %CHP机组出力爬坡速率约束
end
%燃气锅炉热出力以及爬坡约束
for t=1:24
C=[C,0<=H_GB(t)<=200,]; %燃气锅炉的热功率出力约束
end
for t=2:24
C=[C,-100<=H_GB(t)-H_GB(t-1)<=100,]; %燃气锅炉的爬坡速率约束
end
%P2G运行约束
C=[C,0<=P_P2G<=200,]; %P2G运行功率上下限
%电储能和热储能约束
for t=1:24
C=[C,
0<=P_ESC(t)<=miu_ESC(t)*80, %充电上下限约束
0<=P_ESD(t)<=miu_ESD(t)*80, %放电上下限约束
0<=miu_ESC(t)+miu_ESD(t)<=1,%不可能同时出现充放电
10<=S_ES(t)<=120, %储电量上下限约束
0<=H_TSC(t)<=miu_TSC(t)*40, %充热上下限约束
0<=H_TSD(t)<=miu_TSD(t)*40, %放热上下限约束
0<=miu_TSC(t)+miu_TSD(t)<=1,%不可能同时出现充放热
10<=S_TS(t)<=50, %储热量上下限约束
];
end
C=[C,S_ES_init==S_ES(24),S_ES(1)==S_ES_init*(1-0.001)+0.95*P_ESC(1)-P_ESD(1)/0.95,]; %电储能和热储能初始和结束的约束
C=[C,S_TS_init==S_TS(24),S_TS(1)==S_TS_init*(1-0.01)+0.88*H_TSC(1)-H_TSD(1)/0.88,];
% 碳捕集电厂约束
C=[C,sum(P_WI)<=1500,]; % 碳捕集电厂日总出力上限
for t=1:24
C=[C,
60<=P_WI(t)<=100, % 碳捕集电站的出力上下限
0.1*400<=V_WIalpha(t)<=0.9*400, %储气装置容量的上下限
0<=alpha_1(t)<=160,
0<=alpha_2(t)<=160, %流入储气装置的上下限
0<=alpha_3(t)<=160, %流出储气装置的上下限
];
end
for t=2:24
C=[C,-40<=P_WI(t)-P_WI(t-1)<=40,]; % 碳捕集电站的爬坡速率约束(文中未给参数,暂取40)
end
for t=2:24
C=[C,V_WIalpha(t)==V_WIalpha(t-1)+alpha_2(t)-alpha_3(t),]; %储气装置的容量变化约束
end
C=[C,V_WIalpha(1)==0.4*400+alpha_2(1)-alpha_3(1),];
C=[C,Q_CC==Q_CS+Q_P2G,Q_CS>=0,Q_P2G>=0,0<=P_VC,0<=P_WIC,0<=P_GC,];
C=[C,P_EM>=0,Q_P2Gsum-Q_P2G>=0,P_WIN>=0,0<=P_Valpha,0<=P_Walpha,0<=P_Galpha,0<=P_WIalpha,0<=P_WC,];
%% 目标函数
C_F=200*24+17*sum(P_G)+0.04*P_G*(P_G');%碳捕集电厂燃料成本C_F
C_WI=19.8*(0.96-0.76)*sum(P_WI); % 碳捕集电厂成本
V_BUY=V_CHP+V_GB-V_P2G; %天然气购买量
C_H=0.419*sum(V_BUY); %CHP机组和燃气锅炉成本
Q_BUY=Q_P2Gsum-Q_P2G; %购买的二氧化碳的量
C_PG=120*sum(Q_BUY)+20*sum(P_P2G); %P2G成本
C_CS=4.89*sum(Q_CS); %碳封存成本
C_W=21.4*sum(P_W)+14.2*sum(P_V); %系统运行维护成本
C_M=sum(k_EM.*P_EM); %电力市场购电成本
%碳排放权配额模型
E_e_buy=0.68*C_M; %购电配额
E_CHPGB=0.102*3.6*sum(V_BUY); %CHP和GB配额
Q_Q=0.76*P_GN; %碳捕集电厂的碳排放配额
z=sum(V_BUY);
E_IES=E_e_buy+E_CHPGB+sum(Q_Q); %IES总碳排放配额
%实际碳排放模型
E_e_buy_a=0.711*C_M;
E_CHPGB_a=0.249*3.6*sum(V_BUY);
E_IES_a=E_e_buy_a+E_CHPGB_a+sum(Q_N);
E=E_IES-E_IES_a %实际IES总碳排放
%阶梯碳交易成本(分段线性化)
E_v=sdpvar(1,5)%每段区间内的长度,分为5段,每段长度是2000
lamda=19.8;%碳交易基价
C=[C,
E==sum(E_v),%总长度等于E
0<=E_v(1:4)<=2000,%除了最后一段,每段区间长度小于等于1000
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_F-C_CO2+C_WI+C_H+C_PG+C_CS+C_W+C_M;
%% 求解
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 代表求解成功
else
error('求解出错!');
end
%% 数据分许与画图部分
%碳捕集能耗
P_CC=-P_CC;
P_WC=double(P_WC);
P_VC=double(P_VC);
P_WIC=double(P_WIC);
P_CC=double(P_CC);
for t=1:24
Plot_p(1,t)=P_WC(t);
Plot_p(2,t)=P_VC(t);
Plot_p(3,t)=P_WIC(t);
end
Plot_p=Plot_p';
figure
b=bar(Plot_p,'stacked');
b(1).FaceColor = [.9 .8 .1];
b(2).FaceColor = [.0 .9 .0];
b(3).FaceColor = [.9 .1 .3];
hold on
bar(P_CC,'stacked');
xlabel('时段');
ylabel('功率/MW');
legend('风电提供能耗','光伏提供能耗',' 碳捕集电厂提供能耗','碳捕集总能耗');
legend('boxoff');
title('综合能源系统的碳捕集能耗');
box off
%CO2排放量和处理量
Q_N=double(Q_N);
Q_CC=double(Q_CC);
figure
bar([Q_N',Q_CC'])
xlabel('时段');
ylabel('二氧化碳量/t');
legend('二氧化碳排放量','二氧化碳捕集量');
legend('boxoff');
title('CO2排放量和处理量');
box off
% 碳捕集电厂出力和能量市场购电量
P_WI=double(P_WI);
P_WIN=double(P_WIN);
P_EM=double(P_EM);
figure
d=bar(P_EM,'stacked');
d(1).FaceColor = [.9 .1 .3];
hold on
plot(P_WI,'b-.','LineWidth',1.5);
hold on
plot(P_WIN,'m-.o','LineWidth',1.5);
xlabel('时段');
ylabel('功率/MW');
legend('综合能源系统购电量',' 碳捕集电厂等效出力',' 碳捕集电厂净出力');
legend('boxoff');
title('综合能源系统购电和电厂出力情况');
box off
% 电负荷和各供电单元出力
% P_GN(t)+P_WIN(t)+P_CHP(t)+P_WN(t)+P_VN(t)+P_ESD(t)+P_EM(t)==P_EL(t)+P_ESC(t), %电功率平衡约束
P_GN=double(P_GN);
P_WIN=double(P_WIN);
P_CHP=double(P_CHP);
P_WN=double(P_WN);
P_VN=double(P_VN);
P_ESD=double(P_ESD);
P_EM=double(P_EM);
P_EL=double(P_EL);
P_ESC=double(P_ESC);
for t=1:24
Plot_PLoad(1,t)=P_GN(t);
Plot_PLoad(2,t)=P_WIN(t);
Plot_PLoad(3,t)=P_CHP(t);
Plot_PLoad(4,t)=P_WN(t);
Plot_PLoad(5,t)=P_VN(t);
Plot_PLoad(6,t)=P_ESD(t);
Plot_PLoad(7,t)=P_EM(t);
Plot_PLoad(8,t)=-1*P_ESC(t);
Plot_PLoad(9,t)=-1*P_EL(t);
end
Plot_PLoad=Plot_PLoad';
figure
f1=bar(Plot_PLoad,'stacked');
f1(1).FaceColor = [.9 .1 .3];
f1(2).FaceColor = [.2 .7 .2];
f1(3).FaceColor = [.2 .2 .7];
f1(4).FaceColor = [.9 .8 .1];
f1(5).FaceColor = [.6 .8 .1];
hold on
plot(P_EL,'m-.s','LineWidth',2);
xlabel('时间/h');
ylabel('功率/MW');
title('电负荷和各机组供电出力');
hold on
legend('碳捕集电厂净出力','碳捕集电厂的出力','CHP机组电出力','风力机组出力',...
'光伏机组出力','电储能的放电功率','电网的购电量','电储能的充电功率','电负荷功率');
legend('boxoff');
grid on
box off
% 热负荷和各供热单元出力
% H_CHP(t)+H_GB(t)+H_TSD(t)==H_HL(t)+H_TSC(t), %热功率平衡约束
H_CHP=double(H_CHP);
H_GB=double(H_GB);
H_TSD=double(H_TSD);
H_TSC=double(H_TSC);
P_CHP=double(P_CHP);
for t=1:24
Plot_HeatLoad(1,t)=H_CHP(t);
Plot_HeatLoad(2,t)=H_GB(t);
Plot_HeatLoad(3,t)=H_TSD(t);
Plot_HeatLoad(4,t)=-1*H_TSC(t);
Plot_HeatLoad(5,t)=-1*H_HL(t);
end
Plot_HeatLoad=Plot_HeatLoad';
figure
f=bar(Plot_HeatLoad,'stacked');
f(1).FaceColor = [.9 .1 .3];
f(2).FaceColor = [.2 .7 .2];
f(3).FaceColor = [.2 .2 .7];
f(4).FaceColor = [.9 .8 .1];
f(5).FaceColor = [.6 .8 .1];
hold on
plot(H_HL,'m-.s','LineWidth',2);
xlabel('时间/h');
ylabel('功率/MW');
title('热负荷和各机组供热出力');
hold on
legend('CHP机组热出力','燃气锅炉热出力','热储能放热','热储能储热','热负荷功率');
legend('boxoff');
grid on
box off
% 电储能
P_ESD=-P_ESD;
P_ESC1=double(P_ESC);P_ESD1=double(P_ESD);S_ES=double(S_ES);
for t=1:24
P(1,t)=P_ESC1(t);
P(2,t)=S_ES(t);
end
P=P';
figure
a=bar(P,'stacked');
a(1).FaceColor = [.2 .7 .2];
a(2).FaceColor = [.9 .1 .3];
hold on
bar(P_ESD1,'stacked');
hold on
plot(S_ES,'c-.x','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/MW');
title('IES电储能装置功率状态');
hold on
legend('充电功率','电储能蓄电量','放电功率');
legend('boxoff');
grid on
box off
% 热储能
H_TSD=-H_TSD;
S_TS=double(S_TS);H_TSC=double(H_TSC);H_TSD=double(H_TSD);
for t=1:24
R(1,t)=H_TSC(t);
R(2,t)=S_TS(t);
end
R=R';
figure
bar(R,'stacked');
hold on
a=bar(H_TSD,'stacked');
a(1).FaceColor = [.9 .1 .3];
hold on
plot(S_TS,'m-.o','LineWidth',1.5);
xlabel('时间/h');
ylabel('功率/MW');
title('IES热储能装置功率状态');
hold on
legend('充热功率','放热功率','热储能储热量');
legend('boxoff');
grid on
box off