二阶锥规划配电网无功优化：DistFlow模型与33节点Matlab实现

拿到这个“基于二阶锥规划的配电网无功优化算法（33节点系统）”的项目时，我第一反应不是去搜代码，而是先问自己一个实际问题：配电网的无功优化，到底想解决什么？

我上一份工作正好接过一个类似的配网降损分析任务。某片区的馈线末端一到用电高峰就出现低电压，分布式光伏接入后又出现午间电压偏高，传统做法是调度员手动投切电容器，再不行就调主变分接头。可一碰上负荷波动和光伏出力的双重不确定性，手动调根本跟不上节奏。无功能不能自动算、自动调，就成了一个很现实的工程问题。

这套算法做的事情，就是把这个“自动算”变成一套数学优化模型：以网损最小为目标，以潮流方程和运行电压为约束，用二阶锥规划（SOCP）这种凸优化方法去求解，最后在一套公开的33节点配电测试系统上验证效果。整个项目落地的路径是：模型推导 → Matlab建模 → 求解器求解 → 结果分析。

这篇文章我会把模型原理、代码结构和实操过程完整拆开讲。适合两类人看：一类是刚接触配电网优化的在读研究生，想知道SOCP为什么这几年这么火；另一类是做配网自动化、想快速拿一套可靠算法做验证的工程师。看完你至少能独立搭出一版可运行的配电网无功优化Demo，并且知道遇到不收敛、结果不合理时该从哪里排查。

1. 问题背景与方案选型：为什么无功优化偏偏选中SOCP

1.1 配电网无功问题到底长什么样

无功功率不直接做功，但电网离不开它。它是维持电压水平的关键“支撑力”。配电网里负荷大多数是感性负荷，需要消耗无功；分布式光伏逆变器、风电变流器虽然可以发无功，但很多现场早期配置都只按“有功最大化”模式运行，等于白白浪费了调节能力。

无功不足的直接后果是电压跌落，无功过剩则电压抬高。而电压偏移又会反过来增加线路损耗、降低设备寿命，甚至触发保护动作。我在项目现场见过最典型的场景：一条长约十几公里的辐射馈线，末端变压器低压侧电压低到0.90pu，空调压缩机直接启动不了；等午间光伏大发，同一片区电压又能冲到1.08pu以上。两个方向都越限，靠人工已经很难处理。

于是“无功优化”就变成了一个数学问题：在满足潮流约束、电压约束、无功源出力约束的前提下，找到一个最优的无功出力策略，让目标函数达到最小。目标函数多数情况下是网损最小，因为网损直接对应经济性；也有不少场景再加上电压偏差最小的多目标版本。

1.2 为什么淘汰启发式算法，换赛道选凸优化

早些年这类问题的主流解法是粒子群算法、遗传算法、模拟退火等启发式算法，包括现在期刊上依然大量存在。它们的思路很直接：随机生成一批无功出力方案，迭代搜索，保留适应度好的个体。逻辑上没有硬伤，但工程实践中有三个绕不过去的坑：

第一是全局最优性无法证明。启发式算法得到的往往是一个“还不错的解”，没人能保证它是最好的。第二是参数敏感。种群规模、交叉概率、变异概率、惯性权重，每一项都要针对具体网络和负荷水平重新调，换一个算例可能就废了。第三是计算时间不可控。33节点还好，一旦扩展到几百上千节点的实际馈线，单次潮流计算叠加反复迭代，现场可接受的时间窗口往往撑不住。

二阶锥规划（SOCP）是另一条路。它的核心优势在于：约束集合是一个凸锥，目标函数如果是线性的或者凸的，整个问题就是一个凸优化问题。凸优化有一个让人安心的性质——局部最优就是全局最优。而且底层用的是内点法，求解效率在目前的主流算例规模下非常可观。这几年配电网领域大量换用SOCP，本质上就是大家发现：与其花时间调遗传算法的参数，不如把模型转化成凸问题，让成熟求解器去算，又快又有理论保证。

不过要强调一点，SOCP并不是凭空出现的，它依赖一个能被“松弛成凸”的潮流模型。这就引出配电网最常用的DistFlow模型，也是整个算法的最核心部分。

