A Comparative Numerical Study of the Spectral Theory Approach of Nishimura and the Roots Method Based on the Analysis of BDMMAP/G/1 Queue

This paper considers an infinite-buffer queuing system with birth-death modulated Markovian arrival process (BDMMAP) with arbitrary service time distribution. BDMMAP is an excellent representation of the arrival process, where the fractal behavior such as burstiness, correlation, and self-similarity is observed, for example, in ethernet LAN traffic systems. This model was first apprised by Nishimura (2003), and to analyze it, he proposed a twofold spectral theory approach. It seems from the investigations that Nishimura’s approach is tedious and difficult to employ for practical purposes. The objective of this paper is to analyze the same model with an alternative methodology proposed by Chaudhry et al. (2013) (to be referred to as CGGmethod). The CGGmethod appears to be rather simple, mathematically tractable, and easy to implement as compared to Nishimura’s approach. The Achilles tendon of the CGG method is the roots of the characteristic equation associated with the probability generating function (pgf) of the queue length distribution, which absolves any eigenvalue algebra and iterative analysis. Both the methods are presented in stepwise manner for easy accessibility, followed by some illustrative examples in accordance with the context.


Introduction
In telecommunication systems, designing robust and substantial network system has become a very formidable task.Extensive studies on ethernet traffic show that, correlation, burstiness, and self-similarity are very common in the arrival rates of packaged data.Not only packetised traffic exhibits self-similar or fractal characteristics, but also there are several processes like fractional Brownian motion for storage system, autoregressive integrated moving average models, Poisson Pareto burst processes, and so forth admiting self-similarity.Conventional traffic models are not applicable to these types of self-similar traffic models.To tackle these types of communicating queues, several techniques are proposed, for example, modelling a queue with a heavy-tailed service time distribution.Among the available techniques, one of the most exquisite approaches is to consider Markovian arrival process, introduced by Neuts [1] to deal with such correlated and bursty traffic systems; for example, see Lucantoni [2] and Chakravarthy [3].
It is quite evident from the existent literature that the preeminent challenge in analyzing MAP/G/1 queue lies in extracting the boundary probabilities accurately.Insofar as, if these probabilities are known, there are several methods available in the literature to determine the remaining state probabilities and queueing-time distributions; for example, see Ramaswami [4,5], Bini and Meini [6].The most classical approach to get that stationary vector is the matrix-analytic method (MAM), developed by Neuts [7].In this method, the determination of the boundary vector q depends on the evaluation of the vector g.For that we need to find the G matrix.Once G is computed then system-length probabilities and other performance measures can be obtained by exploiting it.Several algorithms have been advocated in the literature for the evaluation of G, for example, cyclic reduction algorithm by Bini et al. [8,9].But all these proposed algorithms have 2 International Journal of Stochastic Analysis some serious computational issues.The reasons behind these difficulties are multidimensional; for example, it involves truncation error as it includes several matrix iterations and approximations; convergence becomes very slow for heavy traffic conditions.The spectral theory approach is one of the most widely used alternative methods available in the literature to obtain the unknown boundary vector; see, for example, Lee et al. [10], Gail et al. [11].Although there are several published algorithms and theorems regarding evaluation of state vectors for these types of queues, its computational difficulties have been a major concern to the practitioners.To get the flavor of the difficulty arising in computation of g as well as state vectors in //1 type queues, we would like to refer to Lee et al. [12], where he says, "Even with these low-order parameter matrices, it is still a formidable task for practitioners, who have completed introductory queueing text, to understand the published theories and write their own computer code even for the simplest MAP/G/1 queues." This paper concerns a queueing system with birth-deathmodulated Markovian arrival process (BDMMAP), introduced by Nishimura [13], where the service time is arbitrarily distributed.BDMMAP is a special type of Markovian arrival process whose underlying process is a birth-death process with a finite state space.This process has several real life applications for details; see Nishimura [13].He proposes a twofold spectral approach to analyze this queue.In his paper, he identifies the zeros of |I − A()| in the unit disk (which are found to be real and distinct) using a modified bisection method and subsequently presents the G matrix in terms of the zeros and associated right null vectors of I − A().Further, to calculate the stationary distribution at the service completion epochs, he has proposed an algorithm based on Fourier inversion formula of the generating function Y + () of the stationary probability distribution.
In this paper, our objective is to compare the CGG approach with that of Nishimura to analyze the present model in terms of computational difficulties.In this regard, we would like to emphasize on the fact that if a practitioner or an engineer who has basic knowledge of spectral theory and understood the method of Nishimura will face extreme difficulty to numerically implement this.Also, in his paper [13], the state vectors are not presented numerically for further implementation of the users.These are the facts which actually motivated us to advocate an alternative approach, namely, CGG method, to analyze the present model, which seems to be easily understandable and rather easy for numerical demonstration.The main advantages of this approach over Nishimura are that one does not need to evaluate any eigenvalues and eigenvectors (especially in terms of ) for the analysis of the present model, which are the essential parts in Nishimura's approach.Here, by CGG method, we first evaluate the vector q, using roots inside the unit circle of the numerator of |I − A()| and then evaluate the rest of the state probabilities by inverting the vector generating function (VGF) using the roots outside the unit circle.Further for the validation of our objective, we present the steady state probabilities at both arbitrary epoch and postdeparture epoch in tabular form, obtained from CGG method for high order arrival matrices, by writing codes in Maple 13 software.
Also to check the accuracy of the two methods, we compare q numerically obtained from both the methods and observe that these values agree up to a desirable degree of accuracy.
This paper is organized as follows: Section 1 consists of introduction of the paper; Section 2 includes the model description; Section 3 consists of the analysis of the present model by both CGG method and Nishimura method; in Section 4, we present the numerical results; in Section 5 we provide conclusion and future scope of this paper, followed by the references.

