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

Authors: A. D. Banik; M. L. Chaudhry; James J. Kim
Journal: Journal of Probability and Statistics
 (2018)

---

## 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 (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.


---

## Body

## 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(C-MAP) and Markovian service process (C-MSP); see Lucantoni et al. [1], Gupta and Banik [2], and references therein. The two processes C-MAP and C-MSP 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 GI/C-MSP/1 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 En/Em/1 queueing model. For further details on the relationships between the characteristics roots of a Markov chain and the corresponding “G” matrix of M/G/1-type Markov chain or “R” matrix of GI/M/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 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 R matrix of the corresponding GI/M/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].


## 2. Description of the Model:GI[X]/C-MSP/1/∞ 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 sizeX whose probability mass function is given by P(X=i)=gi, i=1,2,…, and mean batch size E(X)=∑i=1∞igi=g¯. Let G(z)=∑i=1∞gizi be the probability generating function of gi. The interbatch arrival times are independent identically distributed (i.i.d.) random variables with distribution function (DF) A(x), probability density function a(x) with x>0, and Laplace-Stieltjes transform (LST) A∗(s), where R(s)≥0. Then, 1/λ=-A∗(1)(0); i.e., 1/λ is the mean interbatch arrival time, where A∗(j)(θ) is the j-th (j≥1) derivative of A∗(s) at s=θ.A single server is serving the customers according to a continuous-time Markovian service process (C-MSP) with matrix representation (L0,L1). For a detailed discussion on C-MSP, the readers are referred to Chaudhry et al. [5]. Let us denote L(z)=L0+L1z with L≡L(1)=L0+L1 being an irreducible infinitesimal generator of the underlying Markov chain {J(t)}. The mean stationary service rate μ⋆ of the C-MSP is given by μ⋆=π¯L1e, where π¯=π¯1,π¯2,…,π¯m 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 C-MSP with an average stationary service rate μ∗. The offered load ρ is given by ρ=λg¯/μ∗.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 πj-(n) represent the prearrival epoch probability distribution of an arrival batch when there are n(≥1) customers in the system and the server is busy in phase j. Also, let πj-(0) denote the prearrival epoch probability that the system is empty while the server is idle in phase j. Let π-(n)(n≥1) be the row vector whose j-th (1≤j≤m) component is πj-(n). Also, let us assume π-∗(z)=∑n=0∞π-(n)zn with πj-∗(z) being the j-th 1≤j≤m component of π-∗(z). As presented in [5], πj-∗(z) is an analytic function of z for |z|≤1. Thus, we have (1)πj-∗z=∑i=1mr^kij1-γiz,1≤j≤m,where kij are constants to be determined. Here, γi(1≤i≤mr^) are the mr^ roots inside the unit circle of the equation det⁡Im-Gz-1Sz=0 with r^ being the maximum size of an arrival batch, where S(z)=∑n=0∞Snzn=A∗-Lz with z≤1. The procedure of obtaining Sn is available in literature; e.g., see Appendix 2 of [5]. One may note that the computation of S(z) using above relation may be cumbersome. However, the following numerical scheme may be efficient and is given by Sz=limN→∞⁡∑n=0NSnzn, where Sn may be obtained as proposed in [5]. The equation det⁡Imzr^-g1zr^-1+g2zr^-2+g3zr^-3+⋯+gr^Sz=0 has mr^ roots γi inside the unit circle. That is, the equation det⁡Im-GzSz-1=0 has mr^ roots 1/γi outside the unit circle. Now, comparing the coefficient of zn from both sides of (1), one may have the explicit values of πj-(n)(1≤j≤m,n≥0) in terms of γi(1≤i≤mr^) and krj(1≤i≤mr^,1≤j≤m); 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(i,j)-th entry of ϕk(s) be the LST of the service times of a total of k(≥1) customers and the corresponding phase change of the underlying Markov chain is from i to j. The probability that the service of a customer is completed in the interval x,x+dx along with phase changes is given by the matrix eL0xL1dx+o(dx), where o(x) represents a function of x such that o(x)/x→0 as x→0. Since the LST of total service time of k customers ϕk(s) is the product of their individual LST of each service times ϕ1(s) of those k customers, therefore, (2)ϕ1s=∫0∞e-sxeL0xL1dx=sIm-L0-1L1andϕks=ϕ1ks,k≥2,whereR(s)≥0. Let WF(x)=[WF1(x),WF2(x),…,WFm(x)], x≥0, denote the row vector of order m, where WFr(x)(1≤r≤m) represents the stationary joint probability that the first customer of an arrival batch has to wait up to a maximum of x units of time before getting served with the service process being in phase r. The first customer in an arrival batch waits in the system at most x time units if and only if he sees upon arrival n(≥0) customers ahead of him and so he has to wait for (n+1) service completions including itself in the interval 0,x. Similarly, we denote WA(x)=[WA1(x),WA2(x),…,WAm(x)], x≥0, and WL(x)=[WL1(x),WL2(x),…,WLm(x)], x≥0, as probability vectors whose r-th component WAr(x)(1≤r≤m) of WA(x) denotes the stationary joint probability that an arbitrary customer in an arrival batch has to wait in the system a maximum of x units of time before departing the system with the service process being in phase r. The rth component WLr(x)(1≤r≤m) of WL(x) denotes the stationary joint probability that the last customer in an arrival batch has to wait in the system at most x units of time before getting served and departing the system with the phase of the service process being j. Let WF∗(s)=[WF1∗(s),WF2∗(s),…,WFm∗(s)] be the LST of WF(x). Similarly, let WA∗(s)=[WA1∗(s),WA2∗(s),…,WAm∗(s)] and WL∗(s)=[WL1∗(s),WL2∗(s),…,WLm∗(s)] be the LSTs of WA(x) and WL(x), respectively. Hereafter, we need to define a random variable J which denotes the number of customers before an arbitrary customer in a batch. Following Chaudhry and Templeton [7], one may write gr-=P[J=r]=(1/g¯)∑i=r+1∞gi,r≥0. Therefore, the LST of the waiting-time distribution after a little algebraic manipulation (see Chaudhry et al. [5]) may be derived as(3)WF∗s=∑n=0∞π-nϕ1n+1s=π-∗ϕ1sϕ1s,(4)WA∗s=∑n=0∞π-n∑j=0∞gj-ϕ1n+j+1s=1g¯π-∗ϕ1sIm-Gϕ1sIm-ϕ1s-1ϕ1s,(5)WL∗s=∑n=0∞π-n∑j=1∞gjϕ1n+js=π-∗ϕ1sGϕ1s,with R(s)≥0. One can find component-wise mean waiting times in the system of the above three types( W F = [ W F 1 , W F 2 , … , W F m ], WA=[WA1,WA2,…,WAm], and WL=[WL1,WL2,…,WLm], respectively) by differentiating the above formulae and putting s=0. They are obtained after a little algebraic manipulation; see Chaudhry et al. [5] for explicit formulae of WF, WA, and WL. Using Little’s law, we can obtain mean sojourn time as WA(LL)=Ls/λg¯, where Ls is the mean stationary system-length obtained in [5] and one may note that WA(LL)=WAe while performing numerical computation. Also, one may simplify (3) as stated below. First combine the roots γin(1≤i≤mr^) and [ϕ1(s)]n, yielding a geometric series, which are then summed over n. As a consequence of the formula regarding geometric series, this sum is [Im-γiϕ1(s)]-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).(6)WF∗s=∑n=0∞π-nϕ1snϕ1s=∑n=0∞∑i=1mr^ki1γin,∑i=1mr^ki2γin,…,∑i=1mr^kimγinϕ1snϕ1s=k11,k12,…,k1madjIm-γ1ϕ1sdet⁡Im-γ1ϕ1sϕ1s+k21,k22,…,k2madjIm-γ2ϕ1sdet⁡Im-γ2ϕ1sϕ1s+⋯+kmr^1,kmr^2,…,kmr^madjIm-γmr^ϕ1sdet⁡Im-γmr^ϕ1sϕ1s, with R(s)≥0. In the following, we develop a numerical scheme for inverting the above LSTs component-wise. Let us first look at the r-th component of the LST given by formula (3) which is further simplified as (6). The right-hand side of (3) involves π-∗ϕ1sϕ1(s) which is further simplified as the right-hand side of (6). Both LSTs (3) and (6) are analytic in the region R(s)≥0. Therefore (as ϕ1(s) is the LST of the transition probability matrix of the underlying Markov chain of C-MSP), the r-th (1≤r≤m) component of π-∗(ϕ1(s))ϕ1(s) may be written as a rational function in s, which is evident from the right-hand side of (6). Now, taking into account the zeros of the denominator of each mr^ terms in the right-hand side of (6), we obtain mr^ characteristic equations; namely, (7)det⁡Im-γjϕ1s=0,1≤j≤mr^.Also, each element of the matrix ϕ1(s)=(sIm-L0)-1L1=adj[sIm-L0]L1/det[sIm-L0] has denominator det[sIm-L0] which has m zeros in the region R(s)<0 (see similar proof of Theorem 7.3 of [5]). As a result, det[Im-γjϕ1(s)](1≤j≤mr^) has denominator {det[sIm-L0]}m which has m2 zeros in the region R(s)<0 with each zero having multiplicity m. Similarly, the numerator of det[Im-γjϕ1(s)](1≤j≤mr^) also has m2 zeros in the region R(s)<0 for each j(1≤j≤mr^). After canceling the common zeros of the numerator and the denominator of det[Im-γjϕ1(s)] × det[sIm-L0](1≤j≤mr^), one is left with m zeros of the numerator of det[Im-γjϕ1(s)] × det[sIm-L0](1≤j≤mr^) in the region R(s)<0. Let us define s¯i,j(1≤i≤m,1≤j≤mr^) as the zeros of the denominator of each component of π-∗(ϕ1(s))ϕ1(s) in the region R(s)<0. Also, it may be noted here that none of these zeros s¯i,j(1≤i≤m,1≤j≤mr^) is equal to a zero of det[sIm-L0]. Therefore, one may finally note that the denominator of the rth component of the LST WF∗(s) given by formula (3) which has m2r^ zeros: s¯i,j(1≤i≤m,1≤j≤mr^) with R(s)<0. One may finally obtain the rth (1≤r≤m) component of WF∗(s) using partial fractions as follows: (8)WFr∗s=∑i=1m∑j=1mr^Cijrs-s¯ij,where Cijr(1≤i≤m,1≤j≤mr^,1≤r≤m) are constants to be evaluated and R(s)≥0. It may be remarked here that while making partial fractions, we have assumed that the roots s¯ij are distinct. To get the unknown constants, we use the required number of linear equations after substituting any positive values of s (R(s)>0) in (8). Finally, we equate the component-wise WF∗(s) from (3) with (8) under the same positive values of s. One may note that, as normalizing conditions, the component-wise values of WF∗(0) and the components of WF must be included in the set of linear equations to get the above constants. These components may be obtained by putting s=0 in (8) and we may use WFr=-WFr∗(1)(s)|{s=0}, where WFr∗(1)(s) is the first differentiation of WFr∗(s) with respect to s. After that, taking inverse Laplace transform of (8), component-wise explicit closed-form expression of probability density function (p.d.f.) wFr(x) is given by (9)wFrx=∑i=1m∑j=1mr^Cijres¯ijx,1≤i≤m,1≤j≤mr^,1≤r≤m, with x>0. Finally, the p.d.f. of waiting time of the first customer may be obtained from wF(x)=∑r=1mwFr(x). From this, we can get WFr(u)=∫0uwFr(x)dx, (1≤r≤m), and WF(u)=∫0uwF(x)dx.Similar to the procedure followed earlier, we consider the zeros of the denominator of each element of the matrix(Im-G(ϕ1(s))) in (4). As discussed earlier, the denominator of each element of ϕ1(s) has m zeros in the region R(s)<0 which may be obtained by solving det[sIm-L0]=0. This implies that, given maximum arrival batch size r^, m zeros of the denominator of each element of the matrix (Im-G(ϕ1(s))) are repeated up to r^ times. Let the m roots of det[sIm-L0]=0 be denoted by r¯i(1≤i≤m). Therefore, the denominator of each element of the matrix (Im-G(ϕ1(s))) has repeated roots r¯i(1≤i≤m) with multiplicity r^. Also, we have to investigate the zeros of the denominator of each element of the matrix (Im-ϕ1(s))-1. These zeros may be obtained by using det⁡Im-ϕ1s. It can also be shown that there will be m2 zeros of the numerator of det[Im-ϕ1(s)] in the region R(s)≤0 (since s=0 is root of the equation det[Im-ϕ1(s)]=0). Also, there are m2 zeros of the denominator (which is equal to {det[sIm-L0]}m) of det[Im-ϕ1(s)] in the region R(s)<0 and they are r¯i(1≤i≤m) with multiplicity m each. After canceling the common zeros of the numerator and the denominator of det[Im-ϕ1(s)], one is left with (m-1) zeros of the numerator of det[Im-ϕ1(s)] in the region R(s)<0. Let these (m-1) roots of the numerator of det[Im-ϕ1(s)] be denoted by h¯i(1≤i≤m-1), where R(h¯i)<0. Lastly, the roots of the denominator of each element of the matrix ϕ1(s) in the right-hand side of (4) are also r¯i(1≤i≤m). Finally, if we consider all the roots (with R(s)<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: r¯i’s (1≤i≤m) roots repeated in the region R(s)<0 with multiplicity (r^-1) occur in the denominator of each component of WA∗(s) in (4). On the other hand s¯ij(1≤i≤m,1≤j≤mr^) and h¯i(1≤i≤m-1) are the roots in the region R(s)<0 that occur in the denominator of each component of WA∗(s) in (4). Following this conclusion, after using partial fractions, we may write the r-th (1≤r≤m) component of WA∗(s) (in (4) which also has factor π-∗(ϕ1(s))) as follows: (10)WAr∗s=∑n=1r^-1∑i=1mC¯rins-r¯in+∑i=1m∑j=1mr^Dijrs-s¯ij+∑i=1m-1Hirs-h¯i,1≤r≤m, where C¯rin(1≤i,r≤m,1≤n≤r^-1), Di,jr(1≤i,r≤m,1≤i≤mr^), and Hir(1≤r≤m,1≤i≤m-1) are constants to be evaluated and R(s)≥0. It may be remarked here that while making partial fractions, we have assumed that the roots h¯i and s¯ij 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(s)≥0. The roots of the denominator with R(s)≥0 must be canceled with numerator as the LSTs (3)-(5) are analytic in the region R(s)≥0.Now to obtain the unknown constants, we use the required number of linear equations after substituting any positive values ofs in (10). Finally, we equate the components of WA∗(s) from (4) with the corresponding component of (10) under the same positive values of s. One may note that, as normalizing conditions, the component-wise values of WA∗(0) and the components of WA must be included in the set of linear equations to get the above constants. These components may be obtained by putting s=0 in (10) and we may use WAr=-WAr∗1s{s=0}, where WAr∗(1)(s) is the first differentiation of WAr∗(s) with respect to s. After that, taking inverse Laplace transform of (10), component-wise explicit closed-form expression of probability density function (p.d.f.) wAr(x) is given by (11)wArx=∑n=1r^-1∑i=1mC¯rinxn-1n-1!er¯ix+∑i=1m∑j=1mr^Dijres¯ijx+∑i=1m-1Hireh¯ix,1≤r≤m,with x>0. Finally, the p.d.f. of waiting time of an arbitrary customer may be obtained by wA(x)=∑r=1mwAr(x). 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 WAr(u)=∫0uwAr(x)dx, (1≤r≤m), and WA(u)=∫0uwA(x)dx. It may be finally noted that limu→∞WAr(u) must be equal to the r-th component of WA∗ss=0. 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 byPH-distribution, we take interbatch arrival time distribution to be of PH-type having the representation (α,T), where α and T are of dimension ν. Then using the procedure adopted by Chaudhry et al. [5], S(z) can be derived as follows: (12)Sz=Im⊗αLz⊕T-1Im⊗Teν, with L(z)⊕T=(L(z)⊗Iν)+(Im⊗T), where ⊗ and ⊕ are Kronecker product and sum, respectively. For the derivation of S(z), see Theorem 7.5 in [5]. For more discussion on the derivation of S(z), 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 aPH[X]/C-MSP/1/∞ Queueing System

