
简介面向电力系统调度与优化领域的学习者和研究人员该资源围绕机组组合问题基于混合整数线性规划方法利用MATLAB、YALMIP与Cplex构建求解框架可在满足负荷需求、启停成本、功率上下限等约束条件下制定经济可行的发电计划。压缩包共7个文件包括.m格式的模型代码、3个Excel运算结果表、2个Visio最优出力图表和1个Word基本要求说明文档整体大小仅267KB结构清晰。资源覆盖热备用系数0.05与0.2两种工况下的机组组合求解结果以及对应的机组各时段最优出力曲线同时附有详细说明便于对照理解MILP建模思路、约束设定、求解器调用与结果分析。已有3051人学习浏览适合需要快速掌握机组组合优化方法、复现典型算例或进行相关课程设计的读者参考。1. 机组组合问题为什么需要混合整数线性规划电力系统调度部门每天都要回答同一个问题明天 96 个时段里哪些机组开、哪些机组停、每台带多少负荷。机组组合Unit Commitment, UC就是做这个决策的它和纯经济调度最大的区别在于多了整数变量——机组的开停状态只有 0 和 1不能开半台也不能让机组处于微微启动的中间态。混合整数线性规划MILP正好是处理连续功率 离散状态这类混合决策的数学框架。很多人第一反应是用线性规划不就行了真跑起来才会发现把 0/1 状态变量直接丢进 LP 模型求解器给出的往往是机组半开半停这类物理上不可能的解。MILP 通过分支定界法在整数解空间里搜索虽然计算量比 LP 大得多但能严格保证找到全局最优解或给出最优性间隙gap。对电力系统这种动辄上百台机组、数万个约束的优化问题MILP 在工程上已经是标配——商业求解器如 Gurobi、CPLEX 能处理数万整数变量的问题开源求解器 CBC、SCIP 在中小规模场景下也完全够用。本文就从数学模型、代码实现、求解参数调优到实际排产落地完整走一遍机组组合优化的实现路径。2. 机组组合的 MILP 数学模型目标函数与约束体系2.1 目标函数燃料成本、启动成本与空载成本的分量组合机组组合优化要最小化的总成本通常包含三部分燃料成本、启动成本、空载成本。燃料成本一般用机组出力的二次函数近似但在 MILP 里需要线性化处理。常见做法是分段线性化——把出力范围切成若干段每段用一条直线近似段数越多精度越高但变量和约束也会线性增加。启动成本更特殊它取决于机组停机了多久停机时间短锅炉还是热的启动成本低停机超过一定小时数就是冷启动成本高。工程上通常简化为两段或三段阶梯函数用二进制变量表示热启动和冷启动两种状态。模型里启动成本还要区分启动动作和运行状态——只有机组从 0 变 1 的那一刻才计一次启动成本而不是每个时段都计。一个典型的目标函数写法如下以最小化总成本为例# 目标函数总成本 燃料成本 启动成本 空载成本 # 燃料成本已做分段线性化每段斜率为 c[k]段出力为 p[k] model.setObjective( quicksum( quicksum(c[i][k] * p[i][k][t] for k in range(K)) # 燃料成本 startup_cost[i] * u_start[i][t] # 启动成本 no_load_cost[i] * u[i][t] # 空载成本 for i in range(N) for t in range(T) ), GRB.MINIMIZE )这段代码里u[i][t]是机组 i 在时段 t 的开停状态0/1u_start[i][t]是启动动作变量只有机组从停机变开机的那一个时段为 1。p[i][k][t]是机组 i 在时段 t 第 k 段的出力增量。三部分成本分开建模的好处是可以分别设置权重——比如现货市场中启停频繁启动成本的准确性对报价影响很大。提示启动成本在 MILP 里必须用启动动作变量而不是运行状态变量来建模否则每个运行时段都会被重复计一次启动成本结果会严重偏大。2.2 约束体系系统平衡、备用容量与机组爬坡约束系统约束是机组组合的核心骨架。最基本的两个是功率平衡约束和备用容量约束。功率平衡要求每个时段所有机组出力之和等于该时段的系统负荷备用容量要求开机机组的最大可调出力之和减去当前出力要大于等于系统的旋转备用需求。这两个约束都作用于全系统层面把单台机组的决策耦合在一起是问题难解的根源之一。单个机组层面的约束包括出力上下限、最小运行时间、最小停机时间、爬坡速率约束。其中最小运行/停机时间是典型的整数约束——机组一旦开机就必须至少连续运行若干小时一旦停机就必须至少停够若干小时。这类约束用状态变量 时间窗求和的形式建模# 最小运行时间约束若机组在 t 时段启动则 t ~ tmin_up-1 时段必须为运行状态 for i in range(N): for t in range(T - min_up[i] 1): model.addConstr( quicksum(u[i][tau] for tau in range(t, t min_up[i])) min_up[i] * u_start[i][t] ) # 最小停机时间约束若机组在 t 时段停机则 t ~ tmin_down-1 时段必须为停机状态 for t in range(T - min_down[i] 1): model.addConstr( quicksum(1 - u[i][tau] for tau in range(t, t min_down[i])) min_down[i] * u_stop[i][t], namefmin_down_{i}_{t} )爬坡约束则连接相邻时段的出力变化。这里有一个容易踩的坑爬坡约束不仅要约束同一台机组相邻时段的出力差还要考虑启动/停机时的出力跳变。机组启动的时段出力从 0 直接跳到最小技术出力这个跳变往往超过正常爬坡速率需要单独处理。实务上有两种做法一种是简化模型忽略启动/停机瞬间的爬坡限制另一种是引入启动出力轨迹变量把启动过程拆成几个时段逐步加出力。大规模系统里常用前者小系统精确模拟用后者。2.3 机组组合完整代码用 Python Gurobi 实现一个小型系统下面的代码实现了一个 3 台机组、24 个时段的机组组合模型包含功率平衡、备用、爬坡和最小启停时间约束。数据可以直接替换成实际机组参数运行import gurobipy as gp from gurobipy import GRB # 机组参数最大出力, 最小出力, 爬坡速率, 最小开/停机时间, 启/停成本, 空载成本 gen [ {pmax: 300, pmin: 100, ramp: 60, min_up: 4, min_down: 3, su: 8000, sd: 4000, no_load: 500, c: 0.5}, {pmax: 200, pmin: 50, ramp: 40, min_up: 3, min_down: 2, su: 5000, sd: 2500, no_load: 300, c: 0.7}, {pmax: 100, pmin: 20, ramp: 30, min_up: 2, min_down: 1, su: 3000, sd: 1500, no_load: 200, c: 1.0}, ] T 24 load [280 50 * (1 (t % 12) / 6) for t in range(T)] # 模拟负荷曲线 reserve [0.1 * l for l in load] # 10% 旋转备用 m gp.Model(unit_commitment) N len(gen) u m.addVars(N, T, vtypeGRB.BINARY, nameu) p m.addVars(N, T, vtypeGRB.CONTINUOUS, namep, lb0) u_start m.addVars(N, T, vtypeGRB.BINARY, namestart) u_stop m.addVars(N, T, vtypeGRB.BINARY, namestop) # 出力上下限与爬坡约束 for i in range(N): for t in range(T): m.addConstr(p[i, t] gen[i][pmax] * u[i, t]) m.addConstr(p[i, t] gen[i][pmin] * u[i, t]) for i in range(N): for t in range(1, T): m.addConstr(p[i, t] - p[i, t-1] gen[i][ramp]) m.addConstr(p[i, t-1] - p[i, t] gen[i][ramp]) # 系统功率平衡与备用 for t in range(T): m.addConstr(quicksum(p[i, t] for i in range(N)) load[t]) m.addConstr(quicksum(gen[i][pmax] * u[i, t] - p[i, t] for i in range(N)) reserve[t]) # 目标函数启动与停机成本 obj quicksum(gen[i][c] * p[i, t] gen[i][no_load] * u[i, t] for i in range(N) for t in range(T)) obj quicksum(gen[i][su] * u_start[i, t] gen[i][sd] * u_stop[i, t] for i in range(N) for t in range(T)) m.setObjective(obj, GRB.MINIMIZE) m.optimize()这段代码的运行逻辑是Gurobi 在后台执行分支定界branch and bound和割平面cutting plane算法u[i,t]被当作整数变量处理求解器通过松弛 LP 求解、分支、剪枝逐层逼近整数最优解。参数里ramp是机组每时段能增减的出力上限min_up/min_down是最小连续运行/停机时段数。注意代码中的爬坡约束没有处理启停瞬间的出力跳变这是为了保持模型简洁——实际工程中需要在启动/停机的前几个时段加特殊的出力轨迹约束。3. 求解器选型与混合整数线性规划算法参数调优3.1 商业求解器与开源求解器的实际差异机组组合问题规模差异极大一个省级系统可能有几百台机组、96 个时段整数变量轻松超过十万个一个园区微网可能只有 5~10 台机组、24 个时段几千个整数变量。规模不同求解器选型逻辑完全不同。商业求解器 Gurobi 和 CPLEX 是电力行业事实标准它们经过数十年优化对大规模 MILP 的分支策略、割平面生成、预处理都做得非常成熟。1000 台机组、96 时段的 UC 问题约 10 万二进制变量Gurobi 通常能在几分钟内给出 1% 以内的最优性间隙。开源求解器 CBC 和 SCIP 在中小规模100 台机组以内表现尚可但到大规模问题时求解时间可能差一个数量级。选型建议教学演示和中小规模验证用 CBC安装简单PuLP 直接调用生产系统或大规模调度用 Gurobi 或 CPLEX。另外 HiGHS 是个值得关注的新选择它对 LP 和 MILP 的求解速度在开源里属于第一梯队。3.2 四个必调的 MILP 求解参数MILP 求解器的默认参数覆盖通用场景但对机组组合这类特定问题调参往往能带来 2~5 倍的加速。以下四个参数是我在实际项目中几乎每次都要碰的参数Gurobi 写法作用与建议值MIP 间隙MIPGap0.01提前终止条件1% 对调度够用现货市场出清要求 0.01%~0.1%时间限制TimeLimit300防止求解时间失控到时间返回当前最优解和间隙并行线程Threads0默认用满所有核共享服务器上建议手动限制割平面强度Cuts2网络拓扑和机组特征相关-1到3都试一下其中 MIPGap 是工程上最有用的参数。机组组合问题的整数解非常密集最优解附近往往有成百上千个近似相等的可行解。把间隙从默认的 1e-4 放宽到 1e-2求解时间常常能缩短 80%而成本差异只有百分之零点几——对日前调度来说这点误差完全可以接受。注意MIPGap 是相对间隙定义是|best_bound - best_solution| / |best_solution|。设 0.01 意味着当找到的可行解与理论下界差距在 1% 以内时停止搜索。如果想要严格最优解设 0 或很小的值但要接受求解时间可能呈指数增长。3.3 计算时间失控时的工程处理手段机组组合问题本质是 NP-hard最坏情况下求解时间无上限。工程上必须准备达不到最优怎么办的预案。除了上面提到的 MIPGap 和时间限制还有两个有效手段。第一个是给出好的初始可行解warm start。用上一轮调度结果或启发式算法比如优先顺序法先算出一个可行开机方案作为 MIP 起点的 MIP start。这能让求解器一开始就有可行解之后分支定界可以围绕这个解快速收敛——更早地剪枝、更早地逼近全局最优。Gurobi 里用m.Start属性传入初始状态变量即可。第二个是分解算法。把全系统 UC 按机组分解成单机子问题用拉格朗日松弛或交替方向乘子法ADMM迭代求解每次子问题都是小规模 MILP非常快。代价是收敛到全局最优没有保证——工程上通常用分解法求出高质可行解再用它作为原 MILP 的 MIP start 来收敛到最优。这个先分解、后精修的组合方案是很多实际调度系统的底层逻辑。4. 电力系统机组组合的实际约束与场景落地4.1 网络安全约束把潮流计算嵌入优化模型前面模型只有系统级的功率平衡和备用约束没考虑输电网的传输能力限制。实际系统中机组开在哪里和负荷在哪里同样重要——A 区开机过剩、B 区负荷紧张即使全系统功率平衡也可能因为线路容量不够导致无法送电。这就是网络安全约束机组组合SCUC。常见做法是把直流潮流DC Power Flow方程嵌入 MILP。线性灵敏度因子PTDF矩阵把线路潮流表示为节点注入功率的线性函数这样潮流约束可以写成纯线性约束保留 MILP 结构# 线路潮流约束-limit[l] sum(PTDF[l][n] * (p_gen[n] - load_demand[n] for n in nodes)) limit[l] for l in range(L): m.addConstr( quicksum(PTDF[l][n] * (net_injection[n]) for n in range(N)) line_limit[l], namefline_upper_{l} ) m.addConstr( quicksum(PTDF[l][n] * (net_injection[n]) for n in range(N)) -line_limit[l], namefline_lower_{l} )PTDF 矩阵怎么算对电力系统潮流计算有基础的话可以从节点电纳矩阵 B 求逆得到PTDF Bf * Bbus^{-1}其中 Bf 是支路-节点关联矩阵乘以支路电纳。这个矩阵是常数只取决于网络拓扑和线路参数不随调度结果变化所以可以预先算好、作为常量输入优化模型。这就是线性化潮流的威力——把潮流计算从优化模型里剥离出来预先算成系数矩阵。实际工程中注意 PTDF 方法只适用于交流系统的直流近似忽略无功和网损。对输电网级机组组合几百个节点精度足够对配电网或高阻线路较多的区域可能需要引入网损修正项——把网损近似为出力的线性函数或迭代修正。4.2 新能源接入后的机组组合场景法与鲁棒优化光伏和风电的大规模接入让机组组合从确定性问题变成了不确定性问题。负荷预测还有一定精度风电光伏的出力预测误差可以到 20%~30%。传统确定性 UC 要求备用容量覆盖预测误差但新能源占比越高这个方式的成本越高——要为极端场景预留大量机组空转经济性很差。工程上两类主流方案随机规划场景法和鲁棒优化。随机规划的做法是生成若干新能源出力场景每个场景对应一套机组出力决策但开机决策在所有场景下统一——这叫非预期性约束non-anticipativity。模型变成# 随机UC场景 s 下的出力决策共享同一组开机状态 u[i,t] for s in range(S): for t in range(T): model.addConstr( quicksum(p[i, t, s] for i in range(N)) load[t, s], # 场景s的负荷 namefbalance_s{t}_t{t} ) # 开机状态 u[i,t] 不随场景变化这个模型的规模随场景数线性增长场景数超过 20 个时计算压力已经很大。工程上常用场景削减技术如快速前向选择从几千个原始场景里挑出最典型的 10~20 个保证代表性同时控制计算量。鲁棒优化则不问概率分布只定义一个不确定集合比如风电出力在预测值 ±15% 的区间内波动要求机组组合方案对这个集合内的所有可能出力都可行。求解上经常用到对偶变换把 max-min 问题转化为单层 MILP。代价是解偏保守——它保证的是最坏情况也可行实际运行中大部分场景远没这么恶劣。两条路线各有适用场景系统备用充裕、新能源渗透率还不高的用场景法渗透率高、安全裕度要求极其严格的用鲁棒优化。4.3 现货市场出清场景机组组合与出清价格的耦合国内电力现货市场建设中机组组合不仅是物理调度工具也是市场出清的核心环节。日前市场出清要做安全约束机组组合SCUC和安全约束经济调度SCED——先用 SCUC 算出开机组合再用 SCED 在固定开机的条件下计算每台机组的出清出力并基于节点的边际电价形成节点电价LMP。这个场景下目标函数不再是最小化总成本而是最小化申报价格下的购电费用。机组的申报曲线往往不是单段线性而是分段阶梯报价这正好天然适配 MILP 的分段线性约束。出清价格由对偶变量影子价格给出——功率平衡约束的对偶乘子就是系统能量价格线路阻塞约束的对偶乘子则贡献阻塞分量。在这个场景中MILP 求解的数值稳定性变得极其重要。价格计算依赖对偶变量的精确值任何数值误差都可能放大成电价偏差。需要关注求解器的数值精度参数Gurobi 的NumericFocus并在求解后检查约束违反度。这也解释了为什么现货市场出清对 MIPGap 的要求远高于传统调度——1% 的间隙对物理调度没问题但对市场出清可能意味着几百万元的价差。5. 用热启动与迭代细化把求解时间压缩一个数量级最后落到一个实用技巧大多数机组组合在连续多日运行时前后两天的结果高度相似——昨天的机组组合完全可以直接作为今天的初始解。利用这个特性做热启动warm start求解时间往往能显著下降。Gurobi 里设置热启动的方式是把整数变量的值作为 MIP start 传入。如果上一轮的调度结果是可行的MIP start 可以直接给求解器一个可行下界起点分支定界算法可以据此快速收紧间隙# 假设 prev_schedule 是上一轮求得的开机状态矩阵 for i in range(N): for t in range(T): u[i, t].Start prev_schedule[i][t].X # 从上一轮结果复制初值 m.Params.MIPGap 0.005 # 收紧间隙 m.Params.TimeLimit 120 # 设一个硬性时间上限 m.optimize()配合热启动的另一个手段是迭代细化。先用粗时间粒度求解——比如把 96 个时段聚合成 24 个时段每 4 个时段合并为一段用平均负荷快速度算出近似开机方案。然后把这个方案作为 96 时段精细模型的 MIP start再求解完整模型。粗模型能捕捉机组开停的大致格局精细模型负责修正小时间尺度的爬坡和负荷跟随。我实测过不少场景这个方法比直接求解 96 时段模型快 5~10 倍而且最终解质量几乎不受影响——因为开机格局是问题的主要解结构细化阶段只是在这个格局下找最优出力分配。验证热启动效果有个简单方法对比冷启动和热启动的运行日志看第一个可行解出现的时间H标志行和 MIPGap 收敛曲线的斜率。热启动通常能在前几秒内就输出一个高质量可行解而冷启动可能要经过数次分支才能碰到第一个可行解。如果热启动没有明显加速检查约束是否挂得干净、上一轮结果的机组启停是否有变化——如果负荷或机组检修状态大幅变化热启动的价值会减弱此时可以考虑只对未变化机组做热启动、变化大的机组留空。机组组合是理论和工程结合最紧密的优化问题之一模型完备性和求解效率同样重要。本文还有配套的精品资源点击获取