ESPPRC 定价算法：从基础标号到工程级加速

作者：Sebastilan & Claude（AI 协作）

---

系列定位：BPC 学习系列第四篇。第三篇从零实现了完整的 BPC 求解器，并在小实例上验证了正确性。本篇深入定价子问题（ESPPRC），通过6 层技术逐步叠加 + 消融实验量化每层贡献，最终与开源 SOTA（cspy、RouteOpt 2.0）对标，给出工程级推荐配置。

摘要：同一个定价子问题，用 Gurobi 13（i7-12700）求解需要 10-120 秒（$n=12$ 时 10-101 秒，$n\ge 18$ 超时找不到最优解）；换成基础标号算法，$n=12$ 只需 1ms（快 10000 倍），但 $n\ge 18$ 同样超时；再叠加双向搜索 + Bucket + 并行，$n=18$ 从超时降到 1.4 秒；最后加上 ng-route 松弛，$n=18$ 仅 10ms，$n=31$ 仅 8ms，$n=50$ 以上仍毫秒级。5 个数量级的差距，6 层技术叠加的结果。本文不涉及理论创新，纯粹是现有最佳技术的拆解、组合、实测与对比——每层都有独立消融实验数据，最终比学术主流库 cspy 快 100-31000 倍。

---

1. 引言：定价为什么是 BPC 的瓶颈

在列生成的每轮迭代中，求解 RMP（限制主问题）耗时通常只有毫秒级——它不过是一个规模有限的 LP。真正的时间黑洞在定价子问题：给定 RMP 的对偶值 $\pi$，在所有可行路线中找 reduced cost 最小的那条。

对于中等规模实例（$n=50$ 客户），定价耗时可以占整个列生成过程的 95% 以上。B&B 树中每个节点都要跑完整的列生成循环，这意味着定价函数会被调用数万次。定价的速度直接决定了 BPC 能解多大规模的实例。

定价慢的根本原因：ESPPRC 是 NP-hard 的。精确求解的状态空间为 $O(n\cdot 2^n)$——12 个客户约 5 万状态，18 个客户约 500 万，21 个客户已经超过 3 亿。朴素实现在 $n=18$ 时就会超时（60 秒以上）。

本文的目标是系统回答：如何把 ESPPRC 从"$n=18$ 超时"推进到"$n=50+$ 实用"？技术路线：

层次	策略	效果	出处
算法层	双向标号（Bidir）	TIMEOUT → 25.6s（n=18）	Righini & Salani (2006)
数据结构层	Bucket + RC 排序	25.6s → 2.4s（10.7x）	Sadykov et al. (2021)
并行层	双线程 fwd/bwd	2.4s → 1.4s（1.68x）	工程优化
松弛层	ng-route 松弛	精确→松弛，n=31 从不可行到 8ms	Baldacci et al. (2011)
精确化层	DSSR 迭代收紧	松弛恢复为精确，n≤31 秒级	Righini & Salani (2008)
工程层	标号池预分配	额外 ~1.2x	工程优化

本文结构：第 2 章定义问题与复杂度框架；第 3 章实现基础标号算法（Level 0 基线）；第 4 章对精确型做逐层消融实验；第 5 章引入 ng-route 松弛突破指数壁垒；第 6 章与开源 SOTA 对标；第 7 章给出工程推荐配置。

---

2. 问题定义与复杂度框架

2.1 符号表

符号	含义
$G=(V,A)$	完全有向图， V = { 0 , 1 , … , n } V = \{0, 1, \ldots, n\} V={0,1,…,n}，0 为 depot
$c_{ij}$	弧 $(i,j)$ 的行驶成本（整数，满足三角不等式）
$d_i$	客户 $i$ 的需求量（$d_0=0$）
$Q$	车辆容量约束
${\pi }_i$	RMP 覆盖约束 $i$ 的对偶变量（来自当前列生成迭代）
${\bar{c}}_{ij}$	修正成本：${\bar{c}}_{ij}=c_{ij}-{\pi }_j$（可为负）
${\bar{c}}_r$	路线 $r$ 的 reduced cost：${\bar{c}}_r=c_r-{\sum }_{i\in r}{\pi }_i$
$N=n+1$	节点总数（含 depot）
$\Delta$	ng-route 邻域大小（默认 8）

2.2 ESPPRC 形式化

ESPPRC 本质上是一个最短路径问题，但有三个关键区别使它从多项式跳到 NP-hard：

1. 弧代价可以是负的：修正成本 ${\bar{c}}_{ij}=c_{ij}-{\pi }_j$，当对偶值 ${\pi }_j$ 较大时为负。Dijkstra 要求非负权重，Bellman-Ford 能处理负权但不能处理负环——而 ESPPRC 的 reduced cost 图恰恰充满负环（详见 2.5.2 节）。
2. 有容量约束：路径上所有客户的总需求 $\sum d_i\le Q$，为状态空间增加一个资源维度。
3. 路径必须 elementary：每个客户至多访问一次——这才是复杂度爆炸的根源。

去掉第 3 条，状态只需 $(node,load)$，$O(n^{2}Q)$ 伪多项式可解。加上第 3 条，状态变成 $(node,load,visited\_set)$，$visited\_set$ 有 $2^n$ 种可能——NP-hard。

min ⁡ ∑ ( i , j ) ∈ A c ˉ i j x i j (1) 最小化 reduced cost s.t. ∑ j ∈ V x 0 j = 1 (2) 从 depot 出发 ∑ i ∈ V x i j = ∑ k ∈ V x j k ∀   j ∈ V ∖ { 0 } (3) 流量守恒 ∑ j ∈ V x j 0 = 1 (4) 返回 depot ∑ ( i , j ) ∈ A d j x i j ≤ Q (5) 容量约束 ∑ j ∈ V x i j ≤ 1 ∀   i ∈ V ∖ { 0 } (6) Elementary：每客户至多访问一次 x i j ∈ { 0 , 1 } ∀   ( i , j ) ∈ A (7) 决策变量 \begin{array}{llll} \min & \displaystyle\sum_{(i,j) \in A} \bar{c}_{ij} x_{ij} & & \text{(1) 最小化 reduced cost} \\[0.8em] \text{s.t.} & \displaystyle\sum_{j \in V} x_{0j} = 1 & & \text{(2) 从 depot 出发} \\[0.8em] & \displaystyle\sum_{i \in V} x_{ij} = \sum_{k \in V} x_{jk} & \forall\, j \in V \setminus \{0\} & \text{(3) 流量守恒} \\[0.8em] & \displaystyle\sum_{j \in V} x_{j0} = 1 & & \text{(4) 返回 depot} \\[0.8em] & \displaystyle\sum_{(i,j) \in A} d_j x_{ij} \leq Q & & \text{(5) 容量约束} \\[0.8em] & \displaystyle\sum_{j \in V} x_{ij} \leq 1 & \forall\, i \in V \setminus \{0\} & \text{(6) Elementary：每客户至多访问一次} \\[0.8em] & x_{ij} \in \{0,1\} & \forall\, (i,j) \in A & \text{(7) 决策变量} \end{array} mins.t.​(i,j)∈A∑​cˉij​xij​j∈V∑​x0j​=1i∈V∑​xij​=k∈V∑​xjk​j∈V∑​xj0​=1(i,j)∈A∑​dj​xij​≤Qj∈V∑​xij​≤1xij​∈{0,1}​∀j∈V∖{0}∀i∈V∖{0}∀(i,j)∈A​(1) 最小化 reduced cost(2) 从 depot 出发(3) 流量守恒(4) 返回 depot(5) 容量约束(6) Elementary：每客户至多访问一次(7) 决策变量​

复杂度根源：约束 (6) 是 Elementary 约束——每个客户至多访问一次。这使得 ESPPRC 本质上是带资源约束的 Hamiltonian Path 问题，已知 NP-hard。如果去掉 (6)，问题退化为普通的 SPPRC（Shortest Path with Resource Constraints），可用伪多项式时间 DP 求解。

2.3 复杂度分解框架

实现 ESPPRC 的标准方法是标号算法（Label-Setting DP）。定价总耗时可分解为：

$$T_{pricing}=\frac{|Labels|\times T_{per\_label}}{P}$$

