We consider a fleet of N _bus electrical buses and a set of N _pile piles. Assuming that the electricity consumption of all buses is known during their respective operating hours, the objective of the EV-bus charging scheduling problem is to find an optimal daily schedule for charging the buses by minimizing the cost of electricity while avoiding any breakdown. The EV-bus charging scheduling problem with minimization of price can thus be formulated as follows, with T = ( t 1 , t 2 , … , t k , … , t T ) being a discretization of time periods: min ∑ i , j , k p k ω j x i j k , x i j k = o i k z i j k , (2) where z _ijk is a binary variable indicates whether bus i is charging at pile j for the time period t _k while o _ik = 0 if bus i is in operation and o _ik = 1 otherwise, i.e., if bus i is available for charging at time t _k . ω _j is the power output of pile j, p _k the cost of electricity at time t _k . We consider a time period of 24 h, and need only to find the values z _ijk during non-operating hours, i.e., during charging windows. This reduces the number of time steps to L _time ≤ T. In the sequel, we denote by S _il (respectively, D _il ) the starting (respectively, ending) time of the l-th charging window of bus i. The constraints of the problem can be expressed as: ∀ i , k : ∑ j x i j k ≤ 1 , (3) ∀ k : ∑ i , j x i j k ≤ N pile , (4) ∀ i , l : ∑ j , k = 1 D i l ω j x i j k − ∑ k = 1 D i l c i k + 1 ≤ 1 , (5) ∀ i , l : ∑ j , k = 1 D i l ω j x i j k − ∑ k = 1 S i ( l + 1 ) c i k + 1 ≥ 0.3 , (6) ∀ i , l : ∑ j , k ∈ [ S i l , D i l ] x i j k ( 1 − x i j ( k + 1 ) ) ≤ 1 . (7) Constraint (3) ensures that any bus can only be charged at most one pile at a given time, while (4) means that the total number of buses being charged at the same time should not be greater than the total number of piles N _pile. (5) and (6) control the state-of-charge (SOC) of each bus i between 100 and 30% at all times by adding constraints at each time a bus leaves the charging station, with c _ik being the power consumption of bus i at time t _k . Note here that the power ω _j is also normalized to be consistent with SOC ∈ [0, 1]. Lastly (7) guarantees that a bus can only be charged continuously, and not intermittently, at most once during each charging window. Hence, the original EV-bus charging problem Eqs 2–7 is an integer programming problem with inequality constraints.

Due to hardware constraints, quantum annealing can only be directly sampled from an Ising Hamiltonian. In this section, we reformulate thus the integer programming problem to a quadratic unconstrained binary optimization (QUBO) problem, which by definition does not include any constraint or non-binary variables [31, 33]. For each equality constraint of the form C _i (x) = 0 in the optimization problem, the penalty method consists in augmenting the cost function with a term P i 2 ( x ) which penalizes the solutions that violate the constraint. For the inequality constraints, following [31], we first transform them into equality constraints with slack variables and then apply the penalty method. More specifically, constraints (3) and (7) are equivalent to ∑ j x i j k ≤ 1 ≡ ∑ j x i j k = s 1 , i k , (8) ∑ j , k ∈ [ S i l , D i l ] x i j k ( 1 − x i j ( k + 1 ) ) ≤ 1 ≡ ∑ j , k ∈ [ S i l , D i l ] x i j k ( 1 − x i j ( k + 1 ) ) = s 5 , i l , (9) with s _1,ik , s _5,il ∈ {0, 1}. This yields the corresponding penalty functions: P 1 2 ( x ) = ( ∑ j x i j k − s 1 , i k ) 2 , (10) P 5 2 ( x ) = ( ∑ j , k ∈ [ S i l , D i l ] x i j k ( 1 − x i j ( k + 1 ) ) − s 5 , i l ) 2 . (11) It should be noted that Eq. 11 is of degree 4 in x, i.e, it is not a quadratic form. Therefore, a first limitation of the penalty method is that the already quadratic constraint (7) cannot be included in the augmented cost function.

