配电网故障恢复这个话题，这几年在电力系统里热度一直很高。我自己把“同时考虑孤岛与重构的故障恢复策略”这个方向完整复现了一遍，代码是基于Matlab实现的。坦率说，刚拿到这个题目时我觉得并不复杂，无非就是“断开故障支路、闭合联络开关、让失电负荷恢复供电”那一套经典做法。但真正开始建模才发现，把孤岛运行和网络重构放到同一个优化框架里求解，坑远比想象中多——尤其是拓扑连通性约束、潮流收敛性、二阶锥松弛紧度这几关，处理不好，结果根本无法落地。

这篇文章不是教科书式的理论复述，而是我从零开始复现这套策略的完整记录，包括数学建模、Matlab代码实现、求解器调试，以及我自己踩过的好几个坑。如果你想做的是配电网故障恢复、分布式电源孤岛划分、或者网络重构方向的课题，这篇文章可以直接拿来做技术底稿。

1. 项目定位与核心思路拆解

1.1 为什么要把孤岛和重构放在同一个模型里

传统配电网故障恢复，基本思路是“重构优先”：故障发生后，通过调整分段开关和联络开关的状态，把非故障失电区域的负荷转移到其他馈线上，由主网继续供电。重构解决的是“怎么切、怎么合”的问题。

但随着分布式电源接入比例越来越高，情况发生了变化。故障隔离之后，如果故障区域下游还有分布式电源（比如光伏、储能、小型燃气轮机），那么靠DG给孤岛内的负荷供电是完全可行的——这就是“孤岛运行”。孤岛的好处很明显：少切负荷、减少停电范围、提高供电可靠性。

问题是，以往很多研究把这两件事分开做：先划分孤岛，再在网络中搜索重构方案；或者反过来，先重构再考虑DG是否能形成有效孤岛。分开做的结果经常会打架——你按重构方案合上了一个开关，结果发现某个DG节点被隔离开了，孤岛根本形成不了；或者你先固定了孤岛范围，结果重构空间被锁死，白白丢了很多能恢复的负荷。

所以，这个项目标题里“同时考虑”三个字才是核心痛点。它把孤岛划分和网络重构统一成一个组合优化问题，用一套决策变量同时决定：哪些开关分闸、哪些开关合闸、哪些DG独立成岛、岛内带哪些负荷。这样做的好处是可以在全局最优层面协调“主网供电 + 孤岛供电”两种恢复手段，而不是贪心地分步求解。

1.2 复现这个项目的技术路径

我复现时把整体流程拆成了四段：

1. 故障场景注入：指定某条支路发生永久性故障，断开故障支路两侧开关，得到失电区域；
2. 联合优化建模：以网损最小、失电负荷最小、开关操作次数最小为多目标，建立带辐射状拓扑约束和潮流约束的混合整数优化模型；
3. 二阶锥松弛处理：把配电网潮流方程中的非线性项通过变量替换和锥松弛转化为二阶锥约束，交给商业求解器求解；
4. 结果校验与可视化：求解得到开关状态和孤岛范围后，重新做确定性潮流校验，确认电压、电流在允许范围内。

整个链路中最容易出问题的环节是第三步到第四步之间——优化模型给出的解在理论上“可行”，但物理验算时可能会发现潮流不收敛，或者某些孤岛内根本没有电源节点。这类问题在后面会专门展开讲。

2. 数学模型构建：从物理问题到优化公式

2.1 目标函数：多目标如何加权

配电网故障恢复的目标不是单一的。我在建模时做了三层目标：

• 失电负荷量最小：这是最核心的指标，恢复供电的负荷越多越好；
• 开关操作次数最少：每个开关动作一次都有机械寿命损耗，而且频繁操作会增加误动风险；
• 网损最小：在恢复供电的前提下，让系统运行尽量经济。

数学上，目标函数写成：

min F = α₁ · Σ(λᵢ · P_ᵈ,ᵢ · (1 - xᵢ)) + α₂ · Σ(zₖ) + α₃ · Σ(rᵢⱼ · Iᵢⱼ²)

其中，xᵢ是负荷节点的供电状态变量（1表示已恢复供电，0表示失电），λᵢ表示负荷的重要程度权重（一二级负荷权重更高），zₖ表示开关状态是否发生变化（0/1变量），Iᵢⱼ表示支路电流，rᵢⱼ表示支路电阻。

