A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer

A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer

Sebastian Feld
LMU Munich
Mobile and Distributed Systems Group
Munich, Germany
&Christoph Roch
LMU Munich
Mobile and Distributed Systems Group
Munich, Germany
&Thomas Gabor
LMU Munich
Mobile and Distributed Systems Group
Munich, Germany
&Christian Seidel
Volkswagen Data:Lab
Munich, Germany
&Florian Neukart
Volkswagen Group of America
San Francisco, USA
&Isabella Galter
Volkswagen Data:Lab
Munich, Germany
&Wolfgang Mauerer
OTH Regensburg
Regensburg, Germany
&Claudia Linnhoff-Popien
LMU Munich
Mobile and Distributed Systems Group
Munich, Germany

The Capacitated Vehicle Routing Problem (CVRP) is an NP-optimization problem (NPO) that has been of great interest for decades for both, science and industry. The CVRP is a variant of the vehicle routing problem characterized by capacity constrained vehicles. The aim is to plan tours for vehicles to supply a given number of customers as efficiently as possible. The problem is the combinatorial explosion of possible solutions, which increases superexponentially with the number of customers. Classical solutions provide good approximations to the globally optimal solution. D-Wave’s quantum annealer is a machine designed to solve optimization problems. This machine uses quantum effects to speed up computation time compared to classic computers. The problem on solving the CVRP on the quantum annealer is the particular formulation of the optimization problem. For this, it has to be mapped onto a quadratic unconstrained binary optimization (QUBO) problem. Complex optimization problems such as the CVRP can be translated to smaller subproblems and thus enable a sequential solution of the partitioned problem. This work presents a quantum-classic hybrid solution method for the CVRP. It clarifies whether the implemenation of such a method pays off in comparison to existing classical solution methods regarding computation time and solution quality. Several approaches to solving the CVRP are elaborated, the arising problems are discussed, and the results are evaluated in terms of solution quality and computation time.


A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer

  Sebastian Feldthanks: sebastian.feld@ifi.lmu.de LMU Munich Mobile and Distributed Systems Group Munich, Germany Christoph Roch LMU Munich Mobile and Distributed Systems Group Munich, Germany Thomas Gabor LMU Munich Mobile and Distributed Systems Group Munich, Germany Christian Seidel Volkswagen Data:Lab Munich, Germany Florian Neukart Volkswagen Group of America San Francisco, USA Isabella Galter Volkswagen Data:Lab Munich, Germany Wolfgang Mauerer OTH Regensburg Regensburg, Germany Claudia Linnhoff-Popien LMU Munich Mobile and Distributed Systems Group Munich, Germany


noticebox[b]Preprint. Work in progress.\end@float

1 Introduction

Optimization problems can be found in many different domains of applications, be it economics and finance (Black and Litterman, 1992), logistics (Caunhye et al., 2012), or healthcare (Cabrera et al., 2011). Their high complexity engaged reseachers to develop efficient methods for solving these problems (Papadimitriou and Steiglitz, 1998). With D-Wave Systems releasing the first commercially available quantum annealer in 2011 111https://www.dwavesys.com/news/d-wave-systems-sells-its-first-quantum-computing-system-lockheed-martin-corporation, there is now the possibility to find solutions for such problems in a completely different way than classical computation does. To use D-Wave’s quantum annealer the optimization problem has to be formulated as a quadratic unconstrained binary optimization (QUBO) problem (Boros et al., 2007), which is one of two input types acceptable by the machine (alternative: the Ising model (Glauber, 1963)). Doing this, the metaheuristic quantum annealing seeks for the minimum of the optimization function, i.e., the best solution of the defined configuration space (McGeoch, 2014). There has been recent research about solving real world problems on a quantum annealer, like Volkswagen’s Traffic Flow Optimization (Neukart et al., 2017) or the recently announced Tsunami Evacuation Optimization project by Tohoku University. 222https://www.dwavesys.com/sites/default/files/Ohzeki.pdf

The paper at hand regards the Capacitated Vehicle Routing Problem (CVRP), an NP-optimization problem that plays a major role in common operations research and is excessively studied since its proposal in (Dantzig and Ramser, 1959). The classic CVRP can be described as the problem of designing optimal routes from one depot to a number of geographically scattered customers subject to some side constraints (see Figure 1). It can be formulated as follows:

Let be a graph with being a set of vertices representing customer locations with the depot located at vertex and being a set of undirected edges. With every edge a non-negative cost is associated. This cost may, for instance, represent the (geographical) distance between two customers and . Furthermore, assume there are vehicles stationed at the depot that have the same capacity . In addition, every customer has a certain demand (Laporte, 1992). The CVRP consists of finding a set of vehicle routes such that

  • each customer in is visited exactly once by exactly one vehicle;

  • all routes start and end at the depot;

  • the sum of customer demand within a route does not exceed the vehicles’ capacity;

  • the sum of costs of all routes is minimal given the constraints above;

To solve the CVRP on D-Wave’s quantum annealer, the formulated QUBO problem has to be mapped to the hardware. However, quantum computation compared to classical computation is still in its infancy and one of the major problems is that quantum hardware is limited regarding the number of quantum bits (qubits) and their connectivity on the chip. Generally, this leads to difficulties in mapping large QUBO problems to the hardware.

