企业官网建设流程全解析

1. 从一次热电联产联调事故说起：为什么单算三种网络会翻车

去年我做了一个园区级综合能源系统的规划评估项目，系统里有中压配电网、中压燃气管网、热水供热管网，还装了两台燃气轮机组、三台热泵、两台燃气锅炉。评审前我想省事，用三个成熟的独立程序分别算电、气、热：电网潮流程序假设燃气轮机满发，气网程序按最大供气负荷定容，热网程序把CHP热出力当成定值。结果显示电网无过载、气网容量够、热网温度达标，看起来完美。

现场边界的几条联调管线焊完，DCS系统一通上，问题全暴露了：某一回线夏季高负荷时母线电压跌到0.90 pu以下，气网末端压力比运行规程低了15%，热网最远端用户供水温度只有62℃，设计值可是75℃。三个独立程序算出来的"最优方案"在真实耦合运行下根本兜不住。

问题出在哪？就在于这三个网络根本不应该分开算。电网的燃气轮机烧的是天然气，气网的压力波动会直接影响发电出力；热泵要从电网取电才能向热网供热，电网电压跌了热泵制热效率跟着掉；燃气锅炉向热网供热的同时也在气网里抽气。任何一个网络只要有波动，都会通过耦合设备传导出去。这就是标题里"计及多能耦合"要解决的核心问题——电气热能流计算，本质上是在联立求解三个网络方程和一个耦合设备方程组的数学问题。

所以这篇文章，我把这套计算的物理模型、数学方程、Matlab实现思路以及调试经验从头到尾串一遍。适合正在做综合能源方向的研究生、搞园区能源规划的工程师，以及想从独立潮流过渡到多能流计算但找不到完整参考的人。

2. 电气热三种网络的数学模型：先弄清楚每个网络在解什么未知量

2.1 电力网络：牛拉法潮流是基础，但耦合节点的处理方式不同

电力网络这一步没什么新鲜玩意儿，就是常规牛顿-拉夫逊潮流计算。节点分三类：PQ节点（给定有功和无功）、PV节点（给定有功和电压幅值）、平衡节点（给定电压幅值和相角，吸收全网功率差额）。状态变量是各节点电压幅值V和相角θ，极坐标下的功率不平衡方程：

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

但在综合能源系统里要注意，耦合节点不再是简单的PQ/PV。燃气轮机节点在电网模型里通常当成PV节点，因为原动机有功出力可是由气网侧的天然气供应量决定的，不是独立给定的。这就在不同网络之间引入了相互制约的耦合变量，需要把气网的气压和CHP的气耗量纳入整个联立方程组求解框架。热泵节点更特殊——它既是电网的功率负荷，又是热网的固定热出力源，电功率需求本身由热网侧的制热需求决定。

2.2 天然气网络：Weymouth方程和节点气压，很多人第一次栽在这

天然气网络和电网在数学上很像，都是节点流量平衡加支路特性方程。高压输气管道通常用Weymouth方程描述稳态气流：

F_ij^2 = C_ij^2 * (p_i^2 - p_j^2)

其中F_ij是管道流量（kg/s或m³/h），p是节点绝对压力，C_ij是和管道直径、长度、摩擦系数相关的常数。变量上，流量F相当于电流，压力p相当于电压，但注意是压力的平方差。每个节点满足流量平衡：

Σ F_in - Σ F_out + F_source - F_load = 0

节点类型分两类：压力给定节点（类似电网的平衡节点，比如气源/门站节点）和流量给定节点（类似PQ节点，气负荷节点）。状态变量是全网节点压力，如果初始压力给得离谱，方程的平方项会让雅可比矩阵直接病态。

压缩机的建模是热电联算里最容易漏的一块。电驱动压缩机要从电网取电，燃气驱动压缩机要从气网取气，压缩机耗气/耗电功率普遍按以下经验公式折算：

P_comp = k * F * (p_out/p_in)^m - 1

