A Generic Two-Phase Stochastic Variable Neighborhood Approach for Effectively Solving the Nurse Rostering Problem

In this contribution, a generic two-phase stochastic variable neighborhood approach is applied to nurse rostering problems. The proposed algorithm is used for creating feasible and efficient nurse rosters for many different nurse rostering cases. In order to demonstrate the efficiency and generic applicability of the proposed approach, experiments with real-world input data coming from many different nurse rostering cases have been conducted. The nurse rostering instances used have significant differences in nature, structure, philosophy and the type of hard and soft constraints. Computational results show that the proposed algorithm performs better than six different existing approaches applied to the same nurse rostering input instances using the same evaluation criteria. In addition, in all cases, it manages to reach the best-known fitness achieved in the literature, and in one case, it manages to beat the best-known fitness achieved till now.


Introduction
In this paper, the problem of nurse rostering is faced, which refers to the schedule of the personnel's shift in a hospital.This problem belongs to the wide category of timetabling problems.These problems deal with the allocation of resources to specific timeslots so that some specific constraints are satisfied and the created timetables/rosters are valid and effective.According to each case, the constraints, the sources and the elements defining the effectiveness of each timetable/roster are determined.These timetabling/rostering problems are non-deterministic polynomial time (NP)-complete in their general form [1], as far as their computational complexity is concerned, meaning that the difficulty to find a solution rises exponentially to their size and a deterministic algorithm, giving an acceptable solution in polynomial time, cannot be found [2,3].Therefore, alternative optimization methods, namely metaheuristics, have been developed in order to reach a (near) optimal solution for various kinds of the nurse rostering problem [4,5].Metaheuristics comprise a major class of approaches to solve the nurse rostering problem.They have been designed in order to cope with complex optimization problems in cases where other optimization methods have failed to be either effective or efficient.The main advantages of metaheuristic methods are their effectiveness and general applicability.In the literature, a lot of heuristic methods have been developed for dealing with the nurse rostering problem: genetic algorithms [6][7][8], tabu search [9,10], simulated annealing [11,12], variable neighborhood search [13][14][15], scatter search [16,17], iterated local search [18,19], particle swarm optimization [20], memetic algorithms [21], ant colony optimization [22], etc.
The algorithm presented in this contribution comprises a heuristic method to solve the nurse rostering problem.More precisely, it is a stochastic variable neighborhood approach, which uses three different swap mechanisms, which are different from other swap mechanisms presented in the literature [13].The use of these three swap operators by the proposed algorithm enables it to search in three different neighborhoods of the problem's search space.The reason why we decided to use a variable neighborhood search algorithm in order to solve this specific problem is that variable neighborhood search algorithms have been widely applied in many multi-objective optimization problems having very satisfactory results [23][24][25].The innovation of the proposed algorithm is two-fold.First, although there are plenty of variable neighborhood algorithms applied to scheduling and timetabling problems in the literature [13,[26][27][28][29], there is no two-phase stochastic variable neighborhood approach, to the best of our knowledge, applied to the nurse rostering problem.This was our main motivation in order to design and apply a two-phase stochastic variable neighborhood algorithm, so as to solve effectively and efficiently the nurse rostering problem.The second novelty of the proposed algorithm is the application of a -stochastic moving segment grouping swap‖ (see Subsection 3.2), which is innovative, to our knowledge, and different from other types of swaps presented in the literature [13,26].
Therefore, in this contribution, a new two-phase stochastic algorithm based on variable neighborhood search [30,31] has been designed, developed and applied to the nurse rostering problem.The generic two-phase stochastic variable neighborhood algorithm proposed has been used in order to create feasible and efficient schedules of the personnel's shifts in many different hospitals having different types of constraints.In order to demonstrate the effectiveness, efficiency and generic applicability of the proposed algorithm, its performance is compared with six other very effective algorithms published in the literature that have been applied to the same problem instances [32][33][34][35][36][37].