Model Description
Although Nishimura [13] presented a detailed description about the queue //1 but for the sake of completeness, we here present a brief overview of the model.

Service Time Distribution.
Here the customers are served according to FCFS discipline and by a single server having service time arbitrarily distributed whose cdf, pdf, mean service time, and Laplace-Stieltjes transform are, respectively, (), h(), 1/, and ℎ * () = ∫ ∞ 0  − d().Traffic intensity  is given by /.

Analysis
In the literature (see Lucantoni [14]), the probability generating function (pgf) Y + () of the stationary probability vector just after service completion epochs for the //1 queue is given by ( Here Y + () = ∑ ∞ =0 y +    , where y +  = ( + ,1 , . . .,  + , ) ( ≥ 0) are the -dimensional row vector whose th entry is the stationary joint probability that just after service completion epochs, the arrival phase is  and the number of customers in the system is .
The main steps involved in the analysis of //1 type queue are as follows: (i) one needs to evaluate the vector g (q = (1 − )g) and (ii) after getting g, one needs to invert the generating function Y + () to get the stationary probability vector.In the following sections, we briefly discuss the steps required for the analysis of //1 queue by Nishimura's method and by CGG method, respectively.Nishimura.We here briefly discuss the procedure adopted by Nishimura.For clarity of understanding we divide his approach into two subsequent segments; namely, first one is the evaluation of g and the second one is the evaluation of steady state probability vector, which is as follows.

The Evaluation of g.
In the following, we describe the procedure adopted by Nishimura for the evaluation of g.
(i) Nishimura shows that, for a BDMMAP/G/1 queue, the   's ( = 1, . . ., ) (zeros in the unit disk || ≤ 1 of |I−A()|) are distinct and simple and G can be exhibited as where ^ are the right null vectors of I − A() corresponding to   .
(ii) ^'s are calculated from the following equations: Accordingly,   and ^ are the eigenvalues and right eigenvectors of the G matrix.A() is the matrix probability generating function of arrival numbers, that is, of A  which is an  ×  matrix and whose (, )-entry is the transition probability under the condition that the arrival phase is  at the service starting epoch, the arrival phase at the service completion epoch is , and the number of arriving customers during the service time is .
(iii) Now g is calculated from gG = g and ge = 1.

