# Multiserver Queue with Guard Channel for Priority and Retrial Customers

Authors: Kazuki Kajiwara; Tuan Phung-Duc
Journal: International Journal of Stochastic Analysis
 (2016)

---

## Abstract

This paper considers a retrial queueing model where a group of guard channels is reserved for priority and retrial customers. Priority and normal customers arrive at the system according to two distinct Poisson processes. Priority customers are accepted if there is an idle channel upon arrival while normal customers are accepted if and only if the number of idle channels is larger than the number of guard channels. Blocked customers (priority or normal) join a virtual orbit and repeat their attempts in a later time. Customers from the orbit (retrial customers) are accepted if there is an idle channel available upon arrival. We formulate the queueing system using a level dependent quasi-birth-and-death (QBD) process. We obtain a Taylor series expansion for the nonzero elements of the rate matrices of the level dependent QBD process. Using the expansion results, we obtain an asymptotic upper bound for the joint stationary distribution of the number of busy channels and that of customers in the orbit. Furthermore, we develop an efficient numerical algorithm to calculate the joint stationary distribution.


---

## Body

## 1. Introduction

In this paper, we consider multiserver retrial queues with guard channels for priority and retrial customers. Retrial queues are characterized by the fact that a blocked customer repeats its request after a random time. During retrial intervals, the customer is said to be in the orbit. This type of queueing models is widely used in modelling and performance analysis of communication and service systems, especially in cellular networks [1–6]. For instance, Tran-Gia and Mandjes [6] report the influence of retrials on the performance of cellular networks using retrial queueing models. Marsan et al. [3] carry out a fixed point approximation analysis for retrial queueing models arising from cellular networks while some extension is presented in [4]. Artalejo and Lopez-Herrero analyze a multiserver queue for cellular networks operating under a random environment using a four-dimensional Markov chain.The guard channel concept has been extensively used in communication systems [1–3, 6]. This is also referred to as trunk (or circuit or bandwidth) reservation in teletraffic literature [7, 8]. Tran-Gia and Mandjes [6] propose some multiserver retrial queues with fresh and handover calls and guard channels for a base station in cellular networks. In [6], the orbit size is assumed to be finite which simplifies the analysis.The analysis of multiserver retrial queues with infinite orbit size is challenging due to the fact that the underlying Markov chain is state dependent because the retrial rate is proportional to the number of customers in the orbit. Thus, even for the fundamental model with one type of traffic and without guard channels, the stationary distribution is expressed in terms of simple functions for only some special cases, that is, one or two servers [9]. We refer to [9–11] for some efforts in finding analytical expressions for the joint stationary distribution for the cases of more than two servers. For models with both retrial and guard channels, although some numerical methods [1, 3, 6, 12–14] have been presented for various models, there is no analytical result available.This motivates us to consider a new model with both retrials and guard channels for which we explore new analytical and numerical results. From the modelling point of view, the novelty is the priority given to retrial customers. It should be noted that retrial customers are treated the same as normal customers in [1–3, 5, 6] which is suitable for the context of cellular networks since the base station might not be able to recognize redial calls so as to give them some priority. In [15], retrial customers are given a preemptive priority over waiting customers. To the best of our knowledge, the current paper is the first to consider the priority for retrial customers in the context of queueing models with guard channels. Our model may also be fit for systems with human servers where the service differentiation among two classes of customers is needed. In such a service system, the server can easily recognize retrial customers so as to give them some priority over fresh customers who arrive at the system for the first time. We formulate the queueing system using a level dependent QBD process where the level and the phase are referred to as the number of customers in the orbit and the number of busy channels, respectively.The stationary distribution of a level dependent QBD process can be expressed in terms of a sequence of rate matrices [16]. Thus, we can characterize the stationary distribution through the sequence of rate matrices. The QBD process of our model possesses some special structure; that is, only the last two rows are nonzero allowing us to get some insights into the structure of the stationary distribution. Liu and Zhao [17] use this property to obtain upper and lower asymptotic bounds for the stationary distribution of the fundamental retrial model without guard channels. Liu et al. [18] further extend their analysis to the model with nonpersistent customers. B. Kim and J. Kim [19] and Kim et al. [20] refine the tail asymptotic results in Liu and Zhao [17] and Liu et al. [18], respectively. Phung-Duc [21] presents a perturbation analysis for a multiserver retrial queue with two types of nonpersistent customers. In [21], the author derives the Taylor series expansion formulae for the nonzero elements of the rate matrices. The difference of our model in comparison with the above work is that the last two rows of the rate matrices are nonzero in our model while for those in [17, 21] only the last row is nonzero. This makes the analysis of our model more complex and challenging.The main contribution of our paper is threefold. First, using a censoring technique and a perturbation method, we obtain the Taylor series expansion for the rate matrices in terms of the number of customers in the orbit. Our formula is general in the sense that we can obtain the expansion with arbitrary number of terms. This was not reported in Liu and Zhao [17]. Second, using this result we obtain an asymptotic upper bound for the stationary distribution which is more challenging compared to [17, 21] due to the denseness of the rate matrices. It should be noted that this is the first asymptotic result for multiserver retrial queue with guard channels. Third, we present an efficient method to calculate the stationary distribution of the model whose computational complexity is linear to the number of servers. An earlier version of this paper was presented in [22].The rest of our paper is organized as follows. Section2 presents the model and some preliminary results on the level dependent QBD formulation. Section 3 is devoted to the presentation of the Taylor series expansion for the rate matrices. In Section 4, we show the asymptotic upper bound for the joint stationary distribution while a numerical algorithm is presented in Section 5. Section 6 provides some numerical examples and Section 7 concludes our paper and presents some future directions.


## 2. Model and Formulation

### 2.1. Model

In this paper, we consider a queueing model with two types of customers (types 1 and 2). There arec servers; among them g servers are assigned as guard channels. Customers of type 1 (high priority) and type 2 (low priority) arrive at the system according to two independent Poisson processes with rates λ
1 and λ
2, respectively. Customers of type 1 can use all c servers while those of type 2 cannot use the guard channels. If there are c
-
g busy servers, the rest of g servers automatically become guard channels for customers of type 1. Furthermore, we assume that a blocked customer (both types 1 and 2) retries after some exponentially distributed time with mean 1
/
μ. Upon retrial, if there is an idle channel the retrial customer occupies it immediately; otherwise it enters the orbit again. Thus, a retrial customer has the same priority as that of a high priority one. As a result, we may expect that decreasing the number of retrials by a customer improves the quality of service (QoS). Service times for both types 1 and 2 customers are assumed to follow the same exponential distribution with mean 1
/
ν.Remark 1.
In this paper, we restrict ourselves to the case of one guard channel; that is,g
=
1. This is because the asymptotic analysis for case g
=
1 is complex enough and is essentially different from case g
=
0. Asymptotic analysis for case g
>
1 may need a new technique which will be left for a future study. On the other hand, the stability condition presented in this paper can be extended to case g
>
1 in a straightforward manner while the numerical algorithm can be adapted to case g
>
1.Remark 2.
From a theoretical point of view, the assumption that retrial customers (both normal and priority) have the same priority significantly simplifies the analysis. This is because if retrial customers keep their initial priority, we need to distinguish two types of retrial customers for which we should have two orbits for priority and normal customers.


### 2.2. Level Dependent QBD Process

