Application of Quantum Annealing to Nurse Scheduling Problem

Quantum annealing is a promising heuristic method to solve combinatorial optimization problems, and efforts to quantify performance on real-world problems provide insights into how this approach may be best used in practice. We investigate the empirical performance of quantum annealing to solve the Nurse Scheduling Problem (NSP) with hard constraints using the D-Wave 2000Q quantum annealing device. NSP seeks the optimal assignment for a set of nurses to shifts under an accompanying set of constraints on schedule and personnel. After reducing NSP to a novel Ising-type Hamiltonian, we evaluate the solution quality obtained from the D-Wave 2000Q against the constraint requirements as well as the diversity of solutions. For the test problems explored here, our results indicate that quantum annealing recovers satisfying solutions for NSP and suggests the heuristic method is potentially achievable for practical use. Moreover, we observe that solution quality can be greatly improved through the use of reverse annealing, in which it is possible to refine returned results by using the annealing process a second time. We compare the performance of NSP using both forward and reverse annealing methods and describe how this approach might be used in practice.

Among non-deterministic polynomial time (NP)-hard problems, timetabling or rostering problems represent a number of practically important examples. For example, in operations research, the nurse scheduling problem (NSP) arises when finding the optimal schedule for a set of available nurses over a fixed timetable of shifts. Solutions to NSP are required to respect hard constraints, such as days off and minimum availability, as well as soft constraints, such as minimum shift assignments, for each nurse. Examples of NSP are often cast as linear or quadratic programming problems, depending on the nature of the constraints, but they may also be formulated in terms of unconstrained optimization and solved using search methods, including tabu search.
Quantum annealing (QA) is a metaheuristic method for solving combinatorial optimization problems derived from the principles of quantum mechanics 1 . QA operates by driving the Hamiltonian dynamics of an initial quantum state to a sought-after final state that represents the minimum energy configuration of an encoded optimization problem 2 . In the limit that the dynamics are strictly adiabatic and the Hamiltonian sufficiently complex, this coincides with the universal model of adiabatic quantum computing 3,4 . In practice, however, the adiabatic condition is rarely obtained and the guarantee that the true solution will be recovered is lost. Instead, these quasi-adiabatic dynamics yield the sought-after final state with some non-unit probability and it remains open as to when such behavior can provide a computational advantage 2,5 .
We demonstrate the use of QA to solve NSP and we evaluate its efficiency and accuracy for this problem from empirical results. Our approach uses the commercial quantum annealer available from D-Wave Systems to implement several hard constraints. The D-Wave 2000Q is a commercially available quantum annealing device based on superconducting flux qubits designed to solve quadratic unconstrained binary optimization (QUBO) 6 . It uses the principles of QA operating in the presence of a transverse-Ising Hamiltonian by encoding the problem into a sparsely connected graph expressing the hardware interactions. By reducing NSP to QUBO form and then embedding this problem into the D-Wave processor, we use QA to recover candidate solutions for different problem instances.
Previous applications of QA to unconstrained optimization span a broad variety of topics [7][8][9][10][11][12][13][14] , while some efforts have investigated closely related timetabling problems, such as the job-shop problem 15 . The distinguishing features of NSP include the multiple hard constraints, which make the problem difficult to address for conventional numerical solvers 16 . Equally challenging is the effort to recover approximately optimal solutions, which may fail to find the true minimum with some parameterized tolerance.
Our interest lies in casting NSP as a combinatorial optimization problem in order to validate the capabilities of QA. This may be seen as a first step toward developing quantum-enhanced solvers for other scheduling problems as well. However, the question as to whether QA satisfies the question of quantum computational advantages cannot be addressed here. In addition to limitations on the existing hardware, in both capacity and control 17 , there is poor theoretical understanding of the computational power for QA generally. Prior analyses have suggested that NSP and other classical combinatorial optimization problems belong to the complexity class stoqAQC, for which the relationship to the uniquely quantum complexity class BQP is not yet established 2 . These so-called stoquastic problems may be solved efficiently by strictly classical methods, but this is also unknown 5 . However, there is a clear potential for quantum advantage for non-stoquastic problems, which do not correspond to NSP or other combinatorial optimization problems 18 .
Our presentation is organized as follows. Sec. 2 defines NSP formally in terms of the hard and soft constraints considered here; Sec. 3 introduces QA specialized to the D-Wave 2000Q in terms of a transverse Ising Hamiltonian that encodes an NSP instance. Results from evaluating various problem instances on the 2000Q processor are summarized in Sec. 4, and the influence of reverse annealing to improve these results are discussed. Lastly, we conclude in Sec. 6 with comments for future work.

