多能耦合区域综合能源系统电气热能流计算与Matlab实现

前阵子给一个园区做多能互补方案，我先用常规潮流程序把电力网络算了，结果发现燃气锅炉和CHP的出力根本定不下来——热负荷一变，电网机组的出力就得跟着变，单独算电网等于在打移动靶。后来狠下心把电网、气网、热网放到同一个Matlab框架里联立求解，才算是把整个区域综合能源系统的稳态运行点一次算清楚。

这篇内容就是围绕“计及多能耦合的区域综合能源系统电气热能流计算”这个题目展开的，我会把数学模型、耦合元件建模、Matlab代码实现、算例验证和调试经验全串起来讲。适合正在做多能互补、园区综合能源规划、微能网运行优化的学生和工程师参考，尤其是手头有项目但不知道如何把三个网络方程统一到一套代码里的朋友。

1. 为什么电气热联立求解比想象中复杂得多

1.1 多能耦合的本质：能量转换设备把三个网络拴在一起

所谓多能耦合，本质就是电网、气网和热网之间不再彼此独立，而是通过CHP机组、燃气锅炉、电锅炉、热泵这些能量转换设备产生了强关联。

举个例子，CHP机组一端消耗天然气，另一端同时输出电能和热能。在热负荷较高的季节，CHP被迫多发一些热，电出力也会跟着提高；如果电网侧此时并不缺电，多出来的电就会从平衡节点倒送，或者改变其他发电机的出力。这些连锁反应只有把电、气、热三个网络放在同一个方程体系里解，才能得到自洽的结果。

换句话说，耦合设备就是“桥梁”，它们把电网节点方程、气网节点方程、热网节点方程中的若干项串了起来。你在代码里只要漏掉任何一侧的耦合项，计算结果就会在能量平衡上出现明显缺口，而且你很难凭直觉判断缺口来自哪里。

1.2 三个网络的数学差异：不是改个符号就能复用

很多人第一反应是，电网潮流都能算了，照葫芦画瓢把电阻改成管道阻力，不就能算热网气网了吗？实测下来完全不是这么回事。三个网络的状态变量和方程性质差异非常大，我先用表格梳理一下：

电网的复功率方程是电压、相角的二次非线性；气网的Weymouth方程直接跟压差的平方根打交道；热网更特殊，它分为水力工况和热力工况两层，水力方程里的“压降与流量平方成正比”带有明显的绝对值非线性。你说它们难，倒不是公式本身有多复杂，而是这些非线性特性在牛顿法的雅可比矩阵里表现得很不一样，初值给不好、符号写错、单位不一致，都会导致迭代发散。

1.3 我的方案：统一求解法，而不是顺序迭代

处理多能潮流，业内主要有两条技术路线。

• 顺序迭代法：先假定耦合设备的某一侧出力，比如 CHP 电出力固定，算出电网潮流；再拿着电出力去算热网和气网，得到新的热/气侧状态后返回修正 CHP 参数，如此反复迭代。
• 统一求解法：把所有网络的方程和耦合约束合并成一个大型非线性方程组，用扩展的牛顿-拉夫逊法一次性求解所有变量。

顺序迭代法写起来直观，但碰到强耦合场景容易在几次迭代之间来回震荡，尤其当热负荷和电负荷同时波动时，CHP 的出力在两个解附近反复跳，收敛非常痛苦。我在实际项目里最终选了统一求解法，虽然装配矩阵的代码量更大，但鲁棒性明显好一截，后续如果要加储能、加氢能设备，也只需要在变量向量和雅可比矩阵里“加行加列”，扩展性可控。

2. 电气热网络的数学模型：先把方程组列对

2.1 电网方程：极坐标下的经典潮流

电网部分沿用常规潮流计算的极坐标牛顿-拉夫逊法。对节点 i，有功和无功平衡方程写作：

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。

