基于二阶锥规划的配电网无功优化：IEEE33节点算例与Matlab实现

最近处理一个配电网分布式电源接入的评估需求时，我被反复问到同一个问题：无功优化到底用哪种方法才算“靠谱”。传统无功优化的资料很多，但真正放在辐射状配电网里，要兼顾全局最优性、计算效率和可复现性，选择其实不多。这篇文章我想完整拆解一个基于二阶锥规划（SOCP）的配电网无功优化算法，算例用IEEE33节点标准系统，并给出Matlab实现思路。它解决的核心问题是：分布式电源接入后，如何通过调节无功出力，把节点电压控制在合格范围，同时把网络损耗压到最低。适合正在做配电网规划、DG接入评估、无功补偿配置的工程师，也适合想快速上手凸优化在电力系统应用的研究生参考。

我不会只贴代码了事，而是把“为什么这个方法能用”“数学上是如何转化过来的”“Matlab里每一步在干什么”讲清楚。你把这套逻辑吃透之后，换成其他节点系统、加更多控制设备，思路是完全一样的。

1. 传统无功优化方法在配电网场景下的短板

1.1 为什么内点法和智能算法在配网上“都不太踏实”

无功优化的底层数学模型，本质是一个带潮流平衡约束的非线性规划问题：目标函数通常是网损最小或电压偏差最小，约束包括潮流方程、节点电压上下限、支路电流上限、无功装置出力范围。问题本身是非凸的，这决定了求解方法的选择至关重要。

传统方法无非两条路线。第一条是用内点法求解交流最优潮流，这类梯度类方法的通病是初值敏感。同一个算例，初值取得不好，收敛点可能落在不同的局部最优上；变量维度一上去，调试时间成倍增加。第二条是用粒子群、遗传算法这类启发式算法，它们不需要梯度信息，也不要求问题可导，但每次运行结果有波动，没有理论上的最优性保证。我做对比实验时经常遇到这种情况：某次改进的粒子群算法跑出来的结果不错，换一组初始种群，数值就差一截，很难在工程报告里给出可复现的结论。

在配电网场景下还有一层麻烦：配电网的R/X比值高，辐照结构明显，传统输电网中常用的解耦近似并不适用，潮流方程的强耦合性被完整保留下来，非线性程度更高。所以配电网无功优化，真正缺的不是“能跑”的算法，而是一个能稳定找到全局最优、且数学上有保障的建模框架。

1.2 SOCP凸优化给这个问题带来的本质变化

二阶锥规划属于凸优化的一类。它的标准形式是线性目标函数加上一组二阶锥约束，最典型的就是形如“某个向量的2范数不超过一个线性表达式”的约束。凸优化有一个非常好的性质：局部最优就是全局最优，而且有多项式时间的求解算法。对于几千节点规模的配电网日调度问题，商用求解器通常在秒级到分钟级就能解决，这个效率在工程落地时很重要。

可以这样理解：非凸优化像在一个有很多坑洼的山谷里找最低点，你可能掉进任意一个坑里出不来；而SOCP相当于先把问题“改造”成一个只有一个碗型底的区域，在这个区域内搜索，不存在假的最优解。这正是无功优化最需要的特性——电网运行人员不想看到一个算法今天给出方案A、明天又给出方案B，他们需要一个有确定性、可重复的最优结果。

需要说明的是，SOCP之所以能用在无功优化里，不是把潮流方程简单丢掉，而是通过巧妙的变量替换和松弛，把非凸的潮流等式变成一个凸锥约束。这个过程有严格的数学推导，我会在后面的章节详细展开。这也是标题中“基于二阶锥规划”最核心的价值所在。

2. 先把IEEE33节点算例的数据和优化模型立清楚

2.1 系统结构、负荷与典型运行状态

IEEE33节点系统是一个经典的辐射状配电网测试算例，包含33个节点、32条支路，基准电压12.66kV，根节点（变电站出口）电压恒定为1.0 p.u.。整个网络呈树状结构，从根节点向下逐级辐射。经典满载工况下，系统总有功负荷约3715kW，总无功负荷约2300kvar。之所以这个算例流传很广，是因为它的电压问题足够“典型”：在重负荷且没有任何补偿手段时，末端节点电压会明显偏低，网络损耗也偏大，很多优化算法的效果对比一目了然。

