QBD approximations of a call center queueing model with general patience distribution

Computers & Operations Research 35 (2008) 2463 – 2481 www.elsevier.com/locate/cor

QBD approximations of a call center queueing model with general patience distribution Ken’ichi Kawanishi∗ Department of Computer Science, Gunma University, 1-5-1 Tenjin-cho, Kiryu-shi, Gunma 376-8515, Japan Available online 22 January 2007

Abstract This paper investigates a computationally practical way for analyzing a call center queueing model, i.e., a finite-capacity, multiserver queueing model, where each server goes on a single vacation. Poisson arrival process and exponential service and vacation times are assumed. We also assume that each customer may leave the queue due to impatience. Customers’ patience times are i.i.d. random variables with a general distribution. Level-dependent finite QBD (quasi-birth–death) processes are employed to approximate such a queueing model. Two approaches are considered. The first one uses the phase-type (PH) distribution to approximate the general patience distribution, whereas the second one is based on the idea of replacing the eventual reneging of customers with balking. We find that the first approach is almost impossible to compute numerically due to the exponential growth of the size of the block matrices in a level-dependent finite QBD. We examine the validity and applicability of the approximation based on the second approach and show that it gives us a practical way to obtain performance measures of call center systems in practical scale with sufficiently reasonable accuracy. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Quasi-birth–death process; Call center; Queues with impatient customers; Performance evaluation

1. Introduction Suppose that we have a multi-server queueing model. Customers are served according to the first-come first-serve (FCFS) discipline and may leave the queue before being served. In the context of queueing theory, such queueing models are commonly referred to as multi-server queues with impatient customers. In order to reveal impacts of impatient customers on performance measures, research efforts on queues with impatient customers or simply queues with impatience have been conducted over a long period. It seems that the problem of queues with impatience was first analyzed by Palm in 1937 [1]. Multi-server queues with impatience, however, have attracted much attention in the queueing literature possibly because of explosive demands to efficiently design and manage call or contact centers [2,3]. For the literature on queueing models with impatient customers as well as recent advances, see [4–8] and the references therein. In this paper, we report approximate analyses of a finite-capacity, multi-server queueing model for possible use in call centers. Our queueing model assumes that customers may leave the queue before being served. In addition, we assume that the server goes on a single vacation, which is aimed for constructing a more detailed model of the actual ∗ Tel.: +81 277 30 1838; fax: +81 277 30 1837.

E-mail address: [email protected]. 0305-0548/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2006.12.002

2464

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

service process in call centers based on telephone lines. For details, see Section 2. We assume a Markovian setting for customers’ arrival process, service and vacation times. More precisely, customers arrive according to a Poisson arrival process, and the service time of each customer is an independent and identically distributed (i.i.d.) exponential random variable, and the vacation time is also an i.i.d. exponential random variable. We, however, assume that the patience time for the individual customer in the queue is an i.i.d. random variable with a general distribution. The primary purpose of this paper is to explore a computationally practical way for analyzing such a queueing model. We approximately analyze the queueing model by using two approaches. The first approach is based on the one used by Pla et al. [9], i.e., customers’ patience distribution is approximated by the phase-type (PH) distribution. We show that the approach may induce the exponential growth of the size of the block matrices in a level-dependent finite quasi-birth–death (QBD) process and hence has a potential disadvantage as a practical tool to evaluate performance of the queueing model. The idea of the second approach is to replace the eventual reneging of customers with balking, i.e., we assume that customers whose patience time is longer than the offered virtual waiting time are immediately rejected upon arrival. Such a modeling approach allows us to use another level-dependent finite QBD process with block matrices whose size does not exponentially increase. Note that the central idea of the second approach was applied to a queueing model in cellular networks to analyze the effect of handoff area sojourn time distribution on performance [10]. We investigate how the idea works in our queueing model, which is a call center queueing model and needs an additional analysis to evaluate the balking probability because of a single vacation. Numerical examples for various customers’ patience distributions will be provided to examine the accuracy of the approximations. We show that it gives us a practical way to obtain performance measures of call center systems in practical scale with sufficiently reasonable accuracy. We, however, also find that the accuracy diminishes as the call center capacity increases, and a need to model the vacation becomes unnecessary. The rest of the paper is organized as follows. In Section 2, we briefly review the related work on queues with impatient customers, in particular, with emphasis on applications to call centers. In Section 3, we present a QBD model of a call center which is a basis for our approximations. In Section 4, we show our approximations based on two approaches and discuss the computational complexity for the modeling approaches. In Section 5, we report the results of some numerical examples of the two approximations together with computer simulation to compare the accuracy of the approximations. Finally, we summarize our results in Section 6. 2. Related work Most of papers adopted M/M/c type queues with or without impatient customers and their variants for call center queueing models, i.e., a class of qeueing models that have c identical servers, a flow of customers who arrive according to a Poisson process, and the exponentially distributed service time of each customer. While the Poisson arrival process may be an acceptable assumption if a large number of potential arrival sources can be justified, there still remains discussion about the exponential service time distribution if we consider call center queueing models. First, the service time reportedly fits the lognormal distribution, rather than an exponential one [11–13]. Another issue is the definition of the service time itself. In call centers, an agent (server) is required to finish an additional job (after-call work or wrap-up1 ), such as entering or updating data into the customer database after speaking with the customer via telephone, and the agent cannot serve the waiting customers during the wrap-up time. It may be suggested that the wrap-up time should be included in the service time itself because it remains true that the agent is not available for waiting customers during the wrap-up time. Although such a modeling approach simplifies performance analysis, the duration of the wrap-up time can also be treated as a single vacation [14,15] and is indeed expected to have impacts on performance measures such as the queue length distribution, loss probability, waiting time distribution. Note that all call centers in the real world are capacity-limited in terms of the number of agents as well as trunk lines. During the wrap-up time, a new incoming customer can occupy the line held by the customer who has left the system and hence the new customer may be queued, but will not be lost. Hence, modeling the wrap-up time as a part of the service time (e.g., M/M/c/N queue with an equivalent mean service time) cannot properly capture the effects of the after-call work, in particular, for the case where the trunk lines used in service are comparable to the number of agents. In fact, the modeling of after-call work as a single vacation was used and the effects of the wrap-up time on performance were discussed by several authors [16–19], but without impatient customers. 1 We use the terminology “after-call work” and “wrap-up” interchangeablely in this paper

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2465

