Appointment Scheduling at Outpatient Clinics Using Two-Stage Stochastic Programming Approach

Clinics with a large volume of patients are often burdened with limited resources such as nurses and providers. Also, an efficient health system seeks short wait times for the patients to see the provider (indirect wait time) and within the clinic (direct wait time) during the day of the appointment. Additionally, the appointment duration, volume of patients, no-show behavior are uncertain. The direct and indirect wait times, stochastic parameters, rising treatment costs, and increased demand of patients motivate the need for efficient appointment schedules and clinic operations. In this paper, we develop a two-stage stochastic mixed-integer linear programming model (SMILP) integrated with a simulation model to generate a scheduling template for the providers to schedule individual patient appointments and resources. The model minimizes the expected wait times for the patients with a fair and equitable utilization of the resources. Computational experiments were conducted using a data-driven simulation model, and the results indicate that the proposed approach can significantly decrease patients’ direct and indirect wait times when compared to a deterministic indexing policy used for scheduling appointments.

to overcome the following problems: the no-show behavior of patient arrival, patient/provider check-in delays, overbookings, long wait times, and poor provider/staff utilization. These are pervasive in outpatient clinics ([12]- [21]). On the other hand, appointment scheduling systems, which give patients flexibility in choosing their appointment time, not only lead to satisfied patients but also have outstanding effects on other patient flow metrics such as decreasing the no-show rate as well as patient appointment delays (time between patient desired time and assigned appointment time) and higher patient retention rates. All of these result in better reimbursement rates by payers for providers ([22], [23]).
Many studies have documented the no-show rate in medical practice. The study in [19] reported a 42% average no-show rate which ranges from 6% to 92% in outpatient clinics. Also, the study in [20] reports 13% to 24% no-show rates at endoscopy clinics for different service types, and [24] indicates post intervention no-show rates ranging from 28% to 45%. The results in [25] show the overall no-show rate at Obstetrics and Gynecology (OBGYN) clinics as 30.1%. They investigated the strong relationship between patients' appointment delays and no-show cases in OBGYN specialty clinics, and similarly in psychotherapy appointments, a 21% no-show rate was reported by [21].
Another patient flow metric is patients' appointment delays known as indirect wait time in the literature. A survey was conducted by [26] to determine the average indirect wait time for new patients to see a provider in 15 major and 15 midsized metropolitan areas in different specialty clinics. The survey was organized in 2004, 2009, 2014, and 2017, and the results show an increase in the indirect wait time in 2017 compared to other years. For instance in 2004, it was reported for 15 mid-sized metropolitan markets between 88,000 and 143,000 people including 1,414 medical offices in large metro markets and 494 medical offices in mid-sized metro markets. The report on indirect wait time for cardiology, dermatology, obstetrics and gynecology, orthopedic surgery and family medicine, and average indirect wait time of obstetrics and gynecology was studied in [26]. The average is about 25 days with a standard deviation of 10 days. In this study, we develop a framework to minimize direct and indirect wait times simultaneously while considering fair and equitable utilization of the resources in a clinic. A case study based on the operations in OBYGN clinic is also presented to corroborate the performance of the approach.
In this paper, we present an integrated approach of twostage stochastic and simulation models to generate scheduling template while considering the sources of uncertainties such as no-show behavior, direct and indirect wait times, and fair and equitable assignments of new patients to the providers. Following are the major contributions of this paper: 1) a two-stage stochastic model to generate a weekly scheduling pattern for individual providers that gets updated at regular intervals based on the type and mix of services rendered; 2) a clinic simulation model to determine sequencing rules to maintain the direct wait times within a threshold; 3) extensive data-driven simulated computational experiments to corroborate the effectiveness of the proposed method with a first-call-first-serve policy used in the call center. The overall framework is depicted in Figure 1. This paper is organized as follows. Section II reviews the relevant literature and describes the problem. Section III discusses the development of our model formulation and solution scheme. Section IV presents the computational experiments and insights of the results based on data driven simulation study for an OBGYN clinic. Finally, Section V provides concluding remarks.

