★ FEATURED ARTICLE 简介这份资源面向电力系统优化方向的研究生、科研人员与电力工程师聚焦安全约束机组组合SCUC与经济调度问题提供一套可运行的MATLAB/CPLEX求解方案。压缩包共11个文件约301KB包含6个xls数据表、2个m脚本、2个vsdx流程图和1份docx说明文档分别用于存放节点导纳矩阵、机组出力与求解结果、核心优化程序及建模思路整理。内容围绕IEEE-30节点系统展开涵盖机组启停成本与最小运行时间、负荷平衡、线路潮流与备用约束等关键环节并给出热备用0.05与0.2两种场景下的机组组合结果对比。已有1960人学习适合希望理解SCUC建模流程、掌握混合整数线性规划求解思路并对照结果验证算法的读者参考。1. 从一份压缩包说起SCUC 机组组合到底在算什么如果你手里也躺着一个名为「考虑SCUC的电力系统机组组合(含经济调度)优化」的压缩包点开之前最好先搞清楚它到底在解决什么问题。SCUC 全称 Security-Constrained Unit Commitment安全约束机组组合本质是一道混合整数规划在满足负荷平衡、备用容量、机组爬坡、最小启停时间以及网络安全约束的前提下决定未来若干时段内每台机组「开还是关」「开多少出力」让总运行成本最低。它和纯经济调度ED最大的区别在于ED 默认机组组合状态已知只优化出力分配而 SCUC 要先把 0/1 的启停决策定下来再在给定组合下做经济调度所以它是一个「先离散后连续」的两层结构。这类模型最典型的落地场景是省级或区域电网的日前发电计划编制时间尺度通常覆盖 24 个时段步长 1 小时也有 15 分钟粒度的版本。适合谁看做电力系统优化调度的研究生、电网调度自动化方向的工程师、以及需要把 SCUC 嵌进微电网能量管理系统或虚拟电厂聚合平台的开发者。压缩包里大概率是模型文件、算例数据和求解脚本的组合但无论里面是什么你都得先能自己把这条链路跑通否则拿到别人的包也只是个黑匣子。2. 把 SCUC 拆成可求解的数学模型从目标函数到约束矩阵2.1 目标函数与决策变量的物理含义SCUC 的目标函数由三块构成机组运行成本、启动成本、停机成本。运行成本通常是出力的二次函数工程上为了用 MILP 求解器一般做分段线性化。决策变量分两类连续变量 P[i,t] 表示机组 i 在时段 t 的出力二进制变量 u[i,t]、v[i,t]、w[i,t] 分别表示运行、启动、停机状态。启动成本又分热启动和冷启动取决于机组已经停机了多久这需要引入停机时间计数器或者用启动类型变量来区分。我一般会把目标函数写成min sum_i sum_t [ C_i(P[i,t]) SU_i * v[i,t] SD_i * w[i,t] ]其中 C_i 是分段线性成本SU_i 和 SD_i 是单次启动和停机成本。注意停机成本在很多教材里被忽略但在实际算例中如果停机不花钱优化器会倾向于频繁启停来套利导致结果不可信。2.2 约束条件逐条拆解与矩阵组装约束是 SCUC 最容易翻车的地方漏一条约束结果就可能完全不可行或者成本低得离谱。核心约束包括功率平衡sum_i P[i,t] D[t]D[t] 是时段 t 的系统负荷。旋转备用sum_i u[i,t] * Pmax[i] D[t] R[t]R[t] 是备用需求。出力上下限u[i,t] * Pmin[i] P[i,t] u[i,t] * Pmax[i]。爬坡约束P[i,t] - P[i,t-1] RU[i] * u[i,t-1] SU[i] * v[i,t]下降同理。最小启停时间用 u 和 v/w 的线性组合表达通常写成 sum_{kt}^{tT_on-1} u[i,k] T_on[i] * v[i,t]。网络安全约束直流潮流下用 PTDF 矩阵把线路潮流表示为注入功率的线性函数再限制在线路容量内。把这些约束写成矩阵形式就是 A_eq * x b_eq 和 A_ineq * x b_ineq。下面这段 Python 用 PuLP 搭一个 3 机 4 时段的迷你 SCUC方便你验证约束逻辑import pulp # 3台机组4个时段 n_gen, n_t 3, 4 Pmax [100, 80, 50] Pmin [20, 15, 10] cost_a [0.02, 0.03, 0.04] # 二次项系数 cost_b [30, 28, 35] # 一次项系数 D [120, 150, 140, 130] # 负荷 R [20, 20, 20, 20] # 备用 RU [40, 30, 20] # 上爬坡 RD [40, 30, 20] # 下爬坡 Ton [2, 2, 1] # 最小开机时间 Toff [2, 2, 1] # 最小停机时间 prob pulp.LpProblem(SCUC_mini, pulp.LpMinimize) u pulp.LpVariable.dicts(u, (range(n_gen), range(n_t)), catBinary) v pulp.LpVariable.dicts(v, (range(n_gen), range(n_t)), catBinary) w pulp.LpVariable.dicts(w, (range(n_gen), range(n_t)), catBinary) P pulp.LpVariable.dicts(P, (range(n_gen), range(n_t)), lowBound0) # 目标线性化成本用一次项近似启动成本固定 prob pulp.lpSum(cost_b[i] * P[i][t] 50 * v[i][t] for i in range(n_gen) for t in range(n_t)) for t in range(n_t): prob pulp.lpSum(P[i][t] for i in range(n_gen)) D[t] prob pulp.lpSum(u[i][t] * Pmax[i] for i in range(n_gen)) D[t] R[t] for i in range(n_gen): prob P[i][t] Pmin[i] * u[i][t] prob P[i][t] Pmax[i] * u[i][t] prob u[i][t] - u[i][t-1] v[i][t] - w[i][t] if t 0 else u[i][t] v[i][t] # 最小开机时间约束简化版仅对t0之后 for i in range(n_gen): for t in range(n_t): if t Ton[i] n_t: prob pulp.lpSum(u[i][k] for k in range(t, t Ton[i])) Ton[i] * v[i][t] prob.solve(pulp.PULP_CBC_CMD(msg0)) print(Status:, pulp.LpStatus[prob.status]) for i in range(n_gen): print([u[i][t].varValue for t in range(n_t)], [P[i][t].varValue for t in range(n_t)])这段代码里u[i][t] - u[i][t-1] v[i][t] - w[i][t]是启停逻辑的核心它保证状态切换时启动和停机变量不会同时为 1。最小开机时间约束用滑动窗口求和表达Ton[i]是机组 i 的最小开机时段数。求解器用 CBC 是因为它免费且对中小规模 MILP 够用但如果你要跑 100 台机组 24 时段建议换 Gurobi 或 HiGHSCBC 的分支定界在规模上去后会明显变慢。2.3 经济调度在 SCUC 中的嵌套位置很多人把 SCUC 和 ED 混在一起写结果模型又大又慢。正确的做法是分层SCUC 决定 u/v/wED 在 u 固定的情况下优化 P。但在实际求解时两者往往在同一个 MILP 里同时求解因为 ED 的约束功率平衡、爬坡本身就是 SCUC 约束的一部分。所谓「含经济调度」通常是指模型在确定机组组合后再单独跑一次纯 LP 的 ED 来校验成本或者用 ED 的结果做安全校核。我一般会在 SCUC 求解完成后把 u 固定重新用 LP 求解 P对比两次的 P 是否一致如果不一致说明 SCUC 里的线性化成本或爬坡约束有遗漏。3. 用 Python 搭一套可复现的 SCUC 求解流程3.1 数据准备负荷曲线、机组参数与网络参数数据是 SCUC 的燃料格式不对后面全白搭。我习惯用三个 CSVload.csv存时段和负荷gen.csv存机组参数network.csv存线路和 PTDF。机组参数至少要有 Pmax、Pmin、二次成本系数、启动成本、最小启停时间、初始状态。网络参数如果做直流潮流需要支路电抗和容量以及节点-机组映射。import pandas as pd gen pd.read_csv(gen.csv) # 列gen_id, Pmax, Pmin, a, b, c, SU, SD, Ton, Toff, init_state, init_hours load pd.read_csv(load.csv) # 列hour, demand, reserve print(gen.head()) print(load.head())init_state表示机组在第一个时段之前是开还是关init_hours表示已经持续了多少个时段。这两个字段直接影响最小启停时间约束的边界处理漏掉就会导致第一个时段的启停决策违反物理限制。3.2 用 Pyomo 构建 SCUC 模型并调用求解器PuLP 适合教学但工程上我更推荐 Pyomo因为它对稀疏矩阵和复杂约束的表达更自然而且能无缝切换求解器。下面是一个 Pyomo 版本的骨架import pyomo.environ as pyo model pyo.ConcreteModel() model.G pyo.Set(initializegen[gen_id].tolist()) model.T pyo.Set(initializeload[hour].tolist()) model.u pyo.Var(model.G, model.T, domainpyo.Binary) model.v pyo.Var(model.G, model.T, domainpyo.Binary) model.w pyo.Var(model.G, model.T, domainpyo.Binary) model.P pyo.Var(model.G, model.T, domainpyo.NonNegativeReals) # 参数用字典挂上去 Pmax gen.set_index(gen_id)[Pmax].to_dict() Pmin gen.set_index(gen_id)[Pmin].to_dict() b gen.set_index(gen_id)[b].to_dict() SU gen.set_index(gen_id)[SU].to_dict() D load.set_index(hour)[demand].to_dict() def obj_rule(m): return sum(b[i] * m.P[i,t] SU[i] * m.v[i,t] for i in m.G for t in m.T) model.obj pyo.Objective(ruleobj_rule, sensepyo.minimize) def balance_rule(m, t): return sum(m.P[i,t] for i in m.G) D[t] model.balance pyo.Constraint(model.T, rulebalance_rule) def pmin_rule(m, i, t): return m.P[i,t] Pmin[i] * m.u[i,t] model.pmin pyo.Constraint(model.G, model.T, rulepmin_rule) def pmax_rule(m, i, t): return m.P[i,t] Pmax[i] * m.u[i,t] model.pmax pyo.Constraint(model.G, model.T, rulepmax_rule) solver pyo.SolverFactory(cbc) results solver.solve(model, teeFalse) print(results.solver.termination_condition)Pyomo 的好处是约束可以按集合索引批量生成model.balance一行就覆盖了所有时段。求解器用 CBC 时如果模型规模大记得在solve里加options{sec: 300}限制单次求解时间避免卡死。3.3 求解结果校验功率平衡与爬坡是否真的满足求解器返回 optimal 不代表结果物理上合理。我每次都会做三项校验第一逐时段检查 sum(P) 是否等于负荷误差超过 1e-4 就说明平衡约束没生效第二检查相邻时段 P 的差值是否超过爬坡限值第三检查最小启停时间是否被违反。下面这段校验代码可以直接复用P_val {(i,t): pyo.value(model.P[i,t]) for i in model.G for t in model.T} u_val {(i,t): pyo.value(model.u[i,t]) for i in model.G for t in model.T} for t in model.T: total sum(P_val[i,t] for i in model.G) assert abs(total - D[t]) 1e-4, fBalance violated at {t} for i in model.G: for t in list(model.T)[1:]: ramp P_val[i,t] - P_val[i,t-1] assert ramp RU[i] 1e-4, fRamp up violated gen {i} at {t}如果校验失败先看求解器日志里的约束松弛量再回头检查约束的索引是否漏了边界时段。边界时段是 SCUC 最容易出错的地方因为 t0 没有前一个时段tT 没有后一个时段很多约束在边界处需要特殊处理。4. 避坑与排查SCUC 求解中最容易翻车的五个地方4.1 现象求解器返回 infeasible但负荷和机组容量明明够原因通常是备用约束和爬坡约束在边界时段冲突。比如第一个时段所有机组都处于停机状态但最小开机时间要求它们不能立刻启动导致功率平衡无解。解决方法是检查init_state和init_hours是否合理必要时在第一个时段放宽最小启停时间约束或者增加一台必开机组作为 slack。4.2 现象求解时间随机组数量指数增长24 时段跑几个小时这是 MILP 的典型病。原因可能是分支定界树太大或者模型里有大量对称解。解决手段有三个一是给机组按成本排序加一个「同类型机组启停顺序」的对称破缺约束二是把二次成本分段线性化时不要分太多段3 到 5 段足够三是设置 MIP gap比如 0.5% 就接受不必追求全局最优。4.3 现象结果里机组频繁启停启动成本却很低检查启动成本是否只写了热启动成本冷启动成本被漏掉。另外如果最小停机时间约束写错机组可以停一个时段就重新启动优化器会利用这个漏洞。正确的做法是同时约束最小开机和最小停机时间并且用启动类型变量区分热启动和冷启动。4.4 现象网络安全约束加上去后模型不可行直流潮流约束里 PTDF 矩阵的维度必须和节点注入功率匹配。常见错误是机组挂错节点或者 PTDF 矩阵没有包含平衡节点的处理。我一般会先用一个纯 ED 模型验证 PTDF 矩阵的正确性给定一组注入功率手算线路潮流和模型输出对比误差在 1e-6 以内才继续。4.5 现象经济调度结果和 SCUC 里的出力不一致这通常是因为 SCUC 里的成本函数做了线性化而 ED 用了原始二次成本。两者本来就会有偏差但如果偏差超过 5%说明线性化的分段点选得不好。解决方法是根据机组的经济出力区间重新划分分段或者在 ED 里也用同样的线性化成本。5. 进阶技巧用滚动时域和热启动把 SCUC 跑得更快更稳当你把基础版跑通之后下一步就是让它能应对实际规模。我常用的两个技巧是滚动时域和热启动。滚动时域的思路是把 24 小时切成多个窗口比如每次优化 6 小时但只执行前 1 小时然后窗口向前滑动。这样做的好处是每个子问题规模小求解快而且能自然处理负荷预测误差。代价是全局最优性会打折扣但对于日前计划来说滚动时域的鲁棒性往往比全局最优更重要。热启动则是把上一次求解的 u 变量作为初始可行解传给求解器。Pyomo 里可以通过model.u[i,t].value previous_u[i,t]来设置初始值Gurobi 和 CBC 都支持 MIP start。实测下来热启动能把求解时间缩短 30% 到 50%尤其是当负荷曲线变化不大时。还有一个容易被忽略的技巧是约束松弛。实际系统中备用约束和网络安全约束偶尔会冲突导致模型不可行。与其让整个求解失败不如给关键约束加一个松弛变量并在目标函数里给松弛变量一个高惩罚。这样求解器会告诉你哪条约束最紧而不是直接罢工。# 给备用约束加松弛变量 model.slack_reserve pyo.Var(model.T, domainpyo.NonNegativeReals) def reserve_rule(m, t): return sum(m.u[i,t] * Pmax[i] for i in m.G) D[t] R[t] - m.slack_reserve[t] model.reserve pyo.Constraint(model.T, rulereserve_rule) # 目标函数里加惩罚 model.obj pyo.Objective( exprsum(b[i] * model.P[i,t] for i in model.G for t in model.T) sum(SU[i] * model.v[i,t] for i in model.G for t in model.T) 10000 * sum(model.slack_reserve[t] for t in model.T), sensepyo.minimize )惩罚系数取 10000 是因为它要远大于任何一台机组的运行成本这样求解器会优先消除松弛。但也不能太大否则数值稳定性会变差一般取总成本量级的 100 倍左右。最后说一个我自己的习惯每次改完模型先拿 3 机 4 时段的迷你算例跑一遍确认约束逻辑没写反再上 100 机 24 时段的真实算例。这个习惯帮我省下了无数次「跑了一晚上发现结果全错」的后悔药。SCUC 这玩意儿模型里一个符号写反结果可能看起来还挺合理但成本会低得离谱等你发现的时候已经浪费了半天。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