节点类型仍然是三种：平衡节点提供相位基准，PV节点给定有功和电压幅值，PQ节点给定有功无功负荷。在多能系统里，CHP 机组的电气侧通常建模为 PQ 或 PV 节点，这取决于控制方式；但它的有功出力不是预先给定的常数，而是同时被热负荷和气网侧的天然气供给约束，所以在统一求解时这个“给定值”会在联立方程的过程中被反复修正。

2.2 气网方程：Weymouth管道模型和压缩机

气网稳态计算里有两条核心约束。

第一条是节点流量平衡：对任意气网节点，进气源的注入流量加上从相邻管道流入的流量，等于该节点的气负荷加上泄漏损失（一般忽略泄漏）。写成矩阵形式就是 A * f + L_source - L_load = 0，A 是节点-管道关联矩阵，f 是管道流量向量。

第二条是管道流量方程，工程上最常用 Weymouth 公式：

f_pq = K_pq * sign(π_p - π_q) * sqrt(abs(π_p^2 - π_q^2))

这里的 π 是节点气压的平方，K_pq 是管道参数，跟管径、长度、摩阻系数、气体性质有关。

压缩机在气网里很关键，它负责升压，典型建模方式有定出口压力、定升压比、定耗气量三种。我用的 Matlaab 代码里默认支持定升压比模式，即 p_out = CR * p_in，同时压缩机自身要消耗一部分天然气，可以作为该节点的附加负荷进入流量平衡方程。如果你的实际系统输气距离不长、压力变化不大，可以把压缩机先当成理想元件跳过，但后期做高压输气网络时一定要补回来。

2.3 热网方程：水力工况和热力工况分开建模

热网比电网麻烦在它有两层方程。

水力工况层：每个节点的流量守恒方程，加上每个回路的压力平衡方程。管道压降方程为：

h = K_h * m * abs(m)

其中 m 是管道质量流量，K_h 是管道的阻力特性系数。正因为有 abs(m) 这一项，方程在流量变号附近不可导，牛顿法收敛时会特别敏感。

热力工况层：忽略动态传热，稳态热功率方程为：

φ = C_p * m * (T_supply - T_return)

管道沿程温降采用指数衰减模型：

T_end = T_a + (T_start - T_a) * exp(-λ * L / (C_p * m))

多个支路汇入同一个节点时，混合温度按流量加权平均，这一点在代码里特别容易写错。我刚开始就是直接取了算术平均值，结果热负荷大的一侧温度算高，整体能量平衡始终对不上。

3. 耦合元件建模：CHP、燃气锅炉、电锅炉怎么进入方程组

3.1 CHP机组的热电联产约束

CHP 是电气热三网耦合最核心的设备。最简单的稳态模型是电热比线性约束：

H_chp = c_m * P_chp

c_m 是热电比，抽凝式机组通常在 0.8~1.5 之间，背压式可能更高。CHP 消耗的天然气流量为：

G_chp = P_chp / (η_e * LHV)

或者更精细一些，把电效率和热效率分开折算燃料消耗：

G_chp = P_chp / (η_e * LHV) + H_chp / (η_h * LHV)

在统一求解的残差方程里，CHP 会同时出现在三个位置：电网节点注入项、气网节点负荷项、热网节点热源项。换句话说，你每给 CHP 分配一个变量编号，雅可比矩阵里就有三块对应区域需要填入非零偏导数，这是“多能耦合”在数学上的直观体现。

3.2 燃气锅炉与电锅炉：相对简单的两网耦合

燃气锅炉消耗天然气、输出热水，关联气网和热网：

G_gb = H_gb / (η_gb * LHV)

电锅炉消耗电能、输出热能，关联电网和热网：

P_eb = H_eb / η_eb

热泵也是电转热，但用的是逆卡诺循环，耗电量与供热量之间用 COP 联系：

P_hp = H_hp / COP

为了把这几种设备在代码里统一管理，我维护了一张耦合设备映射表：