II. RELEVANT LITERATURE
In recent years, appointment scheduling in outpatient clinics has attracted much attention in health care delivery systems. The authors in [27] and [28] discussed a variety of methods in modeling, optimization, and future work in appointment scheduling. Designing an effective appointment scheduling system in outpatient clinics helps a smooth flow of patients while considering patients' and physicians' preferences, and matching supply and demand. The research in [29] developed a framework to find the possibility that patients may cancel or no-show at their time of appointments. The authors in [28] divided decision making in outpatient appointment scheduling into static and dynamic models. The work in [28] and [30] define static scheduling as the decisions about appointment times that are made prior to the start of the appointment session while in dynamic case, the schedule may be modified later depending on the state of the system. The work in [31] considered a static appointment scheduling where the number of patients is already known. The work in [29] developed a dynamic scheduling where appointments are assigned to patients depending on the clinic's appointment schedule at the time of the patient's call. In our work, we consider the dynamic model for our computational experiments, and the proposed two-stage model proposes a scheduling template for the call center scheduler.
Another classification on appointment scheduling is with respect to the wait times. The work in [29], [30], [32], [33], and [34] are among the most recent studies on appointment scheduling which consider wait time in model formulations. Other than wait times, the research in [34] addresses crucial characteristics such as the randomness of service time, patient punctuality and no-show behavior. Similar to the research in [34], we minimize indirect and direct wait times in the model by determining appropriate schedule template for the providers. In another study by [33], a weekly schedule pattern in outpatient clinics for an OBGYN specialty is determined while considering only direct wait time for the patients.
A two-stage stochastic integer programming approach is one of the most common methods in appointment scheduling ( [20], [30], [32], [33] and [31]). The authors in [20] presented a two-stage stochastic mixed-integer program considering uncertainty in a system for optimizing booking and appointment times with the objective of maximizing expected profit for outpatient clinics. The research work in [32] formulated the appointment scheduling in two steps: a two-stage stochastic linear program for static appointment scheduling capturing no-show behavior; and subsequently, a dynamic scheduling using a multi-stage stochastic linear program. In research by [30], the stochastic overbooking in outpatient appointment scheduling for clinical use was modeled. Patient wait time, provider over time and patient revenue were considered in the objective function of the model formulation.
In this research, we devise a two-stage stochastic mixedinteger program for appointment scheduling while considering uncertainties in patient's demand and no-show behavior of patients. The model also considers fair and equitable and non-appointment slot assignments for each provider. We consider an OBGYN clinic for our case-study, and different types of services and categories are as shown in Table 1.

