{'fraction_non_alphanumeric': 0.026521378785614325, 'fraction_numerical': 0.009184657817149224, 'mean_word_length': 4.917510853835021, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Surgical Scheduling via Optimization and Machine Learning with Long-Tailed Data

Using data from cardiovascular surgery patients with long and highly variable post-surgical lengths of stay (LOS), we develop a modeling framework to reduce recovery unit congestion. We estimate the LOS and its probability distribution using machine learning models, schedule procedures on a rolling basis using a variety of optimization models, and estimate performance with simulation. The machine learning models achieved only modest LOS prediction accuracy, despite access to a very rich set of patient characteristics. Compared to the current paper-based system used in the hospital, most optimization models failed to reduce congestion without increasing wait times for surgery. A conservative stochastic optimization with sufficient sampling to capture the long tail of the LOS distribution outperformed the current manual process and other stochastic and robust optimization approaches. These results highlight the perils of using oversimplified distributional models of LOS for scheduling procedures and the importance of using optimization methods well-suited to dealing with long-tailed behavior.

spite the use of a wide range of machine learning models and a rich set of patient characteristics.
• Optimization models that rely on machine learning predictions of LOS without accounting for extended LOS did not improve scheduling performance (recovery unit congestion and wait times of patients) relative to current paper-based systems in use. • We show a data-driven conservative stochastic optimization approach that accounts for stochasticity in extended LOS can achieve scheduling performance improvements, outperforming other stochastic and robust optimization approaches. • We apply and evaluate our methodology in the context of a pediatric academic medical center using real medical and operational data.

Introduction
For hospital-based surgical care, the capacity of the intensive care unit (ICU) is often a crucial downstream bottleneck. ICU bed shortages are associated with adverse patient outcomes, lost revenue from cancelled procedures, and a variety of detrimental spillover issues for numerous parts of the hospital. For many types of surgery, patients are sent to an intensive care unit (ICU) to recover until they meet the criteria to be transferred to a step-down unit. The capacity of the ICU, especially the number of beds available, is often a crucial bottleneck of surgical planning. When the ICU is at capacity, staff may transfer patients prematurely to the step-down unit or surgical procedures may be cancelled at the last-minute with adverse impact on patient and family experience, hospital reputation, finances and staff morale [3].
While most patients that require an ICU bed require it urgently, elective surgical procedures are scheduled as far as a year in advance. Our primary goal in this work was to develop a scheduling model that would reduce post-surgical bed congestion in practice in the presence of difficult-to-predict, long-tailed LOS data. Our secondary goals were to address the challenges associated with optimization in the presence of long-tailed empirical data. Numerous works have examined how to optimize surgical scheduling in order to reduce recoverybed congestion. Despite the importance of post-surgical LOS, studies commonly discard empirical LOS data after they have been used to fit the parameters of a distribution and use only synthetic data to evaluate model performance. We illustrate limitations associated with using synthetic, rather than empirical, data.

Setting
This work was performed at the heart center of a high surgical volume pediatric academic medical center (PAMC) in the United States. The heart center uses three operating rooms, 26 cardiovascular ICU (CVICU) beds shared by elective and urgent surgical patients, and an acute care unit. In recent years, growing demand for elective surgical services and fixed CVICU capacity has resulted in numerous surgical cancellations. From September 2019 to May 2020, 84% of cardiovascular surgery cancellations happened on the 40% of days with 23 or more patients in total in the CVICU. In particular, 35% of cancellations happen on the 11% of days with more than 10 elective surgical patients in the CVICU (See Figure 1). This suggests that a significant fraction of cancellations may be prevented if optimized scheduling smooths elective surgical patient CVICU occupancy.

Overview
Numerous theoretical works using optimization to schedule elective surgical cases have demonstrated substantial reductions to ICU congestion in numerical experiments [17,19,31,23]. However, relatively little work has been implemented or has demonstrated measurable improvements in practice. Most existing positive results rely on the use of synthetic data, e.g., a lognormal distribution fit to approximate empirical data from the institution studied, to represent procedure duration and post-operative length of stay. Such common practice, combined with the lack of real-life implementations of these algorithms, leaves many real world challenges of surgical scheduling unaddressed and unidenti- We examine reducing ICU congestion by optimizing the scheduling of complex pediatric cardiovascular surgical procedures, using real data from the hospital for patient flow and post-operation ICU LOS. Our approach combines machine learning for LOS forecasting, optimization for surgical scheduling, and simulation for performance evaluation.
To motivate the problem, we first establish the theoretical upper bound of system performance by running a deterministic optimization formulation with known a priori information on post-operation LOS. The results demonstrate significant reductions in CVICU congestion and patient wait times, compared to the institution's current manual and primarily paper-based process. We then develop predictive models using machine learning for LOS based on patient data available at the time the procedure is scheduled. When perfect information of LOS is replaced with point-predictions using machine learning, however, deterministic optimization is found to be not sufficient to reduce congestion without increasing wait times. This is despite access to a very rich set of patient characteristics for LOS prediction.
The key bottleneck is identified to be the lack of accurate predictions of LOS, especially for a small group of patients with very long LOS. This has a disproportionately large impact on optimization performance due to the long tail of the distribution of LOS. We observe that patients with similar characteristics may vary drastically in realized LOS, and patient populations display temporal non-stationarity, both suggesting inherent un-predictability of LOS at the time of scheduling. This inspired the development of stochastic and robust optimization algorithms that incorporate LOS probability distributions instead of point estimates to address this bottleneck.
We develop a novel data-driven modeling framework for scheduling cardiovascular surgery under uncertain LOS. The framework combines machine-learning predicted LOS probability distributions, rolling information update and optimization methods. Under this framework, three scheduling algorithms are formulated: two based on stochastic optimization, namely Standard-RSO and Conservative-RSO, and one robust optimization formulation, RRO. The optimal schedules given by each algorithm are evaluated through simulation using historical patient arrivals and LOS.
Machine learning models are used to obtain both point-estimates of LOS and its probability distribution characterized by predictive errors. Standard-RSO, which considers both under-and over-estimation errors, fails to reduce ICU congestion versus the status quo. In contrast, both Conservative-RSO and RRO, with carefully tuned parameters, managed to reduce ICU congestion without excessively increasing average patient wait times. We attribute the promising performances of the latter two algorithms to their targeted design focusing on addressing common under-estimation for prolonged LOS by machine learning predictions. Moreover, Conservative-RSO achieves better average and worstcase performance compared to RRO in reducing ICU congestion, revealing that algorithmic complexity does not always equal better performance in practical settings.
The main contributions of our work are two fold. First, the negative results associated with deterministic optimization as well as Standard-RSO offer important lessons that are generally applicable to others seeking to optimize surgical scheduling for complex patients. We highlight the importance of evaluating model performance using empirical data or simulation that fully captures the long-tailed nature of patient LOS. In addition, as hospital operations and medical practices are constantly changing, non-stationarity in LOS should be taken into account when developing simulating and evaluating data-driven models. Our negative results, revealed through simulation with empirical data, provide clear evidence against extrapolating good performance on synthetic data using standard distributional assumptions (e.g. lognormal distribution) into real world scenarios. These findings may also help explain the dearth of studies reporting successful implementation and measured improvements of similar approaches.
Our second contribution is to propose a promising algorithm, Conservative-RSO, that is specifically designed to address the challenge of long-tailed ICU LOS in the surgical scheduling problem. Compared to Standard-RSO, Conservative-RSO demonstrates that good performance can be achieved by incorporating careful design choices without increasing computational complexity; the comparison with RRO further suggests that computationally complicated models do not always translate into better performance in practice. We believe the promising performance of Conservative-RSO helps identify a direction of future research on datadriven optimization methods to address one of the most challenging obstacles in hospital operations.
The remainder of the paper is structured as follows. Section 2 provides a review of related literature on surgical scheduling optimization and LOS predictions. Section 3 presents our modeling framework under deterministic LOS estimations, introducing formulations for both offline optimization and rolling optimization. These deterministic models are used to establish performance upper-bounds given accurate predictions of patient LOS. Section 4 presents the development of machine learning models for LOS predictions and discusses the poor performance and challenges of using deterministic optimization with machine learning prediction. Section 5 presents our modeling framework using predicted LOS distributions, and introduces the three algorithms, Standard-RSO, Conservative-RSO and RRO. The performances of these algorithms under numerical experiments are presented in Section 6. Section 7 discusses insights from the design and performance of algorithms, as well as their limitations and potential extensions. Finally, Section 8 provides concluding remarks on implications of the work, best practices in applying schedule optimization in practice, and directions for future research.
Gurobi 9.0.3 [10] and Python were used to implement and solve all mixed integer optimization algorithms developed in the paper. All optimization problems were solved using the compute facility of Stanford Research Computing Center and Stanford University. On this shared facility, each optimization problem was solved with two CPUs and with either 8GB or 16GB of memory for each CPU.

Surgery Scheduling
Extensive research has been carried out on scheduling of patients to improve operating room (OR) performance.
Most of the literature on elective surgery scheduling focuses on OR room capacities without considering the capacity of the subsequent recovery units. The variability of surgical procedure durations is usually considered as the primary source of stochasticity and as the primary challenge to surgical scheduling. We refer to [32,22] for a comprehensive review and focus on research that studies the surgical scheduling problem under the constraint of limited downstream capacity.
Deterministic formulations of surgical schedule optimization under limited downstream capacity have been examined in [11,9,21,7]. [11] considered optimizing surgical scheduling with limited capacity at the Post Anesthesia Care Unit (PACU) to minimize the number of nurses needed for the PACU. The case duration and the recovery times in the PACU are treated as deterministic. Similarly, [21] considered the surgical scheduling problem with deterministic recovery times in the PACU and the ICU for all patients under the assumption of perfect information. [9] also investigates surgery scheduling with downstream PACU capacity constraints. It proposes a two step optimization procedure where patients are first assigned to ORs and dates. In the second step surgery times are assigned to patients. More recently, [7] developed a combined machine learning and optimization approach to reduce congestion in the PACU, where the recovery time in the PACU is estimated using machine learning predictions. [7] reports good results based on simulations using empirical data, but does not report the results of implementation. Most of the existing deterministic formulations focus on the PACU capacity as the key downstream bottleneck with limited discussions on the ICU capacity. Compared to the recovery time at the PACU (which is typically in the order of hours), ICU LOS may be multiple days or months and involves greater variance than single-day PACU stays. The substantial difference in variability between CVICU long-tailed LOS and PACU recovery time are why the approach in [7] is not applicable in the present setting.
Given the multi-period nature and significant uncertainties involved in ICU capacity planning, a number of mathematical formulations have been proposed in recent years. [17] provided a well-known benchmark instance of surgery scheduling with ICU capacity constraints using a stochastic mixed integer programming model. The work uses sample average approximation and assumes the LOS in the ICU to be random with known distributions (arbitrary with finite support).
[31] proposes a two level time horizon for surgery scheduling. In the first level, patients are selected from a waitlist to be scheduled with a timeframe (e.g. a week) using approximate dynamic programming. In the second level, selected patients are scheduled using a sample average approximation method similar to [17]. Other stochastic programming approaches built on the work in [17]; we refer to [23] for a review of these formulations.
Besides stochastic optimization models, [19,23] develop optimization methods for robust performance of surgery scheduling under worst-case realizations of LOS or LOS distributions. Most relevant to our work is [19]. The authors formulated a two-stage robust optimization approach to reduce congestion in downstream capacities, and developed solution techniques which we adapt and apply to our setting.
All of the above-mentioned papers on stochastic or robust optimization evaluated model performance using synthetic data based on strong distributional assumptions, and none has been implemented in practice. In contrast, we use real LOS data in evaluating our model performance.
In addition, [8] considered a surgery scheduling problem focusing on minimizing downstream costs including overcapacity costs at the ICU. Models for this problem proposed in [8] and others such as [1] are concerned with the tactical problem of allocating OR block of times to surgical specialties to optimize patient flow into the downstream units. In contrast, our work and others' mentioned above focus on the operational problem of assigning individual patients to surgical time blocks.
To our knowledge, we are the first to study a combined machine learning and optimization approach for schedule optimization to reduce congestion at the CVICU: an environment that includes non-stationary, long-tail LOS behavior. We are also the first to propose a stochastic formulation for surgical scheduling that is specifically designed to address significant variability in the tail of the LOS distribution. We use real world postoperation LOS data for LOS prediction, model parameter tuning and performance evaluation, reflecting the difficulties of working in a setting with non-stationary operations and patient volumes and long-tailed LOS distributions.

Post-surgery LOS prediction
In developing our machine learning models for LOS prediction, we refer to a separate line of literature focusing on predictions of post-surgery (in particular cardiac surgeries) LOS in the intensive care units.
Most work suggests that post-surgery LOS prediction at admission time -either using predictive modeling or by expert opinion -is challenging especially in the case of prolonged ICU LOS.
One common way to predict prolonged ICU LOS is through binary classification, such as in [28,6]. However, binary classification does not provide the level of granularity required for optimizing surgical scheduling in our context.
For regression-based models, [29] builds and evaluates multivariate regression models using data from 246 hospitals for heart failure patients. The model achieves a modest R 2 value of 4.8%, where only 1.2% of variation is explained by patient characteristics. Similarly, [4] shows through univariate regression that only 12% of the variation could be explained by patient characteristics and general hospital characteristics in aggregate for patients with a primary diagnosis of acute myocardial infarction. More complicated LOS regression models using machine-learning for cardiac surgery patients are artificial neural networks developed in [25,14] and adaptive neuro-fuzzy systems explored in [16]. Although machine-learning based models generally result in a higher R 2 value, predictive accuracy for patients with prolonged LOS remains low. For instance, the neural network in [25] achieves an overall accuracy of over 60% but is unable to predict any LOS above 15 days, and the model in [14] with R 2 = 0.41 underestimates the LOS for almost all patients with actual LOS of longer than 100 hours.
More generally, [27,18] demonstrate that prolonged LOS predictions are challenging even for experienced physicians. In addition, [13] suggests that information gathered at admission did not have a significant impact on the identification of patients with prolonged ICU LOS. The variables that had the greatest impact on prolonged ICU LOS were those measured on day 5 of ICU stays. Although LOS prediction using patient features collected during and after surgery tends to achieve a higher level of accuracy [13,24,30], most of these data are not available when a procedures is being scheduled weeks to months in advance.
In our attempt to develop a predictive model for post-operation LOS, we identified the same difficulty as observed by the others in dealing with the long-tailed, non-stationary distribution of LOS for complex cardiovascular surgeries. This challenge motivated our design of a stochastic optimization formulation that specifically addresses the unpredictability and the consistent underestimation of prolonged LOS by predictive models. We discuss the formulation in detail in Section 5.

Deterministic Formulations for Schedule Optimization
We first study the problem of surgical scheduling where the LOS in the ICU are treated as deterministic. We de-velop two optimization models -offline and rolling deterministic optimization -given operational constraints at the PAMC and some a priori information on patient LOS.

Offline Deterministic Optimization
We formulate the offline surgery scheduling problem as a mixed integer program (MIP) with the objective of smoothing out CVICU elective census overtime to reduce cancellations without creating excessive delays.
Define set D = {1, 2, . . . , N d } as the set of available dates in the time period of interest, and define set P = {1, 2, . . . , N p } as the set of all elective surgical patients that are to be scheduled within the time period. We define binary decision variables, where x d,p = 1 if patient p is scheduled for surgery on day d and otherwise x d,p = 0.
For the offline problem, it is assumed that both patient arrivals, P , and the LOS of each patient, l p , are known in advance. Given x and {l p : p ∈ P }, the ICU overflow variable, u d for all d ∈ D, counts the number of elective patients on day d that exceed ICU capacity, c.
The offline deterministic optimization problem is formulated as below.
We use the notation (·) + = max(·, 0). In the offline problem, the first term of the objective function (1a) represents the total wait time for all patients. Here, d min p denotes the earliest date by which patient p will be ready for his surgery. The second term, β d∈D f (u d ), is the total cost associated with ICU overflow, weighted by constant β. We set the function f (·) to be a convex, piece-wise linear function of u d such that it approximates the quadratic cost function, The purpose of the convex cost function is to impose higher penalty for greater overflow in order to smooth out and minimize large peaks in ICU occupancy. This also reflects the highly non-linear increase in the possibility of cancellation as elective ICU occupancy increases, as demonstrated in Figure 1. We provide the exact MIP formulation of f (·) in Appendix A.
Constraint (1b) enforces the assignment of every patient to exactly one date of surgery. The second constraint (1c) calculates if a patient needs an ICU bed on a given day based on her assignment and LOS. The third constraint (1d) calculates the ICU overflow, u d , on each day given ICU bed capacity, c.
For the last constraint (1e), we use Q op to represent any remaining institution-specific operational constraints. Scheduling for complex procedures is commonly constrained by single surgeon-patient matches, surgeon and OR room availability, which may differ depending on the context. For our formulation, Q op is tailored to the context of the PAMC as explained below.
Patient availability. When scheduling patients for surgery, there is often a clinically determined upper bound on patient wait time and also a corresponding lower bound that accounts for surgical delays due to patient availability, travel arrangement and necessary insurance or health checks. We incorporate this constraint by imposing windows of availability for individual patients as hard constraints in our formulation, similar to [9]. In practice, the window of availability of patients can be estimated case by case upon arrival and hence used as inputs for optimization.
For the purpose of simulation, since a patient's actual window of availability at the time of arrival is not recorded in historical census and surgery data from the PAMC, we use the following heuristics to obtain an estimation. The earliest available date of patient p for surgery, denoted as d min p ∈ D, is set to be date that is half of the lead time prior to her originally scheduled surgery date. The latest available date, denoted as d max p ∈ D, is set to be 90 days after the actual surgery date. This heuristic ensures that the set of feasible scheduling solutions is not too restricted, and any resultant increase in wait time will be penalized by the first term in the objective function.
Surgeon and OR Room availability. Each procedure requires one OR room and takes a pre-specified amount of time for the patient's assigned surgeon to complete. While many past work focused their attention on the stochastic nature of surgical durations (see, for example, [19]), OR room availability is not a bottleneck resource in our setting. We therefore treat surgical du-rations as deterministic and uses the actual length of procedure in optimization 1 .
Specific to the PAMC, a maximum of 45 hours of surgery in OR is available on every weekday, where each surgeon performs operations for no longer than 15 hours every day 2 . During implementation, we also incorporate surgeon-specific day-offs in the formulation. For example, one of the surgeons does not perform surgeries on Mondays.
Complex surgeries. On top of the regular capacity constraints above, special arrangements are often required for complex, full-day procedures. At the PAMC, Pulmonary Artery Reconstruction (PAR) surgery is a type of highly complex procedure that requires special treatment in the algorithm. The surgeons only perform PAR surgeries on certain weekdays ('PAR days') and not the others. If a PAR surgery is scheduled for a surgeon on any day, no other surgery can be performed by the same surgeon on the same day. In addition, the surgeons need to take at least a one-day break in between PAR surgeries, i.e., PAR surgeries cannot be scheduled on two consecutive days for the same surgeon.
Full mathematical formulations of Q op for our context and implementation details are provided in Appendix A.
Here, the weight for ICU overflow, β, and ICU bed capacity, c are parameters to be tuned. For the value of β, we select a value among β ∈ {5, 10, 25, 50, 100, 200} so that the model achieves the the most reduction on ICU overflow without excessively lengthening the average and median patient wait times compared to the original surgical schedule. Meanwhile, we determine the value of c based on the minimum achievable ICU occupancy upper bound given true historical on LOS, i.e., s.t. Constraints (1b), (1c), (1e), and (1f) Solving (2) yields c = 8 at the PAMC, which is used for all numerical experiments for the rest of the paper.
1 Alternatively, we also estimated each patients' surgical duration using the average surgical duration for each procedure type in the empirical data. This approach produces similar results to that using the actual surgical durations. 2 Surgeries that took longer than 15 hours in reality were treated as 15 hours for feasibility.

Rolling Deterministic Optimization (RDO)
In practice, the arrival and LOS of incoming patients are unknown and performing offline optimization over a one-year horizon is not practical. We thus develop a rolling-horizon alternative to the offline formulation, where patients are scheduled in batches with estimated LOS at the time of scheduling. Meanwhile, we dynamically observe the realization of LOS for patients who have undergone surgery in the previous period and update LOS estimates for these patients.
We define a sequence of scheduling days, s b ∈ D, for b ∈ {1, 2, . . . , B} in chronological order. On each scheduling day s b , scheduling is performed for all patients who arrived between s b−1 and s b − 1. We use P b to denote the batch of patients who are scheduled on day s b . These patients may be scheduled for surgery any time during the scheduling horizon, where N d is sufficiently large to cover the scheduling horizon for all batches. 3 Note that scheduling decisions for patients in P b are also influenced by past patients, P past b = b−1 k=1 P k who are still in the system. This includes all patients who arrived before day s b−1 who are currently staying in the ICU on day s b , as well as those who have their surgery scheduled for days within the current scheduling horizon. Surgical and ICU capacities are adjusted accordingly to account for these past scheduling decisions.
To introduce the mathematical formulation, we use the same set of notations as the offline problem. Note that the binary decision variables x d,p , y d,p are now defined for all d ∈ D and p ∈ P b ∪ P past b . Since scheduling decisions for patients in P past b are already made, we let constantsx = {x d,p : d ∈ D, p ∈ P past b } denote previous scheduling decisions. This constraint implicitly updates surgical and ICU capacities available for incoming patients by accounting for patients who have been previously scheduled. The deterministic batch optimization problem (BOP) on scheduling day s b is formulated as below.
3 Note that we still enforce the constraint that each patient is scheduled no later than d max The objective function (3a) mirrors that of the offline formulation in Equation (1a). Constraint (3b) requires that patients scheduled previously are not rescheduled for a different date. Constraint (3c) requires that every patient in P b is scheduled for one surgery day no earlier than s b . Constraints (3d)-(3g) mirror constraints (1c)-(1f) in the offline problem. Note that constraints (3e) and (3f) are influenced by both patients in P b and those in P past b . This is because the latter may also undergo surgery and/or need ICU beds after the scheduling date s b , thus competing for ICU resources and other resources in Q op (e.g. surgeon and OR room availability) with patients in P b . Full mathematical formulation with our context-specific Q op and implementation details are provided in Appendix A.1.
Surgical schedule optimization can thus be operationalized in practice by solving deterministic BOP sequentially for b = 1, 2, . . . , B, with an appropriately chosen sequence of scheduling days, {s b } and estimates of patient LOS, {l p : p ∈ P b ∪ P past b }. For patients in P b , estimations of l p can be made based on patient features at the time of arrival (see Section 4). For patients in P past b who have received their surgeries, l p estimates can be progressively updated as uncertainties in LOS realize. We formalize this information update procedure below.

Information Update Procedure for deterministic BOP.
At the start of each batch b, the value of l p for all p ∈ P past b is updated as follows.
(a) If the patient has undergone surgery and has been discharged from ICU by scheduling day s b , her realized LOS is observed and we update l p to be the true LOS. For simulation purpose only, we use a simplified heuristic for patients under (b). We assume the true LOS can be evaluated with certainty for those soon to be discharged (within 5 days), and thus update l p to be the true LOS for those patients. For the remaining patients, we set l p ← max{m + 10, l p }. In practice, l p can be reevaluated on a case-by-case basis for each patient based on any new information available. On every scheduling day, patient arrivals P b and any realization of additional LOS information for P past b are observed. An optimal solution of deterministic BOP is then obtained and implemented for all patients in P b . The existing schedule and required surgical resources are updated accordingly. We describe this rolling optimization process in Algorithm 1.

Algorithm 1 Rolling Deterministic Optimization (RDO)
Update LOS estimate l p using the Information Update Procedure.
implement the optimal solution x * for all p ∈ P b .

8:
3. Schedule and Capacity Update. 9: In practice, the sequence of s b is an additional design choice that can be set based on hospital practices and nature of different procedures. In all our experiments (both the deterministic formulation and stochastic or robust formulations in later sections), we set scheduling days to be the start of every month, and schedule patients in monthly batches given our context. For complex, elective cardiovascular procedures, there is usually a significant lead time between patient referral and when a patient undergoes surgery. For our patient cohort, the median lead time is 59 days and the average lead time is 79 days. Therefore, monthly scheduling is chosen as an example and proof of concept given the context: although patients need to wait for being scheduled, the effect of this additional wait time is limited in this case as the total wait time is being minimized under our optimization algorithms. In Appendix C, we also re-run our numerical experiments for biweekly scheduling days and show that the resultant insights are similar. For other types of procedures that require faster scheduling, s b can be adjusted accordingly such that scheduling can happen more frequently in smaller batches. We discuss the implications of this design further in Section 7.

Model Performance under Perfect Information
Before introducing machine learning approaches for estimation of LOS, we first evaluate the performance of our models in the real world setting of the PAMC, using true historical LOS as model inputs, {l p }. The purpose of this exercise is to identify the upper bound of the achievable level of operational performance versus the status quo, which sets the optimal performance benchmark for remaining numerical experiments.
We evaluate the performance of our scheduling model by considering the 596 elective cardiovascular surgeries that were carried out at the PAMC from September 2018 to March 2020. Patient arrivals are modeled using deterministic, historical arrivals. An optimal schedule was generated using the optimization models given the true LOS of each patient, and the resultant daily CVICU census was simulated based on the optimal surgical schedule and the actual LOS.
The performance of the optimization models is assessed using two metrics: 1) percentage reduction in the number of high ICU occupancy days (10 or more elective patients) in the ICU versus the historical status quo, and 2) median and average reduction in patient wait times compared to the original schedule. Good performance corresponds to a measurable reduction in the number of high ICU occupancy days without increasing patient wait times.
When evaluating our models, we treat the first six months  Overall, both models achieve significant reduction in both high occupancy days in the ICU and patient wait times.
The above results show that it is theoretically possible to eliminate the number of days with more than 10 elective patients in the ICU given perfect information on LOS at the time of scheduling. Such performance improvement can be achieved along with a reduction in average and median wait times.
Having established the performance upper bound under perfect information, we next proceed to consider incorporating machine learning predictions of LOS in schedule optimization, and the impact of inaccurate predictions on deterministic optimization models.

