{'fraction_non_alphanumeric': 0.032087148185825705, 'fraction_numerical': 0.008604611732790776, 'mean_word_length': 4.752987076322848, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Quantum Annealing Approach for Dynamic Multi-Depot Capacitated Vehicle Routing Problem.

Quantum annealing (QA) is a quantum computing algorithm that works on the principle of Adiabatic Quantum Computation (AQC), and it has shown significant computational advantages in solving combinatorial optimization problems such as vehicle routing problems (VRP) when compared to classical algorithms. This paper presents a QA approach for solving a variant VRP known as multi-depot capacitated vehicle routing problem (MDCVRP). This is an NP-hard optimization problem with real-world applications in the fields of transportation, logistics, and supply chain management. We consider heterogeneous depots and vehicles with different capacities. Given a set of heterogeneous depots, the number of vehicles in each depot, heterogeneous depot/vehicle capacities, and a set of spatially distributed customer locations, the MDCVRP attempts to identify routes of various vehicles satisfying the capacity constraints such as that all the customers are served. We model MDCVRP as a quadratic unconstrained binary optimization (QUBO) problem, which minimizes the overall distance traveled by all the vehicles across all depots given the capacity constraints. Furthermore, we formulate a QUBO model for dynamic version of MDCVRP known as D-MDCVRP, which involves dynamic rerouting of vehicles to real-time customer requests. We discuss the problem complexity and a solution approach to solving MDCVRP and D-MDCVRP on quantum annealing hardware from D-Wave.


I. INTRODUCTION
The vehicle routing problem (VRP) is a well-known combinatorial optimization problem commonly encountered in logistics, transportation, supply chain management, and scheduling [1]. There are several variants of VRP such as multi-depot vehicle routing problem (MDVRP), vehicle routing and scheduling, vehicle routing problem with time windows (VRPTW), capacitated vehicle routing problem (CVRP) and multi-depot capacitated vehicle routing problem (MDCVRP), which are computationally challenging and are NP-hard problems [2]. In this paper, we discuss a solution approach using the principles of quantum computing to improve the computational performance in solving the MDCVRP and its dynamic version, the D-MDCVRP.
A MDCVRP can be briefly defined as follows. Given a set of spatially distributed depots with heterogeneous material handling capacities, a variable set of service vehicles at each depot with heterogeneous material carrying capacity, and a set of spatially distributed customer locations, with heterogeneous demands that need to be served, the objective in an MDCVRP is to assign routes to various service vehicles to serve the customer locations considering the depot and vehicle capacities. For example, in a supply chain network for a courier collection service, there will be pickup requests originating from spatially distributed customer locations and these requests need be served by different depots with the help of a heterogeneous vehicle fleet. The D-MDCVRP is a dynamic version of MDCVRP, which requires dynamic rerouting of vehicles to respond to real-time customer requests.
The Quantum computing paradigm uses the principles of quantum mechanical systems such as superposition and entanglement to improve the algorithmic computational performance when compared to classical counterparts. Adiabatic quantum computation, introduced by Farhi [3], is based on the the adiabatic theorem [4]. Quantum annealing works on the principle of adiabatic quantum computation. Kadowaki-Nishimori [5] proposed the quantum annealing technique by introducing quantum fluctuations into the classical simulated annealing algorithm for solving combinatorial optimization problems in the traverse Ising model, which is equivalent to a quadratic unconstrained binary optimization (QUBO). Therefore, quantum annealing can be used to solve optimization problems that follow the QUBO formulation.
Quantum annealing was used to solve a wide variety of combinatorial optimization problems in fields such as bioinformatics [6], Quantum chemistry [7], computational biology [8], traffic flow optimization [9], fault diagnosis [10], train-ing of deep neural networks [11], vehicle routing [12], job shop scheduling [13] and nurse scheduling in healthcare [14]. Following the application of quantum annealing for solving combinatorial optimization in various domains, we provide the formulations for solving the MDCVRP and D-MDCVRP in a way that facilitates quantum annealing optimization.
Paper Contributions: The contributions made through this paper are: (1) a quadratic unconstrained binary optimization (QUBO) formulation for a MDCVRP under vehicle/depot capacity constraints; (2) a QUBO formulation for dynamic rerouting of multiple vehicles in response to real-time customer requests; and (3) discussion on problem complexity and a solution framework through quantum annealing.
Paper Organization: Section II provides a background to quantum annealing, QUBO, and D-Wave hardware for quantum annealing. Section III provides a review of existing literature on the use of quantum annealing for solving vehicle routing problems. Section IV provides the MDCVRP formulation with constraints along with associated QUBO formulation. Section V discusses the D-MDCVRP formulation. Section VI discusses on the problem complexity and solution framework for solving MDCVRP and D-MDCVRP followed by concluding remarks and future work in Section VII.

II. BACKGROUND
Quantum annealing is a metaheuristic algorithm for solving combinatorial optimization problems based on the principles of quantum mechanics [5]. Quantum annealing is implemented by evolving the Hamiltonian dynamics from an initial quantum state to a final quantum state, which corresponds to the solution of an optimization problem of interest [15]. In this context, the dynamism is strictly adiabatic and thus relates to the concept of adiabatic quantum computing [16]. If H I and H F represent the initial and final Hamiltonians of a quantum system, then the Hamiltonian at any time, H(t), between the evolution can be written as [16]: In Eq. 1, T is the time taken to evolve the system from the initial state to the final state, and t is any time between 0 and T (0 ≤ t ≤ T ). The quantum system is initially in the ground state of the initial Hamiltonian. When the system is evolved slowly in an adiabatic manner, the system stays in the ground state. The final Hamiltonian, and thus the final ground state, corresponds to the minimization problem of interest.
A quantum system comprises of several individual qubits, and the Hamiltonian of an N-qubit Ising system can be written as: where s i represents the state of i th qubit, which can be either −1 or 1, h i and J ij are the bias and interaction terms relating to individual qubits [17], and s is a vector of all qubit states.
In quadratic unconstrained binary optimization (QUBO), the value taken by each individual decision variable is either 0 or 1. The Ising model can be translated to a binary optimization representation using the following transformation: x i = 1+si 2 by which the −1 and 1 states of s i are mapped to 0 and 1 states of x i . Therefore, Eq. 2 can be written as: In Eq. 3, Q is an N × N symmetric matrix. Therefore, any optimization problem that can be represented as a QUBO can be solved using quantum annealing. If a constrained binary optimization is available, then it needs to be converted to a QUBO. A commonly used method for such a conversion is the penalty method [18]. For example, if f (x) and c(x) ≤ 0 represent the objective function and constraint respectively, then an equivalent unconstrained formulation can be written as q(x) = f (x) + γg(c(x)), where γ is the penalty term, and g is a function defined over the constraint. In this paper, we discuss QUBO formulations of MDCVRP and D-MDCVRP so they can be solved using quantum annealing.

III. LITERATURE REVIEW
As this paper is considering vehicle routing-related problems, we provide a brief review of previous work that considered quantum annealing approaches to solve such problems.
Boros mentioned a set of search algorithms for solving QUBO problems by illustrating simulation results obtained from various computational experiments [19]. According to Choi [20], the adiabatic quantum computation technique can solve QUBO problems that use an Ising-spin Hamiltonian. Furthermore, the adiabatic quantum computation can also solve constrained polynomial optimization problems [21]. Vyskocil [22] proposed how to solve a constrained mixedinteger linear programming problem in the QUBO framework by eliminating large coefficients that often result due to quadratic penalties. The computational techniques such as the adiabatic quantum method, quantum circuits, and quantum walks are suitable to solve Hamiltonian problems but, to address the QUBO problem, we require the application of Dwave architecture platform [23].
Moylette [14] proposed a quadratic speedup quantum algorithm for TSP by using the quantum backtracking algorithm to a classical computation algorithm. The above mentioned previous works show that quantum techniques prove to be beneficial for solving VRP related problems such as MDVRP that motivated towards this work.
Hirotaka [24] proposed a QUBO formulation for capacitated vehicle routing problem (CVRP), by introducing the concept of time, capacity and state of vehicles associated with every departure and destinations locations. Sebastian [25] formulated a hybrid method for solving CVRP using quantum annealer. The heuristic based approach used clustering and routing phase to determine the efficient vehicle routes in each cluster. Clark [26] investigated providing real-time routing for multiple robots, using hybrid classical-quantum approach for generating collision free routing for multiple robots on simulation grid. Christos [27] developed a QUBO model formulation for travelling salesman problem (TSP), related to time windows that can handle a small scale TSP with time windows on D-Wave platform.
From the above literature review, we noticed that previous works considered routing problems with a single entity (such as a vehicle or a traveling salesman); however, in real-world applications, we commonly encounter routing problems with heterogeneous vehicles and depots. Therefore, this paper focuses on QUBO formulations for solving MDCVRP and D-MDCVRP on the quantum annealing platform.

A. Problem Parameters
Here, we define the parameters that will later be used in the optimization formulation.
x ijk : A binary variable, which is equal to 1 if location j is served after location i by vehicle k.
µ ik : A binary variable, which is equal to 1 if location i is the first location served by vehicle k after leaving its depot.
η ik : A binary variable, which is equal to 1 if location i is the last location served by vehicle k before returning to its depot.
γ kd : A binary variable, which is equal to 1 if vehicle k belongs to depot d.   x jik + k∈K µ ik = 1 ∀i ∈ T Constraint 3: Each trip should start and end at a depot. The variable µ ik = 1 indicates that the first customer served by vehicle k after starting from the depot is i. Similarly, η ik = 1 indicates that the last customer served by vehicle k before returning to the depot is i. Each vehicle starting from a depot should visit only one customer location. i∈T µ ik = 1 ∀k ∈ K Constraint 4: Each vehicle should end at its depot from only a single customer location.
Constraint 5: When a vehicle k reaches a location i from j, it should leave i to reach some other location p. This ensures that a vehicles route is continuous and does not terminate at a customer location. This is often referred to as a flow constraint or continuity constraint.
Constraint 6: A vehicle k should not form closed loops with a subset of customer locations. A closed loop (also called a subtour) is formed when the numbers of arcs traveled between any subset of customer locations across all vehicles is equal to the number of customer locations in that subset. Therefore, we constrain the maximum number of arcs to be one less than the number of customer locations. This is known as the subtour elimination constraint [28].
For example, consider a scenario with two depots D 1 and D 2 , and one vehicle at each depot. Let vehicle V 1 be associated with D 1 and V 2 with D 2 . Assume that there are four customer locations, C i , i = 1 . . . 4 to be served. When the subtour elimination constraint is not provided, a potential set of routes for vehicles can be as follows. For V 1 , the route can be Here, the vehicles satisfy all the five constraints above, but provide an infeasible route where V 1 goes in a loop between C 2 and C 3 . To eliminate such loops, we provide the subtour elimination. This constraint needs to be provided for every subset of customer locations.
Let P(T ) denote the power set of T , i.e., the set of all subsets derived from T . Let S represent an element from P(T ), whose cardinality is at least two since formation of a subtour requires atleast two locations. The subtour elimination constraint can be written as Constraint 8: Total customer demand across all vehicles should be less than the depot capacity (assuming that the depot capacity is less than the summation of all the vehicle capacities). k∈K γ kd i∈T j∈T j =i

C. Problem Hamiltonian
Here, we first detail the general approach for the construction of Hamiltonian for any equality and inequality constraint. We later use this approach to construct the Hamiltonian for various constraints in the MDCVRP formulation. First, let us consider an equality constraint of the form, where x i is the i th binary decision variable, A i represents its coefficient, b is an integer constant, and n x is the number of decision variables. The Hamiltonan term corresponding to this equality constraint can be written as ( Now, let us consider an inequality constraint of the form, To represent any optimization formulation as a QUBO, the inequality constraints need to be transformed into equality constraints [22]; this is accomplished by introducing an additional set of binary decision variables called slack variables, and the corresponding Hamiltonian can be written as ( Here, n λ represents the number of slack variables, which can be calculated as n λ = 1 + log 2 b , where . is the ceiling function. The number of slack variables should be sufficient enough to represent all the values from 0 (when all the x i variables are zero) to b. Following this discussion, we will write the Hamiltonian terms of various equality and inequality terms in the MDCVRP formulation. We will first start with writing the Hamiltonian term corresponding to the objective function.
The Hamiltonian term, H O that corresponds to the objective function is simply the objective function in Eq. 4.
The Hamiltonian terms that correspond to various constraints in Section IV are given below. Let C i , i = 1 . . . 8 represent the eight constraints, and let H Ci , i = 1 . . . 8 represent the Hamiltonian terms corresponding to the constraints.
In Eqs. 14-21, B is a large positive constant that corresponds to the penalty incurred when the constraints are violated. After obtaining Hamiltonian terms that correspond to the objective function and various constraints, the overall Hamiltonian is equal to the sum of the individual Hamiltonian terms. Therefore, After static formulation, let us now discuss formulation to facilitate dynamic rerouting to respond to real-time customer requests.

REROUTING
The vehicles across several depots will initially be routed based on the output from the static MDCVRP formulation discussed in Section IV. Assume that at a certain time when are vehicles are in service, a new set of customer requests denoted as R become available. Therefore, the overall set of customer locations, denoted as W = T ∪ R. Let Γ W and Θ W represent the set of locations that were served and yet to be served respectively. Note that Γ W ⊆ T , R ⊆ Θ W , and W = Γ W + Θ W . When a new set of customer requests arrives and rerouting is performed, the positions of vehicles can be anywhere along the obtained from the static formulation.
Let C represent the set of current locations of the vehicles. Let γ kc represent the binary variable that determines the current location of vehicle k. γ kc = 1 when vehicle k is at location c and 0 otherwise. These current locations will now become the initial locations of the vehicles for the rerouting process. In the static formulation, the initial locations of various vehicles are their associated depots. Since the vehicles have limited capacity, and some customer locations were already served, the remaining vehicle capacities and similarly, remaining depot capacities needs to be considered in the rerouting analysis. We provide the revised formulation that needs to be adopted for dynamic rerouting.
The objective function, which is the minimization of total distance to be covered by all the vehicles across all depots can be written as: Constraints 1-6 (Eqs. 5-10) will have the same form as in the static formulation except for a change in the values taken by various indices. The revised set of constraints are given below.
k∈K j∈ΘW j =i k∈K j∈ΘW j =i Here, β ik corresponds to the first location reached by a vehicle from the initial position in the rerouting process, and β ik = 1 if location i is the first location served by vehicle k.
k∈K i,j∈S i =j i∈ΘW j∈ΘW j =i In Eq. 29, γ kw = 1 if vehicle k served a location w in Γ W . k∈K γ kd i∈ΘW j∈ΘW j =i The procedure for writing the Hamiltonian terms corresponding to various objective and constraints in the dynamic rerouting formulation is similar to that used for the static formulation in Section IV-C. Since the objective function (Eq. 22) and the first six constraints (Eqs. [23][24][25][26][27][28] are very similar to the MDCVRP formulation, we will write the Hamiltonian terms for the last two constraints denoted as H d C7 and H d C8 . Following Section IV-C, the number of slack variables required for a given k and d in Eq. 29 and 30 can be calculated as n λk = 1 + log 2 (Q k − w∈ΓW q w γ kw ) and n λd = 1+log 2 (V d − k∈K w∈ΓW q w γ kw γ kd ) respectively. Given n λk and n λd , H d C7 and H d C8 can be written as

VI. DISCUSSION
In this section, we discuss the problem complexity and a solution framework for solving the static and dynamic multidepot capacitated vehicle routing problem.
Problem Complexity: We will use the MDCVRP formulation discussed in Section IV to illustrate the problem complexity. However, the same discussion applied to the dynamic formulation as well. The total number of decision variables is equal to the number of binary decisions that are used to describe various vehicle routes (e.g., x ijk , µ ik , η ik ) and also the number of slack variables (e.g., λ lS , λ lk , λ ld ) introduced to convert the inequality constraints to equality constraints. For simplicity, let us refer to the former set of decision variables as route decision variables. The total number of route decision variables is equal to |T |(|T |−1)|K|+2|T ||K|, where |T |(|T | − 1)|K| is the number of x ijk variables and |T ||K| represents the number of each of µ ik and η ik .
With regard to slack variables, let us first consider the number of slack variables in the subtour elimination constraint (Eq. 10). We will have a subtour elimination constraint for every |S| combinations of customer locations. Let n z represent the number of combinations of obtaining z elements from a set of n elements. Therefore, the total number of subtour elimination constraints will be equal to 2≤|S|≤|T | |T | |S| . For a given value of |S| in the subtour elimination constraint, the number of slack variables is equal to 1+log 2 (|S|−1) . The total number of slack variables considering all the subtour elimination constraints is equal to |T | |S|=2 |T | |S| 1 + log 2 (|S| − 1) . Following the discussion on subtour elimination constraint, the number of slack variables in the vehicle capacity constraint (Eq. 11) and depot capacity constraint (Eq. 12) can be calculated as k∈K 1 + log 2 Q k and d∈D 1 + log 2 V d respectively. Thus, the total number of decision variables (N D ) is equal to (33) Solution framework: In order to solve the problem using quantum annealing, the overall Hamiltonian needs to be written in the QUBO form. This can be accomplished using the PyQUBO package in Python [29]. The PyQUBO compiles the overall Hamiltonian and obtains the coefficient matrix of the binary decision variables, i.e., the Q matrix (in Eq. 3). Depending on the total number of decision variables, the QUBO formulation can be solved either using a quantum solver (such as D-Wave 2000Q) or a hybrid quantum-classical solver [30]. A hybrid solver decomposes the problem into several smaller problems and use a combination of classical and quantum solvers to improve the computational performance. Since D-Wave 2000Q has 2048 qubits, a hybrid solver can handle a higher number of decision variables as it uses a combination of classical and quantum solvers. Also, according to Chancellor [31], the Chimera graph is not completely utilized due to physical constraints.The analysis steps are summarized below: 1) Choose penalty values associated with various constraints; the penalty values may be chosen to be equal for all constraints. 2) Construct the overall Hamiltonian for the objective and constraint functions by converting any inequality constraints to equality constraints through addition of slack variables. 3) Obtain the matrix of coefficients (the Q matrix) of the decision variables using the PyQUBO package. 4) Solve the QUBO problem by choosing an appropriate solver, either a quantum or a hybrid solver, depending on the number of decision variables.