2. 核心模型：DistFlow潮流方程与二阶锥松弛

2.1 DistFlow方程：从辐射状馈线写出的递推关系

配电网跟输电网的最大区别是拓扑结构：输电网通常是环网运行，配电网则多采用辐射状结构，也就是一条馈线从根节点出发，层层分支向末端供电，没有闭环。这个结构给了我们一个简化描述潮流的机会。

DistFlow模型就是针对辐射状配电网的一种准确潮流描述。它不像牛顿-拉夫逊法那样要解一个庞大的节点功率平衡方程组，而是按照支路顺序，用递推关系把父节点和子节点的功率、电压联系起来。

对于一个支路i，设父节点为节点i，子节点为节点i+1，支路电阻为R、电抗为X，那么DistFlow的核心递推关系可以写成：

[ P_{i+1} = P_i - R_i \frac{P_i^2 + Q_i^2}{V_i^2} - p_{i+1}^{load} ]

[ Q_{i+1} = Q_i - X_i \frac{P_i^2 + Q_i^2}{V_i^2} - q_{i+1}^{load} ]

[ V_{i+1}^2 = V_i^2 - 2(R_i P_i + X_i Q_i) + (R_i^2 + X_i^2) \frac{P_i^2 + Q_i^2}{V_i^2} ]

三个式子分别表示有功、无功和电压平方沿馈线的传递关系。这里的( P_i, Q_i )是从支路i首端流向末端的功率，( V_i )是节点i的电压幅值，( p_{i+1}^{load}, q_{i+1}^{load} )是节点i+1的负荷或者分布式电源注入。

这个模型在物理意义上很直观：功率从电源点往负荷点走，每经过一段线路，就有有功和无功的损耗，电压也随之下降。很多文献里也把DistFlow叫BranchFlow模型，英文缩写就是DistFlow或DistFlow方程，指的都是同一族表达式。

但这里藏着一个麻烦：等式右边有一个分式项(\frac{P_i^2 + Q_i^2}{V_i^2})。这个项是非线性的，而且不是凸的。如果直接把它丢给求解器，得到的就是一个非凸优化问题，求解起来跟启发式算法一样会遇到局部最优的困境。

2.2 从非凸到凸：二阶锥松弛到底做了什么

为了把这个非线性约束转成凸约束，标准做法是引入辅助变量。先定义支路电流平方的辅助变量：

[ L_i = \frac{P_i^2 + Q_i^2}{V_i^2} ]

这个等式又等价于：

[ P_i^2 + Q_i^2 = L_i V_i^2 ]

我们把这个约束“放松”成不等式：

[ P_i^2 + Q_i^2 \leq L_i V_i^2 ]

这一步就是整个算法里的“松弛”。看起来只是把等号变成了小于等于号，但其实质是把一个非凸集合扩展成了一个凸集合。配电网无功优化的DistFlow方程就变成了：

[ V_{i+1}^2 = V_i^2 - 2(R_i P_i + X_i Q_i) + (R_i^2 + X_i^2) L_i ]

[ P_{i+1} = P_i - R_i L_i - p_{i+1}^{load} ]

[ Q_{i+1} = Q_i - X_i L_i - q_{i+1}^{load} ]

原本的三个式子仍然保留，只有那个二次等式约束被“放松”成了一个凸锥约束。这就是二阶锥规划的由来。

很多第一次接触SOCP的人会担心：把等号放成小于等于，解出来的结果还是原问题的解吗？答案是：在目标函数合适的前提下，最优解会自然地落在锥边界上。原因非常直观——我们优化的目标是网损最小，而网损恰恰是(\sum R_i L_i)的某种加权和。如果( P_i^2 + Q_i^2 < L_i V_i^2)，说明电流平方被“虚高”了，网损因此被抬高，这样的解不可能最优。所以在网损最小化这个目标的驱动下，优化算法会自动迫使不等式取等号，使松弛是“紧”的，我们就得到了原问题的可行解。

在Matlab的YALMIP或CVX中，二阶锥约束一般写成标准锥形式。一个旋转锥 (P^2 + Q^2 ≤ V L)可以改写成：

