SiouxFalls路网用户均衡交通分配的Frank-Wolfe实现

简介：围绕Frank-Wolfe算法求解UE用户均衡模型，这份Python实现以SiouxFalls网络为案例，面向交通工程、交通网络分析与交通分配方向的学习者，帮助理解如何通过迭代调整路径流量逐步逼近用户均衡状态，并验证没有出行者可单方面改路径降低费用。压缩包共7个文件，整体约7KB，以6个txt和1个py为主：文本文件提供网络拓扑、节点、路段及OD需求等基础数据，Python脚本完成BPR路阻计算、最短路分配与Frank-Wolfe迭代更新，最终输出各路段流量结果。已有2335人学习，适合课程实验、算法复现与结果对照。读者可据此跑通从数据读取、方向搜索到流量更新的完整流程，观察收敛过程与路段流量变化，并为改换网络、调整需求或对比其他交通分配算法提供可复用的代码与数据基础，对课程设计、论文实验与算法入门均有参考价值。

1. 一张 76 条边的路网，为什么值得单独写一套 Frank-Wolfe

SiouxFalls 是交通分配练手绕不开的一张网：24 个节点、76 条无向路段、528 个 OD 对，总需求 360600。规模不大，但它把 UE（User Equilibrium，用户均衡）里容易出错的地方全暴露了一遍——BPR 阻抗、全有全无分配、线搜索、收敛间隙，一个都躲不掉。这套 Python 实现做的事很具体：用 Frank-Wolfe 算法把 Beckmann 目标函数压下去，让相对间隙稳定跑到 1e-4 以下，最后输出每条弧的均衡流量和 V/C 比。做交通网络分析课程设计的人能直接照着跑，手上有别的路网数据、想快速拿一个可信基准解做对照的从业者也能用。下面从数据层开始，一层层把能复现的细节摊开。

2. 数据层建模：把 76 条无向边展开成 152 条有向弧

2.1 三张表：路段表、OD 表、弧表

标准 SiouxFalls 数据集给的是两份东西：一份 76 行的路段表，每行一条无向边；一份 528 行的 OD 表，每行一个起讫点对和需求量。不同来源的字段名不完全一致，我一般先统一成下面这套口径，后面所有代码都按它走，换网络时只改读入部分，算法层一行不动。

为什么不直接在原始路段表上跑？因为 Frank-Wolfe 的每一次迭代都在「弧流量向量」上做加减、点乘、求导，全部是定长数组运算。如果每次都要拿 (from, to) 去字典里查索引，528 个 OD 对乘上几千次迭代，开销会翻好几倍。派生出一张弧表，把 t0 和 cap 提成两个 numpy 数组，后面就只剩向量化操作了。

2.2 无向边转有向弧：一次复制，两个弧号

这一步是整个数据层最容易埋雷的地方，先看代码。

import numpy as np import pandas as pd def load_network(net_csv, od_csv): """读入 SiouxFalls 路段表与 OD 表，展开成有向弧。 net_csv 需含: from_node, to_node, length, free_flow_time, capacity od_csv 需含: o, d, demand """ net = pd.read_csv(net_csv) od = pd.read_csv(od_csv) # 字段名统一小写并去空格，避免不同来源的 CSV 大小写不一致 net.columns = [c.strip().lower() for c in net.columns] od.columns = [c.strip().lower() for c in od.columns] keep = ['from_node', 'to_node', 'length', 'free_flow_time', 'capacity'] fwd = net[keep].copy() # 反向弧：起终点互换，长度/时间/容量不变（SiouxFalls 双向同规格） bwd = fwd.rename(columns={'from_node': 'to_node', 'to_node': 'from_node'}) arcs = pd.concat([fwd, bwd], ignore_index=True) arcs['arc_id'] = np.arange(len(arcs), dtype=int) arcs['t0'] = arcs['free_flow_time'].astype(float) arcs['cap'] = arcs['capacity'].astype(float) # pair_id 用来把正反两条弧归回同一条物理路段，出结果时按路段聚合 arcs['pair_id'] = np.tile(np.arange(len(fwd)), 2) return arcs.reset_index(drop=True), od

pd.concat之后行数从 76 变成 152，前 76 条是正向弧，后 76 条是反向弧，arc_id连续。pair_id是个小技巧：正向弧 i 和反向弧 i+76 共享同一个 pair_id，后面要统计「某条路的双向总流量」或者画路段带宽图时，直接按 pair_id 分组就行，不用再拿 (from, to) 反查。

需要提醒的是，这种「rename 一下就算反向弧」的做法只对对称路网成立。SiouxFalls 双向同规格，所以没问题；如果你的数据里有单向匝道、潮汐车道，或者两个方向容量不同，就必须在原始表里显式写两行，一行正向一行反向，再走同样的展开逻辑。我见过有人拿一份非对称路网硬套这段代码，结果反向弧容量凭空多了一倍，最后 V/C 全都不对，查了两天才发现是数据层的问题。

2.3 BPR 阻抗函数与它的导数

路段走行时间用标准 BPR 函数：

t_a(x_a) = t0_a × (1 + α × (x_a / c_a)^β)

对流量求导：

dt_a/dx_a = t0_a × α × β × (x_a / c_a)^(β-1) / c_a

def bpr(t0, cap, flow, alpha=0.15, beta=4.0): """BPR 阻抗，全部按数组运算，flow 可以是标量也可以是向量""" return t0 * (1.0 + alpha * np.power(flow / cap, beta)) def bpr_grad(t0, cap, flow, alpha=0.15, beta=4.0): """dt/dx，线搜索里判断单调性要用；flow=0 且 beta>1 时结果为 0，安全""" return t0 * alpha * beta * np.power(flow / cap, beta - 1) / cap

参数取 α=0.15、β=4 是 BPR 的经典取值，也是 SiouxFalls 文献里最常用的组合。这两个数不是随便定的，它们决定了阻抗曲线的形状：β=4 意味着流量超过容量之后阻抗会陡峭上升，V/C 到 1.2 时阻抗已经是自由流的 1.15×1.2⁴≈2.4 倍。这个陡峭程度直接影响了 Frank-Wolfe 后期的收敛速度——越陡，线搜索给出的步长越小，迭代次数越多。

bpr_grad里有个边界要留意：当 β<1 时，x=0 处 (x/c)^(β-1) 会发散，得加个 eps 兜住。β=4 的场景下不会遇到，但如果你换成 β=0.5 的简化形式，记得把np.power(flow/cap, beta-1)改成np.power(np.maximum(flow/cap, 1e-12), beta-1)。

2.4 先跑通四个检查再往下写

数据读进来之后别急着写算法，先花五分钟把这四个检查跑一遍，后面能省掉大量排查时间。

def sanity_check(arcs, od, n_nodes=24): # 1) 弧数必须是路段数的两倍 assert len(arcs) == 2 * (len(arcs) // 2), '弧展开异常' # 2) 总需求核对，SiouxFalls 标准值为 360600 total = od['demand'].sum() print(f'OD 对 {len(od)} 个，总需求 {total:.0f}') # 3) 弧端点必须在 1..n_nodes 范围内 assert arcs['from_node'].between(1, n_nodes).all() assert arcs['to_node'].between(1, n_nodes).all() # 4) 单位量级抽查：t0 应为分钟级，cap 应为千级 veh/h print(f"t0 范围 {arcs['t0'].min():.2f}~{arcs['t0'].max():.2f} min") print(f"cap 范围 {arcs['cap'].min():.0f}~{arcs['cap'].max():.0f} veh/h")

四个检查分别拦四类问题：弧数不对说明展开逻辑写错了；总需求对不上说明 OD 表读漏了行或者有重复行；端点越界说明节点编号是 0-based 而不是 1-based；量级不对说明单位混了——比如 t0 给的是秒、demand 给的是日需求，直接跑下去 V/C 能爆到几十，而算法本身不会报任何错，只是给你一个看起来很正常的错结果。

提示：总需求 360600 是 SiouxFalls 高峰小时口径下的标准值，如果你的数据算出来是它的 24 倍或者 1/24，八成是把日需求和小时需求搞混了。

3. 最短路与全有全无分配：528 个 OD 对怎么叠成一条弧流量向量

3.1 邻接表：不要在这层用图库

Frank-Wolfe 每一轮迭代要做一次全有全无分配，也就是对全部 528 个 OD 对各跑一次最短路。收敛到 1e-4 一般要几十到一百多次迭代，算下来是几万次最短路调用。这个量级下，图对象本身的构造和遍历开销会占大头，所以我一般不用通用图库，直接上纯 list 加 heapq。24 个节点、152 条弧的规模，单次 Dijkstra 在几十微秒量级，整个分配跑完不到一秒。

