Staffing for many-server systems facing non-standard arrival processes
SSTAFFING FOR MANY-SERVER SYSTEMSFACING NON-STANDARD ARRIVAL PROCESSES
BY M. HEEMSKERK, M. MANDJES & B. MATHIJSEN
Abstract.
Arrival processes to service systems often display (i) larger than anticipated fluctuations,(ii) a time-varying rate, and (iii) temporal correlation. Motivated by this, we introduce a specific non-homogeneous Poisson process that incorporates these three features. The resulting arrival process isfed into an infinite-server system, which is then used as a proxy for its many-server counterpart.This leads to a staffing rule based on the square-root staffing principle that acknowledges the threefeatures. After a slight rearrangement of servers over the time slots, we succeed to stabilize systemperformance even under highly varying and strongly correlated conditions. We fit the arrival streammodel to real data from an emergency department and demonstrate (by simulation) the performanceof the novel staffing rule.
Keywords.
Queueing; Applied Probability; OR in Health Services.
Affiliations.
Mariska Heemskerk and Michel Mandjes are with Korteweg-de Vries Institute forMathematics, University of Amsterdam; Science Park 904, 1098 XH Amsterdam; The Netherlands( email : j.m.a.heemskerk|[email protected] ).Britt Mathijsen is a former PhD candidate at the Department of Mathematics and ComputerScience; Eindhoven University of Technology; P.O. Box 513, 5600 MB Eindhoven; The Netherlands.( email : [email protected] ). Version: June 2, 2020. Funding.
The research of the first two authors is funded by the NWO Gravitation ProgrammeN etworks (Grant Number 024.002.003). The research of the second author is partly funded by anNWO Top Grant (Grant Number 613.001.352). The research of the third author was funded by anNWO Free Competition Grant (Grant Number 613.001.213). Introduction
The design of staffing algorithms for service systems has been attracting a great deal of interest eversince Erlang published his first papers. The delay experienced by the system’s users is predominantlycaused by the inherent randomness in the arrival stream. From a managerial point of view, it isnatural to address this randomness in such a way that operational costs and customer satisfactionare balanced. An important complication that recently received a lot of attention is that, as hasbeen observed in various empirical studies, the variance of the arrival stream is larger than thecorresponding mean; this phenomenon, called overdispersion, is not captured by standard Poissonprocesses. The challenge that arises is to develop staffing algorithms that are based on moresophisticated, realistic arrival stream models. See [9, 14, 19, 27, 35] for related work on the design ofstaffing rules for service systems with overdispersed input. Such staffing rules have broad application a r X i v : . [ c s . PF ] M a y BY M. HEEMSKERK, M. MANDJES & B. MATHIJSEN potential, in settings that include call center environments [1, 3, 5, 7, 40], cloud computing [34, 36],and health care delivery [28].
Staffing in stochastic service systems.
Due to the intrinsic stochastic variability of most servicesystems, they are naturally described by a stochastic model. More specifically, queueing models haveproven to reveal ‘good’ staffing rules; they determine the number of staff needed to effectively butefficiently cope with the demand imposed on the system. Evidently, any staffing rule should be suchthat the average workload brought in per time unit is smaller than the system capacity, to make surethat the delays experienced by patients remain bounded. Moreover, certain performance targets areset to guarantee the patients a specific ‘Quality-of-Service’ (QoS) level, which is typically expressedin terms of waiting time. The objective of the system operator is, on the other hand, to bring theutilization level as close as possible to , so as to control operational cost by efficiently using theresources. Staffing rules aim to strike a proper balance between the interests of the patients and thesystem operator.The goal that is often strived for is to design a service system in such a way that its patients go intoservice more or less immediately upon arrival [19, 40]. Consequently, a commonly chosen service-levelagreement (SLA) is to bound the probability of delay by some (typically small) QoS parameter ε > [5, 18, 42]. In the typical situation that the arrival rate varies over time, this delay probability isclearly time-dependent too. The manager’s objective would be to set up a staffing schedule thatminimizes the operational costs under the constraint that the delay probability, at any point of time,does not exceed ε . The ideal staffing algorithm is one that stabilizes the probability of delay overtime around some value close to ε , bringing the system in the so-called Quality-and-Efficiency-Driven(QED) regime [13]. Note that, in case there are periods where the algorithm induces a probabilityof delay significantly smaller than ε , it would mean that (at least locally) ‘too many’ resources aredeployed; the variability in the arrival stream is anticipated suboptimally. Realistically modeling the arrival process.
Queueing models have been used intensively to describeand understand congestion phenomena that arise in case of scarce resources. As a first step indesigning a realistic model, it is key to study the arrival process at hand. We will continue thisintroduction with a short recap on the properties that should be present in a realistic arrival streammodel.A common assumption in queueing theory is that of Poissonian arrivals, entailing that the mean andvariance of the number of arrivals (roughly) match. However, it is often observed that service systemsface arrival streams that are highly variable (mean (cid:28) variance; overdispersion), while in specificcases systems have to deal with almost deterministic arrivals (variance (cid:28) mean; underdispersion).As an example of the latter, consider service systems in healthcare with scheduled yet not necessarilypunctual arrivals (so that arrival epochs randomly fluctuate around the appointed arrival time), asstudied in e.g., [21, 22, 25]. In such settings clearly some sort of ‘induced deterministicness’ plays a role, in the sense that arrivals are actively being directed to (or away from) the system. In thisthesis however, we will focus on ‘undirected’ arrival streams only.For ‘undirected’ arrival streams, overdispersion is a phenomenon commonly found in data. Exampleswhere one could expect to encounter overdispersed arrivals include a call center of a bank, an insurancecompany and an emergency department in a hospital; see e.g. [3, 23, 24, 27]. In such settings,arrivals are usually triggered (or inhibited) by occasional events or (un-)favorable circumstanceswhich can cause unforeseen peaks (or dips) on top of the usual daily patterns. This so-called ‘randomenvironment’ gives rise to an effect commonly referred to as parameter uncertainty [3, 40], whichnaturally leads to overdispersion.Speaking of daily patterns: in nearly all practical applications, the mean number of arrivals is notconstant over time (e.g. over the course of the day) and follows a predictable pattern. It must benoted that the variability that causes overdispersion is of a different nature than the variabilityinduced by nonstationarity. Nonstationarity can be modeled by a non-homogeneous Poisson proces,replacing the constant arrival rate of a Poisson process by a (deterministic) time-varying one.However, for non-homogeneous Poisson processes the mean and variance of the number of arrivalsstill match, hence such processes fail to capture the entirety of the desired dynamics observed inarrival processes. Nevertheless, nonstationarity is another important feature of a real-life arrivalprocess [10, 11, 37] and as such should be incorporated in any realistic arrival stream model as well.Besides being overdispersed and having a time-varying rate, a realistic arrival stream might evenhave dependencies between the numbers of arrivals in disjoint time intervals. That is to say: it’shighly unlikely that the random environment affects the arrival stream in an i.i.d. fashion over thedifferent intervals; the effects at hand possibly play a role for a longer period of time. Indeed, arrivaldata often exhibits these kinds of dependencies, e.g. in call centers [16, 17].
Existing staffing methods.
As mentioned, our objective is to develop a staffing rule such that thedelay probability is sufficiently low, uniformly over time. With this rather stringent service-levelrequirement in mind it is fairly natural to approximate the system relying on its infinite-servercounterpart. The famous square-root staffing principle is based on exactly this observation. Inthe classical setting (M/G/ ∞ with arrival rate R and unit-mean service times) it uses that thesteady-state number of busy servers is Poisson distributed with mean R . By asymptotic normalitythe coefficient of variation (i.e., standard deviation divided by the mean) has the approximate form / √ R , and that the steady-state probability of delay in a corresponding finite-server setting with s servers, say p s ( R ) can be approximated by p s ( R ) ≈ − Φ (cid:18) s − R √ R (cid:19) BY M. HEEMSKERK, M. MANDJES & B. MATHIJSEN for large R , with Φ( · ) the distribution function of a standard Normal random variable. For β suchthat − Φ ( β ) = ε we find for the required number of servers [38, 39] s = R + β √ R. (1)Note that Eqn. (1) has an appealing interpretation: the number of servers s should evidently betaken larger than the expected workload R , with the extra term β √ R (‘uncertainty hedge’) beingof the same order as the natural load fluctuations of the workload process. Refinements of ordersmaller than √ R are explored in e.g. [18, 29, 42].Although the excess probability P ( M > s ) corresponding to the infinite-server system (with M denoting the stationary number of busy servers in this M/G/ ∞ queue) is likely to be smaller thanthe probability of delay in its finite-server counterpart, still square-root staffing rules have shown togive accurate results [5, 18]. This can be explained by the fact that as R grows large, the hedge β √ R prevents congestion more and more effectively.So far we discussed the situation of a constant Poissonian arrival rate. For large-scale systemsthe predominant assumption in the literature is that patients arrive according to a time-varyingPoissonian arrival rate. Staffing algorithms for non-homogeneous Poisson processes (NHPPs) havebeen studied for several decades.If the arrival process is an NHPP with nonstationary arrival rate λ ( · ) , then the number of arrivals N ( s, t ) in the interval [ s, t ) , with s < t , is Poisson distributed with parameter (cid:90) ts λ ( r ) d r. Note that such a non-homogeneous arrival process not yet exhibits overdispersion ( E N ( s, t ) = V ar ( N ( s, t )) ). For the resulting model M/G/ ∞ -based staffing rules cannot be applied directly, butvarious approaches have been proposed.In a first approach, the nonstationarity is essentially ignored: one uses a simple stationary approxi-mation (SSA), based on a stationary model in which the arrival rate is chosen equal to the long-runaverage [32]. This method performs poorly in most scenarios [11], for example when the actual rateis slowly varying with respect to the service time or when the relative amplitude (level of nonstation-arity) of the rate is larger than 10%. A second approach is the pointwise stationary approximation(PSA) [10, 11, 37], which considers the system at time t as if it has dealt with an arrival rate λ ( t ) with s t servers from the start (i.e., assuming steady state), thus ignoring nonstationarity in a differentway. This method works well in settings where the arrival rate changes sufficiently slowly [10, 37], soit covers scenarios on the other side of the spectrum.As a comprimise between the two extremes, [37] suggests the average stationary approximation(ASA) that generalizes both SSA and PSA, replacing the arrival rate at time t with ¯ λ t = E Sa (cid:90) tt − a/ E S λ ( r ) d r for some positive constant a and mean service time E S . An alternative to this was proposed by [19],saying that one could replace R in the staffing formula by m ∞ ( t ) , the expected ‘offered load’ in aninfinite-server system with time-dependent arrival rate λ ( t ) at time t : m ∞ ( t ) = E (cid:90) tt − S λ ( r ) d r, where S denotes the service time. They showed that this method stabilizes the probability of delayclose to some target value at all times, independently of the arrival rate being slowly varying ornot. Importantly, in [19] asymptotic normality was used to arrive at the approximation, hence theirmethod follows the tradition of the square-root staffing procedure described above.As mentioned above, NHPP models fail to exhibit overdispersion, a phenomenon observed acrossvarious types of service systems; see e.g., [6, 20, 23, 31, 33]. The parameter uncertainty underlyingoverdispersion potentially jeopardizes the effectiveness of the square-root rule, typically leading tooveroptimistic staffing algorithms. This complication was brought forward in many studies, e.g.in [2, 3, 7, 8, 12, 20, 23, 29, 30, 33, 41]. Different methods were proposed to overcome this issue;typically a Poisson mixture is used to model parameter uncertainty. That is, the deterministicPoissonian arrival rate is replaced by a sampled one; see [8, 40, 6, 20, 4, 28, 26, 29] for examples.Relatively little attention has been paid to staffing rules in the context of arrival processes in whichthe numbers of arrivals in disjoint intervals are dependent.
Contributions and organization.
The contributions of this paper are twofold. In the first place, wepresent a flexible model for the arrival process, based on [15], that can deal with (any level of)overdispersion, nonstationarity and dependencies between arrivals of consecutive time slots. Thechallenge being to come up with a model that remains of practical use/computationally tractable,we believe that the model proposed here is among the simplest models with these three properties.Moreover, fitting data to our model is a relatively straightforward task. The model is presented inSection 2.1In the second place, we develop staffing rules meeting the criteria mentioned in the introduction, togo with this comprehensive yet simple model for the arrival stream. It requires low computationalcost to determine staffing prescriptions based on these rules. In Section 2.2 we present the newstaffing rule. Subsequently, in Section 2.3 we present a case study based on a healthcare-relateddata set to show that the rule succeeds to stabilize the delay probability around the targeted ε . Theobservations here lead to a much improved version of the staffing rule that was introduced in Section2.2. This concludes Section 2.In the rest of the paper we work out the details necessary for implementation and further asses theperformance of the proposed staffing rules. Section 3 presents straightforward statistical proceduresto estimate the modeling parameters. We perform the estimation procedure both for a real dataset from an emergency department, and for a stylized example. In Section 4 we perform extensive BY M. HEEMSKERK, M. MANDJES & B. MATHIJSEN experiments to assess the performance of our staffing rule in settings with overdispersion, a time-varying arrival rate, and temporal correlation. Here, we incorporate impatience into the model (asin reality, customers might abandon the system before their service begins), in order to analyze howthis affects the performance. Finally, Section 5 concludes the paper.2.
Model and staffing rule
In this section we first present our arrival stream model meeting the criteria mentioned in theintroduction (overdispersion, time-varying rate, correlation between disjoint time intervals). Ourmodel is arguably the simplest among all models satisfying these requirements. It is relativelycompact and only requires a few input parameters. We then introduce a suitable staffing rule tomatch such arrival streams. It is new compared to the earlier described methods in the introduction,as it combines all three features of realistic arrival processes while still using the concept of square-root staffing, where the mean under the square-root is replaced by the variance of the number ofcustomers in the approximative infinite-server system. We conclude the section by an illustrativecase study, in which we demonstrate the procedure and its performance.2.1.
Model description.
The model we consider could be termed a mixed M t /G/s t queue withinfinite waiting room. We systematically introduce the components of the model, starting with thearrival process. Arrival process . In our setup the arrival process is a Cox process, i.e., a time-dependent Poissonprocess with random arrival rate. At time t , the arrival rate is Λ( t ) (cid:62) . This Λ( t ) consists of a deterministic trend λ ( t ) (capturing the daily pattern), which is inflated by a stochastic busynessfactor that incorporates the desired overdispersion and temporal correlation. As a consequence, themodel proposed possesses the three desired properties.More specifically, the arrival rate is built up as follows. Following common practice, we assume that λ ( t ) is piecewise constant on time intervals of fixed size ∆ . For t ∈ [ j ∆ , ( j + 1)∆) we can thereforewrite λ ( t ) = λ j . We introduce a sequence of random variables W ≡ { W j } j ∈ Z , which are independentand distributed as a random variable W (cid:62) ; we normalize them such that E W = 1 , and assume V ar ( W ) < ∞ . The busyness factor of slot j is affected by the current value of the W process ( W j ,that is), but also by the previous I values ( W j − I up to W j − ). The parameter I ∈ N reflects theamount of dependence between the stochastic arrival rates pertaining to consecutive disjoint slots.Let N be the total number of time slots of size ∆ in the considered time frame, i.e. N = 24 when ∆ = 1 hour and the considered time frame is a day. The level of dependence from previous values ofthe W process is dealt with in an autoregressive way, with parameter α ∈ (0 , . Concretely, thismeans that for t ∈ [ j ∆ , ( j + 1)∆) the stochastic arrival rate is given by Λ( t ) = λ j · (cid:16) c α I (cid:88) (cid:96) =0 α (cid:96) W j − (cid:96) (cid:17) ; (2) here c α := (1 − α ) / (1 − α I +1 ) is a normalizing constant that ensures that the busyness factor hasmean 1: E (cid:32) I (cid:88) (cid:96) =0 α (cid:96) W j − (cid:96) (cid:33) = 1 − α I +1 − α = 1 c α . It means that the busyness factor gets a new value every ∆ time units; thus, ∆ − can be regarded asthe sampling frequency . Note that the process W is not observable; as we show later, in our staffingformula we just need to know V ar ( W ) .The values λ j reflect the mean arrival rates during the individual time slots. When assumingperiodicity in the data (e.g., daily and weekly patterns), one can estimate these values in astraightforward way from historic data. The value of ∆ is situation-dependent; one often picks , or minutes. This leaves us with estimating α , I , and V ar ( W ) . The procedure we followed is thatwe use standard least-squares tools to estimate α and V ar ( W ) for given I ; this we do for multiplevalues of I , so as to select an ‘optimal’ I . We elaborate on these estimation issues in Section 3. Service times.
The patients’ service times are independent and identically distributed samples froma general non-negative distribution; we denote the underlying random variable by S , and write P ( t ) := P ( S > t ) . In the numerical experiments in Sections 2.3 and 4 we focus on the case ofexponentially distributed service times (with mean µ − ), but in the staffing rule one could pick inprinciple any distribution. Number of servers.
At time t , the number of servers is s t . The value of s t is as determined in Section2.2. We assume that services are always completed, even if s t drops to a value that is insufficient toserve all patients present; as this assumption is fairly natural in practice this is the way the systemdynamics will be modeled in the simulation experiment.2.2. Staffing rule.
The staffing rule we propose is essentially an adaptation of the classical square-root staffing rule in Eqn. (1): for some constant β > , s t = m ∞ ( t ) + β (cid:112) v ∞ ( t ); (3)here m ∞ ( t ) and v ∞ ( t ) are the mean queue length and variance of the mixed M t /G/ ∞ counterpartof the mixed M t /G/s t system introduced in Section 2.1. Note that given an overdispersed arrivalstream, the term β (cid:112) v ∞ ( t ) (the hedge) is larger than in the classical SRS rule, where it would equal β (cid:112) m ∞ ( t ) .The use of such a rule is justified by asymptotic normality, which is backed by the results in [15].The initial choice for the constant β is (with Φ( · ) the normal CDF): β = Φ − ( ε ) . (4)It is expected that this choice is not optimal, given the approximative nature of the procedure. Infact, β is always smaller than optimal, since the actual number of customers in a finite-server systemwill be higher than predicted by an infinite-server proxy (where each customer is served immediately BY M. HEEMSKERK, M. MANDJES & B. MATHIJSEN upon arrival and hence can leave without waiting). The idea is to slightly tweak the value β in orderto more closely attain the desired service level (i.e., P ( delay ) < ε ). More importantly, irrespective ofthe level the shape of s t as a function of time should ensure a delay probability that is relatively flatover time. This depends mostly on the shape of m ∞ ( t ) and v ∞ ( t ) .Hence, the next step is to determine expressions for m ∞ ( t ) and v ∞ ( t ) . Following the approachof [15], we deduce that the queue-length process of an infinite-server queue fed by a Cox processarrival process with arrival rate Λ( t ) is again a Cox process. More specifically, the distribution ofthe number of patients at time t is Poisson with random parameter R ( t ) = (cid:90) t −∞ Λ( s ) P ( S > t − s ) d s. (5)We thus obtain that m ∞ ( t ) = E R ( t ) = E (cid:20)(cid:90) ∞ Λ( t − s ) P ( S > s ) d s (cid:21) (6)and, by the law of total variance (conditioning on Λ( s ) , for s ∈ ( −∞ , t ] ), v ∞ ( t ) = V ar ( R ( t )) = E [ V ar ( M ( t )) | Λ( · ))] + V ar ( E [ M ( t ) | Λ( · )])= m ∞ ( t ) + V ar (cid:18)(cid:90) ∞ Λ( t − s ) P ( S > s ) d s (cid:19) . (7)As an aside, note that indeed m ∞ ( t ) (cid:54) v ∞ ( t ) , which reflects the overdispersion that we introduced.These expressions can easily be simplified using the observation that Λ( t ) is piecewise constant. Forthe case that S is exponentially distributed, µ ∞ ( t ) and v ∞ ( t ) can be evaluated in closed form in t = n ∆ with n ∈ N ; see Appendix A .2.3. Case study: MOL staffing for overdispersed hospital arrival data.
We continue byillustrating our approach and its performance in a case study. The data set used was provided bythe SEElab and originates from the emergency department (ED) of an Israeli hospital. It contains -minute resolution arrival counts of a -year time period ( − ), which covers a total of days. The average arrival volume per day, exceeding arrivals, is sufficienty large and themean length of stay (LOS) is almost hours.We aggregate different weekdays separately, which implies that we have N = 224 observations foreach day of the week (see Table 1). Figure 1 presents the sample mean and variance of the numberof arrivals per slot for a 5, 10, 15, 30 and 60 minute resolution, based on these observations.When considering the hourly mean arrival rates (same timescale as LOS) for each hour of the week itis fairly variable: with a time average of , its minimum is . and its maximum . . In conclusion,the level of nonstationarity is high and the rate is rapidly changing with respect to the LOS. Notethat this means that both of the methods mentioned in the introduction, SSA and PSA, would notbe accurate.An observation from Figure 1 is that different weekdays indeed show different patterns in the arrivalstream and also the level of overdispersion and nonstationarity is visibly different although data has Start day
End day
Total
Total
Mean LOS (min.)
St. dev. LOS (min.)
Table 1.
Summary statistics of the hospital ED.already been averaged over 224 samples. Note for example that Sunday (in Figure 1 the 4th day)is an exceptionally busy day with high peaks in the mean and variance, in contrast to Friday andSaturday (i.e., Israeli weekend).Furthermore, observe that the (positive) difference between the mean and the variance increases withthe chosen resolution, but only when choosing a resolution larger than 30 minutes overdispersionbecomes more apparent in the plots. Comparing the values corresponding to the different subfiguresin Figure 1, it shows that the growth in the variance is still roughly linear. In fact, our model predictswhat we observe in Figure 1: the variance in the number of arrivals (just like the mean) grows roughlylinear in the length of the time slot, but due to the presence of (nonnegative) correlation betweenrates in consecutive time slots, an extra term should be added to the variance when aggregatingdata from smaller time slots. On top of this visualization of the sample means and variances, we canalso compute the empirical covariance matrix to quantify the correlations between the arrival countsin all different time slots, for each of the resolutions.When fitting the model to the arrival data we find that choosing the parameters α , I and V ar ( W ) differently for different weekdays significantly improves the fit. As this modeling decision also affectsthe staffing rule and the subsequent performance analysis, we decided to simplify reality and examineonly isolated Sundays in the rest of this case study. Note that consequently we pretend that Sundayssucceed one another, so that time slots around midnight are correlated in the model although inreality there is a week of (ignored) events in between. It is expected that the error resulting fromthis simplification is small as the arrival volume is small around midnight, and even smaller for smallvalues of I .In the rest of this subsection we will present our rule’s performance when staffing a multi-serversystem where the arrival stream is taken from the data set, restricted to Sundays. In Section 4 wepresent a systematic evaluation using stylized input.Given ∆ = 1 hour and I = 10 (where we intentionally pick a large value of I to be on the safeside), we find by the statistical inference procedure that will be described in Section 3 (cf. Table 2)that α = 0 . and V ar ( W ) = 0 . are the best fit for the data. As this will be the input for the Arrivals . . . . . . .
5 minute resolution time slot
10 minute resolution time slot0 100 200 300 400 500 600
15 minute resolution time slot 0 50 100 150 200 250 300
30 minute resolution time slot
60 minute resolution time slot
Figure 1.
Mean (black) and variance (red) of the number of arrivals in each timeslot for various resolutions. The week starts at Thursday, as the first day arrival datawas recorded was a Thursday. The dashed blue line in the 60 minute resolution plotshows the sum of the variances in the corresponding 30 minute time slots. The gapbetween the dashed blue line and the red line indicate the presence of nonnegativecorrelation between the number of arrivals in consecutive time slots.simulation, we expect that using these parameters for the staffing rule will give the best result. Tocheck this, we will also generate the delay probabilities if we plug in different parameter settingsin the staffing rule, i.e., I = 0 , , (with corresponding α and V ar ( W ) as given in Table 3). Asin Table 1, the mean length of stay is minutes, so the hourly service rate to be used in thesimulation is µ = 60 / ≈ . . theoretical mean ( t ) theoretical vars ( t ) staffing rule 0.01 - no tuningstaffing rule 0.01 - after tuningstaffing rule 0.05 - no tuningstaffing rule 0.05 - after tuningstaffing rule 0.1 - no tuningstaffing rule 0.1 - after tuningsimulated mean ( t ) simulated vars ( t ) Figure 2.
Mean, variance and staffing rules given different service-level agreements.Figure 2 shows the empirical mean and variance obtained from the simulation, which determinesthe probability distribution of the number of patients per time slot. Observe that the theoreticalvalues at the end of the time slot, as given in Eqns. (13) and (14), more or less coincide with theseempirical values; see Figure 2. The prescribed number of servers in time slot n depends on theservice level that was set and is chosen according to the staffing rule in Eqn. (3), with t = n ∆ and β as in Eqn. (4). We compare ε ∈ { . , . , . } ; see the solid gray lines in Figure 2.Next, the empirical probability of exceeding the staffing level in an infinite-server system is computedusing the simulation results. As the staffing rule is based on analysis of infinite-server systems, itcan be expected that this probability behaves well: asymptotic normality predicts that it should beclose to the required level ε in every time slot, which implicitly says that the service level should bemore or less stable. However, Figure 3 shows that the exceedance probability sometimes crosses therequired service level and does not follow a smooth straight line. - server proxy0.05effective delay probabilityinfinite - server proxy0.01effective delay probabilityinfinite - server proxy Figure 3.
Exceedance vs delay probability for different service levels. Each servicelevel is designated by a different color, where the dashed line describes the effectivedelay probability for the finite-server system under study and the dotted line describesthe exceedance probability for the infinite-server proxy.This can partly be explained by rounding errors (note that the ‘tooths’ in the lines are often causedby a difference of server), and moreover it is noted that asymptotic results in the end are just approximations. All in all the performance of this very straightforward and easy-to-use staffing ruleis satisfying. But importantly, the infinite-server results can of course only be used as a proxy. Theactual delay probabilities from the finite-server setting (where the staffing rule dictates the numberof patients that can be served at a time) will be significantly larger due to queueing caused by thewaiting patients. If this queueing bias would result in a uniform shift upwards, the staffing rulewould still prove perfectly useful, as we can easily tune the delay probability down by tweaking β . Unfortunately this is not the case; Figure 3 shows a heavy spike around noon, so the system islocally performing unacceptably poorly. Note that in the (finite-server) setting with abandonmentsthe performance would certainly be better, depending on the abandonment rate. Now, instead ofgoing immediately into service as in the infinite-server setting, all customers initiate an exponentialclock (with a rate that might even be comparable to the service rate) right upon arrival, for potentialabandonment of the system. The infinite-server proxy is way more accurate in such a setting. InSection 4 we will consider this adaptation, but for now we try to further improve the staffing rulefor the basic setting (i.e., the setting without abandonments).Note that, although for most of the day the delay probability seems rather stable, around noon ittakes on values twice the targeted service level. Comparing Figure 2, we find that around noon,which is not incidentally precisely the area where the increase in load is extremely high (due tononstationarity of the arrival stream), the prescribed number of servers follows the same slope asthat of the square-root of the variance. However, apparently this is not enough; the system cannot deal with the backlog that is rapidly building up around noon. Based on this observation, wecook up a heuristic that could potentially overcome this hurdle: the hedge β (cid:112) v ∞ ( t ) is replaced bya more involved one, that accounts for extreme fluctuations in the arrival rate in settings where thelevel of nonstationarity is high. Slope heuristic.
Let v n the ratio between the variance in time slot n + 1 and n . The idea is to scaleup the number of servers when v n (cid:29) while mildly reducing the number of servers when v n < ,without changing the total number of staffed servers over the day. It is important to only makesubtle changes, so that the ‘shape’ (viz., Figure 2) prescribed by the infinite-server proxy staysunaltered. Consequently, we are after an increasing function f ( x ) with the property that | f ( x ) − | ≤ | x − | . Functions f δ ( x ) = x δ with < δ ≤ satisfy these conditions and have the advantage that they caneasily be tuned via the parameter δ . We arrive (for t = n ∆ ) at s mn ∆ = m ∞ ( n ∆) + β ( v n ) /δ (cid:112) v ∞ ( n ∆) , (8)for n = 1 , . . . , . Then δ can be picked such that the variance of the resulting delay probability(given a staffing level according to s δn ∆ for n = 1 , . . . , ) is minimized. Alternatively, a few valuesfor δ are compared to arrive at a value for which this variance is relatively small. ♦ Remark 2.1.
Note that Eqn. (8) simplifies to Eqn. (3) if λ ( t ) is constant. The heuristic introducesan extension to the staffing rule that was originally proposed to account for nonstationarity; if thereis no nonstationarity present ( λ ( t ) ≡ λ ), the slope-adapted rule reduces to the original rule, in whichcase the latter’s performance is satisfactory.Moreover, note that in the infinite-server setting performance would not improve by using the rulein Eqn. (8); in this setting Eqn. (3) is the best we can get. That is to say, implementation of thisheuristic is useful in situations where an extremely steep slope in the arrival rate causes an avalancheof queueing patients once the number of servers is restricted. ♦ Figure 4.
Using the slope-adapted staffing rule improves the stability of the delayprobability. The delay probabilities with no tuning (dashed lines) coincide withthose in Fig. 3. The dotted lines depict the delay probabilities after tuning and areclearly more stable than before. With δ ε denoting the selected δ for service level ε , δ . = 3 / , δ . = 1 / and δ . = 1 / .Figure 4 compares the performance of the slope-adapted staffing rule (cf. Eqn. (8)) with that of theoriginally proposed staffing rule (cf. Eqn. (3)). We observe better stability with (approximately)the same number of servers (in the example of Figure 4 the total number of servers for both rulesdiffers by 2 or 3 servers) and (on average) a slightly smaller delay probability. Nevertheless, thedelay probabilities still exceed the targeted service level ε . Hence, the simple adaptation of choosinga higher value of β in Eqn. (3) would further improve performance.3. Statistical procedures
In this section we describe how to determine, based on historical data, the parameters in the arrivalstream model introduced in Section 2, i.e., λ j for j = 0 , . . . , N − , α , I , and V ar ( W ) . Given a value ∆ > it is enough to know, for each j : ¯Λ j , the average number of arrivals in ∆ j := [ j ∆ , ( j + 1)∆) ,and Σ := Σ( α, V ar ( W )) , the covariance matrix, representing for j = 0 , . . . , N − and k = 1 , . . . , I the (nonnegative) covariance between the number of arrivals in ∆ j and ∆ j + k (note that here and inwhat follows, the indices in the subscripts should be taken modulo N ; for the sake of readability wedo not write that explicitly). Deterministic trend.
The average number of arrivals ¯Λ j should correspond to the average over valuesof the mixed Poisson random variable with random parameter Λ j := λ j · c α I (cid:88) (cid:96) =0 α (cid:96) W j − (cid:96) , with c α := (1 − α ) / (1 − α I +1 ) . Using that E Λ j = λ j , the W j have unit mean and that c α is anormalizing constant, the ¯Λ j are unbiased estimators for the λ j . Covariance matrix.
Recall that the covariance of two independent mixed Poisson random variables(meaning that the enveloping Poisson random variables are independent) with dependent parametersequals the covariance of the parameters. In addition, recall that the variance of a mixed Poissonrandom variable is the sum of the expectation and variance of its parameter.Given that I ≤ (cid:98) ( N − / (cid:99) (cf. Appendix C), we obtain the following expressions for the entries ofthe covariance matrix: Σ j,j = E Λ j + V ar (Λ j ) = λ j + λ j c α − α I +1) − α V ar ( W ) ; (9) Σ j,j + k = Σ j + k,j = C ov(Λ j , Λ j + k )= λ j λ j + k c α C ov (cid:32) I (cid:88) (cid:96) =0 α (cid:96) W j − (cid:96) , I (cid:88) (cid:96) =0 α (cid:96) W j + k − (cid:96) (cid:33) (10) = λ j λ j + k c α α k − α I − k +1) − α V ar ( W ) , (11)where it’s noted that Σ j,j + k = Σ j + k,j = 0 for k > I . Let C k ( α, I ) := c α α k − α I − k +1) − α . Then Eqns.(9) and (11) can be captured by Σ j,j + k = Σ j + k,j = λ j (cid:0) { k =0 } + λ j + k C k ( α, I ) V ar ( W ) (cid:1) , (12)for k = 0 , , . . . , I (and otherwise). Note that by l’Hôpital’s rule lim α ↑ C k ( α, I ) = lim α ↑ c α α k − α I − k +1) − α = I − k + 1( I + 1) . Hence we set C k (1 , I ) := ( I − k + 1) / ( I + 1) . Procedure for α , I and V ar ( W ) . The idea is to vary I in an outer loop and to estimate α and V ar ( W ) (for any given I ); one could then compare how much gain is made by using different I with respect to the base case where I = V ar ( W ) = 0 (standard Poisson) and I = 0 (no correlation).Subsequently, it makes sense to select the largest I that is a significant improvement over I − (orover the standard Poisson case, where I = V ar ( W ) = 0 ). Note that the model only allows for valuesof I ranging from to . I α V ar ( W ) MSE* MSE Gain (%)Poisson - - - 7.435 0.000 %0 - 0.017 5.562 5.974 19.645 %1 1.000 0.041 5.207 4.343 41.583 %2 1.000 0.060 4.267 3.476 53.244 %3 1.000 0.075 3.389 2.998 59.681 %4 1.000 0.089 2.821 2.685 63.881 %5 1.000 0.102 2.569 2.463 66.872 %6 0.907 0.112 2.344 2.342 68.501 %7 0.879 0.121 2.182 2.231 69.992 %8 0.866 0.129 2.246 2.120 71.481 %9 0.867 0.138 2.234 2.015 72.905 %10 0.866 0.146 2.098 1.936 73.963 %11 0.861 0.152 1.949 1.894 74.520 %
Table 2.
Fitted parameters for Wednesday, Hospital 3.
I α V ar ( W ) MSE* MSE Gain (%)Poisson - - - 12.253 0.000 %0 - 0.015 9.818 9.752 20.4 %1 1.00 0.034 11.212 7.485 38.9 %5 1.00 0.084 6.238 4.566 62.7 %10 0.81 0.11 4.414 4.158 66.1 %
Table 3.
Fitted parameters for Sunday, Hospital 3.To be able to determine the values of α and V ar ( W ) given λ j , λ j + k and I , we need the empiricalcovariance matrix Σ derived from the arrival data. Note that an estimate of any two nonzero entriesof Σ provides enough information to solve for α and V ar ( W ) , after having equated them to theexpression in Eqn. (12). However, each nonzero pair leads to a different solution. We wish todetermine values for α and V ar ( W ) such that the theoretical covariance matrix Σ( α, V ar ( W )) asgiven by Eqn. (12) is the ‘best’ approximation for Σ . Therefore, the next step in the procedure is tominimize the average of the entrywise mean squared errors, where we sum over the entries for whichthe theoretical covariance matrix is nonzero (noting that the number of nonzero entries, being equalto N (2 I + 1) , depends on the choice of I ). In Table 2 this value is labeled with MSE*, with a separatecolumn for the exact MSE values where all entries of the empirical covariance matrix are takeninto account. The gain is computed as the relative gain in (exact) MSE compared to the standard Poisson case (where I = V ar ( W ) = 0 ). We observe that from I = 5 , not much improvement is stillto be gained, so I = 5 seems to be a good choice when we aim for moderate complexity and a goodfit. Some more examples.
In Table 3 the same procedure is used to obtain the best fit for I = 0 , , , ,which are used in the case study in Section 2.3. Similar gain percentages in MSE are obtained bychoosing I larger, although in Table 2 I = 5 and I = 10 achieve a better fit than in Table 3. Atthe same time, the variance of W as well as the correlation parameter α is consistently smaller inthe data set that corresponds to Sundays; although total arrival volume is larger here, temporalcorrelation and overdispersion seems to be less prominent. I α V ar ( W ) MSEPoisson - - 0.000 1.00 0.172 87.81 0.871 0.343 18.52 0.561 0.425 3.993 0.518 0.467 0.7754 0.505 0.489 0.1045 0.500 0.500 06 0.499 0.507 0.02407 0.498 0.510 0.04158 0.496 0.510 0.0483
Table 4.
Fitted parameters given I and corresponding MSE.We add a theoretical example to assess the precision of this minimization method. With the λ j as above, set α ∗ = V ar ( W ) ∗ = 0 . and I ∗ = 5 . As these parameters together define the arrivalprocess, this gives a certain covariance matrix Cov ∗ . We use the minimization method to find,given some choice of I ∈ { , , . . . , } , the optimal values for α and V ar ( W ) in terms of the MSEof Cov( α, V ar ( W )) with respect to Cov ∗ . The results can be found in Table 4, together with thecorresponding MSE. Observe that the method recovers the true values α ∗ and V ar ( W ) ∗ in casewe set I = I ∗ = 5 . For lower degree of correlation, it is found that a larger value for α (moredependence) is compensated by a smaller value for V ar ( W ) (less overdispersion), however apartfrom the case I = 4 choosing I too small inevitably leads to a big loss in precision. For I = 4 theMSE is acceptably small. On the other hand, setting I too large leads only to small errors, whichmeans that selecting a value I above the true value leads to marginal differences. Performance
In this section we extend the numerical work on the performance of the presented staffing rules inthe finite-server setting with described arrival stream, based on simulations instead of data. Thegoal is to assess the individual effects of nonstationarity, temporal correlation and overdispersion.Moreover, to reduce the gap between the finite-server system and its infinite-server proxy we addan abandonment rate θ to the model. Note that in any service system that involves waiting, it isnatural to have a (possibly small) positive abandonment rate.We start with a particular stylized instance for the arrival stream, again inspired by the hospitaldata, i.e., the levels of overdispersion, nonstationarity and temporal correlation are comparable. Thedaily pattern is represented by a sine function with a cycle length of hours (with ∆ = 1 hour),having a dip early in the morning (at 4:30) and a peak late in the afternoon (at 16:30). That is, λ j = N + pN · sin( 2 π
24 ( j + 13 . for j = 0 , . . . , , where N is the system size and p the level of nonstationarity. The parameters are set as in Table 5. N = . I = p = . α = µ = . V ar ( W ) = . Table 5.
Parameter setting base case.Note that, as in Section 2.3, the concerning system is fairly small whereas the level of nonstationarityis extremely high. As the service rate is low, this means that the patient ‘sees’ effectively differentarrival rates during its stay and nonstationarity can not be ignored. The correlation structure isabundantly present and the level of overdispersion seems mild (though nonzero). However, theeffect of V ar ( W ) being positive on the size of the hedge (i.e., on √ v ∞ ) is quite large; comparecolumns 1 and 3 in Table . On the other hand, taking I > slightly mitigates this effect; comparecolumns 2 and 3 in Table . This table was included to show the effects of different modeling choicesmade independent of the impact of a (wildly fluctuating) daily pattern, which will be added in theexperiment that follows.In Table 6 the staffing levels for × different settings are given. Here abandonments are incorporatedto different extents: the abandonment rate is θ = a ∗ µ , for a = 0 , . , (’no’, ’mild’ or ’max’). Notethat incorporating abandonments does not affect the prescription for the staffing level, as our staffingrule does not account for it.We find that in this instance with stationary deterministic daily pattern, performance does notchange significantly when a correlation structure is added to the model, where it is noted that weaccount for it in the staffing rule (in this case that means that less servers were used to achieveapproximately the same delay probability). We do however need significantly more servers to attain a comparable level for the delay probability when switching from the ‘no overdispersion’ setting(column 1) to the setting with overdispersion (column 3), which of course is the motivation forthe staffing rule introduced in this paper. Note that performance gets worse anyway, despite thecomplication in the number of servers.From Table 6 it becomes very clear that firm improvement in performance can be achieved byincorporating a positive abandonment rate. Nevertheless, we see that even without a daily pattern,the β constant needs to be tuned somewhat until the delay probabilities match the targetedprobabilities, partly due to the overdispersion (the first column displays better performance) andpartly due to the inaccuracy of the infinite-server proxy (when abandonments are incorporatedperformance gets better until it’s nearly perfect). Strangely, only in the cases where ε = 0 . withno/mild abandonments, the performance got worse when correlation was left out. Of course theproxy is least accurate in this case, but apparently having dependence between arrival rates ‘helps’here. standard I > I = 0 √ v ∞ s . s . s . no abandonments 0.098 0.12 0.13 ε = 0 . mild abandonments 44 0.086 45 0.10 46 0.11max abandonments 0.079 0.093 0.10no abandonments 0.051 0.055 0.068 ε = 0 . mild abandonments 46 0.046 48 0.048 49 0.060max abandonments 0.043 0.044 0.056no abandonments 0.011 0.017 0.016 ε = 0 . mild abandonments 50 0.010 52 0.015 55 0.015max abandonments 0.0099 0.014 0.014 Table 6.
Delay probability obtained through simulation, for the setting withoutnonstationarity (we set p = 0 ). The first column gives the probability when setting V ar ( W ) = 0 , in the third column the correlation structure is ignored ( I = 0 ). Themiddle column is the delay probability as dictated by our model given the statedparameter setting (see Table 5).From the plots in Figure 5 we observe that incorporating only mild abandonments significantlyimproves the performance of our staffing rule, which makes sense as the infinite-server proxy is moreaccurate for finite-server models with abandonments. Note that the finite-server setting endowedwith an abandonment rate θ = µ coincides with the infinite-server setting, the setting in the lastrow of plots. Note that the somewhat erratic nature of the delay probability is due to inevitablerounding errors resulting from the fact that the number of servers needs to be integer. (a) (b) (c)(d) (e) (f)(g) (h) (i) Figure 5.
Delay probability in finite-server setting, staffing level according new rule.In the first column we set V ar ( W ) to zero (Subfigures (a), (d) and (g)), in the arrival stream as wellas in the staffing rule.Given that the delay probabilities in this setting define some sort of baseline for the performance(this should be the easiest setting to handle), it is remarkable that the performance does not get(significantly) worse when taking into account overdispersion and correlation. In that sense it lookslike our staffing rule is prescribing the correct number of servers. The plots even suggest slightimprovement in many settings. Comparing the overdispersed setting where correlation is left out(column 3) with the setting with both overdispersion and correlation (column 2), there is only a veryslight improvement in performance over all plots with different abandonment rates and staffing levels.That is, our rule seems to account well for overdispersion, but it is struggling slightly harder to dealwith the correlation structure. However, it can be concluded that nonstationarity is the main factorthat complicates achieving stable delay probabilities, mostly in the setting without abandonments.In order to make a fair comparison with the case study in Section 2.3, it is necessary to apply theslope-adapted staffing rule (cf. Eqn. (8)) here as well. We will only apply it to the case with noabandonments, hence we mirror Figure 5 (a), (b) and (c): see Figure 6 (a), (b) and (c). From Figure 6it can be concluded that the slope-adapted staffing rule indeed stabilizes the delay probabilities (a) (b) (c) Figure 6.
Delay probability in finite-server setting, staffing level determined withslope-adapted staffing rule. Here δ from Eqn. (8) is of the form δ = k , for the valueof k that maximally stabilizes the delay probability.over the day. However, further improvement could be made by tweaking β , to get the stabilizedprobabilities below the targeted level. The resulting improvement is not shown, as the procedureand its effect are trivial. 5. Conclusion and discussion
In this paper we propose new staffing rules for a specific queueing model with overdispersed andnonstationary input with temporal correlation. The objective is to stabilize the delay probabilitythroughout the day around a fixed target value, which the final staffing rule developed succeeds todo.In the numerical experiments in Section 2.3 the originally proposed rule based on an infinite-serverproxy proves insufficient for staffing purposes. The main complication turns out to be nonstationarity.Considering the same model with abandonments, we observe significantly better performance, due tothe fact that the infinite-server proxy is more accurate for finite-server models with abandonments.Applying the introduced slope heuristic, the adapted staffing rule renders a major improvement(already without abandonments!).The observed performance is robust for the choice of parameters for overdispersion and temporalcorrelation; as long as the combination of parameters results in an accurate estimate for the variancein the number of arriving customers, the prescribed (slope-adapted) staffing level is appropriate. InAppendix B, it is shown that this variance is decreasing in α and I and at the same time increasingin V ar ( W ) , so that different parameter settings can result in the same variance. Although thestatistical procedure in Section 3 does not lead to a unique ‘optimal’ choice for the parameters α , I and V ar ( W ) , because of this robustness it is sufficient to select a reasonable parameter setting.Note that the implementation of nonstationarity in the model is rather straightforward: fitting aconstant arrival rate to fixed time slots is the simplest and also a widely used procedure to implementtime-of-day or time-of-week effects. It is remarked, though, that the discontinuities might causepoor predictions close to the slot boundaries. Therefore, in [43] a slight adaptation is suggested: itis proposed to use piecewise linear (hence continuous) rates. Acknowledgments.
The authors would like to thank Johan van Leeuwaarden (Tilburg University)for useful suggestions and advice.
References [1] Z. Aksin, M. Armony, and V. Mehtotra. The modern call center: A multi-disciplinary perspectiveon operations management research.
Production and Operations Management , 16(6):665–688,2007.[2] A.N. Avramidis, A. Deslauriers, and P. L’Ecuyer. Modeling daily arrivals to a telephone callcenter.
Management Science , 50(7):896–908, 2004.[3] A. Bassamboo, R.S. Randhawa, and A. Zeevi. Capacity sizing under parameter uncertainty:safety staffing principles revisited.
Management Science , 56(10):1668–1686, 2010.[4] A. Bassamboo and A. Zeevi. On a data-driven method for staffing large call centers.
OperationsResearch , 57(3):714–726, 2009.[5] S.C. Borst, A. Mandelbaum, and M.I. Reiman. Dimensioning large call centers.
OperationsResearch , 52(1):17–34, 2004.[6] B.P.K. Chen and S.G. Henderson. Two issues in setting call center staffing levels.
Annals ofOperations Research , 108(1):175–192, 2001.[7] N. Gans, G. Koole, and A. Mandelbaum. Telephone call centers: Tutorial, review, and researchprospects.
Manufacturing & Service Operations Management , 5(2):79–141, 2003.[8] W.K. Grassmann. Finding the right number of servers in real-world queuing systems.
Interfaces ,18(2):94–104, 1988.[9] L. V. Green, P. J. Kolesar, and W. Whitt. Coping with time-varying demand when settingstaffing requirements for a service system.
Production and Operations Management , 16(1):13–39,2007.[10] L.V. Green and P. Kolesar. The pointwise stationary approximation for queues with nonstation-ary arrivals.
Management Science , 37(1):84–97, 1991.[11] L.V. Green, P. Kolesar, and A. Svoronos. Some effects of nonstationarity on multiserverMarkovian queueing systems.
Operations Research , 39(3):502–511, 1991.[12] I. Gurvich, J. Luedtke, and T. Tezcan. Staffing call-centers with uncertain demand forecasts: achance-constrained optimization approach.
Management science , 56(7):1093–1115, 2010.[13] S. Halfin and W. Whitt. Heavy-traffic limits for queues with many exponential servers.
OperationsResearch , 29(3):567–588, 1981.[14] B. He, Y. Liu, and W. Whitt. Staffing a service system with non-Poisson nonstationary arrivals.
Probability in the Engineering and Informational Sciences , 30(4):593–621, 2016.[15] M. Heemskerk, J.S.H. van Leeuwaarden, and M. Mandjes. Scaling limits for infinite-serversystems in a random environment.
Stochastic Systems , 7(1):1–31, 2017. [16] R. Ibrahim, P. L’Ecuyer, N. Regnard, and H. Shen. On the modeling and forecasting call centerarrivals. In
Proceedings of the 2012 Winter Simulation Conference , pages 1–12, 2012.[17] R. Ibrahim, H. Ye, P. L’Ecuyer, and H. Shen. Modeling and forecasting call center arrivals: aliterature survey and a case study.
International Journal of Forecasting , 32(3):865–874, 2016.[18] A.J.E.M Janssen, J.S.H. van Leeuwaarden, and A.P. Zwart. Refining square-root safety staffingby expanding Erlang-C.
Operations Research , 59(6):1512–1522, 2011.[19] O.B. Jennings, A. Mandelbaum, W.A. Massey, and W. Whitt. Server staffing to meet time-varying demand.
Management Science , 42(10):1383–1394, 1996.[20] G. Jongbloed and G. Koole. Managing uncertainty in call centers using Poisson mixtures.
Applied Stochastic Models in Business and Industry , 17(4):307–318, 2001.[21] O. Jouini and S. Benjaafar. Appointment scheduling with non-punctual arrivals.
IFACProceedings Volumes , 42(4):235–239, 2009.[22] S.-H. Kim, V. Vel, W. Whitt, and W.C. Cha. Poisson and non-Poisson properties in appointment-generated arrival processes: The case of an endocrinology clinic.
Operations Research Letters ,43(3):247–235, 2015.[23] S.-H. Kim and W. Whitt. Are call center and hospital arrivals well modeled by nonhomogeneousPoisson processes?
Manufacturing & Service Operations Management , 16(3):464–480, 2014.[24] S.-H. Kim and W. Whitt. Choosing arrival process models for service systems: Tests of anonhomogeneous Poisson process.
Naval Research Logistics (NRL) , 61(1):66–90, 2014.[25] S.-H. Kim, W. Whitt, and W.C. Cha. A data-driven model of an appointment-generated arrivalprocess at an outpatient clinic.
INFORMS Journal on Computing , 30(1):181–199, 2017.[26] Y.L. Koçaga, M. Armony, and A.R. Ward. Staffing call centers with uncertain arrival rate andco-sourcing.
Production and Operations Management , 24(7):1101–1117, 2015.[27] S. Liao, G. Koole, C. van Delft, and O. Jouini. Staffing a call center with uncertain non-stationaryarrival rate and flexibility.
OR Spectrum , 34(3):691–721, 2012.[28] S. Maman. Uncertainty in the demand of service: The case of call centers and emergencydepartments. Master thesis, Technion - Israel Institute of Technology, Haifa, 2009.[29] B.W.J. Mathijsen, A.J.E.M. Janssen, J.S.H. van Leeuwaarden, and A.P. Zwart. Robust heavy-traffic approximations for service systems facing overdispersed demand.
Queueing Systems ,90(3-4):257–289, 2018.[30] V. Mehrotra, O. Ozlük, and R. Saltzmann. Intelligent procedures for intra-day updating of callcenter agent schedules.
Production and operations management , 19(3):353–367, 2010.[31] T.R. Robbins, D.J. Medeiros, and T.P. Harrison. Does the Erlang C model fit in real callcenters? In
Proceedings of the Winter Simulation Conference , pages 2853–2864. WinterSimulation Conference, 2010.[32] M.H. Rothkopf and S.S. Oren. A closure approximation for the nonstationary
M/M/s queue.
Management Science , 25(6):522–534, 1979. [33] S.G. Steckley, S.G. Henderson, and V. Mehrotra. Forecast errors in service systems. Probabilityin the Engineering and Informational Sciences , 23(2):305–332, 2009.[34] J. Tan, H. Feng, X. Meng, and L. Zhang. Heavy-traffic analysis of cloud provisioning. In
Proceedings of the 24th International Teletraffic Congress , pages 1–8, 2012.[35] J.S.H. van Leeuwaarden, B. Mathijsen, and B. Zwart. Economies-of-scale in many-serverqueueing systems: tutorial and partial review of the QED Halfin-Whitt heavy-traffic regime.
SIAM Review , 61(3):403–440, 2019.[36] J.S.H. van Leeuwaarden, B.W.J. Mathijsen, and F. Sloothaak. Cloud provisioning in the QEDregime. In
Proceedings of the 9th EAI International Conference on Performance EvaluationMethodologies and Tools , pages 180–187, 2016.[37] W. Whitt. The pointwise stationary approximation for M t /M t /s queues is asymptoticallycorrect as the rates increase. Management Science , 37(3):307–314, 1991.[38] W. Whitt. Understanding the efficiency of multi-server service systems.
Management Science ,38(5):708–723, 1992.[39] W. Whitt. Approximations for the
GI/G/m queue.
Production and Operations Management ,2(2):114–161, 1993.[40] W. Whitt. Dynamic staffing in a telephone call center aiming to immediately answer all calls.
Operations Research Letters , 24(5):205–212, 1999.[41] J. Zan.
Staffing service centers under arrival-rate uncertainty . PhD thesis, University of Texas,2012.[42] B. Zhang, J.S.H. van Leeuwaarden, and B. Zwart. Staffing call centers with impatient customers:refinements to many-server asymptotics.
Operations Research , 60(2):461–474, 2012.[43] Z. Zheng and P.W. Glynn. Fitting continuous piecewise linear Poisson intensities via maximumlikelihood and least squares. In
Proceedings of the 2017 Winter Simulation Conference , pages1740–1749, 2017.
Appendix A. Computations for infinite-server queue
In this appendix we calculate m ∞ ( t ) and v ∞ ( t ) in terms of λ ( t ) , α and V ar ( W ) , which can beextracted from arrival data by means as proposed in Section 3. Let ¯ F ( s ) := P ( S > s ) . In thisappendix we consider exponentially distributed service times, but similar calculations can be donefor other distributions in a straightforward manner.Let t = n ∆ for some n ∈ Z ≥ . We assume that λ ( s ) is a periodic step function with step size ∆ and cycle length N (i.e., λ (0) = λ ( N ∆) ) and write λ k := λ ( t ) for t ∈ [ k ∆ , ( k + 1)∆) for somenon-negative value λ k . As a consequence, for λ , . . . , λ N − , we have λ k = λ (cid:96) if k mod N = (cid:96) mod N. Let us start with evaluating m ∞ ( t ) for this setting of periodic λ ( · ) and exponential service times(with mean µ − ). In the first place, an elementary calculation reveals that Eqn. (6) simplifies to (with t = n ∆ ), m ∞ ( t ) = E (cid:20)(cid:90) ∞ Λ( t − u ) ¯ F ( u ) d u (cid:21) = 1 − e − µ ∆ µ ∞ (cid:88) j =1 λ n − j (e − µ ∆ ) j − For j = 1 , . . . , N , we introduce (using the periodicity) κ j ( n ) := ∞ (cid:88) (cid:96) =1 λ n − ( (cid:96) − N − j (e − µ ∆ ) ( (cid:96) − N + j − = λ n − j (e − µ ∆ ) j − ∞ (cid:88) (cid:96) =1 (e − µ ∆ N ) (cid:96) − = λ n − j · e − µ ∆( j − − e − µ ∆ N . This leads to an expression for m ∞ ( t ) in terms of a finite sum: m ∞ ( t ) = 1 − e − µ ∆ µ N (cid:88) j =1 κ j ( n ) = 1 − e − µ ∆ − e − µ ∆ N µ N (cid:88) j =1 λ n − j e − µ ∆( j − . (13)We now move on to compute v ∞ ( t ) . To this end, we define γ ( j ) := λ j c α (cid:90) ( n − j )∆( n − j − ¯ F ( u ) d u. The idea is to rearrange the contributions to the random arrival rate due to each of the W j in theexpression for v ∞ ( t ) in Eqn. (7): V ar (cid:18)(cid:90) ∞ Λ( t − u ) ¯ F ( u ) d u (cid:19) = V ar ∞ (cid:88) j =1 λ n − j (cid:16) c α I (cid:88) (cid:96) =0 α (cid:96) W n − j − (cid:96) (cid:17) (cid:90) j ∆( j − ¯ F ( u ) d u = V ar ∞ (cid:88) j =1 (cid:16) λ n − j c α (cid:90) j ∆( j − ¯ F ( u ) d u (cid:17) I (cid:88) (cid:96) =0 α (cid:96) W n − j − (cid:96) = V ar ∞ (cid:88) j =1 γ ( n − j ) I (cid:88) (cid:96) =0 α (cid:96) W n − j − (cid:96) = V ar ∞ (cid:88) j =1 (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 α (cid:96) γ ( n − j + (cid:96) ) (cid:19) W n − j Noting that the W j are independent and identically distributed, the expression in the previousdisplay becomes V ar ( W ) I (cid:88) j =1 (cid:18) j − (cid:88) (cid:96) =0 α (cid:96) γ ( n − j + (cid:96) ) (cid:19) + V ar ( W ) ∞ (cid:88) j = I +1 (cid:18) I (cid:88) (cid:96) =0 α (cid:96) γ ( n − j + (cid:96) ) (cid:19) = V ar ( W ) c α (e µ ∆ − µ ∞ (cid:88) j =1 (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 α (cid:96) λ n − j + (cid:96) e − µ ( j − (cid:96) )∆ (cid:19) , where we use that (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 α (cid:96) γ ( n − j + (cid:96) ) (cid:19) = I ∧ ( j − (cid:88) (cid:96) =0 α (cid:96) λ n − j + (cid:96) c α (cid:90) ( j − (cid:96) )∆( j − (cid:96) − ¯ F ( u ) d u = c α · (e µ ∆ − µ B j with B j := (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 α (cid:96) λ n − j + (cid:96) (e − µ ∆ ) j − (cid:96) (cid:19) . The next step is again to exploit the periodicity. For this we study (cid:80) ∞ j =1 B j , under the assumption I < N (which is fairly natural). Elementary calculus reveals that v ∞ ( t ) can be expressed as a finitesum, due to ∞ (cid:88) j =1 B j = ∞ (cid:88) j =1 (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 α (cid:96) λ n − ( j − (cid:96) ) (e − µ ∆ ) j − (cid:96) (cid:19) = ∞ (cid:88) j =1 α j (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) = N (cid:88) j =1 ∞ (cid:88) k =1 α j (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) − ( k − N (cid:18) e − µ ∆ α (cid:19) j − (cid:96) +( k − N (cid:19) = N (cid:88) j =1 ( α j (cid:18) I ∧ ( j − (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) + ∞ (cid:88) k =1 α j (cid:18) I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) − kN (cid:18) e − µ ∆ α (cid:19) j − (cid:96) + kN (cid:19) )= I (cid:88) j =1 α j (cid:18) ( j − (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) + N (cid:88) j = I +1 α j (cid:18) I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) + N (cid:88) j =1 α j ∞ (cid:88) k =1 (cid:18) e − µ ∆ α (cid:19) kN (cid:18) I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) (cid:19) = I (cid:88) j =1 α j (cid:18) ( j − (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) + N (cid:88) j = I +1 α j (cid:18) I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) + 1 α N e µ ∆ N − N (cid:88) j =1 α j (cid:18) I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) (cid:19) . Note that for convergence of the infinite series, we have to assume that e − µ ∆ < α .The final expression for v ∞ ( t ) is as follows: v ∞ ( t ) = V ar ( W ) c α (e µ ∆ − µ · N (cid:88) j =1 α j · D j , (14)where D j := (cid:18) ( j − ∧ I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) + 1 α N e µ ∆ N − (cid:18) I (cid:88) (cid:96) =0 λ n − ( j − (cid:96) ) (cid:18) e − µ ∆ α (cid:19) j − (cid:96) (cid:19) . Appendix B. Variance of the arrival process as a function of the parameter space
In this appendix we show that the variance of the arrival process is decreasing in α and I andincreasing in V ar ( W ) . With the nonstationary doubly-stochastic arrival rate process Λ( t ) given byEqn. (2), we get a nonstationary mixed Poisson arrival process. This process is overdispersed, whichbecomes visible once we write down its variance: V ar (cid:0)(cid:0) Poisson(Λ( t )) (cid:1)(cid:1) = λ j + λ j (1 − α ) (1 − α I +1 ) − α I +1) − α V ar ( W ) , (15)which is larger than λ j as the second term at the right-hand side of Eqn. (15) is positive for V ar ( W ) > . It is hence directly seen that Eqn. (15) is increasing in V ar ( W ) . Observe that thefactor (1 − α ) (1 − α I +1 ) − α I +1) − α = 1 − α α α I +1 − α I +1 (16) depends both on α and I . We state and prove that Eqn. (16) is (strictly) decreasing in α and I .Note that indeed α I +1 − α I +1 < α I − α I for I = 0 , , . . . , for all α ∈ (0 , . Now consider the function f I ( α ) = − α α α I − α I for some I > (for I = 1 we find thisfunction is constant; note that this corresponds to the case where I = 0 in our model, i.e. there isno correlation between past time slots). After taking the logarithm and differentiating, we end upwith the condition that the function is (strictly) decreasing if Iα I − ( 11 + α I + 11 − α I ) <
11 + α + 11 − α . Rewriting gives Iα I − < I − (cid:88) k =0 α k = α I − + (cid:80) ( I − / − k =0 (cid:0) α k + α I − − k ) (cid:1) if I is odd (cid:80) I/ − k =0 (cid:0) α k + α I − − k ) (cid:1) if I is even , and these two cases are easy to check individually, since α I − < (cid:0) α k + α I − − k ) (cid:1) for all relevant k . Appendix C. Constraint on I In this appendix we explain why we take I , the number of elapsed time slots that affect the busynessfactor, to be at most equal to (cid:98) ( N − / (cid:99) . Note that the number of nonzero entries in the covariancematrix equals N (2 I + 1) , as for each time slot j = 0 , . . . , N − we have a nonzero entry on thediagonal and for k = 1 , . . . , I both Σ j,j + k and Σ j,j − k are nonzero. Of course the dimension of thematrix only allows for N entries. In other words: N (2 I + 1) ≤ N must hold, i.e., I ≤ (cid:98) ( N − / (cid:99) . (17)To be even more precise, strictly it is only required to set I ≤ (cid:98) N/ (cid:99) , however in the case where N is even, for k = I = N/ we should replace Eqn. (11) by Σ j,j + N/ = Σ j + N/ ,j = C ov(Λ j , Λ j + N/ )= λ j λ j + N/ c α C ov (cid:32) I (cid:88) (cid:96) =0 α (cid:96) W j − (cid:96) , I (cid:88) (cid:96) =0 α (cid:96) W j + k − (cid:96) (cid:33) = λ j λ j + N/ c α α N/ V ar ( W ) , (18)for j = 0 , . . . , N − .The following example serves as an illustration of the complication that arises when I is not restrictedas in Eqn. (17). Example 1.
Let N = 24 . Then Eqn. (11) holds for k = 1 , . . . , (when we pick I = 11 ). If wewere to choose I = 12 and used Eqn. (10) for the covariance between arrivals in ∆ and ∆ , wewould get λ λ c α C ov (cid:0) α I W + · · · + αW + W , α W + · · · + αW + W (cid:1) . (19)However, the second occurence of W in Eqn. (19) is incorrect and should be written as W (cid:48) : it isdescribes an i.i.d. copy of W . As a result, the covariance just equals λ λ c α α V ar ( W ) (cf. Eqn.(18)). Note that it’s in fact still possible to write all nonnegative covariances in a × -matrix.Namely, for I = 11 still N entries equal zero; as for I = 12 and k = N/ the entries Σ j,j + k and Σ j,j − k happen to coincide (for j = 0 , . . . , N − ), these N values exactly fill up the ‘previouslyunoccupied’ entries.If however, we had chosen I = 13 , we would need to write both the covariance between arrivals in ∆ and ∆ (where time slot ∆ passes first) and the covariance between arrivals in ∆ and the ∆ after (i.e., ∆ passes first) on the same entry, but their values do not match. ♦ All in all, we see that it makes sense to exclude (for simplicity)
I > (cid:98) ( N − / (cid:99)(cid:99)