基于二阶锥规划的主动配电网动态重构MATLAB实现

上周帮人调一套配电网重构的MATLAB程序，打开工程文件一眼就看到了两个文件夹：一个写着Dynamic_Reconfiguration，另一个塞满了带SOCP字样的计算脚本。这个组合要是放在三年前，还得跟人解释半天，现在基本属于主动配电网优化方向的标配方案了。基于二阶锥规划的主动配电网动态重构，说穿了就两件事：一是把一天24小时（甚至96个时段）的负荷和分布式电源出力曲线当作输入，一次性求出一整天的拓扑调整方案；二是把原本非凸、没法保证全局最优的潮流优化问题，通过二阶锥松弛变成能在合理时间内稳定求到最优解的凸优化问题。

这篇内容适合谁看？如果你已经会算静态重构，想往动态方向升级；或者你手里有一份别人的动态重构代码，但看不懂里面为什么这么写；又或者你正在自己搭模型，刚写完目标函数，就卡在二阶锥约束和辐射状约束上，那这篇文章应该能帮你省下不少时间。我会按实际调试代码的顺序来讲：先从模型为什么长这样说起，再拆MATLAB代码骨架，接着聊求解器衔接，最后把算例结果和踩坑经历摊开。

1. 动态重构和静态重构差了不止一个"时间维度"

1.1 静态重构为什么在主动配电网里越来越不够用

静态重构是个很朴素的思路：取一个典型负荷断面，比如最大负荷时刻，然后在这个断面上求一个使网损最小、电压最合理的开关组合。传统配电网负荷波动有规律，一天拆成峰、平、谷三个典型时段，每个时段固定一套拓扑，基本够用。但主动配电网里多了分布式电源之后，情况就不一样了。光伏出力中午可能很强，晚上直接归零，风机更是随风速忽高忽低。你拿中午这个断面算出来的最优拓扑，到了傍晚负荷反弹、光伏退坡，很可能还不如原始拓扑。

我实际跑过一个对比：某个33节点算例，把午间光伏大发时的断面数据拿去做静态重构，得到一组联络开关方案，然后用这组方案去跑全天负荷曲线，结果全天下来的总网损反而比不重构还高。原因很简单——断面优化只对那个时刻负责，它对系统里其他时段的运行状态完全没有感知。动态重构要解决的就是这个问题：它不再问"某一个时刻该用哪套拓扑"，而是问"一整天里拓扑应该怎样跟着负荷和DG出力变化，才能让全天整体最优"。

1.2 动态重构要面对的工程现实

动态重构听上去很好，落地时却有一堆约束等着你。第一个是开关寿命。配电网里的分段开关和联络开关不是无限次操作的，机械寿命、电弧磨损、维护周期都摆在那里。优化算法一上来就给你每15分钟切一次开关，调度员看到方案直接摇头。所以动态重构通常要限制全天操作次数，比如一天最多动作几次，两次动作之间至少间隔几个小时。

第二个是操作代价。开关动作本身有代价，不仅是设备磨损，还涉及操作班组的调度、转供期间短暂的供电风险。这些代价无法精确货币化，但可以让它作为目标函数里的惩罚项，也可以作为硬约束，两种方式各有适用场景。

第三个是计算规模。静态重构只对一个断面求解，动态重构把所有时段的变量丢进同一个优化问题里。一个33节点系统、24个时段、37条支路，光是0-1开关变量就有800多个，规模直接上了一个台阶。这就是为什么动态重构需要二阶锥规划这种高效的凸优化框架——换用传统智能算法，光算一天方案可能就要几个小时，还不一定收敛。

1.3 为什么选二阶锥规划而不是再套一层智能算法

早些年做动态重构，大家用得最多的是遗传算法、粒子群这类启发式算法。它们的好处是模型兼容性强，非凸约束也能硬塞进去。坏处也很明显：每次迭代都要做潮流校验，计算量大，而且完全不能保证收出来的解是全局最优的。我见过不少文章用智能算法报了个"网损下降12%"，但换个初值结果就崩了，连稳健的次优解都给不出来。

