Next Article in Journal
Forecasting Spatially-Distributed Urban Traffic Volumes via Multi-Target LSTM-Based Neural Network Regressor
Next Article in Special Issue
COVID-19 Spatial Diffusion: A Markovian Agent-Based Model
Previous Article in Journal
Exact Solutions of Bernoulli and Logistic Fractional Differential Equations with Power Law Coefficients
Previous Article in Special Issue
Short-Scale Stochastic Solar Energy Models: A Datacenter Use Case
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Diffusion Limit of Multi-Server Retrial Queue with Setup Time

1
Institute of Applied Mathematics and Computer Science, Tomsk State University, 36 Lenin Ave., 634050 Tomsk, Russia
2
Department of Policy and Planning Sciences, Faculty of Engineering, Information and Systems, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan
3
Public Policy Program, VNU Vietnam Japan University, My Dinh Campus, Nam Tu Liem, Hanoi, Vietnam
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(12), 2232; https://doi.org/10.3390/math8122232
Submission received: 1 October 2020 / Revised: 10 December 2020 / Accepted: 14 December 2020 / Published: 16 December 2020
(This article belongs to the Special Issue Queue and Stochastic Models for Operations Research)

Abstract

:
In the paper, we consider a multi-server retrial queueing system with setup time which is motivated by applications in power-saving data centers with the ON-OFF policy, where an idle server is immediately turned off and an off server is set up upon arrival of a customer. Customers that find all the servers busy join the orbit and retry for service after an exponentially distributed time. For this model, we derive the stability condition which depends on the setup time and turns out to be more strict than that of the corresponding model with an infinite buffer which is independent of the setup time. We propose asymptotic methods to analyze the system under the condition that the delay in the orbit is extremely long. We show that the scaled-number of customers in the orbit converges to a diffusion process. Using this diffusion limit, we obtain approximations for the steady-state probability distribution of the number of busy servers and that of the number of customers in the orbit. We verify the accuracy of the approximations by simulations and numerical analysis. Numerical results show that the retrial system under the limiting condition consumes more energy than that with an infinite buffer in front of the servers.

1. Introduction

Presently, cloud computing is used by many companies and individuals. In our everyday business, we use many kinds of sharing services such as Dropbox, Slack, Overleaf, etc. which are based on cloud computing. Cloud computing is even more important under the current situation where social distancing is encouraged to combat COVID-19, leading to remote work for a large portion of our society. Cloud computing is supported by data centers in which a huge number of servers are available. These servers consume a huge amount of energy and thus saving-energy is crucial in management of data centers and cloud computing. Since user traffic has peak-on and peak-off nature, it is desired that more server resource is allocated in peak-on period and less resource is allocated in peak-off period. Furthermore, because user traffic stochastically varies, we need control policies that add more resources when the workload is large and release resources when the workload is small. The same mechanism is also observed in 5G networks with network functions virtualization (NFV) [1,2]. Motivated by these situations, many power-saving mechanisms were proposed [3,4,5].
One of the most natural mechanisms is the ON-OFF policy which turns OFF a server once it has no job to handle and turns on an OFF server again once a job arrives. However, an OFF server cannot be active immediately to serve the waiting job but needs a setup time, during which the server cannot process the job but consumes energy. Queues with setup time are appropriate models for these situations. A server is turned off when it has no job to process while it is turned on again when waiting jobs are available. These models are challenging because the underlying Markov chain has a non-homogeneous structure, i.e., the number of setup servers may depend on the number of jobs in the system. Gandhi et al. [4] study this problem and propose analysing the model using a renewal reward approach. Phung-Duc [6] analyzes the same model using a generating function approach and a matrix analytic method. Furthermore, structure of the optimal policy was studied by Maccio and Down [5,7].
Furthermore, from the user point of view, the request might be blocked if a resource is not allocated, i.e., all the servers are busy upon arrival. In practice, in case the request is blocked, a protocol will reconnect after some time. This time is random and depends on the number of retries. In some application, the retrial interval is doubled for two consecutive blockings.
Motivated by the above applications, we consider a multiserver retrial queues with setup time. In this model, the server is switched off immediately once it finishes serving a job and no new job is available. An off server will be switched on upon arrival of a job and the server changes to the setup state. After some setup time, the server becomes active so that it can serves the job. If all servers are occupied (serving a job or in setup), the arriving job is blocked and makes a retrial (joins the orbit) after an exponentially distributed time. A retrial job behaves the same as a fresh one. The model was first proposed and numerically analyzed in [8], where a non-trivial sufficient stability condition was derived. Some analytical results for single server case are obtained in [9]. In this paper, we study the model in depth, in which our focus is not a numerical solution but to obtain analytical solution under a specific regime. We consider the situation where the retrial time interval is relatively long. In such a regime, the number of retrial jobs explodes. However, using a proper scaling, we can obtain the explicit distribution for the number of jobs in the orbit. The explicit solution is then used to approximate the distribution of the number of customers in the orbit and that of the states of the servers.
It should be noted that multiserver retrial queues without setup time whose underlying Markov chain is two-dimensional are already very difficult to obtain an analytical solution [10,11]. This paper combines two challenging features, i.e., setup time and retrial. As a result, our model in this paper is formulated by a three-dimensional Markov chain and thus is even more challenging to obtain analytic results. Thus, we consider the system under the slow retrial asymptotic regime.
The slow retrial asymptotic regime was studied by Cohen [12], where the first order asymptotic (a law of large number) for multiserver retrial queue without setup time was obtained in the stationary regime. This result is extended in our paper to the model with setup time in non-stationary regime. Our result states that the scaled number of jobs in the orbit converges to a deterministic process which is characterized by an ordinary differential equation (ODE). Furthermore, we study the second order asymptotics, where the scaled number of jobs in the orbit weakly converges to a diffusion process. The limiting results are then used to obtain the approximation of performance measures. As related works, the second order asymptotic result for multiserver retrial queue without setup was derived in [13] (pp. 55–59). Some recent development of asymptotic analysis of retrial queues can be found in [14,15,16]. Furthermore, we refer to the books [13,17] for basic results of retrial queueing systems. Difussion limits for queueing models have been extensively studied [18,19,20,21,22,23,24]. The methodological difference is that we are based on characteristic function while [18,19,20,21,22,23,24] are based on sample path equations or other techniques.
The rest of our paper is organized as follows. Mathematical description of the problem is presented in Section 2. In Section 3, we obtain main equations for the problem solution. The first stage of the asymptotic analysis is made in Section 4, where we derive function a ( x ) which will be used for further diffusion analysis. Based on the results of the first stage, we obtain the condition for the existence of the steady-state regime for our model in Section 5. The second stage of the asymptotic analysis is made in Section 6, where we derive function b ( x ) which is used together with a ( x ) for the asymptotic diffusion analysis in Section 7. Using obtained continuous distributions, we build approximations for discrete distributions under study in Section 8. In Section 9, we analyze applicability areas of obtained approximations using numerical methods and simulation. Concluding remarks are presented in Section 10.

2. Mathematical Model

