企业官网建设流程全解析

区域综合能源系统的能流计算，这两年是真的火。我年前接了个园区级综合能源规划项目，业主方开口就问：“你们能不能把电、气、热三个网放一个模型里算？”我第一反应是，这不就是潮流计算嘛，电网的潮流算过无数遍了，气网和热网无非就是多几个方程而已。真上手做的时候才发现，三套物理方程、三套单位制、三套初值习惯全部搅在一起，光是让雅可比矩阵不奇异就折腾了快两周。最后还是老老实实把“计及多能耦合的区域综合能源系统电气热能流联算”这套Matlab代码完整写出来，才算把这个事情真正跑通。

这篇文章就把我这次的实际做法完整拆开来讲：从多能耦合的数学模型怎么建，到电气热三个子网络在Matlab里如何统一求解，再到热网水力热力计算那些容易翻车的细节，最后是调试经验。内容适合正在做综合能源系统仿真的研究生、搞园区能源规划的工程师，以及想把“电热气联算”从论文搬进工程代码的人。代码思路可以直接抄作业，我会把为什么这样设计也一并说清楚。

1. 为什么要做电气热能流联算：单一系统计算的边界在哪

1.1 传统电力潮流的“盲区”

传统电力系统潮流计算，本质就是解一组节点功率平衡方程：给定负荷和发电机出力，求各节点电压幅值和相角。这个框架成熟、工具多，MATPOWER一个命令就能跑IEEE算例。但它有一个天然边界：只能回答“电怎么走”，回答不了“气往哪供、热往哪送”。

到了区域综合能源系统这个场景，能源站里燃气轮机发1度电要烧多少方天然气，余热锅炉能带多少平方米采暖负荷，电制氢设备吃多少电、产多少气——这些跨网络的关键关系，传统电潮流的数学模型里根本没有位置。你没法在只含电气量变量的方程里表达“燃气轮机进气量决定发电量、发电量的余热又决定供热能力”这种强耦合关系。

所以做区域综合能源系统的能流计算，第一步就是承认单网络计算已经不够用了，必须把电网、气网、热网放到同一个方程框架里，让它们通过耦合元件在数学上真正“对话”。

1.2 耦合元件：把三个网络绑在一起的关键

把三个网络连起来的是能源转换设备，行业内统一叫耦合元件。常见的四类：热电联产机组（CHP）、电转气设备（P2G）、电锅炉、热泵。它们的特点是，同一个物理设备同时出现在两个甚至三个子网络的方程里，变量在两个网络中都出现。

举个例子，CHP机组蒸汽轮机侧发电、抽汽侧供热，在电网方程里它是发电机节点（有功注入），在气网方程里它是天然气负荷节点（燃料消耗），在热网方程里它是热源节点（供热量）。一个设备，三个网络里的三类角色。传统的序贯解法是：先算电网，把CHP电出力给气网，气网算完把气耗返回来，再算热网，反复迭代直到稳定。这种做法实现简单，但边界条件来回传递时非常容易振荡，尤其是CHP的“以热定电”模式，热网温度一变，电出力跟着变，电潮流变了气耗又变，三层嵌套很有可能不收敛。

所以我在这次项目中直接采用统一求解法：把所有子网络方程和耦合元件方程全部堆成一个非线性方程组，用牛顿-拉夫逊法一次联立求解。这是这篇代码实现的核心路线，也是我认为最值得分享的部分。

2. 核心数学模型：电、气、热三套方程怎么统一起来

2.1 电网络侧：经典的极坐标牛顿法方程

电网部分，我沿用成熟的极坐标牛顿-拉夫逊法。对每个PQ节点，有功和无功注入方程分别是：

ΔPi = Pg_i - Pd_i - Vi * Σ(Vj * (Gij*cos(θi-θj) + Bij*sin(θi-θj))) = 0 ΔQi = Qg_i - Qd_i - Vi * Σ(Vj * (Gij*sin(θi-θj) - Bij*cos(θi-θj))) = 0

我没有在这个部分做任何创新，因为这个环节越标准越好。PV节点只保留有功方程，无功作为待求量；平衡节点电压幅值和相角都给定。这里要提醒一个点：在电热气联算时，电网方程的变量集合里除了电压幅值和相角，还会额外混入耦合设备的有功出力和天然气消耗量，所以电网子矩阵已经不是单纯的电网络雅可比了，需要从全局变量表统一索引。

