Next Article in Journal
Turbulence as a Network of Fourier Modes
Next Article in Special Issue
On Probability Characteristics for a Class of Queueing Models with Impatient Customers
Previous Article in Journal
A Mathematical Pre-Disaster Model with Uncertainty and Multiple Criteria for Facility Location and Network Fortification
Previous Article in Special Issue
Renewal Redundant Systems Under the Marshall–Olkin Failure Model. A Probability Analysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Institute of Applied Mathematics and Computer Science, Tomsk State University, Tomsk 634050, Russia
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(4), 531; https://doi.org/10.3390/math8040531
Submission received: 15 March 2020 / Revised: 31 March 2020 / Accepted: 1 April 2020 / Published: 3 April 2020

Abstract

:
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.

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
B ( x ) = q 1 1 e μ 1 x + q 2 1 e μ 2 x ,
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 the number of busy servers at the moment t by n ( t ) { 0 , 1 , , N } , and the number of customers in the orbit by i ( t ) { 0 , 1 , 2 } . The problem is to find steady-state probability distribution of the two-dimensional process { n ( t ) , i ( t ) } that we denote by
P ( n , i ) = P { n ( t ) = n , i ( t ) = i }
for all n { 0 , 1 , , N } and i { 0 , 1 , 2 } independently of t.
To solve the problem, let us consider the stochastic process { n 1 ( t ) , n 2 ( t ) , i ( t ) } , where n 1 ( t ) and n 2 ( t ) 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
P ( n 1 , n 2 , i , t ) = P { n 1 ( t ) = n 1 , n 2 ( t ) = n 2 , i ( t ) = i } ,
where 0 n 1 + n 2 N , i { 0 , 1 , 2 } . So, we can write differential equations
P ( n 1 , n 2 , i , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 + i σ ) P ( n 1 , n 2 , i , t ) + q 1 λ P ( n 1 1 , n 2 , i , t ) + q 2 λ P ( n 1 , n 2 1 , i , t ) + q 1 ( i + 1 ) σ P ( n 1 1 , n 2 , i + 1 , t ) + q 2 ( i + 1 ) σ P ( n 1 , n 2 1 , i + 1 , t ) + ( n 1 + 1 ) μ 1 P ( ( n 1 + 1 , n 2 , i , t ) ) + ( n 2 + 1 ) μ 2 P ( n 1 , n 2 + 1 , i , t )
for the case n 1 + n 2 N 1 , and
P ( n 1 , n 2 , i , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 ) P ( n 1 , n 2 , i , t ) + q 1 λ P ( n 1 1 , n 2 , i , t ) + q 2 λ P ( n 1 , n 2 1 , i , t ) + λ P ( n 1 , n 2 , i 1 , t ) + q 1 ( i + 1 ) σ P ( n 1 1 , n 2 , i + 1 , t ) + q 2 ( i + 1 ) σ P ( n 1 , n 2 1 , i + 1 , t )
for the case n 1 + n 2 = N .
Let us introduce partial characteristic functions
H ( n 1 , n 2 , u , t ) = i = 0 e j u i P ( n 1 , n 2 , i , t ) ,
where j = 1 . Then we can rewrite the system of balance Equations (2)–(3) as follows:
H ( n 1 , n 2 , u , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + q 1 λ H ( n 1 1 , n 2 , u , t ) + q 2 λ H ( n 1 , n 2 1 , u , t ) + ( n 1 + 1 ) μ 1 H ( n 1 + 1 , n 2 , u , t ) + ( n 2 + 1 ) μ 2 H ( n 1 , n 2 + 1 , u , t ) j σ q 1 e j u H ( n 1 1 , n 2 , u , t ) u j σ q 2 e j u H ( n 1 , n 2 1 , u , t ) u + j σ H ( n 1 , n 2 , u , t ) u
for n 1 + n 2 N 1 , and
H ( n 1 , n 2 , u , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + q 1 λ H ( n 1 1 , n 2 , u , t ) + q 2 λ H ( n 1 , n 2 1 , u , t ) + λ e j u H ( n 1 , n 2 , u , t ) j σ q 1 e j u H ( n 1 1 , n 2 , u , t ) u j σ q 2 e j u H ( n 1 , n 2 1 , u , t ) u
for n 1 + n 2 = N .
Let us use linear finite difference operators A , B , I 0 , I 1 to rewrite the system Equations (4) and (5) in a compact form. Then using notation H ( u , t ) for ( N + 1 ) × ( N + 1 ) matrix function with indexes n 1 , n 2 = 0 , , N and entries equal to H ( n 1 , n 2 , u , t ) for n 1 + n 2 N or zero otherwise, we can rewrite Equations (4) and (5) as follows:
H ( u , t ) t = H ( u , t ) A + e j u B + j σ H ( u , t ) u I 0 e j u I 1 .
Here
H ( u , t ) A = ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + q 1 λ H ( n 1 1 , n 2 , u , t ) + q 2 λ H ( n 1 , n 2 1 , u , t ) + ψ ( n 1 + n 2 , N ) ( n 1 + 1 ) μ 1 H ( n 1 + 1 , n 2 , u , t ) + ( n 2 + 1 ) μ 2 H ( n 1 , n 2 + 1 , u , t ) ,
H ( u , t ) B = 1 ψ ( n 1 + n 2 , N ) λ H ( n 1 , n 2 , u , t ) ,
H ( u , t ) I 0 = ψ ( n 1 + n 2 , N ) H ( n 1 , n 2 , u , t ) ,
H ( u , t ) I 1 = q 1 H ( n 1 1 , n 2 , u , t ) + q 2 H ( n 1 , n 2 1 , u , t ) ,
where indicator ψ ( n 1 + n 2 , N ) is defined as follows:
ψ ( n 1 + n 2 , N ) = 1 , if   0 n 1 + n 2 < N , 0 , otherwise : n 1 + n 2 = N .
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 E . Applying this operator on Equation (6) and taking into account that
( A + B ) E = 0
and
( I 0 I 1 ) E = 0 ,
we obtain the following scalar equation:
H ( u , t ) t E = e j u 1 H ( u , t ) B + j σ e j u H ( u , t ) u I 0 E .
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: σ 0 .
Denoting σ = ε , we make the following substitutions in system Equations (6) and (9):
τ = t ε , u = ε w , H ( u , t ) = F ( w , τ , ε ) .
We obtain
ε F ( w , τ , ε ) τ = F ( w , τ , ε ) A + e j ε w B + j F ( w , τ , ε ) w I 0 e j ε w I 1 ,
ε F ( w , τ , ε ) τ E = e j ε w 1 F ( w , τ , ε ) B E + j e j ε w F ( w , τ , ε ) w I 0 E .
We will solve system Equations (11) and (12) under asymptotic condition ε 0 .
Let x = x ( τ ) = lim ε 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
R ( n 1 , n 2 , x ) = ρ 1 ( x ) n 1 n 1 ! · ρ 2 ( x ) n 2 n 2 ! R ( 0 , 0 , x ) , R ( 0 , 0 , x ) = 1 / n 1 + n 2 N ρ 1 ( x ) n 1 n 1 ! · ρ 2 ( x ) n 2 n 2 ! ,
where
ρ 1 ( x ) = q 1 ( λ + x ) μ 1 , ρ 2 ( x ) = q 2 ( λ + x ) μ 2 ,
and x = x ( τ ) is a solution of the ordinary differential equation
x ( τ ) = R ( x ) B x ( τ ) I 0 E .
Proof. 
Making asymptotic transition ε 0 in Equations (11) and (12) and denoting lim ε 0 F ( w , τ , ε ) = F ( w , τ ) , we derive
F ( w , τ ) ( A + B ) + j F ( w , τ ) w ( I 0 I 1 ) = 0 , F ( w , τ ) τ E = j w F ( w , τ ) B E + j F ( w , τ ) w I 0 E .
Let us find a solution of system Equation (16) in the form
F ( w , τ ) = R ( x ) e j w x ( τ ) .
Substituting Equation (17) into Equation (16), we derive
R ( x ) A + B x ( τ ) I 0 I 1 = 0 ,
x ( τ ) = R ( x ) B x ( τ ) I 0 E .
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:
q 1 ( λ + x ) R ( n 1 1 , n 2 , x ) + q 2 ( λ + x ) R ( n 1 , n 2 1 , x ) ( λ + x + n 1 μ 1 + n 2 μ 2 ) R ( n 1 , n 2 , x ) + ( n 1 + 1 ) μ 1 R ( n 1 + 1 , n 2 , x ) + ( n 2 + 1 ) μ 2 R ( n 1 , n 2 + 1 , x ) = 0
for n 1 + n 2 N 1 , and
q 1 ( λ + x ) R ( n 1 1 , n 2 , x ) + q 2 ( λ + x ) R ( n 1 , n 2 1 , x ) ( n 1 μ 1 + n 2 μ 2 ) R ( n 1 , n 2 , x ) = 0
for n 1 + n 2 = N . Substituting here expressions Equations (13) and (14), we obtain true identities. So, the theorem is proved. □
Let us denote
a ( x ) = R ( x ) B x ( τ ) I 0 E .
Corollary 1.
Let
r ( x ) = n 1 + n 2 = N R ( n 1 , n 2 , x ) ,
then function a ( x ) may be written in the form
a ( x ) = ( λ + x ( τ ) ) r ( x ) x ( τ )
Proof. 
From Equation (20) we derive
a ( x ) = R ( x ) B x ( τ ) I 0 E = a ( R ( x ) B E x ( τ ) R ( x ) I 0 E = λ n 1 + n 2 = N R ( n 1 , n 2 , x ) x n 1 + n 2 N 1 R ( n 1 , n 2 , x ) = λ n 1 + n 2 = N R ( n 1 , n 2 , x ) x n 1 + n 2 N R ( n 1 , n 2 , x ) n 1 + n 2 = N R ( n 1 , n 2 , x ) = ( λ + x ) n 1 + n 2 = N R ( n 1 , n 2 , x ) x n 1 + n 2 N R ( n 1 , n 2 , x ) .
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). □

4. 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 { j w x ( τ ) } (see Equation (17)), let us make the following substitution in system Equation (9):
H ( u , t ) = e j u σ x ( σ t ) H ( 1 ) ( u , t ) .
(here H ( 1 ) ( u , t ) E is characteristic function of centered process i ( t ) 1 σ x ( σ t ) ). Then we obtain the system
H ( 1 ) ( u , t ) t + j u x ( σ t ) H ( 1 ) ( u , t ) = H ( 1 ) ( u , t ) A + e j u B + j σ j σ x ( σ t ) H ( 1 ) ( u , t ) + H ( 1 ) ( u , t ) u I 0 e j u I 1 , H ( 1 ) ( u , t ) u + j u x ( σ t ) H ( 1 ) ( u , t ) E = e j u 1 H ( 1 ) ( u , t ) BE + j σ e j u j σ x ( σ t ) H ( 1 ) ( u , t ) + H ( 1 ) ( u , t ) u I 0 E .
Taking into account notation Equation (20), we can rewrite this system in the form
H ( 1 ) ( u , t ) t + j u a ( x ) H ( 1 ) ( u , t ) = H ( 1 ) ( u , t ) A + e j u B x I 0 e j u I 1 + j σ H ( 1 ) ( u , t ) u I 0 e j u I 1 , H ( 1 ) ( u , t ) u E + j u a ( x ) H ( 1 ) ( u , t ) E = e j u 1 H ( 1 ) ( u , t ) B e j u x I 0 + e j u j σ H ( 1 ) ( u , t ) u I 0 E .
Denoting σ = ε 2 and making the following substitutions in Equation (24)
τ = t ε 2 , u = ε w , H ( 1 ) ( u , t ) = F ( 1 ) ( w , τ , ε ) ,
we obtain system of differential equations
ε 2 F ( 1 ) ( w , τ , ε ) τ + j ε w a ( x ) F ( 1 ) ( w , τ , ε ) = F ( 1 ) ( w , τ , ε ) A + e j ε w B x I 0 e j ε w I 1 + j ε F ( 1 ) ( w , τ , ε ) w I 0 e j ε w I 1 , ε 2 F ( 1 ) ( w , τ , ε ) τ E + j ε w a ( x ) F ( 1 ) ( w , τ , ε ) E = e j ε w 1 F ( 1 ) ( w , τ , ε ) B e j ε w x I 0 + e j ε w j ε F ( 1 ) ( w , τ , ε ) w I 0 E .
Let us denote
lim ε 0 F ( 1 ) ( w , τ , ε ) = F ( 1 ) ( w , τ ) , lim ε 0 F ( 1 ) ( w , τ , ε ) τ = F ( 1 ) ( w , τ ) τ , lim ε 0 F ( 1 ) ( w , τ , ε ) w = F ( 1 ) ( w , τ ) w .
We can prove the following theorem.
Theorem 2.
Function F ( 1 ) ( w , τ ) has the following form:
F ( 1 ) ( w , τ ) = Φ ( w , τ ) R ( x ) ,
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
Φ ( w , τ ) τ = a ( x ) w Φ ( w , τ ) w + b ( x ) ( j w ) 2 2 Φ ( w , τ ) .
Here function a ( x ) is determined by Equation (20), b ( x ) has the form
b ( x ) = a ( x ) + 2 g ( x ) B x I 0 E + 2 x R ( x ) I 0 E ,
and top-left triangle matrix g ( x ) is determined by the following system of equations:
g ( x ) A + B + x I 1 I 0 = a ( x ) R ( x ) + R ( x ) x I 1 B , g ( x ) E = 0 .
Proof. 
Let us write only the members of the first equation of system Equation (26) that contains only zero or the first order of ε :
j ε w a ( x ) F ( 1 ) ( w , τ , ε ) = F ( 1 ) ( w , τ , ε ) A + B + j ε w B x I 0 I 1 j ε w I 1 + j ε F ( 1 ) ( w , τ , ε ) w I 0 I 1 + O ε 2 .
Solution of this equation we write in the form
F ( 1 ) ( w , τ , ε ) = Φ ( w , τ ) R ( x ) + j ε w f ( x ) + O ε 2 ,
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
j ε w a ( x ) Φ ( w , τ ) R ( x ) = Φ ( w , τ ) { R ( x ) A + B x I 0 I 1 + j ε w f ( x ) A + B x ( I 0 I 1 ) } + j ε Φ ( w , τ ) w R ( x ) I 0 I 1 + O ε 2 .
Taking into account Equation (20), dividing Equation (33) by j ε w Φ ( w , τ ) and tending ε 0 , we obtain
a ( x ) R ( x ) = f ( x ) A + B x I 0 I 1 + R ( x ) B x I 1 + Φ ( w , τ ) / w w Φ ( w , τ ) R ( x ) I 0 I 1 .
Let us rewrite this equality in the form of non-homogeneous equation
f ( x ) A + B + x I 1 I 0 = a ( x ) R ( x ) R ( x ) B x I 1 Φ ( w , τ ) / w w Φ ( w , τ ) R ( x ) I 0 I 1 .
According to superposition principle, we can write a solution of this equation in the form
f ( x ) = C R ( x ) + g ( x ) φ ( x ) Φ ( w , τ ) / w w Φ ( w , τ ) .
Substituting it into Equation (34), we obtain equations
φ ( x ) A + B + x I 1 I 0 = R ( x ) I 0 I 1 ,
g ( x ) A + B + x I 1 I 0 = a ( x ) R ( x ) + R ( x ) x I 1 B .
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
d R ( x ) d x A + B x I 0 I 1 R ( x ) I 0 I 1 = 0 .
Taking into account Equation (36), we can conclude that
φ ( x ) = d R ( x ) d x ,
where φ ( x ) E = d R ( x ) d x E = d ( R ( x ) E ) d x = 0 .
Matrix g ( x ) 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) 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 ε :
ε 2 Φ ( w , τ ) τ + j ε w a ( x ) Φ ( w , τ ) 1 + j ε w f ( x ) E = j ε w + ( j ε w ) 2 2 Φ ( w , τ ) R ( x ) + j ε w f ( x ) B x I 0 + j ε w x I 0 + j ε Φ ( w , τ ) w R ( x ) I 0 E + O ε 3 .
Then we derive
ε 2 Φ ( w , τ ) τ + j ε w a ( x ) Φ ( w , τ ) + ( j ε w ) 2 a ( x ) Φ ( w , τ ) f ( x ) E = Φ ( w , τ ) { j ε w R ( x ) B x I 0 + ( j ε w ) 2 f ( x ) B x I 0 + ( j ε w ) 2 R ( x ) x I 0 + ( j ε w ) 2 2 R ( x ) B x I 0 + ( j ε w ) 2 Φ ( w , τ ) / w w R ( x ) I 0 } E + O ε 3 .
Applying Equation (20), we obtain
ε 2 Φ ( w , τ ) τ + ( j ε w ) 2 a ( x ) Φ ( w , τ ) f ( x ) E = Φ ( w , τ ) { ( j ε w ) 2 f ( x ) B x I 0 + ( j ε w ) 2 R ( x ) x I 0 + ( j ε w ) 2 2 a ( x ) + ( j ε w ) 2 Φ ( w , τ ) / w w R ( x ) I 0 } E + O ε 3 .
Let us divide this equation by ε 2 :
Φ ( w , τ ) τ = Φ ( w , τ ) { ( j w ) 2 f ( x ) B x I 0 + ( j w ) 2 R ( x ) x I 0 + ( j w ) 2 2 a ( x ) ( j w ) 2 a ( x ) f ( x ) + ( j w ) 2 Φ ( w , τ ) / w w R ( x ) I 0 } E .
So, we derive
Φ ( w , τ ) / τ Φ ( w , τ ) = ( j w ) 2 2 2 f ( x ) B x I 0 + R ( x ) x I 0 a ( x ) f ( x ) E + a ( x ) w Φ ( w , τ ) / w Φ ( w , τ ) R ( x ) I 0 E .
Substituting Equation (35) here, we obtain
Φ ( w , τ ) / τ Φ ( w , τ ) = ( j w ) 2 2 2 g ( x ) B x I 0 E + 2 R ( x ) x I 0 E + a ( x ) + w Φ ( w , τ ) / w Φ ( w , τ ) φ ( x ) B x I 0 E R ( x ) I 0 E .
Denoting
b ( x ) = a ( x ) + 2 g ( x ) B x I 0 E + 2 x R ( x ) I 0 E ,
we can rewrite Equation (39) in the form
Φ ( w , τ ) τ = w Φ ( w , τ ) w φ ( x ) B x I 0 E R ( x ) I 0 E + ( j w ) 2 2 b ( x ) Φ ( w , τ ) .
Differentiating Equation (20) by x, we obtain
a ( x ) = R ( x ) x B x I 0 E R ( x ) I 0 E .
Then, taking into account Equation (38) and substituting Equation (42) into Equation (41), we derive the equation
Φ ( w , τ ) τ = a ( x ) w Φ ( w , τ ) w + b ( x ) ( j w ) 2 2 Φ ( w , τ )
that coincides with Equation (28). So, the theorem is proved. □
Corollary 2.
Let g 1 ( x ) = n 1 + n 2 = N g ( n 1 , n 2 , x ) , r ( x ) = n 1 + n 2 = N R ( n 1 , n 2 , x ) , then function b ( x ) from Equation (40) may be written in the form
b ( x ) = a ( x ) + 2 x 1 r ( x ) + ( λ + x ) g 1 ( x ) .
Proof. 
Let us rewrite Equation (40) in scalar form:
b ( x ) = a ( x ) + 2 λ n 1 + n 2 = N g ( n 1 , n 2 , x ) + x n 1 + n 2 N 1 g ( n 1 , n 2 , x ) x n 1 + n 2 N 1 R ( n 1 , n 2 , x ) = a ( x ) + 2 λ g 1 ( x ) + x 1 x n 1 + n 2 = N R ( n 1 , n 2 , x ) x n 1 + n 2 N g ( n 1 , n 2 , x ) g 1 ( x ) .
Applying here condition g ( x ) E = 0 , we obtain Equation (43). □
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.

5. Method of Asymptotic Diffusion Analysis

Taking into account notation σ = ε 2 and substitutions Equation (25) from the previous section, we consider stochastic process
σ i ( σ t ) 1 σ x ( σ t ) = σ i ( τ ) 1 σ x ( τ )
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
y ( τ ) = lim σ 0 σ i ( τ ) 1 σ x ( τ )
is a solution of stochastic differential equation
d y ( τ ) = a ( x ) y d τ + b ( x ) d w ( τ )
that depends on continuous parameter x.
Proof. 
Consider Equation (28) with a ( x ) and b ( x ) determined by Equations (22) and (40). Solution of this equation determines asymptotic characteristic function for centered and normalized stochastic process σ i ( τ ) 1 σ x ( τ ) of the number i ( t ) of customers in the orbit under condition σ 0 .
Let us make in this equation an inverse Fourier transform on argument w. Then for probability density function (p.d.f.) P ( y , τ ) = P y ( τ ) < y y of process y ( τ ) , we obtain the equation
P ( y , τ ) τ = y a ( x ) y P ( y , τ ) + 1 2 2 y 2 { b ( x ) P ( y , τ ) } ,
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:
z ( τ ) = x ( τ ) + ε y ( τ )
(here ε = σ as before.)
Lemma 2.
Stochastic process z ( τ ) is a solution to stochastic differential equation
d z ( τ ) = a ( z ) d τ + σ b ( z ) d w ( τ )
with a precision up to an infinitesimal of order ε 2 .
Proof. 
Due to x ( τ ) is a solution of differential equation d x ( τ ) = a ( x ) d τ and process y ( τ ) satisfies Equation (44), the following equality is true:
d z ( τ ) = d x ( τ ) + ε y ( τ ) = a ( x ) + ε y a ( x ) d τ + ε b ( x ) d w ( τ ) .
We can represent its coefficients in the form
a ( x ) + ε y a ( x ) + O ( ε 2 ) = a ( x + ε y ) = a ( z ) , ε b ( x ) + O ( ε 2 ) = ε b ( x ) + O ( ε ) = ε b ( x + ε y ) = ε b ( z ) ,
and, so, we can rewrite Equation (47) as follows:
d z ( τ ) = a ( z ) d τ + σ b ( z ) d w ( τ ) + O ( ε 2 ) .
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 ( τ ) :
s ( z , τ ) = s ( z ) = P z ( τ ) < z z .
Let us prove the following statement.
Theorem 3.
Stationary p.d.f. s ( z ) of stochastic process z ( τ ) has the form
s ( z ) = C b ( z ) exp 2 σ 0 z a ( x ) b ( x ) d x ,
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
z a ( z ) s ( z ) + 1 2 2 z 2 σ b ( z ) s ( z ) = 0 .
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. □

6. 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 approximation P ^ ( 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:
z ( τ ) = x ( τ ) + ε y ( τ ) = x ( τ ) + σ σ i ( τ ) 1 σ x ( τ ) = x ( τ ) + σ i ( τ ) x ( τ ) = σ i ( τ ) .
There are various ways to construct probability distribution P ^ ( 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
G ( i ) = C b ( σ i ) exp 2 σ 0 σ i a ( x ) b ( x ) d x .
Taking into account normalization condition, we can write
P ^ ( i ) = G ( i ) i = 0 G ( i ) .
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 approximation R ^ 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:
R ^ 2 ( n 1 , n 2 ) = i = 0 R ( n 1 , n 2 , σ i ) P ^ ( i ) ,
where R ( n 1 , n 2 , x ) is calculated using Equation (13).
So, basing on these two types of approximations, we can write two types R ^ 1 ( n ) and R ^ 2 ( n ) of approximations for the stationary probability distribution R ( n ) of the number of busy servers as follows:
R ^ k ( n ) = ν = 0 n R k ( ν , n ν ) , k = 1 , 2 .
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
Δ = max 0 i < ν = 0 i P 1 ( ν ) P 2 ( ν )
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 from P ^ ( i ) , R ^ 1 ( n ) or R ^ 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
μ 1 = 0.5 , , μ 2 = 3 , q 1 = 0.4 , q 2 = 0.6 .
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 approximations R ^ 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).

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 μ 2 = m. Some calls (about 20%) require long time to be serviced. We model them using rate μ 1 = 0.2 m for exponential distribution. So, we can generally model service time in our system as hyper-exponential with parameters μ 1 = 0.2 m, μ 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 i ¯ · C E where C E is some constant and i ¯ = 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:
K ( m ) = m · C s + C E i = 0 i · P i + C L i = w + 1 P i .
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.

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 Δ 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.

Author Contributions

Conceptualization, A.N.; Formal analysis, A.M.; Investigation, A.M., A.N. and S.P.; Methodology, A.N.; Software, A.M.; Writing—original draft, S.P.; Writing—review & editing, A.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Aguir, S.; Karaesmen, F.; Askin, O.Z.; Chauvet, F. The impact of retrials on call center performance. OR Spektrum 2004, 26, 353–376. [Google Scholar] [CrossRef] [Green Version]
  2. Stolletz, R. Performance Analysis and Optimization of Inbound Call Centers; Springer: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  3. Avrachenkov, K.; Yechiali, U. Retrial networks with finite buffers and their application to internet data traffic. Probab. Eng. Inf. Sci. 2008, 22, 519–536. [Google Scholar] [CrossRef] [Green Version]
  4. Roszik, J.; Sztrik, J.; Kim, C.S. Retrial queues in the performance modelling of cellular mobile networks using MOSEL. Int. J. Simul. 2005, 6, 38–47. [Google Scholar]
  5. Artalejo, J.R.; Lopez-Herrero, M.J. Cellular mobile networks with repeated calls operating in random environment. Comput. Oper. Res. 2010, 37, 1158–1166. [Google Scholar] [CrossRef]
  6. Ke, X.; Shi, L.; Guo, W.; Chen, D. Multi-Dimensional Traffic Congestion Detection Based on Fusion of Visual Features and Convolutional Neural Network. IEEE Trans. Intell. Trans. Syst. 2019, 20, 2157–2170. [Google Scholar] [CrossRef]
  7. Nazarov, A.A.; Moiseev, A.N. Distributed System of Processing of Data of Physical Experiments. Rus. Phys. J. 2014, 57, 984–990. [Google Scholar] [CrossRef]
  8. Phung-Duc, T. Retrial Queueing Models: A Survey on Theory and Applications. arXiv 2019, arXiv:1906.09560v1. [Google Scholar]
  9. Falin, G.I.; Templeton, J.G.C. Retrial Queues; Chapman & Hall: London, UK, 1997. [Google Scholar]
  10. Falin, G.I. A survey of retrial queues. Queueing Syst. 1990, 7, 127–168. [Google Scholar] [CrossRef]
  11. Artalejo, J.R.; Gomez-Corral, A. Retrial Queueing Systems: A Computational Approach; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  12. Shekhar, C.; Raina, A.A.; Kumar, A. A brief review on retrial queue: Progress in 2010–2015. Int. J. Appl. Sci. Eng. Res. 2016, 5, 324–336. [Google Scholar]
  13. Kim, C.; Mushko, V.; Dudin, A. Computation of the steady state distribution for multi-server retrial queues with phase type service process. Ann. Oper. Res. 2012, 201, 307–323. [Google Scholar] [CrossRef]
  14. Phung-Duc, T.; Masuyama, H.; Kasahara, S.; Takahashi, Y. A matrix continued fraction approach to multiserver retrial queues. Ann. Oper. Res. 2013, 203, 161–183. [Google Scholar] [CrossRef] [Green Version]
  15. Neuts, M.; Rao, B. Numerical investigation of a multiserver retrial mode. Queueing Syst. 2002, 7, 169–189. [Google Scholar] [CrossRef]
  16. Artalejo, J.; Pozo, M. Numerical calculation of the stationary distribution of the main multiserver retrial queue. Ann. Oper. Res. 2002, 116, 41–56. [Google Scholar] [CrossRef]
  17. Ridder, A. Fast simulation of retrial queues. In Proceedings of the 3rd Workshop on Rare Event Simulation and Related Combinatorial Optimization Problems, Pisa, Italy, 5–6 October 2000. [Google Scholar]
  18. Moiseeva, E.; Nazarov, A. Asymptotic analysis of RQ-systems M/M/1 on heavy load condition. In Proceedings of the IV International Conference Problems of Cybernetics and Informatics, Baku, Azerbaijan, 12–14 September 2012; pp. 164–166. [Google Scholar]
  19. Fedorova, E. Quasi-geometric and gamma approximation for retrial queueing systems. Commun. Comput. Inf. Sci. 2014, 487, 123–136. [Google Scholar]
  20. Nazarov, A.; Chernikova, Y. Gaussian approximations of probabilities distribution of states of the retrial queueing system with r-persistent exclusion of alternative customers. Commun. Comput. Inf. Sci. 2015, 564, 200–209. [Google Scholar]
  21. Lapatin, I.; Nazarov, A. Asymptotic Analysis of the Output Process in Retrial Queue with Markov-Modulated Poisson Input Under Low Rate of Retrials Condition. Commun. Comput. Inf. Sci. 2019, 1141, 315–324. [Google Scholar]
  22. Nazarov, A.; Phung-Duc, T.; Paul, S.; Lizura, O. Asymptotic-Diffusion Analysis for Retrial Queue with Batch Poisson Input and Multiple Types of Outgoing Calls. Lect. Notes Comput. Sci. 2019, 11965, 207–222. [Google Scholar]
  23. Moiseev, A.; Demin, A.; Dorofeev, V.; Sorokin, V. Discrete-event approach to simulation of queueing networks. Key Eng. Mater. 2016, 685, 939–942. [Google Scholar] [CrossRef]
  24. Mesheryakov, R.; Moiseev, A.; Demin, A.; Dorofeev, V.; Sorokin, V. Using parallel computing in queueing network simulation. Key Eng. Mater. 2016, 685, 943–947. [Google Scholar] [CrossRef]
Figure 1. Values of the total cost K ( m ) for various values of bound w.
Figure 1. Values of the total cost K ( m ) for various values of bound w.
Mathematics 08 00531 g001
Table 1. Values of Kolmogorov distances Δ k of approximations R ^ k ( n ) from Equation (54) for various values of σ and ρ .
Table 1. Values of Kolmogorov distances Δ k of approximations R ^ k ( n ) from Equation (54) for various values of σ and ρ .
ρ σ 210.50.20.1
0.6 Δ 1 0.0440.0320.0200.0110.006
Δ 2 0.0590.0530.0420.0190.008
0.7 Δ 1 0.0650.0460.0300.0150.009
Δ 2 0.0830.0630.0380.0130.003
0.8 Δ 1 0.0780.0560.0370.0190.011
Δ 2 0.0750.0500.0250.0050.002
0.9 Δ 1 0.0730.0490.0300.0150.008
Δ 2 0.0450.0210.0080.0010.001
Table 2. Values of Kolmogorov distance of approximation P ^ for various values of σ and ρ .
Table 2. Values of Kolmogorov distance of approximation P ^ for various values of σ and ρ .
σ 210.50.20.1
ρ = 0.6 0.0850.1290.1410.0820.039
ρ = 0.7 0.1220.1220.0810.0270.014
ρ = 0.8 0.1040.0700.0320.0110.014
ρ = 0.9 0.0520.0220.0030.0120.010

Share and Cite

MDPI and ACS Style

Moiseev, A.; Nazarov, A.; Paul, S. Asymptotic Diffusion Analysis of Multi-Server Retrial Queue with Hyper-Exponential Service. Mathematics 2020, 8, 531. https://doi.org/10.3390/math8040531

AMA Style

Moiseev A, Nazarov A, Paul S. Asymptotic Diffusion Analysis of Multi-Server Retrial Queue with Hyper-Exponential Service. Mathematics. 2020; 8(4):531. https://doi.org/10.3390/math8040531

Chicago/Turabian Style

Moiseev, Alexander, Anatoly Nazarov, and Svetlana Paul. 2020. "Asymptotic Diffusion Analysis of Multi-Server Retrial Queue with Hyper-Exponential Service" Mathematics 8, no. 4: 531. https://doi.org/10.3390/math8040531

APA Style

Moiseev, A., Nazarov, A., & Paul, S. (2020). Asymptotic Diffusion Analysis of Multi-Server Retrial Queue with Hyper-Exponential Service. Mathematics, 8(4), 531. https://doi.org/10.3390/math8040531

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop