企业官网建设流程全解析

搞电力系统优化的，十有八九会碰到这么个场面：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和Bbus：

Cg = 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 求解器接入失败的排查清单

环境搭建阶段的问题，我整理成了一张速查表，基本上覆盖了绝大多数情况。

接入失败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”的警告。这时候先把所有数据转成标幺值，用baseMVA=100统一缩放，再看结果是否恢复稳定。

如果模型不可行，问题更容易定位。YALMIP里求解后可以直接用check(cons)查看约束残差，数值为正表示满足约束，负数绝对值越大说明违反越严重。把每条约束残差打印出来，通常马上就能找到瓶颈：是发电机下限总和大于负荷，还是某条线路容量设成了0。还有一点，DC-OPF里如果松弛母线的相角没有被固定，平衡约束会有无穷多解，求解器虽然可能给出解，但结果不稳定，所以theta(refIdx) == 0这条一定不能少。

5.3 想上 AC-OPF 和混合整数模型怎么办

CPLEX和Gurobi强大归强大，但AC-OPF这种非凸问题它们并不擅长。交流潮流方程里有电压幅值与相角的乘积项，不是二次规划，更不是线性规划，直接丢给CPLEX/Gurobi会报“不支持的非线性约束”。想继续用这套工具链，有两个方向：一是把AC潮流线性化，用LPAC或者DistFlow近似；二是换一个能处理非线性的求解器，比如IPOPT，YALMIP同样可以调用。不过从工程实践看，绝大多数效率优化场景用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的作用是：机组停机时u=0，强迫出力为0；机组开机时u=1，出力恢复上下限约束。目标函数里的c .* u是空载成本，整个问题变成混合整数二次规划，CPLEX和Gurobi处理这种模型正是拿手好戏。从DC-OPF到UC，核心代码变化很小，这也是我推荐用这套工具链的原因。

5.4 个人实操心得

样板戏演完，最后说几句实在话。我自己的习惯是：凡是涉及IEEE39的优化，第一步永远是跑一个不带网络约束的经济调度，看总成本和机组出力是否合理；第二步再加载DC-OPF，用check(cons)检查约束；第三步才放整数变量或转AC。这样一层层往上加，出了问题永远能定位到是数据、约束还是求解器。

还有一个更偏个人偏爱的技巧：把sdpsettings('solver')做成一个循环，分别跑CPLEX和Gurobi，把两个求解器的目标值、求解时间、迭代次数存下来。导师或合作者问起为什么选某个求解器时，直接甩一张对比表，比任何口舌都有说服力。IEEE39算例规模小，跑一遍往往不到一秒，多跑几次完全没成本。希望这篇内容能让你少走几趟弯路。