With this paper we present an intuitive way to split the CVRP into smaller optimization problems by taking advantage of a classical 2-Phase-Heuristic (Laporte and Semet, 2002), see Figure 1. This heuristic divides the CVRP into two phases, the clustering phase and the routing phase. The clustering phase itself can be mapped to the NP-complete Knapsack Problem (KP) (Karp, 1972), which tries to pack different sized items (here: customers) into capacity restricted knapsacks (here: vehicles). Doing this, the sum of the objective values of the items in a knapsack should be maximized, i.e. the euclidean distance between customers assigned to a vehicle should be minimized. The routing phase can be represented by the NP-hard Travelling Salesman Problem (TSP) (Lawler, 1985). Thus, the minimal tour in which all customers of a cluster are visited once is sought. Doing this, the tour starts and ends in one place, i.e. the depot. Figure 1 shows a CVRP example with the 2-Phase-Heuristic. First the customers are grouped into clusters (b) before efficient vehicle routes in each cluster are searched (c).

(a) (a)
(b) (b)
(c) (c)
Figure 1: Overview of the CVRP and the 2-Phase-Heuristic. (a) Initial state with 9 customers and 1 depot. (b) Clustering phase results in three clusters found. (c) Routing phase determines shortest path inside each cluster.

In this paper we investigate different quantum-classic hybrid approaches to solve the CVRP, expound their difficulties in finding good solutions, and finally propose a hybrid method based on the 2-Phase-Heuristic to solve the CVRP using D-Wave’s quantum annealer. We map the optimization problems to a QUBO problem, and analyze performance from an application-specific perspective by using large benchmark datasets.

The paper is structured as follows: Section 2 gives a brief overview of existing methods for solving the CVRP. An introduction to quantum annealing and the common QUBO problem is given in Section 3. In Section 4, two approaches for solving the CVRP are briefly discussed before the concept of our hybrid method is presented. Section 5 first introduces the test setup, and subsequently presents and discusses the results with regard to solution quality and computational performance on commonly used CVRP datasets. Finally, we conclude this paper in Section 6.

2 Related Work

Over the last decades several families of heuristics have been proposed for solving the CVRP. They can be divided into construction heuristics, improvement heuristics and metaheuristics (Laporte and Semet, 2001).

Construction heuristics try to generate a good solution gradually. In every step, they insert customers into partial tours or combine sub-tours considering some capacities and costs to generate a complete solution. One of the most fundamental construction heuristics is the Clarke and Wright savings algorithm (Clarke and Wright, 1964), which first constructs a single tour for each customer, calculates the saving that can be obtained by merging those single customer tours, and iteratively combines the best sub-tours until no saving can be obtained anymore.

Improvement heuristics try to iteratively enhance a given feasible solution, which is often generated by a construction heuristic. A common methodology is to replace or swap customers between sub-tours taking capacity constraints into account. Popular improvement methods can be found in (Lin, 1965) and (Or, 1976).

Metaheuristics can be thought of as top-level strategies which guide local improvement operators to find a global solution. Groër et al. describe a library of local search heuristics for the (C)VRP (Groër et al., 2010). In addition, Crispin and Syrichas propose a classical quantum annealing metaheuristic for vehicle scheduling (Crispin and Syrichas, 2013). To approximate quantum annealing on a classical computer, they use a stochastic variant called Path-Integral Monte Carlo (PIMC) to simulate the quantum fluctuations of a quantum system. In our work, the quantum annealing hardware is responsible for that. However, the complexity exists in mapping the CVRP to a format readable by the hardware.

One of the most important classical 2-Phase-Heuristics is the Sweep algorithm (Gillett and Miller, 1974), where feasible clusters are formed by rotating a ray centered at the depot. After that the TSP is solved for each cluster. Fisher and Jaikumar also tried to solve the VRP with a cluster-first, route-second algorithm (Fisher and Jaikumar, 1981). They formulated a Generalized Assignment Problem (GAP) instead of using a geometry based method to form the clusters. Bramel and Simchi-Levi described a 2-Phase-Heuristic where the seeds were determined by solving capacitated location problems and the remaining vertices were gradually included into their allotted route in a second stage (Bramel and Simchi-Levi, 1995).

3 Quantum Annealing on D-Wave processor

Quantum annealing in general is a metaheuristic for solving complex optimization problems (Kadowaki and Nishimori, 1998). D-Wave’s quantum annealing algorithm is implemented in hardware using a framework of analog control devices to manipulate a collection of quantum bit (qubit) states according to a time-dependent Hamiltonian, denoted , shown in Equation 1.


The basic process of quantum annealing is to physically interpolate between an initial Hamiltonian with an easy to prepare minimal energy configuration (or ground state), and a problem Hamiltonian , whose minimal energy configuration is sought that corresponds to the best solution of the defined problem. This transition is described by an adiabatic evolution path which is mathematically represented as function and decreases from to (McGeoch, 2014). If this transition is executed sufficiently slow, the propability to find the ground state of the problem Hamiltonian is close to (Albash and Lidar, 2018).

D-Wave’s quantum annealing hardware is known to solve a specific optimization problem called a quadratic unconstrained binary optimization (QUBO) problem (Boros et al., 2007). QUBO is a unifying model which can be used for representing a wide range of combinatorial optimization problems. However, in order to use quantum annealing on D-Wave’s hardware the CVRP has to be formulated as a QUBO problem. The functional form of the QUBO the quantum annealer is designed to minimize is:


