# Analytically Simple and Computationally Efficient Results for theGIX/Geo/c Queues

Authors: Mohan L. Chaudhry; James J. Kim; Abhijit D. Banik
Journal: Journal of Probability and Statistics
 (2019)

---

## Abstract

A simple solution to determine the distributions of queue-lengths at different observation epochs for the modelGIX/Geo/c is presented. In the past, various discrete-time queueing models, particularly the multiserver bulk-arrival queues, have been solved using complicated methods that lead to incomplete results. The purpose of this paper is to use the roots method to solve the model GIX/Geo/c that leads to a result that is analytically elegant and computationally efficient. This method works well even for the case when the inter-batch-arrival times follow heavy-tailed distributions. The roots of the underlying characteristic equation form the basis for all distributions of queue-lengths at different time epochs.


---

## Body

## 1. Introduction

The study of discrete-time queues is relatively recent compared to its continuous-time counterpart. Over the last few decades, the field quickly gained value among queueing theorists and researchers due to the digitization of information technology, particularly in the area of signal processing devices, microcomputers, and computer networks. In analyzing the discrete-time queues, many researchers have recognized that the continuous-time models are no longer suitable to accurately measure the performance of systems in which the basic operational units are digital, such as machine cycle time, bits, and packets (see Miyazawa and Takagi [1] for details). For further details on the discrete-time models and telecommunications, one may see Bruneel and Kim [2]. In this connection, see also Takagi [3].In discrete-time queueing theory, various techniques have been introduced by many researchers to analyze theGI/Geo/1 queues. Some of their techniques include the imbedded Markov chain, supplementary variable, semi-Markov process, birth and death process, matrix-analytic, combinatorial, and numerical methods. What is pervasive across these techniques is in their final product: The GI/Geo/1 queue has three distinct queue-length distributions at three different time epochs. Though each queue-length distribution has its own importance, they all have other purposes as well. For instance, the queue-length distribution at a pre-arrival time epoch is used to compute the actual waiting-time distribution, and in the case of queues with finite capacity, it is used to evaluate the blocking probabilities (see Chaudhry and Gupta [4]). The queue-length distribution at a random time epoch is needed to compute the virtual waiting-time distribution, whereas the queue-length distribution at an outside observer’s time epoch is used to obtain the mean waiting-time-in-queue using Little’s law. Moreover, relations between the queue-length distributions at different time epochs involve interesting mathematical derivations. For some of the earlier work on the discrete-time single-server queueing models, see Zukerman and Rubin [5].In comparing theGI/Geo/1 queue with its continuous-time counterpart GI/M/1, the main difference between the two models is in the measurement of time. While a continuous-time model has a time parameter that is a real number, a discrete-time model has a time parameter that is an integer. In the model GI/Geo/1, the time axis is divided into individual time slots where the duration of one time slot is a single unit of time. In each individual time slot, two events may occur: arrival and departure. When an arrival occurs before a departure, it is called GI/Geo/1 with early arrival system (EAS), and when a departure occurs before an arrival, it is called GI/Geo/1 with late arrival system (LAS). Further, in GI/Geo/1 with LAS, if a server is idle and a customer arrives, then either his service can start immediately or in the next time slot. In the former case, it is known as an immediate access (IA), whereas in the latter case, it is known as a delayed access (DA). When discussing the GI/Geo/1 type queues, there exist six queue-length distributions: pre-arrival, random, and outside observer’s time epochs (three for EAS and three for LAS-DA).The queueing modelGIX/Geo/1 considers the model GI/Geo/1 with a batch arrival. The earliest work on GIX/Geo/1 appears to be that of Vinck and Bruneel [6] who use the complex contour integration technique. Though they provide an analytical solution to GIX/Geo/1 with EAS at different time epochs, they do not provide the corresponding numerical results. Furthermore, in their work, the analysis of GIX/Geo/1 with LAS-DA is missing, yet it is deemed to be an important aspect of GIX/Geo/1 when considering applications in telecommunications (see Takagi [3]). In response to this, Chaudhry and Gupta [7] provide the first complete solution to GIX/Geo/1 using the supplementary variable technique. One of the main contribution of their work is that they do not stop after finding the probability generation function (p.g.f.) which was a common way to conclude the analysis of a queueing model at that time, but perform the inversion of the p.g.f. Doing so enables the finding of steady-state queue-length distributions in terms of the roots of the model’s characteristic equation. This technique is referred to as the roots method.Historically, the roots method was dismissed by some queueing theorists due to perceived difficulties in computing the roots of a model’s characteristic equation. Neuts states (see Stidham [8]) “in discussing matrix-analytic solutions, I had pointed out that when the Rouché roots coincide or are close together, the method of roots could be numerically inaccurate. When I finally got copies of Crommelin’s papers, I was elated to read that, for the same reasons as I, he was concerned about the clustering of roots. In 1932, Crommelin knew; in 1980, many of my colleagues did not….” Readers can refer to Neuts [9] for his other comments on the roots method. In the 1980s, commercial computing software such as MAPLE were not able to find the roots (they do now). To address the issue of root-finding in queueing theory, the QROOT software was developed by M. L. Chaudhry in the early 1990s and demonstrated that such roots can be found (see Chaudhry [10]). The roots method was then successfully adopted to solve a wide variety of queueing problems as noted by Janssen and van Leeuwaarden [11] who write “initially, the potential difficulties of root-finding were considered to be a slur on the unblemished transforms since the determination of the roots can be numerically hazardous and the roots themselves have no probabilistic interpretation. However, Chaudhry et al. [12] have made every effort to dispel the skepticism towards root-finding in queueing theory….”While the roots method is simply another way of solving queueing problems, there are added advantages to it as well. Gouweleeuw [13] states “it is more efficient to use the roots method to get explicit expressions for probabilities from generating functions.” Furthermore, a recent paper by Maity and Gupta [14] compares the spectral theory approach against the roots method. Maity and Gupta [14] identify several difficulties in getting results using the spectral theory approach, an approach which may be simpler than the matrix-geometric approach as stated in several papers by Chakka (see, e.g., Chakka [15]). As well, Daigle and Lucantoni [16] state “whenever the roots method works, it works blindingly fast.”However, while the roots method historically only dealt with queues with light-tailed distributions, more recent research by Harris et al. [17] conclude that the roots method cannot solve queues with heavy-tailed inter-arrival times. The heavy-tailed distributions constitute a class of probability distributions that are characterized by slower decay rate than the exponential or geometric distribution. When considering the heavy-tailed distributions as an inter-arrival time distribution, the consensus among some researchers is that the roots method cannot be applied due to the unique probabilistic properties of the heavy-tailed random variables. In sharing this view, Harris et al. [17] state that “the standard root-finding problem gets complicated particularly when the inter-arrival time distribution possesses a complicated non-closed form or non-analytic Laplace-Stieltjes transform (L-ST).” The same difficulty persists in discrete-time queues since the discrete heavy-tailed probability distributions such as Weibull and Log-Normal distributions do not have a closed form of probability generating function (p.g.f.). In addition, the discrete Pareto distribution, for certain values of its parameter, has an infinite mean just like the continuous Pareto distribution.Nevertheless, the heavy-tailed distributions are useful tools in modeling real-life examples (see Willinger and Paxon [18], Leland et al. [19], Park et al. [20, 21], and Pitkow [22]). In particular, the heavy-tailed distributions (or synonymously referred to as the power, long, or fat-tailed distribution) are particularly useful when modeling the inter-arrival times of network packets and connection sizes under heavy traffic congestion (see Harris et al. [17]). While the use of heavy-tailed distributions has predominantly been in the continuous-time realm, the use of discrete heavy-tailed distributions in probabilistic modeling is evident in a number of synchronous Time Division Multiplexed (TDM) communication systems such as slotted wireless communications channels, asynchronous transfer mode, and optical packet switching (see Hernández et al. [23]). In the field of medical science and biostatistics, Para and Jan [24] have recently discovered that the Burr-type XIII and Lomax distributions are helpful in modeling a discrete data which exhibits heavy tails.In this paper, we show that the modelGIX/Geo/c including the heavy-tailed inter-batch-arrival times can be solved through roots. In doing so, we express all queue-length distributions of GIX/Geo/c in terms of the roots of the model’s characteristic equation.Our analytical proof on the existence of roots (see Appendix A) and the root-finding algorithm (see Section 11) form the foundation of our contribution. In our numerical analysis, we demonstrate how the roots of the characteristic equation behave within the contour of a unit circle. In doing so, several interesting results appear as we vary the traffic intensity, the batch size distribution, inter-batch-arrival times, and the maximum batch size.While solving the modelGIX/Geo/c through roots directly addresses some statements made by several authors, there are added benefits to it as well. Queueing theory, under the scope of algorithmic efficiency, places heavy emphasis on achieving a simpler and more efficient solution procedure as it saves memory usage and computing-time. In our method, the roots are quickly found and the queue-lengths (which are all in terms of roots) can be computed efficiently. Efficiency equates to higher engineering productivity, whereas in a purely mathematical sense, a new approach to GIX/Geo/c through roots reveals something new in the literature. For instance, the real root from the characteristic equation of GIX/Geo/c can also be used to approximate the asymptotic loss probability in GIX/Geo/c/S with a moderately sized S, where S is the maximum capacity of the system. Also, the characteristic equation of GIX/Geo/c can be transformed into that of GIX/M/c. Doing so enables us to compute the roots of the characteristic equation of GIX/M/c with an inter-batch-arrival time distribution that has nonclosed and nonexistent L-STs. Lastly, it is identified that that some of the preexisting derivations in a model with light-tailed distributions are no longer applicable if the mean inter-batch-arrival time is infinite.The roots method significantly simplifies and extends the similar work done by others using different techniques: In solving the non-bulk-arrival queueGI/Geo/c, Chan and Maa [25] use the imbedded Markov chain technique to derive the queue-length distribution for EAS at the pre-arrival epoch only. In this connection, one may also see Gao et al. who use the generating function approach to solve their non-bulk queues GI/Geo/c [26, 27] and GI/D/c [28, 29]. We have avoided that approach by using the difference-equation approach. Besides, none of the above authors consider heavy-tailed inter-arrival times. Similarly, Wittevrongel et al. [30] solve using complex analysis and contour integration while only considering the EAS with no numerical results. Lastly, Chaudhry et al. [31] solve using the supplementary variable technique that consider the light-tailed inter-batch-arrival times only. The results by Chaudhry et al. [31] are in terms of iterative relations between different queue-length distributions which can be analytically laborious and numerically inefficient. Our work also responds to the comments made by several mathematicians such as Kendall, Kleinrock, and Neuts [9]. Kendall [32] made a famous remark that queueing theory wears the Laplacian curtain. Kleinrock [33] states “one of the most difficult parts of this method of spectrum factorization is to solve for the roots.”


