企业官网建设流程全解析

在去年做园区综合能源规划评审的时候，甲方问了我一个挺扎心的问题：“你们给出的电负荷、热负荷、气负荷预测，是三个团队各算各的，还是放在同一个框架里算的？”我承认当时多少有点心虚，因为行业里确实习惯电力潮流归电力、热网水力归热网、天然气管网归燃气管网，最后拼成一张系统图完事。但近几年不行了——燃气轮机同时供电气和热，电锅炉在低谷时段批量启动，P2G设备甚至能反向影响天然气管网的压力。CHP机组的发电计划，直接决定热网有多少热量、气网抽多少气；电负荷波动会通过电锅炉传导到热网；天然气管网压力波动又能限制燃机出力上限。面对这种强耦合场景，还在用解耦方式做区域综合能源系统的电气热能流计算，结果可信度真的要打问号。这篇文章我就从物理模型、耦合设备建模、Matlab求解器实现、算例验证到实际调试踩坑，把这一整套流程完整讲一遍，适合正在做综合能源规划、运行优化或者课题研究，需要自己手写能流计算代码的同学。

1. 为什么电气热能流必须放在同一张图里迭代求解

1.1 “各算各的”到底会差多少

先举一个很典型的场景。某个工业园区里有一台额定电功率3 MW的CHP机组，热电比1.2，也就是说发3 MW电的同时能回收3.6 MW热量。如果按传统的解耦思路，电力潮流计算时会把CHP当成一个固定出力为3 MW的发电机节点，热力计算时把热负荷5 MW当成固定节点负荷，天然气计算时把耗气量当成固定值。问题是，实际运行中CHP大概率是“以热定电”——热网需要多少热，CHP就发多少电。热负荷变化时，注入配电网的电功率跟着变，配电网的电压分布也随之改变。反过来，如果电网侧电压偏低、调度要求CHP降出力，热网瞬时缺的热量就得靠燃气锅炉顶上，燃气锅炉的耗气量又叠加到气网负荷上。这几条链路互相咬合，没有任何一个子系统可以独立求得准确解。

我做过一个不算复杂的测试：一个IEEE 33节点配电网、6节点气网、6节点热网组成的小型综合能源系统，里面只有一台CHP和一台电锅炉。解耦计算和耦合计算结果对比下来，配电网节点电压最大偏差接近1%，气网节点压力最大偏差超过4%，热网节点温度偏差约0.5℃。如果不考虑P2G，这个偏差还不算离谱；一旦加上P2G，气网压力偏差可以直接拉到15%以上，甚至出现方向性错误——比如解耦算出来某节点压力正常，耦合算出来该节点已经跌破最低供气压力。

1.2 耦合的数学本质

从数学上看，三个子系统的稳态模型可以统一写成一组非线性代数方程：

F_elec(x_elec, x_gas, x_heat) = 0 F_gas(x_gas, x_elec, x_heat) = 0 F_heat(x_heat, x_elec, x_gas) = 0

解耦计算等价于把x_gas和x_heat固定在某组常数上去解F_elec，把x_elec和x_heat固定去解F_gas，依次类推。这在数学上等同于完全忽略了交叉偏导项∂F_elec/∂x_gas、∂F_gas/∂x_elec这些非对角块。当耦合设备渗透率低的时候，忽略交叉偏导带来的误差也许在工程接受范围内；但耦合度一旦上去，这些交叉偏导不仅数值可观，而且随着运行点剧烈变化，解耦结果就完全不可控了。

所以，真正要做“计及多能耦合”的能流计算，技术上就两条路：一是把所有方程联立起来做统一牛顿法求解，二是分解迭代但反复交换耦合变量直到收敛。两者没有绝对的优劣，实践里怎么选，我后面会详细讲。

2. 从物理方程到可计算的残差形式

2.1 电力网络潮流方程

电力部分继续沿用经典极坐标牛顿-拉夫逊潮流模型。对每个PQ节点，写有功和无功两个功率平衡残差：

ΔP_i = P_i^sp - V_i * Σ(V_j * (G_ij*cosθ_ij + B_ij*sinθ_ij)) = 0 ΔQ_i = Q_i^sp - V_i * Σ(V_j * (G_ij*sinθ_ij - B_ij*cosθ_ij)) = 0

