Asymptotic Diffusion Analysis of Multi-Server Retrial Queue with Hyper-Exponential Service

: A multi-server retrial queue with a hyper-exponential service time is considered in this paper. The study is performed by the method of asymptotic diffusion analysis under the condition of long delay in orbit. On the basis of the constructed diffusion process, we obtain approximations of stationary probability distributions of the number of customers in orbit and the number of busy servers. Using simulations and numerical analysis, we estimate the accuracy and applicability area of the obtained approximations.


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 Sections 3 and 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.

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 µ 1 , µ 2 > 0, 0 < q 1 < 1, and q 1 + q 2 = 1.
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 summing operator by E. Applying this operator on Equation (6) and taking into account that (A + B)E = 0 and we obtain the following scalar equation: It will be an additional equation to the system Equation (6).

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: σ → 0.
E{εi(τ)} be a scalar function which determines the average number of customers in the orbit normalized by ε = σ asymptotically for ε → 0. Denote probabilities of the number of busy servers working at the first phase (n 1 ) or second phase (n 2 ) of given hyper-exponential distribution of service time by R(n 1 , n 2 , x) for n 1 , n 2 = 0, . . . , N and n 1 + n 2 ≤ N. Let R(x) be (N + 1) × (N + 1) left-top triangle matrix with entries equal to R(n 1 , n 2 , x) for n 1 + n 2 ≤ N and others equal to zero.

Theorem 1.
Under asymptotic condition ε → 0, probabilities R(n 1 , n 2 , x) are determined by expressions where and x = x(τ) is a solution of the ordinary differential equation Proof. Making asymptotic transition ε → 0 in Equations (11) and (12) and denoting lim ε→0 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 R(n 1 , n 2 , x) which solution depends on x(τ). Equation (19) (that coincides with Equation (15)) is an ordinary differential equation of the first order for function x(τ). 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: (13) and (14), we obtain true identities. So, the theorem is proved.
then function a(x) may be written in the form Proof. From Equation (20) we derive Taking into account that ∑ n 1 +n 2 ≤N R(n 1 , n 2 , x) = 1 and using notation Equation (21), we obtain Equation (22).

Second Stage of the Asymptotic Analysis
Taking into account that characteristic function of process 1 σ x(σt) of the normalized number of customers in the orbit has form exp{jwx(τ)} (see Equation (17)), let us make the following substitution in system Equation (9): (here H (1) (u, t)E is characteristic function of centered process i(t) − 1 σ x(σt)). Then we obtain the system Taking into account notation Equation (20), we can rewrite this system in the form Denoting σ = ε 2 and making the following substitutions in Equation (24) we obtain system of differential equations Let us denote We can prove the following theorem.
Theorem 2. Function F (1) (w, τ) has the following form: where non-zero entries R(n 1 , n 2 , x) of matrix R(x) are determined by Theorem 1, and scalar function Φ(w, τ) is a solution of differential equation Here function a(x) is determined by Equation (20) and top-left triangle matrix g(x) 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 Φ(w, τ) and f(x) 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 jεwΦ(w, τ) and tending ε → 0, 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 g(x) 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 is a particular solution to non-homogeneous system Equation (37), so it should satisfy some additional condition which we choose in the form g(x)E = 0 (as was declared in Equation (30)). Then solution g(x) of system Equation (37) Then we derive Applying Equation (20), we obtain Let us divide this equation by ε 2 : Substituting Equation (35) here, we obtain 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.
In the next section, we will use obtained functions a(x) and b(x) as a drift and diffusion coefficients of a diffusion process which we apply to derive approximations for probability distributions under study.

Method of Asymptotic Diffusion Analysis
Taking into account notation σ = ε 2 and substitutions Equation (25) from the previous section, we consider stochastic process under asymptotic condition σ → 0. Here i(τ) is the number of customers in the orbit and function x(τ) is determined in the previous section as a solution of differential equation x (τ) = a(x). We can prove the following statement.

Lemma 1. Asymptotic stochastic process
is a solution of stochastic differential equation that depends on continuous parameter x. (22)

