简介这是一份基于IPOPT求解电力系统经济调度问题的MATLAB实现源码面向电力系统工程师、研究人员及学习非线性优化的学生。资源将经济调度中机组出力分配、功率平衡与成本最小化等约束转化为可直接求解的优化模型涵盖目标函数、雅可比矩阵、非线性约束及初始点定义等核心模块并通过紧凑接口完成IPOPT调用。压缩包共7个文件全部为.m脚本整体仅5KB代码精简、结构清晰适合快速阅读与二次开发。当前已有139人学习可用于理解内点法在电力调度中的实际落地过程也可作为课程设计或论文仿真的参考基底。1. 经济调度遇见 IPOPT这个组合为什么值得你花一晚上调通多时段电力系统经济调度里最磨人的不是把目标函数写出来而是让内层 ED 子问题在成百上千条约束下不翻车地收敛。用 scipy 的 SLSQP 调过的人都知道那种黑匣子式迭代很容易卡在可行域边缘转向商用求解器又要处理授权。IPOPT 是开源的内点法优化器对二次成本、非线性网损、大规模稀疏结构都很对口是电力系统经济调度ED方向常见的求解内核ED_using_IPOPT 这类例程名里挂的核心求解器十有八九就是它。这篇文章面向手里已有模型、想把成本曲线和约束落到可复现 Python 代码的从业者。我按“建模 → 最小代码 → 参数边界 → 踩坑”的顺序讲最后给一套能自检的验证习惯照做基本能跑通而且能判断结果对不对。2. 把 ED 模型喂给 IPOPT 之前目标函数、约束与维度设计在写SolverFactory(ipopt)之前先把数学模型写清楚。这一步能省掉后面八成的问题。2.1 经济调度的数学形式从机组成本到等式约束经典 ED 问题的目标函数是各机组发电成本之和。最常见的成本函数是二次型最小化 Σ_i (a_i · P_i² b_i · P_i c_i)其中 P_i 是机组 i 的有功出力a_i、b_i、c_i 由燃料成本曲线拟合得到。约束分四类功率平衡等式忽略网损时为 Σ P_i DD 为该时段负荷、机组出力上下限P_min ≤ P_i ≤ P_max、爬坡约束多时段时体现P_i^t − P_i^(t−1) ≤ RU_iP_i^(t−1) − P_i^t ≤ RD_i以及可选网损约束。网损的处理是 ED 从线性走向非线性的关键节点。常见的工程近似是 B 系数法L(P) Σ_i Σ_j P_i · B_ij · P_j Σ_i B_0i · P_i B_00功率平衡等式随之变为 Σ P_i D L(P)。B 系数是二次形式和二次成本叠加后整个模型变成真正的非线性规划。这个场景是 IPOPT 的发挥区因为它本来就是为了大规模稀疏 NLP 设计的。还有个容易被忽略的判断忽略网损时ED 是带线性约束的凸二次规划拉格朗日函数的 KKT 条件可以直接写成“等微增率”公式对所有不触界的机组满足 b_i 2a_i·P_i λ其中 λ 就是平衡约束的影子价格也是系统边际成本。这个性质后面当校验用比盯着目标函数值靠谱得多。2.2 为什么这一步选 IPOPT和 SLSQP、LP、商用求解器的对比选求解器要看问题的规模和可导性。小型单时段 ED 直接用 scipy 的 SLSQP 能解但一旦扩到 24 时段、100 台机组SLSQP 每次迭代都要处理稠密近似 Hessian内存和时间都跟不上。线性规划需要预先讲故事成本曲线分段线性化分段多了 LP 规模的起飞而且分段点处的最优解容易贴着拐点经济性解释变得别扭。商用 QCP 求解器自然强但授权和对外替代性常是项目边界条件。IPOPT 走原始-对偶内点法每次迭代求解一个关于 Hessian、Jacobian 和约束右端项的 KKT 线性系统。它依赖的是稀疏线性代数组件而不是稠密矩阵所以对“很多约束、很多变量但非零元素稀疏”的 ED 模型特别合适。ED 场景里 Hessian 是每台机组一个对角块Jacobian 按平衡约束、爬坡约束、上下限约束展开后也高度结构化正好是这种算法的舒适区。边界也要说清楚IPOPT 是连续 NLP 求解器不能直接处理机组启停的 0/1 变量。ED_using_IPOPT 这类例程通常是在给定启停状态的前提下求连续出力机组组合那一层的整数决策要交给外层 MIP 或启发式框架。忽略这个边界后面一定会踩坑。2.3 维度选择和稀疏结构决策变量顺序会影响收敛吗用 Pyomo 建模时变量顺序由索引方式决定一般不需要过度优化。但如果要迁移到 cyipopt 直传梯度变量排列就要讲究。我的习惯是把 P 铺成一维数组索引取 t*G g也就是“先时段、后机组”。原因是功率平衡约束、爬坡约束都是按时段展开的这种排列让 Jacobian 的非零块沿着时段方向形成带状结构稀疏线性求解器的符号分解更顺。维度设计上还有一个常被忽略的点量纲。出力用 MW负荷用 MW成本系数用元/MW²h 和元/MWh这是 ED 内部的统一约定。一旦把负荷写成 kW、价格又按万元计问题数学上等价但目标量级和约束量级差距拉大IPOPT 内部缩放的负担就变重迭代曲线容易震荡。我的做法是先在模型外把所有数据归一到同一个基准最简单的基准就是负荷峰值 MW 数求解完再反标幺回来省掉一堆缩放玄学。约束的顺序也可以稍微留个心。IPOPT 默认的线性求解器 mumps 对行排序有自己的策略一般情况下不需要人工干预。但如果你在 Pyomo 里用ConstraintList动态追加约束注意先加平衡约束、再加爬坡约束不要让同一批约束在列表里穿插出现否则调试时可读性很差出错时定位也麻烦。3. 用 Pyomo IPOPT 跑通最小 ED从 3 机 1 时段到 24 时段这一章给一套可以直接抄的例程。目标是从最小模型跑通再逐步加维度。3.1 环境准备一条命令装出可执行的 ipoptPyomo 只是建模层真正求解的是 ipopt 可执行文件。只pip install pyomo而不装 ipopt 二进制运行时 SolverFactory 会报找不到求解器。我一般用 conda 一次性解决依赖conda create -n ed python3.11 -y conda activate ed conda install -c conda-forge pyomo ipopt -y which ipopt倒数第二行会同时拉入 pyomo 和 ipopt 可执行文件mumps 线性求解器也会作为 ipopt 的依赖一并装好。最后一行which ipopt用来确认二进制在 PATH 里只要能看到路径Pyomo 就能找到。如果看不到优先检查安装源而不是急着改代码。参数说明Python 版本选当前稳定版即可conda-forge 里的 ipopt 包会匹配可用的 mumps 版本不需要自己编译。pip 路线也可以但 pip 上的 ipopt 包不一定自带求解器二进制这是新手最容易遇到的环境坑。3.2 最小 3 机组 1 时段 ED完整 Python 代码先拿出最小模型3 台机组、1 个时段、没有爬坡、没有网损。目标只是验证链路通不通。import pyomo.environ as pyo # 机组参数a(元/MW^2/h), b(元/MWh), c(元/h), Pmin(MW), Pmax(MW) GEN [ {a: 0.006, b: 4.5, c: 150.0, pmin: 50.0, pmax: 300.0}, {a: 0.012, b: 3.8, c: 260.0, pmin: 30.0, pmax: 200.0}, {a: 0.004, b: 6.2, c: 80.0, pmin: 20.0, pmax: 150.0}, ] D 400.0 # 时段负荷MW m pyo.ConcreteModel() m.I pyo.RangeSet(len(GEN)) # 决策变量机组出力带上界和下界 m.P pyo.Var(m.I, boundslambda m, i: (GEN[i-1][pmin], GEN[i-1][pmax])) # 目标函数总发电成本最小 def cost_rule(m, i): g GEN[i-1] return g[a] * m.P[i]**2 g[b] * m.P[i] g[c] m.COST pyo.Objective(exprsum(cost_rule(m, i) for i in m.I), sensepyo.minimize) # 功率平衡约束总出力 负荷 def balance_rule(m): return sum(m.P[i] for i in m.I) D m.BAL pyo.Constraint(rulebalance_rule) # 用 Pyomo 调用 ipopt solver pyo.SolverFactory(ipopt) solver.options { tol: 1e-8, max_iter: 500, acceptable_tol: 1e-6, print_level: 5, } result solver.solve(m, teeTrue) # teeTrue 可以看到迭代日志 m.display()逻辑说明决策变量m.P是一维数组每个元素代表一台机组的出力用 bounds 参数直接约束上下限不需要额外写不等式约束。目标函数是二次成本求和。功率平衡约束写成等式这是 ED 的核心约束。solve(m, teeTrue)会把 IPOPT 的迭代日志打到终端第一次跑建议开 tee能直接看到收敛过程。参数说明tol是 IPOPT 的最优性容忍度控制 KKT 条件的残差大小1e-8 对这个 3 机模型很宽松max_iter是最大迭代次数防止模型不可行时无限跑acceptable_tol是“可接受停止”的阈值配合 print_level 能看到它什么时候提前收住。这个 3 机问题理论最优出力大约在 P1 194.44 MW、P2 126.39 MW、P3 79.17 MW系统边际成本约 6.83 元/MWh跑出来的结果应该和这个对照在容许范围内。3.3 扩展到 24 时段负荷曲线与爬坡约束单时段通了接下来把时间维加上。这个改动会暴露初值、爬坡和索引三个新的问题。import numpy as np T 24 load_profile [400 60 * np.sin(2 * np.pi * (t - 6) / T) for t in range(T)] # 实际工程中用 numpy 或 pandas 读 SCADA 数据这里用简单的正弦近似演示 RU [120.0, 100.0, 80.0] # 每台机组的上爬坡上限 MW/h RD [120.0, 100.0, 80.0] # 下爬坡上限 MW/h m pyo.ConcreteModel() m.G pyo.RangeSet(len(GEN)) m.T pyo.RangeSet(T) # 时段 1..24 def P_bounds(m, g, t): return (GEN[g-1][pmin], GEN[g-1][pmax]) m.P pyo.Var(m.G, m.T, boundsP_bounds) def obj_rule(m): return sum(GEN[g-1][a] * m.P[g, t]**2 GEN[g-1][b] * m.P[g, t] GEN[g-1][c] for g in m.G for t in m.T) m.COST pyo.Objective(ruleobj_rule, sensepyo.minimize) # 每个时段都要满足功率平衡 def bal_rule(m, t): return sum(m.P[g, t] for g in m.G) load_profile[t-1] m.BAL pyo.Constraint(m.T, rulebal_rule) # 爬坡约束相邻时段出力差不能超过限值 def ramp_up_rule(m, g, t): if t 1: return pyo.Constraint.Skip # 第一个时段没有前一时段 return m.P[g, t] - m.P[g, t-1] RU[g-1] m.RAMP_UP pyo.Constraint(m.G, m.T, ruleramp_up_rule) def ramp_dn_rule(m, g, t): if t 1: return pyo.Constraint.Skip return m.P[g, t-1] - m.P[g, t] RD[g-1] m.RAMP_DN pyo.Constraint(m.G, m.T, ruleramp_dn_rule) solver pyo.SolverFactory(ipopt) solver.options {tol: 1e-8, max_iter: 2000} result solver.solve(m, teeFalse)逻辑说明m.P变成二维变量(机组, 时段)目标函数要遍历所有机组和所有时段。负荷平衡约束从一条变成 24 条用Constraint(m.T)一次写入。爬坡约束引用前一时段的出力所以t 1要跳过否则会产生一个不存在的P[g, 0]。参数说明max_iter从 500 提到 2000是要给约束增多后的迭代留余量但也不要盲目放大后面遇到不收敛时要先怀疑可行域而不是迭代次数。模型里没给初值IPOPT 会从 0 开始对爬坡约束比较多的模型0 初值意味着初始点离可行域远这时可以把m.P[g, t].set_value(GEN[g-1][pmin])放在 solve 前往往能少走很多弯路。3.4 结果校验影子价格与等微增率对碰跑通只是第一步还要验证结果合理。ED 的最优性有个硬指标所有不触界的机组边际成本 b_i 2a_i·P_i 应该等于平衡约束的影子价格 λ。这个 λ 正好是功率平衡约束的 dual 值。m.dual pyo.Suffix(directionpyo.Suffix.IMPORT) # 声明要导入对偶值 result solver.solve(m, teeFalse) for t in range(1, 4): # 看前三个时段 lam m.dual[m.BAL[t]] print(ft{t} lambda{lam:.4f} 元/MWh) for g in range(1, len(GEN)1): P_val pyo.value(m.P[g, t]) mgc GEN[g-1][b] 2 * GEN[g-1][a] * P_val print(f 机组{g} P{P_val:.2f} marginal_cost{mgc:.4f})逻辑说明Suffix(directionpyo.Suffix.IMPORT)是让 Pyomo 在求解完成后把对偶乘子读回模型否则m.dual[m.BAL[t]]拿不到值。循环里取前三个时段作抽查把每台机组的边际成本打印出来和 λ 对比。在没触界的情况下两者应该几乎相等偏差在tol数量级内。参数说明dual 值的准确性依赖tol设置tol 越严λ 和边际成本的对碰越接近。如果只关心经济信号而不做高精度校验acceptable_tol1e-6就够用。对碰之后再去读目标函数数值才有意义否则可能模型求解成功但经济性解释全错。4. 参数怎么设IPOPT 的必调参数与边界值IPOPT 的参数很多真正影响 ED 收敛的就那几个。下面按“先用默认、再按现象调”的顺序讲。4.1 收敛判定tol 与 acceptable_tolIPOPT 默认tol1e-8这个值对大多数 ED 模型过于严格。ED 的出力单位是 MW经济信号能精确到小数点后两位就足够调度员用了把收敛残差压在 1e-8 只会让迭代在最后阶段反复逼近白白增加迭代次数。我的通常做法是参数常用值说明tol1e-6 到 1e-8最优性残差门槛越严迭代越多acceptable_tol1e-6 或 1e-7可接受停止阈值acceptable_iter15连续多少步满足可接受条件就停max_iter5000硬性上限防止不可行模型死循环print_level5控制迭代日志详细程度排查时调到 8acceptable_tol是 IPOPT 的一个“后悔药”机制如果连续acceptable_iter步的残差都在可接受范围内就提前结束即使没有达到tol。对大规模 ED肉眼观察历史负荷变化时1e-6 和 1e-8 的解在 MW 层面看不出差别但求解时间可能差好几倍。4.2 线性求解器mumps 与 ma57 的取舍IPOPT 每步迭代都要解一个大型稀疏线性系统线性求解器决定了下限。conda 默认安装的 mumps 对几万变量的问题表现不错我遇到的绝大多数 ED 模型都用它收住。模型到几十万变量以上、或 Hessian 结构非常稠密时mumps 的符号分解会变成瓶颈这时可以尝试solver.options[linear_solver] ma57但注意 ma57 是 HSL 软件包里的求解器需要额外授权不是所有二进制都内置。换成 ma57 前先确认当前 ipopt 编译时是否包含了它否则 SolverFactory 会报错。常见做法是先把问题规模降下来验证再考虑换求解器不要一开始就在线性代数层面较劲。4.3 缩放与初值这就是 IPOPT 的玄学来源同样的模型把成本系数单位从元换成万元目标函数值从几万变成几十收敛速度可能天差地别。IPOPT 默认开启nlpscaling_methodgradient-based自动缩放但自动缩放不是万能的。当目标函数和约束的量级差超过三个数量级时我倾向于手动归一数据而不是依赖求解器内部缩放。初值的影响同样被低估。Pyomo 不设初值时变量默认从 0 开始。对爬坡约束密集的 24 时段模型0 初值意味着每个时段都远离可行域IPOPT 需要先花大量迭代寻找可行点。我的习惯是先把所有出力初值设为各机组 Pmin或者用上一个时段的解作为下一个时段的初值这个习惯在扩展机组组合循环时特别管用。mu_strategy也值得说一下。默认是adaptive适合大多数凸 ED如果遇到目标函数在山谷附近反复震荡可以改成monotone并设置mu_init0.1让障碍参数单调下降换取更可靠的步长代价是迭代次数增加。不要同时调mu_init和tol两个参数一次只改一个否则出问题不知道归因于谁。4.4 遇到 0/1 变量时IPOPT 不是 MIP这是 ED_using_IPOPT 方向最容易产生误解的地方。把机组启停变量松弛成 [0,1] 连续变量直接丢给 IPOPT确实能跑出 0.9999 这种数值但它并不代表组合最优λ 也会被污染。工程上正确姿势是把它当作子问题外层用分支定界、启发式或 MIP 决定启停内层用 IPOPT 求解给定启停下的连续出力。IPOPT 的输出可以作为 MIP 分支节点的下界也可以作为调度计划校验工具但不要指望它直接回答“开几台机组”。5. 避坑清单经济调度用 IPOPT 翻车的 5 个高频场景这一章把我实际排查过的典型问题写出来每条都按现象、原因、解决展开。5.1 Restoration phase failed先怀疑模型可行性还是初值现象IPOPT 迭代到一半日志里报Restoration phase failed或者直接退出码显示EXIT: Restoration Failed。模型不大但它在无数个点上找可行方向都没找到。原因最常见是可行域为空。比如负荷曲线里某时段 D 小于所有机组 Pmin 之和或者爬坡约束与上下限组合后不存在满足条件的出力序列。也有可能是初值离可行域太远内点法在恢复阶段用尽迭代。解决先把等式约束临时改成不等式例如把Σ P D换成Σ P D - slack给每个平衡约束加一个松弛变量让模型总能找到可行点。如果松弛量很大才可行说明原模型数据有矛盾。顺着时段查负荷和机组上限基本能找到问题点。5.2 Jacobian / Hessian 无法求导表达式可导性出事现象报类似Cannot evaluate Jacobian、Error in evaluator的信息或者求解过程在第一二次迭代就中断。原因Pyomo 表达式里掺了不可导的 Python 函数比如math.floor、round、min直接作用在变量上或者在规则里用if m.P[i] 0这类依赖变量值的逻辑分支。IPOPT 要求目标函数和约束在可行域内可导这类写法破坏了光滑性。解决把min(a, b)、max(a, b)改成辅助变量加上限/下限约束。分段线性成本里的绝对值也要拆成正负偏差变量不要直接在表达式里写abs。这是 NLP 建模的基本纪律但每次踩坑都会有人再犯。5.3 分段线性成本让 IPOPT 越跑越慢现象目标函数由二次成本改成分段线性成本后IPOPT 迭代步数明显变多收敛后目标值还在相邻分段点之间抖动。原因分段线性函数在分段点处不可导IPOPT 虽然能通过某种近路径推进但 Hessian 信息在这一带失真内点法步长被迫缩小迭代变成蚂蚁搬家。解决如果成本曲线本来就是二次拟合不要换成 PWL。如果来源数据就是分段报价先用光滑函数拟合再交给 IPOPT或引入 SOS2 约束配合 MIP 求解器处理。非要在 IPOPT 里硬解就把分段点数量减少到最小并且保证初值不落在分段点附近。5.4 同一份模型换数据后结果乱跳现象模型和数据格式都没变只是把负荷从 400 MW 换成 405 MW结果某台机组出力从 100 MW 跳到 180 MW甚至换一台机组顶上去λ 也跟着突变。原因很多 ED 的最优解在机组成本曲线交叉点附近本来就对负荷敏感出力跳变可能是真实的最优响应不一定是 bug。但更常见的原因是量纲没统一负荷写成 405000 kW成本系数还是按 MW 标定的导致 scaling 崩溃IPOPT 的解精度其实不合格。解决统一用 MW 和元/MWh 建模求解前打印一次所有张量数据的 min/max确认量级一致。再不行用上一章的 shadow price 校验方法看触界机组和非触界机组的 λ 是否自洽。数据变更后不要只看目标值先看 λ 和边际成本的对碰结果。5.5 ipopt 找不到环境问题排在第一位现象Pyomo 求解时报No such file or directory: ipopt或Solver (ipopt) did not exit normally代码本身换个环境就能跑。原因最常见是只装了 pyomo没装 ipopt 可执行文件或者用 pip 装的 pyomo 包并没有携带求解器二进制而机器上又没独立安装 ipopt。解决回到 3.1 的which ipopt确认 PATH 里有可执行文件。如果没有用 conda-forge 的 ipopt 包补上。还不行就在SolverFactory(ipopt, executable/path/to/ipopt)里显式指向路径。每次换机器、换 CI 环境第一件事就是检查这条命令别急着去调模型参数。6. 把 IPOPT 用出效率的三个习惯校验、warmstart 与解析解对碰最后讲三个我每天都会用到的操作它们让 IPOPT 从“能跑”变成“跑得可信”。6.1 永远留一条解析可解的对照用例我在做任何 ED 模型前都会先留一个 2 到 3 机、单时段、无网损的用例因为那个问题可以用等微增率手算。模型升级到 24 时段、加网损、加爬坡后把对照用例的输入套进新代码里如果输出和手算结果偏离超过 1e-3 级别说明引入新约束时有 bug。这个对照用例不用删放在测试目录里当回归测试每次改代码都跑一遍。6.2 用上一轮的解做 warmstart机组组合外循环里同一个机组的启停状态往往只有一两个变量变化连续出力解的变化很小。把上一轮 IPOPT 的最优出力值设为下一轮的初值求解时间能降一半以上。for g in m.G: for t in m.T: m.P[g, t].set_value(prev_solution[g-1][t-1])逻辑说明prev_solution是上一轮结果结算后存的二维数组按机组和时段索引。对爬坡约束较多的模型这个初值天然满足爬坡可行域的大部分IPOPT 不需要花阶段一重新找可行点。参数说明如果启停状态发生大变化比如某台机组从关到开这台机组的初值仍要用 Pmin 给不能用旧出力否则会引入不可行初值。6.3 影子价格要写成回归报告ED 的结果不应只存目标值和出力矩阵还要把每个时段平衡约束的 λ、各机组边际成本、约束是否触界一起落盘。拿 λ 和边际成本差值做一个回归指标差值超过 1e-2 就报警。这个指标能同时捕捉量纲错误、约束写错和收敛精度不足三类问题。我自己的教训是有一回换了负荷数据后目标值看着很合理λ 却比边际成本高出一截排查了半天发现是某个时段的负荷被当成 kW 读入模型本身没错是数据预处理漏了换算。自那以后所有 ED 模型里必须有一份 dual 落盘先看 λ 再看出力。以上这组流程从最小模型到 24 时段从参数设置到结果校验够把一个 ED_using_IPOPT 的骨架搭得扎实。希望帮到你。本文还有配套的精品资源点击获取