The Evaluation of Steady State Probability
Vector.The procedure adopted by Nishimura for evaluating the steady state probability vectors is as follows.(ii) With the help of the above facts and assumptions, he has represented the generating function in the form ) U (), International Journal of Stochastic Analysis (b) where and  , ()'s are calculated from the following iterative procedure: (1) (2) where   () = ( −1 ()  ()) which is a symmetric tridiagonal matrix, with   () = √ −1 ()  ().
(iii) After obtaining Y + (), he uses fast Fourier transform (FFT) to get the stationary probabilities.In order to apply FFT he has to determine a number  (preferably  = 2  ) so that after that number the remaining probabilities are negligible and further he obtains the approximated state vector as y +  = (1/) ∑ −1 =0 Y + (  ) − , ( = 0, 1, . . .,  − 1).In order to calculate the eigenvalues   (  ) of (  ),  =  2/ he uses Durand-Kerner-Aberth method for calculating zeros of a polynomial.Note 1.We have tried to implement the above method by writing codes in Maple 13 software, in a PC with Intel CORE(TM) i5, 3.10 Ghz processor.Our experience and the realizations during the computation by Nishimura's method are summarized in the following remarks.
Remark 1.An addressable disadvantage in Nishimura's approach is that we need to find the eigenvalues of D() (as well as of F()) as a function of  only, using a bisection method based on Givens' algorithm, which may be cumbersome in computation for higher order arrival matrices.In this regard, we would like to quote as in Lee et al. [12] that ". ..we have experienced difficulties in finding all the eigenvalues as a function of z when the order is high.In this case, commercially available mathematical packages may help."As far as our experience, commonly available mathematical packages may not always solve the present problem up to the desired level.
Remark 2. For the evaluation of G by Nishimura's method, one has to compute the ^'s corresponding to each of the   's, which evidently needs a thorough computational involvement even for a moderate value of .However, it seems to us that computation up to g or the boundary vector q by this method is not as difficult as the evaluation of the other stationary vectors.
Remark 3. In order to derive Y + () by Nishimura's method, one needs to first find U() & V().For that one has evaluate the eigenvalues and eigenvectors of F(), which is again a tedious task.After getting these eigenvectors, he uses a recursive algorithm to evaluate U() & V(), which may address numerical inaccuracy for high order matrices.
Remark 4. In order to apply FFT on Y + (), he has to calculate the eigenvalues   (  ) of (  ),  =  2/ , using the Durand-Kerner-Aberth method.Although, Durand-Kerner-Aberth method has superlinear convergence for simple zeros but its convergence heavily depends on selecting the initial guess.Also, if the concerned polynomial has zeros with large moduli and small moduli at the same time, selecting the initial points becomes rather delicate.Also, for stopping the iteration, one needs a tedious backward rounding error analysis, which asserts that each computed approximation is the zero by a "nearby" polynomial (see Bini [15]).So, according to us, one may not find it very easy from the computation point of view.
Remark 5.Although FFT algorithm has computational complexity ( log ) instead of ( 2 ) (as in the case of DFT) (if  is of the form 2  , a one-dimensional FFT of length  requires less than 3 log 2  floating-point operations (times a proportionality constant)), it looks computationally very much challenging to the users to invert the pgf by FFT algorithm in Nishimura's case, due to the complex structure of the pgf presented by him.We have tried to invert the pgf by FFT for very low order arrival matrices and found it very difficult.So, as far as the authors' experience, this methodology may not be suitable for the practitioners with  = q 2 = q 3 = q 4 = q 5 = q 6 = q 7 = qe = 1 −  Encountering the difficulty to numerically implement the present queueing model by Nishimura's method, we apply the recently developed CGG method (see [16]) to obtain the steady state vectors.For the sake of completeness, we present the CGG approach in stepwise manner in the following section.

