企业官网建设流程全解析

做区域综合能源系统仿真的朋友，几乎都绕不开“电气热能流计算”这道坎。我接这个课题时，正处于几种状态叠加：Matlab装了不下三次、文献翻了一摞、模型反复推倒重来、代码天天在跑又天天发散。痛定思痛后，我决定把“计及多能耦合的区域综合能源系统电气热能流计算”完整梳理一遍，并用Matlab从零搭出一套能跑通、易扩展、可后续接优化算法的计算程序。这篇文章就围绕这套代码的建模思路、求解方案和实战经验展开，计划做冷热电联供、多能互补、综合能源系统规划方向的同学，以及想快速实现能流计算的电气/热动工程师，都能从这里找到可直接参考的框架。

1. 项目概述与研究思路

1.1 单能源网络能流计算的明显短板

传统电力系统潮流计算只考虑发电机、负荷和线路，天然气管网分析只关心气压与流量，供热管网设计只盯着温度和热负荷。这三套模型独立运行在很多场景下够用，但只要系统里出现跨网络转换设备——比如燃气轮机、燃气锅炉、电锅炉、电转气P2G、余热回收——单一能流计算就失灵了。原因很简单：CHP机组从天然气网取气、向电网送电、向热网供热，电气热三个网络因此被同一台设备强制绑定。你单独算电网时，气网压力变化了，CHP的天然气量就不一样，电出力跟着变，电网潮流必然重新分布。

这就像三个原本各自独立的水池被水管连通了，只盯着一个水池算水位，另外两个水池一旦波动，这个水池的水位根本稳不住。所以必须建立多能耦合的综合能源系统模型，把电网、气网、热网放到同一个计算框架里，用统一的能流计算来反映设备工况变化带来的跨网影响。这个统一能流计算不是把三个程序简单拼在一起，而是要仔细处理好网络边界、设备互连方程和迭代收敛机制。

1.2 整体技术路线与求解方案选型

实现电气热能流计算，业界主流通常有两条路线：统一求解法和顺序迭代法。统一求解法把所有网络方程和耦合约束组合成一个大方程组，同时迭代求解。优点是理论上可以统一处理强耦合、收敛性能更稳定，缺点是雅可比矩阵维度大、初值极其敏感、程序模块化程度低，一旦某个子网络参数变化，整个矩阵的索引和分块逻辑都要跟着改。

顺序迭代法则是把三个网络分别求解，只在耦合设备节点处交换功率、流量等变量，反复迭代到稳定。它牺牲了一部分理论上的强耦合性，但换来的是工程上的灵活性。我最终选了顺序迭代法，原因有三个：一是项目周期紧，顺序迭代可以直接复用成熟的单网能流代码，比如电网潮流就采用电力系统分析中经典的牛顿—拉夫逊法，气网和热网各自写一套迭代求解器，每个子网络单独调试，问题容易定位；二是后续计划做多能流优化调度，顺序迭代的模块化结构更容易嵌入目标函数和约束条件；三是在区域综合能源系统这类中等规模场景下，顺序迭代法的计算时间通常在毫秒到秒级，完全满足科研验证需求。

计算框架的整体思路可以概括为：先把电、气、热各自的网络参数和负荷数据准备好，然后初始化耦合设备变量；随后按照“电网→气网→热网”的顺序循环求解子网络能流；每轮根据耦合设备模型更新耦合变量，包括CHP的发电出力、耗气量、供热功率，以及P2G的耗电量和产气量；最后判断所有耦合变量的变化量是否小于收敛阈值。如果发散，就针对具体原因调整初值、收敛判据或迭代顺序。这个框架在后面代码章节会拆解到函数级，这里先建立整体认知。

2. 核心模型与数学基础

2.1 电网潮流模型：牛顿—拉夫逊法怎么落地

电网子系统的能流计算采用极坐标形式的牛顿—拉夫逊法。节点有功、无功功率平衡方程是核心，任意节点i的有功注入和无功注入分别满足：