We consider a multi-server retrial queue with Poisson arrivals with intensity λ , infinite-capacity orbit and N servers. On the arrival of a customer, if all servers are busy, the customer goes to the orbit where he or she stays for a random time exponentially distributed with parameter σ (mean 1 / σ ) and then tries to get a service again. If there is a free server, the customer occupies it. Each server serves a customer in two phases: the first one is a setup time which is exponentially distributed with parameter μ 1 and the second one is a real service which is exponentially distributed with parameter μ 2 . After the service completed at the second phase, the customer leaves the system.
The service in the system has a specific feature: when one customer completes its service at the second phase and make its server free, and if we have another customer which is served at the first phase at another server, this customer (which was served at the first phase) immediately goes to the server that just becomes free and starts its service at the second phase. Therefore, its total time of servicing becomes less.
Let us denote:
  • n 1 ( t ) is the number of servers that are working at the first phase at the instant t;
  • n 2 ( t ) is the number of servers that are working at the second phase at the instant t;
  • i ( t ) is the number of customers in the orbit at the instant t;
  • P ( n 1 , n 2 , i , t ) = P { n 1 ( t ) = n 1 , n 2 ( t ) = n 2 , i ( t ) = i } is the joint probability distribution of the stochastic process { n 1 ( t ) , n 2 ( t ) , i ( t ) } .
The goal of the study is to obtain the probability distribution P ( n 1 , n 2 , i , t ) P ( n 1 , n 2 , i ) for the process { n 1 ( t ) , n 2 ( t ) , i ( t ) } in the steady-state regime. The problem is solved using the method of asymptotic diffusion analysis [16,25] under an asymptotic condition of long delay in the orbit: σ 0 (mean retrial interval 1 / σ ).

3. Main Equations