Evaluation of the Vector q
Step 1.We write the pgf as (a) Y() = [ 1 (),  1 (), . . .,   ()] is the vector generating function of ,  ,2 , . . .,  , ] is a row vector whose th element denotes the joint probability that there are  ( ≥ 0) customers in the system and the phase of the arriving customer is  (1 ≤  ≤ ) at arbitrary epoch.
(c) q = [ 1 ,  2 , . . .,   ] is a row vector whose th element represents the joint probability that an arriving customer is in phase  and the server is idle.
Step 2. It is seen that, in applications, most of the distributions have rational Laplace-Stieltjes transform.They exploit this fact.They consider the LST of the service time distribution as ℎ() = ()/(), where the degree of the polynomial () is  and that of () is at most .Now the (, )th element of (I − A()) is given by (I − A()) , =  , ()/(), where Step 3. From the system of ( 15 Step 4. Each component Y  () of the VGF () given in being analytic in || ≤ 1 implies that  1 ,  2 , . . .,   must be the zeros of the numerator of Y  () and hence we can determine the unknown vector q by considering any one component of () say   () (1 ≤  ≤ ).This implies that   (  ) = 0, 1 ≤  ≤  and another equation comes from normalizing condition, that is,   (1) =     (1).

Evaluation of the State Probability
Vector at Postdeparture Epoch y +  .If () has  roots, namely,   (1 ≤  ≤ ), outside the unit circle then by partial fraction method   () =  0, + ∑  =1 ( , /( −   )),  = 1, . . ., .We get where, Remark 6.In Nishimura's method, we have to compute all the eigenvalues   and the corresponding eigenvectors ^ of I − A().Also we need to evaluate the eigenvalues and eigenvectors of D() and F(), which is not easy even for a moderate value of , as mentioned earlier.CGG method has an impeccable advantage over this method as it does not involve determination of any eigenvalues and eigenvectors.
In this method, one needs to find the roots of the c.e. in Ω and Ω * .Nowadays, with improved computational aids, the determination of roots of a large degree polynomial, however, close and in any interval is not a big deal.Also, by CGG method, we can extract the vector q from any one component of the vector generating function Y() which is a notable side of this approach.
Remark 7. Nishimura [13] shows that for //1 queue the roots of the characteristic equation |I − A()| inside the unit circle are real and distinct, so we need not to bother about the multiple root case inside Ω in CGG method.This is an added advantage of using CGG method.
Remark 8.In order to apply the CGG method only we have to assume that the service time distribution has rational LST.It is well known that a wide class of distribution produces rational LST.Also, for the distribution which does not possess rational LST, we can use some approximation method which will replace it with a rational function; for example, for deterministic distribution (which has transcendental LST), we can use Padé approximation to rationalize its LST (for references see [17][18][19][20]).Also to justify our claim for rational approximation, we present the boundary vector in Table 3, for the present model in the case of deterministic service time distribution and match it to that obtained by Nishimura's method.

Numerical Results
An extensive numerical effort has been put forth and a considerable number of examples have been studied, in order to validate the objective of this paper.However, in this section, we present few of them in the form of self-explanatory tables.
For examples 1 and 2, considered here, we took  = 3.560994, service time distribution as phase type distribution with parameters  := 2,  := ( −41.2 35 25 −52.5 ),  := ( 1 0 ) with mean 0.067934.In Table 1, we obtain the distribution of stationary probability at arbitrary epoch using CGG method for  = 7.The other input parameters are displayed in the table.We present the stationary probability vector ,   's in the table.It is shown in the table that the numerical value of qe matches 1 − .Here, we would like to mention that we have compared the values of q obtained by both the Nishimura and CGG method and which are found to be the same.As we were unable to invert the pgf  + () even for the low value of  by Nishimura's method due to its complicated structure presented by him, we restrict up to the comparison of q only.The last row and column of Table 1 represent the marginal distributions of the number of the customers in the system and the distribution of the phases, respectively.In Table 2, we present the distribution of stationary probability at service completion epoch using CGG method.The input parameters in this example are the same as in Table 1.The other things are analogous to that of Table 1.
In Table 3, we present the boundary vector q for the present model in case of deterministic service time distribution and further it is observed that it matches up to desired degree of accuracy with that obtained by Nishimura's method.
Here we use Padé approximation with the parameters [5,5] to rationalize the transcendental LST which results from the deterministic distribution.The other particulars about the parameters taken for the example are presented in the table.In Table 4, we present the higher order random epoch probabilities for the deterministic service time distribution for the input parameters the same as in Table 3.Here the c.e. () = 0 admits 510 roots of which 10 roots lie inside Ω and the others lie in Ω * .For this very example we would like to emphasize on the fact that it is observed that of the 500 roots in Ω * several roots are repeated and the number of distinct roots is 100.It is often commented in queueing literature that it is quite difficult to handle the multiple roots although it may be admitted from our observation that it is not always the case.Here we find the probabilities from the inversion of () using the partial fraction technique with simple algebraic manipulation in case of roots with multiplicity more than 1.Also to check the accuracy of the probabilities computed here, we find the tail probabilities using the root with smallest absolute value (here it is  = 1.09213352) of the c.e. that lie in Ω * and these values agree up to satisfaction after  = 18 (for details see the Table 4).
The purpose of the examples presented here is to compare the accuracy of the two methods (up to the level of q) and provide some numerical demonstration of the present queueing model by obtaining steady state distributions both at arbitrary and postdeparture epoch by CGG method for various input parameters.Due to the lack of space, we have not presented the numerical results for higher values of .However, these can be obtained from the author, if needed.

Conclusion and Future Scope
This paper presents a parallel study of the Nishimura and CGG method, on the basis of the numerical implementation of the //1 queue.The main objective of this paper is to inform the readers about the two methods from the computation point of view.A number of varied examples throughout the investigation bring out the fact that implementation of Nishimura's method is very very difficult if not impossible, even for the low order input parameters, whereas the CGG method seems to tackle the present queueing model in a much better way.
Further, the CGG method may be applied to analyze a wide variety of queueing models with practical importance, for example, TSMAP/G/1 queue (which is mainly used to model tree structured LAN traces); for details, see [21].Also, it would be interesting if any recursive algorithm or some analytic computation scheme for the computation of the roots of the c.e. can be developed for the present queueing model.

Table 2 :
Distribution of stationary probability at service completion epoch using CGG method.
basic knowledge of queueing, spectral theory, and computer programming to computationally implement it.

Table 4 :
Distribution of stationary probability at arbitrary epoch using CGG method for deterministic service time distribution.