P_i = Σ_j V_i V_j (G_ij cosθ_ij + B_ij sinθ_ij)

Q_i = Σ_j V_i V_j (G_ij sinθ_ij - B_ij cosθ_ij)

其中θ_ij = θ_i - θ_j，V_i是节点i的电压幅值，G_ij和B_ij是节点导纳矩阵的实部和虚部。程序里需要形成不平衡量向量，包含所有PQ节点的有功、无功不平衡量以及PV节点的有功不平衡量，然后求解雅可比矩阵方程，获得相角和电压幅值的修正量。迭代初值一般选平启动，即所有PQ节点电压幅值1.0标幺、相角0度，平衡节点电压和相角固定。

这个模型本身不复杂，但工程实现时有两个讲究。第一是标幺化处理，负荷、发电机和线路阻抗都要统一到同一个功率基准，通常取100MVA或10MVA。标幺化不是为了让数字看起来顺眼，而是让电压、功率、阻抗的数量级都落在1附近，雅可比矩阵的病态程度大幅降低，Matlab里求解线性方程也更稳。第二是PV节点的无功越限处理，如果某台发电机在计算中无功出力超过上限，就需要把它从PV节点改成PQ节点，或者直接去掉该节点的无功平衡方程再重算，否则迭代过程容易来回振荡。我第一次完整算例时就遇到这个坑，后面排查了很久才发现是光伏节点的无功越限没处理。

2.2 天然气网能流模型：Weymouth方程和节点气压

天然气管网计算的关键是管段流量和节点气压的耦合关系。稳态工况下，管段从节点m流向节点n的流量常用Weymouth方程描述：

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

这里π是节点气压，单位常用MPa，C_mn是管道常数，由管径、摩擦系数、气体温度、压缩因子等决定，sign函数用来表达流量方向。如果π_m大于π_n，流量从m流向n，sign取正；反之取负。直接对气压而不是气压平方求导，方程会产生根号项和方向项，数值上不太稳定，所以我在程序中把求解变量设置成气压平方π²，让方程变成线性平方差形式，牛顿法的雅可比矩阵更规整。

对每个不含气源的节点，需要建立流量平衡方程：流入该节点的流量总和加上气源注入或P2G注入，等于该节点的气负荷。如果管网里存在压缩机，要在压缩机节点单独加一个增压比关系方程，同时把压缩机的耗气量加到对应节点负荷里。压缩机相当于给天然气网络增加了一个“抬高压力”的元件，它自己也要消耗一部分天然气才能工作，这部分消耗往往被初学者忽略。气网求解的初值我推荐取设计压力的0.8倍，而不是1.0倍，因为在Weymouth方程的二次型结构下，初值太靠上会让雅可比矩阵初始斜率过大，迭代步长容易越过真解。

2.3 热力网络能流模型：水力与热力分开算再耦合

热网模型比电网和气网复杂一些，因为它涉及水力工况和热力工况两个层面。水力计算负责求取管段质量流量和节点压力或压头，热力计算负责求取节点供水温度和回水温度。整体采用“水力—热力”解耦迭代：先根据热负荷和假定的供回水温度求质量流量，再固定流量求温度分布，然后根据温度分布修正流量，反复到收敛。

简化起见，我采用节点法。对每个负荷节点，热负荷功率满足 H_i = c_p · m_i · (T_s,i - T_r,i)，其中c_p是水的比热，m_i是流经该负荷的质量流量，T_s和T_r分别是供水温度和回水温度。管道温度损耗采用一阶散热模型：

T_end = (T_start - T_a) × exp(-λL / (c_p m)) + T_a