Next, handling constraints (4) to (6) requires to introduce a slack variable that is either an integer [constraint (4)] or a real number [constraints (5)–(6)]. We choose to encode such non-binary variables s with binary encoding, i.e., s = ( 1 / Γ ) ∑ i = 0 N − 1 2 i s ( i ) , Γ = ∑ i = 0 N − 1 2 i , where N is the number of bits and upper-indexed s ⁽ⁱ⁾ are used to denote the encoding binary variable of s. Other options, not explored in this article, for encoding the continuous variables include one-hot encoding and order encoding, as in [34]. Constraints (4)–(6) become P 2 2 ( x ) = ∑ i , j x i j k − s 2 , k 2 , (12) P 3 2 ( x ) = ∑ j , k = 1 D i l ω j x i j k − ∑ k = 1 D i l c i k + s 3 , i l 2 , (13) P 4 2 ( x ) = ∑ j , k = 1 D i l ω j x i j k − ∑ k = 1 S i ( l + 1 ) c i k − 0.7 s 4 , i l + 0.7 2 , (14) where s _2,jk is a binary-encoded integer slack variable in [0, N _pile] and s _3,il , s _4,il are binary-encoded continuous slack variables in [0, 1]. Furthermore, it is easy to see that the penalties P ₃ and P ₄ can be written into a single one: P 3,4 2 ( x ) = ∑ j , k = 1 D i l ω j x i j k − ∑ k = 1 D i ( l + 1 ) c i k − s 3,4 , i l + 0.7 2 , (15) with s 3,4 , i l ∈ [ 0,0.7 − ∑ k = D i l S i ( l + 1 ) c i k ] .

Finally, it can be observed that the barrier method in classical optimization is a more efficient way of handling inequality constraints [35], yet not applicable here due to the log function required by this method but not feasible in practice on the quantum annealer.

By employing the penalty method and encoding the continuous variables, the original constraints are reformulated into penalty functions P _{1 … 5}. The effective QUBO Hamiltonian reads then: H QUBO ( x ) = H ∗ ( x ) + ∑ i λ i P i 2 ( x ) , λ i ∈ R , (16) where H ^∗ = ∑ _i,j,k p _k ω _j x _ijk is the Hamiltonian corresponding to the original cost function (2). It remains to determine the values of λ _i such that the structure of the original problem can be kept unchanged in the reformulated QUBO problem. More specifically, let C denote the set of all viable solutions x, we need to ensure that: ∀ x 1 ∈ C , x 2 ∉ C : H QUBO ( x 1 ) < H QUBO ( x 2 ) , (17) ∀ x 1 , x 2 ∈ C  and  H ∗ ( x 1 ) < H ∗ ( x 2 ) : H QUBO ( x 1 ) < H QUBO ( x 2 ) . (18) We proceed to an analysis of the different terms in Eq. 16, following the method in [31]. The detailed calculations are described in the Supplementary Material, and we find that in order to preserve the structure of the original problem, it is sufficient that the coefficients λ _i satisfy the following inequalities: min δ p × min ω 2 − 2 N − 2 > λ 3,4 > max p × max ω × L time ε 2 , λ 1,2 > max p × max ω × L time , (19) where   min δp is the minimal price difference between any two given times during the day (including the difference between lowest price and 0), ɛ is a margin to control the state of charge within the boundaries [0.3 + ɛ, 1 −ɛ], typically set to ɛ = 0.1. The non-quadratic penalty P 5 2 is omitted. These inequalities indicate thus how to appropriately set the values of the coefficients λ _i for the penalty functions. Another practical issue is the ratio r between the largest and smallest terms in the QUBO objective function H _QUBO. Indeed, a large ratio will require the hardware device to have high dynamic range and precision, which is challenging for the current technology advancement. More precisely, recalling that the original objective function is H* = ∑ _i,j,k p _k ω _j x _ijk , we obtain the lower and upper bounds, respectively: r m a x > 1 min p × min ω × min δ p × min ω 2 − 2 N − 2 ∼ 1 2 − 2 N − 2 , r m i n < max p × max ω × L time ε 2 × 1 max p × max ω = L time ε 2 , (20) where r _max , r _min give the range of the ratio r for which the structure of the problem stays the same. Hence, in order to guarantee that the structure of the problem is preserved, N, ɛ, L _time should be such that: 1 2 − 2 N − 2 > L time ε 2 . (21) Furthermore, observe that the upper bound L time ε 2 of r _min scales with the number of time steps, L _time but is independent of N _bus and N _pile. Consequently, increasing the size of this problem solved by quantum annealing and penalty method would not represent a higher challenge in terms of dynamic range of the physical device. This is a satisfying result as the QUBO form of many hard problems requires scaling dynamic range [36]. On the contrary, increasing the number of time steps requires higher performances from the QPUs.

