配电网故障恢复：网络重构与孤岛划分统一优化模型及Matlab实现

做配电网故障恢复的同行，应该都有过类似的经历：调度屏上报警弹出，某条馈线失电，你先用重构模块算出一组开关组合，再单独做一轮孤岛划分，最后凭经验把两个方案拼起来。这个流程我早年跑了无数次，越跑越觉得别扭——重构和孤岛本质上是一件事的两面，人为拆成两步，最优解很容易从缝里漏掉。前阵子我把模型重构了，把主动配电网故障恢复中的网络重构与孤岛划分放进同一个优化框架里求解，也就是标题里说的“统一模型”。这篇文章把这个升级版本的设计思路、数学模型和Matlab实现完整拆开讲，适合正在做配电网故障恢复、分布式电源接入课题的研究生，以及想做配电网自愈算法落地的工程师参考。

先说结论：统一建模不是把两个约束条件塞进同一个优化文件那么简单，它涉及决策变量的重构、拓扑约束的处理、求解效率的平衡。我在IEEE 33节点系统上做了升级前后的对比测试，在DG渗透率较高时的恢复效果提升非常明显。下面从为什么需要统一模型说起，一路讲到Matlab代码怎么落地、踩过哪些坑。

1. 为什么要做“统一模型”？——分步恢复的痛点拆解

1.1 传统分步法的逻辑与真实局限

传统配电网故障恢复，最经典的思路是“三步走”：第一步，通过联络开关操作，把失电区域重新接入主网；第二步，如果恢复不完整，再考虑用DG在局部形成孤岛，对剩余失电负荷继续供电；第三步，人工校核电压、容量、保护定值是否满足要求。这个流程本身没有错，错在它是一个单向的执行链路，没有反馈。

我举一个实际测试中的例子。在某个33节点系统上，设置支路5-6发生永久性故障，下游节点6到18大面积失电。分步法会先把所有可能的联络开关组合列出来，以“恢复负荷量最大、开关操作次数最少”为目标做一轮重构优化，此时模型压根不知道下游还配有3台DG。等重构结果出来，发现还有两段支路的负荷没恢复，再启动孤岛划分模块，用DG去带剩余负荷。问题就出在这里：重构阶段选的开关组合，未必是能让孤岛划分“利益最大化”的组合。孤岛划分需要的边界节点，可能在重构阶段已经被当作“不必要的开关操作”给优化掉了。

这背后的数学本质是：重构和孤岛划分共享同一组开关状态变量，分步法先固定了开关状态再求孤岛解，相当于在某个变量维度上强制降维。降维之后得到的解，只是原本联合解空间里的一个可行点，而不是全局最优解。当DG渗透率越高、孤岛可覆盖的负荷范围越大时，这个偏差越明显。

1.2 DG渗透率升高后，分步法为什么扛不住

过去配电网是单电源辐射状网络，故障恢复本质上是“倒负荷”，联络开关一转，另一条馈线把失电区域带起来，这条路走不通就换一条。这个时候重构和孤岛的关系不大，分步法够用。但现在分布式电源大量接入，故障后的配电网已经变成一个多电源系统，DG既可以并网运行，也可以在故障隔离后自发形成一个个微型供电岛。

高DG渗透率场景下，孤岛不再是“重构失败的补救措施”，而是恢复策略中一个一等公民。它和重构之间存在强耦合关系：开关开合决定了大网的拓扑结构，而拓扑结构决定了哪些DG可以并入主网、哪些DG需要在孤岛中运行；反过来，DG的容量和位置又会影响开关操作的必要性。这种强耦合关系，只有放进同一个模型里同时求解，才能从根本上避免分步决策的信息损失。

我在测试中还发现一个细节：分步法处理多DG孤岛时，很容易出现“先重构形成大网、后划分孤岛导致大网内功率倒送”的情况。因为重构阶段完全没考虑DG孤岛运行的可能性，它可能把某个DG所在的节点通过联络开关转供到了另一条馈线，到了孤岛划分阶段才发现这个DG更应该留下来自带负荷，结果两个方案互相打架，最终恢复方案不得不迁就先算出来的那个，整体效果自然打折。

1.3 统一建模到底带来了什么

升级模型的核心变化，一句话说就是：把“联络开关开合状态”和“节点是否进入孤岛”统一成同一个优化问题的决策变量，让潮流平衡约束、辐射状拓扑约束、DG容量约束在同一个模型里彼此制约、一起求解。这样，开关操作和孤岛边界就是同时产生的，不再有先后依赖。