这张表的价值在于，写代码前先想清楚每个设备在三个网络的“支点”在哪里，后面组装配残差和雅可比时就不会漏项。

耦合元件的变量，我建议在一开始就作为全局变量放进总变量向量 x 里，而不是在每次迭代后单独计算。这样做的好处是让牛顿法的搜索方向同时包含所有网络的可行修正量，避免子网络之间出现“互相猜”的现象。

4. Matlab代码实现：从零搭建计算主循环

4.1 总变量编号与数据结构

统一求解的第一步是定义总变量向量。我把所有待求状态量按顺序拼接成一个列向量：

x = [ 电网相角θ; 电网PQ节点电压幅值V; 气网节点气压π; 热网节点供水温度T_s; 热网节点回水温度T_r; 热网管道流量m; 耦合设备内部变量 ]

每个变量都要分配一个全局索引。我习惯用结构体保存三个网络的原始数据：

data.elec.bus = struct('id',{}, 'type',{}, 'V0',{}, 'th0',{}, 'Pd',{}, 'Qd',{}); data.elec.branch = struct('fbus',{}, 'tbus',{}, 'r',{}, 'x',{}, 'b',{}); data.gas.node = struct('id',{}, 'type',{}, 'p0',{}, 'load',{}, 'source',{}); data.gas.pipe = struct('fnode',{}, 'tnode',{}, 'K',{}, 'comp',{}); data.heat.node = struct('id',{}, 'type',{}, 'Ts0',{}, 'Tr0',{}, 'phi',{}, 'flow',{}); data.heat.pipe = struct('fnode',{}, 'tnode',{}, 'Kh',{}, 'L',{}, 'Ta',{});

不建议把变量零散地放在一堆全局变量里，那样后期改网络规模时很容易索引错乱。用结构体统一挂载数据，再用一个索引映射函数把“网络类型+节点编号”转换成主变量 x 中的位置，是维护这套代码最划算的方式。

4.2 残差函数：每行残差对应一个物理方程

主计算循环围绕“残差函数”和“雅可比矩阵”两个核心展开。残差函数 F(x) 按顺序返回所有方程的不平衡量。

电网部分残差是节点有功无功不平衡量；气网部分残差是节点流量不平衡量；热网部分残差包括水力残差（节点流量平衡、回路压降）和热力残差（热功率方程、管道温降方程）。

管道流量的核心函数可以这样写：

function f = weymouth_flow(K, pi_k, pi_m) d = pi_k^2 - pi_m^2; if abs(d) < 1e-12 f = 0; else f = K * sign(d) * sqrt(abs(d)); end end

关键是把 sign() 和 sqrt(abs()) 用对，否则当管道两端气压相等时，流量方向会在迭代中反复乱跳，导致整体发散。

热力节点混合温度的计算也放一个小函数：

function T_mix = mix_temperature(m_in, T_in) % m_in: 流入节点的各支路流量向量 % T_in: 对应支路末端温度向量 T_mix = sum(m_in .* T_in) / sum(m_in); end

这个函数看起来简单，但它对应的物理含义是：供暖管道多路热水汇合时，混合温度等于流量加权平均，而不是算术平均。

4.3 雅可比矩阵组装：先数值差分验证，再上解析式

统一求解法的装配工作量主要集中在雅可比矩阵 J 上。J 中的元素是每个残差方程对每个变量的偏导数。三个网络内部有各自的非零块，耦合设备又在三块区域之间产生耦合非零块。

组装方式我推荐分三步走。

第一步用数值差分验证逻辑：

eps0 = 1e-6; J(:, j) = (F(x + eps0*ej) - F(x - eps0*ej)) / (2*eps0);

先对一个小算例做数值差分雅可比，如果牛顿迭代都能收敛，说明残差方程本身逻辑是对的。

第二步再写解析雅可比替换，提升计算速度。不要一上来就写解析式，那样很难定位到底是残差写错了还是求导求错了。

第三步用稀疏矩阵存储：