2.2 气网络侧：Weymouth方程和节点流量平衡

气网建模的核心是管道流量与两端压力的关系。工程上最常用的是稳态Weymouth方程：

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

其中C_ij是管道常数，与管径、长度、温度、气体组分有关。这个sign函数在写Matlab代码时是个坑：求导时sqrt的导数是1/(2*sqrt(...))，如果压差为0导数直接无穷大，必须在代码里做防零保护。我采用的办法是在雅可比解析求导时，如果abs(p_i²-p_j²) < 1e-8，就用数值导数的eps近似，实测稳定很多。

气网的节点流量平衡方程是所有流向该节点的管道流量代数和，加上气源注入和负荷消耗，整体满足KCL形式：

ΣF_in - ΣF_out + Gs_i - Gd_i = 0

气网节点分为三类：压力已知的气源节点、负荷已知的普通节点、压力流量都不确定的中间节点。求解时至少需要一个压力参考节点，否则雅可比矩阵奇异——这跟电网里必须有平衡节点是一个道理。压缩机模型如果要做，一般把它当成一个消耗少量天然气的增压装置，在节点流量平衡里增加一个自耗气项，这次算例中我暂未包含压缩机，但代码框架里留好了接口位置。

2.3 热网络侧：水力模型和热力模型的双层结构

热网是所有子网络里最麻烦的，因为它不是一组方程能解决的，而是“水力”和“热力”两层耦合在一起。水力模型描述工质流量如何分配：

节点流量守恒：Σm_in = Σm_out 管道压降方程：p_loss = K * m_dot²（沿程阻力） 回路压降为零环方程：ΣΔp_loop = 0

热力模型的未知量是各节点的供水和回水温度，以及热源/负荷节点的热功率。节点温度的混合关系满足质量加权平均，管道温度损失简化为一阶指数模型：

T_out = T_env + (T_in - T_env) * exp(-λ*L / (Cp*m_dot))

这个方程里λ是管道综合传热系数，L是管长，Cp是比热容，m_dot是该管段质量流量。热负荷节点的热功率用下面这个公式计算：

Ф = Cp * m_dot * (T_supply - T_return)

联算时热网的未知量（节点供水温度、回水温度、管道流量）和电网电压、气网压力全部塞进同一个未知向量，热网的雅可比行列比电网和气网都大，这也是为什么很多现成代码里热网部分做得最粗糙的原因。

2.4 耦合元件建模与变量扩展

耦合元件方程是全局方程组的“粘合剂”。我这套代码里内置了三类，用参数配置切换：

CHP以热定电模式：给定热出力Ф_chp，则电出力P_chp = Ф_chp * η_e / η_t，天然气耗量F_gas = Ф_chp / (η_t * LHV_gas)，其中η_e、η_t分别是发电效率和热回收效率，LHV_gas是天然气低位热值。P2G设备：输入电功率P_p2g，产气量F_p2g = P_p2g * η_p2g / LHV_gas，形成一个电网负荷节点和一个气网注入源。电锅炉和热泵：输入电功率P，输出热功率Ф = COP * P，形成电网负荷节点和热网热源节点。

统一求解时，未知变量向量x就是一个混编数组，比如：

x = [ 电网V角; 电网V幅值; 气网节点压力; 热网供水温度; 热网回水温度; 热网管道流量; 耦合环节附加变量 ]

所有子网络和耦合元件的等式残差按顺序堆成一个大列向量F(x)，再对所有变量求偏导形成全局雅可比矩阵。下面是变量和方程对应关系的一个总览：

3. Matlab整体实现框架：从算例数据到统一求解器

3.1 输入数据结构：别用零散变量，用嵌套结构体组织系统

写Matlab代码最忌讳的就是堆一大把命名混乱的全局变量。我这个项目里所有系统数据都用嵌套结构体组织，每一类对象一个数组，字段名统一。

% 电网络节点：1平衡节点，2 PQ节点，3 PV节点 bus(1).type = 1; bus(1).Vm = 1.06; bus(1).Va = 0; bus(1).Pd = 0; bus(1).Qd = 0; bus(1).Pg = 0; bus(1).Qg = 0; bus(2).type = 2; bus(2).Vm = 1.0; bus(2).Va = 0; bus(2).Pd = 1.2; bus(2).Qd = 0.4; bus(2).Pg = 0; bus(2).Qg = 0; % 气网络节点：1气源节点，2负荷节点 gas(1).type = 1; gas(1).p = 4.0e6; gas(1).Gd = 0; % 节点压力给定值，单位Pa gas(2).type = 2; gas(2).p0 = 3.5e6; gas(2).Gd = 0.8; % 负荷m³/h % 热力站：供回水温度设定 heat_node(1).Ts = 90; heat_node(1).Tr = 50; % 一次网供回水温度℃ heat_load(1).Phi = 0.5; % 热负荷MW

