# Queueing System with Heterogeneous Customers as a Model of a Call Center with a Call-Back for Lost Customers

Authors: Sergey Dudin; Chesoong Kim; Olga Dudina; Janghyun Baek
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A multiserver queueing system with infinite and finite buffers, two types of customers,
and two types of servers as a model of a call center with a call-back for lost
customers is investigated. Type 1 customers arrive to the system according to a
Markovian arrival process. All rejected type 1 customers become type 2 customers.
Typer, r=1,2, servers serve type r customers if there are any in the system and
serve type r′, r′=1,2,  r′≠r, customers if there are no type r customers in the system.
The service times of different types of customers have an exponential distribution
with different parameters. The steady-state distribution of the system is analyzed.
Some key performance measures are calculated. The Laplace-Stieltjes transform of
the sojourn time distribution of type 2 customers is derived. The problem of optimal
choice of the number of each type servers is solved numerically.


---

## Body

## 1. Introduction

A call center is a specialized unit of companies that handles voice requests from clients. As call centers provide a communication channel between companies and their customers, the effective call center performance is extremely important for company’s reputation. For modeling and optimizing call centers, the queueing theory is used. For background information and an overview of the present state of the art in the study of call centers, the reader is referred to the survey [1] and the papers [2, 3], as well as references therein.It is not a secret that there is a high competition between companies, and a customer, who does not receive service in a call center of one company, can simply move to competitors’ company. Loss of customers or potential customers may cause damage or loss of profit. This is especially important for such organizations as banks, Internet shops, telemarkets, and other organizations that deal with selling and providing the paid services. For these organizations, loss of calls causes loss of clients and increases the risk of clients dealing with competitors. For example, if a potential client, who wants to get the information about deposits, does not succeed in connecting with an operator of a call center of one bank, this client can contact another bank and open a deposit there. Thus, the issue of effective service of calls with minimum losses of customers is very important. At the same time, to achieve the call center performance without losses is very difficult because(i)
the arrival flow of customers can be correlated (the time intervals in which customers arrive in the system often alternate with intervals in which few customers arrive in the system);(ii)
it is too expensive to have a large number of operators (most of the costs of a call center involve wages).So, it seems to be reasonable to call back a customer, who did not wait for the answer of an operator of a call center, and offer him (her) the services. We assume that if customer leaves the system without service, his (her) phone number is added to some infinite size virtual queue, and an operator calls back him (her) later. Let us refer to a customer who leaves the system without service as a virtual customer and a customer who waits on hold as a real customer. One way is to provide a call-back to virtual customers when there is an idle operator and there are no customers waiting on hold. But in this case, the waiting time of a call-back can be very long and the customer will no longer need the services provided by the call center. In our paper, we consider more reasonable and general way to provide a call-back. We assume that all call center operators are divided into two groups. For the first group of operators (type 1), processing of real customers has a priority, while the operators from the second group (type 2) firstly provide service to virtual customers, if there are any in the system. This approach allows us to solve the issue of determination of the number of different types of operators under constraints on the waiting time of an arbitrary virtual customer but complicates the investigation of the system.Models of call centers with a call-back option were considered in [4–6]. In all these models, a call-back was considered as an option which is proposed to a customer. We consider the model of a call center in which a call-back is not proposed to the customers but operators who call back all customers leaving the system without service. In [4, 5], an asymptotic analysis of the model of a call center with a call-back option and identical servers in the case of heavy load for a large number of operators is presented. In our paper, we consider more general customers’ arrival process and the possibility of customers abandonment. Also, we present an exact analytical analysis of the model without the restriction that the number of operators is sufficiently large. In [6], all operators are assumed to be identical and call-backs to virtual customers are provided only when there is an idle operator and there are no real customers in the system. In our paper, we considered a more general and complicated model where all operators are divided into groups.The paper is organized as follows. In Section2, the mathematical model is described. The ergodicity condition and the stationary distribution of the system states are analyzed in Section 3. The expressions for the main performance measures of the system are given in Section 4. The Laplace-Stieltjes transform of the sojourn time distribution of an arbitrary virtual customer is presented in Section 5. Section 6 contains numerical illustrations and their short discussion. Section 7 concludes the paper.


## 2. Mathematical Model