因子	含义	决定因素
$|Labels|$	生成的标号总数	状态空间大小 + 支配剪枝效率
$T_{per\_label}$	每个标号的平均处理时间	支配检查代价 + 扩展代价 + 内存操作
$P$	并行度	线程数（fwd/bwd 独立）

每种优化策略只打击其中一个因子：

策略	打击因子	理论效果
Bidir（双向）	$|Labels|$：指数深度减半	$2^n\to 2\cdot 2^{n/2}$
Bucket + RC排序	$T_{per\_label}$：支配检查加速	支配检查次数 /16x
Parallel	$P$：双线程	墙钟时间 /2
ng-route松弛	$|Labels|$：状态空间宽度	$2^n\to 2^{\Delta }$
Pool预分配	$T_{per\_label}$：内存操作	实测无显著效果

关键设计原则：这 4 个因子相互独立，策略可以乘性叠加，无交互效应（第 4 章实验验证）。

---

2.5 能不能直接用求解器解？

在介绍标号算法之前，一个自然的问题：ESPPRC 能写成 MIP，为什么不直接交给 Gurobi 求解？

答案是：可以，但远不够用。本章通过实验数据说明 MIP 的极限，并逐步引出 ng-route 松弛的必要性。

2.5.1 ESPPRC 的 MIP 建模

使用 MTZ（Miller-Tucker-Zemlin）子环消除，可以将 ESPPRC 紧凑地建模为 MIP。

决策变量：

变量	类型	含义
x i j ∈ { 0 , 1 } x_{ij} \in \{0,1\} xij​∈{0,1}	二值	是否使用弧 $(i,j)$
$u_i\in [0,Q]$	连续	到达客户 $i$ 时的累计载荷

MTZ 模型：

min ⁡ ∑ i ∈ V ∑ j ∈ V , j ≠ i c ˉ i j x i j (1) 最小化 reduced cost s.t. ∑ j = 1 n x 0 j = 1 (2) 从 depot 出发一次 ∑ i = 1 n x i 0 = 1 (3) 返回 depot 一次 ∑ i ∈ V , i ≠ j x i j = ∑ k ∈ V , k ≠ j x j k ∀   j ∈ { 1 , … , n } (4) 流量守恒 ∑ i ∈ V , i ≠ j x i j ≤ 1 ∀   j ∈ { 1 , … , n } (5) 每客户至多访问一次 u j ≥ u i + d j − ( Q + d j ) ( 1 − x i j ) ∀   i , j ∈ { 1 , … , n } , i ≠ j (6) MTZ 子环消除 + 载荷传播 d j ⋅ ∑ i x i j ≤ u j ≤ Q ⋅ ∑ i x i j ∀   j ∈ { 1 , … , n } (7) 载荷上下界 x i j ∈ { 0 , 1 } ,   u i ≥ 0 (8) 变量域 \begin{array}{llll} \min & \displaystyle\sum_{i \in V} \sum_{j \in V,j \neq i} \bar{c}_{ij} x_{ij} & & \text{(1) 最小化 reduced cost} \\[0.8em] \text{s.t.} & \displaystyle\sum_{j=1}^{n} x_{0j} = 1 & & \text{(2) 从 depot 出发一次} \\[0.8em] & \displaystyle\sum_{i=1}^{n} x_{i0} = 1 & & \text{(3) 返回 depot 一次} \\[0.8em] & \displaystyle\sum_{i \in V, i \neq j} x_{ij} = \sum_{k \in V, k \neq j} x_{jk} & \forall\, j \in \{1,\ldots,n\} & \text{(4) 流量守恒} \\[0.8em] & \displaystyle\sum_{i \in V, i \neq j} x_{ij} \leq 1 & \forall\, j \in \{1,\ldots,n\} & \text{(5) 每客户至多访问一次} \\[0.8em] & u_j \geq u_i + d_j - (Q + d_j)(1 - x_{ij}) & \forall\, i,j \in \{1,\ldots,n\}, i \neq j & \text{(6) MTZ 子环消除 + 载荷传播} \\[0.8em] & d_j \cdot \displaystyle\sum_{i} x_{ij} \leq u_j \leq Q \cdot \displaystyle\sum_{i} x_{ij} & \forall\, j \in \{1,\ldots,n\} & \text{(7) 载荷上下界} \\[0.8em] & x_{ij} \in \{0,1\},\ u_i \geq 0 & & \text{(8) 变量域} \end{array} mins.t.​i∈V∑​j∈V,j=i∑​cˉij​xij​j=1∑n​x0j​=1i=1∑n​xi0​=1i∈V,i=j∑​xij​=k∈V,k=j∑​xjk​i∈V,i=j∑​xij​≤1uj​≥ui​+dj​−(Q+dj​)(1−xij​)dj​⋅i∑​xij​≤uj​≤Q⋅i∑​xij​xij​∈{0,1}, ui​≥0​∀j∈{1,…,n}∀j∈{1,…,n}∀i,j∈{1,…,n},i=j∀j∈{1,…,n}​(1) 最小化 reduced cost(2) 从 depot 出发一次(3) 返回 depot 一次(4) 流量守恒(5) 每客户至多访问一次(6) MTZ 子环消除 + 载荷传播(7) 载荷上下界(8) 变量域​

约束 (6) 是核心：当 $x_{ij}=1$ 时，强制 $u_j\ge u_i+d_j$（载荷单调递增，消除子环）；当 $x_{ij}=0$ 时，$(Q+d_j)$ 的大 M 使约束自动满足。完整 Python 实现见 components/espprc/testset/espprc_oracle.py（mip_oracle 函数）。

MIP 与标号算法的实验对比（i7-12700，超时 120 秒）：

实例	$n$	MIP（Gurobi 13）	状态	精确标号（C8）	ng-BIDIR_POOL（$\Delta =8$）
tiny3	3	0.005s	optimal	<0.1ms	<0.1ms
tiny5	5	0.020s	optimal	<0.1ms	<0.1ms
small8	8	0.403s	optimal	<0.1ms	<0.1ms
en13k4_converged	12	15.7s	optimal	1ms	0.35ms
en13k4_iter1	12	101s	optimal	1ms	0.35ms
en13k4_iter2	12	10.6s	optimal	1ms	0.35ms
pn16k8_converged	15	3.3s	optimal	<1ms	0.5ms
pn16k8_iter1	15	3.4s	optimal	<1ms	0.5ms
pn19k2_initial	18	121s	TIMEOUT	1.4s	10ms
pn19k2_perturbed	18	121s	TIMEOUT	1.4s	10ms
pn22k2_initial	21	121s	TIMEOUT	TIMEOUT	9ms
pn22k2_perturbed	21	121s	TIMEOUT	TIMEOUT	9ms
an32k5_initial	31	121s	TIMEOUT	TIMEOUT	8ms

几点值得注意：

1. MIP 对对偶值极为敏感：同一实例 $n=12$，三个不同对偶点下 MIP 耗时从 10.6s 到 101s 相差近 10 倍，而标号算法始终稳定在 1ms。
2. MIP 超时后解质量差：$n=18$ 时 MIP 不仅超时，还找不到最优解（RC=-542 vs oracle=-545）。BPC 框架中定价必须返回真正的最小 RC 列，MIP 超时给出的次优解会导致 LP 无法收敛。
3. 标号算法快 4-5 个数量级：$n=12$ 的三个实例，MIP 平均约 42 秒，C8 标号算法 1ms——差距约 42000 倍。

为什么标号算法远快于 MIP？ MIP 变量数 $O(n^2)$、约束数 $O(n^2)$，看似规模不大，但 Gurobi 内部的 branch-and-bound 搜索树会指数爆炸——对偶值越"尖锐"（大 ${\pi }_i$ 使某些弧修正成本极负），分支决策越困难。标号算法虽然理论复杂度 $O(n\cdot 2^n)$ 同样指数，但通过支配规则大幅剪枝，实际生成的标号数远小于理论上界。

2.5.2 去掉 elementarity 会怎样？

既然 elementary 约束是 NP-hard 的来源，能不能直接去掉它？

