#### Abstract

This paper deals with a batch arrival infinite-buffer single server queue. The interbatch arrival times are generally distributed and arrivals are occurring in batches of random size. The service process is correlated and its structure is presented through a continuous-time Markovian service process (). We obtain the probability density function (p.d.f.) of actual waiting time for the first and an arbitrary customer of an arrival batch. The proposed analysis is based on the roots of the characteristic equations involved in the Laplace-Stieltjes transform (LST) of waiting times in the system for the first, an arbitrary, and the last customer of an arrival batch. The corresponding mean sojourn times in the system may be obtained using these probability density functions or the above LSTs. Numerical results for some variants of the interbatch arrival distribution (Pareto and phase-type) have been presented to show the influence of model parameters on the waiting-time distribution. Finally, a simple computational procedure (through solving a set of simultaneous linear equations) is proposed to obtain the “” matrix of the corresponding -type Markov chain embedded at a prearrival epoch of a batch.

#### 1. Introduction

In recent years, analysis of queueing processes with nonrenewal arrival and service processes has been carried out to model data transmission of complex computer and communication networks. The traditional queueing analysis using Poisson processes is not powerful enough to capture the dependence involved in the arrival (service) process. The performance analysis of such bursty and correlated type of the arrival (service) processes may be done through some analytically tractable processes, namely, Markovian arrival process and Markovian service process ; see Lucantoni et al. [1], Gupta and Banik [2], and references therein. The two processes and are convenient representations of a versatile Markovian point process or N-process; see Neuts [3] and Ramaswami [4]. Chaudhry et al. [5] obtain stationary system-length distributions for the infinite-buffer batch arrival queues at various epochs using roots. Further, Singh et al. [6] have shown that it is beneficial to use roots rather than the matrix-geometric/matrix-analytic methods in terms of computation time. For further details on the roots method, the readers are referred to the related references in Chaudhry et al. [7].

One may note that Grassmann [8] shows the relation between the characteristics roots and the corresponding eigenvalues of the possible quasi birth-death (QBD) process in connection with the computation of the waiting-time distribution in a queueing model. For further details on the relationships between the characteristics roots of a Markov chain and the corresponding “” matrix of -type Markov chain or “” matrix of -type Markov chain one may see Gail et al. [9].

In this paper we obtain the explicit closed-form expressions for the Laplace-Stieltjes transform (LST) of waiting-time distribution. For analytic purpose, we start our work with associated prearrival epoch probability distribution; see Chaudhry et al. [5] for details. After that we consider characteristic equations involved in the Laplace-Stieltjes transform’s (LST’s) waiting-time distribution for the first and an arbitrary customer of an arrival batch. Once the roots of these characteristic equations are found, it becomes easy to obtain the stationary waiting-time distribution for the first and an arbitrary customer of an arrival batch. Also, a simple computational procedure based on the solution of a set of linear simultaneous equations has been proposed to compute the matrix of the corresponding -type Markov chain embedded at a prearrival epoch of a batch. It is well-known that the determination of the “” matrix can be done through the minimal nonnegative solution of a nonlinear matrix equation; see Neuts [10] for details. Finally, it may be remarked here that the same queueing model has been dealt with in [5]. But in both the queueing models, corresponding probability density functions of actual waiting time for the first/an arbitrary customer have not been derived. It is needless to mention here that probability density functions of the waiting times are important for obtaining statistical measures of Quality of Service (QoS) requirement of a queueing system. Thus, the present work is a nontrivial extension of [5].

#### 2. Description of the Model: Queue and the Analysis of the Stationary Distribution of Waiting Time

We assume a renewal input single server queueing system where the customers arrive in batches of random size whose probability mass function is given by , , and mean batch size . Let be the probability generating function of . The interbatch arrival times are independent identically distributed (i.i.d.) random variables with distribution function (DF) , probability density function with , and Laplace-Stieltjes transform (LST) , where Then, ; i.e., is the mean interbatch arrival time, where is the -th derivative of at .

A single server is serving the customers according to a continuous-time Markovian service process () with matrix representation . For a detailed discussion on , the readers are referred to Chaudhry et al. [5]. Let us denote with being an irreducible infinitesimal generator of the underlying Markov chain . The mean stationary service rate of the is given by , where is the probability row vector which can be computed from with , where is a column vector of ones with an appropriate dimension. The customers are served singly according to a with an average stationary service rate . The offered load is given by .

