Optimal appointment scheduling with a stochastic server : Simulation based K-steps look-ahead selection method

Article history: Received April 16 2017 Received in Revised Format July 22 2017 Accepted July 31 2017 Available online July 31 2017 This paper studies the problem of scheduling a finite set of customers with stochastic service times for a single-server system. The objective is to minimize the waiting time of customers, the idle time of the server, and the lateness of the schedule. Because of the NP-hardness of the problem, the optimal schedule is notoriously hard to derive with reasonable computation times. Therefore, we develop a simulation based K-steps look-ahead selection method which can result in nearly optimal schedules within reasonable computation times. Furthermore, we study the different distributed service times, e.g., Exponential, Weibull and lognormal distribution and the results show that the proposed algorithm can obtain better results than the lag order approximation method proposed by Vink et al. (2015) [Vink, W., Kuiper, A., Kemper, B., & Bhulai, S. (2015). Optimal appointment scheduling in continuous time: The lag order approximation method. European Journal of Operational Research, 240(1), 213-219.]. Finally, a realistic appointment scheduling includes experiments to verify the good performance of the proposed method. © 2018 Growing Science Ltd. All rights reserved


Introduction
The appointment scheduling problem (ASP) to a stochastic server has been well known and widely studied in the literature since 1950s.The problem is consisted of a finite set of customers with stochastic service time and the service time of the customers is commonly assumed to be distributed with a known distribution or at least some historical data.The problem aims to select a deterministic schedule of appointment times to optimize competing performance criteria, such as waiting time of customers, the idle time of the server, and the overtime of the schedule.Commonly, the scheduling is finite, which is limited by the number of customers seen on a particular day.Second, customers must arrive punctually according to a defined schedule of appointment times, which is reasonable because empirical evidence suggests that patients arrive early more than late.
The problem we address is similar to the problem in Vink et al. (2015).Vink et al. (2015) studied a problem of scheduling a finite set of customers with stochastic service time for a single-server system.In addition, they proposed a lag order approximation method in which the lag order refers to the number of predecessors taken into account.Contrary to Vink et al. (2015), we will propose a simulation based K-steps look-ahead selection method in which K successors are taken into account to optimize the current decision.In addition, the biggest difference between these two methods is that the lag order approximation method is a forward method while the simulation based K-steps look-ahead selection method is a kind of backward induction.
In ASP, the service provider is able to schedule arriving customers with the help of an appointment schedule, which is consisted of a series of appointment times.Then the customer will arrive at the appointment time which is the earliest time he can receive service.In order to simplify the problem, an ASP can be divided into two stages: the service provider gives an appointment schedule in the first stage and the server executes the service in the second stage.In reality, we can imagine that the customers present themselves in random order at the first stage, which is reasonable if we assume that the sequence of customers is according to the rule of First Come First Served (FCFS).Then the service provider just need to decide on how to choose the appointment times of customers that are to be scheduled to the server.
The basic setting in this paper belongs to the so-called static ASP, in which the appointment schedules are scheduled prior to the beginning of the actual service (Cayirli & Veral, 2003).Assume that there is one server in the system who works whenever there is a customer in the system and is idle, otherwise.We also assume that there are customers that need to be scheduled at time zero.Furthermore, suppose that the service-time distribution (or at least some historical data) and customer' utility function due to waiting time, as well as the server's utility function, in terms of idle time and possible lateness after the final client (overtime), are known.The objective is then to minimize a convex combination of, possibly weighted, sum of the customer's waiting time, the server's idle time and lateness.
The problem is very important and widely studied, and to the best of our knowledge, there are few clean solution strategies in previous literature.Vink et al. (2015) proposed a lag order approximation method which can solve the problem well.In this paper, we explore another alternative approach.Our approach is able to deal with general service-time distributions, even the case that distribution is unknown with just some historical data.In our approach the optimal appointment times depend only on a limited number of clients that arrived subsequently to the customer's appointment, leading to an optimization problem with reduced dimensionality.For instance, we optimize the schedule for customers , 1 and 2 with consideration of their expected waiting times, corresponding idle times, and lateness of the server.Then according to the optimal solution of these three customers, we decide the appointment time of customer .We refer to this method as K-steps look-ahead selection method in which successors are taken into account to optimize the current decision.In each iteration, the subproblem will be solved by a simulation approach, therefore, the total approach is called simulation based K-steps look-ahead selection method (SBKSLS).
The remainder of this paper is organized as follows.Section 2 is a brief review of the related literature.Section 3 mathematically formulates the problem.In section 4, we present SBKSLS.The performance of the proposed method is evaluated in Section 5 by studying some numerical examples.Finally, Section 6 concludes and discusses directions for further research.

