综合能源系统电气热能流计算：Matlab实现与收敛策略解析

如果你在复现“计及多能耦合的区域综合能源系统电气热能流计算”这类题目搜一圈会发现论文很多但能直接落地到Matlab的工程化经验非常零散。电、气、热三套系统单独建模都很成熟可一旦加了多能耦合问题立刻变成“怎么把三个异质网络放进同一个能流计算框架里”并且还要保证数值收敛。我前段时间完整跑过一套Matlab实现从方程建模到收敛调试都过了一遍这篇文章就把整个思路和代码骨架整理出来适合正在做综合能源系统方向研究、或者接活需要复现这类算例的工程师参考。文章不绕弯子直接从能流方程开始讲到求解策略、Matlab实现细节最后给一个实际算例和踩坑记录。1. 先搞清楚要算什么电、气、热三类网络各自的能流方程能流计算本质是求解一组非线性代数方程让系统达到稳态平衡状态。理解每一步的目的是关键。恰恰是那些基础概念现在想了想早期走弯路的原因就是混淆了各类网络的物理本质。1.1 电网潮流模型从经典牛顿法出发电网潮流采用最常见的极坐标牛顿-拉夫逊模型。对节点 i有功和无功平衡方程分别为P_i - V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij) 0Q_i - V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij) 0其中 θ_ij θ_i - θ_j。电网模型中状态变量是各节点的电压幅值 V 和相角 θ平衡节点slack bus给定电压幅值和相角PV节点给定有功和电压幅值PQ节点给定有功和无功。耦合设备接入后等效为节点注入功率的修正CHP在电网节点表现为有功注入和无功由PV或PQ设定决定电锅炉和P2G则表现为额外有功负荷。这一层的Matlab实现没什么特别关键是把导纳矩阵 Y 用稀疏矩阵存储求解线性修正量时直接列主元分解即可。值得注意的是当电网只是整体能流计算的一个子模块时slack节点的选取会对耦合迭代收敛产生明显影响。我实际测试时把slack选在靠近大负荷区域和选在馈线末端牛顿法的收敛轨迹差异很大因为耦合设备改变了局部功率平衡模式。1.2 气网能流模型Weymouth方程与节点平衡天然气网络的能流计算核心是管道的稳态流量方程。工程上最常用的是Weymouth公式描述管道两端压力与流量的关系。如果直接用压力 p 建模方程为f_ij C_w * sqrt(p_i² - p_j²)但这里有个坑该式只适用于 p_i ≥ p_j即默认正向流动。实际迭代过程中压差可能随时变号直接套公式会让雅可比矩阵在零压差附近出现数值跳变。为了避免这个问题我在Matlab代码里把状态量换成压力平方 π p²并把流量方程改写成带符号的形式f_ij sign(π_i - π_j) * C_w * sqrt(|π_i - π_j|)这个形式在正反向流动时都可用但 sqrt 在 π_i π_j 处不可导。做法是在小邻域内做线性化处理见后面的代码片段。气网节点平衡方程为A * f L_g 0其中 A 是节点-管道关联矩阵L_g 是节点净负荷正值流出。对于气源节点通常给定压力作为边界条件对于纯负荷节点给定流量。在综合能源系统中气网节点并不单纯挂负荷还连接燃气锅炉、CHP和P2G设备。CHP和燃气锅炉消耗气网流量是正值气体负荷P2G则向气网注入天然气表现为负值负荷。因此气网能流计算必须与之耦合迭代不能单独求解完气网再去搭接口。1.3 热网能流模型水力与热力耦合的两个层次热网能流计算是三个网络中最容易让新手懵掉的部分因为它包含两个层次的物理方程水力工况计算和热力工况计算这两者互相耦合。水力工况描述热水在管道中的流量分配。节点流量守恒方程为A_s * m_s m_qA_r * m_r -m_q其中 m_s 是供水管网各管道流量m_r 是回水管网各管道流量m_q 是节点注入/流出的质量流量。管道压降由达西-魏斯巴赫公式简化而来Δp_ij K_ij * m_ij * |m_ij|K_ij 与管长、管径、摩擦系数有关。闭合回路压降必须满足基尔霍夫第二定律任意基本回路中各管道压降的代数和为零。求解水力工况就是联立节点流量守恒、回路压降为零这些约束通常用一个独立的牛顿法解出流量分布 m。热力工况描述温度场。节点热功率方程φ_i C_p * m_i * (T_s,i - T_r,i)其中 C_p 是水的比热容。管道沿程温降方程为T_end T_a (T_start - T_a) * exp(-λ * L / (C_p * m))λ 是管道单位长度传热系数L 是管长m 是流过该管道的质量流量T_a 是环境温度。节点混合温度由流量加权平均决定。热网求解时既可以水力、热力交替迭代也可以联立求解。由于温度方程中管道温降呈指数形式联立求解的雅可比矩阵容易出现病态所以很多实际代码采用“先解水力再对温度做拓扑前推回代”的顺序结构。1.4 耦合设备把三种能源串起来的接口节点电-气-热三个网络并不是各自独立运行的连接它们的设备就是耦合方程所在。这个项目里我主要建模了四类设备热电联产机组CHP消耗天然气同时输出电功率和热功率。运行区域可以用定热电比或可变热电比简化建模 P_e η_e * F_chp * LHVφ_h η_h * F_chp * LHV定义热电比 α φ_h / P_e。这台设备直接将气网天然气消耗、电网电输出和热网热源注入三者耦合起来。燃气锅炉消耗天然气输出热功率 φ_h η_boiler * F_gboiler * LHV它连接气网和热网是CHP热出力不足时的补充热源。电锅炉消耗电功率输出热功率 φ_h η_eb * P_elec它连接电网和热网常作为灵活负荷参与调节。P2G电转气消耗电功率向气网注入天然气 F_p2g η_p2g * P_elec / LHV它连接电网和气网方向是电→气。耦合设备方程是整个系统的“接口”求解策略的关键就在于如何让这些接口变量在电、气、热三套方程之间一致收敛。如果只算单一网络的潮流其实很简单一旦把这些接口变量都打开系统状态变量维度会从几十扩展到上百这才是研究的核心难点。2. 统一求解还是交替迭代初值、雅可比矩阵与收敛性权衡面对多能耦合能流计算主流方案无非两种统一求解法和交替迭代法。这个选择直接决定了Matlab代码的整体结构所以值得单独拿出来讲透。2.1 统一求解法的数学框架统一求解法的思路是把电网、气网、热网的所有状态变量都放进同一个未知向量 x 里把所有平衡方程拼接成一个大残差向量 F(x)然后用牛顿法同时求解x [θ_i, V_i, π_i, m_ij, T_s,i, T_r,i]F(x) [电网有功残差; 电网无功残差; 气网节点流量残差; 热网水力残差; 热网热力残差]J * Δx -F(x)这种方法的优势是理论上耦合最强、收敛阶数最高不存在接口变量交替迭代带来的收敛延迟。缺点也显而易见雅可比矩阵规模大分块结构复杂而且电、气、热三个网络的量纲差异巨大——电压是kV量级压力平方可能到MPa量级温度是几十到上百度。如果不做归一化处理雅可比矩阵的条件数会非常差迭代几步就发散了。实际建模时要在统一牛顿法前先完成“归一化”这一步电压用标幺值压力用基准压力归一温度用设计供回水温差归一。这样把各状态变量都压缩到0.1~10的量级雅可比矩阵的奇异值分布才比较合理。2.2 交替迭代法的实现逻辑交替迭代法采用“顺序求解-接口更新”策略先解电网潮流利用电网结果更新CHP电出力、电锅炉电负荷、P2G电耗再把这些量代入气网解气网能流得到天然气消耗和压力分布接着更新热网热源出力解热网水力-热力工况然后检查接口变量变化量是否小于收敛阈值如果不满足回到电网重新求解。这种方法的代码组织最清晰每个网络保留独立求解器便于调试和扩展设备模型。缺点也很明显强耦合工况下简单顺序迭代可能陷入震荡甚至发散。因为CHP的耗气量影响气网压力气网压力又反过来限制CHP进气量CHP出力变化再影响电网潮流这个循环如果不加阻尼或加速因子收敛非常慢。2.3 我的选择分块统一牛顿法及其理由我在最终代码里没有采用纯粹的统一求解或纯粹的顺序迭代而是用了“分块统一牛顿法”。核心思路是电网和气网联立成一个大牛顿系统热网作为耦合子系统在每次外层迭代内部求解。具体来说外层统一求解的状态变量是电网的 θ、V 和气网的 π外加CHP和锅炉的热出力接口变量。热网子系统给定热源出力和热负荷后内部先解水力方程得到流量分布 m再做温度前推回代得到 T_s、T_r返回热平衡残差给外层牛顿法。这样做的理由有三点第一电网和气网之间的耦合最紧密CHP、燃气锅炉、P2G都直接连接这两个网络统一求解能显著提高收敛性第二热网的温度方程指数非线性较强单独求解更容易控制数值稳定性放在内部循环处理外层维度降低第三从代码复现的角度看热网作为子函数维护后续更换热网拓扑只需要改子函数不需要动主迭代骨架。2.4 初值设计与量纲归一化无论选择哪种求解框架初值问题都绕不开。我踩过最深的坑就是气网压力初值给得太随意导致第一步迭代就撞到Weymouth方程的不可导邻域。电网初值采用平启动全部PQ节点电压幅值给1.0相角给0。气网初值给统一基准压力标幺值1.0对应实际压力2MPa这样管道压差为零点附近的初始流量也为零虽然距离真实工况较远但数值上不会爆。热网初值比较讲究供水温度给设计值90°C回水温度给50°C管道流量先通过线性网络近似分布一次再进入牛顿循环。这样热力前推回代的前几步温度残差已经比较小外层接口变量更新也更平滑。此外所有状态变量在进入统一雅可比矩阵前都必须除以各自的基值。我在代码中用一个结构体保存基值信息base.Sbase 10e6; % 功率基准 10 MVA base.Vbase 12.66e3; % 电网电压基准 12.66 kV base.pbase 2e6; % 气网压力基准 2 MPa base.Tbase 40; % 热网温差基准 40 °C归一化之后牛顿法修正量的限幅阈值才能统一设置例如限制最大修正量为0.5否则压差修正量可能比其他变量大几个数量级导致迭代震荡。3. Matlab代码的骨架数据结构、矩阵组装与主迭代循环这一部分直接给可落地的代码设计和关键片段。我假设你已经具备基本的Matlab编程能力但可能没有大规模稀疏矩阵组装和牛顿法调试经验所以会连原理带代码一起讲。3.1 用struct组织三类网络的拓扑数据Matlab面向过程的编程风格对这类多网络耦合项目非常合适我没有用class而是用一组struct保存网络数据。这样做的好处是代码透明、调试方便、与论文算法伪代码的对应关系直观。电网部分pn.bus [... % bus type Pd(MW) Qd(MVar) ... ]; pn.branch [... % 支路数据 % from to R X B ]; pn.Y makeYbus(pn); % 稀疏导纳矩阵气网部分pg.node [... % 节点编号、类型(0气源/1负荷)、压力基准 ]; pg.pipe [... % 起始节点、终止节点、Weymouth系数C_w、管长L ]; pg.A makeIncidence(pg); % 节点-管道关联矩阵热网部分ph.node [... % 编号、类型(0热源/1热负荷)、热功率需求、设计供回水温度 ]; ph.pipe [... % 起始节点、终止节点、K系数、长度L、传热系数lambda ];把数据和拓扑分开保存后续组装残差函数时直接通过函数句柄传递这些struct可以避免全局变量的滥用。实测下来这种结构在扩展P2G或储能设备模型时非常顺手只要在对应struct里加一个设备字段即可。3.2 残差函数与雅可比矩阵组装统一残差函数是主迭代的核心。我单独写了一个文件residual_full.m输入是所有状态变量拼接的大向量 x 和网络数据输出是三网残差的拼接向量 F。这个文件保持纯函数形式不改变任何外部状态方便在牛顿迭代中反复调用。雅可比矩阵的组装有两类思路数值差分和解析推导。初始版本我用Matlab自带的jacobianest之类的差分工具但多次评估后发现多能流系统的雅可比矩阵差分误差在耦合接口位置被放大导致牛顿法修正方向偏差。后来改为解析雅可比电网子块沿用传统潮流的 H、N、J、L 分块气网子块直接对Weymouth方程求导热网子块内部不对外层提供完整雅可比而是通过隐函数求导得到热出力对接口变量的灵敏度。对于气网Weymouth方程的雅可比我在零压差附近做了小邻域线性化处理这是保证牛顿法稳定过零点的关键function [f, df] weymouth(pi_i, pi_j, C, eps_scale) dpi pi_i - pi_j; if abs(dpi) eps_scale f C * dpi / sqrt(eps_scale); df C / sqrt(eps_scale); else f C * sign(dpi) * sqrt(abs(dpi)); df C / (2 * sqrt(abs(dpi))); end end这里 eps_scale 取 1e-6。实际测试中加入这个小邻域处理后管道两端压差过零不再引起牛顿迭代振荡气网和电网统一求解的第一步修正量就稳定得多。3.3 阻尼牛顿步与收敛判据多能流系统非线性极强直接全步长牛顿法很容易过冲。我采用简单的阻尼策略先计算完整牛顿步 Δx如果 Δx 的无穷范数超过0.5就按比例缩放dx -J \ F; if norm(dx, inf) 0.5 dx dx * (0.5 / norm(dx, inf)); end收敛判据采用残差二范数和接口变量变化量的双重判断。残差二范数要求小于1e-6接口变量CHP电出力、热出力、P2G功率的前后步变化量小于1e-5。单独只判断F范数不够因为在强耦合情况下F范数下降缓慢但接口变量其实已经震荡必须同时监控。主循环简版如下x init_state(pn, pg, ph); for iter 1:30 [F, J] build_full_residual_jacobian(x, pn, pg, ph, device); dx -J \ F; if norm(dx, inf) 0.5 dx dx * (0.5 / norm(dx, inf)); end x x dx; if norm(F, 2) 1e-6 break; end end3.4 热网水力-热力顺序求解片段热网子系统的实现独立封装成一个函数。水力方程用牛顿法解流量热力方程用拓扑前推回代解温度。前推回代是热网计算的高效做法具体是利用管道有向图对节点进行拓扑排序从热源节点开始逐个计算管道末端温度。代码片段如下T_s zeros(n_heat, 1); T_s(heat_source_idx) T_supply_set; % 热源供水温度给定 for k 1:n_heat_pipe head pred(k); % 当前管道首端节点 tail succ(k); % 当前管道末端节点 m_flow m_pipe(k); % 该管道流量来自水力求解结果 T_s(tail) T_env (T_s(head) - T_env) * exp(-lambda * L / (Cp * m_flow)); end需要非常注意的是如果管道流量 m_flow 低于某个下限指数项的分母接近零温降计算会失稳。我在代码里对流量设了下限0.01 kg/s低于该值的管道直接按绝热处理避免温度出现奇异值。这个处理虽然看似不严谨但实际工况下流量很低意味着该管道几乎不参与供热对整体能流分布的影响可以忽略。4. 算例实测33节点配电网耦合小型气网与热网跑一遍模型搭完了没有算例验证一切都是纸上谈兵。我搭了一个中等规模的测试系统把整个求解框架完整跑了一遍。4.1 算例拓扑与参数设置测试系统由三部分组成电网采用IEEE 33节点标准配电网基准功率10MVA基准电压12.66kV。在节点18接入一台电锅炉额定功率0.5MW在节点22接入一台CHP机组额定电出力0.5MW在节点25接入一台P2G装置额定功率0.3MW。电网总负荷保持不变。气网采用自定义的7节点网络节点0为天然气源给定压力基准2MPa标幺值1.0其余6个节点为负荷节点。CHP接在气网节点3燃气锅炉接在气网节点5P2G接在气网节点6。管道Weymouth系数经过标幺换算后在0.2到0.5之间总气负荷约1.2MW折算为天然气热值。热网采用6节点环状供热网络节点1为热源由CHP和燃气锅炉共同供热节点2和节点3为热负荷节点设计供热负荷总量1.0MW。设计供水温度90°C回水温度50°C环境温度20°C。4.2 迭代过程与结果分布先看收敛情况。基准工况下统一牛顿法冷启动的迭代过程如下迭代步残差二范数接口变量最大变化量14.7e-18.2e-236.2e-23.5e-353.1e-35.2e-484.5e-56.3e-6106.1e-78.9e-8前两步残差下降不快主要是气网Weymouth方程在初始压力接近时导数较大修正步被限幅。过了第三个迭代步后进入一致的快速下降区间整体收敛表现比较理想。整个计算时间在Matlab原生环境下约1.8秒主要是热网水力牛顿法的内部迭代占了一部分。从能流结果看基准工况下CHP以500kW电出力和700kW热出力运行气网气源总流量折算热功率约为1.75MW其中CHP耗气占0.95MW燃气锅炉占0.55MW其余为纯气负荷。电锅炉消耗电网功率约480kW含损耗向热网注入约450kW热功率。主网根节点从系统吸收的有功比不接入耦合设备时下降了约350kW说明本地CHP确实承担了一部分负荷减少了向上游电网的取电。4.3 耦合系数对能流分布的影响分析能流计算的价值不只是得到一组稳态解更重要的是观察关键参数变化后系统潮流分布的迁移方向。我做了两组灵敏度分析。第一组是改变CHP的热电比 α从1.4降到1.2。此时CHP在输出相同总燃料热量的情况下电出力增加、热出力减少。电网侧节点22注入的有功增大主网从根节点吸收的功率进一步下降约180kW热网侧CHP热出力不足以满足原热负荷燃气锅炉需要多烧天然气补足约140kW热出力导致气网气源总流量上升约8%。这个现象说明多能耦合系统中某一侧出力调整的影响会通过耦合设备传导到另外两个网络且方向往往是相反的——电少了气就多热少了气也多。第二组是P2G功率从300kW提升到600kW。此时电网在节点25增加了300kW电负荷相当于抬高了该区域的局部功率需求配电网末端电压降落幅度增大气网则多了约600kW折算热功率的天然气注入气源供给量反而下降。实际观察到的结果是主网根节点有功吸收量比基准工况增加约250kW气网节点压力整体抬升0.05标幺值。这和P2G“消纳电能、替代天然气供应”的运行逻辑完全吻合。4.4 收敛行为与计算效率实测我还专门对比了统一求解法和交替迭代法在同一算例上的表现。交替迭代法每个网络独立求解但接口变量更新采用直接替换方法收敛迭代次数残差终值计算时间(秒)统一牛顿法106.1e-71.8交替迭代法218.3e-72.6交替Aitken加速137.4e-72.0直接交替迭代法在CHP热电比为1.4时虽能收敛但接口变量在CHP电出力和热出力之间来回拉扯中间几步出现了小幅震荡。加上Aitken加速后收敛速度几乎追平统一求解法。如果未来系统规模扩大到几十个耦合设备我建议优先考虑统一求解法因为交替迭代法的接口变量数量增加后收敛半径会迅速缩小。5. 调试过程中踩过的坑与工程化建议这个项目从模型搭建到最终稳定跑通我最少调试了三天。很多问题在单网络潮流计算中根本不会出现但一旦多网耦合就会暴露出来。下面这几个坑值得每一个复现者提前避雷。5.1 气网Weymouth方程在零压差处的数值跳变最开始我直接使用 f C * sqrt(p_i² - p_j²)没有考虑流动方向。当收敛迭代中某个管道两端压力变得非常接近时sqrt 函数的导数趋向无穷大雅可比矩阵元素异常迭代直接从残差1e-3跳变到1e3。排查了好久才发现是气网分支的问题。解决办法就是我前面展示的小邻域线性化用 eps_scale 把零压差附近的流量-压差关系改为线性函数。这个修改不影响有效工况下的精度但能让牛顿法安全通过零压差区域。另一处细节是管道流量方向判定。用带符号的Weymouth形式时如果 dpi 和 m 的符号不匹配整个气网平衡方程的方向就不对残差永远降不下来。我建议在初始化阶段检查每根管道的连接方向是否与真实高低压关系一致不一致就交换首末节点。5.2 热网温度方程与流量迭代的耦合震荡热网水力计算给出流量分布热力前推回代再用这个流量算温度。问题在于温度变化会改变供回水温差进而改变节点热功率方程的等效质量流量需求反过来又影响水力方程中的节点注入流量。我在热负荷较大的节点上遇到过流量与温度交替迭代震荡的情况表现为供水温度在89到91°C之间来回跳始终不收敛。解决方法是给温度更新加一个低通滤波每次热力计算得到新的供水温度后不直接替换而是用 T_new 0.5 * T_new 0.5 * T_old 更新。这个阻尼系数看似简单却把热网子系统的收敛稳定性显著提升了。后来我又把水力方程中的回路压降方程残差阈值单独放宽到1e-5避免水力迭代过于苛刻导致温度更新时断时续。5.3 对统一雅可比矩阵的病态征兆提前预警统一雅可比矩阵规模大调试时很难直观看出哪里出了问题。我发现一个非常有效的信号观察行列式的数量级变化。正常迭代中雅可比矩阵行列式应该在同一数量级附近小幅变化如果某次迭代行列式突然跳变几个数量级说明矩阵中某个方程或变量连接出了问题。这个判断方法帮我在20分钟内定位到一个误写的管道关联矩阵行避免了靠打印逐行检查的海量时间。另一个实用技巧是分网络监控残差分量。每次迭代时把电网残差、气网残差、热网残差分别计算并打印出来观察是哪个网络先发散。F build_residual(x, pn, pg, ph); r_e norm(F(1:n_e), 2); r_g norm(F(n_e1:n_en_g), 2); r_h norm(F(n_en_g1:end), 2); fprintf(iter %d: 电网 %.3e, 气网 %.3e, 热网 %.3e\n, iter, r_e, r_g, r_h);三网残差分量几乎总是能告诉你发散源头在哪而不是让你在几百个方程里大海捞针。5.4 代码加速稀疏化、预分解与并行初值扫描工程落地的最后一步是性能优化。Matlab对稀疏矩阵的支持非常好我把节点导纳矩阵、关联矩阵、统一雅可比矩阵全部转换成 sparse 类型后计算耗时降了约60%。统一牛顿法中每次迭代都需要求解 J * dx -F这一步用直接分解需要消耗大量时间。实测表明如果在迭代前对雅可比矩阵做一次 LU 分解迭代内用分解好的 L、U 回代求解10步迭代可以省下约35%的时间。一开始我图省事用inv(J)后来发现不仅慢而且矩阵规模稍大时数值误差明显修正方向都会偏。建议大家无论如何都要用J \ F或者显式分解方式永远不要用逆矩阵。并行初值扫描是我后加的一个小功能由于不同工况下的初值敏感度不同我用parfor对多个初值做粗筛选择残差最小的初值作为正式迭代起点。这个方法在强耦合重负荷工况下特别有用比如当CHP热电比调到1.0以下的工况从标准初值出发可能发散但从几个随机扰动初值中选一个往往能顺利收敛。最后再分享一点个人体会这种多能流Matlab工程的调试周期通常比预期长一倍不是因为模型有多难而是因为三套网络的错误会互相掩盖。电网解不出来可能是气网参数错了热网温度发散可能是电网潮流修正量限幅太紧。我的建议是先断开所有耦合设备分别验证三个独立网络的能流求解器然后再逐个接回耦合设备。每接一个设备就重跑一遍全系统确认收敛后才继续下一个。整个过程虽然慢但每一步的失败原因都清晰可控。等你把CHP、电锅炉、P2G、燃气锅炉全部接入并稳定收敛之后再回头看这个“电气热能流计算”的框架就真的彻底掌握了。