Dynamic Distributed Ambulatory Care Scheduling

Abstract

We investigate an ambulatory care scheduling problem derived from a real case in Ontario, Canada that offers multi-appointment, multi-class, multi-priority treatments in geographically distributed campuses with multiple resources. We consider a dynamic setting with uncertain patient arrival and use of the emergency department. This problem is formulated as an infinite-horizon Markov decision process model. Since we cannot solve large-sized instances via conventional approaches, we hybridize this model with a neural network to simplify feasibility constraints while respecting all assumptions. Given the curse of dimensionality, we use an affine approximation architecture to estimate the value function. An equivalent linear programing model is solved through column generation in order to compute approximate optimal policies and derive two easy-to-implement scheduling policies. Simulation results demonstrate that the approximate optimal policy and heuristics outperform alternative scheduling policies. Finally, we demonstrate that the application of our methodology can enhance performance metrics in a large ambulatory care center in Canada. We show that a template-based scheduling rule can result in high resource utilization but poor scheduling decisions. However, an efficient scheduling policy equips a booking clerk with intelligent scheduling rules that are difficult for her to predict in real-time and work well in comparison to scheduling templates.

1. Introduction

Healthcare providers are under pressure to enhance service quality and reduce costs (Papanicolas et al., 2018). Given the greater emphasis on preventive medicine, outpatient services have become an increasingly important element in our health systems (Cayirli and Veral, 2003). Delays in outpatient services are the consequence of insufficient capacity as well as inefficient scheduling policies (Sauré et al., 2012). Thus, the optimal scheduling of outpatients is crucial to decrease wait times and increase the utilization of limited and expensive personnel.

Ambulatory care (AC) denotes health services offered on an outpatient basis (Comerford and Shah, 2019), including diagnosis, observation, consultation, treatment, and rehabilitation services. AC services aim to reduce avoidable inpatient stays, equipping patients with the knowledge to manage their treatment backed by a multidisciplinary team (Cooper and de Lord, 2018). AC is one of the largest-volume health services in Canada, making it a significant component of the health system (CIHI, 2022).

Historically, complex treatments were administered in an inpatient setting due to their uncertainties (Cooper and de Lord, 2018). This put healthcare providers in unviable positions due to increasing expenditures and ever-rising demand. Thus, AC has become the primary provider of various treatments in recent years (e.g., chemotherapy, infusion therapy) (Comerford and Shah, 2019). Despite these advancements, wait times remain excessive, a direct consequence of earlier diagnoses, an aging population, advances in diagnostics, increased survival rate and more effective care procedures (Cooper and de Lord, 2018; Fung-Kee-Fung et al., 2022). Delayed treatments can have a substantial negative impact on outcomes (Fung-Kee-Fung et al., 2022). Long wait times are one of the main sources of patient (and staff) dissatisfaction, depression and anxiety (Frick et al., 2007). Apart from affecting patient and staff satisfaction, long wait times may adversely affect patient adherence to scheduled appointments (O’Neill et al., 2012).

In this paper, we investigate an AC scheduling problem for a hospital with geographically distributed campuses, where all new arrivals are added to a unique waiting list. Patients on the waiting list may end up either using the Emergency Department (ED) or transferring to another AC center due to long wait times (or medical abnormalities). The daily challenge facing the booking clerk is to assign the available capacity between various patient classes in order to decrease the number of patients whose wait time exceeds a preset threshold with greater importance given to any late bookings of higher-priority demand (Patrick et al., 2008). Thus, our research objective is to design a decision support tool to schedule patients in a distributed center with patients’ arrival and use of the ED being uncertain. While we develop our problem assumptions based on a partnership with a hematology and oncology AC center in Canada, our methodology remains practical for various other AC scheduling applications, simply by relaxing some of our assumptions.

The remainder of this paper is organized as follows. Section 2 presents a review of the literature and specifies the contributions. Section 3 provides the problem description and introduces the chief components of a dynamic program for our distributed AC scheduling problem, and Section 4 develops an approximate dynamic programing (ADP) algorithm to solve large-sized instances of this problem. Section 5 offers computational results and managerial implications. Finally, Section 6 provides conclusions and future research directions.

2. Literature Review

Outpatient scheduling problems have received considerable attention from scholars and practitioners (see Ahmadi-Javid et al. (2017) and Cayirli and Veral (2003) for broad literature reviews). Typically, outpatient scheduling problems are divided into two levels: (i) Advance scheduling (allocating capacity to future demand), and (ii) appointment scheduling (specifying patients’ appointment times) (Pan et al., 2020). In this section, we first review the literature on outpatient scheduling using ADP algorithms. Second, we review chemotherapy appointment scheduling papers. Finally, we discuss the main contributions of our research.

2.1. Dynamic Programming for Outpatient Scheduling

The exponential growth of complexity in dynamic programs is known as the “curse of dimensionality,” one of the main downsides of dynamic programing. ADP is a growing research area that aims to overcome this obstacle primarily through three streams: (i) Linear programing-based algorithms, (ii) simulation-based algorithms, and (iii) aggregation-based algorithms. Linear programing-based ADP algorithms have attracted the most attention in the AC scheduling literature. Patrick et al. (2008) studied a multi-priority advance scheduling problem for a diagnostic facility in a public healthcare setting. They proposed a Markov decision process (MDP) model that minimizes wait time, deferral and diversion costs. The authors developed a linear programing-based ADP algorithm with linear basis functions and solved it through column generation (CG). Sauré et al. (2012) addressed a multi-priority multi-appointment advance scheduling problem where patients can receive treatment across multiple days and for irregular lengths of time. The authors formulated this problem as a discounted infinite-horizon MDP and developed an ADP algorithm using an equivalent linear programing formulation and an affine architecture to deal with the intractable numbers of variables and constraints.

Zhou et al. (2022) studied a multi-class advance scheduling problem for a diagnostic facility considering heterogeneous wait time targets and access equity. They developed a finite-horizon MDP model which minimizes the total expected costs. To deal with the curse of dimensionality, the authors reformulated the MDP as a multi-stage stochastic programing model and proposed a modified Benders decomposition algorithm based on new dual integer cuts. Diamant et al. (2018) studied an advance scheduling problem where patients require multiple appointments without any pre-specified order. Due to patient no-shows, Diamant et al. considered overbooking to ensure full capacity utilization. The authors proposed an infinite-horizon MDP that maximizes the discounted profit and solved it through an equivalent linear programing model with the value function approximation and decision variable aggregation. Likewise, Diamant (2021); Göçgün (2018a); Göçgün and Puterman (2014); Parizi and Ghate (2016); Sauré et al. (2020); Wang et al. (2018); Wang and Fung (2015a, 2015b) developed linear programing-based ADP algorithms to address outpatient scheduling problems.

Another stream in the ADP literature is simulation-based algorithms. Sauré et al. (2015) proposed a simulation-based ADP algorithm for advance patient scheduling. This algorithm was based on policy iteration and a post-decision state formulation, and used a non-linear logistic value function approximation. Sauré et al. demonstrated that policies obtained through their proposed algorithm outperform other scheduling policies in most cases. Göçgün (2018b) studied a multi-appointment advance scheduling problem for radiation therapy with stochastic arrivals and appointment cancelations. The author formulated an MDP model and deployed a simulation-based ADP algorithm that relies on least-squares-based approximate policy iteration. Computational results demonstrated better performance of the proposed ADP algorithm in comparison to the myopic heuristic decision rule. Schuetz and Kolisch (2012) considered a multi-class advance scheduling problem for service industries where the service company can either confirm or reject arriving demand. They let the service duration be stochastic and allowed for no-shows and cancelations. Schuetz and Kolisch designed a simulation-based ADP algorithm based on a post-decision state formulation and showed the efficiency of their methodology regarding the objective function, solution time and memory requirement.

Wang and Fung (2015a) addressed a multi-class sequential appointment scheduling problem while respecting patient preferences for the appointment time and solved it using a simulation-based ADP algorithm. The authors demonstrated that a well-designed appointment system can enhance the resource utilization and the revenue of outpatient departments. Lu et al. (2018) studied an advance scheduling problem where all requests receive estimated appointment dates upon their arrival. The authors formulated an MDP model and solved it through a simulation-based ADP algorithm that aggregates multiple days into a single period. This research showed that informing patients of their appointment dates in advance is a cheap yet effective solution to improve their experiences. Similarly, Geng et al. (2011); Li et al. (2018); Lin et al. (2011) deployed simulation-based ADP algorithms to solve outpatient scheduling problems.

Aggregation-based algorithms are the last stream in the ADP literature for AC scheduling problems. Wang et al. (2018) studied an appointment scheduling problem considering patient preference for the appointment time. While appointment requests arrive sequentially within a booking horizon, patients choose between a number of time slots (specified based on the availability of physicians). The authors proposed a dynamic program that maximizes the expected revenue for a given service day. They applied an aggregation-based approximation method which aggregates several time slots into a single slot and used a complete-set policy to decrease the number of system states as well as decisions. Li et al. (2018) also investigated a sequential appointment scheduling problem where patient choices for physicians and time slots are considered. They proposed an MDP model that maximizes the expected patient satisfaction and solved the problem through an aggregation-based ADP algorithm. Li et al. found that different aggregation levels lead to different results and computation times. Gedik et al. (2017) investigated a patient admission planning problem for a proton therapy facility in order to maximize the total expected number of treatment sessions delivered to patients in a planning period with stochastic patient arrivals. The problem was formulated as an MDP model that was simplified using a state aggregation technique developed based on the fixed-weight method.

2.2. Hematology and Oncology Appointment Scheduling

A hematology and oncology AC center offers various treatments, including—but not limited to—intravenous immunoglobulin, chemotherapy, iron infusion, paracentesis, phlebotomy, blood transfusion, as well as some diagnostic procedures. Patients of these clinics typically go through an initial consultation session, followed by some regular treatment sessions. Haghi et al. (2023) investigated the simultaneous scheduling of consultation and treatment sessions for different patient classes. The authors evaluated decisions regarding the assignment and sequencing of patients for nurses and beds. They assumed that treatment duration is uncertain and consists of setup and infusion times, where each nurse can only set up one patient at a time and can monitor several patients during the infusion. According to Garaix et al. (2020), the daily care in a hematology and oncology AC center includes (i) consultation, (ii) drug preparation, and (iii) injection. In this paper, the authors addressed all three decisions simultaneously, assuming a series of fixed injection plans once patients are booked for their first injection appointment. They allow nurses to work overtime with a penalty while injection duration uncertainty and patient punctuality.

On the other hand, some papers focus only on the drug injection scheduling decisions. Gul (2023) assumed that each patient has a primary nurse who should be served by the same nurse throughout their course of treatment, which can be violated with a penalty. Their approach decides on patient sequencing, appointment time setting, and patient assignment to nurses and beds. Their model minimizes patient clinic wait time, nurse overtime, and violations from the primary nurse assignment under infusion duration uncertainty. Demir et al. (2021) examined chemotherapy scheduling, taking into account the availability of nurses and beds amid service duration uncertainties. The authors considered a two-component service duration (premedication and infusion). Nurses can handle only one patient for premedication, however, they can proctor the infusion process of multiple patients. Their method sequences patients on a daily appointment list, determines appointment times, and assigns patients to nurses and chairs in order to minimize a weighted combination of patient waiting time, chair idle time, and nurse overtime. Similarly, Benzaid et al. (2020); Karakaya et al. (2023) investigated the drug injection scheduling decisions.