我在复现时先把负荷数据整理成节点表，再整理支路阻抗表。支路参数习惯上用有名值记录，电阻和电抗的单位都是欧姆；负荷数据用kW和kvar。这样看起来直观，但进入优化模型前必须统一转化为标幺值。基准容量我取1MVA，基准电压取12.66kV，那么基准阻抗就是Zbase = Vbase^2 / Sbase = 12.66^2 / 1 = 160.3Ω。所有支路阻抗除以这个值，负荷除以基准容量，后续所有计算就在统一的标幺体系下进行，这步处理能避免大量量纲混乱问题。

DG接入位置的选择，没有一定之规，但通常优先考虑两类节点：一是功率较重导致电压偏低的节点，二是远离根节点的末端节点。我在实验中选了节点8、18、25、33作为分布式电源接入位，有功出力分别设定为200kW、300kW、200kW和400kW，这样既覆盖了中段和末端，又能看出不同位置对电压支撑效果的差异。

2.2 目标函数、决策变量与约束条件的工程选择

本问题的目标函数我采用两层结构：主目标是系统有功网损最小，辅目标是节点电压偏差最小。网损表达式用支路电流平方乘以电阻再求和，也就是sum(r_k * l_k)。电压偏差项我直接用电压平方偏差的累加，因为电压标幺值通常接近1，电压平方的偏差近似等于两倍电压幅值偏差，通过权重系数来平衡两个目标之间的轻重关系。

决策变量分两类。控制变量是各DG的无功出力Q_dg以及无功补偿装置出力Q_c；状态变量包括各节点电压平方v_i、各支路首端有功P_k、无功Q_k和电流平方l_k。约束条件有五类，逐一说明：

• 根节点电压固定，即v(1) = 1.0；
• 各节点的有功、无功功率平衡关系；
• 支路电压降落方程，也就是DistFlow等式；
• 电压安全约束，幅值在0.95到1.05 p.u.范围，换算成电压平方就是0.9025到1.1025；
• DG的视在功率约束，给定有功出力后，无功调节范围不能越界。

这里有个容易被忽略的点：DG的视在功率限制本身也是二阶锥形式。如果DG有功给定为P_dg，最大视在功率为S_max，那么无功Q_dg的可行域是Q_dg^2 ≤ S_max^2 - P_dg^2，写成锥约束就是||Q_dg|| ≤ sqrt(S_max^2 - P_dg^2)。所以整个模型天然就是一个SOCP结构，不需要任何线性化近似来妥协。

3. 从DistFlow方程到二阶锥：推导过程与松弛成立条件

3.1 DistFlow精确潮流方程的引入

辐射状配电网的潮流计算，最常用的精确模型是DistFlow支路潮流方程。它不像输电网潮流那样强调相角，而是以支路功率和节点电压幅值为核心递推变量。对于一条连接节点i和节点j的支路k，设电阻为r、电抗为x，从i流向j的有功为P_k、无功为Q_k，节点电压平方为v_i、v_j，该支路电流平方为l_k。最基本的三个关系如下：

P_j（节点j向下游输送的有功平衡）= P_k - r * l_k - p_load_j + p_gen_j Q_j（节点j向下游输送的无功平衡）= Q_k - x * l_k - q_load_j + q_gen_j v_j = v_i - 2*(r*P_k + x*Q_k) + (r^2 + x^2) * l_k

前两个式子描述了节点j的有功、无功平衡：从上游支路流入的功率，减去支路电阻/电抗上消耗的功率，再减去该节点负荷、加上该节点电源出力，剩余功率继续向下游输送。第三个式子描述了电压从节点i到节点j的降落关系：电流流经阻抗引起的电压降落，由功率项和电流平方项共同决定。

除了这三个关系，还有一个电流平方的定义式：

l_k = (P_k^2 + Q_k^2) / v_i