In the following section, we obtain the actual waiting-time distribution for the first customer and an arbitrary customer of an arrival batch. In order to get the actual waiting-time distributions, we first obtain the prearrival epoch probability distribution of system-length at a prearrival epoch. For this, we proceed exactly in the same way as that discussed by Chaudhry et al. [5]. Let represent the prearrival epoch probability distribution of an arrival batch when there are customers in the system and the server is busy in phase . Also, let denote the prearrival epoch probability that the system is empty while the server is idle in phase . Let be the row vector whose -th component is . Also, let us assume with being the -th component of . As presented in [5], is an analytic function of for . Thus, we have where are constants to be determined. Here, are the roots inside the unit circle of the equation with being the maximum size of an arrival batch, where with The procedure of obtaining is available in literature; e.g., see Appendix 2 of [5]. One may note that the computation of using above relation may be cumbersome. However, the following numerical scheme may be efficient and is given by , where may be obtained as proposed in [5]. The equation has roots inside the unit circle. That is, the equation has roots outside the unit circle. Now, comparing the coefficient of from both sides of (1), one may have the explicit values of in terms of and ; see Chaudhry et al. [5] for a detailed procedure.

#### 3. Waiting-Time Analysis

In the following we obtain the probability density function (p.d.f.) of actual waiting time for the first and an arbitrary customer of an arrival batch. Let the -th entry of be the LST of the service times of a total of customers and the corresponding phase change of the underlying Markov chain is from to . The probability that the service of a customer is completed in the interval along with phase changes is given by the matrix , where represents a function of such that as Since the LST of total service time of customers is the product of their individual LST of each service times of those customers, therefore,

where Let , , denote the row vector of order , where represents the stationary joint probability that the first customer of an arrival batch has to wait up to a maximum of units of time before getting served with the service process being in phase . The first customer in an arrival batch waits in the system at most time units if and only if he sees upon arrival customers ahead of him and so he has to wait for service completions including itself in the interval Similarly, we denote , , and , , as probability vectors whose -th component of denotes the stationary joint probability that an arbitrary customer in an arrival batch has to wait in the system a maximum of units of time before departing the system with the service process being in phase . The th component of denotes the stationary joint probability that the last customer in an arrival batch has to wait in the system at most units of time before getting served and departing the system with the phase of the service process being . Let be the LST of . Similarly, let and be the LSTs of and , respectively. Hereafter, we need to define a random variable which denotes the number of customers before an arbitrary customer in a batch. Following Chaudhry and Templeton [7], one may write Therefore, the LST of the waiting-time distribution after a little algebraic manipulation (see Chaudhry et al. [5]) may be derived aswith One can find component-wise mean waiting times in the system of the above three types

, , and , respectively) by differentiating the above formulae and putting They are obtained after a little algebraic manipulation; see Chaudhry et al. [5] for explicit formulae of , , and . Using Little’s law, we can obtain mean sojourn time as , where is the mean stationary system-length obtained in [5] and one may note that while performing numerical computation. Also, one may simplify (3) as stated below. First combine the roots and , yielding a geometric series, which are then summed over . As a consequence of the formula regarding geometric series, this sum is , and this expression is expanded, using the adjoint matrix and determinants (see also [9], (35), and Section A.1 for a proof of formula (6) given below). with In the following, we develop a numerical scheme for inverting the above LSTs component-wise. Let us first look at the -th component of the LST given by formula (3) which is further simplified as (6). The right-hand side of (3) involves which is further simplified as the right-hand side of (6). Both LSTs (3) and (6) are analytic in the region Therefore (as is the LST of the transition probability matrix of the underlying Markov chain of ), the -th component of may be written as a rational function in , which is evident from the right-hand side of (6). Now, taking into account the zeros of the denominator of each terms in the right-hand side of (6), we obtain characteristic equations; namely, Also, each element of the matrix has denominator which has zeros in the region (see similar proof of Theorem 7.3 of [5]). As a result, has denominator which has zeros in the region with each zero having multiplicity . Similarly, the numerator of also has zeros in the region for each . After canceling the common zeros of the numerator and the denominator of × , one is left with zeros of the numerator of × in the region . Let us define as the zeros of the denominator of each component of in the region . Also, it may be noted here that none of these zeros is equal to a zero of . Therefore, one may finally note that the denominator of the th component of the LST given by formula (3) which has zeros: with . One may finally obtain the th component of using partial fractions as follows: where are constants to be evaluated and . It may be remarked here that while making partial fractions, we have assumed that the roots are distinct. To get the unknown constants, we use the required number of linear equations after substituting any positive values of () in (8). Finally, we equate the component-wise from (3) with (8) under the same positive values of . One may note that, as normalizing conditions, the component-wise values of and the components of must be included in the set of linear equations to get the above constants. These components may be obtained by putting in (8) and we may use , where is the first differentiation of with respect to . After that, taking inverse Laplace transform of (8), component-wise explicit closed-form expression of probability density function (p.d.f.) is given by with Finally, the p.d.f. of waiting time of the first customer may be obtained from . From this, we can get , , and .

