说句实在话看到标题里带“SCI一区复现”的时候我第一反应是又双叒叕一篇综合能源系统运行优化。但真把这篇Energy期刊的文章拆开读完之后我发现它跟很多“挂羊头卖狗肉”的IES文章不一样一是它把P2G和碳捕集设备做成了物理耦合的碳循环而不是两个孤立模块二是它用了epsilon约束算法求解碳排放成本和运维成本的双目标优化问题不是简单做个加权和。这两点放在一起足够让一个做调度优化的人忙活好几周。这篇文章就是我复现这个工作的全过程记录包含模型理解、数学建模、YALMIPGurobi代码落地以及我踩过的一些很绕的坑适合正在做综合能源优化、碳捕集/P2G方向、或者刚想学epsilon约束法的同学参考。1. 为什么做这个复现以及复现之前要读懂哪些东西做IES方向的人应该都有体会P2G在模型里通常被当成“电转氢再转气”的能源路由器碳捕集设备则是一个降低净排放的装置。论文里把这两个设备同时写进模型的不少但真正把P2G甲烷化需要的CO2和碳捕集设备捕集到的CO2连成一条流约束的其实不多。这篇论文让我印象最深的地方就是它把“碳循环”做成了物理约束而不是简单地在排放系数上打个折。加上它用的是epsilon约束法而不是加权和法这一点对我来说吸引力很大因为加权和法在处理混合整数非线性或非凸前沿时确实容易漏解。复现之前我花了两三天时间把论文的模型部分读了不止一遍。这里建议所有打算做复现的人动手写代码之前先把下面四件事搞清楚否则后面一定会返工。第一系统拓扑是什么。这篇论文的拓扑大致是燃气轮机供应电和热余热锅炉回收热量燃气锅炉作为热力备用电解槽消耗电能制氢氢气进入甲烷化反应器与CO2反应生成合成天然气合成天然气可以注入气网或者回供给燃气轮机碳捕集设备捕集烟气中的CO2一部分CO2送去甲烷化剩余部分封存或外售。这个拓扑决定了所有变量的耦合关系拓扑错了后面的约束全是错的。第二设备建模假设。论文里有没有把电解槽的效率当成常数碳捕集设备的捕集率是固定比例还是随负荷变化燃气轮机热电比是定值还是变量这些假设直接决定你的模型是LP、QP还是MILP。我复现时按最常见的做法把电解槽和甲烷化效率设为常数碳捕集效率也设为固定比例这样模型规模控制在可接受范围。第三时间分辨率。论文里是24小时典型日优化还是8760小时全年优化。复现典型日优化就够了24个时段步长1小时。如果你想做全年动态省省吧单日模型跑通再扩展别一上来就把自己累死。第四双目标怎么处理。我在这篇论文上做的“增加”就是基于epsilon约束法重新构造求解流程。这个方法需要先分别求两个单目标的最小值然后确定epsilon的取值区间再逐步收紧约束求帕累托前沿。论文里写得比较精简实际操作中边界求不准、区间取不好前沿就会断。还有一点要提醒复现一区论文不是照抄代码。论文一般不公开源码即使公开很多参数也是基于特定算例的。你需要做的是重构模型把论文里未明确的参数用合理典型值补全。比如设备的效率、上下限、碳排放因子、碳价、分时电价这些论文给不全很正常这时候就用行业常用值或者文献典型值并在记录里标注哪些参数是假设的。2. 先把能流图画清楚P2G、碳捕集和热电联供的耦合逻辑代码写不出来的时候十有八九是能流图没想明白。我复现的第一周基本上没碰电脑全在纸上画系统能量流和碳流把每个设备的输入输出列成一张表再对着表检查模型有没有漏掉耦合关系。2.1 热电联供单元燃气轮机和余热锅炉热电联供单元是整个系统的电热源头。燃气轮机消耗天然气发电同时产生高温烟气余热锅炉回收这部分热量供热。燃气锅炉作为备用热源在热负荷高峰或者燃气轮机出力受限时补充供热。在建模时要抓住三个关键参数发电效率、热电比、烟气排放因子。发电效率决定同样天然气能发多少电热电比决定发多少电的同时能供多少热排放因子决定烧同样天然气产生多少CO2。这三个参数抓准了CHP单元的能量流和碳流就清楚了。燃气轮机的运行区间也要注意。它不是一个简单的点模型而是有最小技术出力、最大出力、爬坡速率限制的。尤其是有没有0-1启停变量直接决定了模型是LP还是MILP。论文里如果没有明确说连续运行我建议用连续变量加上下限约束就行先跑通再考虑加启停逻辑。2.2 P2G的双段结构电解槽与甲烷化反应器P2G不是一个设备而是“电解槽甲烷化反应器”两个环节串联。电解槽消耗电能和水产出氢气氢气进入甲烷化反应器与CO2反应生成甲烷合成天然气。两个环节各有各的效率建模时必须分开写。实际建模时电解槽的输入是电功率输出是氢气流量。这里要用氢气的热值换算不能直接用功率乘以一个系数就完事。甲烷化反应需要CO2和H2摩尔配比大约是1比4也就是说想生成多少合成天然气必须同时保证有足够的氢气和CO2。这就是我在后文说的“假绿氢”坑的根源——如果漏写了CO2来源约束模型会让甲烷化凭空消耗CO2碳排直接被严重低估。P2G产出的合成天然气可以注入天然气网络也可以直接回供给燃气轮机和燃气锅炉使用。这一步其实是把电力系统、天然气系统和碳捕集系统连接起来的关键桥梁。也就是说P2G不仅是电负荷同时也是气源它消耗电和CO2产出可用的天然气。这个双角色属性是所有耦合约束的出发点。2.3 碳捕集设备如何嵌进能量流碳捕集设备捕集燃气轮机和锅炉烟气中的CO2捕集过程需要消耗电和热胺法吸收需要再沸器热负荷不过很多简化模型只考虑电耗。捕集到的CO2一部分作为原料送给P2G甲烷化另一部分可以封存或外售都可以算作减少系统净排放。我画能流图时最关注的是这条“碳回路”燃气轮机燃烧天然气发电供热产生CO2碳捕集设备从烟气中捕集CO2一部分CO2送入甲烷化反应器与电解槽产出的H2反应生成合成天然气合成天然气又回到燃气轮机那里。这样就形成一个“天然气→CO2→甲烷化原料→合成天然气→天然气”的碳循环CCS为P2G提供碳源P2G为CHP提供燃料两者互为补充。为了把各设备的输入输出关系看清楚我复现前列了这样一张表建议你也做一张设备/单元输入输出关键耦合关系燃气轮机天然气、空气电、热、烟气CO2出力产生碳排放热电比约束余热锅炉高温烟气热回收燃气轮机余热燃气锅炉天然气、空气热、烟气CO2热备用产生碳排放电解槽电、水氢气电负荷P2G第一段甲烷化反应器氢气、CO2合成天然气P2G第二段消耗CO2碳捕集设备烟气CO2、电浓缩CO2捕集率约束为P2G供碳这张表是我后面写代码时的“宪法”每个变量必须能在这张表里找到归宿每个约束必须能对应到一条物理关系上。如果代码写完发现某个变量在表中找不到明确位置那大概率是建模出了问题。3. 双目标优化建模的难点碳排放成本怎么进目标函数不少文章的双目标优化其实就是“成本最小化”再加一个“碳排放最小化”两个目标各乘一个权重加起来。这篇论文不一样的地方在于用epsilon约束法把其中一个目标变成约束去枚举解集。但在做epsilon约束之前首先得把两个目标函数本身写清楚尤其是碳排放成本绝不是“总购电量乘个因子”这么简单。3.1 运维成本和碳排放成本的目标函数写法运维成本我按常见做法拆成四块购电成本、购气成本、设备运行维护成本、可能还有启停成本。写成数学形式大致是C_op Σ_t [ c_buy_e(t) * P_buy(t) c_buy_g(t) * G_buy(t) ] Σ_t Σ_i [ c_mt_i * P_i(t) ]其中P_buy是购电功率G_buy是购气流量c_mt是设备单位出力维护成本P_i是设备i的出力。P2G的电解槽耗电、碳捕集设备的耗电都体现在购电功率里如果系统内有新能源则需要在电功率平衡里扣掉。碳排放成本稍微复杂一点。系统的总的二氧化碳排放不能只看末端排放要把煤气链条上的碳流算进去E_total E_grid E_gt E_gb - E_cap - E_p2gE_grid是外购电对应的碳排放按电网平均排放因子估算E_gt和E_gb分别是燃气轮机和燃气锅炉燃烧天然气产生的CO2E_cap是碳捕集设备捕集到的CO2量E_p2g是P2G甲烷化过程中消耗掉的CO2量。这里要注意捕集的CO2如果封存了和送到P2G去利用了对净排放的影响是一样的都是让这部分碳不排入大气。所以碳排放成本写为C_em ρ_CO2 * E_totalρ_CO2就是碳价单位一般是元/吨或$/吨需要和CO2量的单位对齐。3.2 约束条件里最容易出错的碳流平衡目标函数写完之后最关键的约束就是碳流平衡。我在这个环节吃过亏所以专门多写几句。首先是捕集量约束碳捕集设备捕集的CO2量不能超过烟气理论可捕集量乘上捕集效率也就是E_cap(t) ≤ η_cap * (E_gt(t) E_gb(t))这里的η_cap就是碳捕集率一般在0.8到0.9之间。然后是P2G利用碳约束这里也是我最容易漏的地方P2G甲烷化消耗的CO2必须来自碳捕集设备提供的CO2也就是说模型里必须有一个约束E_p2g(t) ≤ E_cap(t)甚至可以写成等式假设所有捕集到的CO2都优先用于P2G剩余部分才封存。如果不写这条约束模型会让P2G凭空消耗CO2相当于系统从空气中免费取碳碳排放成本会被严重低估结果和物理实际完全不符。除了碳流还有常规的电、热、气功率平衡约束。电功率平衡购电功率加燃气轮机发电功率减去电解槽耗电和碳捕集设备耗电要等于电负荷。热功率平衡余热锅炉供热加燃气锅炉供热等于热负荷。气功率平衡购气量加P2G产气量要等于燃气轮机和燃气锅炉的耗气量。这三条平衡约束看起来平平无奇但有一个容易被忽略的点P2G产气量如果直接进入气网那么购气量和耗气量之间是可能不同时反号的需要理清是“注入气网”还是“直接供给人端”。一般简化模型是直接供给燃气轮机气功率平衡写成购气产气用气这样就能让P2G替代一部分外购天然气。设备出力约束方面燃气轮机和燃气锅炉有上下限约束电解槽有额定功率约束碳捕集设备有处理能力上限。如果有爬坡约束需要写成相邻时段出力的差值限制。这些相对常规但很繁琐少写一个都可能让结果在某个时段失真。3.3 为什么加权和法在这里容易翻车我知道很多人一看到双目标优化第一反应就是加权和法把两个目标线性加权成一个目标。这个方法简单但在这个场景下有硬伤。因为系统里有设备上下限约束、可能还有0-1启停变量这个问题的帕累托前沿不一定是凸的。加权和法本质上是在目标空间里朝一个固定方向做线性滑动它只能找到帕累托前沿的凸包上的点那些落在凹区域里的解加权和是扫不到的。epsilon约束法的做法完全不同。它把其中一个目标比如运维成本转成不等式约束 C_op ≤ epsilon然后去最小化另一个目标碳排放成本。每给定一个epsilon就得到一个在“运维成本不超过e”条件下的碳排放最低解。不断改变epsilon就能枚举出前沿上的所有点包括非凸部分。这正是这篇论文选择它的原因也是我复现时坚持用它而不是简单加权的原因。在实际实现中epsilon约束法的另一个好处是直观你可以控制枚举多少个点每一个点都有明确的“成本上限”解释向别人解释结果时也好沟通。4. Matlab求解框架与epsilon约束法的落地代码模型理清楚之后就是代码落地了。我用的组合是Matlab YALMIP Gurobi这也是目前做这种混合整数线性规划最顺手的组合之一。4.1 求解器选型YALMIPGurobi的理由为什么不直接用Matlab自带的linprog或intlinprog因为linprog只能处理线性规划如果模型里涉及0-1启停变量就必须用混合整数线性规划求解器intlinprog虽然能用但在大规模多时段、多设备、带大量约束的情况下求解速度和稳定性都不如商业求解器。YALMIP是一个建模层工具它把符号变量和约束写成求解器能理解的形式再交给Gurobi求解整个过程非常自然。我选择YALMIPGurobi的另一个原因是Gurobi求解MILP的速度足够快而且YALMIP提供了方便的optimize接口切换求解器也很容易。节点数量几十上百时没问题24时段模型通常几秒就能出解。当然你需要先装好Gurobi并申请学术license这一步可能会劝退一些人但为了复现效率和稳定性值得。4.2 变量表、约束组和代价函数怎么写写Matlab代码之前我先把每个变量定义好维度通常都是1×T的矩阵T是时段数24。我的变量命名表大致如下变量名含义维度E_gt燃气轮机发电功率1×TH_gt燃气轮机供热功率1×TH_boiler燃气锅炉供热功率1×TP_buy购电功率1×TG_buy购气流量1×TP_el电解槽输入电功率1×TH2_p2g电解槽产氢量1×TSNG_p2g甲烷化产气量1×TCO2_cap碳捕集捕集量1×TCO2_p2gP2G甲烷化消耗的CO2量1×TE_total系统净排放量1×T实际代码里我会把YALMIP的sdpvar变量定义写清楚然后逐条添加约束。这里给一个简化版的约束构建示意方便你理解写法T 24; E_gt sdpvar(1, T, full); H_gt sdpvar(1, T, full); P_el sdpvar(1, T, full); P_buy sdpvar(1, T, full); G_buy sdpvar(1, T, full); CO2_cap sdpvar(1, T, full); CO2_p2g sdpvar(1, T, full); SNG_p2g sdpvar(1, T, full); C []; % 电功率平衡购电燃机发电 负荷电解槽耗电碳捕集耗电 C [C, P_buy E_gt P_load P_el P_ccs]; % 热功率平衡 C [C, H_gt H_boiler H_load]; % 气功率平衡购气P2G产气 燃机耗气锅炉耗气 C [C, G_buy SNG_p2g Gas_gt Gas_boiler]; % 碳捕集量不超过可捕集量 C [C, CO2_cap eta_cap * (Emission_gt Emission_boiler)]; % P2G消耗的CO2不超过捕集量 C [C, CO2_p2g CO2_cap]; % 设备上下限 C [C, 0 E_gt E_gt_max, 0 H_gt H_gt_max]; % 目标函数 op_cost sum(P_buy .* c_buy_e G_buy .* c_buy_g) ... sum(c_mt_gt * (E_gt H_gt) c_mt_el * P_el); em_cost rho_co2 * (sum(E_grid_co2) sum(E_gt_co2) ... sum(E_boiler_co2) - sum(CO2_cap) - sum(CO2_p2g));这只是最简单的一版框架。实际论文里还会有更多约束比如P2G产氢量和产气量之间的效率转换关系、甲烷化反应的物料配比关系、碳捕集设备的能耗特性等。这些约束要一个不少地加进去模型才算完整。4.3 epsilon约束法的循环骨架epsilon约束法的核心代码其实不复杂逻辑就是“先求边界再逐步扫描”。我用的是这个流程第一步分别求两个单目标的最优解得到帕累托前沿的两个端点。这里要注意求运维成本最小时碳排放成本完全自由求碳排放成本最小时运维成本完全自由。得到的结果是运维最小解对应的碳排放值以及碳排放最小解对应的运维值。第二步把碳排放成本作为主目标运维成本作为约束。epsilon的取值范围从运维成本最小值到碳排放最小解对应的运维成本最大值之间选取。将区间等分逐步收紧epsilon每取一个epsilon就求解一次单目标优化得到一个帕累托解。第三步把这些解画在二维坐标系里就是帕累托前沿。代码框架如下% 第1步求单目标边界 optimize(C, op_cost); op_min value(op_cost); op_at_em_min_op value(em_cost); % 运维最优解对应的碳排 optimize(C, em_cost); em_min value(em_cost); op_at_em_min value(op_cost); % 碳排最优解对应的运维成本 % 第2步确定epsilon区间并扫描 eps_max op_at_em_min; eps_min op_min; N 30; pf_op zeros(1, N); pf_em zeros(1, N); for k 1:N eps eps_max - (eps_max - eps_min) * k / N; C_eps [C, op_cost eps]; diagnostics optimize(C_eps, em_cost); if diagnostics.problem 0 pf_op(k) value(op_cost); pf_em(k) value(em_cost); end end注意一点如果扫描点落在不可行区域diagnostics.problem会非零这时需要跳过该点否则会得到一个无效解或者NaN。我一开始没有处理这个情况结果在画图的时候前沿出现了断点后面调整网格才解决。关于扫描点数N我建议先设20到30个点跑一遍看轮廓。如果某些区间变化剧烈比如前沿在某个epsilon附近突然拐弯就在那个区间附近加密比如再加5到10个点。不要一上来就N100后面我在踩坑部分会说为什么。5. 复现踩坑记录从断裂前沿到漏掉的二氧化碳流约束复现工作真正花时间的不是搭框架而是调试各种意想不到的问题。这里把我踩过的一些比较有代表性的坑写出来给后来人省点时间。5.1 坑1epsilon区间取错帕累托前沿第一段就断了我第一次跑epsilon约束法的时候直接把运维成本最小的单目标解当成epsilon下限碳排放成本最小的单目标解对应的运维成本当成epsilon上限结果画出来的前沿中间缺了一整块。检查之后发现求碳排放最小时我确实得到了一个可行解但这个解对应的运维成本并不是“能达到的最大运维成本”而只是在那一组最优解里的一个值。由于碳排放成本函数可能存在多个最优解有的解运维成本低有的解运维成本高直接取其中一个值就会导致epsilon区间偏窄。解决方法是做字典序优化在求碳排放成本最小之后把碳排放成本固定在最优值再最小化或最大化运维成本得到准确的区间端点。这样得到的epsilon最大值才是合理的。这个细节论文里通常不会强调但实际复现时非常关键。5.2 坑2漏了CO2流约束出现“凭空消碳”的假解这个问题是我在验证结果时发现的。一开始我在模型中只写了碳捕集设备的捕集量约束没有把P2G消耗的CO2和碳捕集设备捕集的CO2关联起来。结果优化器给出了一个非常离谱的结果碳排放成本几乎是0P2G大量产气但碳捕集设备却没怎么运行。后来一查问题出在P2G甲烷化消耗的CO2没有来源约束优化器让它凭空消耗CO2来抵消碳排放相当于给系统开了一个无限供碳的作弊器。加上CO2_p2g CO2_cap这条约束之后结果立刻正常了。这个坑提醒我所有物料平衡类的约束尤其是碳流一定不能只看目标和成本的表达要回到能流图上逐条核对“物质的来源和去向”。5.3 坑3单位不一致导致结果偏离原文复现中另一个折腾人的问题是单位。天然气可以用m³、kg、kWh三种单位CO2排放因子可以用kg CO2/kWh、kg CO2/m³、tCO2/MWh热值有高位低位之分。我在最开始建模时把燃气轮机的发电效率写成了百分比但把天然气的能量按m³来算结果所有功率平衡约束全部对不上目标函数数值差了好几个数量级。建议在建模前把所有数据统一到一套单位体系我的做法是全部用kWh作为能量单位CO2量用吨碳价用元/吨。这样虽然需要多做一些换算但至少不会出现数值量级上的诡异偏差。换算关系要写清楚比如1 m³天然气约等于10 kWh按热值粗略换算CO2排放因子要用对应的热值基础。5.4 坑4求解时间随epsilon网格数急剧上升epsilon约束法本质上是把双目标转成多次单目标求解每次都是MILP。第一次我图方便设置N100结果跑了将近一个小时才出全部结果而且很多点实际上重复了。后来我把N降到30并给Gurobi设置了相对对偶间隙约束每次求解时间控制在几秒以内整体时间降到了五分钟以内前沿形状没有明显变化。这里有个经验是先用粗网格跑一遍找到前沿变化最剧烈的区域再在那个区域附近加密。这样既省时间又能保证前沿的细节不丢。另外每次循环求解时可以使用上一次的解作为初始热启动也可以显著加快收敛。我把这4个坑整理成一个表格方便对照问题现象根因解决方式前沿断裂epsilon区间端点求解不准用字典序优化确定端点碳排放被异常压低P2G消耗CO2无来源约束加CO2_p2g CO2_cap目标函数数值量级错乱单位不统一统一为kWh、吨、元求解时间爆炸网格过密、MILP时间未限制N30粗扫局部加密设MIP gap6. 帕累托前沿复现结果与一些验证建议代码跑通只是第一步真正重要的是结果对不对、能不能和论文吻合。我复现出来的帕累托前沿和预期比较一致下面简单说说结果以及我是怎么验证的。6.1 复现出的帕累托前沿长什么样以横轴为运维成本纵轴为碳排放成本画出来的帕累托前沿是一条向右下倾斜的曲线运维成本越低的一端碳排放成本越高碳排放成本越低的一端运维成本越高。前沿整体呈一个比较明显的L形或者凹形这个形状说明一个问题当系统在低排放区间运行的时候每降低一单位碳排放需要付出的运维成本增量越来越大边际减排成本是递增的。如果系统里P2G和碳捕集设备的耦合约束写对了前沿中段通常会出现一个“拐弯”区域。在这个区域里系统开始启动碳捕集设备并把捕集到的CO2供给P2G同时增加P2G出力来替代部分天然气外购。过了这个拐弯点进一步减排就需要大幅提高设备出力或者抑制电负荷成本上升会非常陡峭。6.2 不同运行区间下设备行为的变化我还对比了帕累托前沿两个极端点上的设备出力情况。运维成本最优点上燃气轮机基本处于满发或接近满发状态碳捕集设备可能只维持最低运行或不投运P2G的出力也较少系统主要通过外购电力和天然气来满足负荷。碳排放成本最优点上情况完全反过来碳捕集设备满负荷运行P2G尽可能多地利用系统内的电力和CO2来产气燃气锅炉尽量少用必要时甚至通过降低燃气轮机出力来减少碳排放。这两个极端解之间的折中解可以通过模糊隶属度函数从帕累托前沿上选择。我一般的做法是分别计算每个解对两个目标的隶属度然后选隶属度最大的点作为折中解。这个点对应的系统运行方式往往就是碳捕集和P2G都维持在中等出力水平燃气轮机和燃气锅炉合理分担热负荷整体效果很均衡。6.3 验证复现结果是否可靠的几个方法验证复现结果是否可靠我总结了几个比较实用的方法供你参考。第一单点核验。随机取某一个epsilon下的解检查电、热、气功率平衡是否严格满足检查碳捕集量是否不超过可捕集量检查P2G消耗CO2是否来自捕集量。任何一个约束不满足说明代码里有bug。第二和论文数值对比。如果论文公开了典型日负荷曲线和参数把数据输入后对比上下界端点数值。如果数量级对不上优先检查单位换算和效率参数的选取。第三灵敏度分析。改变碳价看帕累托前沿怎么移动。碳价越高最优解越偏向低碳端这是一种常识性的检验。如果碳价升高反而选择了高碳解那一定是模型或代码里有逻辑错误。第四MIP gap控制。所有MILP求解都有数值误差设置合适的相对MIP gap比如1%以内保证每次求解结果稳定。否则帕累托前沿上会出现一些无规律的噪声点很难看。复现一区论文最大的收获不是跑通代码那一刻的快感而是搞清楚哪些约束决定了结果的合理性。碳流约束、epsilon区间的边界计算、单位统一这三样东西任何一个不仔细都会让结果看起来“很合理”但实际上离物理真相很远。这篇Energy文章的模型确实有它的精巧之处把P2G和碳捕集设备通过CO2流串成了一个闭环再用epsilon约束法把双目标问题切得清清楚楚。如果手头有同学也在做类似的综合能源系统优化建议可以拿这个模型做起点后面再加储能、需求响应或新能源出力不确定性扩展空间很大。