LetC
(
t
) and N
(
t
) denote the number of busy servers and the number of retrial customers in the orbit at time t. Letting X
(
t
)
=
(
C
(
t
)
,
N
(
t
)
)
(
t
≥
0
), the bivariate process {
X
t
,
t
≥
0
} is a Markov chain in the state space S
=
{
0,1
,
…
,
c
}
×
Z
+, where Z
+
=
{
0,1
,
2
,
…
}. We assume that {
X
(
t
)
} is positive recurrent. The necessary and sufficient condition for the positive recurrence of {
X
(
t
)
} is given in Lemma 3.Lemma 3.
{
X
(
t
)
} is positive recurrent if and only if(1)
λ
c
ν
<
1
, where λ
=
λ
1
+
λ
2.Proof.
The proof is presented in AppendixA.It is easy to see that{
X
t
,
t
≥
0
} is a level dependent QBD process whose infinitesimal generator Q is given as follows:(2)
Q
=
Q
1
0
Q
0
0
O
O
⋯
Q
2
1
Q
1
1
Q
0
1
O
⋯
O
Q
2
2
Q
1
2
Q
0
2
⋯
O
O
Q
2
3
Q
1
3
⋯
⋮
⋮
⋮
⋮
⋱
,where O is the zero matrix with an appropriate dimension and {
Q
0
n
,
Q
1
n
,
n
∈
Z
+
} and {
Q
2
n
,
n
∈
N
} are square matrices of size c
+
1. Furthermore, bearing in mind that N
= North,
S
= South,
W
= West,and E
= East, these matrices are given and partitioned as follows:(3)
Q
0
n
=
0
⋯
0
0
0
⋮
⋱
⋮
⋮
⋮
0
⋯
0
0
0
0
⋯
0
λ
2
0
0
⋯
0
0
λ
=
Q
0
,
N
W
n
Q
0
,
N
E
n
Q
0
,
S
W
n
Q
0
,
S
E
n
,
Q
2
n
=
0
n
μ
0
⋯
0
0
0
0
0
n
μ
⋯
0
0
0
0
0
0
⋯
0
0
0
⋮
⋮
⋮
⋱
⋮
⋮
⋮
0
0
0
⋯
0
n
μ
0
0
0
0
⋯
0
0
n
μ
0
0
0
⋯
0
0
0
=
Q
2
,
N
W
n
Q
2
,
N
E
n
Q
2
,
S
W
n
Q
2
,
S
E
n
,
Q
1
n
=
b
0
n
λ
0
⋯
0
0
0
ν
b
1
n
λ
⋯
0
0
0
0
2
ν
b
2
n
⋯
0
0
0
⋮
⋮
⋮
⋱
⋮
⋮
⋮
⋮
⋮
⋮
⋱
b
c
-
2
n
λ
0
0
0
0
⋯
0
b
c
-
1
n
λ
1
0
0
0
⋯
0
c
ν
b
c
n
=
Q
1
,
N
W
n
Q
1
,
N
E
n
Q
1
,
S
W
n
Q
1
,
S
E
n
, where N
=
{
1,2
,
…
}, b
i
(
n
)
=
-
(
λ
+
i
ν
+
n
μ
(
1
-
δ
i
,
c
)
)
(
i
=
0,1
,
2
,
…
,
c
) and δ
i
,
c is the Kronecker symbol; that is, δ
i
,
c
=
1 if i
=
c and 0 otherwise. Let π
i
,
n denote the stationary probability that there are i busy servers and n customers in the orbit; that is,(4)
π
i
,
n
=
lim
t
→
∞
⁡
P
C
t
=
i
,
N
t
=
n
,
i
=
0,1
,
…
,
c
,
n
∈
Z
+
.Furthermore, let(5)
π
n
=
π
0
,
n
,
π
1
,
n
,
…
,
π
c
,
n
,
π
=
π
0
,
π
1
,
…
. We have(6)
π
0
Q
1
0
+
π
1
Q
2
1
=
0
,
n
=
0
,
π
n
-
1
Q
0
n
-
1
+
π
n
Q
1
n
+
π
n
+
1
Q
2
n
+
1
=
0
,
n
∈
N
,
π
e
=
1
,where e and 0 are vectors with appropriate dimensions with all 1 elements and all zero elements, respectively. It is established in [16] that the solution of (6) is given by(7)
π
n
=
π
n
-
1
R
n
,
n
∈
N
, where {
R
n
,
n
∈
N
} is the minimal nonnegative solution of(8)
Q
0
n
-
1
+
R
n
Q
1
n
+
R
n
R
n
+
1
Q
2
n
+
1
=
O
,
n
∈
N
.Furthermore, π
0 is determined by (9)
π
0
Q
1
0
+
R
1
Q
2
1
=
0
,
π
0
I
+
R
1
+
R
1
R
2
+
⋯
e
=
1
. Thus the problem of finding the stationary distribution is equivalent to that of obtaining the rate matrices. However, the rate matrices do not have closed forms in general leading to algorithmic approaches for numerical calculation. To this end, we present Lemmas 4 and 5.Lemma 4 (Proposition2.2 in [23]).
LetM denote the set of square matrices of size c
+
1. Furthermore, let R
n
:
M
→
M denote the following function:(10)
R
n
X
=
-
Q
0
n
-
1
Q
1
n
+
X
Q
2
n
+
1
-
1
,
n
∈
N
. It is easy to see that {
R
n
,
n
∈
N
} satisfies(11)
R
n
=
R
n
R
n
+
1
=
R
n
∘
R
n
+
1
∘
R
n
+
2
∘
⋯
,
n
∈
N
, where f
(
g
(
·
)
)
=
f
∘
g
(
·
).Lemma 5 (Proposition2.4 in [23]).
{
R
k
n
,
k
∈
Z
+
,
n
∈
N
} is defined by the following recursive formulae: (12)
R
0
n
=
O
,
k
=
0
,
R
k
n
=
R
n
R
k
-
1
n
+
1
=
⋯
=
R
n
∘
R
n
+
1
∘
⋯
∘
R
n
+
k
-
1
O
,
n
,
k
∈
N
.One has(13)
lim
k
→
∞
⁡
R
k
n
=
R
n
,
n
∈
N
.Remark 6.
Lemmas4 and 5 allow us to derive a numerical algorithm for calculating the rate matrices. They also imply that the rate matrices are matrix continued fractions. Bright and Taylor [24, 25] propose a recursive algorithm for computing rate matrix R
(
n
). From Lemmas 4 and 5, we observe that the first c
-
1 rows of R
(
n
) are zero. In Section 5, we propose a method for calculating R
(
n
) with computational complexity of O
(
c
) by exploiting this sparsity and Lemma 5. It should be noted that the recursive algorithm in [24, 25] has the computational complexity of O
(
c
3
) due to the denseness of D-matrices (see Bright and Taylor [24, 25]).It is easy to see that the firstc
-
1 rows of R
(
n
) and R
k
(
n
) are zeros. Similar to other block matrices, we also partition R
(
n
) as(14)
R
n
=
R
N
W
n
R
N
E
n
R
S
W
n
R
S
E
n
. It is easy to see that R
N
W
(
n
)
=
O and R
N
E
(
n
)
=
O. Furthermore, we denote the elements of R
S
W
(
n
) and R
S
E
(
n
) as follows:(15)
R
S
W
n
=
r
0
c
-
1
,
n
r
1
c
-
1
,
n
⋯
r
c
-
2
c
-
1
,
n
r
0
c
,
n
r
1
c
,
n
⋯
r
c
-
2
c
,
n
,
R
S
E
n
=
r
c
-
1
c
-
1
,
n
r
c
c
-
1
,
n
r
c
-
1
c
,
n
r
c
c
,
n
.Remark 7.
It should be noted that, in comparison with a previous version [22], some notations have been changed. In particular, r
i
(
0
,
n
) and r
i
(
1
,
n
) (i
=
0,1
,
2
,
…
,
c) in [22] are replaced by r
i
(
c
-
1
,
n
) and r
i
(
c
,
n
) in the current paper.In Section3, we show the Taylor series expansion of R
(
n
) in terms of 1
/
n. Comparing the last two rows on both sides of (8) yields(16)
R
S
W
n
Q
1
,
N
W
n
+
R
S
E
n
Q
1
,
N
W
n
+
R
S
E
n
R
S
W
n
+
1
Q
2
,
N
W
n
+
1
=
O
,
Q
0
,
S
E
n
-
1
+
R
S
W
n
Q
1
,
N
E
n
+
R
S
E
n
Q
1
,
N
E
n
+
R
S
E
n
R
S
W
n
Q
2
,
N
E
n
+
1
+
R
S
E
n
R
S
E
n
Q
2
,
S
E
n
+
1
=
O
.Rewriting (16) in the scalar form, we obtain(17)
b
0
n
r
0
c
-
1
,
n
+
ν
r
1
c
-
1
,
n
=
0
,
(18)
λ
r
i
-
1
c
-
1
,
n
+
b
i
n
r
i
c
-
1
,
n
+
i
+
1
ν
r
i
+
1
c
-
1
,
n
+
r
~
i
c
-
1
,
n
=
0
,
i
=
1,2
,
…
,
c
-
2
,
(19)
λ
r
c
-
2
c
-
1
,
n
+
b
c
-
1
n
r
c
-
1
c
-
1
,
n
+
c
ν
r
c
c
-
1
,
n
+
r
~
c
-
1
c
-
1
,
n
=
-
λ
2
,
(20)
λ
1
r
c
-
1
c
-
1
,
n
+
b
c
n
r
c
c
-
1
,
n
+
r
~
c
c
-
1
,
n
=
0
,
(21)
b
0
n
r
0
c
,
n
+
ν
r
1
c
,
n
=
0
,
(22)
λ
r
i
-
1
c
,
n
+
b
i
n
r
i
c
,
n
+
i
+
1
ν
r
i
+
1
c
,
n
+
r
~
i
c
,
n
=
0
,
i
=
1,2
,
…
,
c
-
1
,
(23)
λ
1
r
c
-
1
c
,
n
+
b
c
n
r
c
c
,
n
+
r
~
c
c
,
n
=
-
λ
,where(24)
r
~
i
c
-
1
,
n
=
n
+
1
μ
r
c
-
1
c
-
1
,
n
r
i
-
1
c
-
1
,
n
+
1
+
r
c
c
-
1
,
n
r
i
-
1
c
,
n
+
1
,
r
~
i
c
,
n
=
n
+
1
μ
r
c
-
1
c
,
n
r
i
-
1
c
-
1
,
n
+
1
+
r
c
c
,
n
r
i
-
1
c
,
n
+
1
.Lemma 8 (Lemma2.1 in [17, 18] and Proposition  3 in [21]).
One has(25)
Q
2
n
-
1
+
Q
1
n
-
1
+
R
n
Q
2
n
e
=
0
,
n
∈
N
.Comparing the last two elements on both sides of (25) yields(26)
∑
i
=
0
c
-
1
r
i
c
-
1
,
n
=
λ
2
n
μ
,
n
∈
N
,
(27)
∑
i
=
0
c
-
1
r
i
c
,
n
=
λ
n
μ
,
n
∈
N
.Proof.
This lemma follows from the fact that the following matrix represents the infinitesimal generator of the ergodic Markov chain{
X
t
,
t
≥
0
} censored in levels {
l
i
,
i
=
0,1
,
…
,
n
-
1
}, where l
(
i
)
=
(
(
0
,
i
)
,
(
1
,
i
)
,
…
,
(
c
,
i
)
). Let Q
≤
n
-
1 denote the infinitesimal generator of the censored Markov chain. We have(28)
Q
≤
n
-
1
=
Q
1
0
Q
0
0
O
⋯
O
Q
2
1
Q
1
1
Q
0
1
⋱
O
O
Q
2
2
Q
1
2
⋱
⋮
⋮
O
⋱
⋱
Q
0
n
-
2
O
⋯
O
Q
2
n
-
1
Q
^
n
-
1
,where(29)
Q
^
n
-
1
=
Q
1
n
-
1
+
R
n
Q
2
n
. Therefore,(30)
Q
2
n
-
1
+
Q
^
n
-
1
e
=
0
. By comparing the last two elements of both sides of this equation, we obtain the announced result.Corollary 9.
For casec
=
2, explicit expressions for nonzero elements of R
(
n
) are given as follows:(31)
r
0
1
,
n
=
λ
2
ν
n
μ
λ
+
ν
+
n
μ
,
r
1
1
,
n
=
λ
2
λ
+
n
μ
n
μ
λ
+
ν
+
n
μ
,due to (17) and (26) with c
=
2. It follows from (27) and (21) with c
=
2 that(32)
r
0
2
,
n
=
λ
ν
n
μ
λ
+
ν
+
n
μ
,
r
1
2
,
n
=
λ
λ
+
n
μ
n
μ
λ
+
ν
+
n
μ
.
Furthermore, substituting these explicit expressions into (20) and (23) and arranging the results, we obtain(33)
r
2
1
,
n
=
λ
2
λ
+
n
μ
λ
1
λ
+
ν
+
n
+
1
μ
+
λ
2
ν
n
μ
λ
+
ν
+
n
μ
3
λ
+
2
ν
+
2
n
+
1
μ
ν
,
r
2
2
,
n
=
λ
ν
λ
+
ν
+
n
+
1
μ
3
λ
+
2
ν
+
2
n
+
1
μ
+
λ
+
n
μ
λ
λ
+
n
+
1
μ
+
λ
1
ν
n
μ
λ
+
ν
+
n
μ
3
λ
+
2
ν
+
2
n
+
1
μ
.Remark 10.
It should be noted that the explicit results in Corollary9 cannot be obtained for case c
≥
3. Thus, in the next section, we present an asymptotic expansion for the rate matrices in the general case with an arbitrary value of c.


## 2.1. Model

In this paper, we consider a queueing model with two types of customers (types 1 and 2). There arec servers; among them g servers are assigned as guard channels. Customers of type 1 (high priority) and type 2 (low priority) arrive at the system according to two independent Poisson processes with rates λ
1 and λ
2, respectively. Customers of type 1 can use all c servers while those of type 2 cannot use the guard channels. If there are c
-
g busy servers, the rest of g servers automatically become guard channels for customers of type 1. Furthermore, we assume that a blocked customer (both types 1 and 2) retries after some exponentially distributed time with mean 1
/
μ. Upon retrial, if there is an idle channel the retrial customer occupies it immediately; otherwise it enters the orbit again. Thus, a retrial customer has the same priority as that of a high priority one. As a result, we may expect that decreasing the number of retrials by a customer improves the quality of service (QoS). Service times for both types 1 and 2 customers are assumed to follow the same exponential distribution with mean 1
/
ν.Remark 1.
In this paper, we restrict ourselves to the case of one guard channel; that is,g
=
1. This is because the asymptotic analysis for case g
=
1 is complex enough and is essentially different from case g
=
0. Asymptotic analysis for case g
>
1 may need a new technique which will be left for a future study. On the other hand, the stability condition presented in this paper can be extended to case g
>
1 in a straightforward manner while the numerical algorithm can be adapted to case g
>
1.Remark 2.
From a theoretical point of view, the assumption that retrial customers (both normal and priority) have the same priority significantly simplifies the analysis. This is because if retrial customers keep their initial priority, we need to distinguish two types of retrial customers for which we should have two orbits for priority and normal customers.


## 2.2. Level Dependent QBD Process

