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

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. Type 𝑟, 𝑟=1,2, servers serve type 𝑟 customers if there are any in the system and serve type 𝑟^′, 𝑟^′=1,2,  𝑟^′≠𝑟, 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.

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 Section 2, 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 an 𝑁-server queuing system with two types of customers, a buffer of capacity 𝐾 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 𝜈_𝑡, 𝑡≥0, which is an irreducible continuous-time Markov chain with the state space {0,1,…,𝑊}. Arrivals occur only at the epochs of transitions in the underlying process 𝜈_𝑡, 𝑡≥0. The intensities of transitions of the process 𝜈_𝑡, 𝑡≥0, that are accompanied (not accompanied) by the arrival of a customer are defined by the square matrix 𝐷₁(𝐷₀) of size 𝑊=𝑊+1.

The matrix 𝐷(1)=𝐷₀+𝐷₁ is an infinitesimal generator of the process 𝜈_𝑡, 𝑡≥0. The stationary distribution vector 𝜽 of this process satisfies the system of equations 𝜽𝐷(1)=𝟎, 𝜽𝐞=1. Here and throughout this paper, 𝟎 is a zero row vector, and e denotes a unit column vector.

The average intensity 𝜆 (fundamental rate) of the MAP is defined by 𝜆=𝜽𝐷₁𝐞. The coefficient of variation 𝑐_var of intervals between customer arrivals is calculated as 𝑐²_var=2𝜆𝜽(−𝐷₀)⁻¹𝐞−1. The coefficient of correlation 𝑐_cor of intervals between successive arrivals is given as 𝑐_cor=(𝜆𝜽(−𝐷₀)⁻¹(𝐷(1)−𝐷₀)(−𝐷₀)⁻¹𝐞−1)/𝑐²_var.

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 are 𝑘, 𝑘=0,𝐾−1, real customers in the finite buffer during an arbitrary real customer arrival epoch, the arriving customer leaves the system with probability 𝑝_𝑘, 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, 𝑝_𝐾=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 is 𝑀, 𝑀=0,𝑁, and the number of type 2 servers is 𝑁−𝑀.

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 𝜇₁ and 𝜇₂, respectively.

3. The Process of System States

Let 𝑖_𝑡 be the number of customers in the system, 𝑘_𝑡 the number of real customers in the finite buffer, 𝑚_𝑡 the number of busy type 1 servers, 𝑛_𝑡 the number of busy type 1 servers that serve the real customers, 𝑙_𝑡 the number of busy type 2 servers that serve the real customers, and 𝜈_𝑡 the state of the underlying process of the MAP during the epoch 𝑡, 𝑡≥0.

The behavior of the system under consideration can be described in terms of the regular irreducible continuous-time Markov chain 𝜉_𝑡={𝑖_𝑡,𝑘_𝑡,𝑚_𝑡,𝑛_𝑡,𝑙_𝑡,𝜈_𝑡}, 𝑡≥0, with the state space {𝑖,𝑚,𝑛,𝑙,𝜈},𝑖<𝑁,𝑚=max{0,𝑖−(𝑁−𝑀)},min{𝑖,𝑀},𝑛=0,𝑚,𝑙=0,min{𝑖−𝑚,𝑁−𝑀},𝜈={}0,𝑊𝑖,𝑘,𝑛,𝑙,𝜈,𝑖≥𝑁,𝑘=0,min{𝑖−𝑁,𝐾},𝑛=0,𝑀,𝑙=0,𝑁−𝑀,𝜈=0,𝑊. (1)

Let us introduce the operator 𝒯(𝐴_{𝑘,𝑘−1},𝑘=𝑘₁,𝐾₁;𝐴_𝑘,𝑘,𝑘=𝑘₂,𝐾₂;𝐴_𝑘,𝑘+1,𝑘=𝑘₃,𝐾₃) which define a block tridiagonal matrix with subdiagonal blocks 𝐴_{𝑘,𝑘−1}, diagonal blocks 𝐴_𝑘,𝑘, and updiagonal blocks 𝐴_𝑘,𝑘+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 𝜉_𝑡 in the lexicographic order of the components 𝑘, 𝑚, 𝑛, 𝑙, and 𝜈. Let 𝑄 be a generator of the Markov chain 𝜉_𝑡,𝑡≥0, consisting of the blocks 𝑄_𝑖,𝑗, which define the transition rates of this chain from the states having value 𝑖 of the first component to the states having value 𝑗 of this component of the chain. The diagonal entries of the matrices 𝑄_𝑖,𝑖 are negative, and the modulus of diagonal entries defines the total intensity of leaving the corresponding state of the Markov chain 𝜉_𝑡,𝑡≥0.

