机组组合的混合整数线性规划建模：0-1变量、约束与Pyomo求解

简介：面向电力系统调度与最优化方向学习者、研究者的机组组合优化资源包，完整演示基于混合整数线性规划（MILP）的机组启停与出力分配建模思路，借助MATLAB、YALMIP与CPLEX实现模型构建与求解，适合电力专业学生、调度工程师快速入门复现。压缩包共7个文件，其中m脚本负责优化模型定义与求解流程，docx文档说明问题基本要求，3个xls分别保存热备用系数0.05与0.2下的组合求解结果，2个vsdx以图表呈现不同备用需求下各时段的机组最优出力，整包仅267KB，小巧而完整。资源已有3051人学习下载，可作为课程设计、科研参考或工程培训素材。通过对照m代码、Excel结果与Visio图表，可直观理解0-1整数变量如何表示机组启停状态、连续变量如何分配出力水平，同时掌握CPLEX分支定界法求解混合整数规划的关键设置与结果分析方法，为后续扩展安全约束机组组合、考虑新能源不确定性等实际问题打下坚实基础。

1. 从凌晨负荷曲线说起：机组组合为什么是混合整数线性规划问题

凌晨三点负荷掉到500多MW，三台机组里留哪台在线、停机那台明天中午还要再启动一次，启动煤耗、寿命损耗都得算进成本。电网调度平时说的机组组合，就是给定负荷预测、备用要求和机组技术参数，先把“哪些机组在哪些时段开机”这个0-1决策定下来，再把有功功率分配到在线机组上。混合整数线性规划恰好把这两层放进同一个优化问题：连续变量负责功率分配，整数变量负责启停，目标是最小化燃料、启动和空载成本。它解决的不是“把曲线填平”，而是让调度员拿到一个带MIP gap的寻优结果，而不是靠候选顺序拍脑袋。这篇笔记适合正在做日前/日内计划、或者想把机组组合换成数学规划选型的人。

2. 机组组合建模：目标函数、约束与0-1变量如何一层层长出来

2.1 为什么纯线性规划不够：整数变量来自三个现实约束

很多刚接触机组组合的人第一反应是：负荷平衡、燃料成本都是线性关系，用线性规划不就行了？对，如果所有机组默认开机，经济调度确实是个线性规划问题。但实际调度不可能让机组随便开。第一，机组有最小技术出力，火电机组低于某个出力就得停机，这个“要么不开，一开就至少带多少负荷”的语义天然是离散的。第二，机组有最小连续运行时间和最小连续停机时间，昨天刚停了一台机组，今天就让它再启动，不仅违反运行规程，还会让设备寿命加速折旧。第三，旋转备用需要在线机组留出可调容量，哪些机组在线决定了备用够不够，这又是一个0-1判断。

这三种约束都没法用“连续变量构造一个很大的惩罚因子”来糊弄。你可以把0-1变量松弛成0到1之间的连续变量，得到的结果可能是某台机组开0.3台，这在物理上没有意义。所以这里必须用混合整数线性规划。把0-1变量和连续变量放进同一个线性模型里，用分支定界法去搜，才是机组组合的标准做法。

2.2 目标函数怎么铺：空载、边际燃料与启动成本不能漏项

机组组合的目标函数一般是最小化总运行成本。我一般把成本拆成三块：空载成本、边际燃料成本和启动成本。空载成本是机组只要并网就产生的固定费用，单位是元/小时，乘以0-1启停变量；边际燃料成本是每多发1MWh多消耗的燃料费用，单位是元/MWh，乘以连续出力变量；启动成本是机组从停机到并网的一次性费用，单位是元/次，乘启动事件变量。写成紧凑形式：

min Σ_t Σ_i ( N_i · u_it + M_i · p_it + S_i · v_it )

