# Analysis of Multiserver Retrial Queueing System with Varying Capacity and Parameters

Authors: Alexander N. Dudin; Olga S. Dudina
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A multiserver queueing system,  the dynamics of which depends on the state of some external continuous-time
Markov chain (random environment, RE), is considered. Change of the state of the RE may cause variation of
the parameters of the arrival process, the service process, the number of available servers, and the available buffer
capacity, as well as the behavior of customers. Evolution of the system states is described by the multidimensional
continuous-time Markov chain. The generator of this Markov chain is derived. The ergodicity condition is presented.
Expressions for the key performance measures are given. Numerical results illustrating the behavior of the system
and showing possibility of formulation and solution of optimization problems are provided. The importance of the
account of correlation in the arrival processes is numerically illustrated.


---

## Body

## 1. Introduction

Queueing theory is widely used for decision making about the resources needed to provide service in a variety of real life systems including contact centers, intelligent transportation systems, telecommunication networks, manufacturing and administrative systems, and banking. Due to the increasing technical and logical complexity of real life systems and diversity of the provided services, queueing models of these systems and their fragments become more and more involved. One of the essential features of modern systems, which should be accounted in queueing models, is that the parameters of the system operation may vary, for example, depending on time of a day or night, available amount of required resources, and possible sharing with some other systems. This fact gave rise to the progress in the study of the so-called queues operating in the random environment (R
E).In contrast to the classical queueing models, where the parameters and the distributions characterizing arrival, service, and other processes describing dynamics of the system are assumed to be fixed in advance, queues operating in theR
E presuppose that some or all parameters may dynamically vary due to influence of some external medium called R
E
. Early research in this topic was done by Gnedenko and Kovalenko [1], Yechiali and Naor [2], Yadin and Syski [3], O’Cinneide and Purdue [4], Purdue [5], Neuts [6, 7], and others. A brief history of the development of theory of queues in the R
E, the reference list, and real life examples of queues in the R
E can be found, for example, in the papers [8–12]. In [8], an unreliable M
/
M
/
1 retrial queue in a Markovian random environment is analyzed via matrix-analytic methods. Ergodicity condition is proved and approximate distribution of the number of customers in the system is computed. Optimization problem of choosing the arrival and service rates for each environment state is considered. In [9, 10], the finite source M
A
P
/
P
H
/
N retrial queue operating in a random environment is studied. The arrival flow is described by the Markov arrival process (M
A
P) (for definition, more details, and usefulness in modelling in telecommunications, see [13–15]), and the service time has a phase-type (P
H) distribution (for definition and more details, see [7]). It is assumed that the parameters defining the M
A
P process and P
H distribution depend on the state of the R
E that is a continuous-time Markov chain with a finite state space. In [9], it is assumed that there is additional M
A
P arrival process of negative customers. The arrival of the negative customer with equal probability goes to any busy server to remove the customer being in service. In [9, 10], the finite state multidimensional Markov chain describing the behavior of the systems is investigated. The algorithms for calculating the stationary state probabilities are elaborated. Main performance measures are obtained and the illustrative numerical examples are presented. In [11], the B
M
A
P
/
P
H
/
N
/
N queue operating in the R
E is investigated. The arrival flow is described by the batch Markov arrival process (B
M
A
P). The system does not have a buffer. An arriving customer who did not succeed to find a free server upon arrival is lost. Due to possibility of batch arrivals, disciplines of partial admission, complete admission, and complete rejection are analyzed. The stationary distribution of the system states and the waiting time distribution are computed. Numerical illustrations are presented. In particular, it is demonstrated that reasonable engineering approximations of performance measures of the system may be very poor. In [12], the B
M
A
P
/
P
H
/
N retrial queue operating in the R
E is investigated. As in the previously described model, the system does not have a buffer. But the customer who did not succeed to find a free server upon arrival is not lost, but he proceeds to the retrial orbit, which is a location from which customers may attempt to gain or regain service. The intensities of repeating the attempts are also assumed to be dependent on the state of the R
E
. In [12], the sufficient condition of stability of the system is proved, the stationary distribution of the system states is computed. The presented numerical examples illustrate a poor quality of the approximation of the main performance measures of the system by means of the simpler queueing models. An effect of possible smoothing the traffic and an impact of the retrial intensity are shown. Analysis presented in [12] is more complicated than the analyses in [9, 10] due to two evident reasons: (i) the state space of the multidimensional Markov chain describing the behavior of the system is finite in [9, 10] and is infinite in [12]; (ii) the generator of the multidimensional Markov chain has tridiagonal structure in [9, 10] (i.e., the chain is level dependent quasi-birth-and-death process) while it has more general upper-Hessenbergian structure in [12].The main advantages of the model considered in this paper, compared to the ones analyzed in [11, 12] and to all the other papers devoted to multiserver queues operating in the R
E, are as follows.(i)
The model under study has more flexible service discipline that combines the features of loss systems, systems with a finite buffer and systems with retrials. The queueing systems that incorporate both normal queues and retrial orbits are called hybrid retrial queues. The importance of their investigation stems from the fact that many modern technologies of customers random access assume the existence of some places where the customers, who did not succeed to get access upon arrival, may be temporarily kept (in registers for handover customers in cells of mobile communication networks, in IVR (Interactive Voice Response) machines in call centers, etc.). Such hybrid systems were considered, for example, in [16, 17]. Queueing systems in R
E with customers loss considered in [11] and with retrials considered in [12] are obtained as special cases of the hybrid retrial queue in R
E considered in this paper.
(ii)
We suppose that the capacity of the service area of the system (the number of servers and places in the buffer) may depend on the state of theR
E, while it is usually assumed that the capacity is constant and the R
E influences only the parameters of the arrival and service processes. Note that the changes in service capacity may be considered as server breakdowns, among other possibilities.
(iii)
Customers balking, impatience (abandonment), and nonpersistence are taken into account. This allows us to account psychology of customers, use of visible queue option, information obsolescence during the waiting time, customer’s mobility, service provider’s competition, and so forth.Due to the generality of the considered queueing model, it has a lot of possible applications for the investigation of a variety of real life systems. Besides the examples of application to the performance evaluation, capacity planning, and optimization of some systems, which coincide with the ones described in [12] (application to modelling hot spot in airport and wireless local area network), it is worth to mention important potential applications for work force management in call-centers. Fluctuation of the intensities of customer arrivals and retrials occurs here, for example, due to the well known existence of hours of low, middle, and peak load of the center. The servers correspond to the operators of the call-center and places in the buffer correspond to IVR machines. Impatience of customers staying in the buffer reflects the possibility to get the requested information directly from IVR without contacting to an operator. It is necessary to create the schedule of the work of the operators and IVRs to dynamically fit the number of the active operators and IVRs to the current level of the load of the call-center in such a way as to guarantee the best quality of the customer’s service under the existing restrictions imposed on the total number of available operators and their working schedule and on energy consumption by IVRs.The rest of the paper is organized as follows. In Section2, the mathematical model is described. The process of the system states is defined; its generator in a block matrix form is presented. In Section 3, the stability condition is derived and the problem of computation of the stationary probabilities of the system states is touched on. Formulas for computation of the steady state performance measures of the system are presented in Section 4. Numerical illustrations are given and briefly discussed in Section 5. Section 6 concludes the paper.


