做区域综合能源系统的能流计算很多人的第一反应是把电网潮流、天然气网、热网分别算一遍最后再想办法拼起来。我最早也是这么干的——结果一个5节点电网加4节点气网加4节点热网的小系统顺序迭代了二十多轮还在震荡最后干脆不收敛。后来我把电、气、热三套方程和所有耦合设备方程放进同一个Matlab框架里做统一牛顿-拉夫逊求解才算把计及多能耦合这件事真正落地。这篇内容适合正在搭区域综合能源系统仿真平台、或者被多能耦合能流计算困扰的电气与能源专业同学参考里面所有坑都是我实际调代码时踩过的不是教科书上那种理论上应该收敛的写法。1. 电网、气网、热网为什么不能各算各的1.1 三套方程完全不是一个物种先得把三个子系统的数学形态看清楚否则后面的代码无从谈起。电网潮流用的是极坐标下的复数功率平衡方程核心未知量是节点电压幅值和相角方程里面全是三角函数。对某个节点i有功和无功平衡可以写成P_i V_i * sum_j V_j * (G_ij * cos(theta_ij) B_ij * sin(theta_ij)) Q_i V_i * sum_j V_j * (G_ij * sin(theta_ij) - B_ij * cos(theta_ij))这里的V和theta是标幺值量级大概在0.9到1.1之间theta在正负几度到十几度之间。这是一个典型的非线性方程组电力系统的人习惯用牛顿-拉夫逊法解。气网不一样。天然气管道的稳态流动通常用Weymouth方程描述管道流量和两端压力平方差的平方根成正比f_ij sign(p_i^2 - p_j^2) * C_ij * sqrt(|p_i^2 - p_j^2|)这个方程有两个特点第一它把压力平方当作基本变量天然是非线性的第二当管道两端压力相等时流量对压力的导数趋向无穷大。这意味着如果初值给得不好牛顿法一步就可能把压力迭代成负值然后整个矩阵崩溃。热网更麻烦它本身是水力-热力双层的。水力层要算管道流量分配满足节点流量连续和回路压降平衡热力层要算节点温度和热功率平衡温度在管道里还有指数衰减T_end (T_start - T_a) * exp(-lambda * L / (Cp * m)) T_a水力方程里有流量的一次项和二次项热力方程里有指数项两者通过流量决定温度衰减、温度决定热功率互相咬合。换句话说热网内部自己就是强耦合的你再把电网和气网叠上去问题复杂程度是乘法不是加法。1.2 耦合设备是把三张网缝起来的针脚区域综合能源系统和单纯三张网并列运行的本质区别在于存在同时连接多张网络的物理设备。这些设备的运行状态一旦变化会同步改变电网节点注入功率、气网节点耗气量、热网节点热功率。我在测试系统里用到的耦合设备有四类每一类的数学关系都很简单但架不住数量一多、耦合关系交错之后整体方程就变得非常硬设备类型能源转换关系典型模型燃气轮机CHP天然气 - 电 热P eta_e * G * LHVH c_m * P燃气锅炉天然气 - 热H eta_b * G * LHV电锅炉电 - 热H eta_eb * P_eb热泵电 - 热COP放大H COP * P_hp以燃气轮机CHP为例它烧的是气网送来的天然气发出来的是电同时通过余热回收产热。这三个量——耗气量G_chp、电出力P_chp、热出力H_chp——分别落在气网负荷方程、电网注入功率方程、热网热源方程里。如果哪一侧的方程变了另外两侧都得跟着动。这就是多能耦合四个字最直白的含义。它不是把三个能流计算结果放在一张表里叫耦合而是这些设备的输入输出变量必须同时在三个子系统的方程中联立求解一点都拆不开。1.3 顺序迭代在我这里的实测表现我第一版代码用的是顺序迭代法流程大致是先固定CHP热出力和耗气量解电网潮流得到电锅炉的用电功率再固定电功率解气网得到各节点压力再解热网更新热源功率然后回到第一步。每轮循环之间把上一轮的结果当作常数传给下一轮。理论上这可行但实际跑起来非常难受。我那个测试系统里电锅炉和热泵同时存在热负荷一增加理想情况下CHP的热出力应该跟着涨耗气量加大气网压力下降同时电网的平衡节点出力变化。然而顺序迭代里每一轮都在用上一轮旧的值去计算这一轮的耦合关系反馈天然慢了一个周期。在耦合强度高的工况下残差会在10^-1量级来回震荡用松弛因子压能勉强收敛但收敛速度慢到没法看。顺序迭代还有一个隐藏问题它没有统一的收敛判据。电网侧用功率偏差判断气网侧用压力偏差判断热网侧用温度偏差判断三者不在同一个尺度上到底以谁为准实际调试的时候经常出现电网已经收敛、热网还差一大截或者反过来。这种看起来收敛了但其实没全收敛的状态是做综合能源系统仿真最头疼的事。2. 统一能流计算的数学模型变量、方程与耦合关系放弃顺序迭代之后我转向了统一求解框架。思路很朴素把电网、气网、热网和所有耦合设备的方程全部写进一个残差向量把所有未知量全部写进一个变量向量然后用牛顿-拉夫逊法一次性解。2.1 电网方程极坐标牛顿-拉夫逊电网部分采用最常见的极坐标潮流方程。节点分为三类平衡节点、PV节点、PQ节点但我在统一框架里做了个简化——把所有节点的电压幅值和相角都放进未知量向量然后对平衡节点额外加两个简单等式约束V 1.0theta 0。这样做虽然变量数量多了一点但代码结构非常统一不需要在迭代过程中反复判断节点类型。电网方程对变量的偏导数很规整和单电网潮流计算完全一样。不同点在于常规潮流里CHP并网节点和电锅炉节点的注入功率是给定的常数在统一能流里这些P和Q本身是未知量由耦合设备方程决定。所以电网方程的雅可比需要额外增加对耦合设备变量比如P_chp、P_eb的偏导数项这部分在后面的矩阵结构里会细说。2.2 气网方程Weymouth 稳态模型气网考虑的是稳态流动节点流量平衡方程为sum_{j in i} f_ij - L_i 0其中L_i是节点负荷包括常规气负荷和所有接在该节点上的耦合设备耗气量。管道流量用Weymouth方程计算方向通过sign函数处理。我没有在首版代码里加入压缩机模型因为区域级的天然气网络通常不需要长距离增压如果你要处理带压缩机的系统需要额外引入压缩比和压缩机耗气量方程原理上不冲突只是变量和残差各多一组。气网方程求导的时候最容易出问题的地方是Weymouth方程在零流量点附近。导数为无穷大意味着雅可比矩阵在这个区域数值上很险我的处理方式是让所有非气源节点压力初值比气源节点低一点比如气源节点p1.0标幺其他节点p0.95从初值上避开零流量点。2.3 热网方程水力平衡叠加热力平衡热网统一采用供热-回水双温度网络模型。水力方程包括两部分节点流量连续A * m m_load 回路压降平衡B * K * m .* |m| 0其中A是节点-管道关联矩阵B是回路关联矩阵K是管道阻力系数矩阵。热力方程包括节点功率平衡、节点混合温度方程和管道温降方程。这里要特别注意一个设计决策供水温度和回水温度必须作为两组独立的变量同时出现在未知量向量里。很多人第一次建模时只算供水侧回水侧用一个固定温度替代这在单一热源系统里勉强可用但CHP和电锅炉共同供热的时候不同热源的供水温度不一样节点混合之后回水温度也会随负荷分布变化固定回水温度会直接算错热功率分配。节点混合温度的方程是T_mix sum_k (m_k * T_k_out) / sum_k m_k管道温降则用指数衰减公式。这两个方程合在一起构成了热网残差的主体。2.4 耦合设备方程如何进全局变量我在变量向量里给每一台耦合设备都分配了专属变量比如CHP的三变量组P_chp、H_chp、G_chp、电锅炉两变量组P_eb、H_eb、燃气锅炉两变量组G_gb、H_gb。对应的残差方程是CHPP_chp eta_e * G_chp * LHV CHPH_chp c_m * P_chp 电锅炉H_eb eta_eb * P_eb 燃气锅炉H_gb eta_gb * G_gb * LHV这些变量被分别填入三个子系统的方程里。电网方程里看到P_chp和P_eb气网方程里看到G_chp和G_gb热网方程里看到H_chp、H_eb和H_gb。一个耦合设备的变量横跨多个子系统方程这才是计及多能耦合在数学上的真正落地。关键原则是变量数和方程数必须严格一致。每引入一组耦合设备变量就必须引入一组对应的耦合设备方程否则雅可比矩阵是奇异的Matlab解线性方程组时会直接报错或者算出离谱的结果。3. 统一牛顿-拉夫逊的工程取舍3.1 两种求解策略的对比把两种方法摆在一起对比差异非常明显对比项顺序迭代法统一求解法编程复杂度低可复用单能流程序高需要统一变量和残差框架收敛性弱耦合系统尚可强耦合易震荡稳健牛顿法二次收敛雅可比矩阵规模三个小矩阵分别求逆一个大矩阵统一求解收敛判据不一致电网/气网/热网各算各的统一残差向量一个判据对初值依赖较低可分系统分别给初值较高但给子系统合理初值即可适用场景弱耦合、快速估算强耦合、高精度、优化调度底座我个人的实践体会是如果只是做一个非常粗略的估算比如热负荷和电负荷只通过CHP关联且耦合强度很低顺序迭代还能凑合用。但一旦涉及电锅炉、热泵、P2G这类强耦合设备顺序迭代基本就是浪费时间。统一求解的编程工作量前期大一些但换来的是稳定收敛和一个可以继续往上面堆各种扩展模块的干净框架。3.2 统一求解的雅可比矩阵长什么样统一求解的核心是构造一个大的非对称雅可比矩阵J它的行对应残差方程列对应未知量。结构上可以分块理解电网变量 气网变量 热网变量 耦合设备变量 电网方程 [Jee] [0] [0] [Jec] 气网方程 [0] [Jgg] [0] [Jgc] 热网方程 [0] [0] [Jhh] [Jhc] 耦合方程 [Jce] [Jcg] [Jch] [Jcc]这个矩阵里对角线上的Jee、Jgg、Jhh分别是电网、气网、热网内部方程对自身变量的偏导数基本就是单能流计算的雅可比矩阵。真正的多能耦合体现在非对角块上Jec电网功率平衡方程对CHP电出力、电锅炉功率等耦合变量的偏导Jgc气网节点流量平衡方程对CHP耗气量、燃气锅炉耗气量的偏导Jhc热网节点热功率方程对CHP热出力、电锅炉热出力等的偏导Jce耦合设备方程比如CHP的P eta_e * G * LHV对电网电压、气网压力等系统变量的偏导。矩阵规模大概是几十到几百阶对Matlab来说直接用小规模稠密矩阵求解没有任何压力。我建议先把J当普通矩阵用反斜杠解等到系统规模超过两百个变量再去考虑sparse存储过早优化反而增加调试难度。3.3 量纲归一化是统一求解的前提这是我在整个项目里踩过最重的一坑。一开始我把电网的标幺电压0.95~1.05、气网节点的真实压力几十万Pa、热网管道流量几kg/s、温度300多K直接堆进同一个变量向量结果牛顿法第一步就算不下去雅可比矩阵的逆基本是垃圾。原因很简单不同变量的量级差了好几个数量级。压力变量是10^5级别电压变量是1级别迭代步长dx -J\F里压力方向的步长天然比电压方向大好几个数量级电压被压力绑架了几乎动不了。解决方法是给每个子系统做标幺化处理。我选的基准是电网功率基准S_base10MVA电压基准取节点额定电压气网压力基准p_base1MPa热网流量基准m_base1kg/s温度基准T_base373.15K。所有变量在进入迭代之前都除以各自的基准值迭代完成后在输出阶段再换算回物理单位。归一化之后雅可比矩阵的条件数直接从10^10级别降到10^3左右牛顿法的迭代次数从十几步降到五六步。这个工作虽然看起来只是把单位换算了一下但实际上是整个程序能不能收敛的分水岭。4. Matlab代码骨架索引映射、分块组装与阻尼牛顿4.1 用struct组织三个子系统的参数Matlab里实现这个框架我强烈建议先用struct组织数据不要一上来就写class。struct的字段访问清晰在调试的时候可以直接breakpoint看数据而且对小规模系统没有性能问题。我最初给三张网各建了一个struct大概长这样% 电网 grid.bus [ % bus type Pload Qload V0 theta0 1 3 0 0 1.0 0 2 1 2.0 1.0 1.0 0 3 1 3.0 1.5 1.0 0 ]; grid.branch [...]; % 线路参数: from, to, R, X, B % 气网 gas.node [...]; % 节点编号, 气源/负荷标记, 负荷流量 gas.pipe [...]; % 管道编号, from, to, C_ij % 热网 heat.node [...]; % 节点编号, 类型(热源/负荷), 热负荷 heat.pipe [...]; % 管道编号, from, to, 长度, 直径, 阻力系数等程序稳定运行之后如果你想继续扩展成面向对象结构把这三个struct包进一个系统类里那也是水到渠成的事。但一开始就上class会引入很多和数据组织无关的复杂度不利于排查数学问题。4.2 全局变量编号与idx映射表统一求解中所有未知量都放在一个列向量x里那么某个变量在x的哪个位置就必须有一个全局映射表。我习惯用Matlab的struct来存这个索引表每次迭代都通过它访问变量nB size(grid.bus, 1); nGas size(gas.node, 1); nPipeH size(heat.pipe, 1); nNodeH size(heat.node, 1); idx.V 1 : nB; idx.theta idx.V(end) 1 : idx.V(end) nB; idx.pGas idx.theta(end) 1 : idx.theta(end) nGas; idx.mH idx.pGas(end) 1 : idx.pGas(end) nPipeH; idx.Ts idx.mH(end) 1 : idx.mH(end) nNodeH; idx.Tr idx.Ts(end) 1 : idx.Ts(end) nNodeH; % 耦合设备变量接着排 idx.Pchp idx.Tr(end) 1 : idx.Tr(end) nChp; idx.Hchp idx.Pchp(end) 1 : idx.Pchp(end) nChp; idx.Gchp idx.Hchp(end) 1 : idx.Hchp(end) nChp; idx.Peb idx.Gchp(end) 1 : idx.Gchp(end) nEb; idx.Heb idx.Peb(end) 1 : idx.Peb(end) nEb;这段代码是示意性的实际系统里还要加上燃气锅炉变量但结构完全一样。有了这个idx映射表后续组装修残差和雅可比矩阵的时候就非常省心再也不会出现这个变量到底在第几列的混乱。我建议把idx映射表的构建单独写成一个函数任何地方要改变量顺序只改这一处。4.3 残差和雅可比怎么分块组装组装的流程是固定套路先开一个全零的残差向量F和一个全零的雅可比矩阵J然后依次填入电网、气网、热网、耦合设备的贡献。我这里用一段伪代码说明function [F, J] assembleSystem(x, idx, net) N length(x); F zeros(N, 1); J zeros(N, N); % 电网残差和电网-电网/电网-耦合分块 [Fe, Je] electricResidual(x, idx, net.grid); F(idx.equation.e) Fe; J(idx.equation.e, idx.all) Je; % 气网残差 [Fg, Jg] gasResidual(x, idx, net.gas); F(idx.equation.g) Fg; J(idx.equation.g, idx.all) Jg; % 热网残差 [Fh, Jh] heatResidual(x, idx, net.heat); F(idx.equation.h) Fh; J(idx.equation.h, idx.all) Jh; % 耦合设备残差 [Fc, Jc] couplingResidual(x, idx, net.couplings); F(idx.equation.c) Fc; J(idx.equation.c, idx.all) Jc; end这里有个细节值得强调每个子系统函数返回的雅可比不光是子系统内部方程对内部变量的导数还包括对耦合设备变量的导数。比如electricResidual里计算电网有功方程对CHP电出力的偏导这个偏导可能是1因为P_chp直接以正项出现在节点功率平衡里也可能是0取决于CHP是否接在这个节点。用解析方式求导很简单但你必须把y——也就是耦合变量出现在哪些方程、以什么系数出现——搞清楚。实在没把握的时候对一个局部的小规模子系统用中心差分验证一次雅可比的正确性是非常值得做的。牛顿迭代主循环就是标准的写法外面套一个阻尼策略x x0; for k 1:maxIter [F, J] assembleSystem(x, idx, net); if norm(F, inf) tol break; end dx -J \ F; % 阻尼牛顿步长减半直到残差下降 alpha 1; while alpha 1e-4 xtry x alpha * dx; [Ftry, ~] assembleSystem(xtry, idx, net); if norm(Ftry, inf) norm(F, inf) break; end alpha alpha * 0.5; end x x alpha * dx; enddamping的关键是保证每次迭代残差都在下降只要残差不降就把步长砍半重试。这个方法对初值不太离谱的系统都能救回来。4.4 初值不用怕阻尼牛顿帮你兜底统一求解的初值并没有想象中那么玄各子系统用各自的合理初值组合起来就行。我常用的初值方案子系统初值设置电网所有节点V1.0theta0平衡节点V1.0气网气源节点p1.0标幺其余节点p0.95热网供水温度所有节点设为设计供水温度回水温度设为设计回水温度管道流量按热负荷估算耦合设备各设备按设计工况初始CHP取额定出力的50%锅炉取额定热负荷的60%燃气轮机CHP的初始耗气量一定不要超过气源供气能力否则牛顿法第一轮就会把气网压力迭代成负值整个矩阵直接崩掉。我曾经因为初值给得太满反复出现压力为负的问题后来把所有非气源节点压力初值统一降到0.95标幺同时把CHP耗气量初值压到额定值的40%问题立刻消失。这不是数学问题是纯经验问题。5. 测试算例与结果5.1 测试系统5节点电网加4节点气网加4节点热网为了验证框架我搭了一个小规模区域系统参数刻意贴近真实工程尺度电网5节点配电网基准容量10MVA节点1作为平衡节点节点2和3接常规电负荷节点4接电锅炉节点5接CHP。气网4节点树状网络节点1是气源节点2接普通气负荷节点3接CHP节点4接燃气锅炉。热网4节点供热网络热源节点是CHP、电锅炉和燃气锅炉的公共接入点节点2和3接热负荷节点4作为回水混合节点。耦合设备参数取典型的工业值燃气轮机CHP电效率0.35、热电比0.6电锅炉电热效率0.98燃气锅炉热效率0.9。常规负荷设定为电网总负荷约8MW热负荷约6MW气负荷约2MW。所有功率以10MVA作为电功率基准气网压力以1MPa为基准温度基准取373.15K。5.2 收敛记录和关键结果统一牛顿法在这个算例里的表现让我印象很深。初值按上面的方案给定第1次迭代残差从10^-1级别跳到10^-3第4次迭代已经到10^-5第6次迭代残差降到10^-7以下满足收敛判据。相比之下我之前用顺序迭代法在这个系统上跑了二十多轮还在10^-1量级震荡差别真的非常大。最终的关键结果如下电网节点电压幅值(pu)相角(deg)1平衡1.0000.0020.994-1.5230.988-2.414电锅炉0.982-2.555CHP0.991-1.78气网节点压力(MPa)1气源1.00020.9283CHP接口0.8724燃气锅炉0.903热网供水温度在88到90℃之间回水温度在48到50℃之间各负荷节点的热功率平衡都精确满足。CHP在收敛点的电出力约3.2MW热出力约1.92MW耗气量折算约0.122kg/s电锅炉和燃气锅炉按经济调度逻辑分别承担剩余热负荷——这些工作点完全是由耦合方程和各子系统方程同时决定的没有人为指定这正是计及多能耦合的意义所在。5.3 和拆开单算的差异对比作为对照我用拆开单算的方式处理同一个算例先固定CHP电出力为3.2MW、热出力为1.92MW、耗气量为0.122kg/s把耦合设备全部当成常数然后三个子系统分别求解。结果与统一求解对比对比量统一求解拆开单算偏差电网节点4电压0.982pu0.979pu0.3%气网节点3压力0.872MPa0.841MPa3.6%热网节点2供水温度89.2℃88.4℃0.8KCHP所在节点电压0.991pu0.996pu0.5%表面上看电压偏差只有0.3%到0.5%好像不算大但气网节点压力偏差到了3.6%这个量级在工程上已经不能忽略了。更关键的是这还只是解耦一次的结果在强耦合工况下拆开单算产生的偏差会通过设备方程逐轮放大最后根本对不上。做区域综合能源系统分析的人常说单能流算出来只能当参考原因就在这里。6. 实操中的四个坑和后续扩展方向6.1 四个最容易被忽略的坑第一个坑是单位制混乱。我算天然气耗量的时候一会儿用标准立方米每小时一会儿用kg/s热值一会儿用MJ/m³一会儿用MJ/kg中间换算错了三次才把量纲对齐。建议所有内部计算统一用SI单位制只在输入输出阶段做换算并且把换算系数集中放在一个配置文件里。第二个坑是热网温降公式里的除零问题。管道流量m等于0的时候exp(-lambdaL/(Cpm))直接崩溃。实际工程中流量为0的管道不应出现在运行工况里但算例里如果某根管道因为负荷波动趋近零流量程序就会挂掉。我的处理是在管道流量初值设定时加一个下限比如0.01kg/s同时避免在迭代中让流量变号。第三个坑是雅可比矩阵的行列索引写反。分块组装的时候如果残差函数里某个变量的索引填错矩阵的非对称性会变得很奇怪收敛速度会急剧下降。我在每个子系统函数里加了assert检查确保返回的雅可比矩阵尺寸和残差长度一致再用一个小算例的数值差分验证一次各分块导数之后就再没出过这类问题。第四个坑是气网初值避开零流量点。Weymouth方程在零流量处导数无穷大如果所有节点压力初值都一样管道流量初值恰好为0雅可比矩阵奇异。我的做法是让非气源节点压力初值比气源节点低0.05标幺从根上避开这个问题。6.2 这个框架还能往哪些方向扩展统一求解框架的好处是任何新的设备或子系统都可以按变量-方程的套路往里面加。储能设备就是给蓄电池/蓄热罐各加一组状态方程在稳态能流里体现为不同时段之间的耦合优化调度可以把能流计算作为内层约束外层用内点法或粒子群寻找经济最优工作点如果要做准动态分析把气网和热网的时间常数加进去本质上还是在这个统一残差框架上扩展状态变量。我个人实际使用中最满意的一点是这个框架让我从挨个设备调参数、反复试错的状态里解放了出来。以前改一个热负荷值要重新跑三个流程再人工核对现在只需要改输入数据统一求解一遍所有工作点一次性到位。如果你也在做区域综合能源系统的能流计算我的建议是不要被统一求解前期建模型的工程量吓退——那段调试期换来的是后续所有分析和优化工作的稳定地基。