搞过综合能源系统仿真的朋友应该都有过这种经历单张网的能流计算早就不是难题电力潮流、热网水力热力分析、燃气管网水力计算分别做各自都有成熟算法但一旦把这些网络放到同一个“区域综合能源系统”里事情马上就变得不一样了。配电网、区域热网和配气网通过热电联产CHP机组、燃气锅炉、电转气P2G等设备紧紧耦合在一起电网的负荷波动会改变CHP的出力CHP的燃气消耗又直接扰动气网的压力分布热网那边因为供水温度调整反过来影响电负荷。这时候传统的“各算各的”已经完全不适用必须在同一个框架内对电、气、热三类能流做联立求解——这就是计及多能耦合的区域综合能源系统电气热能流计算要解决的问题。我最近在Matlab里把这套计算完整实现了一遍电、气、热三网的稳态能流方程组统一用牛顿法迭代求解中间穿插各类耦合设备的能量转换关系。这篇文章把我从建模、编码、调试到算例验证的全过程整理出来希望对正在做综合能源系统仿真、运行优化或规划评估的同学有些实际帮助。1. 从“三张网各算各的”到“一张图联算”多能耦合带来的计算范式变化理解多能耦合为什么让问题变复杂是动手写代码之前必须先想明白的一件事。区域综合能源系统的“区域”二字决定了它的规模特征配电网、区域供热管网、配气管网节点数量通常从几十到几百不等规模不大但耦合关系密。这种系统里能量流的传递路径是环环相扣的。1.1 独立计算原本的适用范围把三类网络各自独立计算的逻辑先摆出来。配电网潮流计算给定负荷和电源注入用牛顿-拉夫逊法求解节点电压幅值和相角核心是节点的有功无功功率平衡方程。区域热网计算分两步走先做水力计算求管道流量和节点压力再做热力计算求供回水温度水力与热力交替迭代。燃气管网水力计算同样是为给定的气负荷求解节点压力和管道流量核心方程是管道流量与两端压力平方差之间的Weymouth关系。这三套方法单独拿出来都已经非常成熟工程上也有大量软件可用。如果三个系统边界上的负荷和注入量都是固定的常数独立计算完全够用这也是传统规划分析里习以为常的做法。问题在于综合能源系统通过耦合设备打破了这种“边界给定”的假设。1.2 一个真实会发生的耦合失效场景举个例子体会一下。冬季夜间园区的电负荷降低但热负荷反而升高。燃气锅炉大概率要顶上热负荷缺口锅炉耗气量增加气网末端压力随之下降。如果末端压力低到CHP机组的最小进气压力要求以下CHP就不得不降低出力它的电出力减少配电网的电压又会进一步走低。如果按固定负荷独立计算电网你完全看不到这个连锁反应独立算气网时你也无法确定CHP和锅炉的耗气量到底是多少因为它们是由电网和热网的状态共同决定的。这种跨网络的“蝴蝶效应”如果不把多能耦合放进计算模型里就只能在事后用经验系数去修正精度和可靠性都难以保证。1.3 三类方程在数学结构上的差异另一个难点来自三类方程的数学结构差异。电网功率方程是强非线性的变量是电压幅值和相角热网的水力方程包含管道流量的平方项热力方程则相对线性但流量的存在让水力计算与热力计算无法彻底解耦气网的Weymouth方程本身是二次的而且天然气管路可能存在双向流动方程的方向性需要特殊处理。把这三类方程统一放进一个非线性方程组F(x)0初值选择、量纲统一、收敛判据设计全都需要重新考虑。单看每一类网络都不难难的是让它们在一个求解框架里稳定地一起收敛。2. 电-气-热网络数学模型的推导从物理定律到可求解方程组统一求解的前提是把三类网络分别用一组可微的代数方程描述出来。这里我按电网、热网、气网的顺序把模型展开最后汇聚成一个统一变量向量。2.1 电力系统潮流方程电网采用最常用的极坐标牛顿-拉夫逊模型。对任意节点i有功和无功平衡方程写成P_i V_i \sum_{j \in N_i} V_j (G_{ij} \cos\theta_{ij} B_{ij} \sin\theta_{ij}) Q_i V_i \sum_{j \in N_i} V_j (G_{ij} \sin\theta_{ij} - B_{ij} \cos\theta_{ij})其中V和θ是节点电压幅值和相角G、B是导纳矩阵对应元素θ_ij θ_i - θ_j。PQ节点同时保留有功和无功两个平衡方程PV节点只保留有功方程平衡节点不参与计算。区域配电网里节点以PQ为主整体规模不大Matlab求解几十个节点的潮流毫无压力。这里需要特别留意如果电网某节点接了CHP或P2G那么该节点的功率注入就不再是常数而是气网变量耗气量或热网变量热出力的函数。这一行节点方程在雅可比矩阵中就会出现跨网络的偏导数项这正是统一求解与独立求解的本质区别。2.2 热力系统模型水力与热力热网模型比电网繁琐因为它有两层结构。水力计算包含两类方程。第一是节点流量连续性管道流入流出流量的代数和等于该节点净注入流量。第二是回路压力平衡闭合回路中所有管道压降的代数和为零。管道压降写成H_p R_p |m_p| m_pR_p是管道阻力系数m_p是管道质量流量。这里特意写成|m|·m而不是m²是为了保留流量方向信息。如果迭代过程中某条管道流量方向发生翻转|m|·m的导数依然连续不会导致雅可比矩阵突变。热力计算同样包含两类方程。管道沿程温度衰减T_{end} T_{ground} (T_{start} - T_{ground}) \exp\left(-\frac{\lambda L}{C_p m_p}\right)节点温度混合\sum_p m_{p,in} T_{p,out} m_{node} T_{node,mix}热负荷功率与供回水温度的关系\Phi_d C_p m_d (T_s - T_o)其中C_p是水的比热容T_s和T_o分别是供水和回水温度。热力方程和水力方程通过管道流量m_p相互耦合这也是为什么一些教材喜欢把热网计算拆成内部交替迭代的原因。但我在后面的实现里选择把所有变量统一放入一个向量而非做嵌套迭代理由会在第6节细说。2.3 天然气网络模型配气网的稳态模型相对简洁。管道流量用Weymouth方程描述f_g K_{wm} \sqrt{p_a^2 - p_b^2}K_wm由管道内径、长度、气体温度、压缩因子等参数决定。节点气体流量平衡方程\sum_{p \in i,out} f_{g,p} - \sum_{p \in i,in} f_{g,p} f_{i,load} - f_{i,source}压力平方差的写法让方程足够光滑但差值为负时开方无物理意义。程序实现中我会把方向从平方差符号中提出具体做法见最后一节。2.4 统一变量向量与统一方程组把三类网络的未知数放进同一个向量x [\theta; V; p_g; T_s; T_o; m]对应的方程F(x)0由三块构成电网节点的有功无功不平衡量、气网节点的气体流量不平衡量、热网节点的水力与热力不平衡量。变量总数在三网齐全的系统里通常在两百左右Matlab的稀疏矩阵可以轻松处理。这套变量的排列顺序一旦确定后续所有代码都必须严格遵循我习惯用一个“变量字典”结构体来记录每个变量在x中的起止索引避免硬编码带来的索引错乱。3. 耦合设备建模CHP、燃气锅炉与电转气的能量转换关系这一节是整个课题的核心也是“计及多能耦合”这几个字真正的落点。耦合设备是连接不同网络的实体接口建模的实质是把一种形式的能量转换为另一种形式并在两个或多个网络的平衡方程中同时体现出来。3.1 热电联产CHP机组CHP是区域综合能源系统里最常见的耦合设备消耗天然气同时输出电功率和热功率。定热电比模式常见于背压式机组P_e \eta_e f_{g,chp} LHV P_h c_m P_e其中η_e是发电效率LHV是天然气低位热值c_m是热电比。这个模式实现最简单两个方程就把电网出力、热网出力和气网耗气量绑定在了一起。变热电比模式常见于抽凝式机组P_e c_h P_h \eta_{total} f_{g,chp} LHV这种模式下CHP的电出力与热出力不再成固定比例可以在一个可行域内调节灵活性更好但建模时需要在方程中显式处理可行域边界。在能流计算里我倾向于把CHP的燃气消耗量f_g,chp作为未知变量而不是给定常数。这样电网中与CHP相连节点的有功注入P_e会随着f_g,chp变化气网节点平衡方程里也会出现与P_e或P_h相关的偏导项这些偏导就是要填入统一雅可比矩阵的耦合块。3.2 燃气锅炉与电转气P2G燃气锅炉比较简单天然气转化为热的效率通常在0.8到0.95之间\Phi_{gb} \eta_{gb} f_{g,gb} LHV它在气网侧制造耗气量在热网侧注入热功率本质上和CHP的一部分行为一致。电转气设备是目前研究热点方向正好相反电能转化为天然气f_{g,p2g} \eta_{p2g} P_{p2g} / LHVP2G在电网侧表现为一个可调电负荷在气网侧表现为一个气源。如果P2G的输入电功率是给定值那它的建模很简单相当于把一部分电能“转换”为气网注入流量如果它的输入功率还需要根据系统状态去优化就需要把P2G消耗的电功率也作为变量纳入统一的优化框架而不是纯能流计算。本文算例里P2G按给定功率处理重点放在多能流统一求解上。3.3 热泵与电锅炉热泵和电锅炉同样是电网到热网的耦合设备。热泵的制热性能系数COP通常在3左右输入1 kW电可以输出约3 kW热\Phi_{hp} COP \cdot P_{hp}电锅炉更简单效率接近1输入多少电就输出多少热。我写程序时把这几种设备统一抽象成“耦合设备对象”每个设备定义两个接口从源侧抽取的功率函数、向目标侧注入的功率函数。这样新增设备类型时网络求解核心完全不用动只需要注册新的接口函数扩展性会好很多。3.4 耦合变量如何进入统一雅可比矩阵这是程序实现中最容易出错的地方。先看一个分块结构示意J \begin{bmatrix} A_{ee} 0 0 A_{eh} A_{eg} \\ 0 A_{hh} A_{hg} 0 0 \\ 0 A_{gh} A_{gg} 0 0 \end{bmatrix}实际结构会因为耦合设备的接入位置而改变但核心思想一致每个耦合设备在组装雅可比矩阵时要做一次“连接登记”。假设CHP接在电网第i节点、气网第j节点、热网第k节点那么电网第i节点的有功方程对f_g,chp的偏导、热网第k节点的热功率方程对f_g,chp的偏导、气网第j节点的流量平衡方程对P_e或P_h的偏导都要在对应位置填进去。我建议把雅可比矩阵做成稀疏矩阵然后用辅助函数把设备参数映射到具体的行列索引而不是手动填一个全矩阵。稀疏矩阵既能提升求解速度也能避免大量零元素的存储浪费。4. Matlab程序实现统一变量字典、分块雅可比矩阵与收敛判据模型清晰之后代码实现就是水到渠成的事。但实现层面依然有几个关键设计点直接影响程序能不能稳定收敛、后续能不能方便扩展。4.1 整体程序架构我最终采用的文件结构如下main_IEGS.m主程序初始化数据、调用求解器、输出结果data_IEGS.m定义电网、热网、气网的节点与支路参数以及耦合设备参数build_residual.m根据当前x计算所有不平衡量Fassemble_jacobian.m组装统一稀疏雅可比矩阵Jnewton_solver.m牛顿迭代主循环utils/idx.m变量字典记录各物理量在x中的起止索引变量字典是关键设计。因为系统变量来自三个网络和多个耦合设备如果直接在代码里硬编码索引后面调试会非常痛苦。我的实现类似function idx init_idx(data) n_bus data.pf.N; % 电网节点数 n_node_g data.gas.N; % 气网节点数 n_node_h data.heat.N; % 热网节点数 n_pipe data.heat.N_pipe; % 热网管道数 idx.theta (1:n_bus); offset n_bus; idx.V offset (1:n_bus); offset offset n_bus; idx.p_g offset (1:n_node_g); offset offset n_node_g; idx.Ts offset (1:n_node_h); offset offset n_node_h; idx.To offset (1:n_node_h); offset offset n_node_h; idx.m offset (1:n_pipe); end有了idx结构体任何函数要修改某个变量时先从字典读出对应索引再对x进行赋值或取值逻辑清晰且不容易错。4.2 核心迭代代码牛顿迭代主循环写出来就是标准的骨架function [x, info] newton_solver(x0, data, opt) x x0(:); F0 build_residual(x, data); for k 1:opt.max_iter F build_residual(x, data); err max(abs(F)); if err opt.tol_abs opt.tol_rel * max(abs(F0)) info.converged true; info.iter k; return; end J assemble_jacobian(x, data); dx -J \ F; x enforce_bounds(x dx, data); if any(isnan(x)) warning(NaN encountered at iteration %d, k); info.converged false; return; end end info.converged false; endenforce_bounds用来保证变量在物理合理范围内比如气网节点压力不能为负、热网温度不能低于环境温度等。这一步看似简单实际能规避大量后期发散问题。4.3 初值设置策略初值对多能流计算的收敛性影响比单网潮流大得多。我的经验是电网侧V设1.0 puθ设0这是惯例。气网侧所有节点压力取气源压力的60%到80%作为初值。比如气源压力1.6 MPa所有节点初值取1.2 MPa。热网侧供水温度按设计值如90℃设置回水温度按设计值如50℃设置管道流量按热负荷除以供回水温差再换算成质量流量来估算。热网初值尤其关键。温度初值如果偏离实际太多管道温度衰减公式里的指数项可能出现数值溢出导致整个牛顿迭代在第一步就直接发散。4.4 收敛判据的设计电、气、热三个网络的物理量量纲完全不同不能共用一个绝对误差阈值。电网的有功不平衡量单位是MW气网流量不平衡量单位是kg/s热网的热功率不平衡量单位是MW或kW三者数量级可能差出好几个量级。我建议用混合判据\max\{ \|F_{elec}\|_\infty, \|F_{gas}\|_\infty, \|F_{heat}\|_\infty \} \epsilon_{abs} \epsilon_{rel} \|F_0\|_\infty这样既保证电网功率不平衡量收敛到了合理水平也保证气网流量不平衡量和热网功率不平衡量各自满足精度要求。单纯用一个统一的绝对阈值很可能出现某一类网络“以为收敛了”实际还差得很远的情况。5. 算例验证典型日场景下的能流分布与耦合效应光说不练没有说服力。我搭了一个小型区域综合能源系统把上面的算法完整跑了一遍并对比了独立计算与联立计算的差异。5.1 算例系统说明测试系统的结构如下电网修改版IEEE 33节点配电网基准电压12.66 kV基准功率100 MVA热网9节点区域供热网络设计供/回水温度90/50℃气网6节点配气管网气源压力1.6 MPa耦合设备1台定热电比CHP机组热电比1.2接在电网18节点、热网热源节点1、气网负荷节点41台燃气锅炉接在热网热源节点2、气网负荷节点41台P2G接在电网31节点、气网节点6设备关键参数整理成表设备参数数值CHP发电效率0.35CHP热电比1.2CHP天然气低位热值LHV38.9 MJ/m³燃气锅炉热效率0.9P2G综合转换效率0.6热网设计供水温度90℃热网设计回水温度50℃气网气源节点压力1.6 MPa5.2 典型场景与仿真结果设置两个典型日场景。场景A冬季白天电负荷较高热负荷中等场景B冬季夜间电负荷低、热负荷高。两个场景下的关键结果如下表指标场景A场景BCHP电出力1.8 MW1.4 MWCHP热出力2.16 MW1.68 MW燃气锅炉热出力4.5 MW7.2 MWP2G输入电功率0.5 MW0.2 MW电网18节点电压0.981 pu0.969 pu气网5节点压力1.21 MPa1.02 MPa统一迭代收敛次数11次14次收敛次数在预期范围内统一求解的迭代次数只比单电网牛顿法多出三四次说明把三网联立起来并没有造成不可接受的计算负担。更重要的是结果本身的物理合理性场景B中热负荷增大锅炉耗气量增加气网末端压力从1.21 MPa下降到1.02 MPa同时CHP的进气压力也受牵连电出力从1.8 MW降到1.4 MW18节点电压从0.981 pu降到0.969 pu。这个连锁反应链完全符合能量流动的物理直觉。5.3 与独立计算结果的对比我还做了另一组对照实验先给定CHP出力分别计算电网、热网、气网再把结果和统一联立求解结果比较。在场景B下独立计算假设CHP电出力固定为2.0 MW、锅炉耗气量固定为0.12 kg/s得到的气网5节点压力偏差约8%电网18节点电压偏差约0.015 pu。这个偏差放到工程里是不能忽略的8%的压力偏差可能直接改变管网压缩机或者调压站的运行点0.015 pu的电压偏差也会影响无功补偿设备的动作判断。这说明计及多能耦合并不是学术上的“锦上添花”。当系统内部存在CHP、P2G这类双向耦合设备且它们的运行状态随负荷动态变化时联立求解才能反映真实的能流分布独立计算只能在耦合较弱或耦合设备出力恒定的特殊情况下使用。6. 调试中的五个坑量纲、初值、索引与平方根符号处理代码写出来能跑通和任何工况下都能稳定跑通中间隔着无数个坑。我把自己调试过程中最折腾人的几个问题整理如下希望能帮你跳过这几段黑暗时光。6.1 量纲统一是最隐蔽的坑热网功率单位与比热容单位不匹配是我吃过最惨的亏。热网的热功率通常用MW但水的比热容C_p是4.182 kJ/(kg·K)如果直接套用Φ C_p·m·ΔT数字会差出1000倍。解决办法是在程序开头就固定一套单位制我采用的是全部SI制功率用W、质量流量用kg/s、温度差用K、压力用Pa。气网的Weymouth方程常数K_wm也必须和压力、流量单位严格配套文献里大量数据是以bar和m³/h为单位的复制过来之前务必重新换算这一条能救你一命。6.2 气网压力平方根问题牛顿迭代过程中节点压力的中间值可能让p_a² - p_b²变成负数。如果直接开方Matlab会返回复数或NaN后面整个线性方程组瞬间爆炸。建议把Weymouth方程改写成f_g K_{wm} \cdot sign(p_a^2 - p_b^2) \cdot \sqrt{|p_a^2 - p_b^2|}方向信息从平方差的符号中提取开方内永远非负。代价是在p_a等于p_b附近导数不光滑但实际稳定运行点极少刚好落在压力相等处工程上完全可接受。6.3 雅可比矩阵索引务必用有限差分校验自己手写解析雅可比矩阵时最常见的错误是耦合设备偏导数填错位置尤其在多设备多网络的场景下。我的做法是先用有限差分法算出数值雅可比矩阵与解析雅可比矩阵对比直到两者相对误差小于1e-5再跑正式迭代。这个方法看起来多花十几分钟实际上能节省几天的调试时间。6.4 热网独立迭代与统一迭代的取舍我最初把热网内部的水力-热力交替迭代放在外层统一迭代里结果经常出现外层迭代已经达到精度热网内部却还在振荡的情况。后来狠下心把所有热网未知量直接放进统一x向量不做任何内层嵌套收敛性反而明显改善。原因是内层迭代使用的容差即使设得很严也会给外层引入一个“数值噪声源”破坏牛顿法对非线性的全局处理能力。统一的牛顿框架虽然让雅可比矩阵更大但数值行为要干净得多。6.5 一遇奇异就检查设备运行模式场景B调试过程中我遇到过雅可比矩阵接近奇异的情况最后定位到是CHP定热电比模式导致电网方程与热网方程对耗气量变量的偏导成比例。解决办法是增加自由度把CHP的燃气消耗量从固定输入改成变量让它同时受电网功率方程、热网热功率方程和气网流量平衡方程共同约束。如果某个工况下设备出力触及上下限就固定它并去掉对应自由度这种处理方式在工程上非常常见。提示调试统一能流程序时如果雅可比矩阵发生奇异优先怀疑耦合设备的运行模式把“定值模式”改成“变量模式”往往能解决问题。最后分享一点实际操作中的体会。我现在做区域综合能源系统的多能流计算已经习惯性地把所有网络变量放在一个x里所有方程放在一个F里统一用牛顿法迭代。只要变量字典和雅可比矩阵的分工清楚这套框架的扩展比想象中容易得多——后续加设备、加约束、甚至改成动态仿真核心数学结构都不会动摇。如果让我给正在做类似工作的朋友一条最中肯的建议先用一个十几节点的小系统把统一迭代跑通再用有限差分验证雅可比矩阵最后再扩展到完整的大系统。能流计算这件事模型对、量纲对、索引对剩下的就只是时间问题了。