简介：基于二阶锥规划的主动配电网最优潮流求解MATLAB代码，面向电气工程、计算机及数学等专业学生，用于课程设计、期末大作业或毕业设计。代码采用参数化编程，注释详尽，便于调整参数与二次开发；压缩包内含20个文件，以两个MATLAB脚本为核心，另附PDF与PPT讲解文档、IEEE33节点拓扑图、参考文献及辅助文本，整体约10.95MB。内容覆盖风电、无功补偿、储能、有载调压变压器等多时段24小时协调优化场景，附带可直接运行的案例数据，并提供数据导入、模型构建、求解器配置与结果分析等完整流程，能够帮助读者理解二阶锥规划如何处理潮流约束，掌握主动配电网优化调度的建模与求解方法。目前已有130人学习下载，适合需要从理论到实践快速入门的进阶学习者，也可为相关课题研究提供可复用的算法框架，具备较强的工程参考价值。

1. 二阶锥规划做主动配电网最优潮流：这套 MATLAB 代码解决的不只是「能算」

二阶锥规划（SOCP）这几年在主动配电网最优潮流里几乎成了默认选项。原因很直接：配电网最优潮流本身是非凸问题，启发式算法只能给你一个无法证明优劣的可行解，而 SOCP 凸松弛能从数学上保证全局最优。这份 MATLAB 代码就是把整套方法落到 IEEE 33 节点算例上的完整实现——从 DistFlow 建模、二阶锥松弛，到求解器调用和结果输出，每一步都有对应可执行代码。适合两类人：做配电网方向毕设或论文的研究生，以及想快速搭出主动配电网优化调度原型的一线工程师。

2. 模型与松弛原理：DistFlow 方程、变量代换与二阶锥松弛到底松在哪

2.1 从交流潮流到 DistFlow：径向配电网怎么描述

配电网和输电网最大区别在于拓扑：输电网是环网，配电网正常运行时是辐射状的。这意味着牛拉法那种对所有节点做复数电压迭代的通用解法虽然能用，但对于优化场景来说太重了——你要优化的不是一组工作点，而是在大量约束下搜索最优工作点，模型里多一个非凸变量，求解难度就上一个台阶。

常见做法是用 DistFlow 支路潮流模型。它其实是交流潮流在辐射状网络下的等价改写：对每条支路(i, j)，定义从 i 流向 j 的有功P_ij和无功Q_ij，同时用电压幅值平方U_i = V_i^2、电流幅值平方L_ij = I_ij^2代替原来的复数电压和相角。节点功率平衡写成都一样，但支路电压降方程变成了这种形式：

U_j = U_i - 2(R_ij * P_ij + X_ij * Q_ij) + (R_ij^2 + X_ij^2) * L_ij

这个方程的好处是它直接描述有功、无功、电压幅值之间的关系，不引入相角。配电网里相角差本来就不大，把这个自由度去掉不会损失工程精度，却把问题的非凸性一下子压缩到一个二次耦合式上，也就是这条支路自身的功率与电压电流关系。

P_ij^2 + Q_ij^2 = U_i * L_ij

注意，这是一个精确关系，不是近似。DistFlow 模型在辐射状网络上和完整交流潮流是等价的，区别只是换了一组变量。这份资源里的代码采用的就是这个模型，你在主程序里会看到U、L、P、Q四个核心变量，后面所有约束都围绕它们展开。

2.2 二次耦合项怎么「松」成锥：变量替换和松弛的几何意义

上面那个等式P^2 + Q^2 = U * L是整个模型里最后一块非凸拼图。等式的可行域是一条双曲曲面，把它直接丢给求解器，问题仍然是 NP-hard 级别的。

二阶锥规划的核心操作，就是把等式「松」成不等式：

P_ij^2 + Q_ij^2 <= U_i * L_ij

这个不等式在数学上是一个旋转二阶锥（rotated second-order cone），是凸集。为了说明这个松弛是有意义的，可以把旋转锥等价地写成标准二阶锥形式：

