Quantum Annealing of Vehicle Routing Problem with Time, State and Capacity

TL;DR

This paper introduces a novel QUBO formulation for the capacitated vehicle routing problem that incorporates time, capacity, and state constraints, enabling quantum annealing solutions with D-Wave 2000Q.

Contribution

It presents a new QUBO model for CVRP that includes time-evolution, capacity adjustments, and vehicle states, expanding the scope of quantum annealing applications.

Findings

Successful formulation of CVRP as QUBO with time and capacity constraints

Demonstration of quantum annealing on D-Wave 2000Q for the problem

Enhanced modeling of vehicle states and traveling rules

Abstract

We propose a brand-new formulation of capacitated vehicle routing problem (CVRP) as quadratic unconstrained binary optimization (QUBO). The formulated CVRP is equipped with time-table which describes time-evolution of each vehicle. Therefore, various constraints associated with time are successfully realized. With a similar method, constraints of capacities are also introduced, where capacitated quantities are allowed to increase and decrease according to the cities which vehicles arrive. As a bonus of capacity-qubits, one also obtains a description of state, which allows us to set a variety of traveling rules, depending on each state of vehicles. As a consistency check, the proposed QUBO formulation is also evaluated by quantum annealing with D-Wave 2000Q.

Equations

$$\displaystyle x_{\tau,a}^{(i)}\qquad\bigl{(}1\leq\tau\leq T,\quad 1\leq a\leq N,\quad 1\leq i\leq k\bigr{)},$$

$$\displaystyle(\tau)\qquad\mapsto\qquad\bigl{[}t_{\tau},t_{\tau+1}\bigr{)}\qquad\bigl{(}1\leq\tau\leq T-1\bigr{)},\qquad t_{\tau+1}-t_{\tau}\equiv\Delta t.$$

$$\displaystyle x_{\tau,a}^{(i)}=1\quad\bigl{(}\text{arriving}\bigr{)},\qquad x_{\tau,a}^{(i)}=0\quad\bigl{(}\text{not arriving}\bigr{)}.$$

$$n_{ab}^{(\tau )}$$

$$d_{ab}^{(\tau )}$$

$$\mathcal{H}$$

$$\displaystyle+\underbrace{\lambda\times\biggl{(}\sum_{\left\{\begin{subarray}{c}1\leq a<b\leq N\cr 1\leq\tau\leq T\cr 1\leq i\leq k\end{subarray}\right\}}x_{\tau,a}^{(i)}\,x_{\tau,b}^{(i)}+\sum_{\left\{\begin{subarray}{c}1\leq a\leq N\cr 1\leq\tau\neq\tau^{\prime}\leq T\cr 1\leq i\leq j\leq k\end{subarray}\right\}}x_{\tau^{\prime},a}^{(i)}\,x_{\tau,a}^{(j)}+\sum_{\left\{\begin{subarray}{c}1\leq a\leq N\cr 1\leq\tau\leq T\cr 1\leq i<j\leq k\end{subarray}\right\}}x_{\tau,a}^{(i)}\,x_{\tau,a}^{(j)}}_{(\ast 3)}\biggr{)}.$$

$$\displaystyle\lambda\times x_{\tau,a}^{(i)}\,x_{\tau,b}^{(i)}\qquad\bigl{(}\forall a\neq\forall b,\forall i,\forall\tau\bigr{)}.$$

$$\displaystyle\lambda\times x_{\tau,a}^{(i)}\,x_{\tau^{\prime},a}^{(i)}\qquad\bigl{(}\forall a,\forall i,\forall\tau\neq\forall\tau^{\prime}\bigr{)}.$$

$$\displaystyle\lambda\times x_{\tau,a}^{(i)}\,x_{\tau^{\prime},a}^{(j)}\qquad\bigl{(}\forall a,\forall i\neq\forall j,\forall\tau,\forall\tau^{\prime}\bigr{)}.$$

$$\mu =d_{\max },\rho =\frac{d_{\max }-d_{\min }}{\lambda }.$$

