A simulation-optimization approach for the surgery scheduling problem : a case study considering stochastic surgical times

Article history: Received September 2


Introduction
Hospital management has changed from a resources optimization approach to a balance between efficiency and quality (Brailsford & Vissers, 2011).Consequently, the decision-making process has become increasingly complex (Hulshof et al., 2012) and the use of operations research has grown steadily in recent years (Dobrzykowski et al., 2014).The hypothesis is that the appropriate allocation of resources, in support processes, positively influences the quality of service (Van De Klundert et al., 2008).It is also expected that a proper scheduling of resources involved with the patient (human, equipment and facilities) leads to the improvement of critical variables such as cost, coverage, and response times (Cardoen et al., 2010;Elkhuizen et al., 2007;Liu et al., 2011;Velasco et al., 2012).Despite the growth of research in this area, there is still a gap regarding the actual practices of hospitals (Brailsford & Vissers, 2011;van Sambeek et al., 2011).One of the main challenges for the research community is to achieve collaborative spaces that allow a greater impact of developments in real systems (Pitt, Monks et al., 2015;Velasco et al., 2012).
In this context, the scheduling of operating rooms is particularly relevant.Surgical services account for 40% of total costs of hospitals and are an important part of their income (Denton et al., 2007).Consequently, this issue has been addressed from different approaches (Cardoen et al., 2010;Hasvold & Scholl, 2011).In particular, the surgery scheduling problem (SSP) aims to define the time and room where a set of surgeries must be performed (Abdelrasol et al., 2014).Several authors have demonstrated the positive impact of the use of formal methodologies and tools to make this decision (Díaz-López et al., 2015;Ghazalbash et al., 2012;Litvak, 2010;Shylo et al., 2013).However, in a regional context, only a few approaches have been made regarding this problem (Alonso et al., 2013;Ceballos-Acevedo et al., 2014;Díaz-López et al., 2015;Velásquez-Restrepo et al., 2012).Moreover, Velasco et al. (2012) concluded that there was a field of interest for the academic community in the study of problems in real environments.
Regarding this, Abdelrasol et al. (2014) reviewed the literature and identified two lines of work that will be of interest to hospitals and researchers in the coming years.The authors conclude that the development of efficient algorithms and the inclusion of the variability of the surgical times should be addressed in order to close the gap between management practices and research.Assuming a deterministic surgical time, as is usual in research, makes the implementation of the results difficult (Velásquez-Restrepo et al., 2012).The same challenge has been addressed in other classical problems (Holm et al. 2013;Juan et al., 2014;Sarin et al., 2013).Consequently, the use of hybrid simulation-optimization techniques has flourished (Juan et al., 2015).
Due to the complexity of relations to be modeled, it is common to use discrete event simulation for operating rooms scheduling (Aringhieri et al., 2015;Baesler et al., 2015;M'Hallah & Al-Roomi, 2014).However, optimization could overcome the limitations of simulation when selecting the evaluated scenarios (Zhang & Xie, 2015).According to Baesler et al. (2015), most of the work in this area is concentrated in the simulation phase and uses dispatching rules for the optimization problem.This situation implies an enhancement field since, in scheduling problems, the use of more elaborate optimization techniques has shown a positive impact on the quality of the solutions (Molina-Sánchez & González-Neira, 2016).In fact, after a revision of 216 papers in literature, published between 2000 to 2014, Samudra et al. (2016) suggest the usage of simulation-optimization methods to solve the complex surgery scheduling problems.
Recently, some studies have explored this approach using search techniques for the optimization phase.For example, Testi et al. (2007) proposed a hierarchical approach of three stages weekly scheduling.The objectives were improve throughput, overtime and reduce waiting list.In the first phase the number of sessions to assign to each specialty is solved ass a bin packing problem.The second phase consists in the assignment of specialties to the operating rooms.And, in the final step, the simulation procedure is implemented to define the best sequence.Saadouli et al. (2015) studied the problem of scheduling surgeries within a planning horizon of one day.The authors solve the optimization problem using exact methods and discrete events simulation is used to evaluate the solution quality.Baesler et al. (2015) use simulated annealing at the optimization phase.The authors solve instances of up to 170 patients, and found that the use of more sophisticated search techniques, other than the dispatching rules, can deliver up to 18% of improvement in utilization levels.Wang et al. (2014) group surgeries to schedule each kind of them in a different operating room to reduce the flow variability and waiting time.They present an alternate simulation-optimization structure where the results of simulation are the inputs of the optimization model and its results are used for the next simulation, and so on, until meeting the stopping criteria.With this approach authors reduced the waiting time in 89% and the overtime in 46%.Beaulieu et al. (2012) proposed a methodology with four phases to deal with surgery scheduling and rescheduling.In the first phase authors proposed a MILP model for the assignment of the surgeries to a given day.In the second step, diverse strategies are implemented to schedule the surgeries in that day.These assignments are evaluated in the third stage with the use of simulation.Then, in the fourth phase decisions of rescheduling are taken if it is required.In this work, a four steps approach is proposed.First, surgery cases are assigned to a given day.Second, they are scheduled according to different strategies, which are evaluated through a simulation tool in the third step.If needed, feedback and rescheduling occur in the fourth step.Authors evaluated the performance of their approach with randomly generated instances based on real data of a large hospital in Montreal.Duma and Aringhieri (2015) designed an online solution method to adjust a patient schedule that have been affected by a delay during its execution.Possible modifications are rescheduling the surgery or assignment of overtime for the surgery.The authors also developed an offline simulation-optimization model to evaluate the impact of the online procedure.Aringhieri and Duma (2015) developed a simulation-optimization model to Schedule surgeries.The method includes a greedy algorithm to construct the solution and a local search that improve it.Authors center their attention on the analysis of a surgical clinical path of a patient to optimize the use of resources and assess the effect of the optimization in different indexes focused both on patients and facility.Molina-Pariente et al. (2016) consider stochastic surgery times and uncertain arrivals of surgeons and emergency surgeries to minimize under and overtime costs of operating the rooms and costs of exceed the capacity constraints.Authors proposed an iterated greedy local search hybridized with a Monte Carlo simulation and compare the results obtained with this method with the deterministic solution of the problem, demonstrating the cost reduction when solving the stochastic problem.Landa et al. (2016) solve the advance and allocation scheduling problems.In the first one authors assigned a date for surgery and an operating room block for the surgeries.They determined the sequence of the surgeries in each operating room in a specific day.They developed a hybrid algorithm that combines Monte Carlo simulation with neighborhood search methods to obtain solutions that maximizes the utilization rate while patient cancellations are reduced.Addis et al. (2016) schedule also elective procedures using the block scheduling strategy and considering stochastic surgery times and stochastic arrivals of new patients.Therefore, authors propose a rolling horizon method for patients selection and assignment where an integer linear programming model (ILP) is used iteratively, minimizing waiting time and tardiness of assignment of patients, the uncertainties that are generated between one period to another are recovered at the next iteration of the ILP model.Finally, Ozcan et al. (2017) applied a simulation-optimization approach, in which simulated annealing metaheuristic deals with the optimization issue, to allocate resources for thyroid surgical treatment without ignoring the other types of surgeries.The implemented approach allows managers to involve hospital and patients objectives by obtaining the opportunity costs of prioritizing one objective to another.
In this context, it is necessary to continue exploring other optimization techniques to improve the quality of the solutions and reduce computational time-a factor that can hamper implementation of sophisticated optimization techniques (Baesler et al., 2015).Additionally, research is focused on optimizing objective functions related to efficiency.However, in environments of high variability, it is interesting to quantify the impact of high rates of utilization in the service level (Holm et al., 2013).Consequently, this paper proposes a simulation-optimization technique to solve the stochastic version of the surgery scheduling problem.
Instead of generating one solution for each instance, the proposed technique generates a set of solutions so that the hospital can decide which one to implement.In general, the proposed technique implies that the surgical time can be modeled by probability density functions.Then, these functions are used to generate a set of realizations corresponding to different percentiles of surgical time.Subsequently, for each percentile, the deterministic version of the problem is solved using GRASP.Thus, there are as many solutions for one instance as test percentiles chosen.Finally, through Monte Carlo simulation, the confidence intervals are calculated for three indicators: i) utilization percentage of rooms, ii) surgeries percentage that are delayed, and iii) average delay time in surgery.In this manner, the proposed technique quantifies the impact of high levels of utilization in the service level.