其中T_a是环境温度，λ是管道单位长度传热系数，L是管长，m是质量流量。算的时候务必注意质量流量单位是kg/s，热量单位是kW，否则常会出现差三个数量级的情况。热网通常是辐射状或弱环网，我的做法是先判断拓扑是否为树状，如果是树状就从热源逐管段推流量，如果是环网就要用图论中的基本回路法列压降平衡方程，程序复杂度会升高一个级别。区域综合能源系统里绝大多数热网按辐射状设计，所以逐支路法够用。

2.4 耦合设备模型：CHP、P2G和燃气锅炉的数学表达

耦合设备是电气热三个网络的“粘合剂”。我的代码里主要放了四类：抽凝式CHP、燃气锅炉、电锅炉和P2G。这里给最常用的简化模型。抽凝式CHP在给定的进气量下，电出力P_e和热出力H_th之间有一个可行域，通常表现为多边形约束。在能流计算中，我先取一个固定热电比β，使H_th = β·P_e，再把CHP的天然气消耗量F_gas按综合效率反算。实际程序里，将F_gas作为气网的负荷增量、P_e作为电网的电源注入、H_th作为热网的热源注入，三个网络就通过一个CHP变量实现了互联。

P2G的模型更直接：电转气过程消耗电功率P_e，产生天然气流量F_gas = η_p2g · P_e / LHV_natgas。在电网里表现为负荷，在气网里表现为气源注入。燃气锅炉则消耗天然气、输出热功率H_gb = η_gb · F_gas · LHV_natgas，它把气网和热网耦合起来。电锅炉则把电网和热网耦合起来。这些耦合设备的模型都必须写成可调用的Matlab函数，外层迭代时反复调用即可。每个设备函数内部至少要包含两条信息：一是正常工况下的变量转换关系，二是出力上下限和可行域判断。能流计算如果算出来的CHP出力超出可行域，外层循环就必须做限制处理，否则算出的结果再收敛也是没有物理意义的。

3. Matlab代码实现架构与关键函数

3.1 数据结构设计与文本参数读取

代码要易读，数据组织别乱。我用struct数组存三类网络的节点和支路信息。电网友节点存id、类型（平衡/PV/PQ）、有功负荷Pd、无功负荷Qd、有功发电Pg、无功发电Qg、电压幅值V、相角theta；支路存from、to、电阻R、电抗X、对地电纳B、变比tap。气网节点存id、类型（气源/负荷/连接点）、压力p、注入量inj、负荷load；支路存from、to、管径、长度、管道常数C、压缩机标志。热网节点存id、类型（热源/负荷/连接点）、供水温度T_supply、回水温度T_return、热负荷heat、质量流量mass_flow；管道存from、to、长度length、传热系数lambda、管径。

所有网络参数和负荷数据都存放在文本文档中，主程序用readtable统一读取。这个习惯帮我省了大量调整参数的时间，因为算例参数经常要改，如果硬编码在脚本里，每改一次就要去翻代码。我推荐数据文件按列命名清晰，比如“case_electric_bus.txt”的列名就是“id,type,Pd,Qd,Pg,Qg,V,theta”，读进来以后直接转成数组，后续都传struct，别在函数里到处用全局变量。下面这段是数据读取的基本写法：

%% 读取电网、气网、热网数据 busE = readtable('case_electric_bus.txt'); branchE = readtable('case_electric_branch.txt'); busG = readtable('case_gas_node.txt'); branchG = readtable('case_gas_pipe.txt'); busH = readtable('case_heat_node.txt'); branchH = readtable('case_heat_pipe.txt'); % 耦合设备参数 dev = readtable('coupling_device.txt');

节点编号最好从1开始连续编号，这样索引映射最简单。如果遇到非连续编号，就要额外加一张映射表，既浪费内存，还容易出索引错位的问题。

3.2 电网能流核心函数的实现要点

核心函数名我用“elecPowerFlow.m”，输入nodeE、branchE和电压初值，输出节点电压和功率分布。函数主体就是牛顿—拉夫逊迭代，下面给出一段经过注释的关键代码：