$$x_{1,a}^{(i)}={\delta }_{a,s_i}$$

$$x_{\tau ,a}^{(i)}=0\text{when τ+nei,a(τ)>T},$$

$$\displaystyle\biggl{(}\frac{d_{ab}^{(\tau)}-\mu}{\rho}\times x_{\tau+n_{ab}^{(\tau)},a}^{(i)}\,x_{\tau,b}^{(i)}+\sum_{1\leq\delta\tau\leq n_{ab}^{(\tau)}-1}\lambda\times x_{\tau+\delta\tau,a}^{(i)}\,x_{\tau,b}^{(i)}\biggr{)}\bigg{|}_{b=e_{i}}$$

$$\displaystyle\quad\underset{\text{replace}}{\longrightarrow}\quad\biggl{(}\sum_{1\leq\delta\tau\leq T-\tau}\lambda\times x_{\tau+\delta\tau,a}^{(i)}x_{\tau,e_{i}}^{(i)}\biggr{)},$$

$$\displaystyle x_{\tau,a|c_{1},c_{2},\cdots,c_{M}}^{(i)}\qquad\bigl{(}1\leq\tau\leq T,\quad 1\leq a\leq N,\quad 1\leq i\leq k\bigr{)},$$

$$\displaystyle q_{m}\leq c_{m}\leq Q_{m}\qquad\bigl{(}q_{m},Q_{m}\in\mathbb{Z},\quad 1\leq m\leq M\bigr{)}.$$

$$\displaystyle q\leq c\leq Q\qquad\Leftrightarrow\qquad q_{m}\leq c_{m}\leq Q_{m}\quad\bigl{(}1\leq m\leq M\bigr{)}.$$

$$\displaystyle x_{\tau,a|c}^{(i)}=1\qquad\bigl{(}\text{arriving with the status $c$ of capacitated variables}\bigr{)}.$$

$$\displaystyle B_{ab}^{(\tau)}\equiv\bigl{(}B_{ab|1}^{(\tau)},B_{ab|2}^{(\tau)},\cdots,B_{ab|M}^{(\tau)}\bigr{)}.$$

$$B_{ab|m}^{(\tau )}=$$

$$\mathcal{H}=$$

$$+$$

$$c\mapsto c^{\prime }=c+B_{ab}^{(\tau )},$$

$$\displaystyle\lambda\times x_{\tau,a|c^{\prime}}^{(i)}\,x_{\tau,a|c}^{(i)}\qquad\bigl{(}\forall a,\forall i,\forall\tau,\forall c^{\prime}\neq\forall c\bigr{)}.$$

$${\hat{x}}_{\tau ,a}^{(i)}:$$

$$1)x\to \hat{x}:$$

$$n_a^{(\tau )}$$

$$\mathcal{H}$$

$$\displaystyle+\sum_{\left\{\begin{subarray}{c}1\leq a,b\leq N\cr 1\leq\tau\leq T-1\cr 1\leq i\leq k\end{subarray}\right\}}^{\text{``staying''}}\left(\begin{array}[]{c}\displaystyle\vspace{0.1cm}\delta_{ab}\biggl{(}\frac{0-\mu}{\rho}\times x_{\tau+n_{a}^{(\tau)},a}\widehat{x}_{\tau,a}^{(i)}+\sum_{1\leq\delta\tau\leq n_{a}^{(\tau)}-1}\lambda\times x_{\tau+\delta\tau,a}^{(i)}\widehat{x}_{\tau,a}^{(i)}\biggr{)}\cr\displaystyle+(1-\delta_{ab})\sum_{1\leq\delta\tau\leq n_{b}^{(\tau)}}\lambda\times\widehat{x}_{\tau+\delta\tau,a}^{(i)}\widehat{x}_{\tau,b}^{(i)}\end{array}\right)+$$