三个目标的权重系数如何设置，直接影响最终恢复策略的形态。我的经验是，α₁必须远大于α₂和α₃，比如取100:1:0.01。理由是：故障恢复场景中“少停电”是压倒性目标，如果权重均衡，求解器会为了省两次开关操作而牺牲几百千瓦的负荷恢复，这种方案工程上完全不可接受。但α₂又必须存在，否则求解器会倾向把馈线上所有开关都拉闸再重新合闸——这样确实能构造出最优拓扑，但实际操作风险和代价太大。

注意：如果代码里用的是线性加权，一定要把α₁、α₂、α₃的量纲考虑进去。失电负荷通常是kW量级（几百到几千），开关操作次数是个位数，网损也是kW量级。直接按数值大小给权重会失真，建议先做归一化。

2.2 约束条件：不只是潮流方程

目标函数决定“往哪个方向优化”，真正让问题变得难解的是约束条件。我梳理下来，这个模型至少要包含五类约束：

（1）潮流约束（DistFlow方程）

配电网是辐射状网络，采用DistFlow模型非常合适。对每条支路(i,j)：

Pⱼ = Pᵢ - rᵢⱼ·(Pᵢ² + Qᵢ²)/Vᵢ² - P_负荷,ⱼ + P_发电,ⱼ

Qⱼ = Qᵢ - xᵢⱼ·(Pᵢ² + Qᵢ²)/Vᵢ² - Q_负荷,ⱼ + Q_发电,ⱼ

Vⱼ² = Vᵢ² - 2·(rᵢⱼ·Pᵢ + xᵢⱼ·Qᵢ) + (rᵢⱼ² + xᵢⱼ²)·(Pᵢ² + Qᵢ²)/Vᵢ²

这里P、Q是支路潮流，V是节点电压。注意这个方程本身有非线性项(P²+Q²)/V²，后面做二阶锥松弛就是从这里入手。

（2）节点电压与支路电流约束

V_min ≤ Vᵢ ≤ V_max，通常取0.95到1.05标幺值；Iᵢⱼ ≤ I_max。这个约束在故障恢复中特别容易触发——因为重构后供电半径变长，线路末端电压很容易越限。

（3）辐射状拓扑约束

配电网正常运行时必须是辐射状，不允许闭环。数学上要求：所有被投入的支路组成的图是连通的，且支路数 = 节点数 - 孤岛数。这里孤岛数是一个随拓扑变化的整数，这是整套模型里最棘手的约束之一，后面会说具体处理办法。

（4）DG出力约束

0 ≤ P_DG ≤ P_DG_max，以及最小运行时间约束（这里简化为出力上下限即可）。如果在孤岛内，DG必须能支撑岛内负荷的功率平衡，否则孤岛无意义。

（5）开关状态与网络拓扑一致性约束

开关变量是0/1的，支路如果断开，则该支路潮流必须为0。这个约束用Big-M法处理：

-P_BigM · (1 - zₖ) ≤ Pₖ ≤ P_BigM · (1 - zₖ)

这个Big-M参数选多大很有讲究。M太大，会让求解器数值稳定性变差；M太小，又会砍掉可行域。我在33节点系统里取M = 10 MW，比系统最大负荷功率高一个量级就够用。

2.3 二阶锥松弛：非线性问题的标准解法

DistFlow方程里的非线性项(P² + Q²)/V²，处理方式有几种——直接非线性求解、线性化近似、二阶锥松弛。目前工程复现最推荐的是二阶锥松弛（Second-Order Cone Relaxation）。

做法是引入辅助变量： lᵢⱼ = (Pᵢⱼ² + Qᵢⱼ²) / Vᵢ² uᵢ = Vᵢ²

经过变量替换，潮流方程变为线性方程： ΣP = Pᵢ - rᵢⱼ·lᵢⱼ - P_负荷,ⱼ + P_发电,ⱼ uⱼ = uᵢ - 2·(rᵢⱼ·Pᵢ + xᵢⱼ·Qᵢ) + (rᵢⱼ² + xᵢⱼ²)·lᵢⱼ

而描述P、Q、l之间关系的约束则松弛为旋转锥约束： lᵢⱼ ≥ (Pᵢⱼ² + Qᵢⱼ²) / uᵢ

等价地写为二阶锥格式： ‖[2·Pᵢⱼ, 2·Qᵢⱼ, lᵢⱼ - uᵢ]‖₂ ≤ lᵢⱼ + uᵢ