norm([2*P_ij; 2*Q_ij; U_i - L_ij]) <= U_i + L_ij

也就是说，原来要求功率点精确落在曲面上，现在允许它落在曲面「上方」的锥体内。问题变成凸优化了，求解器的速度和稳定性都大幅提升。

但这里有个关键的工程问题：松弛后的最优解不一定在原曲面上，可能落在锥体内部，此时解虽然可行，却不满足真实潮流方程。业界把这种现象叫做「松弛不精确」或「锥松弛间隙不为零」。好在绝大多数配电网最优潮流场景下，只要目标函数里包含网损项或电压惩罚项，最优解会被「拉」到锥边界上，松弛是精确的。这也是为什么代码里目标函数通常会显式写出网损最小——不只是为了节能，也是为了保障松弛有效性。

提示：如果目标函数完全不含网损项（比如纯最小化 DG 购电成本），锥约束可能松掉。第 4 章会专门讲怎么检查和补救。

2.3 主动配电网引入 DG 后，为什么 SOCP 比启发式算法更合适

传统配电网是被动网络，优化变量就那么几个：变压器抽头、无功补偿容量。到了主动配电网（ADN）场景，分布式电源出力、储能充放电功率、可调负载、联络开关状态全都成了变量，决策空间维度一下子从十几个涨到上百个。

我见过不少人在这个场景下用遗传算法或粒子群，结果往往很尴尬：同一组参数跑三次出来三个不同的网损值，而且谁都说不清是不是全局最优。启发式算法对初值敏感、对罚函数系数敏感，调试成本全堆在「玄学调参」上。

SOCP 在这个场景下的优势有三条。第一，凸问题的局部最优就是全局最优，结果可复现，写论文、写报告都站得住脚。第二，YALMIP、Mosek、Sedumi 这套工具链非常成熟，配电网规模（几十到几百个节点）对它们来说是小问题，求解时间一般是秒级到分钟级。第三，目标函数无论是网损、购电成本、弃风弃光惩罚，还是这些项的加权组合，都能保持凸性，后续往多时段、网络重构方向扩展也不会破坏模型结构。

有人会问：为什么不直接用非线性求解器（fmincon、ipopt）解原问题？能解，但要面对初值、局部极小值、约束雅可比奇异这些麻烦；而且非线性内点法每次求解结果依赖初值，本质上没有解决「可复现性」问题。SOCP 相当于用「贴上精确性标签的松弛」换来了「可证明的全局最优」，这笔交易在工程上是划算的。

3. 把 MATLAB 代码跑通：YALMIP 建模、求解器调用与结果输出

3.1 代码结构速览：算例数据、主程序、后处理怎么分工

拿到代码后我习惯先不看主程序，先把文件按「数据 / 建模 / 后处理」三组分开。这份代码的典型组织方式是：一个数据文件保存 IEEE 33 节点算例（节点编号、支路阻抗、负荷有功无功、DG 接入位置），一个主程序完成 YALMIP 建模与求解，再加上一两个画图或验证的子函数。

IEEE 33 节点算例的公共参数我列一下，方便你对照数据文件核对基准值：

这套参数在网络公开文献里很常见，属于标准算例。如果数据文件里出现了不认识的列，先对照这张表确认有没有把 R、X、P、Q 的列读对，我见过有人把负荷有功列读成无功列，结果求解器报 infeasible 还不知道哪错了。

3.2 变量定义与目标函数：网损最小的两种写法

建模第一步是把决策变量用 YALMIP 的sdpvar声明出来。这份代码里核心变量是电压平方、电流平方、支路功率和节点注入功率：

n = 33; % 节点数 nBranch = 32; % 支路数，不含联络开关 U = sdpvar(n, 1); % 节点电压幅值平方 U_i L = sdpvar(nBranch, 1); % 支路电流幅值平方 L_ij P = sdpvar(nBranch, 1); % 支路有功，方向由 from 指向 to Q = sdpvar(nBranch, 1); % 支路无功 Pg = sdpvar(n, 1); % 节点注入有功（DG + 变电站） Qg = sdpvar(n, 1); % 节点注入无功