$$\displaystyle+\lambda\times\left(\begin{array}[]{l}\displaystyle\vspace{0.1cm}\sum_{\left\{\begin{subarray}{c}1\leq a\leq N\cr 1\leq\tau\leq T\cr 1\leq i\leq k\end{subarray}\right\}}\widehat{x}_{\tau,a}^{(i)}\,x_{\tau,a}^{(i)}+\sum_{\left\{\begin{subarray}{c}1\leq a<b\leq N\cr 1\leq\tau\leq T\cr 1\leq i\leq k\end{subarray}\right\}}\Bigl{(}x_{\tau,a}^{(i)}\,x_{\tau,b}^{(i)}+\widehat{x}_{\tau,a}^{(i)}\,x_{\tau,b}^{(i)}\Bigr{)}+\cr\displaystyle\vspace{0.1cm}\qquad+\sum_{\left\{\begin{subarray}{c}1\leq a\leq N\cr 1\leq\tau\neq\tau^{\prime}\leq T\cr 1\leq i\leq j\leq k\end{subarray}\right\}}\Bigl{(}x_{\tau^{\prime},a}^{(i)}\,x_{\tau,a}^{(j)}+\widehat{x}_{\tau^{\prime},a}^{(i)}\,x_{\tau,a}^{(j)}\Bigr{)}+\cr\displaystyle\qquad+\sum_{\left\{\begin{subarray}{c}1\leq a\leq N\cr 1\leq\tau\leq T\cr 1\leq i<j\leq k\end{subarray}\right\}}\Bigl{(}x_{\tau,a}^{(i)}\,x_{\tau,a}^{(j)}+\widehat{x}_{\tau,a}^{(i)}\,x_{\tau,a}^{(j)}\Bigr{)}+\cr\qquad\qquad\qquad\qquad\qquad\qquad+\Bigl{(}\text{``hermitian conjugate''}\Bigr{)}\end{array}\right).$$

$$\displaystyle\mathcal{H}=\sum_{A,B=1}^{N}\underbrace{\bigl{(}C_{AB}-\mu\bigr{)}}_{\equiv\epsilon_{AB}}\times x_{A}\,x_{B}-\xi\sum_{A=1}^{N}x_{A}+\sum_{\left\langle A,B\right\rangle\in E_{\text{frb.}}}\lambda\times x_{A}\,x_{B},$$

$$\xi =\lambda (\text{TSP}),\xi \ge \frac{3}{2}\lambda (\text{VRP}),$$

Full text

1. 1 Introduction

2. 2 Time-table in TSP/VRP

3. 3 Multiple-capacitated TSP/VRP

4. 4 State of vehicles

5. 5 Practical constraints with time, state and capacity

6. 6 Validity of the formulation with D-Wave 2000Q

7. 6.1 Negatively-shifted energy method and choice of baselines

8. 7 Conclusion and discussion

1 Introduction

Vehicle Routing Problem (VRP) [1] is a basic mathematical problem related to optimization of planning, logistics, and transportation. Owing to the recent interests in quantum annealing machines, which were first studied theoretically by Kadowaki-Nishimori [2] and made available commercially by D-Wave Systems Inc. [3], the investigation of VRP as a quadratic unconstrained binary optimization (QUBO) has become very important, particularly in an attempt to achieve quantum-mechanical optimization of real-world problems encountered in our daily life concerning various mobility services.

VRP is a generalization of Traveling Salesman Problem (TSP), i.e., a problem to find the traveling path that has the lowest cost. Similarly, the purpose of VRP is to find the best routing scheme for multiple vehicles that achieves the lowest cost under various circumstances; for instance, each city should be visited exactly once. QUBO formulation of the TSP is found in [4]; it was constructed by adding the square of linear-constraint functions to the associated cost functions [5, 6, 7] using Lagrange-like multipliers. A straightforward extension of this QUBO formulation can also be applied to VRP by introducing several copies of QUBO systems for TSP. However, such formulation intrinsically suffers from the strict concept of time. It does have a concept of time-step, but it is generally not equivalent to the concept of time in many practical applications. Obviously, if one of the cities is located far from the rest, then traveling there will take more time than traveling to other cities. In such case, the conventional time-step formulation cannot describe the time, which flows commonly and homogeneously for all vehicles. Therefore, the introduction of time in the conventional QUBO formulations is the main obstacle in formulating various important VRP constraints associated with time such as time-window, (non-)simultaneous arrivals, and chronological variation of cost.

