The Vehicle Routing Problem: Exact and Heuristic Solutions
Understand how to solve complex routing problems with Python
The Vehicle Routing Problem (VRP) aims to determine the best set of routes to be performed by a fleet of vehicles to serve a given set of customers. Due to its several applications and challenging combinatorial aspects, it is one of the most studied problems in Operations Research and Numerical Optimization.
The Capacitated Vehicle Routing Problem (CVRP) is one of the most common variants, as it introduces vehicles with limited load capacity and possibly duration/distance constraints. Other usual variants also introduce multiple depots, heterogeneous fleets, pickup and deliveries, and time-window constraints.
Combinatorial aspects of these problems are such that considering a simple set of 15 points, there are 6 × 10¹¹ possible routes connecting them (Dantzig & Ramser, 1959). Therefore, some real-world applications would remain impractical until computational and algorithm research advances over the last couple of decades. Branch-Cut-and-Price algorithms have been able to prove the optimality of CVRP instances with a few hundred customers (Fukasawa et al., 2006; Pecin et al., 2017), and state-of-the-art metaheuristics combined with local search techniques can provide good quality (sometimes optimal) solutions to these instances within a few seconds (Vidal et al., 2012; Vidal, 2022).
Throughout this article, we will introduce the Capacitated Vehicle Routing Problem with load (and duration) constraints and solve it using Mixed-Integer Programming (MIP) and specialized (meta)heuristic algorithms. In the first part, we will use the Python AML Pyomo, with the HiGHS solver, whereas in the second part, we will use the Google OR Tools package.
The reader who is more interested in real-world applications than in theoretical aspects of the problem might quickly skim through the MIP section and pay more interest to the Specialized (Meta)Heuristics and Useful Extensions sections.
Those interested in understanding in detail the MIP formulation but not yet familiar with numerical optimization might find it useful to take a look at my previous stories about Linear Programming and the Branch & Bound method before proceeding with this one.
As usual, you can find the complete code in this git repository.
Now, let us dive in!
Mixed-Integer Programming
The mathematical formulation presented in this section will use the same equations presented by Toth & Vigo (2002) in the model therein referred to as the “three-index vehicle flow formulation”.
Consider a set V of nodes (demands and depot) and a set K of vehicles. We will use lowercase i and j to indicate node indexes and lowercase k to indicate vehicle indexes. As this model is valid for the asymmetric case, let us assume the nodes are part of a complete directed graph G(V, A) with arcs A. In this problem, there is a single depot node indexed by 0 and all vehicles have the same capacity Q. Consider two groups of decision variables:
- x_{i, j, k}: Is a binary variable that indicates an active arc from node i to node j performed by vehicle k.
- y_{i, k}: Is a binary variable that indicates that the demand from node i is fulfilled by vehicle k.
Consider our objective to minimize the cost value associated with active arcs. Total duration or distance are usual examples. Say the cost of traversing the arc i, j is cᵢⱼ. The objective function can be stated as the following.
We also need to include constraints that ensure:
- Each customer i is visited once, therefore has one active arc which starts from it and one that arrives on it.
- If any arc variable indexed by vehicle k goes into one node i or out of it, the demand q of this node is assigned to vehicle k.
- The total demand assigned to a vehicle must not exceed its capacity Q.
- Exactly |K| nodes start at the depot and arrive at the depot.
- There are no subtours… However, the number of subtours is potentially too large to be enumerated from the start. We will dive into more detail about how to do it.
As usual for Python tutorials, let us start our hands-on part by importing the libraries used in this section:
import time
from itertools import cycle
import numpy as np
from scipy.spatial.distance import pdist, squareform
import matplotlib.pyplot as plt
import matplotlib as mpl
import networkx as nx
import pyomo.environ as pyo
from pyomo.contrib.appsi.solvers.highs import HighsAnd now let us instantiate a random problem with N demand nodes. In this example, the depot node is assumed to be the first node (index 0), so we ensure its corresponding demand is also zero.
np.random.seed(42) # Results should be always the same
N = 10
demands = np.random.randint(1, 10, size=N)
demands[0] = 0
capacity = 15
n_vehicles = 4
coordinates = np.random.rand(N, 2)
distances = squareform(pdist(coordinates, metric="euclidean"))
distances = np.round(distances, decimals=4) # avoid numerical errorsThe number of necessary vehicles can be calculated by using a Bin Packing Problem. An example of how to perform it is also included in the complete source code.
There are two approaches for modeling a problem in pyomo: Abstract and Concrete models. In the first approach, the algebraic expressions of the problem are defined before some data values are supplied, whereas, in the second, the model instance is created immediately as its elements are defined. You can find more about these approaches in the library documentation or in the book by Bynum et al. (2021). Throughout this article, we will adopt the Concrete model formulation.
model = pyo.ConcreteModel()Let us instantiate the sets of demand nodes V, arcs A, and vehicles K. Notice that the depot node is included in the set of nodes V as in the original mathematical formulation.
model.V = pyo.Set(initialize=range(len(demands)))
model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j])
model.K = pyo.Set(initialize=range(n_vehicles))Now our parameters for capacity, demand load, and arc costs.
model.Q = pyo.Param(initialize=capacity)
model.q = pyo.Param(model.V, initialize={i: d for (i, d) in enumerate(demands)})
model.c = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A})And the decision variables indicating active arcs in a given vehicle and that a node visited by the vehicle.
model.x = pyo.Var(model.A, model.K, within=pyo.Binary)
model.y = pyo.Var(model.V, model.K, within=pyo.Binary)Before including the constraints, I will create our objective that computes the total cost of active arcs.
model.obj = pyo.Objective(
expr=sum(
model.x[i, j, k] * model.c[i, j]
for (i, j) in model.A
for k in model.K
),
sense=pyo.minimize,
)We also must include constraints previously listed. First let us implement them using the usual Pyomo signature: function(model, *domain).
def arcs_in(model, i):
if i == model.V.first():
return sum(model.x[:, i, :]) == len(model.K)
else:
return sum(model.x[:, i, :]) == 1.0
def arcs_out(model, i):
if i == model.V.first():
return sum(model.x[i, :, :]) == len(model.K)
else:
return sum(model.x[i, :, :]) == 1.0
def vehicle_assignment(model, i, k):
return sum(model.x[:, i, k]) == model.y[i, k]
def comp_vehicle_assignment(model, i, k):
return sum(model.x[i, :, k]) == model.y[i, k]
def capacity_constraint(model, k):
return sum(model.y[i, k] * model.q[i] for i in model.V) <= model.QAnd then incorporate them into our model.
model.arcs_in = pyo.Constraint(model.V, rule=arcs_in)
model.arcs_out = pyo.Constraint(model.V, rule=arcs_out)
model.vehicle_assignment = pyo.Constraint(model.V, model.K, rule=vehicle_assignment)
model.comp_vehicle_assignment = pyo.Constraint(model.V, model.K, rule=comp_vehicle_assignment)
model.capacity_constraint = pyo.Constraint(model.K, rule=capacity_constraint)Notice that I didn’t include the subtour elimination constraints yet. We should consider all possible permutations of the N nodes taking k at a time and these could be prohibitively large to enumerate even for moderate size instances. Alternatively, in our solution procedure, we will recursively include subtour elimination constraints each time a new solution is found if we verify that this solution produces subtours. In some commercial solvers, these are called “lazy constraints” and can be incorporated directly into the solver via callback.
First let us create a function that, given a subtour, all remaining nodes, a node from the subtour, and a vehicle, returns a Pyomo expression corresponding to the mathematical formulation previously stated. Also, let us include a ConstraintList to which we will include new elements as we proceed with the solution.
def subtour_elimination(model, S, Sout, h, k):
nodes_out = sum(model.x[i, j, k] for i in S for j in Sout)
return model.y[h, k] <= nodes_out
model.subtour_elimination = pyo.ConstraintList()We must create some functions that will, given a solution, return subtours created (if they exist). To do so, we will first create a list of active arcs in the model using the find_arcs function. This list will be used to create an incomplete directed graph using the Networkx DiGraph class. And the find_subtours function should return a list of sets of connected components.
def find_arcs(model):
arcs = []
for i, j in model.A:
for k in model.K:
if np.isclose(model.x[i, j, k].value, 1):
arcs.append((i, j))
return arcs
def find_subtours(arcs):
G = nx.DiGraph(arcs)
subtours = list(nx.strongly_connected_components(G))
return subtoursOur goal is to eliminate groups of connected components that do not include the depot node. So in the next step, we will create functions that iterate over the list of sets, and include new constraints if the set of components does not include the depot node. This will use the method add of the ConstraintList class.
def eliminate_subtours(model, subtours):
proceed = False
for S in subtours:
if 0 not in S:
proceed = True
Sout = {i for i in model.V if i not in S}
for h in S:
for k in model.K:
model.subtour_elimination.add(
subtour_elimination(model, S, Sout, h, k)
)
return proceedAnd now we have everything ready to propose a solution procedure that iteratively solves the MIP, verifies if the current solution has subtours, and, if so, includes new constraints to eliminate them. Otherwise, the solution found is optimal.
def solve_step(model, solver):
sol = solver.solve(model)
arcs = find_arcs(model)
subtours = find_subtours(arcs)
time.sleep(0.1)
proceed = eliminate_subtours(model, subtours)
return sol, proceed
def solve(model, solver):
proceed = True
while proceed:
sol, proceed = solve_step(model, solver)
return solNow let us instantiate the solver and solve our model. The Highs solver is available in Pyomo (check the imports) if the highspy package is installed. So make sure to run a pip install highspy .
solver = Highs()
solver.highs_options = {
"log_file": "Highs.log",
"mip_heuristic_effort": 0.2,
"mip_detect_symmetry": True,
"mip_rel_gap": 1e-6,
}
solution = solve(model, solver)One more function to find the tours created and we are ready to plot results.
def find_tours(model):
tours = []
for k in model.K:
node = 0
tours.append([0])
while True:
for j in model.V:
if (node, j) in model.A:
if np.isclose(model.x[node, j, k].value, 1):
node = j
tours[-1].append(node)
break
if node == 0:
break
return tours
Amazing results for the small instance with a total of 10 nodes. However, even for this small instance, the solver took almost half a minute to obtain the solution and the complexity increases significantly with more demand points. Fortunately, there are specialized algorithms publicly available to find good quality solutions for much larger instances in a short computational time. Let us take a look at them in the next section.
Specialized (Meta)Heuristics
Throughout the years several specialized (meta)heuristics have been proposed for variants of the VRP. Most of them strongly rely on local search algorithms so that, given a solution, different perturbations are tried in order to sequentially improve its cost until no further improvement is possible in the given neighborhood. When using Google Or Tools, we will too use local search methods associated with constructive algorithms.
In this section, instance 150d from Rochat and Taillard (1995) will be used. Its data was obtained from the CVRPLIB. This instance has 150 customers and one depot node, so we surely wouldn't be able to solve it using the MIP strategy previously presented.
Let us once again start by importing the libraries used.
from itertools import cycle
import numpy as np
import pandas as pd
from scipy.spatial.distance import pdist, squareform
import matplotlib.pyplot as plt
import matplotlib as mpl
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcpAnd let us load the problem data from a file that includes the coordinates and demand of each node.
dataset = pd.read_csv("./data/tai150d.csv", index_col=0)
coordinates = dataset.loc[:, ["x", "y"]]
demands = dataset.d.values
capacity = 1874
n_vehicles = 15
N = coordinates.shape[0]
distances = squareform(pdist(coordinates, metric="euclidean"))
distances = np.round(distances, decimals=4)The first step to using OR Tools VRP solver is to instantiate a routing manager and a model.
# Create the routing index manager: number of nodes, number of vehicles, depot node
manager = pywrapcp.RoutingIndexManager(
N, n_vehicles, 0
)
# Create Routing Model
routing = pywrapcp.RoutingModel(manager)Next, we will include callbacks to quantify dimensions related to arcs/edges and nodes. The method RegisterTransitCallback of our routing instance can be used to quantify any dimension related to arcs/edges whereas the method RegisterUnaryTransitCallback can quantify values related to nodes. It might be important to re-scale your parameters depending on their original magnitude because ortools considers only integer values in dimensions.
# Same valid for any callback related to arcs/edges
def distance_callback(from_index, to_index):
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return distances[from_node, to_node]
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
# Same valid for any callback related to nodes
def demand_callback(from_index):
from_node = manager.IndexToNode(from_index)
return demands[from_node]
demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)Now, we will include a capacity constraint using the demand_callback_index previously defined. Notice that duration constraints could have been defined using the same syntax just passing instances of RegisterTransitCallback as the first argument. Furthermore, it is important to remark that the routing model handles heterogeneous fleets, so we must pass a list of values in the third argument.
# Any constraint associated with vehicles can take same arguments
routing.AddDimensionWithVehicleCapacity(
demand_callback_index,
0, # null capacity slack
[capacity,] * n_vehicles, # vehicle maximum capacities (list for each vehicle)
True, # start cumul to zero
'Capacity'
)Similarly, the definition of the objective also takes a callback as the main argument. In this example, let us minimize the distances defined in the transit_callback_index.
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)At last, we must define solver parameters and solve our model.
# Setting heuristic strategies
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
routing_enums_pb2.FirstSolutionStrategy.CHRISTOFIDES
)
search_parameters.local_search_metaheuristic = (
routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH
)
search_parameters.time_limit.FromSeconds(300)
# Solve the problem
solution = routing.SolveWithParameters(search_parameters)The following piece of code can be used to extract the routes used in our solution.
tours = []
for vehicle_id in range(n_vehicles):
index = routing.Start(vehicle_id)
tours.append([])
while not routing.IsEnd(index):
node_index = manager.IndexToNode(index)
previous_index = index
index = solution.Value(routing.NextVar(index))
tours[-1].append(node_index)
else:
node_index = manager.IndexToNode(index)
tours[-1].append(node_index)And one can access the objective value by running the simple line solution.ObjectiveValue() . Using the configuration presented, I got an objective of 2679, which is quite close to the proven optimal value of 2645 (1.2% gap). The routes obtained are represented below.
The complete code (plots included) is available in this git repository.
Useful Extensions
The OR Tools library is fantastic as a general solver for routing problems for it handles several variants of the VRP, such as time windows, heterogeneous fleets, and multiple depots. However, algorithms fit for the canonical CVRP can perform even better. It is totally worth taking a look at the HGS-CVRP package (Vidal, 2022) that combines a state-of-the-art genetic algorithm with several local search moves. The algorithm finds the optimal solution for the instance tai150d within less than 20 seconds.
Regarding some real-world applications, it is likely that one should rely on road distances (or durations) rather than Euclidean distances to connect locations. Google provides a nice paid interface to do so, which you can check in this tutorial. However, if you are looking for open-source alternatives it is worth checking out the OpenStreetMap API. Some useful requests are:
- https://router.project-osrm.org/table/v1/driving/<LOCATIONS>
- http://router.project-osrm.org/route/v1/car/<LOCATIONS>?overview=false&steps=true
In which <LOCATIONS> should be a list of longitude and latitude pairs separated by commas within the pairs and by semicolons between different pairs. You can also specify sources and destinations in the table request, which is useful in case the complete table is too large to handle in a single request.
Besides doing precise routing calculations, visualization can be an important tool. The Python library folium can be quite useful to do it. It is definitely worth taking a look at it.
Further Reading
Earlier in this article we implemented an exact MIP model for the CVRP, which is not suitable for moderate-size instances. However, algorithms that combine column generation to Branch and Cut have been successful in solving instances with up to a few hundred customers. It is worth taking a look at the research papers of Fukasawa et al. (2006) and Pecin et al. (2017).
Those interested in a previous introduction to column generation might find it in my previous Medium article.
Regarding meta-heuristics, the papers of Vidal et al. (2012) and Vidal (2022) are fantastic. Both are also available in the form of technical reports, with links available in the HGS-CVRP repository.
Conclusions
In this article, two approaches for solving the Capacitated Vehicle Routing Problem (CVRP) were presented: Mixed-Integer Programming and (Meta)Heuristics. The first alternative was used to solve a small instance in which it has been successful, although it is not able to handle moderate-size or large instances. The second approach was used to solve a challenging problem from the literature with 150 customers to which the solver found a good quality solution with a 1.2% gap to the known optimal within 300s.
References
Bynum, M. L. et al., 2021. Pyomo-optimization modeling in python. Springer.
Dantzig, G. B., & Ramser, J. H., 1959. The truck dispatching problem. Management science, 6(1), 80–91.
Fukasawa, R., Longo, H., Lysgaard, J., Aragão, M. P. D., Reis, M., Uchoa, E., & Werneck, R. F., 2006. Robust branch-and-cut-and-price for the capacitated vehicle routing problem. Mathematical programming, 106, 491–511.
Pecin, D., Pessoa, A., Poggi, M., & Uchoa, E., 2017. Improved branch-cut-and-price for capacitated vehicle routing. Mathematical Programming Computation, 9, 61–100.
Rochat, Y., & Taillard, É. D., 1995. Probabilistic diversification and intensification in local search for vehicle routing. Journal of heuristics, 1, 147–167.
Toth, P., & Vigo, D., 2002. An overview of vehicle routing problems. The vehicle routing problem, 1–26.
Vidal, T., 2022. Hybrid genetic search for the CVRP: Open-source implementation and SWAP* neighborhood. Computers & Operations Research, 140, 105643.
Vidal, T., Crainic, T. G., Gendreau, M., Lahrichi, N., & Rei, W., 2012. A hybrid genetic algorithm for multidepot and periodic vehicle routing problems. Operations Research, 60(3), 611–624.
