Quantum Walk-Based Vehicle Routing Optimisation

I Introduction

Quantum computers may offer a unique advantage when it comes to many combinatorial optimisation problems, where the search for an optimal solution becomes quickly infeasible due to solution spaces that scale exponentially with increasing problem size [1]. While a classical computer’s central processing unit is restricted to assessing the quality of solutions one after another, a quantum processor is capable of operating on the complete solution space at once via the principle of quantum superposition. However, this capability alone does not provide any significant utility. For example, consider a quantum system initialised in an equal superposition of states, one for each solution in the solution space. A measurement of the system is equally likely to collapse the system into any one of these states, which is equivalent to simply picking a solution at random. Where quantum superposition offers significant utility, is when it is combined with a suitable quantum algorithm. Such an algorithm should be capable of significantly increasing the probability of measuring the system in a state corresponding to the most, or one of the most optimal solutions.

One candidate is the Quantum Walk-based Optimisation Algorithm (QWOA) [2, 3]. This algorithm makes use of alternating quality-dependent phase-shifts and continuous-time quantum walks over a circulant graph connecting all possible solutions in the solutions space. By tuning the applied phase-shifts and quantum walks via a set of classical parameters, it is hoped that quantum interference will result in a concentration of probability density at states corresponding to high-quality solutions. This tuning process is carried out via a classical optimisation procedure which optimises for the expectation value of quality, as measured from the quantum circuit. Its effectiveness has been recently demonstrated in the context of portfolio optimisation problems [4].

In this paper, we show the applicability of the QWOA to the combinatorial optimisation problem of capacitated vehicle routing. The capacitated vehicle routing problem (CRVP) asks which route(s) should be taken by supply vehicles with a limited capacity in order to deliver products from a single depot to multiple locations, each requiring unique quantities of various products. Trips to and from the depot as well as between external locations are all characterised by route-dependent costs. The globally optimal solution is the route or set of routes that successfully delivers all required products whilst minimising the total delivery cost, subject to the further constraint that every vehicle must return to the depot upon completion of delivery. A general form of the CVRP was first studied by Dantzig and Ramser [5] in 1959, and although algorithms have since been developed to solve smaller-scale problems exactly [6], the focus for larger scale problems has been on heuristic-based methods for finding near-optimal solutions. The CVRP has clear real world significance because improving vehicle schedules, even by a tiny proportion, can lead to a large reduction in the transportation costs over time. Such problems have very recently been the subject of quantum algorithm development, in particular, using the quantum annealing approach [7–10]. This paper, however, focuses on the application of QWOA to the Capacitated Vehicle Routing Problem (CVRP), an algorithm that operates within the gate-based model of quantum computation. The purpose of this paper is to demonstrate that the QWOA can be effectively applied to the CRVP, and to present and analyse the results of numerical simulations of the application of the algorithm to an example problem.

The structure of this paper is as follows: In Section II, the CVRP will be formally introduced, along with an illustrative example and a brief discussion of its solution space. In Section III, the QWOA will be introduced, including its theoretical framework and quantum circuit implementation. In Section IV, the CVRP will be shown to satisfy the necessary prerequisite features for application of the QWOA, including possessing efficient processes for computing the cardinality of the solution space, for indexing/unindexing of the solution space, and for computing solution qualities. Finally, in Section V, the numerical results of the simulated application of the QWOA to an example CVRP will be presented and analysed.

II The Capacitated Vehicle Routing Problem

A. Formal Definition

A capacitated vehicle routing problem consisting of n delivery locations shall be referred to as a problem of size n. The delivery network for such a problem shall be characterised by a complete graph with n + 1 vertices. The vertex representing the depot is labelled with a zero, and the delivery locations are represented by vertices labelled 1, 2, …, n. The number of packages required at each location are included in a package vector, P, of dimension n, containing non-negative integers, P_i, where P_i is the number of packages required at location i and i = 1, 2, …, n. The costs associated with the trips between nodes are captured in a cost matrix, C. The cost matrix is square, n + 1 dimensional, and element C_ij is the cost of the trip from node i to node j with i, j = 0, 1, …, n and with C_ij taking positive finite values. Since it makes no sense to talk of a trip from a node to itself, the leading diagonal of the cost matrix consists of zeros. If costs are equal in both directions for all trips between nodes, then the cost matrix will be symmetrical, however, this need not be the case. The vehicle capacity is represented by the natural number variable, V. No consideration is made for time, and any solution which minimises cost can be undertaken by a single vehicle. Thus, while particular solutions (delivery routes) may be undertaken with multiple vehicles, the number of vehicles does not, and need not, make an appearance in this formulation of the CVRP.

Given a particular instance of the CVRP, solving the problem reduces to finding a solution from the space of all possible solutions, $\mathbb{S}$, which minimises total cost. Or to put it more formally, for some instance of the CVRP, specified by a unique set of the parameters, P, C and V, and with a quality function, f(x), which returns the total cost of any given solution x, the problem reduces to finding an optimal solution x*, such that

$f\left(x^{*}\right)=\min \left\{f\left(x\right)\mid x\in \mathbb{S}\right\}.(1)$

The difficulty lies in the total number of possible solutions, M, which grows exponentially with increasing problem size n.

B. Illustrative Example