In addition to the concept of time, it is also important to introduce the concept of capacity (i.e., capacitated VRP or CVRP), which may describe the capacity of carrying passengers or packages. There are some attempts toward a QUBO formulation of CVRP: one tackles a possibility of constraint terms describing inequality [8]; and another utilizes a hybrid cluster algorithm combined with the TSP QUBO systems [9]. In this paper, we propose a different approach to investigate a concise QUBO formulation. In fact, capacity and time are similar as the description of time already implies that vehicles should travel within their own capacity of time. As further explained, these two concepts are similarly implemented by introducing a new kind of interactions that depend on the cities of departure and destination only. Furthermore, by introducing capacity-qubits (See Section 3), one can realize multiple capacitated variables which can be allowed to increase and decrease (i.e. pickup and delivery) during each travel.

However, there are further applications of this approach. As a bonus of capacity-qubits, one can also describe the concept of state. In particular, cost and time-duration can change depending on each state of vehicles. As a simplest example, we demonstrate a two-state model: arrival-state and departure-state. The mean time from arrival to departure is the duration of visits. We can also set up how long each vehicle stays depending on the cities visited by vehicles.

The organization of this paper is as follows: The brand-new QUBO formulation of VRP is introduced with time-table (in Section 2), with capacity-qubits (in Section 3), and with the concept of state (in Section 4). Other related constraints from real-world applications are in Section 5, and validity of our models is discussed in Section 6. Conclusion and discussion are presented in Section 7.

2 Time-table in TSP/VRP

The first new ingredient introduced to the TSP/VRP concept in this paper is time. In particular, we introduce time-table in the TSP/VRP formulation. First, we would prepare binary qubits parametrized by three integers $(\tau,a,i)$:

[TABLE]

where $\tau$ parametrizes each time-interval of the time-table, with the assumption that there are $N$ cities to visit and $k$ vehicles are present. The time-interval means that we divide the total time into several time units as follows:

[TABLE]

Herein $\Delta t$ is the unit of time-division.111Note that the unit of time-division can also depend on time $(\tau)$ and vehicle $(i)$, i.e., $\Delta t_{\tau}^{(i)}\equiv t_{\tau+1}^{(i)}-t_{\tau}^{(i)}$. Although such a generalization is also important for some applications, we keep the uniform value $\Delta t_{\tau}^{(i)}=\Delta t$ for the sake of simplicity. For example, three hours from 9:00 AM to 12:00 AM can be divided into nine intervals with twenty-minute duration each. Suppose a vehicle $(i)$ does (or does not) arrive at a city $(a)$ in a time-interval $(\tau)$, the binary qubits $x_{\tau,a}^{(i)}$ takes the following value:

[TABLE]

Since the conventional QUBO formulation of TSP/VRP always assumes one-step forward, we would invent a new kind of interaction which describes how much time each vehicle spends for each travel. Thus, we would first introduce time-duration matrices \bigl{(}n_{ab}^{(\tau)}\bigr{)}_{1\leq a\neq b\leq N} as well as cost matrices \bigl{(}d_{ab}^{(\tau)}\bigr{)}_{1\leq a\neq b\leq N} for each time-interval $(\tau)$ as222Here $\lceil x\rceil$ ($x\in\mathbb{R}$) is the minimum integer which is greater than or equal to $x$, i.e., $x\leq\lceil x\rceil<x+1$ and $\lceil x\rceil\in\mathbb{Z}$.

[TABLE]

Thus, our proposed Hamiltonian $\mathcal{H}$ can then be written as follows:

[TABLE]

The new form of interaction is introduced in the first line: $(\ast 1)$ and $(\ast 2)$