Literature review
ASP has been widely studied and this topic has attracted the interest of many academicians and practitioners over the last 60 years.Many of the studies have been researched in outpatient clinics and other health-care environments, therefore, we use "patient" and "customers" interchangeably in the following review.Bailey (1952) is the first paper to study ASP.For a comprehensive overview on ASP, we refer to the review works given by Cayirli and Veral (2003) and Gupta and Denton (2008).Cayirli and Veral (2003) provided a comprehensive survey of research on appointment scheduling in outpatient services and identified future research directions that provided opportunities to expand existing knowledge.Gupta and Denton (2008) summarized key issues in designing and managing patient appointment systems for health services and exposed open research areas and opportunities for future work.
One of the main issues in appointment scheduling is the different distributions of the service times.Most of the contributions on appointment scheduling are based on exponential service times, such as in Wang (1999), Wang (1997) and Kuiper et al. (2015) studied the problem with a phase-type distribution for the service times.In the previous papers, it is common to assume independent and identically distributed random variables for the service times.It seems reasonable to assume that the service-time distribution of the customers are independent, since customers call in randomly for an appointment.However, in reality the service times often do not follow an exponential distribution, let alone the service-time distributions of the arriving clients are identical.Thus, more and more papers studied the problem with looser requirements for the distribution of service times.Thus, Chakraborty et al. (2010) studied the problem for patients with general service time distributions.Mak et al. (2014) developed distributionfree models that solve the appointment sequencing and scheduling problem by assuming only moments' information of job durations.Begen et al. (2012) considered the problem of appointment scheduling with discrete random durations but under the more realistic assumption that the duration probability distributions are not known and only a set of independent samples is available, e.g., historical data.In our paper, the method proposed can not only solve the known distribution, but also the distributions are assumed unknown, that only a set of independent samples is available, e.g., historical data.
There are numerical papers studied research methodologies about ASP.Weiss (1990) saw the situation as a D/GI/1 queueing system, where the arrival times of the customers are a decision variables.Wang (1997) used the phase-type distribution functions and matrix algebraic manipulation to obtain recursive expressions for the customer flow-time distributions and developed a computational procedure which can efficiently evaluate the mean flow-times for large number of customers that need to be scheduled.Wang (1999) showed that the service order depends on the order of service rates and the optimal schedule can be obtained by solving a set of nonlinear equations if the service times are exponentially distributed with different rates.Denton and Gupta (2003) presented that the ASP can be expressed as a two-stage stochastic linear program and a standard L-shaped algorithm is employed to obtain optimal solutions.Robinson and Chen (2003) used the structure of the optimal solution as the basis for a simple closedform heuristic for setting appointment times and the heuristic is shown to perform on average within 2% (and generally within 0.5%) of the optimal policy.Bendavid and Golany (2011) employed the Cross-Entropy method to solve the problem, which aimed to set gates so to minimize the sum of the expected holding and shortage costs.Mancilla and Storer (2012) developed a heuristic solution approach based on Benders' decomposition for a single-resource stochastic appointment sequencing and scheduling problem with waiting time, idle time, and overtime costs.Kuiper et al. (2015) demonstrated how to optimally generate appointment schedules and in their procedure, they replaced general service-time distributions by their phase-type counterparts, and then optimize a utility function.Kemper et al. (2014) indicated that rules were needed to assure a good trade-off between quality (in terms of the customer's waiting time) and cost (in terms of the server's idle time) and this paper presented a technique to generate such rules.Kemper et al. (2014) also suggested that one should schedule jobs in the order of increasing variances, for convex loss functions with scale families of service time distributions.
Although there are so many papers studied the approach to solve the ASP, however, there is few clean solution available.Vink et al. (2015) proposed a lag order approximation method in which the optimal appointment times depend only on a limited number of customers that arrived previously to the customer's appointment, leading to an optimization problem with reduced dimensionality.In our approach, the optimal appointment times also depend only on a limited number of customers leading to an optimization problem with reduced dimensionality.Contrary to Vink et al. (2015), we used a backward induction, that the current decision depends on a limited number of customers arrived subsequently to the customer's appointment.In each iteration, the subproblem will solved by a simulation approach, therefore, the total approach is called as SBKSLS.

Problem description
Assume that there is one server in the system who works whenever there is a customer in the system and is idle, otherwise.In addition, we assume that there are customers, denoted by 1,2, . . ., that need to be scheduled at time zero.Each customer has a stochastic service time, which is denoted by the random variable for customer .The service system has a single server and if upon arrival customer finds the server idle, he immediately starts his service.If the server is busy, then customer awaits his turn until all customers that are scheduled before customer have finished their services.We also assume that all the customers are punctual, and no-show and walk-in customers are not allowed.
The vector ( , , . . ., ) is denoted as an appointment schedule for this service system.For a given schedule, denotes the time that server has been idle upon starting the service of customer .In addition, we denote by the waiting time of customer .It is noted that the sojourn time of customer can be defined by Eq. (1).
The planning horizon , in which customers can be scheduled, is finite.However, it can happen that after the planning horizon there are still customers who need to be served.Therefore, we define the lateness as the overtime that the server has to be made in order to finish all services.The interappointment time, which is the time between customer and 1 arrival, is defined by Eq. ( 2).
The idleness is given by Eq. ( 3).
The waiting time is given by Eq. ( 4).

max
, 0 Clearly, it is reasonable to assume that 0, so that both 0 and 0.Moreover, it holds that 0, 1,2, . . ., .The objective of ASP is to find a schedule ( , , . . ., ), or equivalently , . . ., , such that a utility function , which depends on , and , is minimized.Throughout this paper, we assume that has the form as shown in Eq. ( 6).

Utility functions
In this subsection, we present two utility functions that are commonly used in the literature.The utility function includes the expected waiting times, the expected idle times, and the expected lateness with different weighing factors.
One often chooses general polynomial function and sets , and , where , , 0 and 0. However, setting and taking 1,2 gives us two remarkable insights, which we will refer to as the absolute value utility function and the quadratic utility function.Note that in these cases the idle and waiting times are equally weighted.

Absolute value utility function
The absolute value utility function can be obtained by taking and , with , ∈ . We know from Eq. ( 6) that the utility function reduces to This utility function penalizes deviations from the planning (either caused by waiting or by idling) linearly.It has been used (with 0) by, e.g., Wang (1999) and Kuiper et al. (2015).

Quadratic utility function
The quadratic utility function UF penalizes the deviation from the schedule quadratically instead of linearly.This can be achieved by taking.and .Since for 1, the utility function reduces to This utility function has been used ( 0) by Kemper et al. (2014).

Simulation based K-steps look-ahead selection method
In this section, we introduce SBKSLS in its general form.Basically, the optimal schedule is found through the optimization of Eq. ( 6), that is min ,..., , . . ., The waiting time of customer is a random variable depending on , . . ., .The main idea of the proposed method is to neglect part of the customers that influence the waiting time (and idle time and lateness) of the utility function in Eq. ( 6), and obtain the schedule for each customer in terms of successors.

Procedure for SBKSLS
This subsection introduces the procedure on obtaining a good schedule if the sequence of customers is given.This procedure tries to insert all the customers into time axis one by one according to the sequence.The procedure consists of iterations for inserting the customers.
Let the inserting sequence of customers be 1,2, . . ., .For customer , we solve the subproblem just for customers , . . ., with the objective function which is shown in Eq. ( 10).min ,..., , . . ., Then we can determine the value of appointment time for customer which is denoted by * .For customer , based on the obtained appointment times for customers , . . ., , i.e., * , . . ., * , we solve the subproblem just for customers , . . ., with the objective function which is shown in Eq. ( 11).min ,..., * , . . ., * , , . . ., It is noted that, for the th iteration, we can obtain the value of appointment time for customers , . . ., .Therefore, the problem is divided into subproblems which is easily solved.The method for solving the subproblem is presented in Section 4.2.and the procedure for SBKSLS is presented in Algorithm 1.

Algorithm 1
The procedure for SBKSLS Input: The service times of a given sequence customers.Output: The appointment schedule.

Simulation approach for solving subproblem
In this section, we will introduce a simulation approach to solve the subproblem in each iteration.Begen and Queyranne (2011) considered the problem of determining an optimal appointment schedule for a given sequence of jobs on a single processor to minimize the expected total cost and proved that there exists an optimal appointment schedule that is integer and can be found in polynomial time.Here we also divide the time into several identical parts.
Let and denote the minimum and the maximum possible values of appointment times.These two values can be obtained by historical experiences or results of preliminary experiments.Then we divide into parts and per unit length is equal to .Here is a pre-defined number.This is reasonable since that minute is the smallest length of time in reality.Thus the searching space for customer is defined as: | , 0,1, … , .The searching space is a discretized area.The lattice size is pre-defined according to the scale of the problem and the computer's capacity.Therefore, there are 1 possible solutions for each subproblem.For each solution, we use a simulation approach to obtain the value.The procedure is as follows: in a set of random numbers of a particular distribution, one minimizes the utility function over the interappointments * , . . ., * , , . . ., .Here is also a pre-defined number.Then compare the value of 1 possible solutions and select the smallest one as the final solution.Thus, we can solve the subproblem by a simulation approach.

Monte Carlo Simulation for generating UF values
In order to simulate different distributions we apply a similar approach which uses random numbers in a Monte Carlo Simulation study.The procedure is as follows: in a set of 10 random numbers of a particular distribution, one minimizes the utility function over the interappointments , . . ., .We repeat this 100 times so in total 10 random numbers are needed to get a sample mean, and a sample variance, , of system's cost.Next we illustrate the results of the simulation approach for a system with 11 customers.We assume that the service times of the customers are independent and exponentially distributed (i.i.d.) with parameter 1 for 1,2, . . ., .The results are shown in Fig. 1 and the red line represents z .
. We can find that almost all the points appear between two lines.Therefore, the simulation approach can obtain a UF value which can reflect the performance of schedule well.

Complexity analysis
Here we discuss the complexity of the proposed heuristic algorithm.First we analyze the procedure of solving subproblem.The complexity of this procedure is 1 .Based on the above, we return to the general procedure of the proposed heuristic algorithm.Its complexity is 1 .Therefore, the complexity of the proposed algorithm is 1 .

Numerical experiments
In this section, we will determine the key parameters in SBKSLS and apply the proposed method to a realistic appointment scheduling problems.We will determine the key parameters in Section 5.1.Then in Section 5.2, we will study different distributed service times, e.g., Weibull and lognormal distribution.In Section 5.3 we will apply the proposed method to a real-life example of a CT-scan in a hospital.

Determine the parameters of SBKSLS
It is necessary to investigate the computing speed (CPU time) of the proposed algorithm.The CPU time is influenced by many aspects, e.g., the number of customers ( ), the amount of customers for the each iteration ( ) and the amount of random numbers in each iteration ( ).However, the first factor is the character of the problem and it is determined by the scale of the instances.Then and are two important parameters in the proposed approach and it needs to find an appropriate combination of and for the proposed method which can consider both the quality of the solution and the CPU time.
First, we conduct an experiment to determine parameter and the instances are generated which can refer to Vink et al. (2015).We illustrate the results of the proposed method for a system with 11 customers in Fig. 2. We assume that the service times of the customers are independent and exponentially distributed (i.i.d.) with parameter 1 for 1,2, . . ., .The system operates under linear utility and quadratic utility function with 1 and 0, respectively.For the problem, we solve it with five different (20, 50, 100, 200 and 500) with a given 2 and each instance is solved 10 times.We calculate the average values of the instances as to these combinations ( and ), respectively and the results are shown in Table 1.The first two columns under the parameter category display the values of and , respectively.The rest of Table 1 are divided into two categories: the results for linear utility functions and quadratic utility functions.In each case, the value and CPU time are displayed the resulting average (A), largest (L) and smallest (S) values.The variance of the value ( ) is also displayed.
From Table 1, we observe that the value all have a tendency to decrease with the growth in value and the CPU time has an opposite trend.Besides, the CPU time is acceptable for all instances.Thus, we set 500 for the following experiments because of the best values and the smaller variance for 10 times, which means the value has a better stability.
Second, an experiment is conducted to determine parameter .The instances are generated as the same as the previous experiment and for this problem, we solve it with three different (0, 1 and 2) with a given 500.Table 2 summarizes the results of SBKSLS and the other methods.It is noted that, all the results of lag order approximation method and simulation in this paper are from Vink et al. (2015).We can observe that increasing value reduces the expected total cost of the system.In addition, the results obtained by our proposed method are better than these obtained by lag order approximation method, which is proposed by Vink et al. (2015), on both the objective and the CPU time.When 2, the gap is even smaller than 1% between our results and the optimal results obtained by simulation.According to above two experiments, we set 2 and 500 for the following experiments.Fig. 2 displays the optimal sizes per arriving customers for different values.We can find that when the value increases the approach results in larger slot sizes.For example, for the first slot (the time between the first and second appointment) the slot size increases from 1.11 based on 0-step look-ahead to a slot size of about 1.35 based on II-steps look-ahead.Furthermore, for the result of II-steps look-ahead, the slots in the middle of the scheme are larger than the slots in the beginning and the end of the scheme.We can also find that the results obtained by II-steps look-ahead is very close to the optimal results which is obtained by simulation method.

Performance analysis on different distributed service times
In this section, we will study the behavior of different distributed service times, e.g., Weibull and lognormal distribution.According to Cayirli and Veral (2003) service-time distributions have a typical coefficient ranging from 0.35 to 0.85.Given this range, we illustrate our proposed method by using Weibull and lognormal distribution with reasonable coefficient of variations: 0.35, 0.5 and 0.85.In Tables 3 and 4, we display the results of various values to this problem.The results are compared to the lag order method proposed by Vink et al. (2015) and an optimal schedule derived through simulation.Fig. 3 Optimal lot size for SBKSLS compared with lag order method and simulation.11, Weibull distributed service times with mean 1 and CV=0.5 From the Table 3 and Table 4, we conclude that the application of the proposed method results in a small loss of quality of the appointment schedule.In any case, linear or quadratic utility and either for Weibull or lognormal service-time distribution, our approach with 0 step generates schedules that are no more than 12% from the optimal UF value.Furthermore, from the tables we see that the II-steps look-ahead generates schedules that are reasonably close to the optimal UF value (around 0.1-1.0%).In addition, the computation time is also smaller than that of lag order method for all cases.In Fig. 3 and Fig. 4 we illustrate the schedules derived by the SBKSLS, lag order II method and the simulated optimal in both a Weibull and lognormal setting with coefficient of variation equal to 0.5.We can find that the results

Fig. 2 .
Fig. 2. Optimal lot size for SBKSLS compared with lag order method and simulation with quadratic utility.11, i.i.d.exponential ( 1) service times

Fig. 4 .
Fig. 4. Optimal lot size for SBKSLS compared with lag order method and simulation.11,i.i.d.lognormal distributed service times with mean 1 and CV=0.5

Table 1
Optimization results of SBKSLS with different amount of random numbers for 2

Table 2
Optimization results of SBKSLS compared with lag order method and simulation.

Table 4
Optimization results of SBKSLS compared with lag order method and simulation.