计及多能耦合的区域综合能源系统电气热能流计算研究（Matlab代码实现）
做综合能源系统仿真这几年，我接触最多的一个方向就是多能耦合系统的稳态能流计算。很多人一听“电气热能流”就以为是把电网潮流、天然气网水力计算、热网水力热力计算三个东西分别跑一遍然后拼在一起，真上手做一遍就会发现，事情没那么简单。一个区域综合能源系统里，CHP机组同时连着电网和气网，电锅炉同时连着电网和热网，燃气锅炉同时连着气网和热网，这些耦合设备把三条网络硬生生绑成了一个整体。你单独算电网潮流的时候，气网供给CHP的燃料量是未知的；单独算气网的时候，CHP的耗气量又取决于电网的发电需求；而热网的供水温度又反过来影响CHP的抽汽量。
这篇内容就是基于我实际做过的一个Matlab仿真项目，把整套电气热能流计算从数学模型、求解策略到代码实现完整梳理一遍。适合正在做综合能源系统方向的研究生、刚入行做园区能源规划的工程师，以及想把多能流计算模型真正落地成可运行代码的人。我会把每一步为什么这么做、踩过哪些坑、参数怎么定，都写清楚。
1. 问题是先理清楚的：为什么单算"三条网"不算真正意义的多能能流
1.1 单能流系统已经有成熟工具，为什么还要做联合能流计算
先说个最直观的类比：你手头有一个电网潮流程序、一个天然气网水力计算程序、一个热网水力热力计算程序，三个程序各自都跑得很顺。但把它们拼在一起的时候，问题就出来了。电网潮流需要知道气网供给CHP的天然气量，从而确定CHP的发电功率和热功率；气网需要知道CHP的耗气量才能做节点流量平衡；热网需要知道CHP的供热量才能确定节点温度和管道流量。这三个量互为因果，谁都没法先独立确定。
这就是所谓的多能耦合的内在复杂性。严格来说，只有当耦合设备两侧的变量都被当作未知数、放进同一个迭代框架里求解，才能叫计及多能耦合的能流计算。如果只是把CHP当作电网里的一个PQ节点，气网里当作一个固定负荷，那本质上还是三个孤立的单能流问题，谈不上耦合。
另外，单能流系统的求解工具已经非常成熟，比如电网潮流有Matpower，天然气网有管网仿真软件，热网也有专门的热力计算工具。但找遍市面，没有一个现成工具能一次性求解电、气、热三网联立的能流方程。原因其实也简单：这三个网络的物理特性差异太大。电网是代数方程主导，求解讲究收敛性和电压约束；气网是压力驱动的非线性方程，跟电网形式有相似之处但物理含义完全不同；热网更特殊，它被水力方程和热力方程双重绑定，先要解出各管道的流量分布，再在那个基础上求解温度的传播与混合。物理特性差异这么大，就很难用一套通用的求解器统一处理。
1.2 一个典型的区域综合能源系统拓扑长什么样
在动手之前，先把对象系统画清楚。我这次仿真的典型区域综合能源系统由三部分组成：一个IEEE 9节点电网（可以外购电，也可以由CHP发电），一个6节点天然气网（气源经过加压站、管道，供给CHP和燃气锅炉），一个8节点热网（包含若干供热管道、换热站和热负荷节点）。
这三块之间的连接点就是耦合设备的安装位置。CHP机组挂在电网3号节点、气网4号节点、热网5号节点上，燃气锅炉挂在气网5号节点和热网7号节点上，电锅炉挂在电网7号节点和热网6号节点上。这样三个网络通过耦合设备形成了一个闭环：气网给CHP供气，CHP发电并产热，发电送入电网、热量送入热网，电网再给电锅炉供电、电锅炉补充供热。
这个拓扑结构比较典型，涵盖了综合能源系统里最常见的三种耦合设备：CHP（气转电+热）、燃气锅炉（气转热）、电锅炉（电转热）。做代码实现的时候，这几个设备的模型就得分别写清楚，因为它们的工作特性完全不同。我后面的算例、代码结构也都是围绕这套拓扑展开的。
2. 电气热三网联立求解的数学建模：核心方程式逐一拆解
2.1 电网潮流的极坐标形式与节点类型处理
电网部分我采用的是最经典的新ton-Raphson极坐标潮流算法。对每个PQ节点，有两个待求的误差量：有功功率误差和无功功率误差；对每个PV节点，只有一个有功功率误差，电压幅值给定。方程组的核心表达就是功率平衡方程：
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))
这里的G_ij和B_ij对应支路导纳的实部和虚部，θ_ij是节点i和j的相角差。常规潮流里，平衡节点承担全网功率差额，典型做法是设一个外网等效节点，给定1.04的电压幅值和0度相角。
做综合能源系统能流计算时，电网部分有一个必须注意的改动：CHP机组不能简单当作固定的PQ节点或PV节点。因为CHP的发电量是气网供气量和自身热电比共同决定的结果，在气网没有收敛之前，你根本不知道CHP能发多少电。所以，在联合求解框架里CHP的并网节点一般当作P节点处理：有功功率由耦合迭代给出，电压幅值由无功调节维持，无功功率不做限制。这样处理的好处是电、气、热三个网络的变量可以被耦合方程统一调度，坏处是雅可比矩阵的结构和形成方式要跟着改，不能直接照搬标准潮流程序的节点分类逻辑。
2.2 天然气网稳态能流的Weymouth方程与节点压力计算
天然气网的稳态能流计算，核心是管道的Weymouth方程。高压输气管道中，气体的质量流量和管道两端压力的平方差成正比，这一点和电流与电压差的关系有些类似，但又不太一样。一般表达式写作：
F_ij = sign(π_i — π_j) * k_ij * sqrt(|π_i² — π_j²|)
其中，π_i和π_j是节点i、j的压力，k_ij是与管道内径、长度、压缩因子、气体性质有关的常数。节点流量平衡方程为：
Σ F_ij = L_i_gas
该式表示流入节点的天然气净流量等于该节点的负荷（包括天然气负荷和耦合设备耗气量）。气源节点是压力已知的边界节点，其他节点的压力是待求变量。
在Matlab实现里，这个非线性方程组的求解也使用牛顿法。雅可比矩阵的构造逻辑是：对每一条管道支路，求F_ij对π_i、π_j的偏导数；对每个节点，累加所有关联支路的偏导量。因为Weymouth方程的导数分母带有sqrt项，当节点压力差特别小或者接近零时，导数计算容易数值溢出。我后面会专门讲这个问题，这里是初学者最容易忽略的。
2.3 热网"水力-热力"双重耦合的方程体系
热网的建模是三个网络里最容易写乱的一个，因为变量最多、量纲也最容易出问题。热网的能流方程分为两个层次：
水力方程描述管道流量和水泵节点的压力平衡关系。其首要约束是节点流量连续性：流入某节点的体积流量等于流出该节点的流量，与供热负荷的流量需求相平衡。其次是环路压力方程：在闭合的供热环中，沿回路的压力损失之和为零，这一组方程将管道阻力（与流量平方和管径相关）和网络拓扑紧密耦合。
热力方程则关注温度场：每个节点的供水温度根据上游管道混合温度确定，节点处的热量平衡方程可写为：
Φ_i = C_p * m_i * (T_s,i — T_r,i)
其中，C_p是水的比热容，m_i是流经该节点的工质流量，T_s,i和T_r,i分别是该节点供回水温度。管道沿线的温度损失方程采用指数衰减模型：
T_end = T_amb + (T_start — T_amb) * exp(-λ * L / (C_p * m))
这里的λ是管道单位长度的热损失系数，L是管道长度。
水力方程的结果（即各管道流量分布）是热力方程求解的前提。所以在求解热网时，通常先解水力方程得到m的分布，再解热力方程确定温度分布。两者合在一起才算完成了热网的能流计算。这个"先水力、后热力"的顺序不能颠倒，否则流量未知，热量方程根本无从下手。
2.4 耦合设备建模：CHP、燃气锅炉、电锅炉的参数化处理
耦合设备是多能流计算的灵魂所在。我在代码里实现三种设备的模型，分别对应不同数学形式：
CHP机组采用典型的背压式机组模型，数学模型可以表达为：发电功率P_e和供热功率H_heat都与燃料消耗量F_gas呈线性关系。即P_e = η_e * F_gas * LHV_gas，H_heat = η_h * F_gas * LHV_gas。其中η_e是发电效率，η_h是热回收效率，LHV_gas是天然气低位热值。背压式CHP的热电比是固定的，这一特性使得电功率和热功率之间呈严格比例关系。在联立求解过程中，可根据系统的电、热负荷需求直接确定燃料消耗量，进而实现电网、气网和热网之间的信息交换与约束统一。
燃气锅炉模型是最简单的，就是一个转换效率η_gb，天然气输入F_gas，热量输出H_heat = η_gb * F_gas * LHV_gas。
电锅炉模型类似，只是输入变成了电能P_e，输出热量H_heat = η_eb * P_e。
这三个设备在Matlab里分别写成三个函数，输入输出都是结构体变量。这样迭代求解的结构就能统一：耦合设备从电网取/送出的电功率影响电网潮流方程；从气网取用的气功率影响气网流量平衡方程；向热网注入的热功率影响热网温度计算。设备模型越真实，迭代间的变量传递就越复杂，但原则是清晰统一的：每种耦合设备本质上都是三条网络之间的一条数据通道。
3. 求解策略怎么选：统一求解法与分解迭代法的取舍
3.1 两种主流框架的对比与适用场景
多能耦合能流的基本求解策略，学术界和工程上无非两条路：统一求解法和分解迭代法。
统一求解法将电网、气网、热网的所有方程放在一起，形成一个大规模非线性方程组，一次性用牛顿法迭代求解。好处是收敛性好，尤其是当耦合程度很强（比如CHP占系统发电比例很高）时，不容易出现发散。缺点是雅可比矩阵规模大、结构复杂，编程量大，且气网、热网的方程特性差异大，矩阵容易病态。更要命的问题在调试阶段：一旦不收敛，你很难判断问题出在电网方程、气网方程还是热网方程上。
分解迭代法的思路则是三个网络各自独立求解，耦合变量在外层做一个迭代循环。迭代流程是：先给定耦合设备的初始运行状态（比如假设CHP发电功率和热功率），分别求解电网、气网、热网的能流；然后根据三网的结算结果，更新耦合设备的运行状态；再重新求解三网。如此循环，直到前后两次迭代的耦合变量偏差小于收敛阈值。
我在实际项目里用的是分解迭代法改进版本。原因很直接：编程调试灵活，哪一层出问题一目了然；而且每条子网的求解函数可以单独测试，跟已知结果的单能流系统对比验证，大大降低了出错概率。缺点是单纯分解迭代在强耦合场景下可能迭代震荡或者收敛极慢。针对这个问题，我给耦合变量的更新加了松弛因子，效果非常明显：
X_new = X_old + α * (X_calculated — X_old)
这里α取0.5到0.8之间比较稳妥。α太大会导致震荡发散，α太小则收敛速度慢。当CHP的出力占比超过系统总出力的30%时，保守一点取α=0.5；占比较低时取0.8，可以兼顾速度和稳定性。
3.2 初值策略的选择：直接影响能否收敛的关键一环
初值设定是这个项目里我踩过最多坑的地方。电网潮流的初值相对简单，平启动即可：所有PQ节点电压幅值取1.0，相角取0。但气网的初值极为敏感：管道Weymouth方程里有sqrt(|π_i² — π_j²|)项，如果给定初值使得节点间压力差过大或者节点压力过低，迭代过程中平方差出现负值或者开方项数值异常，就直接挂了。
我的做法是：根据气源压力和工作压力范围，先预估一下整个气网的基准压力。假如气源节点压力是2.0MPa，末端节点压力大概在1.2MPa附近，那所有节点的初值统一取1.6MPa，别去费心思逐个节点估算。这样做的原理是，牛顿法的局部收敛特性要求初值不能偏离真解太远，取平均压力能保证初始压差在合理范围内，从而避免开方项出现非法值的概率。
热网的初值也类似。供水温度取设计值的上端比如100°C，回水温度取80°C，所有节点统一取这两个基准值。先解水力方程得到流量分布，再在流量基础上解热力方程，温度场的迭代比较温和。
判断收敛的标准我设为三网的最大不平衡量小于1e-6（标幺值）。注意是三个网络的不平衡量都同时满足，不是只看某一个网络的收敛情况。因为耦合系统的特点就是各网络相互牵制，只收敛其中一个没有意义。
3.3 解耦迭代的收敛判据与松弛更新逻辑
外层耦合迭代的流程，我写成下面的伪代码逻辑：
- 初始化：假设CHP的发电功率P_e、热功率H_heat，燃气锅炉的供热量，电锅炉的耗电功率
- 电网潮流求解：把CHP作为P节点（有功给定），电锅炉作为负荷节点，求解电网模型，得到各节点电压、相角和各耦合节点的结算结果
- 气网能流求解：根据CHP和燃气锅炉的耗气量（由第一步的初始假设或上一次迭代结果折算）作为气网负荷，求解气网模型，得到修正后的节点压力和可用供气量
- 热网能流求解：根据CHP、燃气锅炉、电锅炉的热功率注入与热负荷需求，求解热网模型，得到供热温度和流量
- 根据气网结算得到的流量，反推CHP的发电功率上限并修订热功率值；再根据热网的实际供热需求，修正CHP和电锅炉的分配
- 判断修正量与上次迭代的偏差是否小于收敛阈值，不满足则返回步骤2继续迭代
- 所有偏差满足阈值后，输出全局能流结果
这套流程写出来看着简单，实际操作中容易忽略一个点：步骤5的耦合变量更新要带上松弛因子，而且每个设备的更新松弛因子最好单独设置。比如CHP的更新比较稳定，α可取0.8；电锅炉因为受电网和热网双重约束，振荡风险高，α取0.5。
4. Matlab代码实现全过程：从函数架构到算例验证
4.1 全套代码的文件结构与核心数据流设计
下面是这次项目全套Matlab代码的文件组织方式，每个文件的功能都做了拆分，互不干扰：
system_main.m 主程序入口，控制整体迭代流程 load_system_data.m 读入电、气、热三网的拓扑参数与负荷数据 build_ybus.m 根据支路参数构建电网节点导纳矩阵 electric_powerflow.m 电网潮流求解（牛顿-拉夫逊法） gas_powerflow.m 天然气网能流求解（牛顿法） hydraulic_solve.m 热网水力平衡方程求解 thermal_solve.m 热网热力方程求解 coupling_update.m 耦合设备变量更新与松弛处理 plot_reports.m 三网能流结果的可视化输出在设计数据流的时候，我建议所有中间变量和结果都用结构体变量按网络分类封装。三个大结构体sys_elec、sys_gas、sys_heat分别存放各网络物理模型所需的数据，比如节点参数、支路参数、求解结果。耦合设备的参数单独放在struct_coupling结构体里。
这个架构的好处是，每段代码的输入和输出都是清晰的数据结构，函数之间的耦合度被降到最低。比如electric_powerflow函数只接收sys_elec结构体和耦合节点的注入功率，只返回更新后的电压和相角结果，它不需要关心气网压力是多少、热网供水温度是多少。模块解耦的好处，在问题排查阶段优势尤其明显：电网不收敛时只需要检查电网这部分的网络连接和参数配置就行，不会被气网和热网干扰。
4.2 核心求解函数关键技术点：牛顿法雅可比矩阵的组装
先看电网潮流求解的核心循环。标准牛顿法潮流的核心就是构造雅可比矩阵，并反复求解增量方程：
for iter = 1:iter_max % 计算功率不平衡量 [dP, dQ] = compute_power_mismatch(V, theta, Ybus, S_spec); mismatch = [dP; dQ]; if max(abs(mismatch)) < tol break; end % 构造雅可比矩阵 J = build_jacobian(V, theta, Ybus); % 求解增量 dx = J \ mismatch; % 更新状态变量 theta = theta + dx(1:n_node); V = V + dx(n_node+1:end) .* V; end这段代码看起来简单，但有几个细节需要强调。雅可比矩阵的组装，公式推导可以参考任何一本电力系统分析教材，但在Matlab里有一个程序效率的大坑：如果直接用两个for循环逐个元素填充雅可比矩阵，三节点系统跑起来没问题，一旦节点数超过30，每一次迭代都会因为矩阵组装太慢让人等到怀疑人生。更好的做法是利用稀疏矩阵的索引批量赋值，或者直接使用vectorized方式按节点对数组同时计算。我实测过，同样的IEEE 9节点算例，向量化后的代码比双重循环快至少一个数量级，这个差距在蒙特卡洛批量计算或者考虑多场景规划时会体现得非常明显。
再来看看气网的求解核心。气网求解的难点其实不在牛顿法的框架，而在于Weymouth方程导数的正确处理：
for iter = 1:iter_max % 计算节点流量不平衡量（包括耦合设备负荷） node_mismatch = compute_gas_mismatch(Pressure, gas_loads, pipe_k, pipe_connect); J = build_gas_jacobian(Pressure, pipe_k, pipe_connect); deltaP = J \ node_mismatch; Pressure = Pressure + deltaP; if max(abs(node_mismatch)) < tol break; end end这里build_gas_jacobian函数有一个关键逻辑：对于管道ij，流量F_ij对节点压力P_i的偏导要分为两种情况计算。当P_i大于P_j时，导数为正，值为k_ij * P_i / sqrt(|P_i² — P_j²|)。当P_i小于P_j时，导数为负。由于开方项位于分母中，当压力差趋近于零时，雅可比矩阵元素会趋向无穷大，该数值特性极易导致迭代不收敛。为避免此问题，需在程序中对压力差设置最小限值，例如当|P_i² — P_j²|小于某个阈值（我取1e-6）时，将该导数按上限值截断处理。曾经由于忽略该细节导致程序反复发散，排查了很久才意识到问题所在。
4.3 一个典型算例的测试结果与误差分析
我用的算例配置如下：CHP安装在三网交汇的关键节点，电功率输出为60MW，热功率输出为50MW；燃气锅炉的热功率为30MW，电锅炉的热功率为20MW；电网侧的总电负荷约为110MW，外网联络线提供剩余约30MW功率；气源节点压力设定为2.0MPa，两个气网负荷节点的天然气流量需求分别约为4000和2500 m³/h。
经过分解迭代求解，整个系统在14次外层迭代后收敛（松弛因子取0.6），各网络的不平衡量均达到10⁻⁶级别。电网的平衡节点有功出力为-28.5MW（即从外网购电28.5MW），各节点电压幅值均在0.98到1.03之间，满足常规设计要求。气网的节点压力从气源2.0MPa下降到末端约1.35MPa，压降符合管道长度和流量的物理预期。热网的供水温度经过管网衰减后，最远负荷节点比热源出口下降了约9°C，回水温度分布也比较均匀。
为了验证程序的正确性，我做了两组对照实验。第一组是单一电网潮流对照：把CHP当作固定出力节点，运行电网潮流程序，与Matpower的计算结果对比，节点电压偏差小于1e-6，说明电网模块实现无误。第二组是气网对照：把已知流量和压力参数的6节点气网用商业管网仿真软件校核，结果压力偏差在2%以内。这样逐模块验证过的程序，再去做三网联算，才敢说结果是可信的。
5. 实战中一定会遇到的坑：三个典型问题与排查实录
5.1 气网开方项异常造成不收敛的排查过程
项目最早期调试阶段，气网部分一直报出NaN，程序直接崩溃。当时把错误定位到sqrt处，但百思不得其解，因为从物理概念上，节点压力平方差为负似乎不该出现。后来通过打印每次迭代的节点压力值发现，迭代第一次就出现某个节点压力值跳到接近零，随后它与相邻节点间的压力平方差直接变负，开方返回复数。
根因其实是两个问题叠加。第一个是初值给得太差，我把节点初值都设成0.5MPa，但气源压力是2.0MPa，迭代前期压力调整跨度大，部分节点被牛顿法推到了不合理区域；第二个是没有在Weymouth方程里加上数值保护，当开方根内部出现极小负数时，直接就得到NaN了。修正方案是两个方向同时入手：初值统一取工作压力范围的中间值，同时在Weymouth公式里加入限幅保护，用max(eps, abs(p_i²—p_j²))来保险。
5.2 耦合变量振荡不收敛，为什么需要逐设备设置松弛因子
框架刚跑起来的时候，出现了比较有意思的现象：电压和压力都稳定，但CHP的出力在相邻两次迭代之间来回跳，幅度还不小（±8MW左右），导致外层循环无法收敛。直觉解决方案是把全局松弛因子调小，发现效果有限，因为电网、气网模块本身是精确计算的，真正振荡的来源是耦合变量更新逻辑——我先更新CHP发电量，再根据它来更新耗气量，而耗气量又马上被气网结算结果修正，形成正反馈回路。
后来我把CHP、燃气锅炉、电锅炉的松弛因子分别独立设置：CHP为0.6，燃气锅炉为0.5，电锅炉为0.4。因为电锅炉同时受电网和热网两侧约束，其更新最容易引起波动。结果外层迭代从完全不收敛变成14次收敛，效果立竿见影。这个教训让我意识到，在多能流迭代里，不要贪图代码简单就用一个统一的松弛因子，不同耦合设备的时间常数和灵敏性完全不同，分开设置才是合理做法。
5.3 单位换算的"隐蔽炸弹"：天然气标方与热功率的换算
最后一个要提醒的坑是单位。这是我检查看了半天才发现的低级错误，但后果很严重。天然气网里，气负荷的单位习惯用m³/h（标方），而电气热联立方程里，CHP的耗气量需要换算成MW热功率才能和电网功率平衡方程式统一量纲。
天然气低位热值LHV一般取36 MJ/m³，换算关系是1 m³/h天然气约等于0.01 MW。计算方法是：36 MJ/m³乘以1000除以3600秒，结果正好约等于10 kW，也就是0.01 MW。这个数值关系看起来很直白，但如果代码里直接用一个常数去乘除，很容易在某处漏掉或者写错一个数量级。更隐蔽的是锅炉效率、CHP效率这些百分比数，在单位换算时也要参与运算，一旦顺序弄错，结果就可能出现CHP耗气量比气网总供气量还大的离谱局面。
我的建议是，在所有计算之前，定义一组统一的基准值和转换函数，把天然气质量流量、标准体积流量、热功率三套单位全部转换成统一的标幺值或者统一的MW制后再进入迭代计算。程序内部只用一套单位制，只在输入和输出层做转换。
6. 一些个人经验和可以继续扩展的方向
做这个项目最大的体会是：多能耦合能流计算，"耦合"两个字才是重点。不要把它当成三个单能流程序的简单拼接，而是要在求解框架设计时就充分考虑耦合变量的传递机制。初期我犯的错就是把80%的精力花在三个网络各自的算法实现上，最后联调阶段才发现真正的难点在更新与收敛策略上。
另外一个实用的经验是，Matlab里做这类项目，一定要充分利用结构体和稀疏矩阵。结构体让三网的参数管理变得清晰，而稀疏矩阵让大规模节点的求解效率成倍提升。如果只是为了"能跑通"，用全矩阵也凑合，但如果你以后想把算例规模扩大，或者在项目基础上做优化调度、多场景分析，全矩阵的性能瓶颈就会非常明显。
这套代码后续可以扩展的方向还挺多：一是把稳态能流扩展为准稳态时序仿真，模拟一天24小时或者四季的更替场景；二是在能流计算基础上叠加经济优化目标，比如以运行成本最小为目标函数调整CHP出力比；三是加入储能设备模型，让气网、热网具备调节能力，看看储能对多能互补系统能流分布的影响。如果你已经在做或者准备做综合能源系统的仿真分析，欢迎一起交流具体实现中的细节。