嵌套结构体的好处是，后续写残差函数时可以非常方便地用for循环遍历所有节点组装方程，不需要反复根据名字查变量。这里的热网管道数据结构我单独用了一个管道数组，每条管道定义了起点、终点、长度L、内径D、传热系数λ。

3.2 统一非线性求解器：残差堆叠与雅可比组装

有了数据之后，核心就是写残差函数和雅可比矩阵。这部分我强烈建议不要一上来就推导每个元素的解析导数，先用数值雅可比把整体跑通，验证方程组装正确，再逐步替换成解析导数提升性能。数值雅可比的实现非常简单：

F0 = eval_residual(x, sys_data); % 初值残差 J = zeros(length(F0), length(x)); for k = 1:length(x) xp = x; xp(k) = xp(k) + 1e-7; J(:,k) = (eval_residual(xp, sys_data) - F0) / 1e-7; end

实测下来，对于100个变量以内的算例，数值雅可比每次迭代的耗时在几十毫秒级别，完全可接受。等到后面要做几百节点的大算例，再换解析雅可比不迟。这里还有个细节：统一的牛顿法迭代格式是dx = -J \ F，然后x = x + dx，但Matlab反斜杠对稀疏矩阵更友好，组装J时尽量用sparse函数，能处理更大规模的系统。

3.3 主迭代循环代码骨架

我直接给出这次项目里主求解器的骨架，照着写就能跑通一个基础版本：

function [x, hist] = ries_energy_flow(x0, sys, opt) % 区域综合能源系统统一能流求解器 % 输入：x0初值，sys系统结构体，opt选项 % 输出：x解向量，hist收敛历史 x = x0(:); tol = opt.tol; maxiter = opt.maxiter; hist = zeros(maxiter, 1); for iter = 1:maxiter F = eval_residual(x, sys); J = assemble_jacobian(x, sys); % 稀疏雅可比 dx = -J \ F; x = x + opt.damping * dx; normF = norm(F, inf); hist(iter) = normF; if normF < tol fprintf('收敛于第%d次迭代，残差=%.2e\n', iter, normF); return; end if norm(dx, inf) < opt.xtol fprintf('变量增量已很小，提前终止，残差=%.2e\n', normF); return; end end error('迭代超限，未收敛'); end

关键点在eval_residual这个函数里。我把它按模块拆成四段：电网残差、气网残差、热网残差、耦合设备残差，每段独立计算再拼接成大向量。这样调试时哪一段出问题可以直接定位。对于初值的阻尼系数，我默认设0.9，收敛慢一点但稳定，等跑通了再改回1.0提速。

3.4 关于Matlab版本和工具箱的碎碎念

好多人私信问我用的哪个版本，是不是必须用Simulink。这里说清楚：我用的Matlab 2023b，但这个代码只用到了基础矩阵运算和非线性方程求解，不依赖任何工具箱，从2016b到2026b都能跑。网上热传的那些Matlab安装过程中报类似License Manager Error -8的问题，基本都是激活文件hostid不匹配引起的，跟计算代码无关，别把两个问题混在一起排查。如果你用的是2026b预览版，跑这段代码也没有任何障碍，因为核心语法没变过。另外对于不想配环境的人，直接用Matlab Online网页版把.m文件传上去也能运行，我这个算例规模完全够用。

4. 热网水力-热力计算的实战细节

4.1 水力计算：回路压降方程绝对不能省

热网水力计算是初学最容易偷懒的地方。很多人直接用节点流量守恒把所有管道流量求出来，发现方程不够，就随便把某个流量设成已知。这种处理在简单辐射状热网里勉强能work，但一旦有环网，流量分配会完全不对。

真实的热网水力计算必须包括回路压降方程：对每个独立回路，沿回路一圈的压降代数和为零，否则就该有平面的水流循环出现，物理上不成立。加上这个方程后，管道流量才是完整可解的。这里我分享一个调参经验：管道阻力系数K和直径、粗糙度有关，如果不知道怎么给，先用经验公式K = 8λ_fL/(π² * ρ * D⁵)估算。λ_f是摩擦系数，水的工况下取0.02左右比较安全，实测下来流量分配不会差太远。