PV节点只写有功残差，平衡节点什么都不写，状态变量是除平衡节点外所有节点的电压相角和PQ节点的电压幅值。这块大家都很熟悉，不多展开。需要特别注意的是，节点注入功率P_i^sp里面，如果接入了CHP或者电锅炉，对应的功率项不能预先写成常数。比如CHP节点，P_i^sp = P_load - P_chp，而P_chp本身是耦合设备模型算出来的、随迭代变化的值。忽略这一点，代码写得再漂亮，本质上还是解耦的。

2.2 天然气网络稳态流量方程

稳态天然气网络用节点气压作为状态变量，管段流量用Weymouth方程描述。中高压管网中，管段流量和两端气压的关系近似写成：

f_ij = sign(p_i - p_j) * C_ij * sqrt(|p_i² - p_j²|)

C_ij是管段综合系数，包含管径、长度、温度、压缩因子等参数。注意这里必须带sign符号项，否则管内流量方向固定死了，实际管网中气流方向是随压力分布变化的，忽略sign会导致负值开根号直接报错。

每个节点的流量平衡残差写成：

ΔF_i = Σf_in - Σf_out + F_source - F_load - F_device = 0

其中F_device就是从该节点取气的耦合设备耗气量，比如CHP、燃气锅炉、P2G的耗气（P2G是产气，符号相反）。天然气负荷里，一部分是固定负荷，这部分提前给定；另一部分是随耦合设备出力变化的“动态负荷”，必须在每次迭代里更新。

压缩机模型这里只说一种常用的简化处理：给定升压比或者给定出口压力，压缩机的自耗气量按经验公式估算。如果一开始就上精细的动态压缩机模型，代码复杂度和调试难度会陡增，而且对稳态能流结果的影响通常有限，建议第一版代码先做简化。

2.3 热力网络的水力-热力联合方程

热力网络是所有子系统里最容易写错、也最容易数值翻车的一块。它本身包含两个层次：水力方程和热力方程，两者通过管道流量耦合在一起。

水力方程，一是节点流量连续性：

A * m = q

A是节点-管道关联矩阵，m是管道流量向量，q是节点净注入流量；二是管道压降方程：

Δp_drop = K * m * |m|

K是管道阻力系数，流量与压降呈平方关系。在供热管网里，循环泵提供的压头平衡所有管道压降。

热力方程这边，首先每个热负荷节点的热功率与供回水温差有关：

Φ_load = C_p * m_node * (T_s - T_o)

m_node是流经该负荷节点的流量，T_s和T_o分别是供水温度和回水温度。然后是管道沿程温降，考虑环境散热时管道末端温度不等于首端温度：

T_out = (T_in - T_a) * exp(-λ*L / (C_p*m)) + T_a

最后是节点混合温度方程，多条支路的热水在节点汇合后温度是流量加权平均：

Σ(m_in * T_in) = T_mix * Σ(m_in)

热力网络里我最想强调的一个实操点：状态变量不应该把“每条管道流量、每个节点供回水温度”全塞进去硬解，因为节点混合温度实际上是由代数关系直接决定的。更稳的做法是，把节点混合前后两段变量分离，通过消元把混合温度表达成支路温度的函数，只保留独立变量进雅可比矩阵。直接全塞进去，雅可比矩阵里会出现大量线性相关的行，轻则收敛慢，重则矩阵奇异直接报错。

2.4 耦合设备：CHP、电锅炉、燃气锅炉与P2G

耦合设备的建模决定了整个系统的耦合强度，也是这篇代码的灵魂。我的做法是统一成“功率-能耗转换器”的视角来看：

CHP机组是最经典的耦合元件。简化的背压式机组模型用热电比描述：

Φ_recover = R_h * P_chp F_gas = P_chp / (η_e * LHV)

P_chp是发电机出力，Φ_recover是回收的热功率，R_h是热电比，η_e是发电效率，LHV是天然气低位热值。这里的热电比在中小型机组上通常在1.0～1.5之间，选型时拿厂家样本数据即可。注意，CHP的“以热定电”和“以电定热”两种运行方式在代码里对应不同的变量约束，实现时要留一个控制模式开关。