Methods
The problem studied is the SSP for elective procedures.This is due to some hospitals has separate rooms for elective surgeries than those for emergency and urgency surgical procedures.The set of surgeries to be performed and the availability of a set of homogeneous operating rooms are known.In this context, two decisions must be made: i) the allocation of the procedures to the rooms, and ii) the starting time of each procedure.The duration of surgeries is not known with certainty; however, it can be modeled by using a random variable.The aim is to maximize the occupation of operating rooms.Likewise, the impact of the chosen strategy by means of two indicators of service quality must be quantified: i) the percentage of surgeries that had delays, and ii) the average delay time of scheduled surgeries.
A methodology that combines simulation and optimization techniques to attack this problem was designed.Figueira and Almada-Lobo ( 2014) have proposed a classification of hybridization approaches of these two methods, depending on their role in solving the problem.Specifically, from the point of view of hierarchical structure, the simulation-optimization approach selected for this research is sequential simulation-optimization (SSO) and in particular the evaluation procedure named Solution Completion by Simulation (SCS).The SSO consists in that optimization and simulation modules run sequentially one after another.The SCS allows to complete the initial solution obtained in the optimization phase by giving more accurate values for the different variables and objective function through simulation (Figueira & Almada-Lobo, 2014).Thence, we resolved the deterministic SSP at the optimization phase using GRASP metaheuristic due to the NP-hardness of the problem (Brucker, 2007;Pinedo, 2012).Posteriorly, we applied a Monte Carlo simulation to validate the results and calculate the expected value of the indicators.A graphical representation of the methodology used in this study, based on the proposed algorithm in (Juan et al., 2015), is shown in Fig. 1.Initially, it was assumed that all procedures would last as many minutes as the 50th percentile of its probability density function.Taking the durations as known, with the value of 50 th percentile, GRASP was used to schedule the operating rooms.The resulting schedule was analyzed by using a Monte Carlo simulation model, taking the corresponding probability distribution for each surgery.It allows us to estimate the expected values of utilization and quality of service indicators.Finally, this process was repeated by increasing 5% in the percentile used for surgeries' durations up to the 95th percentile.Thus, a total of ten different scheduling scenarios were considered.The details of the phases of optimization and simulation will be discussed next.