在综合能源系统里如果压缩机是电驱动的，那它就成了电网和气网之间的又一个耦合环节，不建这个模型，气网压力可能算得虚高。

2.3 热力网络的"水力-热力"双层模型

热网比电网和气网都麻烦，因为它有水力过程（水泵、流量、压力）和热力过程（供回水温度、热负荷、管道散热）两个子模型，而且这俩是嵌套耦合的——水力工况决定各管段流量分布，流量分布又决定节点混合温度和管道散热，反过来温度又影响热源侧的出力需求。

水力模型就是管网的基尔霍夫定律。对每个节点流量平衡：

Σ Q_in - Σ Q_out - Q_load = 0

每个环路压降平衡：

Σ Δp_loop = 0

管段压降和流量的关系通常用Darcy-Weisbach或简化幂函数：Δp = R * Q^1.75 到 Q^2。

热力模型是另一套方程。节点供水温度T_s和回水温度T_r是两个核心状态变量，管道末端温度按下式计算：

T_end = T_ground + (T_start - T_ground) * exp(-λ * L / (c_p * Q))

节点混合温度则根据各支路流量加权平均。热负荷节点的用热功率由下式确定：

P_heat = Q * c_p * (T_s - T_r)

注意这里的Q是流过负荷节点的质量流量，T_s和T_r是供回水温度。热网稳态计算中，通常给定热负荷和部分节点的供回水温度，求解全网流量和温度分布。实际算的时候，水力迭代和热力迭代要交替进行，我先固定温度算流量，再用新流量更新温度，反复迭代直到两者都收敛。

2.4 三种网络的时间常数差异与稳态假设

电网的动态过程是毫秒级，气网是秒到分钟级，热网的热惯性是小时级甚至更长。要做稳态多能流计算，就得接受一个前提：假定系统运行在某一稳定工况，忽略动态过程。这在工程规划阶段完全够用——我们算的是"在某个典型运行方式下系统是否能满足供需平衡"，不是暂态稳定性。但理解这个时间尺度差异很重要，它决定了迭代策略：热网通常最后收敛，因为它的响应最慢，初值误差可能要到最外层循环才体现出来。

3. 耦合设备建模：算例能不能收敛，一大半看这几个环节

3.1 CHP机组：电热双输出的"能量路由器"

CHP是多能流计算里最核心的耦合设备。燃气轮机的燃料气通过燃烧做功发电，余热通过换热器供热。理想化模型可以写成：

F_gas * LHV = P_elec / η_e = P_heat / η_h

其中F_gas是燃料消耗量，LHV是天然气低位热值，η_e是发电效率，η_h是热回收效率。实际建模更常用的是电力-热力可行域约束：背压式机组热电比恒定，PQ关系是一个固定斜率；抽凝式机组的热出力可以在一定范围内调节，电出力也会随之变化。多能流计算中通常把CHP当成气网的一个节点负荷（消耗天然气），电网的一个有功源（发出电功率），热网的一个热源（输出热功率），三者通过上面的方程绑定在一起。

建这个模型时要格外小心变量的正方向和量纲。气网里CHP耗气量是负数源，电网里是正功率注入，热网里是热源注入，三个网络各自的方程对同一个CHP的变量符号约定完全不一样。我在实现时踩的坑是电网变量用标幺值，气网用m³/h，热网用kW，导致耦合矩阵里雅可比元素数量级差了好几个数量级，一度收敛很慢。后来统一改写各网络计算函数，让它们都返回国际单位制下的SI核心量纲变量，再在各自函数内部转换成标幺值计算，联立求解才稳定下来。

3.2 热泵和燃气锅炉：两个好算但容易被忽略的环节

热泵模型简单：

Q_heat = COP * P_elec

COP（制热性能系数）通常取固定值或随工况变化的经验曲线。热泵在电网是负荷、在热网是热源，它的电功率需求在联立方程组里不是独立给定的，而是由热网侧需要它输出的热量除以COP反算出来。这个"反算"逻辑要写清楚，不然容易把功率方向弄反。

