去年做某个区域综合能源站的能流计算项目时我最开始的想法很简单电、气、热各自调用一套成熟求解器先把电网潮流算出来把电锅炉功率塞给热网把CHP机组耗气量塞给气网再反复迭代几轮就完事了。这个思路看起来顺理成章实际跑起来就原形毕露——当CHP容量占区域负荷比例超过30%时这种交替迭代方式经常出现振荡甚至发散程序连续跑十几轮都不收敛调试成本极高。后来我彻底切换思路把所有网络的方程统一拼成一个非线性方程组用一套加权的牛顿-拉夫逊求解器同时求解所有变量收敛性和稳定性一下子就上来了。这篇文章就围绕这套“计及多能耦合的区域综合能源系统电气热能流计算”完整建模和Matlab实现来写。我会把能源集线器的耦合矩阵怎么建、统一求解器怎么组装雅可比矩阵、收敛性怎么调试这几个核心问题一次性讲透。内容比较适合正在做综合能源系统仿真、需要自己编写多能流求解程序的研究生和工程师也适合刚接触综合能源系统、想快速理解多能耦合本质的读者。代码思路我已经在Matlab R2022b下完整跑通后面贴出的关键片段可以直接抄。1. 多能耦合能流计算的整体设计思路1.1 为什么传统“分网络依次算”的思路走不通单看电网潮流我们处理的是纯电气量负荷和发电机出力都是“给定的边界条件”。但在综合能源系统里电网的边界条件本身就是内部变量CHP机组消耗天然气来发电它发多少电不是由调度指令独立决定而是由气网能送来多少气、热网需要多少热共同决定。电锅炉消耗电能产热让电网负荷和热网热源互相锁定。燃气锅炉则直接建立气网流量和热网热源之间的映射。三个网络之间形成一条紧密咬合的链式关系。如果坚持“先算电网、再算热网、再算气网”的交替求解每一轮迭代时都必须强制假设另外两个网络的状态不变。这个假设在耦合很弱时还能凑合一旦CHP容量占比升高或者热泵、电锅炉这类电热耦合设备功率较大时网络之间的反馈增益就会被放大交替迭代极易出现数值振荡。我实测在强耦合场景下交替法经常需要上百次迭代也未必收敛而统一求解在同样场景下20次迭代左右就能稳定落到1e-6的精度。因此这个项目的第一设计决策就是把电网、气网、热网的所有方程全部放在同一个非线性方程组里用统一迭代同时求解所有变量。虽然雅可比矩阵规模变大、结构变得复杂但换来的是更强的收敛特性和更清晰的物理一致性。1.2 用能源集线器把设备耦合“翻译”成矩阵把三个网络的方程硬拼在一起还不够必须把耦合设备的关系用一种统一的形式表达出来。这里我用的是能源集线器Energy Hub理论。这个概念理解起来其实很容易它就像一个多端口网络黑箱输入侧接电、气、热等不同能源输出侧也是用户需要的各种能源形式黑箱内部的设备用效率矩阵来描述能量转换关系。举个例子一个典型的区域能源集线器内装有变压器、CHP机组、燃气锅炉和电锅炉。输入侧是电网购电P_el和网购气P_g输出侧是电负荷L_e和热负荷L_h。设备效率参数变压器效率η_T0.95CHP发电效率η_chp_e0.40CHP热效率η_chp_h0.45燃气锅炉效率η_gb0.85电锅炉效率η_eb0.98。那么能量平衡关系可以写成L_e η_T × (P_el - P_eb) η_chp_e × F_chpL_h η_chp_h × F_chp η_gb × F_gb η_eb × P_eb其中P_eb是电锅炉消耗的电功率F_chp和F_gb分别是CHP和燃气锅炉消耗的天然气量。再加上内部能量守恒F_chp F_gb P_g以及电母线平衡P_el P_eb P_grid_in就构成一个完整的Hub内部方程系统。这个抽象最大的价值是让设备级的不统一物理量在Hub层面收敛为一个统一的输入输出关系。后续不管是做能流计算、优化调度还是灵敏度分析都可以把整个Hub当作一个带内部状态的模块来使用不需要每个耦合设备单独写一套与网络交互的逻辑代码结构清晰得多。1.3 统一求解与交替求解一个关键取舍做多能流求解业界大致有两种路线。一种是交替迭代法Sequential Method按照电网→热网→气网的顺序逐个求解每个子网络可以使用自己成熟的求解器代码复用率高。另一种是统一求解法Simultaneous Method把所有平衡方程组成一个整体方程组一次建立统一的雅可比矩阵同时更新所有状态变量。我在项目中对这两种路线做了系统对比。交替迭代在弱耦合下表现不错但存在两个硬伤第一各网络子迭代之间的收敛判定难以统一你很难判断“电网已经收敛但气网还没进入状态此时该不该停止”第二强耦合下振荡严重往往需要设置很强的阻尼或引入松弛因子而这又会显著拖慢收敛速度。统一求解法虽然雅可比阶数高但现在Matlab里用sparse函数构建稀疏矩阵求解线性方程组Axb的成本并不高对于几十到几百节点规模的区域级综合能源系统完全可接受。最终我选择了统一求解法作为主求解器交替法只作为初值生成和结果校验的辅助工具。这是整个项目最核心的技术选型决策直接决定了后面的代码架构。2. 电气热能流建模的核心细节2.1 电网经典牛顿-拉夫逊但注意注入项来源电网子模型采用经典极坐标下的牛顿-拉夫逊潮流方程。对每个PQ节点和PV节点分别写有功和无功不平衡方程ΔP_i P_i^spec - V_i × ΣV_j × (G_ij×cosθ_ij B_ij×sinθ_ij) 0ΔQ_i Q_i^spec - V_i × ΣV_j × (G_ij×sinθ_ij - B_ij×cosθ_ij) 0这里的P_i^spec和Q_i^spec看起来和普通潮流一致但有一个关键区别注入功率里包含了耦合设备注入的功率项。比如CHP机组接入电网节点kP_k^spec中不仅有常规的发电机出力和负荷还要加上CHP电出力P_chp。而P_chp并非固定参数它是气网供气量F_chp的函数。于是电网节点方程里就出现了气网节点压力的“影子变量”这就为统一雅可比矩阵的非对角块引入了非零元素。我在实际建模时还额外考虑了一种特殊情况如果某个节点同时接了CHP和电锅炉那么该节点的净负荷可能变成负值即向电网倒送功率。这本身物理上没问题但会让牛顿潮流方程中的初值设定更加敏感后面我会说初值怎么处理。2.2 气网从Weymouth方程到节点流量平衡天然气网络采用稳态潮流模型管道流量与两端压力差的关系用经典的Weymouth方程描述f_ij sign(p_i - p_j) × C_ij × sqrt(|p_i² - p_j²|)其中p_i、p_j是节点绝对压力C_ij是管道传输系数包含管径、长度、气体组分、温度等参数影响。这个方程强非线性当压力差趋于零时流量对压力的导数趋向无穷大是气网收敛困难的主要来源之一。气网节点流量平衡方程为Δf_i Σ f_ij f_source_i - f_load_i 0其中f_load_i是节点取气负荷。注意这里的取气负荷不仅包括居民和工业用户还包括连接在气网节点上的燃气锅炉和CHP的耗气量。耦合设备的耗气量又要回到能源集线器内部方程来确定不能简单当作常数负荷。气网求解的状态变量一般取节点压力p_i参考节点气源节点压力已知作为平衡节点。实际项目中我还用到了压缩机模型但只在高压输气场站需要考虑区域级气网通常可以忽略——不过在做代码架构时要留好接口别把压缩机完全写死。2.3 热网水力与热力两个层面的方程组热力网络比电网和气网都复杂因为要同时描述水力工况和热力工况。水力工况的核心关系是节点流量连续Σ m_out Σ m_in以及回路压力平衡泵扬程等于沿线压降之和Σ Δp_loop 0模型的状态变量是各管道流量m和节点压力。热力工况则要区分供水管网和回水管网核心关系是节点功率平衡和管道温降节点混合温度由流入该节点的各支路流量加权平均得到T_mix (Σ m_in×T_in) / (Σ m_in)管道沿程温降按指数规律衰减T_end (T_start - T_amb) × exp(-λ×L / (c_p×m)) T_amb热源或热负荷节点向管网注入或取走的功率满足φ c_p × m × (T_supply - T_return)热网求解的状态变量我取了四个管道流量向量、节点供水温度向量、节点回水温度向量以及热源/换热站的供回水压差相关变量。这样做虽然增大了方程组的规模但物理上更完备能同时支持供热和供冷场景。2.4 耦合设备CHP、燃气锅炉、电锅炉、热泵怎么进模型这是整个项目最容易写乱的部分。我把耦合设备分成三类来处理第一类是燃气-电-热三向耦合的CHP第二类是燃气-热双向耦合的燃气锅炉第三类是电-热双向耦合的电锅炉和热泵。CHP机组我采用简化的背压式模型P_chp η_chp_e × F_chpH_chp η_chp_h × F_chp也就是说电出力和热出力都正比于燃料耗量热电比ρ η_chp_h/η_chp_e是固定值。这个模型对大多数区域级的规划研究已经足够。如果你研究的系统是抽凝式CHP热出力可以在一定范围内连续调节那需要在这个基础上加一个热电比可调变量代码架构里预留好接口就行。燃气锅炉和电锅炉相对简单分别是H_gb η_gb × F_gbH_eb η_eb × P_eb热泵则用COP模型H_hp COP × P_hp每一类设备接入Hub时设备耗电/耗气量、发电/发热量要同时反映到对应子网络的节点注入方程里。这是耦合设备建模的核心同一时间一台设备至少要同时修改两个网络的节点注入项并且这两个修改的变量在迭代中是关联的。我在代码里用结构体数组来存每一台设备的类型、接入节点、额定容量和效率参数这样后续无论增删设备都不需要改动求解器本体。3. 基于Matlab的统一能流求解器实现3.1 代码架构把三类网络映射为统一数据结构这套Matlab代码整体遵循“数据驱动”的组织方式。我不建议为电网、气网、热网分别设计完全不同的数据结构因为这会让统一求解器里到处是类型判断。更实用的是把三类网络统一成一个广义网络拓扑描述node所有节点的基态信息附加net_type字段标记该节点属于电网、气网还是热网branch所有支路的基态信息附加branch_type字段标记是线路、管道还是供热管道hub能源集线器信息包含设备清单和耦合矩阵状态变量则统一拼成一个列向量X [theta; % 电网电压相角(Vθ节点除外) V; % 电网电压幅值 p_gas; % 气网节点压力 m_flow; % 热网管道流量 T_supply; % 热网供水节点温度 T_return]; % 热网回水节点温度主程序流程非常清晰读入基础数据文件初始化状态变量进入统一迭代循环迭代收敛后写结果到结构体。function [X, iter, err_his] multiEnergyPowerFlow(data) X initState(data); for iter 1:data.maxIter [F, J] residualAndJacobian(X, data); dX -J \ F; % 带阻尼的线搜索 alpha 1.0; [F_new, ~] residualAndJacobian(X alpha*dX, data); while norm(F_new) norm(F) alpha alpha * 0.5; if alpha 1e-4, break; end [F_new, ~] residualAndJacobian(X alpha*dX, data); end X X alpha * dX; err norm(dX, inf); err_his(iter) err; if err data.tol, break; end end end3.2 统一雅可比矩阵的分块组装统一雅可比矩阵是这个程序里最需要注意的部分。它的结构可以写成下面这个分块形式J [J_EE, J_EG, J_EH J_GE, J_GG, J_GH J_HE, J_HG, J_HH]其中对角块J_EE、J_GG、J_HH分别对应电网、气网、热网自身变量之间的偏导数关系这部分与单独的网络潮流计算基本一致。非对角块J_EG、J_GE、J_EH等才是多能耦合的真正体现——它们记录了某一个网络的变量对另一个网络平衡方程的灵敏度影响。以CHP机组为例它接入电网节点k和气网节点n电功率平衡方程的形式为F_e,k P_k^spec - V_k×ΣV_j×(G_kj×cosθ_kj B_kj×sinθ_kj) - P_chp(X) 0其中P_chp是X中气网变量p_gas,n和Hub内部变量的函数所以J_EG(k, n) ∂F_e,k / ∂p_gas,n -∂P_chp / ∂p_gas,n ≠ 0这就是一个典型的跨网络雅可比非对角元素。组装时用稀疏矩阵最舒服% 组装统一雅可比矩阵 J sparse(nVar, nVar); % 电网对角块 J(ie, ie) dFe_dXe; % 电网自身偏导数 % 气网对角块 J(ig, ig) dFg_dXg; % 气网自身偏导数 % 热网对角块 J(ih, ih) dFh_dXh; % 热网自身偏导数 % 耦合交叉块 J(ie, ig) dFe_dXg; % 气网变量对电网方程的影响 J(ig, ie) dFg_dXe; % 电网变量对气网方程的影响 % 其他交叉块依此类推初学的时候最容易漏掉交叉块。我曾经犯过一个错吉川矩阵只组装了对角块交叉块全为零迭代后电气热三个网络的变量各自收敛但整体能量不平衡——因为CHP的发电量、耗气量和发热量在三个子网络里对应的不是同一条曲线程序“看起来正常”实际结果完全物理失真。检查雅可比交叉块是否完整的一个简单方法算完一次迭代后用有限差分法数值验证每一个非零元素的位置和值比对解析结果。3.3 能量集线器内部方程装载方式Hub内部方程并不是单独求解的而是作为整体残差向量的一部分参与统一迭代。我把每个Hub的能量平衡方程嵌入到统一的残差向量F末尾让Hub内部变量也成为X的一部分% 残差向量组装 F [F_electrical; % 电网节点平衡方程 F_gas; % 气网节点平衡方程 F_hydraulic; % 热网水力方程 F_thermal; % 热网热力方程 F_hub]; % 能源集线器内部平衡方程同时状态向量扩充为X [theta; V; p_gas; m_flow; T_supply; T_return; F_chp; % 各Hub的CHP耗气量 F_gb; % 各Hub的燃气锅炉耗气量 P_eb; % 各Hub的电锅炉电功率 P_hp]; % 各Hub的热泵电功率Hub内部设备的耗气量、耗电量都不再“打包成负荷”而是显式参与迭代。这样做的好处是在牛顿法中设备变量与网络变量同时更新不存在时间差和数值错位。代价是雅可比矩阵规模进一步扩大但对于几十个Hub规模的系统来说稀疏矩阵求解完全不是瓶颈。3.4 收敛判据与阻尼控制收敛判据我采用了混合策略同时检查变量增量范数和方程残差范数。纯电网潮流一般只看一个量但多能流里不同物理量的量纲差异很大单纯的电压误差无法反映气网压力和温度的收敛状态。我的设置是data.tol 1e-6; % 最大变量增量标幺化后 data.tolF 1e-6; % 最大方程残差所有变量进入统一向量前都做了标幺化处理。电网变量用系统基准容量100MVA标幺气网压力用基准压力50bar标幺温度用基准温差50℃标幺。这一步非常重要否则雅可比矩阵中温度相关元素和电压相关元素的数量级可能差出四个量级导致矩阵病态。阻尼控制使用最简单的线搜索策略。每次得到牛顿方向dX后先尝试全步长alpha1如果残差范数不降反升就把步长减半最多减4次。这个简单策略在实际中非常有效能解决绝大多数收敛振荡问题。我更推荐在正式计算前先打印出每次迭代残差范数变化看一眼就知道是初值问题还是方程本身病态问题再决定是调整初值还是调整阻尼。4. 收敛性调试与常见问题排雷4.1 量纲悬殊导致的病态雅可比第一次跑完全耦合算例时我发现雅可比矩阵的条件数达到了1e14几乎接近数值奇异收敛速度也异常慢。排查后发现典型的问题气网方程里压力和流量是几十bar量级的平方关系而电网方程里电压只有幅值0.95~1.05pu两者之间的交叉偏导数数量级差距巨大。如果不做任何缩放极坐标牛顿法在交叉迭代时会把气压的小扰动放大到电压方程里造成虚假的电压波动。解决思路很直接所有非标幺变量进入统一方程组前都做基准标幺缩放不是只缩放变量本身还要同步缩放对应的方程。具体操作时我会为每个方程设置“基准残差”将电功率方程残差除以100MVA将气网流量方程残差除以基准流量将热网功率方程除以基准热功率这样统一雅可比矩阵中各类元素的数量级能控制在1e-3到1e3之间条件数可降到1e8左右牛顿法收敛就顺畅了。4.2 初值选择和耦合环节的收敛抖动多能流对初值的要求比纯电网潮流高得多。我第一次跑综合算例时电压初值给了统一的1.0pu结果电网部分飞快收敛但气网节点的压力初值给得太低导致CHP的耗气量成为负数热量平衡方程直接物理越界迭代一下子就崩了。之后我总结了一套比较稳妥的初值策略电网变量用平启动电压幅值1.0相角0气网压力初值设为气源压力的85%~90%线性递减分布这样Weymouth方程的压力差初值不会太小避免了流量计算初始导数过大热网管道流量初值统一按额定流量的0.7倍给定供水温度给设计工况温度回水温度给设计工况温度减10℃。Hub的设备变量比如F_chp用“该Hub热负荷设计值除以对应设备效率”来做初值。这套策略在三个不同规模的测试系统上都很稳基本都能在20次迭代以内达到1e-6精度。4.3 强耦合场景下统一迭代失稳的处理当系统中CHP容量占比超过一定程度或者热泵装机很大时我遇到过统一迭代在某个中间步突然发散的情况。排查发现问题出在“某个物理量越过了不可行域”。比如温度初值设置不合理某条供热管道的末端温度可能低于环境温度此时指数温降公式会出现绝对值异常增长残差范数急剧变大。处理方法有两个可以叠加使用。第一是变量限幅每一步牛顿更新后检查物理量是否在合理范围内比如节点温度必须大于环境温度、气网压力必须保持正定如果越界就按边界值截断。第二是自适应阻尼增强把线搜索从“残差不降则减半”升级为“残差下降率小于5%就继续减半”虽然略微增加迭代次数但大幅提高稳定性。我最终工程中用的是阻尼上限0.8、下限0.125、每步搜索次数上限5次的方案所有测试工况都能平稳收敛。4.4 收敛性对比数据实录为了验证统一求解器的实际性能我搭了一套约60节点规模的测试系统电网含33个节点气网含11个节点热网含14个节点和2个能源集线器。一台上是CHP容量占比35%另一台上是电锅炉供暖占比20%。用同一套数据分别跑交替法和统一法求解方式初始场景迭代次数收敛精度是否出现振荡交替迭代681e-5是多次振荡统一迭代无阻尼261e-6否但中间有波动统一迭代带线搜索181e-6否平滑收敛5. 常见问题速查与调试建议下面是我在整个项目调试过程中整理出的问题速查表特别适合初学综合能源能流计算的读者出现问题优先对照排查。现象可能原因排查方向处理建议电网部分正常气网压力越变越离谱气网初值偏离平衡点太远检查Weymouth方程初始压力差是否过小用气源压力85%~90%作初值热网温度出现负值或超过100℃温度初值不合理或迭代步长过大检查管道温降公式是否被大流量冲过头加入变量限幅限制温度物理边界非对角雅可比为零但结果异常交叉块漏组装用有限差分验证每个交叉偏导数补全J_EG/J_GE等交叉块迭代到某个点残差范数不降反升变量越过物理不可行域打印中间物理量检查CHP耗气量是否为正使用阻尼线搜索并限幅收敛慢每次迭代只降一个量级各网络变量数量级差异大检查雅可比条件数和各元素数量级做标幺化和残差缩放收敛后总能量不平衡Hub内部方程未完全耦合进统一残差单独算每个Hub内部能量进出是否守恒把设备变量显式加入X而不是打包成负荷大规模系统迭代速度慢雅可比未用稀疏矩阵存储检查是否用了full而非sparse用sparse组装线性求解用反斜杠运算符调试时我最推荐的习惯是写一个数值差分校验函数每次修改模型后都跑一遍验证解析雅可比矩阵和数值雅可比矩阵是否一致。这个习惯救过我很多次尤其是新增一种耦合设备类型时交叉块里很容易漏掉一项偏导数。关于初值判断有个更工程化的技巧先按纯电网牛顿潮流跑一遍把网络状态基本摸清固定电网状态单独迭代一次气网和热网再把三个网络的中间状态作为统一迭代的初值。这个过程看起来多花几毫秒但能避免大量“初值规律不好”的返工。还有个针对热网的细节。热网水力方程里的回路压降约束在不同环路数量较多时容易出现欠约束的情况。我的经验是建模阶段就要明确哪个管道流量是独立的、哪些通过节点流量连续性自然确定。如果回路方程写多了一行或写漏了一行雅可比矩阵会奇异性增强收敛阈值很难达到1e-6。遇到这种情况检查热网支路的流量变量是否和节点压力变量一一对应即可。关于阻尼系数我试过固定0.9的效果没有自适应阻尼好试过纯Newton-CG方法共轭梯度法在小规模上可行但代码复杂度偏高最终稳定采用“阻尼线搜索”的经典组合。对区域级多能流问题这套组合花费最少、收益最直接。6. 一点个人项目的体会与扩展思考最后说一点个人经验。这套代码我在自己的项目里维护了快一年最大的感悟是综合能源能流计算的难点从来不在数学推导而在于工程数据的完备性。热网里一个换热站的实测供水温度、气网里一个调压站的出口压力、CHP机组的实际热电比曲线这些数据只要有一项不准确模型再精巧也白搭。如果你只是需要一个快速可用的落地工具建议把精力花在数据清洗和参数辨识上而不是一味追求更花哨的求解算法。另外有个经常被忽视的细节电网基准容量和热网基准功率的选取如果差太多会直接影响标幺化后的雅可比矩阵条件数。我在项目中试过电网基准100MVA、热网基准50MW、气网基准100kW的组合效果就明显比三者用同一个基准要好。这个需要根据你的具体系统规模去调没有一个固定公式但调试思路是一样的——先打印一次雅可比矩阵的非零元素数量级分布看看哪一块明显偏离1e0~1e2针对性地做缩放。后续如果想扩展可以从稳态能流走向准动态时序仿真把光伏出力曲线、储能充放电策略和需求响应模块加进来。这个框架的扩展性也不错只要在Hub内部新增一个设备类型再补上对应的偏导数解析式统一求解器不需要做结构性改动。如果你正卡在多能流收敛或者Matlab代码实现的问题上可以对照这个项目里的思路重新捋一遍自己的建模逻辑多半能找出问题所在。