LetC
(
t
) and N
(
t
) denote the number of busy servers and the number of retrial customers in the orbit at time t. Letting X
(
t
)
=
(
C
(
t
)
,
N
(
t
)
)
(
t
≥
0
), the bivariate process {
X
t
,
t
≥
0
} is a Markov chain in the state space S
=
{
0,1
,
…
,
c
}
×
Z
+, where Z
+
=
{
0,1
,
2
,
…
}. We assume that {
X
(
t
)
} is positive recurrent. The necessary and sufficient condition for the positive recurrence of {
X
(
t
)
} is given in Lemma 3.Lemma 3.
{
X
(
t
)
} is positive recurrent if and only if(1)
λ
c
ν
<
1
, where λ
=
λ
1
+
λ
2.Proof.
The proof is presented in AppendixA.It is easy to see that{
X
t
,
t
≥
0
} is a level dependent QBD process whose infinitesimal generator Q is given as follows:(2)
Q
=
Q
1
0
Q
0
0
O
O
⋯
Q
2
1
Q
1
1
Q
0
1
O
⋯
O
Q
2
2
Q
1
2
Q
0
2
⋯
O
O
Q
2
3
Q
1
3
⋯
⋮
⋮
⋮
⋮
⋱
,where O is the zero matrix with an appropriate dimension and {
Q
0
n
,
Q
1
n
,
n
∈
Z
+
} and {
Q
2
n
,
n
∈
N
} are square matrices of size c
+
1. Furthermore, bearing in mind that N
= North,
S
= South,
W
= West,and E
= East, these matrices are given and partitioned as follows:(3)
Q
0
n
=
0
⋯
0
0
0
⋮
⋱
⋮
⋮
⋮
0
⋯
0
0
0
0
⋯
0
λ
2
0
0
⋯
0
0
λ
=
Q
0
,
N
W
n
Q
0
,
N
E
n
Q
0
,
S
W
n
Q
0
,
S
E
n
,
Q
2
n
=
0
n
μ
0
⋯
0
0
0
0
0
n
μ
⋯
0
0
0
0
0
0
⋯
0
0
0
⋮
⋮
⋮
⋱
⋮
⋮
⋮
0
0
0
⋯
0
n
μ
0
0
0
0
⋯
0
0
n
μ
0
0
0
⋯
0
0
0
=
Q
2
,
N
W
n
Q
2
,
N
E
n
Q
2
,
S
W
n
Q
2
,
S
E
n
,
Q
1
n
=
b
0
n
λ
0
⋯
0
0
0
ν
b
1
n
λ
⋯
0
0
0
0
2
ν
b
2
n
⋯
0
0
0
⋮
⋮
⋮
⋱
⋮
⋮
⋮
⋮
⋮
⋮
⋱
b
c
-
2
n
λ
0
0
0
0
⋯
0
b
c
-
1
n
λ
1
0
0
0
⋯
0
c
ν
b
c
n
=
Q
1
,
N
W
n
Q
1
,
N
E
n
Q
1
,
S
W
n
Q
1
,
S
E
n
, where N
=
{
1,2
,
…
}, b
i
(
n
)
=
-
(
λ
+
i
ν
+
n
μ
(
1
-
δ
i
,
c
)
)
(
i
=
0,1
,
2
,
…
,
c
) and δ
i
,
c is the Kronecker symbol; that is, δ
i
,
c
=
1 if i
=
c and 0 otherwise. Let π
i
,
n denote the stationary probability that there are i busy servers and n customers in the orbit; that is,(4)
π
i
,
n
=
lim
t
→
∞
⁡
P
C
t
=
i
,
N
t
=
n
,
i
=
0,1
,
…
,
c
,
n
∈
Z
+
.Furthermore, let(5)
π
n
=
π
0
,
n
,
π
1
,
n
,
…
,
π
c
,
n
,
π
=
π
0
,
π
1
,
…
. We have(6)
π
0
Q
1
0
+
π
1
Q
2
1
=
0
,
n
=
0
,
π
n
-
1
Q
0
n
-
1
+
π
n
Q
1
n
+
π
n
+
1
Q
2
n
+
1
=
0
,
n
∈
N
,
π
e
=
1
,where e and 0 are vectors with appropriate dimensions with all 1 elements and all zero elements, respectively. It is established in [16] that the solution of (6) is given by(7)
π
n
=
π
n
-
1
R
n
,
n
∈
N
, where {
R
n
,
n
∈
N
} is the minimal nonnegative solution of(8)
Q
0
n
-
1
+
R
n
Q
1
n
+
R
n
R
n
+
1
Q
2
n
+
1
=
O
,
n
∈
N
.Furthermore, π
0 is determined by (9)
π
0
Q
1
0
+
R
1
Q
2
1
=
0
,
π
0
I
+
R
1
+
R
1
R
2
+
⋯
e
=
1
. Thus the problem of finding the stationary distribution is equivalent to that of obtaining the rate matrices. However, the rate matrices do not have closed forms in general leading to algorithmic approaches for numerical calculation. To this end, we present Lemmas 4 and 5.Lemma 4 (Proposition2.2 in [23]).
LetM denote the set of square matrices of size c
+
1. Furthermore, let R
n
:
M
→
M denote the following function:(10)
R
n
X
=
-
Q
0
n
-
1
Q
1
n
+
X
Q
2
n
+
1
-
1
,
n
∈
N
. It is easy to see that {
R
n
,
n
∈
N
} satisfies(11)
R
n
=
R
n
R
n
+
1
=
R
n
∘
R
n
+
1
∘
R
n
+
2
∘
⋯
,
n
∈
N
, where f
(
g
(
·
)
)
=
f
∘
g
(
·
).Lemma 5 (Proposition2.4 in [23]).
{
R
k
n
,
k
∈
Z
+
,
n
∈
N
} is defined by the following recursive formulae: (12)
R
0
n
=
O
,
k
=
0
,
R
k
n
=
R
n
R
k
-
1
n
+
1
=
⋯
=
R
n
∘
R
n
+
1
∘
⋯
∘
R
n
+
k
-
1
O
,
n
,
k
∈
N
.One has(13)
lim
k
→
∞
⁡
R
k
n
=
R
n
,
n
∈
N
.Remark 6.
Lemmas4 and 5 allow us to derive a numerical algorithm for calculating the rate matrices. They also imply that the rate matrices are matrix continued fractions. Bright and Taylor [24, 25] propose a recursive algorithm for computing rate matrix R
(
n
). From Lemmas 4 and 5, we observe that the first c
-
1 rows of R
(
n
) are zero. In Section 5, we propose a method for calculating R
(
n
) with computational complexity of O
(
c
) by exploiting this sparsity and Lemma 5. It should be noted that the recursive algorithm in [24, 25] has the computational complexity of O
(
c
3
) due to the denseness of D-matrices (see Bright and Taylor [24, 25]).It is easy to see that the firstc
-
1 rows of R
(
n
) and R
k
(
n
) are zeros. Similar to other block matrices, we also partition R
(
n
) as(14)
R
n
=
R
N
W
n
R
N
E
n
R
S
W
n
R
S
E
n
. It is easy to see that R
N
W
(
n
)
=
O and R
N
E
(
n
)
=
O. Furthermore, we denote the elements of R
S
W
(
n
) and R
S
E
(
n
) as follows:(15)
R
S
W
n
=
r
0
c
-
1
,
n
r
1
c
-
1
,
n
⋯
r
c
-
2
c
-
1
,
n
r
0
c
,
n
r
1
c
,
n
⋯
r
c
-
2
c
,
n
,
R
S
E
n
=
r
c
-
1
c
-
1
,
n
r
c
c
-
1
,
n
r
c
-
1
c
,
n
r
c
c
,
n
.Remark 7.
It should be noted that, in comparison with a previous version [22], some notations have been changed. In particular, r
i
(
0
,
n
) and r
i
(
1
,
n
) (i
=
0,1
,
2
,
…
,
c) in [22] are replaced by r
i
(
c
-
1
,
n
) and r
i
(
c
,
n
) in the current paper.In Section3, we show the Taylor series expansion of R
(
n
) in terms of 1
/
n. Comparing the last two rows on both sides of (8) yields(16)
R
S
W
n
Q
1
,
N
W
n
+
R
S
E
n
Q
1
,
N
W
n
+
R
S
E
n
R
S
W
n
+
1
Q
2
,
N
W
n
+
1
=
O
,
Q
0
,
S
E
n
-
1
+
R
S
W
n
Q
1
,
N
E
n
+
R
S
E
n
Q
1
,
N
E
n
+
R
S
E
n
R
S
W
n
Q
2
,
N
E
n
+
1
+
R
S
E
n
R
S
E
n
Q
2
,
S
E
n
+
1
=
O
.Rewriting (16) in the scalar form, we obtain(17)
b
0
n
r
0
c
-
1
,
n
+
ν
r
1
c
-
1
,
n
=
0
,
(18)
λ
r
i
-
1
c
-
1
,
n
+
b
i
n
r
i
c
-
1
,
n
+
i
+
1
ν
r
i
+
1
c
-
1
,
n
+
r
~
i
c
-
1
,
n
=
0
,
i
=
1,2
,
…
,
c
-
2
,
(19)
λ
r
c
-
2
c
-
1
,
n
+
b
c
-
1
n
r
c
-
1
c
-
1
,
n
+
c
ν
r
c
c
-
1
,
n
+
r
~
c
-
1
c
-
1
,
n
=
-
λ
2
,
(20)
λ
1
r
c
-
1
c
-
1
,
n
+
b
c
n
r
c
c
-
1
,
n
+
r
~
c
c
-
1
,
n
=
0
,
(21)
b
0
n
r
0
c
,
n
+
ν
r
1
c
,
n
=
0
,
(22)
λ
r
i
-
1
c
,
n
+
b
i
n
r
i
c
,
n
+
i
+
1
ν
r
i
+
1
c
,
n
+
r
~
i
c
,
n
=
0
,
i
=
1,2
,
…
,
c
-
1
,
(23)
λ
1
r
c
-
1
c
,
n
+
b
c
n
r
c
c
,
n
+
r
~
c
c
,
n
=
-
λ
,where(24)
r
~
i
c
-
1
,
n
=
n
+
1
μ
r
c
-
1
c
-
1
,
n
r
i
-
1
c
-
1
,
n
+
1
+
r
c
c
-
1
,
n
r
i
-
1
c
,
n
+
1
,
r
~
i
c
,
n
=
n
+
1
μ
r
c
-
1
c
,
n
r
i
-
1
c
-
1
,
n
+
1
+
r
c
c
,
n
r
i
-
1
c
,
n
+
1
.Lemma 8 (Lemma2.1 in [17, 18] and Proposition  3 in [21]).
One has(25)
Q
2
n
-
1
+
Q
1
n
-
1
+
R
n
Q
2
n
e
=
0
,
n
∈
N
.Comparing the last two elements on both sides of (25) yields(26)
∑
i
=
0
c
-
1
r
i
c
-
1
,
n
=
λ
2
n
μ
,
n
∈
N
,
(27)
∑
i
=
0
c
-
1
r
i
c
,
n
=
λ
n
μ
,
n
∈
N
.Proof.
This lemma follows from the fact that the following matrix represents the infinitesimal generator of the ergodic Markov chain{
X
t
,
t
≥
0
} censored in levels {
l
i
,
i
=
0,1
,
…
,
n
-
1
}, where l
(
i
)
=
(
(
0
,
i
)
,
(
1
,
i
)
,
…
,
(
c
,
i
)
). Let Q
≤
n
-
1 denote the infinitesimal generator of the censored Markov chain. We have(28)
Q
≤
n
-
1
=
Q
1
0
Q
0
0
O
⋯
O
Q
2
1
Q
1
1
Q
0
1
⋱
O
O
Q
2
2
Q
1
2
⋱
⋮
⋮
O
⋱
⋱
Q
0
n
-
2
O
⋯
O
Q
2
n
-
1
Q
^
n
-
1
,where(29)
Q
^
n
-
1
=
Q
1
n
-
1
+
R
n
Q
2
n
. Therefore,(30)
Q
2
n
-
1
+
Q
^
n
-
1
e
=
0
. By comparing the last two elements of both sides of this equation, we obtain the announced result.Corollary 9.
For casec
=
2, explicit expressions for nonzero elements of R
(
n
) are given as follows:(31)
r
0
1
,
n
=
λ
2
ν
n
μ
λ
+
ν
+
n
μ
,
r
1
1
,
n
=
λ
2
λ
+
n
μ
n
μ
λ
+
ν
+
n
μ
,due to (17) and (26) with c
=
2. It follows from (27) and (21) with c
=
2 that(32)
r
0
2
,
n
=
λ
ν
n
μ
λ
+
ν
+
n
μ
,
r
1
2
,
n
=
λ
λ
+
n
μ
n
μ
λ
+
ν
+
n
μ
.
Furthermore, substituting these explicit expressions into (20) and (23) and arranging the results, we obtain(33)
r
2
1
,
n
=
λ
2
λ
+
n
μ
λ
1
λ
+
ν
+
n
+
1
μ
+
λ
2
ν
n
μ
λ
+
ν
+
n
μ
3
λ
+
2
ν
+
2
n
+
1
μ
ν
,
r
2
2
,
n
=
λ
ν
λ
+
ν
+
n
+
1
μ
3
λ
+
2
ν
+
2
n
+
1
μ
+
λ
+
n
μ
λ
λ
+
n
+
1
μ
+
λ
1
ν
n
μ
λ
+
ν
+
n
μ
3
λ
+
2
ν
+
2
n
+
1
μ
.Remark 10.
It should be noted that the explicit results in Corollary9 cannot be obtained for case c
≥
3. Thus, in the next section, we present an asymptotic expansion for the rate matrices in the general case with an arbitrary value of c.


