A Note on the Waiting-Time Distribution in an Infinite-Buffer GI[X]/C-MSP/1 Queueing System

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-timeMarkovian service process (C-MSP).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 “R” matrix of the corresponding GI/M/1-type Markov chain embedded at a prearrival epoch of a batch.


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 /-/1 queues at various epochs using roots.Further, Singh et al. [6] have shown that it is beneficial to use roots rather than the matrixgeometric/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   /  /1 queueing model.For further details on the relationships between the characteristics roots of a Markov chain and the corresponding "G" matrix of //1type Markov chain or "R" matrix of //1-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 waitingtime 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 R matrix of the corresponding //1-type Markov chain embedded at a prearrival epoch of a batch.It is well-known that the determination of the "R" 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].
A single server is serving the customers according to a continuous-time Markovian service process (C-) with matrix representation (L 0 , L 1 ).For a detailed discussion on -, the readers are referred to Chaudhry et al. [5].Let us denote L() = L 0 + L 1  with L ≡ L(1) = L 0 + L 1 being an irreducible infinitesimal generator of the underlying Markov chain {()}.The mean stationary service rate  ⋆ of the - is given by  ⋆ = L 1 e, where  = [ 1 ,  2 , . . .,   ] is the probability row vector which can be computed from L = 0 with e = 1, where e 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 waitingtime 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  (≥ 1) customers in the system and the server is busy in phase .Also, let  −  (0) denote the prearrival epoch probability that the system is empty while the server is idle in phase .Let  − () ( ≥ 1) be the row vector whose th (1 ≤  ≤ ) component is  −  ().Also, let us assume  − * () = ∑ ∞ =0  − ()  with  − *  () being the -th 1 ≤  ≤  component of  − * ().As presented in [5],  − *  () is an analytic function of  for || ≤ 1.Thus, we have where   are constants to be determined.Here,   (1 ≤  ≤ r) are the r roots inside the unit circle of the equation det[I  −( −1 )S()] = 0 with r being the maximum size of an arrival batch, where S() = ∑ ∞ =0 S    =  * (−L()) with || ≤ 1.The procedure of obtaining S  is available in literature; e.g., see Appendix 2 of [5].One may note that the computation of S() using above relation may be cumbersome.However, the following numerical scheme may be efficient and is given by S() = lim →∞ ∑  =0 S    , where S  may be obtained as proposed in [5].The equation det[I   r − ( 1  r−1 +  2  r−2 +  3  r−3 +⋅ ⋅ ⋅+ r)S()] = 0 has r roots   inside the unit circle.That is, the equation det[I  − ()S( −1 )] = 0 has r roots 1/  outside the unit circle.Now, comparing the coefficient of   from both sides of (1), one may have the explicit values of  −  () (1 ≤  ≤ ,  ≥ 0) in terms of   (1 ≤  ≤ r) and   (1 ≤  ≤ r, 1 ≤  ≤ ); see Chaudhry et al. [5] for a detailed procedure.

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  (≥ 1) 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  L 0  L 1  + (), where () represents a function of  such that ()/ → 0 as  → 0. Since the LST of total service time of  customers   () is the product of their individual LST of each service times  1 () of those  customers, therefore, where R() ≥ 0. Let W  () = [ 1 (),  2 (), . ..,   ()],  ≥ 0, denote the row vector of order , where   () (1 ≤  ≤ ) 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  (≥ 0) customers ahead of him and so he has to wait for ( + 1) service completions including itself in the interval (0, ].Similarly, we denote W  () = [ 1 (),  2 (), . . .,   ()],  ≥ 0, and W  () = [ 1 (),  2 (), . . .,   ()],  ≥ 0, as probability vectors whose -th component   () (1 ≤  ≤ ) of W  () 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   () (1 ≤  ≤ ) of W  () 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 .
] be the LSTs of W  () and W  (), 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 waitingtime distribution after a little algebraic manipulation (see Chaudhry et al. [5]) may be derived as with R() ≥ 0. One can find component-wise mean waiting times in the system of the above three types and W  = [ 1 ,  2 , . . .,   ], respectively) by differentiating the above formulae and putting  = 0.They are obtained after a little algebraic manipulation; see Chaudhry et al. [5] for explicit formulae of W  , W  , and W  .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   () = W  e while performing numerical computation.Also, one may simplify (3) as stated below.First combine the roots    (1 ≤  ≤ r) and [ 1 ()]  , yielding a geometric series, which are then summed over .As a consequence of the formula regarding geometric series, this sum is [I  −    1 ()] −1 , 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 R() ≥ 0. 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  − * ( 1 ()) 1 () which is further simplified as the right-hand side of ( 6).Both LSTs (3) and ( 6) are analytic in the region R() ≥ 0. Therefore (as  1 () is the LST of the transition probability matrix of the underlying Markov chain of -), the -th (1 ≤  ≤ ) component of  − * ( 1 ()) 1 () 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 r terms in the right-hand side of (6), we obtain r characteristic equations; namely, Also, each element of the matrix  with R() < 0. One may finally obtain the th (1 ≤  ≤ ) component of W *  () using partial fractions as follows: where    (1 ≤  ≤ , 1 ≤  ≤ r, 1 ≤  ≤ ) are constants to be evaluated and R() ≥ 0. 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  (R() > 0) in (8).Finally, we equate the component-wise W *  () from ( 3) with ( 8) under the same positive values of .One may note that, as normalizing conditions, the component-wise values of W *  (0) and the components of W  must be included in the set of linear equations to get the above constants.These components may be obtained by putting  = 0 in (8) and we may use   = − * (1)   ()| {=0} , where  * (1)  () is the first differentiation of  *  () with respect to .After that, taking inverse Laplace transform of (8), componentwise explicit closed-form expression of probability density function (p.d.f.)   () is given by with  > 0. Finally, the p.d.f. of waiting time of the first customer may be obtained from   () = ∑  =1   ().From this, we can get   () = ∫  0   (), (1 ≤  ≤ ), and   () = ∫  0   ().Similar to the procedure followed earlier, we consider the zeros of the denominator of each element of the matrix (I  − ( 1 ())) in (4).As discussed earlier, the denominator of each element of  1 () has  zeros in the region R() < 0 which may be obtained by solving det[I  − L 0 ] = 0.This implies that, given maximum arrival batch size r,  zeros of the denominator of each element of the matrix (I  −( 1 ())) are repeated up to r times.Let the  roots of det[I  −L 0 ] = 0 be denoted by   (1 ≤  ≤ ).Therefore, the denominator of each element of the matrix (I  − ( 1 ())) has repeated roots   (1 ≤  ≤ ) with multiplicity r.Also, we have to investigate the zeros of the denominator of each element of the matrix (I  −  1 ()) −1 .These zeros may be obtained by using det[I  −  1 ()].It can also be shown that there will be  2 zeros of the numerator of det[I  −  1 ()] in the region R() ≤ 0 (since  = 0 is root of the equation det[I  −  1 ()] = 0).Also, there are  2 zeros of the denominator (which is equal to {det[I  − L 0 ]}  ) of det[I  −  1 ()] in the region R(s) < 0 and they are   (1 ≤  ≤ ) with multiplicity  each.After canceling the common zeros of the numerator and the denominator of det[I  −  1 ()], one is left with ( − 1) zeros of the numerator of det[I  −  1 ()] in the region R() < 0. Let these ( − 1) roots of the numerator of det[I  −  1 ()] be denoted by ℎ  (1 ≤  ≤  − 1), where R(ℎ  ) < 0. Lastly, the roots of the denominator of each element of the matrix  1 () in the right-hand side of (4) are also   (1 ≤  ≤ ).Finally, if we consider all the roots (with R() < 0) 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 (1 ≤  ≤ ) roots repeated in the region R() < 0 with multiplicity (r − 1) occur in the denominator of each component of W *  () in (4).On the other hand   (1 ≤  ≤ , 1 ≤  ≤ r) and ℎ  (1 ≤  ≤  − 1) are the roots in the region R() < 0 that occur in the denominator of each component of W *  () in (4).Following this conclusion, after using partial fractions, we may write the -th (1 ≤  ≤ ) component of W *  () (in (4) which also has factor  − * ( 1 ())) as follows: where    (1 ≤ ,  ≤ , 1 ≤  ≤ r − 1),   , (1 ≤ ,  ≤ , 1 ≤  ≤ r), and    (1 ≤  ≤ , 1 ≤  ≤  − 1) are constants to be evaluated and R() ≥ 0. 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 R() ≥ 0. The roots of the denominator with R() ≥ 0 must be canceled with numerator as the LSTs (3)-( 5) are analytic in the region R() ≥ 0.
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 W *  () 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 W *  (0) and the components of W  must be included in the set of linear equations to get the above constants.These components may be obtained by putting  = 0 in (10) and we may use   = − * (1)   ()| {=0} , where  * (1)  () is the first differentiation of  *  () with respect to .After that, taking inverse Laplace transform of (10), componentwise explicit closed-form expression of probability density function (p.d.f.)   () is given by with  > 0. Finally, the p.d.f. of waiting time of an arbitrary customer may be obtained by   () = ∑  =1   ().Further, let us denote the component-wise and total distribution function of the waiting-time distribution of an arbitrary 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).

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 (, T), where  and T are of dimension ].Then using the procedure adopted by Chaudhry et al. [5], S() can be derived as follows: with L() ⊕ T = (L() ⊗ I ] ) + (I  ⊗ T), where ⊗ and ⊕ are Kronecker product and sum, respectively.For the derivation of S(), see Theorem 7.5 in [5].For more discussion on the derivation of S(), 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. [] /-/1/∞ 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  [] /-/1/∞ queueing system using the numerical scheme presented in Section 3. The -type interbatch arrival time is given by  = [0.350.65] and T = [ −1.183 2.423 1.367 −2.896 ] with  = 0.032317.The batch size distribution is given by  1 = 0.0,  2 = 0.0,  3 = 0.5, and  4 = 0.5,   = 0.0 ( ≥ 5) with mean batch size  = 3.