Proof. Consider Equation (28) with a(x) and b(x) determined by Equations
which is the Fokker-Planck equation for p.d.f. P(y, τ). Hence, stochastic process y(τ) is a diffusion process with drift coefficient a (x)y and diffusion coefficient b(x). Therefore, process y(τ) is a solution of stochastic differential Equation (44).
Let us consider the following stochastic process: (here ε = √ σ as before.) Lemma 2. Stochastic process z(τ) is a solution to stochastic differential equation with a precision up to an infinitesimal of order ε 2 .
Proof. Due to x(τ) is a solution of differential equation dx(τ) = a(x)dτ and process y(τ) 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 O(ε 2 ).
Suppose that the system is in steady-state regime. Consider stationary p.d.f. s(z) for process z(τ): Let us prove the following statement.
where C is a normalizing constant.
Proof. Due to z(τ) is a solution of stochastic differential Equation (46), the process is diffusion with drift coefficient a(z) and diffusion coefficient b(z). Therefore, its p.d.f. s(z) 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 s(∞) = 0, we obtain Equation (48). The theorem is proved.

Approximations of the Stationary Distributions
Consider matrix characteristic function F (1) (w, τ) (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 Φ(w, τ) of the number of customers in the orbit and matrix R(x) 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 P{n 1 (t) = n 1 , n 2 (t) = n 2 , i(t) = i} = P(n 1 , n 2 , i, t) in the form of production P(n 1 , n 2 , i, t) = R(n 1 , n 2 )P(i), where R(n 1 , n 2 ) 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), P(i) is the stationary distribution of the number of customers in the orbit.
To construct an approximationP(i) for discrete probability distribution P(i), consider p.d.f. s(z) of stochastic process z(τ). In the non-asymptotic case we can represent z(τ) as follows: There are various ways to construct probability distributionP(i) of discrete process σi(τ) on the base of p.d.f. s(z) of continuous stochastic process z(τ). We will use the following one.
Using Equations (48) and (49), we write non-negative function G(i) of discrete argument i ∈ {0, 1, . . . } in the form Taking into account normalization condition, we can writê This probability distribution we will use as an approximation for probability distribution P(i) = P {i(t) = i} 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 R(n) 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 a(x) = 0, which is an equation for variable x. Let us denote its positive solution by κ. Substituting x = κ into Equation (13), we obtain the following expression for calculation of the first-type approximation R 1 (n 1 , n 2 ) of probabilities R(n 1 , n 2 ):R 1 (n 1 , n 2 ) = R(n 1 , n 2 , κ).
Approximation Equation (51) allows us determine another form of approximationR 2 (n 1 , n 2 ) for probability distribution R(n 1 , n 2 ) of the number of servers working in the first and in the second phases: where R(n 1 , n 2 , x) is calculated using Equation (13). So, basing on these two types of approximations, we can write two typesR 1 (n) andR 2 (n) of approximations for the stationary probability distribution R(n) of the number of busy servers as follows:R The accuracy of the constructed approximations and their applicability area will be discussed in the next section.

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 P 1 (ν) and P 2 (ν) are probability distributions of two discrete variables. For our purposes, one of them will be an approximation fromP(i),R 1 (n) orR 2 (n), 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 0.05 as enough to an approximation be enough accurate. We reach ∆ ≤ 0.001 for simulation experiments themselves (when both P 1 (ν) and P 2 (ν) 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 (N = 5) and hyper-exponential service time with parameters Let us denote load of the system by ρ (0 < ρ < 1), and value of arrival process intensity λ we determine as λ = ρN.
In Table 1, you may find Kolmogorov distances ∆ k for approximationsR k 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 (∆ k is decreasing). We highlight values of Kolmogorov distance that satisfy chosen applicability condition ∆ ≤ 0.05 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 ρ ≤ 0.7 and the second one is more preferable for ρ > 0.7. Furthermore, we may notice that for σ < 1 both of them are applicable for all cases ρ ≥ 0.6, and the second approximation is applicable for greater values σ for narrower area (ρ ≥ 0.9).  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 σ < 2 while ρ ≥ 0.9 (boldface is used to highlight applicable cases again).

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 µ 2 = m. Some calls (about 20%) require long time to be serviced. We model them using rate µ 1 = 0.2m for exponential distribution. So, we can generally model service time in our system as hyper-exponential with parameters µ 1 = 0.2m, µ 2 = m, q 1 = 0.2, q 2 = 0.8. 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 m · C s , where C s is some constant. Furthermore, we have some losses due to customers that should make repeated calls. We calculate them asī · C E where C E is some constant andī = ∞ ∑ i=0 i · P i 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 C L · ∞ ∑ i=w+1 P i , where C L 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 C s = 3, C E = 0.04, C L = 2, 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 K(m) for the cases w = 10, w = 20 and w = 50 ( Figure 1). As we see, in all these cases function K(m) 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.

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 ∆ ≤ 0.05, 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.