数据¶
所提供的容量约束车辆路径问题(CVRP)算例来自 Augerat 等人的 Set A 数据集。它们遵循 TSPLib 格式:
- 节点数量由关键字 DIMENSION 指定(其中包括一个配送中心,因此客户数量为节点数减 1)。
- 卡车容量由关键字 CAPACITY 指定。
- 边的类型由 EDGE_WEIGHT_TYPE 指定。请注意,在我们的模型中仅接受 EUC_2D 这一种边类型。
- 在关键字 NODE_COORD_SECTION 之后:每个节点的 ID 以及其 x、y 坐标。
- 在关键字 DEMAND_SECTION 之后:每个节点的 ID 及其需求。
- 配送中心列在关键字 DEPOT_SECTION 之后。请注意,在我们的模型中仅接受一个配送中心。
可用车辆数量等于客户数量。
建模思路¶
容量约束车辆路径问题(CVRP)的 OptAgent 模型使用 列表决策变量。对每辆卡车,我们定义一个列表变量,表示其所访问的客户序列。通过对所有列表施加 partition 约束,我们确保每个客户恰好由一辆卡车服务。
当一辆卡车至少访问一个客户时,它才被车队所使用。借助 count 算子(返回列表中的元素数量),我们可以检查每辆卡车是否被使用,从而计算车队中使用的卡车总数。
我们可以使用需求数组上的 at 算子来访问序列中每个客户的需求。每辆卡车所配送的总量通过一个 lambda 函数计算,该函数将 sum 算子应用于所有被访问的客户。请注意,该 sum 中的项数以及列表的大小在搜索过程中会变化。
从一个客户到下一个客户所行驶的距离同样使用二维距离矩阵上的 at 算子来访问。我们使用另一个 lambda 函数 计算每辆卡车的总行驶距离,该函数对序列中相邻客户之间的距离求和。
最后依次声明两个最小化目标:先最小化使用的车辆数,再最小化总行驶距离。OptAgent 按目标声明顺序执行字典序优化。
Python 实现¶
import math
from pathlib import Path
from optagent import OptModel, solve
def read_cvrp_instance(filename):
"""Read an Augerat-format CVRP instance and return a dict of arrays."""
tokens = Path(filename).read_text(encoding="utf-8").split()
nb_nodes = 0
truck_capacity = 0
depot_x = 0
depot_y = 0
index = 0
while index < len(tokens):
token = tokens[index]
if token == "DIMENSION":
nb_nodes = int(tokens[index + 2])
index += 3
continue
if token == "CAPACITY":
truck_capacity = int(tokens[index + 2])
index += 3
continue
if token == "EDGE_WEIGHT_TYPE":
if tokens[index + 2] != "EUC_2D":
raise ValueError(f"Only EUC_2D is supported, got {tokens[index + 2]}")
index += 3
continue
if token == "NODE_COORD_SECTION":
index += 1
break
index += 1
nb_customers = nb_nodes - 1
customers_x = [0] * nb_customers
customers_y = [0] * nb_customers
for _ in range(nb_nodes):
node_id = int(tokens[index])
index += 1
x = int(tokens[index])
index += 1
y = int(tokens[index])
index += 1
if node_id == 1:
depot_x = x
depot_y = y
else:
customers_x[node_id - 2] = x
customers_y[node_id - 2] = y
while index < len(tokens) and tokens[index] != "DEMAND_SECTION":
index += 1
index += 1
demands = [0] * nb_customers
for _ in range(nb_nodes):
node_id = int(tokens[index])
index += 1
demand = int(tokens[index])
index += 1
if node_id != 1:
demands[node_id - 2] = demand
distance_matrix = compute_distance_matrix(customers_x, customers_y)
distance_depots = compute_distance_depots(
depot_x, depot_y, customers_x, customers_y
)
return {
"nb_customers": nb_customers,
"nb_trucks": nb_customers,
"truck_capacity": truck_capacity,
"distance_matrix": distance_matrix,
"distance_depots": distance_depots,
"demands": demands,
}
def compute_dist(xi, xj, yi, yj):
exact_dist = math.sqrt(math.pow(xi - xj, 2) + math.pow(yi - yj, 2))
return int(math.floor(exact_dist + 0.5))
def compute_distance_matrix(customers_x, customers_y):
nb_customers = len(customers_x)
distance_matrix = [
[None for _ in range(nb_customers)] for _ in range(nb_customers)
]
for i in range(nb_customers):
distance_matrix[i][i] = 0
for j in range(nb_customers):
dist = compute_dist(
customers_x[i], customers_x[j], customers_y[i], customers_y[j]
)
distance_matrix[i][j] = dist
distance_matrix[j][i] = dist
return distance_matrix
def compute_distance_depots(depot_x, depot_y, customers_x, customers_y):
nb_customers = len(customers_x)
distance_depots = [None] * nb_customers
for i in range(nb_customers):
dist = compute_dist(depot_x, customers_x[i], depot_y, customers_y[i])
distance_depots[i] = dist
return distance_depots
def main(instance_file, output_file=None, time_limit=10):
data = read_cvrp_instance(instance_file)
nb_customers = data["nb_customers"]
nb_trucks = data["nb_trucks"]
truck_capacity = data["truck_capacity"]
model = OptModel()
customers_sequences = [
model.list(nb_customers) for k in range(nb_trucks)
]
model.constraint(model.partition(customers_sequences))
demands = model.array(data["demands"])
dist_matrix = model.array(data["distance_matrix"])
dist_depot = model.array(data["distance_depots"])
trucks_used = [model.count(customers_sequences[k]) > 0 for k in range(nb_trucks)]
dist_routes = []
for k in range(nb_trucks):
sequence = customers_sequences[k]
count = model.count(sequence)
demand_lambda = model.lambda_function(lambda customer: demands[customer])
route_quantity = model.sum(sequence, demand_lambda)
model.constraint(route_quantity <= truck_capacity)
dist_lambda = model.lambda_function(
lambda i: dist_matrix[sequence[i - 1], sequence[i]]
)
dist_routes.append(
model.sum(model.range(1, count), dist_lambda)
+ model.iif(
count > 0,
dist_depot[sequence[0]]
+ dist_depot[sequence[count - 1]],
0,
)
)
nb_trucks_used = model.sum(*trucks_used)
total_distance = model.sum(*dist_routes)
# Lexicographic objectives, in the same order as the original Hexaly model.
model.minimize(nb_trucks_used)
model.minimize(total_distance)
solution = solve(model, time_limit_s=float(time_limit))
result_values = {'trucks_used': nb_trucks_used.value, 'total_distance': total_distance.value, **{f'truck_{k}': sequence.value for k, sequence in enumerate(customers_sequences)}}
lines = [
f"Trucks used = {result_values['trucks_used']}; "
f"Total distance = {result_values['total_distance']}; "
f"Status = {solution.feasible}"
]
for truck in range(nb_trucks):
sequence = result_values[f"truck_{truck}"]
if sequence:
customers = " ".join(str(customer + 2) for customer in sequence)
lines.append(f"Truck {truck + 1}: {customers}")
result_text = "\n".join(lines)
print(result_text)
if output_file is not None:
Path(output_file).write_text(result_text + "\n", encoding="utf-8")
return solution
本地运行¶
Notebook 直接调用 main 并显式传入实例路径。以下三个代码格相互独立,可以按需要单独运行;调整 time_limit 可以控制每个实例的求解时间。
INSTANCE_DIR = Path.cwd() / "instances"
print("Instances:", INSTANCE_DIR)
solution_a_n32_k5 = main(
INSTANCE_DIR / "A-n32-k5.vrp",
time_limit=1,
)
solution_a_n33_k5 = main(
INSTANCE_DIR / "A-n33-k5.vrp",
time_limit=10,
)
solution_a_n33_k6 = main(
INSTANCE_DIR / "A-n33-k6.vrp",
time_limit=10,
)