Predicting LOS with Machine Learning
Having established the performance upper bound of schedule optimization using accurate LOS predictions in the previous section, we now present the development of machine learning models for LOS prediction in this section, as well as the schedule optimization results using the predicted LOS.
We discuss both regression and classification models for LOS prediction. 5-fold cross validation is used for hyper-parameter tuning and model selection during the training process. The model with the highest cross-validation score is selected for validation on the test set. The sklearn library [20] in Python was used to train and implement the selected models.
The predicted LOS is then used as input for RDO. Scheduling optimization is only performed on the test set and not on the training set. This simulates the real world setting where the incoming patients are prospectively scheduled based on LOS predictions at the time of admission. We compare the schedule optimization result with the status quo.

The Dataset
Our data set consists of medical records, surgery data and CVICU census data for 2,352 elective cardiovascular surgical patients at the PAMC from 2014 to 2020. It is split into a training set with 1,738 patients from January 2014 to mid-June 2018, and a test set with 614 surgical patients from mid-June 2018 to May 2020.
The patient features are chosen based on clinical knowledge of the experienced surgeons at the unit, as detailed below. The missing values in the data are imputed using the mean for continuous features and the mode for categorical features from the training set.
To compare our data quality with the existing literature, we first train a binary classification model to identify patients with risks of prolonged stays of more than 5 or 10 days, making up 28% and 12.8% of the training set respectively. Candidate machine learning models considered include logistic regression, random forest and gradient boosting machine. Gradient boosting machine is selected based on cross-validation AUC.
We refer to Ettema et al. (2010) [6] for performance benchmarks. The authors of [6] reviewed and validated 14 models for the prediction of prolonged LOS in the ICU after cardiac surgery using data from 11,395 surgeries, and found that the area under the curve (AUC) scores of the best performing models range from 0.71 to 0.75. Our model achieves an AUC score of 0.77 for stays more than 5 days and 0.73 for stays more than 10 days on the test set, which is inline with that of the best-performing models identified in [6] (See Table 1). Although a direct comparison with existing literature is difficult due to different patient populations and different definitions of prolonged stay, our model and features can achieve high levels of distinguishing performance compared to the benchmarks.

