做了大半年区域综合能源系统电气热能流计算，从最初只会单独跑电网潮流，到把气网、热网全部耦合起来在Matlab里写成一套能稳定收敛的代码，中间踩过的坑比想象中多得多。这个项目的完整标题是《计及多能耦合的区域综合能源系统电气热能流计算研究（Matlab代码实现）》，核心就一句话：当电、气、热三种异质能源通过热电联产（CHP）、燃气锅炉、电锅炉这些设备耦合在一起时，怎么用一套可复现的数值算法，算出整个系统稳态运行的各种关键状态量。

很多初学者以为这项任务就是把三个网的潮流程序简单堆在一起，实际上完全不是这么回事。电网算的是节点复功率平衡，气网算的是节点压力和管段流量平衡，热网更麻烦，水力计算和热力计算得交替进行。三种能网的时间常数、量纲、建模思路都不一样，耦合变量又来回反馈，稍不注意就会发散。这篇博文把我从模型搭建、方程推导到Matlab落地实现、调参避坑的完整过程梳理一遍，内容包括电气热网稳态模型、耦合设备建模、顺序迭代求解框架、雅可比矩阵组装、典型算例分析，以及实际调试中遇到的一堆典型问题。适合正在做综合能源系统、园区多能互补、微电网方向研究的研究生和工程师参考。

1. 思路拆解：多能耦合到底“耦”在哪里

1.1 三种能网的计算范式差异

先看基础。电网的稳态潮流分析，核心是节点功率平衡方程，已知负荷和机组注入功率，求解各节点电压幅值和相角。方程是光滑的二次型，牛顿-拉夫逊法收敛性很好，只要初值不离谱，四五次迭代基本都能收住。这是三种网络里最成熟、参考资料最多的一个。

气网的稳态计算就不太一样了。它要满足的是节点气体流量平衡，即流入某个节点的天然气流量减去流出流量，等于该节点的用气负荷。描述管段流量的公式常用Weymouth方程，这个方程带压力平方差以及绝对值符号，在管道两端压力相等的位置导数不存在，容易引发迭代震荡。我实际调试时发现，气网方程整体是单调的，但收敛半径比电网小不少，初值差一点就可能在压力低的位置出现负的压力平方，直接报错。

热网是最麻烦的。它的稳态运行状态需要分两步描述：水力工况决定各管段质量流量，热力工况决定各节点供水和回水温度。水力方程有两大类，一类是节点流量平衡方程，流入流出必须相等；另一类是环路压头平衡方程，也就是沿闭合回路流动的压头损失代数和为零。流量在方程里以平方项出现，所以是非线性的。热力方程又建立在流量基础之上，温度计算时流量待在分母上，初值给不好一样发散。

三张网络合在一起，本质上是一个多物理场耦合的非线性方程组问题。我自己的体会是：如果只把三套程序按顺序各跑一遍，那不叫多能耦合计算，那只是三个独立的工具串行使用。真正的耦合计算要求三个子系统之间通过耦合设备逐轮交换边界条件，一直迭代到整个系统的运行状态达到一致。

1.2 耦合设备是系统的“关节”

区域综合能源系统里，电网、气网、热网并不是孤立存在的，它们通过能量转换设备连在一起。常用的耦合设备有热电联产机组CHP、燃气锅炉、电锅炉、热泵，甚至电转气P2G。这类设备在能流计算里扮演的角色很清晰：上游能网给它们提供输入能量，下游能网承接它们的输出能量，转换关系则由各自效率决定。我用一个表格把最常见的几种设备整理了出来，方便后续建模时直接对照：

从能流角度看，任何一个环节负荷的变化都会沿整个网络链式传导。举个例子：冬季某区域热负荷升高，CHP机组热出力随之上调，天然气耗量增加，气网节点压力就会下降；与此同时CHP电出力也发生变化，电网潮流分布跟着改变。如果热负荷需求继续增长，电锅炉也可能被调度，电网负荷进一步变化。这就是多能耦合的本质——所有变量相互制约、相互传导，任何一个边界条件的改变都可能在另一个网络里引起不可忽略的响应。