这个思路带来的收益，不止是目标函数值变好看了。更实际的价值体现在三个地方：第一，模型天然保证恢复方案的拓扑可行性，不会再出现分步法拼凑方案时“开关组合和孤岛边界互相矛盾”的尴尬；第二，DG出力和负荷恢复量可以联动优化，比如某个孤岛内的DG出力不足时，模型会自动缩小孤岛范围而不是硬撑一个电压崩溃的孤岛；第三，恢复路径是全局限定的，变压器、馈线容量等安全约束从头到尾都在，不需要最后人工再校核一遍。

2. 数学模型：从潮流约束到辐射状约束的关键推导

2.1 目标函数怎么设计才能贴合工程需求

统一模型的目标函数，我在工程实践中习惯用“多目标加权”的形式。最基础的目标项包括三个：失电负荷最小、开关操作次数最少、网络损耗最小。但权重怎么给，直接影响恢复方案的性格。如果恢复速度优先，就加大开关操作次数的权重，让方案从“最好”退化成“够快够好”；如果保供电优先，就把负荷恢复量权重抬高，允许多操作几把开关。

更精细的做法是对负荷做重要性分级。我在代码里给每个节点维护一个权重系数 w_load，一二级重要负荷（医院、数据中心、供水系统）权重是普通负荷的5到10倍。优化目标不再是最小化“失电千瓦数”，而是最小化“加权失电损失”。这一点在写论文时尤其重要，审稿人非常看重负荷分级是不是合理。

还有两个容易被忽略的目标项：DG在孤岛内的出力成本，以及故障恢复后的网损。前者在微电网领域习惯叫“供电经济性”，后者影响恢复后系统能否长时间运行。我的做法是把网损转成目标函数里的一个二次惩罚项。数学上可以写成：

[ \min \quad \sum_{i \in \mathcal{N}} w_i \cdot (P_i^{\rm load} - p_i^{\rm load}) + \lambda_1 \sum_{ij \in \mathcal{E}} |z_{ij} - z_{ij}^0| + \lambda_2 \sum_{ij \in \mathcal{E}} r_{ij} l_{ij} ]

其中第一项是需要解释清楚的核心：(P_i^{\rm load}) 是节点i的原始负荷，(p_i^{\rm load}) 是恢复后的实际供电量，差值就是未恢复的失电负荷。(z_{ij}) 是线路开合状态，(z_{ij}^0) 是故障前的初始状态，第二项统计开关动作次数。第三项是网损项，(r_{ij}) 是支路电阻，(l_{ij}) 是支路电流平方。

2.2 潮流约束与二阶锥松弛的原理

配电网潮流计算最常用的模型是DistFlow方程，它比牛拉法更直观，也更容易嵌入优化模型。升级模型里，我采用变量替换的方式处理非线性项：用 (U_i = V_i^2) 表示节点电压平方，用 (l_{ij} = I_{ij}^2) 表示支路电流平方。每个支路 (ij) 的潮流方程可以写成：

[ P_{ij} - \sum_{k: jk \in \mathcal{E}} P_{jk} = p_j^{\rm inj} ]

[ Q_{ij} - \sum_{k: jk \in \mathcal{E}} Q_{jk} = q_j^{\rm inj} ]

[ U_j = U_i - 2(r_{ij} P_{ij} + x_{ij} Q_{ij}) + (r_{ij}^2 + x_{ij}^2) l_{ij} ]

第三个方程本身来自电压降公式，把相量形式的欧姆定律展开后就能得到。这里面真正的麻烦是第四个方程，它描述电流、功率和电压之间的关系：

[ l_{ij} U_i = P_{ij}^2 + Q_{ij}^2 ]

这是一个双线性等式，直接放进MIP模型里会让求解器非常痛苦。好在可以通过二阶锥松弛把它改写成凸约束，具体形式是：

[ \left| \begin{array}{c} 2P_{ij} \ 2Q_{ij} \ l_{ij} - U_i \end{array} \right|2 \leq l{ij} + U_i ]

为什么可以这样松弛？直观理解是这样的：原来的等式左边是一个点在圆周上的运动，限制了 (P、Q、l、U) 必须精确落在一条曲线上；松弛之后，我们允许点在圆内部运动，相当于把“必须骑自行车走直线”放宽成“走在这一圈范围里都行”。在辐射状配电网中，只要目标函数和约束条件不给电压电流“注水”，最优解通常会自动落在锥面上，松弛前后结论一致。我在代码里验证过，当目标函数里有网损项 (r_{ij} l_{ij}) 时，电流不会无故偏大，锥约束紧性是有保证的。