with being a vector of binary variables of size , and being an real-valued matrix describing the relationship between the variables. Given the matrix , the annealing process tries to find binary variable assignments to minimize the objective function in Equation 2.

The quantum processing unit (QPU) is a physical implementation of an undirected graph with qubits as vertices and couplers as edges between them. These qubits are arranged according in a so-called chimera graph, as illustrated in Figure 2. In relation to the QUBO problem, each qubit on the QPU represents such a QUBO variable and couplers between qubits represent the costs associated with qubit pairs, mathematically described in matrix . If the problem structure can not be embedded directly to the chimera graph, auxiliary qubits may be introduced to augment the available couplings.

Figure 2: Excerpt of the structure of a chimera graph. The full qubit graph extends to a lattice of groups of qubits. Figure in reference to (Biswas et al., 2017).

If – like in this paper – large data sets are used, the size of the resulting QUBO problem may exceed the limited number of available qubits on the QPU and the problem cannot be put on the chip altogether anymore. For this case, D-Wave provides a tool called QBSolv that splits the QUBO into smaller components and solves them sequentially on the D-Wave hardware 333https://github.com/dwavesystems/qbsolv. A detailed view on the QBSolv algorithm is given in Section 5.1.

4 Concept of Hybrid Solution Method

The classical 2-Phase-Heuristic separates the CVRP into a clustering phase and a routing phase. Both phases can be seen as detached optimization problems, the Knapsack Problem (KP) and the Travelling Salesman Problem (TSP), respectively. We have investigated three different approaches for the 2-Phase-Heuristic using a quantum annealer, see Figure 3. The insights of the preliminary exploration will be given in the following paragraphs. The concept of the most suitable approach will be presented in detail in Subsections 4.1 and 4.2.

Figure 3: Assignment of the clustering phase (KP) and the routing phase (TSP) of the 2-Phase-Heuristic to a classical solution method (Classic) or quantum annealing solution method (Quantum).

In the first approach (Q2Q in Figure 3), both optimization problems (KP and TSP) are mapped into two separate QUBO problems to be sequentially solved using the quantum annealer. However, we have faced severe difficulties within the clustering phase. The problem was to find values for an edge-weighting parameter such that customers having a short distance to each other (low cost) are grouped together. Experimental results have shown that the edge-weighting parameter needs to be chosen individually for each dataset. Thus, we perceive this approach as impractical.

The second approach (Q1Q in Figure 3) maps both optimization problems into one integrated QUBO problem, thus, the two phases are solved simultaneously. As a consequence, the QUBO problem includes two optimization functions: one minimizes the distances between customers to form dense clusters, and one tries to find the shortest route between customers inside a cluster. However, we observed that both optimization functions seemed to hinder each other. Neglecting the clustering optimization function led to valid routes inside the clusters, but at the same time the clusters were very sparse. The other way round, i.e. neglecting the routing optimization function, led to dense clusters but also to invalid routes inside the clusters. In summary, both efforts lead to invalid or unusable solutions.

The third approach (HS in Figure 3) as a candidate for a CVRP solution method combines the positive aspects of the previous mentioned approaches. To achieve this, the clustering phase (KP) is solved using a classical algorithm while the routing phase (TSP) is mapped to a QUBO problem in order to solve it on the quantum annealer.

4.1 Phase 1 – Clustering

The clustering phase of the hybrid solution method we propose is inspired by Shin and Han (Shin and Han, 2011). Based on their work, we add a charateristics called clustering core point parameter that will be presented below. The clustering phase can be subdivided: (1) cluster generation and (2) cluster improvement.

Within the cluster generation the core stop of a cluster, i.e. the first customer in a cluster, is chosen. (Koenig, 1995) propose to select the core stop either based on the maximum demand of the customers, or based on the largest distance to the depot. The motivation behind choosing the customer with the highest demand is the assumption that this one is the most critical customer in relation to the vehicles’ capacity constraint. After selecting that particular customer, the vehicle can be filled with goods for customers having a smaller demand. The motivation behind choosing the customer with the largest distance to the depot is the assumption that this one is the most critical customer in relation to the routes’ length constraint and that other customers may be supplied while approaching or receding that particular customer. Once the core stop of a cluster has been selected, the geometric center CC() of cluster is calculated using


with and being the and coordinates of customer and being the number of customers within cluster . Now, the customer with the smallest distance to the cluster center is selected from the set of unclustered customers and added to the cluster. After the cluster center is recalculated the steps are repeated until the demand of a customer to be added would exceed the vehicle’s capacity. If this is the case, a new core stop is selected based on the previously explained criteria and the still unclustered customers are assigned to the new cluster. This procedure stops when each customer has been assigned to a cluster.

Once the clusters have been generated, the cluster improvement is executed to enhance the clusters. This step assigns a customer belonging to cluster to another cluster , if that would reduce the distance to the cluster center, i.e. if the distance to is smaller than the distance to . However, assigning a customer to a new cluster must not violate the capacity limitation. If the reassignment is possible and valid, and are recalculated and the improvement process begins again. The improvement step terminates if it is not possible to assign a customer to another cluster or when a certain stop criterion is reached (e.g., number of iterations).

4.2 Phase 2 – Routing

