简介本资源是一份面向交通工程、运筹优化与智能交通系统方向的Python实践代码包聚焦用户均衡UE问题的算法实现与小规模网络验证适合交通规划初学者、研究生课程设计及算法爱好者动手复现。资源采用经典Frank-Wolfe算法求解SiouxFalls标准测试网络的UE分配结果路阻函数严格遵循美国联邦公路局BPR形式OD需求与路网拓扑均以结构化文本文件提供最终输出各路段流量分布完整覆盖建模—求解—结果解析闭环。压缩包共7个文件6个txt数据文件1个核心Python脚本总大小仅7KB轻量易读txt文件分别定义节点、链接、OD对及出行需求py文件封装了Frank-Wolfe迭代逻辑与收敛判断代码简洁、注释清晰便于理解算法原理与调试修改。目前已有2331人学习下载是掌握交通分配基础算法、开展小规模网络仿真实验的高性价比入门材料。1. Frank-Wolfe 真不是“调包就跑”的黑匣子SiouxFalls 网络上跑通 UE 交通分配才是验证你理解交通均衡本质的第一块试金石很多人一看到“用户均衡UE”就下意识点开某深度学习框架的教程结果在交通流分配场景里栽得莫名其妙——UE 不是分类任务没有标签也没有梯度回传它是一组非线性方程组在路径选择理性假设下的不动点解。Frank-WolfeFW算法之所以被交通工程领域沿用近六十年正因为它不依赖导数、天然适配网络拓扑结构、每步迭代只需求解一个最短路子问题且收敛性有严格保证。这份 Python 实现就是把这套经典理论从教科书里“抠”出来在 SiouxFalls 这个被引超 2000 次的基准网络上完整复现从 OD 矩阵读入、自由流时间加载、BPR 函数建模、路径生成、流量分配到 FW 主循环的步长控制与收敛判定——全部可调试、可打断、可替换模块。适合刚学完《交通分配理论》想动手验证公式的研究生也适合需要嵌入自定义路网逻辑的交通仿真工程师。它不封装成黑盒 API而是暴露每一步的流量张量、路径集、边际成本向量让你看清“均衡”到底是怎么一滴一滴挤出来的。2. UE 模型与 Frank-Wolfe 算法为什么选它为什么必须手写核心循环UEUser Equilibrium模型的本质是假设每个出行者都选择当前网络中对自己而言最短路径而所有被选中的路径在该 OD 对之间具有相等且最小的广义费用。这导致系统状态满足一个变分不等式VI$$ \sum_{a \in A} t_a(v_a^) (v_a - v_a^) \geq 0, \quad \forall \mathbf{v} \in \Omega $$其中 $t_a(v_a)$ 是路段 a 的行程时间函数常用 BPR 形式$\Omega$ 是满足 OD 需求约束的可行流量集合。直接求解这个 VI 非常困难而 Frank-Wolfe 提供了一种基于方向导数的迭代下降策略每步构造一个线性近似求解该近似下的极小点即全网络最短路分配再沿连接当前解与该极小点的方向做一维搜索更新。提示这不是 SGD也不是 Adam。FW 的“方向”是确定性的——它永远指向当前流量下所有 OD 对的系统最优SO分配方向而 SO 分配恰好对应于将所有出行者视为一个整体进行全局优化的结果。UE 和 SO 的差距正是“自私路径选择”带来的系统效率损失即 Braess 悖论的量化基础。所以 FW 迭代过程本身就是在不断逼近“如果所有人合作路网能多快”和“所有人各自为政路网实际多慢”之间的平衡点。2.1 SiouxFalls 网络小而全的交通分配“Hello World”SiouxFalls 是一个仅含 24 个节点、76 条有向路段、24 个 OD 对的合成路网但它完整包含了环形结构、并联路段、瓶颈节点等典型特征。其标准数据集包含三部分SiouxFalls_node.csv节点 ID、坐标x, ySiouxFalls_link.csv路段 ID、起点/终点节点、自由流时间FFTT、容量Cap、BPR 参数 α0.15, β4.0SiouxFalls_trips.csvOD 对origin, destination、出行需求量trips我们不依赖任何交通仿真平台如 VISUM、Aimsun导出中间格式而是直接用 Pandas 读取 CSV构建纯内存图结构。关键在于路段索引必须全局唯一且连续因为后续所有流量向量v、时间向量t、雅可比近似矩阵都以该索引对齐。import pandas as pd import numpy as np # 读取路段数据强制重置索引为 0~75 links pd.read_csv(SiouxFalls_link.csv) links links.reset_index(dropTrue) n_links len(links) # 构建邻接表node - list of (next_node, link_idx) adj {} for idx, row in links.iterrows(): u, v int(row[from]), int(row[to]) if u not in adj: adj[u] [] adj[u].append((v, idx)) # OD 需求矩阵24x24行origin列destination trips_df pd.read_csv(SiouxFalls_trips.csv) od_pairs trips_df[[origin, destination, trips]].values n_nodes 24 demand np.zeros((n_nodes, n_nodes)) for o, d, q in od_pairs: demand[int(o)-1, int(d)-1] q # 节点编号从 1 开始数组索引从 0 开始这段代码看似简单但决定了后续所有路径搜索和流量累加的正确性。adj是 Dijkstra 的输入基础demand是 FW 外层循环的驱动源而n_links将作为所有一维流量向量v的长度——这是整个实现的“空间契约”一旦错位后续所有矩阵运算都会静默出错。2.2 BPR 函数与时间-流量映射别让 α 和 β 成为玄学参数BPRBureau of Public Roads函数是交通工程中最经典的路段阻抗模型$$ t_a(v_a) t_a^0 \left(1 \alpha \left(\frac{v_a}{c_a}\right)^\beta \right) $$其中 $t_a^0$ 是自由流时间$c_a$ 是通行能力$\alpha$ 和 $\beta$ 是经验校准参数。在 SiouxFalls 标准数据中$\alpha0.15$, $\beta4.0$ 是固定值不可更改——这不是调参而是复现基准实验的前提。很多初学者试图“优化”这两个值来拟合实测数据但在本项目中它们是定义问题的一部分改了就不再是 SiouxFalls UE 问题。我们封装一个向量化函数支持对整条路段流量向量v一次性计算时间向量tdef bpr_time(v, fftt, capacity, alpha0.15, beta4.0): 向量化计算 BPR 时间 :param v: (n_links,) 流量向量 :param fftt: (n_links,) 自由流时间向量 :param capacity: (n_links,) 容量向量 :return: (n_links,) 当前流量下的行程时间向量 # 防止除零capacity 为 0 的路段设为极大值实际网络中不会出现 cap_safe np.where(capacity 0, 1e6, capacity) ratio np.clip(v / cap_safe, 0, 100) # 截断极端拥堵比避免浮点溢出 return fftt * (1 alpha * (ratio ** beta)) # 示例用初始零流量计算自由流时间 fftt links[fftt].values.astype(float) capacity links[capacity].values.astype(float) t0 bpr_time(np.zeros(n_links), fftt, capacity) # 所有路段均为自由流时间注意np.clip——当v/c超过 100 时ratio**4会达到 1e8乘上fftt可能触发inf进而污染后续最短路计算。这是真实踩坑点某次调试中因某条路段容量被误读为 0导致v/c无穷大Dijkstra 返回inf路径FW 方向向量全为 0迭代卡死。clip是低成本的“后悔药”。2.3 Frank-Wolfe 主循环四步铁律缺一不可FW 求解 UE 的每一次迭代严格遵循以下四步不是三步不是五步评估当前流量下的路段时间$t(v^k)$求解方向子问题在当前时间下求所有 OD 对的最短路径并按 OD 需求分配流量 → 得到辅助流量 $y^k$一维搜索找最优步长 $\gamma_k \arg\min_{\gamma \in [0,1]} ; C(v^k \gamma (y^k - v^k))$其中 $C(v)$ 是系统总时间SCT更新流量$v^{k1} v^k \gamma_k (y^k - v^k)$关键在于第 2 步的 $y^k$ 不是“新路径流量”而是“若所有人此刻重新选择最短路系统会变成什么样”的瞬时状态第 3 步的 $\gamma_k$ 决定了我们向这个理想状态“走多远”。FW 收敛慢但稳定通常需 30~100 次迭代而 $\gamma_k$ 若用精确线搜索求导令为 0可大幅加速。def fw_iteration(v, demand, adj, fftt, capacity, alpha0.15, beta4.0): 执行一次 Frank-Wolfe 迭代 :return: new_v, gamma, sct_old, sct_new, y # Step 1: 计算当前时间 t bpr_time(v, fftt, capacity, alpha, beta) # Step 2: 求解最短路并分配y^k y all_or_nothing_assignment(demand, adj, t, n_nodes, n_links) # Step 3: 精确线搜索minimize SCT(gamma) sum_i t_i(v gamma*(y-v)) * (v_i gamma*(y_i-v_i)) # SCT(gamma) 是关于 gamma 的四次函数但可解析求导得三次方程用 numpy.roots 求实根 def sct_func(gamma): v_new v gamma * (y - v) t_new bpr_time(v_new, fftt, capacity, alpha, beta) return np.sum(t_new * v_new) # 求导dSCT/dgamma 0 → 三次方程系数推导略见附录公式 # 这里用数值方法在 [0,1] 上采样 100 点取最小值点更鲁棒 gammas np.linspace(0, 1, 101) scts np.array([sct_func(g) for g in gammas]) gamma_opt gammas[np.argmin(scts)] # Step 4: 更新 v_new v gamma_opt * (y - v) sct_old sct_func(0) sct_new sct_func(gamma_opt) return v_new, gamma_opt, sct_old, sct_new, y # 初始化流量为零 v np.zeros(n_links) max_iter 80 convergence_tol 1e-5 sct_history [] for k in range(max_iter): v, gamma, sct_old, sct_new, y fw_iteration(v, demand, adj, fftt, capacity) sct_history.append(sct_new) # 相对变化率作为收敛判据 if k 0 and abs(sct_old - sct_new) / sct_old convergence_tol: print(fConverged at iteration {k}, SCT{sct_new:.2f}) breakall_or_nothing_assignment是核心子程序需对每个 OD 对运行 Dijkstra生成路径再将需求量全量加载到路径上所有路段。它不返回路径列表只返回一个(n_links,)的辅助流量向量y。这是 FW 的精髓你不需要存储所有路径只需要知道“如果重选每条路段会新增多少车”。3. 路径生成与全有全无分配Dijkstra 必须返回路径而不是只返回距离Frank-Wolfe 的“方向子问题”本质是求解一个线性规划$$ \min_{y \in \Omega} ; \sum_{a} t_a(v^k) \cdot y_a $$其解必为极点extreme point即所有 OD 流量都走唯一最短路径All-or-Nothing, AON。因此我们必须为每个 OD 对找到一条最短路径并将全部需求量加载到该路径的各路段上。这就要求 Dijkstra 算法必须返回路径本身而不仅是距离。3.1 带路径重建的 Dijkstra用 predecessor 数组代替 distance 数组标准 Dijkstra 只维护dist[v]但我们还需prev[v]记录到达v的前驱节点。当目标节点d被标记为已访问时通过prev数组反向追溯即可得到完整路径。import heapq def dijkstra_path(adj, t, start, end, n_nodes): 返回从 start 到 end 的最短路径节点序列及总时间 :param adj: 邻接表 {u: [(v, link_idx), ...]} :param t: (n_links,) 当前路段时间向量 :param start, end: 起终点节点 ID0-based :return: path_nodes: [start, ..., end], total_time dist [float(inf)] * n_nodes prev [-1] * n_nodes # 记录前驱节点 dist[start] 0 heap [(0.0, start)] while heap: d, u heapq.heappop(heap) if d dist[u]: continue if u end: break if u not in adj: continue for v, link_idx in adj[u]: # 边权 路段时间 t[link_idx] new_dist d t[link_idx] if new_dist dist[v]: dist[v] new_dist prev[v] u heapq.heappush(heap, (new_dist, v)) # 重建路径 if dist[end] float(inf): return [], float(inf) path [] curr end while curr ! -1: path.append(curr) curr prev[curr] path.reverse() return path, dist[end] # 示例从节点 0 到节点 1 的最短路径 path, time_cost dijkstra_path(adj, t0, 0, 1, n_nodes) print(fPath 0-1: {path}, time{time_cost:.2f})注意adj[u]中的link_idx是路段在全局links表中的索引t[link_idx]即该路段当前时间。这是路径与时间向量对齐的关键纽带。3.2 全有全无分配路径→路段流量的映射必须可逆有了路径节点序列下一步是将其转换为路段 ID 序列再将 OD 需求量q加载到这些路段上。难点在于同一节点对之间可能存在多条最短路径时间相同而标准 Dijkstra 只返回一条。UE 理论允许这种歧义但实践中我们采用“首次发现即采纳”策略确保结果可复现。def path_to_links(path, adj): 将节点路径转换为路段 ID 列表 :param path: [n0, n1, n2, ..., nk] :return: list of link indices links_on_path [] for i in range(len(path) - 1): u, v path[i], path[i1] # 在 adj[u] 中查找指向 v 的路段 found False for next_node, link_idx in adj[u]: if next_node v: links_on_path.append(link_idx) found True break if not found: raise ValueError(fNo link from node {u} to {v} in adjacency list) return links_on_path def all_or_nothing_assignment(demand, adj, t, n_nodes, n_links): 执行全有全无分配返回 (n_links,) 辅助流量向量 y y np.zeros(n_links) # 遍历所有 OD 对 for o in range(n_nodes): for d in range(n_nodes): q demand[o, d] if q 0: continue # 求 o-d 最短路径 path, _ dijkstra_path(adj, t, o, d, n_nodes) if not path: continue # 无连通路径跳过 # 转换为路段 ID link_ids path_to_links(path, adj) # 将需求 q 加载到路径上所有路段 for lid in link_ids: y[lid] q return y这里path_to_links是关键胶水函数。它确保了无论路径多长最终都归结为对y[lid]的原子累加。y向量的维度76与v完全一致为 FW 主循环的向量运算铺平道路。3.3 避坑路径歧义、浮点误差与图连通性三大雷区现象某次迭代后SCT 不降反升或震荡不收敛原因Dijkstra 中heapq.heappush(heap, (new_dist, v))的new_dist是浮点数当多条路径时间极其接近如 12.0000001 vs 12.0000002时浮点精度导致“错误”路径被优先弹出prev[v]被错误覆盖。解决在new_dist dist[v]判断前加入容差1e-8if new_dist dist[v] - 1e-8:。这是交通分配中公认的浮点安全实践。现象all_or_nothing_assignment返回的y中某些路段流量为 0但该路段明明在 OD 路径上原因path_to_links查找u-v时遍历adj[u]未 break导致只取了最后一条匹配边如果有多条平行路段。解决for next_node, link_idx in adj[u]:循环内if next_node v:后立即break确保取第一条即原始 CSV 中顺序。现象dijkstra_path返回空路径all_or_nothing_assignment报ValueError原因SiouxFalls 网络虽小但存在单向路段。adj表只存u-v未存v-u。当 OD 对o10, d5时若adj[10]中无指向 5 的边且adj[5]中也无指向 10 的边因为是单向则无路径。但 SiouxFalls 标准数据中所有 OD 对都是连通的。解决检查SiouxFalls_link.csv中from/to列是否被误读为字符串导致int(row[from])失败或adj构建时漏掉某些from节点如from0但adj键为0字符串。用print(set(links[from]))与set(adj.keys())对比可快速定位。4. 收敛性监控与结果验证别信“跑了 50 步”要看 SCT 曲线和路径一致性Frank-Wolfe 的收敛性不能只看迭代次数。UE 解的正确性体现在三个可验证层面系统总时间SCT单调下降且趋于平稳、各 OD 对最短路径时间相等、路段流量满足守恒律。忽略任一层面都可能得到一个“看起来像解”的错误结果。4.1 SCT 历史曲线下降速率告诉你算法是否健康SCTSystem Cost Total定义为 $\sum_a t_a(v_a) \cdot v_a$即所有路段的“时间×流量”之和。它是 UE 的目标函数虽然 UE 本身不是最小化 SCT但 FW 迭代路径恰好使 SCT 单调递减。绘制sct_history是第一道防线。import matplotlib.pyplot as plt plt.figure(figsize(10, 5)) plt.plot(sct_history, b-o, markersize3) plt.xlabel(Iteration) plt.ylabel(System Total Time (veh·hr)) plt.title(Frank-Wolfe Convergence on SiouxFalls) plt.grid(True) plt.yscale(log) # 对数纵轴看清后期微小变化 plt.show()健康曲线应呈现前 10 步快速下降大步长10~30 步缓慢下降步长缩小30 步后趋于水平相对变化 1e-5。若曲线出现上扬说明bpr_time或all_or_nothing_assignment有 bug若 50 步后仍剧烈波动可能是gamma搜索不精确或浮点误差累积。4.2 OD 对路径时间一致性检验UE 的黄金标准UE 的核心定义是对任意 OD 对 $(o,d)$所有被使用的路径 $p$ 满足 $c_p(v) \min_{p \in P_{od}} c_{p}(v)$即路径广义费用相等。我们选取几个高频 OD 对如 1→2, 1→10, 10→1用最终收敛的流量v_final重新运行 Dijkstra记录所有最短路径的时间并验证是否一致。def check_od_consistency(v_final, adj, fftt, capacity, n_nodes, od_list[(0,1), (0,9), (9,0)]): 检查指定 OD 对的最短路径时间是否一致 t_final bpr_time(v_final, fftt, capacity) for o, d in od_list: # 获取所有最短路径不止一条 paths k_shortest_paths(adj, t_final, o, d, n_nodes, k5) costs [] for path in paths: if not path: continue cost 0 for i in range(len(path)-1): u, v path[i], path[i1] for next_node, link_idx in adj[u]: if next_node v: cost t_final[link_idx] break costs.append(cost) if costs: print(fOD {o1}-{d1}: min cost {min(costs):.4f}, max cost {max(costs):.4f}, frange {max(costs)-min(costs):.6f}) # 注意k_shortest_paths 需另行实现如 Yens algorithm此处为示意 check_od_consistency(v, adj, fftt, capacity, n_nodes)理想情况下range应 ≤ 1e-4。若达 0.1 以上说明流量分配未真正收敛或路径搜索未考虑所有等价路径。4.3 路段流量守恒用 OD 矩阵反推验证你的v是否合法对任意路段a其流量v_a必须等于所有经过a的 OD 路径的流量之和。我们可以用demand和最终路径集反向计算一个v_recon并与v_final对比。def reconstruct_flow_from_paths(v_final, demand, adj, fftt, capacity, n_nodes, n_links): 用最终流量 v_final 重新计算所有 OD 最短路径再累加得到 v_recon t_final bpr_time(v_final, fftt, capacity) v_recon np.zeros(n_links) for o in range(n_nodes): for d in range(n_nodes): q demand[o, d] if q 0: continue path, _ dijkstra_path(adj, t_final, o, d, n_nodes) if not path: continue link_ids path_to_links(path, adj) for lid in link_ids: v_recon[lid] q return v_recon v_recon reconstruct_flow_from_paths(v, demand, adj, fftt, capacity, n_nodes, n_links) error np.max(np.abs(v - v_recon)) print(fMax reconstruction error: {error:.2e})error应 1e-3。这是对all_or_nothing_assignment和dijkstra_path的终极交叉验证——它不依赖任何中间变量只用输入demand和输出v独立重跑一遍分配逻辑。5. 进阶技巧如何用这份代码做你自己的路网分析三个实战习惯这份 SiouxFalls 实现的价值远不止于跑通一个例子。它是一套可迁移的“交通分配脚手架”我把它用在三个真实场景中某高校交通实验室的课程设计、某公司定制化仿真引擎的 UE 模块验证、以及我自己写论文时的 baseline 对比。每次使用我都强制执行以下三个动作从未翻车。5.1 动作一用“零流量初始化 单 OD 测试”隔离路径模块在接入真实路网前我总会先构造一个极简测试网3 个节点 A-B-C两条路段 A→B、B→COD 需求仅 A→C100。此时理论解是v_AB100, v_BC100, t_ABt0_AB*(1α*(100/c_AB)^β)。我把v初始化为零只跑 1 次 FW 迭代然后手动检查dijkstra_path是否返回[A,B,C]path_to_links是否返回[0,1]假设 AB 是 link 0BC 是 link 1y向量是否为[100,100,0,...]v_new是否为[100*γ, 100*γ, 0,...]这 5 分钟的手动验证能提前揪出 80% 的路径映射错误。比对着 76 维向量 debug 高效十倍。5.2 动作二把 BPR 参数做成配置项但绝不用于“调优”我在bpr_time函数签名中保留alpha和beta参数但实际调用时永远传入0.15和4.0。为什么还留着因为当我要分析“如果政策降低 α如推广导航软件提升路径信息质量UE 流量会如何变化”时只需改两个数字重跑即可。但必须清楚这不是模型校准而是what-if 敏感性分析。真正的校准需用实测行程时间反推那属于另一个专业流程不在本脚手架范围内。5.3 动作三保存每轮迭代的v和t用于可视化与调试我修改主循环增加v_history [v.copy()] # 存储每次迭代后的 v t_history [bpr_time(v, fftt, capacity).copy()] for k in range(max_iter): v, gamma, sct_old, sct_new, y fw_iteration(v, demand, adj, fftt, capacity) v_history.append(v.copy()) t_history.append(bpr_time(v, fftt, capacity).copy()) # ... convergence check有了v_history我就能用matplotlib.animation画出流量在路网上的动态演化用np.linalg.norm(v_history[-1] - v_history[-2])替代 SCT 差值作为收敛判据甚至把v_history[10]导出为 CSV导入 GIS 软件看“第 10 步时哪几条路已经堵死”。数据可追溯是工程落地的生命线。从那以后我每次接手新路网都先跑通 SiouxFalls再替换links.csv和trips.csv最后只改三处n_nodes、n_links、adj构建逻辑。其他 95% 的代码原封不动。这份克制省下了我至少 200 小时的重复造轮时间。希望帮到你。本文还有配套的精品资源点击获取