Related Work
Valouxis and Housos presented in [32] a hybrid methodology that utilizes the strengths of operations research and artificial intelligence.In particular, an approximate integer linear programming model is firstly solved, and its solution is further improved using local search techniques.Furthermore, a tabu search strategy is applied in order to construct effective and efficient solutions.Li et al., present a hybrid artificial intelligence approach for a class of over-constrained nurse rostering problems, [33] which comes in two phases.The first phase solves a relaxed version of the problem, which only includes hard rules and part of the nurses' requirements for shifts.In the second phase, adjustments with descend local search and tabu search are applied to improve the solution.The algorithm presented in [34] is a shift sequence-based approach that consists also of two stages: (a) Constructing high quality sequences for nurses by only considering the sequence constraints and (b) Iteratively constructing schedules for nurses and the overall rosters, based on the sequences built and considering the schedule and roster constraints.Greedy local search carried out during and after the roster construction manages to improve the (partial) rosters built.Puente et al., present a genetic algorithm approach to solve the medical doctor rostering problem in a hospital emergency department in [35].More specifically, they intend to automate the creation of timetables by applying genetic algorithms in an actual hospital emergency department.Firstly, a heuristic-schedule builder, designed ad hoc to satisfy the hard constraints, produces an initial population of feasible solutions.Afterwards, iteratively, a genetic algorithm obtains new generations of feasible individuals, thanks to the use of a specific crossover operator, based on the exchange of whole workweeks that operates together with a repair function.Musa and Saxena describe in [36] a single-phase goal-programming algorithm for scheduling nurses in one unit of a hospital.The goals represent the scheduling policies of the hospital and nurses' preferences for weekends on and off.Experiments on one unit with 11 nurses resulted in satisfactory results.Finally, in [37], Weil et al., present the efficiency of constraint programming for solving the nurse rostering problem.Experimental results obtained are very satisfactory regarding response time and flexibility of the approach.
The proposed variable neighborhood search algorithm uses the same formalism for modeling the nurse rostering problem, tries to minimize the same fitness function and uses the same performance criteria in order to evaluate the quality of resulted rosters, as the ones used in [32][33][34][35][36][37].Therefore, a straightforward comparison of their experimental results is fair.Moreover, in order to have a fair comparison with these algorithms, we decided to use the exact same input instances used by these six approaches.Computational results showed that the proposed two-phase variable neighborhood search algorithm achieves better results compared to these six very effective algorithms.The comparison was carried out on the basis of seven instances taken from real world situations that were also used as input by the six published effective approaches mentioned above.In one case, the proposed algorithm manages to beat the best-known fitness achieved till now.In addition, in the rest of the six cases, the proposed algorithm manages to reach, for each different instance, the best-known fitness achieved in the literature and demonstrates, experimentally, that in these instances, there is more than one roster that achieves the best-known fitness.
This paper is organized as follows.Section 2 defines the nurse rostering problem and the constraints used, in most cases, in order to evaluate the resulting shift schedules.Section 3 describes the proposed two-phase variable neighborhood algorithm, while Section 4 describes the input data used.Section 5 assesses and compares the performance of the proposed algorithm to that of existing approaches.Finally, Section 6 provides a summary and future extensions.

Problem Definition
The nurse rostering problem has to satisfy a large number of constraints and is affected by many parameters.The entities that are involved in the construction of a feasible and effective solution of the nurse rostering problem are the nurses, the shifts and the time periods.More precisely, nurses have to make some specific shifts in specific time periods.Therefore, in order to create a feasible timetable, for the nurse-shift couple, the time periods that these shifts will take place must be assigned.Constraints regarding the construction of a nurse roster can be divided into two categories: -hard‖ constraints and -soft‖ constraints.When all hard constraints are satisfied, then a feasible nurse roster is constructed, which is a roster that can actually be used by the hospital it was made for.However, the number of soft constraints satisfied is the main factor that affects the quality of a nurse roster.The final aim, of course, is to create a feasible nurse roster while maximizing its quality, i.e., to create a roster that satisfies all hard constraints and, at the same time, satisfies the maximum possible number of soft constraints.

Constraints
There are many different types of nurse rostering problems found in the literature, each having their own constraints.However, in most cases, the hard constraints that must be satisfied in order to keep the nurse roster valid are the following:  All shift type demands during the planning period must be met  The shift coverage requirements must be fulfilled  Each nurse should work at most one shift per day Also, the soft constraints that should be satisfied, in most types of the nurse rostering problem, in order the nurse roster to be considered of high quality are the following:  Maximum number of shifts that should be assigned to a nurse  Minimum number of shifts that should be assigned to a nurse  Maximum number of consecutive working days  Minimum number of consecutive working days  Maximum number of consecutive free days  Minimum number of consecutive free days  Maximum number of hours worked  Minimum number of hours worked  Maximum number of consecutive working weekends  Maximum number of working weekends in four weeks  Number of days off after a series of night shifts  Complete weekends (i.e., if a nurse has to work only on some days of the weekend, then a penalty occurs)  Identical shift types during the weekend (i.e., assignments of different shift types to the same nurse during a weekend are penalized)  Unwanted patterns (i.e., an unwanted pattern is a sequence of assignments that is not in the preferences of a nurse, based on her contract)  Unwanted patterns not involving specific shift types  Unwanted patterns involving specific shift types  Alternative skill (i.e., if assignments of a nurse to a shift type requiring a skill that she does not have occurs, then the solution is penalized accordingly)  Day on/off request (i.e., requests by nurses to work or not to work on specific days of the week should be respected, otherwise solution quality is compromised)  Shift on/off request (i.e., similar to the previous, but now for specific shifts on certain days) As stated in the Introduction Section, in this contribution, the proposed generic variable neighborhood search algorithm is applied to seven different nurse rostering instances each of which belongs to a different type of the nurse rostering problem.For a more detailed description of each type of the nurse rostering problem faced by the proposed algorithm, the reader can refer to the respective references [32][33][34][35][36][37].A detailed description of the input instances used in the experimental results is presented in Section 4.