电锅炉模型最简单，电功率转热功率：

Φ_eb = η_eb * P_eb

效率通常在0.95以上，基本可以当成纯电阻加热处理。燃气锅炉也一样，天然气化学能转热能：

Φ_gb = η_gb * F_gas * LHV

P2G要稍微费点心思。电解水制氢的效率决定了电转氢的量：

F_H2 = η_p2g * P_p2g / LHV_H2

但氢气能不能直接注入天然气网络，取决于掺氢比例约束。实际工程中掺氢比例通常限制在5%～20%（视管网材质和终端燃气具而定），如果代码里直接把氢气产量全部注入气网，算出来的气压可能明显偏离真实约束。稳妥做法是加一个简单的上限约束，比如注入气网节点的氢气等效流量不超过该节点总流量的10%。

3. Matlab代码总体架构：统一牛顿法为主、分解法兜底

3.1 两种求解框架的取舍

统一求解法把电气热三个子系统的状态变量堆成一个长向量，所有方程堆成一个长残差向量，直接做一次牛顿迭代：

x = [Va; Vm; p_gas; m_pipe; Ts; Tr] F = [F_elec; F_gas; F_hydraulic; F_thermal; F_coupling] x_{k+1} = x_k - J(x_k) \ F(x_k)

这个方案的优点是严格的二阶收敛，交叉耦合项完整保留，是我在强耦合场景下的首选。缺点是编程工程量确实大，尤其是雅可比矩阵的组装，四个子系统加耦合块的偏导要全手动推导和编码，初期调试周期不短。

分解法则是按“电力→天然气→热力”顺序轮流迭代，每个子系统用成熟求解器，子系统之间互相交换耦合变量，直到整体收敛。模块化好、代码复用率高，但强耦合时经常出现迭代震荡甚至发散，要在各子系统之间加松弛因子：

x_new = x_old + α*(x_computed - x_old)

α取0.5～0.8比较常见。实测下来，分解法在耦合度中等（单台CHP）的场景还能勉强工作；一旦P2G和电锅炉都上来，收敛速度会急剧恶化，迭代次数从几次飙到几十次。

我自己的工程建议是两阶段策略：先用分解法迭代5～8轮，得到一个合理的中间初值，然后切换统一牛顿法精解。这个做法既避免了统一法冷启动初值太差导致发散的问题，又能在最后阶段拿到严格耦合解。代码里控制这个切换非常简单，一个if判断即可。

3.2 代码模块划分与数据结构

我自己写这类代码的习惯是，用struct组织全系统数据，不用class。struct的好处是轻量、透明，课题组里其他人接手时一眼就能看明白字段含义，改起来也无压力。整体模块划分如下：

01_load_case.m 读取并检查网络数据 02_build_index.m 生成状态变量全局索引 03_init_guess.m 设定迭代初值 04_ies_residual.m 计算残差向量 05_ies_jacobian.m 计算雅可比矩阵 06_newton_solver.m 牛顿迭代主循环 07_post_process.m 结果统计与可视化 run_ies_pf.m 主脚本，一键运行

这种模块化划分的核心思路是让“残差计算”和“雅可比计算”完全分离，后续想升级模型，比如换更精细的压缩机模型，只需要改04和05两个文件，不用动主循环。

数据结构大体长这样：

% ---- 电气热综合能源系统数据定义 ---- ies = struct(); % 电力系统 ies.elec.baseMVA = 10; ies.elec.bus = [bus_i, type, Pd, Qd, ...]; % 节点数据 ies.elec.branch = [fbus, tbus, r, x, b, ...]; % 支路数据 % 天然气系统 ies.gas.baseP = 1.0; % 基准压力 MPa ies.gas.node = [node_i, load, source, ...]; % 节点数据 ies.gas.pipe = [fnode, tnode, C, ...]; % 管段数据 % 热力系统 ies.heat.baseMW = 10; ies.heat.node = [node_i, phi_load, ...]; % 节点数据 ies.heat.pipe = [fnode, tnode, K, L, d, ...]; % 管道数据 % 耦合设备 ies.coup(1).type = 'chp'; ies.coup(1).ebus = 5; % 接入电网节点 ies.coup(1).gnode = 2; % 接入气网节点 ies.coup(1).hnode = 2; % 接入热网节点 ies.coup(1).P = 3.0; % 额定电功率 MW ies.coup(1).Rh = 1.2; % 热电比 ies.coup(1).eta_e = 0.42; % 发电效率

