%% CCS+CANQING
clc
clear all
P_e_load=sdpvar(1,24); %微网经过需求响应后实际的电负荷
P_h_load=sdpvar(1,24); %微网经过需求响应后实际的热负荷
P_e_cut=sdpvar(1,24); %微网的可削减电负荷
P_e_tran=sdpvar(1,24); %微网的可转移电负荷
P_h_cut=sdpvar(1,24); %微网的可削减热负荷
P_h_tran=sdpvar(1,24); %微网的可转移电负荷
S_bat=sdpvar(1,24); %微网中的储电设备的储电余量
P_bat_cha=sdpvar(1,24); %储电设备的充电功率
P_bat_dis=sdpvar(1,24); %储电设备的放电功率
B_bat_cha=binvar(1,24); %储电设备的放电状态位,取1时为放电,0为未放电
B_bat_dis=binvar(1,24); %储电设备的充电状态位,取1时为充电,0为未充电
S_tes=sdpvar(1,24); %微网中的储电设备的储电余量
P_tes_cha=sdpvar(1,24); %储电设备的充电功率
P_tes_dis=sdpvar(1,24); %储电设备的放电功率
B_tes_cha=binvar(1,24); %储电设备的放电状态位,取1时为放电,0为未放电
B_tes_dis=binvar(1,24); %储电设备的充电状态位,取1时为充电,0为未充电
P_e_buy=sdpvar(1,24); %微网向外电网的购买的电功率
P_e_sell=sdpvar(1,24); %微网向外电网的售出的电功率
B_e_buy=binvar(1,24); %购电标志
B_e_sell=binvar(1,24); %售电标志
P_h_buy=sdpvar(1,24); %微网向外电网的购买的热功率
P_h_sell=sdpvar(1,24); %微网向外电网的售出的热功率
B_h_buy=binvar(1,24); %购热标志
B_h_sell=binvar(1,24); %售re标志
P_e_PV=sdpvar(1,24); %风力的实际出力值
V_GT_CH4=sdpvar(1,24);
V_GT_H2=sdpvar(1,24);
P_e_GT=sdpvar(1,24); %燃气轮机的发电功率
P_e_GT_SUP=sdpvar(1,24); %燃气轮机实际的发电功率
P_h_GT=sdpvar(1,24); %燃气轮机的产热功率
V_GB_CH4=sdpvar(1,24);
V_GB_H2=sdpvar(1,24);
P_h_GB=sdpvar(1,24); %锅炉的产热功率
%两阶段P2G
% WT
P_e_WT=sdpvar(1,24);
P_e_PV=sdpvar(1,24);
% EL
P_g_EL=sdpvar(1,24); %pem制氢功率
P_e_EL=sdpvar(1,24); %pem耗电功率
P_e_ELbuy=sdpvar(1,24);
% P_h_EL=sdpvar(1,24); %
V_EL_H2=sdpvar(1,24);
% CCS
P_e_CCS=sdpvar(1,24); %CCS耗电功率
M_CCS_CO2=sdpvar(1,24); %捕集CO2功率
% 购气
V_BUY_CH4=sdpvar(1,24),
%% 计算参数
M=1000000;
L_CH4=35.86*1000; %KJ/M3
L_H2=120*1000; %KJ/M3
V_CO2=0.5; %在 100℃ 和 1 atm 条件下,1 kg CO₂ 的体积约为 0.696 m³。
V_CH4=1.4; %标况 1kg天然气体积1.4M3
V_H2=11.17; %1 kg 氢气的体积约为 11.17 m³
K_GB_H=0.85;
P_h_GBmin=200;
P_h_GBmax=1000;
P_h_GBdown=-200;
P_h_GBup=200;
K_GT_E=0.45;
P_e_GTmin=1800;
P_e_GTmax=4500;
P_e_GTdown=-900;
P_e_GTup=900;
K_EL_H2=0.75;
P_e_ELbuymax=1000;
P_e_ELmin=0;
P_e_ELmax=2000;
P_e_ELdown=-800;
P_e_ELup=800;
A_CCS_CO2=1.02; %kWh/kg 是捕获 CO2 所消耗电能的转换系数。 kWh/kg
P_e_CCSmin=0;
P_e_CCSmax=600;
P_e_CCSdown=-240;
P_e_CCSup=240;
S_bat_int=400;
S_batmin=200;
S_batmax=1200;
P_bat_chamax=600;
P_bat_dismax=600;
S_tes_int=400;
S_tesmin=200;
S_tesmax=800;
P_tes_chamax=400;
P_tes_dismax=400;
P_e_buymax=1000;
P_e_sellmax=1000;
P_h_buymax=800;
P_h_sellmax=800;
P_e_load0=3000*[0.4476 0.4628 0.4389 0.4778 0.5440 0.7732 0.8071 0.8324 0.8476 0.8645 0.9238 0.9441 1.0000 0.9710 0.8417 0.7736 0.8346 0.8294 0.7630 0.6677 0.5724 0.5315 0.4905 0.4785];
P_h_load0=2500*[0.5920 0.6055 0.5843 0.5959 0.6578 0.7467 0.8202 0.9555 0.9594 1.0000 0.9613 0.9439 0.9865 0.9555 0.9942 0.9265 0.9884 0.9826 0.9227 0.8202 0.7312 0.6868 0.6500 0.6365 ];
P_e_WT_PRE=1500*[0.531959039 0.717892503 0.565004111 0.563108967 0.630595044 0.500832957 0.572302704 0.543549865 0.418886088 0.61296521 0.784283688 0.607356299 0.688252427 0.728935073 0.536267485 0.672290798 0.562845084 0.627322335 0.602754432 0.410946372 0.358698603 0.436885332 0.556948548 0.466038402];
P_e_PV_PRE=800*[0 0 0 0 0 0.096994745 0.268134132 0.558446021 0.555730945 0.405237615 0.599796904 0.694387523 0.639471338 0.604943309 0.531663197 0.547697019 0.446525826 0.259199023 0.125313905 0 0 0 0 0];
C_e=[0.4*ones(1,5),1*ones(1,5),1.4*ones(1,4),1*ones(1,6),0.4*ones(1,4)];
C_e1=[0.25*ones(1,24)];
C_h=[0.27*ones(1,9),0.35*ones(1,8),0.15*ones(1,7)];
C_h1=[0.1*ones(1,24)];
C_g=[1.4*ones(1,9),1.8*ones(1,3),2.2*ones(1,4),1.8*ones(1,4),2.2*ones(1,3),1.8*ones(1,1)];
C_OM_WT=0.03;
C_OM_PV=0.01;
C_OM_GT=0.044;
C_OM_GB=0.035;
C_OM_ESS=0.011;
C_OM_TES=0.018;
C_OM_EL=0.022;
C_OM_CCS=0.031;
%% 约束
%
C=[];
C=[C,
P_e_GT==P_e_GT_SUP+P_e_CCS,
P_e_GT>=0,
];
% GT
for i=1:24
C=[C,
0.1*(V_GT_H2(i)+V_GT_CH4(i))<=V_GT_H2(i)<=0.2*(V_GT_H2(i)+V_GT_CH4(i))
];
end
C=[C,
P_e_GT==(V_GT_H2*L_H2+V_GT_CH4*L_CH4)*K_GT_E,
%0.1<=GT_H2<=0.2,
P_e_GTmin<=P_e_GT<=P_e_GTmax,
P_e_GTdown<=P_e_GT(2:24)-P_e_GT(1:23)<=P_e_GTup,
];
% WT
C=[C,
0<=P_e_WT<=P_e_WT_PRE,
0<=P_e_PV<=P_e_PV_PRE
];
% EL机组
C=[C,
0<=P_e_ELbuy<=P_e_ELbuymax,
P_e_EL*K_EL_H2==V_EL_H2*L_H2,
P_e_ELmin<=P_e_EL<=P_e_ELmax,
P_e_ELdown<=P_e_EL(2:24)-P_e_EL(1:23)<=P_e_ELup,
];
% MR机组(看不太到意义)
% C=[C,
% P_e_MR*K_MR_CH4*3600==V_MR_CH4*L_CH4,
% P_e_MRmin<=P_e_MR<=P_e_MRmax,
% P_e_MRdown<=P_e_MR(2:24)-P_e_MR(1:23)<=P_e_MRup,
% ];
% CCS机组
C=[C,
P_e_CCS==A_CCS_CO2*M_CCS_CO2,
P_e_CCSmin<=P_e_CCS<=P_e_CCSmax,
P_e_CCSdown<=P_e_CCS(2:24)-P_e_CCS(1:23)<=P_e_CCSup,
0<=M_CCS_CO2<=(P_e_GT+0.15*P_h_GT)/A_CCS_CO2,%最大碳捕集约束
];
% % EL,MR与CCS耦合特性
% %
% C=[C,
% V_CCS_CO2==P_e_MR*A_MR_CO2,
% V_MR_H2*L_H2==L_CH4*V_MR_CH4*K_MR_CH4,
%
% %V_CCS_CO2==V_MR_CH4, 不要反复定义可能会导致约束冲突无法求解
% ];
% 热电耦合特性
C=[C,
P_e_GTmin-P_e_CCS<=P_e_GT_SUP<=P_e_GTmax-P_e_CCS,
max((P_e_GTmin-0.15*P_h_GT-P_e_CCS),(0.85*(P_h_GT-450)-P_e_CCS))<=P_e_GT_SUP<=P_e_GTmax-0.2*P_h_GT-P_e_CCS,
max((P_e_GTmin-0.15*P_h_GT)-P_e_CCSmax,(0.85*(P_h_GT-450)-P_e_CCSmax))<=P_e_GT<=P_e_GTmax-0.2*P_h_GT-P_e_CCSmin,
(1/A_CCS_CO2)*max((P_e_GTmin-0.15*P_h_GT-P_e_GT_SUP),(0.85*(P_h_GT-450)-P_e_GT_SUP))<=M_CCS_CO2<=(1/A_CCS_CO2)*(P_e_GTmax-0.20*P_h_GT-P_e_GT_SUP), %产气功率上下限约束 式(17)
];
% GB
for i=1:24
C=[C,
%GB_H2(i)*(V_GB_H2(i)+V_GB_CH4(i))==V_GB_H2(i),
0.02*(V_GB_H2(i)+V_GB_CH4(i))<=V_GB_H2(i)<=0.2*(V_GB_H2(i)+V_GB_CH4(i))
];
end
C=[C,
P_h_GB==(V_GB_H2*L_H2+V_GB_CH4*L_CH4)*K_GB_H,
%GB_H2*(V_GB_H2+V_GB_CH4)==V_GB_H2,
%0.02<=GB_H2<=0.2,
P_h_GBmin<=P_h_GB<=P_h_GBmax,
P_h_GBdown<=P_h_GB(2:24)-P_h_GB(1:23)<=P_h_GBup,
];
% ESS
%上下限约束
C=[C,
S_batmin<=S_bat<=S_batmax,
S_bat(1)==S_bat(24),
];
C=[C,
S_bat(1)==S_bat_int+0.95*P_bat_cha(1)-P_bat_dis(1)/0.95,
];
for t=2:24
C=[C,
S_bat(t)==S_bat(t-1)+0.95*P_bat_cha(t)-P_bat_dis(t)/0.95,
];
end
for t=1:24
C=[C,
0<=P_bat_cha(t)<=P_bat_chamax,
0<=P_bat_cha(t)<=B_bat_cha(t)*M,
0<=P_bat_dis(t)<=P_bat_dismax,
0<=P_bat_dis(t)<=B_bat_dis(t)*M,
B_bat_cha(t)+B_bat_dis(t)<=1,
];
end
% TES
C=[C,
S_tesmin<=S_tes<=S_tesmax,
S_tes(1)==S_tes(24),
];
C=[C,
S_tes(1)==S_tes_int+0.95*P_tes_cha(1)-P_tes_dis(1)/0.95,
];
for t=2:24
C=[C,
S_tes(t)==S_tes(t-1)+0.95*P_tes_cha(t)-P_tes_dis(t)/0.95,
];
end
for t=1:24
C=[C,
0<=P_tes_cha(t)<=P_tes_chamax,
0<=P_tes_cha(t)<=B_tes_cha(t)*M,
0<=P_tes_dis(t)<=P_tes_dismax,
0<=P_tes_dis(t)<=B_tes_dis(t)*M,
B_tes_cha(t)+B_tes_dis(t)<=1,
];
end
% buy sell
for t=1:24
C=[C,
0<=P_e_buy(t)<=P_e_buymax,
0<=P_e_buy(t)<=B_e_buy(t)*M,
0<=P_e_sell(t)<=P_e_sellmax,
0<=P_e_sell(t)<=B_e_sell(t)*M,
B_e_buy(t)+B_e_sell(t)<=1,
0<=P_h_buy(t)<=P_h_buymax,
0<=P_h_buy(t)<=B_h_buy(t)*M,
0<=P_h_sell(t)<=P_h_sellmax,
0<=P_h_sell(t)<=B_h_sell(t)*M,
B_h_buy(t)+B_h_sell(t)<=1,
];
end
% DR
for t=1:24
C=[C,
P_e_load(t)==P_e_load0(t)+P_e_cut(t)+P_e_tran(t), %电负荷功率平衡约束
P_h_load(t)==P_h_load0(t)+P_h_cut(t)+P_h_tran(t), %电负荷功率平衡约束
-0.2*P_e_load0(t)<=P_e_tran(t)<=0.2*P_e_load0(t), %微网的可转移电功率上下限约束
-0.2*P_h_load0(t)<=P_h_tran(t)<=0.2*P_h_load0(t), %微网的可转移电功率上下限约束
0<=P_e_cut(t)<=0.15*P_e_load0(t),
0<=P_h_cut(t)<=0.15*P_h_load0(t),
];
end
C=[C,
sum(P_e_tran)==0,
sum(P_h_tran)==0,
]; %转移的电-热负荷总量为0约束
% 电-热-氢-天然气 平衡
C=[C,
%V_BUY_CH4<=0.9*(V_GT_CH4+V_GB_CH4),
P_e_GT_SUP+P_e_buy+P_bat_dis==P_bat_cha+P_e_sell+P_e_load,
P_h_GT+P_h_GB+P_h_buy+P_tes_dis==P_tes_cha+P_h_sell+P_h_load,
V_EL_H2==V_GT_H2+V_GB_H2,
V_GT_CH4+V_GB_CH4==V_BUY_CH4,
P_e_ELbuy+P_e_WT+P_e_PV==P_e_EL,
];
%% 目标函数
% 购能成本
%碳交易系数
%配额量
% ae_buy=0.823; %购电
ae_gen=0.789;
ah_gen=0.385;
% ah_buy=0.1372;
%实际量
% be_buy=1.08; %购电
be_gen=0.823;
bh_gen=0.1372;
% bh_buy=0.677;
C_e_buy=C_e*P_e_buy'+C_e*P_e_ELbuy'; %购电
C_e_sell=C_e1*P_e_sell';
C_h_buy=C_h*P_h_buy';
C_h_sell=C_h1*P_h_sell'; %购热
C_S=C_e*P_bat_cha'-C_e*P_bat_dis'+C_h*P_tes_cha'-C_h*P_tes_dis'; %储能
C_G=3600*C_g*V_BUY_CH4'; %购气
C_C=-0.25*(sum(P_e_GT)*ae_gen+sum(P_h_GT)*ah_gen+sum(P_h_GB)*ah_gen-sum(3600*V_BUY_CH4/V_CO2)+(sum(P_e_buy)*(ae_gen-be_gen))+(sum(P_h_buy)*(ah_gen-bh_gen)));
C_Ctrade=0.3*sum(M_CCS_CO2);
C_DR=0.15*sum(abs(P_e_tran))+0.1*sum(P_e_cut)+0.08*sum(abs(P_h_tran))+0.04*sum(P_h_cut);
C_OM=C_OM_CCS*sum(P_e_CCS)+C_OM_EL*sum(P_e_EL)+C_OM_ESS*sum(abs(P_bat_cha)+abs(P_bat_dis))+C_OM_TES*sum(abs(P_tes_cha)+abs(P_tes_dis))+C_OM_GB*(sum(P_h_GB))...
+C_OM_GT*(sum(P_e_GT)+sum(P_h_GT))+C_OM_PV*(sum(P_e_PV))+C_OM_WT*(sum(P_e_WT));
obj=C_e_buy-C_e_sell+C_h_buy-C_h_sell+C_S+C_G+C_C+C_DR-C_Ctrade+C_OM,
%% 求解设置
ops=sdpsettings('solver','cplex','verbose',2,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
ops.cplex.exportmodel='a.lp';
result=solvesdp(C,obj,ops);
if result.problem==0
% problem =0 代表求解成功
% do nothing!
else
error('求解出错');
end
%% figure
%最优电负荷
P_e_load=value(P_e_load);
Plot_e=[];
for t=1:24
Plot_e(1,t)=P_e_GT_SUP(t);
Plot_e(2,t)=P_e_buy(t);
Plot_e(3,t)=-P_bat_cha(t)+P_bat_dis(t);
Plot_e(4,t)=-P_e_sell(t);
%Plot_e(6,t)=-P_e_EC(t);
end
Plot_e=Plot_e';
figure
bar(Plot_e,'stacked');
hold on
plot(P_e_load,'r-o','LineWidth',1.0);
xlabel('时间/h');
ylabel('功率/kW');
title('最优负荷情况');
legend('GT机组','购买','电储能','售电','负荷');
P_h_GT+P_h_GB+P_h_buy+P_tes_dis==P_tes_cha+P_h_sell+P_h_load,
%最优热负荷
P_h_load=value(P_h_load);
Plot_h=[];
for t=1:24
Plot_h(1,t)=P_h_GT(t);
Plot_h(2,t)=P_h_GB(t);
Plot_h(3,t)=P_h_buy(t);
Plot_h(4,t)=-P_tes_cha(t)+P_tes_dis(t);
Plot_h(5,t)=-P_h_sell(t);
%Plot_e(6,t)=-P_e_EC(t);
end
Plot_h=Plot_h';
figure
bar(Plot_h,'stacked');
hold on
plot(P_h_load,'r-o','LineWidth',1.0);
xlabel('时间/h');
ylabel('功率/kW');
title('最优负荷情况');
legend('GT机组','GB','购买','储能','售热','负荷');
将上面这段matlab代码进行界面化。