## 2. Mathematical Model

We consider a retrial multiserver queueing system with varying capacity and behavior of the customers. The structure of the system under study is presented in Figure1.Figure 1
Queueing system under study.The dynamics of the system depends on the state of theR
E. The R
E is defined by the stochastic process r
t,  t
≥
0, which is an irreducible continuous-time Markov chain with the state space {
1,2
,
…
,
R
} and the infinitesimal generator H
.Under the fixed stater of the R
E, capacity of the system is equal to K
(
r
) including N
(
r
) servers, 0
≤
N
(
r
)
≤
K
(
r
), and a waiting room (buffer) of size K
(
r
)
-
N
(
r
), r
∈
{
1
,
…
,
R
}
. Without the loss of generality, we assume that the states of the R
E are enumerated in ascending order of the capacity of the system; that is, 0
≤
K
(
1
)
≤
K
(
2
)
≤
⋯
≤
K
(
R
)
. We also suggest that 0
≤
N
(
1
)
≤
N
(
2
)
≤
⋯
≤
N
(
R
)
. We call the servers and the buffer as the service area of the system.Arrival of customers is defined by the Markovian arrival process (M
A
P). The underlying process of this M
A
P is {
r
t
,
ν
t
},  t
≥
0, where r
t is the state of the R
E and the process ν
t has a finite state space {
0,1
,
…
,
W
}. Under the fixed state r of the R
E the process ν
t behaves as irreducible continuous-time Markov chain. The sojourn time of this chain in the state ν is exponentially distributed with the positive finite parameter λ
ν
(
r
)
. When the sojourn time in the state ν expires, with probability p
0
r
ν
,
ν
′, the process ν
t jumps to the state ν
′ without generation of a customer, ν
,
ν
′
∈
0
,
…
,
W,  ν
≠
ν
′, r
∈
1
,
…
,
R
. With probability p
1
r
ν
,
ν
′, the process ν
t jumps to the state ν
′ with generation of a customer, ν
,
ν
′
∈
0
,
…
,
W, r
∈
1
,
…
,
R
.The behavior of the arrival process under the fixed stater of the R
E is completely characterized by the matrices D
0
(
r
) and D
1
(
r
) defined by the entries (1)
D
0
r
ν
,
ν
=
-
λ
ν
r
,
ν
∈
0
,
…
,
W
,
D
0
r
ν
,
ν
′
=
λ
ν
r
p
0
r
ν
,
ν
′
,
ν
,
ν
′
∈
0
,
…
,
W
,
ν
≠
ν
′
,
D
1
r
ν
,
ν
′
=
λ
ν
r
p
1
r
ν
,
ν
′
,
ν
,
ν
′
∈
0
,
…
,
W
,
r
∈
1
,
…
,
R
.The square matrixD
(
r
)
(
1
)
=
D
0
(
r
)
+
D
1
(
r
) of dimension W
¯
=
W
+
1 represents the generator of the process ν
t,  t
≥
0, under the fixed state r, r
∈
{
1
,
…
,
R
}.The average arrival rateλ
(
r
) under the fixed state r of the R
E is given as (2)
λ
r
=
θ
r
D
1
r
e
,where θ
(
r
) is the invariant vector of the stationary distribution of the Markov chain ν
t, t
≥
0, under the fixed state r. The vector θ
(
r
) is the unique solution to the system (3)
θ
r
D
r
1
=
0
,
θ
r
e
=
1
.The squared coefficient of variationc
v
a
r
(
r
) of intervals between successive arrivals under the fixed state r of the R
E is given as (4)
c
v
a
r
r
=
2
λ
r
θ
r
-
D
0
r
-
1
e
-
1
,
r
∈
1
,
…
,
R
.The coefficient of correlationc
c
o
r
(
r
) of two successive intervals between arrivals under the fixed state r of the R
E is given as (5)
c
cor
r
=
λ
r
θ
r
-
D
0
r
-
1
D
r
1
-
D
0
r
-
D
0
r
-
1
e
-
1
c
v
a
r
r
,
r
∈
1
,
…
,
R
.Let us introduce the following matrices:(6)
D
~
1
=
diag
⁡
D
1
r
,
r
∈
1
,
…
,
R
,
D
~
0
=
H
⊗
I
W
¯
+
diag
⁡
D
0
r
,
r
∈
1
,
…
,
R
,
D
~
1
=
D
~
0
+
D
~
1
.The averaged (over all the states of theR
E) intensity λ of input flow of customers is defined as (7)
λ
=
θ
D
~
1
e
,where the vector θ is the unique solution of the system (8)
θ
D
~
1
=
0
,
θ
e
=
1
.The squared coefficient of variation c
v
a
r of intervals between successive arrivals is given as (9)
c
v
a
r
=
2
λ
θ
-
D
~
0
-
1
e
-
1
.The coefficient of correlation c
c
o
r of two successive intervals between arrivals is given as (10)
c
cor
=
λ
θ
-
D
~
0
-
1
D
~
1
-
D
~
0
-
1
e
-
1
c
v
a
r
.We assume that, during the epochs of the transitions of processr
t, t
≥
0, the states of the process ν
t, t
≥
0, do not change, and only the intensities of transitions of this process change.If, at an arbitrary arrival moment, some server is free, an arriving customer starts the service. If all servers are busy, but the buffer is not full, the customer is placed to the buffer. If the buffer is full and the state of theR
E is r, the customer balks (leaves the system permanently) with probability 1
-
p
1
(
r
) or, with the complementary probability, moves to orbit and tries again later on.If theR
E transits from the state r to the state r
′, where r
′
<
r, the capacity of the system (the number of available servers and/or places in the buffer) decreases. If the current number of customers in service and/or in the buffer exceeds the available number of servers and/or places in the buffer under the state r
′ of the R
E, the redundant customers leave the service area. The customers leaving the service area go, independently of each other, into orbit with probability p
2
(
r
) or leave the system permanently with the complementary probability.When the state of theR
E is r, r
∈
1
,
…
,
R, each customer staying in orbit repeats the attempts to reach the service area after an exponentially distributed time described by the parameter α
(
r
), α
(
r
)
≥
0
. Customers retry independent of each other and are not absolutely persistent. If the attempt is successful, the customer leaves orbit and moves to service or to the buffer. Otherwise, the customer returns to orbit with probability p
3
(
r
)
. With the complementary probability, the customer leaves the system permanently.Customers in orbit are impatient. When the state of theR
E is r, r
∈
{
1
,
…
,
R
}, each customer may leave the orbit and the system after an exponentially distributed time described by the parameter γ
(
r
),  γ
(
r
)
>
0
. If the customers staying in the orbit are patient, we put γ
(
r
)
=
0
.We assume that the probability of the service completion of a customer during the time interval(
t
,
t
+
Δ
t
) is equal to μ
(
r
)
Δ
t
+
o
(
t
) when the state of the R
E is r, r
∈
{
1
,
…
,
R
}
.If theR
E transits from the state r to the state r
′, where r
′
>
r, the capacity of the system increases and the corresponding number of customers from the buffer, if any, occupies additional N
(
r
′
)
-
N
(
r
) servers.Mention that in the described queueing model the discipline of customers admission may be different under the various states of theR
E
. For example, if p
1
(
r
)
=
1 and α
(
r
)
=
0, no customers are admitted to the system if the buffer is full; that is, the system operates as the system with a finite buffer. If, additionally, K
(
r
)
=
N
(
r
), the system operates as the system with customers loss (Erlang loss model). If p
1
(
r
)
=
0, the system operates as the usual retrial system where all arriving customers, who meet the service area full, go into orbit. In the presented analysis, it is accounted that N
(
r
) can be equal to 0, that is, for some states of the R
E service to customers are not provided at all.In order to improve the readability of the paper, we collect and already introduced some new notation in Notation.Consider the following:(i)
leti
t, i
t
≥
0, be the number of customers in orbit,
(ii)
letr
t, r
t
∈
{
1
,
…
,
R
}, be the state of the R
E,
(iii)
letn
t, n
t
∈
{
0
,
…
,
K
(
r
)
}, be the number of customers in the service area,
(iv)
letν
t, ν
t
∈
{
0
,
…
,
W
}, be the state of the second component of the underlying process of customers arrivals
at the momentt, t
≥
0
.It is easy to see that the processξ
t
=
{
i
t
,
r
t
,
n
t
,
ν
t
}, t
≥
0, is the four-dimensional irreducible Markov chain. Let us enumerate the states of the Markov chain ξ
t, t
≥
0, in the direct lexicographic order of the components (
i
,
r
,
n
,
ν
). We call the set of the states having value (
i
,
r
) of the two first components of the Markov chain the macrostate (
i
,
r
).LetQ be the generator of the Markov chain ξ
t, t
≥
0
. It is formed by the blocks Q
i
,
j, consisting of the matrices (
Q
i
,
j
)
r
,
r
′ that define the intensities of the transitions of the Markov chain ξ
t, t
≥
0, from the macrostate (
i
,
r
) to the macrostate (
j
,
r
′
), r
,
r
′
∈
{
1
,
…
,
R
}
. The diagonal entries of the matrix Q
i
,
i are negative. The modulus of each entry defines the intensity of departing from the corresponding state of the Markov chain ξ
t,  t
≥
0.Lemma 1.
The generatorQ has the following block structure: (11)
Q
=
Q
0,0
Q
0,1
Q
0,2
…
Q
0
,
K
¯
O
O
⋯
Q
1,0
Q
1,1
Q
1,2
…
Q
1
,
K
¯
Q
1
,
K
¯
+
1
O
⋯
O
Q
2,1
Q
2,2
…
Q
2
,
K
¯
Q
2
,
K
¯
+
1
Q
2
,
K
¯
+
2
⋯
⋮
⋮
⋮
⋱
⋮
⋮
⋮
⋱
,where K
¯
=
max
⁡
{
max
⁡
{
K
(
R
)
-
K
(
1
)
,
N
(
R
)
-
N
(
1
)
}
,
1
}
.
The nonzero blocksQ
i
,
j,  i
,
j
≥
0, are defined as follows. Consider(i)
(12)
Q
i
,
i
=
Q
i
,
i
r
,
r
′
,
r
,
r
′
∈
1
,
…
,
R
,where (13)
Q
i
,
i
r
,
r
=
-
μ
r
C
r
+
i
α
r
+
γ
r
I
K
r
+
1
+
μ
r
C
r
E
K
r
-
+
i
p
3
r
α
r
I
^
K
r
⊗
I
W
¯
+
1
-
p
1
r
I
^
K
r
+
E
K
r
+
⊗
D
1
r
+
I
K
r
+
1
⊗
D
0
r
+
H
r
,
r
I
K
r
+
1
W
¯
,
(14)
Q
i
,
i
r
,
r
′
=
H
r
,
r
′
M
r
,
r
′
0
⊗
I
W
¯
,
r
′
<
r
,
(15)
Q
i
,
i
r
,
r
′
=
H
r
,
r
′
M
r
,
r
′
+
⊗
I
W
¯
,
r
′
>
r
,
i
≥
0
.Hereinafter, the entry (
M
r
,
r
′
n
)
k
,
k
′ of the matrix M
r
,
r
′
n, r
∈
1
,
…
,
R,  r
′
∈
1
,
…
,
r
-
1, of dimension K
r
+
1
×
K
r
′
+
1 defines the probability that when the state of the RE changes from r to r
′ (and the capacity of the service area decreases from K
(
r
) to K
r
′) the number of customers in the service area changes from k to k
′ and n customers move to orbit. The matrices M
r
,
r
′
n have the following nonzero entries. Consider (16)
M
r
,
r
′
0
k
,
k
=
1
,
k
∈
0
,
…
,
N
r
′
,
M
r
,
r
′
n
k
,
N
r
′
=
p
r
n
,
k
-
N
r
′
,
k
∈
N
r
′
+
1
,
…
,
N
r
.
M
r
,
r
′
n
k
,
N
r
′
+
k
-
N
r
=
p
r
n
,
N
R
-
N
r
′
,
k
∈
N
r
+
1
,
…
,
min
⁡
K
r
′
+
N
r
-
N
r
′
,
K
r
⁡
,
M
r
,
r
′
n
k
,
K
r
′
=
p
r
n
,
k
-
K
r
′
,
k
∈
min
⁡
K
r
′
+
N
r
-
N
r
′
,
K
r
+
1
,
…
,
K
r
.The probability p
(
r
)
(
n
,
k
) that n customers go into orbit in the case when k customers leave the service area is given by (17)
p
r
n
,
k
=
C
k
n
1
-
p
2
r
n
p
2
r
k
-
n
,
n
≤
k
,
0
,
n
>
k
.The matrix M
r
,
r
′
+,  r
∈
{
1
,
…
,
R
-
1
},  r
′
∈
{
r
+
1
,
…
,
R
}, of dimension K
r
+
1
×
K
r
′
+
1, which defines the transition probabilities of the component n
t when the state of the RE changes from r to r
′ (and the capacity of the service area increases from K
(
r
) to K
r
′), is equal to the matrix I
K
r
+
1
∣
O
. Consider(ii)
(18)
Q
i
,
i
+
k
=
Z
1,1
k
O
O
…
O
O
Z
2,1
k
Z
2,2
k
O
…
O
O
⋮
⋮
⋱
⋱
⋮
⋮
Z
R
-
1,1
k
Z
R
-
1,2
k
Z
R
-
1,3
k
⋱
Z
R
-
1
,
R
-
1
k
O
Z
R
,
1
k
Z
R
,
2
k
Z
R
,
3
k
…
Z
R
,
R
-
1
k
Z
R
,
R
k
,
i
≥
0
,
k
∈
1
,
…
,
K
¯
,where the entries of the matrices Z
r
,
r
′
k define the intensities of the transitions of components {
r
t
,
n
t
,
ν
t
} of the Markov chain ξ
t leading to the arrival of k customers into the orbit. These matrices are defined by the following formulas: (19)
Z
r
,
r
1
=
p
1
r
I
^
K
r
⊗
D
1
r
,
Z
r
,
r
k
=
O
,
k
>
1
,
Z
r
,
r
′
k
=
H
r
,
r
′
M
r
,
r
′
k
⊗
I
W
¯
,
k
≥
1
,
r
′
<
r
,
r
,
r
′
∈
1
,
…
,
R
.Consider(iii)
(20)
Q
i
,
i
-
1
=
diag
⁡
i
α
r
E
K
r
+
+
1
-
p
3
r
I
^
K
r
+
γ
r
⊗
I
K
r
+
1
⊗
I
W
¯
,
r
∈
1
,
…
,
R
,
i
≥
1
.Proof of the lemma is performed by means of the analysis of the intensities of all possible transitions of the Markov chainξ
t during the time interval having infinitesimal length. The generator Q is block upper-Hessenbergian matrix; that is, all the blocks below the first subdiagonal are equal to zero blocks. This stems from the fact that, at any moment, no more than one customer may leave the orbit, so the blocks Q
i
,
j are equal to zero when j
<
i
-
1
. The number K
¯ defines the maximal number of customers that can simultaneously move to the orbit at an arbitrary moment. We recall that one customer may arrive to the orbit if the buffer is full at the customer’s arrival moment and at most max
⁡
{
K
(
R
)
-
K
(
1
)
,
N
(
R
)
-
N
(
1
)
} customers may arrive to the orbit if the R
E jumps from state R to state 1 when all servers are busy and the buffer is full. So, the blocks Q
i
,
j are equal to zero when j
>
K
¯
. This explains structure (11) of generator Q
.The negative diagonal entries of the matrix(
Q
i
,
i
)
r
,
r define, up to the sign, intensities of the exit of the Markov chain ξ
t from the macrostate (
i
,
r
)
. This exit can happen due to a customer service completion, attempt of a customer from the orbit to reach the service area, departure of a customer from the orbit due to impatience, change of the state of the R
E, and change of the state of the underlying process of M
A
P arrivals. The nondiagonal entries of the matrix (
Q
i
,
i
)
r
,
r define intensities of transition of the Markov chain ξ
t inside of the macrostate (
i
,
r
)
. These transitions can happen due to the customer service completion, change of the state of the R
E, and change of the state of the underlying process of M
A
P arrivals. This explains form (13) of the block (
Q
i
,
i
)
r
,
r
.Transition from the macrostate(
i
,
r
) to the macrostate i
,
r
′,  r
′
<
r, can happen when the state of the R
E changes from r to r
′,  r
′
<
r
. Such a change of the state of the R
E implies the possible reduction of the number of available servers and places in the buffer. This may cause departure of n customers from the service to the orbit. Matrices M
r
,
r
′
n,  r
∈
1
,
…
,
R,  r
′
∈
1
,
…
,
r
-
1, define the probability that when the state of the R
E changes from r to r
′ (and the capacity of the service area decreases from K
(
r
) to K
r
′) the number of customers in the service area changes from k to k
′ and n customers move to the orbit. This explains formula (14). Careful analysis of probabilities, which form the matrices M
r
,
r
′
n, r
∈
1
,
…
,
R, r
′
∈
1
,
…
,
r
-
1, leads to formulas (16).Transition from the macrostate(
i
,
r
) to the macrostate i
,
r
′, r
′
>
r, can happen when the state of the R
E changes from r to r
′,  r
′
>
r
. Such a change of the state of the R
E implies the possible increase of the number of available servers and places in the buffer. This explains formula (15).The increase of the value of the first components of the Markov chainξ
t from i to i
+
k,  k
∈
{
1
,
…
,
K
¯
}, may happen if one customer arrives to the orbit from outside (for k
=
1) or k customers leave the service area due to the reduction of the available space there. This explains formulas (18)-(19).Finally, the decrease of the value of the first components of the Markov chainξ
t from i to i
-
1 may happen if a customer makes an attempt from the orbit when the service area is not full or a customer leaves the orbit due to impatience. This explains formula (20). Lemma is proved.Remark 2.
It is easy to verify that the following limits exist:(21)
Y
0
=
lim
i
→
∞
⁡
R
i
-
1
Q
i
,
i
-
1
,
Y
1
=
lim
i
→
∞
⁡
R
i
-
1
Q
i
,
i
+
I
,
Y
k
=
lim
i
→
∞
⁡
R
i
-
1
Q
i
,
i
+
k
-
1
,
k
∈
2
,
…
,
K
¯
+
1
,where R
i is a diagonal matrix with the diagonal entries which are defined as the moduli of the corresponding diagonal entries of the matrix Q
i
,
i, i
≥
0
. The matrices Y
k, k
∈
{
0
,
…
,
K
¯
+
1
}, have the following form: (22)
Y
0
=
diag
⁡
Ω
~
1
,
…
,
Ω
~
R
,
Y
1
=
Q
~
1,1
Q
~
1,2
…
Q
~
1
,
R
Q
~
2,1
Q
~
2,2
…
Q
~
1
,
R
⋮
⋮
⋱
⋮
Q
~
R
,
1
Q
~
R
,
2
…
Q
~
R
,
R
,
Y
k
=
Z
~
1,1
k
-
1
O
O
…
O
O
Z
~
2,1
k
-
1
Z
~
2,2
k
-
1
O
…
O
O
⋮
⋮
⋱
⋱
⋮
⋮
Z
~
R
-
1,1
k
-
1
Z
~
R
-
1,2
k
-
1
Z
~
R
-
1,3
k
-
1
⋱
Z
~
R
-
1
,
R
-
1
k
-
1
O
Z
~
R
,
1
k
-
1
Z
~
R
,
2
k
-
1
Z
~
R
,
3
k
-
1
…
Z
~
R
,
R
-
1
k
-
1
Z
~
R
,
R
k
-
1
,
k
>
1
,where (23)
Ω
~
r
=
E
K
r
+
⊗
I
W
¯
,
if
p
3
r
=
1
,
α
r
>
0
,
γ
r
=
0
,
Ω
~
r
=
γ
r
γ
r
+
α
r
I
K
r
+
1
-
I
^
K
r
+
α
r
γ
r
+
α
r
E
K
r
+
+
I
^
K
r
⊗
I
W
¯
if
p
3
r
≠
1
,
α
r
>
0
or
γ
r
≠
0
,
Ω
~
r
=
O
,
if
α
r
=
0
,
γ
r
=
0
,
r
∈
1
,
…
,
R
,
Q
~
r
,
r
′
=
R
1
r
Q
0,0
r
,
r
′
+
δ
r
-
r
′
,
0
I
^
K
r
⊗
I
W
¯
,
if
p
3
r
=
1
,
α
r
>
0
,
γ
r
=
0
,
R
2
r
Q
0,0
r
,
r
′
+
δ
r
-
r
′
,
0
I
K
r
⊗
I
W
¯
,
if
α
r
=
0
,
γ
r
=
0
,
O
,
if
p
3
r
≠
1
,
α
r
>
0
or
γ
r
≠
0
,
r
,
r
′
∈
1
,
…
,
R
,
δ
i
,
j indicates the Kronecker delta, (24)
R
1
r
=
I
^
K
r
⊗
μ
r
N
r
-
H
r
,
r
I
W
¯
+
Σ
0
r
-
1
-
p
1
r
Σ
1
r
-
1
,
r
∈
1
,
…
,
R
,
R
2
r
=
μ
r
C
r
I
W
¯
+
I
K
r
⊗
Σ
0
r
-
H
r
,
r
I
W
¯
-
I
^
K
r
⊗
1
-
p
1
r
Σ
1
r
-
1
,
r
∈
1
,
…
,
R
,
Z
~
r
,
r
′
k
=
R
1
r
Z
r
,
r
′
k
,
if
p
3
r
=
1
,
α
r
>
0
,
γ
r
=
0
,
R
2
r
Z
r
,
r
′
k
,
if
α
r
=
0
,
γ
r
=
0
,
O
,
if
p
3
r
≠
1
,
α
r
>
0
or
γ
r
≠
0
,
r
,
r
′
∈
1
,
…
,
R
.Here, Σ
0
(
r
) and Σ
1
(
r
) are the diagonal matrices, the diagonal entries of which are defined as the corresponding diagonal entries of the matrices -
D
0
(
r
) and D
1
(
r
), respectively.
According to the definition given in [18], the existence of the limits Y
k, k
∈
{
0
,
…
,
K
¯
+
1
}, means that the Markov chain ξ
t,  t
≥
0, belongs to the class of continuous-time asymptotically quasi-Toeplitz Markov chains (A
Q
T
M
C). This fact allows us to use the results of [18] for derivation of the sufficient condition of the ergodicity of the Markov chain ξ
t and computation of its steady state distribution.