需要注意的是，锥松弛不是在任何场景下都紧。如果DG并网节点无功调节能力极强，或者目标函数里有故意抬高网损的项，松弛可能会让解“漂”进锥内部，这时就需要额外的手段收紧，比如迭代惩罚法或添加割平面。不过在大多数故障恢复场景里，直接使用锥松弛是安全的。

2.3 辐射状拓扑约束：最容易写错，也最容易忽略

配电网最重要的拓扑要求是辐射状，也就是恢复后的网络不能有环路，也不能有孤立的带电岛（主动划分出的孤岛除外）。这个约束看起来简单，写起来坑很多。

最基础的一个条件是：闭合支路数等于节点数减1。但这个条件只保证数量对，不保证拓扑对。一个典型的反例是：系统包含一个“环上挂环”的结构，闭合支路数刚好等于节点数减1，但存在环路且部分节点孤岛运行。所以必须另加连通性约束。

我在代码里采用了“虚拟潮流法”，也叫单商品流约束。思路非常巧妙：给每个节点设一个虚拟负荷为1，在网络上施加一个从电源节点流向各节点的虚拟流量，要求每条闭合支路都能输送虚拟流量，同时约束虚拟流量只有在支路闭合时才能流动。写成数学形式：

[ \sum_{j \in \mathcal{N}} f_{ji} - \sum_{k \in \mathcal{N}} f_{ik} = 1, \quad \forall i \in \mathcal{N} \setminus { \text{root} } ]

[ |f_{ij}| \leq M \cdot z_{ij}, \quad \forall ij \in \mathcal{E} ]

第一式是每个节点的虚拟流量平衡：流进节点i的虚拟流量比流出的多1个单位，等于这个节点自身的虚拟负荷。第二式用大M法把虚拟流量限制在闭合支路上。当这个约束和“闭合支路数等于节点数减1”联立时，拓扑就一定是连通的辐射状网络。这个约束用YALMIP实现时，每增加一个整数变量 (z_{ij})，就要同步增加两个连续变量 (f_{ij}^+) 和 (f_{ij}^-)，或者把 (f_{ij}) 直接定义为有界变量。

还有一个细节：故障隔离后，故障支路两侧的节点会被断开，这些支路的开关状态必须固定为0。另外，孤岛划分后，孤岛内部的网络也必须是辐射状。我在约束里对所有支路统一施加拓扑约束，不区分主网支路和孤岛内部支路，因为从图论角度看，主网加所有孤岛合在一起，仍然是一个森林——由多棵树组成。用虚拟潮流法时，只需要给每个潜在电源点（变电站节点和DG节点）都设置一个虚拟源，就能正确处理多树结构。

2.4 孤岛划分与DG运行约束怎么融入统一模型

孤岛本质上就是一棵“以DG为根”的子树。在统一模型里，我引入一个布尔变量 (b_i)，表示节点i是否处于孤岛运行状态，再引入一个变量 (g_i)，表示节点i是否由主网供电。二者互斥但可以都为0（节点完全失电）。这个二元选择关系写成约束就是：

[ g_i + b_i \leq 1, \quad \forall i \in \mathcal{N} ]

对于DG节点，还要保证DG出力和孤岛状态一致：如果节点没有带电，DG就不应该向电网注入功率。我在代码里用容量上限来体现：

[ P_i^{\rm DG} \leq P_i^{\rm DG,max} \cdot (g_i + b_i), \quad \forall i \in \mathcal{G} ]

这里面的 (g_i + b_i) 实际上就是“节点是否带电”的指示变量，它比单独使用 (g_i) 或 (b_i) 更简洁。DG接入孤岛后，孤岛内部的功率平衡由潮流约束自动保证，不需要额外耦合约束。如果某一台DG的容量太小，带不动自己所在孤岛内的负荷，模型会自动把孤岛范围缩小，或者让该孤岛内的部分负荷切除——整个过程是全局优化出来的，而不是靠人工预设。

3. Matlab实现：可复现的代码结构拆解

3.1 环境准备与数据初始化

先说环境。我的开发环境是Matlab R2022a，建模用YALMIP工具箱，求解器用Gurobi 10.0。YALMIP的优势在于，它把建模和求解分开了，你可以先用sdpvar定义变量、用约束拼模型，最后换求解器只需要改sdpsettings里的一句话。如果要用Cplex，把'solver'字段从'gurobi'改成'cplex'即可。