Optimization phase
For this stage we firstly proposed the mathematical formulation of the problem and secondly solved it with a simulation-optimization procedure in order to evaluate the performance of the solution method under uncertainty.

Mathematical formulation
The mixed integer linear programming MILP proposed model is: , , , : Binary variable that takes value 1 if surgery is sequenced immediately after surgery on day at room : continuous variable for the finishing time of surgery Objective function:

Sets
, , , Eq. ( 1) calculates the utilization rate in hours.Constraints set (2) ensures that in each room at each day in each position of the sequence at most only one surgery is scheduled.Constraints set (3) guarantees that every surgery is scheduled at most one time in the scheduling horizon.Constraints set (4) makes that at each room on each day the surgeries be sequenced one next to another one.Constraints set (5) controls that surgeries do not overlap and calculates the finalization time of every scheduled surgery and together with constraint set (6) allows to calculate the completion time of a scheduled surgery.Constraints set (7) states that the completion time of a surgery must not exceed the maximum scheduling time and sets in zero the completion time of an unscheduled surgery.Constraints set (9) define the domain of the binary decision variables.

Simulation-optimization procedure
As it was mentioned in Methods section GRASP metaheuristic was selected for the optimization phase of the simulation-optimization procedure.GRASP is an iterative multi-start metaheuristic, designed for solving combinatorial problems.Each iteration of the metaheuristic consists of two steps: construction and local search (Resende & Ribeiro, 2010).For construction, the surgeries were organized according to a greedy function.Subsequently, a list of candidates (i.e.surgeries that can be scheduled) was constructed, and one of them was randomly selected.This process was repeated until all the surgeries were scheduled.Later in the local search phase, the initial solution is taken and all exchanges of surgeries among different days and rooms were evaluated.Finally, it was assessed whether, after all exchanges of local search, it was possible to schedule some of the surgeries that had not initially been sequenced in the construction phase.Fig. 2 presents a graphical representation of the GRASP algorithm proposed.To conclude on quality, the metaheuristic's solutions were compared with exact solutions of an equivalent linear programming model.

Simulation phase
In problems where the parameters have random behavior, it is expected that the performance measures have it also (Li et al., 2012).Consequently, the Monte Carlo simulation was used to validate the results and to estimate the risk (Bayer, 2014;Tako et al., 2014).To calculate the required number of iterations of the simulation, the proposed algorithm in Framinan and Perez-Gonzalez (2015) was used.An amplitude interval for indicators of 2.5% around the mean was established as a convergence threshold to be reached.Thus, the simulation stopped when the indicator with the highest variance reached such accuracy.Three indicators were calculated: occupation rate ( 14), waiting time (15) and waiting rate ( 16).The occupation rate shows the proportion of available time that the rooms are used for surgeries.The waiting time represents the average time in minutes that a patient has to wait if his corresponding surgery is starting late.The waiting time is the proportion of surgeries that begin after its schedule hour.

