Table of Contents
International Journal of Stochastic Analysis
Volume 2015, Article ID 958730, 9 pages
Research Article

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

Department of Mathematics, Indian Institute of Technology, Kharagpur 721302, India

Received 18 September 2014; Revised 31 December 2014; Accepted 5 January 2015

Academic Editor: Lukasz Stettner

Copyright © 2015 Arunava Maity and U. C. Gupta. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 CGG method). The CGG method 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.

1. 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 depends on the evaluation of the vector . For that we need to find the matrix. Once 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 , for example, cyclic reduction algorithm by Bini et al. [8, 9]. But all these proposed algorithms have 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 as well as state vectors in 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-death-modulated 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 in the unit disk (which are found to be real and distinct) using a modified bisection method and subsequently presents the matrix in terms of the zeros and associated right null vectors of . 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 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 , using roots inside the unit circle of the numerator of 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 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.

2. Model Description

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

2.1. Arrival Process

Here the arrival rate matrices corresponding to batch arrivals are tridiagonal and of the following form: where , , , , and the infinitesimal generator of the underlying process is given by ; that is, where , , , . In the region , is assumed to be analytic, so is , , . The stationary probability vector for the underlying process is given by

The fundamental arrival rate is given by , where is column vector of 1’s of appropriate dimension.

2.2. 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, , , , and . Traffic intensity is given by .

3. Analysis

In the literature (see Lucantoni [14]), the probability generating function (pgf) of the stationary probability vector just after service completion epochs for the queue is given by Here , where 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 type queue are as follows: one needs to evaluate the vector and after getting , one needs to invert the generating function to get the stationary probability vector. In the following sections, we briefly discuss the steps required for the analysis of queue by Nishimura’s method and by CGG method, respectively.

3.1. Analysis Done by 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 and the second one is the evaluation of steady state probability vector, which is as follows.

3.1.1. The Evaluation of

In the following, we describe the procedure adopted by Nishimura for the evaluation of .

(i) Nishimura shows that, for a BDMMAP/G/1 queue, the ’s (zeros in the unit disk of ) are distinct and simple and can be exhibited as where are the right null vectors of corresponding to .

(ii) ’s are calculated from the following equations: Accordingly, and are the eigenvalues and right eigenvectors of the matrix. is the matrix probability generating function of arrival numbers, that is, of 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 is calculated from and .

3.1.2. The Evaluation of Steady State Probability Vector

The procedure adopted by Nishimura for evaluating the steady state probability vectors is as follows.

(i) In order to derive the steady state probability vector, he has shown that the zeros of and are the same, simple, and real in (see Nishimura, [13, Theorems 1 and 4]), where . He has further assumed (see Nishimura, [13, page 434, Assumption 1]) that, for , eigenvalues of are simple.

(ii) With the help of the above facts and assumptions, he has represented the generating function in the form where (a)(b)where , with .(c)and ’s are calculated from the following iterative procedure:(1)(2)where is a complex number with and , where is a left eigenvector corresponding to the eigenvalue of which is a symmetric tridiagonal matrix, with .

(iii) After obtaining , he uses fast Fourier transform (FFT) to get the stationary probabilities. In order to apply FFT he has to determine a number (preferably ) so that after that number the remaining probabilities are negligible and further he obtains the approximated state vector as , . In order to calculate the eigenvalues of , 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 (as well as of ) 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 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 or the boundary vector by this method is not as difficult as the evaluation of the other stationary vectors.

Remark 3. In order to derive by Nishimura’s method, one needs to first find & . For that one has evaluate the eigenvalues and eigenvectors of , which is again a tedious task. After getting these eigenvectors, he uses a recursive algorithm to evaluate & , which may address numerical inaccuracy for high order matrices.

Remark 4. In order to apply FFT on , he has to calculate the eigenvalues of , , 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 instead of (as in the case of DFT) (if is of the form , a one-dimensional FFT of length requires less than 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 basic knowledge of queueing, spectral theory, and computer programming to computationally implement it.
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.

3.2. CGG Approach

This method exploits the roots of the equation that lies within the domain and exterior to ; we denote it as . In this method, they first find the arbitrary epoch probabilities and thereby obtain the departure epoch probabilities.

In the following, we briefly describe the entire CGG method in stepwise manner. For the sake of easy understanding, we divide the steps in two segments: first one describes the procedure to evaluate (as well as ) and the second one is about how they evaluate the other state vectors.

3.2.1. Evaluation of the Vector

Step 1. We write the pgf as (a) is the vector generating function of and is the th component of , that is, .(b) is a row vector whose th element denotes the joint probability that there are customers in the system and the phase of the arriving customer is at arbitrary epoch.(c) 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 is given by , where

Step 3. From the system of (15), using Cramer’s rule, we can find , . We denote and by and , respectively, where

Step 4. Each component of the VGF given in being analytic in implies that must be the zeros of the numerator of and hence we can determine the unknown vector by considering any one component of say . This implies that , and another equation comes from normalizing condition, that is, .

3.2.2. Evaluation of the State Probability Vector at Postdeparture Epoch

If has roots, namely, , outside the unit circle then by partial fraction method , . We get where, Now the system-length distribution at postdeparture epoch is obtained by

Remark 6. In Nishimura’s method, we have to compute all the eigenvalues and the corresponding eigenvectors of . Also we need to evaluate the eigenvalues and eigenvectors of and , 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 from any one component of the vector generating function which is a notable side of this approach.

Remark 7. Nishimura [13] shows that for queue the roots of the characteristic equation 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 [1720]). 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.

4. 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 , service time distribution as phase type distribution with parameters , with mean 0.067934. In Table 1, we obtain the distribution of stationary probability at arbitrary epoch using CGG method for . 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 matches . Here, we would like to mention that we have compared the values of 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 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.