其中 u_it 是机组i在时段t的运行状态，p_it是出力，v_it是启动事件。注意 v_it 不能直接用 u_it - u_i,t-1 塞进目标函数，因为差值是可能为负的，需要在约束里单独定义。漏掉空载成本会让模型倾向于让大机组空转开出力下限来“占位”，漏掉启动成本会让模型频繁启停，看起来燃料成本很低，实际账完全不对。

2.3 约束条件逐条写：从负荷平衡到最小启停时间

一个能跑的机组组合模型，至少要有下面这张约束表：

爬坡约束里的 M 是松弛项，我用 (2-u_it-u_i,t-1)·Pmax_i 来替代大M。当机组在相邻时段都在线时，2-u_it-u_i,t-1=0，爬坡限制严格生效；只要有一侧不在线，后面那项就自动放大到机组容量，相当于放过启动/停机瞬间的出力跳变。最小启停时间的约束写法很关键，启动事件v和停机事件w如果不同时做互斥约束，求解器可能把两个事件同时置1，让模型进入自相矛盾的状态。这段后面会专门讲坑。

3. 用一套6机24时段最小数据在Pyomo里把模型建起来

3.1 测试数据先落地：机组参数和负荷曲线怎么配

做机组组合不要一上来就接几百台机组的真实数据，先从6台机组、24个时段开始。这个规模下，Gurobi基本秒解，你可以快速验证约束写得对不对。我用的演示参数如下：

# 机组参数：pmin 最小出力(MW)，pmax 最大出力(MW)，mc 边际燃料成本(元/MWh) # no_load 空载成本(元/h)，startup 单次启动成本(元)， # min_up/min_down 最小连续运行/停机时间(h)，ramp_up/ramp_down 爬坡速率(MW/h) pmin = [150, 150, 20, 20, 25, 10] pmax = [455, 455, 130, 130, 80, 25] mc = [10, 20, 30, 40, 50, 55] no_load = [1000, 1000, 500, 500, 300, 200] startup = [4500, 5000, 550, 560, 250, 230] min_up = [8, 8, 5, 5, 3, 1] min_down = [8, 8, 5, 5, 3, 1] ramp_up = [180, 180, 80, 80, 60, 30] ramp_down = [180, 180, 80, 80, 60, 30] # 24小时负荷预测曲线(MW)，谷段在凌晨4点，峰段在午间12点 demand = [700, 650, 620, 600, 580, 620, 680, 780, 850, 900, 950, 980, 960, 930, 920, 940, 950, 900, 850, 800, 760, 720, 690, 660] # 旋转备用按当日负荷的10%取，实际调度里由调度机构给定 reserve = [round(d * 0.1) for d in demand]

这里所有成本单位统一成元，你可以直接替换成自己的燃料成本和电价系数。重点检查pmax之和要覆盖峰值负荷加备用，我这里的6台机组总容量1275MW，峰值负荷980MW加10%备用在1078MW附近，留有足够空间。如果你自己的数据里总容量比负荷还低，模型直接不可行，先查这一步。

3.2 决策变量与目标函数：代码先搭骨架

用Pyomo建模时，我习惯先定义三组变量：u是启停状态，p是出力，v和w分别是启动和停机事件。v和w都是0-1变量，不是连续变量，这一步别为了省事把它们定义成NonNegativeReals，否则最小启停时间约束会失效。

import pyomo.environ as pyo m = pyo.ConcreteModel() m.G = pyo.Set(initialize=range(6)) m.T = pyo.Set(initialize=range(24)) # u为启停状态，p为出力，v/w为启动/停机事件 m.u = pyo.Var(m.G, m.T, domain=pyo.Binary) m.p = pyo.Var(m.G, m.T, domain=pyo.NonNegativeReals) m.v = pyo.Var(m.G, m.T, domain=pyo.Binary) m.w = pyo.Var(m.G, m.T, domain=pyo.Binary) def objective_rule(m): return sum( no_load[i] * m.u[i, t] + mc[i] * m.p[i, t] + startup[i] * m.v[i, t] for i in range(6) for t in range(24) ) m.obj = pyo.Objective(rule=objective_rule, sense=pyo.minimize)