U、L是平方量纲，整个模型都以它们为变量，不需要再写V^2这类非线性表达式。P、Q的维数是支路数，不是节点数，写约束时注意下标对齐。Pg、Qg里根节点对应的是变电站注入，DG 接入节点对应的是分布式电源出力，没有 DG 的节点该项会被限成 0。

目标函数最直接的写法是网损最小。每条支路损耗是R_ij * L_ij，所以：

% 支路电阻向量，单位 p.u. R = branch(:, 3) / Zbase; Objective = R' * L; % 总网损 = sum(R_ij * L_ij)

另一种写法是把网损展开成「流入功率减去流出功率」的形式，即sum(P(in)) - sum(sum(P(out)))，但数值上两种写法等价，直接R' * L更简洁，而且显式把网损写进目标能压制锥松弛间隙，后续验证环节省不少麻烦。如果要做多目标（网损 + DG 成本），就写成Objective = R'*L + c' * Pg(dgIdx)，注意把 DG 成本系数c的量纲折到和网损同一基准。

3.3 约束逐条落进 YALMIP：潮流、电压、DG 出力与锥约束

约束是这份代码的骨架，我按四类展开。先预计算每条支路的首末端节点索引：

from = branch(:, 1); % 支路首端节点 to = branch(:, 2); % 支路末端节点 Constraints = []; % ---- 节点功率平衡 ---- for i = 1:n inIdx = find(to == i); % 流入节点 i 的支路 outIdx = find(from == i); % 从节点 i 流出的支路 Constraints = [Constraints, ... sum(P(inIdx)) - sum(P(outIdx)) + Pg(i) - Pd(i) == 0]; Constraints = [Constraints, ... sum(Q(inIdx)) - sum(Q(outIdx)) + Qg(i) - Qd(i) == 0]; end % ---- 支路电压降方程 ---- for k = 1:nBranch i = from(k); j = to(k); Constraints = [Constraints, ... U(j) == U(i) - 2*(R(k)*P(k) + X(k)*Q(k)) + (R(k)^2 + X(k)^2)*L(k)]; end % ---- 二阶锥约束：P^2 + Q^2 <= U_i * L_k ---- for k = 1:nBranch i = from(k); Constraints = [Constraints, ... norm([2*P(k); 2*Q(k); U(i) - L(k)]) <= U(i) + L(k)]; end % ---- 根节点电压、电压上下限 ---- U0 = 1.0; % 变电站出口 1.0 p.u. Constraints = [Constraints, U(1) == U0]; % 根节点电压固定 Constraints = [Constraints, 0.95^2 <= U <= 1.05^2]; % ---- DG 出力约束 ---- Constraints = [Constraints, 0 <= Pg(dgIdx) <= PgMax]; Constraints = [Constraints, -0.31*Pg(dgIdx) <= Qg(dgIdx) <= 0.31*Pg(dgIdx)];

第一段功率平衡方程的含义是：流入节点的功率加上本地注入功率，等于流出功率加上本地负荷。Pd、Qd是各节点负荷，单位也要折成 p.u.，否则电压降方程里的R*P量纲全乱。

第二段电压降方程直接对应 DistFlow 模型，(R^2 + X^2)*L(k)这一项在有些教程里会被省略，认为数值太小，但严谨的复算会告诉你它影响第 4 章的锥间隙验证，别省。

第三段锥约束是最关键也是最好抄错的部分。norm形式等价于P^2 + Q^2 <= U_i * L_k，它把旋转锥变成了 YALMIP 能直接识别的标准二阶锥。注意左边必须用norm([...])包起来，右边是标量，不能写成P(k)^2 + Q(k)^2 <= U(i)*L(k)，那样 YALMIP 会当成多项式非凸约束拒绝。

最后一组 DG 无功约束用了功率因数 ±0.95 对应的比值关系（0.31 约等于 tan(acos(0.95))）。实际数据里上下限可以按你的算例改，但结构不要动。

