简介:这份资源面向交通工程、运筹优化方向的学习者与研究者,聚焦用户均衡(UE)交通分配问题的算法实现。它以经典SiouxFalls路网为算例,采用Frank-Wolfe算法迭代求解,最终得到满足Wardrop第一原理的均衡流量——即任何出行者都无法通过单方面改变路径来降低自身出行费用。路阻函数使用美国联邦公路局BPR函数,出行需求由OD矩阵给定,输出为每条路段上的分配流量,适合用作交通网络分析课程实验或算法复现练习。压缩包共7个文件,以6个txt数据文件和1个py脚本为主,txt涵盖路网拓扑、节点、路段及OD需求等基础数据,py脚本承载核心求解流程,整体约7KB,结构精简便于快速上手。目前已有2331人学习下载。读者可借此掌握Frank-Wolfe在UE模型中的完整实现思路,理解BPR阻抗与流量更新的耦合关系,并可直接运行脚本复现SiouxFalls分配结果,为后续扩展至更大规模网络或对比其他算法提供可复用的代码基础。
1. 从 SiouxFalls 网络说起:为什么 UE 模型值得你花两周去啃
如果你正在做交通工程方向的毕设或课设,大概率绕不开一个名字——SiouxFalls 网络。它只有 24 个节点、76 条路段,放在真实城市路网里小得像个玩具,但几乎所有交通分配算法的论文都拿它当基准测试。原因很简单:它足够小,小到你能在笔记本上跑完整个迭代过程;又足够典型,典型到能暴露算法收敛性的所有问题。
UE(User Equilibrium)用户均衡模型的核心假设只有一句话:每个出行者都选对自己最有利的路径,最终没有人能通过单方面换路来缩短自己的出行时间。这个均衡状态用数学语言写出来就是 Wardrop 第一原理,求解它本质上是一个带线性约束的非线性凸优化问题。Frank-Wolfe 算法就是解这类问题最经典的下降方向法之一,它每一步只需要解一个全有全无分配,计算量小、实现门槛低,特别适合作为交通分配入门的第一个完整实现。
这套方案要解决的就是:给定 SiouxFalls 网络的拓扑结构、路段自由流走行时间和通行能力,用 Python 从零实现 Frank-Wolfe 算法,求出 UE 均衡状态下的路段流量和路径走行时间。适合正在做交通规划课设、交通分配算法毕设,或者想从零理解 UE 求解全流程的人。
2. 需求拆解与技术选型:为什么是 Frank-Wolfe 而不是别的
2.1 UE 模型的数学本质与求解难点
UE 模型的数学形式是一个非线性规划问题:
min Z(x) = Σ_a ∫₀^{x_a} t_a(w) dw s.t. Σ_k f_k^{rs} = q_{rs} ∀ r,s f_k^{rs} ≥ 0 x_a = Σ_{rs} Σ_k f_k^{rs} · δ_{a,k}^{rs}其中 t_a(x_a) 是路段走行时间函数,通常用 BPR 函数表示:
t_a(x_a) = t_a^0 · [1 + α · (x_a / c_a)^β]t_a^0 是自由流时间,c_a 是通行能力,α 和 β 是参数,工程上常取 α=0.15、β=4。目标函数 Z(x) 对 x_a 的偏导数恰好就是 t_a(x_a),这个性质是 Frank-Wolfe 算法能用的关键——它让下降方向的搜索变得极其简单。
求解难点在于:目标函数是凸的,但约束是线性的,变量维度等于路段数乘以 OD 对数量。SiouxFalls 有 76 条路段、528 个 OD 对,直接上通用优化器不是不行,但你会失去对算法每一步的理解。Frank-Wolfe 的价值在于它把每一步都拆得清清楚楚:方向搜索、步长确定、流量更新,三个步骤对应三段代码,调试的时候哪里出问题一眼就能看出来。
2.2 为什么选 Frank-Wolfe 而不是梯度下降或内点法
梯度下降需要投影到可行域上,而 UE 的可行域是高维空间里的一个多面体,投影操作本身就是一个子优化问题,实现复杂度直接翻倍。内点法收敛快,但需要计算 Hessian 矩阵,76 条路段还好,一旦扩展到几百条路段,内存和计算量都会成为瓶颈。
Frank-Wolfe 的优势在于:每一步的方向搜索等价于一次全有全无分配(All-or-Nothing Assignment),也就是把所有 OD 需求分配到当前最短路径上。这个操作只需要跑一次 Dijkstra 或 Bellman-Ford,时间复杂度可控,而且物理意义明确——它模拟的是“所有出行者都按当前路况选最短路径”的极端情况。步长用线搜索确定,一维搜索计算量可以忽略。
注意:Frank-Wolfe 的收敛速度是次线性的,通常需要几百到上千次迭代才能达到较高精度。SiouxFalls 网络上一般跑 200 到 500 次迭代就能得到比较稳定的结果,但如果你要写论文做对比实验,建议跑到相对间隙(Relative Gap)小于 1e-4 再停。
2.3 技术栈与依赖选择
Python 实现这套算法,核心依赖只有三个:
| 库 | 用途 | 是否必须 |
|---|---|---|
| NumPy | 矩阵运算、数组操作 | 必须 |
| NetworkX | 图结构构建、最短路径计算 | 推荐 |
| Matplotlib | 收敛曲线绘制、结果可视化 | 可选 |
NetworkX 的shortest_path函数底层是 Dijkstra,对于 SiouxFalls 这种规模的网络足够快。如果你追求极致性能,可以用scipy.sparse.csgraph.dijkstra,但 NetworkX 的接口更友好,调试阶段更省心。
数据文件方面,SiouxFalls 网络的标准输入包括三个文件:SiouxFalls_node.tntp(节点坐标)、SiouxFalls_net.tntp(路段属性)、SiouxFalls_trips.tntp(OD 需求矩阵)。这些文件在交通工程领域是公开的标准测试数据,格式是 TNTP 格式,用文本编辑器就能打开。
3. 总体架构设计:四个模块与一条数据流
3.1 模块划分与职责边界
整个方案拆成四个模块,每个模块只做一件事:
数据加载模块负责读取 TNTP 格式的三个文件,输出统一的内部数据结构。节点表包含节点 ID 和坐标,路段表包含起终点、自由流时间、通行能力、BPR 参数,OD 表包含起点、终点、需求量。这个模块的关键是处理 TNTP 格式里的特殊符号,比如~和;,以及 OD 矩阵的稀疏存储格式。
网络构建模块把路段表转换成邻接字典或邻接矩阵,同时预计算每条路段的 BPR 参数。邻接结构的选择直接影响最短路径计算的效率,SiouxFalls 只有 24 个节点,用邻接字典就够了,但如果你要扩展到更大的网络,建议用 CSR 格式的稀疏矩阵。
Frank-Wolfe 核心模块是整个方案的心脏,包含三个子过程:方向搜索(全有全无分配)、步长计算(线搜索)、流量更新(凸组合)。这个模块的输入是当前路段流量,输出是更新后的路段流量和当前的相对间隙。
结果输出模块负责把最终流量、走行时间、收敛曲线整理成可读的报告。通常的做法是输出一个 CSV 文件记录每次迭代的目标函数值和相对间隙,再用 Matplotlib 画一张收敛曲线图。
3.2 数据流与迭代逻辑
数据流从 OD 需求出发,经过全有全无分配得到初始流量,然后进入 Frank-Wolfe 主循环。每次迭代做三件事:
- 根据当前路段流量计算走行时间,更新网络权重
- 对每个 OD 对跑最短路径,把需求分配到最短路径上,得到辅助流量
- 计算下降方向,做线搜索确定步长,更新路段流量
循环的终止条件通常有两个:迭代次数达到上限,或者相对间隙小于阈值。相对间隙的定义是:
gap = Σ_a x_a · (t_a(x_a) - t_a^*(x_a)) / Σ_a x_a · t_a(x_a)其中 t_a^*(x_a) 是辅助流量对应的走行时间。这个指标衡量的是当前流量离均衡状态还有多远,值越小越接近均衡。
3.3 关键数据结构设计
路段用字典存储,键是路段 ID,值是一个包含起终点、自由流时间、通行能力、当前流量的对象。OD 需求用嵌套字典存储,外层键是起点,内层键是终点,值是需求量。
# 路段数据结构示例 links = { 1: { 'from': 1, 'to': 2, 'free_flow_time': 6.0, 'capacity': 25900.0, 'bpr_alpha': 0.15, 'bpr_beta': 4.0, 'flow': 0.0 }, # ... 其余 75 条路段 } # OD 需求数据结构示例 od_demand = { 1: {2: 100.0, 3: 200.0, ...}, 2: {1: 150.0, 3: 300.0, ...}, # ... 其余起点 }这种结构的优点是直观,调试的时候直接打印就能看懂。缺点是内存占用比 NumPy 数组大,但 SiouxFalls 规模小,这点开销可以忽略。
4. 核心模块实现:从全有全无分配到线搜索
4.1 数据加载与网络构建
TNTP 格式的文件用分号分隔字段,用~表示注释。读取的时候需要跳过注释行,把每行拆成字段列表。路段文件的关键字段是起点、终点、自由流时间、通行能力、BPR 参数,OD 文件的关键字段是起点、终点、需求量。
import numpy as np import networkx as nx def load_network(net_file): """读取 TNTP 格式的路段文件,返回路段字典和 NetworkX 图""" links = {} G = nx.DiGraph() with open(net_file, 'r') as f: for line in f: line = line.strip() if not line or line.startswith('~'): continue parts = line.split() if len(parts) < 6: continue link_id = int(parts[0]) from_node = int(parts[1]) to_node = int(parts[2]) free_flow = float(parts[3]) capacity = float(parts[4]) # BPR 参数通常取默认值,部分文件会显式给出 alpha = float(parts[5]) if len(parts) > 5 else 0.15 beta = float(parts[6]) if len(parts) > 6 else 4.0 links[link_id] = { 'from': from_node, 'to': to_node, 'free_flow_time': free_flow, 'capacity': capacity, 'bpr_alpha': alpha, 'bpr_beta': beta, 'flow': 0.0 } G.add_edge(from_node, to_node, weight=free_flow, link_id=link_id) return links, G这段代码的关键点在于:G.add_edge的时候把自由流时间作为初始权重,后续迭代中会动态更新。link_id作为边的属性存储,方便从最短路径结果反查路段流量。
OD 文件的读取稍微麻烦一点,TNTP 格式的 OD 矩阵通常写成嵌套结构,外层是起点,内层是终点和需求量。解析的时候需要维护当前起点,遇到新的起点就更新状态。
def load_od_demand(od_file): """读取 TNTP 格式的 OD 需求文件,返回嵌套字典""" od = {} current_origin = None with open(od_file, 'r') as f: for line in f: line = line.strip() if not line or line.startswith('~'): continue if line.startswith('Origin'): current_origin = int(line.split()[1]) od[current_origin] = {} elif current_origin is not None: parts = line.split(';') for part in parts: part = part.strip() if not part: continue dest, demand = part.split(':') od[current_origin][int(dest)] = float(demand) return od4.2 BPR 函数与全有全无分配
BPR 函数是路段走行时间的核心计算逻辑,输入是当前流量,输出是走行时间。实现的时候注意参数的单位和量纲,自由流时间通常是分钟,通行能力是辆/小时,流量也是辆/小时。
def bpr_time(link, flow): """计算 BPR 走行时间""" t0 = link['free_flow_time'] c = link['capacity'] alpha = link['bpr_alpha'] beta = link['bpr_beta'] return t0 * (1 + alpha * (flow / c) ** beta)全有全无分配是 Frank-Wolfe 方向搜索的核心。对每个 OD 对,用当前走行时间作为权重跑一次最短路径,然后把该 OD 对的需求全部分配到这条路径上。所有 OD 对处理完后,累加得到辅助流量。
def all_or_nothing(links, G, od_demand): """全有全无分配,返回辅助流量字典""" # 更新图权重为当前走行时间 for link_id, link in links.items(): u, v = link['from'], link['to'] G[u][v]['weight'] = bpr_time(link, link['flow']) aux_flow = {link_id: 0.0 for link_id in links} for origin, destinations in od_demand.items(): for dest, demand in destinations.items(): if demand <= 0: continue try: path = nx.shortest_path(G, origin, dest, weight='weight') except nx.NetworkXNoPath: continue # 把需求累加到路径经过的每条路段上 for i in range(len(path) - 1): u, v = path[i], path[i + 1] link_id = G[u][v]['link_id'] aux_flow[link_id] += demand return aux_flow这段代码的效率瓶颈在nx.shortest_path的调用次数。SiouxFalls 有 528 个 OD 对,每次迭代要跑 528 次 Dijkstra,200 次迭代就是 10 万次以上。实际跑起来大概需要几分钟,可以接受。如果要加速,可以把所有 OD 对的最短路径计算并行化,或者用scipy.sparse.csgraph.dijkstra一次性算出所有起点到所有终点的最短路径树。
4.3 线搜索与流量更新
线搜索的目标是找到最优步长 λ,使得目标函数 Z(x + λ(d - x)) 最小。由于 Z 是凸函数,这个一维搜索可以用黄金分割法或三分法求解。工程上更常用的做法是直接求导令其为零,但 BPR 函数是非线性的,解析解不好求,所以数值搜索更实际。
def line_search(links, aux_flow, direction, tol=1e-6): """黄金分割法线搜索最优步长""" def objective(lam): total = 0.0 for link_id, link in links.items(): x = link['flow'] + lam * direction[link_id] t0 = link['free_flow_time'] c = link['capacity'] alpha = link['bpr_alpha'] beta = link['bpr_beta'] # 目标函数的积分形式 total += t0 * (x + alpha * c / (beta + 1) * (x / c) ** (beta + 1)) return total a, b = 0.0, 1.0 gr = (np.sqrt(5) - 1) / 2 c = b - gr * (b - a) d = a + gr * (b - a) while abs(b - a) > tol: if objective(c) < objective(d): b = d d = c c = b - gr * (b - a) else: a = c c = d d = a + gr * (b - a) return (a + b) / 2黄金分割法的收敛速度是线性的,对于一维搜索来说足够快。tol控制搜索精度,1e-6 对应步长精度大约百万分之一,实际跑的时候可以放宽到 1e-4 节省时间。
流量更新就是凸组合:
def update_flow(links, aux_flow, step_size): """按步长更新路段流量""" for link_id in links: links[link_id]['flow'] = ( links[link_id]['flow'] + step_size * (aux_flow[link_id] - links[link_id]['flow']) )4.4 主循环与收敛判断
主循环把上面几个模块串起来,每次迭代计算相对间隙,判断是否满足终止条件。
def frank_wolfe(links, G, od_demand, max_iter=500, gap_tol=1e-4): """Frank-Wolfe 主循环""" # 初始全有全无分配 aux_flow = all_or_nothing(links, G, od_demand) for link_id in links: links[link_id]['flow'] = aux_flow[link_id] history = [] for iteration in range(max_iter): # 方向搜索 aux_flow = all_or_nothing(links, G, od_demand) # 计算相对间隙 numerator = 0.0 denominator = 0.0 for link_id, link in links.items(): t_current = bpr_time(link, link['flow']) t_aux = bpr_time(link, aux_flow[link_id]) numerator += link['flow'] * (t_current - t_aux) denominator += link['flow'] * t_current gap = numerator / denominator if denominator > 0 else 0.0 history.append({'iteration': iteration, 'gap': gap}) if gap < gap_tol: break # 线搜索与更新 direction = {lid: aux_flow[lid] - links[lid]['flow'] for lid in links} step = line_search(links, aux_flow, direction) update_flow(links, aux_flow, step) return links, history相对间隙的计算逻辑是:分子是当前流量与辅助流量之间的走行时间差加权和,分母是当前总走行时间。这个指标无量纲,值越小说明越接近均衡。通常跑到 1e-4 以下就可以认为收敛了,论文里常报告 1e-6 的结果。
5. 避坑与常见问题:那些让你跑不出结果的地方
5.1 流量初始化为零导致除零错误
现象:第一次调用all_or_nothing时,BPR 函数里flow / capacity没问题,但相对间隙计算时分母Σ x_a · t_a可能为零,导致ZeroDivisionError。
原因:初始流量全为零,所有路段的走行时间等于自由流时间,但总走行时间为零。
解决:在进入主循环之前先做一次全有全无分配,把辅助流量赋给初始流量。这样初始流量就不为零了,相对间隙计算也不会出问题。
5.2 最短路径权重更新不及时
现象:迭代过程中流量在变,但最短路径结果始终不变,收敛曲线一条直线。
原因:all_or_nothing里更新了G[u][v]['weight'],但 NetworkX 的shortest_path可能缓存了之前的结果,或者更新权重的时候漏掉了某些边。
解决:每次调用all_or_nothing之前强制刷新所有边的权重,确保G[u][v]['weight']等于当前流量下的 BPR 时间。可以在更新完权重后打印几条边的权重确认一下。
5.3 线搜索步长始终接近零
现象:迭代几百次后相对间隙下降极慢,步长越来越小,流量几乎不变。
原因:Frank-Wolfe 在接近最优解时会出现“锯齿”现象,下降方向与梯度方向接近正交,步长自然趋近于零。这是算法本身的特性,不是实现错误。
解决:如果精度要求不高,跑到 1e-3 或 1e-4 就可以停。如果要求高精度,可以考虑用共轭方向法或松弛法加速收敛。另外检查一下 BPR 参数是否合理,α=0.15、β=4 是标准值,如果文件里给的是其他值,确认单位是否一致。
5.4 OD 需求矩阵解析错误
现象:程序跑起来但总流量对不上,或者某些 OD 对的需求为零。
原因:TNTP 格式的 OD 文件里,需求值可能用科学计数法表示,或者分隔符不是标准的分号。另外有些文件会在行末加注释,解析时需要跳过。
解决:读取 OD 文件后打印总需求量和 OD 对数量,和已知的 SiouxFalls 标准值对比。SiouxFalls 的总需求量是 360,600,OD 对数量是 528。如果对不上,逐行检查解析逻辑。
5.5 收敛曲线震荡不下降
现象:相对间隙在某个值附近上下波动,不持续下降。
原因:可能是线搜索精度不够,步长没取到最优值;也可能是 BPR 函数参数写反了,比如把 α 和 β 搞混。
解决:先检查 BPR 函数的实现,确认alpha * (flow / capacity) ** beta的顺序正确。然后提高线搜索精度,把tol从 1e-4 降到 1e-6。如果还是震荡,打印每次迭代的步长和目标函数值,看看是不是步长在 0 和 1 之间跳变。
6. 验收与答辩:怎么证明你的 UE 模型真的收敛了
6.1 收敛性验证的三条硬指标
答辩的时候导师不会只看你跑出了结果,他会问你:你怎么知道这就是均衡解?通常的做法是从三个维度验证:
相对间隙曲线:把每次迭代的相对间隙画成折线图,横轴是迭代次数,纵轴是间隙值(对数坐标)。合格的曲线应该单调下降,最终稳定在 1e-4 以下。如果曲线震荡或者下降太慢,说明实现有问题。
流量守恒检验:对所有路段,流入某个节点的总流量应该等于流出该节点的总流量(加上该节点的产生/吸引量)。这个检验能发现分配过程中的累加错误。
路径走行时间一致性:在均衡状态下,同一个 OD 对之间所有被使用的路径走行时间应该相等,且小于等于未被使用路径的走行时间。随机抽几个 OD 对,用最终流量算一下各条路径的走行时间,验证这个性质。
6.2 答辩时怎么讲清楚 Frank-Wolfe 的每一步
答辩陈述建议按这个顺序讲:先讲 UE 模型的数学形式,重点说清楚目标函数为什么是那个积分形式;然后讲 Frank-Wolfe 的三个步骤,每一步对应代码里的哪个函数;最后展示收敛曲线和最终流量分布。
导师最常问的问题是:“为什么用 Frank-Wolfe 而不是直接调 scipy.optimize?”回答的思路是:Frank-Wolfe 的每一步都有明确的物理意义,方向搜索对应全有全无分配,步长对应出行者对路况变化的响应程度。而且它不需要计算 Hessian 矩阵,扩展到大规模网络时内存友好。
另一个高频问题是:“你怎么保证收敛到全局最优?”回答:UE 模型的目标函数是严格凸的,约束是线性的,所以局部最优就是全局最优。Frank-Wolfe 每次迭代都保证目标函数不增,只要步长选得合适,一定收敛。
6.3 结果展示的实用技巧
最终结果建议用两张图展示:一张是收敛曲线,一张是路段流量分布图。流量分布图可以用 NetworkX 的draw函数,边的粗细表示流量大小,颜色表示走行时间。SiouxFalls 网络节点少,画出来很直观。
如果要做对比实验,可以跑不同 BPR 参数下的结果,或者对比 Frank-Wolfe 和梯度下降的收敛速度。表格里列出迭代次数、最终间隙、总走行时间三个指标,一目了然。
提示:答辩前把代码跑一遍,确认没有随机性。Frank-Wolfe 本身是确定性的,但如果你用了随机初始流量或者随机 OD 抽样,结果每次都不一样,答辩现场会很尴尬。
从那以后我每次做交通分配实验,都会在跑完主循环后强制走一遍流量守恒检验和路径走行时间一致性检验,这两个检查花不了五分钟,但能挡住答辩现场百分之八十的追问。希望帮到你。
本文还有配套的精品资源,点击获取