## 3. Taylor Series Expansion

In this section, we derive the Taylor series expansion for all the nonzero elements of the rate matrices. In particular, we find the Taylor series expansion ofr
i
(
c
-
1
,
n
) and r
i
(
c
,
n
)
(
i
=
0,1
,
…
,
c
) in terms of 1
/
n. We use {
θ
m
0
,
k
,
m
∈
Z
+
} and {
θ
m
1
,
k
,
m
∈
Z
+
} to denote the coefficients of the Taylor series expansion, where k is the number of idle servers. We use the convention that if k
<
0 or c
<
k then θ
m
(
0
,
k
)
=
0 and θ
m
(
1
,
k
)
=
0. Furthermore, o
(
x
) implies l
i
m
x
→
0
o
(
x
)
/
x
=
0 and O
(
x
) implies l
i
m
s
u
p
x
→
0
O
(
x
)
/
x
<
∞, respectively.In this section, Lemma12 gives the one-term expansion while Lemma 13 improves Lemma 12 by replacing small order o
(
·
) by big order O
(
·
). Furthermore, Theorem 14 provides the general expansion formulae for the higher order Taylor series expansions of r
i
(
c
-
1
,
n
) and r
i
(
c
,
n
)
(
i
=
0,1
,
…
,
c
).Remark 11.
In this section, we find the Taylor expansion for nonzero components. The basic idea in a perturbation approach is that the coefficient of(
m
+
1
)th term is derived based on the coefficients of the lower term expansions, that is, ith term expansion (i
=
m
,
m
-
1
,
…
,
1). In this paper, we find the coefficients of (
m
+
1
)th term of r
c
-
k
(
c
-
1
,
n
) and r
c
-
k
(
c
,
n
) for k
=
0,1
,
…
,
c in parallel. In principle, Lemmas 12 and 13 and Theorem 14 below could be merged into one theorem; we however present these separately in order for clarity.Lemma 12.
One has the one-term series expansion (n
→
∞) for nonzero elements of R
(
n
) as follows: (34)
r
c
-
k
c
-
1
,
n
=
θ
0
0
,
k
1
n
k
+
o
1
n
k
,
k
=
0,1
,
…
,
c
,
r
c
-
k
c
,
n
=
θ
0
1
,
k
1
n
k
+
o
1
n
k
,
k
=
0,1
,
…
,
c
, where the sequences {
θ
0
0
,
k
,
k
=
0,1
,
…
,
c
} and {
θ
0
1
,
k
,
k
=
0,1
,
…
,
c
} are given as follows:(35)
θ
0
0
,
k
=
0
,
k
=
0
,
λ
2
μ
,
k
=
1
,
λ
2
μ
∏
i
=
1
k
-
1
c
-
i
ν
μ
,
k
=
2
,
…
,
c
,
θ
0
1
,
k
=
λ
c
ν
,
k
=
0
,
λ
μ
,
k
=
1
,
λ
μ
∏
i
=
1
k
-
1
c
-
i
ν
μ
,
k
=
2
,
…
,
c
.Proof.
The technical details are provided in AppendixB.Lemma 13.
The series expansion formulae (n
→
∞) in Lemma 12 can be improved as(36)
r
c
-
k
c
-
1
,
n
=
θ
0
0
,
k
1
n
k
+
O
1
n
k
+
1
,
k
=
0,1
,
…
,
c
,
(37)
r
c
-
k
c
,
n
=
θ
0
1
,
k
1
n
k
+
O
1
n
k
+
1
,
k
=
0,1
,
…
,
c
.Proof.
The technical details are provided in AppendixC.Theorem 14.
The nonzero elements ofR
(
n
) (n
→
∞) can be expanded as follows:(38)
r
c
-
k
c
-
1
,
n
=
∑
i
=
0
m
θ
i
0
,
k
-
1
i
1
n
k
+
i
+
O
1
n
k
+
m
+
1
,
m
∈
N
,
(39)
r
c
-
k
c
,
n
=
∑
i
=
0
m
θ
i
1
,
k
-
1
i
1
n
k
+
i
+
O
1
n
k
+
m
+
1
,
m
∈
N
,where {
θ
m
0
,
k
,
θ
m
1
,
k
,
k
=
0,1
,
…
,
c
,
m
∈
N
} are recursively defined as follows:(40)
θ
m
0,0
=
-
λ
1
c
ν
θ
m
-
1
0,1
+
μ
c
ν
∑
j
=
0
m
-
1
Φ
j
0,0
θ
m
-
j
-
1
0,1
-
1
j
+
1
+
μ
c
ν
∑
j
=
1
m
Φ
~
j
1,0
θ
m
-
j
0,0
-
1
j
,
θ
m
0,1
=
∑
j
=
2
m
i
n
c
,
m
+
1
θ
m
+
1
-
j
0
,
j
-
1
j
,
θ
m
0
,
k
=
c
-
k
+
1
ν
μ
θ
m
0
,
k
-
1
+
λ
μ
θ
m
-
2
0
,
k
+
1
+
λ
+
c
-
k
ν
μ
θ
m
-
1
0
,
k
+
∑
j
=
0
m
-
2
Φ
j
0
,
k
θ
m
-
j
-
2
0,1
-
1
j
+
∑
j
=
0
m
-
1
Φ
j
1
,
k
θ
m
-
j
-
1
0,0
-
1
j
+
1
,
k
=
2,3
,
…
,
c
,
θ
m
1,0
=
-
λ
1
c
ν
θ
m
-
1
1,1
+
μ
c
ν
∑
j
=
0
m
-
1
Φ
j
0,0
θ
m
-
j
-
1
1,1
-
1
j
+
1
+
μ
c
ν
∑
j
=
1
m
Φ
~
j
1,0
θ
m
-
j
1,0
-
1
j
,
θ
m
1,1
=
∑
j
=
2
m
i
n
c
,
m
+
1
θ
m
+
1
-
j
1
,
j
-
1
j
,
θ
m
1
,
k
=
c
-
k
+
1
ν
μ
θ
m
1
,
k
-
1
+
λ
μ
θ
m
-
2
1
,
k
+
1
+
λ
+
c
-
k
ν
μ
θ
m
-
1
1
,
k
+
∑
j
=
0
m
-
2
Φ
j
0
,
k
θ
m
-
j
-
2
1,1
-
1
j
+
∑
j
=
0
m
-
1
Φ
j
1
,
k
θ
m
-
j
-
1
1,0
-
1
j
+
1
,
k
=
2,3
,
…
,
c
. Furthermore,(41)
Φ
j
0
,
k
=
∑
i
=
0
j
θ
i
0
,
k
+
1
-
1
j
k
+
i
j
-
i
j
-
i
!
,
Φ
j
1
,
k
=
∑
i
=
0
j
θ
i
1
,
k
+
1
-
1
j
k
+
i
j
-
i
j
-
i
!
,
Φ
~
j
1,0
=
∑
i
=
1
j
θ
i
1,1
-
1
j
i
j
-
i
j
-
i
!
, where (
ϕ
)
n
(
-
∞
<
ϕ
<
∞
,
n
∈
Z
+
) denotes the Pochhammer symbol defined by (42)
ϕ
n
=
1
,
n
=
0
,
ϕ
ϕ
+
1
⋯
ϕ
+
n
-
1
,
n
∈
N
.Proof.
The technical details are provided in AppendixD.


## 4. Asymptotic Upper Bound

In this section, we present an asymptotic upper bound for the stationary distribution. To this end, we use Lemmas15 and 16.Lemma 15.
For a square matrixA
=
a
1,1
⋯
a
1
,
n
⋮
⋱
⋮
a
n
,
1
⋯
a
n
,
n and a vector x
=
(
x
1
,
x
2
,
…
,
x
n
), one has(43)
x
A
1
≤
x
1
A
∞
, where x
1
=
∑
i
=
1
n
x
i and A
∞
=
m
a
x
1
≤
i
≤
n
∑
j
=
1
n
a
i
j.Lemma 16 (Fact5 in [17]).
For integersn
(
≥
1) and a
^
>
0
,
b
^ satisfying b
^
≠
a
^
m
-
m
2
(
m
=
0,1
,
…
), one has (44)
∏
j
=
N
n
1
+
a
^
j
+
b
^
j
2
=
O
n
a
^
,
n
⟶
∞
.Remark 17.
In [17], only the last row of the rate matrices is nonzero. This fact allows us to evaluate the tail probability using the product of a sequence of scalars. However, since the last two rows of the rate matrices are nonzero in our model, we need to deal with the product of a sequence of matrices. Thus, in order to apply the technique given in [17], that is, Lemma 16, we need to use Lemma 15.Theorem 18.
One definesπ
n
E
=
(
π
c
-
1
,
n
,
π
c
,
n
). One then has (45)
π
n
E
1
=
O
n
a
×
λ
c
ν
n
,
n
⟶
∞
, where a
=
(
c
2
ν
+
λ
)
/
c
μ.Proof.
The proof uses Lemmas15 and 16. We have (46)
R
S
E
n
=
r
c
-
1
c
-
1
,
n
r
c
c
-
1
,
n
r
c
-
1
c
,
n
r
c
c
,
n
, where(47)
r
c
-
1
c
-
1
,
n
=
θ
0
0,1
1
n
-
θ
1
0,1
1
n
2
+
O
1
n
3
,
r
c
c
-
1
,
n
=
-
θ
1
0,0
1
n
+
θ
2
0,0
1
n
2
+
O
1
n
3
,
r
c
-
1
c
,
n
=
θ
0
1,1
1
n
-
θ
1
1,1
1
n
2
+
O
1
n
3
,
r
c
c
,
n
=
θ
0
1,0
-
θ
1
1,0
1
n
+
θ
2
1,0
1
n
2
+
O
1
n
3
.
It follows fromπ
n
=
π
n
-
1
R
(
n
) that π
n
E
=
π
n
-
1
E
R
S
E
(
n
). Thus, applying Lemma 15 repeatedly, we obtain (48)
π
n
E
1
≤
π
0
E
1
R
S
E
1
∞
⋯
R
S
E
n
-
1
∞
R
S
E
n
∞
.
For a sufficiently largen, R
S
E
n
∞ is given by (49)
R
S
E
n
∞
=
r
c
-
1
c
,
n
+
r
c
c
,
n
=
θ
0
1,0
+
θ
0
1,1
-
θ
1
1,0
1
n
+
θ
2
1,0
-
θ
1
1,1
1
n
2
+
O
1
n
3
=
θ
0
1,0
1
+
θ
0
1,1
-
θ
1
1,0
θ
0
1,0
n
+
θ
2
1,0
-
θ
1
1,1
θ
0
1,0
n
2
+
O
1
n
3
=
λ
c
ν
1
+
a
n
+
b
n
2
+
O
1
n
3
, where(50)
a
=
θ
0
1,1
-
θ
1
1,0
θ
0
1,0
,
b
=
θ
2
1,0
-
θ
1
1,1
θ
0
1,0
.
Thus, for the parameters that satisfy Lemma16, we have(51)
π
0
E
1
∏
i
=
1
n
R
S
E
i
∞
=
O
n
a
×
λ
c
ν
n
,
n
⟶
∞
, implying the desired result.Corollary 19.
One has(52)
π
i
,
n
=
O
n
a
-
c
+
i
×
λ
c
ν
n
,
i
=
0,1
,
…
,
c
,
n
⟶
∞
.Proof.
Fromπ
n
=
π
n
-
1
R
(
n
), we have(53)
π
i
,
n
=
π
c
-
1
,
n
-
1
r
i
c
-
1
,
n
+
π
c
,
n
-
1
r
i
c
,
n
,
i
=
0,1
,
…
,
c
. It follows from Theorem 14 that(54)
r
i
c
-
1
,
n
=
O
1
n
c
-
i
,
r
i
c
,
n
=
O
1
n
c
-
i
,
n
⟶
∞
. Theorem 18 yields(55)
π
n
E
1
=
O
n
a
×
λ
c
ν
n
,
n
⟶
∞
. Thus,(56)
π
i
,
n
=
O
n
a
-
c
+
i
×
λ
c
ν
n
,
n
⟶
∞
.


