1. Introduction
Nowadays, retrial queues (RQ) are very popular mathematical models of various real systems. These may be call centers [
1,
2] where a calling customer, if he or she find all operators busy, tries to make a new call again after some (random) time. In the case of modeling of data transmission via TCP [
3] if a data packet is lost or corrupted, the packet should be retransmitted after certain time. In cellular mobile networks, more complex models of retrial queues are used (see [
4,
5]). Furthermore, we see great opportunities in using multi-server retrial queues for modeling and design data processing systems with high-intensity data flows [
6,
7]. These are not only fields where retrial queues may be applied (e.g., survey [
8]).
A retrial queue is a such queueing system with limited number of servers where customers that arrived and do not find a free server are not lost, but go to special ‘waiting’ state, or are moved to special buffer (
orbit) for waiting. In the orbit, the customers are not organizing in a queue (e.g., first come, first served (FCFS)), but wait for a random time and then try to obtain the service again. There are a great number of RQs modifications considered by various authors—From ordinary single-server or multi-server systems to retrial queues with several orbits, priorities, collisions and so on. Detail information about retrial queue models, their analysis and applications you may find in surveys [
8,
9,
10,
11,
12] and others.
Most of the authors consider single-line retrial queues and or retrial queues with exponentially distributed service time. In this paper, we consider a multi-line retrial queue with hyper-exponential service time, but we hold in mind a such model with arbitrary distributed service time as the aim for future work.
Unfortunately, exact analytical solutions for retrial queues may be derived in very rare cases. Due to this, the most authors use numerical methods, approximations or even simulation methods to study more complex retrial queues [
13,
14,
15,
16,
17]. This also applies to multi-server models [
16] including systems with non-exponential service time [
13].
Traditionally, we study retrial queues under various asymptotic conditions such as heavy load or long delay [
18,
19,
20,
21]. In these studies we and our colleagues obtained approximations for stationary probability distributions of the number of customers in the orbit that are applicable in appropriate conditions. Unlike these previous works, in this paper, we apply a new approach—The asymptotic diffusion method [
22]—To perform more detailed and accurate analysis of the model. This will give more precise approximations and wider areas of their applicability.
The rest of the paper is organized as follows. In
Section 2, we describe in detail a mathematical model of the considered retrial queue and derive balance equations. In
Section 3 and
Section 4, we perform asymptotic analysis of obtained equations under condition of long delay in the orbit. In
Section 5, we introduce the method of asymptotic diffusion analysis and construct a diffusion process which is used in
Section 6 to build approximations for steady-state probability distributions of the number of busy servers and the number of customers in the orbit. Using simulations and numerical experiments, we determine the applicability area of suggested approximations in
Section 7.
2. Mathematical Model
Consider a retrial queue with
N servers and Poisson arrivals with intensity
. Service times are independent and identically distributed as hyper-exponential variables with cumulative distribution function
where
,
, and
.
When a customer arrives, if there is a free server the customer moves to it for service, otherwise the customer goes to the orbit and stays there during random period exponentially distributed with parameter . After this period the customer tries to get a service again according to the described algorithm.
Let us denote the number of busy servers at the moment
t by
, and the number of customers in the orbit by
. The problem is to find steady-state probability distribution of the two-dimensional process
that we denote by
for all
and
independently of
t.
To solve the problem, let us consider the stochastic process
, where
and
are the numbers of busy servers that realizes the service time of the first phase or the second phase respectively. The process is Markovian, hence we can write balance equations for its probability distribution
where
,
. So, we can write differential equations
for the case
, and
for the case
.
Let us introduce partial characteristic functions
where
. Then we can rewrite the system of balance Equations (
2)–(
3) as follows:
for
, and
for
.
Let us use linear finite difference operators
to rewrite the system Equations (
4) and (
5) in a compact form. Then using notation
for
matrix function with indexes
and entries equal to
for
or zero otherwise, we can rewrite Equations (
4) and (
5) as follows:
Here
where indicator
is defined as follows:
We write these operators on the right hand of matrix function affected by them similar to notations of matrix operations.
Let us denote summing operator by
. Applying this operator on Equation (
6) and taking into account that
and
we obtain the following scalar equation:
It will be an additional equation to the system Equation (
6).
3. First Stage of the Asymptotic Analysis
System of Equations (
6) and (
9) can not be solved in a direct way, so, we will solve it under asymptotic condition of long delay in the orbit:
.
Denoting
, we make the following substitutions in system Equations (
6) and (
9):
We will solve system Equations (
11) and (
12) under asymptotic condition
.
Let be a scalar function which determines the average number of customers in the orbit normalized by asymptotically for . Denote probabilities of the number of busy servers working at the first phase () or second phase () of given hyper-exponential distribution of service time by for and . Let be left-top triangle matrix with entries equal to for and others equal to zero.
Theorem 1. Under asymptotic condition , probabilities are determined by expressions and is a solution of the ordinary differential equation Proof. Making asymptotic transition
in Equations (
11) and (
12) and denoting
, we derive
Let us find a solution of system Equation (
16) in the form
Substituting Equation (
17) into Equation (
16), we derive
So, Equation (
18) is a homogeneous system of linear equations for
which solution depends on
. Equation (
19) (that coincides with Equation (
15)) is an ordinary differential equation of the first order for function
. So, we can find a solution of Equation (
18) and substitute it into Equation (
19).
Let us solve system Equation (
18). To do this, we rewrite matrix Equation (
18) in the form of the following homogeneous system of linear equations:
for
, and
for
. Substituting here expressions Equations (
13) and (
14), we obtain true identities. So, the theorem is proved. □
Corollary 1. then function may be written in the form Proof. From Equation (
20) we derive
Taking into account that
and using notation Equation (
21), we obtain Equation (
22). □
4. Second Stage of the Asymptotic Analysis
Taking into account that characteristic function of process
of the normalized number of customers in the orbit has form
(see Equation (
17)), let us make the following substitution in system Equation (
9):
(here
is characteristic function of centered process
). Then we obtain the system
Taking into account notation Equation (
20), we can rewrite this system in the form
Denoting
and making the following substitutions in Equation (
24)
we obtain system of differential equations
We can prove the following theorem.
Theorem 2. Function has the following form:where non-zero entries of matrix are determined by Theorem 1, and scalar function is a solution of differential equation Here function is determined by Equation (20), has the formand top-left triangle matrix is determined by the following system of equations: Proof. Let us write only the members of the first equation of system Equation (
26) that contains only zero or the first order of
:
Solution of this equation we write in the form
where
and
are some scalar and matrix functions whose expressions we will obtain later. Substituting this expression into Equation (
31), we derive
Taking into account Equation (
20), dividing Equation (
33) by
and tending
, we obtain
Let us rewrite this equality in the form of non-homogeneous equation
According to superposition principle, we can write a solution of this equation in the form
Substituting it into Equation (
34), we obtain equations
Notice that Equation (
37) for matrix
coincides with the first expression of Equation (
30), therefore statement Equation (
30) is true (the second expression of Equation (
30) will be discussed later).
Now, consider Equation (
18). Let us differentiate it by
x, we derive equation
Taking into account Equation (
36), we can conclude that
where
.
Matrix
is a particular solution to non-homogeneous system Equation (
37), so it should satisfy some additional condition which we choose in the form
(as was declared in Equation (
30)). Then solution
of system Equation (
37) is uniquely determined.
Now, let us consider the second Equation of Equation (
26). We substitute expansion Equation (
32) in it and write the members of derived equation including up to the second order of
:
Applying Equation (
20), we obtain
Let us divide this equation by
:
Substituting Equation (
35) here, we obtain
Denoting
we can rewrite Equation (
39) in the form
Differentiating Equation (
20) by
x, we obtain
Then, taking into account Equation (
38) and substituting Equation (
42) into Equation (
41), we derive the equation
that coincides with Equation (
28). So, the theorem is proved. □
Corollary 2. Let , , then function from Equation (40) may be written in the form Proof. Let us rewrite Equation (
40) in scalar form:
Applying here condition
, we obtain Equation (
43). □
In the next section, we will use obtained functions and as a drift and diffusion coefficients of a diffusion process which we apply to derive approximations for probability distributions under study.
5. Method of Asymptotic Diffusion Analysis
Taking into account notation
and substitutions Equation (
25) from the previous section, we consider stochastic process
under asymptotic condition
. Here
is the number of customers in the orbit and function
is determined in the previous section as a solution of differential equation
. We can prove the following statement.
Lemma 1. Asymptotic stochastic processis a solution of stochastic differential equationthat depends on continuous parameter x. Proof. Consider Equation (
28) with
and
determined by Equations (
22) and (
40). Solution of this equation determines asymptotic characteristic function for centered and normalized stochastic process
of the number
of customers in the orbit under condition
.
Let us make in this equation an inverse Fourier transform on argument
w. Then for probability density function (p.d.f.)
of process
, we obtain the equation
which is the Fokker–Planck equation for p.d.f.
. Hence, stochastic process
is a diffusion process with drift coefficient
and diffusion coefficient
. Therefore, process
is a solution of stochastic differential Equation (
44). □
Let us consider the following stochastic process:
(here as before.)
Lemma 2. Stochastic process is a solution to stochastic differential equationwith a precision up to an infinitesimal of order . Proof. Due to
is a solution of differential equation
and process
satisfies Equation (
44), the following equality is true:
We can represent its coefficients in the form
and, so, we can rewrite Equation (
47) as follows:
It coincides with Equation (
46) with a precision up to infinitesimal
. □
Suppose that the system is in steady-state regime. Consider stationary p.d.f.
for process
:
Let us prove the following statement.
Theorem 3. Stationary p.d.f. of stochastic process has the formwhere C is a normalizing constant. Proof. Due to
is a solution of stochastic differential Equation (
46), the process is diffusion with drift coefficient
and diffusion coefficient
. Therefore, its p.d.f.
is a solution of the Fokker–Planck equation
This equation is an ordinary differential equation of the second order. Solving it and taking into account normalization condition and boundary condition
, we obtain Equation (
48). The theorem is proved. □
6. Approximations of the Stationary Distributions
Consider matrix characteristic function
(Theorem 2) of the joint distribution of the number of customers in the orbit and the number of servers working in the first and at the second phases. Due to Equation (
27), it equals to production of characteristic function
of the number of customers in the orbit and matrix
that represents two-dimensional probability distribution of the number of servers working in the first and at the second phases. Therefore, we can represent stationary probability distribution under study
in the form of production
where
is the joint two-dimensional stationary distribution of the number of servers working in the first and in the second phases of hyper-exponential service time Equation (
1),
is the stationary distribution of the number of customers in the orbit.
To construct an approximation
for discrete probability distribution
, consider p.d.f.
of stochastic process
. In the non-asymptotic case we can represent
as follows:
There are various ways to construct probability distribution of discrete process on the base of p.d.f. of continuous stochastic process . We will use the following one.
Using Equations (
48) and (
49), we write non-negative function
of discrete argument
in the form
Taking into account normalization condition, we can write
This probability distribution we will use as an approximation for probability distribution of the number of customers in the orbit of considered retrial queue in steady-state regime.
We can construct two different approximations for stationary probabilities
of the number of busy servers. The first one can be build on the base of Formulas (
13) and (
14). In stationary regime, differential Equation (
15) transforms to the equality
which is an equation for variable
x. Let us denote its positive solution by
. Substituting
into Equation (
13), we obtain the following expression for calculation of the first-type approximation
of probabilities
:
Approximation Equation (
51) allows us determine another form of approximation
for probability distribution
of the number of servers working in the first and in the second phases:
where
is calculated using Equation (
13).
So, basing on these two types of approximations, we can write two types
and
of approximations for the stationary probability distribution
of the number of busy servers as follows:
The accuracy of the constructed approximations and their applicability area will be discussed in the next section.
7. Applicability Area of Obtained Approximations
To estimate the accuracy of obtained approximations Equations (
51)–(
54) we perform numerical experiments with usage of simulation modeling. We use the Kolmogorov distance
as an error estimation. Here
and
are probability distributions of two discrete variables. For our purposes, one of them will be an approximation from
,
or
, and another will be an empiric distribution of the same parameter constructed on the base of simulation results. Simulation model and simulation parameters (e.g., volume of collected statistical data) are chosen in a such way to error of simulation be not have a significant influence on calculation of value Equation (
55) for approximations under test.
We consider values of Kolmogorov distance Equation (
55) less or equal to
as enough to an approximation be enough accurate. We reach
for simulation experiments themselves (when both
and
in Equation (
55) are taken from simulation results of two experiments with identical parameters).
Special developed software based on
ODIS system [
23,
24] is used to perform simulations. Numerical calculations of probabilities Equations (
51)–(
54) are made using
PTC Mathcad application.
For our purposes of analyzing applicability area of results Equations (
51)–(
54), let us consider a retrial queue with five servers (
) and hyper-exponential service time with parameters
Let us denote load of the system by (), and value of arrival process intensity we determine as .
In
Table 1, you may find Kolmogorov distances
for approximations
from Equation (
54) calculated for various values of load parameter
and intensity of retrials
. As we see, while
is decreasing, the accuracy of the approximations becomes better (
is decreasing). We highlight values of Kolmogorov distance that satisfy chosen applicability condition
by boldface font. We can not conclude what approximation from these is better in all cases but the first one seems more preferable in the case of
and the second one is more preferable for
. Furthermore, we may notice that for
both of them are applicable for all cases
, and the second approximation is applicable for greater values
for narrower area (
).
The similar results for approximation of the number of customers in the orbit Equation (
51) are presented in
Table 2. As in the previous case, this approximation became more accurate while
is decreasing too. In addition, this approximation became more precise almost always while
is growing. We can assume that it is applicable for all
while
(boldface is used to highlight applicable cases again).
8. Call Center Example
Consider an example which illustrates application of obtained results in real system. Some call center has five workstations for operators that service calls from customers. Ordinary call requires time to be serviced equal to 1/m in average, where m is some measure of the average productivity of one operator and its workstation. For such calls we use exponentially distributed service time with rate = m. Some calls (about 20%) require long time to be serviced. We model them using rate m for exponential distribution. So, we can generally model service time in our system as hyper-exponential with parameters m, = m, . Let arrival process be Poisson with intensity = 2. If a customer find the system busy (there is no free operator), he or she tries to make a new call after a random time that we suppose exponentially distributed with parameter = 0.2. So, we can consider a model of the call center in the form of multi-server retrial queue with hyper-exponential service time and mentioned parameters.
Let us consider that the call center operation requires some cost and has some penalties. We suppose that we can improve a performance of the call center by increasing an average speed of calls processing
m (e.g., involving more qualified personal or deploying better equipment) but we should pay a cost equals to
, where
is some constant. Furthermore, we have some losses due to customers that should make repeated calls. We calculate them as
where
is some constant and
is mathematical expectation of the number of calls in the orbit. In addition, if we have too many "waiting" customers we obtain a penalty which we calculate as
, where
is some constant and
w is some upper bound after which calls in the orbit are taken into account for the penalty calculation. Thus we can write the total cost of the call center operation
K in the following form:
Using the following values of constants
and varying value of
m in the interval from 0.75 to 1.2 (with step 0.05) that corresponds to values of the system load
from 0.96 to 0.6 respectively, we obtain values of the total cost
for the cases
,
and
(
Figure 1). As we see, in all these cases function
has a minimum in the considered range of
m. Therefore, when we plan work of the call center we should take into account that increasing of performance of operators and workstations does not always cause the cost be less. Moreover, value of
m for the minimum point depends on value of upper bound
w.
9. Discussion
In this paper, we apply a new approach [
22] for studying multi-server retrial queue with hyper-exponential service time. Analysis of the model is made under asymptotic condition of long delay in the orbit. Using asymptotic analysis technique [
18,
19], we obtain parameters of diffusion process Equation (
45) which is used to build approximations for stationary probability distributions of the number of customers in the orbit Equation (
51) and the number of busy servers Equations (
52)–(
54).
We use numerical experiments and simulations for estimation of the approximations accuracy. Using Kolmogorov distance
Equation (
55) as an error of approximation and criteria of approximation applicability in the form
, we obtain applicability areas of proposed approximations. As it shown, these areas are wide enough. At least in the most of cases, they are wider than applicability areas of approximations obtained using basic asymptotic analysis technique as it was made in [
18,
19,
20,
21].
We suppose to use asymptotic diffusion method and obtained results in future works that will be devoted as to analysis of retrial queues with other configurations both to study of multi-server retrial queue with arbitrary distributed service time. Using hyper-exponential distribution as an approximation for distributions of other types may be the first step in this direction. Further studies may involve supplementary variables technique or using Markov chains based on departure epochs.
We thank the reviewers of this paper whose comments allow us improve the work.