After the clustering phase is completed, the goal is now to find the shortest route inside each cluster. Thus, the Travelling Salesman Problem (TSP) is executed for every generated cluster. The TSP can be reduced to the Hamiltonian Cycle Problem (HPC), which can be formulated as QUBO problem as follows (Lucas, 2014):


The binary variable is if the customer with index is located at position in the Hamiltonian Cycle. The first term (constraint) requires that each customer must occur only once in the cycle, while the second term enforces that each position in the cycle must be assigned to exactly one customer. This defines the order of the customers within the tour. The squared differences of these terms ensure that exactly one customer has a unique position in the tour. Otherwise, a high penalty value would be added to the solution energy, which states the solution itself as suboptimal or rather invalid. The third term ensures that the order of customers found is possible. That is, if and are both and with being the set of edges between the nodes representing the customers, then the penalty value should also state the solution invalid (Lucas, 2014). In this work we have evaluated our algorithm using CVRP/TSP datasets with fully meshed vertices, i.e. customers are connected with undirected edges. That is why the third term can be neglected.

In order to find the Hamiltonian Cycle with the shortest length, the following minimization function is needed:


Here is the euclidean distance between the customers and . The minimization function sums all costs of the edges between successive customers. The total solution for the TSP QUBO problem is then composed of:


The penalty value must be chosen sufficiently small so that it does not violate the constraint . A possible choice would be (Lucas, 2014). With , has to be chosen larger than the greatest distance between two customers.

By multiplying the QUBO formulas, one obtains the coefficients for the QUBO problem matrix which can be written as a lower (or upper) triangular matrix to be mapped to the quantum annealing processor. Figure 4 shows an excerpt of an exemplary QUBO problem in which only the coefficients are entered for simplification. As soon as several coefficients are noted on one cell of the matrix they must be added. In addition, every coefficient is multiplied with the penalty value and the distance is multiplied with penalty value .

Figure 4: Excerpt of a visualized TSP QUBO problem matrix: A1 corresponds to customer A on position 1 in the TSP cycle. The colored entries correspond to the coefficients and distances from Equation 4 and 5 , respectively.

5 Evaluation

In this section we present the results of the hybrid solving method with regards to solution quality and computational results. First the preliminaries for the test setup are given, then the test results are shown.

5.1 Preliminaries

As already mentioned, the D-Wave quantum annealer can be regarded as a discrete optimization machine that accepts problems in QUBO formulation. The QUBO problem matrix increases with the problem size, i.e. with the number of problem variables used. For the TSP logical variables, with being the number of customers, have to be used to describe it as a QUBO problem (see Section 4.2). These variables have to be mapped to the qubits and the logical links between them to the couplers of the physical QPU. Because of the almost fully meshed dependences between the logical variables it is not possible that the logical problem structure matches the physical one. For such issues D-Wave provides a minor embedding technique to find a valid embedding to the hardware. This technique is also included in D-Wave’s QBSolv tool, which we have used to fit our large QUBO problems to the physical hardware.

QBSolv splits the QUBO into smaller components (subQUBOs), which are then solved independently of each other. This process is executed iteratively as long as there is an improvement and it can be defined using the QBSolv parameter num_repeats. This parameter determines the number of times to repeat the splitting of the QUBO problem matrix after finding a better sample. With doing so, the QUBO matrix is split into different components using a classical tabu search heuristic in each iteration. QBSolv can be used in a completely classical way to solve the subQUBOs or as a quantum-classic hybrid method by solving the single subQUBOs on the quantum annealer. For more details about QBSolv, see (Michael Booth and Roy, 2017).

There exist many different benchmark datasets for the CVRP and the TSP, which can be downloaded from (Xavier, 2014), (Reinelt, 2013), (Cook, 2009). In addition, the Best Known Solution (BKS) of each dataset is noted. It gives information about the best solution, i.e., the shortest euclidean distance found by any solution method. In order to test and compare the proposed hybrid solution method with regard to the solution quality, various test datasets of Christofides and Eilon (Xavier, 2014) have been selected. Details about the CVRP datasets can be extracted from the name with the format E-nX-kY. For example, E-n22-k4 stands for a certain dataset , for the number of customers including the depot and for the minimal number of vehicles needed to solve the problem. The TSP datasets have the name format CityX, which just indicates the number of customers which have to be visited in the TSP tour. As already mentioned, the customer and depot coordinates relate to a 2D euclidean space.

5.2 Results

In this subsection the results of our hybrid method are presented. First, we exclusively analyze the TSP which is executed on the quantum annealer to see how the different-sized problem instances are handled.

5.2.1 TSP – Solution Quality

Table 1 shows the results for different-sized TSP datasets. The problem instance and its size, also included in the name, can be read from column one and two. In addition, the BKS, the best solution of 100 runs, and the average deviation of all 100 runs from the BKS, is noted in column three to five, respectively. One can see that for smaller sized problem instances (1)(2) the BKS has been found and the average deviation is generally low (0.00% to 0.31%). With increasing problem size (3)(4)(5) the BKS is not found and the average deviation increases continuously (2.70% to 25.91%). Therefore it can be concluded that the TSP can be solved comparatively well for smaller sized problem instances on the quantum annealer, while with regard to larger sized instances only a good approximation is found.