[ \left| \begin{matrix} 2P \ 2Q \ V - L \end{matrix} \right|_2 \leq V + L ]

这种形式正是求解器直接支持的二阶锥约束。

2.3 松弛不紧怎么办

“松弛紧”听上去很理想，但实际跑代码时偶尔也会翻车，尤其是目标函数不只包含网损，还加了电压偏差项、无功调节代价项的时候。因为这些额外目标可能跟“拉紧松弛”的方向不一致，最终解仍然满足不等式但不等号严格成立，这时得到的解在物理上解释不通：电流平方比实际要求的值更大，电压和功率流对不上。

我的处理习惯是在求解完成后做一次紧性校验：

[ residual = \max \left| P_i^2 + Q_i^2 - L_i V_i^2 \right| ]

如果残差超过1e-4或者1e-3这个量级，就得想办法。常用的补救方法是在目标函数里加一个很小的惩罚项，例如加(\epsilon \sum L_i)或(\epsilon \sum (P_i^2 + Q_i^2))，把解“拽回”锥边界上。权重要取小，比如1e-4到1e-6，否则会干扰原始目标。

另外还有一种情况是电压初值给得太离谱，导致锥约束在可行域边缘反复震荡。这时候把电压初始化成额定值、支路功率初始化成负荷叠加值，通常就能缓解。

3. Matlab实现全流程：从环境准备到参数导入

3.1 环境搭建：Matlab + YALMIP + 求解器

这套算法落地的标准工具链是Matlab加上一个建模工具，再加一个求解器。我推荐用YALMIP做建模层，因为它对二次锥约束的表达非常友好，语法接近数学公式本身，不需要手动把所有约束都改写成标准锥形式。

求解器方面，商业求解器和免费求解器都可以。免费方案里SDPT3、SeDuMi都能解二阶锥问题，对33节点这个小算例绰绰有余。如果想体验更快的求解速度，可以配置主流商用求解器，配好后在YALMIP里只要一行代码切换。

安装顺序上有一点要注意：先把Matlab装好，然后下载YALMIP，把解压后的文件夹加入Matlab路径，再安装求解器。很多人容易卡在“加入了路径但YALMIP命令找不到”这一步，通常是没执行savepath或者当前路径没有保存。正确做法是在Matlab里运行pathtool，把YALMIP和求解器所在目录都加进去，然后点Save，再运行yalmip和solver_installed检查是否就绪。

3.2 33节点标准算例的参数导入与预处理

33节点系统是配电网优化里最常用的公开测试算例之一。它包含33个节点、32条支路、1个主电源点，负荷基本分布在末端附近，网络呈典型的辐射状结构。基准电压是12.66kV，基准功率通常取10MVA或1MVA，做标幺化时务必统一基准值，不然算出来的电压会整体偏移。

拿到参数后的第一件事，不是写优化模型，而是把支路参数表和数据检查清楚。我的做法是：先读入支路首端、末端节点编号、电阻电抗、负荷有功无功，然后用稀疏矩阵把所有参数组装起来。这样写出来的代码更接近实际工程，尤其是在后续扩展成几百个节点的网络时，手工连支路矩阵会后悔的。

预处理阶段还有个容易忽略的点：如果有联络开关，要先把联络开关处理成断开支路，保证网络是纯辐射状。否则DistFlow模型的“父节点→子节点”递推关系就不成立了。实际代码里我会把支路按“从根节点出发，前推后带”的顺序排一下，或者直接使用节点注入功率平衡方程来建模，后者对支路顺序不敏感，实现起来更省心。

4. 核心代码实现与模块拆解

4.1 决策变量与DistFlow约束构建

下面是基于YALMIP建模的核心代码结构，重点看约束怎么构建：