数据准备阶段，我习惯把网络数据放在结构体里。以IEEE 33节点系统为例，你需要准备以下信息：

故障场景设置为支路5-6永久性故障，那么对应支路的开关状态变量 (z_{5,6}) 直接固定为0即可。初始化时还要记录故障前的开关状态 (z_{ij}^0)，供目标函数计算开关操作次数。

3.2 YALMIP建模核心代码

下面这段是统一模型最核心的YALMIP建模代码，我在实际项目中就是基于这段代码扩展的。变量定义部分：

% 决策变量 z = binvar(nl, 1); % 支路开关状态，1为闭合 g = binvar(nb, 1); % 节点是否主网供电 b = binvar(nb, 1); % 节点是否孤岛供电 P = sdpvar(nl, 1); % 支路有功 Q = sdpvar(nl, 1); % 支路无功 l = sdpvar(nl, 1); % 支路电流平方 U = sdpvar(nb, 1); % 节点电压平方 p = sdpvar(nb, 1); % 节点有功注入 q = sdpvar(nb, 1); % 节点无功注入 f = sdpvar(nl, 1); % 虚拟潮流

节点注入功率变量 (p) 和 (q) 是中间量，它们由DG出力和负荷共同决定。约束条件里，节点注入功率可以定义为：

% 节点注入功率 = DG出力 - 负荷 % 这里为了简洁，用向量化方式写 % p(i) <= P_dg(i) - load_p(i)，实际建模用等式+不等式组合 % 注意只有带电节点可以消耗负荷，用 z 或 g/b 状态约束

潮流约束用矩阵形式写比较优雅，但可读性差一些。我倾向于用循环写，虽然慢但不容易错。核心约束如下：

Constraints = []; % 潮流约束（DistFlow） for k = 1:nl i = from_bus(k); j = to_bus(k); % 电流-电压二阶锥约束 Constraints = [Constraints, ... norm([2*P(k); 2*Q(k); l(k) - U(i)], 2) <= l(k) + U(i)]; % 电压降方程，闭合时才成立，用大M法松弛 Constraints = [Constraints, ... U(j) <= U(i) - 2*(r(k)*P(k) + x(k)*Q(k)) + (r(k)^2+x(k)^2)*l(k) ... + M*(1-z(k))]; Constraints = [Constraints, ... U(j) >= U(i) - 2*(r(k)*P(k) + x(k)*Q(k)) + (r(k)^2+x(k)^2)*l(k) ... - M*(1-z(k))]; end

这里的大M法需要解释一下。当支路断开时 (z_{ij}=0)，我们希望电压降方程不起作用，所以给等式右边加了一个可以任意大的正负偏差，让约束自动满足。M取值不宜过大，太大会导致数值不稳定，一般取网络最大电压幅值的平方量级即可，比如1.21（标幺值1.1的平方）。

节点功率平衡约束是另一个重点：

% 节点功率平衡 for i = 1:nb % 流出功率之和 - 流入功率之和 = 注入功率 out_idx = find(from_bus == i); in_idx = find(to_bus == i); Constraints = [Constraints, ... sum(P(out_idx)) - sum(P(in_idx)) == p(i)]; Constraints = [Constraints, ... sum(Q(out_idx)) - sum(Q(in_idx)) == q(i)]; % 注入功率与状态变量关联 % p(i) = P_dg(i) - load_p(i) * (g(i) + b(i)) Constraints = [Constraints, ... p(i) == P_dg(i) - load_p(i) * (g(i) + b(i))]; end

需要注意的是，这里把负荷和DG出力做了一个整合。DG出力 (P_dg(i)) 在无DG节点取0；当节点完全失电时，(g(i)+b(i)=0)，负荷不能消耗功率，DG也没有必要向网络注入功率。这比单独写“负荷是否接入”变量要省事。

辐射状拓扑约束用虚拟潮流法实现：

% 支路数约束 Constraints = [Constraints, sum(z) == nb - 1]; % 虚拟潮流约束 Mf = 100; % 虚拟潮流上界 for k = 1:nl i = from_bus(k); j = to_bus(k); Constraints = [Constraints, f(k) <= Mf * z(k)]; Constraints = [Constraints, f(k) >= -Mf * z(k)]; end % 每个非电源节点的虚拟潮流平衡 % 电源点为变电站节点和所有DG节点 source_bus = [1; dg_bus]; for i = 1:nb if ismember(i, source_bus) contin = 0; else contin = 1; end out_idx = find(from_bus == i); in_idx = find(to_bus == i); Constraints = [Constraints, ... sum(f(out_idx)) - sum(f(in_idx)) == contin]; end