According to the statistical analysis of a call center by Brown et al. [11], the coefficient of variation (CV) of the patience time distribution for the individual customer is not exactly equal to one. For example, customers who are performing stock trading or request regular service have CV greater than one but new and Internet customers have the patience time distribution with less variability. Hence, the patience time distribution may depend on the groups of customers, types of services that the customers request, and possibly others, and it is hard to identify the patience distribution uniquely. In the same literature, the hazard rate of the customers’ patience distribution was also investigated and shown to have two main peaks. The two peaks are possibly due to a “Please wait” message played by the call center system for the first time when customers enter the queue, and at the second time if customers are still in the queue. Moreover, in order to improve the grade-of-service, customers in the queue are automatically routed to another system after the deterministic maximum waiting time [6]. Therefore, the actual random deadline may be given by the minimum of the deterministic maximum waiting time and the patience time for the individual customer. Because of many uncertainties on the patience time distribution for customers as well as controls and management by the call center systems, it is unlikely that we can obtain useful performance insights by assuming an exponential distribution as a random deadline. It is highly preferable to use a general patience time distribution in call center models. We should also take into account that actual call centers are capacity-limited systems, as mentioned early in this section. By using a matrix analytic method [20] and assuming that the patience time in the queue for the individual customer is an i.i.d. PH distribution, Pla et al. [9] analyzed M/M/c/N + P H queueing model with FCFS, last-come first-served (LCFS), and service in random order (SIRO) disciplines for cellular networks. It is possible to apply the same approach to a call center queueing model with a single vacation. One possible shortcoming of the approach, however, is that the size of the matrix describing the M/M/c/N + P H model explosively increases as the buffer size increases. Brandt and Brandt [6] obtained stationary performance measures, including exact stationary queue length distribution, of M(n)/M(n)/c + GI with finite or infinite buffer model when the arrival and service rates depend on the queue length and the customers’ patience time is an i.i.d. random variable. The service mechanism related to the after-call work in call centers, however, seems not to have been considered and still needs further discussion. Relaxing the general distribution of the patience time for the individual customer, but still maintaining a single vacation modeling approach, Kawanishi [21] discussed a call center queueing model with finite capacity when the patience time for the individual customer is an i.i.d exponential random variable. As we have discussed, the exponential patience time distribution may not properly capture the actual behavior of impatient customers in call centers. 3. A call center queueing model with infinitely patient customers Before beginning an approximate analysis, we describe a multi-server queueing model with infinitely patient customers, i.e., customers do not leave the queue before being served. Suppose that we are given a c-server queueing model. The capacity N (c) of the queueing model is assumed to be finite. Poisson arrival customers with rate  are served according to the FCFS discipline. If a customer arrives and finds that the system has no remaining capacity, then that customer is blocked and lost. The service of an accepted customer immediately commences if the customer finds an available server upon arrival. Otherwise, the customer waits for service until a server becomes available. A customer leaves the system when the service of exponential random length with mean 1/ ends. Whenever a customer departs from the system, the server goes on a single vacation of exponential random length with mean 1/ that is independent from the service time and has a different parameter in general. After the single vacation ends, the server becomes available and immediately begins providing service to a customer, if there are any, in the queue. If there is no customer in the queue, the server remains available until a new customer arrives. Let us move on to the stationary analysis of the queueing model. We can identify three distinct states for each server: (1) the busy state, in which the customer is receiving the service from the server; (2) the vacation state, in which the server is on vacation after providing the service; and (3) the idle state, in which the server is idle. We denote by (n, i, j ) the state of the queueing system when the queue length is n, i servers are on vacation, j servers are in busy state, and the remaining servers (if any) are in idle state. If we denote by a ∨ b ≡ max(a, b), then the possible transitions are as follows: • from (0, i, j ) to (0, i, j + 1) with rate  for 0 i < c and 0  j < c − i; • from (0, i, j ) to (0, i + 1, j − 1) with rate j  for 0 i < c and 0 < j < c − i; • from (0, i, j ) to (0, i − 1, j ) with rate i for 0 < i c and 0 j c − i;

2466

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

• from (n, i, c − i) to (n + 1, i, c − i) with rate  for 0 n < N and [0 ∨ (c + n − N + 1)] i c; • from (n, i, c − i) to (n, i + 1, c − i − 1) with rate (c − i) for 0 n N and [0 ∨ (c + n − N − 1)]i < c; • from (n, i, c − i) to (n − 1, i − 1, c − i + 1) with rate i for 0 < n N and [0 ∨ (c + n − N )] < i c. We denote by Nq (t), Iv (t), and Ib (t), the queue length, the number of servers on vacation, and the number of servers busy at time t, respectively. Let us define the stationary probability p(i, j ) by p(i, j ) lim Pr[Nq (t) = 0, Iv (t) = i, Ib (t) = j ] t→∞

(1)

for 0 i < c and 0 j < c − i. Also, we define the stationary probability (n, i) by (n, i) lim Pr[Nq (t) = n, Iv (t) = i, Ib (t) = c − i] t→∞

(2)

for 0 n N, [0∨(c +n−N )] i c. Let p(n)=(p(0, n), p(1, n−1), . . . , p(n, 0)) for 0 n < c and (n)=((n, 0∨ (c + n − N )), . . . , (n, c)) for 0 n N . Then, the stationary probability vector  defined by (p(0), . . . , p(c − 1), (0), . . . , (N ))

(3)