nurse Scheduling problem
Several alternative definitions for NSP exist, yet its fundamental concepts can be summarized as follows. NSP is a problem to create a rotating roster of nurses working at a hospital while respecting constraints on their availability and level of effort. In the simplest example of a two-shift system, a schedule assigns nurses to the day duty and night duty shifts of the roster. We will apply two hard constraints to this example that enforce a minimum number of nurses assigned to each shift while also ensuring a minimal period of rest between shifts for each nurse. We refer to these as shift constraints and nurse constraints, respectively.
Shift constraints require that a sufficient number of nurses be assigned to each shift. However, the necessary number may depend on the experience of the assigned nurses, as more experienced nurses may be capable of performing more work during their shift. Nurses usually work as part of a group with each in charge of multiple jobs. Hence it is important to allocate the minimum number of nurses needed to cover all the jobs during a shift. In addition, sometimes it would be also required to distribute jobs to nurses evenly. By contrast, nurse constraints represent the need to satisfy the appropriate working condition for each nurse. This includes time between shifts to get enough rest as well as days off and scheduled vacations. We will apply the following constraints to our instance of NSP: 1. Upper and lower limit of the number of breaks. 2. The number of nurses in duty for each shift slot. 3. Upper and lower limit of time interval between two days of duty.

Methods
Quantum annealing with the D-Wave 2000Q. The D-Wave 2000Q is based on quantum annealing, which is a derivative of adiabatic quantum optimization. The latter is based on the time-dependent Schrodinger equation where ψ T ( ) denotes the quantum mechanical wave function of an underlying physical system and H(t) is the time-dependent Hamiltonian that drives the dynamics. A generic form of this Hamiltonian is ( ) 1. Therefore, the quantum state ψ(0) evolves under an interpolation from H 0 to H 1 in order to prepare the final state ψ T ( ). Assuming the initial state is an eigenstate of H 0 , then the adiabatic theorem promises that the quantum state will remain an instantaneous eigenstate of H(t) provided the dynamics evolve sufficiently slow. The latter condition may be enforced by choice of the annealing time T or the schedules. Consequently, we may select the final Hamiltonian H 1 to represent a computational problem in which the eigenstates encode a well-defined solution. More precisely, we will focus on the case in which the ground state encodes the computational solution.
Let s i x , s i z be Pauli spin operators at sites. The D-Wave 2000Q implements an initial Hamiltonian expressed as where V is a set of spin sites, and a final Hamiltonian that takes the Ising spin model form where J ij describe the interaction between sites i, j and h i are the weights of the linear terms. The resulting time-dependent Hamiltonian corresponds to the well-known transverse field Ising model which is used widely in statistical physics to describe complex spin systems. Finding the ground state of the Ising model itself is known to www.nature.com/scientificreports www.nature.com/scientificreports/ be NP-Hard and, therefore, this Hamiltonian is capable of expressing a broad variety of combinatorial optimization problems including, as we show below, NSP. By contrast, the ground state of H 0 is easy to deduce.
Notwithstanding the premise of adiabatic quantum optimization, the technical challenges of reliably preparing a quantum physical system in a pure, zero-temperature quantum state prevent these ideals from being realized in practice. Rather, the behavior of the 2000Q is better approximated by a mixed quantum state evolving as an open system. In addition to the limits on controllability, there is more general concern that knowing the optimal annealing schedule and duration require a priori information about the solution itself. Therefore, quantum annealing is most often treated as a heuristic that may be applied with good accuracy in certain situations.
In our examples below, we interface with the quantum annealer by providing a logical representation of the Ising spin model H 1 to be solved. Additional steps address the transformation of the logical input into a physical representation which can be embedded into the hardware 19 . This step, known as minor embedding, depends strongly on the connectivity of the vertex set V and the sparsity of the chimera layout 20,21 . Consequently, additional auxillary spin variables may be introduced during embedding to ensure logical connections are satisfied 22 . In addition, the ability to express the logical parameters, J i,j and h i , is limited by the dynamic range of the hardware control.
The 2000Q offers a variety of controls for modifying the annealing time and schedules that determine the quantum dynamics leading to solution. Forward annealing corresponds to the process of driving the quantum system from the initial Hamiltonian H 0 to the final Hamiltonian H 1 and then performing measurements. We use forward annealing to compute the outcome of a given problem instance and we repeat this process many times to estimate the frequency with which computed solutions are observed. Reverse annealing builds on the result of a forward anneal by first initializing the processor to a previously computed result and then reversing the Hamiltonian dynamics to anneal backward. This process reintroduces the transverse field and potentially prepares a new, intermediate quantum state. The process is then completed by annealing forward in time to the final Hamiltonian. We investigate the relative accuracy of both forward and reverse annealing to solve instances of NSP.
ising model formulation of nSp. A common approach to casting combinatorial optimization as an Ising model is to first express the problem as unconstrained optimization and, more specifically, as quadratic unconstrained binary optimization (QUBO) 23 . The QUBO form offers a direct mapping into the Ising model using a simple change of variable from binary to bipolar representation. We take this approach to formulate NSP as a QUBO problem with respect to minimization and then transform to the equivalent Ising model. Consider a set of N nurses labeled as = … n N 1, , and a schedule consisting of D working days labelled as specify the assignment of nurse n to day d. We then consider specific instances of the shift and nurse constraints discussed above. For the hard shift constraint, we require that the schedule must ensure at least 1 nurse is assigned each working each day. For the hard nurse constraint, the schedule must ensure no nurse works two or more consecutive days, while the soft nurse constraint requires that all nurses should have approximately even work schedules.
We construct objective functions that correspond to each shift and nurse constraint and then use the sum of these terms to express the QUBO form. We introduce composite indices i n d ( , ) and j n d ( , ) as functions of the nurse n and the day d. We construct the hard nurse constraint by introducing a symmetric, real-valued matrix J such that = + J a i n d j n d ( , ), ( , 1) and zero otherwise. The positive correlation constant a enforces the nurse constraint by penalizing a schedule for nurse n to work two consecutive days. The resulting objective function is quadratic, i.e., J q q i j i j , , and takes its minimal when the hard nurse constraint is satisfied. Note that the nurse constraint can be modified by changing the entries of the matrix J.
We express the hard shift constraint in terms of the required workforce W(d) needed on each day d and the level of effort E(n) available from each nurse n. We seek an equality solution for this constraint by introducing a quadratic function that penalizes schedules with too many or too few nurses assigned. We take a similar approach for the soft nurse constraint by introducing a quadratic penalty for failing to account for nurse preferences in the work schedule. We use F(n) to specify the number of work days that each nurse wishes to be scheduled and G n d ( , ) to define a the preference for nurse n to work on day d.
As a simplified example of the soft nurse constraint, we decompose the preference function into the product 1 In addition, they can also have options whether they may work on weekend/night or not by tuning h 2 (d): The formulation can be more sophisticated by including three-shift systems, distinction of weekdays and weekends regarding burden, or day-off request with priority. We simply require the minimum duty days F(n) for all nurses n are equal to or greater than [ where the positive real-valued numbers λ and γ tune the relative significance of each term. The objective function has its minimum when all the constraints are satisfied and takes on a positive value otherwise. We will assume that the functions E(n), F(n) and W(d) are integer-valued functions of n or d but this is not required. We will require the minimum duty days F(n) for all nurses n are equal to or greater than ) means the integer part of x.
We next transform the QUBO expression in Eq. (7) into an equivalent Ising spin model. This requires changing from the binary variables q i to the bipolar spin variable = − s q 2 1 i i . The resulting quadratic terms are then collected to match the form of the Ising spin model in Eq. (4). Notably, the connectivity between spin sites is determined by the relatively sparse nurse constraints J and F(n) as well as the shift constraints set by E(n), W(d), and G(n, d).

