1. 这个项目到底在解决什么问题

1.1 一次典型故障后的两个恢复手段

先说场景。你坐在配网调度台前，某条10kV馈线出口断路器跳闸，故障点定位结果显示中段某个位置发生了永久性故障。这时候下游一大片配电变压器全部失电，用户电话一个接一个打进来。传统处理方式是：派人巡线、隔离故障、然后人工合联络开关，把失电区域“倒”到相邻馈线上。这套做法依赖现场经验，算得慢，而且往往只考虑到“能不能送电”，没仔细想“送过去之后会不会过载、会不会电压越限”。

这几年分布式电源（光伏、储能、柴油机、小水电）越来越多，故障恢复又多了一个手段：让失电区域里那些带DG的支路主动形成孤岛运行，DG自己带着一部分负荷继续供电，不需要从变电站送电过去。这种“主动孤岛”不再是过去防孤岛保护里避之不及的故障孤岛，而是调度员刻意划分出来的供电孤岛，前提是孤岛内部功率能自平衡、电压能稳住。

所以，一个完整的配电网故障恢复策略，应该同时考虑两条路子：重构——通过操作分段开关和联络开关改变网络拓扑，把负荷转移到健康馈线；孤岛——让DG失电后自发组织成独立的小供电网格。这个项目复现的正是这种“孤岛划分+网络重构”联合优化的运行策略，核心载体是Matlab代码。

1.2 分开做为什么不靠谱

很多早期做法是两步走：先用启发式算法找一个重构方案，把能转供的负荷都转走，剩下的实在没辙了再划孤岛。问题在于，重构和孤岛在物理上是耦合的——开关状态决定了网络拓扑，拓扑决定了哪些节点被隔离出来，隔离区域里的DG才能决定能不能形成自平衡孤岛。反过来，如果你想保留某个孤岛，那孤岛边界上的开关就不能动，这会直接影响重构方案能不能执行。

举个我在复现中反复用到的例子。33节点系统末端有一段长馈线，电压本身偏低。如果单纯重构，把所有失电负荷一次性转接到另一条馈线上，末端电压可能掉到0.95以下，这时候要么甩负荷，要么就得开关站里加无功补偿。但如果末端附近有一台光伏DG，把它划成孤岛，让DG带住最末端的几个负荷，转供馈线的负担立刻轻了，电压也抬回来了，整条系统的恢复率反而更高。

这就是“同时考虑”的价值：不是先算重构再补孤岛，也不是先划孤岛再看剩多少负荷能转供，而是把两个决策放进同一个优化模型里，让求解器自己去权衡。复现这个策略，最核心的工作其实不是写代码，而是把物理问题翻译成一组合适的数学约束。

1.3 复现目标与验证环境

我复现时用的是常说的33节点标准测试配电网算例，就是33-bus那个经典算例，节点0作为变电站出口，包含5条联络开关支路。故障设在其中一条主干馈线上，然后让优化程序给出一个包含“哪些开关合、哪些开关分、哪些DG做孤岛根、哪些负荷恢复供电”的完整方案。

程序跑完直接输出：恢复负荷比例、开关操作次数、孤岛数量、最低节点电压。整篇文章会从建模推导讲到Matlab代码实现，最后把我复现时踩过的坑都列出来。如果你是在校学生做故障恢复课题，或者正在做配网自愈/馈线自动化工程，这篇内容应该能帮你省下不少时间。

2. 统一优化模型的搭建思路

2.1 决策变量：三组0-1变量加上连续量

要把这个联合优化问题扔给求解器，首先得把“电网状态”表达成变量。我最终用了三组0-1变量加几组连续变量，全部放在一个MILP（混合整数线性规划）模型里。

第一组是节点变量，记为x(i)。x(i)=1表示节点i处于通电状态，对应负荷也恢复供电；x(i)=0表示该节点断电。注意这里的“节点”是配电网上每个母线节点，不只含负荷节点，因为很多中间节点没负荷，但必须通电才能串起下游。