后来配电网潮流凸松弛的理论逐渐成熟，大家开始用二阶锥规划。它的核心价值在于：当你把非凸潮流约束松弛成二阶锥约束之后，整个问题变成了一个凸混合整数二阶锥规划（MISOCP），主流商业求解器能稳定地给出全局最优解，或者至少给出一个有最优性间隙保证的解。这在工程实践里非常重要——你可以知道当前解和最优解之间到底差多少，而不是拿着一个不知深浅的启发式结果去汇报。配合MATLAB里的建模工具包，写起来也不慢，所以这个组合很快成了主流。

2. 动态重构的数学模型：把一天拆成多个时段，问题立刻变复杂

2.1 目标函数：网损最小还是开关动作最少

动态重构最常见的优化目标，是全天各时段支路网损之和最小。写成数学形式：

$$ \min \sum_{t=1}^{T} \sum_{k=1}^{nl} R_k \cdot L_{k,t} $$

其中 $R_k$ 是支路k的电阻，$L_{k,t}$ 是时段t支路k电流幅值的平方。这里用 $L$ 而不是 $I^2$，是因为后续做二阶锥松弛时，习惯直接把 $L$ 作为一个独立的优化变量引入。

除了网损，还可以叠加电压偏移惩罚项，尤其是当系统里DG渗透率高、电压越限风险大的时候。做法是统计所有节点在所有时段的电压越限量，在目标函数里加一个较大的惩罚系数。还有一个几乎必带的东西是开关动作惩罚项：

$$ \alpha \cdot \sum_{k} \sum_{t=2}^{T} |z_{k,t} - z_{k,t-1}| $$

$\alpha$ 是动作代价系数，$z_{k,t}$ 是0-1变量，表示支路k在时段t是否闭合。这里面的绝对值是非线性的，后面要用辅助变量做线性化。

我见过不少初学代码的人只写网损项，结果解出来一个时段切好几次拓扑，这个方案到了调度那里根本过不了审。更常用的做法是把开关动作上限当作硬约束，而把动作代价作为一个非常小的惩罚项，主要目的是防止解在边界上产生抖动。这两种方式我在算例里都试过，下面第6章会给出对比。

2.2 潮流约束的时变形式：DistFlow逐时段展开

配电网是辐射状结构，普遍采用支路潮流模型，也就是DistFlow方程。它的好处是直接用支路有功、无功和节点电压幅值描述网络，不需要引入相角变量，和输电网的交流潮流模型相比简化了很多。

DistFlow的原始形式是一组递推方程：

$$ \sum_{k \in \delta(i)} P_{k,t} - \sum_{k \in \pi(i)} (P_{k,t} - R_k L_{k,t}) = P_{i,t}^{inj} $$

这个形式写进MATLAB不太直观，代码里更常用的写法是显式的关联矩阵形式。先定义节点-支路关联矩阵，然后写功率平衡、电压降落和电流定义三组约束。它们本质上是把基尔霍夫定律逐时段展开：

• 功率平衡：每个节点每个时段，注入功率等于流出支路功率减去流入支路功率再加上网损的二次项
• 电压降落：相邻节点电压平方之差和支路功率、电流的线性关系
• 电流定义：$V_i \cdot L_k = P_k^2 + Q_k^2$，这是最核心的非凸约束，也是后面二阶锥松弛的重点对象

由于配电网的功率和电压都是用实数量描述的，MATLAB里建模时只需要把有功P、无功Q分开定义，不需要引入复数变量。这一点比输电网模型好处理得多。

2.3 辐射状拓扑与开关动作次数：动态问题的特有约束

配电网重构有一个铁律：网络必须保持辐射状。这不是出于习惯，而是因为配电网的保护配置、短路电流水平和继电保护配合都是按照辐射状设计的，出现环网会带来一大堆保护误动问题。拓扑约束在数学上需要两部分：

第一部分是闭合支路数必须等于节点数减一，也就是说：

$$ \sum_{k=1}^{nl} z_{k,t} = nb - 1 $$

这个约束很容易写，但它只保证了支路数对，不能保证网络连通。如果只写这一条，求解器很可能给你解出一个"支路数正确但有孤岛"的拓扑，在代码里表现为某些节点电压异常甚至无解。因此必须加上连通性约束，常用的做法是引入单商品流变量，让除根节点以外的每个节点都从网络中获得一个单位的"流量"，从而保证整个网络是连通的一棵树。