Results
Using the Ising model for NSP, we solve this combinatorial optimization problem with the D-Wave 2000Q. Our implementations use the randomized embedding algorithm based on work by Cai, Macready and Roy and implemented in the D-Wave software toolchain 24 . We first discuss results obtained using forward annealing and then with reverse annealing. For our studies, we fix the annealing time for a single sample to 200 μs and we collect 1000 samples per problem instance to estimate the solution frequency. We present results for N = 3 and 4 nurses and D = 5-14 days. Throughout our study, we fix the parameters γ = 0.3 and λ = 1.3, and for the soft nurse constraints, we use the simplest example with = h n ( ) 1 ( ) 1 and the penalty = a 7/2. The results of simulated annealing show that ground states were found well when those vales of a, γ and λ were used (see Fig. 2).
A graphical representation of the computed schedules are shown in Fig. 1 of two solutions obtained from forward annealing for an instance with = N 3 and = D 4. The schedule in Fig. 1 (left) satisfies all the constraints of this NSP instance, whereas the result in Fig. 1 (right) does not satisfy the nurse constraint that every nurse has an assignment. Consequently, solution (left) minimizes the energy of the Ising spin model, i.e., it represents a ground state, while the result in (right) is necessarily higher due to the penalty term. Similar results are obtained for all the instances tested here, and in the following, we present statistical estimates for the frequency with which the computed schedules satisfy all the constraints. Of course, the ground state may be degenerate and multiple satisfying schedules are possible.
Forward annealing. We first consider the quality of solutions obtained by forward annealing. Figure 2 presents the probability of finding a ground state with D-Wave 2000Q (left) and simulated annealing (right) for the case of N = 3 and 4 over a schedule of D days. While there is significant probability to recover completely satisfying solutions for smaller schedules, we observe that the probability to satisfy all constraints falls below our sampling level for larger schedules. In particular, we did not observe completely satisfying solutions from D-Wave 2000Q for N = 4, D = 10, 11, 12, 13, 14, whereas some ground states were obtained with D-Wave 2000Q for N = 3, Comparing those results, we may say simulated annealing has a considerable advantage.
We next evaluate the deviation of computed solutions from the completely satisfying solution using the Hamming distance. For these calculations, we express the schedules as ordered binary vectors of size ND and we sum the elements of their inner product to calculate the number of position in which they agree. Subtracting this number from ND yields the Hamming distance, which measure the number of position in which the vectors disagree. In particular, a value of zero indicates exact agreement between the computed schedule and the observed ground state solution. Figure 3 plots the Hamming distance between consecutive solutions obtained by means of the D-wave processor. It is apparent from these statistics that the average Hamming distance is well below the maximum value if ND but significantly greater than 0 for all forward annealing solutions. We will present a comparison with reverse annealing solutions in the next section.