Problem Size
Best Known Solution
Best solution
found of 100 runs
Avg. deviation
of 100 runs from BKS
(1) Burma14 14 3323 3323 0.00%
(2) Ulysses16 16 6859 6859 0.31%
(3) Ulysses22 22 7013 7019 2.70%
(4) WesternSahara29 29 27603 28293 8.16%
(5) Djibouti38 38 6656 7396 25.91%
Table 1: Results for various TSP datasets. Parameter num_repeats was set to 250.

In Figure 5 the deviation of each found solution from the BKS is visualized with a boxplot diagram. In this figure the variance of each test result in relation to the deviation is shown in more detail. It also seems that with increasing problem size the variance of the results expands. Whereas for datasets (1)(2)(3) the variance is comparatively low (0.00%, 0.00%-1.56%, and 0.09%-5.29%), larger datasets (4)(5) show more fluctuations (2.50%-13.77% and 11.12%-36.01%).

Figure 5: Deviation of the results from the BKS for each dataset. Every dataset was run 100 times and num_repeats was set to 250.

Experimental tests showed that the solution quality depends on the QBSolv parameter num_repeats. In Figure 6 the deviation from the BKS for the already known datasets with different settings for the num_repeats parameter are shown. Dataset Burma14 has been neglected since even num_repeats set to 50 finds the BKS in every run. Each setting has been executed 10 times. One can see that with increasing the num_repeats parameter there is a tendency that the solution quality improves, i.e., the deviation from the BKS decreases.

(a) (a) Ulysses16
(b) (b) Ulysses22
(c) (c) WesternSahara29
(d) (d) Djibouti38
Figure 6: Different settings for the num_repeats parameter of QBSolv for various datasets.

5.2.2 CVRP – Solution Quality

Now the results of the hybrid method including its classical clustering phase are presented. One can optimize the clustering of the customers by choosing the core stop of a cluster (max_distance or max_request). Depending on the selected parameter, either the customer with the largest distance to the depot or the customer with the highest demand is set as the seed of a cluster. In Table 2 the best solution found is noted, i.e. the sum of all vehicle routes and the deviation from the BKS for both options of the core stop parameter.

Hybrid Solution Method
Problem BKS
Shortest distance
Shortest distance
Avg. deviation
from BKS
Avg. deviation
from BKS
(1) E-n22-k4 375 385 407 2.66% 8.53%
(2) E-n33-k4 835 965 852 15.07% 2.05%
(3) E-n51-k5 521 618 557 18.75% 6.91%
(4) E-n76-k7 682 748 814 9.70% 19.46%
(5) E-n101-k8 815 892 898 11.16% 11.39%
Table 2: Distances and deviations from the BKS for certain E datasets of Christofides and Eilon. Every dataset was executed 10 times.

The results do not allow a concrete statement about the choice of the core stop parameter of the hybrid method. One can see, that it is independent of the problem size. In fact the core point has to be selected individually for each dataset to obtain a good clustering and a short route length as a consequence. Especially for the first three datasets (1)(2)(3) the hybrid method finds good approximations to the BKS with regard to the solution quality (2.66%-6.91% deviation).

In Table 3 the solution quality of the hybrid method is compared to other fundamental construction and 2-phase-heuristics. The solution quality has been compared based on seven CMT datasets of Christofides, Mingozzi and Toth (Xavier, 2014). Datasets CMT6-CMT10 have been neglected in the evaluation as they are identical to CMT1-CMT5, however, they regard an additional time window, which in turn is ignored in the classic CVRP. The solutions for the selected datasets found by the other heuristics were extracted from (Fisher and Jaikumar, 1981). For every solution method the best solution found (Distance) and the deviation (Dev.) from the BKS is depicted. In addition, the problem size and the BKS for each dataset is noted in column two and three.

Christofides et al.
2 Phase
Problem Size BKS Distance Dev. Distance Dev. Distance Dev. Distance Dev. Distance Dev.
CMT1 50 524.61 585 11.5% 524 0.12% 550 4.84% 532 1.41% 556b 5.98%
CMT2 75 835.26 900 7.75% 857 2.60% 883 5.72% 874 4.64% 926a 10.86%
CMT3 100 826.14 886 7.25% 833 0.83% 851 3.01% 851 3.01% 905a 9.55%
CMT4 150 1028.42 1204 17.07% - - 1093 6.28% 1079 4.92% 1148a 11.63%
CMT5 199 1291.29 1540 19.26% 1420 9.97% 1418 9.81% 1389 7.57% 1429a 10.66%
CMT11 120 1042.12 - - - - - - - - 1084b 4.02%
CMT12 100 819.56 877 7.01% 848 3.47% 876 6.89% 949 15.79% 828a 1.03%
Table 3: Results for selected CMT datasets: Solution found with (a) max_distance or (b) max_request

For problem instances CMT1 and CMT5, the hybrid method was able to keep up or even surpass the Saving-Heuristic of Clarke and Wright. With regard to datasets CMT2, CMT3 and CMT4 major deviations from the BKS are recorded. Here the hybrid method can keep up with or even dominate the Saving-Heuristic. However, with regard to the other heuristics it was not competitive. The last two problem instances are structured problems in which the customers are organized in clusters by the authors of the datasets. According to Fisher and Jaikumar, these datasets resemble real problems rather than problem instances CMT1-CMT5 (Fisher and Jaikumar, 1981). Therefore, with these instances one can recognize that the clustering algorithm of the hybrid method works comparatively well, since each of the other four heuristics was surpassed with regard to the solution quality. However, in relation to the BKS it should be mentioned that both, the hybrid method and the other heuristics, never found the BKS (exception: problem instance CMT1 - Fisher and Jaikumar). This is basically a well-known problem of clustering customers in the CVRP, since finding the BKS or the global optimum is related to the geographical structure of the customers of a problem instance.

