👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆下载资源链接👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆
《《《《《《《《更多资源还请持续关注本专栏》》》》》》》
论文与完整源程序_电网论文源程序的博客-CSDN博客
https://blog.csdn.net/liang674027206/category_12531414.html
附带WORD解读！！！！！！！！！！！
双碳背景下，为实现低碳排放，需从低碳政策和低碳技术两个路径进行协调。为此建立了含 P2G-CCS 耦合和燃气掺氢的虚拟电厂(virtual power plant，VPP)，并提出了基于阶梯碳交易机制的VPP优化调度策略。首先，在低碳技术层面，针对P2G-CCS耦合和燃气掺氢子系统，建立了掺氢燃气轮机、掺氢燃气锅炉、两段式电转气(power to gas，P2G)和碳捕集系统(carbon capture system，CCS)的数学模型；其次，在低碳政策层面，建立了阶梯碳交易模型对系统碳排放进行约束；最后，在建模基础上，提出了以碳交易成本、购气和煤耗成本、碳封存成本、机组启停成本和弃风成本之和最低为目标函数的优化调度策略。对建立的模型线性化处理后，采用MATLAB调用CPLEX和粒子群算法进行求解，通过设置不同的情景进行对比，验证了所提模型的有效性，并分析了不同固定掺氢比、变掺氢比、不同的阶梯碳交易参数对VPP低碳性和经济性的影响。
部分代码展示：
%% 图7 变掺氢比 0%-20%
clc
clear
close all
PriM = 200; %燃煤价格500元/吨
Prico2feng = 50; %CO2封存的单位成本50元/吨
Smaxco2feng = 20; %吨
Pco2Base = 215; %200元/吨 注意单位换算
Eperco2ccs = 0.269; %kg单位碳补集电功率kW
PriG = 3.5; %天然气的价格3.5元/m3
kWind = 1; %风电放大比例（将风电调大，这样才会有弃风弃光，才会有电制氢和甲烷化的经济效益）
%% 常量初始化/变量初始化
%% 风电预测MPPT
EwindMppt = kWind*1000*[255 270 230 210 300 315 200 175 150 20 75 125 200 220 210 255 305 315 305 315 285 275 160 290 ];
Eload=1000*[205 205 210 195 180 190 190 200 255 345 375 380 375 355 350 345 360 360 357 365 365 355 355 347 ];
Hload=1000*[ 300 325 340 350 370 365 365 340 315 300 280 270 250 245 245 245 245 250 252 260 270 290 295 310 ];
T=24;
%%
m2qCH4 = 5.071e7; % 天然气热值5.071*10^7J/kg
m2qH2 = 1.4e8; %氢气热值1.4*10^8 J/kg
q2e = 3.6e6; %1kwh的电能相当于3.6e6J的热能
v2mCH4 = 0.71428; %1立方米天然气质量0.71428千克
v2mh2 = 89.9e-3; %1立方米氢气质量89.9e-3千克
v2mco2 = 1.964; %1立方米co2质量1.964千克
vch42mco2=1.9; %一立方米天然气完全燃烧后可以生成二氧化碳的重量1.9kg
%% 风机
Ewind = sdpvar(1,24);
Ewindcur = sdpvar(1,24);
%% 燃气轮机
ECHPmax = 1000*350; %燃气轮机电出力上限kW
ECHPmin = 0;
HCHPmax = 1000*300; %燃气轮机热出力上限kW
HCHPmin = 0;
ditaEHCHPmax = 1000*150; %燃气轮机总功率爬坡上限kW
ditaEHCHPmin = -1000*150;
nHCHP =0.4; %可利用热能站混合燃气总热能的比例
nECHP =0.35; %可利用电能站混合燃气总热能的比例
ECHP = sdpvar(1,24); %燃气轮机电出力
HCHP = sdpvar(1,24); %燃气轮机热出力
EHCHP = sdpvar(1,24); %燃气轮机电热总出力
mco2CHP = sdpvar(1,24);%燃气轮机碳排放质量kg
vco2CHP = sdpvar(1,24);%燃气轮机碳排放体积m3
mCH4CHP = sdpvar(1,24);%质量kg
mh2CHP = sdpvar(1,24);%质量kg
vch4CHP = sdpvar(1,24);%体积：标准立方米
vH2CHP = sdpvar(1,24);%体积：标准立方米 %掺氢比例在10%-20% %改为0%-20%
%% 燃气锅炉
HGBmax = 1000*80; %kW 燃气锅炉热出力上限
HGBmin = 0;
ditaHGBmax = 1000*25; %kW 燃气锅炉爬坡
ditaHGBmin =-1000*25; %kW
nHGB = 0.92; %混合燃气热量 到 可利用热量 的转化系数
HGB = sdpvar(1,24);
% qCH4GB = sdpvar(1,24); %GB天然气热量J
% qH2GB = sdpvar(1,24); %GB氢气热量J
mCH4GB = sdpvar(1,24); %GB天然气质量kg
vco2CH4 = sdpvar(1,24);
mh2GB = sdpvar(1,24); %GB氢气质量kg
vch4GB = sdpvar(1,24); %GB天然气体积m3
vH2GB = sdpvar(1,24); %GB氢气体积m3 %掺氢比例在2%-20% %改为0%-20%
mco2GB = sdpvar(1,24); %GB二氧化碳质量kg
vco2GB = sdpvar(1,24); %燃气轮机碳排放体积m3
%% 电加热锅炉
EEBmax = 1000*40;%kW
EEBmin = 0;
ditaEEBmax = 1000*10; %kW
ditaEEBmin = -1000*10; %kW
nEEB = 0.9; %这个直接就是kWh电 到 kWh热，很简单
EEB = sdpvar(1,24); %kW电加热锅炉的耗电
HEB = sdpvar(1,24); %kW电加热锅炉的产热
%% 火电机组
EMmax = 1000*162;%kW %火电机组最大发电功率 kW
EMmin = 1000*45;%kW %火电机组最小发电功率 kW
ditaEMmax = 1000*100;%kW %爬坡
ditaEMmin =-1000*100;%kW %爬坡
EM = sdpvar(1,24); %火电机组发电功率kW
YEM = binvar(1,24); %火电机组启停变量（1是运行，0是停止）
YEMqi = binvar(1,24); %0变1
YEMting = binvar(1,24); %1变0
mco2EM = sdpvar(1,24); %火电机组碳排放量kg
vco2EM = sdpvar(1,24); %火电机组碳排放量m3
mEM = sdpvar(1,24); %煤耗 kg
%% 电转气
nP2H = 0.85; %电制氢效率 电能J转氢气热能J
EP2Hmax = 1000*120;%kW 电转气耗电功率上限
EP2Hmin = 0; %
nCH4 = 0.7; %甲烷化效率 氢气热能J转天然气热能J
EP2H = sdpvar(1,24); %耗电kW
mh2P2H =sdpvar(1,24); %制氢kg
mh2CH4 =sdpvar(1,24); %甲烷化耗氢kg
mch4CH4 = sdpvar(1,24); %甲烷化制取天然气 kg
vch4CH4 = sdpvar(1,24); %甲烷化制取天然气体积m3
mco2CH4 = sdpvar(1,24); %甲烷化吸收co2 kg
%% 碳补集系统
ECCSmin = 0;
ECCSmax = 1000*150; %kW %碳补集耗电功率上限
VPFmax = 29200; %m3 %碳补集富液体积上限
VPFmin = 0;
VPFstart = 14600;%m3
ECCS = sdpvar(1,24); %碳补集耗电功率kW
mco2CCSin = sdpvar(1,24); %吸收co2质量kg
mco2CCSout = sdpvar(1,24); %释放co2质量kg
vco2CCSin = sdpvar(1,24); %吸收co2体积m3
vco2CCSout = sdpvar(1,24); %释放co2体积m3
VF = sdpvar(1,24);
VP = sdpvar(1,24);
YCCSin = binvar(1,24);
%% 储热
HSSsocmax = 1000*120; %kWh 容量
HSSsocmin = 1000*20; %kWh
HSSmaxC = 1000*15; %kW 功率
HSSdmaxD = 1000*15; %kW
nHSS = 0.95;
sunHSS=0.01;
HSSsocstart=1000*60; %kWh 初始容量
HSSsoc=sdpvar(1,24); %kW
HSSC=sdpvar(1,24); %充电功率
HSSD=sdpvar(1,24); %放电功率
Y_HSS = binvar(1,T);
%% 储电设备
ESSsocmax = 1000*60; %kWh
ESSsocmin = 1000*10; %kWh
ESSmaxC = 1000*5; %kW
ESSmaxD = 1000*5; %kW
nESS = 0.9;
sunESS=0.001;
ESSsocstart=1000*30; %kWh
ESSsoc=sdpvar(1,24); %kWh
ESSC=sdpvar(1,24);
ESSD=sdpvar(1,24);
Y_ESS = binvar(1,T);
%% 天然气管网
VGgrid = sdpvar(1,24);
%% C02封存
mco2Storage = sdpvar(1,24);
M =1e8;
%阶梯碳交易
Lco2 = 50; %碳交易区间宽度
ditaPco2 = 0.25; %负碳区间的补偿上升率
thitaPco2 = 0.25; %正碳区间的价格增长率
ejy = sdpvar(1,24); %碳交易的功率
fco2 = sdpvar(1,24); %碳交易各时刻的费用
% -Pco2Base*(2+3*ditaPco2)*Lco2+Pco2Base*(1+3*ditaPco2)*(ejy+2*Lco2), -M<=ejy<=-2*Lco2
% -Pco2Base*(1+ditaPco2)*Lco2+Pco2Base*(1+2*ditaPco2)*(ejy+Lco2), -2*Lco2 <=ejy<=Lco2
% Pco2Base*(1+ditaPco2)*ejy, -Lco2 <=ejy<=Lco2<=0
% fco2 = Pco2Base*ejy, 0<=ejy<=Lco2
% Pco2Base*Lco2+Pco2Base*(1+thitaPco2)*(ejy-Lco2), Lco2<=ejy<=2*Lco2
% Pco2Base*(2+thitaPco2)*Lco2+Pco2Base*(1+2*thitaPco2)*(ejy-2*Lco2), 2*Lco2<=ejy<=M
dco2 = binvar(6,24);
zco2 = sdpvar(6,24); %公式（B-3）的y就等同于fco2
%% 设备建模约束
%% 风机
C=[ ];
C=[C,0<=Ewind,Ewind<=EwindMppt,
0<=Ewindcur,Ewindcur<=EwindMppt,
Ewind+Ewindcur==EwindMppt,
];
%C02封存
C=[C, 0<=mco2Storage,mco2Storage<=Smaxco2feng*1000 ];
%燃气轮机建模约束
C=[C, ECHPmin<=ECHP,ECHP<=ECHPmax, %公式(27)
HCHPmin<=HCHP,HCHP<=HCHPmax,
ECHP(2:24)+HCHP(2:24)-ECHP(1:23)-HCHP(1:23)>= ditaEHCHPmin,
ECHP(2:24)+HCHP(2:24)-ECHP(1:23)-HCHP(1:23)<= ditaEHCHPmax,
(m2qCH4*mCH4CHP+m2qH2*mh2CHP)/q2e*nHCHP == HCHP,
(m2qCH4*mCH4CHP+m2qH2*mh2CHP)/q2e*nECHP == ECHP,
mCH4CHP == v2mCH4*vch4CHP,
mh2CHP == v2mh2*vH2CHP,
mco2CHP == vch42mco2*vch4CHP,
mco2CHP == v2mco2*vco2CHP,
vH2CHP>=0,
vH2CHP<=(vH2CHP+vch4CHP)*0.2, %掺氢体积比，最大20%
];
%燃气锅炉建模约束
C=[C, HGBmin<=HGB,HGB<=HGBmax,
ditaHGBmin<=HGB(2:24)-HGB(1:23),HGB(2:24)-HGB(1:23)<=ditaHGBmax,
HGB==nHGB*vch4GB*36e6/3.6e6+nHGB*mh2GB*m2qH2/3.6e6,
mco2GB == vch42mco2*vch4GB,
mco2GB == v2mco2*vco2GB,
mh2GB == v2mh2*vH2GB,
vH2GB>=0,
vH2GB<=(vH2GB+vch4GB)*0.2, %掺氢体积比，最大20%
];
%电加热锅炉
C=[C, EEBmin<=EEB,EEB<=EEBmax,
ditaEEBmin<=EEB(2:24)-EEB(1:23),EEB(2:24)-EEB(1:23)<=ditaEEBmax,
HEB==EEB*nEEB,
];
%火电机组
C=[C,
EMmin*YEM<=EM,EM<=EMmax*YEM,
ditaEMmin<=EM(2:24)-EM(1:23),EM(2:24)-EM(1:23)<=ditaEMmax,
mEM==0.320*EM,%一度电约等于320g标准煤
mco2EM == 0.7976*EM, %1 kW·h电折算标准煤0.32 kg,co2排放0.7976kg
mco2EM == v2mco2*vco2EM,
YEMqi(2:24)-YEMting(2:24)==YEM(2:24)-YEM(1:23),
YEMqi(1)==0,
YEMting(1)==0,
];
%电转气
C=[C,
EP2Hmin<=EP2H,EP2H<=EP2Hmax,
m2qH2*mh2P2H==nP2H*EP2H*3.6e6,
nCH4*m2qH2*mh2CH4==m2qCH4*mch4CH4,
mco2CH4==44/16*mch4CH4, %实际应该比44/16这个系数略小，因为甲烷化生成的天然气里有氢气，
vch4CH4 == mch4CH4/v2mCH4,
mco2CH4 == v2mco2*vco2CH4,
mco2CH4 == mco2CCSout/44*12,
mh2CH4>=0,
];
%碳补集系统
C=[C,
ECCSmin<=ECCS,ECCS<=ECCSmax,
mco2CCSin*Eperco2ccs == ECCS, %co2补给1kg耗电0.269kWh
ECCSmin<=ECCS,ECCS<=ECCSmax,
mco2CCSin == v2mco2*vco2CCSin,
0<=vco2CCSout,vco2CCSout<=ECCSmax/Eperco2ccs,
VF(1)==VPFstart+vco2CCSin(1)/25-vco2CCSout(1)/25,%单位体积的富液能够吸收其 25 倍体积的 CO2
VF(2:24)==VF(1:23)+vco2CCSin(2:24)/25-vco2CCSout(2:24)/25,
VF(24)==VPFstart,
VP+VF == VPFmax,
mco2CCSout == v2mco2*vco2CCSout,
0<=vco2CCSin,vco2CCSin<=YCCSin*M, %新添加
0<=vco2CCSout,vco2CCSout<=(1-YCCSin)*M,
0<=mco2CCSin,mco2CCSin<=(mco2CHP+mco2GB+mco2EM)*0.9,
];
%储热
C=[C,HSSsocmin<=HSSsoc,HSSsoc<=HSSsocmax,
HSSsoc(2:T)==HSSsoc(1:T-1)*(1-sunHSS)+nHSS*HSSC(2:T)-HSSD(2:T)/nHSS,
HSSsoc(1)==HSSsocstart*(1-sunHSS)+nHSS*HSSC(1)-HSSD(1)/nHSS,
HSSsoc(24)== HSSsocstart,
0<=HSSC,HSSC<=HSSmaxC,
0<=HSSD,HSSD<=HSSmaxC,
0<=HSSC,HSSC<=M*Y_HSS,
0<=HSSD,HSSD<=M*(1-Y_HSS),
];
%储电
C=[C,ESSsocmin<=ESSsoc,ESSsoc<=ESSsocmax,
ESSsoc(2:T)==ESSsoc(1:T-1)*(1-sunESS)+nESS*ESSC(2:T)-ESSD(2:T)/nESS,
ESSsoc(1)==ESSsocstart*(1-sunESS)+nESS*ESSC(1)-ESSD(1)/nESS,
ESSsoc(24)== ESSsocstart,
0<=ESSC,ESSC<=ESSmaxC,
0<=ESSD,ESSD<=ESSmaxD,
0<=ESSC,ESSC<=M*Y_ESS,
0<=ESSD,ESSD<=M*(1-Y_ESS),
];
H2SSsoc = sdpvar(1,24);
H2SSC = sdpvar(1,24);
H2SSD = sdpvar(1,24);
Y_H2SS = binvar(1,24);
H2SSsocmin = 1000; %kg
H2SSsocmax = 5000; %100*1000/5*89.3e-3kg
nH2SS = 0.98;
H2SSmaxC = 0.25*H2SSsocmax;
H2SSmaxD = 0.25*H2SSsocmax;
H2SSsocstart = 0.5*H2SSsocmax;
%储氢
C=[C,H2SSsocmin<=H2SSsoc,H2SSsoc<=H2SSsocmax,
H2SSsoc(2:T)==H2SSsoc(1:T-1)+nH2SS*H2SSC(2:T)-H2SSD(2:T)/nH2SS,
H2SSsoc(1)==H2SSsocstart+nH2SS*H2SSC(1)-H2SSD(1)/nH2SS,
H2SSsoc(24)== H2SSsocstart,
0<=H2SSC,H2SSC<=H2SSmaxC,
0<=H2SSD,H2SSD<=H2SSmaxD,
0<=H2SSC,H2SSC<=M*Y_H2SS,
0<=H2SSD,H2SSD<=M*(1-Y_H2SS),
];
%% 平衡约束
C=[C,ECHP+EM+Ewind+ESSD==EP2H+EEB+ESSC+ECCS+Eload ]; %电平衡 要注意没有外网交互
C=[C,HCHP+HGB+HEB+HSSD==HSSC+Hload ];
C=[C,VGgrid+vch4CH4== vch4CHP+vch4GB ];
C=[C,mh2P2H+H2SSD==mh2CH4+mh2CHP+mh2GB+H2SSC ];
%% 阶梯碳交易
%碳交易量的计算
mco2_shiji = sdpvar(1,24);
mco2_peie = sdpvar(1,24);
lamdaECHP = 0.2;
lamdaHCHP = 0.2;
lamdaEM = 0.2;
lamdaGB = 0.2;
C=[C,mco2_shiji == mco2CHP+mco2GB+mco2EM-mco2CCSin,
mco2_shiji >= 0,
mco2_peie == lamdaECHP*ECHP+lamdaHCHP*HCHP+lamdaEM*EM+lamdaGB*HGB, ];
%碳交易分段函数的线性化
a=repmat([Pco2Base*(1+3*ditaPco2),Pco2Base*(1+2*ditaPco2),Pco2Base*(1+ditaPco2),Pco2Base,Pco2Base*(1+thitaPco2),Pco2Base*(1+2*thitaPco2)]',1,24);
b=repmat([-Pco2Base*(2+3*ditaPco2)*Lco2+Pco2Base*(1+3*ditaPco2)*2*Lco2, -Pco2Base*(1+ditaPco2)*Lco2+Pco2Base*(1+2*ditaPco2)*Lco2,0,0,Pco2Base*Lco2+Pco2Base*(1+thitaPco2)*(-Lco2),Pco2Base*(2+thitaPco2)*Lco2+Pco2Base*(1+2*thitaPco2)*(-2*Lco2)]',1,24);
% -Pco2Base*(2+3*ditaPco2)*Lco2+Pco2Base*(1+3*ditaPco2)*(ejy+2*Lco2), -M<=ejy<=-2*Lco2
% -Pco2Base*(1+ditaPco2)*Lco2+Pco2Base*(1+2*ditaPco2)*(ejy+Lco2), -2*Lco2 <=ejy<=-Lco2
% Pco2Base*(1+ditaPco2)*ejy, -Lco2 <=ejy<=Lco2<=0
% fco2 = Pco2Base*ejy, 0<=ejy<=Lco2
% Pco2Base*Lco2+Pco2Base*(1+thitaPco2)*(ejy-Lco2), Lco2<=ejy<=2*Lco2
% Pco2Base*(2+thitaPco2)*Lco2+Pco2Base*(1+2*thitaPco2)*(ejy-2*Lco2), 2*Lco2<=ejy<=M
boundco2l = repmat([ -M,-2*Lco2,-Lco2,0,Lco2,2*Lco2 ]',1,24 );
boundco2u = repmat([ -2*Lco2,-Lco2,0,Lco2,2*Lco2,M ]',1,24 );
C=[C,fco2==sum(a.*zco2+b.*dco2),
(mco2_shiji-mco2_peie)/1000==sum(zco2),
sum(dco2)==1,
dco2.*boundco2l<=zco2,zco2<=dco2.*boundco2u, ];
%% 目标函数
F = sdpvar(1,1); %总成本
Fco2 = sdpvar(1,1); %1）碳交易与碳封存成本。
C=[C,Fco2==sum(fco2)+Prico2feng*sum(mco2Storage/1000)]; %100元封每吨CO2
Fth = sdpvar(1,1); %2）燃煤机组启停+煤耗成本+运维成本
C=[C, Fth== 15e4*sum(YEMqi+YEMting)+PriM*sum(mEM/1000)+0.1*sum(EM)]; %火电机组出力的平方暂且不处理
Fcur = sdpvar(1,1); %弃风成本
C=[C, Fcur == 0.6*sum(Ewindcur) ];
Fbuy = sdpvar(1,1); %购气成本
C=[C, Fbuy == PriG*sum(VGgrid) ];
C=[C, F == Fco2+Fth+Fcur+Fbuy ];
%% 调用cplex求解器求解
% ops=sdpsettings('verbose', 2, 'solver', 'cplex');
ops = sdpsettings('solver', 'cplex+', 'verbose', 2, 'debug', 1);
ops.cplex.MIPGap =0.0001;
ops.cplex.TimeLimit = 200;
ops.cplex.MIPFocus=3;
ops.cplex.Presolve=2;
ops.cplex.DualReductions = 0;
%进行求解计算
result=optimize(C,F,ops);
if result.problem==0
fprintf('求解成功 ！！！\n');
else
error('求解出错！请查找错误来源');
end
%% 绘图
mh2CH4=value(mh2CH4);mh2CHP=value(mh2CHP);mh2GB=value(mh2GB);H2SSC=value(H2SSC);
figure
% BianH2xiaohao = mh2CH4+mh2CHP+mh2GB+H2SSC;
BianH2xiaohao = mh2CH4+mh2CHP+mh2GB;
plot(BianH2xiaohao/1000);
xlabel('时间/h');
ylabel('功率/(吨/h)');
title('掺氢耗氢功率图');
legend('变掺氢耗氢');
xlim([0 30 ]);
set(gca,'FontSize',12 );
save('BianH2xiaohao.mat','BianH2xiaohao');
效果展示：
76号资源-复现源程序：论文可在知网下载《基于阶梯碳交易的含P2G-CCS耦合和燃气掺氢的虚拟电厂优化调度》本人博客有解读资源-CSDN文库https://download.csdn.net/download/LIANG674027206/89139682👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆下载资源链接👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆👆
《《《《《《《《更多资源还请持续关注本专栏》》》》》》》
论文与完整源程序_电网论文源程序的博客-CSDN博客https://blog.csdn.net/liang674027206/category_12531414.html
975
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