Lemma 1. The infinitesimal generator 𝑄 of the Markov chain 𝜉_𝑡,𝑡≥0, has a block tridiagonal structure:𝑄𝑄=^0,0𝑄_0,1𝑄𝑂𝑂⋯𝑂𝑂𝑂𝑂𝑂⋯_1,0𝑄_1,1𝑄_1,2𝑂⋯𝑂𝑂𝑂𝑂𝑂⋯𝑂𝑄_2,1𝑄_2,2𝑄_2,3⋯𝑂𝑂𝑂𝑂𝑂⋯⋮⋮⋱⋱⋱⋮⋮⋮⋮⋮⋯𝑂𝑂𝑂𝑂⋯𝑄_{𝑁+𝐾−1,𝑁+𝐾−2}𝑄_{𝑁+𝐾−1,𝑁+𝐾−1}𝑄_{𝑁+𝐾−1,𝑁+𝐾}𝑂𝑂⋯𝑂𝑂𝑂𝑂⋯𝑂𝑄_{𝑁+𝐾,𝑁+𝐾−1}𝑄₁𝑄₂𝑂⋯𝑂𝑂𝑂𝑂⋯𝑂𝑂𝑄₀𝑄₁𝑄₂⋯⋮⋯⋮⋮⋮⋮⋮⋮⋱⋱⋱. (2)
The nonzero blocks 𝑄_𝑖,𝑗, 𝑖,𝑗≥0, have the following form: 𝑄_𝑖,𝑖=diagdiagdiag𝐴(𝑖,𝑚,𝑛,𝑙),𝑙=0,min{𝑁−𝑀,𝑖−𝑚},𝑛=0,𝑚,𝑚=𝑖^∗,min{𝑖,𝑀},𝑖=𝑄0,𝑁−1,_𝑖,𝑖=𝒯𝐵_{𝑘,𝑘−1},𝑘=1,𝑖−𝑁;𝐵_𝑘,𝑘,𝑘=𝑂0,𝑖−𝑁;_𝑘,𝑘=0,𝑖−𝑁−1,𝑖=𝑄𝑁,𝑁+𝐾−1,₁=𝒯𝐵_{𝑘,𝑘−1},𝑘=1,𝐾;𝐵_𝑘,𝑘,𝑘=𝑂0,𝐾;_𝑘,𝑘=𝑄0,𝐾−1,𝑖≥𝑁+𝐾,_𝑖,𝑖+1=𝒯𝑂_𝑚,𝑚=1,𝑖;𝑂^(𝑖)_𝑚,𝑚,𝑚=𝐶0,𝑖;_{(𝑖)𝑚,𝑚+1},𝑚=𝐶0,𝑖,𝑖<𝑁−𝑀;diag^(𝑖)_𝑚,𝑚+1,𝑚=𝑖^∗,𝑖,𝑖≥𝑁−𝑀,𝑖=𝑄0,𝑀−1,_𝑖,𝑖+1=𝒯𝑂_𝑚,𝑚=1,𝑀;𝐶^(𝑖)_𝑚,𝑚,𝑚=𝐶0,𝑀;_{(𝑖)𝑚,𝑚+1},𝑚=,𝒯𝐶0,𝑀−1𝑖<𝑁−𝑀;_{(𝑖)𝑚,𝑚},𝑚=𝑖^∗𝐶+1,𝑀;_{(𝑖)𝑚,𝑚+1},𝑚=𝑖^∗𝑂,𝑀−1;_𝑚,𝑚=𝑖^∗,𝑀−2,𝑖≥𝑁−𝑀,𝑖=𝑄𝑀,𝑁−1;_𝑖,𝑖+1=𝒯𝑂_𝑘,𝑘=1,𝑖−𝑁;𝐻_𝑘,𝑘,𝑘=𝐻0,𝑖−𝑁;_𝑘,𝑘+1,𝑘=0,𝑖−𝑁,𝑖=𝑄𝑁,𝑁+𝐾−1,₂=𝒯𝑂_𝑘,𝑘=1,𝐾;𝐻_𝑘,𝑘,𝑘=𝐻0,𝐾;_𝑘,𝑘+1,𝑘=𝑄0,𝐾−1,𝑖≥𝑁+𝐾,_{𝑖,𝑖−1}=𝒯𝐸^(𝑖)_{𝑚,𝑚−1},𝑚=1,𝑖;𝐸^(𝑖)_𝑚,𝑚,𝑚=𝑂0,𝑖−1;_𝑚,𝑚=𝒯𝑂0,𝑖−2,𝑖≤𝑁−𝑀;_𝑚,𝑚=𝑖^∗+1,𝑖;𝐸^(𝑖)_{𝑚,𝑚−1},𝑚=𝑖^∗𝐸,𝑖;_{(𝑖)𝑚,𝑚},𝑚=𝑖^∗,𝑖−1,𝑖>𝑁−𝑀,𝑖=𝑄1,𝑀,_{𝑖,𝑖−1}=𝒯𝐸^(𝑖)_{𝑚,𝑚−1},𝑚=𝐸1,𝑀;_{(𝑖)𝑚,𝑚},𝑚=0,𝑀;𝑂_𝑚,𝑚=,𝒯𝑂0,𝑀−1𝑖≤𝑁−𝑀;_𝑚,𝑚=𝑖^∗𝐸+1,𝑀;_{(𝑖)𝑚,𝑚−1},𝑚=𝑖^∗,𝑀;𝐸^(𝑖)_𝑚,𝑚,𝑚=𝑖^∗,𝑀,𝑖>𝑁−𝑀,𝑖=𝑄𝑀+1,𝑁,_{𝑖,𝑖−1}=𝒯𝑍^(𝑖)_{𝑘,𝑘−1},𝑘=1,𝑖−𝑁;𝑍_𝑘,𝑘,𝑘=𝑂0,𝑖−𝑁−1;_𝑘,𝑘=0,𝑖−𝑁−2,𝑖=𝑄𝑁+1,𝑁+𝐾,₀=𝒯𝑉_{𝑘,𝑘−1},𝑘=1,𝐾;𝑉_𝑘,𝑘,𝑘=𝑂0,𝐾;_𝑘,𝑘=0,𝐾−1,𝑖>𝑁+𝐾. (3)
Here,  𝐼 is the identity matrix, 𝑂 is a zero matrix of appropriate dimension, and 𝑂_𝑙, 𝑙=𝑘,𝐾, is the sequence which consists of 𝐾−𝑘+1 zero matrices of appropriate dimension: 𝑖^∗=max{0,𝑖−(𝑁−𝑀)}. (4)⊗indicates the Kronecker product of matrices: 𝐴()=−)𝜇𝑖,𝑚,𝑛,𝑙(𝑛+𝑙₁+(𝑖−𝑛−𝑙)𝜇₂𝐼_𝑊+𝐷₀𝐷,𝑖−𝑚≤𝑁−𝑀,₀,𝑖−𝑚>𝑁−𝑀; (5) 𝐵_𝑘,𝑘=−𝐵⊗𝐼_𝑊+𝐼_{(𝑀+1)(𝑁−𝑀+1)}⊗(−𝑘𝛼𝐼_𝑊+𝐷₀) where 𝐵=diag{diag{(𝑛+𝑙)𝜇₁+(𝑁−𝑛−𝑙)𝜇₂,𝑙=0,𝑁−𝑀},𝑛=0,𝑀}; 𝐵_{𝑘,𝑘−1}=𝑘𝛼𝐼_{(𝑀+1)(𝑁−𝑀+1)𝑊}; (6) 𝐶^(𝑖)_𝑚,𝑚+1=𝒯(𝑂_𝑛,𝑛=1,𝑚;𝑂_𝑛,𝑛=0,𝑚;𝐶^(𝑖,𝑚)_𝑛,𝑛+1,𝑛=0,𝑚) where 𝐶^(𝑖,𝑚)_𝑛,𝑛+1=𝐼_{min{𝑖−𝑚,𝑁−𝑀}+1}⊗𝐷₁;𝐶^(𝑖)_𝑚,𝑚=𝑂_{(𝑖+1)𝑊×(𝑖+2)𝑊}𝑂,𝑚=0,𝑖<𝑁−𝑀;_{(𝑁−𝑀+1)𝑊×(𝑁−𝑀+1)𝑊}𝑂,𝑚=0,𝑖≥𝑁−𝑀;_{(𝑚+1)(min{𝑖−𝑚,𝑁−𝑀}+1)𝑊×(𝑚+1)(min{𝑖−𝑚+1,𝑁−𝑀}+1)𝑊},𝑚=𝐼1,𝑀−1;_𝑀+1⊗𝐼⁺_𝑖⊗𝐷₁,𝑚=𝑀, (7)where 𝐼⁺_𝑖 is a matrix of size (𝑖−𝑀+1)×(𝑖−𝑀+2) with nonzero entries (𝐼⁺_𝑖)_𝑙,𝑙+1,𝑙=0,𝑖−𝑀, which are equal to 1. Consider𝑂^(𝑖)_𝑚,𝑚=𝑂_{(𝑖+1)𝑊×(𝑖+2)𝑊}𝑂,𝑚=0,_{(𝑚+1)(min{𝑖−𝑚,𝑁−𝑀}+1)𝑊×(𝑚+1)(min{𝑖−𝑚+1,𝑁−𝑀}+1)𝑊,}𝑚=𝐻1,𝑖;_𝑘,𝑘=𝐼_{(𝑀+1)(𝑁−𝑀+1)}⊗𝑝_𝑘𝐷₁;𝐻_𝑘,𝑘+1=𝐼_{(𝑀+1)(𝑁−𝑀+1)}⊗1−𝑝_𝑘𝐷₁;𝐸^(𝑖)_𝑚,𝑚=𝐼_𝑚+1⊗𝒯𝜇₁𝑙,𝑙=1,𝑖−𝑚;𝜇₂(𝑖−𝑚−𝑙),𝑙=0,𝑖−𝑚−1;𝑂_𝑙,𝑙=0,𝑖−𝑚−2⊗𝐼_𝑊;𝐸^(𝑖)_{𝑚,𝑚−1}=𝒯𝜇₁𝑙,𝑙=1,𝑚;𝜇₂(𝑚−𝑙),𝑙=0,𝑚−1;𝑂_𝑙,𝑙=0,𝑚−2⊗𝐼_{(𝑖−𝑚+1)𝑊};𝑍_𝑘,𝑘=𝒯𝜇₁𝑙𝐼_{𝑊(𝑁−𝑀+1)},𝑙=1,𝑀;Ψ+𝜇₂(𝑀−𝑙)𝐼_{(𝑁−𝑀+1)𝑊},𝑙=𝑂0,𝑀;_𝑙,𝑙=𝐼0,𝑀−1,𝑘=0;_𝑀+1⊗Ψ,𝑘=𝜇1,𝐾−1,Ψ=𝒯₁𝑙,𝑙=1,𝑁−𝑀;𝜇₂(𝑁−𝑀−𝑙),𝑙=0,𝑁−𝑀;𝑂_𝑙,𝑙=0,𝑁−𝑀−1⊗𝐼_𝑊;𝑍^(𝑖)_{𝑘,𝑘−1}=Φ,𝑘=1,𝑖−𝑁−1;Φ+𝐼_𝑀+1⊗𝒯𝑂_𝑙,𝑙=1,𝑁−𝑀;𝜇₁𝑙,𝑙=𝜇0,𝑁−𝑀;₂(𝑁−𝑀−𝑙),𝑙=0,𝑁−𝑀−1⊗𝐼_𝑊𝑂,𝑘=𝑖−𝑁,Φ=𝒯_𝑛,𝑛=1,𝑀;𝜇₁𝑛,𝑛=0,𝑀;𝜇₂(𝑀−𝑛),𝑛=0,𝑀−1⊗𝐼_{(𝑁−𝑀+1)𝑊};𝑉_𝑘,𝑘=𝒯𝜇₁𝑛𝐼_{𝑁−𝑀+1},𝑛=1,𝑀;𝑉_𝑛,𝑛,𝑛=0,𝑀;𝑂_𝑛,𝑛=0,𝑀−1⊗𝐼_𝑊𝐼,𝑘=0,_𝑀+1⊗Ψ,𝑘=1,𝐾, (8)where 𝑉_𝑛,𝑛=𝒯(𝜇₁𝑙,𝑙=1,𝑁−𝑀;𝜇₂(𝑁−𝑛−𝑙),𝑙=0,𝑁−𝑀;𝑂_𝑙,𝑙=0,𝑁−𝑀). Consider 𝑉_{𝑘,𝑘−1}=𝒯𝑂_𝑛,𝑛=1,𝑀;𝜇₁𝑛,𝑛=𝜇0,𝑀;₂(𝑀−𝑛),𝑛=0,𝑀−1⊗𝐼_{(𝑁−𝑀+1)𝑊}. (9)

The proof of the lemma is implemented by the analysis of all transitions of the Markov chain 𝜉_𝑡,𝑡≥0, during the interval of infinitesimal length and rewriting the intensities of these transitions in block matrix form.

Corollary 2. The Markov chain 𝜉_𝑡, 𝑡≥0, belongs to the class of continuous time quasi-birth-and-death processes; see, for example, [7].

Theorem 3. The Markov chain 𝜉_𝑡,𝑡≥0, is ergodic if the following inequality holds true: 𝐲𝒜𝐵⊗𝐼_𝑊𝐞>𝐲𝒜𝐼_𝑀+1⊗𝐷₁𝐞, (10) where the vector 𝐲 is the unique solution to the system 𝐲𝑄_0,0+𝑇₀𝑄_1,0=𝟎,𝐲(𝐼+𝒜)𝐞=1. (11)
Here, 𝒜=^𝐾𝑘_𝑘=0𝑙=0𝑇_𝑙; (12) the matrices 𝑇_𝑘,𝑘=0,𝐾−1, are calculated using the backward recursion 𝑇_𝐾−1=−𝑄_{𝐾−1,𝐾}𝑄_𝐾,𝐾⁻¹,𝑇_𝑘=−𝑄_𝑘,𝑘+1𝑄_{𝑘+1,𝑘+1}+𝑇_𝑘+1𝑄_{𝑘+2,𝑘+1}⁻¹,𝑘=𝐾−2,𝐾−3,…,0; (13)(14) the matrices 𝑄_𝑘,𝑘, 𝑄_𝑘,𝑘+1, and 𝑄_{𝑘,𝑘−1} are defined as follows: 𝑄_0,0=𝒯𝜇₁𝑛,𝑛=1,𝑀;−𝜇₁𝑛,𝑛=0,𝑀;𝑂_𝑛,𝑛=0,𝑀−1⊗𝐼_𝑊+𝐼_𝑀+1⊗𝐷₀+𝑝₀𝐷₁,𝑄_𝑘,𝑘=diag−𝜇₁𝑛−𝜇₂(𝑀−𝑛),𝑛=0,𝑀⊗𝐼_𝑊−𝑘𝛼𝐼_(𝑀+1)𝑊+𝐼_𝑀+1⊗𝐷₀+𝑝_𝑘𝐷₁,𝑘=1,𝐾,𝑄_𝑘,𝑘+1=1−𝑝_𝑘𝐼_𝑀+1⊗𝐷₁,𝑘=0,𝐾−1,𝑄_{𝑘,𝑘−1}=𝑘𝛼𝐼_(𝑀+1)𝑊+𝒯𝑂_𝑛,𝑛=1,𝑀;𝜇₁𝑛,𝑛=()𝜇0,𝑀;𝑀−𝑛₂,𝑛=0,𝑀−1⊗𝐼_𝑊,𝑘=1,𝐾. (15)

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) where the vector 𝐱 is the unique solution to the following system of linear algebraic equations: 𝐱𝑄₀+𝑄₁+𝑄₂=𝟎,𝐱𝐞=1. (17)
Let the row vector 𝐱 have the form 𝐱=𝐱₀,𝐱₁,…,𝐱_𝐾, (18) where 𝐱_𝑘=𝐱_𝑘⊗𝐞̂⊗𝜽, 𝑘=0,𝐾. Here, 𝐞̂ is a row vector of size 𝑁−𝑀+1 with all zero entries except the entry (𝐞̂)₀ which is equal to 1, 𝜽 is the stationary distribution vector of the process 𝜈_𝑡,𝑡≥0, and 𝐱_𝑘, 𝑘=0,𝐾 is a row vector of size 𝑀+1 such as the vector (𝐱₀,…,𝐱_𝐾) is stochastic.
By direct substitution, we verify that the vector 𝐱 from (18) provides the unique solution of system (17) if the vector 𝐱=(𝐱₀,𝐱₁,…,𝐱_𝐾), where 𝐱_𝑘=𝐱_𝑘⊗𝜽 is the unique solution to the following system of linear algebraic equations: 𝐱𝑄=𝐱,𝐱𝐞=1. (19)
Here, 𝑄 is a block tridiagonal matrix with subdiagonal blocks 𝑄_{𝑘,𝑘−1}, diagonal blocks 𝑄_𝑘,𝑘, and updiagonal blocks 𝑄_𝑘,𝑘+1 defined by formulas (15).
To calculate the vectors 𝐱_𝑘, 𝑘=0,𝐾, we use the numerically stable algorithm which consists of the following steps.
Step 1. Calculate the matrices 𝑇_𝑘, 𝑘=0,𝐾−1, using the backward recursion (14) under the initial condition (13).
Step 2. Calculate the vector 𝐱₀ as the unique solution to the following system: 𝐱₀𝑄_0,0+𝑇₀𝑄_1,0=𝟎,𝐱₀𝐼+^𝐾𝑘_𝑘=0𝑙=0𝑇_𝑙𝐞=1. (20)
Step 3. Calculate the vectors 𝐱_𝑘=𝐱₀∏^𝑘_𝑙=0𝑇_𝑙, 𝑘=1,𝐾.
Substituting the vector 𝐱 from (18) into inequality (16) and performing some algebra, we obtain the inequality ^𝐾_𝑘=0𝐱_𝑘⊗𝜽𝐵⊗𝐼_𝑊𝐞>^𝐾_𝑘=0𝐱_𝑘⊗𝜽𝐼_𝑀+1⊗𝐷₁𝐞 (21) or ^𝐾_𝑘=0𝐱_𝑘𝐵⊗𝐼_𝑊𝐞>^𝐾_𝑘=0𝐱_𝑘𝐼_𝑀+1⊗𝐷₁𝐞. (22)
Let the matrix 𝒜 be given by formulas (12) and denote the vector 𝐱₀ as 𝐲. 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 𝜉_𝑡 is fulfilled, then the stationary probabilities of system states 𝜋(𝑖,𝑚,𝑛,𝑙,𝜈), 𝑖<𝑁, 𝑚=max{0,𝑖−(𝑁−𝑀)},min{𝑖,𝑀},𝑛=0,𝑚,  𝑙=0,min{𝑖−𝑚,𝑁−𝑀}, 𝜈=0,𝑊, and 𝜋(𝑖,𝑘,𝑛,𝑙,𝜈), 𝑖≥𝑁, 𝑘=0,min{𝑖−𝑁,𝐾}, 𝑛=0,𝑀, 𝑙=0,𝑁−𝑀, 𝜈=0,𝑊, exist. Let us form the row vectors 𝝅_𝑖 of these probabilities enumerated in the lexicographic order of the components (𝑚,𝑛,𝑙,𝜈) if 𝑖<𝑁 and in the lexicographic order of the components (𝑘,𝑛,𝑙,𝜈) if 𝑖≥𝑁.

