搞区域综合能源系统的人，几乎都绕不开电气热能流计算这块硬骨头。所谓电气热能流，不是三条独立潮流的简单拼盘，而是考虑多能耦合后，电、气、热三种网络在一个系统里联立求解的运行状态。我之前在做一个园区级综合能源系统方案时，因为要用Matlab实现计及多能耦合的电气热能流计算，前前后后折腾了小两个月。这篇记录我会把建模思路、代码实现细节、调参踩坑经验都摊开来讲，适合正在做综合能源方向研究的学生、刚入行的系统工程师，以及打算自己搭稳态分析工具的朋友参考。

1. 项目背景与核心思路拆解

1.1 为什么区域综合能源系统需要“多能耦合”视角

传统电力系统潮流计算只管电网，天然气管道水力计算只管气网，供热管网水力热力计算只管热网。但在区域综合能源系统里，电网、气网、热网不是孤立的，它们通过热电联产机组（CHP）、燃气锅炉、热泵、电锅炉、P2G设备联系在一起。任何一个网络内部的变化，都会通过耦合设备传导到其他网络。

举个最简单的例子：冬季热负荷增加，燃气锅炉的耗气量上升，气网的节点流量重新分配，气源压力发生变化；如果CHP机组参与供热，它的热出力提升往往带来电出力改变，电网潮流随即变化。反过来，如果电网发生故障导致热泵停机，热网热源出力缺口就需要其他热源填补，又回头影响气网用气量。这种“电-气-热”互相牵制的稳态关系，正是电气热能流计算要解决的问题。

因此，项目标题里强调“计及多能耦合”，不是概念上的装饰，而是实实在在的数学建模要求。我们需要把三个物理网络的方程联立起来求解，而不是分别计算后拼个结果。在Matlab里实现这套计算，核心难点也在于此。

1.2 电气热能流计算的本质：从三个独立网络到一个耦合系统

从数学上看，电气热能流计算本质上是求解一组非线性方程组：电网节点功率平衡方程、气网节点流量平衡方程、热网水力与热力方程，以及耦合设备能量守恒方程。状态变量包括电网节点电压幅值和相角、气网节点压力、热网节点供水和回水温度、各耦合设备的内部运行变量。

求解这组方程有两种主流思路。第一种是统一求解法，把电、气、热所有方程和变量拼成一个大雅可比矩阵，用牛顿-拉夫逊法整体迭代。第二种是解耦迭代法，交替求解电网、气网、热网子问题，每次迭代后更新耦合设备变量，直到收敛。

我自己最终选择的是统一求解法。原因很直接：解耦迭代虽然可以复用成熟的潮流程序，但迭代过程中容易出现震荡，尤其当耦合设备比例较高时，收敛性很难保证。统一求解法把所有耦合关系放在一个雅可比矩阵里，物理上更清晰，数学上收敛性也更好，只是代码实现时需要仔细处理变量索引和雅可比矩阵的组装。

2. 系统建模：电气热网络数学模型与耦合环节

2.1 电网潮流模型与牛顿-拉夫逊解法

电网部分采用经典极坐标形式的牛顿-拉夫逊法。每个节点的有功和无功功率方程分别为：

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的电压相角差。节点类型依然分为PQ节点、PV节点和平衡节点。在区域综合能源系统中，耦合设备通常作为电网的PQ节点或PV节点处理，但如果CHP机组采用“以热定电”模式，它的电出力会随热出力变化，这时必须把CHP相关变量作为额外状态量加入统一求解框架，而不能简单当PQ节点。

在Matlab代码里，我会用稀疏矩阵保存导纳矩阵Ybus，节点功率偏差作为残差分量。计算电网残差时，从状态向量中取出V和θ，计算注入功率，再减去给定值。

% 电网功率残差计算 Vn = V .* exp(1j * theta); Iinj = Ybus * Vn; S = Vn .* conj(Iinj); Pcal = real(S); Qcal = -imag(S); % 注意符号约定，这里以注入为正 res_grid = [Pcal(PQPV) - P_spec(PQPV); Qcal(PQ) - Q_spec(PQ)];

这种写法在Matlab里非常直观，调试时也容易对照教科书公式。

2.2 天然气网络水力方程与气潮流计算

天然气网络稳态计算的核心是节点流量平衡。每根管道流量与两端压力相关，最常用的是Weymouth公式：

f_mn = sign(π_m² - π_n²) · C_mn · sqrt(|π_m² - π_n²|)