为什么要“松弛”？因为原来等式约束是精确物理关系，松弛成不等式后可行域变大了，问题从非凸变成凸问题，再加上0/1变量后是混合整数二阶锥规划（MISOCP），主流求解器可以直接处理。

但松弛有一个隐患：如果最优解里某些支路的锥约束不是“紧”的，即l值明显大于物理实际值，那么这个解虽然在数学可行，物理上却不可行。我在第4节会讲怎么检查和处理。

3. 求解方法与Matlab代码实现

3.1 求解器选型：CPLEX还是Gurobi

这一节先给结论：Matlab环境下的复现，强烈推荐用YALMIP工具箱 + 商业求解器（CPLEX或Gurobi）。

为什么不推荐自己写粒子群、遗传算法？我复现初期也试过用智能算法做联合优化。粒子群的优点是实现简单、不需要凸优化知识；缺点是极容易陷入局部最优，而且每次运行结果还不一样——这对故障恢复策略来说是不可接受的。工程上，故障恢复方案需要在分钟内给出确定性的、可验证的最优解，MISOCP是更稳妥的路线。

YALMIP的安装已经比较成熟，主要步骤是：下载YALMIP源码包，添加到Matlab路径，再安装CPLEX或Gurobi求解器，并在YALMIP中通过yalmip('clear')和solvesdp（或optimize）调用。

具体到本项目，求解器配置核心代码：

% 求解器配置 ops = sdpsettings('solver', 'cplex', ... % 指定求解器 'verbose', 2, ... % 打印求解过程 'showprogress', 1, ... % 显示进度 'cplex.mip.tolerances.mipgap', 0.0001, ... % MIP间隙 'cplex.mip.tolerances.integrality', 1e-5); % 整数容差

这里mipgap设置为1e-4表示求解器在最优解间隙小于0.01%时停止，这个精度对故障恢复问题已经足够。设得太小（比如1e-9）反而会让求解时间急剧上升，尤其在规模更大的配电网场景里。

3.2 数据准备：以33节点测试系统为例

复现时我选了一个经典的33节点配电测试系统，母线电压10kV，总负荷约3.7MW，含5条联络开关支路，我在其中4个节点接入了分布式电源（容量从100kW到300kW不等）。这个系统规模适中，既能体现联合优化的复杂度，又方便调试和可视化。

数据准备阶段最需要注意的是单位统一。我在代码里把功率单位统一成kW，电压用标幺值，阻抗用标幺值：

% 系统数据加载（部分节点数据示意） % 支路矩阵 columns: [起点, 终点, 电阻(pu), 电抗(pu), 容量(A), 初始状态] branch = [ 1, 2, 0.0055, 0.0030, 300, 1; 2, 3, 0.0030, 0.0015, 300, 1; % ... 其余支路 18, 33, 0.0040, 0.0020, 200, 0; % 联络开关, 初始断开 ]; % 节点数据 columns: [节点编号, 有功负荷(kW), 无功负荷(kvar)] bus = [ 1, 0, 0; 2, 100, 60; 3, 90, 40; % ... 其余节点 ]; % DG数据 columns: [节点编号, 有功上限(kW), 无功上限(kvar)] dg = [ 8, 200, 100; 12, 150, 80; 25, 300, 150; 30, 100, 50; ];

这里有个容易忽略的细节：联络支路的初始状态是0（断开），但故障恢复时它可能变成1（闭合）。在定义变量时，我建议直接对所有支路定义二进制变量，只是用初始状态做目标函数中的“开关动作次数”计算基准。

3.3 核心约束的YALMIP实现

定义决策变量：

% 支路状态变量，1表示闭合，0表示断开 z = binvar(n_branch, 1); % 节点供电状态，1表示恢复供电，0表示失电 x = binvar(n_bus, 1); % 支路有功、无功、电流平方项，节点电压平方项 P_branch = sdpvar(n_branch, 1); Q_branch = sdpvar(n_branch, 1); l_branch = sdpvar(n_branch, 1); u_node = sdpvar(n_bus, 1); % DG出力变量 P_dg = sdpvar(n_dg, 1); Q_dg = sdpvar(n_dg, 1);

潮流约束（这里是关键代码，每一行都有讲究）：

% 对每条支路构造DistFlow方程 for k = 1:n_branch i = branch(k, 1); j = branch(k, 2); r = branch(k, 3); x = branch(k, 4); % 节点j的有功平衡 % 流入 = 上一级注入 - 线路损耗 + DG出力 - 负荷（未恢复部分不计） constraints = [constraints, P_branch(k) - r * l_branch(k) + P_dg_at(j) * x(j) - P_load(j) * x(j) == ... sum(P_branch(find(branch(:,1) == j)))]; end

