简介在配电网重构与网架优化中辐射状拓扑是最基本的运行要求但如何将“无环且连通”这一工程直觉转化为可在整数线性规划中表达的数学约束始终是建模难点。针对带环配电网断线解环思想通过选择待断开支路将原始环网还原为树状结构然而仅依靠基本环约束往往无法彻底排除漏网环。为此引入虚拟潮流单商品流约束与支路总数条件共同构成完备的辐射状约束体系。该建模方法以图论中的生成树与环空间为理论基础结合IEEE 33节点算例和matlab整数规划求解可有效支撑配电网重构、故障恢复及分布式电源选址定容等场景。本文从拓扑约束的数学本质出发完整演示约束矩阵构建、求解与结果校验为工程实践和论文复现提供可直接借鉴的实现路径。1. 配电网重构为什么卡在辐射状约束做配电网重构这几年我最深的体会是真正难倒人的往往不是潮流计算而是怎么把“必须是辐射状”这句大白话写成能塞进优化模型的数学约束。很多人一开始写重构模型注意力全放在网损目标、电压约束、DG出力上拓扑约束随手用一句“断开所有环”带过结果求解出来的拓扑要么带环要么有孤岛根本没法用。配电网跟输电网最大的区别就在这里。输电网允许环网运行甚至靠环网提高可靠性配电网如果有合环继电保护会失去选择性故障电流分配混乱可靠性反而下降。所以配电网重构、故障恢复、网架规划这类问题里辐射状是最基本的硬约束。可“无环”和“连通”这两件事天然是组合性质的不是简单的线性等式能表达。最怕的就是模型优化完了还得人工去调开关那还不如直接启发式搜索。最近我复现了一篇基于“断线解环”思想的辐射状拓扑约束建模论文用matlab从零搭了一遍。所谓断线解环理解起来很直白先假设配电网里所有开关全部闭合这时候网络是带环的然后通过“断掉某些线路”把环一个个解掉最终剩下一棵连通的树。这个思路很多文献都用过但真正落到线性规划约束里坑比想象中多。我一开始只加了基本环约束结果解出来的拓扑照样带环折腾了很久才把完整的约束体系补齐。这篇文章就完整复现一遍这个建模过程。我会以IEEE 33节点配电网为例给出matlab代码拆解基本环的生成、约束矩阵的构建、intlinprog求解以及结果校验。适合正在做配电网重构、分布式电源选址定容、或想要复现同类论文的同学参考。废话不多说直接进入正题。2. 断线解环思想的本质把“树”拆成“断哪几条边”要理解断线解环先得清楚配电网原始网络是什么形态。一个典型的配电网把所有分段开关和联络开关都闭合之后就变成了一个连通图G(V,E)其中V是节点集合E是支路集合。这个图不再辐射而是带环的。我们的目标是从E中选出一个子集F使F构成的图仍然连通且没有环也就是生成树。2.1 环数、支路数与基本环根据图论基本公式一个连通的图G如果节点数是|V|支路数是|E|那么它的环数由circuit rank决定如果G是一棵树则|E||V|-1。如果G有环则|E||V|-1。需要断开的支路数量为|E| - |V| 1。这个值就是图中基本环的数量。以IEEE 33节点系统为例总共33个节点37条支路含32个分段开关和5个联络开关所以需要断开的支路数是37-3315条。换句话说为了从环网变成辐射状树恰好得断开5条闭合支路断开多了某个负荷节点就会失电断开少了拓扑就带环。基本环怎么找任意取出一棵生成树树上有|V|-1条支路剩下|E|-|V|1条支路是非树支路。每加一条非树支路就会和树中的唯一路径形成一个环称为基本环。所有这些基本环的集合构成环空间的一组基。断线解环的基本想法就是要求这组基中的每一个基本环都不能被完整保留下来。2.2 基本环约束断线解环的数学表达引入二进制变量x_e表示支路e是否闭合1为闭合0为断开。对于一个基本环C_k如果闭合支路把C_k中所有支路都保留了那么这个环就不可能被解开。因此要让每个环至少断一条线需要满足[ \sum_{e \in C_k} (1 - x_e) \ge 1 ]等价于[ \sum_{e \in C_k} x_e \le |C_k| - 1 ]这就是“断线解环”最核心的线性约束。写成矩阵形式就是Aineq * x bineq每行对应一个基本环。如果你只是随手复现以为加上这个约束再加一个总闭合支路数等于N-1就万事大吉那就错了。2.3 为什么只加基本环约束不够我在第一次复现时就是只加了基本环约束和支路总数约束然后让intlinprog求解一个随机目标结果得到了一组“看起来”满足条件的解。把它画出来一看图上明显还有一个大圈但检测一下发现每个基本环约束都满足。为什么因为基本环的集合虽然是环空间的基但基本环中每条边都断一条并不等于所有环都被破坏。举个例子一个完全图K4有6条边4个节点圈数等于3。如果选择断开的3条边恰好避开某一个三角形那剩下的图可能还是会包含一个由其它边构成的环。更简单地说基本环约束只是必要条件不是充分条件。它能排除大多数环但不能排除由多个基本环的对称差组合出来的新环。配电网网络规模越大、环路越多这种“漏网之环”出现的概率就越大。所以只靠“断线”本身还不够还得把“连通”这个条件显式地建进去。这也是这篇复现里最重要的一个补丁。3. 完整约束构件的设计环约束 虚拟潮流连通性既然基本环约束不够就得想别的办法把连通性约束补上。在整数规划里最常用的手段是引入虚拟潮流也叫单商品流约束。思路是给网络定一个根节点比如配电网的电源母线然后让根节点向所有其他节点发送“虚拟功率”每个非根节点消耗一个单位。因为只有闭合支路才允许虚拟潮流通过所以如果某个节点和根节点之间没有连通的闭合路径那么它的虚拟功率需求就不可能被满足。这样连通性就变成了流平衡方程。3.1 环约束容易踩的坑很多人第一次写这个模型会把所有可能环全部枚举出来每个环都加一条“至少断开一条”。这个思路理论上没问题但配电网支路一多环的数量会爆炸约束矩阵会很大求解效率很低。而用基本环可以减少约束数量但刚才说了基本环约束不是充分的。另一个坑是如果你使用的生成树不同得到的“基本环”集合也不同但环约束的实质效果应该一样。然而在数值求解中不同的环集合对求解器友好程度差别很大有些环大小悬殊约束强弱不同会影响求解时间。后面我会讲到怎么在matlab里用minspantree来获得一个比较规整的基本环集合。3.2 用单商品流补足连通性单商品流的模型长这样。设根节点为r对于每条支路e(i,j)引入两个非负变量f_{ij}和f_{ji}表示从i到j和从j到i的虚拟潮流。每个节点的流量平衡方程为对于根节点r[ \sum_{j \in N(r)} f_{rj} - \sum_{j \in N(r)} f_{jr} N - 1 ]对于非根节点k[ \sum_{j \in N(k)} f_{kj} - \sum_{j \in N(k)} f_{jk} -1 ]这里的N-1表示根节点发出N-1个单位的虚拟功率每个非根节点消耗掉1个单位。因为总发出等于总消耗所以平衡。支路容量约束必须和开关变量x_e耦合[ 0 \le f_{ij} \le M x_e, \quad 0 \le f_{ji} \le M x_e ]这里M是一个足够大的数理论上取N-1就行。因为虚拟潮流总量不会超过N-1取更大的M只会让求解器产生数值灾难没有必要。如果x_e0该支路两个方向的虚拟潮流都为0说明这条支路断开不能传递功率。如果x_e1虚拟潮流可以在两个方向上流动但受容量上限限制。有了这个约束闭合支路图一定是连通的。为什么因为任何非根节点要得到1个单位虚拟潮流必须存在一条从根到它的闭合路径否则流量平衡无法成立。根节点到所有节点都有路径这就是连通性。3.3 与支路总数约束的配合单商品流只保证连通还不能保证无环。因为如果闭合支路数超过N-1图仍然是连通的但一定有环。所以必须再加上总闭合支路数约束[ \sum_{e \in E} x_e N - 1 ]这条约束其实在断线解环里就是“正好断开|E|-N1条线”。结合单商品流连通的图边数等于节点数减1那就必然是一棵树不存在任何环。这是图论里的基本定理连通且边数N-1的图一定是树。三个部分放在一起才算一个完整的辐射状拓扑约束基本环约束确保每个基本环至少断一条线这是断线解环的直接体现。支路总数约束确保闭合支路数等于N-1不会有多余的环。单商品流连通约束确保所有节点都在同一个连通分量里杜绝孤岛。三者缺一不可。我复现的时候把这套约束写进intlinprog才真正拿到一个在任何初始拓扑下都能保持辐射状的通用模型。4. matlab复现从IEEE33节点网络开始下面进入代码部分。我用的系统是IEEE 33节点配电网这算是配电网研究里的“hello world”。整个系统33个节点37条支路其中32条是分段开关5条是联络开关需要断开的支路数就是5。节点编号从1到33支路按标准数据顺序排列。4.1 网络数据与邻接矩阵在matlab里我习惯用一个结构体存网络数据字段包括start_node、end_node、resistance、reactance、is_switch等。这里不贴全部33节点数据了只给出结构示例。% 支路数据 % 每一行: 起点, 终点, 电阻(ohm), 电抗(ohm) branch_data [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; % ... 省略中间支路 32 33 0.3410 0.5302; % 联络开关支路 21 8 2.0 2.0; 9 15 2.0 2.0; 12 22 2.0 2.0; 18 33 2.0 2.0; 25 29 2.0 2.0; ]; start_node branch_data(:,1); end_node branch_data(:,2); E size(branch_data, 1); N max(max(start_node), max(end_node)); % 33注意联络开关支路的电阻电抗在实际IEEE33节点数据里是有具体数值的我这里用2.0示意复现时建议查标准数据表。4.2 用matlab的minspantree求基本环matlab在图论方面有graph对象非常好用。我们可以直接用graph构造拓扑然后求最小生成树再遍历非树边来构造基本环。% 构建无向图 G graph(start_node, end_node); % 求一棵最小生成树 T minspantree(G); % 提取树中的边编号 tree_edges findedge(T); % 找出非树边 all_edges (1:E); non_tree_edges setdiff(all_edges, tree_edges); % 初始化基本环集合用一个cell数组 cycles {}; for k 1:length(non_tree_edges) e non_tree_edges(k); u start_node(e); v end_node(e); % 在树T中找u到v的路径 path shortestpath(T, u, v); % 将路径上的支路转换为原网络支路编号 edges_in_path []; for i 1:length(path)-1 a path(i); b path(i1); % 在原始网络中查找这条边对应的支路编号 edge_id find((start_nodea end_nodeb) | (start_nodeb end_nodea)); edges_in_path [edges_in_path, edge_id]; end % 基本环 树路径 非树边 cycles{k} [edges_in_path, e]; end这段代码的关键是findedge和shortestpath的搭配。由于minspantree返回的是T但T的边和原始G的边不一定一一对应所以用findedge可以拿到树边在原图G中的边索引但上面代码直接用了findedge(T)有一点问题findedge(T)返回的是T对象里的边索引不是原图G的索引。更好一点的做法是直接对G求生成树然后根据T.Edges.EndNodes去匹配原边的编号。下面给一个更稳妥的版本% 求生成树返回树边在原图G中的边索引 T minspantree(G); tree_edges findedge(G, T.Edges.EndNodes(:,1), T.Edges.EndNodes(:,2));这样tree_edges就是原网络中的支路编号。然后按同样的方法遍历非树边。4.3 建立整数线性规划模型有了基本环集合之后就可以搭建整数规划模型了。决策变量分为两部分xE个0-1变量表示支路开合。f2*E个连续非负变量表示每条支路两个方向的虚拟潮流。因此总变量个数是 37 2*37 111。如果不用虚拟流只做基本环约束变量就少很多但为了正确性咱们还是老老实实加上。目标函数可以先随便给一个比如最小化断开支路的权重和。这里设定一个随机权重w_e目的是验证约束不是寻找最优重构。w rand(E,1); % 随机权重 % 完整决策变量 z [x; f1; f2]; nvars E 2*E; % 目标函数只跟x有关f的系数为0 f_obj [w; zeros(2*E, 1)]; % 类型x为二进制f为连续 ctype [repmat(B, 1, E), repmat(C, 1, 2*E)]; lb [zeros(E,1); zeros(2*E,1)]; ub [ones(E,1); (N-1)*ones(2*E,1)]; % 虚拟潮流上界 N-1后面的约束矩阵需要仔细构造。下面章节详细拆解。5. 代码逐段拆解约束矩阵、求解与结果校验这一节是重点也是复现过程中最花时间的地方。约束矩阵一旦拼错intlinprog很可能直接给一个不可行解或者求解器报错。5.1 约束矩阵构建我们先构造不等式约束部分。基本环约束写成[ \sum_{e \in C_k} x_e \le |C_k| - 1 ]由于变量顺序是x在前面所以不等式矩阵Aineq的左半部分就是环矩阵右半部分对应f变量都是0。n_cycles length(cycles); Aineq zeros(n_cycles, nvars); bineq zeros(n_cycles, 1); for k 1:n_cycles cycles_k cycles{k}; Aineq(k, cycles_k) 1; % 对x部分 bineq(k) length(cycles_k) - 1; end接着是另一个不等式约束虚拟潮流的容量限制。对每一条支路e有[ f_{ij} \le (N-1)x_e,\quad f_{ji} \le (N-1)x_e ]这个约束里同时包含x和f所以要在Aineq里补上。列索引x的索引是1:Ef_{ij}的索引是E 2e - 1或我们自己定义的有序编号f_{ji}的索引是E 2e。为了让自己不容易搞混建议一开始就定义一个索引映射x_idx (e) e; f_idx (e, dir) E 2*(e-1) dir; % dir1表示起点到终点dir2表示终点到起点然后循环每条支路for e 1:E row n_cycles 2*(e-1) 1; Aineq(row, x_idx(e)) N-1; Aineq(row, f_idx(e,1)) -1; % 表示 f_ij - (N-1)*x 0 bineq(row) 0; row2 row 1; Aineq(row2, x_idx(e)) N-1; Aineq(row2, f_idx(e,2)) -1; bineq(row2) 0; end然后是等式约束。第一个等式是总闭合支路数等于N-1Aeq zeros(N, nvars); % 先预留后面再填 beq zeros(N, 1); row 1; Aeq(row, x_idx(1:E)) 1; beq(row) N - 1;第二个等式是每个节点的虚拟潮流平衡。N个节点每个节点一条等式根节点特殊处理。对节点k流出减流入等于d_k根节点d_r N-1其他节点d_k -1构造时需要把所有与节点k相连的支路的f都找出来。root 1; for k 1:N row row 1; % 找出所有包含节点k的支路 neighbors edges_in_node{k}; % 提前构造 for idx 1:length(neighbors) e neighbors(idx); u start_node(e); v end_node(e); if u k % 流出 f_{uv}: dir1 Aeq(row, f_idx(e,1)) Aeq(row, f_idx(e,1)) 1; % 流入 f_{vu}: dir2从v到u Aeq(row, f_idx(e,2)) Aeq(row, f_idx(e,2)) - 1; elseif v k % 流出 f_{vu}: dir2 Aeq(row, f_idx(e,2)) Aeq(row, f_idx(e,2)) 1; % 流入 f_{uv}: dir1 Aeq(row, f_idx(e,1)) Aeq(row, f_idx(e,1)) - 1; end end if k root beq(row) N - 1; else beq(row) -1; end end注意这里变量约束里f的下界为0所以不会有负流量。这个模型允许一条闭合支路上两个方向同时有虚拟潮流吗在流量平衡中如果某个支路上双向都有正流量就形成浪费最优解一定会避免因为目标函数不受f影响所以可能出现双向流的“不唯一”情况但依然满足连通性约束。如果想更严谨可以加约束 f_{ij} f_{ji} (N-1)*x_e进一步限制双向同时有流。但在实际求解中单向流总是能表示的因此不加也能跑出正确树形。5.2 目标函数与求解拼好Aineq、bineq、Aeq、beq之后调用intlinprog。x是整数变量f是连续变量可以通过intcon指定整数变量索引。intcon 1:E; % 只有x是整数 options optimoptions(intlinprog, Display, final); [z_sol, fval, exitflag] intlinprog(f_obj, intcon, Aineq, bineq, Aeq, beq, lb, ub, options);如果模型正确exitflag应该为1fval是一个小数。解出来后x部分就是开关状态x_sol z_sol(1:E); closed_edges find(x_sol 0.5);把闭合支路画出来马上能看到拓扑结构。5.3 用图理论校验辐射状结果为了确认结果真的是树我写了一个简单的校验函数。使用matlab内置的graph和conncompG_sol graph(start_node(closed_edges), end_node(closed_edges)); is_connected (length(conncomp(G_sol)) 1); num_edges size(G_sol.Edges, 1); num_nodes N; is_tree is_connected (num_edges num_nodes - 1);如果is_tree为true说明辐射状约束完全满足。我第一次跑通时看到这个true长舒一口气。后来我又在多个随机生成的小型网络上测试全部通过证明这套约束是通用的。6. 复现过程遇到的三类问题与解决复现过程不总是一帆风顺。这里记录我遇到的三个问题希望帮你少走弯路。6.1 基本环约束不充分的实证一开始我偷懒没加虚拟流约束只用了基本环约束和支路总数约束然后跑了一个随机目标。结果求解器很快给出一组解x_sol中闭合支路数是32但画图后发现有一个明显的环而且还不连通。我用conncomp一检测才发现图被分成了两个连通分量其中一个分量里有个小环。这种情况在单纯使用基本环约束时经常发生特别是在环与环共享支路、结构比较复杂的配电网络里。所以别迷信基本环约束它只是“解环”的必要手段不是完整约束。6.2 big-M取值导致数值问题另一个问题是big-M。最初我把虚拟潮流容量写成M*x_eM习惯性地取了10000。结果intlinprog跑得很慢有时候还会出现奇怪的不稳定解。后来分析发现在这个问题里虚拟潮流的最大值不会超过N-1也就是31因为根节点最多发出N-1个单位的流量。把M从10000改成31之后求解速度肉眼可见地提升数值稳定性也好很多。记住big-M不是越大越好只要比理论上界略大一点点就够了。6.3 联络开关编号顺序对结果的影响IEEE33节点系统的联络开关通常放在支路数据的末尾但这会导致生成树时非树支路大概率就是这些联络开关。基本环的大小差异会很大有些环包含十几条支路有些环只有三四条支路。对intlinprog来说约束矩阵的行列规模虽然不大但环集不同收敛速度略有不同。我做了个小实验把支路顺序随机打乱后再求生成树得到的环集不同最终最优解可能会不同因为目标函数随机但辐射状约束都能满足。这说明环集不是唯一的你可以选择任意生成树来构造基本环。如果追求稳定可以选一个节点度数比较均衡的生成树但通常影响不大。7. 这个模型后续还能怎么玩辐射状约束模型一旦搭好就不只是能复现论文了它可以作为配电网重构、故障恢复、网架规划等一系列优化问题的公共底座。7.1 扩展到重构优化最常见的扩展是把目标函数改成网络损耗最小。以配电网重构为例在确定开关状态后流经各支路的电流就是由拓扑决定的。如果要做精确的潮流重构目标函数会变成非线性intlinprog就不再适用得用其他求解器如YALMIP非线性求解器。不过如果把Distflow线性潮流和这个辐射状约束放在一起就能构造一个混合整数二阶锥规划MISOCP或不完全线性化模型用cplex或gurobi求解。7.2 结合Distflow线性潮流Distflow方程是配电网中最常用的支路潮流模型。在辐射状网络下Distflow可以写成三维变量P、Q、U的线性或二阶锥约束。把前面构造的辐射状拓扑约束与Distflow结合就能实现“拓扑与潮流联合优化”。这是当前配电网重构论文的主流建模框架很多高水平论文都是这么写的。断线解环在这里的作用就是提供拓扑约束让潮流方程只在树状网络上生效。7.3 与其他辐射状约束方法的对比除了断线解环虚拟潮流还有几种常用的辐射状约束表达方式。比如“每个非根节点只有一个父节点”的单亲约束它直接把父节点选择纳入决策变量逻辑更直白但需要额外命令节点是有向的。再比如“割集约束”通过枚举所有割集来保证连通约束数量指数级一般用Gomory割平面动态加入。从实际应用看断线解环虚拟流在中小规模配电网中实现简单、求解快、扩展性也不错大规模网络则可以考虑单亲约束配合先进的求解器。我在实际使用中倾向于把断线解环约束作为“快速验证模型逻辑”的第一个版本。它最直观能和图论概念一一对应排查错误也方便。等模型跑通后再根据规模切换到更高级的约束形式。如果你也正在复现类似的配电网约束建模建议先把基本环计算、虚拟流约束、结果校验这三件事做好再往里面叠加潮流和目标函数这样就不会在根上翻车。最后再说一个小技巧每次改约束时都用一个固定的小网络比如把IEEE33节点裁成10个节点先跑通确认无误后再放大网络能省下大量调试时间。本文还有配套的精品资源点击获取 配电网重构的辐射状约束建模：断线解环与虚拟潮流实现
简介在配电网重构与网架优化中辐射状拓扑是最基本的运行要求但如何将“无环且连通”这一工程直觉转化为可在整数线性规划中表达的数学约束始终是建模难点。针对带环配电网断线解环思想通过选择待断开支路将原始环网还原为树状结构然而仅依靠基本环约束往往无法彻底排除漏网环。为此引入虚拟潮流单商品流约束与支路总数条件共同构成完备的辐射状约束体系。该建模方法以图论中的生成树与环空间为理论基础结合IEEE 33节点算例和matlab整数规划求解可有效支撑配电网重构、故障恢复及分布式电源选址定容等场景。本文从拓扑约束的数学本质出发完整演示约束矩阵构建、求解与结果校验为工程实践和论文复现提供可直接借鉴的实现路径。1. 配电网重构为什么卡在辐射状约束做配电网重构这几年我最深的体会是真正难倒人的往往不是潮流计算而是怎么把“必须是辐射状”这句大白话写成能塞进优化模型的数学约束。很多人一开始写重构模型注意力全放在网损目标、电压约束、DG出力上拓扑约束随手用一句“断开所有环”带过结果求解出来的拓扑要么带环要么有孤岛根本没法用。配电网跟输电网最大的区别就在这里。输电网允许环网运行甚至靠环网提高可靠性配电网如果有合环继电保护会失去选择性故障电流分配混乱可靠性反而下降。所以配电网重构、故障恢复、网架规划这类问题里辐射状是最基本的硬约束。可“无环”和“连通”这两件事天然是组合性质的不是简单的线性等式能表达。最怕的就是模型优化完了还得人工去调开关那还不如直接启发式搜索。最近我复现了一篇基于“断线解环”思想的辐射状拓扑约束建模论文用matlab从零搭了一遍。所谓断线解环理解起来很直白先假设配电网里所有开关全部闭合这时候网络是带环的然后通过“断掉某些线路”把环一个个解掉最终剩下一棵连通的树。这个思路很多文献都用过但真正落到线性规划约束里坑比想象中多。我一开始只加了基本环约束结果解出来的拓扑照样带环折腾了很久才把完整的约束体系补齐。这篇文章就完整复现一遍这个建模过程。我会以IEEE 33节点配电网为例给出matlab代码拆解基本环的生成、约束矩阵的构建、intlinprog求解以及结果校验。适合正在做配电网重构、分布式电源选址定容、或想要复现同类论文的同学参考。废话不多说直接进入正题。2. 断线解环思想的本质把“树”拆成“断哪几条边”要理解断线解环先得清楚配电网原始网络是什么形态。一个典型的配电网把所有分段开关和联络开关都闭合之后就变成了一个连通图G(V,E)其中V是节点集合E是支路集合。这个图不再辐射而是带环的。我们的目标是从E中选出一个子集F使F构成的图仍然连通且没有环也就是生成树。2.1 环数、支路数与基本环根据图论基本公式一个连通的图G如果节点数是|V|支路数是|E|那么它的环数由circuit rank决定如果G是一棵树则|E||V|-1。如果G有环则|E||V|-1。需要断开的支路数量为|E| - |V| 1。这个值就是图中基本环的数量。以IEEE 33节点系统为例总共33个节点37条支路含32个分段开关和5个联络开关所以需要断开的支路数是37-3315条。换句话说为了从环网变成辐射状树恰好得断开5条闭合支路断开多了某个负荷节点就会失电断开少了拓扑就带环。基本环怎么找任意取出一棵生成树树上有|V|-1条支路剩下|E|-|V|1条支路是非树支路。每加一条非树支路就会和树中的唯一路径形成一个环称为基本环。所有这些基本环的集合构成环空间的一组基。断线解环的基本想法就是要求这组基中的每一个基本环都不能被完整保留下来。2.2 基本环约束断线解环的数学表达引入二进制变量x_e表示支路e是否闭合1为闭合0为断开。对于一个基本环C_k如果闭合支路把C_k中所有支路都保留了那么这个环就不可能被解开。因此要让每个环至少断一条线需要满足[ \sum_{e \in C_k} (1 - x_e) \ge 1 ]等价于[ \sum_{e \in C_k} x_e \le |C_k| - 1 ]这就是“断线解环”最核心的线性约束。写成矩阵形式就是Aineq * x bineq每行对应一个基本环。如果你只是随手复现以为加上这个约束再加一个总闭合支路数等于N-1就万事大吉那就错了。2.3 为什么只加基本环约束不够我在第一次复现时就是只加了基本环约束和支路总数约束然后让intlinprog求解一个随机目标结果得到了一组“看起来”满足条件的解。把它画出来一看图上明显还有一个大圈但检测一下发现每个基本环约束都满足。为什么因为基本环的集合虽然是环空间的基但基本环中每条边都断一条并不等于所有环都被破坏。举个例子一个完全图K4有6条边4个节点圈数等于3。如果选择断开的3条边恰好避开某一个三角形那剩下的图可能还是会包含一个由其它边构成的环。更简单地说基本环约束只是必要条件不是充分条件。它能排除大多数环但不能排除由多个基本环的对称差组合出来的新环。配电网网络规模越大、环路越多这种“漏网之环”出现的概率就越大。所以只靠“断线”本身还不够还得把“连通”这个条件显式地建进去。这也是这篇复现里最重要的一个补丁。3. 完整约束构件的设计环约束 虚拟潮流连通性既然基本环约束不够就得想别的办法把连通性约束补上。在整数规划里最常用的手段是引入虚拟潮流也叫单商品流约束。思路是给网络定一个根节点比如配电网的电源母线然后让根节点向所有其他节点发送“虚拟功率”每个非根节点消耗一个单位。因为只有闭合支路才允许虚拟潮流通过所以如果某个节点和根节点之间没有连通的闭合路径那么它的虚拟功率需求就不可能被满足。这样连通性就变成了流平衡方程。3.1 环约束容易踩的坑很多人第一次写这个模型会把所有可能环全部枚举出来每个环都加一条“至少断开一条”。这个思路理论上没问题但配电网支路一多环的数量会爆炸约束矩阵会很大求解效率很低。而用基本环可以减少约束数量但刚才说了基本环约束不是充分的。另一个坑是如果你使用的生成树不同得到的“基本环”集合也不同但环约束的实质效果应该一样。然而在数值求解中不同的环集合对求解器友好程度差别很大有些环大小悬殊约束强弱不同会影响求解时间。后面我会讲到怎么在matlab里用minspantree来获得一个比较规整的基本环集合。3.2 用单商品流补足连通性单商品流的模型长这样。设根节点为r对于每条支路e(i,j)引入两个非负变量f_{ij}和f_{ji}表示从i到j和从j到i的虚拟潮流。每个节点的流量平衡方程为对于根节点r[ \sum_{j \in N(r)} f_{rj} - \sum_{j \in N(r)} f_{jr} N - 1 ]对于非根节点k[ \sum_{j \in N(k)} f_{kj} - \sum_{j \in N(k)} f_{jk} -1 ]这里的N-1表示根节点发出N-1个单位的虚拟功率每个非根节点消耗掉1个单位。因为总发出等于总消耗所以平衡。支路容量约束必须和开关变量x_e耦合[ 0 \le f_{ij} \le M x_e, \quad 0 \le f_{ji} \le M x_e ]这里M是一个足够大的数理论上取N-1就行。因为虚拟潮流总量不会超过N-1取更大的M只会让求解器产生数值灾难没有必要。如果x_e0该支路两个方向的虚拟潮流都为0说明这条支路断开不能传递功率。如果x_e1虚拟潮流可以在两个方向上流动但受容量上限限制。有了这个约束闭合支路图一定是连通的。为什么因为任何非根节点要得到1个单位虚拟潮流必须存在一条从根到它的闭合路径否则流量平衡无法成立。根节点到所有节点都有路径这就是连通性。3.3 与支路总数约束的配合单商品流只保证连通还不能保证无环。因为如果闭合支路数超过N-1图仍然是连通的但一定有环。所以必须再加上总闭合支路数约束[ \sum_{e \in E} x_e N - 1 ]这条约束其实在断线解环里就是“正好断开|E|-N1条线”。结合单商品流连通的图边数等于节点数减1那就必然是一棵树不存在任何环。这是图论里的基本定理连通且边数N-1的图一定是树。三个部分放在一起才算一个完整的辐射状拓扑约束基本环约束确保每个基本环至少断一条线这是断线解环的直接体现。支路总数约束确保闭合支路数等于N-1不会有多余的环。单商品流连通约束确保所有节点都在同一个连通分量里杜绝孤岛。三者缺一不可。我复现的时候把这套约束写进intlinprog才真正拿到一个在任何初始拓扑下都能保持辐射状的通用模型。4. matlab复现从IEEE33节点网络开始下面进入代码部分。我用的系统是IEEE 33节点配电网这算是配电网研究里的“hello world”。整个系统33个节点37条支路其中32条是分段开关5条是联络开关需要断开的支路数就是5。节点编号从1到33支路按标准数据顺序排列。4.1 网络数据与邻接矩阵在matlab里我习惯用一个结构体存网络数据字段包括start_node、end_node、resistance、reactance、is_switch等。这里不贴全部33节点数据了只给出结构示例。% 支路数据 % 每一行: 起点, 终点, 电阻(ohm), 电抗(ohm) branch_data [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; % ... 省略中间支路 32 33 0.3410 0.5302; % 联络开关支路 21 8 2.0 2.0; 9 15 2.0 2.0; 12 22 2.0 2.0; 18 33 2.0 2.0; 25 29 2.0 2.0; ]; start_node branch_data(:,1); end_node branch_data(:,2); E size(branch_data, 1); N max(max(start_node), max(end_node)); % 33注意联络开关支路的电阻电抗在实际IEEE33节点数据里是有具体数值的我这里用2.0示意复现时建议查标准数据表。4.2 用matlab的minspantree求基本环matlab在图论方面有graph对象非常好用。我们可以直接用graph构造拓扑然后求最小生成树再遍历非树边来构造基本环。% 构建无向图 G graph(start_node, end_node); % 求一棵最小生成树 T minspantree(G); % 提取树中的边编号 tree_edges findedge(T); % 找出非树边 all_edges (1:E); non_tree_edges setdiff(all_edges, tree_edges); % 初始化基本环集合用一个cell数组 cycles {}; for k 1:length(non_tree_edges) e non_tree_edges(k); u start_node(e); v end_node(e); % 在树T中找u到v的路径 path shortestpath(T, u, v); % 将路径上的支路转换为原网络支路编号 edges_in_path []; for i 1:length(path)-1 a path(i); b path(i1); % 在原始网络中查找这条边对应的支路编号 edge_id find((start_nodea end_nodeb) | (start_nodeb end_nodea)); edges_in_path [edges_in_path, edge_id]; end % 基本环 树路径 非树边 cycles{k} [edges_in_path, e]; end这段代码的关键是findedge和shortestpath的搭配。由于minspantree返回的是T但T的边和原始G的边不一定一一对应所以用findedge可以拿到树边在原图G中的边索引但上面代码直接用了findedge(T)有一点问题findedge(T)返回的是T对象里的边索引不是原图G的索引。更好一点的做法是直接对G求生成树然后根据T.Edges.EndNodes去匹配原边的编号。下面给一个更稳妥的版本% 求生成树返回树边在原图G中的边索引 T minspantree(G); tree_edges findedge(G, T.Edges.EndNodes(:,1), T.Edges.EndNodes(:,2));这样tree_edges就是原网络中的支路编号。然后按同样的方法遍历非树边。4.3 建立整数线性规划模型有了基本环集合之后就可以搭建整数规划模型了。决策变量分为两部分xE个0-1变量表示支路开合。f2*E个连续非负变量表示每条支路两个方向的虚拟潮流。因此总变量个数是 37 2*37 111。如果不用虚拟流只做基本环约束变量就少很多但为了正确性咱们还是老老实实加上。目标函数可以先随便给一个比如最小化断开支路的权重和。这里设定一个随机权重w_e目的是验证约束不是寻找最优重构。w rand(E,1); % 随机权重 % 完整决策变量 z [x; f1; f2]; nvars E 2*E; % 目标函数只跟x有关f的系数为0 f_obj [w; zeros(2*E, 1)]; % 类型x为二进制f为连续 ctype [repmat(B, 1, E), repmat(C, 1, 2*E)]; lb [zeros(E,1); zeros(2*E,1)]; ub [ones(E,1); (N-1)*ones(2*E,1)]; % 虚拟潮流上界 N-1后面的约束矩阵需要仔细构造。下面章节详细拆解。5. 代码逐段拆解约束矩阵、求解与结果校验这一节是重点也是复现过程中最花时间的地方。约束矩阵一旦拼错intlinprog很可能直接给一个不可行解或者求解器报错。5.1 约束矩阵构建我们先构造不等式约束部分。基本环约束写成[ \sum_{e \in C_k} x_e \le |C_k| - 1 ]由于变量顺序是x在前面所以不等式矩阵Aineq的左半部分就是环矩阵右半部分对应f变量都是0。n_cycles length(cycles); Aineq zeros(n_cycles, nvars); bineq zeros(n_cycles, 1); for k 1:n_cycles cycles_k cycles{k}; Aineq(k, cycles_k) 1; % 对x部分 bineq(k) length(cycles_k) - 1; end接着是另一个不等式约束虚拟潮流的容量限制。对每一条支路e有[ f_{ij} \le (N-1)x_e,\quad f_{ji} \le (N-1)x_e ]这个约束里同时包含x和f所以要在Aineq里补上。列索引x的索引是1:Ef_{ij}的索引是E 2e - 1或我们自己定义的有序编号f_{ji}的索引是E 2e。为了让自己不容易搞混建议一开始就定义一个索引映射x_idx (e) e; f_idx (e, dir) E 2*(e-1) dir; % dir1表示起点到终点dir2表示终点到起点然后循环每条支路for e 1:E row n_cycles 2*(e-1) 1; Aineq(row, x_idx(e)) N-1; Aineq(row, f_idx(e,1)) -1; % 表示 f_ij - (N-1)*x 0 bineq(row) 0; row2 row 1; Aineq(row2, x_idx(e)) N-1; Aineq(row2, f_idx(e,2)) -1; bineq(row2) 0; end然后是等式约束。第一个等式是总闭合支路数等于N-1Aeq zeros(N, nvars); % 先预留后面再填 beq zeros(N, 1); row 1; Aeq(row, x_idx(1:E)) 1; beq(row) N - 1;第二个等式是每个节点的虚拟潮流平衡。N个节点每个节点一条等式根节点特殊处理。对节点k流出减流入等于d_k根节点d_r N-1其他节点d_k -1构造时需要把所有与节点k相连的支路的f都找出来。root 1; for k 1:N row row 1; % 找出所有包含节点k的支路 neighbors edges_in_node{k}; % 提前构造 for idx 1:length(neighbors) e neighbors(idx); u start_node(e); v end_node(e); if u k % 流出 f_{uv}: dir1 Aeq(row, f_idx(e,1)) Aeq(row, f_idx(e,1)) 1; % 流入 f_{vu}: dir2从v到u Aeq(row, f_idx(e,2)) Aeq(row, f_idx(e,2)) - 1; elseif v k % 流出 f_{vu}: dir2 Aeq(row, f_idx(e,2)) Aeq(row, f_idx(e,2)) 1; % 流入 f_{uv}: dir1 Aeq(row, f_idx(e,1)) Aeq(row, f_idx(e,1)) - 1; end end if k root beq(row) N - 1; else beq(row) -1; end end注意这里变量约束里f的下界为0所以不会有负流量。这个模型允许一条闭合支路上两个方向同时有虚拟潮流吗在流量平衡中如果某个支路上双向都有正流量就形成浪费最优解一定会避免因为目标函数不受f影响所以可能出现双向流的“不唯一”情况但依然满足连通性约束。如果想更严谨可以加约束 f_{ij} f_{ji} (N-1)*x_e进一步限制双向同时有流。但在实际求解中单向流总是能表示的因此不加也能跑出正确树形。5.2 目标函数与求解拼好Aineq、bineq、Aeq、beq之后调用intlinprog。x是整数变量f是连续变量可以通过intcon指定整数变量索引。intcon 1:E; % 只有x是整数 options optimoptions(intlinprog, Display, final); [z_sol, fval, exitflag] intlinprog(f_obj, intcon, Aineq, bineq, Aeq, beq, lb, ub, options);如果模型正确exitflag应该为1fval是一个小数。解出来后x部分就是开关状态x_sol z_sol(1:E); closed_edges find(x_sol 0.5);把闭合支路画出来马上能看到拓扑结构。5.3 用图理论校验辐射状结果为了确认结果真的是树我写了一个简单的校验函数。使用matlab内置的graph和conncompG_sol graph(start_node(closed_edges), end_node(closed_edges)); is_connected (length(conncomp(G_sol)) 1); num_edges size(G_sol.Edges, 1); num_nodes N; is_tree is_connected (num_edges num_nodes - 1);如果is_tree为true说明辐射状约束完全满足。我第一次跑通时看到这个true长舒一口气。后来我又在多个随机生成的小型网络上测试全部通过证明这套约束是通用的。6. 复现过程遇到的三类问题与解决复现过程不总是一帆风顺。这里记录我遇到的三个问题希望帮你少走弯路。6.1 基本环约束不充分的实证一开始我偷懒没加虚拟流约束只用了基本环约束和支路总数约束然后跑了一个随机目标。结果求解器很快给出一组解x_sol中闭合支路数是32但画图后发现有一个明显的环而且还不连通。我用conncomp一检测才发现图被分成了两个连通分量其中一个分量里有个小环。这种情况在单纯使用基本环约束时经常发生特别是在环与环共享支路、结构比较复杂的配电网络里。所以别迷信基本环约束它只是“解环”的必要手段不是完整约束。6.2 big-M取值导致数值问题另一个问题是big-M。最初我把虚拟潮流容量写成M*x_eM习惯性地取了10000。结果intlinprog跑得很慢有时候还会出现奇怪的不稳定解。后来分析发现在这个问题里虚拟潮流的最大值不会超过N-1也就是31因为根节点最多发出N-1个单位的流量。把M从10000改成31之后求解速度肉眼可见地提升数值稳定性也好很多。记住big-M不是越大越好只要比理论上界略大一点点就够了。6.3 联络开关编号顺序对结果的影响IEEE33节点系统的联络开关通常放在支路数据的末尾但这会导致生成树时非树支路大概率就是这些联络开关。基本环的大小差异会很大有些环包含十几条支路有些环只有三四条支路。对intlinprog来说约束矩阵的行列规模虽然不大但环集不同收敛速度略有不同。我做了个小实验把支路顺序随机打乱后再求生成树得到的环集不同最终最优解可能会不同因为目标函数随机但辐射状约束都能满足。这说明环集不是唯一的你可以选择任意生成树来构造基本环。如果追求稳定可以选一个节点度数比较均衡的生成树但通常影响不大。7. 这个模型后续还能怎么玩辐射状约束模型一旦搭好就不只是能复现论文了它可以作为配电网重构、故障恢复、网架规划等一系列优化问题的公共底座。7.1 扩展到重构优化最常见的扩展是把目标函数改成网络损耗最小。以配电网重构为例在确定开关状态后流经各支路的电流就是由拓扑决定的。如果要做精确的潮流重构目标函数会变成非线性intlinprog就不再适用得用其他求解器如YALMIP非线性求解器。不过如果把Distflow线性潮流和这个辐射状约束放在一起就能构造一个混合整数二阶锥规划MISOCP或不完全线性化模型用cplex或gurobi求解。7.2 结合Distflow线性潮流Distflow方程是配电网中最常用的支路潮流模型。在辐射状网络下Distflow可以写成三维变量P、Q、U的线性或二阶锥约束。把前面构造的辐射状拓扑约束与Distflow结合就能实现“拓扑与潮流联合优化”。这是当前配电网重构论文的主流建模框架很多高水平论文都是这么写的。断线解环在这里的作用就是提供拓扑约束让潮流方程只在树状网络上生效。7.3 与其他辐射状约束方法的对比除了断线解环虚拟潮流还有几种常用的辐射状约束表达方式。比如“每个非根节点只有一个父节点”的单亲约束它直接把父节点选择纳入决策变量逻辑更直白但需要额外命令节点是有向的。再比如“割集约束”通过枚举所有割集来保证连通约束数量指数级一般用Gomory割平面动态加入。从实际应用看断线解环虚拟流在中小规模配电网中实现简单、求解快、扩展性也不错大规模网络则可以考虑单亲约束配合先进的求解器。我在实际使用中倾向于把断线解环约束作为“快速验证模型逻辑”的第一个版本。它最直观能和图论概念一一对应排查错误也方便。等模型跑通后再根据规模切换到更高级的约束形式。如果你也正在复现类似的配电网约束建模建议先把基本环计算、虚拟流约束、结果校验这三件事做好再往里面叠加潮流和目标函数这样就不会在根上翻车。最后再说一个小技巧每次改约束时都用一个固定的小网络比如把IEEE33节点裁成10个节点先跑通确认无误后再放大网络能省下大量调试时间。本文还有配套的精品资源点击获取