2. 电气热数学模型怎么搭

2.1 电网潮流的极坐标牛顿法

电网潮流我采用的是经典的极坐标形式牛顿-拉夫逊法。对第i个节点，有功和无功功率方程分别为：

Pi = Vi * Σ_j Vj * (Gij cosθij + Bij sinθij)

Qi = Vi * Σ_j Vj * (Gij sinθij - Bij cosθij)

其中Gij和Bij是节点导纳矩阵中的电导和电纳，θij等于θi减θj。把给定的节点注入功率与计算功率做差，得到不平衡量ΔP和ΔQ，然后对电压幅值和相角求偏导，形成雅可比矩阵的四个分块：

Hij = ∂Pi / ∂θj，Nij = Vj * ∂Pi / ∂Vj，Jij = ∂Qi / ∂θj，Lij = Vj * ∂Qi / ∂Vj

写成矩阵形式就是：

| ΔP | = | H N | * | Δθ | | ΔQ | | J L | | ΔV/V |

在Matlab里，最核心的是两点：第一个是正确组装节点导纳矩阵Ybus，第二个是雅可比矩阵的更新。Ybus我建议用稀疏矩阵存储，省内存而且后续矩阵求逆速度快。雅可比矩阵的组装可以完全向量化，把同类型节点和支路的贡献批量算出来，避免一层层for循环把速度拖垮。

2.2 气网管流的Weymouth方程

天然气网络稳态计算我建议用节点压力平方π作为状态变量。对连接节点i和节点j的管道，管段流量写成：

Qij = sign(πi - πj) * Cij * sqrt( |πi^2 - πj^2| )

其中Cij是管道常数，由管径、长度、摩擦系数等因素共同决定。用压力平方做变量有一个直接好处：Weymouth方程里出现的是πi和πj的一次差，而不是压力的二次方，求导过程简单很多。而且流量方向由πi和πj的相对大小自动决定，不需要人为指定。

节点流量平衡方程为：

Σ_j Qij - Li = 0

其中Li是节点i的用气负荷，气源节点则作为压力松弛节点处理，给定压力值作为系统的参考。牛顿法迭代时，对每个非松弛节点建立上述平衡方程，求管道流量对π的偏导填充雅可比矩阵。气网的雅可比元素带有分母项，即平方根项，如果初值设置导致某个管道两端压力平方差接近零时，这一项就会变得很大，这是气网数值不稳的主要来源之一。

2.3 热网的“水力+热力”双层结构

热网稳态计算必须分成水力工况和热力工况两步来走，这个顺序不能颠倒，因为热力计算需要的水流量必须由水力计算先确定。

水力工况方面，第一步列出节点流量平衡方程：A * q = L，A是节点-支路关联矩阵，q是各管段质量流量，L是节点注入或流出的质量流量。第二步列出环路压头平衡方程：B * (K * q * |q|) = 0，B是环路-支路关联矩阵，K是管段阻力系数。这里关键是要保证独立环路数量足够，程序里通过环路矩阵自动识别闭合回路。

热力工况方面，对每个负荷节点有热功率方程：

Φi = Cp * mi * (Ts_i - Tr_i)

其中Cp是水的比热容，mi是通过负荷节点的质量流量，Ts和Tr分别表示该负荷节点的供水温度和回水温度。管道的热量损失用指数温降模型描述：

T_end = Ta + (T_start - Ta) * exp(-λ * L / (Cp * m))

λ是管道单位长度热损失系数，L是管长。当多根管道的水汇入同一个节点时，混合后的温度按流量加权平均计算：

T_out = Σ(m_in * T_in) / Σm_in

这一条在供水网和回水网里都有，是热网装配方程组时比较容易出错的地方。我在程序里把供水网络和回水网络分开存储，通过负荷节点把它们串起来，这样逻辑清晰，调试起来也好定位问题。

