Queueing system MAP|PH|N|N+R with impatient heterogeneous customers as a model of call center

Applied Mathematical Modelling 37 (2013) 958–976

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Queueing system MAPjPHjNjN þ R with impatient heterogeneous customers as a model of call center Chesoong Kim a,⇑, Sergey Dudin b, Olga Taramin b, Janghyun Baek c a

Sangji University, Wonju, Kangwon 220-702, Republic of Korea Belarusian State University, 4, Nezavisimosti Ave., Minsk 220030, Belarus c Chonbuk National University, Jeonju, Republic of Korea b

a r t i c l e

i n f o

Article history: Received 4 August 2011 Received in revised form 2 March 2012 Accepted 12 March 2012 Available online 21 March 2012 Keywords: Call center Multi-server queueing system Markovian arrival process Phase type service time distribution Impatient heterogeneous customers

a b s t r a c t A multi-server queueing system with a Markovian arrival process (MAP), a finite buffer and impatient heterogeneous customers useful in modeling a call center is investigated. The servers are identical and independent of each other. The service time of a customer has a Phase type distribution (PH). If all servers are busy at a customer arrival epoch, the customer moves to the buffer with probability, which depends on the number of customers in the system, or, alternatively, leaves the system forever. During a waiting period the customers can be impatient and leave the system. An efficient algorithm for calculating the stationary probabilities of system states is proposed. A special approach for reducing the system dimension is used. Some key performance measures are calculated. The Laplace– Stieltjes transforms of the sojourn and waiting time distribution are derived. Some illustrative numerical results are presented. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction A call center is a centralised office used by companies for receiving and servicing its clients requests by telephone. Call centers are playing an increasingly important role in today’s world. According to [1], all fortune 500 companies have at least one call center. The continued growth of call centers has brought a big interest from practitioners and academic researchers. To describe the operation of call center queueing theory is used. Mathematical models allow us to solve a wide range of optimization problems, e.g. the problem of determining the optimal number of call agents or customer service representatives (CSRs) needed to serve clients with a preassigned level of quality of service. This task has a great practical importance, because most of the costs of call center involves wages with benefits (see, e.g. [2]). The simplest model of a call center is the Erlang-C system (MjMjN queueing system). In this model a stationary Poisson arrival process and exponentially distributed service times are assumed. Also, customers’ abandonment and impatience are ignored. Another simple model of a call center is MjMjNjN queueing system (Erlang-B system). This model assumes no buffer. The customers are assumed to be lost when all servers are busy. In [3], MjMjNjR queueing system with finite buffer and impatient customers as a model of a call center is analyzed. Customers’ retrial in a call center is considered in [4]. In [5], a queueing model with two classes of customers, with and without abandonment is analyzed, and the problem of announcing delays to customers upon their arrival is considered. In [6], tandem queueing models of two-stage call centers are investigated. These models are useful for modelling the situation when a customer after servicing at the first stage of a call center

⇑ Corresponding author. E-mail addresses: [email protected] (C. Kim), [email protected] (S. Dudin), [email protected] (O. Taramin), [email protected] (J. Baek). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.03.021

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

959

can require additional service. For the references and the present state of art in investigation of call centers the reader is referred to the paper [2]. Models of call centers in the overwhelming majority of the existing papers assume that the arrival flow of customers is Poisson and the service time is exponentially distributed. These assumptions greatly simplify the study of such systems, but at the same time reduce the adequacy of the model. Our statistical analysis of real arrival flows and service processes of the call center of one of the largest banks in Belarus has shown that arrival flows are correlated and the service time does not have exponential distribution. Similar conclusions about the arrival process and the service time in real call centers are reported in many papers, for references see, e.g. [7]. In this paper, we consider a queueing model which describes the operation of the call center with different types of impatient customers. This model is much more general than the existing ones in the several direction. In order to study the effect of variation and correlation within the arrival process as well as the correlation of the arrivals of different types of customers, we use Markovian arrival process (MAP) to model the arrivals of customers. Also, another interesting feature of the presented paper is the assumption of phase type (PH) services which generalizes the exponential services assumed in most of the existing papers in the literature. This allows to one take into account the variation of the service time. The fact that the service time distribution in modern telephone networks is poorly approximated by an exponential distribution, but is well described by a PH distribution is reported in [8,9]. The impact of the service time distribution, specifically the variability, on the system performance measures is illustrated numerically in this paper. It is well known that if we consider a multi-server queueing system with PH service process and the number of servers is large, the Markovian process describing the system has a huge state space. This is due to the fact that at a given moment we must know the state of PH underlying process for each busy server. In the existing papers, which deal with the multi-server queues with PH service process, the maximal number of servers does not exceed 10 in the numerical examples. It means that using the standard technique we can not calculate performance measures of a call center with large number of CSRs (greater than 10). To avoid this, we use the special method proposed by Ramaswami and Lucantoni, see [10,11], for reduction the dimension of systems with phase type service time distribution. Using this reduction technique enables one to study the model with large number of CSRs. The model under study, as shown below, will also help practitioners to optimize the call centers economically. The paper is organized as follows. In Section 2, the mathematical model is described. The stationary distribution of the system states is analyzed in Section 3. The expressions for the main performance measures of the system are given in Section 4. In Sections 5 and 6, we focus on the analysis of the sojourn and waiting time distributions of different type customers. Illustrative numerical examples are given in Section 7, and some concluding remarks are given in Section 8.

2. Mathematical model We consider a queueing system to model a call center with N CSRs (servers), a finite waiting space (buffer) of capacity R and two types of calls, referred to as priority and non-priority calls (customers). The customers arrive to the system according to MAP. Generally speaking, MAP is correlated, so it is ideal to model correlated and bursty traffic in modern telecommunication networks and call centers in particular, see, e.g., [12,13]. Different methods of the estimation of the MAP parameters using a finite set of observed data such as a set of customer arrival times recorded at the real call center are presented, e.g., in paper [14]. The arrival of customers in the MAP is directed by the stochastic process mt ; t P 0. The process mt ; t P 0, is an irreducible continuous time Markov chain with the state space f0; 1; . . . ; Wg. The sojourn time of this chain in the state m is exponentially ð0Þ distributed with the positive finite parameter km . When the sojourn time in the state m expires, with probability pm;m0 the proðrÞ cess mt jumps to the state m0 without generation of customers, m; m0 2 f0; 1; . . . ; Wg; m – m0 , and with probability pm;m0 the process mt jumps to the state m0 with generation of type r customer, r ¼ 1; 2; m; m0 2 f0; 1; . . . ; Wg. ðrÞ ðrÞ The behavior of the MAP is completely characterized by the matrices D0 ; D1 ; r ¼ 1; 2, defined by entries ðD1 Þm;m0 ¼ ðrÞ ð0Þ km pm;m0 ; m; m0 2 f0; 1; . . . ; Wg; r ¼ 1; 2, and ðD0 Þm;m ¼ km ; m 2 f0; 1; . . . ; Wg; ðD0 Þm;m0 ¼ km pm;m0 ; m; m0 2 f0; 1; . . . ; Wg; m – m0 . The matrix ð1Þ ð2Þ Dð1Þ ¼ D0 þ D1 þ D1 represents the generator of the process mt ; t P 0. ð1Þ ð2Þ The average total arrival rate k is given as k ¼ hðD1 þ D1 Þe where h is the invariant vector of a stationary distribution of the Markov chain ðmt ; t P 0Þ. The vector h is the unique solution to the system hDð1Þ ¼ 0; he ¼ 1. Here e is a column vector of appropriate size consisting of 1’s and 0 is a row vector of appropriate size consisting of zeroes. The average arrival rate kr of ðrÞ type r customers is given as kr ¼ hD1 e; r ¼ 1; 2. The squared integral (without differentiating the types of customers) coefficient of variation cv ar of intervals between sucðrÞ cessive arrivals is given as cv ar ¼ 2khðD0 Þ1 e  1. The squared coefficient of variation cv ar of inter-arrival times of type r cusðrÞ ðr Þ 1 tomers is given as cv ar ¼ 2kr hðD0  D1 Þ e  1; r – r; r ; r ¼ 1; 2. The integral coefficient of correlation ccor of two successive ðrÞ intervals between arrivals is given as ccor ¼ ðkhðD0 Þ1 ðDð1Þ  D0 ÞðD0 Þ1 e  1Þ=cv ar . The coefficient of correlation ccor of two ðrÞ ðr Þ ðrÞ successive intervals between type r customers’ arrivals is computed by ccor ¼ ðkr hðD0 þ D1 Þ1 D1 ðr Þ 1 ðrÞ ðD0 þ D1 Þ e  1Þ=cv ar ; r – r. More information about the MAP with different types of customers and related research is given, e.g. in [15,16].

960

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