It is well known that the probability vectors 𝝅_𝑖, 𝑖≥0, satisfy the following system of linear algebraic equations: 𝝅₀,𝝅₁,…,𝝅_𝑖,…𝑄=𝟎,𝝅₀,𝝅₁,…,𝝅_𝑖,…𝐞=1, (23) where 𝑄 is the infinitesimal generator of the Markov chain 𝜉_𝑡, 𝑡≥0.

To solve this system, the following numerically stable algorithm can be used.

Theorem 4. The vectors 𝝅_𝑖, 𝑖≥0, are defined as follows: 𝝅_𝑖=𝝅₀ℱ_𝑖,𝑖≥1, (24) where the matrices Φ_𝑖 are calculated using the recurrent formulas: ℱ₀=𝐼_𝑊,ℱ_𝑖=−ℱ_𝑖−1𝑄_{𝑖−1,𝑖}𝑄_𝑖,𝑖+𝑄_𝑖,𝑖+1𝐺_𝑖⁻¹,𝑖=ℱ1,𝑁+𝐾−1,_𝑁+𝐾=−ℱ_{𝑁+𝐾−1}𝑄_{𝑁+𝐾−1,𝑁+𝐾}𝑄₁+𝑄₂𝐺⁻¹,ℱ_𝑖=−ℱ_𝑖−1𝑄₀𝑄₁+𝑄₂𝐺⁻¹,𝑖>𝑁+𝐾, (25) and the vector 𝝅₀ is the unique solution to the system 𝝅₀𝑄_0,0+𝑄_0,1𝐺₀=𝟎,𝝅₀^∞_𝑖=0ℱ_𝑖𝐞=1. (26)
The matrices 𝐺_𝑖 are calculated using the backward recursion 𝐺_𝑖=−𝑄_{𝑖+1,𝑖+1}+𝑄_{𝑖+1,𝑖+2}𝐺_𝑖+1⁻¹𝑄_𝑖+1,𝑖,𝐺𝑖=𝑁+𝐾−2,𝑁+𝐾−3,…,0,_{𝑁+𝐾−1}=−𝑄₁+𝑄₂𝐺_𝑁+𝐾⁻¹𝑄_{𝑁+𝐾,𝑁+𝐾+1},𝐺_𝑁+𝐾=𝐺, (27) where the matrix 𝐺 is the minimal nonnegative solution of the matrix equation 𝑄₂𝐺²+𝑄₁𝐺+𝑄₀=𝑂. (28)

