Python实现SiouxFalls用户均衡与Frank-Wolfe算法

简介：面向交通工程、交通网络分析与交通分配方向的学习者，这份材料以经典SiouxFalls网络为案例，用Python实现Frank-Wolfe算法求解UE用户均衡，帮助理解“无用户可通过单方面改变路径降低出行费用”的均衡条件。包内共7个文件，含6个txt与1个py，压缩包约7KB；txt主要用于存放节点、路段、OD需求及网络基础参数，py文件则对应Frank-Wolfe迭代求解流程，结构轻量但覆盖完整实验链路。已有2335人学习，适合课程实验、算法入门与结果复现。读者可据此走通从读取网络数据、套用BPR路阻函数、加载OD出行需求到输出各路段流量的全过程，并借均衡结果核对分配逻辑与收敛表现，为后续交通网络分析或科研基线实验提供可操作的参考。

1. SiouxFalls 上的用户均衡，为什么值得亲手实现一遍

SiouxFalls 这个 24 节点、76 条路段的算例，是大多数人接触用户均衡（User Equilibrium，UE）时跑的第一个网络。它小到能在一张纸上画出拓扑，又复杂到足以让随手写的分配结果和真正的均衡解差出可观的总阻抗。用 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 的公开算例文件里给的就是 b=0.15、power=4，我一般不去改这两个数——改了之后你算出来的结果没法跟任何公开基准对照，出了问题也分不清是代码错了还是参数错了。

把 BPR 代入 Beckmann 目标函数，积分可以解析算出来：

Z(x) = Σ_a t_a^0 · [ x_a + α · x_a^(β+1) / ((β+1) · c_a^β) ]

这个解析式在后面有两个用处：一是画收敛曲线，看目标函数是不是单调不增；二是当间隙算出来可疑时，用它交叉验算线搜索有没有把解带偏。很多人只算间隙不看目标函数，结果线搜索写错了半个月都没发现，这类问题在后面的排查章里还会细说。

2.2 梯度就是路段阻抗，子问题因此退化成全有全无分配

Frank-Wolfe 是求解“线性约束 + 凸目标”这类问题的经典方法，它的迭代格式是

x_{k+1} = x_k + γ_k · (y_k - x_k), γ_k ∈ [0, 1] y_k = argmin_{y ∈ D} ∇Z(x_k)^T · y

D 是流量守恒约束构成的可行域。关键在于这个子问题怎么解。对 Beckmann 目标函数求偏导：

∂Z/∂x_a = t_a(x_a)

梯度向量就是当前路段流量下的阻抗向量。于是子问题变成

min Σ_a t_a(x_a^k) · y_a

约束是 y 也要满足每个 OD 对的需求守恒。这个式子的含义极其直白：把路段阻抗固定成当前解算出来的那一组值，然后每个 OD 对的需求全部走最短路，这样总阻抗最小。这就是全有全无分配（All-or-Nothing，AON），也叫最短路分配。

这个推导带来一个非常实际的结论：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 需求乘以各自最短路径阻抗之和。这个量叫 SPTT（Shortest 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, alpha=ALPHA_BPR, beta=BETA_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, alpha=ALPHA_BPR, beta=BETA_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, tol=1e-6, max_iter=60): """在 [0,1] 上二分求 phi'(gamma)=0；phi' 单调不减，二分必收敛""" 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)

几个参数要说清楚。tol=1e-6 是步长精度，对应到间隙上的影响大约是同阶甚至更小，对 gap_tol=1e-4 的目标完全够用；如果想跑到 1e-6，把 tol 收到 1e-9 再配合更多迭代次数。max_iter=60 是保险丝，实际二分 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 表。

路段表的核心四列：

OD 表三列：origin、dest、demand，单位是 veh/h。SiouxFalls 的公开 OD 需求总量是 360600，非零 OD 对 528 个（24×23=552 对里有一部分需求为零）。

这里有一个必须提前说清的约定问题：节点编号是 1-based（1 到 24），但 NumPy 数组下标是 0-based。路段我用 0-based 的 link_id 作数组下标，节点号保持 1-based 并给数组开到 n_nodes+1 长度，多出来的一格空着不用。这样 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 href="https://download.csdn.net/download/Libertasss/86734295" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>