燃气锅炉更简单：

P_heat = η_b * F_gas * LHV

本质上是天然气流量和热出力之间的线性/准线性关系。但要注意，锅炉的烟气往往需要引风机、循环泵，这些设备本身还要耗电，严格算的话还得把这部分电耗算进电网负荷。如果是做高精度校核计算这些是不能忽略的，但一般规划阶段可以归到站用电的固定比例里。

3.3 能量枢纽耦合矩阵：把不同设备抽象成统一形式

多能流计算发展到后面，大家喜欢用能量枢纽（Energy Hub）的方式统一描述耦合设备。核心思想就是把一个站点的输入（电、天然气、太阳能等）通过一个耦合矩阵C映射到输出（电、热、冷等）：

[ P_e_out ] [ C_ee C_eh ] [ P_e_in ] [ P_h_out ] = [ C_he C_hh ] * [ P_g_in ]

矩阵中的元素表示耦合设备的效率或分配系数。CHP在能量枢纽里表现为C_he（天然气转热）和C_ee（天然气转电）同时非零；热泵表现为C_eh（电转热）。用能量枢纽建模的好处是，不管站点里装了多少台设备，都可以先综合成汇总矩阵再接入整体计算程序，代码模块复用性很好。我在Matlab里把每个站点的耦合矩阵定义成struct存起来，再统一组装到全网方程组里，调试时只看矩阵元素就能发现问题出在哪个设备上。

4. Matlab实现路线：统一求解和分解迭代我都试过，各有利弊

4.1 统一求解法：大雅可比矩阵怎么拼出来

统一求解法是把电网、气网、热网的所有方程堆成一个大的非线性方程组F(X)=0，牛顿-拉夫逊迭代一次同时修正所有状态变量。整个系统的状态变量向量长这样：

X = [V_e; θ_e; p_g; T_s; T_r; Q_flow]

其中V_e和θ_e是电网每个非平衡节点的电压幅值和相角，p_g是气网非压力给定节点的压力，T_s和T_r是热网节点的供回水温度，Q_flow是热网各管段流量（或者不显式列流量，用水力方程直接求解）。

雅可比矩阵是分块结构：

