简介面向电力系统调度与最优化方向学习者、研究者的机组组合优化资源包完整演示基于混合整数线性规划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 约束条件逐条写从负荷平衡到最小启停时间一个能跑的机组组合模型至少要有下面这张约束表约束表达式作用负荷平衡Σ_i p_it D_t任意时段总出力等于负荷旋转备用Σ_i Pmax_i · u_it ≥ D_t R_t在线机组可调上限覆盖负荷加备用出力上下限Pmin_i · u_it ≤ p_it ≤ Pmax_i · u_it开机时出力受限停机时出力强制为0爬坡约束p_it - p_i,t-1 ≤ RU_i M·(2-u_it-u_i,t-1)限制相邻时段出力增量最小连续运行Σ_{kt-UT_i1}^{t} v_ik ≤ u_it启动后必须持续运行若干小时最小连续停机Σ_{kt-DT_i1}^{t} w_ik ≤ 1-u_it停机后必须保持停机若干小时启停事件逻辑v_it - w_it u_it-u_i,t-1v_itw_it ≤ 1启动/停机事件与状态变化一致爬坡约束里的 M 是松弛项我用 (2-u_it-u_i,t-1)·Pmax_i 来替代大M。当机组在相邻时段都在线时2-u_it-u_i,t-10爬坡限制严格生效只要有一侧不在线后面那项就自动放大到机组容量相当于放过启动/停机瞬间的出力跳变。最小启停时间的约束写法很关键启动事件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(initializerange(6)) m.T pyo.Set(initializerange(24)) # u为启停状态p为出力v/w为启动/停机事件 m.u pyo.Var(m.G, m.T, domainpyo.Binary) m.p pyo.Var(m.G, m.T, domainpyo.NonNegativeReals) m.v pyo.Var(m.G, m.T, domainpyo.Binary) m.w pyo.Var(m.G, m.T, domainpyo.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(ruleobjective_rule, sensepyo.minimize)目标函数里no_load乘umc乘pstartup乘v。注意启动成本只跟启动事件v挂钩不是跟“状态变化”的直接差值挂钩这样每一次冷态启动只计一次费用。如果某台机组初始就在线t0不需要启动事件这个边界会在后面约束里处理。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, ruleload_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, rulereserve_rule) # 出力上下限u0时强迫p0 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, rulep_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, rulep_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, rulestartup_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, ruleshutdown_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, rulestartup_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, ruleramp_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, ruleramp_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, rulemin_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, rulemin_down_rule)这段代码里最值得看的是 t0 分支。如果忽略初始状态模型会默认“第一小时前全部机组停机”这会让启动事件v在第一个时段失去约束机组可以白启动而不计启动成本。这里我用 initial_status 把真实初始状态带进 v[0] 和 w[0] 的边界保证边界不漂。爬坡约束用 (2-u_it-u_i,t-1)*pmax 而不是写死的大M好处是你不用去调M的值pmax本身就能起到“不在线就不约束”的作用缺点是当相邻时段有一侧离线时约束右端会放大这是工程近似不是严格建模。如果你要处理启动过程中的分段出力轨迹需要另外引入启动轨迹变量但6机模型用这个松弛就够了。3.4 求解与结果查看先确认能跑出可行解模型建完剩下的就是用求解器跑。我偏好Gurobi没有授权就用开源的CBCPyomo的接口是统一的。solver pyo.SolverFactory(gurobi) solver.options[MIPGap] 0.001 solver.options[TimeLimit] 120 solver.options[Threads] 4 res solver.solve(m, teeTrue) 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(fG{i1}:{pyo.value(m.p[i, t]):6.1f}MW for i in on) print(ft{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里基本一致只是写入方式略有差异参数典型值什么时候动MIPGap0.00124时段小系统做精度验证大系统可以先放到0.01MIPGapAbs视成本量级而定当相对gap失真时改用绝对gap控制TimeLimit60300秒日内滚动调度要留出足够的通信和校验余量Threads48单机多核并行跑分支MIPFocus1/2/3想快速拿可行解设1想证明最优性设2平衡设3Gurobi里设置方式是 solver.options[MIPGap] 0.001CBC里一般是 solver.options[mipgap] 0.01具体参数名以你装的版本帮助为准。不要一上来就设 TimeLimit3600 和 MIPGap0.0001整数规划的证明最优在机组组合这种规模下可能跑出几小时结果就是调度方案出得太慢反而没人敢用。我的习惯是先放宽到 TimeLimit120、MIPGap0.01拿一个可行解看约束有没有问题确认模型正确后再收紧gap。MIPFocus 也是有用的优化项如果你有大量机组需要快速给出一个次优可行方案MIPFocus1 会让求解器更早剪枝虽然牺牲一点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同时把t0的边界处理好。只要变量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 和 MIPFocus2 去证明最优性。判断标准不是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.001TimeLimit设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里补网络约束、网损和备用分区。这个流程也方便你对照人工排班结果看到底是组合方案的差异在影响成本还是功率分配里出了问题。这些年我做机组组合模型有个习惯任何新版本上线前先用三到五天的历史数据回测成本比上一版低、启停次数比上一版稳定才敢拿去给调度值班用。求解器参数换一个环境就重测一轮不能拿昨天的参数赌明天的运行结果。希望帮到你。本文还有配套的精品资源点击获取