def build_graph(arcs, n_nodes=24): """邻接表：adj[u] = [(v, arc_id, t0), ...]，1-based 节点下标""" adj = [[] for _ in range(n_nodes + 1)] for r in arcs.itertuples(): adj[int(r.from_node)].append((int(r.to_node), int(r.arc_id), float(r.t0))) # 按 arc_id 排序，保证同一节点的出弧顺序固定，结果可复现 for u in range(n_nodes + 1): adj[u].sort(key=lambda x: x[1]) return adj

邻接表长度取n_nodes + 1，是为了让节点编号直接当下标用，省掉一次 -1 的偏移转换。别小看这个细节，分配代码里有大量arcs.at[a, 'from_node']这样的取值，如果下标体系混着 0-based 和 1-based，调试的时候会非常难受。

最后那句排序不是可有可无的。Dijkstra 在遇到等成本路径时，谁先进堆谁就赢，出弧顺序变了，全有全无的结果就会变，虽然最终均衡解一样，但中间迭代的数值会对不上，做结果复现和回归测试时很难受。

3.2 Dijkstra：松弛条件里的容差是调出来的

import heapq def dijkstra(adj, origin, cost): """cost: 当前每条弧的阻抗数组，长度等于弧数 返回 dist（各节点最短成本）和 prev_arc（各节点的前驱弧编号）""" n = len(adj) - 1 INF = float('inf') dist = [INF] * (n + 1) prev_arc = [-1] * (n + 1) dist[origin] = 0.0 pq = [(0.0, origin)] while pq: d, u = heapq.heappop(pq) if d > dist[u] + 1e-15: # 惰性删除：堆里的旧标签直接丢弃 continue for v, a, _ in adj[u]: nd = d + cost[a] # 严格小于才更新，容差 1e-12 用来吸收浮点噪声 if nd < dist[v] - 1e-12: dist[v] = nd prev_arc[v] = a heapq.heappush(pq, (nd, v)) return dist, prev_arc

三个地方值得说。第一，if d > dist[u] + 1e-15: continue是标准的惰性删除写法，堆里可能同时存在同一个节点的多个旧标签，弹出时发现比当前记录差就直接跳过，不用实现 decrease-key。第二，松弛条件用nd < dist[v] - 1e-12而不是nd < dist[v]，是为了避免两条成本几乎相同的路径反复互相覆盖、把堆撑大。第三，prev_arc记录的是「到达该节点的最后一条弧」，回溯路径时从终点一路沿前驱弧往回走，走到起点为止。

容差 1e-12 不是拍脑袋定的。SiouxFalls 的 t0 量级在 1~10 分钟，路径总成本在几十的量级，双精度浮点的相对误差在 1e-16 左右，累加十几条弧之后误差放大到 1e-15 量级，取 1e-12 有足够余量又不会把真正不同的路径误判成平局。换成更大的网络、成本量级到 1e4 以上时，这个容差要跟着放大到 1e-10，否则会因为舍入误差被判成「不更新」，导致实际走的不是最短路。

3.3 全有全无分配：一次遍历叠完 528 个 OD 对

def all_or_nothing(arcs, adj, od, cost): """在当前阻抗下，把每个 OD 对的需求全部压到它的最短路上，返回弧流量向量""" y = np.zeros(len(arcs), dtype=float) for r in od.itertuples(): o, d, q = int(r.o), int(r.d), float(r.demand) if q <= 0: continue dist, prev_arc = dijkstra(adj, o, cost) if dist[d] == float('inf'): raise ValueError(f'OD {o}->{d} 不可达，检查弧方向是否建全') v = d while v != o: # 沿前驱弧回溯到起点 a = prev_arc[v] y[a] += q v = int(arcs.at[a, 'from_node']) return y

外层循环是 528 个 OD 对，内层是 Dijkstra 加一次回溯。注意这里产出的是弧流量 y，不是路径流量——Frank-Wolfe 的方向搜索只需要知道「每条弧上应该有多少流量」，不需要保存具体路径。这个设计让内存占用维持在 O(E) 而不是 O(路径数 × 弧数)。

