做配电网重构、故障恢复或者网架规划研究的同行,应该都有过这样的经历:一篇EI论文里把辐射状约束写得极其漂亮,几行公式就把“断线解环”说得明明白白,可真到用Matlab代码复现的时候,求解器要么给你吐出一堆孤岛,要么干脆散伙不收敛。我当初复现那篇基于断线解环思想的配电网辐射状拓扑约束建模方法时,前面几天基本都在跟“为什么约束都写了,结果还是环网/孤岛”搏斗。这篇就是把那段实操经验完整捋一遍:断线解环思想背后的数学本质是什么,三种常用的辐射状约束建模写法各自坑在哪,以及Matlab里从数据准备到求解验证的完整实现流程。
这篇文章适合正在做配电网拓扑优化、故障恢复策略、分布式电源接入方案的研究生和工程师,也适合那些被“树约束怎么线性化”卡住、想快速用Matlab把模型跑通的人。如果你已经知道辐射状拓扑是什么意思,但不知道约束怎么写才能让求解器真正认账,那这篇文章应该能帮你省下好几天调试时间。
1. 断线解环到底解的是什么?先讲清楚辐射状约束为何难写
1.1 为什么配电网必须“开路运行”:保护、短路电流与调度习惯
很多刚接触配电网优化的人第一反应是:既然电力系统主网是环网运行,配电网为什么非得辐射状?这不是脱裤子放屁吗?答案是:传统配电网的继电保护和运行方式从头到尾都是按“单端供电、树状拓扑”设计的。变电站出口的过流保护、馈线分段开关的定值配合,都假设每一段线路只有一个电源方向。一旦出现环网,短路电流路径不唯一,保护配合关系会彻底错乱,故障隔离范围不可控,这是运行规程层面直接不允许的。
另外,配电网的调度操作习惯也是“开环运行”。实际工作中调度员做负荷转供时,操作票里写的是“先合联络开关、再断开中间某个分段开关”或者反过来,目的就是避免长时间合环运行。所以你去看各种配网重构论文,目标函数可以是网损最小、恢复负荷最多、开关操作次数最少,但几乎每个模型里都塞着同一个硬约束:最终拓扑必须是辐射状,也就是数学上的一棵树。
这就引出了建模时的核心矛盾:物理上要一棵树,但树约束在整数规划里并不是一条简单的线性式子。
1.2 “断开所有环”的直觉建模与隐藏漏洞
断线解环思想的出发点很直观:先把所有候选线路、联络开关都视为闭合,整个网络是一个包含很多环路的连通图;然后我们的任务是选择断开其中一部分线路,把有环的连通图“剪”成一棵辐射状的树。
这个“剪”的过程,最容易写成:每个基本环上至少断开一条线。写成约束就是对于每个基本回路集合 C,有:
[ \sum_{l \in C} (1 - x_l) \ge 1 ]
其中 x_l = 1 表示线路闭合/投运,x_l = 0 表示断开。每个环至少断一条,理论上环确实被破坏了。我第一次看这个约束也觉得已经够了,但拿小算例一跑立刻露馅:这个约束只能保证“无环”,完全不能保证“连通”。
举个极端例子:一个“8”字形网络,两条子环共享一条公共边。只要把中间那条公共边断开,两个基本环的约束都满足了,但网络被切成了两个互不相连的孤岛。每个孤岛内部倒是无环,可这不是配电网要的辐射状——辐射状必须是一整棵连通的树,而不是一片森林。所以断线解环的完整数学条件,必须是“连通 + 边数 = 节点数 - 1”的组合,或者等价地用其他方式同时锁住这两个性质。
1.3 树约束的数学本性:边数、连通性与指数级约束
一个无向图是树的充要条件有两个:一是连通,二是边数等于节点数减一。边数条件是天然线性的:
[ \sum_{k \in L} x_k = N - 1 ]
但连通性条件麻烦得多。通俗地讲,连通意味着“对任意一个节点非空真子集 S,至少有一条线路跨越 S 和 V\S 的割集”。这个约束要枚举所有子集,数量是 2 的指数级,直接写进 MILP 里是不现实的。
所以几乎所有实用方法都在做一件事:用辅助变量把指数级约束转成多项式规模的线性约束,或者用割平面延迟加入。断线解环思想本身并不等于某一条具体公式,而是一类建模策略的总称——它关心的是“怎么把候选的闭环网络变成一棵合法的辐射状树”。下面这节就给出三种我实际验证过能跑通的写法。
2. 三种辐射状约束的数学表达:从生成树到断线解环的完整映射
2.1 边数加连通性:最直观但连通性需要额外变量
最简单的思路就是把树的充要条件逐条翻译成约束。边数约束上面已经写了。连通性我们通常用“每个节点都能从根节点收到流量/物资”来刻画,这就是单商品流思想,我把它单独放在 2.2 节细讲。这一节想强调的是,如果你真的只写 ∑x = N-1,而连通性靠别的约束或者干脆不写,那结果几乎必然是错的。
一个常见翻车现象是:求解器给出的网络边数刚好等于 N-1,但图被分成两个连通分量,其中一个分量的边数恰好等于内部节点数减一。这时 ∑x = N-1 仍然满足,可整体是森林而不是树。我在调试一个配电网故障恢复模型时,就遇到过三角形区域加一个悬空节点的组合,边数对、环也全破了,但负荷节点根本没接到电源上。所以连通性不是锦上添花,而是辐射状约束的半个身子。
2.2 单商品虚拟量法:把连通性写进潮流平衡里
这是工程上最常用、也最容易在Matlab里落地的一种写法。核心思想是:选一个根节点(通常取变电站母线),假设从根节点向全网发出 N-1 个单位的虚拟商品,每个非根节点恰好消费 1 个单位,虚拟商品只能通过“闭合”的线路运输。如果最终拓扑是连通的,那么每个节点都能收到这 1 个单位;如果存在孤岛,孤岛内部不管怎么倒腾,净流入永远凑不出消费量,约束必然被违反。
具体模型如下。设 N 为节点数,L 为候选支路集合,参考方向人为指定(比如从编号小的节点指向编号大的节点)。x_k ∈ {0,1} 表示支路 k 是否闭合;f_k 表示支路 k 上沿参考方向流动的虚拟量,可正可负。约束为:
[ \sum_{k \in L} x_k = N - 1 ]
[ -M x_k \le f_k \le M x_k, \quad \forall k \in L ]
[ \sum_{k \in \delta^+(r)} f_k - \sum_{k \in \delta^-(r)} f_k = N - 1 ]
[ \sum_{k \in \delta^+(v)} f_k - \sum_{k \in \delta^-(v)} f_k = -1, \quad \forall v \ne r ]
其中 δ⁺(v) 表示以 v 为起点的参考方向支路集合,δ⁻(v) 表示以 v 为终点的支路集合。第一条锁边数,第二条锁走向(断开支路不允许有虚拟量),第三、四条相当于每个节点都从根节点“拿到”了消费量,从而锁连通。
这里有两个容易被忽略的细节。第一个是 M 的取值,后面第 4 节细讲。第二个是 f_k 可以正可以负,参考方向只是记账用的,不要把它理解成实际潮流方向,否则后面叠加真实功率平衡约束时容易绕晕。
2.3 基本环割集约束:断线解环思想的理论版
如果把断线解环思想直接公式化,那就是对每个基本环写一条“至少断一条”的约束。前面说了它单独不够,但它也不是完全没有价值——配合连通性约束,它能起到加速作用,在分支切割算法里尤其好用。
我把这种方式称为“环约束 + 连通约束”的混合模式。环约束负责砍掉环,连通约束负责防止砍成孤岛。实际编码时,可以先不加环约束,只跑连通约束和边数约束,求解器也能给出正确的树;但加上环约束之后,可行域被提前剪掉不少坏分支,求解速度往往有肉眼可见的提升。用 MATLAB 的 graph 对象可以直接调 cyclebasis 函数拿到基本环集合,实现成本很低。
三种方法的对比见表:
| 方法 | 辅助变量 | 约束规模 | 求解难度 | 适用场景 |
|---|---|---|---|---|
| 单商品虚拟量法 | L 个连续虚拟量 | O(L) | 中等 | 通用MILP,最推荐首选 |
| 基本环断线约束+连通约束 | 无 | O(环数) | 较快但需防孤岛 | 小规模、拓扑结构简单的网络 |
| 生成树割集约束 | 无 | 指数级,需惰性割平面 | 大规模网络占优 | 节点多、环路多的场景 |
“生成树割集约束”我在这篇里不展开代码了,原理就是前文说的“任意子集至少一条跨割边”,通过分离算法动态加约束。大多数人的项目规模用单商品虚拟量法已经完全够用,先把这条跑通,再考虑高级手段。
3. Matlab建模实操:从支路数据表到可求解的MILP
3.1 数据准备:节点支路表、状态和参考方向
无论你是用手工录入的算例,还是从 Excel/CSV 导入外部数据,第一步都是把网络整理成统一格式。我习惯用矩阵 branch,每一行是一条候选支路,列依次是:首端节点、末端节点、电阻、电抗、初始状态。注意状态这一列必须统一语义:1 表示闭合/投运,0 表示断开。很多网架数据里联络开关的初始状态写作 1 表示“开关在合位”,会因为不同算例习惯不同带坑。
节点编号不一定连续,最好先做一次重映射,让节点集合变成 1 到 N 的连续整数。生成关联矩阵的示例代码如下:
N = max(max(branch(:,1)), max(branch(:,2))); L = size(branch, 1); root = 1; % 变电站母线作为根节点 % 构建关联矩阵 inc:inc(node, k) = 1 表示节点是支路k参考方向的起点,-1为终点 inc = zeros(N, L); for k = 1:L i = branch(k, 1); j = branch(k, 2); if i < j inc(i, k) = 1; inc(j, k) = -1; else inc(j, k) = 1; inc(i, k) = -1; end end这里让参考方向统一按“节点编号从小到大”,是为了后面写节点平衡约束时方便取正负号。关联矩阵当成一个方便的索引工具就行,不需要额外记忆。
3.2 用 YALMIP 写辐射状约束的完整代码
我常用 YALMIP 作为建模层,它的好处是可以先写约束再换求解器。辐射状约束的核心代码段如下:
x = binvar(L, 1); % 1 表示闭合/投运 f = sdpvar(L, 1); % 虚拟商品流,参考方向上的代数量 M = N - 1; % 最大可能虚拟流量,刚好够用 Constraints = []; % 约束1:边数等于节点数减一,保证“树”的数量条件 Constraints = [Constraints, sum(x) == N - 1]; % 约束2:断开支路不允许有虚拟量 Constraints = [Constraints, -M * x <= f <= M * x]; % 约束3:根节点净发出 N-1 个单位虚拟量 rootOut = find(inc(root, :) == 1); rootIn = find(inc(root, :) == -1); Constraints = [Constraints, sum(f(rootOut)) - sum(f(rootIn)) == N - 1]; % 约束4:每个非根节点净消费 1 个单位虚拟量 for v = 1:N if v == root continue; end outIdx = find(inc(v, :) == 1); inIdx = find(inc(v, :) == -1); Constraints = [Constraints, sum(f(outIdx)) - sum(f(inIdx)) == -1]; end % 目标函数示例:最小化闭合线路的电阻总和 obj = sum(branch(:, 3) .* x); % 求解(不指定求解器时 YALMIP 会用默认求解器) optimize(Constraints, obj);目标函数具体是什么不重要,重要的是这组约束把它独自拆出来测试时,应该能稳定输出一棵合法的树。我把这个约束块单独封装成一个函数,任何上游模型需要辐射状拓扑时直接调用它,顺便传一个根节点参数,这样后面叠加潮流约束时不用反复改动。
如果不用 YALMIP,你也可以用 MATLAB 自带的 intlinprog 重写:x 是二进制变量,f 是连续变量,把所有约束写进 Aineq、Aeq 矩阵即可。不过人工拼矩阵比较费眼神,我更建议用小工具箱做建模层,出问题也容易排查。
3.3 结果后处理:用 conncomp 做辐射状复核
求解器输出 x 的数值解后,第一件事不是去看目标函数值,而是验证拓扑是否真的满足辐射状。这一点很多人会偷懒,但我强烈建议写成强制检查项:
selEdgeIdx = find(round(value(x)) > 0.5); G = graph(branch(selEdgeIdx, 1), branch(selEdgeIdx, 2)); nodeBins = conncomp(G); if max(nodeBins) == 1 && numedges(G) == N - 1 disp('拓扑合法:连通且为树'); else error('拓扑非法:存在孤岛或边数不满足'); endMATLAB 的 graph 对象在处理这类检查时非常方便,conncomp 返回每个节点所在的连通分量编号,最大值等于 1 就说明全图只有一个连通分量。配合 numedges 等于 N-1,双重条件同时满足,才叫真正的辐射状。
3.4 一个 5 节点迷你算例:断两条线得到合法树
用一个 5 节点 6 支路的小例子把流程串一遍。支路数据如下(电阻单位省略):
| 支路 | 首端 | 末端 | 电阻 |
|---|---|---|---|
| 1 | 1 | 2 | 0.10 |
| 2 | 2 | 3 | 0.10 |
| 3 | 3 | 4 | 0.10 |
| 4 | 4 | 5 | 0.10 |
| 5 | 5 | 1 | 0.20 |
| 6 | 1 | 3 | 0.30 |
网络本身是一个五边形加一条弦,按照断线解环思想,闭合线路数必须等于 4,也就是要断开 2 条。目标函数取最小化闭合线路的电阻总和,则最优解显然是闭合低阻的 1-2、2-3、3-4、4-5 四条支路,形成一条链状树;断开高阻的 5-1 联络线和 1-3 弦。输出判断里,conncomp 最大值为 1,numedges 为 4,辐射状校验通过。
这个例子虽然小,但足够验证约束写对了没有。我在实际复现时,永远先跑通 5 节点算例,再切换到 33 节点甚至上百节点的系统,因为大算例出了问题很难定位是不收敛还是约束本身写错。
4. 复现时最容易翻车的四个细节(含实测调试记录)
4.1 边数条件满足但出现孤岛:最经典的伪辐射状
前面讲过,只写 sum(x) == N-1 并不保证连通。这里再补一个我实际遇到过的例子:某个 7 节点系统,求解器给出边数为 6,从数字看好得很,但画出来是两个互不相连的链,一个链 4 个节点、3 条边,另一个链 3 个节点、2 条边。两者加起来边数正好 5?不对,应该是 6?这例子不一定准确,但道理是成立的:边数条件只能约束数量,不能约束分布。
调试方法很简单,用 conncomp 看连通分量数。如果大于 1,先检查虚拟商品流约束里有没有把根节点的净发出量写错成 N 而不是 N-1,或者把某个非根节点的净消费量写成了 0。这类低级错误在手工拼约束时很容易发生。
4.2 Big-M 到底取多大才不伤数值
使用单商品虚拟量法时,f 的上界不是随便取的。有人习惯写 M = 1e6,结果模型数值条件很差,整数求解器要么慢,要么干脆给出莫名其妙的结果。这里有一个天然上界:根节点净发出的虚拟量是 N-1,所以任何一条支路上的虚拟流量最大也不可能超过 N-1。直接取 M = N-1 就够,既不会把可行解挡在门外,也把大 M 的数值伤害降到最低。
如果你的网络里某些支路串联了多个节点,单条支路虚拟量上限其实还可以再缩,但没必要,M = N-1 已经足够安全。用代码表示就是M = N - 1;,不要把 M 提升到万级。
4.3 联络开关初始状态的符号陷阱
配电网数据的开关状态经常存在两套习惯:一套用 1 表示闭合,另一套用 0 表示闭合。你从不同的开源算例、论文附录里抄数据时,很容易被这一条暗算。我踩过一次:明明算例里联络开关都写着 1,我理解成“初始闭合”,结果模型天生就是一大片环网,起手就无解。
解决思路是在读数据时统一状态语义,进入模型之前强制转换:
status = branch(:, 5); % 如果原数据用 0 表示闭合,就翻转为 1,同时打印提示同时建议在模型内部只使用 x = 1 表示闭合/投运这唯一语义,所有的“断线”操作都放到结果解释层去做,不要在约束里混用“断开变量”。我见过一些复现代码,用 y = 1 表示断开,目标函数里一会儿用 x 一会儿用 1-y,犯错的概率直接翻倍。
4.4 整数解“看似辐射状”的二次校验方法
即使求解器报告了最优解,也建议做一次独立的拓扑校验,而不是直接信任输出。原因在于第三方求解器在默认参数下可能提前终止,或者因数值问题返回一个次优整数解。我的习惯是:拿到 x 后重新构图,检查三个条件——边数等于 N-1、连通分量数等于 1、基本环集合为空。用 MATLAB 写就是前面那段 conncomp 代码,再加一句:
cyc = cyclebasis(G); assert(isempty(cyc), '存在未断开的环');虽然代数约束理论上已经保证了树,但这一步能把“模型缺陷”和“求解器数值波动”区分开。如果代数解通过了校验,后面加潮流约束时就不用怀疑拓扑问题了,能节省大量排查时间。
5. 故障恢复与网架规划场景下的扩展玩法
5.1 从“最小断线数”到“最大化恢复供电”
前面讲的模型核心是拓扑约束,目标函数随意。实际故障恢复场景中,断线解环思想依然好用,只是把问题的语义换了:故障发生后,故障支路被强制隔离(x = 0 固定),剩余网络通过合上联络开关来恢复失电负荷。这时目标函数变成最大化恢复的负荷量,而辐射状约束照旧,只不过候选支路集合里新增了原本断开的联络开关。
如果你把每个负荷节点的恢复状态设成二进制变量,还要加一条“节点恢复的前提是从根可达”的逻辑约束。这个逻辑约束用虚拟商品流同样能表达:恢复节点才允许消费 1 单位虚拟量,失电节点消费 0。这样辐射状约束和恢复状态就耦合到一起了。
5.2 叠加 DistFlow 潮流约束后的整体模型
纯拓扑模型不需要潮流,但工程上往往还要看电压和负载率。辐测状拓扑真正的优势在于,配电网潮流可以按支路逐级前推回代,优化模型里常用 DistFlow 的线性化形式。忽略网损的一阶线性近似大致长这样:
[ P_v = \sum_{k \in \delta^+(v)} P_k - \sum_{k \in \delta^-(v)} P_k + p_v ]
电压差约束:
[ V_i - V_j = \frac{r_{ij} P_{ij} + x_{ij} Q_{ij}}{V_0} ]
这里的 P、Q 是实际有功无功潮流,和之前虚拟商品流 f 完全无关,两套变量各管各的。很多复现失败的人就是没分清“虚拟流”和“真实潮流”,把两套东西混在一个平衡方程里,导致约束互相打架。
我的建议是先把拓扑约束单独跑通,确认 vaild tree 没问题,再叠加 DistFlow。分步调试能让问题定位清晰得多。
5.3 后续可以怎么改造
- 多时段重构:把决策变量从单断面扩展到时间维度,断线解环约束每个时段都要满足,但联络开关状态变化次数要限制,否则操作成本爆炸。
- 含分布式电源的孤岛运行:如果允许 DG 带孤岛运行,辐射状约束就放宽成“每个连通分量内部是树,且至少包含一个电源”,本质上是多根生成森林问题,根节点集合变化后约束只需要小改。
- 三相不平衡网络:每个开关/线路要按相建模,辐射状约束需要按相别分别检验,复杂度明显上升,但思路一致。
最后再分享一点小经验
我在实际复现那篇 EI 论文的模型时,最大的感受是:论文里的一行“约束(5)满足辐射状拓扑”,到你代码里往往是十几行小心维护的索引操作和符号约定。不要一上来就冲大算例,先把 5 节点模型吃到透,再扩展规模;不要把 M 写太大;不要信任求解器的整数解输出,用 conncomp 复核一次花不了几秒钟。
如果你现在也在做配电网拓扑优化,建议把你的辐射状约束封装成一个独立的函数模块,输入是支路数据和根节点,输出是 x 和 f 以及一套校验标志。这样无论是做故障恢复、网架规划还是后续加灵活性资源,都不用再动这个核心模块。断线解环思想的价值不在于那几条公式,而在于它把“把环网剪成树”这个直觉变成了一个可以交给求解器反复搜索的数学结构。把这个结构吃透了,你再看其他论文里的生成树约束、割集约束,基本都能一眼看穿它到底在约束什么。