电力系统热电联产机组调度优化：从以热定电到热电解耦的代码实战

简介这份资源面向电力系统调度、综合能源优化方向的学习者与研究人员聚焦热电联产CHP机组与火电、风电、变相储热协同的经济调度问题提供基于MATLAB的建模与求解代码帮助读者理解如何在机组启停、负荷需求与排放约束下求取最优运行策略。资源包共8个文件约62KB以4个.m脚本为核心配套2个.xlsx初始数据表、1个程序说明.txt及1个相变储热程序压缩包覆盖锅炉、储热与调度主程序等模块便于对照数据复现实验。目前已有573人学习下载。读者可从中获得完整的调度模型实现思路、储热环节的代码组织方式与数据接口设计并借助说明文档快速理清变量定义、函数调用与结果验证流程适合作为课程设计、论文复现或算法改进的参考底稿。1. 电力系统热电联产机组调度优化从“以热定电”到“热电解耦”的代码实战很多做综合能源或者虚拟电厂的朋友第一次拿到热电联产机组的调度优化任务时最容易踩的坑就是直接把“以热定电”的刚性约束写进模型结果求解器要么报 infeasible要么给出的调度方案在现实中根本没法执行。热电联产机组的调度优化核心难点在于它同时产出电和热而热负荷又具有极强的刚性和时段耦合特性。如果你正在做电力系统经济调度、园区级综合能源管控或者手头有一个含 CHP 机组的微网优化项目那这篇笔记就是把你从“公式看得懂、代码跑不通”的状态里拽出来。我会按“模型怎么建 → 代码怎么写 → 参数怎么调 → 坑怎么避”的顺序把热电联产机组调度优化这件事拆成能直接抄作业的工程路径。2. 热电联产机组调度模型从物理约束到数学表达2.1 为什么“以热定电”是调度优化的第一道坎热电联产机组和纯凝火电机组最大的区别在于它的电出力不是独立可调的而是被热出力“绑”住的。常见做法是采用抽凝式机组其可行运行域是一个由多个不等式围成的多边形。很多新手一上来就写P_elec k * P_heat这种线性关系只适用于背压机组抽凝机组的电热关系是一个区域不是一条线。在调度优化里这个区域通常用四个顶点围成的凸多边形来近似。具体来说抽凝机组的电出力上限和下限都随热出力变化。如果你忽略这个区域直接把电、热当成两个独立变量求解器会给出一个数学上最优、物理上荒谬的解比如热出力为零时电出力却拉满。我一般会先把机组的电热特性曲线整理成一组线性不等式约束。假设热出力为 (H)电出力为 (P)那么约束形式通常是电出力上限(P \leq a_1 H b_1)电出力下限(P \geq a_2 H b_2)热出力上下限(H_{min} \leq H \leq H_{max})电出力上下限(P_{min} \leq P \leq P_{max})这些系数 (a_1, b_1, a_2, b_2) 必须从机组的热力特性说明书里查或者通过历史运行数据拟合。拟合的时候注意不要用全工况数据一把梭要按供暖季和非供暖季分开拟合否则误差会大到让调度方案失去意义。2.2 目标函数与约束条件怎么写成代码能用的形式调度优化的目标函数通常是运行成本最小化包括燃料成本、启停成本、外购电成本。对于热电联产机组燃料成本是电、热出力的二元函数。常见做法是把它线性化或者分段线性化。如果你用 Python 的 Pyomo 或者 scipy 做优化分段线性化可以用 SOS2 约束但更简单的做法是直接用凸组合表示。约束条件除了电热特性还要考虑功率平衡电负荷 电需求 各机组电出力之和 外购电热平衡热负荷 各机组热出力之和 热储能 弃热机组爬坡约束(-R_{down} \leq P_t - P_{t-1} \leq R_{up})最小启停时间约束这里有个血泪经验热负荷的功率平衡里弃热变量一定要加而且要给一个很小的惩罚系数。否则当热负荷低于机组最小热出力时模型会直接无解。现实中可以通过散热或者储热罐来“消纳”多余热量代码里不体现这个求解器就翻车。2.3 用 Python Pyomo 搭一个最小可跑模型下面这段代码是一个单台抽凝式热电联产机组、单时段调度的最小示例。你可以直接复制运行只需要改参数。import pyomo.environ as pyo # 创建模型 m pyo.ConcreteModel() # 参数设置 H_max 100.0 # 最大热出力 MW H_min 20.0 # 最小热出力 MW P_max 200.0 # 最大电出力 MW P_min 50.0 # 最小电出力 MW a1, b1 1.5, 30.0 # 电出力上限系数 a2, b2 0.8, 10.0 # 电出力下限系数 c_heat 15.0 # 单位热出力成本 元/MW c_elec 25.0 # 单位电出力成本 元/MW H_demand 60.0 # 热负荷需求 MW P_demand 120.0 # 电负荷需求 MW # 变量 m.H pyo.Var(bounds(H_min, H_max), doc热出力) m.P pyo.Var(bounds(P_min, P_max), doc电出力) m.P_buy pyo.Var(bounds(0, 500), doc外购电) m.H_curtail pyo.Var(bounds(0, 200), doc弃热) # 约束 m.elec_upper pyo.Constraint(exprm.P a1 * m.H b1) m.elec_lower pyo.Constraint(exprm.P a2 * m.H b2) m.heat_balance pyo.Constraint(exprm.H m.H_curtail H_demand) m.power_balance pyo.Constraint(exprm.P m.P_buy P_demand) # 目标函数 m.cost pyo.Objective( exprc_heat * m.H c_elec * m.P 0.01 * m.H_curtail, sensepyo.minimize ) # 求解 solver pyo.SolverFactory(glpk) results solver.solve(m) print(f热出力: {pyo.value(m.H):.2f} MW) print(f电出力: {pyo.value(m.P):.2f} MW) print(f外购电: {pyo.value(m.P_buy):.2f} MW) print(f弃热: {pyo.value(m.H_curtail):.2f} MW) print(f总成本: {pyo.value(m.cost):.2f} 元)这段代码的逻辑说明elec_upper和elec_lower两个约束把电出力限制在由热出力决定的可行域内。heat_balance里加了H_curtail保证热负荷低于最小热出力时模型仍然可行。目标函数里给弃热加了 0.01 的惩罚系数是为了让求解器优先减少弃热但又不至于因为惩罚太大导致数值问题。参数说明a1, b1和a2, b2是电热特性线的斜率和截距必须根据实际机组填写。c_heat和c_elec是成本系数如果机组是燃煤的可以折算成标煤单价。H_demand和P_demand是当前时段的负荷预测值。求解器我用的 GLPK开源免费适合小规模问题。如果规模大换 CBC 或者 Gurobi。3. 多时段调度与热储能把时间耦合约束写对3.1 热储能模型为什么不能照搬电储能电储能通常只考虑充放电功率和容量热储能还要考虑散热损失和温度分层。很多做电力系统出身的人第一次写热储能模型时直接套用电储能公式结果发现储热罐的“自放电”系数设成 0.01 和 0.1 对调度结果影响巨大。热储能的散热损失跟环境温度、保温层厚度、储热介质都有关常见做法是取一个经验值比如每小时 1% 到 3%。在多时段调度里热储能的约束要写成储热状态(S_t (1 - \sigma) S_{t-1} H_{charge,t} \cdot \eta_{charge} - H_{discharge,t} / \eta_{discharge})储热容量(S_{min} \leq S_t \leq S_{max})充放热功率(0 \leq H_{charge,t} \leq H_{charge,max})(0 \leq H_{discharge,t} \leq H_{discharge,max})不能同时充放通常用二进制变量或者互补约束这里有个玄学问题如果调度周期是 24 小时储热罐的初始状态和终止状态如果不约束求解器会倾向于在第一天把储热放空、最后一天充满导致结果不可持续。我一般会强制 (S_{24} S_0)或者至少让 (S_{24} \geq S_0)。3.2 多时段调度的代码实现与滚动优化下面是一个 24 小时调度的代码框架用 Pyomo 的RangeSet来定义时间索引。为了可读性我只写了核心约束实际项目里还要加机组启停和爬坡。import pyomo.environ as pyo import numpy as np T 24 m pyo.ConcreteModel() m.t pyo.RangeSet(0, T-1) # 负荷预测数据示例 H_demand [60, 55, 50, 48, 45, 50, 65, 80, 85, 80, 75, 70, 70, 72, 75, 80, 85, 90, 88, 82, 75, 70, 65, 62] P_demand [120, 115, 110, 105, 100, 110, 130, 150, 160, 155, 145, 140, 140, 142, 145, 150, 160, 170, 165, 155, 145, 135, 128, 122] # 变量 m.H pyo.Var(m.t, bounds(20, 100)) m.P pyo.Var(m.t, bounds(50, 200)) m.P_buy pyo.Var(m.t, bounds(0, 500)) m.H_curtail pyo.Var(m.t, bounds(0, 200)) m.S pyo.Var(m.t, bounds(0, 300)) # 储热状态 m.H_ch pyo.Var(m.t, bounds(0, 50)) # 充热 m.H_dis pyo.Var(m.t, bounds(0, 50)) # 放热 # 电热特性约束 a1, b1 1.5, 30.0 a2, b2 0.8, 10.0 def elec_upper_rule(m, t): return m.P[t] a1 * m.H[t] b1 m.elec_upper pyo.Constraint(m.t, ruleelec_upper_rule) def elec_lower_rule(m, t): return m.P[t] a2 * m.H[t] b2 m.elec_lower pyo.Constraint(m.t, ruleelec_lower_rule) # 热平衡 def heat_balance_rule(m, t): return m.H[t] m.H_dis[t] m.H_curtail[t] H_demand[t] m.H_ch[t] m.heat_balance pyo.Constraint(m.t, ruleheat_balance_rule) # 电平衡 def power_balance_rule(m, t): return m.P[t] m.P_buy[t] P_demand[t] m.power_balance pyo.Constraint(m.t, rulepower_balance_rule) # 储热动态 sigma 0.02 # 散热损失率 eta_ch 0.95 eta_dis 0.95 def storage_rule(m, t): if t 0: return m.S[t] (1 - sigma) * 100 m.H_ch[t] * eta_ch - m.H_dis[t] / eta_dis return m.S[t] (1 - sigma) * m.S[t-1] m.H_ch[t] * eta_ch - m.H_dis[t] / eta_dis m.storage pyo.Constraint(m.t, rulestorage_rule) # 终端约束 m.terminal pyo.Constraint(exprm.S[T-1] 100) # 目标函数 c_heat, c_elec, c_buy 15.0, 25.0, 80.0 def obj_rule(m): return sum(c_heat * m.H[t] c_elec * m.P[t] c_buy * m.P_buy[t] 0.01 * m.H_curtail[t] for t in m.t) m.cost pyo.Objective(ruleobj_rule, sensepyo.minimize) solver pyo.SolverFactory(glpk) results solver.solve(m) for t in m.t: print(f时段{t}: 热{pyo.value(m.H[t]):.1f}, 电{pyo.value(m.P[t]):.1f}, f储热{pyo.value(m.S[t]):.1f}, 弃热{pyo.value(m.H_curtail[t]):.1f})这段代码的关键点储热动态约束里初始状态我设成了 100终端约束要求最后时段储热不低于 100这样避免求解器把储热“吃干榨净”。sigma取 0.02是常见的小型储热罐散热率。eta_ch和eta_dis取 0.95是充放热效率。如果你的储热罐更大或者保温更好这些参数要相应调整。参数说明c_buy是外购电成本通常比自发电成本高设成 80 元/MW 是为了让模型优先自发电。如果外购电有分时电价可以把c_buy改成按时间索引的数组。H_demand和P_demand是 24 小时预测曲线实际项目里要用预测模型输出。3.3 滚动优化与日前调度的衔接日前调度给出的是 24 小时计划但实际运行中负荷预测会偏差机组状态也会变。我一般会在日前调度之后每 15 分钟做一次滚动优化只优化未来 4 小时但只执行第一个时段的指令。滚动优化的代码框架和上面一样只是时间窗口缩短初始状态用实际测量值。这里有个坑滚动优化如果每次都用 GLPK 重新求解计算时间可能跟不上。对于小规模问题GLPK 毫秒级出解没问题。但如果机组数量超过 5 台时段超过 96 个建议换用 Gurobi 或者用热启动。另外滚动优化的目标函数里要加一个对日前计划的跟踪惩罚项否则调度结果会频繁大幅调整对机组不利。4. 代码优化与求解器调参让模型跑得快、解得稳4.1 从 Pyomo 到矩阵形式什么时候该换写法Pyomo 写模型可读性好但构建大规模问题时Python 层的循环会成为瓶颈。如果你的模型有 100 台机组、96 个时段Pyomo 的变量数会超过 10 万构建时间可能比求解时间还长。这时候常见做法是把模型导出成 LP 文件或者直接用 scipy 的linprog或者pulp的矩阵接口。我一般会先用 Pyomo 做原型验证确认模型逻辑正确后再用pyomo.contrib里的LinearExpression或者直接写矩阵。如果问题规模特别大可以考虑用 JuMPJulia或者直接调 Gurobi 的 C API。但对于大多数园区级项目Pyomo GLPK/CBC 足够了。代码优化的另一个方向是减少二进制变量。机组启停约束如果用二进制变量求解时间会指数上升。如果机组启停成本相对于燃料成本很小可以先用连续变量松弛求解后再取整。或者用“最小启停时间”的线性化方法避免引入二进制。4.2 求解器参数怎么调以 GLPK 和 CBC 为例GLPK 和 CBC 都是开源求解器适合中小规模线性规划和混合整数规划。GLPK 的参数可以通过options传给 Pyomosolver pyo.SolverFactory(glpk) solver.options[msg_lev] GLP_MSG_OFF # 关闭求解日志 solver.options[tm_lim] 60000 # 超时 60 秒 results solver.solve(m)CBC 的参数类似solver pyo.SolverFactory(cbc) solver.options[seconds] 60 solver.options[ratioGap] 0.01 # 相对间隙 1% results solver.solve(m)ratioGap是混合整数规划的重要参数。默认是 0求解器会一直找最优解可能跑几个小时。设成 0.01 意味着找到与最优解差距 1% 以内的可行解就停。对于调度问题1% 的间隙通常可以接受因为负荷预测本身就有误差。如果求解器报 infeasible不要急着改模型。先用results.solver.termination_condition看终止条件。如果是infeasible常见原因是热平衡约束太紧或者电热特性约束的可行域为空。我一般会先把弃热变量上限放大或者把热负荷需求调低看模型能不能跑通再逐步收紧。4.3 用热启动和分解算法加速多时段调度多时段调度如果直接求解变量数随时段线性增长。对于 24 时段、单机组变量数几百个没问题。但如果要优化一周 168 个时段或者多机组变量数上万求解时间会明显增加。这时候可以用拉格朗日松弛或者 Benders 分解把问题拆成主问题和子问题。不过对于大多数工程应用更实用的做法是“滚动优化 热启动”。热启动就是把上一次求解的结果作为下一次求解的初始值。Pyomo 支持通过warmstartTrue传给求解器但 GLPK 不支持热启动CBC 和 Gurobi 支持。如果你用 CBC可以在求解前把变量的value设成上次的解for t in m.t: m.H[t].value previous_H[t] m.P[t].value previous_P[t] solver pyo.SolverFactory(cbc) solver.options[warmstart] True results solver.solve(m)热启动对混合整数规划的效果尤其明显因为求解器可以复用之前的分支定界树。但要注意如果模型约束变了热启动的初始值可能不可行求解器会自己调整不会报错。5. 避坑与排查热电联产调度代码的五个常见翻车点5.1 现象求解器报 infeasible但检查约束感觉没问题原因最常见的是热负荷低于机组最小热出力而代码里没有弃热变量或者弃热上限设得太小。另一个原因是电热特性约束的可行域在某个热出力区间内为空比如a1 * H b1 a2 * H b2。解决先检查H_min和H_max之间是否所有点都满足a1 * H b1 a2 * H b2。如果不满足说明特性线参数拟合错了。然后确认弃热变量存在且上限足够大。5.2 现象调度结果里机组电出力频繁在上下限之间跳变原因目标函数里没有考虑机组爬坡成本或者爬坡约束没加。纯线性规划求解器会倾向于让机组在可行域边界上跳因为这样成本最低。解决加爬坡约束并且在目标函数里加一个很小的出力变化惩罚项比如0.001 * (P[t] - P[t-1])**2但平方项会变成二次规划GLPK 不支持。可以用绝对值线性化或者直接加爬坡约束就够了。5.3 现象储热罐在调度周期结束时被放空第二天没法用原因没有加终端约束或者终端约束太松。解决强制S[T-1] S[0]或者设一个最低储热水平。如果项目允许跨天调度可以把调度周期拉长到 48 小时但计算量翻倍。我一般会设S[T-1] 0.5 * S_max保证第二天有足够的热储能。5.4 现象外购电成本设得很高但模型还是大量外购电原因电热特性约束导致机组电出力上限被热出力限制。如果热负荷很低机组电出力上限也低即使外购电贵模型也只能外购。解决检查热负荷预测是否偏低或者考虑热电解耦方案比如加电锅炉或者热泵。如果机组是抽凝式可以适当放宽电出力下限让机组在低热出力时多发电。5.5 现象代码在本地跑得好部署到服务器上结果不一样原因求解器版本不同或者浮点数精度差异。GLPK 和 CBC 在不同平台上的默认参数可能不同。解决在代码里显式设置求解器参数不要依赖默认值。另外把模型导出成 LP 文件在不同平台上用同一个求解器求解对比结果。如果结果差异大检查是不是有变量没有初始化或者约束里有随机数。6. 进阶技巧用灵敏度分析验证调度方案的鲁棒性调度优化给出的是一条确定性的曲线但现实中热负荷和电负荷都有预测误差。我一般会在模型跑通之后做一轮灵敏度分析把热负荷和电负荷分别上下浮动 5%、10%看总成本变化多少看机组出力是否还在可行域内。如果某个负荷水平下模型直接 infeasible说明调度方案没有留足裕度。具体做法是写一个循环每次改负荷数据重新求解记录成本和关键变量。下面是一个简单的灵敏度分析代码片段import copy base_H H_demand.copy() base_P P_demand.copy() results_list [] for delta in [-0.1, -0.05, 0, 0.05, 0.1]: for t in range(T): H_demand[t] base_H[t] * (1 delta) P_demand[t] base_P[t] * (1 delta) # 重新求解模型需要重新构建或者更新参数 # 这里省略模型重建代码实际项目里可以用 Pyomo 的 Param 来更新 results_list.append({ delta: delta, cost: pyo.value(m.cost), feasible: results.solver.termination_condition optimal }) for r in results_list: print(f负荷偏差 {r[delta]*100:.0f}%: 成本 {r[cost]:.0f} 元, 可行: {r[feasible]})这段代码的关键是不要每次重新构建模型而是把负荷定义成pyo.Param然后更新参数值再求解。这样速度快很多。如果某个偏差下模型不可行说明你的弃热或者外购电上限设得太紧需要放宽。另一个进阶技巧是给热负荷和电负荷加一个“鲁棒裕度”比如在约束里写成H H_curtail H_demand * 1.05这样即使预测偏低 5%模型仍然可行。代价是成本会略高但换来的是调度方案的鲁棒性。我一般会在供暖季初期用 5% 裕度中期用 3%末期用 2%根据实际预测误差调整。最后说一个我自己的习惯每次写完调度代码先不要看成本先看弃热和外购电。如果弃热很大说明热储能或者热电解耦没做好如果外购电很大说明机组电出力被热出力绑得太死。这两个指标比总成本更能反映模型的健康度。希望帮到你。本文还有配套的精品资源点击获取