**SPPRC（允许重复访问）**的复杂度从 $O(n\cdot 2^n)$ 降到 $O(n^{2}Q)$（伪多项式），标号状态退化为 $(node,load)$——去掉了 $2^n$ 的访问 bitmask。$n=50,Q=200$ 时只有约 50 万个状态，毫秒级可解。

但问题立即出现：reduced cost 图里充满负 2-cycle。

$${\bar{c}}_{jk}+{\bar{c}}_{kj}=(c_{jk}-{\pi }_k)+(c_{kj}-{\pi }_j)=c_{jk}+c_{kj}-{\pi }_j-{\pi }_k$$

CG 早期对偶值 $\pi$ 较大时，几乎每对相邻节点都满足 ${\pi }_j+{\pi }_k>c_{jk}+c_{kj}$，构成负 2-cycle。SPPRC 算法会找到 $j\to k\to j\to k\to \cdots$ 的振荡路径，RC 极度负数——但这些路径对 VRP 完全无意义（重复客户无法配送）。

实证：第 4 章的弧消除实验（Bellman-Ford）发现 arcs_eliminated = 0——修正成本图中负环如此普遍，Bellman-Ford 无法给出任何有意义的下界，无法消除任何弧。

因此，必须在某种程度上维护 elementarity，问题是"维护多少"——这一设计空间将在第 4 章精确型实验结束后、第 5 章引入 ng-route 时详细展开。ng-route 松弛不能用 MIP 建模——其原因将在第 5 章介绍 ng-route 原理后详述。

---

3. Level 0：基础标号算法

3.1 标号定义与支配规则

标号是算法状态的完整快照。在节点 $v$ 的一个标号 $L$ 包含：

Label L = {
    node:  当前节点 v
    rc:    从 depot 到 v 的累计 reduced cost
    load:  从 depot 到 v 的累计载重
    vmask: 已访问客户的 bitmask（uint64_t，bit k-1 对应客户 k）
    prev:  父标号索引（用于路径回溯）
}

支配规则：在同一节点 $v$，标号 $L_1$ 支配 $L_2$，当且仅当：

$$r{c}_1\le r{c}_2,loa{d}_1\le loa{d}_2,visite{d}_1\subseteq visite{d}_2$$

三个条件同时成立。直觉：$L_1$ 更便宜、更轻、访问的客户集合更小（未来扩展限制更少）。被支配的 $L_2$ 永远不可能产生比 $L_1$ 更好的完整路线，安全丢弃。

支配规则是整个算法的效率基础——没有它，标号数会随深度指数爆炸。

3.2 算法流程

3.3 核心实现

基础标号算法的核心逻辑约 30 行 C++，包含四个部分：修正边权计算（${\bar{c}}_{ij}=c_{ij}-{\pi }_j$）、BFS 主循环（逐标号扩展）、支配检查（全量扫描 node_labels[v]）、路径收集（回 depot 后 RC < 0）。完整实现见 GitHub 开源库 components/espprc/exact/espprc.h（ESPPRC::solve 函数）。

3.4 Level 0 性能基线

实例	$n$	MIP（Gurobi 13）	标号 Level 0	加速比
tiny3	3	0.005s	<0.1ms	~50x
small8	8	0.4s	<0.1ms	~4000x
en13k4	12	10-101s	1ms	10000-100000x
pn16k8	15	3.3s	<1ms	~6600x
pn19k2	18	TIMEOUT（120s）	TIMEOUT（60s）	均超时

即使是最朴素的标号算法，也比 MIP 求解器快 4-5 个数量级。$n=18$ 时两者都超时，但标号算法的超时上限是 60s，MIP 需要 120s 仍找不到最优解——这是两类算法的代际差异。

$n$ 每增加 3，状态空间约增 $2^3=8$ 倍，从 18 到 21 是不可逾越的鸿沟。

复杂度天花板：$O(n\cdot 2^n)$ 的精确 bitmask，$n=21$ 时状态空间 $\approx 44$ 亿。即使每个状态只需 1 纳秒，也需要 44 秒——不计算支配检查开销。

---

4. 精确型逐层消融实验

实验设置：测试平台：i9-14900K，MINGW64 GCC 15.2，Release 模式，超时 60 秒。4 个 runtime bool 开关对应 4 种优化，构成 9 种配置（$C_0$ 到 $C_8$）。

配置	Pool	Bucket	Bidir	Para	说明
C0	-	-	-	-	裸基线
C1	ON	-	-	-	单测 Pool
C2	-	ON	-	-	单测 Bucket
C3	-	-	ON	-	单测 Bidir
C4	ON	ON	-	-	现有 EXACT_POOL
C5	ON	-	ON	-	Pool+Bidir
C6	-	ON	ON	-	Bucket+Bidir
C7	ON	ON	ON	-	全部无并行
C8	ON	ON	ON	ON	满配

4.1 因子 $|Labels|$：双向搜索

回顾 2.3 节的复杂度分解：$T_{pricing}=|Labels|\times T_{per\_label}/P$。双向搜索打击的是第一个因子 $|Labels|$——通过将搜索深度从 $n$ 减半到 $n/2$，标号总数从 $O(2^n)$ 降到 $O(2\cdot 2^{n/2})$。

原理：传统正向标号从 depot 扩展到所有状态，深度最大为 $n$ 层。双向搜索将状态空间从"单向深度 $n$"压缩为"双向深度 $n/2$"：

• 正向：从 depot 出发，只扩展到 $load\le Q/2$（载重不超过容量一半）
• 反向：从 depot "反向"出发，只扩展到 $load\le Q-Q/2$
• Merge：枚举正向标号 $L_f$（在节点 $v$）和反向标号 $L_b$（在节点 $v$），若 $visite{d}_f\cap visite{d}_b=\varnothing$ 且 $loa{d}_f+loa{d}_b\le Q$，则拼接成完整路线，RC = $r{c}_f+r{c}_b+c_{v_f\to v_b}$

理论状态空间缩减：$2^n\to 2\cdot 2^{n/2}$。$n=18$ 时从 $2^{18}\approx 26$ 万降到 $2\cdot 2^9\approx 1024$，数量级差距。

Merge 的关键实现：正向和反向标号都按 RC 排序，利用双重早停：

// Merge 阶段：正向标号按 RC 升序，反向标号按 RC 升序
// 若 rc_f + min_bwd_rc > best_rc，后续所有 (f, b) 对都不更优，直接跳出
std::sort(fwd_at_v.begin(), fwd_at_v.end(),
          [](const Label* a, const Label* b) { return a->rc < b->rc; });

for (const Label* Lf : fwd_at_v) {
    if (Lf->rc + global_min_bwd_rc + min_return_cost > best_rc) break; // 外层早停

    for (const Label* Lb : bwd_at_v) {
        if (Lf->rc + Lb->rc + arc_cost > best_rc) break;  // 内层早停

        if ((Lf->vmask & Lb->vmask) == 0 &&              // 不重叠
            Lf->load + Lb->load <= capacity) {            // 容量满足
            double rc = Lf->rc + Lb->rc + arc_cost;
            if (rc < best_rc) {
                best_rc = rc;
                // 记录最优合并
            }
        }
    }
}

消融结果：

配置	n=12	n=15	n=18	n=21
C0 裸基线	1ms	0ms	TIMEOUT	SKIP
C3 +Bidir	0ms	0ms	25,574ms	TIMEOUT

指标	n=12	n=18
标号数 C0	2,148	827K*
标号数 C3	1,002	680K
支配检查 C0	343K	176亿*
支配检查 C3	33K	78亿

$n=18$：从不可解变为 25.6s。在我们的消融实验中，没有双向搜索的所有配置（C0、C1、C2、C4）在 $n=18$ 时均超时。理论上，纯 Bucket+RC 排序也能减少支配检查次数，但由于单向搜索的标号总数依然是 $O(2^n)$ 级别，仅靠加速每次支配检查不足以在 60 秒内完成。双向搜索将指数的底数从 $2^n$ 降到 $2^{n/2}$，是突破 $n\ge 18$ 瓶颈的关键。

标号数从 2148 降到 1002（$n=12$），符合 $2^n\to 2\cdot 2^{n/2}$ 的理论预测（比例约 2:1）。

4.2 因子 $T_{per\_label}$：Bucket + RC 排序