## 3. System Stability and Stationary Distribution

As follows from [18], the sufficient condition for the ergodicity of A
Q
T
M
C is the fulfillment of the inequality (25)
y
Y
0
e
>
y
∑
k
=
2
K
¯
+
1
k
-
1
Y
k
e
,where the vector y is the unique solution to the system (26)
y
∑
k
=
0
K
¯
+
1
Y
k
=
y
,
y
e
=
1
.Thus, to check whether or not the Markov chainξ
t is ergodic, it is necessary to substitute the matrices Y
k,  k
∈
{
0
,
…
,
K
¯
+
1
}, to system (26), solve this system, and verify the fulfillment of inequality (25). If this inequality is fulfilled, the Markov chain under study is ergodic. The solution of a finite system of (26) on computer does not meet any essential problems.Mention that there exists an important particular case (the customers are impatient or nonpersistent at least for one state of theR
E) when ergodicity of the Markov chain ξ
t can be established more easy than via solution of system (26) and verification of inequality (25). For this case, the following statement is true.Theorem 3.
If customers are impatient or nonpersistent (γ
(
r
)
≠
0 or p
3
(
r
)
≠
1) at least for one state r of the RE, then the Markov chain ξ
t is ergodic for any set of parameters of the queueing system under study.Proof.
LetL
=
{
l
1
,
l
2
,
…
,
l
S
} be the set of the states of the R
E for which at least one of inequalities p
3
(
r
)
≠
1 or γ
(
r
)
≠
0, is fulfilled, r
∈
L
. It can be verified that, in this case, by means of coordinated perturbations of the block rows and block columns, the matrix Y
=
∑
k
=
0
K
¯
+
1
Y
k can be transformed to the canonical normal form (27)
Y
=
Y
1,1
Y
1,2
O
Y
2,2
.Here, Y
1,1 is the matrix obtained from the matrix Y by means of removing the block rows and block columns with numbers l,  l
∈
L
,
Y
1,2 is the matrix obtained from the matrix Y by means of removing the block rows with numbers l, l
∈
L, and the block columns with numbers r,  r
∈
{
1
,
…
,
R
}
∖
L, and Y
2,2
=
diag
⁡
{
Ω
~
l
,
l
∈
L
}
. This means that the matrix Y is reducible. So, it follows, from Theorem  2 in [18], that the sufficient condition for the ergodicity of the Markov chain ξ
t can be rewritten in the form (28)
z
Y
¯
0
e
>
z
∑
k
=
2
K
¯
+
1
k
-
1
Y
¯
k
e
,where the vector z is the unique solution to the system (29)
z
Y
2,2
=
z
,
z
e
=
1
,and the matrices Y
¯
0 and Y
¯
k are obtained from the matrices Y
0 and Y
k by means of removing the block rows and block columns with numbers r,  r
∈
{
1
,
…
,
R
}
∖
L
. It is easy to check that Y
¯
k
=
O,  k
∈
{
1
,
…
,
K
¯
+
1
} and Y
¯
0 is the stochastic matrix. Consequently, (30)
z
Y
¯
0
e
=
z
e
=
1
>
0
=
z
∑
k
=
2
K
¯
+
1
k
-
1
Y
¯
k
e
.
Theorem is proved.In the sequel, we assume that the ergodicity condition is fulfilled. Then, the following stationary probabilities exist:(31)
π
i
,
r
,
n
,
ν
=
lim
t
→
∞
⁡
P
i
t
=
i
,
r
t
=
r
,
n
t
=
n
,
ν
t
=
ν
,
i
≥
0
,
r
∈
1
,
…
,
R
,
n
∈
0
,
…
,
K
r
,
ν
∈
0
,
…
,
W
.Let us form the row-vectorsπ
i as follows: (32)
π
i
,
r
,
n
=
π
i
,
r
,
n
,
0
,
π
i
,
r
,
n
,
1
,
…
,
π
i
,
r
,
n
,
W
,
n
∈
0
,
…
,
K
r
,
π
i
,
r
=
π
i
,
r
,
0
,
π
i
,
r
,
1
,
…
,
π
i
,
r
,
K
r
,
r
∈
1
,
…
,
R
,
π
i
=
π
i
,
1
,
π
i
,
2
,
…
,
π
i
,
R
,
i
≥
0
.It is well known that the vectorsπ
i, i
≥
0, satisfy the system (33)
π
0
,
π
1
,
…
Q
=
0
,
π
0
,
π
1
,
…
e
=
1
,where Q is the generator of ξ
t, t
≥
0. System (33) is infinite and cannot be directly solved on computer. However, it can be successfully solved by means of the numerically stable algorithms developed in [18, 19].