目标函数里，no_load乘u，mc乘p，startup乘v。注意启动成本只跟启动事件v挂钩，不是跟“状态变化”的直接差值挂钩，这样每一次冷态启动只计一次费用。如果某台机组初始就在线，t=0不需要启动事件，这个边界会在后面约束里处理。

3.3 约束铺开：负荷平衡、备用、爬坡与最小启停时间

这组约束是模型的正文，也是出问题最多的部分。我把每个约束都用Pyomo规则的写法列出来，命名尽量直白：

# 机组初始状态：0表示初始停机，1表示初始在线 initial_status = [0] * 6 # 负荷平衡：任意时刻所有在线机组出力之和等于负荷预测 def load_balance_rule(m, t): return sum(m.p[i, t] for i in range(6)) == demand[t] m.load_balance = pyo.Constraint(m.T, rule=load_balance_rule) # 旋转备用：在线机组最大可调容量之和要覆盖负荷加备用 def reserve_rule(m, t): return sum(pmax[i] * m.u[i, t] for i in range(6)) >= demand[t] + reserve[t] m.reserve = pyo.Constraint(m.T, rule=reserve_rule) # 出力上下限：u=0时强迫p=0 def p_hi_rule(m, i, t): return m.p[i, t] <= pmax[i] * m.u[i, t] m.p_hi = pyo.Constraint(m.G, m.T, rule=p_hi_rule) def p_lo_rule(m, i, t): return m.p[i, t] >= pmin[i] * m.u[i, t] m.p_lo = pyo.Constraint(m.G, m.T, rule=p_lo_rule) # 启动/停机事件：v和w分别捕捉u从0到1、从1到0的变化 def startup_event_rule(m, i, t): if t == 0: return m.v[i, 0] >= m.u[i, 0] - initial_status[i] return m.v[i, t] >= m.u[i, t] - m.u[i, t - 1] m.startup_event = pyo.Constraint(m.G, m.T, rule=startup_event_rule) def shutdown_event_rule(m, i, t): if t == 0: return m.w[i, 0] >= initial_status[i] - m.u[i, 0] return m.w[i, t] >= m.u[i, t - 1] - m.u[i, t] m.shutdown_event = pyo.Constraint(m.G, m.T, rule=shutdown_event_rule) # 互斥约束：同一时段不能既启动又停机 def startup_shutdown_mutex(m, i, t): return m.v[i, t] + m.w[i, t] <= 1 m.mutex = pyo.Constraint(m.G, m.T, rule=startup_shutdown_mutex) # 爬坡约束：相邻时段出力变化受限，启停瞬间用pmax做松弛 def ramp_up_rule(m, i, t): if t == 0: return pyo.Constraint.Skip return m.p[i, t] - m.p[i, t - 1] <= ramp_up[i] + (2 - m.u[i, t] - m.u[i, t - 1]) * pmax[i] m.ramp_up = pyo.Constraint(m.G, m.T, rule=ramp_up_rule) def ramp_down_rule(m, i, t): if t == 0: return pyo.Constraint.Skip return m.p[i, t - 1] - m.p[i, t] <= ramp_down[i] + (2 - m.u[i, t] - m.u[i, t - 1]) * pmax[i] m.ramp_down = pyo.Constraint(m.G, m.T, rule=ramp_down_rule) # 最小连续运行时间：窗口内的启动事件数不超过当前在线状态 def min_up_rule(m, i, t): lo = max(0, t - min_up[i] + 1) return sum(m.v[i, k] for k in range(lo, t + 1)) <= m.u[i, t] m.min_up = pyo.Constraint(m.G, m.T, rule=min_up_rule) # 最小连续停机时间：窗口内的停机事件数不超过当前停机状态 def min_down_rule(m, i, t): lo = max(0, t - min_down[i] + 1) return sum(m.w[i, k] for k in range(lo, t + 1)) <= 1 - m.u[i, t] m.min_down = pyo.Constraint(m.G, m.T, rule=min_down_rule)