这一段初学者很容易写错。故障恢复场景下，失电负荷的“重要性”其实是通过x变量体现的——负荷没恢复（x=0），它就不出现在功率平衡方程里，自然就不消耗功率。如果直接把负荷功率全部放进方程，相当于默认所有负荷都恢复了，那优化问题的可行域会一下子变得非常奇怪。

电压约束部分：

% 电压上下限 for i = 1:n_bus constraints = [constraints, u_min^2 <= u_node(i) <= u_max^2]; end % 电压降落方程 for k = 1:n_branch i = branch(k, 1); j = branch(k, 2); constraints = [constraints, u_node(j) == u_node(i) - 2*(r*P_branch(k) + x*Q_branch(k)) + ... (r^2 + x^2) * l_branch(k)]; end

二阶锥约束（用YALMIP的cone命令）：

% 旋转锥约束: l >= (P^2 + Q^2)/u for k = 1:n_branch constraints = [constraints, cone([2*P_branch(k); 2*Q_branch(k); l_branch(k) - u_node(branch(k,1))], ... l_branch(k) + u_node(branch(k,1)))]; end

辐射状约束是重点。我在代码中采用了一种经典做法——利用“支路功率方向”构造虚拟流。给每个节点分配一个虚拟需求（通常取1），每条闭合支路的虚拟流从电源节点流向负荷节点，从而保证连通性和无环性：

% 虚拟流约束 % 设根节点（变电站）的虚拟需求为 n_bus - 1，其他节点为1 f = sdpvar(n_branch, 1); % 虚拟流 for k = 1:n_branch i = branch(k, 1); j = branch(k, 2); % 虚拟流不能超过支路状态决定的上限（Big-M） constraints = [constraints, -n_bus * z(k) <= f(k) <= n_bus * z(k)]; end % 节点虚拟流平衡 for j = 2:n_bus % 变电站节点除外 inflow = sum(f(branch(:,2) == j)); % 流入j的虚拟流 outflow = sum(f(branch(:,1) == j)); % 流出j的虚拟流 constraints = [constraints, inflow - outflow == 1]; end

这套虚拟流约束是配电网重构建模的经典技巧。它的原理是：如果网络中有环，环上的虚拟流无法同时满足所有节点的平衡方程；如果某个节点与电源不连通，它的虚拟需求无法被满足。这个约束的好处是不需要显式判断孤岛数量，很适合与MISOCP结合。

3.4 目标函数的代码实现

最终的目标函数代码：

% 失电负荷惩罚项（一二级负荷权重更高） load_penalty = 0; for i = 1:n_bus load_penalty = load_penalty + load_weight(i) * P_load(i) * (1 - x(i)); end % 开关动作惩罚项 switch_penalty = 0; for k = 1:n_branch switch_penalty = switch_penalty + abs(z(k) - branch(k, 6)); end % 网损项 loss_penalty = sum(r_branch .* l_branch); % 总目标 objective = 100 * load_penalty + 1 * switch_penalty + 0.01 * loss_penalty; % 求解 optimize(constraints, objective, ops);

abs(z(k) - branch(k,6))的写法很巧妙——branch里第6列存的是开关初始状态，两者相减取绝对值，正好统计出动作次数。这里用小技巧避免了引入额外0/1变量，能显著减小MIP规模，求解速度至少提升30%。

4. 复现过程中常见的坑与排查技巧

4.1 二阶锥松弛不紧：结果“数学可行，物理可疑”

这是我的老实话报告：在使用33节点系统测试时，最优解里会偶发某些远端的支路锥松弛不紧，导致网损被低估。你算出来的“最优”网损是20kW，但实际验潮后可能是30kW，甚至潮流根本不收敛。

检查方法很简单：求解完成后，对所有闭合支路检查锥约束的松弛间隙：

% 计算松弛间隙 gap = zeros(n_branch, 1); for k = 1:n_branch if value(z(k)) > 0.5 P_val = value(P_branch(k)); Q_val = value(Q_branch(k)); l_val = value(l_branch(k)); u_val = value(u_node(branch(k,1))); gap(k) = (l_val * u_val - (P_val^2 + Q_val^2)) / (l_val * u_val); end end % 如果gap显著大于1e-4，说明该支路松弛不紧

