Multiserver Queue with Guard Channel for Priority and Retrial Customers

This paper considers a retrial queueingmodel where a group of guard channels is reserved for priority and retrial customers. Priority and normal customers arrive at the system according to two distinct Poisson processes. Priority customers are accepted if there is an idle channel upon arrival while normal customers are accepted if and only if the number of idle channels is larger than the number of guard channels. Blocked customers (priority or normal) join a virtual orbit and repeat their attempts in a later time. Customers from the orbit (retrial customers) are accepted if there is an idle channel available upon arrival. We formulate the queueing system using a level dependent quasi-birth-and-death (QBD) process. We obtain a Taylor series expansion for the nonzero elements of the rate matrices of the level dependent QBD process. Using the expansion results, we obtain an asymptotic upper bound for the joint stationary distribution of the number of busy channels and that of customers in the orbit. Furthermore, we develop an efficient numerical algorithm to calculate the joint stationary distribution.


Introduction
In this paper, we consider multiserver retrial queues with guard channels for priority and retrial customers.Retrial queues are characterized by the fact that a blocked customer repeats its request after a random time.During retrial intervals, the customer is said to be in the orbit.This type of queueing models is widely used in modelling and performance analysis of communication and service systems, especially in cellular networks [1][2][3][4][5][6].For instance, Tran-Gia and Mandjes [6] report the influence of retrials on the performance of cellular networks using retrial queueing models.Marsan et al. [3] carry out a fixed point approximation analysis for retrial queueing models arising from cellular networks while some extension is presented in [4].Artalejo and Lopez-Herrero analyze a multiserver queue for cellular networks operating under a random environment using a four-dimensional Markov chain.
The guard channel concept has been extensively used in communication systems [1][2][3]6].This is also referred to as trunk (or circuit or bandwidth) reservation in teletraffic literature [7,8].Tran-Gia and Mandjes [6] propose some multiserver retrial queues with fresh and handover calls and guard channels for a base station in cellular networks.In [6], the orbit size is assumed to be finite which simplifies the analysis.
The analysis of multiserver retrial queues with infinite orbit size is challenging due to the fact that the underlying Markov chain is state dependent because the retrial rate is proportional to the number of customers in the orbit.Thus, even for the fundamental model with one type of traffic and without guard channels, the stationary distribution is expressed in terms of simple functions for only some special cases, that is, one or two servers [9].We refer to [9][10][11] for some efforts in finding analytical expressions for the joint stationary distribution for the cases of more than two servers.For models with both retrial and guard channels, although some numerical methods [1,3,6,[12][13][14] have been presented for various models, there is no analytical result available.This motivates us to consider a new model with both retrials and guard channels for which we explore new analytical and numerical results.From the modelling point of 2 International Journal of Stochastic Analysis view, the novelty is the priority given to retrial customers.It should be noted that retrial customers are treated the same as normal customers in [1-3, 5, 6] which is suitable for the context of cellular networks since the base station might not be able to recognize redial calls so as to give them some priority.In [15], retrial customers are given a preemptive priority over waiting customers.To the best of our knowledge, the current paper is the first to consider the priority for retrial customers in the context of queueing models with guard channels.Our model may also be fit for systems with human servers where the service differentiation among two classes of customers is needed.In such a service system, the server can easily recognize retrial customers so as to give them some priority over fresh customers who arrive at the system for the first time.We formulate the queueing system using a level dependent QBD process where the level and the phase are referred to as the number of customers in the orbit and the number of busy channels, respectively.
The stationary distribution of a level dependent QBD process can be expressed in terms of a sequence of rate matrices [16].Thus, we can characterize the stationary distribution through the sequence of rate matrices.The QBD process of our model possesses some special structure; that is, only the last two rows are nonzero allowing us to get some insights into the structure of the stationary distribution.Liu and Zhao [17] use this property to obtain upper and lower asymptotic bounds for the stationary distribution of the fundamental retrial model without guard channels.Liu et al. [18] further extend their analysis to the model with nonpersistent customers.B. Kim and J. Kim [19] and Kim et al. [20] refine the tail asymptotic results in Liu and Zhao [17] and Liu et al. [18], respectively.Phung-Duc [21] presents a perturbation analysis for a multiserver retrial queue with two types of nonpersistent customers.In [21], the author derives the Taylor series expansion formulae for the nonzero elements of the rate matrices.The difference of our model in comparison with the above work is that the last two rows of the rate matrices are nonzero in our model while for those in [17,21] only the last row is nonzero.This makes the analysis of our model more complex and challenging.
The main contribution of our paper is threefold.First, using a censoring technique and a perturbation method, we obtain the Taylor series expansion for the rate matrices in terms of the number of customers in the orbit.Our formula is general in the sense that we can obtain the expansion with arbitrary number of terms.This was not reported in Liu and Zhao [17].Second, using this result we obtain an asymptotic upper bound for the stationary distribution which is more challenging compared to [17,21] due to the denseness of the rate matrices.It should be noted that this is the first asymptotic result for multiserver retrial queue with guard channels.Third, we present an efficient method to calculate the stationary distribution of the model whose computational complexity is linear to the number of servers.An earlier version of this paper was presented in [22].
The rest of our paper is organized as follows.Section 2 presents the model and some preliminary results on the level dependent QBD formulation.Section 3 is devoted to the presentation of the Taylor series expansion for the rate matrices.
In Section 4, we show the asymptotic upper bound for the joint stationary distribution while a numerical algorithm is presented in Section 5. Section 6 provides some numerical examples and Section 7 concludes our paper and presents some future directions.

Model and Formulation
2.1.Model.In this paper, we consider a queueing model with two types of customers (types 1 and 2).There are  servers; among them  servers are assigned as guard channels.Customers of type 1 (high priority) and type 2 (low priority) arrive at the system according to two independent Poisson processes with rates  1 and  2 , respectively.Customers of type 1 can use all  servers while those of type 2 cannot use the guard channels.If there are - busy servers, the rest of  servers automatically become guard channels for customers of type 1.Furthermore, we assume that a blocked customer (both types 1 and 2) retries after some exponentially distributed time with mean 1/.Upon retrial, if there is an idle channel the retrial customer occupies it immediately; otherwise it enters the orbit again.Thus, a retrial customer has the same priority as that of a high priority one.As a result, we may expect that decreasing the number of retrials by a customer improves the quality of service (QoS).Service times for both types 1 and 2 customers are assumed to follow the same exponential distribution with mean 1/].Remark 1.In this paper, we restrict ourselves to the case of one guard channel; that is,  = 1.This is because the asymptotic analysis for case  = 1 is complex enough and is essentially different from case  = 0. Asymptotic analysis for case  > 1 may need a new technique which will be left for a future study.On the other hand, the stability condition presented in this paper can be extended to case  > 1 in a straightforward manner while the numerical algorithm can be adapted to case  > 1.
Remark 2. From a theoretical point of view, the assumption that retrial customers (both normal and priority) have the same priority significantly simplifies the analysis.This is because if retrial customers keep their initial priority, we need to distinguish two types of retrial customers for which we should have two orbits for priority and normal customers.

