搞电力系统优化的十有八九会碰到这么个场面Matlab里装着Matpower手头是IEEE39节点系统想在论文里加一个经济调度或者最优潮流算例结果被各种求解器折腾到深夜。我自己的固定搭配是Matlab YALMIP CPLEX/Gurobi用IEEE39做试验田。这篇文章不绕弯子直接讲清楚怎么把这套链路搭起来以及建模求解过程中那些文档里不会写的坑。适合刚入门优化调度的研究生也适合想快速在标准算例上验证算法的工程师。1. 从需求到方案为什么是 YALMIP CPLEX/Gurobi 的组合1.1 三层分工建模层、求解层和算例层很多朋友第一次接触这套工具链时会觉得混乱YALMIP、CPLEX、Gurobi、Matlab、IEEE39五个词堆在一起到底谁是主角先把这个关系理顺。Matlab是宿主环境负责存数据、写脚本、画图YALMIP是优化建模工具箱提供变量、约束、目标函数的高级描述CPLEX和Gurobi是真正干活的商业求解器负责在约束围成的可行域里找出最优解IEEE39则是用来检验整个流程的标准测试系统。打个比方YALMIP是翻译官把“发电机出力必须在上下限之间”这种人类语言翻译成求解器能读的标准数学形式CPLEX/Gurobi是发动机装上了才能跑出结果。三者加在一起才能把“算一个IEEE39算例”变成一件顺手的事。有人会问直接用Matlab的linprog、quadprog不行吗不是不行而是一旦模型变复杂手工构造系数矩阵的成本会失控。拿DC-OPF来说你要自己拼节点电纳矩阵、发电机关联矩阵、线路潮流约束矩阵每次改补偿母线或者加一个整数变量矩阵维度就全变了。YALMIP的价值就在于把“矩阵维度”这件事隐含掉你只需要声明变量是sdpvar还是binvar再写约束关系剩下的事它来处理。对于IEEE39这种几十个变量的算例看起来没什么但在实际项目中模型复杂度远不止如此早一天切换到YALMIP就早一天摆脱手撸矩阵的噩梦。1.2 为什么拿 IEEE39 当第一个算例IEEE39节点系统也叫新英格兰10机39节点系统是一套被广泛使用的标准测试算例。它包含39条母线、10台发电机和46条交流线路负荷总量在6000MW量级。相比IEEE14、IEEE30它的网络结构更接近真实电网节点和支路数又不会大到让调试变成灾难相比IEEE118它又足够简单适合把模型细节一条条抠清楚。所以我在跑新算法时通常先拿IEEE39验证可行性再往更大的系统迁移。这个算例在Matpower里就叫case39可以用一句话加载进来这也是我推荐新手从这里入手的原因。真正做优化时IEEE39还有一个好处它的机组类型和成本参数是分开的发电机有功上下限、母线负荷、线路容量都齐全既能支撑只关心机组出力的经济调度也能支撑考虑网络约束的最优潮流。你可以在同一套数据上反复改模型不用到处找数据。对于CPLEX和Gurobi来说IEEE39的维度不大LP或QP通常几秒内就能解完所以即使模型写错了也能很快从解的结果中发现问题而不是干等求解器跑十分钟。1.3 你可能在做的三类问题经济调度、最优潮流和机组组合围绕IEEE39可以做的优化问题主要有三类难度依次递增。第一类是经济调度ED只考虑各台发电机的成本函数和出力上下限目标是把总发电成本压到最低网络结构完全忽略适合练手和理解成本函数。第二类是直流最优潮流DC-OPF在ED基础上加入直流潮流方程和线路容量约束让节点注入功率必须满足“流入等于流出”的网络规律这是电力系统优化最常用的线性模型。第三类是机组组合UC再叠加整数变量用0/1表示机组开停变成混合整数二次规划CPLEX和Gurobi的整数求解能力在这里就能真正体现出来。要注意的是CPLEX和Gurobi擅长的是线性规划、二次规划、混合整数线性/二次规划这类凸问题。AC-OPF里包含电压幅值与相角相乘的非线性项属于非凸问题这两个求解器并不能直接高效处理这也是很多新手踩坑的地方。所以我倾向于用DC-OPF作为主线讲解把工具链和建模方法讲透之后要扩展成SCUC安全约束机组组合也只需要修改约束和整数变量不会动摇整个流程。2. 环境搭建让 Matlab 真正“认得出”两个求解器2.1 YALMIP 安装与路径配置YALMIP本身是一堆m文件不需要编译安装流程非常简单。从GitHub或官网下载最新的yalmip.zip解压到某个固定目录比如D:\toolbox\yalmip然后在Matlab命令行里执行addpath(genpath(D:\toolbox\yalmip)); savepath;genpath会递归把这个目录下所有子目录加入搜索路径savepath会把路径保存到startup文件这样下次启动Matlab不用重复设置。这里有一点要提醒不同版本的Matlab对YALMIP的兼容性略有差别如果你用的是很老的YALMIP版本在较新的Matlab上可能报一些莫名其妙的数组兼容错误所以尽量保持YALMIP为最近一年内的版本。不要嫌更新麻烦求解器接不上很多时候不是求解器问题而是YALMIP太老。2.2 CPLEX 接入从安装目录到 MATLAB 路径CPLEX现在归属于IBM ILOG CPLEX Optimization Studio。商业用户用正式许可学生或教学场景可以下载社区版社区版免费但有变量和约束数量限制大概是1000个对于IEEE39算例完全够用。安装完成后关键是找到Matlab接口目录。不同版本的目录名规律是C:\Program Files\IBM\ILOG\CPLEX_Studio2210\cplex\matlab\x64_win64其中2210对应2022.1版本后面的x64_win64对应64位Windows。找到目录后在Matlab中执行addpath(C:\Program Files\IBM\ILOG\CPLEX_Studio2210\cplex\matlab\x64_win64); savepath;如果使用网络版许可证还需要设置环境变量ILOG_LICENSE_FILE指向许可文件位置否则求解器虽然能识别但一调用就会报许可证找不到。这里最常见的问题是有人只装了CPLEX的Python接口或通用可执行文件没在Matlab路径中加入cplex\matlab目录导致YALMIP始终提示找不到CPLEX。检查路径是否加对其实是低垂的果实却经常被忽略。2.3 Gurobi 接入gurobi_setup 和许可证Gurobi的安装流程比CPLEX更顺手。安装包会默认放到如C:\gurobi1100\win64这样的目录其中1100代表版本号11.0.0。以我常用的Gurobi为例进入win64\matlab目录里面有一个现成的gurobi_setup.m脚本直接在Matlab里运行cd(C:\gurobi1100\win64\matlab); gurobi_setup这个脚本会把Matlab接口路径自动加入搜索路径。接下来需要处理许可证。学术用户可以注册Gurobi账号申请免费学术许可然后按官网提示用grbgetkey把许可写到本机。如果是单机浮动许可最好手动配置环境变量GRB_LICENSE_FILE指向gurobi.lic文件所在目录。好多人在这一步卡住原因是许可文件在服务器上而本地没有告诉Gurobi去哪找。记住Gurobi安装成功不等于许可可用gurobi_setup只负责路径不负责许可证。2.4 用 yalmiptest 做一次体检环境是否配置完成不需要自己写复杂的测试模型YALMIP自带一个体检命令yalmiptest。在Matlab命令窗输入它YALMIP会依次测试支持的求解器并在结果列表里显示CPLEX、Gurobi等是否可用。我一般只看对应的行是不是“Success”如果是说明求解器可以被正常调用如果是“Failure”或“Solver not found”说明路径或许可出了问题。这种体检比较耗时间但首次配置时非常值。除了yalmiptest也可以写一个三行的快速测试x sdpvar(1,1); result optimize(x 1, x^2, sdpsettings(solver,gurobi)); if result.problem 0, disp(Gurobi OK); end把solver换成cplex就是CPLEX的测试。注意sdpsettings是YALMIP里设置求解器选项的入口后面所有调用都会用到。这里有个经验如果两个求解器都接了在真正求解IEEE39前先跑一次这个最小测试能帮你确认不是环境的问题再开始怀疑模型。YALMIP的求解器调度逻辑是如果你不指定solver它会自动选择一个已注册的可用求解器但为了可控我建议每次都显式指定。3. IEEE39 算例建模把电路图纸翻译成优化约束3.1 case39 数据结构速览Matpower里加载IEEE39最常用的命令是loadcase(case39)。返回的mpc是一个结构体里面最关键的是四个矩阵bus存放母线数据每一行是一根母线包括母线编号、有功负荷、无功负荷、电压上下限等gen存放发电机数据每一行是一台发电机包括挂接母线编号、当前出力、有功上下限等branch存放线路和变压器支路数据包括首末端母线、电阻、电抗、长期容量等gencost存放成本函数系数。以gen矩阵为例gen(:,10)是有功出力下限gen(:,9)是有功出力上限这两列是DC-OPF必然要用的。bus矩阵里的bus(:,3)是母线有功负荷bus(:,2)是母线类型其中数值3通常代表松弛母线。我看很多初学者一拿到case39就着急写模型结果对数据列的含义全靠猜。我的建议是先把这些矩阵打出来看一遍用disp(mpc.bus(1:5,:))肉眼确认格式再用代码把发电机所在的母线号提取出来genBus mpc.gen(:,1); isSlack mpc.bus(:,2) 3; refIdx find(isSlack);这么做能避免后面建模时出现母线编号和矩阵行号对不上的低级错误。DC-OPF里松弛母线必须存在一个参考角度通常取case39的31号母线作为平衡节点如果你的case39数据版本不同可能会稍有差异所以用bus(:,2)3动态找出来要比写死31更稳妥。3.2 直流潮流DC-OPF的建模推导DC-OPF要做的事情一句话概括在所有发电机出力和线路传输容量都满足限制的前提下找到让总发电成本最小的潮流分布。这里的“潮流”用一个线性方程近似描述即忽略无功、电阻和电压变化认为线路有功功率只取决于两端相角差和线路电抗P_ij (theta_i - theta_j) / x_ij那么每个节点上注入的有功功率发电机出力减去负荷必须等于所有从该节点流出的线路功率之和写成矩阵形式Cg * Pg - Pd B * theta其中Cg是发电机关联矩阵把每台发电机的出力放到它所在的母线位置上B是节点电纳矩阵由所有线路电抗取倒数后拼装而成。因为B是一个奇异矩阵如果没有参考角度方程会有无穷多解所以必须额外加一个约束theta(refIdx) 0固定松弛母线的相角。线路容量约束也很直观任意一条线路上的有功潮流必须在额定容量范围内-Fmax_ij (theta_i - theta_j) / x_ij Fmax_ij目标函数是最小化所有发电机的总成本。如果每台机组成本是二次函数a_i*Pg_i^2 b_i*Pg_i c_i那么总成本是凸二次函数CPLEX和Gurobi可以精确求解如果二次项系数为0问题退化为线性规划。这一点对求解器非常友好也是为什么DC-OPF被广泛应用的原因。3.3 YALMIP 代码实现变量、约束、目标一次写清在YALMIP里建模最大的感受是“想做收敛代码时先把几个核心对象声明明白。第一步定义决策变量Pg sdpvar(ng, 1); theta sdpvar(n, 1);Pg是一个ng维向量表示每台发电机的有功出力theta是一个n维向量表示每根母线的相角。sdpvar是YALMIP声明连续变量的函数如果遇到机组组合问题把某几个变量声明为binvar即可变成整数变量。第二步构建节点电纳矩阵和发电机关联矩阵这一步是纯数据预处理建议写成独立函数便于复用。第三步写约束直接使用[]连接cons [Cg * Pg - Pd Bbus * theta, Pmin Pg Pmax, theta(refIdx) 0];YALMIP会把这种向量化约束拆成逐条标量约束比你手写for循环拼矩阵省太多。最后设置目标函数cost sum(a .* Pg.^2 b .* Pg); optimize(cons, cost, sdpsettings(solver,cplex));注意YALMIP的optimize输入顺序是“约束、目标、选项”很多新人在初学时容易把顺序写成“目标、约束、选项”这一行放错就会报错。写完后用value(Pg)取回数值解整个DC-OPF的核心工作就完成了。4. 完整实操从读取 matpower 数据到打印调度结果4.1 数据读取与矩阵预处理现在把前面讲的串成一个能直接跑的脚本。第一步读取case39数据并提取维度信息mpc loadcase(case39); n size(mpc.bus, 1); ng size(mpc.gen, 1); nl size(mpc.branch, 1); Pd mpc.bus(:, 3); Pmin mpc.gen(:, 10); Pmax mpc.gen(:, 9); branch mpc.branch; active branch(:, 11) 0; branch branch(active, :); nl size(branch, 1); Fmax branch(:, 6); Fmax(Fmax 0) 1e6; % 母线编号到矩阵行号的映射 [~, fb] ismember(branch(:,1), mpc.bus(:,1)); [~, tb] ismember(branch(:,2), mpc.bus(:,1)); [~, genBusIdx] ismember(mpc.gen(:,1), mpc.bus(:,1));这里我把branch里非运行支路过滤掉branch(:,11)是状态位置1才代表投入。Fmax取rateA有些数据里rateA为0表示不限容量所以我直接把它换成一个足够大的数。ismember返回的是母线编号在bus矩阵中的行号后面用它索引theta就不会出错。接下来构建Cg和BbusCg zeros(n, ng); for k 1:ng Cg(genBusIdx(k), k) 1; end Bbus zeros(n, n); for k 1:nl i fb(k); j tb(k); x branch(k,4); if x 0 Bbus(i,i) Bbus(i,i) 1/x; Bbus(j,j) Bbus(j,j) 1/x; Bbus(i,j) Bbus(i,j) - 1/x; Bbus(j,i) Bbus(j,i) - 1/x; end end这一段看着繁琐但逻辑很简单电抗取倒数就是线路电纳组装进节点电纳矩阵。为什么要这么写因为YALMIP只负责优化部分矩阵预处理还得自己来。如果把四舍五入错误、母线编号错位这些坑留到优化阶段排查起来会非常痛苦。我习惯把这两段预处理单独放进build_bdc.m函数等模型越来越大时能省下不少时间。4.2 完整 DC-OPF 求解脚本的主体数据准备好后正式的YALMIP建模就简洁多了refIdx find(mpc.bus(:,2) 3); refIdx refIdx(1); theta sdpvar(n, 1); Pg sdpvar(ng, 1); a 0.001 * ones(ng, 1); b 20 * ones(ng, 1); cost sum(a .* Pg.^2 b .* Pg); cons [Cg * Pg - Pd Bbus * theta, Pmin Pg Pmax, theta(refIdx) 0]; for k 1:nl f (theta(fb(k)) - theta(tb(k))) / branch(k,4); cons [cons, -Fmax(k) f Fmax(k)]; end ops sdpsettings(solver, gurobi, verbose, 2); sol optimize(cons, cost, ops);这里refIdx是通过母线类型找到的在case39中通常就是31号母线。目标成本系数a、b是简化的示例值实际使用时应从gencost或你自己的数据中读取否则结果只具有教学意义。如果你愿意也可以改为solver, cplex只要路径和许可证没问题两个求解器都能处理这个线性二次规划。注意solver参数要在sdpsettings里指定YALMIP对大小写不敏感但为了保险建议统一小写。verbose是输出详细程度设为2可以看到求解日志出问题时比设0更容易定位。4.3 结果解析与可视化求解完成后第一件事是检查sol.problem它返回0表示求解成功。然后取出结果if sol.problem 0 Pg_opt value(Pg); theta_opt value(theta); total_cost value(cost); line_flow (theta_opt(fb) - theta_opt(tb)) ./ branch(:,4); loading abs(line_flow) ./ Fmax; disp([总出力(MW): , num2str(sum(Pg_opt))]); disp([总负荷(MW): , num2str(sum(Pd))]); disp([总成本: , num2str(total_cost)]); end在DC-OPF模型里因为没有线路损耗总出力会等于总负荷。如果这个等式不成立说明平衡约束写错或某些支路被过滤错了。我通常还会画两张图一张是发电机出力柱状图用于和机组上下限对比另一张是线路负载率图用于检查有没有线路越限。bar(Pg_opt); ylabel(有功出力/MW); figure; plot(loading, o); hold on; yline(1, r--); ylabel(线路负载率);画出来之后一眼就能看出哪些发电机接近上限、哪些线路重载这对验证模型可靠性很有帮助。千万不要只盯着总成本一个数看那是最容易掩盖错误的方式。我在实际算例中踩过很多次总成本算出来低得离谱一查才发现某台发电机的下限被写成负数或者线路容量全被设成了Inf。5. 常见问题速查表 个人避坑心得5.1 求解器接入失败的排查清单环境搭建阶段的问题我整理成了一张速查表基本上覆盖了绝大多数情况。现象可能原因解决办法yalmiptest显示CPLEX或Gurobi Not foundMatlab路径里没有加入接口目录重新执行addpath确认目录存在optimize报No suitable solver没有在sdpsettings指定solver或指定名称写错显式写solver, gurobi检查拼写调用求解器时报许可证错误许可证没有安装或环境变量没设置检查GRB_LICENSE_FILE或ILOG_LICENSE_FILE解出来的值全是NaN模型不可行或数值严重病态检查约束是否矛盾缩放单位查看sol.info求解器非常慢整数变量过多或模型非凸先求解不含整数的松弛版本看看是否秒出接入失败90%是路径问题10%是许可问题。先跑yalmiptest再跑最小测试。不要随意把求解器升级到最新版YALMIP可能跟不上最新接口我遇到过Gurobi新版本发布时旧版YALMIP调用报内部错误升级YALMIP后解决。这里有个小技巧在optimize之前用ops sdpsettings(solver,gurobi,gurobi.NumericFocus,3)这个参数会告诉Gurobi提高数值精度代价是增加一点求解时间当模型出现无意义的微小扰动时很管用。5.2 数值病态与求解卡死怎么办DC-OPF本身是一个线性问题但数据预处理里有一个很隐蔽的坑线路电抗太小会让节点电纳矩阵里出现数值非常大的元素矩阵条件数变差求解器在迭代时容易出现数值错误。Matpower的case39数据本身是标幺值一般问题不大但如果你自己改成有名值或者把某条线路的参数填错解出来的相角可能离谱、sol.problem报一些“numerical issues”的警告。这时候先把所有数据转成标幺值用baseMVA100统一缩放再看结果是否恢复稳定。如果模型不可行问题更容易定位。YALMIP里求解后可以直接用check(cons)查看约束残差数值为正表示满足约束负数绝对值越大说明违反越严重。把每条约束残差打印出来通常马上就能找到瓶颈是发电机下限总和大于负荷还是某条线路容量设成了0。还有一点DC-OPF里如果松弛母线的相角没有被固定平衡约束会有无穷多解求解器虽然可能给出解但结果不稳定所以theta(refIdx) 0这条一定不能少。5.3 想上 AC-OPF 和混合整数模型怎么办CPLEX和Gurobi强大归强大但AC-OPF这种非凸问题它们并不擅长。交流潮流方程里有电压幅值与相角的乘积项不是二次规划更不是线性规划直接丢给CPLEX/Gurobi会报“不支持的非线性约束”。想继续用这套工具链有两个方向一是把AC潮流线性化用LPAC或者DistFlow近似二是换一个能处理非线性的求解器比如IPOPTYALMIP同样可以调用。不过从工程实践看绝大多数效率优化场景用DC-OPF就够了AC-OPF更多用于校核和详细潮流分析。如果要做机组组合反而是在DC-OPF上很自然的扩展。把Pg对应的开关量用binvar声明再增加开机状态约束u binvar(ng, 1); cons [Pmin .* u Pg Pmax .* u, Cg * Pg - Pd Bbus * theta, theta(refIdx) 0]; cost sum(a .* Pg.^2 b .* Pg c .* u);这里Pmin .* u Pg Pmax .* u的作用是机组停机时u0强迫出力为0机组开机时u1出力恢复上下限约束。目标函数里的c .* u是空载成本整个问题变成混合整数二次规划CPLEX和Gurobi处理这种模型正是拿手好戏。从DC-OPF到UC核心代码变化很小这也是我推荐用这套工具链的原因。5.4 个人实操心得样板戏演完最后说几句实在话。我自己的习惯是凡是涉及IEEE39的优化第一步永远是跑一个不带网络约束的经济调度看总成本和机组出力是否合理第二步再加载DC-OPF用check(cons)检查约束第三步才放整数变量或转AC。这样一层层往上加出了问题永远能定位到是数据、约束还是求解器。还有一个更偏个人偏爱的技巧把sdpsettings(solver)做成一个循环分别跑CPLEX和Gurobi把两个求解器的目标值、求解时间、迭代次数存下来。导师或合作者问起为什么选某个求解器时直接甩一张对比表比任何口舌都有说服力。IEEE39算例规模小跑一遍往往不到一秒多跑几次完全没成本。希望这篇内容能让你少走几趟弯路。