Case Study
The proposed methodology was tested to schedule surgeries for one week at the Méderi University Hospital in Bogotá, Colombia.The surgery service at the hospital has 19 operating rooms and performs approximately 35,000 procedures a year.According to a previous study, the use of formal tools for scheduling rooms has high potential for impact at this hospital (Estupiñán et al., 2016).The authors found that by shifting from manual to algorithmic programming with the use of dispatching rules, the occupation indicator can be improved considerably.These results confirm the hypothesis that a lowquality scheduling can explain poor performances in service indicators (Litvak, 2010;van Sambeek et al., 2011;Velasco et al., 2012).
Currently, the scheduling process is made manually and there are no specialty blocks, this implies an open surgery scheduling.Historical data corresponding to all surgeries performed during May 2014 were analyzed.According to the method proposed by Baesler et al. (2015), by using tests of homogeneity and goodness of fit, the probability density function that best describes the duration of surgeries was determined.In total, 26 probability density functions for surgeries were generated, which explain 89% of the occupation of operating rooms.The goodness of fit tests was run on Expertfit (Flexim 7).
In order to parametrize GRASP, random instances were generated for the parameterization of the metaheuristic through a 3 factorial design.Factors were greedy parameter of GRASP (with levels 0.2, 0.25 and 0.3), way in surgeries are sorted (ascending order of duration, ascending order of variability, descending order of duration) and number of iterations of GRASP (100, 500 and 900).According to ANOVA, the interaction of the three factors is statistically significant (p-value < 0.0001 for a 95% confidence).Additionally, a Tukey test was performed-to determine the best treatment.It was found that the best combination of parameters is as follows: sorting surgeries in ascending order according to their duration, selecting the 20% of non-scheduled procedures for the list of candidates, and running 900 iterations.Assumptions of normality (p -value> 0.05 in the Kolmogorov-Smirnov test) and homogeneity of variance (p-value = 0.056 for the Levene test) were tested.
Expected results of the case study are divided into two categories: i) The performance GRASP-based approach results were measured with the outcomes obtained by a mixed integer linear programming model that was run for six hours.ii) The real instance obtained with the hospital database was run with the simulationoptimization approach to establish utilization and service levels for ten different scenarios.

GRASP performance in comparison with MILP model
A first aspect of interest is the comparison of GRASP against the results obtained with an integer programming model.The problem of scheduling operating rooms is highly related to the scheduling problem for parallel machines; therefore, it can be formulated in that way (Pinedo, 2012).However, since the objective of the mathematical model is to determine the quality of the solutions given by the metaheuristic, a simplification of the formulation through a multiple knapsack problem approach is proposed (Chen et al., 2016).
To compare the results, ten deterministic small instances were generated, one for each percentile of surgical time, by using the probability distributions previously fitted with data of original surgeries at the Hospital.The random generated instances consisted in 27 types of surgeries, 10 operating rooms and a total of 63 surgeries to be scheduled in only one day.Those instances were tested with the MILP model and with GRASP metaheuristic.To define these ten instances, many tests were carried out with instances of different sizes.The process of creating them begun with big instances similar to the Hospital ones, and we were reducing the size of them by decrementing the number of surgeries, days and rooms until we achieve to solved it in GAMS.With big and small instances, we obtained always an out of memory response.The reason for creating only small instances is that the big and medium instances designed, were tested with GAMS and the response was always out of memory.GRASP was executed 10 times (with 10 different seeds) for each Scenario with the parameters selected in section 3, and the outcomes of worst, best and average results of occupation rate and execution time were recorded (Tables Table 1  and Table 2 respectively).The model was solved using the solver CPLEX 12.7 of GAMS and the GRASP approach was coded using Visual Basic for Applications of Excel.It was found that the use of the algorithm involves an average reduction of 4.26% in the utilization compared to the solutions found with the model.However, the computation time savings is on average around 93%. Furthermore, because the comparison is being done with a simplified model, this underestimates the computational saving of the algorithm.The result is important for two reasons.On one hand, although there is recent evidence of good performance using GRASP in other operating room scheduling problems (Cartes Rubilar & Medina Duran, 2016), to the best of our knowledge there are no previous studies that test the performance of GRASP in the surgery scheduling problem.The reduction in the found objective function is good, compared to what was reported by other authors in production scheduling problems (Molina-Sánchez & González-Neira, 2016).On the other hand, to solve stochastic problems it is vital to find efficient algorithms from a computational time point of view (Figueira & Almada-Lobo, 2014).Since deterministic solutions present difficulties in implementation, in real systems, it is often preferable to find good solutions in reasonable computational times (Sarin et al., 2013).