处理办法有两个。第一个是对目标函数加一个小的惩罚项λ·Σl_branch（λ取1e-4左右），把松弛“压紧”；第二个更稳妥——把MISOCP求得的开关状态固定下来，重新用传统牛顿-拉夫逊潮流验算一遍，确认所有节点电压和支路潮流都满足约束。

我最终在代码里两种方案都做了：优化阶段用惩罚项收紧松弛，后处理阶段做潮流复核。这样的“双保险”让我能放心地输出恢复方案。

4.2 孤岛内部没有电源：连通性约束失效

有一次跑了半天，结果发现一个孤岛里压根没有DG节点——所有负荷都靠主网供电，理论上不可能成立。原因是我在虚拟流约束里只用了变电站作为唯一虚拟流源，孤岛内的DG节点虽然物理上有电源，但虚拟流模型中并没有体现出来，导致孤岛的连通性虽然满足，但它是个“无源之岛”。

解决方法是：虚拟流约束的源节点不只有变电站，还要包括所有DG节点。改写节点虚拟流平衡方程，让DG节点也满足“注入等于负荷”：

% 节点虚拟流平衡（修正版） for j = 2:n_bus inflow = sum(f(branch(:,2) == j)); outflow = sum(f(branch(:,1) == j)); % DG节点作为虚拟流源，多一个注入 if is_dg_node(j) constraints = [constraints, inflow + is_dg_node(j) - outflow == 1]; else constraints = [constraints, inflow - outflow == 1]; end end

实际上更系统的做法是：在网络重构问题中规定每个连通分量必须包含至少一个电源节点（变电站或DG）。这可以显式添加一组线性约束，也可以用虚拟流法隐含保证。如果用的是智能算法，就没法方便地表达这种约束，必须在解码时做连通性校验——这也是我坚持用MISOCP的另一个原因。

4.3 求解时间爆炸：从四个小时到四分钟

33节点系统的MISOCP规模其实不大，但如果约束写成循环嵌套的写法不当，YALMIP展开变量时会非常慢。我第一版代码直接跑了四个小时都没出结果。

后来定位到两个瓶颈：

• Big-M参数过大。多个约束里M都取了10000，导致MIP的线性松弛非常松，分支定界要探索的节点数暴增。把M统一收紧到与系统规模匹配的量级后，求解时间大幅下降。
• 约束循环里重复调用find。YALMIP在构造约束时，每一条constraints = [constraints, ...]都会触发一次约束合并。循环几百次之后，这部分时间开销不比求解器本身少。优化方式是把全部约束先放在元胞数组里，最后一次性合并：

c_list = cell(n_cons, 1); idx = 0; % ... 构造过程中逐一填入 c_list{idx} = expr; constraints = [c_list{:}];

优化后整体求解时间压缩到4分钟左右（CPLEX，mipgap默认值）。对故障恢复场景来说，这个速度完全可用。

4.4 常见问题速查表

5. 复现后的核心体会与扩展建议

整套复现跑通之后，我个人最深的感受是：配电网故障恢复的关键难点不在“恢复”而在“约束”。开关组合的数量虽然庞大，但真正难啃的是辐射状约束和孤岛连通性约束的数学表达——这一步建模质量直接决定了求解器的性能和结果的工程可接受度。

如果看完这篇文章后你也想复现，我的建议是：不要一上来就追求大系统。先把33节点系统、固定故障点、单DG场景跑通，确认二阶锥松弛、虚拟流约束、权重系数逻辑都正确，再逐步扩展到多故障、多DG、动态重构。每一步都先做潮流验证，确认优化结果在物理上经得起推敲。

此外，这套框架本身有很好的扩展空间。比如把故障恢复与储能调度联合优化，引入可移动式应急电源（比如发电车）的接入位置决策，或者在目标函数中加入关键负荷（医院、数据中心）的差异化供电保障系数——这些都能在当前模型基础上通过增加变量和约束实现，底层求解框架不用动。

由于我复现时在连通性约束上吃过亏，最后再说一个经验：每调一次虚拟流约束、每加一个新的DG节点，都要重新跑一次“无故障状态下求解是否给出初始网络”的自检。如果连无故障场景都给出荒诞的拓扑，那一定是约束逻辑出了问题，这时候去调权重、改参数都没有意义。先保证模型的“零故障正确性”，再去追求故障场景的最优性，整个复现过程才会走得比较顺。