function [nodeE, iter] = elecPowerFlow(nodeE, branchE, tol) % 牛顿-拉夫逊法求解电网潮流 % 节点类型：1平衡，2PV，3PQ，这里默认平衡节点为1号节点 maxIter = 30; V = nodeE.V; theta = nodeE.theta; for iter = 1:maxIter [Pcal, Qcal] = calcInject(nodeE, branchE, V, theta); dP = nodeE.Pg - nodeE.Pd - Pcal; dQ = nodeE.Qg - nodeE.Qd - Qcal; % 组装不平衡量，去掉平衡节点行，PV节点去掉无功方程 dPQ = assembleMismatch(dP, dQ, nodeE.type); if max(abs(dPQ)) < tol break; end J = buildJacobian(nodeE, branchE, V, theta); dx = J \ dPQ; % 更新角度与电压，平衡节点不更新 idx = find(nodeE.type ~= 1); theta(idx) = theta(idx) + dx(1:length(idx)); V(idx) = V(idx) + dx(length(idx)+1:end); end nodeE.V = V; nodeE.theta = theta; end

这里需要强调两个细节。第一，雅可比矩阵的分块索引必须与节点编号保持一致，建议用一个单独的索引映射函数来管理，不要直接硬编码行列号。第二，如果某个PV节点在迭代过程中无功越限，要在每次迭代前检查一次，越限了就把它降级为PQ节点，并把无功定值设为限值。这段代码里的“calcInject”和“buildJacobian”是子函数，建议单独写成m文件，方便后续做其它算例时复用。

3.3 气网能流核心函数的实现要点

气网求解函数“gasPowerFlow.m”的核心是用牛顿法迭代节点气压平方。首先根据管道参数和初始气压计算所有管段流量，然后对非气源节点列出流量残差，用雅可比更新气压平方。核心代码片段如下：

function [nodeG, iter] = gasPowerFlow(nodeG, branchG, tol) % 天然气网牛顿法能流计算，求解变量取节点气压平方p2 N = height(nodeG); p2 = nodeG.p.^2; maxIter = 30; for iter = 1:maxIter f = calcPipeFlow(branchG, nodeG, p2); % 计算所有管段流量 F = zeros(N, 1); for n = 1:N if nodeG.type(n) == 1 % 气源节点，压力给定，跳过 continue; end F(n) = sum(f(branchG.to == n)) - sum(f(branchG.from == n)) ... - nodeG.load(n) + nodeG.inj(n); end if max(abs(F)) < tol break; end J = buildGasJacobian(branchG, nodeG, p2); % 仅更新非气源节点 idx = find(nodeG.type ~= 1); p2(idx) = p2(idx) - J \ F(idx); end nodeG.p = sqrt(p2); end

气网程序最容易出问题的就是方向符号。Weymouth方程里如果平方差是负数，直接开方会出错，所以在“calcPipeFlow”里必须先计算dp2 = π_m² - π_n²，再用sign(dp2)*sqrt(abs(dp2))处理，这样管段反向流动也能正常工作。压缩机节点我单独处理，把它看成一个提升比，出口压力固定在给定值，压缩机消耗的燃料按比例加到所在节点的负荷中。这个处理方式虽然忽略了一些动态特性，但对于稳态能流计算已经足够精确。

3.4 热力网络能流函数的实现要点

热力计算我写了两个函数：“heatHydraulic.m”负责水力计算，“heatThermal.m”负责温度计算。热网通常是辐射状或弱环网，水力计算从热源节点开始逐段推算出各管段质量流量，热力计算从源到负荷沿管段传递温度。关键的温度计算代码片段如下：

