1. 从一次热电联产联调事故说起为什么单算三种网络会翻车去年我做了一个园区级综合能源系统的规划评估项目系统里有中压配电网、中压燃气管网、热水供热管网还装了两台燃气轮机组、三台热泵、两台燃气锅炉。评审前我想省事用三个成熟的独立程序分别算电、气、热电网潮流程序假设燃气轮机满发气网程序按最大供气负荷定容热网程序把CHP热出力当成定值。结果显示电网无过载、气网容量够、热网温度达标看起来完美。现场边界的几条联调管线焊完DCS系统一通上问题全暴露了某一回线夏季高负荷时母线电压跌到0.90 pu以下气网末端压力比运行规程低了15%热网最远端用户供水温度只有62℃设计值可是75℃。三个独立程序算出来的最优方案在真实耦合运行下根本兜不住。问题出在哪就在于这三个网络根本不应该分开算。电网的燃气轮机烧的是天然气气网的压力波动会直接影响发电出力热泵要从电网取电才能向热网供热电网电压跌了热泵制热效率跟着掉燃气锅炉向热网供热的同时也在气网里抽气。任何一个网络只要有波动都会通过耦合设备传导出去。这就是标题里计及多能耦合要解决的核心问题——电气热能流计算本质上是在联立求解三个网络方程和一个耦合设备方程组的数学问题。所以这篇文章我把这套计算的物理模型、数学方程、Matlab实现思路以及调试经验从头到尾串一遍。适合正在做综合能源方向的研究生、搞园区能源规划的工程师以及想从独立潮流过渡到多能流计算但找不到完整参考的人。2. 电气热三种网络的数学模型先弄清楚每个网络在解什么未知量2.1 电力网络牛拉法潮流是基础但耦合节点的处理方式不同电力网络这一步没什么新鲜玩意儿就是常规牛顿-拉夫逊潮流计算。节点分三类PQ节点给定有功和无功、PV节点给定有功和电压幅值、平衡节点给定电压幅值和相角吸收全网功率差额。状态变量是各节点电压幅值V和相角θ极坐标下的功率不平衡方程ΔP_i P_i_spec - V_i * Σ V_j * (G_ij * cosθ_ij B_ij * sinθ_ij) 0 ΔQ_i Q_i_spec - V_i * Σ V_j * (G_ij * sinθ_ij - B_ij * cosθ_ij) 0但在综合能源系统里要注意耦合节点不再是简单的PQ/PV。燃气轮机节点在电网模型里通常当成PV节点因为原动机有功出力可是由气网侧的天然气供应量决定的不是独立给定的。这就在不同网络之间引入了相互制约的耦合变量需要把气网的气压和CHP的气耗量纳入整个联立方程组求解框架。热泵节点更特殊——它既是电网的功率负荷又是热网的固定热出力源电功率需求本身由热网侧的制热需求决定。2.2 天然气网络Weymouth方程和节点气压很多人第一次栽在这天然气网络和电网在数学上很像都是节点流量平衡加支路特性方程。高压输气管道通常用Weymouth方程描述稳态气流F_ij^2 C_ij^2 * (p_i^2 - p_j^2)其中F_ij是管道流量kg/s或m³/hp是节点绝对压力C_ij是和管道直径、长度、摩擦系数相关的常数。变量上流量F相当于电流压力p相当于电压但注意是压力的平方差。每个节点满足流量平衡Σ F_in - Σ F_out F_source - F_load 0节点类型分两类压力给定节点类似电网的平衡节点比如气源/门站节点和流量给定节点类似PQ节点气负荷节点。状态变量是全网节点压力如果初始压力给得离谱方程的平方项会让雅可比矩阵直接病态。压缩机的建模是热电联算里最容易漏的一块。电驱动压缩机要从电网取电燃气驱动压缩机要从气网取气压缩机耗气/耗电功率普遍按以下经验公式折算P_comp k * F * (p_out/p_in)^m - 1在综合能源系统里如果压缩机是电驱动的那它就成了电网和气网之间的又一个耦合环节不建这个模型气网压力可能算得虚高。2.3 热力网络的水力-热力双层模型热网比电网和气网都麻烦因为它有水力过程水泵、流量、压力和热力过程供回水温度、热负荷、管道散热两个子模型而且这俩是嵌套耦合的——水力工况决定各管段流量分布流量分布又决定节点混合温度和管道散热反过来温度又影响热源侧的出力需求。水力模型就是管网的基尔霍夫定律。对每个节点流量平衡Σ Q_in - Σ Q_out - Q_load 0每个环路压降平衡Σ Δp_loop 0管段压降和流量的关系通常用Darcy-Weisbach或简化幂函数Δp R * Q^1.75 到 Q^2。热力模型是另一套方程。节点供水温度T_s和回水温度T_r是两个核心状态变量管道末端温度按下式计算T_end T_ground (T_start - T_ground) * exp(-λ * L / (c_p * Q))节点混合温度则根据各支路流量加权平均。热负荷节点的用热功率由下式确定P_heat Q * c_p * (T_s - T_r)注意这里的Q是流过负荷节点的质量流量T_s和T_r是供回水温度。热网稳态计算中通常给定热负荷和部分节点的供回水温度求解全网流量和温度分布。实际算的时候水力迭代和热力迭代要交替进行我先固定温度算流量再用新流量更新温度反复迭代直到两者都收敛。2.4 三种网络的时间常数差异与稳态假设电网的动态过程是毫秒级气网是秒到分钟级热网的热惯性是小时级甚至更长。要做稳态多能流计算就得接受一个前提假定系统运行在某一稳定工况忽略动态过程。这在工程规划阶段完全够用——我们算的是在某个典型运行方式下系统是否能满足供需平衡不是暂态稳定性。但理解这个时间尺度差异很重要它决定了迭代策略热网通常最后收敛因为它的响应最慢初值误差可能要到最外层循环才体现出来。3. 耦合设备建模算例能不能收敛一大半看这几个环节3.1 CHP机组电热双输出的能量路由器CHP是多能流计算里最核心的耦合设备。燃气轮机的燃料气通过燃烧做功发电余热通过换热器供热。理想化模型可以写成F_gas * LHV P_elec / η_e P_heat / η_h其中F_gas是燃料消耗量LHV是天然气低位热值η_e是发电效率η_h是热回收效率。实际建模更常用的是电力-热力可行域约束背压式机组热电比恒定PQ关系是一个固定斜率抽凝式机组的热出力可以在一定范围内调节电出力也会随之变化。多能流计算中通常把CHP当成气网的一个节点负荷消耗天然气电网的一个有功源发出电功率热网的一个热源输出热功率三者通过上面的方程绑定在一起。建这个模型时要格外小心变量的正方向和量纲。气网里CHP耗气量是负数源电网里是正功率注入热网里是热源注入三个网络各自的方程对同一个CHP的变量符号约定完全不一样。我在实现时踩的坑是电网变量用标幺值气网用m³/h热网用kW导致耦合矩阵里雅可比元素数量级差了好几个数量级一度收敛很慢。后来统一改写各网络计算函数让它们都返回国际单位制下的SI核心量纲变量再在各自函数内部转换成标幺值计算联立求解才稳定下来。3.2 热泵和燃气锅炉两个好算但容易被忽略的环节热泵模型简单Q_heat COP * P_elecCOP制热性能系数通常取固定值或随工况变化的经验曲线。热泵在电网是负荷、在热网是热源它的电功率需求在联立方程组里不是独立给定的而是由热网侧需要它输出的热量除以COP反算出来。这个反算逻辑要写清楚不然容易把功率方向弄反。燃气锅炉更简单P_heat η_b * F_gas * LHV本质上是天然气流量和热出力之间的线性/准线性关系。但要注意锅炉的烟气往往需要引风机、循环泵这些设备本身还要耗电严格算的话还得把这部分电耗算进电网负荷。如果是做高精度校核计算这些是不能忽略的但一般规划阶段可以归到站用电的固定比例里。3.3 能量枢纽耦合矩阵把不同设备抽象成统一形式多能流计算发展到后面大家喜欢用能量枢纽Energy Hub的方式统一描述耦合设备。核心思想就是把一个站点的输入电、天然气、太阳能等通过一个耦合矩阵C映射到输出电、热、冷等[ P_e_out ] [ C_ee C_eh ] [ P_e_in ] [ P_h_out ] [ C_he C_hh ] * [ P_g_in ]矩阵中的元素表示耦合设备的效率或分配系数。CHP在能量枢纽里表现为C_he天然气转热和C_ee天然气转电同时非零热泵表现为C_eh电转热。用能量枢纽建模的好处是不管站点里装了多少台设备都可以先综合成汇总矩阵再接入整体计算程序代码模块复用性很好。我在Matlab里把每个站点的耦合矩阵定义成struct存起来再统一组装到全网方程组里调试时只看矩阵元素就能发现问题出在哪个设备上。4. Matlab实现路线统一求解和分解迭代我都试过各有利弊4.1 统一求解法大雅可比矩阵怎么拼出来统一求解法是把电网、气网、热网的所有方程堆成一个大的非线性方程组F(X)0牛顿-拉夫逊迭代一次同时修正所有状态变量。整个系统的状态变量向量长这样X [V_e; θ_e; p_g; T_s; T_r; Q_flow]其中V_e和θ_e是电网每个非平衡节点的电压幅值和相角p_g是气网非压力给定节点的压力T_s和T_r是热网节点的供回水温度Q_flow是热网各管段流量或者不显式列流量用水力方程直接求解。雅可比矩阵是分块结构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、J_ge等就是耦合设备带来的交叉导数。CHP的存在意味着气网的气压会影响电网注入功率导致J_ge非零热泵的存在意味着热网负荷会影响电网功率导致J_eh非零。统一求解最大的优点是强耦合工况下收敛性好电网气网热网的变化在一个迭代步里同时传递不会出现分解迭代那种一个网络已经算完了另一个还没反应的滞后问题。缺点也很明显编程量大、雅可比矩阵维度大一个中等规模系统动辄几百上千维、初值要求高、程序耦合度高没法复用现有代码。调试时我通常把JOFF标志位打开单独用解析雅可比和数值差分雅可比做对比交叉导数错误的概率很高。4.2 分解迭代法模块化清晰但收敛条件要小心分解迭代法是网与网之间用松弛法迭代。基本流程给定CHP、热泵等耦合设备出力初值。把CHP电出力作为电网注入功率把热泵当成电网负荷单独求解电网潮流。把CHP耗气量、锅炉耗气量作为气网负荷求解气网。把CHP热出力、锅炉热出力、热泵热出力作为热源求解热网。根据气网算出的CHP可用气量重新修正CHP的电出力根据热网算出的热负荷需求重新修正热泵和锅炉的热出力。返回第二步直到耦合变量前后两次迭代的差值小于允许误差。分解迭代的优点是各网络可以分别用成熟程序求解程序组织清晰。缺点是收敛性完全取决于耦合强度——如果耦合设备功率占比大、或者热电比很高简单的高斯迭代就会发散。我在一个热电比0.95的算例里第一次跑分解迭代出现了典型的往返震荡CHP电出力一会儿0.8 MW一会儿1.3 MW气网压力也跟着来回跳最后只能降低松弛因子到0.4才勉强压住。两种方法的关键对比如下对比项统一求解法分解迭代法收敛速度快二次收敛慢线性收敛/可能震荡初值要求高各网络都要给合理初值相对宽松各网络独立初值编程复杂高大雅可比矩阵易出错中模块可复用强耦合适应性好差需松弛因子配合程序可维护性低高我个人建议如果做单次规划校验、不需要频繁改系统拓扑用统一求解法如果要做随机多场景分析、需要在不同案例间频繁切换网络参数用分解迭代法更现实配上阻尼机制能稳住大半场景。4.3 我在Matlab里搭建的实现框架和关键代码我基于统一求解法实现了一套代码核心框架分三层网络参数层、耦合设备层、求解核心层。网络参数层用struct存电网的bus/branch数据、气网的node/pipe数据、热网的node/pipe数据耦合设备层定义CHP、热泵、锅炉、压缩机各自的参数、连接关系和模型函数求解核心层负责组装雅可比矩阵、迭代求解、结果整理。主程序流程可以概括为%% 主流程 % 1. 定义网络和设备参数示例 grid define_electric_network(); % 电网bus/branch gas define_gas_network(); % 气网node/pipe heat define_heat_network(); % 热网node/pipe units define_coupling_units(); % CHP/热泵/锅炉/压缩机 % 2. 初始化状态变量 X0 initialization(grid, gas, heat, units); % 3. 牛顿-拉夫逊迭代 for iter 1:max_iter [F, J] build_residual_and_jacobian(X0, grid, gas, heat, units); delta_X -J \ F; % 阻尼处理 alpha damping_factor(F, delta_X, residual_norm); X0 X0 alpha * delta_X; if norm(F, inf) tol break; end end组装残差和雅可比矩阵是整个程序最核心也最易错的部分。我建议把残差矢量组织成清晰的分段形式电网残差放前面、气网残差居中、热网残差放后面这样调试时看哪个分量不收敛一目了然。雅可比矩阵用稀疏矩阵表示避免大矩阵运算内存爆炸。一个重要技巧让耦合设备的残差方程尽量写成显式关系。比如CHP与其在耦合矩阵里写复杂的隐式方程不如直接在残差函数里写成wasted_heat 0; % 说明此处仅作为示意 R_chp_e unit.P_elec - unit.eta_e * unit.F_gas * unit.LHV; R_chp_h unit.P_heat - unit.eta_h * unit.F_gas * unit.LHV;这样残差语义清晰雅可比矩阵的元素也能手算校核。另一个我后来才掌握的技巧先跑纯电网、纯气网、纯热网三个子程序全部收敛了再联立。如果联立后不收敛基本可以确定是耦合变量或交叉导数写错了。可以加一个开关把耦合项的雅可比设置为零矩阵验证残差是否对应独立的三个网络。5. 收敛性调试实录初值、阻尼与气网的非凸难题5.1 气网Weymouth方程的非凸性迭代为什么会在两个解之间跳Weymouth方程F² C²(p_i² - p_j²)里的流量和压差是平方关系它的解域是凸的吗不是。同一个F可以对应两个符号相反的压差方程在F0处不光滑。牛顿法在高非线性区域很容易在多个根之间跳跃甚至震荡。我在调一个含压缩机的气网算例时出现了压力在1.2 MPa和0.8 MPa之间来回跳的现象雅可比残差总是降不到1e-6。解决方法是给流量变量一个初值符号约束根据管道连接关系预判方向同时将方程改写为R_ij F_ij - C_ij * sqrt(p_i^2 - p_j^2) * sign(p_i - p_j)这样能避免负值平方根带来的NaN问题也能大幅改善迭代稳定性。5.2 热网冷启动困难温度初值怎么给才合理热网计算里最容易犯的错误是一上来就零初值。给定热负荷为1.5 MW的节点如果初始温度设成0℃绝对温度管道的散热模型和节点混合方程会产生巨大的负雅可比元素导致迭代直接发散。我的经验是用略高于环境的温度值作为初值比如统一设供温75℃、回温45℃、环境温度15℃然后跑两层固定点迭代把它稳定下来再进入牛顿主循环。另外热网的流量和温度量级差异大流量几百kg/h温度几十度不归一化的话雅可比矩阵的条件数会很差。我给每个变量指定了一个基准值比如流量基准1000 kg/h、温度基准100℃在雅可比组装后统一缩放条件数立刻改善了两个数量级。5.3 我常用的四个工程化收敛辅助手段如果加阻尼之后还是收敛不理想我一般按以下顺序排查第一个手段是残差归一化。电网的有功不平衡量级在MW、无功在MVar气网的流量不平衡在kg/s热网的功率不平衡在kW直接比较它们的范数会掩盖小量级分量的问题。我通常把各类残差除以其所在网络的总负荷量级统一变成无量纲相对误差再统一看最大值。第二个手段是迭代步长限制。每次修正量δX不得超过该状态变量初值的某个百分比比如电压不超过0.1 pu、压力不超过0.1 MPa、温度不超过5℃。防止牛顿法大步长跳过可行域。这个可以用以下代码简单实现max_step 0.1 * abs(X0 eps); delta_X sign(delta_X) .* min(abs(delta_X), max_step);第三个手段是迭代顺序调整。分解迭代法的顺序直接影响收敛速度。我试过几种顺序后认为先算气网、再算电网、最后算热网的顺序最稳原因是气网为CHP和锅炉提供了燃料边界确定了这个边界后电网热网的源项就相对确定了再算后两者不像先算电网那么盲目。第四个手段是中等规模系统下改用信赖域算法替代线搜索。Matlab的fsolve有信赖域选项多能流这类高度非线性的问题上比自带的高斯牛顿法稳得多。我在统一求解的主循环里没有自己写信赖域而是直接调用优化工具箱的fsolve配合自带的雅可比函数句柄省了很多麻烦。5.4 三类常见的假收敛以及排查思路有些情况下迭代停下来了但结果明显不对。我整理的排查思路是第一类是孤立网络平衡但耦合节点失衡。可能因为耦合设备残差没有加入整体残差向量各网络独自平衡了但CHP电出力和气耗量之间却不满足设备模型。排查办法是算完后把每个耦合设备的输入输出都打印出来人工代入换算效率看和设定值是否吻合。第二类是压力/温度出现物理量纲不合理但迭代仍然继续。可能原因是初值方向错误迭代走到了另一个数学解。比如热网远端回水温度算出来比供水还高这种解满足方程组但从物理上看不可能。排查办法是对状态变量加边界约束让残差函数在越界时返回一个大数惩罚或者直接在一次牛顿步后判断上下界。第三类是迭代收敛但残差极小但流量不平衡。因为在管网模型中某些循环管路环网里的流量解不唯一水力方程在某些构型下存在多解。解决办法是检查环路独立方程是否被从雅可比矩阵中错误地删除了或者对环网添加基准流量约束。6. 一个典型算例的结果解读从迭代残差到工程判断6.1 算例配置与初始条件我用一个中等规模的测试系统演示总体效果33节点配电网、20节点气网含2个气源节点、1个压缩机站、16节点热网含1个热源站。耦合设备配置为2台CHP每台额定电出力1.5 MW热电比0.9、1台热泵COP 3.5额定热出力1.2 MW、2台燃气锅炉每台热出力2 MW。总电负荷约8.5 MW总气负荷约4.2 MW等效总热负荷约6.8 MW。迭代初值电网各节点电压1.0 pu、相角0气网各节点压力统一设2.0 MPa热网供水温度75℃、回水温度45℃各管段流量按负荷近似预分配。迭代目标最大残差小于1e-7无量纲归一后。6.2 迭代过程与收敛特征这个算例用统一求解法第1次迭代残差降到3.7e-1第3次降到2.1e-3第5次落到8.6e-6第6次达标。对比一下同样算例用分解迭代法需要24次迭代且由于耦合设备功率占比偏高松弛因子设到0.45才不震荡。收敛后各网络的关键指标指标数值电网最大电压偏差5.7%节点12距CHP接入节点较远气网最低节点压力1.33 MPa末端气负荷节点热网供回水温差最大偏差7.8℃平衡节点净注入有功4.9 MWCHP1电出力/热出力1.32 MW / 1.19 MW热泵耗电功率/热出力0.38 MW / 1.33 MW气源总供气量3890 m³/h标况6.3 怎么判断这个结果算对了三类平衡校验多能流计算有一个优势是天然存在多组可交叉验证的平衡关系只要有一个对不上基本可以确定算错了电力平衡平衡节点注入功率 CHP发电 - 常规用电负荷 - 热泵电耗 - 电驱动压缩机功耗 0网损约束。天然气平衡气源供气量 CHP耗气量 锅炉耗气量 压缩机能耗 气负荷这个用燃气负荷总量和Weymouth方程算出的管段流量差值验证。热力平衡CHP热出力 锅炉热出力 热泵热出力 热网负荷总和 管道散热损失。这个特别容易忽略管道散热我在一个算例里发现热源总供热比热负荷多出12%一开始以为是计算错误后来单独算了散热才发现是管道损失占的。除了三个网络的平衡还要检查耦合设备自身的模型一致性。CHP的发电效率η_e P_elec / (F_gas * LHV)要落在合理区间0.25~0.45热泵的实际COP Q_heat / P_elec要接近设定值偏差不超过5%如果偏差过大说明耦合方程在联立解中被带偏了优先检查雅可比矩阵中对应块的交叉导数是否正确。6.4 联立计算后发现了哪些独立计算看不到的现象这个算例最有意思的地方在于联立计算后才暴露出来的问题当热网负荷从6.8 MW增加到8.0 MW时热泵热出力需求增加电耗从0.38 MW增加到0.43 MW电网潮流重新分布CHP热出力需求增加气耗量上升气网末端压力从1.33 MPa跌到1.20 MPa气网压力下降后燃气轮机的出力能力受限平衡节点被迫多注入0.6 MW。这串连锁反应如果还用三个网络独立计算完全看不到——独立计算会假设热泵电耗固定、CHP气耗固定结果必然低估了高负荷时电网和气网的压力。这也是为什么我一直认为多能流计算不是学术上的炫技它解决的是工程实际中的账目对不上问题。做园区能源规划、做设备选型、做运行方式校核之前先做一次全耦合的电气热能流计算比任何经验估算都靠谱。7. 给入门者的三个实用建议最后分享几点实际做项目过程中的经验。第一不要一上来就追求大规模高精度。先把一个最简单的三节点电网、两节点气网、两个节点热网跑通用线性方程组手解的结果验证程序正确再逐步加节点、加设备。我一开始就奔着33节点去结果调了一个月后来用5节点小系统三天就通了。第二Matlab里做这个研究工具箱够用就好。本质上核心代码就是牛顿-拉夫逊迭代加上稀疏矩阵运算不需要一堆专业工具箱。我在项目里只用了基础的矩阵运算\操作符、稀疏矩阵函数sparse、nnz其他都是自己写的函数。优化工具箱的fsolve可以作为备选解法器。第三变量和量纲管理是整个程序的命门。建议用一套统一的单位体系贯穿始终比如电网用标幺值、气网用MPa和kg/s、热网用kW和kg/h在接口函数里统一换算。每写一个设备模型就写一个单元测试函数验证输入输出是否正确别等联立了再排查。多能流计算这块内容能聊的细节远远不止这些但把上述这套模型和实现方法吃透常见的电气热联算需求基本都能覆盖。后续如果想深入可以做动态多能流、考虑不确定性场景、结合优化算法做日前调度底层框架完全可以直接复用。