General Overview of the Stochastic Variable Neighborhood Algorithm
The flowchart describing the general overview of the proposed algorithm is presented in Figure 1.As shown there, the proposed algorithm is a hybrid one consisting of two phases: (a) The first phase deals with the assignment of nurses to working days (b) The second phase deals with the assignment of nurses to shift types At first, the values of the algorithm's parameters are set, namely, the population size, the maximum number of repetition cycles of first phase, the maximum number of generations of second phase and, finally, the swapping probability (see Subsection 3.2).Due to the differences in nature, structure, philosophy and type of hard and soft constraints among input instances, there is a small difference in the swapping probability value used for each instance.Note that the value of swapping probability determines the searching behavior of the algorithm.A high value will cause an exhausting cell swap, while a low value will cause skipping of cell swaps (see Subsection 3.2).Exhaustive experiments showed that for some instances, a lower value of the swapping probability is beneficial to the algorithm, while the opposite holds for others.Except for swapping probability, the values used for the other algorithm's parameters are the same.Table 1 lists the parameters' values for each instance, as well as the time consumed in order to find the optimal parameters of the algorithm.The effect of parameter setting to algorithm performance is investigated in Section 5.1.The experimentation procedure in order to find the optimal parameters of the algorithm is described as follows.At first, the number of cycles in the first phase was determined.For each input instance, we ran five experiments setting the number of cycles equal to 1, 2, 3, 4 and 5, respectively.For all instances, experimental results showed that a number of cycles equal to 1 suffices in order to achieve the best possible results.Next, the population size was determined.For each input instance, we ran five experiments, setting the population size equal to 1 to 5, respectively.For all instances, experimental results showed that a population size equal to 2 achieves the best possible results.The maximum number of generations in the second phase was set arbitrarily equal to 100, which is a very big value, since our main purpose was to reach or even beat the best ever reported roster for each input instance.However, as mentioned in Section 3.3, the user is able to choose the termination criterion he/she likes to apply between the maximum number of generations and the number of generations for which the fitness remains the same; that is, no improvement is reported.
Finally, we determined the value of the swapping probability (p swap ) that leads the algorithm to the best possible results.At first, for each input instance, we ran seven experiments, setting p swap equal to 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0.For the sixth instance, since p swap = 0.4 was the value giving the best results, we further experimented, setting p swap equal to 0.35, 0.36, 0.37, 0.38, 0.39, 0.41, 0.42, 0.43, 0.44 and 0.45.Since p swap = 0.45 was the value giving the best results, we further experimented, setting p swap equal to 0.445, 0.446, 0.447, 0.448, 0.449, 0.451, 0.452, 0.453, 0.454 and 0.455.Since again p swap = 0.45 was the value giving the best results, we decided to use this value.For the fourth instance, since p swap = 0.5 was the value giving the best results, we further experimented, setting p swap equal to 0.45, 0.46, 0.47, 0.48, 0.49, 0.51, 0.52, 0.53, 0.54 and 0.55.Since again p swap = 0.5 was the value giving the best results, we decided to use this value.For the seventh instance, since p swap = 0.8 was the value giving the best results, we further experimented, setting p swap equal to 0.75, 0.76, 0.77, 0.78, 0.79, 0.81, 0.82, 0.83, 0.84 and 0.85.Since p swap = 0.85 was the value giving the best results, we further experimented, setting p swap equal to 0.845, 0.846, 0.847, 0.848, 0.849, 0.851, 0.852, 0.853, 0.854 and 0.855.Since again p swap = 0.85 was the value giving the best results, we decided to use this value.For the second, the third and the fifth instances, since p swap = 1.0 was the value giving the best results, we further experimented, setting p swap equal to 0.95, 0.96, 0.97, 0.98 and 0.99.Since p swap = 0.97 was the value giving the best results, we further experimented, setting p swap equal to 0.965, 0.966, 0.967, 0.968, 0.969, 0.971, 0.972, 0.973, 0.974 and 0.975.Since again p swap = 0.97 was the value giving the best results, we decided to use this value for these two instances.For the rest instances, namely, the first instance, since initially p swap = 1.0 was the value giving the best results, we further experimented, setting p swap equal to 0.95, 0.96, 0.97, 0.98 and 0.99.Since again p swap = 1.0 was the value giving the best results, we further experimented, setting p swap equal to 0.995, 0.996, 0.997, 0.998 and 0.999.Since again p swap = 1.0 was the value giving the best results, we further experimented, setting p swap equal to 0.9995, 0.9996, 0.9997, 0.9998 and 0.9999.Although, once again p swap = 1.0 was the value giving the best results, we decided to use p swap = 0.99995, since a value of p swap = 1.0 would make the execution of swaps deterministic and not stochastic.Using p swap = 1.0 would lead the algorithm to poor diversification and big intensification.To conclude, we noticed that for all instances, setting p swap to a big value causes big intensification, while setting p swap to a small value causes big diversification.
After setting the parameters' values, the initialization of each individual of the proposed algorithm's population takes place.This is done at random with respect to the working requirements of each day.Next, the first phase of the algorithm, which deals with the assignment of nurses to working days, is executed.After the first phase is completed, the best individual found is copied to all individuals of the population, and shift types are randomly assigned to them.Therefore, the input to the second phase of the algorithm is a population of individuals, all of them being equivalent, by means of workings, to the best individual found by the first phase, with shift types randomly assigned to them.Finally, the second phase of the algorithm, which deals with the assignment of nurses to shift types, is executed, leading to the near optimal solution found by the proposed algorithm.
Since the seven nurse rostering input instances, to which the proposed algorithm is applied, are totally different cases of nurse rostering problems, there are significant differences, between them, in nature, structure, philosophy and type of hard and soft constraints.As a result, a different evaluation function is used for each different input instance taking into account different constraints and having different weight values.The specific constraints and the weight values used for each different constraint for each different input instance are the ones presented in [38].However, a general form of the evaluation function applied to all instances can be presented as follows: where the ith Constraint is different for each different input instance, Weight_of_ith_Constraint is the weight of the ith Constraint (which is different for each different input instance) and Times_ith_Constraint_is_violated is the number of times the ith Constraint is violated.In addition, Weight_function(Weight_of_ith_Constraint) equals Weight_of_ith_Constraint, if the ith Constraint is assumed linear, while Weight_function(Weight_of_ith_Constraint) equals Weight_of_ith_Constraint × Weight_of_ith_Constraint if the ith Constraint is assumed quadratic.Whether a constraint is assumed linear or quadratic depends, once again, on the input instance, and it is explicitly stated.At this point, we have to mention that for all input instances, the evaluation function takes into account only constraints concerning working days and day offs in the first phase of the algorithm, while it takes into account all kinds of constraints in the second phase of the algorithm.