## 5. Numerical Algorithm

In this section, we propose a computational algorithm for the stationary distribution of our model extending that proposed by Phung-Duc et al. [26] for the fundamental M/M/c/c retrial queues without guard channels. In Section 5.1, we show some results which are the basis for the algorithm. Section 5.2 presents the algorithms for the rate matrices and the stationary distribution. Section 5.3 proposes a simple method for determining the truncation point used in Algorithm 2 in Section 5.2 for the stationary distribution. Section 5.4 derives some performance measures such as the blocking probability for type 2 customers (low priority) and that for type 1 and retrial calls.


### 5.1. Efficient Computation

Due to Lemma4, we need to compute k inverse matrices in order to obtain R
k
(
n
). It may take a long time when the number of servers is large. Thus, instead of computing the inverse matrices, we propose a new method exploiting the fact that only the last two rows are nonzero. The computational complexity of our new method is only O
(
c
). In particular, the computational complexity in all the theorems and lemmas below is O
(
c
).It should be noted that the computation ofR
(
n
) and R
k
(
n
) is equivalent to that of their last two rows r
(
n
) and r
k
(
n
); that is,(57)
r
n
=
r
c
-
1
,
n
r
c
,
n
,
r
k
n
=
r
k
c
-
1
,
n
r
k
c
,
n
,where r
(
i
,
n
) and r
k
(
i
,
n
) (i
=
0,1) are the row vectors of c
+
1 elements.Definition 20.
We define the functionr
n as follows. Let(58)
X
x
,
y
=
O
x
y
,
L
r
Y
=
y
c
-
1
y
c
, where y
c
-
1 and y
c are the second last and the last rows of Y. Furthermore, (59)
r
n
x
y
=
L
r
R
n
X
x
,
y
, where x and y are vectors with an appropriate dimension.It is easy to see thatr
(
n
) and r
k
(
n
) satisfy the following equations: (60)
r
n
=
r
n
r
n
+
1
,
r
k
n
=
r
n
r
k
-
1
n
+
1
=
r
n
∘
r
n
+
1
∘
⋯
∘
r
n
+
k
-
1
O
, for n
,
k
∈
N. Lemmas 21 and 22 compute r
k
(
c
-
1
,
n
) and r
k
(
c
,
n
) using r
k
-
1
(
n
+
1
), respectively. Furthermore, Lemma 23 computes the stationary distribution of the censored Markov chain on level 0 using r
(
1
).Lemma 21.
For arbitraryn, k, one has(61)
r
k
,
i
c
-
1
,
n
=
α
i
+
β
i
r
k
,
c
c
-
1
,
n
,
i
=
0,1
,
…
,
c
-
1
, where {
α
i
,
β
i
,
i
=
0,1
,
…
,
c
} and r
k
,
c
(
c
-
1
,
n
) are given as follows: (62)
α
c
=
0
,
β
c
=
1
,
α
c
-
1
=
0
,
β
c
-
1
=
-
b
c
n
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
,
n
+
1
λ
1
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
-
1
,
n
+
1
,
α
c
-
2
=
-
λ
2
λ
,
β
c
-
2
=
-
b
c
-
1
n
β
c
-
1
+
c
ν
+
n
+
1
μ
r
k
-
1
,
c
-
2
c
-
1
,
n
+
1
β
c
-
1
+
n
+
1
μ
r
k
-
1
,
c
-
2
c
,
n
+
1
λ
,
α
i
-
1
=
-
b
i
n
α
i
+
i
+
1
ν
α
i
+
1
λ
,
i
=
c
-
2
,
c
-
3
,
…
,
1
,
β
i
-
1
=
-
b
i
n
β
i
+
i
+
1
ν
β
i
+
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
β
c
-
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
,
n
+
1
λ
,
i
=
c
-
2
,
c
-
3
,
…
,
1
,
r
k
,
c
c
-
1
,
n
=
-
b
0
n
α
0
+
ν
α
1
b
0
n
β
0
+
ν
β
1
.Proof.
The technical details are provided in AppendixE.Lemma 22.
For arbitraryn and k, one has (63)
r
k
,
i
c
,
n
=
α
i
+
β
i
r
k
,
c
c
,
n
,
i
=
0,1
,
…
,
c
-
1
, where {
α
i
,
β
i
,
i
=
0,1
,
…
,
c
} and r
k
,
c
(
c
,
n
) are given as follows: (64)
α
c
=
0
,
β
c
=
1
,
α
c
-
1
=
-
λ
λ
1
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
-
1
,
n
+
1
,
β
c
-
1
=
-
b
c
n
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
,
n
+
1
λ
1
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
-
1
,
n
+
1
,
α
i
-
1
=
-
b
i
n
α
i
+
i
+
1
ν
α
i
+
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
α
c
-
1
λ
,
i
=
c
-
1
,
c
-
2
,
…
,
1
,
β
i
-
1
=
-
b
i
n
β
i
+
i
+
1
ν
β
i
+
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
β
c
-
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
,
n
+
1
λ
,
i
=
c
-
1
,
c
-
2
,
…
,
1
. Furthermore,(65)
r
k
,
c
c
,
n
=
-
b
0
n
α
0
+
ν
α
1
b
0
n
β
0
+
ν
β
1
.Proof.
This lemma can be proved using the same technique as in Lemma21.Lemma 23.
Solutionx
0
=
(
x
0
,
x
1
,
…
,
x
c
) to (66)
x
0
Q
1
0
+
R
1
Q
2
1
=
0
,
x
0
e
=
1
, is given by x
i
=
β
i
x
c
(
i
=
0,1
,
…
,
c
), where {
β
i
,
i
=
0,1
,
…
,
c
-
1
,
c
} is recursively defined as (67)
β
c
=
1
,
β
c
-
1
=
λ
+
c
ν
-
μ
r
c
-
1
1,1
λ
1
+
μ
r
c
-
1
0,1
,
β
i
-
1
=
λ
+
i
ν
β
i
-
i
+
1
ν
β
i
+
1
-
μ
r
i
-
1
0,1
β
c
-
1
+
r
i
-
1
1,1
λ
,
i
=
c
-
1
,
c
-
2
,
…
, and then(68)
x
c
=
1
β
0
+
β
1
+
β
2
+
⋯
+
β
c
.Remark 24.
x
0 is proportional to π
0.Remark 25.
Computation ofr
k
(
c
-
1
,
n
) and r
k
(
c
,
n
) using Lemmas 21 and 22 might be numerically unstable due to overflow. Thus, we use recursive formulae in Theorem 26 to obtain a numerically stable scheme.Theorem 26.
Sequence{
x
i
,
i
=
0,1
,
…
,
c
} represents either {
r
i
,
k
c
-
1
,
n
,
i
=
0,1
,
…
,
c
} or {
r
i
,
k
c
,
n
,
i
=
0,1
,
…
,
c
}. {
x
i
,
i
=
0,1
,
…
,
c
-
2
} is calculated in terms of x
c
-
1 and x
c as follows:(69)
x
i
=
i
+
1
ν
x
i
+
1
+
D
i
B
i
,
i
=
0,1
,
…
,
c
-
2
, where {
B
i
,
D
i
,
i
=
0,1
,
…
,
c
-
2
} are given as follows: (70)
B
0
=
λ
+
n
μ
,
D
0
=
0
,
B
i
=
λ
+
i
ν
+
n
μ
-
λ
i
ν
B
i
-
1
,
i
=
1,2
,
…
,
c
-
2
,
D
i
=
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
x
c
-
1
+
r
k
-
1
,
i
-
1
c
,
n
+
1
x
c
+
λ
D
i
-
1
B
i
-
1
,
i
=
1,2
,
…
,
c
-
2
. In addition,(71)
B
i
>
λ
,
D
i
>
0
.Proof.
We prove Theorem26 using mathematical induction. Let {
x
i
,
i
=
0,1
,
…
,
c
} denote {
r
i
,
k
c
-
1
,
n
,
i
=
0,1
,
…
,
c
} defined in (57). We have (72)
x
0
=
ν
λ
+
n
μ
x
1
. Thus, B
0
=
λ
+
n
μ and D
0
=
0. For i
=
1,2
,
…
,
c
-
2, we prove by mathematical induction. For j
=
1,2
,
…
,
i
-
1, assuming that (73)
B
j
=
λ
+
j
ν
+
n
μ
-
λ
j
ν
B
j
-
1
,
D
j
=
n
+
1
μ
r
k
-
1
,
j
-
1
c
-
1
,
n
+
1
x
c
-
1
+
r
k
-
1
,
j
-
1
c
,
n
+
1
x
c
+
λ
D
j
-
1
B
j
-
1 are true, we show that it is also true for j
=
i. Indeed, using the assumption of mathematical induction, we have (74)
λ
i
ν
x
i
+
D
i
-
1
B
i
-
1
-
λ
+
i
ν
+
n
μ
x
i
+
i
+
1
ν
x
i
+
1
+
x
~
i
=
0
, where x
~
i
=
(
n
+
1
)
μ
(
r
k
-
1
,
i
-
1
(
c
-
1
,
n
+
1
)
x
c
-
1
+
r
k
-
1
,
i
-
1
(
c
,
n
+
1
)
x
c
). Arranging this formula yields (75)
x
i
=
i
+
1
ν
x
i
+
1
+
λ
D
i
-
1
/
B
i
-
1
+
x
~
i
λ
+
i
ν
+
n
μ
-
λ
i
ν
/
B
i
-
1
=
i
+
1
ν
x
i
+
1
+
D
i
B
i
, implying that case j
=
i is also true. Thus, for any i
=
1,2
,
…
,
c
-
2, the desired result is established. We can show similar result for {
r
i
,
k
c
,
n
,
i
=
0,1
,
…
,
c
}.Remark 27.
Using Theorem26, we can calculate r
k
,
i
(
c
-
1
,
n
) (i
=
0,1
,
…
,
c
-
2) in terms of r
k
,
c
-
1
(
c
-
1
,
n
) and r
k
,
c
(
c
-
1
,
n
) and r
k
,
i
(
c
,
n
) (i
=
0,1
,
…
,
c
-
2) in terms of r
k
,
c
-
1
(
c
,
n
) and r
k
,
c
(
c
,
n
). Furthermore, r
k
,
c
-
1
(
c
-
1
,
n
), r
k
,
c
(
c
-
1
,
n
), r
k
,
c
-
1
(
c
,
n
), and r
k
,
c
(
c
,
n
) are obtained from Lemmas 21 and 22.Remark 28.
Tran-Gia and Mandjes [6] propose some models where blocked handover calls do not retry but are lost. The results in Section 5.1 are easily adapted to these models.