3.3 状态变量全局索引：代码里最容易乱的地方

统一求解法的代码里最核心也最容易写乱的地方就是状态变量的全局索引。我强烈建议在02_build_index.m里把所有变量的起始位置一次性算好，存成一张索引表，后续残差和雅可比组装都从索引表取位置，不要在函数里反复计算。

function idx = build_index(ies) n_elec_bus = size(ies.elec.bus, 1); n_gas_node = size(ies.gas.node, 1); n_heat_node = size(ies.heat.node, 1); n_pipe = size(ies.heat.pipe, 1); % 电力状态变量：除平衡节点外所有Va + 所有PQ节点Vm idx.Va = (1:n_elec_bus-1)'; idx.Vm = (n_elec_bus : n_elec_bus + n_pq - 1)'; % 天然气状态变量：全部节点压力 idx.p_gas = (max(idx.Vm)+1 : max(idx.Vm)+n_gas_node)'; % 热力状态变量：管道流量 + 节点供水温度 + 节点回水温度（经代数消元后） idx.m_pipe = ...; idx.Ts = ...; idx.Tr = ...; % 耦合设备的状态变量（如有需要额外引入） idx.coup = ...; % 残差向量的索引与之一一对应 idx.F = idx; % 实际工程中常用相同索引结构 end

初值设定同样在03_init_guess.m里集中处理。电压相角统一设0，电压幅值设1.0。气网节点压力统一设成气源出口压力附近，比如0.8 MPa。热力部分是最需要谨慎的：供水温度设设计值90℃，回水温度设设计值50℃。千万不要把回水温度初值设成环境温度20℃，那样供回水温差巨大，热功率方程里残差极小但雅可比元素也很小，数值上会产生严重的尺度问题。

4. 核心求解器实现：组装雅可比矩阵的关键细节

4.1 残差函数怎么写才高效

残差函数是整个求解器的核心，我见过很多初学者把它写得又慢又难调试。最基础的几个原则，首先是所有映射关系在build_index阶段就固化好，残差函数里只查表；其次要尽量矢量化，能用矩阵运算就绝不用for循环。比如电力残差的计算，从节点导纳矩阵直接切片就能一次算出全部节点的注入功率，不需要逐节点循环。

function F = ies_residual(x, ies, idx) F = zeros(size(x)); % 1. 取出各子系统的状态变量 Va = x(idx.Va); Vm = x(idx.Vm); p_gas = x(idx.p_gas); m_pipe = x(idx.m_pipe); Ts = x(idx.Ts); Tr = x(idx.Tr); % 2. 电力残差（此处仅示意核心逻辑） V = Vm .* exp(1j * Va); % 构建复电压向量（需完整还原全节点） S_calc = V .* conj(Ybus * V); % 节点注入复功率 F(idx.F_P) = real(S_calc) - P_spec; F(idx.F_Q) = imag(S_calc) - Q_spec; % 3. 天然气残差 % 根据Weymouth方程计算各管段流量，再按节点累加 % 耦合设备的耗气量由当前状态变量计算后计入F_load_device % 4. 热力残差 % 水力方程、节点温度平衡方程分别组装 % 5. 耦合设备方程 % CHP、电锅炉等自身的等式约束 end

写残差函数有个非常实用的调试技巧：先单独调用一次能流计算程序，把某个子系统节点上已知的解析结果代入残差函数，检查该部分残差是否为零向量。如果该部分不为零，说明这一块的公式或索引写错了。按子系统逐个验证，能把调试时间缩短一半以上。

4.2 雅可比矩阵的分块组装与验证

雅可比矩阵是统一求解法里工作量最大、最容易被细节坑到的部分。整体结构是一个分块矩阵：

电力变量 气网变量 热力变量 电力方程 J_EE J_EG J_EH 气网方程 J_GE J_GG J_GH 热力方程 J_HE J_HG J_HH

