简介面向交通工程、交通网络分析与交通分配方向的学习者这份材料以经典SiouxFalls网络为案例用Python实现Frank-Wolfe算法求解UE用户均衡帮助理解“无用户可通过单方面改变路径降低出行费用”的均衡条件。包内共7个文件含6个txt与1个py压缩包约7KBtxt主要用于存放节点、路段、OD需求及网络基础参数py文件则对应Frank-Wolfe迭代求解流程结构轻量但覆盖完整实验链路。已有2335人学习适合课程实验、算法入门与结果复现。读者可据此走通从读取网络数据、套用BPR路阻函数、加载OD出行需求到输出各路段流量的全过程并借均衡结果核对分配逻辑与收敛表现为后续交通网络分析或科研基线实验提供可操作的参考。1. SiouxFalls 上的用户均衡为什么值得亲手实现一遍SiouxFalls 这个 24 节点、76 条路段的算例是大多数人接触用户均衡User EquilibriumUE时跑的第一个网络。它小到能在一张纸上画出拓扑又复杂到足以让随手写的分配结果和真正的均衡解差出可观的总阻抗。用 Frank-Wolfe 算法在 Python 里把它跑通最直接的收益是你能亲眼看到“没人后悔自己选的路”这种流量分布是怎么从一次全有全无分配里长出来的也能摸清收敛间隙、线搜索步长、BPR 阻抗这几个量之间的拉扯关系。做交通分配、路网评估或者要给路径规划找基准解的开发者都用得上。下面按模型为什么长这样、数据怎么组织、主循环怎么写、哪里会出问题的顺序走一遍。2. UE 与 Frank-Wolfe为什么每轮只需要做一次最短路2.1 Wardrop 第一原理和 Beckmann 目标函数是同一件事用户均衡的原始表述很朴素每个 OD 对之间所有被使用的路径阻抗相等且不大于任何未被使用路径的阻抗。这句话写成数学就是一组互补松弛条件——路径流量乘上“该路径阻抗减去最小阻抗”等于零。互补松弛直接解非常麻烦因为路径集合的规模随网络指数增长SiouxFalls 才 24 个节点路径数就已经上千。Beckmann 在 1956 年给出的转换是把路段流量 x 当成决策变量构造min Z(x) Σ_a ∫_0^{x_a} t_a(w) dw s.t. Σ_{a∈A} δ_{a,r} f_r q_od 对每个 OD 对 f_r ≥ 0其中 t_a(·) 是路段阻抗函数δ 是路径与路段的关联矩阵f_r 是路径流量q_od 是 OD 需求。这个规划的目标函数是严格凸的因为 t_a 单调递增可行域是线性多面体所以它有唯一的最优路段流量解而且 KKT 条件展开后恰好就是 Wardrop 第一原理。换句话说解这个凸规划等于解用户均衡。这里必须把 BPR 阻抗函数摆出来因为它是整个模型里唯一一个“可以有不同意见”的地方t_a(x_a) t_a^0 · (1 α · (x_a / c_a)^β)t_a^0 是自由流走行时间c_a 是通行能力α 通常取 0.15β 通常取 4。SiouxFalls 的公开算例文件里给的就是 b0.15、power4我一般不去改这两个数——改了之后你算出来的结果没法跟任何公开基准对照出了问题也分不清是代码错了还是参数错了。把 BPR 代入 Beckmann 目标函数积分可以解析算出来Z(x) Σ_a t_a^0 · [ x_a α · x_a^(β1) / ((β1) · c_a^β) ]这个解析式在后面有两个用处一是画收敛曲线看目标函数是不是单调不增二是当间隙算出来可疑时用它交叉验算线搜索有没有把解带偏。很多人只算间隙不看目标函数结果线搜索写错了半个月都没发现这类问题在后面的排查章里还会细说。符号含义SiouxFalls 取值t_a^0自由流走行时间分钟1 到 10 不等c_a通行能力veh/h4432 到 25900α, βBPR 参数0.15, 4q_odOD 需求非零对 528 个路段数有向边7638 对双向节点数—242.2 梯度就是路段阻抗子问题因此退化成全有全无分配Frank-Wolfe 是求解“线性约束 凸目标”这类问题的经典方法它的迭代格式是x_{k1} x_k γ_k · (y_k - x_k), γ_k ∈ [0, 1] y_k argmin_{y ∈ D} ∇Z(x_k)^T · yD 是流量守恒约束构成的可行域。关键在于这个子问题怎么解。对 Beckmann 目标函数求偏导∂Z/∂x_a t_a(x_a)梯度向量就是当前路段流量下的阻抗向量。于是子问题变成min Σ_a t_a(x_a^k) · y_a约束是 y 也要满足每个 OD 对的需求守恒。这个式子的含义极其直白把路段阻抗固定成当前解算出来的那一组值然后每个 OD 对的需求全部走最短路这样总阻抗最小。这就是全有全无分配All-or-NothingAON也叫最短路分配。这个推导带来一个非常实际的结论Frank-Wolfe 每迭代一轮只需要跑一遍“每个起点一次 Dijkstra 沿最短路树加载需求”不需要存路径集合不需要解线性方程组内存占用只跟路段数和节点数有关。对 SiouxFalls 这种规模一轮迭代在普通笔记本上是毫秒级对几万个路段的路网一轮也还在秒级。这是它至今仍被当作交通分配基准算法的原因。代价也很明确。AON 解永远是可行域的顶点——所有需求都压在最短路上其他路段流量为零。F-W 是在这些顶点之间做凸组合迭代点会沿着顶点连成的折线逼近最优解收敛速度是 O(1/k) 的次线性。翻译成人话前 10 轮间隙掉得飞快从 1e-1 掉到 1e-2 可能只要几轮但从 1e-3 掉到 1e-4 可能要几十轮再往下就是漫长的尾部。这个特性不是代码写错了是算法本身的性质后面讲收敛判据时还会回来。还有一点容易被忽略F-W 的子问题解 y_k 顺带给出了当前 OD 需求下的最短路径总阻抗。因为 y 就是 AON 解Σ_a t_a(x_a^k) · y_a 恰好等于所有 OD 需求乘以各自最短路径阻抗之和。这个量叫 SPTTShortest Path Travel Time它会在算收敛间隙时直接用到不需要额外算一遍最短路。有些实现每轮跑两次 AON一次求方向一次求间隙纯属浪费翻倍的运行时间换不来任何精度。2.3 步长 γ 从哪来一维凸函数上的二分求根方向定了剩下的是在 x_k 到 y_k 这条线段上找最优点φ(γ) Z(x_k γ · (y_k - x_k)), γ ∈ [0, 1]把 φ 对 γ 求导φ(γ) Σ_a (y_a - x_a) · t_a(x_a γ · (y_a - x_a))再求一次导φ(γ) Σ_a (y_a - x_a)^2 · t_a(x_a γ(y_a - x_a))BPR 函数单调递增所以 t_a ≥ 0于是 φ(γ) ≥ 0φ 是凸函数φ 单调不减。这就意味着找 φ 的最小值点等价于找 φ(γ) 0 的根而单调函数的求根用二分最稳不需要导数信息也不会像牛顿法那样在接近 0 的时候炸掉。两个端点值得先判断。φ(0) Σ_a (y_a - x_a) t_a(x_a) SPTT - TSTT其中 TSTT Σ_a x_a t_a(x_a) 是当前解的总阻抗。因为 y 是在 t(x) 下使总阻抗最小的解所以 SPTT ≤ TSTTφ(0) ≤ 0 恒成立。如果 φ(0) 恰好等于 0说明当前解已经均衡直接停。另一端 φ(1) 可正可负如果它 ≤ 0说明 φ 在 [0,1] 上单调递减直接取 γ 1。下面这段是阻抗、目标函数和线搜索的最小实现可以直接跑import numpy as np ALPHA_BPR, BETA_BPR 0.15, 4.0 def link_cost(t0, cap, x, alphaALPHA_BPR, betaBETA_BPR): BPR 阻抗t0 * (1 alpha * (v/c)^beta)全向量化 ratio np.maximum(x, 0.0) / cap # 防浮点 -1e-18 导致 power 出 nan return t0 * (1.0 alpha * np.power(ratio, beta)) def beckmann(t0, cap, x, alphaALPHA_BPR, betaBETA_BPR): Beckmann 目标的解析积分用来画收敛曲线和交叉验算 return float(np.sum( t0 * (x alpha * np.power(x, beta 1) / ((beta 1) * np.power(cap, beta))) )) def line_search(t0, cap, x, y, tol1e-6, max_iter60): 在 [0,1] 上二分求 phi(gamma)0phi 单调不减二分必收敛 d y - x def dphi(g): return float(np.dot(d, link_cost(t0, cap, x g * d))) if dphi(0.0) 0.0: return 0.0 # 已经均衡方向不下降 if dphi(1.0) 0.0: return 1.0 # 端点最优走满步 lo, hi 0.0, 1.0 for _ in range(max_iter): mid 0.5 * (lo hi) if dphi(mid) 0.0: hi mid else: lo mid if hi - lo tol: break return 0.5 * (lo hi)几个参数要说清楚。tol1e-6 是步长精度对应到间隙上的影响大约是同阶甚至更小对 gap_tol1e-4 的目标完全够用如果想跑到 1e-6把 tol 收到 1e-9 再配合更多迭代次数。max_iter60 是保险丝实际二分 20 次左右就已经把区间压到 1e-6 以下正常不会跑满。link_cost 里那个 np.maximum 看着多余因为 γ∈[0,1] 且 x、y 都非负理论上 xγd 不会小于 0但浮点误差会给出 -1e-17 这种值np.power 遇到负数加非整数指数会返回 nan一个 nan 进去整条流量向量就废了这种错误在日志里表现得极其隐晦。另外提醒一句φ(0) 的计算里已经包含了 SPTT如果线搜索函数内部顺手把它返回出来主循环就不必重算一遍能省掉一次全向量乘法。SiouxFalls 上无所谓几万路段时这点开销会积少成多。3. SiouxFalls 的数据组织路段表、OD 表和邻接表3.1 三张表的字段约定以及两套编号必须分清SiouxFalls 的数据可以拆成三块节点表、路段表、OD 表。节点表在纯分配问题里其实用不上因为 F-W 不需要坐标真正要喂给算法的是路段表和 OD 表。路段表的核心四列字段类型说明tailint1 到 24路段起点节点号headint1 到 24路段终点节点号capacityfloatveh/h通行能力必须和流量同单位free_flow_timefloatmin自由流走行时间OD 表三列origin、dest、demand单位是 veh/h。SiouxFalls 的公开 OD 需求总量是 360600非零 OD 对 528 个24×23552 对里有一部分需求为零。这里有一个必须提前说清的约定问题节点编号是 1-based1 到 24但 NumPy 数组下标是 0-based。路段我用 0-based 的 link_id 作数组下标节点号保持 1-based 并给数组开到 n_nodes1 长度多出来的一格空着不用。这样 Dijkstra 里的 dist[v]、prev[v] 可以直接用节点号索引写起来最不容易错。最怕的是混着用——一半代码 0-based 一半 1-based跑出来的结果要么越界要么悄悄错位而且错位之后流量守恒还是对的查起来极其费劲。路段表我直接内嵌进代码避免依赖外部文件# SiouxFalls 76 条有向路段(tail, head, capacity, free_flow_time) # 容量保留 1 位小数对结果的影响在 1e-9 量级 SIOUX_LINKS [ (1, 2, 25900.2, 6), (1, 3, 23403.5, 4), (2, 1, 25900.2, 6), (2, 6, 4958.2, 5), (3, 1, 23403.5, 4), (3, 4, 17110.5, p a hrefhttps://download.csdn.net/download/Libertasss/86734295 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p