5.2.3 CVRP – Computational Results

The computation time for the executed test instances must be considered differentiated. As already mentioned, QBSolv can be used as a pure classical solver as well as a hybrid solver for large QUBO problems. With the classic version, the subQUBOs are solved locally, while with the hybrid version the subQUBOs are solved sequentially on the D-Wave hardware. Doing this, the QUBO is split locally and the subQUBOs are sent to the hardware via a remote connection and the individual jobs are placed in a queue. As a result of this process, additional latency and possibly waiting times may occur.

To demonstrate the difference in relation to the computational time, the CMT1 dataset has been solved locally as well as on the quantum annealer using QBSolv. The dataset has been run on a Dell 2.8 GHz i7 with 16 GB RAM Notebook. For measuring the CPU processing time the Python module cProfile has been used. In Listing 1 and Listing 2 the measured CPU times for the locally executed QBSolv are shown, while Listing 3 and Listing 4 show the same for the remotely executed QBSolv. While the classical run took seconds to complete (see first line in Listing 2), the hybrid version needed seconds (first line in Listing 4). A distance of has been found in the classical way (last line in Listing 1) and the hybrid version found a distance of (last line in Listing 3).

1number of clusters: 5
4-4922.00000 Energy of solution
54 Number of Partitioned calls, 1 output sample
60.01600 seconds of classic cpu time
9-7646.00000 Energy of solution
104 Number of Partitioned calls, 1 output sample
110.04600 seconds of classic cpu time
14-12552.00000 Energy of solution
152 Number of Partitioned calls, 1 output sample
160.01600 seconds of classic cpu time
19-6903.00000 Energy of solution
208 Number of Partitioned calls, 1 output sample
210.03100 seconds of classic cpu time
24-16253.00000 Energy of solution
2532 Number of Partitioned calls, 1 output sample
260.12500 seconds of classic cpu time
28Total distance 557
Listing 1: Output of locally used QBSolv.
1738180 function calls (661690 primitive calls) in 1.474 seconds
3Ordered by: internal time
5ncalls  tottime  percall  cumtime  percall      filename:lineno(function)
65       0.262    0.052    0.262    0.052        {dwave_qbsolv.qbsolv_binding.run_qbsolv}
751      0.066    0.001    1.258    0.025        __init__.py:1(<module>)
8364     0.058    0.000    0.058    0.000        {imp.find_module}
937752   0.056    0.000    0.056    0.000        binary_quadratic_model.py:484(add_interaction)
105       0.043    0.009    0.044    0.009        TSPPreparer.py:56(generate_adjacency_matrix)
1178703   0.040    0.000    0.060    0.000        records.py:438(__getattribute__)
125       0.033    0.007    0.033    0.007        TSPPreparer.py:25(generate_edge_list_tsp)
135       0.026    0.005    0.028    0.006        TSPPreparer.py:98(generate_edge_matrix)
14705     0.025    0.000    0.025    0.000        {compile}
155       0.023    0.005    0.078    0.016        binary_quadratic_model.py:595(add_interactions_from)
164       0.020    0.005    0.202    0.050        __init__.py:4(<module>)
1711195   0.020    0.000    0.093    0.000        response.py:675(__getitem__)
185       0.016    0.003    0.018    0.004        TSPSolver.py:11(generate_tsp_qubo)
Listing 2: Output of cProfile for locally used QBSolv.
1number of clusters: 5
4-4922.00000 Energy of solution
54 Number of Partitioned calls, 1 output sample
63.31200 seconds of classic cpu time
9-7646.00000 Energy of solution
102 Number of Partitioned calls, 1 output sample
112.03200 seconds of classic cpu time
14-12551.00000 Energy of solution
152 Number of Partitioned calls, 1 output sample
161.25100 seconds of classic cpu time
19-6903.00000 Energy of solution
204 Number of Partitioned calls, 1 output sample
213.56600 seconds of classic cpu time
24-16246.00000 Energy of solution
252 Number of Partitioned calls, 1 output sample
261.20000 seconds of classic cpu time
28Total distance 566
Listing 3: Output of remotely used QBSolv.
1882498 function calls (821251 primitive calls) in 15.792 seconds
3Ordered by: internal time
5ncalls  tottime  percall  cumtime  percall      filename:lineno(function)
6241     11.205   0.046    11.205   0.046        {method ’acquire’ of ’thread.lock’ objects}
75023    1.028    0.000    1.028    0.000        {built-in method __new__ of type object at 0x53A271E0}
855      0.926    0.017    0.926    0.017        {method ’read’ of ’_ssl._SSLSocket’ objects}
91       0.618    0.618    0.618    0.618        {minorminer.find_embedding}
101       0.368    0.368    0.368    0.368        {method ’do_handshake’ of ’_ssl._SSLSocket’ objects}
111       0.181    0.181    0.181    0.181        {method ’connect’ of ’_socket.socket’ objects}
125       0.077    0.015    11.369   2.274        {dwave_qbsolv.qbsolv_binding.run_qbsolv}
1351      0.068    0.001    1.287    0.025        __init__.py:1(<module>)
1441476   0.065    0.000    0.065    0.000        binary_quadratic_model.py:484(add_interaction)
155       0.044    0.009    0.045    0.009        TSPPreparer.py:56(generate_adjacency_matrix)
165       0.035    0.007    0.035    0.007        TSPPreparer.py:25(generate_edge_list_tsp)
1761672   0.033    0.000    0.049    0.000        records.py:438(__getattribute__)
183       0.031    0.010    0.032    0.011        decoder.py:370(raw_decode)
19701     0.026    0.000    0.091    0.000        binary_quadratic_model.py:595(add_interactions_from)
205       0.026    0.005    0.028    0.006        TSPPreparer.py:98(generate_edge_matrix)
21705     0.025    0.000    0.025    0.000        {compile}
225       0.017    0.003    0.019    0.004        TSPSolver.py:11(generate_tsp_qubo)
Listing 4: Output of cProfile for remotely used QBSolv.