The First Phase of the Stochastic Variable Neighborhood Algorithm
The first phase of the proposed algorithm, which deals with the assignment of nurses to working days, is presented in Figure 2. As shown, this phase consists of the execution of procedure Successive_Segment_Swap_Mutation().This procedure is applied to each individual sequentially and is repeated for a specified number of cycles.A detailed description of this procedure along with procedure Selective_Partial_Swap(), which is used by Successive_Segment_Swap_Mutation() in the first phase of the algorithm, is given in the following paragraphs.The procedure, Successive_Segment_Swap_Mutation(), is applied as follows.At first, a list of all nurses L 1 is created at random.Next, for each nurse n 1 in L 1 and for each nurse, n 2 , next to n 1 (i.e., after n 1 in L 1 ), procedure Selective_Partial_Swap() is applied.This procedure, which is described in the next paragraph, is applied between nurse n 1 and other nurses, until no other nurse n 2 exists in list L 1 and is repeated from the beginning for each nurse, n 1 , in list, L 1 .The structure of Successive_Segment _Swap_Mutation() is presented in Figure 3.  Procedure Selective_Partial_Swap(), which in fact implements a -stochastic moving segment grouping swap‖, is applied for each day, d 1 , of the scheduling period as follows.At first, the left extreme (this is always equal to d 1 ) and the right extreme (this is equal to d 2 ) of the cell segment, in which swaps will be performed, are defined.Next, swaps are performed between cells included in a cell segment defined previously for rosters belonging to nurse n 1 and nurse n 2 under a certain probability.After all swaps in the selected cell segment have been performed, the fitness of the roster is computed.If the fitness of the created roster (i.e., after the swaps) is improved, then the swaps are accepted; otherwise, the swaps are discarded, and the roster sustains the structure it had before the swaps.After that, d 2 is increased by one, and swapping cells between d 1 and d 2 for rosters belonging to nurses n 1 and n 2 is repeated as long as d 2 is less or equal to the last day of the scheduling period.The structure of this procedure is presented in Figure 4.At this point, we have to mention that the application of a -stochastic moving segment grouping swap‖ to the nurse rostering problem is innovative, to our knowledge.

The Second Phase of the Algorithm
The second phase of the proposed algorithm, which deals with the assignment of nurses to shift types, is presented in Figure 5.As shown, this phase consists of the sequential execution of the following procedures: Procedure Successive_Segment_Swap_Mutation() is the same with the one executed in the first phase of the algorithm; however, this time, it is executed in order to assign nurses to shift types and not to working days as performed in the first phase.This is the goal of the other two procedures, too.The execution of these three procedures is repeated for each individual of the population for a number of times, i.e., generations.The execution of the second phase of the algorithm is performed, until a specified termination criterion is met.We have implemented two termination criteria, which are the following:  The total number of generations  The number of generations for which the fitness remains the same, that is, no improvement is reported More precisely, at the beginning of the algorithm, the user is asked to select the termination criterion he/she prefers to use: 1.If he/she chooses -the total number of generations‖, next, he/she has to insert this number.2. If he/she chooses -the total number of generations for which the fitness remains the same‖, next, he/she has to insert this number.
In the next sections, we present a detailed description of procedures, Selective_Day_Swap_Mutation() and Random_Segment_Swap_Mutation().