is uniquely determined by solving the system of equations Q=0 with the normalization condition 1 =1 numerically [22], where 0 is a row vector of zeros, 1 is a column vector of ones, and Q is the infinitesimal generator which has the block tridiagonal form as ⎛ ⎞ A0,0 B0,0 O ··· ··· O .. ⎟ ⎜ .. ⎜ D0,1 A0,1 B0,1 . . ⎟ ⎜ ⎟ ⎜ ⎟ .. ⎜ O ⎟ . D0,2 A0,2 ⎜ ⎟ ⎜ . ⎟ . . . .. .. .. ⎜ .. ⎟ B0,c−2 ⎜ ⎟ ⎜ ⎟ Q=⎜ (4) D0,c−1 A0,c−1 B0,c−1 ⎟ ⎜ ⎟ . . .. .. ⎟ ⎜ B0 D0,c A0 ⎜ ⎟ ⎜ ⎟ .. ⎜ . D1 A1 O ⎟ ⎜ ⎟ ⎜ . ⎟ .. .. .. ⎝ .. . . . BN−1 ⎠ AN O ··· ··· O DN for N 1. In the following, we explicitly show the matrices in Q. First, let us consider the matrices that have two subscripts. When i servers are on vacation and j servers are busy for 0 i < c and 0 j < c − i, the services for arriving customers immediately commence. Hence, we have ⎛ 0 · · · 0 0⎞  . . ⎜ 0  . . . .. .. ⎟  0 0 B0,0 = ( 0), B0,1 = , . . . , B0,c−1 = ⎜ (5) .⎟ .. ⎝ .. . . ⎠, 0  0 . . 0 .. . 0 ··· 0  0 where B0,k−1 is a k × (k + 1) matrix for 1k c. Recall that the server immediately goes on a single vacation after providing a service. When the single vacation time expires, the server becomes idle if the queue remains empty. Because the vacation time is exponentially distributed with mean of 1/, we have ⎛ ⎞ 0 ··· ··· 0 .. ⎟ ⎜ ...

 . ⎟ ⎜  0 0 ⎜ 0 .. ⎟ .. D0,1 = , D0,2 =  0 , . . . , D0,c = ⎜ (6) . 0 2 . ⎟ ⎜ ⎟,  ⎜. . ⎟ 0 2 .. ... 0 ⎠ ⎝ .. 0 · · · 0 c

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2467

where D0,k is a (k + 1) × k rectangle matrix for 1 k c. Because the service time for each customer is a random variable of an exponential distribution with mean of 1/, A0,k for each 0 k c − 1 is a (k + 1) × (k + 1) square matrix given by A0,0 = (−) and ⎛  A0,1 =

 0





 , 

A0,2 =

 0 0

2  0



⎜  ⎜0 0 ⎜. .  , . . . , A0,c−1 = ⎜ ⎜. ⎜  ⎝ ... 0

(c − 1) .. . .. . ···

⎞ 0 ··· 0 .⎟ .. .. . . .. ⎟ ⎟ .. . 2 0 ⎟ ⎟, ⎟ .. .  ⎠ ··· 0 

(7)

where the diagonal elements shown by the symbol “” are arranged so that the row sums of Q are zero. We move on now to the case where none of the servers are in idle state. To simplify the description of the remaining matrices, let us consider a queueing model with infinitecapacity that corresponds to the original finite-capacity queueing ∞ ∞ model for the time being. We can show that the matrix sequences {Bn }∞ n=0 , {An }n=0 , and {Dn }n=1 in the infinite-capacity version of the queueing model all consist of (c + 1) × (c + 1) matrices. Moreover, we can see that the matrices are homogeneous in the sense that they do not depend on n. If we denote by (X)i,j the (i, j ) element of a matrix X, where the row and column indexes of a matrix start at 0, not 1, then it follows that the matrices Bn , An , and Dn are given by (Bn )i,j = i,j , (An )i,j = (c − i)i,j −1 − [(c − i) + i]i,j − i,j ,

(8) (9)

for n0, and (Dn )i,j = ii,j +1 ,

(10)

for n 1, where i,j is the Kronecker delta defined by i,j 

i = j, i  = j.

1, 0,

(11)

Transforming these matrices to those of the original queueing model with finite capacity is trivial. We give the N N N structure of the matrices in {Bn }N−1 n=0 , {An }n=0 , and {Dn }n=1 , except the diagonal elements in {An }n=0 , for c = 2 and N = 5 as an example that may be helpful for general cases

 Bn = 0 0

 An = 0 0

0 Dn =  0

0  0 2  0 0 0 2

0 0 





  0 0 for n = 0, 1, 2, B3 = , 0 , B4 =     0   , A5 = (),  for n = 0, 1, 2, 3, A4 = 0     0  0 0 , D5 = (2 0). 0 for n = 1, 2, 3, D4 = 0 2 0 0 0  0

(12)

(13)

(14)

Now we consider the steady-state performance measures. Let us denote by pB , EW q , and Wqc (0), respectively, the probability that an arriving customer is blocked and lost because no capacity remains in the system, the mean waiting time in queue, and the probability of delay for an accepted customer. Then, 1 − pB is given by the probability that a customer finds a server in idle state or a empty waiting space in the queue upon arrival; i.e., 1 − pB =

c−1

n=0

p(n)1 +

N−1

c

n=0 k=[0∨(c+n−N+1)]

(n, k).

(15)

2468

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

Using the total arrival intensity  =

c−1

n=0 p(n)B0,n 1



+

N−1 n=0

(n)Bn 1 together with Little’s formula, we have

N

1 EW q = n(n)1 . 

(16)

n=0

We can obtain 1 − Wqc (0) as the ratio of the arrival intensity of customers who are immediately served (and hence experience no delay) to the total arrival intensity; i.e., 1 − Wqc (0) =

c−1

1 p(n)1 . 

(17)

n=0

4. QBD approximations of a call center queueing model with impatient customers In this section, we discuss approximate approaches of our call center queueing model when customers may leave the queue due to impatience. Two approaches are considered. 4.1. Method 1 The first one is to approximate customers’ patience time by a PH random variable. A (continuous-time) phase-type random variable is defined as the time to the absorbing state 0 of a Markov chain on state space {0, 1, . . . , m} with initial probability vector (0 , ) and infinitesimal generator  0 0 , (18) t T where  is a row vector of size m, T is an m × m matrix, t is a column vector of size m and 0 is a row vector of zeros. We have the relation T 1 + t = 0 . We also have 0 + 1 = 1. In this paper, we assume that the patience distribution should not put probability mass at 0 and hence 0 = 0 as the service immediately commences if an arriving customer finds an available server. In order to fulfill the conditions for the infinitesimal generator of a Markov chain, T has negative diagonal elements and non-negative off-diagonal elements and t has non-negative elements. If we denote by (T )i,j as element (i, j ) of T and by ti as element i of t respectively, we have that (T )i,i < 0, Ti,j 0, ti 0,