朴素支配检查的痛点：标号算法每生成一个新标号，需要与同一节点上所有已有标号做支配检查。朴素实现是 $O(|labels\_at\_node|)$ 全量扫描。当标号数很大时（$n=18$ 时节点上可能累积数千标号），支配检查成为主要时间瓶颈——4.1 节的 Bidir 虽然把标号总数压下来了，但每次支配检查的成本仍是线性的。

Bucket 方案：按 load 分桶（桶数 $\theta =20$），新标号只需检查 load ≤ 自己的桶。桶内按 RC 排序 + min_rc 跳过，进一步减少无效比较。具体实现：

按 load 分桶（桶数 $\theta =20$）：node_buckets[v][b] 存放节点 $v$、载重桶 $b$ 内的所有标号。对载重为 $load$ 的新标号，支配者的载重不超过 $load$，只需检查桶 $0$ 到 $\lfloor load/step\rfloor$：

bool is_dominated(int node, double rc, int load, uint64_t vmask) {
    int b_max = std::min(load / step, n_buckets - 1);
    for (int b = 0; b <= b_max; b++) {
        if (min_rc[node][b] > rc + 1e-10) {
            bucket_skipped++;
            continue;   // 整桶跳过！该桶最小 RC 都大于 rc
        }
        for (int li : node_buckets[node][b]) {
            const Label& el = labels[li];
            if (!el.active) continue;
            if (el.rc <= rc + 1e-10 && el.load <= load
                && (el.vmask & vmask) == el.vmask)
                return true;
        }
    }
    return false;
}

RC 排序：每个桶内的标号按 RC 升序维护。扫描时若遇到 $el.rc>rc$，后续标号的 RC 更大，不可能支配新标号（支配要求 $r{c}_1\le r{c}_2$），立即跳出：

// 桶内 RC 升序排序 + 早停
for (int li : node_buckets[node][b]) {
    const Label& el = labels[li];
    if (el.rc > rc + 1e-10) break;   // RC 排序早停
    if (el.load <= load && (el.vmask & vmask) == el.vmask)
        return true;
}

效果数据：

对比	n=12 支配检查次数	n=15 支配检查次数	n=18 时间
C3 (Bidir only)	33K 次	4.5K 次	25.6s
C6 (Bidir+Bucket)	5.3K 次	536 次	2,378ms
加速比	6x	8x	10.7x

支配检查次数从 78 亿降到 4.9 亿（$n=18$），时间从 25.6s 降到 2.4s，加速 10.7 倍。

Bucket 的盈亏平衡点：桶的维护（插入排序、min_rc 更新）有固定开销。$n=12/15$ 时标号总数仅数百，开销可能抵消收益——实测这两个规模时间都是 0-1ms，差异在测量噪声内。Bucket 的收益只在标号数达到数万以上时才显著。

4.3 因子 $P$：并行 fwd/bwd

原理：双向搜索的正向阶段和反向阶段之间没有共享可变状态——正向生成的标号集合对反向阶段完全透明（反向不依赖正向结果，只依赖图结构）。这是天然适合 2 线程的并行结构：

// 两个线程各自跑一个方向
PhaseResult fwd_result, bwd_result;

std::thread fwd_thread([&]() {
    fwd_result = run_phase(
        /*direction=*/true,   // 正向
        /*load_limit=*/Q / 2, // 正向载重上限
        capacity, n_cust, c_bar, demand, ...);
});

// 主线程跑反向
bwd_result = run_phase(
    /*direction=*/false,  // 反向
    /*load_limit=*/Q - Q/2,
    capacity, n_cust, c_bar, demand, ...);

fwd_thread.join();   // 等待正向完成后再 merge
merge(fwd_result, bwd_result);  // 单线程合并

效果数据：

配置	n=18 时间	加速比
C7 全部无并行	2,377ms	1x
C8 满配	1,415ms	1.68x

接近理论上限 2x（线程调度开销 + Merge 阶段单线程）。$n\le 15$ 时线程创建/同步开销可能大于收益（实测 n=12/15 两者均为 0-1ms）。

为什么只用 2 线程？ 双向搜索天然分为正向和反向两个独立阶段，这是并行度的上限。如果要进一步并行（例如将正向阶段内部的节点扩展并行化），需要处理标号池的并发写入和支配检查的线程安全问题——标号之间的支配关系是全局的，不能简单分区。在我们的测试规模（$n\le 50$）下，线程同步开销会抵消并行收益。RouteOpt 2.0 在处理 $n=1000+$ 的超大实例时确实采用了更复杂的并行策略，但这超出了本文的范围。

累积加速路径（$n=18$）：

C0 裸基线:     TIMEOUT (>60s)
C3 +Bidir:    25.6s   (关键突破)
C6 +Bucket:   2.4s    (10.7x 加速)
C8 +Parallel: 1.4s    (1.68x 加速)

总计：Bidir × Bucket × Parallel ≈ 18x (从 25.6s 到 1.4s)

4.4 无效策略复盘

以下 3 个策略经过严格 A/B 实测，均无效。

4.4.1 Pool 预分配内存（$T_{per\_label}$）

预期效果：std::vector 的动态扩容（每次满载后容量翻倍、重新分配）有额外开销。预先 reserve() 大块内存可以消除 reallocation。

实测结果：

对比	n=12	n=15	n=18
C0 vs C1	1ms vs 1ms	0ms vs 0ms	均超时
C3 vs C5	0ms vs 0ms	0ms vs 0ms	25.6s vs 25.5s
C6 vs C7	0ms vs 0ms	0ms vs 0ms	2,378ms vs 2,377ms

所有对照组完全一致。原因：现代 std::vector 的分摊 O(1) 扩容成本极低；标号总数不超过 200 万时，内存分配不是瓶颈。Pool 是"正确的工程习惯"（零代价保留），但不计入有效加速策略。

需要指出的是，Pool 预分配的收益取决于标号总数——当标号数不足以触发 vector 多次 reallocation 时（我们的测试中 $n\le 18$ 最大约 40 万标号），其效果自然不可观测。在更大规模（百万级标号）的场景下，Pool 可能会有可观测的收益，但我们的测试集无法验证这一点。

4.4.2 Completion Bounds 剪枝（$|Labels|$）

原理：预计算 0-1 背包：max_reward[q] = 在剩余容量 $q$ 内能获得的最大对偶收益。若 new_rc ≥ max_reward[Q - new_load]，即使后续完美利用剩余容量也无法获得负 RC，直接剪枝。

// Completion Bounds 预计算（0-1 背包）
std::vector<double> max_reward(capacity + 1, 0.0);
for (int i = 0; i < n_cust; i++) {
    if (cover_duals[i] <= 0) continue;
    for (int q = capacity; q >= demand[i+1]; q--)
        max_reward[q] = std::max(max_reward[q], max_reward[q - demand[i+1]] + cover_duals[i]);
}

// 扩展时剪枝
if (new_rc - max_reward[capacity - new_load] >= -1e-6) {
    labels_pruned++;
    continue;  // 剪枝：不可能产生负 RC 路径
}

pruning 统计（看似很好）：

实例	CG 阶段	pruned 比例
en13k4_iter1	早期	0%
en13k4_iter2	中期	~10%
pn16k8_converged	收敛	~48%

A/B 实测结果（精确满配 + CB vs 精确满配）：

实例	n	+CB	无CB	加速比
en13k4_iter1	12	0.70ms	0.62ms	0.87x（慢）
en13k4_iter2	12	0.52ms	0.31ms	0.59x（慢）
en13k4_converged	12	0.75ms	0.46ms	0.61x（慢）
pn19k2_initial	18	1541ms	1522ms	0.98x（无效）

核心教训：pruning 统计 ≠ 实际加速。Bucket + RC 排序体系下，被 CB 剪掉的标号大概率已经被支配消除——两种机制剪的是同一批"弱标号"，CB 只是多做了一次检查（浮点比较 + 查表），反而多了开销。大实例（$n\ge 18$）CB 完全无效（pruned=0）。

4.4.3 弧消除（$T_{per\_label}$）

原理：对每条弧 $(i,j)$，用 Bellman-Ford 算法计算对偶图中的最短路下界 $h_i,h_j$。若 ${\bar{c}}_{ij}+h_j-h_i\ge 0$，则弧 $(i,j)$ 不可能参与负 RC 路径，删除之以缩小图。