We assume that type 1 customers are priority customers and type 2 customers are non-priority customers. Priority customers are more important for the call center, so they have the advantages of access speed to CSR, but the behavior of different type customers after beginning the talk is assumed to be statistically identical, i.e., each CSR serves incoming calls with the same efficiency independently of customers priority. This assumption is usual in call center literature, see motivation, e.g. in [5]. We also suggest that the call center queue is ‘‘visible’’ (see, e.g. [17]), which means the following. An arriving customer, who cannot enter into service immediately, is informed about the queue length. The customer then decides either to leave immediately due to the length of the queue or join the queue first and possibly leave after some random time by being impatient. So, in mathematical model we assume that if at an arbitrary customer arrival epoch there is a free server, a customer is admitted to the system and occupies the free server. If at a non-priority customer arrival epoch all servers are busy with i; i 2 f0; 1; . . . ; R  1g, customers are presenting in the buffer then this non-priority customer leaves the system with probability qi or moves to the buffer with supplementary probability. If the non-priority customer moves to the buffer, it is placed at the position number i þ 1 in the buffer. If at a priority customer arrival epoch all servers are busy, there are i; i 2 f0; 1; . . . ; R  1g, customers in the buffer, and the number of priority customers in the buffer is k; k 2 f0; 1; . . . ; ig, then this priority customer leaves the system with probability pk or moves to the buffer with supplementary probability. If the priority customer moves to the buffer, it is placed at the position number k þ 1 in the buffer, i.e., the priority customer is placed after priority customers and ahead of non-priority customers already presenting in the buffer. If at an arbitrary customer arrival epoch the buffer is full, the customer leaves the system forever. The priority and non-priority customers can be impatient, i.e., the type r customer leaves the system after arrival in random time, which is exponentially distributed with the parameter ar ; 0 < ar < 1, if this customer is not served in between. The service time (talk time and after talk work time) of a customer by each server has a PH distribution with an irreducible representation ðb; SÞ. This service time can be interpreted as a time until the underlying Markov process gt ; t P 0, with a finite state space f1; . . . ; M; M þ 1g reaches the single absorbing state M þ 1 conditioned on the fact that the initial state of this process is selected among the states f1; . . . ; Mg according to the probabilistic row vector b ¼ ðb1 ; . . . ; bM Þ. Transition rates of the process gt within the set f1; . . . ; Mg are defined by the sub-generator S and transition rates into the absorbing state (which lead to the service completion) are given by the entries of the column vector S0 ¼ Se. R1 The service time distribution function has the form AðxÞ ¼ 1  beSx e, Laplace–Stieltjes transform (LST) 0 esx dAðxÞ of this 1 1 distribution is bðsI  SÞ S0 ; Res > 0. The mean service time is calculated by b1 ¼ bðSÞ e. The squared coefficient of varia2 tion is given by cv ar ¼ b2 =b1  1 where b2 ¼ 2bðSÞ2 e. For more information about PH distribution and its usefulness see, e.g. [18]. The methods of modelling of PH process using a set of service times obtained at the real system, and in particular call center, can be founded in paper [19]. Note that the main goal of the paper is to develop the algorithms for computing the key performance measures of the system under consideration. Once these measures are calculated, we can solve various optimization problems relevant to call center management in practice. 3. The process of system states Let it be the number of customers in the system, it 2 f0; 1; . . . ; N þ Rg; kt be the number of priority customers in the buffer, ðmÞ kt 2 f0; 1; . . . ; maxf0; it  Ngg; mt be the state of the directing process of the MAP; mt 2 f0; 1; . . . ; Wg; gt be the number of P ðmÞ ðmÞ servers at the phase m of service, m 2 f1; 2; . . . ; Mg; gt 2 f0; 1; . . . ; minfit ; Ngg; M g ¼ minfi t ; Ng, at the epoch m¼1 t ðmÞ t; t P 0. Note that the meaning of the components gt is chosen according to the approach by Ramaswami and Lucantoni, see [10,11]. So, the behavior of the system under consideration can be described in terms of the regular irreducible continuous-time Markov chain ð1Þ

ðMÞ

nt ¼ fit ; kt ; mt ; gt ; . . . ; gt g;

t P 0;