J = sparse(I_idx, J_idx, val_list, n_var, n_var);

三个网络拼起来之后变量数量很容易超过几百个，稠密矩阵求逆在 matlab 里虽然也能跑，但每迭代一次的时间会随变量数快速增长，尤其是做参数扫描或运行优化时，速度差距非常明显。

4.4 牛顿迭代主循环

主循环可以压得很精简：

x = x0; for k = 1:maxIter F = compute_residual(x, data); J = assemble_jacobian(x, data); dx = -(J \ F); x = x + dx; if norm(F, inf) < tol && norm(dx, inf) < tol break; end end

这里我用了双重收敛判据：一边看残差范数，一边看变量增量范数。只盯变量增量会出现一种假象——步长小不代表方程满足，尤其是在热网这种温度变量数量级很小的场景里，ΔT=1e-6 看起来已经很“收敛”了，但对应的热功率不平衡可能还很大。

5. 算例验证：一个小型区域系统跑通全流程

5.1 算例拓扑与参数

为了验证代码，我搭了一个小型区域综合能源系统：

• 电网：3个节点，节点1为平衡节点，节点2接入CHP机组，节点3接电锅炉和常规电负荷，额定容量100 kVA，基准电压10 kV；
• 气网：4个节点，节点1为气源，节点2给CHP供气，节点3给燃气锅炉供气，节点4带常规气负荷；
• 热网：4个节点，节点1连接CHP作为热源1，节点2连接燃气锅炉作为热源2，节点3和4分别带热负荷100 kW和120 kW，热力管网包含3条管道。

这部分参数我在建模时故意选得很接近园区实际情况：热网供水温度初值80℃，回水初值50℃，电锅炉额定热出力60 kW，热泵COP取3.2，CHP热电比取0.93。

5.2 迭代过程与收敛性观察

以“平启动”初值出发，统一求解法在本算例中的迭代过程大致如下：

这里我解释一下，该数据是在特定初值、特定参数下得到的结果，不具备普遍性，但能直观说明一个现象：统一法前期收敛很快，后期进入线性收敛区，整体11步以内就压到了浮点精度附近。

收敛后的关键运行状态如下：

• 电网侧：节点1电压1.00∠0°，节点2电压0.982∠-2.1°，节点3电压0.976∠-3.5°；
• 气网侧：气源节点1气压1000 kPa，节点2约996 kPa，节点3约993 kPa；
• 热网侧：节点1供水温度80.0℃，节点3供水温度76.2℃，节点4供水温度74.5℃；回水温度分别约49.3℃、48.6℃、47.9℃；
• 耦合设备运行点：CHP电出力300 kW，热出力279 kW；燃气锅炉热出力150 kW；电锅炉热出力60 kW，电功率约67 kW；热泵热出力45 kW，电功率约14 kW。

5.3 结果合理性分析

收敛之后不能直接宣告胜利，还要做校核。我习惯检查三类能量平衡：

• 电网有功平衡：平衡节点实际出力 + CHP电出力 ≈ 电负荷 + 电锅炉耗电 + 热泵耗电 + 网损；
• 气网流量平衡：气源总注入 ≈ CHP燃气耗量 + 燃气锅炉燃气耗量 + 常规气负荷 + 管损；
• 热网热平衡：热源总供热量 ≈ 热负荷 + 管道散热损失。

我的Matlab代码在校核后，三项平衡的最大相对误差都在1e-6量级以内，说明联立解是自洽的。这里补一句：如果你在跑完自己的算例后能量不平衡误差很大，九成以上是某个耦合设备漏在了某个网络的残差里，而不是牛顿法迭代次数不够。

6. 踩坑记录与收敛性调优经验

6.1 热网初值：供水温度必须大于回水温度

这句话看着像废话，但在我调试初期真的翻过车。我给热网全部节点初值设置为一样的80℃，回水也为80℃，结果供水回水温差为0，热功率方程残差对温度的导数出现异常，雅可比矩阵条件数急剧恶化。

