用IPOPT求解电力系统经济调度：从建模到调参实战指南

简介：这是一份基于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 = D，D 为该时段负荷）、机组出力上下限（P_min ≤ P_i ≤ P_max）、爬坡约束（多时段时体现：P_i^t − P_i^(t−1) ≤ RU_i，P_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 环境准备：一条命令装出可执行的 ipopt

Pyomo 只是建模层，真正求解的是 ipopt 可执行文件。只pip install pyomo而不装 ipopt 二进制，运行时 SolverFactory 会报找不到求解器。我一般用 conda 一次性解决依赖：

conda create -n ed python=3.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, bounds=lambda 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(expr=sum(cost_rule(m, i) for i in m.I), sense=pyo.minimize) # 功率平衡约束：总出力 = 负荷 def balance_rule(m): return sum(m.P[i] for i in m.I) == D m.BAL = pyo.Constraint(rule=balance_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, tee=True) # tee=True 可以看到迭代日志 m.display()

逻辑说明：决策变量m.P是一维数组，每个元素代表一台机组的出力，用 bounds 参数直接约束上下限，不需要额外写不等式约束。目标函数是二次成本求和。功率平衡约束写成等式，这是 ED 的核心约束。solve(m, tee=True)会把 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, bounds=P_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(rule=obj_rule, sense=pyo.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, rule=bal_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, rule=ramp_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, rule=ramp_dn_rule) solver = pyo.SolverFactory("ipopt") solver.options = {"tol": 1e-8, "max_iter": 2000} result = solver.solve(m, tee=False)

逻辑说明：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(direction=pyo.Suffix.IMPORT) # 声明要导入对偶值 result = solver.solve(m, tee=False) for t in range(1, 4): # 看前三个时段 lam = m.dual[m.BAL[t]] print(f"t={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(direction=pyo.Suffix.IMPORT)是让 Pyomo 在求解完成后把对偶乘子读回模型，否则m.dual[m.BAL[t]]拿不到值。循环里取前三个时段作抽查，把每台机组的边际成本打印出来和 λ 对比。在没触界的情况下，两者应该几乎相等，偏差在tol数量级内。

参数说明：dual 值的准确性依赖tol设置，tol 越严，λ 和边际成本的对碰越接近。如果只关心经济信号而不做高精度校验，acceptable_tol=1e-6就够用。对碰之后再去读目标函数数值才有意义，否则可能模型求解成功但经济性解释全错。

4. 参数怎么设：IPOPT 的必调参数与边界值

IPOPT 的参数很多，真正影响 ED 收敛的就那几个。下面按“先用默认、再按现象调”的顺序讲。

4.1 收敛判定：tol 与 acceptable_tol

IPOPT 默认tol=1e-8，这个值对大多数 ED 模型过于严格。ED 的出力单位是 MW，经济信号能精确到小数点后两位就足够调度员用了，把收敛残差压在 1e-8 只会让迭代在最后阶段反复逼近，白白增加迭代次数。我的通常做法是：

acceptable_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_method=gradient-based自动缩放，但自动缩放不是万能的。当目标函数和约束的量级差超过三个数量级时，我倾向于手动归一数据，而不是依赖求解器内部缩放。

初值的影响同样被低估。Pyomo 不设初值时，变量默认从 0 开始。对爬坡约束密集的 24 时段模型，0 初值意味着每个时段都远离可行域，IPOPT 需要先花大量迭代寻找可行点。我的习惯是先把所有出力初值设为各机组 Pmin，或者用上一个时段的解作为下一个时段的初值，这个习惯在扩展机组组合循环时特别管用。

mu_strategy也值得说一下。默认是adaptive，适合大多数凸 ED；如果遇到目标函数在山谷附近反复震荡，可以改成monotone并设置mu_init=0.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 的骨架搭得扎实。希望帮到你。

本文还有配套的精品资源，点击获取