简介这是一套面向电力系统工程师与科研人员的IPOPT经济调度实现示例围绕电力系统经济调度问题提供完整建模与求解方案。压缩包内含7个MATLAB脚本文件整体仅5KB涵盖目标函数、约束条件、雅可比矩阵、初始点设置及全局变量定义等模块便于理解代码结构与快速复用。项目通过内点法求解机组出力最优分配在满足功率平衡、机组出力限制等约束下最小化总发电成本适合具备一定优化基础和MATLAB使用经验的读者入门实践。已有139人学习下载可作为学习IPOPT接口调用、掌握电力经济调度建模思路的参考范例。1. 电力系统经济调度的 IPOPT 解法先跑通一个最小模型再谈别的负荷预测曲线每小时刷新15 分钟后下一时段的机组出力方案就要落到调度台。这就是电力系统经济调度ED的典型场景给定负荷算出一组机组出力让总发电成本最低同时满足功率平衡和机组自身限制。这个优化问题一旦带上火电机组的非线性成本特性、网损修正和爬坡约束就成了一个地道的非线性规划NLP而 IPOPT 正是开源求解器里处理这类问题最常用的那一个。标题里那一串看起来像仓库文件夹名的字符本质上就是用 IPOPT 求解电力系统经济调度这个最小工程骨架。如果你正在接手调度系统的算法模块、在做求解器选型或者想把手里的商业求解器换掉一部分这篇笔记能让你从模型写法一路走到参数调优和排错。全文代码围绕一个三机系统展开单母线模型几秒钟出结果足够把原理和坑都讲清楚。2. 把 ED 写成 IPOPT 能解的 NLP目标函数、约束与网损模型2.1 为什么选 IPOPT从 ED 的非线性结构看求解器选型经济调度的经典目标函数来自火电机组的煤耗特性曲线。工程上最常用的是二次函数形式F_i(P_i) a_i P_i^2 b_i P_i c_i其中 P_i 是第 i 台机组的出力MWa_i、b_i、c_i 是机组成本系数。如果只保留 b_i P_i 项模型退化成线性规划LP用单纯形法就能解但二次项描述的是机组在高低负荷区间的煤耗率差异忽略它会让高负荷机组的边际成本被低估调度结果在经济性上明显偏差。加入二次项之后问题就从 LP 变成二次规划QP如果再考虑阀点效应在二次项上加一个正弦修正项目标函数直接变成非凸 NLP。这就是选求解器的关键分岔点。Gurobi、CPLEX 这类商业求解器处理凸 QP 很快但遇到非凸目标或者非线性约束比如网损 B 系数公式它们要么不支持要么性能暴跌。而 IPOPT 用的是内点法对 NLP 的支持是原生且完备的目标函数和约束只需要能求一阶、二阶导数凸不凸它都能迭代只是收敛到局部最优还是全局最优的区别。对实际 ED 场景来说机组成本函数大多是凸的即便加阀点项也只是轻度非凸IPOPT 配合多初始点尝试已经足够。在开源阵营里IPOPT 的生态也最成熟。Pyomo 里一行 SolverFactory(ipopt) 就能调起来和 AMPL、GAMS、JuMP 都有接口。如果你在学术论文里看到用 IPOPT 求解 ED那基本是默认动作。2.2 目标函数、功率平衡与机组限制ED 的三层标准结构ED 的标准数学形式可以拆成三层。第一层是目标函数最小化所有机组的发电成本之和min Σ (a_i P_i² b_i P_i c_i)第二层是等式约束也就是功率平衡Σ P_i P_load P_lossP_load 是系统总负荷P_loss 是输电网损。这个约束必须严格满足哪怕相差 0.1 MW在实际运行中都会造成频率偏差。IPOPT 处理这类等式约束的手段是将其纳入拉格朗日函数通过内点障碍项和牛顿迭代逐步逼近可行解最终收敛时约束残差会小于你设定的 tol 参数。第三层是不等式约束。最基本的包括机组出力上下限P_min ≤ P_i ≤ P_max如果做的是多时段调度还要加上爬坡约束-R_i ≤ P_i(t) - P_i(t-1) ≤ R_i以及系统备用约束Σ P_max_i ≥ P_load P_reserve在 IPOPT 内部不等式约束不会显式参与迭代而是通过引入松弛变量和障碍项把一个带不等式约束的 NLP 转成一系列只含等式约束的子问题。这也是为什么 IPOPT 对不等式约束的初始点比较敏感——如果初始点让某个不等式越界太多障碍项会给出很大的对偶值迭代前期会花大量步数把它拉回来。这一点在第 5 章避坑部分会再展开。2.3 网损的 B 系数近似从忽略到精确的迭代路径实际系统里的网损不能直接忽略。常见做法是 B 系数法把网损写成机组出力的二次型P_loss Σᵢ Σⱼ P_i B_ij P_jB 矩阵通过潮流计算离线求得对 3 机系统是一个 3x3 系数矩阵对上百台机组的系统就是上百阶的稠密矩阵。把 P_loss 代入功率平衡约束后约束本身变成二次等式这对 IPOPT 没有任何额外负担——它本来就是非线性求解器。不过在实现上有一个经验问题B 系数是按某个基准运行点线性化得到的如果你直接把它写进约束得到的调度结果在真实潮流下网损会有偏差。我一般会做 2 到 3 次迭代先按 P_loss 0 求解一组出力用这组出力算网损再代回约束重新求解重复两轮后网损偏差通常可以压到 0.5% 以内。这在工程上比一次性解准得多而且每轮 IPOPT 求解只要几秒。2.4 数据组织与量纲写模型前先准备好这几张表ED 建模翻车最常见的原因不是优化模型本身而是数据和量纲混乱。我习惯把机组参数整理成下表这样的结构再开始写代码。机组a (元/MW²·h)b (元/MWh)P_min (MW)P_max (MW)爬坡 (MW/h)G10.01214.0108520G20.01513.5108020G30.01015.057025这个表格里有三个量纲陷阱成本系数 a 的量级是十的负二次方b 是十的一次方而功率是几十 MW。三者相乘后目标函数值在几千的量级而约束的残差在几百的量级IPOPT 内部的自动缩放scaling通常能处理但如果你的机组规模到了上百台某些 a 系数会变成十的负五次方我建议手动做一次缩放把目标函数整体除以 1000或者把功率基准改成 100 MW。量纲问题不解决后面第 5 章说的Restoration failed和迭代停滞就会频繁来找你。3. 用 Pyomo 调用 IPOPT 求解经济调度最小可复现代码3.1 环境准备安装 IPOPT 与 PyomoPyomo 是 Python 生态里最成熟的优化建模框架它支持将模型自动转换为 NLP 格式并调用 IPOPT。安装用 conda 最省事因为 conda-forge 里已经打包好了 IPOPT 的可执行文件不需要自己编译源码。conda create -n ed python3.11 -y conda activate ed conda install -c conda-forge ipopt pyomo -y安装完成后验证一下 IPOPT 是否能正常调用ipopt -v能打印出版本信息就说明求解器就绪。如果你用的是 Windowsconda-forge 也提供了可用的 ipopt 包但要注意从源码编译的 IPOPT 在 Windows 上容易遇到缺少 libgfortran 的问题能装 conda 版就不要自己编译——这算是血泪经验。3.2 最小 ED 模型三个机组、一个功率平衡约束下面的代码是一个完整的最小 ED 模型。三台机组单母线目标函数为二次成本约束为功率平衡和出力上下限。负荷取 150 MW。import pyomo.environ as pyo # 机组数据a 为二次成本系数b 为一次成本系数Pmin/Pmax 为出力上下限 gen_data { 1: {a: 0.012, b: 14.0, Pmin: 10, Pmax: 85}, 2: {a: 0.015, b: 13.5, Pmin: 10, Pmax: 80}, 3: {a: 0.010, b: 15.0, Pmin: 5, Pmax: 70}, } P_load 150 # MW m pyo.ConcreteModel() # 机组集合 m.G pyo.Set(initializegen_data.keys()) # 出力变量直接绑定上下界避免额外写不等式约束 m.P pyo.Var(m.G, boundslambda m, i: (gen_data[i][Pmin], gen_data[i][Pmax])) # 目标函数最小化总发电成本 m.cost pyo.Objective( exprsum(gen_data[i][a] * m.P[i]**2 gen_data[i][b] * m.P[i] for i in m.G) ) # 功率平衡约束总出力 总负荷 m.balance pyo.Constraint( exprsum(m.P[i] for i in m.G) P_load ) # 调用 IPOPT 求解 solver pyo.SolverFactory(ipopt) solver.options[tol] 1e-8 result solver.solve(m, teeTrue) # 打印结果 for i in sorted(m.G): print(f机组 {i}: P {pyo.value(m.P[i]):.2f} MW) print(f总出力 {sum(pyo.value(m.P[i]) for i in m.G):.2f} MW) print(f总成本 {pyo.value(m.cost):.2f} 元/h)这段代码有几个值得说明的设计。第一变量 m.P 的 bounds 参数直接传入了区间Pyomo 会将其生成为变量边界IPOPT 内部用障碍项处理比显式写 Inequality 约束效率更高。第二功率平衡用的是 而且是单条约束IPOPT 对等式约束的处理比不等式更直接。第三tol 设到 1e-8 而不是默认的 1e-6原因后面会说但先记住这个值在 ED 场景下更好用。运行这段代码你会看到 G1 和 G3 承担大部分出力G2 因为成本系数偏高而略少这是符合直觉的结果。3.3 加网损后的模型改动B 系数约束怎么加进去实际系统不能忽略网损在上一节模型的基础上扩展把 B 系数矩阵和网损表达式加进去import numpy as np # B 系数矩阵三机系统的网损二次型系数 # 注意这里的 B 按 100 MVA 基准折算为标幺值实际需根据潮流计算得到 B np.array([ [0.021, 0.001, 0.002], [0.001, 0.017, 0.001], [0.002, 0.001, 0.023], ]) # 定义网损表达式P_loss P^T B P def loss_expr(m): idx sorted(m.G) P [m.P[i] for i in idx] return sum(P[i] * B[i, j] * P[j] for i in range(len(idx)) for j in range(len(idx))) # 修改功率平衡约束总出力 负荷 网损 m.balance pyo.Constraint( exprsum(m.P[i] for i in m.G) P_load loss_expr(m) ) solver pyo.SolverFactory(ipopt) solver.options[tol] 1e-8 solver.solve(m, teeTrue) for i in sorted(m.G): print(f机组 {i}: P {pyo.value(m.P[i]):.2f} MW) print(f网损 {pyo.value(loss_expr(m)):.2f} MW)关键点在于balance 约束的右侧从常数变成了变量 P 的二次型整个约束变成非线性的。IPOPT 的处理方式是把这个二次型纳入约束雅可比矩阵通过牛顿法迭代求解。相比上一节的纯线性约束这会让每轮迭代的线性系统求解成本变高但对三机系统来说时间可以忽略。还有一个工程细节不要手动把网损固定成一个常数再求解那样的话功率平衡可能无法满足。正确做法是让网损表达式完整地参与建模IPOPT 会在迭代中自动协调各机组出力使得总出力 - 负荷 - 网损 0在收敛时成立。4. IPOPT 参数怎么设6 个影响 ED 结果的旋钮与推荐基准4.1 IPOPT 参数地图先认识这几个旋钮IPOPT 的可调参数超过 100 个但做 ED 求解真正需要关心的就 6 个。下面这张表是实际项目里常用的配置基准。参数名作用ED 场景推荐值tol收敛容忍度KKT 误差小于该值判定收敛1e-6 到 1e-8acceptable_tol可接受的低精度收敛阈值1e-5acceptable_iter连续多少次迭代满足 acceptable_tol 后提前终止5max_iter最大迭代次数5000linear_solver求解 KKT 线性系统所用的线性代数库mumps / ma57 / ma97mu_strategy障碍参数 μ 的更新策略adaptivebound_push初始点往可行域内部推的距离1e-2print_level控制台输出详细程度5调试时这些参数里对 ED 结果影响最大的是 tol、linear_solver 和 mu_strategy。tol 直接决定功率平衡约束最终残差设 1e-6 时150 MW 负荷可能残留 0.00015 MW 的不平衡量设 1e-8 时残差降到 0.0000015 MW。对调度系统来说这个量级差异看似微乎其微但如果你把出力结果送到潮流计算里做验证1e-6 的残差有时会造成节点电压的可见偏差。4.2 按 ED 场景调参精度优先还是速度优先调参的推荐顺序是先定 tol再定 linear_solver最后调 mu_strategy。下面是实际项目里的一套起点配置opt pyo.SolverFactory(ipopt) # 精度ED 对功率平衡残差敏感tol 偏严 opt.options[tol] 1e-8 # 可接受解如果在迭代中途连续 5 次达到 1e-5也算收敛 opt.options[acceptable_tol] 1e-6 opt.options[acceptable_iter] 5 # 线性求解器mumps 开源且稳健追求性能可换 ma57/ma97需另外安装 opt.options[linear_solver] mumps # 障碍参数adaptive 在目标函数比较光滑时收敛更快 opt.options[mu_strategy] adaptive # 迭代上限防止个别病态案例无限迭代 opt.options[max_iter] 5000 opt.options[print_level] 5 res opt.solve(m, teeTrue)这里解释三个关键选择。linear_solver 选 mumps 是最省心的默认项因为它是 IPOPT 自带的开源稀疏线性直接求解器不需要额外的许可证如果你的系统到了几百台机组、约束几万个mumps 的内存占用和求解时间会明显上升这时候换成 ma57 或 ma97HSL 库的求解器能快 2 到 3 倍但需要单独向 HSL 申请学术或商业授权。tol 从 1e-6 收紧到 1e-8 通常会让迭代次数增加 10% 到 30%但对小规模 ED 也就多几十毫秒值得。还要重点提一下不可接受但对结果有微妙影响的参数bound_push。它控制初始点距离变量边界有多远。IPOPT 默认值是 1e-8意味着初始点几乎贴在边界上。如果你的某个机组在解里处于 Pmin 边界初始点贴边界会造成障碍项数值过大前面几十轮迭代都在把变量往可行域内部推。把这个参数调成 1e-3 或 1e-2等于告诉求解器初始点请离边界远一点对 ED 这种大量机组运行在边界附近的场景有奇效。如果遇到迭代不收敛我还习惯把 print_level 调到 5 查看每次迭代的 inf_pr原始残差和 inf_du对偶残差。inf_pr 不降说明是可行性问题inf_du 不降说明是目标函数梯度或 Hessian 的问题。这个判断方向在排错时非常关键。5. 经济调度求解避坑指南IPOPT 常见报错与排查记录5.1 Restoration failed不可行问题还是初始点问题现象求解器运行几十轮迭代后突然报错 Restoration failed!退出码显示 Infeasible_Problem_Detected结果文件是空的。原因有两种可能。第一种是约束真的不可行比如你把负荷设成 200 MW但机组 Pmax 之和只有 170 MW功率平衡方程永远无法满足。第二种是初始点离可行域太远IPOPT 的恢复阶段restoration phase找不到一个可行点——它试图最小化约束违反度但步长被障碍项卡住。解决先检查约束本身是否自洽sum(Pmax) 是否 负荷 网损再用一个简单技巧给变量一个靠近问题中心的初值而不是用默认的 0。# 给所有机组一个初始出力让功率平衡残差小一些 for i in m.G: m.P[i].set_value((gen_data[i][Pmin] gen_data[i][Pmax]) / 2)这个做法的原理是IPOPT 的恢复阶段受障碍项影响初始点越接近可行域恢复阶段越容易成功。对 ED 这种边界约束居多的模型从每台机组的Pmin Pmax/ 2 出发几乎不会触发 Restoration failed。5.2 求解器报 NaN成本函数或约束中有不可导点现象日志里出现 NaN in ... 或 Invalid number 字样迭代直接中断。原因目标函数或约束表达式里出现了除零、负数开方、或绝对值在零点不可导。ED 里最常见的是阀点效应写成 P_i 的正弦函数再加绝对值项在 P_i 0 附近正弦函数的梯度方向和绝对值项的次梯度冲突IPOPT 的有限差分或自动求导会算出 NaN。解决把非光滑项做光滑化处理。比如用 |x| 近似为 sqrt(x^2 epsilon)epsilon 取 1e-4 量级对于在零处不可导的项改成在零附近用二次函数过渡。另一个检查点是确保 B 矩阵正定如果 B 不是正定矩阵网损表达式可能在某个方向上为负导致约束右侧比实际更小产生看似可行但数值不稳定的迭代路径。5.3 初始点敏感同一种负荷两个初值两个结果现象同一份数据、同一个求解器把初值从各机组平均分配改成全部置为 Pmin得到的最优解不同机组组合和总成本都变了。原因这通常是目标函数非凸造成的。如果模型里加了阀点效应的正弦项或者机组成本曲线在某个区间呈凹形IPOPT 作为局部优化算法只能保证找到初值附近的局部最优。你看到的多解现象是真实的数学性质不是 bug。解决多初始点策略是标准做法。随机生成 10 组初值分别求解取成本最低的解import random best_obj float(inf) best_solution None for _ in range(10): for i in m.G: lo gen_data[i][Pmin] hi gen_data[i][Pmax] m.P[i].set_value(random.uniform(lo, hi)) res solver.solve(m, teeFalse) if res.solver.termination_condition optimal: obj_val pyo.value(m.cost) if obj_val best_obj: best_obj obj_val best_solution [pyo.value(m.P[i]) for i in sorted(m.G)]这套做法在工程上比试图证明全局最优性实在得多。电力调度系统的运行规程也不要求证明全局最优只要求方案可行、成本可比。5.4 跨平台结果不一致macOS、Windows、Linux 上的微妙差异现象在 Windows 笔记本上跑出来的调度方案和 Linux 服务器上跑出来的方案不同机组出力差零点几 MW目标成本差几块钱。原因IPOPT 内部依赖 BLAS/LAPACK 线性代数库不同平台的浮点运算顺序不同导致舍入误差有差异。此外 conda 在不同系统上打包的 IPOPT 可能依赖不同版本的 MUMPS线性求解器的数值行为不完全一致。解决第一把 tol 从 1e-6 改成 1e-8可以有效压缩这种平台间差异——越高的收敛精度越接近数学上的真实最优平台间的微小舍入被迭代过程吸收。第二在项目验收时固定运行环境用同一台 Linux 服务器作为基准所有对比结果都在同一环境上产生。这不算玄学只是浮点运算的本质属性接受它并控制它就是成熟的做法。6. 最后一步热启动、对偶价格与交叉验证经济调度在实际运行中不是孤立的单次求解而是每 15 分钟滚动计算一次。相邻两个时段的数据变化很小上一时段的解就是下一时段绝佳的初始点。用热启动替代冷启动迭代次数通常可以减少一半以上。# 上一时段的解 last_P {1: 72.3, 2: 40.1, 3: 37.8} # 直接把上一时段结果作为本轮初始点 for i in m.G: m.P[i].set_value(last_P[i]) solver.solve(m, teeTrue) # 读功率平衡约束的对偶变量即系统边际成本 m.dual pyo.Suffix(directionpyo.Suffix.IMPORT) lm m.dual[m.balance] print(f系统边际成本 {lm:.4f} 元/MWh)这个对偶变量值得多说两句。IPOPT 在求解过程中已经把功率平衡约束的拉格朗日乘子算出来了它恰好就是经济学上的系统边际成本也就是现货市场出清价格的基础。很多项目额外跑一个 LP 来算这个价格其实没必要——IPOPT 的乘子精度在小规模 ED 上足够直接使用。我早期在项目里犯过的错是只取机组出力不管对偶变量后来做市场结算模块时又不得不把模型重写一遍。建议一开始就把乘子读出来它会让你的系统多一个有用的输出。最后做一次交叉验证。小规模 ED 的解可以用 scipy.optimize.minimize 的 SLSQP 方法对比验证也可以在不等约束整数机组组合要求不高的情况下用网格搜索或蒙特卡洛枚举近似验证最优解的合理性。验证时重点看两组指标功率平衡残差是否达到 1e-6 以下所有机组的出力是否严格在上下界内。这两点通过之后IPOPT 的结果就可以放心进入后续的潮流计算或者下发流程。说到底求解器只是调度系统里的一环它的解必须能在下一步的校验中站得住脚。希望这些经验和踩坑记录能帮到你少走几段弯路。本文还有配套的精品资源点击获取 IPOPT求解电力系统经济调度：从最小模型到参数调优
简介这是一套面向电力系统工程师与科研人员的IPOPT经济调度实现示例围绕电力系统经济调度问题提供完整建模与求解方案。压缩包内含7个MATLAB脚本文件整体仅5KB涵盖目标函数、约束条件、雅可比矩阵、初始点设置及全局变量定义等模块便于理解代码结构与快速复用。项目通过内点法求解机组出力最优分配在满足功率平衡、机组出力限制等约束下最小化总发电成本适合具备一定优化基础和MATLAB使用经验的读者入门实践。已有139人学习下载可作为学习IPOPT接口调用、掌握电力经济调度建模思路的参考范例。1. 电力系统经济调度的 IPOPT 解法先跑通一个最小模型再谈别的负荷预测曲线每小时刷新15 分钟后下一时段的机组出力方案就要落到调度台。这就是电力系统经济调度ED的典型场景给定负荷算出一组机组出力让总发电成本最低同时满足功率平衡和机组自身限制。这个优化问题一旦带上火电机组的非线性成本特性、网损修正和爬坡约束就成了一个地道的非线性规划NLP而 IPOPT 正是开源求解器里处理这类问题最常用的那一个。标题里那一串看起来像仓库文件夹名的字符本质上就是用 IPOPT 求解电力系统经济调度这个最小工程骨架。如果你正在接手调度系统的算法模块、在做求解器选型或者想把手里的商业求解器换掉一部分这篇笔记能让你从模型写法一路走到参数调优和排错。全文代码围绕一个三机系统展开单母线模型几秒钟出结果足够把原理和坑都讲清楚。2. 把 ED 写成 IPOPT 能解的 NLP目标函数、约束与网损模型2.1 为什么选 IPOPT从 ED 的非线性结构看求解器选型经济调度的经典目标函数来自火电机组的煤耗特性曲线。工程上最常用的是二次函数形式F_i(P_i) a_i P_i^2 b_i P_i c_i其中 P_i 是第 i 台机组的出力MWa_i、b_i、c_i 是机组成本系数。如果只保留 b_i P_i 项模型退化成线性规划LP用单纯形法就能解但二次项描述的是机组在高低负荷区间的煤耗率差异忽略它会让高负荷机组的边际成本被低估调度结果在经济性上明显偏差。加入二次项之后问题就从 LP 变成二次规划QP如果再考虑阀点效应在二次项上加一个正弦修正项目标函数直接变成非凸 NLP。这就是选求解器的关键分岔点。Gurobi、CPLEX 这类商业求解器处理凸 QP 很快但遇到非凸目标或者非线性约束比如网损 B 系数公式它们要么不支持要么性能暴跌。而 IPOPT 用的是内点法对 NLP 的支持是原生且完备的目标函数和约束只需要能求一阶、二阶导数凸不凸它都能迭代只是收敛到局部最优还是全局最优的区别。对实际 ED 场景来说机组成本函数大多是凸的即便加阀点项也只是轻度非凸IPOPT 配合多初始点尝试已经足够。在开源阵营里IPOPT 的生态也最成熟。Pyomo 里一行 SolverFactory(ipopt) 就能调起来和 AMPL、GAMS、JuMP 都有接口。如果你在学术论文里看到用 IPOPT 求解 ED那基本是默认动作。2.2 目标函数、功率平衡与机组限制ED 的三层标准结构ED 的标准数学形式可以拆成三层。第一层是目标函数最小化所有机组的发电成本之和min Σ (a_i P_i² b_i P_i c_i)第二层是等式约束也就是功率平衡Σ P_i P_load P_lossP_load 是系统总负荷P_loss 是输电网损。这个约束必须严格满足哪怕相差 0.1 MW在实际运行中都会造成频率偏差。IPOPT 处理这类等式约束的手段是将其纳入拉格朗日函数通过内点障碍项和牛顿迭代逐步逼近可行解最终收敛时约束残差会小于你设定的 tol 参数。第三层是不等式约束。最基本的包括机组出力上下限P_min ≤ P_i ≤ P_max如果做的是多时段调度还要加上爬坡约束-R_i ≤ P_i(t) - P_i(t-1) ≤ R_i以及系统备用约束Σ P_max_i ≥ P_load P_reserve在 IPOPT 内部不等式约束不会显式参与迭代而是通过引入松弛变量和障碍项把一个带不等式约束的 NLP 转成一系列只含等式约束的子问题。这也是为什么 IPOPT 对不等式约束的初始点比较敏感——如果初始点让某个不等式越界太多障碍项会给出很大的对偶值迭代前期会花大量步数把它拉回来。这一点在第 5 章避坑部分会再展开。2.3 网损的 B 系数近似从忽略到精确的迭代路径实际系统里的网损不能直接忽略。常见做法是 B 系数法把网损写成机组出力的二次型P_loss Σᵢ Σⱼ P_i B_ij P_jB 矩阵通过潮流计算离线求得对 3 机系统是一个 3x3 系数矩阵对上百台机组的系统就是上百阶的稠密矩阵。把 P_loss 代入功率平衡约束后约束本身变成二次等式这对 IPOPT 没有任何额外负担——它本来就是非线性求解器。不过在实现上有一个经验问题B 系数是按某个基准运行点线性化得到的如果你直接把它写进约束得到的调度结果在真实潮流下网损会有偏差。我一般会做 2 到 3 次迭代先按 P_loss 0 求解一组出力用这组出力算网损再代回约束重新求解重复两轮后网损偏差通常可以压到 0.5% 以内。这在工程上比一次性解准得多而且每轮 IPOPT 求解只要几秒。2.4 数据组织与量纲写模型前先准备好这几张表ED 建模翻车最常见的原因不是优化模型本身而是数据和量纲混乱。我习惯把机组参数整理成下表这样的结构再开始写代码。机组a (元/MW²·h)b (元/MWh)P_min (MW)P_max (MW)爬坡 (MW/h)G10.01214.0108520G20.01513.5108020G30.01015.057025这个表格里有三个量纲陷阱成本系数 a 的量级是十的负二次方b 是十的一次方而功率是几十 MW。三者相乘后目标函数值在几千的量级而约束的残差在几百的量级IPOPT 内部的自动缩放scaling通常能处理但如果你的机组规模到了上百台某些 a 系数会变成十的负五次方我建议手动做一次缩放把目标函数整体除以 1000或者把功率基准改成 100 MW。量纲问题不解决后面第 5 章说的Restoration failed和迭代停滞就会频繁来找你。3. 用 Pyomo 调用 IPOPT 求解经济调度最小可复现代码3.1 环境准备安装 IPOPT 与 PyomoPyomo 是 Python 生态里最成熟的优化建模框架它支持将模型自动转换为 NLP 格式并调用 IPOPT。安装用 conda 最省事因为 conda-forge 里已经打包好了 IPOPT 的可执行文件不需要自己编译源码。conda create -n ed python3.11 -y conda activate ed conda install -c conda-forge ipopt pyomo -y安装完成后验证一下 IPOPT 是否能正常调用ipopt -v能打印出版本信息就说明求解器就绪。如果你用的是 Windowsconda-forge 也提供了可用的 ipopt 包但要注意从源码编译的 IPOPT 在 Windows 上容易遇到缺少 libgfortran 的问题能装 conda 版就不要自己编译——这算是血泪经验。3.2 最小 ED 模型三个机组、一个功率平衡约束下面的代码是一个完整的最小 ED 模型。三台机组单母线目标函数为二次成本约束为功率平衡和出力上下限。负荷取 150 MW。import pyomo.environ as pyo # 机组数据a 为二次成本系数b 为一次成本系数Pmin/Pmax 为出力上下限 gen_data { 1: {a: 0.012, b: 14.0, Pmin: 10, Pmax: 85}, 2: {a: 0.015, b: 13.5, Pmin: 10, Pmax: 80}, 3: {a: 0.010, b: 15.0, Pmin: 5, Pmax: 70}, } P_load 150 # MW m pyo.ConcreteModel() # 机组集合 m.G pyo.Set(initializegen_data.keys()) # 出力变量直接绑定上下界避免额外写不等式约束 m.P pyo.Var(m.G, boundslambda m, i: (gen_data[i][Pmin], gen_data[i][Pmax])) # 目标函数最小化总发电成本 m.cost pyo.Objective( exprsum(gen_data[i][a] * m.P[i]**2 gen_data[i][b] * m.P[i] for i in m.G) ) # 功率平衡约束总出力 总负荷 m.balance pyo.Constraint( exprsum(m.P[i] for i in m.G) P_load ) # 调用 IPOPT 求解 solver pyo.SolverFactory(ipopt) solver.options[tol] 1e-8 result solver.solve(m, teeTrue) # 打印结果 for i in sorted(m.G): print(f机组 {i}: P {pyo.value(m.P[i]):.2f} MW) print(f总出力 {sum(pyo.value(m.P[i]) for i in m.G):.2f} MW) print(f总成本 {pyo.value(m.cost):.2f} 元/h)这段代码有几个值得说明的设计。第一变量 m.P 的 bounds 参数直接传入了区间Pyomo 会将其生成为变量边界IPOPT 内部用障碍项处理比显式写 Inequality 约束效率更高。第二功率平衡用的是 而且是单条约束IPOPT 对等式约束的处理比不等式更直接。第三tol 设到 1e-8 而不是默认的 1e-6原因后面会说但先记住这个值在 ED 场景下更好用。运行这段代码你会看到 G1 和 G3 承担大部分出力G2 因为成本系数偏高而略少这是符合直觉的结果。3.3 加网损后的模型改动B 系数约束怎么加进去实际系统不能忽略网损在上一节模型的基础上扩展把 B 系数矩阵和网损表达式加进去import numpy as np # B 系数矩阵三机系统的网损二次型系数 # 注意这里的 B 按 100 MVA 基准折算为标幺值实际需根据潮流计算得到 B np.array([ [0.021, 0.001, 0.002], [0.001, 0.017, 0.001], [0.002, 0.001, 0.023], ]) # 定义网损表达式P_loss P^T B P def loss_expr(m): idx sorted(m.G) P [m.P[i] for i in idx] return sum(P[i] * B[i, j] * P[j] for i in range(len(idx)) for j in range(len(idx))) # 修改功率平衡约束总出力 负荷 网损 m.balance pyo.Constraint( exprsum(m.P[i] for i in m.G) P_load loss_expr(m) ) solver pyo.SolverFactory(ipopt) solver.options[tol] 1e-8 solver.solve(m, teeTrue) for i in sorted(m.G): print(f机组 {i}: P {pyo.value(m.P[i]):.2f} MW) print(f网损 {pyo.value(loss_expr(m)):.2f} MW)关键点在于balance 约束的右侧从常数变成了变量 P 的二次型整个约束变成非线性的。IPOPT 的处理方式是把这个二次型纳入约束雅可比矩阵通过牛顿法迭代求解。相比上一节的纯线性约束这会让每轮迭代的线性系统求解成本变高但对三机系统来说时间可以忽略。还有一个工程细节不要手动把网损固定成一个常数再求解那样的话功率平衡可能无法满足。正确做法是让网损表达式完整地参与建模IPOPT 会在迭代中自动协调各机组出力使得总出力 - 负荷 - 网损 0在收敛时成立。4. IPOPT 参数怎么设6 个影响 ED 结果的旋钮与推荐基准4.1 IPOPT 参数地图先认识这几个旋钮IPOPT 的可调参数超过 100 个但做 ED 求解真正需要关心的就 6 个。下面这张表是实际项目里常用的配置基准。参数名作用ED 场景推荐值tol收敛容忍度KKT 误差小于该值判定收敛1e-6 到 1e-8acceptable_tol可接受的低精度收敛阈值1e-5acceptable_iter连续多少次迭代满足 acceptable_tol 后提前终止5max_iter最大迭代次数5000linear_solver求解 KKT 线性系统所用的线性代数库mumps / ma57 / ma97mu_strategy障碍参数 μ 的更新策略adaptivebound_push初始点往可行域内部推的距离1e-2print_level控制台输出详细程度5调试时这些参数里对 ED 结果影响最大的是 tol、linear_solver 和 mu_strategy。tol 直接决定功率平衡约束最终残差设 1e-6 时150 MW 负荷可能残留 0.00015 MW 的不平衡量设 1e-8 时残差降到 0.0000015 MW。对调度系统来说这个量级差异看似微乎其微但如果你把出力结果送到潮流计算里做验证1e-6 的残差有时会造成节点电压的可见偏差。4.2 按 ED 场景调参精度优先还是速度优先调参的推荐顺序是先定 tol再定 linear_solver最后调 mu_strategy。下面是实际项目里的一套起点配置opt pyo.SolverFactory(ipopt) # 精度ED 对功率平衡残差敏感tol 偏严 opt.options[tol] 1e-8 # 可接受解如果在迭代中途连续 5 次达到 1e-5也算收敛 opt.options[acceptable_tol] 1e-6 opt.options[acceptable_iter] 5 # 线性求解器mumps 开源且稳健追求性能可换 ma57/ma97需另外安装 opt.options[linear_solver] mumps # 障碍参数adaptive 在目标函数比较光滑时收敛更快 opt.options[mu_strategy] adaptive # 迭代上限防止个别病态案例无限迭代 opt.options[max_iter] 5000 opt.options[print_level] 5 res opt.solve(m, teeTrue)这里解释三个关键选择。linear_solver 选 mumps 是最省心的默认项因为它是 IPOPT 自带的开源稀疏线性直接求解器不需要额外的许可证如果你的系统到了几百台机组、约束几万个mumps 的内存占用和求解时间会明显上升这时候换成 ma57 或 ma97HSL 库的求解器能快 2 到 3 倍但需要单独向 HSL 申请学术或商业授权。tol 从 1e-6 收紧到 1e-8 通常会让迭代次数增加 10% 到 30%但对小规模 ED 也就多几十毫秒值得。还要重点提一下不可接受但对结果有微妙影响的参数bound_push。它控制初始点距离变量边界有多远。IPOPT 默认值是 1e-8意味着初始点几乎贴在边界上。如果你的某个机组在解里处于 Pmin 边界初始点贴边界会造成障碍项数值过大前面几十轮迭代都在把变量往可行域内部推。把这个参数调成 1e-3 或 1e-2等于告诉求解器初始点请离边界远一点对 ED 这种大量机组运行在边界附近的场景有奇效。如果遇到迭代不收敛我还习惯把 print_level 调到 5 查看每次迭代的 inf_pr原始残差和 inf_du对偶残差。inf_pr 不降说明是可行性问题inf_du 不降说明是目标函数梯度或 Hessian 的问题。这个判断方向在排错时非常关键。5. 经济调度求解避坑指南IPOPT 常见报错与排查记录5.1 Restoration failed不可行问题还是初始点问题现象求解器运行几十轮迭代后突然报错 Restoration failed!退出码显示 Infeasible_Problem_Detected结果文件是空的。原因有两种可能。第一种是约束真的不可行比如你把负荷设成 200 MW但机组 Pmax 之和只有 170 MW功率平衡方程永远无法满足。第二种是初始点离可行域太远IPOPT 的恢复阶段restoration phase找不到一个可行点——它试图最小化约束违反度但步长被障碍项卡住。解决先检查约束本身是否自洽sum(Pmax) 是否 负荷 网损再用一个简单技巧给变量一个靠近问题中心的初值而不是用默认的 0。# 给所有机组一个初始出力让功率平衡残差小一些 for i in m.G: m.P[i].set_value((gen_data[i][Pmin] gen_data[i][Pmax]) / 2)这个做法的原理是IPOPT 的恢复阶段受障碍项影响初始点越接近可行域恢复阶段越容易成功。对 ED 这种边界约束居多的模型从每台机组的Pmin Pmax/ 2 出发几乎不会触发 Restoration failed。5.2 求解器报 NaN成本函数或约束中有不可导点现象日志里出现 NaN in ... 或 Invalid number 字样迭代直接中断。原因目标函数或约束表达式里出现了除零、负数开方、或绝对值在零点不可导。ED 里最常见的是阀点效应写成 P_i 的正弦函数再加绝对值项在 P_i 0 附近正弦函数的梯度方向和绝对值项的次梯度冲突IPOPT 的有限差分或自动求导会算出 NaN。解决把非光滑项做光滑化处理。比如用 |x| 近似为 sqrt(x^2 epsilon)epsilon 取 1e-4 量级对于在零处不可导的项改成在零附近用二次函数过渡。另一个检查点是确保 B 矩阵正定如果 B 不是正定矩阵网损表达式可能在某个方向上为负导致约束右侧比实际更小产生看似可行但数值不稳定的迭代路径。5.3 初始点敏感同一种负荷两个初值两个结果现象同一份数据、同一个求解器把初值从各机组平均分配改成全部置为 Pmin得到的最优解不同机组组合和总成本都变了。原因这通常是目标函数非凸造成的。如果模型里加了阀点效应的正弦项或者机组成本曲线在某个区间呈凹形IPOPT 作为局部优化算法只能保证找到初值附近的局部最优。你看到的多解现象是真实的数学性质不是 bug。解决多初始点策略是标准做法。随机生成 10 组初值分别求解取成本最低的解import random best_obj float(inf) best_solution None for _ in range(10): for i in m.G: lo gen_data[i][Pmin] hi gen_data[i][Pmax] m.P[i].set_value(random.uniform(lo, hi)) res solver.solve(m, teeFalse) if res.solver.termination_condition optimal: obj_val pyo.value(m.cost) if obj_val best_obj: best_obj obj_val best_solution [pyo.value(m.P[i]) for i in sorted(m.G)]这套做法在工程上比试图证明全局最优性实在得多。电力调度系统的运行规程也不要求证明全局最优只要求方案可行、成本可比。5.4 跨平台结果不一致macOS、Windows、Linux 上的微妙差异现象在 Windows 笔记本上跑出来的调度方案和 Linux 服务器上跑出来的方案不同机组出力差零点几 MW目标成本差几块钱。原因IPOPT 内部依赖 BLAS/LAPACK 线性代数库不同平台的浮点运算顺序不同导致舍入误差有差异。此外 conda 在不同系统上打包的 IPOPT 可能依赖不同版本的 MUMPS线性求解器的数值行为不完全一致。解决第一把 tol 从 1e-6 改成 1e-8可以有效压缩这种平台间差异——越高的收敛精度越接近数学上的真实最优平台间的微小舍入被迭代过程吸收。第二在项目验收时固定运行环境用同一台 Linux 服务器作为基准所有对比结果都在同一环境上产生。这不算玄学只是浮点运算的本质属性接受它并控制它就是成熟的做法。6. 最后一步热启动、对偶价格与交叉验证经济调度在实际运行中不是孤立的单次求解而是每 15 分钟滚动计算一次。相邻两个时段的数据变化很小上一时段的解就是下一时段绝佳的初始点。用热启动替代冷启动迭代次数通常可以减少一半以上。# 上一时段的解 last_P {1: 72.3, 2: 40.1, 3: 37.8} # 直接把上一时段结果作为本轮初始点 for i in m.G: m.P[i].set_value(last_P[i]) solver.solve(m, teeTrue) # 读功率平衡约束的对偶变量即系统边际成本 m.dual pyo.Suffix(directionpyo.Suffix.IMPORT) lm m.dual[m.balance] print(f系统边际成本 {lm:.4f} 元/MWh)这个对偶变量值得多说两句。IPOPT 在求解过程中已经把功率平衡约束的拉格朗日乘子算出来了它恰好就是经济学上的系统边际成本也就是现货市场出清价格的基础。很多项目额外跑一个 LP 来算这个价格其实没必要——IPOPT 的乘子精度在小规模 ED 上足够直接使用。我早期在项目里犯过的错是只取机组出力不管对偶变量后来做市场结算模块时又不得不把模型重写一遍。建议一开始就把乘子读出来它会让你的系统多一个有用的输出。最后做一次交叉验证。小规模 ED 的解可以用 scipy.optimize.minimize 的 SLSQP 方法对比验证也可以在不等约束整数机组组合要求不高的情况下用网格搜索或蒙特卡洛枚举近似验证最优解的合理性。验证时重点看两组指标功率平衡残差是否达到 1e-6 以下所有机组的出力是否严格在上下界内。这两点通过之后IPOPT 的结果就可以放心进入后续的潮流计算或者下发流程。说到底求解器只是调度系统里的一环它的解必须能在下一步的校验中站得住脚。希望这些经验和踩坑记录能帮到你少走几段弯路。本文还有配套的精品资源点击获取