## 2. Model Description

Consider the modelGIX/Geo/c in which customers arrive in batches of size X with a maximum size r0<r<∞. The probability mass function (p.m.f.) of X is bh=PX=h1≤h≤r with mean b¯=∑h=1rhbh and p.g.f. Bz=∑h=1rbhzhz≤1. The n-th inter-batch-arrival time, say Tnn≥1, is a discrete-time period that is measured from the moment just before the n-th batch-arrival (say at time t to the moment just before the n+1-th batch-arrival (say at time t+m). Inter-batch-arrival times are independent identically distributed random variables (i.i.d.r.v.’s) distributed as T which is divided into m time slots. The random variable T has a p.m.f. am=PT=ma0=0,m≥1, mean a¯=∑m=1∞mam=1/λ, and p.g.f. Az=∑m=1∞amzmz≤1. The batch size and the inter-batch-arrival times are independent. There are c servers in the model where each has service time Y that is independent to one another and geometrically distributed as(1)PY=y=1−μy−1μ,0<μ<1,y≥1.The mean service time of one server isEY=1/μ and the traffic intensity of the model is ρ=b¯λ/cμ<1. In addition, let ωji be the probability of an event that j customers are served out of i customers within a single time slot. The conditional probability ωji then follows the binomial distribution(2)ωji=ijμj1−μi−j,0≤j≤i<c,cjμj1−μc−j,0≤j≤c≤i,with ω00=1 and a/b=0 for a<b or a<0. Lastly, given that GIX/Geo/c is a discrete-time queueing model, it has two different but related aspects (EAS and LAS-DA). Since the model is solved under the steady-state condition, the queue-length distribution of the Late Arrival System with Immediate Access (LAS-IA) is identical to that of the EAS.


## 3.GIX/Geo/c with EAS at a Pre-Arrival Time Epoch

InGIX/Geo/c with EAS, the n-th batch arrival occurs before the departures in the n-th inter-batch arrival time. This is illustrated in Figure 1.Figure 1
GI
X/Geo/c with EAS.LetNn− be the number of customers in the system, including those, if any, receiving service just before the n-th batch arrival (i.e., at the n-th pre-arrival time epoch). The stochastic process Nn−,n≥1 forms a homogenous Markov chain:(3)Nn+1−=Nn−+Xn−Dn,Nn−,Dn≥0,Xn≥1,where Dn is the number of customer departures during Tn and Xn is the size of the n-th batch arrival. In considering the steady-state aspect of GIX/Geo/c with EAS, Nn− becomes N− as n⟶∞ with a p.m.f. Qj−=limn⟶∞PNn−=jj≥0 and Xn becomes X as n⟶∞. Let Pi,jn=PNn−=jN0−=ii,j≥0,n≥1 be the n-step transition probabilities of Nn−,n≥1. We define Pi,j=limn⟶∞Pi,jn and then using the imbedded Markov chain analysis, we have the following transition probabilities:(4)Pi,j=∑h=1rbhki+h−j,i≥0,j≥1,1−∑k=1∞Pi,k,i≥0,j=0,where bh=0 for h≤0. As n⟶∞,Dn becomes D with p.m.f. kj=limn⟶∞PDn=jj≥0, which is defined as(5)ki+h−j=0,i+h<j,∑m=1∞amcmi+h−j1−μcm−i+h−jμi+h−j,c≤j≤i+h,∑m=1∞ami+hj1−μmj1−1−μmi+h−j,1≤j≤i+h≤c,∑m=1∞∑w=1m∑s=jc∑l=c−s+1ccw−1i+h−l−sμi+h−l−s1−μcw−1−i+h−l−s×clμl1−μc−lsj1−μm−wj1−1−μm−ws−jam,1≤j<c<i+h,for 1≤h≤r. The Pi,j’s given above can be derived by replacing the single-arrival notion with a batch-arrival condition h1≤h≤r in the one-step transition probabilities of GI/Geo/c which are available in Chan and Maa [25]. Furthermore, the same results are also available in Chaudhry et al. [31]. Lastly, in expressing Qj−,j≥0 in terms of roots, the Chapman-Kolmogorov equation,(6)Qj−=∑i=0∞Qi−Pi,j,j≥0,is used extensively. Since (6) is a j set of first-order linear difference equations with Qj− being the unknown functions to be determined, we can assume a solution of a general form(7)Qj−=Czj,j≥maxc,r,where C is a nonzero constant. Note that Qj−,j≥maxc,r is the pre-arrival queue-length distribution for the overloaded servers of GIX/Geo/c with EAS since j≥maxc,r indicates that all c servers are busy throughout T. Substituting the above into (6) yields the characteristic equation of GIX/Geo/c(8)1=Bz−1Kz,where Kz=∑j=0∞kjzj=∑m=1∞amμz+1−μcm=A1−μ+μzc. Since (8) has r roots inside the unit circle z=1 (see Appendix A for proof), let these roots be z1,z2,…,zr. Hence, Qj− for j≥maxc,r becomes r-fold and is expressed as(9)Qj−=∑h=1rChzhj,j≥maxc,r,where Ch (yet to be evaluated) for 1≤h≤r are the nonzero constants. In completely finding Qj−,j≥0, we replace Qj−,j≥maxc,r in (6) with (9) such that(10)Qj−=∑i=0maxc−1,r−1Qi−Pi,j+∑i=maxc,r∞∑h=1rChzhiPi,j,j≥0.The normalizing condition(11)∑j=0maxc−1,r−1Qj−+∑j=maxc,r∞∑h=1rChzhj=1,is used in conjunction with (10) such that when c≤r we let j=1,2,…,2r−1 in (10) to have 2r equations that are required to solve for Qj−,0≤j≤r−1 and Ch,1≤h≤r. On the contrary, when c>r we let j=1,2,…,c+r−1 in (10), and with the normalizing condition we have c+r equations that are required to solve for Qj−,0≤j≤c−1 and Ch,1≤h≤r (the Ch’s when c>r are different from those when c≤r). Solving these equations leads to the solution(12)Qj−=determined  above,0≤j≤maxc−1,r−1,∑h=1rChzhj,j≥maxc,r.


## 4.GIX/Geo/c with EAS at a Random Time Epoch

In the steady-stateGIX/Geo/c with EAS, a pre-arrival time epoch falls just before a batch-arrival, whereas a random time epoch falls at any instant between two consecutive batch arrivals. Let Rn be the discrete-time period measured from t to t+i1≤i≤m. In addition, Rn are i.i.d.r.v.’s distributed as R with p.m.f. ai∗=PR=ia0∗=0,i≥1. Visual illustration of Rn and Tn for GIX/Geo/c with EAS are provided in Figure 2, where Rn is the duration of time measured from the n-th pre-arrival time epoch to the n+1-th random time epoch, whereas Tn is the duration of time measured from the n-th pre-arrival time epoch to the n+1-th pre-arrival time epoch.Figure 2
GI
X/Geo/c with EAS at pre-arrival and random time epochs.Based on discrete renewal theory and the length-biased sampling phenomenon (see Chaudhry and Templeton [34] for definition and Kim [35] for similar application in continuous-time queues), if R consists of i time slots, then its p.m.f. ai∗,1≤i≤m can be expressed in terms of am,m≥1 such that(13)ai∗=1a¯ai+ai+1+⋯+am−2+am−1+am,1≤i≤m.Asm⟶∞, the above expression becomes(14)ai∗=1a¯limm⟶∞∑l=imal,i≥1.In addition, letQj be the steady-state p.m.f. of the number of customers in the system at a random time epoch for GIX/Geo/c with EAS. Using the above relation, we can determine Qj with the following Chapman-Kolmogorov equation:(15)Qj=∑i=0∞Qi−Pi,j∗,j≥1,where Qi−,i≥0 are already known and Pi,j∗’s are the transition probabilities from Section 3, except that am is replaced with am∗. In addition, Q0=1−∑j=1∞∑i=0∞Qi−Pi,j∗ holds. The GASTA property (Geometric Arrivals See Time Averages; see Takagi [3]) holds if al=λ1−λl−1l≥1 is substituted into the relation between am and ai∗. Doing so leads to ai∗=aii≥1; hence, due to this property, Qj=Qj−j≥0 in GeoX/Geo/c.


## 5.GIX/Geo/c with EAS at an Outside Observer’s Time Epoch

By definition, an outside observer’s time epoch inGIX/Geo/c with EAS falls anywhere just after a batch-arrival and immediately before an event of departures. Let Qjo,j≥0 be the steady-state queue-length distribution of GIX/Geo/c with EAS at an outside observer’s time epoch.In Figure3, the n-th inter-batch-arrival time period is a single time slot i.e.,m=1; hence, the n+1-th random time epoch and the n+1-th pre-arrival time epoch coincide at t+1 (this can be checked by letting m=1 in Figure 2).In Figure3, if there are j customers in the system at an outside observer’s time epoch (with probability Qjo) and there are l customers in the system at a random time epoch (with probability Ql), then there must be an event of j−l customer departures (with probability ωj−lj) between an outside observer’s time epoch and a random time epoch. Based on this notion, we form the relation(16)Ql=∑j=ll+cQjoωj−lj,l≥0.Figure 3
Outside observer’s and random time epochs inGIX/Geo/c with EAS.Since (16) is an l set of first-order linear difference equations with Qjo being the unknown functions to be determined, we assume that Qjo for j≥maxc,r is a geometric sum that consists of the same roots of (8) but with different constant coefficients (hence they are unknown). Let these unknown constant terms be Eh1≤h≤r, then we have(17)Qjo=∑h=1rEhzhj,j≥maxc,r,which is the p.m.f. of the steady-state number of customers in the system for the overloaded servers of GIX/Geo/c with EAS at an outside observer’s time epoch. By replacing Qjo,j≥maxc,r in (16) with the above geometric sum, we have(18)Ql=∑j=lmaxc−1,r−1Qjoωj−lj+∑j=maxc,rl+c∑h=1rEhzhjωj−lj,l≥0,with the normalizing condition(19)∑j=0maxc−1,r−1Qjo+∑j=maxc,r∞∑h=1rEhzhj=1.Whenc≤r we let l=1,2,…,2r−1 in (18), and with (19) it leads to 2r equations that are required to solve for Qjo,0≤j≤r−1 and Eh,1≤h≤r. On the contrary, when c>r we let l=1,2,…,c+r−1 in (18), and with (19) it leads to c+r equations that are required to solve for Qjo,0≤j≤c−1 and Eh,1≤h≤r (the Eh’s when c>r are different from those when c≤r). Solving these equations leads to an expression of Qjo for j≥0 as(20)Qjo=determined  above,0≤j≤maxc−1,r−1,∑h=1rEhzhj,j≥maxc,r.


## 6.GIX/Geo/c with LAS-DA at a Pre-Arrival Time Epoch

InGIX/Geo/c with LAS-DA, the n-th batch arrival occurs after the departures in the n-th inter-batch arrival time. This is illustrated in Figure 4.Figure 4
GI
X/Geo/c with LAS-DA.LetMn− be the number of customers in the system, including the ones, if any, in service just before the n-th batch-arrival. The random variable Mn− becomes M− and has a p.m.f. Pj−=limn⟶∞PMn−=jj≥0 as we assume the steady-state condition n⟶∞. Before proceeding to find Pj−, it is worth noting that Pj− can be found independently of Qk−0≤j≤k (as done by Chaudhry and Kim [36] in solving the model GIX/Geo/1). However, we leverage Qk− from Section 3 to determine Pj− for the purpose of demonstrating continuity in our solution procedure. Let Pj−s be the p.m.f. of there being j customers in the system GIX/Geo/c with LAS-DA at time s, and let Qk−t be the p.m.f. of there being k0≤k≤j customers in the system GIX/Geo/c with EAS at time t. Figure 5 illustrates that in EAS, the n-th batch arrives at the beginning of the time slot t,t+1, whereas in LAS-DA, the same n-th batch would arrive at the end of the time slot s,s+1.Figure 5
A comparison between EAS and LAS-DA at then-th batch-arrival.If up toc customers depart during s,t, then we can develop a relation between Qk−t and Pj−s as(21)Qk−t=∑j=kn+cPj−sωj−kj,k≥0.As a remark,0≤k≤j indicates that if no customers depart during s,t, then k=j with probability ω0j, whereas any customer departures during s,t result in k<j with probability ωj−kj. By letting Qn−=limn⟶∞Qn−tn≥0 and Pj−=limn⟶∞Pj−s j≥0, so that we have the steady-state queue-length distributions, we have(22)Qk−=∑j=kn+cPj−ωj−kj,k≥0.Since (22) is a k set of first-order linear difference equations with Pj− being the unknown functions to be determined, we can make a similar assumption as the one made in Sections 3 and 5, and let(23)Pj−=∑h=1rFhzhj,j≥maxc,r+1be the p.m.f. of the steady-state number of customers in the system for the overloaded servers of GIX/Geo/c with LAS-DA at a pre-arrival time epoch. The Fh’s in (23) are the unknown nonzero constants. As a remark, (23) differs from (9) since we have j≥maxc,r+1 in (23) but j≥maxc,r in (9). This difference is due to the property of GIX/Geo/c with LAS-DA: assume that there are initially no customers in the system and then a batch of customers arrives after some inter-batch-arrival time. If the next inter-batch-arrival time is one time slot (with probability a1), then no customers access any of the c servers, whereas if the inter-batch-arrival time is at least two time slots (with probability am,m≥2), then customers start accessing the idle servers at the beginning of the second time slot. This property results in the definition of the overloaded servers of GIX/Geo/c with LAS-DA to come into effect when j≥maxc,r+1. To solve for Pj−,0≤j≤maxc−1,r and Fh′s, (22) is alternatively expressed as(24)Qk−=∑j=kmaxc−1,rPj−ωj−kj+∑j=maxc,r+1k+c∑h=1rFhzhjωj−kj,k≥0.We also consider the normalizing condition(25)∑j=0maxc−1,rPj−+∑j=maxc,r+1∞∑h=1rFhzhj=1,such that when c≤r+1 we let k=1,2,…,2r in (24), and with (25) it leads to 2r+1 equations that are required to solve for Pj−,0≤j≤r and Fh,1≤h≤r. On the contrary, when c>r+1 we let k=1,2,…,c+r−1 in (24), and with (25) it leads to c+r equations that are required to solve for Pj−,0≤j≤c−1 and Fh,1≤h≤r (the Fh’s when c>r+1 are different from those when c≤r+1). Solving these equations leads to Pj−,j≥0 as(26)Pj−=determined above,0≤j≤maxc−1,r,∑h=1rFhzhj,j≥maxc,r+1.


## 7.GIX/Geo/c with LAS-DA at Random and Outside-Observer’s Time Epochs

LetPj and Pjo be the p.m.f.’s of the number of customers in system in GIX/Geo/c with LAS-DA at random and outside observer’s time epochs, respectively. In LAS-DA, an outside observer’s time epoch falls in a time interval that begins just after a departure and immediately before a batch-arrival. Based on this notion, Pjo is the same as Pj. In addition, we have(27)Qjo=Pjo=Pj,where Qjo is available in (20).


## 8. p.m.f. of Waiting-Time-in-Queue and Performance Measures

LetTq be the amount of time an arbitrary customer of an incoming batch spends in queue upon batch-arrival and until entering service. As mentioned by Chaudhry and Templeton [34], if the position of an arbitrary customer of an incoming batch is F, then its p.m.f. is PF=f≡rf=∑h=f∞bh/b¯¯1≤f≤r. Since the system capacity of the GIX/Geo/c queues is infinite, the waiting time distribution remains the same in both LAS-DA and EAS (see Hunter [37] and Chaudhry and Gupta [7]). Let the p.m.f. of Tq be wk=PTq=k,k≥0. As stated by Chaudhry et al. [31], wk is(28)wk=∑i=1cQ^i,k=0,∑d=1ck∑j=max0,d−cd−1ck−1jμj1−μck−1−j∑l=d−jcclμl1−μc−lQ^c+d,k≥1,where Q^i=∑l=0i−1Ql−ri−li≥1.


## 9. Performance Measures

The modelGIX/Geo/c has various performance measures. Since the distributions are known, various moments can be calculated. In particular, denote Lso as the mean number of customers in the system and Lqo as the mean number of customers in queue, both at an outside observer’s time epoch. Lso and Lqo can be found in terms of Qjo such that(29)Lso=∑n=1∞nQno,Lqo=∑n=c+1∞n−cQno.The mean waiting-time in the system and the mean waiting-time-in-queue can also be found using Little’s law:(30)Ws=a¯Lsob¯,(31)Wq=a¯Lqob¯.In addition,Wq can be found separately by using the result in Section 8, which is Wq=∑k=1∞kwk. Furthermore, it can be shown that the average number of idle servers can be found using the expression(32)∑j=0c−1c−jQjo=c1−ρ.The left-hand side of (32) determines the average number of idle servers using the queue-length-distribution at an outside observer’s time epoch. The right-hand side of (32) determines the same number, except that it is independent of Qj−,j≥0 (hence independent of roots).


## 10. Approximation of Loss Probability inGIX/Geo/c/S Using the Loss Formula

In finite-buffer queues, the loss probability is defined as the probability of a customer being rejected from the system. In the case of finite-buffer queues with bulk-arrivals, an incoming batch of customers may be partially or totally rejected, depending on how close the number of customers in the system is to the finite-buffer immediately prior to a batch-arrival. Gouweleeuw [13] and Gouweleeuw and Tijms [38] approximate the loss probability of a wide range of queueing models with general input such as BMAP/G/1/S and GI/G/1/S. In particular, while discussing the discrete-time model GIX/G/c/S, Gouweleeuw [13] states that “… an enormous amount of methods to calculate or approximate the loss probability of a customer in such a queueing system is present in the literature…. An overview is presented of those methods which explore the relation between the finite-buffer queue and its infinite-buffer counterpart.” From this statement, it is understandable that Gouweleeuw’s [13] main effort is in presenting an efficient and accurate approximation of loss probabilities using the steady-state queue-length distribution of the infinite-buffer counterpart.The asymptotic loss probability is the loss probability asS becomes a large integer. In discussing the asymptotic loss probability of GIX/M/c/S, Kim and Choi [39] express the asymptotic loss probabilities of GIX/M/c/S in terms of the roots of the model’s characteristic equation. They prove that there is exactly one positive real root (say z1) that has the largest modulus among the r roots. Kim and Choi [39] go further to claim that z1 with its corresponding constant coefficient C1 can accurately approximate the asymptotic loss probability, even for a moderately sized S. By this claim, we can also accurately approximate the loss probability of GIX/Geo/c/S (both total and partial rejections). Depending on whether the model involves EAS or LAS-DA, we can easily identify z1 from (8) as well as C1 and F1 from Sections 3 and 6, respectively, and approximate the asymptotic loss probabilities of GIX/Geo/c/S.


## 11.GIX/Geo/c Involving Heavy-Tailed Inter-Batch-Arrival Times

The heavy-tailed distributions distinguish themselves from the light-tailed distributions by having a significantly slower rate of decay. InGIX/Geo/c, a slow decaying arrival-time distribution renders a very lengthy (or an infinite) mean inter-batch-arrival time that equates to a very small (or a zero) arrival rate λ. The arrival rate is directly proportional to ρ, and as ρ⟶0, the roots of the characteristic equation of GIX/Geo/c tend towards the origin (see, e.g., Table 1) which was a concern to some researchers.Table 1
Discrete probability distributions.
G
I
a
m
m
≥
1
A
z
z
≤
1
Parameter(s)
Mean
G
e
o
m
e
t
r
i
c
λ
λ
1
−
λ
m
−
1
λ
z
/
1
−
1
−
λ
z
λ
>
0
1
/
λ
P
o
i
s
s
o
n
λ
λ
m
−
1
/
m
−
1
!
e
−
λ
z
e
−
λ
1
−
z
λ
>
0
1
+
λ
W
e
i
b
u
l
l
β
β
−
m
/
β
′
−
1
1
/
β
′
−
1
∑
m
=
1
∞
β
−
m
z
m
β
′
=
∑
j
=
0
∞
β
−
j
,
β
>
0
∑
m
=
1
∞
m
β
−
m
/
β
′
−
1
S
L
N
φ
e
−
l
n
2
m
/
2
φ
∑
m
=
1
∞
e
−
ln
2
m
/
2
z
m
φ
=
∑
j
=
1
∞
e
−
ln
2
j
/
2
φ
∑
m
=
1
∞
m
e
−
l
n
2
m
/
2
P
a
r
e
t
o
M
δ
/
m
M
+
1
Does not exist
M
>
0
,
δ
=
φ
M
+
1
−
1
∞
,
M
≤
k
φ
M
/
φ
M
+
1
,
M
>
kThe heavy-tailed distributions are believed to impose another challenge on the roots method due to the nonclosed form of p.g.f.’s. While the light-tailed distributions have closed form ofAz (see Table 2), this is not the case for Weibull, Standard Log-Normal (SLN), and Pareto distributions. Moreover, directly solving (8) with an am that follows the Pareto distribution with M≤1 will consume a lengthy computing time (or not compute at all) due to the infinite series of a p.m.f. that decays at an extremely slow rate.Table 2
Binomial batch size distribution withbh=r/h−1ph−1qr−h+1p=0.45,q=0.55,r=21 and ρ=0.8.
z
h
G
I
GI
X/Geo/5
GI
X/M/5
S
L
N
W
e
i
b
u
l
l
2
P
a
r
e
t
o
4.5
,
2
I
G
2,3
z
1
0.9656
0.9786
0.8031
0.9453
z
2
0.5625
−
0.4433
i
0.0977
−
0.5157
i
0.3845
−
0.3232
i
0.5494
−
0.4076
i
z
3
0.5625
+
0.4433
i
0.0977
+
0.5157
i
0.3845
+
0.3232
i
05494
+
0.4076
i
z
4
02578
−
0.5013
i
0.3205
−
0.5482
i
0.1938
−
0.3462
i
0.2830
−
0.4770
i
z
5
0.2578
+
0.5013
i
0.3205
+
0.5482
i
0.1938
+
0.3462
i
0.2830
+
0.4770
i
z
6
0.0669
−
0.4518
i
0.6354
−
0.4463
i
0.0655
−
0.3256
i
0.0970
−
0.4537
i
z
7
0.0669
+
0.4518
i
0.6354
+
0.4463
i
0.0655
+
0.3256
i
0.0970
+
0.4537
i
z
8
−
0.0364
−
0.2412
i
−
0.3959
−
0.0553
i
−
0.0337
−
0.2760
i
−
0.0364
−
0.2412
i
z
9
−
0.0364
+
0.2412
i
−
0.3959
+
0.0553
i
−
0.0337
+
0.2760
i
−
0.0364
+
0.2412
i
z
10
−
0.0496
−
0.3859
i
−
0.3663
−
0.1622
i
−
0.2536
−
0.1213
i
−
0.3530
−
0.0487
i
z
11
−
0.0496
+
0.3859
i
−
0.3663
+
0.1622
i
−
0.2536
+
0.1213
i
−
0.3530
+
0.0487
i
z
12
−
0.1391
−
0.3290
i
−
0.3087
−
0.2578
i
−
0.2789
−
0.0417
i
−
0.3275
−
0.1428
i
z
13
−
0.1391
+
0.3290
i
−
0.3087
+
0.2578
i
−
0.2789
+
0.0417
i
−
0.3275
+
0.1428
i
z
14
−
0.2149
−
0.2670
i
−
0.2270
−
0.3348
i
−
0.1367
−
0.2398
i
−
0.2071
−
0.2973
i
z
15
−
0.2149
+
0.2670
i
−
0.2270
+
0.3348
i
−
0.1367
+
0.2390
i
−
0.2071
+
0.2973
i
z
16
−
0.3160
−
0.0388
i
−
0.1330
−
0.3877
i
−
0.2052
−
0.1893
i
−
0.2779
−
0.2275
i
z
17
−
0.3160
+
0.0388
i
−
0.1330
+
0.3877
i
−
0.2052
+
0.1893
i
−
0.2779
+
0.2275
i
z
18
−
0.3012
−
0.1166
i
−
0.0414
−
0.4447
i
−
0.0332
−
0.1820
i
−
0.1219
−
0.3499
i
z
19
−
0.3012
+
0.1166
i
−
0.0414
+
0.4447
i
−
0.0332
+
0.1820
i
−
0.1219
+
0.3499
i
z
20
−
0.2688
−
0.1941
i
−
0.0364
−
0.2412
i
−
0.0406
−
0.2424
i
−
0.0282
−
0.4006
i
z
21
−
0.2688
+
0.1941
i
−
0.0364
+
0.2412
i
−
0.0406
+
0.2424
i
−
0.0282
+
0.4006
i
N
91
144
49
81
L
C
10
−
110
10
−
520
10
−
4
10
−
4
z
=
1As a remark, inParetoM the function φ· is the Riemann zeta function. The k-th moment of ParetoM exists as ETk=φM−k+1/φM+1 if and only if 1≤k<M is true. On the contrary, the moments become infinite if k≥M.Solving theGIX/Geo/c queueing model involving the heavy-tailed inter-batch-arrival times requires a different numerical approach from the one traditionally used to deal with the light-tailed inter-batch-arrival time distributions. To accomplish this, we express Kz in (8) as(33)KNz=∑m=1Namμz+1−μcm.Modern computing software such as MAPLE can easily determine the roots of (8) when using the above expression of Kz with N being a positive integer. One must choose an adequately sized N given that the value of N is indirectly proportional to the rate of decay of am (larger N is required for am with a slower decay). To offset this balance, we have implemented a simple algorithm in MAPLE that determines N as well as the roots of (8) (Algorithm 1).Algorithm 1
(1)
with
p
l
o
t
s:
(2)
L
C
≔
10
−
400:
(3)
μ
≔
(4)
a
≔
m
⟶
(5)
b
≔
h
⟶
(6)
c
≔
(7)
r
≔
(8)
N
≔
1:
(9)
whileμc∗N∗aN>LC do
(10)
N≔N+1:
(11)
end do:
(12)
ifμc∗N∗aN≤LC then
(13)
N≔N−1:
(14)
end if:
(15)
y≔zr−expandzr∗∑h=1rbh∗z−h∑m=1Nμ∗z+1−μc∗m∗am
(16)
AllRoots≔fsolvey,z,z=−1−I..1+I,complex:
(17)
multiplecomplexplot,AllRoots,symbol=diagonalcross,symbolsize=10,style=point,labels=″Re″,″Im″,color=″Red″,cos+sinI,−π..π,color=″SteelBlue″To run the above program,μ,am,bh,c and r must be defined. In line (8) of our program the LC of 10−400 is used as the default value but it can be scaled depending on the desired degree of accuracy. In addition, the term μc∗N in lines (9) and (25) of our program were chosen as the demarcation point to determine N since it is the constant coefficient of the last term of z in KNz; hence, LC stands for the Last Coefficient. Line (44) plots all r roots that are found within the contour of a unit circle z=1.In comparing our root-finding method with the one by Harris et al. [17], we first briefly summarize their analytical foundation: Harris et al. [17] solve the continuous-time single-arrival (i.e., single root) models using the transform approximation method (TAM). TAM approximates the L-ST of a heavy-tailed inter-arrival distribution by discrete summation consisting of “N” terms such that(34)KNz=1N∑k=1Ne−zXk,where Xk,k=1,2,…,N, are chosen to cover the outcome space of the original inter-arrival random variable with KXj=j−0.5/N. Ultimately, N is manually picked to satisfy limn⟶∞KNz=Kz, where Kz is the L-ST of a heavy-tailed inter-arrival distribution such that KNz=a¯Nμ1−z for GI/M/1 and KNz=a¯Ncμ1−z for GI/M/c. Therefore, once N is known, solving the equation(35)z=KNz,results in the root that is needed to solve the single-arrival model. In explaining TAM, Harris et al. [17] place caveat on the fact that TAM results in N as large as 106. Moreover, Harris et al. [17] lay the analytical foundation of TAM based on the single-arrival models such as GI/M/1 and GI/M/c. Though their method could be extended to bulk-arrival models (i.e., multiple roots), it may not only lead to a laborious analytical foundation, but will require even larger N to approximate the L-ST of a heavy-tailed inter-batch arrival distribution.In our numerical results, we have employedN=333 to determine 1,000 roots of the characteristic equation in GIX/M/c when the inter-batch-arrival times follow a heavy-tailed distribution with an infinite mean. In addition, we also compute the roots of the characteristic equation for GIX/M/c with heavy-tailed inter-batch-arrival times (see Section 12); hence, a direct comparison between the TAM by Harris et al. [17] and our method can be made in terms of the size of N: finding the roots at a smaller N equates to higher efficiency and shorter computing time.


## 12. Deriving the Root Equation: FromGIX/Geo/c to GIX/M/c

Here we derive the root equation ofGIX/M/c. We also explain how to compute the roots of the characteristic equation of GIX/M/c involving heavy-tailed inter-batch-arrival times. Let t be the inter-batch-arrival time of GIX/M/c such that it is a nonzero real-number t∈R∗. It has a distribution function At=PT^≤tt>0 with the mean ET^ and arrival rate λ^ such that ET^=1/λ^. Similarly, let Y^ be the exponential service time of a server with the mean EY^ and exponential service rate μ^0<μ^<1 such that EY^=1/μ^. In GIX/Geo/c, T is divided into m time slots, and multiplying m by a very small, real number (say Δ) equates to a continuous inter-batch-arrival time t. Based on this notion, we have t=Δm (similarly, ET^=ΔET and EY^=ΔEY).To derive the root equation ofGIX/M/c, we consider the ki+h−j from Section 2 in a form that is prior to being conditioned on am and bh which is(36)cmi+h−jμi+h−j1−μcm−i+h−j,m≥1,1≤h≤r,i≥0,j≥1,i+h−j≥cm.Substitutingm=t/Δ and μ=Δμ^ into the above equation gives(37)ct/Δ!i+h−j!ct/Δ−i+h−j!Δμ^i+h−j1−Δμ^ct/Δ−i+h−j.To transform theki+h−j of GIX/Geo/c into the k^i+h−j of GIX/M/c, taking the limit of the above result as Δ⟶0 (so that the discrete time parameter becomes the continuous time parameter) gives(38)limΔ⟶0ct/Δ!i+h−j!ct/Δ−i+h−j!Δμ^i+h−j1−Δμ^ct/Δ−i+h−j.Multiplying and then dividing the right-hand-side bycti+h−jt>0,c≥1,i+h≥j, it becomes(39)limΔ⟶0ct/Δ!ct/Δi+h−jct/Δ−i+h−j!ctμ^i+h−ji+h−j!1−Δμ^1/Δct1−Δμ^−i+h−j.Sincei+h≥j and using the fact that lima⟶∞a/b1/ab=1/b!, the above simplifies to(40)cμ^ti+h−ji+h−j!e−cμ^t.We now condition this by the p.d.f.att>0 over 0,∞ such that(41)k^i+h−j=∫0∞cμ^ti+h−ji+h−j!e−cμ^tatdt,1≤h≤r,t>0,or(42)k^i+h−j=∫0∞cμ^ti+h−ji+h−j!e−cμ^tdAt,1≤h≤r,t>0,where At=PT^≤tt>0. The p.m.f. k^i+h−j has a p.g.f. K^z(43)K^z=∑i=j−h∞k^i+h−jzi+h−j,1≤h≤r,z≤1,which can also be expressed as K^z=∫0∞e−cμ^1−ztdAt=a˜cμ^1−z, where a˜s is the L-ST a˜s=∫0∞e−stdAt. Now substituting K^z into (8) leads to(44)1=∑h=1rbhz−h∑i=j−h∞∫0∞e−cμ^tcμ^tzi+h−ji+h−j!dAt,or(45)1=Bz−1K^z.The equation (45) coincides with the root equation provided by Chaudhry and Kim [40] who consider the GIX/M/c with light-tailed inter-batch-arrival times. To compute the roots of (45) when the inter-batch-arrival times follow a continuous heavy-tailed distribution, we modify K^z in (45) into K^Nz, where K^Nz=∑n=0N∫0∞e−cμ^tcμ^tzn/n!dAt. We replace μc∗N in lines (22) and (25) of Algorithm 1 with Re∫0∞e−cμ^tcμ^tN/N!dAt, and replace line (31) entirely with y≔zr−expandzr∗∑h=1rbh∗z−h∑n=0N∫0∞e−cμtcμtn/n!atdtzn in Section 11. Doing so enables us to easily compute the roots of (45) when the inter-batch-arrival times of GIX/M/c follow a distribution such as Pareto (Paretoα,k):(46)at=kααtα+1,t≥k>0,α>0,or inverse-Gaussian (IGα,k):(47)at=k2πt3e−kt−α2/2α2t,t,α,k>0.


## 13. Numerical Results

The numerical results ofGIX/Geo/c involving light-tailed inter-batch-arrival times are published by Chaudhry et al. [31]. Though we can recreate their results using the roots method, we instead focus on the numerical results of GIX/Geo/c as well as GIX/M/c involving heavy-tailed inter-batch-arrival times. Our numerical results are largely divided into three parts.First, we consider different input parameters (inter-batch-arrival times, maximum batch size, and traffic intensity), and illustrate how the roots behave within the bounds of a unit circle. We then deduce some numerical analysis.We then consider an extreme case where we have plotted 1,000 roots with an infinite mean inter-batch-arrival time. In doing so, we indicate how to adjust the setting in MAPLE that will otherwise not compute those 1,000 roots in default setting.Lastly, we present the queue-length and the waiting-time-in-queue distributions ofGIX/Geo/c with heavy-tailed inter-batch-arrival times. Since all queue-length distributions are in terms of roots of the model’s characteristic equation, once the roots are found the solution procedure is deemed complete. Hence, this part of our numerical results serves as a proof of concept.


### 13.1. Computing the Roots of the Characteristic Equation

In Tables1, 3, and 4, we consider the binomial, (1,10), and normalized Poisson batch size distributions, respectively. In the same tables, the traffic intensities are in a descending order, and the maximum batch size includes both even and odd numbers. We find the roots of the characteristic equation for both GIX/Geo/c and GIX/M/c in Tables 3 and 4, whereas in Table 1 we focus on GIX/Geo/c with a larger r. We have purposely used c=5 in all tables in order to isolate and study the effect of other input parameters. All computations were performed in MAPLE using the program in Section 11. All roots were found up to the tenth decimal place and rounded to four decimal places. We have also embedded the figures of plotted roots in each table for visual illustration. At the end of Table 1, we provide our analysis in bullet points.Table 3
(1,10) batch size distribution withb1=0.9,b10=0.1, and ρ=0.25.
z
h
G
I
GI
X/Geo/5
GI
X/M/5
S
L
N
W
e
i
b
u
l
l
2
P
a
r
e
t
o
4.5
,
2
I
G
2,3
z
1
0.7408
0.7820
0.6065
0.7086
z
2
0.0974
−
0.5902
i
0.5448
−
0.4547
i
0.0209
−
0.4062
i
0.4544
−
0.4278
i
z
3
0.0974
+
0.5902
i
0.5448
+
0.4547
i
0.0209
+
0.4062
i
0.4544
+
0.4278
i
z
4
0.4798
−
0.4518
i
0.1663
−
0.6423
i
0.3137
−
0.3790
i
0.1048
−
0.5568
i
z
5
0.4798
+
0.4518
i
0.1663
+
0.6423
i
0.3137
+
0.3790
i
0.1048
+
0.5568
i
z
6
−
0.4953
−
0.6026
−
0.3262
−
0.5033
z
7
−
0.4344
−
0.2646
i
−
0.5021
−
0.3460
i
−
0.1765
−
0.3123
i
−
0.2121
−
0.4860
i
z
8
−
0.4344
+
0.2646
i
−
0.5021
+
0.3460
i
−
0.1765
+
0.3123
i
−
0.2121
+
0.4860
i
z
9
−
0.2377
−
0.4915
i
−
0.2221
−
0.5901
i
−
0.2895
−
0.1665
i
−
0.4275
−
0.2780
i
z
10
−
0.2377
+
0.4915
i
−
0.2221
+
0.5901
i
−
0.2895
+
0.1665
i
−
0.4275
+
0.2780
i
N
15
49
22
34
L
C
10
−
40
10
−
240
10
−
3
10
−
3
z
=
1Table 4
Normalized Poisson batch size distribution withbh=1/∑l=1rPl/l!Ph/h!P=0.3,r=50 and ρ=0.
Pareto
X[1]/Geo/5
z
1
0.3004
z
20
0.0072
−
0.0069
i
z
2
0.0007
−
0.0146
i
z
21
0.0072
+
0.0069
i
z
3
0.0007
+
0.0146
i
z
22
0.0073
−
0.0051
i
z
4
0.0012
−
0.0158
i
z
23
0.0073
+
0.0051
i
z
5
0.0012
+
0.0158
i
z
24
0.0101
−
0.0449
i
z
6
0.0025
−
0.0138
i
z
25
0.0101
+
0.0449
i
z
7
0.0025
+
0.0138
i
z
26
−
0.0146
z
8
0.0027
−
0.0235
i
z
27
−
0.0145
−
0.0018
i
z
9
0.0027
+
0.0235
i
z
28
−
0.0145
+
0.0018
i
z
10
0.0039
−
0.0127
i
z
29
−
0.0141
−
0.0036
i
z
11
0.0039
+
0.0127
i
z
30
−
0.0141
+
0.0036
i
z
12
0.0050
−
0.0115
i
z
31
−
0.0136
−
0.0053
i
z
13
0.0050
+
0.0115
i
z
32
−
0.0136
+
0.0053
i
z
14
0.0060
−
0.0101
i
z
33
−
0.0128
−
0.0069
i
z
39
−
0.0095
−
0.0111
i
z
45
−
0.0048
−
0.0138
i
z
15
0.0060
+
0.0101
i
z
34
−
0.0128
+
0.0069
i
z
40
−
0.0095
+
0.0111
i
z
46
−
0.0048
+
0.0138
i
z
16
0.0067
−
0.0086
i
z
35
−
0.0119
−
0.0085
i
z
41
−
0.0080
−
0.0122
i
z
47
−
0.0030
−
0.0143
i
z
17
0.0067
+
0.0086
i
z
36
−
0.0119
+
0.0085
i
z
42
−
0.0080
+
0.0122
i
z
48
−
0.0030
+
0.0143
i
z
18
0.0069
−
0.0032
i
z
37
−
0.0108
−
0.0099
i
z
43
−
0.0065
−
0.0131
i
z
49
−
0.0012
−
0.0146
i
z
19
0.0069
+
0.0032
i
z
38
−
0.0108
+
0.0099
i
z
44
−
0.0065
+
0.0131
i
z
50
−
0.0012
+
0.0146
i
L
C
=
10
−
55
N
=
19From the above tables and through additional numerical tests, we have deduced the following points:Point 1: whilez1 is always a positive real root (as proved by Kim and Choi [39]), when r is an odd number, there will be r−1 imaginary roots, and when r is an even number, there is also a negative real root and r−2 imaginary roots. In addition, z1 has the largest modulus out of all r roots inside z=1.Point 2: the imaginary roots of the characteristic equations always exist in complex conjugate pairs. Ifr is an odd number, there are r−1/2 pairs of complex conjugates, whereas if r is an even number, there are r−2/2 pairs of complex conjugates.Point 3: theLC forGIX/M/c can be significantly larger than the LC for GIX/Geo/c even though the numerical accuracy of roots is preserved. This is due to the magnitude of μc∗N being much smaller than that of ∫0∞e−cμ^tcμ^tN/N!dAt when μ=μ^.Point 4: we can verify the accuracy of roots by back-substituting any of ther roots into (8) or (45). While substituting into the characteristic equation with a root that is found at a higher N leads to a value that is closer to 1, a moderately sized N will find roots that are just as effective. The accuracy of roots can be verified in another way: compute the queue-length distributions with the r roots found at a moderately sized N and then use them to compute the left hand side of (32). Since the right hand side of (32) is independent of roots, one can choose to use a higher N to make the left-hand-side of (32) match the right-hand-side up to the desired decimal place.Point 5: as seen in Table1, as a¯⟶∞, the roots tend to converge towards the origin. However, the value of μ,0<μ<1 also influences the plotting pattern of roots. Given the relation ρ=b¯/a¯cμ, having a¯=∞ leads to ρ=0. However, making μ⟶0 will have a counter effect and may withhold the roots from converging towards the origin (see later in Figure 6, where a¯=∞ and μ=10−6). On the other hand, a large μ (and c) coupled with ρ=0 would even further the clustering of roots towards the origin. In either case, roots can be found.Point 6: the rootz1 nears 1 as ρ approaches 1. This behavior can be seen by observing the decrease in z1 in the tables as we make ρ=0.8, 0.25, and 0. However, similar to Point 5, when ρ=0, making μ⟶0 will have a counter effect and result in z1⟶1 (see later in Figure 6, where a¯=∞ and μ=10−6).Point 7: the magnitude ofμ0<μ<1 is directly proportional to the size of N. While the size of N is primarily determined by the size of LC, if LC is fixed at a certain value, a smaller μ will lead to smaller N while a larger μ will lead to a larger N. This is true since given the same power, an exponent with a smaller base will always be smaller than the other exponent with a larger base. In GIX/Geo/c and GIX/M/c with heavy-tailed inter-batch-arrival times, the am and at decay at an extremely slow rate as m,t⟶∞ while μcm and e−cμ^tcμ^tn/n! decay faster with smaller values of μ and μ^, thus requiring less number of iterations N to reach LC. We have tested this concept on ParetoX[1]/Geo/5 by finding the N at various values of μ while letting LC=10−100,bh=1/∑l=1rPl/l!Ph/h!P=0.1,r=10 and ρ=0. Readers can see Table 5 for results.Point 8: different batch size distributions lead to different plotting patterns of roots. In Table4, the (1,10) batch size distribution plotted equal number of roots on each side of the imaginary axis (with the exception of ParetoX[4.5, 2]/M/5), whereas this was not the case in other tables. However, there is always an equal number of roots plotted on each side of the real axis (this is due to the complex conjugate pairing as per Point 2).Point 9: whena¯=∞ in GIX/Geo/c, it results in Q0=Q0o=P0o=P0=1. Intuitively, this can be understood as the system being empty (in steady-state) at the outside observer’s and random time epochs since there are no arrivals while remaining customers continue to get served. While this concept can be numerically verified, it can be analytically proven as follows: since a¯=∞, the relation Qj=∑i=0∞Qi−Pi,j∗j≥1 becomes Qj=0j≥1 which results in Q0=1. Another approach is by letting ρ=0 in (32) which leads to ∑j=0c−1c−jQjo=c. This indicates that the mean number of idle servers is c (i.e., an empty system). Therefore, it can be simplified to cQ0o=c or Q0o=1. When a¯=∞ in GIX/M/c, the queue-length distribution at a random epoch (say pj,j≥0) becomes p0=1. As a remark, several relations between the performance measures of GIX/M/c are derived by Yao et al. [41]. If the inter-batch-arrival times follow the heavy-tailed distributions with an infinite mean, some of their relations that involve ρ become invalid. As an example, Yao et al. [41] derive the relation(48)Wq1−=ρλ^b¯Ls−+1−c+∑k=0c−1c−k−1pk−,where Wq1− is the mean waiting-time-in-queue of the first customer within an incoming batch, λ^ is the arrival rate, and pk− is the p.m.f. of the number of customers in the system at a pre-arrival time epoch with mean Ls−. In this relation, when ρ=0, the right-hand-side becomes undefined while the left-hand-side may not. The same phenomena can be observed when ρ=0 in another relation that determines Wq− which is the waiting-time-in-queue of a random customer within an incoming batch.Figure 6
Plotting 1,000 roots.Table 5
Relation betweenμ and N.
Pareto
X[1]/Geo/5
μ
N
0.1
19
0.2
27
0.3
36
0.4
48
0.5
63
0.6
86
0.7
123
0.8
196
0.9
413
0.99
4239


### 13.2. Extreme case

In Figure6 we have plotted the roots of the characteristic equation for ParetoX[1]/Geo/5. Doing so leads to a¯=∞ which results in ρ=0. The batch size distribution is a binomial distribution p=0.45,q=0.55 with r=1,000. We have used μ=10−6 and LC=10−10,000 which led to N=333.Using our MAPLE program from Section 11, we have successfully plotted all 1,000 roots inside the unit circle withz1=0.9999959545 (a smaller LC will lead to a z1 that is closer to 1). As a remark, one may encounter the error message in MAPLE “Length of output exceeds limit of 1000000” as z1 becomes very close to 1. This can be overcome by changing the default setting in MAPLE in the following manner: go to “Tools,” “Options,” “Precision,” and change the “Limit expression length” from 1,000,000 to 90,000,000 or greater. Doing so prevents MAPLE from rounding z1 to 1 which is what is needed to compute the queue-length distributions in terms of roots. Computation of the results for the ParetoX[α, k]/M/c queue for various parameters, readers may see Figure 7.Figure 7
Plotting the roots of the characteristic equation for theParetoX[α, k]/M/7 queue with b1=0.5,b25=0.5,μ=1, and N=50.


### 13.3. Computing the Distributions ofGIX/Geo/c Involving Heavy-Tailed Distributions

We have computed the queue-lengths and the waiting-time-in-queue distributions, as well as the performance measures for the queueGIX/Geo/c. We have considered the Weibullβ and SLN distributions as the inter-batch-arrival time distributions. All computations were performed in MAPLE (Tables 6 and 7).Table 6
Inter-batch-arrival pattern isWeibullM with am=M−m/V−1, where V=∑j=0∞M−j. Parameters taken are M=2,b1=0.425,b2=0.325,b3=0.075,b4=0.05,b5=0.125,μ=0.11,andc=5. This gives b¯=2.125,a¯=13.7054,andρ=0.2819.
Weibull
X[2]/Geo/5
j
Q
j
−
Q
j
Q
j
o
P
j
−
P
j
P
j
o
m
w
k
0
0.3300
0.4172
0.4338
0.3041
0.4338
0.4338
0
0.1996
1
0.2279
0.2097
0.2033
0.2161
0.2033
0.2033
1
0.0382
2
0.1606
0.1291
0.1323
0.1623
0.1323
0.1323
2
0.0314
3
0.1002
0.0756
0.0805
0.1053
0.0805
0.0805
3
0.0258
4
0.0635
0.0458
0.0505
0.0690
0.0505
0.0505
4
0.0212
5
0.0411
0.0286
0.0347
0.0497
0.0347
0.0347
5
0.0174
6
0.0270
0.0188
0.0228
0.0327
0.0228
0.0228
6
0.0143
7
0.0176
0.0122
0.0149
0.0215
0.0149
0.0149
7
0.0117
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
Sum
0.9999
0.9999
0.9999
0.9999
0.9999
0.9999
Sum
0.9999
∑
j
=
0
c
−
1
Q
j
o
c
−
n
=
3.5905
c
1
−
ρ
=
3.5905
L
q
o
=
0.1847
W
q
=
1.1912
b
¯
W
q
/
a
¯
=
0.1847
L
s
o
=
1.5942
W
s
=
10.2820
b
¯
W
s
/
a
¯
=
1.5942Table 7
Inter-batch-arrival pattern isSLN with am=Ve−lnm2/2 with V=1/∑j=1∞e−lnj2/2 Parameters taken are b1=0.425,b2=0.325,b3=0.075,b4=0.05,b5=0.125,μ=0.2,andc=5. This gives b¯=2.125,a¯=4.6519,andρ=0.4568.
SLN
X/Geo/5
j
Q
j
−
Q
j
Q
j
o
P
j
−
P
j
P
j
o
m
w
k
0
0.2080
0.2765
0.2339
0.1630
0.2339
0.2339
0
0.4175
1
0.2188
0.2141
0.1860
0.1844
0.1860
0.1860
1
0.0805
2
0.1811
0.1646
0.1602
0.1737
0.1602
0.1602
2
0.0591
3
0.1253
0.1107
0.1190
0.1335
0.1190
0.1190
3
0.0432
4
0.0819
0.0715
0.0837
0.0955
0.0837
0.0837
4
0.0315
5
0.0549
0.0477
0.0630
0.0725
0.0630
0.0630
5
0.0229
6
0.0396
0.0344
0.0464
0.0535
0.0464
0.0464
6
0.0167
7
0.0278
0.0241
0.0331
0.0382
0.0331
0.0331
7
0.0122
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
Sum
0.9999
0.9999
0.9999
0.9999
0.9999
0.9999
Sum
0.9999
∑
j
=
0
c
−
1
Q
j
o
c
−
n
=
2.7160
c
1
−
ρ
=
2.7160
L
q
o
=
0.5059
W
q
=
1.1075
b
¯
W
q
/
a
¯
=
0.5059
L
s
o
=
2.7899
W
s
=
6.1075
b
¯
W
s
/
a
¯
=
2.7899


## 13.1. Computing the Roots of the Characteristic Equation

In Tables1, 3, and 4, we consider the binomial, (1,10), and normalized Poisson batch size distributions, respectively. In the same tables, the traffic intensities are in a descending order, and the maximum batch size includes both even and odd numbers. We find the roots of the characteristic equation for both GIX/Geo/c and GIX/M/c in Tables 3 and 4, whereas in Table 1 we focus on GIX/Geo/c with a larger r. We have purposely used c=5 in all tables in order to isolate and study the effect of other input parameters. All computations were performed in MAPLE using the program in Section 11. All roots were found up to the tenth decimal place and rounded to four decimal places. We have also embedded the figures of plotted roots in each table for visual illustration. At the end of Table 1, we provide our analysis in bullet points.Table 3
(1,10) batch size distribution withb1=0.9,b10=0.1, and ρ=0.25.
z
h
G
I
GI
X/Geo/5
GI
X/M/5
S
L
N
W
e
i
b
u
l
l
2
P
a
r
e
t
o
4.5
,
2
I
G
2,3
z
1
0.7408
0.7820
0.6065
0.7086
z
2
0.0974
−
0.5902
i
0.5448
−
0.4547
i
0.0209
−
0.4062
i
0.4544
−
0.4278
i
z
3
0.0974
+
0.5902
i
0.5448
+
0.4547
i
0.0209
+
0.4062
i
0.4544
+
0.4278
i
z
4
0.4798
−
0.4518
i
0.1663
−
0.6423
i
0.3137
−
0.3790
i
0.1048
−
0.5568
i
z
5
0.4798
+
0.4518
i
0.1663
+
0.6423
i
0.3137
+
0.3790
i
0.1048
+
0.5568
i
z
6
−
0.4953
−
0.6026
−
0.3262
−
0.5033
z
7
−
0.4344
−
0.2646
i
−
0.5021
−
0.3460
i
−
0.1765
−
0.3123
i
−
0.2121
−
0.4860
i
z
8
−
0.4344
+
0.2646
i
−
0.5021
+
0.3460
i
−
0.1765
+
0.3123
i
−
0.2121
+
0.4860
i
z
9
−
0.2377
−
0.4915
i
−
0.2221
−
0.5901
i
−
0.2895
−
0.1665
i
−
0.4275
−
0.2780
i
z
10
−
0.2377
+
0.4915
i
−
0.2221
+
0.5901
i
−
0.2895
+
0.1665
i
−
0.4275
+
0.2780
i
N
15
49
22
34
L
C
10
−
40
10
−
240
10
−
3
10
−
3
z
=
1Table 4
Normalized Poisson batch size distribution withbh=1/∑l=1rPl/l!Ph/h!P=0.3,r=50 and ρ=0.
Pareto
X[1]/Geo/5
z
1
0.3004
z
20
0.0072
−
0.0069
i
z
2
0.0007
−
0.0146
i
z
21
0.0072
+
0.0069
i
z
3
0.0007
+
0.0146
i
z
22
0.0073
−
0.0051
i
z
4
0.0012
−
0.0158
i
z
23
0.0073
+
0.0051
i
z
5
0.0012
+
0.0158
i
z
24
0.0101
−
0.0449
i
z
6
0.0025
−
0.0138
i
z
25
0.0101
+
0.0449
i
z
7
0.0025
+
0.0138
i
z
26
−
0.0146
z
8
0.0027
−
0.0235
i
z
27
−
0.0145
−
0.0018
i
z
9
0.0027
+
0.0235
i
z
28
−
0.0145
+
0.0018
i
z
10
0.0039
−
0.0127
i
z
29
−
0.0141
−
0.0036
i
z
11
0.0039
+
0.0127
i
z
30
−
0.0141
+
0.0036
i
z
12
0.0050
−
0.0115
i
z
31
−
0.0136
−
0.0053
i
z
13
0.0050
+
0.0115
i
z
32
−
0.0136
+
0.0053
i
z
14
0.0060
−
0.0101
i
z
33
−
0.0128
−
0.0069
i
z
39
−
0.0095
−
0.0111
i
z
45
−
0.0048
−
0.0138
i
z
15
0.0060
+
0.0101
i
z
34
−
0.0128
+
0.0069
i
z
40
−
0.0095
+
0.0111
i
z
46
−
0.0048
+
0.0138
i
z
16
0.0067
−
0.0086
i
z
35
−
0.0119
−
0.0085
i
z
41
−
0.0080
−
0.0122
i
z
47
−
0.0030
−
0.0143
i
z
17
0.0067
+
0.0086
i
z
36
−
0.0119
+
0.0085
i
z
42
−
0.0080
+
0.0122
i
z
48
−
0.0030
+
0.0143
i
z
18
0.0069
−
0.0032
i
z
37
−
0.0108
−
0.0099
i
z
43
−
0.0065
−
0.0131
i
z
49
−
0.0012
−
0.0146
i
z
19
0.0069
+
0.0032
i
z
38
−
0.0108
+
0.0099
i
z
44
−
0.0065
+
0.0131
i
z
50
−
0.0012
+
0.0146
i
L
C
=
10
−
55
N
=
19From the above tables and through additional numerical tests, we have deduced the following points:Point 1: whilez1 is always a positive real root (as proved by Kim and Choi [39]), when r is an odd number, there will be r−1 imaginary roots, and when r is an even number, there is also a negative real root and r−2 imaginary roots. In addition, z1 has the largest modulus out of all r roots inside z=1.Point 2: the imaginary roots of the characteristic equations always exist in complex conjugate pairs. Ifr is an odd number, there are r−1/2 pairs of complex conjugates, whereas if r is an even number, there are r−2/2 pairs of complex conjugates.Point 3: theLC forGIX/M/c can be significantly larger than the LC for GIX/Geo/c even though the numerical accuracy of roots is preserved. This is due to the magnitude of μc∗N being much smaller than that of ∫0∞e−cμ^tcμ^tN/N!dAt when μ=μ^.Point 4: we can verify the accuracy of roots by back-substituting any of ther roots into (8) or (45). While substituting into the characteristic equation with a root that is found at a higher N leads to a value that is closer to 1, a moderately sized N will find roots that are just as effective. The accuracy of roots can be verified in another way: compute the queue-length distributions with the r roots found at a moderately sized N and then use them to compute the left hand side of (32). Since the right hand side of (32) is independent of roots, one can choose to use a higher N to make the left-hand-side of (32) match the right-hand-side up to the desired decimal place.Point 5: as seen in Table1, as a¯⟶∞, the roots tend to converge towards the origin. However, the value of μ,0<μ<1 also influences the plotting pattern of roots. Given the relation ρ=b¯/a¯cμ, having a¯=∞ leads to ρ=0. However, making μ⟶0 will have a counter effect and may withhold the roots from converging towards the origin (see later in Figure 6, where a¯=∞ and μ=10−6). On the other hand, a large μ (and c) coupled with ρ=0 would even further the clustering of roots towards the origin. In either case, roots can be found.Point 6: the rootz1 nears 1 as ρ approaches 1. This behavior can be seen by observing the decrease in z1 in the tables as we make ρ=0.8, 0.25, and 0. However, similar to Point 5, when ρ=0, making μ⟶0 will have a counter effect and result in z1⟶1 (see later in Figure 6, where a¯=∞ and μ=10−6).Point 7: the magnitude ofμ0<μ<1 is directly proportional to the size of N. While the size of N is primarily determined by the size of LC, if LC is fixed at a certain value, a smaller μ will lead to smaller N while a larger μ will lead to a larger N. This is true since given the same power, an exponent with a smaller base will always be smaller than the other exponent with a larger base. In GIX/Geo/c and GIX/M/c with heavy-tailed inter-batch-arrival times, the am and at decay at an extremely slow rate as m,t⟶∞ while μcm and e−cμ^tcμ^tn/n! decay faster with smaller values of μ and μ^, thus requiring less number of iterations N to reach LC. We have tested this concept on ParetoX[1]/Geo/5 by finding the N at various values of μ while letting LC=10−100,bh=1/∑l=1rPl/l!Ph/h!P=0.1,r=10 and ρ=0. Readers can see Table 5 for results.Point 8: different batch size distributions lead to different plotting patterns of roots. In Table4, the (1,10) batch size distribution plotted equal number of roots on each side of the imaginary axis (with the exception of ParetoX[4.5, 2]/M/5), whereas this was not the case in other tables. However, there is always an equal number of roots plotted on each side of the real axis (this is due to the complex conjugate pairing as per Point 2).Point 9: whena¯=∞ in GIX/Geo/c, it results in Q0=Q0o=P0o=P0=1. Intuitively, this can be understood as the system being empty (in steady-state) at the outside observer’s and random time epochs since there are no arrivals while remaining customers continue to get served. While this concept can be numerically verified, it can be analytically proven as follows: since a¯=∞, the relation Qj=∑i=0∞Qi−Pi,j∗j≥1 becomes Qj=0j≥1 which results in Q0=1. Another approach is by letting ρ=0 in (32) which leads to ∑j=0c−1c−jQjo=c. This indicates that the mean number of idle servers is c (i.e., an empty system). Therefore, it can be simplified to cQ0o=c or Q0o=1. When a¯=∞ in GIX/M/c, the queue-length distribution at a random epoch (say pj,j≥0) becomes p0=1. As a remark, several relations between the performance measures of GIX/M/c are derived by Yao et al. [41]. If the inter-batch-arrival times follow the heavy-tailed distributions with an infinite mean, some of their relations that involve ρ become invalid. As an example, Yao et al. [41] derive the relation(48)Wq1−=ρλ^b¯Ls−+1−c+∑k=0c−1c−k−1pk−,where Wq1− is the mean waiting-time-in-queue of the first customer within an incoming batch, λ^ is the arrival rate, and pk− is the p.m.f. of the number of customers in the system at a pre-arrival time epoch with mean Ls−. In this relation, when ρ=0, the right-hand-side becomes undefined while the left-hand-side may not. The same phenomena can be observed when ρ=0 in another relation that determines Wq− which is the waiting-time-in-queue of a random customer within an incoming batch.Figure 6
Plotting 1,000 roots.Table 5
Relation betweenμ and N.
Pareto
X[1]/Geo/5
μ
N
0.1
19
0.2
27
0.3
36
0.4
48
0.5
63
0.6
86
0.7
123
0.8
196
0.9
413
0.99
4239


## 13.2. Extreme case

In Figure6 we have plotted the roots of the characteristic equation for ParetoX[1]/Geo/5. Doing so leads to a¯=∞ which results in ρ=0. The batch size distribution is a binomial distribution p=0.45,q=0.55 with r=1,000. We have used μ=10−6 and LC=10−10,000 which led to N=333.Using our MAPLE program from Section 11, we have successfully plotted all 1,000 roots inside the unit circle withz1=0.9999959545 (a smaller LC will lead to a z1 that is closer to 1). As a remark, one may encounter the error message in MAPLE “Length of output exceeds limit of 1000000” as z1 becomes very close to 1. This can be overcome by changing the default setting in MAPLE in the following manner: go to “Tools,” “Options,” “Precision,” and change the “Limit expression length” from 1,000,000 to 90,000,000 or greater. Doing so prevents MAPLE from rounding z1 to 1 which is what is needed to compute the queue-length distributions in terms of roots. Computation of the results for the ParetoX[α, k]/M/c queue for various parameters, readers may see Figure 7.Figure 7
Plotting the roots of the characteristic equation for theParetoX[α, k]/M/7 queue with b1=0.5,b25=0.5,μ=1, and N=50.


## 13.3. Computing the Distributions ofGIX/Geo/c Involving Heavy-Tailed Distributions

We have computed the queue-lengths and the waiting-time-in-queue distributions, as well as the performance measures for the queueGIX/Geo/c. We have considered the Weibullβ and SLN distributions as the inter-batch-arrival time distributions. All computations were performed in MAPLE (Tables 6 and 7).Table 6
Inter-batch-arrival pattern isWeibullM with am=M−m/V−1, where V=∑j=0∞M−j. Parameters taken are M=2,b1=0.425,b2=0.325,b3=0.075,b4=0.05,b5=0.125,μ=0.11,andc=5. This gives b¯=2.125,a¯=13.7054,andρ=0.2819.
Weibull
X[2]/Geo/5
j
Q
j
−
Q
j
Q
j
o
P
j
−
P
j
P
j
o
m
w
k
0
0.3300
0.4172
0.4338
0.3041
0.4338
0.4338
0
0.1996
1
0.2279
0.2097
0.2033
0.2161
0.2033
0.2033
1
0.0382
2
0.1606
0.1291
0.1323
0.1623
0.1323
0.1323
2
0.0314
3
0.1002
0.0756
0.0805
0.1053
0.0805
0.0805
3
0.0258
4
0.0635
0.0458
0.0505
0.0690
0.0505
0.0505
4
0.0212
5
0.0411
0.0286
0.0347
0.0497
0.0347
0.0347
5
0.0174
6
0.0270
0.0188
0.0228
0.0327
0.0228
0.0228
6
0.0143
7
0.0176
0.0122
0.0149
0.0215
0.0149
0.0149
7
0.0117
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
Sum
0.9999
0.9999
0.9999
0.9999
0.9999
0.9999
Sum
0.9999
∑
j
=
0
c
−
1
Q
j
o
c
−
n
=
3.5905
c
1
−
ρ
=
3.5905
L
q
o
=
0.1847
W
q
=
1.1912
b
¯
W
q
/
a
¯
=
0.1847
L
s
o
=
1.5942
W
s
=
10.2820
b
¯
W
s
/
a
¯
=
1.5942Table 7
Inter-batch-arrival pattern isSLN with am=Ve−lnm2/2 with V=1/∑j=1∞e−lnj2/2 Parameters taken are b1=0.425,b2=0.325,b3=0.075,b4=0.05,b5=0.125,μ=0.2,andc=5. This gives b¯=2.125,a¯=4.6519,andρ=0.4568.
SLN
X/Geo/5
j
Q
j
−
Q
j
Q
j
o
P
j
−
P
j
P
j
o
m
w
k
0
0.2080
0.2765
0.2339
0.1630
0.2339
0.2339
0
0.4175
1
0.2188
0.2141
0.1860
0.1844
0.1860
0.1860
1
0.0805
2
0.1811
0.1646
0.1602
0.1737
0.1602
0.1602
2
0.0591
3
0.1253
0.1107
0.1190
0.1335
0.1190
0.1190
3
0.0432
4
0.0819
0.0715
0.0837
0.0955
0.0837
0.0837
4
0.0315
5
0.0549
0.0477
0.0630
0.0725
0.0630
0.0630
5
0.0229
6
0.0396
0.0344
0.0464
0.0535
0.0464
0.0464
6
0.0167
7
0.0278
0.0241
0.0331
0.0382
0.0331
0.0331
7
0.0122
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
⋮
Sum
0.9999
0.9999
0.9999
0.9999
0.9999
0.9999
Sum
0.9999
∑
j
=
0
c
−
1
Q
j
o
c
−
n
=
2.7160
c
1
−
ρ
=
2.7160
L
q
o
=
0.5059
W
q
=
1.1075
b
¯
W
q
/
a
¯
=
0.5059
L
s
o
=
2.7899
W
s
=
6.1075
b
¯
W
s
/
a
¯
=
2.7899


## 14. Conclusion

Our complete solution procedure to the modelGIX/Geo/c begins with solving for the queue-length distribution of GIX/Geo/c with EAS at a pre-arrival time epoch. In doing so, we express the result in terms of the roots of the characteristic equation. We then develop several relations between the queue-length distributions in order to find the queue-length distributions at other time epochs in GIX/Geo/c with EAS and LAS-DA. The significance of our solution procedure is that we have simplified the expression of all queue-length distributions in terms of roots. The roots are quickly and accurately found, even for the heavy-tailed inter-batch-arrival times, and this becomes advantageous for researchers and practitioners who are looking for an analytically simple and computationally efficient way to compute the steady-state queue-length distributions of the model GIX/Geo/c.In our numerical results, we focused on theGIX/Geo/c queue as well as the GIX/M/c queue involving heavy-tailed inter-batch-arrival times. Because the heavy-tailed distributions do not have the closed form of p.g.f.’s or L-STs and, in certain cases, have an infinite mean (i.e., nonexistent p.g.f.’s or L-STs), applying the roots method was considered to be not possible and numerically hazardous by some authors. However, we have demonstrated that the roots method is deemed effective and efficient even if the inter-batch-arrival times follow the heavy-tailed distributions.



---

*Source: 6480139-2019-09-03.xml*