3.4 求解器配置与结果输出：从解状态到电压分布图

求解器配置直接用sdpsettings：

ops = sdpsettings('solver', 'mosek', ... 'verbose', 2, ... 'savesolveroutput', 1); sol = optimize(Constraints, Objective, ops); if sol.problem == 0 U_val = value(U); P_val = value(P); Q_val = value(Q); L_val = value(L); netLoss = R' * L_val; % 最优网损 fprintf('总网损: %.4f kW\n', netLoss * baseMVA * 1000); else sol.info end

sol.problem == 0表示求解成功。如果有求解器装好了但 YALMIP 识别不到，ops里显式指定solver是最快的排查办法。免费环境首选sedumi，有授权的话mosek速度优势明显，尤其你要算 123 节点以上规模时，mosek 的内点法对二阶锥问题处理得比 sedumi 稳得多。verbose=2会打印每次迭代的残差信息，第一次跑建议开着，能直观看到求解器卡在哪一步。

提示：MATLAB 装好后，第一件事是把 YALMIP 和求解器目录addpath进搜索路径，再跑一次yalmiptest确认能识别。否则程序还没开始建模就死在求解器检测上。

结果输出阶段，电压分布曲线是检查重点。sqrt(U_val)就是各节点电压幅值 p.u.，画plot(1:33, sqrt(U_val))后，看曲线有没有超过 0.95~1.05 的平线。最优网损要在结果里打印出来，它和第 4 章的回代验证配合，能交叉确认松弛精度。

4. 避坑指南：二阶锥求解在 MATLAB 里最常见的五个翻车现场

4.1 求解器报 Numerical issues：先查数据尺度而不是怀疑模型

现象：optimize返回Numerical problems，解出来的电压是几千甚至几万，sol.problem不是 0。

原因：几乎都是数据尺度问题。线路阻抗用欧姆、负荷用 kW，而电压平方U用 p.u.，两者同时出现在方程里，数值差四五个数量级。内点法对这种病态矩阵非常敏感，迭代直接发散。

解决：所有参数量纲先折成 p.u.。功率基准baseMVA = 10，电压基准baseKV = 12.66，阻抗基准Zbase = baseKV^2 / baseMVA（欧姆），负荷功率除以baseMVA再除以 1000（因为原单位是 kW）。我一般会在数据加载后加一条断言：assert(max(U_val) < 10 && min(U_val) > 0.01)，用电压数量级来判断前期处理有没有做完。

4.2 松弛不精确：锥间隙怎么查、惩罚项怎么加

现象：目标函数值低得不对劲，比如网损只有零点几 kW，但把解回带到普通潮流程序一算，电压全对不上。

原因：锥约束被「松」掉了。当目标函数不含网损项，或者网损项权重很小时，求解器没有动力把解推到锥边界上，最优解落在锥体内部，此时它满足二阶锥不等式，但不是真实潮流的解。

解决：先量化间隙，再决定要不要修。用求解结果算：

U_branch_val = U_val(from); % 支路首端电压平方 gap = max((P_val.^2 + Q_val.^2) ./ (U_branch_val .* L_val));

gap在 1e-6 量级说明松弛是紧的，大于 1e-3 就要处理。最常见补救是目标里显式加入网损项R'*L，并把它的权重加到不影响原目标量级（比如 1e-3 到 1 之间）；如果还是松，就在锥约束上加一个小的线性化惩罚，但那是下策——改目标权重通常就够。

4.3 电流平方约束写成等式：非凸等式直接毁掉松弛

现象：代码里写了L == (P^2 + Q^2)./U，YALMIP 不报错，但求解时间异常长，或者结果明显不对。

原因：这个等式是非凸的，它把凸的锥约束退化回了非凸曲面。YALMIP 对等式约束默认不会帮你检查凸性，遇到这种非凸等式，内部要么转成非线性规划交给 fmincon，要么硬着头皮内点法，结果都很差。

