热力管道热惯性建模与虚拟储能优化调度实践

前阵子做一个园区综合能源系统调度的仿真课题一开始图省事把供热网络当成一个简单的“热力传送带”热量从热源出发当天发当天到完全忽略管道里那几百吨水的存在。结果系统级仿真一跑问题全露出来了某些时刻热源明明已经把出力压了下来用户侧回水温度还在往上飘某些时刻热源加出力用户侧升温却姗姗来迟。跟实际运行数据一对偏差大得没法看。后来老老实实把热力管道的热惯性纳进模型用有限差分法求解管道的瞬态温度场再把热网蓄热能力折算成虚拟储能加进调度模型结果才终于对得上。这篇就把整个实现过程复盘一遍从控制方程、离散格式、虚拟储能量化到Matlab代码架构和调试踩坑完整走一条可复现的技术路线。这个工作适合两类人看一类是做综合能源系统优化调度研究的需要把热网动态特性纳入调度模型另一类是搞热网仿真或供热运行调度的工程师想给传统“以热定电”的思路加一点柔性空间。文章里涉及的物理模型、公式和代码都以可落地为首要目标不会停留在概念分析层面。1. 热惯性不是“可选项”它是综合能源系统调度里绕不开的物理约束1.1 供热网络到底“重”在哪先算一笔简单的账。一根DN500的供热管道内径大概0.5米左右跑1公里长里面装的水量是V π × (0.25)² × 1000 ≈ 196 m³也就是约196吨水。水的比热容是4.186 kJ/(kg·K)这根管道每升高1℃就要吸收大约820 MJ的热量。如果是热水管网温差动辄三四十度一条十公里长的管线整体热容量轻松上到几百吉焦。这种数量级的能量足够一个中型热源站满负荷运行好几个小时。这意味着什么意味着供热网络本质上是一个巨大的、随着温度变化不断吸放热的蓄能体。它的时间常数不是秒级、分钟级而是小时级。在调度时间尺度通常1小时或15分钟内管网的充放热过程会直接改变热源出力与用户负荷之间的时序匹配关系。1.2 忽略热惯性导致的典型调度错误如果把热网当静态模型默认“热源供了多少热用户侧立刻就用掉多少”调度模型会犯两类典型错误。第一类是过拟合尖峰负荷。比如早晨供热负荷陡增静态模型会认为热源必须立刻跟着加出力于是调度方案里就会出现一个很高的出力尖峰。但实际上由于管网里存着大量热水热源可以提前半小时开始缓慢提温利用管网蓄热来“熨平”这个尖峰做到更小的热源装机或更平稳的出力曲线。第二类是热量“错位”。调度模型在某个时段刻意压低热源出力以为用户侧温度会立刻下降但实际上管网和建筑的热惯性会把这个影响延迟一两个小时。这就会导致实际运行中出现“热源已经减了、用户侧还在热”“热源已经加了、用户侧还没暖起来”的现象。很多做综合能源调度的人把精力全放在设备模型和优化算法上结果误差的根源其实在热网这一环。这也是为什么说热惯性不是锦上添花的精度提升而是影响调度结果正确性的基础约束。2. 热力管道动态建模先把控制方程、边界条件和离散格式定明白2.1 一维对流-扩散方程与关键参数供热管道内热水的温度变化在轴向方向占主导径向温度梯度通过集总参数的方式折算到散热损失里。工程上用得最多的是下面这个一维对流-扩散方程ρ c_p A (∂T/∂t) -ρ c_p A u (∂T/∂x) λ A (∂²T/∂x²) - K_loss (T - T_env)其中ρ为水的密度取1000 kg/m³c_p为水的比热容取4186 J/(kg·K)A为管道内截面积单位m²u为管道内流速单位m/s由流量除以截面积得到λ为水的导热系数约0.6 W/(m·K)K_loss为管道单位长度综合散热系数单位W/(m·K)这里面对流项第二项是主导项它决定了热水从热源到用户侧的“传输延迟”扩散项通常很小但在某些低流速工况下不能完全无视散热损失项则决定了长距离输送时的沿途温降对回水温度预测影响很大。2.2 为什么选显式迎风差分格式对这样的方程做数值求解有限差分法是最直接的路径。差分格式有很多种我在这个项目里选了显式一阶迎风差分原因是它在满足稳定条件的前提下编程最简单、物理意义最直观、调试最省心。迎风差分的核心思想是热水是沿x轴正方向流动的那么k节点下一时刻的温度应该受“上游”k-1节点而不是下游k1节点的影响这符合热水的物理输运特性。离散后每个内部节点的递推式是T_k^(n1) T_k^(n) - (u Δt / Δx) (T_k^(n) - T_(k-1)^(n)) (λ Δt / (ρ c_p Δx²)) (T_(k1)^(n) - 2T_k^(n) T_(k-1)^(n)) - (K_loss Δt / (ρ c_p A)) (T_k^(n) - T_env)三个修正项分别对应对流输运、导热扩散、沿途散热损失。实际写Matlab循环时用向量化写法可以避免显式for循环带来的效率问题但为了体现清晰的物理过程先按最基本的for循环写% 初始化 Nx 101; % 空间节点数 dx L / (Nx - 1); % 空间步长 Nt round(T_total / dt); % 时间步数 T T_init * ones(1, Nx); % 初始温度场 T_new T; T_history zeros(Nt, Nx); T_history(1, :) T; for n 1 : Nt - 1 for k 2 : Nx - 1 conv u * (T(k) - T(k - 1)) / dx; diffu thermal_diff * (T(k 1) - 2 * T(k) T(k - 1)) / dx^2; loss K_loss / (A * rho * cp) * (T(k) - T_env); T_new(k) T(k) - dt * conv dt * diffu - dt * loss; end % 边界条件 T_new(1) T_supply; % 入口给定供水温度 T_new(Nx) T_new(Nx - 1); % 出口自由出流近似 T T_new; T_history(n 1, :) T; end2.3 边界条件和初值设置里的门道管道入口是最容易处理的边界因为热源侧的供水温度是可控量直接给第一类边界条件T(0,t) T_supply(t)即可。如果你需要模拟“给定入口流量和入口温度”的情况也完全可以按上式把T_supply改成随时间变化的曲线。管道出口处的边界处理要小心。如果直接假设出口温度不变会人为引入一个虚假的热量堆积如果自由出流∂T/∂x 0则用户侧的回水温度会略微偏高但整体误差在工程可接受范围内。更严格的处理是把出口接到回水管道继续建模形成完整的供回水管网但这对模型复杂度提升很大且对调度问题的边际收益有限我的做法是保留供回水双管但出口采用自由出流条件。初值问题更隐蔽。第一次跑仿真如果直接给一个常数初值比如所有管道初始温度都是90℃那么前几个时间步会出现明显的“虚假瞬态”——管道的蓄热效应会把真实的热波给抹掉一部分。解决方法是先用恒定入口温度把管网跑到一个准稳态再把准稳态温度场作为初始值或者直接读入前一天该时段的历史温度场。2.4 稳定条件不是“建议”是硬约束显式格式最大的短板就是稳定性限制。有两个条件必须同时满足对流通项满足CFL条件C u Δt / Δx ≤ 1导热项满足傅里叶数限制Fo λ Δt / (ρ c_p Δx²) ≤ 0.5冷却水的热扩散系数约1.4×10⁻⁷ m²/s这个数值非常小所以傅里叶数基本不会成为限制条件。真正的约束来自CFL条件。如果流速是1 m/s管道用100个节点来剖分Δx 10 m那么最大允许时间步长就是10秒。这意味着模拟24小时需要8640步计算量完全可接受但如果你把空间步长取得太小时间步长就会被拖累整段仿真时间会指数级上升。我在实践中得到的经验是空间步长不需要取得太细10米到50米都足够保证工程精度关键是时间步长要同时满足CFL条件和调度模型的时间分辨率要求。如果调度步长是1小时但CFL条件要求最大计算步长为30秒那就在每个调度时段内部嵌套一个子循环把动态求解和调度决策解耦开。3. 虚拟储能量化把热管网折合成一台“温度储能电池”3.1 热网蓄热能力如何映射为储能参数热惯性在物理上表现为管网的蓄热能力在调度模型里就可以把它等价成一个虚拟储能装置。虚拟储能不需要额外建储热罐它利用的是供热管网本身已有的水容量和温度允许波动范围。这个概念跟电池储能高度相似温度升高相当于充电管网吸热温度降低相当于放电管网放热允许的温度波动范围决定了储能容量。具体量化上把供热管网等效为一个“虚拟储能电池”用三个参数描述:虚拟储能量E_vess当前管网热状态相对基准热状态的偏差虚拟储能功率P_vess单位时间内管网热量的增减速率荷电状态SOC_vess当前可用储能量与最大储能量的比值管网中第i段管道温度从T_base升高到T(t)时吸收的热量为E_vess(t) Σ M_i c_p (T_i(t) - T_base)其中M_i ρ A_i L_i是该段管道内的水质量。所有管段累加起来就是整个供热网络的虚拟储能量。3.2 一段DN500管道到底能“存”多少热用前面提到的那根1公里DN500管道来算一笔账。管内水质量M ≈ 196吨。如果用户侧允许的供水温度波动范围是±5℃也就是10℃的可利用温差那么单根管道的可利用虚拟储能容量是E_vess_max M c_p ΔT 196000 × 4186 × 10 ≈ 8.2 × 10⁹ J 8.2 GJ换算成更方便理解的电量单位1 GJ ≈ 277.78 kWh8.2 GJ ≈ 2278 kWh。这还只是一根1公里的供回水管路中供水管单向的容量。实际上供回水双管都有蓄热能力再算上整个供热半径内的全部管线一个中等规模园区的热网虚拟储能量很容易达到几十吉焦。相比之下同样要储存这么多能量需要定制一个几百立方米的常压储热水罐无论投资还是占地都不是一个量级。这就是热网虚拟储能在经济性上的天然优势。为了更直观我做了一个不同管径、不同温差下的虚拟储能容量对比表管道规格长度(km)可用温差(℃)可利用热容量(GJ)等效电量(MWh)DN300151.60.44DN500154.11.14DN50021016.44.56DN80031057.816.06这张表的意义在于做调度方案时一眼就能看出热网虚拟储能的调节潜力大概在什么量级从而判断它对削峰填谷能起到多大作用。3.3 从温度和流量变化推算可调功率虚拟储能的充放电功率本质是管网整体温度水平的变化率。如果把一段管道抽象成集总参数模型它的蓄放热功率可以写成P_vess(t) M c_p dT(t)/dt离散化之后在一个调度时段Δt内虚拟储能平均充放功率近似为P_vess(t) ≈ M c_p (T(t1) - T(t)) / Δt这个式子的物理含义很清晰调度决策让管网平均温度提升K度就等于让虚拟储能以某个功率持续充电K度对应的时间。如果Δt取1小时M c_p取管道水容量×比热那么每提1℃对应一个固定的充热功率。这个功率不能无限大。受限于管道的热交换速率、热源侧的升温能力和用户侧的温度容忍度虚拟储能的充放功率需要加上限约束P_min ≤ P_vess(t) ≤ P_max同时因为供回水温度都有运行上限和下限虚拟储能SOC也需要限制在安全区间E_vess_min ≤ E_vess(t) ≤ E_vess_max到这里虚拟储能就完成了从“热管网热状态”到“调度模型储能约束”的映射。4. 调度模型与Matlab代码架构热网动态约束是怎么进优化器的4.1 调度目标与约束设计虚拟储能进入调度模型之后优化问题的结构就清晰了。以典型的热电联产燃气锅炉电锅炉综合能源系统为例调度目标是在满足热负荷需求的前提下最小化整个调度周期内的运行成本min Σ_t ( c_gas · F_chp(t) c_gas · F_boiler(t) c_elec · P_eb(t) )其中F_chp(t)为燃气轮机在t时段的燃料耗量F_boiler(t)为燃气锅炉燃料耗量P_eb(t)为电锅炉耗电量。c_gas和c_elec分别为气价和电价。这里可以根据实际问题再加入碳排放成本、启停成本等。约束条件需要覆盖四个方面热源侧约束常规机组出力上下限、爬坡约束、启停逻辑热网热平衡约束任一时刻热源总供热量 用户热负荷 管网散热损失 虚拟储能充放功率虚拟储能状态转移方程E_vess(t1) E_vess(t) P_ch(t)ΔT - P_dis(t)ΔT温度边界约束各节点供回水温度必须在允许范围内注意到一点传统调度模型里的热平衡是瞬时平衡即“发多少热用多少热”而引入虚拟储能后热平衡变成了带存储项的平衡本质上是微分方程离散化后的形式。这一步是整个建模思路转变的核心。4.2 代码模块怎么拆Matlab实现时我把整个项目拆成五个模块各自独立。好处是每个模块可以单独测试改参数和换算例时不需要牵一发动全身。main.m 主脚本参数总控、调用各模块、输出结果 get_network_data.m 管网拓扑与参数读取 heat_load_profile.m 热负荷曲线生成或用真实历史数据 pipe_dynamics.m 管道有限差分求解器输出各节点温度场 virtual_storage.m 虚拟储能量化计算E_vess、P_vess、SOC optimize_schedule.m 调用linprog/intlinprog求解调度模型 plot_results.m 结果可视化温度曲线、出力曲线、SOC曲线pipe_dynamics.m是核心求解器。它接收热源供水温度曲线、流量曲线和室外温度用第二节的有限差分算法计算管道出口温度动态响应返回各节点在各时刻的温度矩阵和虚拟储能状态。optimize_schedule.m是调度求解器。它采用“动态模拟-优化迭代”的解耦结构先用上一轮的调度结果跑一遍热网动态模拟得到虚拟储能的可行域再把可行域带入优化模型求解更新调度结果如此迭代2-3轮即可收敛到稳定的调度方案。4.3 求解器选型linprog还是intlinprog还是YALMIPCplex如果整个优化模型是线性的热网管道动态方程在给定流量下是线性的只有温度一个状态变量直接调用Matlab优化工具箱自带的linprog就够了。我一开始的版本就是这样代码干净部署没有额外依赖适合快速验证算法。如果加入了机组启停的0-1变量、分时电价下的设备启停策略问题就变成混合整数线性规划MILP需要换成intlinprog。Matlab自带的intlinprog对小规模问题几十个0-1变量求解没有问题求解速度可以接受。但如果你把整个热网的每个节点温度都作为优化变量而不只是把虚拟储能聚合量作为变量那么优化问题的维度会急剧膨胀自带的intlinprog可能就要跑很久了。这种情况下建议用YALMIP作为建模层后端接Cplex或者Gurobi。YALMIP的好处是建模灵活约束表达式写起来直观切换求解器只需改一行代码ops sdpsettings(solver, cplex, verbose, 0); optimize(constraints, objective, ops);我这里选择的是把动态热网模型在优化外部做迭代降维优化器内部只保留虚拟储能状态变量这样一个24小时、15分钟分辨率的调度问题决策变量控制在几百个以内直接用linprog就能秒解。4.4 一个典型调度结果的长相以冬季典型日为例热负荷在早上7点出现早高峰晚上18点到22点出现晚高峰。设置分时电价谷电0.3元/kWh、平电0.6元/kWh、峰电1.1元/kWh。不引入虚拟储能的调度方案燃气锅炉会在早晚高峰全出力运行电锅炉基本只在谷电时段运行热源出力曲线跟热负荷曲线严丝合缝地“贴”在一起。引入虚拟储能后调度结果发生了两个明显变化第一热源出力曲线变得平缓。早高峰来临前一两个小时燃气锅炉就开始缓慢加出力把管网温度提前“顶”上去一部分用户侧实际感受到的升温时间并没有延迟但热源的爬坡速率大大降低了。第二电锅炉的谷电利用更充分。夜间谷电时段电锅炉不仅满足即时热负荷还额外把管网温度往上抬相当于把热量“存”在管网里白天峰电时段再通过管网自然降温把热量放出来。这一充一放等于用便宜的谷电替代了昂贵的峰电供热运行成本下降幅度在15%到25%之间具体取决于电价差和管网规模。输出结果时我会画三张图第一张是热源各设备出力曲线第二张是管网供回水温度动态曲线第三张是虚拟储能SOC曲线。这三张图放在一起整个系统“什么时候充热、什么时候放热、什么时候直接供热”的逻辑一目了然。5. 我踩过的坑和调参经验步长、初值、验证缺一不可5.1 时间步长和空间步长怎么配对才不会抖显式迎风差分格式有个典型症状当空间步长不变你一味减小时间步长CFL条件明明满足得很好但温度曲线反而出现锯齿状波动。这不是发散是数值耗散和数值频散在作祟。我自己调试时踩过一次很深的坑。管道流速1 m/s我把Δx取到2米CFL条件要求Δt ≤ 2秒实际取1.5秒照理说很安全。结果一跑管道出口温度曲线在升温阶段出现了明显的高频振荡。排查了很久最后发现问题出在空间步长与管道总长的比例关系上2米步长让1000米的管道分成500段边界条件的数值反射在细网格下反而更容易被激发。最终调整方案是不要盲目追求空间分辨率。管道分段数控制在20到100段之间然后用CFL条件倒推时间步长。以1000米管道为例Δx取10米100段流速1 m/sΔt取9秒CFL 0.9既稳定又不会出现数值振荡。仿真结果跟20米步长对比出口温度误差不超过0.3℃但计算时间节省了将近一半。5.2 初值给不好前十几小时全是“假动态”热网动态仿真最容易被忽视的就是初值。直接给T_init 90℃常数初值跑前几个小时仿真结果跟实际情况根本对不上。因为真实管网在运行中从热源到末端本来就有自然温降越靠近用户侧温度越低入口90℃末端可能只有65℃。用一个平的初值相当于人为给管网“充”了一波热量仿真初期管网会把这部分虚假热量慢慢放出来温度曲线看起来像在下降实际上是数值过程在自愈。我的做法先用恒定入口温度取当天预测的平均供水温度把管网模型单独跑24小时达到准稳态然后把这个准稳态温度场作为调度仿真的初始值。这样初始误差几乎可以完全消除。如果手头有历史运行数据直接采用历史对应时刻的温度场作为初值更准。5.3 虚拟储能SOC的初值必须跟管网热状态对齐这是调度模型和热网动态模型耦合时最隐蔽的一个坑。调度模型里的虚拟储能它的SOC对应的是管网当前温度相对基准温度的偏差。如果SOC初值偏低而管网实际温度很高优化器就会倾向于在调度周期里“赚”这波热量导致供热不足、用户侧温度跌出约束下限。解决办法是在优化循环启动前先用实时温度数据计算虚拟储能的初始SOCE_vess_init Σ M_i c_p (T_i(0) - T_base)把E_vess_init作为优化问题的初始状态约束写进去。这一步做对了调度结果才是在真实物理状态基础上的最优解而不是在一个虚构的“空电池”状态上做规划。5.4 验证模型的两板斧阶跃响应和能量守恒很多人在Matlab里跑通代码、画出曲线就觉得完事了。模型验证这一步反而是整个项目里最不该省的部分。我的习惯是做两个检查。第一是阶跃响应检查。给入口温度施加一个从80℃到90℃的阶跃观察出口温度的响应曲线。管道长度1000米、流速0.5 m/s时理论传输延迟应该是2000秒约33分钟然后出口温度应该以一个接近一阶惯性的曲线逐渐逼近90℃。如果数值仿真出来的延迟时间跟理论值对不上或者稳态值有偏差说明方程里对流项或者散热损失项的系数有问题。第二是能量守恒检查。统计一个完整调度周期内热源总供热量、用户侧总用热量、管网散热损失、虚拟储能终始状态变化量它们应该满足Q_supply Q_demand Q_loss ΔE_vess这个等式左右两边的偏差如果超过2%说明某些环节的能量计算有疏漏。我在实际项目中发现散热损失项的计算是最容易出偏差的因为K_loss同时包含了保温层传导和管道埋深处土壤的传热不同季节、不同土壤湿度下数值会变。我的处理方式是用夏冬两季的实际运行数据分别标定冬季和夏季的K_loss值调度模型中按季节切换精度明显提升。5.5 一套可复用的调试顺序最后分享一套我经过多次项目验证的调试顺序按这个顺序来可以少走很多弯路。第一步先用单根管道验证动态求解器。固定入口温度、固定流量跑一个简单算例对照解析解或阶跃响应理论值。第二步再扩展到一个小型放射状热网两三根管道、一个热源、几个负荷节点验证拓扑连接的边界条件传递是否正确。第三步接着把虚拟储能模块接入单独验证虚拟储能量的计算是否跟管网平均温度的变化一致。第四步最后才把优化调度模型接进来先跑无分时电价工况确认调度结果就是一个满足热平衡的可行解再逐步加入分时电价、机组启停等复杂约束。每步都验证通过再进下一步比直接搭一个完整系统再回头debug要高效得多。这跟我早期直接把所有模块一股脑写完、结果花了一周时间定位一个边界条件传参错误相比效率差距是数量级的。从热网动态模型构建、有限差分求解、虚拟储能量化到调度优化和Matlab实现整条技术链路到这里就完整了。这套方法本身有很好的扩展性比如把建筑围护结构热惯性也纳入虚拟储能范畴、在管网模型中增加多个分支节点的压力流量耦合、或者把调度模型从确定性规划扩展到考虑负荷预测不确定性的鲁棒优化。不过那都是后续延伸的话题了。就当前这套框架来说它的价值在于用一套不算复杂的数学模型把热网从调度问题里的“静态负担”变成了“动态调节资源”让综合能源系统的热侧真正具备跟电侧同等的调度灵活性。