用Python手把手实现SPPRC标签算法：从理论到代码实战（附Solomon数据集）
·
用Python手把手实现SPPRC标签算法：从理论到代码实战（附Solomon数据集）
如果你正在学习运筹优化，或者是一名Python开发者，想要深入理解并实现一个能解决实际工业问题的算法，那么这篇文章就是为你准备的。我们不会停留在枯燥的理论推导上，而是直接动手，用Python从零开始构建一个解决带资源约束的最短路径问题的标签算法。我会带你走过从问题定义、算法核心思想，到代码实现的每一个关键步骤，最后用经典的Solomon数据集进行实战测试，并可视化我们的求解结果。整个过程，你会看到如何将复杂的运筹学模型，转化为清晰、高效、可复用的Python代码。
1. 问题拆解：什么是SPPRC与ESPPRC？
在开始写代码之前，我们必须先弄清楚要解决什么问题。你可能会看到两个相似的名词：SPPRC和ESPPRC。它们都是“最短路径问题”这个大家族中，加了“紧箍咒”的成员。
- SPPRC：全称是 Shortest Path Problem with Resource Constraints，即带资源约束的最短路径问题。这里的“资源”可以是时间、车辆载重、电量等等。路径必须满足这些资源的总消耗不超过上限。关键点在于，它允许路径重复访问同一个节点。这听起来有点反直觉，但在数学上，这种松弛让问题变得“友好”一些——存在伪多项式时间的精确算法。
- ESPPRC：全称是 Elementary Shortest Path Problem with Resource Constraints，即带资源约束的基本最短路径问题。它比SPPRC多了一个“Elementary”约束，意思是路径中每个节点最多只能访问一次。别小看这个限制，它瞬间将问题的难度提升到了NP-hard级别。
提示：在实际应用中，比如用列生成法求解车辆路径问题（VRP）时，子问题往往就是ESPPRC。但因为ESPPRC太难解，我们通常会先把它松弛成SPPRC来求一个下界。所以，掌握SPPRC的解法是攻克更复杂问题的重要基石。
为了更直观地理解它们的区别和联系，我们用一个简单的表格来对比：
| 特性 | SPPRC | ESPPRC |
|---|---|---|
| 全称 | Shortest Path Problem with Resource Constraints | Elementary Shortest Path Problem with Resource Constraints |
| 节点访问 | 允许重复访问 | 禁止重复访问（基本路径） |
| 计算复杂度 | 伪多项式时间可解 | NP-hard问题 |
| 典型应用 | 列生成中的松弛子问题 | 列生成中的精确子问题（理论上） |
| 关系 | ESPPRC松弛掉“基本路径”约束后得到 | SPPRC加上“基本路径”约束后得到 |
我们的实战将从SPPRC入手，因为它的算法框架更清晰，同时也是理解ESPPRC算法的基础。理解了如何用标签法处理资源约束，后续再加入“禁止重复访问”的逻辑就会容易得多。
2. 算法核心：标签法是如何工作的？
标签法，特别是标签校正算法，是解决SPPRC的一把利器。你可以把它想象成在一个迷宫中，派出了许多支探险队，每支队伍都记录着自己的“履历”（标签），然后根据一套规则决定哪些队伍值得继续探索，哪些队伍可以提前收队。
一个标签本质上是一个状态记录，它描述了某条部分路径抵达某个节点时的全部重要信息。通常，一个标签会包含以下核心字段：
# 一个标签类的简化结构示意
class Label:
def __init__(self):
self.current_node = None # 当前所在的节点
self.previous_label_index = None # 上一个标签的索引，用于回溯完整路径
self.cost = 0.0 # 从起点到当前节点的累积成本（如距离）
self.resource_consumption = [] # 资源消耗向量，如 [时间, 载重]
self.dominated = False # 是否被“支配”（即是否无效）
self.visited_nodes = set() # 已经访问过的节点集合（用于ESPPRC）
算法的核心流程是一个循环，我把它总结为以下几个步骤：
- 初始化：创建起点的初始标签（成本为0，资源消耗为0），放入“未处理标签集合”。
- 选择与扩展：从“未处理集合”中取出一个标签，检查它是否有效（未被支配）。如果有效，就尝试从它的
current_node出发，访问所有可行的下一个节点。 - 可行性检查：对于每一个下一个节点，计算延伸后的新资源消耗（如到达时间、总载重）。检查是否满足资源约束（如时间窗、容量限制）。只有全部满足，才能生成新标签。
- 支配规则检查：这是算法的加速关键。新标签生成后，需要和同一节点上已有的所有标签进行比较。如果存在一个旧标签，它在所有资源消耗和成本上都不比新标签差（即旧标签的成本 ≤ 新标签成本，且旧标签的每种资源消耗 ≤ 新标签的对应消耗），那么新标签就被“支配”了，是无效的，可以直接丢弃。反之，如果新标签支配了某个旧标签，则旧标签被标记为
dominated并移除。 - 收集结果：如果被扩展的标签已经到达终点，并且成本为负（在列生成中，我们寻找的是检验数为负的路径），则将其加入“候选路径集合”。
- 循环与终止：重复步骤2-5，直到“未处理标签集合”为空，或找到足够数量的候选路径。
这个过程中，支配规则就像是一个严格的面试官，它确保只有“帕累托最优”的潜力股（即无法在所有维度上同时被其他路径超越的路径）才能进入下一轮。这极大地剪掉了搜索空间的分支。
3. 代码实战：构建Python版SPPRC求解器
理论说得再多，不如一行代码。我们现在就基于上面的思路，用Python实现一个结构清晰、可读性强的SPPRC求解器。我们会使用面向对象的设计，让代码模块化，便于理解和扩展。
首先，我们需要一个数据类来读取和存储Solomon数据集的信息。
3.1 数据处理：解析Solomon数据集
Solomon基准数据集是车辆路径问题领域的经典测试集，包含了客户点的坐标、需求、时间窗等信息。我们写一个工具函数来读取它。
# utils/data_loader.py
import re
import math
import numpy as np
class SolomonData:
def __init__(self):
self.customer_num = 0 # 顾客数量
self.node_num = 0 # 节点总数 (起点+顾客点+终点)
self.vehicle_capacity = 0 # 车辆容量
self.coord_x = [] # 节点x坐标
self.coord_y = [] # 节点y坐标
self.demand = [] # 节点需求（起点终点为0）
self.ready_time = [] # 时间窗开始时间
self.due_time = [] # 时间窗结束时间
self.service_time = [] # 服务时间
self.distance_matrix = None # 距离矩阵
def load_solomon_data(filepath, max_customers=None):
"""
读取Solomon格式的数据文件。
:param filepath: 数据文件路径
:param max_customers: 最大读取的顾客数，用于快速测试，为None时读取全部
:return: SolomonData 对象
"""
data = SolomonData()
with open(filepath, 'r') as f:
lines = f.readlines()
# 解析文件头信息
for i, line in enumerate(lines):
if i == 4: # 第5行通常是车辆数和容量信息
parts = re.split(r'\s+', line.strip())
data.vehicle_capacity = float(parts[2])
elif i >= 9: # 第10行开始是客户数据
if line.strip() == '':
continue
parts = re.split(r'\s+', line.strip())
# 列顺序：ID, XCOORD, YCOORD, DEMAND, READY_TIME, DUE_TIME, SERVICE_TIME
data.coord_x.append(float(parts[1]))
data.coord_y.append(float(parts[2]))
data.demand.append(float(parts[3]))
data.ready_time.append(float(parts[4]))
data.due_time.append(float(parts[5]))
data.service_time.append(float(parts[6]))
# 如果限制了最大顾客数，提前结束读取
if max_customers is not None and len(data.demand) - 1 >= max_customers:
# -1是因为第一个读入的是起点（ID为0）
break
# 确定节点数量：起点 + 实际读取的顾客点 + 终点（我们通常将终点设为与起点相同或最后一个点）
# 这里我们简单地将最后一个点复制一份作为终点，形成开放路径
data.node_num = len(data.coord_x)
data.customer_num = data.node_num - 1 # 假设第一个点是起点
# 构建距离矩阵（欧几里得距离）
data.distance_matrix = np.zeros((data.node_num, data.node_num))
for i in range(data.node_num):
for j in range(data.node_num):
if i != j:
dx = data.coord_x[i] - data.coord_x[j]
dy = data.coord_y[i] - data.coord_y[j]
data.distance_matrix[i][j] = math.sqrt(dx*dx + dy*dy)
else:
data.distance_matrix[i][j] = 0.0
# 为了方便，我们假设旅行时间等于距离
# 在实际VRPTW中，可能需要单独的时间矩阵
return data
3.2 核心引擎：实现标签类与算法
接下来是重头戏——标签算法的实现。我们将定义一个Label类和一个SPPRC_Solver类。
# core/spprc_solver.py
import copy
from typing import List, Optional
class Label:
"""标签类，代表一条从起点到当前节点的部分路径及其状态。"""
__slots__ = ('current_node', 'prev_label_idx', 'cost', 'resources', 'dominated', 'visited_bitmask')
def __init__(self, current_node: int, prev_label_idx: int, cost: float, resources: List[float]):
"""
初始化一个标签。
:param current_node: 当前节点ID
:param prev_label_idx: 前驱标签在列表中的索引
:param cost: 累积成本
:param resources: 资源消耗列表，如 [累计时间, 累计载重]
"""
self.current_node = current_node
self.prev_label_idx = prev_label_idx
self.cost = cost
self.resources = resources[:] # 复制列表，避免引用问题
self.dominated = False
# 使用位掩码高效记录访问过的节点（适用于节点数少于64的情况，否则可用位数组）
self.visited_bitmask = 0
if current_node >= 0:
self.visited_bitmask |= (1 << current_node)
def is_visited(self, node: int) -> bool:
"""检查某个节点是否已被访问过。"""
if node < 0 or node > 63: # 简单示例，假设节点数<64
# 实际项目中对于更多节点，应使用Python的int无限位或位数组
raise ValueError("Node index out of bitmask range for this simple example.")
return (self.visited_bitmask & (1 << node)) != 0
def add_visited(self, node: int):
"""标记某个节点为已访问。"""
if node < 0 or node > 63:
raise ValueError("Node index out of bitmask range.")
self.visited_bitmask |= (1 << node)
class SPPRC_Solver:
"""SPPRC标签校正算法求解器。"""
def __init__(self, data: SolomonData, resource_limits: List[float]):
"""
初始化求解器。
:param data: Solomon数据对象
:param resource_limits: 资源上限列表，如 [最大时间, 最大载重]
"""
self.data = data
self.resource_limits = resource_limits
self.labels: List[Label] = [] # 存储所有生成的标签
# 记录每个节点关联的标签索引（用于快速查找和支配检查）
self.node_to_label_indices: List[List[int]] = [[] for _ in range(data.node_num)]
def _is_feasible_extension(self, from_label: Label, to_node: int) -> Optional[List[float]]:
"""
检查从给定标签扩展到to_node是否可行。
如果可行，返回新的资源消耗向量；否则返回None。
"""
# 1. 检查是否重复访问（对于SPPRC，这一步通常跳过或作为可选项）
# 如果我们求解的是ESPPRC，这里需要强制检查：
# if from_label.is_visited(to_node):
# return None
# 2. 计算新的资源消耗
new_resources = from_label.resources.copy()
# 资源1: 时间
travel_time = self.data.distance_matrix[from_label.current_node][to_node]
arrival_time = new_resources[0] + travel_time
# 如果早于时间窗开始，需要等待
if arrival_time < self.data.ready_time[to_node]:
arrival_time = self.data.ready_time[to_node]
# 检查是否晚于时间窗结束
if arrival_time > self.data.due_time[to_node]:
return None
new_resources[0] = arrival_time + self.data.service_time[to_node] # 加上服务时间
# 资源2: 载重
new_resources[1] += self.data.demand[to_node]
if new_resources[1] > self.resource_limits[1]: # 超过容量限制
return None
# 可以继续添加其他资源约束...
return new_resources
def _dominates(self, label_a: Label, label_b: Label) -> bool:
"""
判断标签A是否支配标签B。
支配条件：A在所有维度上都不比B差，且至少在一个维度上严格更好。
"""
if label_a.current_node != label_b.current_node:
return False
# 比较成本和所有资源
all_not_worse = (label_a.cost <= label_b.cost + 1e-7)
for ra, rb in zip(label_a.resources, label_b.resources):
all_not_worse = all_not_worse and (ra <= rb + 1e-7)
# 必须至少有一个维度严格更好
strictly_better = (label_a.cost < label_b.cost - 1e-7)
for ra, rb in zip(label_a.resources, label_b.resources):
strictly_better = strictly_better or (ra < rb - 1e-7)
return all_not_worse and strictly_better
def _clean_dominated_labels_at_node(self, node: int):
"""清理某个节点上被支配的标签。"""
label_indices = self.node_to_label_indices[node]
non_dominated = []
for i in range(len(label_indices)):
idx_i = label_indices[i]
label_i = self.labels[idx_i]
if label_i.dominated:
continue
dominated = False
# 与列表中其他未支配标签比较
for j in range(len(label_indices)):
if i == j:
continue
idx_j = label_indices[j]
label_j = self.labels[idx_j]
if label_j.dominated:
continue
if self._dominates(label_j, label_i):
label_i.dominated = True
dominated = True
break
if not dominated:
non_dominated.append(idx_i)
self.node_to_label_indices[node] = non_dominated
def solve(self, start_node: int = 0, end_node: Optional[int] = None, max_paths: int = 10):
"""
求解从start_node到end_node的SPPRC。
:param start_node: 起点索引
:param end_node: 终点索引，默认为最后一个节点
:param max_paths: 希望找到的最大路径数量（用于提前终止）
:return: 找到的优化路径列表，每个元素为(成本, 节点列表)
"""
if end_node is None:
end_node = self.data.node_num - 1
# 初始化：创建起点标签
initial_resources = [self.data.ready_time[start_node], 0.0] # 开始时间，载重
initial_label = Label(start_node, -1, 0.0, initial_resources)
self.labels.append(initial_label)
self.node_to_label_indices[start_node].append(0)
# 未处理标签集合（这里用队列，但按特定顺序弹出可能更高效）
unprocessed_indices = [0]
optimal_paths = []
while unprocessed_indices and len(optimal_paths) < max_paths:
label_idx = unprocessed_indices.pop(0)
current_label = self.labels[label_idx]
if current_label.dominated:
continue
# 如果当前标签已经到达终点，记录为候选路径
if current_label.current_node == end_node:
# 在列生成中，我们通常寻找成本为负的路径
if current_label.cost < -1e-7:
path = self._backtrack_path(label_idx)
optimal_paths.append((current_label.cost, path))
continue
# 扩展当前标签
current_node = current_label.current_node
for next_node in range(self.data.node_num):
if next_node == current_node:
continue
# 检查扩展可行性
new_resources = self._is_feasible_extension(current_label, next_node)
if new_resources is None:
continue
# 计算新成本（这里成本就是距离，在列生成中会是 reduced cost）
new_cost = current_label.cost + self.data.distance_matrix[current_node][next_node]
# 创建新标签
new_label = Label(next_node, label_idx, new_cost, new_resources)
# 复制访问位掩码（如果实现ESPPRC，需要设置 visited）
new_label.visited_bitmask = current_label.visited_bitmask
new_label.add_visited(next_node) # 标记新节点为已访问
# 检查新标签是否被该节点现有标签支配
dominated = False
for existing_idx in self.node_to_label_indices[next_node]:
existing_label = self.labels[existing_idx]
if existing_label.dominated:
continue
if self._dominates(existing_label, new_label):
dominated = True
break
if not dominated:
# 新标签不被支配，加入列表
new_label_idx = len(self.labels)
self.labels.append(new_label)
self.node_to_label_indices[next_node].append(new_label_idx)
unprocessed_indices.append(new_label_idx)
# 新标签可能支配该节点原有的标签
to_remove = []
for existing_idx in self.node_to_label_indices[next_node]:
if existing_idx == new_label_idx:
continue
existing_label = self.labels[existing_idx]
if existing_label.dominated:
continue
if self._dominates(new_label, existing_label):
existing_label.dominated = True
to_remove.append(existing_idx)
# 清理被支配的标签索引
self.node_to_label_indices[next_node] = [
idx for idx in self.node_to_label_indices[next_node]
if idx not in to_remove
]
# 按成本排序返回路径
optimal_paths.sort(key=lambda x: x[0])
return optimal_paths
def _backtrack_path(self, label_idx: int) -> List[int]:
"""从终点标签回溯，重建完整路径。"""
path = []
while label_idx != -1:
label = self.labels[label_idx]
path.append(label.current_node)
label_idx = label.prev_label_idx
return path[::-1] # 反转得到从起点到终点的顺序
3.3 可视化与结果分析
算法跑出了结果，我们当然要直观地看看这些路径长什么样。用matplotlib可以轻松实现。
# visualization/plotter.py
import matplotlib.pyplot as plt
def plot_solution(data: SolomonData, path: List[int], title: str = "SPPRC Solution"):
"""
绘制求解出的路径。
:param data: 包含坐标的数据对象
:param path: 节点索引列表，表示路径
:param title: 图表标题
"""
plt.figure(figsize=(10, 8))
# 绘制所有节点
plt.scatter(data.coord_x[0], data.coord_y[0], c='green', s=200, marker='s', edgecolors='black', label='Depot', zorder=5)
plt.scatter(data.coord_x[1:-1], data.coord_y[1:-1], c='lightblue', s=100, edgecolors='black', label='Customers', zorder=4)
if data.node_num > 1:
plt.scatter(data.coord_x[-1], data.coord_y[-1], c='red', s=200, marker='s', edgecolors='black', label='End Depot', zorder=5)
# 绘制路径
for i in range(len(path) - 1):
start_node = path[i]
end_node = path[i + 1]
plt.plot([data.coord_x[start_node], data.coord_x[end_node]],
[data.coord_y[start_node], data.coord_y[end_node]],
'b-', linewidth=1.5, alpha=0.7, zorder=3)
# 在连线中间加个箭头表示方向
mid_x = (data.coord_x[start_node] + data.coord_x[end_node]) / 2
mid_y = (data.coord_y[start_node] + data.coord_y[end_node]) / 2
plt.annotate('', xy=(data.coord_x[end_node], data.coord_y[end_node]),
xytext=(mid_x, mid_y),
arrowprops=dict(arrowstyle='->', color='blue', lw=1, alpha=0.7))
# 标注节点ID
for i in range(data.node_num):
plt.annotate(str(i), (data.coord_x[i], data.coord_y[i]),
textcoords="offset points", xytext=(0,5), ha='center', fontsize=8)
plt.xlabel("X Coordinate")
plt.ylabel("Y Coordinate")
plt.title(title)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.axis('equal')
plt.tight_layout()
plt.show()
def print_path_stats(data: SolomonData, path: List[int], cost: float):
"""打印路径的详细统计信息。"""
total_distance = 0.0
total_load = 0.0
current_time = data.ready_time[path[0]]
print(f"Path: {' -> '.join(map(str, path))}")
print(f"Total Cost (Distance): {cost:.2f}")
print("\nStep-by-step breakdown:")
for i in range(len(path) - 1):
from_node = path[i]
to_node = path[i + 1]
dist = data.distance_matrix[from_node][to_node]
total_distance += dist
total_load += data.demand[to_node]
arrival = current_time + dist
wait = max(0.0, data.ready_time[to_node] - arrival)
service_start = arrival + wait
service_end = service_start + data.service_time[to_node]
print(f" {from_node} -> {to_node}: Dist={dist:.2f}, "
f"Arrive={arrival:.2f}, Wait={wait:.2f}, "
f"Service=[{service_start:.2f}, {service_end:.2f}], "
f"CumulLoad={total_load:.2f}")
current_time = service_end
print(f"\nSummary: Total Distance = {total_distance:.2f}, Total Load = {total_load:.2f}")
print(f"Feasibility Check: Load <= Capacity? {total_load <= data.vehicle_capacity}")
4. 实战演练：在Solomon数据集上运行
现在，让我们把所有的部件组装起来，在一个实际的Solomon数据集上运行我们的算法。我们选用经典的c101.txt数据集，并先只取前20个客户点进行快速测试。
# main.py
import time
from utils.data_loader import load_solomon_data
from core.spprc_solver import SPPRC_Solver
from visualization.plotter import plot_solution, print_path_stats
def main():
# 1. 加载数据
data_file = "./data/c101.txt"
customer_limit = 20 # 为了快速演示，先处理20个客户
print(f"Loading Solomon data from {data_file} (first {customer_limit} customers)...")
data = load_solomon_data(data_file, max_customers=customer_limit)
print(f"Data loaded: {data.node_num} nodes, capacity={data.vehicle_capacity}")
# 2. 设置资源约束
# 假设资源1是时间（最大为最后一个时间窗的结束时间），资源2是载重（车辆容量）
max_time = max(data.due_time) # 一个宽松的时间上限
resource_limits = [max_time, data.vehicle_capacity]
# 3. 创建求解器并求解
solver = SPPRC_Solver(data, resource_limits)
start_node = 0
end_node = data.node_num - 1 # 最后一个节点作为终点
print("\nStarting SPPRC labeling algorithm...")
start_time = time.time()
optimal_paths = solver.solve(start_node, end_node, max_paths=5)
elapsed_time = time.time() - start_time
print(f"Algorithm finished in {elapsed_time:.3f} seconds.")
print(f"Number of labels generated: {len(solver.labels)}")
# 4. 输出结果
if optimal_paths:
best_cost, best_path = optimal_paths[0]
print(f"\n*** Best Path Found ***")
print_path_stats(data, best_path, best_cost)
plot_solution(data, best_path, title=f"SPPRC Solution (Cost: {best_cost:.2f})")
# 如果有其他优化路径，也显示一下
if len(optimal_paths) > 1:
print(f"\nOther feasible paths found:")
for i, (cost, path) in enumerate(optimal_paths[1:], start=2):
print(f" Path {i}: Cost={cost:.2f}, Nodes={len(path)}, Route={path}")
else:
print("No feasible path found under the given constraints.")
if __name__ == "__main__":
main()
运行这段代码，你会看到终端输出详细的路径信息和统计，同时弹出一个窗口展示路径的可视化结果。从起点（绿色方块）出发，经过一系列客户点（蓝色圆点），最终到达终点（红色方块）。连线上的箭头指示了行驶方向。
5. 性能优化与进阶思考
我们实现了一个正确但基础的版本。在处理更大规模的问题时（比如100个客户点），你可能会发现计算时间变长了。这时，可以考虑以下优化方向：
- 更高效的支配规则实现：我们现在的支配检查是O(n²)的。可以尝试按资源字典序对标签排序，这样可能只需与相邻标签比较。
- 双向标签算法：同时从起点和终点生成标签，在中间汇合，可以大幅减少搜索空间。
- 启发式支配规则：使用一些松弛的支配条件（如NG-path松弛）来更激进地剪枝，虽然可能损失最优性，但能极大加速，适合在列生成中快速获得负检验数路径。
- 资源离散化：将连续的资源（如时间）离散化为多个状态，用动态规划的思想来合并状态，适用于资源约束范围不大的情况。
- 并行化：标签的扩展和支配检查可以并行处理。
从SPPRC到ESPPRC，核心的修改在于强制要求路径是基本的（无重复节点）。这需要在标签中维护已访问节点集合（我们代码中的visited_bitmask已经为此做了准备），并在扩展时严格检查。虽然问题变难了，但算法的基本骨架不变。
把这个项目作为你运筹优化算法实践的起点，试着去修改代码，实现ESPPRC，或者尝试上述的某种优化技巧。真正的理解，来自于动手实现和不断调试的过程中。
这里是“一人公司”的成长家园。我们提供从产品曝光、技术变现到法律财税的全栈内容，并连接云服务、办公空间等稀缺资源，助你专注创造，无忧运营。
更多推荐
11
0- 0
已为社区贡献1条内容
所有评论(0)