拿到这个课题时我第一反应是配电网故障重构听上去像个调度问题但真正动手用 Matlab 做起来核心反而落在数学优化上。你手里有一张已经失电的下游网络故障支路被保护切开剩下那些分段开关和联络开关怎么组合才能既把所有能救的负荷救回来又不违反电压、电流和拓扑约束。这本身是个组合爆炸问题而二阶锥SOCP是近几年把这类问题从“启发式碰运气”变成“可证明最优”的最实用工具之一。这篇东西不是教科书复述是我自己用 Matlab YALMIP 搭这套模型时踩过的坑和验证过的路适合刚接触配电网重构、或者会用 Matlab 但没碰过凸优化的同学参考。1. 配电网故障重构到底在解决什么问题1.1 从一张失电负荷表说起配电网和输电网最大的区别在于一个“配”字网络是闭环设计、开环运行的。正常情况下所有分段开关闭合、联络开关断开系统呈辐射状一旦某条支路发生永久性故障断路器或馈线自动化动作把故障段隔离故障点下游的一大片负荷就全部失电。故障重构要做的就是在这个“局部网络被切开”的状态下重新组合剩余可操作的开关——闭合部分联络开关、断开部分分段开关——把失电负荷转移到其他正常馈线上。这里的关键点在于不是把所有能合的开关全合上就行因为配电网不允许闭环运行合上所有开关会形成环网保护整定、短路电流、潮流分布全部乱套。所以故障重构本质上是一个带约束的组合优化问题。变量是每个开关的开/合状态约束包括辐射状拓扑、节点电压上下限、支路容量目标是尽可能恢复失电负荷同时兼顾网损、开关操作次数、重要负荷优先恢复等工程指标。1.2 为什么不能“暴力枚举所有开关组合”有人会问一个配电网的分段开关和联络开关加起来也就几十个暴力枚举不行吗行但只限于很小的系统。IEEE 33 节点配电网算例有 5 个联络开关和 32 个分段开关全部开关组合是 2 的 37 次方约 1370 亿种组合。就算你根据拓扑可行性筛掉大部分剩下的可行解仍然数量庞大而实际配电网动辄上百个节点故障后还要考虑不同故障位置、不同负荷水平枚举法根本撑不住。这也是为什么故障重构领域从早期的启发式算法支路交换法、最优流模式法发展到遗传算法、粒子群等智能算法再到现在主流的混合整数凸优化方法。前两类方法的问题是要么依赖初始拓扑和局部寻优容易陷入局部最优要么收敛性没有严格保证跑十次可能得到十个不同结果。而基于二阶锥松弛的混合整数凸规划把非线性潮流的难点用数学工具化解掉求解器给出的是带最优性边界的全局解这在工程上意义完全不一样。2. 为什么偏偏选二阶锥SOCP2.1 DistFlow 潮流方程的非凸性来源配电网潮流计算有个经典简化模型叫 DistFlow它利用辐射状网络的递推特性用支路有功 P、无功 Q、末端电压幅值平方 U 来描述潮流U_j U_i - 2*(r*P x*Q) (r^2 x^2) * (P^2 Q^2) / U_i这个式子里(P^2 Q^2) / U_i 那一项就是网损项。问题出在哪P、Q、U 都是变量P^2、Q^2 是非线性二次项再除以 U_i 就变成非凸项。如果直接把这一项放进优化模型整个可行域就是非凸的求解器只能靠启发式或者局部搜索没法保证找到全局最优。对辐射状配电网来说这个非凸性主要就是“二次项除变量”造成的。如果能把它转成凸约束剩下的二进制开关变量虽然还是混合整数问题但连续部分凸了整体就可以交给成熟的分支定界框架去求解——这正是现代商业求解器最擅长的领域。2.2 旋转锥松弛的数学思路与精确性条件二阶锥松弛的思路非常直接引入一个辅助变量 s令s_k (P_k^2 Q_k^2) / U_i然后把上面的等式潮流U_j U_i - 2*(r*P x*Q) (r^2 x^2) * s中的 s 替进去。这个不等式在数学上等价于一个旋转二阶锥约束(2P)^2 (2Q)^2 (U_i - s)^2 (U_i s)^2旋转锥是凸的于是非凸的 DistFlow 等式被“松弛”成了凸约束。为什么叫松弛因为原本的严格等式要求 s 恰好等于 (P^2Q^2)/U_i现在只要求 s 大于等于它。如果最优解里这个不等式是紧的即等号成立那么这个松弛就是精确的求出来的结果就是原问题的全局最优解。什么时候松弛是精确的理论上已经有比较成熟的研究在辐射状网络、负荷有界且没有逆向潮流等条件下DistFlow 的二阶锥松弛通常是紧的。实际工程中我做 IEEE 33、IEEE 69 节点算例时基本每次都能得到紧解松弛间隙在 1e-6 量级。但注意这并不意味着可以完全不检查——后面我会专门说怎么验证松弛是否精确。2.3 和遗传算法、MILP 线性化方案横向对比我在做方案选型的时候其实还对比过另外两条路线这里直接说结论。第一条路线是智能算法遗传算法/粒子群是配电网重构论文里的老面孔了。优点是建模门槛低不用理解凸优化的数学细节把开关状态编码成染色体、目标函数写成适应度就行。缺点是第一每次搜索都要反复调用潮流计算33 节点还好几百个节点就很慢第二收敛结果不稳定同样的参数跑几次结果可能不一样这对工程验收和论文复现都很致命第三很难处理约束——辐射状约束和电压约束通常只能用罚函数罚系数怎么定是个玄学。第二条路线是 MILP 线性化。把 DistFlow 里的二次项通过分段线性化或大 M 法等手段变成线性约束好处是模型简单、求解快坏处是精度受线性化分段数影响分段多了变量爆炸分段少了误差大。尤其网损项对电压幅值敏感时线性化误差直接传导到重构结果里。SOCP 路线在这三者中属于“精度和效率的平衡点”二次项是精确建模的只是把等式松弛成锥约束不需要做近似求解速度上Cplex/Gurobi 对二阶锥混合整数规划的支撑已经很成熟33 节点算例通常几秒到几十秒就能收敛到全局最优。对比维度遗传算法/粒子群MILP 线性化二阶锥松弛全局最优性无保证有保证线性化误差内有保证松弛精确时求解速度慢反复潮流计算快快建模精度依赖潮流计算精度依赖线性化分段数二次项精确表达代码复杂度低中中高适合场景教学演示、小系统中大规模快速估算工程级精确分析3. Matlab YALMIP 的实现框架3.1 准备工作工具箱与算例数据Matlab 里做二阶锥混合整数规划我强烈建议不要手写求解算法直接上 YALMIP 商业求解器的组合。YALMIP 是一个建模语言它的作用是把你的数学约束“翻译”成求解器能识别的标准形式你负责写模型它负责对接求解器。求解器我优先推荐 Cplex 或 Gurobi 这两个商业软件学术 license 免费申请支持二阶锥约束和混合整数规划如果申请不到也可以用开源的 SCIP 或 ECOS 先跑通模型但大规模时性能差距比较明显。安装流程不展开只说两个坑第一YALMIP 要加到 Matlab 路径里pathtool里添加文件夹后记得savepath第二求解器要确保 Matlab 能调用到可执行文件yalmiptest命令能直接验证工具箱是否配置成功。配电网算例数据我推荐从 IEEE 33 节点开始。标准参数是基准电压 12.66 kV基准容量 10 MVA总负荷约 3715 kW 2300 kvar5 个联络开关分布在 8-21、9-15、12-22、18-33、25-29不同版本编号略有差异。数据组织建议用一个结构体存节点和支路信息支路数组每行是[首端节点, 末端节点, 电阻(标幺), 电抗(标幺), 初始开关状态]。3.2 模型变量与约束的代码骨架下面是核心代码框架我整理成可以直接改着用的结构。先定义变量% 节点数 nb支路数 nl z binvar(nl, 1); % 支路开关状态1闭合0断开 U sdpvar(nb, 1); % 各节点电压幅值平方 P sdpvar(nl, 1); % 支路首端有功 Q sdpvar(nl, 1); % 支路首端无功 s sdpvar(nl, 1); % 辅助变量表示 (P^2Q^2)/U_i然后是 DistFlow 约束。这里有个关键细节等式潮流只在开关闭合时成立开关断开时该支路潮流必须为 0节点电压也不需要满足递推关系。我在实现中用 Big-M 法把等式拆成两个不等式乘上开关状态M 10; % Big-M 取值见 3.3 节说明 Constraints []; for k 1:nl i branch(k, 1); j branch(k, 2); r branch(k, 3); x branch(k, 4); % 旋转锥约束s (P^2 Q^2) / U_i % 用标准二阶锥形式表示旋转锥 Constraints [Constraints, ... cone([2*P(k); 2*Q(k); U(i)-s(k)], U(i)s(k))]; % DistFlow 潮流等式用 Big-M 解耦 Constraints [Constraints, ... U(j) - U(i) 2*(r*P(k) x*Q(k)) - (r^2x^2)*s(k) M*(1-z(k))]; Constraints [Constraints, ... U(j) - U(i) 2*(r*P(k) x*Q(k)) - (r^2x^2)*s(k) -M*(1-z(k))]; % 开关断开时支路潮流清 0 Constraints [Constraints, ... -M*z(k) P(k) M*z(k)]; Constraints [Constraints, ... -M*z(k) Q(k) M*z(k)]; end节点功率平衡也要分情况。对每个节点注入功率等于该节点所连支路潮流之和for j 1:nb % 节点 j 的净注入功率DG 出力 - 负荷 P_inj P_dg(j) - P_load(j); Q_inj Q_dg(j) - Q_load(j); % 与该节点相连的支路 from_idx find(branch(:,1) j); to_idx find(branch(:,2) j); Constraints [Constraints, ... sum(P(from_idx)) - sum(P(to_idx)) P_inj 0]; Constraints [Constraints, ... sum(Q(from_idx)) - sum(Q(to_idx)) Q_inj 0]; end这里要注意方向定义我约定支路潮流方向是从首端流向末端所以对节点 j 来说以 j 为首端的支路是流出以 j 为末端的支路是流入。负荷用“负注入”表示DG 用“正注入”。电压约束直接加边界Umin 0.95^2; % 电压下限 0.95 pu Umax 1.05^2; % 电压上限 1.05 pu Constraints [Constraints, Umin U Umax];辐射状约束是最容易写错的部分。完整的辐射状约束 闭合支路数 节点数 - 1 网络连通。只加前者不够因为可能形成“孤岛环网”的组合只加连通性也没用可能多条支路成环。我采用单商品流约束来同时保证连通性% 闭合支路数约束 Constraints [Constraints, sum(z) nb - 1]; % 单商品流以节点 1 为根节点 f sdpvar(nl, 1); % 虚拟流量 for k 1:nl i branch(k, 1); j branch(k, 2); Constraints [Constraints, 0 f(k) nb * z(k)]; end % 根节点净流出为 -(nb-1)其他节点净流出为 1 for j 1:nb from_idx find(branch(:,1) j); to_idx find(branch(:,2) j); if j 1 Constraints [Constraints, ... sum(f(from_idx)) - sum(f(to_idx)) -(nb-1)]; else Constraints [Constraints, ... sum(f(from_idx)) - sum(f(to_idx)) 1]; end end单商品流的含义很直观每个非根节点向根节点发送 1 单位虚拟流量根节点总共接收 nb-1 单位。如果某个节点不连通流量约束必然矛盾如果存在环网闭合支路数约束会强制多断开一条支路。两者配合得到的拓扑一定是辐射状的。3.3 Big-M 参数与边界设置的经验Big-M 的取值是我调试时翻车最多的地方。太大会导致数值病态Cplex 容易报“numerical trouble”或者收敛到错误解太小则会把可行域错误截断明明有解却报 infeasible。我的经验是M 的取值要和潮流方程的量纲匹配。把功率、电压都换成标幺值后DistFlow 方程里 rP 的量级大约是 0.01~0.1xQ 同理(r^2x^2)*s 的量级更小所以 M 取 1~10 就足够了。不要因为“担心不够大”就随手填 1000那会让约束矩阵的条件数急剧恶化。另外节点电压平方 U 的量级在 0.9~1.1 之间所以 Umin 0.95^2 这种写法不要写成 Umin 0.95。很多新手在这里踩坑电压下限明明是 0.95 pu平方后是 0.9025写错的话模型直接无解或者解出来的电压偏低。故障场景的设置也不难比如支路 12-13 发生永久故障直接强制z(12) 0具体支路编号取决于你的数据结构并把这个支路从可重构开关集合里剔除即可。如果你想模拟“故障隔离后下游失电”的场景可以把故障点下游节点负荷标记为失电负荷在目标函数里加上恢复奖励项。3.4 目标函数怎么设计才贴近工程目标函数是重构问题的核心我建议至少包含两部分网损项和恢复项。网损项的表达式是sum(r .* (P.^2 Q.^2) ./ U)但注意这里用的是辅助变量 s所以可以直接写成objective sum(r .* s); % 网损之和恢复项用于故障场景。对失电负荷节点失电惩罚权重设高模型会优先闭合路径把该节点恢复供电。权重怎么设我一般把重要负荷如医院、数据中心权重设为普通负荷的 10~20 倍这样即使不能恢复全部负荷也会优先恢复关键负荷更贴近配电调度的实际需求。% 恢复惩罚未恢复的失电负荷权重累加 % w_load 为各节点负荷权重nodal_load 为有功负荷 recovery_penalty sum(w_load .* nodal_load .* (1 - restored_flag)); objective objective 1000 * recovery_penalty;restored_flag 怎么定义一个节点只要有任意一条连通的潮流路径从电源点供电它就算恢复。在 DistFlow 框架里判断节点是否恢复的一个常用近似是检查该节点电压是否在正常范围内。更严格的做法是引入恢复二进制变量并和拓扑连通绑定——但这样会显著增加模型复杂度。工程上我通常用后处理验证的方式求解完成后对失电节点逐个检查拓扑连通路径确认恢复状态。4. 仿真结果怎么解读4.1 以 IEEE 33 节点为例看故障重构前后对比我用标准 IEEE 33 节点算例做过一组完整测试故障场景设置为支路 12-13 故障下游约 4 个节点失电。重构前的状态是故障支路断开下游负荷全部失电联络开关全部保持断开变压器到网络末端电压明显偏低。重构后的结果非常直观系统闭合了联络开关 12-22具体编号以数据文件为准断开分段开关 11-12失电负荷全部恢复供电网络最低电压从 0.928 pu 恢复到了 0.95 pu 以上系统网损从故障后的某个高值降回到接近正常运行水平。这里要提醒一点重构结果不是唯一的。同一个故障场景如果目标函数里网损权重和恢复权重比例不同最优开关组合可能不一样。比如网损权重更高时模型会选择让网络更接近“均衡分配”的拓扑恢复权重更高时模型会优先把失电负荷接回来哪怕网损稍高。典型结果对比可以参考指标故障隔离后重构前故障重构后失电负荷约 300 kW0最低电压0.928 pu0.952 pu系统网损偏高网络潮流不均降低 20%~30%开关操作次数02~3 次4.2 收敛性、松弛间隙与求解时间SOCP 模型求解完成后一定要检查两件事求解器报告的 gap以及二阶锥约束的“紧度”。YALMIP 求解结束后用optimize返回的诊断信息可以直接看求解状态。gap 一般要求控制在 0.1% 以内33 节点算例通常能轻松达到 0.01% 以下。如果 gap 降不下去先检查是不是 M 值过大导致数值问题。二阶锥紧度的检查方法是读回value(s)和value(U)、value(P)、value(Q)计算s - (P.^2 Q.^2)./U看是否接近 0。我在运行中得到的结果这个差值基本都在 1e-6 量级说明松弛是精确的。如果发现差值很大说明锥约束被松弛掉了需要排查负荷设置、电压边界等条件。求解时间方面IEEE 33 节点用 Cplex 跑完整混合整数二阶锥模型包含 37 个二进制变量和 100 多个连续变量在我自己的笔记本上大约 3~10 秒收敛到全局最优。IEEE 69 节点大约 15~30 秒。这个速度对离线重构分析完全够用但如果要做实时故障恢复秒级响应还需要进一步降维比如只把故障影响范围内的开关纳入候选集。5. 常见问题与排查技巧实录5.1 求解器报 Infeasible模型无解这是最常遇到的问题。第一次遇到别慌按这个顺序排查先检查电压上下限是否写成了标幺值而不是标幺值平方再检查单商品流约束里根节点编号是否和数据一致最后检查是不是负荷数据方向搞反了——把负荷当成注入、DG 当成流出潮流直接就不可行。我自己的排错技巧是先把所有二进制变量固定为故障前的正常状态跑一个纯二阶锥可行性问题。如果能求解说明连续约束没问题问题出在开关组合如果连固定拓扑都无解那必然是模型写错了跟重构无关。5.2 解出来的拓扑带环网或孤岛拓扑是辐射状但代码跑出来却有环十有八九是辐射状约束没写全。再次强调闭合支路数等于节点数减一这只是必要条件。如果没有单商品流约束模型完全可以把两个环和一个孤岛组合在一起支路数照样满足。加上单商品流之后只要 f 的上下界和 z 绑定正确环网和孤岛都会被排除。另一个隐蔽问题是单商品流约束里 f 的整数性。f 不需要是整数变量连续变量就够——但要注意f 的流量单位必须和支路数同一量纲。f(k) nb * z(k)这个约束里nb 是节点数f 的量级必须是“每节点 1 单位流量”不能把 f 设置为兆瓦级功率。5.3 数值病态问题Big-M 太大模型“假不可行”Big-M 取值过大会导致 Cplex/Gurobi 在分支定界的过程中出现大量数值不可行现象是log 里出现 numerical trouble求解器反复尝试预求解但 gap 一直降不下来甚至报 infeasible但你把 M 改小后立刻就能求解。我的建议是把 M 的最小值设为潮流方程中所有项系数量级总和的 10 倍即可。用标幺值建模时M 10 是安全的如果系统特别大最多取到 50~100不要再高。另一个配套操作是给 Cplex 设置NumericalEmphasis参数为 1能显著提升数值稳定性。5.4 求解时间过长怎么办33 节点跑几十秒属于正常但如果你算的是 200 节点以上的实际馈线完整模型可能要好几分钟甚至更久。三个降维技巧第一缩小开关候选集。故障重构只需要操作“和失电区域有电气路径关联”的开关那些离故障点很远的联络开关基本不需要动作。人为把候选集缩小到故障区域周边 2~3 层二进制变量数量能降一半以上。第二设置求解时间上限和 gap 容忍。对于工程分析来说gap 到 1% 就可以接受设置Cplex的timelimit和mipgaptol参数让模型在可接受时间内给出近似最优解而不是死磕全局最优。第三用两阶段法先用一个简化 MILP 快速得到一个较好的初始解再把这个解作为 warm start 喂给 SOCP 模型能明显加速收敛。5.5 松弛不精确时怎么处理虽然理论条件和实际算例都说明旋转锥松弛通常是紧的但总有例外——特别是在计及 DG 逆变器无功调节、负荷模型比较独特的场景下。松弛不精确的表现是最优解里 s 明显大于 (P^2Q^2)/U_i对应支路的网损被高估重构结果可能偏离真实最优。处理办法有两个。第一加约束强制锥紧对某些关键支路直接限制s (1epsilon) * (P^2Q^2)/U_i但这种约束本身是非凸的实际效果有限。第二后验修正SOCP 求解后固定开关组合用标准潮流计算器重新计算精确网损和电压验证结果可行性。如果差异在可接受范围直接用潮流结果作为最终方案。我个人经验是配电网重构场景下松弛不精确属于偶发情况不会影响整体方案但一定要有检查意识。每次跑完都看一眼 gaptightness形成一个固定习惯比事后发现模型不可信要踏实得多。5.6 常见问题速查表异常现象可能原因处理办法报 Infeasible电压边界写成 pu 而非 pu 平方检查 Umin 0.95^2报 Infeasible负荷/ DG 功率方向写反检查节点注入功率表达式拓扑有环只加了支路数N-1缺连通性加单商品流约束M 太大导致数值病态Big-M 设置不合理M 取 10~50加数值强调求解时间过长候选开关集太大缩小候选集、设 gap 容忍松弛不精确锥约束不紧检查紧度后验潮流验证写在最后这套基于 Matlab 与二阶锥的配电网故障重构模型我是从“看着论文公式发懵”到“能独立跑出结果”一路摸索过来的。说实话最难的不是把代码写出来而是理解每个约束为什么这么写、每个参数为什么这么设。尤其是辐射状约束和 Big-M 这两个点我反复调试了至少一个礼拜。如果你也正好卡在这两个地方耐心调把模型拆开逐个检查一定能跑通。