III. NOTATION AND MODEL FORMULATION
We design an appointment scheduling model capturing multitype patient channeling to different provider levels in an OBGYN clinic specialty. Based on the research on OBGYN specialty clinics by [33] and [35], we divide patients into three categories, and consequently, seven patient types.
Let N , R and T be the sets (indexed by n, r, t) denoting patient categories, providers and planning horizon, respectively. The set H represents the total available slots for a provider within the planning horizon, and N is the set of new patient categories which should have fair and equitable assignments to the providers. Each day in the planning horizon is divided into 16 fifteen-minutes slots. Service time duration for each patient type is based on the literature on OBGYN clinic ( [33], [35]). The authors in [35] collected data from the West Little Rock (WLR) clinic operated under the University of Arkansas for Medical Sciences (UAM). The parameters a r and ρ r denote the target for new patient category and penalty for the deviation from the target for a provider r. The average service time for a new patient is defined as penalty. Parameters c i and F represent the criticality for a patient category i and overall criticality factor allowed in a slot, respectively, and M is a large constant. The parameter v r represents the total number of available slots for a provider r within the entire planning horizon H . Letω be a multi-variate random variable representing the uncertainties in call-volume and no-show behaviors of the patients.
Two-stage stochastic programming is a common approach for modeling decision making problems that involve uncertainty. First-stage decision variables represent 'here-andnow' decisions that are determined before the realization of randomness, and second-stage decisions are determined after scenarios representing uncertainties are presented. In this two-stage formulation, the first-stage decides the scheduling template for the providers, which indicates the number of patient categories to be assigned for each slot. The firststage scheduling template is exposed to multiple scenarios of patient call volumes in the second-stage where the indirect wait time is minimized. We minimize the deviation of newpatient type of each provider from the given target in the firststage, and expected indirect wait time in the second-stage. In the first-stage, the integer variable x irt denotes the number of patient category i assigned to provider r per time-slot t, variable e r measures the deviation of the assignments for a provider r from its target, and z rt is a binary variable indicating whether a slot has been used by the provider r or not.
In the second-stage, let P be the set of patients, indexed by p, and p(i) and p(d) represent patient's category and desired slot of appointment, respectively. The parameter w i represents a penalty weight for the patient category i ∈ N , and is a parameter in the range (0,1). The no-show parameter for a patient p(i) for a scenario ω is represented as a binary vector s p(i) (ω), where a value of one represents that the patient will not show up for the appointment, and zero for otherwise. The second-stage binary variable y p(i)rt (ω) represents the assignment of patient p(i) with a provider r for the time slot t in a scenario ω. The first-stage objective function minimizes the deviation from the proportion of new patients allotted to each provider, and second-stage objective function minimizes the expected indirect wait times for the patients.
The first-stage formulation is represented as follows: For a given first-stage solution x and a realization ω of the random variableω, the second-stage recourse function is defined as follows: t∈T :p(d)≤t r∈R In the above formulation, constraints (2) and (3) capture the deviation for total number of appointments for the new patients for each provider r. The lower or upper deviation from the target a r is captured using the variable e r and is penalized in the objective function (1). Constraint (4) is a capacity restriction which limits the number of different patient categories assigned in a slot, and constraint (5) defines the indicating relationship between the variables x and z.
Considering other non-clinic activities, the total availability for a provider is defined in the constraint (6). Finally, restrictions for the variables are defined in the constraints (7). In the second-stage formulation, constraints (9) state that individual assignment of patients depends upon the scheduling template generated in the first-stage. Constraints (10) assure that a patient gets an appointment within the planning horizon. For the objective, a super-linear function (t − p(d)) (1+ ) is used so that indirect delays among the patients are fair and equitable, and the function is weighted by the parameter w i based on the patient category i. The secondstage objective function (8) represents a weighted indirect wait time for the patients for a given scenario ω.

A. CLINIC PATIENT FLOW SIMULATION
The patient flow simulation for the clinic is executed for each day based on the scheduling template proposed by the optimization model described in the previous section. Patients' direct and indirect wait times are measured. The clinic patient flow involves two stages: time with the nurse and time with the provider. Once a patient walks into the clinic, the patient waits in the lobby for the nurse to become available. After being visited by the nurse, the patient waits for the provider in the exam room. We assume that patients who are scheduled for a particular day must be served before the end of that day even if the provider has to work overtime. Patients are assumed not to leave the exam room until the provider finishes all required tasks. Also, patients may arrive late for their appointments. By arriving late, patients may increase the wait times for the patients that follow them. Therefore, we assume that patients are called in the order of their arrival time.
The first-stage solution x is evaluated in a simulation model for a clinic, and the average direct wait time for the patients in each slot is noted. Whenever the direct wait time for a slot is beyond a given user threshold of twenty minutes, the assignment sequence for the x variables is prohibited by adding the constraints (12)-(15) for the next iteration of branch and cut procedure. Let L be the set of integers {1, 2, 3, 4}, and a parameter m l ∈ H . Similarly, G is the set of consecutive time-periods where the direct wait time from the simulation is beyond the user's threshold. The binary variablex irtl is used to discretize the variable x. The constraints added for the direct wait time are as follows: i∈N ,l∈L,t⊆Gx x irtl ∈ {0, 1}.
Constraints (12) -(13) discretize the variable x to identify the exact combination of number of patients in different types that are scheduled for a set of consecutive time slots G. Constraints (14) denote that such combination should be avoided in the next iteration of branch and cut algorithm. It should be noted that the constraints (12)-(15) are dynamically added, i.e, every feasible x is evaluated in a clinic simulation model during the branch and cut procedure, and whenever the threshold for direct wait time is violated, the respective constraints are added. For every combination (denoted by set G) that exceeds the given threshold, the corresponding set of constraints (12) - (14) are generated and added to the model for the next iteration of branch and cut procedure. Algorithm 1 presents pseudocode for our approach based on the two-stage optimization model and simulation. The two-stage stochastic programming model without direct wait time is solved, and a scheduling template is obtained from the first-stage solution. Using the scheduling template, each day in the planning horizon is simulated using given input with 200 replications, and specific percentiles of patient flow measures are calculated. If any of the measures violate the patients' direct wait time thresholds, a new set of constraints (12), (13), and (14) are added to the model, which is then re-optimized. The process is repeated until the optimization Algorithm 1 Clinic Patient Flow Simulation Input: Output of the two-stage stochastic model (2) -(11) Output: Constraints (12)- (14) Initialization : 1: Optimize the two-stage stochastic model (2) -(11); 2: for t in T do 3: Run simulation model for each day in the planning horizon; 4: if (estimates ≥ any of the thresholds) then 5: construct the corresponding constraints (12)-(14); 6: end if 7: end for 8: return constraints (12)-(14) model reaches a state where the indirect wait time is minimized without violating any patient flow thresholds for direct wait time.