arcs.at[a, 'from_node']能直接取到值，前提是 arcs 的 DataFrame 索引刚好等于 arc_id（0 到 151 连续）。如果中途做过dropna或者sort_values打乱过索引，这里会静默取到错的行，流量加到错误的弧上——不会报错，只会让结果慢慢偏掉。稳妥做法是在load_network返回前加一句arcs = arcs.reset_index(drop=True)，把索引钉死。

累加顺序不影响最终结果，浮点加法的舍入误差在 1e-10 量级，对 gap 的影响可以忽略。但如果要跟别人对数，建议按 OD 表原始顺序遍历，别用groupby之类会改变顺序的操作。

3.4 平局路径与结果可复现性

SiouxFalls 里存在不少等成本路径，尤其是在自由流阻抗下，很多 OD 对的第二条、第三条最短路和第一条几乎同价。这意味着全有全无的解本身不唯一：同一个网络、同一个阻抗，两次运行可能把某个 OD 的需求压到不同的弧上。

这不是 bug。Frank-Wolfe 的每一步方向都是可行域的顶点，选哪个顶点只影响路径，不影响收敛性，因为下一步的线搜索会把流量重新按阻抗分配回去。真正需要担心的是结果不可复现：如果某次跑出来 gap 到了 1e-5，另一次同一个脚本跑出来卡在 3e-4，先别怀疑算法，先看两次的出弧顺序是不是一样。

固定顺序的办法在 3.1 里已经做了：建图时按 arc_id 排序。再配合固定的浮点容差，同一份数据、同一个 Python 版本，两次运行的结果能对到小数点后十位。跨机器跨版本时浮点库可能有极微小差异，但 gap 曲线形状应该完全一致。

还有一种情况要区分：如果两条路径成本只差 1e-13，代码把它们判成平局，实际流量分配差异可能到几千 veh/h。这种量级的流量差异会传导到 BPR 函数上，让下一步的阻抗出现可见差别。遇到这种现象，把 Dijkstra 的容差从 1e-12 收紧到 1e-14 试试，如果结果还是飘，那说明网络里确实存在严格等价的路径，属于正常现象，看 gap 和总阻抗就行，不用纠结单条弧的具体数值。

4. Frank-Wolfe 主循环：用二分线搜索把相对间隙压到 1e-4

4.1 四个变量、一个方向、一次线搜索

UE 问题的数学形式是 Beckmann 目标函数：

min Z(x) = Σ_a ∫₀^x_a t_a(w) dw s.t. 流量守恒约束，x ≥ 0

Frank-Wolfe 是求解这个凸规划最经典的一阶方法，每一轮干四件事：算当前阻抗、找下降方向、定步长、更新流量。四个变量对应四行代码，但要跑得稳，每一行都有讲究。

先把主循环摆出来，再逐个拆。

def frank_wolfe(arcs, adj, od, max_iter=500, gap_tol=1e-4, alpha=0.15, beta=4.0, verbose=True): t0 = arcs['t0'].to_numpy() cap = arcs['cap'].to_numpy() # 1) 初始化：自由流阻抗下的全有全无，得到初始可行解 x0 x = all_or_nothing(arcs, adj, od, t0.copy()) hist = [] for k in range(max_iter): # 2) 用当前流量更新弧阻抗 t = bpr(t0, cap, x, alpha, beta) # 3) 方向搜索：在当前阻抗下再做一次全有全无 y = all_or_nothing(arcs, adj, od, t) # 4) 收敛判据：相对间隙 tstt_x = float(np.dot(x, t)) tstt_y = float(np.dot(y, t)) gap = (tstt_x - tstt_y) / tstt_x hist.append((k, tstt_x, gap)) if verbose and k % 10 == 0: print(f'iter {k:3d} TSTT {tstt_x:14.2f} gap {gap:.3e}') if gap <= gap_tol: break # 5) 线搜索定步长，6) 凸组合更新流量 lam = line_search(x, y, t0, cap, alpha, beta) x = x + lam * (y - x) return x, np.array(hist, dtype=float)

初始化那一步用的是自由流阻抗 t0，不是零流量下的阻抗——两者的区别在于 BPR 在 x=0 时阻抗恰好等于 t0，所以数值上一样，但语义上要说清楚：初始解是「所有车都走自由流最短路径」，这是个极端拥堵的解，正好作为迭代起点。

