在去年做园区综合能源规划评审的时候甲方问了我一个挺扎心的问题“你们给出的电负荷、热负荷、气负荷预测是三个团队各算各的还是放在同一个框架里算的”我承认当时多少有点心虚因为行业里确实习惯电力潮流归电力、热网水力归热网、天然气管网归燃气管网最后拼成一张系统图完事。但近几年不行了——燃气轮机同时供电气和热电锅炉在低谷时段批量启动P2G设备甚至能反向影响天然气管网的压力。CHP机组的发电计划直接决定热网有多少热量、气网抽多少气电负荷波动会通过电锅炉传导到热网天然气管网压力波动又能限制燃机出力上限。面对这种强耦合场景还在用解耦方式做区域综合能源系统的电气热能流计算结果可信度真的要打问号。这篇文章我就从物理模型、耦合设备建模、Matlab求解器实现、算例验证到实际调试踩坑把这一整套流程完整讲一遍适合正在做综合能源规划、运行优化或者课题研究需要自己手写能流计算代码的同学。1. 为什么电气热能流必须放在同一张图里迭代求解1.1 “各算各的”到底会差多少先举一个很典型的场景。某个工业园区里有一台额定电功率3 MW的CHP机组热电比1.2也就是说发3 MW电的同时能回收3.6 MW热量。如果按传统的解耦思路电力潮流计算时会把CHP当成一个固定出力为3 MW的发电机节点热力计算时把热负荷5 MW当成固定节点负荷天然气计算时把耗气量当成固定值。问题是实际运行中CHP大概率是“以热定电”——热网需要多少热CHP就发多少电。热负荷变化时注入配电网的电功率跟着变配电网的电压分布也随之改变。反过来如果电网侧电压偏低、调度要求CHP降出力热网瞬时缺的热量就得靠燃气锅炉顶上燃气锅炉的耗气量又叠加到气网负荷上。这几条链路互相咬合没有任何一个子系统可以独立求得准确解。我做过一个不算复杂的测试一个IEEE 33节点配电网、6节点气网、6节点热网组成的小型综合能源系统里面只有一台CHP和一台电锅炉。解耦计算和耦合计算结果对比下来配电网节点电压最大偏差接近1%气网节点压力最大偏差超过4%热网节点温度偏差约0.5℃。如果不考虑P2G这个偏差还不算离谱一旦加上P2G气网压力偏差可以直接拉到15%以上甚至出现方向性错误——比如解耦算出来某节点压力正常耦合算出来该节点已经跌破最低供气压力。1.2 耦合的数学本质从数学上看三个子系统的稳态模型可以统一写成一组非线性代数方程F_elec(x_elec, x_gas, x_heat) 0 F_gas(x_gas, x_elec, x_heat) 0 F_heat(x_heat, x_elec, x_gas) 0解耦计算等价于把x_gas和x_heat固定在某组常数上去解F_elec把x_elec和x_heat固定去解F_gas依次类推。这在数学上等同于完全忽略了交叉偏导项∂F_elec/∂x_gas、∂F_gas/∂x_elec这些非对角块。当耦合设备渗透率低的时候忽略交叉偏导带来的误差也许在工程接受范围内但耦合度一旦上去这些交叉偏导不仅数值可观而且随着运行点剧烈变化解耦结果就完全不可控了。所以真正要做“计及多能耦合”的能流计算技术上就两条路一是把所有方程联立起来做统一牛顿法求解二是分解迭代但反复交换耦合变量直到收敛。两者没有绝对的优劣实践里怎么选我后面会详细讲。2. 从物理方程到可计算的残差形式2.1 电力网络潮流方程电力部分继续沿用经典极坐标牛顿-拉夫逊潮流模型。对每个PQ节点写有功和无功两个功率平衡残差ΔP_i P_i^sp - V_i * Σ(V_j * (G_ij*cosθ_ij B_ij*sinθ_ij)) 0 ΔQ_i Q_i^sp - V_i * Σ(V_j * (G_ij*sinθ_ij - B_ij*cosθ_ij)) 0PV节点只写有功残差平衡节点什么都不写状态变量是除平衡节点外所有节点的电压相角和PQ节点的电压幅值。这块大家都很熟悉不多展开。需要特别注意的是节点注入功率P_i^sp里面如果接入了CHP或者电锅炉对应的功率项不能预先写成常数。比如CHP节点P_i^sp P_load - P_chp而P_chp本身是耦合设备模型算出来的、随迭代变化的值。忽略这一点代码写得再漂亮本质上还是解耦的。2.2 天然气网络稳态流量方程稳态天然气网络用节点气压作为状态变量管段流量用Weymouth方程描述。中高压管网中管段流量和两端气压的关系近似写成f_ij sign(p_i - p_j) * C_ij * sqrt(|p_i² - p_j²|)C_ij是管段综合系数包含管径、长度、温度、压缩因子等参数。注意这里必须带sign符号项否则管内流量方向固定死了实际管网中气流方向是随压力分布变化的忽略sign会导致负值开根号直接报错。每个节点的流量平衡残差写成ΔF_i Σf_in - Σf_out F_source - F_load - F_device 0其中F_device就是从该节点取气的耦合设备耗气量比如CHP、燃气锅炉、P2G的耗气P2G是产气符号相反。天然气负荷里一部分是固定负荷这部分提前给定另一部分是随耦合设备出力变化的“动态负荷”必须在每次迭代里更新。压缩机模型这里只说一种常用的简化处理给定升压比或者给定出口压力压缩机的自耗气量按经验公式估算。如果一开始就上精细的动态压缩机模型代码复杂度和调试难度会陡增而且对稳态能流结果的影响通常有限建议第一版代码先做简化。2.3 热力网络的水力-热力联合方程热力网络是所有子系统里最容易写错、也最容易数值翻车的一块。它本身包含两个层次水力方程和热力方程两者通过管道流量耦合在一起。水力方程一是节点流量连续性A * m qA是节点-管道关联矩阵m是管道流量向量q是节点净注入流量二是管道压降方程Δp_drop K * m * |m|K是管道阻力系数流量与压降呈平方关系。在供热管网里循环泵提供的压头平衡所有管道压降。热力方程这边首先每个热负荷节点的热功率与供回水温差有关Φ_load C_p * m_node * (T_s - T_o)m_node是流经该负荷节点的流量T_s和T_o分别是供水温度和回水温度。然后是管道沿程温降考虑环境散热时管道末端温度不等于首端温度T_out (T_in - T_a) * exp(-λ*L / (C_p*m)) T_a最后是节点混合温度方程多条支路的热水在节点汇合后温度是流量加权平均Σ(m_in * T_in) T_mix * Σ(m_in)热力网络里我最想强调的一个实操点状态变量不应该把“每条管道流量、每个节点供回水温度”全塞进去硬解因为节点混合温度实际上是由代数关系直接决定的。更稳的做法是把节点混合前后两段变量分离通过消元把混合温度表达成支路温度的函数只保留独立变量进雅可比矩阵。直接全塞进去雅可比矩阵里会出现大量线性相关的行轻则收敛慢重则矩阵奇异直接报错。2.4 耦合设备CHP、电锅炉、燃气锅炉与P2G耦合设备的建模决定了整个系统的耦合强度也是这篇代码的灵魂。我的做法是统一成“功率-能耗转换器”的视角来看CHP机组是最经典的耦合元件。简化的背压式机组模型用热电比描述Φ_recover R_h * P_chp F_gas P_chp / (η_e * LHV)P_chp是发电机出力Φ_recover是回收的热功率R_h是热电比η_e是发电效率LHV是天然气低位热值。这里的热电比在中小型机组上通常在1.01.5之间选型时拿厂家样本数据即可。注意CHP的“以热定电”和“以电定热”两种运行方式在代码里对应不同的变量约束实现时要留一个控制模式开关。电锅炉模型最简单电功率转热功率Φ_eb η_eb * P_eb效率通常在0.95以上基本可以当成纯电阻加热处理。燃气锅炉也一样天然气化学能转热能Φ_gb η_gb * F_gas * LHVP2G要稍微费点心思。电解水制氢的效率决定了电转氢的量F_H2 η_p2g * P_p2g / LHV_H2但氢气能不能直接注入天然气网络取决于掺氢比例约束。实际工程中掺氢比例通常限制在5%20%视管网材质和终端燃气具而定如果代码里直接把氢气产量全部注入气网算出来的气压可能明显偏离真实约束。稳妥做法是加一个简单的上限约束比如注入气网节点的氢气等效流量不超过该节点总流量的10%。3. Matlab代码总体架构统一牛顿法为主、分解法兜底3.1 两种求解框架的取舍统一求解法把电气热三个子系统的状态变量堆成一个长向量所有方程堆成一个长残差向量直接做一次牛顿迭代x [Va; Vm; p_gas; m_pipe; Ts; Tr] F [F_elec; F_gas; F_hydraulic; F_thermal; F_coupling] x_{k1} x_k - J(x_k) \ F(x_k)这个方案的优点是严格的二阶收敛交叉耦合项完整保留是我在强耦合场景下的首选。缺点是编程工程量确实大尤其是雅可比矩阵的组装四个子系统加耦合块的偏导要全手动推导和编码初期调试周期不短。分解法则是按“电力→天然气→热力”顺序轮流迭代每个子系统用成熟求解器子系统之间互相交换耦合变量直到整体收敛。模块化好、代码复用率高但强耦合时经常出现迭代震荡甚至发散要在各子系统之间加松弛因子x_new x_old α*(x_computed - x_old)α取0.50.8比较常见。实测下来分解法在耦合度中等单台CHP的场景还能勉强工作一旦P2G和电锅炉都上来收敛速度会急剧恶化迭代次数从几次飙到几十次。我自己的工程建议是两阶段策略先用分解法迭代58轮得到一个合理的中间初值然后切换统一牛顿法精解。这个做法既避免了统一法冷启动初值太差导致发散的问题又能在最后阶段拿到严格耦合解。代码里控制这个切换非常简单一个if判断即可。3.2 代码模块划分与数据结构我自己写这类代码的习惯是用struct组织全系统数据不用class。struct的好处是轻量、透明课题组里其他人接手时一眼就能看明白字段含义改起来也无压力。整体模块划分如下01_load_case.m 读取并检查网络数据 02_build_index.m 生成状态变量全局索引 03_init_guess.m 设定迭代初值 04_ies_residual.m 计算残差向量 05_ies_jacobian.m 计算雅可比矩阵 06_newton_solver.m 牛顿迭代主循环 07_post_process.m 结果统计与可视化 run_ies_pf.m 主脚本一键运行这种模块化划分的核心思路是让“残差计算”和“雅可比计算”完全分离后续想升级模型比如换更精细的压缩机模型只需要改04和05两个文件不用动主循环。数据结构大体长这样% ---- 电气热综合能源系统数据定义 ---- ies struct(); % 电力系统 ies.elec.baseMVA 10; ies.elec.bus [bus_i, type, Pd, Qd, ...]; % 节点数据 ies.elec.branch [fbus, tbus, r, x, b, ...]; % 支路数据 % 天然气系统 ies.gas.baseP 1.0; % 基准压力 MPa ies.gas.node [node_i, load, source, ...]; % 节点数据 ies.gas.pipe [fnode, tnode, C, ...]; % 管段数据 % 热力系统 ies.heat.baseMW 10; ies.heat.node [node_i, phi_load, ...]; % 节点数据 ies.heat.pipe [fnode, tnode, K, L, d, ...]; % 管道数据 % 耦合设备 ies.coup(1).type chp; ies.coup(1).ebus 5; % 接入电网节点 ies.coup(1).gnode 2; % 接入气网节点 ies.coup(1).hnode 2; % 接入热网节点 ies.coup(1).P 3.0; % 额定电功率 MW ies.coup(1).Rh 1.2; % 热电比 ies.coup(1).eta_e 0.42; % 发电效率3.3 状态变量全局索引代码里最容易乱的地方统一求解法的代码里最核心也最容易写乱的地方就是状态变量的全局索引。我强烈建议在02_build_index.m里把所有变量的起始位置一次性算好存成一张索引表后续残差和雅可比组装都从索引表取位置不要在函数里反复计算。function idx build_index(ies) n_elec_bus size(ies.elec.bus, 1); n_gas_node size(ies.gas.node, 1); n_heat_node size(ies.heat.node, 1); n_pipe size(ies.heat.pipe, 1); % 电力状态变量除平衡节点外所有Va 所有PQ节点Vm idx.Va (1:n_elec_bus-1); idx.Vm (n_elec_bus : n_elec_bus n_pq - 1); % 天然气状态变量全部节点压力 idx.p_gas (max(idx.Vm)1 : max(idx.Vm)n_gas_node); % 热力状态变量管道流量 节点供水温度 节点回水温度经代数消元后 idx.m_pipe ...; idx.Ts ...; idx.Tr ...; % 耦合设备的状态变量如有需要额外引入 idx.coup ...; % 残差向量的索引与之一一对应 idx.F idx; % 实际工程中常用相同索引结构 end初值设定同样在03_init_guess.m里集中处理。电压相角统一设0电压幅值设1.0。气网节点压力统一设成气源出口压力附近比如0.8 MPa。热力部分是最需要谨慎的供水温度设设计值90℃回水温度设设计值50℃。千万不要把回水温度初值设成环境温度20℃那样供回水温差巨大热功率方程里残差极小但雅可比元素也很小数值上会产生严重的尺度问题。4. 核心求解器实现组装雅可比矩阵的关键细节4.1 残差函数怎么写才高效残差函数是整个求解器的核心我见过很多初学者把它写得又慢又难调试。最基础的几个原则首先是所有映射关系在build_index阶段就固化好残差函数里只查表其次要尽量矢量化能用矩阵运算就绝不用for循环。比如电力残差的计算从节点导纳矩阵直接切片就能一次算出全部节点的注入功率不需要逐节点循环。function F ies_residual(x, ies, idx) F zeros(size(x)); % 1. 取出各子系统的状态变量 Va x(idx.Va); Vm x(idx.Vm); p_gas x(idx.p_gas); m_pipe x(idx.m_pipe); Ts x(idx.Ts); Tr x(idx.Tr); % 2. 电力残差此处仅示意核心逻辑 V Vm .* exp(1j * Va); % 构建复电压向量需完整还原全节点 S_calc V .* conj(Ybus * V); % 节点注入复功率 F(idx.F_P) real(S_calc) - P_spec; F(idx.F_Q) imag(S_calc) - Q_spec; % 3. 天然气残差 % 根据Weymouth方程计算各管段流量再按节点累加 % 耦合设备的耗气量由当前状态变量计算后计入F_load_device % 4. 热力残差 % 水力方程、节点温度平衡方程分别组装 % 5. 耦合设备方程 % CHP、电锅炉等自身的等式约束 end写残差函数有个非常实用的调试技巧先单独调用一次能流计算程序把某个子系统节点上已知的解析结果代入残差函数检查该部分残差是否为零向量。如果该部分不为零说明这一块的公式或索引写错了。按子系统逐个验证能把调试时间缩短一半以上。4.2 雅可比矩阵的分块组装与验证雅可比矩阵是统一求解法里工作量最大、最容易被细节坑到的部分。整体结构是一个分块矩阵电力变量 气网变量 热力变量 电力方程 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、J_GE这些就是耦合贡献全部来自耦合设备。比如CHP接入后电力方程的残差对气网压力变量没有直接偏导因为CHP的耗气量影响的是气网方程但气网方程的残差对电力变量有偏导因为CHP耗气量是发电出力的函数。组装时用稀疏矩阵累加即可不需要等一个稠密矩阵。function J ies_jacobian(x, ies, idx) N length(x); J sparse(N, N); % 电力块 J(idx.F_P, idx.Va) dFp_dVa; J(idx.F_P, idx.Vm) dFp_dVm; % ... 其他子系统块同理 % 耦合块 % CHP气网方程对电出力的偏导 J(idx.F_gas, idx.Va) J(idx.F_gas, idx.Va) dFgas_dVa_chp; % P2G气网方程对电功率的偏导负号因为P2G消耗电产生气 % 电锅炉热网方程对电功率的偏导 end关于雅可比矩阵的验证我最推荐先用数值差分做一次对照J_num zeros(N, N); h 1e-6; for i 1:N xp x; xp(i) xp(i) h; xm x; xm(i) xm(i) - h; J_num(:, i) (F(xp) - F(xm)) / (2*h); end数值雅可比不用跑很多次就取一个典型运行点和解析雅可比做差看到最大偏差在1e-5量级即可认为解析推导正确。之后正式求解就用解析雅可比速度快且精度稳定。还有一个工程建议如果不想手推所有偏导第一版可以全程用数值雅可比跑通逻辑跑通后再逐步替换解析块。以前端工程思维来看这属于典型的“先让它跑起来再优化”策略。4.3 迭代收敛判据的坑收敛判据我建议同时盯两个指标残差范数max(|F|)和修正量范数max(|dx|)。只盯一个指标会出问题我实际遇到过残差已经降到1e-10量级、但状态变量还在小幅变飘的情况——这通常是因为方程之间存在弱线性相关多解或近多解导致的。所以代码里两个判据同时满足才认为收敛tol_F 1e-8; tol_dx 1e-10; for iter 1:max_iter F ies_residual(x, ies, idx); J ies_jacobian(x, ies, idx); if max(abs(F)) tol_F % 记录残差已达标 flag_F_ok true; else flag_F_ok false; end dx -J \ F; x x dx; if max(abs(dx)) tol_dx flag_F_ok disp([收敛于第 , num2str(iter), 次迭代]); break; end % 不收敛时输出当前状态方便追踪 if mod(iter, 5) 0 fprintf(iter%d, maxF%.2e, maxdx%.2e\n, iter, max(abs(F)), max(abs(dx))); end end另外不同物理量纲的变量混在一个向量里x的修正量范数会被大数值量级的变量主导。比如节点温度是几十的量级而电压幅值标幺值是1左右直接比max(|dx|)相当于默认温度修正更重要。严谨的写法是加“伪标幺化”——把每个分量的修正量除以该变量的基准值再取范数这样各变量能公平参与收敛判断。5. 算例设计与结果验证怎么证明代码算对了5.1 一个可复现的小型电气热耦合测试系统我自建了一套测试系统数据原型都来自公开文献改出来的方便复现。电网采用IEEE 33节点配电网经过改造后系统包括一个33节点配电网、一个6节点天然气网和一个6节点热力网。核心参数如下表子系统基准值主要参数配电网10 MVA12.66 kV33节点单回辐射网天然气网1.0 MPa6节点1个气源含2条长输管道热力网10 MW6节点3个热负荷设计供/回水温度90/50℃耦合设备接入位置和参数设备接入电网节点接入气网节点接入热网节点额定功率效率/热电比CHP5223 MW发电效率0.42热电比1.2电锅炉12-51.5 MW效率0.95燃气锅炉-442 MW效率0.90这套系统规模不大但麻雀虽小五脏俱全电气热三个网络的耦合要素CHP、电锅炉、燃气锅炉全部覆盖。如果想测试P2G可以在节点12加一组500 kW的P2G设备气网侧接入节点5。5.2 耦合与解耦的结果对比我做了两组计算一组是严格耦合求解另一组是解耦计算把CHP电出力、耗气量、热出力全部预设为固定值。结果对比如下指标解耦计算耦合计算偏差节点5电压幅值0.9820.9900.8%气网节点3压力0.812 MPa0.777 MPa4.3%热网节点3供水温度88.6℃89.1℃0.5℃系统总购电总购气折合一次能源21.3 MW20.6 MW3.3%这个结果说明了什么在只含一台CHP和一台电锅炉的“中等耦合”场景下解耦计算电压偏差1%左右也许还能忍受但气网压力偏差4.3%已经能影响管网的运行安全性判断了——如果解耦算出来节点压力在安全范围内实际耦合运行却可能低于最低供气压力。而P2G加入后偏差量级会进一步放大。能流计算如果连压力是否越限都判断错后续的优化调度、可靠性评估全都建立在错误的基础上。5.3 收敛性分析和异常表现同一套数据我用三种策略分别测试统一牛顿法冷启动6次迭代收敛到1e-8过程平稳。分解法17次迭代勉强收敛中间出现两次明显的震荡如果耦合设备再多一点大概率直接发散。两阶段法先分解5轮再切统一法总计7次迭代收敛比冷启动统一法更稳定。冷启动统一牛顿法之所以也能收敛是因为这套测试系统的初值给得比较合理尤其是热力网络回水温度设了50℃而不是20℃。初值一差收敛性和稳定性立刻变脸。调初值的过程我总结成一句话电网初值随便给气网初值靠气源热网初值靠设计——热力部分最敏感。求解过程中还有一个常见异常迭代发散的征兆通常是max(F)先降后升反复震荡。遇到这种情况先别急着改阻尼参数优先检查雅可比矩阵的正确性尤其检查耦合设备相关的非对角块符号。我调试时踩过一次最大的坑就是CHP耗气量对气网方程贡献的符号写反了导致整个系统变成一个“正反馈”回路怎么迭代都发散。6. 工程实践中的单位、初值与数值病态6.1 单位制统一是第一优先级电气热混合计算里单位制是最大的隐性坑。电力系统通常用标幺值气网用MPa、m³/h热网用MW、kg/s、℃。三个子系统的残差和雅可比元素混合在一起时数值量级可以差出10^6倍以上。直接组装成一个大矩阵条件数会非常难看求解精度大幅下降。我的做法是各子系统内部保留自己的习惯单位但在组装全局残差和雅可比之前先把所有残差除以本子系统的基准值做一个“伪标幺化”。比如气体残差除以该气网节点的基准流量kg/s热力残差除以基准热功率MW。这样全局残差各分量都在同一数量级后续迭代收敛判据也更合理。6.2 阻尼策略和数值病态的应对即使初值给得不错强耦合系统在牛顿迭代中途也可能迈出过大的步子。我的代码里加了一个简单的回溯线搜索当||F(x αdx)|| (1-ρ)||F(x)||时把步长α减半重试。ρ取0.10.3加这个阻尼后整体迭代次数会多几次但换来的是大幅提升的稳定性在工程里完全值得。雅可比矩阵用sparse存储、用左除求解千万别用矩阵求逆。另外如果条件数实在太大可以在求解前做行列缩放、改善条件数这个操作对大型系统效果明显。6.3 模型简化的边界要清楚最后聊聊模型简化的尺度。我见过不少新手一上来就堆复杂模型——压缩机动态方程、管道暂态过程、P2G多步化学反应全塞进去代码跑不动不说稳态能流的本质反而被淹没了。我的建议是第一版代码里做这样几个简化天然气压缩机用定升压比模型暂不引入变转速、变工况性能曲线。热力管道温降先用平均温度近似管段较短时可忽略温降只有长输管道才计入λL/(C_p m)项。P2G产氢按效率直接折算但必须带掺氢比例约束。供热管网忽略网络热惯性稳态工况下这不影响能流结果。这些简化不是偷懒而是稳态能流计算本身就要回答的核心问题——在给定的运行条件下电气热三个网络各自的状态量是多少。动态特性、变工况性能是后续动态仿真和优化研究的范畴混在一起只会让问题失去焦点。把稳态能流的计算逻辑先走通后面再逐步细化这条路比我一开始追求模型完备要高效得多。