Because the process { n 1 ( t ) , n 2 ( t ) , i ( t ) } is a three-dimensional Markov chain, we can write the following equalities. For n 1 + n 2 = 0 :
P ( 0 , 0 , i , t + Δ t ) = P ( 0 , 0 , i , t ) ( 1 λ Δ t ) ( 1 i σ Δ t ) + P ( 0 , 1 , i , t ) μ 2 Δ t + o ( Δ t ) .
For 1 n 1 + n 2 N 1 :
P ( n 1 , n 2 , i , t + Δ t ) = P ( n 1 , n 2 , i , t ) ( 1 λ Δ t ) ( 1 n 1 μ 1 Δ t ) ( 1 n 2 μ 2 Δ t ) ( 1 i σ Δ t ) + P ( n 1 1 , n 2 , i , t ) λ Δ t + P ( n 1 + 1 , n 2 1 , i , t ) ( n 1 + 1 ) μ 1 Δ t + P ( n 1 + 1 , n 2 , i , t ) n 2 μ 2 Δ t + P ( n 1 1 , n 2 , i + 1 , t ) ( i + 1 ) σ Δ t + o ( Δ t ) .
For n 1 + n 2 = N :
P ( n 1 , n 2 , i , t + Δ t ) = P ( n 1 , n 2 , i , t ) ( 1 λ Δ t ) ( 1 n 1 μ 1 Δ t ) ( 1 n 2 μ 2 Δ t ) + P ( n 1 1 , n 2 , i , t ) λ Δ t + P ( n 1 , n 2 , i 1 , t ) λ Δ t + P ( n 1 + 1 , n 2 1 , i , t ) ( n 1 + 1 ) μ 1 Δ t + P ( n 1 1 , n 2 , i + 1 , t ) ( i + 1 ) σ Δ t + o ( Δ t ) .
We derive the following system of differential balance equations from here. For n 1 + n 2 = 0 :
P ( 0 , 0 , i , t ) t = ( λ + i σ ) P ( 0 , 0 , i , t ) + μ 2 P ( 0 , 1 , i , t ) .
For 1 n 1 + n 2 N 1 :
P ( n 1 , n 2 , i , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 + i σ ) P ( n 1 , n 2 , i , t ) + λ P ( n 1 1 , n 2 , i , t ) + ( n 1 + 1 ) μ 1 P ( n 1 + 1 , n 2 1 , i , t ) + n 2 μ 2 P ( n 1 + 1 , n 2 , i , t ) + ( i + 1 ) σ P ( n 1 1 , n 2 , i + 1 , t ) .
For n 1 + n 2 = N :
P ( n 1 , n 2 , i , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 ) P ( n 1 , n 2 , i , t ) + λ P ( n 1 1 , n 2 , i , t ) + λ P ( n 1 , n 2 , i 1 , t ) + ( n 1 + 1 ) μ 1 P ( n 1 + 1 , n 2 1 , i , t ) + ( i + 1 ) σ P ( n 1 1 , n 2 , i + 1 , t ) .
Denoting j = 1 and using partial characteristic functions
H ( n 1 , n 2 , u , t ) = i = 0 e j u i P ( n 1 , n 2 , i , t ) ,
we obtain the following system of differential equations.
For n 1 + n 2 = 0 :
H ( 0 , 0 , u , t ) t = λ H ( 0 , 0 , u , t ) + μ 2 H ( 0 , 1 , u , t ) + j σ H ( 0 , 0 , u , t ) u .
For 1 n 1 + n 2 N 1 :
H ( n 1 , n 2 , u , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + λ H ( n 1 1 , n 2 , u , t ) + ( n 1 + 1 ) μ 1 H ( n 1 + 1 , n 2 1 , u , t ) + n 2 μ 2 H ( n 1 + 1 , n 2 , u , t ) + j σ H ( n 1 , n 2 , u , t ) u e j u j σ H ( n 1 1 , n 2 , u , t ) u .
For n 1 + n 2 = N :
H ( n 1 , n 2 , u , t ) t = ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + e j u λ H ( n 1 , n 2 , u , t ) + λ H ( n 1 1 , n 2 , u , t ) + ( n 1 + 1 ) μ 1 H ( n 1 + 1 , n 2 1 , u , t ) e j u j σ H ( n 1 1 , n 2 , u , t ) u .
Let us use linear finite difference operators A , B , I 0 , I 1 to rewrite the system (1)–(3) in the following compact form:
H ( u , t ) t = A + e j u λ B H ( u , t ) + I 0 e j u I 1 j σ H ( u , t ) u .
Here H ( u , t ) is ( N + 1 ) × ( N + 1 ) top-left triangle matrix with entries equal to H ( n 1 , n 2 , u , t ) for indices n 1 0 , n 2 0 , n 1 + n 2 N and equal to zero for other numbers n 1 , n 2 . Operators in (4) are defined as follows:
A H ( u , t ) n 1 , n 2 = λ H ( 0 , 0 , u , t ) + μ 2 H ( 0 , 1 , u , t ) for   n 1 + n 2 = 0 , ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + λ H ( n 1 1 , n 2 , u , t ) + ( n 1 + 1 ) μ 1 H ( n 1 + 1 , n 2 1 , u , t ) + n 2 μ 2 H ( n 1 + 1 , n 2 , u , t ) for   1 n 1 + n 2 N 1 , ( λ + n 1 μ 1 + n 2 μ 2 ) H ( n 1 , n 2 , u , t ) + λ H ( n 1 1 , n 2 , u , t ) + ( n 1 + 1 ) μ 1 H ( n 1 + 1 , n 2 1 , u , t ) for   n 1 + n 2 = N ,
B H ( u , t ) n 1 , n 2 = H ( n 1 , n 2 , u , t ) for   n 1 + n 2 = N , 0 otherwise   ,
I 0 H ( u , t ) n 1 , n 2 = H ( n 1 , n 2 , u , t ) for   n 1 + n 2 N 1 , 0 otherwise   ,
I 1 H ( u , t ) n 1 , n 2 = H ( n 1 1 , n 2 , u , t ) .
Summing up all equations (1)–(3), we derive
t n 1 + n 2 N H ( n 1 , n 2 , u , t ) = e j u 1 λ n 1 + n 2 = N H ( n 1 , n 2 , u , t ) + e j u n 1 + n 2 N 1 j σ H ( n 1 , n 2 , u , t ) u .
Let us denote the summing operator for indices n 1 + n 2 = N by E 1 , the summing operator for indices n 1 + n 2 N 1 by E 2 , and the total summing operator by E = E 1 + E 2 . Then Equation (9) together with (4) may be rewritten in the form
H ( u , t ) t = A + e j u λ B H ( u , t ) + I 0 e j u I 1 j σ H ( u , t ) u , t E H ( u , t ) = e j u 1 λ E 1 H ( u , t ) + e j u j σ u E 2 H ( u , t ) .
This system of equations is the main for the study of the current paper. We will solve it using the method of asymptotic diffusion analysis [16,25] in the next sections.

4. First Stage of Asymptotic Analysis

Let us find the solution of system (10) using the method of asymptotic analysis [16,25] under asymptotic condition of long delay in the orbit: σ 0 . To do this, let us make the following substitutions:
σ = ε , τ = ε t , u = ε w , H ( u , t ) = F ( w , τ , ε ) .
Then we obtain the following system:
ε F ( w , τ , ε ) τ = A + e j ε w λ B F ( w , τ , ε ) + I 0 e j ε w I 1 j F ( w , τ , ε ) w , ε τ E F ( w , τ , ε ) = e j ε w 1 λ E 1 F ( w , τ , ε ) + j e j ε w w E 2 F ( w , τ , ε ) .
We can prove the following statement.
Theorem 1.
Under the asymptotic condition σ 0 , the following equality holds.
lim σ 0 E e j w σ i τ σ = e j w x ( τ ) .
Here the scalar function x ( τ ) is a solution of ordinary differential equation
x ( τ ) = a ( x ) = λ E 1 R x E 2 R ,
where R = R ( x ) is a left-top triangle matrix which is a solution of homogeneous linear system
A + λ B + x ( I 1 I 0 ) R = 0
and satisfies the normalization condition
ER = n 1 + n 2 N R ( n 1 , n 2 ) = 1 .
Proof. 
Denoting lim ε 0 F ( w , τ , ε ) = F ( w , τ ) and making asymptotic transition ε 0 in (11), we derive
( A + λ B ) F ( w , τ ) + ( I 0 I 1 ) j F ( w , τ ) w = 0 , τ E F ( w , τ ) = j w λ E 1 F ( w , τ ) + j w E 2 F ( w , τ ) .
We find the solution of this system in the form
F ( w , τ ) = e j w x ( τ ) R ,
where R is a top-left triangle matrix with non-zero entries R ( n 1 , n 2 ) for n 1 + n 2 N , and x ( τ ) is a scalar function which has a meaning of asymptotic (while σ 0 ) value of the normalized number of customers in the orbit σ i τ σ . Entries R ( n 1 , n 2 ) have a meaning of asymptotic (while σ 0 ) probabilities that n 1 servers are working at the first phase and n 2 servers are working at the second phase. Substituting (17) into (16), we obtain equalities (14) and (13). Because entries R ( n 1 , n 2 ) of matrix R are probabilities, they satisfy (15). Since the scalar function x ( τ ) is an asymptotic value of the normalized number of customers in the orbit σ i τ σ , equality (12) is true. The theorem is proved. □
Probability distribution R is a solution of system of Equation (14) whose coefficients depend on x, then R is a matrix function: R = R ( x ) . Due to this, we can rewrite (13) as the following expression:
a ( x ) = λ E 1 R ( x ) x E 2 R ( x ) ,
which determines the function a ( x ) . This function is very important for the study due to the following two reasons. The first one is that we showed that x ( τ ) = a ( x ) . Therefore, function a ( x ) characterizes dynamic properties of the process x ( τ ) . The second reason is that a ( x ) will be used in Section 7 as a drift coefficient for the diffusion process which determines the asymptotic number of customers in the orbit.
Remark 1.
For σ 0 (the mean retrial time 1 / σ ), the number of customers in the orbit explosively increases, Theorem 1 shows that i ( t ) / ( 1 / σ ) converges to a deterministic process independent of σ and thus has the same order as 1 / σ . Theorem 1 can be interpreted as the law of large numbers.

5. Existence of Steady-State Regime

Considering function a ( x ) which is determined as (18), let us derive a condition of steady-state regime existence.
Due to (18), a ( 0 ) > 0 and values of the process x ( τ ) are growing in the neighborhood of point x = 0 . If a ( x ) > 0 for all x, then steady-state regime does not exist in the considered retrial queue. The inequality lim x a ( x ) < 0 is the necessary condition for the existence of the steady-state regime.
Let us prove the following statement.
Theorem 2.
Steady-state regime in the considered retrial queue exists if and only if
λ < N μ 2 ( N 1 ) μ 2 + μ 1 N μ 2 + μ 1 .
Proof. 
Sufficiency of condition (19) was proved in [8]. Let us prove its necessity. We rewrite (18) in the form:
a ( x ) = λ n 1 + n 2 = N R ( n 1 , n 2 , x ) x n 1 + n 2 N 1 R ( n 1 , n 2 , x ) .
To determine value of lim x a ( x ) , let us find limit properties of probabilities R ( n 1 , n 2 , x ) under x . Transition intensities for states ( n 1 , n 2 ) are depicted in Figure 1 (for simplicity, we draw the graph for partial case of N = 4 ). In the considered system, we have the following possible transitions:
  • from state ( n 1 , n 2 ) to state ( n 1 + 1 , n 2 ) with intensity ( λ + x ) , where λ is an intensity of arrivals and x is an asymptotic intensity of retrials from orbit,
  • from state ( n 1 , n 2 ) to state ( n 1 1 , n 2 + 1 ) with intensity n 1 μ 1 when one of n 1 customers completes service at the first phase and goes to the second one,
  • from state ( n 1 , n 2 ) to state ( n 1 1 , n 2 ) with intensity n 2 μ 2 for n 1 > 0 when one of n 2 customers completes its service at the second phase and leaves the system but one customer served at the first phase instantly goes to the second phase,
  • from state ( 0 , n 2 ) to state ( 0 , n 2 1 ) with intensity n 2 μ 2 when one of n 2 customers completes its service at the second phase and leaves the system.
Using the method cut of graph, we derive equations for probabilities R ( n 1 , n 2 , x ) . We can derive them from system of (14) and (15) but it is clearer how they are obtained if using cut of graph. It should be noted that the system of (14) and (15) is equivalent to the system of equations for finding the stationary distribution of the Markov chain with the transition diagram as in Figure 1. This Markov chain represents the M/M/N/N queueing system with setup time with the arrival rate given by λ + x .
For diagonal cuts when n 1 + n 2 = n N 1 we can write
( λ + x ) n 1 + n 2 = n R ( n 1 , n 2 , x ) = μ 2 n 1 = 1 n + 1 n 2 R ( n + 1 n 2 , n 2 , x ) .
Therefore, under the limit condition x all the probabilities R ( n 1 , n 2 , x ) are equal to zero when n 1 + n 2 N 1 . Then due to the normalization requirement, we have
lim x n 1 + n 2 = N R ( n 1 , n 2 , x ) = 1 ,
and, so, all non-zero probabilities R ( n 1 , n 2 ) are located on the diagonal n 1 + n 2 = N . Moreover, from (21), it follows that all the sums n 1 + n 2 = n R ( n 1 , n 2 , x ) are infinitesimals of order 1 x N n for n N 1 . For horizontal cuts we can write the following equalities:
N μ 2 R ( 0 , N , x ) = μ 1 R ( 1 , N 1 , x ) , n μ 2 R ( 0 , n 2 , x ) = μ 1 n 1 = 1 N ( n 2 1 ) n 1 R ( n 1 , n 2 1 , x ) for   n 2 N 1 .
Points ( 0 , n 2 ) for n 2 N 1 lay on diagonals of the graph, therefore, all R ( 0 , n 2 , x ) for n 2 N 1 are infinitesimals of order 1 x N n 2 . Then all the sums n 1 = 1 N ( n 2 1 ) R ( n 1 , n 2 1 , x ) are infinitesimals of order 1 x N n 2 too. Hence, for n 2 = N 1 , the sum R ( 1 , N 2 , x ) + R ( 2 , N 2 , x ) and probability R ( 2 , N 2 , x ) have the infinitesimal order of 1 x . For n 2 N 2 , all the sums n 1 = 1 N ( n 2 1 ) R ( n 1 , n 2 1 , x ) and probabilities R ( N ( n 2 1 ) , n 2 1 , x ) are infinitesimals of order 1 x 2 or higher. Therefore, under the limit x , only two probabilities R ( 0 , N ) and R ( 1 , N 1 ) are not equal to zero. Due to the normalization condition, we can write the following:
N μ 2 R ( 0 , N ) = μ 1 R ( 1 , N 1 ) , R ( 0 , N ) + R ( 1 , N 1 ) = 1 .
Solution of this system is as follows:
R ( 0 , N ) = μ 1 N μ 2 + μ 1 , R ( 1 , N 1 ) = N μ 2 N μ 2 + μ 1 .
Let us consider the sum R ( 1 , N 2 , x ) + R ( 2 , N 2 , x ) . As we find before, it is an infinitesimal of order 1 x . Let us write the equation for the probability R ( 2 , N 2 , x ) :
2 μ 1 + ( N 2 ) μ 2 R ( 2 , N 2 , x ) + ( λ + x ) R ( 1 , N 2 , x ) + 3 μ 1 R ( 3 , N 3 , x ) = 0 .
Here R ( 3 , N 3 , x ) is an infinitesimal of order 1 x 2 , hence, the probability R ( 1 , N 2 , x ) , which is multiplied here by x, also, is an infinitesimal of order 1 x 2 because only in this case equality (24) is true and sum R ( 1 , N 2 , x ) + R ( 2 , N 2 , x ) is an infinitesimal of order 1 x . Therefore, for n 1 + n 2 N 1 , only one probability R ( 0 , N 1 , x ) has infinitesimal order 1 x , all other probabilities R ( n 1 , n 2 , x ) have higher infinitesimal order. Let us write equation for the probability R ( 0 , N 1 , x ) as follows:
λ + x + ( N 1 ) μ 2 R ( 0 , N 1 , x ) = N μ 2 R ( 0 , N , x ) + ( N 1 ) μ 2 R ( 1 , N 1 , x ) + μ 1 R ( 1 , N 2 , x ) .
Taking here the limit x , we obtain
lim x x R ( 0 , N 1 , x ) = N μ 2 R ( 0 , N ) + ( N 1 ) μ 2 R ( 1 , N 1 ) .
Due to (23), we can rewrite it as follows:
lim x x R ( 0 , N 1 , x ) = N μ 2 ( N 1 ) μ 2 + μ 1 N μ 2 + μ 1 .
Let us come back to the function a ( x ) = x ( τ ) which determines derivative of the process x ( τ ) = σ i ( τ σ ) under asymptotic condition σ 0 . Let us find limit of the function a ( x ) from (20) under condition x . Due to (22) and (25), we can write
lim x a ( x ) = λ N μ 2 ( N 1 ) μ 2 + μ 1 N μ 2 + μ 1 .
Because lim x a ( x ) < 0 is necessary for the existence of steady-state regime, the condition (19) is true. The theorem is proved. □
Remark 2.
In [6], the M/M/N/Setup model with infinite buffer is analyzed and the stability condition is given λ < N μ 2 . Theorem 2 shows that the stability condition for the corresponding retrial model is more strict than that of the infinite buffer counterpart.

6. Second Stage of Asymptotic Analysis

In Section 4, we made the first stage of the asymptotic analysis and, similar to the law of large numbers, we obtain equality (12) which determines the convergence of characteristic function of the process σ i τ σ to the deterministic function x ( τ ) under condition σ 0 . Now, let us perform the second stage of the asymptotic analysis to obtain more detailed parameters of the convergence.
Let us make the following substitution in system (10)
H ( u , t ) = H ( 1 ) ( u , t ) e j u σ x ( σ t ) .
We obtain
H ( 1 ) ( u , t ) t + j u a ( x ) H ( 1 ) ( u , t ) = A + e j u λ B + x e j u I 1 I 0 H ( 1 ) ( u , t ) + I 0 e j u I 1 j σ H ( 1 ) ( u , t ) u , t E H ( 1 ) ( u , t ) + j u a ( x ) E H ( 1 ) ( u , t ) = e j u 1 λ E 1 x e j u E 2 H ( 1 ) ( u , t ) + j σ e j u u E 2 H ( 1 ) ( u , t ) .
Substitution (26) is made to go to a centered process for i ( t ) : the function H ( 1 ) ( u , t ) is a matrix characteristic function of the centered process i ( t ) 1 σ x ( σ t ) , where function x ( τ ) have been obtained at the first stage of the asymptotic analysis in Section 4.
Denoting σ = ε 2 and making in (27) substitutions
τ = ε 2 t , u = ε w , H ( 1 ) ( u , t ) = F ( 1 ) ( w , τ , ε ) ,
we derive
ε 2 F ( 1 ) ( w , τ , ε ) τ + j ε w a ( x ) F ( 1 ) ( w , τ , ε ) = A + e j ε w λ B + x e j ε w I 1 I 0 F ( 1 ) ( w , τ , ε ) + j ε I 0 e j ε w I 1 F ( 1 ) ( w , τ , ε ) w , ε 2 τ E F ( 1 ) ( w , τ , ε ) + j ε w a ( x ) E F ( 1 ) ( w , τ , ε ) = e j ε w 1 λ E 1 x e j ε w E 2 F ( 1 ) ( w , τ , ε ) + e j ε w j ε w E 2 F ( 1 ) ( w , τ , ε ) .
We can prove the following statement.
Theorem 3.
Let Φ ( w , τ ) be an asymptotic (as σ 0 ) characteristic function of normalized and centered process σ i τ σ 1 σ x ( τ ) :
Φ ( w , τ ) = lim σ 0 E exp j w σ i τ σ 1 σ x ( τ ) .
Then it satisfies the equation
Φ ( w , τ ) τ = a ( x ) w Φ ( w , τ ) w + b ( x ) ( j w ) 2 2 Φ ( w , τ ) ,
where function a ( x ) is determined by (18), b ( x ) has the form
b ( x ) = a ( x ) + 2 ( λ + x ) E 2 g ( x ) + x E 2 R ( x ) ,
and top-left triangle matrix g ( x ) is a particular solution of the system of equations
A + λ B + x I 1 I 0 g ( x ) = a ( x ) R ( x ) + x I 1 λ B R ( x ) ,
and satisfies additional condition
E g ( x ) = 0 .
Proof. 
We rewrite the first equation of (29) with precision up to infinitesimals of order ε 2 :
j ε w a ( x ) F ( 1 ) ( w , τ , ε ) = A + λ B + x I 1 I 0 + j ε w λ B x I 1 F ( 1 ) ( w , τ , ε ) + j ε I 0 I 1 F ( 1 ) ( w , τ , ε ) w + O ε 2 ,
and let us write its solution in the form of expansion:
F ( 1 ) ( w , τ , ε ) = Φ ( w , τ ) R ( x ) + j ε w f ( x ) + O ε 2 ,
where f ( x ) is some matrix function whose expression we will obtain later. Substituting (36) into (35), we obtain
j ε w a ( x ) R ( x ) = A + λ B + x I 0 I 1 R ( x ) + j ε w f ( x ) + j ε w λ B x I 1 R ( x ) + I 0 I 1 R ( x ) j ε Φ ( w , τ ) / w Φ ( w , τ ) + O ε 2 .
Taking into account (14) and dividing by j ε w , under the asymptotic condition ε 0 , we obtain the following equation for the matrix function f ( x ) :
a ( x ) R ( x ) = A + λ B + x I 1 I 0 f ( x ) + λ B x I 1 R ( x ) + I 0 I 1 R ( x ) Φ ( w , τ ) / w w Φ ( w , τ ) .
Applying the 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 (37), we obtain equations
A + λ B + x I 1 I 0 g ( x ) = a ( x ) R ( x ) + x I 1 λ B R ( x ) ,
A + λ B + x I 1 I 0 φ ( x ) = I 0 I 1 R ( x ) .
Differentiating (14) by x, we obtain equation
A + λ B + x I 1 I 0 d R ( x ) d x + I 1 I 0 R ( x ) = 0
which coincides with (40). Therefore, its solution φ ( x ) may be presented in the form
φ ( x ) = d R ( x ) d x .
Notice that an additional condition E φ ( x ) 0 is true due to the normalization condition. Because matrix function g ( x ) is a particular solution of non-homogeneous system (39), we choose an additional equation for it in the form (34) so that g ( x ) is uniquely defined. Let us write the second equation of system (29) with the precision up to O ε 3 :
ε 2 τ E F ( 1 ) ( w , τ , ε ) + j ε w a ( x ) E F ( 1 ) ( w , τ , ε ) = j ε w λ E 1 x E 2 + j ε w x E 2 F ( 1 ) ( w , τ , ε ) + j ε w E 2 F ( 1 ) ( w , τ , ε ) + ( j ε w ) 2 2 λ E 1 x E 2 F ( 1 ) ( w , τ , ε ) + O ε 3 .
Substituting (36) here, we obtain
ε 2 Φ ( w , τ ) τ + j ε w a ( x ) Φ ( w , τ ) 1 + j ε w E f ( x ) = j ε w { λ E 1 x E 2 Φ ( w , τ ) R ( x ) + j ε w f ( x ) + j ε w x E 2 Φ ( w , τ ) R ( x ) + j ε Φ ( w , τ ) w E 2 R ( x ) } + ( j ε w ) 2 2 Φ ( w , τ ) λ E 1 x E 2 R ( x ) + O ε 3 .
Taking into account (13), combining similar terms, dividing by ε and taking the limit ε 0 , we derive
Φ ( w , τ ) τ + ( j w ) 2 a ( x ) Φ ( w , τ ) E f ( x ) = ( j w ) 2 2 Φ ( w , τ ) a ( x ) + ( j w ) 2 Φ ( w , τ ) λ E 1 x E 2 f ( x ) + x E 2 R ( x ) + Φ ( w , τ ) / w w Φ ( w , τ ) E 2 R ( x ) .
Substituting (38) here and taking into account that E R ( x ) 1 , E g ( x ) 0 , E φ ( x ) 0 and E f ( x ) 0 , we obtain the following equation:
Φ ( w , τ ) τ = w Φ ( w , τ ) w λ E 1 x E 2 φ ( x ) E 2 R ( x ) + ( j w ) 2 2 Φ ( w , τ ) a ( x ) + 2 λ E 1 x E 2 g ( x ) + x E 2 R ( x ) .
Considering (41) and differentiating (13), we obtain
a ( x ) = λ E 1 x E 2 R ( x ) x E 2 R ( x ) = λ E 1 x E 2 φ ( x ) E 2 R ( x ) .
This is the coefficient in the first term in the right part of (42). As for the coefficient in the second term, we denote as b ( x ) :
b ( x ) = a ( x ) + 2 λ E 1 x E 2 g ( x ) + x E 2 R ( x ) .
We rewrite equality (42) in the form
Φ ( w , τ ) τ = a ( x ) w Φ ( w , τ ) w + b ( x ) ( j w ) 2 2 Φ ( w , τ )
which coincides with (31). Due to E g ( x ) 0 , we can write the following:
λ E 1 x E 2 g ( x ) = λ E 1 g ( x ) + x E 1 g ( x ) = ( λ + x ) E 1 g ( x ) .
So, (43) is coincident with (32). Equality (30) is true due to substitutions (26), (28), expansion (36), the limit ε = σ 0 and E R ( x ) 1 . Thus, the theorem is proved. □
Remark 3.
Theorem 3 shows that i τ σ 1 σ x ( τ ) / 1 / σ converges to a diffusion process. This can be regarded as a central limit theorem for the process i ( t ) , where the mean is 1 σ x ( τ ) and the variance is 1 / σ .

7. Method of Asymptotic Diffusion Analysis

In this section, we consider an implementation of the method of asymptotic diffusion analysis for obtaining probability distribution of the process i ( t ) of the number of customers in the orbit under asymptotic condition σ 0 in considered retrial queue. In what follows, we use y and y ( τ ) interchangeably.
Let us make the following inverse Fourier transform in (31):
1 2 π e j w y Φ ( w , τ ) d w = P ( y , τ ) .
We obtain the equation
P ( y , τ ) τ = y a ( x ) y P ( y , τ ) + 1 2 2 y 2 b ( x ) P ( y , τ ) .
Here y ( τ ) = lim σ 0 σ i ( τ σ ) 1 σ x ( τ ) is a limit for the normalized and centered process σ i τ σ 1 σ x ( τ ) ; Φ ( w , τ ) is its characteristic function and P ( y , τ ) is its probability density function (p.d.f.).
Equation (44) is the Fokker – Planck equation for p.d.f. P ( y , τ ) , therefore, the process y ( τ ) is a diffusion process with a drift coefficient equal to a ( x ) y and a diffusion coefficient equal to b ( x ) .
Diffusion process y ( τ ) is a solution of the stochastic differential equation
d y ( τ ) = a ( x ) y d τ + b ( x ) d w ( τ ) ,
where w ( τ ) is the Wiener process. This equation is difficult to solve directly.
We rewrite the ordinary differential Equation (13) in the form:
d x ( τ ) = a ( x ) d τ .
Consider the following stochastic process:
z ( τ ) = x ( τ ) + ε y ( τ ) ,
where ε = σ as in the previous section. It is easy to confirm that the process z ( τ ) is the limit of the normalized process σ i τ σ under condition σ 0 . Thus, z ( τ ) has a more direct relation with i ( t ) which is the quantity of interest.
We derive a stochastic differential equation for z ( τ ) .
Differentiating (47) and taking into account (45), (46), we obtain
d z ( τ ) = d x ( τ ) + ε d y ( τ ) = a ( x ) + ε y a ( x ) d τ + ε b ( x ) d w ( τ ) .
Let us rewrite coefficients in this equality in the form using the first order Taylor series expansion with respect to ϵ .
a ( x ) + ε y a ( x ) = a ( x + ε y ) + O ε 2 = a ( z ) + O ε 2 , ε b ( x ) = ε b ( x + ε y ) + O ( ε ) = ε b ( z ) + O ε 2 .
So, we can rewrite the equality with precision up to O ε 2 as follows:
d z ( τ ) = a ( z ) d τ + ε b ( z ) d w ( τ ) .
Due to ε = σ , we can write the equation
d z ( τ ) = a ( z ) d τ + σ b ( z ) d w ( τ )
which is a stochastic differential equation and its solution z ( τ ) is a diffusion process with drift coefficient a ( z ) and diffusion coefficient σ b ( z ) .
Stationary p.d.f. π ( z ) of the process z ( τ ) is a solution of the Fokker–Planck equation
a ( z ) π ( z ) + σ 2 b ( z ) π ( z ) = 0 .
We find a solution of this differential equation in the form.
a ( z ) π ( z ) + σ 2 b ( z ) π ( z ) = 0 .
The general solution for this differential equation is given in the form.
π ( z ) = C b ( z ) exp 2 σ 0 z a ( x ) b ( x ) d x ,
where C is an arbitrary constant. Because, we are looking for a solution which is a probability density function with support in [ 0 , ] , we choose C as follows.
C = 0 1 b ( z ) exp 2 σ 0 z a ( x ) b ( x ) d x d z 1 .
Therefore, we found asymptotic probability density function π ( z ) of the normalized number of customers in the orbit in the steady state, i.e., σ i ( τ σ ) as σ 0 and then τ . In the next section, we further use it to build a discrete distribution of the number of customers in the orbit.

8. Approximations of Steady-State Distributions

The goal of the work is to find enough precise approximations for discrete probability distribution P ( n 1 , n 2 ) of the number of servers working in the first and in the second phases and for probability distribution P ( i ) of the number of customers in the orbit in steady-state regime of considered multi-server retrial queue with setup time.
Two-dimensional probability distribution P ( n 1 , n 2 ) can be approximated by entries R ( n 1 , n 2 , x ) of the top-left triangle matrix R ( x ) if we choose appropriate value x of function x ( τ ) which is determined by the differential equation x ( τ ) = a ( x ) . To do this, we choose x = κ where κ is a positive root of the equation
a ( x ) = 0 .
For such value of x we have x ( τ ) = 0 that means the steady-state regime.
Algorithm for constructing of the approximations R ( n 1 , n 2 ) and R ( n ) for R ( n 1 , n 2 , x ) and R ( n , x ) is as follows. Here R ( n ) and R ( n , x ) are the probabilities that the number of busy servers (in both phases) is n.
1. Let us rewrite Equation (14) for entries R ( n 1 , n 2 , x ) of the matrix R ( x ) :
( λ + x ) R ( 0 , 0 , x ) + μ 2 R ( 0 , 1 , x ) = 0 for   n 1 + n 2 = 0 ,
( λ + x + n 1 μ 1 + n 2 μ 2 ) R ( n 1 , n 2 , x ) + ( λ + x ) R ( n 1 1 , n 2 , x ) + ( n 1 + 1 ) μ 1 R ( n 1 + 1 , n 2 1 , x ) + n 2 μ 2 R ( n 1 + 1 , n 2 , x ) = 0 for   1 n 1 + n 2 N ,
( n 1 μ 1 + n 2 μ 2 ) R ( n 1 , n 2 , x ) + ( λ + x ) R ( n 1 1 , n 2 , x ) + ( n 1 + 1 ) μ 1 R ( n 1 + 1 , n 2 1 , x ) = 0 for   n 1 + n 2 = N
and let us find a solution R ( n 1 , n 2 , x ) of the system which satisfies the normalization requirement n 1 + n 2 N R ( n 1 , n 2 , x ) 1 for given N , λ , μ 1 and μ 2 .
2. From equality (13), we derive
a ( x ) = λ n 1 + n 2 = N R ( n 1 , n 2 , x ) x n 1 + n 2 N 1 R ( n 1 , n 2 , x ) .
Using notation
r ( x ) = n 1 + n 2 = N R ( n 1 , n 2 , x )
and applying the normalization requirement
1 = n 1 + n 2 N R ( n 1 , n 2 , x ) = n 1 + n 2 = N R ( n 1 , n 2 , x ) + n 1 + n 2 N 1 R ( n 1 , n 2 , x ) ,
we can rewrite (49) in the form
a ( x ) = ( λ + x ) r ( x ) x .
3. Using numerical methods, we find a solution x = κ of the equation
a ( x ) = 0 .
4. Substituting x = κ into R ( n 1 , n 2 , x ) , we obtain approximation R ( n 1 , n 2 ) = R ( n 1 , n 2 , κ ) of two-dimensional probability distribution of the number of servers working in the first and in the second phases of service.
5. Using expression
R ( n ) = n 1 = 0 n R ( n 1 , n n 1 ) ,
we find approximation R ( n ) of the probability distribution of the number of busy servers in the considered retrial queue.
To find an approximation for discrete steady-state probability distribution P ( i ) of the number of customers in the orbit, we apply p.d.f. π ( z ) of continuous limit process normalized by σ . To do this, we apply expression (48) and find values of the following function:
G ( i ) = 1 b ( σ i ) exp 2 σ 0 σ i a ( x ) b ( x ) d x .
We obtain an approximation P A ( i ) for the distribution P ( i ) in the form
P A ( i ) = G ( i ) / i = 0 G ( i ) .
It should be noted that this approximation is not unique and we might choose other alternative distribution also [26].
Algorithm of constructing of the approximation P A ( i ) is as follows.
1. Let us rewrite Equation (33) for each non-zero entry g ( n 1 , n 2 , x ) of the matrix g ( x ) :
( λ + x ) g ( 0 , 0 , x ) + μ 2 g ( 0 , 1 , x ) = a ( x ) R ( 0 , 0 , x ) for   n 1 + n 2 = 0 ,
( λ + x + n 1 μ 1 + n 2 μ 2 ) g ( n 1 , n 2 , x ) + ( λ + x ) g ( n 1 1 , n 2 , x ) + ( n 1 + 1 ) μ 1 g ( n 1 + 1 , n 2 1 , x ) + n 2 μ 2 g ( n 1 + 1 , n 2 , x ) = a ( x ) R ( n 1 , n 2 , x ) + x R ( n 1 1 , n 2 , x ) for   1 n 1 + n 2 N 1 ,
( n 1 μ 1 + n 2 μ 2 ) g ( n 1 , n 2 , x ) + ( λ + x ) g ( n 1 1 , n 2 , x ) + ( n 1 + 1 ) μ 1 g ( n 1 + 1 , n 2 1 , x ) = a ( x ) λ R ( n 1 , n 2 , x ) + x R ( n 1 1 , n 2 , x ) for   n 1 + n 2 = N .
Then we find solution g ( n 1 , n 2 , x ) of this system which satisfies additional requirement n 1 + n 2 N g ( n 1 , n 2 , x ) 0 for given N , λ , μ 1 and μ 2 .
2. From (32), we can write
b ( x ) = a ( x ) + 2 ( λ + x ) n 1 + n 2 = N g ( n 1 , n 2 , x ) + x n 1 + n 2 N 1 R ( n 1 , n 2 , x ) .
Here n 1 + n 2 N 1 R ( n 1 , n 2 , x ) = 1 n 1 + n 2 = N R ( n 1 , n 2 , x ) = 1 r ( x ) , therefore, we can write function b ( x ) in the form
b ( x ) = a ( x ) + 2 ( λ + x ) n 1 + n 2 = N g ( n 1 , n 2 , x ) + x ( 1 r ( x ) ) .
3. We substitute obtained functions a ( x ) and b ( x ) into (52) to calculate values G ( i ) .
4. Applying (53), we obtain probability distribution P A ( i ) which approximates the goal distribution P ( i ) for small values of σ .

9. Applicability Area of Obtained Results

To estimate precision and applicability area of obtained approximations, we compare the results with empiric probabilities obtained on the base of simulation results. For numerical estimation of the error, we traditionally use Kolmogorov distance
Δ = max 0 i < k = 0 i P 1 ( k ) P 2 ( k ) ,
where P 1 ( k ) and P 2 ( k ) are discrete distributions that should be compared. In our case, one of them is approximation (53) and another one is an empiric distribution of the number of customers in the orbit obtained from simulation results.
We use
Δ 0.05
as a precision criterion for the approximation to be applicable. To avoid simulation error influence, we choose such sample sizes in the simulation experiments to empiric distributions have error Δ 0.001 (to estimate it, we use comparisons between several simulation results). You may choose any level of the error in the form of the Kolmogorov distance. Due to the definition of the Kolmogorov distance, the value 0.05 may be interpreted as a 5% error and this fact has no other meaning.
For simulation experiments, the following parameters of the retrial queue with setup time were chosen: the number of servers N = 5 , setup time distribution parameter μ 1 = 2 , main service distribution parameter μ 2 = 1 . Intensity of arrivals λ were chosen to satisfy steady-state existence requirement (19). Therefore, if we denote the system load parameter by ρ [ 0 ; 1 ) , then we can choose value of λ as follows:
λ = ρ N μ 2 ( N 1 ) μ 2 + μ 1 N μ 2 + μ 1 .
This can be rewritten as
ρ = λ N μ 2 N μ 2 + μ 1 ( N 1 ) μ 2 + μ 1 ,
representing the traffic intensity of our system.
We will vary parameters ρ and σ such that (55) is satisfied. Results of such analysis for various values of parameters ρ and σ are presented in Table 1. We highlight by boldface font values that satisfy condition (55). Therefore, we can conclude that obtained approximation (53) becomes more precise when system load parameter ρ grows and delay in the orbit increases ( σ becomes less). For values ρ 0.9 and for values σ 0.2 it becomes applicable in a wide range.
Distributions of the number of customers in the orbit obtained by simulation experiments and by the approximation for ρ = 0.8 are shown in Figure 2. As one can see, these distributions are very close to each other and especially for the cases ( σ = 1 , 0.5 , 0.2 ) where Δ 0.05 , and this fact proves the high accuracy of the obtained approximation. Just for large values of σ , we have a jump at the starting point n = 0 in real distribution which we cannot approximate properly.
An important observation is that the approximation becomes more accurate as the traffic intensity increases. This is because, we approximate the scaled number of customers in the orbit by a truncated distribution of the one (stationary distribution of the diffusion process) which might have negative support. When the traffic intensity, the mass for negative part is small and thus this approximation is appropriate.
Results on errors in mean and variance for these distributions are shown in Table 2.
In Figure 3, we compare the mean power consumption of our model in the limitting regime σ 0 (denoted by Asymptotic) with that of the buffered model [6], i.e., instead of retrial, blocked customers join the buffer in front of the servers ( σ ). Figure 3 shows the total power consumption by setup servers (phase 1) and active servers (phase 2) against the arrival rate λ . We assume that both the setup server and the active one consume an energy unit per unit time. Thus, the power consumption is equal to the mean of the sum of the number of setup servers and that of active ones. The mean number of setup servers and that of active ones in the limit regime is calculated using the distribution R ( n 1 , n 2 ) and is given by n 1 + n 2 = 1 N R ( n 1 , n 2 ) × ( n 1 + n 2 ) . For a reference, we also plot the power consumption of the ON-IDLE model (M/M/c queue), where the mean number of active servers and that of idle servers are given by λ / μ 2 and ( c λ / μ 2 ) , respectively. Here we assume that power concumption of an idle server is of 60% that of an active one. As a result, power consmption for hte ON-IDLE model is simply given by λ / μ 2 + 0.6 × ( c λ / μ 2 ) . Fixed parameters are given by μ 2 = 1 , c = 5 for both models. We observe that the power consumption of the asymptotic regime is larger than that of the buffered model. This is intuitive because in the buffered model, once a server completes a service (phase 2), a waiting customer can immediately occupy the server without setup (phase 1). Furthermore, ON-OFF models are more power-saving than the ON-IDLE model when the arrival rate is relatively small, and the opposite trend is observed when the arrival rate is large enough.

10. Conclusions

In the paper, we considered a multiserver retrial queue with setup time. Using methods of asymptotic analysis under the condition of long delay in the orbit, we derived necessary condition for the existence of the steady-state regime. Using the method of asymptotic diffusion analysis, we obtained diffusion process whose probability density function is used an approximation for the probability distribution of the number of customers in the orbit. Using simulations and numerical experiments, we showed that the approximation is good for a wide enough range of parameters. Further studies may be made for more complex models with setup time, e.g., models with impatient customers and/or with outgoing calls. The asymptotic regime where the number of servers tends to infinity [18,27] might also be worth for investigation. In our paper, although we proved that the scaled version of the number of customers in the orbit converges to a diffusion process which has a stationary distribution. In the future work, we plan to give a rigorous proof of the convergence of the stationary scaled number of customers in the orbit to the stationary distribution of the diffusion process.

Author Contributions

Conceptualization, A.N., T.P.-D.; formal analysis, A.M., T.P.-D.; investigation, A.M., A.N., S.P., T.P.-D.; methodology, A.N.; software, A.M.; writing–original draft, S.P.; writing–review & editing, A.M., T.P.-D. All authors have read and agreed to the published version of the manuscript.

Funding

The third author T.P. was supported in part by JSPS Kakenhi No. 18K18006.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ren, Y.; Phung-Duc, T.; Chen, J.C.; Yu, Z.W. Dynamic auto scaling algorithm (dasa) for 5g mobile networks. In Proceedings of the 2016 IEEE Global Communications Conference (GLOBECOM), Washington, DC, USA, 4–8 December 2016; pp. 1–6. [Google Scholar]
  2. Phung-Duc, T.; Ren, Y.; Chen, J.C.; Yu, Z.W. Design and analysis of deadline and budget constrained autoscaling (DBCA) algorithm for 5G mobile networks. In Proceedings of the 2016 IEEE international conference on cloud computing technology and science (CloudCom), Luxembourg, 12–15 December 2016; pp. 94–101. [Google Scholar]
  3. Gandhi, A.; Harchol-Balter, M.; Adan, I. Server farms with setup costs. Perform. Eval. 2010, 67, 1123–1138. [Google Scholar] [CrossRef]
  4. Gandhi, A.; Doroudi, S.; Harchol-Balter, M.; Scheller-Wolf, A. Exact analysis of the M/M/k/setup class of Markov chains via recursive renewal reward. ACM Sigmetrics Perform. Eval. Rev. 2013, 41, 153–166. [Google Scholar] [CrossRef]
  5. Maccio, V.J.; Down, D.G. Structural properties and exact analysis of energy-aware multiserver queueing systems with setup times. Perform. Eval. 2018, 121, 48–66. [Google Scholar] [CrossRef]
  6. Phung-Duc, T. Exact solutions for M/M/c/setup queues. Telecommun. Syst. 2017, 64, 309–324. [Google Scholar] [CrossRef] [Green Version]
  7. Maccio, V.J.; Down, D.G. On optimal control for energy-aware queueing systems. In Proceedings of the 2015 27th International Teletraffic Congress, Ghent, Belgium, 8–10 September 2015; pp. 98–106. [Google Scholar]
  8. Phung-Duc, T.; Kawanishi, K. Multiserver retrial queue with setup time and its application to data centers. J. Ind. Manag. Optim. 2019, 15, 15–35. [Google Scholar] [CrossRef] [Green Version]
  9. Phung-Duc, T. Single server retrial queues with setup time. J. Ind. Manag. Optim. 2017, 13, 1329–1345. [Google Scholar] [CrossRef]
  10. Phung-Duc, T.; Masuyama, H.; Kasahara, S.; Takahashi, Y. M/M/3/3 and M/M/4/4 retrial queues. J. Ind. Manag. Optim. 2009, 5, 431–451. [Google Scholar] [CrossRef]
  11. Phung-Duc, T.; Masuyama, H.; Kasahara, S.; Takahashi, Y. State-dependent M/M/c/c+ r retrial queues with Bernoulli abandonment. J. Ind. Manag. Optim. 2010, 6, 517–540. [Google Scholar] [CrossRef]
  12. Cohen, J.W. Basic problems of telephone traffic theory and the influence of repeated calls. Philips Telecommun. Rev. 1957, 18, 49–100. [Google Scholar]
  13. Falin, G.I.; Templeton, J.G.C. Retrial Queues; Chapman & Hall: London, UK, 1997. [Google Scholar]
  14. Danilyuk, E.Y.; Moiseeva, S.P.; Sztrik, J. Asymptotic Analysis of Retrial Queueing System M/M/1 with Impatient Customers, Collisions and Unreliable Server. J. Sib. Fed. Univ. Math. Phys. 2020, 13, 218–230. [Google Scholar] [CrossRef]
  15. Fedorova, E.; Nazarov, A.; Moiseev, A. Asymptotic Analysis Methods for Multi-Server Retrial Queueing Systems. In Applied Probability and Stochastic Processes; Springer: Singapore, 2020; pp. 159–177. [Google Scholar]
  16. Moiseev, A.; Nazarov, N.; Paul, S. Asymptotic Diffusion Analysis of Multi-Server Retrial Queue with Hyper-Exponential Service. Mathematics 2020, 8, 531. [Google Scholar] [CrossRef] [Green Version]
  17. Artalejo, J.R.; Gomez-Corral, A. Retrial Queueing Systems; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
  18. Halfin, S.; Whitt, W. Heavy-traffic limits for queues with many exponential servers. Oper. Res. 1981, 29, 567–588. [Google Scholar] [CrossRef]
  19. Mandelbaum, A.; Massey, W.A.; Reiman, M.I. Strong approximations for Markovian service networks. Queueing Syst. 1998, 30, 149–201. [Google Scholar] [CrossRef]
  20. Robert, P. Stochastic Networks and Queues; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013; Volume 52. [Google Scholar]
  21. Whitt, W. On the heavy-traffic limit theorem for GI/G/∞ queues. Adv. Appl. Probab. 1982, 14, 171–190. [Google Scholar] [CrossRef]
  22. Whitt, W. Stochastic-Process Limits: An Introduction to Stochastic-Process Limits and Their Application to Queues; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2002. [Google Scholar]
  23. Whitt, W. Efficiency-driven heavy-traffic approximations for many-server queues with abandonments. Manag. Sci. 2004, 50, 1449–1461. [Google Scholar] [CrossRef] [Green Version]
  24. Pender, J.; Phung-Duc, T. A law of large numbers for M/M/c/delayoff-setup queues with nonstationary arrivals. In Analytical and Stochastic Modelling Techniques and Applications; LNCS 9845; Springer: Cham, Switzerland, 2016; pp. 253–268. [Google Scholar]
  25. 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. In International Conference on Distributed Computer and Communication Networks; LNCS 11965; Springer: Cham, Switzerland, 2019; pp. 207–222. [Google Scholar]
  26. Jantschi, L. A Test Detecting the Outliers for Continuous Distributions Based on the Cumulative Distribution Function of the Data Being Tested. Symmetry 2019, 11, 835. [Google Scholar] [CrossRef] [Green Version]
  27. Maccio, V.J.; Down, D.G. Asymptotic performance of an energy-aware G/G/C queue with general setup times. Queueing Model. Service Manag. 2020, 3, 111–135. [Google Scholar]
Figure 1. Transition for states ( n 1 , n 2 ) .
Figure 1. Transition for states ( n 1 , n 2 ) .
Mathematics 08 02232 g001
Figure 2. Comparisons of the approximation (dashed line) and the simulation results (solid line) obtained for the probability distribution of the number of customers in the orbit for ρ = 0.8 and various values of σ .
Figure 2. Comparisons of the approximation (dashed line) and the simulation results (solid line) obtained for the probability distribution of the number of customers in the orbit for ρ = 0.8 and various values of σ .
Mathematics 08 02232 g002
Figure 3. Power consumption against the arrival rate λ .
Figure 3. Power consumption against the arrival rate λ .
Mathematics 08 02232 g003
Table 1. Values of Kolmogorov distance ( Δ ) of approximation distribution P A given by (53) and its exact one (by simulation) for various values of σ (retrial rate) and ρ (traffic intensity). The cases with Δ 0.05 are in boldface.
Table 1. Values of Kolmogorov distance ( Δ ) of approximation distribution P A given by (53) and its exact one (by simulation) for various values of σ (retrial rate) and ρ (traffic intensity). The cases with Δ 0.05 are in boldface.
σ 210.50.20.1
ρ = 0.6  0.0578 0.0962 0.0907 0.04030.0143 
ρ = 0.7  0.0842 0.0810 0.0471 0.0122 0.0070 
ρ = 0.8  0.0650 0.0396 0.0154 0.0035 0.0027 
ρ = 0.9  0.0247 0.0103 0.0020 0.0012 0.0010 
Table 2. Errors in mean and variance of approximation distribution P A given by (53) for ρ = 0.8 and various values of σ .
Table 2. Errors in mean and variance of approximation distribution P A given by (53) for ρ = 0.8 and various values of σ .
σ 210.50.2
Mean error8% 5% 1.5% 0.002%
Variance error0.5% 0.3% 0.04% <0.0001%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Nazarov, A.; Moiseev, A.; Phung-Duc, T.; Paul, S. Diffusion Limit of Multi-Server Retrial Queue with Setup Time. Mathematics 2020, 8, 2232. https://doi.org/10.3390/math8122232

AMA Style

Nazarov A, Moiseev A, Phung-Duc T, Paul S. Diffusion Limit of Multi-Server Retrial Queue with Setup Time. Mathematics. 2020; 8(12):2232. https://doi.org/10.3390/math8122232

Chicago/Turabian Style

Nazarov, Anatoly, Alexander Moiseev, Tuan Phung-Duc, and Svetlana Paul. 2020. "Diffusion Limit of Multi-Server Retrial Queue with Setup Time" Mathematics 8, no. 12: 2232. https://doi.org/10.3390/math8122232

APA Style

Nazarov, A., Moiseev, A., Phung-Duc, T., & Paul, S. (2020). Diffusion Limit of Multi-Server Retrial Queue with Setup Time. Mathematics, 8(12), 2232. https://doi.org/10.3390/math8122232

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