第二组是支路变量，记为z(e)。z(e)=1表示支路e上的开关闭合，z(e)=0表示断开。支路集合里同时包含正常的分段开关支路和几条联络开关支路。故障支路直接从支路集合里剔除掉，或者在模型中强制锁死为断开。

第三组变量是这个方案里最特殊的，记为s(j)。s(j)=1表示节点j上的DG被选为“孤岛根”，也就是说，这个DG所在的那棵树，在拓扑上不连到变电站，而是由它独立供电。正常重构方案里s全为0，有孤岛方案时至少有一个s为1。

连续变量方面，每条支路上有两路功率变量P(e)、Q(e)，一个电压变量v(i)，还有一个专门用于生成辐射状拓扑的虚拟流变量F(e)。所有负荷、DG出力、支路阻抗参数都提前转换为标幺值，基准功率取10MVA，基准电压取12.66kV。这套标幺体系选得好，后面MILP的系数不会差太多量级，求解也稳定。

2.2 辐射状约束：闭合开关数、虚拟流与“根”

配电网和输电网最大的不同就是要求辐射状运行，也就是闭环设计、开环运行。故障恢复后，系统可能同时存在一个主网树和几个孤岛树，但整体必须仍然是一组互不环网的树。这是整个模型里最绕、也最容易写错的地方。

首先有一个拓扑学上的简单公式：闭合支路数 = 通电节点数 - 树的棵数。如果系统只有一个变电站根，没有孤岛，那闭合支路数就等于通电节点数减1。如果存在m个孤岛，那通电节点被分成了1+m棵树，闭合支路数就要再减m。用变量写出来就是：

sum(zvar) == sum(xvar) - 1 - sum(svar);

这条约束能保证“支路数对得上”，但它不保证连通。比如，我可以拿出三根支路接成一个三角形环，节点通电数4，支路数3，等式满足了，但环网出现了。所以必须再补一个“连通性”约束，也就是让每一棵树内部有路，树和树之间没有联系。

工程上最稳妥的做法是引入虚拟流。把每个通电节点想象成一张“报到单”的消费者，每棵树的根是报到单的生产者。每个通电节点必须恰好消费1份报到单，根节点生产若干份。只有闭合支路上才能传递报到单。这样，只要约束满足，系统就自动是连通且无环的。

这个虚拟流变量的约束写法，常见有两种。一种是按“节点净流出 = 1”写，所有非根节点都等于1；另一种用不等式加Big-M放开。我在代码里用的是双向不等式版本，因为处理断电节点（x=0）时更干净：

% F变量上限与支路闭合状态绑定 prob.Constraints.flowCap = Fvar <= nNode .* zvar; % 每个节点的虚拟净流出：消费1份，除非是根提供了额外供给 for i = 1:nNode inIdx = find(toVec == i); outIdx = find(fromVec == i); expr = sum(Fvar(inIdx)) - sum(Fvar(outIdx)); prob.Constraints.flowLow{i} = expr >= 1 - supply(i) - nNode .* (1 - xvar(i)); prob.Constraints.flowUpp{i} = expr <= 1 - supply(i) + nNode .* (1 - xvar(i)); end

supply(i)这个概念很关键。普通节点supply=0，它就必须净消费1份虚拟流。变电站根节点supply是一个范围在1到nNode之间的连续变量S0，代表变电站那棵树里有几个通电节点。DG孤岛根节点j的supply则和s(j)绑定：s(j)=1时supply可以是1到nNode之间的数，s(j)=0时supply强制为0。这样约束一闭合，模型自然就知道每个孤岛根养活了几个节点，不需要我们手工指定“孤岛范围”。这算是我复现里觉得最妙的一处设计。

2.3 潮流与电压的线性化写法

故障恢复本质上是个潮流问题，但MILP不能直接塞进非线性潮流方程。工程上常用的折中方案是DistFlow线性化。DistFlow是一种面向辐射状配电网的支路潮流方程，忽略网损高阶项后可以写成线性形式。

首先是有功和无功在每个节点的功率平衡。把每条支路的方向定好，比如支路e从节点from指向节点to，那么节点i的功率平衡就是：