2.4 耦合设备的数学表达

耦合设备建模其实没有太多花样，核心就是把输入能量按效率映射为输出能量。以CHP机组为例：

Pe = ηe * Fgas Qh = ηh * Fgas

其中Fgas是天然气输入的化学能功率，ηe和ηh分别是发电效率和产热效率。实际运行中CHP存在一个电热比概念r = ηe / ηh，调度人员在算例中通常会给定这个比例或者给定热出力让程序反推电出力。燃气锅炉和电锅炉的建模更简单，燃气锅炉输出热功率等于天然气输入乘以热效率，电锅炉输出热功率等于电输入乘以转换效率。热泵则用COP来描述，输入1kW电可以输出几kW热。

耦合设备在顺序迭代里承担的角色就是变量转换的桥梁：CHP消耗的天然气量变成气网负荷节点，CHP电出力变成电网电源节点，热出力变成热网热源节点。每个设备内部还要考虑出力上下限约束，我在代码里加了简单的限幅逻辑，迭代过程中如果计算出的出力越过界限，就把它钳制在边界上，防止不现实的结果污染整个潮流计算。

3. 求解策略与Matlab实现细节

3.1 统一求解还是顺序迭代

针对多能流计算，学术界和工程界主要有两种求解思路：统一求解法和顺序迭代法。

统一求解法是把电网、气网、热网的所有方程组装进一个超大雅可比矩阵，所有的状态变量同时迭代更新。理论上收敛性最好，能处理强的耦合关系，而且各网络之间的交互可以在一次迭代中完全体现。但代价是代码复杂度飙升，矩阵规模大，一旦不收敛，问题定位非常痛苦，新手不建议直接上手。

顺序迭代法是按电网→热网→气网的顺序逐个求解，每个网络内部用成熟的牛顿法，网络之间通过耦合设备更新边界条件，外循环反复进行直到所有网络的变量不再变化。这种方法实现简单，模块化程度高，能复用成熟可靠的电网潮流代码。缺点是如果耦合特别强，外循环可能出现震荡甚至不收敛。

我做这个项目时用的就是顺序迭代法。原因很直接：项目目标是先算得对、算得稳，再谈算得快。实测下来，对于节点规模几十个左右的区域综合能源系统，顺序迭代的收敛速度完全够用，总计算时间在普通笔记本上也就是一两秒的数量级。如果后续要处理大量CHP占比很高的强耦合场景，再考虑升级成统一求解法也不迟。

3.2 数据结构与雅可比矩阵组装

Matlab代码里我采用结构体数组来管理电网、气网和热网的拓扑与参数，每个网络保持独立的数据文件，这样换算例时只需要替换数据，不需要动求解核心。电网节点结构体包含节点编号、节点类型、有功负荷、无功负荷、初始电压和相角；气网节点包含编号、类型、用气负荷、初始压力；热网节点包含编号、类型、热负荷、水流量等。

%% 电网节点数据示例 grid.bus = struct('id', 1, 'type', 1, 'Pd', 0.1, 'Qd', 0.05, ... 'V0', 1.0, 'theta0', 0, 'Pg', 0, 'Qg', 0); grid.branch = struct('from', 1, 'to', 2, 'r', 0.02, 'x', 0.06, 'b', 0.01); %% 气网节点数据示例 gas.node = struct('id', 1, 'type', 1, 'L', 0.3, 'p0', 4.0); gas.pipe = struct('from', 1, 'to', 2, 'C', 0.08); %% 热网节点数据示例 heat.node = struct('id', 1, 'type', 1, 'Phi', 0.8, 'Cp', 4.182); heat.pipe = struct('from', 1, 'to', 2, 'K', 8000, 'L', 200, 'lambda', 0.25); heat.loop = [1, 2, 3]; % 独立环路包含的支路编号

