基于Frank-Wolfe算法的SiouxFalls网络用户均衡交通分配实战
简介这份资源面向交通工程、运筹优化方向的学习者与研究者聚焦用户均衡UE交通分配问题的算法实现。它以经典SiouxFalls网络为算例采用Frank-Wolfe算法迭代求解最终得到满足Wardrop第一原理的均衡路段流量即任何出行者都无法通过单方面改变路径来降低自身出行费用。路阻函数使用美国联邦公路局BPR函数出行需求由OD矩阵给定适合作为交通网络分析与算法编程的入门实践素材。压缩包共7个文件以6个txt数据文件和1个py脚本为主txt分别承载网络拓扑、节点、路段与OD需求等基础数据py文件为完整的Frank-Wolfe求解程序整体约7KB结构精简、便于直接运行与二次修改。目前已有2331人学习下载。读者可据此掌握UE模型的建模思路、BPR阻抗计算、方向搜索与流量更新的完整流程并对照输出结果验证算法收敛性为后续研究或课程作业提供可复用的代码框架。1. 从 SiouxFalls 说起为什么交通分配值得你花一个周末跑通如果你手头有一份路网数据和 OD 需求表想知道每个路段上到底会分到多少车流那 UE用户均衡模型就是绕不开的经典方法。SiouxFalls 网络是交通分配领域最常用的测试网络之一24 个节点、76 条有向边、OD 对数量适中跑起来不费劲但足够验证算法是否正确。这个项目做的事情很明确用 Frank-Wolfe 算法求解 UE 用户均衡模型输入 SiouxFalls 网络拓扑和 OD 需求输出每条路段上的均衡流量和旅行时间。适合交通工程方向的学生、做交通仿真的开发者以及想理解「Wardrop 第一原理怎么变成代码」的工程师。下面我从目录结构开始一路拆到核心迭代逻辑和二次开发接口。2. 项目结构与技术选型从数据文件到求解器的分层设计2.1 目录结构数据、算法、可视化三层分离拿到源码包后我习惯先看目录树。这个项目的组织方式很清晰常见做法是按功能分层ue-frank-wolfe/ ├── data/ │ ├── SiouxFalls_network.csv # 路段拓扑起终点、自由流时间、通行能力 │ └── SiouxFalls_od.csv # OD 需求矩阵 ├── src/ │ ├── network.py # 路网加载与邻接表构建 │ ├── assignment.py # Frank-Wolfe 主循环 │ ├── shortest_path.py # Dijkstra 最短路 │ └── utils.py # BPR 函数、收敛指标 ├── main.py # 入口脚本 ├── requirements.txt └── README.mddata/放的是 SiouxFalls 标准测试数据CSV 格式方便替换成自己的路网。src/里network.py负责把边表转成邻接结构shortest_path.py实现最短路搜索assignment.py是 Frank-Wolfe 的核心迭代utils.py放 BPR 阻抗函数和相对间隙计算。main.py串起整个流程。这种分层的好处是换网络只改data/换算法只改assignment.py互不干扰。2.2 技术选型为什么用纯 Python NumPy 而不是调库项目没有依赖任何交通仿真专用库核心计算用 NumPy 做向量化。原因很实际Frank-Wolfe 的每一步都需要全有全无分配all-or-nothing也就是把所有 OD 需求一次性加载到最短路上NumPy 的数组操作比纯循环快一个量级。requirements.txt里通常只有numpy、pandas、matplotlib三个包。pandas 读 CSVnumpy 做矩阵运算matplotlib 画收敛曲线和路段流量热力图。如果你要扩展到更大网络可以把最短路换成scipy.sparse.csgraph.dijkstra但 SiouxFalls 这个规模用纯 Python 堆实现完全够用。2.3 数据格式边表和 OD 矩阵的字段含义SiouxFalls_network.csv的列一般是from_node, to_node, free_flow_time, capacity, length。其中free_flow_time是零流量下的行驶时间capacity是路段通行能力BPR 函数会用到这两个参数。SiouxFalls_od.csv是 OD 矩阵行和列都是节点编号值是对应的出行需求。加载时要注意节点编号从 1 开始还是从 0 开始SiouxFalls 标准数据通常从 1 开始代码里做索引映射时要减 1否则最短路会找不到节点。这个坑我在第一次跑的时候踩过后面会细说。3. Frank-Wolfe 核心迭代从全有全无分配到收敛判定3.1 BPR 阻抗函数路段走行时间怎么随流量变化UE 模型的核心是路段阻抗函数项目里用的是 BPR 函数import numpy as np def bpr_time(flow, free_flow_time, capacity, alpha0.15, beta4.0): 计算路段在给定流量下的走行时间 flow: 路段流量 (veh/h) free_flow_time: 自由流时间 (min) capacity: 通行能力 (veh/h) alpha, beta: BPR 参数默认 0.15 和 4.0 ratio flow / capacity return free_flow_time * (1 alpha * np.power(ratio, beta))这个函数的逻辑很直接流量越大走行时间越长。alpha控制拥堵增加的幅度beta控制曲线的陡峭程度。标准 BPR 用 0.15 和 4.0但你可以根据实际路网标定。注意capacity不能为 0否则会除零。如果某条边容量缺失我一般会填一个默认值比如 1800 veh/h或者直接报错让用户补数据。3.2 全有全无分配把 OD 需求一次性压到最短路上Frank-Wolfe 每次迭代的第一步是找当前阻抗下的最短路然后把所有 OD 需求分配到最短路上import heapq def dijkstra(adj, source, impedance): 单源最短路返回前驱节点和距离 adj: 邻接表 {node: [(neighbor, edge_id), ...]} source: 起点 impedance: 每条边的当前阻抗数组 dist {node: float(inf) for node in adj} prev {node: None for node in adj} dist[source] 0 pq [(0, source)] while pq: d, u heapq.heappop(pq) if d dist[u]: continue for v, eid in adj[u]: nd d impedance[eid] if nd dist[v]: dist[v] nd prev[v] (u, eid) heapq.heappush(pq, (nd, v)) return dist, prev def all_or_nothing(network, od_matrix, impedance): 全有全无分配每个 OD 对走最短路流量累加到路段上 返回辅助流量向量 aux_flow np.zeros(len(network.edges)) for origin in od_matrix.index: dist, prev dijkstra(network.adj, origin, impedance) for dest in od_matrix.columns: demand od_matrix.loc[origin, dest] if demand 0: continue # 回溯路径累加流量 node dest while prev[node] is not None: u, eid prev[node] aux_flow[eid] demand node u return aux_flowdijkstra用堆实现返回每个节点的最短距离和前驱边。all_or_nothing遍历所有 OD 对把需求累加到最短路径经过的边上。这里有个性能点如果 OD 对很多可以只对需求大于 0 的对做回溯。SiouxFalls 的 OD 矩阵比较稠密但实际路网往往稀疏加个判断能省不少时间。3.3 线搜索与流量更新步长怎么选才不会震荡Frank-Wolfe 的第二步是找最优步长让目标函数在辅助流量和当前流量之间最小化。项目里用的是二分法或者黄金分割法def line_search(current_flow, aux_flow, network, tol1e-6, max_iter50): 在 current_flow 和 aux_flow 之间找最优步长 lambda 目标函数sum over edges of integral of bpr_time def objective(lam): flow current_flow lam * (aux_flow - current_flow) times bpr_time(flow, network.free_flow, network.capacity) # 数值积分近似用梯形法 return np.sum(times * flow) lo, hi 0.0, 1.0 for _ in range(max_iter): mid1 lo (hi - lo) / 3 mid2 hi - (hi - lo) / 3 if objective(mid1) objective(mid2): hi mid2 else: lo mid1 if hi - lo tol: break return (lo hi) / 2这里的目标函数是 Beckmann 变换的积分形式实际代码里可以用解析积分代替数值积分精度更高。线搜索的步长lambda在 0 到 1 之间更新公式是new_flow current_flow lambda * (aux_flow - current_flow)。如果步长选得不好流量会在两个极端之间震荡收敛很慢。黄金分割法比固定步长稳但每轮多几次函数评估SiouxFalls 规模下完全可接受。3.4 收敛判定相对间隙降到多少算收敛Frank-Wolfe 的收敛指标是相对间隙relative gapdef relative_gap(current_flow, aux_flow, impedance): 相对间隙 (当前总阻抗 - 辅助总阻抗) / 当前总阻抗 current_total np.sum(current_flow * impedance) aux_total np.sum(aux_flow * impedance) if current_total 0: return float(inf) return (current_total - aux_total) / current_total相对间隙衡量的是当前流量离均衡状态还有多远。一般降到 1e-4 以下就可以停SiouxFalls 通常迭代 20 到 40 次就能达到。如果迭代 100 次还没收敛大概率是 BPR 参数或者 OD 需求有问题比如需求超过了网络总容量导致阻抗爆炸。4. 跑通与验证从命令行到收敛曲线4.1 环境准备与依赖安装项目依赖很少Python 3.8 以上就行python -m venv venv source venv/bin/activate # Windows 用 venv\Scripts\activate pip install -r requirements.txtrequirements.txt里一般写numpy1.21、pandas1.3、matplotlib3.4。如果装不上先升级 pip。我习惯用虚拟环境避免和系统包冲突。4.2 运行入口main.py 的参数与输出main.py通常接受几个命令行参数python main.py --network data/SiouxFalls_network.csv \ --od data/SiouxFalls_od.csv \ --max-iter 100 \ --gap 1e-4 \ --output results/--network和--od指定数据文件--max-iter是最大迭代次数--gap是收敛阈值--output是结果目录。运行后会生成link_flows.csv路段流量、convergence.csv每轮相对间隙和convergence.png收敛曲线。如果没加--output默认输出到当前目录。4.3 结果验证流量守恒与收敛曲线怎么看跑完之后先检查流量守恒所有从起点出发的流量之和应该等于 OD 需求总和。然后看收敛曲线正常的曲线应该单调下降最后趋于平缓。如果曲线震荡或者上升说明线搜索有问题。路段流量可以画成热力图SiouxFalls 的流量分布应该集中在几条主干道上和实际路网的拥堵模式一致。我一般会对比自由流时间和均衡时间的比值如果某条边比值超过 2说明拥堵严重需要检查容量设置。5. 避坑与排查那些让我熬夜的报错和异常5.1 节点编号从 1 开始导致 KeyError现象运行dijkstra时报KeyError: 0或者找不到节点。原因SiouxFalls 数据节点编号从 1 开始但代码里邻接表用 0 基索引。解决加载数据时统一减 1或者用字典映射。我一般会在network.py里加一行node_map {n: i for i, n in enumerate(sorted(nodes))}后面全部用映射后的索引。5.2 OD 需求超过网络容量导致阻抗爆炸现象迭代几次后流量变成inf或者nan。原因某条边的容量设得太小而 OD 需求又很大BPR 函数里flow/capacity超过 1 很多beta4时时间会爆炸。解决检查容量字段确保单位一致veh/h必要时调大容量或者降低需求。也可以给 BPR 加一个上限比如min(ratio, 10)防止数值溢出。5.3 线搜索步长震荡不收敛现象相对间隙在 1e-2 附近来回跳降不下去。原因线搜索精度不够或者目标函数计算有误。解决把黄金分割的max_iter调到 100tol调到 1e-8。如果还不行检查目标函数的积分是否用了正确的 BPR 原函数。我一般会打印每轮的步长和间隙看是不是步长在 0 和 1 之间乱跳。5.4 最短路回溯时路径成环现象回溯路径时陷入死循环。原因Dijkstra 的前驱更新有误或者图里有负权边。解决确保所有边的阻抗为正BPR 函数在流量非负时不会产生负值。如果用了自定义阻抗检查有没有负数。另外回溯时加一个访问标记防止环。5.5 结果流量不守恒现象路段流量之和和 OD 需求对不上。原因全有全无分配时漏了某些 OD 对或者需求矩阵里有负数。解决检查 OD 矩阵是否所有值非负遍历时用if demand 0: continue跳过零需求。另外确认每个 OD 对都做了回溯没有提前退出。6. 二次开发换网络、换算法、加可视化6.1 替换成自己的路网数据换网络是最常见的需求。你只需要准备两个 CSV边表和 OD 矩阵。边表至少包含from_node, to_node, free_flow_time, capacityOD 矩阵行和列是节点编号。把文件放到data/下改main.py的参数即可。如果节点编号不是连续的network.py里的映射逻辑会自动处理。我一般会先用小网络比如 4 个节点测试确认流程通了再换大网络。6.2 把 Frank-Wolfe 换成其他算法Frank-Wolfe 的收敛速度是次线性的大网络可能需要很多轮。如果你想换算法assignment.py里的主循环是唯一需要改的地方。比如换成梯度投影法只需要把线搜索和流量更新替换成投影步骤。或者用 Origin-Based 算法收敛更快但实现复杂。我的建议是先把 Frank-Wolfe 跑通理解每一步的输入输出再动手改。接口上只要保持solve(network, od_matrix)的签名不变上层代码不用动。6.3 加一个路段流量热力图项目自带收敛曲线但路段流量可视化更直观。用 matplotlib 画热力图import matplotlib.pyplot as plt import numpy as np def plot_link_flows(network, flows, output_path): 画路段流量热力图边宽表示流量大小 fig, ax plt.subplots(figsize(10, 8)) for i, (u, v) in enumerate(network.edges): x [network.nodes[u][x], network.nodes[v][x]] y [network.nodes[u][y], network.nodes[v][y]] width flows[i] / max(flows) * 5 ax.plot(x, y, linewidthwidth, colorsteelblue, alpha0.7) ax.set_title(Link Flows (UE)) plt.savefig(output_path, dpi150) plt.close()这段代码假设节点有坐标SiouxFalls 数据里通常有x, y字段。如果没有可以用 networkx 的 spring layout 自动布局。画出来的图能直观看到哪些路段流量大方便和实际观察对比。6.4 一个让我长记性的教训第一次跑这个项目时我没检查 OD 矩阵的索引直接用了默认的 0 基结果所有需求都加载到了错误的节点上流量图看起来完全乱套。从那以后我每次加载数据都强制走一遍assert set(od.index) set(network.nodes)确认节点集合一致再往下跑。这个习惯帮我省了很多调试时间。希望这个项目能帮你少走弯路顺利跑通自己的第一个 UE 分配。本文还有配套的精品资源点击获取