1. 项目概述与论文复现的核心价值

先说说这个题目本身。P2G（Power to Gas，电转气）是电气综合能源系统里这两年绕不开的一个关键环节，它把富余风电、光伏的电能转化成氢气甚至合成天然气，让电力系统和气网系统产生双向耦合。我当初选择复现这篇论文，核心目的是想把这套"计及P2G厂站"的规划模型彻底吃透，并且在Matlab环境里跑通完整的优化求解流程——从天然气网络建模、电力系统约束、P2G厂站运行特性，到规划方案的迭代寻优，每一步都要落地成可执行的代码。

为什么值得复现？这类论文最大的价值不在于最终的规划结果，而在于建模思路：它把电-气耦合从单一的电转气购能问题，升级成了考虑厂站内部运行约束（电解槽效率、储氢罐容量、甲烷化反应热量平衡）的一体化规划问题。对于准备做综合能源系统方向的学生，或者刚接触能源互联网优化的工程师，把这篇论文的代码复现出来，等于一次性打通了电力系统优化调度、天然气网络稳态分析、混合整数线性规划（MILP）求解三条技术线。

我在复现过程中的体会是，这类论文的难点不在数学推导，而在工程落地。学校给的论文材料往往只给最终模型和算例参数，中间过程——比如天然气管道怎么线性化、网损怎么处理、NLP和MILP怎么切换——需要自己大量补课。下面我把整个复现过程拆开讲，从建模到代码实现再到踩坑记录，给你一条相对平整的路。

2. 系统建模：电网、气网和P2G厂站的协同表达

2.1 电力系统部分怎么建

电力系统在规划模型里一般保留节点功率平衡、机组出力上下限、爬坡约束这几个核心模块。复现时我确认了原论文的假设：规划层不考虑暂态过程，用直流潮流模型近似交流潮流，这么做既保留了网络拓扑约束对规划结果的影响，又让模型整体保持线性。

需要特别留心的是P2G厂站作为负荷接入时怎么建模。P2G的本质是一个大功率电负荷，所以任何含P2G的节点，负荷不再是固定值而是决策变量的一部分。这在实现上要改动现有的节点功率平衡方程：

% 节点功率平衡，P_P2G是决策变量 % 原方程: sum(Pg) - sum(Pload) = sum(Pij) % 修改后: sum(Pg) - sum(Pload) - P_P2G = sum(Pij) balance_eq = sum(Pg_idx) - sum(Pload_idx) - sum(P2G_idx) == sum(flow_idx);

很多同学复现时容易在功率平衡里忘掉P2G这一项，直接导致规划结果里P2G容量越大，系统不平衡量越大，仿真结果完全失真。

再就是网损。直流潮流模型一般忽略网损，但在规划问题里全网总有功平衡往往对计算结果影响很大。我看的这篇论文采用的是迭代修正法：先求解不计网损的模型，再根据潮流结果计算网损，作为已知量回调到下一个迭代中。这种做法实现起来简单，但要注意迭代的收敛条件设置。我在实测中用的是网损前后两次变化小于0.1%作为停止条件，通常在3到5轮就能收敛。

2.2 天然气网络稳态模型

气网部分是这个领域公认的难点。天然气管网的核心变量是节点气压和管道流量，两者之间的非线性关系让问题变得很棘手。常见的处理手段有两种：一是分段线性化（piecewise linearization），把Weimouth方程按流量区间分段逼近；二是直接在大规模MILP里用增量线性化方法嵌入。

原论文采用的是增量线性化方法，这是我复现时印象最深的一块。Weimouth方程本身是：

% 管道流量与两端压力满足非线性关系 % F_ij = C_ij * sqrt(pi^2 - pj^2) % 令 H_ij = pi^2 - pj^2，则流量关于H_ij是根号关系 % 线性化: 将H_ij取值范围分成N段，每段内流量线性逼近

分段数量是精度和计算量的折中。我试过5段、8段、10段，发现对规划问题来说8段已经足够，分段数超过10段之后计算时间几乎翻倍，而结果差异不到0.5%。如果论文没给具体分段数，我建议直接采用8段作为默认参数，这也是该领域文献里最常见的选择。

气源、储气罐、负荷这三类节点的约束相对直接，主要注意气源出力的上限和节点气压的上下限。节点气压范围在规划问题里经常被误设为固定值，实际上气网节点气压允许在合理区间内浮动，比如0.9到1.2倍的基准值。如果气压约束过紧，会过度限制气网消纳P2G产气的能力。

2.3 P2G厂站内部的"黑箱"变"白箱"