function [nodeH] = heatThermal(nodeH, branchH, env) % 节点供水温度计算，从热源出发沿流动方向求解 % env.Ta为环境温度，env.cp为比热 for k = 1:height(branchH) fromN = branchH.from(k); toN = branchH.to(k); T_start = nodeH.T_supply(fromN); lam = branchH.lambda(k); L = branchH.length(k); mdot = branchH.mass_flow(k); if mdot < 1e-6 continue; % 防止零流量导致指数异常 end T_end = (T_start - env.Ta) * exp(-lam * L / (env.cp * mdot)) + env.Ta; nodeH.T_supply(toN) = T_end; end % 根据热负荷和流量反算回水温度 nodeH.T_return = nodeH.T_supply - nodeH.heat ./ (env.cp * nodeH.mass_flow); end

注意：供水网络和回水网络通常拓扑对称但方向相反。回水温度要通过热负荷平衡方程逐点反算，如果某个负荷节点有多条来水管道，就应该按流量加权混合温度，不能简单平均。还有一个小细节：管段质量流量不能出现接近0的值，否则散热方程里的指数项会溢出，算出的温度会变成负数。所以水力计算完成后要先检查流量，低于阈值的管段要单独处理。

3.5 多能耦合外层迭代与收敛控制

有了三个子网求解器，剩下就是把它们“缝合”到一起。我的外循环主函数结构如下：

function [res, iterOut] = IES_PowerFlow(casefile) % 初始化耦合设备 dev = initCouplingDevice(casefile); for kOut = 1:20 % 1) 根据当前CHP/P2G状态，更新电网边界 nodeE.Pg = basePg + dev.Pe_chp; nodeE.Pd = basePd + dev.Pe_p2g; % P2G耗电 % 2) 更新气网边界 nodeG.load = baseGLoad + dev.Fg_chp + dev.Fg_gb; nodeG.inj = baseGInj + dev.Fg_p2g; % 3) 更新热网边界 nodeH.heatSource = baseHeat + dev.Hth_chp + dev.Hth_gb + dev.Hth_eb; % 4) 分别求解子网络 [nodeE, ~] = elecPowerFlow(nodeE, branchE, 1e-6); [nodeG, ~] = gasPowerFlow(nodeG, branchG, 1e-6); [nodeH, ~] = heatPowerFlow(nodeH, branchH, 1e-6); % 5) 根据子网新状态更新耦合设备参数 devNew = updateCouplingDevice(dev, nodeE, nodeG, nodeH); hist(kOut) = max(abs([devNew.Pe_chp - dev.Pe_chp; devNew.Fg_chp - dev.Fg_chp; devNew.Hth_chp - dev.Hth_chp])); if hist(kOut) < 1e-5 break; end dev = dev + 0.5 * (devNew - dev); end end

注意外循环更新策略很关键。我踩过的坑是直接把新设备出力整体替换旧值，遇到强耦合场景经常发散。后来给更新加了一个阻尼系数alpha=0.5，相当于给迭代加低通滤波，抵掉高频振荡，收敛性大幅改善。每个子网求解器的输出反过来会成为另一个子网的边界，因此子网内部也要做保护处理，比如气网节点气压越界时及时返回错误码，不要让外层循环背着错误继续跑。

4. 算例验证与结果分析

4.1 测试系统参数与耦合设备接入

为验证代码正确性，我搭了一个中等规模算例：电网采用修改的13节点辐射配电网，气网为6节点环状结构，热网为6节点辐射状结构，三个网络通过2台CHP、1台燃气锅炉、1台P2G耦合。网络参数和负荷数据都放在文本文件中。设备参数整理如下表：

约束条件为：电网平衡节点电压标幺值1.0，气网气源压力1.0MPa，热网供水温度90℃。因为这个算例主要用来验证算法框架，所以负荷我都取了比较常规的冬季度典型值，没有刻意拉太高。

4.2 电/气/热流计算结果与收敛性分析

