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
 $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_{j_1\to j_2}$: Real number representing the distance from city $j_1$ to city $j_2$ (city $j = 0$ denotes depot) ($j_1, j_2 \in \{0,1,\dots,N \}$)
Variables¶
 $x_{i,j,k} \in \{0, 1\} \quad (i \in \{0, \dots , L+1\}, j \in \{0, \dots, N\}, k \in \{0, \dots,
K1\})$
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 $j = 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_{0,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ $$ $$ x_{L+1,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 \{0, \dots, K  1\}) $$ 
City $j$ (excluding the depot) is visited exactly once by any vehicle
$$ \sum_{i=1}^L \sum_{k=0}^{K1} 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 \{0, \dots, K  1\}) $$
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=0}^{K1}\sum_{i=0}^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=0$)
$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=1$)
$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=2$)
$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 \{0, \dots, K1\}) $$
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=0}^{K1} 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. If the solution is supposed to stay in the depot except at the start of delivery ($i = 0$) or at the end of delivery ($i=L+1$), we ignore the situation $i$.
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 15, and the transportation demand for each city is randomly generated.
Also,
adjust capacity
to be a solvable condition.
import numpy as np
from geopy.distance import geodesic
# The number of vehicles
nvehicle = 2
# The number of cities to visit
ncity = 15
avg_cities_per_vehicle = ncity // nvehicle
# Fix the random seed
seed = 12345
rng = np.random.default_rng(seed)
# Randomly determin the delivery demand (weight) in each city
demand = rng.integers(1, 100, size=ncity)
demand_max = np.max(demand)
demand_mean = demand.mean()
# Set an appropriate vehicle payload Q for the above problem setting
capacity = int(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 = [
(
rng.uniform(lon_range[0], lon_range[1]),
rng.uniform(lat_range[0], lat_range[1]),
)
for i in range(ncity + 1)
]
# Coordinate distance matrix between two cities
distance_matrix = np.array(
[
[geodesic(coord_i[::1], coord_j[::1]).m for coord_j in ind2coord]
for coord_i in ind2coord
]
)
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: list[tuple], title: str, best_tour: dict = dict()):
l = len(coord)
center = [
sum(lat for _, lat in coord) / l,
sum(lon for lon, _ in coord) / 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( # type: ignore
locations=[coord[city][::1] for city in tour], color=_color, weight=3
).add_to(m)
else:
for k, node in enumerate(coord):
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)) # type: ignore
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 binary polynomial formulation¶
Next, we define the necessary decision variables, using Amplify's VariableGenerator
. For
a
total of $N+1$ cities and depots, we will visit the depot at the beginning, a maximum of $L$ cities,
and
the depot again at the end with $K$ vehicles, so we define the binary 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 VariableGenerator
gen = VariableGenerator()
# Maximum number of cities one vehicle can visit based on payload capacity
max_tourable_cities = upperbound_of_tour(capacity, demand)
x = gen.array("Binary", shape=(max_tourable_cities + 2, ncity + 1, nvehicle))
As mentioned in the section, Formulation, the constraints are:
 All vehicles start from and end at the depot (city $j = 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_{0,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ & x_{L+1,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 \{0, \dots, K1\})\\ (3)\quad & \sum_{i=1}^L \sum_{k=0}^{K1} 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 \{0, \dots, K1\})\\ \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_{0,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ x_{L+1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \end{align*} $$
This constraint is achieved by setting values to an array of variables. This reduces the number of variables in the problem.
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 \{0, \dots, K1\}) $$
This constraint equation can be implemented by one_hot
function.
When axis
argument is set, firstly the variable array is summed up along the axis, and
then
the equality constraints for those summations are generated.
In this case, we set axis=1
to sum over $j$.
from amplify import one_hot
one_trip_constraints = one_hot(x[1:1, :, :], axis=1)
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=0}^{K1} x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}) $$
This also can be implemented by the one_hot
function. In the same way,
axis=(0, 2)
is specified to sum over $i$ and $k$ as follows:
one_visit_constraints = one_hot(x[1:1, 1:, :], axis=(0, 2))
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 \{0, \dots, K1\}) $$
This can be implemented using less_equal
and einsum
.
First, we sum $w_j x_{i,j,k}$ over $j$ by using einsum
function.
For the details of einsum
, refer to this page.
from amplify import einsum
weight_sums = einsum("j,ijk>ik", demand, x[1:1, 1:, :])
Then, we construct the inequality constraints for each $k$ by using less_equal
function.
As well as one_hot
described above, less_equal
also has axis
argument.
In this case, axis=0
and 1
of weight_sums
correspond to $i$ and
$k$, respectively (note that the axis of $j$ is already summed up), and therefore we can sum over $i$
with
axis=0
.
Also, although not necessary in the present problem setup, we specify penalty_formation
argument as "Relaxation"
of 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/amplify/v1/penalty.html#ineqpenalty for further information.
from amplify import less_equal, ConstraintList
capacity_constraints: ConstraintList = less_equal(
weight_sums, # type: ignore
capacity,
axis=0,
penalty_formulation="Relaxation",
)
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=0}^{K1}\sum_{i=0}^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. To represent
$x_{i+1,j_2,k}$,
we take a slice of variable array x
shifted by one along $i$.
from amplify import Poly, 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: Poly = einsum("pq,ipk,iqk>", distance_matrix, x[:1], x[1:]) # type: ignore
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)
.
from amplify import Model
constraints = one_trip_constraints + one_visit_constraints + capacity_constraints
constraints *= np.max(distance_matrix) # set weight
model = Model(objective, constraints)
Client configuration¶
Now, we create a client for the combinatorial optimization solver (Fixstars Amplify Annealing Engine).
from amplify import FixstarsClient
from datetime import timedelta
client = FixstarsClient()
# Timeout is set to be 2 seconds
client.parameters.timeout = timedelta(milliseconds=2000)
# client.token = "Please put your Amplify API Token"
By passing the model and client to solve
function, we can execute the solver based on
the
QUBO model.
from amplify import solve
result = solve(model, client)
if len(result) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
x_values = result.best.values # solution with the lowest objective function value
Postprocessing and solution visualization¶
The next step is to postprocess and visualize the results.
In the following we define the function onehot2sequence
to format the obtained solution
and
the function process_sequence
to perform postprocessing. The
onehot2sequence
function creates a numbered list in visit order from the onehot binary vector returned by `solve 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.
# Converting onehot variable table to a dictionary. key: vehicle index, value: list containing the order of cities visited by each vehicle
def onehot2sequence(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
onehot2sequence
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 evaluate
method of the decision
variable
array x
and assign the result to the variable array.
Then, postprocess the result using onehot2sequence
and process_sequence
that
we defined earlier.
solution = x.evaluate(x_values)
sequence = onehot2sequence(
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(f"Cost: {result.solutions[0].objective}") # 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].objective:.2f}"
plot_solution(ind2coord, title, best_tour)