配电网辐射状拓扑约束的断线解环建模与MATLAB实现

做配电网优化的人第一关往往不是目标函数恰恰是“怎么把网络必须是辐射状的这一条写成约束”。这几天我在用 MATLAB 复现一篇 EI 源论文里“基于断线解环思想”的配电网辐射状拓扑约束建模方法前前后后踩了不少坑也把不少理论上的老套路重新梳理了一遍。这篇就把我的完整实现过程、数学模型的来龙去脉、代码层面的坑都写出来适合正在做网络重构、网架规划、故障自愈、分布式电源接入优化的工程师和研究生参考。先剧透一个结论只用“闭合支路数等于节点数减一”这条约束是绝对不够的你还需要解决连通性和无环问题而“断线解环”正好提供了一套非常工程化的建模语言。1. 为什么配电网必须辐射状先搞清楚“断线解环”到底在解什么1.1 配电网拓扑的硬性要求输电网一般是环网运行因为潮流可以沿着多条路径流动供电可靠性高。但配电网恰恰相反绝大部分中低压配电网在正常运行方式下必须保持“辐射状”也就是从变电站母线出发每条线路像一棵树一样展开任何一个负荷节点都只有唯一一条供电路径。原因很直接配电网保护基本不依赖方向元件辐射状结构可以保证故障电流方向唯一继电保护配合简单可靠同时短路电流也容易控制不会因为环网导致保护越级跳闸。但实际电网硬件并不是天然树的。很多城市电缆网、园区电网在设计时搞的是“环网接线”平时通过打开联络开关、分段开关把环网拆成几条辐射状馈线来运行。这个“平时把环断开”的操作就是电网里常说的“解环”。调度员解环靠的是现场经验而优化建模里要解的“环”是把所有开关状态统一当作决策变量在数学上保证最后选出来的闭合支路集合是一个无环、连通的树。1.2 “断线解环”如何翻译成数学问题如果你只给优化模型下一条指令“请把网络搞成辐射状”模型是听不懂的。它需要几个明确的数量和结构条件连通负荷节点不能脱离电源也就是电网不能出现电力孤岛无环闭合支路不能构成任何回路基数正确树的性质之一是边数等于点数减一。从图论角度看配电网是图 G(V,E)其中 V 是节点集合E 是联络开关所在的可控支路集合。我们需要从 E 中选出一个子集 T使得 G_T(V,T) 是无向树。但纯图论条件“T 是树”很难直接用线性不等式表达因为“不存在一个子图构成环”是一个全局拓扑性质直接写全是枚举指数级的子图约束。“断线解环”思想就是反过来的既然环网里每一个基本回路都必须被至少断开一条支路那我直接给每个基本回路写一条约束要求这个回路里闭合支路的数量不超过回路支路总数减一。本质上就是在告诉求解器每个可能形成的环都不许闭合。这个思想来自调度员的操作经验但在数学上它等价于对回路空间做割集约束。不过实际操作要注意如果配电网原始拓扑里环的数量很多基本回路数量也会很多。以某公开的33节点配网算例为例初始环网结构下基本回路大概有5个约束数量还不算夸张但如果是一个密集网格型网络基本回路数量可能呈指数增长这时候直接枚举所有回路约束就不可取了。所以目前工程实现里更稳的做法是采用“有向生成树”约束把整棵树拆成“每个节点有且仅有一个父节点”这种局部条件再配合断线解环的基数约束兼顾全局连通性和求解效率。1.3 别再迷信“闭合支路数节点数-1”这个约束几乎是所有初学重构模型的人最先写出来的东西sum(x) n - 1其中 x 是每条支路是否闭合的 0-1 变量n 是节点数。这个约束是对的但单独使用远远不够。因为“边数等于点数减一”只能说明支路数量对的但不能排除两种情况一是网络分成两个甚至多个互不相连的子网整体加起来还是 n-1 条边二是某个子网内部有环同时另一个子网缺边两边凑在一起边数刚好满足。我举一个特别直观的例子。三个节点 1、2、3三条支路分别构成三角形另有根节点0通过两条支路连接1和2。如果我们随便选4条支路闭合可能出现节点3挂在支路1-3上、但节点1和2又通过另外一条支路闭合导致一个三角形子环同时根到1、根到2的支路里有断开。这种状态下边数可能刚好等于节点数减一但拓扑显然不是辐射状。所以要真正实现辐射状约束必须同时把“根节点存在”“每个非根节点有一条入边”“闭合支路与有向边对应”这些条件叠加上去。这也是我复现那篇 EI 论文时感受最深的一点好模型不会只靠一条巧妙的约束而是几条约束的配合才把拓扑结构锁死。2. 核心建模思路拆解断线解环如何写成可求解的约束2.1 决策变量怎么写先定义配电网的可控支路集合。每条支路 k 有首端节点 i、末端节点 j、电阻 r、电抗 x、容量上限 cap 和初始开关状态。在 MATLAB 代码里我习惯用一个结构体存储branch struct(... from, [1 1 2 3 4 5],... to, [2 3 4 5 6 6],... r, [0.1 0.2 0.15 0.2 0.1 0.1],... x, [0.05 0.1 0.08 0.1 0.05 0.05],... cap, [200 200 200 200 200 200]);优化变量包括变量类型含义x(k)二元支路 k 是否闭合1 为闭合y_f(i,j)二元节点 i 是节点 j 的父节点y_t(j,i)二元节点 j 是节点 i 的父节点P(k), Q(k)连续支路 k 上流过的有功、无功U(i)连续节点 i 的电压幅值平方注意我这里把每一条无向支路拆成了两个方向的有向候选父子关系。y_f 和 y_t 并不是独立变量它们必须与 x 关联如果支路闭合那么要么 i 是 j 的父节点要么 j 是 i 的父节点二者只能有一个如果支路断开则两个方向都不能成立。写成约束就是y_f(i,j) y_t(j,i) x(k)这个等式非常关键它把“支路闭合”这个开关变量的含义翻译成了图结构上的“方向选择”。很多初学者直接把 y 和 x 写成“小于等于关系”结果模型在可行域上放松了优化器可能会选中一个没有方向、支路却闭合的拓扑导致后面潮流约束大量报错。2.2 父节点约束把树结构局部化断线解环思想的严格数学表达不需要枚举全部回路只需要为每个节点引入“父节点唯一”约束。设根节点为0要求根节点没有任何父节点每个非根节点有且仅有一个父节点两个节点之间的父节点指向必须借助闭合支路。写成约束组就是sum(y_f(i,root)) 0 sum(y_t(j,i)) 1 对于任意非根节点 i再加上方向变量与开关变量之间的关联约束整个图就自然形成了一个“有向树”。为什么这能同时保证连通和无环推理很直接每个非根节点都有一条入边那就意味着从任意一个非根节点出发沿着父节点方向不断向上回溯必然要么到达根节点要么陷入循环。但如果陷入了循环循环里每一个节点都会有一条来自环内节点的入边而“每个节点有且仅有一条入边”的条件会让这条循环链条之外再找不到新增边来接入根节点整体会出现某个子网没有根的情况这与“根节点无父节点”和“闭合支路数等于 n-1”联合起来是矛盾的。因此整个图最终一定是一棵以根节点为唯一根的有向树。如果你看过一些文献里的表述会发现有人把这套约束叫“单商品流”约束也有人叫“有向生成树”约束。其实本质是同一个东西只是在实现细节上有的地方用 f(i) 表示节点 i 是否从根可达然后写 f(j) f(i) 来传播连通性。我复现的这篇论文里用的就是上述“每个节点有唯一父节点”的版本因为它不需要额外引进一个连续的“虚拟潮流”变量求解规模相对更小数值上也不太容易出现尺度问题。2.3 环网基础回路约束要不要加有人会问既然有向生成树约束已经能保证辐射状那论文标题里的“断线解环思想”体现在哪里答案在加速技巧上。虽然基础约束已经足够但大多数 MILP 求解器在分支定界过程中并不会主动识别“回路闭合”这种全局违反它只能通过一步步分支去试。为了提高求解效率通常还会补一组“基本回路割约束”对每个通过初始闭合开关形成的环路 L要求sum(x(k), k in L) length(L) - 1含义就是每条基本回路至少要有 1 条支路断开也就是“断线解环”。有了这条约束求解器在很早的阶段就能剪掉大量环网分支特别是对 33、69 节点这种中等规模的算例效果非常明显。不过要注意基础回路集合需要在预处理阶段用深度优先搜索或并查集方法生成而且生成的回路应该尽量选择较短的回路。长回路约束对剪枝帮助不大反而增加约束数量。在实际中我的做法是先用图论函数生成所有基本回路然后按回路长度从小到大加入模型只保留长度不超过 6 的短回路约束。这样既保留了断线解环的加速效果也不会因为约束过多拖慢求解器。3. MATLAB 实现要点从数据准备到求解配置3.1 潮流平衡方程的线性化处理拓扑约束只是模型的一半另一半是支路功率和节点电压。如果把精确交流潮流直接写进混合整数模型那就是一个高度非凸的 MINLP普通求解器根本跑不动。论文里常见的做法是采用线性化 DistFlow 模型。对任意闭合支路 k首端节点 i末端节点 j线性化 DistFlow 写为U_j U_i - 2 * ( r_ij * P_ij x_ij * Q_ij )注意这里 U 是电压幅值的平方不是电压本身。为什么取平方而不是直接用电压因为支路功率方程里与电压相关的项是 U^2 形式的直接以 U 为变量会引入非线性取平方后线路压降约束就变成了线性关系。但支路开断时节点 i 和 j 之间不存在电气联系所以这个等式不能强制成立。必须用 big-M 方法将它配成可松弛版本-U_j U_i - 2*(r*P x*Q) -M*(1-x_k) -U_j U_i - 2*(r*P x*Q) M*(1-x_k)当 x_k1 时两边收紧压降等式生效当 x_k0 时压降关系被 M 项放松不影响拓扑。这里的 M 选取极其讲究我后文专门讲踩坑。节点功率平衡约束也要按支路开关状态处理如果支路断开P_ij 和 Q_ij 必须强制为 0如果支路闭合支路功率应该参与节点功率平衡。所以我通常额外加-M*x_k P_k M*x_k -M*x_k Q_k M*x_k这一步很多初学代码的人会漏掉。他们只写了支路断开时的压降放宽却忘了把断开支路上的功率清为零结果求解完以后断开支路仍然在偷偷转移功率整个“辐射状网络”的电气上却是环网。这个 bug 非常隐蔽因为目标函数可能看起来正常但你去检查每条支路的功率时会发现断开支路上有非零潮流。3.2 目标函数如何构建如果目标是网损最小精确表达式是loss sum( r_k * (P_k^2 Q_k^2) / U_k^2 )这个表达式用于 MILP 是不行的因为里面有变量平方和比值。复现论文里一般分两种情况第一种情况是用线性目标做重构寻优例如最小化总网损的线性近似版或最小化开关动作次数。开关动作次数的目标可以直接写Objective sum(x ~ x0);其中 x0 是初始开关状态。不过 Yalmip 里不能直接把不等于关系放进去要拆成两个 0-1 辅助变量。我一般写成switch_on binvar(n_branch,1); switch_off binvar(n_branch,1); Constraints [Constraints, x x0 switch_on - switch_off]; Constraints [Constraints, switch_on switch_off 1]; Objective sum(switch_on) sum(switch_off);这样模型干净求解器也容易处理。第二种情况是真正把网损作为目标那需要引入辅助变量做二阶锥松弛或分段线性近似。这里 Yalmip 支持二阶锥约束可以直接写loss_term sdpvar(n_branch,1); Constraints [Constraints, loss_term(k) r_k*(P_k^2 Q_k^2)]; Objective sum(loss_term);但要注意标准 intlinprog 不支持二阶锥必须调用支持二阶锥的商业求解器。我实测下来模型规模在几十到几百个节点时用这类求解器计算效率尚可但每次求解的耗时波动非常大特别依赖热启动和初始可行解。3.3 Yalmip 代码骨架我习惯用 Yalmip 做建模层因为它的约束拼装非常直观调试也方便。核心代码结构如下% 定义变量 nb size(branch,1); n length(bus(:,1)); x binvar(nb,1); y_f binvar(nb,2); % 第一列表示 from-to, 第二列表示 to-from P sdpvar(nb,1); Q sdpvar(nb,1); U sdpvar(n,1); Constraints []; % 1. 闭合支路数约束 Constraints [Constraints, sum(x) n - 1]; % 2. 每个非根节点有唯一父节点根节点无父节点 root 1; for k 1:nb f branch.from(k); t branch.to(k); Constraints [Constraints, y_f(k,1) y_f(k,2) x(k)]; Constraints [Constraints, y_f(k,1) x(k)]; Constraints [Constraints, y_f(k,2) x(k)]; end % 根节点入边为零 root_in_edges find(branch.to root); Constraints [Constraints, sum(y_f(root_in_edges,2)) 0]; % 其他节点入边为一 for nod 1:n if nod ~ root incoming [find(branch.to nod); find(branch.from nod)]; % 注意branch.to nod 对应的是 to-from 方向的入边 % 需要根据方向变量列数判断 end end这里代码里面的细节容易出错我先说清楚方向对应关系我定义 y_f(k,1) 表示 from 节点是 to 节点的父节点y_f(k,2) 表示 to 节点是 from 节点的父节点。那么节点 nod 的入边来源包括两类支路一类是 to 端等于 nod 的支路此时入边方向为 from-to对应 y_f(k,1)另一类是 from 端等于 nod 的支路此时入边方向为 to-from对应 y_f(k,2)。写成代码时不要搞混否则模型会在不知不觉中被多加了“根节点必须有父节点”或“某节点不能有父节点”这类错误约束。潮流和电压约束再按前文公式逐条加进去最后调用求解器options sdpsettings(verbose, 2, solver, ); sol optimize(Constraints, Objective, options);solver留空时 Yalmip 会自动挑可用求解器。如果机器上没有商业求解器会退回 MATLAB 内置的 intlinprog但节点数超过 100 时intlinprog 的性能会比较吃力。我通常建议至少用支持 MILP 的开源求解器或者直接在另一台机器上装商业求解器并配好 Yalmip 路径。4. 实操复现中的踩坑记录与排查方法4.1 必踩的坑连通性约束丢失导致孤岛我第一次复现时为了简化模型只写了“闭合支路数 n-1”和“每个节点有唯一父节点”但没有把根节点无入边这个约束写清楚。结果优化器输出了一个非常离谱的拓扑根节点本身变成了一个孤岛网络被分成了两个部分其中一部分有环另一部分是完全断开的。由于整体闭合支路数量还是 n-1这个问题从数据上完全看不出来只有把解画成图才能发现。调试这类问题一个特别实用的小工具是在求解完成后把闭合支路抽出来用 DFS 或 MATLAB 的图论函数检查连通分量数量close_branch branch(find(value(x) 0.5), :); G graph(close_branch.from, close_branch.to); components conncomp(G); if length(unique(components)) 1 disp(存在孤岛); end我每次运行完都会执行这个检查这也是复现论文模型时最可靠的“自检”环节。后来我排查发现根节点无入边约束虽然写了但写错了方向变量列号导致实际约束变成“根节点不能有出边”和我的意图正好相反。这类方向性 bug 很难通过看 log 发现必须画网络拓扑确认。4.2 电压约束里 M 值选取不当big-M 在 MILP 里是双刃剑。M 取太小会误伤可行解M 取太大求解器数值条件恶化产生大量数值问题。我一开始为了省事全部用 10000 作为 M结果求解器给出“模型可能不可行”的提示但手动检查发现可行解明明存在。问题的根源是电压平方 U 的量纲一般是 kV^2节点电压在 1.05 到 0.95 之间波动时U 大约是 0.9 到 1.1。而支路功率 P 的量纲可能是 kW 或 MW线路电阻的平方加起来可能只有几百。如果线路是 0.5 欧功率 1000 kW那么压降项也就是 1000 左右。M 取 10000 以后断开支路的电压松弛范围特别大对求解器来说很多“环网”拓扑会暂时看起来可行分支定界就要花大量时间去排除这些错误的解。我现在的做法是逐条支路计算 M不搞全局统一值。对每条支路M_k 取该支路在极限负载下的最大压降M_k 2 * (r_k * cap_k x_k * cap_k) 0.1;这样既保证断开支路不约束压降又不至于让 M 过大。实测下来求解时间能缩短 20% 到 50%数值稳定性也明显改善。4.3 目标函数非线性导致求解器报错很多人复现 EI 论文时发现论文摘要里写着“网损最小”就不假思索地把目标函数写成精确网损公式然后丢给 intlinprog。intlinprog 只支持线性目标遇到二次项直接报错。解决办法有两个一是改用最小化“支路功率平方之和”的 SOC 形式前提是用支持二次目标或二阶锥的求解器二是用分段线性化把二次曲线近似成线性函数。如果你只是做网络重构、不关心网损的精确数值线性近似通常也能得到相近的开关方案。我的经验是先跑一个目标函数为“开关动作次数最小”的版本验证所有约束都可行再换回非线性目标。这样可以把拓扑约束的问题和非线性求解的问题分开排查不会一上来就遇到一堆报错根本分不清是模型错还是求解器不支持。4.4 常见问题速查表现象可能原因处理办法求解器提示不可行根节点无入边约束方向写错检查方向变量列号画拓扑确认结果存在孤岛连通性约束缺失补上唯一父节点约束和根节点无入边约束断开支路上功率非零缺少功率置零约束加 -Mx P Mx求解时间过长缺少环网基本回路约束或 M 值过大加入断线解环约束逐支路设置 M目标函数报错使用非线性目标但求解器不支持换支持二次目标的求解器或线性化电压越限但模型显示可行big-M 太大掩盖了电压约束将 M 按支路容量和阻抗重新计算同一个模型多次求解结果不同求解器未设置随机种子或初始可行解设置固定随机种子采用热启动节点功率平衡总是不满足忽略了根节点所在的平衡节点约束检查根节点的功率平衡方程是否写完整5. 求解性能优化与模型扩展5.1 热启动带来的巨大收益辐射状拓扑约束的 MILP 模型求解难度在于 0-1 变量之间的组合关系。有一部分变量如父节点变量和开关变量强耦合求解器的分支定界会陷入比较长时间的试探。热启动在这种情况下效果非常显著。你可以用图论方法快速生成一棵初始生成树作为可行解比如用 Prim 或 Kruskal 算法按支路容量从大到小选边得到一组 x 初值然后把它作为整数变量的初始解传给求解器。我记得对某个几十节点的算例无热启动时求解器花了 300 多秒还在 65% 的 gap 附近徘徊给了热启动以后几十秒内就收敛到 1% 以内的 gap。这个差距完全是实际工程效率的差距做研究可能还能等做在线自愈或规划方案比选等不起。Yalmip 里可以这样设置初始解assign(x, x0); options sdpsettings(usex0, 1); sol optimize(Constraints, Objective, options);但要注意usex0 并不是所有求解器都支持。intlinprog 支持某些求解器也支持但有的开源求解器会忽略这个选项。最稳妥的办法是直接在求解器层面设置初始解结构。5.2 环网基本回路约束的工程价值前面提到的断线解环约束实际价值是在求解过程起始阶段就给求解器打下“不许成环”的底子。单纯的父节点约束虽然数学上已经能表达辐射状但属于一种“隐式约束”求解器要等变量分支到一定程度才能看出环被闭死。而基本回路约束是“显式约束”每一条回路都被直接限制早期剪枝效果非常明显。生成基本回路时我一般先用初始环网生成一个生成树然后每加入一条非树边就得到一个基本回路。例如对 33 节点系统原始环网中有 37 条支路和 33 个节点生成树有 32 条支路非树边有 5 条因此对应 5 个基本回路。这 5 个回路的约束数量很少却能把求解器的搜索空间切割掉一大半。如果你遇到的是大规模网格网络基本回路数量会多得多。这时候就不要全加进去了而是挑最短的若干条回路加入。实际经验是回路长度超过 6 的约束基本帮不上忙还会增加 LP 松弛的时间。控制在 10 条以内比较合适。5.3 从网络重构扩展到分布式电源优化辐射状拓扑约束建模不仅是网络重构的专利很多配电网规划、分布式电源选址定容、微电网并网方案问题里都要用到同一套拓扑约束。区别只在目标函数和额外变量不同。例如做分布式光伏接入优化时需要在节点上增加安装容量变量目标函数变成光伏容量收益减去网损同时还要考虑电压越限。这时拓扑约束可以保持原封不动只需把节点功率平衡方程改成包含光伏注入P_in - P_out P_load - P_pv这种模块化扩展能力正是“断线解环”建模思路的重要价值把拓扑结构约束变成一个独立的模块后续不管换什么应用场景都不需要重新推导拓扑条件。5.4 双层规划与鲁棒优化的天然接口论文题目里出现的“EI复现”往往不只是跑通单层模型还涉及双层规划或者多场景鲁棒优化。拓扑变量在外层连续功率变量在内层。如果每一层都完整求解 MILP计算量会非常大。我个人的实践建议是外层用断线解环约束和基本回路约束生成候选拓扑内层用连续潮流求解器计算该拓扑下的网损和电压然后通过割平面方式把不可行场景反馈回外层。这种方法不需要每次迭代都求解一个完整 MILP实际工程可操作性更强。复现时如果发现论文模型和解法步骤对不上先检查是不是用了 Benders 分解或列约束生成这种框架不用硬套单层 MILP。另外给网上找的代码做清理时最好先把所有散落的约束写成一个函数文件并且手动构造一个简单的三节点系统做单元测试。三节点系统只需要几条约束可以用穷举法验证模型选择出的拓扑是否合理。如果三节点都过不了就不要急着跑 33 节点或更大的算例。收尾时的一点私人经验我这次复现最大的体会是辐射状拓扑约束的难点从来不在某个单一公式上而在于多个结构条件之间的配合。父节点唯一性、根节点无父节点、支路闭合方向一致、基本回路断开的数量约束这些条件看起来重复实际上各自承担了不同的职责有的保证无环有的保证连通有的负责与潮流耦合。少写一条结果就会变成看似正确但一画图就露馅的伪辐射状。如果你接下来要自己动手实现类似模型我建议按这个顺序来先跑通三节点小系统验证拓扑约束再切换到标准中压算例加入基本回路约束加速最后再把潮流和功率平衡逐条加上去。每个阶段都把求解结果画成网络图检查一次不要轻信求解器的状态为 solved 就认为模型没问题。最后再分享一个小技巧在 MATLAB 里把求解结果导出为节点邻接矩阵后可以顺便计算一下所有负荷节点的供电距离。辐射状网络中任意节点到根节点的路径是唯一的供电距离应该是一个确定值。如果某个节点计算结果距离为无穷大不用怀疑一定是约束模型出了问题。这个方法比盯着数值日志排查快得多。