% 关联矩阵A：A(i,e)=1表示支路e从i流出，-1表示流入 A * Pvar == (loadP - Pgen) .* xvar; A * Qvar == (loadQ - Qgen) .* xvar;

这里有个细节值得强调：DG出力在节点不通电时是不该存在的，所以DG出力必须乘以x(i)，而不能写成常数。我一开始偷懒把DG出力写成固定常数，结果求解器给出的方案里出现了“断电节点还在向网络注入功率”的荒谬结果，后来才反应过来必须乘上节点状态变量。

变电站根节点特殊处理。它本质上是松弛节点，有功无功注入可以变化，所以我把变电站节点的功率平衡单独写成注入变量Psub和Qsub，再加一个容量约束：

prob.Constraints.subBalanceP = A(subIdx,:) * Pvar == -Psub; prob.Constraints.subBalanceQ = A(subIdx,:) * Qvar == -Qsub; prob.Constraints.subLimit = Psub.^2 / 10 + Qsub.^2 / 10 <= 0.5; % 示例参数

电压方面，线性化DistFlow给出：如果支路e闭合，那么下游节点电压v_to近似等于上游节点电压v_from减去2(rP+xQ)。写成大M形式，闭合时等式生效，断开时约束自动放松：

for e = 1:nBranch f = fromVec(e); t = toVec(e); expr = vvar(t) - vvar(f) + 2 .* (r(e) .* Pvar(e) + x(e) .* Qvar(e)); prob.Constraints.volLow{e} = expr >= -Mv .* (1 - zvar(e)); prob.Constraints.volUpp{e} = expr <= Mv .* (1 - zvar(e)); end

Mv取多大很有讲究。电压差标幺值最多0.1左右，2(rP+xQ)的量级通常不超过0.05，所以我取Mv=2，既保证足够大不会误伤可行解，又不会大到引起数值病态。如果Mv取100、1000这种夸张值，intlinprog在分支定界时很容易出现虚假不可行，后面会详细讲。

2.4 目标函数与开关操作惩罚

目标函数直接决定恢复策略的偏好。我采用“恢复负荷最大化为主、开关操作次数最小化为辅”的写法：

% 最大化加权恢复负荷量，减去开关操作代价 obj = sum(loadP .* weight .* xvar) - lambda .* sum(opvar); prob.Objective = obj; prob.ObjectiveSense = 'max';

weight是节点负荷的重要程度，普通负荷取1，一级负荷可以取10、20。开关操作次数这里需要一点小技巧：故障前的开关状态z0是已知的，每个开关e的操作次数用绝对值 |z(e)-z0(e)| 表示，但绝对值不能直接进线性规划，所以引入辅助变量op(e)，用两条不等式linearize：

prob.Constraints.op1 = opvar >= zvar - z0vec; prob.Constraints.op2 = opvar >= z0vec - zvar;

lambda的取值也有讲究。如果lambda太大，求解器宁可少恢复部分负荷也不肯多扳一个开关，这不符合调度直觉；太小则开关操作次数可能失控。我在33节点算例里测试，lambda取0.1到0.3之间比较合适。比如恢复一个0.1MW负荷收益是0.1，操作一个开关惩罚0.1，这样一次额外操作必须有等量负荷收益才划得来。实际工程里如果更看重操作安全，可以把这个系数调大到0.5以上。

另外给个提醒，目标函数里负荷、开关操作惩罚的单位必须一致。负荷是标幺功率（MW），开关是次数，两者直接相加在量纲上不严格，但工程上大多数论文都这么干。如果你在写学术论文，建议把开关操作惩罚写成“相当于损失多少负荷”的语义，比如一个开关=0.1MW损失，这样审稿人不会挑刺。

2.5 孤岛其实是被“约束”自动刻画出来的

很多入门玩家看到“考虑孤岛”第一反应是：需要在约束里明确写出“哪些节点属于哪个孤岛”，然后给每个孤岛加功率平衡约束。这种思路不是不行，但候选孤岛集合会爆炸，而且孤岛边界和重构开关相互影响，非常难枚举。