4.2 热力计算：温度是跟着工质“跑”的

热力计算里最重要的思维转变是，温度不是“节点属性”，而是跟随质量流量在管道里传播的。供水温度从热源出发，经过各段管道衰减，到达换热站后把热量交给二级网，自己变成回水温度；回水再沿着回水管网流回热源重新加热。

这个过程中有两点必须注意。第一，节点混合温度按质量流量加权平均，不能简单算术平均。比如两条管道分别以2kg/s和5kg/s流入同一个节点，温度分别是85°C和70°C，混合后是(285 + 570)/(2+5)，约74.3°C，而不是77.5°C。第二，管道温降公式里的m_dot在分母上，流量越小温度损耗越大，所以热网在低负荷工况下温度和流量需要联立求解，不能把流量算完再套温度。

我在代码里用一个纯函数计算管道出口温度：

function Tout = pipe_outlet_temp(Tin, m_dot, L, D, lambda, Tenv, Cp) % 管道出口温度，单位°C % 传热系数K_heat = lambda*L/D近似 exp_term = exp(-lambda * L / (Cp * m_dot * 1000)); % 换算kW Tout = Tenv + (Tin - Tenv) * exp_term; end

这里的Cp单位是kJ/(kg·°C)，乘以1000换算成J，避免单位错乱。这个函数我实测了很多次，只要lambda给得不是太离谱，温度和实际运行数据能对上。

4.3 单位的坑：MW、kPa、kg/s、°C混用是最大的错误来源

我必须专门用一节来讲单位，因为这是联算中最容易出错、却最容易被忽视的地方。电网里习惯用标幺值，电压是p.u.，功率是MW；气网里压力是kPa或MPa，流量是m³/h或kg/s；热网里流量是kg/s，热功率是MW，温度是°C。三个系统的单位习惯完全不一样，直接塞进同一个方程组，数值尺度差异会大到让雅可比矩阵病态。

我建议的做法是统一采用国际单位制作为“计算基准单位”：功率MW、压力Pa、质量流量kg/s、温度°C换算成K但温差保持用数值。温度使用绝对温度K作为内部变量，因为热力公式涉及环境温度相减时，用°C做单位没有问题，但在计算比热容相关的能量平衡时，统一K更安全。整个算例中我只允许在输入输出层做单位转换，计算层全部使用同一套单位系统，这个习惯帮我避免了至少三次“看起来收敛了但结果物理上不对”的尴尬。

4.4 热网迭代策略：内层可以先自解，外层再联立

虽然整体用统一求解法，但热网内部我还是做了特殊处理：在每一轮全局牛顿迭代中，水力方程和热力方程各自有自己的收敛逻辑。全局迭代负责把电、气、热拉通，热网内部通过一个内层循环先完成水力-热力自洽，再返回给全局残差函数。这样做的原因是热网的方程具有很强的“局部性”，如果完全平铺到全局方程里，变量之间的相互依赖链条太长，数值收敛性反而变差。

具体做法是：在eval_residual函数里，当计算热网残差时，先调用一个solve_hydraulic_thermal_local函数，输入各耦合节点的边界约束，输出热网内部一致的温度和流量分布，再把它代入热网残差方程。这种“全局联立、局部自洽”的混合策略，从收敛性和计算效率两个角度看，都比纯平铺好不少，也算是我这次实践的一个原创优化。

5. 常见问题与排查技巧实录

5.1 初值设置：纯零初值在气网里必炸

这是所有刚接触气网计算的人都要踩的坑。给气网节点压力初值设为0，代入Weymouth方程后sqrt(|0-0|)这一项以及它的导数都是0，整个气网子矩阵全零，雅可比矩阵直接奇异。正确做法是：所有气网节点压力初值给在气压站出口压力附近，比如4MPa的管网，初值给3.5~4.0MPa之间，让sqrt内部的压差至少有0.5MPa的正量，雅可比矩阵才不会出现0列。

热网的温度初值也一样，别给环境温度（比如20°C），否则管道温降公式里(Tin - Tenv)是0，雅可比矩阵对应位置全零。我给的是一个“线性斜坡”初值：热源出口温度给的规整值90°C，沿线节点按照离热源距离递减，回水管网从50°C反向递增回来。这样初始雅可比矩阵大致反映了真实物理方向，牛顿法的收敛球半径能大不少。

