搞配电网故障恢复的人大概都有同感重构和孤岛划分这俩问题分开弄各有各的成熟办法一旦要放一起同时优化就各种打架。前段时间我把这两块重新整理塞进同一个优化模型用Matlab重新实现了一个升级版本。这次重构把动态孤岛、多时段耦合、负荷分级都加了进去实测下来比之前的基础版本稳定不少。这篇文章把我搭建这个统一模型的过程、Matlab实现要点、以及调试时踩过的坑完整写出来。适合正在做主动配电网故障恢复方向的研究生、做配网调度算法开发的工程师也适合刚接触配电网优化、想找个完整范例入门的同学。看完你能知道为什么重构和孤岛划分必须统一建模也能直接照着代码结构复现一版能跑的模型。1. 重构与孤岛划分统一建模的必要性1.1 把两个决策放进一个模型的原因配电网故障恢复简单说就是线路故障后通过开关操作把非故障失电区域的负荷重新送上电。传统配网只有分段开关和联络开关恢复手段就是倒闸操作也就是网络重构。但现在分布式电源接入比例越来越高故障后可以让一部分区域主动脱离主网靠本地分布式电源继续供电这就是孤岛运行。问题是这两个决策并不是独立的。你选择哪些联络开关闭合、哪些分段开关断开本质上就是在重新画配电网的拓扑边界。而孤岛的边界同样是由一组断开的支路来标记的。换句话说重构和孤岛划分都在做同一件事确定配电网最终运行拓扑中的开关状态组合。如果分开建模先做孤岛划分再做重构很可能出现孤岛划分方案在重构约束下根本不可行的情况反过来先重构再划分孤岛又会忽略一些孤岛运行的潜力。只有把开关状态和孤岛归属放到一个模型里同时决策得到的结果才是全局最优的恢复方案。1.2 升级版本相比基础版改进的地方基础版本一般只做一个时间断面的静态恢复给定故障位置、给定负荷大小、给定分布式电源出力求一组开关状态使得恢复负荷量最大。这种模型用来理解原理没问题但离工程实践有距离。实际调度中故障恢复持续几个小时这段时间里负荷在变、光伏出力在变、储能电量在变一个静态断面的最优方案放到下一个时段可能就不再满足约束。这个升级版本主要做了三件事把单时段模型扩展成多时段滚动恢复模型时段间考虑开关动作次数限制和储能荷电状态转移。引入负荷分级一、二级负荷优先恢复普通负荷在功率不足时允许部分切除。加入了孤岛内功率平衡的显式约束保证每个孤岛在恢复期间都能自维持而不是只在某一个时刻能撑住。这三项改进带来的直接效果是恢复方案不再只是“纸面上能算出来”而是真正能在调度系统里执行下去的方案。1.3 这个模型的适用范围统一模型适合的场景很明确配电网发生永久性故障需要进行故障隔离和恢复重构系统里有一定比例的分布式电源和储能具备孤岛运行条件调度部门希望一次性给出从故障发生到恢复完成的完整开关操作序列。对于只有单一电源的传统辐射状配网孤岛划分意义不大单独做重构就够了没必要用这个模型。2. 数学模型目标函数与约束条件怎么搭2.1 目标函数怎么设计才合理统一模型的目标函数我建议不要只写一个“恢复负荷最大”因为单独这一项会引发一个问题模型为了多恢复一点点负荷可能让开关来回倒换十几次。实际配网开关是有动作次数寿命限制的操作越频繁故障隐患越大。常用的做法是双目标加权$$ \max \sum_{t}\sum_{i} w_i \cdot P_{i,t}^{load} \cdot x_{i,t} - \lambda \sum_{(i,j)} |s_{ij,t} - s_{ij,t-1}| $$第一项是恢复负荷收益$w_i$是节点负荷权重$P_{i,t}^{load}$是节点有功负荷$x_{i,t}$是节点在时段t是否带电。第二项是开关动作惩罚$s_{ij,t}$表示支路开关状态整个求和项统计所有开关的累计动作次数。权重$w_i$怎么定最简单的是按负荷等级给一、二、三级负荷分别取100、10、1。注意权重差距不要太大太大相当于强制要求优先恢复一级负荷可能引发无法恢复三级负荷时把一级负荷以外的全部切掉差距太小又体现不出重要性差别。工程上一般取10倍到100倍之间的比例比较合适。惩罚系数$\lambda$也要小心。我见过有人把$\lambda$设得很大结果模型为了少动一次开关宁可少恢复几千千瓦负荷。正确做法是先跑一版不带惩罚项的模型拿到最优恢复量作为基准然后设置一个远小于负荷权重量级的$\lambda$让开关动作惩罚只在恢复量相同或相近时起到“打破平局”的作用。2.2 潮流约束DistFlow模型的妙处配电网潮流计算有各种方法但做恢复重构这类开关状态优化问题最合适的是DistFlow方程。它专门针对辐射状配电网推导把支路潮流写成两个递推方程$$ P_{ij} - \sum_{k:(j,k)} P_{jk} p_j^d - p_j^g $$ $$ Q_{ij} - \sum_{k:(j,k)} Q_{jk} q_j^d - q_j^g $$ $$ V_j^2 V_i^2 - 2(R_{ij}P_{ij} X_{ij}Q_{ij}) (R_{ij}^2X_{ij}^2)\frac{P_{ij}^2Q_{ij}^2}{V_i^2} $$前两个方程是节点有功、无功平衡第三个方程描述支路两端电压关系。这里有个麻烦第三个方程里有电流平方项是非凸的直接扔给求解器解不动。解决办法有两种。第一种是用二阶锥松弛把电压平方和电流平方定义成新变量将第三个方程松弛为二阶锥约束。这是目前学术界的主流做法松弛后的问题在满足一定条件下能保证全局最优。第二种是用线性DistFlow直接忽略掉损耗项$(R_{ij}^2X_{ij}^2)\frac{P_{ij}^2Q_{ij}^2}{V_i^2}$把方程变成线性的求解速度飞快但精度稍差。我在升级版本里用的是二阶锥形式因为故障恢复场景下线路潮流变化大损耗项不能简单忽略。如果你只是想做快速评估线性DistFlow也够用但要注意电压结果在重载情况下会偏乐观。2.3 辐射状拓扑约束最容易写错的地方配电网正常运行必须保持辐射状结构这是硬约束。数学上辐射状有两个充要条件所有节点连通支路数等于节点数减去子图数量。如果已知故障隔离后系统分成m个独立子网包含主网和各个孤岛那么支路数等于节点数减m。但“支路数等于节点数减子图数”只是必要条件还得加上连通性约束才能保证每个子图内部连通。我调试时踩过这个坑只写了支路数约束模型算出来一个“断开但无法从任何电源供电”的孤立节点离谱的是它完全满足支路数条件。连通性约束的常见实现方法是单商品流法。在每条支路上定义一个虚拟流变量$f_{ij}$从电源节点向负荷节点流动约束条件如下$$ \sum_{j:(i,j)} f_{ij} - \sum_{j:(j,i)} f_{ji} \le 1, \quad \forall i $$ $$ |f_{ij}| \le M \cdot s_{ij}, \quad \forall (i,j) $$第一个约束保证每个节点最多产生一个单位的虚拟流表示每个负荷最多由一个电源供电第二个约束保证支路断开时虚拟流也为零。加上这个约束后所有带电节点就必然能连接到某个电源节点拓扑连通性有了保证。注意$M$的取值。$M$太大比如$10^6$会破坏求解器的数值稳定性太小又会错误地截断可行解。实际调试时取节点数加1就够了33节点系统取34左右就行。2.4 孤岛划分约束功率平衡是底线孤岛能不能成立最关键的是功率平衡。孤岛内所有负荷的功率之和必须小于等于孤岛内所有分布式电源和储能的出力之和。写成约束$$ \sum_{i \in \Omega_g} P_{i,t}^g \sum_{i \in \Omega_s} P_{i,t}^{dis} \ge \sum_{i \in \Omega_g} (w_i \cdot P_{i,t}^{load} \cdot x_{i,t}) P_{loss,t} $$这里$\Omega_g$表示属于该孤岛的节点集合$\Omega_s$表示孤岛内储能节点集合。$P_{loss,t}$是孤岛网络损耗如果不方便精确计算可以用一个比例系数粗略估计比如取孤岛总负荷的5%。这里有个实现细节孤岛的功率平衡约束是“条件约束”只有孤岛成立时才生效。怎么处理引入孤岛状态变量$z_{g,t}$1表示节点g属于某孤岛然后用大M法把条件约束写成$$ \sum_{i \in \Omega_g} P_{i,t}^g \sum_{i \in \Omega_s} P_{i,t}^{dis} M(1-z_{g,t}) \ge \sum_{i \in \Omega_g} P_{i,t}^{load} \cdot x_{i,t} $$当$z_{g,t}1$时这个约束退化为正常的功率平衡当$z_{g,t}0$时$M$项让约束自动松弛不产生影响。还要注意孤岛内必须至少有一个分布式电源或储能否则该孤岛无法建立电压和频率。这个条件在辐射状约束和连通性约束下其实能自动满足但加上显式约束会更稳妥$$ \sum_{i \in \Omega_g} u_{i}^{DG} \ge z_{g,t} $$其中$u_i^{DG}$是分布式电源是否投入运行的变量。3. Matlab实现从公式到能跑的代码3.1 算例数据怎么组织做配电网重构研究最常用的基准算例就是33节点系统这个系统常被称为“标准33节点算例”在文献里出现频率极高。它包含33个节点、37条支路其中有5条联络开关支路系统额定电压12.66kV总负荷约3.7MW。基础版本可以直接用这个系统的原始数据升级版本需要改造加入几个分布式电源和储能节点。数据组织方式建议用结构体数组每个节点和每条支路都是结构体的一个元素。这样写代码时语义清晰调试也方便。基本字段如下% 节点数据结构 node(id).bus_type % 1-PQ节点 2-PV节点 3-平衡节点 node(id).P_load % 有功负荷 kW node(id).Q_load % 无功负荷 kvar node(id).level % 负荷等级 1/2/3 % 支路数据结构 branch(id).from % 支路首端节点 branch(id).to % 支路末端节点 branch(id).R % 电阻 Ohm branch(id).X % 电抗 Ohm branch(id).s0 % 初始开关状态 0-断开 1-闭合 % 分布式电源数据结构 dg(id).node % 接入节点 dg(id).P_max % 最大有功出力 kW dg(id).P_min % 最小有功出力 kW % 储能数据结构 ess(id).node % 接入节点 ess(id).E_init % 初始电量 kWh ess(id).E_max % 额定容量 kWh ess(id).P_charge_max % 最大充电功率 kW ess(id).P_discharge_max % 最大放电功率 kW用结构体数组的好处是传递方便、查看某个字段时一目了然。我用表格存储过一段时间写约束时访问列数据确实方便但每个字段都要用表格变量名引用代码会显得冗长。对比下来结构体更合适。3.2 建模工具箱与求解器配置Matlab里做混合整数规划建模最常用的工具箱是YALMIP其次是直接调优化工具箱自带的intlinprog。我个人的经验是除非模型特别简单否则别用纯intlinprog方式写因为你把自己所有的约束都手动转成矩阵形式那工作量太大了而且容易出错。YALMIP的好处是可以用接近数学公式的方式写约束代码可读性强排错方便。变量定义用binvar或sdpvar约束用方括号组合求解调optimize接口。下面这段是我定义决策变量和核心约束的代码框架% 决策变量定义 % s_branch: 支路开关状态变量 (N_branch x T) 二进制变量 % x_node: 节点带电状态变量 (N_node x T) 二进制变量 % P_dg: 分布式电源有功出力 (N_dg x T) 连续变量 % P_ess_in / P_ess_out: 储能充放电功率 (N_ess x T) 连续变量 % soc: 储能荷电状态 (N_ess x T1) 连续变量 s_branch binvar(N_branch, T, full); x_node binvar(N_node, T, full); P_dg sdpvar(N_dg, T, full); P_ess_in sdpvar(N_ess, T, full); P_ess_out sdpvar(N_ess, T, full); soc sdpvar(N_ess, T1, full); % 辐射状约束支路数 节点数 - 子图数 % 简化处理假设故障隔离后系统分为 M_sub 个子网 Constraints [Constraints, sum(s_branch(:, t), 1) N_node - M_sub];这一段只展示了变量的定义和一个拓扑约束骨架。完整的约束集合还包括潮流约束的逐条展开、电压上下限约束、DG出力上下限、储能充放电功率上下限和SOC递推式等。求解器方面YALMIP本身不求解问题需要一个后端求解器。Matlab自带的intlinprog可以求解中小规模的混合整数线性规划问题但如果你用了二阶锥松弛intlinprog就不够用了。这种情况下推荐使用支持混合整数二阶锥规划的求解器。这里有个经验在配置求解器时务必设置MIPGap容差和最长求解时间。故障恢复问题一旦扩大到多时段二进制变量数量很容易超过一千个如果不设置求解时间上限求解器可能几个小时都跑不完。我一般在求解前做三件事设置MIPGap为0.01%保证解的质量接近最优即可不必强求0。设置求解时间上限为600秒超过时限取当前最优可行解。打开求解日志输出方便观察求解进度判断模型是否存在数值问题。3.3 拓扑连通性约束的YALMIP写法单商品流法连通性约束在YALMIP里写起来很直接。引入虚拟流变量f_branch和源节点变量source_node然后逐条添加约束% 虚拟流变量 f_branch sdpvar(N_branch, T, full); % 连通性约束虚拟流不能通过断开支路 Constraints [Constraints, -M * s_branch(:, t) f_branch(:, t) M * s_branch(:, t)]; % 每个节点虚拟流守恒约束 % A_node_branch是节点-支路关联矩阵A_node_branch[i,j]1表示支路j从节点i流出-1表示流入 for i 1:N_node inflow_expr A_node_branch(i, :) * f_branch(:, t); Constraints [Constraints, inflow_expr 1]; end关联矩阵A_node_branch的构造有个小技巧遍历每条支路把支路首端节点对应的元素设为1末端节点对应的元素设为-1。不需要考虑左右端方向与实际潮流方向一致因为配电网重构问题中的潮流方向本来就不固定只要保证关联矩阵一致即可。写约束时还要注意故障隔离后主网电源节点和孤岛的分布式电源节点都可以作为节点产生虚拟流的“源”。最简单的办法是让所有电源节点都能产生至多一个单位的虚拟流这样每个负荷节点最多被一个电源覆盖与孤岛划分的物理意义一致。3.4 结果后处理与可视化怎么做求解完成后严格来说先别急着分析结果第一步是验证解是否满足所有约束。我习惯用YALMIP提供的check函数% 检查约束违背情况 diagnostics check(Constraints); % 输出最大约束违背量 max_violation max(diagnostics);如果max_violation大于1e-5说明解有问题即使是求解器返回了“成功”标志也要警惕。有时候求解器因数值问题返回了不严格满足约束的“伪最优解”这一步能帮你及时发现。拓扑可视化用Matlab的Graph对象最方便% 根据最优开关状态绘制恢复后拓扑 G graph([], []); for j 1:N_branch if value(s_branch(j, 1)) 0.5 % 开关闭合 G addedge(G, branch(j).from, branch(j).to); end end figure; plot(G, Layout, force);把恢复后的网络结构画出来能直观地看出孤岛怎么划分的、哪些支路断开了。我这个升级版本还额外输出两张图各时段电压剖面图和恢复负荷柱状图。电压剖面图能直观看出有没有电压越限的节点恢复负荷柱状图则展示不同时段各类负荷的恢复比例写论文时这两张图都很有说服力。4. 升级版本的难点动态孤岛与多时段耦合4.1 多时段模型中的储能和动作次数约束多时段模型和单时段的本质区别在于引入了“状态转移”约束。储能是最典型的例子。储能的荷电状态变化满足递推关系$$ SOC_{t1} SOC_t \eta_c \cdot P_{in,t} \cdot \Delta t - \frac{P_{out,t} \cdot \Delta t}{\eta_d} $$这个式子本身是线性的但要小心两个细节。一是充放电功率不能同时为正否则模型会利用这个漏洞“免费”充电放电导致SOC虚增。处理办法是引入两个二进制变量表示充放电状态或者用互补约束$P_{in,t} \cdot P_{out,t} 0$的线性化形式。二是SOC本身有上下限约束。故障恢复期间储能一般不充到100%也不放到0%因为要保留一定的功率调节能力。升级版本里我把SOC范围限制在20%到90%这样既保护了储能设备寿命也保留了响应裕度。开关动作次数约束是多时段模型另一个关键点。直接统计所有开关在所有时段的总动作次数会得到一个二次表达式$|s_{t}-s_{t-1}|$之和其中绝对值是非线性的。处理办法是引入辅助二进制变量$action_{ij,t}$表示开关是否动作$$ action_{ij,t} \ge s_{ij,t} - s_{ij,t-1} $$ $$ action_{ij,t} \ge s_{ij,t-1} - s_{ij,t} $$ $$ \sum_{t} action_{ij,t} \le Action_max $$每次开关动作都会改变潮流路径必须重新计算潮流。这在多时段模型里意味着潮流约束要配着开关状态变量“切换”。我在实现时用了大M法把潮流方程写成条件等式开关闭合时潮流方程生效断开时对应的潮流全部为0。这一步是整个模型里最容易出错的地方一定要反复验证。4.2 分布式电源出力的不确定性怎么考虑升级版里如果要更贴近实际还应该考虑光伏和风电出力的不确定性。最简单的办法是场景法生成一组典型出力场景每个场景赋予一个概率目标函数变成所有场景下的期望恢复量。这样做模型规模会成倍增长求解时间也急剧上升。另一个思路是鲁棒优化假设分布式电源出力落在某个区间内要求在所有可能出力下恢复方案都可行。这种模型更保守恢复效果比期望值模型差但方案的安全性最高。实际工程调度更看重鲁棒性因为恢复方案要保证“一定能执行”而不是“大概率能执行”。我个人的建议是如果你主要做学术研究先用确定性模型把机理讲清楚再扩展不确定性如果你是做工程部署建议直接上鲁棒模型或者两阶段鲁棒模型牺牲一点恢复量换可靠性。4.3 求解效率优化方法多时段统一模型的二进制变量会迅速膨胀。33节点系统、16个时段、37条支路支路状态变量就有592个再加上节点带电状态变量528个总计超过1000个二进制变量。直接求解很容易陷入长时间迭代。我实测有效的几个优化手段固定故障隔离区的开关状态。故障发生后故障点附近的开关必须断开隔离这部分变量可以直接设为固定值不在优化变量范围内能减少几十个二进制变量。预选候选开关集合。故障恢复中真正需要调整的开关往往只在失电区域边界附近。把所有联络开关和失电区边界的分段开关作为候选集其他开关保持原状态可以显著压缩解空间。先求一个可行解作为热启动。可以先用简化线性模型快速求一组开关状态作为完整模型的初始可行解这样求解器的MIP求解过程能更快收敛。这些方法每个都能带来成倍的求解速度提升加在一起一个16时段模型从跑不动变成几分钟内收敛效果非常明显。5. 常见问题与排查技巧实录5.1 求解结果出现环网这个是最常见的问题。明明加了辐射状约束模型输出结果里还是出现了环状结构。大概率是你只写了支路数等于节点数减子图数的必要条件但没有加连通性约束。我前面说过这两个约束缺一不可。排查方法很简单把最优解对应的闭合开关拓扑画出来数一数闭合支路数量。如果闭合支路数等于节点数减子图数但出现了环那一定是连通性约束出了问题重点检查单商品流约束里的关联矩阵符号对不对。5.2 出现“孤立节点”或正常负荷没恢复如果某些负荷节点没有带电但你检查开关状态发现它和电源节点之间的支路都闭合着那问题大概率出在潮流约束里。我遇到过一个情况把所有支路都闭合了但潮流解不出来因为分布式电源出力被设成了固定值0而主网电源又因为某个隔离开关误设为断开无法供电。排查思路先看这个节点的“父节点”方向有没有闭合路径连到电源再检查对应支路的开关变量最优值是不是1还要看该节点所在的孤岛是否满足功率平衡约束。这个三个检查跑一遍大部分问题都能定位。5.3 求解器提示模型不可行模型不可行先别怀疑求解器先怀疑约束之间互相矛盾。最常见的矛盾是负荷恢复量太大、分布式电源容量不足。解决方法是把目标函数里的“强制恢复”改成“尽量恢复”——引入表示负荷是否恢复的松弛变量让模型允许在电源不足时放弃部分三级负荷。另一个容易忽视的原因是孤岛功率平衡约束里的$M$取值不合理。上面说过$M10^6$这种取值在求解器眼里经常真的就是10的6次方会导致数值灾难。我之前用了一个1e6的$M$模型怎么调都无解后来改成$M1000$马上就出结果了。大M参数的选法我总结一个原则M取“该物理量在实际系统中的最大可能值”的1.1倍不要更多。比如支路潮流$P_{ij}$的M取该条线路最大传输容量的1.1倍虚拟流的M取节点总数加1。这样既保证约束松弛效果又不会破坏数值稳定性。5.4 最终恢复方案电压越限恢复后的拓扑满足所有结构约束但某些节点电压低于0.95p.u.甚至低于0.9p.u.这在工程上不可接受。原因是重构和孤岛划分模型里只做了静态潮流校验没有对电压约束设置严格的上下限。解决办法有两个方向。一是把电压约束$V_{min} \le V_i \le V_{max}$直接写进潮流约束中让模型在搜索过程中自动排除电压越限的方案。二是用灵敏度方法做后校验先求拓扑再计算电压如果越限就在目标函数中惩罚越限节点迭代几次也能收敛到可行方案。第二种方法在模型规模大时更实用。5.5 常见问题速查表现象可能原因排查方向求解结果出现环网缺少连通性约束或关联矩阵错误检查单商品流约束和A_node_branch矩阵出现孤立带电节点电源节点设置错误或隔离开关误开检查电源节点集合和故障支路开关状态模型无解大M参数过大、负荷恢复要求过强调整M取值、引入负荷恢复松弛变量电压越限电压约束未入模或DG容量不足添加上下限约束、增加无功或储能支撑求解时间过长二进制变量过多、MIPGap太小预选候选开关、设置时限和MIPGap恢复方案开关动作过多目标中开关惩罚权重太小增大开关动作惩罚系数λ6. 实操心得与复盘把标准33节点算例配上分布式电源和储能改造成含DG版本升级版的模型从单时段扩展到16个时段整体代码量大概在800行左右。这个规模对于个人研究者来说不算小但结构清晰的话调试起来反而比一堆散乱的脚本更省心。我个人实际跑下来最大的体会是配电网故障恢复这种混合整数问题千万别一上来就上完整商业模型。先把单时段的33节点系统调通所有约束全部验证无误再往上加多时段、储能SOC、不确定性一步一步来。每个阶段都保存一版能跑的代码出了问题可以随时退回上一版对比。还有一个容易被忽略的点求解器的日志一定要打开看。很多时候结果不对不是因为模型错了而是求解器在中途因为数值问题提前终止或者找到了一个“次优可行解”但质量很差。看日志能发现这些隐藏信息节省大量排查时间。最后分享一个小技巧如果你用的是YALMIP写完约束后先别急着求解可以用simplify合并同类项再用coefficients提取约束矩阵看看约束里有没有明显的零行或重复行。这个习惯帮我抓出过不少因为变量引用错误导致的冗余约束。把这套流程走通你可以顺手把模型扩展到更大规模的配网系统或者更细的负荷模型核心框架不用大改收益却会成倍放大。