with the state space

    fi; 0; m; gð1Þ ; . . . ; gðMÞ g; i 2 f0; 1; . . . ; Ng [ fi; k; m; gð1Þ ; . . . ; gðMÞ g; i 2 fN þ 1; N þ 2; . . . ; R þ Ng ; k 2 f0; 1; . . . ; i  Ng;

m 2 f0; 1; . . . ; Wg; gðmÞ 2 f0; 1; . . . ; minfi; Ngg; m 2 f1; 2; . . . ; Mg:

 NþM1 . In the overM1 whelming majority of the existing papers, the behavior of the queues with PH service time distribution is described by means ð1Þ

ðMÞ

The dimension of the process gt ¼ fgt ; . . . ; gt g; t P 0, when all servers are busy is equal to K ¼ ð1Þ

ðNÞ



ðmÞ

of stochastic process including components ft ; . . . ; ft ; t P 0, where ft ; m 2 f1; 2; . . . ; Ng, is the state of the PH underlying process on the mth busy server. Note that the dimension K is significantly less than the dimension of the process fft ; . . . ; ft g; t P 0, which is equal to K ¼ M N . For example, if we fix N ¼ 50 and M ¼ 2, the dimension K ¼ 250 while the dimension K is only 51. ð1Þ

ðNÞ

961

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

Since the Markov chain nt is regular irreducible and has a finite state space, then for any choice of the system parameters there exist stationary probabilities of the system states which are defined as follows:



n



o

ð1Þ ðMÞ p i; 0; m; gð1Þ ; . . . ; gðMÞ ¼ t!1 limP it ¼ i; kt ¼ 0; mt ¼ m; gt ¼ gð1Þ ; . . . ; gt ¼ gðMÞ ;

i 2 f0; 1; . . . ; Ng;

m 2 f0; 1; . . . ; Wg; gðmÞ 2 f0; 1; . . . ; ig; m 2 f1; 2; . . . ; Mg;

M P

gðmÞ ¼ i;

m¼1



n



o

ð1Þ ðMÞ p i; k; m; gð1Þ ; . . . ; gðMÞ ¼ t!1 limP it ¼ i; kt ¼ k; mt ¼ m; gt ¼ gð1Þ ; . . . ; gt ¼ gðMÞ ;

i 2 fN þ 1; N þ 2; . . . ; N þ Rg;

k 2 f0; 1; . . . ; i  Ng; M P

gðmÞ 2 f0; 1; . . . ; Ng; m 2 f1; 2; . . . ; Mg;

m 2 f0; 1; . . . ; Wg;

gðmÞ ¼ N:

m¼1

Let us form the row vectors pði; 0; mÞ of the probabilities pði; 0; m; gð1Þ ; . . . ; gðMÞ Þ; i 2 f0; 1; . . . ; Ng, and the row vectors pði; k; mÞ of the probabilities pði; k; m; gð1Þ ; . . . ; gðMÞ Þ; i 2 fN þ 1; N þ 2; . . . ; N þ Rg; k 2 f0; 1; . . . ; i  Ng, enumerated in the reverse lexicographic order of the components gð1Þ ; . . . ; gðMÞ . Then let us form the row vectors pi :

pi ¼ ðpði; 0; 0Þ; pði; 0; 1Þ; . . . ; pði; 0; WÞÞ; i 2 f0; 1; . . . ; Ng; pði; kÞ ¼ ðpði; k; 0Þ; pði; k; 1Þ; . . . ; pði; k; WÞÞ; k 2 f0; 1; . . . ; i  Ng; pi ¼ ðpði; 0Þ; pði; 1Þ; . . . ; pði; i  NÞÞ; i 2 fN þ 1; N þ 2; . . . ; N þ Rg: It is well-known that the probability vectors pi ; i 2 f0; 1; . . . ; N þ Rg, satisfy the following system of linear algebraic equations:

ðp0 ; p1 ; . . . ; pNþR ÞQ ¼ 0;

ðp0 ; p1 ; . . . ; pNþR Þe ¼ 1;

ð1Þ

where Q is the infinitesimal generator of the Markov chain nt ; t P 0. Lemma 1. The infinitesimal generator Q of the Markov chain nt ; t P 0, has the block-three-diagonal structure:

0Q

0;0

B Q 1;0 B B B O B Q ¼B . B .. B B @ O O

Q 0;1

O

...

O

Q 1;1

Q 1;2

...

O

Q 2;1 .. . O

Q 2;2 .. . O

... O .. .. . . . . . Q NþR1;NþR1

O

O

...

Q NþR;NþR1

O

1

C C C O C C C: .. C . C C Q NþR1;NþR A O

Q NþR;NþR

The non-zero blocks Q i;j ; i; j P 0, have the following form:

Q i;i ¼ D0  Ai ðN; SÞ þ IW  DðiÞ ;

i 2 f0; 1; . . . ; N  1g;

  ð1Þ ð2Þ Q N;N ¼ D0  AN ðN; SÞ þ p0 D1 þ q0 D1  IK N þ IW  DðNÞ ;     ð2Þ ð1Þ e iNþ1 þ a2 C b iNþ1  I Q i;i ¼ IiNþ1  D0  AN ðN; SÞ þ qiN D1  IK N þ IW  DðNÞ þ Q iNþ1  D1  IK N  a1 C WK N ; i 2 fN þ 1; N þ 2; . . . ; N þ R  1g;     e Rþ1 þ a2 C b Rþ1  I Q NþR;NþR ¼ IRþ1  Dð1Þ  AN ðN; SÞ þ IW  DðNÞ  a1 C ; WK N   Q i;i1 ¼ IW  LNi N; e S ;

i 2 f1; 2; . . . ; Ng;

    e iNþ1 E b b Q i;i1 ¼ EiNþ1  IW  L0 ðN; e SÞPN1 ðbÞ þ a1 C iNþ1 þ a2 C iNþ1 E iNþ1  IWK N ; Q i;iþ1 ¼ D1  Pi ðbÞ; Q i;iþ1 ¼

i 2 fN þ 1; N þ 2; . . . ; N þ Rg;

i 2 f0; 1; . . . ; N  1g;

   ð1Þ ð2Þ IiNþ1  Q iNþ1 EþiNþ1  D1 þ ð1  qiN Þ b E þiNþ1  D1  IK N ;

i 2 fN; N þ 1; . . . ; N þ R  1g;

962

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

where  I is an identity matrix, O is a zero matrix of appropriate dimension;   and  are symbols of Kronecker’s sum and product respectively, see, e.g. [20];  W ¼W þ1 ;   0 0  e S¼ ; S0 S ð1Þ ð1Þ  D1 ¼ D1 þ D2 ;   

e l ¼ diagf0; 1; . . . ; l  1g; C b l ¼ diagfl  1; l  2; . . . ; 0g; l 2 f2; 3; . . . ; R þ 1g; C Q l ¼ diagfp0 ; p1 ; . . . ; pl1 g; l 2 f1; 2; . . . ; Rg; b E l ; E l ; l 2 f2; 3; . . . ; R þ 1g, are the matrices of size l  ðl  1Þ which are defined as

1 0 1 1 0 0  0 1 0  0 B0 1  0C B1 0 0  0C C B C B C  B. . . B .. C B .. .. 0 1 0    0 C; b . El ¼ B ¼ E . .C C; B C l B .. C C B B .. .. .. . . @0 0  1A @. . . . .A 0

0 

Eþ l ;

b Eþ l ;

0

0 0  1

0  0

l 2 f1; 2; . . . ; Rg, are the matrices of size l  ðl þ 1Þ which are defined as

0

0 1 0

B0 0 B Eþl ¼ B B .. .. @. . 

0 0

 0 0

1

0

1 0  0 0

1

C B 1  0 0C C þ B0 1  0 0C b C; B. . . ¼ ; E .. . . . .. .. C . C l B. . . . .. .. C . . . .A A @. . 0

0

 0 1

0  1 0

 iþM1  Ki ¼ ; i 2 f0; 1; . . . ; Ng ; M1 ðiÞ  D ¼ diagfAi ðN; SÞe þ LNi ðN; e SÞeg; i 2 f1; 2; . . . ; Ng; Dð0Þ ¼ O11 . SÞ; i 2 f0; . . . ; Ng, The algorithms for calculation of the matrices P i ðbÞ; i 2 f0; . . . ; N  1g; Ai ðN; SÞ; i 2 f0; . . . ; Ng, and LNi ðN; e are presented in Appendix A. The proof of the lemma is implemented by means of the analysis of all transitions of the Markov chain nt ; t P 0, during the interval of an infinitesimal length and rewriting intensities of these transitions in the block matrix form. Note that the matrix Pi ðbÞ defines the transition probabilities of the process gt ; t P 0, at the epoch of starting the new service given that i servers are busy at this epoch. The matrix LNi ðN; e SÞ defines the intensities of transitions of this process at the service completion epoch given that i servers are busy at this epoch. The matrix Ai ðN; SÞ defines the intensities of transitions of the process gt ; t P 0, which do not lead to the service completion given that i servers are busy. The modules of the diagonal entries of the matrix DðiÞ define the total intensity of leaving the corresponding states of the process gt ; t P 0, given that i servers are busy. If the dimension of the system (1) is small, it can be easily solved on a computer by standard methods. Otherwise, to solve this system the following numerically stable algorithm can be used. Theorem 1. The vectors pi ; i 2 f0; 1; . . . ; N þ Rg, are given as follows

pi ¼ pi1 T i1 ¼ p0 F i ; i 2 f1; 2; . . . ; N þ Rg; where the matrices F i are calculated using the recurrent formulas:

F 0 ¼ I;

F i ¼ F i1 T i1 ;

i 2 f1; 2; . . . ; N þ Rg;

the matrices T i ; i 2 f0; 1; . . . ; N þ R  1g, are calculated using the backward recurrence

T i ¼ Q i;iþ1 ðQ iþ1;iþ1 þ T iþ1 Q iþ2;iþ1 Þ1 ;

i ¼ N þ R  2; N þ R  3; . . . ; 0;

under the initial condition

T NþR1 ¼ Q NþR1;NþR ðQ NþR;NþR Þ1 ; the vector p0 is the unique solution to the system

p0 ðQ 0;0 þ T 0 Q 1;0 Þ ¼ 0; p0

NþR P

F l e ¼ 1:

l¼0

The numerical stability of the proposed algorithm follows from the fact, that all inverted matrices computed in this algorithm are irreducible sub-generators. So, as it is well known from the matrix theory, these inverse matrices exist and are nonnegative.

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

963

4. Performance measures As soon as the vectors pi ; i 2 f0; 1; . . . ; N þ Rg, have been calculated, we are able to find various performance measures of the system (call center) under consideration. The stationary distribution of the number of customers in the system lim Pfit ¼ ig ¼ pi e; i 2 f0; 1; . . . ; N þ Rg. t!1 P The average number of customers in the system e L ¼ NþR i¼1 ipi e. P The average number of customers in the buffer N buffer ¼ NþR i¼Nþ1 ði  NÞpi e. P PiN ¼ NþR The average number of priority customers in the buffer N buffer i¼Nþ1 k¼1 kpði; kÞe. 1 P NþR PiN1 buffer ¼ i¼Nþ1 k¼0 ði  k  NÞpði; kÞe. The average number of non-priority customers in the buffer N 2 P The average number of busy servers N serv er ¼ NþR i¼1 minfi; Ngpi e. The loss probability of an arbitrary customer at the entrance to the call center due to buffer overflow is given as

Pentloss ¼ k1 pNþR ðIRþ1  D1  IK N Þe: The loss probability of an arbitrary type r customer at the entrance to the system due to buffer overflow is calculated as follows

  ðrÞ Pentloss ¼ k1 r pNþR IRþ1  D1  I K N e; r

r ¼ 1; 2:

The probability P escloss that an arbitrary customer arrives when all servers are busy, buffer is not full, and the customer does not join the buffer and leaves the system is given as

Pescloss ¼ k1

NþR1 P i¼N

h



i

pi Q iNþ1  D1ð1Þ þ qiN IiNþ1  Dð2Þ  IK N e: 1

Pescloss , r

The probabilities r ¼ 1; 2, that an arbitrary type r customer arrives when all servers are busy, the buffer is not full, and the customer does not join the buffer and leaves the system are given as

Pescloss ¼ k1 1 1

¼ k1 Pescloss 2 2

NþR1 P i¼N NþR1 P i¼N





pi Q iNþ1  D1ð1Þ  IK N e;   ð2Þ qiN pi IiNþ1  D1  IK N e:

The intensity of flow of customers, which get the service in the system, is calculated as

kout ¼

i P h Nserv er NþR ¼ pi Imaxf1;iNþ1g  IW  LmaxfNi;0g ðN; eSÞ e: b1 i¼1

The loss probability of an arbitrary customer is calculated as

Ploss ¼ 1 

kout : k

The probability P imploss that an arbitrary customer after arrival will go to the buffer and leave it due to impatience

Pimploss ¼ Ploss  Pentloss  P escloss : Let us introduce the probabilities zðk; gð1Þ ; . . . ; gðMÞ Þ; k 2 f1; 2; . . . ; Rg, that during the waiting time of a priority customer in the buffer this customer does not leave the call center due to impatience conditioned on the fact that at its arrival epoch ð1Þ ðMÞ there are k  1 priority customers in the buffer and the states of the processes gt ; . . . ; gt are gð1Þ ; . . . ; gðMÞ correspondingly. The column vectors zðkÞ consisting of the probabilities zðk; gð1Þ ; . . . ; gðMÞ Þ enumerated in the reverse lexicographic order of the components gð1Þ ; . . . ; gðMÞ are given as follows:

zð1Þ ¼ ða1 I  AÞ1 LeK N ; zðkÞ ¼ ðka1 I  AÞ1 ðL þ ðk  1Þa1 IÞzðk  1Þ;

k 2 f2; 3; . . . ; Rg;

where

A ¼ AN ðN; SÞ þ DðNÞ ;

  L ¼ L0 N; e S PN1 ðbÞ:

The probability P imploss that an arbitrary priority customer will go to the buffer and leave it due to impatience is given by 1

964

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

Pimploss ¼ k1 1 1 The probability lated as

NþR1 P iN P i¼N k¼0

Pimploss 2

  ð1Þ ð1  pk Þpði; kÞ D1 e  IK N ðeK N  zðk þ 1ÞÞ:

that an arbitrary non-priority customer will go to the buffer and leave it due to impatience is calcu-

imploss Pimploss ¼ k1  k1 Pimploss Þ: 2 ðkP 2 1

The proof of this formula follows from the expression kPimploss ¼ k1 P imploss þ k2 P imploss . 1 2 The loss probability of an arbitrary priority customer is given as

Ploss ¼ Pimploss þ Pentloss þ P1escloss : 1 1 1 The loss probability of an arbitrary non-priority customer is given by loss Ploss ¼ Pimploss þ Pentloss þ P2escloss ¼ k1  k1 Ploss 2 ðkP 2 2 1 Þ: 2

The intensity of flow of type r customers, which get the service in the system, is calculated as

kout ¼ kr ð1  Ploss r r Þ;

r ¼ 1; 2:

5. Distribution of the sojourn and waiting times of a priority customer Let V 1 ðxÞ be the distribution function of the sojourn time of an arbitrary priority customer in the system under study and R 1 sx e dV 1 ðxÞ; Re s > 0, be its LST. 0

v 1 ðsÞ ¼

Theorem 2. The LST

v 1 ðsÞ is calculated as follows:

v 1 ðsÞ ¼ Pentloss þ Pescloss þ k1 1 1 1 where the vectors

v 1 ðs; kÞ;

N1     NþR1 N P P iP ð1Þ pi D1ð1Þ  IK i ebðsI  SÞ1 S0 þ ð1  pk Þpði; kÞ D1 e  IK N v 1 ðs; k þ 1Þ ; i¼0

i¼N k¼0

k 2 f1; 2; . . . ; Rg, are given as

v1 ðs; 1Þ ¼ ððs þ a1 ÞI  AÞ1 ½LeK

N

bðsI  SÞ1 S0 þ a1 eK N ;

v1 ðs; kÞ ¼ ððs þ ka1 ÞI  AÞ1 ½ðL þ ðk  1Þa1 IÞv1 ðs; k  1Þ þ a1 eK

ð2Þ N

;

k 2 f2; 3; . . . ; Rg:

ð3Þ

Proof. Let us tag an arbitrary priority customer and keep track of its service process in the system. We will derive the expression for the LST v 1 ðsÞ by means of the method of collective marks (method of additional event, method of catastrophes) for references, see, e.g., [21,22]. To this end, we interpret the variable s as the intensity of some virtual stationary Poisson flow of catastrophes. So, v 1 ðsÞ has the meaning of probability that no one catastrophe arrives during the sojourn time of an arbitrary priority customer. The following situations are possible at the arrival epoch of the tagged customer:  The customer leaves the system immediately upon arrival. The probability of this event can be calculated by the formula P1entloss þ P escloss . In this case, the probability that no one catastrophe arrives during the sojourn time is equal to one. 1  There is an idle server at the arrival epoch of the tagged customer and this customer immediately goes to the service. The PN1  ð1Þ probability of this event is k1 1 i¼0 pi D1  IK i e. In this case, the probability that no one catastrophe arrives during the sojourn time is equal to the probability that no one catastrophe arrives during the service time and given as bðsI  SÞ1 S0 .  All servers are busy at the arrival epoch and the tagged customer attends the buffer. Let v 1 ðs; k; gð1Þ ; . . . ; gðMÞ Þ be the LST of distribution of the tagged priority customer’s sojourn time conditioned on the fact that, at the given moment, the position of the tagged customer in the buffer is equal to k; k 2 f1; 2; . . . ; Rg, and the states of the ð1Þ ðMÞ processes gt ; . . . ; gt ; t P 0, are gð1Þ ; . . . ; gðMÞ correspondingly. Let us form the vectors v 1 ðs; kÞ of these LSTs enumerated in the reverse lexicographic order of the components gð1Þ ; . . . ; gðMÞ . If we find the vectors v 1 ðs; kÞ then the probability that the tagged customer will be admitted to the system at the epoch when all server are busy and no one catastrophe arrives during its sojourn time will be written as:

k1 1

NþR1 N P iP i¼N k¼0

  ð1Þ ð1  pk Þpði; kÞ D1 e  IK N v 1 ðs; k þ 1Þ:

Based on a probabilistic sense of the LST, we obtain the system for calculation of the vectors

v 1 ðs; kÞ:

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

965

ððs þ ka1 ÞI þ AÞv 1 ðs; kÞ þ ð1  dk;1 ÞðL þ ðk  1Þa1 IÞv 1 ðs; k  1Þ þ a1 eK N þ dk;1 LeK N bðsI  SÞ1 S0 ¼ 0T ; k 2 f1; 2; . . . ; Rg; ð4Þ where di;j is symbol of Kronecker’s delta. It is easy to verify that the solution of (4) can be written as (2) and (3). Using the formula of total probability we prove the theorem. h Having the LST v 1 ðsÞ been computed we can determine the distribution of an arbitrary priority customer’s sojourn time using numerical methods, see, e.g. [23]. From the practical point of view it is interesting to know the average sojourn time of an arbitrary priority customer in the system, which can be easily found by means of the LST v 1 ðsÞ. Corollary 1. The average sojourn time of an arbitrary priority customer in the system is given as:

 N1     NþR1 P P P iN ð1Þ ð1Þ 1 V soj b1 pi D1  IK i e  ð1  pk Þpði; kÞ D1 e  IK N v 01 ð0; k þ 1Þ ; 1 ¼ k1 i¼0

ð5Þ

i¼N k¼0

where

v01 ð0; 1Þ ¼ ða1 I  AÞ1 ðI þ b1 LÞeK v01 ð0; kÞ ¼ ðka1 I  AÞ1



N

;

 eK N  ðL þ ðk  1Þa1 IÞv 01 ð0; k  1Þ ;

k 2 f2; 3; . . . ; Rg:

Proof. Formula (5) for calculation of the average sojourn time of an arbitrary priority customer in the system is based on 0 definition V soj 1 ¼ v 1 ðsÞjs¼0 taking into account that v ðs; kÞjs¼0 ¼ eK N ; k 2 f1; 2; . . . ; Rg. By analogy with Theorem 2 we can derive expressions for LSTs of distribution of an arbitrary priority customer’s waiting time (Theorem 3) and of an arbitrary serviced priority customer’s sojourn time (Theorem 4). h Theorem 3. The LST y1 ðsÞ of distribution of an arbitrary priority customer’s waiting time can be found as follows:

y1 ðsÞ ¼ Pentloss þ Pescloss þ k1 1 1 1

N1     NþR1 P iN P P ð1Þ pi D1ð1Þ  IK i e þ ð1  pk Þpði; kÞ D1 e  IK N y1 ðs; k þ 1Þ i¼0

i¼N k¼0

where the vectors y1 ðs; kÞ; k 2 f1; 2; . . . ; Rg, are calculated as

y1 ðs; 1Þ ¼ ððs þ a1 ÞI  AÞ1 ðL þ a1 IÞeK N ; y1 ðs; kÞ ¼ ððs þ ka1 ÞI  AÞ1 ½ðL þ ðk  1Þa1 IÞy1 ðs; k  1Þ þ a1 eK N ;

k 2 f2; 3; . . . ; Rg:

Corollary 2. The average waiting time of an arbitrary priority customer is calculated as follows:

V wait ¼ k1 1 1

NþR1 P iN P

  ð1Þ ð1  pk Þpði; kÞ D1 e  IK N y01 ð0; k þ 1Þ;

i¼N k¼0

where

y01 ð0; 1Þ ¼ ða1 I  AÞ1 eK N ; y01 ð0; kÞ ¼ ðka1 I  AÞ1 ½eK N  ðL þ ðk  1Þa1 IÞy01 ð0; k  1Þ;

k 2 f2; 3; . . . ; Rg:

Theorem 4. The LST w1 ðsÞ of distribution of the sojourn time of an arbitrary serviced priority customer in the system is calculated as 1 w1 ðsÞ ¼ ðkout 1 Þ

N1     NþR1 P iN P P ð1Þ 1 pi Dð1Þ  I S þ ð1  p Þ p ði; kÞ D e  I ðs; k þ 1Þ ; ebðsI  SÞ w Ki 0 KN 1 k 1 1 i¼0

i¼N k¼0

where the vectors w1 ðs; kÞ; k 2 f1; 2; . . . ; Rg, are computed by formulas

w1 ðs; 1Þ ¼ ððs þ a1 ÞI  AÞ1 LeK N bðsI  SÞ1 S0 ; w1 ðs; kÞ ¼ ððs þ ka1 ÞI  AÞ1 ðL þ ðk  1Þa1 IÞw1 ðs; k  1Þ;

k 2 f2; 3; . . . ; Rg:

966

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

Corollary 3. The average sojourn time of an arbitrary serviced priority customer in the system is given as follows:

 N1   NþR1 P P P iN ð1Þ ð1Þ v out 1 0 V ser ¼ ðk Þ b p D  I ð1  p Þ p ði; kÞðD e  I Þw ð0; k þ 1Þ e  1 i K K k 1 1 1 N 1 1 i i¼0

i¼N k¼0

where

w01 ð0; 1Þ ¼ ða1 I  AÞ1 ½ða1 I  AÞ1 þ b1 ILeK N ; w01 ð0; kÞ ¼ ðka1 I  AÞ1 ½ðka1 I  AÞ1 ðL þ ðk  1Þa1 IÞw1 ð0; k  1Þ  ðL þ ðk  1Þa1 IÞw01 ð0; k  1Þ;

k 2 f2; 3; . . . ; Rg:

v v Corollary 4. The average waiting time of an arbitrary serviced priority customer is calculated as V waitser ¼ V ser  b1 . 1 1

6. Distribution of the sojourn and waiting times of a non-priority customer Derivation of the formula for calculation of the LST of distribution of an arbitrary non-priority customer’s sojourn time in the system v 2 ðsÞ is more involved. We will derive expression for the LST v 2 ðsÞ by means of the method of collective marks. So, v 2 ðsÞ has the meaning of probability that no one catastrophe arrives during the sojourn time of an arbitrary non-priority customer. We will tag an arbitrary non-priority customer and will keep track of its staying in the system. Let v 2 ðs; i; l; k; m; gð1Þ ; . . . ; gðMÞ Þ be the probability that catastrophe will not arrive during the rest of the tagged customer sojourn time in the system conditioned on the fact that, at the given moment, the number of customers in the buffer is equal to i; i 2 f1; 2; . . . ; Rg, the tagged customer has position number l; l 2 f1; 2; . . . ; ig, in the buffer, the number of priority customers in the buffer is equal ð1Þ ðMÞ to k; k 2 f0; 1; . . . ; i  1g, and the states of the processes mt ; gt ; . . . ; gt ; t P 0, are m; gð1Þ ; . . . ; gðMÞ ; m 2 f0; 1; . . . ; Wg; gðmÞ 2 PM ðmÞ f0; 1; . . . ; Ng; m¼1 g ¼ N; m 2 f1; 2; . . . ; Mg. Let us enumerate probabilities v 2 ðs; i; l; k; m; gð1Þ ; . . . ; gðMÞ Þ in the reverse lexicographic order of the components gð1Þ ; . . . ; gðMÞ and direct lexicographic order of the component m and form from these probabilities the column vectors v 2 ðs; i; l; kÞ. It follows from the formula of total probability that if we will have functions v 2 ðs; i; l; kÞ been calculated the LST v 2 ðsÞ is computed by

v 2 ðsÞ ¼ Pentloss þ Pescloss þ k1 2 2 2 NþR1 P

þ

i¼N

ð1  qiN Þ

iP N k¼0

N1   P pi D1ð2Þ  IK i ebðsI  SÞ1 S0 i¼0



pði; kÞðD1ð2Þ e  IK N Þv 2 ðs; i  N þ 1; i  N þ 1; kÞ :

The system of linear algebraic equations for vectors

v2 ðs; i; l; kÞ ¼ ½ðs þ ka1 þ ði  kÞa2 ÞI  D0  A

v2 ðs; i; l; kÞ is derived by means of the formula of total probability:

1 

dk;0 ðð1  dl;1 ÞIW  Lv 2 ðs; i  1; l  1; 0Þ   þ dl;1 IW  LbðsI  SÞ S0 Þ þ ð1  dk;0 Þ IW  L þ ka1 I v 2 ðs; i  1; l  1; k  1Þ 1

þ ði  lÞa2 v 2 ðs; i  1; l; kÞ þ ðl  1  kÞa2 v 2 ðs; i  1; l  1; kÞ þ a2 eWK N h ð1Þ ð2Þ þ ð1  di;R Þ ð1  pk ÞD1  IK N v 2 ðs; i þ 1; l þ 1; k þ 1Þ þ ð1  qi ÞD1  IK N v 2 ðs; i þ 1; l; kÞ i  ð1Þ ð2Þ þðpk D1 þ qi D1 Þ  IK N v 2 ðs; i; l; kÞ þ di;R D1  IK N v 2 ðs; R; l; kÞ ; i 2 f1; 2; . . . ; Rg; l 2 f1; 2; . . . ; ig; k 2 f0; 1; . . . ; l  1g:

ð6Þ

Let us explain formula (6) in brief. The diagonal entries of the matrix in the first square brackets in (6) are equal to the total intensity of the events which can happen after an arbitrary epoch: catastrophe arrival, transition of the directing process of the MAP, transition of any directing process of the PH service processes, and abandon of priority and non-priority customers from the buffer. The non-diagonal entries of the matrix in the first square brackets in (6) are equal to the intensities of the events which do not impact on the components s; i; l and k (transition of the MAP underlying process without generation of customer, transition of one of PH underlying processes which does not lead to service completion). The first two terms in the round brackets in (6) correspond to the case when the number of priority customers in the buffer is equal to 0 and one service process is completed. In this case the number of customers and position of the tagged customer in the buffer are decreased by one. The number bðsI  SÞ1 S0 defines the probability that catastrophe will not arrive during the service time of the tagged customer. Third term corresponds to the case when the number of priority customers in the buffer is greater than 0 and service process in any server is completed or one priority customer leaves the buffer due to impatience. In this case not only the number of customers and position of the tagged customer, but also the number of priority customers in the buffer is decreased by one. Fourth, fifth and sixth terms correspond to the case when a non-priority customer leaves the buffer due to impatience. The rest of terms corresponds to the case when a new customer arrives to the system. It is worth to note, that if a priority customer arrives and is admitted to the system, the position of the tagged customer is increased by one.

967

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

Let us introduce the column vectors

v2 ðs; i; lÞ ¼ ðv 2 ðs; i; l; 0Þ; . . . ; v2 ðs; i; l; l  1ÞÞT ; v 2 ðs; iÞ ¼ ðv2 ðs; i; 1Þ; . . . ; v 2 ðs; i; iÞÞT ; v2 ðsÞ ¼ ðv2 ðs; 1Þ; . . . ; v 2 ðs; RÞÞT : System (6) of the linear algebraic equations can be rewritten into the matrix form as

b lþ1 v 2 ðs; i þ 1; l þ 1Þ þ ð1  dl;1 Þ Q b l1 v 2 ðs; i  1; l  1Þ  ðsI  Q li;i Þv 2 ðs; i; lÞ þ Q li;iþ1 v 2 ðs; i þ 1; lÞ þ Q i;i1 i;iþ1 b l v 2 ðs; i  1; lÞ þ aðlÞ ðsÞ ¼ 0T ; þ ð1  dl;i Þ Q i;i1

i 2 f1; 2; . . . ; Rg; l 2 f1; 2; . . . ; ig;

ð7Þ

where

   ð1Þ ð2Þ el  I  IK N ; Q li;i ¼ Il  ðD0  AÞ  a1 C  a2 C il  IWK N þ di;R Il  D1 þ ð1  di;R Þ Q l  D1 þ qi Il  D1 WK N

l

2 f1; 2; . . . ; ig; i 2 f1; 2; . . . ; Rg; Q li;iþ1 ¼ ð1  qi ÞIl  D1  IK N ; ð2Þ

b lþ1 ¼ ðIl  Q l ÞEþ  Dð1Þ  IK ; Q l i;iþ1 N 1

l 2 f1; 2; . . . ; ig; i 2 f1; 2; . . . ; R  1g; l 2 f1; 2; . . . ; ig; i 2 f1; 2; . . . ; R  1g;

  e l E þ a2 C blb b l1 ¼ E  I  L þ a1 C E l  IWK N ; Q i;i1 l l W Q li;i1 ¼ ði  lÞa2 Il  IWK N ;

b lþ1 ¼ O; l 2 f1; 2; . . . ; Rg; Q R;Rþ1

b l1 ¼ O; l 2 f2; 3; . . . ; ig; i 2 f2; 3; . . . ; Rg; Q 1;0

i 2 f2; 3; . . . ; Rg; Q l1;0 ¼ O; l 2 f1; 2; . . . ; i  1g;

C il ¼ diagfi; i  1; . . . ; i  l þ 1g;

l 2 f1; 2; . . . ; ig;

aðlÞ ðsÞ ¼ a2 elWK N þ dl;1 IW  LeWK N bðsI  SÞ1 S0 ;

l 2 f1; 2; . . . ; ig; i 2 f1; 2; . . . ; Rg:

Let us introduce the following notation: 0 1 X1;1 X1;2 O ... O O B X2;1 X2;2 X2;3 . . . O O C B C B O X X . . . O O C 3;2 3;3 B C  X¼B . where . . . . . .. C .. .. .. .. B .. C B C @ O O O . . . XR1;R1 XR1;R A O O O . . . XR;R1 XR;R

Xi;i ¼ diagfQ 1i;i ; . . . ; Q ii;i g; 0

Xi;i1

B BQ b1 B i;i1 B B O B ¼B B . B .. B B B O @ O 0

Xi;iþ1

Q 1i;i1

Q 1i;iþ1

B B O B ¼B B . B .. @ O

i 2 f1; 2; . . . ; Rg;

O

O



O

Q 2i;i1

O



O

b2 Q i;i1

Q 3i;i1



O

.. . O

..

. O

..

. 

.. . b i2 Q i;i1

O

O



O

b2 Q i;iþ1

O



O

Q 2i;iþ1

b3 Q i;iþ1



O

.. .

..

..

.. .

O

O

.

.

   Q ii;iþ1

O

1

C O C C C O C C C; .. C . C C C C Q i1 i;i1 A i1 b Q i;i1 O

i 2 f2; 3; . . . ; Rg;

1

C O C C C; .. C . C A

i 2 f1; 2; . . . ; R  1g:

b iþ1 Q i;iþ1

 Vector aðsÞ is given as aðsÞ ¼ ða1 ðsÞ; . . . ; aR ðsÞÞT where ai ðsÞ ¼ ðað1Þ ðsÞ; . . . ; aðiÞ ðsÞÞT ; i 2 f1; 2; . . . ; Rg.Using this notation we can rewrite the system (7) into the form

ðX  sIÞv 2 ðsÞ þ aðsÞ ¼ 0T :

ð8Þ

968

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

It can be verified that the diagonal entries of the matrix X  sI dominate in all rows of this matrix. So the inverse matrix exists. Thus we have proved the following assertion. Theorem 5. The vector v 2 ðsÞ consisting of the conditional LST v 2 ðs; i; l; k; m; gð1Þ ; . . . ; gðMÞ Þ; i 2 f1; 2; . . . ; Rg; l 2 f1; 2; . . . ; ig; k 2 f0; 1; . . . ; l  1g; m 2 f0; 1; . . . ; Wg; gðmÞ 2 f0; 1; . . . ; ig; m 2 f1; 2; . . . ; Mg, is calculated by

v2 ðsÞ ¼ ðX  sIÞ1 aðsÞ:

ð9Þ

Formula (9) gives the explicit form of the vector v 2 ðsÞ, but in practice the matrix X  sI usually has a big dimension. Using the fact that this matrix has block tridiagonal form the vector v 2 ðsÞ can be easily found from system (8) using the tridiagonal matrix algorithm or Thomas algorithm, see [24], that is modified for block tridiagonal matrices, see Appendix. Having the probabilities

v 2 ðs; i; l; k; m; gð1Þ ; . . . ; gðMÞ Þ been computed we can find LST v 2 ðsÞ.

Corollary 5. The average sojourn time V soj 2 of an arbitrary non-priority customer is calculated by

 N1     @ v ðs; i  N þ 1; i  N þ 1; kÞ NþR1 iN P P P 2 ð2Þ ð2Þ 1 V soj b1 pi D1  IK i e  ð1  qiN Þ pði; kÞ D1  IK N js¼0 : 2 ¼ k2 @s i¼0 i¼N k¼0 Here the column vectors @ v2 ðs;iNþ1;iNþ1;kÞ js¼0 are calculated as the blocks of the vector dvds2 ðsÞ js¼0 ¼ X1 ½a0 ð0Þ  v 2 ð0Þ. @s By analogy with Theorem 5 we can derive expressions for LSTs of distribution of an arbitrary non-priority customer’s waiting time (Theorem 6) and arbitrary serviced non-priority customer’s sojourn time (Theorem 7). Theorem 6. The LST y2 ðsÞ of distribution of an arbitrary non-priority customer’s waiting time can be found as follows:

y2 ðsÞ ¼ Pentloss þ P escloss þ k1 2 2 2

N1  P i¼0



pi D1ð2Þ  IK i e þ

NþR1 P

iN P

i¼N

k¼0

ð1  qiN Þ





pði; kÞ D1ð2Þ  IK N y2 ðs; i  N þ 1; i  N þ 1; kÞ



~. where the column vectors y2 ðs;i  N þ 1; i  N þ 1; kÞ are calculated as the blocks of the vector y2 ðsÞ ¼ ðsI  XÞ1 a ~ R ÞT ; a ~i ¼ a ~ð1Þ ; . . . ; a ~ðiÞ T ; a ~ðlÞ ¼ a2 elWK þ dl;1 IW  LeWK ; l 2 f1; 2; . . . ; ig; i 2 f1; 2; . . . ; Rg. ~ ¼ ða ~1 ; . . . ; a Here a N N of an arbitrary non-priority customer is calculated by Corollary 6. The average waiting time V wait 2

V wait ¼ k1 2 2

NþR1 P

ð1  qiN Þ

i¼N

iN P

k¼0



pði; kÞ Dð2Þ 1  IK N

 @y ðs; i  N þ 1; i  N þ 1; kÞ 2 js¼0 : @s

Theorem 7. The LST w2 ðsÞ of distribution of the sojourn time of an arbitrary serviced non-priority customer in the system is calculated as 1 w2 ðsÞ ¼ ðkout 2 Þ

N1     NþR1 iN P P P ð2Þ pi D1ð2Þ  IK i ebðsI  SÞ1 S0 þ ð1  qiN Þ pði; kÞ D1  IK N w2 ðs; i  N þ 1; i  N þ 1; kÞ i¼0

i¼N

k¼0

^ðsÞ. where the column vectors w2 ðs; i  N þ 1; i  N 1; kÞ are calculated of the vector w2 ðsÞ ¼ ðsI  XÞ1 a þ T ðlÞas the blocks T T ð1Þ ðiÞ ^ ^ ^ ^ ^ ^ ^ Here aðsÞ ¼ ða1 ðsÞ; . . . ; aR ðsÞÞ ; ai ðsÞ ¼ a ðsÞ; . . . ; a ðsÞ ; a ðsÞ ¼ 0lWK þ dl;1 IW  LeWK N bðsI  SÞ1 S0 ; l 2 f1; 2; . . . ; ig; N i 2 f1; 2; . . . ; Rg:

v Corollary 7. The average sojourn time V ser of an arbitrary serviced non-priority customer is computed by 2

 N1     @w ðs; i  N þ 1; i  N þ 1; kÞ NþR1 iN P P P 2 ð2Þ ð2Þ 1 v js¼0 : V ser ¼ ðkout b1 pi D1  IK i e  ð1  qiN Þ pði; kÞ D1 e  IK N 2 Þ 2 @s i¼0 i¼N k¼0 v v Corollary 8. The average waiting time of an arbitrary serviced non-priority customer is calculated as V waitser ¼ V ser  b1 . 2 2

7. Numerical Examples As it was mentioned in Section 1, analysis of the data obtained from the call center of one of the largest banks in Belarus has shown that an arrival flow of calls is correlated. In Experiment 1, we intend to show that the correlation in the input process has a great impact on the performance measures of the system and traditional in literature assumption that the arrival flow is described by the stationary Poisson process can lead to pure prediction of performance measures of the call center operation.

969

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976 ð1Þ

ð2Þ

Below we introduce four MAPs defined by the matrices D0 , D1 and D1 . All these MAPs have the same average total arrival rate k ¼ 11:5, the average intensity of priority customers k1 ¼ 1:5, and the average intensity of non-priority customers k2 ¼ 10, but different integral coefficients of correlation. ð1Þ ð2Þ The first process coded as MAP0 is defined by the matrices D0 ¼ 11:5; D1 ¼ 1:5 and D1 ¼ 10. It has the coefficients of correlation ccor ¼ ccor ¼ 0 and the coefficients of variation cv ar ¼ cv ar ¼ 1; r ¼ 1; 2. The second process coded as MAP0:2 is defined by the matrices ðrÞ

D0 ¼



39:29454 0:01169

ðrÞ

     5:08073 0:02314 33:80460 0:1548 ð1Þ ð2Þ ; D1 ¼ ; D1 ¼ : 1:27440 0:01829 0:14639 0:12199 0:97604 0:23127

ð1Þ

ð2Þ

ð1Þ

ð2Þ

It has the coefficients of correlation ccor ¼ 0:20; ccor ¼ 0:39 and the coefficients of variation cv ar ¼ 7:11; cv ar ¼ 12:13. The fourth process coded as MAP0:4 2 also has the coefficient of correlation c cor ¼ 0:4, but is defined by the matrices

D0 ¼



39:25439 0:01168

     1:84191 0:31141 37:00106 0:10001 ð1Þ ð2Þ ; D1 ¼ ; D1 ¼ : 1:27343 0:13905 1:12116 0:00102 0:00052 0

ð1Þ

ð2Þ

ð1Þ

ð2Þ

It has the coefficients of correlation ccor ¼ 0:037; ccor ¼ 0:001 and the coefficients of variation cv ar ¼ 1:12; cv ar ¼ 96:56. 0:4   MAP0:2 ; MAP0:4 1 and MAP 2 have the same integral coefficient of variation c v ar ¼ 12:3. 0:43680 0:40768 PH service process is characterized by the vector b ¼ ð0:2; 0:8Þ and the matrix S ¼ . 0:42608 1:71809 The mean service time b1 is equal to 1.934, coefficient of variation is equal to 1.936. This PH process was modeled using statistical analysis of service time intervals (including talk time and after call work time) at the call center of the real bank. We assume that the buffer capacity R ¼ 10, the intensities a1 ¼ 1; a2 ¼ 0:5, the probabilities qk ¼ pk ¼ 0:02ðk þ 1Þ; k 2 f0; 1; . . . ; 9g. Let us vary the number of servers N in the interval ½1; 50. Figs. 1–3 illustrate the dependence of the average number of customers in the system e L, the average number of busy servers N serv er , the average number of customers N buffer in the buffer, the average number of priority and non-priority customers denoted by N buffer and N buffer , respectively, the intensity of output flow of customers kout , the intensity of output flow of pri1 2 out ority and non-priority customers denoted by kout 1 and k2 , respectively, on the number of servers for different MAPs presented above. From these figures, one can conclude that coefficients of correlation and variation in the arrival process have the profound impact on the system performance measures. The average number of customers in the system e L, average number of busy servers N serv er and the intensity of output flow kout decrease with increasing of correlation in the arrival process. So, the

Fig. 1. The average number of customers in the system e L and the average number of busy servers N serv er as functions of the number of servers for different MAPs.

Fig. 2. The average number of customers N buffer , priority customers N buffer and non-priority customers N buffer in the buffer as functions of the number of 1 2 servers for different MAPs.

970

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

out Fig. 3. The intensity of flow of customers kout , priority customers kout 1 and non-priority customers k2 which get service as functions of the number of servers for different MAPs.

increase of the coefficient of correlation in the arrival process has negative impact on the intensity of the flow of successfully served customers, which is one of the most important performance measures of the system. Note that besides the presented results we considered analogous dependencies for some MAP arrival processes with negative coefficient of correlation. We can conclude that the values of all characteristics in this case are very close to the ones in case of MAP0 . Note, that on Figs. 1–3 the curves, which characterize some integral measures of the system (without differentiating the types of customers) and correspond to the MAP0:4 and MAP0:4 1 2 , practically coincide. However the curves which separately characterize processing of two types of customers are quite different for these MAPs. The dependence of the average number of customers in the buffer N buffer on the coefficient of correlation in the arrival flow can be explained as follows. For small N, the average number of customers in the buffer decreases with increasing of the coefficient of correlation due to higher loss probability for flows with large coefficient of correlation. When the number of servers N increases, the average number of customers in the buffer is larger for arrival processes with a large coefficient of correlation, because the probability to meet a free server at the arrival epoch of the customer from arrival process with low coefficient of correlation is larger than for the customer from arrival process with large coefficient of correlation. The behavior of the intensity of output flow of priority customers kout 1 is not so predictable. It is surprising that for small N 0:4 the intensity kout can be even greater for the correlated flow (MAP ) than for the stationary Poisson arrival process (MAP 0 ). 2 1 0:4 Explanation of this phenomenon is the following. The arrival process MAP 2 has large coefficient of variation of inter-arrival times of the non-priority customers and due to that the loss probability of non-priority customers from MAP0:4 2 is higher than the same probability for other considered flows. So, high loss probability of non-priority customers makes better conditions for priority customers. Due to that the intensity of output flow of non-priority customers kout is less for MAP0:4 than for 2 2 0:4 MAP1 . For large N, the loss probability of a priority customer from arrival process with low coefficient of correlation tends to zero, and the intensity kout 1 tends to the arrival intensity of priority customers k1 . At the same time, for large N the loss probability of a customer from arrival process with large coefficient of correlation does not tend to zero. So, for large N the intensity kout is greater for the arrival processes with low coefficient of correlation. 1 The dependence of the average number of priority customers N buffer and the average number of non-priority customers 1 N buffer in the buffer on the number of servers can be interpreted by analogy with dependence of N buffer on the number of serv2 ers. But due to the large coefficient of variation of the non-priority customers in MAP0:4 2 for small N average number of pri0:4 ority customers N buffer can be greater for MAP than for the stationary Poisson arrival process (MAP0 ). 2 1 entloss entloss entloss loss loss loss The dependence of the loss probabilities P ; P1 ; P2 ; P ; P 1 and P2 on the number of servers for different MAPs is illustrated on Figs. 4 and 5. One can see from these figures that the loss probabilities Ploss and Pentloss increase with increasing of the correlation in the 0:4 arrival process. Loss probability of priority customers for processes MAP0:4 1 and MAP 2 , which have the same integral coefficients of correlation and variation, but different coefficients of correlation and variation in the flow of priority and non-pri-

Fig. 4. Loss probabilities P entloss ; P entloss and P entloss as functions of the number of servers for different MAPs. 1 2

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

971

Fig. 5. Loss probabilities P loss ; P loss and P loss as functions of the number of servers for different MAPs. 1 2

ority customers, are substantially different due to mentioned above reason. Because of the large coefficient of variation of intervals between arrivals of non-priority customers, the loss probabilities P loss and P2entloss for the process MAP0:4 2 is higher 2 0:4 than for MAP1 . As the loss probability of non-priority customer is greater, the loss probabilities of priority customers P loss 1 0:4 and P 1entloss for MAP0:4 2 is lower than for MAP 1 . imploss imploss imploss The dependence of the loss probabilities P ; P1 and P2 on the number of servers for different MAPs is illustrated on Fig. 6. One can see that the dependencies of the probabilities Pimploss ; P imploss and P imploss on the number of servers are very sim1 2 buffer ilar to the same dependencies of N buffer , N buffer and N correspondingly. It can be explained by the formula 1 2 kr Pimploss ¼ ar N buffer ; r ¼ 1; 2. This formula holds true because the right and left sides of this equation determine the intensity r r of leaving the system due to impatience by type r customer. So, the probability Pimploss is proportional to the average number r of type r customers in the buffer N buffer . r Figs. 7 and 8 show the dependence of the average sojourn time, waiting time of an arbitrary priority and non-priotity customer in the system, and the average sojourn time of an arbitrary serviced priority and non-priotity customer in the system on the number of servers for different MAPs. It is worth to note that the results of this experiment show that the famous Little’s formula holds true for the system under study, i.e.

V wait ¼ r

Nbuffer r ; kr

which also equals to

imploss

Pr

ar

; r ¼ 1; 2.

Fig. 6. Loss probabilities P imploss ; P imploss and P imploss as functions of the number of servers for different MAPs. 1 2

wait v Fig. 7. V soj and V ser as functions of the number of servers for different MAPs. 1 1 ; V1

972

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

wait v Fig. 8. V soj and V ser as functions of the number of servers for different MAPs. 2 2 ; V2

Experiment 2. The problem of optimal choice of the number of CSRs is very important for effective organization of operation of a call center. In this experiment we solve numerically the problem of optimal choice of the number of servers N. The aim of optimization is maximization of one of possible cost criteria (average profit of the system per unit time):

loss loss out JðNÞ ¼ a1 kout 1 þ a2 k2  c 1 k1 P 1  c 2 k2 P 2  dN:

Here ar ; r ¼ 1; 2 is the average profit obtained by the system from the service of one type r customer, kout is the intensity of r the flow of type r customers, which get successful service in the system, cr ; r ¼ 1; 2 is the charge of the system when a type r customer is lost, P loss is the loss probability of a type r customer, d is the charge paid for maintenance of one server per unit r time. Our goal is to find the optimal value N which provides the maximal value of this cost criterion. It is worth to note, that the problem of the suitable choice of the cost coefficients (in our case coefficients ar ; cr ; r ¼ 1; 2) in the cost criterion always play crucial role in successful implementation of optimization. We assume here that in our model the cost coefficients can come from the experts in the real call center to which the model will be applied. So, we fixed the following cost coefficients in the cost criterion: a1 ¼ 20; a2 ¼ 15; c1 ¼ 10; c2 ¼ 5; d ¼ 3. The rest of parameters are the same as ones presented in the first experiment. The dependence of the cost criterion J on the number N of servers for different MAPs is presented on Fig. 9. Let us analyze the result of optimal design of the number of servers. As it is seen from Fig. 9, the optimal value of cost criterion J is 93.5128 and the optimal number of servers N ¼ 25 for MAP0 , the optimal value of cost criterion J is 77:1070 and the optimal number of servers N ¼ 30 for MAP0:2 , in case MAP0:4 J ¼ 12:1619 under N ¼ 37, and in case 1 MAP0:4 J ¼ 10:5863 under N ¼ 38. 2 As a conclusion to this experiment we should say that the optimal number of servers N is very sensitive to the coefficient of correlation in the arrival process. Using the models with stationary Poisson arrival process for performance evaluation of the call centers with correlated arrival flows of customers can lead to huge error in evaluation of the profit. For example, if we assume that the arrival flow of priority and non-priority customers to the call center is described by the stationary Pois son arrival process, while actually the correlated flow MAP0:4 2 arrives to the call center, we will hire N ¼ 25 call center agents and expect the profit J ¼ 77:1070. But actually if we will hire N ¼ 25 agents, we will get zero profit. So, prediction of the system performance under assumption that the customers arrive according to the stationary Poisson process, which is used as the descriptor of arrival processes in call centers literature, can be essentially wrong and too optimistic. Experiment 3. This experiment has two goals. The one is to show that the buffer capacity R is an important parameter of the system and effectiveness of the system operation can be improved by means of the proper choice of R. The second goal is to illustrate that characteristics of the system depend not only on the coefficient of correlation in the arrival process, but also on variation in the PH service process.

Fig. 9. The cost criterion as a function of the number of servers for different MAPs.

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

973

Fig. 10. The cost criterion as a function of the buffer capacity for different variation in the service process.

So, we formulate and numerically solve the following problem of optimal choice of the buffer capacity. The aim of optimization is maximization of the following cost criterion loss loss out IðRÞ ¼ a1 kout 1 þ a2 k2  c 1 k1 P 1  c 2 k2 P 2  bR

which determines the average profit of the system per unit time. Here b is the charge paid for maintenance of one place in the buffer per unit time. Our goal is to find the optimal value R which provides the maximum value to the cost criterion IðRÞ. To investigate the impact of the variation in the service process we consider three PH service processes with the same average service time b1 ¼ 2:98, but different coefficients of variation. The first PH process is characterized by the vector   0:669798 0:669798 b ¼ ð1; 0Þ and the matrix S ¼ . It corresponds to the Erlang distribution. 0 0:669798 The second PH process is characterized by the vector b ¼ ð1Þ and the matrix S ¼ ð0:334899Þ. It corresponds to the expo  0:1 0:01 nential distribution. The third PH process is characterized by the vector b ¼ ð0:2; 0:8Þ and the matrix S ¼ . 0:02 1 The squared coefficients of variation of these service processes are cv ar ¼ 0:5; cv ar ¼ 1, and cv ar ¼ 4, correspondingly. We consider the MAP0:4 presented in the first experiment. Let us fix the following parameters of the system: 2 N ¼ 10; qk ¼ pk ¼ 0:05; k ¼ 0; R  1. The rest of parameters are the same as ones presented in the first experiment. Let us vary the buffer capacity R in the interval [1,20]. The cost coefficients in the cost criterion are assumed to be as follows: a1 ¼ 40; a2 ¼ 20; c1 ¼ 10; c2 ¼ 5; b ¼ 0:5. The dependence of the cost criterion I on the buffer capacity R for different PH service processes is presented on Fig. 10. One can conclude from Figure 10 that if the coefficient of variation increases than the optimal buffer capacity R decreases. The optimal value of cost criterion I is equal to 14:71 and the optimal buffer capacity R ¼ 11 for Erlang service time distribution with coefficient of variation cv ar ¼ 0:5, the optimal value of cost criterion I is 15.94 and the optimal buffer capacity R ¼ 10 for exponential service process which has coefficient of variation cv ar ¼ 1, and for PH service process with coefficient of variation cv ar ¼ 4 the optimal value I ¼ 19:49 under the optimal buffer capacity R ¼ 9. As a conclusion from this experiment we can say that the coefficient of variation of the service time distribution impacts on the system performance measures, but not so essential as the coefficient of correlation in the arrival process. One can see from this numerical example that for correct prediction of the system operation it is necessary to take into account the variation of the service time distribution. 8. Conclusion In this paper we consider a novel multi-server queueing system with impatient heterogeneous customers which can be suited for modeling of a real-life call center. The special approach for reducing the system dimension is used. The efficient algorithm for calculating the stationary probabilities of system states is presented. The main performance measures and the Laplace–Stieltjes transforms of the sojourn and waiting time distribution are derived. The numerical results confirm the importance of the analysis of the queueing models with the MAP arrival process because correlation in the arrival process essentially impacts on the system operation. Importance of taking into account the coefficient of variation of the service time distribution is also illustrated numerically. Possibility of effective use of the presented results for solving the problems of optimal design of modern call centers is demonstrated. Acknowledgements This research was supported by Basic Science ResearchProgram through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant No. 2011-0015214).