Lemma 3. {𝑋(𝑡)} is positive recurrent if and only if
where Proof.The proof is presented in Appendix A.
(4) Furthermore, let   = ( 0, ,  1, , . . .,  , ) ,  = ( 0 ,  1 , . ..) . ( We have where e and 0 are vectors with appropriate dimensions with all 1 elements and all zero elements, respectively.It is established in [16] that the solution of ( 6) is given by where {R () ,  ∈ N} is the minimal nonnegative solution of Furthermore,  0 is determined by 2 ) = 0,  0 (I + R (1) + R (1) R (2) Thus the problem of finding the stationary distribution is equivalent to that of obtaining the rate matrices.However, the rate matrices do not have closed forms in general leading to algorithmic approaches for numerical calculation.To this end, we present Lemmas 4 and 5.
Lemma 5 (Proposition 2.4 in [23]).{R ()  ,  ∈ Z + ,  ∈ N} is defined by the following recursive formulae: One has International Journal of Stochastic Analysis Remark 6. Lemmas 4 and 5 allow us to derive a numerical algorithm for calculating the rate matrices.They also imply that the rate matrices are matrix continued fractions.Bright and Taylor [24,25] propose a recursive algorithm for computing rate matrix R () .From Lemmas 4 and 5, we observe that the first  − 1 rows of R () are zero.In Section 5, we propose a method for calculating R () with computational complexity of () by exploiting this sparsity and Lemma 5.
It should be noted that the recursive algorithm in [24,25] has the computational complexity of ( 3 ) due to the denseness of D-matrices (see Bright and Taylor [24,25]).
It is easy to see that the first  − 1 rows of R () and R ()   are zeros.Similar to other block matrices, we also partition R ()  as It is easy to see that R ()  = O and R ()  = O.Furthermore, we denote the elements of R ()   and R ()  as follows: Remark 7. It should be noted that, in comparison with a previous version [22], some notations have been changed.In particular,  (0,)  and  (1,)    ( = 0, 1, 2, . . ., ) in [22] are replaced by  (−1,)  and  (,)  in the current paper.

Corollary 9.
For case  = 2, explicit expressions for nonzero elements of R () are given as follows: due to (17) and ( 26) with  = 2.It follows from ( 27) and ( 21) with  = 2 that Furthermore, substituting these explicit expressions into ( 20) and ( 23) and arranging the results, we obtain ] . ( Remark 10.It should be noted that the explicit results in Corollary 9 cannot be obtained for case  ≥ 3. Thus, in the next section, we present an asymptotic expansion for the rate matrices in the general case with an arbitrary value of .
Remark 11.In this section, we find the Taylor expansion for nonzero components.The basic idea in a perturbation approach is that the coefficient of ( + 1)th term is derived based on the coefficients of the lower term expansions, that is, th term expansion ( = ,  − 1, . . ., 1).In this paper, we find the coefficients of ( + 1)th term of  (−1,) − and  (,) − for  = 0, 1, . . .,  in parallel.In principle, Lemmas 12 and 13 and Theorem 14 below could be merged into one theorem; we however present these separately in order for clarity.

Asymptotic Upper Bound
In this section, we present an asymptotic upper bound for the stationary distribution.To this end, we use Lemmas 15 and 16.
Lemma 15.For a square matrix A = ( ) and a vector x = ( 1 ,  2 , . . .,   ), one has where Lemma 16 (Fact 5 in [17]).For integers (≥ 1) and â > 0, b Remark 17.In [17], only the last row of the rate matrices is nonzero.This fact allows us to evaluate the tail probability using the product of a sequence of scalars.However, since the last two rows of the rate matrices are nonzero in our model, we need to deal with the product of a sequence of matrices.Thus, in order to apply the technique given in [17], that is, Lemma 16, we need to use Lemma 15.
Proof.The proof uses Lemmas 15 and 16.We have where where Thus, for the parameters that satisfy Lemma 16, we have implying the desired result.

Numerical Algorithm
In this section, we propose a computational algorithm for the stationary distribution of our model extending that proposed by Phung-Duc et al. [26] for the fundamental M/M// retrial queues without guard channels.In Section 5.1, we show some results which are the basis for the algorithm.Section 5.2 presents the algorithms for the rate matrices and the stationary distribution.Section 5.3 proposes a simple method for determining the truncation point used in Algorithm 2 in Section 5.2 for the stationary distribution.Section 5.4 derives some performance measures such as the blocking probability for type 2 customers (low priority) and that for type 1 and retrial calls.

Efficient Computation.
Due to Lemma 4, we need to compute  inverse matrices in order to obtain R ()  .It may take a long time when the number of servers is large.Thus, instead of computing the inverse matrices, we propose a new method exploiting the fact that only the last two rows are nonzero.The computational complexity of our new method is only ().In particular, the computational complexity in all the theorems and lemmas below is ().
It should be noted that the computation of R () and R ()    is equivalent to that of their last two rows r () and r ()  ; that is, where r (,) and r (,)  ( = 0, 1) are the row vectors of  + 1 elements.
Definition 20.We define the function   as follows.Let where y −1 and y  are the second last and the last rows of Y. Furthermore, where x and y are vectors with an appropriate dimension.

Lemma 22.
For arbitrary  and , one has where {  ,   ,  = 0, 1, . . ., } and  (,) , are given as follows: International Journal of Stochastic Analysis 9 Furthermore, Proof.This lemma can be proved using the same technique as in Lemma 21.
Remark 25.Computation of r (−1,)  and r (,)  using Lemmas 21 and 22 might be numerically unstable due to overflow.Thus, we use recursive formulae in Theorem 26 to obtain a numerically stable scheme.
Remark 28.Tran-Gia and Mandjes [6] propose some models where blocked handover calls do not retry but are lost.The results in Section 5.1 are easily adapted to these models.

Computational Algorithm.
In this section, we present an algorithm for computing the rate matrices and then a procedure for the computation of the stationary distribution.Algorithm 1 shows a method for r () while Algorithm 2 computes approximation π = (π 0 , π1 , . . ., π ) to the stationary distribution, where {  ,  ∈ Z + } is an arbitrary increasing sequence and  is the truncation point given in advance.We will discuss how to choose the truncation point in Section 5.3.

Determination of Truncation Point 𝑁.
In Algorithm 2, the truncation point is given in advance and it should be large 1 and r ()  0 using Lemmas 21 and 22 and Theorem 26. while ‖r ()    − r ()  −1 ‖ ∞ >  do  fl  + 1; Compute r ()    and r ()  −1 using Lemmas 21 and 22 and Theorem 26.r() fl r ()    ; end while Algorithm 1: Computation of r () .enough such that the tail probability is sufficiently small; that is,

Input
where  is given in advance.However, since   is not explicitly obtained for general M/M// retrial queues, a direct determination of such  is difficult.In this paper, we use the explicit results for M/M/1/1 retrial queue to determine this truncation point.In particular, we consider M/M/1/1 retrial queue with arrival rate /, retrial rate , and service rate ].This queue is stable since  = /(]) < 1 due to the stability condition of our original model.
Let  , ( = 0, 1,  ∈ Z + ) denote the probability that the number of busy servers is  and the number of customers in the orbit is  in M/M/1/1 retrial queue.It is shown in [26] that where  ∈ Z + and ()  (−∞ <  < ∞,  ∈ Z + ) denotes the Pochhammer symbol.
Using this result, we set the truncation point as follows: We verify the accuracy of this choice using numerical results.

Blocking Probability.
We derive blocking probabilities as performance measures.In our model, priority (type 1) and retrial customers are blocked when all the servers are occupied while normal (type 2) customers are blocked when at least  − 1 servers are occupied.Thus, the blocking probability of normal customers is given by and the blocking probability of priority and retrial customers is given by

Numerical Results
In this section, we show some numerical examples.In particular, in Section 6.1 we confirm the effectiveness of the Taylor series expansion for the rate matrices.Section 6.2 is devoted to the numerical investigation of the asymptotic behavior for the joint stationary distribution.Section 6.3 presents the blocking probabilities for priority and normal customers against the number of channels.

Accuracy of the Taylor Series Expansion.
The rate matrix is calculated using Algorithm 1 with  = 10 −10 and   = 2  .We call the rate matrix obtained by Algorithm 1 under this setting exact result.First, we present some numerical examples to show the effectiveness of the Taylor series expansion.Tables 1 and 2 show numerical results of R () ( = 100) for  = 1 and 10, respectively.Tables 3 and 4 show numerical results of R ()   ( = 1000) for  = 1 and 10, respectively.Other parameters are given by  = 5, ] = 1, and  2 / 1 = 4 and  is calculated from traffic intensity (= /]).We obtain exact value for the rate matrices using the matrix continued fraction approach [23] with enough accuracy (relative error of the order of 10 −10 ).The one-, two-, and three-term expansions ( = 1, 2, 3) are expressed by R (,1) , R (,2) , and R (,3) , respectively.In these tables, we show the relative errors, that is, We observe that the Taylor series expansion gives a good approximation in the sense that the relative error is quite small.The relative errors for case  = 1000 are smaller than those for case  = 100.This fact agrees with the Taylor series expansion formulae.We also observe that the relative error increases with the traffic intensity.This suggests that we need more computational effort for the cases of relatively heavy load in comparison with those of relatively light load.
Furthermore, we observe that relative error in Tables 2  and 4 ( = 10) is smaller than the corresponding one in Tables 1 and 2 ( = 1), respectively.This implies that the Taylor series expansion gives good approximation for the case of a relatively large retrial rate.This is the case of interest in practice where customers are impatient.
Figure 1 represent  (−1,)  against the number of expansion terms.The parameters are given by  = 1000,  = 100,  = 1, ] = 1,  2 / 1 = 24, and  = 0.9.We observe that the Taylor series expansion converges to the exact value after about 5 terms.Interestingly, we observe that the Taylor series expansions for  (−1,)  oscillate and converge to the exact values.