这段代码里最值得看的是 t=0 分支。如果忽略初始状态，模型会默认“第一小时前全部机组停机”，这会让启动事件v在第一个时段失去约束，机组可以白启动而不计启动成本。这里我用 initial_status 把真实初始状态带进 v[0] 和 w[0] 的边界，保证边界不漂。爬坡约束用 (2-u_it-u_i,t-1)*pmax 而不是写死的大M，好处是你不用去调M的值，pmax本身就能起到“不在线就不约束”的作用；缺点是当相邻时段有一侧离线时，约束右端会放大，这是工程近似，不是严格建模。如果你要处理启动过程中的分段出力轨迹，需要另外引入启动轨迹变量，但6机模型用这个松弛就够了。

3.4 求解与结果查看：先确认能跑出可行解

模型建完，剩下的就是用求解器跑。我偏好Gurobi，没有授权就用开源的CBC，Pyomo的接口是统一的。

solver = pyo.SolverFactory('gurobi') solver.options['MIPGap'] = 0.001 solver.options['TimeLimit'] = 120 solver.options['Threads'] = 4 res = solver.solve(m, tee=True) assert res.solver.termination_condition == pyo.TerminationCondition.optimal, f"求解失败: {res.solver.termination_condition}" total_cost = pyo.value(m.obj) print("最优总成本: {:.0f} 元".format(total_cost)) for t in range(24): on = [i for i in range(6) if pyo.value(m.u[i, t]) > 0.5] out = " ".join(f"G{i+1}:{pyo.value(m.p[i, t]):6.1f}MW" for i in on) print(f"t={t:02d} demand={demand[t]:3.0f} reserve={reserve[t]:3.0f} | {out}")

我一般把 MIPGap 设为 0.001，也就是0.1%的相对gap。这个数值对24时段6机组的规模很容易达到，但对几百台机组的大系统不能照抄，后面会讲怎么调。打印结果时不要只打印成本，要把每个时段的在线机组和出力都打出来，肉眼扫一遍：凌晨有没有该停的大机组在线，午间高峰有没有启动数量不够导致备用不足，这些都比看目标值更直观。

4. 求解器参数怎么设：MIP gap、时间限制与热启动技巧

4.1 分支定界在做什么：为什么MIP gap比目标值更值得看

混合整数线性规划算法的求解过程不是一个“黑匣子”一头扎到底。Gurobi这类求解器先把整数变量松弛成连续变量，得到一个松弛解，然后通过分支定界不断把整数解的空间切小。每找到一个可行整数解，它就有一个上界；满载的松弛解是下界。上界和下界之间的相对差距就是MIP gap。所以你看到的“gap 0.1%”的意思是：当前找到的这个整数解，离理论上最优整数解的距离不超过0.1%。这在工程上是可接受的收敛标准，而且这个数字比单纯看“目标值是多少”要诚实得多。很多人的模型跑出几百万成本，但没看gap，可能那个解跟最优解差了几万块，尤其是启停成本大的机组，一个错误的启停决策就能让成本差出一大截。

4.2 三个必调参数：MIPGap、TimeLimit与Threads

求解器参数里，我先动MIPGap、TimeLimit和Threads三个。它们的含义在Gurobi和CBC里基本一致，只是写入方式略有差异：

Gurobi里设置方式是 solver.options['MIPGap'] = 0.001，CBC里一般是 solver.options['mipgap'] = 0.01，具体参数名以你装的版本帮助为准。不要一上来就设 TimeLimit=3600 和 MIPGap=0.0001，整数规划的证明最优在机组组合这种规模下可能跑出几小时，结果就是调度方案出得太慢，反而没人敢用。我的习惯是先放宽到 TimeLimit=120、MIPGap=0.01，拿一个可行解看约束有没有问题，确认模型正确后再收紧gap。MIPFocus 也是有用的优化项：如果你有大量机组需要快速给出一个次优可行方案，MIPFocus=1 会让求解器更早剪枝，虽然牺牲一点gap，但能保证方案先出来。