974

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

Appendix A. Computation process A.1. Calculation of the matrices LNi ðN; e SÞ and Ai ðN; SÞ; i 2 f0; . . . ; Ng 1. Calculate the matrices sðkÞ ðSÞ; k 2 f0; . . . ; M  1g, which are obtained by removing of k first rows and k first columns from the matrix S. 2. Calculate the matrices T j ¼ sðM2jÞ ðSÞ; j 2 f1; . . . ; M  2g. ðwÞ 3. Calculate the matrices Li ðT j Þ using the recurrent formulas: ð0Þ

Li ðT j Þ ¼ ðN  iÞtjrj ;1 ; 0 B B B B B B ðwÞ Li ðT j Þ ¼ B B B B B B @

i 2 f0; . . . ; N  1g;

ðN  iÞt jrj w;1 I

j 2 f1; . . . ; M  2g;

O



LN1 ðT j Þ

ðN  i  1Þt jrj w;1 I



O .. .

LN2 ðT j Þ ...

ðw1Þ

 .. .

O

O



O

O

   Li

ðw1Þ

w 2 f1; . . . ; rj  2g;

C C C C C C O C; C .. C . C C j trj w;1 I C A O

ðw1Þ

i 2 f0; . . . ; N  1g;

1

O

ðT j Þ