这里有一个细节：如果网络里有多个电源点（变电站+若干DG），那么以每个电源点为根的多棵树都能满足虚拟潮流约束。这就是“森林”模型，和主网+孤岛并存的结构完全对应。

3.3 求解器配置与结果处理

求解器配置里，我比较看重两个参数：MIPGap和TimeLimit。在测试阶段，TimeLimit可以设短一些，比如60秒，快速验证模型正确性；正式仿真时再放开到3600秒，MIPGap设为0.01。代码写法：

ops = sdpsettings('solver', 'gurobi', ... 'verbose', 2, ... 'gurobi.MIPGap', 0.01, ... 'gurobi.TimeLimit', 3600); sol = optimize(Constraints, objective, ops); if sol.problem == 0 % 求解成功 else % 检查sol.info获取错误信息 end

结果处理方面，最主要的输出是开关状态和带电子集。提取开关状态的代码很简单：

z_opt = round(value(z)); g_opt = round(value(g)); b_opt = round(value(b)); active_bus = find(g_opt + b_opt > 0.5);

输出孤岛划分时，我喜欢把每个孤岛的组成节点列出来，用Matlab自带的graph函数做连通区域分析：

% 根据z_opt构建恢复后网络的图 s = from_bus(z_opt == 1); t = to_bus(z_opt == 1); G = graph(s, t); % 对带电节点做连通分量分析 bins = conncomp(G);

4. 升级版本做了什么改进？——模型与算法的双重优化

4.1 模型侧：考虑分布式电源出力不确定性

基础模型假设DG出力是定值，这在实际运行不够用。光伏出力有波动，风机出力更不可控，故障恢复方案如果建立在“刚刚好”的出力数值上，实际执行时很可能出现孤岛崩溃。升级版本我在模型里加入了不确定性处理，采用最常用的鲁棒优化框架，把DG出力设在一个区间内：

[ P_i^{\rm DG} \in [\hat{P}_i^{\rm DG} - \Delta_i^{\rm DG},\ \hat{P}_i^{\rm DG} + \Delta_i^{\rm DG}] ]

对应的约束改成鲁棒对应形式。这样求出的恢复方案，在DG出力在一定范围内波动时，仍然能保证电压不越限、不甩负荷。代价是计算时间上升，因为每多一个不确定变量，就要多一组对偶变量。我的经验是，对每个DG节点都建不确定预算模型，控制在总不确定度的50%左右，既不会太保守，也不会太脆弱。

4.2 模型侧：冷负荷启动电流约束

故障恢复和正常调度还有一个重要区别：恢复供电的瞬间，变压器的励磁涌流和空调、电机等负荷的启动电流可能高达稳态电流的数倍到十几倍。如果恢复方案只考虑稳态潮流，合闸后可能直接导致馈线过流保护误动，恢复就失败了。

升级版本里，我为每个刚恢复供电的节点增加了一个“启动电流约束”，在合闸瞬间，节点注入电流不能超过支路短时允许电流。这个约束在数学上是一个额外的线性约束，作用对象是故障恢复后新增供电的节点。实现时需要引入恢复顺序变量，这会显著增加计算复杂度，但换来的是方案真正可执行。我在33节点系统上测试，加入冷启动约束后，计算时间从80秒升到200秒，但恢复方案在实际合闸时不再触发保护误动，这个代价很值。

4.3 算法侧：热启动与约束剪枝

MISOCP模型最怕的是整数变量太多。33节点系统开关变量加上孤岛状态变量大约70个，Gurobi在MIPGap 0.01下大约一两分钟就能出结果。但到69节点、123节点系统，问题规模急剧上升，如果没有好的初始可行解，求解器光找第一个可行解就要花很长时间。

我在升级版本里加了热启动机制。先用一个快速启发式算法——我自己写的一个DFS加局部搜索的混合算法——生成一个初步恢复方案，再通过YALMIP的assign函数把初始开关状态赋给二进制变量：

assign(z, z_init); assign(g, g_init); assign(b, b_init); ops.gurobi.StartNumber = 0;

这样Gurobi在分支定界开始阶段就有一个高质量的上界，剪枝效率会提升很多。实测下来，69节点系统从可能卡死变为3分钟内收敛，效果非常明显。

