(PDF) Operating room scheduling problem under uncertainty: Application of continuous phase-type distributions

IISE Transactions ISSN: 2472-5854 (Print) 2472-5862 (Online) Journal homepage: https://www.tandfonline.com/loi/uiie21 Operating Room Scheduling Problem under Uncertainty: Application of Continuous Phase-Type Distributions Mohsen Varmazyar, Raha Akhavan-Tabatabaei, Nasser Salmasi & Mohammad Modarres To cite this article: Mohsen Varmazyar, Raha Akhavan-Tabatabaei, Nasser Salmasi & Mohammad Modarres (2019): Operating Room Scheduling Problem under Uncertainty: Application of Continuous Phase-Type Distributions, IISE Transactions To link to this article: https://doi.org/10.1080/24725854.2019.1628372 Accepted author version posted online: 10 Jun 2019. Submit your article to this journal View Crossmark data Full Terms & Conditions of access and use can be found at https://www.tandfonline.com/action/journalInformation?journalCode=uiie21 Operating Room Scheduling Problem under Uncertainty: Application of Continuous Phase- Type Distributions Mohsen Varmazyara, Raha Akhavan-Tabatabaeib, Nasser Salmasic*, and Mohammad Modarresa a Department of Industrial Engineering, Shairf University of Technology, Tehran, Iran b School of Management, Sabanci University, Istanbul, Turkey c Corning Incorporated, Wilmington, NC, USA t Corresponding author Nasser Salmasi

