T R R ECHNICAL ESEARCH EPORT Interpolation approximations for $M|G|infty$ arrival processes by Konstantinos P. Tsoukatos, Armand M. Makowski CSHCN T.R. 99-35 (ISR T.R. 99-69) The Center for Satellite and Hybrid Communication Networks is a NASA-sponsored Commercial Space Center also supported by the Department of Defense (DOD), industry, the State of Maryland, the University of Maryland and the Institute for Systems Research. This document is a technical report in the CSHCN series originating at the University of Maryland. Web site http://www.isr.umd.edu/CSHCN/ Report Documentation Page Form Approved OMB No. 0704-0188 Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 1. REPORT DATE 3. DATES COVERED 1999 2. REPORT TYPE - 4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER Interpolation approximations for $M|G|infty$ arrival processes 5b. GRANT NUMBER 5c. PROGRAM ELEMENT NUMBER 6. AUTHOR(S) 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION Army Research Laboratory,2800 Powder Mill Road,Adelphi,MD,20783 REPORT NUMBER 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 11. SPONSOR/MONITOR’S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release; distribution unlimited 13. SUPPLEMENTARY NOTES The original document contains color images. 14. ABSTRACT see report 15. SUBJECT TERMS 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF ABSTRACT OF PAGES RESPONSIBLE PERSON a. REPORT b. ABSTRACT c. THIS PAGE 16 unclassified unclassified unclassified Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std Z39-18 Sponsored by: ARL, NASA, NSF Interpolation approximations j j1 for M G arrival processes (cid:3) y Konstantinos P. Tsoukatos Armand M. Makowski Electrical Engineering Department and Institute for Systems Research University of Maryland, College Park, MD 20742 Abstract We present an approximate analysis of a discrete{time queue with correlated ar- rival processes of the so{called MjGj1 type. The proposed heuristic approximations are developed around asymptotic results in the heavy and light tra(cid:14)c regimes. Inves- tigation of the system behavior in light tra(cid:14)c quanti(cid:12)es the di(cid:11)erences between the gradual MjGj1 inputs and the point arrivals of a classical GIjGIj1 queue. In heavy tra(cid:14)c, salient features are e(cid:11)ectively captured by the exponential distribution and the Mittag{Le(cid:15)er special function, under short{ and long{range dependence respectively. By interpolating between the heavy and light tra(cid:14)c extremes we derive approxima- tions to thequeuesize distribution, applicable to all tra(cid:14)c intensities. We examine the accuracy of these expressions and discuss possible extensions of our results in several numerical examples. 1 Introduction The conclusions of a series of measurement studies demonstrating that network tra(cid:14)c ex- hibits persistent long term correlations have spurred recent activity in the study of queueing systems with correlated arrival processes. Analytical results reveal that, when strong de- pendencies are present, diverse queueing patterns may arise, in contrast to the familiar exponential decay encountered in traditional tra(cid:14)c models with bounded exponential mo- ments. (cid:3)TheworkofthisauthorwassupportedthroughNSFGrantNSFDCDR-88-03012andtheArmyResearch Laboratory under Cooperative Agreement No. DAAL01-96-2-0002. yTheworkofthisauthorwassupportedpartiallythroughNSFGrantNSFDCDR-88-03012,NASAGrant NAGW277S and the Army Research Laboratory under Cooperative Agreement No. DAAL01-96-2-0002. 1 In this paper we deal with a discrete time queue, viewed as a surrogate for a network multiplexer, driven by an MjGj1 arrival stream. Both the discrete time MjGj1 process considered here and its continuous time variant are among the tra(cid:14)c models arising from large aggregations of on/o(cid:11) sources, that have attracted a great deal of attention. Reasons for this, such as flexibility in representing correlation functions of actual tra(cid:14)c traces, are discussed in [9]. A fluid queue fed by the continuous time version of the process has been studied at least as early as 1974 [3]. Later, Cox notices that the MjGj1 busy server process with heavy tailed G is a second order asymptotically self{similar process [4]. Here, we are interested in the entire steady{state queue length distribution at a multi- plexer fed by the MjGj1 arrival process. For arbitrary pmf G, the system lacks the desired Markovian structure and a calculation using numerical inversion techniques, e.g. [13], is not possible, since no z-transform expressions are available. An exception is provided by the ge- ometric case, where a two dimensional Markov chain formulation and a functional equation for the z-transform is given in [2]. To circumvent the di(cid:14)culties of an exact analysis, one may rely on information gleaned from various asymptotic regimes. A promising approach consists of deriving approximations from the analysis of large bu(cid:11)er asymptotics [7, 8, 10]; these estimates are exact in the limit as the bu(cid:11)er level goes to in(cid:12)nity. Our objective is to explore alternative approximations to the queue length probabilities, developed arounda combination of lightandheavy tra(cid:14)casymptotics. Such approximations become exact in the limit as the tra(cid:14)c intensity goes to zero and one respectively. In light tra(cid:14)c we take advantage of the fact that the MjGj1 arrival process is obviously \Poisson driven", sothattheReiman{Simontheory[12]applies, underaboundedexponentialmoment assumption. The resulting light tra(cid:14)c limits of the queue with MjGj1 arrivals di(cid:11)er from those of a classical GIjGIj1 queue. This is a manifestation of the fact that work that joins the system gradually, as is the case with MjGj1 inputs, generates less queueing than work that arrives instantaneously. From the heavy tra(cid:14)c regime, we collect the associated limit distribution of the queue size: This is given through the exponential function in the standard short{range dependent setup, and the Mittag{Le(cid:15)er special function in the case where the MjGj1 process is long{range dependent. The approximation to the queue size distribution is subsequently generated by interpolating between the heavy and light tra(cid:14)c extremes. For some common pmfs G the approximant assumes a simple (cid:12)nal form. More interestingly, it has the potential of capturing accurately the queue size distribution at small bu(cid:11)er sizes, for which approximations based on large bu(cid:11)er asymptotics are usually ill (cid:12)tted. On the other hand, when G has (cid:12)nite exponential moment, we do not expect the heavy{light tra(cid:14)c interpolation to be accurate for bu(cid:11)er sizes much larger than the maximum burst length: It simply does not possess the correct decay rate { it does so only as the tra(cid:14)c intensity tends to one, i.e., in the heavy tra(cid:14)c limit. Surprisingly, this drawback is often absent under long{ range dependence, since there are cases where the queue size distribution has hyperbolic asymptotics with the same exponent for all tra(cid:14)c intensities! Then an approximation is more valuable, especially when considering that, in the presence of heavy tails, alternative estimates by means of simulation take an unreasonably long time to obtain. Yet, this is 2 somewhat compromised by the unaivalability of rigorously established light tra(cid:14)c limits, under long{range dependence. In Section 6 we rely on a postulated relationship, but the problem is still unresolved. The paper is organized asfollows: The description of thesystem model is given in Section 2. Section 3 contains the main conclusions of the light and heavy tra(cid:14)c analyses. These are the ingredients for the approximation, which is presented in Section 4 and discussed through numerical examples in Section 5. Further extensions of the results are suggested in Section 6. 2 The System Model We introduce the queueing model of interest, together with the required notation. We start by presenting the MjGj1 arrival processes and several of their properties; additional facts can be found in [4, 10]. Consider a population of in(cid:12)nitely many information sources, operating in discrete{time. Sources can be in one of two states, active or idle. During time slot [n,n+1), n = 0,1,..., β new sources become active. Source j, j = 1,...,β begins generating information by n+1 n+1 thestartofslot[n+1,n+2),itsactivityperiodhasdurationσ (innumberofslots). While n+1,j active, each source emits information at a constant rate of one information unit (packet) per time slot. After its activity period expires, each source switches o(cid:11) permanently, never to generate packets again. Let b denote the number of active sources, or equivalently, the n number of packets generated by the active sources at the beginning of time slot [n,n+1). If initially (i.e., at time n = 0) there were already b active sources, we denote by σ the 0,j residual activity duration (in time slots) for the jth active source, j = 1,...,b. Throughout, the IN{valued rvs b, fβ , n = 0,1,...g, fσ , n = 1,2,...; j = 1,2,...g n+1 n,j and fσ , j = 1,2,...g satisfy the following assumptions: (i) These rvs are mutually inde- 0,j pendent; (ii) The rvs fβ , n = 0,1,...g are i.i.d. Poisson rvs with parameter λ > 0; (iii) n+1 The rvs fσ , n = 1,...; j = 1,2,...g are i.i.d. with common pmf G on f1,2,...g. Let n,j σ be a generic IN{valued rv distributed according to the pmf G, assume throughout that E[σ] < 1; (iv) The rvs fσ , j = 1,2,...g are i.i.d. IN{valued rvs distributed accord- 0,j ing to the equilibrium pmf G associated with G, i.e., if σ denotes a generic IN{valued rv e e distributed according to the pmf G , then e P[σ (cid:21) n] P[σ = n] = , n = 1,2,... (1) e E[σ] In summary, the process fb , n = 0,1,...g results from discrete{time Poisson(λ) arrivals n of information sessions, where the session duration is distributed according to the pmf G and the packet generation rate of an ongoing session is one packet per time slot. Under the 3 enforced assumptions fb , n = 0,1,...g can be identi(cid:12)ed as the busy server process of a n discrete{time MjGj1 queue; for this reason the packet arrival process fb , n = 0,1,...g is n referred to as the MjGj1 arrival process. The following proposition shows that fb , n = n 0,1,...g is a correlated process, with time dependencies controlled by the tail of σ [10]. Proposition 2.1 If b is taken to be a Poisson rv with parameter λE[σ], then the process fb , n = 0,1,...g is a (strictly) stationary ergodic process with the properties: n (a) For each n = 0,1,..., the rv b is a Poisson rv with parameter λE[σ]; n (b) Its covariance function is given by h i cov(b ,b ) = λE [σ −j]+ = λE[σ]P[σ > j], n,j = 0,1,... n+j n e (c) Its index of dispersion of counts (IDC) is given by X1 X1 λ IDC (cid:17) cov(b ,b ) = λE[σ] P[σ > j] = E[σ(σ +1)], n+j n e 2 j=0 j=0 and the process is short{range dependent (i.e., IDC (cid:12)nite) if and only if E[σ2] is (cid:12)nite. We now feed this MjGj1 arrival stream fb , n = 0,1,...g into a discrete{time single n server queue with in(cid:12)nite bu(cid:11)er capacity. Such a queueing system routinely serves as a model for a network multiplexer: If q denotes the number of packets remaining in the n multiplexer bu(cid:11)er by the end of slot [n−1,n), and the multiplexer output link can transmit c packets/slot, then the bu(cid:11)er content sequence fq , n = 0,1,...g evolves according to the n Lindley recursion q = 0; q = [q +b −c]+, n = 0,1,... (2) 0 n+1 n n+1 From Part (a) of Proposition 2.1 the average input rate to the multiplexer is E[b ] = n λE[σ], and the system is stable if the tra(cid:14)c intensity ρ (cid:17) λE[σ]/c satis(cid:12)es ρ < 1. In that case q =) q, where the IR{valued rv q is the stationary queue size in the multiplexer n n bu(cid:11)er. We are interested in evaluating P(b,ρ) (cid:17) P [q > b], b (cid:21) 0 and Q (ρ) (cid:17) E [qm], 0 (cid:20) ρ < 1, m = 1,2, ρ m ρ i.e., the probability that the stationary queue size exceeds b, and the queue size (cid:12)rst and sec- ond moments, when the tra(cid:14)c intensity is ρ. To that end we develop simple approximations that all flow from asymptotic results under heavy and light tra(cid:14)c conditions. 3 Heavy and Light Tra(cid:14)c The interpolation approximation we have in mind hinges on the availability of explicit ex- pressions for limits of system quantities as ρ ! 1 (heavy tra(cid:14)c limits), and derivatives with respect to ρ as ρ ! 0 (light tra(cid:14)c derivatives). It thus requires examination of the behavior of the queue with MjGj1 arrivals under each one of these two asymptotic regimes. 4 Light Tra(cid:14)c We start with the light tra(cid:14)c regime. The right-hand derivatives at ρ = 0 of the various metrics of interest are evaluated using the Reiman{Simon technique [12]. For the system to be in the domain of applicability of the Reiman{Simon results, an assumption on (cid:12)niteness of the exponential moment of σ is needed. h i Assumption (A) There exists θ(cid:63) > 0 such that E eθσ < 1 for θ < θ(cid:63). The detailed light tra(cid:14)c analysis of the queue with MjGj1 arrivals is provided in [15]. We summarize the conclusions in the following. Proposition 3.1 Consider the Lindley recursion (2) with integer release rate c = 1,2,..., and let b = 0,1, .... Under Assumption (A) it holds that (a) For each n = 0,1,...,c ∂n dn P(b,0+) = 0 and Q (0+) = 0, m = 1,2,.... (3) ∂ρn dρn m (b) In addition, for c = 1, E[σ]2 ∂2 P(b,0+) = Eh[σ −b]+2iP[σ > b]+2 Eh[σ −b]+i2 ∂ρ2 h i −3 E [σ −b]+ P[σ > b]+P[σ > b]2 (4) d2 E[σ2] Q (0+) = (5) dρ2 1 E[σ] and ! d2 1 E[σ2]2 Q (0+) = 1+ . (6) dρ2 2 2 E[σ]2 Proposition (3.1) delineates a light tra(cid:14)c behavior for the queue with MjGj1 arrivals that is certainly di(cid:11)erent from the one of a classical GIjGIj1 queue. As seen from (3), when c = 1, in which case the multiplexer can serve no more thanone source per time slot, the (cid:12)rst derivative of the tail probability is zero. Hence, in a Taylor expansion of P [q > b] around ρ ρ = 0 the linear term in ρ would o(cid:11)er no contribution. In contrast, the stationary workload W in a single server MjGj1 queue is known to satisfy P [W > x] (cid:24) ρ(1−G (x)) (ρ ! 0), ρ e that is, in the classical queueing setup the corresponding expansion starts with a non-zero ρ term. For MjGj1 arrivals it is the second derivative (4) which is the most informative. This highlights the role of the activity duration rv σ, through both its distribution and its (cid:12)rst two moments. Notice that even if Assumption (A) were to be relaxed, (4) shows that for P [q > b] to decay like ρ2 for small ρ it is necessary that E[σ2] be (cid:12)nite. If E[σ2] = 1, as ρ is the case for long{range dependent MjGj1 arrivals, expression (4) yields in(cid:12)nity and ρ2 is no longer the correct order of decay. A di(cid:11)erent, perhaps smaller exponent should be sought 5 in the long{range dependent case (see (21)). Finally, relations (3) reflect (though in a rough manner) the statistical multiplexing gain: Since the (cid:12)rst non-zero contribution to a the tail probability is no lower than ρc+1, (3) implies that increasing the multiplexer release rate c while maintaining the same tra(cid:14)c intensity ρ would result in a decreasing tail probability P [q > b], as could be expected. ρ Heavy tra(cid:14)c Considered next is the behavior of the queue with MjGj1 arrivals in heavy tra(cid:14)c, that is, as the arrival rate λE[σ] tends to the multiplexer release rate c from below. Clearly, as the tra(cid:14)c intensity ρ converges to one the system becomes unstable and the queue length grows unbounded. It is thus necessary to seek a suitable normalizer for the queue length process, so that its normalized version has a non trivial heavy tra(cid:14)c limit. This problem is typically addressed in a setup where the system of interest is embedded into a family of queueing systems, parametrized by an integer, say l = 1,2,..., ensuring that, as l " 1, the appropriate trend to instability is established. Such an approach was pursued in [16], providing a complete characterization of the arising heavy tra(cid:14)c limits. We tacitly assume here that the heavy tra(cid:14)c limit of the stationary distribution coincides with the stationary distribution of the heavy tra(cid:14)c limit, and gather the required results from [16] in a convenient form: Proposition 3.2 The heavy tra(cid:14)c limits of the stationary queue length distribution associ- ated with (2) can be classi(cid:12)ed as follows: (a) If E[σ2] < 1 then ! 2E[σ] limP [(1−ρ) q > x] = exp − x , x (cid:21) 0. (7) ρ!1 ρ E[σ2] (b) If P[σ > n] = n−α, n = 1,2,..., with 1 < α < 2, then ! h i (α−1)E[σ] limP (1−ρ)1/(α−1)q > x = E − xα−1 , x (cid:21) 0, (8) ρ!1 ρ α−1 Γ(2−α) where X1 xn E (x) (cid:17) , ν > 0, x 2 IR, (9) ν Γ(νn+1) n=0 is the Mittag{Le(cid:15)er special function [5]. Part (a) of Proposition 3.2 addresses the classical short{range dependent case, for which theheavytra(cid:14)cnormalizeris(1−ρ)andthelimitingheavytra(cid:14)cdistributionisexponential. Part (b) deals with a long{range dependent MjGj1 arrival process. Under long{range dependence, the heavy tra(cid:14)cqueue length distribution is expressed througha Mittag{Le(cid:15)er function with hyperbolic decay, while the power{law behavior of the heavy tra(cid:14)c normalizer is (1 − ρ)1/(α−1). We mention that this result can be stated in a more general manner to cover the situation where the tail of σ is regularly varying of order α, 1 < α < 2. 6 The results under the light and heavy tra(cid:14)c regimes are subsequently combined into approximations for all values of the tra(cid:14)c intensity. 4 Interpolation Approximations Whenever Assumption (A) is satis(cid:12)ed, P [q > b] is in(cid:12)nitely di(cid:11)erentiable with respect to ρ ρ at ρ = 0, hence it can be approximated by bringing together heavy tra(cid:14)c limits and light tra(cid:14)c derivatives into a Taylor series{like expansion. To this end we follow the approach proposed in [6]. In passing, we also discuss approximations for Q (ρ), m = 1,2. The details m are as follows: The interpolation method Consider the normalized queue length rv (1−ρ) q and de(cid:12)ne F(x,ρ) (cid:17) P [(1−ρ) q > x], 0 (cid:20) ρ < 1, x (cid:21) 0 (10) ρ and F(x,1) (cid:17) limP [(1−ρ) q > x]. (11) ρ ρ!1 Assume that partial derivatives of F(x,ρ) with respect to ρ, up to order n, at ρ = 0+, b are available. Construct F (x,ρ), the nth order interpolation approximation to F(x,ρ), by n means of the polynomial ! Xn ρi ∂i Xn 1 ∂i Fb (x,ρ) (cid:17) F(x,0+)+ F(x,1)− F(x,0+) ρn+1. (12) n i! ∂ρi i!∂ρi i=0 i=0 Observe that ∂i ∂i b b F (x,1) = F(x,1) and F (x,0+) = F(x,0+), i = 0,1,...,n, n ∂ρi n ∂ρi b that is, F (x,ρ) is precisely that unique n + 1 degree polynomial in ρ which matches the n n+1 partial derivatives of F(x,ρ) at ρ = 0+ and its heavy tra(cid:14)c limit. Now, by reversing the (1−ρ) normalization in Fb (x,ρ) we generate the nth order interpolation approximation n to P [q > b] as ρ P [q > b] (cid:25) Fb ((1−ρ) b,ρ). (13) ρ n Note that, in principle, this may lie outside [0,1], in which case it is obviously a poor approximation. To calculate the quantities associated with (13) it remains to express the partial derivatives appearing in (12) in terms of the light tra(cid:14)c derivatives of P [q > b]. We ρ have ∂ ∂ ∂ F(x,0+) = P(x,0+)+x P(x,0+) (14) ∂ρ ∂ρ ∂x 7 and ∂2 ∂2 ∂2 F(x,0+) = P(x,0+)+2x P(x,0+) ∂ρ2 ∂ρ2 ∂ρ∂x ∂ ∂2 +2x P(x,0+)+x2 P(x,0+). (15) ∂x ∂x2 In case additional light tra(cid:14)c information is available, repeated application of the chain rule will yield higher order derivatives, as needed. Approximate expressions We are now ready to write approximate expressions anchored on the heavy and light tra(cid:14)c results of Section 3. Proposition 3.2(a) provides the limit (11) that should be inserted in (12). Proposition 3.1(a) can be used to substitute for the partials in (14) and then in (12). Thus, if the multiplexer release rate is c = 1,2,..., the cth order interpolation approximation to P [q > b] is simply ρ ! 2E[σ] P [q > b] (cid:25) Fb ((1−ρ) b,ρ) = ρc+1exp − (1−ρ) b . (16) ρ c E[σ2] More can be accomplished in the case c = 1, since Proposition 3.1(b) a(cid:11)ords us a promising 2nd order interpolation approximation. Insertion of (15) in (12) yields ! 1 ∂2 2E[σ] Fb (b,ρ) = ρ2(1−ρ) P(b,0+)+ρ3exp − b (17) 2 2 ∂ρ2 E[σ2] and the latter leads to the 2nd order approximant P [q > b] (cid:25) Fb ((1−ρ) b,ρ), c = 1, (18) ρ 2 where Proposition 3.1(b) is used to supply the second partial derivative in (17). Next, we briefly deal with moment approximations. We restrict attention to the case c = 1 and consider only the queue length (cid:12)rst and second moment. The relevant light tra(cid:14)c limits are given by (3), (5) and (6). In heavy tra(cid:14)c it can be inferred from (7) that ! E[σ2] m lim(1−ρ)m Q (ρ) = m! , m = 1,2,.... m ρ!1 2E[σ] Moment interpolations are then developed in very much the same manner as distribution interpolations. We skip the details of the derivation and list the (cid:12)nal expressions E[σ2] ρ2 Q (ρ) = , c = 1 (19) 1 2E[σ] 1−ρ and ! ! ρ2 E[σ2]2 E[σ2]2 Q (ρ) (cid:25) ρ −1 + +1 , c = 1. (20) 2 4(1−ρ)2 E[σ]2 E[σ]2 8