实测结果：所有实例 arcs_eliminated = 0。

根本原因：CVRP 覆盖约束对偶使修正边权图普遍存在负环。例如：
$${\bar{c}}_{12}+{\bar{c}}_{21}=(c_{12}-{\pi }_2)+(c_{21}-{\pi }_1)=2c_{12}-{\pi }_1-{\pi }_2$$
当 ${\pi }_1+{\pi }_2>2c_{12}$（对偶值较大时），这就是负数。Bellman-Ford 遇到负环时将相关节点的 bound 设为 $-\infty$，保守地保留所有弧。即使 CG 收敛时的对偶也有负环——这是 CVRP 修正边权图的结构性特征。

结论：ESPPRC 的障碍不在于"有负权弧"，而在于状态空间指数级。已有策略（ng-route 压缩 mask、双向减半深度、Bucket 加速支配）是正确方向。

4.5 精确型小结

完整消融表

配置	n=12 时间	n=15 时间	n=18 时间	n=21 时间
C0 裸基线	1ms	0ms	TIMEOUT	SKIP
C1 +Pool	1ms	0ms	TIMEOUT	SKIP
C2 +Bucket	1ms	0ms	TIMEOUT	SKIP
C3 +Bidir	0ms	0ms	25,574ms	TIMEOUT
C4 Pool+Bucket	1ms	1ms	TIMEOUT	SKIP
C5 Pool+Bidir	0ms	0ms	25,481ms	TIMEOUT
C6 Bucket+Bidir	0ms	0ms	2,378ms	TIMEOUT
C7 全部无并行	0ms	1ms	2,377ms	TIMEOUT
C8 满配	1ms	1ms	1,415ms	TIMEOUT

策略独立性验证

各策略加速比（$n=18$）：

• Bidir：TIMEOUT → 25.6s（必要，不量化）
• Bucket（Bidir 开启）：25.6s → 2.4s = 10.7x
• Parallel：2.4s → 1.4s = 1.68x
• Pool：无可观测差异 = 1.0x

乘性组合：$10.7\times 1.68\approx 18x$（从 25.6s 到 1.4s），与实测完全吻合。策略间无交互效应，符合"打在独立因子"的预期。

精确型极限

$n=21$，满配 C8 仍然 TIMEOUT。从 $n=18$ 到 $n=21$ 只增加 3 个客户，但理论状态空间增加 $2^{(21-18)/2}\approx 2.8$ 倍，加上图密度增加，实际计算量可能 5-10 倍。精确 bitmask 算法的 $O(n\cdot 2^n)$ 天花板无法通过工程优化突破——需要从根本上缩减状态空间。

精确型三因子极限分析

回顾复杂度分解 $T_{pricing}=|Labels|\times T_{per\_label}/P$，精确型三因子的优化状态：

因子	策略	效果	到极限了吗？
$|Labels|$	Bidir（深度 $n\to n/2$）	TIMEOUT → 25.6s	✅ 理论极限，深度无法更短
$|Labels|$	Completion Bounds	无效	❌ 剪枝目标与支配重叠
$T_{per\_label}$	Bucket + RC 排序	10.7x	✅ min_rc 已实现 O(1) 跳过
$T_{per\_label}$	Pool 预分配	无效	— 当前规模不触发
$T_{per\_label}$	弧消除 (Bellman-Ford)	无效	❌ CVRP 结构性负环
$P$	并行 fwd/bwd	1.68x	✅ 2 线程 = 实际上限

结论：精确型的三个因子均已压至工程极限。瓶颈不在任何单一因子的优化不足，而在状态空间本身——$O(n\cdot 2^n)$ 的指数底数无法通过工程手段消除。必须从根本上缩减状态空间，这正是下一节松弛型的出发点。

---

衔接：从精确到松弛的设计空间

精确型消融已经做到工程极限：双向 + Bucket + 并行，$n=18$ 需要 1.4 秒，$n=21$ 仍然超时。问题的根源是状态空间 $O(n\cdot 2^n)$，工程优化无法改变指数的底数。

如果不能再压缩单标号处理代价，就必须减少标号的数量——也就是缩小状态空间。从完全精确（记住所有 $2^n$ 种访问组合）到完全松弛（SPPRC，不记访问状态），中间存在一个连续的设计空间：

松弛程度	方法	状态空间	路径性质
完全精确	Full bitmask	$O(n\cdot 2^n)$	严格 elementary
k-cycle 消除	禁止长度 $\le k$ 的环	$O(n^k)$（多项式）	禁止短环
ng-route 松弛	只在邻域内保持 elementary	$O(n\cdot 2^{\Delta })$	近似 elementary
完全松弛	SPPRC	$O(n^{2}Q)$（伪多项式）	允许重访

k-cycle 消除（禁止长度 $\le k$ 的环）是一种中间策略：$k=2$ 时禁止 2-cycle（消除振荡），$k=n$ 时退化为完全 elementary。理论上 $k$ 越大约束越紧，但状态空间随 $k$ 多项式增长，计算代价也相应增加。实践中 $k$ 很少超过 5，而且单纯 k-cycle 消除并不足以充分利用空间结构信息。

我们实现了 k-cycle 消除（模板参数 $K$，BFS 标号扩展，hash map 按 $(node,recent[])$ 分组做支配检查）并与 ng-route 实测对比：

配置	n=12 时间	n=12 RC	n=15 时间	n=15 RC	n=31 时间	n=60 时间
k-cycle $k=2$	0.14ms	-296	0.08ms	-230	545ms	1,202ms
k-cycle $k=3$	0.59ms	-296	0.15ms	-196	33s	超时
k-cycle $k=5$	1.38ms	-296	0.20ms	-196	SKIP	SKIP
ng-route $\Delta =5$	0.14ms	-240	0.05ms	-148	52ms	99ms
ng-route $\Delta =8$	0.59ms	-240	0.09ms	-131 ✓	86ms	366ms
精确	1.02ms	-240	0.11ms	-131	SKIP	SKIP

（n=12 精确 RC = -240；n=15 精确 RC = -131；✓ 表示与精确值吻合）

实测揭示了两个反直觉的事实：

k-cycle 既慢又松。对比 $k=3$ 与 ng-route $\Delta =8$：在 $n=15$ 实例上，k-cycle $k=3$ 耗时 0.15ms 但 RC = -196，而 ng-route $\Delta =8$ 耗时 0.09ms 且 RC = -131（精确值）。k-cycle 用了更多时间，却得到了更差的 LP bound（更负的 RC = 更松的松弛）。在 $n=31$ 时这一差距暴增：k-cycle $k=3$ 需要 33 秒，ng-route $\Delta =8$ 仅需 86 毫秒，而且 ng-route 的 bound 质量更高。

$k$ 增大不一定收紧 bound。k=3 与 k=5 在 $n=12$ 和 $n=15$ 上给出完全相同的 RC，说明这些实例中最优路径的 cycle 长度超过 5——禁止更短的 cycle 对 bound 无影响，但计算量仍在上升。

根本原因在于：k-cycle 按"时序邻近度"（最近 $k-1$ 步）约束，而 CVRP 的危险重访是"空间邻近度"（最近的节点之间）。ng-route 直接针对后者，用更少的状态获得更紧的约束。

ng-route（Baldacci et al., 2011） 是这个设计空间中最受工程欢迎的选择：它不按环长度限制，而是按空间邻近度维护 elementary 约束——只在每个节点的 $\Delta$ 个最近邻范围内禁止重访。状态空间从 $2^n$ 压缩到 $2^{\Delta }$，$\Delta =8$ 时缩减因子超过 $2^{n-8}$。更关键的是，ng-route 松弛的质量极好：在实践中，最近邻之间的负 2-cycle 才是真正的危险，远处节点的重访几乎不出现在最优路线中。第 5 章将详细说明这一点。

---

5. ng-route 松弛：突破指数壁垒

5.1 因子 $|Labels|$：ng-route 邻域遗忘

核心思想：精确算法必须记住"曾经访问过哪些客户"（完整 bitmask），以保证不重复。ng-route 的创新在于：只记忆每个节点的 $\Delta$ 个最近邻的访问状态。远处的节点可以被"遗忘"，从而允许再次访问。