The enormous difference in computational time can be seen in the method listing of cProfile in Listings 2 and 4, respectively. Here it can clearly be seen that the hybrid version of QBSolv mainly stays in the method method ’acquire’ of ’thread.lock’ objects, which is not listed in the locally QBSolv run. Therefore we assume that here the main QBSolv thread waits for the child threads to find a valid embedding for the respective SubQUBOs to the D-Wave hardware. This process is not needed using QBSolv locally. However, the actual annealing time for solving a QUBO problem remotely is in the range of microseconds.

The hybrid QBSolv method splits the problem QUBOs into , , , and subQUBOs (Number of Partitioned Calls, see Listing 3) and solved them one after the other on the quantum annealer. In addition, Listing 5 contains the individual qpu_access_time per subQUBO in microseconds. These add up to seconds. The classic QBSolv version requires seconds to solve the QUBOs.

1qpu_access_time 7783            qpu_access_time 7789
2qpu_access_time 7787            qpu_access_time 7783
3qpu_access_time 7808            qpu_access_time 7807
4qpu_access_time 7814            qpu_access_time 7778
5qpu_access_time 7768            qpu_access_time 7795
6qpu_access_time 7802            qpu_access_time 7788
7qpu_access_time 7807            qpu_access_time 7782
Listing 5: qpu_access_time of using QBSolv remotely. The qpu_access_time includes the whole time range from programming to post processing overhead time.

In summary, it can be stated that the real computation time to solve a QUBO problem on D-Wave’s quantum annealer is comparatively shorter than with the classic QBSolv version ( speedup). However, methods like finding a valid embedding of the QUBO problem to the hardware, which is not needed for using QBSolv locally, generates overhead. For this reason, the classic version of QBSolv is currently more advantageous in practice regarding the total calculation time.

6 Conclusion

To the best of our knowledge, this work presents the first study that solves the capacitated vehicle routing problem (CVRP) using quantum annealing hardware. The most important step was to find an intuitive way to map this optimization problem to a QUBO problem. Doing this, the classical 2-phase-heuristic has been used, which enables to divide the complex problem into a clustering phase as well as a routing phase and solves them sequentially or simultaneously (see approaches Q2Q and Q1Q in Figure 3). Due to various problems within the clustering phase, a hybrid method proved to be the best. We showed that our hybrid method was able to compete with other classical construction and 2-phase-heuristics and in some cases even surpass them with regard to solution quality. However, it should be mentioned that there are other solution methods like metaheuristics, which provide a better solution with regard to the used benchmark datasets. Especially when using datasets whose BKS contains overlapping routes, the hybrid method – and in general heuristics with clustering methods that work distance-based – has got difficulties in finding the BKS. However, the hybrid method usually provided a good approximation.

The computational time of the hybrid solution method must be considered differentiated. Due to the currently limited number of available qubits on the D-Wave QPU, the tool QBSolv must be used for large QUBO problems which are not directly embeddable on the D-Wave chip. This tool makes it possible to split the QUBO into smaller subQUBOs and place them one after the other on the quantum annealer. However, this hybrid solution option involves certain latency and waiting times which lacks the hoped acceleration of the computational time compared to the classical option. Then again, it should be mentioned that the real solution time for an embeddable QUBO problem on the D-Wave quantum annealer is located in the range of microseconds.

At the 2018 D-Wave Qubits Europe users conference D-Wave provided an outlook about the future hardware directions of quantum annealing. They stated that the connectivity and the number of qubits on D-Wave machines will immensely rise over the next years 444https://www.dwavesys.com/sites/default/files/mwj_dwave_qubits2018.pdf. These news give hope that in the future D-Wave’s quantum annealers are more suitable for these kind of optimization problems and a shorter total computation time can be achieved.