Another stream of studies is the use of scheduling templates for hematology and oncology AC appointment scheduling. Given a fixed scheduling template for a long planning horizon, Hesaraki et al. (2023) focused on the online, daily appointment scheduling for a hematology and oncology AC center. The main purpose of this paper is to minimally adjust the input scheduling template in order to accommodate new and returning patients, optimizing the interests of multiple stakeholders. Hahn-Goldberg et al. (2014) similarly focused on daily appointment scheduling for a hematology and oncology AC center. New appointment requests are scheduled according to a predefined scheduling template. If a request does not fit the template, their proposed approach updates the template using an optimization model. To handle last-minute additions and cancelations, the authors proposed a shuffling algorithm that adjusts appointment start times within a predefined time limit. Their findings show that dynamic template scheduling can improve makespan by up to 20% compared to current practices. Likewise, Condotta and Shakhlevich (2014); Faridimehr et al. (2021); Huang et al. (2019) studied the scheduling templates for hematology and oncology AC patients.

2.3. Main Contributions

Our main contributions are four-fold. To the best of our knowledge, our paper is the first to develop an ADP algorithm for a distributed multi-appointment, multi-class, multi-priority, multi-resource AC scheduling problem. Not only do we assume heterogeneity of resources (beds and nurses) and different resource requirements of patient classes, but we also consider that service durations associated with beds and nurses have different settings (i.e., a bed is fully busy during a patient appointment while a nurse may not be). Second, nearly all the literature focuses on either advance scheduling problems or appointment scheduling problems because handling both at the same time is complex (Sauré et al., 2020). Considering both decision levels requires binary variables to formulate the problem, introducing computational challenges (i.e., we cannot solve large-scale instances of the problem within a reasonable time via conventional approaches). To address this, we propose a hybrid methodology that combines ADP with a neural network (NN). This approach effectively incorporates patient allocation and sequencing considerations, enabling us to solve large-scale advance scheduling problems that are more constrained than those typically addressed by infinite-horizon ADP algorithms in the literature. Moreover, the methodology is versatile and can be applied broadly to other optimization problems under complex constraints (e.g., machine scheduling, storage location assignment, crew scheduling, facility layout), where a component of the problem can be partitioned in a similar fashion. Third, we assume multi-component service durations as described in Section 3. Such an assumption is rare in the literature, and to the best of our knowledge, our research is the first to consider three unique components associated with hematology and oncology treatments in the service duration. This assumption is applicable to a wide range of other problems, such as (i) operating room scheduling where surgical procedures have pre- and post-operative phases, (ii) manufacturing scheduling problems with setup, processing and breakdown phases, (iii) home healthcare visits where some procedures involve initial setup by a nurse when they can perform other tasks, (iv) logistics problems with the notion of loading and unloading in freight or transportation. Fourth, we propose two easy-to-implement heuristics for our problem that significantly outperform benchmark policies and perform well in comparison to our primary solution approach. The policies provide generalizable managerial insights that could benefit similar scheduling problems.

3. Dynamic Programming

This section describes all the assumptions and objectives of the problem under consideration and, then formulates the MDP model.

3.1. Problem Description

We focus on a distributed hematology and oncology AC (i.e., the center offers services on multiple campuses), that aims to avoid hospital admissions during the course of treatment by providing timely care. We investigate a setting with nearby campuses and assume patients are indifferent to the campus they are allocated to, as long as they receive all their treatment at the same campus (continuity of care). Each campus in the hematology and oncology AC has a different but known availability of human and physical resources. All campuses have the resources (nurses and beds) to serve benign treatments, and each is also equipped to serve a specific group of malignant treatments. Each campus has fixed operating hours with no planned overtime, meaning our planned schedule includes no built-in overtime. However, if the solutions are implemented, some overtime may still occur due to uncertainties involved in practice.

Demand is prioritized (on a scale from one to five—the lower the priority, the higher the urgency) and may require a single- or multi-appointment treatment (each treatment on a particular day) with a specific treatment frequency. Each patient is associated with a wait time threshold considering her treatment class (or type) and priority level. The wait time is an important performance metric for hematology and oncology AC centers given its impact on health outcomes and patient satisfaction. A deferral cost occurs each time a patient on the waiting list is left unscheduled. The deferral cost ensures that scheduling decisions are not delayed unnecessarily. For example, if a patient arrives on day 25 and has a three-day wait time threshold, scheduling her for day 28 on either day 26 or day 27 will incur no wait time cost. However, making scheduling decisions earlier is preferable and is incentivized by adding the deferral cost. Patients may utilize the ED due to their health condition, which is costly for the health center. Similar to the wait time, each class-priority is associated with a deferral threshold. If a patient is not scheduled within the deferral threshold, she will be transferred to another hematology and oncology AC center (or rejected), incurring a large cost. Therefore, the objective function in this study minimizes a weighted sum of wait time, deferral, ED use, and transfer costs.

Each patient has a ready-to-treat date or earliest permissible day for scheduling, which is determined based on factors such as readiness of test results and personal preparation. Once the first session of a multi-appointment treatment is scheduled, the follow-up appointments are automatically scheduled based on the treatment frequency. Treatments in our hematology and oncology AC may have different service durations but each includes three components: Initialization, passivity and finalization. In the initialization period, a nurse serves the patient and cannot do any other activity (e.g., the nurse prepares the equipment and attaches the blood transfusion set to the patient to start the treatment). The passivity period follows the initialization period and does not require the full involvement of the nurse (e.g., the time blood transfusion takes). In the finalization period, the same nurse is again busy completing the treatment (e.g., the nurse detaches the blood transfusion set, discharges the patient and cleans the bed). The patient occupies a bed during her entire treatment. Note that the nurse can handle multiple patients at the same time, limited to a maximum, but she cannot serve the initialization and/or finalization periods of two patients at the same time.

Finally, we consider two sources of uncertainty. First, we assume that patient demand is uncertain and follows a truncated Poisson distribution for each class-priority. New patients are added to the waiting list until they are scheduled or transferred to another hematology and oncology AC center. Second, we consider the uncertain use of the ED, which follows a truncated Poisson distribution for each class-priority. For patients on the waiting list, the probability of using the ED increases by their priority and wait time. For patients who are already scheduled, whether they have started treatment or not, we assume that this probability increases according to their priority. Patients on the waiting list or already scheduled may have to use the ED due to their health condition. The higher the priority level and the larger the wait time, the higher the chance of using the ED.

3.2. Markov Decision Process Model

This sub-section defines the components of an infinite-horizon MDP model for the advance scheduling problem.

- Indices:
$c$	Index of campuses.
$b$	Index of beds.
$k$	Index of patient classes.
$o$	Index of priority levels.
$i,i^{\prime }$	Indices of the remaining appointments.
$e$	Index of the number of deferrals.
$t,t^{\prime }$	Indices of days.
- Sets:
$\mathcal{C}$	Set of all campuses ( $\{1,2,\dots ,C\}$ ).
$\mathcal{B}$	Set of all beds ( $\{1,2,\dots ,B\}$ ).
${\mathcal{B}}_c$	Set of beds available in campus $c$ ( $\{1,2,\dots ,B_c\}$ ).
$\mathcal{K}$	Set of all patient classes ( $\{1,2,\dots ,K\}$ ).
$\mathcal{O}$	Set of all priority levels ( $\{1,2,\dots ,O\}$ ).
${\mathcal{I}}_k$	Set of all possible remaining appointments for class $k$ ( $\{1,2,\dots ,I_k\}$ ).
${\mathcal{E}}_k$	Set of all possible deferrals for class $k$ ( $\{1,2,\dots ,E_k\}$ ).
$\mathcal{T}$	Set of all decision days ( $\{1,2,\dots ,T\}$ ).
${\mathcal{T}}^{\prime }$	Set of all planning days ( $\{1,2,\dots ,T,T+1,\dots ,T^{\prime }\}$ ).
- Parameters:
$S_k^1$ , $S_k^2$ , $S_k^3$	Service duration of initialization, passivity and finalization periods for class $k$ , respectively.
$G_{kb}$	If patients with class $k$ can be served on bed $b$ , 1; otherwise, 0.
$V_{kot}^1$	The discounted wait time cost of scheduling a patient on the waiting list with class $k$ and priority $o$ on day $t$ .
$V_{koe}^2$	The deferral cost of a patient on the waiting list to the next decision epoch with class $k$ , priority $o$ and having waited $e$ days.
$V_k^3$	The cost of using the emergency department (ED) for patient class $k$ .
$V_{koi}^4$	The transfer cost of a patient on the waiting list with class $k$ , priority $o$ and $i$ remaining appointments.
$F{R}_k$	The frequency of class $k$ .
$A{V}_{ko}$	The earliest permissible day for scheduling a patient with class $k$ and priority $o$ .
- Decision variables:
$z_{koietc}^1$	The number of patients with class $k$ , priority $o$ , $i$ remaining appointments and having waited $e$ days that are allocated to campus $c$ on day $t$ for their first appointment.
$z_{ktc}^2$	The number of patients with class $k$ that are allocated to campus $c$ on day $t$ for their follow-up appointment(s).

3.2.1. Decision Epochs, Decision Horizon and Planning Horizon

The decision epochs correspond to the end of each day. At each decision epoch, we make a decision for the next $T$ days (a $T$ -day decision horizon). We extend the planning horizon to $T^{\prime }$ days ( $T^{\prime }\ge T$ ), ensuring it is large enough that patients scheduled within the decision horizon can complete their treatments. Note that we investigate a rolling planning horizon, where day $t$ of the current planning horizon is associated with day $t-1$ of the next planning horizon.

3.2.2. State Variable

The state variable takes the form

$$\begin{array}{rl}\overrightarrow{s}& =(p_{1111}^w,p_{1112}^w,\dots ,p_{koie}^w,\dots ,p_{KO{I}_K{E}_k}^w;\\ & p_{111}^m,p_{112}^m,\dots ,p_{ktc}^m,\dots ,p_{K{T}^{\prime }C}^m)\end{array}$$

where $p_{koie}^w$ refers to the number of patients on the waiting list with class $k$ , priority $o$ , $i$ remaining appointments and $e$ days wait time, and $p_{ktc}^m$ denotes the number of patients allocated to campus $c$ on day $t$ with class $k$ . Note that patients are transferred to another hospital if they are not scheduled after $E_k$ days. The state space is restricted to a $(\sum _{k\in \mathcal{K}}[I_k\times O\times E_k]+K\times T^{\prime }\times C)-$ dimension vector.

3.2.3. Decision Variables

The primary decisions are the allocation of patients to campuses such that the total costs are minimized over the decision horizon. Thus, $z_{koietc}^1\in {\mathbb{Z}}_{\ge 0}$ determines the number of patients from class $k$ , priority $o$ , with $i$ remaining appointments and $e$ days wait time that are allocated to campus $c$ on day $t$ for their first appointment. And $z_{ktc}^2\in {\mathbb{Z}}_{\ge 0}$ specifies the number of patients from class $k$ that are allocated to campus $c$ on day $t$ for their follow-up appointments. For simplicity reasons, we use $\overrightarrow{z}$ and ${\mathrm{\Xi }}_{\overrightarrow{s}}$ to refer to the decision vector and all possible values of $\overrightarrow{z}$ for state variable $\overrightarrow{s}$ , respectively.

3.2.4. Transition Probabilities

Two sources of uncertainty occur at each decision epoch: (i) daily arrival of new patients, and (ii) use of the ED. The former source of uncertainty does not depend on the state of the system while the latter does. The state of the system for the next decision epoch, denoted by ${\overrightarrow{s}}^{\prime }=({\overrightarrow{p}}^{{}^{\prime }w},{\overrightarrow{p}}^{{}^{\prime }m})$ , relies on three random variables: (i) ${\theta }_{koi}$ , the number of new patients from class $k$ , priority $o$ and with $i$ remaining appointments, (ii) ${\zeta }_{koie}$ , the number of patients on the waiting list with class $k$ , priority $o$ , $i$ remaining appointments and $e$ days wait time that used the ED, and (iii) ${\varphi }_{ktc}$ , the number of patients with class $k$ scheduled in campus $c$ on day $t$ that used the ED. The transition probabilities that determine ${\overrightarrow{s}}^{\prime }$ based on $\overrightarrow{s}$ and $\overrightarrow{z}$ can be formulated as below:

$$pr({\overrightarrow{s}}^{\prime }|\overrightarrow{s},\overrightarrow{z})=\{\begin{array}{ll}{\displaystyle \prod _{k\in \mathcal{K}}\prod _{o\in \mathcal{O}}\prod _{i\in {\mathcal{I}}_k}\prod _{e\in {\mathcal{E}}_k}\prod _{t\in {\mathcal{T}}^{\prime }}\prod _{c\in \mathcal{C}}\mathbb{P}({\theta }_{koi})\mathbb{P}({\zeta }_{koie})\mathbb{P}({\varphi }_{ktc})}& \text{If }{\overrightarrow{s}}^{\prime }\text{ satisfies Equations }(1)-(5),\\ 0& \text{Otherwise.}\end{array}$$

Transitions are governed by the following equations.

$$\begin{array}{rl}{p}_{koi0}^{{}^{\prime }w}& ={\theta }_{koi}\mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k\end{array}$$

(1)

$$\begin{array}{rl}{p}_{koie}^{{}^{\prime }w}& =p_{koi(e-1)}^w-\sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koi(e-1)tc}^1-{\zeta }_{koi(e-1)}\\ & +\sum _{\begin{array}{c}{i}^{\prime }\in {\mathcal{I}}_k\\ i^{\prime }=i+1\end{array}}{\zeta }_{ko{i}^{\prime }(e-1)}\\ & \mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k;e>0\end{array}$$

(2)

$$\begin{array}{rl}{\zeta }_{koie}& \le p_{koie}^w-\sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koietc}^1\\ & \mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k;e>0\end{array}$$

(3)

$$\begin{array}{rl}{p}_{ktc}^{{}^{\prime }m}& =p_{k(t+1)c}^m+\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{\begin{array}{c}{t}^{\prime }\in \mathcal{T};t^{\prime }=t\\ t^{\prime }\ge A{V}_{ko}\end{array}}\sum _{e\in {\mathcal{E}}_k}{z}_{koie(t^{\prime }+1)c}^1\\ & +\sum _{\begin{array}{c}{t}^{\prime }\in {\mathcal{T}}^{\prime };t^{\prime }=t\\ t^{\prime }\ge \underset{o\in \mathcal{O}}{min}A{V}_{ko}\end{array}}{z}_{k(t^{\prime }+1)c}^2-{\varphi }_{k(t+1)c}\\ & \mathrm{\forall }k\in \mathcal{K};t\in {\mathcal{T}}^{\prime }/T^{\prime };c\in \mathcal{C};\sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}$$

(4)

$$\begin{array}{rl}{\varphi }_{ktc}& \le p_{ktc}^m+\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{\begin{array}{c}{t}^{\prime }\in \mathcal{T};t^{\prime }=t\\ t^{\prime }\ge A{V}_{ko}\end{array}}\sum _{e\in {\mathcal{E}}_k}{z}_{koie{t}^{\prime }c}^1+\sum _{\begin{array}{c}{t}^{\prime }\in {\mathcal{T}}^{\prime };t^{\prime }=t\\ t^{\prime }\ge \underset{o\in \mathcal{O}}{min}A{V}_{ko}\end{array}}{z}_{k{t}^{\prime }c}^2\\ & \mathrm{\forall }k\in \mathcal{K};i\in {\mathcal{I}}_k;t\in {\mathcal{T}}^{\prime };c\in \mathcal{C};\sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}$$

(5)

$$\begin{array}{rl}& \sum _{k\in \mathcal{K}}\sum _{c\in \mathcal{C}}{p}_{k{T}^{\prime }c}^{{}^{\prime }m}=0\end{array}$$

(6)

New arrivals are added to the waiting list ( $e=0$ ) based on equation (1). Equation (2) ensures that the new waiting list associated with class $k$ , priority $o$ , $i$ remaining appointments, and $e$ days wait time ( $e>0$ ) consists of patients on the current waiting list (excluding recently scheduled ones and those with $i$ remaining appointments that used the ED) plus patients with $i+1$ remaining appointments that used the ED (visiting the ED reduces the number of appointments by one). Equation (4) ensures that the new list of scheduled patients on day $t$ and campus $c$ includes previously and recently scheduled patients (for both first and follow-up appointments) minus those that were scheduled on the same day and campus but used the ED. Equations (3) and (5) specify upper bounds on ${\zeta }_{koie}$ and ${\varphi }_{ktc}$ , respectively. Finally, equation (6) ensures no patient is scheduled on the last day of the next planning horizon (since it is rolling).

3.2.5. Feasible Decision Set

The set of feasible decisions compatible with state $(\overrightarrow{p^w},\overrightarrow{p^m})\in \overrightarrow{s}$ must satisfy the following constraint sets:

$$\begin{array}{rl}& \sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koietc}^1\le p_{koie}^w\\ & \mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k\end{array}$$

(7)

$$\begin{array}{rl}& z_{ktc}^2=\sum _{o\in \mathcal{O}}\sum _{\begin{array}{c}{t}^{\prime }\in \mathcal{T}\\ t^{\prime }\ge A{V}_{ko};t>t^{\prime }\\ (t-t^{\prime })(\mathrm{mod}F{R}_k)=0\end{array}}\sum _{\begin{array}{c}i\in {\mathcal{I}}_k\\ i>\frac{t-t^{\prime }}{F{R}_k}\end{array}}\sum _{e\in {\mathcal{E}}_k}{z}_{koie{t}^{\prime }c}^1\\ & \mathrm{\forall }k\in \mathcal{K};t\in {\mathcal{T}}^{\prime };c\in \mathcal{C};t\ge \underset{o\in \mathcal{O}}{min}A{V}_{ko};\sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}$$

(8)

$$\begin{array}{rlr}& F(p_{ktc}^m+\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}\sum _{\begin{array}{c}{t}^{\prime }\in \mathcal{T};t^{\prime }=t\\ t^{\prime }\ge A{V}_{ko}\end{array}}{z}_{koie{t}^{\prime }c}^1& \\ & +\sum _{\begin{array}{c}{t}^{\prime }\in {\mathcal{T}}^{\prime };t^{\prime }=t\\ t^{\prime }\ge \underset{o\in \mathcal{O}}{min}A{V}_{ko}\end{array}}{z}_{k{t}^{\prime }c}^2,\mathrm{\forall }k\in \mathcal{K})\le 0\mathrm{\forall }t\in {\mathcal{T}}^{\prime };c\in \mathcal{C}\end{array}$$

(9)

$$\begin{array}{rl}& z_{koietc}^1\in {\mathbb{Z}}_{\ge 0}\mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k;t\in \mathcal{T};\\ & c\in \mathcal{C};t\ge A{V}_{ko};\sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}$$

(10)

$$\begin{array}{rl}& z_{ktc}^2\in {\mathbb{Z}}_{\ge 0}\mathrm{\forall }k\in \mathcal{K};t\in {\mathcal{T}}^{\prime };c\in \mathcal{C};\\ & t\ge \underset{o\in \mathcal{O}}{min}A{V}_{ko};\sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}$$

(11)

Constraint set (7) restricts the total number of patients on the waiting list with class $k$ , priority $o$ , $i$ remaining appointments and $e$ days wait time that are allocated to campuses for their first appointment during the decision horizon to be less than or equal to $p_{koie}^w$ . Constraint set (8) ensures that follow-up appointments are scheduled on time in an appropriate campus. If at least one patient is scheduled on day $t^{\prime }$ in an eligible campus $c$ for the first appointment ( $z_{koie{t}^{\prime }c}^1>0$ and $\sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0$ ) and requires at least one follow-up appointment, this constraint set schedules the corresponding patients on an appropriate day $t$ on campus $c$ ( $t>t^{\prime }$ , $i>\frac{t-t^{\prime }}{F{R}_k}$ and $(t-t^{\prime })(\mathrm{mod}F{R}_k)=0$ ). Given a patient allocation, Constraint set (9) ensures that there exists a feasible patient schedule with no overtime (at the design phase) and that respects the limited number of beds and nurses. Here, a patient allocation includes both first and follow-up appointments. In this constraint set, $F(\cdot )$ returns a continuous value. If the feasible patient schedule has overtime equal to zero, $F(\cdot )\le 0$ ; otherwise, $F(\cdot )>0$ . Readers are referred to Constraint sets (EC.12) to (EC.18) in the E-companion for a detailed mathematical representation of function $F(\cdot )$ . Finally, Constraint sets (10) and (11) specify the valid domains of non-negative integer decision variables.

It is noteworthy to mention that we had initially modeled the MDP feasibility constraints via binary variables to make both patient allocation and scheduling decisions. The use of binary variables enables us to directly ensure the feasibility of solutions with regard to the limited availability of beds and nurses, using no overtime (instead of using function $F(\cdot )$ ). However, the MDP resulted in an unreasonably large model in terms of both constraints and variables. Based on our experiments, commercial solvers are unable to handle such models. Therefore, we reformulated the MDP via integer variables and introduced on function $F(\cdot )$ in Constraint set (9) in order to ensure that there exists a feasible patient schedule. Unlike conventional methods, our approach enables us to solve a complex, large-sized healthcare scheduling problem. Function $F(\cdot )$ is flexible enough that it can be adapted for a wide range of problem types, such as (i) operating room scheduling where surgical procedures have pre- and post-operative phases, and (ii) manufacturing scheduling problems with setup, processing and breakdown phases.

To build a predictive model for function $F(\cdot )$ , we could use logistic regression models for each campus-day to predict if there is overtime. However, Sauré et al. (2015) demonstrated a non-linear trend in a simpler scheduling problem, suggesting that logistic regression may not be ideal. Our experiments also indicated that simpler modeling choices, such as logistic regression, significantly underperformed NN (see Section 5.1.2).

We build function $F(\cdot )$ in three steps (Figure 1). First, we generate a data set of input and output values. Data set generation relies on a simulation model that randomly generates patient allocations (input values) and two heuristics that compute nearly minimum overtime values for the generated patient allocations while respecting the limited availability of beds and nurses (output values). Second, we specify the network architecture and train an NN using the input and output values. Third, the trained NN is linearized and used to build function $F(\cdot )$ .

Figure 1.

Illustration of steps for building function $F(\cdot )$ .

Data set generation: To be generalizable, NNs should be trained over large-sized and balanced data sets. We create a large data set of inputs and outputs through simulation. The simulation generates random patient allocations for each campus on a single day. For each patient allocation, heuristics are deployed to find feasible patient schedules such that overtime is minimized across all campuses. The simulation may include thousands of patient allocations and finding the minimum overtime value using mathematical programing for each scenario requires significant CPU time. Therefore, we use two different heuristics to increase robustness in computing high-quality solutions. Algorithm 1 outlines this simulation model.

In Heuristic 1, for each campus, we sort eligible patient classes based on the number of beds available (the fewer beds available, the higher the priority). Then, we schedule all allocated patients of each class to feasible “bed-nurse” combinations with the earliest service start times. Here, a “bed-nurse” for a patient is a feasible combination when (i) the bed is appropriate for the patient and is free for her entire service duration, (ii) the nurse is not busy with the initialization or finalization period(s) of other patients during her initialization and finalization periods, (iii) the nurse must not serve more than a predefined number of patients at the same time. Once patient scheduling is complete for a campus, the algorithm calculates the overtime/under-utilization of all beds. If at least one of the beds experiences overtime, the heuristic returns one; otherwise, zero. Algorithm 2 and Figure EC.1 in the E-companion provide a pseudocode and an illustrative example for Heuristic 1, respectively.

Unlike Heuristic 1, Heuristic 2 does not necessarily schedule patients of one class in a row. For campus $c$ , this heuristic computes the proportion of bed time available for each class based on equation (12).