### 5.2. Computational Algorithm

In this section, we present an algorithm for computing the rate matrices and then a procedure for the computation of the stationary distribution. Algorithm1 shows a method for r
(
n
) while Algorithm 2 computes approximation π
^
=
(
π
^
0
,
π
^
1
,
…
,
π
^
N
) to the stationary distribution, where {
k
l
,
l
∈
Z
+
} is an arbitrary increasing sequence and N is the truncation point given in advance. We will discuss how to choose the truncation point in Section 5.3.Algorithm 1:Computation of r
(
n
).
Input: Q
0
(
n
)
,
Q
1
(
n
)
,
Q
2
(
n
)
,
k
n
;
n
∈
Z
+
,
ϵ
Output: r
^
(
n
)
l
≔
1;
Computer
k
1
(
n
) and r
k
0
(
n
) using Lemmas 21 and 22 and Theorem 26.
while
r
k
l
(
n
)
-
r
k
l
-
1
(
n
)
∞
>
ϵ
do
l
≔
l
+
1;
Compute r
k
l
(
n
) and r
k
l
-
1
(
n
) using Lemmas 21 and 22 and Theorem 26.
r
^
(
n
)
≔
r
k
l
(
n
);
end whileAlgorithm 2:Stationary distribution.
Input: λ
,
μ
,
ν
,
c
,
k
n
;
n
∈
Z
+
,
ϵ
,
N
Output: π
^
n
;
n
=
0,1
,
…
,
N
Computer
^
(
N
) using Algorithm 1.
for
n
=
1 to N
-
1
do
r
^
N
-
n
≔
r
N
-
n
r
^
(
N
-
n
+
1
);
end for
Computex
0 using Lemma 23.
for
n
=
1 to N
do
x
n
≔
x
c
-
1
,
n
-
1
r
^
(
c
-
1
,
n
)
+
x
c
,
n
-
1
r
^
(
c
,
n
);
end for
for
n
=
0 to N
do
π
^
n
≔
x
n
∑
n
=
0
N
x
n
e;
end for


### 5.3. Determination of Truncation PointN

In Algorithm2, the truncation point is given in advance and it should be large enough such that the tail probability is sufficiently small; that is,(76)
∑
n
=
N
+
1
∞
π
n
e
<
ϵ
, where ϵ is given in advance.However, sinceπ
n is not explicitly obtained for general M/M/c/c retrial queues, a direct determination of such N is difficult. In this paper, we use the explicit results for M/M/1/1 retrial queue to determine this truncation point. In particular, we consider M/M/1/1 retrial queue with arrival rate λ
/
c, retrial rate μ, and service rate ν. This queue is stable since ρ
=
λ
/
(
c
ν
)
<
1 due to the stability condition of our original model.Letp
i
,
n
(
i
=
0,1
,
n
∈
Z
+
) denote the probability that the number of busy servers is i and the number of customers in the orbit is n in M/M/1/1 retrial queue. It is shown in [26] that (77)
p
0
,
n
=
ρ
n
n
!
1
-
ρ
λ
/
c
μ
+
1
λ
c
μ
n
,
p
1
,
n
=
ρ
n
+
1
n
!
1
-
ρ
λ
/
c
μ
+
1
1
+
λ
c
μ
n
, where n
∈
Z
+ and (
ϕ
)
n
(
-
∞
<
ϕ
<
∞
,
n
∈
Z
+
) denotes the Pochhammer symbol.Using this result, we set the truncation point as follows:(78)
N
=
inf
⁡
n
∣
∑
i
=
0
n
p
0
,
i
+
p
1
,
i
>
1
-
ϵ
,
ϵ
>
0
.We verify the accuracy of this choice using numerical results.


### 5.4. Blocking Probability

We derive blocking probabilities as performance measures. In our model, priority (type 1) and retrial customers are blocked when all the servers are occupied while normal (type 2) customers are blocked when at leastc
-
1 servers are occupied. Thus, the blocking probability of normal customers is given by(79)
π
c
-
1
+
π
c
≔
∑
n
=
0
∞
π
c
-
1
,
n
+
∑
n
=
0
∞
π
c
,
n
, and the blocking probability of priority and retrial customers is given by(80)
π
c
≔
∑
n
=
0
∞
π
c
,
n
.


## 5.1. Efficient Computation

Due to Lemma4, we need to compute k inverse matrices in order to obtain R
k
(
n
). It may take a long time when the number of servers is large. Thus, instead of computing the inverse matrices, we propose a new method exploiting the fact that only the last two rows are nonzero. The computational complexity of our new method is only O
(
c
). In particular, the computational complexity in all the theorems and lemmas below is O
(
c
).It should be noted that the computation ofR
(
n
) and R
k
(
n
) is equivalent to that of their last two rows r
(
n
) and r
k
(
n
); that is,(57)
r
n
=
r
c
-
1
,
n
r
c
,
n
,
r
k
n
=
r
k
c
-
1
,
n
r
k
c
,
n
,where r
(
i
,
n
) and r
k
(
i
,
n
) (i
=
0,1) are the row vectors of c
+
1 elements.Definition 20.
We define the functionr
n as follows. Let(58)
X
x
,
y
=
O
x
y
,
L
r
Y
=
y
c
-
1
y
c
, where y
c
-
1 and y
c are the second last and the last rows of Y. Furthermore, (59)
r
n
x
y
=
L
r
R
n
X
x
,
y
, where x and y are vectors with an appropriate dimension.It is easy to see thatr
(
n
) and r
k
(
n
) satisfy the following equations: (60)
r
n
=
r
n
r
n
+
1
,
r
k
n
=
r
n
r
k
-
1
n
+
1
=
r
n
∘
r
n
+
1
∘
⋯
∘
r
n
+
k
-
1
O
, for n
,
k
∈
N. Lemmas 21 and 22 compute r
k
(
c
-
1
,
n
) and r
k
(
c
,
n
) using r
k
-
1
(
n
+
1
), respectively. Furthermore, Lemma 23 computes the stationary distribution of the censored Markov chain on level 0 using r
(
1
).Lemma 21.
For arbitraryn, k, one has(61)
r
k
,
i
c
-
1
,
n
=
α
i
+
β
i
r
k
,
c
c
-
1
,
n
,
i
=
0,1
,
…
,
c
-
1
, where {
α
i
,
β
i
,
i
=
0,1
,
…
,
c
} and r
k
,
c
(
c
-
1
,
n
) are given as follows: (62)
α
c
=
0
,
β
c
=
1
,
α
c
-
1
=
0
,
β
c
-
1
=
-
b
c
n
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
,
n
+
1
λ
1
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
-
1
,
n
+
1
,
α
c
-
2
=
-
λ
2
λ
,
β
c
-
2
=
-
b
c
-
1
n
β
c
-
1
+
c
ν
+
n
+
1
μ
r
k
-
1
,
c
-
2
c
-
1
,
n
+
1
β
c
-
1
+
n
+
1
μ
r
k
-
1
,
c
-
2
c
,
n
+
1
λ
,
α
i
-
1
=
-
b
i
n
α
i
+
i
+
1
ν
α
i
+
1
λ
,
i
=
c
-
2
,
c
-
3
,
…
,
1
,
β
i
-
1
=
-
b
i
n
β
i
+
i
+
1
ν
β
i
+
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
β
c
-
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
,
n
+
1
λ
,
i
=
c
-
2
,
c
-
3
,
…
,
1
,
r
k
,
c
c
-
1
,
n
=
-
b
0
n
α
0
+
ν
α
1
b
0
n
β
0
+
ν
β
1
.Proof.
The technical details are provided in AppendixE.Lemma 22.
For arbitraryn and k, one has (63)
r
k
,
i
c
,
n
=
α
i
+
β
i
r
k
,
c
c
,
n
,
i
=
0,1
,
…
,
c
-
1
, where {
α
i
,
β
i
,
i
=
0,1
,
…
,
c
} and r
k
,
c
(
c
,
n
) are given as follows: (64)
α
c
=
0
,
β
c
=
1
,
α
c
-
1
=
-
λ
λ
1
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
-
1
,
n
+
1
,
β
c
-
1
=
-
b
c
n
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
,
n
+
1
λ
1
+
n
+
1
μ
r
k
-
1
,
c
-
1
c
-
1
,
n
+
1
,
α
i
-
1
=
-
b
i
n
α
i
+
i
+
1
ν
α
i
+
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
α
c
-
1
λ
,
i
=
c
-
1
,
c
-
2
,
…
,
1
,
β
i
-
1
=
-
b
i
n
β
i
+
i
+
1
ν
β
i
+
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
β
c
-
1
+
n
+
1
μ
r
k
-
1
,
i
-
1
c
,
n
+
1
λ
,
i
=
c
-
1
,
c
-
2
,
…
,
1
. Furthermore,(65)
r
k
,
c
c
,
n
=
-
b
0
n
α
0
+
ν
α
1
b
0
n
β
0
+
ν
β
1
.Proof.
This lemma can be proved using the same technique as in Lemma21.Lemma 23.
Solutionx
0
=
(
x
0
,
x
1
,
…
,
x
c
) to (66)
x
0
Q
1
0
+
R
1
Q
2
1
=
0
,
x
0
e
=
1
, is given by x
i
=
β
i
x
c
(
i
=
0,1
,
…
,
c
), where {
β
i
,
i
=
0,1
,
…
,
c
-
1
,
c
} is recursively defined as (67)
β
c
=
1
,
β
c
-
1
=
λ
+
c
ν
-
μ
r
c
-
1
1,1
λ
1
+
μ
r
c
-
1
0,1
,
β
i
-
1
=
λ
+
i
ν
β
i
-
i
+
1
ν
β
i
+
1
-
μ
r
i
-
1
0,1
β
c
-
1
+
r
i
-
1
1,1
λ
,
i
=
c
-
1
,
c
-
2
,
…
, and then(68)
x
c
=
1
β
0
+
β
1
+
β
2
+
⋯
+
β
c
.Remark 24.
x
0 is proportional to π
0.Remark 25.
Computation ofr
k
(
c
-
1
,
n
) and r
k
(
c
,
n
) using Lemmas 21 and 22 might be numerically unstable due to overflow. Thus, we use recursive formulae in Theorem 26 to obtain a numerically stable scheme.Theorem 26.
Sequence{
x
i
,
i
=
0,1
,
…
,
c
} represents either {
r
i
,
k
c
-
1
,
n
,
i
=
0,1
,
…
,
c
} or {
r
i
,
k
c
,
n
,
i
=
0,1
,
…
,
c
}. {
x
i
,
i
=
0,1
,
…
,
c
-
2
} is calculated in terms of x
c
-
1 and x
c as follows:(69)
x
i
=
i
+
1
ν
x
i
+
1
+
D
i
B
i
,
i
=
0,1
,
…
,
c
-
2
, where {
B
i
,
D
i
,
i
=
0,1
,
…
,
c
-
2
} are given as follows: (70)
B
0
=
λ
+
n
μ
,
D
0
=
0
,
B
i
=
λ
+
i
ν
+
n
μ
-
λ
i
ν
B
i
-
1
,
i
=
1,2
,
…
,
c
-
2
,
D
i
=
n
+
1
μ
r
k
-
1
,
i
-
1
c
-
1
,
n
+
1
x
c
-
1
+
r
k
-
1
,
i
-
1
c
,
n
+
1
x
c
+
λ
D
i
-
1
B
i
-
1
,
i
=
1,2
,
…
,
c
-
2
. In addition,(71)
B
i
>
λ
,
D
i
>
0
.Proof.
We prove Theorem26 using mathematical induction. Let {
x
i
,
i
=
0,1
,
…
,
c
} denote {
r
i
,
k
c
-
1
,
n
,
i
=
0,1
,
…
,
c
} defined in (57). We have (72)
x
0
=
ν
λ
+
n
μ
x
1
. Thus, B
0
=
λ
+
n
μ and D
0
=
0. For i
=
1,2
,
…
,
c
-
2, we prove by mathematical induction. For j
=
1,2
,
…
,
i
-
1, assuming that (73)
B
j
=
λ
+
j
ν
+
n
μ
-
λ
j
ν
B
j
-
1
,
D
j
=
n
+
1
μ
r
k
-
1
,
j
-
1
c
-
1
,
n
+
1
x
c
-
1
+
r
k
-
1
,
j
-
1
c
,
n
+
1
x
c
+
λ
D
j
-
1
B
j
-
1 are true, we show that it is also true for j
=
i. Indeed, using the assumption of mathematical induction, we have (74)
λ
i
ν
x
i
+
D
i
-
1
B
i
-
1
-
λ
+
i
ν
+
n
μ
x
i
+
i
+
1
ν
x
i
+
1
+
x
~
i
=
0
, where x
~
i
=
(
n
+
1
)
μ
(
r
k
-
1
,
i
-
1
(
c
-
1
,
n
+
1
)
x
c
-
1
+
r
k
-
1
,
i
-
1
(
c
,
n
+
1
)
x
c
). Arranging this formula yields (75)
x
i
=
i
+
1
ν
x
i
+
1
+
λ
D
i
-
1
/
B
i
-
1
+
x
~
i
λ
+
i
ν
+
n
μ
-
λ
i
ν
/
B
i
-
1
=
i
+
1
ν
x
i
+
1
+
D
i
B
i
, implying that case j
=
i is also true. Thus, for any i
=
1,2
,
…
,
c
-
2, the desired result is established. We can show similar result for {
r
i
,
k
c
,
n
,
i
=
0,1
,
…
,
c
}.Remark 27.
Using Theorem26, we can calculate r
k
,
i
(
c
-
1
,
n
) (i
=
0,1
,
…
,
c
-
2) in terms of r
k
,
c
-
1
(
c
-
1
,
n
) and r
k
,
c
(
c
-
1
,
n
) and r
k
,
i
(
c
,
n
) (i
=
0,1
,
…
,
c
-
2) in terms of r
k
,
c
-
1
(
c
,
n
) and r
k
,
c
(
c
,
n
). Furthermore, r
k
,
c
-
1
(
c
-
1
,
n
), r
k
,
c
(
c
-
1
,
n
), r
k
,
c
-
1
(
c
,
n
), and r
k
,
c
(
c
,
n
) are obtained from Lemmas 21 and 22.Remark 28.
Tran-Gia and Mandjes [6] propose some models where blocked handover calls do not retry but are lost. The results in Section 5.1 are easily adapted to these models.


