前阵子给一个园区做多能互补方案我先用常规潮流程序把电力网络算了结果发现燃气锅炉和CHP的出力根本定不下来——热负荷一变电网机组的出力就得跟着变单独算电网等于在打移动靶。后来狠下心把电网、气网、热网放到同一个Matlab框架里联立求解才算是把整个区域综合能源系统的稳态运行点一次算清楚。这篇内容就是围绕“计及多能耦合的区域综合能源系统电气热能流计算”这个题目展开的我会把数学模型、耦合元件建模、Matlab代码实现、算例验证和调试经验全串起来讲。适合正在做多能互补、园区综合能源规划、微能网运行优化的学生和工程师参考尤其是手头有项目但不知道如何把三个网络方程统一到一套代码里的朋友。1. 为什么电气热联立求解比想象中复杂得多1.1 多能耦合的本质能量转换设备把三个网络拴在一起所谓多能耦合本质就是电网、气网和热网之间不再彼此独立而是通过CHP机组、燃气锅炉、电锅炉、热泵这些能量转换设备产生了强关联。举个例子CHP机组一端消耗天然气另一端同时输出电能和热能。在热负荷较高的季节CHP被迫多发一些热电出力也会跟着提高如果电网侧此时并不缺电多出来的电就会从平衡节点倒送或者改变其他发电机的出力。这些连锁反应只有把电、气、热三个网络放在同一个方程体系里解才能得到自洽的结果。换句话说耦合设备就是“桥梁”它们把电网节点方程、气网节点方程、热网节点方程中的若干项串了起来。你在代码里只要漏掉任何一侧的耦合项计算结果就会在能量平衡上出现明显缺口而且你很难凭直觉判断缺口来自哪里。1.2 三个网络的数学差异不是改个符号就能复用很多人第一反应是电网潮流都能算了照葫芦画瓢把电阻改成管道阻力不就能算热网气网了吗实测下来完全不是这么回事。三个网络的状态变量和方程性质差异非常大我先用表格梳理一下网络状态变量核心方程最容易出问题的点电网电压幅值 V、相角 θ节点复功率平衡方程功率与电压乘积项带来的强非线性气网节点气压 p或气压平方 π节点流量平衡、管道Weymouth方程流量与压差平方根相关低压差附近导数突变热网节点供水温度 T_s、回水温度 T_r、管道流量 m水力方程流量平衡、回路压降、热力方程热功率、温降压降与流量平方呈非线性温度混合需要按流量加权电网的复功率方程是电压、相角的二次非线性气网的Weymouth方程直接跟压差的平方根打交道热网更特殊它分为水力工况和热力工况两层水力方程里的“压降与流量平方成正比”带有明显的绝对值非线性。你说它们难倒不是公式本身有多复杂而是这些非线性特性在牛顿法的雅可比矩阵里表现得很不一样初值给不好、符号写错、单位不一致都会导致迭代发散。1.3 我的方案统一求解法而不是顺序迭代处理多能潮流业内主要有两条技术路线。顺序迭代法先假定耦合设备的某一侧出力比如 CHP 电出力固定算出电网潮流再拿着电出力去算热网和气网得到新的热/气侧状态后返回修正 CHP 参数如此反复迭代。统一求解法把所有网络的方程和耦合约束合并成一个大型非线性方程组用扩展的牛顿-拉夫逊法一次性求解所有变量。顺序迭代法写起来直观但碰到强耦合场景容易在几次迭代之间来回震荡尤其当热负荷和电负荷同时波动时CHP 的出力在两个解附近反复跳收敛非常痛苦。我在实际项目里最终选了统一求解法虽然装配矩阵的代码量更大但鲁棒性明显好一截后续如果要加储能、加氢能设备也只需要在变量向量和雅可比矩阵里“加行加列”扩展性可控。2. 电气热网络的数学模型先把方程组列对2.1 电网方程极坐标下的经典潮流电网部分沿用常规潮流计算的极坐标牛顿-拉夫逊法。对节点 i有功和无功平衡方程写作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。节点类型仍然是三种平衡节点提供相位基准PV节点给定有功和电压幅值PQ节点给定有功无功负荷。在多能系统里CHP 机组的电气侧通常建模为 PQ 或 PV 节点这取决于控制方式但它的有功出力不是预先给定的常数而是同时被热负荷和气网侧的天然气供给约束所以在统一求解时这个“给定值”会在联立方程的过程中被反复修正。2.2 气网方程Weymouth管道模型和压缩机气网稳态计算里有两条核心约束。第一条是节点流量平衡对任意气网节点进气源的注入流量加上从相邻管道流入的流量等于该节点的气负荷加上泄漏损失一般忽略泄漏。写成矩阵形式就是 A * f L_source - L_load 0A 是节点-管道关联矩阵f 是管道流量向量。第二条是管道流量方程工程上最常用 Weymouth 公式f_pq K_pq * sign(π_p - π_q) * sqrt(abs(π_p^2 - π_q^2))这里的 π 是节点气压的平方K_pq 是管道参数跟管径、长度、摩阻系数、气体性质有关。压缩机在气网里很关键它负责升压典型建模方式有定出口压力、定升压比、定耗气量三种。我用的 Matlaab 代码里默认支持定升压比模式即 p_out CR * p_in同时压缩机自身要消耗一部分天然气可以作为该节点的附加负荷进入流量平衡方程。如果你的实际系统输气距离不长、压力变化不大可以把压缩机先当成理想元件跳过但后期做高压输气网络时一定要补回来。2.3 热网方程水力工况和热力工况分开建模热网比电网麻烦在它有两层方程。水力工况层每个节点的流量守恒方程加上每个回路的压力平衡方程。管道压降方程为h K_h * m * abs(m)其中 m 是管道质量流量K_h 是管道的阻力特性系数。正因为有 abs(m) 这一项方程在流量变号附近不可导牛顿法收敛时会特别敏感。热力工况层忽略动态传热稳态热功率方程为φ C_p * m * (T_supply - T_return)管道沿程温降采用指数衰减模型T_end T_a (T_start - T_a) * exp(-λ * L / (C_p * m))多个支路汇入同一个节点时混合温度按流量加权平均这一点在代码里特别容易写错。我刚开始就是直接取了算术平均值结果热负荷大的一侧温度算高整体能量平衡始终对不上。3. 耦合元件建模CHP、燃气锅炉、电锅炉怎么进入方程组3.1 CHP机组的热电联产约束CHP 是电气热三网耦合最核心的设备。最简单的稳态模型是电热比线性约束H_chp c_m * P_chpc_m 是热电比抽凝式机组通常在 0.8~1.5 之间背压式可能更高。CHP 消耗的天然气流量为G_chp P_chp / (η_e * LHV)或者更精细一些把电效率和热效率分开折算燃料消耗G_chp P_chp / (η_e * LHV) H_chp / (η_h * LHV)在统一求解的残差方程里CHP 会同时出现在三个位置电网节点注入项、气网节点负荷项、热网节点热源项。换句话说你每给 CHP 分配一个变量编号雅可比矩阵里就有三块对应区域需要填入非零偏导数这是“多能耦合”在数学上的直观体现。3.2 燃气锅炉与电锅炉相对简单的两网耦合燃气锅炉消耗天然气、输出热水关联气网和热网G_gb H_gb / (η_gb * LHV)电锅炉消耗电能、输出热能关联电网和热网P_eb H_eb / η_eb热泵也是电转热但用的是逆卡诺循环耗电量与供热量之间用 COP 联系P_hp H_hp / COP为了把这几种设备在代码里统一管理我维护了一张耦合设备映射表耦合设备电网侧表现气网侧表现热网侧表现CHP机组注入有功 P_chp消耗燃气 G_chp注入热功率 H_chp燃气锅炉无直接关联消耗燃气 G_gb注入热功率 H_gb电锅炉消耗有功 P_eb无直接关联注入热功率 H_eb热泵消耗有功 P_hp无直接关联注入热功率 H_hp这张表的价值在于写代码前先想清楚每个设备在三个网络的“支点”在哪里后面组装配残差和雅可比时就不会漏项。耦合元件的变量我建议在一开始就作为全局变量放进总变量向量 x 里而不是在每次迭代后单独计算。这样做的好处是让牛顿法的搜索方向同时包含所有网络的可行修正量避免子网络之间出现“互相猜”的现象。4. Matlab代码实现从零搭建计算主循环4.1 总变量编号与数据结构统一求解的第一步是定义总变量向量。我把所有待求状态量按顺序拼接成一个列向量x [ 电网相角θ; 电网PQ节点电压幅值V; 气网节点气压π; 热网节点供水温度T_s; 热网节点回水温度T_r; 热网管道流量m; 耦合设备内部变量 ]每个变量都要分配一个全局索引。我习惯用结构体保存三个网络的原始数据data.elec.bus struct(id,{}, type,{}, V0,{}, th0,{}, Pd,{}, Qd,{}); data.elec.branch struct(fbus,{}, tbus,{}, r,{}, x,{}, b,{}); data.gas.node struct(id,{}, type,{}, p0,{}, load,{}, source,{}); data.gas.pipe struct(fnode,{}, tnode,{}, K,{}, comp,{}); data.heat.node struct(id,{}, type,{}, Ts0,{}, Tr0,{}, phi,{}, flow,{}); data.heat.pipe struct(fnode,{}, tnode,{}, Kh,{}, L,{}, Ta,{});不建议把变量零散地放在一堆全局变量里那样后期改网络规模时很容易索引错乱。用结构体统一挂载数据再用一个索引映射函数把“网络类型节点编号”转换成主变量 x 中的位置是维护这套代码最划算的方式。4.2 残差函数每行残差对应一个物理方程主计算循环围绕“残差函数”和“雅可比矩阵”两个核心展开。残差函数 F(x) 按顺序返回所有方程的不平衡量。电网部分残差是节点有功无功不平衡量气网部分残差是节点流量不平衡量热网部分残差包括水力残差节点流量平衡、回路压降和热力残差热功率方程、管道温降方程。管道流量的核心函数可以这样写function f weymouth_flow(K, pi_k, pi_m) d pi_k^2 - pi_m^2; if abs(d) 1e-12 f 0; else f K * sign(d) * sqrt(abs(d)); end end关键是把 sign() 和 sqrt(abs()) 用对否则当管道两端气压相等时流量方向会在迭代中反复乱跳导致整体发散。热力节点混合温度的计算也放一个小函数function T_mix mix_temperature(m_in, T_in) % m_in: 流入节点的各支路流量向量 % T_in: 对应支路末端温度向量 T_mix sum(m_in .* T_in) / sum(m_in); end这个函数看起来简单但它对应的物理含义是供暖管道多路热水汇合时混合温度等于流量加权平均而不是算术平均。4.3 雅可比矩阵组装先数值差分验证再上解析式统一求解法的装配工作量主要集中在雅可比矩阵 J 上。J 中的元素是每个残差方程对每个变量的偏导数。三个网络内部有各自的非零块耦合设备又在三块区域之间产生耦合非零块。组装方式我推荐分三步走。第一步用数值差分验证逻辑eps0 1e-6; J(:, j) (F(x eps0*ej) - F(x - eps0*ej)) / (2*eps0);先对一个小算例做数值差分雅可比如果牛顿迭代都能收敛说明残差方程本身逻辑是对的。第二步再写解析雅可比替换提升计算速度。不要一上来就写解析式那样很难定位到底是残差写错了还是求导求错了。第三步用稀疏矩阵存储J sparse(I_idx, J_idx, val_list, n_var, n_var);三个网络拼起来之后变量数量很容易超过几百个稠密矩阵求逆在 matlab 里虽然也能跑但每迭代一次的时间会随变量数快速增长尤其是做参数扫描或运行优化时速度差距非常明显。4.4 牛顿迭代主循环主循环可以压得很精简x x0; for k 1:maxIter F compute_residual(x, data); J assemble_jacobian(x, data); dx -(J \ F); x x dx; if norm(F, inf) tol norm(dx, inf) tol break; end end这里我用了双重收敛判据一边看残差范数一边看变量增量范数。只盯变量增量会出现一种假象——步长小不代表方程满足尤其是在热网这种温度变量数量级很小的场景里ΔT1e-6 看起来已经很“收敛”了但对应的热功率不平衡可能还很大。5. 算例验证一个小型区域系统跑通全流程5.1 算例拓扑与参数为了验证代码我搭了一个小型区域综合能源系统电网3个节点节点1为平衡节点节点2接入CHP机组节点3接电锅炉和常规电负荷额定容量100 kVA基准电压10 kV气网4个节点节点1为气源节点2给CHP供气节点3给燃气锅炉供气节点4带常规气负荷热网4个节点节点1连接CHP作为热源1节点2连接燃气锅炉作为热源2节点3和4分别带热负荷100 kW和120 kW热力管网包含3条管道。这部分参数我在建模时故意选得很接近园区实际情况热网供水温度初值80℃回水初值50℃电锅炉额定热出力60 kW热泵COP取3.2CHP热电比取0.93。5.2 迭代过程与收敛性观察以“平启动”初值出发统一求解法在本算例中的迭代过程大致如下迭代次数最大残差范数 max|F|最大变量增量 max|Δx|16.3e-013.6e-0138.9e-021.7e-0154.7e-032.3e-0276.2e-058.5e-0493.9e-084.2e-06112.8e-111.9e-08这里我解释一下该数据是在特定初值、特定参数下得到的结果不具备普遍性但能直观说明一个现象统一法前期收敛很快后期进入线性收敛区整体11步以内就压到了浮点精度附近。收敛后的关键运行状态如下电网侧节点1电压1.00∠0°节点2电压0.982∠-2.1°节点3电压0.976∠-3.5°气网侧气源节点1气压1000 kPa节点2约996 kPa节点3约993 kPa热网侧节点1供水温度80.0℃节点3供水温度76.2℃节点4供水温度74.5℃回水温度分别约49.3℃、48.6℃、47.9℃耦合设备运行点CHP电出力300 kW热出力279 kW燃气锅炉热出力150 kW电锅炉热出力60 kW电功率约67 kW热泵热出力45 kW电功率约14 kW。5.3 结果合理性分析收敛之后不能直接宣告胜利还要做校核。我习惯检查三类能量平衡电网有功平衡平衡节点实际出力 CHP电出力 ≈ 电负荷 电锅炉耗电 热泵耗电 网损气网流量平衡气源总注入 ≈ CHP燃气耗量 燃气锅炉燃气耗量 常规气负荷 管损热网热平衡热源总供热量 ≈ 热负荷 管道散热损失。我的Matlab代码在校核后三项平衡的最大相对误差都在1e-6量级以内说明联立解是自洽的。这里补一句如果你在跑完自己的算例后能量不平衡误差很大九成以上是某个耦合设备漏在了某个网络的残差里而不是牛顿法迭代次数不够。6. 踩坑记录与收敛性调优经验6.1 热网初值供水温度必须大于回水温度这句话看着像废话但在我调试初期真的翻过车。我给热网全部节点初值设置为一样的80℃回水也为80℃结果供水回水温差为0热功率方程残差对温度的导数出现异常雅可比矩阵条件数急剧恶化。后来我把初值策略改成供水各节点统一80℃回水各节点统一50℃管道流量初值取一个与热负荷匹配的估计值再让牛顿法自己去修正。这一步改动非常小收敛稳定性却提升了一个档次。6.2 气网符号函数与低压差附近的处理Weymouth方程的绝对值符号是头号天坑。在迭代早期管道两端气压很可能接近此时如果 sign(d) 写得不对或者 sqrt(abs(d)) 的导数没有平滑牛顿步长会猛冲一下然后第二、三迭代步直接发散。我的处理是在计算雅可比矩阵中 Weymouth 对气压的偏导时给 dπ_k^2-π_m^2 加一个很小的死区或平滑让导数不会在原点处变成无穷大。工程上可以理解为微小压差下流量对压差的变化率不可能真的无穷大做一点约束既符合物理常识又能救回迭代稳定。6.3 三网单位不统一导致的病态矩阵电网习惯用标幺值电压1.0左右气网用kPa几百到上千热网温度用℃几十到一百管道流量用 kg/s可能只有个位数。把这些直接塞进同一个变量向量雅可比矩阵非零元素的量级会从1e-2到1e3跨度巨大矩阵条件数非常差迭代后期精度会很尴尬。我的做法是对变量做基准化温度除以100气压除以1000流量除以各自的基准值让所有变量基本落在0.01~10区间内。这一步做完收敛速度明显变快而且不会出现“前半段正常后半段始终差一点达不到容差”的怪现象。6.4 顺序迭代法的教训强耦合时不要死磕我在做第一版代码时其实是先写的顺序迭代法因为觉得简单。结果是只要热负荷和 CHP 热电比稍微敏感一点外层循环就在“CHP多发电→热侧偏热→调低电出力→电网不平衡→提高CHP出力”的循环里反复横跳设定最大迭代次数也救不回来。后来我统计了一下震荡点根本原因在于 CHP 同时被电、气、热三个约束牵制顺序法的每一层都在“临时固定”另外两个网络的变量信息传递存在一步延迟。统一求解法之所以稳健是因为它在一个非线性方程组里同时处理了所有物理约束三类残差共享同一个变量向量各变量在同一个迭代步内同步更新不需要人为设定收敛次序。6.5 建议的调试路径按网络逐层叠加如果一开始就搭完整的三网联立系统出问题之后几乎没法定位。有一个非常实用的调试策略先在纯电网纯热网的小系统上把耦合设备只有电锅炉/热泵的情况跑通再加入气网此时只需要固定燃气锅炉和CHP的气侧消耗把气网作为“跟随网络”接入最后才把CHP的电-热-气三联互动打开做完整统一求解。每加一层网络之前保留一个可运行的存档版本方便对照调试。这套思路帮我节省了大量时间也建议你复现类似项目时采用。一点个人的实操体会最后单独说一个我反复测试中认为最有用的技巧不要迷信“全解析雅可比”。只要条件允许先用数值差分雅可比把模型逻辑验证清楚再根据矩阵稀疏模式用 symbolic 工具辅助求导最后人工校正关键耦合块的偏导数。这样虽然前期麻烦一点但能大幅降低“模型对但代码错”的隐蔽 bug 概率。多能耦合潮流计算在数学上并没有比常规电力潮流高深太多难的是三套方程、三类变量、多种单位体系要在同一个架构里稳定共处。我的体验是只要建模一致性守住了收敛性基本不会成为瓶颈如果你的迭代总是发散先回头检查耦合设备的能量流方向是不是统一了。把功夫下在建模和数据结构上Matlab 这套代码完全可以在几十个节点的园区级规模上稳定运行。希望这篇文章能帮你少走一些我走过的弯路。