ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

机组组合优化实战:基于YALMIP+CPLEX的MILP建模与热备用分析

机组组合优化实战:基于YALMIP+CPLEX的MILP建模与热备用分析 简介面向电力系统调度与优化领域的研究人员和工程师这份资料以混合整数线性规划MILP为核心系统讲解并实现机组组合优化问题。模型用整数变量表示机组启停状态连续变量表示出力水平通过YALMIP工具箱在MATLAB中完成建模并调用Cplex求解器在负荷需求、旋转备用和电网安全等约束下寻求燃料成本最小化的最优调度方案。压缩包共7个文件以MATLAB源码.m为核心配套Excel结果表格、Visio最优出力图表以及需求说明文档整体大小仅267KB结构清晰便于对照学习0.05与0.2两种热备用水平下的完整求解与结果展示。已有3051人学习该资源。读者可从中获得可运行程序、典型场景计算结果和可视化图表既能掌握MILP建模、YALMIP与Cplex的配合使用也能为实际电力系统机组组合方案设计与教学提供参考。1. 机组组合不是“能开就开”而是 0/1 变量的取舍把几十台火电、气电机组排进 96 个时段表面看是“够用就开不够再补”可一旦机组数量超过十台单纯凭经验排出的方案和最优方案之间燃料成本差距能到几个百分点。机组组合问题的本质是混合整数线性规划每台机组必须先用 0/1 变量决定是否在线再在连续空间里分配出力。jizuzuheyouhua.m是一套直接可跑的 MILP 建模代码配套的热备用 0.05 和 0.2 两种求解结果 Excel 文件完整展示了负荷平衡、热备用、出力上下限和启动成本之间怎么互相牵制。适合正在做电力系统机组组合课程设计或想把经济调度从经验表格升级成优化模型的工程师。读完可以照着代码改参数也能从结果文件里反推每一组约束的工程含义。如果对 YALMIP 不熟建议先运行一遍再回头看模型求解器返回接近 240 个整数变量解的过程比任何公式都直观。2. MILP 建模启停变量、出力区间与热备用约束怎么进 YALMIP2.1 决策变量binvar 控制启停sdpvar 控制出力机组组合的难点不在目标函数而在决策变量的维度和语义必须与实际调度动作一一对应。用z(i,t)表示第i台机组在第t个时段是否在线0 代表停机1 代表运行用p(i,t)表示对应时段的出力。如果不给变量加整数约束纯线性规划会把每台机组都安排成“开一半”比如 50% 出力这在物理上根本不存在要么并入电网要么解列。所以机组状态的 0/1 化是混合整数线性规划在电力系统里最典型的用途。YALMIP 中定义这两类变量各自只需要一行N 10; % 机组数量 T 24; % 时段数 z binvar(N, T, full); % 启停状态0/1 变量 p sdpvar(N, T, full); % 出力水平连续变量binvar生成二进制变量矩阵sdpvar生成连续变量矩阵。full说明变量不是按稀疏结构存储而是完整矩阵当 N 到 50、T 到 96 时这样定义能减少 YALMIP 内部索引换算的开销。很多课程设计模板会漏掉full在小规模测试时没影响遇到大规模算例就会明显变慢。2.2 目标函数与约束条件一次把模型写完整目标函数一般取最小化燃料成本和启动成本之和。燃料成本与出力近似线性时直接用系数乘出力如果燃料曲线是二次的MILP 里更常见的做法是分段线性化把出力区间切成几段每段对应一个连续变量和一个 0/1 激活变量。下面的代码使用线性系数c并把启动成本放到辅助变量su上su sdpvar(N, T, full); % 启动辅助变量 startup_cost [300;150;80;60;40;30;20;20;15;15]; obj sum(sum(c .* p)) sum(sum(startup_cost .* su)); F []; F [F, sum(p, 1) Demand]; % 负荷平衡 F [F, sum(z .* Pmax, 1) - sum(p, 1) reserve]; % 热备用 F [F, Pmin .* z p Pmax .* z]; % 出力上下限 F [F, su(:, 1) z(:, 1)]; F [F, su(:, 2:end) z(:, 2:end) - z(:, 1:end-1)];这段代码里obj是 YALMIP 的目标函数对象F是约束集合。负荷平衡要求每个时段总出力等于需求热备用约束要求在线机组的可用容量之和减去实际出力后仍然大于等于预定的备用容量。Pmin .* z p Pmax .* z把出力变量限制在机组运行范围内z0时强制p0。启动逻辑用相邻时段状态变化来触发只有当前时段在线而上一时段停机的点才会产生启动成本。约束类型YALMIP 写法物理含义负荷平衡sum(p,1) Demand所有在线机组出力之和等于系统负荷热备用sum(z.*Pmax,1) - sum(p,1) reserve可用容量与当前出力之差必须覆盖热备用出力区间Pmin.*z p Pmax.*z停机出力为 0运行出力不能越限启动触发su z_t - z_{t-1}停机转运行的那个时段计入启动成本reserve的取值可以按最大单机容量或负荷比例设定热备用 0.05 和 0.2 就是两组不同的负荷比例。这个参数越高系统里要留的空载容量越多机组提前启动的时段也会提前整体成本随之上升。2.3 爬坡约束和最小运行时间约束不能省只写上一节的四组约束MILP 就能解出经济性不错的组合但它不满足实际机组的调节能力。烟煤机组从 30% 负荷爬到 80% 需要几十分钟出力变化有爬坡速率上限。爬坡约束用相邻时段出力差来写- RD_i p(i,t) - p(i,t-1) RU_iYALMIP 中直接作用在矩阵切片上Rup repmat(50, N, T-1); % 上升爬坡速率单位 MW/h Rdown repmat(50, N, T-1); F [F, -Rdown p(:, 2:end) - p(:, 1:end-1) Rup];这里把时段之间的出力变化限制在上下爬坡速率之间。另一个容易被漏掉的是最小运行时间约束一台大型联合循环机组启动后至少要稳定运行 4 小时以上不能为了省半小时成本就停机。这类约束可以用累计运行状态或者启动/停机前后事件来线性化YALMIP 里用for循环按机组展开是实际项目中最常见的处理方式。2.4 模型规模估算整数变量数量决定求解难度机组组合的求解时间和整数变量的个数不是线性关系。假设 10 台机组24 个时段仅启停状态z就有 240 个 0/1 变量如果时段数拉长到 96整数变量接近 1000 个CPLEX 的分支定界搜索空间会呈指数级膨胀。这也是为什么很多课程设计的算例只保留 6 到 10 台机组每天 24 个时段就是为了让求解器能在几分钟内给出全局最优解。遇到更大规模问题时通常要先做机组聚合或备用约束聚合把整数变量数量压下去再交给 MILP 求解器处理。3. YALMIP 与 CPLEX 求解从 optimize 到 gap 控制3.1 optimize 调用与结果读取模型定义完成后YALMIP 只是把目标函数和约束打包真正求解工作由 CPLEX 完成。optimize会先做预求解再用分支定界和割平面方法搜索整数解。代码骨架如下ops sdpsettings(solver, cplex, verbose, 2, ... cplex.mip.tolerances.mipgap, 1e-4); sol optimize(F, obj, ops); if sol.problem 0 P_opt value(p); Z_opt value(z); else disp(sol.info); endsdpsettings是 YALMIP 中唯一的求解器配置入口。verbose设为 2 会输出分支定界日志包括当前上下界和节点数日志出现长时间不动的上下界说明求解进入了瓶颈。sol.problem是求解状态码0 表示成功非 0 时不要急着改参数先看sol.info它往往比目标值更有诊断价值。sol.problem 0 表示可行且已返回最优或满足 gap 的解 sol.problem 1 表示求解时间达到限制返回当前最好可行解 sol.problem 2 表示模型不可行这里problem1在实际调度里很常见。时间限制到了但还有可行解调度系统可以先用当前解执行同时启动下一轮滚动优化而不是把程序卡死。把判断分支写进自动化脚本比只检查problem0要稳健。3.2 CPLEX 参数表哪些参数值得手动调CPLEX 参数很多但处理机组组合问题时值得手动调的只有那几个。下面的表是我在多个工程项目里的起点值YALMIP 中的设置项作用常用起点cplex.mip.tolerances.mipgap相对 MIP 间隙1e-4cplex.mip.limits.timelimit最大求解时间600秒cplex.mip.limits.nodes分支节点上限默认cplex.mip.display求解日志输出级别2mipgap是收敛条件表示当前可行解与最优目标值之间允许的最大相对差距。做课程设计时为了证明“全局最优”可以改成1e-6现场运行建议5e-3因为 0.5% 的成本差在燃料价格波动面前可以忽略但求解时间能减少一个量级。timelimit控制最长求解时间机组组合通常是 24 小时切片滚动优化每个窗口留给求解器的逻辑时间不宜超过 10 到 15 分钟。nodes上限适合放在批处理里防止内存爆掉但在单算例调试时保持默认即可。3.3 先跑 LP 松弛再跑 MILP可以快速判断模型是否合理拿到新模型时可以先做一次 LP 松弛求解把所有binvar临时改成连续变量看目标函数下界。这一步不费时间却很容易发现约束中的符号错误。如果 LP 松弛最优值异常低说明约束太松如果 LP 松弛直接不可行那 MILP 必然不可行。YALMIP 里最简单的方式是重新用一个sdpvar变量替换z求解后再把结果和目标函数值对比。常见做法是保留两套变量定义调试时切换不修改约束代码。z_lp sdpvar(N, T, full); % 先当连续变量 p_lp sdpvar(N, T, full); ops_lp sdpsettings(solver, cplex, verbose, 0); olp optimize(F, obj, ops_lp); % 实际使用时需将 z 替换成 z_lpLP 松弛给出的目标值可以作为 MILP 的下界参考一旦分支定界出来的目标值低于这个下界说明模型写错或数据单位不一致。3.4 结果导出与 Excel 文件对照求解完成后value(p)和value(z)是两个应该最先看的变量。把它们按时段整理成表格就是资源包里热备用0.05状态下的机组组合问题求解结果.xls这类文件的来源。导出代码建议保留原始列名out table((1:T), Demand, value(p(1,:)), value(z(1,:)), ... VariableNames, {时段, 负荷, 机组1出力, 机组1状态}); writetable(out, 机组组合结果.xls);读取旧版 xls 文件时要注意 sheet 名中的中文问题。readtable默认读第一个 sheet指定 sheet 时用Sheet, 2更稳妥。先运行head确认列名再做列选择比直接按位置索引要安全。4. 热备用 0.05 与 0.2 结果对比从 xls 里读出的调度逻辑4.1 用 readtable 读取两个场景结果资源包里有两个结果文件文件名对应两个热备用比例。直接用 MATLAB 读取避免手工比对。中文文件名用单引号括起来即可。读取后先看表头和尺寸r05 readtable(热备用0.05状态下的机组组合问题求解结果.xls); r20 readtable(热备用0.2状态下的机组组合问题求解结果.xls); height(r05) height(r20)如果两个文件高度一致说明时段数量相同可以按行对比。列名如果带单位后缀比如“机组1出力(MW)”在后续计算里要通过VariableNames重命名否则用中文索引容易出错。更稳妥的做法是读取前先打印r05.Properties.VariableNames手动确认后把列名映射成自己习惯的英文变量名。4.2 启停状态对比在线机组数与启停切换次数把启停状态按时段做行求和可以得到每个时段的在线机组数state05 r05{:, 4:end}; % 假设第1列是时段第2列是负荷第3列后是机组状态 state20 r20{:, 4:end}; numOn05 sum(state05, 2); numOn20 sum(state20, 2); bar([numOn05, numOn20]);这段代码用在线机组数量做快速对比。备用比例低时在线机组数在低谷时段可能只有 3 台备用比例升高后低谷时段可能仍要保留 5 台在线。另一个更精细的指标是启停切换次数对state矩阵做列方向差分非零位置就是机组状态变化点。切换越频繁说明模型越缺少最小运行时间约束也说明热备用提高后系统更倾向于让机组在低负载区间运行而不是完全停机。对比项热备用 0.05 的一般表现热备用 0.2 的一般表现在线机组数量低谷时段更少低谷时段明显增加单机平均负载率偏高偏低启停切换总次数可能更频繁相对平缓总燃料成本较低较高高峰时段可调容量紧张充裕这张表描述的是常见规律不是结果文件的绝对数值。拿到具体 xls 后应该用自己的数据和表里的方向做对比如果趋势完全相反优先检查模型里的热备用约束是否乘错了z.*Pmax。4.3 vsdx 图里隐藏的约束一致性信息两个.vsdx文件是 Visio 格式的出力图画的是热备用 0.05 和 0.2 下各台机组的最优出力曲线。这类图比 Excel 数字更直观适合做约束反查。看到出力曲线在某个时段出现阶梯跳变时先核对是否撞上爬坡速率限制看到状态在相邻时段反复跳变时检查最小启停约束是否被注释掉看到某台机组长期卡在Pmin附近要确认是不是热备用约束把它硬留在线上。每一类图形特征背后都对应一组约束能看出模型边界在哪里。4.4 把机组组合结果接到潮流计算里做校核机组组合给出的是各时段机组启停与出力真正下发前还要验证线路和断面不过载。常见做法是把优化结果写入电力系统潮流计算 MATLAB 脚本用runpf或自己写的牛顿-拉夫逊迭代做交流潮流校验。热备用 0.2 场景下更多机组在线潮流分布通常更分散但局部线路可能因为机组出力重分配而出现反向潮流。如果潮流脚本里发现线路过载不需要重新求解整个 MILP而是手动调低该区域机组的出力上限再加一轮迭代。这个流程放在课程设计里能从“能算”提升到“能解释”的层次。5. 求解完成后先看原始残差而不是先看目标值MILP 求解器返回最优解时很多人习惯先看总成本但调试模型时一定先执行约束残差检查。YALMIP 提供了一个被低估的函数check它返回约束集合中每条约束的残差在旧版模板里这个功能由checkset承担。先更新变量值再检查约束满足程度比直接调参数更快暴露问题assign(z, value(z)); assign(p, value(p)); assign(su, value(su)); res check(F); [max_res, idx] max(res); if max_res 1e-6 fprintf(最大残差 %f出现在第 %d 条约束\n, max_res, idx); endassign把求解值写回到 YALMIP 变量check才会基于这些数值计算原始残差。残差非正表示满足正数越大说明这条约束在求解结果里被违反。最常见的违反来源有两个一是负荷平衡写成结果总和高于需求二是爬坡约束的矩阵维度写错导致某些时段根本没有被约束。idx可以指向具体约束但不是所有约束都自带可读标签这时需要在构建F时给每条约束命名比如F [F, sum(p,1) Demand : power_balance]让定位从“第几条”变成直接对应“负荷平衡”。检查残差通过后再看 CPLEX 日志里的 gap 和节点数。如果 gap 长时间不降问题通常出在整数变量过多或备用约束把可行域切得太碎可以试试把热备用约束先放宽 10%看目标值是否显著下降。显著下降说明备用成本过高本次热备用场景的经济代价比预想大不显著下降说明模型冗余约束多可以保留原参数。把check结果和 CPLEX 日志放在一起判断才能区分“模型写错”和“求解器没算够”这也是 MILP 类优化问题最值得建立的调试顺序。本文还有配套的精品资源点击获取
返回列表