Procedure Selective_Day_Swap_Mutation()
This procedure is applied as follows.At first, a nurse, n 1 , is selected at random.Next, the order of all combinations between nurse n 1 and all other nurses is created randomly.Then, for each day of the scheduling period and for each pair of nurses created at random, a swap is performed between the cells of the current day of each pair of nurses.If the fitness of the created roster (i.e., after the swap) is improved, then the swap is accepted; otherwise, the swap is discarded, and the roster sustains the structure it had before the swap.This procedure is performed for each nurse of the roster.The structure of this procedure is presented in Figure 6.This procedure is applied as follows.At first, a list of all nurses in random order, L 1 , and a list of all nurses in a random order, too, L 2 , are created.Next, for each nurse, n 1 , in list L 1 and for each nurse, n 2 , in list L 2 , procedure Selective_Partial_Swap() is applied.This procedure, which is described in the previous section, is applied between nurse n 1 and each nurse, n 2 , belonging to list, L 2 , until no other nurse n 2 exists in list L 2 and is repeated from the beginning for each nurse n 1 in list L 1 .The structure of this procedure is presented in Figure 7.

Input Data
As stated in the Introduction Section, the proposed two-phase variable neighborhood algorithm is applied to seven different nurse rostering input instances.These instances, which have significant differences with each other, are presented in the following sections.They comprise a set of benchmark data that represents a wide variety of nurse rostering problems with non-trivial properties, which are derived from complete real world complex instances.

Input Instance, Valouxis-1
In this input instance, there are 16 nurses working three daily work shift types, and the planning horizon is 28 days long.The demand is assumed to be the same every week, and the daily requirement for personnel from Monday to Friday is 4-4-2 for the Day, Evening and Night work shifts, respectively, and the demand for Saturday and Sunday is 3-3-2, respectively.The legal work stretches are 2-4 days long, while the least time break between two work stretches is two calendar days long.In the planning horizon, all individuals must have at least one Sunday rest, while the total Day, Evening and Night work shifts must be 5-8, 5-8 and 2-5, respectively.The total work shifts of all types must be 15-18 work shifts.A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [39].For a more detailed description of this input instance, the reader can refer to [32].

Input Instance, BCV3-46.2
This input instance was collected from a rather small department of a real hospital, using the nurse rostering model and algorithms developed at KaHo Sint-Lieven [40].It comprises a non-cyclic problem.The number of nurses equals 46, the number of shift types equals 3 (Day, Early, Late and Night), the scheduling period is 26 days long and there is only one skill level (Nurse).A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [41].For a more detailed description of this input instance, the reader can refer to [34].

Input Instance, MUSA
In this input instance, there are 11 nurses, there is only one shift type (Day), the scheduling period is 14 days long and there is only one skill level (Nurse).A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [42].For a more detailed description of this input instance, the reader can refer to [36].

Input Instance, LLR
This input instance belongs to a class of over-constrained nurse rostering problems.The number of nurses equals 27, the number of shift types equals 3 (Morning, Afternoon and Night), the scheduling period is seven days long and there is only one skill level (Nurse).A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [43].For a more detailed description of this input instance, the reader can refer to [33].

Input Instance, BCV4-13.1
This input instance was also collected from a rather small department at a real hospital, using the nurse rostering model and algorithms developed at KaHo Sint-Lieven [40].It comprises a non-cyclic problem.The number of nurses equals 13, the number of shift types equals 4 (Day, Early, Late and Night), the scheduling period is 29 days long and there are two skill levels (Nurse and Head nurse).A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [44].For a more detailed description of this input instance, the reader can refer to [34].

Input Instance, WHPP
In this instance, the number of nurses equals 30, the number of shift types equals 3 (Day, Early and Night), the scheduling period is 14 days long and there is only one skill level (Nurse).A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [45].For a more detailed description of this input instance, the reader can refer to [37].

Input Instance, HED01
This input instance stems form the needs of an actual hospital emergency department (HED) in Spain.An actual situation set up by the HED's management is described.HED's permanent staff consists of 16 workers, while there are also four temporary workers.The number of shift types equals 5 (Weekday morning-M; Weekday afternoon-A; Weekday night-N; Weekday stand-by duty-D; and Holiday stand-by duty-H), the scheduling period is 31 days long and there are two skill levels (Permanent staff and Temporary staff).Moreover, a minimum number of doctors must be assigned to each working shift.As a general rule for working days, different members of staff will be assigned to the different existing shifts: Four members will be assigned to the morning (M) and afternoon (A) shifts; two will be assigned to the night (N) shift; and one person will be on 24 h weekday stand-by duty (D).On Saturdays, Sundays and holidays, the HED's medical staff will work on duty, whereby four members of staff will generally work for an uninterrupted 24 h holiday stand-by duty (H).A specific description of the precise definition of the evaluation function used, the hard and soft constraints and their respective weight values is given in [46].For a more detailed description of this input instance, the reader can refer to [35].