Model Development
We now train machine learning models for LOS with the objective of using the predictions for surgical scheduling.
We first consider classification models that classify LOS into ordinal buckets: 0-1 day, 2-5 days, 6-10 days and > 10 days. For each LOS bucket, we define two metrics to evaluate model performance. 'Accuracy 1' is calculated as the exact number of true predictions in the bucket divided by the actual number of patients in the bucket. 'Accuracy 2' extends the definition of true prediction to include the scenario where the patient is predicted to be in the LOS bucket adjacent to his or her true bucket.
Candidate models include gradient boosting machine (GBM), random forest and ordinal regression models. GBM achieves the best cross-validated accuracy score on LOS classification with multiple LOS buckets. The model achieves an overall accuracy of 53%, and 88% of the predicted LOS either fall in the true bucket or the adjacent buckets. The most significant predictors according to the impurity-based feature importance are procedure types (0.47), weight (0.22), height (0.17) and the surgeon (0.07).
A breakdown of the predictive accuracy on each buckets of the test set is presented in Table 2. While predictive accuracy for the first three patient groups with LOS ≤ 10 days is relatively high based on Accuracy 2, predictive accuracy drops sharply for the group of patients with longer than 10 days of LOS, where the majority of the patients' LOS is being under-estimated. Similar behavior is observed under other choices of bucketing, such as using finer buckets for the LOS> 10 patient group.  We use the predicted buckets of LOS are as inputs for RDO in Algorithm 1. Given each patient's predicted LOS bucket, the upper bound of the first three buckets is used as the LOS point estimates l p , and l p = 30 is used for the > 10 bucket because 95% of the patients left the CVICU in less than or equal to 30 days. Although the machine learning model uses month of surgery from historical data, this feature was studied and had little effect on LOS prediction. Optimization performance on ICU occupancy and wait times is shown in red in Figure 3, in comparison with performance upper bounds in grey assuming the classification algorithm achieves 100% accuracy. Here, we consider a range of weights on the cost of ICU overflow relative to the cost of total wait times: β ∈ {200, 100, 50, 25, 10, 5}. With perfect predictions of each patient's LOS bucket (grey bars), the optimization model significantly reduces the number of high ICU occupancy days without pushing back surgeries when β = 10 and 5. In contrast, the machine learning predicted outcomes show only slight improvement in both objectives for β = 50 (red bars). The trade off between the two objectives becomes significant as β decreases. When the weight on the cost of ICU overflow is sufficiently high, optimization reduces the number of high ICU occupancy days at the cost of increasing patient wait times and vice versa.