## 5.2. Computational Algorithm

In this section, we present an algorithm for computing the rate matrices and then a procedure for the computation of the stationary distribution. Algorithm1 shows a method for r
(
n
) while Algorithm 2 computes approximation π
^
=
(
π
^
0
,
π
^
1
,
…
,
π
^
N
) to the stationary distribution, where {
k
l
,
l
∈
Z
+
} is an arbitrary increasing sequence and N is the truncation point given in advance. We will discuss how to choose the truncation point in Section 5.3.Algorithm 1:Computation of r
(
n
).
Input: Q
0
(
n
)
,
Q
1
(
n
)
,
Q
2
(
n
)
,
k
n
;
n
∈
Z
+
,
ϵ
Output: r
^
(
n
)
l
≔
1;
Computer
k
1
(
n
) and r
k
0
(
n
) using Lemmas 21 and 22 and Theorem 26.
while
r
k
l
(
n
)
-
r
k
l
-
1
(
n
)
∞
>
ϵ
do
l
≔
l
+
1;
Compute r
k
l
(
n
) and r
k
l
-
1
(
n
) using Lemmas 21 and 22 and Theorem 26.
r
^
(
n
)
≔
r
k
l
(
n
);
end whileAlgorithm 2:Stationary distribution.
Input: λ
,
μ
,
ν
,
c
,
k
n
;
n
∈
Z
+
,
ϵ
,
N
Output: π
^
n
;
n
=
0,1
,
…
,
N
Computer
^
(
N
) using Algorithm 1.
for
n
=
1 to N
-
1
do
r
^
N
-
n
≔
r
N
-
n
r
^
(
N
-
n
+
1
);
end for
Computex
0 using Lemma 23.
for
n
=
1 to N
do
x
n
≔
x
c
-
1
,
n
-
1
r
^
(
c
-
1
,
n
)
+
x
c
,
n
-
1
r
^
(
c
,
n
);
end for
for
n
=
0 to N
do
π
^
n
≔
x
n
∑
n
=
0
N
x
n
e;
end for


## 5.3. Determination of Truncation PointN

In Algorithm2, the truncation point is given in advance and it should be large enough such that the tail probability is sufficiently small; that is,(76)
∑
n
=
N
+
1
∞
π
n
e
<
ϵ
, where ϵ is given in advance.However, sinceπ
n is not explicitly obtained for general M/M/c/c retrial queues, a direct determination of such N is difficult. In this paper, we use the explicit results for M/M/1/1 retrial queue to determine this truncation point. In particular, we consider M/M/1/1 retrial queue with arrival rate λ
/
c, retrial rate μ, and service rate ν. This queue is stable since ρ
=
λ
/
(
c
ν
)
<
1 due to the stability condition of our original model.Letp
i
,
n
(
i
=
0,1
,
n
∈
Z
+
) denote the probability that the number of busy servers is i and the number of customers in the orbit is n in M/M/1/1 retrial queue. It is shown in [26] that (77)
p
0
,
n
=
ρ
n
n
!
1
-
ρ
λ
/
c
μ
+
1
λ
c
μ
n
,
p
1
,
n
=
ρ
n
+
1
n
!
1
-
ρ
λ
/
c
μ
+
1
1
+
λ
c
μ
n
, where n
∈
Z
+ and (
ϕ
)
n
(
-
∞
<
ϕ
<
∞
,
n
∈
Z
+
) denotes the Pochhammer symbol.Using this result, we set the truncation point as follows:(78)
N
=
inf
⁡
n
∣
∑
i
=
0
n
p
0
,
i
+
p
1
,
i
>
1
-
ϵ
,
ϵ
>
0
.We verify the accuracy of this choice using numerical results.


## 5.4. Blocking Probability

We derive blocking probabilities as performance measures. In our model, priority (type 1) and retrial customers are blocked when all the servers are occupied while normal (type 2) customers are blocked when at leastc
-
1 servers are occupied. Thus, the blocking probability of normal customers is given by(79)
π
c
-
1
+
π
c
≔
∑
n
=
0
∞
π
c
-
1
,
n
+
∑
n
=
0
∞
π
c
,
n
, and the blocking probability of priority and retrial customers is given by(80)
π
c
≔
∑
n
=
0
∞
π
c
,
n
.


## 6. Numerical Results

In this section, we show some numerical examples. In particular, in Section6.1 we confirm the effectiveness of the Taylor series expansion for the rate matrices. Section 6.2 is devoted to the numerical investigation of the asymptotic behavior for the joint stationary distribution. Section 6.3 presents the blocking probabilities for priority and normal customers against the number of channels.


### 6.1. Accuracy of the Taylor Series Expansion

The rate matrix is calculated using Algorithm1 with ϵ
=
1
0
-
10 and k
n
=
2
n. We call the rate matrix obtained by Algorithm 1 under this setting exact result.First, we present some numerical examples to show the effectiveness of the Taylor series expansion. Tables1 and 2 show numerical results of R
(
n
) (n
=
100) for μ
=
1 and 10, respectively. Tables 3 and 4 show numerical results of R
(
n
) (n
=
1000) for μ
=
1 and 10, respectively. Other parameters are given by c
=
5, ν
=
1, and λ
2
/
λ
1
=
4 and λ is calculated from traffic intensity ρ
(
=
λ
/
c
ν
). We obtain exact value for the rate matrices using the matrix continued fraction approach [23] with enough accuracy (relative error of the order of 1
0
-
10). The one-, two-, and three-term expansions (m
=
1,2
,
3) are expressed by R
(
n
,
1
), R
(
n
,
2
), and R
(
n
,
3
), respectively. In these tables, we show the relative errors, that is, R
n
,
1
-
R
n
∞
/
R
n
∞, R
n
,
2
-
R
n
∞
/
R
n
∞, and R
n
,
3
-
R
n
∞
/
R
n
∞.Table 1
Relative error forR
(
n
) (n
=
100
,
μ
=
1).
(ρ)
One term
Two terms
Three terms
0.1
0.0051053401
0.0003425140
0.0000228094
0.2
0.0086100661
0.0006446694
0.0000491957
0.3
0.0120849796
0.0009702635
0.0000821267
0.4
0.0155304303
0.0013188638
0.0001219509
0.5
0.0189467632
0.0016900430
0.0001690102
0.6
0.0223343192
0.0020833798
0.0002236397
0.7
0.0256934342
0.0024984580
0.0002861679
0.8
0.0290244403
0.0029348670
0.0003569166
0.9
0.0323276648
0.0033922015
0.0004362009Table 2
Relative error forR
(
n
) (n
=
100
,
μ
=
10).
(ρ)
One term
Two terms
Three terms
0.1
0.0004109454
0.0000030629
0.0000000407
0.2
0.0008055339
0.0000063695
0.0000000804
0.3
0.0011997344
0.0000099832
0.0000001323
0.4
0.0015935474
0.0000139034
0.0000001977
0.5
0.0019869735
0.0000181294
0.0000002747
0.6
0.0023800133
0.0000226607
0.0000003638
0.7
0.0027726674
0.0000274965
0.0000004657
0.8
0.0031649363
0.0000326363
0.0000005807
0.9
0.0035568205
0.0000380794
0.0000007095Table 3
Relative error forR
(
n
) (n
=
1000
,
μ
=
1).
(ρ)
One term
Two terms
Three terms
0.1
0.0004109342
0.0000030754
0.0000000215
0.2
0.0008055116
0.0000063974
0.0000000500
0.3
0.0011997010
0.0000100293
0.0000000863
0.4
0.0015935030
0.0000139704
0.0000001309
0.5
0.0019869182
0.0000182201
0.0000001843
0.6
0.0023799470
0.0000227778
0.0000002472
0.7
0.0027725901
0.0000276429
0.0000003200
0.8
0.0031648480
0.0000328146
0.0000004033
0.9
0.0035567214
0.0000382924
0.0000004976Table 4
Relative error forR
(
n
) (n
=
1000
,
μ
=
10).
(ρ)
One term
Two terms
Three terms
0.1
0.0000401092
0.0000000304
0.0000000000
0.2
0.0000800545
0.0000000640
0.0000000001
0.3
0.0001199958
0.0000001007
0.0000000001
0.4
0.0001599331
0.0000001406
0.0000000002
0.5
0.0001998664
0.0000001837
0.0000000003
0.6
0.0002397957
0.0000002300
0.0000000004
0.7
0.0002797210
0.0000002795
0.0000000005
0.8
0.0003196424
0.0000003321
0.0000000006
0.9
0.0003595598
0.0000003880
0.0000000007We observe that the Taylor series expansion gives a good approximation in the sense that the relative error is quite small. The relative errors for casen
=
1000 are smaller than those for case n
=
100. This fact agrees with the Taylor series expansion formulae. We also observe that the relative error increases with the traffic intensity. This suggests that we need more computational effort for the cases of relatively heavy load in comparison with those of relatively light load.Furthermore, we observe that relative error in Tables2 and 4 (μ
=
10) is smaller than the corresponding one in Tables 1 and 2 (μ
=
1), respectively. This implies that the Taylor series expansion gives good approximation for the case of a relatively large retrial rate. This is the case of interest in practice where customers are impatient.Figure1 represent r
c
(
c
-
1
,
n
) against the number of expansion terms. The parameters are given by n
=
1000, c
=
100, μ
=
1, ν
=
1, λ
2
/
λ
1
=
24, and ρ
=
0.9. We observe that the Taylor series expansion converges to the exact value after about 5 terms. Interestingly, we observe that the Taylor series expansions for r
c
(
c
-
1
,
n
) oscillate and converge to the exact values.Figure 1
r
c
(
c
-
1
,
n
) versus # of expansion terms.