Research was funded by Volkswagen Group, department Group IT.


  • [1]
  • Albash and Lidar [2018] Albash, T. and Lidar, D. A. (2018). Adiabatic quantum computation. Rev. Mod. Phys. 90, 015002. doi:10.1103/RevModPhys.90.015002
  • Biswas et al. [2017] Biswas, R., Jiang, Z., Kechezhi, K., Knysh, S., Mandrà, S., O’Gorman, B., et al. (2017). A nasa perspective on quantum computing: Opportunities and challenges. Parallel Computing 64, 81 – 98. doi:https://doi.org/10.1016/j.parco.2016.11.002. High-End Computing for Next-Generation Scientific Discovery
  • Black and Litterman [1992] Black, F. and Litterman, R. (1992). Global portfolio optimization. Financial analysts journal , 28–43
  • Boros et al. [2007] Boros, E., Hammer, P. L., and Tavares, G. (2007). Local search heuristics for quadratic unconstrained binary optimization (qubo). Journal of Heuristics 13, 99–132
  • Bramel and Simchi-Levi [1995] Bramel, J. and Simchi-Levi, D. (1995). A location based heuristic for general routing problems. Operations Research 43, 649–660
  • Cabrera et al. [2011] Cabrera, E., Taboada, M., Iglesias, M. L., Epelde, F., and Luque, E. (2011). Optimization of healthcare emergency departments by agent-based simulation. Procedia computer science 4, 1880–1889
  • Caunhye et al. [2012] Caunhye, A. M., Nie, X., and Pokharel, S. (2012). Optimization models in emergency logistics: A literature review. Socio-economic planning sciences 46, 4–13
  • Clarke and Wright [1964] Clarke, G. and Wright, J. W. (1964). Scheduling of vehicles from a central depot to a number of delivery points. Operations Research 12, 568–581
  • Cook [2009] [Dataset] Cook, W. (2009). The traveling salesman problem. http://www.math.uwaterloo.ca/tsp/world/countries.html
  • Crispin and Syrichas [2013] Crispin, A. and Syrichas, A. (2013). Quantum annealing algorithm for vehicle scheduling. In Systems, Man, and Cybernetics (SMC), 2013 IEEE International Conference on (IEEE), 3523–3528
  • Dantzig and Ramser [1959] Dantzig, G. B. and Ramser, J. H. (1959). The truck dispatching problem. Management science 6, 80–91
  • Fisher and Jaikumar [1981] Fisher, M. L. and Jaikumar, R. (1981). A generalized assignment heuristic for vehicle routing. Networks 11, 109–124. doi:10.1002/net.3230110205
  • Gillett and Miller [1974] Gillett, B. E. and Miller, L. R. (1974). A heuristic algorithm for the vehicle-dispatch problem. Operations Research 22, 340–349. doi:10.1287/opre.22.2.340
  • Glauber [1963] Glauber, R. J. (1963). Time-dependent statistics of the ising model. Journal of mathematical physics 4, 294–307
  • Groër et al. [2010] Groër, C., Golden, B., and Wasil, E. (2010). A library of local search heuristics for the vehicle routing problem. Mathematical Programming Computation 2, 79–101. doi:10.1007/s12532-010-0013-5
  • Kadowaki and Nishimori [1998] Kadowaki, T. and Nishimori, H. (1998). Quantum annealing in the transverse ising model. Physical Review E 58, 5355
  • Karp [1972] Karp, R. M. (1972). Reducibility among combinatorial problems. In Complexity of computer computations (Springer). 85–103
  • Koenig [1995] Koenig, B. (1995). Heuristiken zur ein-depot-tourenplanung
  • Laporte [1992] Laporte, G. (1992). The vehicle routing problem: An overview of exact and approximate algorithms. European journal of operational research 59, 345–358
  • Laporte and Semet [2001] Laporte, G. and Semet, F. (2001). The vehicle routing problem (Philadelphia, PA, USA: Society for Industrial and Applied Mathematics), chap. Classical Heuristics for the Capacitated VRP. 109–128
  • Laporte and Semet [2002] Laporte, G. and Semet, F. (2002). Classical Heuristics for the Capacitated VRP, chap. 5. 109–128. doi:10.1137/1.9780898718515.ch5
  • Lawler [1985] Lawler, E. L. (1985). The traveling salesman problem: a guided tour of combinatorial optimization. Wiley-Interscience Series in Discrete Mathematics
  • Lin [1965] Lin, S. (1965). Computer solution of the traveling salesman problem 10, 2245–
  • Lucas [2014] Lucas, A. (2014). Ising formulations of many np problems. Frontiers in Physics 2, 5
  • McGeoch [2014] McGeoch, C. C. (2014). Adiabatic quantum computation and quantum annealing: Theory and practice. Synthesis Lectures on Quantum Computing 5, 1–93
  • Michael Booth and Roy [2017] Michael Booth, S. P. R. and Roy, A. (2017). Partitioning optimization problems for hybrid classical/quantum execution
  • Neukart et al. [2017] Neukart, F., Compostella, G., Seidel, C., von Dollen, D., Yarkoni, S., and Parney, B. (2017). Traffic flow optimization using a quantum annealer. Frontiers in ICT 4, 29
  • Or [1976] Or, I. (1976). Traveling salesman-type combinatorial problems and their relation to the logistics of regional blood banking UMI order no. 77-10076
  • Papadimitriou and Steiglitz [1998] Papadimitriou, C. H. and Steiglitz, K. (1998). Combinatorial optimization: algorithms and complexity (Courier Corporation)
  • Reinelt [2013] [Dataset] Reinelt, G. (2013). TSPLIB. https://wwwproxy.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/
  • Shin and Han [2011] Shin, K. and Han, S. (2011). A centroid-based heuristic algorithm for the capacitated vehicle routing problem 30, 721–732
  • Xavier [2014] [Dataset] Xavier, I. (2014). CVRP LIB. http://vrp.atd-lab.inf.puc-rio.br/index.php/en/
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description