(19)

for 1 i  = j m. We can formulate a level-dependent finite QBD [9] by keeping track of the phase variables of each waiting customer. Suppose that none of the servers are in idle state, i.e., all servers are in either busy or vacation states. Then, we can represent the state of the queueing model by (n, i, {i1 , i2 , . . . , in }) when i servers are in vacation state (thereby the remaining c − i servers are in busy state) and the phase variable for the kth customer in the queue is ik for 1 k n, given that the queue length is n > 0. By “the kth customer,” we mean that the customer is at the kth position in the queue from the head. Lastly, we represent the state of the queueing model by (0, i) when i servers are in vacation state and the remaining c − i servers are in busy state but the queue length is zero. N N ˜ N−1 ˜ N ˜ N The matrices {Bn }N−1 n=0 , {An }n=0 , and {Dn }n=1 should be replaced by the matrices {Bn }n=0 , {An }n=0 , and {Dn }n=1 as follows. Because an arriving customer who has to wait for service draws the phase-type random patience time which starts in phase i with probability i for 1i m, B˜ n for each n 0 is given by n

   B˜ n = Bn ⊗Im ⊗ Im · · · ⊗ Im ⊗,

(20)

where Ik is a k × k identity matrix and ⊗ denotes the Kronecker product operation. For D˜ n , we have to consider two possibilities. One of the two is an event of customer departure from the queue due to impatience. Let us define Un by i−1

n   

Im ⊗ Im ⊗ · · · Im ⊗ t ⊗Im ⊗ Im · · · ⊗ Im , Un     i=1

n−i

(21)

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2469

for n > 0. Then, we can see that the transition of one customer departure from the queue due to impatience, given that n customers are in the queue for n > 0, can be described by Un . The other is also an event of customer departure from the queue due to the end of a single vacation. Because the service discipline is FCFS, the first customer leaves the queue and then the service of the customer commences immediately if the queue is not empty. When the queue is empty, one of the servers in vacation state induces a transition by the end of a single vacation but the queue remains empty. Hence, Dn for n > 0 is given by n−1

   D˜ n = Dn ⊗ 1 ⊗Im ⊗ Im · · · ⊗ Im +En ⊗ Un , 

(22)

where En = Ic+1 for 1 nN − c and En = [0 IN −n+1 ] for N − c < nN . A˜ n has transition rates that do not accompany events of customer arrival or departure for n waiting customers in the queue. They describe changes in the phase variables of the phase-type patience time distributions and the number of servers in busy or vacation states. If we denote by Tn for n > 0 n−1

   Tn  T ⊕ T ⊕ · · · T ⊕ T , where ⊕ denotes the Kronecker sum operation, then A˜ n for n 0 is given by An , n = 0, ˜ An = An ⊕ Tn , 0 < nN.

(23)

(24)

While the phase-type distribution can be used to approximate a general distribution and is analytically tractable because of the Markovian properties, one phase variable is not sufficient to obtain an accurate approximation of a general distribution. Consequently, the size m of matrix T has to be greater than or equal to two. Thus, it may be hard to solve a large capacity call center queueing model because we have to deal with the matrices whose size grows exponentially as N becomes large. For example, computing the stationary probability vector  for a queueing model with capacity N by using the numerical algorithm (e.g., the linear level reduction algorithm [22] or the folding algorithm [23]) may require to compute the inverse of matrices whose size are comparable to O(mN ) such as AN . Computing the inverse of matrices with large size is costly and may happen to be numerically erroneous results, and hence it prevents us from obtaining performance measures for large value of N and m > 1. 4.2. Method 2 The second approach is based on the modeling assumption that customers who eventually leave the queue due to impatience are immediately rejected upon arrival. Because customers whose patience time is less than the offered virtual waiting time are not served anyway, such a modeling approach is expected to yield an approximate queue length process of the actual system. In this approach, we assume that an accepted customer will not leave the queue due to impatience and remains in the queue until it is served. Then, the queue length process for the approximate queueing model with impatient customers can be formulated by a level-dependent finite QBD process, i.e., a queueing model with infinitely patient customers by replacing the arrival rate  in the ith row of Bn with (1 − bn,i ), where bn,i is the probability of a customer who is immediately rejected, given that it finds n customers are in the queue and i servers are in vacation state upon arrival. It is obvious that the probability bn,i , which can be regarded as the balking probability of a customer, depends on the waiting time of the queueing model with infinitely patient customers. We can find that the waiting time is greater than t if no more than n departure events from the queue occur during the interval (0, t], given that an accepted customer has to wait for service and finds n waiting customers in the queue upon arrival for 0 nN − 1 [24,25]. Then, the conditional probability wn,i (t) of an accepted customer whose waiting time in the queue is greater than t, given that the customer finds n customers ahead in the queue and i servers are in vacation state upon arrival, can be expressed as wn,i (t) =

n

k=0

(P (k, t)1 )i ,

(25)

2470

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

where (P (·, t)1 )i is the ith element of vector P (·, t)1 and P (·, t) is a (c + 1) × (c + 1) matrix related to the number of events in terms of customer departure from the queue because of the completion of a server vacation. In view of the building blocks in Q, it is clear that P (n, t) can be identified with the probability matrix of n arrival events in (0, t] for the Markovian arrival process (MAP) with representing matrices C and D given by [19,21] (C)i,j = (c − i)i,j −1 − [(c − i) + i]i,j , (D)i,j = ii,j +1 . Using Eq. (25), we approximately evaluate bn,i for 0 n N − 1 and 0 i c by  ∞ wn,i (t) dF (t), bn,i =

(26) (27)

(28)

0

where F (t) is the distribution function of the patience time for the individual customer with mean 1/ . Note that we do not suffer from increasing the size of the matrices of the level-dependent finite QBD process in this approach, which is desirable to compute stationary performance measures numerically. We need to compute P (k, t)1 (0 k N − 1) to obtain bn,i . Computing the vector P (·, t)1 is a non-trivial task but we can reduce it to the problem of computing the exponential function of a matrix. Recall that P (n, t) satisfies the (forward) Chapman–Kolmogorov system of differential equations [26,27] d P (0, t) = P (0, t)C, dt d P (n, t) = P (n, t)C + P (n − 1, t)D, dt

(29) n1,