Although the quadratic penalty model maps the constrained optimization problem to an unconstrained one while preserving the relation between solutions, there are several practical limitations while implementing the quadratic penalty model on a real quantum annealer. Firstly, the Hamiltonian that can be realized onto the real physical device is an Ising Hamiltonian which contains at most two-body interactions, i.e., quadratic terms in the objective function. In the previous Section 2.3, we have shown that some penalty functions can contain higher-order terms that can not be written into the QUBO objective function, cf. quadratic constraint (7) and the corresponding penalty (11).

The second practical limitation of the penalty method relates to the connectivity of the hardware architecture. We refer to the quadratic penalty term of the first constraint P 1 , i k 2 ( x ) = ( ∑ j x i j k − s 1 , i k ) 2 as an example. To implement this constraint, we are requiring an all-to-all connectivity between the qubits representing x _ijk for some fixed i, j and an ancillary qubit s _1,ik . It is beyond current technology to realize arbitrary connectivity between qubits on a physical device, as the interactions are often restricted between neighbors. Specifically, the D-Wave 2000q architecture is based on a Chimera graph, which has only limited connectivity; each qubit being coupled to at most 6 other qubits. The next generation of D-Wave devices (D-Wave Advantage) has an improved connectivity with the Pegasus graph of degree 15 [37], but was not available for testing as of writing time. Nevertheless, it is still very challenging to solve a highly connected problem on either of these architectures. Besides, as shown in Section 2.3, the coefficient λ is typically much larger than the terms in the original objective function, which could out-range the practical dynamic range of current device.

In practice, before the logical Ising Hamiltonian can be sampled, the problem has to be mapped to the physical hardware graph via minor embedding [38, 39]. Minor embedding maps one logical qubit into a chain of physical qubits with negative neighboring interaction strength J _{chain_strength}. This process creates significantly more physical qubits, for example the number of physical qubits needed to embed a N-clique onto the Chimera graph of D-Wave 2000q is O(N ²), i.e., only up to 64 fully connected qubits can be embedded onto over 2000 qubits. Besides, the interaction J _{chain_strength} must be of similar magnitude or greater than the penalties to prevent breaking of chains, which would further address the problem of effective dynamic range. We refer to [38, 40] for more discussion around the optimal value of J _{chain_strength}.

Several previous contributions have tried to find more efficient ways of embedding constraints with Quantum Annealing. Hen et al. propose to design a specific driver Hamiltonian H _D that commutes with the penalty Hamiltonian H _penalty: [H _D , H _penalty] = 0 but not with the original problem Hamilotnian H*: [H _D , H*] ≠ 0 [41, 42]. In this case, the state would only evolve in a sub-space of the Hilbert space that satisfies the constraints. This approach is however not yet implementable on the current version of quantum annealer devices. In another work, Vyskočil et al. suggest to solve the problem of embedding constraints by using ancillary variables and mixed-integer linear programming to find the optimal combinatorial design of these qubits on the hardware graph [43, 44]. Ajagekar et al. apply a decomposition of the problem into a constrained MILP problem and unconstrained QUBO sub-problems that would be solved by quantum annealer [45]. However, these methods either could not be implemented with the current hardware or very hard to generalize to realistic optimization problems with several constraints.