Proof. Proof is based on the results of the paper [8], taking into account the structure of the generator 𝑄.

4. Performance Measures

As soon as the vectors 𝝅_𝑖,𝑖≥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 is 𝐿=∑^∞_𝑖=1𝑖𝝅_𝑖𝐞.

The average number of busy servers is 𝑁_busy=∑^∞_𝑖=1min{𝑖,𝑁}𝝅_𝑖𝐞.

The average number of busy type 1 servers is 𝑁^busy₁=^𝑁−1_𝑖=1^{min{𝑖,𝑀}}_{𝑚=max{0,𝑖−(𝑁−𝑀)}}𝑚𝝅(𝑖,𝑚)𝐞+𝑀^∞_𝑖=𝑁𝝅_𝑖𝐞. (29)

The average number of busy type 2 servers is 𝑁^busy₂=^𝑁−1_𝑖=1^{min{𝑖,𝑀}}_{𝑚=max{0,𝑖−(𝑁−𝑀)}}()𝝅()𝐞+()𝑖−𝑚𝑖,𝑚𝑁−𝑀^∞_𝑖=𝑁𝝅_𝑖𝐞. (30)

The average number of busy servers serving the real customers is 𝑁^busy_real=^𝑁−1_𝑖=1^{min{𝑖,𝑀}𝑚}_{𝑚=max{0,𝑖−(𝑁−𝑀)}𝑛=0}^{min{𝑖−𝑚,𝑁−𝑀}}_𝑙=0()𝝅()𝐞+𝑛+𝑙𝑖,𝑚,𝑛,𝑙_{∞𝑖=𝑁min{𝑖−𝑁,𝐾}𝑀𝑘=0𝑛=0𝑁−𝑀𝑙=0}(𝑛+𝑙)𝝅(𝑖,𝑘,𝑛,𝑙)𝐞. (31)