状态转移（对比精确 vs ng-route）：

// 精确 ESPPRC：全集记忆
new_vmask = old_vmask | (1ULL << (j - 1));   // 永久记住访问了 j

// ng-route 松弛：邻域遗忘
new_vmask = (old_vmask & ng_mask[j]) | (1ULL << (j - 1));
//              ↑ 只保留 j 邻域内的访问记录        ↑ 记住 j 自身

ng_mask[j] 是节点 $j$ 的 $\Delta$ 个最近邻的 bitmask。访问 $j$ 时，不在 $j$ 邻域内的客户访问记录被清零——这些客户可以再次访问，但实际上它们离 $j$ 很远，在优质路线中一般不会出现绕路。

ng_mask 的构建：

// build_ng_masks: 按距离选最近的 Δ 个邻居
std::vector<uint64_t> build_ng_masks(
    const std::vector<std::vector<int>>& dist,
    int n_cust, int delta)
{
    std::vector<uint64_t> ng_mask(N, 0u);
    for (int j = 1; j <= n_cust; j++) {
        // 按距离排序，取最近的 delta 个
        std::vector<int> others(n_cust);
        std::iota(others.begin(), others.end(), 1);
        std::sort(others.begin(), others.end(),
                  [&](int a, int b) { return dist[j][a] < dist[j][b]; });

        uint64_t mask = (1ULL << (j - 1));  // 始终包含 j 自身
        int cnt = 0;
        for (int k : others) {
            if (k == j) continue;
            mask |= (1ULL << (k - 1));
            if (++cnt >= delta) break;
        }
        ng_mask[j] = mask;
    }
    return ng_mask;
}

复杂度：$O(n\cdot 2^{\Delta })$。$\Delta =8$ 时状态空间每个节点最多 $2^8=256$ 种 vmask，总状态 $\le n\cdot 256$。对比精确的 $n\cdot 2^n$，$n=25$ 时缩减因子超过 $2^{17}$。

理论代价：ng-route 得到的是 ESPPRC 的下界（RC 可能比精确更小），路径可能不满足 elementary 约束（有重复节点）。但实践中 $\Delta =8$ 的精度已非常接近精确（见 5.3 节）。

为什么 ng-route 有效？直觉解释

ng-route 有效的关键在于：负 2-cycle 的空间局部性。

回顾 2.5.2 节，负 2-cycle 的条件是：
$${\bar{c}}_{jk}+{\bar{c}}_{kj}=c_{jk}+c_{kj}-{\pi }_j-{\pi }_k<0$$

由于距离满足三角不等式，$c_{jk}+c_{kj}\ge 2\cdot {\min }_{路径}$，空间越近的节点对 $c_{jk}+c_{kj}$ 越小，越容易满足负 2-cycle 条件。因此，最危险的重复访问模式集中在空间相邻的节点对。

ng-route 恰好在每个节点的 $\Delta$ 个最近邻（空间最近的邻居）之间维护 elementary 约束。这意味着：

• 近邻节点对（最容易形成负环）：被 ng_mask 约束禁止重访，消除了最有害的振荡
• 远距离节点对（很难形成负环）：允许重访，但 $c_{jk}+c_{kj}$ 很大，算法自然不会选择这些绕路

这就是为什么 $\Delta =8$——仅 8 个最近邻——就能让 LP bound 与完全精确达到相同质量：精确定价所做的额外工作（记住远处节点的访问状态），在优质路线中几乎从不被用到。

5.1.1 ng-route 为何无法用 MIP 建模

既然 ng-route 松弛的效果如此接近精确，一个自然的问题：能否将 ng-route 的约束也写成 MIP，交给 Gurobi 求解？

答案是否定的。ng-route 的约束本质上是路径依赖的：

new_vmask = (old_vmask & ng_mask[j]) | {j}   // 路径依赖

new_vmask 的值不仅取决于"当前访问了 $j$"，还取决于"沿这条路径到达 $j$ 之前积累的访问历史"。MIP 约束是变量之间的全局线性关系，无法直接表达"沿某条路径的访问历史"。

要在 MIP 中捕获 ng-route 的语义，需要为每个（节点, vmask）状态对引入变量——变量数达到 $O(n\cdot 2^{\Delta })$，并且约束仍然是指数级的。这与直接跑标号算法相比没有任何优势，反而丢失了标号算法中的支配消除机制。

结论：ng-route 只能用标号算法实现。这也是为什么高性能 ESPPRC 库都以 C++ 标号算法为核心，而不是调用通用 MIP 求解器。

5.2 全栈组合：BIDIR_POOL

将 ng-route 与前述所有精确型优化组合，形成 BIDIR_POOL 求解器：

BIDIR_POOL = ng-route + 双向标号 + PoolEngine(RC排序桶 + min-RC跳过 + 预分配) + 并行

PoolEngine 的 is_dominated（关键数据结构）：

struct PoolEngine {
    std::vector<Label> labels;                          // 预分配标号池
    std::vector<std::vector<std::vector<int>>> node_buckets;  // [node][bucket][label_idx]
    std::vector<std::vector<double>> min_rc;             // [node][bucket] 最小 RC

    bool is_dominated(int node, double rc, int load, uint64_t vmask) {
        const Label* ldata = labels.data();
        int b_max = std::min(load / step, n_buckets - 1);
        for (int b = 0; b <= b_max; b++) {
            if (min_rc[node][b] > rc + 1e-10) {
                bucket_skipped++;
                continue;   // 整桶跳过
            }
            for (int li : node_buckets[node][b]) {
                const Label& el = ldata[li];
                if (!el.active) continue;
                if (el.rc > rc + 1e-10) break;    // RC 排序早停
                if (el.load <= load && (el.vmask & vmask) == el.vmask)
                    return true;
            }
        }
        return false;
    }
};

与精确满配（C8）相比，ng-route 在不降低支配规则的情况下将 vmask 的位数从 $n$ 压缩到 $\Delta$，使得更多标号相互支配、更快被消除。

5.3 性能爆发

$n$	精确满配 C8	松弛 BIDIR_POOL
12	0.5ms	0.35ms
15	0.8ms	0.5ms
18	1,415ms	10ms
21	TIMEOUT	9ms
31	TIMEOUT	8ms

$n=18$：从 1.4s 到 10ms，140 倍加速。$n=21,31$：精确不可行，松弛 10ms 以内。

这不是渐进改进，而是数量级跳跃。

5.4 LP 质量评估：接近免费的午餐

松弛定价得到的是 LP 下界——但下界有多宽松？用 4 个 BPC 实例对比松弛 LP 和精确 LP 的差距：

实例	$n$	松弛 LP（$\Delta =8$）	精确 LP	差值	精确代价
E-n13-k4	13	264.00	264.00	0	1.0x
E-n22-k4	22	373.50	373.71	0.21	1.9x
A-n32-k5	32	758.43	758.43	0	4.4x
E-n51-k5	51	517.06	517.06	0	12.5x

4 个实例中 3 个 LP 值完全相同，1 个仅差 0.21（<0.06%）。而精确定价的代价是 1.9-12.5 倍以上。

结论：$\Delta =8$ 的 ng-route 是"几乎免费的午餐"——以极低的 LP 精度损失换取数量级的速度提升。BPC 列生成的默认配置：全程使用 BIDIR_POOL（$\Delta =8$）。

5.5 因子 $|Labels|$：DSSR 从松弛渐进精确

当需要精确的 elementary 解时（如 CG 收敛验证或 L3 定价），使用 DSSR（Decremental State Space Relaxation）。

算法逻辑：

初始化：critical_mask = 0（所有节点均使用 ng-route 松弛）

循环：
    1. 运行 ng-route（携带 critical_mask）：
       new_vmask = (old_vmask & (ng_mask[j] | critical_mask)) | {j}
    2. 检查最优路径是否有重复节点
    3. 若有重复节点 v：将 v 加入 critical_mask
    4. 重复，直到无重复节点（elementary）或无新增 critical 节点

DSSR 状态转移的含义：critical 节点被所有路线强制"记忆"（永不遗忘），相当于逐步收紧到精确。通常 1-3 轮迭代后收敛。