First of all, we generate several instances of the EV-bus scheduling problem with different values of L _time, N _bus, N _pile based on real operational data. The problems are generated by randomly sampling from 15 typical bus operational schedules and randomly setting the number N _pile and power ω _j of charging piles in a reasonable range. We refer to Supplementary Figure S1 in the Supplementary Material for a visualization of these typical bus schedules.

To solve the problem via a quantum annealer, the logical problem (QUBO) graph needs first to be mapped to the hardware graph. In this case, the D-Wave-2000q has a Chimera architecture, which possesses a limited connectivity of degree 6. This implies that the number of physical qubits after the mapping is larger than the number of logical qubits in its original QUBO form, due to the minor embedding. We note that the minor embedding problem is NP-hard itself hence a heuristic algorithm is employed here [49]. Indeed, as analyzed in Section 2.4 and Section 3, the highly connected QUBO graph caused by the quadratic penalties increases significantly the number of qubits needed on the hardware to solve the EV-bus scheduling problem, even at a small scale.

The cost of embedding is illustrated in Figures 1A–D, where a problem with ∼ 150 qubits will cost ∼ 1,500 qubits on the hardware. We further observe a quadratic increase of physical qubits in Figure 1A when fixing N _bus = 1, N _pile = 1 and increasing L _time from 24 to 120, and a linear increase in Figures 1B–D where L _time is fixed. Consistently with the analysis in Section 2.3, we find thus that increasing the number of time steps L _time is more challenging in terms of hardware requirement than increasing N _bus, N _pile. This result further emphasizes that the extra qubits needed are mostly due to the quadratic penalty terms in P 3,4 2 ( x ) = ( ∑ j , k = 1 D i n ω j x i j k + … ) 2 , where the number of variables x in the penalty increases with L _time.

After minor embedding, the embedded problem is then sampled using the D-Wave quantum annealer as an efficient solver for QUBO problems. As discussed in Sections 2.3 and 2.4, the penalty coefficients λ _i for the penalties P _i need to be set appropriately to guarantee that the structure of the original problem is preserved. Additionally, the parameter J _{chain_strength} is introduced during minor-embedding to ensure that the physical qubits can be correctly mapped back to the logical qubits after the sampling. Although theoretical bounds, i.e., minimum theoretical values to ensure that the penalized optimal solution is the original one, are provided for λ _i in Section 2.3, in practice, the optimal values for P _i can not be directly calculated. Instead, we explore the trade-off between penalty and chain strength by computing the time to 99%-success, as defined in [14], and chain break fraction for an instance of the problem N _bus = 1, N _pile = 1, L _time = 24 solved with different settings of λ _3,4 and J _{chain_strength}. The coefficients λ _1,2 are constant to 10, since according to our previous results, these values are above their theoretical bounds and are of no significant effect on the efficiency of the solver.

As shown in Figure 2A, the time to 99%-success for this problem takes similar values along the lines J _{chain_strength} ∼ λ _3,4. However, this conclusion does not hold for other problems, as shown in Supplementary Figure S2 in the Supplementary Material, where the optimal λ _3,4 seems to be in the range [50, 800] for problems of larger sizes. We also observe that for λ _3,4 ≫ J _{chain_strength}, almost no solution could be found, which can be explained by an increased chain break fraction, as shown in Figure 2B. Equally for λ _3,4 ≪ J _{chain_strength}: although the quality of mapping is guaranteed by an increased chain strength, yet the efficiency of sampling by D-Wave QPUs is reduced. We conclude that choosing the optimal values for λ _3,4 and J _{chain_strength} is largely an empirical question that is problem dependent, but it is safer to set a penalty λ _3,4 close to the theoretical bounds and a chain strength J _{chain_strength} of the same order of magnitude as λ _3,4.