这篇论文区别于普通电-气耦合研究的最大亮点，是把P2G厂站内部过程展开了。P2G不是一个简单的输入电、输出气的黑箱，而是一条由电解槽、储氢罐、甲烷化反应器、气体压缩机四部分组成的工艺链条。

电解槽环节的关键约束是额定容量和运行范围。电解槽的输入功率不能低于某个比例，否则电解效率急剧下降，所以一般设一个最小运行功率约束（比如额定功率的20%）。储氢罐的作用是缓冲电解产氢和甲烷化耗氢之间的时间不匹配，它的状态方程是一个离散时间递推式，容量约束和初始/终态储量约束都必须加进去。

甲烷化反应环节存在热量平衡问题——甲烷化是强放热反应，建模时要在P2G功率输出和运行温度之间做一个温度约束的简化处理。不过大多数规划类论文并不真正求解热平衡，而是直接用氢转甲烷的转换效率乘以输入氢量得到产气量。如果复现时想更严谨，可以加一个温度惩罚项，但我个人认为对于年度规划问题来说意义不大，徒增非线性。

气体压缩机在P2G厂站模型里往往是最容易被忽略的一块。P2G产气压力通常低于天然气输气管网的压力等级，必须经过压缩机升压才能注入气网。压缩机本身消耗的功率虽然占比不大（约占P2G总耗电的2%到5%），但在规划模型里如果不计这部分自耗电，P2G的净效率会被高估。我建议至少按压缩比和流量做一个线性化的功耗估算，别完全省略。

3. 规划模型构建与求解器选型

3.1 目标函数的三层结构

原论文的目标函数是典型的多层规划架构：投资成本+运行成本+环境成本。我复现时把目标函数拆成了三层，方便后续做敏感性分析：

• 投资成本层：P2G厂站各设备的单位投资成本乘容量再乘年值系数，注意设备寿命不同，折算系数也不同。电解槽寿命一般按10到15年算，甲烷化设备按20年算，别统一套一个系数。
• 运行成本层：包括购电成本、购气成本、机组启停成本。这里购电成本要区分分时电价，论文算例里通常给的是峰平谷三段电价。
• 环境成本层：按碳排放量折算成惩罚费用。需注意碳价参数设置，原论文一般会说明基准碳价是多少，复现时务必核对单位——是元/吨还是元/千克，弄错一个量级整个结果全乱。

这里有个心得：论文的原始算例数据不同，成本项权重差异很大。复现前先把目标函数各成本项的数量级算一遍，如果发现某一项比另一个项小几个数量级，大概率是单位问题而非真实差异。

3.2 约束条件的层级拆分

规划模型本质上是双层问题：投资决策（长期）与运行决策（短期）。复现时如果直接用一个大规模MILP求解器硬解整个模型，计算规模会非常恐怖——因为运行层要模拟365天×24小时的调度过程，变量总量轻松上万。

原论文这里的处理思路值得学习：把规划问题拆成主问题和子问题的迭代式。主问题是投资决策，输出P2G厂站建设方案；子问题是给定投资方案后的年度运行优化，输出运行成本和可行域反馈。主-子问题之间通过Benders分解的思路交互。

我在复现中没有写完整的Benders分解，而是采用了更工程化的启发式迭代方案：先给一个初始P2G容量猜测值，求运行子问题得到该方案下的最优运行成本，把第一轮结果里被触发的容量瓶颈约束提取出来，用于修正下一轮的投资方案，如此迭代三到四轮便可收敛。

必须说明的是，这个简化方案牺牲了严格的全局最优性。如果审稿要求严格的最优解，正版的Benders或直接MILP求解器是必需的；但如果只是做工程方案分析，这个迭代法的结果已经足够可靠，而且速度快一个数量级。

3.3 求解器选型与性能表现

Matlab环境下求解MILP问题，我比较过几套方案：自带的intlinprog、 YALMIP+Gurobi、 YALMIP+CPLEX。实际测试结果让我有点意外：intlinprog在中小规模算例（节点数少于30）表现尚可，但一旦进入IEEE 39节点或118节点级别的气电耦合系统，intlinprog的求解时间和数值稳定性都肉眼可见地变差。

Gurobi在MILP求解上的性能优势非常明显，尤其是大量二元变量的场景。以我复现的30节点电网加20节点气网算例为例，Gurobi求解时间约120秒，intlinprog则耗了近800秒，而且Gurobi的解质量（目标值更优）更好。

给一个小建议：复现这类论文，如果资金宽裕，优先用Gurobi。如果只有Matlab基础工具箱，也完全可以跑通，只是要把算例规模控制在合理范围内，并设置合适的求解精度和最大迭代次数。

4. Matlab代码实现：从框架到核心函数

4.1 代码整体架构设计