代码跑通后，我记录了一组典型结果：电网所有节点电压幅值在0.95~1.05pu之间，最低电压出现在P2G接入的10号节点附近，这说明电转气负荷对馈线末端电压有明显影响；气网节点压力在0.82~1.00MPa之间，CHP2所在节点耗气量大，压力相对较低；热网供水温度从热源到末端下降约6℃，回水温度基本维持在60℃左右，符合设计预期。外层迭代到第9次收敛，最大耦合变量偏差小于1e-5，总耗时约0.8秒。

作为对比，把相同算例中的P2G停运后再算，外层迭代第5次就收敛了，可见强耦合会让收敛速度明显变慢。这个结果对做系统规划很有参考意义：P2G和CHP配置过多，影响的不只是能量平衡，连稳态能流都会变得“更硬”，初值稍微差一点就发散。同时我也验证了阻尼系数的影响：把外循环阻尼从0.5改为1.0，也就是无阻尼直接替换，算例的迭代次数从第9次变成第15次，而且中间有两次耦合变量大幅振荡。这说明顺序迭代里阻尼不是可有可无的调参项，而是保证数值稳定的必要手段。

5. 工程实现中的坑与排查技巧

5.1 初值选择影响：从冷启动到热启动

我第一次跑气网时，用的是全节点压力标幺值1.0，结果牛顿法死活不收敛。后来发现Weymouth方程是二次型，初值取1.0时气压平方差很大，雅可比矩阵容易把迭代带飞。改成取设计压力附近的值后，问题立刻消失。电网初值用平启动基本没问题，但热网的泵压和质量流量不能乱设，流量初值如果给得太小，管道散热方程里的指数项会变得特别大，温度直接算成负数。我建议先用一次粗略的水力计算估计质量流量，再进入完整“水力—热力”迭代，不要一上来就直接联合算。

5.2 单位统一与标幺化处理：最容易错的地方

多能流计算最容易出问题的不是算法，而是单位。电网习惯用标幺值，气网习惯用MPa和立方米每小时，热网习惯用kW和kg/s，三者混在一起，稍不注意就会出现“万”和“千”差三个数量级的迷惑。我的做法是内部计算全部统一到国际单位：功率用kW、气压用MPa、流量用kg/s，最后显示时才转换。特别提醒天然气的热值LHV常见单位是kJ/m³，换算成kW时一定要乘流量再除以3600。如果程序里算出的设备耗气量明显偏离物理直觉，第一件事就是去查单位换算表。

5.3 耦合变量更新策略：发散时的处理技巧

如果外循环迭代震荡，不要一味缩小收敛阈值，先看耦合设备出力变化历史。常见发散原因有三个：一是CHP的热电比不可行，设备模型超出可行域；二是P2G消耗电力和产气量折算失误；三是气网或热网求解失败但外层没感知。我采取的排查手段是每次外迭代打印耦合变量变化表，超过5次没收敛时，把阻尼系数降到0.2甚至0.1，并限制CHP出力变化步长不超过20%。另一个屡试不爽的小技巧是：先算不含P2G的工况，让CHP先稳定，再加入P2G，避免一上来就全耦合。这种逐步增加耦合度的启动方式，能让问题定位清晰很多。

5.4 常见问题速查表

下面这张表是我在整个项目周期里积累的排查经验，直接照着查能省大量调试时间：

最后再分享一个我后来一直在用的小技巧：不要等到整个程序跑完才画曲线，每一步迭代都把关键变量存到数组里，尤其是耦合设备的历史序列。有一次我怀疑外循环不收敛，画出CHP1的电出力序列后才发现它在一百多和一百六十之间来回跳，阻尼系数加少了。这种动态曲线比任何日志都直观，能帮你快速判断到底是数值问题还是模型问题。Matlab做这个特别顺手，一行plot就能解决，但前提是你在写循环时就把每一步的数据留下来。整个项目做到后面，你会发现能流计算本身不是障碍，障碍永远是边界条件、初值和单位这三个老熟人。把这套框架吃透，后面再往上叠加优化调度、故障分析或者多场景不确定性计算，都能少走一大半弯路。