做区域综合能源系统仿真的朋友几乎都绕不开“电气热能流计算”这道坎。我接这个课题时正处于几种状态叠加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 - θ_jV_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²|)这里π是节点气压单位常用MPaC_mn是管道常数由管径、摩擦系数、气体温度、压缩因子等决定sign函数用来表达流量方向。如果π_m大于π_n流量从m流向nsign取正反之取负。直接对气压而不是气压平方求导方程会产生根号项和方向项数值上不太稳定所以我在程序中把求解变量设置成气压平方π²让方程变成线性平方差形式牛顿法的雅可比矩阵更规整。对每个不含气源的节点需要建立流量平衡方程流入该节点的流量总和加上气源注入或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平衡2PV3PQ这里默认平衡节点为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注意外循环更新策略很关键。我踩过的坑是直接把新设备出力整体替换旧值遇到强耦合场景经常发散。后来给更新加了一个阻尼系数alpha0.5相当于给迭代加低通滤波抵掉高频振荡收敛性大幅改善。每个子网求解器的输出反过来会成为另一个子网的边界因此子网内部也要做保护处理比如气网节点气压越界时及时返回错误码不要让外层循环背着错误继续跑。4. 算例验证与结果分析4.1 测试系统参数与耦合设备接入为验证代码正确性我搭了一个中等规模算例电网采用修改的13节点辐射配电网气网为6节点环状结构热网为6节点辐射状结构三个网络通过2台CHP、1台燃气锅炉、1台P2G耦合。网络参数和负荷数据都放在文本文件中。设备参数整理如下表设备额定电功率/kW热电比效率/%连接说明CHP15001.35电36、热49电网5节点、气网3节点、热网2节点CHP23001.20电34、热41电网8节点、气网4节点、热网4节点燃气锅炉600—89气网5节点、热网5节点P2G200—62电网10节点、气网2节点约束条件为电网平衡节点电压标幺值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 常见问题速查表下面这张表是我在整个项目周期里积累的排查经验直接照着查能省大量调试时间现象可能原因处理方式电网潮流发散雅可比索引错位、PV节点越限未处理用符号微分校验雅可比越限节点转PQ气网压力出现负值初值太差或管道方向计算错误检查sign处理初值取0.8倍设计压力热网温度异常低或为负质量流量初值不合理、散热指数溢出先做水力预计算限制温度更新范围外循环振荡耦合变量更新增益过大减小阻尼系数限制单步变化幅度收敛但结果不合理单位换算有误或设备可行域未校验逐项核对参数单位增加设备出力约束最后再分享一个我后来一直在用的小技巧不要等到整个程序跑完才画曲线每一步迭代都把关键变量存到数组里尤其是耦合设备的历史序列。有一次我怀疑外循环不收敛画出CHP1的电出力序列后才发现它在一百多和一百六十之间来回跳阻尼系数加少了。这种动态曲线比任何日志都直观能帮你快速判断到底是数值问题还是模型问题。Matlab做这个特别顺手一行plot就能解决但前提是你在写循环时就把每一步的数据留下来。整个项目做到后面你会发现能流计算本身不是障碍障碍永远是边界条件、初值和单位这三个老熟人。把这套框架吃透后面再往上叠加优化调度、故障分析或者多场景不确定性计算都能少走一大半弯路。