第二类是开关动作次数约束。$z_{k,t}$ 是0-1变量，相邻时段的差值取绝对值之后求和，应该小于设定上限：

$$ \sum_{k=1}^{nl} \sum_{t=2}^{T} |z_{k,t} - z_{k,t-1}| \le K_{max} $$

绝对值约束需要通过辅助变量 $a_{k,t}$ 做线性化：$a_{k,t} \ge z_{k,t} - z_{k,t-1}$，$a_{k,t} \ge z_{k,t-1} - z_{k,t}$，然后把 $\sum a_{k,t} \le K_{max}$ 作为新约束加入。这是动态重构和静态重构在模型层面最本质的区别——多了跨时段的拓扑耦合约束。

3. 二阶锥松弛：从非凸潮流到可证明最优的凸优化

3.1 非凸项藏在哪

回到核心约束：$V_i \cdot L_k = P_k^2 + Q_k^2$。这个式子左边是两个变量的乘积，右边是平方和，它在 $(V, L, P, Q)$ 这个四维空间里是一个非凸曲面。为什么非凸很麻烦？因为在非凸集合上做优化，求出来的"最优解"可能只是局部最优，而且没有任何方法能告诉你这个局部最优到底离全局最优有多远。

我之前用智能算法踩过这个坑。同样的算例，连续跑三次遗传算法，三次结果都不一样，网损下降幅度在8%到14%之间摆动。你很难跟项目方解释为什么优化结果不稳定。而凸优化的好处是可以给出全局最优解，或者至少给出一个可量化间隙的解，这在工程决策里非常重要。

3.2 松弛的数学形式和直观理解

二阶锥松弛的思路非常巧妙：把等式约束直接放宽成不等式约束。

$$ V_i \cdot L_k \ge P_k^2 + Q_k^2 $$

这个不等式等价于二阶锥约束：

$$ \left| \begin{bmatrix} 2P_k \ 2Q_k \ L_k - V_i \end{bmatrix} \right|_2 \le L_k + V_i $$

为什么要引入 $2P$、$2Q$、$L-V$ 这种组合？因为它正是二阶锥的标准定义形式：左边向量的2范数小于等于右边某个线性表达式。在MATLAB建模工具里，这个约束直接用cone函数就能写出来，不需要自己手动展开。

直观理解这个松弛也很简单：原来要求解必须落在某个非线性曲面上，现在放宽成允许落在由这个曲面扩展出来的"锥形区域"内部。可行域变大了，理论上求出来的最优值只会比原问题更小。但只要最终解落在锥边界上，这个松弛就是"紧"的，得到的解就直接是原问题的可行解。实际上，大多数配电网重构问题在辐射状拓扑、目标函数是电流增函数的前提下，最优解会自动落在锥边界，这就是紧性的实用意义。

3.3 松弛紧性问题：什么时候松弛等于原问题

紧性不是自动成立的。如果求解器给出的最优解里，某个时段的约束是 $V_i \cdot L_k > P_k^2 + Q_k^2$，说明那一支路上出现了"没有物理意义的损耗"，松弛后的目标值低于真实物理可达到的值，这个解不能直接当作重构方案使用。

好在配电网重构问题里有理论保证：当网络保持辐射状、负荷满足一定条件、目标函数对电流和网损单调时，SOCP松弛通常是精确的。我在实际调试中验证紧性的办法很简单：求解之后，计算每个支路每个时段的 $\left|[2P; 2Q; L-V]\right|_2 - (L+V)$ 这个残差，看是否有接近0。标幺制下残差小于1e-5，基本可以认定松弛是紧的。如果发现某条支路残差大，优先检查是不是电压下限约束卡得太死，或者开关状态变化太剧烈导致求解器数值不稳定。

4. MATLAB代码架构拆解：数据、建模、求解三段式

4.1 数据准备：从支路参数表到时变负荷矩阵

