简介这份资源面向交通工程、运筹优化方向的学习者与研究者聚焦用户均衡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 条路段放在真实城市路网里小得像个玩具但几乎所有交通分配算法的论文都拿它当基准测试。原因很简单它足够小小到你能在笔记本上跑完整个迭代过程又足够典型典型到能暴露算法收敛性的所有问题。UEUser 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.tntpOD 需求矩阵。这些文件在交通工程领域是公开的标准测试数据格式是 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, weightfree_flow, link_idlink_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, weightweight) 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 次 Dijkstra200 次迭代就是 10 万次以上。实际跑起来大概需要几分钟可以接受。如果要加速可以把所有 OD 对的最短路径计算并行化或者用scipy.sparse.csgraph.dijkstra一次性算出所有起点到所有终点的最短路径树。4.3 线搜索与流量更新线搜索的目标是找到最优步长 λ使得目标函数 Z(x λ(d - x)) 最小。由于 Z 是凸函数这个一维搜索可以用黄金分割法或三分法求解。工程上更常用的做法是直接求导令其为零但 BPR 函数是非线性的解析解不好求所以数值搜索更实际。def line_search(links, aux_flow, direction, tol1e-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_iter500, gap_tol1e-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,600OD 对数量是 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 抽样结果每次都不一样答辩现场会很尴尬。从那以后我每次做交通分配实验都会在跑完主循环后强制走一遍流量守恒检验和路径走行时间一致性检验这两个检查花不了五分钟但能挡住答辩现场百分之八十的追问。希望帮到你。本文还有配套的精品资源点击获取