All optimization problems in the complexity class NP can be reformulated as the problem of finding the ground state (lowest-energy configuration) of a Hamiltonian48. This is also the method we use for the HVRP in this work. Since the HVRP combines two distinct problems, a routing problem and a capacity problem, we have to derive an Ising Hamiltonian that captures both these problems simultaneously.
Routing problem
For the routing problem, we start from the TSP formulation in Ref.48 with the Hamiltonian \(H^\textrmTSP = H_A^\textrmTSP + H_B^\textrmTSP\), where
$$\beginaligned H_A^\textrmTSP&= A \sum _i=1^N \left( 1 – \sum _\alpha =1^N y_i \alpha \right) ^2 + A \sum _\alpha =1^N \left( 1 – \sum _i=1^N y_i \alpha \right) ^2 + A \sum _(i, j) \notin \mathscr E \sum _\alpha =1^N y_i \alpha y_j \alpha +1 \, , \endaligned$$
(10)
$$\beginaligned H_B^\textrmTSP&= B \sum _(i, j) \in \mathscr E W_i j \sum _\alpha =1^N y_i \alpha y_j \alpha +1 \, , \endaligned$$
(11)
encoding the total cost. Here \(N = |\mathscr N|\) is the number of nodes including the depot, A and B are positive constants, and W encodes the distances between the nodes. The index i represents the nodes and \(\alpha\) the order in a prospective cycle. The binary variables \(y_i \alpha \) can be referred to as ’routing variables’ indicating in which order of the cycle node i is visited. There are \(N^2\) variables, with \(y_i,N+1\equiv y_i,1\) for all i, such that the route ends where it starts. Minimizing the first term in Eq. (10) leads to each node being visited exactly once and minimizing the second term in Eq. (10) ensures that each cycle has exactly one entry for each order. The last term in Eq. (10), which ensures that non-existent arcs are not used, can be neglected for the problem we investigate because we assume a fully connected graph. The Hamiltonian \(H_B^\textrmTSP\) in Eq. (11) adds up the total distance travelled.
Above, we used \(y_i \alpha ^v\) as a routing variable that tracks the order \(\alpha\) in which node i is visited within a cycle. To blend this framework with the mathematical structure of the HVRP, formulated in Eqs. (1)–(9), we introduce a transformation from the decision variable y to x as follows:
$$\beginaligned x_ij^v&= \sum _\alpha =1^N_0-1 y_i\alpha ^v y_j\alpha +1^v \, , \endaligned$$
(12)
$$\beginaligned x_0i^v&= y_i1^v + \sum _\alpha =2^N_0 \left( 1 – \sum _\beginarrayc j=1 \\ j \ne i \endarray^N_0 y_j \alpha -1^v \right) y_i \alpha ^v \, , \qquad x_i0^v = y_i N_0^v + \sum _\alpha =1^N_0 – 1 y_i \alpha ^v \left( 1 – \sum _\beginarrayc j=1 \\ j \ne i \endarray^N_0 y_j \alpha +1^v \right) \, . \endaligned$$
(13)
The summation in Eq. (12) is not equal 0 if and only if i and j are subsequent stops on the same route. Equation (13) ensures that the first and last stops are automatically connected to the depot (assuming a single depot). By using the decision variables \(y_i \alpha ^v\) that indicate the position in a prospective cycle instead of \(x_ij^v\) that is 1 if and only if a vehicle traverses from customer i to j38, we take care of the subtour-elimination constraint encoded in Eq. (6), which otherwise allows for closed loops between customers, without returning to the depot.
The transformation from \(y_i \alpha ^v\) to \(x_ij^v\) as defined in Eqs. (12)–(13) is unique. However, the inverse transformation, i.e., from the HVRP model variables \(x_ij^v\) to the quantum model variables \(y_i \alpha ^v\), is not unique. This becomes evident if we consider that any permutation of the sequence \(y_i \alpha ^v\) for a fixed i will yield the same \(x_ij^v\) values. Hence, a single solution in the HVRP model could correspond to multiple valid solutions in the Ising model. Despite this non-uniqueness, every feasible solution in the HVRP model has a corresponding set of solutions in the Ising model, ensuring the robustness and validity of our approach.
We can now write the Ising formulation for the routing problem. Here we extend the formulation compared to previous works to capture different types of trucks, not just multiple trucks of the same type (having the same capacity)48,51. Let \(V = |\mathscr V|\) be the number of trucks, where \(\mathscr V\) now is the set of vehicles chosen for the optimization [instead of the set of vehicle types, as in Eqs. (1)–(9)], and denote by \(N_0 = |\mathscr N_0|\) the number of customers to visit. The indices v now represent a specific truck of a specific type [instead of just a specific type, as in Eqs. (1)–(9)]. The Ising Hamiltonian \(H = H_A + H_B + H_C\) we derive is then merging Eqs. (12)–(13) with Eqs. (1)–(9):
$$\beginaligned H_A&= A \sum _v=1^V \sum _i=1^N_0 \sum _j=1^N_0 c^v_ij \sum _\alpha =1^N_0 – 1 y_i \alpha ^v y_j \alpha + 1^v + A \sum _v=1^V \sum _i=1^N_0 (t_v + c_0i^v) \left[ y_i1^v + \sum _\alpha =2^N_0 \left( 1 – \sum _\beginarrayc j=1 \\ j \ne i \endarray^N_0 y_j \alpha -1^v \right) y_i \alpha ^v \right] \nonumber \\&\quad + A \sum _v=1^V \sum _i=1^N_0 c_i0^v \left[ y_i N_0^v + \sum _\alpha =1^N_0 – 1 y_i \alpha ^v \left( 1 – \sum _\beginarrayc j=1 \\ j \ne i \endarray^N_0 y_j \alpha +1^v \right) \right] \, , \endaligned$$
(14)
$$\beginaligned H_B&= B \sum _i=1^N_0 \left( 1 – \sum _\alpha =1^N_0 \sum _v=1^V y_i \alpha ^v \right) ^2 \, , \endaligned$$
(15)
$$\beginaligned H_C&= C \sum _\alpha =1^N_0 \left( 1- \sum _i=1^N_0 \sum _v=1^V y_i\alpha ^v \right) ^2 \, . \endaligned$$
(16)
The Hamiltonian is composed of different parts. \(H_A\) in Eq. (14) captures the total cost of the original mathematical formulation [Eq. (1)], i.e., the minimization of the variable and fixed cost. The first term estimates the variable cost for traveling between the different customers/cities, while the second and third terms measure the cost of leaving and arriving at the depot. Note that the fixed cost is included in the second term as well, since all the vehicles leaving the depot generate the fixed cost.
For this particular mapping it is necessary to define the set of vehicles that are used for the optimization beforehand. Therefore, we can neglect the inequality constraint defined in Eq. (2) from the original formulation, which ensures that the number of vehicles of a specific type does not exceed the number of available vehicles.
The constraint given by \(H_B\) in Eq. (15) ensures that each city is visited exactly once [i.e., satisfies the constraints in Eqs. (3) and (4) in the original HVRP formulation]. Furthermore, \(H_C\) in Eq. (16) ensures that each city has a unique position in the cycle and that not more than one city can be travelled to at the same time. To make sure that the constraints are not violated, we require \(0< \text max(|H_A|) < B, C\). If these inequalities are not fulfilled, we risk obtaining unphysical solutions that have lower energies than valid ones.
Visualization of a problem instance with a suggested solution (other solutions are possible; we assume a fully connected graph). The first truck \(v_1\) visits cities \(c_2\) and \(c_3\) before returning to the depot. The second truck \(v_2\) only visits city \(c_1\).
One notable technicality about the formulation is that certain solutions that may be considered valid are excluded by the constraint in Eq. (15). However, the excluded solutions are physically equivalent to some allowed solution, as illustrated by the example in Fig. 1 with two trucks and three cities:
$$\beginaligned y = \left[ \beginpmatrix 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \endpmatrix, \beginpmatrix 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \endpmatrix \right] \qquad \text and \qquad y = \left[ \beginpmatrix 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \endpmatrix, \beginpmatrix 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \endpmatrix \right] , \endaligned$$
(17)
where the matrix rows correspond to the index i (the node visited), the matrix columns correspond to the index \(\alpha\) (the order in which the nodes are visited), and each matrix corresponds to a vehicle v.
In both cases, the two different matrices describe the routes of the two different vehicles and the order in which they visit the different customers \(c_i\). They both represent a physically valid solution where the first truck visits first the second and then the third customer (\(v_1: c_2 \rightarrow c_3\)), while the second truck goes from the depot to the first customer and then back to the depot (\(v_2: c_1\)). The constraint in Eq. (15), however, rules out the latter representation as it has two non-zero entries for \(\alpha =2\). If we want to allow this larger set of viable representations, physically equivalent to solutions already allowed by Eq. (15), we can replace that constraint by a reformulated one,
$$\beginaligned H’_B = B \sum _v=1^V \sum _\alpha =1^N_0 \left( u_\alpha ^v – \sum _i=1^N_0 y_i \alpha ^v \right) ^2, \endaligned$$
(18)
where we have introduced \(N_0^2\) additional auxiliary qubits \(u_\alpha ^v\). Especially in the NISQ era, where quantum resources are scarce, it is important to encode the problem with as few qubits as possible. Thus we do not consider Eq. (18) a viable route for implementations, but use Eq. (15) for the simulations later in this article.
Capacity problem
The capacity constraint is of a similar nature as the constraints for the knapsack problem—both are described by an inequality constraint, which for the knapsack problem is to not add too many items to the knapsack and for the trucks to not overload the vehicles. Therefore, we can use the formulation given in Ref.48 to model the inequality constraint introduced by the capacities.
The knapsack problem with integer weights is the following. We have a list of N objects, labeled by i, with the weight of each object given by \(w_i\) and its value by \(c_i\). The knapsack has a limited capacity of W. The binary decision variable \(x_i\) denotes whether an item is contained (1) in the knapsack or not (0). The total weight of the knapsack is \(\mathscr W = \sum _i=1^N w_i x_i\) with a total value of \(\mathscr C = \sum _i=1^N c_i x_i\). The knapsack problem is to maximize \(\mathscr C\) under the constraint \(\mathscr W \le W\). This problem belongs to the complexity class NP-hard48.
We can write an Ising formulation of the knapsack problem as follows. Let \(z_n\) for \(1 \le n \le W\) be a binary variable which is 1 if the final weight of the knapsack is n and 0 otherwise48. The Hamiltonian \(H^\textrmK = H_A^\textrmK + H_B^\textrmK\) whose energy we seek to minimize is then
$$\beginaligned H_A^\textrmK&= A \left( 1 – \sum _n=1^W z_n \right) ^2 + A \left( \sum _n=1^W n z_n – \sum _i=1^N w_i x_i \right) ^2 \, , \endaligned$$
(19)
$$\beginaligned H_B^\textrmK&= -B \sum _i=1^N c_i x_i \, . \endaligned$$
(20)
To ensure that the hard constraint is respected, we require \(0< \text max(|H_B^\textrmK|) < A\).
Reducing the number of auxiliary qubits
It is possible to reduce the number of variables required for the auxiliary variable \(z_n\). We want to encode a variable which can take the values from 0 to W. Let \(M \equiv \lfloor \text log_2(W) \rfloor\). We then require \(M+1\) binary variables instead of W binary variables:
$$\beginaligned \sum _n=1^W n z_n \rightarrow \sum _m = 0^M-1 2^m z_m + \left( W + 1 – 2^M \right) z_M. \endaligned$$
(21)
It is important to note that in this encoding method, multiple bitstrings \(z_m\) can represent the same number. This situation arises due to the transition from one-hot encoding to a binary representation. In one-hot encoding, each integer from 0 to W would have a unique bitstring representation. However, in our binary representation, where we encode numbers using fewer bits, different bitstrings can be interpreted as the same integer. This degeneracy is possible if \(W \ne 2^M+1 -1\).
We can make use of the inequality constraint given in the knapsack formulation [Eq. (19)] to encode the capacity constraints for the HVRP. Therefore, we can neglect \(H_B^\textrmK\) [Eq. (20)] and only consider \(H_A^\textrmK\) [Eq. (19)]. Let \(Q^v\) be the maximum capacity of vehicle v. The Hamiltonian then becomes
$$\beginaligned H^\textrmK&= A \sum _v \left( 1 – \sum _k=0^Q^v z_k^v \right) ^2 + A \sum _v \left( \sum _k=0^Q^v k z_k^v – \sum _\alpha ,i q_i y_i \alpha ^v \right) ^2, \endaligned$$
(22)
or equivalently, using the log formulation,
$$\beginaligned H^\textrmK&= A \sum _v \left( \sum _k=0^M^v-1 2^k z_k^v + (Q^v + 1 – 2^M^v) z_M^v^v – \sum _\alpha ,i q_i y_i \alpha ^v \right) ^2, \endaligned$$
(23)
where \(M^v \equiv \lfloor \text log_2(Q^v) \rfloor\). Note that by using the log trick, the decision variable \(z_k^v\) switches from a one-hot encoding to a binary representation. The lowest energy possible for the Hamiltonian in Eq. (23) is zero; this minimum is reached when the \(z_k^v\) are chosen such that each vechicle v carries an amount of goods, given by the first inner sum in Eq. (23), not exceeding its maximum capacity, and that amount exactly corresponds to the demand of the customers it visits along its route, given by the second inner sum in Eq. (23).
The full Ising Hamiltonian for the HVRP
We are now ready to write down the full Hamiltonian for the HVRP. The full Ising Hamiltonian \(H_\mathscr C\) contains four terms, where the first term \(H_A\) captures the actual optimization problem and the other terms are penalty terms to ensure that invalid configurations are penalized with a high energy:
$$\beginaligned H_\mathscr C&= H_A + H_B + H_C + H_D \, , \endaligned$$
(24)
$$\beginaligned&\beginaligned H_D&= D \sum _v=1^V \left( \sum _k=0^M^v-1 2^k z_k^v + (Q^v + 1 – 2^M^v) z_M^v^v – \sum _\alpha =1^N_0 \sum _i=1^N_0 q_i y_i \alpha ^v \right) ^2 \endaligned \endaligned$$
(25)
For the terms \(H_A\) to \(H_C\), see Eqs. (14)–(16). Note that \(H_D\), together with \(H_B\), ensures that we satisfy the commodity-flow constraints in Eqs. (6) and (7) in the original HVRP formulation.
The formulation presented in this paper combines the capacity problem and the routing problem in one Ising Hamiltonian. Similarly, a unified approach is also attempted in Ref.31 with the difference that the authors add a constraint for clustering the customers as well. Here, by using the decision variables \(y_i \alpha ^v\) that indicate the position in a prospective cycle instead of \(x^v\) that is 1 if and only if a vehicle traverses from customer i to j, we circumvent the subtour-elimination constraint. This constraint needs to loop through all possible subtours as it is presented in Ref.38. Moreover, a solution obtained with our mapping does not necessarily use all the vehicles that are available. It can find the most cost-efficient subset of vehicles needed to solve the task.
Resources
The required resources (qubits) for solving the HVRP with our approach can be separated into two parts. The first part comes from representing all the possible routes between the customers in the routing part of the problem and scales with \(N_0^2 \cdot V\). Additional qubits are required for the constraining term \(H_D\). The total number of qubits required is
$$\beginaligned \#q = N_0^2 \cdot V + \underbrace\sum _v=1^V \lfloor \log _2 Q^v \rfloor + 1_H_D. \endaligned$$
(26)
Using the alternative formulation for \(H_B\) [see Eq. (18)] adds more auxiliary qubits (\(N_0 \cdot V\)), yielding
$$\beginaligned \#q =(N_0 + 1) N_0 V + \sum _v=1^V \lfloor \log _2 Q^v \rfloor + 1. \endaligned$$
(27)
Furthermore, we can estimate the number of single- and two-qubit gates required to apply the QAOA (see next section) to the problem in this formulation. For the routing part of the problem, each layer of the QAOA needs \(O (N_0^2 \cdot V)\) single-qubit gates (i.e., proportional to the number of qubits) and \(O (N_0^3 \cdot V)\) two-qubit gates (i.e., proportional to \(N_0\) times the number of qubits). For the capacity part of the problem, the number of single-qubit gates is the same, while the number of two-qubit gates is \(O (N_0^2 \cdot Q)\). What this translates to in circuit depth will depend on the connectivity of the quantum hardware.
As a comparison, we note that modern high-performance optimizers (classical computers) for the HVRP can find near-optimal solutions for problem instances with more than 1,000 customers52,53. Typically, these optimizers exhibit polynomial or pseudo-polynomial time complexity. For a quantum computer to solve problem instances of this size, it would need at least millions of controllable qubits. Systems of this size are likely still several years away from being realized.
link