4.2 方向搜索：为什么 y 永远落在可行域的顶点上

Frank-Wolfe 的子问题是把非线性的 Z(x) 在当前点做一阶近似，变成线性规划：

min ∇Z(x)ᵀ y，s.t. y ∈ 可行域

而 ∇Z(x) 的第 a 个分量恰好就是 t_a(x_a)（Beckmann 函数求导的结果）。所以子问题变成「在固定阻抗下最小化总走行时间」，这正是全有全无分配要干的事。

关键在于这个线性规划的最优解一定在可行域多面体的顶点上，而交通分配的可行域顶点有明确的物理含义：每个 OD 对的需求全部集中在一条路径上。所以 y 的每一维要么是 0，要么是若干 OD 对需求的和，绝不会有「一半走这条路一半走那条」的中间状态。

这个性质带来两个直接后果。好的一面是，算 y 不需要解大规模线性规划，跑一次全有全无就精确得到顶点解，这是 Frank-Wolfe 在交通分配里长盛不衰的核心原因。不好的一面是，顶点解非常「极端」，每次迭代只能在顶点之间跳跃，越接近最优解，两个顶点之间的连线越难贴合真正的最优方向，这就是收敛速度在后期急剧下降的根源。工程上表现为：前 20 次迭代 gap 从 1 掉到 1e-2，后 80 次才从 1e-2 磨到 1e-4。

4.3 二分法线搜索：dZ/dλ 单调递增这件事救了一切

定步长这步，目标函数沿方向 d = y - x 是一元函数：

φ(λ) = Z(x + λd)，λ ∈ [0, 1]

求导：

φ'(λ) = Σ_a d_a × t_a(x_a + λ d_a)

因为每条弧的 t_a(·) 是单调递增的，而 d_a 是常数，所以 φ'(λ) 是若干个单调递增函数乘常数再求和——也是单调递增的。单调递增意味着一分法能解：找到 φ'(λ)=0 的那个 λ 就行。

这个性质是 BPR 函数的「单调性」直接给的，也是 Frank-Wolfe 能配精确线搜索的前提。如果阻抗函数不是单调的（比如某些考虑排队的模型会出现非单调段），二分法就会失效，得换成黄金分割或者直接扫描。

def line_search(x, y, t0, cap, alpha=0.15, beta=4.0, tol=1e-10, max_iter=60): """二分求 φ'(λ)=0；φ' 单调递增，所以可以用二分""" d = y - x def dphi(lam): f = x + lam * d # 试探点处的弧流量 return float(np.dot(d, bpr(t0, cap, f, alpha, beta))) if dphi(0.0) >= 0.0: return 0.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)

函数开头的两个边界判断不是可有可无的防御性代码，是必须的。dphi(0) >= 0意味着沿着方向走一步目标函数就上升，说明当前解已经是最优（此时 gap 会同步给出极小的值，两个判据互相印证）。dphi(1) <= 0意味着走到方向端点目标函数还在下降，直接取 λ=1。这两个分支如果不写，二分会在 [0,1] 内部收敛到一个假的根，步长算出来是错的。

二分迭代 60 次能把区间压到 2⁻⁶⁰ ≈ 1e-18，远超浮点精度，实际上 40 次左右就够了。tol=1e-10是给一个提前退出的口子，避免每次都跑满。

4.4 relative gap 的正确写法，以及三种常见错法

收敛判据用的是相对间隙：

gap = (Σ_a x_a t_a(x) - Σ_a y_a t_a(x)) / Σ_a x_a t_a(x)

分子是「当前解的总阻抗」减去「线性子问题下界」，分母用当前解的总阻抗做归一化。这个量恒大于等于 0（因为 y 是子问题最优解），且随迭代单调下降，是所有交通分配软件里的标准输出。

第一种错法：分母用 Σ y_a t_a。这样算出来的 gap 会比正确值小一到两个数量级，看起来收敛很快，实际上提前停在了 1e-2 的水平。第二种错法：用相邻两次迭代的 TSTT 变化率 (TSTT_k - TSTT_{k-1}) / TSTT_k 当判据。这个量在后期确实很小，但流量还在弧之间大幅漂移——总阻抗趋于稳定不代表每条弧的流量稳定，因为阻抗曲线在均衡点附近是平的。第三种错法：用绝对间隙不除分母。换个需求量大的网络，绝对间隙的量级跟着变，同一个阈值在不同网络上代表完全不同的收敛程度。