j 2 f1; . . . ; M  2g;

t jk;l

where is the ðk; lÞ th entry of the matrix T j and rj is the number of rows of the matrix T j . ðwÞ 4. Calculate the matrices U i ðT j Þ using the recurrent formulas: ð0Þ

U i ðT j Þ ¼ t j1;rj ; 0 ðwÞ U i ðT j Þ

B B B ¼B B B @

i 2 f1; . . . ; Ng;

tj1;rj w I

ðw1Þ

UN

O



O

O

U N1 ðT j Þ   

O

O

ðT j Þ

O

t j1;rj w I

.. .

.. .

O

j 2 f1; . . . ; M  2g;

ðw1Þ

.. .

O

w 2 f1; . . . ; rj  2g;

..

O

i 2 f1; . . . ; Ng;

.. .

.. .

tj1;rj w I

ðw1Þ Ui ðT j Þ

.



1 C C C C; C C A

j 2 f1; . . . ; M  2g:

ðr 2Þ Li j ðT j Þ;

ðr j 2Þ

5. Calculate the matrices Li ðN; T j Þ ¼ i 2 f0; . . . ; N  1g, and U i ðN; T j Þ ¼ iU i ðwÞ 6. Calculate the matrices Ai using the recurrent formulas:

0

ð0Þ