雅可比矩阵组装是整个数值计算的心脏，这里有一个非常重要的工程建议：组装过程和求解过程分离。先把雅可比矩阵通过向量化计算一次性填好，再用\运算符求解线性方程组，不要让求解器在内部重复计算矩阵。电网的雅可比矩阵我这里直接列一个简化的组装示意：

Ybus = makeYbus(grid.branch); % 构建节点导纳矩阵 [H, N, J, L] = calcJacobian(V, theta, Ybus); dP = Pspec - Pcalc(V, theta, Ybus); dQ = Qspec - Qcalc(V, theta, Ybus); dtheta = H \ (dP - N * dV_over_V); dV = V .* (L \ (dQ - J * dtheta));

热网水力雅可比矩阵要特别注意流体阻力项K * q * |q|对q的导数，实际导数形式是K * 2 * |q|。如果某个管段流量初值恰好设成零，这一项直接为零，雅可比矩阵奇异，整个求解直接失败。这个问题我在后面排查章节会详细展开。

3.3 牛顿法迭代与收敛判据

顺序迭代的外循环流程分为五步：

1. 初始化所有网络的状态变量，包括电网电压幅值和相角、气网节点压力平方、热网管段流量和节点温度。
2. 求解电网潮流，计算CHP等电源节点的新电出力，更新电网节点注入功率。
3. 求解热网水力方程得到各管段流量，再解热力方程得到各节点供回水温度，更新CHP热出力。
4. 把CHP和燃气锅炉消耗的天然气量叠加到气网负荷节点上，求解气网节点压力。
5. 检查前后两次外循环各耦合变量的变化量，如果都小于容差则收敛，否则回到第二步继续迭代。

收敛判据我设置为两个条件必须同时满足：一是每个网络内部的功率或流量不平衡量小于设定容差，电网取1e-6标幺值，气网取1e-4相对误差，热网温度变化小于1e-4摄氏度；二是耦合变量两天循环之间的相对变化小于1e-5。两个条件同时满足才宣告整个多能流计算收敛，只满足一个很容易出现假收敛。

耦合变量的更新不要一步到位，我实际调参发现加一个阻尼因子很有用。具体做法是不直接采用本轮求解出的CHP热出力数值，而是取上一轮数值和本轮数值的加权平均，权重在0.6到0.8之间。这个操作和电网潮流里的步长控制思路类似，用牺牲一点收敛速度的代价换取更强的稳定性，尤其在热负荷在迭代过程中变化剧烈时效果明显。

3.4 代码整体框架概览

整个项目的主函数我设计成模块化结构，每个子系统对应一个独立求解函数，外循环只负责耦合变量传递和收敛判断。这样做最大的好处是，把电网潮流换成别的算法或者把气网模型升级为非等温模型时，不用动整个框架，只需要替换对应的子系统函数。

function [grid_res, gas_res, heat_res] = ries_energyflow(data) % 初始化 grid_state = init_grid(data.grid); gas_state = init_gas(data.gas); heat_state = init_heat(data.heat); % 耦合设备参数 chp = data.chp; gb = data.gas_boiler; eb = data.electric_boiler; % 外循环 for outer = 1:max_outer % Step 1: 电网潮流 grid_state = solve_grid(grid_state, data.grid, chp.Pe, eb.Pe); % Step 2: 热网计算 chp.Qh = chp_calc(chp.Fgas, chp.eta_e, chp.eta_h); heat_state = solve_heat(heat_state, data.heat, chp.Qh, gb.Qh); % Step 3: 气网计算 chp.Fgas_new = (chp.Pe / chp.eta_e + chp.Qh / chp.eta_h) / 2; gas_state = solve_gas(gas_state, data.gas, chp.Fgas_new, gb.Fgas); % Step 4: 阻尼更新耦合变量 chp.Fgas = 0.7 * chp.Fgas + 0.3 * chp.Fgas_new; % 收敛判断 if check_convergence(grid_state, gas_state, heat_state, tol) break; end end end

4. 算例设计与结果分析