The average number of busy servers serving the virtual customers is 𝑁^busy_virt=𝑁_busy−𝑁^busy_real. (32)

The average number of real customers in the finite buffer is 𝑁^real_buf=^∞_𝑖=𝑁^{min{𝑖−𝑁,𝐾}}_𝑘=1𝑘𝝅(𝑖,𝑘)𝐞. (33)

The average number of virtual customers in the infinite buffer is 𝑁^virt_buf=^∞_𝑖=𝑁^{min{𝑖−𝑁,𝐾}}_𝑘=0(𝑖−𝑘−𝑁)𝝅(𝑖,𝑘)𝐞. (34)

The intensity of flow of real customers, who receive service in the system, is 𝜆^real_out=𝜇₁𝑁^real_busy. (35)

The intensity of flow of virtual customers, who receive service in the system, is calculated as 𝜆^virt_out=𝜇₂𝑁^virt_busy. (36)

The intensity of flow of customers, who receive service in the system, is calculated as 𝜆_out=𝜆^real_out+𝜆^virt_out. (37)

The probability 𝑃_esc that an arbitrary real customer arrives when all servers are busy, 𝑘, 𝑘<𝐾, real customers present in the buffer, and the customer does not join the buffer and becomes virtual is defined as 𝑃_esc=𝜆^∞−1_𝑖=𝑁^{min{𝑖−𝑁,𝐾−1}}_𝑘=0𝑝_𝑘𝝅(𝑖,𝑘)𝐼⊗𝐷₁𝐞. (38)

The probability 𝑃_ent that an arbitrary real customer becomes a virtual customer at the entrance to the system due to the presence of 𝐾 real customers in the buffer is given by 𝑃_ent=𝜆^∞−1_{𝑖=𝑁+𝐾}𝝅(𝑖,𝐾)𝐼⊗𝐷₁𝐞. (39)

The probability that an arbitrary real customer becomes a virtual customer is calculated as 𝑃_to-virt𝜆=1−^realout𝜆. (40)

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 𝑃_imp=𝑃_to-virt−𝑃_ent−𝑃_esc. (41)

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 LST 𝑣(𝑠) 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, 𝑣(𝑠) has the meaning of the probability that no catastrophe arrives during the sojourn time of the tagged virtual customer.

Let 𝑣(𝑠,𝑟,𝑘,𝑛,𝑙,𝜈) 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 𝑟,𝑟≥1, in the buffer for virtual customers, the number of real customers in the buffer is equal to 𝑘,𝑘=0,𝐾, the number of busy type 1 servers that serve real customers is equal to 𝑛,𝑛=0,𝑀, the number of busy type 2 servers that serve real customers is equal to 𝑙,𝑙=0,𝑁−𝑀, and the state of the process 𝜈_𝑡,𝑡≥0, is 𝜈.

Let us enumerate the probabilities 𝑣(𝑠,𝑟,𝑘,𝑛,𝑙,𝜈) in the lexicographic order of the components 𝑘,𝑛,𝑙,𝜈 and form the column vectors 𝐯(𝑠,𝑟,𝑘,𝑛,𝑙)=(𝑣(𝑠,𝑟,𝑘,𝑛,𝑙,0),…,𝑣(𝑠,𝑟,𝑘,𝑛,𝑙,𝑊))^𝑇,𝐯(𝑠,𝑟,𝑘,𝑛)=(𝐯(𝑠,𝑟,𝑘,𝑛,0),…,𝐯(𝑠,𝑟,𝑘,𝑛,𝑁−𝑀))^𝑇,𝐯(𝑠,𝑟,𝑘)=(𝐯(𝑠,𝑟,𝑘,0),…,𝐯(𝑠,𝑟,𝑘,𝑀))^𝑇,𝐯(𝑠,𝑟)=(𝐯(𝑠,𝑟,0),…,𝐯(𝑠,𝑟,𝐾))^𝑇. (42)