Reverse annealing.
We study the performance of reverse annealing when using the schedules computed by forward annealing. Reverse annealing is a heuristic methods to improve the frequency with which a computed solution satisfies the problem constraints. The procedure basically consists of the following three steps: (1)  An example of a reverse annealing schedule used in our study is shown in Fig. 4. The schedule reverses to = . s 0 6 at time 2 μs, holds for h t = 10 μs and then anneals forward at 12 μs to end at 14 μs. Though there are various choices of the target s and hold time h t , we have found the heuristic to work for our problem instances when ∈ .
.   www.nature.com/scientificreports www.nature.com/scientificreports/ 12, 13, 14 are exactly zero. While it is noted that solutions were successfully improved for most of the cases by reverse annealing, the average success rate decreased. By this methodology, the probability of finding ground states is highly related to the frequency distribution of forward annealing, since initial inputs of reverse annealing are randomly chosen. Comparing Fig. 5 with the probability of forward annealing given in Fig. 2, we see that the accuracy is not always improved by reverse annealing. This is understood as follows. When accuracy is improved by reverse annealing, this is because the relatively many forward annealing solutions which are close to ground states are obtained and those solutions are frequently picked up as initial inputs of reverse annealing.
We examine the corresponding Hamming distance between consecutive solutions computed with reverse annealing in Fig. 6, which shows the mean Hamming distance and their standard deviation. Comparing Fig. 6 with Fig. 3, we notice that overall reverse annealing decrease the mean and the standard deviation of the   www.nature.com/scientificreports www.nature.com/scientificreports/ Hamming distance becomes close to 0, which implies that reverse annealing successfully refined solutions and grouped them together.
We now discuss application of reverse annealing to the low-lying energy states of the forward annealing solution distribution. These states represent the best solutions obtained from forward annealing and, therefore, may be the most likely to be effected by the process. This approach can be useful to estimate the essential effect of reverse annealing, since it is natural to ask whether reverse annealing works well when an initial state already close to the true ground state.
We used the lowest energy solution from each forward annealing distribution as the initial state of reverse annealing. Figure 7 presents the probability of finding a solution that completely satisfies the constraints of the NSP instance from reverse annealing with these initial states. In a process of analyzing data, we noticed that even if a satisfying solution is used as the initial input, it is not guaranteed that computed solution will be the same state. There are many local minima and reverse annealing may get trapped in such local minima depending on the schedule parameters. We find that the accuracy is close to 100% when the ground state is used for an input, otherwise the probability decreases. Figure 7 is based on different lowest energy states, some of which are ground states, of forward annealing.
Some ground states were successfully obtained by reverse annealing even though the initial input was not a ground state. This shows that it is possible to improve solution quality by reverse annealing for the case of using the low-lying energy states as input. However for the = N 4 case, ground states were not observed for D = 10, 11, 12, 13, 14. Figure 8 show the mean Hamming distance and the standard deviation for these results. Both statistics are greatly reduced compared to the forward annealing results as expected from the improved accuracy.
A lower value of s_target corresponds to a larger transverse field, which leads to a broader search of the solution space. A greater value of hold_time (microseconds) implies that the reverse anneal is given more time to explore, which leads to a broader and deeper search of the solution space.