Computational Results
The proposed two-phase stochastic variable neighborhood search algorithm approach is coded in C++ and is run on Intel ® Core™ 2 Duo CPU E7500 2.93 GHz under the Windows 7 OS.The algorithm parameters' values used are the ones presented in Table 1, Section 3.1.In order to demonstrate its efficiency and very good performance, the proposed algorithm is compared with six very effective algorithms for solving the nurse rostering problem issued in the literature [32][33][34][35][36][37] in solving the same seven input instances.
In Table 2, the performance and efficiency of the proposed two-phase stochastic variable neighborhood search algorithm is shown by comparing the best timetables constructed by it with the best timetables created by the other six algorithms for each different input instance.Since in [32][33][34][35][36][37] the best rosters created are presented, we also decided in the current contribution to present and compare the best rosters constructed by the proposed approach in order to have a fair comparison between the algorithms.Note that, for each different input instance, in order to compare the roster constructed by the proposed variable neighborhood search algorithm and the roster constructed by the respective published nurse rostering algorithm, we used the same fitness function, the same hard and soft constraints and the same constraint weights as the ones used by the respective algorithms.Table 2 demonstrates that the proposed algorithm outperforms other published approaches in 6/7 cases (85.7%), considering the best roster per instance, while it achieves the same result in 1/7 cases (14.3%).From experimental results presented in Table 2, one can easily come to the conclusion that the proposed algorithm is very efficient and achieves better results compared to the other six techniques issued in the literature that have been applied to the same instances of the nurse rostering problem.
Moreover, in order to demonstrate the efficiency and very satisfactory performance of the proposed two-phase variable neighborhood search algorithm, the best rosters found by it are compared with the best-known timetables ever reported for the same seven different input nurse rostering instances [38].
Table 3 demonstrates that the proposed two-phase algorithm manages to reach the best known fitness ever reported in the literature in 6/7 cases (85.7%), while it manages to beat the best known fitness ever reported in the literature in 1/7 cases (14.3%).From experimental results presented in Table 3, one can easily come to the conclusion that the proposed algorithm is very efficient and achieves results equal to the best known ever reported for the majority of these quite different nurse rostering instances, while it manages to beat the best known ever reported result in one case.Except for that, the proposed algorithm has demonstrated experimentally that the best result ever reported for these six instances is not unique, since the best nurse rosters that the proposed algorithm created are in all six cases different from the best ones ever reported [38].This means that in these six cases, there are at least two different best ever reported rosters and that maybe the global optimum is yet to be found.The executables implementing the proposed stochastic variable neighborhood search algorithm, as well as the best rosters achieved for each input instance, can be accessed in [47].The superiority of the proposed algorithm compared to other approaches comes mainly from the fact that the algorithm succeeds in searching the search space using a new variable neighborhood search approach.In the literature, there are three commonly used swaps, that is, simple move, simple swap and Kempe swap, each one of them leading to a search algorithm to investigate a different neighborhood [23].The proposed algorithm uses only simple swaps, that is, swaps between cells without considering what these cells contain in order to perform the swap (i.e., whether they are empty or not).However, each one of the three swap procedures used by the proposed algorithm (see Sections 3.2, 3.3.1,3.3.2) applies a different swap mechanism, which leads the algorithm to search in a different neighborhood of the search space.The application of the proposed mutation operators that implement the algorithm's swapping mechanisms, applied in this specific way and order, is innovative to our knowledge and different from the swapping mechanisms already investigated in [13].The strongest point of this variable neighborhood search approach is the combination of a classic neighborhood search with a -stochastic moving segment grouping swap‖ (see Subsection 3.2).The -stochastic moving segment grouping swap‖ achieves exhaustive local search, since the segment is not stable, and as a result, it ensures intensification.On the other hand, since this swap is a stochastic one, it ensures diversification.To conclude, combining and applying these three swap procedures (see Sections 3.2, 3.3.1,3.3.2) enriches the variable neighborhood search approach, since it enhances the classic neighborhood search with a -stochastic moving segment grouping swap‖.This combination ensures both intensification and diversification of the search space.
Since the nature of the proposed two-phase algorithm is stochastic, different computational results may be obtained in different runs.So, in order to demonstrate its efficiency, in Table 4, we present not only the best, but also the worst and the average results (and the respective standard deviations-STDs), considering the fitness function value achieved and the execution time of the algorithm.Additionally, we present the respective coefficient of variation (CV) and the success rate for each input instance.All results presented in Table 4 concern the execution of the proposed algorithm to the seven aforementioned nurse rostering input instances for 100 Monte Carlo runs.Experimental results presented in Table 4 show that the average fitness reached by the proposed two-phase variable neighborhood search approach is, in most cases, very close to the best one achieved for each input instance.This demonstrates that the proposed algorithm is stable and efficient.We also notice that CV of fitness function value ranges from 0% to 41.7%, with the great majority of CV values being below 10%.More specifically, in four out of seven cases, CV equals 0, which means that the algorithm is totally homogenous.This observation leads to the conclusion that the behavior of the algorithm concerning the resulted fitness function value for 100 Monte Carlo runs per input instance is quite homogenous.Note that we have intentionally avoided calculating CV values for the execution time.This is done because STD and average values are close to each other, so the CV value for the execution time would be misleading.Moreover, in Table 4, we present the success rate, i.e., the percentage of cases that the proposed algorithm achieves the best fitness function value among 100 Monte Carlo runs.The fact that, in most cases, the success rate achieved is bigger than 60%, demonstrates the efficiency of the proposed algorithm.In addition, in four cases, the success rate is 100%.
Finally, we present convergence results (maximum, average and minimum evaluation function value versus number of function evaluations) in order to illustrate the evolutionary behavior of the proposed algorithm.Figure 8 illustrates the convergence behavior of the proposed algorithm for four input instances, namely, MUSA, LLR, BCV4-13.1 and HED01.In all cases, the algorithm's convergence behavior is very satisfactory, since it avoids falling into local optima very quickly and shows significant improvement during the evolutionary process.In these experiments, because our main concern was to reach or even beat the best ever reported result for each different input instance, we used as the determination criterion -the total number of generations‖.