ip Authors’ Biography Mohsen Varmazyar has received his Ph.D. in industrial engineering from cr Sharif University of Technology (SUT) in IRAN. He holds Master degree in industrial engineering at SUT. His research interests focus on applied operations research, scheduling and stochastic processes. us an Raha Akhavan-Tabatabaei is associate professor of Operations Management and co-director of Masters in Business Analytics at Sabanci School of Management. Prior to this position, she was associate professor of Industrial M Engineering and the founding director of Masters in Analytics at Universidad de los Andes in Bogota, Colombia. Before that, she worked as senior ed industrial engineer at Intel Corporation in Arizona, USA. She has received her Ph.D. and M.Sc. in Industrial Engineering and Operations Research from North Carolina State University, and her B.Sc. from Sharif University of pt Technology. Her research is focused on stochastic modeling and data-driven decision making with applications in healthcare, logistics, revenue ce management and reliability. Nasser Salmasi has received his PhD in the area of Industrial Engineering Ac from Oregon State University, USA. He was working at the Department of Industrial Engineering at Sharif University of Technology, Tehran, Iran as an Associate Professor for nine years. Since October 2015 he has joined Corning Incorporated in USA as an operations research analyst. His primary areas of research interests are applied operations research, sequencing and scheduling, and simulation. 1 Mohammad Modarres is a professor at Dept. of Industrial Engineering, Sharif University of Technology, Iran. He received his PhD in Systems Engineering and Operations Research from University of California, Los Angles (UCLA) in 1975. His research interests are operations research, revenue management and robust optimization. Abstract t This paper studies the stochastic operating room (OR) scheduling problem integrated by Post-Anesthesia Care ip Unit (PACU), named the operating theater room (OTR) problem. Due to the inherent uncertainty in surgery cr duration and its consecutive PACU time, the completion time of a patient should be modeled as the sum of a us number of random variables. Some researchers have proposed the use of normal distribution for its well-known additive property, but there are questions regarding its fitting adequacy to real OTR data, which tends to be an asymmetric with a long tail. We propose to estimate the surgery and PACU times with the family of continuous phase-type (CPH) distributions, which provides both fitting adequacy and additive property. We first compute M the completion time of each patient analytically and compare the results with normal and lognormal distributions on a series of real OTR datasets. Then, we develop a search algorithm embedding constructive ed heuristic and meta-heuristic algorithms as a sequence generator engine for the patients, and apply the CPH distribution as a chance constraint to eventually find the schedule of each sequence in OTR problems. The best pt algorithm among several tested constructive heuristic algorithms is used as the neighborhood structure of meta- heuristic algorithms. We finally construct a numerical example of OTR problem to illustrate the application of ce the proposed algorithm. Ac Keywords: operating theater room scheduling; Post-Anesthesia Care Unit (PACU); minimization of makespan; continuous phase-type (CPH) distributions; heuristics and meta- heuristic algorithms. 1 Introduction Healthcare expenditures have been increasing noticeably over the last couple of decades due to the growth and aging of the population (Saadouli et al. 2015). Thus, hospital managers have been 2 interested in improving the efficiency and effectiveness of healthcare systems by developing modern planning and scheduling tools for their resources. The operating room (OR) is one of the most critical resources, and is considered as the engine of hospitals since it is the largest cost as well as the largest revenue center (Zhao and Li, 2014). However, research indicates that in some situations, hospital managers have not reached their target utilization levels for these highly valuable resources (Erdogan and Denton, 2011; Zhao and Li, 2014). Hence, it is essential to improve the patient flow and increase the efficiency of ORs in order to provide timely treatments for the patients and to maximize the utilization of available OR resources (Abdelrasol et al. 2013; Lee and Yih, 2014; Zhao and Li, 2014). As a required step, OR sequencing and scheduling plays a crucial role in the decision making of t ip hospital managers. cr OR sequencing and scheduling determines the start time of each surgery to be performed in different surgical groups, as well as the resources assigned to each surgery over a scheduled period. The overall us surgery process contains three stages in a predetermined sequence, i.e., the pre-operative stage, the an peri-operative and the post-operative stage (Pham and Klinkert, 2008). In the pre-operative stage, patients are transported from nursing units to the preoperative holding unit (PHU) and then moved to M an OR. In peri-operative stage, the patient is anaesthetized for surgery by an anesthetist and then operated on by one or several surgeon(s). In the post-operative stage, the patient may be transported to ed several different destinations. Most patients are taken to the postanaesthesia care unit (PACU) where pt they recover from residual effects of anesthesia under the care of PACU nurses. The problem of OR scheduling including the PACU is commonly referred to as the operating theater room (OTR) ce problem (Guerriero and Guido, 2011; Wang et al. 2015). Ac Due to the complexity of OTR problem, so far no exact method has been introduced to tackle it within a reasonable amount of time (Denton et al. 2007; Mancilla and Storer, 2012). The complexity of OTR problem results from multiple stages, coordination of many resources, and uncertainty inherent to surgical services (Cardoen et al. 2010b). The OTR problem with multiple stages and one resource is treated as a hybrid flowshop problem which is NP-hard (Dekhici and Belkadi, 2014; Wang et al. 2015; Xiang et al. 2015). Uncertainty in surgery duration is because of its variance and skewness and its dependence on other sources of variability such as gender, age, and the experience of the surgeon 3 (da Silva Godinho, 2014; Stepaniak et al. 2010). The PACU time could also be uncertain because the type of anesthesia used in the surgery or other issues such as patient conditions or age may prolong the duration of this stage (Baesler et al. 2015). Considering uncertainty in surgical services, the patients’ OTR completion time is uncertain which in turn renders the starting time of the following patients uncertain as well. Accurate estimation of completion time for each patient at each stage is critical in constructing an efficient OTR scheduling algorithm (Samudra et al. 2016; Shylo et al. 2012). The OTR completion time of patients requires the computation of the convolution or minimum/maximum of several random variables that represent surgery and PACU times. t ip On the other hand, the number of hours that ORs will be operational on each day is a strategic cr decision in hospital management, and hence developing a schedule to ensure the timely completion of surgery for the last patient is of high importance (Tànfani and Testi, 2012). Therefore, it is essential to us compute the probability distribution of the OTR completion time during the available hours of the an ORs. Depending on the probability distributions assumed for the OTR random variables, computing convolutions, minimum/maximum, and the probability of timely completion time could become a M complex operation (Samudra et al. 2016). Thus, choosing a distribution to accurately model the surgery, PACU, and the convolution of these stochastic times (i.e., the total OTR time) is a major ed challenge in the OTR scheduling problem. pt The OTR scheduling problem has been intensively studied in the literature from the decision scientists’ perspective (Abe et al. 2016; Cardoen et al. 2010b; Ferrand et al. 2014; Guerriero and ce Guido, 2011). Dexter et al. (1999) used on-line and off-line bin-packing techniques to plan elective Ac patients and evaluated their performances using the lognormal distribution in a simulation model. Lamiri et al. (2007) addressed the elective surgery planning under uncertain surgery times, estimated by a lognormal distribution. They first modeled the planning problem as a stochastic integer program and then used Monte Carlo simulation to approximate the problem by a mixed- integer program. They assigned elective patients to a particular OR in a particular period by a column-oriented formulation and solved the linear relaxation of the latter formulation via column generation. The duration of surgeries on the basis of a number of different hazard models was predicted by Joustra et al. (2013) to 4 plan the surgeries in the Academic Medical Center. They used several hazard models to compare the results with the predictions provided by surgeons. Bruni et al. (2015) suggested a comprehensive stochastic programming modeling framework, which handles the inherent uncertainty characterizing the arrival of emergency patients and the duration of surgery. In their research, surgery durations were sampled from lognormal distributions. The reviewed literature shows that fitting with popular distributions, especially normal and lognormal, has been extensively studied in OTR problem. Despite the widespread use of normal and lognormal distributions, they have several shortcomings when applied to OTR scheduling problems. For instance, the normal distribution could generate negative values, and symmetry renders it unrealistic to fit surgical durations. Also, in order to t ip efficiently schedule the OTR, computing the convolution of several random variables that represent cr the duration of a series of activities is often required. Although computing these convolutions can be a straightforward procedure for the normal distribution, it is a complex operation for the lognormal us distribution (Hans et al. 2008). Moreover, computation of objective functions such as minimization of an makespan requires the calculation of the minimum of the maximum completion times of all the surgeries. For these reasons, it is difficult to work with distributions that do not admit a closed form M function for the distribution of the convolution or minimum/maximum, such as the normal and lognormal distributions. ed Modeling stochastic surgery and PACU times with the aforementioned limitations constitutes several pt challenges when tackling the OTR scheduling problem:  Developing a general solution strategy used in different real-life situations is hampered by the ce lack of a generally acceptable distribution for modeling surgery and PACU times.  Ac Finding a single distribution that can accurately fit all the surgery and PACU times in an instance may be impossible.  Computing the convolutions of the chosen distributions in general is not a trivial task, especially if the surgery and PACU times follow different families of distributions. To tackle these challenges, we propose a new approach using continuous phase-type (CPH) distributions to model the stochastic surgery and PACU times and apply it in the OTR scheduling 5 problem. By using this family of distributions, we can closely approximate any positive and continuous distribution with good precision and compute their convolutions in an exact manner. In order to show the application of CPH distribution in OTR problem, we define our research problem as el pat; time; room;int( PACU ) stoch single ; Cmax based on scheduling notations used for OTR problems introduced by (Cardoen et al. 2010a). The reason for choosing this notation and the assumptions of this research are the following:  The problem only deals with elective patients for whom the surgery can be well planned in advance (el). t  The decisions apply to the patient level since they are made for individual patients (pat); ip  For each patient selected to be scheduled, decisions related to the assignment of a date, a time, cr and an OR have to be made using an open scheduling strategy (time, room).   us The surgical suite is integrated by post anesthesia care unit (int(PACU)). The duration of the surgeries and PACU times are addressed in a stochastic fashion (stoch). an  The objective function considered in this research is a single criterion (single), as M minimization of makespan ( Cmax ), which implies high resource utilization. The rest of the paper is organized as follows. In Section ‎2, we discuss the stochastic model of surgery ed and PACU times using CPH distributions and compare our model with normal and lognormal distributions. In Section ‎3, we propose a search algorithm by defining sequence generator engines and pt a schedule evaluator and combine the CPH distribution with heuristics algorithms to find better ce solutions. We provide a small example and show the application of CPH distribution on a small-size OTR problem. In Sections ‎4 and ‎5, we present several constructive heuristics and meta-heuristic Ac algorithms and discuss the computational results over a set of test problems. In Section ‎6, we provide the conclusions and show directions for future research. 2 Modeling stochastic times with CPH distributions To the best of our knowledge, there is no available research in the literature on the application of CPH family of distributions in the context of OTR problems. However, the CPH distributions’ density and 6 closure properties make them an attractive choice for modeling surgery and PACU in real life scenarios. In this section, we first introduce the CPH distributions and their distinguishing properties. Then, we apply their properties to fit real sample datasets of surgery times and compare the results of CPH fit with normal and lognormal distributions in terms of fitting quality. Afterwards, we consider the simplest form of OTR problem with one patient and model OTR times for this patient with CPH, normal, and lognormal distributions. 2.1 Definition and properties of the CPH distributions t CPH distributions have been introduced and formalized by Neuts (1975, 1981), as the distribution of ip time until absorption in a discrete-state continuous-time Markov chain (CTMC) with n transient states cr and one absorbing state. We assume that { X (n)}n0 denote the CTMC with finite state space us S  {0,1, 2,..., n} , where the absorbing state is numbered 0 and the transient states are numbered 1, 2, …, n. CPH distribution is shown by Z ~ PH C (, T) and the infinitesimal generator is: an T t  Q , M (1)  0 0 where T is a square matrix of dimension n, t is a column vector and 0 is a row vector of dimension n. ed Based on the infinitesimal generator matrix (Q), we have that Tij  0, Tii  0 for i  j and ti  0 and T1  t  0 where 1 is the column vector of 1s of the appropriate dimension n. The probability pt distribution of the initial states is denoted with the row vector (   0 ) and  0  1  1 . The cumulative ce distribution function, the probability density function and the ith moments of CPH distribution are Ac presented in Appendix A. Figure 1 gives matricial and graphical representations of several popular classes of CPH distributions, including the exponential, Erlang, Hyper-exponential, and Hyper- Erlang families. More detailed explanation on the representation of CPH distributions is presented in Buchholz et al. (2014). CPH distributions are dense in the set of continuous density function with support on [0, ∞). This means that there is a CPH distribution arbitrarily close to any positive distribution. Such a CPH 7 distribution is found by a process known as distribution fitting, for which there exist different efficient algorithms, like those proposed by Bobbio et al. (2005) and Thümmler et al. (2006). Some properties of CPH distributions make them flexible and tractable (Neuts, 1981). Firstly, a special form of the distribution function, the expected value, and higher-order moments for a CPH distribution can be represented in closed formulae in terms of  and T . Secondly, the CPH distributions are closed under convolution, meaning that the convolution of CPH distributions is also a CPH distribution. The initial vector and generator matrix are easily calculated from the original parameters of the distributions. Formal statements of these properties can be found in Appendix A. Based on density and closure properties of CPH distributions, the modeling of surgery and PACU t ip times by this family of distributions is an attractive choice. On one hand, surgery and PACU times can cr be modeled theoretically as accurately as needed using CPH distributions, since these times are always positive. On the other hand, regardless of the CPH distribution applied for each surgery or PACU time, the convolution is calculated by the same formula, as the CPH family is closed under us an convolution (see Appendix A). CPH distributions are thus able to accurately model different surgery or PACU time while allowing the computation of the required convolutions in an efficient manner. M Erlang Erlang Distribution Distribution with with Parameters Parameters n, and λλ n, and Exponential Exponential Distribution Distribution with Parameter λλ with Parameter    0 0 0  0  ed  0 λ     λ λ   λ λ 1 1 0 Q ; 1 1 2 ... ... n 0 T   0 0   0 0     0 0 0   pt Hyper-Exponential Hyper-Exponential Distribution Distribution Hyper- Hyper- Erlang Erlang Distributions Distributions π1 Q1 0 0  1   λ1 λ1 λ1 λ1 0 Q1 0   1 0  π1 ... Q ; ce 0 ... 0 ...    0 λ1  2 ... 0 0    Qm  λ2 0 0 2 0 T   ... ... ... ...  ;  i i 0  π2   0  0 0 ... n 1 0   0 i 0 0   0  λn 0 ... 0 n  λm λm λm   πm ... λm Qi   ...    Ac  0 0 i i  πn n    0 0 0 i  Figure 1. Examples of CPH distributions 2.2 Fitting stochastic times with CPH, Normal and Lognormal distributions In practice, surgery or PACU time distributions are often obtained from empirical data with no explicit probability density function. In fact, closed formulae for the probability density function may not exist at all. Moreover, surgery or PACU times in reality have different behaviors such as high 8 variability, long and sharp tails (Baesler et al. 2015; da Silva Godinho, 2014; Stepaniak et al. 2009). Thus, it is difficult to work with distributions that do not support such behaviors, like the normal distribution. In contrast, CPH distributions provide a good alternative under these conditions (Perez and Riano, 2007). Different algorithms are developed to fit the most appropriate CPH distribution to empirical data. A wide range of CPH fitting algorithms is studied by Perez and Riano (2007). They show that the Expectation-Maximization for Hyper-Erlang (EMHE) algorithm (Thümmler et al. 2006) is the best performer in terms of goodness-of-fit for maximum likelihood estimation. We apply the EMHE algorithm to fit three real datasets of Start Anesthesia to Start Operation t ip (SASO), Surgery, and PACU times. SASO time is defined as the duration of time from starting cr anesthesia procedure until starting the surgery. Surgery time is the time of surgery and the PACU time is the time of patients recovering from anesthesia after the surgery. The three datasets are collected us from orthopedic operating theater in the Scottish NHS hospital from 1998 to 1999, and are used in an other studies (Bowers and Mould, 2004). We assume that Surgery time is independent of SASO and PACU times for two reasons: (1) the correlation between Surgery time with SASO and PACU times M in our real datasets are at most 0.09 and -0.046, respectively; (2) most researchers assume that Surgery time is independent of SASO and PACU times (Baesler et al. 2015; Cardoen et al. 2010b; Lee and ed Yih, 2014; Saadouli et al. 2015). pt We also fit the same set of data by normal and lognormal distributions, as the two most commonly used distributions in the OTR literature. In all three cases, i.e., Hyper-Erlang, normal and lognormal ce we use 80% of the data for training and the remaining 20% for testing out-of-sample error (Cruz et al., Ac 2011; Epprecht et al., 2017; Saeed et al., 2017). The estimated parameters of Hyper-Erlang, normal, and lognormal distributions are shown in Table 1. In this table, the mean and the variance of normal and lognormal distributions are presented by  and  2 and shown by N (  , ) and ln N ( , ) , 2 2 respectively. Hyper-Erlang distribution (Figure 1) is a mixture of m mutually independent Erlang distributions (with parameters of the number of phases and rate parameter ( i )) weighted with the initial probabilities 1 , . . .,  m , where  i  0 , and the vector  is stochastic, i.e.,   i  1 . We m i 1 9 consider the Hyper-Erlang distribution with 10 states to fit the set of data, due to its well-balanced fitting precision and manageable computational time. Parameters and representations of Hyper-Erlang distributions are given in Table 1 and Figure 2. We compare the results of fitting the same set of data by normal and lognormal distributions with the EMHE algorithm. In Figure 3 to Figure 5, the empirical probability functions for the traces of SASO, Surgery, and PACU as well as the fitted density functions of Hyper-Erlang distribution with 10 states, normal and lognormal distribution are shown. To evaluate the goodness-of-fit, we apply the chi- square test to statistically compare the result of each fitting distribution with the frequency of data at a significance level of 5%. Table 2 to Table 4 show the results of chi-square test for each fitting t ip distribution and test data (the 20% out-of-sample data). Based on the results of the chi-square test, cr there is a statistically significant difference among the fitted normal distribution and the data frequency (P-value= 0.00) while there is no evidence of statistical difference between the fitting of us Hyper-Erlang and lognormal distributions with the data frequency (P-value= 1.00). an To evaluate the dispersion of fitting error, we bootstrap from the data (Chernick, 2008) and report the mean squared error of fitting for SASO, Surgery, and PACU on 1000 equally sized datasets. The M results along with the 95% confidence intervals around the mean squared error for each real dataset and the corresponding fitted distributions are given in Table 5. In almost all cases, the Hyper-Erlang ed distribution shows a statistically smaller mean squared error compared to that of the normal and pt lognormal distributions. The 95% confidence intervals around the mean squared error for Hyper- Erlang do not overlap with those of the normal and lognormal distributions, indicating that the error of ce Hyper-Erlang distribution fitting is significantly smaller than the other two distributions. Ac Table 1. Parameters of Hyper-Erlang, Normal, and Lognormal distribution for trace data of SASO, Surgery and PACU times Time SASO Surgery PACU Distribution Number of phases =[2,4,4] Number of phases =[1,3,6] Number of phases =[1,2,3,4] Hyper-Erlang   [0.01; 0.02; 0.87]   [0.08; 0.58; 0.34]   [0.03; 0.13; 0.70; 0.14]   [0.56;1.78; 0.21]   [0.02; 0.04; 0.26]   [0.00; 0.15; 0.57; 2.76] Normal N(20.59,13.48) N(55.50,43.08) N(10.68,42.81) Lognormal LnN(2.82,0.68) LnN(3.71,0.84) LnN(1.51,0.98) 10 Hyper- Hyper- Erlang Erlang Distribution Distribution with with 10 10 States States Hyper- Hyper- Erlang Erlang Distribution Distribution with with 10 10 States States Hyper- Hyper- Erlang Erlang Distribution Distribution with with 10 10 States States Trace Trace Data Data of of SASO SASO Trace Trace Data Data of of Surgery Surgery Trace Trace Data Data of of PACU PACU λλ11== 0.21 λλ11== 0.21 λλ11== 0.21 0.21 λλ11== 0.26 λλ11== 0.26 λλ11== 0.26 0.26 λλ11== 0.26 λλ11== 0.26 λλ11== 2.76 ππ11 =0.87 0.21 0.21 ππ11 =0.34 =0.34 0.26 0.26 0.26 0.26 λλ11=2.76 =2.76 λλ11== 2.76 2.76 2.76 =0.87 ππ11 =0.14 =0.14 λλ11= λλ11= = = λλ11== 00.2.2 22.7.7 11 66 λλ22== 0.57 λλ22== 0.57 00.2.2 λλ22== 1.78 λλ22== 1.78 ππ22 =0.70 =0.70 0.57 0.57 λλ22== 0 66 λλ22== 1.78 1.78 1.78 1.78 0.5 .577 ππ22 =0.02 =0.02 λλ22== 1.78 1.78 ππ22 =0.58 =0.58 λλ22== 0.04 0.04 λλ22== 0.04 0.04 λλ22== 0.04 0.04 .155 ππ33 =0.13 λλ33== 0.15 λλ33== 00.1 =0.13 0.15 ..5566 ==00 λλ33 ππ33 =0.1 =0.1 λλ33== 0.56 0.56 ππ33 =0.08 =0.08 λλ33== 0.02 ππ44 =0.03 λλ44== 0.005 0.005 0.02 =0.03 Figure 2. Hyper-Erlang representation for trace data of SASO, Surgery and PACU times 0.14 SASO Surgery 0.26 0.24 0.12 Hyper-Erlang with 10 Hyper-Erlang with 10 state 0.22 state 0.2 0.1 Normal t Normal 0.18 ip Lognormal 0.08 0.16 lognormal PDF 0.14 PDF 0.06 0.12 0.1 cr 0.04 0.08 0.06 0.02 0.04 0.02 us 0 0 1 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 128 136 144 152 160 168 176 184 192 200 208 216 1 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 SASO time (min) Surgery time (min) Figure 3. Trace data of SASO and the approximating Figure 4. Trace data of Surgery and the approximating an distributions distributions 0.32 0.3 0.28 PACU M 0.26 0.24 Hyper-Erlang with 10 state Table 2- The results of chi-square test for SASO data 0.22 0.2 Normal Distribution Hyper-Erlang Normal Lognormal 0.18 χ2 5.24 549.19 7.44 PDF 0.16 2 lognormal χ -crit 31.41 41.34 41.34 ed 0.14 0.12 0.1 p-value 1.00 0.00 1.00 0.08 0.06 significant No Yes No 0.04 0.02 pt 0 1 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 PACU time (min) ce Figure 5. Trace data of PACU and the approximating distributions Ac Table 3- The results of chi-square test for Surgery Table 4- The results of chi-square test for PACU data data Hyper- Hyper- Distribution Normal Lognormal Distribution Normal Lognormal Erlang Erlang χ2 20.57 85.75 20.37 χ2 8.08 140.86 14.69 2 2 χ -crit 61.66 70.99 70.99 χ -crit 31.41 41.34 41.34 p-value 1.00 0.00 1.00 p-value 0.99 0.00 0.98 significant No Yes No significant No Yes No 11 Table 5- The results of Mean Squared Error, Standard Error and Confidence Interval of SASO, Surgery and PACU data Distribution Hyper-Erlang Normal Lognormal Mean Squared Error 24.79 63.62 49.81 SASO Standard Error 0.19 0.40 0.33 Confidence Interval [24.40, 25.16] [62.82, 64.41] [49.15, 50.46] Mean Squared Error 139.19 309.04 156.84 Surgery Standard Error 0.57 0.83 0.59 Confidence Interval [138.04, 140.33] [307.37, 310.70] [155.66, 158.01] Mean Squared Error 905.02 17881.44 1707.81 PACU Standard Error 11.66 134.42 19.86 Confidence Interval [881.69, 928.35] [17612.60, 18150.27] [1668.09, 1747.52] t ip Through these examples we can observe that the CPH (Hyper-Erlang) distribution outperforms the other two distributions in terms of fitting quality for each dataset. This higher performance is due to cr the special properties of the CPH distributions as discussed earlier in this section, which renders them an attractive choice for fitting surgery procedure times. us In the next subsection, we discuss the steps to calculate the convolution of three independent an stochastic times (i.e., SASO + surgery + PACU = OTR time) based on CPH, normal and lognormal M distributions and compare the outcome. 2.3 Calculating OTR time by CPH, Normal and Lognormal distributions ed The process of surgery in the OTR problem consists of a sequence of several stages, including SASO, pt surgery and PACU. We assume that one patient from the set surgery and goes through the successive stages as illustrated in Figure 6. The total duration of SASO, ce surgery and PACU for this patient is denoted by TS (T S  t SASO  t Surgery  t PACU ) and referred to as the Ac OTR time. The hours that ORs are open and operational in each day is a strategic decision and hospital managers consider definite thresholds that we denote by T, in each day for the surgery of elective patients. In order to design an efficient schedule to ensure the timely completion of surgery for the last patient, we assume that every patient p  finishes his/her operation before a threshold T with a probability greater than  : 12 P (T S  T )   p  . (2) In all calculations, we set the β = 0.85 as an arbitrary probability value and assume that the threshold T is equal to the available OTR service time in the hospital reported in the test problem. Our aim in this subsection is to compare the performance of CPH distribution with normal and lognormal distributions in modeling the OTR times, i.e., TS (T S  t SASO  t Surgery  t PACU ). We use the same examples presented in Section ‎2.2 for each stage in Figure 6. Figure 6. Patient flow for surgery t ip Depending on the probability distributions assumed for the duration of each stage in Figure 6, cr computing the OTR time can be a straightforward procedure (e.g., CPH distributions and normal distributions due to their additive nature), or a complex operation for those distributions that lack this us property (e.g., lognormal distributions). Based on Moschopoulos (1985) and Shylo et al. (2012), the an analytical expression for the convolution of lognormal distributions is complicated and there is no explicit formula for it. To overcome such a situation for the lognormal distribution, we use Monte M Carlo (MC) simulation algorithm, as one of the most popular approaches in the OTR problem (Shylo et al. 2012). ed In order to compare the CPH distribution with normal and lognormal distributions in modeling the OTR time, we approximate SASO, Surgery, and PACU times with these distributions and calculate pt their sum in the cases of CPH and normal distributions and apply MC simulation for the lognormal ce distribution. In what follows, we demonstrate how to model OTR times using a CPH, normal and lognormal approximations. Ac 2.3.1 CPH approximation We use the Hyper- Erlang distribution with 10 states to model the distribution of each stage (SASO, Surgery, and PACU time) as shown in Table 1. Based on the closure properties of CPH distribution in Appendix A, the convolution of k CPH random variables results in a CPH distribution with distinct parameters. Therefore, the OTR time (i.e., T S  t SASO  t Surgery  t PACU ) also follows the CPH 13 Table 6. Initial probability of OTR time distribution Element of Element of Element of Value Value Value Matrix (i, j) Matrix (i, j) Matrix (i, j) (1,1) 0.09 (1,3) 0.87 (1,7) 0.02 Table 7. Square matrix of OTR time distribution Element Element Element Element Element Element of Matrix Value of Matrix Value of Matrix Value of Matrix Value of Matrix Value of Matrix Value (i, j) (i, j) (i, j) (i, j) (i, j) (i, j) (1,1) -0.06 (9,9) -1.77 (14,14) -0.04 (20,20) -0.26 (22,22) -0.15 (2,27) 0.00 (1,2) 0.06 (9,10) 1.77 (2,15) 0.02 (2,21) 0.00 (22,23) 0.15 (6,27) 0.00 (2,2) -0.06 (10,10) -1.77 (6,15) 0.07 (6,21) 0.00 (23,23) -0.15 (10,27) 0.00 (3,3) -0.21 (2,11) 0.00 (10,15) 0.60 (10,21) 0.00 (2,24) 0.00 (11,27) 0.00 (3,4) 0.21 (6,11) 0.02 (15,15) -0.26 (11,21) 0.00 (6,24) 0.00 (14,27) 0.01 (4,4) -0.21 (10,11) 0.15 (15,16) 0.26 (14,21) 0.00 (10,24) 0.00 (20,27) 0.04 (4,5) 0.21 (11,11) -0.02 (16,16) -0.26 (20,21) 0.01 (11,24) 0.02 (27,27) -2.76 (5,5) -0.21 (2,12) 0.03 (16,17) 0.26 (21,21) -0.01 (14,24) 0.03 (27,28) 2.76 t (5,6) 0.21 (6,12) 0.12 (17,17) -0.26 (2,22) 0.00 (20,24) 0.18 (28,28) -2.76 ip (6,6) -0.21 (10,12) 1.02 (17,18) 0.26 (6,22) 0.00 (24,24) -0.57 (28,29) 2.76 (7,7) -1.77 (12,12) -0.04 (18,18) -0.26 (10,22) 0.00 (24,25) 0.57 (29,29) -2.76 cr (7,8) 1.77 (12,13) 0.04 (18,19) 0.26 (11,22) 0.00 (25,25) -0.57 (29,30) 2.76 (8,8) -1.77 (13,13) -0.04 (19,19) -0.26 (14,22) 0.01 (25,26) 0.57 (30,30) -2.76 (8,9) 1.77 (13,14) 0.04 (19,20) 0.26 (20,22) 0.03 (26,26) -0.57 us distribution. We calculate the distribution of OTR time by the convolution of the three fitted Hyper- an Erlang distributions. The result of OTR time distribution based on its initial probability ( T ) and S M square matrix (TT ) are presented in form of a sparse matrix in Table 6 and Table 7. Given the S parameters of the resulting CPH distribution and using its cumulative distribution function presented ed in Equation A.1, we compute P (T S  T ) . pt 2.3.2 Normal approximation ce The convolution of k normal random variables is normally distributed as N (i 1 i , i 1 i2 ) . The k k estimated parameters of the normal distribution for each stage are shown in Table 1. Therefore, we Ac calculate the OTR time distribution as the convolution of three fitted normal distributions resulted as N(86.7, 3870.3). Given these parameters we can now use the normal distribution function to compute P (T S  T ) . 14 2.3.3 Lognormal approximation Since the lognormal distribution lacks the additive property, we resort to MC simulation to estimate P (T S  T ) . The algorithm of MC simulation to estimate the P (T S  T ) is shown in Figure 7. In this algorithm, N 0 independent random variates are generated for each random variable (i.e., SASO, Surgery and PACU time) from their corresponding fitted lognormal distributions with parameters shown in Table 1. More specifically, for each stage depicted in Figure 6, random variable t s , describing time in each stage s  S and S  {SASO,OR,PACU} , is sampled N 0 times, and then the OTR time ( TS ) is calculated as the sum of these three random variates for N 0 samples. The t ip P (T S  T ) is then estimated based on the number of OTR times that are less than T out of the N 0 cr samples. us To correctly interpret the inherent variability of the simulated samples we define a confidence interval (CI) as shown in Inequality (3) based on Oberle (2015). an pˆ (1  pˆ ) pˆ (1  pˆ ) pˆ  z   p  pˆ  z  (3) 2 N0 2 N0 M In Inequality (3), p̂ is the natural estimator of P (T S  T ) , N 0 is the number of generated random  variates,  is the level of confidence, and z  is the z-score associated with ed . 2 2 Input: S, sequence of stages in surgery process; {𝑡𝑠̃ }𝑠∈𝑆 random variable describing time in each stage 𝑠 ∈ 𝑆; T, pt Threshold; 𝛽, required probability; 𝑁0 , number of random variates to generate. Output: true if 𝑃(𝑇̃𝑆 ≤ 𝑇) ≥ 𝛽, false otherwise. 1: 𝑗 ← 0 % j represents number of replications in which 𝑇̃𝑆 ≤ 𝑇 ce 2: for i=1 to 𝑁0 do 3: 𝑇̃𝑆𝑖 ← 0 % 𝑇̃𝑆𝑖 represents the total duration in one replication 4: for all 𝑠 ∈ 𝑆 do 𝑇̃𝑆𝑖 ← 𝑇̃𝑆𝑖 + generate(𝑡̃𝑠 ) % generates one sample from distribution 𝑡̃𝑠 Ac 5: 6: end for 7: if 𝑇̃𝑆𝑖 ≤ 𝑇 then 8: 𝑗 ←𝑗+1 9: end if 10: end for 11: 𝑗 12: if ≥ 𝛽 then 𝑁0 13: return true 14: else 15: return false 16: end Figure 7. Monte Carlo (MC) simulation algorithm 15 In all cases, we use N 0  10000 as suggested by Shylo et al. (2012). Three levels of confidence, i.e., 95%, 97.5%, and 99%, are considered to compare the results. We compare this CI with the probabilities resulted from the two other approximations, i.e., CPH and normal. 2.3.4 Results of comparing the three approximations In order to compare the CPH approximation with real data, normal and lognormal approximations, we present their results by illustrating the distribution of OTR time ( TS ) for each case as shown in Figure 8. Based on this figure, the lognormal and CPH approximations perform similarly in mimicking the real data, while the normal approximation performs significantly worse in comparison with two other t ip approximations. cr 0.25 Real Dara 0.2 CPH us Normal Probability 0.15 Lognormal 0.1 an 0.05 M 0 1 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 100 Time (min) Figure 8. Probability density function of TS for CPH, lognormal, and Normal approximations ed To compare the result of P (T S  T ) for CPH and normal approximations with the CI of P (T S  T ) pt for the lognormal approximation, we show the numerical example for P (T S  T ) with different ce possible threshold times ( T ) for three defined levels of confidence, i.e., 95%, 97.5% and 99% (Table 8). The results show that the P (T S  T ) for CPH approximation falls within the CI of lognormal Ac approximation while this probability for the normal approximation is out of CI. For instance, the P (T S  T ) for CPH approximation in T=150 is 0.8943 and this probability is in the CI of lognormal approximation ([0.8678, 0.9121]) while this probability for normal approximation (0.8453) is out of the CI. Therefore, by comparing the three approximations presented above, we conclude that there is 16 no statistical difference between the CPH and lognormal approximations, and they both perform better than the normal distribution approximation. Although MC simulation is easy to implement, it needs a large number of samples to get an acceptable precision in the results. Also, the unavailability of an explicit probability distribution function after simulation is another disadvantage of the MC simulation algorithm (Bandyopadhyay and Bhattacharya, 2014). Therefore, we can conclude that the CPH distribution is the most attractive choice in terms of precision and ease of implementation in the OTR scheduling problems. ̃ 𝑺 ≤ 𝑻) of CPH, lognormal, and normal approximations and CI of lognormal approximation Table 8. 𝑷(𝑻 Approximation CPH Normal Lognormal CI of Lognormal Approximation t Probability Approximation Approximation Approximation   0.975   0.95   0.9 ip P (T S  150) 0.8943 0.8453 0.8900 [0.8678,0.9121] [0.8706,0.9093] [0.8737,0.9062] P (T S  140) 0.8687 0.8039 0.8630 [0.8386,0.8873] [0.8416,0.8843] [0.8451,0.8808] cr P (T S  130) 0.8375 0.7565 0.8440 [0.8182,0.8697] [0.8215,0.8664] [0.8251,0.8628] us P (T S  120) 0.7998 0.7034 0.8130 [0.7853,0.8406] [0.7888,0.8371] [0.7927,0.8332] P (T S  110) 0.7552 0.6456 0.7750 [0.7454,0.8045] [0.7491,0.8008] [0.7532,0.7967] P (T S  100) an 0.7032 0.5842 0.7280 [0.6964,0.7595] [0.7004,0.7555] [0.7028,0.7511] P (T S  90) 0.6435 0.5207 0.6630 [0.6294,0.6965] [0.6337,0.6922] [0.6384,0.6875] P (T S  80) 0.5750 0.4567 0.5940 [0.5591,0.6288] [0.5635,0.6244] [0.5684,0.6195] M P (T S  70) 0.4946 0.3938 0.4880 [0.4525,0.5234] [0.4570,0.5189] [0.4620,0.5139] P (T S  60) 0.3962 0.3335 0.3890 [0.3544,0.4235] [0.3587,0.4192] [0.3636,0.4143] ed 3 Stochastic Model for the OTR problem pt Since the OTR scheduling problem, defined in Section 1 as ce el pat; time; room;int( PACU ) stoch single ; Cmax , is NP-Hard, (Arnaout and Kulbashian, 2008; Ac Baesler et al. 2015; Latorre-Núñez et al. 2016), we use heuristic algorithms to find near optimal solutions. Our strategy is to connect the CPH distribution with constructive heuristic and meta- heuristic algorithms, for searching better solutions. Our proposed search algorithm for OTR problem, inspired by Baesler et al. (2015) and Gómez et al. (2015), has two components. The first component is the sequence generator engine which is responsible for exploring the solution space, generating a new sequence for patients to build a 17 solution. The OTR literature includes a number of efficient algorithms that can play the role of the sequence generator engine. The most common algorithms are simulated annealing, tabu search, genetic algorithms, longest expected processing time (LEPT), among others (Baesler et al. 2015; Cardoen et al. 2010b). The second component is the schedule evaluator that computes the performance metrics for a sequence of patients such as the probability of completing the surgery before the set threshold T, P (T S  T ) . The interaction between these two components is depicted in Figure 9. To better illustrate our proposed algorithm, we first apply it on a small-size problem in this section. t Then, in Section 4, we present several heuristic algorithms that can be paired with our proposed ip schedule evaluator. cr New New Sequence Sequence  Initial Initial condition  Set Set of condition of patients patients Sequence Generator Heuristic Heuristic Algorithm Algorithm us Schedule Evaluator Performance Measure Final Final schedule schedule an Meta-heuristic Meta-heuristic Algorithm Algorithm P (TS  T ) OTR OTR Schedule Schedule M Figure 9. Search algorithm for OTR problem Consider two operating rooms and one recovery bed as depicted in Figure 10 with stochastic surgery ed ( t Surgery ) and PACU ( tPACU ) times. We assume that for t Surgery two distributions scenarios can happen: pt the first scenario happens with probability p1 and has low variance and the second scenario has long ce surgery durations and relatively higher variance with probability p2 . Let us use an Erlang distribution to model the low-variance scenario and the exponential distribution to model the high-variance Ac surgery time. We also assume that the PACU time is modeled with Erlang distributions. 18 OTR OTR Environment Environment OR PACU OR1 PACU1 OR2 Figure 10. OTR environment In this example, we first compute the convolution of surgery and PACU times (i.e., the OTR time) for one patient and show the PH representation of the OTR time. Then, we consider four patients, two t patients with type '1' surgery and two patients with type '2' surgery, and determine the schedule of ip each patient by considering P (T S  T ) . cr 3.1 The OTR problem for one patient us To compute the OTR time for one patient, we assume that the surgery time ( t Surgery ) includes two an distributions scenarios: Erlang distribution with probability p1 and an exponential distribution with probability p2 . The Erlang distribution is the sum of two exponential phases with the same intensity M 1 . The exponential distribution is parameterized with 2 . Let p1  0.6 , p2  0.4 , 1  0.04 , ed 2  0.01 . Moreover, assume that the PACU time ( tPACU ) follows an Erlang distribution with two pt phases. The t Surgery and tPACU are PH C (Surgery , TSurgery ) and PH C (PACU , TPACU ) respectively, and ce follow the PH representation with the structure shown in Figure 11, Equations (4) and (5).  0.04 0.04 0  Surgery  [0.6,0,0.4] , TSurgery  0  Ac  0 0.04 (4)  0 0 0.01  0.01 0  (5) PACU  [0.3,0.7] , TPACU   .  0 0.5  19 Surgery Surgery PACU PACU 0.04 0.04 0.6 0.01 0.3 0.4 0.7 0.5 0.01 Figure 11. Structure of CPH distribution for surgery and PACU time The OTR time of a patient (TOTR  t Surgery  t PACU ) is the convolution of CPH distributions for surgery and PACU times as shown in Figure 11, which is distributed PH C (OTR , TOTR ) with: t ip  0.04 0.04 0.00 0.00 0.00   0.00 0.04 0.00 0.01 0.03  cr  OTR  [0.6,0,0.4,0,0] , TOTR   0.00 0.00 0.01 0.00 0.01  . (6)   us  0.00 0.00 0.00 0.01 0.00   0.00 0.00 0.00 0.00 0.50  Using the formula introduced in Appendix A, we compute the required performance measures for an PH C (OTR , TOTR ) . The expected duration of OTR time is E [TOTR ]  101.4 , its variance is M Var [TOTR ]  10370 , the probability that the OTR time is completed before 150 minutes and 270 minutes as arbitrary threshold times are P (TOTR  150)  0.79 and P (TOTR  270)  0.93 , respectively. ed 3.2 The OTR problem with four patients pt To present the schedule of four patients, assume that the surgery time of type '1' is ce PH C (Surgery 1 , TSurgery 1 ) following Equation (4) and the surgery time of type '2' is Ac PH C (Surgery 2 , TSurgery 2 ) with:  0.05 0.05 0.00  Surgery 2  [0.3,0,0.7] , TSurgery 2   0.00 0.05 0.00  .  (7)  0.00 0.00 0.01 The PH representation of PACU time is given in Equation (5). We assume that the sequence of patients is pa2- pa1- pa1- pa2 where pai is a patient with type 'i' surgery, i {1, 2} . We also consider 20 that the ORs can start the first surgery at 8:00 AM and should complete the last one by 12:30 PM. i.e., the ORs are only available for 4:30 hours or 270 minutes. In this instance, the goal is to assign these four patients to two ORs to make sure that all surgeries are performed within the ORs’ available time with a probability ( P (T S  T ) ) greater than or equal to an arbitrary probability value such as 0.85. Figure 12 depicts the schedule of assigning four patients to two ORs. Based on the sequence of pa2- pa1- pa1- pa2, one patient wih type '2' surgery and one wih type '1' surgery start their operation at the beginning of the planning horizon simultaneously. The expected surgery completion times of these two surgeries are 82 and 70 minutes, respectively. The expected surgery completion time of patient t with type '1' surgery is shorter than patient with type '2' surgery, then patient with type '1' surgery is ip assigned to PACU and patient with type '2' surgery has to remain on operating room (OR2) after the cr completion of the surgery until PACU becomes available (i.e., blocking OR2). This patient leaves us OR2 when the recovery of patient with type '1' surgery is completed. The P (T S  T ) for each patient is calculated by cumulative distribution function of the CPH presented in Appendix A, Equation A.1. an For instance, the P (T S  T ) for patient with type '1' surgery in OR1 is calculated by Equation (8). M 0.04 0.04 0     0 0.04 0 270 (8) TSurgery1 270  0 0.01 Pr(TS  270)  1   Surgery1 e 1  1  [0.6,0,0.4] e 1  0.9730 0 ed Since the P (T S  T ) for both patients is greater than 0.85, we can operate the patients in both ORs. pt However, the second patient with type '2' surgery (red box in Figure 12) should be operated in the next day because the P (T S  T ) is less than 0.85. ce 82 82 183.4 183.4 Ac Pr(T≤270)= Pr(T≤270)=0.9530 0.9530 Pr(T≤270)= Pr(T≤270)=0.7961 0.7961 OR2 OR2 2 2 OR1 OR1 1 1 70 70 140 140 Pr(T≤270)= Pr(T≤270)=0.9730 0.9730 Pr(T≤270)= Pr(T≤270)=0.9009 0.9009 PACU PACU 1 2 1 101.4 101.4 132.8 132.8 171.4 171.4 Figure 12. Scheduling of four patients 21 Hospitals usually use the average surgery time (AST) to schedule the surgeries, disregarding the inherent variability of the surgery durations which results in the low utilization ORs (M. J. Campbell and Swinscow, 2011; Master et al. 2016). Based on this example and Figure 12, when we consider the AST criterion instead of P (T S  T ) for scheduling these four patients, the second patient with type '2' surgery (red box in Figure 12) can still be operated on the same day because based on AST the completion time of surgery for the second patient with type '2' surgery (183.4) is still less than the threshold time of 270 minutes. However, considering the variability and skewness inherent in the OTR time, the AST method renders an t infeasible schedule with respect to the probability of completion before 270 minutes. ip 4 Heuristic algorithms for OTR problem cr Heuristic algorithms have been used for the OTR scheduling problem. Based on Cardoen et al. us (2010b), we consider two categories of such algorithms: (1) constructive heuristics and (2) meta- an heuristics. In this section, we first discuss several constructive heuristics and then meta-heuristic algorithms based on simulated annealing (SA) and tabu search (TS). We apply them in the search M algorithm (Figure 9) as a "sequencing generator engine" and then deploy CPH methods to evaluate the generated schedules. Constructive heuristics, SA and TS are well known heuristics in the OTR ed scheduling problem with reported high performance. At the end of this section, we present a theoretical upper bound on the performance of this search algorithm. pt 4.1 Constructive heuristic algorithm for OTR problem ce We choose 12 known constructive heuristic algorithms in flowshop scheduling problem (Baker and Ac Trietsch, 2013; Pinedo, 2012) to determine the sequence of patients. We use the following constructive heuristics: (a) increasing expected mean of surgery time (Shortest Expected Surgery Processing Time, SESPT); (b) increasing expected mean of OTR time (Shortest Expected OTR Processing Time, SEOPT); (c) decreasing expected mean of surgery time (Longest Expected Surgery Processing Time, LESPT); (d) decreasing expected mean of OTR time (Longest Expected OTR Processing Time, LEOPT); (e) half increasing in surgery time and half decreasing in surgery time 22 (SESPT-LESPT); (f) half increasing in OTR time and half decreasing in OTR time (SEOPT-LEOPT); (g) half decreasing in surgery time and half increasing in surgery time (LESPT- SESPT); (h) half decreasing in OTR time and half increasing in OTR time (LEOPT- SEOPT); (i) mixed surgery time (MST) based on work by Marcon and Dexter (2006); (j) mixed OTR time (MOT) based on work by Marcon and Dexter (2006); (k) Campbell, Dudek, and Smith (CDS) algorithm based on work by H. G. Campbell et al. (1970); (l) Nawaz, Enscor, and Ham (NEH) algorithm based on work by Nawaz et al. (1983). 4.2 Tabu Search algorithm for the OTR problem t The TS algorithm, developed by Glover (1989, 1990), uses sequence generator techniques to find a ip near-optimal solution for optimization problems. The algorithm needs an initial solution, which is cr called the initial seed. Then, the neighborhood of the seed is examined to turn the search to the best us neighborhood solution as the new seed. The process of moving from the current seed towards the new selected seed is called a move. In order to avoid cycling, the moves that would take the new seed back an to the recently visited solutions are forbidden or penalized for a certain number of iterations. The features of the forbidden moves are kept in a list, called the tabu list (TL). A forbidden move is M ignored when a generated neighborhood based on a tabu move provides a solution with an objective ed function value better than all observed solutions so far. In this case, the seed with the best objective function value is chosen as the new seed. The best found objective function value of the search so far pt is called the aspiration level. The stopping criteria used in TS algorithm is based on time of executing ce the algorithm. In our research problem, a solution represents a particular patient sequence that can be represented as Ac a list of patients who must be served in order of the list. We define two kinds of neighborhoods: (1) pairwise exchange (PE), which means that two patients of the same sequence are selected randomly to exchange places, which generates a new neighbor solution. (2) Best constructive heuristic (BCH), which mean that two slots of patients are selected randomly and the best constructive heuristic is applied for the patients sequence between two selected slots. 23 In this work, we propose two different meta-heuristic algorithms based on TS to solve our research problem: (1) TS with PE (TS-PE); (2) TS with BCH (TS- BCH). In TS-PE algorithm, the TL keeps the type of two patients exchanged positions. For instance, consider four-patient and assume that the current seed, i.e., pa1 – pa2 – pa3 – pa4, is switched to pa1 – pa4 – pa3 – pa2. In this case, since the position of patients in the second and fourth slots are changed, for a finite number of moves (the number of these moves is determined by the TL size), the patient with surgery of type two and four cannot be changed their position. In TS- BCH algorithm, the TL keeps track of the slots to which the patients are assigned, i.e., the slots of patients are memorized in the TL. In mentioned instance, the position of patients in the second and the fourth slots cannot be changed. t ip 4.3 Simulated Annealing algorithm for the OTR problem cr SA, introduced by Kirkpatrick et al. (1983), is an iterative random search technique that has been us widely applied to solve various combinatorial optimization problems including the OTR problem ' (Cardoen et al. 2010a; Sier et al. 1997). At each iteration, a new solution ( S ) is obtained by the an movement into the neighborhood of the current solution ( S ) . If the new solution is better, i.e., M Cmax (S ' )  Cmax (S ) , then the move is admitted and the algorithm continues. Otherwise, the Cmax ( S ) Cmax ( S ' ) ed degrading solution is admitted as the current solution with the probability exp T , where T is the temperature. At first, the initial temperature (T0 ) is set high and gradually reduced through pt multiplication by a cooling factor  . The algorithm carries on until convergence by appropriately ce reducing the temperature or when the time of executing the algorithm is finished. Ac In this work, we propose two different meta-heuristic algorithms based on SA to solve our research problem based on two kinds of neighborhoods defined in Section ‎4.2: (1) SA with PE (SA-PE); (2) SA with BCH (SA- BCH). 4.4 Upper bound on the search algorithm In this section, we present the worst-case ratio bound (ρ) for every heuristic H we have proposed in our search algorithm. Let SH be a schedule generated by a heuristic H and let S* be the optimal 24 schedule. Heuristic H is defined as the worst-case ratio bound ρ if for any problem instance the following equation holds (Koulamas and Kyparisis, 2000): C max (S H )  (9) C max (S * ) Our proposed research problem can be considered as a flexible flowshop scheduling problem with two stages. Koulamas and Kyparisis (2000) presented the worst-case ratio bound for flexible flowshop scheduling problem with m stages as following when the processing times are stochastic. C max (S H ) Pmax  1  mf ( m ) l max (10) C max (S * ) Psum t ip where cr - Pmax  max{ pi , j } and Psum   pi , j , pi,j is the process time of job j in stage i, i,j i j - us l max  max{l1 ,..., l m } , li is the number of machine in stage i=1,..,m, an - f (m )  m 2 1. We assume that the surgery and PACU times are to be independent, identically distributed CPH M random variables drawn from PH C (Surgery , TSurgery ) and PH C (PACU , TPACU ) . Since our times are ed CPH distributed, we use the properties of “Convolution of PHC” and “Maximum of PHC”, shown in Appendix A, to calculate the distribution of Pmax and Psum. PH C ( Pmax , TPmax ) and PH C ( Psum , TPsum ) pt are distributions of Pmax and Psum. The worst-case ratio bound for our research problem by using ce Pmax E [Pmax ] 1 E[ ] and E [X ]  (- T) 1 is as follows: Psum E [Psum ] Ac C max (S H ) Pmax  Pmax (TPmax )11 E[ ]  1  6l max E [ ]  1  6l max . (11) C max (S * ) Psum  Psum (TPsum )11 5 Computational results In this section, we first compare the performance of constructive heuristic and meta-heuristic algorithms in terms of efficiency in subsection 5.1. Then, we compare the best meta-heuristic 25 algorithm found in subsection 5.1 with the only existing algorithm similar to our research problem in literature, to the best of our knowledge. At the end of this section, we demonstrate the performance of our proposed search algorithm (Figure 9) by comparing the CPH and AST methods in the schedule evaluator component and applying the best meta-heuristic algorithm as a sequence evaluator found in subsection 5.1. 5.1 Results of heuristic algorithms In this research, 12 constructive heuristic and four meta-heuristic algorithms are proposed for the el pat; time; room;int( PACU ) stoch single ; Cmax problem. In order to evaluate the performance of t ip these algorithms, we generate the test problem by assuming the number of patients between 5 to 500 cr in three different categories: small (problems with 5–100 patients for scheduling of one week), medium (problems with 100 to 200 patients for scheduling of two weeks), and large (problems with us 200 to 500 patients for scheduling of more than two weeks to one month). Each category is divided an into subcategories as intervals with five units. For example, in the first category, we consider [1- 5], [6- 10], …, [96-100] as subcategories. Then we generate 20 cases randomly from a uniform discrete M distribution for each subcategory to assign surgeries to patients. Generally, we generate 2000 different test problems, 400 for small, 400 for medium, and 1200 different test problems for large categories. ed Finally, to test the performance on real data we bootstrap from our datasets presented in Subsection 2.2 and generate 23 new instances of surgery types, motivated by Baesler et al. (2015). In each 2000 pt different test problem, the surgery type of each patient is identified by the number between 1 to 23. ce We assume that there are 3 operating rooms and 6 recovery beds and the ORs are available only for 8 hours, from 8:00 AM to 16:00 PM on each day. Ac All proposed algorithms are coded in MATLAB 2014a (MATLAB Release 2014a, 2015) and executed on a personal computer with Intel(R) Core(TM) i3-2120, running at 3.30 GHz with 8 GB RAM. To compare the performance of the proposed constructive heuristics, we consider 2000 generated test problems for all categories and solve each one by each proposed constructive heuristic algorithm. We calculate the objective function of each test problem (minimization of makespan) and then compute 26 the preference percentage of each constructive heuristic. The results of preference percentage for each constructive heuristic are presented Figure 13 and the average, standard error (SE), and range of preference percentage for each heuristic in different categories of test problems are shown in Table 9. For example, the average of preference percentage for CDS algorithm and category of [1-100] is 33.81%. It means that in 33.81% of the test problems in the [1-100] category, the objective function of CDS algorithm is better than other constructive heuristic algorithms. Based on the results, the average of preference percentage for both CDS and LEOPT algorithms are the same in 1 to 200 patients and both of them perform better than other heuristics for the number of mentioned patients. For 200 to 500 patients, both LESPT and LEOPT algorithms outperform other heuristics. The running time and the t ip complexities of the proposed constructive heuristics are presented in Table 10 where n is the number cr of patients and m is the number of stages (m=2 in our research problem). Based on Table 10, the running time of CDS algorithm is smaller than LEOPT algorithm in 1 to 200 patients. For 200 to 500 us patients, the running time of LESPT algorithm is better than LEOPT algorithm. We apply the results an of these experiments as a neighborhood structure (BCH) of SA- BCH and TS- BCH algorithms. Therefore, we consider the BCH by selecting two slots of patients randomly and if the number of M patients between two slots are less than 200 patients, then we apply the CDS algorithm to order these special patients. Otherwise, we use the LESPT rule to order these special patients. ed SESPT LESPT SESPT-LESPT LESPT-SESPT MST CDS pt NEH SEOPT LEOPT SEOPT-LEOPT LEOPT-SEOPT MOT 0.6 ce 0.5 0.4 Ac Percentage 0.3 0.2 0.1 0 5 35 65 95 125 155 185 215 245 275 305 335 365 395 425 455 485 Number of Patients Figure 13. The preference percentage of each heuristic 27 Table 9. The average, standard error (SE), and range of preferencem percentage for each heuristic in different test problem categories [1-100] [100-200] [200-300] [300-400] [400-500] Algorit Avera Ran Avera Ran Avera Ran Avera Ran Avera Ran hm SE SE SE SE SE ge ge ge ge ge ge ge ge ge ge 0.00 0.0 0.0 0.0 0.0 SESPT 0.64% 0.05 2.64% 0.19 1.94% 0.10 3.41% 0.20 2.66% 0.10 4 1 1 1 1 13.75 0.02 16.73 0.0 29.77 0.0 27.76 0.0 28.52 0.0 LESPT 0.30 0.30 0.31 0.45 0.25 % 3 % 2 % 2 % 2 % 2 SESPT- 0.00 0.0 0.0 0.0 0.0 1.51% 0.10 1.45% 0.10 4.33% 0.15 6.83% 0.15 5.68% 0.15 LESPT 6 1 1 1 1 LESPT- 0.00 0.0 0.0 0.0 0.0 0.91% 0.08 5.63% 0.14 0.99% 0.10 2.68% 0.15 4.10% 0.15 SESPT 5 1 1 1 1 0.00 0.0 0.0 0.0 0.0 MST 4.15% 0.15 0.50% 0.15 0.54% 0.06 0.00% 0.00 0.00% 0.00 9 1 0 0 0 33.81 0.01 23.70 0.0 0.0 0.0 0.0 CDS 0.28 0.36 8.71% 0.20 6.54% 0.15 6.44% 0.25 % 9 % 2 1 1 1 0.00 0.0 0.0 0.0 0.0 NEH 0.64% 0.05 5.12% 0.15 5.43% 0.15 6.87% 0.20 2.99% 0.15 4 1 1 1 1 0.00 0.0 0.0 0.0 0.0 t SEOPT 0.25% 0.05 0.74% 0.05 0.75% 0.05 1.95% 0.10 2.25% 0.10 3 0 0 1 1 ip 32.48 0.01 23.66 0.0 28.04 0.0 26.71 0.0 28.08 0.0 LEOPT 0.26 0.21 0.35 0.47 0.35 % 9 % 1 % 2 % 3 % 2 SEOPT 10.10 0.01 15.43 0.0 15.83 0.0 13.84 0.0 15.05 0.0 cr - 0.30 0.30 0.30 0.24 0.25 % 7 % 2 % 2 % 2 % 1 LEOPT LEOPT 0.00 0.0 0.0 0.0 0.0 - 0.76% 0.05 2.41% 0.10 1.95% 0.10 2.19% 0.10 3.96% 0.10 us 4 1 1 1 1 SEOPT 0.00 0.0 0.0 0.0 0.0 MOT 1.01% 0.05 1.98% 0.10 1.70% 0.10 1.23% 0.10 0.25% 0.05 5 1 1 1 0 Table 10. Running Time and Complexity of proposed constructive heuristics an Running Time Complexity of Algorithm Constructive [1-200] [200-500] Heuristics M SESPT 0.64 1.41 LESPT 0.56 1.26 SESPT-LESPT 0.56 1.29 LESPT-SESPT 0.59 1.36 ed MST 0.54 1.10 O(n*log(n)) SEOPT 0.55 1.22 LEOPT 0.66 3.13 pt SEOPT-LEOPT 0.66 1.44 LEOPT-SEOPT 0.58 1.25 MOT 0.55 1.29 ce CDS 0.57 1.31 O(n2) NEH 0.55 1.15 O((m-1)*n2) Ac To present the appropriate meta-heuristic algorithms, we use 2000 test problems for three different categories of patients (small, medium, and large). The solution time for the test problems is considered based on the level of the test problem, and each problem is solved by 120s, 300s, and 600s for small, medium, and large problem, respectively. In order to tune the algorithms to solve different sized problems, extensive experiments were performed to find suitable values for the parameters. The empirical formulae for the value of the SA and TS parameters are presented in Appendix B. 28 To statistically compare the performance of the four mentioned meta-heuristic algorithms, a single- factor (algorithm effect) experiment with four levels (number of algorithms) is performed. The hypothesis to be tested for this purpose is: H0: SA-PE = SA- BCH = TS-PE= TS- BCH H1: Otherwise (12) The single factor experiment is conducted as follows: Y ab    methoda  blockb  ab (13) where Yab: The response variable; μ: The overall mean; methoda: The meta-heuristic algorithms (effect factor), a = 1, 2, 3, 4. t blockb: The test problem’s number factor (random factor), b = 1,…, 2000. ip  ab : The error term cr This is a mixed model, since it includes fixed factors (algorithms) as well as random factor (problem us instances). The normal probability plots of the residuals and equality of variance are demonstrated in Appendix C, Figure C1 and Tables C1 to C3. The residual plots are close to a straight line and there is an no statistically significant difference among the variances of algorithms. The results show that M ANOVA can be used to analyze the results. In our experimental design, the four meta-heuristic algorithms are considered as the effect factor and the test problem’s number factor are considered as ed blocks. The experimental design is coded in the Statistical Analysis System (SAS) software package, release 9.1 (SAS Release, 2003) with a significance level of 5%. The ANOVA tables are provided in pt Appendix C. Based on the results of the ANOVA, there is a statistically significant difference among ce the algorithms (the p-values are less than 0.05). In order to find the meta-heuristic algorithm with the best performance, we also perform the Tukey Ac test (as detailed in Appendix C). The results show that there are significant differences in their performance of all proposed meta-heuristic algorithms. Since TS- BCH algorithm provided the minimum average (Table C6 in Appendix C), it can be considered as the best algorithm among the proposed ones for our research problem. 29 5.2 Comparing the TS-BCH and Baesler’s algorithms In this section, we compare the performance of TS-BCH algorithm presented in the previous subsection with the performance of the only other applicable algorithm to our problem presented by Baesler et al. (2015). We use 2000 test problems for three different categories of patients (small, medium, and large). We also consider 23 surgery types based on Baesler et al. (2015) and assume that there are 3 operating rooms and 6 recovery beds. The distributions for surgery and PACU times of the 23 surgery types are presented in Table 11. We employ the EMHE algorithm (Thümmler et al. 2006) to fit the presented distributions by Hyper-Erlang distribution when we apply the TS-BCH algorithm. Since the assumption of approximate normality and equality of variance are satisfied, we apply paired t ip t-test experiment to compare the performances of the two algorithms (TS-BCH and Baesler’s cr algorithm) for 2000 test problems. The result of normal probability plots of the residuals for the two algorithms, equality of variances and the paired t-test are shown in Appendix D. Since the p-value for us the experiment is almost equal to zero, we can conclude that there is a statistically significant an difference between the results of the two algorithms. Since the average of the objective function of the TS-BCH algorithm is lower than the Baesler’s algorithm, we can conclude that the TS-BCH has a M better performance compared to the Baesler’s algorithm. ed Table 11. Distributions for Surgery and PACU times of 23 surgery types pt Surger PACU time Surgery Surgery time PACU time Surgery time [min] y Type [min] Type [min] [min] 1 TRIA(25, 57.7, 71) 55+EXPO(17.8) 13 5+EXPO(23.6) 5+EXPO(25) ce 2 TRIA(30, 60, 150) NORM(67.5, 24.6) 14 20+EXPO(9.51) 5+WEIB(3.85, 0.241) 10+61*BETA(0.98, 3 NORM(65.6, 26.3) 15 TRIA(8, 13.7, 41) WEIB(10.5, 0.491) 0.641) 4 45+WEIB(47.8, 1.29) TRIA(10, 84, 115) 16 TRIA(10, 14, 76) 5+WEIB(0.687, 0.256) Ac 5 5+WEIB(26.4, 1.22) 5+EXPO(21) 17 TRIA(30, 46.5, 96) 30+EXPO(55) 30+155*BETA(0.633, 6 35+EXPO(67.2) 18 20+EXPO(9.51) 5+WEIB(4.01, 0.244) 0.917) 56*BETA(0.526, 7 UNIF(24.5, 70.5) TRIA(5, 82.7, 105) 19 5+EXPO(22.5) 0.946) 20+81*BETA(0.899, 8 NORM(93.3, 26.4) NORM(105, 32.3) 20 TRIA(19.5, 43.5, 101) 1.26) 9 5+WEIB(59.6, 1.72) NORM(101, 42.1) 21 10+EXPO(52) TRIA(5, 68, 95) 40+96*BETA(0.693, 15+140*BETA(0.494, 10+121*BETA(1.44, 10 22 TRIA(35, 53, 71) 1.16) 0.481) 0.825) 10+306*BETA(1.66, 11 NORM(107, 47.4) 23 25+EXPO(24.1) 5+WEIB(40, 0.573) 5.12) 12 25+EXPO(36) 10+WEIB(9.27, 0.295) 30 5.3 Comparing the results of TS-BCH algorithm with CPH and AST methods The TS-BCH algorithm presented in the previous subsection uses CPH distribution to compute the probability of completing the surgery before the set threshold T, P (T S  T ) . To show the efficiency of applying CPH distribution within the TS-BCH algorithm, we compare the results with usual methods used in hospitals such as AST. In both methods, we find the best sequence of patients by applying our proposed search algorithm (Figure 9) and schedule them for all ORs and required days. In order to compare CPH and AST methods, we need a method to show the superiority one against another. We choose MC simulation, presented in Section ‎2.3, as a method to evaluate CPH and AST t methods. To this end, we use one sample from each category of defined test problem: small i.e., 50 ip patients, medium i.e., 150 patients and large i.e., 250 patients. We also consider 23 surgery types cr based on Baesler et al. (2015), 3 operating rooms, 6 recovery beds and 8 hours available OR time for each day, from 8:00 AM to 16:00 PM. us The results of TS-BCH algorithm with CPH and AST methods for three different experiments are an depicted in Table 12 to Table 14. In these tables, "the probabilities of completing the surgery before  480) ) in each day", "mean of probabilities" and "a number of days M the defined threshold time ( P(TS all patients scheduled in three ORs" are shown. In Table 12, we first find the best sequence of 50 ed patients by TS-BCH algorithm and schedule them by CPH or AST method for three ORs and required days. We use these results as an input for the MC simulation and determine the P(TS  480) for each pt OR and day. For example, the probability of 0.97 in the first day and first OR of CPH method shows ce the P(TS  480) and 0.93 is the mean of P(TS  480) for all ORs and days in CPH method. P(TS  480) and Ac In almost all results depicted in Table 12 to Table 14, the values of mean of P(TS  480) for AST method are less than 0.85 while they are more than 0.85 for CPH one. Therefore, considering the variability and skewness inherent in the OTR times, the AST method renders an infeasible schedule with respect to the probability of completion before defined threshold time in almost all cases. Table 12. Probability of completing the surgery before the threshold- 50 patients Method ORs Day 1 Day 2 Day 3 Mean 31 OR1 0.97 0.90 1.00 CPH OR2 0.88 0.87 1.00 0.93 OR3 0.85 0.88 1.00 OR1 0.59 0.86 - AST OR2 0.65 0.79 - 0.71 OR3 0.58 0.77 - Table 13. Probability of completing the surgery before the threshold- 150 patients Method ORs Day 1 Day 2 Day 3 Day 4 Day 5 Day 6 Day 7 Mean OR1 0.86 0.92 0.91 0.95 0.90 0.92 0.97 CPH OR2 0.92 0.85 0.97 0.95 0.92 0.93 0.97 0.93 OR3 0.89 0.91 0.92 0.95 0.88 0.98 0.98 OR1 0.57 0.51 0.69 0.63 0.65 1.00 - AST OR2 0.58 0.68 0.66 0.58 0.59 1.00 - 0.67 OR3 0.54 0.64 0.64 0.58 0.60 1.00 - Table 14. Probability of completing the surgery before the threshold- 250 patients ORs Day Day Day Day Day Day Day Day Day Day Day Day Method Mean 1 2 3 4 5 6 7 8 9 10 11 12 OR1 0.89 0.88 0.95 0.94 0.91 0.83 0.90 0.93 0.92 0.90 0.90 1.00 t CPH OR2 0.89 0.90 0.93 0.97 0.92 0.87 0.89 0.93 0.92 0.89 0.87 1.00 0.92 ip OR3 0.86 0.98 0.94 0.97 0.93 0.88 0.89 0.97 0.92 0.89 0.91 1.00 OR1 0.70 0.62 0.59 0.59 0.55 0.64 0.57 0.61 0.58 1.00 - - AST OR2 0.70 0.60 0.61 0.54 0.56 0.66 0.57 0.56 0.57 1.00 - - 0.64 OR3 0.56 0.58 0.63 0.55 0.55 0.67 0.58 0.64 0.54 1.00 - - cr 6 Conclusions and Suggestions for Future Research us To the best of our knowledge, this is the first research on operating theater room (OTR) problem that an applies continuous phase-type (CPH) distributions in order to measure the probability of completion time of surgery before a given threshold. The OTR problem is more challenging because of M considering multiple stages, coordination of many resources, and uncertainty inherent to surgical services. Our main contribution in this research is proposing a stochastic modeling of surgery and ed PACU times with CPH distribution for the first time. We present a model for the OTR scheduling pt problem based on the formalism of the CPH distributions, using the matrix-geometric properties of the CPH distributions to evaluate the convolutions that naturally appear in summing the stochastic ce time of consecutive stages in the OTR. Ac We demonstrate that the CPH distribution outperforms the commonly used distributions such as normal and lognormal for this purpose. We also show that the CPH distribution is the most attractive choice against Monte Carlo (MC) simulation in terms of precision and ease of implementation in the OTR scheduling problems. We develop a search algorithm consisting of constructive heuristic and meta-heuristic algorithms as the patient sequence generator engine and embed in it the CPH evaluator as a chance constraint. The results of experiments for the search algorithm with constructive heuristic show that the CDS and 32 LEOPT algorithms outperform other heuristics on 1 to 200 patients. The LESPT and LEOPT methods give the same results in 200 to 500 patients and both of them perform better than the other heuristics for the mentioned number of patients. We propose four meta-heuristic algorithms based on SA and TS algorithm to consider as a sequence generator engine for the search algorithm. The results of experiments demonstrate that the TS algorithm considering CDS and LESPT methods as a neighborhood search provides better performance than the other proposed meta-heuristic algorithms. We also compare the best meta- heuristic algorithm (TS-BCH) with Baesler’s algorithm, the only algorithm similar to our research problem that we could find in the literature. We show that the TS-BCH has a better performance t ip compared to the Baesler’s algorithm. cr We also show the performance of our proposed search algorithm by comparing the CPH and common methods used in hospitals such as AST to determine the schedule of patients. We conclude that the us AST method renders an infeasible schedule because of disregarding the inherent variability of the an surgery durations and violating the completion time threshold in most of the cases, whereas the CPH method produces feasible schedules that meet the completion time threshold with the designated M probability. Future research can focus on considering other optimization criteria such as minimization of the ed tardiness, waiting time and overtime of patients. Extending the proposed approach to OTR problem pt with downstream resources such as CCU and ICU can be developed for the proposed research problem. ce Acknowledgment Ac The authors thankfully acknowledge professor John Bowers from University of Stirling, for sharing the datasets of Scottish NHS hospital from 1998 to 1999. 7 References 33 Abdelrasol, Z. Y., Harraz, N., and Eltawil, A. (2013). A proposed solution framework for the operating room scheduling problems. In Proceedings of the world congress on engineering and computer science (Vol. 2, pp. 23–25). Abe, T. K., Beamon, B. M., Storch, R. L., and Agus, J. (2016). Operations research applications in hospital operations: Part III. IIE Transactions on Healthcare Systems Engineering, 6(3), 175– 191. Arnaout, J.-P. M., and Kulbashian, S. (2008). Maximizing the utilization of operating rooms with t stochastic times using simulation. In Proceedings of the 40th conference on winter simulation ip (pp. 1617–1623). cr Baesler, F., Gatica, J., and Correa, R. (2015). Simulation optimisation for operating room scheduling. Int J Simul Model (IJSIMM), 14(2), 215–226. us an Baker, K. R., and Trietsch, D. (2013). Principles of sequencing and scheduling. John Wiley & Sons. M Bandyopadhyay, S., and Bhattacharya, R. (2014). Discrete and Continuous Simulation: Theory and Practice. CRC Press. ed Bobbio, A., Horváth, A., and Telek, M. (2005). Matching three moments with minimal acyclic phase pt type distributions. Stochastic Models, 21(2-3), 303–326. ce Bowers, J., and Mould, G. (2004). Managing uncertainty in orthopaedic trauma theatres. European Journal of Operational Research, 154(3), 599–608. Ac Bruni, M. E., Beraldi, P., and Conforti, D. (2015). A stochastic programming approach for operating theatre scheduling under uncertainty. IMA Journal of Management Mathematics, 26(1), 99–119. Buchholz, P., Kriege, J., and Felko, I. (2014). Input modeling with phase-type distributions and Markov models: theory and applications. Springer. 34 Campbell, H. G., Dudek, R. A., and Smith, M. L. (1970). A heuristic algorithm for the n job, m machine sequencing problem. Management Science, 16(10), B–630. Campbell, M. J., and Swinscow, T. D. V. (2011). Statistics at square one. John Wiley & Sons. Cardoen, B., Demeulemeester, E., and Beliën, J. (2010a). Operating room planning and scheduling: A classification scheme. International Journal of Social Health Information Management, 1(1), 71–83. Cardoen, B., Demeulemeester, E., and Beliën, J. (2010b). Operating room planning and scheduling: A t ip literature review. European Journal of Operational Research, 201(3), 921–932. cr Chernick, M. R. (2008). Bootstrap methods: A guide for practitioners and researchers. Hoboken. Nj: us Wiley. an Cruz, A., Muñoz, A., Zamora, J. L., and Espínola, R. (2011). The effect of wind generation and weekday on Spanish electricity spot price forecasting. Electric Power Systems Research, 81(10), M 1924–1935. ed Da Silva Godinho, C. A. P. (2014). Optimizing Operating Theater Planning-A Data Mining And Optimization Appoach. pt Dekhici, L., and Belkadi, K. (2014). Bi-objective Operating Theater scheduling case of the paediatric ce hospital of Oran. In Logistics and Operations Management (GOL), 2014 International Conference on (pp. 181–187). Ac Denton, B., Viapiano, J., and Vogl, A. (2007). Optimization of surgery sequencing and scheduling decisions under uncertainty. Health Care Management Science, 10(1), 13–24. Dexter, F., Macario, A., Traub, R. D., Hopwood, M., and Lubarsky, D. A. (1999). An operating room scheduling strategy to maximize the use of operating room block time: computer simulation of 35 patient scheduling and survey of patients’ preferences for surgical waiting time. Anesthesia & Analgesia, 89(1), 7–20. Epprecht, C., Guegan, D., Veiga, Á., da Rosa, J. C., and others. (2017). Variable selection and forecasting via automated methods for linear models: LASSO/adaLASSO and Autometrics. Erdogan, S. A., and Denton, B. T. (2011). Surgery planning and scheduling. Wiley Encyclopedia of Operations Research and Management Science. Ferrand, Y. B., Magazine, M. J., and Rao, U. S. (2014). Managing operating room efficiency and t ip responsiveness for emergency and elective surgeries—a literature survey. IIE Transactions on Healthcare Systems Engineering, 4(1), 49–64. cr us Glover, F. (1989). Tabu search—part I. ORSA Journal on Computing, 1(3), 190–206. an Glover, F. (1990). Tabu search—part II. ORSA Journal on Computing, 2(1), 4–32. Gómez, A., Mariño, R., Akhavan-Tabatabaei, R., Medaglia, A. L., and Mendoza, J. E. (2015). On M modeling stochastic travel and service times in vehicle routing. Transportation Science, 627 – ed 641. Guerriero, F., and Guido, R. (2011). Operational research in the management of the operating theatre: pt a survey. Health Care Management Science, 14(1), 89–114. ce Hans, E., Wullink, G., Van Houdenhoven, M., and Kazemier, G. (2008). Robust surgery loading. Ac European Journal of Operational Research, 185(3), 1038–1050. Joustra, P., Meester, R., and van Ophem, H. (2013). Can statisticians beat surgeons at the planning of operations? Empirical Economics, 44(3), 1697–1718. Kirkpatrick, S., Gelatt, C. D., Vecchi, M. P., and others. (1983). Optimization by simulated annealing. Science, 220(4598), 671–680. 36 Koulamas, C., and Kyparisis, G. J. (2000). Asymptotically optimal linear time algorithms for two- stage and three-stage flexible flow shops. Naval Research Logistics (NRL), 47(3), 259–268. Lamiri, M., Dreo, J., and Xie, X. (2007). Operating room planning with random surgery times. In Automation Science and Engineering, 2007. CASE 2007. IEEE International Conference on (pp. 521–526). Latorre-Núñez, G., Lüer-Villagra, A., Marianov, V., Obreque, C., Ramis, F., and Neriz, L. (2016). Scheduling operating rooms with consideration of all resources, post anesthesia beds and t emergency surgeries. Computers & Industrial Engineering, 97, 248–257. ip Lee, S., and Yih, Y. (2014). Reducing patient-flow delays in surgical suites through determining start- cr times of surgical cases. European Journal of Operational Research, 238(2), 620–629. us Mancilla, C., and Storer, R. (2012). A sample average approximation approach to stochastic an appointment sequencing and scheduling. IIE Transactions, 44(8), 655–670. M Marcon, E., and Dexter, F. (2006). Impact of surgical sequencing on post anesthesia care unit staffing. Health Care Management Science, 9(1), 87–98. ed Master, N., Scheinker, D., and Bambos, N. (2016). Predicting pediatric surgical durations. arXiv pt Preprint arXiv:1605.04574. ce MATLAB Release 2014a. (2015). MathWorks, Inc. Natick, Massachusetts, USA. Ac Moschopoulos, P. G. (1985). The distribution of the sum of independent gamma random variables. Annals of the Institute of Statistical Mathematics, 37(1), 541–544. Nawaz, M., Enscore, E. E., and Ham, I. (1983). A heuristic algorithm for the m-machine, n-job flow- shop sequencing problem. Omega, 11(1), 91–95. 37 Neuts, M. F. (1975). Computational uses of the method of phases in the theory of queues. Computers & Mathematics with Applications, 1(2), 151–166. Neuts, M. F. (1981). Matrix-geometric solutions in stochastic models: an algorithmic approach. Johns Hopkins University, Baltimore. Oberle, W. (2015). Monte Carlo Simulations: Number of Iterations and Accuracy. Perez, J., and Riano, G. (2007). Benchmarking of fitting algorithms for continuous phase-type distributions. Working Paper, COPA Universidad de Los Andes, 1–20. t ip Pham, D.-N., and Klinkert, A. (2008). Surgical case scheduling as a generalized job shop scheduling cr problem. European Journal of Operational Research, 185(3), 1011–1025. us Pinedo, M. L. (2012). Scheduling: theory, algorithms, and systems. Springer Science & Business an Media. Saadouli, H., Jerbi, B., Dammak, A., Masmoudi, L., and Bouaziz, A. (2015). A stochastic M optimization and simulation approach for scheduling operating rooms and recovery beds in an ed orthopedic surgery department. Computers & Industrial Engineering, 80, 72–79. Saeed, F., Gazem, N., Patnaik, S., Balaid, A. S. S., and Mohammed, F. (2017). Recent Trends in pt Information and Communication Technology: Proceedings of the 2nd International Conference ce of Reliable Information and Communication Technology (IRICT 2017) (Vol. 5). Springer. Ac Samudra, M., Van Riet, C., Demeulemeester, E., Cardoen, B., Vansteenkiste, N., and Rademakers, F. E. (2016). Scheduling operating rooms: achievements, challenges and pitfalls. Journal of Scheduling, 19(5), 493–525. SAS Release, 9.1. (2003). SAS Institute Inc. Cary, North Carolina, USA. Retrieved from http://www.sas.com/ 38 Shylo, O. V, Prokopyev, O. A., and Schaefer, A. J. (2012). Stochastic operating room scheduling for high-volume specialties under block booking. INFORMS Journal on Computing, 25(4), 682– 692. Sier, D., Tobin, P., and McGurk, C. (1997). Scheduling surgical procedures. Journal of the Operational Research Society, 48(9), 884–891. Stepaniak, P. S., Heij, C., and De Vries, G. (2010). Modeling and prediction of surgical procedure times. Statistica Neerlandica, 64(1), 1–18. t ip Stepaniak, P. S., Heij, C., Mannaerts, G. H. H., de Quelerij, M., and de Vries, G. (2009). Modeling procedure and surgical times for current procedural terminology-anesthesia-surgeon cr combinations and evaluation in terms of case-duration prediction and operating room efficiency: us a multicenter study. Anesthesia & Analgesia, 109(4), 1232–1245. an Tànfani, E., and Testi, A. (2012). Advanced Decision Making Methods Applied to Health Care (Vol. 173). Springer Science & Business Media. M Thümmler, A., Buchholz, P., and Telek, M. (2006). A novel approach for phase-type fitting with the ed EM algorithm. IEEE Transactions on Dependable and Secure Computing, 3(3), 245–258. pt Wang, Y., Tang, J., Pan, Z., and Yan, C. (2015). Particle swarm optimization-based planning and ce scheduling for a laminar-flow operating room with downstream resources. Soft Computing, 19(10), 2913–2926. Ac Xiang, W., Yin, J., and Lim, G. (2015). A short-term operating room surgery scheduling problem integrating multiple nurses roster constraints. Artificial Intelligence in Medicine, 63(2), 91–106. Zhao, Z., and Li, X. (2014). Scheduling elective surgeries with sequence-dependent setup times to multiple operating rooms using constraint programming. Operations Research for Health Care, 3(3), 160–167. 39 Appendix A. Characterization of Continuous Phase-Type Distribution and Its Properties A.1. Characterization The distribution of time Z until the process reaches the absorbing state is said to be CPH distributed and is denoted Z ~ PH C (, T) . The cumulative distribution function of the CPH distribution Z ~ t PH C (, T) is ip FZ (x )  1  e Tx 1 for x  0 cr (A.1) the density function is where e Tx f Z (x )  e Tx t denotes the matrix exponential, defined as: for x  0 us (A.2) an  (Tx )k e Tx  (A.3) k 0 k! M and the ith moments is i  E [X i ]  i !(- T)i 1 . (A.4) ed A.2. Closure properties pt One of the appealing features of CPH distributions is that the class is closed under a number of operations. The closure properties are the main contributing factor to the popularity of these ce distributions in stochastic modeling. Ac Assume that Z i ~ PH C ((i ) , T(i ) ) for i = 1, 2 are two independent CPH distributed random variables of order ni . (i) Convolution of PHC: the sum Z  Z 1  Z 2 ~ PH C (, T) has a CPH distribution of order n  n1  n2 with representation 40   (1) t (1)  (2)    ((1) , 0(1) (2) ) and   . (A.5)  0  (2)  Proof. See Neuts (1981), Theorem 2.2.2. (ii) Mixture of PHC: the convex mixture Z   Z 1  (1   )Z 2 ~ PH C (, T) has a CPH distribution of order n  n1  n2 with representation   (1) 0    ((1) ,(1   )(2) ) and   . (A.6)  0  (2)  Proof. See Neuts (1981), Theorem 2.2.4. (iii) Minimum of PHC: The minimum Z  min(Z 1 , Z 2 ) ~ PH C (, T) has a CPH distribution t ip of order n  n1.n2 with representation cr   (1)  (2) and   (1)  (2) , (A.7) where  is the Kronecker product and  is Kronecker sum. Proof. See Neuts (1981), Theorem 2.2.9. us an (iv) Maximum of PHC: The maximum Z  max(Z 1 , Z 2 ) ~ PH C (, T) has a CPH distribution of M order n  n1.n2  n1  n2 with representation ed   (1)   (2)  (1)  t (2) t (1)   (2)      ((1)  (2) , 0(2) (1) , 0(1) (2) ) and    0  (1) 0 . (A.8)  0 0  (2)    pt Proof. See Neuts (1981), Theorem 2.2.9. ce Appendix B. Parameter tuning for meta-heuristic algorithms Ac Table B1. Simulated annealing parameter Initial temperature (T 0 ) 1000 Cooling temperature ( ) 1 Table B2. Tabu Search parameter Tabu list size Number of patients From To Parameter formula  Number of Patients  1 100  10  41 101 150 10  Number of Patients  151 500  20  Appendix C. The results of residual plot, equality of variances, ANOVA table and Tukey test for the OTR problem t ip cr us an Figure C1. Residual plot for four meta- heuristic algorithms Table C1. Test for Equal Variances: SA-PE, SA- BCH, TS-PE, TS- BCH Method M Null hypothesis All variances are equal Alternative hypothesis At least one variance is different Significance level α = 0.05 ed Bartlett’s method is used. This method is accurate for normal data only. Table C2. 95% Bonferroni Confidence Intervals for Standard Deviations pt Sample N StDev CI SA-PE 3999 10180.3 (9903.20, 10472.3) ce SA- BCH 3999 10196.3 (9918.75, 10488.7) TS-PE 3999 10133.3 (9857.50, 10423.9) TS- BCH 3999 10079.0 (9804.71, 10368.1) Ac Individual confidence level = 98.75% Table C3. Result of Test for Equal Variances: four meta-heuristic algorithms Method Test Statistic P-Value Bartlett 0.65 0.885 Table C4. The ANOVA table for four meta-heuristic algorithms Source DF Sum of squares Mean squares F-value Pr > F Dependent variable: Y Model 2002 820862591684 410021274 4831.86 <.0001 Error 5997 508892986 84857 42 Corrected total 7999 821371484671 R-Square Coeff Var Root MSE Objective Mean 0.999380 1.703614 291.3038 17099.17 Source DF Type III SS Mean square F-value Pr > F Method 3 47058377 15686125 184.85 <.0001 Block 1999 820815533307 410613073 4838.83 <.0001 Table C5. Tukey test for four meta-heuristic algorithms. Differences of least squares means. Algorithms Comparison* Difference Between Means 95% Confidence Limits Algorithms Comparison* 1-2 11.37 [-6.68 29.42] 1-3 50.51 [32.45 68.56] *** 1-4 192.38 [174.32 210.44] *** 2-1 -11.37 [-29.42 6.68] 2-3 39.14 [21.08 57.19] *** 2-4 181.01 [162.95 199.07] *** 3-1 -50.51 [-68.56 -32.45] *** 3-2 -39.14 [-57.19 -21.08] *** t 3-4 141.87 [123.81 159.93] *** ip 4-1 -192.38 [-210.44 -174.32] *** 4-2 -181.01 [-199.07 -162.95] *** 4-3 -141.87 [-159.93 -123.81] *** cr 1 is the SA-PE algorithm. 2 is the SA- BCH algorithm. 3 is the TS-PE algorithm.4 is the TS- BCH algorithm. us Table C6. The average of the objective function values for the test problems by using four meta-heuristic algorithms Meta-heuristic Algorithm Average SA-PE 17162 an SA- BCH 17151 TS-PE 17112 TS- BCH 16970 M ed Appendix D. The results of residual plot, equality of variances and paired t-test experiment for the TS- pt BCH algorithm and Baesler’s algorithm ce A l g o r i t h m= T S - B CH A l g o r i t h m= B a e s l e r 40000 40000 30000 30000 Ac R R e e s 20000 s u 20000 u l l t t 10000 10000 0 0 - 4 - 3 - 2 - 1 0 1 2 3 4 - 4 - 3 - 2 - 1 0 1 2 3 4 No r ma l Qu a n t i l e s No r ma l Qu a n t i l e s Figure D1. Residual plot for TS-BCH algorithm Figure D2. Residual plot for Baesler’s algorithm Table D1. Test for Equal Variances: TS-BCH algorithm and Baesler’s algorithm Null hypothesis All variances are equal 43 Alternative hypothesis At least one variance is different Significance level α = 0.05 Table D2. 95% Bonferroni Confidence Intervals for Standard Deviations Sample N StDev CI TS-BCH 2000 10095.8 (9874.4, 10333.7) Baesler’s algorithm 2000 10235.6 (10014.2, 10473.6) Individual confidence level = 97.5% Table D3. Result of Test for Equal Variances: TS-BCH algorithm and Baesler’s algorithm Method Test Statistic P-Value Multiple comparisons 0.93 0.334 Levene 0.70 0.404 Table D4. Paired t-test for TS-BCH and Baesler’s algorithm t ip The TTEST Procedure Difference: TS-BCH - Baesler’s algorithm cr N Mean Std Dev Std Err 2000 -228.5 516.36 11.546 Mean -228.5 95% CL Mean -251.2 -205.9 Std Dev 516.36 us 95% CL Std Dev 500.84 532.88 an DF t Value Pr > |t| 1999 -19.79 <.0001 M Table D5. The average of the objective function values for the test problems by using TS-BCH and Baesler’s algorithm Algorithm Average ed TS-BCH 33981005 Baesler’s algorithm 34438089 pt ce Ac 44