J = [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的存在意味着气网的气压会影响电网注入功率，导致J_ge非零；热泵的存在意味着热网负荷会影响电网功率，导致J_eh非零。

统一求解最大的优点是强耦合工况下收敛性好，电网气网热网的变化在一个迭代步里同时传递，不会出现分解迭代那种"一个网络已经算完了另一个还没反应"的滞后问题。缺点也很明显：编程量大、雅可比矩阵维度大（一个中等规模系统动辄几百上千维）、初值要求高、程序耦合度高没法复用现有代码。调试时我通常把JOFF标志位打开，单独用解析雅可比和数值差分雅可比做对比，交叉导数错误的概率很高。

4.2 分解迭代法：模块化清晰，但收敛条件要小心

分解迭代法是网与网之间用松弛法迭代。基本流程：

1. 给定CHP、热泵等耦合设备出力初值。
2. 把CHP电出力作为电网注入功率，把热泵当成电网负荷，单独求解电网潮流。
3. 把CHP耗气量、锅炉耗气量作为气网负荷，求解气网。
4. 把CHP热出力、锅炉热出力、热泵热出力作为热源，求解热网。
5. 根据气网算出的CHP可用气量，重新修正CHP的电出力；根据热网算出的热负荷需求，重新修正热泵和锅炉的热出力。
6. 返回第二步，直到耦合变量前后两次迭代的差值小于允许误差。

分解迭代的优点是各网络可以分别用成熟程序求解，程序组织清晰。缺点是收敛性完全取决于耦合强度——如果耦合设备功率占比大、或者热电比很高，简单的高斯迭代就会发散。我在一个热电比0.95的算例里第一次跑分解迭代，出现了典型的往返震荡：CHP电出力一会儿0.8 MW一会儿1.3 MW，气网压力也跟着来回跳，最后只能降低松弛因子到0.4才勉强压住。

两种方法的关键对比如下：

我个人建议：如果做单次规划校验、不需要频繁改系统拓扑，用统一求解法；如果要做随机多场景分析、需要在不同案例间频繁切换网络参数，用分解迭代法更现实，配上阻尼机制能稳住大半场景。

4.3 我在Matlab里搭建的实现框架和关键代码

我基于统一求解法实现了一套代码，核心框架分三层：网络参数层、耦合设备层、求解核心层。网络参数层用struct存电网的bus/branch数据、气网的node/pipe数据、热网的node/pipe数据；耦合设备层定义CHP、热泵、锅炉、压缩机各自的参数、连接关系和模型函数；求解核心层负责组装雅可比矩阵、迭代求解、结果整理。

主程序流程可以概括为：

%% 主流程 % 1. 定义网络和设备参数（示例） grid = define_electric_network(); % 电网bus/branch gas = define_gas_network(); % 气网node/pipe heat = define_heat_network(); % 热网node/pipe units = define_coupling_units(); % CHP/热泵/锅炉/压缩机 % 2. 初始化状态变量 X0 = initialization(grid, gas, heat, units); % 3. 牛顿-拉夫逊迭代 for iter = 1:max_iter [F, J] = build_residual_and_jacobian(X0, grid, gas, heat, units); delta_X = -J \ F; % 阻尼处理 alpha = damping_factor(F, delta_X, @residual_norm); X0 = X0 + alpha * delta_X; if norm(F, inf) < tol break; end end

组装残差和雅可比矩阵是整个程序最核心也最易错的部分。我建议把残差矢量组织成清晰的分段形式，电网残差放前面、气网残差居中、热网残差放后面，这样调试时看哪个分量不收敛一目了然。雅可比矩阵用稀疏矩阵表示，避免大矩阵运算内存爆炸。

一个重要技巧：让耦合设备的残差方程尽量写成显式关系。比如CHP，与其在耦合矩阵里写复杂的隐式方程，不如直接在残差函数里写成：

wasted_heat = 0; % 说明：此处仅作为示意 R_chp_e = unit.P_elec - unit.eta_e * unit.F_gas * unit.LHV; R_chp_h = unit.P_heat - unit.eta_h * unit.F_gas * unit.LHV;

这样残差语义清晰，雅可比矩阵的元素也能手算校核。

另一个我后来才掌握的技巧：先跑纯电网、纯气网、纯热网三个子程序，全部收敛了再联立。如果联立后不收敛，基本可以确定是耦合变量或交叉导数写错了。可以加一个开关把耦合项的雅可比设置为零矩阵，验证残差是否对应独立的三个网络。

5. 收敛性调试实录：初值、阻尼与气网的非凸难题

5.1 气网Weymouth方程的非凸性：迭代为什么会在两个解之间跳

Weymouth方程F² = C²(p_i² - p_j²)里的流量和压差是平方关系，它的解域是凸的吗？不是。同一个F可以对应两个符号相反的压差，方程在F=0处不光滑。牛顿法在高非线性区域很容易在多个根之间跳跃甚至震荡。我在调一个含压缩机的气网算例时，出现了压力在1.2 MPa和0.8 MPa之间来回跳的现象，雅可比残差总是降不到1e-6。解决方法是给流量变量一个初值符号约束（根据管道连接关系预判方向），同时将方程改写为：

R_ij = F_ij - C_ij * sqrt(p_i^2 - p_j^2) * sign(p_i - p_j)

这样能避免负值平方根带来的NaN问题，也能大幅改善迭代稳定性。

5.2 热网"冷启动"困难：温度初值怎么给才合理

热网计算里最容易犯的错误是一上来就零初值。给定热负荷为1.5 MW的节点，如果初始温度设成0℃（绝对温度），管道的散热模型和节点混合方程会产生巨大的负雅可比元素，导致迭代直接发散。我的经验是：用略高于环境的温度值作为初值，比如统一设供温75℃、回温45℃、环境温度15℃，然后跑两层固定点迭代把它稳定下来再进入牛顿主循环。另外，热网的流量和温度量级差异大（流量几百kg/h，温度几十度），不归一化的话雅可比矩阵的条件数会很差。我给每个变量指定了一个基准值，比如流量基准1000 kg/h、温度基准100℃，在雅可比组装后统一缩放，条件数立刻改善了两个数量级。

5.3 我常用的四个工程化收敛辅助手段

如果加阻尼之后还是收敛不理想，我一般按以下顺序排查：

第一个手段是残差归一化。电网的有功不平衡量级在MW、无功在MVar，气网的流量不平衡在kg/s，热网的功率不平衡在kW，直接比较它们的范数会掩盖小量级分量的问题。我通常把各类残差除以其所在网络的总负荷量级，统一变成无量纲相对误差，再统一看最大值。

第二个手段是迭代步长限制。每次修正量δX不得超过该状态变量初值的某个百分比，比如电压不超过0.1 pu、压力不超过0.1 MPa、温度不超过5℃。防止牛顿法大步长跳过可行域。这个可以用以下代码简单实现：

max_step = 0.1 * abs(X0 + eps); delta_X = sign(delta_X) .* min(abs(delta_X), max_step);

第三个手段是迭代顺序调整。分解迭代法的顺序直接影响收敛速度。我试过几种顺序后，认为先算气网、再算电网、最后算热网的顺序最稳，原因是气网为CHP和锅炉提供了燃料边界，确定了这个边界后电网热网的源项就相对确定了，再算后两者不像先算电网那么盲目。

第四个手段是中等规模系统下改用信赖域算法替代线搜索。Matlab的fsolve有信赖域选项，多能流这类高度非线性的问题上比自带的高斯牛顿法稳得多。我在统一求解的主循环里没有自己写信赖域，而是直接调用优化工具箱的fsolve配合自带的雅可比函数句柄，省了很多麻烦。

5.4 三类常见的"假收敛"以及排查思路

有些情况下迭代停下来了，但结果明显不对。我整理的排查思路是：

第一类是孤立网络平衡但耦合节点失衡。可能因为耦合设备残差没有加入整体残差向量，各网络独自平衡了，但CHP电出力和气耗量之间却不满足设备模型。排查办法是算完后把每个耦合设备的输入输出都打印出来，人工代入换算效率，看和设定值是否吻合。

第二类是压力/温度出现物理量纲不合理但迭代仍然继续。可能原因是初值方向错误，迭代走到了另一个数学解。比如热网远端回水温度算出来比供水还高，这种解满足方程组但从物理上看不可能。排查办法是对状态变量加边界约束，让残差函数在越界时返回一个大数惩罚，或者直接在一次牛顿步后判断上下界。

第三类是迭代收敛但残差极小但流量不平衡。因为在管网模型中，某些循环管路（环网）里的流量解不唯一，水力方程在某些构型下存在多解。解决办法是检查环路独立方程是否被从雅可比矩阵中错误地删除了，或者对环网添加基准流量约束。

6. 一个典型算例的结果解读：从迭代残差到工程判断

6.1 算例配置与初始条件

我用一个中等规模的测试系统演示总体效果：33节点配电网、20节点气网（含2个气源节点、1个压缩机站）、16节点热网（含1个热源站）。耦合设备配置为2台CHP（每台额定电出力1.5 MW，热电比0.9）、1台热泵（COP 3.5，额定热出力1.2 MW）、2台燃气锅炉（每台热出力2 MW）。总电负荷约8.5 MW，总气负荷约4.2 MW（等效），总热负荷约6.8 MW。

迭代初值：电网各节点电压1.0 pu、相角0；气网各节点压力统一设2.0 MPa；热网供水温度75℃、回水温度45℃，各管段流量按负荷近似预分配。迭代目标：最大残差小于1e-7（无量纲归一后）。

6.2 迭代过程与收敛特征

这个算例用统一求解法，第1次迭代残差降到3.7e-1，第3次降到2.1e-3，第5次落到8.6e-6，第6次达标。对比一下：同样算例用分解迭代法，需要24次迭代，且由于耦合设备功率占比偏高，松弛因子设到0.45才不震荡。

收敛后各网络的关键指标：

6.3 怎么判断这个结果算对了：三类平衡校验

多能流计算有一个优势是天然存在多组可交叉验证的平衡关系，只要有一个对不上，基本可以确定算错了：

电力平衡：平衡节点注入功率 + CHP发电 - 常规用电负荷 - 热泵电耗 - 电驱动压缩机功耗 = 0（网损约束）。

天然气平衡：气源供气量 = CHP耗气量 + 锅炉耗气量 + 压缩机能耗 + 气负荷，这个用燃气负荷总量和Weymouth方程算出的管段流量差值验证。

热力平衡：CHP热出力 + 锅炉热出力 + 热泵热出力 = 热网负荷总和 + 管道散热损失。这个特别容易忽略管道散热，我在一个算例里发现热源总供热比热负荷多出12%，一开始以为是计算错误，后来单独算了散热才发现是管道损失占的。

除了三个网络的平衡，还要检查耦合设备自身的模型一致性。CHP的发电效率η_e = P_elec / (F_gas * LHV)要落在合理区间（0.25~0.45），热泵的实际COP = Q_heat / P_elec要接近设定值（偏差不超过5%），如果偏差过大，说明耦合方程在联立解中被带偏了，优先检查雅可比矩阵中对应块的交叉导数是否正确。

6.4 联立计算后发现了哪些独立计算看不到的现象

这个算例最有意思的地方在于联立计算后才暴露出来的问题：当热网负荷从6.8 MW增加到8.0 MW时，热泵热出力需求增加，电耗从0.38 MW增加到0.43 MW，电网潮流重新分布；CHP热出力需求增加，气耗量上升，气网末端压力从1.33 MPa跌到1.20 MPa；气网压力下降后，燃气轮机的出力能力受限，平衡节点被迫多注入0.6 MW。这串连锁反应，如果还用三个网络独立计算，完全看不到——独立计算会假设热泵电耗固定、CHP气耗固定，结果必然低估了高负荷时电网和气网的压力。

这也是为什么我一直认为：多能流计算不是学术上的炫技，它解决的是工程实际中的"账目对不上"问题。做园区能源规划、做设备选型、做运行方式校核之前，先做一次全耦合的电气热能流计算，比任何经验估算都靠谱。

7. 给入门者的三个实用建议

最后分享几点实际做项目过程中的经验。

第一，不要一上来就追求大规模高精度。先把一个最简单的三节点电网、两节点气网、两个节点热网跑通，用线性方程组手解的结果验证程序正确，再逐步加节点、加设备。我一开始就奔着33节点去，结果调了一个月，后来用5节点小系统三天就通了。

第二，Matlab里做这个研究，工具箱够用就好。本质上核心代码就是牛顿-拉夫逊迭代加上稀疏矩阵运算，不需要一堆专业工具箱。我在项目里只用了基础的矩阵运算（\操作符）、稀疏矩阵函数（sparse、nnz），其他都是自己写的函数。优化工具箱的fsolve可以作为备选解法器。

第三，变量和量纲管理是整个程序的命门。建议用一套统一的单位体系贯穿始终，比如电网用标幺值、气网用MPa和kg/s、热网用kW和kg/h，在接口函数里统一换算。每写一个设备模型，就写一个单元测试函数验证输入输出是否正确，别等联立了再排查。

多能流计算这块内容，能聊的细节远远不止这些，但把上述这套模型和实现方法吃透，常见的电气热联算需求基本都能覆盖。后续如果想深入，可以做动态多能流、考虑不确定性场景、结合优化算法做日前调度，底层框架完全可以直接复用。