Abstract
A simple solution to determine the distributions of queuelengths at different observation epochs for the model GI^{X}/Geo/c is presented. In the past, various discretetime queueing models, particularly the multiserver bulkarrival 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 GI^{X}/Geo/c that leads to a result that is analytically elegant and computationally efficient. This method works well even for the case when the interbatcharrival times follow heavytailed distributions. The roots of the underlying characteristic equation form the basis for all distributions of queuelengths at different time epochs.
1. Introduction
The study of discretetime queues is relatively recent compared to its continuoustime 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 discretetime queues, many researchers have recognized that the continuoustime 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 discretetime models and telecommunications, one may see Bruneel and Kim [2]. In this connection, see also Takagi [3].
In discretetime queueing theory, various techniques have been introduced by many researchers to analyze the GI/Geo/1 queues. Some of their techniques include the imbedded Markov chain, supplementary variable, semiMarkov process, birth and death process, matrixanalytic, combinatorial, and numerical methods. What is pervasive across these techniques is in their final product: The GI/Geo/1 queue has three distinct queuelength distributions at three different time epochs. Though each queuelength distribution has its own importance, they all have other purposes as well. For instance, the queuelength distribution at a prearrival time epoch is used to compute the actual waitingtime distribution, and in the case of queues with finite capacity, it is used to evaluate the blocking probabilities (see Chaudhry and Gupta [4]). The queuelength distribution at a random time epoch is needed to compute the virtual waitingtime distribution, whereas the queuelength distribution at an outside observer’s time epoch is used to obtain the mean waitingtimeinqueue using Little’s law. Moreover, relations between the queuelength distributions at different time epochs involve interesting mathematical derivations. For some of the earlier work on the discretetime singleserver queueing models, see Zukerman and Rubin [5].
In comparing the GI/Geo/1 queue with its continuoustime counterpart GI/M/1, the main difference between the two models is in the measurement of time. While a continuoustime model has a time parameter that is a real number, a discretetime 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 queuelength distributions: prearrival, random, and outside observer’s time epochs (three for EAS and three for LASDA).
The queueing model GI^{X}/Geo/1 considers the model GI/Geo/1 with a batch arrival. The earliest work on GI^{X}/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 GI^{X}/Geo/1 with EAS at different time epochs, they do not provide the corresponding numerical results. Furthermore, in their work, the analysis of GI^{X}/Geo/1 with LASDA is missing, yet it is deemed to be an important aspect of GI^{X}/Geo/1 when considering applications in telecommunications (see Takagi [3]). In response to this, Chaudhry and Gupta [7] provide the first complete solution to GI^{X}/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 steadystate queuelength 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 matrixanalytic 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 rootfinding 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 rootfinding 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 rootfinding 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 matrixgeometric 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 lighttailed distributions, more recent research by Harris et al. [17] conclude that the roots method cannot solve queues with heavytailed interarrival times. The heavytailed distributions constitute a class of probability distributions that are characterized by slower decay rate than the exponential or geometric distribution. When considering the heavytailed distributions as an interarrival time distribution, the consensus among some researchers is that the roots method cannot be applied due to the unique probabilistic properties of the heavytailed random variables. In sharing this view, Harris et al. [17] state that “the standard rootfinding problem gets complicated particularly when the interarrival time distribution possesses a complicated nonclosed form or nonanalytic LaplaceStieltjes transform (LST).” The same difficulty persists in discretetime queues since the discrete heavytailed probability distributions such as Weibull and LogNormal 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 heavytailed distributions are useful tools in modeling reallife examples (see Willinger and Paxon [18], Leland et al. [19], Park et al. [20, 21], and Pitkow [22]). In particular, the heavytailed distributions (or synonymously referred to as the power, long, or fattailed distribution) are particularly useful when modeling the interarrival times of network packets and connection sizes under heavy traffic congestion (see Harris et al. [17]). While the use of heavytailed distributions has predominantly been in the continuoustime realm, the use of discrete heavytailed 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 Burrtype XIII and Lomax distributions are helpful in modeling a discrete data which exhibits heavy tails.
In this paper, we show that the model GI^{X}/Geo/c including the heavytailed interbatcharrival times can be solved through roots. In doing so, we express all queuelength distributions of GI^{X}/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 rootfinding 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, interbatcharrival times, and the maximum batch size.
While solving the model GI^{X}/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 computingtime. In our method, the roots are quickly found and the queuelengths (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 GI^{X}/Geo/c through roots reveals something new in the literature. For instance, the real root from the characteristic equation of GI^{X}/Geo/c can also be used to approximate the asymptotic loss probability in GI^{X}/Geo/c/S with a moderately sized , where is the maximum capacity of the system. Also, the characteristic equation of GI^{X}/Geo/c can be transformed into that of GI^{X}/M/c. Doing so enables us to compute the roots of the characteristic equation of GI^{X}/M/c with an interbatcharrival time distribution that has nonclosed and nonexistent LSTs. Lastly, it is identified that that some of the preexisting derivations in a model with lighttailed distributions are no longer applicable if the mean interbatcharrival time is infinite.
The roots method significantly simplifies and extends the similar work done by others using different techniques: In solving the nonbulkarrival queue GI/Geo/c, Chan and Maa [25] use the imbedded Markov chain technique to derive the queuelength distribution for EAS at the prearrival epoch only. In this connection, one may also see Gao et al. who use the generating function approach to solve their nonbulk queues GI/Geo/c [26, 27] and GI/D/c [28, 29]. We have avoided that approach by using the differenceequation approach. Besides, none of the above authors consider heavytailed interarrival 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 lighttailed interbatcharrival times only. The results by Chaudhry et al. [31] are in terms of iterative relations between different queuelength 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 model GI^{X}/Geo/c in which customers arrive in batches of size with a maximum size . The probability mass function (p.m.f.) of is with mean and p.g.f. . The th interbatcharrival time, say , is a discretetime period that is measured from the moment just before the nth batcharrival (say at time to the moment just before the th batcharrival (say at time ). Interbatcharrival times are independent identically distributed random variables (i.i.d.r.v.’s) distributed as which is divided into time slots. The random variable has a p.m.f. , mean , and p.g.f. . The batch size and the interbatcharrival times are independent. There are servers in the model where each has service time that is independent to one another and geometrically distributed as
The mean service time of one server is and the traffic intensity of the model is . In addition, let be the probability of an event that customers are served out of customers within a single time slot. The conditional probability then follows the binomial distributionwith and for or . Lastly, given that GI^{X}/Geo/c is a discretetime queueing model, it has two different but related aspects (EAS and LASDA). Since the model is solved under the steadystate condition, the queuelength distribution of the Late Arrival System with Immediate Access (LASIA) is identical to that of the EAS.
3. GI^{X}/Geo/c with EAS at a PreArrival Time Epoch
In GI^{X}/Geo/c with EAS, the th batch arrival occurs before the departures in the th interbatch arrival time. This is illustrated in Figure 1.
Let be the number of customers in the system, including those, if any, receiving service just before the th batch arrival (i.e., at the th prearrival time epoch). The stochastic process forms a homogenous Markov chain:where is the number of customer departures during and is the size of the th batch arrival. In considering the steadystate aspect of GI^{X}/Geo/c with EAS, becomes as with a p.m.f. and becomes as . Let be the nstep transition probabilities of . We define and then using the imbedded Markov chain analysis, we have the following transition probabilities:where for . As becomes with p.m.f. which is defined asfor The ’s given above can be derived by replacing the singlearrival notion with a batcharrival condition in the onestep 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 in terms of roots, the ChapmanKolmogorov equation,is used extensively. Since (6) is a set of firstorder linear difference equations with being the unknown functions to be determined, we can assume a solution of a general formwhere is a nonzero constant. Note that is the prearrival queuelength distribution for the overloaded servers of GI^{X}/Geo/c with EAS since indicates that all servers are busy throughout . Substituting the above into (6) yields the characteristic equation of GI^{X}/Geo/cwhere . Since (8) has roots inside the unit circle (see Appendix A for proof), let these roots be . Hence, for becomes fold and is expressed aswhere (yet to be evaluated) for are the nonzero constants. In completely finding , we replace in (6) with (9) such that
The normalizing conditionis used in conjunction with (10) such that when we let in (10) to have equations that are required to solve for and . On the contrary, when we let in (10), and with the normalizing condition we have equations that are required to solve for and (the ’s when are different from those when ). Solving these equations leads to the solution
4. GI^{X}/Geo/c with EAS at a Random Time Epoch
In the steadystate GI^{X}/Geo/c with EAS, a prearrival time epoch falls just before a batcharrival, whereas a random time epoch falls at any instant between two consecutive batch arrivals. Let be the discretetime period measured from to . In addition, are i.i.d.r.v.’s distributed as with p.m.f. . Visual illustration of and for GI^{X}/Geo/c with EAS are provided in Figure 2, where is the duration of time measured from the th prearrival time epoch to the th random time epoch, whereas is the duration of time measured from the th prearrival time epoch to the th prearrival time epoch.
Based on discrete renewal theory and the lengthbiased sampling phenomenon (see Chaudhry and Templeton [34] for definition and Kim [35] for similar application in continuoustime queues), if consists of time slots, then its p.m.f. can be expressed in terms of such that
As , the above expression becomes
In addition, let be the steadystate p.m.f. of the number of customers in the system at a random time epoch for GI^{X}/Geo/c with EAS. Using the above relation, we can determine with the following ChapmanKolmogorov equation:where are already known and ’s are the transition probabilities from Section 3, except that is replaced with . In addition, holds. The GASTA property (Geometric Arrivals See Time Averages; see Takagi [3]) holds if is substituted into the relation between and . Doing so leads to ; hence, due to this property, in Geo^{X}/Geo/c.
5. GI^{X}/Geo/c with EAS at an Outside Observer’s Time Epoch
By definition, an outside observer’s time epoch in GI^{X}/Geo/c with EAS falls anywhere just after a batcharrival and immediately before an event of departures. Let be the steadystate queuelength distribution of GI^{X}/Geo/c with EAS at an outside observer’s time epoch.
In Figure 3, the th interbatcharrival time period is a single time slot ; hence, the th random time epoch and the th prearrival time epoch coincide at (this can be checked by letting in Figure 2).
In Figure 3, if there are customers in the system at an outside observer’s time epoch (with probability ) and there are customers in the system at a random time epoch (with probability ), then there must be an event of customer departures (with probability ) between an outside observer’s time epoch and a random time epoch. Based on this notion, we form the relation
Since (16) is an set of firstorder linear difference equations with being the unknown functions to be determined, we assume that for 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 , then we havewhich is the p.m.f. of the steadystate number of customers in the system for the overloaded servers of GI^{X}/Geo/c with EAS at an outside observer’s time epoch. By replacing in (16) with the above geometric sum, we havewith the normalizing condition
When we let in (18), and with (19) it leads to equations that are required to solve for and . On the contrary, when we let in (18), and with (19) it leads to equations that are required to solve for and (the ’s when are different from those when ). Solving these equations leads to an expression of for as
6. GI^{X}/Geo/c with LASDA at a PreArrival Time Epoch
In GI^{X}/Geo/c with LASDA, the th batch arrival occurs after the departures in the th interbatch arrival time. This is illustrated in Figure 4.
Let be the number of customers in the system, including the ones, if any, in service just before the th batcharrival. The random variable becomes and has a p.m.f. as we assume the steadystate condition . Before proceeding to find , it is worth noting that can be found independently of (as done by Chaudhry and Kim [36] in solving the model GI^{X}/Geo/1). However, we leverage from Section 3 to determine for the purpose of demonstrating continuity in our solution procedure. Let be the p.m.f. of there being customers in the system GI^{X}/Geo/c with LASDA at time , and let be the p.m.f. of there being customers in the system GI^{X}/Geo/c with EAS at time . Figure 5 illustrates that in EAS, the th batch arrives at the beginning of the time slot , whereas in LASDA, the same th batch would arrive at the end of the time slot .
If up to customers depart during , then we can develop a relation between and as
As a remark, indicates that if no customers depart during , then with probability , whereas any customer departures during result in with probability . By letting and , so that we have the steadystate queuelength distributions, we have
Since (22) is a set of firstorder linear difference equations with being the unknown functions to be determined, we can make a similar assumption as the one made in Sections 3 and 5, and letbe the p.m.f. of the steadystate number of customers in the system for the overloaded servers of GI^{X}/Geo/c with LASDA at a prearrival time epoch. The ’s in (23) are the unknown nonzero constants. As a remark, (23) differs from (9) since we have in (23) but in (9). This difference is due to the property of GI^{X}/Geo/c with LASDA: assume that there are initially no customers in the system and then a batch of customers arrives after some interbatcharrival time. If the next interbatcharrival time is one time slot (with probability ), then no customers access any of the servers, whereas if the interbatcharrival time is at least two time slots (with probability ), 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 GI^{X}/Geo/c with LASDA to come into effect when . To solve for and s, (22) is alternatively expressed as
We also consider the normalizing conditionsuch that when we let in (24), and with (25) it leads to equations that are required to solve for and . On the contrary, when we let in (24), and with (25) it leads to equations that are required to solve for and (the ’s when are different from those when ). Solving these equations leads to as
7. GI^{X}/Geo/c with LASDA at Random and OutsideObserver’s Time Epochs
Let and be the p.m.f.’s of the number of customers in system in GI^{X}/Geo/c with LASDA at random and outside observer’s time epochs, respectively. In LASDA, an outside observer’s time epoch falls in a time interval that begins just after a departure and immediately before a batcharrival. Based on this notion, is the same as . In addition, we havewhere is available in (20).
8. p.m.f. of WaitingTimeinQueue and Performance Measures
Let be the amount of time an arbitrary customer of an incoming batch spends in queue upon batcharrival and until entering service. As mentioned by Chaudhry and Templeton [34], if the position of an arbitrary customer of an incoming batch is , then its p.m.f. is . Since the system capacity of the GI^{X}/Geo/c queues is infinite, the waiting time distribution remains the same in both LASDA and EAS (see Hunter [37] and Chaudhry and Gupta [7]). Let the p.m.f. of be . As stated by Chaudhry et al. [31], iswhere .
9. Performance Measures
The model GI^{X}/Geo/c has various performance measures. Since the distributions are known, various moments can be calculated. In particular, denote as the mean number of customers in the system and as the mean number of customers in queue, both at an outside observer’s time epoch. and can be found in terms of such that
The mean waitingtime in the system and the mean waitingtimeinqueue can also be found using Little’s law:
In addition, can be found separately by using the result in Section 8, which is . Furthermore, it can be shown that the average number of idle servers can be found using the expression
The lefthand side of (32) determines the average number of idle servers using the queuelengthdistribution at an outside observer’s time epoch. The righthand side of (32) determines the same number, except that it is independent of (hence independent of roots).
10. Approximation of Loss Probability in GI^{X}/Geo/c/S Using the Loss Formula
In finitebuffer queues, the loss probability is defined as the probability of a customer being rejected from the system. In the case of finitebuffer queues with bulkarrivals, 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 finitebuffer immediately prior to a batcharrival. 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 discretetime model GI^{X}/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 finitebuffer queue and its infinitebuffer 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 steadystate queuelength distribution of the infinitebuffer counterpart.
The asymptotic loss probability is the loss probability as becomes a large integer. In discussing the asymptotic loss probability of GI^{X}/M/c/S, Kim and Choi [39] express the asymptotic loss probabilities of GI^{X}/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 ) that has the largest modulus among the roots. Kim and Choi [39] go further to claim that with its corresponding constant coefficient can accurately approximate the asymptotic loss probability, even for a moderately sized . By this claim, we can also accurately approximate the loss probability of GI^{X}/Geo/c/S (both total and partial rejections). Depending on whether the model involves EAS or LASDA, we can easily identify from (8) as well as and from Sections 3 and 6, respectively, and approximate the asymptotic loss probabilities of GI^{X}/Geo/c/S.
11. GI^{X}/Geo/c Involving HeavyTailed InterBatchArrival Times
The heavytailed distributions distinguish themselves from the lighttailed distributions by having a significantly slower rate of decay. In GI^{X}/Geo/c, a slow decaying arrivaltime distribution renders a very lengthy (or an infinite) mean interbatcharrival time that equates to a very small (or a zero) arrival rate . The arrival rate is directly proportional to , and as , the roots of the characteristic equation of GI^{X}/Geo/c tend towards the origin (see, e.g., Table 1) which was a concern to some researchers.
The heavytailed distributions are believed to impose another challenge on the roots method due to the nonclosed form of p.g.f.’s. While the lighttailed distributions have closed form of (see Table 2), this is not the case for Weibull, Standard LogNormal (SLN), and Pareto distributions. Moreover, directly solving (8) with an that follows the Pareto distribution with 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.
As a remark, in the function is the Riemann zeta function The th moment of exists as if and only if is true. On the contrary, the moments become infinite if
Solving the GI^{X}/Geo/c queueing model involving the heavytailed interbatcharrival times requires a different numerical approach from the one traditionally used to deal with the lighttailed interbatcharrival time distributions. To accomplish this, we express in (8) as
Modern computing software such as MAPLE can easily determine the roots of (8) when using the above expression of with being a positive integer. One must choose an adequately sized given that the value of is indirectly proportional to the rate of decay of (larger is required for with a slower decay). To offset this balance, we have implemented a simple algorithm in MAPLE that determines as well as the roots of (8) (Algorithm 1).

To run the above program, and must be defined. In line (8) of our program the LC of is used as the default value but it can be scaled depending on the desired degree of accuracy. In addition, the term in lines (9) and (25) of our program were chosen as the demarcation point to determine since it is the constant coefficient of the last term of in ; hence, LC stands for the Last Coefficient. Line (44) plots all roots that are found within the contour of a unit circle .
In comparing our rootfinding method with the one by Harris et al. [17], we first briefly summarize their analytical foundation: Harris et al. [17] solve the continuoustime singlearrival (i.e., single root) models using the transform approximation method (TAM). TAM approximates the LST of a heavytailed interarrival distribution by discrete summation consisting of “N” terms such thatwhere , are chosen to cover the outcome space of the original interarrival random variable with . Ultimately, is manually picked to satisfy , where is the LST of a heavytailed interarrival distribution such that for GI/M/1 and for GI/M/c. Therefore, once is known, solving the equationresults in the root that is needed to solve the singlearrival model. In explaining TAM, Harris et al. [17] place caveat on the fact that TAM results in as large as . Moreover, Harris et al. [17] lay the analytical foundation of TAM based on the singlearrival models such as GI/M/1 and GI/M/c. Though their method could be extended to bulkarrival models (i.e., multiple roots), it may not only lead to a laborious analytical foundation, but will require even larger to approximate the LST of a heavytailed interbatch arrival distribution.
In our numerical results, we have employed to determine 1,000 roots of the characteristic equation in GI^{X}/M/c when the interbatcharrival times follow a heavytailed distribution with an infinite mean. In addition, we also compute the roots of the characteristic equation for GI^{X}/M/c with heavytailed interbatcharrival 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 : finding the roots at a smaller equates to higher efficiency and shorter computing time.
12. Deriving the Root Equation: From GI^{X}/Geo/c to GI^{X}/M/c
Here we derive the root equation of GI^{X}/M/c. We also explain how to compute the roots of the characteristic equation of GI^{X}/M/c involving heavytailed interbatcharrival times. Let be the interbatcharrival time of GI^{X}/M/c such that it is a nonzero realnumber . It has a distribution function with the mean and arrival rate such that . Similarly, let be the exponential service time of a server with the mean and exponential service rate such that . In GI^{X}/Geo/c, is divided into time slots, and multiplying by a very small, real number (say ) equates to a continuous interbatcharrival time . Based on this notion, we have (similarly, and ).
To derive the root equation of GI^{X}/M/c, we consider the from Section 2 in a form that is prior to being conditioned on and which is
Substituting and into the above equation gives
To transform the of GI^{X}/Geo/c into the of GI^{X}/M/c, taking the limit of the above result as (so that the discrete time parameter becomes the continuous time parameter) gives
Multiplying and then dividing the righthandside by , it becomes
Since and using the fact that , the above simplifies to
We now condition this by the p.d.f. over such thatorwhere . The p.m.f. has a p.g.f. which can also be expressed as , where is the LST . Now substituting into (8) leads toor
The equation (45) coincides with the root equation provided by Chaudhry and Kim [40] who consider the GI^{X}/M/c with lighttailed interbatcharrival times. To compute the roots of (45) when the interbatcharrival times follow a continuous heavytailed distribution, we modify in (45) into , where . We replace in lines (22) and (25) of Algorithm 1 with , and replace line (31) entirely with in Section 11. Doing so enables us to easily compute the roots of (45) when the interbatcharrival times of GI^{X}/M/c follow a distribution such as Pareto ():or inverseGaussian ():
13. Numerical Results
The numerical results of GI^{X}/Geo/c involving lighttailed interbatcharrival 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 GI^{X}/Geo/c as well as GI^{X}/M/c involving heavytailed interbatcharrival times. Our numerical results are largely divided into three parts.
First, we consider different input parameters (interbatcharrival 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 interbatcharrival 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 queuelength and the waitingtimeinqueue distributions of GI^{X}/Geo/c with heavytailed interbatcharrival times. Since all queuelength 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 Tables 1, 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 GI^{X}/Geo/c and GI^{X}/M/c in Tables 3 and 4, whereas in Table 1 we focus on GI^{X}/Geo/c with a larger . We have purposely used 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.
From the above tables and through additional numerical tests, we have deduced the following points: Point 1: while is always a positive real root (as proved by Kim and Choi [39]), when is an odd number, there will be imaginary roots, and when is an even number, there is also a negative real root and imaginary roots. In addition, has the largest modulus out of all roots inside . Point 2: the imaginary roots of the characteristic equations always exist in complex conjugate pairs. If is an odd number, there are pairs of complex conjugates, whereas if is an even number, there are pairs of complex conjugates. Point 3: the LC forGI^{X}/M/c can be significantly larger than the LC for GI^{X}/Geo/c even though the numerical accuracy of roots is preserved. This is due to the magnitude of being much smaller than that of when . Point 4: we can verify the accuracy of roots by backsubstituting any of the roots into (8) or (45). While substituting into the characteristic equation with a root that is found at a higher leads to a value that is closer to , a moderately sized will find roots that are just as effective. The accuracy of roots can be verified in another way: compute the queuelength distributions with the roots found at a moderately sized 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 to make the lefthandside of (32) match the righthandside up to the desired decimal place. Point 5: as seen in Table 1, as , the roots tend to converge towards the origin. However, the value of also influences the plotting pattern of roots. Given the relation , having leads to . However, making will have a counter effect and may withhold the roots from converging towards the origin (see later in Figure 6, where and ). On the other hand, a large (and ) coupled with would even further the clustering of roots towards the origin. In either case, roots can be found. Point 6: the root nears 1 as approaches 1. This behavior can be seen by observing the decrease in in the tables as we make , , and . However, similar to Point 5, when , making will have a counter effect and result in (see later in Figure 6, where and ). Point 7: the magnitude of is directly proportional to the size of . While the size of is primarily determined by the size of LC, if LC is fixed at a certain value, a smaller will lead to smaller while a larger will lead to a larger . 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 GI^{X}/Geo/c and GI^{X}/M/c with heavytailed interbatcharrival times, the and decay at an extremely slow rate as while and decay faster with smaller values of and , thus requiring less number of iterations to reach LC. We have tested this concept on Pareto^{X}[1]/Geo/5 by finding the at various values of while letting and . Readers can see Table 5 for results. Point 8: different batch size distributions lead to different plotting patterns of roots. In Table 4, the (1,10) batch size distribution plotted equal number of roots on each side of the imaginary axis (with the exception of Pareto^{X}[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: when in GI^{X}/Geo/c, it results in . Intuitively, this can be understood as the system being empty (in steadystate) 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 , the relation becomes which results in . Another approach is by letting in (32) which leads to . This indicates that the mean number of idle servers is (i.e., an empty system). Therefore, it can be simplified to or . When in GI^{X}/M/c, the queuelength distribution at a random epoch (say ) becomes . As a remark, several relations between the performance measures of GI^{X}/M/c are derived by Yao et al. [41]. If the interbatcharrival times follow the heavytailed distributions with an infinite mean, some of their relations that involve become invalid. As an example, Yao et al. [41] derive the relationwhere is the mean waitingtimeinqueue of the first customer within an incoming batch, is the arrival rate, and is the p.m.f. of the number of customers in the system at a prearrival time epoch with mean . In this relation, when , the righthandside becomes undefined while the lefthandside may not. The same phenomena can be observed when in another relation that determines which is the waitingtimeinqueue of a random customer within an incoming batch.
13.2. Extreme case
In Figure 6 we have plotted the roots of the characteristic equation for Pareto^{X}[1]/Geo/5. Doing so leads to which results in . The batch size distribution is a binomial distribution with We have used and which led to .
Using our MAPLE program from Section 11, we have successfully plotted all 1,000 roots inside the unit circle with (a smaller will lead to a that is closer to ). As a remark, one may encounter the error message in MAPLE “Length of output exceeds limit of 1000000” as 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 to which is what is needed to compute the queuelength distributions in terms of roots. Computation of the results for the Pareto^{X}[α, k]/M/c queue for various parameters, readers may see Figure 7.
13.3. Computing the Distributions of GI^{X}/Geo/c Involving HeavyTailed Distributions
We have computed the queuelengths and the waitingtimeinqueue distributions, as well as the performance measures for the queue GI^{X}/Geo/c. We have considered the and distributions as the interbatcharrival time distributions. All computations were performed in MAPLE (Tables 6 and 7).
14. Conclusion
Our complete solution procedure to the model GI^{X}/Geo/c begins with solving for the queuelength distribution of GI^{X}/Geo/c with EAS at a prearrival 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 queuelength distributions in order to find the queuelength distributions at other time epochs in GI^{X}/Geo/c with EAS and LASDA. The significance of our solution procedure is that we have simplified the expression of all queuelength distributions in terms of roots. The roots are quickly and accurately found, even for the heavytailed interbatcharrival times, and this becomes advantageous for researchers and practitioners who are looking for an analytically simple and computationally efficient way to compute the steadystate queuelength distributions of the model GI^{X}/Geo/c.
In our numerical results, we focused on the GI^{X}/Geo/c queue as well as the GI^{X}/M/c queue involving heavytailed interbatcharrival times. Because the heavytailed distributions do not have the closed form of p.g.f.’s or LSTs and, in certain cases, have an infinite mean (i.e., nonexistent p.g.f.’s or LSTs), 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 interbatcharrival times follow the heavytailed distributions.
Appendix
Linear difference equations have frequently been used in the theory of queues (see, e.g., Chaudhry and Templeton [34]). The application of linear difference equations in solving the GI^{X}/Geo/c queues is explained below. First, by rearranging (6), we have
Now substituting into above expression leads towhich is the underlying characteristic equation of GI^{X}/Geo/c. To prove that (8) has roots inside the unit circle , let us rewrite this equation as
Now let
Consider absolute values of and on the circle , where is positive and sufficiently small. This gives
Finally, we havewhere . Thus, for and sufficiently small, on . Since and satisfy the conditions of Rouché’s theorem, it follows that (8) has roots inside the unit circle, , since has roots inside . Extending on the work of rootfinding in queueing theory by Chaudhry and Templeton [34] and using complex analysis, each of these roots can be individually computed using the following expression:
Data Availability
The data used to support the findings of this study are included within the article.
Disclosure
This paper is based on the second author’s doctoral thesis titled “A Unified Approach to MultiServer BulkArrival Queues Using Roots” completed on 2 April 2018 at the Royal Military College of Canada (see Kim [42]).
Conflicts of Interest
The authors declare that they have no conflict of interest.
Acknowledgments
This work was supported by the Canadian Defence Applied Research Program (grant no. GRC0000B1638).