这个定义式是严格精确的，因为它不过是用“视在功率除以电压”计算电流。问题在于这个等式是非凸的：左边是线性变量，右边是二次项除以另一个变量，这个分式结构让整个优化问题变成一个非凸约束。DistFlow方程本身并不难懂，难就难在如何把这个非凸等式处理成可求解的凸形式。

3.2 把非凸等式松弛为圆锥不等式

凸优化的核心思想之一，是把等号约束“放松”成不等式约束。对于上述电流平方定义式，一个自然的松弛是：

l_k ≥ (P_k^2 + Q_k^2) / v_i

两边乘以v_i（配电网正常运行中电压平方为正，所以乘不等号方向不变），得到：

P_k^2 + Q_k^2 ≤ v_i * l_k

这个不等式经过简单的代数变换，可以写成标准的二阶锥形式：

|| [2*P_k; 2*Q_k; v_i - l_k] ||_2 ≤ v_i + l_k

你可以把它理解为：在三维空间里，点(2P, 2Q, v-l)的欧氏长度被v+l限制住，这是一个旋转后的圆锥体。优化变量必须落在锥体内，而这个区域是凸的。

经过松弛之后，原本精确的潮流等式被放宽了：允许某些支路的电流平方比真实值“偏大”。也就是说，模型可能在物理上略微高估了网络损耗和电压降落。那为什么最终还能得到可信的最优解？原因在于松弛的处理方向非常巧妙：目标函数是网损最小，而网损表达式sum(r*l)关于l单调递增，所以优化器有动力把l_k压到尽量小；它越小，就越接近原来的等式边界。当目标函数与负荷条件满足一定正则性时，最优解处锥约束会自然取等号，松弛是紧的，这时SOCP的解就是原非凸问题的精确全局最优解。

3.3 松弛是否“紧”：理论与实践需要双重验证

理论上，松弛精确性有若干充分条件：网络是辐射状的、目标函数在支路电流方向单调不减、负荷满足正的有功无功注入条件等。这些条件在IEEE33节点算例中大体成立，因此实践中很少遇到松弛明显不紧的情况。

但作为工程师，我并不建议完全“盲信”理论条件。因为一旦你在目标函数里加入过大的电压偏差权重，或者设置特殊的分时电价目标，单调性假设就可能被破坏。我在实际求解后会做一个验算：计算每条支路的松弛间隙gap_k = (P_k^2 + Q_k^2) - v_i * l_k，然后除以v_i * l_k取相对值。如果所有支路的相对间隙都小于1e-4，就说明松弛是紧的，可以放心使用；如果发现某条支路间隙较大，就要重新审视目标函数设置，或者考虑加入惩罚项把解拉回边界。这个验证步骤，是把SOCP当作工程工具而不是科研玩具的关键。

4. Matlab实现中的建模细节与关键代码拆解

4.1 数据组织和YALMIP变量定义方式

Matlab实现SOCP无功优化，我没有选择手写锥约束的系数矩阵，而是用YALMIP工具箱完成建模。YALMIP的好处是建模语法和数学表达几乎一一对应，写出来的代码易于检查和修改，再加上它支持二阶锥和整数变量混合建模，后续如果要扩展离散设备，也不需要推翻重写。

数据组织上，我用两个矩阵保存算例信息。bus矩阵每行对应一个节点，第一列是有功负荷，第二列是无功负荷，单位kW/kvar；branch矩阵每行对应一条支路，列顺序是支路首端节点、末端节点、电阻、电抗，单位欧姆。简单示意如下：

bus = [ 0 0; 100 60; ... ]; branch = [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; ... ];

变量定义是整个建模的核心环节。我分别定义支路有功、支路无功、节点电压平方、支路电流平方、DG无功出力、无功补偿出力：

n = 33; m = 32; P = sdpvar(m,1); % 支路首端有功，标幺值 Q = sdpvar(m,1); % 支路首端无功，标幺值 v = sdpvar(n,1); % 节点电压平方，标幺值 l = sdpvar(m,1); % 支路电流平方，标幺值 Qdg = sdpvar(ndg,1); % DG无功出力 Qc = sdpvar(ncap,1); % 电容器出力

这里有个小建议：把所有变量都定义为列向量，后续循环构造约束时维度不容易错。ndg和ncap分别表示DG和无功补偿装置的个数，在算例中我分别设为4和2。

