搞区域综合能源系统的人几乎都绕不开电气热能流计算这块硬骨头。所谓电气热能流不是三条独立潮流的简单拼盘而是考虑多能耦合后电、气、热三种网络在一个系统里联立求解的运行状态。我之前在做一个园区级综合能源系统方案时因为要用Matlab实现计及多能耦合的电气热能流计算前前后后折腾了小两个月。这篇记录我会把建模思路、代码实现细节、调参踩坑经验都摊开来讲适合正在做综合能源方向研究的学生、刚入行的系统工程师以及打算自己搭稳态分析工具的朋友参考。1. 项目背景与核心思路拆解1.1 为什么区域综合能源系统需要“多能耦合”视角传统电力系统潮流计算只管电网天然气管道水力计算只管气网供热管网水力热力计算只管热网。但在区域综合能源系统里电网、气网、热网不是孤立的它们通过热电联产机组CHP、燃气锅炉、热泵、电锅炉、P2G设备联系在一起。任何一个网络内部的变化都会通过耦合设备传导到其他网络。举个最简单的例子冬季热负荷增加燃气锅炉的耗气量上升气网的节点流量重新分配气源压力发生变化如果CHP机组参与供热它的热出力提升往往带来电出力改变电网潮流随即变化。反过来如果电网发生故障导致热泵停机热网热源出力缺口就需要其他热源填补又回头影响气网用气量。这种“电-气-热”互相牵制的稳态关系正是电气热能流计算要解决的问题。因此项目标题里强调“计及多能耦合”不是概念上的装饰而是实实在在的数学建模要求。我们需要把三个物理网络的方程联立起来求解而不是分别计算后拼个结果。在Matlab里实现这套计算核心难点也在于此。1.2 电气热能流计算的本质从三个独立网络到一个耦合系统从数学上看电气热能流计算本质上是求解一组非线性方程组电网节点功率平衡方程、气网节点流量平衡方程、热网水力与热力方程以及耦合设备能量守恒方程。状态变量包括电网节点电压幅值和相角、气网节点压力、热网节点供水和回水温度、各耦合设备的内部运行变量。求解这组方程有两种主流思路。第一种是统一求解法把电、气、热所有方程和变量拼成一个大雅可比矩阵用牛顿-拉夫逊法整体迭代。第二种是解耦迭代法交替求解电网、气网、热网子问题每次迭代后更新耦合设备变量直到收敛。我自己最终选择的是统一求解法。原因很直接解耦迭代虽然可以复用成熟的潮流程序但迭代过程中容易出现震荡尤其当耦合设备比例较高时收敛性很难保证。统一求解法把所有耦合关系放在一个雅可比矩阵里物理上更清晰数学上收敛性也更好只是代码实现时需要仔细处理变量索引和雅可比矩阵的组装。2. 系统建模电气热网络数学模型与耦合环节2.1 电网潮流模型与牛顿-拉夫逊解法电网部分采用经典极坐标形式的牛顿-拉夫逊法。每个节点的有功和无功功率方程分别为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的电压相角差。节点类型依然分为PQ节点、PV节点和平衡节点。在区域综合能源系统中耦合设备通常作为电网的PQ节点或PV节点处理但如果CHP机组采用“以热定电”模式它的电出力会随热出力变化这时必须把CHP相关变量作为额外状态量加入统一求解框架而不能简单当PQ节点。在Matlab代码里我会用稀疏矩阵保存导纳矩阵Ybus节点功率偏差作为残差分量。计算电网残差时从状态向量中取出V和θ计算注入功率再减去给定值。% 电网功率残差计算 Vn V .* exp(1j * theta); Iinj Ybus * Vn; S Vn .* conj(Iinj); Pcal real(S); Qcal -imag(S); % 注意符号约定这里以注入为正 res_grid [Pcal(PQPV) - P_spec(PQPV); Qcal(PQ) - Q_spec(PQ)];这种写法在Matlab里非常直观调试时也容易对照教科书公式。2.2 天然气网络水力方程与气潮流计算天然气网络稳态计算的核心是节点流量平衡。每根管道流量与两端压力相关最常用的是Weymouth公式f_mn sign(π_m² - π_n²) · C_mn · sqrt(|π_m² - π_n²|)其中π_m和π_n为节点绝对压力C_mn为管道常数由管径、长度、摩擦因数、气体组分决定。sign项保证流量方向从高压流向低压。气网节点类型分为三类气源节点压力给定、负荷节点流量给定、以及平衡节点用来平衡全网用气误差。计算气潮流时有一点容易被忽略天然气是可压缩流体流量对压力呈平方关系方程非线性程度比电网功率方程还厉害。初值给不好牛顿法在第一轮就可能发散。我的做法是所有非气源节点初始压力设为气源压力的0.9倍并配合阻尼牛顿法来启动。气网残差写入统一求解器时还需要考虑压缩机等设备。简化处理中可以暂时忽略压缩机或者把压缩机等效为流量约束。对于研究型的能流计算先忽略压缩机也能验证核心算法是否跑通。% 气网管道流量计算Weymouth公式 dp_sq pi_gas(m) .^ 2 - pi_gas(n) .^ 2; flow(m,n) sign(dp_sq) .* C(m,n) .* sqrt(abs(dp_sq)); res_gas node_injection - node_demand - sum(flow, 2);2.3 热力网络水力-热力联合计算热网是最容易让人头疼的部分因为它包含水力工况和热力工况两个层次。水力工况求解各管道质量和节点压力热力工况求解节点供水和回水温度。水力模型的基本约束是节点质量流量守恒流入节点的质量流量等于流出节点的质量流量。另外每个供热环路的压力降代数和必须为零。简化计算时可以忽略环网只考虑枝状管网这样只需要质量守恒方程就能求出各管段流量。热力模型的核心是热功率方程节点热负荷φ_i等于该节点供回水温差乘以质量流量和比热容φ_i C_p · m_i · (T_s,i - T_o,i)管道温度沿程会下降考虑环境散热时用指数衰减公式T_end (T_start - T_amb) · exp(-λL / (C_p · m)) T_amb多根管道汇入同一个节点时混合温度按照能量加权平均T_mix ∑(m_j · T_j) / ∑m_j在热网中热源和热负荷都表现为对温度或热功率的约束。例如CHP机组的供水温度固定而热用户按照热负荷给定。热力方程既有线性也有强非线性和电网、气网耦合后对雅可比矩阵的数值精度要求更高。2.4 耦合设备的耦合方程与状态变量扩展多能耦合是通过耦合设备体现的需要为每一类设备建立数学方程。最典型的几类设备如下。CHP联产机组输入天然气输出电和热。简化建模可以写成F_gas P_elec / η_eQ_heat F_gas · η_h · LHV_ratio其中η_e是发电效率η_h是热回收效率。如果背压式CHP热电比固定抽汽式CHP则热电比可以调节需要增加一个抽汽系数变量。燃气锅炉天然气消耗与热出力成正比F_gas Q_heat / (η_gb · LHV)热泵消耗电功率产生热功率Q_heat COP · P_elecP2G设备消耗电能产生天然气作为气网的附加气源和电网新增负荷F_gas η_p2g · P_elec这些方程需要在统一求解时作为额外的残差方程加入同时把设备内部变量加入状态向量。例如CHP机组的热出力Q_heat和电出力P_elec热泵的耗电量P_elec都需要参与迭代。这样做的优势是耦合设备对三个网络的边界影响可以被同时感知不会出现解耦迭代时那种“气网算完才想起热网”的滞后效应。3. Matlab代码实现从数据到算法的落地细节3.1 代码整体架构与数据结构设计写Matlab程序最怕把所有代码堆在一个脚本里后面根本没法改。我做这个项目时把代码拆成四层结构数据定义层、网络模型层、求解器层、结果输出层。数据定义层用结构体组织每个物理网络一个struct。电网保存导纳矩阵Ybus、节点编号、负荷、发电机数据气网保存管道首末端节点编号、管道常数C、节点类型、压力初值热网保存管道参数、节点负荷、比热容、环境温度耦合设备保存连接拓扑和效率参数。data.grid.Ybus Ybus; data.grid.bus bus; % bus(:,1)编号, bus(:,2)类型, bus(:,3)V, bus(:,4)theta data.gas.pipe pipe; % pipe(:,1),pipe(:,2)首末端, pipe(:,3)C常数 data.gas.node gas_node; % [编号,类型,压力,负荷] data.heat.pipe heat_pipe; % [首端,末端,长度,直径,散热系数] data.heat.node heat_node; data.coupling.c 0; % 依次定义CHP、热泵、锅炉等 data.coupling.chp chp_data;数据结构定了后面的残差函数、雅可比函数、牛顿迭代主循环都能统一处理。我最想强调的一点是变量索引必须单独写一个函数负责把全局状态向量x拆成电网相角、电压幅值、气网压力、热网供水和回水温度、耦合变量。千万别在残差函数里靠魔法数字取位置否则改一个设备类型索引全部乱套。function [theta, V, pi_g, Ts, Tr, x_cpl] split_state(x, idx) theta x(idx.theta); V x(idx.V); pi_g x(idx.pi_g); Ts x(idx.Ts); Tr x(idx.Tr); x_cpl x(idx.cpl); end3.2 统一求解法与解耦求解法的实现选择之前说了我选统一求解法。实际编码时核心就是构建统一的残差向量R(x)和雅可比矩阵J。残差向量按顺序排列电网有功偏差、电网无功偏差、气网流量偏差、热网水力偏差、热网热力偏差、耦合设备方程偏差。状态向量也按同样顺序排列。统一求解法的雅可比矩阵维度是所有变量数之和。对一个中小型系统比如电力5节点、气网4节点、热网3节点、CHP和热泵各1台状态变量数大概在24个左右雅可比矩阵24×24数值差分求解非常快。如果扩展到几十节点就应该用解析雅可比或者稀疏矩阵。我用Matlab的sparse函数组装J配合左除\求解线性方程组速度和内存都很好。dx -J_full \ res; x x dx;为了提升迭代稳定性我加入了阻尼策略如果新的残差范数没有下降就减小步长这有点像Levenberg-Marquardt的思路。实际测试中阻尼因子取0.5能有效避免震荡。解耦迭代法也有它的应用场景比如系统规模很大气网和热网都有成熟的商业计算软件外接替代统一求解是合理的。但从学术复现和代码透明性角度统一求解法更友好这也是我在博客里推荐统一法的原因。3.3 关键函数与核心代码片段讲解残差函数residual(x)是整个程序的心脏。以我搭建的简易算例为例简要展示代码结构。function res residual(x) [theta, V, pi_g, Ts, Tr, x_cpl] split_state(x, idx); % 1. 电网残差 res_grid grid_residual(theta, V, data.grid); % 2. 气网残差 res_gas gas_residual(pi_g, data.gas); % 3. 热网水力残差 res_heat_hyd heat_hydraulic_residual(data.heat); % 4. 热网站热力残差 res_heat_therm heat_thermal_residual(Ts, Tr, data.heat); % 5. 耦合设备残差 res_cpl coupling_residual(x_cpl, theta, V, pi_g, Ts, Tr, data.coupling); res [res_grid; res_gas; res_heat_hyd; res_heat_therm; res_cpl]; end耦合设备残差需要把电、气、热变量都串起来。比如CHP机组的残差方程为function res_cpl coupling_residual(x_cpl, V, theta, pi_g, Ts, Tr, cpl) % 提取CHP变量 P_chp x_cpl(1); Q_heat x_cpl(2); F_gas x_cpl(3); % CHP耦合方程 res_cpl(1) F_gas - P_chp / cpl.chp.eta_e; res_cpl(2) Q_heat - cpl.chp.eta_h * F_gas * cpl.chp.LHV_ratio; res_cpl(3) node_gas_injection(1) F_gas - gas_load(1); % 气网节点1的流量平衡计入F_gas end这个示例省略了很多细节但表达了核心逻辑。雅可比矩阵可以用jacobianest之类的数值差分函数自动生成或者手动写解析式。数值差分的优点是代码量小、不易出错缺点是步长选择不当会降低精度。我建议第一步先用数值差分跑通再考虑将关键方程替换为解析导数。3.4 算例设置与结果可视化为了验证算法我构造了一个小型区域综合能源系统电网采用5节点系统气网4节点热网3节点。电网节点1接平衡电源节点2接CHP机组节点4接热泵。气网节点1为气源连接到节点2的负荷和CHP的耗气。热网节点1为CHP热源节点2为热泵热源节点3为热负荷。给定负荷条件后Matlab程序在初始值接近真实解时5次迭代内收敛到1e-8精度。如果不调整初值直接采用电网标幺电压1.0、气网统一压力3 MPa、热网统一温度75℃作为初值迭代次数大约8次。整体运行时间不到0.2秒说明统一求解法在小规模系统中性能完全没问题。结果可视化方面我会输出四类图电网节点电压幅值和相角分布图、气网节点压力分布图、热网供回水温度分布图、耦合设备功率流向图。Matlab的bar和stairs画压力分布很直观quiver可以画能流方向但实际写代码时我只会输出数值表格让用户自行决定绘图方式。4. 常见问题与调试经验实录4.1 初值设置与收敛性调试我遇到过最多的收敛性问题是气网初值太差。曾有一段时间我把所有气网节点压力初值设为2 MPa但平衡节点实际压力是4.2 MPa结果牛顿法前两步的雅可比矩阵接近奇异迭代直接发散。后来我加了一段“预置初值”逻辑气网非气源节点压力按气源压力乘0.9初始化热网供水温度统一设为85℃回水温度统一设为45℃电网电压标幺值设为1.0。这个初值策略在多个算例里都表现稳定。如果遇到发散可以尝试先用小负荷系数运行。例如将所有负荷乘以0.5求解完成后再把负荷系数提高到0.75用上一个解作为初值继续求解直到全负荷。这种“逐步加载”策略在处理强非线性系统时非常有效。4.2 耦合设备雅可比矩阵的数值处理数值雅可比虽然方便但步长h的选取很敏感。如果h取得太小差分时出现数值噪声取得太大斜率的线性近似误差也会偏大。我的经验是对不同量纲的变量分开设步长电压和温度用1e-6压力和相角用1e-5。更稳妥的做法是解析求导。比如CHP方程的导数非常简单F对P的导数是1/η_eQ对F的导数是η_h·LHV_ratio。电网潮流对电压和角度的导数公式在教科书里都有照着写成稀疏矩阵即可。两者结合计算速度快一倍以上。4.3 结果验证与误差分析技巧代码写完之后必须验证正确性。我习惯用“退化测试”来检查将耦合设备出力全部设置为0此时系统退化为独立电网、气网、热网。分别用传统的单一网络潮流程序计算结果然后和统一求解结果对比误差应在1e-6级别。如果偏差大十有八九是耦合方程或者节点连接关系写错了。还可以做能量守恒检验总输入能量天然气功率加上网侧功率应该等于总负荷功率加上各网络损耗。热网损耗包括管道散热气网损耗在绝热稳态下为零电网损耗就是线路有功损耗。这个测试每次计算完都做能快速暴露符号方向错误。以下是我项目中最常遇到的几个问题汇总问题现象可能原因解决办法网格雅可比奇异气网压力初值相差过大预置初值采用逐步加载迭代发散但不报错残差方程量纲不一致统一单位避免压力用MPa而流量用kg/sCHP耦合不收敛热出力和电出力变量耦合顺序错误检查耦合方程确认变量排列顺序结果独立网络对不上符号约定不一致统一采用“注入为正”再比对教科书公式热网温度出现负值混合节点温度公式写错检查质量流量方向确认分母不为零5. 扩展方向与个人实操体会5.1 从稳态能流走向优化调度能流计算稳态求解器是综合能源系统优化调度的基石。我之前在后面做优化时把能流计算封装成一个函数用遗传算法搜索各设备出力每个个体通过调用能流求解器来校验可行性和运行代价。这样优化算法不用关心内部网络约束只需要给定设备出力组合就能得到系统是否越限、损耗有多大等指标。如果你后续要做电-气-热联合优化建议把能流计算模块做得“无状态”化输入设备出力和负荷输出节点状态和损耗。这样各种优化算法都可以直接调用不需要重复建模。5.2 我踩过的坑与个人建议最后分享几个实际教训。第一Matlab代码里不要用大量全局变量尽量用结构体传递数据否则调试的时候一个变量被覆盖根本找不出来。第二一开始不要追求大算例先搭一个3节点电网3节点气网2节点热网的最小系统跑通后再扩规模。第三所有物理量统一用国际单位制气网压力用Pa流量用kg/s这个细节能省去大量排错时间。我现在写能流代码的习惯是每个算例单独放在一个m文件中文件开头列出所有假设条件和单位文件结尾输出关键节点状态。这样三个月后回来看还能快速捡起来。希望这篇记录能帮你避开我踩过的那些坑。如果你也是做综合能源方向的欢迎按这个框架动手试一遍亲手把收敛过程跑出来才算真正把电气热能流计算吃透。