论文复现这件事十个有九个是倒在“理论上很简单”这六个字上。这次要说的这篇《计及碳排放成本的电-气-热综合能源系统节点能价计算方法研究》标题四平八稳但等你真动手做复现就会知道什么叫“真正做到了电热气潮流耦合”的分量——它不是把电网、气网、热网三个模型拼在一篇论文里而是让它们在同一套方程里互相咬合、彼此约束。我前后花了两周多时间从读模型到推公式再到把三个网络的潮流联立起来跑通最后成功提取出每个节点的能价信号。这篇就把整个复现过程和盘托出从建模思路、工具选型、关键代码实现到热网收敛困难这种坑位都给你捋一遍。无论你是做综合能源系统优化的研究生还是搞碳排放成本核算的工程师里面都有可以直接抄作业的东西。1. 项目概述这篇论文到底在算什么1.1 从标题里读出来的三个关键信号先拆标题。“计及碳排放成本”说明目标函数里不只有购能成本和设备运行成本还得把碳排放的代价装进去这直接决定了节点能价的量级和分布。中间的“电-气-热综合能源系统”是对象电网、气网、热网三条网络同时出现不是三个独立系统算完再叠个报告而是它们通过耦合设备热电联产机组、燃气锅炉、电锅炉、电转气装置发生能量交换一方的潮流改变会让另外两方的运行点跟着变。最关键的收尾词是“节点能价计算方法”。很多人会把节点能价理解成“平均成本除以电量”这是最大的误区。这里的节点能价对应的是电力系统中“节点边际电价LMP”的概念推广到多能网络——它是在满足所有物理约束的前提下某个节点新增一单位能量需求时系统总成本含碳排放成本的边际增量数学上就是对节点能量平衡约束的拉格朗日乘子。换句话说能价是影子价格反映了这个节点能量的真实稀缺程度。“真正做到了电热气潮流耦合”这句话看似感性其实是在说三个网络的潮流方程是同时联立求解的电力网络用交流潮流方程天然气网络用带Weymouth方程的稳态气流模型热力网络用水力-热力耦合模型耦合设备卡在中间做变量衔接。这样求出来的能价才是自洽的而不是先算完电网电价再回头去修气价。1.2 整体建模框架三层网络加耦合设备加能价提取整个复现我把它拆成四块理解。第一块是电力网络模型。节点注入有功、无功平衡支路潮流用极坐标下的交流潮流方程描述这部分有成熟的MATPOWER体系可以参考。第二块是天然气网络模型用稳态天然气潮流描述管道流量和两端气压满足Weymouth方程气源、加压站、天然气负荷都有对应的变量和边界条件。第三块是热力网络模型它和电网气网最大的不同在于“工质”会流动你得同时描述水力方程供回水管道流量和热力方程节点温度混合与管道热损失。第四块才是重点——耦合设备把三条网络串起来热电联产机组CHP消耗天然气、输出电和热燃气锅炉烧气产热电锅炉用电产热电转气P2G把多余电力转成天然气。每个耦合设备都是三网之间能量转换的事实通道也是能价传导的物理基础。在这个框架上优化模型是在满足电网、气网、热网所有潮流约束、设备出力约束、碳排放约束的前提下最小化系统总运行成本加碳排放成本。求解出最优运行点后再通过KKT条件提取各节点能量平衡约束的对偶乘子就是我需要的节点能价。整体流程用一句话概括就是先建网络模型再写耦合约束然后求解最优能量流最后从乘子里“解读”出价格信号。2. 复现准备工作工具选型与模型参数梳理2.1 工具选型为什么我选MATLAB加YALMIP复现这种三层网络耦合优化工具选错后面全是罪。我试过Python的Pyomo框架建模灵活但是处理天然气网Weymouth方程这种强非线性约束时调试起来不如MATLAB直观。PyPSA是电力系统专用工具箱气网和热网得完全自己搭等于半套重写。GAMS也有人用但报错信息对新手太不友好加上授权问题劝退。最终我定的是MATLAB加YALMIP加IPOPT求解器。选这个组合的理由有三个。第一YALMIP的约束建模方式非常接近数学公式我能把自己纸面上推的每个方程直接对应成一行代码排查模型错误时可以逐条核对心里有底。第二IPOPT是内点法求解器处理非线性规划NLP的收敛性和鲁棒性都够用这对Weymouth方程和热网温度混合方程这种强非凸约束来说很重要。第三MATLAB处理矩阵运算天然顺畅能在求解前快速检查雅可比矩阵和约束梯度的结构。事后证明这个选择很值。YALMIP内置了dual命令直接返回约束对应的拉格朗日乘子省去了我手写KKT系统再自己解乘子的工作量。要是自编求解器硬抠乘子光是乘子符号方向和缩放问题就能卡你一周。2.2 测试系统设计与数据准备复现论文时最尴尬的事情是论文里没公开完整的节点数据。参数只能从图和表格里反推反推不出来的就得按工程经验合理补全。我先参考经典的多能系统测试算例自己拼了一个中等规模系统电网取5节点环形结构网架参数参考IEEE 5节点系统修改气网用6节点含2个气源和1台加压站热网设计成6节点环状加枝状混合结构节点热负荷按居民区和工业区两类设置。这套系统不算大但足以验证“耦合效应能传导多远”。我自己设计数据时严守三条原则各网络基准值统一。功率基准取100 MVA气网基准值换算成标准状态下的体积流量热力网络功率基准也折算到同一套单位体系。耦合设备的容量和所在节点要合理。CHP放在电网节点2、气网节点3、热网节点4的交叉位置数据上保证它满发时电能和热能都有消纳空间不会让某个网络直接越限。碳排放参数要可调。排放因子和碳价先给一个基准值再写一个系数方便后面对比“碳价高低对节点能价的影响”。另外我把所有单位统一成国际制再进入计算。电网功率用MW、气网流量用m³/h标准状态、热网功率用MWth、价格用元/MWh碳排放用tCO₂/MWh。单位乱是这类复现最容易阴沟翻船的地方后面会专门讲。3. 核心实现三网耦合潮流的建模与能价计算全流程3.1 电力网络潮流建模电网部分我用了标准极坐标交流潮流。每个节点有两个变量电压幅值V和相角θ。节点注入有功功率P_i和无功功率Q_i满足P_i V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij)Q_i V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij)这套方程在综合能源系统里依然适用但要注意一点电网节点除了原有负荷和电源还多了两类耦合设备注入——CHP的发电功率、电锅炉/P2G的用电功率。所以在节点功率平衡方程里耦合设备的功率是作为变量出现的它们会随气网和热网的状态改变。电网不再是“给固定注入算潮流分布”的被动网络而是整个耦合系统里的一个环节。我建议复现时电网部分先用MATPOWER做个单网验证确认网架数据和潮流求解没问题再接入YALMIP。直接一上来就写大优化模型出了问题根本分不清是电网方程写错了还是气网约束写错了。这种“分层验证”的思路整个复现过程帮我省了一半的排查时间。3.2 天然气网络建模Weymouth方程与节点气压约束天然气网络的稳态潮流比电网复杂因为管道流量和气压是非线性关系。标准做法是用Weymouth方程描述管道流量f_ij s_ij C_ij √(s_ij (p_i² - p_j²))其中s_ij是流向标识正数表示从i流向jC_ij是管道常数。方程里出现的是气压平方差不是气压差这个非线性会让求解器非常吃初始点。管道两头气压要是初值给得离可行域太远IPOPT直接报“Restoration failed”。天然气管网还有两个特有的设备要建模。一个是加压站它靠压缩机消耗一部分天然气来提升下游气压消耗量跟压缩比和流量相关。另一个是气源气源有最大产气量约束和单位产气成本。天然气负荷的构成也得注意除了传统的燃气用户CHP和燃气锅炉的耗气量本身就是变量这部分把气网和电网热网绑在了同一套优化里。Weymouth方程还有个坑它经常写成Φ_ij C_ij sign(Δp) Δp²这种带绝对值和符号函数的形式符号函数在0附近不可导。我的处理办法是把流向固定为正值方向通过增加一个二进制变量做流向上限约束避开不可导点。这在YALMIP里可以借助binvar实现但会增加求解难度。如果系统规模不大也可以用“假设初始流向已知然后迭代修正”的方法先用初始流向求最优解再看是否出现逆流有逆流就更新流向再求一次。实测迭代两三次基本稳定。3.3 热力网络建模水力-热力模型热力网络是三层网络里最容易把人绕晕的部分因为它和电网气网有个本质区别能量是靠热水流动携带的你不能只定义“节点功率”还得描述配水管网的水流量和温度分布。水力模型部分节点有流量连续约束流入某节点的流量等于流出加节点负荷消耗。管路压降用Darcy-Weisbach方程描述流量越大压降越猛。热力模型部分每个管道有供水温度T_s和回水温度T_r节点热功率由流量和温差共同决定H_link c_p m_link (T_s - T_r)热水在管道流动时有热损失温度沿程衰减我按指数衰减模型处理衰减系数取0.01~0.05 /km具体看保温状况。而节点处多股水流汇合时混水温度按流量加权平均计算。这套模型的难点在于它有一堆非线性项——流量乘温差、温度衰减指数函数、混水温度加权。热网方程和非线性电网方程、气网方程叠加在一起整条约束雅可比矩阵的结构变得非常稠密且病态这也为后面的收敛问题埋了伏笔。复现热网时有个通用技巧如果热网节点数不多先把回水温度当作固定值求一版解再放开回水温度做精修。3.4 耦合设备建模CHP、燃气锅炉、P2G怎么接入耦合设备是“电热气潮流耦合”的物理核心建模时最忌单打独斗。我分别处理了四种设备。CHP热电联产机组是其中最重要的。它的输入是天然气输出是电和热三个量之间有耦合关系。最简单的模型是定电热比P_e η_e F_gasP_h η_h F_gas。但实际CHP的电出力可调范围受热出力约束我再加了一个可行域包络约束保证电热出力组合在三角或多边形可行域内。CHP运行时还有一个特性为了满足热负荷它可能被迫发出超过单纯经济调度所需电量的电这个特性会让电网节点的能价发生变化是后面结果分析的重点。燃气锅炉模型简单烧气产热效率常数P_h η_b F_gas。电锅炉是耗电产热P_h η_eb P_e。P2G是电转气消耗电力和水产出天然气主要是氢气/甲烷F_gas η_p2g P_e这个方向通常只在风电出力大、气价贵的场景下才被启用。这些设备接入时还有一个关键细节耦合设备连接不同网络节点时所在节点要明确。CHP在电网侧连节点2、在气网侧连节点3、在热网侧连节点4这意味着CHP的发电出力影响电网节点2的功率平衡耗气量是气网节点3的负荷产热量是热网节点4的注入。能价传导的路径正是沿着这些“端口节点”铺开的。3.5 碳排放成本建模与节点能价提取碳排放成本建模可以走两条路一是硬约束路线给系统设定总碳排放上限超过上限就不可行碳价是约束的对偶乘子二是成本化路线给每吨碳排放一个价格直接加进目标函数。论文标题写的是“计及碳排放成本”所以我采用成本化路线为主、排放上限约束为辅的组合方式。目标函数变成min ∑ 购能成本 ∑ 设备运行成本 ∑ 碳排放量 × 碳价其中碳排放量包括外购电力折算的上游排放、天然气燃烧的排放CHP和燃气锅炉、P2G如果消耗的是高碳电力也有间接排放。这部分要写得足够细体现在结果上就是碳价一涨气源节点的能价和CHP所在节点的能价会联动上涨但上涨幅度会因为节点在气网中的位置不同而不同——离气源远、管压低的节点气价升得更狠。求解这一个NLP模型后节点能价的提取方式就非常关键了。YALMIP里求解模型会返回一个sol结构对每个约束调用dual函数就能拿到对偶乘子。具体来说电网节点有功平衡约束的乘子对应节点电价气网节点流量平衡约束的乘子对应节点气价热网节点热功率平衡约束的乘子对应节点热价。提取后要做一个符号和量纲校验乘子单位应该是元/MWh如果出现量纲异常九成是目标函数里成本单位和功率基准没对齐。我个人建议提取之后顺便算一下每个节点的系统边际成本做对比如果能价和边际成本偏差太大说明乘子取错了约束或符号反了。% 伪代码示意YALMIP求解与乘子提取 % 定义变量略 Constraints [ ... ]; % 所有网络约束 耦合设备约束 Objective ...; % 购能成本 运行成本 碳成本 ops sdpsettings(solver,ipopt,usefull,1); sol optimize(Constraints, Objective, ops); % 提取能价节点功率平衡约束的拉格朗日乘子 bus_price dual(Constraints(1)); % 电网节点约束 gas_price dual(Constraints(2)); % 气网节点约束 heat_price dual(Constraints(3)); % 热网节点约束4. 实操过程一个最小三网耦合系统跑通全流程4.1 最小案例场景设计为了人肉验证模型正确性我设计了一个最小系统电网3节点、气网3节点、热网3节点三条网络各有1个源、1个负荷、1段传输通道中间用1台CHP把三段能量串起来。主气源在气网节点1CHP在气网节点2取气电力网络里CHP接在电网节点2热力网络的CHP接在热网节点2另外热网节点3设了一个燃气锅炉做备用热源。这个系统小到什么程度小到我可以手算平衡给定负荷、CHP电热比和锅炉效率先手工求解一个粗略解再用这个解当初值喂给IPOPT。这个手动初值对非线性规划收敛帮助巨大——热网温度方程和Weymouth方程如果初值远离可行域IPOPT经常在没有可行域的空间里挣扎最后给你一个“Converge to a point of local infeasibility”。我还额外做了四个版本的对照模型纯电网版无气网热网CHP当成固定电源、电气耦合版没有热网CHP只发电不用热、电热气全耦合版三种网络完整接入不考虑碳成本、电热气耦合加碳成本版。后面对比这四个版本的结果能清楚看出每一层耦合和碳成本对节点能价的影响。4.2 结果解读耦合效应怎么改变节点能价先看第一个对照结果。纯电网版里电网节点2接入CHP作为固定电源节点能价主要由外购电边际成本和网损决定大概在320~340元/MWh。电气耦合版里CHP的电出力不再固定而是由气网气价决定结果节点2能价掉了近10%因为天然气价折算成单位电能成本比外购电便宜。这说明一个道理多能系统的节点电价不单由电网自己决定气源价格、管道压力约束都会通过CHP的“耗气边界”传导进来。电热气全耦合时变化更明显。由于CHP必须满足热网节点2的热负荷它的电出力被“捆绑”调度某些时段电出力被迫上升电网节点2的能价反而因为本地供给增加而下降。而热网节点2的热价因为CHP承担了大部分热负荷边际热源成了备用燃气锅炉所以热价主要被锅炉的天然气效率和气价锚定。这个结果正好验证了耦合效应的传导逻辑热负荷通过CHP“绑架”了发电出力进而压制了局部节点电价同时“抬高”了气价——因为CHP和锅炉都在抢气。再看加了碳成本的版本。碳价按100元/tCO₂计CHP和燃气锅炉的单位碳排放量折算进目标函数后气网节点2和热网节点3的能价同步上浮。这里最有意思的是热价上涨幅度大于气价上涨幅度——因为燃气锅炉效率比CHP低单位热输出的碳排放更多碳成本分摊到每MWh热量的涨幅自然更大。用通俗的话说碳成本不是均匀摊到大伙头上的谁的效率低、谁的排放因子高谁所在的就能价承担更多碳成本。这恰恰是计及碳排放成本的节点能价计算方法的独特价值也是复现这篇论文最值得讲清楚的地方。4.3 收敛性与求解参数设置三网耦合的NLP求解收敛性是复现过程最大的硬骨头。我试过直接上全耦合模型IPOPT默认参数跑出来的结果十次有七八次不收敛。后来摸索出一套可行的参数和策略。IPOPT的tol设1e-6每轮迭代步长限制可以放宽但线搜索的接受阈值不要放松否则容易在热网温度方程上胡乱跳。变量缩放很重要。电网功率变量动辄上百MW热网流量变量可能是几kg/s管道气压是几MPa混在一起求解器数值条件极差。我统一做了标幺化功率用100 MVA基准流量用基准流量气压用基准气压把所有变量都压到1e-2到1e2的量级范围。这一步做完收敛率直接从不到三成提到了七成以上。耦合设备的变量初始点注意不要全部设成0。CHP电出力初值不要设0给一个接近当前热负荷对应出力的值热网温度初值在85℃附近气网气压初值在3~5 MPa区间。这些初值不要求准但要求“位置对”让求解器从一开始就在物理可行域附近搜索。我强烈建议复现这种论文时把收敛问题拆成三个层次解决先解决变量缩放再解决初始点最后才调节求解器参数。顺序反了你会陷入一个调参的泥潭——今天调好了明天换个负荷又炸。5. 踩坑实录与问题排查速查表5.1 最大的坑耦合变量单位不统一导致雅可比矩阵病态这个坑我栽得最惨。第一次跑全耦合模型IPOPT直接报“The maximal condition number is too large matrix is singular”一度以为方程写错了。后来定位到是单位基准问题电功率变量是MW气网流量是kg/s热网功率是MWth但Weymouth方程里的管道常数C_ij用的又是m³/h和bar这些量级差了好几个数量级导致雅可比矩阵条件数爆炸。修法不复杂所有模型输入先做标幺化气网流量、气压、管道常数全部折算到统一基准热网流量同理。我建议你写代码之前先把所有物理量归一化到一套标幺体系里再写约束方程。别偷懒这个步骡省不了。5.2 热网温度方程带来的强非凸振荡热网节点水温和回水温度混合方程是强非凸的特别是多股不同温度的水流在节点混合时“温度乘流量”耦合项很容易让迭代过程在两个局部解之间来回振荡表现为IPOPT迭代数到几百步目标函数还在锯齿状跳动。我的经验是先用固定回水温度的简化模型求一个可行解拿到这个解之后把其中的回水温度作为初始化值再放开回水温度变量重新求解。这个方法本质上是一种同伦/延拓思想用一个更简单的模型引导求解器进入凸性较好的区域再切入完整模型。放在工程里说就是“先用傻瓜模型把车开上正路再切到专业驾驶模式”。5.3 KKT乘子提取不出来或数值异常YALMIP的dual命令提取乘子但很多人不知道乘子只在约束边界处有意义如果节点约束不起作用没有到达平衡边界乘子会是0这不是bug是数学性质。所以提取能价前先确认你关注的那个节点约束是否是起作用约束——检查该节点是否有注入变化引起目标函数变化。如果节点能价提取到0或异常大多半原因如下约束符号反了乘子出现镜像但数值没问题校正方向即可。目标函数单位是万元功率基准是MW乘子单位变成万元/MWh要做单位换算。求解器返回的不是全局最优而是局部最优乘子值也会跟着偏建议换多个初始点多试几次。5.4 问题速查表现象可能原因解决办法IPOPT报矩阵奇异变量量级悬殊单位未标幺全网变量标幺化统一基准迭代振荡不收敛热网温度非凸耦合初始点偏离物理可行域先固定回水温度求解再放开复算节点能价提取到0该节点约束不起作用或乘子提取了非平衡约束确认节点约束处于活动边界检查约束索引能价出现负值目标函数方向或符号设置错误检查优化方向复核成本项符号气网流量莫名逆流未固定流向符号函数不可导导致求解器乱跳用二进制变量约束流向或迭代修正流向碳成本对能价无影响碳成本项没写进目标函数或排放因子为0检查目标函数是否包含碳项核对排放因子6. 复现过程的经验总结与下一步扩展方向最后分享几条从这次复现里沉淀下来的体会。做论文复现尤其是这种多网耦合类的最忌一上来就堆代码。我的习惯是先在纸上把每个网络的变量、约束、耦合关系画清楚——画成一张类似“能流图约束表”的东西标注哪些变量是跨网络共享的、哪些约束是连接两个网络的。代码只是这张图的“翻译”而已。这次CHP设备的建模如果当初没有先在纸上画清楚“电网节点2、气网节点3、热网节点4之间通过CHP形成能量闭环”的关系后面调试时十有八九要找半天bug。还有一点要强调复现论文不是目的理解算法背后的经济学含义才是目的。节点能价这个东西换个角度看就是多能系统的“路标”——它告诉你哪个节点能量稀缺稀缺到什么程度碳成本又让哪些节点涨价。把这些机制吃透了你才有能力去做扩展。我个人计划把这个模型往三个方向扩展一是从稳态走向多时段动态考虑热网管道的蓄热特性这样热价会出现时间平移效应二是加入阶梯碳价机制让碳成本对能价的影响是非线性的三是把求解规模放大到几十个节点的算例系统测试这个框架的工程实用边界。如果这些方向跑通了再来分享第二篇。