Ai

0

B SM;M1 B B B 0 B ¼B . B .. B B @ 0 0

iSM1;M

0



0

0 2SM;M1 .. .

ði  1ÞSM1;M 0 .. .

  .. .

0 0 .. .

0 0

0 0

0

ðjÞ

Ai

i U ðN; T j Þ O N N B ðj1Þ B LN1 ðN; T j Þ A1 B B B O LN2 ðN; T j Þ B ¼B .. .. B . . B B B O O @

O i 2 f1; . . . ; Ng;

 0    iSM;M1 O

i1 U ðN; T j Þ N1 N1 ðj1Þ A2

O

.. . O O

0

1

C C C C C C; C C C SM1;M A 0 0 .. .

ðT j Þ; i 2 f1; . . . ; Ng; j 2 f1; . . . ;M  2g.

i 2 f1; . . . ; Ng;

0



O

O



O

O

1

C C C C C  O O C C; .. .. .. C . . . C C ðj1Þ 1 U ðN; T j Þ C  Ai1 A Niþ1 Niþ1 ðj1Þ    LNi ðN; T j Þ Ai

j 2 f1; . . . ; M  2g: ðM2Þ

7. Calculate the matrices Ai ðN; SÞ as follows A0 ðN; SÞ ¼ O11 ; Ai ðN; SÞ ¼ Ai ; i 2 f1; . . . ; Ng. ðM1Þ e 8. Calculate the matrices Li ðN; e SÞ as follows Li ðN; e SÞ ¼ Li ð SÞ; i 2 f0; . . . ; N  1g; LN ðN; e SÞ ¼ O11 .

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