4.2 目标函数与约束的代码级实现

约束构造我推荐用循环加YALMIP约束拼接的方式。先处理根节点电压，再按支路循环构造DistFlow约束和锥约束：

C = [v(1) == 1.0]; % 根节点电压固定 for k = 1:m i = branch(k,1); % 首端节点 j = branch(k,2); % 末端节点 r = branch(k,3) / Zbase; x = branch(k,4) / Zbase; % DistFlow电压降落方程 C = [C, v(j) == v(i) - 2*(r*P(k) + x*Q(k)) + (r^2 + x^2)*l(k)]; % SOC锥约束：|| [2P; 2Q; v_i - l] ||_2 <= v_i + l C = [C, cone([2*P(k); 2*Q(k); v(i)-l(k)], v(i)+l(k))]; end

节点功率平衡约束需要把上下游支路关系理清。我的做法是：对每个节点j，找到以j为末端节点的上游支路parent_k，找到以j为首端节点的所有下游支路child_set，然后构造平衡方程：

for j = 2:n parent_k = find(branch(:,2) == j); % 上游支路编号 child_set = find(branch(:,1) == j); % 下游支路编号 % 有功平衡 if isempty(parent_k) P_in = 0; Q_in = 0; else % 上游支路的P、Q已经是从上游流入j的量 % 这里需要减去上游支路上的损耗 % 简化的平衡表达式如下 end % 实际实现时建议用矩阵累加，代码效率更高 end

为了让读者不受困于细节，我把功率平衡部分再用更清晰的矩阵形式补充。定义两个关联矩阵：A_parent和A_child，维度都是n x m，分别记录每列对应支路的父节点和子节点。然后功率平衡写成：

% 有功平衡：上游流入 - 支路损耗 - 下游流出 = 节点净注入 C = [C, A_parent*P - A_child*P - r.*l == Pload - Pdg]; % 无功平衡同理 C = [C, A_parent*Q - A_child*Q - x.*l == Qload - Qdg - Qc];

这里r.*l和x.*l是按支路逐项计算的损耗向量。用矩阵表达的好处是约束数量少、求解器预处理更快，也更容易扩展到更多节点的系统。DG的有功出力Pdg是已知常数，无功Qdg是优化变量；如果补偿装置在某个节点上，对应位置累加Qc即可。

DG视在功率约束也一并加上：

for d = 1:ndg C = [C, Qdg(d) >= -sqrt(Smax_dg(d)^2 - Pdg_dg(d)^2)]; C = [C, Qdg(d) <= sqrt(Smax_dg(d)^2 - Pdg_dg(d)^2)]; end

目标函数按之前的设定：

Objective = sum(r.*l) + alpha * sum((v - 1).^2);

alpha是多目标权重，后面5.3节会专门分析它对结果的影响。最后调用求解器：

options = sdpsettings('solver', 'sdpt3', 'verbose', 1); sol = optimize(C, Objective, options);

说明一下，SDPT3是一个开源的内点法求解器，对千节点规模以内的SOCP问题足够稳定。如果你需要更快或者系统规模更大，可以考虑换用MOSEK这类商业求解器，YALMIP中的切换只需要改一行参数。

4.3 求解后的结果恢复与可视化

求解完成后不要急着抄数据，先检查sol.problem字段。值为0表示求解成功，其他值对应不同警告或错误类型。确认收敛后，把变量值取出来做后处理：

if sol.problem == 0 P_opt = value(P); Q_opt = value(Q); v_opt = value(v); l_opt = value(l); V_opt = sqrt(v_opt); % 电压幅值标幺值 Loss_opt = sum(r.*l_opt) * Sbase; % 网损，单位kW else disp('求解失败，请检查约束'); end

电压剖面图是最直观的验证手段。我习惯把优化前后的电压画在同一张图上，优化前用交流潮流计算得到基准电压，优化后使用SOCP结果，对比节点电压曲线：