其中π_m和π_n为节点绝对压力，C_mn为管道常数，由管径、长度、摩擦因数、气体组分决定。sign项保证流量方向从高压流向低压。气网节点类型分为三类：气源节点压力给定、负荷节点流量给定、以及平衡节点（用来平衡全网用气误差）。

计算气潮流时有一点容易被忽略：天然气是可压缩流体，流量对压力呈平方关系，方程非线性程度比电网功率方程还厉害。初值给不好，牛顿法在第一轮就可能发散。我的做法是，所有非气源节点初始压力设为气源压力的0.9倍，并配合阻尼牛顿法来启动。

气网残差写入统一求解器时，还需要考虑压缩机等设备。简化处理中，可以暂时忽略压缩机，或者把压缩机等效为流量约束。对于研究型的能流计算，先忽略压缩机也能验证核心算法是否跑通。

% 气网管道流量计算（Weymouth公式） dp_sq = pi_gas(m) .^ 2 - pi_gas(n) .^ 2; flow(m,n) = sign(dp_sq) .* C(m,n) .* sqrt(abs(dp_sq)); res_gas = node_injection - node_demand - sum(flow, 2);

2.3 热力网络水力-热力联合计算

热网是最容易让人头疼的部分，因为它包含水力工况和热力工况两个层次。水力工况求解各管道质量和节点压力，热力工况求解节点供水和回水温度。

水力模型的基本约束是节点质量流量守恒：流入节点的质量流量等于流出节点的质量流量。另外，每个供热环路的压力降代数和必须为零。简化计算时，可以忽略环网，只考虑枝状管网，这样只需要质量守恒方程就能求出各管段流量。

热力模型的核心是热功率方程：节点热负荷φ_i等于该节点供回水温差乘以质量流量和比热容：

φ_i = C_p · m_i · (T_s,i - T_o,i)

管道温度沿程会下降，考虑环境散热时用指数衰减公式：

T_end = (T_start - T_amb) · exp(-λL / (C_p · m)) + T_amb

多根管道汇入同一个节点时，混合温度按照能量加权平均：

T_mix = ∑(m_j · T_j) / ∑m_j

在热网中，热源和热负荷都表现为对温度或热功率的约束。例如，CHP机组的供水温度固定，而热用户按照热负荷给定。热力方程既有线性也有强非线性，和电网、气网耦合后，对雅可比矩阵的数值精度要求更高。

2.4 耦合设备的耦合方程与状态变量扩展

多能耦合是通过耦合设备体现的，需要为每一类设备建立数学方程。最典型的几类设备如下。

CHP联产机组：输入天然气，输出电和热。简化建模可以写成：

F_gas = P_elec / η_e
Q_heat = F_gas · η_h · LHV_ratio

其中η_e是发电效率，η_h是热回收效率。如果背压式CHP，热电比固定；抽汽式CHP则热电比可以调节，需要增加一个抽汽系数变量。

燃气锅炉：天然气消耗与热出力成正比：

F_gas = Q_heat / (η_gb · LHV)

热泵：消耗电功率产生热功率：

Q_heat = COP · P_elec

P2G设备：消耗电能产生天然气，作为气网的附加气源和电网新增负荷：

F_gas = η_p2g · P_elec

这些方程需要在统一求解时作为额外的残差方程加入，同时把设备内部变量加入状态向量。例如，CHP机组的热出力Q_heat和电出力P_elec，热泵的耗电量P_elec，都需要参与迭代。这样做的优势是耦合设备对三个网络的边界影响可以被同时感知，不会出现解耦迭代时那种“气网算完才想起热网”的滞后效应。

3. Matlab代码实现：从数据到算法的落地细节

3.1 代码整体架构与数据结构设计

写Matlab程序最怕把所有代码堆在一个脚本里，后面根本没法改。我做这个项目时，把代码拆成四层结构：数据定义层、网络模型层、求解器层、结果输出层。

数据定义层用结构体组织，每个物理网络一个struct。电网保存导纳矩阵Ybus、节点编号、负荷、发电机数据；气网保存管道首末端节点编号、管道常数C、节点类型、压力初值；热网保存管道参数、节点负荷、比热容、环境温度；耦合设备保存连接拓扑和效率参数。