Three-Shift System and Additional Constraints
Our example of a two-shift NSP can be extended to address a three-shift system in which each day is divided into shifts labeled as daytime, early nighttime and late nighttime. Both the two-shift and three-shift systems are popular in medical clinic settings. For the three-shift scheduled, the Hamiltonian is modified to include the new as follows:  www.nature.com/scientificreports www.nature.com/scientificreports/ and The number of duty on late nighttime shift should be finely tuned to sustain a healthy work environment of nurses. In addition, the requests of day-off for individual nurse can be incorporated to the Hamiltonian as an external magnetic field term, where g(n, d) is defined by the priority of the day-off request for a specific nurse and working shift high priority 1 5 middle priority low priority (11) By tuning the positive parameter η, the day-off request is highly considered within scheduling. This constraint should be treated as a soft constraint, whereas the first two terms correspond to hard constraints.
To demonstrate that this formulation is capable of solving a three-shift NSP, we optimized the NSP with 3 nurses and 6 working slots in a three-shift system by taking the day-off request of nurse into account. In Fig. 9, three different types of day-off request are incorporated for each optimization. While Fig. 9(a) has no day-off request, Fig. 9(b,c) show the job shifts according to the day-off request. In Fig. 9(b), all day-off request (n, d) = (1,4), (2,6) and (3,5) are satisfied. On the other hand, in Fig. 9(c), not all requests are satisfied and some constraints are given priority over the day-off requests of lowest priority. As in Table 1, setting the parameter η lowest among other parameters leads day-off request as most soft constraints. The optimization is achieved by QUBO and its solutions are guaranteed to be a ground state of the Hamiltonian by the brute force search method.

Discussion and Conclusion
We have presented a formulation of the nurse scheduling problem (NSP) as an Ising spin Hamiltonian that can be solved using quantum annealing with the commercially available quantum processor. Our approach casts the constraint satisfaction problem into QUBO form and then translates directly to the Ising model. We estimate the frequency for satisfying all constraints using a variety of quantum annealing techniques performed on the D-Wave 2000Q quantum   www.nature.com/scientificreports www.nature.com/scientificreports/ annealer. Our results demonstrate that for a fixed sample size and fixed annealing duration T the probability of successfully matching all constraints decays with increasing size of the schedule D and increasing size of the roster N. Reverse annealing methods do improve the probability of success but not uniformly. This is an encouraging sign that the quantum annealer is suited to solve a complementary set of practical problems but we do not address the question of performance relative to conventional solver method. Current realization of quatnum annealers are insufficient in capacity to explore realistic sizes of N and D. However, it may prove possible to decompose larger problems into a set of smaller optimization problems that could then be solved 25 . For example, we have tested an open-source software package called qbsolv, which is based on a variant of tabu search. We found that qbsolv was capable of find satisfying solutions to NSP instances of size (N, D) = (4, 160) using the same parameters used in this work.
Though quantum computation is in an early stage, the approach we have presented may be extended to more general scheduling problems, and anticipate several future directions to improve this methodology. While we have confirmed that reverse annealing can be useful to obtain ground state solutions, it is unclear the necessary parameters to ensure this improvement and further empirical investigations are warranted. In general, reverse annealing depends on several schedule variables and it is natural to try to investigate more on their dependence on those parameters. While our study has been limited to the schedule shown in Fig. 4, different schedule parameters may influence the results. An example of such a study is given in 26 . Our work is a first step toward solving scheduling problems with a quantum annealer. In addition, a way of accommodating the hard constraints should be improved. Though simulate annealing can find solutions which satisfy both soft and hard constraints quite nicely, we cannot deny the possibility that the hard constrains were not treated well due to quantum tunnelling. This problem in addition to noise of the processor would be responsible for the low probability of finding ground states. This requires further investigation (see 27 for example). There would be various ways to extend and generalize our method on the NSP to different scheduling problems.