其中对角块J_EE、J_GG、J_HH分别是各子系统自身的偏导矩阵，非对角块J_EG、J_GE这些就是耦合贡献，全部来自耦合设备。比如CHP接入后，电力方程的残差对气网压力变量没有直接偏导（因为CHP的耗气量影响的是气网方程），但气网方程的残差对电力变量有偏导（因为CHP耗气量是发电出力的函数）。组装时用稀疏矩阵累加即可，不需要等一个稠密矩阵。

function J = ies_jacobian(x, ies, idx) N = length(x); J = sparse(N, N); % 电力块 J(idx.F_P, idx.Va) = dFp_dVa; J(idx.F_P, idx.Vm) = dFp_dVm; % ... 其他子系统块同理 % 耦合块 % CHP：气网方程对电出力的偏导 J(idx.F_gas, idx.Va) = J(idx.F_gas, idx.Va) + dFgas_dVa_chp; % P2G：气网方程对电功率的偏导（负号，因为P2G消耗电产生气） % 电锅炉：热网方程对电功率的偏导 end

关于雅可比矩阵的验证，我最推荐先用数值差分做一次对照：

J_num = zeros(N, N); h = 1e-6; for i = 1:N xp = x; xp(i) = xp(i) + h; xm = x; xm(i) = xm(i) - h; J_num(:, i) = (F(xp) - F(xm)) / (2*h); end

数值雅可比不用跑很多次，就取一个典型运行点，和解析雅可比做差，看到最大偏差在1e-5量级即可认为解析推导正确。之后正式求解就用解析雅可比，速度快且精度稳定。

还有一个工程建议：如果不想手推所有偏导，第一版可以全程用数值雅可比跑通逻辑，跑通后再逐步替换解析块。以前端工程思维来看，这属于典型的“先让它跑起来，再优化”策略。

4.3 迭代收敛判据的坑

收敛判据我建议同时盯两个指标：残差范数max(|F|)和修正量范数max(|dx|)。只盯一个指标会出问题，我实际遇到过残差已经降到1e-10量级、但状态变量还在小幅变飘的情况——这通常是因为方程之间存在弱线性相关，多解或近多解导致的。所以代码里两个判据同时满足才认为收敛：

tol_F = 1e-8; tol_dx = 1e-10; for iter = 1:max_iter F = ies_residual(x, ies, idx); J = ies_jacobian(x, ies, idx); if max(abs(F)) < tol_F % 记录残差已达标 flag_F_ok = true; else flag_F_ok = false; end dx = -J \ F; x = x + dx; if max(abs(dx)) < tol_dx && flag_F_ok disp(['收敛于第 ', num2str(iter), ' 次迭代']); break; end % 不收敛时输出当前状态，方便追踪 if mod(iter, 5) == 0 fprintf('iter=%d, maxF=%.2e, maxdx=%.2e\n', iter, max(abs(F)), max(abs(dx))); end end

另外，不同物理量纲的变量混在一个向量里，x的修正量范数会被大数值量级的变量主导。比如节点温度是几十的量级，而电压幅值标幺值是1左右，直接比max(|dx|)相当于默认温度修正更重要。严谨的写法是加“伪标幺化”——把每个分量的修正量除以该变量的基准值再取范数，这样各变量能公平参与收敛判断。

5. 算例设计与结果验证：怎么证明代码算对了

5.1 一个可复现的小型电气热耦合测试系统

我自建了一套测试系统，数据原型都来自公开文献改出来的，方便复现。电网采用IEEE 33节点配电网，经过改造后系统包括一个33节点配电网、一个6节点天然气网和一个6节点热力网。

核心参数如下表：

耦合设备接入位置和参数：

这套系统规模不大，但麻雀虽小五脏俱全，电气热三个网络的耦合要素（CHP、电锅炉、燃气锅炉）全部覆盖。如果想测试P2G，可以在节点12加一组500 kW的P2G设备，气网侧接入节点5。

5.2 耦合与解耦的结果对比

我做了两组计算，一组是严格耦合求解，另一组是解耦计算（把CHP电出力、耗气量、热出力全部预设为固定值）。结果对比如下：