% 基本参数 nb = 33; % 节点数 nl = 32; % 支路数 % 假设branch_tbl每行: [首端节点, 末端节点, R, X, P_load, Q_load] branch = branch_tbl; R = branch(:, 3); X = branch(:, 4); % 决策变量 V2 = sdpvar(nb, 1, 'full'); % 电压平方 P = sdpvar(nl, 1, 'full'); % 支路有功 Q = sdpvar(nl, 1, 'full'); % 支路无功 L = sdpvar(nl, 1, 'full'); % 支路电流平方 Constraints = []; % 根节点电压固定为额定值 Constraints = [Constraints, V2(1) == 1.0]; % 支路容量约束（可选） MaxS = 0.2; % 标幺 Constraints = [Constraints, P.^2 + Q.^2 <= MaxS^2]; % DistFlow 约束 for k = 1:nl i = branch(k, 1); j = branch(k, 2); % 电压递推 Constraints = [Constraints, ... V2(j) == V2(i) - 2*(R(k)*P(k) + X(k)*Q(k)) + (R(k)^2 + X(k)^2)*L(k)]; % 二阶锥松弛约束: P^2 + Q^2 <= V2(i) * L(k) Constraints = [Constraints, cone([2*P(k); 2*Q(k); V2(i)-L(k)], V2(i)+L(k))]; end % 节点注入功率平衡 for k = 1:nb idx_out = find(branch(:, 1) == k); % 以k为首端的支路 idx_in = find(branch(:, 2) == k); % 以k为末端的支路 P_balance = sum(P(idx_out)) - sum(P(idx_in)); Q_balance = sum(Q(idx_out)) - sum(Q(idx_in)); % 负荷和DG注入 P_inj = P_load(k) - P_gen(k); Q_inj = Q_load(k) - Q_gen(k); Constraints = [Constraints, P_balance == P_inj]; Constraints = [Constraints, Q_balance == Q_inj]; end

这一段代码里，cone([2*P; 2*Q; V2(i)-L(k)], V2(i)+L(k))就是关键表达。YALMIP里cone的第一个参数是向量，第二个参数是标量，表示( |v|_2 \leq t)。这个语法一定要写对，很多新手把cone的写反了，程序直接报维度错误。

另外要注意：节点注入约束用“出支路之和减入支路之和”来表达，这样不需要支路排序，更稳健。我在初学阶段按DistFlow递推法逐步写约束，结果支路顺序一乱就全错了；换成这种节点平衡写法之后，同样的模型可以平滑迁移到其他网络结构。

4.2 目标函数与求解配置

目标函数里最常用的是网损最小：

Objective = sum(R .* L); optimize(Constraints, Objective, ops);

注意这里的L是电流平方，乘以支路电阻R得到了功率损耗。在标幺值体系下，优化结果再乘以基准功率就换算成了kW。

如果是多目标，可以加电压偏差项。例如希望电压尽量靠近额定值，就在目标里加一个二次项：

Objective = sum(R .* L) + 0.5 * sum((V2 - 1.0).^2);

但这个二次项会破坏“目标函数线性+二阶锥”的纯粹性。虽然YALMIP仍会把问题归类为二次锥规划（SOCP/QP）的变体，但选择权重的逻辑要变，尤其要关注2.3节说的松弛紧性问题。

求解配置上，我一般这样设置：

ops = sdpsettings('verbose', 2, 'solver', 'sdpt3', 'savesolveroutput', 1);

这里的solver可以换成你实际安装的求解器名字。需要留意的是，不同求解器对初始点的处理方式不同。SDPT3通常能自己找到可行初始点，但给出一个合理初值往往能加快收敛：

assign(V2, ones(nb, 1)); assign(P, P_init); assign(Q, Q_init); assign(L, zeros(nl, 1));

4.3 结果输出：网损、电压分布与紧性校验

求解完成后，用value()提取变量：

V2_opt = value(V2); P_opt = value(P); Q_opt = value(Q); L_opt = value(L); V_opt = sqrt(V2_opt); loss_opt = sum(R .* L_opt); % 标幺值 loss_kw = loss_opt * Sbase * 1000; % kW % 紧性校验 residual = max(abs(P_opt.^2 + Q_opt.^2 - V2_opt(branch(:,1)).*L_opt)); disp(['最大松弛残差: ', num2str(residual)]);

这步千万别省。残差大说明模型虽然求出了“数学最优”，但物理上可能是伪解。我在实际项目中有一次残差到了0.01量级，找了好久才发现是无功约束方向写反了，导致锥松弛过度宽松。