The Waiting-Time Analysis of a 𝑃𝐻
with stationary mean service rate  * = 0.500000 and  = [0.4649780.428426 0.106596] which gives traffic intensity (offered load)  = 0.226216.The characteristic roots to obtain the system-length distribution at a prearrival epoch are calculated through the nonlinear algebraic equation: where  = 3, r = 4 and   ( ≥ 1) are given above and S() may be obtained by (12); see Theorem 7.5 of [5].The r = 12 roots of (14) inside || < 1 are evaluated and the corresponding   (1 ≤  ≤ 12, 1 ≤  ≤ 3) 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 W *  () 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 componentwise densities along with the total waiting-time density function   () against . Figure 1 shows that componentwise 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. [] /-/1/∞ 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  [] /-/1/∞
(17) Now using the above procedure as discussed in previous section, we have computed the component-wise densities of waiting time for the  [] /-/1/∞ 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

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 /-/1 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 with  = 3.198147 which gives  = 0.593567.It may be noted here that for an infinite-buffer /-/1 queue the expression for LST of waiting-time distribution of an arriving customer is the same as W *  () = W *  () which is given in (3) and for the sake of notational convenience we may denote this as W *  ().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.

Conclusion
In this paper, we analyze the waiting-time distribution of the  [] /-/1/∞ 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  2 = 16 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 "R" matrix of the corresponding //1type Markov chain with a desired level of accuracy; see Section A.2.This procedure may serve as an alternative method for computation of the R 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 /-P/1-type queue with set-up time with or without vacation(s).These problems may attract researchers' attention in the future.

Table 1 :
Actual waiting-time (in the system) distribution of a customer.
[12]e2of[12].In Table1, 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 as