电力系统机组组合优化调度:Matlab+Yalmip建模与求解实战 简介面向电力系统调度研究者的电子资料以IEEE14节点、30节点、118节点标准测试系统为对象系统阐述机组组合优化调度中的机组启停、功率分配、成本最小化与多约束处理等核心问题适用于使用Matlab进行算法验证和科研仿真的工程师及高校师生。资源包为单个PDF文档约610KB内容涵盖三类节点系统的规模特点如14节点含5台发电机、30节点含6台发电机、118节点含54台发电机、调度方案运行结果展示以及部分Matlab代码示例便于对照学习与复现。目前已有35人浏览学习。文档提炼了不同网络规模下优化目标与约束条件的建模思路展示了从简单到复杂电网的求解效果能帮助读者理解算法适应性与计算要求并为搭建自己的机组组合仿真框架、扩展经济调度等后续研究提供参考。1. 电力系统机组组合优化调度不是“多算一次经济调度”而是一道先选人再上场的双层决策题电力系统机组组合优化调度UC为什么难难就难在它要同时回答两个问题明天哪些机组开机以及开着的机组每台发多少电。前者是离散决策后者是连续分配两者互相咬合。假如某电网明天最大负荷要涨到某一水平你光靠经济调度把开着的机组出力调来调去是不行的——机组有最小技术出力有爬坡极限还有从冷态到并网需要的时间。我最早在某公司做日内调度时真正让我熬夜的从来不是潮流计算而是开机方案风功率预测一上来火电开多了要承担启停成本开少了晚上高峰又顶不住。这篇文章就把这个问题的完整解法拆开从模型到Matlab代码再到我在IEEE14、30、118节点三个标准系统上反复调试换来的参数经验全都梳理给你。适合正在做电力市场出清、新能源并网运行、或准备用UC作为算法验证题目的从业者和研究者。2. 机组组合数学模型怎么建目标函数、约束与求解器选型2.1 UC建模的核心目标函数拆成燃料成本与启动成本先明确一个观点不要一上来就写约束先想清楚目标函数。UC的目标函数由两部分构成——运行燃料成本和开机启动成本。燃料成本通常是机组出力的二次函数启动成本则按机组从停机到开机的一次性费用计。有些实现还会加上空载费用但最小化模型里空载费用已经隐含在燃料成本的常数项里不需要单独列。数学化之后目标函数是燃料成本sum( a_i * p_i(t)^2 b_i * p_i(t) c_i * u_i(t) )启动成本sum( S_i * v_i(t) )其中 v_i(t) 是开机信号当机组在 t 时刻由 0 变为 1 时取 1否则取 0这里的 u_i(t) 是 0-1 启停变量p_i(t) 是机组 i 在时段 t 的出力。整个问题就变成在所有可行的开机组合与出力组合里挑一组让总成本最小。你和纯经济调度的区别就在这一层经济调度里 u_i(t) 是固定已知的参数而 UC 里它是决策变量。这个“自由度”正是求解难度从连续优化跳到混合整数规划的原因。2.2 四类必须写对的约束负荷平衡、备用、容量与爬坡、最小启停时间UC 的约束可以归为四类写代码前先在纸面上列清楚后面排查会遇到大量“模型看着对求解器报无解”的案例根子多半在这四类的某一条里。第一类是系统功率平衡。每个时段所有开机机组的出力之和要等于该时段的系统负荷减去新能源预测出力。这个约束写成等式sum(p_i(t)) D(t) - W(t)。第二类是旋转备用约束。它保障系统在最大出力机组跳闸或负荷突增时的安全。要注意备用约束的正确写法是“可用容量之和”不是“当前出力之和”sum( Pmax_i * u_i(t) ) D(t) R(t)。很多初学者的代码错在写成 sum(p_i(t)) D(t) R(t)一旦负荷低、机组已压到最小出力这个约束会让模型强行多开机成本暴涨。第三类是机组容量与爬坡约束。容量约束是 Pmin_i * u_i(t) p_i(t) Pmax_i * u_i(t)这里用 u_i(t) 把停运机组的出力锁死为 0。爬坡约束绑定相邻时段p_i(t) - p_i(t-1) RU_ip_i(t-1) - p_i(t) RD_i。爬坡量纲问题会在后文单独说这是最经典的坑。第四类是最小启停时间约束。机组一旦从停机转为开机必须连续运行至少 MinUp 个小时反过来从开机转为停机也必须维持停机至少 MinDown 个小时。这类约束用线性不等式累加形式表达Matlab 里用循环生成比向量化更好读也更容易调。在写代码之前我建议你先用纸笔把 14 节点系统的一个单时段 UC 手动解一遍。你会发现开机组合一变最优出力分配整个就变了把这两个问题拆开解再用循环迭代就是很多启发式算法的原型。2.3 求解器选型为什么我不建议从纯启发式开始工程的现实是程序写完之后要能交代最优性。我见过不少同行一上来写遗传算法或粒子群调了两周参数结果某天负荷曲线稍一变化算法给出的解比 MILP 求解器贵了几个百分点还没有办法解释gap。除非系统规模大到上千台机组或者你需要处理非线性很强的约束否则混合整数线性规划是更可靠的选择。Matlab 链路里有三个选型方向。最省事的是直接用 intlinprog 写 0-1 整数规划优点是零依赖缺点是大规模问题求解慢最常用的是 Yalmip Gurobi 或 Cplex 的组合建模用 Yalmip 描述约束求解交给商业求解器这也是目前大部分论文和工业原型的标准配置如果你预算有限可以 Yalmip 开源的 SCIP 或 HiGHS。我一般在 30 节点以下用 intlinprog 做快速验证118 节点以上果断转 Yalmip Gurobi并把 MIP gap 放到 0.1% 到 1% 之间既能保证可解释性又不至于让求解器为了最后一个 0.01% 的收益空转几小时。这里顺便讲一个判断依据如果目标函数里只含一次项、约束都是线性的建模为 MILP如果燃料成本保留二次项Yalmip 会自动把它转成二阶锥或交给 MIQP 求解器118 节点下这种模型会让求解时间成倍增加。所以工程实现里经常把二次成本曲线做分段线性化用若干段直线逼近原曲线把问题拉回 MILP 框架。这条“把问题拉回线性”的思路在第三节代码里会落地。3. 用Matlab把UC跑通从数据准备到Yalmip求解的最小闭环3.1 数据准备机组参数表的列顺序就是你的调试点UC 输入不需要节点导纳矩阵只需要机组参数和负荷曲线。最常见的做法是准备一个文本文件每行代表一台机组列依次为编号、最小出力 Pmin、最大出力 Pmax、爬坡上限 RU、爬坡下限 RD、初始启停状态、燃料成本二次项 a、一次项 b、常数项 c、启动成本 S。下面这个代码片段演示读取% gen14.txt 每行一台机组列顺序见注释 % 列: 编号 Pmin Pmax RU RD u0 a b c S gen load(gen14.txt); nG size(gen, 1); % 机组台数 T 24; % 调度时段数按 1 小时步长 Pmin gen(:,2); Pmax gen(:,3); RU gen(:,4); RD gen(:,5); u0 gen(:,6); % 第 0 时段状态用于最小启停约束的初始衔接 a gen(:,7); b gen(:,8); c gen(:,9); S gen(:,10);这里有一个关键设计把 u0 单独拎出来当作初始状态而不把第一时刻当作上一时刻。从物理角度看机组在调度期开始之前是开着还是停着会直接影响第一个时段能否开机或停机忘记接 u0 会让最小启停约束在边界处失效。负荷曲线建议也放在单独文件里不要为了省事把负荷硬编码进脚本。因为后面你要做 24 点、48 点、96 点的切换试验负荷文件改成可变长度读取能让你少改很多代码。3.2 定义决策变量Yalmip里0-1变量与连续变量分开声明Yalmip 建模的代码风格很固定。启停状态用 binvar出力用 sdpvar两者都按 nG 行、T 列的矩阵声明。注意第三个参数 full 是必须的如果省略Yalmip 默认生成对称结构这里的矩阵必须显式声明为非对称的。u binvar(nG, T, full); % 启停状态u(i,t) 1 表示开机 p sdpvar(nG, T, full); % 出力p(i,t) 在 [Pmin, Pmax] 内 % 引入开机信号变量 v用于计算启动成本 v binvar(nG, T, full); % v(i,t) 1 表示在 t 时刻完成启动 Cstr []; % 开机信号的定义: v 至少等于 u(t)-u(t-1) Cstr [Cstr, v u - [zeros(nG,1), u(:,1:end-1)]]; Cstr [Cstr, v 0];说明一下这个开机信号的做法u(:,1:end-1) 是上一时段的启停矩阵两者相减得到差分开机时差分等于 1。把差值用约束限定为 v 的下界配合目标函数中启动成本为正求解器为了最小化成本会自动把多余的 v 压到 0所以不需要额外限制 v 的上界。这样写比用 max 函数优雅得多——max 在 Yalmip 里会引入非光滑算子某些求解器处理起来反而慢。3.3 目标函数与约束组装搭建一个可扩展的约束结构目标函数建议按行拆开写分别累加燃料成本和启动成本。燃料成本默认是二次函数这里先保留二次项第四节的 118 节点案例再讲分段线性化。objective 0; % 燃料成本: sum_t sum_i a(i)*p(i,t)^2 b(i)*p(i,t) c(i)*u(i,t) for t 1:T objective objective sum( a .* p(:,t).^2 b .* p(:,t) c .* u(:,t) ); end % 启动成本: sum_t sum_i S(i) * v(i,t) objective objective sum( sum( S .* v ) );写循环累加的好处是便于后续扩展多时段耦合约束代价是建模速度略慢。14 节点系统完全无所谓118 节点系统如果你嫌慢可以把循环改成矩阵运算但工程上我更偏好清晰可调试而不是为了省几十毫秒牺牲可读性。代码里的 a、b、c 是列向量p(:,t) 也是列向量点乘保证逐机组运算。约束部分照前面分类逐个追加% 1) 出力上下限约束停运机组出力为 0 Cstr [Cstr, Pmin .* u p Pmax .* u]; % 2) 系统功率平衡每个时段成立 Cstr [Cstr, sum(p, 1) load_w(1:T)]; % load_w 是净负荷曲线 % 3) 旋转备用约束可用容量 负荷 备用 Cstr [Cstr, sum(Pmax .* u, 1) load_w(1:T) reserve(1:T)]; % 4) 机组爬坡约束从 t2 到 T for t 2:T Cstr [Cstr, p(:,t) - p(:,t-1) RU]; Cstr [Cstr, p(:,t-1) - p(:,t) RD]; end % 5) 最小启停时间约束这里用最直观的累加形式 for i 1:nG for t 1:T % 最小开机时间若 u(i,t)0 且 u(i,t-1)1则接下来 MinUp-1 个时段必须为 1 if t T - MinUp(i) 1 Cstr [Cstr, u(i,t) - u(i,t-1) u(i,t1) - u(i,t) ... 这个写法暂不推荐最小启停时间的累加写法要严谨推荐下面这种被广泛使用的四段式形式。假设最小开机时间为 Mu最小停机时间为 Md引入上一时段状态 u0 辅助边界for i 1:nG for t 1:T if t 1 % 首时段与初始状态衔接 Cstr [Cstr, u(i,1) - u0(i) max(0, 1 - u0(i))]; % 简化形式 else % 开机后 Mu 个时段内的后效约束 uptime_idx t 1 : min(T, t Mu(i) - 1); if ~isempty(uptime_idx) Cstr [Cstr, u(i,t) - u(i,t-1) u(i, uptime_idx(1))]; Cstr [Cstr, u(i,t) - u(i,t-1) u(i, uptime_idx(end))]; end end end end必须坦白说最小启停时间约束的正确线性化是 UC 建模中最多人卡壳的地方之一上面这段为了控制篇幅做了化简。实践中我建议直接把“u0 之后连续 M 个时段为 0”的约束用 yalmip 的 约束数组生成写成下面这种更直白的形式哪怕多几行也不容易错for i 1:nG for t 1:T % 停机约束如果本时段从开机变为停机则后续 Md(i)-1 个时段必须维持停机 if t1 T for k t1 : min(T, t Md(i) - 1) Cstr [Cstr, u(i,t-1) - u(i,t) 1 - u(i,k)]; end end end end这种双层循环生成的约束数量多但每一行都有清晰的物理含义排查的时候你能直接指出这条约束是哪台机组、哪个时刻、哪类限制。等到问题调通再考虑改成压缩形式。3.4 求解与结果落盘Gurobi的MIP gap要主动放松约束与目标都组装好之后求解配置决定了你的项目是“秒级返回”还是“泡面级等待”。下面是标准求解配置options sdpsettings(solver, gurobi, ... % 可换 cplex 或 intlinprog verbose, 1, ... % 打印求解日志 mipRelGap, 0.001); % 相对最优间隙 0.1% sol optimize(Cstr, objective, options); if sol.problem 0 u_sol value(u); p_sol value(p); total_cost value(objective); disp([UC 求解成功总成本 , num2str(total_cost)]); else disp([求解失败Yalmip 错误码: , num2str(sol.problem)]); % 排查: 用 sol.info 看求解器原始返回信息 endmipRelGap 是工程里最值得设置的参数。0.001 意味着求解器找到的解与最优解之间的差距不超过 0.1% 就收手。做研究时有人习惯追求 1e-6 的严格最优但工程上负荷预测本身就有 3% 以上的误差0.1% 的求解精度已经完全够用而计算时间可能相差一个数量级。求解成功之后务必把结果保存到 mat 文件而不是直接关掉工作区。后面做 96 点时段的滚动调度、或者对比不同目标函数参数时你会发现自己最需要的是历史结果存档。这个小习惯帮我避免过很多次“方案A和方案B哪个好”造成的争论——因为文件里存着每次求解的全部变量。4. IEEE14/30/118节点怎么切换规模差异、参数修正与提速策略4.1 三种规模到底差在哪不是节点数而是机组数IEEE 14、30、118 节点系统在 UC 里真正起作用的不是母线数量而是发电机台数和负荷量级。14 节点常见配置只有 3 到 5 台机组跑起来几乎是瞬时的适合验证模型正确性30 节点通常有 6 台左右机组可以开始体现爬坡约束的价值118 节点接入五十余台机组负荷和备用需求也大了一个量级这才是真正检验建模效率和求解器性能的场景。三套系统的差异用一张表列清楚项目IEEE 14节点IEEE 30节点IEEE 118节点常见机组台数3~56 左右54 左右负荷量级百 MW 级百 MW 级数千 MW 级是否建议加线路潮流约束可加可不加建议加必须加建议求解器intlinprog / GurobiGurobi / CplexGurobi / CplexMIP gap0.1%0.1%~0.5%0.5%~1%注意节点数本身不是 UC 问题的维度真正决定规模的是机组数乘以时段数。118 节点系统 54 台机组乘以 24 个时段0-1 变量就有 1296 个连续变量同样数量这已经是中等规模的 MILP 问题。4.2 从14节点到30节点什么时候该加线路约束14 节点系统上如果数据只给了机组参数和负荷曲线没有线路阻抗数据UC 就可以完全忽略网络约束。30 节点开始市面上能拿到的数据包里通常带线路参数这时候我会建议加一套直流潮流约束否则 UC 的结果可能在某个线路上严重越限物理上根本不可行。直流潮流约束的引入方式很标准把支路功率表示成节点相角和电纳的函数再加支路容量上下限。Yalmip 里这需要额外引入相角变量 theta约束写成 P_line B * theta 的形式。但这里有一个工程提速技巧大部分 UC 论文里用的是“传输约束的凸松弛”先把线路容量约束加上跑一遍如果没有任何一条支路达到容量边界说明网络约束不起作用可以从模型里删掉以提速。这个技巧在 30 节点上非常实用因为 6 台机组的 UC 很容易满足所有潮流约束删了之后求解速度几乎不变但代码大幅简化。4.3 上118节点求解时间失控的破解顺序118 节点是检验工程能力的分水岭。我遇到过的典型情况是同一条代码14 节点 2 秒出结果30 节点 10 秒出结果换成 118 节点直接半小时没跑完最后被迫中断。这不是求解器不行而是模型里有三个拖慢速度的元凶。第一个元凶是成本函数保留了二次项。前面提过二次项会让问题变成 MIQP 或需要外近似处理。Gurobi 虽然能直接解 MIQP但求解时间比同规模的 MILP 高一个量级以上。破解做法是把二次成本分段线性化比如把每台机组的出力区间切成 4 到 6 段用分段线性函数近似原曲线。分段数越多精度越高但约束也越多工程上 4 段已经很够用。% 对机组 i 的二次成本做 5 段线性近似示例 p_seg_size (Pmax(i) - Pmin(i)) / 5; seg_boundaries Pmin(i) : p_seg_size : Pmax(i); % 引入分段出力变量 pseg, size 1 x 5 pseg sdpvar(1, 5); Cstr [Cstr, pseg 0, pseg p_seg_size]; Cstr [Cstr, p(i,t) sum(pseg) Pmin(i) * u(i,t)]; % 目标中对应项的系数为各段中点的斜率注意最后那条等式约束里 pseg 只是区间内的增量机组总出力等于 Pmin 乘以开机状态加上各段增量之和。这样目标函数里的二次成本就变成各段斜率的线性加权整个模型退化为纯 MILP。第二个元凶是备用约束被写成全时段强约束。118 节点的高峰负荷时段的备用需求很大如果每个时段都要求全系统可用容量大于等于负荷加固定备用模型会被无谓地撑大。常见做法是把备用需求改成按负荷比例计算的时变值比如 R(t) 5% * D(t) 2% * W(t)这样低谷时段约束强度降低求解更快。第三个元凶是求解精度设置过高。直接把 mipRelGap 从 0.001 放到 0.01即 1% 的相对间隙118 节点的求解时间经常从半小时降到两分钟以内。工程上 1% 的成本差通常远小于负荷预测误差带来的影响代价完全值得。5. 机组组合求解的五个高频翻车现场与排查清单5.1 从24点改到96点后模型突然无解爬坡约束的量纲陷阱现象24 个时段的模型运行正常把负荷数据改成 15 分钟一个点T96之后求解器立刻报 INFEASIBLE。原因爬坡约束的数值是用“每小时内机组能增加多少 MW”定义的比如 RU 50 MW/h。调度时段改成 15 分钟之后每个时段的爬坡能力实际只剩 50 / 4 12.5 MW/时段但代码里 RU 仍然按 50 参与约束导致相邻时段出力跳变被允许得过大而后面的最小启停时间约束又限制机组不能随意启停互相拧巴。如果负荷曲线变化平缓还好一旦早晚高峰负荷陡升陡降必无解。解决把爬坡值统一换算成“每个调度时段的爬坡能力”。正确写法是 RU_时段 RU * (时段长度/60)。如果 T96 就除以 4T48 就除以 2。类似地最小启停时间约束中的 MinUp / MinDown 也要从“小时数”换算成“时段数”。这个换算的最佳位置在数据读取之后用一个 config 变量单独定义调度步长后续所有约束参数都基于它计算。5.2 求解器报无解但约束“看起来都满足”备用约束方向性错误现象模型在峰值负荷时段无解。检查出力上下限、负荷平衡都找不到冲突点最后逐条翻约束才发现旋转备用约束的数值异常偏大。原因代码把备用错误写成 sum(p) D R而不是 sum(Pmax * u) D R。用当前出力之和做备用意味着备用容量被计算成“正在发出的功率”而实际上备用指的是“还能再多发多少”差值等于 (Pmax - p) * u。我用 sum(p) 写错过一次结果是模型为了满足备用约束不得不让所有机组满发然后负荷平衡约束又要求削减出力两个约束互相矛盾无解是必然。解决把备用约束改成 sum(Pmax .* u) D R同时加一条备用容量上限约束 sum((Pmax - p) .* u) R。如果数据里有部分机组不参与备用调节还要在参与备用的机组子集上计算。5.3 118节点求解一小时不收敛二次成本没有分段线性化现象程序没有崩溃也没有报错但 Gurobi 日志里的 MIP gap 一直卡在 20% 左右下不去跑了半个多小时没有任何改善迹象。原因目标函数里的二次项让每个连续变量对应的梯度随出力变化分支定界过程中下界估计非常松散。说白了就是求解器在尝试逼近一个非线性曲面线性松弛给不出紧的下界。解决第 4.3 节的分段线性化代码套用上去同时把 MIP gap 设为 1%。这两步做完同样的 118 节点问题通常在几分钟内收敛。我自己的习惯是先加分段线性化再看求解时间如果仍然太慢才放 gap。5.4 大M系数拍脑袋乱设解出来的出力是“假的”现象某次求解正常结束但检查结果时发现有一台机组在停机状态下出力为 0.8 MW既不在 [Pmin, Pmax] 区间内也明显小于 Pmin。排查很久才发现问题出在一个关联约束里的大M系数。原因为了表示“只有开机状态下才能有出力”这类关联关系代码里用了 p M * u 这个约束其中 M 1e6。由于求解器的数值容差有限当 u 0 时约束 p 1e6 实际上不起限制作用p 可以残留一个极小值。更糟的是如果其他约束里有类似量级的 M 参与数值稳定性会明显变差。解决工程上 M 系数永远取该约束物理上界的 1.1 倍不要为了省事统一用大数。对于 p M * u直接写 p Pmax * uPmax 本身就是出力上界这才是符合物理的大M。同理备用约束里表示可用容量时直接用 Pmax * u 而不是大M。把 M 系数全部换成物理边界之后这类“伪解”问题会消失求解器数值稳定性也会改善。5.5 结果总成本很低但负荷曲线回放对不上忘了网损现象UC 出来的总成本很有竞争力但把出力曲线画出来叠加负荷曲线发现总有那么几个时段出力之和明显小于负荷或者大于负荷不过功率平衡约束又是等式约束理论上不该出现这个偏差。原因代码里用的负荷数据可能是系统总负荷但 UC 没有计入网损。14 节点系统网损只有百分之几察觉不出来118 节点网损可以达到负荷的 3% 到 5%忽略后等式的平衡点在物理上根本不成立。解决第一档用历史数据求一个平均网损系数把负荷乘上 1/(1-η) 再参与平衡约束第二档在直流潮流模型里显式加入网损表达式。工程上先用第一档观察净负荷曲线在所有时段都满足出力上下限再考虑升级到第二档。6. 让结果可信时域回放、随机扩展与验收习惯6.1 让UC结果过“时域回放”这一关求解不报错只是第一步结果可信才是交付的底线。我每次跑完会把 u_sol 和 p_sol 画到一张图上横轴是 24 个时段纵轴是每台机组的出力用不同颜色区分再叠加负荷曲线。人眼扫一遍主要看三点有没有机组频繁启停时间上像锯齿、有没有机组出力长期压在 Pmin 或 Pmax 边界上、负荷曲线和出力总面积是否吻合。这些任何数值指标都替代不了因为它能暴露那些“约束满足但物理诡异”的解。6.2 从确定性UC向随机UC扩展的思路确定性 UC 假设负荷和新能源出力已知实际工程里这两者都是预测值。当你把确定性 UC 跑顺之后下一步是把不确定性引入常见做法之一是场景法把风电的预测误差按正态分布采样出多个场景在每个场景下都要求满足备用约束目标函数改为所有场景成本的期望。做法之二是鲁棒优化用不确定区间描述预测误差保证最坏场景下也能满足负荷和备用。两种做法都是在现有模型上增加一层维度而不是推翻模型所以把确定性 UC 的代码写清晰、把约束结构做模块化是在为后面的扩展省时间。6.3 一个让我少踩很多坑的固定习惯我现在拿到任何新的 UC 数据包第一件事不是直接跑 118而是先裁一个小系统取 6 台机组、24 个时段跑通一遍完整流程再逐步扩规模。每改一个参数比如爬坡上下限或备用比例我都会把结果存一份带日期的版本。某天某导师来问“上周那组结果负载率为什么那么高”时我能在五分钟内调出当时的参数文件复盘而不是对着已经清空的工作区发愣。希望这组模型、代码和调试顺序也能帮你把这个方向上的坑提前填平。本文还有配套的精品资源点击获取