4.1 测试系统与计算条件

我用的测试系统是自己拼装的区域综合能源系统算例，规模不算大，但足以说明问题。电网部分参考了IEEE 33节点配电系统的拓扑骨架，加了两台分布式燃气轮机和若干负荷；天然气网络取7个节点、8条管道，带一个气源节点和两条环状结构；热网由16个节点组成，包含3个独立环路，负荷节点覆盖供热区域中的主要热力站。

基准值统一取为：电功率100MW、气流量100kg/s、热功率100MW；压力用MPa作单位，温度用摄氏度。CHP机组共3台，燃气锅炉1台，电锅炉1台。算例场景设定为典型冬季供暖日：热负荷32MW，电负荷18MW，CHP承担约七成热负荷，其余由燃气锅炉和电锅炉分担，这是北方园区供能系统很常见的一种运行方式。

4.2 迭代收敛过程实录

在这个算例下，顺序迭代外循环跑到第14次才收敛。单次外循环内，电网潮流平均迭代3次收敛，热网水力方程迭代10到12次、热力方程迭代6到8次，气网方程迭代4到5次。把各子系统独立的迭代次数和整个外循环迭代次数乘起来，总计算时间在普通笔记本上约为0.9秒。对于规划阶段的分析计算，这个性能完全够用。

我特意记录了第一次和第二次外循环之间各个状态变量的变化幅度。刚开始时，各网络状态量变化很大，因为耦合变量初值只是我拍脑袋给的估算值，和实际收敛值差别明显。到了第五次外循环以后，各耦合变量的变化幅度就稳定下降，呈现典型的近线性收敛趋势。最终收敛时，电网最大功率不平衡量约8×10⁻⁷标幺值，气网最大节点流量不平衡量约5×10⁻⁶，热网温度最大变化约0.0002摄氏度，都在容差范围内。

4.3 耦合效应对结果的影响

为了说清楚耦合计算到底和多能流独立计算有多大差别，我把同一个算例里CHP的天然气耗量固定住，不让它随气网压力和气负荷反馈变化，相当于做一次解耦近似计算。两种方式的结果对比非常能说明问题：解耦后气网最低节点压力是3.71MPa，耦合计算后同一节点压力降到3.58MPa。原因不难理解，CHP实际消耗的天然气比固定值更大，重负荷条件下气网管道压降因此明显增加。

电网侧的差异同样存在。耦合计算后，电网最低节点电压从解耦近似时的0.9628标幺值降为0.9547标幺值。幅度虽然只有不到一个百分点，但如果后续要在同一框架内评估无功补偿容量或者电压合格率，这个差异就不能再忽略了。这也说明了一个重要的工程结论：系统的运行状态是多能耦合共同决定的，单独计算任一能源网络都可能高估运行裕度。

5. 实操中的坑与排查清单

5.1 初值敏感性是最大的敌人

这项任务做久了会发现，算法本身不复杂，真正让人崩溃的是初值。电网潮流相对好办，平启动电压1.0∠0°基本能覆盖大多数情况。气网初值统一取标称压力附近即可，通常不容易出问题。热网是重灾区，我踩过最惨的一次坑是把管段流量初值全部设为0，心想反正迭代会自己收敛，结果牛顿法直接在第一次雅可比求解时崩溃——前面说过，阻力项导数中含有2乘以|q|，流量为零时这项等于零，雅可比矩阵奇异，线性方程组无解。

后来我养成了一个习惯：热网各管段流量初值按负荷功率反推。方法是m = Φ / (Cp * ΔT)，其中ΔT取设计供回水温差，比如95/70摄氏度对应的25摄氏度温差。这样即使用户给的负荷分布不太均匀，初值也至少处在正确的数量级上，水力方程从一开始就有合理的梯度信息，迭代收敛速度明显变快。

5.2 量纲统一问题

