这几年区域综合能源系统Integrated Energy System, IES的课题越来越多园区级的冷热电三联供、区域级的电-气-热联合规划几乎成了电力方向毕业设计和横向项目里的热门方向。而不管具体场景怎么变第一道绕不过去的坎基本都是同一个多能流计算。很多人上手直接用电力潮流的思路去套天然气网和热力网结果要么算出来数值很诡异要么迭代半天根本不收敛最后卡在仿真这一步后面优化调度、可靠性分析全都没法往下走。这篇文章把我自己用 Matlab 实现区域综合能源系统电气热能流计算的完整过程梳理出来从数学模型搭建、耦合设备建模到交替迭代求解器的代码结构再到实际跑算例时踩过的坑一次性讲清楚。适合正在做 IES 方向毕业设计、需要调论文仿真的硕士生也适合刚接触园区综合能源规划的工程技术人员参考。1. 项目核心思路先把网系和耦合分清楚1.1 单一能流计算的舒适区与盲区传统电力潮流、天然气稳态水力分析、供热管网水力计算单独拿出来都有成熟的算法和商业软件。电力潮流用牛拉法或快速解耦法天然气网用 Weymouth 稳态方程热力管网用图论基本回路法这些都是教科书级别的成熟内容。但综合能源系统的难点不在任何一个单独的子系统而在多个能源网络之间的边界条件不再独立。热电联产机组把天然气变成电和热燃气锅炉把天然气变成热电转气P2G把电变成天然气电锅炉和热泵把电变成热。任意一个负荷波动都会跨网传播。如果分开做单能流计算就会在边界条件上犯一个很隐蔽的错误电网调度以为 CHP 能发出 3 MW 电但气网实际供给的天然气只够发 2.5 MW整个系统根本达不到计算假设的运行点。所以计及多能耦合这几个字不是修饰而是问题的本质。多能流计算的最终目标是给整个 IES 找一个一致的稳态运行点——在这个运行点上电网潮流方程满足气网节点压力方程满足热网温度与流量方程满足而且所有耦合设备的输入输出能量关系也满足。1.2 典型耦合设备及其能量关系既然要写 Matlab 代码第一步就是把耦合设备抽象成数学关系。我在项目里用的是最常用的一批设备模型CHP热电联产输入天然气输出电和热。背压式机组可以简化成固定电热比抽凝式机组则需要用可行域描述。简化模型写成F_gas · LHV P_e / η_e Q_h / η_h Q_h c_m · P_e (背压式固定电热比)燃气锅炉GB输入天然气输出热。Q_h η_gb · F_gas · LHV电锅炉EB输入电输出热。Q_h η_eb · P_e热泵HP输入电输出热从低温热源取热。Q_h COP · P_eP2G电转气输入电输出天然气。F_gas η_p2g · P_e / LHV这些公式很简单但它们是连接三个网络的钩子。在交替迭代法中耦合设备的输入侧是一个网络的边界比如 CHP 耗气量是气网的负荷输出侧又是另一个网络的电源或热源CHP 电出力是电网的注入功率热出力是热网的热源功率。搞清楚每个设备的能量流向和效率矩阵比背潮流公式重要得多。1.3 统一求解 vs 交替求解我为什么先交替多能流计算的求解路线分两大派统一求解法把电、气、热所有方程联立构造一个大型稀疏雅可比矩阵一次牛顿迭代同时更新所有状态量。优点是二阶收敛解的一致性最好缺点是矩阵规模大、方程结构复杂尤其热网要考虑供水/回水双网络导致雅可比矩阵的分块结构和常规电力潮流差异很大代码调试门槛高。交替求解法分解协调法每轮迭代分别求解电网、气网、热网三个子系统的能流方程通过耦合设备变量交换信息。优点是模块化强电网潮流、气网方程、热网水力/热力模型可以分别写成独立的函数排错方便缺点是迭代次数通常比统一法多强耦合场景下可能振荡。我在实际做这个课题时强烈建议先用交替法把整体框架跑通验证各网络之间的耦合关系是否正确再视情况往统一法方向扩展。交替法的代码结构和论文里的系统框图天然对应审稿人读起来也舒服。下面整个项目就是按交替法展开的。2. 电气热能流的核心数学模型2.1 电网潮流牛拉法极坐标方程电网部分不需要重新开发算法直接用经典牛顿-拉夫逊法。极坐标下每个节点的注入功率不平衡量为ΔPi Pi_spec - Vi · Σj Vj [ Gij·cos(θij) Bij·sin(θij) ] ΔQi Qi_spec - Vi · Σj Vj [ Gij·sin(θij) - Bij·cos(θij) ]节点分成 PQ、PV、平衡三类雅可比矩阵由四个分块构成∂P/∂θ、∂P/∂V、∂Q/∂θ、∂Q/∂V。迭代格式是J · Δx -ΔF状态量更新为电压幅值和相角。这里有个细节在多能流框架里电网节点类型需要根据耦合设备重新设定。比如 CHP 并网节点如果按给定电出力处理可以视为 PQ 节点或 PV 节点取决于是否带无功调节能力电锅炉和 P2G 节点则是在原有负荷上增加一个电负荷项这个电负荷大小又取决于热网/气网的返回结果。电网求解器的输入不能是固定负荷向量必须是一个由耦合变量动态拼接的边界条件。2.2 天然气网Weymouth 方程与压缩机天然气稳态流动最常用的是 Weymouth 方程Fij sign(πi^2 - πj^2) · Cij · sqrt( |πi^2 - πj^2| )其中 π 为节点压力Fij 为管段流量Cij 为管道常数由管径、长度、摩擦系数、气体特性决定。节点流量平衡方程是Σj∈i Fij S_i L_iS_i 为气源注入量L_i 为节点负荷取流入节点为正。气网节点分类与电网类似气源节点压力给定类似平衡节点负荷节点流量给定。如果管网上有压缩机还需要加入压缩机模型。最简单的处理是给定升压比 Rc使得出口压力 Rc × 进口压力同时压缩机消耗一部分天然气可用固定比例或线性函数近似。压缩机模型让气网雅可比矩阵多了一组约束我在初期实现时先不加入压缩机跑通后再扩展这样能有效降低调试难度。2.3 热网潮流水力模型与热力模型分开写热网计算比电气两网都麻烦因为要分成水力和热力两步而且两者互相耦合。水力模型求管道流量 q节点流量连续方程A · q GG 是节点净流出流量回路压降方程Bf · h 0管道压降方程h_l R_l · q_l · |q_l|求解方法用基本回路法给定一个初始流量分布计算每个基本回路的压降代数和用牛顿法修正回路流量直到所有回路压降为 0。热力模型求温度 T节点功率方程Φ_i Cp · q_i · (T_s,i - T_o,i)供水管道温度沿程损失T_end Ta (T_start - Ta) · exp(-λ·L / (Cp·q))多股热水汇合节点温度按流量加权平均T_mix Σ(q_in · T_in) / Σ(q_in)热网的关键难点在于水力方程中的节点流量 G Φ / (Cp·ΔT)而 ΔT 是热力计算的结果热力方程计算温度分布时又需要已知管道流量。因此热网内部本身就是一个小型迭代我把它设计成一个嵌套循环固定负荷热功率 → 假设一个 ΔT 初值 → 求水流流量 → 求温度分布 → 反推新的 ΔT → 再修正流量。这个细节很关键很多人在这一步翻车直接用热负荷除以固定温差算流量然后温度分布怎么算都不收敛就是因为忽略了水力-热力的双向耦合。2.4 耦合设备的统一抽象为了让主程序更干净我把所有耦合设备抽象成一个统一的coupling_model函数。它的本质是一个多输入多输出的映射关系设备输入侧网络输出侧网络关键参数CHP气网耗气电网电、热网热η_e、η_h / 电热比燃气锅炉气网耗气热网热η_gb电锅炉电网耗电热网热η_eb热泵电网耗电热网热COPP2G电网耗电气网产气η_p2g初始化时给每个设备设一组初值每一轮外部迭代后由这个函数根据效率公式更新设备的输入/输出量再传给对应子网络求解器。耦合设备的变量更新方向绝对不能搞反比如 P2G 在电网上是负荷消耗电在气网上是气源注入天然气方向反了迭代一定振荡发散。3. Matlab 代码实现框架与关键模块3.1 数据结构设计用 struct 组织三个网络Matlab 做这类计算最自然的数据组织方式是结构体struct。我建议在项目根目录下建下面这种结构IES_MultiFlow/ ├── main_IES_Flow.m % 主程序 ├── data/ │ └── case_IES_5g4h6.m % 算例数据电网5节点气网4节点热网6节点 └── funcs/ ├── pf_grid.m % 电网牛拉法潮流 ├── gas_nodal.m % 气网节点压力求解 ├── heat_hydraulic.m % 热网水力求解 ├── heat_thermal.m % 热网热力求解 ├── coupling_model.m % 耦合设备映射 └── calc_jacobian_grid.m % 电网雅可比矩阵数据文件里用结构体定义各个网络的拓扑和参数例如% 电网5节点含1平衡节点、1个CHP并网节点、1个电锅炉节点 grid.nb 5; grid.Ybus [...]; % 节点导纳矩阵 grid.type [3; 1; 1; 1; 1]; % 1PQ, 2PV, 3平衡 grid.Pload [0; 0.2; 0.4; 0.3; 0.5]; % MW grid.Qload [0; 0.1; 0.15; 0.12; 0.2]; % MVar % 气网4节点1个气源 gas.nn 4; gas.C [...]; % 管道Weymouth常数 gas.pi0 [7; 6.5; 6.2; 6.0]; % 节点压力初值 (MPa) gas.L [0; 0; 0; 0]; % 负荷在迭代中由耦合设备更新 % 热网6节点含2个热源CHP、燃气锅炉3个热负荷 heat.nh 6; heat.A [...]; % 节点-支路关联矩阵 heat.Bf [...]; % 基本回路矩阵 heat.R [...]; % 管道阻力系数 heat.Cp 4.182; % 水的比热容 kJ/(kg·K) heat.Ta 25; % 环境温度 ℃提示数据文件与求解逻辑分离后续换算例只需要新增一个 case 文件主程序一行不用改。这个习惯在做多能流系列仿真时能节省大量时间。3.2 电网牛拉法求解器的核心实现电网求解器我有现成的模板核心部分长这样function [V, theta, iter] pf_grid(Ybus, Sbus, type, V0, theta0) % 极坐标牛顿-拉夫逊潮流 % Sbus: 复数节点注入功率发电机-负荷PQ/PV/平衡节点均计入 % type: 节点类型向量1PQ 2PV 3平衡 n length(Sbus); V V0; theta theta0; tol 1e-8; maxIter 50; for iter 1:maxIter % 计算注入功率 Vc V .* exp(1j * theta); Ibus Ybus * Vc; Scalc Vc .* conj(Ibus); % 不平衡量PQ节点取P和QPV节点只取P平衡节点不取 dP real(Sbus) - real(Scalc); dQ imag(Sbus) - imag(Scalc); % 组装dF根据type向量选择) ... % 雅可比矩阵省略具体分块可参考标准教科书 J calc_jacobian_grid(Ybus, V, theta, type); % 解方程并更新 dx J \ dF; theta theta dx(1:n); V V dx(n1:end); if max(abs(dx)) tol, break; end end end实际调试时要注意Sbus 是多能流框架给出的当前轮耦合边界下的净注入功率在每次外部迭代里都要重新拼装。我习惯把 CHP 并网节点的注入功率、电锅炉节点的负荷、P2G 节点的负荷都写成一个build_grid_boundary(coup)函数返回当前轮次电网需要的 Sbus这样主循环会清爽很多。3.3 气网求解牛顿法处理 Weymouth 方程气网的实现也是牛顿法只是状态量是节点压力方程是节点流量不平衡量。对每个流量给定节点偏差量为ΔF_i Σj∈i Fij S_i - L_i雅可比矩阵里需要 Weymouth 方程的导数。我直接用解析形式Fij sign(Δπ) · Cij · sqrt(|Δπ|), 其中 Δπ πi^2 - πj^2 ∂Fij/∂πi Cij · πi / sqrt(|Δπ|) ∂Fij/∂πj -Cij · πj / sqrt(|Δπ|)注意当 Δπ 接近 0 时导数趋于无穷实际计算时给分母加一个很小的常数eps兜底。气源节点压力给定不参与不平衡量计算并在迭代后固定压力值。核心函数骨架function [pi, iter] gas_nodal(gas, L) % gas: 气网结构体, L: 节点天然气负荷正值为耗气 p gas.pi0; for iter 1:50 F weymouth_flow(p, gas.C, gas.branch); dev sum(F, 2) gas.S - L; % 节点流量不平衡量 J weymouth_jacobian(p, gas.C, gas.branch, gas.S_type); dp J \ dev; p p dp; if max(abs(dp)) 1e-6, break; end end pi p; end3.4 热网水力与热力求解内部嵌套迭代热网水力求解我用基本回路法先用节点流量方程得到一个可行的初始支路流量分布再以回路压降为零作为修正方程用牛顿法迭代修正回路流量。热力求解则在已知流量的基础上自热源节点向负荷节点逐管道推演供水温度再从负荷侧回推回水温度。由于水力与热力耦合我在模块内部套了一层循环伪代码如下function [q_h, Ts, To] heat_solve(heat, Phi_load, Ts_supply_hot) % 初始假设供回水温差为 30℃ dT 30 * ones(heat.nh, 1); for inner 1:20 % 由热负荷和温差求节点流量需求 G Phi_load ./ (heat.Cp * dT / 1000); % 注意单位换算Cp用kJ/(kg·K) % 水力求解基本回路法 q_h heat_hydraulic(heat, G); % 热力求解供水管温度推演 回水管混合温度 [Ts, To] heat_thermal(heat, q_h, Ts_supply_hot); % 更新节点温差 dT_new Ts - To; if norm(dT_new - dT, inf) 1e-4, break; end dT 0.5 * dT 0.5 * dT_new; % 阻尼更新防止振荡 end end注意热量单位的换算非常容易出错。热功率 Φ 如果单位是 MW而 Cp 是 kJ/(kg·K)流量 q 是 kg/s则需要满足Φ (MW) (Cp×q×ΔT) / 1000。很多人算出来流量大几个数量级基本都是这个换算错误。3.5 交替迭代主循环与收敛判据这些都是子模块真正把三个网络拧在一起的是主循环。主程序框架如下% 初始化耦合设备变量 coup.P_chp 3.0; % CHP电出力 MW coup.Q_chp 3.5; % CHP热出力 MWth coup.Q_gb 1.0; % 燃气锅炉热出力 MWth coup.P_eb 0.8; % 电锅炉电功率 MW coup.F_gas_total 0; % 气网总供气 alpha 0.7; % 松弛因子 maxOuter 30; for k 1:maxOuter % 1. 由当前耦合变量构建电网边界求电网潮流 Sbus build_grid_boundary(grid, coup); [V, theta] pf_grid(grid.Ybus, Sbus, grid.type, V0, theta0); % 从电网结果更新以电为输入的设备电锅炉、P2G、热泵 % 电锅炉/热泵耗电是给定的分布式控制量时也可不更新 coup.P_eb_new ...; coup.P_p2g_new ...; % 2. 求解热网此时热源为CHP热出力燃气锅炉热出力 Phi_source [coup.Q_chp; coup.Q_gb; 0; 0; 0; 0]; Phi_load [0; 0; 1.2; 0.8; 1.0; 0]; % 热负荷 MWth [q_h, Ts, To] heat_solve(heat, Phi_load, T_supply_set); % 3. 由热平衡检查CHP/燃气锅炉当前热出力是否满足负荷 % 若热源是可调度的则实际热出力由热网潮流结果决定 coup.Q_chp_new ...; coup.Q_gb_new ...; % 4. 求解气网负荷为CHP耗气燃气锅炉耗气-P2G产气 L_gas [coup.F_chp_gas; coup.F_gb_gas; 0; 0] - coup.F_p2g; [pi_gas] gas_nodal(gas, L_gas); % 5. 由气网结果更新CHP的实际供电出力气源约束下 coup.P_chp_new coup.F_chp_gas * LHV * coup.eta_e / 3600; % 6. 松弛更新耦合变量检查收敛 coup.P_chp (1-alpha)*coup.P_chp alpha*coup.P_chp_new; coup.Q_chp (1-alpha)*coup.Q_chp alpha*coup.Q_chp_new; err abs(coup.P_chp_new - coup.P_chp) ... abs(coup.Q_chp_new - coup.Q_chp); if err 1e-5, break; end end外层迭代的收敛判据我同时检查两个量一是耦合变量更新量的绝对值二是各子网络求解器返回的最大不平衡量残差。只查耦合变量有时会漏掉电压压力温度虽然变了但设备功率没变的假收敛。4. 算例测试与结果分析4.1 测试系统设计为了验证代码我搭了一个小规模测试算例拓扑关系如下电网5 节点。节点 1 为平衡节点节点 2 为 CHP 并网节点PQ给定注入功率节点 3 为电锅炉节点节点 4、5 为纯电负荷节点。气网4 节点。节点 1 为气源压力给定 7 MPa节点 2 为中间输气节点节点 3 为 CHP 耗气节点节点 4 为燃气锅炉耗气节点。热网6 节点。节点 1 为 CHP 热源节点 2 为燃气锅炉热源节点 3、4、5 为热负荷节点节点 6 为管网中间汇合节点管道按树干式连接含一个基本回路。设计这个系统的好处是三个网络的耦合节点全部拉开便于观察每个设备在迭代中的行为。如果是自己想快速验证也可以从电网气网电转气的两网耦合开始跑通之后再加热网排查范围会小很多。4.2 典型计算结果在一个轮次的收敛解上典型的运行状态如下电网电压幅值在 0.975 ~ 1.02 pu 之间CHP 节点电压略低于平衡节点电锅炉节点因为有集中电负荷电压相对偏低。气网气源节点 7 MPaCHP 节点压力约 6.2 MPa燃气锅炉节点约 6.0 MPa基本符合输气压力沿程递减的趋势。热网供水温度设定 70 ℃远端负荷节点供水温度降到 64 ℃左右回水温度约 45 ℃供回水温差约 20 ℃。整个热网的热功率平衡误差控制在 1e-4 MW 以内。需要说明具体数值与网络参数强相关我这里给的是合理的数量级不是通用标准。真正能说明问题的是收敛行为交替法在松弛因子 α0.7 时大约 8~12 次外迭代内达到 1e-5 的收敛精度。同一算例改用统一求解牛顿法4~6 次迭代即可收敛。交替法的迭代次数对松弛因子非常敏感。α0.3 时迭代 20 次以上才收敛α1.0无松弛时CHP 强耦合场景下经常出现等幅振荡怎么都不收敛。4.3 收敛性能与初值敏感性我专门做了一组敏感性子测试结论很有代表性气网压力初值影响最大。压力初值偏离真实值太远比如把 7 MPa 的管道初始化成 3 MPaWeymouth 方程的 sqrt 项在迭代初期会出现很大的导数牛顿法容易走飞。热网流量初值不能随便拍脑袋。我用按热负荷占总负荷比例分配总流量的方式初始化支路流量比均匀初始化收敛快接近一倍。压缩机模型如果加入收敛速度会明显变慢因为压缩机节点引入了不等式约束和额外的耗气量方程。建议先不加压缩机组算出一个基准结果再加进来对比。注意交替法的本质是固定点迭代耦合强度越高、设备电热比越固定越容易发散。这时候不要一味降低松弛因子整体迭代次数会增加很多更好的做法是改用分层迭代内层气-电迭代外层再耦合热网或者直接上统一求解法。5. 常见问题与避坑指南5.1 收敛性差先查方向再查初值最后查松弛我在帮别人调试时见过最多的三类不收敛原因耦合变量方向反了。P2G 节点在电网是负荷在气网是气源如果两边符号设置一致必炸。初值离解太远。特别是热网供水温度初值如果设成和环境温度一样温降指数项算出来几乎为零温度推演直接失真。没有阻尼更新。直接用新值替换旧值在 CHP 电热比固定的系统里几乎必然振荡。排查顺序建议先固定其他两个网络单独做电-气两网迭代确认收敛后加入电-热耦合最后把气-热通过 CHP 合并。分步验证能帮你快速定位是哪一组耦合关系的更新逻辑出错。5.2 单位与基值的统一多能流系统最坑的地方是单位混用。我整理了一张自用速查表物理量常用单位与功率的折算关系电功率MW直接参与潮流计算热功率MWth与电功率同量纲注意区分天然气流量Nm³/h1 Nm³ 天然气热值约 35.6 MJ即 9.9 kWh1 Nm³/h ≈ 0.0099 MW水的流量kg/s热功率 4.182 × 流量 × 温差 / 1000 (MW)压力MPa 或 Pa注意 Weymouth 方程中压力单位统一为 Pa 或 MPa 保持一致有一个非常隐蔽的问题天然气热值有高位热值HHV和低位热值LHV之分CHP 效率书里常写 η_e 0.35基于 LHV。如果代码里误用 HHV气网算出的耗气量明显偏小电网侧 CHP 出力却不变最终气网压力会高得离谱。5.3 耦合设备模型的简化陷阱论文里用固定效率模型很快工程应用就要小心。固定效率模型的问题在于CHP 的电出力和热出力实际上是强耦合的抽凝式机组的运行域是一个梯形或多边形可行域固定电热比只是背压式机组的特例。我在后期扩展时把 CHP 模型改成了线性可行域模型P_e_min ≤ P_e ≤ P_e_max Q_h_min ≤ Q_h ≤ Q_h_max Q_h c1 · P_e ≤ c2 (线性约束近似)这类约束不直接进牛顿法而是在多能流收敛后做一次可行域投影修正效果比硬加进雅可比矩阵稳定得多。5.4 Matlab 实现细节几个提高效率、缩短调试周期的习惯雅可比矩阵用 sparse 存储。Matlab 对稀疏矩阵的\运算有专门优化稠密矩阵在节点数超过 50 时速度差异非常明显。避免在循环里用 eval 或者动态变量名。多能流代码需要反复跑上百次参数敏感性分析脚本性能很重要。网络数据统一用数组和稀疏矩阵不要用eval([node, num2str(i)])这种写法。牛顿迭代不必每轮都更新雅可比矩阵。气网和热网的雅可比矩阵在迭代中变化不大可以每隔 3~5 次更新一次用简化牛顿法能省不少时间。初期调试用完整雅可比后期性能优化再切换。用 R2020b 以上版本。本文这套代码在 R2020b、R2023b、R2025 上我都跑过逻辑完全兼容。新版本对稀疏矩阵求解和隐式扩展有优化但老版本跑小算例也够用。把中间结果打印成表格。我在主循环里每隔一轮输出当前 CHP 电出力、热出力、气负荷、电压最大偏差用fprintf格式化输出排查问题时比看断点直观得多。6. 一些扩展思路多能流计算的代码跑通之后相当于给后续研究打了一个底座。我自己的经验是后续这几个方向都能直接复用这套框架优化调度多能流作为校核层把调度方案映射到运行点验证可行性。概率多能流将负荷、风光出力设为随机变量用蒙特卡洛或多点估计法跑几百次多能流求电压、压力、温度的均值与置信区间。代码只需要在最外层加一个抽样循环。碳流分析在多能流收敛的基础上沿能量流追踪碳排放的网间转移做园区碳足迹核算。热网动态扩展稳态热力模型换成准动态模型管道热惯性可以进一步做热负荷柔性分析。如果研究周期允许我建议把统一求解法和交替求解法都实现一版对比它们的收敛性和计算时间这也是一篇论文里很扎实的对比内容。交替法实现简单、排错容易统一法收敛快、理论完整两者互补。我个人的体会是做多能流计算最大的障碍不是理论看不懂而是三个成熟算法在接口处的细节处理。只要把每个网络当前轮次的边界条件在哪、由哪个耦合设备决定这个逻辑理清楚代码写起来其实非常快。最后再分享一个小技巧第一版先别追求通用性针对一个固定算例硬编码拓扑和参数跑通结果后再花半天把数据结构抽象成通用框架。反向操作——一边设计通用框架一边调 bug——大概率会陷入无穷无尽的调试循环。希望这套 MatLab 实现思路能帮你在综合能源系统研究里少走几段弯路。