The first term $(\ast 1)$ gives the cost of travel from a city $(b)$ (departed at $\tau$) to a city $(a)$ (the arrival is at $\tau^{\prime}=\tau+n_{ab}^{(\tau)}$).

• 2.

The second term $(\ast 2)$ forbids any early arrivals at a city $(a)$ (i.e., at $\tau^{\prime}$ in range $1\leq\tau^{\prime}-\tau<n_{ab}^{(\tau)}$). This is realized by repulsive interactions set forward from the departure city $(b)$ to the arrival city $(a)$.

The other terms $(\ast 3)$ are obtained from the basic constraints of penalty terms:

Any vehicle $(i)$ in any time-interval $(\tau)$ will not arrive at two different cities $(a)$ and $(b)$ simultaneously:

[TABLE]

• 2.

If a vehicle $(i)$ arrives at a city $(a)$ in a time-interval $(\tau)$, then that vehicle $(i)$ will not arrive at a city $(a)$ in any other time-interval $(\tau^{\prime})$:

[TABLE]

• 3.

If a vehicle $(i)$ arrives at a city $(a)$ in a time-interval $(\tau)$, then the other vehicles $(j)$ will not arrive at a city $(a)$ in any time-interval $(\tau^{\prime})$:

[TABLE]

At this point, we no longer employ the conventional method of constraint functions (i.e., square of linear-constraint functions, applied in [4]). Instead, we have introduced an additional parameter $\mu\,(>0)$ around the traveling cost $d_{ab}^{(\tau)}$, as well as the standard parameters $\rho\,(>0)$ and $\lambda\,(>0)$.

• •

Overall negative shift of costs (delivered by $\mu$) replaces the role of the negative linear terms generated in the conventional method of constraint functions (See Section 6).

• •

Up to the overall scaling, at least two parameters ($\mu,\rho$) should be adjusted properly to optimize the performance of each Ising machine. The remaining $\lambda$ is adjusted by total scaling to the maximum value suited for each Ising machine.

• •

However, as discussed in Section 6, the standard value of the parameters $\mu$ and $\rho$ are selected as follows:

[TABLE]

Generally, the initial starting points of the vehicles can be implemented using boundary condition of the binary qubits. For instance, one can choose

[TABLE]

Conversely, to set the final destination, we turn off several binary qubits from which one cannot reach the final destination before exceeding the upper time limit $T$ (where the final destination of vehicle $(i)$ is the city $(e_{i})$):

[TABLE]

Additionally, the forward interactions from the final destination should be replaced as follows:

[TABLE]

to make sure that any travel after arriving at the final destination is forbidden.

As this TSP/VRP formulation can deal with time-scheduling, we refer it to as time-scheduled TSP/VRP or TS-TSP/VRP.

3 Multiple-capacitated TSP/VRP

In the last section, we have introduced the concept of time and, as is mentioned in Introduction, time is a kind of capacitated variable whose consumption is accumulated until full of its “capacity”. Therefore, one can replace the concept of time by the concept of capacity for “a monotonically increasing/decreasing capacitated variable”. This replacement is also useful in some practical applications, especially in a case of saving the number of qubits. However, such a concept of time/capacity only resolves the problem of single-capacitated TSP/VRP.

In this section, we shall introduce the concept of multiple capacities in addition to our time-scheduled TSP/VRP formulation, which is referred to as time-scheduled multiple-capacitated VRP or TS-mCVRP. Addition of multiple-capacity can be achieved by adding capacitated variables $(c_{1},c_{2},\cdots,c_{M})$ in the binary qubits of Eq. (2.1) as follows:

[TABLE]

with their capacity bounds given as

[TABLE]

In many cases, we use a shorter notation, $c=(c_{1},c_{2},\cdots,c_{M})$, as a vector of $M$-dimensional capacitated variables.333We also use the following “array” notation:

\displaystyle q\leq c\leq Q\qquad\Leftrightarrow\qquad q_{m}\leq c_{m}\leq Q_{m}\quad\bigl{(}1\leq m\leq M\bigr{)}.