Bucket (% of Test
Compared to classification models, generating pointforecast for individual patients using regression models is even more challenging especially for those with prolonged LOS (see, for example, [26]). We follow the same procedure to select and train a variety of regression models including OLS regression, Lasso regression, gradient boosting machine, quantile regression and more. However, when the generated point-predictions are combined with deterministic optimization, there is still no model that achieves performance improvement on both metrics of interest. Similar to the classification model above, the main difficulty comes from significant underestimation of prolonged LOS. Figure 4 plots the density distribution of relative predictive errors, defined as the ratio between true LOS and predicted LOS, for training and testing set. The distribution is heavily right-skewed with a long-tail of very large relative errors (i.e. true LOS >> predicted LOS). We explore this issue further in the following section.   to be inferior in the tail of the distribution to that in the middle due to imbalanced data. Furthermore, [13] shows that after 5 days in the ICU, the most significant predictors of the remaining LOS are features collected on day 5 instead of those at the time of admission.
Poor predictions of prolonged LOS have a significant impact on optimization performance, especially in the context of cardiovascular surgeries where LOS in the ICU tends to have a long tail distribution. As shown in Figure 5b, 5.7% of the patients with LOS of greater than 30 days in the ICU makes up 42% of the total CVICU daily census count in the test set aggregated from 2018 to 2020. This small group of patients thus has a disproportionately large influence on ICU resource use and operational performance.
A closer examination of our data also reveals the non-stationary nature of hospital operations, another challenge to developing and evaluating predictive models in real life. A comparison between Figure 5b and 5a highlights the shift in patient population in terms of LOS distributions from 2014-2018 to 2018-2020 at the PAMC. The most notable change is an increase in cases with extended LOS in the CVICU, driven by the hospital's decision to admit a larger fraction of complex cases in recent years. Such shifts in hospital operations make prospective predictions of LOS based on past data even less reliable. We split the training and test set in chronological order instead of randomly to reflect the impact of such inherent non-stationarity on model development and performance in practice.

Data-Driven Surgical Scheduling with Machine Learning Under Uncertain LOS
While combining machine learning with a deterministic formulation of optimization can provide significant performance improvement in theory, we have shown that the unpredictability of the long-tailed LOS along with other challenges can lead to poor performance in practice. A natural next step is to consider optimization methods that directly incorporate the uncertainties in LOS.
In this section, we propose a data-driven surgical scheduling framework, designed to address the identified challenges in the previous sections. Our framework combines optimization under uncertainty with machinelearning predicted distribution of LOS, and rolling information update.
Under this general framework, we develop three algorithms: Standard-RSO, Conservative-RSO and RRO. Standard-RSO and Conservative-RSO apply stochastic optimization using the predicted LOS distributions to minimize the expected cost, where the latter further adjusts the predicted LOS distribution to target underestimations for prolonged LOS. RRO uses robust optimization to produce solutions that are robust against under-estimations for prolonged LOS.

Estimating LOS Probability Distribution using Machine Learning
Instead of relying on deterministic predictions of LOS for optimization, we now use machine learning models to predict the probability distribution of LOS. We start by generating point-wise prediction of LOS for each patient using conventional machine learning models discussed in Section 4. For each patient p, let v p be the value of LOS point prediction from the machine learning model for the true LOS, L p , and define the relative prediction error L p and R p are random variables and v p is known. We assume that the random variables R p 's are independent and identically distributed across patients from a distributions S, i.e., Note that R p > 1 represents an under-prediction of LOS and R p < 1 means an over-prediction. We choose to model relative prediction errors so that the absolute prediction errors, |L p − v p |, tend to be larger for patients with prolonged realized LOS. This is in accordance with our finding in Section 4 that prediction is more challenging for prolonged LOS. The distribution S is approximated using the empirical distribution of predictive errors of the machine learning model during training and validation. Figure 4 provides an example of this distribution.
Given estimated S, the distribution of L p for individual patient p is estimated using

Rolling Stochastic Optimization (RSO)
Building on the rolling deterministic optimization framework introduced in section 3.2, we first explore stochastic optimization approaches that utilize estimated probability distributions of patient LOS.
Let the random variable L p denote the LOS of each patient p, and a LOS realization trace, denote a sequence of realized LOS for all patients who have arrived at the system at the start of period b. We use µ b to denote the discrete probability distribution of ω over the set of all possible traces for batch b, denoted as Ω b . We assume that each patient's LOS is independent of other patients, and hence µ b is the joint distribution of independent random variables L p , for all p ∈ P b ∪ P past b .
Using the same set of notations as deterministic BOP introduced in Section 3.2, stochastic optimization formulations aim to minimize the expected cost taken over the distribution µ b : Equation (5) parallels the deterministic objective function of deterministic BOP in Equation (3a). The distinction is that the ICU overflow variable, u d is now a random variable that is a function of x and ω.
Following the approach used in [17], we approximate the objective function (5) The objective function (6a) formulates the SAA approximation of Equation (5). Constraints (6d) and (6e) formulate the stochastic parallel to constraints (3d) and (3e) in deterministic BOP, creating a replica of variables y d,p and u d,p for every sample trace ω (n) . The other constraints remain identical to their deterministic counterpart. Note that we focus on surgical scheduling problems where different realizations of ω only affect daily ICU occupancy and overflow, and do not affect feasibility of a schedule, x. In practice, ICU overflow can be accommodated by setting up temporary ICU beds or utilizing spare resources from other ICUs. Feasibility is thus determined exclusively by operational constraints such as OR room capacity, surgeon and patient availability, which remain unchanged. We combine stochastic BOP with the rolling optimization framework introduced previously for RDO. Under this framework, stochastic BOP is solved sequentially for a sequence of schedule days, {s b : b = 1, . . . , B}. Meanwhile, the probability distributions of LOS L p are progressively updated for every patients as uncertainty realizes. We introduce the information update procedure for stochastic BOP below.
Information Update Procedure for Stochastic BOP. At the start of each batch b, the distribution of L p for all p ∈ P past b is updated as follows.
(a) If the patient has undergone surgery and has been discharged from ICU by scheduling day s b , her realized LOS is observed and we update L p to be a constant equal to the true LOS. (b) If the patient is in the ICU on day s b having stayed for m days, then we update the distribution of L p using the conditional distribution of L p |L p ≥ m.
To obtain the resultant distribution, we update the distribution of prediction errors, R p ∼ S by conditioning on round(v p · R p ) ≥ m while fixing v p . (c) If the procedure of the patient is scheduled after s b , there is no change to the distribution of L p .
We develop two rolling stochastic optimization algorithms, Standard-RSO and Conservative-RSO. Both algorithms adopt the same framework combining stochastic optimization, machine-learning predicted distribution and rolling information update.
The two algorithms differ from each other on how LOS realization traces ω are sampled for SAA. In Standard-RSO, realizations of l (n) p are generated by sampling the relative prediction error r (n) p directly from S and setting In contrast, in Conservative-RSO, we generate conservative realizations of l (n) p by rounding any r (n) p < 1 drawn from S to one, i.e., In other words, Standard-RSO draws LOS realizations from a distribution derived from two-way predictive errors, S; Conservative-RSO draws LOS realizations from an adjusted distribution of S, where only one-way errors from LOS under-predictions (i.e. r p > 1) are retained, and any over-prediction errors r p < 1 are rounded up to 1. Conservative-RSO is designed such that the algorithm focuses on stochasticity arising from difficultto-predict extended LOS for cardiovascular surgeries, where machine learning models consistently yield underpredictions.
We formalize the two rolling stochastic optimization algorithms in Algorithms 2 and 3 below.  Update LOS distribution L p and S using the Information Update Procedure for stochastic BOP.  4. Schedule and Capacity Update. 13: When implementing Standard-RSO and Conservative-RSO, we set N Ω = 10 to obtain our numerical results. Large values of N Ω can lead to very large numbers of variables and constraints under Equations (6d) and (6e). We thus limit the value of N Ω due to computational constraints. We discuss the computational limitations further in Section 7.

Rolling Robust Optimization (RRO)
Robust optimization formulations provide an alternative to stochastic optimization for scheduling under uncertainty. Robust optimization aims to minimize the worst-case cost defined over an uncertainty set, U b , which Update LOS distribution L p and S using the Information Update Procedure for stochastic BOP.   }. Using the same notations as before, the formulation of robust BOP for a pre-defined uncertainty set U b can be written as follows.
where θ(x) is the worst-case cost associated with ICU overflow for all possible LOS realization traces in the uncertainty set, Apart from the uncertainty set U b , the remaining constraints of robust BOP are similar to those of deterministic BOP in Equations (3). It is worth noting that in Equation ( and, for some chosen constant 0 ≤ η ≤ 1, Constraint (9) enforces a budget of uncertainty, η · |A b |, which limits the total possible normalized deviations (i.e., extended ICU days) from the lower bound l min p . Following the terminology in robust optimization literature (see, e.g., [2]), we refer to η as the uncertainty budget, and tune the value of η among {0.5, 0.75, 1.0} based on performance of numerical experiments.
Instead of using stylized assumptions to set values for l min p and l max p and construct the uncertainty set (see, e.g., [19] where F −1 S denotes the inverse cumulative distribution function of S. For our numerical experiments, we present results with α = 0.2 (i.e., the 80th percentile), chosen among {0.1, 0.15, 0.2, 0.25} based on simulated robust optimization performance for our data set.
As uncertainties of LOS realize overtime for previously scheduled patients, P past b , we update the distributions of LOS and parameters used for the uncertainty set dynamically as follows. We develop an adapted column and constraint generation approach (AC&CG) to solve robust BOP in (7) with the above definition of U b . The method of AC&CG for surgical scheduling with ICU capacity constraints was first proposed in [19], which was designed for linear cost functions f (u d ) only. Our AC&CG approach extends the algorithm in [19] to accommodate convex, piece-wise linear formulations of f (u d ).

Information
At a high level, the AC&CG approach involves iteratively solving a Main problem and a Recourse problem. Every iteration t, the Main problem solves for a temporary optimal scheduling decision, x t * , that minimizes the worst-case cost over a subset of the uncertainty set, Ω t ⊂ U b . The optimal objective value provides a lower bound to that of the original BOP. The Recourse problem then finds a worst-case LOS realization trace, ω t * ∈ U b , that maximizes the cost of ICU overflow under the temporary scheduling decision, x t * . The optimal objective value provides an upper bound to that of the original BOP. When the upper and lower bounds obtained are equal, the temporary scheduling decision x t * obtained from the Main problem is also optimal to robust BOP, and the algorithm terminates. Otherwise, we update the set Ω t to include ω t * obtained from the Recourse problem, Ω t+1 ← Ω t ∪ {ω t * }, and re-solve the updated the Main problem with Ω t+1 .
We formalize the above AC&CG algorithm for solving robust BOP in Algorithm 5 in Appendix B. Note that each iteration of AC&CG adds a LOS realization trace, ω t * to the set Ω t . This process adds a significant number of new variables and related constraints to the Main problem. 5 As the number of iterations increases, 5 See Equations (23d) and (23e) in Appendix B.
the size of the Main problem can therefore grow quickly. When implementing AC&CG, we thus terminate the algorithm after at most 10 iterations due to computational constraints. In Appendix C, we show that the optimally gaps, U B−LB LB , are reduced to under 1% for almost all batches after 10 iterations. We discuss the computational challenges further in Section 7.
The full rolling robust optimization algorithm (RRO) combining information update and AC&CG is presented in Algorithm 4. the test set 6 , are used as the estimated distribution S. The LOS predictions and S are then used as inputs to Algorithm 2, 3 and 4 to generate optimal surgical schedules. Performances of the optimal schedules are simulated and evaluated using historical patient arrivals and actual LOS data in the test set from September 2018 to March 2020.
Our numerical results indicate that both Conservative-RSO and RRO outperform the status quo, while Standard-SP fails to achieve consistent performance improvement. Moreover, Conservative-RSO achieves the best overall performance despite its relative computational simplicity compared to RRO. The results highlight the importance of tailoring algorithm for long-tailed distributions and reveals shortcomings of complex algorithms in practical settings. The best-performing algorithm, Conservative-RSO, provides a promising direction for designing efficient CVICU scheduling algorithms.

Performance on Historical True LOS
We first simulate the performance of Standard RSO, Conservative-RSO and RRO on historical patient arrivals.
Each run of either stochastic optimization algorithms may produce a different scheduling policy, because the objective function involves random sampling. For this reason, both Standard-RSO and Conservative-RSO are run 90 times each on the testing period for each values of β ∈ {1, 5, 10, 20}, where greater values of β mean more weight on the cost of ICU overflow compared to wait time.
In contrast, the solution of RRO is deterministic for the same set of model parameters. When evaluating RRO, we thus only run the algorithm once for each set of parameters. For RRO, we tune both values of β ∈ {1, 5, 10, 20} and the uncertainty budget, η ∈ {0.5, 0.75, 1.0}. Greater values of η mean a larger uncertainty set containing longer LOS realizations.
Performance of each optimal policy is simulated using the historical true LOS trace, i.e., the actual realized values of LOS for each patient in the testing period. The change in wait times and the number of high ICU occupancy days in the ICU for each experiments are 6 In practice, prediction errors of the test set will not be available at the time of scheduling. The test errors are included in constructing S so that, if the LOS of an incoming patient in the test set surpasses all cases in the training set, conditional sampling described in the previous section (see Step 2(b) of Method 1 still works as intended. Although this is not ideal, the scheduling process should not benefit much from it and it helps simplify our simulations. In practice, if the LOS of a patient surpasses all previous cases, one approach is to consult health providers. evaluated in comparison to those of the original schedule (i.e. the status quo). Good performance corresponds to an improvement on both metrics.
Performance of the Standard-RSO, Conservative-RSO and RRO are presented in Figure 7, 6 and 8 respectively. Each run of experiments is plotted using the average or median change in wait time (the y coordinate) and the simulated reduction in high ICU occupancy days in the ICU (the x coordinate). The quadrant shaded in blue in each plot indicates performance improvement on both metrics. Different values of β indicated by color, and different values of η (for RRO only) are indicated by marker type.
Standard-RSO shows poor performance in reducing ICU overflow regardless of the value of β, 7 as illustrated in Figure 6a and 6b. Despite considering the stochasticity in LOS, Standard-RSO tends to aggressively schedule most patients much earlier and is unable to effectively handle the trade off between ICU overflow and patient wait times.
In contrast, Conservative-RSO demonstrates consistent improvement of performance in both metrics compared to the status quo in Figure 7a in terms of the average change in wait time and the reduction in high ICU occupancy days for β ∈ {5, 10}. Relatively weaker performance is observed when the median change in wait times is examined in Figure 7b. This is likely because the formulation of the objective function specifically minimizes the sum of wait times instead of the median. As expected, higher values of β lead to less high ICU occupancy days but longer wait times.
RRO also achieves good performance shown in Figure 8. In particular, performance improvement on both metrics parameter are achieved under parameter combinations (η = 1, β = 10) and (η = 0.75, β ∈ {5, 10, 20}). Higher values of β and uncertainty budget η tend to result in less high ICU occupancy days but longer wait times, in line with our expectation.
Although Conservative-RSO and RRO are more conservative in its ICU occupancy forecasts and lead to longer patient wait times, there is enough slack in the original system so that both are able to effectively reduce ICU congestion without excessively pushing back surgeries compared to the original schedule. In addition, the contrast in performance between Standard-RSO and Conservative-RSO further reflects that, with careful design choices targeting the challenges specific to the problem, performance improvement can be achieved 7 For consistency, we present performance for β ∈ {1, 5, 10, 20} for all algorithms. It is worth noting that increasing the value of β for Standard-RSO does not lead to meaningful improvement in performance. One example is provided in Figure 12 of Appendix C. in practical settings without increasing complexity of the algorithm

Performance on Bootstrapped Traces of LOS
In order to obtain confidence intervals for the relative performances of the three methods, we now use bootstrapping [5]   Each algorithm produces one optimal schedule x * for the sequence.  • Step 4. Simulate the performance of the optimal schedule by each algorithm under the LOS realization scenario described by ω (i) . Compare relative performance on both metrics. Performance of RDO is used as a benchmark for comparison.
Similar to Section 6.1, Step 1-4 are repeated for 100 iterations (N E = 100) on each value of β for each algorithm. For RRO, η = 0.75 is used because of its good performance on the historical trace shown in Figure 8. We compare the change in wait times relative to the original schedule, and compare the number of high ICU occupancy days of the three algorithms relative to that of deterministic optimization, RDO. Figure 9 presents the performance of Standard-RSO, Conservative-RSO and RRO in reducing ICU congestion on bootstrapped patient arrivals in comparison to the performance benchmark set by RDO.
Standard-RSO manages to achieve some overflow reduction for β ≥ 10. When β is small, however, Standard-RSO tends to under-perform compared to the deterministic benchmark. In contrast, the interquartile ranges of Conservative-RSO (blue) and RRO (red) both show more significant reduction in the number of high ICU occupancy days compared to RDO for all values of chosen β.
The median of changes in patient wait time compared to the status quo are presented in Figure 10. Standard-RSO demonstrates similar behaviors to its deterministic counterpart, with the tendency to aggressively schedule patients earlier. In contrast, Conservative-RSO strategically delays procedures so that ICU overflow can be effectively avoided. For β ∈ {1, 5, 10}, the Conservative-RSO is able to consistently reduce ICU overflow without pushing back a majority of procedures compared to the original calendar. RRO, on the other hand, achieves comparable performance to Conservative-RSO and is able to reduce ICU overflow without pushing back procedures for all values of β.
Comparing Conservative-RSO and RRO, both algorithms utilize conservative LOS estimates to address under-predictions for long-tailed LOS. Both achieves promising reductions in ICU congestion without excessive increase in wait times. However, Conservative-RSO demonstrates better average and worst-case performance (upper caps of the box plots) in ICU overflow for all selected values of β.
The relative poorer performance of RRO in reducing ICU overflow may be a result of its sensitivity towards the choice of other parameters, such as the uncertainty budget η, the choice of l min p and the percentile 1 − α used for setting l max p . During our effort to calibrate the uncertainty set, the performance of RRO is observed to be drastically different under different parameter choices. While there remains room for potential performance improvement with better parameter tuning, finding the optimal combination of all parameters can be a very difficult task in practice given the large number of possibilities and the long computation time required for each run of the AC&CG algorithm. This reveals a shortcoming of using robust optimization approaches that involve a large number of parameters, and demonstrates that increasing algorithmic complexity does not always leads to better performance in practice. To improve practicality of robust optimization, further research is needed to develop new formulations less sensitive to parameter choices.
Lastly, we note that there is still a notable performance gap between our best-performing model, Conservative-RSO, and the theoretical upper bound established in Section 3 with rolling deterministic optimization under perfect information. This implies significant room for further development of more efficient data-driven optimization methods, with Conservative-RSO serving as a stepping tone towards a promising direction.

Summary of numerical results and additional approaches
We examined a variety of optimization formulations in conjunction with data-driven LOS predictions. All approaches are developed under similar frameworks combining machine learning, rolling information update and optimization methods, but they differ significantly in their simulation performances, strengths and weaknesses. We provide a summary of these algorithms in Table 3.   [15] to search for other promising models from its wide range of builtin machine learning and deep learning algorithms.

RDO
None of the above alternatives result in consistent performance improvements compared to results presented in this section.

Insights and Discussions
Among all the practical formulations considered, Conservative-RSO and RRO both achieved reduction in both ICU congestion and patient wait times. The key driver of their promising performance is the focus on under-estimates of long-tailed, extended LOS of cardiovascular surgery patients, through design of either SAA sampling of LOS realizations, or the uncertainty set of LOS. Good performance also relies on precise choices of the trade-off coefficient β, and several other model parameters in the case of robust optimization. To achieve the desired optimization performance in practice, all model parameters need to be carefully tuned to balance reduced ICU congestion with longer wait times.
Moreover, on each scheduling day, possible LOS outcomes for patients currently in the ICU are updated given how long they have stayed in the ICU. In particular, if some patients are already staying longer than expected, their projected probability distribution of LOS will be updated accordingly based on the new information. This information update procedure further corrects for potential under-estimates of prolonged LOS at the time of scheduling and improves optimization performance.
The fact that the optimization outcomes are highly sensitive to these nuanced design choices underscores the operational challenges of surgical scheduling with long-tailed LOS. Any solution will need to be tailored and fine-tuned to suit the context of specific institutions in order to reach the desired outcome.
Although our proposed scheduling framework has demonstrate promising performance, there remains room for potential extensions and future work to address some of its limitations.
First, under our rolling optimization framework, patients who arrived in-between two scheduling days need to wait until the next scheduling day to be scheduled. In our simulations, we have chosen the scheduling days s b to be one month apart for elective procedures based on our context, as explained at the end of section 3.2. However, for other types of procedures, scheduling may need to happen on a daily bases. While our modeling framework naturally extends to the case of daily scheduling (by setting s b as consecutive days), managing ICU overflow in this case can become more difficult, since scheduling is performed with less patient arrival information. This is a limitation of our framework to be addressed in future works. For instance, authors of [12] propose an alternative rolling horizon framework for surgical scheduling that takes into account uncertain future patient arrivals. With a forward-looking, rolling arrival horizon, scheduling decisions can be made more frequently for realized patient arrivals, in anticipation of stochastic future arrivals. Incorporating designs of this kind with our current framework has the potential to improve practicality and performance for procedures that require more urgent scheduling.
Second, a challenge of both Conservative-RSO and RRO that has not been addressed is their computational limitations. In our numerical experiments, we have limited the number of sampled traces, N Ω for SAA to 10 and the number of iterations of AC&CG to 11. As mentioned in Section 5.2 and 5.3, the memory and run time required to solve the optimal schedule for each batch of patients can grow quickly when N Ω and the number of AC&CG iterations increase. In the case of SAA for Standard and Conservative-RSO, it typically took 4 to 8 hours to run each numerical experiment for N Ω = 10 with 19 batches with two CPUs and 16GB memory per CPU. For N Ω = 30, at least 21 hours were required to run one such numerical experiment with four CPUs and 16GB of memory per CPU. Meanwhile, AC&CG for RRO requires even longer run time: it took typically 7 to 12 hours to run each experiment for 11 iterations with two CPUs and 16GB memory per CPU. The authors of [19] also point out that the proposed algorithm may be computationally inefficient for larger uncertainty sets. Developing a better understanding of how the choice of N Ω and the maximum number of AC&CG iterations impacts optimization performance is a meaningful subject of further research.

Concluding Remarks
Using seven years of data from an academic medical center, we developed machine learning models to predict post-surgical length of stay, optimization models to schedule procedures to minimize post-surgical bed congestion, and simulated the results of the use of these models on the most recent 18 months of held-out empirical test data. We established an upper bound on performance with an offline optimization model of historical data with access to actual LOS.
A conservative stochastic program with sufficient sampling of the LOS distribution to capture tail behavior managed to achieve the best overall performance in reducing ICU congestion without increasing wait times for surgery. Compared to the hospital's current paperbased system, the deterministic and standard stochastic optimization approaches, along with numerous variants of machine learning and optimization formulations, failed to improve performance. The failure of most models to improve over the status quo, illustrate the importance of using empirical data, rather than synthetic parametric data, especially for long-tailed distributions. Several lessons for dealing with such data are offered by the contrast between the models that did and did not improve performance.
The prevalent null results highlight the challenge of the unpredictability and non-stationarity of prolonged LOS in practical settings. This challenge is rarely addressed in previous works, most of which follow the common practice of using synthetic LOS data generated based on strong distributional assumptions. For example, in [23] and [17], the distribution of the LOS in the ICU is assumed to have a mean and standard deviation of no more than 3.5 days. In the robust optimization approach discussed in [19], the worst-case deviation of the uncertainty set is assumed to be 4 days, which the authors admit can be far from reality. Despite the positive results obtained in simulation, the actual performance of these optimization models can be drastically different in practice in the presence of significant LOS estimation errors.
We proposed a novel formulation of stochastic optimization, Conservative-RSO, which demonstrated most consistent improvements with empirical data despite of its relative algorithmic simplicity. Our results show that Conservative-RSO is a promising first step towards addressing the unpredictability and long-tailed behavior of LOS to improve surgical scheduling, and they provide meaningful directions for future work. First, much work is needed to explore and evaluate alternative designs of data-driven optimization formulations, calibrated and evaluated using empirical data of LOS. When empirical data are not available, any distributional assumptions should incorporate the long-tailed, non-stationary nature of LOS distributions to minimize the bias present in using standard distributional assumptions. Furthermore, learning over long-tailed distributions has not been fully examined in prior research. We highlight its central importance in tackling the general problem of surgical scheduling with constrained downstream capacities. Obtaining accurate and precise LOS predictions is one way to push the Pareto frontier of optimization closer to the theoretical performance upper bound using full information. Settings where patient characteristics available at the time of scheduling fail to explain the majority of variation in LOS offer exciting challenges to be tackled by innovative approaches that combine the power of prediction, optimization and information update. Solving the quadratic objective, f (u d ) = u 2 d of the mixedinteger program can be slow. To reduce the runtime required to solve the problem, we implement the quadratic term in the objective function, f (u d ) = u d , using a piece-wise linear approximation, f (u d ) = e1u (1) d + e2u (2) d + e3u (3) d + e4u (4) d + e5u and add additional constraints In our formulation, we set e1 = 1, e2 = e3 = e4 = e5 = 2. Since u d only takes integer values, this coefficient choice leads to f (u d ) ≡ u 2 d for any u d ≤ 5. The piece-wise linear approximation is 20 times faster than the quadratic form in our numerical experiments.
The above approximation is used in all our algorithms. The formulations for stochastic and robust algorithms are analogous (e.g. copies of u ( m) d are created for different LOS realizations), and we thus omit the details here.
Next, we introduce the full mathematical formulation of Q op in constraint (1e) for offline optimization below.

Sets
Next we describe the optimization constraints (12)-(16) that are equivalent to x ∈ Q op .

A.1 Q op for Batch Optimization Problem
When solving BOP for rolling scheduling, Q op is adjusted accordingly. In addition to the Sets and Parameters specified in Section A, we use the following sets and parameters. Sets