简介：这份资源面向交通工程、运筹优化方向的学习者与研究者，聚焦用户均衡（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 读 CSV，numpy 做矩阵运算，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, alpha=0.15, beta=4.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, tol=1e-6, max_iter=50): """ 在 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 gap）：
def 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里一般写numpy>=1.21、pandas>=1.3、matplotlib>=3.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 很多，beta=4时时间会爆炸。解决：检查容量字段，确保单位一致（veh/h），必要时调大容量或者降低需求。也可以给 BPR 加一个上限，比如min(ratio, 10)，防止数值溢出。
5.3 线搜索步长震荡不收敛
现象：相对间隙在 1e-2 附近来回跳，降不下去。原因：线搜索精度不够，或者目标函数计算有误。解决：把黄金分割的max_iter调到 100，tol调到 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, capacity，OD 矩阵行和列是节点编号。把文件放到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, linewidth=width, color='steelblue', alpha=0.7) ax.set_title('Link Flows (UE)') plt.savefig(output_path, dpi=150) plt.close()这段代码假设节点有坐标，SiouxFalls 数据里通常有x, y字段。如果没有，可以用 networkx 的 spring layout 自动布局。画出来的图能直观看到哪些路段流量大，方便和实际观察对比。
6.4 一个让我长记性的教训
第一次跑这个项目时，我没检查 OD 矩阵的索引，直接用了默认的 0 基，结果所有需求都加载到了错误的节点上，流量图看起来完全乱套。从那以后我每次加载数据都强制走一遍assert set(od.index) == set(network.nodes)，确认节点集合一致再往下跑。这个习惯帮我省了很多调试时间。希望这个项目能帮你少走弯路，顺利跑通自己的第一个 UE 分配。
本文还有配套的精品资源，点击获取