4.3 热启动：把上一轮解当后悔药塞进去

滚动调度里，今天跑完明天的日前计划，过了四个小时又要重新跑一次修正计划。这时候如果从零开始求解，是对算力的浪费。常见做法是把上一轮解作为MIP start传给新模型，相当于告诉求解器“上个方案长这样，你从这个附近开始搜”。在Pyomo里实现是直接给变量赋初值：

# last_u / last_p 是上一轮模型求解后取出的变量值 for i in range(6): for t in range(24): m.u[i, t].value = float(last_u[i][t]) m.p[i, t].value = float(last_p[i][t])

赋初值不能保证这个解一定可行，比如负荷曲线变了，原来的启停方案可能不满足新的旋转备用要求。但Gurobi会把MIP start作为一条热路径去引导分支，不可行部分会在预处理阶段被拒绝，不会拖累求解太多。实际项目里，热启动往往能让滚动求解时间从几十秒压到几秒。这个技巧在单次离线优化里没用，但在重复跑的批处理场景里很值钱。

5. 机组组合MILP的5个常见坑与排查

5.1 现象：启动和停机事件同时为1

有次我从一份开源模型里扒约束代码，只抄了 v_it >= u_it-u_i,t-1 和 w_it >= u_i,t-1-u_it，没抄互斥约束。求解结果里，某台机组同一个时段既启动又停机，目标值还很漂亮。原因是停机变量w在目标函数里没有惩罚项，模型可以把w随意置1而不付出代价，而w一旦为1，最小连续停机约束就会变松，等于给模型开了一个“逃逸通道”。解决方法是补上 v_it + w_it <= 1，同时把t=0的边界处理好。只要变量w参与约束，就必须给它明确的物理语义，不能只在上限不等式里用。

5.2 现象：爬坡约束让启动瞬间变成不可行

在一个初始状态全停机的案例里，求解器直接报不可行，日志里指向爬坡约束。排查下来发现，模型里写的是 p_it - p_i,t-1 <= RU_i，没有考虑 t-1 时段机组根本没在线。机组从停机到并网，出力从0跳到最小技术出力，这个跳变不是爬坡约束应该管的，但硬写会让启动瞬间被限制在爬坡速率以内，导致模型找不到任何可行解。解决方法是改用带松弛项的写法，把不在线一侧的爬坡约束放开。如果你要更严格的建模，可以引入启动轨迹变量，把机组在启动过程中的分段出力路径单独刻画；但对大多数24时段机组组合，用松弛项是工程惯例。

5.3 现象：备用约束看着满足，实际可用备用不够

某次模型跑通了，检查约束发现旋转备用确实大于等于负荷加备用，但调度员说方案没法执行。后来定位到原因：备用约束只写成了 sum(pmax*u) >= demand + reserve，这相当于只保证“在线机组的容量之和足够”，没有考虑机组爬坡速率能不能在调度时间尺度内把备用顶上去。比如一台爬坡速率只有30MW/h的机组，它在线但没法在10分钟内增加30MW备用容量。解决方法是把备用约束里的可用容量按爬坡速率折算，或者把备用需求拆到“快速响应机组”这一类上。真正要落地的时候，调度机构会给一个明确的时间窗口，你用那个窗口去折算机组能贡献的备用值，而不是拿容量硬扛。

5.4 现象：MIPGap设太松，结果跑两次差很多