我复现的这套模型绕开了这个问题。思路很直接：不预设任何孤岛集合，只让模型决定哪些DG节点当“根”。拓扑约束保证了每个通电节点最终都归属于某棵树，这棵树的根要么是变电站，要么是一个被选中为孤岛根的DG。树的根确定后，孤岛范围自然就是这棵树覆盖的节点集合。

功率平衡约束再兜底：如果某个DG被选为孤岛根，它那棵树上所有负荷的有功无功消耗，都必须由这台DG的出力兜住。DG出力不够？那对不起，这棵树上肯定有一些x(i)变量会变成0，也就是甩负荷。整个过程是全自动的，不需要额外告诉模型“哪里该划孤岛”。这就是我认为“同时考虑孤岛与重构”在建模层面最优雅的解法——孤岛不是特殊元素，它只是“根不在变电站的树”。

3. Matlab实现流程与核心代码

3.1 文件分工与数据准备

我习惯把代码分成四个文件，各管一段，调试起来不糊涂：

• restore_main.m：主脚本，负责组装模型、求解、输出结果。
• case33_data.m：返回算例参数结构体，包含节点、支路、负荷、DG接入点。
• build_restore_model.m：根据参数结构体生成优化问题对象。
• check_restore_result.m：对求解结果做一次独立校验，防止模型bug“静默出错”。

33节点算例的参数整理是第一步。节点参数包括有功负荷、无功负荷、权重；支路参数包括起始节点、终止节点、电阻、电抗、初始开关状态、是否为联络开关。DG我放在节点22和节点25，两个都是经典的光伏接入位置，额定出力分别取0.4MW和0.3MW，功率因数固定0.95。这个配置不是随便拍的，因为故障点设置在节点7到节点8之间时，节点22和节点25正好落在失电区，DG能不能顶上去直接决定了末端负荷能否保住。

故障支路的处理建议是：直接从候选支路表里删掉，而不是靠zvar强制为0。原因很简单，删掉可以减少一个整数变量，求解更快。程序里我会在读取支路表之后加一行过滤：

br = br(br.include == 1, :); % include=0的支路表示已故障隔离

3.2 核心约束的代码组装

用optimproblem建模比手搓intlinprog的A矩阵省心太多，自带的optimvar就能定义变量类型和边界。变量定义部分长这样：

xvar = optimvar('x', nNode, 1, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); zvar = optimvar('z', nBranch, 1, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); svar = optimvar('s', nDG, 1, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); opvar = optimvar('op', nBranch, 1, 'Type', 'continuous', 'LowerBound', 0); Pvar = optimvar('P', nBranch, 1, 'LowerBound', -Mpq, 'UpperBound', Mpq); Qvar = optimvar('Q', nBranch, 1, 'LowerBound', -Mpq, 'UpperBound', Mpq); Fvar = optimvar('F', nBranch, 1, 'LowerBound', 0, 'UpperBound', nNode); vvar = optimvar('v', nNode, 1, 'LowerBound', 0.95, 'UpperBound', 1.05); vvar0 = vvar(1); % 变电站根节点电压

这里svar的长度是DG数量，而不是节点数量。如果某个DG节点没被选为孤岛根，s=0；选上则s=1。把svar和节点编号的映射关系提前存好，后面加约束时直接用映射数组。

孤岛根与通电状态的联系我额外加了两条，避免出现“断开节点当孤岛根”的脏方案：

% 被选为孤岛根的DG节点必须通电 prob.Constraints.dgOn = svar <= xvar(dgNodeIdx); % DG根必须至少带一个节点（自己） prob.Constraints.dgSupplyLow = Svar >= svar;

供电变量Svar是连续变量，范围0到nNode。注意，Svar和svar绑定后，如果svar=0，Svar被迫等于0；如果svar=1，Svar可以在1到nNode之间浮动。它具体取多少，完全由虚拟流约束自动决定——它等于这台DG孤岛里通电节点的总数，不需要我们预先指定。

支路闭合与两端节点通电的一致性也别漏，否则会出现“开关合着但两端节点标成断电”的怪事：

prob.Constraints.zLinkFrom = zvar <= xvar(fromVec); prob.Constraints.zLinkTo = zvar <= xvar(toVec);