Table 1: Distribution of stationary probability at arbitrary epoch using CGG method.
Table 2: Distribution of stationary probability at service completion epoch using CGG method.
Table 3: Comparison of the vector obtained by Nishimura and CGG method for the queue.

In Table 3, we present the boundary vector 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 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. 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 ) of the c.e. that lie in and these values agree up to satisfaction after (for details see the Table 4).

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

The purpose of the examples presented here is to compare the accuracy of the two methods (up to the level of ) 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.

5. 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 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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of the paper.


The authors thank the anonymous referees for their valuable comments and suggestions that helped to improve the presentation of the paper in the current form. The second author acknowledges the Department of Science and Technology, Govt. of India, for the financial support under the Project Grant SR/S4/MS: 789/12.


  1. M. F. Neuts, “A versatile Markovian point process,” Journal of Applied Probability, vol. 16, no. 4, pp. 764–779, 1979. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  2. D. M. Lucantoni, “The BMAP/G/1 queue: a tutorial,” in Performance Evaluation of Computer and Communication Systems, vol. 729 of Lecture Notes in Computer Science, pp. 330–358, Springer, Berlin, Germany, 1993. View at Publisher · View at Google Scholar
  3. S. R. Chakravarthy, “The batch Markovian arrival process: a review and future work,” in Advances in Probability Theory and Stochastic Processes, vol. 1, pp. 21–49, 2001. View at Google Scholar
  4. V. Ramaswami, “A stable recursion for the steady state vector in Markov chains of M/G/1 type,” Stochastic Models, vol. 4, no. 1, pp. 183–188, 1988. View at Publisher · View at Google Scholar
  5. V. Ramaswami, “Nonlinear matrix equations in applied probability—solution techniques and open problems,” SIAM Review, vol. 30, no. 2, pp. 256–263, 1988. View at Publisher · View at Google Scholar · View at MathSciNet
  6. D. Bini and B. Meini, “On the solution of a nonlinear matrix equation arising in queueing problems,” SIAM Journal on Matrix Analysis and Applications, vol. 17, no. 4, pp. 906–926, 1996. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  7. M. F. Neuts, Structured Stochastic Matrices of M/G/1 Type and Their Applications, vol. 5 of Probability: Pure and Applied, Marcel Dekker, New York, NY, USA, 1989. View at MathSciNet
  8. D. A. Bini, B. Meini, S. Steffé, and B. van Houdt, “Structured Markov chains solver: software tools,” in Proceedings of the 2006 Workshop on Tools for Solving Structured Markov Chains, p. 14, ACM, 2006.
  9. D. A. Bini and B. Meini, “Improved cyclic reduction for solving queueing problems,” Numerical Algorithms, vol. 15, no. 1, pp. 57–74, 1997. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  10. H. W. Lee, J. M. Moon, J. K. Park, and B. K. Kim, “A spectral approach to compute the mean performance measures of the queue with low-order BMAP input,” Journal of Applied Mathematics and Stochastic Analysis, vol. 16, no. 4, pp. 349–360, 2003. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  11. H. R. Gail, S. L. Hantler, and B. A. Taylor, “Spectral analysis of M/G/1 and G/M/1 type Markov chains,” Advances in Applied Probability, vol. 28, no. 1, pp. 114–165, 1996. View at Google Scholar
  12. H. W. Lee, J. M. Moon, B. K. Kim, J. G. Park, and S. W. Lee, “A simple eigenvalue method for low-order D-BMAP/G/1 queues,” Applied Mathematical Modelling, vol. 29, no. 3, pp. 277–288, 2005. View at Publisher · View at Google Scholar
  13. S. Nishimura, “A MAP/G/1 queue with an underlying birth-death process,” Stochastic Models, vol. 19, no. 4, pp. 425–447, 2003. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  14. D. M. Lucantoni, “New results on the single server queue with a batch Markovian arrival process,” Communications in Statistics: Stochastic Models, vol. 7, no. 1, pp. 1–46, 1991. View at Publisher · View at Google Scholar · View at MathSciNet
  15. D. A. Bini, “Numerical computation of polynomial zeros by means of Aberth's method,” Numerical Algorithms, vol. 13, no. 2, pp. 179–200, 1996. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  16. M. L. Chaudhry, G. Singh, and U. C. Gupta, “A simple and complete computational analysis of MAP/R/1 queue using roots,” Methodology and Computing in Applied Probability, vol. 15, no. 3, pp. 563–582, 2013. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  17. B. Fisher and M. Chaudhry, “Computing the distribution for the number of renewals with bulk arrivals,” INFORMS Journal on Computing, vol. 26, no. 4, pp. 885–892, 2014. View at Publisher · View at Google Scholar · View at MathSciNet
  18. N. Akar, “Solving the ME/ME/1 Queue with State-space Methods and the Matrix Sign Function,” Queueing Systems, vol. 63, no. 2, pp. 131–145, 2006. View at Google Scholar
  19. N. Akar and E. Arikan, “A numerically efficient method for the MAP/D/1/K queue via rational approximations,” Queueing Systems, vol. 22, no. 1-2, pp. 97–120, 1996. View at Publisher · View at Google Scholar · View at MathSciNet
  20. S. Winitzki, “Uniform approximations for transcendental functions,” in Computational Science and Its Applications—ICCSA, vol. 2667 of Lecture Notes in Computer Science, pp. 780–789, Springer, Berlin, Germany, 2003. View at Publisher · View at Google Scholar · View at MathSciNet
  21. S. Nishimura and M. Shinno, “Modeling burst traffic using a MAP with a tree structure,” Journal of the Operations Research Society of Japan, vol. 47, no. 3, pp. 129–144, 2004. View at Google Scholar · View at MathSciNet · View at Scopus