## 4. Performance Measures of the System

Having computed the vectors of the stationary probabilitiesπ
i, i
≥
0, it is possible to compute a variety of the steady state performance parameters of the system.The distribution of the number of the customers in orbit is(34)
lim
t
→
∞
⁡
P
i
t
=
i
=
π
i
e
,
i
≥
0
.The average number of customers in the service area is(35)
L
=
∑
i
=
0
∞
∑
r
=
1
R
∑
k
=
0
K
r
k
π
i
,
r
,
k
e
.Here and throughout this section, formulas for the main performance measures contain the infinite sums. However, this does not create essential difficulties in computer implementation. It is well known that if the ergodicity condition is fulfilled, the stationary probability vectorsπ
i converge in norm to zero vector when i approaches infinity. So, computation of an infinite sum may be stopped if the summand becomes less than some preassigned value ε (e.g., ε
=
1
0
-
10).The average number of customers in the buffer is(36)
N
buffer
=
∑
i
=
0
∞
∑
r
=
1
R
∑
k
=
N
r
+
1
K
r
k
-
N
r
π
i
,
r
,
k
e
.The average number of busy servers is(37)
N
server
=
∑
i
=
0
∞
∑
r
=
1
R
∑
k
=
1
K
r
min
⁡
k
,
N
r
π
i
,
r
,
k
e
.The average number of customers in orbit is(38)
L
orbit
=
∑
i
=
1
∞
i
π
i
e
.The intensity of output of customers is(39)
λ
out
=
∑
i
=
0
∞
∑
r
=
1
R
∑
k
=
1
K
r
min
⁡
k
,
N
r
μ
r
π
i
,
r
,
k
e
.The probability that a customer arrives at the system when the buffer is full and leaves the system is(40)
P
loss-ent
=
λ
-
1
∑
i
=
0
∞
∑
r
=
1
R
1
-
p
1
r
π
i
,
r
,
K
r
D
1
r
e
.The probability that a customer arrives at the system when the buffer is full and goes into orbit is(41)
P
orb-ent
=
λ
-
1
∑
i
=
0
∞
∑
r
=
1
R
p
1
r
π
i
,
r
,
K
r
D
1
r
e
.The loss probability of a customer is(42)
P
loss
=
1
-
λ
out
λ
.The probability of customers loss due to the decrease of the number of servers caused by change of the state of theR
E is (43)
P
RE-loss
=
1
λ
∑
i
=
0
∞
∑
r
=
2
R
∑
r
′
=
1
r
-
1
1
-
p
2
r
H
r
,
r
′
·
∑
k
=
N
r
′
+
1
N
r
k
-
N
r
′
π
i
,
r
,
k
+
N
r
-
N
r
′
∑
k
=
N
r
+
1
min
⁡
K
r
′
+
N
r
-
N
r
′
,
K
r
π
i
,
r
,
k
+
∑
k
=
K
r
′
+
N
r
-
N
r
′
+
1
K
r
k
-
K
r
′
π
i
,
r
,
k
·
e
.This formula is quite transparent. The right hand side of this formula represents the fraction. The denominator of this fraction is the average rate λ of customers arrival to the system. The numerator is the average rate of customers loss due to the decrease of the number of servers caused by the change of the state of the R
E
. Such a loss occurs, with probability 1
-
p
2
(
r
), every time when the state r,  r
∈
{
2
,
…
,
R
}, of the R
E changes to the state r
′,  r
′
∈
1
,
…
,
r
-
1, while the number k of customers in service area is greater than N
r
′
. The number of the simultaneously lost customers is equal to k
-
N
r
′ if k
∈
N
r
′
+
1
,
…
,
N
r, equal to N
r
-
N
r
′ if k
∈
N
r
+
1
,
…
,
min
⁡
K
r
′
+
N
r
-
N
r
′
,
K
r, and equal to k
-
K
r
′ if k
∈
K
r
′
+
N
r
-
N
r
′
+
1
,
…
,
K
r
. Applying the formula of total probability, we get formula (43).The probability of an arbitrary customer loss from orbit is(44)
P
(loss-from-orbit)
=
P
loss
-
P
loss-ent
-
P
RE-loss
.The probability that an arbitrary customer from orbit will make an attempt to receive service when the system is full and return to orbit is(45)
P
return-to-orbit
=
α
-
1
∑
i
=
1
∞
∑
r
=
1
R
i
α
r
p
3
r
π
i
,
r
,
K
r
e
,where α
=
∑
i
=
1
∞
∑
r
=
1
R
i
α
(
r
)
π
(
i
,
r
)
e
.The probability that an arbitrary customer from orbit makes an attempt to receive service when the system area is full and leaves the system without service is(46)
P
loss-non-persistence
=
α
-
1
∑
i
=
1
∞
∑
r
=
1
R
i
α
r
1
-
p
3
r
π
i
,
r
,
K
r
e
.


