做了大半年区域综合能源系统电气热能流计算从最初只会单独跑电网潮流到把气网、热网全部耦合起来在Matlab里写成一套能稳定收敛的代码中间踩过的坑比想象中多得多。这个项目的完整标题是《计及多能耦合的区域综合能源系统电气热能流计算研究Matlab代码实现》核心就一句话当电、气、热三种异质能源通过热电联产CHP、燃气锅炉、电锅炉这些设备耦合在一起时怎么用一套可复现的数值算法算出整个系统稳态运行的各种关键状态量。很多初学者以为这项任务就是把三个网的潮流程序简单堆在一起实际上完全不是这么回事。电网算的是节点复功率平衡气网算的是节点压力和管段流量平衡热网更麻烦水力计算和热力计算得交替进行。三种能网的时间常数、量纲、建模思路都不一样耦合变量又来回反馈稍不注意就会发散。这篇博文把我从模型搭建、方程推导到Matlab落地实现、调参避坑的完整过程梳理一遍内容包括电气热网稳态模型、耦合设备建模、顺序迭代求解框架、雅可比矩阵组装、典型算例分析以及实际调试中遇到的一堆典型问题。适合正在做综合能源系统、园区多能互补、微电网方向研究的研究生和工程师参考。1. 思路拆解多能耦合到底“耦”在哪里1.1 三种能网的计算范式差异先看基础。电网的稳态潮流分析核心是节点功率平衡方程已知负荷和机组注入功率求解各节点电压幅值和相角。方程是光滑的二次型牛顿-拉夫逊法收敛性很好只要初值不离谱四五次迭代基本都能收住。这是三种网络里最成熟、参考资料最多的一个。气网的稳态计算就不太一样了。它要满足的是节点气体流量平衡即流入某个节点的天然气流量减去流出流量等于该节点的用气负荷。描述管段流量的公式常用Weymouth方程这个方程带压力平方差以及绝对值符号在管道两端压力相等的位置导数不存在容易引发迭代震荡。我实际调试时发现气网方程整体是单调的但收敛半径比电网小不少初值差一点就可能在压力低的位置出现负的压力平方直接报错。热网是最麻烦的。它的稳态运行状态需要分两步描述水力工况决定各管段质量流量热力工况决定各节点供水和回水温度。水力方程有两大类一类是节点流量平衡方程流入流出必须相等另一类是环路压头平衡方程也就是沿闭合回路流动的压头损失代数和为零。流量在方程里以平方项出现所以是非线性的。热力方程又建立在流量基础之上温度计算时流量待在分母上初值给不好一样发散。三张网络合在一起本质上是一个多物理场耦合的非线性方程组问题。我自己的体会是如果只把三套程序按顺序各跑一遍那不叫多能耦合计算那只是三个独立的工具串行使用。真正的耦合计算要求三个子系统之间通过耦合设备逐轮交换边界条件一直迭代到整个系统的运行状态达到一致。1.2 耦合设备是系统的“关节”区域综合能源系统里电网、气网、热网并不是孤立存在的它们通过能量转换设备连在一起。常用的耦合设备有热电联产机组CHP、燃气锅炉、电锅炉、热泵甚至电转气P2G。这类设备在能流计算里扮演的角色很清晰上游能网给它们提供输入能量下游能网承接它们的输出能量转换关系则由各自效率决定。我用一个表格把最常见的几种设备整理了出来方便后续建模时直接对照耦合设备输入能网输出能网关键转换参数热电联产CHP天然气电力 热力电效率ηe、热效率ηh、电热比r燃气锅炉天然气热力热效率ηb电锅炉电力热力电热转换效率ηe2h热泵电力热力制热性能系数COP电转气P2G电力天然气转换效率η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 / ∂θjNij Vj * ∂Pi / ∂VjJij ∂Qi / ∂θjLij 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 LA是节点-支路关联矩阵q是各管段质量流量L是节点注入或流出的质量流量。第二步列出环路压头平衡方程B * (K * q * |q|) 0B是环路-支路关联矩阵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 牛顿法迭代与收敛判据顺序迭代的外循环流程分为五步初始化所有网络的状态变量包括电网电压幅值和相角、气网节点压力平方、热网管段流量和节点温度。求解电网潮流计算CHP等电源节点的新电出力更新电网节点注入功率。求解热网水力方程得到各管段流量再解热力方程得到各节点供回水温度更新CHP热出力。把CHP和燃气锅炉消耗的天然气量叠加到气网负荷节点上求解气网节点压力。检查前后两次外循环各耦合变量的变化量如果都小于容差则收敛否则回到第二步继续迭代。收敛判据我设置为两个条件必须同时满足一是每个网络内部的功率或流量不平衡量小于设定容差电网取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 end4. 算例设计与结果分析4.1 测试系统与计算条件我用的测试系统是自己拼装的区域综合能源系统算例规模不算大但足以说明问题。电网部分参考了IEEE 33节点配电系统的拓扑骨架加了两台分布式燃气轮机和若干负荷天然气网络取7个节点、8条管道带一个气源节点和两条环状结构热网由16个节点组成包含3个独立环路负荷节点覆盖供热区域中的主要热力站。基准值统一取为电功率100MW、气流量100kg/s、热功率100MW压力用MPa作单位温度用摄氏度。CHP机组共3台燃气锅炉1台电锅炉1台。算例场景设定为典型冬季供暖日热负荷32MW电负荷18MWCHP承担约七成热负荷其余由燃气锅炉和电锅炉分担这是北方园区供能系统很常见的一种运行方式。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 常见故障速查表我把调试过程中遇到的典型问题整理成了一张速查表新接触这套代码的人遇到状况可以先对照排查故障现象可能原因处理办法热网雅可比矩阵奇异管段流量初值为零阻力项导数失效按负荷反推非零流量初值外循环震荡不收敛耦合变量一步更新过大加入0.6~0.8阻尼因子气网压力平方出现负值初值压力过低或气源压力偏小提高气源参考压力检查负荷量级电网潮流迭代发散PV节点无功越限未处理将越限PV节点转成PQ节点重新迭代CHP出力持续越限CHP可行域内的功率分配不合理调整机组间出力分配比例收敛判据同时满足但结果不合理各网络收敛容差过松收紧热网温度容差到1e-4℃计算耗时突然增大矩阵以满阵形式存储改用sparse稀疏存储6. 一些个人经验体会这个项目做下来我最大的感触是多能耦合能量流计算真正的技术含量不在算法推导本身而在于工程调试的耐性和系统性。我花在调初值、查矩阵组装、核对量纲上的时间远远超过写方程的时间。如果你也要做类似的工作我强烈建议按这个顺序推进先把电网潮流跑通并和标准算例比对再把气网单独算通把热网的水力工况单独调稳定最后才做三者耦合。每个环节都先用一个非常小的两节点系统验证结果确认和手算值一致以后再进行下一步。这套流程听起来慢实际上是最快的一条路它能帮你在源头上堵住大部分bug。另外还有一个值得分享的小技巧调试耦合迭代时可以把每一轮外循环的耦合变量打印出来比如CHP的天然气耗量、热出力、电出力观察它们的演化轨迹。如果发现变量在长期震荡而不是单调收敛说明阻尼因子太小或者某个子系统内部的容差设置太松如果变量直接从初值跳到一个不合理的边界值说明耦合设备限幅逻辑可能有问题。这比单纯盯着发散信息排查效率高得多。这套电气热耦合能流计算框架目前的扩展空间还很大。比如把气网模型升级成非等温模型把热网加入管道热惯量动态特性又或者在耦合设备层加入储能装置和需求侧响应逻辑都能在这个框架上做下去。代码本身也写了完整的子函数模块后续做参数优化或者多场景扫描时只需要在外层套一个循环就行。如果你正在做园区多能互补、分布式能源系统规划或者综合能源调度相关的课题这套方法和代码可以作为你的第一版基础工具。