Blocking Probability versus Number of Servers.
We use the following fixed parameters ] = 1,  = 0.7, and  2 / 1 = 24 while varying the number of servers from 1 to 100.Truncation point  is determined using the method in Section 5.3 with  = 10 −10 .Blocking probabilities are   and  −1 +   for high and low priority customers, respectively.Figure 3 represents the blocking probabilities of two types of customers for three values of  (0.1, 1, and 10).Obviously, for the same , the blocking probability for low priority customers is higher than that of high priority customers.
From Figure 3, we can observe a large difference between the curve for type 1 customers and the corresponding one for type 2 customers.This implies that one guard channel is enough  to guarantee the QoS of type 1 customers.Furthermore, the blocking probabilities increase with  because customers who retry in a short interval may suffer from the same congested situation.
An important observation is that all the curves are asymptotically linear when the number of servers is large.Asymptotic analysis for the case of large number of servers may be the topic of any future research.In this direction, Avram et al. [27] consider the blocking probability under slow retrials and Halfin-Whitt regime.

Concluding Remarks
In this paper, we have introduced a new queueing model with a guard channel for retrial and priority customers.The new queueing model is formulated using QBD process which possesses a sparse structure allowing an efficient numerical algorithm and the Taylor series expansion for all the nonzero elements of the rate matrices.We have also derived an asymptotic upper bound for the joint stationary distribution.Numerical results have revealed that the upper bound can be further improved.Future work includes finding the exact asymptotic formulae for the joint stationary distribution.

A. Proof of Lemma 3
We prove the sufficient condition in Lemma 3 using Proposition A.1.

C. Proof of Lemma 13
Proof.We prove Lemma 13 using mathematical induction.

D. Proof of Theorem 14
Proof.We prove Theorem 14 using mathematical induction.First, we show that Theorem 14 is true for  = 1.(E.10)

1 Figure 3 :
Figure 3: Blocking probability versus the number of servers.