简介面向电力系统调度与优化研究人员提供基于混合整数线性规划的机组组合优化完整实现方案覆盖MILP建模、YALMIP建模与CPLEX求解等核心环节。压缩包内共7个文件大小仅267KB包含MATLAB主程序.m用于模型搭建与求解3个Excel表格.xls存放热备用0.05和0.2状态下的求解结果2个Visio图表.vsdx展示不同热备用场景各时段最优出力曲线另有Word说明文档.docx交代基本要求。文件类型覆盖源码、数据、图表与说明结构清晰。已有3051人学习下载。通过研读代码与结果可掌握将离散开停机变量和连续出力变量统一纳入MILP框架的方法理解目标函数、约束条件及求解器参数设置并快速迁移到实际电力系统调度场景中。1. 电力系统机组组合问题的本质混合整数线性规划为什么能解它调度员手里最难的从来不是“今天要发多少电”而是“早上八点该把哪几台机组转起来”。火电机组从冷态到并网要几十个小时提前启动意味着在没有负荷的时候就开始烧煤启动晚了高峰负荷来了却顶不上去。这就是机组组合Unit Commitment, UC问题在 24 小时甚至更长的调度周期内决定每台机组在哪些时段开机、哪些时段停机再分配各时段的有功出力让总成本最低同时满足负荷、备用、爬坡、最小启停时间等约束。真正麻烦的地方在于机组的开停机状态是一个“开/关”的离散决策而每台机组的出力是一个连续变量。两类变量搅在一起问题就变成了混合整数线性规划MILP。商业求解器Gurobi、CPLEX和开源求解器CBC、HiGHS都内置了分支定界、割平面这些专门对付整数变量的算法所以今天几乎所有的 UC 工程实现都落到 MILP 上把调度规则写成线性约束把成本写成线性目标剩下的交给求解器去搜索整数组合。这套思路对电力系统调度、电力市场出清、微电网能量管理都适用。如果你正在做新能源消纳、储能调度或者多能互补优化UC 的 MILP 建模是绕不开的基线方案。下文把最常用的机组组合 MILP 建模和实现路径完整过一遍从约束设计到求解器参数再到验证结果该看什么。搜索“混合整数线性规划算法”的时候通常关心的就是变量怎么定、约束怎么列、参数怎么调——这篇就按这个顺序展开。2. 机组组合混合整数线性规划建模先立变量再列约束2.1 决策变量、参数表和成本函数UC 问题的输入是预测负荷曲线、机组技术参数、燃料成本输出是每个时段的机组启停计划和出力分配。建模第一步是定义决策变量这是后续所有约束和代码的地基。一个最小可用的 UC 模型只需要下面两组变量变量类型含义(u_{t,i})二元变量1/0机组 i 在时段 t 是否开机(p_{t,i})连续变量可以设下限机组 i 在时段 t 的有功出力MW目标函数是总运行成本最小化。只有开机状态的机组会产生燃料成本燃料成本对出力近似一条二次曲线但在线性规划框架里只能分段线性化或者直接用线性近似。更标准的做法是引入启动成本和空载成本。$$ \min \sum_{t1}^{T} \sum_{i1}^{N} \left[ f_i(p_{t,i}) C_i^u \cdot v_{t,i} \right] $$其中 (f_i(p)) 是出力-成本函数(C_i^u) 是启动成本。这里 (v_{t,i}) 需要额外引入“启动动作变量”来正确计费直接使用 (u_{t,i}) 会多算稳态运行成本但很多简化模型里就用 (u_{t,i}) 近似这在快速复用场景下是可接受的。随后在实现中会采用更精确的启动变量写法。2.2 系统约束负荷平衡、备用容量、线路潮流系统约束是每一层 UC 模型都有的“硬约束”。负荷平衡约束最直接任意时段所有机组出力之和必须等于该时段的负荷需求 (D_t)用式 (1) 表示。这是等式约束如果负荷和机组不可调比如存在强制运行的机组可能会造成无解这时候要引入松弛变量来软化约束——这是实际工程里很常见的技巧。备用容量约束要求系统有足够的可调容量应对负荷预测误差和机组跳闸。表述为任意时段开机机组的最大出力之和必须大于负荷加备用需求。即$$ \sum_i p_{t,i}^{\max} \cdot u_{t,i} \geq D_t R_t $$如果电网内部存在输电阻塞还要加入线路潮流约束。在直流潮流假设下支路潮流是节点注入功率的线性函数因此约束仍然是线性的(-\text{limit}l \leq \sum_i GSF{l,i} \cdot p_{t,i} \leq \text{limit}_l)。这里 GSF 是发电转移因子矩阵由电网拓扑和线路阻抗决定。这个约束不算 UC 的必选项但一旦系统存在阻塞忽略它得出的调度计划在物理上就无法执行。运行备用约束和潮流约束都必须线性化因为 MILP 的求解依赖线性松弛的凸性。非线性约束会让分支定界的下界估计失真求解效率急剧下降。这也是 UC 必须用 MILP 而不是直接用非线性模型的一个核心原因——保留线性结构才能用大规模整数规划算法。2.3 机组约束出力上下限、爬坡、最小启停时间每台机组的技术约束是让 UC 模型区别于单纯经济调度的关键。出力上下限非常直接(p_{t,i}^{\min} \cdot u_{t,i} \leq p_{t,i} \leq p_{t,i}^{\max} \cdot u_{t,i})。注意这里上下限乘了开停机变量 (u)所以停机时机组出力强制为 0开机时出力被限制在技术出力范围内。注意如果这么做当机组在启动或停机过程中出力不能瞬间跨越整个区间需要通过爬坡约束来表达。爬坡约束限制出力从一个时段到下一个时段的变化速度出力上升不能超过爬坡上限 (RU_i)下降不能超过 (RD_i)。数学上要写成同时约束正向和负向变化的两行不等式$$ p_{t,i} - p_{t-1,i} \leq RU_i \cdot u_{t-1,i} p_{t,i}^{\max} \cdot (1 - u_{t-1,i}) $$这个写法的含义是如果上一时段机组已开机则正常限制爬坡速率如果上一时段机组是停机状态则这个约束不起限制作用否则从 0 爬升到最低出力会违反约束。这一行是实际建模中容易踩坑的地方很多刚上手的人漏掉 ((1-u_{t-1,i})) 这一项导致模型几乎无解。最小启停时间是让求解变的困难的主要因素。最小运行时间 (UT_i) 表示机组一旦开机至少持续运行 (UT_i) 个时段才能停机最小停机时间 (DT_i) 类似。这类约束在数学上需要引入额外的启动/停机指示变量同时需要处理调度周期边界的初始化状态。常见的线性写法如下其中 (y_{t,i}) 和 (z_{t,i}) 分别是启动和停机动作变量。$$ u_{t,i} - u_{t-1,i} y_{t,i} - z_{t,i} $$$$ \sum_{k\max(1, t-UT_i1)}^{t} y_{k,i} \leq u_{t,i} $$$$ \sum_{k\max(1, t-DT_i1)}^{t} z_{k,i} \leq 1 - u_{t,i} $$逻辑说明第一个式子是状态转移方程描述开停机状态与动作变量的关系状态上跳说明启动状态下降说明停机。第二式强制在任意时刻的往前 (UT_i) 个时段窗口内若本时段为开机状态则窗口内的启动动作不能超过 1 次——这在逻辑上防止了启动后立刻停机。第三式同理约束停机状态。2.4 冷启动和热启动的决策变量设计机组启动有热启动与冷启动之别。停机时间短锅炉仍有余温启动成本低停机时间长了之后机组需要重新加热启动成本会显著上升。这个问题可以通过分段定义启动成本来处理把启动成本写成停机时长的阶梯函数。数学上MILP 可以用多个启动变量表示不同温度段的启停成本。简化模型常把启动成本简化为固定常数这在燃气轮机场景可以接受但燃煤机组必须在模型中区分热态和冷态否则启动成本估算偏差可达 3 到 5 倍。工程实践中我会优先采用“状态转移 最小启停约束”的建模方式再把冷热启动成本用二元变量分段表达。这样模型还是纯线性约束但决策变量会比最小形式多出一组“冷启动判别变量”。这类模型的求解规模大约在每台机组每天增加 24 到 48 个二元变量对当前求解器来说几乎无感。约束类型典型写法常见误用负荷平衡等式忘记加入松弛变量导致无解备用不等式备用在日前计划里计入但未留够爬坡空间最小运行时间窗口内启动次数 ≤ 1只用运行累计时间约束导致线性化困难冷启动分段常数成本用单一平均值导致计划偏差过大3. 用 Python Pyomo 实现混合整数线性规划的机组组合求解3.1 最小算例的数据准备建模层面已经清晰现在用 Python 的 Pyomo 包做一个最小可运行的 UC 求解器。Pyomo 是开源优化建模语言语法接近数学公式写起来快而且同一套模型可以无缝切换到 Gurobi、CBC、HiGHS 等求解器后端。假设有 4 台机组、24 个时段。机组参数用一个 Python dict 表示。import pandas as pd from pyomo.environ import ( ConcreteModel, Var, Binary, NonNegativeReals, Objective, Constraint, minimize, SolverFactory ) # 机组参数max出力(MW)min出力(MW)爬坡速率(MW/h) # 最小运行/停机时间(h)启动成本($)空载成本($/h) gen_params { G1: {pmax: 300, pmin: 60, ramp: 120, min_up: 4, min_down: 3, startup_cost: 500, no_load: 1200}, G2: {pmax: 250, pmin: 50, ramp: 100, min_up: 3, min_down: 2, startup_cost: 400, no_load: 1000}, G3: {pmax: 180, pmin: 40, ramp: 80, min_up: 2, min_down: 2, startup_cost: 300, no_load: 800}, G4: {pmax: 120, pmin: 30, ramp: 60, min_up: 1, min_down: 1, startup_cost: 200, no_load: 500}, } # 24小时负荷曲线单位 MW load [280, 260, 250, 255, 270, 300, 360, 420, 480, 520, 550, 560, 540, 530, 545, 555, 560, 550, 530, 500, 460, 410, 350, 300] # 备用容量需求取每时段负荷的 5%最小不低于 20 MW reserve [max(20, l * 0.05) for l in load]注意燃料成本一般用二次曲线表示这里做了一阶线性近似空载成本代表最低负荷对应的燃料费用。机组 G1 是最贵的大机组G4 是最便宜的小机组在真实的日前调度中半夜负荷低时大机组不应满发这就是 UC 要算出来的结果。3.2 构建模型并求解状态变量与启动变量的配合T 24 G list(gen_params.keys()) U range(T) model ConcreteModel() # 决策变量运行状态二元、出力连续 model.u Var(U, G, withinBinary) model.p Var(U, G, withinNonNegativeReals) # 启动/停机动作变量二元用于精确计费与最小启停时间约束 model.v_start Var(U, G, withinBinary) model.v_stop Var(U, G, withinBinary) # 目标启动成本 运行成本 def obj_rule(m): startup sum(m.v_start[t, g] * gen_params[g][startup_cost] for t in U for g in G) fuel sum(m.p[t, g] * 10 gen_params[g][no_load] * m.u[t, g] for t in U for g in G) return startup fuel model.obj Objective(ruleobj_rule, senseminimize) # 约束1负荷平衡每时段 def balance_rule(m, t): return sum(m.p[t, g] for g in G) load[t] model.balance Constraint(U, rulebalance_rule) # 约束2备用容量 def reserve_rule(m, t): return sum(m.u[t, g] * gen_params[g][pmax] for g in G) load[t] reserve[t] model.reserve Constraint(U, rulereserve_rule) # 约束3出力上下限 def p_limits(m, t, g): return (gen_params[g][pmin] * m.u[t, g] m.p[t, g] gen_params[g][pmax] * m.u[t, g]) model.plim Constraint(U, G, rulep_limits) # 约束4爬坡约束 def ramp_rule(m, t, g): if t 0: return Constraint.Skip pmax, ramp gen_params[g][pmax], gen_params[g][ramp] return (m.p[t, g] - m.p[t-1, g] ramp * m.u[t-1, g] pmax * (1 - m.u[t-1, g])) model.ramp Constraint(U, G, ruleramp_rule) # 约束5状态转移与最小启停 def transition_rule(m, t, g): if t 0: return Constraint.Skip return m.u[t, g] - m.u[t-1, g] m.v_start[t, g] - m.v_stop[t, g] model.trans Constraint(U, G, ruletransition_rule) def min_up_rule(m, t, g): if t 0: return Constraint.Skip min_up gen_params[g][min_up] # 窗口内启动次数不能超过当前开机状态之和 return sum(m.v_start[k, g] for k in range(max(0, t-min_up1), t1)) m.u[t, g] model.min_up Constraint(U, G, rulemin_up_rule) solver SolverFactory(gurobi) result solver.solve(model, teeFalse)逻辑说明目标函数里的v_start精确捕获启动动作u变量只在机组运行时计空载成本balance_rule保证任意时刻总出力等于负荷这个等式约束是所有其他约束的基准plim把出力钳制在技术范围内并通过乘以u让停机时出力归零ramp_rule中用(1-u)项解除停机边界约束是这段代码里最容易出错也最关键的一行。最小启停约束用了“窗口内启动次数 ≤ 当前运行状态”的写法——当窗口时段内有启动动作但此刻处于停机状态时约束直接被违反这就阻止了机组启动后立刻停机的方案。3.3 结果解析取出启动计划和出力曲线求解完成后需要把结果转成可读的调度计划。# 取出每时段每组状态与出力 schedule [] for t in U: row {hour: t1} for g in G: row[f{g}_u] int(model.u[t, g].value) row[f{g}_p] round(model.p[t, g].value, 1) schedule.append(row) df pd.DataFrame(schedule) print(df[[hour, G1_u, G1_p, G2_u, G2_p, G3_u, G3_p, G4_u, G4_p]].to_string()) # 统计启动次数 starts {(t1, g): int(model.v_start[t, g].value) for t in U for g in G if model.v_start[t, g].value 0.5} print(启动事件:, starts) # 验证完全时段出力等于负荷 p_sum df[[c for c in df.columns if _p in c]].sum(axis1) assert max(abs(p_sum - load)) 1e-6, 出力不平衡逻辑说明逐时段打印机组开停机状态和出力值能直接看出半夜哪些机组停了、早高峰哪些机组转起来starts从启动变量里取所有非零启动事件用于核对启动成本和人工检查机组最小运行时间是否满足最后的assert是一个简单的可行性验证确认求解器返回的解没有违反负荷平衡——这一步是每个工程落地前必须做的检查算力越大的模型越不能跳过。3.4 用开源求解器 CBC 替换 Gurobi 的成本对比公司没有商业求解器授权时可以直接把SolverFactory(gurobi)换成SolverFactory(cbc)或SolverFactory(appsi_highs)。HiGHS 是新近开源的高性能求解器对中等规模 UC数百台机组、96 时段以内表现不错。代价是求解时间变长4 台机组 24 时段用例Gurobi 秒级求解CBC 可能要十几秒到几十秒HiGHS 介于两者之间。同样的模型代码更换求解器后要重新验证上下界差MIP gap因为不同求解器内部的分支策略和割平面机制不一样最终找到的整数解可能不同。一个稳妥的流程是先用 CBC 或 HiGHS 快速得到可行解做流程验证再用 Gurobi 做最终的计划文件下发。如果不换求解器只改solver.options调参数也能提高性能下一章展开讲参数设置。4. 混合整数线性规划求解器参数调优与收敛性排查4.1 三个刚需参数MIPGap、TimeLimit、Threads无论规模大小生产环境里求解器都不会放任默认参数跑。UC 模型一般要求按既定的时间窗出结果比如日前调度必须在 20 分钟内完成计算因此参数设置的目标不是“找到全局最优”而是在给定时间内找到工程可用的次优解。参数名作用域含义典型取值MIPGapGurobi / CPLEX / HiGHS当前最优整数解与最优松弛下界的相对差0.011%生产0.00010.01%研究TimeLimit全局求解器运行时间上限秒600 日前调度300 日内滚动Threads全局并行计算线程数物理核数设置为 4 / 8NodeLimit全局最大分支节点数50000应急设置参数设置代码如下solver SolverFactory(gurobi) solver.options[MIPGap] 0.01 # 允许 1% 的最优性差距 solver.options[TimeLimit] 600 # 10 分钟后停止给出当前最优整数解 solver.options[Threads] 8 result solver.solve(model, teeTrue)逻辑说明MIPGap是最重要的生产参数它表示当前找到的可行解与最优解候选下界之间收敛得有多近1% 表示最多差 1%对 UC 来说完全够用TimeLimit防止模型陷入“差一个分支但下界就是上不去”的僵局teeTrue会打印求解日志能实时看到上界和下界的收敛情况。4.2 求解日志怎么看MIP 上界、下界和 gap 的收敛Gurobi 求解 UC 模型时日志里每行会显示Incumbent、BestBd、Gap三列。Incumbent是当前找到的最优整数解的客观值BestBd是最优松弛问题得到的下界。核心思路是Incumbent 是整数变量的可行解而 BestBd 是线性松弛的连续解两者之间的差距才算 MIP gap。Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time 0 0 61000.0 0 14 61000.0000 61000.0000 0.00% - 1s逻辑说明如果Gap长时间不缩小说明整数解空间搜索停滞。优先收紧松弛约束而不是盲目扩大NodeLimit。日志中没有可行解Incumbent 保持无限大且下界也很大大概率是约束矛盾导致无解。这时先检查负荷平衡和备用容量约束是否冲突再看爬坡约束的边界条件。4.3 常见无解原因的排查与修复UC 模型出现“infeasible”结果时的定位优先级如下先检查负荷平衡不可满足再看备用约束是否过紧最后查爬坡约束。Gurobi 的computeIIS()可以直接找出一组最小不可行约束子集这是排错的第一步。from pyomo.opt import TerminationCondition if result.solver.termination_condition TerminationCondition.infeasible: # Gurobi 闭源求解器可直接用 IIS 功能定位冲突约束子集 try: model.dual Suffix(directionSuffix.IMPORT) iis solver.options # 需按求解器 API 调用 computeIIS print(计算 IIS 定位不可行约束) except Exception: print(不可行手动检查负荷/备用/爬坡约束) else: print(求解完成, result.solver.termination_condition)如果不可行源是备用容量在早高峰时段所有开机机组的最大出力之和小于“负荷 备用”那就必须在更早的时段提前启动某台机组。最小运行时间约束会阻止这台机组只在高负荷时段开一小会儿——这正是 UC 模型要自动做出的权衡。如果改成MIPGap 0.05后模型突然变得可行且成本上升很多说明上次求出的可行解质量其实很差需要降低 gap 阈值再算不要直接上线。4.4 大规模问题的两个降维手段机组聚类与时段聚合超过 100 台机组、96 个时段时模型可能有数十万个二元变量直接求解会卡在内存和 CPU 上。两个工程常规手段一是把同类型同容量的机组聚合二是将日内调度时段按负荷特征合并。机组聚类是把完全相同的机组合并为一台“聚合机组”把出力上限乘以台数启动成本乘以台数解出聚合机组的启停后再按贪心算法分配回各物理机组。这种方式在日前市场出清中很常见代价是丢失部分启停次数和费用的精确描述。时段聚合是把 24 个时段按负荷曲率压缩为 8 到 12 个代表性时段求解结束后再插值回原 24 时段做可行性校核。这两种手段互补聚类降空间复杂度压缩降时间复杂度。5. 机组组合模型扩展与热启动技巧5.1 最小启停时间的边界处理首时段已有开机状态的建模模型中t0时段跳过了最小启停约束但现实是调度开始时机组已经有运行状态。常见做法是把前 (UT_i-1) 个时段的约束修正为“直到满足剩余运行时间之前不允许停机”。工程实现上我会在数据文件里增加一个字段init_state取值正数表示已连续运行小时数负数表示已停机小时数。假设 G1 在调度开始前已运行 2 小时且最小运行 4 小时那么前 2 个时段内不允许它停机。这个逻辑通过额外约束表达为u[t,g] 1对 ( t UT_i - |init_state_i| ) 的所有时段生效。5.2 冷启动费用的线性化建模冷启动费用分段函数最常见的处理方式利用两台机组启停判别变量把启动成本写成 ( C^{\text{hot}} \cdot y (C^{\text{cold}} - C^{\text{hot}}) \cdot w )其中w为冷启动判别变量并加约束 ( w \leq y ) 和停机持续时间判别式。线性化的核心是不引入非线性项的同时捕捉启动成本随停机时长变化的特性。一个更细化的模型会按停机时长分成三到五段每段对应一个二元变量。这个方法在燃煤机组为主的系统中效果显著能直接反映冷态启动比热态启动多出的燃料与维护成本让 UC 结果在冬季负荷尖峰前更倾向于保持热备用。5.3 把上一轮的整数解作为 MIP start 二次求解机组组合最常用的热启动技巧是传递上一轮调度计划。在新的 96 时段模型里把上一日 96 时段方案平移 24 时段作为初始可行解传入求解器。Gurobi 通过model.solutions.store_to(model)加solver.options[MIPStart] 1实现。这个技巧在滚动调度里能显著加速收敛——因为相邻日的负荷曲线高度相似上一轮的整数解与当前最优解往往只差几台机组的启停时刻。给一个可行的初始整数解之后求解器可以立即开启割平面分支数量大幅减少MIPGap 收敛到 1% 的时间通常可以缩短一半以上。5.4 验证输出检验备用与爬坡约束的可行性最后一步不能省。无论求解器报告多漂亮都要在结果层面再跑一遍自动校核。脚本读取调度计划后按三个规则检查每时段开机容量与负荷加备用之差相邻时段出力差与爬坡能力之差每台机组连续运行时长不低于最小运行时间。这三项是 UC 计划下发给电厂前必须通过的门禁。任何一项失败都先回去检查对应约束在模型里有没有写对而不是直接调求解器参数。把检查脚本固化到调度流程中会比任何求解器日志都可靠。本文还有配套的精品资源点击获取