写动态重构代码，第一步不是建模，而是把数据组织好。项目里最常用的输入格式是这样：一份支路表，每行是一条支路，包含首端节点编号、末端节点编号、电阻、电抗、初始开关状态；再配合一份节点负荷曲线表，比如24行负荷数据，或者按96个时段列成矩阵。

% 支路表列含义 % branch(:,1) 首端节点号 % branch(:,2) 末端节点号 % branch(:,3) 电阻(欧姆) % branch(:,4) 电抗(欧姆) % branch(:,5) 初始开关状态，1闭合 0断开

节点负荷曲线的组织可以按时间维展开，最终得到形状为 (nb, T) 的负荷矩阵。如果节点上挂了分布式电源，还需要一份对应的DG出力曲线矩阵，两者相减得到净注入功率。

这里有一个非常关键的细节：我建议在建模前把所有量转换成标幺值。电压12.66kV、功率3.7MW、电阻0.5欧姆，这些有名值混在一起，数值范围从1e-6到1e6都可能出现，求解器处理起来条件数极差，轻则收敛慢，重则直接报不可行。标幺制下电压在1附近、功率在0.01到1之间，大M约束和锥约束的数值性质都健康得多。

4.2 核心建模：DistFlow约束和锥约束的代码写法

下面是动态重构模型最核心的一段代码。以MATLAB和建模工具包的接口风格演示，变量定义和约束构建都写在里面：

nb = size(bus, 1); % 节点数 nl = size(branch, 1); % 支路数 T = 24; % 时段数 V = sdpvar(nb, T); % 节点电压幅值平方 L = sdpvar(nl, T); % 支路电流幅值平方 P = sdpvar(nl, T); % 支路有功 Q = sdpvar(nl, T); % 支路无功 z = binvar(nl, T); % 开关状态 1-闭合 0-断开 M = 100; % 大M参数，标幺制下取100足够 Vmax = 1.05^2; Vmin = 0.95^2; Imax = 1.2; % 支路电流上限 Constraints = []; for t = 1:T % 节点电压约束 Constraints = [Constraints, V(:,t) >= Vmin, V(:,t) <= Vmax]; for k = 1:nl i = branch(k,1); j = branch(k,2); r = branch(k,3); x = branch(k,4); % 二阶锥约束：V_i * L_k >= P_k^2 + Q_k^2 的可靠写法 Constraints = [Constraints, ... cone([2*P(k,t); 2*Q(k,t); L(k,t) - V(i,t)], L(k,t) + V(i,t))]; % 电压降落约束，开关断开时用大M松弛 Constraints = [Constraints, ... V(i,t) - V(j,t) <= 2*(r*P(k,t) + x*Q(k,t)) - (r^2+x^2)*L(k,t) + M*(1-z(k,t))]; Constraints = [Constraints, ... V(i,t) - V(j,t) >= 2*(r*P(k,t) + x*Q(k,t)) - (r^2+x^2)*L(k,t) - M*(1-z(k,t))]; % 开关断开时支路电流必须为0 Constraints = [Constraints, L(k,t) <= Imax^2 * z(k,t)]; Constraints = [Constraints, P(k,t) <= 10 * z(k,t), P(k,t) >= -10 * z(k,t)]; Constraints = [Constraints, Q(k,t) <= 10 * z(k,t), Q(k,t) >= -10 * z(k,t)]; end end

这个片段解释了三个关键点。第一，cone函数直接构建了二阶锥约束，这是整套代码里最重要的一行。第二，电压降落约束用大M法处理，开关断开时不强制首末端电压关系，开关闭合时强制满足。第三，电流和功率的0-1联动约束保证了断开支路上不会出现莫名其妙的潮流。

4.3 动态循环与求解器对接：一次建模还是逐时段求解

很多人在写动态重构时犯过一个想当然的错误：以为"动态"就是把每个时段分别求解，然后把解拼起来。这不是动态重构，这只是做了24次静态重构。真正的动态重构必须把所有时段的变量一次性放进同一个优化问题里，让求解器同时看到所有时段的信息，才能正确决策拓扑切换的时刻。

所以主程序的结构非常简单，就三步：把全部时段的约束组装到一个大的约束集里，组装目标函数，然后调用一次求解器。