It is clear from the previous discussion that directly embedding and sampling the problem obtained with the penalty method suffers from several drawbacks and does not enable to solve large-scale problems. An alternative method discussed in Section 3 is proposed to obtain a better scaling with the quantum annealer. Using the Hubbard-Stratonovich transformation, the QUBO Hamiltonian H _QUBO(x) is transformed into an effective Hamiltonian H _QUBO(x, ν) with an additional multiplier ν. Then the remaining problem is to find the saddle point of the effective Hamiltonian, which is the original optimal solution. This method is similar to the Lagrangian dual method in classic optimization, and can reduce the quadratic terms into linear ones in the penalized QUBO energy.

To demonstrate the effectiveness of this method, we solve a bus charging scheduling problem with N _bus = 3, N _pile = 2 and L _time = 48. The QUBO formulation of this problem with quadratic penalty constraints could not be directly embedded onto the 16 × 16 Chimera graph of the D-Wave 2000q sampler. For reference, to embed such a problem onto a 24 × 24 Chimera graph would cost 2,499 physical qubits and the maximum chain length is 31. On the other hand, the dual Hamiltonian as defined in Eq. 25 only requires 320 physical qubits to embed with a maximum chain length of 4. This example illustrates thus clearly that the limitation for connectivity is well mitigated after the Hubbard-Stratonovich transformation.

We investigate then whether the iterative approach can lead to the ground state of the original Hamiltonian, i.e., the optimal schedule with minimum cost and satisfying all constraints. As shown in Figures 3A,B, the optimal solution is obtained after 96 iterations of updating the multipliers ν based on the sampling of the dual Hamiltonian with the D-Wave 2000q annealer. In Figure 3A, we visualize the state-of-charge SOC of the optimal solutions found while sampling the dual Hamiltonian at different iterations. It can be seen that while updating the multipliers ν according to Algorithm 1, the state of charge gradually falls into the interval [0.3, 1], where the constraints are satisfied. Figure 3B displays the distribution of different samples returned by both the D-Wave 2000q annealer and simulated annealing. We note that the multipliers ν are updated according to the lowest energy samples returned by D-Wave’s machine. We observe that the sampled lowest energy of the dual Hamiltonian does not decrease monotonically. This might be caused by multiple reasons. Firstly, the sampling, quantum or classical, is a heuristic process so that the actual ground states of the dual Hamiltonian are not necessarily found at each iteration, as shown in Figure 3B. Secondly, the energy landscape of the dual Hamiltonian might be non-convex, as also discussed by Ohzeki when this method is proposed [32]. Thirdly, the learning rate might be too large.

From an optimization perspective, the aforementioned three points imply that in practice, the task of solving for the multipliers ν is carried out in a complex landscape with noisy evaluations of the gradient. In this context, finding the optimal multipliers, hence the saddle point could be very challenging, especially with a large number of multipliers. To resolve this problem, instead of using a fixed learning rate η = 0.1, we have employed an adaptive gradient method to update the multipliers ν, with ADAM solver for stochastic optimization [50]. To further confirm the effectiveness of the adaptive gradient methods, we have solved the problem multiple times with the two different strategies, ADAM or fixed learning rate. It is observed in Figure 3C that the adaptive gradient approach for updating ν has a significant advantage over fixing the learning rate, requiring about 10 times less iterations to converge. This emphasizes thus that the practical constraint caused by the transformation, namely performing a classical gradient descent, requires careful tuning in order to reduce the additional computational overhead, which ADAM achieves. Yet, we note that it has been proven that ADAM as described in [50] could possibly diverge, even in the convex case [51]. Hence, the convergence of ADAM-type methods has been analyzed theoretically in several recent contributions. For instance, Chen et al. derive a set of sufficient conditions to guarantee the convergence of AmsGrad, a variant of Adam, with a rate of O ( log ⁡ T + d 2 T ) for non-convex stochastic optimization [52], where T is the number of iterations of the algorithm and d the dimension of the problem. Similarly, Zhou et al. achieve a rate of O ( d / T + d / T ) for AmsGrad [53]. The interested reader is also referred to [54–57] for more in-depth analyses of adaptive gradient algorithms. Consequently, although ADAM has shown to be an appropriate choice for solving the bus charging scheduling problem, further work could explore the impact of the adaptive gradient method on the convergence and computational cost of this iterative approach.