注意：gap 的符号也要看一眼。如果算出来是负数，别怀疑算法，先查 y 是不是真的在当前阻抗下做的全有全无——最常见的原因是算 gap 用的 t 和算 y 用的 t 不是同一个数组，中间被覆盖过。

4.5 拼起来跑一次

arcs, od = load_network('siouxfalls_net.csv', 'siouxfalls_trips.csv') sanity_check(arcs, od, n_nodes=24) adj = build_graph(arcs, n_nodes=24) x, hist = frank_wolfe(arcs, adj, od, max_iter=500, gap_tol=1e-4) print(f'迭代 {int(hist[-1, 0])} 次，最终 gap {hist[-1, 2]:.3e}') # 把结果挂回弧表，方便看 V/C arcs['flow'] = x arcs['voc'] = arcs['flow'] / arcs['cap'] arcs['time'] = bpr(arcs['t0'].to_numpy(), arcs['cap'].to_numpy(), x) print(arcs[['from_node', 'to_node', 'flow', 'voc', 'time']].head(10))

在 SiouxFalls 上，gap 从初始的 0.5 左右压到 1e-4 通常需要几十到一百多次迭代，具体次数取决于容差设定和 BPR 参数的陡峭程度。如果跑到 300 次还没下去，问题不在算法，在第 5 章列的那几个坑里。

5. 避坑与排查：5 个让 gap 卡住的翻车现场

这一章全部是血泪经验，按「现象 → 原因 → 解决」写，遇到对应症状直接对号入座。

5.1 回程流量凭空消失：弧只建了一半

现象：总需求 360600，但把弧流量加起来只有 18 万左右；部分路段的 V/C 高得离谱，另一些路段几乎没流量；OD 表里某些方向压根找不到路径，Dijkstra 抛「不可达」。

原因：路段表按无向边组织，一条边只写了一行。代码直接把它当成有向弧用，只保留了 from→to 一个方向。SiouxFalls 的 76 条路段里有很多是主干道，只建单向之后，反向的 OD 对无路可走，全有全无分配时要么报错，要么把需求挤到绕行路径上。

解决：老老实实做正反两条弧展开，也就是 2.2 那段代码。展开后弧数应该是 152。加一条断言assert len(arcs) == 2 * n_edges，跑之前就拦住。判断一个网络有没有展开正确，最快的办法是看「每个节点的出度是否非零」——如果某个节点只有入弧没有出弧，十有八九就是方向没建全。

5.2 用 1/n 步长替代线搜索：gap 卡在 1e-2

现象：前 20 次迭代 gap 掉得飞快，从 0.5 降到 0.01，然后就开始原地打转，跑到 1000 次也就到 5e-3，怎么都不下去。换成精确线搜索之后，同样的网络 80 次就到了 1e-4。

原因：把 λ = 1/(k+1) 当成步长，也就是常说的 MSA（Method of Successive Averages）。MSA 确实保证收敛，但收敛速率是 O(1/√k) 量级，而且步长不可和——后期每次只移动一点点，流量分布已经接近均衡了，但 gap 的下降速度被步长拖死。很多人图省事用 MSA，结果误以为 Frank-Wolfe 只能到这个精度。

解决：换成 4.3 的二分线搜索。如果因为阻抗函数不单调而没法用二分，退一步用黄金分割或者三分搜索也行，效果比 MSA 好一个数量级。真要用 MSA，至少把判据放宽到 1e-3，别指望它跑到 1e-4。

5.3 直接覆盖流量：把凸组合写成 x = y

现象：迭代过程正常，每一次的 TSTT 都在下降，但弧流量在两次迭代之间大幅震荡，某条路的流量这次是 8000、下次是 200，再下次又回到 7000。gap 曲线呈锯齿状。

原因：把更新写成了x = y，也就是每一轮直接把方向解当成新解。这在数学上退化成反复做全有全无分配，等价于在顶点之间无限跳，永远不收敛到内点。全有全无是「所有车走同一条路」，而均衡解通常是「车分散在多条路径上」，两者不在一个空间里。