(3.3)

The introduction of capacitated variables is understood as the direct product of capacity-qubits and the original VRP-qubits (See Fig. 1). The interpretation of the binary values, Eq. (2.3), should be now incorporated with the status of capacitated variables as follows:

[TABLE]

Furthermore, to describe variation of capacitated variables, we introduce variation matrices \bigl{(}B_{ab|m}^{(\tau)}\bigr{)}_{1\leq a\neq b\leq N} for each time-interval $(\tau)$:444It is also convenient to use a vector notation of the collection of capacity-variation matrices:

\displaystyle B_{ab}^{(\tau)}\equiv\bigl{(}B_{ab|1}^{(\tau)},B_{ab|2}^{(\tau)},\cdots,B_{ab|M}^{(\tau)}\bigr{)}.

(3.5)

[TABLE]

We can write the Hamiltonian as follows:

[TABLE]

As is in Eq. (2.6), the first line includes the interaction associated with consumption of time and variation of capacitated variables. In particular,

• •

The second term $(\ast 4)$ in the first line ensures that transmission of capacitated variables in each travel satisfies the expected variation relation,

[TABLE]

• •

The third term $(\ast 5)$ forbids early arrivals at a city $(a)$ from a city $(b)$ with any changes of capacitated variables, $c\to c^{\prime}$.

The second part represents the basic constraints of penalty terms:

The newly introduced term is the first term $(\ast 6)$ and it represents that a vehicle $(i)$ arriving at a city $(a)$ in a time-interval $(\tau)$ cannot be assigned more than one capacitated-variable status, i.e., $c^{\prime}\neq c$.

[TABLE]

• 2.

The other terms are essentially the same as in the previous section.

As is in the previous section as well as Section 6, the parameters $(\mu,\rho,\lambda)$ should be chosen as the same standard values as stated in Eq. (2.10).

A new feature obtained by our capacity formulation is that the variation matrices $B_{ab|m}^{(\tau)}$ can take any integer number. Therefore, the capacitated variables can increase and/or decrease (i.e., pickup and/or delivery), strictly satisfying the capacity bounds.

4 State of vehicles

In addition to the simple concept of capacity, the capacitated variables can also be interpreted as the states of vehicles. Through the use of “state”, we can introduce various traveling rules depending on each state of vehicles. For the sake of simplicity, we shall demonstrate two-state TSP/VRP here.

As a simple example of a two-state TSP/VRP, we shall consider the arrival-state and departure-state of the vehicles. In terms of qubits, the states of vehicle are now denoted by the hat “$\widehat{\ \ }$” above the binary variables:

[TABLE]

A new ingredient of introducing two states is that we can now describe traveling/staying phases of vehicle:

[TABLE]

and other transitions, say $x\to x$ and $\widehat{x}\to\widehat{x}$, should be forbidden. As an additional phase of motion is introduced, we can further add a new time-duration matrix $\widehat{n}_{a}^{(\tau)}$ to describe how long the vehicle stays at the city $(a)$:

[TABLE]

Therefore, various traveling rules are additionally introduced as follows:555“Hermitian conjugate” of the variables is defined as: $\widehat{x_{\tau,a}^{(i)}}=\widehat{x}_{\tau,a}^{(i)}$ and $\widehat{\widehat{x}_{\tau,a}^{(i)}}=x_{\tau,a}^{(i)}$.

[TABLE]

Construction of the Hamiltonian is essentially the same as that of capacity (particularly the last line, which was obtained by re-interpreting capacitated variables $c$ of Section 3 as states). A distinguishing point is that we can introduce different cost matrices $d_{ab}^{(\tau)}$ and duration matrices $n_{ab}^{(\tau)}$ depending on the state (i.e., phase) of each vehicle.

This system (including “states” of vehicle) is referred to as time-scheduled state-vehicle routing problem or TS-SVRP. As a trivial extension, one can also consider multiple-state models and also the full system of time-scheduled multiple-capacitated state-vehicle routing problem or TS-mCSVRP.