3.3 求解设置与结果校验

所有约束组装完成后，直接调用solve。Matlab自带的intlinprog对33节点这样的中小规模问题完全够用，不需要额外装商业求解器。求解器选项里我特别关注两个：TimeLimit必须设置，不然大算例可能跑一晚上不收敛；IntegerTolerance默认值比较松，建议压紧一档：

opts = optimoptions('intlinprog', ... 'Display', 'iter', ... 'TimeLimit', 120, ... 'IntegerTolerance', 1e-5, ... 'LPOptimalityTolerance', 1e-4); [sol, fval, exitflag] = solve(prob, 'Options', opts);

求解完成后，我强烈建议做一次“独立校验”，不要直接信MILP的结果。我自己写了一个check_restore_result.m，逻辑是：

1. 根据求解得到的z和s，对所有通电节点做DFS遍历，看是否能从变电站根或DG孤岛根到达全部通电节点。
2. 对每个孤岛根检查功率是否自平衡：孤岛内DG总出力是否不小于负荷总消耗。
3. 对每条闭合支路手动算一遍DistFlow电压降，看看和模型给的v是否接近。
4. 统计开关操作次数、恢复负荷比例。

这个校验脚本救了我好几次。有一次模型因为Big-M取值不当，求解器返回的“最优解”里孤岛根节点编号是错的，但目标函数值看起来还挺合理。要不是独立校验发现了树访问不到，我可能就把错误结果当成正确结果用了。

3.4 一个典型算例的输出结果

以节点7-8之间故障为例，我跑出来的最优方案具有很典型的两段式结构：

• 分段开关：断开故障点两侧的几个分段开关，隔离故障。
• 联络开关：合上联络开关18-33，使主网可以向失电区域中部转供。
• 孤岛划分：因为DG22所在分支离故障点近，转供路径上容量紧张，求解器选择让DG25在末端形成孤岛，带节点26、27、29、30、31、32、33的一部分负荷。

典型结果统计如下：

如果不允许孤岛，只做重构，同一故障下最优方案的恢复负荷比例大约只有87%，末端最低电压还压到0.95的边界上。两者的差距就是“同时考虑孤岛”带来的真实收益。换一句更直白的话：多一个孤岛决策，等于多一个供电电源，末端电压也会因为负荷被本地消化而变得更健康。

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

4.1 第一次求解就不可行

我踩的第一个坑就是模型一上来报“No feasible solution found”。排查思路是按层剥离：先把约束全部注释掉，只留变量边界和目标函数，确认问题框架能跑；然后逐步加约束，每加一组就重新求解，看到底是哪一组把可行域砍没了。

结果发现元凶就是树约束。sum(z) == sum(x) - 1 - sum(s)这条式子，如果两边单位对不上就会出问题，比如svar里面包含所有节点而不仅仅是DG节点，或者sum(x)把变电站根节点漏算了。这类问题光靠看代码很难发现，最快的方法是为一个极简单的小网络单独测拓扑约束：比如3节点2支路、一个变电站、一个DG，穷举所有开关状态，验证约束放行和禁止的状态集合是否和你物理直觉一致。

第二个常见的不可行原因是变电站注入容量太紧。如果Psub上限设得不够大，而负荷权重分布又迫使求解器必须恢复某块区域的负荷，模型也会无解。这时候适当放松变电站容量约束，看恢复方案会不会变，就能定位问题。

4.2 Big-M的病态与符号反了

大M法最大的敌人是自己。M值太大，求解器预处理的边界质量变差，分支定界树膨胀，求解时间成倍增加，甚至出现虚假不可行；M值太小，又会把真实可行解切掉。

我的经验法则：功率类大M用该支路可能出现的最大功率的2到3倍，电压类大M用1到2就够了。以33节点算例为例，基准功率10MVA，最大支路功率撑死不超过0.6pu，所以我取Mpq=5已经是富余了。之前试过Mpq=100，求解时间从几秒涨到了几分钟，而且低电压节点数明显变多——这就是数值病态的典型症状。

