计及多能耦合的区域综合能源系统电气热能流计算研究Matlab代码实现做综合能源系统仿真这几年我接触最多的一个方向就是多能耦合系统的稳态能流计算。很多人一听“电气热能流”就以为是把电网潮流、天然气网水力计算、热网水力热力计算三个东西分别跑一遍然后拼在一起真上手做一遍就会发现事情没那么简单。一个区域综合能源系统里CHP机组同时连着电网和气网电锅炉同时连着电网和热网燃气锅炉同时连着气网和热网这些耦合设备把三条网络硬生生绑成了一个整体。你单独算电网潮流的时候气网供给CHP的燃料量是未知的单独算气网的时候CHP的耗气量又取决于电网的发电需求而热网的供水温度又反过来影响CHP的抽汽量。这篇内容就是基于我实际做过的一个Matlab仿真项目把整套电气热能流计算从数学模型、求解策略到代码实现完整梳理一遍。适合正在做综合能源系统方向的研究生、刚入行做园区能源规划的工程师以及想把多能流计算模型真正落地成可运行代码的人。我会把每一步为什么这么做、踩过哪些坑、参数怎么定都写清楚。1. 问题是先理清楚的为什么单算三条网不算真正意义的多能能流1.1 单能流系统已经有成熟工具为什么还要做联合能流计算先说个最直观的类比你手头有一个电网潮流程序、一个天然气网水力计算程序、一个热网水力热力计算程序三个程序各自都跑得很顺。但把它们拼在一起的时候问题就出来了。电网潮流需要知道气网供给CHP的天然气量从而确定CHP的发电功率和热功率气网需要知道CHP的耗气量才能做节点流量平衡热网需要知道CHP的供热量才能确定节点温度和管道流量。这三个量互为因果谁都没法先独立确定。这就是所谓的多能耦合的内在复杂性。严格来说只有当耦合设备两侧的变量都被当作未知数、放进同一个迭代框架里求解才能叫计及多能耦合的能流计算。如果只是把CHP当作电网里的一个PQ节点气网里当作一个固定负荷那本质上还是三个孤立的单能流问题谈不上耦合。另外单能流系统的求解工具已经非常成熟比如电网潮流有Matpower天然气网有管网仿真软件热网也有专门的热力计算工具。但找遍市面没有一个现成工具能一次性求解电、气、热三网联立的能流方程。原因其实也简单这三个网络的物理特性差异太大。电网是代数方程主导求解讲究收敛性和电压约束气网是压力驱动的非线性方程跟电网形式有相似之处但物理含义完全不同热网更特殊它被水力方程和热力方程双重绑定先要解出各管道的流量分布再在那个基础上求解温度的传播与混合。物理特性差异这么大就很难用一套通用的求解器统一处理。1.2 一个典型的区域综合能源系统拓扑长什么样在动手之前先把对象系统画清楚。我这次仿真的典型区域综合能源系统由三部分组成一个IEEE 9节点电网可以外购电也可以由CHP发电一个6节点天然气网气源经过加压站、管道供给CHP和燃气锅炉一个8节点热网包含若干供热管道、换热站和热负荷节点。这三块之间的连接点就是耦合设备的安装位置。CHP机组挂在电网3号节点、气网4号节点、热网5号节点上燃气锅炉挂在气网5号节点和热网7号节点上电锅炉挂在电网7号节点和热网6号节点上。这样三个网络通过耦合设备形成了一个闭环气网给CHP供气CHP发电并产热发电送入电网、热量送入热网电网再给电锅炉供电、电锅炉补充供热。这个拓扑结构比较典型涵盖了综合能源系统里最常见的三种耦合设备CHP气转电热、燃气锅炉气转热、电锅炉电转热。做代码实现的时候这几个设备的模型就得分别写清楚因为它们的工作特性完全不同。我后面的算例、代码结构也都是围绕这套拓扑展开的。2. 电气热三网联立求解的数学建模核心方程式逐一拆解2.1 电网潮流的极坐标形式与节点类型处理电网部分我采用的是最经典的新ton-Raphson极坐标潮流算法。对每个PQ节点有两个待求的误差量有功功率误差和无功功率误差对每个PV节点只有一个有功功率误差电压幅值给定。方程组的核心表达就是功率平衡方程P_i V_i * Σ V_j * (G_ij * cos(θ_ij) B_ij * sin(θ_ij)) Q_i V_i * Σ V_j * (G_ij * sin(θ_ij) — B_ij * cos(θ_ij))这里的G_ij和B_ij对应支路导纳的实部和虚部θ_ij是节点i和j的相角差。常规潮流里平衡节点承担全网功率差额典型做法是设一个外网等效节点给定1.04的电压幅值和0度相角。做综合能源系统能流计算时电网部分有一个必须注意的改动CHP机组不能简单当作固定的PQ节点或PV节点。因为CHP的发电量是气网供气量和自身热电比共同决定的结果在气网没有收敛之前你根本不知道CHP能发多少电。所以在联合求解框架里CHP的并网节点一般当作P节点处理有功功率由耦合迭代给出电压幅值由无功调节维持无功功率不做限制。这样处理的好处是电、气、热三个网络的变量可以被耦合方程统一调度坏处是雅可比矩阵的结构和形成方式要跟着改不能直接照搬标准潮流程序的节点分类逻辑。2.2 天然气网稳态能流的Weymouth方程与节点压力计算天然气网的稳态能流计算核心是管道的Weymouth方程。高压输气管道中气体的质量流量和管道两端压力的平方差成正比这一点和电流与电压差的关系有些类似但又不太一样。一般表达式写作F_ij sign(π_i — π_j) * k_ij * sqrt(|π_i² — π_j²|)其中π_i和π_j是节点i、j的压力k_ij是与管道内径、长度、压缩因子、气体性质有关的常数。节点流量平衡方程为Σ F_ij L_i_gas该式表示流入节点的天然气净流量等于该节点的负荷包括天然气负荷和耦合设备耗气量。气源节点是压力已知的边界节点其他节点的压力是待求变量。在Matlab实现里这个非线性方程组的求解也使用牛顿法。雅可比矩阵的构造逻辑是对每一条管道支路求F_ij对π_i、π_j的偏导数对每个节点累加所有关联支路的偏导量。因为Weymouth方程的导数分母带有sqrt项当节点压力差特别小或者接近零时导数计算容易数值溢出。我后面会专门讲这个问题这里是初学者最容易忽略的。2.3 热网水力-热力双重耦合的方程体系热网的建模是三个网络里最容易写乱的一个因为变量最多、量纲也最容易出问题。热网的能流方程分为两个层次水力方程描述管道流量和水泵节点的压力平衡关系。其首要约束是节点流量连续性流入某节点的体积流量等于流出该节点的流量与供热负荷的流量需求相平衡。其次是环路压力方程在闭合的供热环中沿回路的压力损失之和为零这一组方程将管道阻力与流量平方和管径相关和网络拓扑紧密耦合。热力方程则关注温度场每个节点的供水温度根据上游管道混合温度确定节点处的热量平衡方程可写为Φ_i C_p * m_i * (T_s,i — T_r,i)其中C_p是水的比热容m_i是流经该节点的工质流量T_s,i和T_r,i分别是该节点供回水温度。管道沿线的温度损失方程采用指数衰减模型T_end T_amb (T_start — T_amb) * exp(-λ * L / (C_p * m))这里的λ是管道单位长度的热损失系数L是管道长度。水力方程的结果即各管道流量分布是热力方程求解的前提。所以在求解热网时通常先解水力方程得到m的分布再解热力方程确定温度分布。两者合在一起才算完成了热网的能流计算。这个先水力、后热力的顺序不能颠倒否则流量未知热量方程根本无从下手。2.4 耦合设备建模CHP、燃气锅炉、电锅炉的参数化处理耦合设备是多能流计算的灵魂所在。我在代码里实现三种设备的模型分别对应不同数学形式CHP机组采用典型的背压式机组模型数学模型可以表达为发电功率P_e和供热功率H_heat都与燃料消耗量F_gas呈线性关系。即P_e η_e * F_gas * LHV_gasH_heat η_h * F_gas * LHV_gas。其中η_e是发电效率η_h是热回收效率LHV_gas是天然气低位热值。背压式CHP的热电比是固定的这一特性使得电功率和热功率之间呈严格比例关系。在联立求解过程中可根据系统的电、热负荷需求直接确定燃料消耗量进而实现电网、气网和热网之间的信息交换与约束统一。燃气锅炉模型是最简单的就是一个转换效率η_gb天然气输入F_gas热量输出H_heat η_gb * F_gas * LHV_gas。电锅炉模型类似只是输入变成了电能P_e输出热量H_heat η_eb * P_e。这三个设备在Matlab里分别写成三个函数输入输出都是结构体变量。这样迭代求解的结构就能统一耦合设备从电网取/送出的电功率影响电网潮流方程从气网取用的气功率影响气网流量平衡方程向热网注入的热功率影响热网温度计算。设备模型越真实迭代间的变量传递就越复杂但原则是清晰统一的每种耦合设备本质上都是三条网络之间的一条数据通道。3. 求解策略怎么选统一求解法与分解迭代法的取舍3.1 两种主流框架的对比与适用场景多能耦合能流的基本求解策略学术界和工程上无非两条路统一求解法和分解迭代法。统一求解法将电网、气网、热网的所有方程放在一起形成一个大规模非线性方程组一次性用牛顿法迭代求解。好处是收敛性好尤其是当耦合程度很强比如CHP占系统发电比例很高时不容易出现发散。缺点是雅可比矩阵规模大、结构复杂编程量大且气网、热网的方程特性差异大矩阵容易病态。更要命的问题在调试阶段一旦不收敛你很难判断问题出在电网方程、气网方程还是热网方程上。分解迭代法的思路则是三个网络各自独立求解耦合变量在外层做一个迭代循环。迭代流程是先给定耦合设备的初始运行状态比如假设CHP发电功率和热功率分别求解电网、气网、热网的能流然后根据三网的结算结果更新耦合设备的运行状态再重新求解三网。如此循环直到前后两次迭代的耦合变量偏差小于收敛阈值。我在实际项目里用的是分解迭代法改进版本。原因很直接编程调试灵活哪一层出问题一目了然而且每条子网的求解函数可以单独测试跟已知结果的单能流系统对比验证大大降低了出错概率。缺点是单纯分解迭代在强耦合场景下可能迭代震荡或者收敛极慢。针对这个问题我给耦合变量的更新加了松弛因子效果非常明显X_new X_old α * (X_calculated — X_old)这里α取0.5到0.8之间比较稳妥。α太大会导致震荡发散α太小则收敛速度慢。当CHP的出力占比超过系统总出力的30%时保守一点取α0.5占比较低时取0.8可以兼顾速度和稳定性。3.2 初值策略的选择直接影响能否收敛的关键一环初值设定是这个项目里我踩过最多坑的地方。电网潮流的初值相对简单平启动即可所有PQ节点电压幅值取1.0相角取0。但气网的初值极为敏感管道Weymouth方程里有sqrt(|π_i² — π_j²|)项如果给定初值使得节点间压力差过大或者节点压力过低迭代过程中平方差出现负值或者开方项数值异常就直接挂了。我的做法是根据气源压力和工作压力范围先预估一下整个气网的基准压力。假如气源节点压力是2.0MPa末端节点压力大概在1.2MPa附近那所有节点的初值统一取1.6MPa别去费心思逐个节点估算。这样做的原理是牛顿法的局部收敛特性要求初值不能偏离真解太远取平均压力能保证初始压差在合理范围内从而避免开方项出现非法值的概率。热网的初值也类似。供水温度取设计值的上端比如100°C回水温度取80°C所有节点统一取这两个基准值。先解水力方程得到流量分布再在流量基础上解热力方程温度场的迭代比较温和。判断收敛的标准我设为三网的最大不平衡量小于1e-6标幺值。注意是三个网络的不平衡量都同时满足不是只看某一个网络的收敛情况。因为耦合系统的特点就是各网络相互牵制只收敛其中一个没有意义。3.3 解耦迭代的收敛判据与松弛更新逻辑外层耦合迭代的流程我写成下面的伪代码逻辑初始化假设CHP的发电功率P_e、热功率H_heat燃气锅炉的供热量电锅炉的耗电功率电网潮流求解把CHP作为P节点有功给定电锅炉作为负荷节点求解电网模型得到各节点电压、相角和各耦合节点的结算结果气网能流求解根据CHP和燃气锅炉的耗气量由第一步的初始假设或上一次迭代结果折算作为气网负荷求解气网模型得到修正后的节点压力和可用供气量热网能流求解根据CHP、燃气锅炉、电锅炉的热功率注入与热负荷需求求解热网模型得到供热温度和流量根据气网结算得到的流量反推CHP的发电功率上限并修订热功率值再根据热网的实际供热需求修正CHP和电锅炉的分配判断修正量与上次迭代的偏差是否小于收敛阈值不满足则返回步骤2继续迭代所有偏差满足阈值后输出全局能流结果这套流程写出来看着简单实际操作中容易忽略一个点步骤5的耦合变量更新要带上松弛因子而且每个设备的更新松弛因子最好单独设置。比如CHP的更新比较稳定α可取0.8电锅炉因为受电网和热网双重约束振荡风险高α取0.5。4. Matlab代码实现全过程从函数架构到算例验证4.1 全套代码的文件结构与核心数据流设计下面是这次项目全套Matlab代码的文件组织方式每个文件的功能都做了拆分互不干扰system_main.m 主程序入口控制整体迭代流程 load_system_data.m 读入电、气、热三网的拓扑参数与负荷数据 build_ybus.m 根据支路参数构建电网节点导纳矩阵 electric_powerflow.m 电网潮流求解牛顿-拉夫逊法 gas_powerflow.m 天然气网能流求解牛顿法 hydraulic_solve.m 热网水力平衡方程求解 thermal_solve.m 热网热力方程求解 coupling_update.m 耦合设备变量更新与松弛处理 plot_reports.m 三网能流结果的可视化输出在设计数据流的时候我建议所有中间变量和结果都用结构体变量按网络分类封装。三个大结构体sys_elec、sys_gas、sys_heat分别存放各网络物理模型所需的数据比如节点参数、支路参数、求解结果。耦合设备的参数单独放在struct_coupling结构体里。这个架构的好处是每段代码的输入和输出都是清晰的数据结构函数之间的耦合度被降到最低。比如electric_powerflow函数只接收sys_elec结构体和耦合节点的注入功率只返回更新后的电压和相角结果它不需要关心气网压力是多少、热网供水温度是多少。模块解耦的好处在问题排查阶段优势尤其明显电网不收敛时只需要检查电网这部分的网络连接和参数配置就行不会被气网和热网干扰。4.2 核心求解函数关键技术点牛顿法雅可比矩阵的组装先看电网潮流求解的核心循环。标准牛顿法潮流的核心就是构造雅可比矩阵并反复求解增量方程for iter 1:iter_max % 计算功率不平衡量 [dP, dQ] compute_power_mismatch(V, theta, Ybus, S_spec); mismatch [dP; dQ]; if max(abs(mismatch)) tol break; end % 构造雅可比矩阵 J build_jacobian(V, theta, Ybus); % 求解增量 dx J \ mismatch; % 更新状态变量 theta theta dx(1:n_node); V V dx(n_node1:end) .* V; end这段代码看起来简单但有几个细节需要强调。雅可比矩阵的组装公式推导可以参考任何一本电力系统分析教材但在Matlab里有一个程序效率的大坑如果直接用两个for循环逐个元素填充雅可比矩阵三节点系统跑起来没问题一旦节点数超过30每一次迭代都会因为矩阵组装太慢让人等到怀疑人生。更好的做法是利用稀疏矩阵的索引批量赋值或者直接使用vectorized方式按节点对数组同时计算。我实测过同样的IEEE 9节点算例向量化后的代码比双重循环快至少一个数量级这个差距在蒙特卡洛批量计算或者考虑多场景规划时会体现得非常明显。再来看看气网的求解核心。气网求解的难点其实不在牛顿法的框架而在于Weymouth方程导数的正确处理for iter 1:iter_max % 计算节点流量不平衡量包括耦合设备负荷 node_mismatch compute_gas_mismatch(Pressure, gas_loads, pipe_k, pipe_connect); J build_gas_jacobian(Pressure, pipe_k, pipe_connect); deltaP J \ node_mismatch; Pressure Pressure deltaP; if max(abs(node_mismatch)) tol break; end end这里build_gas_jacobian函数有一个关键逻辑对于管道ij流量F_ij对节点压力P_i的偏导要分为两种情况计算。当P_i大于P_j时导数为正值为k_ij * P_i / sqrt(|P_i² — P_j²|)。当P_i小于P_j时导数为负。由于开方项位于分母中当压力差趋近于零时雅可比矩阵元素会趋向无穷大该数值特性极易导致迭代不收敛。为避免此问题需在程序中对压力差设置最小限值例如当|P_i² — P_j²|小于某个阈值我取1e-6时将该导数按上限值截断处理。曾经由于忽略该细节导致程序反复发散排查了很久才意识到问题所在。4.3 一个典型算例的测试结果与误差分析我用的算例配置如下CHP安装在三网交汇的关键节点电功率输出为60MW热功率输出为50MW燃气锅炉的热功率为30MW电锅炉的热功率为20MW电网侧的总电负荷约为110MW外网联络线提供剩余约30MW功率气源节点压力设定为2.0MPa两个气网负荷节点的天然气流量需求分别约为4000和2500 m³/h。经过分解迭代求解整个系统在14次外层迭代后收敛松弛因子取0.6各网络的不平衡量均达到10⁻⁶级别。电网的平衡节点有功出力为-28.5MW即从外网购电28.5MW各节点电压幅值均在0.98到1.03之间满足常规设计要求。气网的节点压力从气源2.0MPa下降到末端约1.35MPa压降符合管道长度和流量的物理预期。热网的供水温度经过管网衰减后最远负荷节点比热源出口下降了约9°C回水温度分布也比较均匀。为了验证程序的正确性我做了两组对照实验。第一组是单一电网潮流对照把CHP当作固定出力节点运行电网潮流程序与Matpower的计算结果对比节点电压偏差小于1e-6说明电网模块实现无误。第二组是气网对照把已知流量和压力参数的6节点气网用商业管网仿真软件校核结果压力偏差在2%以内。这样逐模块验证过的程序再去做三网联算才敢说结果是可信的。5. 实战中一定会遇到的坑三个典型问题与排查实录5.1 气网开方项异常造成不收敛的排查过程项目最早期调试阶段气网部分一直报出NaN程序直接崩溃。当时把错误定位到sqrt处但百思不得其解因为从物理概念上节点压力平方差为负似乎不该出现。后来通过打印每次迭代的节点压力值发现迭代第一次就出现某个节点压力值跳到接近零随后它与相邻节点间的压力平方差直接变负开方返回复数。根因其实是两个问题叠加。第一个是初值给得太差我把节点初值都设成0.5MPa但气源压力是2.0MPa迭代前期压力调整跨度大部分节点被牛顿法推到了不合理区域第二个是没有在Weymouth方程里加上数值保护当开方根内部出现极小负数时直接就得到NaN了。修正方案是两个方向同时入手初值统一取工作压力范围的中间值同时在Weymouth公式里加入限幅保护用max(eps, abs(p_i²—p_j²))来保险。5.2 耦合变量振荡不收敛为什么需要逐设备设置松弛因子框架刚跑起来的时候出现了比较有意思的现象电压和压力都稳定但CHP的出力在相邻两次迭代之间来回跳幅度还不小±8MW左右导致外层循环无法收敛。直觉解决方案是把全局松弛因子调小发现效果有限因为电网、气网模块本身是精确计算的真正振荡的来源是耦合变量更新逻辑——我先更新CHP发电量再根据它来更新耗气量而耗气量又马上被气网结算结果修正形成正反馈回路。后来我把CHP、燃气锅炉、电锅炉的松弛因子分别独立设置CHP为0.6燃气锅炉为0.5电锅炉为0.4。因为电锅炉同时受电网和热网两侧约束其更新最容易引起波动。结果外层迭代从完全不收敛变成14次收敛效果立竿见影。这个教训让我意识到在多能流迭代里不要贪图代码简单就用一个统一的松弛因子不同耦合设备的时间常数和灵敏性完全不同分开设置才是合理做法。5.3 单位换算的隐蔽炸弹天然气标方与热功率的换算最后一个要提醒的坑是单位。这是我检查看了半天才发现的低级错误但后果很严重。天然气网里气负荷的单位习惯用m³/h标方而电气热联立方程里CHP的耗气量需要换算成MW热功率才能和电网功率平衡方程式统一量纲。天然气低位热值LHV一般取36 MJ/m³换算关系是1 m³/h天然气约等于0.01 MW。计算方法是36 MJ/m³乘以1000除以3600秒结果正好约等于10 kW也就是0.01 MW。这个数值关系看起来很直白但如果代码里直接用一个常数去乘除很容易在某处漏掉或者写错一个数量级。更隐蔽的是锅炉效率、CHP效率这些百分比数在单位换算时也要参与运算一旦顺序弄错结果就可能出现CHP耗气量比气网总供气量还大的离谱局面。我的建议是在所有计算之前定义一组统一的基准值和转换函数把天然气质量流量、标准体积流量、热功率三套单位全部转换成统一的标幺值或者统一的MW制后再进入迭代计算。程序内部只用一套单位制只在输入和输出层做转换。6. 一些个人经验和可以继续扩展的方向做这个项目最大的体会是多能耦合能流计算耦合两个字才是重点。不要把它当成三个单能流程序的简单拼接而是要在求解框架设计时就充分考虑耦合变量的传递机制。初期我犯的错就是把80%的精力花在三个网络各自的算法实现上最后联调阶段才发现真正的难点在更新与收敛策略上。另外一个实用的经验是Matlab里做这类项目一定要充分利用结构体和稀疏矩阵。结构体让三网的参数管理变得清晰而稀疏矩阵让大规模节点的求解效率成倍提升。如果只是为了能跑通用全矩阵也凑合但如果你以后想把算例规模扩大或者在项目基础上做优化调度、多场景分析全矩阵的性能瓶颈就会非常明显。这套代码后续可以扩展的方向还挺多一是把稳态能流扩展为准稳态时序仿真模拟一天24小时或者四季的更替场景二是在能流计算基础上叠加经济优化目标比如以运行成本最小为目标函数调整CHP出力比三是加入储能设备模型让气网、热网具备调节能力看看储能对多能互补系统能流分布的影响。如果你已经在做或者准备做综合能源系统的仿真分析欢迎一起交流具体实现中的细节。
多能耦合下区域综合能源系统电气热能流联合计算与Matlab实现
计及多能耦合的区域综合能源系统电气热能流计算研究Matlab代码实现做综合能源系统仿真这几年我接触最多的一个方向就是多能耦合系统的稳态能流计算。很多人一听“电气热能流”就以为是把电网潮流、天然气网水力计算、热网水力热力计算三个东西分别跑一遍然后拼在一起真上手做一遍就会发现事情没那么简单。一个区域综合能源系统里CHP机组同时连着电网和气网电锅炉同时连着电网和热网燃气锅炉同时连着气网和热网这些耦合设备把三条网络硬生生绑成了一个整体。你单独算电网潮流的时候气网供给CHP的燃料量是未知的单独算气网的时候CHP的耗气量又取决于电网的发电需求而热网的供水温度又反过来影响CHP的抽汽量。这篇内容就是基于我实际做过的一个Matlab仿真项目把整套电气热能流计算从数学模型、求解策略到代码实现完整梳理一遍。适合正在做综合能源系统方向的研究生、刚入行做园区能源规划的工程师以及想把多能流计算模型真正落地成可运行代码的人。我会把每一步为什么这么做、踩过哪些坑、参数怎么定都写清楚。1. 问题是先理清楚的为什么单算三条网不算真正意义的多能能流1.1 单能流系统已经有成熟工具为什么还要做联合能流计算先说个最直观的类比你手头有一个电网潮流程序、一个天然气网水力计算程序、一个热网水力热力计算程序三个程序各自都跑得很顺。但把它们拼在一起的时候问题就出来了。电网潮流需要知道气网供给CHP的天然气量从而确定CHP的发电功率和热功率气网需要知道CHP的耗气量才能做节点流量平衡热网需要知道CHP的供热量才能确定节点温度和管道流量。这三个量互为因果谁都没法先独立确定。这就是所谓的多能耦合的内在复杂性。严格来说只有当耦合设备两侧的变量都被当作未知数、放进同一个迭代框架里求解才能叫计及多能耦合的能流计算。如果只是把CHP当作电网里的一个PQ节点气网里当作一个固定负荷那本质上还是三个孤立的单能流问题谈不上耦合。另外单能流系统的求解工具已经非常成熟比如电网潮流有Matpower天然气网有管网仿真软件热网也有专门的热力计算工具。但找遍市面没有一个现成工具能一次性求解电、气、热三网联立的能流方程。原因其实也简单这三个网络的物理特性差异太大。电网是代数方程主导求解讲究收敛性和电压约束气网是压力驱动的非线性方程跟电网形式有相似之处但物理含义完全不同热网更特殊它被水力方程和热力方程双重绑定先要解出各管道的流量分布再在那个基础上求解温度的传播与混合。物理特性差异这么大就很难用一套通用的求解器统一处理。1.2 一个典型的区域综合能源系统拓扑长什么样在动手之前先把对象系统画清楚。我这次仿真的典型区域综合能源系统由三部分组成一个IEEE 9节点电网可以外购电也可以由CHP发电一个6节点天然气网气源经过加压站、管道供给CHP和燃气锅炉一个8节点热网包含若干供热管道、换热站和热负荷节点。这三块之间的连接点就是耦合设备的安装位置。CHP机组挂在电网3号节点、气网4号节点、热网5号节点上燃气锅炉挂在气网5号节点和热网7号节点上电锅炉挂在电网7号节点和热网6号节点上。这样三个网络通过耦合设备形成了一个闭环气网给CHP供气CHP发电并产热发电送入电网、热量送入热网电网再给电锅炉供电、电锅炉补充供热。这个拓扑结构比较典型涵盖了综合能源系统里最常见的三种耦合设备CHP气转电热、燃气锅炉气转热、电锅炉电转热。做代码实现的时候这几个设备的模型就得分别写清楚因为它们的工作特性完全不同。我后面的算例、代码结构也都是围绕这套拓扑展开的。2. 电气热三网联立求解的数学建模核心方程式逐一拆解2.1 电网潮流的极坐标形式与节点类型处理电网部分我采用的是最经典的新ton-Raphson极坐标潮流算法。对每个PQ节点有两个待求的误差量有功功率误差和无功功率误差对每个PV节点只有一个有功功率误差电压幅值给定。方程组的核心表达就是功率平衡方程P_i V_i * Σ V_j * (G_ij * cos(θ_ij) B_ij * sin(θ_ij)) Q_i V_i * Σ V_j * (G_ij * sin(θ_ij) — B_ij * cos(θ_ij))这里的G_ij和B_ij对应支路导纳的实部和虚部θ_ij是节点i和j的相角差。常规潮流里平衡节点承担全网功率差额典型做法是设一个外网等效节点给定1.04的电压幅值和0度相角。做综合能源系统能流计算时电网部分有一个必须注意的改动CHP机组不能简单当作固定的PQ节点或PV节点。因为CHP的发电量是气网供气量和自身热电比共同决定的结果在气网没有收敛之前你根本不知道CHP能发多少电。所以在联合求解框架里CHP的并网节点一般当作P节点处理有功功率由耦合迭代给出电压幅值由无功调节维持无功功率不做限制。这样处理的好处是电、气、热三个网络的变量可以被耦合方程统一调度坏处是雅可比矩阵的结构和形成方式要跟着改不能直接照搬标准潮流程序的节点分类逻辑。2.2 天然气网稳态能流的Weymouth方程与节点压力计算天然气网的稳态能流计算核心是管道的Weymouth方程。高压输气管道中气体的质量流量和管道两端压力的平方差成正比这一点和电流与电压差的关系有些类似但又不太一样。一般表达式写作F_ij sign(π_i — π_j) * k_ij * sqrt(|π_i² — π_j²|)其中π_i和π_j是节点i、j的压力k_ij是与管道内径、长度、压缩因子、气体性质有关的常数。节点流量平衡方程为Σ F_ij L_i_gas该式表示流入节点的天然气净流量等于该节点的负荷包括天然气负荷和耦合设备耗气量。气源节点是压力已知的边界节点其他节点的压力是待求变量。在Matlab实现里这个非线性方程组的求解也使用牛顿法。雅可比矩阵的构造逻辑是对每一条管道支路求F_ij对π_i、π_j的偏导数对每个节点累加所有关联支路的偏导量。因为Weymouth方程的导数分母带有sqrt项当节点压力差特别小或者接近零时导数计算容易数值溢出。我后面会专门讲这个问题这里是初学者最容易忽略的。2.3 热网水力-热力双重耦合的方程体系热网的建模是三个网络里最容易写乱的一个因为变量最多、量纲也最容易出问题。热网的能流方程分为两个层次水力方程描述管道流量和水泵节点的压力平衡关系。其首要约束是节点流量连续性流入某节点的体积流量等于流出该节点的流量与供热负荷的流量需求相平衡。其次是环路压力方程在闭合的供热环中沿回路的压力损失之和为零这一组方程将管道阻力与流量平方和管径相关和网络拓扑紧密耦合。热力方程则关注温度场每个节点的供水温度根据上游管道混合温度确定节点处的热量平衡方程可写为Φ_i C_p * m_i * (T_s,i — T_r,i)其中C_p是水的比热容m_i是流经该节点的工质流量T_s,i和T_r,i分别是该节点供回水温度。管道沿线的温度损失方程采用指数衰减模型T_end T_amb (T_start — T_amb) * exp(-λ * L / (C_p * m))这里的λ是管道单位长度的热损失系数L是管道长度。水力方程的结果即各管道流量分布是热力方程求解的前提。所以在求解热网时通常先解水力方程得到m的分布再解热力方程确定温度分布。两者合在一起才算完成了热网的能流计算。这个先水力、后热力的顺序不能颠倒否则流量未知热量方程根本无从下手。2.4 耦合设备建模CHP、燃气锅炉、电锅炉的参数化处理耦合设备是多能流计算的灵魂所在。我在代码里实现三种设备的模型分别对应不同数学形式CHP机组采用典型的背压式机组模型数学模型可以表达为发电功率P_e和供热功率H_heat都与燃料消耗量F_gas呈线性关系。即P_e η_e * F_gas * LHV_gasH_heat η_h * F_gas * LHV_gas。其中η_e是发电效率η_h是热回收效率LHV_gas是天然气低位热值。背压式CHP的热电比是固定的这一特性使得电功率和热功率之间呈严格比例关系。在联立求解过程中可根据系统的电、热负荷需求直接确定燃料消耗量进而实现电网、气网和热网之间的信息交换与约束统一。燃气锅炉模型是最简单的就是一个转换效率η_gb天然气输入F_gas热量输出H_heat η_gb * F_gas * LHV_gas。电锅炉模型类似只是输入变成了电能P_e输出热量H_heat η_eb * P_e。这三个设备在Matlab里分别写成三个函数输入输出都是结构体变量。这样迭代求解的结构就能统一耦合设备从电网取/送出的电功率影响电网潮流方程从气网取用的气功率影响气网流量平衡方程向热网注入的热功率影响热网温度计算。设备模型越真实迭代间的变量传递就越复杂但原则是清晰统一的每种耦合设备本质上都是三条网络之间的一条数据通道。3. 求解策略怎么选统一求解法与分解迭代法的取舍3.1 两种主流框架的对比与适用场景多能耦合能流的基本求解策略学术界和工程上无非两条路统一求解法和分解迭代法。统一求解法将电网、气网、热网的所有方程放在一起形成一个大规模非线性方程组一次性用牛顿法迭代求解。好处是收敛性好尤其是当耦合程度很强比如CHP占系统发电比例很高时不容易出现发散。缺点是雅可比矩阵规模大、结构复杂编程量大且气网、热网的方程特性差异大矩阵容易病态。更要命的问题在调试阶段一旦不收敛你很难判断问题出在电网方程、气网方程还是热网方程上。分解迭代法的思路则是三个网络各自独立求解耦合变量在外层做一个迭代循环。迭代流程是先给定耦合设备的初始运行状态比如假设CHP发电功率和热功率分别求解电网、气网、热网的能流然后根据三网的结算结果更新耦合设备的运行状态再重新求解三网。如此循环直到前后两次迭代的耦合变量偏差小于收敛阈值。我在实际项目里用的是分解迭代法改进版本。原因很直接编程调试灵活哪一层出问题一目了然而且每条子网的求解函数可以单独测试跟已知结果的单能流系统对比验证大大降低了出错概率。缺点是单纯分解迭代在强耦合场景下可能迭代震荡或者收敛极慢。针对这个问题我给耦合变量的更新加了松弛因子效果非常明显X_new X_old α * (X_calculated — X_old)这里α取0.5到0.8之间比较稳妥。α太大会导致震荡发散α太小则收敛速度慢。当CHP的出力占比超过系统总出力的30%时保守一点取α0.5占比较低时取0.8可以兼顾速度和稳定性。3.2 初值策略的选择直接影响能否收敛的关键一环初值设定是这个项目里我踩过最多坑的地方。电网潮流的初值相对简单平启动即可所有PQ节点电压幅值取1.0相角取0。但气网的初值极为敏感管道Weymouth方程里有sqrt(|π_i² — π_j²|)项如果给定初值使得节点间压力差过大或者节点压力过低迭代过程中平方差出现负值或者开方项数值异常就直接挂了。我的做法是根据气源压力和工作压力范围先预估一下整个气网的基准压力。假如气源节点压力是2.0MPa末端节点压力大概在1.2MPa附近那所有节点的初值统一取1.6MPa别去费心思逐个节点估算。这样做的原理是牛顿法的局部收敛特性要求初值不能偏离真解太远取平均压力能保证初始压差在合理范围内从而避免开方项出现非法值的概率。热网的初值也类似。供水温度取设计值的上端比如100°C回水温度取80°C所有节点统一取这两个基准值。先解水力方程得到流量分布再在流量基础上解热力方程温度场的迭代比较温和。判断收敛的标准我设为三网的最大不平衡量小于1e-6标幺值。注意是三个网络的不平衡量都同时满足不是只看某一个网络的收敛情况。因为耦合系统的特点就是各网络相互牵制只收敛其中一个没有意义。3.3 解耦迭代的收敛判据与松弛更新逻辑外层耦合迭代的流程我写成下面的伪代码逻辑初始化假设CHP的发电功率P_e、热功率H_heat燃气锅炉的供热量电锅炉的耗电功率电网潮流求解把CHP作为P节点有功给定电锅炉作为负荷节点求解电网模型得到各节点电压、相角和各耦合节点的结算结果气网能流求解根据CHP和燃气锅炉的耗气量由第一步的初始假设或上一次迭代结果折算作为气网负荷求解气网模型得到修正后的节点压力和可用供气量热网能流求解根据CHP、燃气锅炉、电锅炉的热功率注入与热负荷需求求解热网模型得到供热温度和流量根据气网结算得到的流量反推CHP的发电功率上限并修订热功率值再根据热网的实际供热需求修正CHP和电锅炉的分配判断修正量与上次迭代的偏差是否小于收敛阈值不满足则返回步骤2继续迭代所有偏差满足阈值后输出全局能流结果这套流程写出来看着简单实际操作中容易忽略一个点步骤5的耦合变量更新要带上松弛因子而且每个设备的更新松弛因子最好单独设置。比如CHP的更新比较稳定α可取0.8电锅炉因为受电网和热网双重约束振荡风险高α取0.5。4. Matlab代码实现全过程从函数架构到算例验证4.1 全套代码的文件结构与核心数据流设计下面是这次项目全套Matlab代码的文件组织方式每个文件的功能都做了拆分互不干扰system_main.m 主程序入口控制整体迭代流程 load_system_data.m 读入电、气、热三网的拓扑参数与负荷数据 build_ybus.m 根据支路参数构建电网节点导纳矩阵 electric_powerflow.m 电网潮流求解牛顿-拉夫逊法 gas_powerflow.m 天然气网能流求解牛顿法 hydraulic_solve.m 热网水力平衡方程求解 thermal_solve.m 热网热力方程求解 coupling_update.m 耦合设备变量更新与松弛处理 plot_reports.m 三网能流结果的可视化输出在设计数据流的时候我建议所有中间变量和结果都用结构体变量按网络分类封装。三个大结构体sys_elec、sys_gas、sys_heat分别存放各网络物理模型所需的数据比如节点参数、支路参数、求解结果。耦合设备的参数单独放在struct_coupling结构体里。这个架构的好处是每段代码的输入和输出都是清晰的数据结构函数之间的耦合度被降到最低。比如electric_powerflow函数只接收sys_elec结构体和耦合节点的注入功率只返回更新后的电压和相角结果它不需要关心气网压力是多少、热网供水温度是多少。模块解耦的好处在问题排查阶段优势尤其明显电网不收敛时只需要检查电网这部分的网络连接和参数配置就行不会被气网和热网干扰。4.2 核心求解函数关键技术点牛顿法雅可比矩阵的组装先看电网潮流求解的核心循环。标准牛顿法潮流的核心就是构造雅可比矩阵并反复求解增量方程for iter 1:iter_max % 计算功率不平衡量 [dP, dQ] compute_power_mismatch(V, theta, Ybus, S_spec); mismatch [dP; dQ]; if max(abs(mismatch)) tol break; end % 构造雅可比矩阵 J build_jacobian(V, theta, Ybus); % 求解增量 dx J \ mismatch; % 更新状态变量 theta theta dx(1:n_node); V V dx(n_node1:end) .* V; end这段代码看起来简单但有几个细节需要强调。雅可比矩阵的组装公式推导可以参考任何一本电力系统分析教材但在Matlab里有一个程序效率的大坑如果直接用两个for循环逐个元素填充雅可比矩阵三节点系统跑起来没问题一旦节点数超过30每一次迭代都会因为矩阵组装太慢让人等到怀疑人生。更好的做法是利用稀疏矩阵的索引批量赋值或者直接使用vectorized方式按节点对数组同时计算。我实测过同样的IEEE 9节点算例向量化后的代码比双重循环快至少一个数量级这个差距在蒙特卡洛批量计算或者考虑多场景规划时会体现得非常明显。再来看看气网的求解核心。气网求解的难点其实不在牛顿法的框架而在于Weymouth方程导数的正确处理for iter 1:iter_max % 计算节点流量不平衡量包括耦合设备负荷 node_mismatch compute_gas_mismatch(Pressure, gas_loads, pipe_k, pipe_connect); J build_gas_jacobian(Pressure, pipe_k, pipe_connect); deltaP J \ node_mismatch; Pressure Pressure deltaP; if max(abs(node_mismatch)) tol break; end end这里build_gas_jacobian函数有一个关键逻辑对于管道ij流量F_ij对节点压力P_i的偏导要分为两种情况计算。当P_i大于P_j时导数为正值为k_ij * P_i / sqrt(|P_i² — P_j²|)。当P_i小于P_j时导数为负。由于开方项位于分母中当压力差趋近于零时雅可比矩阵元素会趋向无穷大该数值特性极易导致迭代不收敛。为避免此问题需在程序中对压力差设置最小限值例如当|P_i² — P_j²|小于某个阈值我取1e-6时将该导数按上限值截断处理。曾经由于忽略该细节导致程序反复发散排查了很久才意识到问题所在。4.3 一个典型算例的测试结果与误差分析我用的算例配置如下CHP安装在三网交汇的关键节点电功率输出为60MW热功率输出为50MW燃气锅炉的热功率为30MW电锅炉的热功率为20MW电网侧的总电负荷约为110MW外网联络线提供剩余约30MW功率气源节点压力设定为2.0MPa两个气网负荷节点的天然气流量需求分别约为4000和2500 m³/h。经过分解迭代求解整个系统在14次外层迭代后收敛松弛因子取0.6各网络的不平衡量均达到10⁻⁶级别。电网的平衡节点有功出力为-28.5MW即从外网购电28.5MW各节点电压幅值均在0.98到1.03之间满足常规设计要求。气网的节点压力从气源2.0MPa下降到末端约1.35MPa压降符合管道长度和流量的物理预期。热网的供水温度经过管网衰减后最远负荷节点比热源出口下降了约9°C回水温度分布也比较均匀。为了验证程序的正确性我做了两组对照实验。第一组是单一电网潮流对照把CHP当作固定出力节点运行电网潮流程序与Matpower的计算结果对比节点电压偏差小于1e-6说明电网模块实现无误。第二组是气网对照把已知流量和压力参数的6节点气网用商业管网仿真软件校核结果压力偏差在2%以内。这样逐模块验证过的程序再去做三网联算才敢说结果是可信的。5. 实战中一定会遇到的坑三个典型问题与排查实录5.1 气网开方项异常造成不收敛的排查过程项目最早期调试阶段气网部分一直报出NaN程序直接崩溃。当时把错误定位到sqrt处但百思不得其解因为从物理概念上节点压力平方差为负似乎不该出现。后来通过打印每次迭代的节点压力值发现迭代第一次就出现某个节点压力值跳到接近零随后它与相邻节点间的压力平方差直接变负开方返回复数。根因其实是两个问题叠加。第一个是初值给得太差我把节点初值都设成0.5MPa但气源压力是2.0MPa迭代前期压力调整跨度大部分节点被牛顿法推到了不合理区域第二个是没有在Weymouth方程里加上数值保护当开方根内部出现极小负数时直接就得到NaN了。修正方案是两个方向同时入手初值统一取工作压力范围的中间值同时在Weymouth公式里加入限幅保护用max(eps, abs(p_i²—p_j²))来保险。5.2 耦合变量振荡不收敛为什么需要逐设备设置松弛因子框架刚跑起来的时候出现了比较有意思的现象电压和压力都稳定但CHP的出力在相邻两次迭代之间来回跳幅度还不小±8MW左右导致外层循环无法收敛。直觉解决方案是把全局松弛因子调小发现效果有限因为电网、气网模块本身是精确计算的真正振荡的来源是耦合变量更新逻辑——我先更新CHP发电量再根据它来更新耗气量而耗气量又马上被气网结算结果修正形成正反馈回路。后来我把CHP、燃气锅炉、电锅炉的松弛因子分别独立设置CHP为0.6燃气锅炉为0.5电锅炉为0.4。因为电锅炉同时受电网和热网两侧约束其更新最容易引起波动。结果外层迭代从完全不收敛变成14次收敛效果立竿见影。这个教训让我意识到在多能流迭代里不要贪图代码简单就用一个统一的松弛因子不同耦合设备的时间常数和灵敏性完全不同分开设置才是合理做法。5.3 单位换算的隐蔽炸弹天然气标方与热功率的换算最后一个要提醒的坑是单位。这是我检查看了半天才发现的低级错误但后果很严重。天然气网里气负荷的单位习惯用m³/h标方而电气热联立方程里CHP的耗气量需要换算成MW热功率才能和电网功率平衡方程式统一量纲。天然气低位热值LHV一般取36 MJ/m³换算关系是1 m³/h天然气约等于0.01 MW。计算方法是36 MJ/m³乘以1000除以3600秒结果正好约等于10 kW也就是0.01 MW。这个数值关系看起来很直白但如果代码里直接用一个常数去乘除很容易在某处漏掉或者写错一个数量级。更隐蔽的是锅炉效率、CHP效率这些百分比数在单位换算时也要参与运算一旦顺序弄错结果就可能出现CHP耗气量比气网总供气量还大的离谱局面。我的建议是在所有计算之前定义一组统一的基准值和转换函数把天然气质量流量、标准体积流量、热功率三套单位全部转换成统一的标幺值或者统一的MW制后再进入迭代计算。程序内部只用一套单位制只在输入和输出层做转换。6. 一些个人经验和可以继续扩展的方向做这个项目最大的体会是多能耦合能流计算耦合两个字才是重点。不要把它当成三个单能流程序的简单拼接而是要在求解框架设计时就充分考虑耦合变量的传递机制。初期我犯的错就是把80%的精力花在三个网络各自的算法实现上最后联调阶段才发现真正的难点在更新与收敛策略上。另外一个实用的经验是Matlab里做这类项目一定要充分利用结构体和稀疏矩阵。结构体让三网的参数管理变得清晰而稀疏矩阵让大规模节点的求解效率成倍提升。如果只是为了能跑通用全矩阵也凑合但如果你以后想把算例规模扩大或者在项目基础上做优化调度、多场景分析全矩阵的性能瓶颈就会非常明显。这套代码后续可以扩展的方向还挺多一是把稳态能流扩展为准稳态时序仿真模拟一天24小时或者四季的更替场景二是在能流计算基础上叠加经济优化目标比如以运行成本最小为目标函数调整CHP出力比三是加入储能设备模型让气网、热网具备调节能力看看储能对多能互补系统能流分布的影响。如果你已经在做或者准备做综合能源系统的仿真分析欢迎一起交流具体实现中的细节。