Similar to the procedure followed earlier, we consider the zeros of the denominator of each element of the matrix in (4). As discussed earlier, the denominator of each element of has zeros in the region which may be obtained by solving . This implies that, given maximum arrival batch size , zeros of the denominator of each element of the matrix are repeated up to times. Let the roots of be denoted by Therefore, the denominator of each element of the matrix has repeated roots with multiplicity Also, we have to investigate the zeros of the denominator of each element of the matrix . These zeros may be obtained by using . It can also be shown that there will be zeros of the numerator of in the region (since is root of the equation ). Also, there are zeros of the denominator (which is equal to of in the region and they are with multiplicity each. After canceling the common zeros of the numerator and the denominator of , one is left with zeros of the numerator of in the region . Let these roots of the numerator of be denoted by , where Lastly, the roots of the denominator of each element of the matrix in the right-hand side of (4) are also . Finally, if we consider all the roots (with ) of denominators involved in each element of the vector and matrices occurring in the right-hand side of (4), we reach the following conclusion: ’s roots repeated in the region with multiplicity occur in the denominator of each component of in (4). On the other hand and are the roots in the region that occur in the denominator of each component of in (4). Following this conclusion, after using partial fractions, we may write the -th component of (in (4) which also has factor ) as follows: where , , and are constants to be evaluated and . It may be remarked here that while making partial fractions, we have assumed that the roots and are distinct. If they are repeated, a slight modification in partial fraction will be needed. Also, throughout the above discussions of roots, we discard any roots of the denominator of the above LSTs (3)-(5) in the region . The roots of the denominator with must be canceled with numerator as the LSTs (3)-(5) are analytic in the region .

Now to obtain the unknown constants, we use the required number of linear equations after substituting any positive values of in (10). Finally, we equate the components of from (4) with the corresponding component of (10) under the same positive values of . One may note that, as normalizing conditions, the component-wise values of and the components of must be included in the set of linear equations to get the above constants. These components may be obtained by putting in (10) and we may use , where is the first differentiation of with respect to . After that, taking inverse Laplace transform of (10), component-wise explicit closed-form expression of probability density function (p.d.f.) is given by with Finally, the p.d.f. of waiting time of an arbitrary customer may be obtained by . Further, let us denote the component-wise and total distribution function of the waiting-time distribution of an arbitrary customer of an arrival batch by , , and . It may be finally noted that must be equal to the -th component of . Exactly similar analysis may be carried out to obtain the component-wise p.d.f. of the waiting time of the last customer in an arrival batch using the LST given by (5).

#### 4. Numerical Results

Using the results formulated in the previous section, some numerical results are presented in this section. Although numerical calculations were carried out with high precision, due to lack of space the results have been reported to 6 decimal places. Since various distributions can be either represented or approximated by -distribution, we take interbatch arrival time distribution to be of -type having the representation , where and are of dimension . Then using the procedure adopted by Chaudhry et al. [5], can be derived as follows: with , where and are Kronecker product and sum, respectively. For the derivation of , see Theorem 7.5 in [5]. For more discussion on the derivation of , see [5]. Numerical computations were performed on a PC having Intel(R) Core i5 processor @1.8 GHz with 4 GB DDR3 RAM using MAPLE 2015.

##### 4.1. The Waiting-Time Analysis of a Queueing System

In Figure 1, we have computed and plotted the component-wise p.d.f. of waiting times in the system of an arbitrary customer in a queueing system using the numerical scheme presented in Section 3. The -type interbatch arrival time is given by and with . The batch size distribution is given by , and with mean batch size The matrices having positive lag-1 correlation coefficient 0.271908 as with stationary mean service rate and which gives traffic intensity (offered load) The characteristic roots to obtain the system-length distribution at a prearrival epoch are calculated through the nonlinear algebraic equation: where and are given above and may be obtained by (12); see Theorem 7.5 of [5]. The roots of (14) inside are evaluated and the corresponding values are calculated using Section 2 where the prearrival epoch probabilities are derived. Details of the characteristic equations and their roots used for the evaluation of waiting-time densities are available with the authors and may be supplied upon specific request to the authors. After computation of the required roots, we construct partial fractions for each component of as given in (10). Finally, using inverse Laplace-transforms we obtain component-wise p.d.f. of waiting times in the system. In Figure 1, we have plotted the above three component-wise densities along with the total waiting-time density function against Figure 1 shows that component-wise density of waiting time behaves significantly different from that of . Also, in Figure 2, we have plotted the corresponding component-wise distribution function as well as total distribution function against time. Similar behaviour may be observed in the case of DFs, which is what is expected.

##### 4.2. The Waiting-Time Analysis of a Queueing System

In Figure 3, we have repeated similar experiment as given in above figures; i.e., we plotted component-wise p.d.f.’s against time in a queue, where stands for Pareto distribution whose probability density function and distribution functions are given by , , , , and , , respectively, with shape parameter and scale parameter . Batch size distribution is taken arbitrary with probability generating function and mean batch size The service matrices having positive lag-1 correlation coefficient 0.018690 are taken as with and . One may note that the moments do not exist for this case of Pareto interbatch arrival time distribution. From this Pareto distribution function (through equating ), we calculate median which is equal . Here we assume so that We use geometric transform approximation method (G-TAM); see Shortle et al. [11]. For this, consider probabilities Substituting in the distribution function, i.e., , leads to The weights are , , and . A binary search on to match (which is the median of an interbatch arrival time) gives Thus, the approximate Laplace transform of interbatch arrival time may be taken as The Padé approximation [4/5] of (16) is given as follows:

Now using the above procedure as discussed in previous section, we have computed the component-wise densities of waiting time for the using exactly the same procedure as described in Figure 1. Similar to Figure 1, in Figure 3, we have plotted the above three component-wise densities as well as versus . Here also we have seen significant difference among the component-wise densities. Also, we have plotted the corresponding distribution function in component-wise as well as total distribution function in Figure 4.

##### 4.3. Comparison of the Numerical Results Obtained through Roots and Matrix-Geometric Method

In this section we made a thorough numerical comparison of our results with that of Samanta [12] and our computed results are appended in Table 1. For this comparison, we have obtained component-wise waiting-time distribution function (using the procedure discussed in Section 3) of an infinite-buffer - queueing model with parameters given in Table 2 of [12]. For the same queueing model, we have obtained waiting-time (in system) distribution function (using the procedure discussed in Section 3) which was presented in Table 2 of [12] using matrix-geometric method. As expected, in both the above cases, the results were matched up to six decimal places against Table 2 of [12]. In Table 1, we have presented first seven and the last values of the distribution function. For this experiment, - service matrices with positive lag-1 correlation coefficient 0.5381025 are taken aswith stationary mean service rate and = . The interarrival time is phase-type whose representation was taken aswith which gives . It may be noted here that for an infinite-buffer queue the expression for LST of waiting-time distribution of an arriving customer is the same as which is given in (3) and for the sake of notational convenience we may denote this as .

In Figures 5 and 6, we have plotted the probability densities and distribution functions, respectively, against time for the queueing model discussed in Table 1.

#### 5. Conclusion

In this paper, we analyze the waiting-time distribution of the queue and obtain steady-state probability distribution of the waiting time of the first customer and an arbitrary customer of an arrival batch. The proposed method is based on the roots of a characteristic equation involved in finding probability distribution at a prearrival epoch. After that, we use prearrival epoch probabilities and the roots of a set of characteristic equations involved in the LSTs of the waiting time for the first customer and an arbitrary customer of an arrival batch. Besides the computational advantages in terms of time, this method has other advantages over the matrix-geometric method in the sense that it is accurate even if one does not use all the roots. For example, the values of the distribution function presented in Table 1 of the previous section may be obtained (with an accuracy of up to 6 decimal places) by just using the first 8 roots from the list of roots. As a consequence, we may obtain good approximate waiting-time distribution using less number of roots. Finally, we discuss a simple computational procedure (based on the solution of a set of simultaneous linear equations) which can compute the “” matrix of the corresponding -type Markov chain with a desired level of accuracy; see Section A.2. This procedure may serve as an alternative method for computation of the matrix which is usually computed through the nonnegative solution of a nonlinear matrix equation. Similar procedure may be used to obtain the waiting-time distribution of a --type queue with set-up time with or without vacation(s). These problems may attract researchers’ attention in the future.

#### Appendix

#### A. A Proof and Some Additional Numerical Results Related to the Paper Using Roots

##### A.1. Proof of Formula (6) of the Paper

which yields the result (6).

##### A.2. Demonstration of a Simple Computational Procedure for “” Matrix (Used in Matrix-Geometric Method) Related to the Queueing Model **-** of the Paper Using Roots

As discussed by Neuts [10], the stationary probabilities at an embedded prearrival epoch in a -type Markov chain may be obtained in matrix-geometric forms. That is, for a -type Markov chain Neuts [10] showed that there exists a nonnegative matrix “” with appropriate dimension such that the stationary probability vectors at an embedded prearrival epoch are given by , , where the matrix “” has spectral radius in case of a positive recurrent Markov chain and the matrix “” may be computed as a nonnegative solution of a nonlinear matrix equation (see below). Gail et al. [9] proved that the eigenvalues of the “” matrix are related to the characteristic roots of the corresponding -type Markov chain. In the following a simple computation procedure for obtaining the “” matrix as a solution of a set of simultaneous linear equations has been proposed for and queueing models.

First, we explain this computational procedure for computing the matrix in the case of queueing model presented in Table 1. Following [10], the stationary probability vectors at an embedded Markov prearrival epoch are given below. Since we have assumed and the Markov chain is irreducible, there exists the minimal nonnegative solution of the matrix equation , where is strictly positive and of spectral radius less than one. The matrix can be calculated by an iterative scheme given below: where is the value of at the iteration. The invariant probability vector of the Markov chain is obtained by where is obtained by solving the left eigenvector of the matrix with the normalizing condition Now we obtain the roots of the characteristics equation , as , and Thereafter, we get the stationary probability vectors , at an embedded prearrival epoch using the roots , and . One may note that in this case the has dimension , i.e., unknown elements, and they may be obtained by solving the following simultaneous linear equations derived from (A.3) as follows: Equating , , and (A.4)-(A.7) component-wise, we obtain the simultaneous equations with elements of the matrix as variables. After solving these simultaneous linear equations with unknown elements of the matrix , we obtain the matrix as follows: The correctness of the matrix is checked by computing the eigenvalues of given in (A.8) and the eigenvalues are obtained as , , , and which are exactly the same as the characteristic roots. Also, we obtained matrix using iterative scheme (A.2) and found that the matrix is exactly the same as that given in (A.8) after iterations.

Next, we investigate the same computational procedure for computing the matrix in the case of queueing model presented in Section 4.1. It may be noted here that for the case of queueing model as discussed by Neuts [10] (Remarks on bottom of pg. 189), the may be solved by the following set of simultaneous linear equations. One may note that for this case the matrix computed by the above procedure (presented below) has spectral radius less than Here the eigenvalues of the matrix are not the same as the characteristics roots in this case. The eigenvalues of are given by , , , , , , , and . Also, we have checked that this satisfies any of the following equations: Therefore, by the above simple computational procedure, one may compute the matrix of the corresponding -type using the stationary probabilities (obtained through characteristic roots) at an embedded prearrival epoch.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The second author was supported partially by NSERC under research grant number RGPIN-2014-06604.