5 Practical constraints with time, state and capacity

In practical applications for mobility service, there are various kinds of constraints that should be simultaneously implemented. The following list describes examples of such constraints, accepted by our TS-mCSVRP formulation:

setting: Consider an optimization problem for a delivery service. Cost of delivery is the total delivery time, and there are 50 customers, which are shared by five vehicles.

• –

Five vehicles comprises three middle-size cars and two trucks.

• –

Delivery is served in eight hours (e.g., from 9:00 to 17:00).

• –

Delivery schedule is set by a twenty-minute unit.

• –

The cost of optimization would be the total consumption of delivery time.

• 2)

time-variation:  Traffic conditions can change depending on the delivery time. Therefore, the cost $d_{ab}^{(\tau)}$ depends on time $\tau$. Similarly, the cost of inbound and outbound routes is different in general, i.e., $d_{ab}^{(\tau)}\neq d_{ba}^{(\tau)}$.

• 3)

priority for delivery:  Some customers may request a priority for delivery, which can be achieved by adding some extra weight factor to the cost function, $d_{ab}^{(\tau)}\to\varphi*d_{ab}^{(\tau)}$  ($0<\varphi<1$).

• 4)

time-window and type-window:  There are three kinds of constraints associated with window. These constraints can be realized by turning off the associated binary qubits (i.e., $x_{\tau,a|c}^{(i)}=0$) by hand:

• –

Each customer has a request of delivery time, scheduled by twenty minutes. The delivery-time request can be multiple (i.e., disconnected) time-window for each customer.

• –

Some of the vehicles cannot deliver to some customers; e.g. the roads are too narrow for trucks and some of the packages are too large or heavy for small cars to deliver. This requires type-window for delivery service.

• –

Each vehicle should serve within their own working hours. In particular, some of the drivers work in a short time. This requires time-window for vehicles.

Notably, only the formulation of time-window has a discrepancy with our capacity and state description as it induces overtime traveling sometimes (See Section 6). Except for time-window, one can further implement the following two constraints:

capacity-constraints:  Capacity constraints are realized using capacity-qubits:

• –

Each vehicle has its own volume and weight limitations for their capacity.

• –

For some vehicles that are used for both pickup and delivery, they need to be scheduled without over-capacity for volume and weight in delivering their service.

• 6)

scheduling as state:  Some detail about the schedule can be described using state-description.

• –

After arriving at each customer’s premise, drivers should spend twenty minutes on their customer services.

• –

Drivers can take a one-hour rest for lunch. This can be formulated using the time-dependence of duration matrices $n_{ab}^{(\tau)}$ or $n_{a}^{(\tau)}$.

6 Validity of the formulation with D-Wave 2000Q

Thus far, we have discussed the proposed new QUBO formulation of TSP/VRP. Conversely, it is also important to discuss the validity of the proposed formulation using a quantum annealing machine, D-Wave 2000Q.

Fig. 2 shows a primitive instance of a delivery service (TS-VRP with windows) for six customers with three vehicles in two hours. The unit of time-scheduling is chosen to be fifteen minutes. Interestingly, a travel from the city (4) to the city (5) automatically chooses overtime traveling, because time-window occasionally prevents the shortest travel.

Note that the parameters $(\rho,\mu,\lambda)$ are chosen as the standard value, Eq. (2.10), and therefore further optimization of the parameters would improve the performance. With taking into account this result, this instance shows that quantum annealing of our formulation works properly as far as small-size QUBO systems ($\simeq 83$ logical qubits for this instance) are considered.

6.1 Negatively-shifted energy method and choice of baselines

Further clarification is needed about our proposed method for the basic constraints (discussed around Eq. (2.10)). As mentioned earlier, we did not apply the conventional method of “square of linear-constraint functions,” as in [4]. Instead, we apply a negatively-shifted energy method; therefore we will discuss how our proposed method can replace the conventional method.

We consider a typical Hamiltonian, which can be generally expressed as follows:

[TABLE]