电压约束的符号反了也是重灾区。DistFlow的线性化是v_to = v_from - 2(rP+xQ)，如果代码里写成v_to = v_from + 2(rP+xQ)，那么求解器为了最小化电压降会故意把支路功率方向设置成反的，结果就是恢复方案看起来“全对”，但实际潮流算出来电压全红。这里有个很笨但很有效的自检办法：把网络设成单馈线简单辐射状，所有开关全合、全负荷投入，然后把模型退化成纯线性潮流，对比手算或商业潮流软件的结果，几秒钟就能看出符号问题。

4.3 孤岛划出来但负荷带不动

如果求解结果里出现了孤岛根，但看功率数据，孤岛内的总负荷明显超过DG出力，那大概率不是模型的问题，而是目标函数里负荷权重的锅。比如我把二级负荷的权重设成10，一级负荷设成1，求解器有可能为了恢复一个二级大负荷而让DG过载，但功率平衡约束是死的，它过载不了，所以更常见的情况是：某个重要大负荷被恢复了，但孤岛内其他普通负荷全被牺牲掉，看起来就像“孤岛根带不动”。

优化器没有感情，它只会按权重干活。如果你希望孤岛内的负荷尽量均衡，就在目标函数里把权重差异调小；如果你更看中保障某几个重要用户，那就要接受周围普通负荷被甩掉。这个取舍本质上是一个调度策略问题，不是bug。我在复现时把节点18、22、33这几个相对重要的负荷权重设为2，其余都设为1，跑出来的方案就比较均衡，孤岛内不会出现“保了一家大工厂、整个村子黑到底”的奇葩结果。

4.4 求解速度与求解器选型

33节点跑几秒，但如果你换到69节点或者118节点，Matlab自带的intlinprog可能就要几十秒甚至几分钟了。这时候有几个立竿见影的加速手段：

• 把故障支路直接从候选集合删掉，而不是留着当0-1变量。
• 把那些不存在DG的节点上的s变量从模型里彻底移除。s变量的数量对求解复杂度影响很大，毕竟每个s都可能引入一棵新树。
• 给x变量加合理的上界修剪：如果某个节点在故障隔离后没有任何供电路径，它的x直接置0，而不是交给求解器去试。
• 用TimeLimit限制求解时间，并让intlinprog返回当前最好可行解。工程上“次优但可行”通常比“最优但跑不完”更有价值。

如果算例规模更大，或者你想做在线滚动优化，可以考虑换商业MILP求解器。Matlab的optimproblem接口可以直接设置求解器，不需要改模型代码，只是不同求解器对Big-M的容忍度不太一样，换完之后建议重新调一遍M值和求解器容差。

4.5 给后来者的调参建议

复盘整个复现过程，我认为最值得花时间打磨的不是代码本身，而是一组参数的质量：大M值、DG位置与容量、负荷权重、开关操作惩罚系数。这些参数决定了模型的“性格”，同样的代码，参数不同，恢复策略完全不同。

我给后来者三条非常具体的建议：

第一，先把“无孤岛纯重构”跑通，再把s变量加进去。分步调试能让你在出问题时清楚地知道是重构约束的问题还是孤岛约束的问题。第二，完成一次求解后，一定要用独立潮流校验环节确认结果，不要迷信求解器返回的“最优”二字。MILP的数值解在有大量Big-M约束时可能出现轻微违背，独立校验能把这种风险暴露出来。第三，把故障场景参数化，多测几个不同位置的故障，不要只盯着一个场景调参数。有些参数组合在故障点A上表现很好，换到故障点B就完全不行，这种鲁棒性测试在工程交付中是躲不掉的。

最后分享一个我自己的习惯：我通常会在跑完所有算例后，把每个故障场景的结果表格导出来，包括恢复率、操作次数、孤岛数量、计算时间，然后对着表格调整策略参数。数据摆在面前的时候，很多原本模糊的取舍都会变得一目了然。这个项目本身是个很不错的起点，如果你打算继续往深做，可以考虑把网损迭代、分布式光伏出力不确定性、储能参与黑启动这些因素逐步加进去，每一步都能让策略离工程落地更近一点。