计及多能耦合的区域综合能源系统电气热能流联合求解与Matlab实现

1. 为什么区域综合能源系统的热能流必须“电热气联算”我最初拿到“计及多能耦合的区域综合能源系统电气热能流计算研究Matlab代码实现”这个题目时第一反应是这不过就是把电网潮流、气网水力、热网水力三个模块各算各的最后拼在一起出个报告。真正动手以后才发现如果只是这样拼算出来的结果根本没法用。原因在于区域综合能源系统里的关键设备往往横跨多个网络。燃气轮机一面从天然气网络取气一面发出电功率余热还要进热力管网电锅炉把电功率转成热功率直接改变电网注入功率的同时又给热网增加了一个热源P2G设备则把富余电功率转化为天然气。这些设备就是所谓的“耦合元件”它们的存在使得三个网络的状态量相互牵制。最典型的矛盾场景是电网调度要求燃气轮机满发但气网节点压力在高峰时段跌到下限机组实际进气量不足出力根本顶不上去。如果电网潮流和气网水力分开计算这个问题永远不会暴露。所以这里的核心关键词是“计及多能耦合”。它的意思是三个子网络不再各自独立求解而是把电、气、热统一放进同一个状态变量向量里通过统一迭代求解得到全系统各节点的电压、相角、气压、温度、热功率等运行状态。这也是“综合能源系统”区别于“各能源系统简单叠加”的真正价值所在。做这个项目的人通常是这么几类研究综合能源系统规划与运行的硕士研究生、做园区级能源站设计的技术人员、以及需要评估多能互补方案可行性的工程师。如果你属于其中之一那么本文介绍的建模思路、求解框架和Matlab代码实现路径可以直接作为你自主实现的热能流计算程序的起点。2. 电气热能流联合求解的物理模型搭建2.1 电网模型极坐标潮流是公认的锚点电网部分采用传统极坐标潮流方程。对每个节点i需要满足有功注入方程 [ P_i U_i \sum_{j \in i} U_j (G_{ij}\cos\theta_{ij} B_{ij}\sin\theta_{ij}) ]无功注入方程 [ Q_i U_i \sum_{j \in i} U_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]节点类型沿用电网潮流经典分类平衡节点Slack、PV节点、PQ节点。在综合能源系统场景里燃气轮机和发电机通常设为PV节点负荷和电锅炉设为PQ节点。平衡节点承担全系统有功误差这个角色一般由与大电网的联络节点或者容量最大的常规机组扮演。需要特别强调的是电网模型的时间尺度选取。电气热能流计算属于稳态能流计算范畴所以电网、气网、热网全部采用稳态模型不考虑惯性环节和动态过程。这一点和机电暂态仿真有本质区别后者要求解微分方程而这里只需要解代数方程组。2.2 天然气网模型Weymouth方程是管流计算的基石天然气网络稳态计算的核心是管道流量方程。对一条连接节点i和节点j的管道气体流量通常用Weymouth方程描述[ f_{ij} C_{ij} \sqrt{p_i^2 - p_j^2} ]其中 (C_{ij}) 是与管径、长度、气体性质、摩擦系数有关的常数。当 (p_i p_j) 时流量方向反转所以更严谨的写法是[ f_{ij} C_{ij} \cdot \mathrm{sign}(p_i^2 - p_j^2) \cdot \sqrt{|p_i^2 - p_j^2|} ]每个气网节点的流量守恒方程为[ \sum_{j \in i} f_{ij} F_{Gi} - F_{Li} - F_{CHP,i} 0 ]其中 (F_{Gi}) 是气源注入流量(F_{Li}) 是普通气负荷(F_{CHP,i}) 是燃气轮机耗气流量。气网节点也分两类一类是气源节点类似于电网的平衡节点压力给定如2.5 MPa负责平衡全网流量另一类是负荷节点注入流量给定压力待求。实际天然气网络中通常还有压缩机站压缩机两端压力的关系一般表示为 (p_{out} k_c \cdot p_{in})其中 (k_c) 为压缩比。加入压缩机后气网雅可比矩阵的稀疏结构会发生变化这一点在后面的代码实现部分我会重点讲。2.3 热网模型水力与热力是两套必须同时满足的方程热力网络比电、气网络多一层复杂性——它要同时求解水力方程和热力方程。水力方程包括两类。第一类是节点流量连续性方程对每个热网节点流入流量等于流出流量即 (\sum A_{ij} m_j m_{ext,i})其中 (m_j) 代表管道质量流量(A_{ij}) 是节点-支路关联矩阵(m_{ext,i}) 是节点注入或流出的外部流量。第二类是回路压降方程沿热网中每个闭合回路所有管道压降代数和为零即 (\sum B_{ij} h_j 0)。管道压降 (h_j) 与流量 (m_j) 的关系一般用 Darcy-Weisbach 公式[ h_j K_j m_j |m_j| ]其中 (K_j) 是与管径、长度、摩擦系数相关的阻力系数。热力方程部分需要区分供水网络和回水网络。热源从回水网络取低温水加热后送入供水网络热用户从供水网络取高温水换热后回流到回水网络。因此每个节点有两个温度状态变量供水温度 (T_s) 和回水温度 (T_r)。节点热功率平衡方程[ \phi_i C_p m_{ext,i} (T_{s,i} - T_{r,i}) ]其中 (\phi_i) 是节点注入或消耗的热功率(C_p) 是水的比热容。热网节点的温度混合方程和管道温降方程也需要一并纳入[ T_{end} (T_{start} - T_a) e^{-\frac{\lambda L}{C_p m}} T_a ]这里的 (T_a) 是环境温度(\lambda) 是管道单位长度传热系数(L) 是管长。温降方程让热网方程组的非线性程度比电网还高对初值的选择也更敏感。3. 耦合环节建模把三个网络“缝”在一起的关键3.1 能源集线器视角下的设备建模耦合环节的建模方式有很多种我在项目中采用的是“能源集线器 设备效率矩阵”的方法。核心思路是把燃气轮机、燃气锅炉、电锅炉、P2G等设备视为能源转换环节用效率系数描述输入-输出关系再用这些关系把子网络之间的耦合变量显式表达出来。以最常见的热电联产机组CHP为例。设其消耗的天然气流量为 (F_{CHP})发出的电功率 (P_{CHP}) 和热功率 (\Phi_{CHP}) 分别为[ P_{CHP} \eta_e \cdot LHV \cdot F_{CHP} ][ \Phi_{CHP} \eta_h \cdot LHV \cdot F_{CHP} ]其中 (LHV) 是天然气低位热值(\eta_e)、(\eta_h) 分别为发电效率和热回收效率。在气网的节点方程里(F_{CHP}) 是一个负荷项在电网的节点方程里(P_{CHP}) 是一个电源注入项在热网的节点方程里(\Phi_{CHP}) 是一个热源项。同一个变量 (F_{CHP}) 同时出现在三个子网络的方程中这就是“耦合”的数学本质。电锅炉的建模更简单(\Phi_{EB} \eta_{EB} \cdot P_{EB})(P_{EB}) 是电网节点的电负荷增量(\Phi_{EB}) 是热网节点的热源增量。P2G设备则是逆过程(F_{P2G} \eta_{P2G} \cdot P_{P2G} / LHV)。3.2 耦合变量在统一迭代中的处理方式在统一求解框架下耦合变量不需要单独迭代更新而是作为中间变量嵌入各子网络方程中。具体做法是每次迭代开始时根据当前各网络状态变量计算所有耦合设备的输入输出功率然后带入电网、气网、热网的残差方程最后组装成统一雅可比矩阵进行修正。这样做有一个明显的好处——不需要为耦合设备单独设计迭代环节既避免了嵌套迭代的慢收敛也减少了代码实现的复杂度。代价是雅可比矩阵的维度和非零元数量显著增加对矩阵组装和求逆的计算资源要求更高。对于区域级系统电网几十个节点、气网十几个节点、热网十几个节点Matlab的稀疏矩阵运算完全能够胜任。3.3 耦合环节的可行域约束除了能量转换关系耦合设备还需要考虑运行约束。比如燃气轮机有最大和最小出力限制电锅炉有爬坡限制在稳态计算里不做动态考虑但可以设上下限P2G设备通常有容量限制。这些约束在热能流计算中一般作为节点类型判断的依据当某节点类型对应的控制量越过边界时该节点需要从PV转为PQ或者从给定出力的节点转为边界受限节点。这个处理方式和电力系统潮流计算中发电机无功越限的处理思路完全一致。我在代码里预留了“节点类型切换”的函数接口实测下来包含边界条件的算例会增加几次迭代但不会影响最终收敛。4. Matlab代码实现核心逻辑从数据结构到统一牛顿法4.1 数据结构设计面向对象还是结构体数组做综合能源系统仿真Matlab的数据结构设计直接决定后期扩展的难易程度。我尝试过面向对象classdef也尝试过纯结构体数组。对于这种规模的项目我的建议是不要一上来就搞面向对象用结构体数组足够而且后期调试更直观。推荐的数据组织方式如下% 电网数据 grid.bus struct(id, {1;2;3}, ... type, {slack;PV;PQ}, ... P, {0; 0.5; -0.3}, ... Q, {0; 0.2; -0.1}); grid.branch struct(from, {1;1;2}, ... to, {2;3;3}, ... R, {0.01; 0.02; 0.015}, ... X, {0.05; 0.08; 0.06}); % 气网数据 gas.node struct(id, {1;2;3}, ... type, {source;load;load}, ... p, {2.5; 1.8; 1.6}, ... % p是节点气压MPasource节点为给定值 F_load, {0; 0.1; 0.08}); % 气负荷kg/s % 热网数据 heat.node struct(id, {1;2;3}, ... type, {source;load;load}, ... Ts, {110; 95; 90}, ... Tr, {70; 65; 60}, ... phi, {1.0; -0.4; -0.5}); % 热功率MW正为注入这里的关键设计思路是每个网络的数据包单独存放耦合设备CHP、电锅炉、P2G单独建一个结构体数组coupling.chp struct(bus_id, {2}, gas_node, {1}, heat_node, {1}, ... eta_e, {0.35}, eta_h, {0.45}, ... F_max, {0.2}, F_min, {0.02}); coupling.eb struct(bus_id, {5}, heat_node, {2}, ... eta, {0.95}, P_max, {0.5});这样做的原因是统一迭代时需要频繁读取耦合设备的当前工作点结构体数组的字段访问速度比嵌套cell数组快一个数量级而且代码可读性好后续加新设备类型只需要在coupling里加一个字段块。4.2 统一残差方程组的组装统一牛顿法的核心是把所有子网络的残差方程放进同一个向量 (R(x)) 中对状态变量 (x) 求雅可比矩阵 (J)然后求解修正方程 (J\Delta x -R(x))。状态变量向量定义为[ x [\theta, U, p, m, T_s, T_r]^T ]其中 (\theta) 和 (U) 是电网相角与电压幅值不含平衡节点电压相角(p) 是气网非气源节点压力(m) 是热网管道质量流量(T_s) 和 (T_r) 是热网节点供回水温度。残差向量按以下顺序组织function R unifiedResidual(x, grid, gas, heat, coupling) % 第一步从x中拆出各子网络状态量 n_gbus length(grid.bus); idx_theta 2:n_gbus; % 平衡节点相角不参与迭代 idx_U n_gbus1:2*n_gbus; % 电压幅值PV和PQ节点 % …… 类似方式切出气压、流量、温度区间 % 第二步计算耦合设备当前工作点 [P_chp, Q_chp, Phi_chp, F_chp] calcCoupling(coupling, x); % 第三步计算各子网络残差 R_grid calcGridResidual(x_theta, x_U, P_over, Q_over); R_gas calcGasResidual(x_p, F_load_total); R_heat_hyd calcHeatHydraulicResidual(x_m); R_heat_ther calcHeatThermalResidual(x_Ts, x_Tr, x_m); R [R_grid; R_gas; R_heat_hyd; R_heat_ther]; end这里最容易被忽视的是“边界条件穿越”电网的平衡节点、气网的气源节点、热网的平衡节点通常取一个已知压力或温度的热源不参与状态变量求解它们的变量值恒定对应残差方程也要从方程组中移除。否则雅可比矩阵奇异迭代必发散。4.3 雅可比矩阵组装解析求导优于数值微分雅可比矩阵的组装是整个程序里最影响性能和收敛性的部分。两种方案第一种是数值微分即用有限差分近似雅可比矩阵。实现简单每列需要多算一次残差计算量大而且差分步长选择不当会严重影响收敛精度。我的实测结果是对小系统能收敛但迭代次数明显增加大系统容易出现震荡。第二种是解析求导。对每个残差方程直接推导对各状态变量的偏导数。虽然推导过程繁琐但组装出来的雅可比矩阵精度高、收敛速度快。比如电网潮流残差对电压相角的偏导、气网Weymouth方程对气压的偏导、热网温降方程对温度与流量的偏导这些都是标准公式可以在文献里直接查到。这里我给一个示例说明如何用Matlab稀疏矩阵组装雅可比矩阵function J assembleJacobian(x, grid, gas, heat, coupling) n length(x); J sparse(n, n); % 电网部分的雅可比分块 [dP_dtheta, dP_dU, dQ_dtheta, dQ_dU] ... calcGridJacobian(grid, x_theta, x_U); J(grid_rows, grid_cols) ... J(grid_rows, grid_cols) [dP_dtheta, dP_dU; dQ_dtheta, dQ_dU]; % 气网部分的雅可比分块 [dF_dp] calcGasJacobian(gas, x_p); J(gas_rows, gas_cols) J(gas_rows, gas_cols) dF_dp; % 热网部分…… % 关键耦合设备对雅可比矩阵的贡献 [dR_dF_chp] calcCouplingJacobian(coupling, x); J(chp_rows, F_chp_cols) J(chp_rows, F_chp_cols) dR_dF_chp; end耦合设备的雅可比贡献是最容易漏掉的部分。以CHP为例气网残差方程里包含 (-F_{CHP}) 项而 (F_{CHP}) 又依赖于CHP的进气压力实际中因为流量与压力差有关严格说还受气网压力影响但稳态计算中通常简化为给定功率对应给定进气量不考虑气压对进气量的反作用。在简化模型中 (F_{CHP}) 是常数雅可比贡献为零在考虑气压影响的扩展模型中(F_{CHP}) 对气网节点压力的偏导就不是零必须填入对应位置否则迭代会产生错误的修正方向。4.4 初值选择与迭代收敛控制个人经验中最影响综合能源系统热流计算成败的是初值选择。电网部分还好电压初值取1.0 p.u.、相角取0是常规做法气网压力初值取气源压力值的80%左右往往比直接取1.0更稳妥热网管道流量初值需要根据节点热负荷估算否则温度方程的收敛很容易发散。迭代循环建议用如下框架x initState(grid, gas, heat); tol 1e-6; max_iter 60; for iter 1:max_iter R unifiedResidual(x, grid, gas, heat, coupling); if max(abs(R)) tol fprintf(已收敛迭代次数%d\n, iter); break; end J assembleJacobian(x, grid, gas, heat, coupling); dx -J \ R; % 阻尼因子防止修正过大导致震荡 alpha 1.0; while 1 x_new x alpha * dx; R_new unifiedResidual(x_new, grid, gas, heat, coupling); if norm(R_new) norm(R) break; end alpha alpha * 0.5; if alpha 1e-4 break; end end x x_new; end阻尼因子这一段是我实际调试中加上的。电气热能流耦合之后比纯电网潮流的非线性强很多尤其是当初值偏差较大时全步长修正极易导致温度或气压变量飞出物理可行域比如压力变成负数。加入阻尼后虽然极端情况下会增加迭代次数但收敛稳定性明显提升。5. 算例验证一个小型区域综合能源系统5.1 测试系统结构与参数为了验证程序正确性我搭了一个小型测试系统电网6节点、气网5节点、热网6节点。电网部分包含1台常规机组、1个CHP和3个负荷节点气网部分包含1个气源节点、1个CHP用气节点、1个燃气锅炉用气节点和1个普通气负荷节点热网部分包含CHP热源、燃气锅炉热源、2个热负荷节点和1个电锅炉。需要说明的是测试系统规模虽小但包含了全部三种耦合设备类型CHP、燃气锅炉、电锅炉足以检验统一迭代框架的正确性。5.2 计算结果对比以下是其中一个工况的计算结果摘要状态量独立计算各网络分别求解统一联合求解状态量独立计算统一联合求解电网节点3电压/p.u.0.9820.961气网节点3压力/MPa1.721.58电网节点5电压/p.u.1.0040.997热网供水温度/°C103.697.2热网节点4流量/(kg/s)2.312.09CHP电功率/MW1.501.32两组结果在部分状态量上差异明显。原因在于独立计算时CHP电功率1.50 MW是按调度指令给定的不需要反推气网进气量而统一求解时气网在高峰负荷下节点压力下降CHP实际进气量对应发电能力只有1.32 MW两者之间的差值最终由常规机组增发补上。这个现象恰恰说明了“计及多能耦合”的价值不考虑气网约束的调度指令在实际运行中是落实不了的。迭代次数方面纯电网潮流收敛只需4次迭代加入气网和热网后统一迭代收敛需要11次最大残差下降过程平稳没有出现震荡。这说明统一牛顿法的收敛特性总体上是可靠的。6. 调试与避坑记录我在实现过程中踩过的几个关键坑6.1 单位制不统一是最隐蔽的错误来源电网里功率用MW气网里流量可能用kg/s、m³/h、标准立方米/天热网里热功率又是MW而质量流量又是kg/s。如果程序里不做单位归一化计算结果会出现量级错误而且很难排查。我在代码入口处做了严格约定所有功率量统一为MW气网流量统一为kg/s热网质量流量统一为kg/s温度统一为°C压力统一为MPa。所有的转换函数集中在unitConvert.m里后续加设备类型时只需要调用这些函数不要再在计算逻辑里临时换算。这个习惯帮我省了大量的调试时间。6.2 热网回水温度初值会给收敛“埋雷”热网方程中供水网络和回水网络的温度是分开求解的。回水温度初值如果取得过低比如直接取环境温度20°C热源节点的热功率平衡方程残差会非常大导致第一次迭代修正量超出物理范围。我的做法是用热负荷功率估算各节点回水温度初值T_r_init(i) T_s_init(i) - phi(i) / (C_p * m_init(i));其中 (m_init) 用节点热负荷除以供回水温差估算。实测这样给初值后热网方程基本不会发散。6.3 雅可比矩阵中的“零块”不要强行填满统一雅可比矩阵天然具有分块结构电网残差对气网变量、热网变量的偏导为零气网残差对热网状态变量偏导也为零。这些零块不需要计算直接让稀疏矩阵对应位置保持0即可。但有一点必须警惕虽然大多数子网络之间的偏导为零耦合设备的偏导会填到跨网络的位置上。比如CHP的热功率残差对电电网电压的偏导只通过电功率间接影响正常情况下为零但CHP的燃气流量 (F_{CHP}) 如果被建模为与电网出力强相关则在气网残差对电网相角的偏导位置上就会有一个非零元素。这个元素非常容易被漏掉我在调试时就是因为漏了这个非零元素导致气网压力一直震荡不收敛。6.4 压缩机节点与PV节点的“类型切换”逻辑要预留气网的压缩机通常维持出口压力恒定类似于电网的PV节点但当气源压力不足时压缩比会达到上限此时出口压力无法维持压缩机实际上变成了“流量给定”节点。我在代码中预留了类型切换逻辑参考电网无功越限的处理方式每次迭代后检查压缩机压缩比是否越限若越限则将节点类型切换为流量给定节点并重新组装雅可比矩阵。实测表明加入这一逻辑后系统在气源供应紧张工况下的收敛性明显改善。6.5 用“分步调试”代替“一步到位”我的调试顺序是先单独验证电网潮流程序与Matpower结果对比再单独验证气网水力计算与文献算例对比然后单独验证热网水力热力计算与已知稳态解对比最后才拼装耦合迭代。每个环节都单独测试通过后再联调。千万不要一上来就拼完整程序否则出错后根本没有参考基准只能大海捞针。具体来说可以先用小步骤验证耦合环节% 第一步固定热网和气网状态只让CHP的电网注入功率参与迭代 % 第二步固定电网和热网状态只让CHP的耗气量参与气网迭代 % 第三步三者同步迭代这种“坐标下降”式的过渡调试法能让你在每一步都能判断是哪个子网络出了问题。最后再分享一个小技巧在统一迭代的主循环里每一轮都把max(abs(R))打印出来观察残差下降曲线。如果出现残差先降后升再降的锯齿状先考虑阻尼因子是否过小如果是长期不降但震荡优先检查雅可比矩阵跨网络非零元素是否漏填。我见过不少人在这两个问题上耗掉大量时间其实问题往往都不复杂只是排查顺序不对。