最近在复现一篇Energy上一区的综合能源系统运行优化文章，模型里考虑了P2G（电转气）和碳捕集设备，配合热电联供机组，目标函数是碳排放成本加运维成本。原文章用的是加权法，我试着把求解器换成epsilon约束算法，直接跑双目标帕累托前沿。整套代码用Matlab实现，从中踩了不少坑，也理清了很多细节。这篇内容就把整个复现过程和扩展思路完整写出来，给同样做IES优化调度、多目标优化或者需要复现SCI论文代码的朋友一个能直接照做的参考。
1. 项目概述：复现任务到底在做什么
1.1 原文章核心与我的扩展点
先说清楚原文章的基本盘。它做的是含P2G和碳捕集设备的热电联供综合能源系统（IES）运行优化，核心是用一个典型日（或者几个典型场景）的时序数据，决定各个设备的出力、储能充放、P2G产氢/产气、碳捕集吸收CO2的量，使得整个系统在满足电、热、气负荷需求的前提下运行成本最低。文章里把碳排放成本也纳入目标，所以本质是一个“经济性+环保性”的双目标问题。
我复现的时候没有完全照抄原文的目标函数处理方式，而是增加了一个epsilon算法（也就是epsilon约束法）来求解双目标优化。为什么这么改？因为原文如果是单目标加权或者约束转化，那得到的只是一个“可行解”或者若干个权重下的折中解，很难真正看到帕累托前沿的形态。而epsilon约束法通过把一个目标转成约束，然后逐步收紧这个约束，可以系统性地扫描出完整的帕累托前沿，这是它相比线性加权法最实用的地方。
复现过程中，我把整个系统拆成了几个部件模型：热电联产机组（CHP）、燃气锅炉、电转气装置（电解槽+甲烷化）、碳捕集设备（吸收塔+再生塔抽象模型）、储能装置（电储能和热储能）、以及电网和气网的交互。这些部件的能量流关系是建模的关键，P2G消耗电能产生氢气和氧气，氢气可以和CO2反应生成甲烷，而CO2恰好来自碳捕集设备，这样就把P2G和碳捕集在物质流上耦合起来。这个耦合关系处理不好，求解器很容易报无解。
1.2 为什么选epsilon约束法处理双目标
很多刚接触优化的人会习惯性地用加权法，把碳排放成本和运维成本乘上权重加起来，变成一个单目标。这样做不是不行，但有一个很现实的问题：两个目标的数量级很可能不一样，比如运维成本可能是几十万，碳排放成本可能是几万，权重如果拍脑袋定，最终解会严重偏向某个目标。就算你用归一化处理，权重不同得到的不同解，也无法保证能覆盖整个帕累托前沿，尤其是目标空间非凸的时候，加权法会漏掉中间的帕累托点。
epsilon约束法的思路完全不同。它先把其中一个目标（比如碳排放成本）设为一个约束上限epsilon，然后优化另一个目标（比如运维成本）。通过不断改变epsilon的取值，比如从碳排放成本的最大值逐步缩小到最小值，每跑一次就得到一个帕累托最优解，最后把所有解连起来就是完整的帕累托前沿。这种方法对目标空间的凸性没有严格要求，能比加权法更完整地找回非凸部分的帕累托点。
用Matlab实现时，epsilon约束法的代码结构也不复杂。核心就是写一个循环，每次修改约束条件中epsilon的值，然后调用求解器（我用的Yalmip+Gurobi，或者直接用LinProg），把结果记录下来。我后面会详细讲这个循环怎么写，以及怎么避免循环过程中出现的数值问题。
2. 系统建模：P2G、碳捕集与热电联供的数学描述
2.1 能源hub结构梳理
做IES优化之前，第一步一定要画清能源hub的拓扑图。我们处理的是电、热、气三种能源品类的耦合系统，输入侧有电网购电、天然气网购气、可再生能源（比如风电/光伏）出力，中间经过CHP机组、燃气锅炉、P2G设备、碳捕集设备、储能设备，最后输出电负荷、热负荷和气负荷。
我把模型简化成一个节点式结构：电网节点、热网节点、气网节点、CO2流节点。电负荷由CHP发电、可再生能源、电网购电以及电储能放电共同满足；热负荷由CHP余热回收、燃气锅炉产热、热储能放电共同满足；气负荷由天然气网购气、P2G产生的天然气共同满足。CO2流则是碳捕集设备从CHP烟气中捕获CO2，一部分供给P2G甲烷化反应，另一部分可以封存或者出售。
这里要注意：P2G的甲烷化需要H2和CO2，H2来自电解水，CO2来自碳捕集，所以如果没有碳捕集，P2G还得额外购买CO2，经济性就差很多。原文章大概率也是看重了“碳捕集+P2G”这条闭环的价值。我建模时把CO2流也作为变量，而不是固定比值，这样优化器可以自己权衡“捕集多少CO2、用多少CO2产气、封存多少”。
2.2 设备建模与能量流方程
设备建模我用的都是线性化模型，因为要保证整个优化问题是MILP或LP，这样求解快、收敛稳。先列几个关键设备的简化方程。
CHP机组：输入天然气，输出电和热。我采用定热电比的抽凝式模型，功率关系为：
[ P_{chp}(t) = \eta_{chp}^e \cdot F_{chp}(t) ] [ H_{chp}(t) = \eta_{chp}^h \cdot F_{chp}(t) ]
其中F_chp是天然气消耗量（kW），eta_e和eta_h分别是发电效率和热回收效率。为了更灵活，可以引入热电比可调区间，但那样会引入二进制变量，复杂度上升。我复现时先用的定热电比，后续如果想扩展再改成可行域约束。
燃气锅炉就简单了，输入天然气产热：
[ H_{gb}(t) = \eta_{gb} \cdot F_{gb}(t) ]
P2G设备包含电解槽和甲烷化反应器。电解槽消耗电，产出氢气和氧气；甲烷化反应中，H2和CO2按2:1（实际化学计量1:4，但简化用2:1？需要准确：CO2 + 4H2 -> CH4 + 2H2O，所以H2与CO2摩尔比4:1。但为了简化，很多文章直接用P2G效率，认为P2G输入电能和CO2，输出天然气。我采用的公式是：
[ F_{p2g}(t) = \eta_{p2g} \cdot P_{p2g}(t) ]
同时CO2消耗量与产气量成正比：
[ C_{p2g}^{co2}(t) = \alpha_{co2} \cdot F_{p2g}(t) ]
alpha_co2是单位天然气产出消耗的CO2质量系数，由反应计量和分子量算出来。这里如果没搞清楚，后面约束容易写错。
碳捕集设备捕集CHP烟气中的CO2：
[ C_{ccs}(t) = \beta_{ccs} \cdot \gamma_{ccs} \cdot P_{chp}(t) ]
其中beta_ccs是烟气中CO2浓度相关系数，gamma_ccs是捕集率。捕集率可以作为0-1变量或者连续变量，若作为连续变量，则模型可动态调整捕集力度，对应不同的能耗成本。碳捕集本身需要消耗能量（热或者电），我简化为消耗一部分电，增加一个附加电负荷：
[ P_{ccs}(t) = k_{ccs} \cdot C_{ccs}(t) ]
k_ccs是捕集单位CO2需要的电耗。
储能模型用通用形式，电储能和热储能都一样：
[ E(t+1) = E(t) \cdot (1 - \sigma) + P_{ch}^{eff}(t) \cdot \eta_{ch} - \frac{P_{dis}(t)}{\eta_{dis}} ]
注意充放不能同时进行，需要引入二进制变量，或者用互补约束线性化。我用的Yalmip可以直接用binvar。
能量平衡约束是核心，每个能源节点都要满足供需平衡。电功率平衡：
[ P_{chp}(t) + P_{ren}(t) + P_{grid}(t) + P_{dis}^{e}(t) = P_{load}(t) + P_{p2g}(t) + P_{ccs}(t) + P_{ch}^{e}(t) ]
热功率平衡：
[ H_{chp}(t) + H_{gb}(t) + H_{dis}(t) = H_{load}(t) + H_{ch}(t) ]
气功率平衡：
[ F_{grid}(t) + F_{p2g}(t) = F_{chp}(t) + F_{gb}(t) + F_{load}(t) ]
CO2平衡：
[ C_{ccs}(t) = C_{p2g}^{co2}(t) + C_{store}(t) ]
这些方程写出来之后，整个系统的基本骨架就有了。
2.3 碳排放成本与运维成本的目标函数
原文章的目标函数通常是总成本最小，包含购能成本、运维成本、碳排放成本。我这次要拿出来做双目标的两个成本分别是：
运维成本：所有设备的运行维护费用，通常可以用单位出力的运维系数乘以出力来算。
[ C_{om} = \sum_{t} \left( \sum_{i \in devices} c_{om,i} \cdot P_i(t) \right) ]
碳排放成本：系统从电网购电隐含的碳排放和天然气燃烧直接排放的CO2，减去碳捕集量，再乘以碳价。表达式：
[ C_{co2} = c_{co2} \cdot \sum_{t} \left( \phi_{grid} \cdot P_{grid}(t) + \phi_{gas} \cdot (F_{chp}(t) + F_{gb}(t)) - C_{ccs}(t) \right) ]
这里的phi是碳排放因子，C_ccs是捕集量，被减去。如果碳捕集封存还得给封存成本，可以再加上。为了让epsilon约束法好用，我会把这两个成本都写成线性表达式。
注意：运维成本里面要包含P2G、碳捕集这些设备的启停或者变载成本吗？如果原文章不考虑启停，只考虑单位出力的线性成本，那就简单了。我做的时候加了几个二元变量表示设备启停状态，但后来发现求解时间变长不少。如果没有严格要求，建议先不做启停，线性出力的优化已经足够说明问题。
3. 双目标优化与epsilon算法的实现思路
3.1 多目标优化的两种常见套路
对双目标优化，最常用的是加权法和epsilon约束法。加权法的数学形式是min w1C1 + w2C2，通过变化权重得到一系列解。缺陷很明显：如果帕累托前沿非凸，加权法能得到的解不会落在非凸部分的凹陷处。而且权重如何选取没有统一标准，很多人就取0.5和0.5，一跑完发现解是在两个极端的中间，但对决策者来说，可能更想知道“如果碳排放限制在某个水平，最低运维成本是多少”。
epsilon约束法的形式是：
[ \min C_{om} ] [ s.t. C_{co2} \le \epsilon ]
然后让epsilon在C_co2的可能范围内变化。这样每一个epsilon对应一个优化问题，求解一次得到一个帕累托点。关键是epsilon的取值集要覆盖从“最环保”到“最经济”的范围。通常我先求两个极端解：不考虑碳排放约束时C_co2的最大值，以及强制碳排放为零（或极小）时C_co2的最小值。然后在这个区间里等分N个点，依次求解。
除了这两种方法，还可以用多目标进化算法如NSGA-II，但那是启发式算法，求出的解不保证精确帕累托最优，而且Matlab实现起来要自己写非支配排序，代码量大。对于这个规模不大的MILP模型，epsilon约束法配合精确求解器是更稳妥的选择。
3.2 epsilon约束法的原理与实现细节
具体的迭代流程分三步。
第一步，求解两个端点。先不加碳排放约束，只优化运维成本，得到解1，记录此时的碳排放成本E_max。再强制碳排放成本E_min（可以设为一个很小的值或者通过约束Cco2 <= E_min），优化运维成本，但这个可能无解，需要逐步放松下面的E_min，直到有解并得到最小的可行碳排放成本E_min。
第二步，在[E_min, E_max]之间取K个等分点，生成epsilon序列：epsilon_k = E_min + (E_max - E_min) * k / K，k=0,1,...,K。K通常取20~50，视计算时间而定。
第三步，逐个求解带Cco2 <= epsilon_k约束的运维成本最小化问题，得到运维成本和对应的碳排放成本，画图就是帕累托前沿。
这里有个细节：epsilon约束法在约束右端项变化时，需要保证每个子问题都有可行解，而且为了避免重复解，可以在增加epsilon约束的同时加上一个小的目标修正项，比如在目标函数中减去一个非常小的倍数的Cco2（比如-0.001*Cco2），来优先选择碳排放更低的解，从而让帕累托前沿上的点分布更均匀。这个技巧在文献里叫“epsilon约束法的增强版”，我也在代码里用了。
另一个细节是数值尺度。epsilon的值是碳排放成本，可能是几千到几万，而运维成本可能是几十万，如果直接混在一起，约束条件中的阈值不敏感时，求解器会报数值问题。我会在建模时把两个成本都除以一个基准值，比如总负荷均值，变成无单位的相对量，或者统一单位到“元”后，只在约束里使用epsilon，并且让epsilon的步长选择足够大，避免因为浮点误差导致相邻两个epsilon解相同。
3.3 如何生成帕累托前沿
生成帕累托前沿不是只把点画出来就完了，还需要做后处理。我一般会把每个优化子问题得到的C_om和C_co2存成一个两列矩阵，然后剔除重复解和劣解。重复解是因为epsilon间隔太小，或者目标函数修正项没起效，导致连续两个epsilon产生相同的解。劣解则是因为某些epsilon下求解器数值误差产生了被支配的点。剔除后，把剩下的点按C_co2从小到大排列，就是最终的帕累托前沿。
如果还要对比原文章的结果，可以把原文章加权法得到的解也画在一个坐标里，看看是否落在帕累托前沿上。我复现时发现加权法等权重解确实落在我求出的前沿上，但非凸区域的解在加权法下是缺失的，这就证明了epsilon约束法的优势。
4. Matlab代码实现：从框架到细节
4.1 数据准备与参数初始化
整个代码我用的是Yalmip建模，求解器选Gurobi。Yalmip的好处是可以用符号化方式写约束，变量多了也不会乱。如果你没有Gurobi，用Matlab自带的linprog或者intlinprog也可以，但求解速度会慢不少，而且处理二进制变量时intlinprog有时候会不稳定。我强烈建议装一个Gurobi或者Cplex，学生版免费。
参数初始化我是用结构体存储的，比如para.chp_e = 0.35，para.chp_h = 0.45，para.eta_p2g = 0.6，para.carbon_price = 50。这样读起来清楚，后期调参也方便。数据部分，我加载了一个24小时的负荷曲线、风电出力和电价、气价。这些数据可以自己造，也可以用现成的测试系统数据。原文章如果是用某地区实际数据，复现时没有原始数据，我只能用类似量级的合成数据，然后说明结果趋势一致。
在写Yalmip模型之前，先把时间长度T设为24，所有变量定义成T维向量，或者定义成矩阵。比如Pchp = sdpvar(1,T)，表示每个时段的CHP电出力。
4.2 约束条件与变量的矩阵化处理
我习惯先把所有变量定义在一个struct里，方便引用。例如：
T = 24; x.Pchp = sdpvar(1, T); x.Pgb = sdpvar(1, T); % 燃气锅炉耗气 x.Hgb = sdpvar(1, T); x.Pp2g = sdpvar(1, T); x.Pccs_elec = sdpvar(1, T); x.Cccs = sdpvar(1, T); x.Pgrid = sdpvar(1, T); x.Fgas = sdpvar(1, T); % 购气 x.SODe = sdpvar(1, T+1); % 电储能容量 % ... 热储能类似约束用for循环或者矩阵形式。Yalmip支持直接在约束里用循环，但矩阵化更紧凑。比如电平衡：
C = []; C = [C, x.Pchp + P_re - x.Pp2g - x.Pccs_elec + x.Pgrid == P_load];如果不等式约束要加上下限：
C = [C, 0 <= x.Pchp <= para.Pchp_max];对于储能的充放互斥约束，我引入两个二进制变量，然后加约束：
on_ch = binvar(1,T); on_dis = binvar(1,T); C = [C, on_ch + on_dis <= 1]; C = [C, Pch_storage >= 0, Pch_storage <= on_ch * Pch_max]; C = [C, Pdis_storage >= 0, Pdis_storage <= on_dis * Pdis_max];这样能保证充放不同时进行。求解时加上binvar变量会让模型变成MILP，但规模不大，Gurobi处理起来很快。
P2G和碳捕集的耦合约束比较关键，我在后面专门讲。
4.3 调用求解器与循环epsilon求解
我先写一个函数run_ies(para, data, epsilon, is_feasibility)，输入参数包括epsilon（碳排放上限），输出优化的成本和变量。函数内部用Yalmip建模，设置目标为运维成本（如果epsilon是-1，表示不加碳排放约束，求两个极端点）。
核心调用代码大致是这样的：
% 第一次优化：不加碳约束，求最大碳排放 res_max = run_ies(para, data, -1); E_max = res_max.Cco2; % 第二次优化：设一个很小的碳上限，如果不可行，逐步放松 eps_lo = 0; while true try res_min = run_ies(para, data, eps_lo); if res_min.problem == 0 break; end end eps_lo = eps_lo + delta; end E_min = res_min.Cco2; % 生成epsilon序列 K = 30; eps_list = linspace(E_min, E_max, K); % 循环求解 pareto = []; for k = 1:K res = run_ies(para, data, eps_list(k)); if res.problem == 0 pareto = [pareto; res.Cco2, res.Com]; end end在run_ies内部，目标函数我设置成：
obj = res.Com - 1e-4 * res.Cco2;因为当epsilon约束Cco2 <= eps_k时，目标本身是min Com，但加上一个极小的负Cco2项后，求解器会在满足Com最小的前提下尽量让Cco2更小，这样有利于得到均匀的前沿点。这个系数不能太大，否则就变成优先环保了；我实测1e-4在数值尺度下很安全。
注意：由于每次循环都会新建Yalmip模型并求解，时间会随着K线性增长。K取30时，每个子问题大概0.2~0.5秒，总共也就十几秒，完全可以接受。如果你的模型加入了启停等太多整数变量，可能一个子问题就要几十秒，这时K取小一点，比如15，或者改用热启动，把上一个解作为初始点，加速收敛。
4.4 结果可视化与帕累托前沿绘图
绘图是Matlab的强项。我先画出帕累托前沿散点图，再画折线，加上坐标轴标签。通常横坐标是碳排放成本，纵坐标是运维成本，单位都是“万元”或者“元”。
figure; plot(pareto(:,1)/1e4, pareto(:,2)/1e4, '-o', 'LineWidth', 1.5); xlabel('碳排放成本（万元）'); ylabel('运维成本（万元）'); title('epsilon约束法得到的帕累托前沿'); grid on;如果想更直观，可以把每个帕累托点对应的P2G出力和碳捕集量也画出来，看看随着碳排放限制变严，P2G和CCS的出力是如何变化的。我做的实验显示，当碳排放上限比较宽松时，系统更多依赖电网购电，设备出力低；当碳排放上限收紧，系统不得不增加P2G消耗风电、开大碳捕集，运维成本上升。这个趋势和物理直觉一致，也是对模型正确性的一个验证。
我还会输出两个极端方案的具体调度结果，比如调用stairs画电功率平衡图，看CHP、风电、电网、P2G、碳捕集在各个时段的平衡情况。
5. 常见问题与排查经验
5.1 求解无可行解，先查耦合变量一致性
这个问题我踩了好几次。症状是加入P2G和碳捕集的CO2耦合后，求解器提示infeasible。排查方法：先把CO2流约束关闭，单独优化P2G而不考虑CO2来源，看是否有可行解，再单独优化CCS，看是否有可行解。如果单独都可行，加耦合就不可行，那多半是CO2平衡约束写反了，或者产量系数和消耗系数对不上。
比如甲烷化反应需要的CO2量应该等于P2G产气量乘以一个系数，但很多人会把系数设为0.2或者0.5，实际上一立方甲烷需要多少CO2，要从分子量和标况体积去换算。我用的是简化系数，但需要保证在P2G最大出力时，CCS的最大捕集量能够满足CO2需求，否则会在某个时段出现“用气比产气还快”的无解情况。
另外，P2G和CCS同时消耗电能，如果电网购电上限加上可再生能源出力，不足以同时供负荷、P2G、CCS，也会无解。要检查电平衡约束中Pccs_elec和Pp2g是不是都被加到了负荷侧。
5.2 目标函数标幺化与数值尺度问题
Yalmip+Gurobi对数值范围比较敏感。如果你的购电成本是几百元，碳排放成本是几十万元，两者放在一个目标函数里，Gurobi默认的容差下可能忽略小项，导致解不精确。我建议对成本和功率都做标幺化处理，比如选取一个基准功率base=1000kW，所有功率变量除以base，成本除以base或除以一个基准价格。这样数量级都在0.1~几十之间，求解稳定很多。
我做epsilon约束时，直接对原始成本求epsilon区间，步长之间的差距可能只有几十元，但Gurobi的整数可行域解可能对这么小的步长不敏感，导致重复解。后来我把epsilon序列取成对数间隔或者先归一化到[0,1]再到原始值，效果更好。
5.3 epsilon序列生成的两个误区
误区一：epsilon最大值用“不约束时碳排放成本”，最小值用“0”。如果最小值设置为0，很可能系统根本做不到零碳排放，于是前面一大半epsilon无解，白白浪费计算时间。正确做法是动态求最小碳排放，或者从E_max往下搜索，找到第一个可行解作为E_min。
误区二：epsilon间隔均匀划分就一定好。实际上帕累托前沿在某些区域变化剧烈，均匀划分会导致剧烈区域点很稀疏，平坦区域点很密。如果想更有针对性，可以先用粗网格跑一遍，找到拐点区域，再在拐点区域细化epsilon步长。我在复现时用了两阶段：第一阶段K=20均匀，第二阶段在相邻点变化斜率最大的区间内部再插入点，这样画出来的前沿更顺滑。
5.4 一个容易忽略的坑：P2G与碳捕集的时序耦合
P2G和碳捕集不仅有同一个时段的耦合，还有跨时段耦合吗？如果碳捕集捕集的CO2储存在储气罐里，供后续时段使用，那就要加一个CO2存储的动态约束。很多简化模型不做储存，假设即时用完，这样每个时段的CO2捕集量必须等于CO2消耗量加封存量。但实际中P2G可能夜间风电多的时候多出力，而白天CHP满发时碳捕集才多，两者峰值不同步，即时耦合会限制系统灵活性。
我建议在建模时给碳捕集加一个CO2储罐，变量是储罐容量，约束类似储能。这样P2G不一定要和CCS同时运行，产甲烷时可以从储罐取CO2。这个扩展在原文章里可能没有，但加上之后系统运行成本会下降不少，而且更符合实际。如果你想跟原文章结果严格对比，就先不加；如果想做改进分析，加上很有意思。
5.5 求解速度慢的优化技巧
如果模型较大，比如T=24且有多台CHP、P2G、储能，整数变量增多后，Gurobi要反复分支定界。我这篇文章因为是复现，设备数量不大，所以还好。但如果以后扩展成多场景随机优化，需要考虑以下技巧：一是用Yalmip的assign和initial为每个子问题设置初值；二是把epsilon约束法的子问题改成增量式，上一个解作为热启动；三是尽量消除多余的二进制变量，比如储能用分段线性化代替二进制；四是使用big-M法时M不要给太大，够用就行，否则数值条件数很差。
6. 我的实操心得与后续扩展
复现这个题目，最大的收获不是套用公式，而是把“P2G+碳捕集”这个耦合系统的能量与物质流真正想明白了。最开始我照着原文的图直接写代码，结果总是对不上，后来自己画了一张完整的能量和物质流图，把每一个输入输出和单位都标出来，再把方程一条条写出来，问题瞬间清楚。建议看到这里的朋友，不要急着写代码，先动手画图，这会节省你一天时间。
epsilon约束法的代码实现比我预想的简单，但用起来有几个细节值得记住：一是一定要先求两个端点确定epsilon范围，二是目标中加极小正则项避免重复解，三是注意数值标幺化。做到这三点，帕累托前沿基本一次跑通。
后续如果要扩展，可以往这几个方向走：考虑不确定性，用鲁棒优化或者场景随机规划替代确定性模型；加入需求侧响应，让负荷也变成可调变量；或者把碳排放成本改成阶梯碳价，更贴近实际碳交易市场的机制。我自己已经在尝试把模型改成两阶段鲁棒优化，第一阶段决定设备启停，第二阶段在不确定性场景下做经济调度，计算量明显增大，但结果也更有说服力。
最后分享一个小技巧：当你调试一个复杂的IES优化模型时，从24小时缩成1小时先跑通。把T改成1，负荷取一个典型值，所有时序约束变成单点约束，这样如果模型本身有什么低级错误，立刻就能暴露。我每次写新模型都用这个方法，几乎没被大模型的“神秘bug”耗到半夜。
代码的骨架我已经放在文中的章节里了，具体数据因版权原因我就不贴完整原始数据了，你可以根据自己研究的区域设定负荷曲线和能源价格。跑通之后，再逐步把T恢复到24，验证时序储能的动态约束。这个过程虽然有点繁琐，但每一步都是扎实的。