where $C_{AB}\,(>0)$ is the general cost matrix and $E_{\text{fbd}}$ is a set of forbidden pairs (i.e., edges) of configurations. We usually choose $\xi=0$.

In the conventional method (See [4]), the parameter $\xi$ inevitably exists and is strictly correlated with $\lambda\,(>0)$ as follows:

[TABLE]

and $\mu=0$. If one forgets about the correlation and chooses $\xi=0$ and $\lambda>0$ (and $\mu=0$ in the conventional method), the ground state becomes trivial, $x_{A}=0$ ($A=1,2,\cdots,N$). Therefore, the basic role of $\xi$-term is to enhance spontaneous popping-up of qubits: $x_{A}\to 1$ ($A=1,2,\cdots,N$).

However, such effect can also be generated using the overall negative shift of cost energy obtained from $\mu$. This is possible since optimization only cares about the relative values of cost. Hence we can adjust $\mu$ such that

[TABLE]

for all the pairs $\left<A,B\right>$ focused.666This means that $\epsilon_{AB}\leq C_{\rm max}^{(\rm focused)}-\mu=0\leq C_{\rm max}-\mu$ for all the focused pairs $\left\langle A,B\right\rangle$. Therefore, the parameter $\mu$ plays a role of “cut-off scale” of cost energy. In the proposed method, we can set $\xi=0$ or utilize it for other purposes.

In a sequence of configurations generated using the negative energy, Eq. (6.3), there are forbidden configurations, which receive penalty caused by $\lambda$. As these configurations are forbidden, we impose that such configurations receive relatively positive energy. In our proposed VRP, this condition is given by

[TABLE]

for all pairs $\left<A,B\right>$. We put a baseline (like a coastline) as energy $=0$ to separate feasible configurations (inside sea) from forbidden configurations (on the continent). This is referred to as the baseline condition.

From these two conditions, we choose the maximum range of cost energy,

[TABLE]

as the standard values. The larger range of cost causes larger resolution in quantum annealing machines. The solution of these conditions is given by Eq. (2.10).

It is also possible to strategically choose higher baseline in Eq. (6.4) and/or to select some focused energy range of Eq. (6.3) as

[TABLE]

which may improve the performance. Despite increase in forbidden configurations, it can further improve the range of energy input. In contrast, too large $\delta(>0)$ reduces performance, because the probability of ground states evaporates into the forbidden configurations, which possess much larger entropy.777This phenomenon also occurs in the conventional method: if the parameter $\xi$ becomes too large, it induces a large number of forbidden configurations. Therefore, too large $\xi$ also reduces the performance of quantum annealing.

7 Conclusion and discussion

In this paper, we proposed a new QUBO formulation of CVRP with time, state, and capacity. Introduction of the strict concept of time allowed us to formulate various constraints associated with time, while the introduction of capacity-qubits allowed us to formulate pickup and delivery during each travel of the vehicles. Introduction of state allowed us to describe various traveling rules depending on the state of the vehicles. We evaluated the proposed QUBO formulation using a quantum annealing machine, D-Wave 2000Q and the results show that our formulation properly works in small-size QUBO systems (less than $90$ logical qubits $\simeq$ $6\sim 7$ customers), which can be directly embedded in the current D-Wave machines.

For real-world applications of this formulation, on the other hand, at least more than $2000$ logical qubits ($\simeq$ more than $30$ customers for 20-minute scheduling of a half day) are required. In this sense, it is also interesting to evaluate our proposed TS-mCSVRP on digital Ising machines. It is also important to develop an efficient quantum/classical hybrid algorithm for our TS-mCSVRP, which should hasten the practical usage of our formulation. This point shall be reported in the future communication.

Acknowledgment

This work is based on a joint project of DENSO Corporation, Toyota Tsusho Corporation and Toyota Tsusho Nexty Electronics (Thailand) Co., LTD. The authors greatly thank Shunsuke Takahashi, Masakazu Gomi and Toru Awashima for valuable discussions which make this work possible.