IMPROVING THE SHIFT-SCHEDULING PROBLEM USING NON-STATIONARY QUEUEING MODELS WITH LOCAL HEURISTIC AND GENETIC ALGORITHM

. We improve the shift-scheduling process by using nonstationary queueing models to evaluate schedules and two heuristics to generate schedules. Firstly, we improved the ﬁtness function and the initial population generation method for a benchmark genetic algorithm in the literature. We also proposed a simple local search heuristic. The improved genetic algorithm found solutions that obey the delay probability constraint more often. The proposed local search heuristic also ﬁnds feasible solutions with a much lower computational expense, especially under low arrival rates. Differently from a genetic algorithm, the local search heuristic does not rely on random choices. Furthermore, it ﬁnds one ﬁnal solution from one initial solution, rather than from a population of solutions. The developed local search heuristic works with only one well-deﬁned goal, making it simple and straightforward to implement. Nevertheless, the code for the heuristic is simple enough to accept changes and cope with multiple objectives.


INTRODUCTION
Emergency service systems (ESSs) still see workforce scheduling as a challenge and have even been studied within operations research/management science for more than fifty years (Green & Kolesar, 2004;Simpson & Hancock, 2009). The main objective of workforce scheduling is to minimize operational costs while maintaining adequate service level. However, the daily operations of ESSs are surrounded by randomness and require studies to consider variations in the morning, afternoon, and night operations (Schmid, 2012; Souza et al., 2015).

