1. 复现项目的本质搞懂这篇工作到底在优化什么做综合能源系统方向的人对“考虑P2G和碳捕集设备的热电联供综合能源系统运行优化”这个题目应该不陌生。这是一篇发表在Energy上的SCI一区论文的复现任务原文章本身的工作量集中在两个点上一是引入了P2G电转气与碳捕集设备的耦合二是在传统单目标运行优化基础上把它升级成了碳排放成本运维成本的双目标优化问题。复现时还要额外加上epsilon约束算法来求解这个双目标问题这也是整个复现项目最有意思、也最踩坑的地方。先说清楚这个系统是干什么的。热电联供CHP机组本身就存在“以热定电”或者“以电定热”的刚性问题风电、光伏大量接入后弃电问题突出而系统碳排放又主要来源于燃气轮机和燃气锅炉。P2G的作用是把富余电力转化为天然气或者氢气再回注到燃气设备里使用碳捕集设备则是把烟气里的CO2抓回来一部分提供给P2G作为合成天然气的原料一部分压缩封存或者外售。这样一来电、气、热、碳四条链路就耦合在了一起。你要做的运行优化就是在这套多能流耦合系统里决定每台设备每个时段的出力、启停、储气/储热充放策略让总成本尽量低。复现这类工作最容易犯的错误是“只求跑通代码不去抠物理模型”。因为论文里的公式符号多设备模型描述分散如果你不先把能量流程图和平衡关系理清楚后面改约束、调参数、加epsilon约束循环的时候完全不知道bug是出在求解器上还是出在模型本身。我建议复现的第一步不是写代码而是画一张完整的能量枢纽拓扑图把这套系统的“能流骨架”定下来。这也是本文要讲的第一个重点。原文章的核心工作我理解下来是两件事第一构建包含P2G与碳捕集设备的热电联供IES优化模型第二用双目标优化框架对比单目标框架下系统运行策略与碳排放量的差异。由于原文章是以传统加权法处理双目标复现要求在它的基础上改用epsilon约束法重新求解本质上是一次“方法升级型复现”。这就要求你不仅会读代码还得真正理解原模型的目标函数与决策变量关系否则epsilon约束法的两个目标函数谁当主目标、谁当约束都会选错。2. 优化建模的第一关把能量枢纽拓扑和设备模型吃透2.1 拓扑结构是地基先画图再建模先说说这套系统的拓扑。典型的架构是这样的外部电网、风电场、光伏电站作为电能输入端CHP机组燃气轮机余热锅炉承担主要的热电联供燃气锅炉作为热力备用电储能、热储能作为时间转移装置P2G设备和碳捕集设备是新增的两个核心环节。天然气网既向CHP和燃气锅炉供气也接收P2G产出的天然气或富氢天然气碳捕集设备从CHP烟气和外部碳源中捕集CO2精制后一部分送到P2G作为原料一部分进入储碳罐或外售。关键的一点是P2G和碳捕集设备的耦合关系。P2G制甲烷的反应是CO2 4H2 → CH4 2H2O也就是说P2G合成天然气需要消耗CO2和H2。H2来自电解水需要耗电CO2来自碳捕集装置。这个耦合意味着两个设备的运行状态必须联立建模如果P2G在满负荷产气碳捕集设备提供的CO2不够那P2G的产气量就要受限反过来如果碳捕集设备捕集了大量CO2而P2G没有同步开启捕集的CO2就只能在储碳罐里积累产生存储成本或者封存成本。很多复现代码跑出来的结果不合理大概率是这里解耦处理了。我复现时采用的拓扑里风电、光伏、外部购电、CHP发电都汇到电母线上电负荷、P2G用电、电储能充电是电母线的输出。CHP余热、燃气锅炉产热、热储能充放都汇到热母线上供应热负荷。气网供给CHP和燃气锅炉P2G产气可以进入气网缓存罐。这里要注意一个细节多数论文里P2G产气是“替代天然气购入”的所以模型里要有一个天然气购买项的削减逻辑而不是简单地把电和气的转换效率乘起来。2.2 设备模型的关键公式与参数选择CHP机组模型是IES优化的传统难点主要在于它的热电比可行域是个多边形区域。理想化处理是定热电比模型发热量 发电量 × 热电比。但实际运行中CHP的电出力区间受限且在不同负载率区间的热电比不同。原论文用的是线性化的可行域约束方法即把电力出力和热出力限制在一个凸多边形里。我在复现时建议采用这个方案CHP运行时满足最小/最大电出力约束热出力由“电出力-热出力可行域”多组线性不等式描述。这比定热电比模型更贴近原文也是审稿人比较认可的处理方式。P2G设备的模型相对简单输入是电力输出是天然气热值效率一般在50%~65%。约束条件包括输入功率上下限、爬坡约束以及产气量的连续/离散状态耦合。但有个容易忽略的变量——CO2消耗量它必须和碳捕集设备的产碳量联动( Q_{CO2,P2G}(t) \gamma_{P2G} \cdot P_{P2G}(t) )这个系数 (\gamma_{P2G}) 通常根据P2G制甲烷的化学计量比和效率计算。一旦漏掉这个约束P2G几乎等于凭空产气碳排放成本计算结果会明显偏小。碳捕集设备的模型核心是捕集率与能耗的权衡。设备从烟气中分离CO2需要消耗一定的电和热再生能耗捕集率越高单位捕集成本越高通常用线性或分段线性关系描述( Q_{CCS}(t) \eta_{CCS}(t) \cdot E_{in}(t) )其中 (\eta_{CCS}) 是捕集效率与再沸器热耗相关。这里的变量设置要特别注意捕集到的CO2有两个去向——供给P2G的 ( Q_{CO2,toP2G}(t) ) 和外送封存/利用的 ( Q_{CO2,storage}(t) )两者之和等于捕集总量。储碳罐的容量约束和初始/末端状态约束也别漏否则容易让模型“白嫖”初始值。储能设备建模储电、储热、储碳的基本逻辑一致状态量是一个时段的容量更新方程是容量 上一时段容量 充入量 × 效率 - 放出量 / 效率。注意这里充放效率不要写成同一个方向否则能量守恒会出现偏差。储碳罐的充放对象是CO2质量充放效率通常按100%处理但储存有容量上限和日始末平衡约束。2.3 目标函数的构成成本项分解与标准化双目标问题的两个目标函数一个是碳排放成本一个是运维成本二者量纲不同、量级可能差一个数量级所以epsilon约束法比加权法更合适。碳排放成本这个目标函数要仔细区分“碳排放”和“碳排放成本”两个口径。直接排放包括向外部电网购电在发电侧产生的等效碳排放用区域电网排放因子折算、CHP和燃气锅炉燃烧天然气产生的排放、碳捕集设备未捕集到的残余排放。同时要扣除捕集并封存/利用的部分。于是净碳排放量的表达式是( E_{total} \sum_{t} ( \sigma_{grid} \cdot P_{buy}(t) \lambda_{gas} \cdot F_{CHP}(t) \lambda_{gas} \cdot F_{boiler}(t) - Q_{CCS}(t) ) )碳排放成本就是把上述量乘以碳价tCO2价格。注意很多复现代码只算了CHP和燃气锅炉的排放忘了折算购电的间接排放这个误差在系统风电渗透率低的时候会非常大。我自己排查过这类问题算出来碳成本比原文结果低20%多后来追根溯源才发现是电网侧等效排放系数没有计入。运维成本的目标函数则包括购电成本、购气成本扣除P2G产气的替代收益、CHP和燃气锅炉的燃料成本、P2G与碳捕集设备的运行维护单位成本、储能设备充放的老化折算成本以及启停成本如果模型是含0-1变量的混合整数规划。运维成本里最容易被忽略的是碳捕集设备的再生热耗对应的燃料成本——它不开炉子但要从热母线上“分走”一部分蒸汽。如果热母线平衡里没体现这部分热耗运维成本就会偏低而且热功率平衡也会失衡。在搭建目标函数之前我建议先把决策变量清单列出来包括连续的功率变量各设备出力、储能充放功率、整数变量设备启停、状态变量储能容量。决策变量数量决定了模型的规模。本文这套系统典型24小时调度、1小时步长连续变量大体在200~300个0-1整数变量在50个左右属于中小型混合整数规划MILP可以轻松处理的规模。3. 从单目标到双目标epsilon约束法的求解思路与代码落地3.1 为什么不用线性加权法原作者在论文里用的是线性加权法即把两个目标函数乘以权重后相加得到一个总目标函数。这样做的好处是简单一次求解就能得到一组解但严重的问题是当两个目标的量纲和量级相差较大时权重系数 (\omega_1)、(\omega_2) 的物理意义不直观需要做归一化处理更关键的是线性加权法求解多目标问题只能得到Pareto前沿的凸包部分对于非凸问题无法获得完整的Pareto前沿。IES运行优化问题由于存在设备启停的0-1变量和分段线性约束恰好在某些区间表现出非凸性单纯加权法会漏掉部分有意义的运行策略。epsilon约束法的思路完全不同保留一个主目标比如总运行成本即运维成本把另一个目标碳排放成本由目标函数降级为不等式约束( E_{total} \leq \epsilon )。之后连续改变epsilon的值从单目标最优碳排放值逐步增加到碳排放最宽松的值反复求解新的优化问题就能得到一系列Pareto最优解。这个方法对凸与非凸问题都适用而且epsilon的取值有明确的物理含义——就是碳排放的上限决策者可以直观地判断“我最多允许排这么多碳”。3.2 epsilon约束法的三个关键预处理步骤第一步求解单目标问题只最小化运维成本记此时的最优碳排放量为 ( E_{min}^* )。这就是Pareto前沿的最低碳排放点之一。第二步求解另一个单目标问题只最小化碳排放成本记此时的最优运维成本。注意最小化碳排放时得到的运维成本通常非常高这两个端点就是epsilon的范围边界。第三步把 ( E_{min}^* ) 到第一个单目标下自然产生的最大排放量其实是运维成本最小时对应的 ( E_{max}^* )区间等分成N份我建议20份起步30份效果已经很好得到N1个epsilon断点逐一求解。这个流程有一个问题容易踩坑只跑两个单目标端点拿到的碳排范围往往不是Pareto前沿的真实范围。比如最小化碳排放时系统会大量启用P2G和碳捕集而P2G耗电来自外部电网时间接排放大幅增加反而可能让净碳排放不降反升。因此epsilon区间的上界建议用“最小运维成本下对应的碳排放值”下界用“全局最小碳排放值”而不是简单地拿两个单目标极值。相比之下加权法的权重区间没有这种天然的物理解释。3.3 Matlab中epsilon约束法的主循环框架整个主循环的本质是对每个epsilon值调用一次求解器。% 关键流程示例 % 1. 先求解两个端点单目标问题获得epsilon范围 [cost0, emit0] solveSingleObjective(cost); % 最小化运维成本 [cost1, emit1] solveSingleObjective(emission); % 最小化碳排放成本 eps_min emit1; % 取较小端 eps_max emit0; % 取较大端 N 25; eps_values linspace(eps_min, eps_max, N); pareto_solutions []; for i 1:length(eps_values) epsilon eps_values(i); % 把碳排约束塞进原问题 [x, cost_val, emit_val] solveWithEpsilon(epsilon); pareto_solutions [pareto_solutions; cost_val, emit_val, x]; end这里要特别提醒epsilon约束写法有两种常见形式。第一种是把碳排放约束加入原约束集保持目标函数为运维成本第二种是把碳排放约束作为软约束在目标函数里加惩罚项。复现时我强烈建议用第一种硬约束。因为软约束的罚因子需要反复调参而且在某些epsilon区间惩罚项会与运维成本产生“共振”导致解不稳定。硬约束虽然可能增加求解时间但在逻辑上是最干净的。还要注意一个细节epsilon值应该按从紧到松或从松到紧的顺序遍历但每次求解之间没有热启动依赖。为了提高速度你可以用上一次求得的决策变量值作为本次求解的初始点尤其在处理启停变量时热启动能减少分支定界树的搜索深度。不过热启动也有风险如果epsilon变化太大之前的解离最优解太远热启动反而带来误导。我实践下来epsilon变化幅度小时热启动有效幅度大时冷启动更稳。4. Matlab代码实现文件结构、关键函数与参数设计4.1 代码的整体文件架构这类复现项目最适合的代码结构是“分模块脚本主入口”方式而不是把所有代码堆在一个.m文件里否则等你需要在epsilon循环里反复改动参数时会非常痛苦。我建议的文件结构如下main.m主入口脚本定义系统参数、调用单目标求解函数、执行epsilon循环、绘图data_para.m包含所有设备参数、电价、气价、负荷数据的初始化脚本build_problem.m构建优化问题的约束矩阵与目标函数表达式输出模型结构体modelsolve_milp.m调用求解器求解混合整数线性规划返回结果结构体calc_emission.m给定决策变量计算碳排放量的独立函数calc_cost.m给定决策变量计算运维成本的独立函数plot_results.m绘制电/热/气/碳平衡图与Pareto前端。为什么要把碳排放量和运维成本单独拆成函数因为epsilon约束法的核心就是反复评估这两个指标单独成函数方便你打印、核对、画图。尤其调试时我需要每次求解完立刻对比calc_cost和calc_emission的输出与求解器内部目标函数值是否一致如果存在偏差说明目标函数里有项写漏了。4.2 核心建模代码段的写法与细节以Yalmip工具箱为例如果没有Yalmip用独立求解器的接口也可以但Yalmip的约束表达更直观% 决策变量 P_chp sdpvar(1, 24); % CHP电出力 H_chp sdpvar(1, 24); % CHP热出力 u_chp binvar(1, 24); % CHP启停 P_p2g sdpvar(1, 24); % P2G输入电功率 Q_ccs sdpvar(1, 24); % 碳捕集量 E_cap sdpvar(1, 25); % 储碳容量注意有25个点 ... % 目标函数运维成本 Objective sum(price_e .* P_buy) sum(price_g .* F_buy) ... sum(om_chp .* P_chp) sum(om_ccs .* Q_ccs); % 碳排放约束epsilon约束循环内生成了再添加 Constraint [Constraint, total_emission epsilon]; % 求解 ops sdpsettings(solver, gurobi, verbose, 2, showprogress, 1); result optimize(Constraint, Objective, ops);这个代码段里有个很重要的变量E_cap的索引长度是25而不是24。因为储能容量是连接t时刻和t1时刻之间的状态变量你需要用第1个点作为初始容量、第25个点作为末端容量一天24个时刻正好需要25个断点。直接用24维变量就会出现容量循环引用错位很多复现代码的储能平衡约束看起来正确但实际逻辑是错的多半就是这个索引对齐问题。P2G和碳捕集设备的耦合约束正确写法是% CO2物料平衡碳捕集量 P2G消耗量 储碳净注入量 Q_ccs(t) - alpha_p2g * P_p2g(t) - (E_cap(t1) - E_cap(t)) - Q_ccs_sale(t) 0;这个约束里的 ( Q_{ccs_sale} ) 是指外售或封存的CO2可以产生收益也可以当作纯变量不约束去向。在文章里它通常被放进碳排放成本目标里作为负项碳利用收益。如果复现代码里没有这一项P2G所需CO2被强制全部由储碳罐提供可能造成储碳罐容量约束被不合理突破需要特别注意。4.3 参数整定要注意的几个数值电价和气价序列直接决定成本曲线形态。要注意的是我国现行电价是分时峰谷电价原论文里用的可能是两段式或者三段式。复现时建议按原文献价格表填如果找不到可以从典型场景里取日峰平谷三时段电价峰平谷比值在1.6:1:0.4左右。气价按热值计价元/kWh一般在0.3~0.4元/kWh之间这样购电和购气两种能源输入在成本和碳排放上都处于合理区间模型结果不会一边倒。P2G效率是所有参数里对结果影响最大的一个。电解水制氢效率大约70%甲烷化反应效率大约80%总效率在55%~60%。但注意总效率的分子分母口径不同——分母是输入电能分子是产出天然气的热值。所以P2G每消耗1kWh电力大约只能产出0.55~0.6kWh热值的天然气。有的论文用了0.6以上甚至0.7结果P2G的产气收益会显得过高导致优化结果严重偏向P2G而偏离实际。碳捕集设备的能耗比能耗/捕集量通常在0.3~0.6 kWh/kgCO2左右加上热耗约2~4 MJ/kgCO2。如果按电耗热耗折合一次能源计算单位捕集成本大约在40~120元/tCO2。设置这个参数时要注意和碳价匹配如果碳价是50元/t而捕集成本高于碳价那系统的碳捕集设备在经济上就是“亏损”的除非有政策约束或P2G需要CO2原料否则优化结果会倾向于关闭碳捕集。这一点对理解Pareto前沿形态很关键。4.4 求解器选型与工具链配置Matlab下求解MILP的主流方案是YalmipGurobi或YalmipCPLEX再退一步可以用Matlab内置的intlinprog。但对于带大量0-1变量的问题内置求解器的性能明显不够。我之前在一台普通工作站上测试同样规模的模型内置求解器大约要跑30秒到1分钟Gurobi只需要2~3秒。在epsilon约束循环里这个差距会被放大25倍以上。如果不想买商业求解器也可以用开源的CBC求解器Yalmip支持得很好。CBC对中小规模MILP完全够用但要注意MIP gap设置默认gap可能设在1e-4对运行优化问题来说解的质量已经足够但速度太慢我一般把它放宽到1e-3精度几乎没有损失计算时间可以减少一半以上。追求速度时不妨试试。还有一个实用建议是先把模型写成完全线性的。原模型里的碳排成本和运维成本必须是线性的。如果你把P2G效率写成随负载率变化的非线性函数求解器就要走非线性规划epsilon约束法在非凸情况下容易卡出不理想解复现时尽量保持线性化处理。5. 复现过程中踩过的坑与排查方法5.1 epsilon循环求解结果不稳定的排查我在跑epsilon约束时遇到过一种情况epsilon值变化很小但解出来的运行策略完全不同Pareto前端折断成好几段看起来像“异常断点”。排查下来问题出在两处。第一个原因是某些设备启停变量的耗时状态与设备出力的耦合约束出了问题。原模型里设备启动需要满足最小开机时间或最小停机时间但我在初始代码里把启停变量和出力下限约束写松了导致某些解里设备出现1小时级别的频繁启停。Gurobi求解时会利用状态变量的整数性质去优化但一旦某一段epsilon下启停模式改变整条曲线的成本结构就跳变。解决办法是把设备的最小启停时间约束加上至少2~3小时运行策略的连续性明显改善。第二个原因是碳排放量的计算口径里有非线性项。比如电网购电排放因子在峰时和谷时不一样时 ( E_{grid}(t) \sigma(t) \cdot P_{buy}(t) )这是分段线性函数可以引入辅助变量做线性化。如果没做线性化而用到了非线性表达式求解器内部可能走的是MIQP混合整数二次规划路径结果震荡明显。把这些项全部线性化之后Pareto前端就平滑多了。5.2 碳排放量与成本量纲不匹配导致epsilon范围异常如果你设置的epsilon范围是从-50到1000基本都是单位问题。比如碳排放计算用的是吨碳排放约束的epsilon设的是千克或者碳排放成本目标里用的是元而运维成本单位是万元那很可能某个目标的量级完全不匹配导致epsilon约束要么永远不满足要么永远松弛Pareto前端退化成一个点。我踩过一次碳排放量按kgCO2计算但碳价是按元/tCO2设定的结果算出来的碳排放成本目标函数值比预期大了整整1000倍。最后在预算约束里发现运维成本反而被跳过了Pareto前端完全变形。所以建议在代码开头就做一次单位自检打印两个端点单目标求解的碳排量和成本值目测量级是否合理。IES系统典型24小时的碳排放量一般在一百到几百吨CO2量级运维成本在几万到几十万元量级。看到这个范围内的数值再开始跑epsilon循环。5.3 求解时间过长与MIP gap的关系epsilon循环N25时25次MILP求解单次时间如果超过5秒总时间就会拉长到几分钟如果是更复杂的场景比如多场景随机优化时间会涨到小时级。提升速度的几个方向第一减少整数变量数量。CHP启停变量可以按“只在负荷变化时段候选”的方式缩减但这一般不推荐因为容易漏掉最优策略。更可靠的是把储能设备的充放状态用连续变量表示因为储能的充放效率折算直接通过目标函数调节这样省掉大量0-1变量。只要储能的充放不同时发生由成本方向自动约束即可不需要显式0-1互斥变量。第二设置Gurobi的MIPFocus参数和Heuristics参数。MIPFocus设为1可以更侧重寻找可行解Heuristics越高越早找到好解但可能牺牲最优性证明速度。批量跑epsilon循环时我一般把MIPFocus设为1Heuristics设为0.2~0.3能在总时间内拿到质量足够好的Pareto前端。第三并行化。如果你的Gurobi授权支持多线程Matlab的parfor可以并行求解多个epsilon点的优化问题。注意parfor循环里不能有数据竞争每轮求解用的模型参数要复制一份。这样25个点可以一次同时求解总时间基本等于单个点的求解时间。5.4 结果绘图时的常见错误与解读方法复现完成后画Pareto前端是检验结果正确性的第一道关卡。正常情况下Pareto曲线应该是一条单调递减或递增的凸型曲线碳排放成本越高运维成本越低碳排放成本越低运维成本越高。如果曲线出现上升段或明显的凹折说明epsilon约束的某个区间内模型不可行或求解器没找到全局最优。画各类能量平衡图时我的习惯是把电、热、气、碳四条链路的“源-荷-储-网”区分开来各画一张堆叠柱状图图上同时叠加对应设备的功率曲线。这样做的好处是能快速发现功率平衡约束中的“幽灵漏洞”。比如如果某个时刻热平衡图上CHP热出力燃气锅炉热出力远远大于热负荷热耗那一定是有热负荷数据或P2G热耗项漏写如果电平衡图上某时刻总发电比总负荷多很多却不多见弃电项那也是平衡约束有问题。6. 复现这个项目的几点真实体会做完这个复现最深的体会是原论文的难点不在“运行优化”本身而在P2G和碳捕集设备之间的耦合关系以及把物理过程翻译为线性约束的建模能力。一个看似简单的 ( CO_2 \rightarrow CH_4 ) 化学路径落到数学上就是碳平衡、氢平衡、热平衡三条等式同时作用。如果只盯着设备功率上下限和电热平衡代码能跑通但结果和论文对不上审稿人一眼就能看出模型简化过头了。从方法角度看epsilon约束法相比加权法的优势在实际程序里体现得非常明显。我可以在epsilon从紧到松的过程中观察系统的行为变化刚开始的时候CHP电出力被压到很低碳捕集设备几乎满负荷运行P2G启动但受限于CO2供应随着epsilon放宽CHP的出力逐步升高碳捕集负荷下降P2G仍然运行但更多是为了风电消纳而不是碳捕集。这种“行为模式切换”是加权法很难直观呈现的也是epsilson约束法在IES领域被广泛使用的原因。最后分享一个小的实用技巧epsilon循环结束后记得把所有Pareto点对应的决策变量保存下来不要只保存目标函数值。因为后处理时需要从这些解里提取“某epsilon水平下P2G的投入率”或“碳捕集设备的运行时长”这类信息用来做敏感性分析。很多复现代码只画了Pareto前端就交差但论文审稿意见里经常要求补充两个目标对应特定运行模式的解释性图表。保存全变量你的回旋余地就大多了。从实际体验来说调通模型不难难的是把每个解的运行策略看明白、讲清楚。如果跑出了正确的Pareto前端并且能解释清楚为什么曲线在某个区间斜率变缓通常是碳捕集设备能耗占比上升带来的边际效益递减这篇Energy一区文章的复现就算真正到位了。