EI顶刊复现：考虑源荷不确定性的综合能源运行调度与容量配置两阶段随机优化模型程序代码！

本程序参考EI论文《计及源荷不确定性的综合能源生产单元运行调度与容量配置两阶段随机优化》，文章中建立电转气（P2G）、碳捕集利用与封存（CCS）等关键设备及生产环节的数学模型，利用蒙特卡洛抽样随机生成新能源出力场景。程序中算例丰富、注释清晰、干货满满，创新性很高！下面对文章和程序做简要介绍！

文章亮点：

1. 源荷不确定性；

2. 蒙特卡洛抽样法生成新能源出力场景；

3. 综合能源系统优化调度；

4. 综合能源系统容量配置模型；

5. 两阶段随机优化方法（包含确定性优化方法程序，形成对比）；

6. 电价、气价的灵敏度分析；

文章主要工作：

对融合P2G与CCS的IEPU系统框架，建立各关键设备及生产环节数学模型，基于混合整数线性规划(mixed integer linear programming，MILP)算法，以全生命周期内经济成本最低为优化目标，考虑物料及能量平衡约束，实现典型周内各设备功率的最优逐时调度优化；以运行调度为基础，MILP运行优化算法作为子函数在蒙特卡洛抽样生成随机场景的过程中被不断调用，并返回最低经济成本的期望到上层的遗传算法中寻找最优容量配置，实现两阶段一体化的随机优化策略。最后，本文结合优化计算结果，考察两阶段随机优化和确定性优化的区别，分析讨论不确定性因素对系统优化结果的影响，并进行了灵敏度分析。

文中结果：

程序结果：

部分程序：

clc
clear
close all
%%  初始化设备参数及运行变量
addpath('..\光照强度与电负荷生成');
load('IPV1');  load('IPV2'); load('IPV3'); load('IPV4'); load('IPV5'); load('IPV6');   
load('Eload1');  load('Eload2'); load('Eload3'); load('Eload4'); load('Eload5'); load('Eload6');
It = [IPV1,IPV2,IPV3,IPV4,IPV5,IPV6];
Edemand = [Eload1,Eload2,Eload3,Eload4,Eload5,Eload6]; %电负荷
load('Nday'); %各典型日频次   
nday = [Nday(1)*ones(1,24),Nday(2)*ones(1,24),Nday(3)*ones(1,24),Nday(4)*ones(1,24),Nday(5)*ones(1,24),Nday(6)*ones(1,24)];
%% 以下，注意是把六个典型日的约束一起写，维度是24*6
T = 24*6;   
%% 1.1.1光伏设备模型 
E_PVmppt = sdpvar(1,T);  %光伏板mppt发电功率
A_PV = sdpvar(1,1);  %光伏板面积/m2
k = 0.200;  %1平方米的光伏板1000w/m2的标准电功率为200w
E_PVr = sdpvar(1,1);  %光伏板额定发电功率
ita_PV = 0.200/1000;
%文章内写了两个E_PV,有错位，本代码将其改为E_PVmppt与E_PV
E_PV = sdpvar(1,T);  %光伏板有效发电功率
E_PV_cur = sdpvar(1,T);  %弃光功率
%后文算例中出现135MW的光伏容量配置结果，那么这里的限值就算用300MW吧，即300 000kW.
E_PVr_max = 300000;  %光伏板额定发电功率.kW
%之后，这里直接将约束也写上，省的再回头来写约束了。
C=[];
C=[C, E_PVr ==  A_PV*k,
      E_PVmppt == E_PVr*ita_PV/k*It,
      E_PVmppt == E_PV + E_PV_cur,
      0<=E_PVr,E_PVr<=E_PVr_max,
      %补充
      E_PV >= 0,
      E_PV_cur >= 0,
      A_PV >= 0, 
];
%% 1.1.2 CCS 模型   
V_CO2_PGU = sdpvar(1,T);  %火电机组的二氧化碳排放量
E_PGU = sdpvar(1,T);  %火电机组发电功率
e_PGU = 0.46;   %火电机组的二氧化碳排放强度，见表1的 0.46 N.m3CO2/kW.h
ita_CCS_max = 0.65;%碳捕集效率最大值  0.65 
V_CO2_CCSmax = sdpvar(1,T); %碳捕集最大功率（体积）
V_CO2_CCS = sdpvar(1,T); %实际碳捕集功率（体积）
V_CO2_cur = sdpvar(1,T); %碳捕集功率耗散部分功率（体积）
lamdaCO2 = 0.1937;   %碳捕集功率耗电系数 kW.h/N.m3CO2
E_CCS = sdpvar(1,T); %碳捕集耗电功率
%从图5可以找出火电机组的最大出力功率180MW，最大爬坡常出现在119时刻与162时刻的正负50MW
%火电机组的最小出力功率90MW，
%表1中给出的火电机组容量为300000kW
E_PGUmax = 300000;  %kW
E_PGUmin = 90000;   %kW
dita_E_PGUmax = 50000;%kW
dita_E_PGUmin = -50000;%kW
%从图7可知CO2捕集的最大功率是23000m3每小时
%由此计算碳捕集的最大电功率为 0.1937*23000 = 4.4551e+03 kW
E_CCSmax = 4.4551e+03; %kW
C=[C,V_CO2_PGU == e_PGU*E_PGU,
     V_CO2_CCSmax == ita_CCS_max*V_CO2_PGU,
     V_CO2_CCSmax == V_CO2_CCS + V_CO2_cur,
     E_CCS == lamdaCO2*V_CO2_CCS,
     E_PGUmin<=E_PGU,E_PGU<=E_PGUmax,
     dita_E_PGUmin<=E_PGU(2:T)-E_PGU(1:T-1),E_PGU(2:T)-E_PGU(1:T-1)<=dita_E_PGUmax,
     dita_E_PGUmin<=E_PGU(1:T-1)-E_PGU(2:T),E_PGU(1:T-1)-E_PGU(2:T)<=dita_E_PGUmax,
     0<=E_CCS,E_CCS<=E_CCSmax,
     %补充
     V_CO2_CCS >= 0,
     V_CO2_cur >= 0,
     ];
