Diffusion Limit of Multi-Server Retrial Queue with Setup Time

: 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 ﬁnd 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 inﬁnite 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 inﬁnite buffer in front of the servers.


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.

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/σ → ∞).
Summing up all equations (1)-(3), we derive 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 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.

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: Then we obtain the following system: We can prove the following statement.

Theorem 1.
Under the asymptotic condition σ → 0, the following equality holds.
Here the scalar function x(τ) is a solution of ordinary differential equation where R = R(x) is a left-top triangle matrix which is a solution of homogeneous linear system and satisfies the normalization condition Proof. Denoting lim ε→0 F(w, τ, ε) = F(w, τ) and making asymptotic transition ε → 0 in (11), we derive We find the solution of this system in the form 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: 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.

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
Proof. Sufficiency of condition (19) was proved in [8]. Let us prove its necessity. We rewrite (18) in the form: 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 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 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: 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: Solution of this system is as follows: 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): 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: Taking here the limit x → ∞, we obtain Due to (23), we can rewrite it as follows: 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 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.

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) We obtain 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 we derive We can prove the following statement.

Theorem 3.
Let Φ(w, τ) be an asymptotic (as σ → 0) characteristic function of normalized and centered process Then it satisfies the equation where function a(x) is determined by (18) and top-left triangle matrix g(x) is a particular solution of the system of equations and satisfies additional condition Proof. We rewrite the first equation of (29) with precision up to infinitesimals of order ε 2 : and let us write its solution in the form of expansion: where f(x) is some matrix function whose expression we will obtain later. Substituting (36) into (35), we obtain 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): Applying the superposition principle, we can write a solution of this equation in the form Substituting it into (37), we obtain equations Differentiating (14) by x, we obtain equation which coincides with (40). Therefore, its solution ϕ ϕ ϕ(x) may be presented in the form 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 : Substituting (36) here, we obtain Taking into account (13), combining similar terms, dividing by ε and taking the limit ε → 0, we derive ∂Φ(w, τ) ∂τ Substituting (38) here and taking into account that ER(x) ≡ 1, Eg(x) ≡ 0, Eϕ ϕ ϕ(x) ≡ 0 and Ef(x) ≡ 0, we obtain the following equation: Considering (41) and differentiating (13), we obtain 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): We rewrite equality (42) in the form which coincides with (31). Due to Eg(x) ≡ 0, we can write the following: So, (43) is coincident with (32). Equality (30) is true due to substitutions (26), (28), expansion (36), the limit ε = √ σ → 0 and ER(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

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): We obtain the equation Here 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 where w(τ) is the Wiener process. This equation is difficult to solve directly. We rewrite the ordinary differential Equation (13) in the form: Consider the following stochastic process: 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 Let us rewrite coefficients in this equality in the form using the first order Taylor series expansion with respect to .
So, we can rewrite the equality with precision up to O ε 2 as follows: Due to ε = √ σ, we can write the equation 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 We find a solution of this differential equation in the form.
− [a(z)π(z)] The general solution for this differential equation is given in the form.
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.
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.

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.
5. Using expression 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: We obtain an approximation P A (i) for the distribution P(i) in the form 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.

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 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 (55) 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: This can be rewritten as 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 n 1 +n 2 =1 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.

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.

Conflicts of Interest:
The authors declare no conflict of interest.