B. CALL CENTER SIMULATION
To evaluate the efficiency of the two-stage scheduling template, we simulate the practice at the appointment call center. The call center simulation uses either ''first-call-first-service scheduling policy'' (FCFS) or the scheduling template's allocation from the two-stage stochastic programming model described in the earlier section. A FCFS is used to allot or cancel patients' appointments. The scheduler at the call center estimates the priority of available appointment slots that should be offered to a patient based on her desired date and patient type. When a patient requests an appointment for a desired date and slot, the policy allocates a slot based on its feasibility and proximity to the desired date. For feasibility check, the current total criticality factor is calculated with the new assignment, and it should be less than the allowed threshold. Similarly, when the scheduling template is used, an assignment is made based on the capacity allocated for each patient type in the template. Feasibility check is not required for scheduling template as it is accounted in the two-stage stochastic programming model via constraint (4). To generate a patient's choice regarding accepting an appointment, a random number from the uniform distribution U(0, 1) is generated and compared to an ''acceptance'' threshold. If the random number is lesser than this threshold, the patient accepts the corresponding appointment slot; otherwise, another slot is offered and the process continues until the patient accepts. Appointment cancellations are handled similarly: a random number is generated from U(0, 1), and if the random number is less than the clinic's cancellation rate, the patient or clinic cancels the appointment and the patient is removed from the scheduling grid.

IV. CASE STUDY
In this section, we present a case study to demonstrate the performance of the proposed two-stage stochastic programming and simulation models. We consider a multi-category outpatient appointment scheduling for a women's OBGYN clinic. The characteristics and demand for different patient types for the clinic used in the case study are acquired from the literature of women's OBGYN specialty clinics.

A. DATA AND STUDY DESIGN
Based on the literature in OBGYN clinics, the common issue is related to time, equipment, and exam rooms scheduled for several service categories. Each patient type needs specific services such as prenatal and follow up care of pregnancy and miscarriage. The study in [33] has divided patient types with respect to the needed service type into three categories ( Table 1). In this case study, each clinic session is defined as a day and is divided into 16 time-slots with identical service time duration of 15 minutes. There are two providers available on all days of the week who can provide all service categories for different patient types. Patients are scheduled with any available provider in each clinic session (morning/ afternoon). Service time duration for different servers in the clinic such as time spent with nurse and provider are based on [33] and [35].
We assume there are two sessions: morning and afternoon, and each session has eight time-slots. OBGYN clinics services rendered for different patient types are considered in different sessions of a week day based on the study in [33]. Also, typically women clinics consider appointment scheduling for all providers in a clinic and not for specific ones. Therefore, patients can be seen by any available provider upon their appointment time. Since multiple providers are assigned for each day, overbooking is allowed for each time slot. The parameters used for the study are shown in Table 2. Uniform distribution was used to obtain the service time for different patient categories. The seven patient types are organized under three categories, their respective distributions for the service times are provided in the table. Similarly, a different no-show rate was used different patient types and this was used to populate the no-show rate parameter s p(i) (ω) in the optimization model. The total number of patients were divided into three categories of patients in the proportion of 30%, 30%, and 40% for Low Risk OB, High Risk OB, and Gynecology, respectively. The segregation of categories is used to determine the number of patient types within each category. Similarly, the criticality factors for the patient types are 1.669, 0.39, and 0.67 for the first provider, and 1.2, 0.869, and 0.669 for the second provider. The total criticality factor for any time slot used for the computational experiments is 1.99. Based on the capacity, the total number of calls is 350 for the study. Also, for the study 'Low Risk OB', patient category was categorized as new patient type. The desired day for an appointment by the incoming calls of the patients were in proportion of 20% each for Monday, Tuesday and Thursday, 30% for Wednesday, and finally, 10% for Friday.
The assumption is that the most of the appointment requests fall during the middle of the week and it is substantially lower for the last day of the week. For the distribution of the calls, a decision is made for the day of the call based on the proportion. Then a uniform distribution U(1, 16) is used to assign it to a slot within a day. We used six scenarios for the second-stage, and a standard deviation of three slots U(−3, 3) around the average is used as distribution for creating the scenarios. The task of the first-stage model is to determine the scheduling template for the providers while minimizing the deviation and expected costs from first-and secondstages, respectively.