$$Q_k=\frac{(S_k^1+S_k^2+S_k^3)\times H_k}{\sum _{b\in {\mathcal{B}}_c}(U_c-L_b)\times G_{kb}}\mathrm{\forall }k\in \mathcal{K}$$

(12)

where $S_k^1$ , $S_k^2$ and $S_k^3$ refer to the lengths of the initialization, passivity and finalization periods for class $k$ , respectively, $H_k$ denotes the number of unscheduled patients of class $k$ , and $L_b$ refers to the last time slot that bed $b$ is occupied. Heuristic 2 identifies class $k$ with the maximum $Q_k$ (higher priority). It schedules one patient from this class to a feasible bed-nurse combination with the earliest service start time, computes $Q_k$ again and tries to schedule another patient. Heuristic 2 determines overtime/under-utilization of beds similar to Heuristic 1 when all patients are scheduled. Algorithm 3 and Figure EC.1 in the E-companion provide a pseudocode and an illustrative example for Heuristic 2, respectively. Based on our numerical results, the heuristics compute high-quality solutions, with an average gap of less than 4 $\mathrm{\%}$ in overtime value and around 99% accuracy in identifying feasible/infeasible solutions compared to a mathematical program (see the E-companion for more details).

Using heuristics and the NN introduces a potential source of approximation in our algorithm, which may misclassify feasible solutions as infeasible and vice versa. In other words, the algorithm may assume a solution is feasible (no overtime), while in fact it would require overtime. However, the effectiveness of our approach remains intact if the error at the design phase is small. It is worth noting that the primary purpose of our approach is to derive scheduling policies rather than directly implementable solutions. While approximation errors may impact the quality of these policies, smaller errors yield policies that are closer to the optimal. Once these approximate policies are developed, they can either be integrated into a mathematical model or used to create simple scheduling rules, ensuring that all overtime constraints are explicitly accounted for (to guarantee strict feasibility).

Network architecture and linearization: Let us assume a fully connected feed-forward NN with $nL+1$ layers, numbered from $0$ to $nL$ , where Layer $0$ corresponds to the input layer, Layers $1$ to $nL-1$ refer to hidden layers and layer $nL$ denotes the output layer. Each layer $l\in \mathcal{L}=\{0,1,2,\dots ,nL\}$ includes $n{N}_l$ neurons. Note that $n{N}_{nL}=1$ . If $g_l^n\in \mathbb{R}$ specifies the output of layer $l$ and neuron $n\in {\mathcal{N}}_l=\{1,2,\dots ,n{N}_L\}$ , $g_0^n$ is determined based on the value of neuron $n$ in the input layer (i.e., the number of allocated patients with class $n$ ). The outputs of subsequent layers ( $l\ge 1$ ) are computed using equation (13).

$$\begin{array}{r}{g}_l^{n^{\prime }}=\sigma (\sum _{n\in {\mathcal{N}}_{l-1}}{W}_{l-1}^{n{n}^{\prime }}\times g_{l-1}^n+{\mathrm{\Omega }}_l^{n^{\prime }})\mathrm{\forall }l\in \mathcal{L}/\{0\};n^{\prime }\in {\mathcal{N}}_l\end{array}$$

(13)

where $\sigma (\cdot )$ is an activation function, and $W_{l-1}^{n{n}^{\prime }}$ and ${\mathrm{\Omega }}_l^n$ are given matrices of weights and biases, respectively. Note that the weights and biases matrices could include negative entries. As in many applications (Fischetti and Jo, 2018), we could use the rectified linear unit (ReLU) activation function for the hidden layers in the NN, where $ReLU(g)=max\{0,g\}$ . However, this activation function would create several dead neurons in the NN. If the neurons are not activated initially, then, they are always in the off-state as zero gradients flow through them (Maas et al., 2013). Thus, we also consider a newer variant, the leaky ReLU (LReLU) activation function, which allows a small negative gradient flow through the network ( $LReLU(g)=max\{\alpha \times g,g\}$ where typically $0<\alpha \le 0.5$ ). Note that our choice of activation function is limited to piecewise linear activation functions, such as ReLU and LReLU, for the hidden layers because the implementation of other activation functions (e.g., tanh) into linear programing is not straightforward. ReLU and LReLU are nonlinear activation functions, due to which their functions must be linearized to incorporate the NN into a linear program. We adapt the approach proposed by Fischetti and Jo (2018). We redefine $g_l^n\ge 0$ , introduce two new auxiliary variables $q_l^n\ge 0$ and $u_l^n\in \{0,1\}$ , and use linear Constraint sets (14) to (19).

$$\begin{array}{rl}& \sum _{n^{\prime }\in {\mathcal{N}}_0}{W}_0^{n^{\prime }n}\times (g_0^{n^{\prime }}-q_0^{n^{\prime }})+{\mathrm{\Omega }}_1^n=g_1^n-q_1^n\mathrm{\forall }n\in {\mathcal{N}}_1\end{array}$$

(14)

$$\begin{array}{rl}& \sum _{n^{\prime }\in {\mathcal{N}}_{l-1}}{W}_{l-1}^{n^{\prime }n}\times (g_{l-1}^{n^{\prime }}-q_{l-1}^{n^{\prime }}\times \alpha )+{\mathrm{\Omega }}_l^n=g_l^n-q_l^n\\ & \mathrm{\forall }l\in \mathcal{L}/\{0,1\};n\in {\mathcal{N}}_l\end{array}$$

(15)

$$\begin{array}{rl}& \mathbb{1}\{g_l^n\le 0\}\text{ if }{u}_l^n=1\mathrm{\forall }l\in \mathcal{L};n\in {\mathcal{N}}_l\end{array}$$

(16)

$$\begin{array}{rl}& \mathbb{1}\{q_l^n\le 0\}\text{ if }{u}_l^n=0\mathrm{\forall }l\in \mathcal{L};n\in {\mathcal{N}}_l\end{array}$$

(17)

$$\begin{array}{rl}& g_l^n,q_l^n\ge 0\mathrm{\forall }l\in \mathcal{L};n\in {\mathcal{N}}_l\end{array}$$

(18)

$$\begin{array}{rl}& u_l^n\in \{0,1\}\mathrm{\forall }l\in \mathcal{L};n\in {\mathcal{N}}_l\end{array}$$

(19)

where $\mathbb{1}\{\cdot \}$ denotes the indicator function. In the output layer, we use the Sigmoid activation function for our classifier NN, where $S(g)=\frac{1}{1+e^{-g}}$ ( $S(g)\in (0,1),\mathrm{\forall }g\in (-\mathrm{\infty },+\mathrm{\infty })$ and $S^{-1}(g)=\ln \left(\frac{g}{1-g}\right)$ - the Logit function). If $S(g)\le 0.5$ , the solution is a feasible schedule with no overtime. Otherwise, the solution has overtime. To implement this activation function into a linear mathematical program, we do the following:

$$\begin{array}{rl}S(g_{nL}^1-q_{nL}^1)& \gtrless 0.5\Rightarrow S^{-1}\left[S(g_{nL}^1-q_{nL}^1)\right]\\ & \gtrless S^{-1}(0.5)\Rightarrow g_{nL}^1-q_{nL}^1\gtrless 0\end{array}$$

where $g_{nL}^1-q_{nL}^1$ is the value of the output layer before applying the activation function. Therefore, the solution is feasible if $g_{nL}^1-q_{nL}^1\le 0$ (or $F(\cdot )\le 0$ ). Figure EC.2 in the E-companion provides an illustrative example of how the NN classifies feasible and infeasible patient allocations.

3.2.6. Costs

Since costs of earlier stages have more importance for our problem, we use a discounted value function based on pre-decision state variables. As mentioned earlier, we aim to minimize a weighted sum of the wait time, deferral, ED use and transfer costs. The total immediate cost of being in state $\overrightarrow{s}$ and making decision $\overrightarrow{z}$ is calculated using equation (20), where ${\hat{\zeta }}_{koe}$ represents the probability that a patient of class $k$ , priority $o$ and $e$ days of wait time will use the ED.

$$\begin{array}{rl}cst(\overrightarrow{s},\overrightarrow{z})& =\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}\sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koietc}^1\times V_{ko(e+t)}^1\\ & +\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{\begin{array}{c}e\in {\mathcal{E}}_k\\ e<E_k\end{array}}(p_{koie}^w-\sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koietc}^1)\times V_{koe}^2\\ & +\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{\begin{array}{c}e\in {\mathcal{E}}_k\\ e<E_k\end{array}}([\sum _{i\in {\mathcal{I}}_k}{p}_{koie}^w-\sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koietc}^1]\times {\hat{\zeta }}_{koe})\times V_k^3\\ & +\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}(p_{koi{E}_k}^w-\sum _{\begin{array}{c}t\in \mathcal{T}\\ t\ge A{V}_{ko}\end{array}}\sum _{\begin{array}{c}c\in \mathcal{C}\\ \sum _{b\in {\mathcal{B}}_c}{G}_{kb}>0\end{array}}{z}_{koi{E}_{k}tc}^1)\times V_{koi}^4\end{array}$$

(20)

It should be noted that we compute the discounted wait time cost $V_{kot}^1$ similar to Patrick et al. (2008), which is decreasing in $o$ and is zero if a patient is scheduled within her wait time threshold. The immediate cost function does not capture the ED cost for scheduled patients given that the model cannot impact the use of the ED for scheduled patients (it captures the ED cost for patients on the waitlist). Thus, there is no need to include their ED cost in the immediate cost function. The optimality equation of the system at pre-decision state variable $\overrightarrow{s}$ , $v:\overrightarrow{s}\to {\mathbb{R}}_{\ge 0}$ , can be formulated as equation (21).

$$v(\overrightarrow{s})=\underset{\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}}{min}\{cst(\overrightarrow{s},\overrightarrow{z})+\lambda \times \sum _{{\overrightarrow{s}}^{\prime }\in S}pr({\overrightarrow{s}}^{\prime }|\overrightarrow{s},\overrightarrow{z})\times v({\overrightarrow{s}}^{\prime })\}\mathrm{\forall }\overrightarrow{s}\in S$$

(21)

where $0<\lambda \le 1$ is the discount factor.

4. Approximate Dynamic Programming

The state, decision and outcome spaces grow exponentially in MDP models as the problem size increases, which inevitably makes such models intractable. This section provides an approximate methodology to find computationally feasible solutions for large-sized instances of the proposed MDP model. An infinite-horizon ADP algorithm is proposed through which we determine appropriate coefficients for basis functions. Once the coefficients of the basis functions are computed (Phase 1), we can use them to determine patient scheduling (Phase 2). Note that we do not need to run the first phase each time we run the second phase. However, we can use the first phase and update the coefficients if the system configuration changes (e.g., demand mix variations).

4.1. Policy Determination

To implement our ADP algorithm, we reformulate the optimality equation as the following equivalent linear dynamic programing (LDP) model (Powell, 2007):

$$\begin{array}{rl}& \mathbf{\text{LDP:}}max\{\sum _{\overrightarrow{s}\in S}\eta (\overrightarrow{s})\times v(\overrightarrow{s})\}\end{array}$$

(22)

$$\begin{array}{rl}& \text{subjected to:}\\ & cst(\overrightarrow{s},\overrightarrow{z})+\lambda \times \sum _{{\overrightarrow{s}}^{\prime }\in S}pr({\overrightarrow{s}}^{\prime }|\overrightarrow{s},\overrightarrow{z})\times v({\overrightarrow{s}}^{\prime })\ge v(\overrightarrow{s})\mathrm{\forall }\overrightarrow{s}\in S;\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}\\ & v(\overrightarrow{s})\in {\mathbb{R}}_{\ge 0}\end{array}$$

