Capacitated Vehicle Routing Problem (CVRP)¶
In this tutorial, we consider the capacitated vehicle routing problem (CVRP), a type of vehicle routing problem (VRP).
Specific applications include
 Delivery scheduling in the postal and other transportation industries
 Scheduling in garbage collection and street cleaning
A vehicle routing problem is a problem of determining efficient delivery routes from a depot to multiple cities. More specifically, the problem determines the allocation of delivery vehicles and cities, and the order in which the cities are visited, to minimize the total distance traveled by the delivery vehicles, the total cost, and so on. The routing problem can be interpreted as a generalization of the traveling salesperson problem.
The capacitated vehicle routing problem addressed in this tutorial is the above routing problem with an additional load constraint on each vehicle. In other words, each delivery vehicle must satisfy the load constraints while making deliveries  they cannot carry loads heavier than their payload. This problem can be interpreted as a "traveling salesperson problem + knapsack problem".
Here, we consider the case where there is only one delivery location (depot), and the demand in each city and the vehicle capacity take only integer values.
Formulation¶
First, define the constants and variables needed for the formulation.
Constants¶
 $N$: Number of cities
 $L$: Maximum number of cities that can be visited without exceeding the carrying capacity + 2 (origin + destination)
 $K$: Number of vehicles
 $Q$: Maximum payload of each vehicle, integer
 $W_i$: Integer representing the weight of the load to be delivered to city $i$.
 $D_{i\to j}$: Real number representing the distance from city $i$ to city $j$ (city $0$ = depot) ($i, j \in \{0,1,\dots,N \}$)
Variables¶

$x_{i,j,k} \in \{0, 1\} \quad (i \in \{1, \dots , L\}, j \in \{0, \dots, N\}, k \in \{1, \dots, K\})$
Binary variable indicating whether vehicle $k$ chooses city $j$ as $i$th city to visit ($x_{i,j,k}=1\Leftrightarrow$ vehicle $k$ chooses city $j$ as $i$th city to visit)We then consider the constraints that the variables must satisfy.
Constraints¶
We need the following four constraints.
 All vehicles start from and end at the depot (city $0$)
 Vehicle $k$ visits only one $i$th location
 City $j$ (excluding the depot) is visited exactly once by any vehicle
 The total weight of the loads carried by vehicle $k$ is less than or equal to $Q$
These conditions can be expressed as follows:

All vehicles start from and end at the depot (city $0$)
$$ x_{1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ $$ $$ x_{L,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} $$ 
Vehicle $k$ visits only one $i$th location
$$ \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\}) $$ 
City $j$ (excluding the depot) is visited exactly once by any vehicle
$$ \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}) $$ 
The total weight of the loads carried by the vehicle $k$ is less than or equal to $Q$
$$ \sum_{i=1}^L \sum_{j=1}^N w_j x_{i,j,k} \leq Q \quad (k \in \{1, \dots, K\}) $$
Objective Function¶
In this CVRP, we minimize the total distance traveled by the vehicle. The total distance traveled by a vehicle can be expressed in the following equation.
$$ \sum_{k=1}^K\sum_{i=1}^L\sum_{j_1=0}^N\sum_{j_2=0}^N D_{j_1\to j_2}x_{i,j_1,k}x_{i+1,j_2,k} $$
Explanation of Constraint Expressions¶
This section provides a more intuitive explanation using concrete examples.
The following is an example of a variable assignment that satisfies the constraint conditions, summarized in a table for each vehicle. The horizontal axis represents the alphabetical name of the city (the depot is the starting point), and the vertical axis represents the order of the cities to be visited (0 means departing from the depot, 9 means arriving at the depot at the end). If a square in the table is $1$, it means that the vehicle will visit the city.
 Vehicle 1 ($k=1$)
$x_{i,j,1}$  depot  A  B  C  D  E  F  G  H  I 

