- Email: [email protected]
Waiting list management through master surgical schedules: A case study
Operations Research for Health Care 10 (2016) 49–64
Contents lists available at ScienceDirect
Operations Research for Health Care journal homepage: www.elsevier.com/locate/orhc
Waiting list management through master surgical schedules: A case study Belinda Spratt, Erhan Kozan ∗ Queensland University of Technology, 2 George St, Brisbane, QLD 4000, Australia
article
info
Article history: Received 24 January 2016 Accepted 11 July 2016 Available online 15 July 2016 Keywords: MSSP SCAP Waiting list management Operating theatre planning Hybrid metaheuristics
abstract In this paper we address the problem of generating master surgical schedules (MSSs) that adhere to staff and equipment restrictions whilst ensuring patients are treated in a timely manner. We simultaneously address the master surgical scheduling problem (MSSP) and the surgical case assignment problem (SCAP). Stochastic surgical durations are considered in order to produce more robust schedules and reduce unexpected overtime. Also incorporated into the model are several constraints regarding patient wait targets that are set by the Australian government. The problem is formulated using a mixed integer nonlinear programming (MINLP) approach and solved using a variety of hybrid metaheuristics. The metaheuristics implemented are inspired by simulated annealing (SA) and reduced variable neighbourhood search (RVNS). In particular, we present an adaptive SA and hybridise SA and RVNS to greatly improve solution quality. The solution neighbourhoods used by the metaheuristics are based on the hierarchical structure present in the combined MSSP SCAP. We consider a case study of an Australian public hospital with a large surgical department and compare the performance of our model to historical data. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction The surgical department is often one of the largest sources of hospital costs and patient flow. Inefficient use of the operating theatres (OTs) and related resources can result in increased costs and a decrease in patient and staff welfare. The schedule used by hospital administrators can have a massive impact on patient outcomes. Although the scheduling of surgical procedures is a complex process, many surgical departments still rely on simple heuristics and general intuition to produce weekly or fortnightly OT schedules and some daily updates. Such simple techniques can result in overtime, bed shortages, unused resources or large waiting lists due to unnecessary improper utilisation of operating rooms (ORs). The administrative practices of surgical departments can have a large impact on hospital costs, patient outcomes and the overall efficiency of a hospital. As such, a large number of papers exist which investigate and analyse some of the most significant issues in OT planning. For a thorough review of the existing literature see [1–3].
∗
Corresponding author. E-mail address: [email protected] (E. Kozan).
http://dx.doi.org/10.1016/j.orhc.2016.07.002 2211-6923/© 2016 Elsevier Ltd. All rights reserved.
Van Riet and Demeulemeester [1] review efforts to plan ORs in the case of both elective and emergency surgeries with an emphasis on the trade-offs encountered between costs and waiting times. Ferrand et al. [2] also consider both elective and emergency cases in their survey of methodologies used to ensure ORs are both responsive and efficient. In their literature survey, Guerriero and Guido [3] discuss both optimisation and simulation models used for OT planning and scheduling problems. OT planning problems can exist at the strategic, tactical or operational levels. Decisions at the strategic level include determining how much time to allocate to each surgical specialty in order to optimise the mix of patients processed [4]. Tactical problems include the development of a MSS [5], or determining which time blocks to allocate to a particular surgeon [6] or surgical specialty [7]. Problems at the operational level are advanced and allocation scheduling of elective patients. Advanced scheduling is the assignment of a patient to a time block on a particular day. This is also known as the SCAP. Allocation scheduling is the scheduling of patients within those time blocks, also known as the surgical case sequencing problem (SCSP). There are three main policies used in OT scheduling at the tactical level. These are open scheduling, block scheduling and modified block scheduling [6]. An open scheduling policy allows surgeons to assign a surgical case to any day in the planning horizon. Block scheduling allocates time blocks to surgeons or
50
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
surgical specialties. These time blocks may only be used by the surgeon or surgical specialty that they are assigned to. In a modified block scheduling policy, time blocks may be left open for use or unused time can be reallocated to a different surgeon or specialty. A MSS is often used under a block scheduling strategy in order to assign time blocks to surgical specialties. There are many benefits that arise from the implementation of a MSS [8]. The cyclic nature of the schedule allows the coordination of staff and resources and can improve patient flow and throughput. The predictability of a repeated MSS can be incredibly useful when planning admissions and can help stabilise bed occupancy. It is possible to generate an optimal or near optimal MSS for a hospital, based on their particular requirements and goals, using a variety of different operations research techniques. Many hospitals repeat the same MSS over a period of several months. Agnetis et al. [9] analyse the trade-off between the organisational benefits of repeating a single weekly MSS against a schedule updated weekly to adapt to the current state of the surgical waiting lists. The authors consider three policies for the development of MSSs, each permitting different levels of variability from an original schedule. The authors found that allowing even a small number of deviations from a reference schedule greatly improves the performance of the MSS in terms of lateness and waiting time. Throughout the literature, it is common to incorporate constraints on the suitability of the ORs for each specialty, the availability of surgical teams, and post-surgical capacity (e.g. [7]). Other authors include constraints on the patient mix to match hospital requirements. Most of the models used to solve the MSSP have a single objective which is usually related to maximising throughput or utilisation. Some models throughout the literature have multiple objectives, commonly incorporated in a single weighted objective function. Whilst a large number of mathematical formulations and solution techniques exist to produce optimal OT schedules with respect to a variety of objectives, it is necessary to make these solutions applicable and accessible to hospital administration staff. In order to allow surgical staff and scheduling nurses to utilise schedules generated using operations research techniques, a number of decision support systems are presented in the literature (e.g. [10]). The advanced scheduling problem (the OT planning problem or SCAP) is the problem of assigning patients to time blocks over a one to two week time period. This can be done using an open scheduling, block scheduling, or modified block scheduling policy [11]. Advanced scheduling is frequently addressed in the literature and solved using a variety of both exact and heuristic approaches. Whilst some authors rely on exact approaches to solve small test instances, as the problem is NP-hard [12] it is necessary to use heuristics, meta-heuristics and hybrid heuristics to provide good feasible solutions in reasonable amounts of computational time. Column-generation based heuristics are often used to solve OT scheduling problems at the operational level (cf. Wang et al. [13] and Fei et al. [14]) Fei et al. [15] apply a branch-and-price approach to solve the SCAP. Many authors rely on simulation based models to investigate OT planning problems (e.g. [16,17]). The constraints and objective functions used to define OT practice can differ significantly between models. Many models present in the literature minimise expenses related to OT use. Other models prioritise patient welfare over expenses. Roland et al. [18] places particular significance on human resource constraints. The tactical and operational levels of OT planning and scheduling are often treated separately. A more effective procedure involves scheduling the tactical and operational levels
simultaneously. Agnetis et al. [19] integrate the two levels of planning and compare an integrated approach to the use of a decomposition approach. Whilst the decomposition approach allowed the authors to produce comparable solutions in far less time, the model had several simplifying assumptions that make this method unsuitable here. The authors did not consider specific surgeon availability or suitability, and incorporated uncertainty in surgical durations using added slack time in conjunction with expected durations. These assumptions make this approach unsuitable for implementation on the case study considered in Section 2. Aringhieri et al. [20] solve a binary integer programming (BIP) formulation of the MSSP and the advanced scheduling problem (SCAP) simultaneously by making use of the hierarchical structure present in OT schedules in a metaheuristic. We also consider the hierarchical structure when solving the MINLP formulation presented in Section 3. The model presented by Aringhieri et al. [20] considers deterministic surgical durations, whereas we consider stochastic durations. Aringhieri et al. [20] assume that the strategic decisions can be used as an input to the problem, particularly the number of time blocks assigned to each surgical specialty. We implicitly consider the surgical case mix planning problem when solving the combined MSSP SCAP to allow for variations in the MSS based on the current state of waiting lists. One of the major concerns currently facing Australian hospitals is the length of elective surgery waiting lists. The longer a patient waits for surgery, the more their quality of life is reduced. Long surgery waits can also worsen patient outcomes. Occasionally, patients are not seen on time as they are not fit for surgery at that stage, however most overdue patients are late for surgery due to insufficient resources or poor planning. A limited number of operations research techniques have been used to manage waiting lists and improve patient outcomes whilst adhering to resource availability constraints. Everett [21] provides a decision support tool based on a simulation model of elective surgery waiting lists for Australian public hospitals. The model developed by the author can be used in order to monitor the system, assist in strategic decisions or even surgical scheduling at the tactical and operational levels. The simulation model is quite simplistic to allow for the results to be interpreted by hospital staff with ease and does not include as much detail as we incorporate into the model presented in Section 3. Patient priority is considered by Min and Yih [22] whilst assigning elective patients to available surgical capacity. In their model, the authors assign appropriate patients in order of urgency. A common issue identified in waiting list management is the discrepancy in waiting times between urgency categories which often leads to excessive waiting for non-urgent elective patients. Min and Yih [23] use a patient priority which is a combination of urgency category and waiting time by formulating the problem as an infinite horizon Markov decision process. Marques et al. [24] address the SCAP for a small Portuguese public hospital under an open scheduling strategy. The authors aim to maximise surgical suite utilisation and the number of surgeries scheduled. The authors ensure that all deferred urgent and high priority surgeries are scheduled within the planning horizon and use the remainder of the time to schedule priority and normal surgeries according to a weighted objective function. In this paper we provide a MINLP formulation of the combined MSSP SCAP that maximises the number of surgeries scheduled. Whilst the above waiting list management policies are worth considering, in this paper we consider the policies imposed by the Australian government to ensure the model presented is applicable to Australian public hospitals. We ensure that certain government dictated patient-wait targets are satisfied on the final day of the planning horizon and adhere to resource availability and suitability
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
restrictions (including those on surgical staff). We also incorporate the government’s ‘‘treat in turn’’ policy in which 60% of eligible patients are treated in the order in which they arrived on the waiting list. The model presented in Section 3 contributes to the literature in a number of ways. We tackle all three levels of OT planning and scheduling problems by solving the MSSP and SCAP and, in doing so, implicitly solve the surgical case mix planning problem. We also use an integrated approach when solving the combined MSSP SCAP to produce better solutions than those obtained under a decomposition approach. The model we present is more detailed than those seen in the literature through inclusion of waiting list management policies, stochastic surgical durations, surgeon availabilities and capabilities. This allows for its implementation in a real hospital setting. The metaheuristics presented in Section 4 can be used to produce good feasible solutions to the combined MSSP SCAP in a reasonable amount of computational time. In particular, hybridised SA–RVNS has not previously been used on hospital planning problems. We are also able to maintain feasibility of a highly non-linear formulation through use of the solution neighbourhoods presented in Section 4.1. The surgical department considered in the case study is one of the largest in Australia. Due to the size of the case study, the real-life instances solved in Section 5 are far more computationally complex than the majority of instances previously considered in OT planning and scheduling. In Section 2 we present the case study, including details on the surgical department, length of surgery, recent waiting list data, and current scheduling practice. In Section 3 we present a MINLP formulation of the problem described within the case study. We discuss the metaheuristics considered, solution neighbourhoods, and upper bounds in Section 4. Computational results including metaheuristic parameter tuning and model performance are presented in Section 5. In Section 6 we present conclusions and discuss possible future work. 2. Case study In this section we consider the case study of a large Australian public hospital. This adult hospital has one of the largest surgical departments of public Australian hospital and as such, it is incredibly important that the surgical department runs efficiently. The throughput of the hospital is bound by a number of fixed resources. The hospital has 21 ORs, with one reserved for emergency surgeries. Elective surgeries are only scheduled on weekdays and never after 8pm. On the weekend, five emergency ORs are open. There are around 800 beds in the hospital, including 300 surgical beds. The surgical care unit (SCU) has around twenty beds available, as does the post anaesthesia care unit (PACU). As the hospital is a teaching hospital it is necessary to have surgeons of an appropriate skill level available to supervise the surgeries. Variable resources include the availability of surgeons, nurses, anaesthetists, etc. We only consider the availability of the assigned surgeons as there are sufficient numbers of nurses and anaesthetists available. The scheduling of nurses and anaesthetists is performed based on the MSS, but not simultaneously. The surgical nurses often specialise, whereas the anaesthetists are more generalised, although they may have specialty preferences. For a patient to enter the elective surgery waiting list, certain procedures must be followed. Initially, a patient must be referred to a specialist. After the patient has had a consultation the case is reviewed by the surgeon to ensure the patient does require the surgery. If these conditions have been met, the patient is put on the elective surgery waiting list and assigned an urgency category. There are three urgency categories used for classifying elective patients. A patient’s category is determined by the recommended
51
time before his/her elective procedure (30, 90 and 365 days). Another measure that the Queensland government has taken in order to effectively manage waiting lists is the use of a ‘‘treat in turn’’ model for patient care. The treat in turn model ensures that at least 60% of eligible patients treated in the planning horizon are among the set of patients who have been waiting the longest. A patient is treat in turn eligible if he/she is in category two or three (less urgent than a category one patient), a public patient (private hospitals do not have to adhere to this policy), and ready for care at the start of the planning horizon. The set of all treat in turn eligible patients is known as the treat in turn cohort. The following method is used by Queensland Health [25] to determine the proportion of patients treated in turn. First, the longest waiting cohort is calculated as the subset of treat in turn eligible patients who have been waiting longest. The number of patients in the longest waiting cohort is equal to the number of treat in turn eligible patients treated during the planning horizon. The treat in turn policy is satisfied as long as the number of longwait patients treated during the planning horizon is at least 60% of the number of treat in turn eligible patients treated, regardless of sequence within the planning horizon. Another policy implemented by the hospital is that all category one patients must be given a surgery date at the time they are added to the waiting list. To manage this, each surgical specialty has a case manager who liaises with surgical staff to ensure that patients are seen in a timely manner. Currently there are around 2900 patients on the elective surgery waiting list. A breakdown of a recent waiting list can be seen in Table 1. Staff at the elective surgery booking office, in conjunction with surgeon preferences, assign elective booking dates to patients whilst satisfying the MSS, patient wait targets and the treat in turn model. Currently the hospital operates under a four week rotating MSS with two blocks per OR per day (or a full day block). The use of a MSS ensures that there is consistency in staff schedules and predictable patient flow patterns for the downstream wards. Very rarely are major changes made to the MSS. The schedule has remained largely the same since around 2000. Specialties are able to request extra OR time in order to keep up with changes in demand. If an OR is available a cost analysis is performed and, if approved, the schedule is updated. In many cases, if there is a change in demand for a particular specialty it may mean that another specialty must reduce their hours or the extra surgeries may be outsourced to hospitals in the surrounding region. Minor modifications are often made as surgeon availabilities change. When updating or modifying the MSS, preference is given to full day blocks over half day blocks as there is no need to stop half way through the day. The length of time allocated to each time block is dependent on the specialty that will be using the OR and the availability of visiting medical officers (VMOs). Some specialties require longer working days. Cardiology often requires ten hour days whereas ENT may require a twelve to fourteen hour block for a major reconstructive case. Each Thursday the hospital has an elective bookings meeting to discuss the patients planned for the following week. Some surgeons request that a time block be left empty as they predict that there will be admissions over the weekend that require OR time. If a surgical team is unavailable or does not need their usual OR time, other specialties are able to request the OR during that time. VMOs are also able to request OR time. One difficulty in creating the MSS is the availability of consultants. Many consultants spend the majority of their time working at surrounding private hospitals and are only available for a few days each month. It is important to be able to incorporate the required OR time for each of the consultants whilst satisfying their availabilities. Some specialties do not require that a consultant is present for surgeries, and instead a registrar can perform a list of surgeries for several consultants.
52
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 1 Recent elective surgery waiting list. Urgency category
1
Specialty
On time
Overdue
2 On time
Overdue
On time
Overdue
Cardiology (CARD) Ear, Nose and Throat (ENT) General Surgery (GEN) Gynaecology (GYN) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Urology (UROL) Vascular Surgery (VASC) Total
4 10 80 0 12 12 22 58 58 23 279
0 0 0 0 0 0 0 0 0 0 0
46 46 183 1 79 194 163 213 143 70 1138
0 17 29 0 10 10 26 50 2 5 149
1 12 74 0 4 760 310 70 118 13 1362
0 0 1 0 0 6 2 3 3 1 16
Leading up to surgery, there are a number of patient pathways. The patient must attend the preadmission clinic at the hospital before his/her surgical date. On the day of surgery, elective patients arrive at the hospital and inpatients are allocated to a ward. Depending on the time until surgery and any other procedures required beforehand, the patient is either sent to the ward or to the SCU to prepare for surgery. After surgery, all patients are sent to the PACU to recover; this usually takes four to five hours. After the initial recovery period, day patients recover further in the SCU and inpatients are sent back to the ward or intensive care unit (ICU). Any surgical transfers from the ICU are returned directly to ICU after surgery. Whilst there are a limited number of beds in the SCU and PACU, as the surgical department has priority in booking these beds, there are rarely bottlenecks. Occasionally downstream wards are put under pressure due to the large number of surgical patients being transferred. Although downstream capacity is out of the scope of this project, it would be extremely useful to incorporate into future work. The hospital has over 50 surgical specialties, however when surgeries are booked they are usually booked under nine to twelve main specialties. Table 2 shows the count, mean, and variance of the surgical durations (in min) for elective surgeries performed over 2012 and 2013. More recent data is unavailable. From Table 2 it can be seen that specialties that perform more complex procedures have a larger variance than specialties that perform less complicated procedures. Whilst nine to twelve specialties are usually used when booking at this hospital; surgical data is reported under other specialties as well. This difference in reporting systems is limiting as we are only able to schedule under the specialties that the elective surgery booking office uses on the waiting lists. The MSS used by the hospital is able to be more specific. Given that the hospital under study is a teaching hospital, there tend to be variations in surgical durations based on the experience of the surgeon or the severity of the case. When surgeries run overtime it can result in large overtime costs or even the cancellation of surgeries. Whilst surgeons can provide estimate surgical durations on a case by case basis this data is not available for use in the model. Kayış et al. [26] present a statistical model to estimate surgical durations based procedure type, staff, patient type and when the procedure occurs. The authors found that by hybridising their work with that of Dexter et al. [27] they are able to predict surgical durations without bias and with a high accuracy. Due to a lack of appropriate data, we must consider a more simple approximation when incorporating stochasticity in the model presented in Section 4. Strum et al. [28] showed that surgical durations can be modelled using a lognormal distribution. Stuart and Kozan [29] use a lognormal approximation to the sum of lognormally distributed
3
random variables in order to ensure a 15.87% level of accuracy. The authors consider only a single OR. We tighten that bound to ensure that there is a probability of less than 5% that the surgeries assigned to a particular block run overtime and extend the work by considering multiple ORs. The lack of stochasticity is not a limitation in the case of non-elective surgeries as there is an OR reserved for these cases, although these cases often overflow into underutilised OR blocks. Non-elective surgeries consist of both emergency arrivals and other urgent surgeries from within the hospital. Heng and Wright [30] found that the use of a dedicated OR significantly reduces the number of elective cancellations and improves the waiting times of non-elective cases. Like many hospitals, certain ORs are more suited for use by certain specialties. For example, there is extra equipment in the vascular OR; however this does not prohibit other specialties from using the OR. It is a matter of convenience that the hospital usually assigns certain specialties to each OR. This reduces the movement of valuable equipment between sessions. Another consideration that administrative staff make when adjusting the MSS is the availability of equipment. The hospital has a single machine that can be used by several specialties, it is important that two specialties are not simultaneously scheduled to use the same equipment. Administrative staff also attempt to have the equipment used as frequently as possible. Fig. 1 (see Appendix) shows the surgeries performed in the twenty ORs under study in the week from which the elective surgery waiting list was used in Section 5.1. The figure shows the specialties that used the time block, the number of elective and emergency surgeries performed, and the total time required for the surgeries. For example, in OR one on the Monday morning there was a single elective surgery, no emergency surgeries, and the total time was 3.08 h. In this week, a total of 266 surgeries were performed of which 218 were elective surgeries. Future work includes incorporating emergency arrivals in the model as the hospital is unable to complete all emergency procedures in the one reserved OR. To test whether the assumption of lognormal surgical durations holds, the data is first transformed by taking the natural logarithm of the surgical durations. We then test whether the transformed data is likely to belong to a normal distribution. R© version 3.1.1 is used during statistical analysis. Data is visually analysed using density plots and normal Q–Q plots. Shapiro–Wilk, Anderson–Darling and Lilliefors normality tests are used to determine the suitability of the distribution. The results from these tests are shown in Table 3. From Table 3 it can be seen that the assumption of lognormal surgical durations is inappropriate for most specialties according to the goodness-of-fit tests used. This may be due to the large (n > 100) and heterogeneous populations. Within the model, the lognormal approximation is only used to determine how many patients can be scheduled within a
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
53
Table 2 Reported elective length of surgery data. Specialty
2012 Time (min) Count
Average
Variance
Count
Average
Variance
Acute Surgical Unit (ASU) Breast and Endocrine (B&E) Cardiology (CARD) Colorectal (COLO) Cardiac Surgical Unit (CSU) Dental Surgery (DENT) Ear, Nose and Throat (ENT) Faciomaxillary (FMAX) Gastroenterology (GAS) Gynaecology (GYN) Hepato-Pancreato-Biliary (HPB) Liver Transplant (LTPT) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Persistent Pain Management Service (PPMS) Psychiatry (PSY) Respiratory (RESP) Renal Transplant (RTPT) Upper GI and Soft Tissue (UGI) Urology (UROL) Vascular Surgery (VASC)
–
– 122.31 0.00 290.91 260.41 75.92 336.44 150.29 – 34.34 223.10 53.03 214.85 54.70 172.70 129.90 – – – 114.72 149.79 95.28 179.64
– 2996.82 0.00 55 744.56 7426.43 76.93 105 600.21 16 927.13 – NA 14 516.63 848.91 15 944.61 1477.97 15 385.07 21 533.34 – – – 5467.35 15 159.66 6826.57 8819.81
23 540 – 359 843 4 471 338 2 4 355 6 451 1988 1735 1622 16 130 91 572 518 1472 528
136.31 108.54 – 173.19 266.83 71.46 254.35 99.13 53.53 33.84 182.60 64.81 190.77 47.55 130.59 85.77 34.46 13.51 45.41 93.76 133.20 86.53 152.31
4446.71 2709.94 – 20 950.38 8375.50 1486.97 76 440.92 5574.57 72.72 36.03 12 292.92 900.69 15 277.17 1318.96 7559.47 7736.56 327.36 12.56 582.69 3502.18 12 647.82 6898.36 7149.21
422 88 524 1163 6 806 228 – 1 484 2 610 1941 1507 1933 – – – 401 541 1884 437
2013 Time (min)
Table 3 Statistical analysis on surgical durations. Specialty
Observations
S.W. p-val
A.D. p-val
Lillie p-val
Acute Surgical Unit (ASU) Breast and Endocrine (B&E) Colorectal (COLO) Cardiac Surgical Unit (CSU) Dental Surgery (DENT) Ear, Nose and Throat 1 (ENT1) Ear, Nose and Throat 2 (ENT2) Ear, Nose and Throat 3 (ENT3) Faciomaxillary (FMAX) Gastroenterology (GAS) General (GEN) Gynaecology (GYN) Hepato-Pancreato-Biliary (HPB) Liver Transplant (LTPT) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Renal Transplant (RTPT) Upper GI and Soft Tissue (UGI) Urology (UROL) Vascular Surgery (VASC)
23 540 359 843 4 34 29 3 338 2 1740 4 355 6 451 1998 1735 1622 572 518 1472 528
0.517 0.370 0.000 0.000 0.323 0.013 0.283 0.400 0.061 NA 0.000 0.864 0.000 0.830 0.001 0.000 0.000 0.000 0.001 0.000 0.000 0.005
0.352 0.131 0.000 0.000 NA 0.014 0.328 NA 0.139 NA 0.000 NA 0.000 NA 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.730 0.063 0.000 0.000 NA 0.037 0.365 NA 0.259 NA 0.000 NA 0.000 0.684 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001
time block such that in 95% of cases the total duration of their surgeries does not exceed the length of the time block. Although the lognormal distribution appears to be a poor fit to the model, random sampling and simulation from actual data can be used to verify the maximum number of patients that can be scheduled in each time block, under each specialty. We verify the suitability of the lognormal approximation by randomly sampling (according to a uniform distribution) 50,000 sets of patients for each combination of set size (between one and twenty patients) and specialty. The sum of surgical durations is calculated for each of the sets and the 95th percentile is found. Based on this we then determined the number of patients that can be scheduled in each type of time block. Table 4 shows the maximum number of patients that can be scheduled in each time block under each specialty as well as the lognormal distribution parameters calculated using historical data. Note that an asterisk beside a surgical specialty indicates that the lognormal distribution parameters from the General
specialty were used as there was insufficient data. The values within this table are the same whether calculated using lognormal approximations or random sampling. Thus, we are confident in the use of lognormal approximations for surgical duration in the model, despite poor goodness-of-fit test results. 3. Model formulation The combined MSSP SCAP is formulated using MINLP to allocate surgeries, surgeons and specialties to time blocks. 3.1. Scalar parameters Bmax : number of time blocks in a working day Cmax : number of patient urgency categories Dmax : number of working days in the planning horizon Hmax : number of surgeons that practice at the hospital (including VMOs)
54
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 4 Surgical duration model parameters. Specialty *
Anaesthetics (ANAE) Acute Surgical Unit (ASU) Breast and Endocrine (B&E) Cardiothoracic (CARD)* Colorectal (COLO) Cardiac Surgical Unit (CSU) Dental Surgery (DENT) Ear, Nose and Throat 1 (ENT1) Ear, Nose and Throat 2 (ENT2) Ear, Nose and Throat 3 (ENT3) Faciomaxillary (FMAX) Gastroenterology (GAS) General (GEN) Gynaecology (GYN) Hepato-Pancreato-Biliary (HPB) Liver Transplant (LTPT) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Renal Transplant (RTPT) Trauma (TRMA)* Upper GI and Soft Tissue (UGI) Urology (UROL) Vascular Surgery (VASC) Visiting Medical Officers (VMOs)* *
Mu
Sigma2
Patients per half day
Patients per full day
0.788 0.702 0.478 0.788 0.783 1.426 0.036 1.200 0.498 0.585 0.266 −0.137 0.788 −0.598 0.945 −0.031 0.970 −0.474 0.583 −0.015 0.267 0.788 0.516 0.027 0.786 0.788
0.395 0.217 0.209 0.395 0.534 0.112 0.258 0.513 0.779 0.182 0.453 0.025 0.395 0.031 0.317 0.196 0.353 0.463 0.370 0.724 0.338 0.395 0.542 0.658 0.271 0.395
0 1 1 0 0 0 2 0 0 1 1 5 0 8 0 3 0 3 1 1 1 0 0 1 0 0
1 2 3 1 1 1 6 0 1 3 3 10 1 16 1 7 1 8 2 3 4 1 2 3 2 1
Indicates lognormal distribution parameters from the General specialty were used as there was insufficient data.
Pmax : number of patients on the waiting list at the start of the planning horizon Smax : number of surgical specialties (where surgeons of specialty Smax are VMOs) Tmax : number of operating rooms M: the largest possible number of days overdue for any patient at the end of the planning horizon Mp : the maximum possible number of patients that can be scheduled in a full day time block π : the proportion of eligible patients that must be treated in turn. 3.2. Index sets B: the set of time blocks in a working day. B = {1, . . . , Bmax } C : the set of patient urgency categories. C = {1, . . . , Cmax } D: the set of working days in the time horizon. D = {1, . . . , Dmax } H: the set of surgeons that practice at the hospital (including VMOs). H = {1, . . . , Hmax } S: the set of surgical specialties (where specialty Smax is for VMOs). S = {1, . . . , Smax } T : the set of operating rooms. T = {1, . . . , Tmax } P: the set of patients on the waiting list at the start of the planning horizon. P = {1, . . . , Pmax }. 3.3. Vector parameters
αc : the proportion of patients in urgency category c that must be seen on time, ∀c ∈ C δp : the number of days until patient p is due at the start of the period, ∀p ∈ P Wp : the number of days that patient p will have been on the waiting list at the end of the period, ∀p ∈ P τb : the length of time block b in hours, ∀b ∈ B Ns− : the minimum number of OR hours assigned to specialty s over the time horizon, ∀s ∈ S Ns+ : the maximum number of OR hours assigned to specialty s over the time horizon, ∀s ∈ S
Lp : the of patient p’s surgery, ∀p duration ln N µp , σp2
βpc =
∈
1 if patient p is in urgency category c 0 otherwise,
∀p ∈ P , c ∈ C 1 if patient p is in the treat in turn cohort λp = 0
otherwise,
∀p ∈ P 1 if patient p can be treated by surgeon h Eph = 0
otherwise,
∀p ∈ P , h ∈ H 1
Fhdb =
if surgeon h is available on day d during block b 0 otherwise,
∀h ∈ H , d ∈ D, b ∈ {1, 2} 1 if surgeon h is a member of specialty s Ghs = 0
otherwise,
∀h ∈ H , s ∈ S 1
if operating room t is equipped for surgeries by specialty s 0 otherwise,
Rts =
∀t ∈ T , s ∈ S . 3.4. Decision variables 1
if room t is used by specialty s during time block b on day d 0 otherwise,
Xtsdb =
∀t ∈ T , s ∈ S , d ∈ D , b ∈ B 1
Ytpdb =
if room t is used for patient p during time block b on day d 0 otherwise,
∀t ∈ T , p ∈ P , d ∈ D, b ∈ B
P. Lp
∼
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
1
if room t is used by surgeon h during time block b on day d 0 otherwise,
Zthdb =
A surgeon can only be assigned an OR if their specialty is also assigned to that OR during the time block.
∀t ∈ T , h ∈ H , d ∈ D, b ∈ B
Zthdb ≤
if patient p has on time surgery, or is on time on the final day 0 otherwise,
Xtsdb ≤ Rts ,
otherwise,
∀p ∈ P . Utdb : the expectation of the sum of surgical durations assigned to room t on day d in block b Vtdb : the variance of the sum of surgical durations assigned to room t on day d in block b.
∀t ∈ T , s ∈ S , d ∈ D, b ∈ B.
(8)
Ytpdb eµp
eσp , 2
∀t ∈ T , d ∈ D, b ∈ B.
(9)
Constraint (10) is used to calculate the variance of the sum of surgical durations assigned to each time block. Vtdb =
Ytpdb eσp − 1 e2µp +σp , 2
2
p∈P
3.6. Constraints
∀t ∈ T , d ∈ D, b ∈ B.
It is not possible to share a block between different specialties. There must be at most one specialty assigned to each block.
∀t ∈ T , d ∈ D, b ∈ B.
(1)
Each OR can be assigned in both morning and the afternoon or for a full day only. Xtsd1 + Xtsd2 + 2Xtsd3 ≤ 2,
∀t ∈ T , d ∈ D.
(2)
s∈S
(10)
Constraint (11) ensures a 95% confidence level that surgeries assigned to each time block will be completed on time, given the above approximations.
s∈S
(7)
p∈P
The objective of the static model is to maximise the number of surgeries performed planning horizon. during the Maximise t ∈T p∈P d∈D b∈B Ytpdb .
Xtsdb ≤ 1,
∀t ∈ T , h ∈ H , d ∈ D, b ∈ B.
The duration of each surgery is stochastic and can vary due to a number of different factors. Based on the work by Strum et al. [28], Stuart and Kozan [29,31] approximate the sum of lognormally distributed surgical durations using a lognormal distribution. Constraint (9) is used to calculate the expectation of the sum of surgical durations assigned to each block. Utdb =
3.5. Objective function
Xtsdb Ghs ,
An OR can only be used by a specialty if that OR is properly equipped for use by that specialty.
∀p ∈ P 1 patient p is a member of the longest waiting cohort θp = 0
s∈S
1
γp =
55
2 Utdb
2 Vtdb + Utdb
Vtdb +U 2
1.64485×ln
e
U2 tdb
tdb
≤ τb ,
∀t ∈ T , d ∈ D, b ∈ B.
(11)
A patient’s surgery can be performed at most once during the planning horizon.
Each specialty must not be allocated more than its maximum number of OR hours over the planning horizon.
Ytpdb ≤ 1,
∀p ∈ P .
(3)
t ∈T d∈D b∈B
∀s ∈ S .
(12)
t ∈T d∈D b∈B
If a surgeon is a member of a specialty at the hospital then they cannot share a time block with any other surgeons. It is possible for VMOs to share a time block amongst themselves. This reflects current hospital practices.
τb X tsdb ≤ Ns+ ,
Zthdb ≤ 1 + XtSmax db
h∈H
τb X tsdb ≥ Ns− ,
∀s ∈ S .
(13)
t ∈T d∈D b∈B
GhSmax − 1 ,
h∈H
∀t ∈ T , d ∈ D, b ∈ B.
Each specialty must be allocated at least its minimum number of OR hours over the planning horizon.
(4)
A surgeon can only be assigned a time block if he/she is available during that time block. A surgeon cannot be assigned to a morning and/or afternoon block and a full day block on the same day. A surgeon can only be assigned to a full day block if they are available in both the morning and the afternoon.
Constraint (14) determines whether a patient’s surgery was performed on time, or if the patient is not overdue at the end of the planning horizon.
δp − (Dmax − 1) 1 −
Ytpdb
t ∈T d∈D b∈B
−
(d − 1)Ytpdb ≥ M γp − 1 ,
∀p ∈ P .
(14)
t ∈T d∈D b∈B
Zthdb + Zthd3 ≤ Fhdb ,
∀h ∈ H , d ∈ D, b ∈ {1, 2}.
(5)
t ∈T
A patient can only be treated by a surgeon if the surgeon is capable of performing that particular surgery. This constraint is especially necessary in teaching hospitals. Ytpdb ≤
h∈H
Zthdb Eph ,
∀t ∈ T , p ∈ P , d ∈ D, b ∈ B.
A certain proportion of patients in each category must be seen on time. This proportion is in line with Australian government standards.
p∈P
(6)
γp βpc ≥ αc
βpc ,
∀c ∈ C .
(15)
p∈P
A patient can only be in the longest waiting cohort if they are eligible for the treat in turn policy. This means a patient must be
56
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
ready for care at the start of the month, a public patient, in urgency category two or three and assigned to a surgical specialty.
θp ≤ λp ,
∀p ∈ P .
(16)
The size of the longest waiting cohort is equal to the number of treated in turn eligible patients treated as elective patients during the planning horizon.
θp =
p∈P
Ytpdb λp .
(17)
t ∈T p∈P d∈D b∈B
The longest waiting cohort must be composed of the patients that have been on the waiting list for the longest amount of time. max p∈p
1 − θp Wp λp ≤ θp Wp + 1 − θp M .
(18)
The proportion of longest waiting patients that were treated in turn must meet government standards.
θp Ytpdb ≥ π
t ∈T p∈P d∈D b∈B
θp .
(19)
p∈P
4. Solution approaches As the MSSP and SCAP are both NP-hard, the computational time required increases exponentially with an increase in patients, surgeons, rooms, etc. It is necessary to use heuristics to obtain good, practical solutions in a reasonable amount of time. The use of heuristics often allows for simple modifications to adapt to changes in problem specifications. Construction heuristics are used to find an initial, feasible solution. Hans et al. [32] present a robust formulation of the SCAP in which the sum of surgical durations is assumed to be normally distributed. The authors implement a number of construction heuristics for the SCAP including a first fit approach, longest processing time, and random sampling. Once an initial solution is found, improvement heuristics are used to improve initial solutions. Hans et al. [32] use SA with random swaps based on either inserting a random patient in a suitable random day, or swapping two suitable surgeries. The authors found that SA performed well, but required large amounts of computational time. Agnetis et al. [19] implement a decomposition approach to the combined MSSP SCAP in which they formulate and solve the SCAP as a number of knapsack problems. The MSSP formulation considered by the authors is somewhat simplified as the authors do not consider stochastic surgical durations. This simplifies the implementation of the knapsack formulation. Riise and Burke [33] present metaheuristics to concurrently solve the SCAP and SCSP. The authors use steepest variable neighbourhood search (SVNS) as a comparative baseline for their iterative local search method which utilises variable neighbourhood descent. Landa et al. [34] also consider the SCAP and SCSP. The authors implement a hybrid optimisation algorithm, combining both Monte Carlo simulation and neighbourhood search. The algorithm used consists of two phases, the first to solve the SCAP and the second to address the SCSP. The authors also consider stochastic surgical durations. Other metaheuristics frequently implemented to solve OT planning problems include ant colony optimisation (ACO) [35], tabu search (TS) [20] and genetic algorithms (GA) [18]. The choice of heuristic is largely dependent on the formulation of the problem and the size of the problem instances. As the public hospital used in the case study has a large surgical department with a long waiting list it is inappropriate to evaluate all neighbouring solutions or to use population based metaheuristics (ACO, GA, etc.) as this would require excessive
amounts of memory on a real-life instance. For this reason, we must choose metaheuristics that can be implemented without large memory requirements. Not only is SA simple to implement, it has been used on a variety of hospital planning problems (e.g. [32]) and frequently outperforms other metaheuristics including ACO, TS and GA (cf. Birattari [36]). SA has also presented some of the best known solutions to a number of combinatorial optimisation problems [37]. Due to its very good performance on a wide variety of problems and its appropriateness with regards to the problem structure, SA (in conjunction with an improved version) is suitable for use as a benchmark metaheuristic on this problem. RVNS is well suited to large nonlinear optimisation problems. RVNS has been used on a number of hospital planning problems (e.g. [38]). Although the optimality gap increases with increasing problem size, RVNS can be used to find solutions to very large MINLP problems whereas commercial solvers often fail at medium instances. Both SA and RVNS are suited to the problem structure and size (in terms of memory requirements), easy to implement, and have been shown to outperform other common metaheuristics on a wide variety of optimisation problems. As such, we use these metaheuristics to compare performance with the hybridised metaheuristics presented. Future work includes the incorporation of other suitable metaheuristics including TS. We also present a modified SA in which different neighbourhood structures are used with certain probabilities over time. SA and RVNS are hybridised to utilise the best properties of each of the metaheuristics. To the best of our knowledge hybridised SA–RVNS has not been implemented on OT planning/scheduling problems or MINLP of this magnitude. As such, the use of hybridised SA–RVNS on this large MINLP OT planning and scheduling problem is quite innovative. In this section we present a number of metaheuristics designed to produce good and feasible solutions to the MINLP formulation given in Section 3. We begin by presenting the solution neighbourhoods that are implemented in conjunction with the metaheuristics. 4.1. Solution neighbourhoods The difficulty in implementing certain metaheuristics on other, more general, MINLP models lies in the definition of solution neighbourhoods specific to those problems. Here, we consider three levels of neighbourhood structure. These levels are specialty, surgeon and patient. In order to maintain feasibility, we must consider the implications when we change one of these levels. If the specialty assigned to a block is changed, we must also change the surgeon and patients assigned to that block to ensure constraints (6) and (7) still hold. Similarly, if we only change the surgeon assigned to the block, we may also be required to change the patients assigned. Switching patients does not affect either of the upstream levels. In assigning new specialties, surgeons and patients, we must ensure that availability and suitability constraints are upheld. We also consider neighbourhood changes that swap a specialty from a full day block to two half day blocks and vice versa. These neighbourhood structures ensure that the solution produced at each of the iterations satisfies constraints (1)–(12), (14) and (16)–(19). Depending on the composition of the waiting list, availability of surgeons and suitability of ORs, it may not be possible to satisfy constraints (13) and (15). As such, we cannot guarantee that the following neighbourhood swaps will produce solutions feasible with respect to constraints (13) and (15). We only ever accept a move if the new solution is feasible with respect to all
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
constraints presented in Section 3. Thus, if it is possible to satisfy all constraints, the solutions found using the metaheuristics will be feasible. When finding a random neighbouring solution we perform one of five neighbourhood swaps described below. 1. Swap specialty, surgeon and patients: – Randomly select an OR, day and block. – Remove the patients, surgeon(s) and specialty (if any) currently assigned to that OR during the selected time block. – Randomly select a suitable specialty and surgeon(s) from the new specialty. Assign these to the selected time block. – Randomly assign appropriate patients to the time block without violating the constraints presented in Section 3. 2. Swap surgeon and patients: – Randomly select an OR, day and block. – Remove the patients and surgeon(s) (if any) currently assigned to that OR during the selected time block. – Randomly select a suitable surgeon. Assign the new surgeon to the selected time block. – Randomly assign appropriate patients to the time block without violating the constraints presented in Section 3. 3. Swap patients: – Randomly select an OR, day and block. – Remove the patients (if any) currently assigned to that OR during the selected time block. – Randomly assign appropriate patients to the time block without violating the constraints presented in Section 3. 4. Select a full day block and swap to two half day blocks: – Randomly select an OR and day that is currently assigned to a single specialty and surgeon for the entire day. – Remove the patients, surgeon and specialty currently assigned to that OR during the time block. Reassign the same surgeon and specialty to both the morning and afternoon blocks on the same day in the same OR. For both the morning and afternoon time block, randomly assign appropriate patients to the time blocks without violating the constraints presented in Section 3. 5. Select two half day blocks of the same specialty and swap to a full day block: – Randomly select an OR and day that is currently assigned to the same surgeon across both the morning and afternoon time blocks. – Remove the patients, surgeon and specialty assigned to the half day blocks. Reassign the same surgeon and specialty to a full day block in the same OR on the same day. – Randomly assign appropriate patients to the new full day block without violating the constraints presented in Section 3. In addition to the above neighbourhood swaps, after every N iterations the solution is ‘cleaned’ by grouping scheduled patients that can be seen by the same surgeon and reassigning these patients where appropriate in order to reduce the number of time blocks used. No patients are added or removed from the schedule; the assignment of patients is simply rearranged. As the current subset of patients is feasible with respect to all constraints, the ‘cleaned’ solution will also be feasible. This assists in improving the quality of the solution obtained by each of the metaheuristics. Surgeons and specialties to be placed in the schedule (swap type one or two) are not selected according to a uniform distribution. Suitable specialties are selected according to a probability proportional to the maximum number of patients of that specialty that could fit within the time block. Suitable surgeons are selected according to a probability proportional to the number of available patients that can be scheduled at that iteration.
57
4.2. Simulated annealing SA is a global optimisation metaheuristic that was introduced by Kirkpatrick et al. [39]. Based on the idea of temperature control during the cooling of metals, SA allows a worse solution to be accepted with a probability that decreases over time. The acceptance of worse solutions allows SA to escape local optimum in the search for the global optimum. SA is often used on hospital planning problems due to its ease of implementation and solution quality. Ceschia and Schaerf [40] address the patient admission planning problem with the inclusion of OR utilisation constraints. The authors use a metaheuristic based on SA with particular focus on the use of appropriate structures. On a similar problem Ceschia and Schaerf [37] used this method produced best known solutions to a number of sample instances. Fügener et al. [41] implement a number of metaheuristics including SA to produce a MSS for a small Dutch hospital. The authors consider the impact of the MSS on the occupancy of a number of downstream wards. The authors found that SA performed well in terms of both solution quality and time required. Birattari [36] discusses a number of parameter tuning algorithms including a brute force approach and the use of the F-Race algorithm. The author also investigates the use of a variety of metaheuristics on several instances of well-known combinatorial optimisation problems. After careful parameter tuning, SA tends to outperform other metaheuristics including ACO, GA and TS. Due to its ease of implementation, suitability for large problem instances and quality of solutions, as seen in the work referenced above, we use SA as a baseline metaheuristics to compare with the performance of the other metaheuristics on the model presented R in Section 3. We implement SA in MATLAB⃝ through the use of the solution neighbourhoods presented in Section 4.1. At each iteration we randomly select a neighbourhood structure and find a neighbouring solution. We accept solutions according to the acceptance probability; the best solution is only updated if the solution is both feasible and improved. The acceptance probability is given by the expression exp((new obj − old obj)/T ), where T is the metaheuristic temperature at that iteration. 4.3. Adaptive simulated annealing In order to produce improved solutions, we modify our SA based metaheuristic to a more adaptive form. The acceptance probability of a new solution is dependent on both the change from the previous solution and on the time at which the solution is found. This ensures a smoother transition between acceptance probabilities as time progresses. The acceptance probability is given by the expression exp((new obj − old obj) × toc kp /T ), where T is the metaheuristic temperature, toc is the time since starting the heuristic and kp is the exponent parameter. Whilst the selection of neighbourhood structure in SA is done using a uniform distribution, we increase the probability that swaps of type four and five occur over time as these swaps are more useful in a mostly full schedule. Initially, swaps of type four and five cannot be selected. After just over 75% of the iterations have elapsed the swap type is selected according to a uniform distribution. The probability of selecting swaps one, two, or three is given by Eq. (20) and the probability of selecting swaps four or five is given by Eq. (21).
−1
(20)
.
(21)
p1,2,3 (t ) = 3 + 2 (1 + exp (−t + 0.75 × Imax ))−1
p4,5 (t ) = (3 (1 + exp (−t + 0.75 × Imax )) + 2)
−1
4.4. Reduced variable neighbourhood search Variable neighbourhood search (VNS) is a metaheuristic proposed by Mladenović and Hansen [42]. Originally used to solve
58
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
combinatorial optimisation problems, VNS has since been successfully extended to MILPs and some MINLPs. VNS utilises neighbourhood structures, local search, and perturbation to produce good feasible solutions to optimisation problems. A number of variations of VNS exist which are more suited to certain problem structures. A detailed introduction to VNS can be found in [43]. A number of different hospital planning problems including the nurse rostering problem [44] and the scheduling of patients within a hospital [38] have been solved using VNS inspired metaheuristics. Recently, Aringhieri et al. [45] solve the SCAP under a block scheduling policy with the aim of levelling bed occupancy in downstream wards. The authors implement VNS, RVNS and an adaptive version of the algorithm. Depending on problem formulation and metaheuristic parameters used, VNS can produce solutions as good as exact MINLP solvers on small to medium instances (e.g. [46]). Whilst the optimality gap increases with the size of the problem, VNS can provide good feasible solutions to realistic instances in short amounts of time. On the other hand, commercial MINLP solvers often fail on realistic instances. Due to the size of our problem, we consider RVNS in which a random neighbour is selected and chosen if the objective is improved. This vastly reduces the computational time and storage required and makes the work accessible to those without high performance computers. 4.5. Hybridised simulated annealing and reduced variable neighbourhood search There are limitations associated with both SA and RVNS that can reduce the average quality of solutions found. By randomly selecting solution neighbourhoods, SA is not able to make use of the neighbourhood structures present in the combined MSSP SCAP. It is a common occurrence for RVNS to become ‘stuck’ in a local optima due to its inability to accept worsening moves. To avoid each of these issues, we combine the beneficial features of both adaptive SA (cf. Section 4.3) and RVNS (cf. Section 4.4). In order to make use of the neighbourhood structures present in the problem, we continue to use the framework provided by the RVNS metaheuristic. In particular, we do not randomly select a neighbourhood type, but begin at the first neighbourhood structure and move down the list of structures if the neighbouring solution is not accepted. If the neighbouring solution is accepted then we move back up the list of structures. Unlike RVNS, hybridised SA–RVNS occasionally accepts nonimproving moves in an attempt to escape local maxima. As in the adaptive SA metaheuristic, the acceptance probability is both change and time dependent. Throughout the literature there are a few instances of hybridised SA–RVNS being implemented to solve a variety of optimisation problems. To the best of the authors’ knowledge, hybridised SA–RVNS has not been implemented on a MINLP model of this magnitude. The hybridised SA–RVNS algorithm is as follows: SA–RVNS AlgorithmStep 1 Initialise the algorithm: 1.1 Start with an initial solution. 1.2 Check feasibility. 1.3 If the initial solution is feasible, update the best solution to be the initial solution. Otherwise, let the best objective be zero. 1.4 k ← 1. Step 2 While the iteration limit has not been exceeded, iterate over the following steps:
2.1 Find a random neighbour, from neighbourhood k, of the previous solution. 2.2 Calculate the number of patients scheduled in the neighbouring solution. 2.3 Check the feasibility of the neighbouring solution. 2.4 If the acceptance probability is greater than some uniform random number and the neighbouring solution is feasible then update the latest solution to the neighbour found in step 2.1. 2.4.1 If the latest solution is also the best so far, update the best solution. 2.5 If the acceptance probability was too small, or the solution was not feasible, k ← k + 1 2.6 If k > Kmax, k ← 1. 4.6. Upper bounds In this section we present the manner in which upper bounds for the maximisation problem are found. We apply some simplifications to the model that allow us to solve a relaxed version using IBM ILOG CPLEX Optimization Studio 12.6. Due to the availability of data we allocate the same lognormal distribution parameters to the duration of all types of surgery within a specialty. As we only consider nine specialties in this section, it is possible to calculate the maximum number of patients scheduled in each type of block, whilst satisfying constraints (9)– (11). We consider both half and full day blocks which are five and ten hours respectively. This data is shown in Table 4. By utilising the data in Table 4 we are able to remove (highly nonlinear) constraints (9)–(11) from the model. Constraint (22) can be used as a simple replacement. Note that this substitution of constraints is only simple due to the assumption that surgical durations within the same specialty are independent and identically distributed. Asb : the maximum number of patients of specialty s that can be seen in block type b with less than 5% chance of exceeding the length of the time block. ∀s ∈ S , b ∈ B
p∈P
Ytpdb ≤
Asb Xtsdb ,
∀t ∈ T , d ∈ D, b ∈ B.
(22)
s∈S
The constraints used to address the treat in turn policy (constraints (16) through (19)) result in a quadratic constraint that does not satisfy the positive semi-definite property required by CPLEX. It is possible to introduce a set of constraints to linearise (19); however this results in an additional Pmax binary decision variables. Instead, we remove constraints (16) through (19) as it is too computationally expensive to include these when solving realistic instances. Due to the complexity of the problem it is also necessary to relax the binary decision variable Ytpdb such that Ytpdb ∈ [0, 1]. In summary, to calculate the upper bound on the randomly generated instance considered in Section 5.2 we allow the variable Ytpdb to be a real number on the interval [0, 1] and we include only constraints (1)–(8), (12)–(15), and (22). This may result in noninteger solutions that do not satisfy the treat in turn policy. After implementing these strategies we are able to solve each instance in CPLEX to produce an upper bound on the optimal solution to the model presented in Section 3. 5. Computational results In the following section we present computational results including metaheuristic parameter tuning, the performance of each of the metaheuristics is compared and the performance of the model presented in Section 3 is compared to recent historical data. All computational experiments in Sections 5.1 and 5.2 were R R performed using MATLAB⃝ on an Intel⃝ CoreTM i7 processor at 3.40 GHz with 16.0 GB of RAM. Computational experiments in Section 5.3 were performed using the university’s high performance computing facility.
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 5 Problem parameters. Parameter
Description
Value
Pmax Hmax Smax Tmax Dmax Bmax Cmax
Number of patients Number of surgeons Number of specialties Number of ORs Number of days Number of blocks Number of urgency categories
2094 107 10 20 5 3 3
5.1. Parameter tuning In this subsection we analyse the performance of SA, adaptive SA and hybridised SA–RVNS under different metaheuristic parameters. The performance of the above metaheuristics are analysed on a partially random dataset. The number of surgeons, specialties, ORs, days and blocks are kept constant. The surgeons and patients were taken from the elective surgery waiting list. The surgeon availabilities used were randomly generated. The same availability schedule was used for each run. The main scalar problem parameters used in the model can be found in Table 5. These values result in a problem with 667,488 binary decision variables, 600 non-negative continuous decision variables and 672,278 constraints. In this section we investigate all suitable combinations of metaheuristic parameters. We further investigate metaheuristic parameters based on promising behaviour. In each case we run the metaheuristics until variance of the best objective function value changes by less than 0.1 units after a further ten runs, or until 300 runs have been performed. 5.1.1. Simulated annealing A number of parameters are required by the SA metaheuristic. In particular, these are the initial temperature (T ), the number of iterations between temperature changes (imax ), the factor by which the temperature changes (α), and the minimum (stopping) temperature (Tmin ). In order to understand the full capabilities of SA on the model presented in Section 3 we must determine the most appropriate values for each of the metaheuristic parameters. After preliminary investigation, we consider all combinations of the SA parameter values given in Table 6. The five best performing SA parameter sets are shown in Table 7. The table also shows the average number of patients scheduled, lower and upper bounds of a 95% confidence interval and the average time required for each run. There is significant evidence that the average objective function obtained under T = 10, α = 09, imax = 150, Tmin = 0.001 is greater than that obtained under each of the other SA parameter combinations.
59
5.1.2. Adaptive simulated annealing The metaheuristic parameters required by adaptive SA are the initial temperature (T ), the maximum number of iterations and kp , where the probability of accepting a swap is given by exp 1cost × iter kp /T . Whilst it is also possible to change the probability function associated with selecting each type of neighbourhood swap, the current form performs well and will be used for the remainder of the investigation. To determine the best performance of adaptive SA on the model presented in Section 3 we consider suitable combinations of the adaptive SA parameters presented in Table 8. Table 9 shows the five best performing adaptive SA parameter combinations. The table shows the average objective function obtained, lower and upper bounds on a 95% confidence interval, and the average time required for a single run of the metaheuristic. 5.1.3. Reduced variable neighbourhood search The only metaheuristic parameter required by RVNS is the number of iterations performed. Here we consider between 400 and 6000 iterations. Table 10 shows the five best performing RVNS parameters along with the average number of patients scheduled, the lower and upper bounds of a 95% confidence interval, and the average time required for a single run. After 2200 iterations, an increase in the maximum number of iterations does not result in an increase in average objective function value. This is as RVNS tends to become ‘stuck’ at local optima and the objective function value plateaus. Note that RVNS does not perform better with only 2200 iterations, it merely performs as well as it does with a higher number of iterations. This is reflected in the overlapping confidence intervals. 5.1.4. Hybridised SA–RVNS The parameters required by hybridised SA–RVNS are the same as those required by adaptive SA. After initial analysis we consider the metaheuristic parameters shown in Table 11. Table 12 shows results of the SA–RVNS parameter tuning. The table contains the best performing hybridised SA–RVNS parameters, the average objective function, lower and upper bounds of a 95% confidence interval, and the average time required for a single run. At a 95% confidence level there is insufficient evidence of a difference in mean objective function obtained by each of these top five hybridised SA–RVNS parameter combinations. 5.2. Metaheuristic performance on random data Applying the results of the metaheuristic parameter tuning, any further analysis of the metaheuristics will be performed under the metaheuristic parameters shown in Table 13.
Table 6 SA parameters. Parameter
Description
Values
T
Initial temperature Temperature decrease factor Iterations between decreases Minimum temperature
1000, 100, 10 0.9, 0.7, 0.5 150, 100, 50 T × 10−1 , T × 10−2 , T × 10−3
α imax Tmin
Table 7 Results of SA parameter variation. T
α
imax
Tmin
Runs
Average
Lower
Upper
Time (s)
10 10 100 10 100
0.9 0.9 0.9 0.9 0.9
150 150 150 100 100
0.001 0.010 0.010 0.001 0.010
40 110 70 190 210
270.13 268.14 267.94 267.85 264.45
269.21 267.60 267.26 267.37 263.84
271.04 268.67 268.63 268.33 265.06
202.05 150.31 151.71 131.86 96.90
60
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
Table 8 Adaptive SA parameters. Parameter
Description
Values
T kp Imax
Initial temperature Exponent in acceptance probability function Maximum total number of iterations
1000, 100, 10 6, 4, 3, 2.5, 2 5000, 6000, 7000, 8000, 9000
Table 9 Results of adaptive SA parameter variation. T
kp
Imax
Runs
Average
Lower
Upper
Time (s)
100 100 100 100 100
2.5 3 2.5 6 3.5
9000 9000 8000 9000 9000
80 130 110 100 70
270.81 270.78 270.61 270.43 270.40
270.33 270.34 270.04 269.97 269.81
271.30 271.23 271.18 270.89 270.99
155.86 156.87 137.96 156.99 157.13
This reduces the variance of historical surgical durations and allows for more accurate estimation of surgical durations. The MINLP formulation provided in Section 3 does not allow a surgery to be scheduled if the 95th percentile of that surgery’s duration is greater than τBmax . Here, we replace constraint (11) with constraint (23) and include the constraints (24) and (25) that allows overtime when scheduling a single surgery in a full day block. if there is more than one patient scheduled in room t , on day d, block Bmax 0 otherwise,
1
Table 10 Results of RVNS parameter variation.
φtd =
Imax
Runs
Average
Lower
Upper
Time (s)
2200 2400 5000 3200 2300
60 50 40 270 130
265.58 265.18 264.85 264.80 264.70
264.51 263.91 263.40 264.18 263.78
266.65 266.45 266.30 265.41 265.62
174.12 186.07 388.37 247.76 178.17
∀t ∈ T , d ∈ D . Mt : an upper bound on the time that a surgery can take (in hours).
2 Utdb
Summary statistics are shown in Table 14, including the lower and upper bounds of a 95% confidence interval. From these statistics it can be seen that whilst at a 95% confidence level there is insufficient evidence that the mean objective found by Adaptive SA is different to that found by SA, there is very significant evidence that hybridised SA–RVNS outperforms each of the other metaheuristics. It is also worth noting the variance of the results. The objective value found by RVNS has a large variance as RVNS displays inconsistent performance. Much less variation is seen in the objective found using hybridised SA–RVNS which shows the metaheuristic to consistently perform well. Whilst there is insufficient evidence of a difference between the mean objective function value found by SA and Adaptive SA, Adaptive SA requires less time and has a lower variance. One limitation of the model presented in Section 3 is that no ENT surgeries can be scheduled. This is as, under the data shown in Table 4, it is not possible to schedule any ENT surgeries as the 95th percentile of the duration of a single surgery is greater than ten hours. Modifications required to address this issue are presented in Section 5.3.
2 Vtdb + Utdb
× exp 1.64485 ×
ln
2 Vtdb + Utdb
2 Utdb
≤ τb ,
∀t ∈ T , d ∈ D, b ∈ {1, 2} YtpdBmax , ∀t ∈ T , d ∈ D Mp φtd + 1 ≥
(23) (24)
p∈P
2 UtdB max
2 VtdBmax + UtdB max
2 VtdBmax + UtdB max × exp 1.64485 × ln 2
UtdBmax
≤ τBmax + (1 − φtd ) Mt ,
∀t ∈ T , d ∈ D.
(25)
Since any surgeons in ENT are now considered members of ENT1, ENT2, and ENT3 we need to include additional constraints to ensure patients from each subspecialty are scheduled appropriately. We include a binary decision variable χtdb and constraints (26) and (27). Note that ENTc is used to denote ENT subspecialty c.
5.3. Model performance
if room t on day d, block b is scheduled for an ENT subspecialty 0 otherwise, 1
In this section we investigate the performance of our model on the case study presented in Section 2. We use real data from June to September 2015. As we do not have staff availabilities, we use historical data to approximate actual availabilities. When conducting the metaheuristic parameter tuning no surgeries were scheduled for ENT as the 95th percentile on surgical duration was greater than ten hours. In this section, we split ENT into three subgroups separated by a patient’s urgency category.
χtdb =
∀t ∈ T , d ∈ D , b ∈ B χtdb ≥ Xtsdb , ∀t ∈ T , d ∈ D, b ∈ B
(26)
s∈ENT
Ytpdb βpc ≤ Xt (ENTc )db − χtdb + 1,
∀t ∈ T , p ∈ P , d ∈ D, b ∈ B, c ∈ C .
Table 11 Hybridised SA–RVNS parameters. Parameter
Description
Values
T kp Imax
Initial temperature Exponent in acceptance probability function The maximum total number of iterations
1000, 100, 10 1.5, 2, 2.5, 3 4000, 5000, 6000, 7000, 8000, 9000
(27)
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 12 Results of hybridised SA–RVNS parameter variation. T
kp
Imax
Runs
Average
Lower
Upper
Time (s)
1000 1000 10 1000 10
2.5 3.0 2.5 3.0 2.5
9000 7000 9000 8000 8000
90 60 90 130 80
272.19 272.08 271.92 271.88 271.85
271.76 271.62 271.54 271.52 271.40
272.61 272.55 272.30 272.25 272.30
153.68 119.67 158.54 137.07 140.39
Hence, the MINLP formulation used in this section consists of constraints (1)–(10), (12)–(19), and (23)–(27). Very little modification is required to calculate the new upper bound. Note that the metaheuristic solutions reported are feasible with respect to constraints (1)–(10), (12)–(19) and (23)–(27), i.e. all relevant constraints. The upper bound is calculated using constraints (1)–(8), (12)–(15), (22), (26) and (27). The upper bound is not guaranteed to be feasible with respect to the treat in turn policy. The upper bound also relies on a continuous relaxation of Ytpdb . Queensland University of Technology’s high performance computing facility was used to run CPLEX. Each instance was run across sixteen cores, with 128 GB of memory and a wall time of 168 h. Note that the time reported in Table 16 is the CPU time used. Table 15 shows the number of patients on the elective waiting list over the seven weeks under study. The table is split by urgency category and shows the number of patients on time and overdue at the time the waiting list was taken. The number of surgeries performed in each of these weeks is shown in Table 16. This includes both elective and emergency surgeries. The ratio of total number of surgeries waiting in each category is approximately 1:4:5 in categories one, two, and three respectively. It should be noted that the number of patients to be scheduled each week is now roughly 2900. In this section we consider 26 specialties and approximately 230 surgeons. As such, the problem size is significantly larger than that Section 5.1. The computational experiments performed in this section use the metaheuristic parameters given in Table 13. As in Section 5.1 we run each metaheuristic either until the variance of the best objective function value changes by less than 0.1 units after a further ten runs. We impose a maximum of 300 runs of each metaheuristic. Table 16 shows the results of computational experiments using the seven weeks of data available. The table shows the mean, variance, worst and best objective function value obtained by each metaheuristic after at most 300 runs. These results confirm the trends seen in Section 5.2. Two-sample t-tests were performed
61
on each combination of metaheuristics. For each week, the results display the same pattern. RVNS is the worst performing metaheuristic in terms of both lowest mean objective function value and highest objective function value variance. Whilst the mean objective function value obtained by SA is larger than that obtained by RVNS in six of the seven weeks (p ≪ 0.001), SA produces significantly worse solutions than Adaptive SA (p ≪ 0.001). When performing metaheuristic parameter tuning a maximum of 9000 iterations were considered for both Adaptive SA and hybridised SA–RVNS as the metaheuristics performed as well as or better than SA in considerably less time. In future it may be worth considering a larger total number of iterations. In this case we consider computational time as a limited resource. The best performing metaheuristic for each week was hybridised SA–RVNS. This heuristic provides solutions with a low variance and highest mean objective function value. Hybridised SA–RVNS required the least (or second least) amount of computational time and produced significantly better results than the next best metaheuristic, Adaptive SA (p ≪ 0.001). In several instances, the best objective function value obtained using metaheuristics was within ten units of the best integer solution obtained by CPLEX on the relaxed problem. The range of objective function values obtained through the hybridised SA–RVNS metaheuristic is very promising, especially considering the size of the real life instances addressed. When comparing the results produced by the metaheuristics using the model presented in Section 3 to historical data we find that the metaheuristics are able to schedule around 50–100 additional surgeries each week. This result should be considered with caution as surgeon availabilities were merely estimates based on historical data (due to the lack of available data); however these estimates are as accurate as possible. There is also the possibility that downstream wards may be unable to satisfy the increased demand. 6. Conclusions and future work In this paper we addressed the combined MSSP SCAP in Australian public hospitals using a MINLP formulation and tackled all three levels of OT planning/scheduling problems simultaneously. Whilst commercial solvers are unable to cope with the size of the case study considered the metaheuristics and hybridised metaheuristics produce good feasible solutions in reasonable amounts of computational time. In particular
Table 13 Metaheuristic parameters. Parameter Description
SA
RVNS Adaptive SA
SA–RVNS
T
10 0.9 150
NA NA NA
1000 NA NA
α imax kp Tmin Imax
Starting temperature Temperature decrease Number of iterations between temperature decrease Acceptance probability function exponent The stopping temperature The maximum total number of iterations
100 NA NA
NA NA 2.5 0.001 NA NA NA 2200 9000
2.5 NA 9000
Table 14 Summary statistics of computational experiments. Method
Runs
Average
Variance
Lower
Upper
Time (s)
RVNS SA Adaptive SA SA–RVNS Upper bound
60 40 80 90 –
265.58 270.13 270.81 272.19 305.00
17.91 8.63 4.64 4.24 –
264.51 269.21 270.33 271.76 –
266.65 271.04 271.30 272.61 –
174.12 202.05 155.86 153.68 50.01
62
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 15 Case study weekly waiting list data by urgency category. Week
1 2 3 4 5 6 7
1
2
3
Total
On time
Overdue
On time
Overdue
On time
Overdue
337 323 313 292 302 280 294
0 0 0 0 0 22 1
1127 1118 1117 1131 1102 1077 1130
74 70 63 69 75 104 76
1326 1378 1382 1372 1367 1364 1364
16 11 13 15 14 13 16
2880 2900 2888 2879 2860 2860 2881
Table 16 Case study weekly results—metaheuristic performance. Week
Method
Mean
Variance
Worst
Best
Runs
Mean time (s)
1
Actual RVNS SA Adaptive SA–RVNS Upper bound
297 308.78 320.64 328.06 330.77 383.01
– 394.34 14.50 26.63 20.47 –
– 148 309 315 316 –
– 328 331 343 339 –
– 300 160 150 110 –
– 198.74 226.67 182.16 183.27 685 938.00
2
Actual RVNS SA Adaptive SA–RVNS Upper bound
303 329.84 332.34 344.76 349.53 391.68
– 68.46 28.72 30.42 22.79 –
– 301 321 327 335 –
– 349 349 357 359 –
– 300 90 170 270 –
– 222.44 270.48 200.24 192.51 4 582 191.00
3
Actual RVNS SA Adaptive SA–RVNS Upper bound
266 333.70 335.05 346.31 351.84 391.53
– 48.44 22.62 20.40 15.95 –
– 308 323 331 333 –
– 353 347 357 362 –
– 220 190 130 170 –
– 219.48 238.59 194.85 190.24 5 008 999.00
4
Actual RVNS SA Adaptive SA–RVNS Upper bound
306 330.09 334.86 343.13 348.57 393.40
– 50.47 27.35 19.54 20.17 –
– 305 318 333 336 –
– 345 348 354 359 –
– 300 110 200 220 –
– 217.42 243.05 193.02 192.12 4 386 662.00
292 325.81 328.26 340.52 345.36 395.31
– 33.89 23.97 32.56 30.40 –
– 312 314 322 326 –
– 341 340 356 357 –
–
5
Actual RVNS SA Adaptive SA–RVNS Upper bound
70 120 300 190 –
– 217.20 235.06 244.90 224.78 8 772 716.00
6
Actual RVNS SA Adaptive SA–RVNS Upper bound
232 312.22 314.91 322.43 325.58 392.54
– 150.22 11.05 17.41 19.09 –
– 147 308 314 317 –
– 331 326 331 336 –
– 300 230 70 60 –
– 199.65 227.15 217.10 173.45 3 364 606.00
7
Actual RVNS SA Adaptive SA–RVNS Upper bound
256 304.01 309.56 316.89 319.36 387.82
– 227.02 15.89 25.14 28.35 –
– 137 301 305 304 –
– 321 320 330 331 –
– 300 90 140 160 –
– 215.10 236.76 202.39 210.66 6 935 316.00
hybridised SA–RVNS, which has not previously been implemented on large MINLP problems or OT planning/scheduling problems, outperformed our chosen benchmark metaheuristics and has potential to be applied to other complex hospital scheduling problems on a large scale. When modelling the combined MSSP SCAP we considered the Australian government’s current waiting list management policies: the treat in turn policy and patient wait targets. Future work includes the consideration of other possible policies that may better manage the elective surgery waiting lists. Possible extensions to the model include incorporating length of stay, emergency arrivals, and capacity of downstream wards.
Acknowledgements This research was funded by the Australian Research Council (ARC) Linkage Grant LP 140100394 and supported by the Princess Alexandra Hospital, Brisbane, Australia. We would like to thank Dr. Andy Wong for his assistance in obtaining data for our case study. We would also like to thank Dr. Michael Sinnot, Dr. David Cook, Sean Birgan and other staff at the PAH for their considerable feedback and their time. Appendix (See Fig. 1.)
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
63
Fig. 1. A typical MSS.
References [1] C. Van Riet, E. Demeulemeester, Trade-offs in operating room planning for electives and emergencies: A review, Oper. Res. Health Care (2015). [2] Y.B. Ferrand, M.J. Magazine, U.S. Rao, Managing operating room efficiency and responsiveness for emergency and elective surgeries — a literature survey, IIE Trans. Healthc. Syst. Eng. 4 (1) (2014) 49–64. [3] F. Guerriero, R. Guido, Operational research in the management of the operating theatre: a survey, Health Care Manage. Sci. 14 (2011) 89–114. [4] I. Adan, J. Bekkers, N. Dellaert, J. Vissers, X. Yu, Patient mix optimisation and stochastic resource requirements: a case study in cardiothoracic surgery planning, Health Care Manage. Sci. 12 (2) (2009) 129–141. [5] J.M. Van Oostrum, M.V. Van Houdenhoven, J.L. Hurink, E.W. Hans, G. Wullink, G. Kazemier, A master surgical scheduling approach for cyclic scheduling in operating room departments, OR Spectrum 30 (2) (2008) 355–374. [6] H. Fei, N. Meskens, C. Chu, A planning and scheduling problem for an operating theatre using an open scheduling strategy, Comput. Ind. Eng. (2010) 221–230. [7] P. Santibáñez, M. Begen, D. Atkins, Surgical block scheduling in a system of hospitals: an application to resource and wait list management in a British Columbia health authority, Health Care Manage. Sci. 10 (3) (2007) 269–282. [8] J.M. Van Oostrum, E. Bredenhoff, E.W. Hans, Suitability and managerial implications of a master surgical scheduling approach, Ann. Oper. Res. 178 (1) (2010) 91–104. [9] A. Agnetis, A. Coppi, M. Corsini, G. Dellino, C. Meloni, M. Pranzo, Long term evaluation of operating theater planning policies, Oper. Res. Health Care 1 (2012) 95–104. [10] J. Beliën, E. Demeulemeester, B. Cardoen, A decision support system for cyclic master surgery scheduling with multiple objectives, J. Sched. 12 (2) (2009) 147–161. [11] S. Chaabane, N. Meskens, A. Guinet, M. Laurent, Comparison of two methods of operating theatre planning: Application in Belgian hospital, J. Syst. Sci. Syst. Eng. 17 (2) (2008) 171–186. [12] A. Testi, E. Tànfani, Tactical and operational decisions for operating room planning: Efficiency and welfare implications, Health Care Manage. Sci. 12 (4) (2009) 363–373. [13] Y. Wang, J. Tang, R.Y.K. Fung, A column-generation-based heuristic algorithm for solving operating theater planning problem under stochastic demand and surgery cancellation risk, Int. J. Prod. Econ. 158 (2014) 28–36. [14] H. Fei, C. Chu, N. Meskens, Solving a tactical operating room planning problem by a column-generation-based heuristic procedure with four criteria, Ann. Oper. Res. 166 (1) (2009) 91–108. [15] H. Fei, C. Chu, N. Meskens, A. Artiba, Solving surgical cases assignment problem by a branch-and-price approach, Int. J. Prod. Econ. 112 (1) (2008) 96–108. [16] H. Saadouli, B. Jerbi, A. Dammak, L. Masmoudi, A. Bouaziz, A stochastic optimization and simulation approach for scheduling operating rooms and recovery beds in an orthopedic surgery department, Comput. Ind. Eng. 80 (2015) 72–79.
[17] K. Steins, F. Persson, M. Holmer, Increasing utilization in hospital operating department using simulation modeling, Simulation 86 (8–9) (2010) 463–480. [18] B. Roland, C. Di-Martinelly, F. Riane, Y. Pochet, Scheduling an operating theatre under human resource constraints, Comput. Ind. Eng. 58 (2) (2010) 212–220. [19] A. Agnetis, A. Coppi, M. Corsini, G. Dellino, C. Meloni, M. Pranzo, A decomposition approach for the combined master surgical schedule and surgical case assignment problems, Health Care Manage. Sci. 17 (2014) 49–59. [20] R. Aringhieri, P. Landa, P. Soriano, E. Tànfani, A. Testi, A two level metaheuristic for the operating room scheduling and assignment problem, Comput. Oper. Res. 54 (2015) 21–34. [21] J.E. Everett, A decision support simulation model for the management of an elective surgery waiting system, Health Care Manage. Sci. 5 (2002) 89–95. [22] D. Min, Y. Yih, An elective surgery scheduling problem considering patient priority, Comput. Oper. Res. 37 (2010) 1091–1099. [23] D. Min, Y. Yih, Managing a patient waiting list with time-dependent priority and adverse events, RAIRO-Oper. Res. 48 (2014) 53–74. [24] I. Marques, M.E. Captivo, M. Vaz Pato, A bicriteria heuristic for an elective surgery scheduling problem, Health Care Manage. Sci. (2014) 1–16. [25] Queensland Health, Treating in turn — calculating the treat in turn percent, in: Department of Health (ed.), 2013. [26] E. Kayış, T. Khaniyev, J. Suermondt, K. Sylvester, A robust estimation model for surgery durations with temporal, operational, and surgery team effects, Health Care Manage. Sci. 18 (2015) 222–233. [27] F. Dexter, J. Ledolter, V. Tiwari, R.H. Epstein, Value of a scheduled duration quantified in terms of equivalent numbers of historical cases, Anesth. Analg. 117 (1) (2013) 205–210. [28] D. Strum, J. May, L. Vargas, Surgical procedure times are well modeled by the lognormal distribution, Anesth. Analg. 86 (1998) 47S. [29] K. Stuart, E. Kozan, Online scheduling in the operating theatre, in: Industrial Engineering and Management Systems Conference, 2009 Kitakyushu, Japan, pp. 802–807. [30] M. Heng, J. Wright, Dedicated operating room for emergency surgery improves access and efficiency, Can. J. Surg. 56 (3) (2013) 167–174. [31] K. Stuart, E. Kozan, Reactive scheduling model for the operating theatre, Flex. Serv. Manuf. J. 24 (4) (2012) 400–421. [32] E. Hans, G. Wullink, M. Van Houdenhoven, G. Kazemier, Robust surgery loading, European J. Oper. Res. 185 (3) (2008) 1038–1050. [33] A. Riise, E. Burke, Local search for the surgery admission planning problem, J. Heuristics 17 (2011) 389–414. [34] P. Landa, R. Aringhieri, P. Soriano, E. Tànfani, A. Testi, A hybrid optimization algorithm for surgeries scheduling, Oper. Res. Health Care 8 (2016) 103–114. [35] C. Rizk, J.-P. Arnaout, ACO for the surgical cases assignment problem, J. Med. Syst. 36 (2012) 1891–1899. [36] M. Birattari, Tuning Metaheuristics: A Machine Learning Perspective, Springer, Berlin, Heidelberg, 2009. [37] S. Ceschia, A. Schaerf, Local search and lower bounds for the patient admission scheduling problem, Comput. Oper. Res. 38 (2011) 1452–1463.
64
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
[38] S. Vlah, Z. Lukač, J. Pacheco, Use of VNS heuristics for scheduling of patients in hospital, J. Oper. Res. Soc. 62 (2011) 1227–1238. [39] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680. [40] S. Ceschia, A. Schaerf, Dynamic patient admission scheduling with operating room constraints, flexible horizons, and patient delays, J. Sched. (2014) 1–13. [41] A. Fügener, E.W. Hans, R. Kolisch, N. Kortbeek, P.T. Vanberkel, Master surgery scheduling with consideration of multiple downstream units, European J. Oper. Res. 239 (2014) 227–236. [42] N. Mladenović, P. Hansen, Variable neighborhood search, Comput. Oper. Res. 24 (1997) 1097–1100.
[43] P. Hansen, N. Mladenović, J. Brimberg, J.M. Pérez, Variable neighborhood search, in: M. Gendreau, J.-Y. Potvin (Eds.), Handbook of Metaheuristics, Springer, US, 2010. [44] E.K. Burke, J. Li, R. Qu, A hybrid model of integer programming and variable neighbourhood search for highly-constrained nurse rostering problems, European J. Oper. Res. 203 (2010) 484–493. [45] R. Aringhieri, P. Landa, E. Tànfani, Assigning surgery cases to operating rooms: A VNS approach for leveling ward beds occupancies, Electron. Notes Discrete Math. 47 (2015) 173–180. [46] A. Alonso-Ayuso, L. Escudero, F.J. Martín-Campo, N. Mladenović, A VNS metaheuristic for solving the aircraft conflict detection and resolution problem by performing turn changes, J. Global Optim. 63 (2015) 583–596.
Contents lists available at ScienceDirect
Operations Research for Health Care journal homepage: www.elsevier.com/locate/orhc
Waiting list management through master surgical schedules: A case study Belinda Spratt, Erhan Kozan ∗ Queensland University of Technology, 2 George St, Brisbane, QLD 4000, Australia
article
info
Article history: Received 24 January 2016 Accepted 11 July 2016 Available online 15 July 2016 Keywords: MSSP SCAP Waiting list management Operating theatre planning Hybrid metaheuristics
abstract In this paper we address the problem of generating master surgical schedules (MSSs) that adhere to staff and equipment restrictions whilst ensuring patients are treated in a timely manner. We simultaneously address the master surgical scheduling problem (MSSP) and the surgical case assignment problem (SCAP). Stochastic surgical durations are considered in order to produce more robust schedules and reduce unexpected overtime. Also incorporated into the model are several constraints regarding patient wait targets that are set by the Australian government. The problem is formulated using a mixed integer nonlinear programming (MINLP) approach and solved using a variety of hybrid metaheuristics. The metaheuristics implemented are inspired by simulated annealing (SA) and reduced variable neighbourhood search (RVNS). In particular, we present an adaptive SA and hybridise SA and RVNS to greatly improve solution quality. The solution neighbourhoods used by the metaheuristics are based on the hierarchical structure present in the combined MSSP SCAP. We consider a case study of an Australian public hospital with a large surgical department and compare the performance of our model to historical data. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction The surgical department is often one of the largest sources of hospital costs and patient flow. Inefficient use of the operating theatres (OTs) and related resources can result in increased costs and a decrease in patient and staff welfare. The schedule used by hospital administrators can have a massive impact on patient outcomes. Although the scheduling of surgical procedures is a complex process, many surgical departments still rely on simple heuristics and general intuition to produce weekly or fortnightly OT schedules and some daily updates. Such simple techniques can result in overtime, bed shortages, unused resources or large waiting lists due to unnecessary improper utilisation of operating rooms (ORs). The administrative practices of surgical departments can have a large impact on hospital costs, patient outcomes and the overall efficiency of a hospital. As such, a large number of papers exist which investigate and analyse some of the most significant issues in OT planning. For a thorough review of the existing literature see [1–3].
∗
Corresponding author. E-mail address: [email protected] (E. Kozan).
http://dx.doi.org/10.1016/j.orhc.2016.07.002 2211-6923/© 2016 Elsevier Ltd. All rights reserved.
Van Riet and Demeulemeester [1] review efforts to plan ORs in the case of both elective and emergency surgeries with an emphasis on the trade-offs encountered between costs and waiting times. Ferrand et al. [2] also consider both elective and emergency cases in their survey of methodologies used to ensure ORs are both responsive and efficient. In their literature survey, Guerriero and Guido [3] discuss both optimisation and simulation models used for OT planning and scheduling problems. OT planning problems can exist at the strategic, tactical or operational levels. Decisions at the strategic level include determining how much time to allocate to each surgical specialty in order to optimise the mix of patients processed [4]. Tactical problems include the development of a MSS [5], or determining which time blocks to allocate to a particular surgeon [6] or surgical specialty [7]. Problems at the operational level are advanced and allocation scheduling of elective patients. Advanced scheduling is the assignment of a patient to a time block on a particular day. This is also known as the SCAP. Allocation scheduling is the scheduling of patients within those time blocks, also known as the surgical case sequencing problem (SCSP). There are three main policies used in OT scheduling at the tactical level. These are open scheduling, block scheduling and modified block scheduling [6]. An open scheduling policy allows surgeons to assign a surgical case to any day in the planning horizon. Block scheduling allocates time blocks to surgeons or
50
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
surgical specialties. These time blocks may only be used by the surgeon or surgical specialty that they are assigned to. In a modified block scheduling policy, time blocks may be left open for use or unused time can be reallocated to a different surgeon or specialty. A MSS is often used under a block scheduling strategy in order to assign time blocks to surgical specialties. There are many benefits that arise from the implementation of a MSS [8]. The cyclic nature of the schedule allows the coordination of staff and resources and can improve patient flow and throughput. The predictability of a repeated MSS can be incredibly useful when planning admissions and can help stabilise bed occupancy. It is possible to generate an optimal or near optimal MSS for a hospital, based on their particular requirements and goals, using a variety of different operations research techniques. Many hospitals repeat the same MSS over a period of several months. Agnetis et al. [9] analyse the trade-off between the organisational benefits of repeating a single weekly MSS against a schedule updated weekly to adapt to the current state of the surgical waiting lists. The authors consider three policies for the development of MSSs, each permitting different levels of variability from an original schedule. The authors found that allowing even a small number of deviations from a reference schedule greatly improves the performance of the MSS in terms of lateness and waiting time. Throughout the literature, it is common to incorporate constraints on the suitability of the ORs for each specialty, the availability of surgical teams, and post-surgical capacity (e.g. [7]). Other authors include constraints on the patient mix to match hospital requirements. Most of the models used to solve the MSSP have a single objective which is usually related to maximising throughput or utilisation. Some models throughout the literature have multiple objectives, commonly incorporated in a single weighted objective function. Whilst a large number of mathematical formulations and solution techniques exist to produce optimal OT schedules with respect to a variety of objectives, it is necessary to make these solutions applicable and accessible to hospital administration staff. In order to allow surgical staff and scheduling nurses to utilise schedules generated using operations research techniques, a number of decision support systems are presented in the literature (e.g. [10]). The advanced scheduling problem (the OT planning problem or SCAP) is the problem of assigning patients to time blocks over a one to two week time period. This can be done using an open scheduling, block scheduling, or modified block scheduling policy [11]. Advanced scheduling is frequently addressed in the literature and solved using a variety of both exact and heuristic approaches. Whilst some authors rely on exact approaches to solve small test instances, as the problem is NP-hard [12] it is necessary to use heuristics, meta-heuristics and hybrid heuristics to provide good feasible solutions in reasonable amounts of computational time. Column-generation based heuristics are often used to solve OT scheduling problems at the operational level (cf. Wang et al. [13] and Fei et al. [14]) Fei et al. [15] apply a branch-and-price approach to solve the SCAP. Many authors rely on simulation based models to investigate OT planning problems (e.g. [16,17]). The constraints and objective functions used to define OT practice can differ significantly between models. Many models present in the literature minimise expenses related to OT use. Other models prioritise patient welfare over expenses. Roland et al. [18] places particular significance on human resource constraints. The tactical and operational levels of OT planning and scheduling are often treated separately. A more effective procedure involves scheduling the tactical and operational levels
simultaneously. Agnetis et al. [19] integrate the two levels of planning and compare an integrated approach to the use of a decomposition approach. Whilst the decomposition approach allowed the authors to produce comparable solutions in far less time, the model had several simplifying assumptions that make this method unsuitable here. The authors did not consider specific surgeon availability or suitability, and incorporated uncertainty in surgical durations using added slack time in conjunction with expected durations. These assumptions make this approach unsuitable for implementation on the case study considered in Section 2. Aringhieri et al. [20] solve a binary integer programming (BIP) formulation of the MSSP and the advanced scheduling problem (SCAP) simultaneously by making use of the hierarchical structure present in OT schedules in a metaheuristic. We also consider the hierarchical structure when solving the MINLP formulation presented in Section 3. The model presented by Aringhieri et al. [20] considers deterministic surgical durations, whereas we consider stochastic durations. Aringhieri et al. [20] assume that the strategic decisions can be used as an input to the problem, particularly the number of time blocks assigned to each surgical specialty. We implicitly consider the surgical case mix planning problem when solving the combined MSSP SCAP to allow for variations in the MSS based on the current state of waiting lists. One of the major concerns currently facing Australian hospitals is the length of elective surgery waiting lists. The longer a patient waits for surgery, the more their quality of life is reduced. Long surgery waits can also worsen patient outcomes. Occasionally, patients are not seen on time as they are not fit for surgery at that stage, however most overdue patients are late for surgery due to insufficient resources or poor planning. A limited number of operations research techniques have been used to manage waiting lists and improve patient outcomes whilst adhering to resource availability constraints. Everett [21] provides a decision support tool based on a simulation model of elective surgery waiting lists for Australian public hospitals. The model developed by the author can be used in order to monitor the system, assist in strategic decisions or even surgical scheduling at the tactical and operational levels. The simulation model is quite simplistic to allow for the results to be interpreted by hospital staff with ease and does not include as much detail as we incorporate into the model presented in Section 3. Patient priority is considered by Min and Yih [22] whilst assigning elective patients to available surgical capacity. In their model, the authors assign appropriate patients in order of urgency. A common issue identified in waiting list management is the discrepancy in waiting times between urgency categories which often leads to excessive waiting for non-urgent elective patients. Min and Yih [23] use a patient priority which is a combination of urgency category and waiting time by formulating the problem as an infinite horizon Markov decision process. Marques et al. [24] address the SCAP for a small Portuguese public hospital under an open scheduling strategy. The authors aim to maximise surgical suite utilisation and the number of surgeries scheduled. The authors ensure that all deferred urgent and high priority surgeries are scheduled within the planning horizon and use the remainder of the time to schedule priority and normal surgeries according to a weighted objective function. In this paper we provide a MINLP formulation of the combined MSSP SCAP that maximises the number of surgeries scheduled. Whilst the above waiting list management policies are worth considering, in this paper we consider the policies imposed by the Australian government to ensure the model presented is applicable to Australian public hospitals. We ensure that certain government dictated patient-wait targets are satisfied on the final day of the planning horizon and adhere to resource availability and suitability
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
restrictions (including those on surgical staff). We also incorporate the government’s ‘‘treat in turn’’ policy in which 60% of eligible patients are treated in the order in which they arrived on the waiting list. The model presented in Section 3 contributes to the literature in a number of ways. We tackle all three levels of OT planning and scheduling problems by solving the MSSP and SCAP and, in doing so, implicitly solve the surgical case mix planning problem. We also use an integrated approach when solving the combined MSSP SCAP to produce better solutions than those obtained under a decomposition approach. The model we present is more detailed than those seen in the literature through inclusion of waiting list management policies, stochastic surgical durations, surgeon availabilities and capabilities. This allows for its implementation in a real hospital setting. The metaheuristics presented in Section 4 can be used to produce good feasible solutions to the combined MSSP SCAP in a reasonable amount of computational time. In particular, hybridised SA–RVNS has not previously been used on hospital planning problems. We are also able to maintain feasibility of a highly non-linear formulation through use of the solution neighbourhoods presented in Section 4.1. The surgical department considered in the case study is one of the largest in Australia. Due to the size of the case study, the real-life instances solved in Section 5 are far more computationally complex than the majority of instances previously considered in OT planning and scheduling. In Section 2 we present the case study, including details on the surgical department, length of surgery, recent waiting list data, and current scheduling practice. In Section 3 we present a MINLP formulation of the problem described within the case study. We discuss the metaheuristics considered, solution neighbourhoods, and upper bounds in Section 4. Computational results including metaheuristic parameter tuning and model performance are presented in Section 5. In Section 6 we present conclusions and discuss possible future work. 2. Case study In this section we consider the case study of a large Australian public hospital. This adult hospital has one of the largest surgical departments of public Australian hospital and as such, it is incredibly important that the surgical department runs efficiently. The throughput of the hospital is bound by a number of fixed resources. The hospital has 21 ORs, with one reserved for emergency surgeries. Elective surgeries are only scheduled on weekdays and never after 8pm. On the weekend, five emergency ORs are open. There are around 800 beds in the hospital, including 300 surgical beds. The surgical care unit (SCU) has around twenty beds available, as does the post anaesthesia care unit (PACU). As the hospital is a teaching hospital it is necessary to have surgeons of an appropriate skill level available to supervise the surgeries. Variable resources include the availability of surgeons, nurses, anaesthetists, etc. We only consider the availability of the assigned surgeons as there are sufficient numbers of nurses and anaesthetists available. The scheduling of nurses and anaesthetists is performed based on the MSS, but not simultaneously. The surgical nurses often specialise, whereas the anaesthetists are more generalised, although they may have specialty preferences. For a patient to enter the elective surgery waiting list, certain procedures must be followed. Initially, a patient must be referred to a specialist. After the patient has had a consultation the case is reviewed by the surgeon to ensure the patient does require the surgery. If these conditions have been met, the patient is put on the elective surgery waiting list and assigned an urgency category. There are three urgency categories used for classifying elective patients. A patient’s category is determined by the recommended
51
time before his/her elective procedure (30, 90 and 365 days). Another measure that the Queensland government has taken in order to effectively manage waiting lists is the use of a ‘‘treat in turn’’ model for patient care. The treat in turn model ensures that at least 60% of eligible patients treated in the planning horizon are among the set of patients who have been waiting the longest. A patient is treat in turn eligible if he/she is in category two or three (less urgent than a category one patient), a public patient (private hospitals do not have to adhere to this policy), and ready for care at the start of the planning horizon. The set of all treat in turn eligible patients is known as the treat in turn cohort. The following method is used by Queensland Health [25] to determine the proportion of patients treated in turn. First, the longest waiting cohort is calculated as the subset of treat in turn eligible patients who have been waiting longest. The number of patients in the longest waiting cohort is equal to the number of treat in turn eligible patients treated during the planning horizon. The treat in turn policy is satisfied as long as the number of longwait patients treated during the planning horizon is at least 60% of the number of treat in turn eligible patients treated, regardless of sequence within the planning horizon. Another policy implemented by the hospital is that all category one patients must be given a surgery date at the time they are added to the waiting list. To manage this, each surgical specialty has a case manager who liaises with surgical staff to ensure that patients are seen in a timely manner. Currently there are around 2900 patients on the elective surgery waiting list. A breakdown of a recent waiting list can be seen in Table 1. Staff at the elective surgery booking office, in conjunction with surgeon preferences, assign elective booking dates to patients whilst satisfying the MSS, patient wait targets and the treat in turn model. Currently the hospital operates under a four week rotating MSS with two blocks per OR per day (or a full day block). The use of a MSS ensures that there is consistency in staff schedules and predictable patient flow patterns for the downstream wards. Very rarely are major changes made to the MSS. The schedule has remained largely the same since around 2000. Specialties are able to request extra OR time in order to keep up with changes in demand. If an OR is available a cost analysis is performed and, if approved, the schedule is updated. In many cases, if there is a change in demand for a particular specialty it may mean that another specialty must reduce their hours or the extra surgeries may be outsourced to hospitals in the surrounding region. Minor modifications are often made as surgeon availabilities change. When updating or modifying the MSS, preference is given to full day blocks over half day blocks as there is no need to stop half way through the day. The length of time allocated to each time block is dependent on the specialty that will be using the OR and the availability of visiting medical officers (VMOs). Some specialties require longer working days. Cardiology often requires ten hour days whereas ENT may require a twelve to fourteen hour block for a major reconstructive case. Each Thursday the hospital has an elective bookings meeting to discuss the patients planned for the following week. Some surgeons request that a time block be left empty as they predict that there will be admissions over the weekend that require OR time. If a surgical team is unavailable or does not need their usual OR time, other specialties are able to request the OR during that time. VMOs are also able to request OR time. One difficulty in creating the MSS is the availability of consultants. Many consultants spend the majority of their time working at surrounding private hospitals and are only available for a few days each month. It is important to be able to incorporate the required OR time for each of the consultants whilst satisfying their availabilities. Some specialties do not require that a consultant is present for surgeries, and instead a registrar can perform a list of surgeries for several consultants.
52
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 1 Recent elective surgery waiting list. Urgency category
1
Specialty
On time
Overdue
2 On time
Overdue
On time
Overdue
Cardiology (CARD) Ear, Nose and Throat (ENT) General Surgery (GEN) Gynaecology (GYN) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Urology (UROL) Vascular Surgery (VASC) Total
4 10 80 0 12 12 22 58 58 23 279
0 0 0 0 0 0 0 0 0 0 0
46 46 183 1 79 194 163 213 143 70 1138
0 17 29 0 10 10 26 50 2 5 149
1 12 74 0 4 760 310 70 118 13 1362
0 0 1 0 0 6 2 3 3 1 16
Leading up to surgery, there are a number of patient pathways. The patient must attend the preadmission clinic at the hospital before his/her surgical date. On the day of surgery, elective patients arrive at the hospital and inpatients are allocated to a ward. Depending on the time until surgery and any other procedures required beforehand, the patient is either sent to the ward or to the SCU to prepare for surgery. After surgery, all patients are sent to the PACU to recover; this usually takes four to five hours. After the initial recovery period, day patients recover further in the SCU and inpatients are sent back to the ward or intensive care unit (ICU). Any surgical transfers from the ICU are returned directly to ICU after surgery. Whilst there are a limited number of beds in the SCU and PACU, as the surgical department has priority in booking these beds, there are rarely bottlenecks. Occasionally downstream wards are put under pressure due to the large number of surgical patients being transferred. Although downstream capacity is out of the scope of this project, it would be extremely useful to incorporate into future work. The hospital has over 50 surgical specialties, however when surgeries are booked they are usually booked under nine to twelve main specialties. Table 2 shows the count, mean, and variance of the surgical durations (in min) for elective surgeries performed over 2012 and 2013. More recent data is unavailable. From Table 2 it can be seen that specialties that perform more complex procedures have a larger variance than specialties that perform less complicated procedures. Whilst nine to twelve specialties are usually used when booking at this hospital; surgical data is reported under other specialties as well. This difference in reporting systems is limiting as we are only able to schedule under the specialties that the elective surgery booking office uses on the waiting lists. The MSS used by the hospital is able to be more specific. Given that the hospital under study is a teaching hospital, there tend to be variations in surgical durations based on the experience of the surgeon or the severity of the case. When surgeries run overtime it can result in large overtime costs or even the cancellation of surgeries. Whilst surgeons can provide estimate surgical durations on a case by case basis this data is not available for use in the model. Kayış et al. [26] present a statistical model to estimate surgical durations based procedure type, staff, patient type and when the procedure occurs. The authors found that by hybridising their work with that of Dexter et al. [27] they are able to predict surgical durations without bias and with a high accuracy. Due to a lack of appropriate data, we must consider a more simple approximation when incorporating stochasticity in the model presented in Section 4. Strum et al. [28] showed that surgical durations can be modelled using a lognormal distribution. Stuart and Kozan [29] use a lognormal approximation to the sum of lognormally distributed
3
random variables in order to ensure a 15.87% level of accuracy. The authors consider only a single OR. We tighten that bound to ensure that there is a probability of less than 5% that the surgeries assigned to a particular block run overtime and extend the work by considering multiple ORs. The lack of stochasticity is not a limitation in the case of non-elective surgeries as there is an OR reserved for these cases, although these cases often overflow into underutilised OR blocks. Non-elective surgeries consist of both emergency arrivals and other urgent surgeries from within the hospital. Heng and Wright [30] found that the use of a dedicated OR significantly reduces the number of elective cancellations and improves the waiting times of non-elective cases. Like many hospitals, certain ORs are more suited for use by certain specialties. For example, there is extra equipment in the vascular OR; however this does not prohibit other specialties from using the OR. It is a matter of convenience that the hospital usually assigns certain specialties to each OR. This reduces the movement of valuable equipment between sessions. Another consideration that administrative staff make when adjusting the MSS is the availability of equipment. The hospital has a single machine that can be used by several specialties, it is important that two specialties are not simultaneously scheduled to use the same equipment. Administrative staff also attempt to have the equipment used as frequently as possible. Fig. 1 (see Appendix) shows the surgeries performed in the twenty ORs under study in the week from which the elective surgery waiting list was used in Section 5.1. The figure shows the specialties that used the time block, the number of elective and emergency surgeries performed, and the total time required for the surgeries. For example, in OR one on the Monday morning there was a single elective surgery, no emergency surgeries, and the total time was 3.08 h. In this week, a total of 266 surgeries were performed of which 218 were elective surgeries. Future work includes incorporating emergency arrivals in the model as the hospital is unable to complete all emergency procedures in the one reserved OR. To test whether the assumption of lognormal surgical durations holds, the data is first transformed by taking the natural logarithm of the surgical durations. We then test whether the transformed data is likely to belong to a normal distribution. R© version 3.1.1 is used during statistical analysis. Data is visually analysed using density plots and normal Q–Q plots. Shapiro–Wilk, Anderson–Darling and Lilliefors normality tests are used to determine the suitability of the distribution. The results from these tests are shown in Table 3. From Table 3 it can be seen that the assumption of lognormal surgical durations is inappropriate for most specialties according to the goodness-of-fit tests used. This may be due to the large (n > 100) and heterogeneous populations. Within the model, the lognormal approximation is only used to determine how many patients can be scheduled within a
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
53
Table 2 Reported elective length of surgery data. Specialty
2012 Time (min) Count
Average
Variance
Count
Average
Variance
Acute Surgical Unit (ASU) Breast and Endocrine (B&E) Cardiology (CARD) Colorectal (COLO) Cardiac Surgical Unit (CSU) Dental Surgery (DENT) Ear, Nose and Throat (ENT) Faciomaxillary (FMAX) Gastroenterology (GAS) Gynaecology (GYN) Hepato-Pancreato-Biliary (HPB) Liver Transplant (LTPT) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Persistent Pain Management Service (PPMS) Psychiatry (PSY) Respiratory (RESP) Renal Transplant (RTPT) Upper GI and Soft Tissue (UGI) Urology (UROL) Vascular Surgery (VASC)
–
– 122.31 0.00 290.91 260.41 75.92 336.44 150.29 – 34.34 223.10 53.03 214.85 54.70 172.70 129.90 – – – 114.72 149.79 95.28 179.64
– 2996.82 0.00 55 744.56 7426.43 76.93 105 600.21 16 927.13 – NA 14 516.63 848.91 15 944.61 1477.97 15 385.07 21 533.34 – – – 5467.35 15 159.66 6826.57 8819.81
23 540 – 359 843 4 471 338 2 4 355 6 451 1988 1735 1622 16 130 91 572 518 1472 528
136.31 108.54 – 173.19 266.83 71.46 254.35 99.13 53.53 33.84 182.60 64.81 190.77 47.55 130.59 85.77 34.46 13.51 45.41 93.76 133.20 86.53 152.31
4446.71 2709.94 – 20 950.38 8375.50 1486.97 76 440.92 5574.57 72.72 36.03 12 292.92 900.69 15 277.17 1318.96 7559.47 7736.56 327.36 12.56 582.69 3502.18 12 647.82 6898.36 7149.21
422 88 524 1163 6 806 228 – 1 484 2 610 1941 1507 1933 – – – 401 541 1884 437
2013 Time (min)
Table 3 Statistical analysis on surgical durations. Specialty
Observations
S.W. p-val
A.D. p-val
Lillie p-val
Acute Surgical Unit (ASU) Breast and Endocrine (B&E) Colorectal (COLO) Cardiac Surgical Unit (CSU) Dental Surgery (DENT) Ear, Nose and Throat 1 (ENT1) Ear, Nose and Throat 2 (ENT2) Ear, Nose and Throat 3 (ENT3) Faciomaxillary (FMAX) Gastroenterology (GAS) General (GEN) Gynaecology (GYN) Hepato-Pancreato-Biliary (HPB) Liver Transplant (LTPT) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Renal Transplant (RTPT) Upper GI and Soft Tissue (UGI) Urology (UROL) Vascular Surgery (VASC)
23 540 359 843 4 34 29 3 338 2 1740 4 355 6 451 1998 1735 1622 572 518 1472 528
0.517 0.370 0.000 0.000 0.323 0.013 0.283 0.400 0.061 NA 0.000 0.864 0.000 0.830 0.001 0.000 0.000 0.000 0.001 0.000 0.000 0.005
0.352 0.131 0.000 0.000 NA 0.014 0.328 NA 0.139 NA 0.000 NA 0.000 NA 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.730 0.063 0.000 0.000 NA 0.037 0.365 NA 0.259 NA 0.000 NA 0.000 0.684 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.001
time block such that in 95% of cases the total duration of their surgeries does not exceed the length of the time block. Although the lognormal distribution appears to be a poor fit to the model, random sampling and simulation from actual data can be used to verify the maximum number of patients that can be scheduled in each time block, under each specialty. We verify the suitability of the lognormal approximation by randomly sampling (according to a uniform distribution) 50,000 sets of patients for each combination of set size (between one and twenty patients) and specialty. The sum of surgical durations is calculated for each of the sets and the 95th percentile is found. Based on this we then determined the number of patients that can be scheduled in each type of time block. Table 4 shows the maximum number of patients that can be scheduled in each time block under each specialty as well as the lognormal distribution parameters calculated using historical data. Note that an asterisk beside a surgical specialty indicates that the lognormal distribution parameters from the General
specialty were used as there was insufficient data. The values within this table are the same whether calculated using lognormal approximations or random sampling. Thus, we are confident in the use of lognormal approximations for surgical duration in the model, despite poor goodness-of-fit test results. 3. Model formulation The combined MSSP SCAP is formulated using MINLP to allocate surgeries, surgeons and specialties to time blocks. 3.1. Scalar parameters Bmax : number of time blocks in a working day Cmax : number of patient urgency categories Dmax : number of working days in the planning horizon Hmax : number of surgeons that practice at the hospital (including VMOs)
54
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 4 Surgical duration model parameters. Specialty *
Anaesthetics (ANAE) Acute Surgical Unit (ASU) Breast and Endocrine (B&E) Cardiothoracic (CARD)* Colorectal (COLO) Cardiac Surgical Unit (CSU) Dental Surgery (DENT) Ear, Nose and Throat 1 (ENT1) Ear, Nose and Throat 2 (ENT2) Ear, Nose and Throat 3 (ENT3) Faciomaxillary (FMAX) Gastroenterology (GAS) General (GEN) Gynaecology (GYN) Hepato-Pancreato-Biliary (HPB) Liver Transplant (LTPT) Neurosurgery (NSUR) Ophthalmology (OPHT) Orthopaedic (ORTH) Plastic Surgery (PLAS) Renal Transplant (RTPT) Trauma (TRMA)* Upper GI and Soft Tissue (UGI) Urology (UROL) Vascular Surgery (VASC) Visiting Medical Officers (VMOs)* *
Mu
Sigma2
Patients per half day
Patients per full day
0.788 0.702 0.478 0.788 0.783 1.426 0.036 1.200 0.498 0.585 0.266 −0.137 0.788 −0.598 0.945 −0.031 0.970 −0.474 0.583 −0.015 0.267 0.788 0.516 0.027 0.786 0.788
0.395 0.217 0.209 0.395 0.534 0.112 0.258 0.513 0.779 0.182 0.453 0.025 0.395 0.031 0.317 0.196 0.353 0.463 0.370 0.724 0.338 0.395 0.542 0.658 0.271 0.395
0 1 1 0 0 0 2 0 0 1 1 5 0 8 0 3 0 3 1 1 1 0 0 1 0 0
1 2 3 1 1 1 6 0 1 3 3 10 1 16 1 7 1 8 2 3 4 1 2 3 2 1
Indicates lognormal distribution parameters from the General specialty were used as there was insufficient data.
Pmax : number of patients on the waiting list at the start of the planning horizon Smax : number of surgical specialties (where surgeons of specialty Smax are VMOs) Tmax : number of operating rooms M: the largest possible number of days overdue for any patient at the end of the planning horizon Mp : the maximum possible number of patients that can be scheduled in a full day time block π : the proportion of eligible patients that must be treated in turn. 3.2. Index sets B: the set of time blocks in a working day. B = {1, . . . , Bmax } C : the set of patient urgency categories. C = {1, . . . , Cmax } D: the set of working days in the time horizon. D = {1, . . . , Dmax } H: the set of surgeons that practice at the hospital (including VMOs). H = {1, . . . , Hmax } S: the set of surgical specialties (where specialty Smax is for VMOs). S = {1, . . . , Smax } T : the set of operating rooms. T = {1, . . . , Tmax } P: the set of patients on the waiting list at the start of the planning horizon. P = {1, . . . , Pmax }. 3.3. Vector parameters
αc : the proportion of patients in urgency category c that must be seen on time, ∀c ∈ C δp : the number of days until patient p is due at the start of the period, ∀p ∈ P Wp : the number of days that patient p will have been on the waiting list at the end of the period, ∀p ∈ P τb : the length of time block b in hours, ∀b ∈ B Ns− : the minimum number of OR hours assigned to specialty s over the time horizon, ∀s ∈ S Ns+ : the maximum number of OR hours assigned to specialty s over the time horizon, ∀s ∈ S
Lp : the of patient p’s surgery, ∀p duration ln N µp , σp2
βpc =
∈
1 if patient p is in urgency category c 0 otherwise,
∀p ∈ P , c ∈ C 1 if patient p is in the treat in turn cohort λp = 0
otherwise,
∀p ∈ P 1 if patient p can be treated by surgeon h Eph = 0
otherwise,
∀p ∈ P , h ∈ H 1
Fhdb =
if surgeon h is available on day d during block b 0 otherwise,
∀h ∈ H , d ∈ D, b ∈ {1, 2} 1 if surgeon h is a member of specialty s Ghs = 0
otherwise,
∀h ∈ H , s ∈ S 1
if operating room t is equipped for surgeries by specialty s 0 otherwise,
Rts =
∀t ∈ T , s ∈ S . 3.4. Decision variables 1
if room t is used by specialty s during time block b on day d 0 otherwise,
Xtsdb =
∀t ∈ T , s ∈ S , d ∈ D , b ∈ B 1
Ytpdb =
if room t is used for patient p during time block b on day d 0 otherwise,
∀t ∈ T , p ∈ P , d ∈ D, b ∈ B
P. Lp
∼
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
1
if room t is used by surgeon h during time block b on day d 0 otherwise,
Zthdb =
A surgeon can only be assigned an OR if their specialty is also assigned to that OR during the time block.
∀t ∈ T , h ∈ H , d ∈ D, b ∈ B
Zthdb ≤
if patient p has on time surgery, or is on time on the final day 0 otherwise,
Xtsdb ≤ Rts ,
otherwise,
∀p ∈ P . Utdb : the expectation of the sum of surgical durations assigned to room t on day d in block b Vtdb : the variance of the sum of surgical durations assigned to room t on day d in block b.
∀t ∈ T , s ∈ S , d ∈ D, b ∈ B.
(8)
Ytpdb eµp
eσp , 2
∀t ∈ T , d ∈ D, b ∈ B.
(9)
Constraint (10) is used to calculate the variance of the sum of surgical durations assigned to each time block. Vtdb =
Ytpdb eσp − 1 e2µp +σp , 2
2
p∈P
3.6. Constraints
∀t ∈ T , d ∈ D, b ∈ B.
It is not possible to share a block between different specialties. There must be at most one specialty assigned to each block.
∀t ∈ T , d ∈ D, b ∈ B.
(1)
Each OR can be assigned in both morning and the afternoon or for a full day only. Xtsd1 + Xtsd2 + 2Xtsd3 ≤ 2,
∀t ∈ T , d ∈ D.
(2)
s∈S
(10)
Constraint (11) ensures a 95% confidence level that surgeries assigned to each time block will be completed on time, given the above approximations.
s∈S
(7)
p∈P
The objective of the static model is to maximise the number of surgeries performed planning horizon. during the Maximise t ∈T p∈P d∈D b∈B Ytpdb .
Xtsdb ≤ 1,
∀t ∈ T , h ∈ H , d ∈ D, b ∈ B.
The duration of each surgery is stochastic and can vary due to a number of different factors. Based on the work by Strum et al. [28], Stuart and Kozan [29,31] approximate the sum of lognormally distributed surgical durations using a lognormal distribution. Constraint (9) is used to calculate the expectation of the sum of surgical durations assigned to each block. Utdb =
3.5. Objective function
Xtsdb Ghs ,
An OR can only be used by a specialty if that OR is properly equipped for use by that specialty.
∀p ∈ P 1 patient p is a member of the longest waiting cohort θp = 0
s∈S
1
γp =
55
2 Utdb
2 Vtdb + Utdb
Vtdb +U 2
1.64485×ln
e
U2 tdb
tdb
≤ τb ,
∀t ∈ T , d ∈ D, b ∈ B.
(11)
A patient’s surgery can be performed at most once during the planning horizon.
Each specialty must not be allocated more than its maximum number of OR hours over the planning horizon.
Ytpdb ≤ 1,
∀p ∈ P .
(3)
t ∈T d∈D b∈B
∀s ∈ S .
(12)
t ∈T d∈D b∈B
If a surgeon is a member of a specialty at the hospital then they cannot share a time block with any other surgeons. It is possible for VMOs to share a time block amongst themselves. This reflects current hospital practices.
τb X tsdb ≤ Ns+ ,
Zthdb ≤ 1 + XtSmax db
h∈H
τb X tsdb ≥ Ns− ,
∀s ∈ S .
(13)
t ∈T d∈D b∈B
GhSmax − 1 ,
h∈H
∀t ∈ T , d ∈ D, b ∈ B.
Each specialty must be allocated at least its minimum number of OR hours over the planning horizon.
(4)
A surgeon can only be assigned a time block if he/she is available during that time block. A surgeon cannot be assigned to a morning and/or afternoon block and a full day block on the same day. A surgeon can only be assigned to a full day block if they are available in both the morning and the afternoon.
Constraint (14) determines whether a patient’s surgery was performed on time, or if the patient is not overdue at the end of the planning horizon.
δp − (Dmax − 1) 1 −
Ytpdb
t ∈T d∈D b∈B
−
(d − 1)Ytpdb ≥ M γp − 1 ,
∀p ∈ P .
(14)
t ∈T d∈D b∈B
Zthdb + Zthd3 ≤ Fhdb ,
∀h ∈ H , d ∈ D, b ∈ {1, 2}.
(5)
t ∈T
A patient can only be treated by a surgeon if the surgeon is capable of performing that particular surgery. This constraint is especially necessary in teaching hospitals. Ytpdb ≤
h∈H
Zthdb Eph ,
∀t ∈ T , p ∈ P , d ∈ D, b ∈ B.
A certain proportion of patients in each category must be seen on time. This proportion is in line with Australian government standards.
p∈P
(6)
γp βpc ≥ αc
βpc ,
∀c ∈ C .
(15)
p∈P
A patient can only be in the longest waiting cohort if they are eligible for the treat in turn policy. This means a patient must be
56
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
ready for care at the start of the month, a public patient, in urgency category two or three and assigned to a surgical specialty.
θp ≤ λp ,
∀p ∈ P .
(16)
The size of the longest waiting cohort is equal to the number of treated in turn eligible patients treated as elective patients during the planning horizon.
θp =
p∈P
Ytpdb λp .
(17)
t ∈T p∈P d∈D b∈B
The longest waiting cohort must be composed of the patients that have been on the waiting list for the longest amount of time. max p∈p
1 − θp Wp λp ≤ θp Wp + 1 − θp M .
(18)
The proportion of longest waiting patients that were treated in turn must meet government standards.
θp Ytpdb ≥ π
t ∈T p∈P d∈D b∈B
θp .
(19)
p∈P
4. Solution approaches As the MSSP and SCAP are both NP-hard, the computational time required increases exponentially with an increase in patients, surgeons, rooms, etc. It is necessary to use heuristics to obtain good, practical solutions in a reasonable amount of time. The use of heuristics often allows for simple modifications to adapt to changes in problem specifications. Construction heuristics are used to find an initial, feasible solution. Hans et al. [32] present a robust formulation of the SCAP in which the sum of surgical durations is assumed to be normally distributed. The authors implement a number of construction heuristics for the SCAP including a first fit approach, longest processing time, and random sampling. Once an initial solution is found, improvement heuristics are used to improve initial solutions. Hans et al. [32] use SA with random swaps based on either inserting a random patient in a suitable random day, or swapping two suitable surgeries. The authors found that SA performed well, but required large amounts of computational time. Agnetis et al. [19] implement a decomposition approach to the combined MSSP SCAP in which they formulate and solve the SCAP as a number of knapsack problems. The MSSP formulation considered by the authors is somewhat simplified as the authors do not consider stochastic surgical durations. This simplifies the implementation of the knapsack formulation. Riise and Burke [33] present metaheuristics to concurrently solve the SCAP and SCSP. The authors use steepest variable neighbourhood search (SVNS) as a comparative baseline for their iterative local search method which utilises variable neighbourhood descent. Landa et al. [34] also consider the SCAP and SCSP. The authors implement a hybrid optimisation algorithm, combining both Monte Carlo simulation and neighbourhood search. The algorithm used consists of two phases, the first to solve the SCAP and the second to address the SCSP. The authors also consider stochastic surgical durations. Other metaheuristics frequently implemented to solve OT planning problems include ant colony optimisation (ACO) [35], tabu search (TS) [20] and genetic algorithms (GA) [18]. The choice of heuristic is largely dependent on the formulation of the problem and the size of the problem instances. As the public hospital used in the case study has a large surgical department with a long waiting list it is inappropriate to evaluate all neighbouring solutions or to use population based metaheuristics (ACO, GA, etc.) as this would require excessive
amounts of memory on a real-life instance. For this reason, we must choose metaheuristics that can be implemented without large memory requirements. Not only is SA simple to implement, it has been used on a variety of hospital planning problems (e.g. [32]) and frequently outperforms other metaheuristics including ACO, TS and GA (cf. Birattari [36]). SA has also presented some of the best known solutions to a number of combinatorial optimisation problems [37]. Due to its very good performance on a wide variety of problems and its appropriateness with regards to the problem structure, SA (in conjunction with an improved version) is suitable for use as a benchmark metaheuristic on this problem. RVNS is well suited to large nonlinear optimisation problems. RVNS has been used on a number of hospital planning problems (e.g. [38]). Although the optimality gap increases with increasing problem size, RVNS can be used to find solutions to very large MINLP problems whereas commercial solvers often fail at medium instances. Both SA and RVNS are suited to the problem structure and size (in terms of memory requirements), easy to implement, and have been shown to outperform other common metaheuristics on a wide variety of optimisation problems. As such, we use these metaheuristics to compare performance with the hybridised metaheuristics presented. Future work includes the incorporation of other suitable metaheuristics including TS. We also present a modified SA in which different neighbourhood structures are used with certain probabilities over time. SA and RVNS are hybridised to utilise the best properties of each of the metaheuristics. To the best of our knowledge hybridised SA–RVNS has not been implemented on OT planning/scheduling problems or MINLP of this magnitude. As such, the use of hybridised SA–RVNS on this large MINLP OT planning and scheduling problem is quite innovative. In this section we present a number of metaheuristics designed to produce good and feasible solutions to the MINLP formulation given in Section 3. We begin by presenting the solution neighbourhoods that are implemented in conjunction with the metaheuristics. 4.1. Solution neighbourhoods The difficulty in implementing certain metaheuristics on other, more general, MINLP models lies in the definition of solution neighbourhoods specific to those problems. Here, we consider three levels of neighbourhood structure. These levels are specialty, surgeon and patient. In order to maintain feasibility, we must consider the implications when we change one of these levels. If the specialty assigned to a block is changed, we must also change the surgeon and patients assigned to that block to ensure constraints (6) and (7) still hold. Similarly, if we only change the surgeon assigned to the block, we may also be required to change the patients assigned. Switching patients does not affect either of the upstream levels. In assigning new specialties, surgeons and patients, we must ensure that availability and suitability constraints are upheld. We also consider neighbourhood changes that swap a specialty from a full day block to two half day blocks and vice versa. These neighbourhood structures ensure that the solution produced at each of the iterations satisfies constraints (1)–(12), (14) and (16)–(19). Depending on the composition of the waiting list, availability of surgeons and suitability of ORs, it may not be possible to satisfy constraints (13) and (15). As such, we cannot guarantee that the following neighbourhood swaps will produce solutions feasible with respect to constraints (13) and (15). We only ever accept a move if the new solution is feasible with respect to all
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
constraints presented in Section 3. Thus, if it is possible to satisfy all constraints, the solutions found using the metaheuristics will be feasible. When finding a random neighbouring solution we perform one of five neighbourhood swaps described below. 1. Swap specialty, surgeon and patients: – Randomly select an OR, day and block. – Remove the patients, surgeon(s) and specialty (if any) currently assigned to that OR during the selected time block. – Randomly select a suitable specialty and surgeon(s) from the new specialty. Assign these to the selected time block. – Randomly assign appropriate patients to the time block without violating the constraints presented in Section 3. 2. Swap surgeon and patients: – Randomly select an OR, day and block. – Remove the patients and surgeon(s) (if any) currently assigned to that OR during the selected time block. – Randomly select a suitable surgeon. Assign the new surgeon to the selected time block. – Randomly assign appropriate patients to the time block without violating the constraints presented in Section 3. 3. Swap patients: – Randomly select an OR, day and block. – Remove the patients (if any) currently assigned to that OR during the selected time block. – Randomly assign appropriate patients to the time block without violating the constraints presented in Section 3. 4. Select a full day block and swap to two half day blocks: – Randomly select an OR and day that is currently assigned to a single specialty and surgeon for the entire day. – Remove the patients, surgeon and specialty currently assigned to that OR during the time block. Reassign the same surgeon and specialty to both the morning and afternoon blocks on the same day in the same OR. For both the morning and afternoon time block, randomly assign appropriate patients to the time blocks without violating the constraints presented in Section 3. 5. Select two half day blocks of the same specialty and swap to a full day block: – Randomly select an OR and day that is currently assigned to the same surgeon across both the morning and afternoon time blocks. – Remove the patients, surgeon and specialty assigned to the half day blocks. Reassign the same surgeon and specialty to a full day block in the same OR on the same day. – Randomly assign appropriate patients to the new full day block without violating the constraints presented in Section 3. In addition to the above neighbourhood swaps, after every N iterations the solution is ‘cleaned’ by grouping scheduled patients that can be seen by the same surgeon and reassigning these patients where appropriate in order to reduce the number of time blocks used. No patients are added or removed from the schedule; the assignment of patients is simply rearranged. As the current subset of patients is feasible with respect to all constraints, the ‘cleaned’ solution will also be feasible. This assists in improving the quality of the solution obtained by each of the metaheuristics. Surgeons and specialties to be placed in the schedule (swap type one or two) are not selected according to a uniform distribution. Suitable specialties are selected according to a probability proportional to the maximum number of patients of that specialty that could fit within the time block. Suitable surgeons are selected according to a probability proportional to the number of available patients that can be scheduled at that iteration.
57
4.2. Simulated annealing SA is a global optimisation metaheuristic that was introduced by Kirkpatrick et al. [39]. Based on the idea of temperature control during the cooling of metals, SA allows a worse solution to be accepted with a probability that decreases over time. The acceptance of worse solutions allows SA to escape local optimum in the search for the global optimum. SA is often used on hospital planning problems due to its ease of implementation and solution quality. Ceschia and Schaerf [40] address the patient admission planning problem with the inclusion of OR utilisation constraints. The authors use a metaheuristic based on SA with particular focus on the use of appropriate structures. On a similar problem Ceschia and Schaerf [37] used this method produced best known solutions to a number of sample instances. Fügener et al. [41] implement a number of metaheuristics including SA to produce a MSS for a small Dutch hospital. The authors consider the impact of the MSS on the occupancy of a number of downstream wards. The authors found that SA performed well in terms of both solution quality and time required. Birattari [36] discusses a number of parameter tuning algorithms including a brute force approach and the use of the F-Race algorithm. The author also investigates the use of a variety of metaheuristics on several instances of well-known combinatorial optimisation problems. After careful parameter tuning, SA tends to outperform other metaheuristics including ACO, GA and TS. Due to its ease of implementation, suitability for large problem instances and quality of solutions, as seen in the work referenced above, we use SA as a baseline metaheuristics to compare with the performance of the other metaheuristics on the model presented R in Section 3. We implement SA in MATLAB⃝ through the use of the solution neighbourhoods presented in Section 4.1. At each iteration we randomly select a neighbourhood structure and find a neighbouring solution. We accept solutions according to the acceptance probability; the best solution is only updated if the solution is both feasible and improved. The acceptance probability is given by the expression exp((new obj − old obj)/T ), where T is the metaheuristic temperature at that iteration. 4.3. Adaptive simulated annealing In order to produce improved solutions, we modify our SA based metaheuristic to a more adaptive form. The acceptance probability of a new solution is dependent on both the change from the previous solution and on the time at which the solution is found. This ensures a smoother transition between acceptance probabilities as time progresses. The acceptance probability is given by the expression exp((new obj − old obj) × toc kp /T ), where T is the metaheuristic temperature, toc is the time since starting the heuristic and kp is the exponent parameter. Whilst the selection of neighbourhood structure in SA is done using a uniform distribution, we increase the probability that swaps of type four and five occur over time as these swaps are more useful in a mostly full schedule. Initially, swaps of type four and five cannot be selected. After just over 75% of the iterations have elapsed the swap type is selected according to a uniform distribution. The probability of selecting swaps one, two, or three is given by Eq. (20) and the probability of selecting swaps four or five is given by Eq. (21).
−1
(20)
.
(21)
p1,2,3 (t ) = 3 + 2 (1 + exp (−t + 0.75 × Imax ))−1
p4,5 (t ) = (3 (1 + exp (−t + 0.75 × Imax )) + 2)
−1
4.4. Reduced variable neighbourhood search Variable neighbourhood search (VNS) is a metaheuristic proposed by Mladenović and Hansen [42]. Originally used to solve
58
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
combinatorial optimisation problems, VNS has since been successfully extended to MILPs and some MINLPs. VNS utilises neighbourhood structures, local search, and perturbation to produce good feasible solutions to optimisation problems. A number of variations of VNS exist which are more suited to certain problem structures. A detailed introduction to VNS can be found in [43]. A number of different hospital planning problems including the nurse rostering problem [44] and the scheduling of patients within a hospital [38] have been solved using VNS inspired metaheuristics. Recently, Aringhieri et al. [45] solve the SCAP under a block scheduling policy with the aim of levelling bed occupancy in downstream wards. The authors implement VNS, RVNS and an adaptive version of the algorithm. Depending on problem formulation and metaheuristic parameters used, VNS can produce solutions as good as exact MINLP solvers on small to medium instances (e.g. [46]). Whilst the optimality gap increases with the size of the problem, VNS can provide good feasible solutions to realistic instances in short amounts of time. On the other hand, commercial MINLP solvers often fail on realistic instances. Due to the size of our problem, we consider RVNS in which a random neighbour is selected and chosen if the objective is improved. This vastly reduces the computational time and storage required and makes the work accessible to those without high performance computers. 4.5. Hybridised simulated annealing and reduced variable neighbourhood search There are limitations associated with both SA and RVNS that can reduce the average quality of solutions found. By randomly selecting solution neighbourhoods, SA is not able to make use of the neighbourhood structures present in the combined MSSP SCAP. It is a common occurrence for RVNS to become ‘stuck’ in a local optima due to its inability to accept worsening moves. To avoid each of these issues, we combine the beneficial features of both adaptive SA (cf. Section 4.3) and RVNS (cf. Section 4.4). In order to make use of the neighbourhood structures present in the problem, we continue to use the framework provided by the RVNS metaheuristic. In particular, we do not randomly select a neighbourhood type, but begin at the first neighbourhood structure and move down the list of structures if the neighbouring solution is not accepted. If the neighbouring solution is accepted then we move back up the list of structures. Unlike RVNS, hybridised SA–RVNS occasionally accepts nonimproving moves in an attempt to escape local maxima. As in the adaptive SA metaheuristic, the acceptance probability is both change and time dependent. Throughout the literature there are a few instances of hybridised SA–RVNS being implemented to solve a variety of optimisation problems. To the best of the authors’ knowledge, hybridised SA–RVNS has not been implemented on a MINLP model of this magnitude. The hybridised SA–RVNS algorithm is as follows: SA–RVNS AlgorithmStep 1 Initialise the algorithm: 1.1 Start with an initial solution. 1.2 Check feasibility. 1.3 If the initial solution is feasible, update the best solution to be the initial solution. Otherwise, let the best objective be zero. 1.4 k ← 1. Step 2 While the iteration limit has not been exceeded, iterate over the following steps:
2.1 Find a random neighbour, from neighbourhood k, of the previous solution. 2.2 Calculate the number of patients scheduled in the neighbouring solution. 2.3 Check the feasibility of the neighbouring solution. 2.4 If the acceptance probability is greater than some uniform random number and the neighbouring solution is feasible then update the latest solution to the neighbour found in step 2.1. 2.4.1 If the latest solution is also the best so far, update the best solution. 2.5 If the acceptance probability was too small, or the solution was not feasible, k ← k + 1 2.6 If k > Kmax, k ← 1. 4.6. Upper bounds In this section we present the manner in which upper bounds for the maximisation problem are found. We apply some simplifications to the model that allow us to solve a relaxed version using IBM ILOG CPLEX Optimization Studio 12.6. Due to the availability of data we allocate the same lognormal distribution parameters to the duration of all types of surgery within a specialty. As we only consider nine specialties in this section, it is possible to calculate the maximum number of patients scheduled in each type of block, whilst satisfying constraints (9)– (11). We consider both half and full day blocks which are five and ten hours respectively. This data is shown in Table 4. By utilising the data in Table 4 we are able to remove (highly nonlinear) constraints (9)–(11) from the model. Constraint (22) can be used as a simple replacement. Note that this substitution of constraints is only simple due to the assumption that surgical durations within the same specialty are independent and identically distributed. Asb : the maximum number of patients of specialty s that can be seen in block type b with less than 5% chance of exceeding the length of the time block. ∀s ∈ S , b ∈ B
p∈P
Ytpdb ≤
Asb Xtsdb ,
∀t ∈ T , d ∈ D, b ∈ B.
(22)
s∈S
The constraints used to address the treat in turn policy (constraints (16) through (19)) result in a quadratic constraint that does not satisfy the positive semi-definite property required by CPLEX. It is possible to introduce a set of constraints to linearise (19); however this results in an additional Pmax binary decision variables. Instead, we remove constraints (16) through (19) as it is too computationally expensive to include these when solving realistic instances. Due to the complexity of the problem it is also necessary to relax the binary decision variable Ytpdb such that Ytpdb ∈ [0, 1]. In summary, to calculate the upper bound on the randomly generated instance considered in Section 5.2 we allow the variable Ytpdb to be a real number on the interval [0, 1] and we include only constraints (1)–(8), (12)–(15), and (22). This may result in noninteger solutions that do not satisfy the treat in turn policy. After implementing these strategies we are able to solve each instance in CPLEX to produce an upper bound on the optimal solution to the model presented in Section 3. 5. Computational results In the following section we present computational results including metaheuristic parameter tuning, the performance of each of the metaheuristics is compared and the performance of the model presented in Section 3 is compared to recent historical data. All computational experiments in Sections 5.1 and 5.2 were R R performed using MATLAB⃝ on an Intel⃝ CoreTM i7 processor at 3.40 GHz with 16.0 GB of RAM. Computational experiments in Section 5.3 were performed using the university’s high performance computing facility.
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 5 Problem parameters. Parameter
Description
Value
Pmax Hmax Smax Tmax Dmax Bmax Cmax
Number of patients Number of surgeons Number of specialties Number of ORs Number of days Number of blocks Number of urgency categories
2094 107 10 20 5 3 3
5.1. Parameter tuning In this subsection we analyse the performance of SA, adaptive SA and hybridised SA–RVNS under different metaheuristic parameters. The performance of the above metaheuristics are analysed on a partially random dataset. The number of surgeons, specialties, ORs, days and blocks are kept constant. The surgeons and patients were taken from the elective surgery waiting list. The surgeon availabilities used were randomly generated. The same availability schedule was used for each run. The main scalar problem parameters used in the model can be found in Table 5. These values result in a problem with 667,488 binary decision variables, 600 non-negative continuous decision variables and 672,278 constraints. In this section we investigate all suitable combinations of metaheuristic parameters. We further investigate metaheuristic parameters based on promising behaviour. In each case we run the metaheuristics until variance of the best objective function value changes by less than 0.1 units after a further ten runs, or until 300 runs have been performed. 5.1.1. Simulated annealing A number of parameters are required by the SA metaheuristic. In particular, these are the initial temperature (T ), the number of iterations between temperature changes (imax ), the factor by which the temperature changes (α), and the minimum (stopping) temperature (Tmin ). In order to understand the full capabilities of SA on the model presented in Section 3 we must determine the most appropriate values for each of the metaheuristic parameters. After preliminary investigation, we consider all combinations of the SA parameter values given in Table 6. The five best performing SA parameter sets are shown in Table 7. The table also shows the average number of patients scheduled, lower and upper bounds of a 95% confidence interval and the average time required for each run. There is significant evidence that the average objective function obtained under T = 10, α = 09, imax = 150, Tmin = 0.001 is greater than that obtained under each of the other SA parameter combinations.
59
5.1.2. Adaptive simulated annealing The metaheuristic parameters required by adaptive SA are the initial temperature (T ), the maximum number of iterations and kp , where the probability of accepting a swap is given by exp 1cost × iter kp /T . Whilst it is also possible to change the probability function associated with selecting each type of neighbourhood swap, the current form performs well and will be used for the remainder of the investigation. To determine the best performance of adaptive SA on the model presented in Section 3 we consider suitable combinations of the adaptive SA parameters presented in Table 8. Table 9 shows the five best performing adaptive SA parameter combinations. The table shows the average objective function obtained, lower and upper bounds on a 95% confidence interval, and the average time required for a single run of the metaheuristic. 5.1.3. Reduced variable neighbourhood search The only metaheuristic parameter required by RVNS is the number of iterations performed. Here we consider between 400 and 6000 iterations. Table 10 shows the five best performing RVNS parameters along with the average number of patients scheduled, the lower and upper bounds of a 95% confidence interval, and the average time required for a single run. After 2200 iterations, an increase in the maximum number of iterations does not result in an increase in average objective function value. This is as RVNS tends to become ‘stuck’ at local optima and the objective function value plateaus. Note that RVNS does not perform better with only 2200 iterations, it merely performs as well as it does with a higher number of iterations. This is reflected in the overlapping confidence intervals. 5.1.4. Hybridised SA–RVNS The parameters required by hybridised SA–RVNS are the same as those required by adaptive SA. After initial analysis we consider the metaheuristic parameters shown in Table 11. Table 12 shows results of the SA–RVNS parameter tuning. The table contains the best performing hybridised SA–RVNS parameters, the average objective function, lower and upper bounds of a 95% confidence interval, and the average time required for a single run. At a 95% confidence level there is insufficient evidence of a difference in mean objective function obtained by each of these top five hybridised SA–RVNS parameter combinations. 5.2. Metaheuristic performance on random data Applying the results of the metaheuristic parameter tuning, any further analysis of the metaheuristics will be performed under the metaheuristic parameters shown in Table 13.
Table 6 SA parameters. Parameter
Description
Values
T
Initial temperature Temperature decrease factor Iterations between decreases Minimum temperature
1000, 100, 10 0.9, 0.7, 0.5 150, 100, 50 T × 10−1 , T × 10−2 , T × 10−3
α imax Tmin
Table 7 Results of SA parameter variation. T
α
imax
Tmin
Runs
Average
Lower
Upper
Time (s)
10 10 100 10 100
0.9 0.9 0.9 0.9 0.9
150 150 150 100 100
0.001 0.010 0.010 0.001 0.010
40 110 70 190 210
270.13 268.14 267.94 267.85 264.45
269.21 267.60 267.26 267.37 263.84
271.04 268.67 268.63 268.33 265.06
202.05 150.31 151.71 131.86 96.90
60
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
Table 8 Adaptive SA parameters. Parameter
Description
Values
T kp Imax
Initial temperature Exponent in acceptance probability function Maximum total number of iterations
1000, 100, 10 6, 4, 3, 2.5, 2 5000, 6000, 7000, 8000, 9000
Table 9 Results of adaptive SA parameter variation. T
kp
Imax
Runs
Average
Lower
Upper
Time (s)
100 100 100 100 100
2.5 3 2.5 6 3.5
9000 9000 8000 9000 9000
80 130 110 100 70
270.81 270.78 270.61 270.43 270.40
270.33 270.34 270.04 269.97 269.81
271.30 271.23 271.18 270.89 270.99
155.86 156.87 137.96 156.99 157.13
This reduces the variance of historical surgical durations and allows for more accurate estimation of surgical durations. The MINLP formulation provided in Section 3 does not allow a surgery to be scheduled if the 95th percentile of that surgery’s duration is greater than τBmax . Here, we replace constraint (11) with constraint (23) and include the constraints (24) and (25) that allows overtime when scheduling a single surgery in a full day block. if there is more than one patient scheduled in room t , on day d, block Bmax 0 otherwise,
1
Table 10 Results of RVNS parameter variation.
φtd =
Imax
Runs
Average
Lower
Upper
Time (s)
2200 2400 5000 3200 2300
60 50 40 270 130
265.58 265.18 264.85 264.80 264.70
264.51 263.91 263.40 264.18 263.78
266.65 266.45 266.30 265.41 265.62
174.12 186.07 388.37 247.76 178.17
∀t ∈ T , d ∈ D . Mt : an upper bound on the time that a surgery can take (in hours).
2 Utdb
Summary statistics are shown in Table 14, including the lower and upper bounds of a 95% confidence interval. From these statistics it can be seen that whilst at a 95% confidence level there is insufficient evidence that the mean objective found by Adaptive SA is different to that found by SA, there is very significant evidence that hybridised SA–RVNS outperforms each of the other metaheuristics. It is also worth noting the variance of the results. The objective value found by RVNS has a large variance as RVNS displays inconsistent performance. Much less variation is seen in the objective found using hybridised SA–RVNS which shows the metaheuristic to consistently perform well. Whilst there is insufficient evidence of a difference between the mean objective function value found by SA and Adaptive SA, Adaptive SA requires less time and has a lower variance. One limitation of the model presented in Section 3 is that no ENT surgeries can be scheduled. This is as, under the data shown in Table 4, it is not possible to schedule any ENT surgeries as the 95th percentile of the duration of a single surgery is greater than ten hours. Modifications required to address this issue are presented in Section 5.3.
2 Vtdb + Utdb
× exp 1.64485 ×
ln
2 Vtdb + Utdb
2 Utdb
≤ τb ,
∀t ∈ T , d ∈ D, b ∈ {1, 2} YtpdBmax , ∀t ∈ T , d ∈ D Mp φtd + 1 ≥
(23) (24)
p∈P
2 UtdB max
2 VtdBmax + UtdB max
2 VtdBmax + UtdB max × exp 1.64485 × ln 2
UtdBmax
≤ τBmax + (1 − φtd ) Mt ,
∀t ∈ T , d ∈ D.
(25)
Since any surgeons in ENT are now considered members of ENT1, ENT2, and ENT3 we need to include additional constraints to ensure patients from each subspecialty are scheduled appropriately. We include a binary decision variable χtdb and constraints (26) and (27). Note that ENTc is used to denote ENT subspecialty c.
5.3. Model performance
if room t on day d, block b is scheduled for an ENT subspecialty 0 otherwise, 1
In this section we investigate the performance of our model on the case study presented in Section 2. We use real data from June to September 2015. As we do not have staff availabilities, we use historical data to approximate actual availabilities. When conducting the metaheuristic parameter tuning no surgeries were scheduled for ENT as the 95th percentile on surgical duration was greater than ten hours. In this section, we split ENT into three subgroups separated by a patient’s urgency category.
χtdb =
∀t ∈ T , d ∈ D , b ∈ B χtdb ≥ Xtsdb , ∀t ∈ T , d ∈ D, b ∈ B
(26)
s∈ENT
Ytpdb βpc ≤ Xt (ENTc )db − χtdb + 1,
∀t ∈ T , p ∈ P , d ∈ D, b ∈ B, c ∈ C .
Table 11 Hybridised SA–RVNS parameters. Parameter
Description
Values
T kp Imax
Initial temperature Exponent in acceptance probability function The maximum total number of iterations
1000, 100, 10 1.5, 2, 2.5, 3 4000, 5000, 6000, 7000, 8000, 9000
(27)
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 12 Results of hybridised SA–RVNS parameter variation. T
kp
Imax
Runs
Average
Lower
Upper
Time (s)
1000 1000 10 1000 10
2.5 3.0 2.5 3.0 2.5
9000 7000 9000 8000 8000
90 60 90 130 80
272.19 272.08 271.92 271.88 271.85
271.76 271.62 271.54 271.52 271.40
272.61 272.55 272.30 272.25 272.30
153.68 119.67 158.54 137.07 140.39
Hence, the MINLP formulation used in this section consists of constraints (1)–(10), (12)–(19), and (23)–(27). Very little modification is required to calculate the new upper bound. Note that the metaheuristic solutions reported are feasible with respect to constraints (1)–(10), (12)–(19) and (23)–(27), i.e. all relevant constraints. The upper bound is calculated using constraints (1)–(8), (12)–(15), (22), (26) and (27). The upper bound is not guaranteed to be feasible with respect to the treat in turn policy. The upper bound also relies on a continuous relaxation of Ytpdb . Queensland University of Technology’s high performance computing facility was used to run CPLEX. Each instance was run across sixteen cores, with 128 GB of memory and a wall time of 168 h. Note that the time reported in Table 16 is the CPU time used. Table 15 shows the number of patients on the elective waiting list over the seven weeks under study. The table is split by urgency category and shows the number of patients on time and overdue at the time the waiting list was taken. The number of surgeries performed in each of these weeks is shown in Table 16. This includes both elective and emergency surgeries. The ratio of total number of surgeries waiting in each category is approximately 1:4:5 in categories one, two, and three respectively. It should be noted that the number of patients to be scheduled each week is now roughly 2900. In this section we consider 26 specialties and approximately 230 surgeons. As such, the problem size is significantly larger than that Section 5.1. The computational experiments performed in this section use the metaheuristic parameters given in Table 13. As in Section 5.1 we run each metaheuristic either until the variance of the best objective function value changes by less than 0.1 units after a further ten runs. We impose a maximum of 300 runs of each metaheuristic. Table 16 shows the results of computational experiments using the seven weeks of data available. The table shows the mean, variance, worst and best objective function value obtained by each metaheuristic after at most 300 runs. These results confirm the trends seen in Section 5.2. Two-sample t-tests were performed
61
on each combination of metaheuristics. For each week, the results display the same pattern. RVNS is the worst performing metaheuristic in terms of both lowest mean objective function value and highest objective function value variance. Whilst the mean objective function value obtained by SA is larger than that obtained by RVNS in six of the seven weeks (p ≪ 0.001), SA produces significantly worse solutions than Adaptive SA (p ≪ 0.001). When performing metaheuristic parameter tuning a maximum of 9000 iterations were considered for both Adaptive SA and hybridised SA–RVNS as the metaheuristics performed as well as or better than SA in considerably less time. In future it may be worth considering a larger total number of iterations. In this case we consider computational time as a limited resource. The best performing metaheuristic for each week was hybridised SA–RVNS. This heuristic provides solutions with a low variance and highest mean objective function value. Hybridised SA–RVNS required the least (or second least) amount of computational time and produced significantly better results than the next best metaheuristic, Adaptive SA (p ≪ 0.001). In several instances, the best objective function value obtained using metaheuristics was within ten units of the best integer solution obtained by CPLEX on the relaxed problem. The range of objective function values obtained through the hybridised SA–RVNS metaheuristic is very promising, especially considering the size of the real life instances addressed. When comparing the results produced by the metaheuristics using the model presented in Section 3 to historical data we find that the metaheuristics are able to schedule around 50–100 additional surgeries each week. This result should be considered with caution as surgeon availabilities were merely estimates based on historical data (due to the lack of available data); however these estimates are as accurate as possible. There is also the possibility that downstream wards may be unable to satisfy the increased demand. 6. Conclusions and future work In this paper we addressed the combined MSSP SCAP in Australian public hospitals using a MINLP formulation and tackled all three levels of OT planning/scheduling problems simultaneously. Whilst commercial solvers are unable to cope with the size of the case study considered the metaheuristics and hybridised metaheuristics produce good feasible solutions in reasonable amounts of computational time. In particular
Table 13 Metaheuristic parameters. Parameter Description
SA
RVNS Adaptive SA
SA–RVNS
T
10 0.9 150
NA NA NA
1000 NA NA
α imax kp Tmin Imax
Starting temperature Temperature decrease Number of iterations between temperature decrease Acceptance probability function exponent The stopping temperature The maximum total number of iterations
100 NA NA
NA NA 2.5 0.001 NA NA NA 2200 9000
2.5 NA 9000
Table 14 Summary statistics of computational experiments. Method
Runs
Average
Variance
Lower
Upper
Time (s)
RVNS SA Adaptive SA SA–RVNS Upper bound
60 40 80 90 –
265.58 270.13 270.81 272.19 305.00
17.91 8.63 4.64 4.24 –
264.51 269.21 270.33 271.76 –
266.65 271.04 271.30 272.61 –
174.12 202.05 155.86 153.68 50.01
62
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64 Table 15 Case study weekly waiting list data by urgency category. Week
1 2 3 4 5 6 7
1
2
3
Total
On time
Overdue
On time
Overdue
On time
Overdue
337 323 313 292 302 280 294
0 0 0 0 0 22 1
1127 1118 1117 1131 1102 1077 1130
74 70 63 69 75 104 76
1326 1378 1382 1372 1367 1364 1364
16 11 13 15 14 13 16
2880 2900 2888 2879 2860 2860 2881
Table 16 Case study weekly results—metaheuristic performance. Week
Method
Mean
Variance
Worst
Best
Runs
Mean time (s)
1
Actual RVNS SA Adaptive SA–RVNS Upper bound
297 308.78 320.64 328.06 330.77 383.01
– 394.34 14.50 26.63 20.47 –
– 148 309 315 316 –
– 328 331 343 339 –
– 300 160 150 110 –
– 198.74 226.67 182.16 183.27 685 938.00
2
Actual RVNS SA Adaptive SA–RVNS Upper bound
303 329.84 332.34 344.76 349.53 391.68
– 68.46 28.72 30.42 22.79 –
– 301 321 327 335 –
– 349 349 357 359 –
– 300 90 170 270 –
– 222.44 270.48 200.24 192.51 4 582 191.00
3
Actual RVNS SA Adaptive SA–RVNS Upper bound
266 333.70 335.05 346.31 351.84 391.53
– 48.44 22.62 20.40 15.95 –
– 308 323 331 333 –
– 353 347 357 362 –
– 220 190 130 170 –
– 219.48 238.59 194.85 190.24 5 008 999.00
4
Actual RVNS SA Adaptive SA–RVNS Upper bound
306 330.09 334.86 343.13 348.57 393.40
– 50.47 27.35 19.54 20.17 –
– 305 318 333 336 –
– 345 348 354 359 –
– 300 110 200 220 –
– 217.42 243.05 193.02 192.12 4 386 662.00
292 325.81 328.26 340.52 345.36 395.31
– 33.89 23.97 32.56 30.40 –
– 312 314 322 326 –
– 341 340 356 357 –
–
5
Actual RVNS SA Adaptive SA–RVNS Upper bound
70 120 300 190 –
– 217.20 235.06 244.90 224.78 8 772 716.00
6
Actual RVNS SA Adaptive SA–RVNS Upper bound
232 312.22 314.91 322.43 325.58 392.54
– 150.22 11.05 17.41 19.09 –
– 147 308 314 317 –
– 331 326 331 336 –
– 300 230 70 60 –
– 199.65 227.15 217.10 173.45 3 364 606.00
7
Actual RVNS SA Adaptive SA–RVNS Upper bound
256 304.01 309.56 316.89 319.36 387.82
– 227.02 15.89 25.14 28.35 –
– 137 301 305 304 –
– 321 320 330 331 –
– 300 90 140 160 –
– 215.10 236.76 202.39 210.66 6 935 316.00
hybridised SA–RVNS, which has not previously been implemented on large MINLP problems or OT planning/scheduling problems, outperformed our chosen benchmark metaheuristics and has potential to be applied to other complex hospital scheduling problems on a large scale. When modelling the combined MSSP SCAP we considered the Australian government’s current waiting list management policies: the treat in turn policy and patient wait targets. Future work includes the consideration of other possible policies that may better manage the elective surgery waiting lists. Possible extensions to the model include incorporating length of stay, emergency arrivals, and capacity of downstream wards.
Acknowledgements This research was funded by the Australian Research Council (ARC) Linkage Grant LP 140100394 and supported by the Princess Alexandra Hospital, Brisbane, Australia. We would like to thank Dr. Andy Wong for his assistance in obtaining data for our case study. We would also like to thank Dr. Michael Sinnot, Dr. David Cook, Sean Birgan and other staff at the PAH for their considerable feedback and their time. Appendix (See Fig. 1.)
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
63
Fig. 1. A typical MSS.
References [1] C. Van Riet, E. Demeulemeester, Trade-offs in operating room planning for electives and emergencies: A review, Oper. Res. Health Care (2015). [2] Y.B. Ferrand, M.J. Magazine, U.S. Rao, Managing operating room efficiency and responsiveness for emergency and elective surgeries — a literature survey, IIE Trans. Healthc. Syst. Eng. 4 (1) (2014) 49–64. [3] F. Guerriero, R. Guido, Operational research in the management of the operating theatre: a survey, Health Care Manage. Sci. 14 (2011) 89–114. [4] I. Adan, J. Bekkers, N. Dellaert, J. Vissers, X. Yu, Patient mix optimisation and stochastic resource requirements: a case study in cardiothoracic surgery planning, Health Care Manage. Sci. 12 (2) (2009) 129–141. [5] J.M. Van Oostrum, M.V. Van Houdenhoven, J.L. Hurink, E.W. Hans, G. Wullink, G. Kazemier, A master surgical scheduling approach for cyclic scheduling in operating room departments, OR Spectrum 30 (2) (2008) 355–374. [6] H. Fei, N. Meskens, C. Chu, A planning and scheduling problem for an operating theatre using an open scheduling strategy, Comput. Ind. Eng. (2010) 221–230. [7] P. Santibáñez, M. Begen, D. Atkins, Surgical block scheduling in a system of hospitals: an application to resource and wait list management in a British Columbia health authority, Health Care Manage. Sci. 10 (3) (2007) 269–282. [8] J.M. Van Oostrum, E. Bredenhoff, E.W. Hans, Suitability and managerial implications of a master surgical scheduling approach, Ann. Oper. Res. 178 (1) (2010) 91–104. [9] A. Agnetis, A. Coppi, M. Corsini, G. Dellino, C. Meloni, M. Pranzo, Long term evaluation of operating theater planning policies, Oper. Res. Health Care 1 (2012) 95–104. [10] J. Beliën, E. Demeulemeester, B. Cardoen, A decision support system for cyclic master surgery scheduling with multiple objectives, J. Sched. 12 (2) (2009) 147–161. [11] S. Chaabane, N. Meskens, A. Guinet, M. Laurent, Comparison of two methods of operating theatre planning: Application in Belgian hospital, J. Syst. Sci. Syst. Eng. 17 (2) (2008) 171–186. [12] A. Testi, E. Tànfani, Tactical and operational decisions for operating room planning: Efficiency and welfare implications, Health Care Manage. Sci. 12 (4) (2009) 363–373. [13] Y. Wang, J. Tang, R.Y.K. Fung, A column-generation-based heuristic algorithm for solving operating theater planning problem under stochastic demand and surgery cancellation risk, Int. J. Prod. Econ. 158 (2014) 28–36. [14] H. Fei, C. Chu, N. Meskens, Solving a tactical operating room planning problem by a column-generation-based heuristic procedure with four criteria, Ann. Oper. Res. 166 (1) (2009) 91–108. [15] H. Fei, C. Chu, N. Meskens, A. Artiba, Solving surgical cases assignment problem by a branch-and-price approach, Int. J. Prod. Econ. 112 (1) (2008) 96–108. [16] H. Saadouli, B. Jerbi, A. Dammak, L. Masmoudi, A. Bouaziz, A stochastic optimization and simulation approach for scheduling operating rooms and recovery beds in an orthopedic surgery department, Comput. Ind. Eng. 80 (2015) 72–79.
[17] K. Steins, F. Persson, M. Holmer, Increasing utilization in hospital operating department using simulation modeling, Simulation 86 (8–9) (2010) 463–480. [18] B. Roland, C. Di-Martinelly, F. Riane, Y. Pochet, Scheduling an operating theatre under human resource constraints, Comput. Ind. Eng. 58 (2) (2010) 212–220. [19] A. Agnetis, A. Coppi, M. Corsini, G. Dellino, C. Meloni, M. Pranzo, A decomposition approach for the combined master surgical schedule and surgical case assignment problems, Health Care Manage. Sci. 17 (2014) 49–59. [20] R. Aringhieri, P. Landa, P. Soriano, E. Tànfani, A. Testi, A two level metaheuristic for the operating room scheduling and assignment problem, Comput. Oper. Res. 54 (2015) 21–34. [21] J.E. Everett, A decision support simulation model for the management of an elective surgery waiting system, Health Care Manage. Sci. 5 (2002) 89–95. [22] D. Min, Y. Yih, An elective surgery scheduling problem considering patient priority, Comput. Oper. Res. 37 (2010) 1091–1099. [23] D. Min, Y. Yih, Managing a patient waiting list with time-dependent priority and adverse events, RAIRO-Oper. Res. 48 (2014) 53–74. [24] I. Marques, M.E. Captivo, M. Vaz Pato, A bicriteria heuristic for an elective surgery scheduling problem, Health Care Manage. Sci. (2014) 1–16. [25] Queensland Health, Treating in turn — calculating the treat in turn percent, in: Department of Health (ed.), 2013. [26] E. Kayış, T. Khaniyev, J. Suermondt, K. Sylvester, A robust estimation model for surgery durations with temporal, operational, and surgery team effects, Health Care Manage. Sci. 18 (2015) 222–233. [27] F. Dexter, J. Ledolter, V. Tiwari, R.H. Epstein, Value of a scheduled duration quantified in terms of equivalent numbers of historical cases, Anesth. Analg. 117 (1) (2013) 205–210. [28] D. Strum, J. May, L. Vargas, Surgical procedure times are well modeled by the lognormal distribution, Anesth. Analg. 86 (1998) 47S. [29] K. Stuart, E. Kozan, Online scheduling in the operating theatre, in: Industrial Engineering and Management Systems Conference, 2009 Kitakyushu, Japan, pp. 802–807. [30] M. Heng, J. Wright, Dedicated operating room for emergency surgery improves access and efficiency, Can. J. Surg. 56 (3) (2013) 167–174. [31] K. Stuart, E. Kozan, Reactive scheduling model for the operating theatre, Flex. Serv. Manuf. J. 24 (4) (2012) 400–421. [32] E. Hans, G. Wullink, M. Van Houdenhoven, G. Kazemier, Robust surgery loading, European J. Oper. Res. 185 (3) (2008) 1038–1050. [33] A. Riise, E. Burke, Local search for the surgery admission planning problem, J. Heuristics 17 (2011) 389–414. [34] P. Landa, R. Aringhieri, P. Soriano, E. Tànfani, A. Testi, A hybrid optimization algorithm for surgeries scheduling, Oper. Res. Health Care 8 (2016) 103–114. [35] C. Rizk, J.-P. Arnaout, ACO for the surgical cases assignment problem, J. Med. Syst. 36 (2012) 1891–1899. [36] M. Birattari, Tuning Metaheuristics: A Machine Learning Perspective, Springer, Berlin, Heidelberg, 2009. [37] S. Ceschia, A. Schaerf, Local search and lower bounds for the patient admission scheduling problem, Comput. Oper. Res. 38 (2011) 1452–1463.
64
B. Spratt, E. Kozan / Operations Research for Health Care 10 (2016) 49–64
[38] S. Vlah, Z. Lukač, J. Pacheco, Use of VNS heuristics for scheduling of patients in hospital, J. Oper. Res. Soc. 62 (2011) 1227–1238. [39] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680. [40] S. Ceschia, A. Schaerf, Dynamic patient admission scheduling with operating room constraints, flexible horizons, and patient delays, J. Sched. (2014) 1–13. [41] A. Fügener, E.W. Hans, R. Kolisch, N. Kortbeek, P.T. Vanberkel, Master surgery scheduling with consideration of multiple downstream units, European J. Oper. Res. 239 (2014) 227–236. [42] N. Mladenović, P. Hansen, Variable neighborhood search, Comput. Oper. Res. 24 (1997) 1097–1100.
[43] P. Hansen, N. Mladenović, J. Brimberg, J.M. Pérez, Variable neighborhood search, in: M. Gendreau, J.-Y. Potvin (Eds.), Handbook of Metaheuristics, Springer, US, 2010. [44] E.K. Burke, J. Li, R. Qu, A hybrid model of integer programming and variable neighbourhood search for highly-constrained nurse rostering problems, European J. Oper. Res. 203 (2010) 484–493. [45] R. Aringhieri, P. Landa, E. Tànfani, Assigning surgery cases to operating rooms: A VNS approach for leveling ward beds occupancies, Electron. Notes Discrete Math. 47 (2015) 173–180. [46] A. Alonso-Ayuso, L. Escudero, F.J. Martín-Campo, N. Mladenović, A VNS metaheuristic for solving the aircraft conflict detection and resolution problem by performing turn changes, J. Global Optim. 63 (2015) 583–596.