Investigating the Effect of Parameter Setting to Algorithms' Performance
In this section, the effect of parameter setting to algorithms' performance is investigated and the results of different specific parameter values are demonstrated.Except for that, an indication of the contribution of each component of the algorithm is presented.More specifically, the contribution of the two-phase approach compared to a single-phase approach is investigated.As stated in Section 3.1, the first phase of the algorithm deals with the assignment of nurses to working days, while the second phase of the algorithm deals with the assignment of nurses to shift types.Therefore, the first phase of the algorithm cannot solve the nurse rostering problem alone, while the second phase can be applied to solve the nurse rostering problem as a single-phase approach.Experimental results presented below, show that using together the first and the second phase of the algorithm, as a two-phase approach, achieves better results to using only the second phase of the algorithm as a single-phase approach.
Due to the fact that there are no obvious criteria for defining specific parameter values of the proposed algorithm for all instances of the problem, we have selected these values by trial and error.More precisely, we have conducted exhaustive experiments and selected the values that achieved the best simulation results and the best algorithm's behavior.Tables 5-25, presented in the next paragraphs, show the effect of the value of first phase's number of cycles, the effect of the value of the population size and the effect of the value of swapping probability to algorithms' performance and behavior.
In Tables 5-7, experimental results that investigate the effect of parameter setting to algorithms' performance and behavior regarding input instance Valouxis-1 are presented.
As shown in Table 5, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.Except for that, if we do not use the first phase at all (first phase cycles = 0), experimental results are rather worse.As shown in Table 6, setting the value of population size equal to 2 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.As shown in Table 7, setting the value of swapping probability equal to 0.99995 assists the algorithm in achieving both the lowest best fitness value and the highest success rate.In Tables 8-10, experimental results that investigate the effect of parameter setting to algorithms' performance and behavior regarding input instance BCV3-46.2 are presented.
As shown in Table 8, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.Except for that, if we do not use the first phase at all (first phase cycles = 0), execution times are rather worse.As shown in Table 9, setting the value of population size equal to 2 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.As shown in Table 10, setting the value of swapping probability equal to 0.97 assists the algorithm in achieving the lowest best fitness value in the lowest execution time.In Tables 11-13, experimental results that investigate the effect of parameter setting to the algorithms' performance and behavior regarding input instance, MUSA, are presented.
As shown in Table 11, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.Except for that, if we do not use the first phase at all (first phase cycles = 0), execution times and success rate are rather worse.As shown in Table 13, setting the value of swapping probability equal to 0.97 assists the algorithm in achieving the lowest best and average fitness value and the highest success rate.In Tables 14-16, experimental results that investigate the effect of parameter setting to algorithms' performance and behavior regarding input instance, LLR, are presented.
As shown in Table 14, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best and average fitness values and the highest success rate.Except for that, if we do not use the first phase at all (first phase cycles = 0), execution times and success rate are rather worse.As shown in Table 16, setting the value of swapping probability equal to 0.5 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively, as well as the highest success rate.In Tables 17-19, experimental results that investigate the effect of parameter setting to algorithms' performance and behavior regarding input instance, BCV4-13.1, are presented.
As shown in Table 17, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times.For this specific input instance, if we do not use the first phase at all (first phase cycles = 0), experimental results are much less the same.In Tables 20-22, experimental results that investigate the effect of parameter setting to algorithms' performance and behavior regarding input instance, WHPP, are presented.
As shown in Table 20, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times.Except for that, if we do not use the first phase at all (first phase cycles = 0), execution times and success rate are rather worse.As shown in Table 22, setting the value of swapping probability equal to 0.45 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.In Tables 23-25, experimental results that investigate the effect of parameter setting to algorithms' performance and behavior regarding input instance, HED01, are presented.
As shown in Table 23, setting the first phase's number of cycles equal to 1 assists the algorithm in achieving the lowest best fitness value = which is the lowest fitness value ever reported in the literature, in the lowest execution time.Except for that, if we do not use the first phase at all (first phase cycles = 0), the lowest best fitness value achieved is rather worse.As shown in Table 24, setting the value of population size equal to 2 assists the algorithm in achieving the lowest best fitness value, which is the lowest fitness value ever reported in the literature.As shown in Table 25, setting the value of swapping probability equal to 0.85 assists the algorithm in achieving the lowest best fitness value, which is the lowest fitness value ever reported in the literature, in the lowest execution time.

