简介面向电力系统调度、运筹优化和Matlab编程学习者该程序以两台发电机组为简化对象完整演示了基于Matlab与CPLEX联合求解机组组合问题UCP的流程。包体共57个文件主要包括45个.mod模型文件对应不同机组容量与场景参数、10个.m脚本覆盖主程序、数据读取、求解调用、结果打印与成本计算等环节、1个.lp问题实例和1个.asv自动保存文件整个压缩包仅274KB轻量便于快速部署。资源已有1199人学习适合作为教学演示或入门实践素材。通过阅读代码可掌握如何定义目标函数如最小化发电与启停成本、设置负荷平衡与机组出力上下限等约束并体会非线性因素与整数变量在模型中的处理方式同时理解Matlab如何通过接口调用CPLEX以及.mod模型文件与.m主程序的分工协作为后续扩展更多机组或复杂约束提供可直接修改的模板。理解这些建模思路对解决实际机组组合问题乃至其他领域的优化问题也有借鉴意义。1. 两个机组就能讲透的东西Matlab调用CPLEX做机组组合调度做发电调度的人第一次碰机组组合Unit Commitment, UC时最常见的翻车方式不是模型不会写而是面对二三十台机组、96个时段、爬坡和启停约束堆在一起的MIP模型连“报错是数据问题还是约束写错”都判断不了。这个标题对应的“2变量机组组合调度程序”其实就是把UC压到最小可运行的骨架2台发电机、24个时段用Matlab作为建模壳调用CPLEX求解混合整数规划跑通“建模型—组装约束—求解—回读结果”的整条链路。它非常适合两类人想找个matlab调用cplex例题入门的人以及需要把UC从论文公式改成能排产的Matlab代码的从业者。2台机组暴露的问题足够替你排掉后面80%的坑。2. 2变量机组组合数学模型决策变量、约束排布与目标函数的取舍2.1 为什么2变量是入门骨架变量规模与物理含义机组组合问题要回答两个问题未来每个时段哪些机组处于开机状态各台机组发多少有功功率。所以“2变量”不是指只有两个标量变量而是指2台机组的启停状态变量加出力变量合起来才是完整的调度模型。以24时段、2台机组为例决策变量总数为2台×24时段×2类 96个其中48个是二进制启停变量48个是连续出力变量。工业现场的UC模型动辄几万台机组几十个约束块你是没法靠调试器一行行看的。把规模压到2台机组的价值在于所有变量的索引可以用眼睛直接验算所有约束的物理含义可以逐一对应CPLEX求解器返回的任何异常都能立刻定位是建模问题还是数值问题。这套最小骨架包含的物理对象是完整的火电机组有最小技术出力、最大出力、燃料成本曲线、启动成本系统有逐时段负荷需求和备用需求。它和实际工程的差异只是约束数量不是问题本质。从这个意义上说2变量不是玩具是“可解释的最小可行模型”。2.2 目标函数怎么定燃料二次成本与启机成本线性化机组组合的目标是让调度期内的总运行成本最小。每台机组在时段t的成本一般写成燃料成本C_i(p) a_i * p^2 b_i * p c_i其中p是出力a、b、c是机组的燃料成本系数这是典型的二次曲线。启动成本SU_i只要机组从停机变为开机就产生一笔固定费用简化处理时不区分热启动和冷启动。于是总目标函数是min sum_t sum_i ( C_i( p_i,t ) SU_i * u_i,t )这里u_i,t是二进制启停变量p_i,t是连续出力变量。用CPLEX的Matlab接口实现时目标函数要拆成线性部分和二次部分线性部分是b_i * p_i,t c_i * u_i,t?。注意一个容易搞错的细节常数项c_i只有机组开才发生所以实际线性目标是b_i * p_i,t (c_i SU_i) * u_i,t而不是把c_i单独放在一边。二次项是a_i * p_i,t^2放进CPLEX目标矩阵Q。有的教材会把二次成本分段线性化那是为了照顾老式线性求解器。CPLEX本身支持MIQP混合整数二次规划2机组规模直接用二次目标结果是精确的代码也更短。大规模工程模型里大量使用分段线性逼近那是因为模型规模到了几十万变量、且需要保持纯LP结构这个我们在最后一章再谈。2.3 约束排布负荷平衡、备用容量与出力上下限约束是机组组合里最值得花时间的部分。2变量骨架里至少要有三类约束每类都对应明确的物理要求负荷平衡约束每个时段所有开机机组的出力之和等于系统负荷。这是等式约束写成sum_i p_i,t D_t它也是最容易把模型搞成不可行的一处比如负荷数据大于2台机组最大出力之和或者错写成了不等式模型立刻无解。备用容量约束系统需要额外备用容量要求当前已开机机组的最大出力之和要大于负荷加备用写成sum_i Pmax_i * u_i,t D_t R_t。注意这里起作用的是开机机组的容量上限不是实际出力。出力上下限约束机组在开机状态下出力必须在[Pmin_i, Pmax_i]内停机状态下出力必然为0。这是典型的“状态与出力联动约束”写成两个不等式p_i,t Pmax_i * u_i,t上限联动p_i,t Pmin_i * u_i,t下限联动这三类约束构成了UC的核心最小启停时间、爬坡速率约束在2变量骨架里先不引入。原因很实际最小启停时间需要引入“上一时段状态”形成跨时段耦合爬坡约束需要时段间差分这些在N台机组版本里再展开第一版跑通时加进去只会让排查复杂化。约束在Matlab里的组织方式决定了后续调试的体验。常见做法是统一构造一个稀疏矩阵A配合lhs和rhs两个边界向量。每行表示一个约束lhs和rhs相等时是等式lhs为负无穷时是小于等于约束rhs为正无穷时大于等于约束。这种统一格式和CPLEX的Matlab接口天然对齐也方便后续用cplex.writeModel()输出LP文件检查。2.4 参数设置机组参数表与24时段负荷曲线2台机组虽然规模小但参数不能太理想化否则约束根本不起作用。我用一组典型参数单位是传统的火电标幺风格功率用MW成本系数用$/MW^2h、$/MWh、$/h。两台机组刻意设计成不对称的一台是节能但启动慢的大机组一台是灵活但成本较高的小机组这样才能看到调度结果是“怎么组合”的。参数机组1机组2Pmin (MW)100150Pmax (MW)300400a ($/MW^2h)0.0200.015b ($/MWh)2017.5c ($/h)120100SU ($/次)500700负荷曲线取24个时段的典型日负荷晚高峰出现在19-20点低谷在凌晨4-5点范围大致在380MW到700MW之间。这样设计有两个用意凌晨低谷时段2台机组能不能“停机一台”由经济性决定晚高峰时段2台机组的最大出力总和700MW刚好覆盖负荷逼迫备用约束产生作用让1号或2号机组必须保持开机。备用系数取10%即R_t 0.1 * D_t。若总负荷高峰700MW加备用70MW而两台机组的启用上限之和只有700MW那么备用约束会约束到必须把大机组开机这正是模型的核心行为。3. 用Matlab调用CPLEX跑通最小算例完整代码与参数细节3.1 环境准备CPLEX版本与Matlab路径匹配先把环境说清楚这一节不处理好后面代码一个都跑不起来。CPLEX提供Matlab工具箱常见安装路径是C:\Program Files\IBM\ILOG\CPLEX_Studio...\cplex\matlab。在Matlab里把这个目录加入路径addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio2211\cplex\matlab); savepath; % 保存到默认路径避免每次启动重新添加代码说明你有自己的安装版本时把路径里的版本号改成实际目录。savepath能把路径存进Matlab的pathdef文件这样下次启动就不用重新addpath。如果不执行savepath当前会话有效重启后失效。验证是否安装成功运行Cplex构造函数看是否报错try cplex Cplex(test); cplex.solve(); disp([CPLEX版本可用退出码: num2str(cplex.Solution.status)]); catch err disp([Cplex类初始化失败: err.message]); end另一个常见的坑是Matlab版本与CPLEX版本不匹配比如较新的Matlab版本不在CPLEX支持列表里Cplex的动态库加载会报“找不到指定模块”。最直接的判断方法是在CPLEX安装目录的cplex\matlab下有.mexw64文件如果缺失对应的Matlab版本文件基本就说明自带库不兼容。常见做法是用CPLEX配套的make.m命令行在本地重新编译mex文件但我遇到的问题是系统缺少相应版本的C运行时库装上之后再编译就通了。如果你只是想快速跑通示例优先在CPLEX官方支持列表里对应的Matlab版本上运行别在这个环节浪费时间。3.2 数据准备生成2台机组与24时段负荷数据构造放在脚本最前面单独成一个分区方便后续替换成真实机组数据。下面这段把机组参数、负荷、备用全部向量化T 24; % 调度时段数取24整点小时 nUnit 2; % 机组台数2变量骨架 % 机组参数Pmin、Pmax、成本二次/一次/常数、启动成本 Pmin [100; 150]; % MW Pmax [300; 400]; % MW a [0.020; 0.015]; % $/MW^2h b [20; 17.5]; % $/MWh c [120; 100]; % $/h SU [500; 700]; % $/启动次数 % 24时段负荷曲线典型日单位MW D [450 420 400 390 380 390 430 510 570 ... 610 640 660 630 610 590 600 630 650 ... 680 700 650 560 500 460]; % 备用容量取负荷的10% R 0.10 * D;代码说明所有参数都是列向量后面组装索引时按(i-1)*T t来取这个索引规律和Matlab的列优先顺序完全一致。负荷曲线D是行向量长度24。备用系数0.1在工程上表示需要10%的旋转备用实际电网中这个值和系统规模、新能源占比强相关可以后续改成逐时段变化的序列。参数是怎么影响行为的2号机组最大出力400MW大于1号机组且成本曲线更平缓a0.015b17.5所以晚高峰时2号机组倾向于承担基荷和峰荷1号机组负责低谷和补充。如果想让模型出现“低谷停一台”的行为需要确保负荷低谷时两台机组的最小技术出力之和大于负荷需求这里最小出力总和是250MW低谷负荷390MW所以模型大概率会保持两台都在线只是让其中一台贴近最小出力。想看到“停机”行为可以手动把低谷负荷改到240以下。这个敏感性值得你亲手试。3.3 组装模型并求解Cplex对象的标准流程核心部分来了。按照CPLEX的Matlab接口变量顺序约定为先排所有u变量再排所有p变量在每个类型内部按机组编号外循环、时段内循环。即u_i,t的索引(i-1)*T tp_i,t的索引nU (i-1)*T t其中nU 2 * T完整组装和求解代码如下nU nUnit * T; % 启停变量个数 nP nUnit * T; % 出力变量个数 nVar nU nP; % 总变量数 % 线性目标系数 f f zeros(nVar, 1); for i 1:nUnit for t 1:T ui (i-1)*T t; pi nU (i-1)*T t; f(ui) c(i) SU(i); % 开机线性成本常数成本启动成本 f(pi) b(i); % 一次项燃料成本 end end % 变量边界 lb zeros(nVar, 1); ub ones(nVar, 1); % 先默认所有变量都在[0,1]p的ub后面按Pmax修正 for i 1:nUnit for t 1:T pi nU (i-1)*T t; lb(pi) 0; % 出力下限由Pmin约束保证这里不写死 ub(pi) Pmax(i); % 出力上限 end end % 变量类型前nU个二进制后nP个连续 ctype [repmat(B, 1, nU), repmat(C, 1, nP)]; % 二次目标矩阵 QCPLEX目标为 0.5*x*Q*x f*x Q sparse(nVar, nVar); for i 1:nUnit for t 1:T pi nU (i-1)*T t; Q(pi, pi) 2 * a(i); % 因为目标里是 a*p^2CPLEX要乘0.5 end end % 约束矩阵统一组装sparse矩阵A lhs/rhs边界 nConstr T T 2*nUnit*T; A sparse(nConstr, nVar); lhs -inf * ones(nConstr, 1); rhs zeros(nConstr, 1); row 0; % 第1类负荷平衡 sum_i p_i,t D_t for t 1:T row row 1; for i 1:nUnit pi nU (i-1)*T t; A(row, pi) 1; end lhs(row) D(t); rhs(row) D(t); end % 第2类备用约束 sum_i Pmax_i*u_i,t D_t R_t for t 1:T row row 1; for i 1:nUnit ui (i-1)*T t; A(row, ui) -Pmax(i); % 挪到小于等于侧 end lhs(row) -inf; rhs(row) -(D(t) R(t)); end % 第3类p_i,t Pmax_i * u_i,t for i 1:nUnit for t 1:T row row 1; ui (i-1)*T t; pi nU (i-1)*T t; A(row, ui) -Pmax(i); A(row, pi) 1; lhs(row) -inf; rhs(row) 0; end end % 第4类p_i,t Pmin_i * u_i,t for i 1:nUnit for t 1:T row row 1; ui (i-1)*T t; pi nU (i-1)*T t; A(row, ui) Pmin(i); A(row, pi) -1; lhs(row) -inf; rhs(row) 0; end end % 创建Cplex对象并求解 cplex Cplex(uc_2units); cplex.Model.sense minimize; cplex.Model.obj f; cplex.Model.Q Q; cplex.Model.lb lb; cplex.Model.ub ub; cplex.Model.ctype ctype; cplex.Model.A A; cplex.Model.lhs lhs; cplex.Model.rhs rhs; % 求解并检查状态 cplex.solve(); if strcmp(cplex.Solution.status, optimal) x cplex.Solution.x; fprintf(最优成本: %.2f\n, cplex.Solution.objval); else disp([求解失败状态: cplex.Solution.status]); end代码说明这套组装方式把“数学约束”和“数据约束”统一成lhs A*x rhs格式写起来冗余但排查价值很高。备用约束用负号把转成是为了不对lhs单独处理正负无穷如果你愿意也可以用lhs为D(t)R(t)、rhs为正无穷的正向写法我一般习惯统一到“大数在右、负无穷在左”这样LP文件里看起来更规律。参数说明ctype里B和C必须严格从前到后对应变量顺序一旦错位CPLEX不会报错但会把整数变量当成连续变量结果完全失真。Q矩阵对角线的2倍系数是这个环节最经典的一个坑我们在避坑章节专门讲。cplex.Model.A用稀疏矩阵存储2变量规模看不出差异但N台机组时稀疏和稠密的内存差距能差几十倍。3.4 结果回读与校验开机组合和出力分配怎么看求解成功后把变量拆回u和p并且做一次“事后校验”这是我在跑任何调度程序时的固定习惯。不要只盯着objval要重新验证约束数学上是否成立x cplex.Solution.x; u reshape(x(1:nU), T, nUnit); % 按机组分行时段分列 p reshape(x(nU1:end), T, nUnit); % 校验负荷平衡 load_check sum(p, 2); % 每时段出力求和 balance_err max(abs(load_check - D)); fprintf(负荷平衡最大误差: %.4f MW\n, balance_err); % 校验备用约束 reserve_check sum(Pmax .* u, 1); % 每时段开机容量上限 violation reserve_check - (D R); fprintf(备用约束最小裕度: %.2f MW\n, min(violation)); % 打印机组1的开机状态 fprintf(机组1开机时段: ); fprintf(%d , find(u(1,:) 0.5)); fprintf(\n);代码说明reshape必须是(T, nUnit)而不是(nUnit, T)原因是变量索引按列优先组装时第i台机组的所有时段连续排列。方向写反输出结果会错位这个我踩过。find(u(1,:) 0.5)处理二进制变量浮点误差造成的0.9999比直接判1稳。输出来自2机组模型的经济调度含义是成本低的机组优先承担基荷成本高的机组负责高峰和备用容量既看“谁开着”也看“发多少”。4. Matlab调用CPLEX避坑指南5个高频问题与排查路径4.1 Cplex类找不到版本不匹配和路径遗漏现象Cplex(test)直接报“未定义变量Cplex”或“找不到指定的模块”。原因就三类路径没加对、Matlab启动时没有保存路径、CPLEX的Matlab接口版本和当前Matlab不兼容。解决% 第一优先确认实际安装目录 dir(C:\Program Files\IBM\ILOG\CPLEX_Studio*\cplex\matlab) % 如果能看到cplex.m、Cplex.m等文件则路径无误 addpath(genpath(C:\Program Files\IBM\ILOG\CPLEX_Studio2211\cplex\matlab));如果路径正确却报“找不到指定的模块”大概率是Matlab位数与CPLEX位数不一致或缺少运行时依赖。检查cplex\matlab下是否存在对应Matlab大版本如R2023a的文件。我在某次用新版本MatlabR2023b时CPLEX老版本里根本没有对应mexw64只能换一个支持的Matlab版本。玄学的是有时多装一个Visual C Redistributable就能好这也是我建议“先补运行时再考虑重装”的原因。4.2 约束矩阵行列错误导致模型与预期不符现象求解不报错但结果里负荷平衡误差大、机组出力和预期完全不同、备用约束看起来失效。原因几乎都出在sparse矩阵组装时行索引或列索引写错最常见的两种行号row没有递增导致约束互相覆盖ui和pi索引计算时忘了加nU偏移。解决在组装循环里加一句断言或者直接小规模打印检查% 只打印第1个时段的行检查非零元素位置是否符合预期 t 1; row 之后直接用full(A(row, :))查看 fprintf(负荷行非零列索引: %s\n, mat2str(find(A(1, :))));正常逻辑下负荷平衡行在第1行非零列应该是nU1机组1的p1,t和nUT1机组2的p1,t。如果这里出现的是1或T开头的索引说明构造时偏移量丢了。这类问题用LP文件更能一眼看出来最后一章专门讲。另一个实用习惯是构造完A后做一次nnz(A)检查规模固定时非零元素数应该完全可预期。4.3 二次目标Q矩阵系数偏差成本对不上现象求解器返回的最优成本比我手算的期望值低一半左右出力分配也偏向大出力。原因十有八九是CPLEX的目标函数约定。CPLEX文档里二次目标写的是minimize 0.5 * x * Q * x f * x你对目标里是a * p^2所以Q对角元素必须填2 * a。填成a时求解器实际优化的目标是0.5 * a * p^2相当于把二次成本砍半低成本的大出力机组当然会被过度使用。解决Q(pi, pi) 2 * a(i); % 手写注释CPLEX目标是0.5*xQx所以补系数2这个坑的隐蔽性在于模型不报错、解“看起来合理”但成本数字和数据对不上。我的排查习惯是先固定所有机组为开机状态把模型降成纯QP用已知解析解或手算的小算例比如单机组单时段对比目标值能快速确认Q系数。4.4 模型报Infeasible从等式约束到数值尺度逐项查现象cplex.Solution.status返回infeasible。这是入门者最崩溃的一个状态其实排查顺序是固定的。第一步检查负荷平衡等式是否前后一致比如所有时段D(t) R(t)大于最大开机容量组合或者某时段D(t)小于所有机组的Pmin之和这在数学上直接无解。第二步逐个放开约束定位最小冲突集常见做法是暂时删除备用约束跑一次如果变可行问题出在备用约束再删除出力上下限约束跑一次如果变可行问题在上下限联动。第三步检查数值尺度。成本系数从0.015到700幅度接近5个数量级CPLEX的默认数值容差在某些条件下会遇到困难。我一般先把所有成本统一再归一化或者调大cplex.Param.mip.tolerances.integrality.Cur从1e-6改到1e-5。调容差是最后的兜底手段不是第一选择任何工业级应用都必须以“模型本身正确”为前提。4.5 求解时间爆炸用mipgap和timelimit控场现象2台机组24时段跑得飞快但把T改成168一周后求解时间从秒级跳到分钟级甚至不回。这让很多初学的人误以为是“代码效率低”其实是MIP天然的组合复杂性每台机组每个时段都有开/关两个状态24时段48个二进制变量已经有2^48量级的组合空间。工业模型靠分支定界剪枝和启发式不靠“优化代码”。控制求解时间的方法是设置两个参数cplex.Param.mip.tolerances.mipgap.Cur 0.01; % 允许1%最优性间隙 cplex.Param.timelimit.Cur 60; % 最多60秒 cplex.Param.threads.Cur 8; % 用8线程并行参数说明mipgap表示当找到的解与下界差距小于1%时就停。对调度场景1%的成本偏差完全可接受但求解时间可能减少80%以上。timelimit是硬性上限任何情况下不超时。我的习惯是调试阶段用4个时段验证算法正式跑用完整时段配0.5%~1%的gap。记住无脑等待“完美最优解”不是工程思维调度问题里“足够好的解给定时间内拿出来”才是核心诉求。5. 从2变量扩展到实际调度N台机组、热启动与一个诊断技巧5.1 把2台机组改成N台通用数据结构写法2台机组所有索引都是手写的扩展到N台最忌讳的是继续手写三层循环。第一步是建立统一的机组参数表结构常见做法是用结构体数组或表格table替代标量向量。把所有和机组相关的参数放进一个struct索引规律用变量名固化下来uc struct(); uc.Pmin [100; 150; 80]; % 新增机组3 uc.Pmax [300; 400; 180]; uc.a [0.020; 0.015; 0.025]; uc.b [20; 17.5; 24]; uc.c [120; 100; 60]; uc.SU [500; 700; 300];然后所有组装代码里的Pmin(i)都改成uc.Pmin(i)循环上限从nUnit2改成nUnitlength(uc.Pmin)。这样增加机组时只改数据结构不动组装逻辑。更工程化的做法是把变量索引和机组索引用函数封装起来比如ui (i,t) (i-1)*Tt避免在四五处重复写(i-1)*Tt导致改一处漏三处的惨剧。5.2 用MIP Start热启动加速重复求解实际调度经常是一个滚动优化的过程上一时段算完下一时段在新负荷数据下重新求解。每次从零开始分支定界很浪费时间。CPLEX的Matlab接口支持MIP start把上一次的整数解作为本次搜索的初始解try % x_last是上一轮求解的最优解维度与本轮完全一致 cplex.addMipStart(x_last, rolling_start); catch % 老版本不支持时跳过不影响正确性 endMIP Start的作用机理是分支定界算法一开始就有上界剪枝效率大幅提升尤其当日负荷曲线相似时上一轮解距离最优解很近搜索空间可以大幅缩小。注意前提是两次建模的变量顺序、数量和机组参数完全一致否则初始解会被判定非法。如果你用的是YALMIP做建模层可以用assign给变量赋初值然后指定cplex求解效果类似但要求YALMIP版本较新。滚动场景下我能拿到的实际收益是24时段模型从平均40秒降到12秒左右代价只是多传一个向量这买卖划算。5.3 诊断技巧任何可疑模型先写成LP文件再说话最后说一个我每次排查都会做的习惯无论模型报错还是结果可疑先把模型写入LP文件用文本编辑器打开逐行读约束。CPLEX的Matlab接口里一行命令的事cplex.writeModel(uc_diagnose.lp);LP文件是人类可读的纯文本它会把你组装的稀疏矩阵还原成带变量名比如x1、x2的约束表达式。排查诀窍是找你看得懂的锚点负荷平衡行应写成x49 x73 450这样的形式x1和x2是启停变量时它们的Bin声明必须出现在文件头顶Q矩阵系数会显示在“Quadratic Objective”段落。任何索引错位、系数乘2问题、边界写反在LP文件里都会以最原始的代数形式暴露出来。从这个意义上说writemodel比任何调试器都诚实。我自己的流程是数据构造完先写出LP文件跑一遍检查约束形状然后求解结果不对再写一遍LP文件对比差异。这个习惯陪我排掉了至少三分之二的“玄学”问题。如果你想把2变量骨架做成一个能长期复用的调度模型把“跑通”到“可信”的这一步补上固定一个已知负荷场景手算或按比例缩放一个预期成本区间模型输出落在区间内才叫验证通过。这个验证步骤做一次后面换数据、换机组参数都会更有底气希望帮到你。本文还有配套的精品资源点击获取