5) Go to
Step 1 if the obtained solution is not desirable.

VII. CONCLUSION
This paper discussed quadratic unconstrained binary optimization (QUBO) formulations of Multi-Depot Vehicle Routing Problem (MDVRP) and its variant, the Dynamic MDCVRP (D-MDCVRP) so they can be solved on the quantum annealing hardware such as D-Wave 2000Q. MDCVRP and D-MDCVRP seek to serve a set of spatially distributed customer requests through a heterogeneous vehicle fleet operated from multiple depots with heterogeneous capacities. The MDCVRP and D-MDCVRP formulations contain both equality and inequality constraints. To derive the QUBO formulations, the inequality constraints need to be converted to equality constraints. In this paper, we introduced an additional of decision variables called the slack variables to convert inequality constraints to equality constraints. We have also discussed the problem complexity in terms of the number of total decision variables and also a step-by-step solution framework to solve the associated QUBO formulations on the quantum annealing hardware. The methods that we discussed in this paper are applicable to combinatorial optimization problems in other domains such as resource management in data centers, warehouses, and microgrids.
As part of our future work, we will solve the QUBO formulations on the D-Wave 2000Q quantum annealer, and compare the accuracy and computational time of the results against classical heuristic algorithms such as Tabu search. We will also investigate the scalability of the quantum annealing approach by comparing the variation of computational time with respect to the number of vehicles in each depot, and number of depots.