(30)

for t 0 with initial conditions (P (0, 0))i,j =i,j and (P (n, 0))i,j =0 for n 1. If we denote by Qn the (n+1)×(n+1) block matrix defined by ⎛C D ⎞ . .. ⎜ ⎟ C ⎟ Qn  ⎜ (31) .. ⎝ ⎠ . D C for each n0, then it follows that the matrix exponential exp[Qn t] is given by ⎛ ⎞ P (0, t) P (1, t) · · · P (n, t) .. ⎟ . ⎜ P (0, t) . . . ⎟ ⎜ exp[Qn t] = ⎜ ⎟. .. ⎝ . P (1, t) ⎠ P (0, t)

(32)

Note that we omitted the zero matrices. By exploiting the specific structure of the block matrix Qn , we can compute P (n, t) in the numerically stable way [22,28] as P (n, t) =

∞

k=0

e− t

( t)k Vk,n , k!

(33)

where  max{(−C)i,i : 0 i c} and the matrix Vk,n is computed recursively by using the non-negative matrices W C/ + Ic+1 and ZD/ as Vk,n = O, k < n, Vk,0 = W k , Vk+1,n+1 = W V k,n+1 + ZV k,n .

(34) (35) (36)

Using the algorithm to compute P (n, t) numerically, we can obtain Eq. (29). In fact, if we denote by F˜ (s) the Laplace–Stieltjes transform (LST) of a non-negative random variable with distribution function F (·), then it

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

follows that   ∞ P (n, t)1 dF (t) = 0

∞ ∞ 0

e− t

k=0

∞  ∞

( t)k Vk,n 1 dF (t) k!

( t)k Vk,n 1 dF (t) k! k=0 0  ∞ k  k

 k d ˜ = F (s) (−1)  Vk,n 1 ,  k! ds k =

2471

e− t

k=0

(37)

s=

provided that the derivatives of the LST can be computed. We can also obtain Eq. (28) by numerical integration. The approach is one of the alternative methods to compute bn,i numerically when the LST of F (·) cannot be obtained in closed form. In this case, bn,i can be approximately evaluated based on the idea that P (n, t) ≈ P (m) (n, t), =

m

e− t

k=0

( t)k Vk,n , k!

(38)

where P (m) (n, t) is the order-m approximation of P (n, t) obtained by truncating the Taylor series expansion of the matrix exponential at order m. In practice, we set  > 0 as a tolerable error in advance and choose the minimum m such that 1−

m

( t)k k=0

k!

e− t ,

(39)

for each given t. In general, it is observed that the minimum m depends on t and increases as t, which implies that computing P (m) (n, t) for given  is time consuming for large value of t. Without relying on the numerical algorithm for computing P (m) (n, t), however, we can use an explicit expression of P (n, t), e.g., spectral expression based on eigenvalues and eigenvectors. Thanks to the specific structure of C and D given by equations (26) and (27), yet another explicit expression of P (n, t) has been recently reported in [19] which may be useful to compute P (n, t) with less CPU time than when using Eq. (38), but only for small size of n. 5. Numerical examples In this section, we validate our approximate models (Methods 1 and 2) by showing numerical examples together with computer simulation comparisons. Computer simulations, which consist of 10 runs and one million call arrival events in each run, were performed to obtain the blocking probability pB of a customer who finds the queue full and hence is immediately lost upon arrival, the delay probability Wqc (0) of a typical customer accepted by the system, and the mean waiting time EW q in the queue. Note that Wqc (0) and EW q are unconditional, i.e., Wqc (0) and EW q are calculated for not only served customers but customers who leave the queue (if any) due to impatience. More explicitly, Method 1 (and computer simulation) evaluate Wqc (0) and EW q for customers who are served or may leave the queue before being served, whereas Method 2 evaluates Wqc (0) and EW q only for customers served by the system because none of the customers leave the queue before being served once they are accepted by the system. We considered several types of distribution functions F (x) under two parameter settings listed below. • • • • •

number of servers c: 60, maximum capacity N: from 60 up to 70, mean arrival time 1/: 5 and 6, mean service time 1/: 280, mean vacation time 1/: 20.

2472

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

Fig. 1. pB versus N for the patience distribution function of Case 1.1. The mean arrival time of customers is 1/ = 5.

Fig. 2. pB versus N for the patience distribution function of Case 1.1. The mean arrival time of customers is 1/ = 6.

We assume the parameter setting for a moderate scale call center (60 agents and up to 70 lines) or one site of a large-scale multi-site call center, offered by customers with average arrival time of 5 and 6 s, and average handling time (AHT) of 300 s [29]. In order to take into account the vacation time in our approximate model, we partition the AHT into the mean service time (talk length) of 280 s and mean wrap-up time (after-call work time) of 20 s. These parameter setting gives us the blocking probability (the probability of running out of space) of order one percent or less (see figures later), which is a practically reasonable provisioning of quality of service in terms of the loss probability for actual call centers. Note that Method 1 becomes impossible to evaluate numerically for c = 60, N 60 due to the exponential growth of sizes of block matrices, if the number of phases of the PH distribution is more than one.

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2473

Fig. 3. pB versus N for the patience distribution function of Case 1.2. The mean arrival time of customers is 1/ = 5.

Fig. 4. pB versus N for the patience distribution function of Case 1.2. The mean arrival time of customers is 1/ = 6.

5.1. Case 1: patience distributions with common mean The first set of our numerical examples is obtained for two distributions but the distributions share the same mean deadline 1/ = 90, i.e., Case 1.1, 1 − 1{x<1/ } (40) F (x) = 1 − 1{x<100} e−x Case 1.2, where 1{·} is the indicator function and  is a unique solution of the non-linear equation 1 − e−100 = 90, which gives  ≈ 0.00214556. The possible scenario for Case 1.1 is for a customer who is sufficiently patient but controlled by the

2474

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

Fig. 5. EW q versus N for the patience distribution function of Case 1.1. The mean arrival time of customers is 1/ = 5.

Fig. 6. EW q versus N for the patience distribution function of Case 1.1. The mean arrival time of customers is 1/ = 6.

