简介:一份基于Busacker-Gowan迭代算法的最小费用流MATLAB求解代码,面向网络流理论学习者、算法研究人员,以及物流调度、通信路由等需要成本优化的工程场景。压缩包内仅含一个程序文件,即BGf算法的MATLAB源程序,包体仅1KB,代码量小、结构紧凑,便于直接阅读和二次开发。该程序完整覆盖最小费用流求解流程:从网络数据读取、初始流量与费用初始化,到逐次搜索负增益增广路径并调整流量,最终收敛得到最小费用流与总费用。对照算法步骤可清晰理解迭代终止条件与残差网络更新方式,同时掌握用MATLAB实现图论优化算法的编程思路。代码以函数形式组织,适合课程设计、论文复现或实际项目中的成本优化参考;已有七十二人学习浏览,对希望快速入手经典最小费用流算法的读者是一份轻量且聚焦的实用样例。
1. 基于Busacker-Gowan迭代算法最小费用流代码.zip:这份代码包解决什么问题、适合谁
调拨单压到桌上时,最怕的不是路远,而是“每条路都能走,但加起来总成本不对”:几十个供应点、几百条线路,每条线路有单位运费和运输上限,要凑够总调拨量还让总运费最低。这类问题的标准建模就是最小费用流,而Busacker-Gowan迭代算法是求解它最经典的一条路线:每次在当前网络上找一条单位费用最小的源到汇路径,把能塞的流量全塞进去,直到塞满目标流量F。这份基于Busacker-Gowan迭代算法最小费用流代码.zip,把整套思路封装成解压就能用的代码包,适合刚吃透最大流、正在补费用流短板的工程师和学生,也适合需要一个可信基线去验证自己实现的人——先跑通它,再谈优化。
2. 最小费用流问题与Busacker-Gowan迭代:先把算法立住再动代码
2.1 最小费用流到底在算什么:从容量、费用到目标流量F
先把问题说清楚:最小费用流的输入是五件套——有向图、每条边的容量上限、每条边的单位费用、源点、汇点,再加一个目标流量F。输出是一组边上的流量分配,让从源到汇的总流量恰好等于F,同时总费用(流量乘以单位费用之和)最小。这里有两个容易混的点需要先说破。第一,费用是按单位流量算的线性成本,不是“用一次这条路就付多少钱”的固定成本;第二,目标是“恰好送到F”,不是“尽量多送”。如果需要“流量尽量多且费用尽量小”,那是最小费用最大流,只需要把F设成这个网络的最大流值即可,算法本身不用换。
业务对应关系是这样的:运输调拨里,节点是仓库,边是线路,容量是单次最大调拨量,费用是单位运费;流量计费场景里,节点是路由设备,边是链路,容量是带宽,费用是单位带宽价格。这套模型也出现在生产排产、任务分配、数据中心带宽调度里。带着这些业务看参数表,后面对照代码就不会晕。
2.2 Busacker-Gowan迭代的一句话版本:反复找最短费用路并增广
Busacker-Gowan迭代算法在文献里常被写作“连续最短路增广”,因为它做的事情非常朴素:只要还没送够F,就在当前的残量网络上找一条从s到t的单位费用之和最小的路径,沿着这条路尽可能多地塞流量,塞不动了再重新找。每找一次就是一次迭代。算法主循环可以用下面这段伪代码概括:
while flow < F: dist, prev = shortest_path(residual_graph) # 边权取单位费用 if dist[t] == INF: break # 没有可用路径,已达最大流 d = min(F - flow, 沿路径的最小剩余容量) flow += d total_cost += d * dist[t] 沿路径更新容量,并给每条边补一条反向边(容量=d, 费用=-原费用)这段伪代码就是代码包主函数的骨架。需要重点理解的是最后一行的“反向边加负费用”:残量网络里每条原始边除了剩余容量,还伴生一条反向调整边,容量等于当前已送流量,费用取负。含义是,如果后续迭代发现另一条路径整体更便宜,可以把之前分配的流量退回去重走。整个算法的“反悔”机制就靠这条反向边,也是它能收敛到全局最优的关键。
2.3 反向边与残量网络:这套代码能“反悔”的机制
很多第一次接触的人会在反向边这里栽跟头。原始业务图中只有真实的有向边,残量网络则是算法内部的视角:每条边(u,v)上当前流量为f,那么残量网络上有一条从v到u的反向边,容量为f,费用为-w(u,v)。反向边不是真实线路,不能出现在业务成本里,它只是算法虚构的调整通道。Busacker-Gowan每次迭代都是在包含这些反向边的网络上找最短路,所以它不需要一开始就看穿全局最优,而是靠一轮轮迭代修正之前的分配决策。
调试时我一般会打印残量网络,重点看反向边的费用是否确实取负、容量是否与正向已用流量对称。如果哪条反向边容量对不上,那问题基本出在增广更新代码上。这部分和第5章“总费用对不上手算”的翻车现场直接相关,后面细说。
2.4 SPFA vs Bellman-Ford vs Dijkstra+势能:三种最短路选型
最短路实现是这套算法唯一的性能开关,工程上就三种选择。Bellman-Ford不挑边权,能直接处理负费用边,复杂度O(VE)稳定但慢,适合V很小或对最坏情况有严格要求的场景。SPFA是Bellman-Ford的队列优化,期望速度快、代码量短,但最坏可能退化,是算法竞赛和代码包里最常见的版本。Dijkstra+势能理论上最快,但要求边权非负,而残量网络里天然有负的反向边,所以要用势能变换消掉负边,改造细节在第6章。
代码包里通用的实现一般先给SPFA版,理由很直接:实现最短,而且在有负费用边时不会出错。拿到包后先别急着换算法,先用小数据验证正确性,确认业务结果对得上,再按第6章的改造切到Dijkstra+势能提速。
3. 跑通代码包:从zip解压到第一个算例出数
3.1 拿到包先看什么:用代码解耦的思路判断包质量
拿到zip我会先做一次“盲拆”:解压后不看说明文档,直接按文件大小列目录,先找四样东西——主算法文件、建图/数据加载文件、示例入口文件、测试样例数据。好的代码包一定把“算法”和“业务建模”分开,这就是代码解耦:主算法函数只认节点编号和边表,不认“上海仓”“北京机房”这些业务名,业务名统一在数据加载层映射成编号。如果看到包把业务逻辑和数据解析全拼在同一个函数里,后面换数据时极容易把算法改坏。
阅读顺序我一般固定为:先跑示例入口,再读主算法函数,最后看数据格式。示例入口告诉你“这个函数怎么调用”,主算法告诉你“它到底在算什么”,数据格式告诉你“怎么换成自己的问题”。按这个顺序走,十分钟内基本能判断这套代码能不能直接用于你的场景。
3.2 最小依赖与数据边表格式:把图读进内存
SPFA版的最小费用流只依赖Python标准库的deque,不需要numpy,所以环境准备很轻。我一般建一个干净的虚拟环境再跑:
python -m venv venv source venv/bin/activate # Windows 用 venv\Scripts\activate python --version数据文件推荐用边表,每行四个整数“from to cap cost”,支持#注释。示例文件如下,节点编号从0开始:
# 边表: 起点 终点 容量 单位费用 0 1 10 2 0 2 8 4 1 2 5 1 1 3 7 3 2 3 9 2读文件建图的函数可以写得非常简单,核心就是逐行解析后调add_edge:
def load_graph(path, n): mcf = MinCostFlow(n) with open(path, "r", encoding="utf-8") as fp: for line in fp: line = line.strip() if not line or line.startswith("#"): continue fr, to, cap, cost = map(int, line.split()) mcf.add_edge(fr, to, cap, cost) return mcf这里保持整数解析是有意的:第4章会详细讲,整数容量是这套算法不翻车的前提之一。如果只有浮点容量,也建议先放大成整数跑通,再换算回原来的单位。
3.3 核心实现:Busacker-Gowan迭代的Python代码与逐段说明
这是这类代码包的核心,也是你拿到zip后最该反复读的部分。下面给出一个可直接运行的SPFA版本,变量名保持工程可读性:
# min_cost_flow.py from collections import deque class MinCostFlow: def __init__(self, n): self.n = n self.g = [[] for _ in range(n)] def add_edge(self, fr, to, cap, cost): # 正向边: [终点, 剩余容量, 单位费用, 反向边下标] self.g[fr].append([to, cap, cost, len(self.g[to])]) # 反向边: 初始容量0, 费用取负, 用于回退流量 self.g[to].append([fr, 0, -cost, len(self.g[fr]) - 1]) def min_cost_flow(self, s, t, maxf): n, INF = self.n, 10 ** 18 prev_v, prev_e = [0] * n, [0] * n flow, total = 0, 0 while flow < maxf: # 一次迭代: 在残量网络上跑SPFA找最短费用路 dist = [INF] * n inq = [False] * n dist[s], inq[s] = 0, True q = deque([s]) while q: v = q.popleft() inq[v] = False for i, e in enumerate(self.g[v]): to, cap, w, _ = e if cap > 0 and dist[to] > dist[v] + w: dist[to] = dist[v] + w prev_v[to] = v prev_e[to] = i if not inq[to]: q.append(to) inq[to] = True if dist[t] == INF: break # 残量网络里没有s->t路径,说明已达最大流 # 沿找到的最短路确定增广量: 路径上最小剩余容量与剩余需求的较小值 d = maxf - flow v = t while v != s: e = self.g[prev_v[v]][prev_e[v]] d = min(d, e[1]) v = prev_v[v] # 增广: 累加流量与费用, 更新正向/反向边容量 flow += d total += d * dist[t] v = t while v != s: e = self.g[prev_v[v]][prev_e[v]] e[1] -= d # 正向剩余容量扣减 self.g[v][e[3]][1] += d # 反向容量对称增加 v = prev_v[v] return flow, total这段代码的逻辑说明:SPFA部分只有“剩余容量大于0”的边参与松弛,这保证了找的一定是残量网络里的可行路径。add_edge里每条正向边都配一条费用取负的反向边,二者通过rev下标互指。增广时先扣正向容量再加反向容量,顺序不能反。dist[t]存的是这条路径的单位费用之和,所以费用增量是d乘以dist[t],而不是直接用dist[t]。最后返回的flow如果小于maxf,说明网络最大流不够,调用方必须检查这个信号。
配套的示例入口如下,节点0是源,节点3是汇:
# example.py from min_cost_flow import MinCostFlow if __name__ == "__main__": mcf = MinCostFlow(4) mcf.add_edge(0, 1, 10, 2) mcf.add_edge(0, 2, 8, 4) mcf.add_edge(1, 2, 5, 1) mcf.add_edge(1, 3, 7, 3) mcf.add_edge(2, 3, 9, 2) f, c = mcf.min_cost_flow(0, 3, 6) print("流量:", f, "总费用:", c)运行命令和预期输出:
python example.py # 输出: 流量: 6 总费用: 30手算验证:最短费用路是0→1→3,单位费用2加3等于5,路径容量是min(10,7)=7,目标流量只要6,所以直接送满,总费用6×5=30。把maxf改成10再跑,算法会在送完6后继续找0→2→3路,单位费用4加2等于6,两条路叠加构成最小费用流。这个例子能让读者先相信“输出是对的”,再动手换自己的数据。
3.4 对接自己的数据:四要素建模清单
业务数据进代码前,先回答四个问题:源点编号是多少,汇点编号是多少,目标流量F是多少,边表里每条边的容量和单位费用是否已核对。最常见的翻车是把有向边的方向搞反:有向边只能沿业务方向运输;如果实际是双向运输,要显式加两条方向相反的有向边,不能用一条无向边代替。另一个高频错误是费用填成“整批费用”而不是“单位费用”,单位费用乘以流量才是总费用,填错的话输出在数量级上直接对不上。建模阶段多花十分钟核对这四要素,后面能少踩一半坑。
4. 参数怎么调:流量上限、整数容量与费用符号的三处细节
4.1 目标流量F:先用最大流兜底
拿到一个新网络,我第一件事不是直接跑费用流,而是先把F设成极大值跑一遍,只看返回的flow。这个flow就是当前网络能承受的最大流量。如果业务要求的F大于这个值,网络无可行解,后面费用再小也没意义。代码写起来很直接:
f_upper, _ = mcf.min_cost_flow(s, t, 10 ** 18) print("最大流上限:", f_upper) if demand > f_upper: raise ValueError("目标流量超过网络最大流,先扩容或减量")这里利用了这个包的一个特性:当残量网络里没有s到t路径时,break返回当前flow。它顺手扮演了一次最大流计算器,省得再引一套Dinic实现。注意只看返回的flow,费用忽略。这个上限值后面还能用来算最小费用最大流,一鱼两吃。
4.2 整数容量为什么是后悔药:迭代次数可预测
Busacker-Gowan的迭代次数由“每次增广量至少为1”保证。如果容量和F都是整数,增广量d是沿路径上整数剩余容量的最小值,必然不小于1,所以循环最多F次,最坏情况可预估。这就是代码包全部使用整数数据的核心原因。浮点容量下,d可能小到1e-12,迭代次数会爆炸,表现为程序怎么跑都跑不完。
工程里的常见做法是把浮点容量和费用统一放大成整数:乘以100或1000再取整,跑完把总费用除以放大倍数换算回来。这样既保留精度,又让迭代次数回到可控范围。费用项也同样处理,避免浮点比较误差污染最短路判断。
4.3 负费用边的使用边界
建模里可能存在“越运越赚”的边,例如带处置补贴的线路,单位费用就是负数。SPFA能处理负边,但前提是残量网络里不能有负费用环。一旦存在负环,算法会沿环无限刷低费用,循环停不下来。检测方法也很简单:在SPFA里记录每个节点入队次数,超过节点数n就判定有负环,直接报错。具体代码片段在第5章的翻车现场里给出。遇到负环,先把环上的边打印出来人工核对费用,确认是哪条异常便宜的边在作怪,再决定是去掉它还是加业务约束。
4.4 数据规模与最短路选型
下表是我在不同规模下的选择习惯,供你拿到新算例时对照:
| 规模特征 | 推荐最短路 | 预期迭代次数 | 说明 |
|---|---|---|---|
| 节点数≤200,F≤1000 | SPFA | 不超过F | 代码最短,调试友好 |
| 节点数≤2000,F≤10^5 | 势能Dijkstra | 不超过F | 每次O(E log V),快一个量级 |
| 节点数≥10^4且容量连续 | 先整数化,再势能Dijkstra | 不可控时换费用缩放算法 | 普通逐路增广会慢到不可用 |
当F很大(比如10^7)而网络又稀疏时,逐条增广的迭代次数成为瓶颈。这时不建议硬扛Busacker-Gowan,直接换网络单纯形或费用缩放实现。这也是这套代码的边界:它适合中小流量、路径清晰的问题,不适合每秒上亿规模的调度系统。认清这个边界,才知道什么时候该投入重构。
5. 避坑:Busacker-Gowan迭代算法的5个翻车现场与排查方法
5.1 翻车现场一:程序不结束,像死循环
现象:设好F和算例后程序长时间不退,CPU占用拉满。原因:残量网络里存在负费用环。SPFA在负环上停不下来,最短路距离不断变小。解决:在SPFA的入队逻辑里加入队计数,超过节点数n就抛出“负环”异常。参考片段如下:
cnt = [0] * n ... if not inq[to]: q.append(to) inq[to] = True cnt[to] += 1 if cnt[to] > n: raise RuntimeError("检测到负费用环,请检查费用为负的边")出异常后,把费用为负的边全部打印出来,逐条核对业务含义。绝大多数情况是某条边的补贴费用被重复录入,或者反向边的费用符号写成了正数。
5.2 翻车现场二:总费用对不上手算
现象:同样的图,手算是30,程序输出却是别的值。原因:大概率是反向边更新写错。add_edge里正向边和反向边的rev下标必须互指;增广时更新的是self.g[v][e[3]][1],如果rev下标记错,就会把容量加到不相关的边上。解决:写一个校验函数,检查所有边正向剩余容量与反向剩余容量之和是否恒等于初始容量。如果哪条边不满足,顺着增广回溯代码逐行查,重点看反向边的下标是不是取到了正向边的len值。
5.3 翻车现场三:流量没送满就停了
现象:返回的flow小于F,但业务预期要求必须送满。原因:代码在dist[t]等于INF时break,这意味着残量网络里已没有s到t的可行路径,当前流达到该网络的最大流;如果这个最大流值本来就小于F,说明网络容量不够或建图漏了边。解决:先用min_cost_flow(s, t, 10**18)看最大流上限,再对比F。如果上限足够,检查是不是把F设成了极大值而不是业务需求值;如果上限不够,回去核对边表里哪些边漏了或容量偏小。
5.4 翻车现场四:浮点容量导致结果抖动或卡死
现象:把算例改成浮点容量后,输出带一堆奇怪的尾数,程序偶尔卡住。原因:浮点容量下增广量d不断变小,加上最短路比较用了精确相等,可能在接近最优解时出现轻微抖动或微小负环。解决:容量和费用统一转成整数再跑,算完再除回放大倍数;最短路松弛判断用“严格大于”,即dist[to] > dist[v] + w,不写等号。如果必须保留浮点,则给比较加一个1e-9的容差,但迭代次数仍然不可控,不建议依赖这个方案。
5.5 翻车现场五:把最小费用流跑成了最小费用最大流
现象:F设得很大,程序输出一个巨大的流量,费用确实很低,但业务根本用不了那么多。原因:把“恰好送F”和“尽量多送且费用小”搞混了。这个包实现的是给定F的最小费用流,不是自动追求最大流。解决:调F之前先按第4.1节算出最大流上限,再按业务需求设定F。要跑最小费用最大流,就把F设成上限值,一次跑完即可。如果直接拿10^18当F,费用最优没错,但流量完全脱离业务控制。
6. 把代码用出价值:势能Dijkstra改造、对拍验证与一个工程习惯
6.1 用势能Dijkstra替换SPFA:提速改造只动三处
当图上到几千个节点以后,SPFA每次迭代都可能退化,我一般会切到Dijkstra加势能。改造就三步:先用SPFA跑一遍得到每个点的初始势能h;Dijkstra里边的权用“原费用+h[u]-h[v]”;每轮得到dist后,把所有点的h累加上dist。关键代码差异如下:
# Dijkstra 用势能消掉负边: 缩减费用 cost + h[u] - h[to] 保证非负 if cap > 0: nd = d + cost + h[u] - h[to] if dist[to] > nd: dist[to] = nd prev_v[to] = u prev_e[to] = i # 每轮结束后更新势能 for v in range(n): if dist[v] < INF: h[v] += dist[v]这段改造不改变计算结果,只改变速度。改完之后,增广路径的真实费用可以直接回溯原边权累加,不必依赖h的差值。每次重构后务必用对拍脚本确认改造前后的flow和total完全一致。
6.2 验证方法:用线性规划对拍当标准答案
我不太相信一次写对的算法代码,所以任何实现都会配一个对拍脚本。最小费用流本质是线性规划,网络流矩阵又是幺模矩阵,LP得到的最优解就是整数最优解,可以放心拿它当标准答案。示意代码如下:
# 示意:用 scipy.optimize.linprog 做小规模对拍 # 变量 x_e: 每条边流量; 约束: 流量守恒 + 容量上限 + 源点净流出=F # 目标: min Σ cost_e * x_e from scipy.optimize import linprog res = linprog(c, A_eq=A, b_eq=b, bounds=[(0, cap) for cap in caps]) # 对比 res.fun 与 min_cost_flow 返回的 total对拍时用小图,比如节点数不超过10,随机生成几十组容量和费用,每组都对比flow和cost。如果cost有差异,优先查反向边更新和增广路径回溯这两段代码。这套对拍脚本留着,以后每次改代码都是后悔药,避免把原本跑通的正确性改丢。
6.3 我的一个工程习惯
最后说一个长期养成的习惯:我会把网络的节点名直接写成业务名,在数据加载层映射成编号,而不是在算法文件里写数字0、1、2。出问题时能直接把结果拉回业务表里人工核对。算法代码本身永远不拼业务判断。另外,每次重构后跑一遍对拍脚本,这相当于给自己的实现留一份保险。Busacker-Gowan这类代码,坑大多数不在数学公式,而在边的更新和数据的符号上。把这两处看住,这个方向就值得投入。希望帮到你。
本文还有配套的精品资源,点击获取