DSSR vs Full Bitmask 精确对比：

实例	$n$	DSSR	轮数	critical节点	精确 Full Bitmask
en13k4_iter1	12	1.1ms	1	0	2.8ms
pn16k8_converged	15	0.2ms	1	0	1.3ms
pn16k8_iter1	15	1.0ms	2	2	1.6ms
pn19k2_perturbed	18	12,574ms	3	8	TIMEOUT
pn22k2_initial	21	TIMEOUT	4	12	TIMEOUT

DSSR 在纯 ESPPRC 场景全面优于 Full Bitmask（1.5-7x 加速）。$n=18$ 时 DSSR 可解（虽然慢）而精确不行。

使用场景：CG 末期验证收敛（L3 定价）。RouteOpt 2.0 在含 R1C cuts 时选 Full Bitmask，因为 cuts 与 bitmask 更容易集成。

---

5.6 松弛型微优化 A/B 测试

BIDIR_POOL 已是主力松弛求解器，我们对 3 个"常识性"优化做了编译时开关 A/B 实测，结论在意料之外。

实验设置：4 个版本编译自同一 bench_optimization.cpp，只通过 -D 标志切换：baseline（全关）/ BIN_LIMIT=128 / ADAPTIVE_MEET=1 / DIST_SORT=1。各版本运行全部测试集实例，3 遍取中位数。

5.6.1 DistSort：邻居按距离排序（$T_{per\_label}$，有效）

原理：构建前向/后向邻居表时，按距离升序排列。近邻先被扩展，更快产生低 RC 标号，有望提早触发支配剪枝。

A/B 结果：

实例	n	Baseline	DistSort	加速比	Baseline 标号数	DistSort 标号数
pn19k2_initial	18	12.7ms	8.7ms	1.45x	19,086	12,455 (-35%)
an32k5_initial	31	11.9ms	9.9ms	1.20x	34,710	30,622 (-12%)
an32k5_perturbed	31	31.0ms	23.8ms	1.30x	62,648	54,722 (-13%)
large50_random	50	7.3ms	6.8ms	1.08x	5,668	4,888 (-14%)
large60_random	60	23.6ms	21.4ms	1.10x	15,024	12,102 (-19%)

所有实例 RC 值完全一致——算法正确性不受影响。中等规模实例（n=18-31）标号数减少 12-35%，加速 1.2-1.5x；大实例（n=50-60）效果收窄至 1.1x。机制符合预期：近邻先扩展 → 低 RC 标号更早进入桶 → 远邻标号被更早支配消除。

工程建议：默认开启。一行排序代码换来约 1.2x 中位加速，零风险。

5.6.2 BIN_LIMIT：每桶标号数上限（$|Labels|$，无效）

原理：仿照 RouteOpt 实现，每个 (node, load_bin) 桶最多保留 128 个 RC 最小的标号，超额丢弃 RC 最大的。减少标号数同时可能损失最优性（以启发式换速度）。

A/B 结果：label 数与 baseline 完全一致，时间也无可观测差异。

原因：我们测试集中所有实例的 (node, load_bin) 桶最大规模均未超过 128——上限从未触发。BIN_LIMIT 是一个"大实例才生效"的策略，在 $n\le 60$ 的测试集上无法评估。RouteOpt 在有 R1C cuts 后标号数暴增才需要此机制。

5.6.3 AdaptiveMeet：自适应双向汇合点（$|Labels|/P$，无效）

原理：固定 half = Q/2 假设前后向容量消耗对称。实际 demand 分布不对称时，调整汇合容量可平衡标号数。

实现：基于均值 demand 预估前后向深度，若比值超过 1.5 则微调 half。

A/B 结果：与 baseline 完全一致（标号数、RC、时间均无差异）。

原因：简化实现的启发式调整幅度（$\pm Q/10$）过于保守，对典型 CVRP 实例（demand 集中在 Q/5-Q/3 区间）几乎不触发调整。真正有效的自适应汇合点需要在线监控首次扩展的标号深度分布，实现复杂度显著提升。

5.7 松弛型三因子极限分析

因子	策略	效果	到极限了吗？
$|Labels|$	ng-route $\Delta =8$	$2^n\to 2^8$	✅ 实用甜点；$\Delta$ 增大边际递减
$|Labels|$	DSSR 迭代	1-3 轮收敛	✅ 已是该策略最优表现
$T_{per\_label}$	Bucket + RC 排序（同精确型）	10.7x	✅ 直接复用
$T_{per\_label}$	DistSort	1.2-1.5x	✅ 已默认启用
$T_{per\_label}$	BIN_LIMIT=128	未触发	— 大实例 + R1C 后有效
$P$	并行 fwd/bwd（同精确型）	1.68x	✅ 直接复用

结论：松弛型的核心突破在 $|Labels|$ 因子——ng-route 将状态空间指数从 $2^n$ 压缩到 $2^{\Delta }$，这是从"不可行"到"毫秒级"的关键。$T_{per\_label}$ 和 $P$ 的策略与精确型完全复用，无需额外优化。三因子均已到达实用极限，进一步提升需要 CG 级策略（动态 ng 邻域、对偶稳定化）。

---

6. 对标开源 SOTA

6.1 vs cspy（学术界最常用开源库）

cspy 是目前学术界最广泛使用的开源 ESPPRC/SPPRC 求解库，被众多 BPC 论文用作定价基准。cspy 虽以 Python 接口为主，但核心算法以 C++ 实现并通过绑定层暴露给 Python/C#，因此性能差距并非语言层面的。

对标结果：

对比场景	我们的求解器	cspy	加速倍数
松弛定价（BIDIR_POOL vs cspy SPPRC）	8-176ms	800ms-176s	100-1,000x
精确定价（EXACT_POOL vs cspy elementary）	1-28ms	0.4-31s	100-31,000x

加速来源分析（按复杂度因子）：

因子	我们的策略	cspy	估计贡献
$|Labels|$	Bidir + ng-route ($\Delta =8$)	单向 + 朴素 elementary	主要，标号数减少 1-2 个数量级
$T_{per\_label}$	Bucket + RC 排序 + DistSort	链表遍历式支配检查	显著，单标号代价降低 5-20x
$P$	双线程 fwd/bwd	单线程	1.5-2x
接口开销	纯 C++ 直接调用	Python→C++ 绑定层	较小，~2-5x

差距的主因不是语言，而是算法和数据结构：cspy 提供的是通用框架（支持自定义资源和回调），而我们针对 CVRP 定价做了专用优化。

6.2 vs RouteOpt 2.0（INFORMS JoC 2025）

RouteOpt（You & Yang, 2025）是目前 CVRP/VRPTW 精确求解的学术 SOTA，发表于 INFORMS Journal on Computing（运筹学 Tier 1 期刊，仅次于 Management Science 和 Operations Research；IJoC 设有"Software Tools"专栏，是计算型求解器论文的最佳归宿），在标准测试集上全面超越 VRPSolver，并首次证明了三个 CVRP 开放实例的最优性。

技术对标：

ESPPRC 单次定价技术（本文研究范围）：

技术	RouteOpt 2.0	我们的库	备注
ng-route 松弛（$\Delta =8$）	✓	✓	状态空间压缩
双向标号	✓	✓	深度减半
Pool / RC 排序	✓（Bucket Graph）	✓（PoolEngine）	不同架构，效果相近
并行 fwd/bwd	✓	✓	双线程
DSSR 精确化	✓	✓	迭代收紧
DistSort（距离排序扩展）	未提及	✓	1.2-1.5x 加速
Jump Arcs	✓（VRPTW）	❌ 不适用	需 2+ 资源维度

结论：在纯 ESPPRC 单次定价层面，我们的库已覆盖 RouteOpt 的全部相关技术。Jump Arcs 仅适用于多资源问题（VRPTW），不适用于 CVRP。

架构差异：RouteOpt 使用 Bucket Graph（按资源值分桶，桶间有转移弧），我们使用 PoolEngine（按载重分桶 + RC 排序）。两者在工程效果上接近，但 Bucket Graph 与 Jump Arcs 的集成更自然，且更容易扩展到多资源维度（VRPTW）。