975

A.2. Calculation of the matrices P i ðbÞ; i 2 f0; . . . ; N  1g ðjÞ

1. Calculate the matrices P i of size ði þ 1Þ  ði þ 2Þ using the recurrent formulas:

0 ð0Þ

Pi

bM1

B 0 B ¼B B .. @ . 0 0

ðjÞ

Pi

B B B B ¼B B B B @

bM

0



0

0

1

bM1 .. .

bM .. .

 .. .

0 .. .

0 .. .

C C C; i 2 f1; . . . ; N  1g; C A

0

0

   bM1

bMj1

zðjÞ

0T

bMj1 I

0T .. .

O .. .

bMj1 I .. .

0T

O

O

j 2 f1; . . . ; M  2g;

0 ðj1Þ

P1

bM 0



0

O



O

ðj1Þ

.. .

 .. .

O .. .

O

   bMj1 I

P2

0

1

C O C C C O C C; .. C C . A

ðj1Þ

Pi

i 2 f1; . . . ; N  1g;

ðjÞ

where the vectors z ¼ ðbMj ; bMjþ1 ; . . . ; bM Þ; j 2 f1; . . . ; M  2g. ðM2Þ 2. Calculate the matrices P i ðbÞ as follows P 0 ðbÞ ¼ b; Pi ðbÞ ¼ P i ; i 2 f1; . . . ; N  1g. A.3. Algorithm for calculation of the vector

