# Generalized Truncated Methods for an Efficient Solution of Retrial Systems

Authors: Ma Jose Domenech-Benlloch; Jose Manuel Gimenez-Guzman; Vicent Pla; Jorge Martinez-Bauset; Vicente Casares-Giner
Journal: Mathematical Problems in Engineering
 (2008)

---

## Abstract

We are concerned with the analytic solution of multiserver retrial queues including the impatience phenomenon. As there are not closed-form solutions to these systems, approximate methods are required. We propose two different generalized truncated methods to effectively solve this type of systems. The methods proposed are based on the homogenization of the state space beyond a given number of users in the retrial orbit. We compare the proposed methods with the most well-known methods appeared in the literature in a wide range of scenarios. We conclude that the proposed methods generally outperform previous proposals in terms of accuracy for the most common performance parameters used in retrial systems with a moderated growth in the computational cost.


---

## Body

## 1. Introduction

A common assumption when evaluating the performance of
communication networks is that users that do not obtain an immediate service
leave the system without retrying. However,  due to
the increasing number of customers and network complexity, the customer
behavior, in general, and the retrial phenomenon, in particular, may have a
significant impact on the network performance. We can find numerous examples of
retrial queues in communication networks. In addition to classical telephony
networks [1], the
effect of retrials in wireless networks [2, 3] is usually a more delicate
matter. We can also find the effect of retrials in data transfers along
Internet [4]. However,
it must be noted that retrial queues can also be observed in a wider range of
systems, such as communication protocols or in many queues in daily life where
a customer retries its access to a queue hoping to experience a lower delay in
the access to a free server. Regarding the appearance of retrials in communication
protocols, we can highlight its common use in the performance evaluation of MAC
protocols like CSMA/CD [5]. The retrial phenomenon is also essential in modeling
call centers (see [6]
for a wealth of references on this issue).The modeling of repeated attempts has been the subject
of numerous investigations. Two functional blocks are typically distinguished
in models which consider retrials: a block that accommodates the servers and
possibly a waiting queue, and a block where users that retry are accommodated,
usually called retrial orbit. More concretely, we consider anM/M/C queueing system (following
Kendall's notation M/M/x stands for a queue with exponentially
distributed interarrival and service times and an infinite capacity waiting
room) with retrials where, following the human behavior, users
can abandon the system with certain probability after an unsuccessful
retry, what is usually referred to as impatience. Therefore, the
resulting model under study will have a nonhomogeneous (because the retrial
rate depends on the number of users in the retrial orbit) and infinite state
space. It is known that the classical theory [7] is developed for random
walks on the semistrip {0,…,C}×ℤ+ with infinitesimal transitions subject to
conditions of space homogeneity. When the space-homogeneity condition does not
hold, that is, in the case described in this paper, the problem of calculating
the equilibrium distribution has not been solved beyond approximate methods
[8, 9]. Indeed, if we focus on the case of multiserver
retrial queues, the absence of closed-form solutions for the main performance
characteristics when C>2 [10] can be emphasized. Therefore, it is clear that
in this case it is necessary to resort to approximate methods. These methods
are usually grouped into three categories: approximations, finite truncated
methods, and generalized truncated methods. Although all the mentioned
categories are in fact approximations and, therefore, sometimes it is not clear
which category a method belongs to, the first category is usually devoted to
methods that can be useful in a certain domain of the system parameters or in
special extreme cases (low retrial intensity, high blocking probabilities,
etc.). Both the finite truncated and generalized truncated methods replace the
original infinite state space by a solvable state space, that is, a model where
steady state probabilities can be computed. However, while for the finite
truncated methods the initial state space is replaced by a finite state space,
in generalized truncated methods it is replaced by another infinite but
solvable state space.The main contribution of this paper is the development
of two new generalized truncated methods that are able to effectively solve
retrial systems with user impatience. These novel methods are inspired in the
model proposed by Neuts and Rao [11]. Comparing the proposed methods with the most
well-known approaches appeared in the literature, we conclude that the new
methods generally outperform the previous proposals in terms of accuracy for
the most common performance parameters used in retrial systems and under a wide
range of scenarios. Moreover, we show that their computational cost increases
moderately compared to the simpler methods. We also show that all the
generalized truncated methods clearly outperform the finite truncated methods
in terms of accuracy, as it was predicted in [10].The rest of the paper is structured as follows.
Section2 describes the model of the system under study and Section 3
summarizes previous approaches appeared in the literature. In Section 4 we
describe the two proposed methods, and their performance is evaluated and
compared to previous approaches in Section 5. Final remarks and a summary of
the results are provided in Section 6.


## 2. System Model

In the system under study (shown in Figure1),
users arrive following a Poisson stream with rate λ and contend for access to a system with C servers, requesting an exponentially
distributed service time with rate μ.
Furthermore, we define the load of the system by ρ=λ/(Cμ).
Without loss of generality, we consider that each user occupies one resource
unit. When a new request finds all servers occupied, it joins the retrial orbit
with probability 1,
having considered an infinite capacity for the retrial orbit. After an
exponentially distributed time of rate μr this user retries, being a successful retrial
if it finds a free server. Otherwise, the user leaves the system with
probability Pi or returns to the retrial orbit with
probability (1−Pi),
starting the retrial procedure again.Figure 1
System model.The model considered can be represented as a
bidimensional continuous-time Markov chain (CTMC) whose state space is defined
by(2.1)S:={s=(k,m):k≤C;m∈ℤ+},being k the number of users being served and m the number of users in the retrial orbit,
constituting a level dependent Quasi Birth and Death Process (QBD) [7] whose transition diagram is
depicted in Figure 2. In QBD related literature, the term level refers to
a set of states with the same second coordinate. The infinitesimal generator of
this process has the following infinite block tridiagonal
structure:(2.2)Q=[A1(0)A0(0)000⋯A2(1)A1(1)A0(1)00⋯0A2(2)A1(2)A0(2)0⋯00⋱⋱⋱⋯].Figure 2
Retrial model under study.The main characteristics of this model are its
infinite state space ({0,…,C}×ℤ+) and also its space heterogeneity produced by
the fact that the retrial rate depends on the number of customers in the
retrial orbit.