%%  1.1.3 P2G 模型
lamda_H2  = 4.2; %kW.h/N.m电解槽产氢电耗系数
lamda_CH4 = 0.3;%kW.h/N.m合成甲烷电耗系数
w = 1 ;%ω为反应平衡系数   
E_EL = sdpvar(1,T); %电解槽耗电功率
E_CH4 = sdpvar(1,T); %甲烷化耗电功率
E_P2G = sdpvar(1,T); %电转气耗电功率
V_H2_EL = sdpvar(1,T); %电解槽产氢气功率
V_CH4_P2G = sdpvar(1,T); %P2G生成甲烷功率
V_CO2_P2G = sdpvar(1,T); %电转气消耗的CO2功率
%补充
V_H2_P2G = sdpvar(1,T); %电转气消耗的H2功率
%同温同压同体积所以摩尔数相同
C=[C,E_EL == lamda_H2*V_H2_EL,
     E_CH4 == lamda_CH4*V_CH4_P2G,
     E_P2G == E_EL + E_CH4,
     V_CH4_P2G == w*V_CO2_P2G,
     V_CH4_P2G == V_H2_P2G/4,
     ];
%从图5可看出电解槽的最大功率在50MW附近，最小在15MW附近，最大下坡在162时刻的20MW附近，最大上坡在10MW附近
%中国碱性电解槽投资成本将从2020年的每千瓦2000元，可以得知表2中的单位
%从表2可知，电解槽的容量是需要配置的
%因此，此处选用可配置的电解槽的最大容量为100MW，最下容量为10MW，
%上坡速率系数为10/50=0.2
%下坡速率系数为-20/50=-0.4
%最小电功率为标准容量的15/50 = 0.3
%基于此的电解槽约束如下：
E_EL_S = sdpvar(1,1);
E_EL_Smax = 100000;
E_EL_Smin = 10000;
C=[C,E_EL_Smin<=E_EL_S,E_EL_S<=E_EL_Smax,
     0.3*E_EL_S<=E_EL,E_EL<=E_EL_S,
     -0.4*E_EL_S<=E_EL(2:T)-E_EL(1:T-1),E_EL(2:T)-E_EL(1:T-1)<=0.2*E_EL_S,
     -0.4*E_EL_S<=E_EL(1:T-1)-E_EL(2:T),E_EL(1:T-1)-E_EL(2:T)<=0.2*E_EL_S,
     ];