We consider anN-server queuing system with two types of customers, a buffer of capacity K for real customers, and an infinite buffer for virtual customers to model a call center with a call-back option. The structure of the system under study is presented in Figure 1.Figure 1
Structure of the system.Real (type 1) customers arrive at the system according to a Markovian arrival process (MAP). The MAP is defined by the underlying processνt, t≥0, which is an irreducible continuous-time Markov chain with the state space {0,1,…,W}. Arrivals occur only at the epochs of transitions in the underlying process νt, t≥0. The intensities of transitions of the process νt, t≥0, that are accompanied (not accompanied) by the arrival of a customer are defined by the square matrix D1(D0) of size W-=W+1.The matrixD(1)=D0+D1 is an infinitesimal generator of the process νt, t≥0. The stationary distribution vector θ of this process satisfies the system of equations θD(1)=0, θe=1. Here and throughout this paper, 0 is a zero row vector, and e denotes a unit column vector.The average intensityλ (fundamental rate) of the MAP is defined by λ=θD1e. The coefficient of variation cvar of intervals between customer arrivals is calculated as cvar2=2λθ(-D0)-1e-1. The coefficient of correlation ccor of intervals between successive arrivals is given as ccor=(λθ(-D0)-1(D(1)-D0)(-D0)-1e-1)/cvar2.If there is an available server during an arbitrary real customer arrival epoch, this customer is admitted into the system and occupies a free type 1 server if there is, or a type 2 server, otherwise. If all servers are busy during a real customer arrival epoch, this customer becomes aware of the current queue length (“visible” queue) and, based on the information provided, decides whether to wait in line or leave the system (balk). Thus, if all servers are busy and there arek, k=0,K-1¯, real customers in the finite buffer during an arbitrary real customer arrival epoch, the arriving customer leaves the system with probability pk, and with the complementary probability the customer chooses not to wait for service. If the buffer for type l customers is full during an arbitrary customer arrival epoch, this customer always leaves the system; that is, pK=1.Real customers are impatient; that is, a real customer leaves the buffer after an exponentially distributed time described by the parameterα, 0<α<∞, due to a lack of service.We assume that, if a real customer leaves the system without service, this customer changes its type (becomes a virtual customer) and joins a virtual waiting place of infinite size (infinite buffer). Virtual customers are always patient. We assume that both types of customers are served according to the First In-First Out rule. The difference between two types of customers is that a real customer physically presents in the system and holds a line, while a virtual customer is not presented in the system.All servers (operators) are divided into two types:(1)
type 1 operators who serve the real customers;(2)
type 2 operators who call back the customers who left the system (virtual customers).The number of type 1 servers isM, M=0,N¯, and the number of type 2 servers is N-M.If a type 1 operator becomes idle and there are no real customers in the queue but there are virtual customers who wait for a call-back, this operator starts to serve a virtual customer to avoid the idle time. Similarly, if a type 2 operator becomes idle and there are no virtual customers in the queue but there are real customers waiting for service, then this operator starts to serve a real customer.The service times of real and virtual customers have an exponential distribution with the parametersμ1 and μ2, respectively.


## 3. The Process of System States