0 (start)  1  0  0  0  0  0  0  0  0  0 
1  0  1  0  0  0  0  0  0  0  0 
2  0  0  0  1  0  0  0  0  0  0 
3  0  0  0  0  0  1  0  0  0  0 
4  0  0  0  0  0  0  1  0  0  0 
5  1  0  0  0  0  0  0  0  0  0 
6  1  0  0  0  0  0  0  0  0  0 
7  1  0  0  0  0  0  0  0  0  0 
8  1  0  0  0  0  0  0  0  0  0 
9 (end)  1  0  0  0  0  0  0  0  0  0 
 Vehicle 2 ($k=2$)
$x_{i,j,2}$  depot  A  B  C  D  E  F  G  H  I 

0 (start)  1  0  0  0  0  0  0  0  0  0 
1  1  0  0  0  0  0  0  0  0  0 
2  1  0  0  0  0  0  0  0  0  0 
3  1  0  0  0  0  0  0  0  0  0 
4  0  0  1  0  0  0  0  0  0  0 
5  0  0  0  0  0  0  0  1  0  0 
6  0  0  0  0  0  0  0  0  0  1 
7  1  0  0  0  0  0  0  0  0  0 
8  1  0  0  0  0  0  0  0  0  0 
9 (end)  1  0  0  0  0  0  0  0  0  0 
 Vehicle 3 ($k=3$)
$x_{i,j,3}$  depot  A  B  C  D  E  F  G  H  I 