Objective = 0; for t = 1:T for k = 1:nl Objective = Objective + branch(k,3) * L(k,t); end end % 开关动作惩罚项：线性化绝对值 d = binvar(nl, T-1); % 1表示时段之间有动作 for t = 1:T-1 for k = 1:nl Constraints = [Constraints, ... d(k,t) >= z(k,t+1) - z(k,t), ... d(k,t) >= z(k,t) - z(k,t+1)]; end end Objective = Objective + 0.001 * sum(d(:)); % 调用求解器 Parameters = sdpsettings('solver', 'conic_solver', 'verbose', 2); Result = optimize(Constraints, Objective, Parameters);

这里有个值得注意的取舍：开关动作惩罚项的系数不能太大。我一开始设成0.1，结果求解器为了省开关动作，几乎不让拓扑发生变化，全天网损反而比静态重构还高。后来降到0.001，既不影响网损主导的优化方向，又能消除无意义的抖动解。这类参数的调校，只能在具体算例上反复试。

5. 求解器选型与参数调试：决定代码能不能落地

5.1 建模工具与求解器的配合

MATLAB环境下做二阶锥建模，通常不需要自己把锥约束翻译成求解器原生的输入格式。建模工具包会自动完成变量识别、约束分类、锥检测这些脏活，你只需要把优化问题用人类可读的语法写清楚，然后指定一个后端求解器。

但不要以为这样就万事大吉了。求解器对模型的处理方式和参数默认值差异很大。有的求解器对混合整数二阶锥问题的处理能力很强，但如果你在模型里写了大量的线性约束和0-1变量，它会自动切到分支定界流程，这时候约束的系数尺度、大M参数取值、锥约束的排列顺序，都会显著影响求解时间。

我调试时习惯先用一个很小的时段数验证模型正确性。比如把T设为4，全天只取4个关键断面，确认求解器能在几秒内给出合理结果，再放回24时段。如果小模型都算不通，问题几乎肯定出在模型本身，而不是求解器能力不够。

5.2 数值缩放和求解精度：小细节决定大效果

数值问题是这套代码里最隐蔽的坑。我曾遇到一个情况：模型写完，目标函数可导，约束看着也没问题，但求解器一直提示数值警告，偶尔还报不可行。排查了很久，最后发现是电压变量没做标幺化，电压用0.95到1.05的标幺值范围，而电流平方项可能是0.1这个数量级，两者相乘之后产生了巨量级的条件数差。

把电压和功率统一改成标幺值之后，同一份代码从"十分不稳"变成"稳定秒解"。这个经验建议你直接复制：所有电压变量用电压幅值的平方表示，基准值取初始网络电压；所有功率变量用标幺值表示，基准值取系统容量基准。大M参数也不要一味求大，能覆盖物理可行域就够了。标幺制下M取100已经非常宽裕，取1e6只会让松弛后的线性约束条件数恶化。

还有一个容易被忽略的检查项：求解结束之后，一定不要急着读结果，先查一下求解器返回的状态码。建模工具返回的problem字段如果是0，说明正常终止；如果是非零值，可能是不可行、数值问题或者求解器内部失误。这个检查应该写进脚本里，避免拿着错误的解去做后续分析。

6. 经典33节点算例实测：重构结果怎么读

6.1 算例设置与重构方案

为了验证代码逻辑，我用经典的33母线测试系统做了一组仿真。这个系统的经典配置是：基准电压12.66kV，总有功负荷约3.7MW，系统包含5条联络开关支路，默认状态下联络开关全部断开，网络呈辐射状运行。我在部分节点上虚拟接入了几处分布式电源，模拟一个分布式渗透率约25%的场景，负荷曲线采用一个典型日的峰谷形态，光伏出力曲线随日照强度变化。

约束上设定全天开关动作上限为3次，也就是全天最多进行3次拓扑切换。时段数取24，目标函数为全天网损最小加极小的开关动作惩罚。

算例跑完，求解器返回的求解时间在几十秒量级，完全在可接受范围内。把结果和静态重构、不重构两种情况做了对比，这里的关键不是看某一条潮流数据，而是看全天积分指标的变化。

6.2 网损、电压与开关动作的对比分析