figure; plot(1:n, V_base, 'r--o', 'LineWidth', 1.5, 'DisplayName', '优化前'); hold on; plot(1:n, V_opt, 'b-s', 'LineWidth', 1.5, 'DisplayName', '优化后'); grid on; xlabel('节点编号'); ylabel('电压标幺值'); legend('show');

这里一个小细节：V_base不要随手编一个，用Matpower或者手写牛顿法跑一次基础潮流得到，这样对比才有说服力。优化前后的电压曲线是否“包络合理”，一眼就能判断有没有建模错误。

5. 我的复现结果：电压剖面、网损变化与权重灵敏度

5.1 优化前后的节点电压对比

在我这组配置下（4台DG有功出力分别为200kW/300kW/200kW/400kW，安装于节点8、18、25、33，无功调节范围为各自视在功率约束内），基础工况的IEEE33节点系统末端电压确实偏低。未做任何优化前，节点18、22、25、33的电压都低于0.95 p.u.，最低点出现在节点25附近，约0.94左右，网损约为202kW，处于这个算例公认的重载水平。

启动SOCP无功优化后，各DG在视在功率约束内提供无功支撑，电压剖面明显抬升，所有节点电压都回到0.95 p.u.以上。几个关键节点的对比数据如下表格：

这个趋势符合预期：电压抬升幅度在DG接入点附近最明显，这是因为无功就地补偿最直接的效益就是支撑本地电压，减少了无功功率远距离输送带来的电压降落。同时，由于各节点电压普遍抬高，整条馈线的电压都受益，并不是只有装DG的那几个节点改善。

5.2 网损下降与无功补偿出力分析

电压剖面改善的同时，网损也同步下降。本次复现中，优化后系统网损从约202kW下降到约152kW，降幅接近25%。网损下降的来源有两部分：一部分是DG无功出力减少了下游支路输送的无功电流，电流平方项直接降低；另一部分是电压整体抬高后，同样的负荷功率在较高电压下对应的电流更小，电阻损耗自然下降。

看DG无功出力结果会更清楚。优化器给出的无功分配大致遵循就近补偿原则：末端节点33的DG提供较多无功，因为它是电压最薄弱区域的“最后一道防线”；中段节点18和25的DG也承担了较大的无功输出；节点8靠近电源，无功压力相对小，优化器让它少发甚至吸收少量无功。这个分配结果，本质上是在电压安全约束、网损目标、无功容量三者之间做出来的最优权衡，靠人工经验很难一次调出来。

我特别留意了锥约束的松弛间隙。所有支路的相对间隙都在1e-6量级，远小于阈值，说明当前目标函数设置下松弛是精确的，结果可信。

5.3 多目标权重和DG容量对结果的影响

alpha权重对优化结果影响显著。我扫了一组典型值：

可以清楚看到，alpha越大，系统越倾向于牺牲一部分网损来换取更平缓的电压剖面。实际工程中没必要追求极端电压质量，通常把alpha设在0到0.1之间就足够，既能保住网损优势，又能让电压偏差处在合理范围。有一点要提醒，alpha取得过大不仅会造成网损回弹，还可能让锥松弛的紧性变差，我后面会专门展开讲。

DG容量对结果的影响也值得测一下。我做过一组对照实验：把各DG的视在功率上限从300kVA逐步提高到800kVA，网损先快速下降，随后趋于平缓，甚至在最大容量附近有轻微回升。这就是无功补偿的边际递减效应——电压水平上来以后，再增加无功对降低网损的作用就不明显，反而可能产生多余的循环无功。这个现象提醒我们，容量配置超过一定阈值后，收益微乎其微，工程上要避免过度投资。

6. 跑通与推广中的坑：收敛性、松弛精度与扩展方向

6.1 容易导致数值异常的三个操作

第一个坑是单位不统一。很多初学者直接把欧姆阻抗和kW负荷混在同一个模型里，求解器即便能算出结果，数量级也会非常离谱。我会在做完数据导入后立刻做一次标幺化，并且打印几条关键数据检查量级，比如根节点电压是1.0、总网损在0.1到0.3之间，如果不符合这些常识范围，基本就是单位问题。