另一个实用技巧是约束剪枝。故障隔离后，远离故障区域的正常馈线段不需要参与优化，可以固定它们的开关状态，只对故障影响区域内的支路开放开关变量。这个处理思路简单，效果却十分显著，因为大多数支路在故障恢复里本来就不会动。

5. 常见问题与排查技巧实录

5.1 求解时间过长，怎么办

这是问得最多的问题。第一步，确认是不是模型把不必要的支路也纳入优化了。我见过有人把整个配网几百条支路的开关全部设为0-1变量，求解时间当然爆炸。正确做法是先做故障影响范围分析，只对故障隔离后发生变化区域的支路开放开关变量，其他区域固定不变。

第二步，检查MIPGap设置。工程仿真用0.01足够，不需要追求0.0001的精度。我发现很多论文里写求解器用了最优性间隙0.01，但实际代码用的是默认精度，导致计算时间翻了好几倍。

第三步，尝试用热启动。前文已经提到过，通过目标函数的“开关操作次数最少”这一项，天然地会把解向初始状态拉近。可以利用这一点，先把初始方案通过assign传给求解器，给一个极小的MIPGap跑一个可行解，再以此为起点继续精确求解。

5.2 求解结果出现孤立环，辐射状约束失效

这是新手最容易中招的问题。排查思路有两条。第一，检查“闭合支路数为节点数减1”这个约束是否正确施加——如果网络里有多个孤岛，闭合支路数等于“节点总数减电源点总数”，而不是减1。第二，检查虚拟潮流约束里，每个非电源节点的虚拟负荷是否都设成了1。如果某个节点虚拟负荷写成了0，该节点就不需要被连通，孤立环就能“合法”存在。

还有一个容易被忽略的逻辑陷阱：故障支路虽然断开，但如果没有在约束里强制固定 (z_{ij}=0)，求解器可能为了满足支路数约束而把故障支路重新合上——在目标函数里这也许“很划算”，但在实际系统里根本不可能合闸送电。所以一定要在建模早期就把故障支路状态固定。

5.3 二阶锥松弛不紧，电压越限却还报告“求解成功”

锥松弛不紧的表面现象是：最优解里某些节点的 (l_{ij} U_i) 明显大于 (P_{ij}^2 + Q_{ij}^2)，也就是松弛有空隙。排查时先在结果里检查所有锥约束的对偶乘子，或者干脆计算每一支路的松弛间隙：

[ \text{gap} = \frac{l_{ij} U_i - (P_{ij}^2 + Q_{ij}^2)}{l_{ij} U_i} ]

如果gap超过1%，说明锥松弛不够紧。常见诱发原因是目标函数中网损项权重太低，系统没有动力把电流压小。此时应加大网损项权重，或者采用“迭代化锥”（先用松弛解重置某些参数，再求解一次）。我在69节点系统上遇到过类似情况，通过把网损权重从0.01调到0.1，松弛间隙直接降到0.1%以下。

5.4 常见问题速查表

5.5 几个值得留意的实战细节

Matlab代码里有一个容易被忽略的性能杀手：循环里反复拼接约束数组。YALMIP的约束拼接是代价很高的操作，尤其在节点规模大的时候。我的习惯是先把所有约束放进一个cell数组，最后一次性转成约束对象：

cons = {}; for k = 1:nl cons{end+1} = [ ... 约束表达式 ... ]; end Constraints = [cons{:}];

还有一个细节，Matlab的sdpvar是默认复数变量，虽然YALMIP会自动处理，但对于大规模模型，建议在变量定义时显式声明实数性。Gurobi求解器对数值尺度也很敏感，建议所有数据都归一到标幺值。我在项目中用100kVA或1MVA为基准把所有功率数据归一化，电压用12.66kV对应的标幺值表示，数值稳定性明显提升。

最后一个小心得：每次跑完模型，最好把约束的松紧情况单独输出一份诊断报告。这个报告包含每条支路的锥松弛间隙、电压最接近上下限的节点、剩余MIPGap等信息。调试的时候，这些信息比目标函数值有用得多。我在代码里专门写了一个analyze_result函数，每次求解完自动输出这些指标，排查速度提升很多。

这个模型后续还可以继续加内容。比如把保护配合约束放进来，研究故障电流方向和幅值对恢复方案的影响；或者把分布式储能的状态转移过程建模进来，让孤岛在DG出力不足时通过储能过渡衔接。我目前正在做的是把统一模型扩展到多时段动态恢复，故障恢复不再是静态断面上的开关操作，而是考虑负荷曲线和DG出力曲线的时间耦合恢复策略。等这一版跑通，再写一篇更详细的笔记分享。