电压分布可以画出来，直接展示优化前后各个节点的电压曲线。或者把结果存成表格，便于写进汇报。我个人习惯是把数据导出成结构化表格，然后交给后续绘图或报表程序。

5. 33节点系统上的结果分析与扩展方向

5.1 优化前后的网损与电压对比实测

在标准33节点系统的典型负荷数据下，我用这套SOCP流程跑了一版，结果如下表所示：

说明一下，这里的绝对数值依赖基准功率选取、负荷数据和求解器精度，不同实现版本之间会有少量浮动，但趋势是一致的：通过无功优化，最低电压被抬高了0.04pu以上，网损下降20%以上。

这个结果其实很好解释。原先电网里大量无功从电源点长途输送到末端，在线路上产生了额外的电流和发热损耗。无功优化让末端附近的无功源去就近补偿，线路上的无功潮流减少了，网损自然降下来，电压支撑也变好了。

我自己拿到这个结果时，会再检查一遍所有变压器分接头位置和电容器投切状态是否在合理范围。因为有些优化器会“投机取巧”，把某个节点电压推到接近上限来降低网损，这在单一断面下看似最优，但在实际运行中却是不可取的“尖峰策略”。所以更稳妥的做法是在约束里把电压目标限定在0.95到1.05pu之间，而不是一个很宽的可接受区间。

5.2 进阶场景：分布式电源接入后的扩展

很多人拿到这套33节点代码后，会立刻想加分布式光伏或风电。这是一个很自然的扩展，也是当前工程现场的刚需。

接入分布式电源，本质上是给某些节点增加一个可以调节的无功注入源。模型里只需把节点的注入功率改一下：

% 光伏有功出力固定为某值 P_gen(k) = P_pv; % 逆变器无功容量约束 Q_gen(k) = sdpvar(1, 1); Constraints = [Constraints, Q_min <= Q_gen(k) <= Q_max];

注意逆变器的无功容量通常与有功出力有关，比如功率因数限制在0.95，即( |Q| \leq P \cdot \tan(\arccos(0.95)) )。写约束时不要把Q的上限设成固定值，那样会忽略有功出力变化对无功空间的挤压。

还有一层的应用是“多时段协同优化”：把一天分成24个断面，每个断面单独算一套SOCP，然后把相邻时段的调节量作为额外约束，比如电容器投切次数限制、有载调压变压器档位动作次数限制。这种扩展会让模型从“静态优化”变成“动态优化”，工程量会大不少，但底层仍然是SOCP，求解能力可以继续复用。

6. 实操中常见的坑与排查经验

6.1 高频报错与解决方案速查

6.2 几个降低踩坑概率的习惯

第一，永远从“纯负荷、无DG、无调节”的基准模型开始验证。33节点系统如果连基础的潮流约束都过不去，说明代码里一定有写错的地方。这时候不要急着调求解器参数。

第二，设置电压上下限时，不要把上限设到1.2pu这种“绝对安全”的宽松值。宽松的电压约束表面上降低了不可行风险，实际上会让求解器选择一些极端电压工况配合网损最小化，结果非常失真。更符合实际的区间是0.94~1.06pu。

第三，建模之前先手算一个小规模的3节点或5节点示例，把约束写成矩阵形式，验证过一遍再上33节点。很多人直接跳到33节点，一报错就不知道是模型问题还是数据问题。从简到繁，是我屡试不爽的排障路径。

第四，求解完成后别只盯着网损数字，一定要看一眼电压分布曲线。电压曲线的形状往往比总网损更能说明优化质量。比如出现某个节点电压明显凹陷，大概率是那里的无功源没有参与调节。

最后说一个我自己的习惯：每当锥松弛残差超标时，我会先检查目标函数的每一项权重——很多时候不是模型错了，而是权重把求解器“带偏”了。只要目标里始终保留一个足够大的网损分量，松弛紧性很少会出问题。这套方法我在多个配电网算例上验证过，无论是最基础的33节点，还是后续扩展了分布式电源、储能接入的版本，稳得住。