v 2 ðsÞ

1. Calculate the matrix functions Z i ðsÞ; i 2 f2; . . . ; Rg, using the recursion

Z i ðsÞ ¼ ðXi1;i1  sI þ Xi1;i2 Z i1 ðsÞÞ1 Xi1;i ;

i 2 f3; . . . ; Rg;

1

under the initial condition Z 2 ðsÞ ¼ ðX1;1  sIÞ X1;2 . 2. Calculate the vectors ci ðsÞ; i 2 f2; . . . ; R þ 1g, as follows

c2 ðsÞ ¼ ðX1;1  sIÞ1 a1 ðsÞ; ci ðsÞ ¼ ðXi1;i1  sI þ Xi1;i2 Z i1 ðsÞÞ1 ðai1 ðsÞ  Xi1;i2 ci1 ðsÞÞ; i 2 f3; . . . ; R þ 1g: 3. The components

v2 ðs; iÞ;

i 2 f1; . . . ; Rg, of the vector

v2 ðsÞ are calculated as

v2 ðs; RÞ ¼ cRþ1 ðsÞ; v 2 ðs; iÞ ¼ Ziþ1 ðsÞv 2 ðs; i þ 1Þ þ ciþ1 ðsÞ;

i ¼ R  1; R  2; . . . ; 1:

References [1] K.A. Gilson, D.K. Khandelwal, Getting more from call centers, McKinsey Quarterly: The Journal of McKinsey and Co, 2005. [2] O.Z. Aksin, M. Armony, V. Mehrotra, The modern call centers: a multi-disciplinary perspective on operations management research, Working paper, Koc University, Istanbul, Turkey, 2007. [3] O. Garnett, A. Mandelbaum, M. Reiman, Designing a call center with impatient customers, Manufact. Service Oper. Manag. 4 (2002) 208–227. [4] M.S. Aguir, O.Z. Aksin, F. Karaesmen, Y. Dallery, On the interaction between retrials and sizing of call centers, Euro. J. Oper. Res. 191 (2008) 398–408. [5] O. Jouini, Y. Dallery, Z. Aksin, Queuing models for flexible multi-class call centers with real-time anticipated delays, Int. J. Product. Econom. 120 (2009) 389–399. [6] J.W. Kim, S.C. Park, Outsourcing strategy in two-stage call centers, Comput. Oper. Res. 37 (2010) 790–805. [7] T. Aktekin, R. Soyer, Call center arrival modelling: a Bayesian state-space approach, Naval Res. Log. 58 (2011) 28–42. [8] A. Pattavina, A. Parini, Modelling voice call inter-arrival and holding time distributions in mobile networks. Performance challenges for efficient next generation networks in: Proceedings of 19th International Teletraffic Congress, 2005, pp. 729–738. [9] A. Riska, V. Diev, E. Smirni, Efficient fitting of long-tailed data sets into hyperexponential distributions, in: Global Telecommunications Conference (GLOBALCOM’02, IEEE), 2002, pp. 2513–2517. [10] V. Ramaswami, Independent Markov processes in parallel, Commun. Statist. Stochast. Models 1 (1985) 19–432. [11] V. Ramaswami, D.M. Lucantoni, Algorithms for the multi-server queue with phase-type service, Commun. Statist. Stochast. Models 1 (1985) 393–417. [12] S.R. Chakravarthy, The batch Markovian arrival process: a review and future work, in: A. Krishnamoorthy, et al.(Eds), Advances in Probability Theory and Stochastic Processes, Notable Publications, NJ, 2001, pp. 21–49. [13] D. Lucantoni, New results on the single server queue with a batch Markovian arrival process, Commun. Statist. Stochast. Models 7 (1991) 1–46. [14] T. Ryden, Statistical inference for Markov-modulated Poisson processes and Markovian arrival processes, Third International Conference on Matrix Analytic Methods in Stochastic Models, Leuven, July 2000 (plenary talk), in: G. Latouche and P. Taylor (Eds.), Advances in Algorithmic Methods for Stochastic Models, Notable Publications, Neshanic Station, NJ, pp. 329–350. [15] Q.M. He, Queues with marked customers, Adv. Appl. Prob. 28 (1996) 567–587. [16] C.S. Kim, S.A. Dudin, Priority tandem queueing model with admission control, Comput. Indust. Eng. 61 (2011) 131–140. [17] B. Cleveland, Call center management on fast forward: succeeding in today’s dynamic customer contact environment, ICMI, 2006. [18] M. Neuts, Matrix-geometric solutions in stochastic models – an algorithmic approach, John Hopkins University Press, 1981. [19] A. Panchenko, A. Thummler, Efficient phase-type fitting with aggregated traffic traces, Perform. Eval. 64 (2007) 629–645.

976

C. Kim et al. / Applied Mathematical Modelling 37 (2013) 958–976

[20] A. Graham, Kronecker Products and Matrix Calculus with Applications, Ellis Horwood, Cichester, 1981. [21] H. Kesten, J.Th. Runnenburg, Priority in waiting line problems, Amsterdam, Mathematisch Centrum, 1956. [22] D. van Danzig, Chaines de Markof dans les ensembles abstraits et applications aux processus avec regions absorbantes et au probleme des boucles, Ann. de l’Inst. H. Pioncare. 14 (1995) 145–199. [23] J. Abate, W. Whitt, Numerical inversion of the Laplace transforms of probability distributions, ORSA J. Comput. 7 (1995) 36–43. [24] S.D. Conte, C. deBoor, Elementary Numerical Analysis – An Algorithmic Approach, McGraw-Hill, New York, 1972.