## 5. Numerical Results

Let us consider the following set of the system parameters. The number of the states of theR
E is R
=
2
. The generator of the R
E is given by (47)
H
=
-
0.03
0.03
0.07
-
0.07
,so the stationary probability of state 1 is φ
1
=
0.7 and the stationary probability of state 2 is φ
2
=
0.3
. These probabilities are the components of the vector φ
=
(
φ
1
,
φ
2
) which is computed as the unique solution to the system φ
H
=
0, φ
e
=
1
.We assume that the arrival flow under state 1 of theR
E is defined by the matrices (48)
D
0
1
=
-
2.7032
0
0
-
0.0877
,
D
1
1
=
2.6853
0.0179
0.0488
0.0389
,and under state 2 of the R
E it is defined by the matrices (49)
D
0
2
=
-
6.7582
0
0
-
0.2193
,
D
1
2
=
6.7133
0.0449
0.1221
0.0972
.For both arrival processes, the coefficient of correlation is c
c
o
r
(
r
)
=
0.2 and the coefficient of variation is c
v
a
r
(
r
)
=
12.34, r
=
1,2
. The average arrival rate λ
(
1
) of customers under state 1 of the R
E is 2, and the average arrival rate of customers λ
(
2
) under state 2 of the R
E is 5
.The rest of the system parameters under state 1 of theR
E are as follows: (50)
K
1
=
25
,
p
1
1
=
0.95
,
p
2
1
=
1
,
p
3
1
=
0.9
,
γ
1
=
0.1
,
α
1
=
0.3
,
μ
1
=
0.4
.Under state 2 of the R
E, the parameters are as follows: (51)
K
2
=
25
,
p
1
1
=
0.9
,
p
2
1
=
1
,
p
3
1
=
0.9
,
γ
2
=
0
,
α
2
=
0.5
,
μ
2
=
0.5
.Because arrival rateλ
(
2
) is 2.5 times higher than λ
(
1
) while the service rate μ
(
2
) is only 1.25 times higher than μ
(
1
), we may interpret state 1 as a normal mode of the system operation while we interpret state 2 as a congestion mode. To provide good quality of customers service, we should properly adjust the number of available servers to the mode of system. From economical considerations, let us assume that the average number of servers operating at an arbitrary moment of time should not exceed the predefined number N
~
. The average number of servers operating at an arbitrary moment evidently is computed by formula N
(
1
)
φ
1
+
N
(
2
)
φ
2
. Suppose that we have an opportunity to arbitrarily assign the numbers of available servers N
(
1
) and N
(
2
) that satisfy inequality (52)
N
1
φ
1
+
N
2
φ
2
≤
N
~aiming to provide good quality of customers service. It is intuitively obvious that the number N
(
2
) of available servers under operation in the congestion mode should be not less than the number of available servers N
(
1
) in the normal mode.Let us vary the number of available servers under the different states of theR
E, N
(
1
) and N
(
2
), in such a way that N
(
1
)
≤
N
(
2
) and inequality (52) holds good.Figures2–6 illustrate the dependence of some performance measures of the system for values N
(
1
) and N
(
2
) varying in the described set when N
~
=
9
.Figure 2
Dependence of the average number of customers in orbit on the numbers of serversN
(
1
) and N
(
2
).Figure 3
Dependence of the average number of customers in the buffer on the numbers of serversN
(
1
) and N
(
2
).Figure 4
Dependence of the average number of busy servers on the numbers of serversN
(
1
) and N
(
2
).Figure 5
Dependence of the probability that a customer arrives at the system when the buffer is full and leaves the system on the numbers of serversN
(
1
) and N
(
2
).Figure 6
Dependence of the loss probability of a customer on the numbers of serversN
(
1
) and N
(
2
).It is seen, from Figure2, that the average number of customers in orbit significantly changes (from about 0 to about 30) when N
(
1
) and N
(
2
) vary, and this number is maximal when both N
(
1
) and N
(
2
) are small and is minimal when about 9 servers are active in both modes.Similar observations can be made about the average numberN
b
u
f
f
e
r of busy servers based on Figure 3. However, the value of N
b
u
f
f
e
r is less sensitive with respect to the change of N
(
2
) compared to L
o
r
b
i
t
.It is seen, from Figure4, that the average number of busy servers also significantly changes (from about 0 to about 7) when N
(
1
) and N
(
2
) vary, and this number is minimal when both N
(
1
) and N
(
2
) are small and maximal when about 9 servers are active in both modes.Because the capacity of service area (the number of servers plus the number of places in a buffer) under both the states of theR
E is pretty large, 25, the probability that a customer arrives at the system when the buffer is full and leaves the system is not very high (about 0.08) even when the number of active servers is small; see Figure 5. This may seem surprising. But, this effect is easily explained by Figure 6. Due to the impatience of customers, the loss probability of customers is very high (close to 1) when the number of active servers is small. So, the system is never overcrowded. It is worth noting, based on Figure 6, that the loss probability of customers due to impatience essentially decreases with growth of N
(
2
) and, especially, growth of N
(
1
)
.All the mentioned performance measures of the system are quite important. However, likely, the most important measure from the point of view of the system manager is the loss probabilityP
(
l
o
s
s
) of an arbitrary customer because the loss of a customer implies the loss of a potential profit from his/her service. So, let us consider, namely, the loss probability P
(
l
o
s
s
) as criterion of the quality of the system operation and find the optimal number of the servers under the different states of the R
E, N
(
1
), and N
(
2
), under the condition that inequality (52) should be fulfilled.Because, to build Figure6, we have computed the probability P
(
l
o
s
s
) for all feasible N
(
1
) and N
(
2
), the solution of optimization problem is trivial. The optimal (minimal) value of the loss probability is P
(
l
o
s
s
)
=
0.02519 when N
(
1
)
=
7 and N
(
2
)
=
13
. The values of stationary probabilities π
(
i
,
r
) under the optimal values N
(
1
)
=
7 and N
(
2
)
=
13 are presented in Table 1. Probabilities π
(
i
,
r
) for i
>
15 are less than 0.001 and they are omitted here.Table 1
Stationary probabilities of the system states underN
(
1
)
=
7 and N
(
2
)
=
13.
i∖r
r
=
1
r
=
2
i
=
0
0.56478
0.23043
i
=
1
0.04455
0.01657
i
=
2
0.02590
0.01007
i
=
3
0.01818
0.00749
i
=
4
0.01352
0.00597
i
=
5
0.01027
0.00491
i
=
6
0.00786
0.00409
i
=
7
0.00414
0.00344
i
=
8
0.00296
0.00290
i
=
9
0.00214
0.00245
i
=
10
0.00154
0.00206
i
=
11
0.00111
0.00173
i
=
12
8.06E − 4
0.00145
i
=
13
5.86E − 4
0.00121
i
=
14
4.28E − 4
0.00100
i
=
15
3.15E − 4
8.25E − 4It is worth noting that ifN
(
1
)
=
N
(
2
)
=
9, that is, the number of active servers is constant and is not adjusted to the alternation of the normal and congestion mode, the value of the loss probability is about three times higher: P
(
l
o
s
s
)
=
0.067389
. So, the adjustment of the number of active servers to the current load of the system definitely makes sense.To briefly illustrate the profound effect of correlation and variance of interarrival times in the arrival process, let us assume now that, instead of theM
A
P flows with the coefficient of correlation 0.2 considered above, arrival flows are described by the stationary Poisson processes with the same mean arrival rates. These processes are defined by (53)
D
0
1
=
-
2
,
D
1
1
=
2
,
D
0
2
=
-
5
,
D
1
2
=
5
.Let us repeat the experiment described above with these arrival processes having zero correlation and coefficient of variation equal to 1. We obtain that the minimal value of the loss probability (P
(
l
o
s
s
)
=
0.005102) is achieved when N
(
1
)
=
8 and N
(
2
)
=
11
. At the same time, under these values of N
(
1
) and N
(
2
), P
(
l
o
s
s
) is equal to =0.037619 for arrival flows with c
c
o
r
(
r
)
=
0.2, r
=
1,2
. So, the ignorance of the possible correlation in the arrival process gives relative error in the prediction of the value of P
(
l
o
s
s
) equal to 0.037619
-
0.005102
/
0.005102
×
100
=
637 percent. This is not admissible in the performance evaluation and capacity planning of real life systems. Note that our experience of investigation of arrival flows in contact centers of several banks shows that these flows are indeed correlated.


## 6. Conclusion

We considered a multiserver queuing system, the dynamics of which depends on the state of theR
E. Change of the state of the R
E may cause variation of the parameters of arrival, service, retrial, impatience processes, number of available servers, and available buffer capacity, as well as the behavior of the customers. Evolution of the system is described by the multidimensional continuous-time Markov chain. The generator of this Markov chain is derived in a block matrix form. The ergodicity condition is presented. In particular, it is shown that if the customers are impatient or nonpersistent at least in one state of the R
E, then the Markov chain is ergodic for any set of the system parameters. Expressions for the key performance measures are given. Numerical results illustrating the behavior of the system and showing the possibility of formulation and solution of optimization problems are provided. Positive effect of adjusting the number of active servers to the current load of the system is demonstrated. The importance of the account of correlation in the arrival processes is numerically illustrated.



---

*Source: 180481-2015-07-21.xml*