In Figure1, we have computed and plotted the component-wise p.d.f. of waiting times in the system of an arbitrary customer in a PH[X]/C-MSP/1/∞ queueing system using the numerical scheme presented in Section 3. The PH-type interbatch arrival time is given by α=0.350.65 and T=-1.1832.4231.367-2.896 with λ=0.032317. The batch size distribution is given by g1=0.0,g2=0.0,g3=0.5, and g4=0.5,gi=0.0(i≥5) with mean batch size g¯=3.5. The C-MSP matrices having positive lag-1 correlation coefficient 0.271908 as (13)L0=-0.54240951900.00372790000.00000000000.0043492170-0.02298872000.00062131700.00000000000.0012426330-2.2696700720,L1=0.02050345300.00000000000.51817816600.00000000000.01739686900.00062131702.25910768800.00497053400.0043492170with stationary mean service rate μ∗=0.500000 and π¯=0.4649780.4284260.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: (14)det⁡Imzr^-g1zr^-1+g2zr^-2+g3zr^-3+gr^Sz=0,where m=3,r^=4 and gi(i≥1) are given above and S(z) may be obtained by (12); see Theorem 7.5 of [5]. The mr^=12 roots of (14) inside |z|<1 are evaluated and the corresponding kij(1≤i≤12,1≤j≤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 WA∗(s) 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 wA(x) against x. Figure 1 shows that component-wise density of waiting time behaves significantly different from that of wA(x). 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.Figure 1
x versus p.d.f.Figure 2
u versus distribution function.


### 4.2. The Waiting-Time Analysis of aPareto[X]/C-MSP/1/∞ Queueing System

In Figure3, we have repeated similar experiment as given in above figures; i.e., we plotted component-wise p.d.f.’s against time in a Pareto[X]/C-MSP/1/∞ queue, where Pareto stands for Pareto distribution whose probability density function and distribution functions are given by a(x)=αaα/(a+x)α+1, α=0.3, a=2.3, x≥0, and A(x)=1-aα/(x+a)α, x>0, respectively, with shape parameter α and scale parameter a. Batch size distribution is taken arbitrary with probability generating function G(z)=0.1z2+0.1z3+0.8z4 and mean batch size g¯=3.7. The service matrices having positive lag-1 correlation coefficient 0.018690 are taken as (15)L0=-2.1567500.0071250.000.149625-0.1710000.0071250.0071250.007125-2.443000,L1=2.0071250.142500.000.000.0071250.0071250.000.002.42875with μ∗=1.540543 and π¯=0.3437500.3281250.328125. One may note that the moments do not exist for this case of Pareto interbatch arrival time distribution. From this Pareto distribution function A(x) (through equating A(x)=0.5), we calculate median which is equal 20.882547. Here we assume λ=0.047887 so that ρ=λg¯/μ∗=0.115012. We use geometric transform approximation method (G-TAM); see Shortle et al. [11]. For this, consider N=100 probabilities yi=1-qi(1≤i≤100,0<q<1). Substituting yi in the distribution function, i.e., yi=1-qi=A(xi)=1-aα/(xi+a)α, leads to xi=aq-i/α-a. The weights are p¯i=(yi+1-yi-1)/2=(qi-1-qi+1)/2, p¯1=y1=1-(q+q2)/2, and p¯100=y100=(q99+q100)/2. A binary search on q to match ∑i=1100p¯ixi=20.882547 (which is the median of an interbatch arrival time) gives q=0.991498. Thus, the approximate Laplace transform of interbatch arrival time may be taken as (16)A∗s=∑i=1100p¯ie-sxi.The Padé approximation [4/5] of (16) is given as follows: (17)A∗s≃1.0-0.303449s+104.611371s2+16.029743s3+790.04091s41.0+20.503315s+189.254931s2+998.209224s3+3069.401520s4+4521.877843s5.Figure 3
x versus p.d.f.Now using the above procedure as discussed in previous section, we have computed the component-wise densities of waiting time for thePareto[X]/C-MSP/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 wA(x) versus x. 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.Figure 4
x versus distribution function.


### 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 GI/C-MSP/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 Table 2 of [12]. In Table 1, we have presented first seven and the last values of the distribution function. For this experiment, C-MSP service matrices with positive lag-1 correlation coefficient 0.5381025 are taken as(18)L0=-15.4490000.1160000.0520000.0000000.092000-5.0690000.0000000.0620000.0000000.213000-4.8320004.4300000.0590005.0330000.000000-5.321000,L1=15.0680000.1570000.0000000.0560000.1230000.0470004.7450000.0000000.0840000.0000000.0000000.1050000.0510000.0000000.0770000.101000with stationary mean service rate μ∗=5.388016 and π¯ = 0.2646450.2530460.2549610.227348. The interarrival time is phase-type whose representation was taken as(19)α=0.220.330.45,T=-7.4230.0003.8128.542-9.9421.0005.2100.000-8.640with λ=3.198147 which gives ρ=0.593567. It may be noted here that for an infinite-buffer GI/C-MSP/1 queue the expression for LST of waiting-time distribution of an arriving customer is the same as WF∗(s)=WA∗(s) which is given in (3) and for the sake of notational convenience we may denote this as WA∗(s).Table 1
Actual waiting-time (in the system) distribution of a customer.
Waiting-time distribution function u W 1 ( u ) W 2 ( u ) W 3 ( u ) W 4 ( u ) W ( u ) 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.5 0.497562 0.005344 0.017378 0.002648 0.522932 1.0 0.521051 0.005759 0.035850 0.003488 0.566148 2.5 0.577662 0.006778 0.082479 0.005570 0.672489 3.0 0.593261 0.007059 0.095329 0.006143 0.701792 5.0 0.642893 0.007953 0.136214 0.007968 0.795029 10.5 0.713091 0.009218 0.194041 0.010550 0.926900 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 90.5 0.752004 0.009919 0.226096 0.011981 0.999999In Figures5 and 6, we have plotted the probability densities and distribution functions, respectively, against time for the queueing model discussed in Table 1.Figure 5
x versus p.d.f.Figure 6
u versus distribution function.


## 4.1. The Waiting-Time Analysis of aPH[X]/C-MSP/1/∞ Queueing System

In Figure1, we have computed and plotted the component-wise p.d.f. of waiting times in the system of an arbitrary customer in a PH[X]/C-MSP/1/∞ queueing system using the numerical scheme presented in Section 3. The PH-type interbatch arrival time is given by α=0.350.65 and T=-1.1832.4231.367-2.896 with λ=0.032317. The batch size distribution is given by g1=0.0,g2=0.0,g3=0.5, and g4=0.5,gi=0.0(i≥5) with mean batch size g¯=3.5. The C-MSP matrices having positive lag-1 correlation coefficient 0.271908 as (13)L0=-0.54240951900.00372790000.00000000000.0043492170-0.02298872000.00062131700.00000000000.0012426330-2.2696700720,L1=0.02050345300.00000000000.51817816600.00000000000.01739686900.00062131702.25910768800.00497053400.0043492170with stationary mean service rate μ∗=0.500000 and π¯=0.4649780.4284260.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: (14)det⁡Imzr^-g1zr^-1+g2zr^-2+g3zr^-3+gr^Sz=0,where m=3,r^=4 and gi(i≥1) are given above and S(z) may be obtained by (12); see Theorem 7.5 of [5]. The mr^=12 roots of (14) inside |z|<1 are evaluated and the corresponding kij(1≤i≤12,1≤j≤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 WA∗(s) 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 wA(x) against x. Figure 1 shows that component-wise density of waiting time behaves significantly different from that of wA(x). 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.Figure 1
x versus p.d.f.Figure 2
u versus distribution function.


## 4.2. The Waiting-Time Analysis of aPareto[X]/C-MSP/1/∞ Queueing System

In Figure3, we have repeated similar experiment as given in above figures; i.e., we plotted component-wise p.d.f.’s against time in a Pareto[X]/C-MSP/1/∞ queue, where Pareto stands for Pareto distribution whose probability density function and distribution functions are given by a(x)=αaα/(a+x)α+1, α=0.3, a=2.3, x≥0, and A(x)=1-aα/(x+a)α, x>0, respectively, with shape parameter α and scale parameter a. Batch size distribution is taken arbitrary with probability generating function G(z)=0.1z2+0.1z3+0.8z4 and mean batch size g¯=3.7. The service matrices having positive lag-1 correlation coefficient 0.018690 are taken as (15)L0=-2.1567500.0071250.000.149625-0.1710000.0071250.0071250.007125-2.443000,L1=2.0071250.142500.000.000.0071250.0071250.000.002.42875with μ∗=1.540543 and π¯=0.3437500.3281250.328125. One may note that the moments do not exist for this case of Pareto interbatch arrival time distribution. From this Pareto distribution function A(x) (through equating A(x)=0.5), we calculate median which is equal 20.882547. Here we assume λ=0.047887 so that ρ=λg¯/μ∗=0.115012. We use geometric transform approximation method (G-TAM); see Shortle et al. [11]. For this, consider N=100 probabilities yi=1-qi(1≤i≤100,0<q<1). Substituting yi in the distribution function, i.e., yi=1-qi=A(xi)=1-aα/(xi+a)α, leads to xi=aq-i/α-a. The weights are p¯i=(yi+1-yi-1)/2=(qi-1-qi+1)/2, p¯1=y1=1-(q+q2)/2, and p¯100=y100=(q99+q100)/2. A binary search on q to match ∑i=1100p¯ixi=20.882547 (which is the median of an interbatch arrival time) gives q=0.991498. Thus, the approximate Laplace transform of interbatch arrival time may be taken as (16)A∗s=∑i=1100p¯ie-sxi.The Padé approximation [4/5] of (16) is given as follows: (17)A∗s≃1.0-0.303449s+104.611371s2+16.029743s3+790.04091s41.0+20.503315s+189.254931s2+998.209224s3+3069.401520s4+4521.877843s5.Figure 3
x versus p.d.f.Now using the above procedure as discussed in previous section, we have computed the component-wise densities of waiting time for thePareto[X]/C-MSP/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 wA(x) versus x. 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.Figure 4
x versus distribution function.


## 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 GI/C-MSP/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 Table 2 of [12]. In Table 1, we have presented first seven and the last values of the distribution function. For this experiment, C-MSP service matrices with positive lag-1 correlation coefficient 0.5381025 are taken as(18)L0=-15.4490000.1160000.0520000.0000000.092000-5.0690000.0000000.0620000.0000000.213000-4.8320004.4300000.0590005.0330000.000000-5.321000,L1=15.0680000.1570000.0000000.0560000.1230000.0470004.7450000.0000000.0840000.0000000.0000000.1050000.0510000.0000000.0770000.101000with stationary mean service rate μ∗=5.388016 and π¯ = 0.2646450.2530460.2549610.227348. The interarrival time is phase-type whose representation was taken as(19)α=0.220.330.45,T=-7.4230.0003.8128.542-9.9421.0005.2100.000-8.640with λ=3.198147 which gives ρ=0.593567. It may be noted here that for an infinite-buffer GI/C-MSP/1 queue the expression for LST of waiting-time distribution of an arriving customer is the same as WF∗(s)=WA∗(s) which is given in (3) and for the sake of notational convenience we may denote this as WA∗(s).Table 1
Actual waiting-time (in the system) distribution of a customer.
Waiting-time distribution function u W 1 ( u ) W 2 ( u ) W 3 ( u ) W 4 ( u ) W ( u ) 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.5 0.497562 0.005344 0.017378 0.002648 0.522932 1.0 0.521051 0.005759 0.035850 0.003488 0.566148 2.5 0.577662 0.006778 0.082479 0.005570 0.672489 3.0 0.593261 0.007059 0.095329 0.006143 0.701792 5.0 0.642893 0.007953 0.136214 0.007968 0.795029 10.5 0.713091 0.009218 0.194041 0.010550 0.926900 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 90.5 0.752004 0.009919 0.226096 0.011981 0.999999In Figures5 and 6, we have plotted the probability densities and distribution functions, respectively, against time for the queueing model discussed in Table 1.Figure 5
x versus p.d.f.Figure 6
u versus distribution function.


## 5. Conclusion

In this paper, we analyze the waiting-time distribution of theGI[X]/C-MSP/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 m2=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 GI/M/1-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 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 GI/C-MSP/1-type queue with set-up time with or without vacation(s). These problems may attract researchers’ attention in the future.



---

*Source: 7462439-2018-09-02.xml*