A CVRP is fully defined by the three aforementioned variables, P, C, and V. As an example, the following variables

$P=\left[\begin{array}{c}\hfill 14\hfill \\ \hfill 24\hfill \\ \hfill 8\hfill \end{array}\right],C=\left[\begin{array}{cccc}\hfill 0\hfill & \hfill 16\hfill & \hfill 19\hfill & \hfill 12\hfill \\ \hfill 16\hfill & \hfill 0\hfill & \hfill 12\hfill & \hfill 17\hfill \\ \hfill 19\hfill & \hfill 12\hfill & \hfill 0\hfill & \hfill 10\hfill \\ \hfill 12\hfill & \hfill 17\hfill & \hfill 10\hfill & \hfill 0\hfill \end{array}\right],V=20$

fully specify a CVRP of size n = 3, as illustrated in Figure 1.

FIGURE 1

FIGURE 1. Shown on the left is an example of a CVRP of size n = 3 and with vehicle capacity V = 20. In blue are the required package numbers P_i, and each trip/edge has its cost C_ij shown in red. On the right is the route for this particular problem corresponding to the solution (1, 2), (3).

The example solution shown in Figure 1 can be described as follows: A delivery is made to location 1, the leftover stock from this delivery is taken to location 2. The vehicle returns to the depot to restock before completing delivery of the remaining packages to location 2. Rather than taking the small number of leftover packages to location 3, the vehicle returns to the depot to restock before completing the round trip to location 3. Since the vehicle returned to the depot between locations 2 and 3, the solution is effectively split into two independent delivery groups/routes, one shown with blue arrows and the other with red. We can represent this solution as a set of subsets: {(1, 2), (3)}. The first subset corresponds to the route shown with blue arrows and the second corresponds to the route shown with red arrows. Note that the order of these subsets does not affect the quality of the solution, only the order of elements within each subset.

C. Solution Space

The representation of the above example solution as an unordered set of ordered subsets can be extended to capture the entire solution space. Separation of the locations into independent delivery groups in all possible ways is akin to taking all possible partitions of the full set {1, 2, …, n}. The order in which locations can be visited within each delivery group must also be taken into account. So for each partition, the full set of solutions it represents can be acquired by taking all combinations of the permutations of each subset. In other words, by first expanding each subset in the partition into a set containing all permutations of its elements, then by taking the Cartesian product of the resulting sets of permutations. The resulting solution space is the set of all possible partitions of the n elements into nonempty and totally ordered subsets.

III The Quantum Walk Optimisation Algorithm

A. Theoretical Framework

The Quantum Walk Optimisation Algorithm (QWOA) [2, 3] was designed to identify optimal, or near-optimal, solutions to combinatorial optimisation problems. Formally, we consider a mapping $f:\mathbb{S}\to \mathbb{R}$, which returns a measure of the cost (or “quality”) associated with each possible solution in the solution space $\mathbb{S}$, where $\mathbb{S}$ has cardinality M.

The starting point of the QWOA is a quantum system with M basis states, one for each solution in $\mathbb{S}$, initialised in an equal superposition,

$\left|s\right.⟩=\frac{1}{\sqrt{M}}\sum _{x\in \mathbb{S}}\left|x\right.⟩.(2)$

This initial state is then evolved through repeated application of the quality-dependent phase-shift and quantum-walk-mixing unitaries. The quality-dependent phase-shift unitary

$U_Q\left({\gamma }_j\right)=\exp \left(-\text{i}{\gamma }_{j}Q\right),(3)$

where ${\gamma }_j\in \mathbb{R}$ and Q is a diagonal operator such that $Q\left|x\right.⟩=f(x)\left|x\right.⟩$. The quantum-walk-mixing unitary is defined as

$U_W\left(t_j\right)=\exp \left(-\text{i}{t}_j\mathcal{L}\right),(4)$

where t_j ≥ 0, and $\mathcal{L}$ is the Laplacian matrix of a circulant graph that connects the feasible solutions to the problem. For the purpose of this work, we choose a complete graph that connects computational basis states that correspond to a valid solution to the problem. Specifically, the dimension of the Hilbert space is always a power of two whilst the total number of valid routing solutions is not, in general, a power of two. In order to ensure a valid solution to the problem, it is necessary to restrict the connectivity to only the computational basis states that correspond to valid solutions to the problem. Thus, the Laplacian is defined as