5.2 变量尺度差距导致的病态雅可比矩阵

统一求解的第一个版本，我直接把实际单位的变量堆在一起求解，结果迭代全程在“振荡-发散-振荡-发散”中循环，MATLAB报矩阵接近奇异。排查后发现，电网电压值是1.0左右（标幺值），气网压力是4e6级别，热网流量是几十，变量尺度差了6个数量级，雅可比矩阵条件数直接破10的8次方。

解决办法有两个，推荐第二个。第一个是每个变量除以该变量的基准值实现无量纲化（类似电力系统的标幺值思想）；第二个是在迭代更新时做对角缩放：dx = - (D \ J) \ (D \ F)，其中D是对角元素绝对值开方构成的对角阵。用缩放后的矩阵求逆，再转换回原变量空间，实测条件数能下降几个量级。这一段我建议直接把变量基准值写在代码开头的注释里，方便其他接手人理解。

5.3 收敛判据只看残差范数还不够

我第一版收斂判据只检查norm(F, inf) < 1e-6，结果遇到过两次“残差很小但变量还在缓慢漂移”的情况。原因是个别变量的量纲很大（如压力），它在残差表达式里天然占的权重小，残差范数达标时它其实还差很远。

正确的做法是双重判据：残差向量和变量增量向量同时检查。代码骨架里我已经写了这两条。此外还要注意，牛顿法在联算中如果阻尼系数取1.0，偶尔会震荡，我会在检测到本轮残差比上轮更大时自动砍半阻尼，代码里加一个简单的回溯线搜索就是三行的事，但对收敛性帮助极大。

5.4 求解失败排查顺序速查表

我把这次调试中遇到最多的几类“不收敛”问题整理成一张表，每一种都直接对应解决方案：

6. 算例验证与后续扩展

6.1 怎么验证自己的计算结果是可信的

写完代码，千万不能看它收敛了就认为结果对。我惯用的验证方法是从简到繁分三层。第一层，把气负荷和热负荷全设为0，系统退化为纯电网潮流，用MATPOWER跑同一个电网结构，两边结果对比，电压幅值和相角的误差要小于1e-6。第二层，把电负荷和热负荷设为0，退化为纯气网稳态计算，用气体管网公式手算一个简单节点对，核对压力值。第三层，设置一个单热源、单负荷的小型热网，手算管道温降和节点混合温度，跟代码输出对比。

这三层能过，再检查全局能量守恒：系统输入的总能量（气源消耗+电网注入+热源燃料）必须等于总负荷加上各个网络损耗。这一步用一句话总结就是“能量既不会凭空产生，也不会凭空消失”，算完以后所有項加起来如果超出误差范围，肯定哪里有问题。我把这个能量校验直接写成脚本，每调一次代码就跑一遍，发现问题立即定位。

6.2 这个基础版本还能往哪些方向扩展

跑通基础能流计算之后，往上扩展的方向很多。按我目前的经验，最适合在这个框架上做的是下面三类。

第一，把确定性能流升级为不确定性能流。光伏出力、风电出力和负荷随机性，可以通过蒙特卡洛抽样反复调用这个能流求解器，然后统计电压、压力、温度的分布区间。我实测一个100节点规模的小系统，抽500次样本，Matlab单线程跑大约需要20秒，完全可接受。第二，在能流计算基础上加优化层，做区域综合能源系统的最优能流（OEF）求解，把CHP的运行区间、气网压力上下限、热网温差不越限作为约束，目标函数用运行成本和碳排放的线性组合。第三，把稳态模型改造成准动态模型，以15分钟为步长做全天的时序仿真，能评估储热罐、储气罐的调节作用。

这些扩展都只需要在现有残差函数和雅可比组装基础上做加法，不需要推翻重构，这也是我当初坚持用统一求解法而不是序贯法的原因之一。

最后再分享一个调试的小技巧。如果你在联算时反复找不到问题，试着把耦合系数强制设为0，让系统退回“电网独立运行，气网独立运行，热网独立运行”，分别在各自子网络求解器下跑通，然后再一次性恢复耦合系数。这个“先切开后拼上”的方式，能帮你快速分清问题出在某个子网络的方程里，还是出在耦合环节的建模上。我在这次项目里至少用这个方法定位了三处BUG，希望对你也有用。