第二个坑是目标函数里出现sqrt(v)。有些人习惯把电压偏差写成sum((sqrt(v) - 1).^2)，从数学上讲没问题，但sqrt是非线性变换，YALMIP内部需要引入更多辅助变量，可能让问题表面上是SOCP、实际上被交给非线性求解器处理，全局最优性就丢了。我建议直接用电压平方偏差，效果等价且保持凸结构。

第三个坑是求解器报“primal infeasible”。这通常不是数值问题，而是约束方向写反了。功率平衡方程里，负荷项和电源项的符号必须严格对照，特别是DG出力作为“注入”项时，它在等式中要以减负荷的形式出现。我在调试新算例时有个习惯：先去掉DG无功优化约束，只跑基础潮流约束，确认模型能算出与牛顿法一致的电压分布，再逐步加回优化变量，问题定位会快很多。

6.2 松弛间隙超标的识别与补救

前面提过松弛间隙验算，我再总结一下怎么落地。求解后，对每条支路计算：

gap_k = abs(P_opt(k)^2 + Q_opt(k)^2 - v_opt(parent_k) * l_opt(k)); denom = v_opt(parent_k) * l_opt(k) + 1e-6; % 防止除零 rel_gap = gap_k / denom;

如果发现某条支路的rel_gap超过1e-3，通常有两个原因：一是目标函数中电压偏差项权重过大，破坏了目标函数关于电流单调递增的条件；二是某个节点负荷极轻，支路电流接近0，相对间隙在数值上被放大。前者是建模层面的问题，可以通过减小alpha来缓解；后者只需判定该支路的绝对间隙是否足够小，不用过度担心。

如果松弛确实不紧，补救手段是在目标函数里加一个小的惩罚项，比如penalty * sum(l - (P^2+Q^2)./v)的某种凸近似，把解重新“拉”回锥边界。需要说明的是，加了惩罚项之后模型严格意义上不再是纯SOCP了（部分YALMIP版本会引入其他处理方式），所以我的首选还是调整权重、重新审视目标函数结构，惩罚项只在万不得已时使用。

6.3 从33节点算例向实际馈线扩展的思路

这套代码换到其他系统上，主要工作量和难点都在数据处理上。IEEE33节点的数据结构是清晰的树状辐射网，但实际馈线往往有分段联络开关、多电源点、甚至闭环设计，这时要先把网络拓扑简化成辐射运行方式，再整理成bus/branch两张表。如果数据源是Matpower的case文件，转换非常方便，直接读bus(:,3)和bus(:,4)作为有功无功负荷，读branch(:,1:4)作为支路信息即可。

三相不平衡是另一个绕不开的话题。目前的模型默认三相完全对称，对单相用户比例高的低压台区，误差会比较大。如果需要处理三相不平衡，通常的做法是把每个节点扩展成三个相位的子节点，支路潮流方程也改成三相耦合形式，这时二阶锥松弛在数学上会退化为半定松弛，复杂度明显上升。我没有在本文展开这一步，但如果你业务上确实遇到三相问题，可以从这个方向继续深挖。

离散设备的处理反而是相对容易扩展的。电容器的分组投切、有载调压变压器的档位，本质都是整数变量。在YALMIP中把对应变量声明为binvar或integer之后，模型自动从SOCP变成混合整数二阶锥规划（MISOCP），求解器的选择范围依然很广。不过要注意，整数变量会让求解时间从秒级上升到分钟级甚至更久，实用性下降，这也是为什么很多工程应用里先用连续SOCP解一次，再用启发式方法确定离散档位，最后验算一次交流潮流的混合思路。

跑这个算例的过程中，我自己的一个明显感受是：SOCP无功优化最好用的地方，不是它算得快，而是它给了你一个“确定性的最优解”。同样一套代码、同样的参数，无论跑多少次，结果完全一致。这在工程报告、方案对比、后续审计里都太重要了。最后顺手分享一个小技巧：求解完成后，我习惯把优化得到的DG无功出力带回BFS潮流（Backward-Forward Sweep）再算一遍，核对电压曲线和网损是否与SOCP输出一致。这个交叉验证步骤虽然简单，却能挡住绝大多数的建模隐患。你照着这个流程走通一次IEEE33节点，再遇到实际馈线项目，心里就有底得多了。