### 2.1. Computation of Performance Parameters

The most common performance parameter used in retrial
systems is the blocking probability (Pb), which is defined as the probability of the
system being in a state where arrivals are not accepted. Notwithstanding, other
performance parameters can describe the behavior of retrial systems more
accurately. Those performance parameters are the
mean number of users in the retrial orbit  (Nret),
the immediate service probability (Pis), the delayed service probability (Pds), and the nonservice probability (Pns). Pis is defined as the probability of a user being
served in its first attempt, Pds as the probability of obtaining service but
not in the first attempt, and Pns as the probability of leaving the system due
to impatience without having been served. Obviously, it must be met that Pis+Pds+Pns=1.


## 2.1. Computation of Performance Parameters

The most common performance parameter used in retrial
systems is the blocking probability (Pb), which is defined as the probability of the
system being in a state where arrivals are not accepted. Notwithstanding, other
performance parameters can describe the behavior of retrial systems more
accurately. Those performance parameters are the
mean number of users in the retrial orbit  (Nret),
the immediate service probability (Pis), the delayed service probability (Pds), and the nonservice probability (Pns). Pis is defined as the probability of a user being
served in its first attempt, Pds as the probability of obtaining service but
not in the first attempt, and Pns as the probability of leaving the system due
to impatience without having been served. Obviously, it must be met that Pis+Pds+Pns=1.


## 3. Overview of Approximate Methods

As stated in Section1, it is known [10] that neither a closed-form
analytical solution nor a direct algorithmic computation of the steady state
probabilities has been obtained for the model under study, so we must resort to
approximate methods. Those approximations can be classified into three
categories [10].(i)
Approximations: this category includes those
solutions in which the originalintractable model is simplified making some
assumptions that usually suppose a good approximation only in a certain domain
of the system parameters or in special extreme cases.(ii)
Finite truncated methods: these methods consist
in replacing the original infinite state space by a finite one, so the
resulting model is solvable.(iii)
Generalized truncated methods: in this case,
the original infinite state space is replaced by another infinite but
simplified and solvable state space.It must be noted that the difference between the
categoryapproximations and the rest is sometimes not clear because all
the methods are in fact approximations, and usually the accuracy of a method is
better in a certain domain of the problem.


### 3.1. Approximations

An intuitive first approach to solve retrial systems
is to consider those retrials as new arrivals requesting the access to a
server. With this approximation, known asloss method [12, Section 2.8], the system
becomes a simple loss model whose solution can be obtained using Erlang's B
formula with an arrival rate that is the sum of new users and retrials.
Although this solution is very simple, intuitively its main drawback is the
expected inaccuracy when the retrial rate is high.On the other side, for high values of the retrial rate
the system can be approximated by anM/M/C queueing system, so that the blocking
probability can be computed using Erlang's C formula. As the loss method and
the M/M/C approach are useful for low and high retrial
rates, respectively, in
[12, Section 2.8], an interpolation of the above two
methods is proposed to obtain a more accurate
solution. Hereafter, we will refer to this method as Int.In [13], Greenberg and Wolff propose another approximation
(denoted by GW) based on the assumption that the returning customers see time
averages. Obviously, the approximation made in GW is expected to work more accurately
when the retrial rate decreases, because when a customer is blocked and
reattempts quickly, the probability of finding the system busy is high.


### 3.2. Finite Truncated Methods

The first finite truncated method was proposed by
Wilkinson in [14].
This method, denoted by Wil, is based on the truncation of the state
space of the QBD beyond a level Q,
that is, the method restricts the maximum number of users in the retrial orbit
to Q.
Obviously, the method is expected to be more accurate as we use larger values
for the truncation level Q,
but its computational cost will also increase. Other authors present more
efficient truncated methods, for example, in [15] the truncation is based on
the exclusion of those states with negligible stationary probabilities.There also exist other finite truncated methods that
modify the truncated state space to introduce, in some sense, the effect of
retrials. The model presented by Fredericks and Reisner [16], and called hereafter FR,
reduces the dimensionality of the model to a one dimensional state space,
eliminating the dimension corresponding to the number of users in the retrial
orbit and introducing its effect as a new arrival rate that depends on the
state of the system.In a similar way, Marsan et al. [17] propose a method (denoted
by Mar) that reduces the infinite state space to a finite one grouping all
the levels of the QBD where there are users in the retrial orbit into a single
level. In [18] we
propose a generalized version of Mar method, called FM. Comparing FM and Mar,
FM shows a substantial improvement in the accuracy at the expense of only a
marginal increase in the computational cost. These two methods include the
effect of the eliminated states by modifying the transition rates of the last
level of the resulting QBD. The difference between both methods lies in the
number of levels taken into account in the finite truncated model. While Mar
groups all the states with retrials into a single level, FM defines a value Q so the aggregation is performed beyond level Q.
Mar can be seen as a particular case of FM, where Q=1.
Note also that by increasing the value of Q both the accuracy and the computational cost
will also increase. Due to that aggregation two new parameters are introduced: Θ and p.
The parameter Θ denotes the mean number of users in the
retrial orbit conditioned to those states where there are at least Q users in the orbit, that is, Θ=E(m∣m≥Q).
The probability that after a successful retrial the number of users in the
retrial orbit does not drop below Q is represented by p.
The expressions needed to compute those parameters are [18](3.1)p=π(C,Q)π(C,Q−1)+π(C,Q),Θ=λ(π(C,Q−1)+π(C,Q))μr[∑k=0k=C−1π(k,Q)+Piπ(C,Q)].As there is a mutual dependence betweenΘ and p and the steady state probabilities (π(k,m)), to find these values we followed a
fixed-point iterative procedure starting with p=0 and Θ=Q.