B. OPTIMAL WEEKLY SCHEDULING TEMPLATE
The output of two-stage stochastic model presented in the earlier section is a scheduling template (SC-T). The scheduling template is used as a guideline in the call center to assign appointments for incoming calls from the patients. The assignment pattern using scheduling template is compared with the 'FCFS' policy where the assignments are made solely based on the capacity of a desired slot of the patient. The SC-T is derived based on the costs of the secondstage scenarios where the potential call pattern is taken into account. However, FCFS policy is short sighted where it looks only the immediate requirement. Sample weekly scheduling templates are represented in the Tables 3 and 4. In the tables, NL, FL, FH, NG, MG, EG, and RG represent New Low-Risk OB, Follow Up Low-Risk OB, Follow Up High-Risk OB, New GYN, MAU GYN, Established GYN, and Results GYN, respectively. When the direct wait time is lenient, there are more patients scheduled for a given time slot. This is depicted in the tables, where Table 4 is more lenient in direct wait time for the patients in clinic compared to the threshold used for Table 3. During the call center simulation, appointments for the incoming calls from the patients are made based on the number and type of patient types assigned in a slot. If a particular type is not available then the nearest unused slot with the requested patient type is assigned. However with the 'FCFS' policy, the assignments are solely based on the capacity of a slot, i.e, whether by adding another patient will exceed the total criticality factor or not.

C. COMPUTATIONAL RESULTS AND INSIGHTS
The computational experiments were performed on a Dell computer with Xeon 2.60 GHz and 80 GB RAM. The algorithm was implemented using Java programming language, and CPLEX 12.10 is used as a solver. The two-stage model presented in the constraints (1)-(11) along with the dynamic direct wait time constraints (12)-(15) is referred as 'SMILP', and the scheduling template obtained by solving the model is referred as 'SC-T.' The deterministic equivalent problem (DEP) for the two-stage formulation was solved using CPLEX 12.10 as optimization solver. Dynamic direct wait time constraints are added to the first-stage model via 'Lazy-ConstraintCallback' feature of CPLEX. Whenever there is a first-stage feasible integer solutionx for SMILP, the callback function runs the clinic simulation withx. A threshold of 20 minutes is used for the study. Whenever the direct wait time in a slot is beyond 20 minutes in clinic simulation, the constraints (12)- (15) are generated for the slot and the one which is adjoining it, and added to the model. This way, the combination inx for the slot will be avoided in the next iteration. In a conventional comparison, the objective function values from the model and 'FCFS' policy are calculated and a value of stochastic solution is estimated. Instead of this, we simulated the clinic and call center operations, and used the policies to compare the utilization of providers, direct and indirect wait times. This provides a more realistic way to measure the performance of the policies. We conducted eight simulation runs, and each run was repeated five times. Figures 2 and 3 compare the performance of direct and indirect wait times between two-stage model and FCFS    time from SMILP compared to ten units for using FCFS. Similary, for the direct wait time, the SMILP model had insignificant wait time (less than a slot) compared to an average of 1.5 slots for the policy. This significant improvement is due to the dynamic addition of constraints (12)- (15) in the SMILP model. The number of such constraints added for each run is shown in Figure 4. In the figure, 'DW Cuts' refer to number of direct wait constraints. On an average, there were 115 such constraints added dynamically to the first-stage of SMILP model for each run.
In terms of distribution of indirect wait time for the patients in a replication, the averages were three and ten slots for SC-T and FCFS policy, respectively. Especially, it should be noted that only few patients had non-zero indirect wait time using the scheduling template compared to more than 70% using FCFS policy. Table 5 summarizes the results from the simulation study. The reported average values are for each patient. AIWT C-1 refers to average indirect wait time for patient category-1, similarly AIWT C-2 and C-3 refer to average indirect wait time for patient category-2 and category-3. The average direct wait time had the biggest benefit by using SC-T VOLUME 8, 2020 as the constraints were generated by the simulation model and added directly to the SMILP model.