IMPROVING THE SHIFT-SCHEDULING PROBLEM USING NON-STATIONARY QUEUEING MODELS
The requirements often respect a service level measure, which we consider as the probability that there is going to be at least one idle server when a new user arrives in the system. Gillard and Knight (2014) uses Singular Spectrum Analysis (a technique to decompose a time series into oscillatory components and noise) to generate demand forecasts and requirements because they consider that other means do not capture seasonal variations properly. Defraeye Thompson (1993Thompson ( , 1997, and Jennings et al. (1996). Among papers about ESS, one can cite Kolesar et al. (1975), using data from a police precinct.
As the problems themselves, finding a solution (schedule) that obeys all requirements and optimizes the goal is a great challenge and many methods arise in the literature. Alfares (2004) categorized solutions of scheduling problems in 10 groups: manual solutions, integer programing, implicit modeling, decomposition, goal programming, working set generation, linear-programing based solution, construction/improvement, metaheuristics, and others. what happens the moment a busy server is scheduled to leave. The literature identified two possibilities for this moment, which are called end-of-shift disciplines. Firstly, the service is preempted and the user rejoins the head of the queue (preemptive discipline). Secondly, the server completes the job before leaving (exhaustive discipline). Ingolfsson (2005) points that the exhaustive discipline is often more realistic when users are humans. However, the exhaustive discipline still requires more advances to represent mealtime breaks (Beojone, 2017). Those interested on nonstationary queueing models may find detailed explanation for them and their numerical solutions in Schwarz et al. (2016), Defraeye and Van Nieuwenhuyse (2016) and the references therein.

METHODOLOGY
This section presents the studied case, a benchmark problem from the literature. It also presents the schedule evaluator (non-stationary queueing model). Further in this section we present the genetic algorithm and the proposed improvements for a shift-scheduling problem. Finally, we present the proposed local search heuristic.
A portion of our study is inspired on Ingolfsson et al. (2002). Therefore, it is imperative to clarify the points of convergence between this paper and Ingolfsson et al. (2002). Firstly, we used the same benchmark problem from New York's police. Secondly, we used the same queueing model. Finally, we use a similar construction for the genetic algorithm, which we improve in the next subsections. Besides the improvements on the genetic algorithm, this paper also contributes with a new local search heuristic, compared to Ingolfsson et al. (2002).

Data and studied case
Kolesar et al. (1975) studied a police precinct in New York to develop a shift-scheduling problem based on stationary queueing models. Ingolfsson et al. (2002) revisited this case to implement a shift-schedule problem using a genetic algorithm, which is explained in more detail ahead. Table  1 summarizes the main information about the New York police case seen in Kolesar et al. (1975) that will be used in our study, such as the arrival rates and the working hours for each possible tour. Note that there are three major shifts, with four mealtime breaks between 2 and 5 hours after the shift's beginning, for a total of 12 possible shifts, as used in Ingolfsson et al. (2002). Service time is considered constant throughout the day, with a constant rate µ equal to 2 calls/hour. This information is the same as used in both Kolesar et al. (1975) and Ingolfsson et al. (2002).

Schedule evaluator (non-stationary queueing model)
According to Gillard and Knight (2014), it is usual to ESSs to operate under stochastic conditions, so it is absolutely important to be able to plan for such situations. Also, these stochastic conditions are time-varying, and it is reasonable to approximate an ESS with a M(t)/M/s(t) queueing model. For those unfamiliar with Kendall's notation, M(t)/M/s(t) represents a queueing system with a nonhomogeneous Poisson arrival process M(t), a negative exponential service time M, and a time-varying number of working servers s(t). Arrival and service processes are represented by the λ (t) and µ rates, respectively, in calls/hour.
As servers have defined shifts, it is necessary to define what happens in case a server is providing service the moment he is scheduled to leave. The behavior assumed by servers in this circunstance is called end-of-shift discipline. As Ingolfsson et al. (2010) did, we adopted preemptive end-ofshift discipline to facilitate the comparison of our computational results. Hence the non-stationary queueing model used in this paper is the same as the one seen in Ingolfsson et al.    Let π n (t) denote the transient state probability that n customers are in the system at time t. We can write the Chapman-Kolmogorov forward differential equations for a time-dependent multiserver queue as seen in Equation (1) to describe the exact behavior of such a system. The Chapman-Kolmogorov forward differential equations were solved, numerically, using the Runge-Kutta method.
In this paper, nonstationary queueing models are firstly aimed at calculating instantaneous delay probabilities and hourly mean delay probabilities. To do so, we define delay probability π D (t) as the probability that an arriving new user must wait before being served. Our objective is to guarantee a lower π D (t) than a given probability P for all values of time t. Equation (2) shows the delay probability as calculated in nonstationary queueing models.
We also need to calculate the hourly mean service level, defined as follows in Equation (3): As λ (t) is constant on the interval [t,t + 1), t∈ {0, 1, 2, 3, . . ., 23}, we can simplify the Equation We use the schedule evaluator (non-stationary queueing model) to support the decision process of both heuristics from subsections 3.3 (genetic algorithm) and 3.4 (local search).

Genetic Algorithm
Many characteristics of shift-scheduling problems lead us to the use of heuristics to solve them, such as the size and complexity of the search space and the non-linear behavior of delay probabilities, along with its discontinuities. We started this application with a standard genetic algorithm.
Goldberg (1989) presented the genetic algorithm. Its main idea is based on Darwin's evolutionary theory, so it needs a population of individuals that will compete to select the most adapted individual. Each individual might be selected to breed with another one. The most adapted individuals will have higher chances to breed and will thus pass its characteristics to the next generations. The least adapted individuals will have smaller chances to breed and will thus tend to disappear. Individuals are also susceptible to random mutations, altering their chances to survive and breed (for more information, see Goldberg, 1989).
The algorithm uses chromosomes to represent individuals, and each individual is a potential solution for the problem. The most adapted individual to have passed through the population at any time will be chosen as the final solution. These solutions are encoded in the form of binary vectors, as shown in Figure 1. With a fixed number of individuals (called the population), a new generation is created as an offspring of the previous population (no members of the previous generation are included in the next generation). To select individuals (chromosomes - Figure 1) for breeding, we measure their adaptation by calculating a fitness function. Then, two distinct individuals are selected to breed with replacement. Higher fitness leads to higher chances of being selected. Regarding the breeding process, selected individuals have their chromosomes pass through a process called crossover, in which their chromosomes are combined to create two new chromosomes (individuals). After crossover, these new chromosomes may suffer "mutation" at any point, where a bit may change from 0 to 1 and vice-versa. When a new generation emerges, we have an iteration, and the selection process begins once more until a number of iterations is reached or another stop criterion becomes true.  The seed for the initial population is the solution of the scheduling problem using SIPP method and integer programing (Dantzig, 1954). This solution method is referred in the remaining of the paper as "SIPP-IP method" and the solution itself (shift schedule for servers) from this method is referred as "SIPP-IP solution".
The process for generating the initial population used here was proposed by Ingolfsson et al. (2002). The initial population size is 60 individuals. The initial population was generated starting from the SIPP-IP solution. The SIPP-IP was one of the members of the initial population, and the remaining 59 individuals were generated randomly as follows: (1) we calculated the mean and variance of the SIPP-IP solution (number of servers in each shift is the piece of data for this calculation); (2) we sampled from a normal distribution with the calculated mean and variance for each shift (defining new shift schedule for all shifts); and (3) sampled numbers were rounded to the nearest integer between 0 and 7 (limits of the adopted chromosome binary encoding; see Figure 1). The mutation probability for each bit was set to 1%, and the maximum number of iterations was set to 100.
Minimizing labor costs here means fewer vehicles, and maintaining an adequate service level means to keep the delay probability lower than 10%. Thus, our fitness functions were built using the same three parameters as in Ingolfsson et al. , and it is presented as Equation (5). M indicates a number large enough that the fitness function will not become negative. This fitness function is referred as FITNESS0 in the remainder of the paper. The second stop criterion (besides the number of iterations) for FITNESS0 was a difference between the fitness of the current population's best individual and average fitness of the population lower than 0.5.
The fitness values of individuals were linearly scaled to ensure that the average fitness value would not change (Ingolfsson et al., 2002). Hence, the selection of individuals for breeding used the scaled value. The purpose of scaling fitness values was to avoid premature convergence during the first iterations and to avoid a random walk among the mediocre individuals in late iterations (Goldberg, 1989). Scaling process is also useful to avoid some bias caused by using a value of M in the fitness function.

FITNESS function choice
However, FITNESS0 puts too much attention to MAXPROB, which cannot indicate whether the service level constraint is respected or not. In the other hand, NUMHOURS is an interesting option, because it only assumes values different from 0 if the service level is higher than the constraint limit. Therefore, as an option for FITNESS0, we worked on a fitness function that tries to obey the service-level constraint at all costs. Equation (6) shows FITNESS1, built using the same parameters as before, except for MAXPROB.

Initial population generation
When evaluating the population along the iterations, the breeding process comes across several individuals that are certainly not optimal. However, due to scaling and the nature of genetic al-gorithms, the chromosomes of these individuals might pass through generations. In a Markovian process, as the queueing models we use, it is very unlikely that a situation in the beginning of the first shift will interfere in the third shift, and even in the second shift. For this reason, we propose another method to generate the initial population, focused in the shifts that do not obey the service level constraint. Before generating any initial population members, we ran the nonstationary queueing model to find the tours in which the delay probability constraint was not obeyed. For example, if the constraint was not obeyed somewhere between 2 and 3 a.m., then all tours beginning at midnight and ending at 8 a.m. were chosen because police officers are in the middle of their shift (mealtime breaks were considered as part of their shifts). This strategy was used to modify only tours that did not obey the service level constraint. Equation (7) shows the proposed method, where CARS 0 i is the number of vehicles assigned to tour i for iteration 0 (initial population); U i and U j are random numbers generated according to a uniform distribution (0 as minimum and 1 as maximum); l is the set of tours chosen; and CARS SIPP−IP j is the number of vehicles assigned to chosen tours in the SIPP-IP solution:

Local Search Heuristic
If we consider that the SIPP-IP solution is a reasonably good initial guess to the problem, it is straightforward to think about using a local search heuristic to find a solution. The module uses the SIPP-IP solution and tries to satisfy the delay probability constraints through minor changes to the current solution (for an introduction to local search algorithms, see Luke, 2013). The main idea is to try to change servers' tours to satisfy the delay probability constraint. In case it is not possible to fulfill the constraints only changing servers' tours, we add one new server to a chosen tour and try to manage it until the delay probability constraint is obeyed. In the other hand, if the SIPP-IP solution fulfills the delay probability constraint, it will be replaced by a poorer SIPP-IP solution, which is calculated after relaxing the delay probability constraint. Table 2 shows the notation used for the local search heuristic. Table 3 illustrates the proposed local search heuristic through a pseudocode. It starts on iteration k ← 0 and calculates a solution from the SIPP-IP. At this point, the SIPP-IP considers a rC D to ensure that C D is still active on the nonstationary queueing model (measured on the schedule evaluator). Then, the code set maxt as 24 hours. Note that the primary stop condition is to make C D inactive. The code verifies whether maxt is acceptable (capable of identifying a mt s w inside [0, maxt)). Then, it seeks the mt s w in which π k D (mt s w ) > C D , defining the tour w as one reference. In the next step, the code looks for y, where y∈W S w and with the lowest hourly mean delay probability (π D ) on its mealtime break, y k sched . Then, it takes one server from tour w and adds one server to tour y, if possible. If this exchange turns π D t k y or π k D t k w worse than π k−1 D t k−1 w , then the script will cancel the server exchange and will not allow the same shift w to be chosen, π k D (t) Delay probability on iteration k at epoch t.
π k D (t) Mean delay probability on iteration k at period [t,t + 1).
by limiting the maxt to mt s w − 1. If there is no acceptable time horizon, the tour y receives a new server and starts a new iteration, resetting maxt to 24 hours.

COMPUTATIONAL RESULTS
This section presents the results computed for the aforementioned studied case with the genetic algorithm and the proposed improvements. Ahead in this section we present the results for the local search and some comparisons between both solution methods in a sensitivity analysis. All computations ran in a Windows 10 Home, with a processor Intel Core® i5 7300HQ 2.50GHz and 8GB of RAM memory. Figure 2 shows the first results for configuration explained here. To obtain more consistent data, we ran the genetic algorithm and calculated an average value for each measure. Hence, Figure  2 does not show actual values for any run of the algorithm but an average value for all of them. Max., min., and avg. are the highest, lowest, and average FITNESS0 values for the population in each iteration, respectively. In observing the reference values of Figure 2, one can see a slight improvement to the average and minimum values for FITNESS0. However, the highest FITNESS0 value is in the initial population (iteration 0) and then falls and increases again, but it never reaches the initial solution value for FITNESS0. On 25 out of 30 runs, the best solution found was the SIPP-IP solution. In four runs, the genetic algorithm found a solution with 28 CARS, and in one run, it found a solution with 27 CARS. None of the 30 solutions obeyed the delay probability constraint. It is important to note that the SIPP-IP solution found (number of servers for each shift) was different from that shown in Ingolfsson et al.  while maxπ D (t) < rC D , do k ← k + 1; rC D ← rC D + 0.05; Run SIPP-IP using rC D as delay probability constraint to obtain server requirements; Schedule Evaluator; end while while maxπ k D (t) > rC D , do maxt ← 24 hours; while maxt > mt e 1 , do Identify last mt s w , such that mt s w ≤maxt and π k D (mt s w ) > C D ; Identify tour mt s y such thatπ k D mt s y ≤π k D (mt s x ) , ∀x∈W S w ; k ← k + 1; n k x ← n k−1 x , ∀x∈W ; if n k w − 1≥0, then n k w ← n k−1 w − 1; n k y ← n k−1 y + 1; Schedule Evaluator; if π k D (mt s w ) > π k−1 D (mt s w ) Or π k D mt s y > π k−1 D (mt s w ), then n k w ← n k−1 w ; n k y ← n k−1 y ;      It is also important to note that changing the fitness function had an effect on runtimes. While FITNESS0 had an average runtime close to 400 seconds, FITNESS1 had 600-second runtimes.

Genetic algorithm
There was also a difference in variation, with FITNESS0 having greater variance and FITNESS1 having a smaller one. We stored previous solutions in a database, which was erased at each new run. FITNESS0 found its solutions after looking at 2,500 solutions, while FITNESS1 looked at 3,200 solutions, on average.  Figure 6 shows the average reference results for FITNESS1 after using the alternative method for generating the initial population. Note that the initial population had good reference values, even for the poorest solution. In 20 out of 30 runs, the best solution was found among the initial population and the first iteration. This solution uses 30 CARS and has 0.09 MAXPROB. For comparison, using the previous initial population generation method, this solution was found only eight times, versus 26 times using the new method. Also, no solutions were found after the eighth iteration. The upper solutions decreased their FITNESS1 along the iterations, even though the average solutions had only slight improvements.   Larger markers indicate the final solution, and the numbers linked to them denote at which iteration the final solution was found. Hence, no more than 6 iterations separate the SIPP-IP solution (iteration 1) from the final solutions. This solution is final, and the local search will not find any other solution for this problem because there were no random choices or steps. Note that, although the algorithm tries to change servers from tours in the best way (taking a server with a mealtime break during a busy period and put it into a tour with a calm mealtime break), it always found worse solutions than the SIPP-IP solution. In the other hand, after adding new servers, rearranging them among tours decreased MAXPROB. Moreover, it is interesting to note that the simple idea of adding servers based on their mealtime breaks is a good strategy to reduce delays.

Comparisons
This final set of tests was aimed at comparing the results from the genetic algorithm and the proposed local search algorithm. The comparisons focused on increases and decreases in demand. The service times remained the same.    Figure 9 explores the delay probability calculated for final solutions from the genetic algorithm, the proposed local search heuristic, and the SIPP-IP. In none of the cases did the SIPP-IP solution have acceptable performance. Also, as the demand increased, the genetic algorithm found more non-acceptable solutions (delay probability higher than 0.10).  Another important point regards runtimes. Figure 10 shows the correlation between the average runtimes and the number of iterations for the proposed local search heuristic. The proposed local search heuristic had almost constant runtimes, much smaller than those seen for the genetic algorithm (about 600 seconds, as seen in Figure 5).

FINAL CONSIDERATIONS
We revisited the police precinct study from Kolesar et al. (1975) and Ingolfsson et al. (2002) to improve the shift-scheduling process using nonstationary queueing models to evaluate schedules. The genetic algorithm used in Ingolfsson et al. (2002) received improvements over the fitness function and initial population generation method. We also proposed a local search heuristic, in the form of a hill-climb heuristic. This simple heuristic passed through an initial test, and its results were compared to those of the genetic algorithm through a sensitivity analysis by decreasing and increasing arrival rates on a small-scale system.
Our work has important implications for managers to improve the service quality and efficiency of an ESS. Firstly, the improved genetic algorithm found solutions that obey the delay probability constraint more often and have close to the optimum number of vehicles. Two important features were responsible for these improvements: the new fitness function, which forced the choice of schedules that obey the delay probability constraint, and the new initial population generation method. We illustrated through a series of tests how these new features affect convergence and the genetic algorithm's final solutions. In addition, results from the improved genetic algorithm led us to the proposed local search heuristic. This heuristic also finds close to optimum schedules with much lower computational expense, especially under low arrival rates. For comparison, the improved genetic algorithm found solutions after 600 seconds, on average, while the proposed local search heuristic found them in less than 5 seconds (even less than 1 second, in some tests). Both the genetic algorithm and local search heuristic determined the appropriate number of vehicles allocated in each possible tour, while considering 1-hour-long mealtime breaks.
Differently from a genetic algorithm, the developed heuristic will always come to the same results, as it does not rely on random choices. Furthermore, it finds one final solution from one initial solution (the SIPP-IP solution), rather than from a population of solutions. As the proposed local search heuristic does not use chromosomes, its solutions are not limited to the limits that the chromosome encoding might impose on the problem.
The developed heuristic works with only one well-defined goal, making it simple and straightforward to implement. It also relies on the intuitive idea of managing servers' tours based on their mealtime breaks. Nevertheless, the code for the heuristic is simple enough to accept changes and cope with multiple objectives.
There is a wide field for further research. Moreover, using Chapman-Kolmogorov equations is computationally expensive, so we encourage the use of approximate methods like randomization/uniformization (Grassmann, 1977) and the stationary backlog-carryover approach (Stolletz, 2008). The search for better and more efficient heuristics and algorithms is encouraged. Also, there is still room for improvement to the genetic algorithm, such as by using approximations instead of exact models or better fitness functions using aggregate measures. We suggest the use of the proposed local search heuristic for small-scale systems, when considering the use of complex performance measures, to take advantage of the simplicity of our method.