简介面向电力系统调度与优化研究者这份资源覆盖机组组合Unit Commitment这一经典难题聚焦如何在满足负荷平衡、备用容量、启停时间与爬坡等约束下通过混合整数线性规划MILP统一优化机组启停状态与出力水平。资源基于MATLAB平台完整演示借助YALMIP建模并调用Cplex求解的流程适合需要掌握经济调度建模与求解方法的进阶学习者。压缩包内共有7个文件包含MATLAB源码主程序、需求说明文档、Visio图表和Excel数据表呈现热备用0.05与0.2两种场景下的机组最优出力与求解结果整体大小仅267KB。已有3051人学习。通过该资源可快速理解MILP在机组组合中的建模关键获得可直接运行的代码框架与结果数据深入掌握YALMIPCplex的调用流程和不同热备用约束对调度结果的影响。其代码数据结构清晰便于替换负荷曲线与机组参数开展扩展实验为后续研究或课程设计提供良好起点。1. 机组组合优化为什么绕不开混合整数线性规划机组组合Unit CommitmentUC的工业级解法现在基本统一到混合整数线性规划MILP上未来 24 到 96 个时段决定每台机组开还是停、出力多少在满足负荷和备用约束下让总成本最低。难点在“开/停”这种 0/1 决策可行域不是凸集普通线性规划解不了。用二进制变量表达机组状态、把燃料成本分段线性化后剩下的就是让求解器分支定界。为什么是 MILP 而不是启发式算法因为求解器跑混合整数线性规划算法时会同时给出解和 MIP gap出清要对成本负责最优性证明几乎不可替代爬坡、启停时间这类约束写起来也自然后续扩展备用、安全约束都在同一框架里。下面按建模、求解到工程部署展开适合做电力现货出清、新能源调度和园区能源管理的工程师。2. 机组组合的混合整数线性规划建模决策变量与约束线性化建 UC 模型是有固定套路的先定变量再写目标最后把运行约束一条条线性化。大多数失败的模型不是约束漏了而是变量没选对导致启动成本、最小启停时间根本没法表达。这一章从最小变量集讲起给出一套可以直接落笔的建模清单。2.1 最小变量集状态、出力、启动和停机每台机组每个时段至少需要四类变量变量类型含义u_{g,t}Binary1 表示机组 g 在时段 t 开机p_{g,t}Continuous机组 g 在时段 t 的有功出力MWv_{g,t}Binary1 表示时段 t 发生启动动作w_{g,t}Binary1 表示时段 t 发生停机动作很多初学者只写u和p看到启动成本就犯难因为u的差分u_t - u_{t-1}虽然能表达状态变化却不能区分“启动”和“停机”的动作方向。引入v和w之后状态转移用一条等式约束即可锁定u_{g,t} - u_{g,t-1} v_{g,t} - w_{g,t}再补一条v_{g,t} w_{g,t} 1避免同一时段既启动又停机的逻辑矛盾。注意这里t0需要给定初始状态u_{-1}通常取机组的实际运行状态或统一设为 0。2.2 目标函数燃料成本分段线性化与启动成本目标函数一般由三部分组成燃料成本、空载成本和启动成本。燃气或燃煤机组的燃料成本是出力的二次函数a * p^2 b * p c直接写进 MILP 不行最常见的处理是分段线性化把[Pmin, Pmax]切成 K 段每段端点(P_k, C_k)用插值权重表示要求权重和等于u_t同时加入 SOS2 约束或引入 K-1 个二进制变量强制权重相邻段非零。段数取 4 到 6 段时精度和求解速度的平衡最好。启动成本SC_g直接乘v_{g,t}进目标即可minimize sum_g sum_t ( b_g * p_{g,t} c_g * u_{g,t} SC_g * v_{g,t} )如果机组有热态、温态、冷态之分可以把启动成本按停机时长分成两三档用额外二进制变量选择档位。小规模项目固定一个启动成本值已经够用别一上来就上冷热态模型。2.3 最小启停时间、爬坡和备用约束的线性化写法2.3.1 最小运行与最小停机时间最小运行时间约束要求机组一旦启动必须保持运行至少 UT 个时段。标准写法是sum_{it}^{tUT-1} u_{g,i} UT_g * v_{g,t}对每个满足t T - UT_g 1的时段生效。逻辑很直接只要v1后面连续 UT 个时段的开机状态和必须不小于 UT也就是全为 1。最小停机时间对称地写成sum_{it}^{tDT-1} (1 - u_{g,i}) DT_g * w_{g,t}w1时后面 DT 个时段必须全部停机。注意这两条约束对窗口末尾不生效因为tUT-1会超出时段范围所以模型允许在最后几个时段启动而不强制运行满时长。工程上会额外延长调度窗口或加惩罚项处理这个尾部效应。2.3.2 爬坡约束启动和停机当段的特殊处理正常运行跨时段的爬坡约束为p_{g,t} - p_{g,t-1} RU_g * u_{g,t-1} SU_g * v_{g,t}p_{g,t-1} - p_{g,t} RD_g * u_{g,t} SD_g * w_{g,t}这就是 UC 里最容易写错的约束。RU * u_{t-1}表示只有上一时段开机时才有正常爬坡限制SU * v_t是启动当段的额外爬坡空间。机组启动瞬间出力要从 0 跳到至少 Pmin所以SU_g必须不小于Pmin_g否则启动当段直接不可行。停机方向同理SD_g必须能在最后一个运行时段把出力降到 0。有的教材把这组约束简化为统一爬坡率规模小、机组同质时问题不大异构机组一多就会出反直觉的不可行解。出力上下限和系统约束相对简单Pmin_g * u p Pmax_g * u负荷平衡sum p_{g,t} D_t旋转备用sum Pmax_g * u_{g,t} D_t R_t。备用约束是整个问题产生整数耦合的主要来源之一后面求解效率部分还会回来说它。2.4 三机组 24 时段算例参数为了后面跑代码先给一组可直接用的参数。三台机组G1 为重载基荷机组G3 为快速启停的峰荷机组。机组Pmin (MW)Pmax (MW)爬坡 RU/RD (MW/h)启动/停机爬坡 SU/SD (MW)最小开/停 时间 (h)b (元/MWh)c (元/h)启动成本 (元)G11005002002504 / 424060008000G2503001501502 / 326040005000G3301501001001 / 130025003000负荷序列取典型日峰谷形态峰值 785 MW、谷值 380 MW备用统一取 60 MW。这个规模用开源求解器几秒内能出最优解适合验证模型逻辑也适合做后续参数实验的基线。3. 用 Python 和 Pyomo 跑通机组组合最小模型模型建好之后落到代码层需要考虑两件事求解器选哪个、建模语言怎么组织约束。这一章先给选型建议再给一段完整可运行的 Pyomo 代码最后说明怎么读结果。3.1 求解器怎么选开源 CBC/HiGHS 与商业求解器的边界UC 问题的规模和求解器能力直接挂钩。教学和中小规模场景用开源求解器完全够生产出清基本都要商业求解器原因不只是速度还有数值稳定性和技术支持。求解器许可适用场景备注GLPKGPL教学示例几十台机组以内可用MILP 性能一般CBCEclipse中小规模单机本文示例使用安装 coinor-cbc 即可HiGHSMIT中规模近年新增了较完整的 MILP 支持活跃度高SCIP学术免费研究用途约束处理器可定制适合论文复现Gurobi / CPLEX商业生产出清数值稳定分布式并行MIP 启发式强选型有个经验法则如果 30 台机组、96 时段的问题在商业求解器上要跑超过 10 分钟先怀疑模型而不是求解器如果开源求解器 5 分钟出不了可行解再考虑换求解器。换求解器永远比改模型容易但模型本身的整数变量结构决定了性能上限。3.2 用 Pyomo 把三机组模型写成可运行代码Pyomo 是电力行业最常见的建模语言之一模型定义和求解器解耦换求解器只改一行。下面这段代码对应 2.4 节的参数直接保存为uc.py运行即可。# 3 台机组、24 时段的机组组合最小 Pyomo 示例 # 依赖pyomo以及可执行文件在 PATH 中的 cbcapt install coinor-cbc from pyomo.environ import ( ConcreteModel, Set, Param, Var, Objective, Constraint, Binary, NonNegativeReals, minimize, SolverFactory ) # 机组参数pmin/pmax(MW)、爬坡 ru/rd、启动/停机爬坡 su/sd、 # 最小开/停机 ut/dt(h)、燃料斜率 b(元/MWh)、空载成本 c(元/h)、启动成本 sc(元) gens { G1: dict(pmin100, pmax500, ru200, rd200, su250, sd250, ut4, dt4, b240, c6000, sc8000), G2: dict(pmin50, pmax300, ru150, rd150, su150, sd150, ut2, dt3, b260, c4000, sc5000), G3: dict(pmin30, pmax150, ru100, rd100, su100, sd100, ut1, dt1, b300, c2500, sc3000), } T 24 load [420, 405, 390, 380, 395, 430, 505, 620, 685, 715, 735, 745, 705, 695, 725, 745, 765, 785, 745, 705, 685, 625, 545, 465] reserve 60 m ConcreteModel() m.G Set(initializegens.keys()) m.T Set(initializerange(T)) m.u Var(m.G, m.T, withinBinary) # 开机状态 m.v Var(m.G, m.T, withinBinary) # 启动动作 m.w Var(m.G, m.T, withinBinary) # 停机动作 m.p Var(m.G, m.T, withinNonNegativeReals) # 出力 MW # 全部机组初始为停机初始开机时把对应键改为 1 m.u_init Param(m.G, initialize{g: 0 for g in gens}) # 目标燃料成本 空载成本 启动成本 def obj_rule(m): return sum( gens[g][b] * m.p[g, t] gens[g][c] * m.u[g, t] gens[g][sc] * m.v[g, t] for g in m.G for t in m.T ) m.obj Objective(ruleobj_rule, senseminimize) # 状态转移u[t] - u[t-1] v[t] - w[t]t0 时用初始状态 def logic_rule(m, g, t): prev m.u[g, t - 1] if t 0 else m.u_init[g] return m.u[g, t] - prev m.v[g, t] - m.w[g, t] m.logic Constraint(m.G, m.T, rulelogic_rule) # 同一时段不允许同时启动和停机 def one_action_rule(m, g, t): return m.v[g, t] m.w[g, t] 1 m.one_action Constraint(m.G, m.T, ruleone_action_rule) # 出力上下限停机时出力必须为 0 def pmin_rule(m, g, t): return gens[g][pmin] * m.u[g, t] m.p[g, t] def pmax_rule(m, g, t): return m.p[g, t] gens[g][pmax] * m.u[g, t] m.pmin Constraint(m.G, m.T, rulepmin_rule) m.pmax Constraint(m.G, m.T, rulepmax_rule) # 负荷平衡每个时段总出力等于负荷 def balance_rule(m, t): return sum(m.p[g, t] for g in m.G) load[t] m.balance Constraint(m.T, rulebalance_rule) # 旋转备用开机机组的容量之和 负荷 备用 def reserve_rule(m, t): return sum(gens[g][pmax] * m.u[g, t] for g in m.G) load[t] reserve m.reserve Constraint(m.T, rulereserve_rule) # 最小运行时间从启动时刻起连续运行 ut 个时段 def min_up_rule(m, g, t): if t T - gens[g][ut]: return Constraint.Skip # 窗口末尾不约束 return sum(m.u[g, i] for i in range(t, t gens[g][ut])) gens[g][ut] * m.v[g, t] m.min_up Constraint(m.G, m.T, rulemin_up_rule) # 最小停机时间从停机时刻起连续停 dt 个时段 def min_down_rule(m, g, t): if t T - gens[g][dt]: return Constraint.Skip return sum(1 - m.u[g, i] for i in range(t, t gens[g][dt])) gens[g][dt] * m.w[g, t] m.min_down Constraint(m.G, m.T, rulemin_down_rule) # 爬坡约束正常运行段用 ru/rd启动和停机当段放开 su/sd def ramp_up_rule(m, g, t): if t 0: return Constraint.Skip return m.p[g, t] - m.p[g, t - 1] gens[g][ru] * m.u[g, t - 1] gens[g][su] * m.v[g, t] m.ramp_up Constraint(m.G, m.T, ruleramp_up_rule) def ramp_down_rule(m, g, t): if t 0: return Constraint.Skip return m.p[g, t - 1] - m.p[g, t] gens[g][rd] * m.u[g, t] gens[g][sd] * m.w[g, t] m.ramp_down Constraint(m.G, m.T, ruleramp_down_rule) # 求解 solver SolverFactory(cbc) result solver.solve(m, teeFalse) print(求解状态:, result.solver.status, result.solver.termination_condition) print(最优成本(元):, round(m.obj(), 1)) # 输出每个时段的开机组合和出力 for t in range(T): line ft{t:2d} load{load[t]:3d} | .join( f{g}:{开 if m.u[g, t].value 0.5 else 停} {m.p[g, t].value:5.1f}MW for g in m.G ) print(line)这段代码把 2.3 节的所有约束原样翻译成了 Pyomo 表达式。u_init参数处理 t0 时的状态衔接避免第一时段启动不付启动成本的漏洞one_action_rule虽然理论上有转移方程后冗余但保留它可以让求解器在剪枝时更快排除无意义分支。最小启停约束里的Constraint.Skip对应 2.3.1 节说的窗口尾部问题结果里最后几个时段的启动需要人工复核。3.3 代码逻辑拆解与结果怎么读运行后重点看三件事。第一每个时段开的机组数和出力之和是否等于负荷备用约束是否起作用——低谷时段 G3 大概率停机高峰时段三台全开。第二G1 因为启动成本高、燃料便宜通常全程开机或只在凌晨短停如果看到 G3 频繁启停检查是不是sc设置得太低导致模型用启动成本换燃料成本。第三求解器输出的termination_condition必须是optimal如果出现infeasible优先检查爬坡约束里的su/sd是否小于对应机组的pmin。这段代码还能直接改造成实验台把c和sc各乘以一个系数观察最小启停约束如何改变开机组合或者把reserve从 60 调到 150看哪些时段因为备用不足被迫增开机组。这些实验比空读理论更能建立对 UC 结构的直觉。4. 混合整数线性规划算法的求解效率松弛、参数与热启动模型写对只是第一步真正让混合整数线性规划算法在电力系统里跑起来的是求解效率。UC 的整数变量规模随着机组数和时段数线性增长30 台机组 96 时段的模型就有几千个二进制变量不加处理直接丢给求解器往往在给你最优解之前先耗尽预算。这一章从三个方面把求解速度提上去。4.1 线性松弛解为什么是分数分支定界的起点求解器处理 MILP 的第一步是解线性松弛也就是把所有二进制变量u、v、w放宽到[0,1]连续区间。松弛解给出了问题的下界最小化问题但它天然是分数解。举个典型例子三台相同的机组容量都是 300 MW负荷 400 MW空载成本 6000 元线性燃料成本相同。整数最优解是开两台每台出力 200 MW成本 400b 12000 元。但松弛解会把负荷平摊到三台机上每台 133.3 MW开机变量被放松到约 0.44成本变成 400b 8000 元——比真实整数解少了 4000 元因为空载成本只付了 44%。这就是分支定界必须存在的原因以松弛解为根节点不断对分数变量分支成 0 或 1同时用松弛下界剪掉不可能优于当前最优解的枝。UC 问题里最伤求解性能的就是对称性多台同参数机组会让分支树出现大量等价的分数节点剪枝效率骤降。解法之一是加“同质机组按编号顺序开机”的对称消减约束虽然理论上冗余实践中经常把求解时间砍掉一半。4.2 影响分支定界效率的三个参数位置4.2.1 爬坡约束里的启动/停机 M 值2.3.2 节的SU_g * v_{g,t}本质是一个大 M 结构机组启动当段临时放大出力增量上限。SU取太小会错误地切掉可行解取太大则让松弛解的空间虚胖分支时难收敛。工程上的取值规则是SU不小于该机组 Pmin也不超过其最大启动爬坡能力SD同理。不要为了“保险”把 SU/SD 放大到 Pmax 水平那等于告诉松弛解“启动当段可以随便冲”整数解的质量会明显变差。4.2.2 分段线性成本的分段数燃料成本分段越多线性逼近越准但每一段在一个时段里都对应额外二进制变量或 SOS2 变量整数变量数量随之膨胀。按经验24 时段问题用 6 段以内96 时段问题用 4 段以内先跑通再加密验差距。很多生产模型干脆用线性成本误差可以通过重新标定 b、c 参数吸收。记住 UC 的主要矛盾是启停组合的整数决策燃料成本曲线多 0.5% 的精度对组合选择影响很小。4.2.3 网络安全约束中的注入上限 M接入线路潮流约束时要写“机组停机则不注入功率”的形式常见写法是P_g,t M * u_g,t。这里的 M 必须取该机组的 Pmax而不是一个统一的全局大数。M 取 1e6 这种值松弛解会以为停机机组还能注入巨大功率潮流约束形同虚设节点处还会出现矩阵系数范围过大的数值警告。把 M 收紧到物理上界是解决这类数值问题最快的手段。4.3 求解器参数表与热启动写法求解器暴露的参数里对 UC 最敏感的是下面几个。参数名以 CBC 为例Gurobi/CPLEX 有对应写法。参数常用设定作用与理由mipgap0.001 ~ 0.01达到相对间隙即停日前出清 0.1% 足够timelimit300 ~ 600 秒出清有时间窗到点返回当前最优可行解cutsaggressive 或默认割平面加强下界节点数减少但单节点变慢threads4 ~ 8分支定界并行加速不是越大越快presolveon强连通分量分解和约束消除默认开启热启动是滚动调度场景里收益最高的技巧。当天和昨天的机组组合高度相似把上一轮计划作为初值传给求解器可以显著压缩首次可行解的搜索时间。写法是把变量值直接 set 到模型上再求解# 用上一轮计划做热启动减少重复搜索 for g in m.G: for t in m.T: m.u[g, t].set_value(prev_u[t][g]) # prev_u 是上一轮解出的状态表 m.p[g, t].set_value(prev_p[t][g]) # 由状态差分反推启动、停机动作 m.v[g, t].set_value( 1 if t 0 and prev_u[t][g] 0.5 and prev_u[t-1][g] 0.5 else 0 ) m.w[g, t].set_value( 1 if t 0 and prev_u[t][g] 0.5 and prev_u[t-1][g] 0.5 else 0 ) solver.options[mipgap] 0.001 solver.options[timelimit] 300 result solver.solve(m, warmstartTrue, teeFalse)热启动只对可行的初值有效。负荷序列突发跳变、检修计划漏改导致上一轮组合越界时求解器要花额外时间修整初值此时热启动反而更慢。判断依据很简单如果连续几次滚动计算里热启动的效果都在退化先查负荷预测修正量再查初始参数不要盲目怀疑求解器。5. 滚动调度里混合整数线性规划模型的实际调法5.1 时段数、滚动窗口与固定时段96 时段模型每 15 分钟重算一次如果每次都把未来 96 时段全部优化一遍求解器压力大而且计划在相邻轮次之间会抖动。生产系统里更常见的做法是滚动窗口每轮只优化未来 24 到 48 个时段求解完成后只执行前 4 到 8 个时段然后滚动重算。窗口前方固定着上次求解得到的“已下发组合”这些时段里u变量直接赋固定值不进优化窗口末尾可以把最后几个时段的机组状态放宽因为它们只是为下一轮提供热启动初始值。固定时段的实现是把对应u变量临时设为固定参数或者直接把该变量的 bounds 设为相同值。有一处容易踩坑固定标志要在下一轮开始前全部清除。如果“头部固定”的逻辑在多次滚动中叠加固定段会跟着窗口一直往前滑计划会逐渐偏离真实负荷且不易察觉。5.2 出解以后的三层校验脚本求解器返回optimal不代表模型一定正确。解出之后立刻跑三层校验能拦住绝大多数模型低级错误。第一层检查最小启停时间是否被满足第二层检查爬坡约束是否在启停边界上越限第三层对比相邻两轮计划的启停动作次数。下面这段是针对前两层的可复用脚本def check_min_up(u_seq, ut): 校验最小运行时间u_seq 是该机组 T 个时段状态列表 T len(u_seq) starts [t for t in range(T) if u_seq[t] 1 and (t 0 or u_seq[t - 1] 0)] for s in starts: if s ut T: return f窗口末尾启动时段不足, s{s}需人工复核 if sum(u_seq[s:s ut]) ut: return f启动后未满 {ut} 时段, s{s} return ok def check_ramp(p_seq, u_seq, ru, su, sd): 校验爬坡区分正常运行段与启动/停机当段 for t in range(1, len(p_seq)): du p_seq[t] - p_seq[t - 1] if u_seq[t - 1] 1 and u_seq[t] 1 and abs(du) ru: return ft{t} 正常运行爬坡越限: {du}MW if u_seq[t - 1] 0 and u_seq[t] 1 and p_seq[t] su: return ft{t} 启动当段出力越限: {p_seq[t]}MW if u_seq[t - 1] 1 and u_seq[t] 0 and p_seq[t - 1] sd: return ft{t} 停机前出力越限: {p_seq[t - 1]}MW return ok把这两段脚本接在solve之后每次出清自动执行并记录日志。校验函数返回非ok时优先怀疑对应机组的约束生成逻辑而不是先调求解器参数。窗口末尾的“不足”提示属于模型固有尾部效应记录即可不需要修改约束但同一个机组连续多轮都报此类问题说明调度窗口长度设置得太短应把窗口向后延长几个时段。本文还有配套的精品资源点击获取