Conclusions and Future Work
In this contribution, a generic two-phase stochastic variable neighborhood search algorithm has been designed, implemented and applied to the nurse rostering problem in order to create feasible and efficient rosters.The algorithm has been tested with seven different real-world nurse rostering instances in order to demonstrate its quality and efficiency.Computational results showed that the proposed algorithm achieves better results compared to six other very effective algorithms published in the literature that have been applied to the same nurse rostering input instances using the same evaluation criteria.Moreover, the proposed algorithm manages in one case to beat the best-known fitness achieved in the literature till now.In addition, in the other six cases, it manages to reach the best-known fitness achieved in the literature and prove experimentally that there are at least two different best ever reported rosters for these instances.Finally, the application of the proposed algorithm and its verifications to other newly published nurse rostering instances and to problems belonging to other timetabling or scheduling domains will be one of the main issues of our future work.

Figure 1 .
Figure 1.The structure of the proposed stochastic variable neighborhood search algorithm.

Figure 2 .
Figure 2. The structure of the first phase of the proposed stochastic variable neighborhood search algorithm.

Figure 5 .
Figure 5.The structure of the second phase of the proposed stochastic variable neighborhood search algorithm.

Set algorithm parameters' values Main Algorithm Initialize d individuals with working days Proceed to and conclude First Phase Copy best individual to other individuals Assign shift types randomly to any individual, driven by working days Proceed to and conclude Second Phase Near optimal solution produced End of AlgorithmTable 1 .
The parameters' values used for each input instance.

Table 3 .
Comparing the best timetables constructed by the proposed algorithm with the best timetables ever reported for these specific instances.

Table 4 .
Computational experiments demonstrating the efficiency, stability and homogeneity of the proposed algorithm.STD, standard deviation; CV, coefficient of variation.

Table 5 .
Investigating the effect of first phase's number of cycles for input instance, Valouxis-1.

Table 6 .
Investigating the effect of population size for input instance, Valouxis-1.

Table 7 .
Investigating the effect of swapping probability for input instance, Valouxis-1.

Table 8 .
Investigating the effect of first phase's number of cycles for input instance, BCV3-46.2.

Table 9 .
Investigating the effect of population size for input instance, BCV3-46.2.

Table 10 .
Investigating the effect of swapping probability for input instance, BCV3-46.2.

Table 11 .
Investigating the effect of first phase's number of cycles for input instance, MUSA.As shown in Table12, setting the value of population size equal to 2 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.

Table 12 .
Investigating the effect of population size for input instance, MUSA.

Table 13 .
Investigating the effect of swapping probability for input instance, MUSA.

Table 14 .
Investigating the effect of first phase's number of cycles for input instance, LLR.As shown in Table15, setting the value of population size equal to 2 assists the algorithm to achieve the lowest best and average fitness values and the highest success rate.

Table 15 .
Investigating the effect of population size for input instance, LLR.

Table 16 .
Investigating the effect of swapping probability for input instance, LLR.

Table 17 .
Investigating the effect of first phase's number of cycles for input instance, BCV4-13.1.As shown in Table18, setting the value of population size equal to 2 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.

Table 18 .
Investigating the effect of population size for input instance, BCV4-13.1.

Table 19 ,
setting the value of swapping probability equal to 0.97 assists the algorithm to achieve the lowest best and average fitness values in the lowest best and average execution times, respectively.

Table 19 .
Investigating the effect of swapping probability for input instance, BCV4-13.1.

Table 20 .
Investigating the effect of first phase's number of cycles for input instance, WHPP.

Average STD CV (%) Best Worst Average STD
Table 21, setting the value of population size equal to 2 assists the algorithm in achieving the lowest best and average fitness values in the lowest best and average execution times, respectively.

Table 21 .
Investigating the effect of population size for input instance, WHPP.

Table 22 .
Investigating the effect of swapping probability for input instance, WHPP.

Table 23 .
Investigating the effect of first phase's number of cycles for input instance, HED01.

Table 24 .
Investigating the effect of population size for input instance, HED01.

Table 25 .
Investigating the effect of swapping probability for input instance, HED01.