1. 为什么区域综合能源系统的潮流计算绕不开多能耦合这几年做综合能源系统方向的朋友应该都有同感单独算电网潮流、天然气网水力分析或者热网热力工况的文章和代码已经非常成熟但一旦把三者放在同一个系统里联立求解问题立刻变得微妙起来。原因很简单——电、气、热三个网络并不是各自封闭运行的它们通过燃气轮机、电锅炉、热电联产机组、电制冷机、吸收式制冷机等耦合设备互相咬合。电网的负荷变化会影响燃气机组的出力进而改变天然气网的流量天然气网的压力和流量约束反过来又限制发电机的进气量最终影响电网的节点电压和功率分布。热网那边同样是这个逻辑热电联产机组的电出力与热出力之间存在运行区间约束热负荷波动会牵制电出力而电锅炉的投切又同时改变电网和热网两侧的平衡状态。正是这种强耦合特性让先分别算好再拼接的传统思路在综合能源系统里失去了可解释性。也正因为如此计及多能耦合的电气热能流计算才成为区域综合能源系统规划、运行、调度乃至可靠性分析的基础工具。这篇博文想分享的是一套基于Matlab实现的区域综合能源系统电气热能流计算代码重点讲清楚建模思路、耦合机制、解算方法以及我在调试过程中踩过的坑方便打算做这个方向或者正在做相关课题的同学直接参考甚至复用。按惯例先交代一下这套代码的适用范围它面向的是包含电力网络、天然气网络和热力网络的区域级综合能源系统耦合设备涵盖热电联产机组、燃气锅炉、电锅炉、电制冷机等常见类型。解算目标是求取稳态工况下的电网节点电压幅值与相角、气网节点压力与流量、热网节点供回水温度及管道流量。程序的整体思路是采用统一求解法又称全同时求解法即把电、气、热三个网络的方程连同耦合设备的约束一起组成一个大规模非线性方程组用牛顿-拉夫逊法迭代求解这也是目前多能流计算文献里精度和收敛性最稳妥的主流方案。2. 三个子网络单独建模时最容易忽略的物理边界2.1 电网部分其实就是带耦合项的标准交流潮流电网子模型沿用的是经典交流潮流模型。对于节点 i有功和无功功率平衡方程写为[ P_{i} U_{i}\sum_{j \in i}U_{j}(G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ][ Q_{i} U_{i}\sum_{j \in i}U_{j}(G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]这里的 (P_{i}) 和 (Q_{i}) 不是简单的负荷值而是要写成电源注入减去负荷消耗再减去耦合设备消耗的形式。例如节点上连着燃气轮机和电锅炉时[ P_{i} P_{G,i} P_{CHP,i} - P_{L,i} - P_{EB,i} ][ Q_{i} Q_{G,i} Q_{CHP,i} - Q_{L,i} - Q_{EB,i} ]其中下标 (G) 表示常规发电机组(CHP) 表示热电联产机组(EB) 表示电锅炉家 (L) 表示常规电负荷。这里有个细节很多初学者会漏掉电锅炉、电制冷机这类纯电耗设备通常以恒定有功功率形式建模但它们的无功功率往往被忽略。实际上电锅炉通过电力电子变换装置接入电网时会有一定的无功消耗如果系统规模不大、耦合设备占比高忽略无功可能会让迭代结果出现偏差。我在代码里把这类设备的无功按功率因数 0.95 估算虽然严格说要根据实际设备参数来定但至少比直接设零更合理。2.2 天然气网的管道方程不能只记一个公式天然气网络的稳态模型核心是管道流量方程。对于连接节点 (m) 和 (n) 的管道标准形式的稳态流量公式是[ f_{mn} C_{mn} \cdot s_{mn} \cdot \sqrt{s_{mn} \cdot (\pi_m^2 - \pi_n^2)} ]其中 (C_{mn}) 是管道常数由管径、长度、摩擦系数、气体温度、压缩因子等共同决定(\pi_{m}) 和 (\pi_{n}) 是节点压力符号系数 (s_{mn}) 的取值为[ s_{mn} \begin{cases} 1 \pi_m \geq \pi_n \ -1 \pi_m \pi_n \end{cases} ]这个式子的本质是压差平方与流量成正比而不是压差本身与流量成正比。实际编程中有几个地方要注意。第一如果直接用 (\pi_{m}^{2}) 和 (\pi_{n}^{2}) 作为状态变量方程形式会简单很多牛顿法的雅可比矩阵也更好求。第二管道常数 (C_{mn}) 的单位和数值范围很关键不同文献里有的用标准立方米每天有的用千克每秒混用会导致流量偏差几个数量级。第三天然气网中节点分为压力已知的平衡节点和压力待求的一般节点平衡节点通常设在气源或管网上游的门站这和电网中的平衡节点Vθ节点逻辑上是对应的。气网节点流量平衡方程按下式写[ \sum_{n \in m}f_{mn} G_{s,m} - G_{L,m} - G_{CHP,m} - G_{GB,m} 0 ]其中 (G_{s,m}) 是气源注入流量(G_{L,m}) 是常规气负荷(G_{CHP,m}) 和 (G_{GB,m}) 分别是热电联产机组和燃气锅炉消耗的天然气流量。这里要注意燃气机组消耗的天然气流量在气网侧是负荷但在电网侧是电源这种角色转换正是多能耦合在方程层面的体现。2.3 热网的质调节与量调节对应不同的方程形态热力网络的稳态模型比电网和气网都要复杂一些关键在于它涉及两个物理过程水力工况和热力工况。水力工况描述热水在管道中的流量分配由节点流量平衡和环路压降约束决定热力工况描述温度沿管道的分布以及节点供回水温度的混合。在实际工程中热网绝大多数时间运行在质调节模式即各管道流量保持不变通过调节供回水温度来适应热负荷变化。这种模式下可以先把水力工况算好甚至在程序中当成已知量输入只对热力工况迭代求解计算量会大幅下降。但如果是量调节模式流量本身也在变就必须联立求解水力与热力方程收敛难度会明显上升。我的代码默认采用质调节模式热力工况用节点法建模。对于热力管网节点温度分为供水温度和回水温度。热源节点处供水温度为设定值负荷节点处回水温度为设定值其他节点则满足流量加权温度混合方程[ \sum_{k \in S_{m}}(\dot{m}{k}c{p}T_{k}^{out}) T_{m}^{mix}c_{p}\sum_{k \in S_{m}}\dot{m}_{k} ]式子表示若干条汇入同一节点的管道热水混合后的温度等于节点混合温度其中 (c_{p}) 是水的比热容(\dot{m}_{k}) 是各支路的质量流量。管道出口温度则按热损失方程衰减[ T_{j}^{out} (T_{j}^{in} - T_{a}) \cdot \exp\left(-\frac{\lambda L}{\dot{m} c_{p}}\right) T_{a} ]这里 (\lambda) 是管道单位长度传热系数(L) 是管长(T_{a}) 是环境温度。如果管道保温效果好或者系统规模不大也可以暂时忽略热损失但实际工程算出来的温度衰减对后续设备效率计算是有影响的所以代码里保留了这个项方便后续扩展到动态热网模型。热力节点方程最终要满足两类约束每个节点的功率平衡供热量等于负荷侧吸收的热量和各环路的压降为零如果不做水力计算则可以跳过。质调节模式下给定总循环流量和节点热负荷每个节点的供回水温度就是待求状态变量负荷节点的回水温度由热负荷和流量共同决定[ T_{r,i} T_{s,i} - \frac{\Phi_{L,i}}{c_{p}\dot{m}_{i}} ]其中 (\Phi_{L,i}) 是该节点的热负荷功率。这个式子简单但它是热网与电网耦合的另一个入口——如果该节点的热量来自热电联产机组热源侧还要再建一个热出力与电出力对应的耦合方程。3. 耦合设备模型整个计算中最容易出错的环节3.1 热电联产机组的热电比约束热电联产CHP机组是综合能源系统里最核心的耦合设备它的运行特性通常用热电比约束来描述。最常见的是定热电比模型[ \Phi_{CHP} c_{m} \cdot P_{CHP} ]其中 (c_{m}) 是热电比。定热电比模型虽然简单但在实际系统里过于理想化——真实的热电联产机组的工作区间是一个多边形电出力和热出力之间存在可行域边界。更精确的建模方式是用线性可行域来约束典型的一组约束写为[ P_{CHP,min} \leq P_{CHP} \leq P_{CHP,max} ][ 0 \leq \Phi_{CHP} \leq \Phi_{CHP,max} ][ P_{CHP} \geq a_{1}\Phi_{CHP} b_{1} ][ P_{CHP} \leq a_{2}\Phi_{CHP} b_{2} ]这四条约束合在一起刻画了背压式和抽汽式机组的不同运行区间。在潮流计算里处理这种不等式约束的方式通常是先不管它迭代收敛后检查是否越界如果越界则把耦合设备出力定格在边界值上再重新迭代一两次。后面我会专门讲这个处理逻辑。3.2 燃气轮机和燃气锅炉的气耗特性燃气轮机消耗的天然气和电出力之间的关系可以表示为[ G_{GT} \frac{P_{GT}}{\eta_{GT} \cdot LHV_{gas}} ](\eta_{GT}) 是燃气轮机的发电效率(LHV_{gas}) 是天然气低位热值通常取 9.7 kWh/m³ 左右。实际计算中要注意单位统一问题——如果电网侧功率单位是 MW气网侧流量单位是 m³/h那么低位热值必须统一成 MWh/m³ 才能保证等式两边的量纲一致。我最初写代码时在这里栽过跟头算出来的气耗比文献值大了一倍查了半天最后发现是热值单位换算出了错。这个看似不起眼的问题其实是很多同学在做多能流计算时收敛错误甚至结果荒诞的根源。燃气锅炉的气耗模型类似[ G_{GB} \frac{\Phi_{GB}}{\eta_{GB} \cdot LHV_{gas}} ]唯一不同的是输入是热出力 (\Phi_{GB}) 而不是电出力。燃气锅炉可以作为分布式热源配置在热网的任意节点上在热负荷高峰时投入。3.3 电锅炉与电制冷机电转热/冷的单向耦合电锅炉和电制冷机的拓扑逻辑相对简单它们只体现电 → 热/冷的单向转换不存在热或冷回馈到电网的通路。电锅炉的模型就是一个效率系数[ \Phi_{EB} \eta_{EB} \cdot P_{EB} ]电制冷机的模型是制冷系数COP[ Q_{cool} COP_{EC} \cdot P_{EC} ]这两个设备本身没有额外状态变量但在电网方程里要计入它们的电功率消耗在热网/冷网方程里要作为热源/冷源来处理。从编程角度说它们的耦合强度不如 CHP 和燃气轮机但数量多、分布广在整体雅可比矩阵里依然会形成非零交叉块不能省略。4. 统一求解法如何把电、气、热三个方程塞进同一个牛顿迭代框架4.1 状态变量统一排队与初值选取策略统一求解法的核心思想很简单把电网、气网、热网的全部未知量放在同一个状态向量里所有网络方程和耦合方程构成同一个残差向量然后整体用牛顿法迭代。以我实现的代码为例状态变量与残差方程的数量和排列决定了雅可比矩阵的结构也直接决定了程序的通用性和扩展难度。我采用的状态向量排列顺序为电网部分PQ 节点的电压幅值 (U) 和相角 (\theta)PV 节点的相角 (\theta)不包括平衡节点其电压幅值和相角已知气网部分所有非平衡气节点的压力平方 (\pi^{2})热网部分所有非给定温度节点的供水和回水温度质调节模式下对应的残差方程排列顺序与状态变量一一对应这样雅可比矩阵的构造最简单调试时也便于定位问题。整个系统的雅可比矩阵从结构上看是一个分块矩阵[ J \begin{bmatrix} J_{EE} J_{EG} J_{EH} \ J_{GE} J_{GG} J_{GH} \ J_{HE} J_{HG} J_{HH} \end{bmatrix} ]其中 (J_{EE})、(J_{GG})、(J_{HH}) 是三个网络各自的自雅可比块非对角块 (J_{EG})、(J_{GH}) 等是网络间的耦合雅可比块。如果忽略多能耦合这些非对角块全为零问题退化为三个独立网络分别求潮流只有耦合设备存在时非对角块才真正不为零这也是多能耦合在数学上的本质。初值选取方面电网部分采用平启动flat start——PQ节点电压幅值取 1.0 p.u.相角取 0气网节点压力初值取气源压力的 80% 左右热网供回水温度初值取设计工况的典型值。平启动对牛顿法而言是最常见的策略但对于强耦合的大规模系统一旦初值离解太远迭代很容易发散。我的经验是给雅可比矩阵加一个对角阻尼因子——在迭代初期把阻尼因子设得大一些比如 0.5随着残差下降逐步减小到 1.0这套阻尼牛顿法在综合能源系统多能流计算里表现相当稳定。4.2 雅可比矩阵分块计算的程序实现思路要在 Matlab 里高效实现分块雅可比矩阵最高效的做法是分别对电、气、热三个子网络模块写各自的解析雅可比函数再由主程序把分块矩阵拼装起来。每个子模块的输入都是整个状态向量和该子网络自身的节点数据输出是该子网络对应残差对该子网络状态变量的偏导块。耦合块的偏导藏在耦合设备方程里需要针对每一类耦合设备单独写。以 CHP 机组为例假设它连接电网的节点 i 和气网的节点 m同时向热网节点 k 供热。它引入的偏导关系包括[ \frac{\partial P_{i}}{\partial \pi_{m}^2}, \quad \frac{\partial G_{m}}{\partial P_{i}}, \quad \frac{\partial \Phi_{k}}{\partial P_{i}} ]这三条链式关系会让整个雅可比矩阵出现三个非零的交叉块。具体数值计算时如果用解析法需要根据耦合方程组逐项求偏导这对代码结构的模块化要求比较高如果嫌麻烦也可以用数值差分近似这些交叉项即对每个耦合变量施加一个微小扰动观察各网络残差的变化量来估计偏导值。数值差分的精度和步长选择比较敏感但胜在实现简单、不易出错。我的代码里提供了解析法和数值差分法两种模式默认打开解析法方便对照验证。在 Matlab 的实现层面程序主体结构可以这样组织% 主函数入口IES_PowerFlow.m % 输入电网数据、气网数据、热网数据、耦合设备数据 % 输出各网络状态变量、收敛信息 state initState(elec, gas, heat); % 状态变量初始化 residual computeResidual(state, elec, gas, heat, coupling); J assembleJacobian(state, elec, gas, heat, coupling); iter 0; while norm(residual, inf) tol iter maxIter delta -J \ residual; % 阻尼因子 lambda 1.0; while norm(computeResidual(state lambda * delta, ...), inf) norm(residual, inf) lambda 0.5 * lambda; if lambda 1e-4, break; end end state state lambda * delta; residual computeResidual(state, elec, gas, heat, coupling); J assembleJacobian(state, elec, gas, heat, coupling); iter iter 1; end这个流程和常规电力潮流计算的牛顿法是同构的多出来的部分只有残差函数和雅可比函数需要同时涵盖气网和热网的方程。程序的骨架清晰之后真正麻烦的其实在于耦合方程接入残差向量时必须保证状态变量与方程的数目严格相等——多一个未知量而少一个方程矩阵就奇异迭代必然发散。这是统一求解法最考验工程细节的地方。4.3 迭代收敛判据与常见发散原因排查收敛判据我采用的是无穷范数[ \max\left( \max_i |\Delta U_i|, \max_i |\Delta\theta_i|, \max_i |\Delta\pi_i^2|, \max_i |\Delta T_{s,i}|, \max_i |\Delta T_{r,i}| \right) \varepsilon ]这个判据比单独看残差要严格因为残差小并不一定代表状态变量修正量小——在病态系统里两者可能会出现数量级不一致的情况。建议同时监控两个指标修正量无穷范数和残差无穷范数以修正量为主判据残差作为辅助参考。收敛容差设置在 (10^{-6}) 左右比较合理太严格会拖慢速度太宽松结果不够准。在调试过程中发散的情况大多可以归结为以下几类原因第一变量-方程数量不匹配。这是新手最容易犯的错。电网部分除了平衡节点外每个节点都有一个有功方程PQ 节点还有一个无功方程气网每个非平衡节点有一个流量平衡方程热网每个节点有一个温度混合方程。算清楚这些之后再和状态变量核对确保一一对应。第二初值偏离实际解太远。气网压力初值尤其敏感因为管道流量方程里有 (\pi^{2}) 项如果初始压力设得太低迭代过程中平方项产生巨大的梯度变化会导致震荡发散。建议气源压力较高的节点附近初值取最高压力的 80% 以上。第三耦合设备的效率参数或热电比设置不合理。比如 CHP 效率设为 0.9再叠加燃气轮机效率 0.4气耗算出来的结果就会让气网负荷暴涨导致气网压力跌破下限整个迭代过程不收敛。这不是算法问题而是模型参数问题排查时要先看物理上合不合理。5. 一个典型算例从仿真结果看耦合效应的实际影响5.1 算例配置与边界条件用一个 6 节点电网、6 节点气网、6 节点热网的小型区域综合能源系统来演示程序效果。电网包含 1 个平衡节点、2 个 PV 节点常规发电机组和 3 个 PQ 节点。气网包含 1 个气源节点和 5 个负荷节点。热网包含 1 个热源节点由 CHP 供热的换热站和 5 个热负荷节点。耦合设备配置为1 台 CHP 机组连接电节点 2、气节点 3、热节点 1、1 台燃气锅炉连接气节点 5 和热节点 4、1 台电锅炉连接电节点 5 和热节点 5。各网络的关键参数如下表所示参数数值电网基准容量100 MVA气源节点压力4.0 MPa对应平方 16 MPa²热网供水温度110 ℃热网回水温度设计值60 ℃CHP 热电比1.2CHP 发电效率0.42燃气锅炉效率0.9电锅炉效率0.98负荷水平方面电网总负荷设定为 150 MW气网总气负荷设为 80 m³/h换算成热值约 776 MW 等价值热网总热负荷为 60 MW。CHP 承担其中 40 MW 的供热燃气锅炉承担 15 MW电锅炉承担 5 MW。5.2 迭代收敛过程程序从平启动初值开始迭代逐次残差变化如下迭代次数修正量无穷范数残差无穷范数01.0初始36.7210.42319.8720.08752.1430.00690.17240.00020.00653.8×10⁻⁶1.2×10⁻⁴62.3×10⁻⁸7.8×10⁻⁷第 6 次迭代达到收敛迭代次数在牛顿法的正常范围内。这个收敛模式说明系统参数设置合理统一求解法的效率与纯电力潮流相当没有出现因耦合引入导致的明显迭代次数增加。5.3 耦合效应量化对比单网络解耦 vs 统一求解为了直观展示多能耦合对计算结果的实际影响我把代码里耦合设备的交互全部断开相当于把所有耦合方程替换为恒定功率注入分别跑三个独立的常规潮流再将结果与统一求解法做对比。差异最明显的以下几个量状态量解耦计算统一求解偏差CHP 所在电节点电压幅值 (p.u.)1.0451.0212.3%电网平衡节点有功出力 (MW)71.568.23.3 MW气源节点输出流量 (m³/h)105.496.78.7 m³/h热网负荷节点回水温度 (℃)62.365.83.5 ℃这个对比说明解耦计算在很多情况下虽然不至于得到完全荒谬的结果但偏差已经足以影响工程决策。尤其是当 CHP 热出力受到热负荷变化牵制时电网侧的有功平衡会受到显著影响——如果按解耦结果做调度电网平衡节点的出力安排会偏高气网气源采购量也会产生非忽略的偏差。对运行优化或者可靠性分析来说这种偏差积累到最后可能让结论反向。从计算效率上看解耦计算确实比统一求解快那么一点少迭代一两次但考虑到底层物理的完整性这个成本是完全值得的。这也是为什么当前多能流领域的主流文献几乎都推荐统一求解法。6. 代码架构设计模块化拆分比一个脚本全塞进去靠谱得多6.1 文件结构安排与接口约定很多同学拿到一个题目就习惯把所有逻辑写在一个大脚本里跑通这样在几十行的小程序里没有问题但一旦系统规模拉大或者需要换算例代码的可维护性会快速恶化。这套代码做成了模块化结构每个网络一个接口文件、一个数据定义文件和一个残差/雅可比函数文件核心文件结构如下IES_PowerFlow/ ├── main_IES_PF.m % 主程序数据加载、初始化、迭代 ├── data/ % 数据目录 │ ├── case6_elec.m % 6节点电网数据 │ ├── case6_gas.m % 6节点气网数据 │ └── case6_heat.m % 6节点热网数据 ├── elec/ │ ├── elec_residual.m % 电网残差 │ └── elec_jacobian.m % 电网雅可比 ├── gas/ │ ├── gas_residual.m % 气网残差 │ └── gas_jacobian.m % 气网雅可比 ├── heat/ │ ├── heat_residual.m % 热网残差 │ └── heat_jacobian.m % 热网雅可比 └── coupling/ ├── coupling_model.m % 耦合设备参数读取与状态映射 ├── coupling_residual.m % 耦合方程残差 └── coupling_jacobian.m % 耦合方程雅可比接口约定上每个子网络函数都采用统一的输入输出格式。输入参数是完整的状态向量 state 以及该子网络自身的结构体数据 data输出是残差向量及其雅可比。这样主程序不需要知道每个函数内部是怎么实现的只要按约定调用即可。想换成别的算例时只需要在 data 目录下新增一组数据文件或者把现有算例的节点和管道参数改掉。6.2 统一状态向量与数据索引防止变量错位的关键统一求解法的难点在于把三个网络的状态变量放在同一个向量里时如何保证每个残差方程对应的偏导存储在雅可比矩阵的正确位置。我的做法是维护一张状态变量-索引映射表在主程序中预先定义idx.U 1:n_elec_PQ; % PQ节点电压幅值索引 idx.Theta idx.U(end)1:idx.U(end)n_elec_PVn_elec_PQ; % 相角索引 idx.Pi idx.Theta(end)1:idx.Theta(end)n_gas_slackless; % 气网压力平方索引 idx.Ts idx.Pi(end)1:idx.Pi(end)n_heat_unknown; % 供水温度索引 idx.Tr idx.Ts(end)1:idx.Ts(end)n_heat_unknown; % 回水温度索引在残差函数里每个子网络函数只需要按照映射表取自己需要的状态变量子集计算完残差后输出到与映射表对应的位置即可。雅可比矩阵的组装则通过一个稀疏矩阵来累积每个子网络函数返回的是它自己的那部分行和列的稀疏矩阵主程序累加到一起。稀疏矩阵的重要性怎么强调都不过分。如果按全稠密矩阵存储6节点系统还能忍受但扩展到几十几百个节点时矩阵规模迅速膨胀存储和求解效率都会崩。Matlab 里用 sparse 命令把雅可比矩阵声明为稀疏格式配合反斜杠求解计算速度可以提高一两个数量级。6.3 网络拓扑编号与物理连接关系的映射计算中另一个需要谨慎处理的是电网、气网、热网三套编号体系之间的对应关系。电网的节点编号是 1 到 6气网和热网各自也有自己的一套编号它们之间唯一的联系是耦合设备表——每台耦合设备记录了自己挂在电网哪个节点、气网哪个节点、热网哪个节点。程序运行时耦合设备模型要根据这张设备表把不同网络的变量提取出来再组装成耦合方程。这里我采用的是一种面向对象的思路虽然 Matlab 的结构体数组也能实现但用 classdef 定义耦合设备类会让代码可读性更好。每个耦合设备实例包含自己的类型、效率参数、接入节点信息、状态变量索引。在组装雅可比时主程序遍历所有耦合设备实例每个实例把自己对应的耦合残差和耦合雅可比块填入全局矩阵的相应位置。7. 那段改了一个参数整个系统发散的排查经历7.1 问题现象明明是同一套程序算例一换就崩代码写完之后我先拿 6 节点系统调通了一切正常。后来想验证通用性把电网换成一个 33 节点的配电网络气网和热网也在对应规模上做了扩充结果一跑就发散——残差在前几步一直下降到了第 4 次迭代突然往上跳之后直接飞掉。我首先怀疑的是初值问题把气网初值调到最高压力的 90%热网温度初值也做了调整但依然发散。随后我逐个检查状态变量和残差的对应关系确认数量没错又用数值差分法对比了解析雅可比两者的结果在单个子网络上是一致的。7.2 排查链路从变量数量检查到耦合块数值病态真正发现问题是在检查雅可比矩阵条件数的时候。我用 condest 估算了整体雅可比矩阵的条件数发现已经超过 (10^{15})几乎处于奇异状态。进一步拆解每个分块发现电网部分的自雅可比块条件数正常气网部分的条件数在 (10^{7}) 左右也还可以接受。问题出在耦合块和热网自块上——热网部分的自雅可比块条件数高达 (10^{12})因为 33 节点系统里热网管道长度差异过大最短的管道和最长的管道差了三个数量级导致温度衰减指数项差异极大雅可比矩阵中的元素出现了严重的量级失衡。7.3 解决方案量纲归一化与支路参数预处理解决办法分两步。第一步是对热网管道参数做预处理把传热系数和管长合并成无量纲的温度衰减因子[ f_{loss} \exp\left(-\frac{\lambda L}{\dot{m} c_{p}}\right) ]这个因子本身就是 0 到 1 之间的数直接用它替代原方程里的指数表达式避免指数项在变量微扰下产生异常敏感的变化。第二步是对整个热网做量纲归一化处理把温度全部表示为相对于环境温度的归一化值使所有方程里涉及的变量都在 0.1 到 10 这个量级范围内。修改之后同一套程序跑 33 节点系统顺利收敛迭代次数 8 次残差降到 (10^{-7}) 以下。这个经历给我最大的启发是多能流计算的收敛性问题很多时候不是算法本身的缺陷而是数值问题的表现。物理模型没问题方程推导没问题但数值尺度不匹配时再好的算法也会失效。做这类程序一定要养成检查雅可比矩阵条件数的习惯以及在做大规模算例之前先对各网络参数做一次量纲审计。8. 这套代码的边界与可以继续加装模块的方向8.1 目前适用的范围与明确不涉及的环节先说清楚边界免得同学拿去用的时候期望过高。这套代码目前面向的是三相平衡的稳态多能流计算因此它适用于输电网或者配电网的对称稳态分析。如果你要做三相不平衡配电网的电气热耦合计算需要把电网部分改成三相潮流模型工作量会明显上升。其次热网部分只考虑了质调节模式和稳态温度分布如果系统包含相变蓄热、季节性储热或者动态热惯性环节现有模型是不覆盖的。气网部分也没有处理压缩机站、储气库的动态充放气过程以及管网瞬变流动。最后一个限制是冷负荷只能通过电制冷机的 COP 模型做简单折算如果要严格做电-气-热-冷四网耦合的能流计算需要额外加一套冷网模型以及吸收式制冷机的工质循环方程。8.2 几个值得做的扩展方向代码的可扩展性其实在架构上已经预留了空间。结合我自己的经验以下几个方向是最常用到的一是增加热网水力工况的联立求解。目前的质调节模式把水力过程当作已知参数如果需要模拟变流量运行量调节则要把每个管道流量也作为状态变量加入迭代并补充环路压降方程。这个扩展会让状态变量数量显著增加但收敛框架不用变。二是耦合设备多运行模式的自动切换。真实的热电联产机组工作在以热定电以电定热最小凝汽等多个模式之间不同模式下热电比约束的形态不同。可以把模式切换逻辑做成一个外循环每次多能流收敛后检查设备运行点是否越出可行域越界则锁到边界值重新迭代一两次。这个思路在最优潮流、机组组合等上层优化问题里也非常有用。三是与上层优化算法嵌套。多能流计算是内层问题外层接遗传算法、粒子群或启发式搜索来优化设备容量配置、能源价格或碳排放约束这套代码的数据接口可以作为适应度函数的核心计算模块。因为整个程序是模块化的外层优化只需要修改耦合设备参数再调用主程序获得能流结果就能算出目标函数值。四是动态仿真方向的延伸。把稳态能流扩展为动态能流的最大难点在于热网的时间常数远大于电网和气网需要采用多时间常数技术multi-rate integration来协调不同网络的动态响应速度。目前稳态版代码中的雅可比矩阵结构可以直接沿用只需在时间迭代层加一个刚性问题求解器和事件检测逻辑。8.3 分享一点个人实操体会最后聊几句我自己的实操感受。做完这一整套多能流计算代码我最大的体会是相比于单独掌握电网潮流或者气网水力分析多能耦合的真正难点在于系统思维。每一个耦合设备的参数都会同时扰动两到三个网络的平衡方程调试时不能只看单个网络的残差而是要学会在整体状态空间里定位问题。建议刚开始做这个方向的同学先拿一个足够小的系统比如 3 节点电网加 3 节点气网再加 3 节点热网跑通统一求解法的完整流程然后把耦合设备一个一个加进去观察每一步收敛行为的改变。这种从简到繁、逐级验证的习惯会帮你省下大量在复杂算例里盲目试错的精力。这套代码还有挺多可以打磨的空间比如增加可视化交互界面、自动生成潮流分布图等但这些属于锦上添花的功能核心的物理模型和数值方法才是最花时间也最值得钻透的部分。希望这篇分享能给你带来一些有价值的参考。 
长河君 · 2026/9/30 18:07:54 · 阅读 1.2k 区域综合能源系统电气热能流计算：基于Matlab的统一求解法实现
1. 为什么区域综合能源系统的潮流计算绕不开多能耦合这几年做综合能源系统方向的朋友应该都有同感单独算电网潮流、天然气网水力分析或者热网热力工况的文章和代码已经非常成熟但一旦把三者放在同一个系统里联立求解问题立刻变得微妙起来。原因很简单——电、气、热三个网络并不是各自封闭运行的它们通过燃气轮机、电锅炉、热电联产机组、电制冷机、吸收式制冷机等耦合设备互相咬合。电网的负荷变化会影响燃气机组的出力进而改变天然气网的流量天然气网的压力和流量约束反过来又限制发电机的进气量最终影响电网的节点电压和功率分布。热网那边同样是这个逻辑热电联产机组的电出力与热出力之间存在运行区间约束热负荷波动会牵制电出力而电锅炉的投切又同时改变电网和热网两侧的平衡状态。正是这种强耦合特性让先分别算好再拼接的传统思路在综合能源系统里失去了可解释性。也正因为如此计及多能耦合的电气热能流计算才成为区域综合能源系统规划、运行、调度乃至可靠性分析的基础工具。这篇博文想分享的是一套基于Matlab实现的区域综合能源系统电气热能流计算代码重点讲清楚建模思路、耦合机制、解算方法以及我在调试过程中踩过的坑方便打算做这个方向或者正在做相关课题的同学直接参考甚至复用。按惯例先交代一下这套代码的适用范围它面向的是包含电力网络、天然气网络和热力网络的区域级综合能源系统耦合设备涵盖热电联产机组、燃气锅炉、电锅炉、电制冷机等常见类型。解算目标是求取稳态工况下的电网节点电压幅值与相角、气网节点压力与流量、热网节点供回水温度及管道流量。程序的整体思路是采用统一求解法又称全同时求解法即把电、气、热三个网络的方程连同耦合设备的约束一起组成一个大规模非线性方程组用牛顿-拉夫逊法迭代求解这也是目前多能流计算文献里精度和收敛性最稳妥的主流方案。2. 三个子网络单独建模时最容易忽略的物理边界2.1 电网部分其实就是带耦合项的标准交流潮流电网子模型沿用的是经典交流潮流模型。对于节点 i有功和无功功率平衡方程写为[ P_{i} U_{i}\sum_{j \in i}U_{j}(G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ][ Q_{i} U_{i}\sum_{j \in i}U_{j}(G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]这里的 (P_{i}) 和 (Q_{i}) 不是简单的负荷值而是要写成电源注入减去负荷消耗再减去耦合设备消耗的形式。例如节点上连着燃气轮机和电锅炉时[ P_{i} P_{G,i} P_{CHP,i} - P_{L,i} - P_{EB,i} ][ Q_{i} Q_{G,i} Q_{CHP,i} - Q_{L,i} - Q_{EB,i} ]其中下标 (G) 表示常规发电机组(CHP) 表示热电联产机组(EB) 表示电锅炉家 (L) 表示常规电负荷。这里有个细节很多初学者会漏掉电锅炉、电制冷机这类纯电耗设备通常以恒定有功功率形式建模但它们的无功功率往往被忽略。实际上电锅炉通过电力电子变换装置接入电网时会有一定的无功消耗如果系统规模不大、耦合设备占比高忽略无功可能会让迭代结果出现偏差。我在代码里把这类设备的无功按功率因数 0.95 估算虽然严格说要根据实际设备参数来定但至少比直接设零更合理。2.2 天然气网的管道方程不能只记一个公式天然气网络的稳态模型核心是管道流量方程。对于连接节点 (m) 和 (n) 的管道标准形式的稳态流量公式是[ f_{mn} C_{mn} \cdot s_{mn} \cdot \sqrt{s_{mn} \cdot (\pi_m^2 - \pi_n^2)} ]其中 (C_{mn}) 是管道常数由管径、长度、摩擦系数、气体温度、压缩因子等共同决定(\pi_{m}) 和 (\pi_{n}) 是节点压力符号系数 (s_{mn}) 的取值为[ s_{mn} \begin{cases} 1 \pi_m \geq \pi_n \ -1 \pi_m \pi_n \end{cases} ]这个式子的本质是压差平方与流量成正比而不是压差本身与流量成正比。实际编程中有几个地方要注意。第一如果直接用 (\pi_{m}^{2}) 和 (\pi_{n}^{2}) 作为状态变量方程形式会简单很多牛顿法的雅可比矩阵也更好求。第二管道常数 (C_{mn}) 的单位和数值范围很关键不同文献里有的用标准立方米每天有的用千克每秒混用会导致流量偏差几个数量级。第三天然气网中节点分为压力已知的平衡节点和压力待求的一般节点平衡节点通常设在气源或管网上游的门站这和电网中的平衡节点Vθ节点逻辑上是对应的。气网节点流量平衡方程按下式写[ \sum_{n \in m}f_{mn} G_{s,m} - G_{L,m} - G_{CHP,m} - G_{GB,m} 0 ]其中 (G_{s,m}) 是气源注入流量(G_{L,m}) 是常规气负荷(G_{CHP,m}) 和 (G_{GB,m}) 分别是热电联产机组和燃气锅炉消耗的天然气流量。这里要注意燃气机组消耗的天然气流量在气网侧是负荷但在电网侧是电源这种角色转换正是多能耦合在方程层面的体现。2.3 热网的质调节与量调节对应不同的方程形态热力网络的稳态模型比电网和气网都要复杂一些关键在于它涉及两个物理过程水力工况和热力工况。水力工况描述热水在管道中的流量分配由节点流量平衡和环路压降约束决定热力工况描述温度沿管道的分布以及节点供回水温度的混合。在实际工程中热网绝大多数时间运行在质调节模式即各管道流量保持不变通过调节供回水温度来适应热负荷变化。这种模式下可以先把水力工况算好甚至在程序中当成已知量输入只对热力工况迭代求解计算量会大幅下降。但如果是量调节模式流量本身也在变就必须联立求解水力与热力方程收敛难度会明显上升。我的代码默认采用质调节模式热力工况用节点法建模。对于热力管网节点温度分为供水温度和回水温度。热源节点处供水温度为设定值负荷节点处回水温度为设定值其他节点则满足流量加权温度混合方程[ \sum_{k \in S_{m}}(\dot{m}{k}c{p}T_{k}^{out}) T_{m}^{mix}c_{p}\sum_{k \in S_{m}}\dot{m}_{k} ]式子表示若干条汇入同一节点的管道热水混合后的温度等于节点混合温度其中 (c_{p}) 是水的比热容(\dot{m}_{k}) 是各支路的质量流量。管道出口温度则按热损失方程衰减[ T_{j}^{out} (T_{j}^{in} - T_{a}) \cdot \exp\left(-\frac{\lambda L}{\dot{m} c_{p}}\right) T_{a} ]这里 (\lambda) 是管道单位长度传热系数(L) 是管长(T_{a}) 是环境温度。如果管道保温效果好或者系统规模不大也可以暂时忽略热损失但实际工程算出来的温度衰减对后续设备效率计算是有影响的所以代码里保留了这个项方便后续扩展到动态热网模型。热力节点方程最终要满足两类约束每个节点的功率平衡供热量等于负荷侧吸收的热量和各环路的压降为零如果不做水力计算则可以跳过。质调节模式下给定总循环流量和节点热负荷每个节点的供回水温度就是待求状态变量负荷节点的回水温度由热负荷和流量共同决定[ T_{r,i} T_{s,i} - \frac{\Phi_{L,i}}{c_{p}\dot{m}_{i}} ]其中 (\Phi_{L,i}) 是该节点的热负荷功率。这个式子简单但它是热网与电网耦合的另一个入口——如果该节点的热量来自热电联产机组热源侧还要再建一个热出力与电出力对应的耦合方程。3. 耦合设备模型整个计算中最容易出错的环节3.1 热电联产机组的热电比约束热电联产CHP机组是综合能源系统里最核心的耦合设备它的运行特性通常用热电比约束来描述。最常见的是定热电比模型[ \Phi_{CHP} c_{m} \cdot P_{CHP} ]其中 (c_{m}) 是热电比。定热电比模型虽然简单但在实际系统里过于理想化——真实的热电联产机组的工作区间是一个多边形电出力和热出力之间存在可行域边界。更精确的建模方式是用线性可行域来约束典型的一组约束写为[ P_{CHP,min} \leq P_{CHP} \leq P_{CHP,max} ][ 0 \leq \Phi_{CHP} \leq \Phi_{CHP,max} ][ P_{CHP} \geq a_{1}\Phi_{CHP} b_{1} ][ P_{CHP} \leq a_{2}\Phi_{CHP} b_{2} ]这四条约束合在一起刻画了背压式和抽汽式机组的不同运行区间。在潮流计算里处理这种不等式约束的方式通常是先不管它迭代收敛后检查是否越界如果越界则把耦合设备出力定格在边界值上再重新迭代一两次。后面我会专门讲这个处理逻辑。3.2 燃气轮机和燃气锅炉的气耗特性燃气轮机消耗的天然气和电出力之间的关系可以表示为[ G_{GT} \frac{P_{GT}}{\eta_{GT} \cdot LHV_{gas}} ](\eta_{GT}) 是燃气轮机的发电效率(LHV_{gas}) 是天然气低位热值通常取 9.7 kWh/m³ 左右。实际计算中要注意单位统一问题——如果电网侧功率单位是 MW气网侧流量单位是 m³/h那么低位热值必须统一成 MWh/m³ 才能保证等式两边的量纲一致。我最初写代码时在这里栽过跟头算出来的气耗比文献值大了一倍查了半天最后发现是热值单位换算出了错。这个看似不起眼的问题其实是很多同学在做多能流计算时收敛错误甚至结果荒诞的根源。燃气锅炉的气耗模型类似[ G_{GB} \frac{\Phi_{GB}}{\eta_{GB} \cdot LHV_{gas}} ]唯一不同的是输入是热出力 (\Phi_{GB}) 而不是电出力。燃气锅炉可以作为分布式热源配置在热网的任意节点上在热负荷高峰时投入。3.3 电锅炉与电制冷机电转热/冷的单向耦合电锅炉和电制冷机的拓扑逻辑相对简单它们只体现电 → 热/冷的单向转换不存在热或冷回馈到电网的通路。电锅炉的模型就是一个效率系数[ \Phi_{EB} \eta_{EB} \cdot P_{EB} ]电制冷机的模型是制冷系数COP[ Q_{cool} COP_{EC} \cdot P_{EC} ]这两个设备本身没有额外状态变量但在电网方程里要计入它们的电功率消耗在热网/冷网方程里要作为热源/冷源来处理。从编程角度说它们的耦合强度不如 CHP 和燃气轮机但数量多、分布广在整体雅可比矩阵里依然会形成非零交叉块不能省略。4. 统一求解法如何把电、气、热三个方程塞进同一个牛顿迭代框架4.1 状态变量统一排队与初值选取策略统一求解法的核心思想很简单把电网、气网、热网的全部未知量放在同一个状态向量里所有网络方程和耦合方程构成同一个残差向量然后整体用牛顿法迭代。以我实现的代码为例状态变量与残差方程的数量和排列决定了雅可比矩阵的结构也直接决定了程序的通用性和扩展难度。我采用的状态向量排列顺序为电网部分PQ 节点的电压幅值 (U) 和相角 (\theta)PV 节点的相角 (\theta)不包括平衡节点其电压幅值和相角已知气网部分所有非平衡气节点的压力平方 (\pi^{2})热网部分所有非给定温度节点的供水和回水温度质调节模式下对应的残差方程排列顺序与状态变量一一对应这样雅可比矩阵的构造最简单调试时也便于定位问题。整个系统的雅可比矩阵从结构上看是一个分块矩阵[ J \begin{bmatrix} J_{EE} J_{EG} J_{EH} \ J_{GE} J_{GG} J_{GH} \ J_{HE} J_{HG} J_{HH} \end{bmatrix} ]其中 (J_{EE})、(J_{GG})、(J_{HH}) 是三个网络各自的自雅可比块非对角块 (J_{EG})、(J_{GH}) 等是网络间的耦合雅可比块。如果忽略多能耦合这些非对角块全为零问题退化为三个独立网络分别求潮流只有耦合设备存在时非对角块才真正不为零这也是多能耦合在数学上的本质。初值选取方面电网部分采用平启动flat start——PQ节点电压幅值取 1.0 p.u.相角取 0气网节点压力初值取气源压力的 80% 左右热网供回水温度初值取设计工况的典型值。平启动对牛顿法而言是最常见的策略但对于强耦合的大规模系统一旦初值离解太远迭代很容易发散。我的经验是给雅可比矩阵加一个对角阻尼因子——在迭代初期把阻尼因子设得大一些比如 0.5随着残差下降逐步减小到 1.0这套阻尼牛顿法在综合能源系统多能流计算里表现相当稳定。4.2 雅可比矩阵分块计算的程序实现思路要在 Matlab 里高效实现分块雅可比矩阵最高效的做法是分别对电、气、热三个子网络模块写各自的解析雅可比函数再由主程序把分块矩阵拼装起来。每个子模块的输入都是整个状态向量和该子网络自身的节点数据输出是该子网络对应残差对该子网络状态变量的偏导块。耦合块的偏导藏在耦合设备方程里需要针对每一类耦合设备单独写。以 CHP 机组为例假设它连接电网的节点 i 和气网的节点 m同时向热网节点 k 供热。它引入的偏导关系包括[ \frac{\partial P_{i}}{\partial \pi_{m}^2}, \quad \frac{\partial G_{m}}{\partial P_{i}}, \quad \frac{\partial \Phi_{k}}{\partial P_{i}} ]这三条链式关系会让整个雅可比矩阵出现三个非零的交叉块。具体数值计算时如果用解析法需要根据耦合方程组逐项求偏导这对代码结构的模块化要求比较高如果嫌麻烦也可以用数值差分近似这些交叉项即对每个耦合变量施加一个微小扰动观察各网络残差的变化量来估计偏导值。数值差分的精度和步长选择比较敏感但胜在实现简单、不易出错。我的代码里提供了解析法和数值差分法两种模式默认打开解析法方便对照验证。在 Matlab 的实现层面程序主体结构可以这样组织% 主函数入口IES_PowerFlow.m % 输入电网数据、气网数据、热网数据、耦合设备数据 % 输出各网络状态变量、收敛信息 state initState(elec, gas, heat); % 状态变量初始化 residual computeResidual(state, elec, gas, heat, coupling); J assembleJacobian(state, elec, gas, heat, coupling); iter 0; while norm(residual, inf) tol iter maxIter delta -J \ residual; % 阻尼因子 lambda 1.0; while norm(computeResidual(state lambda * delta, ...), inf) norm(residual, inf) lambda 0.5 * lambda; if lambda 1e-4, break; end end state state lambda * delta; residual computeResidual(state, elec, gas, heat, coupling); J assembleJacobian(state, elec, gas, heat, coupling); iter iter 1; end这个流程和常规电力潮流计算的牛顿法是同构的多出来的部分只有残差函数和雅可比函数需要同时涵盖气网和热网的方程。程序的骨架清晰之后真正麻烦的其实在于耦合方程接入残差向量时必须保证状态变量与方程的数目严格相等——多一个未知量而少一个方程矩阵就奇异迭代必然发散。这是统一求解法最考验工程细节的地方。4.3 迭代收敛判据与常见发散原因排查收敛判据我采用的是无穷范数[ \max\left( \max_i |\Delta U_i|, \max_i |\Delta\theta_i|, \max_i |\Delta\pi_i^2|, \max_i |\Delta T_{s,i}|, \max_i |\Delta T_{r,i}| \right) \varepsilon ]这个判据比单独看残差要严格因为残差小并不一定代表状态变量修正量小——在病态系统里两者可能会出现数量级不一致的情况。建议同时监控两个指标修正量无穷范数和残差无穷范数以修正量为主判据残差作为辅助参考。收敛容差设置在 (10^{-6}) 左右比较合理太严格会拖慢速度太宽松结果不够准。在调试过程中发散的情况大多可以归结为以下几类原因第一变量-方程数量不匹配。这是新手最容易犯的错。电网部分除了平衡节点外每个节点都有一个有功方程PQ 节点还有一个无功方程气网每个非平衡节点有一个流量平衡方程热网每个节点有一个温度混合方程。算清楚这些之后再和状态变量核对确保一一对应。第二初值偏离实际解太远。气网压力初值尤其敏感因为管道流量方程里有 (\pi^{2}) 项如果初始压力设得太低迭代过程中平方项产生巨大的梯度变化会导致震荡发散。建议气源压力较高的节点附近初值取最高压力的 80% 以上。第三耦合设备的效率参数或热电比设置不合理。比如 CHP 效率设为 0.9再叠加燃气轮机效率 0.4气耗算出来的结果就会让气网负荷暴涨导致气网压力跌破下限整个迭代过程不收敛。这不是算法问题而是模型参数问题排查时要先看物理上合不合理。5. 一个典型算例从仿真结果看耦合效应的实际影响5.1 算例配置与边界条件用一个 6 节点电网、6 节点气网、6 节点热网的小型区域综合能源系统来演示程序效果。电网包含 1 个平衡节点、2 个 PV 节点常规发电机组和 3 个 PQ 节点。气网包含 1 个气源节点和 5 个负荷节点。热网包含 1 个热源节点由 CHP 供热的换热站和 5 个热负荷节点。耦合设备配置为1 台 CHP 机组连接电节点 2、气节点 3、热节点 1、1 台燃气锅炉连接气节点 5 和热节点 4、1 台电锅炉连接电节点 5 和热节点 5。各网络的关键参数如下表所示参数数值电网基准容量100 MVA气源节点压力4.0 MPa对应平方 16 MPa²热网供水温度110 ℃热网回水温度设计值60 ℃CHP 热电比1.2CHP 发电效率0.42燃气锅炉效率0.9电锅炉效率0.98负荷水平方面电网总负荷设定为 150 MW气网总气负荷设为 80 m³/h换算成热值约 776 MW 等价值热网总热负荷为 60 MW。CHP 承担其中 40 MW 的供热燃气锅炉承担 15 MW电锅炉承担 5 MW。5.2 迭代收敛过程程序从平启动初值开始迭代逐次残差变化如下迭代次数修正量无穷范数残差无穷范数01.0初始36.7210.42319.8720.08752.1430.00690.17240.00020.00653.8×10⁻⁶1.2×10⁻⁴62.3×10⁻⁸7.8×10⁻⁷第 6 次迭代达到收敛迭代次数在牛顿法的正常范围内。这个收敛模式说明系统参数设置合理统一求解法的效率与纯电力潮流相当没有出现因耦合引入导致的明显迭代次数增加。5.3 耦合效应量化对比单网络解耦 vs 统一求解为了直观展示多能耦合对计算结果的实际影响我把代码里耦合设备的交互全部断开相当于把所有耦合方程替换为恒定功率注入分别跑三个独立的常规潮流再将结果与统一求解法做对比。差异最明显的以下几个量状态量解耦计算统一求解偏差CHP 所在电节点电压幅值 (p.u.)1.0451.0212.3%电网平衡节点有功出力 (MW)71.568.23.3 MW气源节点输出流量 (m³/h)105.496.78.7 m³/h热网负荷节点回水温度 (℃)62.365.83.5 ℃这个对比说明解耦计算在很多情况下虽然不至于得到完全荒谬的结果但偏差已经足以影响工程决策。尤其是当 CHP 热出力受到热负荷变化牵制时电网侧的有功平衡会受到显著影响——如果按解耦结果做调度电网平衡节点的出力安排会偏高气网气源采购量也会产生非忽略的偏差。对运行优化或者可靠性分析来说这种偏差积累到最后可能让结论反向。从计算效率上看解耦计算确实比统一求解快那么一点少迭代一两次但考虑到底层物理的完整性这个成本是完全值得的。这也是为什么当前多能流领域的主流文献几乎都推荐统一求解法。6. 代码架构设计模块化拆分比一个脚本全塞进去靠谱得多6.1 文件结构安排与接口约定很多同学拿到一个题目就习惯把所有逻辑写在一个大脚本里跑通这样在几十行的小程序里没有问题但一旦系统规模拉大或者需要换算例代码的可维护性会快速恶化。这套代码做成了模块化结构每个网络一个接口文件、一个数据定义文件和一个残差/雅可比函数文件核心文件结构如下IES_PowerFlow/ ├── main_IES_PF.m % 主程序数据加载、初始化、迭代 ├── data/ % 数据目录 │ ├── case6_elec.m % 6节点电网数据 │ ├── case6_gas.m % 6节点气网数据 │ └── case6_heat.m % 6节点热网数据 ├── elec/ │ ├── elec_residual.m % 电网残差 │ └── elec_jacobian.m % 电网雅可比 ├── gas/ │ ├── gas_residual.m % 气网残差 │ └── gas_jacobian.m % 气网雅可比 ├── heat/ │ ├── heat_residual.m % 热网残差 │ └── heat_jacobian.m % 热网雅可比 └── coupling/ ├── coupling_model.m % 耦合设备参数读取与状态映射 ├── coupling_residual.m % 耦合方程残差 └── coupling_jacobian.m % 耦合方程雅可比接口约定上每个子网络函数都采用统一的输入输出格式。输入参数是完整的状态向量 state 以及该子网络自身的结构体数据 data输出是残差向量及其雅可比。这样主程序不需要知道每个函数内部是怎么实现的只要按约定调用即可。想换成别的算例时只需要在 data 目录下新增一组数据文件或者把现有算例的节点和管道参数改掉。6.2 统一状态向量与数据索引防止变量错位的关键统一求解法的难点在于把三个网络的状态变量放在同一个向量里时如何保证每个残差方程对应的偏导存储在雅可比矩阵的正确位置。我的做法是维护一张状态变量-索引映射表在主程序中预先定义idx.U 1:n_elec_PQ; % PQ节点电压幅值索引 idx.Theta idx.U(end)1:idx.U(end)n_elec_PVn_elec_PQ; % 相角索引 idx.Pi idx.Theta(end)1:idx.Theta(end)n_gas_slackless; % 气网压力平方索引 idx.Ts idx.Pi(end)1:idx.Pi(end)n_heat_unknown; % 供水温度索引 idx.Tr idx.Ts(end)1:idx.Ts(end)n_heat_unknown; % 回水温度索引在残差函数里每个子网络函数只需要按照映射表取自己需要的状态变量子集计算完残差后输出到与映射表对应的位置即可。雅可比矩阵的组装则通过一个稀疏矩阵来累积每个子网络函数返回的是它自己的那部分行和列的稀疏矩阵主程序累加到一起。稀疏矩阵的重要性怎么强调都不过分。如果按全稠密矩阵存储6节点系统还能忍受但扩展到几十几百个节点时矩阵规模迅速膨胀存储和求解效率都会崩。Matlab 里用 sparse 命令把雅可比矩阵声明为稀疏格式配合反斜杠求解计算速度可以提高一两个数量级。6.3 网络拓扑编号与物理连接关系的映射计算中另一个需要谨慎处理的是电网、气网、热网三套编号体系之间的对应关系。电网的节点编号是 1 到 6气网和热网各自也有自己的一套编号它们之间唯一的联系是耦合设备表——每台耦合设备记录了自己挂在电网哪个节点、气网哪个节点、热网哪个节点。程序运行时耦合设备模型要根据这张设备表把不同网络的变量提取出来再组装成耦合方程。这里我采用的是一种面向对象的思路虽然 Matlab 的结构体数组也能实现但用 classdef 定义耦合设备类会让代码可读性更好。每个耦合设备实例包含自己的类型、效率参数、接入节点信息、状态变量索引。在组装雅可比时主程序遍历所有耦合设备实例每个实例把自己对应的耦合残差和耦合雅可比块填入全局矩阵的相应位置。7. 那段改了一个参数整个系统发散的排查经历7.1 问题现象明明是同一套程序算例一换就崩代码写完之后我先拿 6 节点系统调通了一切正常。后来想验证通用性把电网换成一个 33 节点的配电网络气网和热网也在对应规模上做了扩充结果一跑就发散——残差在前几步一直下降到了第 4 次迭代突然往上跳之后直接飞掉。我首先怀疑的是初值问题把气网初值调到最高压力的 90%热网温度初值也做了调整但依然发散。随后我逐个检查状态变量和残差的对应关系确认数量没错又用数值差分法对比了解析雅可比两者的结果在单个子网络上是一致的。7.2 排查链路从变量数量检查到耦合块数值病态真正发现问题是在检查雅可比矩阵条件数的时候。我用 condest 估算了整体雅可比矩阵的条件数发现已经超过 (10^{15})几乎处于奇异状态。进一步拆解每个分块发现电网部分的自雅可比块条件数正常气网部分的条件数在 (10^{7}) 左右也还可以接受。问题出在耦合块和热网自块上——热网部分的自雅可比块条件数高达 (10^{12})因为 33 节点系统里热网管道长度差异过大最短的管道和最长的管道差了三个数量级导致温度衰减指数项差异极大雅可比矩阵中的元素出现了严重的量级失衡。7.3 解决方案量纲归一化与支路参数预处理解决办法分两步。第一步是对热网管道参数做预处理把传热系数和管长合并成无量纲的温度衰减因子[ f_{loss} \exp\left(-\frac{\lambda L}{\dot{m} c_{p}}\right) ]这个因子本身就是 0 到 1 之间的数直接用它替代原方程里的指数表达式避免指数项在变量微扰下产生异常敏感的变化。第二步是对整个热网做量纲归一化处理把温度全部表示为相对于环境温度的归一化值使所有方程里涉及的变量都在 0.1 到 10 这个量级范围内。修改之后同一套程序跑 33 节点系统顺利收敛迭代次数 8 次残差降到 (10^{-7}) 以下。这个经历给我最大的启发是多能流计算的收敛性问题很多时候不是算法本身的缺陷而是数值问题的表现。物理模型没问题方程推导没问题但数值尺度不匹配时再好的算法也会失效。做这类程序一定要养成检查雅可比矩阵条件数的习惯以及在做大规模算例之前先对各网络参数做一次量纲审计。8. 这套代码的边界与可以继续加装模块的方向8.1 目前适用的范围与明确不涉及的环节先说清楚边界免得同学拿去用的时候期望过高。这套代码目前面向的是三相平衡的稳态多能流计算因此它适用于输电网或者配电网的对称稳态分析。如果你要做三相不平衡配电网的电气热耦合计算需要把电网部分改成三相潮流模型工作量会明显上升。其次热网部分只考虑了质调节模式和稳态温度分布如果系统包含相变蓄热、季节性储热或者动态热惯性环节现有模型是不覆盖的。气网部分也没有处理压缩机站、储气库的动态充放气过程以及管网瞬变流动。最后一个限制是冷负荷只能通过电制冷机的 COP 模型做简单折算如果要严格做电-气-热-冷四网耦合的能流计算需要额外加一套冷网模型以及吸收式制冷机的工质循环方程。8.2 几个值得做的扩展方向代码的可扩展性其实在架构上已经预留了空间。结合我自己的经验以下几个方向是最常用到的一是增加热网水力工况的联立求解。目前的质调节模式把水力过程当作已知参数如果需要模拟变流量运行量调节则要把每个管道流量也作为状态变量加入迭代并补充环路压降方程。这个扩展会让状态变量数量显著增加但收敛框架不用变。二是耦合设备多运行模式的自动切换。真实的热电联产机组工作在以热定电以电定热最小凝汽等多个模式之间不同模式下热电比约束的形态不同。可以把模式切换逻辑做成一个外循环每次多能流收敛后检查设备运行点是否越出可行域越界则锁到边界值重新迭代一两次。这个思路在最优潮流、机组组合等上层优化问题里也非常有用。三是与上层优化算法嵌套。多能流计算是内层问题外层接遗传算法、粒子群或启发式搜索来优化设备容量配置、能源价格或碳排放约束这套代码的数据接口可以作为适应度函数的核心计算模块。因为整个程序是模块化的外层优化只需要修改耦合设备参数再调用主程序获得能流结果就能算出目标函数值。四是动态仿真方向的延伸。把稳态能流扩展为动态能流的最大难点在于热网的时间常数远大于电网和气网需要采用多时间常数技术multi-rate integration来协调不同网络的动态响应速度。目前稳态版代码中的雅可比矩阵结构可以直接沿用只需在时间迭代层加一个刚性问题求解器和事件检测逻辑。8.3 分享一点个人实操体会最后聊几句我自己的实操感受。做完这一整套多能流计算代码我最大的体会是相比于单独掌握电网潮流或者气网水力分析多能耦合的真正难点在于系统思维。每一个耦合设备的参数都会同时扰动两到三个网络的平衡方程调试时不能只看单个网络的残差而是要学会在整体状态空间里定位问题。建议刚开始做这个方向的同学先拿一个足够小的系统比如 3 节点电网加 3 节点气网再加 3 节点热网跑通统一求解法的完整流程然后把耦合设备一个一个加进去观察每一步收敛行为的改变。这种从简到繁、逐级验证的习惯会帮你省下大量在复杂算例里盲目试错的精力。这套代码还有挺多可以打磨的空间比如增加可视化交互界面、自动生成潮流分布图等但这些属于锦上添花的功能核心的物理模型和数值方法才是最花时间也最值得钻透的部分。希望这篇分享能给你带来一些有价值的参考。