解决：凸组合x = x + lam * (y - x)，λ 必须来自线搜索且在 [0,1] 内。可以在更新后加一句检查：assert np.all(x >= -1e-9)和assert lam >= 0 and lam <= 1，一旦 λ 越界立刻抓出来。

5.4 gap 分母写错：收敛判据永远"达标"

现象：第 0 次迭代打印出来的 gap 就是 1e-5，直接跳出循环，输出的结果看起来像自由流状态，V/C 全都很低，但稍微验算就发现流量守恒根本不满足。

原因：分母用了 Σ y_a t_a，而 y 是子问题的解，Σ y t 天然小于等于 Σ x t，这个「gap」在某些网络里会算出负数或者极小值。另一种情况是把 gap 定义成 (TSTT_x - TSTT_y) / TSTT_y，同样是错的。

解决：严格按 (Σ x t - Σ y t) / Σ x t 写，并且加两条断言：gap 必须 ≥ 0；初始 gap 在 0.3~0.8 这个量级（SiouxFalls 自由流初始解对应的 gap 大约在这个范围）。如果初始 gap 不在这个区间，先别往下跑，回去查数据。

5.5 单位没对齐：V/C 全部爆表

现象：算法正常收敛，gap 到了 1e-6，但输出的 V/C 里有一半路段超过 3，走行时间是最小值的十几倍，总阻抗算出来和文献里的量级差一个数量级。

原因：t0、cap、demand 三个量的时间粒度不一致。常见组合是 demand 用的是日需求（SiouxFalls 的日需求约 360600×24 量级），而 cap 用的是 veh/h。需求直接放大 24 倍，容量没变，V/C 自然全部爆掉。另一种是 t0 用秒、长度用英里，算法照样收敛，但输出的时间全都差 60 倍。

解决：跑之前打印三个量的量级做交叉验证。t0 应该在 1~10 分钟，cap 在几千 veh/h，单个 OD 对的需求在几百到几千 veh/h，总需求 360600。三个数只要有一个量级不对，立刻停下来查。另外，收敛之后再做一个流量守恒检查（第 6 章会写），守恒不满足说明数据层有问题，跟算法无关。

6. 收敛验证与调参：怎么证明你算出来的是一份均衡解

gap 小不等于结果对。gap 只能说明「当前解和它的线性近似下界的距离小」，如果数据层或者阻抗函数有错，gap 照样能收敛到 1e-6，只是收敛到一个错答案。所以我每次跑完都会固定做三项体检。

第一项是流量守恒。对每个节点，流出减流入应该等于该节点的发生量减吸引量，残差应该在 1e-6 以内。

def check_conservation(arcs, od, x, n_nodes=24): gen = od.groupby('o')['demand'].sum() attr = od.groupby('d')['demand'].sum() inflow = np.zeros(n_nodes + 1) outflow = np.zeros(n_nodes + 1) for a in range(len(arcs)): outflow[int(arcs.at[a, 'from_node'])] += x[a] inflow[int(arcs.at[a, 'to_node'])] += x[a] worst = 0.0 for n in range(1, n_nodes + 1): # 守恒：流出 - 流入 = 发生 - 吸引 resid = outflow[n] - inflow[n] - (gen.get(n, 0.0) - attr.get(n, 0.0)) worst = max(worst, abs(resid)) print(f'最大守恒残差 {worst:.3e}') return worst # 正常应该 < 1e-6

第二项是 gap 曲线的形状。在 1e-4 这个量级上，Frank-Wolfe 的 gap 应该是一条接近直线的下降段（半对数坐标下），如果曲线出现明显的台阶或者反弹，说明某次线搜索出了问题——最常见的是 dphi 的边界判断被绕过，λ 取到了 0 之外的值。

第三项是用 MSA 做交叉验证。用同一份数据跑 500 次 MSA，它收敛得慢，但收敛方向是对的。FW 的最终总阻抗应该 ≤ MSA 的总阻抗，而且两者的差距随迭代次数增加而缩小。如果 FW 的结果比 MSA 还差，那一定是更新公式或者判据写错了，这是最省时间的自检手段。

调参这块，真正需要动的只有下面几个：

还有一个小技巧值得单独说：如果只是想看路网拥堵的空间分布，跑到 gap 1e-3 就够了，迭代次数能省掉三分之二；但如果要拿这份结果去做

本文还有配套的精品资源，点击获取