简介这是一份面向电力系统调度员、电气工程学生及优化算法研究者的机组组合优化资源聚焦混合整数线性规划建模与求解。模型以整数变量表示机组启停、连续变量表示出力以最小化发电成本为目标并纳入功率上下限、热备用、负荷平衡等约束借助YALMIP建模并调用CPLEX求解器。整个压缩包共7个文件体积仅267KB包含可运行的MATLAB优化脚本、题目说明文档、两种热备用情景0.05与0.2下的求解结果表以及机组最优出力可视化图表覆盖从模型定义、数据输入、求解器配置到结果解读的完整流程。已有3051人学习/下载适合希望快速上手混合整数规划在机组组合中实际落地的读者。通过学习代码与结果可掌握目标函数与约束条件设置、求解器调用方法并对比不同热备用水平对启停方案和总成本的影响便于迁移到自身电力系统优化场景中。1. 为什么机组组合必须交给混合整数线性规划调度员的“开机决策”不是连续问题电力系统调度员每天面对的最棘手问题不是“某台机发多少电”而是“明天哪些机组要开机、哪些要停机”。这个决策一旦做错后面经济调度算得再准也得为错误买单。机组组合优化Unit Commitment, UC就是解决这个问题的建模框架而混合整数线性规划MILP是目前把机组组合问题落到求解器里最主流、最可靠的数学形式。它把“开还是不开”的0-1决策和“发多少”的连续出力统一进同一个优化模型既能描述机组最小启停时间这类离散逻辑又能用成熟的分支定界算法保证收敛到全局最优解。对刚入门的工程师这个标题意味着你要掌握三件事怎么把调度规则翻译成线性约束怎么选择求解器并调参以及怎么在真实数据上把结果验证到调度员敢用。对有经验的从业者这篇笔记的价值在边界条件哪些约束可以线性化、哪些必须妥协、MIP gap设到多少才不会让调度员翻车。数据、代码、参数、踩坑我按一次完整落地的顺序写下来。2. 把UC问题写成MILP目标函数、关键约束与数学化表达2.1 目标函数怎么拆发电成本与启停成本的不同量纲机组组合的目标函数通常是最小化系统总运行成本它不是一个“总费用”那么简单而是两笔不同性质的费用叠加。第一笔是燃料成本只和机组出力水平有关一般用二次曲线近似但在MILP里二次函数会让求解变慢工程上普遍做法是分段线性化。第二笔是启停成本包括锅炉点火、汽轮机冲转、额外燃料消耗等它是固定成本只要机组状态从0变1就要付一次和出力多少无关。这两笔钱必须分开建模原因很直接如果不单独建模启停成本求解器可能会让机组频繁启停来“优化”燃料成本因为停机时段不发电就不产生燃料费用实际上启停一次的成本远高于省下的那点燃料费。所以0-1变量的引入不只是数学上的需要它直接对应着调度规则。启动成本和停机成本在有些模型里分开计费在有些模型里简化为只算启动成本我一般建议至少保留启动成本项停机成本视机组类型而定——燃煤机组停机成本不可忽略燃气机组往往可以略掉。目标函数的另一个细节是机组燃料成本曲线的分段处理。二次函数f(p)a·p²b·pc在MILP里只能通过分段线性近似实现每段引入一个连续变量和一个0-1标志变量。段数太少精度差段数太多求解时间暴涨。对大多数燃煤机组3到5段足够对燃气轮机2到3段通常合理。不要为了“精确”把每台机组都切到10段求解时间会教你做人。2.2 功率平衡、备用与爬坡三类核心约束的线性化写法UC模型的核心约束按功能分三类每一类都有固定的线性化套路。功率平衡约束是硬约束它要求每个时段所有机组出力之和等于系统负荷属于等式约束Σ pᵢ,ₜ Dₜ。这里要注意如果做的是多节点模型还得把网络潮流约束加进来但那是另一个量级的复杂度先按下不表。备用约束是保证系统安全的约束常见写法是正备用和负备用分开写。正备用约束要求开机机组的最大出力之和至少覆盖负荷加备用需求Σ (Pᵢ,max · uᵢ,ₜ) ≥ Dₜ Rₜ。这个约束的线性化陷阱在于备用容量只统计“开机”的机组停机机组哪怕最大出力再大也不能算数所以必须乘上0-1变量uᵢ,ₜ。负备用约束类似要求开机机组的最小出力之和不高于负荷减去下限备用这在实际运行中同样重要否则负荷低谷时段所有机组都压不下去系统会被迫弃风弃光。爬坡约束描述的是机组出力在相邻时段之间的变化限制它的线性化写法是三个不等式上限方向、下限方向和启停耦合。最常见的是只写两个不等式pᵢ,ₜ − pᵢ,ₜ₋₁ ≤ RUᵢ 和 pᵢ,ₜ₋₁ − pᵢ,ₜ ≤ RDᵢ。但这有个隐藏问题如果机组在时段t刚刚启动pᵢ,ₜ₋₁是0那么第一个约束要求pᵢ,ₜ ≤ RUᵢ这恰好是对的但如果机组在时段t停机pᵢ,ₜ是0第二个约束变成pᵢ,ₜ₋₁ ≤ RDᵢ这可能会导致不可行。所以严谨的做法是引入启停状态变量把爬坡约束写成带big-M的形式后面避坑章再细讲。2.3 最小启停时间最容易被低估的约束机组一旦启动不能立刻停机一旦停机也不能立刻启动。这是物理规律也是调度安全底线。最小启停时间约束是UC模型里最容易写错、也最容易导致不可行的一类约束。它的标准线性化写法是基于机组在过去若干时段的启停历史来限制当前决策。最小运行时间约束的写法是如果机组从停机状态切换到启动状态那么接下来至少T_on个时段必须保持开机。用0-1变量yᵢ,ₜ表示启动动作从0到1的转变约束为Σₖ₌ₜ to ᵗ⁺ᵀᵒⁿ⁻¹ uᵢ,ₖ ≥ T_on · yᵢ,ₜ。这个约束的含义是在启动后的连续T_on个时段内机组必须全程开机。最小停机时间约束对称Σₖ₌ₜ to ᵗ⁺ᵀᵒᶠᶠ⁻¹ (1 − uᵢ,ₖ) ≥ T_off · zᵢ,ₜ其中zᵢ,ₜ表示停机动作。这里有两个细节容易被忽略。第一如果机组在初始时刻已经运行了若干小时那么初始时段的最小运行时间约束要把历史计入否则求解器可能在第一天就让机组停机违反物理规律。第二启动变量和停机变量之间要加互斥约束和启停一致性约束uᵢ,ₜ − uᵢ,ₜ₋₁ yᵢ,ₜ − zᵢ,ₜ同时yᵢ,ₜ zᵢ,ₜ ≤ 1。漏掉这两个约束模型会出现“既启动又停机”的自相矛盾结果目标函数值看起来合理实际不可执行。3. 实现选型与求解器参数从Gurobi到COPT的实战选择3.1 求解器选型商业求解器与开源求解器的边界模型写好了接下来要选求解器。国内电力系统调度机构里Gurobi、CPLEX和COPT是三大主流商业求解器三者对MILP的处理能力在UC这类问题上差距不大差距主要在授权价格、服务响应和国产化要求。如果你的项目有信创要求或者调度系统在国网南网体系内COPT是常见选择如果是高校科研或市场化项目Gurobi的生态和文档更友好。开源求解器里CBC和SCIP可以用但对大规模UC问题求解速度通常差一个数量级以上我一般只在教学或验证小模型时用。3.2 参数怎么调MIPGap、TimeLimit与数值容差MILP求解器的参数调节是一门“玄学”但熟手手里有固定套路。第一个必调参数是MIPGap它控制着“次优解与最优解之间的相对差距阈值”。调度场景里MIPGap设到0.01%万分之一的相对gap已经是极致实际工程0.1%到0.5%完全够用。先把MIPGap设为0.005再结合TimeLimit兜底让求解器到限时就把当前最好可行解交出来这是UC问题的黄金组合。第二个必调参数是TimeLimit它决定了调度员能不能在要求的时间窗口内拿到结果。日前机组组合的典型时间窗口是10到30分钟实时滚动修正只有1到5分钟。我一般把TimeLimit设为窗口上限的60%到70%留出时间给后处理和数据检查。第三个容易被忽视的参数是数值容差FeasibilityTol、OptimalityTol默认值通常够用但如果模型里有量纲差异很大的系数比如成本单位用“元”而功率单位用“MW”求解器会报numerical trouble这时候不是调容差而是先做单位归一化。3.3 用PyomoGurobi搭建最小可运行模型最常见的实现方式是用Python的Pyomo或PuLP建模然后调用求解器。Pyomo的优势是建模语法接近数学表达式支持随机规划扩展社区资料也多适合UC这类复杂约束较多的模型。PuLP更轻量适合快速原型验证但表达min-up/min-down这类复杂约束时不如Pyomo直观。我倾向于用Pyomo Gurobi的组合原因是后续如果要加场景集、加网络约束Pyomo的扩展性更好不会推倒重来。4. 用PyomoGurobi搭出最小可运行模型代码与参数逐一拆解4.1 数据准备从CSV到可建模的格式UC模型需要两类输入数据机组参数表和负荷预测曲线。机组参数包括最大/最小出力、爬坡速率、最小启停时间、启动成本、燃料成本系数。负荷预测数据一般来自调度系统的超短期或短期预测模块按15分钟或1小时的粒度给出。下面是一个机组参数表的CSV结构示例注意我这里把燃料成本预先做了分段线性化每段用斜率和截距表示。实际工程项目里这些参数由设备厂商提供或从历史运行数据回归得到参数质量直接影响优化结果的可信度。# generator_data.csv # gen_id, p_max, p_min, ramp_up, ramp_dn, min_on, min_off, # start_cost, stop_cost, seg_points(出力分段点), seg_slopes(分段斜率) G1, 500, 150, 200, 200, 4, 2, 20000, 5000, [150, 250, 400, 500], [0.28, 0.30, 0.32] G2, 300, 80, 120, 120, 2, 1, 8000, 2000, [80, 150, 220, 300], [0.35, 0.38, 0.40] G3, 100, 30, 50, 50, 1, 1, 3000, 800, [30, 60, 100], [0.45, 0.48] G4, 200, 60, 80, 80, 2, 1, 5000, 1500, [60, 100, 150, 200], [0.40, 0.42, 0.44]这里的段点seg_points和斜率seg_slopes是把二次成本曲线做一阶近似得到的。每段内的成本函数是线性的成本 截距 斜率 × 出力。截距怎么算用 fitted quadratic 在段点处的值和斜率反推。我在实际项目中会用一段脚本自动从原始机组数据生成这个表避免手算截距出错。4.2 模型定义变量、约束、目标import pyomo.environ as pyo model pyo.ConcreteModel() TIMES list(range(24)) # 24小时粒度1小时 GENS [G1, G2, G3, G4] # 决策变量 model.u pyo.Var(GENS, TIMES, domainpyo.Binary, doc机组启停状态) model.p pyo.Var(GENS, TIMES, domainpyo.NonNegativeReals, bounds(0, 500), doc出力) model.y pyo.Var(GENS, TIMES, domainpyo.Binary, doc启动动作) model.z pyo.Var(GENS, TIMES, domainpyo.Binary, doc停机动作) # 功率平衡约束每个时段总出力 负荷 model.balance pyo.Constraint(TIMES, rulelambda m, t: sum(m.p[g,t] for g in GENS) load[t])上面这段代码先定义了核心决策变量三组变量分别表示状态、出力和启停动作。注意p的bounds我粗设为0到500这是一个临时值实际应该从机组参数表读入下面代码会覆盖。功率平衡约束是等式约束负荷load是一个按时间索引的列表或字典。接下来是备用约束和启停一致性约束# 正备用约束开机机组的最大出力之和 负荷 备用 model.reserve pyo.Constraint(TIMES, rulelambda m, t: sum(p_max[g] * m.u[g,t] for g in GENS) load[t] reserve_req[t]) # 启停一致性约束状态变化 启动动作 - 停机动作 model.status_track pyo.Constraint(GENS, TIMES, rulelambda m, g, t: m.u[g,t] - (m.u[g,t-1] if t 0 else init_status[g]) m.y[g,t] - m.z[g,t]) # 启动/停机动作互斥 model.exclusive pyo.Constraint(GENS, TIMES, rulelambda m, g, t: m.y[g,t] m.z[g,t] 1)备用约束的本质是“预留空间”而不是“实际出力”所以用最大出力乘状态变量而不是用出力变量。启停一致性约束是模型里最容易写错的一条特别是t0时段需要处理初始状态这里用init_status字典保存机组在优化周期开始时的状态0或1。互斥约束保证一个时段内不会同时出现启动和停机两个动作。爬坡约束按前面说的带big-M方式写# 爬坡约束用big-M处理启停边界 M 600 # 大于任何机组最大出力的常数 model.ramp_up pyo.Constraint(GENS, TIMES, rulelambda m, g, t: m.p[g,t] - m.p[g,t-1] ramp_up[g] * m.u[g,t-1] M * (1 - m.u[g,t-1])) model.ramp_dn pyo.Constraint(GENS, TIMES, rulelambda m, g, t: m.p[g,t-1] - m.p[g,t] ramp_dn[g] * m.u[g,t] M * (1 - m.u[g,t]))这段爬坡约束是避坑的关键。当机组在t-1时段处于停机状态时m.u[g,t-1]0第一个约束右边的第一项消失但M×(1-0)M实际上把上升爬坡限制放宽到了无效状态——这是合理的因为一个停机机组启动后的第一段出力本来就允许从0跳到最小出力以上。当机组持续开机时u1约束退化为标准的相邻时段爬坡限制。这套写法比朴素写法更鲁棒代价是多了一个需要调整的M值M要足够大但不能太离谱否则数值困难。4.3 求解与结果解析solver pyo.SolverFactory(gurobi) solver.options[MIPGap] 0.005 solver.options[TimeLimit] 60 solver.options[LogToConsole] 1 results solver.solve(model, teeTrue) # 结果输出每台机组每时段启停状态和出力 for g in GENS: for t in TIMES: if pyo.value(model.u[g,t]) 0.5: print(f小时{t1}: {g} 开机, 出力 {pyo.value(model.p[g,t]):.1f} MW)MIPGap设置0.005表示求解器在找到相对gap小于0.5%的解后就停止这个阈值在日前机组组合场景下足够让调度员接受也能显著缩短求解时间。TimeLimit设为60秒是兜底策略实际生产环境如果周期是96点15分钟粒度TimeLimit可能要放大到300秒以上。运行结束后pyo.value抽取出数值注意二元变量在求解器返回时可能是0.9999或0.0001这样的浮点数判断状态要用阈值0.5而不是直接比较等于1。5. 机组组合建模避坑指南五条血泪经验5.1 现象求解器报“数值困难”或警告“unscaled”求解器日志里出现numerical trouble、dual infeasible这类警告或者同一个模型换个数据规模结果就剧烈变化先别怀疑算法大概率是你的模型数值尺度出了问题。我把成本单位设为“元”、出力单位设为“MW”、爬坡速率单位设为“MW/分钟”结果一个约束里出现了0.0001和5000这样跨七个数量级的系数求解器矩阵条件数爆炸。原因目标函数里燃料成本几十万备用约束里功率又是几百这些系数直接进约束矩阵导致线性规划子问题病态。解决把出力统一为标幺值或者以“百MW”“十MW”为单位成本以“万元”为单位。另外一个更省事的做法是在建模之前对原始数据做一遍标准化让所有参数落在10⁻²到10²这个量级区间。做完这一步95%的数值困难问题消失。5.2 现象MIPGap已经很小但机组启停方案前后两天差异巨大模型跑出来的总成本差距不到0.1%但机组的启停方案却完全不同——今天让G1开机明天让G2开机第三天又回到G1。调度员看着结果一脸茫然因为对这些方案来说总成本差异极小求解器被MIPGap阈值“放过”了。原因多台机组参数相近时它们在目标函数上是近似对称的对称性导致大量等价解分支定界的过程会在这堆等价解里乱撞先找到哪个是哪个。解决加对称性破缺约束或优先规则。常见做法是给机组按效率排序增加一组“如果G1和G2参数相同则倾向先开G1”的约束或者在目标函数里加一个极小的人工偏好系数。还有一种做法是给机组的热启动成本加微小差异打破对称性。这个操作不会影响最优解的有效性只是让求解器有稳定方向。5.3 现象解出来机组频繁启停一天内开停三次目标函数里明明有启动成本但优化结果还是让机组频繁启停。检查发现启动成本确实被计入了但启动成本的数值远小于机组空转或低出力时的燃料成本差导致求解器认为“停机更划算”。原因启动成本参数偏小或者忽略了停机成本导致频繁启停的总成本反而低于保持开机状态。解决核对机组参数表的启动成本单位是否正确在模型里同时计入停机成本对燃煤机组要特别注意它们的启停成本往往占比很高。另外可以考虑在目标函数中引入最小运行时间约束的拉格朗日惩罚项但根治方法还是把成本参数校准到设备厂商的实测数据。5.4 现象模型不可行日志提示是爬坡约束冲突一个含燃煤和燃气机组的系统在早高峰爬坡时段出现infeasible。检查约束发现燃气机组爬坡快但容量小燃煤机组容量大但爬坡不够某个时段的总爬坡能力不满足负荷变化需求。原因这可能是模型写错的不可行——爬坡约束没有和启停状态解耦出现了停机机组还在“爬坡”的荒谬约束也可能是数据真的不可行——系统爬坡能力不足但这种情况在实际调度中较少因为备用配置会兜底。解决先把刚才介绍过的big-M爬坡写法替换掉朴素写法排除模型本身的问题。如果仍然不可行把不可行时段对应的约束单独拉出来检查——用Gurobi的computeIIS功能找不可行子集这个功能在诊断阶段非常有用。确认真是爬坡能力不足的时候需要调整的是备用约束或增加快速响应机组而不是硬调爬坡参数。5.5 现象求解时间随着机组数量增加呈指数级爆炸从10台机组扩到50台求解时间从几十秒暴涨到几个小时。这是MILP最核心的痛点分支定界的搜索树节点数量在最坏情况下是2的指数级增长。原因模型里二元变量数量大约等于机组数×时段数50台机组×96时段就是4800个二元变量每个分支节点都要解一个大规模LP松弛问题。解决不是只有换更大机器一条路。先用MIPGap TimeLimit组合策略接受一个次优解再启用求解器的启发式初始解功能给Gurobi的MIPStart传入一个热启动方案从调度员手工方案或上一轮优化结果出发能大幅削减搜索空间。最后可以尝试用“滚动优化”替代一次性全时段优化把24或96时段切成长度为8到12时段的子问题重叠边界做缓冲这是工程上应对大规模UC的标准策略。6. 从“能跑”到“敢用”解的可信度验证与落地完善6.1 对可行解做重放验证优化结果看起来完美但调度员要求的是“可执行”。我的习惯是写一个独立的约束校验器把求解结果逐条回放到原始约束里验证功率平衡是否满足、爬坡是否超限、最小启停时间是否从第一天开始就合法。这条工序听起来多余但MILP求解器在某些数值边缘情况下会返回近似整数解比如0.99被当作1或者约束松弛后的解通过了求解器的可行性检查但实际违反物理规律。重放验证不需要重新求解只是逐条不等式检查毫秒级完成但能拦住大量低级错误。for t in TIMES: total_power sum(dispatch[g][t] for g in GENS) assert abs(total_power - load[t]) 1e-2, f时段{t}功率不平衡这段校验代码的原理极其简单用断言检查每个时段的总出力和负荷差值是否小于阈值。真正的工程校验器会比这个复杂得多要检查爬坡、备用、启停时间等所有约束但核心思想一致把优化结果当黑盒做回归测试而不是相信求解器“可行”标记。6.2 MIPGap到底设多大合适以及给调度留裕度MIPGap值反映的是“你愿意为更快求解牺牲多少最优性”。0.5%的gap意味着最多比最优解贵0.5%对日前机组组合这种日成本百万级别的体量绝对误差可能在几千到几万元但换来的是求解时间从小时级降到分钟级。我一般在日常运行场景用0.5%到1%在方式计算或规划研究中才收紧到0.01%到0.05%。另一个技巧是先用宽松gap快速得到一个可行解缓存的MIPStart传给下一轮收紧gap的求解。6.3 向两阶段模型演进日前机组组合和实时经济调度衔接机组组合只是电网调度链条的第一环。日内实时阶段负荷预测更新后机组启停状态已经锁定只需要做经济调度来微调出力。这个两阶段结构意味着UC模型的前瞻性非常重要如果你在日前阶段做出一个总成本最优但爬坡余量极小的方案实时阶段一旦负荷偏离预测值一两百MW调度员就得紧急启停机组来救场。所以有些现场系统会在备用水准上设保守系数或者在目标函数里加一个“不平衡惩罚项”来间接预留调节容量。这些方法不应偏离你的精确最优解初衷但它反映了工程落地的真相优化结果要给运行留余量而不是把系统压到极限。我最早做UC项目时只盯着目标函数值有没有优化到位被调度员问一句“你这个方案实时阶段怎么调整”就哑火了。后来养成的习惯是每次做完日前优化都把预测负荷上下偏移2%、5%各跑一遍灵敏度分析确认启停方案在这些偏移下仍然基本可行。这个步骤花不了多少时间但能让你交付的结果真正变成调度台可用的东西希望帮到你。本文还有配套的精品资源点击获取