做综合能源的硕士生看到“计及P2G厂站的电-气综合能源系统规划”这种题目第一反应应该是这又是个要啃硬骨头的复现。P2GPower to Gas的本质是电转气把电网里富余的电能变成天然气让电力系统和天然气系统真正意义上双向互动。这篇复现的核心就是用Matlab搭建一套“电网气网P2G选址定容”的联合规划模型输出可跑的代码和清晰的求解逻辑。我当时拿到这个题目第一反应是找作者要源码结果跟大多数情况一样——没要到。自己吭哧吭哧复现了一周多踩了不少坑才把模型跑透。这篇记录就把整个过程拆开讲讲重点是代码怎么组织、非线性约束怎么处理、以及那些容易让模型“跑飞”的细节。适合正在做电-气耦合规划课题的硕士生、刚入门综合能源方向的研究生也适合想了解P2G建模要点的工程师。1. 先理解P2G为什么是电-气系统规划的关键1.1 电-气耦合的现实背景为什么规划要绑定在一起传统规划思路里电网和天然气网是各管各的电网规划只关心电源和线路气网规划只关心气源和管道。但现实早就不是这样了燃气轮机大规模接入后电网的出力上限很多时候不是由机组决定的而是由气网能供多少气决定的。反过来气网里加压站、调压站的驱动能源也可能来自电网。两个系统在运行层面已经深度耦合如果规划阶段还分开做很容易出现“电网规划落地后气网供不上气”或者“气网扩建后电网送不出电”这种尴尬局面。P2G在中间扮演的角色很特殊。它是目前唯一能让能量从电力系统反向流入天然气系统的设备相当于在两个系统之间装了一个双向阀门。电力系统富余时把多余的电变成甲烷注入气网储存气网缺气时这个“人造气源”顶上去。对规划来说P2G的存在意味着电网和气网可以互为备用比如某个区域电网输电阻塞与其扩建线路不如就地建一个P2G厂站把电转成气送走这往往更经济。所以规划阶段就必须把P2G作为一种可选方案放进优化框架里和线路扩容、气源扩建一起比选这才是“计及P2G厂站”的真正含义。另一个现实驱动力是可再生能源消纳。风电光伏大发时段电网消化不了弃风弃光成了常态。P2G恰好提供了一个消纳出口把本来要丢弃的电能转成气存起来既缓解了弃电压力又给气网增加了气源。我在自己的算例里配了30%比例的风电基准场景弃风率大概12%引入P2G之后可以降到5%以下。这就是这个课题最核心的应用场景高比例新能源下的多能互补。1.2 P2G厂站的内部结构与效率边界很多刚接触P2G的人以为它是一台“电变气”的单一设备其实不是。P2G厂站内部是两段工艺的串联第一段是电解水制氢碱性电解槽或者PEM电解槽靠电把水分子拆成氢气和氧气第二段是甲烷化反应氢气在催化剂作用下与二氧化碳反应生成甲烷和水。之所以非要费劲做成甲烷而不是直接用氢气是因为现有天然气管道、储气库和燃气设备对氢含量有严格限制氢气比例过高会腐蚀管材、影响调压器寿命而甲烷可以直接以一定比例掺混甚至完全替代天然气。两段工艺耦合下来整体能量转换效率大约在0.55到0.65之间我代码里取0.6。这个数字很关键它几乎决定了P2G在经济上什么时候划算。做个简单换算1 MWh电力按0.6效率转换得到0.6 MWh热值的天然气。标准天然气热值约36 MJ/m30.6 MWh对应2160 MJ折算下来就是大约60标准立方米的天然气。换句话说P2G厂站每消耗1 MWh电大概能产出60 Nm3的甲烷。如果电价正常水平0.5元/kWh光电费成本就相当于每立方米天然气8元多远高于常规气价所以P2G只有在弃电价格极低或者有额外政策补偿时才具备经济性。这也是为什么几乎所有的P2G规划文献都把场景设在“高比例新能源高弃风”背景下。这个效率边界还会直接影响规划结果P2G容量配置大了弃风消纳多但投资和运维成本也高容量配置小了又起不到明显作用。所以模型里效率、投资成本、弃电价格这些参数都要做敏感性分析不能只跑一组数据就下结论。2. 规划建模的完整拆解从物理网络到数学表达2.1 电网侧建模为什么规划场景下一定要用直流潮流电网模型在规划问题里没有太多悬念几乎清一色用直流潮流。原因很简单规划问题本身就带有0-1整数变量比如某个候选节点建不建P2G这已经让模型变成混合整数规划了。如果电网侧再用交流潮流支路潮流是电压相角的非线性函数整个问题就成了混合整数非线性规划MINLP求解难度和耗时直接上升几个数量级而且很难保证收敛到全局最优。直流潮流的数学形式是线性的支路有功潮流近似等于两端相角差除以支路电抗节点功率平衡也变成线性等式约束。整个电网侧约束组就变成了一组线性方程可以顺利嵌入MILP框架。直流潮流牺牲掉的精度主要在网络损耗和电压支撑方面对规划阶段来说误差通常在5%到10%区间完全够用。但要知道它的边界算出来的线路潮流和相角是近似值电压不是变量而是被忽略的假设。如果算例里含大量风电接入电压支撑问题不能靠直流潮流评估。我自己的做法是MILP先求出规划方案再用MATPOWER跑一次交流潮流做事后校验确认节点电压不越限。这个“规划粗算运行精算”的双校验流程写论文时也很加分。2.2 气网侧建模Weymouth方程是第一个真正的坑气网的建模比电网粗暴得多但也复杂得多。核心原因在于它的物理约束是强非线性的。节点流量平衡是线性的每个气节点的进气量加气源注入量等于负荷量。但管道潮流与节点压力的关系遵循Weymouth方程这是输气管道在稳态下的经典经验公式F_ij K_ij · sign(π_i - π_j) · sqrt(|π_i - π_j|)这里π_i是节点i气压的平方值K_ij是管道常数由管径、管长、摩擦系数决定。式子里那个sign符号和根号组合导致函数既非线性又不平滑。π_i本身又是变量所以约束里有变量相乘、开方、符号判断直接扔给求解器基本无解。初学者最常犯的错误是忽略符号处理。当节点i压力低于节点j时气压平方差为负流量方向要反向如果直接把根号里的数当成非负处理模型逻辑就错了。稳妥的做法是引入辅助变量Δπ|π_i-π_j|再配一个0-1变量表示流向但这样会增加整数变量数量。处理Weymouth方程的标准做法是分段线性化。因为π的取值范围是已知的——气压有上下限所以Δπ的运行区间也是已知的。把sqrt函数在这个区间内切成若干线段用增量模型或SOS2约束把它近似成线性关系。我在代码里每根管道切6段精度误差大约1%对规划结果的影响很小。分段数太少精度差太多则整数变量爆炸6到8段是实践中不错的平衡点。2.3 P2G厂站建模核心就四个约束P2G厂站在规划模型里的数学表达不复杂但细节容易出错。核心约束可以归成四类。第一选址与建设约束。候选节点i要么建设P2G厂站要么不建用0-1变量x_i表示。建设容量P_cap不能超过上限并且只有x_i1时才能有正容量写成约束就是0 ≤ P_cap ≤ x_i · P_cap_max。这里P_cap_max可以取一个合理上限比如该节点最大可接入容量。第二运行约束。P2G厂站在每个时段消耗的电功率不能超过建设容量。注意这里有个很多文献会忽略的细节P2G设备存在最小运行负荷比比如电解槽一般不能低于20%负荷运行否则效率急剧恶化甚至损坏设备。于是约束变成P_min · x_i ≤ P_p2g(t) ≤ P_cap其中P_min是额定容量的某个比例。第三能量转换约束。P2G的产气量与耗电功率成正比G_p2g(t) η · P_p2g(t) · GAS_CONVERT。这里GAS_CONVERT是单位换算系数把MWh电量转换成标准立方米天然气。按照前面的换算η0.6时1 MWh电对应60 Nm3气所以转换系数的数值要仔细核对。第四耦合位置约束。P2G厂站在电网侧是电负荷在气网侧是气源。建模时要明确它接入哪个电网节点、接入哪个气网节点两个接口的功率必须由同一个决策变量耦合起来。这四类约束都是线性的。只要把气网Weymouth方程线性化处理好整个模型就可以保持MILP结构。2.4 目标函数与约束体系汇总目标函数通常是总成本最小化包含三块P2G厂站投资成本、系统运行成本、其他扩展成本。投资成本按年化处理。P2G厂站一次性投资C_inv_total C_fix · x_i C_unit · P_capC_fix是固定建设成本场地、并联设备等C_unit是单位容量成本。年化要把寿命期内总成本折成等年值公式涉及贴现率。不同文献贴现率取6%到10%我代码里取8%整个目标函数里贴现率是影响投资决策的一个敏感参数改起来要谨慎。运行成本包括向上级电网购电成本、向气网购气成本、P2G厂站运维成本可能还有弃风惩罚成本。这里有个时间尺度统一的问题规划模型用典型日代表全年运行成本要乘以典型日代表的天数再年化。比如冬季典型日和夏季典型日分别代表90天和275天那么对应运行成本要乘以对应天数不能直接加总。约束体系汇总起来就是电网节点功率平衡、线路潮流上下限、发电机出力上下限、气网节点流量平衡、管道流量约束、节点气压上下限、气源出力范围、P2G厂站建设与运行约束加上所有变量的取值范围。整个约束体系写出来大约有上千行但借助Matlab的YALMIP工具箱构建起来比想象中省力。3. Matlab代码实现从模型到可运行代码3.1 非线性环节的处理思路对比在写代码之前先得把求解策略想清楚。气网Weymouth非线性的处理主流有三种思路各有适用场景。第一种是MILP分段线性化这也是我代码最终采用的方式。把Weymouth方程分段逼近之后整个模型变成混合整数线性规划能保证全局最优求解器可以直接给出结果。缺点是需要引入额外整数变量分段越细整数变量越多求解时间相应上升。但对规划问题来说这个代价是完全可以接受的。第二种是顺序迭代法也叫解耦迭代。先假设P2G规划方案固定电网问题与气网问题交替求解通过耦合变量迭代修正。这个方法的好处是每一侧都可以用很成熟的商业软件比如电网侧用MATPOWER、气网侧用专门的气网潮流工具但缺点是整个过程无法保证收敛到全局最优甚至可能出现振荡不收敛。第三种是锥松弛。把气压的平方作为新变量对非凸项做二阶锥松弛原问题变成二阶锥规划SOCP。精度比线性化高但YALMIP和Gurobi对SOCP的支持虽然很好遇到大规模算例时数值稳定性有时候不如纯MILP。而且锥松弛在边界条件比较紧的情况下可能松弛间隙较大结果需要通过后验收紧。三类方法对比如下方法优点缺点适用场景MILP分段线性化全局最优、实现简单分段误差、整数变量多规划问题主流方案顺序迭代可复用成熟软件无法保证全局最优运行模拟、在线调度锥松弛精度较高数据要求高、稳定性一般学术研究、特殊网络刚复现这个题目时不用纠结直接用哪个先走MILP线性化把完整链路跑通后面再对比其他方法。3.2 代码目录结构与运行环境准备我的代码组织很清晰五个部分main.m —— 主入口定义算例并调用建模与求解流程data/ —— 电网节点数据、气网节点数据、负荷曲线、风电曲线、P2G候选集model_build.m —— 用YALMIP定义全部决策变量和约束solve_and_report.m —— 求解并输出报告与图表utils/ —— 通用子函数比如数据读取、单位换算、Weymouth系数计算。环境配置方面Matlab版本建议R2020a以上YALMIP必须装好。YALMIP是Matlab环境下的建模工具箱能把优化问题用接近自然语言的方式写出来翻译成求解器能识别的标准形式。求解器推荐GurobiMILP性能目前没有任何替代品能比学校没有Gurobi许可证的话也可以用Cplex学术版或者先用Matlab自带的intlinprog跑小算例验证逻辑。安装YALMIP之后强烈建议先跑一遍yalmiptest确认建模环境与求解器能正常通信。很多“模型没问题但求解报错”的怪事最终查出来都是YALMIP路径配置或者求解器识别的问题。3.3 核心建模代码实现细节下面把我代码里P2G相关的核心建模部分贴出来这段是整个模型最有代表性的地方% 决策变量定义 x_p2g binvar(n_cand, 1); % 候选厂站是否建设0-1 P2G_cap sdpvar(n_cand, 1); % 建设容量单位MW P2G_pin sdpvar(n_cand, T); % 各时段输入电功率 P2G_gout sdpvar(n_cand, T); % 各时段输出天然气流量Nm3/h % 选址-容量耦合约束 M_cap 120; % 容量上限取120MW C [C, 0 P2G_cap x_p2g * M_cap]; C [C, 0 P2G_pin repmat(P2G_cap, 1, T)]; % 最小运行负荷约束20%额定容量 C [C, repmat(0.20 * P2G_cap, 1, T) P2G_pin]; % 能量转换约束效率0.6换算系数60 Nm3/MWh GAS_CONVERT 60; % 0.6效率下1MWh电对应60Nm3气 C [C, P2G_gout P2G_eff * P2G_pin * GAS_CONVERT];注意repmat(P2G_cap, 1, T)这种写法目的是把容量向量扩展到每个时段。很多人会在这里遇到维度报错其实是因为YALMIP对矩阵隐式扩展支持不如一些其他语言老老实实用repmat最稳妥。电网侧直流潮流约束% 节点功率平衡 C [C, sum(Pgen, 2) - P_load - P_p2g_elec B_theta * theta]; % 线路潮流约束 flow B_line * theta; % 直流潮流线性表达式 C [C, -Flow_max flow Flow_max];这里B_theta是节点导纳矩阵的虚部部分B_line是支路-节点关联矩阵与电抗对角阵的组合。构建之前先单独验证这两个系数矩阵可以用简单两节点系统对比手算结果。气网管道流量分段线性化的核心代码% 管道流量-压力关系的分段线性化增量模型 pieces 6; bp linspace(0, DeltaPi_max, pieces 1); flow_at_bp K_pipe * sqrt(bp); % 各分段端点对应流量 delta sdpvar(pieces, 1, full); % 分段增量 zbin binvar(pieces - 1, 1); % 分段激活指示变量 DeltaPi 0; Flow 0; for k 1:pieces DeltaPi DeltaPi delta(k) * (bp(k1) - bp(k)); Flow Flow delta(k) * (flow_at_bp(k1) - flow_at_bp(k)); end这段是示意结构真正的完整实现还要加分段激活约束防止增量变量“跳过”中间段直接填满后面的分段。缺失这条约束时线性化的分段可能被错位填充流量数值会出现物理上不合理的跳跃。我用implies配合zbin来保证分段顺序填充这是整个MILP实现里最容易被忽略的地方。3.4 求解器参数与收敛控制模型构建完成之后求解设置有几个关键参数直接影响成败。Gurobi的MIPGap默认是1e-4若规划算例规模大整数变量数量超过几千个可以把MIPGap放宽到1e-2求解时间能缩短一个数量级而规划结果的实用性差异微乎其微。除了MIPGap还要设置一个合理的时间上限比如7200秒防止求解器无限“钻”下去。代码里求解部分的骨架ops sdpsettings(solver, gurobi, gurobi.MIPGap, 0.01, ... gurobi.TimeLimit, 7200, verbose, 2); sol optimize(C, Objective, ops); if sol.problem 0 % 求解成功 recon_error check(C); else % 分析不可行原因 disp(sol.info); endcheck(C)是对所有约束进行残差检查的关键工具。求解结束后如果有约束残差明显偏离0说明数值稳定性有问题需要回到M值或者参数精度上去排查。这套“求解校验输出”三连流程能让绝大部分隐蔽问题现形。4. 算例复现与结果分析怎么判断代码对不对4.1 测试系统设计与参数设置复现论文时算例选择很讲究。太大跑不动太小说明不了问题。我默认算例用的是IEEE 14节点电网系统改进版加了一个风电场气网用6节点输气系统两个系统通过两个候选P2G厂站位置互联。风电场装机占系统总负荷约30%基准场景弃风率约12%。P2G效率0.6单位投资成本1.2万元/kW生命周期按20年贴现率8%。这个中等规模算例的优点在于MILP求解大概30秒内一定能跑完特别适合调试和复现。等逻辑验证通过后我再把算例替换成大型系统比如IEEE 118节点电网配合20节点气网代码框架完全不用动只需要替换数据文件。关于候选节点选择不是所有节点都能建P2G。我的经验是选三类位置一是风电送出受阻的节点二是电网负荷中心附近且靠近气网的节点三是气网气压偏低的节点。这三类位置P2G的价值体现最明显也最容易跑出有说服力的结果。4.2 从三个层次解读规划结果算例跑完怎么判断结果是否合理我习惯分三个层次看。第一层看规划方案本身P2G落在哪些候选节点容量配置多大。合理的结果应该落在我前面说的三类节点附近。如果P2G出现在莫名其妙的位置先检查候选集设置是否合理再检查目标函数是否写错——比如投资成本忘记年化会导致P2G容量虚高。第二层看系统运行指标弃风率、购电成本、购气成本、总成本的变化。引入P2G后弃风率应当显著下降但购气成本可能上升因为P2G产出的气也要计入系统气源。总成本是否下降取决于弃电价格与气价的相对高低。第三层看气网压力分布P2G注入点附近节点压力应当升高但绝对不能越上限。如果压力越限说明P2G容量配置超过了气网承受能力需要调整容量或增加压缩机。这一层校验能直接暴露“电侧OK但气侧崩了”的问题很多人漏掉这一步导致结果被审稿人质疑。为了验证P2G的规划效果我通常跑三个场景做对比不建设P2G的基准场景、允许P2G规划的场景、P2G强制配置高容量的场景。三个场景总成本和弃风率放一起规律一目了然。4.3 敏感性分析与合理性校验复现论文最怕结果和文献对不上。我的建议是先跑一个无P2G的基准场景把电网购电成本、气网购气成本分别与单独规划结果核对这一步能确认网络参数没偏。然后再开P2G逐步增加候选容量上限观察目标函数变化是否合理——总成本应该随容量上限先下降后平稳或略微上升如果出现先上升后下降的“凹”形多半是约束写错了。敏感性分析重点扫两个参数气价和P2G效率。气价从2元/m3扫到5元/m3时P2G经济性逐渐变差规划容量应该单调下降P2G效率从0.5提到0.7时规划容量应该增加。这种单调性检验虽然基础但能过滤掉绝大部分建模错误。我甚至建议把敏感性分析结果画成曲线放论文里比单一算例结果更有说服力。还有一个容易被质疑的细节P2G产气与气网负荷的时间匹配。P2G在风电大发时段大量产气此时气网负荷可能处于低谷多出来的气如果无法存储实际消纳效果就要打折扣。所以算例里最好加入储气罐或者允许P2G产气替代其他气源的比例约束这样结果更贴近工程实际。5. 复现一路踩过的坑问题排查速查5.1 模型不可行的排查路径复现过程中最高频的问题就是“模型无解”。第一次遇到时直接懵了怀疑是不是求解器坏了。后来发现90%的情况是约束维度写错或者耦合变量符号反了另外10%才是数值问题。我的排查流程很固定先用sol.info看求解器返回的状态如果是“Infeasible”就用YALMIP的check(C)定位残差最大的约束然后把Constraints列表逐个注释用二分法定位哪一组约束导致不可行。一般十分钟内能找到问题。最容易出问题的就三处电网节点功率平衡里的P2G耗电项符号、气网节点平衡里的P2G产气项符号、容量约束里矩阵维度不匹配。5.2 大M值的选择教训很深刻P2G容量约束里需要一个M值很多初学者的习惯是越大越好干脆写1e6。我在某次调试时就吃过这个亏——模型求解特别慢偶尔报数值警告后来把M值从1e6缩小到容量上限的1.2倍模型几秒钟就解完了数值稳定性也恢复正常。原因是求解器处理线性规划时M值过大会导致约束的数值尺度差异太大预求解器的缩放机制容易失效进而引发舍入误差。规范做法是对所有含M的约束M取该变量实际业务边界的1.2到2倍。比如容量上限120MWM取150以内足够。这条经验同样适用于其他0-1变量相关的“big-M”约束。5.3 线性化精度校验这步不能省分段线性化完成之后如果你直接拿结果写论文风险很大。我第一版代码就没做校验后来手工核对发现某条管道分段线性化后的流量与真实Weymouth方程计算值偏差达到6%。虽然选址结论没变但敏感性分析曲线的形状完全错了。从那以后我固定会在代码里加一个校验子函数把规划得到的P2G产气和气负荷回代到真实Weymouth方程里对比线性化流量与真实流量的相对偏差要求控制在2%以内。如果超过就增加分段数或调整分段区间。这步校验的成本很低但能让结果的可信度上一个台阶。写论文时还可以把这部分作为“模型精度验证”小节直接加分。5.4 环境配置与软件版本相关坑最后说几个环境层面的坑。第一YALMIP和Gurobi的版本匹配问题。新版Gurobi发布后旧版YALMIP有时无法正确识别求解器调用会出幺蛾子。建议先跑一个小的随机线性规划确认求解器链路通畅再跑大模型。第二MATLAB版本太旧也可能导致YALMIP某些函数行为异常。我遇到过一次莫名其妙的check报错升级MATLAB版本后自动消失。第三如果学校没有Gurobi用Matlab自带的intlinprog也不是不行。但要注意它的性能极限本文中等算例大约600多个0-1变量用intlinprog能跑但明显慢如果算例规模更大建议还是装个商业求解器。学术许可证对学生免费申请流程也不复杂别在这上面耽误时间。关于常见问题的汇总整理成速查表现象可能原因处理办法模型无解约束维度错、耦合符号反check(C)定位残差二分排查约束求解特别慢大M值过大、分段太密M取边界1.2倍分段控制在6-8段结果物理上不合理候选集设置问题、P2G符号反校验候选节点与平衡方程符号线性化误差偏大分段数不足或区间设置不当增加分段数并回代校验Weymouth求解器报数值警告大M值过大、参数尺度差异大缩小M值、检查变量量纲写在最后复现论文这件事10%的精力花在研究模型公式上90%都在“伺候”数据和代码。模型跑不通别慌把问题拆小、逐项验证比到处找别人的代码模板有用得多。尤其是气网的Weymouth线性化、P2G的耦合符号、M值选择这三个点这个模型里数值层面的坑基本都集中在这几个位置。把这几关过了整个代码框架就变成你自己的了。后面你想加储能、加氢能、加碳交易都能在这个底座上快速扩展这也是复现类工作最大的价值——不是拿到一份能跑的代码而是真正理解这个系统是怎么被数学表达出来的。