建模思路¶
容量约束弧路径问题(CARP)的 OptAgent 模型使用列表变量来表示分配给每辆卡车的边序列。这些边是卡车访问并提供服务的边。卡车在其路径上也可能仅经过其他边而不提供服务。
边可以沿两个方向被访问,但其需求只能被满足一次。OptAgent 的 disjoint 约束确保每条边(沿任一方向)至多出现在一个列表变量中,而 contains(edges_sequences, 2i) + contains(edges_sequences, 2i+1) == 1 保证每对相反方向中恰好有一条被服务。
contains(edges_sequences, edge) 判断某个方向的弧是否出现在任意车辆路线中。
我们可以使用需求数组上的 at 算子来访问序列中每条边的需求。我们使用一个 lambda 函数，通过 sum 算子对所有被访问边的需求求和，从而计算每辆卡车的总配送量。该总配送量必须不超过卡车的容量。
类似地，我们使用二维距离矩阵上的 at 算子来访问从一条边到下一条边所行驶的距离。
Python 实现¶
from pathlib import Path
import re
from optagent import OptModel, solve
class CarpInstance:
"""CARP instance parsed from the DIMACS challenge format."""
EDGE_PATTERN = re.compile(
r"\(\s*(\d+)\s*,\s*(\d+)\s*\)\s+coste\s+(\d+)"
r"(?:\s+demanda\s+(\d+))?"
)
def __init__(self, filename):
lines = Path(filename).read_text(encoding="utf-8").splitlines()
nb_nodes = self._header_int(lines, "VERTICES")
self.nb_required_edges = self._header_int(lines, "ARISTAS_REQ")
nb_not_required_edges = self._header_int(lines, "ARISTAS_NOREQ")
self.nb_trucks = self._header_int(lines, "VEHICULOS")
self.truck_capacity = self._header_int(lines, "CAPACIDAD")
depot = self._header_int(lines, "DEPOSITO")
required_edges = self._read_edges(
lines, "LISTA_ARISTAS_REQ", self.nb_required_edges, required=True
)
other_edges = self._read_edges(
lines, "LISTA_ARISTAS_NOREQ", nb_not_required_edges, required=False
)
# Match the Hexaly model: an undirected adjacency matrix with 0 for
# absent edges, followed by the same O(V^2) Dijkstra preprocessing.
graph = [[0] * (nb_nodes + 1) for _ in range(nb_nodes + 1)]
for origin, destination, cost, _ in required_edges + other_edges:
graph[origin][destination] = cost
graph[destination][origin] = cost
self.costs_data = []
self.demands_data = []
self.origins_data = []
self.destinations_data = []
for origin, destination, cost, demand in required_edges:
self.costs_data.extend((cost, cost))
self.demands_data.extend((demand, demand))
self.origins_data.extend((origin, destination))
self.destinations_data.extend((destination, origin))
relevant_nodes = set(self.origins_data + self.destinations_data + [depot])
shortest_paths = {
node: self._shortest_paths(node, graph) for node in relevant_nodes
}
self.edges_dist_data = [
[
shortest_paths[self.destinations_data[i]][self.origins_data[j]]
for j in range(2 * self.nb_required_edges)
]
for i in range(2 * self.nb_required_edges)
]
self.dist_from_depot_data = [
shortest_paths[depot][origin] for origin in self.origins_data
]
self.dist_to_depot_data = [
shortest_paths[destination][depot] for destination in self.destinations_data
]
@staticmethod
def _header_int(lines, key):
return int(
next(
line.split(":", 1)[1] for line in lines if line.strip().startswith(key)
)
)
@classmethod
def _read_edges(cls, lines, section, count, required):
if count == 0:
return []
start = (
next(
index
for index, line in enumerate(lines)
if line.strip().startswith(section)
)
+ 1
)
edges = []
for line in lines[start : start + count]:
match = cls.EDGE_PATTERN.search(line)
if match is None:
raise ValueError(f"Invalid edge line: {line}")
origin, destination, cost = map(int, match.group(1, 2, 3))
demand = int(match.group(4)) if required else 0
edges.append((origin, destination, cost, demand))
return edges
@staticmethod
def _shortest_paths(start, graph):
distances = [float("inf")] * len(graph)
visited = [False] * len(graph)
distances[start] = 0
for _ in range(len(graph)):
unvisited = (
index for index, distance in enumerate(distances)
if not visited[index]
)
node = min(unvisited, key=distances.__getitem__, default=None)
if node is None or distances[node] == float("inf"):
break
visited[node] = True
for neighbor, edge_cost in enumerate(graph[node]):
if edge_cost and not visited[neighbor]:
candidate = distances[node] + edge_cost
if candidate < distances[neighbor]:
distances[neighbor] = candidate
return distances
def build_carp_model(
nb_required_edges,
nb_trucks,
truck_capacity,
costs_data,
demands_data,
edges_dist_data,
dist_from_depot_data,
dist_to_depot_data,
):
"""Construct the CARP model with directed-edge list variables."""
model = OptModel()
# Seed the search with a capacity-feasible assignment, as in the
# original Hexaly model's initial route construction.
default_sequences = [[] for _ in range(nb_trucks)]
default_loads = [0] * nb_trucks
edge_order = sorted(
range(nb_required_edges),
key=lambda edge: demands_data[2 * edge],
reverse=True,
)
for edge in edge_order:
demand = demands_data[2 * edge]
candidates = [
truck
for truck in range(nb_trucks)
if default_loads[truck] + demand <= truck_capacity
]
if not candidates:
raise ValueError("Unable to build a capacity-feasible initial assignment")
# Best-fit keeps the remaining capacity usable for tight instances.
truck = min(
candidates,
key=lambda candidate: truck_capacity
- default_loads[candidate]
- demand,
)
default_sequences[truck].append(2 * edge)
default_loads[truck] += demand
# Each item is one orientation of a required edge: 2*i is direct, 2*i+1 reverse.
edges_sequences_vars = [
model.list(
2 * nb_required_edges,
default=tuple(default_sequences[truck]),
)
for truck in range(nb_trucks)
]
edges_sequences = model.array(edges_sequences_vars)
model.constraint(model.disjoint(edges_sequences))
for edge in range(nb_required_edges):
model.constraint(
model.contains(edges_sequences, 2 * edge)
+ model.contains(edges_sequences, 2 * edge + 1)
== 1,
)
costs_array = model.array(costs_data)
demands_array = model.array(demands_data)
dist_from_depot_array = model.array(dist_from_depot_data)
dist_to_depot_array = model.array(dist_to_depot_data)
edges_dist_array = model.array(edges_dist_data)
route_distances = []
for k in range(nb_trucks):
sequence = edges_sequences_vars[k]
c = model.count(sequence)
# Truck capacity.
demand_lambda = model.lambda_function(lambda edge: demands_array[edge])
route_quantity = model.sum(sequence, demand_lambda)
model.constraint(route_quantity <= truck_capacity)
# Service cost plus travel from the preceding serviced directed edge.
dist_lambda = model.lambda_function(
lambda i: costs_array[sequence[i]]
+ edges_dist_array[sequence[i - 1], sequence[i]]
)
route_dist = model.sum(model.range(1, c), dist_lambda) + model.iif(
c > 0,
costs_array[sequence[0]]
+ dist_from_depot_array[sequence[0]]
+ dist_to_depot_array[sequence[c - 1]],
0,
)
route_distances.append(route_dist)
total_distance = model.sum(*route_distances)
model.minimize(total_distance)
return model, edges_sequences_vars, total_distance
def main(input_file, output_file=None, time_limit=10):
instance = CarpInstance(input_file)
model, edges_sequences_vars, total_distance = build_carp_model(
nb_required_edges=instance.nb_required_edges,
nb_trucks=instance.nb_trucks,
truck_capacity=instance.truck_capacity,
costs_data=instance.costs_data,
demands_data=instance.demands_data,
edges_dist_data=instance.edges_dist_data,
dist_from_depot_data=instance.dist_from_depot_data,
dist_to_depot_data=instance.dist_to_depot_data,
)
solution = solve(model, time_limit_s=float(time_limit))
if not solution.feasible:
print(
f"Required edges = {instance.nb_required_edges}; Trucks = {instance.nb_trucks}; "
f"Status = {solution.feasible}; no feasible route was found"
)
return solution
result_values = {'total_distance': total_distance.value, **{f'truck_{truck}': sequence.value for truck, sequence in enumerate(edges_sequences_vars)}}
lines = [
f"Required edges = {instance.nb_required_edges}; Trucks = {instance.nb_trucks}; "
f"Total distance = {result_values['total_distance']}; "
f"Status = {solution.feasible}"
]
for truck in range(instance.nb_trucks):
sequence = result_values[f"truck_{truck}"]
if sequence:
serviced_edges = [
(instance.origins_data[edge], instance.destinations_data[edge])
for edge in sequence
]
lines.append(f"Truck {truck + 1}: {serviced_edges}")
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 可以控制每个实例的求解时间。
from pathlib import Path
INSTANCE_DIR = Path.cwd() / "instances"
print("Instances:", INSTANCE_DIR)
solution_egl_e1_a = main(INSTANCE_DIR / "egl-e1-A.dat", time_limit=1)
solution_egl_e1_b = main(INSTANCE_DIR / "egl-e1-B.dat", time_limit=1)
solution_egl_e1_c = main(INSTANCE_DIR / "egl-e1-C.dat", time_limit=1)