system, e.g., automatically routed to another system, after the deterministic maximum waiting time. Case 1.2 is for a customer whose patience time distribution can be assumed to be exponential but he/she is automatically controlled by the system after the deterministic maximum waiting time. Figures 1–4 (5–8) show pB (EW q ) versus the maximum capacity N, and Tables 1 and 2 show Wqc (0) for Case 1.1 and Case 1.2. We applied Method 2 and also performed computer simulation with 95% confidence interval.2 Together with these results, we computed “non-vacation model,” i.e., M/M/c/N + G queueing model in which the mean of the single exponential service time is tuned to 1/ + 1/, and Method 1 in case of the exponential patience distribution 2 The error bars are too narrow and hence almost invisible. It is the same for the subsequent figure examples.

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2475

Fig. 7. EW q versus N for the patience distribution function of Case 1.2. The mean arrival time of customers is 1/ = 5.

Fig. 8. EW q versus N for the patience distribution function of Case 1.2. The mean arrival time of customers is 1/ = 6.

with mean 1/ = 90. Note the capacity N includes the servers as well as the waiting spaces. Hence, the non-vacation model reduces to the loss system (Erlang-B model) if N = c, for which case the patience time distribution is irrelevant to the performance because no customers can wait for service. On the other hand, it is possible for customers to wait in our queueing model even if N = c by virtue of modeling the after-call work as a single vacation. The “balking probability” in Method 2 is computed by using Eq. (38) in Case 1.1 and by numerical integration in Case 1.2. The third, fourth, fifth, and sixth columns in Tables 1 and 2, respectively, show the numerical results of Method 2, and the sample mean, the upper and lower limits of the 95% confidence interval based on computer simulation, for Wqc (0). We listed numerical values of Wqc (0) for the non-vacation model and the exponential patience distribution in the seventh and eighth columns of Tables 1 and 2, respectively.

2476

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

Table 1 Wqc (0) for the patience distribution function of Case 1.1. c = 60 is assumed Simulation 1/

N

Method 2

Mean

Upper

Lower

Non-vacation

Exponential

5 5 5 5 5 5 5 5 5 5 5

60 61 62 63 64 65 66 67 68 69 70

0.2956 0.3507 0.3964 0.4356 0.4697 0.4998 0.5265 0.5502 0.5711 0.5891 0.6042

0.29579 0.35053 0.39542 0.43431 0.46927 0.49853 0.52663 0.55109 0.57010 0.58933 0.60290

0.29582 0.35056 0.39548 0.43438 0.46932 0.49862 0.52669 0.55114 0.57019 0.58946 0.60303

0.29576 0.35050 0.39536 0.43424 0.46923 0.49844 0.52657 0.55104 0.57001 0.58920 0.60277

N/A 0.0963 0.1756 0.2422 0.2988 0.3475 0.3899 0.4271 0.4600 0.4893 0.5155

0.2360 0.3000 0.3258 0.3439 0.3565 0.3652 0.3709 0.3747 0.3770 0.3785 0.3794

6 6 6 6 6 6 6 6 6 6 6

60 61 62 63 64 65 66 67 68 69 70

0.0604 0.0713 0.0798 0.0866 0.0920 0.0965 0.1001 0.1030 0.1053 0.1071 0.1084

0.06079 0.07074 0.07929 0.08664 0.09252 0.09687 0.10008 0.10366 0.10506 0.10746 0.10826

0.06082 0.07076 0.07932 0.08670 0.09257 0.09692 0.10014 0.10385 0.10513 0.10749 0.10833

0.06077 0.07071 0.07925 0.08659 0.09247 0.09681 0.10001 0.10347 0.10500 0.10742 0.10819

N/A 0.0217 0.0390 0.0530 0.0644 0.0736 0.0812 0.0874 0.0925 0.0967 0.1002

0.0542 0.0615 0.0662 0.0691 0.0709 0.0720 0.0726 0.0729 0.0731 0.0732 0.0733

Table 2 Wqc (0) for the patience distribution function of Case 1.2. c = 60 is assumed Simulation 1/

N

Method 2

Mean

Upper

Lower

Non-vacation

Exponential

5 5 5 5 5 5 5 5 5 5 5

60 61 62 63 64 65 66 67 68 69 70

0.2851 0.3349 0.3747 0.4074 0.4348 0.4581 0.4780 0.4951 0.5098 0.5227 0.5338

0.28750 0.33818 0.37873 0.41279 0.43981 0.46233 0.48177 0.49973 0.51074 0.52391 0.53055

0.28754 0.33823 0.37877 0.41284 0.43987 0.46239 0.48182 0.49979 0.51079 0.52395 0.53069

0.28745 0.33813 0.37868 0.41274 0.43975 0.46226 0.48171 0.49966 0.51069 0.52387 0.53041

N/A 0.0963 0.1749 0.2396 0.2933 0.3381 0.3757 0.4073 0.4340 0.4565 0.4755

0.2360 0.3000 0.3258 0.3439 0.3565 0.3652 0.3709 0.3747 0.3770 0.3785 0.3794

6 6 6 6 6 6 6 6 6 6 6

60 61 62 63 64 65 66 67 68 69 70

0.0581 0.0678 0.0749 0.0803 0.0844 0.0875 0.0898 0.0916 0.0929 0.0939 0.0947

0.05916 0.06951 0.07753 0.08276 0.08550 0.08976 0.09158 0.09329 0.09535 0.09518 0.09685

0.05919 0.06954 0.07755 0.08280 0.08554 0.08981 0.09161 0.09333 0.09539 0.09522 0.09692

0.05914 0.06948 0.07750 0.08272 0.08546 0.08972 0.09155 0.09325 0.09530 0.09514 0.09679

N/A 0.0217 0.0388 0.0524 0.0631 0.0714 0.0780 0.0830 0.0869 0.0898 0.0920

0.0542 0.0615 0.0662 0.0691 0.0709 0.0720 0.0726 0.0729 0.0731 0.0732 0.0733

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2477

Fig. 9. pB versus cX for the lognormal patience distribution. The mean arrival time of customers is 1/ = 5. The maximum capacity N is equal to c (=60) and pB in case of the non-vacation model is computed by using the Erlang-B formula.

Fig. 10. pB versus cX for the lognormal patience distribution. The mean arrival time of customers is 1/ = 6. The maximum capacity N is equal to c (=60) and pB in case of the non-vacation model is computed by using the Erlang-B formula.