Theorem 5. The vectors 𝐯(𝑠,𝑟),𝑟≥1, are calculated using the recurrent formulas 𝐯(𝑠,1)=(𝑠𝐼−Ω)⁻¹𝐚(𝑠),𝐯(𝑠,𝑟)=(𝑠𝐼−Ω)⁻¹Ω𝐯(𝑠,𝑟−1),𝑟>1, (43) where Ω=−𝑁𝜇₂𝐼_𝐾+1+𝛼𝐶_𝐾+1𝐼_𝐾+1+𝐴_𝐾+1⊗𝐼_{(𝑀+1)(𝑁−𝑀+1)𝑊}+𝜇₂−𝜇₁𝐼_𝐾+1⊗𝐶_𝑀+1⊗𝐼_{𝑁−𝑀+1}+𝐼_𝑀+1⊗𝐶_{𝑁−𝑀+1}⊗𝐼_𝑊+𝐼_{(𝐾+1)(𝑀+1)(𝑁−𝑀+1)}⊗𝐷₀+𝑃+𝐼−𝑃𝐸⁺_𝐾+1⊗𝐼_{(𝑀+1)(𝑁−𝑀+1)}⊗𝐷₁+𝐴_𝐾+1⊗𝜇₁𝐶_𝑀+1+𝜇₂𝐶_𝑀+1𝐸⁺_𝑀+1⊗𝐼_{(𝑁−𝑀+1)𝑊},Ω=𝐼_{(𝐾+1)(𝑀+1)}⊗𝜇₁𝐶_{𝑁−𝑀+1}𝐴_{𝑁−𝑀+1}+𝜇₂𝐶_{𝑁−𝑀+1}⊗𝐼_𝑊+𝐼⊗𝜇₁𝐶_𝑀+1𝐴_𝑀+1+𝜇₂𝐶_𝑀+1⊗𝐼_{(𝑁−𝑀+1)𝑊},𝐚(𝑠)=𝜇₁−𝜇₂𝐞_{(𝐾+1)(𝑀+1)}⊗𝐶_{𝑁−𝑀+1}𝐞⊗𝐞_𝑊+𝜇₁−𝜇₂𝐞̃⊗𝐶_𝑀+1𝐞⊗𝐞_{(𝑁−𝑀+1)𝑊}+𝑁𝜇₂𝐞̃+(𝑁−𝑀)𝜇₂𝐞_𝐾+1−𝐞̃⊗𝐞_{(𝑀+1)(𝑁−𝑀+1)𝑊}𝜇₂𝑠+𝜇₂. (44)
Here, 𝐶_𝑙=diag{0,1,…,𝑙−1}, 𝐶_𝑙=diag{𝑙−1,𝑙−2,…,0}, 𝑃=diag{𝑝₀,𝑝₁,…,𝑝_𝐾}; 𝐴_𝑙 is a square matrix of size 𝑙 with all zero entries except the entries (𝐴_𝑙)_{𝑖,𝑖−1},𝑖=1,𝑙−1, which are equal to 1; 𝐼 is a square matrix of size 𝐾+1 with all zero entries except the entry (𝐼)_0,0=1; 𝐸⁺_𝑙 is a square matrix of size 𝑙 with all zero entries except the entries (𝐸⁺_𝑙)_𝑖,𝑖+1, 𝑖=0,𝑙−1, which are equal to 1; 𝐞̃ is a column vector of size 𝐾+1 with all zero entries except the entry (𝐞̃)₀=1.