0 (start)  1  0  0  0  0  0  0  0  0  0 
1  0  0  0  0  1  0  0  0  0  0 
2  0  0  0  0  0  0  0  0  1  0 
3  1  0  0  0  0  0  0  0  0  0 
4  1  0  0  0  0  0  0  0  0  0 
5  1  0  0  0  0  0  0  0  0  0 
6  1  0  0  0  0  0  0  0  0  0 
7  1  0  0  0  0  0  0  0  0  0 
8  1  0  0  0  0  0  0  0  0  0 
9 (end)  1  0  0  0  0  0  0  0  0  0 
The abovementioned constraint "2: Vehicle $k$ visits only one $i$th location" corresponds to a single "$1$" appearing in each row (horizontal line) of the above decision variable table and can be expressed by the following equation:
$$ \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\}) $$
Also, the constraint "3: City $j$ (excluding the depot) is visited exactly once by any vehicle" corresponds to a single "$1$" appearing in each column (vertical line) except the depot for all vehicles combined, and can be expressed as follows:
$$ \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}) $$
While these constraints ensure that each city is visited exactly once by any vehicle, the formulation allows that the delivery depot to be visited any number of times. Therefore, postprocessing is performed so that each vehicle visits a depot only twice, once at the start of the delivery and once at the end of the delivery. Specifically, if $1$ appears in the "depot" column except for rows 0 and 9 in the table, it is ignored (considered as $0$). The travel distance for the route obtained after this postprocessing is equal to or shorter than that before the postprocessing, which can be proven by using the triangular inequality.
Implementation¶
We will now implement the above formulation and constraints using Fixstars Amplify.
Problem definition and visualization¶
First, let us define the variables needed to optimize the present problem. The number of vehicles is
set
to 2, the number of cities to 10, and the transportation demand for each city is randomly generated.
Also,
adjust capacity
to be a solvable condition. The dimension
is [number of
cities +
number of depots = number of cities + 1].
import numpy as np
from geopy.distance import geodesic
# The number of vehicles
nvehicle = 2
# The number of cities to visit
ncity = 10
dimension = ncity + 1
avg_cities_per_vehicle = dimension // nvehicle
# Randomly determin the delivery demand (weight) in each city
demand = np.random.randint(1, 100, size=dimension)
demand[0] = 0
demand_max = max(demand)
demand_mean = sum(demand) / len(demand)
# Set an appropriate vehicle payload Q for the above problem setting
capacity = demand_max + int(demand_mean) * avg_cities_per_vehicle
Next, the coordinates of each city are randomly generated. Generate $D_{i \to j}$ by computing the coordinate distances between all cities $i,j$.
# Ranges of the coordinate
lat_range = [51, 53]
lon_range = [2, 0]
# Coordinates of each city
ind2coord = {
i: (
np.random.uniform(lon_range[0], lon_range[1]),
np.random.uniform(lat_range[0], lat_range[1]),
)
for i in range(dimension)
}
# Coordinate distance matrix between two cities
distance_matrix = np.array(
[
[geodesic(coord_i[::1], coord_j[::1]).m for coord_j in ind2coord.values()]
for coord_i in ind2coord.values()
]
)
Create a list of colors to be used when plotting the problem.
_colors = [
"green",
"orange",
"blue",
"red",
"purple",
"pink",
"darkblue",
"cadetblue",
"darkred",
"lightred",
"darkgreen",
"lightgreen",
"lightblue",
"darkpurple",
]
Define the function plot_solution
to plot the coordinates of the cities.
import folium
def plot_solution(coord: dict, title: str, best_tour: dict = dict()):
l = len(coord)
center = [
sum(lat for _, lat in coord.values()) / l,
sum(lon for lon, _ in coord.values()) / l,
]
m = folium.Map(center, tiles="OpenStreetMap", zoom_start=7)
folium.Marker(
location=coord[0][::1],
popup=f"depot",
icon=folium.Icon(icon="car", prefix="fa"),
).add_to(m)
_color = _colors[1]
if best_tour:
for k, tour in best_tour.items():
_color = _colors[k % len(_colors)]
for city in tour:
if city == 0:
continue
folium.Marker(
location=coord[city][::1],
popup=f"person{k}",
icon=folium.Icon(
icon="user", prefix="fa", color="white", icon_color=_color
),
).add_to(m)
folium.vector_layers.PolyLine(
locations=[coord[city][::1] for city in tour], color=_color, weight=3
).add_to(m)
else:
for k, node in coord.items():
if k == 0:
continue
folium.Marker(
location=node[::1],
popup=f"customer{k}",
icon=folium.Icon(
icon="user", prefix="fa", color="white", icon_color=_color
),
).add_to(m)
title = f"<h4>{title}</h4>"
m.get_root().html.add_child(folium.Element(title))
return m
The plot_solution
function displays the defined problem. The carshaped pins represent
the
location of the depot (start and end points), and the humanshaped pins represent the points to be
visited.
# England, UK
title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}"
plot_solution(ind2coord, title)
Upper bound of cities to visit¶
To reduce the size of the decision variable, an upper bound on the number of cities that can be visited by a single car is estimated in advance.
The following upperbound_of_tour
function estimates the maximum number of cities that
can be
visited without exceeding the load capacity by the following steps:
 Selecting the city with the lowest delivery demand (weight) among the cities not yet selected, and subtracting that city's demand from the maximum load capacity
 Repeating step 1 until the maximum load capacity becomes negative, and the number of iterations is the upper limit of the number of cities to visit
def upperbound_of_tour(capacity: int, demand: np.ndarray) > int:
max_tourable_cities = 1
for w in sorted(demand):
capacity = w
if capacity >= 0:
max_tourable_cities += 1
else:
return max_tourable_cities
return max_tourable_cities
Constructing QUBO formulation¶
Next, we define the necessary decision variables, using Amplify's BinarySymbolGenerator
.
For
a total of $N+1$ cities and depots, we will visit at most $L$ cities with $K$ vehicles, so we define
the
QUBO variable as a threedimensional array of $L \times (N+1) \times K$ as follows (see the
definitions of
constants in the Formulation section).
from amplify import BinaryPoly, BinarySymbolGenerator
gen = BinarySymbolGenerator()
# Maximum number of cities one vehicle can visit based on payload capacity
max_tourable_cities = upperbound_of_tour(capacity, demand)
x = gen.array(max_tourable_cities, dimension, nvehicle)
As mentioned in the section, Formulation, the constraints are:
 All vehicles start from and end at the depot (city $0$)
 Vehicle $k$ visits only one $i$th location
 City $j$ (excluding the depot) is visited exactly once by any vehicle
 The total weight of the loads carried by vehicle $k$ is less than or equal to $Q$