Solving hospital instance
A hospital instance, obtained from the May 2014 records, has been analyzed.Fig. 3 shows the solution frontier of the expected value of the three indicators.On the left side, occupation and the percentage of patients waiting are compared.On the right side, occupation and the average waiting time for patients are plotted.Additionally, the figure shows the performance of the deterministic solution proposed by Estupiñán et al. (2016) who used the same dataset of surgeries, allowing us to make a straightforward comparison.It can be seen that, compared with the result of a deterministic approach, the proposed methodology is a solution that improves the three indicators.Also, this border can be used to support the decision-making process.

Fig. 3. Pareto frontier between waiting time, percentage of waiting patients and occupation rate
On the other hand, Fig. 4 compares the confidence intervals of the solutions obtained with the deterministic approach proposed by Estupiñán et al. (2016) and those obtained with the simulationoptimization methodology.It is important to note that the deterministic approach projected an 80% of occupancy.However, with a confidence level of 95%, this indicator will be between 75.6% and 76.3%.This is because the authors did not take into account the variability of surgical times.considering the resource constraints.On the one hand, a deviation of 5% compared to the projected solution may be problematic in some systems.Hence, the variability of surgical times should be included.On the other hand, setting a target for the occupancy level without considering the reductions in the service might lead to patient dissatisfaction.In this context, the proposed methodology delivers a set of solutions to the hospital.From this set, the one that is more aligned with hospital objectives will be selected.
In general, the behavior of the indicators is as expected: high levels of utilization involving low levels of service.What is important, however, is to quantify the reductions to be made in an objective to enhance the other, as their relationships are not linear.Thus, for example, if the hospital is interested in having an average occupancy of 80%, it must be willing to let 44% of their patients wait for surgery delays and for the average waiting time to be 85 minutes.Table 3 presents the confidence intervals for the percentage of patients that have to wait, the occupation rate and the average time (in minutes) that a patient have to wait if the waiting occurs, all of them depending on the percentile used for the duration of the surgeries.
It can be seen that for percentage of patients that have to wait and occupation rate indicators none of the confidence intervals overlap.That is there here are statistically significant differences in these indicators as the percentile used for surgeries duration changes.The average waiting time for patients is the opposite.It cannot be concluded that the indicator to weaken with the increase in the occupancy.The reduction made by the increase in occupation is mainly associated with having more patients waiting.

Conclusions
This work presents a simulation-optimization approach to room surgery scheduling with real data from an important hospital in Bogota, Colombia.The Monte Carlo simulation dealt with stochastic duration times of surgeries and GRASP metaheursitic for optimization purposes.Real data of surgeries' duration times of a normal month was analyzed to determine its probability distribution and simulate real durations in the proposed scheduling method.
Results show that there is an inversely proportional relation between the service level and utilization rate.While utilization decreases, the percentage of people that must wait increases.In fact, reductions in utilization levels can be clustered in three groups depending on the level of service.The sacrifice, in terms of utilization, needed to decrease delay times is explained as follows.The percentage of patients who will not have to wait for a delay in starting their surgery is measured.The reductions of 5% in the percentage of patients that have to wait correspond to increases of the same size in the percentile used for the duration of the procedure.The sacrifices in occupation rate can be classified into three groups.
The first group (circle points in Fig. 4) corresponds to increments in the percentage of delayed surgeries between 34% and 49%, which represents reductions of about less than 3% of occupation (2.52% on average).In the second group, triangle points on the graph, there are increases of up to a 4% of the utilization level, implying an increment in the percentage of patients that have to wait between 9% and 34%.Finally, the third group is the scenario where, in order to ensure that less than 5% of patients must wait, a reduction in the utilization is around 9% (see the squares points on graph).Thus, it implies that there is no one unique optimal solution, but there is a wide range of possibilities that the institution can analyze to choose the one that adjusts better to its policies.The simulation-optimization approach obtained better results than the implementation of dispatching rules.
Future work can incorporate more constraints to the analyzed problem, such as the personnel and instrumental availability.Also, the blocks of time in operating rooms can be added.In terms of the surgical procedures scheduled, emergency and urgency surgeries can be included.Finally, other simulation-optimization approaches should be implemented not only for validation of the solution but for the construction of the schedule.

:
Position of the sequence in which a surgery is scheduled in a specific room and day Parameters: : duration [hours] of surgery : maximum scheduling time [hours] for a day : big positive number Decision variables:

Fig. 4 .
Fig. 4. Confidence intervals of 95% for the indicators using the Tukey methodCompared to the deterministic approach, the proposed methodology represents two advantages: i) it allows managing the risk associated with the indicators estimation, and ii) it allows setting realistic targets

Table 1
Occupation rate of MILP and GRASP, and GAP of GRASP for each scenario

Table 2
Computational time of MILP and GRASP for each scenario

Table 3
Confidence intervals (95%) for the indicators