### 6.2. Asymptotic Behavior ofπ
i
,
n
/
ρ
n

Figure2 shows π
i
,
n
/
ρ
n
(
n
∈
Z
+
) against n for some fixed i. Parameters are given by c
=
100, N
=
1000, μ
=
1, ν
=
1
/
70, λ
1
=
1
/
25, and λ
2
=
24
/
25. Joint probability π
i
,
n is computed using the algorithm presented in [23] (see [22] for details). We observe that the five curves for i
=
100,75,50,25, and 0 have a negative slope. This implies that there should exist positive C
1
,
i, C
2
,
i, and b
i such that(81)
C
1
,
i
ρ
n
n
-
b
i
≤
π
i
,
n
≤
C
2
,
i
ρ
n
n
-
b
i
,
n
⟶
∞
. Thus, the asymptotic results obtained in this paper can be further refined to be tighter.Figure 2
π
i
,
n
/
ρ
n versus the number of customers in orbit (n).


### 6.3. Blocking Probability versus Number of Servers

We use the following fixed parametersν
=
1
,
ρ
=
0.7, and λ
2
/
λ
1
=
24 while varying the number of servers from 1 to 100. Truncation point N is determined using the method in Section 5.3 with ϵ
=
1
0
-
10. Blocking probabilities are π
c and π
c
-
1
+
π
c for high and low priority customers, respectively. Figure 3 represents the blocking probabilities of two types of customers for three values of μ (0.1, 1, and 10). Obviously, for the same μ, the blocking probability for low priority customers is higher than that of high priority customers. From Figure 3, we can observe a large difference between the curve for type 1 customers and the corresponding one for type 2 customers. This implies that one guard channel is enough to guarantee the QoS of type 1 customers. Furthermore, the blocking probabilities increase with μ because customers who retry in a short interval may suffer from the same congested situation.Figure 3
Blocking probability versus the number of servers.An important observation is that all the curves are asymptotically linear when the number of servers is large. Asymptotic analysis for the case of large number of servers may be the topic of any future research. In this direction, Avram et al. [27] consider the blocking probability under slow retrials and Halfin-Whitt regime.


## 6.1. Accuracy of the Taylor Series Expansion

The rate matrix is calculated using Algorithm1 with ϵ
=
1
0
-
10 and k
n
=
2
n. We call the rate matrix obtained by Algorithm 1 under this setting exact result.First, we present some numerical examples to show the effectiveness of the Taylor series expansion. Tables1 and 2 show numerical results of R
(
n
) (n
=
100) for μ
=
1 and 10, respectively. Tables 3 and 4 show numerical results of R
(
n
) (n
=
1000) for μ
=
1 and 10, respectively. Other parameters are given by c
=
5, ν
=
1, and λ
2
/
λ
1
=
4 and λ is calculated from traffic intensity ρ
(
=
λ
/
c
ν
). We obtain exact value for the rate matrices using the matrix continued fraction approach [23] with enough accuracy (relative error of the order of 1
0
-
10). The one-, two-, and three-term expansions (m
=
1,2
,
3) are expressed by R
(
n
,
1
), R
(
n
,
2
), and R
(
n
,
3
), respectively. In these tables, we show the relative errors, that is, R
n
,
1
-
R
n
∞
/
R
n
∞, R
n
,
2
-
R
n
∞
/
R
n
∞, and R
n
,
3
-
R
n
∞
/
R
n
∞.Table 1
Relative error forR
(
n
) (n
=
100
,
μ
=
1).
(ρ)
One term
Two terms
Three terms
0.1
0.0051053401
0.0003425140
0.0000228094
0.2
0.0086100661
0.0006446694
0.0000491957
0.3
0.0120849796
0.0009702635
0.0000821267
0.4
0.0155304303
0.0013188638
0.0001219509
0.5
0.0189467632
0.0016900430
0.0001690102
0.6
0.0223343192
0.0020833798
0.0002236397
0.7
0.0256934342
0.0024984580
0.0002861679
0.8
0.0290244403
0.0029348670
0.0003569166
0.9
0.0323276648
0.0033922015
0.0004362009Table 2
Relative error forR
(
n
) (n
=
100
,
μ
=
10).
(ρ)
One term
Two terms
Three terms
0.1
0.0004109454
0.0000030629
0.0000000407
0.2
0.0008055339
0.0000063695
0.0000000804
0.3
0.0011997344
0.0000099832
0.0000001323
0.4
0.0015935474
0.0000139034
0.0000001977
0.5
0.0019869735
0.0000181294
0.0000002747
0.6
0.0023800133
0.0000226607
0.0000003638
0.7
0.0027726674
0.0000274965
0.0000004657
0.8
0.0031649363
0.0000326363
0.0000005807
0.9
0.0035568205
0.0000380794
0.0000007095Table 3
Relative error forR
(
n
) (n
=
1000
,
μ
=
1).
(ρ)
One term
Two terms
Three terms
0.1
0.0004109342
0.0000030754
0.0000000215
0.2
0.0008055116
0.0000063974
0.0000000500
0.3
0.0011997010
0.0000100293
0.0000000863
0.4
0.0015935030
0.0000139704
0.0000001309
0.5
0.0019869182
0.0000182201
0.0000001843
0.6
0.0023799470
0.0000227778
0.0000002472
0.7
0.0027725901
0.0000276429
0.0000003200
0.8
0.0031648480
0.0000328146
0.0000004033
0.9
0.0035567214
0.0000382924
0.0000004976Table 4
Relative error forR
(
n
) (n
=
1000
,
μ
=
10).
(ρ)
One term
Two terms
Three terms
0.1
0.0000401092
0.0000000304
0.0000000000
0.2
0.0000800545
0.0000000640
0.0000000001
0.3
0.0001199958
0.0000001007
0.0000000001
0.4
0.0001599331
0.0000001406
0.0000000002
0.5
0.0001998664
0.0000001837
0.0000000003
0.6
0.0002397957
0.0000002300
0.0000000004
0.7
0.0002797210
0.0000002795
0.0000000005
0.8
0.0003196424
0.0000003321
0.0000000006
0.9
0.0003595598
0.0000003880
0.0000000007We observe that the Taylor series expansion gives a good approximation in the sense that the relative error is quite small. The relative errors for casen
=
1000 are smaller than those for case n
=
100. This fact agrees with the Taylor series expansion formulae. We also observe that the relative error increases with the traffic intensity. This suggests that we need more computational effort for the cases of relatively heavy load in comparison with those of relatively light load.Furthermore, we observe that relative error in Tables2 and 4 (μ
=
10) is smaller than the corresponding one in Tables 1 and 2 (μ
=
1), respectively. This implies that the Taylor series expansion gives good approximation for the case of a relatively large retrial rate. This is the case of interest in practice where customers are impatient.Figure1 represent r
c
(
c
-
1
,
n
) against the number of expansion terms. The parameters are given by n
=
1000, c
=
100, μ
=
1, ν
=
1, λ
2
/
λ
1
=
24, and ρ
=
0.9. We observe that the Taylor series expansion converges to the exact value after about 5 terms. Interestingly, we observe that the Taylor series expansions for r
c
(
c
-
1
,
n
) oscillate and converge to the exact values.Figure 1
r
c
(
c
-
1
,
n
) versus # of expansion terms.


## 6.2. Asymptotic Behavior ofπ
i
,
n
/
ρ
n

Figure2 shows π
i
,
n
/
ρ
n
(
n
∈
Z
+
) against n for some fixed i. Parameters are given by c
=
100, N
=
1000, μ
=
1, ν
=
1
/
70, λ
1
=
1
/
25, and λ
2
=
24
/
25. Joint probability π
i
,
n is computed using the algorithm presented in [23] (see [22] for details). We observe that the five curves for i
=
100,75,50,25, and 0 have a negative slope. This implies that there should exist positive C
1
,
i, C
2
,
i, and b
i such that(81)
C
1
,
i
ρ
n
n
-
b
i
≤
π
i
,
n
≤
C
2
,
i
ρ
n
n
-
b
i
,
n
⟶
∞
. Thus, the asymptotic results obtained in this paper can be further refined to be tighter.Figure 2
π
i
,
n
/
ρ
n versus the number of customers in orbit (n).


## 6.3. Blocking Probability versus Number of Servers

We use the following fixed parametersν
=
1
,
ρ
=
0.7, and λ
2
/
λ
1
=
24 while varying the number of servers from 1 to 100. Truncation point N is determined using the method in Section 5.3 with ϵ
=
1
0
-
10. Blocking probabilities are π
c and π
c
-
1
+
π
c for high and low priority customers, respectively. Figure 3 represents the blocking probabilities of two types of customers for three values of μ (0.1, 1, and 10). Obviously, for the same μ, the blocking probability for low priority customers is higher than that of high priority customers. From Figure 3, we can observe a large difference between the curve for type 1 customers and the corresponding one for type 2 customers. This implies that one guard channel is enough to guarantee the QoS of type 1 customers. Furthermore, the blocking probabilities increase with μ because customers who retry in a short interval may suffer from the same congested situation.Figure 3
Blocking probability versus the number of servers.An important observation is that all the curves are asymptotically linear when the number of servers is large. Asymptotic analysis for the case of large number of servers may be the topic of any future research. In this direction, Avram et al. [27] consider the blocking probability under slow retrials and Halfin-Whitt regime.


## 7. Concluding Remarks

In this paper, we have introduced a new queueing model with a guard channel for retrial and priority customers. The new queueing model is formulated using QBD process which possesses a sparse structure allowing an efficient numerical algorithm and the Taylor series expansion for all the nonzero elements of the rate matrices. We have also derived an asymptotic upper bound for the joint stationary distribution. Numerical results have revealed that the upper bound can be further improved. Future work includes finding the exact asymptotic formulae for the joint stationary distribution.



---

*Source: 7168359-2016-03-03.xml*