解决：永远不要写L == (P^2 + Q^2)./U。正确写法只有两种：锥不等式norm([2*P; 2*Q; U - L]) <= U + L，或者用旋转锥函数cone([P; Q], sqrt(U*L))。前者对 YALMIP 各版本兼容性最好，也是这份代码采用的写法。想验证等式是否成立，不靠约束，靠的是第 4.2 节的gap指标。

4.4 DG 出力加大后局部电压反升：缺了功率因数与变压器约束

现象：在 18 号等远端节点接入 DG 后，电压曲线出现明显的「鼓包」，局部电压逼近 1.08 p.u.，超过上限。

原因：DG 有功注入改变了原来的单向潮流特征，功率从接入点向根节点倒送，远端电压被抬起来。如果只约束 DG 有功上限、不管无功范围，DG 在高功率因数运行时可能吸收或注入大量无功，进一步加剧电压偏差。

解决：给 DG 加无功上下限，按功率因数约束写，比如 ±0.31 倍有功（对应功率因数 0.95）。同时给变电站根节点注入功率加上限sqrt(P_sub^2 + Q_sub^2) <= S_sub_max，这个约束本身也是二阶锥，正好和模型结构统一。加了这两个约束后，DG 出力增加时模型会在「多发有功」和「电压越限」之间找到更合理的折中。

4.5 YALMIP 识别不到求解器：版本与路径碰出来的坑

现象：代码原封不动，换一台电脑就报No suitable solver for model，或者 mosek 直接提示找不到 license。

原因：YALMIP 的求解器检测依赖 MATLAB 搜索路径和求解器自身的 mex 文件版本。mosek 的 mex 要和 MATLAB 版本严格匹配，跨版本安装经常出现「装了等于没装」的假象；sedumi 则容易出现 addpath 了但没执行install_sedumi的情况。

解决：装好 MATLAB 后，先把 YALMIP 目录、sedumi 目录都加进路径，然后运行yalmiptest，看输出里哪些 solver 显示Success。这一步能过滤掉大部分「玄学」问题。我自己的习惯是免费环境固定用 sedumi，并在sdpsettings里显式写'solver', 'sedumi'，不依赖 YALMIP 自动选择；正式跑大规模算例时再切 mosek，并在调用前先检查exist('mosekdiag','file')判断 license 是否可用。

5. 验证松弛精度与两个顺手扩展：你的结果能不能让别人信服

跑出最优网损只是第一步，能让别人信服才算完成。我每次跑完 SOCP 都会强制走一遍三件套验证流程。第一件是算锥松弛间隙，方法就是 4.2 节那段代码，gap小于 1e-6 才算「紧」；第二件是把求得的节点注入功率回代到牛拉法潮流（用 matpower 的runpf或自写潮流都行），对比每个节点电压幅值，偏差在 1e-4 p.u. 以内可以接受；第三件看约束残差，U的值域是否落在 0.95~1.05 之间，DG 出力有没有超限。这三项过了，结果才算真正落地。

这个验证流程顺手还能扩展出两个实用玩法。第一个是网络重构：把联络开关用binvar建模，锥约束和电压约束里乘上开关状态变量，SOCP 就变成 MISOCP，mosek 能直接求解，目标函数里加开关动作惩罚项就能限制重构次数。第二个是多时段储能调度：把U、P、Q扩展成带时间维的矩阵，储能节点加一个 SOC 递推等式E(t+1) == E(t) + eta*P_cha(t) - P_dis(t)/eta，就变成了典型的日前优化调度模型。这两步扩展之后，你手里的这份代码就不再只是一个静态算例，而是一个能支撑硕士论文核心章节的框架。

从那以后，我每次拿到一套二阶锥最优潮流代码，第一件事都不是看目标函数写得多漂亮，而是先算锥间隙、再回代牛拉法——这两步过了才敢拿结果去写报告。这个习惯帮我挡掉了好几次把错误的松弛解当成最优解的翻车事故，也希望帮到你。

本文还有配套的精品资源，点击获取