1. 复现背景与核心价值拆解
先说结论：这篇论文的复现难度在同类综合能源系统优化文章里算中上，但它值得做。为什么？因为它把两个容易被人忽略的约束同时摆上了桌——碳排放成本和运维成本，而且用的是求解双目标问题的epsilon约束算法，不是那种简单加权求和就完事的做法。
我一开始看到标题里“P2G”和“碳捕集”这两个词时，第一反应是“又来一个堆设备的”。但真正把代码跑起来、把结果图复现出来之后，才发现这篇文章的核心价值不在设备堆砌，而在“如何把碳排放成本这个非线性、难处理的量，用一个能工程落地的算法写进优化模型里”。很多文章喜欢把碳捕集写得很花哨，但一做仿真就露馅——不是模型不收敛，就是求解时间长得离谱。这篇文章的epsilon算法把双目标问题转成单目标参数化求解，思路干净、代码量适中，非常适合用来练手多目标优化的落地实现。
适合谁参考？三种人：一是正在做综合能源系统、微电网、热电机组优化方向的硕士博士，尤其是需要复现SCI文章做对比实验的；二是对epsilon约束算法感兴趣、想看看除了NSGA-II之外还有什么实用多目标解法的研究者；三是企业里做园区能源调度、想算清楚“减碳到底要花多少钱”的工程师。这篇博文我会从模型结构、约束处理、目标函数线性化、epsilon算法实现、Matlab代码结构几个角度展开，最后附上我踩过的坑。
需要提前说明的是，本文中涉及的具体参数、公式推导细节，我会结合“综合能源系统优化”这一领域的常见实践进行补充说明。原文未明确给出的部分，我会标注“基于常见实践补充”，供你参考。
2. 系统架构与数学模型拆解
2.1 热电联供系统的能量流结构
先梳理一下这个系统在干什么。所谓“考虑P2G和碳捕集设备的热电联供综合能源系统”，本质上是一个典型的多能源互补园区：电网、气网作为外部输入，内部有燃气轮机或热电联产机组（CHP）、P2G（电转气）设备、碳捕集装置、储能设备等，输出端是电负荷和热负荷。
我画过很多次这个系统的能量流图，最核心的一条主线是：CHP机组燃烧天然气发电产热，发电的同时产生大量二氧化碳，这部分二氧化碳不直接排放，而是送进碳捕集装置，捕集下来的CO2再和P2G设备产的氢气（来自电解水）合成甲烷，重新送回气网或者给CHP当燃料。这就是“碳循环”的闭环思路。
用白话翻译一下就是：以前燃煤燃气机组是“发电-排碳-完事”，现在是“发电-捕碳-制气-再烧”，把一部分碳截留在系统内部循环利用，外购天然气和实际碳排放量都降下来。代价是什么？设备投资和运行成本上去了，P2G电解槽耗电，碳捕集装置耗热耗电，而且这几台设备之间互相耦合，调度难度明显加大。
这条能量流结构里有个细节容易被忽略：P2G和碳捕集设备不是独立的“减排装置”，它们和CHP、储能之间是强耦合关系。P2G要用电，电从哪来？如果电网电价高，P2G运行不划算；如果让CHP多发电供给P2G，CHP的碳排放又增加，碳捕集又得多开。所以这本质上是一个“牵一发动全身”的协调优化问题。
2.2 决策变量与目标函数解析
这个优化问题的决策变量可以归成三类，我建议你在复现的时候也按这个分类去组织代码：
- 第一类是设备出力变量：CHP的电出力、热出力、燃气轮机消耗的天然气量、P2G的产氢量/产甲烷量、碳捕集装置的捕集量、各类储能的充放电/充放热功率。
- 第二类是交互变量：从电网购电/售电功率，从气网购气量。
- 第三类是衍生变量：捕集下来的CO2去向、甲烷合成量、系统总碳排放量等。
目标函数是本文的关键。题目里明确说了是“碳排放成本+运维成本的双目标优化”，不是单一经济成本。这里有个概念要先掰扯清楚：很多文章写“双目标”，实际上是“经济成本最小化”和“碳排放量最小化”，但一边是钱一边是吨，单位不统一，最后往往就变成线性加权。而这篇文章的处理方式是——把碳排放量乘以一个碳价（单位碳排放成本），折算成碳排放成本，然后和运维成本放在同一个数量级上做双目标。
这个处理方式的好处非常明显：两个目标都是“钱”，物理意义清晰，决策者可以直接看到“减排目标值对应多大的成本代价”。而且用epsilon约束法求解时，可以自然地扫描碳排放成本的上限值（即epsilon参数），得到完整的Pareto前沿。我实测下来的体会是，这种“同量纲双目标”的Pareto前沿图非常直观，审稿人和项目汇报时都更容易理解。
运维成本部分，通常包括各设备的运行维护费用，常用单位出力运维系数乘以出力量来近似。设备之间的启停成本、碳排放惩罚/交易成本也可以根据模型需要纳入，属于常见实践补充范围。
3. 碳捕集与P2G的建模要点
3.1 碳捕集设备的运行约束
碳捕集设备建模，最大的坑在于它的“能耗”和“捕集量”不是简单的比例关系。常见的工程近似是：捕集单位质量的CO2需要消耗固定比例的电能和热能，同时捕集率有上下限约束。
我用一个简化的模型来展开说明，这个线性化方式也是我在复现时验证过最稳的：
- 捕集量约束：(0 \le C_{capt} \le \eta_{cap} \cdot C_{total})，其中(C_{total})是CHP机组理论总排放量，(\eta_{cap})是捕集率上限。
- 能耗约束：捕集系统消耗的电/热功率正比于捕集量，即(E_{cc} = k_{elec} \cdot C_{capt})，(H_{cc} = k_{heat} \cdot C_{capt})。
注意一个隐含问题：碳捕集的能耗加在电负荷和热负荷上，相当于在能量平衡方程里多挖了一块“内部消耗”。很多初学者复现这类模型时，最容易犯的错误就是忘记把这部分能耗写入电/热功率平衡方程，结果系统的总负荷被低估，优化结果偏乐观。
我在实际复现时建议把碳捕集“三件套”——捕集量、耗电量、耗热量——全部设为独立变量，它们之间的线性关系用等式约束表达。这样做的好处是：一方面便于后续扩展成非线性的捕集效率曲线，另一方面调试时也能更清晰地定位是哪个约束出了问题。
3.2 P2G设备的能量转换效率
P2G设备的核心是电解水制氢，然后氢气和捕集到的CO2反应生成甲烷。这个链路的能量转换效率是建模重点。
从能量流角度来看：1标方氢气对应约3.54 kWh的电能消耗（电解槽效率按70-80%算），氢气和CO2反应合成甲烷时又有约10%-15%的能量损失。所以P2G整体从“电”到“天然气”的能量效率通常在50%-65%之间。
这个效率直接决定了P2G在经济上是否划算。简单算一笔账：如果电价为0.6元/kWh，生产1 kWh当量的天然气需要消耗约1.7 kWh电能，光电力成本就是1.02元——市售天然气才3元/m³（约0.3元/kWh当量）。这么一看P2G经济性极差，那它为什么还有研究价值？因为它的价值不在经济性，而在“消纳弃风弃光”和“碳循环利用”。所以在双目标模型里，P2G能否被调度启用，完全取决于碳排放成本在目标函数里占多大权重——这就是双目标优化在工程上的意义所在。
我复现时对P2G采用的是线性效率模型：产甲烷量 = 电输入功率 × 综合转换效率。如果原文有更细的阶梯效率或者变工况效率曲线，需要做分段线性化，这个我放在后面epsilon算法部分展开。
4. 双目标优化与epsilon约束算法实现
4.1 为什么选epsilon约束而不是加权和
大多数入门教材教的多目标处理方式是线性加权求和——给每个目标一个权重，然后求单目标最优解。但用在这个问题上，你会遇到一个尴尬场景：权重怎么定？定0.5/0.5？人工拍脑袋定出来的权重往往只能得到Pareto前沿上的一个点，想画出完整前沿曲线就得反复试权重，而且线性加权在非凸Pareto前沿上会漏掉一部分解。
epsilon约束法的思路完全不同：保留一个目标作为主目标（比如运维成本最小化），把另一个目标（碳排放成本）放进约束条件限制它的上限——(C_{carbon} \le \epsilon)，然后通过连续扫描epsilon的取值，得到一系列单目标优化问题，每个问题求一个最优解，拼起来就是完整的Pareto前沿。
这一步是整个复现工作的灵魂。我不止一次看到有人复现文章时，把epsilon约束法理解成“加一个上界约束就行”，结果epsilon取值区间没设计好，要么约束从不生效，要么问题直接无解。下面我会给出一个可复现的epsilon扫描方案。
4.2 epsilon取值的设计与扫描技巧
epsilon取值的核心原则是：先放开约束求两个单目标极端解，再用极端解来确定epsilon的范围。
我在复现时按下面三步走，供你参考：
- 第一步：求“纯最小运维成本”问题，不约束碳排放，得到最优运维成本(f_{1}^{min})和对应的最大碳排放成本(f_{2}^{max})。
- 第二步：求“纯最小碳排放成本”问题，不约束运维成本，得到最优碳排放成本(f_{2}^{min})和对应的较大运维成本(f_{1}^{max})。这里注意：每次调目标函数时，另一个目标对应的值要记录下来。
- 第三步：在区间([f_{2}^{min}, f_{2}^{max}])里等间距取N个epsilon值（我一般取10-15个点），每个点求一次“运维成本最小化 + 碳排放成本 ≤ epsilon”的单目标问题，得到N个Pareto解。
epsilon数组建议用linspace生成，但要注意一点：epsilon=最小值时约束最紧，很可能无解；epsilon=最大值时约束最松，退化成单目标。所以实际绘制Pareto前沿时，首尾两个点需要单独判断可行解状态，不要直接丢进求解器不管。
关于求解器选型，YALMIP+GUROBI/CPLEX是常规做法，属于该领域通用工具，MATLAB环境下的求解器选型可按需补充说明。如果你不想依赖商业求解器，也可以用MATLAB自带的linprog，但大规模问题求解会很吃力，尤其是加了整数变量之后速度慢得让人崩溃。
4.3 非线性项的线性化处理
这个模型里有个非常容易卡住的地方：CHP机组的热电出力可行域通常是非凸的，碳排放量可能是出力的二次函数，P2G的效率也可能随负荷变化。这些东西一旦直接扔进求解器，轻则求解时间爆炸，重则模型不收敛。
我用的处理方式是“分段线性化 + 二进制变量”。以碳排放函数为例：把CHP电出力区间切成3-5段，每段用线性函数逼近碳排放量曲线。引入二进制变量表示“当前出力落在哪一段”，再用大M法保证只能激活一段。这样做会让模型从纯LP变成MILP，但epsilon算法本身要反复求解多次，MILP在中小规模系统下配合GUROBI其实完全能接受。
一个必须提醒你的细节：分段线性化和大M法处理不当很容易产生不可行解或者次优解。关键是把分段断点选合理（不要选太密，3-5段足够），同时大M的取值不能太大，过大的M会被坏数值稳定性。我习惯给M设一个比约束右侧最大值大一个数量级的数，同时做一期预求解来验证模型不是“伪可行”。
5. Matlab代码实现与运行结果解读
5.1 参数初始化的注意事项
参数初始化直接决定模型是否合理。我复现时最常被坑的是“单位混乱”——有的参数用kW，有的用MW，有的能量单位用kWh，有的用GJ，换算错了结果图看起来特别妖。
建议在代码顶部建立一个统一的参数结构体，注释里标明单位。下面是我推荐的组织方式，基于常见实践补充：
- 负荷数据：电负荷、热负荷，单位kW，时间分辨率取1小时。
- CHP参数：最大/最小电出力、热电比上限/下限、发电效率、耗气系数。
- P2G参数：电解槽额定功率、电转气综合效率、最大输入功率。
- 碳捕集参数：捕集率上限、单位捕集电耗、单位捕集热耗。
- 储能参数：容量上限、充放功率上限、充放效率、初始SOC。
- 价格参数：上网电价、购电电价、天然气价格、碳价。
特别注意：分时电价是一天24个点的序列，碳捕集和P2G的启停策略会强烈依赖电价曲线。如果电价是平的，P2G基本不会启动；只有电价低谷时段，P2G才会被调度。这也是模型合理性的一个直观检验标准——便宜的电拿来制气，贵的电拿去卖。
5.2 核心约束的Matlab实现示例
以电功率平衡为例，我在YALMIP中通常这样写：
C = [C, P_chp + P_dis + P_grid_buy + P_p2g_out == P_load + P_p2g_in + P_cc_elec + P_chg]; % P_dis: 储能放电功率 % P_grid_buy: 购电功率 % P_p2g_out: P2G输出的电（实际为0，P2G是耗电设备） % P_p2g_in: P2G消耗的电功率 % P_cc_elec: 碳捕集消耗电功率 % P_chg: 储能充电功率这里有一个容易看懵的地方：P2G那个“out”是什么意思？其实在能量流建模中，P2G是一个纯负荷，消耗电产生气，所以在电平衡里它只出现在等式左边（作为负荷项），不会出现在右边（作为电源项）。我上面代码里写P2G_out=0是提示作用，实际不需要这个变量。
热功率平衡同理：
C = [C, H_chp + H_dis + H_cc_heat == H_load + H_chg]; % H_cc_heat: 碳捕集消耗的热功率注意：碳捕集设备消耗热功率这个点在很多简化模型里会被省略。如果原文对捕集能耗做了简化，你要么跟随原文，要么在复现说明里指出该简化对结果的影响。完整模型的目标是贴近实际，简化模型的目标是复现原文，两者要在“与原文结果对比”的大前提下取舍。
碳排放约束的实现：
C = [C, C_total == C_chp + C_grid - C_capt]; C = [C, C_total <= epsilon]; % epsilon约束法的核心这里(C_{grid})指的是外购电对应的上游碳排放（间接排放）。如果不考虑这部分，就把它设成0。加了这部分之后模型更完整，但参数需要额外设置电网碳排放因子。
5.3 运行结果与Pareto前沿的可视化
求解完成后，每个epsilon值会得到一个解，包含设备出力序列、成本值、碳排放值。收集这些数据后，用scatter绘制Pareto前沿：
figure; scatter(carbon_obj, cost_obj, 60, 'filled'); xlabel('碳排放成本/元'); ylabel('运维成本/元'); title('Pareto前沿'); grid on;我复现时得到的典型Pareto前沿是一个下降曲线：碳排放成本低的方向，运维成本高；碳排放成本高的方向，运维成本低。曲线形状类似指数衰减，拐点通常出现在碳排放约束从“宽松”转“收紧”的区间。这个拐点的工程意义是——再往下压碳排放，运维成本会急剧上升，性价比很低，是决策者最关心的“甜点区间”。
除了Pareto前沿，我还会画几个关键场景下的设备出力堆叠图：碳排放约束最宽松时的调度方案、约束最严格时的调度方案、以及折中点的调度方案。三张图并排对比，能非常直观地看到“碳捕集和P2G是什么时候被启用的”。
我的实测结果是：碳价高/碳约束紧时，P2G在夜间低谷电价段启动频率明显增加，碳捕集设备的捕集率拉满；碳约束松时，P2G基本不启动，CHP直发直排，系统回到传统热电联供模式。这个行为完全符合工程直觉，也是验证模型正确性的一个重要指标。
6. 调试心得与常见问题速查
6.1 模型不可行的排查流程
epsilon约束法最常见的坑就是“约束设太紧导致模型无解”。如果你遇到YALMIP报infeasible，我的排查顺序是这样：
第一步，检查epsilon是不是设到了(f_{2}^{min})以下。如果epsilon比碳排放理论最小值还低，那就是人为制造了矛盾约束，必须把epsilon下界放松。
第二步，检查碳捕集捕集量上限是否足够低。如果CHP总排放量本来就不大，而碳捕集又被要求捕集到某个量，可能捕集上限无法满足约束。这种情况下要么调低碳排放目标下限，要么放宽捕集率上限。
第三步，检查电热功率平衡是否满足CHP可行域约束。CHP机组的热电比是有限制范围的，热出力大时电出力不能太小，反之亦然。当碳约束变紧时，调度器可能试图让CHP降出力，但热负荷还在高位，这时候就可能出现“电够了但热不够”的冲突——本质是CHP热电耦合约束和热负荷刚性约束打架。
第四步，如果还找不到原因，把epsilon直接设成inf，重新求解。如果设成inf仍然无解，那就是模型本身写错了，问题不在epsilon，而在于约束或变量定义有冲突。这个时候需要逐条约束注释排查，不要偷懒。
6.2 求解时间过长的优化技巧
epsilon算法要跑N次优化，如果每次都求解很慢，你会等到怀疑人生。几个可行技巧：
一是关闭不必要的输出。设置YALMIP和求解器的显示等级，别让求解器每个节点都往控制台刷日志。
二是把完全相同的公共约束提到循环外面，只用循环内的循环变量更新epsilon约束。我见过有人把整个模型放在for循环里重建，N次迭代就是N次模型构建，每次还得重新解析，浪费时间不说，还可能因为数值误差累积出问题。
三是如果用了二进制变量做分段线性化，可以尝试把分段数从5段减到3段，先验证模型逻辑是否正确，最后再用完整分段去精算。分段数对结果的影响通常很小，但对求解速度的影响是巨大的。
6.3 我对这套模型局限性的一点看法
说句实在话，这个模型和代码虽然能复现出漂亮的结果，但它离工程应用还有距离。最大的问题在于：碳捕集和P2G的变工况特性被简化成了线性模型，实际的催化剂活性衰减、电解槽启停寿命损耗、碳捕集溶剂的再生能耗变化，这些动态特性都没办法在这个框架里体现。
但从学术复现和算法学习角度来说，这个模型是极好的教学素材。你把epsilon约束法吃透了，之后看任何“多目标+复杂约束”的优化问题都会觉得清晰很多——因为工具是通用的，变的只是具体约束形式。
7. 写在最后的实操建议
按这篇文章的配置，完整复现一遍，从读论文到跑出Pareto前沿，我建议你预留两到三周时间。第一周搭模型和写代码，第二周调试和验证结果，第三周做对比实验和画图。不用焦虑，踩坑是正常的，重要的是把每一个坑的原因想清楚。
最后给你一个我实际用下来很顺手的检查方法：拿到求解结果后，先别急着看Pareto前沿，去检查最基础的能量平衡是否满足——把每个时段的电功率、热功率的供给项和需求项分别求和，相减，误差应该在1e-6量级。如果对不上，那一定是代码某处约束写错了。能量守恒是这套模型的试金石，任何时候都有效。
另外，原文如果提供了结果图，你在复现时不要苛求和原文的每个点都一模一样，因为参数设置略有差异、求解器版本不同都会影响数值结果。你要复现的是“趋势”和“形状”——Pareto前沿的走势、设备调度的规律、碳排放约束收紧时成本的变化斜率，这些对上了，复现就是成功的。