We can see that the approximate queueing model by Method 2 gives pB and EW q close to the one obtained by computer simulation for both Case 1.1 and Case 1.2, while the non-vacation model cannot capture the effects of the after-call work in call centers. In particular, the non-vacation model overestimates pB , underestimates EW q , and the deviations are not minor. We can also observe in Tables 1 and 2 that Wqc (0) can be computed with high accuracy by the approximate model. In contrast, the non-vacation model underestimates Wqc (0) compared with computer simulation. Note that the non-vacation model is unable to compute the delay-related performance measures in case of N = c because it reduces the Erlang-loss system, while our approximate queueing model can evaluate customers’ delayrelated performance measures even in case of N = c. It should also be emphasized that the mean duration of the after-call work (1/ = 20) is less than one-tenth of the mean talk length (1/ = 280) in our numerical examples. These findings suggest that such a small fraction of the after-call work impacts the performance measures that cannot be

2478

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

Fig. 11. EW q versus cX for the lognormal patience distribution. The mean arrival time of customers is 1/ = 5. The maximum capacity N is equal to c (=60) and hence EW q is not available in case of the non-vacation model.

Fig. 12. EW q versus cX for the lognormal patience distribution. The mean arrival time of customers is 1/ = 6. The maximum capacity N is equal to c (=60) and hence EW q is not available in case of the non-vacation model.

ignored in some cases. Consequently, the after-call work should be incorporated as a vacation, not equivalent service time, into the queueing model for dimensioning call centers accurately. We can also find that the distribution of the random patience time significantly affects the performance measures possibly due to the influence on the queue length process. Numerical examples suggest that using the exponential distribution based on one-parameter fitting (by tuning to the same mean patience time) cannot properly capture the system performance. Hence, the patience distribution significantly affects the system performance of call centers. The final remark is that Method 2 tends to deviate from computer simulation as N increases and its impacts become larger for Case 1.2 than Case 1.1. We need to take into account (at least) two things: the size of the system capacity and the variability of the patience distribution. In case of large value of N, we can expect that the queue length process is well approximated by M/M/c/N + G with the equivalent service time, not our approximate queueing model by Method 2,

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2479

Table 3 Wqc (0) for the lognormal patience distribution with mean 1/ = 90 in case of c = N = 60 Simulation 1/

cX

Method 2

Mean

Upper

Lower

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

0.2956 0.2956 0.2956 0.2955 0.2951 0.2942 0.2926 0.2902 0.2871 0.2835 0.2794 0.2751 0.2707 0.2663 0.2619

0.29517 0.29473 0.29620 0.29485 0.29384 0.29517 0.29402 0.29280 0.28919 0.28745 0.28267 0.28128 0.27836 0.27434 0.27160

0.29524 0.29476 0.29626 0.29491 0.29388 0.29522 0.29407 0.29285 0.28926 0.28752 0.28272 0.28133 0.27840 0.27441 0.27163

0.29511 0.29469 0.29615 0.29479 0.29381 0.29513 0.29398 0.29276 0.28913 0.28738 0.28261 0.28122 0.27832 0.27428 0.27156

6 6 6 6 6 6 6 6 6 6 6 6 6 6 6

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

0.0604 0.0604 0.0604 0.0603 0.0603 0.0601 0.0598 0.0593 0.0586 0.0579 0.0570 0.0561 0.0552 0.0542 0.0533

0.06100 0.06120 0.06090 0.06062 0.06029 0.06042 0.05951 0.05934 0.05891 0.05818 0.05821 0.05736 0.05654 0.05669 0.05597

0.06103 0.06123 0.06092 0.06067 0.06031 0.06043 0.05952 0.05936 0.05893 0.05820 0.05823 0.05737 0.05656 0.05671 0.05599

0.06098 0.06118 0.06088 0.06058 0.06028 0.06040 0.05949 0.05931 0.05889 0.05816 0.05819 0.05734 0.05652 0.05667 0.05595

because the vacation time will become almost irrelevant for the system performance when the capacity is sufficiently large. We discuss the issue of the variability of the patience time distribution in the next examples. 5.2. Case 2: lognormal patience distribution We assume that the patience time for the individual is a random variable with lognormal distribution given by √ customer x F (x) = ((log(x) − )/), where (x) = (1/ 2) −∞ exp(−u2 /2) du is the standard normal probability distribution function. We use the lognormal distribution as an example for a class of customers who are sufficiently patient but 2 the mean and the squared CV of the lognormal distribution, then not infinitely patient. If we denote by E[X] and cX 2 2 2 we have E[X] = exp( +  /2) and cX = exp( ) − 1 [25]. The two parameters of the lognormal distribution were adjusted so that E[X] = 90 and given cX . It is well known that the lognormal distribution can be approximated by Erlang or (two-stage) hyperexponential distributions [25], and we can obtain numerical results by using a level-dependent finite QBD modeling approach of Method 1. We, however, should notice again that we need more than one phase to approximate the lognormal distribution, and it is practically difficult to obtain performance measures in a realistic system size of call center like our numerical parameter setting. Thus, we evaluated numerical examples of approximation by using Method 2 only. Because the lognormal distribution does not seem to have the LST in closed form, we evaluated bn,i by numerical integration only in case N = c. We show numerical examples of pB and EW q versus cX in Figures 9–12, and Wqc (0) in Table 3. Recall that EW q and Wqc (0) are not available for the non-vacation model in case of N =c, and we compute pB by using the Erlang-B formula.

2480

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