where $\eta (\overrightarrow{s})$ is the strictly positive weight of state variable $\overrightarrow{s}$ in the objective function. We assume that $\sum _{\overrightarrow{s}\in S}\eta (\overrightarrow{s})=1$ . Model (22) has a variable for each state and a constraint for each feasible state-decision pair and thus remains intractable. To address this issue, we rely on an affine approximation of the value function $v(\overrightarrow{s})$ (or $v({\overrightarrow{p}}^w,{\overrightarrow{p}}^m)$ ) in the state variables as follows:

$$\begin{array}{rl}\tilde{v}({\overrightarrow{p}}^w,{\overrightarrow{p}}^m)& =V^0+\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}{p}_{koie}^w\times V_{koie}^w\\ & +\sum _{k\in \mathcal{K}}\sum _{t\in {\mathcal{T}}^{\prime }}\sum _{c\in \mathcal{C}}{p}_{ktc}^m\times V_{ktc}^m\mathrm{\forall }\overrightarrow{s}=({\overrightarrow{p}}^w,{\overrightarrow{p}}^m)\in S\end{array}$$

(23)

where $V^0\in \mathbb{R}$ , $V_{koie}^w\in {\mathbb{R}}_{\ge 0}$ and $V_{ktc}^m\in {\mathbb{R}}_{\ge 0}$ . Here, $V_{koie}^w$ can be interpreted as the marginal expected discounted cost of having an additional patient of class $k$ with priority $o$ , $i$ remaining appointments and $e$ days wait time on the waiting list, and $V_{ktc}^m$ refers to the marginal expected discounted cost of having an additional patient of class $k$ scheduled on day $t$ in campus $c$ . To determine these coefficients (known as the policy), we substitute equation (23) into Model (22) and formulate the linear approximate dynamic programing (LADP) model:

$$\begin{array}{rl}\mathbf{\text{LADP:}}& max\{V^0+\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}{\mathbb{E}}_{\eta }[p_{koie}^w]\times V_{koie}^w\\ & +\sum _{k\in \mathcal{K}}\sum _{t\in {\mathcal{T}}^{\prime }}\sum _{c\in \mathcal{C}}{\mathbb{E}}_{\eta }[p_{ktc}^m]\times V_{ktc}^m\}\end{array}$$

(24)

$$\begin{array}{rl}& \text{subjected to:}\\ & (1-\lambda )\times V^0+\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}{\pi }_{koie}(\overrightarrow{s},\overrightarrow{z})\times V_{koie}^w\\ & +\sum _{k\in \mathcal{K}}\sum _{t\in {\mathcal{T}}^{\prime }}\sum _{c\in \mathcal{C}}{\delta }_{ktc}(\overrightarrow{s},\overrightarrow{z})\times V_{ktc}^m\le cst(\overrightarrow{s},\overrightarrow{z})\\ & \mathrm{\forall }\overrightarrow{s}\in S;\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}\\ & V^0\in \mathbb{R};V_{koie}^w\in {\mathbb{R}}_{\ge 0};V_{ktc}^m\in {\mathbb{R}}_{\ge 0}\end{array}$$

where

$$\begin{array}{rl}& {\mathbb{E}}_{\eta }[p_{koie}^w]=\sum _{\overrightarrow{s}\in S}\eta (\overrightarrow{s})\times p_{koie}^w(\overrightarrow{s})\\ & \mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k\\ & {\mathbb{E}}_{\eta }[p_{ktc}^m]=\sum _{\overrightarrow{s}\in S}\eta (\overrightarrow{s})\times p_{ktc}^m(\overrightarrow{s})\mathrm{\forall }k\in \mathcal{K};t\in {\mathcal{T}}^{\prime };c\in \mathcal{C}\\ & {\pi }_{koi0}(\overrightarrow{s},\overrightarrow{z})=p_{koi0}^w(\overrightarrow{s})-\lambda \times {\delta }_{koi}\mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k\\ & {\pi }_{koie}(\overrightarrow{s},\overrightarrow{z})=p_{koie}^w(\overrightarrow{s})-\lambda \times \sum _{{\overrightarrow{s}}^{\prime }\in S}pr({\overrightarrow{s}}^{\prime }|\overrightarrow{s},\overrightarrow{z})\times p_{koie}^{{}^{\prime }w}(\overrightarrow{s},\overrightarrow{z})\\ & \mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k;e>0\\ & {\delta }_{ktc}(\overrightarrow{s},\overrightarrow{z})=p_{ktc}^m(\overrightarrow{s})-\lambda \times \sum _{{\overrightarrow{s}}^{\prime }\in S}pr({\overrightarrow{s}}^{\prime }|\overrightarrow{s},\overrightarrow{z})\times p_{ktc}^{{}^{\prime }m}(\overrightarrow{s},\overrightarrow{z})\\ & \mathrm{\forall }k\in \mathcal{K};t\in {\mathcal{T}}^{\prime }/T^{\prime };c\in \mathcal{C}\\ & {\delta }_{k{T}^{\prime }c}(\overrightarrow{s},\overrightarrow{z})=p_{k{T}^{\prime }c}^m(\overrightarrow{s})\mathrm{\forall }k\in \mathcal{K};c\in \mathcal{C}\end{array}$$

Model (24) has a tractable number of variables, but is still intractable due to the number of constraints (one for each state-decision pair). We therefore formulate the dual of Model (24), as the restricted master problem (RMP), and solve it using CG.

$$\begin{array}{rl}\mathbf{\text{RMP:}}& min\sum _{\overrightarrow{s}\in S}\sum _{\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}}cst(\overrightarrow{s},\overrightarrow{z})\times \chi (\overrightarrow{s},\overrightarrow{z})\end{array}$$

(25)

$$\begin{array}{rl}& \text{subjected to:}\\ & (1-\lambda )\times \sum _{\overrightarrow{s}\in S}\sum _{\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}}\chi (\overrightarrow{s},\overrightarrow{z})=1\\ & \sum _{\overrightarrow{s}\in S}\sum _{\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}}{\pi }_{koie}(\overrightarrow{s},\overrightarrow{z})\times \chi (\overrightarrow{s},\overrightarrow{z})\ge {\mathbb{E}}_{\eta }[p_{koie}^w]\\ & \mathrm{\forall }k\in \mathcal{K};o\in \mathcal{O};i\in {\mathcal{I}}_k;e\in {\mathcal{E}}_k\\ & \sum _{\overrightarrow{s}\in S}\sum _{\overrightarrow{z}\in {\mathrm{\Xi }}_{\overrightarrow{s}}}{\delta }_{ktc}(\overrightarrow{s},\overrightarrow{z})\times \chi (\overrightarrow{s},\overrightarrow{z})\ge {\mathbb{E}}_{\eta }[p_{ktc}^m]\\ & \mathrm{\forall }k\in \mathcal{K};t\in {\mathcal{T}}^{\prime };c\in \mathcal{C}\\ & \chi (\overrightarrow{s},\overrightarrow{z})\in {\mathbb{R}}_{\ge 0}\end{array}$$

where $\chi (\overrightarrow{s},\overrightarrow{z})$ represents the dual variable related to the state-decision pair $(\overrightarrow{s},\overrightarrow{z})$ . We use the dual values associated with constraints in Model (25) to determine $V^0$ , $V_{koie}^w$ and $V_{ktc}^m$ . In order to determine the next state-decision pair to include in Model (25) (RMP), we formulate the sub-problem (SP) (or the pricing problem) as follows:

$$\begin{array}{rl}\mathbf{\text{SP:}}& max\{(1-\lambda )\times V^0+\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}{\pi }_{koie}(\overrightarrow{s},\overrightarrow{z})\times V_{koie}^w\\ & +\sum _{k\in \mathcal{K}}\sum _{t\in {\mathcal{T}}^{\prime }}\sum _{c\in \mathcal{C}}{\delta }_{ktc}(\overrightarrow{s},\overrightarrow{z})\times V_{ktc}^m-cst(\overrightarrow{s},\overrightarrow{z})\}\end{array}$$

(26)

$$\begin{array}{rl}& \text{subjected to:}\\ & \text{Constraint sets }(7)-(11)\end{array}$$

Model (26) looks for a state-decision pair with the most violated primal constraint. We add this state-decision pair to Model (25) in order to re-optimize the master problem and update the approximation coefficients. The CG procedure terminates when the algorithm is sufficiently close to the non-violation of primal constraints (Model (26) with objective function value of ${10}^{-5}$ ). In order to start the CG procedure, we apply the Phase 1 method described by Amor et al. (2009) to find an initial feasible solution.

Once the approximate optimal policy (AOP) is determined, we use it to compute the cost of scheduling a patient with class $k$ , priority $o$ , $i$ remaining appointments and $e$ days wait time as follows:

