先说两句题外的。“区域综合能源系统电气热能流计算”这个名字看着吓人拆开之后做的事情其实很单纯一个网络里同时跑着电、天然气、热三种能量我们要把它们的稳态分布一次性算出来。传统电力潮流算的是电压和相角气网算的是节点压力和管道流量热网算的是温度和质量流量现在要把它们放在同一个方程组里联立求解。我在实际项目里最深的感受是搞明白三个网络各自的模型大概花三成精力搞清楚耦合设备怎么把三个网络“焊”起来再花三成精力剩下四成时间全在跟雅可比矩阵、初值、收敛性这些“具体执行问题”搏斗。这篇博文不打算堆理论而是按我自己写Matlab程序的完整过程来讲从方程建模到代码实现再到调试踩坑尽量讲透。1. 从电力潮流到多能流为什么传统方法在综合能源系统里失效1.1 三个独立程序“摞”不在一起很多刚接触这个方向的人会有一个直觉想法电网用牛拉法算潮流气网用牛顿法算节点压力热网用节点法算温度这不是都有现成程序吗三个程序迭代着互相传数据不就行了这个方法叫“顺序求解法”Sequential Method确实有人这么做但工程上很容易出问题。最典型的情况是电网迭代到第3次需要气网给出天然气供给量气网迭代到第5次需要热网给出热负荷热网又需要电网给出电负荷——三个网络迭代节奏不一致来回传数据会非常混乱。更麻烦的是CHP热电联产机组同时决定发电量、产热量和耗气量你根本说不清这个耦合量应该由哪个网络“说了算”。所以现在的主流做法是“统一求解法”Unified Method也叫联立求解法把电网的电压、相角气网的压力热网的温度、流量全部塞进同一个状态向量里构造一组统一的残差方程用一个大的牛顿-拉夫逊迭代一次性解出来。1.2 三类网络的状态变量“代沟”很大这个“统一”说起来容易做起来第一个坎就是状态变量差异太大。电网状态量是电压幅值V和相角θ量级通常在0.9到1.1之间标幺值气网状态量是节点压力p工程里常用MPa范围可能是0.3到4.0热网状态量是供水温度Ts、回水温度Tr和质量流量m温度在70到110摄氏度左右。这三个量放在同一套矩阵里雅可比矩阵的元素数值就会横跨好几个数量级。比如∂(电网有功不平衡量)/∂(电压相角)可能在1附近而∂(气网流量不平衡量)/∂(管道压力)可能是0.001∂(热网热力不平衡量)/∂(供水温度)可能是几百。如果不做处理迭代过程中的修正量会被大数值项主导小数值项对应的状态变量几乎得不到有效更新收敛性极差。我在第一次写程序时就吃了这个亏电网部分很快收敛气网和热网却一直震荡。后来才意识到这不光是初值问题更是量纲问题。解决办法后面细说核心思想就是“标准化”——把各网络的物理量都折算到基准值附近让雅可比矩阵不会出现极端病态。1.3 耦合设备是“胶水”也是“麻烦源”单网络潮流里每个节点要么是PQ节点、PV节点要么是平衡节点性质很明确。但在综合能源系统里一个CHP节点同时连着电网、气网、热网它既是电网的发电机节点又是气网的负荷节点还是热网的热源节点。关键变量之间的依赖关系是交叉的CHP发电量增加 → 耗气量增加 → 气网节点压力下降CHP产热量取决于发电量热电比→ 热网供水温度变化气网压力变化 → CHP的出力受限 → 电网潮流重新分布这种“交叉耦合”意味着雅可比矩阵里会出现大量非对角块元素。统一求解法的本质就是把这些交叉偏导数全部算出来并与对角块放在一起迭代让每个网络的变化能同时“看见”其他网络的变化。2. 电、气、热三网方程建模状态变量与残差方程的对应关系2.1 电网极坐标潮流方程电网部分用的是最经典极坐标牛顿-拉夫逊潮流模型这点没什么可创新的。唯一需要注意的是在统一求解框架下电网的节点类型可能发生变化CHP节点在有热负荷时不能简单当成PV节点处理——它的发电量受热负荷和天然气供给的约束有功出力不是独立变量。对每个PQ节点和PV节点有功不平衡方程ΔP_i P_g,i - P_load,i - V_i * Σ_j V_j (G_ij * cos(θ_ij) B_ij * sin(θ_ij))对PQ节点外加无功不平衡方程ΔQ_i Q_g,i - Q_load,i - V_i * Σ_j V_j (G_ij * sin(θ_ij) - B_ij * cos(θ_ij))对平衡节点一般选主电网接入点电压幅值和相角固定不参与迭代更新。2.2 天然气网Weymouth方程与节点流量平衡天然气管道稳态模型常用Weymouth方程描述。对连接节点i和j的管道f_ij C_ij * sign(p_i^2 - p_j^2) * sqrt(|p_i^2 - p_j^2|)其中f_ij是管道流量标况下体积流量C_ij是与管道直径、长度、摩擦系数、气体性质有关的常数sign保证流量方向与压差方向一致。每个节点的流量平衡方程Σ_j f_ij f_source,i - f_load,i 0除了管道天然气网络里还有压缩机和阀门。压缩机的作用是给气体增压补偿管道沿程压力损失。压缩机分定压比、定流量、定出口压力等几种控制模式数学模型也不一样。我在代码里用的是最常用的定出口压力模式p_out p_set这样压缩机就相当于在方程里加了一个固定压力约束实现起来比定压比模式容易得多。定压比模式为什么麻烦后面在第6章详细讲。2.3 热网水力与热力方程必须同时求解热网是这三大网络里最容易建模出错的因为它不是一套方程而是水力模型和热力模型两套方程同时成立。水力模型描述质量流量和压力的关系包括两个方程节点质量流量守恒每个节点的流入流量之和等于流出流量之和环路压头守恒任何一个闭合回路的压头降之和为零包含水泵扬程对应状态量是各管道的质量流量m。水力方程本质上是线性的如果不考虑管路阻力非线性的话但阻力损失与流量的平方成正比所以在牛顿法迭代中仍需处理非线性。热力模型描述温度和热功率的关系核心方程有两个管道温度衰减方程考虑热损失T_end T_start * exp(-λ * L / (Cp * m))这个公式的意思是热水在管道里流动沿途向环境散热温度沿管长呈指数衰减。λ是管道热损失系数L是管长Cp是水的比热容。节点温度混合方程T_mix (Σ m_in * T_in) / Σ m_in当不同温度的热水在节点处汇合时混合后的温度按质量流量加权平均。供热站向热网注入的热功率Φ Cp * m * (Ts_supply - Tr_return)在这个模型里常见做法是给负荷节点指定热功率需求Φ然后由热力方程反推节点供/回水温度。这里面有一个容易混乱的点热力方程里质量流量m和温度T是强耦合的因为管道温度衰减的指数项里有m节点热功率方程里也有m。你不能先算水力再算热力必须把m和T放在一起联合迭代。这也是初学者最容易犯的错误之一——很多教材把水力热力分开讲但工程实现时忽略耦合会导致结果完全不对。下面对三类网络的状态量和方程做个总结网络状态变量残差方程对应物理意义电网电压幅值V、电压相角θΔP、ΔQ节点不平衡量节点有功/无功注入平衡气网节点压力p节点流量不平衡量天然气流量注入平衡热网供水温度Ts、回水温度Tr、管道流量m节点热功率、温度混合、环路压降方程组热力学第一定律与质量守恒3. 耦合设备与统一雅可比多能耦合项如何“焊死”在迭代矩阵里3.1 CHP机组综合能源系统的“神经中枢”CHPCombined Heat and Power是区域综合能源系统里最核心的耦合设备因为它同时牵涉三种能量输入天然气输出电和热。我用的CHP模型分两部分。第一部分是发电部分的燃料消耗特性F_gas a * P_elec^2 b * P_elec cF_gas是消耗的天然气体积流量标况下P_elec是发电功率a、b、c是机组特性系数。这个式子说明发电越多耗气量越大而且不是线性关系二次项让机组在低负荷时效率更低。第二部分是热电比模型。把产热量写成发电量的函数H_heat r_heat2power * P_elecr_heat2power是热电比heat-to-power ratio常见燃气内燃机热电比在1.0到2.5之间。注意这个热电比在不同负荷区间可能是变化的工程上可以取分段线性函数但为了程序稳定我更推荐用定热电比模型做第一版跑通了再上复杂模型。CHP在统一迭代里起的作用可以用下图理解它在电网里是发电机有功注入在气网里是负荷天然气消耗在热网里是热源热功率注入。这三个角色的取值不是独立的而是通过上面的两个方程锁死的。体现在雅可比矩阵里就是在电网有功残差对CHP发电量的偏导、气网流量残差对CHP耗气量的偏导、热网热功率残差对CHP产热量的偏导之间建立交叉引用关系。3.2 电锅炉与P2G两种不同方向的双向转换除了CHP常见的耦合设备还有电锅炉电→热、燃气锅炉气→热、电转气P2G电→气和热泵电→热但效率更高。电锅炉的耦合模型最简单Φ_heat η_EB * P_elecη_EB是电锅炉的电热转换效率通常在0.95到0.99之间。电锅炉在系统中常作为CHP的补充热源当CHP产热不足时补足热缺口。P2GPower to Gas是电→气的设备用富余电力电解水制氢再通过甲烷化反应生成天然气主要成分是甲烷。耦合方程F_gas_prod η_P2G * P_elec / HHV_gasHHV_gas是天然气高位热值单位体积的能量η_P2G是整体转换效率。P2G有意思的点在于它是“反向耦合”它把电网的电力需求转变成气网的天然气供应相当于一台“可以随时启停、功率可调的人造气源”。在电气热联立的雅可比矩阵里P2G节点的电网有功残差和气网流量残差互为镜像——电网多消耗一份电气网就多注入一份气。3.3 耦合偏导数的推导以CHP为例统一求解法最关键的一步是正确求出耦合设备涉及的偏导数填入雅可比矩阵的交叉位置。我直接以CHP为例推导一遍。假设CHP节点是电网中的节点i注入有功P_CHP、气网中的节点k消耗天然气F_CHP、热网中的节点j注入热功率H_CHP。电网有功残差对CHP发电量的偏导∂ΔP_i / ∂P_CHP 1气网流量残差对CHP发电量的偏导经由耗气量传递∂Δf_k / ∂P_CHP ∂Δf_k / ∂F_CHP * ∂F_CHP / ∂P_CHP其中∂F_CHP / ∂P_CHP 2*a*P_CHP b热网热功率残差对CHP发电量的偏导经由热电比传递∂ΔΦ_j / ∂P_CHP r_heat2power看到没这里的关键操作是CHP发电量P_CHP作为中间变量同时出现在电网残差、气网残差和热网残差里。统一求解时P_CHP实际上被消去转换为电网节点注入约束变化对气网压力、热网温度的直接偏导关系。这就是“多能耦合项焊进雅可比矩阵”的含义——矩阵里这些交叉元素不是人为拼凑的而是由耦合设备的物理方程通过链式法则自然导出的。4. 牛顿-拉夫逊统一求解分块矩阵装配与Matlab实现细节4.1 状态向量与残差向量的统一构造先把所有要解的状态量按顺序排成一个列向量。我的排序习惯是X [电网的θ; 电网的V; 气网的p; 热网的m; 热网的Ts; 热网的Tr]注意这里没有把CHP发电量作为额外状态变量加入X因为它是中间变量由耦合方程消去了。如果你把P_CHP也塞进X方程数就变了会导致雅可比矩阵不是方阵没法直接求逆。对应的残差向量F构造如下F [电网有功不平衡量ΔP; 电网无功不平衡量ΔQ; 气网节点流量不平衡量; 热网水力残差; 热网热力残差]每个状态量对应一个方程向量维数一致。4.2 雅可比矩阵的分块结构雅可比矩阵J由各子网络和耦合项的偏导数构成。记作分块形式就是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表示电网残差对气网状态量的偏导其他类似。在只有CHP耦合时J_EG的构成逻辑是电网有功残差ΔP对热网供水温度Ts的偏导经由CHP传热过程传递。如果节点i是CHP节点CHP产热量H_heat随发电量P_CHP变化而ΔP中含P_CHP所以∂ΔP_i/∂Ts_j需要用到∂P_CHP/∂H_heat热电比的倒数以及热网方程中H_heat对Ts的关系。这条链路比较绕我建议推导时一步步用链式法则写完整不要跳步。代码实现时我不用符号求导而是直接手动推导偏导公式写函数。原因很简单统一求解的雅可比矩阵规模不大手动推导一次后面所有算例都能复用符号求导在Matlab里虽然可以做但表达式复杂时生成的速度慢而且调试起来不直观。给一个CHP耦合偏导在Matlab中填入的示意代码% 构造CHP耦合偏导 % P_CHP CHP发电量r_h2p 热电比 % 电网侧CHP节点的发电功率随热负荷需求变化 % 实际编程时J矩阵稀疏填充 % 电网残差对CHP发电量的偏导 dF_dPchp 1; % ΔP对P_CHP dQ_dPchp 0; % ΔQ不含P_CHP假设无功给定 % 气网残差对CHP发电量的偏导 dGas_dPchp -(2 * a * P_chp b); % 耗气量F_CHP的导数负号因为残差是流入减流出 % 热网残差对CHP发电量的偏导NN等于 输出热功率对发电量的导数 减去 热网节点热功率对温度的导数乘以温度对发电量的变化 dHeat_dPchp r_h2p; % 产热量对发电量导数但是这里的偏导填入矩阵的位置有讲究CHP节点电网残差出现在电网块的位置它对应的影响会通过“热负荷变化→CHP产热变化→发电量变化→电网注入变化”这条链路传递所以要在J矩阵的电网残差行、热网温度列的位置填入非零元素。4.3 Matlab主迭代循环实现实际写代码时主流程大概是这样的伪代码风格% 初始化 X init_state(); % 设置初值 tol 1e-6; max_iter 30; for k 1:max_iter % 解耦计算各网络残差 [dP, dQ] elec_residual(X); [dGas] gas_residual(X); [dHyd, dTherm] heat_residual(X); F [dP; dQ; dGas; dHyd; dTherm]; % 检查收敛 if norm(F, inf) tol fprintf(收敛于第%d次迭代\n, k); break; end % 计算雅可比矩阵稀疏存储 J jacobian_unified(X); % 求解修正量 dX -J \ F; % 更新状态量 X X dX; % 打印迭代信息 fprintf(迭代次数: %d, 残差范数: %.6e\n, k, norm(F, inf)); end这里最核心的工程技巧是雅可比矩阵一定要用稀疏矩阵存储和计算。如果不加任何处理一个几十节点的综合能源系统雅可比矩阵维度可能在100到200之间直接用满阵存储也花不了太多内存。但一旦系统规模到几百个节点满阵的求逆效率会急剧下降。用Matlab的sparse函数构建稀疏雅可比再用“\”运算求解速度能提升一个数量级。另一个容易出错的地方是联合雅可比矩阵未必是方阵。电网有平衡节点电压相角不迭代气网也要设定源节点压力为定值不迭代热网同样要有参考节点温度作为基准。这三个“松弛条件”缺一不可否则矩阵是奇异的Matlab会直接报错或给出一堆NaN。5. 算例调试初值选取、量纲标准化与收敛过程验证5.1 测试系统结构与参数我调试用的系统是一个小型区域综合能源系统规模大概是13节点电网含平衡节点、6节点气网、5节点热网。耦合设备有一台CHP、一台电锅炉和一个P2G设备。关键参数如下设备参数数值CHP额定发电功率5 MWCHP热电比1.5CHP燃料系数a/b/c0.002 / 0.12 / 0.03电锅炉电转热效率0.97P2G电转气效率0.63天然气高位热值HHV38 MJ/m³5.2 初值策略与迭代收敛过程初值选取对统一求解法来说重要性不亚于方程本身。我用了一段时间才摸出一套靠得住的经验电网采用平启动所有PQ节点电压幅值1.0相角0PV节点电压给定。气网所有非源节点压力初始化为源节点压力的0.8到0.9倍。这个非常关键如果你把气网节点压力初值设为0流量残差会直接给出巨大数值迭代直接崩。我第一次跑程序时气网初值用了0.5 MPa实际应在1.2 MPa附近雅可比矩阵里有大量sqrt项初值偏差太大导致sqrt里出现负数程序直接报复数错误。热网供水温度初值设为系统设计值比如95摄氏度回水温度设为设计值比如65摄氏度。管道流量初值按负荷热功率除以比热容乘以温差估算不要随便给0。耦合设备CHP发电功率初值设为其额定发电量的60%热负荷按热电比反推。P2G和电锅炉的功率初值按其设计功率的50%起。用这套初值典型收敛过程如下表截取部分迭代记录迭代次数电网最大不平衡量气网最大不平衡量热网最大不平衡量总残差范数11.2e-013.4e-024.8e014.9e0122.8e-031.7e-035.2e005.3e0039.1e-058.2e-054.1e-014.1e-0142.3e-063.5e-062.6e-022.6e-0254.2e-089.8e-089.2e-049.2e-0462.9e-093.7e-107.1e-067.1e-0671.1e-102.2e-112.3e-082.3e-08可以看到前两次迭代总残差范数下降很快但第三次以后热网残差下降速度明显慢于电网和气网。这是因为热网方程里温度衰减指数项的非线性更强在所有网络中收敛最慢。我通常以总残差范数小于1e-6作为收敛判据大概7次左右能完成。5.3 量纲标准化的操作细节前面提到的量纲问题具体操作办法有两种。我两种都试过各有利弊。方法一标幺值化。给电网设定基准功率和基准电压给气网设定基准压力和基准流量给热网设定基准温度和基准质量流量。各网络残差方程在标幺值下求解好处是各物理量的数值都在1附近雅可比矩阵数值尺度比较均匀。缺点是耦合设备的方程要跟着换算比如CHP的燃料曲线系数a、b、c都需要按新基准重新标定写程序时容易搞混。方法二物理单位 残差归一化。方程保持原始物理单位但在构造总残差向量时把各网络的残差除以各自的最大量纲基准。比如电网残差除以基准功率100 MW气网残差除以基准流量100 m³/min热网残差除以基准热功率50 MW。这样总残差向量各分量的数量级就比较接近收敛判断也好写。我在常规项目里更偏向方法二因为耦合设备的方程不用动每类方程都保留清晰的物理意义排查错误时直接看残差就能定位问题。只在做科研对比、需要和文献结果对标的时才换用方法一写标幺值版本。6. 实际工程中踩过的坑与处理对策6.1 初值发散热网回水温度初值的坑我最开始跑统一求解时程序在第三、四次迭代就出现NaN。排查了很久发现问题是回水温度初值给得太随意。热网的供水温度和回水温度必须满足负荷侧的热平衡Φ_load Cp * m * (Ts - Tr)如果初值里Tr给得过大接近Ts甚至超过Ts热功率方程里出现负值对应管道温度衰减指数和高次项就会爆炸。后来我的做法是先用设计工况求出一组理性的回水温度再用它做初值。再不行就干脆把回水温度初值设置为供水温度减去一个固定经验差值比如25到30摄氏度保证初始热功率为正。6.2 压缩机雅可比奇异从定压比到定出口压力气网压缩机是最容易让雅可比矩阵翻车的设备。定压比模式的意思是p_out / p_in R_cR_c是固定压缩比。这个约束的雅可比元素设计起来很麻烦特别是在多级压缩、复杂环网结构里压缩比约束会让雅可比矩阵出现近线性相关的行因为多台压缩机串联时出口压力等于入口压力乘以一系列压缩比的乘积方程组间存在强相关Matlab解线性方程时容易报矩阵奇异。我的解决办法是把定压比压缩机改成定出口压力模式。给定每台压缩机的出口压力设定值这样方程就变成一个简单的固定压力约束雅可比矩阵结构更稳定数学上不太会出现奇异性。代价是压缩比会随入口压力波动而变化当进气压力下降时压缩机可能超过最大压缩比限制。解决方式是潮流计算的每次迭代后检查压缩机压比是否超限超限则把该压缩机切换为定压比模式重新计算。这种“模式切换”机制虽然要加一点判断逻辑但程序稳定性提升非常明显。6.3 稀疏矩阵与非零元设置别忽视数值精度雅可比矩阵里大量元素其实数值非常小理论上应该为零因为电网残差对气网压力的偏导本来就是0只有通过耦合设备传递才非零。如果你直接用满阵存储这些小数值会被保留在矩阵里影响求解精度。更稳妥的做法是只在耦合设备对应的行索引和列索引位置填入非零交叉偏导其余位置全部保持为零然后用Matlab的sparse函数构建矩阵。我还遇到过一个隐蔽问题Matlab的sparse矩阵默认按双精度存储但当你手动填入很多理论值为零但数值上因为浮点误差而不严格为零的元素时矩阵会变得比预期更稠密求解效率下降。所以我在代码里加了一个判断数值绝对值小于某个阈值比如1e-10就截断为零。% 雅可比矩阵填充时跳过数值接近零的元素 tol_zero 1e-10; if abs(val) tol_zero J(row, col) J(row, col) val; end这个小技巧看起来不起眼但对大型系统来说可以显著减少稀疏矩阵的非零元数量让“\”求解速度快三倍以上。6.4 设备约束不能硬塞进潮流方程最后要提醒一个理念上的问题很多设备约束不适合直接写进潮流方程里。比如CHP发电量的上下限、热泵的启动最小负荷、P2G的爬坡速率这些约束本质上是运行优化问题里的不等式约束如果你一股脑塞进潮流方程方程数大于变量数雅可比矩阵就不再是方阵根本无法求解。我目前的处理方法是先按解潮流的方式算出一个“不受限”的稳态解再检查所有耦合设备是否越限。如果有设备越限就把该设备的有功出力或产热量固定到边界值把它对应的方程从残差方程组里移除重新解一次。这套“先算后罚”的策略在工程上稳妥虽然迭代次数可能增多但每一步的数学性质都是良定义的。写在最后我之前做这个方向时最大的体会是多能流计算本身并不复杂难点全在“耦合”二字——耦合设备怎么建模、耦合偏导怎么推、耦合对收敛性的影响怎么压住。Matlab的优势在于矩阵操作非常方便稀疏矩阵、左除求解都是内置的写起统一求解法来比手写C省太多事。如果你要在这个基础上扩展我个人觉得有几个方向值得做一是把热网的动态特性管道热惯性加进去变成“准动态多能流”二是引入不确定性用场景法或区间法描述风电、光伏和负荷波动三是把多能流计算嵌进优化模型里做最优调度那是另一个更有工程价值的大坑。先跑通统一的、稳态的、确定性的版本后面一切都好说。