A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer
Abstract
The Capacitated Vehicle Routing Problem (CVRP) is an NPoptimization 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. DWave’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 quantumclassic 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 Feld^{†}^{†}thanks: 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 LinnhoffPopien 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 DWave Systems releasing the first commercially available quantum annealer in 2011 ^{1}^{1}1https://www.dwavesys.com/news/dwavesystemssellsitsfirstquantumcomputingsystemlockheedmartincorporation, there is now the possibility to find solutions for such problems in a completely different way than classical computation does. To use DWave’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. ^{2}^{2}2https://www.dwavesys.com/sites/default/files/Ohzeki.pdf
The paper at hand regards the Capacitated Vehicle Routing Problem (CVRP), an NPoptimization 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 nonnegative 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 DWave’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 2PhaseHeuristic (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 NPcomplete 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 NPhard 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 2PhaseHeuristic. First the customers are grouped into clusters (b) before efficient vehicle routes in each cluster are searched (c).
In this paper we investigate different quantumclassic hybrid approaches to solve the CVRP, expound their difficulties in finding good solutions, and finally propose a hybrid method based on the 2PhaseHeuristic to solve the CVRP using DWave’s quantum annealer. We map the optimization problems to a QUBO problem, and analyze performance from an applicationspecific 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 subtours 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 subtours 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 subtours taking capacity constraints into account. Popular improvement methods can be found in (Lin, 1965) and (Or, 1976).
Metaheuristics can be thought of as toplevel 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 PathIntegral 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 2PhaseHeuristics 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 clusterfirst, routesecond 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 SimchiLevi described a 2PhaseHeuristic 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 SimchiLevi, 1995).
3 Quantum Annealing on DWave processor
Quantum annealing in general is a metaheuristic for solving complex optimization problems (Kadowaki and Nishimori, 1998). DWave’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 timedependent Hamiltonian, denoted , shown in Equation 1.
(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).
DWave’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 DWave’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:
(2) 
with being a vector of binary variables of size , and being an realvalued 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 socalled 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.
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, DWave provides a tool called QBSolv that splits the QUBO into smaller components and solves them sequentially on the DWave hardware ^{3}^{3}3https://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 2PhaseHeuristic 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 2PhaseHeuristic 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.
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 edgeweighting parameter such that customers having a short distance to each other (low cost) are grouped together. Experimental results have shown that the edgeweighting 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
(3) 
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):
(4) 
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:
(5) 
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:
(6) 
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 .
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 DWave 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 DWave provides a minor embedding technique to find a valid embedding to the hardware. This technique is also included in DWave’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 quantumclassic 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 EnXkY. For example, En22k4 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 differentsized problem instances are handled.
5.2.1 TSP – Solution Quality
Table 1 shows the results for differentsized 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 




(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% 
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%).
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.
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.


Problem  BKS 





(1) En22k4  375  385  407  2.66%  8.53%  
(2) En33k4  835  965  852  15.07%  2.05%  
(3) En51k5  521  618  557  18.75%  6.91%  
(4) En76k7  682  748  814  9.70%  19.46%  
(5) En101k8  815  892  898  11.16%  11.39% 
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 2phaseheuristics. The solution quality has been compared based on seven CMT datasets of Christofides, Mingozzi and Toth (Xavier, 2014). Datasets CMT6CMT10 have been neglected in the evaluation as they are identical to CMT1CMT5, 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.



Sweep 


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%  556^{b}  5.98%  
CMT2  75  835.26  900  7.75%  857  2.60%  883  5.72%  874  4.64%  926^{a}  10.86%  
CMT3  100  826.14  886  7.25%  833  0.83%  851  3.01%  851  3.01%  905^{a}  9.55%  
CMT4  150  1028.42  1204  17.07%      1093  6.28%  1079  4.92%  1148^{a}  11.63%  
CMT5  199  1291.29  1540  19.26%  1420  9.97%  1418  9.81%  1389  7.57%  1429^{a}  10.66%  
CMT11  120  1042.12                  1084^{b}  4.02%  
CMT12  100  819.56  877  7.01%  848  3.47%  876  6.89%  949  15.79%  828^{a}  1.03% 
For problem instances CMT1 and CMT5, the hybrid method was able to keep up or even surpass the SavingHeuristic 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 SavingHeuristic. 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 CMT1CMT5 (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 wellknown 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 DWave 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).
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 DWave 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.
In summary, it can be stated that the real computation time to solve a QUBO problem on DWave’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 2phaseheuristic 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 2phaseheuristics 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 distancebased – 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 DWave QPU, the tool QBSolv must be used for large QUBO problems which are not directly embeddable on the DWave 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 DWave quantum annealer is located in the range of microseconds.
At the 2018 DWave Qubits Europe users conference DWave provided an outlook about the future hardware directions of quantum annealing. They stated that the connectivity and the number of qubits on DWave machines will immensely rise over the next years ^{4}^{4}4https://www.dwavesys.com/sites/default/files/mwj_dwave_qubits2018.pdf. These news give hope that in the future DWave’s quantum annealers are more suitable for these kind of optimization problems and a shorter total computation time can be achieved.
Acknowledgement
Research was funded by Volkswagen Group, department Group IT.
References
 [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. HighEnd Computing for NextGeneration 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 SimchiLevi [1995] Bramel, J. and SimchiLevi, 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 agentbased 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. Socioeconomic 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 vehicledispatch problem. Operations Research 22, 340–349. doi:10.1287/opre.22.2.340
 Glauber [1963] Glauber, R. J. (1963). Timedependent 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/s1253201000135
 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 eindepottourenplanung
 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. WileyInterscience 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 salesmantype combinatorial problems and their relation to the logistics of regional blood banking UMI order no. 7710076
 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.uniheidelberg.de/groups/comopt/software/TSPLIB95/
 Shin and Han [2011] Shin, K. and Han, S. (2011). A centroidbased heuristic algorithm for the capacitated vehicle routing problem 30, 721–732
 Xavier [2014] [Dataset] Xavier, I. (2014). CVRP LIB. http://vrp.atdlab.inf.pucrio.br/index.php/en/