$$V_{koietc}=\{\begin{array}{ll}{V}_{ko(e+t)}^1+\lambda \times {\displaystyle \sum _{\begin{array}{c}{i}^{\prime }\in \{1,2,\dots ,i\}\\ t^{\prime }=(i-1)\times F{R}_k+t-1\end{array}}{V}_{k{t}^{\prime }c}^m-(V_{koe}^2+\mathbb{P}({\zeta }_{koie})\times V_k^3+\lambda \times V_{koi(e+1)}^w)}& \mathrm{\forall }e<E_k;t\in \mathcal{T};c\in \mathcal{C}\\ V_{ko(e+t)}^1+\lambda \times {\displaystyle \sum _{\begin{array}{c}{i}^{\prime }\in \{1,2,\dots ,i\}\\ t^{\prime }=(i-1)\times F{R}_k+t-1\end{array}}{V}_{k{t}^{\prime }c}^m-V_{koi}^4}& \mathrm{\forall }e=E_k;t\in \mathcal{T};c\in \mathcal{C}\end{array}$$

Then, we can use the action coefficient $V_{koietc}$ in model (7) to (11) to conduct advance scheduling. This coefficient balances (i) the cost of the wait time and loss of available treatment capacity in the future against the cost of the deferral and having an extra patient on the waiting list, $\mathrm{\forall }e<E_k;t\in \mathcal{T};c\in \mathcal{C}$ , and (ii) the cost of the wait time and loss of available treatment capacity in the future against the transfer cost, $\mathrm{\forall }e=E_k;t\in \mathcal{T};c\in \mathcal{C}$ . For a patient on the waiting list with class $k$ , priority $o$ , $i$ remaining appointments, and $e$ days wait time, we transfer her if $V_{koi{e}^{\prime }tc}>0$ , $\mathrm{\forall }{e}^{\prime }\in \{e,e+1,\dots ,E_k\};t\in \mathcal{T};c\in \mathcal{C}$ . Otherwise, she is scheduled only when $V_{koietc}\le 0$ or deferred.

4.2. Acceleration Techniques

Despite the use of the NN, our original ADP algorithm is unable to derive policies for medium- and large-sized instances ( $K\ge 3$ , $O\ge 3$ , $I\ge 3$ ) in a reasonable CPU time. Thus, we adapt two acceleration techniques to alleviate the computational complexity. First, we implement the dual stabilization method offered by Addis et al. (2012); Wang et al. (2017), wherein the dual multipliers are stabilized by combining them with a centroid dual vector. Second, inspired by Adelman and Mersereau (2008), we propose a heuristic logic for the CG algorithm by relaxing the original SP via decomposition. We randomly sort campuses at each iteration and solve the SP for each one separately. The acceleration techniques are discussed in depth in the E-companion.

5. Computational Results

This section presents comprehensive numerical analyses to evaluate the proposed solution approach and derive managerial implications. We explain the settings in Section 5.1, introduce new scheduling rules in Section 5.2, and assess the AOP in Section 5.3. Finally, Section 5.4 investigates a case study, and Section 5.5 provides managerial insights. Both master- and SPs of the CG in the ADP algorithm are solved via IBM Ilog CPLEX. The ADP algorithm and all other programs are coded in PyCharm 2020.3.5, in which the interpreter is Python 3.7, and TensorFlow 2.8.0 is used to train the NN. All experiments are run on an Intel(R) Core(TM) i9-10900X CPU @ 3.70GHz with 64 GB of DDR4 RAM.

5.1. Analysis Settings

This sub-section describes instance generation, hyperparameter tuning, and benchmark scheduling rules.

5.1.1. Instance Configurations

We assume that a day includes a six-/eight-hour shift (no overtime available), and time is discretized into time slots of 15, 30, or 60 minutes. The decision horizon includes 7 or 14 days, and the planning horizon includes 7 to 70 days. The service durations vary between three and eight time slots. Treatments are associated with 1 to 12 classes and one to four priority levels, and may require one to nine treatment sessions. The treatment frequency for all classes is one week, two weeks, or a month. Similar to patient safety guidelines in our partner hospital, we do not allow a nurse to handle more than four patients at the same time. The wait time and deferral thresholds may be 2 to 12 days and 1 to 24 days, respectively.

We let both the transfer and ED utilization costs be greater than the maximum possible sum of wait time and deferral costs (i.e., for both $V_{koi}^4$ and $V_k^3$ , $V>V_{ko(E_k+T+1)}^1+\sum _{e\in {\mathcal{E}}_k}{V}_{koe}^2$ ), with the former being the most expensive (to discourage patient transfers). Through sensitivity analyses, we also investigate the scenario where the transfer cost is significantly greater than the other costs (i.e., $V_{koi}^4\gg V_{ko(E_k+T+1)}^1+\sum _{e\in {\mathcal{E}}_k}{V}_{koe}^2$ ). We assume that a transfer is costlier for patients with a larger number of appointments ( $V_{ko(i+1)}^4>V_{koi}^4$ ) while it is relatively more expensive for patients with fewer number of appointments $(V_{ko(i+1)}^4<V_{koi}^4\times \left[\frac{i+1}{i}\right])$ . Thus, $V_{koi}^4<V_{ko(i+1)}^4<V_{koi}^4\times \left[\frac{i+1}{i}\right]$ . However, we analyze an alternative scenario where the transfer cost is relatively equal for multi-appointment treatments $(\text{i.e., }{V}_{koi}^4\ge V_{ko(i+1)}^4\times \left[\frac{i}{i+1}\right])$ . We consider two demand settings: medium and high. For the medium demand setting, we adjust arrival rates so that the transfer rate for all class-priorities is 10% or less, based on simulations using a Myopic policy presented in Section 5.1.3. For the high demand setting, we increase the arrival rates by 20%. Ten test instances are generated according to the above-mentioned settings, and their configurations are summarized in Table 1. We have provided all the information on test instances in a data set repository (https://doi.org/10.17632/bnvf9jdmyb.2).

Table 1.

Configurations of test instances.

Inst.	$C$	$N$	$B$	$U$	$T$	$K$	$O$	$\underset{k\in \mathcal{K}}{max}{I}_k$	$\underset{k\in \mathcal{K}}{max}{E}_k$
1	1	2	4	8	7	1	1	1	1
2	1	2	4	8	7	1	1	3	1
3	1	2	4	8	7	1	1	1	3
4	1	4	8	12	7	3	2	3	6
5	1	4	8	12	7	4	2	4	9
6	2	4	12	16	14	5	3	5	12
7	2	6	12	16	14	6	3	6	15
8	2	6	16	20	14	8	4	7	18
9	2	8	16	20	14	10	4	8	21
10	2	8	20	24	14	12	4	9	24

5.1.2. Hyperparameter Tuning and Predictive Model Performance

The proposed ADP algorithm relies on a NN to classify patient allocations based on overtime. To train the NN, our loss function penalizes errors equally for both “predicting feasibility when infeasible” and “predicting infeasibility when feasible.” While this balanced approach may not suit all applications, our method is flexible and can be easily adapted by modifying the loss function to better align with specific preferences. Furthermore, the NN contains hyperparameters that must be tuned before use, including the numbers of hidden layers and their neurons, number of epochs, activation function, optimizer and $\alpha$ for the LReLU activation function (Elakkiya and Selvakumar, 2022) (their potential values are presented in Table 2). Tuning the aforementioned hyperparameters is important because the performance of our NN might depend on their values. It is typical to evaluate the performance of NNs by a $K$ -fold cross-validation technique (Wong and Yeh, 2019). We apply this evaluation method for hyperparameter tuning, where $K=10$ , and use a weighted evaluation metric that includes (i) accuracy of the test fold $k$ ( $ac{c}_k$ ), (ii) number of activated hidden layers ( $nL$ ) and their neurons ( $n{N}_l$ ), and (iii) number of epochs ( $nE$ )—length of training. This function is computed based on equation (27). Note that we apply the commonly used Binary Cross-Entropy as the loss function for our classifier NN.

$$\begin{array}{rl}FF& =w{g}_1\times \frac{\sum _{k\in \{1,2,\dots ,K\}}ac{c}_k}{K}-w{g}_2\times \frac{nL}{10}-w{g}_3\\ & \times \frac{\sum _{l\in \mathcal{L}}n{N}_l}{1000}-w{g}_4\times \frac{nE}{5000}\end{array}$$

(27)

where $w{g}_{\kappa }>0$ are constants representing the weights of different components in the fitness function (we set $w{g}_1=0.7$ and $w{g}_2=w{g}_3=w{g}_4=0.1$ ). Such a multi-component fitness function is beneficial for feed-forward NNs to achieve generalization (Ojha et al., 2017). The iterated racing package proposed by López-Ibánez et al. (2016) for automatic algorithm configuration (available on R as the iRace package) is used for hyperparameter tuning, where the termination criterion is set equal to 10,000 fitness function evaluations. After running this package, the elite configurations are found and reported in Table 2.

Table 2.

Parameter tuning using the iRace package.

Parameters	Penitential values	Elite configuration
# hidden layers	{0, 1, …, 10}	1
# neurons	{5, 10, …, 100}	20
# epochs	{50, 100, …, 5000}	600
Activation function	ReLU, or LReLU	LReLU
Optimizer	rmsprop, adam, adagrad, or adadelta	adam
$\alpha$	{0.01, 0.05, 0.1, 0.2, 0.3}	0.01

ReLU: rectified linear unit; LReLU: leaky rectified linear unit.

To evaluate the NN using the elite configuration, we employed the 10-fold cross-validation method ten times across Instances 6–10. Instead of using the entire validation set, we focus on patient allocations for a campus-day that utilizes between 75% and 125% of the total available bed capacity $(\text{i.e., }75\mathrm{\%}\le \left(\frac{\sum _{k\in \mathcal{K}}\sum _{o\in \mathcal{O}}\sum _{i\in {\mathcal{I}}_k}\sum _{e\in {\mathcal{E}}_k}(z_{koietc}^1+z_{ktc}^2)\times (S_k^1+S_k^2+S_k^3)}{|{\mathcal{B}}_c|\times U_c}\right)\times 100\le 125\mathrm{\%})$ . The NN achieved an average validation accuracy of approximately 95.40%. In comparison, logistic regression achieved an average accuracy of around 84.20%, significantly lower than the NN with a confidence level of 99% (see Figure 2(a)). Moreover, the NN achieves an area under the receiver operating characteristic curve of around 0.94, surpassing the 0.85 achieved by logistic regression (see Figure 2(b)). Thus, the NN consistently outperforms logistic regression across all instances with higher true positive rates and lower false positive rates. The NN’s performance is more stable and reliable, indicating its superior ability to accurately classify feasible/infeasible solutions. Further analysis in the E-companion demonstrates the NN’s solid performance in the downstream optimization problem.

Figure 2.

Overall performance of the neural network (NN) against logistic regression. (a) Accuracy; (b) area under the receiver operating characteristic curve.

5.1.3. Benchmark Scheduling Rules

We compare the performance of the AOP with three well-known heuristic scheduling rules.

Myopic policy: This policy ignores the impact of today’s decision on the future and books patients as soon as possible based on the immediate cost function (20). The Myopic policy defers a patient only when no capacity is available within the decision horizon and transfers a patient if the maximum deferral threshold is reached.

Adapted Patrick, Puterman, and Queyranne (APPQ) policy: Patrick et al. (2008) suggests booking patients in priority order within their wait time thresholds. We adapt this policy to our problem as follows. For each priority level, we sort classes based on the ascending order of their wait time thresholds. Then, for each class-priority, we prioritize scheduling patients with a greater wait time and fewer treatment sessions required. Then, we book as many patients as possible into days $W{T}^{\prime }=\{A{V}_{ko},A{V}_{ko}+1,\dots ,W{T}_{ko}-e\}$ for the highest priority level. For lower priority levels, book as many patients as possible into days $W{T}^{\prime }=\{A{V}_{ko},W{T}_{ko}-e,W{T}_{ko}-e-1,\dots ,A{V}_{ko}+1\}$ or $W{T}^{\prime }=\{W{T}_{ko}-e,W{T}_{ko}-e-1,\dots ,A{V}_{ko}\}$ when $A{V}_{ko}=1$ or $A{V}_{ko}>1$ , respectively. A deferral occurs when there is no available capacity within days $W{T}^{\prime }$ and $|W{T}^{\prime }|>0$ , and a transfer occurs once $|W{T}^{\prime }|=0$ .

Adapted days with minimum number of booking (ADMB) policy: In the order of priority, this policy books as many patients as possible on the first day and then the day with the minimum number of bookings within the wait time threshold (Sauré et al., 2015). We weigh treatments by their service duration to identify the day with the minimum number of bookings and adapt this policy to our problem similar to the APPQ policy.

5.2. New Scheduling Rules

One approach for creating schedules is to solve a mathematical model incorporating AOP in its objective function. However, implementing this approach might not be straightforward. Observing the AOP for various instances and based on insights derived in Section 5.5, we found two generalizable mechanisms that could be deployed in easy-to-implement scheduling rules to improve their performance for our problem. Unlike the APPQ and ADMB policies, which book patients strictly within the wait time threshold, we first extend the booking horizon for the highest-priority patients (e.g., by one day in our instances). Second, we limit capacity for lower-priority patients when the transfer rate for higher-priority patients exceeds a specific threshold (e.g., 5% in our instances). We set maximum allocation limits for each class-priority in our numerical analyses to implement the second mechanism. If the average rejection rate for higher-priority patients exceeds this threshold, we reduce the maximum allocation for lower-priority patients by a specified percentage (e.g., 10% reduction) until either their maximum capacity reaches zero or the rejection rate for higher-priority patients falls below the threshold.

We deploy these two mechanisms in the APPQ and ADMB policies in order to generate two new scheduling policies, called Enhanced APPQ (EPPQ) and ADMB (EDMB). These policies are also evaluated in the following sub-sections.

5.3. Theoretical Instances

In this sub-section, we interpret the AOP and compare its performance with established scheduling rules in the literature.

5.3.1. Comparison of Scheduling Rules

We use a simulation model to evaluate scheduling rules over Instances 1–10 with two demand settings. The simulation model consists of 200 warm-up iterations and 1000 main iterations. We run the simulation model 100 times for each instance and present the results in Table 3. Each cell under “Heuristic scheduling rules” includes a 99% confidence interval for the average cost. Highlighted confidence intervals are significantly worse than the ones related to the AOP. Each cell under “Gap (%)” demonstrates the minimum and maximum gaps of the scheduling policies in comparison to the AOP, $([F_{policy}-F_{AOP}]/F_{AOP})\times 100$ , respectively. The Myopic policy consistently provides the poorest performance, with a 169.19% gap compared to AOP. In Instances 1-3, all other policies perform similarly because these instances include only one class-priority (i.e., no flexibility to manage excessive demand).

Table 3.

Comparison of performance metrics for heuristic scheduling rules.

		Heuristic scheduling rules
Inst.	Demand	Myopic	APPQ	ADMB	EPPQ	EDMB	AOP	Gap (%)
1	Normal	{\columncolor{tabgray}}l}{\phantom{00,}?>17.97 $\pm$ 0.76	8.47 $\pm$ 0.30	8.48 $\pm$ 0.31	10.07 $\pm$ 0.38	10.06 $\pm$ 0.41	8.56 $\pm$ 0.31	[ $-$ 1.11, 109.88]
	High	{\columncolor{tabgray}}l}{\phantom{00,}?>28.38 $\pm$ 0.43	17.73 $\pm$ 0.27	17.72 $\pm$ 0.26	20.03 $\pm$ 0.32	20.04 $\pm$ 0.27	18 $\pm$ 0.26	[ $-$ 1.56, 57.72]
2	Normal	{\columncolor{tabgray}}l}{\phantom{00,}?>14.23 $\pm$ 0.68	9.19 $\pm$ 0.39	9.21 $\pm$ 0.37	9.46 $\pm$ 0.40	9.50 $\pm$ 0.40	9.24 $\pm$ 0.38	[ $-$ 0.50, 54.02]
	High	{\columncolor{tabgray}}l}{\phantom{00,}?>22.47 $\pm$ 0.51	14.41 $\pm$ 0.32	14.41 $\pm$ 0.32	15.47 $\pm$ 0.39	15.47 $\pm$ 0.37	14.54 $\pm$ 0.32	[ $-$ 0.94, 54.41]
3	Normal	{\columncolor{tabgray}}l}{\phantom{00,}?>34.67 $\pm$ 0.61	25.56 $\pm$ 0.43	25.52 $\pm$ 0.42	27.45 $\pm$ 0.47	26.24 $\pm$ 0.46	25.71 $\pm$ 0.42	[ $-$ 0.71, 34.87]
	High	{\columncolor{tabgray}}l}{\phantom{00,}?>45.90 $\pm$ 0.81	35.66 $\pm$ 0.54	35.58 $\pm$ 0.56	37.48 $\pm$ 0.56	36.37 $\pm$ 0.53	35.80 $\pm$ 0.56	[ $-$ 0.63, 28.19]
4	Normal	{\columncolor{tabgray}}l}{\phantom{0,}?>198.81 $\pm$ 3.34	{\columncolor{tabgray}}l}{\phantom{0,}?>114.17 $\pm$ 2.47	{\columncolor{tabgray}}l}{\phantom{0,}?>120.28 $\pm$ 2.74	{\columncolor{tabgray}}l}{\phantom{0,}?>100.75 $\pm$ 6.53	{\columncolor{tabgray}}l}{\phantom{00,}?>84.38 $\pm$ 0.82	80.29 $\pm$ 0.78	[5.10, 147.62]
	High	{\columncolor{tabgray}}l}{\phantom{0,}?>361.84 $\pm$ 4.37	{\columncolor{tabgray}}l}{\phantom{0,}?>228.01 $\pm$ 3.09	{\columncolor{tabgray}}l}{\phantom{0,}?>243.50 $\pm$ 3.52	{\columncolor{tabgray}}l}{\phantom{0,}?>181.35 $\pm$ 1.05	177.13 $\pm$ 6.29	173.04 $\pm$ 4.78	[2.36, 109.11]
5	Normal	{\columncolor{tabgray}}l}{\phantom{0,}?>250.09 $\pm$ 5.71	{\columncolor{tabgray}}l}{\phantom{0,}?>123.45 $\pm$ 3.33	{\columncolor{tabgray}}l}{\phantom{0,}?>135.43 $\pm$ 4.03	{\columncolor{tabgray}}l}{\phantom{00,}?>78.03 $\pm$ 3.49	{\columncolor{tabgray}}l}{\phantom{00,}?>71.40 $\pm$ 1.05	67.79 $\pm$ 1	[5.33, 271.67]
	High	{\columncolor{tabgray}}l}{\phantom{0,}?>360.68 $\pm$ 5.01	{\columncolor{tabgray}}l}{\phantom{0,}?>182.61 $\pm$ 3.82	{\columncolor{tabgray}}l}{\phantom{0,}?>201.30 $\pm$ 3.21	{\columncolor{tabgray}}l}{\phantom{0,}?>150.97 $\pm$ 12.96	116.50 $\pm$ 2.32	113.16 $\pm$ 12.25	[2.96, 218.75]
6	Normal	{\columncolor{tabgray}}l}{\phantom{0,}?>974.13 $\pm$ 23.69	{\columncolor{tabgray}}l}{\phantom{0,}?>551.99 $\pm$ 19.17	{\columncolor{tabgray}}l}{\phantom{0,}?>579.93 $\pm$ 19.97	{\columncolor{tabgray}}l}{\phantom{0,}?>527.37 $\pm$ 67.97	{\columncolor{tabgray}}l}{\phantom{0,}?>317.35 $\pm$ 8.63	304.30 $\pm$ 8.28	[4.29, 220.12]
	High	{\columncolor{tabgray}}l}{?>1,367.08 $\pm$ 18.03	{\columncolor{tabgray}}l}{\phantom{0,}?>866.96 $\pm$ 13.93	{\columncolor{tabgray}}l}{\phantom{0,}?>928.27 $\pm$ 14.05	{\columncolor{tabgray}}l}{\phantom{0,}?>631.30 $\pm$ 3.33	{\columncolor{tabgray}}l}{\phantom{0,}?>629.67 $\pm$ 3.60	617.56 $\pm$ 11.27	[1.96, 121.37]
7	Normal	{\columncolor{tabgray}}l}{?>1,087.97 $\pm$ 27.40	{\columncolor{tabgray}}l}{\phantom{0,}?>660.25 $\pm$ 20.13	{\columncolor{tabgray}}l}{\phantom{0,}?>725.89 $\pm$ 24.68	{\columncolor{tabgray}}l}{\phantom{0,}?>476.71 $\pm$ 4.72	{\columncolor{tabgray}}l}{\phantom{0,}?>475.09 $\pm$ 4.65	455.93 $\pm$ 2.30	[4.20, 138.63]
	High	{\columncolor{tabgray}}l}{?>2,011.16 $\pm$ 31.81	{\columncolor{tabgray}}l}{?>1,344.52 $\pm$ 28.35	{\columncolor{tabgray}}l}{?>1,506.14 $\pm$ 29.89	{\columncolor{tabgray}}l}{?>1,372.18 $\pm$ 177.39	1,106.32 $\pm$ 63.26	1,088.29 $\pm$ 61.53	[1.66, 84.80]
8	Normal	{\columncolor{tabgray}}l}{?>1,847.70 $\pm$ 80.82	{\columncolor{tabgray}}l}{?>1,063.65 $\pm$ 59.82	{\columncolor{tabgray}}l}{?>1,134.04 $\pm$ 71.39	{\columncolor{tabgray}}l}{\phantom{0,}?>790.89 $\pm$ 154.63	516.88 $\pm$ 81.94	492.15 $\pm$ 78.02	[5.02, 275.44]
	High	{\columncolor{tabgray}}l}{?>3,913.71 $\pm$ 58.05	{\columncolor{tabgray}}l}{?>2,754.82 $\pm$ 52.83	{\columncolor{tabgray}}l}{?>2,967.77 $\pm$ 57.71	{\columncolor{tabgray}}l}{?>1,434.88 $\pm$ 13.61	1,382.96 $\pm$ 69.70	1,365.76 $\pm$ 46.93	[1.26, 186.56]
9	Normal	{\columncolor{tabgray}}l}{?>1,276.45 $\pm$ 39.75	{\columncolor{tabgray}}l}{\phantom{0,}?>625.81 $\pm$ 25.88	{\columncolor{tabgray}}l}{\phantom{0,}?>804.53 $\pm$ 33.18	{\columncolor{tabgray}}l}{\phantom{0,}?>786.41 $\pm$ 20.52	{\columncolor{tabgray}}l}{\phantom{0,}?>770.28 $\pm$ 8.12	591.99 $\pm$ 24.48	[5.71, 115.62]
	High	{\columncolor{tabgray}}l}{?>2,563.95 $\pm$ 42.32	{\columncolor{tabgray}}l}{?>1,584.51 $\pm$ 34.34	{\columncolor{tabgray}}l}{?>1,926.06 $\pm$ 37.41	{\columncolor{tabgray}}l}{?>1,152.92 $\pm$ 16.23	{\columncolor{tabgray}}l}{?>1,105.96 $\pm$ 10.62	1,087.61 $\pm$ 10.45	[1.69, 135.74]
10	Normal	{\columncolor{tabgray}}l}{?>1,685.19 $\pm$ 39.93	{\columncolor{tabgray}}l}{\phantom{0,}?>973.84 $\pm$ 36.02	{\columncolor{tabgray}}l}{?>1,029.96 $\pm$ 38.55	{\columncolor{tabgray}}l}{\phantom{0,}?>937.92 $\pm$ 65.20	919.93 $\pm$ 80.58	865.24 $\pm$ 19.38	[6.32, 92.17]
	High	{\columncolor{tabgray}}l}{?>3,878.89 $\pm$ 56.83	{\columncolor{tabgray}}l}{?>2,626.46 $\pm$ 53.59	{\columncolor{tabgray}}l}{?>2,998.90 $\pm$ 53.38	{\columncolor{tabgray}}l}{?>1,512.49 $\pm$ 42.21	1,469.27 $\pm$ 37.75	1,440.94 $\pm$ 35.09	[1.97, 169.19]

Each cells under “Heuristic scheduling rules” includes a 99% confidence interval for cost. Highlighted intervals are significantly worse than the one related to the AOP. Each cell under “Gap (%)” demonstrates the smallest and largest gaps of the scheduling policies in comparison to the AOP, respectively. AOP: approximate optimal policy; APPQ: adapted Patrick, Puterman, and Queyranne; ADMB: adapted days with minimum number of booking; EPPQ: enhanced APPQ; EDMB: enhanced ADMB; .

The APPQ policy usually outperforms the ADMB policy. Unlike APPQ, which schedules patients from the end of the wait time threshold backward (except for the highest priority levels), ADMB books patients on the day with the fewest bookings within their wait time threshold. This may result in booking lower-priority patients earlier (even when later days have free capacity), using space needed for higher-priority patients and potentially leading to their transfer in the next epochs. The EPPQ and EDMB policies outperform their APPQ and ADMB counterparts (over Instances 4–10) with performance gaps of 30.13% and 63.59%, respectively, as they proactively alleviate excessive demand pressure. Notably, EDMB slightly outperforms EPPQ when extra demand is intelligently handled (as explained in Section 5.2), similar to the observations for the original version of these policies (Sauré et al., 2012).

According to Table 3, the AOP consistently outperforms other policies across Instances 4–10, with the EDMB policy providing the closest performance (a gap of 5.30% and no significant difference in 7 out of 14 instances). In particular, EDMB maintains its performance in larger instances, with a maximum gap of 6.32% compared to AOP, and a gap of only 1.97% in Instance 10, a high-demand setting. Comparing normal and high-demand settings, we observe that performance gaps decrease under high demand (e.g., for EDMB, 8.63% gap in normal demand vs. 1.98% in high demand). This is because higher demand pressures policies to handle transfers more effectively, reducing AOP’s flexibility and bringing the more conservative EPPQ and EDMB policies closer to the AOP’s performance.

5.4. Case Study

Our partner hospital offers treatments on three campuses: Campus 1 with four nurses and eight beds, Campus 2 with three nurses and nine beds, and Campus 3 with three nurses and 13 beds. These campuses serve a total of 15 classes, each with five priority levels. Priority level 1 has a two-day wait time threshold and a five-day deferral limit while priority level 5 has a 28-day wait time threshold and a 45-day deferral limit. All campuses are open for eight hours or 32 time slots (time is discretized into 15 minutes time slots) with no overtime. Service duration varies between four and 28 time slots. For all classes, the initialization/finalization period is one time slot, and the passivity period ranges from two to 26 time slots. The booking clerks of our collaborator hospital suggested assuming a 14-day decision horizon and a 56-day planning horizon (no patient should be scheduled later than two months). Treatments can be provided on a weekly, bi-weekly, or monthly basis (varies based on class), therefore, they may need one to six appointments. The rest of the parameters are defined in our data set repository.

It should be noted that for some treatment classes, the length of treatment may exceed the 56-day planning horizon. In such cases, if the total treatment length is known in advance for a patient, the booking clerk can readily allocate her as soon as the rolling planning horizon includes her future appointment dates. On the other hand, if the total treatment length is unknown, and the patient needs to consult with their physician to assess the necessity of additional appointments, the patient can be added to the list of scheduled patients based on their previous appointments. If this approach is infeasible (e.g., due to lack of capacity), the booking clerk has the flexibility to manually schedule her or add her to the waitlist of high-priority arrivals.

We construct an instance for our case study and compute its AOP. Then, we run a simulation model to compare it against the Myopic, APPQ, ADMB, EPPQ and EDMB policies. We also compare all policies against a template-based scheduling (TBS) approach (i.e., fixed allocation of patient classes to beds over the planning horizon) currently used in our partner hospital. A mathematical program is used to build the scheduling templates (see the E-companion for further details). We utilize the scheduling template and perform advance scheduling similar to the APPQ policy, but with two differences: (i) the capacity for each class is fixed on each day-campus, (ii) new arrivals are scheduled within the decision horizon regardless of their wait time thresholds. Those who cannot be scheduled within the deferral threshold are transferred. Summary results for our case study are presented in Table 4.

Table 4.

Comparison of the scheduling policies for the case study.

	Cost
Policies	Total	Deferral (avg.)	Wait time (avg.)	ED use (avg.)	Transfer (avg.)	On-time treatments (%)
Myopic	{\columncolor{tabgray}}l}{?>2,241.23 $\pm$ 92.67	439.57	514.28	280.07	1,007.28	66.68
TBS	{\columncolor{tabgray}}l}{?>2,041.93 $\pm$ 56.57	319.06	273.20	247.18	1,202.49	80.96
APPQ	{\columncolor{tabgray}}l}{?>1,166.61 $\pm$ 57.36	8.53	0	165.93	992.16	94.20
ADMB	{\columncolor{tabgray}}l}{?>1,368.22 $\pm$ 74.97	6.02	0	177.04	1,185.17	89.56
EPPQ	{\columncolor{tabgray}}l}{\phantom{0,}?>881.32 $\pm$ 69.63	2.02	0.85	188.74	689.71	95.13
EDMB	708.26 $\pm$ 57.16	3.01	0.57	211.24	493.44	94
AOP	679.71 $\pm$ 35.39	8.12	1.40	217.08	453.11	94.76

Each cells under “Total cost” includes a 99% confidence interval for average cost. Highlighted confidence intervals are significantly worse than the one related to the AOP.AOP: approximate optimal policy; ED: emergency department; APPQ: adapted Patrick, Puterman, and Queyranne; ADMB: adapted days with minimum number of booking; EPPQ: enhanced APPQ; EDMB: enhanced ADMB; TBS: template-based scheduling.

In line with our earlier findings, Table 4 shows that the Myopic policy performs the worst by a large margin, followed by the TBS policy. These policies incur the highest deferral and wait time costs, leading to increased ED use. Additionally, they incur the highest transfer costs due to their failure to account for the dynamic nature of new arrivals, greedily utilizing capacity that could have been reserved for higher-priority patients. The EPPQ and EDMB policies outperform their counterparts and perform significantly better than the Myopic and TBS policies. While the EPPQ and EDMB handle patient transfers better, they incur higher ED costs compared to APPQ and ADMB policies as they more frequently schedule higher-priority patients who utilize the ED more often. The AOP policy performs significantly better than other policies, except for EDMB. It rarely defers patients or books them beyond their wait time thresholds. Under this policy, 94.76% of new arrivals are served on-time, with the EPPQ and Myopic policies offering the best and worst rate at 95.13% and 66.68%, respectively.

5.5. Managerial Implications

This sub-section presents managerial insights on (i) the higher performance of the AOP, (ii) the reasons for the TBS policy’s inferiority, (iii) the impact of considering patients’ campus preferences, and (iv) the effect of objective function coefficients on the performance metrics. Further insights from the AOP are provided in the E-companion.

5.5.1. Superiority of the AOP

We analyze performance metrics for all scheduling policies in Instance 5 with normal demand (where the AOP exhibits the largest gap compared to EDMB policy for instances with one campus—see Table 3). This analysis is depicted in Figure 3. Here, wait time refers to the number of days patients waited to receive their first treatment (excluding transferred patients) and deferral denotes the number of times patients are deferred until they are either scheduled or transferred. We calculate the transfer rate as patient transfers divided by all arrivals. For each priority, the metrics represent average values over all treatment classes. Note that the bed utilization rate is limited to $\approx$ 93% for all scheduling policies given the limited availability of nurses.

Figure 3.

Performance metrics associated with the scheduling policies for Instance 5 with normal demand. (a) Cost; (b) Utilization; (c) Waiting time; (d) Deferral; (e) Transfer.

According to Figure 3, the Myopic policy performs the worst on most metrics, except for average utilization. The APPQ and ADMB policies show similar performance, never booking late but varying across different priority levels (e.g., waiting time). These policies struggle with transfer rates, transferring around 3% and 4% of the highest priority patients, respectively. They are designed for settings that allow overtime, enabling patient service in overtime to prevent transfers or late bookings. Without overtime, these policies provide lower-priority patients with more booking opportunities due to their longer wait time thresholds, reducing future capacity for higher-priority patients. This leads to higher transfer rates for high-priority patients. The AOP, on the other hand, defers and allows slightly late bookings, especially for higher-priority patients, to reduce transfer rates and total costs. Although late bookings might not be obvious from Figure 3(c) as the figure depicts average values, the AOP offers a slightly longer average wait time for Priority 1 than other policies. Surge capacity is a crucial resource for healthcare centers facing demand that exceeds their capacity (Demirbilek et al., 2019; Kaji et al., 2006). In our setting, transferring patients serves as a surge capacity measure, with lower-priority patients incurring smaller costs. The AOP intelligently uses this surge capacity for lower-priority patients, resulting in much smaller transfer costs.

Based on such observations, we developed the EPPQ and EDMB policies. These policies significantly outperform their counterparts, with EDMB performing very close to the AOP. However, they are more conservative in transferring lower-priority patients than the AOP, leading to unnecessary transfers, lower utilization rates, and higher total costs. A critical parameter affecting their performance is the transfer threshold. A larger transfer threshold may lead to inadequate transfers of lower-priority patients, causing high-priority patient transfers. Conversely, a smaller transfer threshold can reserve too much capacity for high-priority patients. Optimizing the transfer threshold could improve their performance.

5.5.2. Inferiority of the TBS Policy

Despite our partner hospital utilizing the TBS policy, our experiments revealed its poor performance. The TBS policy allocates fixed spots to different patient classes, alleviating the greedy use of resources seen with the Myopic policy. However, it fails to adapt to the dynamic nature of patient arrivals and disaggregates the decision-making process by separating capacity allocation from actual patient allocations. This separation reduces the flexibility of the optimization process and lowers solution quality (Aringhieri et al., 2015; Moosavi and Ebrahimnejad, 2020). Furthermore, the performance of the TBS policy heavily relies on two bounds that specify the number of permissible spots for each patient class (Constraint set EC.11). Evaluating the impact of these bounds is challenging due to the vast number of combinations and the computational expense involved. Therefore, selecting appropriate bounds is not straightforward.

5.5.3. Patient Preference for Campus

This research examines a hospital with nearby campuses, initially assuming that patients are indifferent to their allocated campus. Alternatively, we consider a scenario where patients have campus preferences. We adapt the heuristic scheduling rules to accommodate these preferences, allocating patients to a campus upon arrival based on a probability distribution reflecting their preference. As a result, all of a patient’s appointments must be scheduled at their chosen campus. In the E-companion, we compare the scheduling policies for scenarios where patients can choose their campus versus scenarios where they cannot.

According to our observation, all heuristic scheduling rules perform worse when patients can choose their campus. This is expected since allowing patient choice reduces the flexibility of the solution process. Although the percentage difference in performance between the best and worst policies narrows when patients can choose their campus, the absolute total increase in costs is more significant for underperforming policies. This means that while the relative performance gap decreases, the overall expenses rise more when patients have the option to choose their campus. In summary, it is preferable to make both campus allocation and appointment scheduling decisions simultaneously in a multi-campus setting, provided there are no strict constraints on campus allocation for patients, such as their preferences or medical eligibility.

5.5.4. Objective Function Coefficients

Objective function (20) incorporates four coefficients ranked by ascending magnitude: Deferral, waiting, ED use, and transfer costs. Initially, we assumed that the ED use and transfer costs were substantially larger compared to the deferral and waiting costs. However, in the E-companion, we explore a scenario where the ED use and transfer costs are significantly reduced. The results reveal that our proposed policies—EPPQ, EDMB, and AOP—maintain superior performance despite the narrowing of the gap difference between the policies. This finding underscores the robustness of our policies across different cost structures and highlights the critical importance of a strategic approach to managing high demand even when transfer costs are comparable to waiting time costs.

6. Conclusions and Remarks

This research studies a distributed multi-appointment AC advance scheduling problem with heterogeneous resources. We consider multiple patient classes and priority levels, where each class-priority is associated with a wait time threshold and a deferral threshold. Taking into account uncertain patient arrivals and use of the ED, we formulate the problem as an infinite-horizon MDP model. Since we cannot solve this model via conventional methods, we hybridize this model with a feed-forward NN to simplify feasibility constraints while respecting all assumptions. To alleviate the curse of dimensionality, we rely on an affine approximation architecture to approximate the value function and solve an equivalent linear programing model through CG to compute AOPs. We implement two necessary acceleration techniques, including dual stabilization and decomposition in the CG component.

Experiments demonstrated that the NN achieves over 95% accuracy on validation datasets for theoretical instances. We compare the AOP against three benchmark scheduling policies and two proposed heuristics (EPPQ and EDMB). The AOP performs the best, while the myopic policy performs worst, with EDMB closely matching AOP. Using realistic data from an Ontario AC center, we found that while their current TBS policy outperforms the myopic policy, it lags behind all others. EDMB and AOP met wait time thresholds for about 95% of patients, compared to 80% for TBS.

In summary, our work primarily contributes to the literature by introducing a hybridized ADP and NN algorithm. The methodology is versatile, with broad applicability to other optimization problems with complex constraints (e.g., machine scheduling, storage location assignment, crew scheduling, facility layout), where a component of the problem can be similarly partitioned. Moreover, the proposed heuristic scheduling rules not only outperform benchmark policies but also offer generalizable managerial insights that could instruct similar scheduling challenges.

An extension of our problem could explore integrating medical and operational decisions, such as determining oncology and hematology dosages to enhance health outcomes and patient convenience while managing limited operational resources. We utilize a fully-connected classifier NN to solve large-sized instances; however, this model may be unnecessarily complex for certain problems. Future research could examine alternative predictive models and assess the value of interpretability in these contexts. Additionally, our assumptions may not fully capture specific healthcare settings. Future work could address other uncertainties, such as service duration and punctuality, and investigate more flexible appointment frequency schedule beyond fixed weekly, bi-weekly, or monthly intervals.

Supplemental Material

sj-pdf-1-pao-10.1177_10591478251331143 - Supplemental material for Dynamic Distributed Ambulatory Care Scheduling

Supplemental material, sj-pdf-1-pao-10.1177_10591478251331143 for Dynamic Distributed Ambulatory Care Scheduling by Amirhossein Moosavi, Onur Ozturk and Jonathan Patrick in Production and Operations Management