Asymptotic Analysis of Finite-Source M/GI/1 Retrial Queueing Systems with Collisions and Server Subject to Breakdowns and Repairs

This paper deals with a retrial queuing system with a finite number of sources and collision of the customers, where the server is subject to random breakdowns and repairs depending on whether it is idle or busy. A significant difference of this system from the previous ones is that the service time is assumed to follow a general distribution while the server’s lifetime and repair time is supposed to be exponentially distributed. The considered system is investigated by the method of asymptotic analysis under the condition of an unlimited growing number of sources. As a result, it is proved that the limiting probability distribution of the number of customers in the system follows a Gaussian distribution with given parameters. The Gaussian approximation and the estimations obtained by stochastic simulations of the prelimit probability distribution are compared to each other and measured by the Kolmogorov distance. Several examples are treated and figures show the accuracy and area of applicability of the proposed asymptotic method.


Introduction
Finite-source retrial queues are very useful and effective stochastic systems to model several problems arising in telephone switching systems, telecommunication networks, computer networks and computer systems, call centers, wireless communication systems, etc. To see their importance the interested reader is referred to the following works and references cited in them, for example Artalejo and Corral (2008); Dragieva and Phung-Duc (2017); Falin and Artalejo (1998); Falin and Templeton (1997); Gómez-Corral and Phung-Duc (2016) and Kim and Kim (2016).
Since in practice some components of the systems are subject to random breakdowns it is of basic importance to study reliability of retrial queues with server's breakdowns and repairs. Finite-source retrial queues with an unreliable server have been investigated in several recent papers, see for example Almási et al. (2005); Dragieva (2014); Gharbi and Dutheillet (2011); Ikhlef et al. (2016); Roszik (2004); Wang et al. (2010Wang et al. ( , 2011 and Zhang and Wang (2013).
Searching the scientific databases we have noticed that relatively just a small number of papers have been devoted to systems when the arriving calls (primary or secondary) causes collisions to the request under service and both go to the orbit, see for example Ali and Wei (2015); Choi et al. (1992); Kim (2010); Kumar et al. (2010) and Peng et al. (2014). It should be noted that collisions decrease the effectiveness of the system performance. That is why new protocols should be developed to avoid the collision but unfortunately, it cannot be neglected, see Cao et al. (2018); Jinsoo et al. (2018); Kwak et al. (2017Kwak et al. ( , 2018, Wentink (2017) and Yeo et al. (2017). This fact shows the importance of the mathematical modeling of such systems.
A research group in Tomsk State University (Russian Federation) has developed a very effective asymptotic method (Nazarov and Moiseeva 2006) by the help of which various systems have been investigated. Concerning to finite-source retrial systems with collision we should mention the following papers (Kvach 2014;Kvach and Nazarov 2015a, b, c;Nazarov et al. 2014).
A research group in Univesity of Debrecen (Hungary) has been dealing with systems with unreliable server/s as can be seen, for example in Almási et al. (2005); Sztrik (2005); Sztrik et al. (2006) and Wüchner et al. (2010) and that is why it was understandable that the two research groups started cooperation in 2017.
The present model is a generalization of the M/G/1//N retrial system treated in Kvach and Nazarov (2015a) where the server was reliable and the M/M/1//N system with unreliable server analyzed in . The generalization has its own importance since the resulting new Kolmogorov-equations are more complicated and hence the solution needs more investigations. Even the method is the same but the system, equations, and the solution are different. The situation is similar when we use generating function, Laplace-transform, characteristic function approach for the solution of a given problem. In Tóth et al. (2017) the present model has been analyzed employing stochastic simulation and in  by an algorithmic approach. The examples treated in these papers are the test results for comparisons to the present paper.
In this paper, a finite-source M/GI /1 retrial queuing system with collisions of customers and an unreliable server is considered. Applying the method of asymptotic analysis under the condition of an unlimited growing number of sources it is proved that the limiting probability distribution of the central and the normalized number of customers in the system follows a normal law with given parameters. Based on the first two asymptotic moments a Gaussian approximation is constructed. Accuracy and the range of applicability of the proposed approximation are analyzed. Based on these results in a consecutive paper by plan the asymptotic distribution of the waiting time distribution will be analyzed.
The rest of the paper is organized as follows. In Sections 2 and 3 the full description of the model is given with the help of the corresponding two-dimensional Markov process. Using a supplementary variable of residual service time the corresponding Kolmogorovequations are derived in steady-state. Sections 4 and 5 are devoted to the first and secondorder asymptotic. In Section 6 a Gaussian approximation is proposed to get the prelimit distribution of the number of customers in the system. We deal with several examples and comparisons showing the advantage of the asymptotic methods, and some comments are made. Finally, the paper ends with a Conclusion and some plans are highlighted.

Model Description and Notation
Let us consider ( Fig. 1) a retrial queuing system of type M/GI /1//N with the collision of the customers and an unreliable server. The number of sources is N and each of them can generate a primary request with the rate λ/N. A source cannot generate a new call until the end of the successful service of this customer. If a primary customer finds the server idle and not failed he enters into service immediately in which the required service time has a probability distribution function B(x). Let us denote the service rate function by μ(y) = B (y)(1 − B(y)) −1 and it's Laplace -Stieltjes transform by B * (y), respectively. Otherwise, if the server is busy an arriving (primary or repeated) customer involves in collision with the customer under service and they both move into the orbit. The retrial time of requests is exponentially distributed with the rate σ/N. We assume that the server is unreliable. The lifetime is supposed to be exponentially distributed with failure rate γ 0 if the server is idle and with the rate γ 1 if it is busy. When the server breaks down it is immediately sent for repair and the recovery time is assumed to be exponentially distributed with the rate γ 2 . We deal with the case when the server is down all sources continue the generation of customers and send it to the orbit. Similarly, customers may retry from the orbit to the server but all arriving customers immediately go into the orbit. Furthermore, in this unreliable model, we Fig. 1 Finite-source M/GI /1 retrial queueing system with collision of the customers and unreliable server suppose that the interrupted request goes to the orbit immediately and its next service is independent of the interrupted one. All the above mentioned random variables are assumed to be independent of each other.
Let Q(t) be the number of customers in the system at time t, that is, the total number of customers in the orbit and in service. Similarly, let C(t) be the server state at time t, that is 1, if the server is busy at time t, 2, if the server is down (under repair) at time t.
Thus, we will investigate the process {C(t), Q(t)} which is not Markov. Using the supplementary variable method let us introduce the random process z(t) equals to the residual service time, that is the time interval from the moment t until the end of the successful service of the customer.
As we can see {C(t), Q(t), z(t)} is multi-dimensional Markov process, the component z(t) is determined only in those moments when the server is busy, that is when C(t) = 1.

Kolmogorov-Equations for the Steady-Steady Probability Distribution of the Number of Customers in the System
The present paper is a continuation of investigations published in Nazarov et al. (2017a). However, for a better understanding and compactness of our research, we repeat some parts of the above-mentioned paper.
Let us define the following probabilities It is easy to see that P 1 (j, t) = lim In a standard way, we can write the corresponding steady-state Kolmogorov-equations as follows Denoting the partial characteristic functions by where i = √ −1 is the imaginary unit then system (1) can be rewritten as Summarizing the equations of the system (3) and letting z → ∞ we obtain (4) We aim to solve systems (3)-(4) by the method of asymptotic analysis under the limiting condition of an unlimited growing number of sources. Since for finite values N the authors do not know any method for solving systems of linear differential equations with variable coefficients we apply N → ∞. Besides, to the best knowledge of the authors, such systems do not have an analytical solution.

Asymptotic of the First Order
Theorem 1 Let Q(∞) be the number of customers in the system in steady-state then where the value of parameter κ 1 is the positive solution of the equation and the stationary distributions of probabilities R k [κ 1 ] of the service state k are defined as follows The proof of the theorem was published earlier in Nazarov et al. (2017a) therefore we give only the system of equations from the proof of Theorem 1 which we will need in the proof of Theorem 2, namely Besides, R 1 (z) denotes the steady-state probability that the server is busy and the residual service time is less than z which was obtained in the following form If we denote by R * 1 (w) its Laplace-Stieltjes transform then it is not difficult to show that Furthermore, R 0 , R 1 stand for the steady-state probability that the server is idle, busy, respectively. Here for the simplified notation, we do not show the dependence on κ 1 of the corresponding quantities, since as soon as κ 1 is obtained they are defined. We will return to this system later in the proof of Theorem 2.
Here we should mention that the solution of Eq. 6 is not unique as we could see in the paper (Nazarov et al. 2017b). In that work an example was given for the reliable case with gamma-distributed service time and very special input parameters. The bi-stability phenomenon is very seldom was treated and explained.
In the following, we assume that the solution is unique. This could be validated, for example with the help of simulation by drawing the steady-state distribution of the number of customers in the system.
Hence the prelimit value of the average number of customers in the system, that is when N is fixed can be approximated by

Asymptotic of the Second Order
The novelty of the paper is in this section since we will get a second-order approximation for a system with an unreliable server. The method is the same as in Kvach and Nazarov (2015a) but in this case, the system of Kolmogorov-equations is more complicated due to the unreliability of the server. To get the solution we need more steps and the structure of the solution is different. To get the variance is not obvious even in a special case when the failure rates of the server are zero which is the model treated in Kvach and Nazarov (2015a). To obtain a better approximation we need the variance that is the reason we state the following theorem the value of parameter κ 2 is defined by expression where and R * 1 (a + γ 1 ) is the Laplace-Stieltjes transform of R 1 (z) given by eq.(10).
system (14) can be rewritten as Our method is constructive. We show that the solution of this system can be written in the form of decomposition Substituting this form (17) into (16) keeping in mind the notion of o(ε) after simplifications we receive the following system of equations After some calculations it can be seen that as ε → 0 this system can be rewritten in the form From these systems we get that function Φ 2 (w) can be written as coinciding with (11).

From (19) it follows that the expression
is constant and has a form Thus, system (18) can be rewritten as Let us consider in more detail the second equation of the system (20). Using standard methods it is not difficult to see that the solution of this equation is Execute the limiting transition at z → ∞ and take into account that the first factor of the right part of equality (21) in a limiting condition tends to infinity. Thus we can conclude that the second factor is equal to zero, that is Performing simple calculations we obtain Taking into account the obtained equality (22) let us rewrite system (20) in the following form Solution of the system (23) can be obtained in several steps, namely Step 1. Add the first and third equations of system (23) and subtract the second equation of system (23) from their sum then we obtain From the fourth equation of system (9) it follows that R 1 (0) = λ(1 − κ 1 ). Also it is not difficult to get that a 1 − B * (a + γ 1 ) = R 1 R 0 (a + γ 1 ). Then we can rewrite (24) in the form from which we obtain Step 2. Let us consider the third equation of the system (23). Multiplying term by term on λ(1 − κ 1 ) = aR 0 B * (a + γ 1 ) which follows from the second and fourth equations of the system (9) we obtain an equation of the form Adding the received equation to the second and fourth equations of system (23) we have Thus it is not difficult to see that γ 0 + γ 2 γ 2 R 0 − 1 = − γ 1 + γ 2 γ 2 R 1 . Substituting it into equality (26) we obtain equality in the form from which it follows that Step 3. As a result, it can be seen that the left-hand side of equalities (25) and (27) coincide. Thus let us equate the right-hand side and after simple calculations, we receive the value of parameter κ 2 in the form (12).
Due to the structure and functional coefficients of the system of equations (3)-(4) and then (16) we cannot show formally that (17) is the asymptotic solution. It was a procedure using constructive steps from our previous experiments to get somehow a solution. In this sense, it is not rigorous proof. As we mentioned earlier to the best of our knowledge there is no analytical solution to this kind of system thus we cannot expect rigorous proof. All in our previous papers the goodness of the proposed asymptotic solutions were tested by either numerical or simulation results which justified our constructiveness concerning the form of a solution, see for example Nazarov et al. (2018. In this paper, we follow this method and in the next section, the comparisons are made. The theorem is proved. From the proved theorem it follows that if N → ∞ the limiting distributions for the centered and normalized number of customers in the system has a Gaussian distribution with variance κ 2 defined by the expression (12).
Owing to replacements (13) and (15) in prelimit case the average and variance of the number of customers in the system can be approximated by Nκ 1 and Nκ 2 , respectively. As a consequence the number of customers in the system follows a normal distribution with parameters Nκ 1 and Nκ 2 .
It is not difficult to see, that in the case when γ 0 = γ 1 = 0, that is with a reliable server our results coincide with the ones published in Kvach and Nazarov (2015c). Furthermore, after elementary calculations we can obtain the formulas for the case when the service time is exponentially distributed with intensity μ, that is for the M/M/1//N system with unreliable server analyzed in .

Examples and Comparisons
In this section we would like to show the effectiveness of the proposed approximation comparing the results to the estimations obtained by numerical and to the simulation calculations, see Tóth et al. (2017), respectively.

Gaussian Approximation
Previously it was obtained that the limiting distribution for the centered and normalized number of customers in the system has a Gaussian distribution with mean Nκ 1 and variance Nκ 2 .
We propose to approximate the prelimit discrete probability distribution Π(j) of the number of customers in the system obtained with the help of either numerical algorithm or simulation, by a normal distribution with the above parameters.
Approximation of discrete distribution by a continuous Gaussian distribution can be executed in various ways. We will use the following form.
Let us denote by G(x) the distribution function of the Gaussian distribution with mean Nκ 1 and variance Nκ 2 .
Furthermore, let us denote by P as (i) the asymptotic discrete distribution obtained with the help of Gaussian approximation, that is which is called the Gaussian approximation of the prelimit distribution Π(i).

Numerical Results and Comparative Analysis
For the considered retrial queuing system we choose gamma-distributed service time S with shape parameter α and scale parameter β with Laplace-Stieltjes transform B * (δ) of the form in the case when α = β, that is when the average service time is equal to unit. It can be shown that where V 2 S denotes the squared coefficient of variation of S. This distribution allows us to show the effect of the distribution on the main performance measures because dealing with the same mean we can see the impact of the variance, too.
To compare the prelimit probabilities Π(j) of the number of customers in the system and its Gaussian approximation P as (j ) at various N we give some graphical illustrations. As an example, the following parameters were chosen λ = 1 , σ = 5 , γ 0 = 0.1 , γ 1 = 0.2 , γ 2 = 1 , α = β = 1.778 . Figures 2, 3, 4a demonstrate the difference between the prelimit distribution Π(j) (solid line) and its Gaussian approximation P as (j ) (dashed line) at various values of N . As we can see for the case N = 10 the difference is noticeable. With increasing N the difference begins to decrease for the case N = 30 it is still distinguishable but for N = 50 it is practically invisible. From the presented Figures it is possible to draw an obvious conclusion that with an increase of N the effectiveness of the proposed approximation improves noticeably.
In Figs. 2-4 (b) the difference between the distribution function F (j) (solid line) and its normal approximation F as (j ) (dashed line) is presented as the values of N increases. On each Figure we can see that the curves are very close to each other and already at N = 30 the difference is practically not noticeable.   We think that visual comparison is not enough to determine the range of applicability of the approximation. To determine the accuracy and area of applicability of the Gaussian approximation we will use the Kolmogorov-distance which can be defined as follows Δ = max 0≤j ≤N j n=0 (Π(n) − P as (n)) . Table 1 represents the values of Kolmogorov-distance for the same set of parameters as for graphic illustrations, except for N, α.
We choose the criterion for the applicability of the approximation as follows: the acceptable error must not exceed the value 0.05, that is Δ ≤ 0.05.
From Table 1 we can see that practically in the whole spectrum of values of the variables α = β and N the Kolmogorov-distance does not exceed the set value 0.05 and therefore in this situation the application of Gaussian approximation is preferable.

Conclusion
In this paper a closed retrial queuing system of type M/GI /1 with collisions of customers and an unreliable server was considered. Applying the method of asymptotic analysis under the condition of an unlimited growing number of sources it was proved that the limiting probability distribution of the central and the normalized number of customers in the system follows a normal law with given parameters. Accuracy and the range of applicability of the proposed approximation were treated. In the future, for the considered system we plan to investigate the number of transitions of the customer into the orbit, the sojourn time of the customer in the orbit/system, and other system performance measures. article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.