Proof. Based on a probabilistic sense of the LST, we can obtain the following system of linear algebraic equations for calculation of the vectors 𝐯(𝑠,𝑟,𝑘,𝑛,𝑙): 𝐯()=()𝜇𝑠,𝑟,𝑘,𝑛,𝑙𝑠+𝑘𝛼+𝑛+𝑙₁+(𝑁−(𝑛+𝑙))𝜇₂𝐼_𝑊−𝐷₀⁻¹×𝛿_𝑟,1𝛿_𝑘,0(𝑛+𝑙)𝜇₁+(𝑁−(𝑛+𝑙))𝜇₂+1−𝛿_𝑘,0𝑙𝜇₁+(𝑁−𝑀−𝑙)𝜇₂×𝜇₂𝑠+𝜇₂𝐞_𝑊+1−𝛿_𝑟,1×𝛿_𝑘,0𝑛𝜇₁𝐯()+()𝜇𝑠,𝑟−1,0,𝑛−1,𝑙𝑀−𝑛₂𝐯(𝑠,𝑟−1,0,𝑛,𝑙)+𝑙𝜇₁𝐯()+()𝜇𝑠,𝑟−1,𝑘,𝑛,𝑙−1𝑁−𝑀−𝑙₂𝐯(+𝑠,𝑟−1,𝑘,𝑛,𝑙)1−𝛿_𝑘,0×𝑘𝛼+𝑛𝜇₁𝐯()+()𝜇𝑠,𝑟,𝑘−1,𝑛,𝑙𝑀−𝑛₂𝐯(+𝑠,𝑟,𝑘−1,𝑛+1,𝑙)1−𝑝_𝑘𝐷₁𝐯(𝑠,𝑟,𝑘+1,𝑛,𝑙)+𝑝_𝑘𝐷₁𝐯(𝑠,𝑟,𝑘,𝑛,𝑙),𝑙=0,𝑁−𝑀,𝑛=0,𝑀,𝑘=0,𝐾,𝑟≥1, (45) where 𝛿_𝑖,𝑗 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 𝑠,𝑟,𝑘,𝑛, and 𝑙 (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 𝜇₂/(𝑠+𝜇₂) 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 𝑠, 𝑟, 𝑘, 𝑛, and 𝑙 do not change.
Using the notation of the matrices given in Theorem 5, let us rewrite system (45) into the matrix form as follows: −𝑠+𝑘𝛼−𝑛𝜇₂−𝜇₁+𝑁𝜇₂𝐼_{(𝑁−𝑀+1)𝑊}+𝜇₂−𝜇₁𝐶_{𝑁−𝑀+1}⊗𝐼_𝑊+𝐼_{𝑁−𝑀+1}⊗𝐷₀+𝑝_𝑘𝐷₁𝐯(𝑠,𝑟,𝑘,𝑛)+𝛿_𝑟,1𝜇₁−𝜇₂𝐶_{𝑁−𝑀+1}𝐞⊗𝐞_𝑊+𝛿_𝑘,0𝑛𝜇₁−𝜇₂+𝑁𝜇₂+1−𝛿_𝑘,0(𝑁−𝑀)𝜇₂𝐞_{(𝑁−𝑀+1)𝑊}×𝜇₂𝑠+𝜇₂+1−𝛿_𝑟,1×𝜇₁𝐶_{𝑁−𝑀+1}𝐴_{𝑁−𝑀+1}+𝜇₂𝐶_{𝑁−𝑀+1}⊗𝐼_𝑊𝐯(𝑠,𝑟−1,𝑘,𝑛)+𝛿_𝑘,0𝑛𝜇₁𝐯()+()𝜇𝑠,𝑟−1,0,𝑛−1𝑀−𝑛₂𝐯()+𝑠,𝑟−1,0,𝑛1−𝛿_𝑘,0𝑘𝛼+𝑛𝜇₁𝐯()+()𝜇𝑠,𝑟,𝑘−1,𝑛𝑀−𝑛₂𝐯(+𝑠,𝑟,𝑘−1,𝑛+1)1−𝑝_𝑘𝐼_{𝑁−𝑀+1}⊗𝐷₁𝐯(𝑠,𝑟,𝑘+1,𝑛)=𝟎^𝑇. (46)
System (46) can be further rewritten into −𝑠+𝑘𝛼+𝑁𝜇₂𝐼_{(𝑀+1)(𝑁−𝑀+1)𝑊}+𝜇₂−𝜇₁𝐶_𝑀+1⊗𝐼_{𝑁−𝑀+1}+𝐼_𝑀+1⊗𝐶_{𝑁−𝑀+1}⊗𝐼_𝑊+𝐼_{(𝑀+1)(𝑁−𝑀+1)}⊗𝐷₀+𝑝_𝑘𝐷₁𝐯(𝑠,𝑟,𝑘)+𝛿_𝑟,1𝜇₁−𝜇₂𝐞_𝑀+1⊗𝐶_{𝑁−𝑀+1}𝐞⊗𝐞_𝑊+𝛿_𝑘,0𝜇₁−𝜇₂𝐶_𝑀+1𝐞⊗𝐞_{(𝑁−𝑀+1)𝑊}+𝑁𝜇₂𝐞_{(𝑀+1)(𝑁−𝑀+1)𝑊}+1−𝛿_𝑘,0𝜇₂(𝑁−𝑀)𝐞_{(𝑀+1)(𝑁−𝑀+1)𝑊}𝜇₂𝑠+𝜇₂+1−𝛿_𝑟,1𝐼_𝑀+1⊗𝜇₁𝐶_{𝑁−𝑀+1}𝐴_{𝑁−𝑀+1}+𝜇₂𝐶_{𝑁−𝑀+1}⊗𝐼_𝑊𝐯(𝑠,𝑟−1,𝑘)+𝛿_𝑘,0𝜇₁𝐶_𝑀+1𝐴_𝑀+1+𝜇₂𝐶_𝑀+1⊗𝐼_{(𝑁−𝑀+1)𝑊}𝐯(+𝑠,𝑟−1,0)1−𝛿_𝑘,0𝑘𝛼𝐼_𝑀+1+𝜇₁𝐶_𝑀+1+𝜇₂𝐶_𝑀+1𝐸⁺_𝑀+1⊗𝐼_{(𝑁+𝑀−1)𝑊}𝐯(𝑠,𝑟,𝑘−1)+1−𝑝_𝑘𝐼_{(𝑀+1)(𝑁−𝑀+1)}⊗𝐷₁𝐯(𝑠,𝑟,𝑘+1)=𝟎^𝑇. (47)
Taking into account notations (44), performing algebra, we obtain the following system: (−𝑠𝐼+Ω)𝐯(𝑠,𝑟)+1−𝛿_𝑟,1Ω𝐯(𝑠,𝑟−1)+𝛿_𝑟,1𝐚(𝑠)=𝟎^𝑇,𝑟≥1. (48)
Based on system (48), it is easy to verify that the vectors 𝐯(𝑠,𝑟),𝑟≥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 LST 𝑣(𝑠) of the distribution of an arbitrary virtual customer’s sojourn time in the system is computed as 𝑣(𝑠)=1𝜆^{𝑣𝑖𝑟𝑡}_𝑜𝑢𝑡^∞_𝑖=𝑁^{min{𝑖−𝑁,𝐾}𝑀}_𝑘=0𝑛=0^𝑁−𝑀_𝑙=0𝑝_𝑘𝝅(𝑖,𝑘,𝑛,𝑙)𝐷₁()+×𝐯𝑠,𝑖−𝑁−𝑘+1,𝑘,𝑛,𝑙_{∞𝑖=𝑁min{𝑖−𝑁,𝐾}𝑀𝑘=1𝑛=0𝑁−𝑀𝑙=0}()(.𝑘𝛼𝝅𝑖,𝑘,𝑛,𝑙×𝐯𝑠,𝑖−𝑁−𝑘+1,𝑘−1,𝑛,𝑙) (49)

Corollary 7. The average sojourn time 𝑉^𝑠𝑜𝑗 of an arbitrary virtual customer is computed by 𝑉^𝑠𝑜𝑗1=−𝜆^{𝑣𝑖𝑟𝑡}_𝑜𝑢𝑡^∞_𝑖=𝑁^{min{𝑖−𝑁,𝐾}𝑀}_𝑘=0𝑛=0^𝑁−𝑀_𝑙=0𝑝_𝑘𝝅(𝑖,𝑘,𝑛,𝑙)𝐷₁×𝐯^′(𝑠,𝑖−𝑁−𝑘+1,𝑘,𝑛,𝑙)|_𝑠=0+^∞_𝑖=𝑁^{min{𝑖−𝑁,𝐾}𝑀}_𝑘=1𝑛=0^𝑁−𝑀_𝑙=0𝑘𝛼𝝅(𝑖,𝑘,𝑛,𝑙)×𝐯^′(𝑠,𝑖−𝑁−𝑘+1,𝑘−1,𝑛,𝑙)|_𝑠=0. (50)
Here, the column vectors 𝐯^′(𝑠,𝑟,𝑘,𝑛,𝑙)|_𝑠=0 are calculated as the blocks of the vector 𝐯^′(𝑠,𝑟)|_𝑠=0,𝑟≥1, which can be computed as follows: 𝐯^′(𝑠,1)_𝑠=0=Ω⁻¹𝐞−𝐚^′(𝑠)_𝑠=0,𝐯^′(𝑠,𝑟)_𝑠=0=Ω⁻¹𝐞−Ω𝐯^′(𝑠,𝑟−1)_𝑠=0,𝑟>1. (51)

Proof. Formula (50) is based on the definition 𝑉^soj=−𝑣^′(𝑠)|_𝑠=0, taking into account that 𝐯(𝑠,𝑟)|_𝑠=0=𝐞,𝑟≥1.

Corollary 8. The average waiting time 𝑉^{𝑤𝑎𝑖𝑡} of an arbitrary virtual customer is defined as 𝑉^{𝑤𝑎𝑖𝑡}=𝑉^𝑠𝑜𝑗−1𝜇₂. (52)

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 matrices 𝐷₀ and 𝐷₁, with the same average arrival rate 𝜆=9 but different coefficients of correlation and variation.
The first process, coded as 𝑀⁰, is a stationary Poisson process. It is defined by the matrices 𝐷₀=(−9) and 𝐷₁=(9) and has the coefficient of correlation 𝑐_cor=0 and the coefficient of variation 𝑐_var=1.
The second process, coded as IPP⁰ (interrupted Poisson process), is defined by the matrices 𝐷₀=,𝐷−28.645400.275840.97606−0.97606₁=.26.544741.8248200 (53) It has the coefficient of correlation 𝑐_cor=0 and the coefficient of variation 𝑐_var=12.4.
The third process, coded as MAP^0.2, defined by the matrices 𝐷₀=,𝐷−12.1647600−0.39483₁=12.083880.080880.219910.17492 (54) and has 𝑐_cor=0.2 and 𝑐_var=12.34.
The fourth process, coded as MAP^0.4, is defined by the matrices 𝐷₀=,𝐷−30.5791800.00913−0.99199₁=.30.260680.318500.109190.87367 (55) It has 𝑐_cor=0.4 and 𝑐_var=12.39.
We assume that the number of servers 𝑁 is equal to 15, the service intensities of real and virtual customers 𝜇₁=0.7 and 𝜇₂=0.8, respectively, the buffer capacity 𝐾=15, the intensity of impatience 𝛼=0.3, and the probabilities 𝑝_𝑘=0.02(𝑘+1),𝑘=0,14.
Let us vary the number of type 1 servers 𝑀 in the interval [0,15].

Figure 2 illustrates the dependence of the average number of customers in the system 𝐿, the average waiting time of a virtual customer 𝑉^wait, the average number of real customers 𝑁^real_buf, and the average number of virtual customers 𝑁^virt_buf in the buffer on the number of type 1 servers 𝑀 for the four MAP arrival processes presented above.

(a)

(b)

(c)

(d)

Figure 2 

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.

The dependence of the intensity of flow of served real customers 𝜆^real_out, the intensity of flow of served virtual customers 𝜆^virt_out, the probability 𝑃_ent that an arbitrary real customer becomes a virtual customer at the entrance to the system due to the presence of 𝐾 real customers in the buffer and the probability 𝑃_imp 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 𝑀 for different MAP arrival processes is presented in Figure 3.

(a)

(b)

(c)

(d)

Figure 3 

Intensities of output flow of real and virtual customers, who received service in the system, and the loss probabilities 𝑃ent and 𝑃imp as functions of the number of type 1 servers 𝑀 for different MAPs.

From Figures 2 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 𝜆^virt_out and 𝜆^real_out decrease essentially, while the average number of customers in the system 𝐿, the average waiting time of an arbitrary virtual customer 𝑉_wait, and the loss probability 𝑃_ent 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 servers 𝑀 the intensity 𝜆^real_out increases, while intensity 𝜆^virt_out decreases. However, at the same time, the average waiting time of an arbitrary virtual customer 𝑉_wait increases essentially with increasing 𝑀, 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, 𝑉^wait=𝑁^buffer_virt𝜆𝑃_to-virt=𝑁^buffer_virt𝜆^virt_out. (56)

Additionally, the numerical results show that another equality holds true: 𝑁^real_buf/𝜆=𝑃_imp/𝛼. To prove this equality, let us rewrite it in the form 𝛼𝑁^real_buf=𝜆𝑃_imp. 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 servers 𝑀. The aim of optimization is the maximization of a cost criterion (one of the possible criteria): Λ(𝑀)=𝜆^real_out+1−𝑞₁𝑞₂𝜆^virt_out. (57) Here, 𝑞₁ is the probability that an operator will not succeed to connect to a customer (customer's phone is busy or does not answer) and 𝑞₂ is the probability that a virtual customer will be interested in the services proposed by an operator. It is evident that the probability 𝑞₂ depends on the waiting time. Thus, Λ(𝑀) determines the effective intensity of the flow of customers, who successfully receive service in the system.
Let us fix 𝑞₁=0.1, and 𝑞₂=0,𝑉^wait𝑉≥40,1−^wait40,𝑉^wait<40. (58)
Our goal is to find the optimal value 𝑀^∗ that provides the maximum value of the cost criterion Λ(𝑀). The value of the cost criterion Λ(𝑀) for different values of the number of type 1 servers 𝑀 and different arrival processes is presented in Table 1. The optimal values of Λ(𝑀) for each arrival process are printed in bold type. Additionally, Table 1 contains the values of the cost criterion Λ(𝑀) for the call center that does not provide a call-back option. In this case, Λ(𝑀)=𝜆_out and can be found based on the results of the paper [11].
As it is seen from Table 1, for the stationary Poisson arrival process, the maximum value of Λ(𝑀) is equal to 8.91163 for 𝑀=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 Λ(𝑀) is equal to 8.32646, when the system has seven type 1 operators and eight type 2 operators. For MAP^0.2 arrival process, the maximum value of Λ(𝑀) is equal to 8.39144, when the system has three type 1 operators and twelve type 2 operators and for MAP^0.4 arrival process the maximum value of Λ(𝑀) is equal to 5.84454 when the system has only type 2 operators.
From Table 1, 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 Λ(𝑀) 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 Λ(𝑀) 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 servers 𝑀 and different arrival processes.

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.

Acknowledgment

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant no. 2011-0015214).