### 3.3. Generalized Truncated Methods

In this section, we discuss the use of generalized truncated methods which approximate theM/M/C retrial queue by some infinite although
solvable system. We comment separately on the method proposed by Neuts and Rao in [11] as the methods proposed in
this paper are based on it.The method introduced by Falin in [19] (denoted by Fal) is
based on the assumption that the retrial rate becomes infinite when the number
of customers in orbit exceeds a certain level Q.
Therefore, the system is approximated beyond level Q by a standard M/M/1 queue, and the retrial rate will depend on the
number of users in the retrial orbit as follows:(3.2)μr(m)={mμrif⁢m<Q,∞if⁢m≥Q.Artalejo and Pozo [10] proposed a method (denoted by AP) based on the M/M/2 queue, instead of the M/M/1 queue. Therefore, the retrial rate will depend
not only on the number of users in the retrial orbit (m) but also on the number of busy servers (k) as follows:(3.3)μr(m)={∞if⁢k<C−1⁢,m≥Q,mμrotherwise.


#### 3.3.1. NR Method

The method proposed by Neuts and Rao in [11], and called hereafter NR, which was proved to
converge to the original model in [20], is based on the homogenization of the model beyond a
given level Q of the QBD, which supposes to restrict the
maximum retrial rate, that is,(3.4)μr(m)={mμrif⁢m<Q,Qμrif⁢m≥Q.Figure3 depicts the system model used in NR method.
The resulting approximated model can be solved by the matrix-geometric approach
[7].
Let(3.5)Q^=[A1(0)A0(0)0⋯00A2(1)A1(1)A0(1)⋯0⋮⋮⋮⋯⋮⋮000⋯A1(Q−1)A0(Q−1)000⋯A2A1+RA2],where R is the rate matrix defined as the unique minimal
nonnegative solution of the quadratic matrix equation R2A2+RA1+A0=0.
In order to compute matrix R, we have used the logarithmic reduction
algorithm as proposed in
[9, Section 8.4], using a precision of 10−8 for the iterative procedure. The steady state
probabilities for states where m≤Q  are obtained by solving [π(0)⋯π(Q)]Q^=0 combined with the normalization
condition(3.6)∑l=0Q−1π(l)e+π(Q)(I−R)−1e=1.Figure 3
NR model.AsQ^ is a finite matrix, this system can be solved
by any of the standard methods defined in the classical linear algebra.
Finally, the steady state probabilities of the states in the homogeneous part
of the model are easily computed by π(Q+n)=π(Q)Rn.


## 3.1. Approximations

An intuitive first approach to solve retrial systems
is to consider those retrials as new arrivals requesting the access to a
server. With this approximation, known asloss method [12, Section 2.8], the system
becomes a simple loss model whose solution can be obtained using Erlang's B
formula with an arrival rate that is the sum of new users and retrials.
Although this solution is very simple, intuitively its main drawback is the
expected inaccuracy when the retrial rate is high.On the other side, for high values of the retrial rate
the system can be approximated by anM/M/C queueing system, so that the blocking
probability can be computed using Erlang's C formula. As the loss method and
the M/M/C approach are useful for low and high retrial
rates, respectively, in
[12, Section 2.8], an interpolation of the above two
methods is proposed to obtain a more accurate
solution. Hereafter, we will refer to this method as Int.In [13], Greenberg and Wolff propose another approximation
(denoted by GW) based on the assumption that the returning customers see time
averages. Obviously, the approximation made in GW is expected to work more accurately
when the retrial rate decreases, because when a customer is blocked and
reattempts quickly, the probability of finding the system busy is high.


## 3.2. Finite Truncated Methods

The first finite truncated method was proposed by
Wilkinson in [14].
This method, denoted by Wil, is based on the truncation of the state
space of the QBD beyond a level Q,
that is, the method restricts the maximum number of users in the retrial orbit
to Q.
Obviously, the method is expected to be more accurate as we use larger values
for the truncation level Q,
but its computational cost will also increase. Other authors present more
efficient truncated methods, for example, in [15] the truncation is based on
the exclusion of those states with negligible stationary probabilities.There also exist other finite truncated methods that
modify the truncated state space to introduce, in some sense, the effect of
retrials. The model presented by Fredericks and Reisner [16], and called hereafter FR,
reduces the dimensionality of the model to a one dimensional state space,
eliminating the dimension corresponding to the number of users in the retrial
orbit and introducing its effect as a new arrival rate that depends on the
state of the system.In a similar way, Marsan et al. [17] propose a method (denoted
by Mar) that reduces the infinite state space to a finite one grouping all
the levels of the QBD where there are users in the retrial orbit into a single
level. In [18] we
propose a generalized version of Mar method, called FM. Comparing FM and Mar,
FM shows a substantial improvement in the accuracy at the expense of only a
marginal increase in the computational cost. These two methods include the
effect of the eliminated states by modifying the transition rates of the last
level of the resulting QBD. The difference between both methods lies in the
number of levels taken into account in the finite truncated model. While Mar
groups all the states with retrials into a single level, FM defines a value Q so the aggregation is performed beyond level Q.
Mar can be seen as a particular case of FM, where Q=1.
Note also that by increasing the value of Q both the accuracy and the computational cost
will also increase. Due to that aggregation two new parameters are introduced: Θ and p.
The parameter Θ denotes the mean number of users in the
retrial orbit conditioned to those states where there are at least Q users in the orbit, that is, Θ=E(m∣m≥Q).
The probability that after a successful retrial the number of users in the
retrial orbit does not drop below Q is represented by p.
The expressions needed to compute those parameters are [18](3.1)p=π(C,Q)π(C,Q−1)+π(C,Q),Θ=λ(π(C,Q−1)+π(C,Q))μr[∑k=0k=C−1π(k,Q)+Piπ(C,Q)].As there is a mutual dependence betweenΘ and p and the steady state probabilities (π(k,m)), to find these values we followed a
fixed-point iterative procedure starting with p=0 and Θ=Q.


## 3.3. Generalized Truncated Methods

In this section, we discuss the use of generalized truncated methods which approximate theM/M/C retrial queue by some infinite although
solvable system. We comment separately on the method proposed by Neuts and Rao in [11] as the methods proposed in
this paper are based on it.The method introduced by Falin in [19] (denoted by Fal) is
based on the assumption that the retrial rate becomes infinite when the number
of customers in orbit exceeds a certain level Q.
Therefore, the system is approximated beyond level Q by a standard M/M/1 queue, and the retrial rate will depend on the
number of users in the retrial orbit as follows:(3.2)μr(m)={mμrif⁢m<Q,∞if⁢m≥Q.Artalejo and Pozo [10] proposed a method (denoted by AP) based on the M/M/2 queue, instead of the M/M/1 queue. Therefore, the retrial rate will depend
not only on the number of users in the retrial orbit (m) but also on the number of busy servers (k) as follows:(3.3)μr(m)={∞if⁢k<C−1⁢,m≥Q,mμrotherwise.


### 3.3.1. NR Method

The method proposed by Neuts and Rao in [11], and called hereafter NR, which was proved to
converge to the original model in [20], is based on the homogenization of the model beyond a
given level Q of the QBD, which supposes to restrict the
maximum retrial rate, that is,(3.4)μr(m)={mμrif⁢m<Q,Qμrif⁢m≥Q.Figure3 depicts the system model used in NR method.
The resulting approximated model can be solved by the matrix-geometric approach
[7].
Let(3.5)Q^=[A1(0)A0(0)0⋯00A2(1)A1(1)A0(1)⋯0⋮⋮⋮⋯⋮⋮000⋯A1(Q−1)A0(Q−1)000⋯A2A1+RA2],where R is the rate matrix defined as the unique minimal
nonnegative solution of the quadratic matrix equation R2A2+RA1+A0=0.
In order to compute matrix R, we have used the logarithmic reduction
algorithm as proposed in
[9, Section 8.4], using a precision of 10−8 for the iterative procedure. The steady state
probabilities for states where m≤Q  are obtained by solving [π(0)⋯π(Q)]Q^=0 combined with the normalization
condition(3.6)∑l=0Q−1π(l)e+π(Q)(I−R)−1e=1.Figure 3
NR model.AsQ^ is a finite matrix, this system can be solved
by any of the standard methods defined in the classical linear algebra.
Finally, the steady state probabilities of the states in the homogeneous part
of the model are easily computed by π(Q+n)=π(Q)Rn.


## 3.3.1. NR Method

The method proposed by Neuts and Rao in [11], and called hereafter NR, which was proved to
converge to the original model in [20], is based on the homogenization of the model beyond a
given level Q of the QBD, which supposes to restrict the
maximum retrial rate, that is,(3.4)μr(m)={mμrif⁢m<Q,Qμrif⁢m≥Q.Figure3 depicts the system model used in NR method.
The resulting approximated model can be solved by the matrix-geometric approach
[7].
Let(3.5)Q^=[A1(0)A0(0)0⋯00A2(1)A1(1)A0(1)⋯0⋮⋮⋮⋯⋮⋮000⋯A1(Q−1)A0(Q−1)000⋯A2A1+RA2],where R is the rate matrix defined as the unique minimal
nonnegative solution of the quadratic matrix equation R2A2+RA1+A0=0.
In order to compute matrix R, we have used the logarithmic reduction
algorithm as proposed in
[9, Section 8.4], using a precision of 10−8 for the iterative procedure. The steady state
probabilities for states where m≤Q  are obtained by solving [π(0)⋯π(Q)]Q^=0 combined with the normalization
condition(3.6)∑l=0Q−1π(l)e+π(Q)(I−R)−1e=1.Figure 3
NR model.AsQ^ is a finite matrix, this system can be solved
by any of the standard methods defined in the classical linear algebra.
Finally, the steady state probabilities of the states in the homogeneous part
of the model are easily computed by π(Q+n)=π(Q)Rn.


## 4. Proposed Methods

### 4.1. HM1 Method

When the number of users retrying is higher thanQ,
the NR method considers the retrial rate to be Qμr.
Obviously, this is a rough approximation that can be inaccurate in many cases.
Therefore, we propose to approximate the retrial rate beyond level Q by an estimation of its mean value, that is, M=E[m∣m≥Q],
keeping the idea of the homogenization:(4.1)μr(m)={mμrif⁢m<Q,Mμrif⁢m≥Q.
The physical meaning of
the proposed approximation is founded on assigning to the retrial rate beyond
level Q a value that is the mean number of users
retrying, when there are at least Q users retrying. This value can be computed by
using(4.2)M=E[m∣m≥Q]=∑r≥Qrπ(r)e∑r≥Qπ(r)e=π(Q)[R(I−R)−1+QI](I−R)−1eπ(Q)(I−R)−1e.However, it must be
noticed thatM depends on the value of Q,
which is a configuration parameter, and more importantly, on the steady state
probabilities. Therefore, we obtain a set of nonlinear equations relating M, R, and the probability vectors π(0),…,π(Q).To solve this system, we use an approximation based on
the assumption that when the number of customers in the orbit is sufficiently
high then it is very likely that all servers are being used. Hence,π(Q)≈π(C,Q)ψ,
where ψ=[00⋯01]t,
and thus M and R no longer depend on π(Q):(4.3)M≈ψ[R(I−R)−1+QI](I−R)−1eψ(I−R)−1e.However, as the relationship betweenR and M still holds, we deploy an iterative procedure
as follows: (1) set M0=Q;
(2) given Mn compute R using the logarithmic reduction algorithm,
then compute Mn+1 using (4.3); (3) if |Mn+1−Mn|/Mn≥ε=10−3 go to step 2;  (4) with the final values of M and R compute the steady state probabilities solving
the finite QBD.


### 4.2. HM2 Method

We must notice that the approximation that performs
HM1 (π(Q)≈π(C,Q)ψ) may not be very accurate when the blocking
probabilty is low, although it
is expected to outperform NR. For that reason, we propose the
HM2 method, which computes the value of M using (4.2). More specifically, HM2 method uses
the iterative procedure sketched next: (1) set M0=Q;
(2) given Mn compute R and the steady state probabilities solving the
finite QBD, then compute Mn+1 using (4.2); (3) if |Mn+1−Mn|/Mn≥ε=10−3 go to step 2.Obviously, the performance of the HM2 method is
expected to be higher than that of the HM1 method, but the computational cost of
HM2 is also expected to be higher due to the fact that in HM1 the values ofπ(0),…,π(Q) are computed only once, and this is not true
in HM2.


### 4.3. Computation of Performance Parameters

For the methods HM1 and HM2, we have used the next
expressions to compute the performance parameters described in
Section2.1.
(4.4)Pb=∑m=0Q−1π(m)z+π(Q)(I−R)−1zwithz=[0,0,…,0,1],Pis=1−Pb,Pds=λ−1μr[∑m=0Q−1mπ(m)o+Mπ(Q)(I−R)−1o]witho=[1,1,…,1,0],Pns=λ−1Piμr[∑m=0Q−1mπ(m)z+Mπ(Q)(I−R)−1z]withz=[0,0,…,0,1],Nret=∑m=0Q−1mπ(m)e+π(Q)(R(I−R)−1+QI)(I−R)−1ewithe=[1,1,…,1].


## 4.1. HM1 Method

When the number of users retrying is higher thanQ,
the NR method considers the retrial rate to be Qμr.
Obviously, this is a rough approximation that can be inaccurate in many cases.
Therefore, we propose to approximate the retrial rate beyond level Q by an estimation of its mean value, that is, M=E[m∣m≥Q],
keeping the idea of the homogenization:(4.1)μr(m)={mμrif⁢m<Q,Mμrif⁢m≥Q.
The physical meaning of
the proposed approximation is founded on assigning to the retrial rate beyond
level Q a value that is the mean number of users
retrying, when there are at least Q users retrying. This value can be computed by
using(4.2)M=E[m∣m≥Q]=∑r≥Qrπ(r)e∑r≥Qπ(r)e=π(Q)[R(I−R)−1+QI](I−R)−1eπ(Q)(I−R)−1e.However, it must be
noticed thatM depends on the value of Q,
which is a configuration parameter, and more importantly, on the steady state
probabilities. Therefore, we obtain a set of nonlinear equations relating M, R, and the probability vectors π(0),…,π(Q).To solve this system, we use an approximation based on
the assumption that when the number of customers in the orbit is sufficiently
high then it is very likely that all servers are being used. Hence,π(Q)≈π(C,Q)ψ,
where ψ=[00⋯01]t,
and thus M and R no longer depend on π(Q):(4.3)M≈ψ[R(I−R)−1+QI](I−R)−1eψ(I−R)−1e.However, as the relationship betweenR and M still holds, we deploy an iterative procedure
as follows: (1) set M0=Q;
(2) given Mn compute R using the logarithmic reduction algorithm,
then compute Mn+1 using (4.3); (3) if |Mn+1−Mn|/Mn≥ε=10−3 go to step 2;  (4) with the final values of M and R compute the steady state probabilities solving
the finite QBD.


## 4.2. HM2 Method

We must notice that the approximation that performs
HM1 (π(Q)≈π(C,Q)ψ) may not be very accurate when the blocking
probabilty is low, although it
is expected to outperform NR. For that reason, we propose the
HM2 method, which computes the value of M using (4.2). More specifically, HM2 method uses
the iterative procedure sketched next: (1) set M0=Q;
(2) given Mn compute R and the steady state probabilities solving the
finite QBD, then compute Mn+1 using (4.2); (3) if |Mn+1−Mn|/Mn≥ε=10−3 go to step 2.Obviously, the performance of the HM2 method is
expected to be higher than that of the HM1 method, but the computational cost of
HM2 is also expected to be higher due to the fact that in HM1 the values ofπ(0),…,π(Q) are computed only once, and this is not true
in HM2.


## 4.3. Computation of Performance Parameters

For the methods HM1 and HM2, we have used the next
expressions to compute the performance parameters described in
Section2.1.
(4.4)Pb=∑m=0Q−1π(m)z+π(Q)(I−R)−1zwithz=[0,0,…,0,1],Pis=1−Pb,Pds=λ−1μr[∑m=0Q−1mπ(m)o+Mπ(Q)(I−R)−1o]witho=[1,1,…,1,0],Pns=λ−1Piμr[∑m=0Q−1mπ(m)z+Mπ(Q)(I−R)−1z]withz=[0,0,…,0,1],Nret=∑m=0Q−1mπ(m)e+π(Q)(R(I−R)−1+QI)(I−R)−1ewithe=[1,1,…,1].


## 5. Results

This section is aimed to compare the proposed methods
HM1 and HM2 with the most well-known approaches appeared in the literature so
far. In order to evaluate and compare the proposed algorithms with previous approaches,
we have studied a wide range of scenarios. Lettingρ=λ/(Cμ),
we have studied different system loads by modifying λ and keeping C=50 and μ−1=180s fixed. It must be noted that, due to the
introduction of the impatience phenomenon, we will be able to consider values
of ρ>1.
We have also studied different configurations of the retrial behavior by
varying μr using in all cases Pi=0.2.
Note that it is not necessary to modify Pi to have different loads in the retrial orbit,
as we use μr for that purpose.The comparison of HM1 and HM2 is done with FM and NR
methods, which can be considered representative examples of a finite truncated
method and a generalized truncated method, respectively. It must be noted that
the resulting QBD has been solved using the algorithm proposed in [21] in all cases, although
other algorithms have been proposed in the literature [22, 23]. For the iterative
procedures in FM, HM1, and HM2 the accuracy has been chosen in all cases to be ε=10−3.For comparing the results, we have used the relative
error of different performance parameters, defined as|Ψapprox−Ψexact|/Ψexact for a generic performance parameter Ψ.
In order to obtain an accurate enough estimate of Ψ which can be used as Ψexact,
we ran all methods with increasing and sufficiently high values of Q so that the value of Ψ had stabilized up to the 14th decimal digit.
As an example and for the particular case ρ=0.8 and μr=0.01,
Figure 4 depicts that the general trend is that the higher the value of Q the lower the relative error is. This is due
to the fact that the system under consideration becomes more similar to the
exact model as Q increases. As expected all methods converged
to the same value in the performance parameters under study, Ψ∈{Pb,Pds,Pns,Nret}.
In Table 1 we show the values of Pb and Nret that are going to be considered as the exact
values, obtained as described above. Obviously, to obtain those values we need
higher values of Q as Pb and Nret increase. In general, when we use values of ρ>1,
the value of Q needed to obtain a tight estimation of Ψ increases significantly, so the computational
complexity makes the resolution to be unfeasible from a practical point of view
when ρ>10.
For that reason and since blocking probabilities between 50% and 75% are unacceptable for virtually any practical
application, the comparison among different methods is done in a range of 0.4≤ρ≤1.4.Table 1
Estimations ofPb and Nret.
ρμr=0.001μr=0.01μr=0.1μr=1.00.4Pb=7.664⋅10−9Pb=7.951⋅10−9Pb=9.299⋅10−9Pb=9.360⋅10−9Nret=8.558⋅10−7Nret=9.276⋅10−8Nret=1.487⋅10−8Nret=3.017⋅10−90.6Pb=2.262⋅10−4Pb=2.566⋅10−4Pb=3.378⋅10−4Pb=3.050⋅10−4Nret=3.800⋅10−2Nret=4.627⋅10−3Nret=8.947⋅10−4Nret=1.531⋅10−40.8Pb=2.548⋅10−2Pb=3.318⋅10−2Pb=3.994⋅10−2Pb=2.897⋅10−2Nret=5.883Nret=0.876Nret=0.160Nret=2.009⋅10−21.0Pb=0.347Pb=0.333Pb=0.273Pb=0.172Nret=1.359⋅102Nret=14.332Nret=1.573Nret=0.1541.2Pb=0.6482Pb=0.6387Pb=0.5436Pb=0.3541Nret=450.19Nret=44.74Nret=4.300Nret=0.39281.4Pb=0.7563Pb=0.7524Pb=0.7028Pb=0.5032Nret=745.68Nret=74.44Nret=7.290Nret=0.67002.0Pb=0.8713Pb=0.8706Pb=0.8619Pb=0.7527Nret=1598.46Nret=159.81Nret=15.935Nret=1.5385.0Pb=0.9613Pb=0.9612Pb=0.9607Pb=0.9539Nret=5780.45Nret=578.04Nret=57.80Nret=5.77010.0Pb=0.9821Pb=0.9821Pb=0.9820Pb=0.9809Nret=12728.44Nret=1272.84Nret=127.28Nret=12.725Evolution of the relative error for the different methods deployed.
(a)
Pb(b)
NretThe different methods are compared by using a metric
based on the minimum truncation level needed to obtain a certain relative
error, which has been widely used in the literature [10, 11]. Table 2 shows the minimum
value of Q needed to obtain a relative error lower than 10−4 for Pb, Pds, Pns and Nret.
Note that a wide range of operation points for the occupancy of the servers and
the retrial orbit has been chosen. The number in bold indicates the lowest
value of Q for all the models studied. Finally, it is
important to note that immediate service probability has been omitted as it is
the complementary value of Pb in the system under study.Table 2
Minimum value ofQ to obtain relative errors lower than 10−4.
μr=0.001μr=0.01μr=0.1μr=1.0ρPbNretPdsPnsPbNretPdsPnsPbNretPdsPnsPbNretPdsPns0.4FM6861481071358894454NR594117961147684544HM1594117891056783345HM212264538454633440.6FM1718162715181322101111125555NR11161120131412188109104555HM110151720111315185810114456HM24421368412675844440.8FM5661537630323237131214146465NR4349425821292432121312125655HM13339505720222832101113144566HM216231839192019229991044541.0FM26425026626658545958171517166575NR22625423023650565150151715146765HM113721124925641465457131316175567HM2144190182197383540401210121154541.2FM56655257153587838877201821177675NR51656452348776867868182018156765HM145349856056263738486151619205577HM2454500499500626263641411151354631.4FM85683786479111511011899232124187685NR79286180375010311710587202421166875HM174879285685884101115117171923246678HM274879579479583919291171418155564Results show that the best performance is obtained by
HM2, followed by HM1, NR, and FM. This trend is obtained in almost every
scenario studied and for all the performance parameters. It is not unexpected
that HM2 is more accurate than NR and HM1 as these can be considered as
approximations of HM2. As it could be expected generalized truncated methods
outperform finite truncated methods, as predicted in [12]. More specifically, the NR
method improves FM, as it requires a lower value of Q for obtaining a concrete relative error. HM1
method obtains better accuracy than NR method for the same value of Q for Pb and Nret,
but this is not true for Pds and Pns,
in which the inverse behavior is obtained. Finally, HM2 method clearly
outperforms all methods, as it needs much lower values of Q to obtain the desired accuracy for all the
performance parameters under study.Nonetheless, as we are dealing with numerical methods
some attention must be paid to the computational cost. Note that the same
methodology for solving the QBD and the same precision for the iterative
procedures have been used in all cases. In Figure5 we show the time
needed to solve the system when ρ=0.8 and μr=10−2 for different values of Q. Results have been obtained running Matlab in an
Intel Pentium Core 2.
As observed, the generalized truncated methods are not only more accurate than
FM, but also they are computationally more efficient. However, we must note
that times are negligible from a human point of view in all the generalized
truncated methods. If we compare the results with those obtained in Figure 4, we
conclude that we can obtain high accuracies with small computation times,
specially for the generalized truncated methods.Figure 5
Computation time for the different methods deployed.


### 5.1. Scenario with no Impatience

Although the impatience phenomenon is part of the
human behavior, some of the most well-known methods to solve retrial systems do
not consider it. For that reason, in this section we consider the particular
case of the proposed model in whichPi=0,
that is, we do not take into account user impatience. In this way we will be
able to make a global comparison of all methods. Note that in this scenario,
the load offered to the system (ρ) has to be less than one. Moreover, we have
decreased the number of servers to C=10,
because the computational complexity of AP makes it infeasible to solve the
system when C≥50 servers. We have focused only in the blocking
probability because when we do not consider the impatience phenomenon Pns=0 and, therefore, Pds=Pb.
We can classify all the methods described in Section 3 into two categories:
those in which the truncation level Q is a configurable parameter and those that
offer a unique solution. In the former the precision can be adjusted through
the parameter Q while in the latter precision is an intrinsic
feature of the method and hence it is fixed. In the first group we find Wil,
FM, Fal, NR, AP, HM1, and HM2 while in the latter we find FR, GW, Loss, Int and
Mar. In a first part, we study the performance of the latter, showing Pb and the relative error for each method in
Table 3. As we can see, errors are usually lower when we increase the time
between retrials but, in any case, errors are unacceptable in most cases. As
Mar is the best method in all of them and it is a particular case of FM in the
sequel only the FM method will be considered as it yields an upper bound for
the accuracy of fixed precision methods.Table 3
Pb and its relative error obtained with fixed
precision methods.
ρμr=0.001μr=0.01μr=0.10.6FR1.3987⋅10−25.4617⋅10−24.6031⋅10−2GW4.3889⋅10−20.21870.4022Loss3.4677⋅10−20.21120.3964Int3.4484⋅10−20.20960.3846Mar8.6043⋅10−48.385⋅10−32.066⋅10−2Pb =5.683⋅10−2Pb =6.954⋅10−2Pb =9.089⋅10−20.8FR2.8113⋅10−20.10227.9129⋅10−2GW4.3607⋅10−20.20790.3702Loss4.1717⋅10−20.20640.3690Int4.1451⋅10−20.20420.3521Mar1.1167⋅10−22.3249⋅10−22.3980⋅10−2Pb=0.2469Pb=0.2982Pb=0.3750Focusing on the configurable methods, Table4 shows
the minimum value of Q needed to obtain a relative error in the
blocking probability lower than 10−4.
Note that the number in bold represents the best choice for each scenario. As
we can see, in scenarios without impatience, generalized truncated methods
(Fal, NR, AP, HM1, and HM2) also outperform finite truncated methods (Wil and
FM). Among the generalized truncated methods, the best performance is obtained
by HM2 and AP. While the first is specially useful for low values of μr,
the latter is the best for high values of μr (μr/μ>10).Table 4
Minimum value ofQ to obtain relative errors for Pb lower than 10−4 when Pi=0.
ρμr=0.001μr=0.01μr=0.1μr=1.00.4Wil9101010FM9974Fal11974NR7754AP10741HM17653HM23543Pb =5.633⋅10−3Pb= 6.442⋅10−3Pb =7.998⋅10−3Pb =8.696⋅10−30.6Wil24191817FM2316115Fal2516117NR181396AP221361HM1151074HM25974Pb =5.683⋅10−2Pb =6.954⋅10−2Pb =9.089⋅10−2Pb =9.979⋅10−20.8Wil74454039FM68341721Fal70352314NR53271710AP572491HM14121137HM23620137Pb=0.2469Pb=0.2982Pb=0.3750Pb=0.4043In Figure6 we show the computational cost
associated to all the presented methods when ρ=0.8 and μr=0.01 for different values of Q.
Results show that the lowest computation costs are for the simplest methods
(Wil, Fal, NR) followed by HM1 and HM2. Finally, and computationally speaking,
the worst methods are FM and AP.Figure 6
Computation time for the different methods deployed.Pi=0.


## 5.1. Scenario with no Impatience

Although the impatience phenomenon is part of the
human behavior, some of the most well-known methods to solve retrial systems do
not consider it. For that reason, in this section we consider the particular
case of the proposed model in whichPi=0,
that is, we do not take into account user impatience. In this way we will be
able to make a global comparison of all methods. Note that in this scenario,
the load offered to the system (ρ) has to be less than one. Moreover, we have
decreased the number of servers to C=10,
because the computational complexity of AP makes it infeasible to solve the
system when C≥50 servers. We have focused only in the blocking
probability because when we do not consider the impatience phenomenon Pns=0 and, therefore, Pds=Pb.
We can classify all the methods described in Section 3 into two categories:
those in which the truncation level Q is a configurable parameter and those that
offer a unique solution. In the former the precision can be adjusted through
the parameter Q while in the latter precision is an intrinsic
feature of the method and hence it is fixed. In the first group we find Wil,
FM, Fal, NR, AP, HM1, and HM2 while in the latter we find FR, GW, Loss, Int and
Mar. In a first part, we study the performance of the latter, showing Pb and the relative error for each method in
Table 3. As we can see, errors are usually lower when we increase the time
between retrials but, in any case, errors are unacceptable in most cases. As
Mar is the best method in all of them and it is a particular case of FM in the
sequel only the FM method will be considered as it yields an upper bound for
the accuracy of fixed precision methods.Table 3
Pb and its relative error obtained with fixed
precision methods.
ρμr=0.001μr=0.01μr=0.10.6FR1.3987⋅10−25.4617⋅10−24.6031⋅10−2GW4.3889⋅10−20.21870.4022Loss3.4677⋅10−20.21120.3964Int3.4484⋅10−20.20960.3846Mar8.6043⋅10−48.385⋅10−32.066⋅10−2Pb =5.683⋅10−2Pb =6.954⋅10−2Pb =9.089⋅10−20.8FR2.8113⋅10−20.10227.9129⋅10−2GW4.3607⋅10−20.20790.3702Loss4.1717⋅10−20.20640.3690Int4.1451⋅10−20.20420.3521Mar1.1167⋅10−22.3249⋅10−22.3980⋅10−2Pb=0.2469Pb=0.2982Pb=0.3750Focusing on the configurable methods, Table4 shows
the minimum value of Q needed to obtain a relative error in the
blocking probability lower than 10−4.
Note that the number in bold represents the best choice for each scenario. As
we can see, in scenarios without impatience, generalized truncated methods
(Fal, NR, AP, HM1, and HM2) also outperform finite truncated methods (Wil and
FM). Among the generalized truncated methods, the best performance is obtained
by HM2 and AP. While the first is specially useful for low values of μr,
the latter is the best for high values of μr (μr/μ>10).Table 4
Minimum value ofQ to obtain relative errors for Pb lower than 10−4 when Pi=0.
ρμr=0.001μr=0.01μr=0.1μr=1.00.4Wil9101010FM9974Fal11974NR7754AP10741HM17653HM23543Pb =5.633⋅10−3Pb= 6.442⋅10−3Pb =7.998⋅10−3Pb =8.696⋅10−30.6Wil24191817FM2316115Fal2516117NR181396AP221361HM1151074HM25974Pb =5.683⋅10−2Pb =6.954⋅10−2Pb =9.089⋅10−2Pb =9.979⋅10−20.8Wil74454039FM68341721Fal70352314NR53271710AP572491HM14121137HM23620137Pb=0.2469Pb=0.2982Pb=0.3750Pb=0.4043In Figure6 we show the computational cost
associated to all the presented methods when ρ=0.8 and μr=0.01 for different values of Q.
Results show that the lowest computation costs are for the simplest methods
(Wil, Fal, NR) followed by HM1 and HM2. Finally, and computationally speaking,
the worst methods are FM and AP.Figure 6
Computation time for the different methods deployed.Pi=0.


## 6. Conclusions

We proposed and compared two novel methods that
effectively obtain the value of typical performance parameters in retrial
systems with user impatience. As there are not closed-form solutions to these
systems when there are more than two servers, approximate methods are required
to solve such systems. The proposed methods are improvements of the method
proposed by Neuts and Rao in [11]
being also based on the homogenization of the state space beyond a certain
level. The results show the better performance of generalized truncated methods
compared to finite truncated methods in terms of accuracy. Comparing the generalized
truncated methods among them, we conclude that the proposed HM2 method
outperforms, in almost all cases, previous proposals in terms of accuracy for a
wide range of scenarios and for the performance parameters studied, with
moderated computation cost growths compared to the simpler methods.



---

*Source: 183089-2008-05-18.xml*