In other words,
$$ \begin{align*} (1) \quad & x_{1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ & x_{L,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ (2)\quad & \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\})\\ (3)\quad & \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\})\\ (4)\quad & \sum_{i=1}^L \sum_{j=1}^N w_j x_{i,j,k} \leq Q \quad (k \in \{1, \dots, K\})\\ \end{align*}. $$
Note that $x_{i,j,k}$ corresponds to vehicle $k$ visits city $j$ as $i$th location.
We will use the less_equal
function to express inequality constraints and the
equal_to
function to express equality constraints. Also, the coefficient that determines
the
weight of the constraint term is implemented as $\max (d_{ij})$. This at least avoids the situation
where
the objective function is favored over the constraint term.
Constraint (1): all vehicles start from and end at the depot (city 0)¶
The constraint equation is expressed as follows.
$$ \begin{align*} x_{1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ x_{L,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \end{align*} $$
We achieve this constraint by prefixing. Prefixing allows us to reduce the number of variables to be solved. Create a function that directly assigns values as follows.
To prefix the values of the variables, do x[i][j][k] = 1
or x[i][j][k] = 0
.
The
assignment operation to the decision variables must be performed before defining the objective
functions
or constraint equations.
x[0][1:][:] = 0
x[1][1:][:] = 0
x[0][0][:] = 1
x[1][0][:] = 1
Constraint (2): vehicle $k$ visits only one $i$th location¶
This constraint is already expressed by the following equation in the Formulation section.
$$ \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\}) $$
This constraint equation is implemented as follows.
from amplify.constraint import one_hot
from amplify import sum_poly
def make_onetrip_constraints(x) > list:
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
constraints = [
one_hot(sum_poly(dimension, lambda j: x[i][j][k]))
for i in range(max_tourable_cities)
for k in range(nvehicle)
]
return constraints
Constraint (3): city $j$ (excluding the depot) is visited exactly once by any vehicle¶
As already mentioned in the Formulation section, this constraint is expressed as:
$$ \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}).$$
This equation is implemented as follows:
def make_onevisit_constraints(x) > list:
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
constraints = [
one_hot(
sum_poly(
max_tourable_cities, lambda i: sum_poly(nvehicle, lambda k: x[i][j][k])
)
)
for j in range(1, dimension)
]
return constraints
Constraint (4): the total weight of the loads carried by vehicle $k$ is less than or equal to $Q$¶
This constraint is expressed by the following equation.
$$ \sum_{i=1}^L \sum_{j=1}^N w_j x_{i,j,k} \leq Q \quad (k \in \{1, \dots, K\}) $$
This constraint expression can be implemented using less_equal
. Also, although not
necessary
in the present problem setup, we will specify InequalityFormulation.RelaxationQuadra
as
the
method
argument to less_equal
. While this method is not a strict penalty
function for this constraint, it keeps from generating many auxiliary variables internally, can handle
real number coefficients, and thus increases solvability under large problem size conditions.
See https://amplify.fixstars.com/en/docs/constraint.html#constraintinequalityrelaxation for further information.
from amplify.constraint import less_equal
from amplify import InequalityFormulation
def make_capacity_constraints(x, demand: np.ndarray, capacity: int) > list:
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
constraints = [
less_equal(
demand * x[:, :, k],
capacity,
method=InequalityFormulation.RelaxationQuadra,
)
/ capacity
/ capacity
for k in range(nvehicle)
]
return constraints
Objective function¶
The objective this time is to minimize the total distance traveled by the vehicles, and the objective function may be written as follows:
$$ \sum_{k=1}^K\sum_{i=1}^L\sum_{j_1=0}^N\sum_{j_2=0}^N D_{j_1\to j_2}x_{i,j_1,k}x_{i+1,j_2,k} $$
The einsum
function is used to compute the objective function. More information about
the
einsum
function can be found here.
from amplify import einsum
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
xr = x.roll(1, axis=0)
# Total distance of the route
objective = einsum("pq,ipk,iqk", distance_matrix, x, xr)
Based on the preassignmentbased constraint, the three constraint equations, and the objective
function
implemented so far, we can now construct a QUBO model for CVRP as follows. Here, we set the weight for
the
constraint equations as np.max(distance_matrix)
.
constraints1 = make_onetrip_constraints(x)
constraints2 = make_onevisit_constraints(x)
constraints3 = make_capacity_constraints(x, demand, capacity)
constraints = sum(constraints1) + sum(constraints2) + sum(constraints3)
model = objective + constraints * np.max(distance_matrix)
Setting up Amplify client¶
Let us use the Amplify Annealing Engine (AE) client.
from amplify.client import FixstarsClient
client = FixstarsClient()
client.parameters.timeout = 2000 # Timeout is set to be 2 seconds
# client.token = "Please put your Amplify API Token"
Then we will instantiate the Solver
class with the client we just created, and we will
execute the solution based on the QUBO model.
from amplify import Solver
solver = Solver(client)
result = solver.solve(model)
if len(result.solutions) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
x_values = result.solutions[0].values # solution
Postprocessing and solution visualization¶
The next step is to postprocess and visualize the results.
In the following we define the function onehot2_sequence
to format the obtained solution
and
the function process_sequence
to perform postprocessing as described above. The onehot2_sequence
function creates a numbered list
in
visit order from the onehot binary vector returned by Solver
as follows.
Note that the following type hints are compatible with Python 3.9 or later. If you are using an older
environment (Python 3.8), such as in a local environment, run
from __future__ import annotations
additionally.
# from __future__ import annotations # For Python 3.8, uncomment this
# Converting onehot variable table to a dictionary. key: vehicle index, value: list containing the order of cities visited by each vehicle
def onehot2_sequence(x_values) > dict[int, list]:
nvehicle = x_values.shape[2]
sequence = dict()
for k in range(nvehicle):
sequence[k] = np.where(x_values[:, :, k])[1]
return sequence
Create the process_sequence
function, which will return the list created by the
onehot2_sequence
function with the extra visits to the depot removed as described in described. Using the triangular inequality, we can prove that the distance
traveled after postprocessing is less than or equal to the distance traveled before postprocessing.
def process_sequence(sequence: dict[int, list]) > dict[int, list]:
new_seq = dict()
for k, v in sequence.items():
v = np.append(v, v[0])
mask = np.concatenate(([True], np.diff(v) != 0))
new_seq[k] = v[mask]
return new_seq
Now let's extract and visualize the result.
We pass the solution x_values
to the decode
method of the decision variable
array x
and assign the result to the variable array.
Then, postprocess the result using onehot2_sequence
and process_sequence
that
we defined earlier.
solution = x.decode(x_values)
sequence = onehot2_sequence(
solution
) # Converting onehot variable table to a dictionary. key: vehicle index, value: list containing the order of cities visited by each vehicle
best_tour = process_sequence(
sequence
) # Remove extra visits to the depot from the dictionary above.
print(result.solutions[0].energy) # Print objective function value
print(*best_tour.items(), sep="\n") # Print the obtained route
The plot_solution
function, defined earlier in the Implementation section, visualizes
the
postprocessed solution. The optimized route that each vehicle should take is displayed on the map.
title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}, cost={result.solutions[0].energy:.2f}"
plot_solution(ind2coord, title, best_tour)