做区域综合能源系统课题的朋友应该都有这种体会单独算电网潮流、单独算天然气网络、单独算热力管网各自都有成熟工具和现成程序但一旦系统里加了热电联产机组、电锅炉、燃气锅炉这些“跨界设备”三个网络就再也拆不开了。我在这类项目里反复折腾过很久最直观的感受是多能流计算的难点不在于某一个网络本身而在于耦合元件怎么建模、三个网络的方程怎么组织在一起迭代收敛。这篇文章就把这块内容完整讲透从数学模型到Matlab代码架构再到算例验证和调试经验做一个偏实战向的拆解。这个内容适合正在做综合能源系统潮流计算、规划优化、运行调度方向的研究生和工程师尤其是手头需要自己编写代码实现电气热多能流计算的人。文章核心聚焦“计及多能耦合”这个关键点不会刻意去讲Matlab基础语法而是把精力放在“怎么把物理模型翻译成可运行的代码”这件事上。1. 多能流计算到底在解决什么问题1.1 三个网络不再是“各算各的”区域综合能源系统里电网、天然气网、热网原本是三个独立运行的能源输配网络。电网用交流潮流描述天然气网用节点压力方程和管道流量方程描述热网用水力-热力联合方程描述。单独看任意一个网络求解方法都很成熟电网有牛顿-拉夫逊法气网有节点法热网有水力/热力迭代法。问题出在“多能耦合”上。一个典型的园区级综合能源系统往往包含燃气轮机组加余热锅炉组成的CHP机组、电锅炉、燃气锅炉、P2G装置等设备。燃气轮机组消耗天然气发出电力的同时余热锅炉又把烟气余热转化为热网的热源电锅炉从电网取电却向热网供热燃气锅炉消耗天然气供热。这样一来电网的电力潮流会改变电锅炉的取电功率进而影响热网的注入热量热网的负荷变化又会影响CHP机组的热出力指令进而改变燃气轮机组的发电功率和天然气消耗量P2G的制气量又会影响气网的节点流量平衡。三个网络被这些耦合设备捆绑成了一个真正的“联立方程组”。我记得第一次自己搭这个模型的时候想省事先把电网潮流跑收敛然后把CHP电出力代到热网里再把热网热负荷代回电网结果循环迭代了三十几次都不收敛。后来才想明白这不是程序写法问题而是物理上电网、气网、热网之间存在双向反馈必须把耦合变量纳入统一的迭代框架或者把所有方程联立起来一起解。这个理解是这个项目里所有代码设计的前提。1.2 统一求解和分解求解的取舍实现多能流计算主要两条路线统一求解法和分解迭代法。统一求解法把所有子系统的方程组装成一个大型雅可比矩阵状态变量同时包含电网节点电压幅值、电压相角、气网各节点压力、热网各节点温度和流量一次牛顿迭代同时修正所有变量。这个方法的优点是收敛性通常更好不存在“内外层迭代”之间的误差传递缺点是雅可比矩阵规模大、稀疏结构复杂编程实现难度高尤其耦合元件对雅可比矩阵的贡献项很容易写错。分解迭代法是把电网、气网、热网分别用各自的求解器计算耦合元件在每个子系统里都表现为边界条件。比如电锅炉在电网侧是有功负荷在热网侧是热源二者通过 η_eb × P_eb Q_eb 这个关系式在外部迭代中建立联系。每轮迭代先算电网更新耦合设备功率再算气网再算热网检查耦合变量前后变化量是否小于阈值。这种方法模块化强调试方便可以复用现成的单网络潮流程序缺点是收敛速度可能慢某些情况下会振荡甚至发散。我的建议是初期代码实现优先走分解迭代法。原因是你能清楚地看到每一个物理量的变化过程出了问题容易定位是哪个网络、哪个耦合环节导致的。等把物理机制摸熟了想追求计算效率或者发高质量论文需要做大规模系统仿真时再升级到统一求解法。两条路线的核心区别可以用下面这个表来理解对比维度统一求解法分解迭代法状态变量维度所有网络变量同时求解维度大每个网络分别求解维度小雅可比矩阵需要组装跨网络耦合块各网络雅可比独立编程难度高耦合项易出错低模块独立性强收敛性较强但依赖好的初值可能慢或振荡调试友好度不友好出错难定位直观可逐步验证适用场景大规模系统、高耦合度、研究型工程方案快速评估、学习框架1.3 这套代码能做什么、不能做什么明确边界很重要。基于稳态模型的多能流计算解决的是“给定各网络的负荷和耦合设备运行状态求解系统稳态运行点”的问题。它能够用来做区域综合能源系统的规划方案比选、运行工况分析、耦合设备容量灵敏度研究也能给后续的优化调度模型提供初值或可行性校验。它不能做的是动态过程分析。电网的频率动态、天然气管道的气体波动传播、热网的慢热惯性过程都需要微分方程或偏微分方程描述不在稳态多能流的范围内。如果课题需要做短期动态仿真那应该基于时域仿真工具或者求解DAE方程组和这里讨论的稳态潮流是两套东西。2. 电气热子系统的模型怎么建才不踩坑2.1 电网部分极坐标NR法加功率不平衡量电网模型不必多解释重点讲代码实现里容易出问题的两个地方。第一采用极坐标牛顿-拉夫逊法时状态变量是节点电压幅值U_i和相角θ_i。对PQ节点同时列写有功和无功不平衡方程对PV节点只列有功方程松弛节点作为参考节点不参与迭代。功率不平衡方程为ΔP_i P_i_spec - U_i ∑_j U_j (G_ij cosθ_ij B_ij sinθ_ij)ΔQ_i Q_i_spec - U_i ∑_j U_j (G_ij sinθ_ij - B_ij cosθ_ij)角度差θ_ij θ_i - θ_j。Matlab代码实现时角度单位必须统一用弧度。有个很隐蔽的坑很多初学者在构造节点导纳矩阵的时候用复数矩阵直接计算没问题但在计算不平衡量时如果角度用了度整个雅可比矩阵数值会差57倍左右程序表现就是迷宫一样的乱跳。第二耦合设备参与电网潮流时它们在电网侧的表现要区分类型。CHP机组如果是定功率控制在电网里就是给定有功和无功的传统发电机组节点如果CHP参与电压调节就要按PV节点处理。电锅炉、P2G这类设备本质上是可控负荷电网建模里作为恒功率负荷节点处理即可。核心思路是耦合设备对电网的注入功率必须先更新再进入电网牛顿迭代不能滞后一个迭代周期。电网求解函数的主体结构大致是function [U, theta, iter] NRPowerFlow(bus, branch, gen, load) % bus: 节点参数结构体数组 % branch: 支路参数结构体数组 % 返回节点电压幅值U和相角theta Y buildYbus(bus, branch); % 构建节点导纳矩阵 % 初始化平启动 U ones(nBus, 1); theta zeros(nBus, 1); for iter 1:50 [dP, dQ] calcMismatch(U, theta, Y, bus, gen, load); if max(abs([dP; dQ])) 1e-6 break; end J buildJacobian(U, theta, Y); % 形成极坐标雅可比 dX J \ [dP; dQ]; % 稀疏LU分解求解修正量 % 更新状态变量 end end2.2 气网模型Weymouth方程与节点压力迭代天然气稳态管网模型比电网复杂一点核心方程是管道流量与两端压力的关系。工程中常用Weymouth方程描述高压输气管道中的稳态流量q_mn sign(p_m - p_n) × C_mn × sqrt(|p_m² - p_n²|)其中C_mn是管道常数由管径、长度、摩擦系数、气体组成决定。这个方程写代码时注意两点。第一取平方差要加绝对值符号由两端压力高低决定。第二压力必须用绝对压力不能用表压否则低压管网计算会出现开方负数的荒谬结果。气网的节点方程是流量平衡每个节点的流入流量等于流出流量加负荷消耗量。已知压力和未知压力节点分开处理通常气源节点压力已知负荷节点压力未知。未知压力节点用牛顿法求解修正方程的核心雅可比元素来自Weymouth流量对端点压力的偏导数∂q_mn / ∂p_m C_mn × p_m / sqrt(|p_m² - p_n²|)这个导数在压差趋近于零的时候会趋向无穷大这是气网迭代不收敛的常见原因之一。解决办法是给流量方程加一个很小的正则化项或者限制压力差的下限避免迭代过程中出现两个节点压力几乎相等的情况。管道流量的计算函数可以封装成这样function q gasPipeFlow(p_in, p_out, C) % 输入入口压力、出口压力、管道常数 % 输出标准体积流量正的表示从in流向out dp_sq p_in^2 - p_out^2; q sign(p_in - p_out) * C * sqrt(abs(dp_sq)); end2.3 热网模型水力与热力联合求解热网分成水力模型和热力模型两部分。水力模型描述管道流量分布节点满足质量流量守恒即流入节点的流量等于流出节点的流量加节点负荷消耗流量。环状管网还需要列写回路压降方程。实际项目中区域供热管网大多是环状结构但做多能流计算研究时往往先采用辐射状或树状拓扑这样可以显著降低水力求解难度。热力模型包含节点功率方程和管道温降方程。节点功率由供回水温差和质量流量决定Q_node c_p × G_node × (T_s - T_r)其中c_p是水的比热容G_node是流经负荷节点的质量流量T_s是供水温度T_r是回水温度。多个支路流入同一节点时混合温度按下式计算T_mix Σ(G_j × T_out_j) / ΣG_j这是热网模型里最容易忽略的一项。如果多条供热管道汇入同一个节点而你把每条管道的末端温度直接拿来当节点温度结果会偏差很大。管道温降描述采用苏霍夫公式T_end T_amb (T_start - T_amb) × exp(-λL / (c_p G))其中λ是管道单位长度热损系数L是管道长度T_amb是环境温度。这个公式说明流量越小温降越明显流量很小时管道末端温度会无限接近环境温度。数值实现上用这个公式直接计算末端温度即可但要注意G不能取零。热网求解的典型做法是“先水力,后热力再迭代”。原因是水力方程的解管道流量决定了热力方程中的流量参数而热力方程中的节点温度又会影响回水温度进而影响负荷侧的供回水温差计算。如果只做一次水力求解就固定流量再算热力对于不含回水温度迭代的简化系统是够用的但严格建模时建议做一轮外部迭代把供水温度和回水温度反复更新到收敛。3. 耦合元件建模多能流计算的灵魂环节3.1 CHP机组的两类运行模式CHP热电联产机组是综合能源系统里最重要的耦合设备。它消耗天然气同时输出电功率和热功率把电网和气网、热网三张网联系在一起。建模时先分清运行模式。背压式CHP机组的电热出力之比是固定的称为定热电比模式H_CHP C_m × E_CHP其中C_m是热电比常数。此时CHP的电力出力和热力出力完全耦合给定其中一个另一个就确定了。抽凝式CHP机组可以调整热电比运行范围是一个可行域需要满足电出力上下限、热出力上下限以及发电与供热之间满足燃料消耗约束。多能流计算中通常给定某个运行点然后校验是否落在可行域内。燃气消耗量是联系气网的纽带。CHP消耗的天然气流量F_CHP E_CHP / (η_e × LHV)其中η_e是CHP的发电效率LHV是天然气低位热值。热效率η_h和发电效率之和不等于总效率因为余热回收有损失通常总效率在0.75~0.9之间。在分解迭代法中CHP的作用流程是电网潮流算完后得到CHP的电出力E_CHP按热电比或设定运行模式更新热出力H_CHP然后在热网中把这个热出力作为对应节点的热源注入量同时按效率公式更新气网中该机组节点的天然气消耗流量。一步滞后带来的误差通过外部迭代消除。在统一求解法中CHP机组的建模要复杂一些电网方程、热网方程、气网方程之间的耦合项直接出现在统一的雅可比矩阵中。具体说热网中CHP节点的热源项对电网中E_CHP求偏导气网中CHP节点的气负荷项也对E_CHP求偏导这些偏导数都要作为非零元素填到雅可比矩阵的跨网络分块里。3.2 电锅炉、燃气锅炉、P2G三种耦合的简化处理电锅炉是把电转成热的设备Q_eb η_eb × P_eb在电网里它是一个有功负荷取电功率P_eb在热网里它是一个热源注入热功率η_eb × P_eb。这里效率η_eb通常在0.95~0.99之间电转热的损失远小于燃气锅炉。电锅炉的优势是响应速度快但在电网负荷高峰时段会加剧电力供需矛盾所以多能流计算中它的取电功率往往与电网运行状态互相制约。燃气锅炉只耦合气网和热网Q_gb η_gb × F_gb × LHV它从气网取气向热网供热建模最简单本质是两个网络之间单向耦合。燃气锅炉通常作为调峰热源热网热负荷超过CHP供热量时启用。P2G设备是实现电气双向耦合的典型装置消耗电能制取天然气氢气或甲烷。模型上简化为F_p2g η_p2g × P_p2g / LHV在电网里是电负荷在气网里是气源注入。P2G通常用于可再生能源消纳场景在风电、光伏出力大的时段把多余电力转化成天然气储存。多能流计算中P2G的功率是给定的运行指令不需要额外求解。这几种耦合设备的关系可以归成一张表设备输入侧输出侧效率典型值关键方程CHP机组天然气电热发电效率0.35~0.45热效率0.4~0.5H C_m×E, F E/(η_e×LHV)电锅炉电热0.95~0.99Q η_eb×P燃气锅炉天然气热0.85~0.93Q η_gb×F×LHVP2G设备电天然气0.50~0.70F η_p2g×P/LHV热泵电环境热热COP 3~5Q COP×P3.3 储能和冷负荷的处理方法完整的区域综合能源系统里通常还有蓄热罐、蓄电、吸收式制冷等设备。蓄热罐本质上是能量时移装置在稳态多能流计算中如果给定蓄热罐在某个计算时刻的充放热功率可以直接把它当作热网中一个负的或正的负荷节点处理并不复杂。冷负荷可以近似看成“负热负荷”或通过吸收式制冷机组转换。在热电联产园区里典型的做法是把冷负荷转化为热负荷需求因为吸收式制冷机以热源为输入Q_cold COP_c × Q_heat_in这样处理之后整个系统的热网方程结构不需要改变只是负荷节点多了一些等效热功率。如果课题明确要求电制冷和吸收式制冷并存那就需要在电网中增加电制冷机组的电负荷在热网中增加吸收式制冷的热负荷模型扩展方向是清晰的。4. Matlab代码实现架构与迭代流程4.1 数据结构和输入文件怎么组织写多能流代码最大的教训就是不要把所有变量都塞进脚本里硬编码。系统规模一扩大脚本方式会让你改一个管道参数找半天而且极易出错。我采用的方案是结构化定义三类网络和一个设备结构体。电网部分用bus和branch结构体数组每个节点包含编号、类型松弛节点/PV节点/PQ节点、有功负荷、无功负荷、电压幅值初值等字段每条支路包含首末端节点号、电阻、电抗、对地导纳。天然气网用gasNode和gasPipe结构体节点字段包括编号、类型气源/负荷、压力初值、节点流量注入管道字段包括两端节点编号、管道常数C。热网部分用heatNode和heatPipe结构体节点字段包括类型热源/热负荷/普通节点、供热温度、回水温度、质量流量注入管道字段包括两端节点号、长度、热损系数。耦合设备单独用unit结构体数组定义每个设备的字段包含设备类型、所在电网节点号、所在气网节点号、所在热网节点号、效率、热电比等。这样做的好处是主程序迭代时只需要遍历unit数组就能更新所有耦合边界完全不用改核心求解函数。4.2 主程序骨架分解迭代法的完整流程分解迭代法的核心逻辑是循环执行“电网求解→更新耦合设备热/气出力和消耗→热网求解→更新耦合设备电/气消耗→气网求解→检查耦合变量收敛”。下面给出一套可以直接参照的主循环框架function x solveMEF(d) % d: 系统数据结构包含电网、气网、热网、耦合设备定义 x.E initElectric(d.E); % 电网初值U1, theta0 x.G initGas(d.G); % 气网初值所有节点压力基准压力 x.H initHeat(d.H); % 热网初值供水温度T_s0, 回水温度T_r0, 流量初值 for k 1:100 x.E_old x.E; x.G_old x.G; x.H_old x.H; % ---- 步骤1求解电网潮流 ---- % 根据上一轮的气网/热网状态更新耦合设备在电网中的注入/消耗功率 [P_inj, Q_inj] updateElectricBoundary(x, d); x.E NRPowerFlow(d.E, P_inj, Q_inj, x.E); % ---- 步骤2更新热网边界并求解热网 ---- % CHP热出力、电锅炉热出力根据最新电网结果更新 [Q_heat_source, G_load] updateHeatBoundary(x, d); x.H HeatFlow(d.H, Q_heat_source, x.H); % ---- 步骤3更新气网边界并求解气网 ---- % 燃气锅炉、CHP的用气量以及P2G的产气量都要更新 [F_gas_load, F_gas_source] updateGasBoundary(x, d); x.G GasFlow(d.G, F_gas_load, F_gas_source, x.G); % ---- 步骤4收敛判据 ---- dE max(abs([x.E.U - x.E_old.U; x.E.theta - x.E_old.theta])); dG max(abs(x.G.p - x.G_old.p)); dH max(abs([x.H.T_s - x.H_old.T_s; x.H.T_r - x.H_old.T_r])); if max([dE, dG, dH]) 1e-6 fprintf(多能流计算收敛迭代次数%d\n, k); return; end end error(多能流计算不收敛请检查初值或系统参数); end这个主循环里有几个细节值得注意。迭代上限不要设太大100次以内足够判断收敛性如果100次还不收敛大概率是模型或初值的问题继续迭代也没有意义。收敛判据用状态变量的变化量而不是不平衡量更稳妥。耦合变量的更新用最新值即前一个网络刚算出来的结果立刻用于下一个网络的边界更新这种风格叫“Gauss-Seidel式”更新比全部用旧值更新的“Jacobi式”收敛更快。4.3 热网求解函数的关键实现热网求解函数内部要处理水力迭代和热力迭代的嵌套关系。一个简化的实现思路是先假设各节点供水和回水温度已知求解节点质量流量然后固定流量更新各管道末端温度和各节点混合温度检查温度变化量如果不收敛则用新温度再次求解水力方程。节点质量流量求解的核心是节点流量守恒方程写成矩阵形式就是A × G G_load其中A是节点-支路关联矩阵G是管道质量流量向量G_load是节点注入流量向量。对树状管网这个方程直接就解出来了环路管网需要额外列写回路压降方程联合求解。温度迭代部分的核心代码如下function T_node updateNodeTemp(G_matrix, T_in, T_amb, pipes, cp) % G_matrix: 各管道质量流量向量 % T_in: 管道入口温度向量 % 先用苏霍夫公式计算各管道末端温度 for i 1:length(pipes) T_out(i) T_amb(i) (T_in(i) - T_amb(i)) * ... exp(-pipes(i).lambda * pipes(i).L / (cp * G_matrix(i))); end % 再按混合定律求各节点的混合温度 for k 1:nNode inflow find(tailNode k); % 流入节点k的管道集合 numerator sum(G_matrix(inflow) .* T_out(inflow)); denominator sum(G_matrix(inflow)); if denominator 0 T_node(k) numerator / denominator; end end end这里有个经验性的处理分母流量很小时混合温度数值会不稳建议对denominator加一个很小的保护阈值或者直接将流量小于某个下限的管道从混合计算中剔除否则温度会出现不合理的跳变。5. 算例设计与结果验证5.1 测试系统的搭建思路为了验证代码的正确性需要一个不复杂但包含典型耦合元素的测试系统。我用的是下面这种结构电网采用5节点系统其中节点1是松弛节点外电网等值节点2接入一台CHP机组节点3接入电锅炉节点4和节点5是常规电负荷。天然气网采用4节点结构节点1是气源压力给定节点2连接CHP机组气负荷节点3连接燃气锅炉气负荷节点4连接P2G设备气源。热网采用6节点结构节点1是CHP热源节点节点2是电锅炉热源节点节点3是燃气锅炉热源节点节点4和节点5是热负荷节点节点6是回水汇集节点所有管道采用辐射状连接。测试系统的运行场景设定为冬季典型工况CHP机组带基础电负荷和热负荷电锅炉作为补充热源燃气锅炉调峰。P2G设备以小功率运行。这样能保证三个网络之间至少有四个耦合元件在同时工作耦合关系够复杂。仿真参数的选取也要符合实际。电网基准容量取100 MVA电压基准取10 kV。气网基准压力取1 MPa天然气低位热值取35.8 MJ/Nm³。热网供回水温度基准取130℃/70℃工程中高温水系统常见参数水的比热容取4.2 kJ/(kg·K)。所有效率和热电比参数设置要保证系统能量平衡有意义不能让某个设备效率大于100%。5.2 收敛过程和结果评估实际运行这套代码时典型工况下分解迭代法大约需要15到30次外部迭代收敛到1e-6的精度。若从平启动初值开始前五轮迭代耦合变量变化量较大之后快速下降基本呈现线性收敛趋势这也符合Gauss-Seidel型迭代的特征。收敛后可以校核几个总量指标全系统电功率平衡即发电机出力与外电网交换功率之和等于负荷加电锅炉和P2G耗电热网供热量等于热负荷加管道散热损失气源总供气量等于CHP、燃气锅炉耗气量减去P2G产气量。能量平衡校验中每一项的误差在1e-6量级说明方程组装和迭代过程是正确的。我还建议做一步“解耦一致性验证”将电锅炉、燃气锅炉、P2G的运行功率全部设为0热网只有CHP供热这时把CHP的电出力和热出力固定将CHP节点当作电网的PQ节点气网节点当作普通负荷此时多能流程序的电网部分结果应该和单独运行电网潮流程序完全一致。这一步验证过关代码的电网模块基本可信。同理也可以分别验证气网、热网模块。5.3 耦合强度对收敛性的影响在代码跑通之后可以进一步做耦合强度的灵敏度分析这也是论文里常用到的手段。方法是固定其他参数逐步改变CHP热电比或者电锅炉容量观察外部迭代次数的变化。一般情况下热电比越大CHP电出力到热出力的传递越强电网和热网之间的耦合越紧密分解迭代法需要的迭代次数会上升。如果热电比超过某个阈值甚至可能出现振荡不收敛的情况。我实测遇到过一个案例把电锅炉功率从1 MW逐步增加到8 MW时迭代次数先是缓慢增加超过6 MW后突然发散。原因是电锅炉功率增大导致电网节点电压偏低电网潮流求解精度下降反过来又影响热网边界形成正反馈。解决办法是调整电网无功补偿容量或者将电锅炉连入的节点改为更靠近电源的节点。这个现象本身也是一个有价值的研究点很多文献专门讨论过耦合强度与收敛性的关系。6. 常见问题与排查技巧6.1 不收敛先看初值再看雅可比最后看单位多能流计算不收敛的排查顺序很重要我自己的经验是按“初值→雅可比→单位”的顺序来切勿一开始就怀疑算法框架。初值问题是新手最容易犯的错。气网压力初值如果设得过高或过低Weymouth方程的流量初值会离谱导致气网牛顿迭代第一步就失败。一个稳妥的做法是把所有未知压力节点初值设为气源压力的0.95倍左右不要平启动设成1 MPa以下一大堆低压节点。热网温度初值也不要取环境温度供水温度初值可以直接用设计供温130℃回水温度用70℃这样热力迭代的温度差不会太大。雅可比矩阵的问题在统一求解法中比较常见。排查手段是有限差分验证将状态变量整体加一个1e-6的小扰动重新计算不平衡量与解析雅可比矩阵对比。如果某个非零元素对不上定位到具体方程后重点检查对应耦合元件的偏导数推导。单位问题在气网和热网中尤其明显。天然气流量可以用标准体积流量Nm³/h也可以用质量流量kg/s两者之间要乘标准密度0.717 kg/Nm³左右。功率和流量换算时还要除以3600转换秒和小时。如果不统一单位能量平衡校验一定对不上。6.2 一个实际排查案例气网压力振荡有一次跑算例气网迭代在第四轮之后开始振荡压力值在1.0 MPa和0.92 MPa之间来回跳动。排查发现是P2G设备在电网侧取电功率过大导致电网潮流计算出的P2G电耗偏高按效率换算成产气流量后这个流量加到气网节点上气源节点压力被推高下一轮气网求解又把该节点负荷调低造成振荡。处理办法是在外部迭代中引入阻尼系数即本轮使用的气负荷值取上一轮值和本轮更新值的加权平均F_new α × F_calculated (1-α) × F_oldα取0.6左右阻尼效果明显迭代很快恢复收敛。这种方法虽然只是工程技巧但在很多耦合度较高的系统里非常管用。6.3 从“能跑”到“可信”的三个验证手段程序写完能出结果并不代表可信。要证明代码正确至少要过三关。第一关是有限差分验证雅可比矩阵这保证方程组的求解方向是正确的。第二关是能量平衡校验保证所有功率、流量在物理上是守恒的。第三关是对比验证找一个同类文献的算例系统用你的代码重跑一遍看结果能否和文献给出的电压、压力、温度分布对上。还有一个非常实用的小技巧把热网的管道热损系数暂时设成0环境温度设成不影响温降的数值此时热网供水温度应该基本等于热源出口温度这可以快速验证热力方程有没有写错。我个人做这类项目时始终坚持“从小系统跑通再到全系统验证”的开发习惯避免把问题积累到无法定位的程度。另外建议所有中间结果都存成结构体或表格式数据方便随时画收敛曲线和状态分布图。多能流计算迭代过程中的中间量画成曲线后很多时候一眼就能看出是哪个网络、哪个变量在“捣乱”这比盯着数值数组去排查效率高得多。最后分享一个我自己总结的开发顺序先把所有耦合系数设成0或固定运行值让三个网络完全解耦各自算各自的标准潮流确认三个子系统全部正确后再逐个打开耦合设备每加一个就重新验证一次收敛性和能量平衡。这个习惯帮我省下了大量排查时间强烈建议各位参考。 综合能源系统多能流计算：电气热耦合建模与Matlab代码实战
做区域综合能源系统课题的朋友应该都有这种体会单独算电网潮流、单独算天然气网络、单独算热力管网各自都有成熟工具和现成程序但一旦系统里加了热电联产机组、电锅炉、燃气锅炉这些“跨界设备”三个网络就再也拆不开了。我在这类项目里反复折腾过很久最直观的感受是多能流计算的难点不在于某一个网络本身而在于耦合元件怎么建模、三个网络的方程怎么组织在一起迭代收敛。这篇文章就把这块内容完整讲透从数学模型到Matlab代码架构再到算例验证和调试经验做一个偏实战向的拆解。这个内容适合正在做综合能源系统潮流计算、规划优化、运行调度方向的研究生和工程师尤其是手头需要自己编写代码实现电气热多能流计算的人。文章核心聚焦“计及多能耦合”这个关键点不会刻意去讲Matlab基础语法而是把精力放在“怎么把物理模型翻译成可运行的代码”这件事上。1. 多能流计算到底在解决什么问题1.1 三个网络不再是“各算各的”区域综合能源系统里电网、天然气网、热网原本是三个独立运行的能源输配网络。电网用交流潮流描述天然气网用节点压力方程和管道流量方程描述热网用水力-热力联合方程描述。单独看任意一个网络求解方法都很成熟电网有牛顿-拉夫逊法气网有节点法热网有水力/热力迭代法。问题出在“多能耦合”上。一个典型的园区级综合能源系统往往包含燃气轮机组加余热锅炉组成的CHP机组、电锅炉、燃气锅炉、P2G装置等设备。燃气轮机组消耗天然气发出电力的同时余热锅炉又把烟气余热转化为热网的热源电锅炉从电网取电却向热网供热燃气锅炉消耗天然气供热。这样一来电网的电力潮流会改变电锅炉的取电功率进而影响热网的注入热量热网的负荷变化又会影响CHP机组的热出力指令进而改变燃气轮机组的发电功率和天然气消耗量P2G的制气量又会影响气网的节点流量平衡。三个网络被这些耦合设备捆绑成了一个真正的“联立方程组”。我记得第一次自己搭这个模型的时候想省事先把电网潮流跑收敛然后把CHP电出力代到热网里再把热网热负荷代回电网结果循环迭代了三十几次都不收敛。后来才想明白这不是程序写法问题而是物理上电网、气网、热网之间存在双向反馈必须把耦合变量纳入统一的迭代框架或者把所有方程联立起来一起解。这个理解是这个项目里所有代码设计的前提。1.2 统一求解和分解求解的取舍实现多能流计算主要两条路线统一求解法和分解迭代法。统一求解法把所有子系统的方程组装成一个大型雅可比矩阵状态变量同时包含电网节点电压幅值、电压相角、气网各节点压力、热网各节点温度和流量一次牛顿迭代同时修正所有变量。这个方法的优点是收敛性通常更好不存在“内外层迭代”之间的误差传递缺点是雅可比矩阵规模大、稀疏结构复杂编程实现难度高尤其耦合元件对雅可比矩阵的贡献项很容易写错。分解迭代法是把电网、气网、热网分别用各自的求解器计算耦合元件在每个子系统里都表现为边界条件。比如电锅炉在电网侧是有功负荷在热网侧是热源二者通过 η_eb × P_eb Q_eb 这个关系式在外部迭代中建立联系。每轮迭代先算电网更新耦合设备功率再算气网再算热网检查耦合变量前后变化量是否小于阈值。这种方法模块化强调试方便可以复用现成的单网络潮流程序缺点是收敛速度可能慢某些情况下会振荡甚至发散。我的建议是初期代码实现优先走分解迭代法。原因是你能清楚地看到每一个物理量的变化过程出了问题容易定位是哪个网络、哪个耦合环节导致的。等把物理机制摸熟了想追求计算效率或者发高质量论文需要做大规模系统仿真时再升级到统一求解法。两条路线的核心区别可以用下面这个表来理解对比维度统一求解法分解迭代法状态变量维度所有网络变量同时求解维度大每个网络分别求解维度小雅可比矩阵需要组装跨网络耦合块各网络雅可比独立编程难度高耦合项易出错低模块独立性强收敛性较强但依赖好的初值可能慢或振荡调试友好度不友好出错难定位直观可逐步验证适用场景大规模系统、高耦合度、研究型工程方案快速评估、学习框架1.3 这套代码能做什么、不能做什么明确边界很重要。基于稳态模型的多能流计算解决的是“给定各网络的负荷和耦合设备运行状态求解系统稳态运行点”的问题。它能够用来做区域综合能源系统的规划方案比选、运行工况分析、耦合设备容量灵敏度研究也能给后续的优化调度模型提供初值或可行性校验。它不能做的是动态过程分析。电网的频率动态、天然气管道的气体波动传播、热网的慢热惯性过程都需要微分方程或偏微分方程描述不在稳态多能流的范围内。如果课题需要做短期动态仿真那应该基于时域仿真工具或者求解DAE方程组和这里讨论的稳态潮流是两套东西。2. 电气热子系统的模型怎么建才不踩坑2.1 电网部分极坐标NR法加功率不平衡量电网模型不必多解释重点讲代码实现里容易出问题的两个地方。第一采用极坐标牛顿-拉夫逊法时状态变量是节点电压幅值U_i和相角θ_i。对PQ节点同时列写有功和无功不平衡方程对PV节点只列有功方程松弛节点作为参考节点不参与迭代。功率不平衡方程为ΔP_i P_i_spec - U_i ∑_j U_j (G_ij cosθ_ij B_ij sinθ_ij)ΔQ_i Q_i_spec - U_i ∑_j U_j (G_ij sinθ_ij - B_ij cosθ_ij)角度差θ_ij θ_i - θ_j。Matlab代码实现时角度单位必须统一用弧度。有个很隐蔽的坑很多初学者在构造节点导纳矩阵的时候用复数矩阵直接计算没问题但在计算不平衡量时如果角度用了度整个雅可比矩阵数值会差57倍左右程序表现就是迷宫一样的乱跳。第二耦合设备参与电网潮流时它们在电网侧的表现要区分类型。CHP机组如果是定功率控制在电网里就是给定有功和无功的传统发电机组节点如果CHP参与电压调节就要按PV节点处理。电锅炉、P2G这类设备本质上是可控负荷电网建模里作为恒功率负荷节点处理即可。核心思路是耦合设备对电网的注入功率必须先更新再进入电网牛顿迭代不能滞后一个迭代周期。电网求解函数的主体结构大致是function [U, theta, iter] NRPowerFlow(bus, branch, gen, load) % bus: 节点参数结构体数组 % branch: 支路参数结构体数组 % 返回节点电压幅值U和相角theta Y buildYbus(bus, branch); % 构建节点导纳矩阵 % 初始化平启动 U ones(nBus, 1); theta zeros(nBus, 1); for iter 1:50 [dP, dQ] calcMismatch(U, theta, Y, bus, gen, load); if max(abs([dP; dQ])) 1e-6 break; end J buildJacobian(U, theta, Y); % 形成极坐标雅可比 dX J \ [dP; dQ]; % 稀疏LU分解求解修正量 % 更新状态变量 end end2.2 气网模型Weymouth方程与节点压力迭代天然气稳态管网模型比电网复杂一点核心方程是管道流量与两端压力的关系。工程中常用Weymouth方程描述高压输气管道中的稳态流量q_mn sign(p_m - p_n) × C_mn × sqrt(|p_m² - p_n²|)其中C_mn是管道常数由管径、长度、摩擦系数、气体组成决定。这个方程写代码时注意两点。第一取平方差要加绝对值符号由两端压力高低决定。第二压力必须用绝对压力不能用表压否则低压管网计算会出现开方负数的荒谬结果。气网的节点方程是流量平衡每个节点的流入流量等于流出流量加负荷消耗量。已知压力和未知压力节点分开处理通常气源节点压力已知负荷节点压力未知。未知压力节点用牛顿法求解修正方程的核心雅可比元素来自Weymouth流量对端点压力的偏导数∂q_mn / ∂p_m C_mn × p_m / sqrt(|p_m² - p_n²|)这个导数在压差趋近于零的时候会趋向无穷大这是气网迭代不收敛的常见原因之一。解决办法是给流量方程加一个很小的正则化项或者限制压力差的下限避免迭代过程中出现两个节点压力几乎相等的情况。管道流量的计算函数可以封装成这样function q gasPipeFlow(p_in, p_out, C) % 输入入口压力、出口压力、管道常数 % 输出标准体积流量正的表示从in流向out dp_sq p_in^2 - p_out^2; q sign(p_in - p_out) * C * sqrt(abs(dp_sq)); end2.3 热网模型水力与热力联合求解热网分成水力模型和热力模型两部分。水力模型描述管道流量分布节点满足质量流量守恒即流入节点的流量等于流出节点的流量加节点负荷消耗流量。环状管网还需要列写回路压降方程。实际项目中区域供热管网大多是环状结构但做多能流计算研究时往往先采用辐射状或树状拓扑这样可以显著降低水力求解难度。热力模型包含节点功率方程和管道温降方程。节点功率由供回水温差和质量流量决定Q_node c_p × G_node × (T_s - T_r)其中c_p是水的比热容G_node是流经负荷节点的质量流量T_s是供水温度T_r是回水温度。多个支路流入同一节点时混合温度按下式计算T_mix Σ(G_j × T_out_j) / ΣG_j这是热网模型里最容易忽略的一项。如果多条供热管道汇入同一个节点而你把每条管道的末端温度直接拿来当节点温度结果会偏差很大。管道温降描述采用苏霍夫公式T_end T_amb (T_start - T_amb) × exp(-λL / (c_p G))其中λ是管道单位长度热损系数L是管道长度T_amb是环境温度。这个公式说明流量越小温降越明显流量很小时管道末端温度会无限接近环境温度。数值实现上用这个公式直接计算末端温度即可但要注意G不能取零。热网求解的典型做法是“先水力,后热力再迭代”。原因是水力方程的解管道流量决定了热力方程中的流量参数而热力方程中的节点温度又会影响回水温度进而影响负荷侧的供回水温差计算。如果只做一次水力求解就固定流量再算热力对于不含回水温度迭代的简化系统是够用的但严格建模时建议做一轮外部迭代把供水温度和回水温度反复更新到收敛。3. 耦合元件建模多能流计算的灵魂环节3.1 CHP机组的两类运行模式CHP热电联产机组是综合能源系统里最重要的耦合设备。它消耗天然气同时输出电功率和热功率把电网和气网、热网三张网联系在一起。建模时先分清运行模式。背压式CHP机组的电热出力之比是固定的称为定热电比模式H_CHP C_m × E_CHP其中C_m是热电比常数。此时CHP的电力出力和热力出力完全耦合给定其中一个另一个就确定了。抽凝式CHP机组可以调整热电比运行范围是一个可行域需要满足电出力上下限、热出力上下限以及发电与供热之间满足燃料消耗约束。多能流计算中通常给定某个运行点然后校验是否落在可行域内。燃气消耗量是联系气网的纽带。CHP消耗的天然气流量F_CHP E_CHP / (η_e × LHV)其中η_e是CHP的发电效率LHV是天然气低位热值。热效率η_h和发电效率之和不等于总效率因为余热回收有损失通常总效率在0.75~0.9之间。在分解迭代法中CHP的作用流程是电网潮流算完后得到CHP的电出力E_CHP按热电比或设定运行模式更新热出力H_CHP然后在热网中把这个热出力作为对应节点的热源注入量同时按效率公式更新气网中该机组节点的天然气消耗流量。一步滞后带来的误差通过外部迭代消除。在统一求解法中CHP机组的建模要复杂一些电网方程、热网方程、气网方程之间的耦合项直接出现在统一的雅可比矩阵中。具体说热网中CHP节点的热源项对电网中E_CHP求偏导气网中CHP节点的气负荷项也对E_CHP求偏导这些偏导数都要作为非零元素填到雅可比矩阵的跨网络分块里。3.2 电锅炉、燃气锅炉、P2G三种耦合的简化处理电锅炉是把电转成热的设备Q_eb η_eb × P_eb在电网里它是一个有功负荷取电功率P_eb在热网里它是一个热源注入热功率η_eb × P_eb。这里效率η_eb通常在0.95~0.99之间电转热的损失远小于燃气锅炉。电锅炉的优势是响应速度快但在电网负荷高峰时段会加剧电力供需矛盾所以多能流计算中它的取电功率往往与电网运行状态互相制约。燃气锅炉只耦合气网和热网Q_gb η_gb × F_gb × LHV它从气网取气向热网供热建模最简单本质是两个网络之间单向耦合。燃气锅炉通常作为调峰热源热网热负荷超过CHP供热量时启用。P2G设备是实现电气双向耦合的典型装置消耗电能制取天然气氢气或甲烷。模型上简化为F_p2g η_p2g × P_p2g / LHV在电网里是电负荷在气网里是气源注入。P2G通常用于可再生能源消纳场景在风电、光伏出力大的时段把多余电力转化成天然气储存。多能流计算中P2G的功率是给定的运行指令不需要额外求解。这几种耦合设备的关系可以归成一张表设备输入侧输出侧效率典型值关键方程CHP机组天然气电热发电效率0.35~0.45热效率0.4~0.5H C_m×E, F E/(η_e×LHV)电锅炉电热0.95~0.99Q η_eb×P燃气锅炉天然气热0.85~0.93Q η_gb×F×LHVP2G设备电天然气0.50~0.70F η_p2g×P/LHV热泵电环境热热COP 3~5Q COP×P3.3 储能和冷负荷的处理方法完整的区域综合能源系统里通常还有蓄热罐、蓄电、吸收式制冷等设备。蓄热罐本质上是能量时移装置在稳态多能流计算中如果给定蓄热罐在某个计算时刻的充放热功率可以直接把它当作热网中一个负的或正的负荷节点处理并不复杂。冷负荷可以近似看成“负热负荷”或通过吸收式制冷机组转换。在热电联产园区里典型的做法是把冷负荷转化为热负荷需求因为吸收式制冷机以热源为输入Q_cold COP_c × Q_heat_in这样处理之后整个系统的热网方程结构不需要改变只是负荷节点多了一些等效热功率。如果课题明确要求电制冷和吸收式制冷并存那就需要在电网中增加电制冷机组的电负荷在热网中增加吸收式制冷的热负荷模型扩展方向是清晰的。4. Matlab代码实现架构与迭代流程4.1 数据结构和输入文件怎么组织写多能流代码最大的教训就是不要把所有变量都塞进脚本里硬编码。系统规模一扩大脚本方式会让你改一个管道参数找半天而且极易出错。我采用的方案是结构化定义三类网络和一个设备结构体。电网部分用bus和branch结构体数组每个节点包含编号、类型松弛节点/PV节点/PQ节点、有功负荷、无功负荷、电压幅值初值等字段每条支路包含首末端节点号、电阻、电抗、对地导纳。天然气网用gasNode和gasPipe结构体节点字段包括编号、类型气源/负荷、压力初值、节点流量注入管道字段包括两端节点编号、管道常数C。热网部分用heatNode和heatPipe结构体节点字段包括类型热源/热负荷/普通节点、供热温度、回水温度、质量流量注入管道字段包括两端节点号、长度、热损系数。耦合设备单独用unit结构体数组定义每个设备的字段包含设备类型、所在电网节点号、所在气网节点号、所在热网节点号、效率、热电比等。这样做的好处是主程序迭代时只需要遍历unit数组就能更新所有耦合边界完全不用改核心求解函数。4.2 主程序骨架分解迭代法的完整流程分解迭代法的核心逻辑是循环执行“电网求解→更新耦合设备热/气出力和消耗→热网求解→更新耦合设备电/气消耗→气网求解→检查耦合变量收敛”。下面给出一套可以直接参照的主循环框架function x solveMEF(d) % d: 系统数据结构包含电网、气网、热网、耦合设备定义 x.E initElectric(d.E); % 电网初值U1, theta0 x.G initGas(d.G); % 气网初值所有节点压力基准压力 x.H initHeat(d.H); % 热网初值供水温度T_s0, 回水温度T_r0, 流量初值 for k 1:100 x.E_old x.E; x.G_old x.G; x.H_old x.H; % ---- 步骤1求解电网潮流 ---- % 根据上一轮的气网/热网状态更新耦合设备在电网中的注入/消耗功率 [P_inj, Q_inj] updateElectricBoundary(x, d); x.E NRPowerFlow(d.E, P_inj, Q_inj, x.E); % ---- 步骤2更新热网边界并求解热网 ---- % CHP热出力、电锅炉热出力根据最新电网结果更新 [Q_heat_source, G_load] updateHeatBoundary(x, d); x.H HeatFlow(d.H, Q_heat_source, x.H); % ---- 步骤3更新气网边界并求解气网 ---- % 燃气锅炉、CHP的用气量以及P2G的产气量都要更新 [F_gas_load, F_gas_source] updateGasBoundary(x, d); x.G GasFlow(d.G, F_gas_load, F_gas_source, x.G); % ---- 步骤4收敛判据 ---- dE max(abs([x.E.U - x.E_old.U; x.E.theta - x.E_old.theta])); dG max(abs(x.G.p - x.G_old.p)); dH max(abs([x.H.T_s - x.H_old.T_s; x.H.T_r - x.H_old.T_r])); if max([dE, dG, dH]) 1e-6 fprintf(多能流计算收敛迭代次数%d\n, k); return; end end error(多能流计算不收敛请检查初值或系统参数); end这个主循环里有几个细节值得注意。迭代上限不要设太大100次以内足够判断收敛性如果100次还不收敛大概率是模型或初值的问题继续迭代也没有意义。收敛判据用状态变量的变化量而不是不平衡量更稳妥。耦合变量的更新用最新值即前一个网络刚算出来的结果立刻用于下一个网络的边界更新这种风格叫“Gauss-Seidel式”更新比全部用旧值更新的“Jacobi式”收敛更快。4.3 热网求解函数的关键实现热网求解函数内部要处理水力迭代和热力迭代的嵌套关系。一个简化的实现思路是先假设各节点供水和回水温度已知求解节点质量流量然后固定流量更新各管道末端温度和各节点混合温度检查温度变化量如果不收敛则用新温度再次求解水力方程。节点质量流量求解的核心是节点流量守恒方程写成矩阵形式就是A × G G_load其中A是节点-支路关联矩阵G是管道质量流量向量G_load是节点注入流量向量。对树状管网这个方程直接就解出来了环路管网需要额外列写回路压降方程联合求解。温度迭代部分的核心代码如下function T_node updateNodeTemp(G_matrix, T_in, T_amb, pipes, cp) % G_matrix: 各管道质量流量向量 % T_in: 管道入口温度向量 % 先用苏霍夫公式计算各管道末端温度 for i 1:length(pipes) T_out(i) T_amb(i) (T_in(i) - T_amb(i)) * ... exp(-pipes(i).lambda * pipes(i).L / (cp * G_matrix(i))); end % 再按混合定律求各节点的混合温度 for k 1:nNode inflow find(tailNode k); % 流入节点k的管道集合 numerator sum(G_matrix(inflow) .* T_out(inflow)); denominator sum(G_matrix(inflow)); if denominator 0 T_node(k) numerator / denominator; end end end这里有个经验性的处理分母流量很小时混合温度数值会不稳建议对denominator加一个很小的保护阈值或者直接将流量小于某个下限的管道从混合计算中剔除否则温度会出现不合理的跳变。5. 算例设计与结果验证5.1 测试系统的搭建思路为了验证代码的正确性需要一个不复杂但包含典型耦合元素的测试系统。我用的是下面这种结构电网采用5节点系统其中节点1是松弛节点外电网等值节点2接入一台CHP机组节点3接入电锅炉节点4和节点5是常规电负荷。天然气网采用4节点结构节点1是气源压力给定节点2连接CHP机组气负荷节点3连接燃气锅炉气负荷节点4连接P2G设备气源。热网采用6节点结构节点1是CHP热源节点节点2是电锅炉热源节点节点3是燃气锅炉热源节点节点4和节点5是热负荷节点节点6是回水汇集节点所有管道采用辐射状连接。测试系统的运行场景设定为冬季典型工况CHP机组带基础电负荷和热负荷电锅炉作为补充热源燃气锅炉调峰。P2G设备以小功率运行。这样能保证三个网络之间至少有四个耦合元件在同时工作耦合关系够复杂。仿真参数的选取也要符合实际。电网基准容量取100 MVA电压基准取10 kV。气网基准压力取1 MPa天然气低位热值取35.8 MJ/Nm³。热网供回水温度基准取130℃/70℃工程中高温水系统常见参数水的比热容取4.2 kJ/(kg·K)。所有效率和热电比参数设置要保证系统能量平衡有意义不能让某个设备效率大于100%。5.2 收敛过程和结果评估实际运行这套代码时典型工况下分解迭代法大约需要15到30次外部迭代收敛到1e-6的精度。若从平启动初值开始前五轮迭代耦合变量变化量较大之后快速下降基本呈现线性收敛趋势这也符合Gauss-Seidel型迭代的特征。收敛后可以校核几个总量指标全系统电功率平衡即发电机出力与外电网交换功率之和等于负荷加电锅炉和P2G耗电热网供热量等于热负荷加管道散热损失气源总供气量等于CHP、燃气锅炉耗气量减去P2G产气量。能量平衡校验中每一项的误差在1e-6量级说明方程组装和迭代过程是正确的。我还建议做一步“解耦一致性验证”将电锅炉、燃气锅炉、P2G的运行功率全部设为0热网只有CHP供热这时把CHP的电出力和热出力固定将CHP节点当作电网的PQ节点气网节点当作普通负荷此时多能流程序的电网部分结果应该和单独运行电网潮流程序完全一致。这一步验证过关代码的电网模块基本可信。同理也可以分别验证气网、热网模块。5.3 耦合强度对收敛性的影响在代码跑通之后可以进一步做耦合强度的灵敏度分析这也是论文里常用到的手段。方法是固定其他参数逐步改变CHP热电比或者电锅炉容量观察外部迭代次数的变化。一般情况下热电比越大CHP电出力到热出力的传递越强电网和热网之间的耦合越紧密分解迭代法需要的迭代次数会上升。如果热电比超过某个阈值甚至可能出现振荡不收敛的情况。我实测遇到过一个案例把电锅炉功率从1 MW逐步增加到8 MW时迭代次数先是缓慢增加超过6 MW后突然发散。原因是电锅炉功率增大导致电网节点电压偏低电网潮流求解精度下降反过来又影响热网边界形成正反馈。解决办法是调整电网无功补偿容量或者将电锅炉连入的节点改为更靠近电源的节点。这个现象本身也是一个有价值的研究点很多文献专门讨论过耦合强度与收敛性的关系。6. 常见问题与排查技巧6.1 不收敛先看初值再看雅可比最后看单位多能流计算不收敛的排查顺序很重要我自己的经验是按“初值→雅可比→单位”的顺序来切勿一开始就怀疑算法框架。初值问题是新手最容易犯的错。气网压力初值如果设得过高或过低Weymouth方程的流量初值会离谱导致气网牛顿迭代第一步就失败。一个稳妥的做法是把所有未知压力节点初值设为气源压力的0.95倍左右不要平启动设成1 MPa以下一大堆低压节点。热网温度初值也不要取环境温度供水温度初值可以直接用设计供温130℃回水温度用70℃这样热力迭代的温度差不会太大。雅可比矩阵的问题在统一求解法中比较常见。排查手段是有限差分验证将状态变量整体加一个1e-6的小扰动重新计算不平衡量与解析雅可比矩阵对比。如果某个非零元素对不上定位到具体方程后重点检查对应耦合元件的偏导数推导。单位问题在气网和热网中尤其明显。天然气流量可以用标准体积流量Nm³/h也可以用质量流量kg/s两者之间要乘标准密度0.717 kg/Nm³左右。功率和流量换算时还要除以3600转换秒和小时。如果不统一单位能量平衡校验一定对不上。6.2 一个实际排查案例气网压力振荡有一次跑算例气网迭代在第四轮之后开始振荡压力值在1.0 MPa和0.92 MPa之间来回跳动。排查发现是P2G设备在电网侧取电功率过大导致电网潮流计算出的P2G电耗偏高按效率换算成产气流量后这个流量加到气网节点上气源节点压力被推高下一轮气网求解又把该节点负荷调低造成振荡。处理办法是在外部迭代中引入阻尼系数即本轮使用的气负荷值取上一轮值和本轮更新值的加权平均F_new α × F_calculated (1-α) × F_oldα取0.6左右阻尼效果明显迭代很快恢复收敛。这种方法虽然只是工程技巧但在很多耦合度较高的系统里非常管用。6.3 从“能跑”到“可信”的三个验证手段程序写完能出结果并不代表可信。要证明代码正确至少要过三关。第一关是有限差分验证雅可比矩阵这保证方程组的求解方向是正确的。第二关是能量平衡校验保证所有功率、流量在物理上是守恒的。第三关是对比验证找一个同类文献的算例系统用你的代码重跑一遍看结果能否和文献给出的电压、压力、温度分布对上。还有一个非常实用的小技巧把热网的管道热损系数暂时设成0环境温度设成不影响温降的数值此时热网供水温度应该基本等于热源出口温度这可以快速验证热力方程有没有写错。我个人做这类项目时始终坚持“从小系统跑通再到全系统验证”的开发习惯避免把问题积累到无法定位的程度。另外建议所有中间结果都存成结构体或表格式数据方便随时画收敛曲线和状态分布图。多能流计算迭代过程中的中间量画成曲线后很多时候一眼就能看出是哪个网络、哪个变量在“捣乱”这比盯着数值数组去排查效率高得多。最后分享一个我自己总结的开发顺序先把所有耦合系数设成0或固定运行值让三个网络完全解耦各自算各自的标准潮流确认三个子系统全部正确后再逐个打开耦合设备每加一个就重新验证一次收敛性和能量平衡。这个习惯帮我省下了大量排查时间强烈建议各位参考。