The first observation is a satisfactory agreement between the approximate queueing model by Method 2 and computer simulation for pB within our numerical examples of cX . In contrast, the non-vacation model (Erlang-B formula) largely overestimates pB and independent of cX because of the irrelevance of the patience distribution, which is not observed by computer simulation. For EW q and Wqc (0), Method 2 also gives a satisfactory agreement with computer simulation. We, however, find that the deviations between the approximation and computer simulation tend to grow as cX becomes large. These findings may imply that our approximate queueing model by Method 2 gives accurate results for small value, but not for large value, of CV. Besides the accuracy with respect to the variability of the patience distribution, however, we can find again that incorporating the after-call work as a single vacation allows us to evaluate the system performance of call centers in a unified and satisfactory way. 6. Conclusion In this paper, we have reported approximate analyses of the effect of impatient customers in a finite-capacity, multiserver queueing model, where each server goes on a single vacation. We have assumed that the customers’ patience time in the queue is an i.i.d. random variable with a general distribution. Two approximate methods have been analyzed. The first one (Method 1) is based on using the phase-type distribution to approximate the general patience distribution. The idea of the approximation of the second approach (Method 2) is that customers whose random patience time eventually expires are immediately rejected upon arrival, which was applied to a queueing model in cellular networks [11] and proved to be a successful tool for performance analysis. We have shown that such a modeling approach can be applied to a call center queueing model that has a more complex service mechanism than that of the cellular networks. Method 1 may induce the exponential growth of the size of the block matrices in the level-dependent finite QBD process and hence has a potential disadvantage as a practical tool to evaluate performance of the queueing model. By comparison, Method 2 allows us to use a level-dependent finite QBD process with block matrices whose size does not exponentially increase. We have investigated how the approximate queueing models work and provided numerical examples to show the accuracy of the approximations. Our numerical examples suggest that the approximate analysis by the level-dependent finite QBD approach (Method 2) gives us a practical way to obtain performance measures of call center systems in practical scale with sufficiently reasonable accuracy. We, however, have also found that the accuracy diminishes as the call center capacity increases, and a need to model the vacation becomes unnecessary. Acknowledgment The author thanks the referees and editor for their critical reading of the manuscript and useful suggestions and comments that have significantly improved the quality of the paper. References [1] Palm C. Etude des delais d’attente. Ericsson Technics 1937;5(2):39–56. [2] Garnett O, Mandelbaum A, Reiman M. Designing a call center with impatient customers. Manufacturing & Service Operations Management 2002;4(3):208–27. [3] Koole G, Mandelbaum A. Queueing models of call centers: an introduction. Annals of Operations Research 2002;113:41–59. [4] Harris CM, Hoffman KL, Saunders PB. Evaluation and modeling of the IRS telephone taxpayer information system. Operations Research 1987;35:504–23. [5] Boxma O, de Waal P. Multiserver queues with impatient customers. Proceedings of the 14th international teletraffic congress, Antibes Juan-les-Pins; 1994 p. 743–56. [6] Brandt A, Brandt M. On the M(n)/M(n)/s queue with impatient calls. Performance Evaluation 1999;35:1–18. [7] Aksin OZ, Harker PT. Modeling a phone center: analysis of a multichannel, multiresource processor shared loss system. Management Science 2001;47(2):324–36. [8] Brandt A, Brandt M. Asymptotic results and a Markovian approximation for the M(n)/M(n)/s + GI system. Queueing Systems 2002;41: 73–94. [9] Pla V, Casares-Giner V, Matinez J. On a multiserver finite buffer queue with impatient customers. Proceedings of ITC specialist seminar, Antwerp; 2004. [10] Pla V, Casares-Giner V. Effect of the handoff area sojourn time distribution on the performance of cellular networks. Proceedings of IEEE mobile and wireless communications network (MWCN), City; 2002. p. 401–405.

K. Kawanishi / Computers & Operations Research 35 (2008) 2463 – 2481

2481

[11] Brown L, Gans N, Mandelbaum A, Sakov A, Shen H, Zeltyn S, Zhao L. Statistical analysis of a telephone call center: a queueing-science perspective. http://econpapers.repec.org/paper/woppennin/03-12.htm. [12] Gans N, Koole G, Mandelbaum A. Telephone call centers: tutorial, review, and research prospects. Manufacturing and Service Operations Management 2003;5:79–141. [13] Mandelbaum A, Sakov A, Zeltyn S. Empirical analysis of a call center. Technical Report, 2000, available at http://iew3.technion. ac.il/serveng/References/references.html. [14] Doshi BT. Queueing system with vacation—a survey. Queueing Systems 1986;1:29–66. [15] Takagi H. Queueing analysis: foundation of performance evaluation, vacation and priority systems, Part 1. Amsterdam: North Holland; 1991. [16] Feinberg MA. Performance characteristics of automated call distribution systems. Proceedings of the IEEE conference on global telecommunications (GLOBECOM), San Diego; 1990. 402A.4. [17] Fischer MJ, Garbin DA, Gharakhanian A. Performance modeling of distributed automatic call distribution systems. Telecommunications Systems 1998;9(2):133–52. [18] Jolley WM, Harris RJ. Analysis of post-call activity in queueing systems. Proceedings of the 9th international teletraffic congress, Torremolinos; 1979. p. 1–9. [19] Kawanishi K. On the counting process for a class of Markovian arrival processes with an application to a queueing system. Queueing Systems 2005;49(2):93–122. [20] Neuts MF. Matrix-geometric solutions in stochastic models: an algorithmic approach. Baltimore: Johns Hopkins University Press; 1981. [21] Kawanishi K. A closed-form solution on a level-dependent Markovian arrival process with queueing application. Calcolo 2004;41(3):153–75. [22] Latouche G, Ramaswami V. Introduction to matrix analytic methods in stochastic modeling. Philadelphia: SIAM Society for Industrial & Applied Mathematics; 1999. [23] Ye J, Li S-Q. Folding algorithm: a computational method for finite QBD processes with level-dependent transitions. IEEE Transactions on Communications 1994;43:625–39. [24] Gross D, Harris CM. Fundamentals of queueing theory. New York: Wiley; 1998. [25] Tijms HK. Stochastic models: an algorithmic approach. England: Wiley; 1994. [26] Lucantoni DM, Meier-Hellstern KS, Neuts MF. A single-server queue with server vacations and a class of nonrenewal arrival processes. Advances in Applied Probability 1990;22:676–705. [27] Neuts MF. Structured stochastic matrices of M/G/1 type and their applications. New York: Marcel Dekker; 1989. [28] Neuts MF, Li J-M. An algorithm for the P (n, t) matrices of a continuous BMAP. In: Chakravarthy SR, Alfa AS, editors. Matrix-analytic methods in stochastic models. New York: Marcel Dekker; 1997. p. 7–19. [29] Zeltyn S, Mandelbaum A. Call centers with impatient customers: many-server asymptotics of the M/M/n + G queue. Queueing Systems 2005;51(3–4):361–402.