我见过不少同学复现代码时喜欢把所有逻辑写在一个几百行的主脚本里，面向过程的写法虽然直白，但一旦需要调节参数或换算例，就变得寸步难行。我这次复现采用了模块化设计，简单说就是数据、模型、求解、结果四层分开：

项目根目录/ ├── data/ % 算例数据，按系统分类存放 ├── models/ % 模型构建函数 ├── solver/ % 求解器调用封装 ├── results/ % 结果输出与图表生成 └── main.m % 主入口，整个流程编排

这个架构的收益是在换算例时显现的——只需新增一个data子目录中的数据文件，代码零改动即可运行新的系统。如果你手头没有现成的电-气耦合标准算例，用MATPOWER提供的电力系统数据，搭配一个自建的气网数据文件也能拼凑出可用算例，不必非要获得原论文的配套数据。

4.2 关键数据结构设计

Matlab编程里最容易被忽视的是数据结构设计。我发现用结构体数组按"对象"组织数据，比散落的命名变量清晰得多。比如电网数据可以这样组织：

grid.bus = struct('id', [], 'type', [], 'Pg', [], 'Pd', [], 'Vmin', [], 'Vmax', []); grid.line = struct('from', [], 'to', [], 'R', [], 'X', [], 'capacity', []); grid.gen = struct('id', [], 'bus', [], 'Pmin', [], 'Pmax', [], 'ramp', [], 'cost', []);

P2G厂站的数据结构更复杂一些，需要包含四个子模块的参数：

p2g.electrolyzer = struct('capacity', [], 'eta_elec', [], 'Pmin_ratio', [], 'inv_cost', []); p2g.h2storage = struct('capacity', [], 'init_level', [], 'final_level', [], 'inv_cost', []); p2g.methanation = struct('capacity', [], 'eta_meth', [], 'Q_consume', [], 'inv_cost', []); p2g.compressor = struct('elevation_ratio', [], 'power_coef', [], 'inv_cost', []);

我踩过的一个坑是：初期把所有数据分散在不同变量里，结果跑大规模算例时内存管理混乱，经常出现变量名写错但程序不报错的情况（因为Matlab对变量名检查不严，只是逻辑错了）。用结构体之后，至少能在视觉层面对数据关系一目了然。

4.3 模型构建核心流程

建模型的过程我分为五步走，每一步都通过测试函数验证正确性后再进入下一步：

第一步，读取算例数据并预处理。这一步的关键是节点编号的对齐——电网和气网的节点编号体系不同，必须在数据层建立映射关系，论文里通常给的算例图数据已经标好，直接用即可。

第二步，构建电网模型。直流潮流中关键是导纳矩阵B的形成和节点分类（平衡节点、PV节点、PQ节点）。给定P2G厂站接入的新节点编号后，在B矩阵中插入对应行列，注意新增节点的基准电压和基准功率要与全系统一致。

第三步，构建气网模型。这部分是代码里最容易写错的。管道流量线性化需要预先计算分段点和斜率，我写成了一个独立函数：

function [seg_points, slopes, intercepts] = linearize_weimouth(pipe_C, p_max, p_min, n_seg) % 输入: 管道常数C, 压力上下限, 分段数 % 输出: 每个分段的端点、斜率、截距 H_max = pipe_C^2 * (p_max^2 - p_min^2); H_seg = linspace(0, H_max^0.5, n_seg+1).^2; F_seg = sqrt(H_seg); % 计算各段斜率 slopes = diff(F_seg) ./ diff(H_seg); intercepts = F_seg(1:end-1) - slopes .* H_seg(1:end-1); seg_points = H_seg; end

第四步，把P2G厂站四个模块全部纳入模型。这一块我建议画一个简易的能流框图（在纸上画即可，不必画进代码），明确每一级转换的能量流方向，再逐一写约束。

第五步，组装目标函数和全部约束，交给求解器。

第五步是整个过程中最耗时的环节。我第一次组装模型时，光约束数量就遇到上百条，Matlab命令行窗口里报错信息满天飞，后来学会了一个技巧：每加一组约束后，立即求解一个简化版模型（只有该约束相关变量的固定值），验证可行性。这样做问题定位很快，基本不用从头调试。

4.4 求解封装与结果后处理

求解器封装我写成了通用接口，这样可以在intlinprog和Gurobi之间无缝切换：

function [x, fval, exitflag] = solve_milp(model, use_gurobi) if use_gurobi % 转换为Gurobi的输入格式 result = gurobi_optimize(model); x = result.x; fval = result.objval; exitflag = result.status; else options = optimoptions('intlinprog', 'Display', 'final', ... 'MaxTime', 1800, 'RelativeGapTolerance', 0.01); [x, fval, exitflag] = intlinprog(model.f, model.intcon, ... model.Aineq, model.bineq, model.Aeq, model.beq, ... model.lb, model.ub, options); end end