有人习惯把 MIPGap 设成0.05，想着“5%以内都能接受”。在小负荷场景下，目标成本里启停成本占比很高，5%的相对gap意味着可能差出好几台机组的启停方案。更麻烦的是，不同日期跑出来的结果，如果gap都比较松，目标成本的可比性就很差，你没法判断是模型改进了还是求解器随机碰到的解更好。解决方法是把 MIPGap 收紧到一个运行稳定的水平，比如0.001；如果某天大系统确实跑不满，再结合 TimeLimit 和 MIPFocus=2 去证明最优性。判断标准不是gap数字好不好看，而是连续跑几个案例，启停方案不出现大的抖动。

5.5 现象：终态把所有机组全停，凌晨又一片黑启动

单日优化很容易在最后一个时段把所有机组全停，因为模型看不到第25小时的需求，它觉得“未来不再有负荷”。放到连续滚动调度里，这个终态会让次日清晨的启动任务非常繁重，可能违反最小连续停机时间。解决方法是做滚动调度而不是孤立单日优化：把上一轮的第24时段状态作为下一轮 initial_status 传下去，并且在目标函数后面对最后若干时段加“最小在线机组数”约束，或者引入惩罚项让最后时段仍然保持合理组合。这个坑在纯学术算例里几乎没人提，但真实系统里非常常见。

6. 验证方法：一周回溯测试与两个可落地的扩展方向

6.1 一周回溯脚本：用历史负荷批量算指标

模型写出来不能只跑一天就算成功。我会拿至少一周的历史负荷数据做回溯测试，逐日跑24时段模型，并且把前一天的终态作为次日的初始状态。这个回测能暴露出很多单日算例看不见的问题：启停是否频繁、备用是否全天满足、求解时间是否稳定。

def build_uc_model(demand_day, reserve_day, init_status): # 复用第3章的建模代码，把initial_status替换为入参init_status pass def run_backtest(load_curves): last_status = [0] * 6 metrics = [] for date in sorted(load_curves): demand_day = load_curves[date] model = build_uc_model(demand_day, [round(d * 0.1) for d in demand_day], last_status) # 调求解器，MIPGap设0.001，TimeLimit设120 solve(model) last_status = [int(pyo.value(model.u[i, 23]) > 0.5) for i in range(6)] startup_count = sum( 1 for i in range(6) for t in range(24) if pyo.value(model.v[i, t]) > 0.5 ) metrics.append({ "date": date, "cost": pyo.value(model.obj), "startup_count": startup_count, "solution_time": result.solver.wallclock_time }) return metrics

回测之后我主要看四项指标：总成本、启停次数、最小启停时间违规次数、求解耗时。总成本要和同一天的其他方法对比才有意义；启停次数用来观察模型是不是在频繁启停；最小启停时间违规次数理论上必须是0，一旦出现先查v/w的边界；求解耗时则决定这套模型能不能塞进日内滚动的时间窗。

6.2 两个可落地的扩展方向：滚动窗口与SCED联动

第一个扩展是滚动窗口。把24小时模型改成“每4小时滚动一次、每次优化未来24小时”，能显著缓解终态效应，也让检修计划和负荷预测更新及时进入模型。滚动窗口的最大改动是把初始状态和终端约束都参数化，终端约束不是随便加的，而是把最后8小时的最小在线机组数抬上去，给下一轮调度留余地。

第二个扩展是和SCED联动。机组组合负责u，安全约束经济调度负责p，两步走是工程上最常见的解耦方式。先让MILP把所有时段的最小启停方案算出来，锁定u，然后交给SCED做时段内经济调度，这样既拿到混合整数线性规划在组合层面的全局性，又能在SCED里补网络约束、网损和备用分区。这个流程也方便你对照人工排班结果，看到底是组合方案的差异在影响成本，还是功率分配里出了问题。

这些年我做机组组合模型有个习惯：任何新版本上线前，先用三到五天的历史数据回测，成本比上一版低、启停次数比上一版稳定，才敢拿去给调度值班用。求解器参数换一个环境就重测一轮，不能拿昨天的参数赌明天的运行结果。希望帮到你。

本文还有配套的精品资源，点击获取