data.grid.Ybus = Ybus; data.grid.bus = bus; % bus(:,1)编号, bus(:,2)类型, bus(:,3)V, bus(:,4)theta data.gas.pipe = pipe; % pipe(:,1),pipe(:,2)首末端, pipe(:,3)C常数 data.gas.node = gas_node; % [编号,类型,压力,负荷] data.heat.pipe = heat_pipe; % [首端,末端,长度,直径,散热系数] data.heat.node = heat_node; data.coupling.c = 0; % 依次定义CHP、热泵、锅炉等 data.coupling.chp = chp_data;

数据结构定了，后面的残差函数、雅可比函数、牛顿迭代主循环都能统一处理。我最想强调的一点是：变量索引必须单独写一个函数，负责把全局状态向量x拆成电网相角、电压幅值、气网压力、热网供水和回水温度、耦合变量。千万别在残差函数里靠魔法数字取位置，否则改一个设备类型，索引全部乱套。

function [theta, V, pi_g, Ts, Tr, x_cpl] = split_state(x, idx) theta = x(idx.theta); V = x(idx.V); pi_g = x(idx.pi_g); Ts = x(idx.Ts); Tr = x(idx.Tr); x_cpl = x(idx.cpl); end

3.2 统一求解法与解耦求解法的实现选择

之前说了我选统一求解法。实际编码时，核心就是构建统一的残差向量R(x)和雅可比矩阵J。残差向量按顺序排列：电网有功偏差、电网无功偏差、气网流量偏差、热网水力偏差、热网热力偏差、耦合设备方程偏差。状态向量也按同样顺序排列。

统一求解法的雅可比矩阵维度是所有变量数之和。对一个中小型系统，比如电力5节点、气网4节点、热网3节点、CHP和热泵各1台，状态变量数大概在24个左右，雅可比矩阵24×24，数值差分求解非常快。如果扩展到几十节点，就应该用解析雅可比或者稀疏矩阵。我用Matlab的sparse函数组装J，配合左除\求解线性方程组，速度和内存都很好。

dx = -J_full \ res; x = x + dx;

为了提升迭代稳定性，我加入了阻尼策略：如果新的残差范数没有下降，就减小步长，这有点像Levenberg-Marquardt的思路。实际测试中，阻尼因子取0.5能有效避免震荡。

解耦迭代法也有它的应用场景，比如系统规模很大，气网和热网都有成熟的商业计算软件，外接替代统一求解是合理的。但从学术复现和代码透明性角度，统一求解法更友好，这也是我在博客里推荐统一法的原因。

3.3 关键函数与核心代码片段讲解

残差函数residual(x)是整个程序的心脏。以我搭建的简易算例为例，简要展示代码结构。

function res = residual(x) [theta, V, pi_g, Ts, Tr, x_cpl] = split_state(x, idx); % 1. 电网残差 res_grid = grid_residual(theta, V, data.grid); % 2. 气网残差 res_gas = gas_residual(pi_g, data.gas); % 3. 热网水力残差 res_heat_hyd = heat_hydraulic_residual(data.heat); % 4. 热网站热力残差 res_heat_therm = heat_thermal_residual(Ts, Tr, data.heat); % 5. 耦合设备残差 res_cpl = coupling_residual(x_cpl, theta, V, pi_g, Ts, Tr, data.coupling); res = [res_grid; res_gas; res_heat_hyd; res_heat_therm; res_cpl]; end

耦合设备残差需要把电、气、热变量都串起来。比如CHP机组的残差方程为：

function res_cpl = coupling_residual(x_cpl, V, theta, pi_g, Ts, Tr, cpl) % 提取CHP变量 P_chp = x_cpl(1); Q_heat = x_cpl(2); F_gas = x_cpl(3); % CHP耦合方程 res_cpl(1) = F_gas - P_chp / cpl.chp.eta_e; res_cpl(2) = Q_heat - cpl.chp.eta_h * F_gas * cpl.chp.LHV_ratio; res_cpl(3) = node_gas_injection(1) + F_gas - gas_load(1); % 气网节点1的流量平衡计入F_gas end

这个示例省略了很多细节，但表达了核心逻辑。雅可比矩阵可以用jacobianest之类的数值差分函数自动生成，或者手动写解析式。数值差分的优点是代码量小、不易出错，缺点是步长选择不当会降低精度。我建议第一步先用数值差分跑通，再考虑将关键方程替换为解析导数。

3.4 算例设置与结果可视化

为了验证算法，我构造了一个小型区域综合能源系统：电网采用5节点系统，气网4节点，热网3节点。电网节点1接平衡电源，节点2接CHP机组，节点4接热泵。气网节点1为气源，连接到节点2的负荷和CHP的耗气。热网节点1为CHP热源，节点2为热泵热源，节点3为热负荷。

