1. 为什么故障重构让很多研究者在第一步就卡壳先说说我自己的经历。早几年刚开始接触主动配电网故障重构时我在网上找了一圈资料发现绝大多数文章都在讲理论推导一上来就是配电网潮流方程、混合整数非线性规划看得懂每个字合在一起就不知道该怎么落地。真正动手写MATLAB代码时更是完全不知道从哪一行开始。后来做的时间长了回头再看发现这个方向的门槛根本不在数学公式上而在于你还没想明白三件事问题到底长什么样、为什么非用二阶锥松弛不可、以及“求解器告诉你不可行”之后该往哪个方向排查。这篇博文就是围绕这三件事展开。我会用一套完整的MATLAB实现——从配电网故障场景建模、二阶锥松弛的核心推导到YALMIP调用商业求解器求解再到结果的可视化出图——把整个过程拆开讲清楚。无论你是刚入门的硕士生还是已经在做配电网优化方向的工程师这篇文章的目标都只有一个让你照着写能跑通跑通了以后敢改敢调。先说结论层面的事二阶锥松弛Second-Order Cone RelaxationSOCR在主动配电网故障重构中的应用本质上就是解决一个“非凸找最优”的问题。配电网在正常运行时是辐射状网络但故障发生后为了恢复失电区域的供电调度人员要决定哪些联络开关闭合、哪些分段开关断开这个决策过程天然带有0-1离散变量。与此同时配电网的潮流方程本身是非线性的也就是说功率、电压、电流之间存在乘积关系这又引入了连续层面的非凸性。“离散非凸”叠加在一起就是标准的混合整数非线性规划MINLP这种问题在数学上很难得到全局最优解。业界和学术界比较统一的思路是把这个硬骨头拆成两步啃先用二阶锥松弛把非凸的潮流方程变成凸约束再用Big-M法处理离散变量和连续变量之间的耦合。这样一来原来的MINLP就变成了混合整数二阶锥规划MISOCP商业求解器比如Gurobi、Cplex、Mosek可以直接处理并且能保证在有限时间内找到全局最优解。这套思路不是某个团队独有的魔法而是过去十年主动配电网优化领域的主流范式。我个人最看重的其实是整个过程里“可视化”的价值。算完一堆状态变量之后如果没有图你很难直观地看到“哪条线路失电了”“重构后哪些节点的电压恢复了”“DG出力被削减了多少”。做工程和研究可视化不是锦上添花而是你判断结果是否合理的第一道关卡。下面我按实际上手顺序把整套方案拆成五个部分讲问题建模、数学原理、代码实现、可视化、避坑经验。每个部分我都会尽量给出能直接用的东西。2. 主动配电网故障重构的问题本质一个“开关组合潮流寻优”的双层结构2.1 主动配电网与传统配电网的差异在哪里传统配电网是无源的电能从变电站单向流向负荷调度员能控制的手段非常有限——无非是投切电容器、调节有载调压变压器。但主动配电网Active Distribution NetworkADN不一样它接入了大量分布式电源DG、储能系统和柔性负荷这些资源既可以向电网供电也可以参与调节。说白了传统配电网的故障重构是“被动恢复”——把故障隔离掉合上联络开关让电网自己重新带上负荷。而主动配电网的故障重构是“主动优化”——不仅要恢复供电还要在恢复的过程中同时考虑DG出力的调整、储能是否放电、哪些可中断负荷可以暂时切除甚至要比较“多恢复一个节点但网损变大”和“少恢复一个节点但网损更优”之间的取舍。这就是为什么题目里特别强调“主动配电网”而不是简单说“配电网”。在主动配电网的故障重构中DG出力往往被视为决策变量而不是恒定不变的负荷注入。这带来一个问题DG的接入位置、出力上限、运行方式恒功率因数或恒电压都会直接影响重构结果。2.2 故障重构的决策变量与目标函数从数学上看故障重构要回答两个问题第一哪些开关应该闭合哪些应该断开第二在确定的网络拓扑下各节点注入功率是多少。先说拓扑层的变量。配电网是由分段开关和联络开关构成的网状结构正常运行时通过断开联络开关保持辐射状。故障时某个或某几个分段开关因故障需要断开原本正常的供电路径被切断这时就要闭合一些联络开关来转供负荷。所以决策变量包括每条支路的“开/关”状态通常用0-1变量表示。再说运行层的变量。支路有功、无功、节点电压幅值平方、支路电流平方这些是连续变量。DG的有功出力在某些模型里也是连续决策变量。目标函数的设计不同研究有不同偏好但基本都围绕三个指标做加权失电负荷量最小这是故障重构的第一目标能多恢复一个用户的供电比什么都重要网络损耗最小在满足供电恢复的前提下尽量让有功损耗小一点开关操作次数最少频繁倒闸操作会降低开关寿命工程上不现实。我采用的典型目标函数形式是这样的[ \min ; \sum_{i \in N} \omega_i^{LS} P_i^{LS} \alpha \sum_{(i,j) \in E} r_{ij} I_{ij}^2 \beta \sum_{(i,j) \in E} \left| z_{ij} - z_{ij}^0 \right| ]三个部分分别对应失电惩罚、网损项和开关动作惩罚。其中 (\omega_i^{LS}) 是节点i的失电权重(P_i^{LS}) 是切负荷量(z_{ij}^0) 是故障前的开关状态。这样把多目标合并成单目标工程上最好实现也方便调权重做对比分析。2.3 辐射状拓扑约束故障重构里最容易写错的一环配电网必须保持辐射状运行这是硬约束。怎么在模型里表达“一个连通图恰好是一棵树”最常见的做法是“单商品流约束”single-commodity flow constraints或者“虚拟潮流法”。思路是这样的给网络注入一个虚拟的“单位流量”根节点变电站出口是唯一源点所有其他节点都是汇点。如果网络是连通的并且节点数为N、支路数为N-1那它就是一棵树。用数学表达就是[ \sum_{j \in \Omega_i} F_{ij} - \sum_{k \in \Pi_i} F_{ki} 1, \quad \forall i \neq root ]其中 (F_{ij}) 是支路(i,j)上的虚拟流量(\Omega_i) 是节点i的子节点集合(\Pi_i) 是父节点集合。再配合[ F_{ij} \leq z_{ij} \cdot M ]当支路断开(z_{ij}0)时虚拟流量必须为0这样就保证了所有节点必须通过闭合支路从根节点获取流量从而强制网络保持连通且无环。这个约束我第一次写的时候犯了很蠢的错误——没有把虚拟流量约束和支路状态变量关联起来结果求解器返回的“最优拓扑”是一堆不连通的岛屿每个岛屿各自满足潮流方程完全失去意义。后来才知道必须用 (F_{ij} \leq z_{ij} \cdot M) 这种“大M约束”把离散状态和连续流量拧到一起。2.4 为什么不能直接求解非凸性的根源到这里模型的基本框架已经清楚了。但有一个致命问题配电网潮流方程本身是非凸的。经典的支路潮流方程DistFlow模型长这样[ P_{ij} \sum_{k \in \Omega_j} P_{jk} r_{ij} l_{ij} P_j^L - P_j^{DG} P_j^{LS} ][ Q_{ij} \sum_{k \in \Omega_j} Q_{jk} x_{ij} l_{ij} Q_j^L - Q_j^{DG} ][ v_j v_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) l_{ij} ]其中[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]这个 (l_{ij}) 的表达式就是非凸的根源。它把支路的有功、无功和节点电压耦合在一起形成一个等式约束这个等式约束对应的可行域是一个非凸的曲面。在离散开关变量的影响下整个问题就变成了混合整数非线性规划MINLP。MINLP虽然在理论上已经被研究了很多年但商业求解器的处理能力仍然有限特别是当配电网节点数增加到上百个时直接求解MINLP往往需要极长的计算时间而且无法保证在可接受时间内找到全局最优解。这就是二阶锥松弛登场的理由。3. 二阶锥松弛的精髓把非凸等式变成凸不等式3.1 变量代换先让问题看起来“顺眼”一点二阶锥松弛的第一步是做一个变量代换把潮流方程中的二次项用新变量表示。传统DistFlow方程里的 (v_i) 和 (l_{ij}) 都是原始变量电压是 (V_i) 的平方电流是 (I_{ij}) 的平方[ v_i V_i^2, \quad l_{ij} I_{ij}^2 ]这个代换本身不改变凸性但它把潮流方程的框架变得统一了所有方程中出现的都是节点的电压平方和支路的电流平方而不再是原始电压、电流的乘积。注意(v_i) 对节点i是唯一的但 (l_{ij}) 是每个支路特有的。如果支路断开对应的 (l_{ij}) 应该被强制为0这同样需要Big-M约束来处理。3.2 核心松弛把等号改成不等号代换之后DistFlow方程中的最后一个关系式[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]这是一个非凸等式因为它等价于 (l_{ij} v_i P_{ij}^2 Q_{ij}^2)左边的乘积项 (l_{ij} v_i) 导致非凸。二阶锥松弛的关键操作就是把这个等式放宽为不等式[ l_{ij} \geq \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]再改写一下得到[ l_{ij} v_i \geq P_{ij}^2 Q_{ij}^2 ]这个不等式可以用一个标准的旋转二阶锥来表示[ \left\lVert \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ v_i - l_{ij} \end{bmatrix} \right\rVert_2 \leq v_i l_{ij} ]写成更直观的矩阵形式就是[ \begin{bmatrix} v_i P_{ij} Q_{ij} \ P_{ij} l_{ij} 0 \ Q_{ij} 0 l_{ij} \end{bmatrix} \succeq 0 ]这个矩阵半正定的条件正是二阶锥规划的典型表达。YALMIP里直接写锥约束写法不同但数学上是一回事。3.3 松弛为什么是“紧”的目标函数帮你把松弛压回来到这里一定有读者会问把一个等式放宽成不等式那不是改变了物理意义吗松弛后的解会不会明明不满足真实潮流方程却被当成最优解答案是不一定但在大多数情况下最优解恰好落在锥的边界上也就是说松弛是紧的此时 (l_{ij} v_i P_{ij}^2 Q_{ij}^2) 严格成立松弛前后等价。为什么会有这种现象关键在于目标函数的结构。如果你以网损最小为优化目标也就是 (\min \sum r_{ij} l_{ij})那么由于 (l_{ij}) 是目标函数里的正项求解器有天然动力把它往小了压。而 (l_{ij}) 又不是可以无限小的——它必须满足潮流方程的下界约束 (l_{ij} \geq (P_{ij}^2 Q_{ij}^2)/v_i)——所以最优解处(l_{ij}) 会被压到锥的边界上。此时不等式等号成立松弛解就是物理上可行的真实解。这正是二阶锥松弛在配电网优化里“好用”的根本原因它把一个非凸问题变成了凸问题同时又不牺牲解的物理可实施性。当然松弛不紧的情况也是存在的特别是在目标函数不含网损项、或者网损项权重很小时。这时你会发现求解结果中的 (l_{ij}) 明显大于潮流下界说明松弛产生了间隙需要额外的手段处理比如收紧节点电压范围或者加入割平面。这部分我会在避坑章节里展开讲。3.4 从“难解”到“可解”计算复杂度上的实测对比我自己用IEEE 33节点系统做过对比。同一个故障场景直接用非线性求解器比如fmincon配合分支定界算法去解MINLP跑了一个小时还在原地打转迭代进度几乎不前进。换成MISOCP之后Gurobi在大约两三秒内就给出了全局最优解。这个差距是数量级的。当然节点数增加到IEEE 123节点甚至更大的系统时求解时间会明显上升因为0-1变量的数量变多了分支定界的搜索空间呈指数扩张。但即便如此MISOCP仍然是可以工程落地的而MINLP在同等规模下基本处于“理论可行、实际绝望”的状态。4. MATLABYALMIP建模从数学模型到可运行代码的完整路径4.1 环境准备装好YALMIP和求解器这套代码的核心依赖是YALMIP工具箱——它在MATLAB里扮演“翻译官”的角色把人类友好的数学模型自动转换成求解器能理解的内部表示。YALMIP不是求解器它本身不负责计算真正算东西的是Gurobi或Cplex。环境配置我列一下避免新手在这里耽误时间MATLAB R2019a及以上版本都可以我用的R2022bYALMIP到GitHub上下载最新版把文件夹添加到MATLAB路径即可不需要安装Gurobi需要申请学术许可证装好后在MATLAB里运行gurobi_setup完成对接备选方案如果没有Gurobi许可证可以用Cplex或者Mosek也能解MISOCP只是求解速度会有差异。我自己最常用Gurobi实测在配电网重构这种中等规模的MISOCP问题上它的分支定界性能比Cplex明显要好尤其是处理大量Big-M约束和0-1变量的时候。4.2 节点与支路数据准备IEEE 33节点系统为了让读者能直接跑通我用IEEE 33节点标准算例来演示。这个系统是国际通用的配电网测试系统包含33个节点、32条分段支路、5条联络支路基准电压12.66 kV基准功率10 MVA总负荷3715 kW 2300 kvar。数据准备阶段我会把节点负荷和支路阻抗分别整理成矩阵% 节点负荷数据: [节点编号, 有功负荷(kW), 无功负荷(kvar)] load_data [ 1, 0, 0; 2, 100, 60; 3, 90, 40; 4, 120, 80; 5, 60, 30; ... ]; % 支路数据: [首端, 末端, 电阻(ohm), 电抗(ohm), 开关类型] % 开关类型: 1表示分段开关(常闭), 0表示联络开关(常开) branch_data [ 1, 2, 0.0922, 0.0470, 1; 2, 3, 0.4930, 0.2511, 1; 3, 4, 0.3660, 0.1864, 1; ... 8, 21, 2.0000, 2.0000, 0; 9, 15, 2.0000, 2.0000, 0; ... ];这里有个容易出错的地方后续编程中节点编号最好从一开始不要混入0索引否则下标映射出错的概率很高而且调试起来非常痛苦因为报错往往不在源头而是在后续某个矩阵索引处。4.3 变量定义sdpvar还是binvar这里分清楚YALMIP里连续决策变量用sdpvar0-1决策变量用binvar。对于故障重构模型需要定义的核心变量如下% 节点数 n 33, 支路数 m 37 n 33; % 节点 m 37; % 全部支路(含联络开关) % 0-1变量: 支路开关状态 z binvar(m, 1); % 1闭合 0断开 % 连续变量: 支路有功、无功、电流平方 Pij sdpvar(m, 1); Qij sdpvar(m, 1); Iij_sq sdpvar(m, 1); % 连续变量: 节点电压平方 Vi_sq sdpvar(n, 1); % 连续变量: 节点失电负荷(故障期间的切负荷量) P_LS sdpvar(n, 1);这里有一个细节电流平方变量 (l_{ij}) 是每个支路独有的而电压平方变量 (v_i) 是每个节点独有的。在写约束时下标对应关系要前后一致支路编号和节点编号不能错位。4.4 约束构建逐条写入注意Big-M的取值接下来是约束条件。我按照五大类来组织代码第一类电压平方范围约束每个节点的电压平方必须保持在允许范围内例如0.95到1.05标幺值Vmin_sq 0.95^2; Vmax_sq 1.05^2; Constraints [Constraints, Vmin_sq Vi_sq Vmax_sq];第二类端口功率平衡约束每个节点必须满足功率平衡方程。这里要特别小心故障场景下的负荷削减量 (P_i^{LS})它是一个非负连续变量表示该节点被削减的负荷功率for i 1:n inflow_P 0; inflow_Q 0; outgoing_P 0; outgoing_Q 0; % 找到所有以i为首端或末端的支路 % (这里用关联矩阵A_in和A_out分别表示支路潮流流入/流出方向) if DG_active(i) 0 inflow_P inflow_P P_DG(i); inflow_Q inflow_Q Q_DG(i); end % 节点i的功率平衡: % 流入功率 DG出力 负荷消耗 - 削减量 流出功率 网损相关项 Constraints [Constraints, ... inflow_P sum(Pij(A_in(i,:))) ... load_P(i) - P_LS(i) sum(Pij(A_out(i,:))) sum(r .* Iij_sq(A_out(i,:)))]; end每个节点的功率平衡都得写漏掉任何一个节点最终的“最优解”都会在物理上毫无意义。我第一次跑通33节点系统时发现第18节点的功率平衡约束写错了方向导致那个节点的电压极低但目标函数值依然“最优”后来检查才知道是索引映射出了问题。第三类DistFlow电压递推约束DistFlow方程中如果支路闭合末端节点电压平方必须满足电压递推关系[ v_j v_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) l_{ij} ]如果支路断开则不需要强制该等式成立所以用Big-M松弛处理for k 1:m i branch_data(k, 1); j branch_data(k, 2); r branch_data(k, 3); x branch_data(k, 4); expression Vi_sq(j) - Vi_sq(i) 2*(r*Pij(k) x*Qij(k)) - (r^2 x^2)*Iij_sq(k); Constraints [Constraints, expression (1 - z(k)) * M]; Constraints [Constraints, expression -(1 - z(k)) * M]; endBig-M的取值需要仔细斟酌。M太小会把可行域限制得比真实情况更小甚至让原本可行的重构方案被误判为不可行M太大会导致数值病态求解器在分支定界时区分度变差。我实测过当电压标幺值在0.95~1.05范围时M取10~20比较合理。取值 (M Vmax_sq - Vmin_sq 2*(r_max x_max) * S_base) 可以作为一个保守估算。第四类二阶锥松弛约束这是整个模型的核心。在YALMIP里可以直接写成矩阵不等式也可以写成范数形式。推荐范数形式for k 1:m Constraints [Constraints, ... norm([2*Pij(k); 2*Qij(k); Vi_sq(branch_data(k,1)) - Iij_sq(k)], 2) ... Vi_sq(branch_data(k,1)) Iij_sq(k)]; end这个约束必须对每条支路加上包括联络开关支路。YALMIP会自动识别这种二次范数不等式并归类为二阶锥约束。第五类辐射状拓扑约束用虚拟流量法实现。定义虚拟流量变量 (F_{ij})并附加Big-M条件F sdpvar(m, 1); % 虚拟流量 % 节点虚拟流量平衡: 每个节点需要1单位虚拟流量 for i 2:n % 根节点除外 inflow_F sum(F(A_in(i,:))); outflow_F sum(F(A_out(i,:))); Constraints [Constraints, inflow_F - outflow_F 1]; end % 支路断开时虚拟流量为0 for k 1:m Constraints [Constraints, F(k) z(k) * M_flow]; Constraints [Constraints, F(k) -z(k) * M_flow]; end虚拟流量的M_flow可以取支路数量的级别因为每条支路的虚拟流量最多不超过节点总数。这里如果设置得太大会造成数值精度浪费影响求解效率。4.5 目标函数三项加权权重要说明白目标函数的代码实现Objective 0; % 第一项: 失电负荷惩罚 weight_LS 10; % 失电惩罚权重(单位/kW) Objective Objective weight_LS * sum(P_LS); % 第二项: 网损项 rho 1; % 网损权重 Objective Objective rho * sum(branch_data(:,3) .* Iij_sq); % 第三项: 开关操作惩罚 beta 0.1; % 开关操作权重 z_initial ones(m, 1); % 故障前开关状态(分段闭合,联络断开) Objective Objective beta * sum(abs(z - z_initial));三个权重之间的数量级关系决定了解的风格偏向哪一边。失电惩罚权重最大这是对的因为恢复供电是第一优先级。网损权重次之开关操作权重最小——因为倒一次闸的代价远小于多运行一天的高损耗。这里每个系统的具体数值需要根据你对运行要求的判断手动调。4.6 求解与结果提取模型封装好后直接调用optimize函数ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 1e-4, ... gurobi.TimeLimit, 300); result optimize(Constraints, Objective, ops); if result.problem 0 % 求解成功 z_opt value(z); Pij_opt value(Pij); Qij_opt value(Qij); Vi_opt sqrt(value(Vi_sq)); P_LS_opt value(P_LS); else disp(求解失败, 错误代码: result.info); endMIPGap设成1e-4表示相对最优间隙控制在0.01%以内对于故障重构这种工程问题这个精度足够了。追求过高的间隙比如1e-6会明显增加求解时间而实际得到的方案在物理意义上几乎没有差别。这里有一个很重要的实践体会如果求解器报infeasible不要先去怀疑求解器先检查约束的索引是否正确。我遇到过好几次“不可行”问题最终都是因为某个节点的功率平衡约束里DG出力方向写反或者支路关联矩阵的行列映射错位。顺序上先查关联矩阵再查电压递推约束最后才查锥约束和拓扑约束。5. 结果可视化把配电网“画”出来5.1 为什么要可视化以及可视化哪些内容故障重构计算结束后YALMIP会给你一堆数值每个节点的电压、每条支路的潮流、每个开关的状态。如果你只看表格很难在十秒内判断“这个方案好不好”。而可视化可以让你一眼看出重构后的网络长什么样、失电节点有哪些、电压最低点在哪里。我习惯在求解之后做四张图重构前后的网络拓扑对比图节点电压分布曲线重构前 vs 重构后各节点有功负载和DG出力柱状图叠加失电信息支路功率流向图可选但信息量很大。5.2 拓扑可视化用MATLAB的Graph对象MATLAB从R2016b开始引入了graph对象可以很方便地画网络拓扑。配合自定义节点坐标可以画出与网架结构一致的图形% 定义所有节点坐标(手动整理, 近似IEEE 33节点几何布局) pos [ 0, 0; % 节点1 1, 0; % 节点2 ... ]; figure; hold on; for k 1:m i branch_data(k, 1); j branch_data(k, 2); if z_opt(k) 1 plot([pos(i,1), pos(j,1)], [pos(i,2), pos(j,2)], b-, LineWidth, 1.5); else % 断开支路用浅色虚线 plot([pos(i,1), pos(j,1)], [pos(i,2), pos(j,2)], r--, LineWidth, 1.0); end end % 标记失电节点 idx_loss find(P_LS_opt 0.1); plot(pos(idx_loss, 1), pos(idx_loss, 2), ro, MarkerSize, 8, MarkerFaceColor, r); % 标记DG节点 idx_dg find(DG_active 0); plot(pos(idx_dg, 1), pos(idx_dg, 2), gs, MarkerSize, 8, MarkerFaceColor, g);这种画法不需要引入额外工具箱也不需要复杂的坐标计算。如果你希望更精细可以用GraphPlot对象设置NodeLabel和EdgeLabel把开关编号和状态直接标在图上。5.3 节点电压曲线重构效果最直接的证据故障发生时失电区域的节点电压会跌到很低甚至为0重构之后应当恢复到正常范围。画电压曲线时把所有节点按拓扑顺序排好画出重构前后的电压对比figure; plot(1:n, Vi_base, b-o, LineWidth, 1.2, MarkerSize, 5); hold on; plot(1:n, Vi_opt, r-s, LineWidth, 1.2, MarkerSize, 5); legend(故障时, 重构后, Location, best); xlabel(节点编号); ylabel(电压幅值 (p.u.)); grid on;从电压曲线上你能非常清楚地看到重构后有多少节点的电压被抬高了有没有节点仍然处于低电压状态。我调参数时经常只看这张图——如果某处电压还是低于0.95说明重构方案并没有完全恢复该区域的供电能力可能需要调整失电惩罚权重或增加闭合开关。5.4 收敛过程可视化判断求解器“紧不紧”的附加手段除了网络状态我还会画二阶锥松弛的“间隙图”——就是每条支路上 (l_{ij} v_i) 和 (P_{ij}^2Q_{ij}^2) 之间的差距。这个间隙是判断松弛是否紧的关键指标。gap zeros(m, 1); for k 1:m i branch_data(k, 1); lhs Vi_sq_opt(i) * Iij_sq_opt(k); rhs Pij_opt(k)^2 Qij_opt(k)^2; gap(k) abs(lhs - rhs) / max(lhs, 1e-6); % 相对间隙 end figure; bar(gap); xlabel(支路编号); ylabel(松弛间隙 (相对值));如果所有支路的相对间隙都在1e-3以下说明松弛解本身就足够接近真实潮流解结果可信。如果某些支路的间隙明显偏大说明最优解不满足物理潮流方程你需要考虑修正模型或者增加约束。5.5 配色与出图好图是给评审看的“第一张脸”MATLAB默认的配色比较“老旧”我在投稿和报告时通常会自定义一套配色保证图和论文里的其他图风格统一% 自定义配色方案 color_blue [0.00, 0.45, 0.74]; color_red [0.85, 0.33, 0.10]; color_green [0.47, 0.67, 0.19]; color_gray [0.49, 0.49, 0.49]; set(0, DefaultAxesColorOrder, [color_blue; color_red; color_green; color_gray]); set(0, DefaultLineLineWidth, 1.5); set(0, DefaultAxesFontName, Helvetica); set(0, DefaultAxesFontSize, 12); % 导出高分辨率图 exportgraphics(gcf, figure_reconfig.png, Resolution, 300);exportgraphics是R2020a之后推荐的导出函数比旧版print更稳定不会出现图片周围多余留白的问题。6. 调试经验与常见坑跑通和跑对之间隔着一整座冰山6.1 先跑小系统再上大系统我的调试顺序是先用6节点或9节点的极简配电网系统把整套代码跑通确认模型逻辑正确再切换到33节点。小系统里你可以手工验算最优解是否合理比如故障支路已知、重构路径确定那么解出来的开关状态应该和你手推的一致。这一步验证通过再放大系统才有意义。如果一上来就跑33节点一旦结果不合理你根本不知道是自己建模错了还是求解器数值问题调试会非常痛苦。6.2 Big-M不是越大越好很多新手误以为M取1e6最安全实际上这会严重恶化求解器的数值稳定性。Gurobi内部会做预求解处理但过大M仍然会带来精度问题——分支定界时上下界收敛变慢甚至误判可行解。我在33节点系统里实测M取15左右效果最好。判断依据很简单断开支路的电压递推约束被Big-M放松后最坏情况下电压差量级在M量级如果M远大于正常电压范围标幺值下大约2就属于过度放大。当然不同系统电压等级不同M的取值需要根据基准值重新估算不要盲目套用。6.3 锥约束到底写没写对用“间隙图”自检很多代码跑完求解器返回problem 0你就以为万事大吉了。但实际上如果锥约束写法有误YALMIP可能不会报错它会把一个不是凸的问题当成非凸问题处理或者生成错误的MPS文件。所以我有固定习惯求解完成后立即画出前面提到的“松弛间隙图”看看是否所有闭合支路的间隙都趋近于0。只要有一条支路间隙大于0.01就说明这个结果不能直接使用。拿这个作为质检环节比盲目信任求解器返回的“最优”靠谱得多。6.4 辐射状约束漏了会怎样把单商品流约束注释掉再跑一遍你会得到一个“网络由多个互不相连的岛屿组成”的解。从目标函数值看它可能比正确解更小因为断开更多支路减少了电流平方项使得网损项下降。看着更“优秀”的结果实际上完全不可运行。这种解法失败在物理意义上必须在模型层面避免。如果你遇到“最优解里网络不连通”的问题优先排查辐射状约束是否遗漏以及节点虚拟流量平衡等式是否写全。注意根节点的虚拟流量平衡是恒等式不需要额外约束其他每个节点都要写流入减流出等于1。6.5 DG出力约束主动配电网特有的坑在主动配电网中DG出力不是完全自由的它有最大输出功率限制[ 0 \leq P_i^{DG} \leq P_i^{DG,max} ]这个约束如果漏写DG出力会被求解器当成“无限余额”的免费电源最优解可能变成把所有负荷都塞给一台DG配电网潮流完全失真。我听过不少同学跑重构结果电压全线偏高排查到最后发现DG约束根本就没写——别笑这种事非常高频。另外DG的功率因数通常被假定为恒定即[ Q_i^{DG} P_i^{DG} \cdot \tan(\phi_i) ]如果要考虑DG无功调节能力还需要把无功也设成决策变量。这个取舍要看项目需求。如果不做无功优化建议直接固定功率因数这样问题更小、收敛更快。6.6 求解时间太长怎么办如果你是做大系统仿真节点数超过100个MISOCP的求解时间可能会从秒级涨到分钟级。实用手段有几个第一把MIPGap从1e-4放宽到1e-2求解时间通常可以缩短一个数量级而重构开关方案往往不会变化因为不同拓扑之间的目标函数差距远大于这个间隙。第二给开关操作次数加一个上限。比如要求最多只允许5次开关操作大幅压缩搜索空间。第三用初始解引导。先用启发式算法比如遗传算法粗跑一个不错的可行解传给Gurobi做初始解让分支定界树从一开始就有良好的上下界剪枝效率大幅提升。Gurobi支持通过x0参数传入初始解YALMIP里可以在sdpsettings里配置。6.7 报错信息怎么看YALMIP的报错信息里最常出现的是Converged to an infeasible point和Numerical issues。前者说明约束冲突通常排查顺序是关联矩阵索引 → 节点功率平衡方向 → 电压递推Big-M → 辐射状约束 → 锥约束。后者说明数值尺度差异太大常见的修法是把绝对单位换成标幺值。我在最初写代码时单位没有统一负荷用kW、阻抗用欧姆、电压用kV结果矩阵里数值相差几个数量级求解器直接报数值问题。后来统一成标幺值之后问题消失。这是个很低级的坑但实际项目里非常常见。7. 完整代码框架与后续扩展方向7.1 主程序框架为了方便读者对照我把主程序的代码结构整理如下你可以直接照着搭框架填充数据后运行%% 主动配电网故障重构主程序 clear; clc; close all; %% 1. 数据加载 load_branch_data; load_node_data; load_dg_data; %% 2. 参数初始化 n length(load_data); % 节点数 m size(branch_data, 1); % 支路数 S_base 10; % 基准功率 MVA V_base_sq 12.66^2; % 基准电压平方 %% 3. 构建关联矩阵 A_in zeros(n, m); % A_in(i,k)1 表示支路k的末端是节点i A_out zeros(n, m); % A_out(i,k)1 表示支路k的首端是节点i %% 4. 定义决策变量 z binvar(m, 1); Pij sdpvar(m, 1); Qij sdpvar(m, 1); Iij_sq sdpvar(m, 1); Vi_sq sdpvar(n, 1); P_LS sdpvar(n, 1); F_virtual sdpvar(m, 1); %% 5. 约束构建 Constraints []; ... %% 6. 目标函数 Objective ... %% 7. 求解 ops sdpsettings(solver, gurobi, verbose, 2); result optimize(Constraints, Objective, ops); %% 8. 后处理与可视化 post_process; plot_network; plot_voltage; plot_relaxation_gap;7.2 这套模型还能往哪些方向扩展二阶锥松弛在主动配电网故障重构上的框架可以直接迁移到几个很有价值的子问题上故障恢复后的网络重构与DG孤岛划分联合优化。在极端天气下配电网可能完全失电这时候可以把系统划分为若干孤岛每个孤岛由一台或几台DG支撑运行。这个问题的数学模型和故障重构几乎同构只是约束条件从“必须与根节点连通”变成了“每个孤岛内部必须连通”。三相不平衡配电网的故障重构。真实配电网往往是三相四线制负荷在相间存在差异。三相模型里DistFlow方程要扩展成三相版本二阶锥松弛依然适用——只是变量数量变成三倍求解难度会同步上升。考虑储能系统参与的故障恢复。储能可以看成“可移动的DG”既能放电支撑负荷也能充电吸收盈余功率。在故障重构中加入储能充放电状态变量约束条件会多出储能SOC递推方程和时间耦合约束但整体仍然能套进MISOCP的框架里。多时段滚动重构。故障恢复不是一锤子买卖负荷随时间变化DG出力也随光照、风速波动一个静态的重构方案可能只在某个时段最优。将问题扩展为多时段模型后决策变量数量会成倍增加但二阶锥松弛处理潮流方程的核心思想不变只是求解时间和内存占用会明显上升。分布式求解。如果配电网规模太大单机内存不够可以考虑把原问题分解成若干子问题用ADMM或Benders分解交替求解。每个子问题内部仍然是MISOCP。二阶锥松弛在凸问题上的分解相对成熟这也是近年来主动配电网优化研究的热点方向之一。这些扩展方向起点都是这套“DistFlow建模 二阶锥松弛 YALMIP求解 可视化验证”的基本框架。把这篇博文里的代码吃透上面任何一条路你都有能力快速上手。我从接触这个方向到现在最深的一个体会是二阶锥松弛不是你学完理论就“会了”的东西而是你要在代码里反复打磨才能找到手感的东西。同一个模型M取多大、权重怎么配、计算顺序怎么组织都会显著影响求解速度和结果质量。希望这篇博文能帮你省掉我当年绕的那些弯路让你把精力放在真正值得研究的那些问题上。 1. 为什么故障重构让很多研究者在第一步就卡壳先说说我自己的经历。早几年刚开始接触主动配电网故障重构时我在网上找了一圈资料发现绝大多数文章都在讲理论推导一上来就是配电网潮流方程、混合整数非线性规划看得懂每个字合在一起就不知道该怎么落地。真正动手写MATLAB代码时更是完全不知道从哪一行开始。后来做的时间长了回头再看发现这个方向的门槛根本不在数学公式上而在于你还没想明白三件事问题到底长什么样、为什么非用二阶锥松弛不可、以及“求解器告诉你不可行”之后该往哪个方向排查。这篇博文就是围绕这三件事展开。我会用一套完整的MATLAB实现——从配电网故障场景建模、二阶锥松弛的核心推导到YALMIP调用商业求解器求解再到结果的可视化出图——把整个过程拆开讲清楚。无论你是刚入门的硕士生还是已经在做配电网优化方向的工程师这篇文章的目标都只有一个让你照着写能跑通跑通了以后敢改敢调。先说结论层面的事二阶锥松弛Second-Order Cone RelaxationSOCR在主动配电网故障重构中的应用本质上就是解决一个“非凸找最优”的问题。配电网在正常运行时是辐射状网络但故障发生后为了恢复失电区域的供电调度人员要决定哪些联络开关闭合、哪些分段开关断开这个决策过程天然带有0-1离散变量。与此同时配电网的潮流方程本身是非线性的也就是说功率、电压、电流之间存在乘积关系这又引入了连续层面的非凸性。“离散非凸”叠加在一起就是标准的混合整数非线性规划MINLP这种问题在数学上很难得到全局最优解。业界和学术界比较统一的思路是把这个硬骨头拆成两步啃先用二阶锥松弛把非凸的潮流方程变成凸约束再用Big-M法处理离散变量和连续变量之间的耦合。这样一来原来的MINLP就变成了混合整数二阶锥规划MISOCP商业求解器比如Gurobi、Cplex、Mosek可以直接处理并且能保证在有限时间内找到全局最优解。这套思路不是某个团队独有的魔法而是过去十年主动配电网优化领域的主流范式。我个人最看重的其实是整个过程里“可视化”的价值。算完一堆状态变量之后如果没有图你很难直观地看到“哪条线路失电了”“重构后哪些节点的电压恢复了”“DG出力被削减了多少”。做工程和研究可视化不是锦上添花而是你判断结果是否合理的第一道关卡。下面我按实际上手顺序把整套方案拆成五个部分讲问题建模、数学原理、代码实现、可视化、避坑经验。每个部分我都会尽量给出能直接用的东西。2. 主动配电网故障重构的问题本质一个“开关组合潮流寻优”的双层结构2.1 主动配电网与传统配电网的差异在哪里传统配电网是无源的电能从变电站单向流向负荷调度员能控制的手段非常有限——无非是投切电容器、调节有载调压变压器。但主动配电网Active Distribution NetworkADN不一样它接入了大量分布式电源DG、储能系统和柔性负荷这些资源既可以向电网供电也可以参与调节。说白了传统配电网的故障重构是“被动恢复”——把故障隔离掉合上联络开关让电网自己重新带上负荷。而主动配电网的故障重构是“主动优化”——不仅要恢复供电还要在恢复的过程中同时考虑DG出力的调整、储能是否放电、哪些可中断负荷可以暂时切除甚至要比较“多恢复一个节点但网损变大”和“少恢复一个节点但网损更优”之间的取舍。这就是为什么题目里特别强调“主动配电网”而不是简单说“配电网”。在主动配电网的故障重构中DG出力往往被视为决策变量而不是恒定不变的负荷注入。这带来一个问题DG的接入位置、出力上限、运行方式恒功率因数或恒电压都会直接影响重构结果。2.2 故障重构的决策变量与目标函数从数学上看故障重构要回答两个问题第一哪些开关应该闭合哪些应该断开第二在确定的网络拓扑下各节点注入功率是多少。先说拓扑层的变量。配电网是由分段开关和联络开关构成的网状结构正常运行时通过断开联络开关保持辐射状。故障时某个或某几个分段开关因故障需要断开原本正常的供电路径被切断这时就要闭合一些联络开关来转供负荷。所以决策变量包括每条支路的“开/关”状态通常用0-1变量表示。再说运行层的变量。支路有功、无功、节点电压幅值平方、支路电流平方这些是连续变量。DG的有功出力在某些模型里也是连续决策变量。目标函数的设计不同研究有不同偏好但基本都围绕三个指标做加权失电负荷量最小这是故障重构的第一目标能多恢复一个用户的供电比什么都重要网络损耗最小在满足供电恢复的前提下尽量让有功损耗小一点开关操作次数最少频繁倒闸操作会降低开关寿命工程上不现实。我采用的典型目标函数形式是这样的[ \min ; \sum_{i \in N} \omega_i^{LS} P_i^{LS} \alpha \sum_{(i,j) \in E} r_{ij} I_{ij}^2 \beta \sum_{(i,j) \in E} \left| z_{ij} - z_{ij}^0 \right| ]三个部分分别对应失电惩罚、网损项和开关动作惩罚。其中 (\omega_i^{LS}) 是节点i的失电权重(P_i^{LS}) 是切负荷量(z_{ij}^0) 是故障前的开关状态。这样把多目标合并成单目标工程上最好实现也方便调权重做对比分析。2.3 辐射状拓扑约束故障重构里最容易写错的一环配电网必须保持辐射状运行这是硬约束。怎么在模型里表达“一个连通图恰好是一棵树”最常见的做法是“单商品流约束”single-commodity flow constraints或者“虚拟潮流法”。思路是这样的给网络注入一个虚拟的“单位流量”根节点变电站出口是唯一源点所有其他节点都是汇点。如果网络是连通的并且节点数为N、支路数为N-1那它就是一棵树。用数学表达就是[ \sum_{j \in \Omega_i} F_{ij} - \sum_{k \in \Pi_i} F_{ki} 1, \quad \forall i \neq root ]其中 (F_{ij}) 是支路(i,j)上的虚拟流量(\Omega_i) 是节点i的子节点集合(\Pi_i) 是父节点集合。再配合[ F_{ij} \leq z_{ij} \cdot M ]当支路断开(z_{ij}0)时虚拟流量必须为0这样就保证了所有节点必须通过闭合支路从根节点获取流量从而强制网络保持连通且无环。这个约束我第一次写的时候犯了很蠢的错误——没有把虚拟流量约束和支路状态变量关联起来结果求解器返回的“最优拓扑”是一堆不连通的岛屿每个岛屿各自满足潮流方程完全失去意义。后来才知道必须用 (F_{ij} \leq z_{ij} \cdot M) 这种“大M约束”把离散状态和连续流量拧到一起。2.4 为什么不能直接求解非凸性的根源到这里模型的基本框架已经清楚了。但有一个致命问题配电网潮流方程本身是非凸的。经典的支路潮流方程DistFlow模型长这样[ P_{ij} \sum_{k \in \Omega_j} P_{jk} r_{ij} l_{ij} P_j^L - P_j^{DG} P_j^{LS} ][ Q_{ij} \sum_{k \in \Omega_j} Q_{jk} x_{ij} l_{ij} Q_j^L - Q_j^{DG} ][ v_j v_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) l_{ij} ]其中[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]这个 (l_{ij}) 的表达式就是非凸的根源。它把支路的有功、无功和节点电压耦合在一起形成一个等式约束这个等式约束对应的可行域是一个非凸的曲面。在离散开关变量的影响下整个问题就变成了混合整数非线性规划MINLP。MINLP虽然在理论上已经被研究了很多年但商业求解器的处理能力仍然有限特别是当配电网节点数增加到上百个时直接求解MINLP往往需要极长的计算时间而且无法保证在可接受时间内找到全局最优解。这就是二阶锥松弛登场的理由。3. 二阶锥松弛的精髓把非凸等式变成凸不等式3.1 变量代换先让问题看起来“顺眼”一点二阶锥松弛的第一步是做一个变量代换把潮流方程中的二次项用新变量表示。传统DistFlow方程里的 (v_i) 和 (l_{ij}) 都是原始变量电压是 (V_i) 的平方电流是 (I_{ij}) 的平方[ v_i V_i^2, \quad l_{ij} I_{ij}^2 ]这个代换本身不改变凸性但它把潮流方程的框架变得统一了所有方程中出现的都是节点的电压平方和支路的电流平方而不再是原始电压、电流的乘积。注意(v_i) 对节点i是唯一的但 (l_{ij}) 是每个支路特有的。如果支路断开对应的 (l_{ij}) 应该被强制为0这同样需要Big-M约束来处理。3.2 核心松弛把等号改成不等号代换之后DistFlow方程中的最后一个关系式[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]这是一个非凸等式因为它等价于 (l_{ij} v_i P_{ij}^2 Q_{ij}^2)左边的乘积项 (l_{ij} v_i) 导致非凸。二阶锥松弛的关键操作就是把这个等式放宽为不等式[ l_{ij} \geq \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]再改写一下得到[ l_{ij} v_i \geq P_{ij}^2 Q_{ij}^2 ]这个不等式可以用一个标准的旋转二阶锥来表示[ \left\lVert \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ v_i - l_{ij} \end{bmatrix} \right\rVert_2 \leq v_i l_{ij} ]写成更直观的矩阵形式就是[ \begin{bmatrix} v_i P_{ij} Q_{ij} \ P_{ij} l_{ij} 0 \ Q_{ij} 0 l_{ij} \end{bmatrix} \succeq 0 ]这个矩阵半正定的条件正是二阶锥规划的典型表达。YALMIP里直接写锥约束写法不同但数学上是一回事。3.3 松弛为什么是“紧”的目标函数帮你把松弛压回来到这里一定有读者会问把一个等式放宽成不等式那不是改变了物理意义吗松弛后的解会不会明明不满足真实潮流方程却被当成最优解答案是不一定但在大多数情况下最优解恰好落在锥的边界上也就是说松弛是紧的此时 (l_{ij} v_i P_{ij}^2 Q_{ij}^2) 严格成立松弛前后等价。为什么会有这种现象关键在于目标函数的结构。如果你以网损最小为优化目标也就是 (\min \sum r_{ij} l_{ij})那么由于 (l_{ij}) 是目标函数里的正项求解器有天然动力把它往小了压。而 (l_{ij}) 又不是可以无限小的——它必须满足潮流方程的下界约束 (l_{ij} \geq (P_{ij}^2 Q_{ij}^2)/v_i)——所以最优解处(l_{ij}) 会被压到锥的边界上。此时不等式等号成立松弛解就是物理上可行的真实解。这正是二阶锥松弛在配电网优化里“好用”的根本原因它把一个非凸问题变成了凸问题同时又不牺牲解的物理可实施性。当然松弛不紧的情况也是存在的特别是在目标函数不含网损项、或者网损项权重很小时。这时你会发现求解结果中的 (l_{ij}) 明显大于潮流下界说明松弛产生了间隙需要额外的手段处理比如收紧节点电压范围或者加入割平面。这部分我会在避坑章节里展开讲。3.4 从“难解”到“可解”计算复杂度上的实测对比我自己用IEEE 33节点系统做过对比。同一个故障场景直接用非线性求解器比如fmincon配合分支定界算法去解MINLP跑了一个小时还在原地打转迭代进度几乎不前进。换成MISOCP之后Gurobi在大约两三秒内就给出了全局最优解。这个差距是数量级的。当然节点数增加到IEEE 123节点甚至更大的系统时求解时间会明显上升因为0-1变量的数量变多了分支定界的搜索空间呈指数扩张。但即便如此MISOCP仍然是可以工程落地的而MINLP在同等规模下基本处于“理论可行、实际绝望”的状态。4. MATLABYALMIP建模从数学模型到可运行代码的完整路径4.1 环境准备装好YALMIP和求解器这套代码的核心依赖是YALMIP工具箱——它在MATLAB里扮演“翻译官”的角色把人类友好的数学模型自动转换成求解器能理解的内部表示。YALMIP不是求解器它本身不负责计算真正算东西的是Gurobi或Cplex。环境配置我列一下避免新手在这里耽误时间MATLAB R2019a及以上版本都可以我用的R2022bYALMIP到GitHub上下载最新版把文件夹添加到MATLAB路径即可不需要安装Gurobi需要申请学术许可证装好后在MATLAB里运行gurobi_setup完成对接备选方案如果没有Gurobi许可证可以用Cplex或者Mosek也能解MISOCP只是求解速度会有差异。我自己最常用Gurobi实测在配电网重构这种中等规模的MISOCP问题上它的分支定界性能比Cplex明显要好尤其是处理大量Big-M约束和0-1变量的时候。4.2 节点与支路数据准备IEEE 33节点系统为了让读者能直接跑通我用IEEE 33节点标准算例来演示。这个系统是国际通用的配电网测试系统包含33个节点、32条分段支路、5条联络支路基准电压12.66 kV基准功率10 MVA总负荷3715 kW 2300 kvar。数据准备阶段我会把节点负荷和支路阻抗分别整理成矩阵% 节点负荷数据: [节点编号, 有功负荷(kW), 无功负荷(kvar)] load_data [ 1, 0, 0; 2, 100, 60; 3, 90, 40; 4, 120, 80; 5, 60, 30; ... ]; % 支路数据: [首端, 末端, 电阻(ohm), 电抗(ohm), 开关类型] % 开关类型: 1表示分段开关(常闭), 0表示联络开关(常开) branch_data [ 1, 2, 0.0922, 0.0470, 1; 2, 3, 0.4930, 0.2511, 1; 3, 4, 0.3660, 0.1864, 1; ... 8, 21, 2.0000, 2.0000, 0; 9, 15, 2.0000, 2.0000, 0; ... ];这里有个容易出错的地方后续编程中节点编号最好从一开始不要混入0索引否则下标映射出错的概率很高而且调试起来非常痛苦因为报错往往不在源头而是在后续某个矩阵索引处。4.3 变量定义sdpvar还是binvar这里分清楚YALMIP里连续决策变量用sdpvar0-1决策变量用binvar。对于故障重构模型需要定义的核心变量如下% 节点数 n 33, 支路数 m 37 n 33; % 节点 m 37; % 全部支路(含联络开关) % 0-1变量: 支路开关状态 z binvar(m, 1); % 1闭合 0断开 % 连续变量: 支路有功、无功、电流平方 Pij sdpvar(m, 1); Qij sdpvar(m, 1); Iij_sq sdpvar(m, 1); % 连续变量: 节点电压平方 Vi_sq sdpvar(n, 1); % 连续变量: 节点失电负荷(故障期间的切负荷量) P_LS sdpvar(n, 1);这里有一个细节电流平方变量 (l_{ij}) 是每个支路独有的而电压平方变量 (v_i) 是每个节点独有的。在写约束时下标对应关系要前后一致支路编号和节点编号不能错位。4.4 约束构建逐条写入注意Big-M的取值接下来是约束条件。我按照五大类来组织代码第一类电压平方范围约束每个节点的电压平方必须保持在允许范围内例如0.95到1.05标幺值Vmin_sq 0.95^2; Vmax_sq 1.05^2; Constraints [Constraints, Vmin_sq Vi_sq Vmax_sq];第二类端口功率平衡约束每个节点必须满足功率平衡方程。这里要特别小心故障场景下的负荷削减量 (P_i^{LS})它是一个非负连续变量表示该节点被削减的负荷功率for i 1:n inflow_P 0; inflow_Q 0; outgoing_P 0; outgoing_Q 0; % 找到所有以i为首端或末端的支路 % (这里用关联矩阵A_in和A_out分别表示支路潮流流入/流出方向) if DG_active(i) 0 inflow_P inflow_P P_DG(i); inflow_Q inflow_Q Q_DG(i); end % 节点i的功率平衡: % 流入功率 DG出力 负荷消耗 - 削减量 流出功率 网损相关项 Constraints [Constraints, ... inflow_P sum(Pij(A_in(i,:))) ... load_P(i) - P_LS(i) sum(Pij(A_out(i,:))) sum(r .* Iij_sq(A_out(i,:)))]; end每个节点的功率平衡都得写漏掉任何一个节点最终的“最优解”都会在物理上毫无意义。我第一次跑通33节点系统时发现第18节点的功率平衡约束写错了方向导致那个节点的电压极低但目标函数值依然“最优”后来检查才知道是索引映射出了问题。第三类DistFlow电压递推约束DistFlow方程中如果支路闭合末端节点电压平方必须满足电压递推关系[ v_j v_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) l_{ij} ]如果支路断开则不需要强制该等式成立所以用Big-M松弛处理for k 1:m i branch_data(k, 1); j branch_data(k, 2); r branch_data(k, 3); x branch_data(k, 4); expression Vi_sq(j) - Vi_sq(i) 2*(r*Pij(k) x*Qij(k)) - (r^2 x^2)*Iij_sq(k); Constraints [Constraints, expression (1 - z(k)) * M]; Constraints [Constraints, expression -(1 - z(k)) * M]; endBig-M的取值需要仔细斟酌。M太小会把可行域限制得比真实情况更小甚至让原本可行的重构方案被误判为不可行M太大会导致数值病态求解器在分支定界时区分度变差。我实测过当电压标幺值在0.95~1.05范围时M取10~20比较合理。取值 (M Vmax_sq - Vmin_sq 2*(r_max x_max) * S_base) 可以作为一个保守估算。第四类二阶锥松弛约束这是整个模型的核心。在YALMIP里可以直接写成矩阵不等式也可以写成范数形式。推荐范数形式for k 1:m Constraints [Constraints, ... norm([2*Pij(k); 2*Qij(k); Vi_sq(branch_data(k,1)) - Iij_sq(k)], 2) ... Vi_sq(branch_data(k,1)) Iij_sq(k)]; end这个约束必须对每条支路加上包括联络开关支路。YALMIP会自动识别这种二次范数不等式并归类为二阶锥约束。第五类辐射状拓扑约束用虚拟流量法实现。定义虚拟流量变量 (F_{ij})并附加Big-M条件F sdpvar(m, 1); % 虚拟流量 % 节点虚拟流量平衡: 每个节点需要1单位虚拟流量 for i 2:n % 根节点除外 inflow_F sum(F(A_in(i,:))); outflow_F sum(F(A_out(i,:))); Constraints [Constraints, inflow_F - outflow_F 1]; end % 支路断开时虚拟流量为0 for k 1:m Constraints [Constraints, F(k) z(k) * M_flow]; Constraints [Constraints, F(k) -z(k) * M_flow]; end虚拟流量的M_flow可以取支路数量的级别因为每条支路的虚拟流量最多不超过节点总数。这里如果设置得太大会造成数值精度浪费影响求解效率。4.5 目标函数三项加权权重要说明白目标函数的代码实现Objective 0; % 第一项: 失电负荷惩罚 weight_LS 10; % 失电惩罚权重(单位/kW) Objective Objective weight_LS * sum(P_LS); % 第二项: 网损项 rho 1; % 网损权重 Objective Objective rho * sum(branch_data(:,3) .* Iij_sq); % 第三项: 开关操作惩罚 beta 0.1; % 开关操作权重 z_initial ones(m, 1); % 故障前开关状态(分段闭合,联络断开) Objective Objective beta * sum(abs(z - z_initial));三个权重之间的数量级关系决定了解的风格偏向哪一边。失电惩罚权重最大这是对的因为恢复供电是第一优先级。网损权重次之开关操作权重最小——因为倒一次闸的代价远小于多运行一天的高损耗。这里每个系统的具体数值需要根据你对运行要求的判断手动调。4.6 求解与结果提取模型封装好后直接调用optimize函数ops sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 1e-4, ... gurobi.TimeLimit, 300); result optimize(Constraints, Objective, ops); if result.problem 0 % 求解成功 z_opt value(z); Pij_opt value(Pij); Qij_opt value(Qij); Vi_opt sqrt(value(Vi_sq)); P_LS_opt value(P_LS); else disp(求解失败, 错误代码: result.info); endMIPGap设成1e-4表示相对最优间隙控制在0.01%以内对于故障重构这种工程问题这个精度足够了。追求过高的间隙比如1e-6会明显增加求解时间而实际得到的方案在物理意义上几乎没有差别。这里有一个很重要的实践体会如果求解器报infeasible不要先去怀疑求解器先检查约束的索引是否正确。我遇到过好几次“不可行”问题最终都是因为某个节点的功率平衡约束里DG出力方向写反或者支路关联矩阵的行列映射错位。顺序上先查关联矩阵再查电压递推约束最后才查锥约束和拓扑约束。5. 结果可视化把配电网“画”出来5.1 为什么要可视化以及可视化哪些内容故障重构计算结束后YALMIP会给你一堆数值每个节点的电压、每条支路的潮流、每个开关的状态。如果你只看表格很难在十秒内判断“这个方案好不好”。而可视化可以让你一眼看出重构后的网络长什么样、失电节点有哪些、电压最低点在哪里。我习惯在求解之后做四张图重构前后的网络拓扑对比图节点电压分布曲线重构前 vs 重构后各节点有功负载和DG出力柱状图叠加失电信息支路功率流向图可选但信息量很大。5.2 拓扑可视化用MATLAB的Graph对象MATLAB从R2016b开始引入了graph对象可以很方便地画网络拓扑。配合自定义节点坐标可以画出与网架结构一致的图形% 定义所有节点坐标(手动整理, 近似IEEE 33节点几何布局) pos [ 0, 0; % 节点1 1, 0; % 节点2 ... ]; figure; hold on; for k 1:m i branch_data(k, 1); j branch_data(k, 2); if z_opt(k) 1 plot([pos(i,1), pos(j,1)], [pos(i,2), pos(j,2)], b-, LineWidth, 1.5); else % 断开支路用浅色虚线 plot([pos(i,1), pos(j,1)], [pos(i,2), pos(j,2)], r--, LineWidth, 1.0); end end % 标记失电节点 idx_loss find(P_LS_opt 0.1); plot(pos(idx_loss, 1), pos(idx_loss, 2), ro, MarkerSize, 8, MarkerFaceColor, r); % 标记DG节点 idx_dg find(DG_active 0); plot(pos(idx_dg, 1), pos(idx_dg, 2), gs, MarkerSize, 8, MarkerFaceColor, g);这种画法不需要引入额外工具箱也不需要复杂的坐标计算。如果你希望更精细可以用GraphPlot对象设置NodeLabel和EdgeLabel把开关编号和状态直接标在图上。5.3 节点电压曲线重构效果最直接的证据故障发生时失电区域的节点电压会跌到很低甚至为0重构之后应当恢复到正常范围。画电压曲线时把所有节点按拓扑顺序排好画出重构前后的电压对比figure; plot(1:n, Vi_base, b-o, LineWidth, 1.2, MarkerSize, 5); hold on; plot(1:n, Vi_opt, r-s, LineWidth, 1.2, MarkerSize, 5); legend(故障时, 重构后, Location, best); xlabel(节点编号); ylabel(电压幅值 (p.u.)); grid on;从电压曲线上你能非常清楚地看到重构后有多少节点的电压被抬高了有没有节点仍然处于低电压状态。我调参数时经常只看这张图——如果某处电压还是低于0.95说明重构方案并没有完全恢复该区域的供电能力可能需要调整失电惩罚权重或增加闭合开关。5.4 收敛过程可视化判断求解器“紧不紧”的附加手段除了网络状态我还会画二阶锥松弛的“间隙图”——就是每条支路上 (l_{ij} v_i) 和 (P_{ij}^2Q_{ij}^2) 之间的差距。这个间隙是判断松弛是否紧的关键指标。gap zeros(m, 1); for k 1:m i branch_data(k, 1); lhs Vi_sq_opt(i) * Iij_sq_opt(k); rhs Pij_opt(k)^2 Qij_opt(k)^2; gap(k) abs(lhs - rhs) / max(lhs, 1e-6); % 相对间隙 end figure; bar(gap); xlabel(支路编号); ylabel(松弛间隙 (相对值));如果所有支路的相对间隙都在1e-3以下说明松弛解本身就足够接近真实潮流解结果可信。如果某些支路的间隙明显偏大说明最优解不满足物理潮流方程你需要考虑修正模型或者增加约束。5.5 配色与出图好图是给评审看的“第一张脸”MATLAB默认的配色比较“老旧”我在投稿和报告时通常会自定义一套配色保证图和论文里的其他图风格统一% 自定义配色方案 color_blue [0.00, 0.45, 0.74]; color_red [0.85, 0.33, 0.10]; color_green [0.47, 0.67, 0.19]; color_gray [0.49, 0.49, 0.49]; set(0, DefaultAxesColorOrder, [color_blue; color_red; color_green; color_gray]); set(0, DefaultLineLineWidth, 1.5); set(0, DefaultAxesFontName, Helvetica); set(0, DefaultAxesFontSize, 12); % 导出高分辨率图 exportgraphics(gcf, figure_reconfig.png, Resolution, 300);exportgraphics是R2020a之后推荐的导出函数比旧版print更稳定不会出现图片周围多余留白的问题。6. 调试经验与常见坑跑通和跑对之间隔着一整座冰山6.1 先跑小系统再上大系统我的调试顺序是先用6节点或9节点的极简配电网系统把整套代码跑通确认模型逻辑正确再切换到33节点。小系统里你可以手工验算最优解是否合理比如故障支路已知、重构路径确定那么解出来的开关状态应该和你手推的一致。这一步验证通过再放大系统才有意义。如果一上来就跑33节点一旦结果不合理你根本不知道是自己建模错了还是求解器数值问题调试会非常痛苦。6.2 Big-M不是越大越好很多新手误以为M取1e6最安全实际上这会严重恶化求解器的数值稳定性。Gurobi内部会做预求解处理但过大M仍然会带来精度问题——分支定界时上下界收敛变慢甚至误判可行解。我在33节点系统里实测M取15左右效果最好。判断依据很简单断开支路的电压递推约束被Big-M放松后最坏情况下电压差量级在M量级如果M远大于正常电压范围标幺值下大约2就属于过度放大。当然不同系统电压等级不同M的取值需要根据基准值重新估算不要盲目套用。6.3 锥约束到底写没写对用“间隙图”自检很多代码跑完求解器返回problem 0你就以为万事大吉了。但实际上如果锥约束写法有误YALMIP可能不会报错它会把一个不是凸的问题当成非凸问题处理或者生成错误的MPS文件。所以我有固定习惯求解完成后立即画出前面提到的“松弛间隙图”看看是否所有闭合支路的间隙都趋近于0。只要有一条支路间隙大于0.01就说明这个结果不能直接使用。拿这个作为质检环节比盲目信任求解器返回的“最优”靠谱得多。6.4 辐射状约束漏了会怎样把单商品流约束注释掉再跑一遍你会得到一个“网络由多个互不相连的岛屿组成”的解。从目标函数值看它可能比正确解更小因为断开更多支路减少了电流平方项使得网损项下降。看着更“优秀”的结果实际上完全不可运行。这种解法失败在物理意义上必须在模型层面避免。如果你遇到“最优解里网络不连通”的问题优先排查辐射状约束是否遗漏以及节点虚拟流量平衡等式是否写全。注意根节点的虚拟流量平衡是恒等式不需要额外约束其他每个节点都要写流入减流出等于1。6.5 DG出力约束主动配电网特有的坑在主动配电网中DG出力不是完全自由的它有最大输出功率限制[ 0 \leq P_i^{DG} \leq P_i^{DG,max} ]这个约束如果漏写DG出力会被求解器当成“无限余额”的免费电源最优解可能变成把所有负荷都塞给一台DG配电网潮流完全失真。我听过不少同学跑重构结果电压全线偏高排查到最后发现DG约束根本就没写——别笑这种事非常高频。另外DG的功率因数通常被假定为恒定即[ Q_i^{DG} P_i^{DG} \cdot \tan(\phi_i) ]如果要考虑DG无功调节能力还需要把无功也设成决策变量。这个取舍要看项目需求。如果不做无功优化建议直接固定功率因数这样问题更小、收敛更快。6.6 求解时间太长怎么办如果你是做大系统仿真节点数超过100个MISOCP的求解时间可能会从秒级涨到分钟级。实用手段有几个第一把MIPGap从1e-4放宽到1e-2求解时间通常可以缩短一个数量级而重构开关方案往往不会变化因为不同拓扑之间的目标函数差距远大于这个间隙。第二给开关操作次数加一个上限。比如要求最多只允许5次开关操作大幅压缩搜索空间。第三用初始解引导。先用启发式算法比如遗传算法粗跑一个不错的可行解传给Gurobi做初始解让分支定界树从一开始就有良好的上下界剪枝效率大幅提升。Gurobi支持通过x0参数传入初始解YALMIP里可以在sdpsettings里配置。6.7 报错信息怎么看YALMIP的报错信息里最常出现的是Converged to an infeasible point和Numerical issues。前者说明约束冲突通常排查顺序是关联矩阵索引 → 节点功率平衡方向 → 电压递推Big-M → 辐射状约束 → 锥约束。后者说明数值尺度差异太大常见的修法是把绝对单位换成标幺值。我在最初写代码时单位没有统一负荷用kW、阻抗用欧姆、电压用kV结果矩阵里数值相差几个数量级求解器直接报数值问题。后来统一成标幺值之后问题消失。这是个很低级的坑但实际项目里非常常见。7. 完整代码框架与后续扩展方向7.1 主程序框架为了方便读者对照我把主程序的代码结构整理如下你可以直接照着搭框架填充数据后运行%% 主动配电网故障重构主程序 clear; clc; close all; %% 1. 数据加载 load_branch_data; load_node_data; load_dg_data; %% 2. 参数初始化 n length(load_data); % 节点数 m size(branch_data, 1); % 支路数 S_base 10; % 基准功率 MVA V_base_sq 12.66^2; % 基准电压平方 %% 3. 构建关联矩阵 A_in zeros(n, m); % A_in(i,k)1 表示支路k的末端是节点i A_out zeros(n, m); % A_out(i,k)1 表示支路k的首端是节点i %% 4. 定义决策变量 z binvar(m, 1); Pij sdpvar(m, 1); Qij sdpvar(m, 1); Iij_sq sdpvar(m, 1); Vi_sq sdpvar(n, 1); P_LS sdpvar(n, 1); F_virtual sdpvar(m, 1); %% 5. 约束构建 Constraints []; ... %% 6. 目标函数 Objective ... %% 7. 求解 ops sdpsettings(solver, gurobi, verbose, 2); result optimize(Constraints, Objective, ops); %% 8. 后处理与可视化 post_process; plot_network; plot_voltage; plot_relaxation_gap;7.2 这套模型还能往哪些方向扩展二阶锥松弛在主动配电网故障重构上的框架可以直接迁移到几个很有价值的子问题上故障恢复后的网络重构与DG孤岛划分联合优化。在极端天气下配电网可能完全失电这时候可以把系统划分为若干孤岛每个孤岛由一台或几台DG支撑运行。这个问题的数学模型和故障重构几乎同构只是约束条件从“必须与根节点连通”变成了“每个孤岛内部必须连通”。三相不平衡配电网的故障重构。真实配电网往往是三相四线制负荷在相间存在差异。三相模型里DistFlow方程要扩展成三相版本二阶锥松弛依然适用——只是变量数量变成三倍求解难度会同步上升。考虑储能系统参与的故障恢复。储能可以看成“可移动的DG”既能放电支撑负荷也能充电吸收盈余功率。在故障重构中加入储能充放电状态变量约束条件会多出储能SOC递推方程和时间耦合约束但整体仍然能套进MISOCP的框架里。多时段滚动重构。故障恢复不是一锤子买卖负荷随时间变化DG出力也随光照、风速波动一个静态的重构方案可能只在某个时段最优。将问题扩展为多时段模型后决策变量数量会成倍增加但二阶锥松弛处理潮流方程的核心思想不变只是求解时间和内存占用会明显上升。分布式求解。如果配电网规模太大单机内存不够可以考虑把原问题分解成若干子问题用ADMM或Benders分解交替求解。每个子问题内部仍然是MISOCP。二阶锥松弛在凸问题上的分解相对成熟这也是近年来主动配电网优化研究的热点方向之一。这些扩展方向起点都是这套“DistFlow建模 二阶锥松弛 YALMIP求解 可视化验证”的基本框架。把这篇博文里的代码吃透上面任何一条路你都有能力快速上手。我从接触这个方向到现在最深的一个体会是二阶锥松弛不是你学完理论就“会了”的东西而是你要在代码里反复打磨才能找到手感的东西。同一个模型M取多大、权重怎么配、计算顺序怎么组织都会显著影响求解速度和结果质量。希望这篇博文能帮你省掉我当年绕的那些弯路让你把精力放在真正值得研究的那些问题上。 