V. CONCLUSION AND FUTURE DIRECTIONS
Outpatient access to the facilities is controlled through appointment scheduling. Patients often suffer from significant delays in obtaining appointments and wait time within the clinic. A well-established appointment scheduling system can help clinics reduce patients' indirect wait time while also improving patient flow within the clinic. Clinic managers have to handle multiple issues when scheduling patient appointments. While different patient types have different complexities, there are uncertainties in the pattern of calls for appointments and patients' willingness to wait for an appointment. Some patients may call in advance to book appointments, while others ask for same-day appointments. This study proposes solving the appointment scheduling problem with a two-stage stochastic programming model integrated with a simulation model to minimize patients' indirect wait time for an appointment while maintaining a patient flow in the clinic that is within the constraints of acceptable patient and practitioner performance thresholds. The model proposes a patient scheduling template for the call center in order to help the clinic manager to reduce the patient appointment delays, scheduling errors and improve the efficiency of resource allocation.
We use a numerical case study inspired by a real-world healthcare system to validate our proposed approach instead of a first-call-first-serve approach. Our model also performs better when the provider is assigned greater capacity for other clinical activities. Determining the optimal schedule for re-running the optimization model based on changes in the uncertainties is an important avenue for future research. Moreover, while our two-stage stochastic programming model determines the number of each patient type that can be scheduled in each appointment slot, it does not provide any guidance regarding when the call center should offer each open appointment to a patient. Integrating the call center scheduling process with our two-stage stochastic programming approach could result in further improving patients' indirect wait times. In the proposed risk-neutral twostage stochastic model, we consider expected value in the second-stage function. For a future extension, a mean-risk performance measure like Conditional-Value-at-Risk (CVaR) can be utilized to estimate the robustness of the scheduling template.
SAMIRA FAZEL ANVARYAZDI received the M.S. degrees in pure mathematics and mathematical statistics and the Ph.D. degree in industrial and systems engineering from Wayne State University. She is currently a Lecturer of data analytics with the Olin Business School and a Faculty Scholar with the Institute for Public Health, Washington University in Saint Louis. Her research interests include stochastic programming, robust optimization techniques, data science, and compartmental models with applications in healthcare systems engineering and business problems.
SARAVANAN VENKATACHALAM received the M.S. and Ph.D. degrees in industrial and systems engineering from Texas A&M University, College Station, in 2003 and 2014, respectively. He is currently an Assistant Professor of industrial and systems engineering with Wayne State University. His research interests include stochastic programming, large-scale optimization, and discrete event modeling and simulation, and applications of interest include supply chain management, healthcare, pricing and revenue management, and energy management.
RATNA BABU CHINNAM received the M.S. and Ph.D. degrees in industrial engineering from Texas Tech University, Lubbock, in 1990 and 1994, respectively. He is currently a Professor of industrial and systems engineering and the Director of the Big Data and Business Analytics Group, Wayne State University. His research interests include business analytics, big data, supply chain management, freight logistics, operations management, sustainability, healthcare systems engineering, and smart engineering systems.