Letit be the number of customers in the system, kt the number of real customers in the finite buffer, mt the number of busy type 1 servers, nt the number of busy type 1 servers that serve the real customers, lt the number of busy type 2 servers that serve the real customers, and νt the state of the underlying process of the MAP during the epoch t, t≥0.The behavior of the system under consideration can be described in terms of the regular irreducible continuous-time Markov chainξt={it,kt,mt,nt,lt,νt}, t≥0, with the state space
(1)(0,min{i-m,N-M}¯{i,m,n,l,ν},i<N,m=max{0,i-(N-M)},min{i,M}¯,n=0,m¯,l=0,min{i-m,N-M}¯,ν=0,W¯)⋃({i,k,n,l,ν},i≥N,k=0,min{i-N,K}¯,n=0,M¯,l=0,N-M¯,ν=0,W¯).Let us introduce the operator𝒯(Ak,k-1,k=k1,K1¯;Ak,k,k=k2,K2¯;Ak,k+1,k=k3,K3¯) which define a block tridiagonal matrix with subdiagonal blocks Ak,k-1, diagonal blocks Ak,k, and updiagonal blocks Ak,k+1. By default, we assume that this block matrix is square or the number of its block rows is equal to the number of block columns ±1.Let us enumerate states of Markov chainξt in the lexicographic order of the components k, m, n, l, and ν. Let Q be a generator of the Markov chain ξt,t≥0, consisting of the blocks Qi,j, which define the transition rates of this chain from the states having value i of the first component to the states having value j of this component of the chain. The diagonal entries of the matrices Qi,i are negative, and the modulus of diagonal entries defines the total intensity of leaving the corresponding state of the Markov chain ξt,t≥0.Lemma 1.
The infinitesimal generatorQ of the Markov chain ξt,t≥0, has a block tridiagonal structure:(2)Q=(Q0,0Q0,1OO⋯OOOOO⋯Q1,0Q1,1Q1,2O⋯OOOOO⋯OQ2,1Q2,2Q2,3⋯OOOOO⋯⋮⋮⋱⋱⋱⋮⋮⋮⋮⋮⋯OOOO⋯QN+K-1,N+K-2QN+K-1,N+K-1QN+K-1,N+KOO⋯OOOO⋯OQN+K,N+K-1Q1Q2O⋯OOOO⋯OOQ0Q1Q2⋯⋮⋯⋮⋮⋮⋮⋮⋮⋱⋱⋱).
The nonzero blocksQi,j, i,j≥0, have the following form:
(3)Qi,i=diag{i*,min{i,M}¯diag{diag{0,min{N-M,i-m}¯A(i,m,n,l),l=0,min{N-M,i-m}¯},n=0,m¯},m=i*,min{i,M}¯},i=0,N-1¯,Qi,i=𝒯(0,i-N-1¯Bk,k-1,k=1,i-N¯;Bk,k,k=0,i-N¯;Ok,k=0,i-N-1¯),i=N,N+K-1¯,Q1=𝒯(Bk,k-1,k=1,K¯;Bk,k,k=0,K¯;Ok,k=0,K-1¯),i≥N+K,Qi,i+1={𝒯(Om,m=1,i¯;Om,m(i),m=0,i¯;Cm,m+1(i),m=0,i¯),i<N-M;diag{Cm,m+1(i),m=i*,i¯},i≥N-M,i=0,M-1¯,Qi,i+1={𝒯(Om,m=1,M¯;Cm,m(i),m=0,M¯;Cm,m+1(i),m=0,M-1¯),i<N-M;𝒯(Cm,m(i),m=i*+1,M¯;Cm,m+1(i),m=i*,M-1¯;Om,m=i*,M-2¯),i≥N-M,i=M,N-1¯;Qi,i+1=𝒯(Ok,k=1,i-N¯;Hk,k,k=0,i-N¯;Hk,k+1,k=0,i-N¯),i=N,N+K-1¯,Q2=𝒯(Ok,k=1,K¯;Hk,k,k=0,K¯;Hk,k+1,k=0,K-1¯),i≥N+K,Qi,i-1={𝒯(Em,m-1(i),m=1,i¯;Em,m(i),m=0,i-1¯;Om,m=0,i-2¯),i≤N-M;𝒯(Om,m=i*+1,i¯;Em,m-1(i),m=i*,i¯;Em,m(i),m=i*,i-1¯),i>N-M,i=1,M¯,Qi,i-1={𝒯(Em,m-1(i),m=1,M¯;Em,m(i),m=0,M¯;Om,m=0,M-1¯),i≤N-M;𝒯(i*,M¯Om,m=i*+1,M¯;Em,m-1(i),m=i*,M¯;Em,m(i),m=i*,M¯),i>N-M,i=M+1,N¯,Qi,i-1=𝒯(Zk,k-1(i),k=1,i-N¯;Zk,k,k=0,i-N-1¯;Ok,k=0,i-N-2¯),i=N+1,N+K¯,Q0=𝒯(Vk,k-1,k=1,K¯;Vk,k,k=0,K¯;Ok,k=0,K-1¯),i>N+K.
Here,I is the identity matrix, O is a zero matrix of appropriate dimension, and Ol, l=k,K¯, is the sequence which consists of K-k+1 zero matrices of appropriate dimension:
(4)i*=max{0,i-(N-M)}.⊗
indicates the Kronecker product of matrices:(5)A(i,m,n,l)={-[(n+l)μ1+(i-n-l)μ2]IW-+D0,i-m≤N-M,D0,i-m>N-M;Bk,k=-B⊗IW-+I(M+1)(N-M+1)⊗(-kαIW-+D0) where B=diag{diag{(n+l)μ1+(N-n-l)μ2,l=0,N-M¯},n=0,M¯};
(6)Bk,k-1=kαI(M+1)(N-M+1)W-;Cm,m+1(i)=𝒯(On,n=1,m¯;On,n=0,m¯;Cn,n+1(i,m),n=0,m¯) where Cn,n+1(i,m)=Imin{i-m,N-M}+1⊗D1;(7)Cm,m(i)={O(i+1)W-×(i+2)W-,m=0,i<N-M;O(N-M+1)W-×(N-M+1)W-,m=0,i≥N-M;O(m+1)(min{i-m,N-M}+1)W-×(m+1)(min{i-m+1,N-M}+1)W-,m=1,M-1¯;IM+1⊗Ii+⊗D1,m=M,where Ii+ is a matrix of size (i-M+1)×(i-M+2) with nonzero entries (Ii+)l,l+1,l=0,i-M¯, which are equal to 1. Consider(8)Om,m(i)={O(i+1)W-×(i+2)W-,m=0,O(m+1)(min{i-m,N-M}+1)W-×(m+1)(min{i-m+1,N-M}+1)W,-m=1,i¯;Hk,k=I(M+1)(N-M+1)⊗pkD1;Hk,k+1=I(M+1)(N-M+1)⊗(1-pk)D1;Em,m(i)=Im+1⊗𝒯(μ1l,l=1,i-m¯;μ2(i-m-l),l=0,i-m-1¯;Ol,l=0,i-m-2¯)⊗IW-;Em,m-1(i)=𝒯(μ1l,l=1,m¯;μ2(m-l),l=0,m-1¯;Ol,l=0,m-2¯)⊗I(i-m+1)W-;Zk,k={𝒯(μ1lIW-(N-M+1),l=1,M¯;Ψ+μ2(M-l)I(N-M+1)W-,l=0,M¯;dddOl,l=0,M-1¯),k=0;IM+1⊗Ψ,k=1,K-1¯,Ψ=𝒯(μ1l,l=1,N-M¯;μ2(N-M-l),l=0,N-M¯;Ol,l=0,N-M-1¯)⊗IW-;Zk,k-1(i)={Φ,k=1,i-N-1¯;Φ+IM+1⊗𝒯(Ol,l=1,N-M¯;μ1l,l=0,N-M¯;dddddddddddddμ2(N-M-l),l=0,N-M-1¯)⊗IW-,k=i-N,Φ=𝒯(On,n=1,M¯;μ1n,n=0,M¯;μ2(M-n),n=0,M-1¯)⊗I(N-M+1)W-;Vk,k={𝒯(μ1nIN-M+1,n=1,M¯;V-n,n,n=0,M¯;On,n=0,M-1¯)⊗IW-,k=0,IM+1⊗Ψ,k=1,K¯,where V-n,n=𝒯(μ1l,l=1,N-M¯;μ2(N-n-l),l=0,N-M¯;Ol,l=0,N-M¯). Consider
(9)Vk,k-1=𝒯(On,n=1,M¯;μ1n,n=0,M¯;μ2(M-n),n=0,M-1¯)⊗I(N-M+1)W-.The proof of the lemma is implemented by the analysis of all transitions of the Markov chainξt,t≥0, during the interval of infinitesimal length and rewriting the intensities of these transitions in block matrix form.Corollary 2.
The Markov chainξt, t≥0, belongs to the class of continuous time quasi-birth-and-death processes; see, for example, [7].Theorem 3.
The Markov chainξt,t≥0, is ergodic if the following inequality holds true:
(10)y𝒜(B⊗IW-)e>y𝒜(IM+1⊗D1)e,
where the vector y is the unique solution to the system
(11)y(Q-0,0+T0Q-1,0)=0,y(I+𝒜)e=1.
Here,(12)𝒜=∑k=0K‍∏l=0k‍Tl;
the matrices Tk,k=0,K-1¯, are calculated using the backward recursion
(13)TK-1=-Q-K-1,K(Q-K,K)-1,(14)Tk=-Q-k,k+1(Q-k+1,k+1+Tk+1Q-k+2,k+1)-1,k=K-2,K-3,…,0;
the matrices Q-k,k, Q-k,k+1, and Q-k,k-1 are defined as follows:
(15)Q-0,0=𝒯(μ1n,n=1,M¯;-μ1n,n=0,M¯;On,n=0,M-1¯)⊗IW-+IM+1⊗(D0+p0D1),Q-k,k=diag{-μ1n-μ2(M-n),n=0,M¯}⊗IW--kαI(M+1)W-+IM+1⊗(D0+pkD1),k=1,K¯,Q-k,k+1=(1-pk)(IM+1⊗D1),k=0,K-1¯,Q-k,k-1=kαI(M+1)W-+𝒯(On,n=1,M¯;μ1n,n=0,M¯;(M-n)μ2,n=0,M-1¯)⊗IW-,k=1,K¯.Proof.
It follows from [7] that the necessary and sufficient ergodicity condition of the quasi-birth-and-death process is the fulfillment of the inequality
(16)xQ0e>xQ2e,
where the vector x is the unique solution to the following system of linear algebraic equations:
(17)x(Q0+Q1+Q2)=0,xe=1.
Let the row vectorx have the form
(18)x=(x0,x1,…,xK),
where xk=x^k⊗e^⊗θ, k=0,K¯. Here, e^ is a row vector of size N-M+1 with all zero entries except the entry (e^)0 which is equal to 1, θ is the stationary distribution vector of the process νt,t≥0, and x^k, k=0,K¯ is a row vector of size M+1 such as the vector (x^0,…,x^K) is stochastic.
By direct substitution, we verify that the vectorx from (18) provides the unique solution of system (17) if the vector x-=(x-0,x-1,…,x-K), where x-k=x^k⊗θ is the unique solution to the following system of linear algebraic equations:
(19)x-Q-=x-,x-e=1.
Here,Q- is a block tridiagonal matrix with subdiagonal blocks Q-k,k-1, diagonal blocks Q-k,k, and updiagonal blocks Q-k,k+1 defined by formulas (15).
To calculate the vectorsx-k, k=0,K¯, we use the numerically stable algorithm which consists of the following steps.
Step 1. Calculate the matrices Tk, k=0,K-1¯, using the backward recursion (14) under the initial condition (13).
Step 2. Calculate the vector x-0 as the unique solution to the following system:
(20)x-0(Q-0,0+T0Q-1,0)=0,x-0(I+∑k=0K‍∏l=0k‍Tl)e=1.
Step 3. Calculate the vectors x-k=x-0∏l=0kTl, k=1,K¯.
Substituting the vectorx from (18) into inequality (16) and performing some algebra, we obtain the inequality
(21)∑k=0K‍(x^k⊗θ)(B⊗IW-)e>∑k=0K‍(x^k⊗θ)(IM+1⊗D1)e
or
(22)∑k=0K‍x-k(B⊗IW-)e>∑k=0K‍x-k(IM+1⊗D1)e.
Let the matrix𝒜 be given by formulas (12) and denote the vector x-0 as y. It is evident that inequality (22) takes form (10) and system (20) is reduced to system (11).If ergodicity condition (10) of the Markov chain ξt is fulfilled, then the stationary probabilities of system states π(i,m,n,l,ν), i<N, m=max{0,i-(N-M)},min{i,M}¯,n=0,m¯,  l=0,min{i-m,N-M}¯, ν=0,W¯, and π(i,k,n,l,ν), i≥N, k=0,min{i-N,K}¯, n=0,M¯, l=0,N-M¯, ν=0,W¯, exist. Let us form the row vectors πi of these probabilities enumerated in the lexicographic order of the components (m,n,l,ν) if i<N and in the lexicographic order of the components (k,n,l,ν) if i≥N.It is well known that the probability vectorsπi, i≥0, satisfy the following system of linear algebraic equations:
(23)(π0,π1,…,πi,…)Q=0,(π0,π1,…,πi,…)e=1,
where Q is the infinitesimal generator of the Markov chain ξt, t≥0.To solve this system, the following numerically stable algorithm can be used.Theorem 4.
The vectorsπi, i≥0, are defined as follows:
(24)πi=π0ℱi,i≥1,
where the matrices Φi are calculated using the recurrent formulas:
(25)ℱ0=IW-,ℱi=-ℱi-1Qi-1,i(Qi,i+Qi,i+1Gi)-1,i=1,N+K-1¯,ℱN+K=-ℱN+K-1QN+K-1,N+K(Q1+Q2G)-1,ℱi=-ℱi-1Q0(Q1+Q2G)-1,i>N+K,
and the vector π0 is the unique solution to the system
(26)π0(Q0,0+Q0,1G0)=0,π0∑i=0∞‍ℱie=1.
The matricesGi are calculated using the backward recursion
(27)Gi=-(Qi+1,i+1+Qi+1,i+2Gi+1)-1Qi+1,i,i=N+K-2,N+K-3,…,0,GN+K-1=-(Q1+Q2GN+K)-1QN+K,N+K+1,GN+K=G,
where the matrix G is the minimal nonnegative solution of the matrix equation
(28)Q2G2+Q1G+Q0=O.Proof.
Proof is based on the results of the paper [8], taking into account the structure of the generator Q.


## 4. Performance Measures

As soon as the vectorsπi,i≥0, have been calculated, we are able to find various performance measures of the system under consideration.The average number of customers in the system isL=∑i=1∞iπie.The average number of busy servers isNbusy=∑i=1∞min{i,N}πie.The average number of busy type 1 servers is(29)N1busy=∑i=1N-1‍∑m=max{0,i-(N-M)}min{i,M}‍mπ(i,m)e+M∑i=N∞‍πie.The average number of busy type 2 servers is(30)N2busy=∑i=1N-1‍∑m=max{0,i-(N-M)}min{i,M}‍(i-m)π(i,m)e+(N-M)∑i=N∞‍πie.The average number of busy servers serving the real customers is(31)Nrealbusy=∑i=1N-1‍∑m=max{0,i-(N-M)}min{i,M}‍∑n=0m‍∑l=0min{i-m,N-M}‍(n+l)π(i,m,n,l)e+∑i=N∞‍∑k=0min{i-N,K}‍∑n=0M‍∑l=0N-M‍(n+l)π(i,k,n,l)e.The average number of busy servers serving the virtual customers is(32)Nvirtbusy=Nbusy-Nrealbusy.The average number of real customers in the finite buffer is(33)Nbufreal=∑i=N∞‍∑k=1min{i-N,K}‍kπ(i,k)e.The average number of virtual customers in the infinite buffer is(34)Nbufvirt=∑i=N∞‍∑k=0min{i-N,K}‍(i-k-N)π(i,k)e.The intensity of flow of real customers, who receive service in the system, is(35)λoutreal=μ1Nbusyreal.The intensity of flow of virtual customers, who receive service in the system, is calculated as(36)λoutvirt=μ2Nbusyvirt.The intensity of flow of customers, who receive service in the system, is calculated as(37)λout=λoutreal+λoutvirt.The probabilityPesc that an arbitrary real customer arrives when all servers are busy, k, k<K, real customers present in the buffer, and the customer does not join the buffer and becomes virtual is defined as
(38)Pesc=λ-1∑i=N∞‍∑k=0min{i-N,K-1}‍pkπ(i,k)(I⊗D1)e.The probabilityPent that an arbitrary real customer becomes a virtual customer at the entrance to the system due to the presence of K real customers in the buffer is given by
(39)Pent=λ-1∑i=N+K∞‍π(i,K)(I⊗D1)e.The probability that an arbitrary real customer becomes a virtual customer is calculated as(40)Pto-virt=1-λoutrealλ.The probability that, after arrival, an arbitrary real customer will go to the buffer and leave it due to impatience (becoming a virtual customer) is computed as(41)Pimp=Pto-virt-Pent-Pesc.


## 5. Distribution of the Sojourn Time of a Virtual Customer

To optimize call center operation and to solve the problem of determination of the number of different types operators, it is necessary to know the average sojourn and waiting times of an arbitrary virtual customer. We will derive the distribution of an arbitrary virtual customer’s sojourn time in terms of the Laplace-Stieltjes transform (LST).To derive the expression for the LSTv(s) of the distribution of an arbitrary virtual customer's sojourn time in the system, we use the method of collective marks (method of catastrophes; see, e.g., [9, 10]). Let us tag an arbitrary virtual customer and we will keep tracking its stay in the system. Thus, according to the method of catastrophes, v(s) has the meaning of the probability that no catastrophe arrives during the sojourn time of the tagged virtual customer.Letv(s,r,k,n,l,ν) be the probability that a catastrophe will not arrive during the rest of the tagged customer’s sojourn time in the system conditioned on the fact that, at the given moment, the position of tagged customer is r,r≥1, in the buffer for virtual customers, the number of real customers in the buffer is equal to k,k=0,K¯, the number of busy type 1 servers that serve real customers is equal to n,n=0,M¯, the number of busy type 2 servers that serve real customers is equal to l,l=0,N-M¯, and the state of the process νt,t≥0, is ν.Let us enumerate the probabilitiesv(s,r,k,n,l,ν) in the lexicographic order of the components k,n,l,ν and form the column vectors
(42)v(s,r,k,n,l)=(v(s,r,k,n,l,0),…,v(s,r,k,n,l,W))T,v(s,r,k,n)=(v(s,r,k,n,0),…,v(s,r,k,n,N-M))T,v(s,r,k)=(v(s,r,k,0),…,v(s,r,k,M))T,v(s,r)=(v(s,r,0),…,v(s,r,K))T.Theorem 5.
The vectorsv(s,r),r≥1, are calculated using the recurrent formulas
(43)v(s,1)=(sI-Ω)-1a(s),v(s,r)=(sI-Ω)-1Ω~v(s,r-1),r>1,
where
(44)Ω=-[Nμ2IK+1+αCK+1(IK+1+AK+1)]⊗I(M+1)(N-M+1)W-+(μ2-μ1)IK+1⊗[CM+1⊗IN-M+1+IM+1⊗CN-M+1]⊗IW-+I(K+1)(M+1)(N-M+1)⊗D0+[P-+(I-P-)EK+1+]⊗I(M+1)(N-M+1)⊗D1+AK+1⊗[μ1CM+1+μ2C~M+1EM+1+]⊗I(N-M+1)W-,Ω~=I(K+1)(M+1)⊗[μ1CN-M+1AN-M+1+μ2C~N-M+1]⊗IW-+I-⊗[μ1CM+1AM+1+μ2C~M+1]⊗I(N-M+1)W-,a(s)=[⊗e(M+1)(N-M+1)W-(μ1-μ2)e(K+1)(M+1)⊗CN-M+1e⊗eW-+(μ1-μ2)e~⊗CM+1e⊗e(N-M+1)W-+(Nμ2e~+(N-M)μ2(eK+1-e~))⊗e(M+1)(N-M+1)W-]μ2s+μ2.
Here,Cl=diag{0,1,…,l-1}, C~l=diag{l-1,l-2,…,0}, P-=diag{p0,p1,…,pK};Al is a square matrix of size l with all zero entries except the entries (Al)i,i-1,i=1,l-1¯, which are equal to 1;I- is a square matrix of size K+1 with all zero entries except the entry (I-)0,0=1;El+ is a square matrix of size l with all zero entries except the entries (El+)i,i+1, i=0,l-1¯, which are equal to 1;e~ is a column vector of size K+1 with all zero entries except the entry (e~)0=1.Proof.
Based on a probabilistic sense of the LST, we can obtain the following system of linear algebraic equations for calculation of the vectorsv(s,r,k,n,l):
(45)v(s,r,k,n,l)=[(s+kα+(n+l)μ1+(N-(n+l))μ2)IW--D0]-1×(δr,1[δk,0((n+l)μ1+(N-(n+l))μ2)+(1-δk,0)(lμ1+(N-M-l)μ2)]×μ2s+μ2eW-+(1-δr,1)×[δk,0(nμ1v(s,r-1,0,n-1,l)+(M-n)μ2v(s,r-1,0,n,l))+lμ1v(s,r-1,k,n,l-1)+(N-M-l)μ2v(s,r-1,k,n,l)]+(1-δk,0)×((kα+nμ1)v(s,r,k-1,n,l)+(M-n)μ2v(s,r,k-1,n+1,l))+(1-pk)D1v(s,r,k+1,n,l)+pkD1v(s,r,k,n,l)),l=0,N-M¯,n=0,M¯,k=0,K¯,r≥1,
where δi,j indicates the Kronecker delta.
Let us explain formula (45) in brief. The diagonal entries of the matrix in the first square brackets in (45) 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, and abandon of real customer from the buffer, the service completion of real or virtual customers. The nondiagonal entries of this matrix are equal to the intensities of the events which do not impact the components s,r,k,n, and l (transition of the MAP underlying process without generation of a customer). The first term in the round brackets in (45) corresponds to the case when the tagged customer has the first position in the buffer, one service process is completed and a virtual customer is chosen for service. In this case, the tagged customer is chosen for service, and the number μ2/(s+μ2) defines the probability that a catastrophe will not arrive during the service time of the tagged customer. The second term corresponds to the case when the position of the tagged customer in the buffer is greater than 1, one service process is completed, and a virtual customer is chosen for service. In this case, the position of the tagged customer in the buffer decreases by 1. The third term corresponds to the cases when the service process is completed and a real customer is chosen for service or a real customer leaves the buffer due to impatience. In this case, the number of real customers in the buffer decreases by 1. The fourth term corresponds to the case when a new customer arrives and is admitted to the system. In this case, the number of real customers in the buffer increases by 1. The fifth term corresponds to the case when a real customer arrives to the system and does not attend the buffer. In this case, the values of the components s, r, k, n, and l do not change.
Using the notation of the matrices given in Theorem5, let us rewrite system (45) into the matrix form as follows:
(46)[-(s+kα-n(μ2-μ1)+Nμ2)I(N-M+1)W-+(μ2-μ1)CN-M+1⊗IW-+IN-M+1⊗(D0+pkD1)-(s+kα-n(μ2-μ1)+Nμ2)I(N-M+1)W-]v(s,r,k,n)+δr,1[+(1-δk,0)(N-M)μ2)e(N-M+1)W-(μ1-μ2)CN-M+1e⊗eW-+(δk,0(n(μ1-μ2)+Nμ2)+(1-δk,0)(N-M)μ2)e(N-M+1)W-]×μ2s+μ2+(1-δr,1)×[(C~N-M+1μ1CN-M+1AN-M+1+μ2C~N-M+1)⊗IW-v(s,r-1,k,n)+δk,0(nμ1v(s,r-1,0,n-1)+(M-n)μ2v(s,r-1,0,n))C~N-M+1]+(1-δk,0)[(kα+nμ1)v(s,r,k-1,n)+(M-n)μ2v(s,r,k-1,n+1)(kα+nμ1)v(s,r,k-1,n)]+(1-pk)IN-M+1⊗D1v(s,r,k+1,n)=0T.
System (46) can be further rewritten into
(47){-(s+kα+Nμ2)I(M+1)(N-M+1)W-+(μ2-μ1)[CM+1⊗IN-M+1+IM+1⊗CN-M+1]⊗IW-+I(M+1)(N-M+1)⊗(D0+pkD1)-(s+kα+Nμ2)I(M+1)(N-M+1)W-}v(s,r,k)+δr,1{+(1-δk,0)μ2(N-M)e(M+1)(N-M+1)W-(μ1-μ2)eM+1⊗CN-M+1e⊗eW-+δk,0((μ1-μ2)CM+1e⊗e(N-M+1)W-+Nμ2e(M+1)(N-M+1)W-)+(1-δk,0)μ2(N-M)e(M+1)(N-M+1)W-}μ2s+μ2+(1-δr,1){IM+1⊗[μ1CN-M+1AN-M+1+μ2C~N-M+1]⊗IW-v(s,r-1,k)+δk,0(μ1CM+1AM+1+μ2C~M+1)⊗I(N-M+1)W-v(s,r-1,0)}+(1-δk,0)(kαIM+1+μ1CM+1+μ2C~M+1EM+1+)⊗I(N+M-1)W-v(s,r,k-1)+(1-pk)I(M+1)(N-M+1)⊗D1v(s,r,k+1)=0T.
Taking into account notations (44), performing algebra, we obtain the following system:
(48)(-sI+Ω)v(s,r)+(1-δr,1)Ω~v(s,r-1)+δr,1a(s)=0T,r≥1.
Based on system (48), it is easy to verify that the vectors v(s,r),r≥1, can be calculated using formulas (43).Based on a probabilistic sense of the LST and the law of total probability, one can prove the following assertion.Theorem 6.
The LSTv(s) of the distribution of an arbitrary virtual customer’s sojourn time in the system is computed as
(49)v(s)=1λoutvirt{∑i=N∞‍∑k=0min{i-N,K}‍∑n=0M‍∑l=0N-M‍pkπ(i,k,n,l)D1×v(s,i-N-k+1,k,n,l)+∑i=N∞‍∑k=1min{i-N,K}‍∑n=0M‍∑l=0N-M‍kαπ(i,k,n,l)×v(s,i-N-k+1,k-1,n,l)∑i=N∞‍∑k=0min{i-N,K}‍∑n=0M‍∑l=0N-M‍}.Corollary 7.
The average sojourn timeVsoj of an arbitrary virtual customer is computed by
(50)Vsoj=-1λoutvirt{∑i=N∞‍∑k=0min{i-N,K}‍∑n=0M‍∑l=0N-M‍pkπ(i,k,n,l)D1×v′(s,i-N-k+1,k,n,l)|s=0+∑i=N∞‍∑k=1min{i-N,K}‍∑n=0M‍∑l=0N-M‍kαπ(i,k,n,l)×v′(s,i-N-k+1,k-1,n,l)|s=0∑i=N∞‍∑k=0min{i-N,K}‍∑n=0M‍∑l=0N-M‍}.
Here, the column vectorsv′(s,r,k,n,l)|s=0 are calculated as the blocks of the vector v′(s,r)|s=0,r≥1, which can be computed as follows:
(51)v′(s,1)|s=0=Ω-1[e-a′(s)|s=0],v′(s,r)|s=0=Ω-1[e-Ω~v′(s,r-1)|s=0],r>1.Proof.
Formula (50) is based on the definition Vsoj=-v′(s)|s=0, taking into account that v(s,r)|s=0=e,r≥1.Corollary 8.
The average waiting timeVwait of an arbitrary virtual customer is defined as
(52)Vwait=Vsoj-1μ2.


## 6. Numerical Examples

Experiment 1.
In this experiment we investigate the impact of the correlation and variation in the arrival flow on the system performance measures.
For this purpose, let us introduce four MAPs, defined by the matricesD0 and D1, with the same average arrival rate λ=9 but different coefficients of correlation and variation.
The first process, coded asM0, is a stationary Poisson process. It is defined by the matrices D0=(-9) and D1=(9) and has the coefficient of correlation ccor=0 and the coefficient of variation cvar=1.
The second process, coded asIPP0 (interrupted Poisson process), is defined by the matrices
(53)D0=(-28.645400.275840.97606-0.97606),D1=(26.544741.8248200).
It has the coefficient of correlation ccor=0 and the coefficient of variation cvar=12.4.
The third process, coded asMAP0.2, defined by the matrices
(54)D0=(-12.1647600-0.39483),D1=(12.083880.080880.219910.17492)
and has ccor=0.2 and cvar=12.34.
The fourth process, coded asMAP0.4, is defined by the matrices
(55)D0=(-30.5791800.00913-0.99199),D1=(30.260680.318500.109190.87367).
It has ccor=0.4 and cvar=12.39.
We assume that the number of serversN is equal to 15, the service intensities of real and virtual customers μ1=0.7 and μ2=0.8, respectively, the buffer capacity K=15, the intensity of impatience α=0.3, and the probabilities pk=0.02(k+1),k=0,14¯.
Let us vary the number of type 1 serversM in the interval [0,15].Figure2 illustrates the dependence of the average number of customers in the system L, the average waiting time of a virtual customer Vwait, the average number of real customers Nbufreal, and the average number of virtual customers Nbufvirt in the buffer on the number of type 1 servers M for the four MAP arrival processes presented above.Average number of customers in the system, the average waiting time of a virtual customer, and the average number of real and virtual customers in the buffer as functions of the number of type 1 servers for different MAPs.
(a)
(b)
(c)
(d)The dependence of the intensity of flow of served real customersλoutreal, the intensity of flow of served virtual customers λoutvirt, the probability Pent that an arbitrary real customer becomes a virtual customer at the entrance to the system due to the presence of K real customers in the buffer and the probability Pimp that an arbitrary real customer after arrival will go to the buffer and leave it due to impatience on the number of type 1 servers M for different MAP arrival processes is presented in Figure 3.Intensities of output flow of real and virtual customers, who received service in the system, and the loss probabilitiesPent and Pimp as functions of the number of type 1 servers M for different MAPs.
(a)
(b)
(c)
(d)From Figures2 and 3, one can conclude that the coefficients of correlation and variation in the arrival process have a profound impact on the system performance measures. As the coefficients of correlation and variation increase, the intensities λoutvirt and λoutreal decrease essentially, while the average number of customers in the system L, the average waiting time of an arbitrary virtual customer Vwait, and the loss probability Pent increase. Thus, increasing the coefficients of correlation and variation has a negative effect on system operation.Additionally, one can conclude, that with increasing the number of type 1 serversM the intensity λoutreal increases, while intensity λoutvirt decreases. However, at the same time, the average waiting time of an arbitrary virtual customer Vwait increases essentially with increasing M, that has a negative impact on the quality of service of virtual customers.The results of this experiment also show that Little's formula holds true for the system under study; that is,(56)Vwait=NvirtbufferλPto-virt=Nvirtbufferλoutvirt.Additionally, the numerical results show that another equality holds true:Nbufreal/λ=Pimp/α. To prove this equality, let us rewrite it in the form αNbufreal=λPimp. The left and right sides of this equality define the intensity of real customers' leaving the buffer due to impatience, so the equality is correct.Experiment 2.
The issue of the optimal choice of the number of each type operator is very important for the effective organization of the call center performance. In this experiment, we numerically solve the problem of optimal choice of the number of type 1 serversM. The aim of optimization is the maximization of a cost criterion (one of the possible criteria):
(57)Λ(M)=λoutreal+(1-q1)q2λoutvirt.
Here, q1 is the probability that an operator will not succeed to connect to a customer (customer's phone is busy or does not answer) and q2 is the probability that a virtual customer will be interested in the services proposed by an operator. It is evident that the probability q2 depends on the waiting time. Thus, Λ(M) determines the effective intensity of the flow of customers, who successfully receive service in the system.
Let us fixq1=0.1, and
(58)q2={0,Vwait≥40,1-Vwait40,Vwait<40.
Our goal is to find the optimal valueM* that provides the maximum value of the cost criterion Λ(M). The value of the cost criterion Λ(M) for different values of the number of type 1 servers M and different arrival processes is presented in Table 1. The optimal values of Λ(M) for each arrival process are printed in bold type. Additionally, Table 1 contains the values of the cost criterion Λ(M) for the call center that does not provide a call-back option. In this case, Λ(M)=λout and can be found based on the results of the paper [11].
As it is seen from Table1, for the stationary Poisson arrival process, the maximum value of Λ(M) is equal to 8.91163 for M=15; that is, when all operators are of type 1 and there are no type 2 operators. For the interrupted Poisson process, the maximum value of Λ(M) is equal to 8.32646, when the system has seven type 1 operators and eight type 2 operators. For MAP0.2 arrival process, the maximum value of Λ(M) is equal to 8.39144, when the system has three type 1 operators and twelve type 2 operators and for MAP0.4 arrival process the maximum value of Λ(M) is equal to 5.84454 when the system has only type 2 operators.
From Table1, one can conclude that providing of a call-back for lost customers makes sense and delivers benefits to a call center. For any choice of the number of type 1 operators, the values of Λ(M) are essentially higher for the system with call-backs than without call-backs. Another important conclusion is that the coefficients of correlation and variation in the arrival process strongly affect the values of cost criterion Λ(M) and the values of the optimal number of type 1 operators. So, for the correct prediction of system performance measures it is necessary to take into account the coefficients of correlation and variation in the arrival process.Table 1
Value of the flow intensity of customers, who receive service in the system, for different values of the number of type 1 serversM and different arrival processes.
Λ
(
M
)
M
0
IP
P
0
MA
P
0.2
MA
P
0.4
Without call-back
8.59707
7.22160
7.76157
4.89795
M
=
0
8.89987
8.31396
8.38843
5.84454
M
=
1
8.90027
8.31731
8.39022
5.77657
M
=
2
8.90068
8.32023
8.39130
5.71393
M
=
3
8.90112
8.32264
8.39144
5.65655
M
=
4
8.90158
8.32449
8.39039
5.60255
M
=
5
8.90208
8.32575
8.38787
5.55049
M
=
6
8.90263
8.32640
8.38354
5.49944
M
=
7
8.90325
8.32646
8.37706
5.44876
M
=
8
8.90396
8.32595
8.36812
5.39803
M
=
9
8.9048
8.32490
8.35651
5.34692
M
=
10
8.90577
8.32335
8.3422
5.29522
M
=
11
8.90688
8.32135
8.3255
5.24273
M
=
12
8.90809
8.31894
8.30702
5.18933
M
=
13
8.90932
8.31617
8.28766
5.13492
M
=
14
8.91052
8.31309
8.26835
5.07942
M
=
15
8.91163
8.30977
8.24991
5.02276


## 7. Conclusion

In this paper, a novel multiserver queueing system with impatient heterogeneous customers and two types of servers, which can be suited for modeling a real-life call center with call-backs for lost customers, is investigated. The ergodicity condition is derived. An efficient algorithm for calculating the stationary probabilities of the system states is presented. The main performance measures and the Laplace-Stieltjes transform of the distribution of the sojourn time of an arbitrary virtual customer are derived. The numerical results confirm the importance of the analysis of queueing models with the Markovian arrival process because correlation and variation in the arrival process essentially impact system operation. The problem of optimal choice of the number of each type servers is solved numerically.



---

*Source: 983723-2013-11-28.xml*