结果后处理这块，我建议一定做三张核心图：第一张是不同P2G容量方案下的总成本柱状图，第二张是典型日的电功率和气功率平衡曲线，第三张是P2G厂站内部能量流桑基图（用Matlab绘图函数手动实现）。这三张图是论文复现成果最直观的呈现方式。特别是"不同P2G容量下的成本曲线"这张图，几乎可以一眼看出规划方案的最优点——总成本最低处对应的P2G容量就是最优容量。

5. 常见问题与排查技巧实录

5.1 求解不收敛或收敛极慢

这是复现此类论文最普遍的问题。我在调试时遇到过不少次模型不收敛的情况，原因是多方面的，但绝大多数指向同一个核心——约束条件过于激进。比如P2G年最大利用小时数设得过高，导致投资容量在运行层面无法收回成本，模型就会反复尝试调整投资方案，始终找不到可行解。

如果遇到收敛慢的问题，我强烈建议先检查以下几处：

5.2 线性化误差过大

增量线性化方法的误差主要源于分段点的选取。我最初用等间距分段，发现管道流量较大的情况下误差可以达到5%以上，在系统层面产生可感知的偏差。改进方式是让分段间距不再均匀，而是让单位区间内的流量残差保持一致，即"误差等分法"。具体实现思路是：先初步分段求出各段最大残差，再调整分段点使各段最大残差相等。

这种方法实现较为复杂，但在求解效率和精度之间能取得更好平衡。如果只是复现论文结果，用均匀分段并取8到10段一般就够用了。

5.3 算例结果与论文不一致的排查

这也是复现代码时非常常见的情况：代码能跑通，结果却和论文对不上。通常问题不在代码逻辑而在数据。建议按以下顺序排查：

第一步，检查单位是否一致。论文中天然气流量可能是立方米/小时、千克/小时、Mbtu/小时三种不同单位混用，换算关系搞错会直接导致几十倍的偏差。

第二步，检查基准功率和基准电压设置。ylq9综合能源系统研究中，电力和天然气的基准值通常设为100MVA和1.0MPa，如果基准值取错，整个标幺值体系都会偏移。

第三步，检查典型日的选取。原论文的规划周期是8760小时，但算例里的典型日可能是按季节或按峰谷时段浓缩出来的。复现时要明确原论文是用了完整的8760小时建模，还是用了若干个典型日乘以权重系数。这两者的结果会有不小差异。

第四步，检查P2G效率的计算方式。是"低热值效率"还是"高热值效率"？两者数值相差约10%左右，论文如果没有明确说明，复现时很容易出现几不可见的差异。

5.4 求解器数值稳定性问题

Matlab的intlinprog在处理具有不同量级数值的混合整数问题时，容易陷入数值病态问题——比如投资成本是千万量级，而运行成本是百万量级，二元变量的目标系数很小，导致MIP启发式搜索表现不佳。

解决办法是对模型中的系数进行归一化，把投资成本除以一个基准值（比如总投资上限），让所有目标项的量级集中在1到100之间。这个操作的原理是避免求解器内部处理跨量级数值时产生舍入误差。

别小看这一步，我遇到过一组算例归一化前后目标值相差4%的情况，而这个差异完全来自数值误差而非模型变化。

6. 复现过程的经验总结与后续扩展方向

整个项目复现下来，我最大的体会是：论文复现不是简单的"翻译代码"，而是对建模思路的再发现。原论文里一笔带过的许多假设——比如管道流量线性化的具体分段数、P2G厂站内部压缩机自耗电的处理方式——恰恰是工程实现中最需要斟酌的细节。真正跑通一轮完整流程之后，你对综合能源系统规划的理解深度，会远远超过只看理论推导时的水平。

如果后续想在这个基础上扩展，我个人建议几个方向：一是把P2G厂站模型换成更精细的电制氢全链条模型，加入电解槽的启停成本和动态效率曲线；二是在规划模型中加入不确定性因素，比如风电出力和电价的随机场景；三是把天然气网由稳态模型扩展为动态模型，考虑管道储气效应。这三个方向在当前的研究中都很热门，而且都是能在复现代码基础上做增量式修改完成的。

最后再分享一个经验：复现过程中请务必保存好每一个能运行的版本，并做好注释。你可能现在觉得某个中间版本没用，但三周后当你发现新模型的求解器行为变得诡异时，那个"旧版本"就是你回溯排查的最佳参照物。每次改动前，先跑一遍当前版本，记录目标值和关键约束的满足情况，再做修改。这套好习惯可以帮你省掉大量调试时间。