后来我把初值策略改成：供水各节点统一80℃，回水各节点统一50℃，管道流量初值取一个与热负荷匹配的估计值，再让牛顿法自己去修正。这一步改动非常小，收敛稳定性却提升了一个档次。

6.2 气网符号函数与低压差附近的处理

Weymouth方程的绝对值符号是头号天坑。在迭代早期，管道两端气压很可能接近，此时如果 sign(d) 写得不对，或者 sqrt(abs(d)) 的导数没有平滑，牛顿步长会猛冲一下，然后第二、三迭代步直接发散。

我的处理是：在计算雅可比矩阵中 Weymouth 对气压的偏导时，给 d=π_k^2-π_m^2 加一个很小的死区或平滑，让导数不会在原点处变成无穷大。工程上可以理解为，微小压差下流量对压差的变化率不可能真的无穷大，做一点约束既符合物理常识，又能救回迭代稳定。

6.3 三网单位不统一导致的病态矩阵

电网习惯用标幺值（电压1.0左右），气网用kPa（几百到上千），热网温度用℃（几十到一百），管道流量用 kg/s（可能只有个位数）。把这些直接塞进同一个变量向量，雅可比矩阵非零元素的量级会从1e-2到1e3跨度巨大，矩阵条件数非常差，迭代后期精度会很尴尬。

我的做法是对变量做基准化：温度除以100，气压除以1000，流量除以各自的基准值，让所有变量基本落在0.01~10区间内。这一步做完，收敛速度明显变快，而且不会出现“前半段正常，后半段始终差一点达不到容差”的怪现象。

6.4 顺序迭代法的教训：强耦合时不要死磕

我在做第一版代码时其实是先写的顺序迭代法，因为觉得简单。结果是：只要热负荷和 CHP 热电比稍微敏感一点，外层循环就在“CHP多发电→热侧偏热→调低电出力→电网不平衡→提高CHP出力”的循环里反复横跳，设定最大迭代次数也救不回来。

后来我统计了一下震荡点，根本原因在于 CHP 同时被电、气、热三个约束牵制，顺序法的每一层都在“临时固定”另外两个网络的变量，信息传递存在一步延迟。统一求解法之所以稳健，是因为它在一个非线性方程组里同时处理了所有物理约束，三类残差共享同一个变量向量，各变量在同一个迭代步内同步更新，不需要人为设定收敛次序。

6.5 建议的调试路径：按网络逐层叠加

如果一开始就搭完整的三网联立系统，出问题之后几乎没法定位。有一个非常实用的调试策略：先在纯电网+纯热网的小系统上把耦合设备只有电锅炉/热泵的情况跑通；再加入气网，此时只需要固定燃气锅炉和CHP的气侧消耗，把气网作为“跟随网络”接入；最后才把CHP的电-热-气三联互动打开，做完整统一求解。

每加一层网络之前，保留一个可运行的存档版本，方便对照调试。这套思路帮我节省了大量时间，也建议你复现类似项目时采用。

一点个人的实操体会

最后单独说一个我反复测试中认为最有用的技巧：不要迷信“全解析雅可比”。只要条件允许，先用数值差分雅可比把模型逻辑验证清楚，再根据矩阵稀疏模式用 symbolic 工具辅助求导，最后人工校正关键耦合块的偏导数。这样虽然前期麻烦一点，但能大幅降低“模型对但代码错”的隐蔽 bug 概率。

多能耦合潮流计算在数学上并没有比常规电力潮流高深太多，难的是三套方程、三类变量、多种单位体系要在同一个架构里稳定共处。我的体验是，只要建模一致性守住了，收敛性基本不会成为瓶颈；如果你的迭代总是发散，先回头检查耦合设备的能量流方向是不是统一了。把功夫下在建模和数据结构上，Matlab 这套代码完全可以在几十个节点的园区级规模上稳定运行。希望这篇文章能帮你少走一些我走过的弯路。