超出 ESPPRC 范围的 BPC 框架级技术（不在本文研究范围内）：

技术	层级	说明
三阶段定价调度（L1/L2/L3）	CG 级	自动升降级，减少总定价次数
对偶稳定化	CG 级	减少对偶震荡，加速收敛
大规模弧消除（RC fixing）	CG 级	利用 LP bound 预删弧
R1C cuts（Rank-1 Cuts）	B&B 级	核心 LP 紧化，节点数减少 80-95%

这些技术作用于列生成循环或 B&B 树层面，而非单次定价调用内部。它们是 BPC 求解器的整体竞争力来源，但不改变 ESPPRC 定价器本身的实现质量。

---

7. 总结与工程推荐

7.1 技术叠加总表

Level	技术	打击因子	$n=18$ 效果	适用范围
0	基础标号（bitmask）	基线	TIMEOUT	$n\le 15$
1	+ 双向搜索	$|L|$ 深度	25.6s	$n\le 18$
2	+ Bucket + RC 排序	$T_{per\_label}$	2.4s	$n\le 18$
3	+ 并行	$P$	1.4s	$n\le 18$
4	ng-route（$\Delta =8$）	$|L|$ 宽度	10ms	$n\le 100+$
4+	BIDIR_POOL（全栈）	全部	10ms	实用推荐
5	+ DSSR（L3精确）	正确性	秒级（偶尔）	CG 末期验证

7.2 三因子极限总览

贯穿全文的复杂度分解 $T_{pricing}=|Labels|\times T_{per\_label}/P$，三个因子的最终状态：

因子	精确型最佳策略	松弛型最佳策略	极限判定
$|Labels|$	Bidir（深度 $n\to n/2$）	ng-route（宽度 $2^n\to 2^{\Delta }$）	✅ 均已到极限
$T_{per\_label}$	Bucket + RC 排序（10.7x）	+ DistSort（额外 1.2x）	✅ 接近理论下界
$P$	2-thread fwd/bwd（1.68x）	同左	✅ 实际上限

三个因子均已压至工程极限。ESPPRC 单次定价的优化空间已基本穷尽——下一步提升不在定价器内部，而在 CG 循环层面：动态 ng 邻域（减少定价调用次数）、对偶稳定化（加速收敛）、三阶段调度（自动选择松弛级别）。

7.3 工程推荐配置

BPC 列生成的推荐配置如下：

三层定价架构：

┌─────────────────────────────────────────────────────┐
│  L1 Light   BIDIR_POOL(Δ=4~6)                       │
│             CG 早期，负 RC 列多时使用，极快出列        │
├─────────────────────────────────────────────────────┤
│  L2 Heavy   BIDIR_POOL(Δ=8~12)                      │
│             L1 找不到负 RC 列时升级，高质量列          │
├─────────────────────────────────────────────────────┤
│  L3 Exact   DSSR_POOL (纯 ESPPRC)                   │
│             espprc_exact_pool (含 R1C cuts)          │
│             L2 也找不到时，精确确认 CG 收敛            │
└─────────────────────────────────────────────────────┘

特定场景推荐：

场景	推荐配置	理由
$n\le 15$，验证正确性	精确 bitmask（Level 0）	最简单，作为 oracle
$n\le 18$，需要精确解	EXACT_POOL（C8 满配）	工程上限
$n\le 31$，精确 LP	DSSR_POOL	1-3 轮迭代，通常秒级
BPC 列生成主力定价	BIDIR_POOL（$\Delta =8$）	速度与质量最优均衡
CG 收敛验证	DSSR_POOL 或 EXACT_POOL	必须精确 elementary

不推荐的配置：

• Completion Bounds（pruning 有但不加速）
• Pool 预分配（无可观测效果）
• Bellman-Ford 弧消除（CVRP 负环导致无效）

7.4 进一步研究方向

在 ESPPRC 定价器层面，仍有若干值得探索的方向：

• 多资源扩展：当前仅支持容量约束（单资源）。扩展到时间窗（VRPTW）后，双向搜索的 meet-point 选择、Jump Arcs 等技术将变得有价值
• 动态 ng 邻域：在 CG 迭代中根据对偶值变化自适应调整 $\Delta$ 和邻域构成，平衡松弛质量与求解速度
• 标号池压缩：大规模实例（$n>100$）下标号池内存可能成为瓶颈，label-per-bin limit 等策略有待在更大实例上验证

此外，ESPPRC 定价器作为 BPC 框架的核心组件，其性能最终需要在完整的列生成循环中评估。三阶段定价调度（L1→L2→L3 自动升降级）、对偶稳定化等 CG 级策略将进一步释放定价器的潜力，这些属于 BPC 框架研究的范畴。

---

参考文献

1. Righini, G., & Salani, M. (2006). Symmetry helps: Bounded bi-directional dynamic programming for the elementary shortest path problem with resource constraints. Discrete Optimization, 3(3), 255-273.

2. Righini, G., & Salani, M. (2008). New dynamic programming algorithms for the resource constrained elementary shortest path problem. Networks, 51(3), 155-170.

3. Sadykov, R., Uchoa, E., & Pessoa, A. (2021). A bucket graph-based labeling algorithm with application to vehicle routing. Transportation Science, 55(1), 4-28.

4. Irnich, S., & Desaulniers, G. (2005). Shortest path problems with resource constraints. In Column Generation (pp. 33-65). Springer.

5. Baldacci, R., Mingozzi, A., & Roberti, R. (2011). New route relaxation and pricing strategies for the vehicle routing problem. Operations Research, 59(5), 1269-1283.

6. Bulhões, T., Hà, M. H., Martinelli, R., & Vidal, T. (2018). The vehicle routing problem with service level constraints. European Journal of Operational Research, 265(2), 544-558.

7. You, Z., & Yang, Y. (2025). RouteOpt: An open-source modular exact solver for vehicle routing problems. INFORMS Journal on Computing.

8. Torretta, C. J. (2022). cspy: A Python package with a collection of algorithms for the (resource constrained) shortest path problem. Journal of Open Source Software, 7(76), 4014.

9. Desaulniers, G., Desrosiers, J., & Solomon, M. M. (2005). Column Generation. Springer.

---

完整代码与使用指南

本文实验的 ESPPRC 求解器库已开源，完整代码位于 GitHub 仓库 的 espprc/ 目录。

快速使用

#include "espprc_solver.h"

// 输入
std::vector<std::vector<int>> dist = {...};  // N×N 距离矩阵, N = n_customers + 1
std::vector<int> demand = {...};              // 需求量, demand[0] = 0 (depot)
int capacity = 100;
std::vector<double> duals = {...};            // 覆盖约束对偶值 π[1..n_cust]

// 求解（默认 BIDIR_POOL，推荐配置）
auto r = espprc::solve(dist, demand, capacity, duals);

if (r.has_negative_rc()) {
    // r.optimal_path = [0, v1, v2, ..., vk, 0]
    // r.optimal_rc   = 最小 reduced cost
    // 将 r.optimal_path 作为新列加入 RMP
}

编译

依赖：C++17 编译器（GCC/Clang/MSVC），无需外部库。

cd espprc && mkdir build && cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
cmake --build . --config Release

文件索引

文件	说明
espprc_solver.h	统一接口（推荐入口）
relaxed/espprc_ng_bidir_pool.h	BIDIR_POOL 全栈（推荐生产配置）
relaxed/espprc_ng.h	ng-mask 构建工具
exact/espprc.h	Level 0 基础标号（精确 bitmask oracle）
exact/espprc_exact_modular.h	模块化精确求解器（4 个 runtime 开关，消融实验用）
exact/espprc_dssr_pool.h	DSSR + Pool（L3 精确定价首选）
exact/espprc_exact_pool.h	Full Bitmask + Pool（含 cuts 时 L3 备选）
bench_exact_ablation.cpp	9 配置 × 4 实例消融实验
bench_espprc.cpp	综合基准测试（11 个 Part + 700 实例压力测试）

---

一起交流

如果你觉得有帮助，欢迎：

• 📌 订阅本专栏，跟进后续更新
• 💬 评论区留言，交流你的想法
• ⭐ 点赞收藏，让更多朋友看到

---

— Sebastilan & Claude