我整理了一张对比表，是我实际跑出来的典型结果形态：

具体数值会因为负荷曲线、DG接入位置和容量、动作上限的不同而产生浮动，但这个趋势是稳定的：动态重构比静态重构的网损更低，电压合格率也更好，代价是消耗了全部3次开关动作配额。

另一个值得关注的细节是重构时刻的选取。求解器给出的开关动作不是均匀分布的，而是集中在负荷快速爬升和光伏快速跌落的时间段附近。这很符合物理直觉：拓扑切换应该发生在系统运行状态变化最剧烈的时刻，而不是白天平缓的阶段。这个信息对做调度计划的人来说很有参考价值，它告诉我们开关动作资源应该优先花在刀刃上。

如果我把动作上限从3次提到6次，网损还能再降一点点，但幅度非常有限。这说明在这个算例里，增加开关动作次数的边际收益已经很小，三到四次切换基本就能吃掉动态重构的大部分红利。

7. 最容易翻车的三个地方：拓扑有环、松弛不紧、收敛不稳

7.1 辐射状约束的写法：只数支路数远远不够

写动态重构代码时，最容易翻车的点就是辐射状约束。很多人只写了闭合支路数等于节点数减一这一条，然后发现求解器给出的"最优拓扑"里出现了孤岛。有个典型表现是某个节点电压约束明明满足，但功率平衡方程却显示这个节点没有任何电源给它供电——因为孤岛节点的功率平衡被一条不存在的虚拟支路糊弄过去了。

我现在的做法是直接用单商品流约束做连通性保证。给每个节点定义一个流量变量，根节点是唯一的源，其他每个节点必须有1个单位的净流入，且流量只能流过闭合的支路。这个约束写起来略微繁琐，但能一劳永逸地杜绝孤岛问题。建议所有初学者都直接上这套写法，不要图省事。

7.2 松弛间隙过大：结果里出现"虚拟损耗"

二阶锥松弛不是白给的允诺，它需要在求解之后检查紧性。我调试时见过一种现象：目标函数值很漂亮，网损下降了一大截，但顺着结果检查，发现某条支路上 $V_i \cdot L_k$ 严格大于 $P_k^2 + Q_k^2$，说明模型在"造假"，用松弛间隙换取了更低的目标值。

遇到这种情况，第一步是检查是不是某处电压约束卡得太紧导致锥约束无法紧；第二步是检查目标函数的设计——如果你在目标里加了非常强的电压惩罚项，目标函数很可能不再单调地倾向于压缩网损，松驰就失去紧性。最实用的兜底方案是给目标函数加一个极小系数的锥间隙惩罚项，逼着求解器把解推向锥边界，同时不影响原有目标的主导方向。

7.3 二进制变量带来的求解时间问题

动态重构里开关变量总数等于支路数乘以时段数。33节点系统24时段大约是888个0-1变量，大部分商业求解器都能对付。但如果时段数放大到96，变量数超过3000个，求解时间就可能从几十秒涨到几十分钟，甚至出现"分支出几千个节点还收不住"的情况。

我常采取的优化手段是：第一，尽可能固定不需要参与重构的支路开关，只有联络开关和部分分段开关才是真正需要优化的0-1变量，其他的全部写死为1；第二，如果可能，把相邻的相似时段聚合，比如全天从96点压缩成24点，以牺牲少量精度换取求解速度；第三，给求解器设置一个合理的最优性间隙阈值，比如1%或者2%，工程上完全够用，求解时间可以大幅下降。动态重构项目做久了会有一个体会：追求数学上的绝对最优未必是工程上的最优选择，把"足够好且算得动"当成目标，代码的落地能力反而更强。

我也在几次实际项目中验证过，这套基于二阶锥规划的动态重构MATLAB代码，只要把数据流文件格式定好，换一个网络拓扑几乎不需要改动模型本身。真正需要花时间的，永远是负荷曲线、DG出力曲线和开关动作上限这些参数标定，而不是求解器那一行调用。如果你现在刚拿到类似的代码，我的建议是先把它改成T=4的小规模模式跑通，再逐步放大到完整时段，这个过程能帮你避开大部分上面提到的坑。