这个结果说明了什么？在只含一台CHP和一台电锅炉的“中等耦合”场景下，解耦计算电压偏差1%左右也许还能忍受，但气网压力偏差4.3%已经能影响管网的运行安全性判断了——如果解耦算出来节点压力在安全范围内，实际耦合运行却可能低于最低供气压力。而P2G加入后，偏差量级会进一步放大。能流计算如果连压力是否越限都判断错，后续的优化调度、可靠性评估全都建立在错误的基础上。

5.3 收敛性分析和异常表现

同一套数据，我用三种策略分别测试：

• 统一牛顿法冷启动：6次迭代收敛到1e-8，过程平稳。
• 分解法：17次迭代勉强收敛，中间出现两次明显的震荡，如果耦合设备再多一点，大概率直接发散。
• 两阶段法（先分解5轮再切统一法）：总计7次迭代收敛，比冷启动统一法更稳定。

冷启动统一牛顿法之所以也能收敛，是因为这套测试系统的初值给得比较合理，尤其是热力网络回水温度设了50℃而不是20℃。初值一差，收敛性和稳定性立刻变脸。调初值的过程我总结成一句话：电网初值随便给，气网初值靠气源，热网初值靠设计——热力部分最敏感。

求解过程中还有一个常见异常：迭代发散的征兆通常是max(F)先降后升，反复震荡。遇到这种情况先别急着改阻尼参数，优先检查雅可比矩阵的正确性，尤其检查耦合设备相关的非对角块符号。我调试时踩过一次最大的坑就是CHP耗气量对气网方程贡献的符号写反了，导致整个系统变成一个“正反馈”回路，怎么迭代都发散。

6. 工程实践中的单位、初值与数值病态

6.1 单位制统一是第一优先级

电气热混合计算里单位制是最大的隐性坑。电力系统通常用标幺值，气网用MPa、m³/h，热网用MW、kg/s、℃。三个子系统的残差和雅可比元素混合在一起时，数值量级可以差出10^6倍以上。直接组装成一个大矩阵，条件数会非常难看，求解精度大幅下降。

我的做法是各子系统内部保留自己的习惯单位，但在组装全局残差和雅可比之前，先把所有残差除以本子系统的基准值，做一个“伪标幺化”。比如气体残差除以该气网节点的基准流量（kg/s），热力残差除以基准热功率（MW）。这样全局残差各分量都在同一数量级，后续迭代收敛判据也更合理。

6.2 阻尼策略和数值病态的应对

即使初值给得不错，强耦合系统在牛顿迭代中途也可能迈出过大的步子。我的代码里加了一个简单的回溯线搜索，当||F(x + αdx)|| > (1-ρ)||F(x)||时，把步长α减半重试。ρ取0.1～0.3，加这个阻尼后整体迭代次数会多几次，但换来的是大幅提升的稳定性，在工程里完全值得。

雅可比矩阵用sparse存储、用左除求解，千万别用矩阵求逆。另外，如果条件数实在太大，可以在求解前做行列缩放、改善条件数，这个操作对大型系统效果明显。

6.3 模型简化的边界要清楚

最后聊聊模型简化的尺度。我见过不少新手一上来就堆复杂模型——压缩机动态方程、管道暂态过程、P2G多步化学反应全塞进去，代码跑不动不说，稳态能流的本质反而被淹没了。我的建议是第一版代码里做这样几个简化：

• 天然气压缩机用定升压比模型，暂不引入变转速、变工况性能曲线。
• 热力管道温降先用平均温度近似，管段较短时可忽略温降；只有长输管道才计入λL/(C_p m)项。
• P2G产氢按效率直接折算，但必须带掺氢比例约束。
• 供热管网忽略网络热惯性，稳态工况下这不影响能流结果。

这些简化不是偷懒，而是稳态能流计算本身就要回答的核心问题——在给定的运行条件下，电气热三个网络各自的状态量是多少。动态特性、变工况性能是后续动态仿真和优化研究的范畴，混在一起只会让问题失去焦点。把稳态能流的计算逻辑先走通，后面再逐步细化，这条路比我一开始追求模型完备要高效得多。