给定负荷条件后，Matlab程序在初始值接近真实解时，5次迭代内收敛到1e-8精度。如果不调整初值，直接采用电网标幺电压1.0、气网统一压力3 MPa、热网统一温度75℃作为初值，迭代次数大约8次。整体运行时间不到0.2秒，说明统一求解法在小规模系统中性能完全没问题。

结果可视化方面，我会输出四类图：电网节点电压幅值和相角分布图、气网节点压力分布图、热网供回水温度分布图、耦合设备功率流向图。Matlab的bar和stairs画压力分布很直观，quiver可以画能流方向，但实际写代码时我只会输出数值表格，让用户自行决定绘图方式。

4. 常见问题与调试经验实录

4.1 初值设置与收敛性调试

我遇到过最多的收敛性问题是气网初值太差。曾有一段时间，我把所有气网节点压力初值设为2 MPa，但平衡节点实际压力是4.2 MPa，结果牛顿法前两步的雅可比矩阵接近奇异，迭代直接发散。后来我加了一段“预置初值”逻辑：气网非气源节点压力按气源压力乘0.9初始化，热网供水温度统一设为85℃，回水温度统一设为45℃，电网电压标幺值设为1.0。这个初值策略在多个算例里都表现稳定。

如果遇到发散，可以尝试先用小负荷系数运行。例如将所有负荷乘以0.5，求解完成后，再把负荷系数提高到0.75，用上一个解作为初值继续求解，直到全负荷。这种“逐步加载”策略在处理强非线性系统时非常有效。

4.2 耦合设备雅可比矩阵的数值处理

数值雅可比虽然方便，但步长h的选取很敏感。如果h取得太小，差分时出现数值噪声；取得太大，斜率的线性近似误差也会偏大。我的经验是对不同量纲的变量分开设步长：电压和温度用1e-6，压力和相角用1e-5。

更稳妥的做法是解析求导。比如CHP方程的导数非常简单，F对P的导数是1/η_e，Q对F的导数是η_h·LHV_ratio。电网潮流对电压和角度的导数公式在教科书里都有，照着写成稀疏矩阵即可。两者结合，计算速度快一倍以上。

4.3 结果验证与误差分析技巧

代码写完之后，必须验证正确性。我习惯用“退化测试”来检查：将耦合设备出力全部设置为0，此时系统退化为独立电网、气网、热网。分别用传统的单一网络潮流程序计算结果，然后和统一求解结果对比，误差应在1e-6级别。如果偏差大，十有八九是耦合方程或者节点连接关系写错了。

还可以做能量守恒检验：总输入能量（天然气功率加上网侧功率）应该等于总负荷功率加上各网络损耗。热网损耗包括管道散热，气网损耗在绝热稳态下为零，电网损耗就是线路有功损耗。这个测试每次计算完都做，能快速暴露符号方向错误。

以下是我项目中最常遇到的几个问题汇总：

5. 扩展方向与个人实操体会

5.1 从稳态能流走向优化调度

能流计算稳态求解器是综合能源系统优化调度的基石。我之前在后面做优化时，把能流计算封装成一个函数，用遗传算法搜索各设备出力，每个个体通过调用能流求解器来校验可行性和运行代价。这样优化算法不用关心内部网络约束，只需要给定设备出力组合，就能得到系统是否越限、损耗有多大等指标。

如果你后续要做电-气-热联合优化，建议把能流计算模块做得“无状态”化：输入设备出力和负荷，输出节点状态和损耗。这样各种优化算法都可以直接调用，不需要重复建模。

5.2 我踩过的坑与个人建议

最后分享几个实际教训。第一，Matlab代码里不要用大量全局变量，尽量用结构体传递数据，否则调试的时候一个变量被覆盖根本找不出来。第二，一开始不要追求大算例，先搭一个3节点电网+3节点气网+2节点热网的最小系统，跑通后再扩规模。第三，所有物理量统一用国际单位制，气网压力用Pa，流量用kg/s，这个细节能省去大量排错时间。

我现在写能流代码的习惯是，每个算例单独放在一个m文件中，文件开头列出所有假设条件和单位，文件结尾输出关键节点状态。这样三个月后回来看，还能快速捡起来。希望这篇记录能帮你避开我踩过的那些坑。如果你也是做综合能源方向的，欢迎按这个框架动手试一遍，亲手把收敛过程跑出来，才算真正把电气热能流计算吃透。