$\left.⟨x\right|\mathcal{L}\left|y\right.⟩=\left\{\begin{array}{cc}0\hfill & x\text{or}y\text{are\hspace{0.17em}not\hspace{0.17em}solutions},\hfill \\ M-1,\hfill & x=y,\hfill \\ -1,\hfill & x\ne y.\hfill \end{array}\right.(5)$

Note that the initial state, $\left|s\right.⟩$, can be understood as the complete graph connecting all solutions in the solution space, where each solution/node is occupied with equal probability. The first unitary U_Q applies a phase-shift at each node proportional to the cost/quality of the solution at that node, with the proportionality constant given by the parameter, γ_j. The second unitary U_W can be understood as performing a quantum walk over the complete graph for time t_j, mixing the amplitudes across nodes. Following the mixing of phase-shifted amplitudes across the nodes of the complete graph, constructive and destructive interference will result in quality-dependent amplitude amplification, controlled by the parameters γ_j and t_j. Application of U_Q and U_W is repeated r times such that the final state of the system is given by,

$\left|\mathit{\gamma },\mathit{t}\right.⟩=U_{QWOA}\left(\mathit{\gamma },\mathit{t}\right)\left|s\right.⟩=U_W\left(t_r\right)U_Q\left({\gamma }_r\right)\dots ,U_W\left(t_1\right)U_Q\left({\gamma }_1\right)\left|s\right.⟩,(6)$

where t = (t₁, t₂, …, t_r) and γ = (γ₁, γ₂, …, γ_r).

By tuning the parameters γ and t, it is possible to amplify the amplitudes corresponding with low cost solutions, and therefore increase the probability of a measurement of the system collapsing it into a low cost solution. The process of tuning the parameters is conducted iteratively through the use of a classical optimisation algorithm (e.g. Nelder-Mead) which takes as its objective function the expectation value of the Q operator:

$c\left(\mathit{\gamma },\mathit{t}\right)=\left.⟨\mathit{\gamma },\mathit{t}\right|Q\left|\mathit{\gamma },\mathit{t}\right.⟩.(7)$

The QWOA is therefore a hybrid quantum (amplitude amplification) and classical (variational) approach. The extent of amplitude amplification possible is restricted by the number of iterations, so an increasing r presents the opportunity for greater amplification of the amplitudes corresponding to the optimal or sub-optimal solutions in the final state $\left|\mathit{\gamma },\mathit{t}\right.⟩$.

B. Quantum Circuits

Figure 2 illustrates the overall quantum circuit layout for the QWOA, with each iteration applying first the quality-dependent phase-shift and then the quantum-walk-mixing unitary. This circuit will complete r iterations for a depth-r QWOA, leading to 2r independent classical variational parameters. The expectation value of the solution quality is used to tune the parameters and is obtained by repeated sampling of the output state. This sampling process is efficient in obtaining a precise estimation of the expectation value for any problem in the NPO-PB complexity class as discussed in [2].

FIGURE 2

FIGURE 2. Schematic diagram of the QWOA circuit paired with a classical optimiser. The classical optimiser is used to tune the phase-shift and walk-time parameters, γ and t, in order to produce an optimally amplified state, $\left|\mathit{\gamma },\mathit{t}\right.⟩$, with a low expectation value for cost c(γ, t). The circuit performs $\left|\mathit{\gamma },\mathit{t}\right.⟩\left|0\right.⟩=U_W(t_r)U_Q({\gamma }_r)\dots U_W(t_1)U_Q({\gamma }_1)\left|s\right.⟩\left|0\right.⟩$. The graph Laplacian $\mathcal{L}$ encodes the connectivity between feasible solutions. The operator f evaluates the quality of a solution to k bits of precision and implements the costing algorithm presented in Section IV.

The QWOA utilises a quantum walk over the feasible combinatorial solutions, connected using an arbitrary choice of graph. The requirement on the graph is that there exists an efficient quantum circuit to simulate the quantum walk. We choose to connect the M valid solutions using the complete graph ${\mathbb{K}}_M$, where M is not necessarily a power of 2. The arbitrary-modulus Quantum Fourier Transform (QFT) was shown to be an efficient approach to perform quantum walk on this graph, and further was applicable to any other choice of circulant graph [3]. However, the QFT is known to be sensitive to noise [11], which is an undesirable property for NISQ computation. In terms of the QWOA this implies potential degradation of the quantum walk, causing leakage of amplitude to infeasible solutions. Here, we present an improved highly efficient quantum circuit for fast-forwarded quantum walk simulation on complete graphs connecting the feasible solutions.

The Laplacian described in Eq. 5 can be expressed as $\mathcal{L}=M(\mathbb{I}-\left|s\right.⟩\left.⟨s\right|)$, where $\left|s\right.⟩$ is the equal superposition over the feasible solutions as per Eq. 2. Thus the quantum walk can be expressed as, up to a global phase,

$e^{-it\mathcal{L}}=\mathbb{I}+\left(-1+e^{iMt}\right)\left|s\right.⟩\left.⟨s\right|.(8)$

This expresses the rotation of the system about a specific state $\left|s\right.⟩$. The rotation about a specific state is well-known in the context of deterministic quantum search and optimal fixed-point amplitude amplification [12–14]. We can therefore directly utilise these results to design an exact circuit for a quantum walk on the complete graph (which connects only the feasible solutions). The first step is to prepare an equal superposition state over all feasible solutions. We do this by preparing the equal superposition over the first M integers $(1/\sqrt{M}){\sum }_{y=0}^{M-1}\left|y\right.⟩$, and then use a specific un-indexing unitary $U_{\#}^{†}$ to map each unique index state $\left|i\right.⟩$ to the ith feasible solution state $\left|i{d}^{-1}(y)\right.⟩$, thus producing the equal superposition over the feasible solution states $\left|s\right.⟩$.

We label this unitary to prepare the desired state $\left|s\right.⟩$ starting in the $\left|0\right.⟩$ state as $\mathcal{G}$, as shown in Figure 3. A similar approach was used by Chiew et al. [15] for preparing superpositions over set permutations. Let m = ⌈ log ₂M⌉ be the number of qubits, and 2^m the dimension of the corresponding Hilbert space. Clearly if M is a power of 2, we can set $\mathcal{G}=H^{\otimes n}$. Otherwise, define

$\mathcal{G}=U_{\#}^{†}{H}^{\otimes m}{S}_0\left(\theta \right)H^{\otimes m}{S}_{\chi }\left(\theta \right)H^{\otimes m},(9)$

where

$S_0\left(\theta \right)=\left\{\begin{array}{cc}{e}^{i\theta }\left|j\right.⟩\hfill & j=0,\hfill \\ \left|j\right.⟩\hfill & 1\le j<2^m,\hfill \end{array}\right.(10)$

and

$S_{\chi }\left(\theta \right)=\left\{\begin{array}{cc}{e}^{i\theta }\left|j\right.⟩\hfill & 0\le j<M,\hfill \\ \left|j\right.⟩\hfill & M\le j<2^m.\hfill \end{array}\right.(11)$

FIGURE 3

FIGURE 3. (A) Circuit for preparing the equal superposition over feasible solutions such that $\mathcal{G}\left|0\right.⟩\left|0\right.⟩=\left|s\right.⟩\left|0\right.⟩$, where $\theta =2\arcsin \sqrt{\frac{2^m}{4M}}$ with m = ⌈ log ₂M⌉. Here the $<M$ gate is an integer comparator that conditionally toggles the ancilla if the input register holds a value less than M. (B) Circuit for exact simulation of quantum walk on the complete graph connecting feasible solutions, where ${\varphi }_t=-\frac{Mt}{2}$, by rotating about the $\left|s\right.⟩$ state as suggested by Yoder et al. [14].

It can be directly verified that by choosing $\theta =2\arcsin \sqrt{\frac{2^m}{4M}}$, we have $\left|s\right.⟩=\mathcal{G}\left|0\right.⟩$ as required. Thus, utilising the fixed-point search quantum circuit given in [14], we give a quantum circuit in Figure 3 that exactly simulates quantum walk on the complete graph connecting the feasible solutions.

Note that the spatial complexity of this QWOA circuit is dominated by that of the input register, which will require m = ⌈ log ₂M⌉ qubits, where M, the cardinality of the solution space, is given in Section IVA. The time complexity of the quantum walk section of the circuit (besides the indexing process) is dominated by the use of the quantum comparator, leading to a gate complexity of $\mathcal{O}(m\log m)$ [16]. However, the indexing and un-indexing algorithms will dominate the overall circuit complexity, with a gate complexity of $\mathcal{O}(poly(n))$. Note that these algorithms are discussed in Section IVB and have classical time complexity of $\mathcal{O}(n^2)$. It is assumed that a quantum circuit implementation would therefor have gate complexity of at least $\mathcal{O}(n^2)$ but not greater than $\mathcal{O}(poly(n))$

IV Prerequisite Features for QWOA

There are four primary criteria for a problem to be considered as an adequate candidate for QWOA [3, 4]. Specifically,

1. The cardinality of the solution space (number of unique solutions) M must be efficiently computable, given a specific problem instance.

2. There must exist an efficient indexing and un-indexing algorithm for the solution space. The indexing algorithm should be capable of taking any arbitrary solution and returning its location in an ordered list of the full solution space. The un-indexing algorithm should be able to do the reverse of this process, returning the unique solution corresponding to a given index. Both of these processes should be efficient, in other words, they can not simply involve the construction of the entire solution space and their complexity should be at most polynomial relative to the size of the problem.

3. The quality of any arbitrary solution must be efficiently computable.

4. As per [2], the problem must lie in the NPO-PB class (the class of optimisation problems where the objective function is bounded polynomially in the problem size). This requirement comes from sampling the quantum circuit to estimate the expectation value, to ensure the number of samples required does not grow exponentially.

The CVRP satisfies these criteria as discussed below.

A. Cardinality of the Solution Space

As described earlier, the full solution space for a problem of size n can be represented by the set of all possible partitions of n elements into nonempty and totally ordered subsets. This is closely related to the unsigned Lah numbers, L(n, k), which count the number of ways a set of n elements can be partitioned into k nonempty and totally ordered subsets [17]. It follows logically that the number of solutions, M, in the solution space for a problem of size n would be given by the sum of unsigned Lah numbers:

$M=L\left(n\right)=\sum _{k=1}^{n}L\left(n,k\right).(12)$

The unsigned Lah Numbers can be computed with the equation:

$L\left(n,k\right)=\left(\genfrac{}{}{0ex}{}{n-1}{k-1}\right)\frac{n!}{k!}.(13)$

Alternatively, the unsigned Lah numbers can be computed using the recursive relation [18]:

$L\left(n,k\right)=\left\{\begin{array}{cc}1,\hfill & \text{if\hspace{0.17em}}n=k\hfill \\ 0,\hfill & \text{if\hspace{0.17em}}n,k=0\text{\hspace{0.17em}or\hspace{0.17em}}k>n\hfill \\ L\left(n-1,k-1\right)+\left(n+k-1\right)L\left(n-1,k\right),\hfill & \text{otherwise}.\hfill \end{array}\right.(14)$

As such, the cardinality of the solution space for the CVRP is efficiently computable for a problem of any size.

B. Indexing and Unindexing Algorithms

The application of the QWOA to any particular problem requires that its solution space be efficiently indexable. The logic required to produce such an indexing algorithm for the solution space of a size n CVRP can be found by closely analysing the recursive relation depicted in Eq. 14.

Solutions can be classified according to k, the number of delivery groups (or subsets) that they contain, with k taking values: 1, 2, …, n. As such, the solution space can be divided into n distinct sub-spaces, one for each of these classifications. Each of these sub-spaces can be further divided into those solutions where the n ^th element exists in a singleton, i.e. in a subset by itself, and those where it exists in a subset with other elements.

For a given solution sub-space with k = K, the number of solutions contained is given by L(n, K) = L(n − 1, K − 1) + (n + K − 1)L(n − 1, K) (see Eq. 14). With careful consideration, it can be seen that the first term in this expression corresponds to the group of solutions for which the n ^th element is in a singleton: by removing the n ^th element from every solution, a subset is lost, and the solution sub-space which remains consists of all ordered partitions of n − 1 elements into K − 1 subsets, of which there are L(n − 1, K − 1).

Similarly, the second term corresponds to the group of solutions for which the n ^th element is not in a singleton. This group of solutions can be acquired by taking the solution sub-space consisting of all partitions of n − 1 elements into K ordered subsets, of which there are L(n − 1, K) solutions, and for each of these solutions, inserting the n ^th element into each of the (n + K − 1) available locations.

This is the logic upon which the algorithm is built, which consists of asking several questions for any particular solution:

• How many subsets exist within the solution?

• Is the n ^th element in a singleton?

• If not, in which of the (n + K − 1) available locations is it situated?

• After removing the n ^th element, where is the next element located?

The process of indexing or unindexing the solution space can also be understood as navigating a decision tree that is constructed from this same logic [19]. An example decision tree for the problem size n = 4 is shown in Figure 4.

FIGURE 4

FIGURE 4. A compact decision tree for the CVRP of size n = 4. The lower portion of the diagram shows the general structure and conditional logic which forms the basis for the construction of the tree. The branches labelled with a number represent multiple branches, with p = 0 corresponding with the lower branch and so on.

Each terminal node (those at the right-hand side in Figure 4) represents a solution in the solution space. The solution located at the bottom has index zero, and the indices increase incrementally up the tree, such that the index of any particular solution is equal to the number of solutions that exist below it. Note that the tree has been included in a compact form, though there appear to be only eight terminal nodes, there are in fact 73. This is because the branches labelled with a number actually represent multiple branches, with the lowest branch being associated with p = 0 and the branch above with p = 1 and so on. Navigating through the decision tree from the origin to a particular terminal node involves choosing particular branches, each of which represents making decisions about the number of subsets in a solution, or the location of each element within the ordered subsets of a solution. By summing the number of solutions that are passed over with each decision made or branch taken, the resulting sum will be equal to the index of the final solution. In this way, the indexing process can be understood as navigating through the tree. Similarly, the unindexing process can be understood as choosing the appropriate branches such that the number of solutions passed over is equal to the index, with the terminal node corresponding to the required solution.

The indexing and unindexing algorithms are shown in Algorithm 1 and Algorithm 2. The time complexities of both the indexing and unindexing algorithms scale with n², $\mathcal{O}(n^2)$, which satisfies the requirement that both of these algorithms should be efficiently computable. Since efficient classical algorithms for indexing and un-indexing exist, these algorithms can be directly “quantised” by translation to a reversible circuit. This defines the indexing unitary U_# used in Figure 3.

Algorithm 1 This indexing algorithm returns the index of any solution, soln, from the solution space. It is composed of two parts, the first relates to placing the solution within a particular sub-space, determined by the number of subsets it contains, and the second is responsible for indexing the solution within this sub-space. As a clarification, on line 14, p is assigned the location of the element e within the solution/partition. This is the location where e would have been placed in the partition not containing e, where the left-most location corresponds to p = 0. This can also be calculated from the sum of the following: the number of subsets prior to the one containing e, the number of elements contained within these prior subsets, and the location of e within its own subset.

Algorithm 2 This unindexing algorithm returns the solution to which a given index within the solution space of particular problem size, n, refers. The algorithm is composed of two parts, the first is responsible for selecting the sub-space to which the indexed solution belongs, and its corresponding sub-index within that space and the second is responsible for returning the solution to which the sub-index and hence index refers.

C. Computation of Solution Qualities

With the solution indexing now well established, the next challenge is to ensure the quality of an arbitrary solution can be efficiently computed and that the objective function is polynomially bounded in the number of locations. Computing the quality/cost corresponding to any given solution involves navigating the delivery network in the order specified by the solution while comparing the number of packages on the vehicle and its capacity with the number of packages required at each location/node. Each time an edge is traversed between two nodes, the corresponding cost from the cost matrix is incurred, and the total cost accumulated by the time the vehicle returns to the depot for the final time is the cost/quality of the solution. This process and the conditional logic involved is shown in detail in Algorithm 3, Algorithm 4, Algorithm 5. The complexity of this computation scales linearly with increasing problem size, n, so the quality of any arbitrary solution is efficiently computable, as required. In addition, for constant bounds on elements in the cost matrix, package vector, and the vehicle capacity, the maximum cost scales linearly with the number of locations and hence the objective function is sufficiently well-behaved for expectation sampling.

Algorithm 3 Computing Solution Qualities: This algorithm returns the cost of a particular solution, soln, for a particular instance of the CVRP, defined by the parameters: vehicle capacity, V, package vector, P, and cost matrix, C. Note that the variable LO is short for “leftovers,” and tracks the number of packages remaining on the vehicle. Also, the variable perm is short for “permutation,” and RH is short for “return home.” RH is a Boolean variable that tracks whether the vehicle is due to return to the depot or has just returned to the depot.

Algorithm 4 Computing Solution Qualities

Algorithm 5 Computing Solution Qualities

V Numerical Results

An example CVRP of size n = 8 was created from randomly generated values and is captured by the following paramaters:

$P=\left[\begin{array}{c}\hfill 23\hfill \\ \hfill 18\hfill \\ \hfill 28\hfill \\ \hfill 7\hfill \\ \hfill 23\hfill \\ \hfill 27\hfill \\ \hfill 9\hfill \\ \hfill 22\hfill \end{array}\right],C=\left[\begin{array}{ccccccccc}\hfill 0\hfill & \hfill 10\hfill & \hfill 16\hfill & \hfill 10\hfill & \hfill 14\hfill & \hfill 17\hfill & \hfill 12\hfill & \hfill 11\hfill & \hfill 17\hfill \\ \hfill 10\hfill & \hfill 0\hfill & \hfill 7\hfill & \hfill 8\hfill & \hfill 14\hfill & \hfill 9\hfill & \hfill 4\hfill & \hfill 1\hfill & \hfill 5\hfill \\ \hfill 16\hfill & \hfill 7\hfill & \hfill 0\hfill & \hfill 15\hfill & \hfill 10\hfill & \hfill 10\hfill & \hfill 5\hfill & \hfill 2\hfill & \hfill 11\hfill \\ \hfill 10\hfill & \hfill 8\hfill & \hfill 15\hfill & \hfill 0\hfill & \hfill 5\hfill & \hfill 15\hfill & \hfill 13\hfill & \hfill 15\hfill & \hfill 15\hfill \\ \hfill 14\hfill & \hfill 14\hfill & \hfill 10\hfill & \hfill 5\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 4\hfill & \hfill 15\hfill & \hfill 4\hfill \\ \hfill 17\hfill & \hfill 9\hfill & \hfill 10\hfill & \hfill 15\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 13\hfill & \hfill 5\hfill & \hfill 3\hfill \\ \hfill 12\hfill & \hfill 4\hfill & \hfill 5\hfill & \hfill 13\hfill & \hfill 4\hfill & \hfill 13\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 7\hfill \\ \hfill 11\hfill & \hfill 1\hfill & \hfill 2\hfill & \hfill 15\hfill & \hfill 15\hfill & \hfill 5\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill 2\hfill \\ \hfill 17\hfill & \hfill 5\hfill & \hfill 11\hfill & \hfill 15\hfill & \hfill 4\hfill & \hfill 3\hfill & \hfill 7\hfill & \hfill 2\hfill & \hfill 0\hfill \end{array}\right],V=20.$

The inter-location costs consist of randomly generated integers from the uniform distribution from 1 to 15. Similarly, the depot to location costs are integers between 10 and 20, and the package counts between 5 and 30.

The cardinality of the solution space for an n = 8 problem, calculated as per Eq. 13, is L(8) = 394, 353. Assessing the quality of these solutions reveals that there exist only 148 distinct solution qualities, where the lowest values correspond to the highest quality solutions. There is therefore a large amount of degeneracy in the qualities, which can be seen in the distribution of qualities shown in Figure 5A.

FIGURE 5

FIGURE 5. (A) Initial quality distribution of the example problem-solution space and, (B), the evolution of the probability distribution, relative to solution quality, with increasing iterations, r, of phase-shifts and walks.

The distribution of qualities shown in Figure 5A is similar in nature to a normal distribution, though the distribution is skewed slightly towards the higher values (lower quality solutions). The solution space is populated predominantly by solutions of mid-tier quality, with only a small number of solutions of high quality.

In order to assess the capability of the QWOA to handle the example problem outlined above, the relevant quantum system and algorithm is simulated numerically. This is carried out using the software package, QuOP_MPI [20] in accordance with the framework laid out in Section IIIA. The phase-shift and walk-time parameters are optimised with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, implemented via the SciPy open-source Python library.

Shown in Figure 5A is the initial quality distribution of the example problem described above, prior to the application of any phase-shifts or quantum walks. Even though the initialisation of the system as an equal superposition means each node on the complete graph is associated with an equal part of the overall probability density, the mid-tier qualities are over-represented in the initial quality distribution because there are a large number of nodes corresponding to these qualities. The success of the QWOA will be measured by its ability to evolve the state of the system, via phase-shifts and quantum walks, such that the probability densities concentrate at the nodes corresponding with the highest quality solutions.

The QWOA objective function is the expectation value of the measured solution quality after r QWOA iterations, given by Eq. 7. As shown in Figure 5B, with an increasing number of repeated iterations, r, the probability distribution within the system does in fact concentrate towards higher quality solutions. Upon closer inspection, however, for the given r values, the resulting quality distribution is not concentrating at the most optimal solutions as much as it is at the near-optimal solutions. With significantly larger r, it is expected that convergence towards the most optimal solutions would become more complete, but as shown in Figure 6, the rate of convergence decreases at large r. More interesting is a direct comparison with classical random sampling of equivalent computational effort, where the classical data represents the expected best quality measured from 2r classical random samples of the solution space, computed from 100,000 trials. Note that we compare with 2r classical samples, because this represents the same number of calls to the quality function when compared with r iterations of the QWOA process, where we quantify computational effort by the number of calls to the quality function. It should also be clarified that no allowance has been made for the computational effort involved in the classical optimisation procedure to arrive at an optimal set of parameters, γ and t, instead, the focus has been on the final optimally amplified state. This serves as a proof of concept in that it shows that the QWOA procedure is capable of providing speed up via sampling of the amplified state, given that there exists a computationally efficient method to arrive at a set of optimal parameters. Figure 6 shows that the expected quality/cost measured from a QWOA amplified state converges towards the target/minimum value as $\frac{1}{r^{0.45}}$, while the classical random sampling scales as $\frac{1}{r^{0.27}}$, when considering equivalent computational effort.

FIGURE 6

FIGURE 6. The convergence of the QWOA objective function towards the minimum/target value with increasing r as compared to a classical sampling method of equivalent computational effort. The QWOA curve scales as 1/r^0.45, while the classical curve scales as 1/r^0.27.

Taking a closer look at how the probability densities are being amplified relative to quality may reveal part of the reason why the system seems to concentrate primarily towards sub-optimal solutions. The graph in Figure 7 shows how the probability densities are amplified at any particular node as a function of its solution quality for the r = 10 and r = 35 case. It is perhaps not surprising that the graph in Figure 5A shows more probability density at the near-optimal solutions, because there are significantly more of these solutions present in the solution space. The probability amplification shown in Figure 7, on the other hand, is given by the ratio of final to initial probability at any given node, so the number of solutions at a given quality does not have a multiplicative effect on the data shown. It is therefore quite telling that here the near-optimal solutions are being amplified more strongly than the most optimal solutions. A potential explanation for this is that amplification of near-optimal solutions is favoured by the optimisation process over the most optimal solutions, because although they are not as optimal, they have a larger influence in minimising the objective function due to their superior numbers.

FIGURE 7

FIGURE 7. This graph quantifies the probability amplification at any particular node/solution as a function of its quality for the cases where r = 10 and r = 35.

Any way in which the most optimal solutions can be weighted more heavily with regards to their effect on the objective function would be beneficial in counteracting this potential bias towards sub-optimal solutions. However, this may not be effective in a practical sense, because the expectation value is calculated from repeated observations of the system, and in order for the increased weighting to take effect, the most-optimal solutions must actually be measured, which is not likely in the pre-amplified state.

In the cannon of currently discovered quantum algorithms, QWOA is most naturally compared to the Quantum Approximation Optimisation Algorithm (QAOA). This algorithm seeks to solve combinatorial optimisation problems following a similar alternating ansatz that is capable of solving the CVRP as formulated in this study. However, QAOA is limited to a mixing operation across the entire Hilbert space [4], with invalid states being differentiated via a penalty function. In contrast, the QWOA indexing unitary restricts its equivalent walk operation to the subspace of only valid solutions. A previous numerical study demonstrated that this allows for the QWOA to produce more effective convergence to high-quality solutions than QAOA and QAOAz at the same circuit depth [4].

As QWOA seeks approximate solutions to the CVRP combinatorial optimisation problem, it is comparable in intent to classical heuristic and metaheuristic algorithms. Such methods are numerous and, as with QWOA, exhibit effectiveness that often depends on the structure of the underlying CVRP instance [21]. Nevertheless, current methods can identify solutions within 0.5–1% of the optimum for problems with hundreds to thousands of delivery locations [22]. The effectiveness of these algorithms is determined by benchmarking their performance against established data sets [21, 22], which exceed the scale afforded by classical simulation of quantum dynamics. For these reasons, a comparison between this QWOA application and current classical heuristic and metaheuristic methods lies beyond the scope of this study.

VI Conclusion

To conclude, the quantum-walk-based optimization algorithm (QWOA) can be effectively applied to the Capacitated Vehicle Routing Problem. Efficient algorithms have been developed for the indexing and unindexing of the Solution space for this problem. It has been demonstrated for a randomly generated eight location problem, that the QWOA is capable of significantly amplifying probabilities for optimal and near-optimal solutions, and in doing so, achieves expected qualities with less computational effort than that required by classical random sampling. The rate of convergence towards the most optimal solutions is however still somewhat limited, and the speed-up relative to classical sampling is contingent on navigating the variational process efficiently. Further work will be focused on quantifying and minimising the computational effort involved in the variational process and in finding a way to improve the algorithm further to provide faster convergence towards the most optimal solutions.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

TB is the first author of this paper; he performed the initial literature survey, developed the main algorithms, performed the numerical simulations, and produced the first draft. EM developed the QuOp package that was used to perform the calculations. SM provided much of the theoretical framework and detailed quantum circuits. JW initiated the project and provided overall guidance and supervision. All authors were involved in the editing of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. EM and SM acknowledge the support of the Australian Government Research Training Program Scholarship. SM’s research was also supported by a Hackett Postgraduate Research Scholarship at the University of Western Australia. JW would like to thank Shakib Daryanoosh for their useful discussions.

References

1. Aaronson S Guest Column: NP-Complete Problems and Physical Reality. SIGACT News (2005) 36:30–52. doi:10.1145/1052796.1052804

CrossRef Full Text | Google Scholar

2. Marsh S, Wang JB A Quantum Walk-Assisted Approximate Algorithm for Bounded NP Optimisation Problems. Quan Inf Process (2019) 18:61. doi:10.1007/s11128-019-2171-3

CrossRef Full Text | Google Scholar

3. Marsh S, Wang JB Physical Review Research, 2 (2020). p. 023302. doi:10.1103/physrevresearch.2.023302

CrossRef Full Text | Google Scholar

4. Slate N, Matwiejew E, Marsh S, Wang JB Quantum Walk-Based Portfolio Optimisation. Quantum (2021) 5:513. doi:10.22331/q-2021-07-28-513

CrossRef Full Text | Google Scholar

5. Dantzig GB, Ramser JH The Truck Dispatching Problem. Manage Sci (1959) 6:80–91. doi:10.1287/mnsc.6.1.80

CrossRef Full Text | Google Scholar

6. Toth P, Vigo D Models, Relaxations and Exact Approaches for the Capacitated Vehicle Routing Problem. Discrete Appl Maths (2002) 123:487–512. doi:10.1016/s0166-218x(01)00351-1

CrossRef Full Text | Google Scholar

7. Harikrishnafkumar R, Nannapaneni S, Nguyen NH, Steck JE, Behrman EC arXiv:2005.12478 (2020).

Google Scholar

8. Borowski M, Gora P, Karnas K, Błajda M, Król K, Matyjasek A, et al. International Conference on Computational Science. In: New Hybrid Quantum Annealing Algorithms for Solving Vehicle Routing Problem. Springer (2020). p. 546–61. doi:10.1007/978-3-030-50433-5_42

CrossRef Full Text | Google Scholar

9. Syrichas A, Crispin A Large-scale Vehicle Routing Problems: Quantum Annealing, Tunings and Results. Comput Operations Res (2017) 87:52–62. doi:10.1016/j.cor.2017.05.014

CrossRef Full Text | Google Scholar

10. Feld S, Roch C, Gabor T, Seidel C, Neukart F, Galter I, et al. A Hybrid Solution Method for the Capacitated Vehicle Routing Problem Using a Quantum Annealer. Front ICT (2019) 6:13. doi:10.3389/fict.2019.00013

CrossRef Full Text | Google Scholar

11. Shor PW Scheme for Reducing Decoherence in Quantum Computer Memory. Phys Rev A (1995) 52:R2493–R2496. doi:10.1103/physreva.52.r2493

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Long GL Phys. Rev. A, 64 (2001). p. 022307. doi:10.1103/physreva.64.022307

CrossRef Full Text | Google Scholar

13. Brassard G, Høyer P, Mosca M, Tapp A Quantum Amplitude Amplification and Estimation. Contemp Math (2002) 305:53–74. doi:10.1090/conm/305/05215

CrossRef Full Text | Google Scholar

14. Yoder TJ, Low GH, Chuang IL Fixed-Point Quantum Search with an Optimal Number of Queries. Phys Rev Lett (2014) 113:210501. doi:10.1103/physrevlett.113.210501

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Chiew M, de Lacy K, Yu CH, Marsh S, Wang JB Graph Comparison via Nonlinear Quantum Search. Quan Inf. Process. (2019) 18:302. doi:10.1007/s11128-019-2407-2

CrossRef Full Text | Google Scholar

16. Gidney C Factoring with N+2 Clean Qubits and N-1 Dirty Qubits. (2018). arXiv [Preprint]. Available at: https://arxiv.org/abs/1706.07884 (Accessed January 19, 2018).

Google Scholar

17. Petkovšek M, Pisanski T Pi Mu Epsilon Journal, 12 (2007). p. 417.

Google Scholar

18. Bényi B, Méndez M, Ramirez JL. Australas. J. Comb (2020) 77:157–179.

19. Wilf HS A Unified Setting for Sequencing, Ranking, and Selection Algorithms for Combinatorial Objects. Adv Maths (1977) 24:281–91. doi:10.1016/s0001-8708(77)80046-7

CrossRef Full Text | Google Scholar

20. Matwiejew E Edric-Matwiejew/QuOp_mpi: v0.0.2 (2020) arXiv:2110.03963. doi:10.5281/zenodo.5540364

CrossRef Full Text | Google Scholar

21. Sánchez FFS, Lazo CAL, Quiñónez FYS Comparative Study of Algorithms Metaheuristics Based Applied to the Solution of the Capacitated Vehicle Routing Problem (IntechOpen, 2020) Publication Title: Novel Trends in the Traveling Salesman Problem. D Davendra, and BD Magdalena. Editor (2020) 63–88. doi:10.5772/intechopen.91972

CrossRef Full Text | Google Scholar

22. Vidal T, Crainic TG, Gendreau M, Prins C A Unified Solution Framework for Multi-Attribute Vehicle Routing Problems. Eur J Oper Res (2014) 234:658–73. doi:10.1016/j.ejor.2013.09.045

CrossRef Full Text | Google Scholar