%从图7可以看出，甲烷化耗碳量最大在4000m3，对应的电功率在1200kW
%因此，甲烷化上限可以定在2000kW
%表2可知，每kW的甲烷化设备容量投资是3000元
%如图7所示，爬坡速度比率可以取+-0.3
E_CH4_S = sdpvar(1,1);
C=[C,0<=E_CH4_S,E_CH4_S<=2000,
     0<=E_CH4,E_CH4<=E_CH4_S,
      -0.3*E_CH4_S <= E_CH4(1:T-1)-E_CH4(2:T),E_CH4(1:T-1)-E_CH4(2:T)<= 0.3*E_CH4_S,
     ];
%%  1.1.4 储氢设备模型
%储氢设备的投资成本是7.76 N.m3
%从图6可看出算例中的设备储存氢气量上限是20 000N.m3,存与放的功率上限是2500N.m3，比例是0.125
%因此，可以选择储氢容量配置的上限为 30 000 N.m3,存与放的功率比例上限是0.125
V_H2_S = sdpvar(1,1);
V_H2 = sdpvar(1,T);
V_H2_in = sdpvar(1,T);
V_H2_out = sdpvar(1,T);
Y_H2 = binvar(1,T); %存气状态标识位
M=1e6;
C=[C,0<=V_H2_S,V_H2_S<=20000,
     V_H2(2:T)==V_H2(1:T-1)+V_H2_in(1:T-1)-V_H2_out(1:T-1),
     0<=V_H2_in,V_H2_in<=0.125*V_H2_S,
     0<=V_H2_out,V_H2_out<=0.125*V_H2_S,
     0<=V_H2_in,V_H2_in<=M*Y_H2,
     0<=V_H2_out,V_H2_out<=M*(1-Y_H2),
     0<=V_H2,V_H2<=V_H2_S,
     %补充
     V_H2(1)== V_H2(T),
    ];    
%%  1.1.5 储碳设备模型
%储碳设备的投资成本是7.76 N.m3
%从图7可看出算例中的设备储存CO2量上限是20 000N.m3,存与放的功率上限是2500N.m3，比例是0.125
%因此，可以选择储CO2容量配置的上限为 30 000 N.m3,存与放的功率比例上限是0.125   
V_CO2_S = sdpvar(1,1);
V_CO2 = sdpvar(1,T);
V_CO2_in = sdpvar(1,T);
V_CO2_out = sdpvar(1,T);
Y_CO2 = binvar(1,T); %存气状态标识位
M=1e6;
C=[C,0<=V_CO2_S,V_CO2_S<=20000,
     V_CO2(2:T)==V_CO2(1:T-1)+V_CO2_in(1:T-1)-V_CO2_out(1:T-1),
     0<=V_CO2_in,V_CO2_in<=0.125*V_CO2_S,
     0<=V_CO2_out,V_CO2_out<=0.125*V_CO2_S,
     0<=V_CO2_in,V_CO2_in<=M*Y_CO2,
     0<=V_CO2_out,V_CO2_out<=M*(1-Y_CO2),
     0<=V_CO2,V_CO2<=V_CO2_S,
     %补充
     V_CO2(1)== V_CO2(T),
    ];

欢迎感兴趣的小伙伴关注下方公众号获取完整版代码，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！