电气热三种网络的量纲差异在耦合时是个大坑。电网用标幺值很顺，气网用MPa压力很自然，热网用kg/s和摄氏度很直观，但三套量纲混在一起后，耦合接口处必须做严格的基准转换。比如CHP的天然气耗量在电网模型里用MW表示化学能功率，但进入气网节点负荷时必须换算成kg/s。我处理的方式是定义统一的基准转换函数，在每个耦合接口调用，绝不在设备代码里手写换算系数，防止不同工况下弄混单位。

如果不做这个转换，雅可比矩阵中不同方程的量纲相差好几个数量级，数值条件数很大，求解过程极易出现精度丢失。实际表现就是迭代曲线震荡、收敛变慢，严重时直接报矩阵奇异。检查这类问题有一个快速方法：把雅可比矩阵对角元素打印出来，如果同一个矩阵里元素数量级差异超过10的6次方，基本可以断定量纲或基准没有统一。

5.3 Matlab环境与性能优化

代码我最早在R2023b上开发调试，后来用R2026b的兼容模式也跑通了，Matlab官方版本间基础语法这部分兼容性做得还可以。开发阶段有个建议是先把算例在串行模式下跑通，再做性能优化。我一开始直接开了parpool并行池，结果调试时根本分不清是算法问题还是并行任务调度的干扰，排查效率非常低。

性能优化方面，最有价值的操作是把能向量化的循环全部向量化，同时用sparse函数存储稀疏矩阵。电网节点导纳矩阵、气网节点关联矩阵、热网节点-支路关联矩阵本身都是高度稀疏的，用稀疏存储后，矩阵求逆和乘法的耗时几乎可以忽略不计。我在32个电网节点、7个气网节点、16个热网节点的规模下，整个多能流计算保持在1秒以内完成，这对后续做参数分析或者优化调度外层循环来讲非常关键。

还有一个容易被忽略的点是Matlab环境本身。如果启动时出现license manager error -8之类的授权异常，这类报错通常是许可证路径没有被正确识别，或者网络浮动授权服务器不可达。优先检查系统环境变量里的许可证路径配置和授权文件位置，不要急着重装。矩阵规模不大时CPU单核性能比多核并行重要，所以不要一上来就依赖并行计算。

5.4 常见故障速查表

我把调试过程中遇到的典型问题整理成了一张速查表，新接触这套代码的人遇到状况可以先对照排查：

6. 一些个人经验体会

这个项目做下来，我最大的感触是，多能耦合能量流计算，真正的技术含量不在算法推导本身，而在于工程调试的耐性和系统性。我花在调初值、查矩阵组装、核对量纲上的时间，远远超过写方程的时间。如果你也要做类似的工作，我强烈建议按这个顺序推进：先把电网潮流跑通并和标准算例比对，再把气网单独算通，把热网的水力工况单独调稳定，最后才做三者耦合。每个环节都先用一个非常小的两节点系统验证结果，确认和手算值一致以后再进行下一步。这套流程听起来慢，实际上是最快的一条路，它能帮你在源头上堵住大部分bug。

另外还有一个值得分享的小技巧：调试耦合迭代时，可以把每一轮外循环的耦合变量打印出来，比如CHP的天然气耗量、热出力、电出力，观察它们的演化轨迹。如果发现变量在长期震荡而不是单调收敛，说明阻尼因子太小或者某个子系统内部的容差设置太松；如果变量直接从初值跳到一个不合理的边界值，说明耦合设备限幅逻辑可能有问题。这比单纯盯着发散信息排查效率高得多。

这套电气热耦合能流计算框架目前的扩展空间还很大。比如把气网模型升级成非等温模型，把热网加入管道热惯量动态特性，又或者在耦合设备层加入储能装置和需求侧响应逻辑，都能在这个框架上做下去。代码本身也写了完整的子函数模块，后续做参数优化或者多场景扫描时，只需要在外层套一个循环就行。如果你正在做园区多能互补、分布式能源系统规划或者综合能源调度相关的课题，这套方法和代码可以作为你的第一版基础工具。