HEAVY-TRAFFIC LIMIT OF THE GI/GI/1 STATIONARY DEPARTURE PROCESS AND ITS VARIANCE FUNCTION

Heavy-traffic limits are established for the stationary departure process from a GI/GI/1 queue and its variance function. The limit process is a function of the Brownian motion limits of the arrival and service processes plus the stationary reflected Brownian motion (RBM) limit of the queue-length process. An explicit expression is given for the variance function, which depends only on the first two moments of the interarrival times and service times plus the previously-determined correlation function of canonical (drift −1, diffusion coefficient 1) RBM. The limit for the variance function here is used to show that the approximation for the index of dispersion for counts of the departure process used in our new robust queueing network analyzer is asymptotically correct in the heavy-traffic limit.

1. Introduction.In this paper we establish a heavy-traffic (HT) limit for the stationary departure process from the stable GI/GI/1 queue and its variance function.In doing so, we are primarily motivated by our desire to develop a new robust queueing network analyzer (RQNA) for open networks of single-server queues exploiting the index of dispersion for counts (IDC) for all arrival processes, which we refer to as the RQNA-IDC.
The HT limit here is also of considerable interest more generally, because the stationary departure process from a GI/GI/1 queue is remarkably complicated; e.g., it is only a stationary renewal process in the special case of an M/M/1 model, when it is Poisson, by Burke's [13] theorem; also see [40] and references therein.Indeed, explicit transform expressions for the variance function of the stationary departure process are evidently only available for the M/GI/1 and GI/M/1 models, due to Takacs [45] and Daley [18,19]; see [10] and [31] for related results on GI/GI/1.We exploit the GI/M/1 and M/GI/1/ results here to directly establish HT limits for the variance function in §3 and §4.
The HT limit for the departure process starting empty in the GI/GI/1 model and more general multi-channel models is an old result, being contained in Theorem 2 of [33], but we evidently derive the first HT limits for the stationary departure process and its variance function for any GI/GI/1 model except M/M/1.The key to the process limit here is the recent HT limits for the stationary vector of queue lengths in a generalized Jackson network in Gamarnik and Zeevi [25] and the extension in Budhiraja and Lee [12].The existence of the HT limit for the scaled variance function is a relatively direct consequence.
In particular, Gamarnik and Zeevi [25] created a Markov process representation by appending the residual interarrival times and service times to the queue length at each node as supplementary variables.Their key step is to identify an appropriate Lyapunov function in that formidable highdimensional setting.With that approach, they establish HT limits for the sequence of steady-state distributions of the vector of queue lengths.For the relatively simple GI/GI/1 model that we consider, the HT limit in [25] only recovers a variant of the very early direct HT limit for the steady-state waiting time distribution from the early paper by Kingman [37] and its extension to general G/G/1 queues in Szczotka [44].For GI/GI/1, we exploit the HT process limit for the entire sequence of steady-state queue-length processes that follows from the Markov process representation.That HT process limit for the stationary queue length process is important to establish our associated HT process limit for the stationary departure process.(To quickly see that the process part of the HT limit is missing in [37] and [44], note that the HT space scaling is used without any time scaling.The importance of both time and space scaling is essential for HT process limits, as discussed in [50].)The extension by Budhiraja and Lee in [12] is important too, because they replace a strong moment-generating-function condition on the interarrival-time and service-time distributions in [25] by the usual 2 + ǫ moment conditions, as in Example 3.1 of [33].
An interesting and useful contribution here is the explicit form of the limiting variance function and its connection to basic functions of reflected Brownian motion (RBM) in [1,2]; see (5.20), which draws upon (3.5) and (3.9).For this step, we exploit the transform limits in the M/GI/1 and GI/M/1 special cases.We apply time-space transformations of the underlying Brownian motions in the GI/GI/1 limit to show that the limit with adjusted parameters is the same as for the M/GI/1 or GI/M/1 special case.Thus, we can identify the explicit form of the limiting variance function for the GI/GI/1 model by exploiting the results for the special cases; see the proof of Theorem 5.3.This same formula serves as an approximation more generally; see §7.
1.1.The New RQNA-IDC.We now explain a bit more about the new RQNA-IDC.It is a parametric-decomposition approximation (i.e., it treats the individual queues separately) like the queueing network analyzer (QNA) in [47].Instead of approximately characterizing each flow by its rate and a single variability parameter, we approximately characterize the flow by its rate and the IDC, which is a real-valued function on the positive halfline, in particular the scaled variance function; see §6.The need for such an enhanced version of QNA has long been recognized, as can be seen from [35,36,49,55].
This paper is a sequel to [51], which developed a robust queueing (RQ) approximation for the steady-state workload in a G/G/1 queue; see §6.1 here.That RQ algorithm extends an earlier RQ algorithm by Bandi, Bertsimas and Youssef [8].Indeed, a full RQNA is developed in [8], but the RQ approximation in [8] is a parametric RQ formulation, based on a single variability parameter, as in [47].In contrast, the new RQ algorithm in [51] involves a functional formulation that incorporates the variance of the input over time.The advantage of the new formulation is demonstrated by simulation comparisons in [51].
A full RQNA-IDC for open networks of G/GI/1 single-server queues is outlined in §6 of [51].There it is noted that the main challenge in developing an effective RQNA that can better capture dependence in the arrival processes at the different queues is to be able to approximate the IDC of the stationary departure process from a G/GI/1 queue.Important theoretical support for the new RQNA-IDC developed and studied in [53] is provided by the present paper.Corollary 6.1 here shows that the proposed approximation for the IDC of the stationary departure process is asymptotically correct for the GI/GI/1 queue in the HT limit.
1.2.On the Heavy-Traffic Scaling.We consider the standard HT scaling of both time and space (discussed in [32,33] and Chapters 5 and 9 of [50]) applied to the stationary departure process.That scaling involves the HT scaling after the usual time limit to obtain the stationary departure process.
To explain, let D ρ,k (t) denote the number of departures in the interval [0, t] for the GI/GI/1 model with traffic intensity ρ < 1, starting with k customers in the system at time 0 and fully specified remaining times until the first new arrival and the service completion for the customer in service.Then, for each s > 0, let D ρ,k,s (t) ≡ D ρ,k (s + t) − D ρ,k (s) : t ≥ 0 be that same departure process after time s.Under general regularity conditions (which we do not discuss), for each k, 0 ≤ k < ∞, and ρ, 0 < ρ < 1, there is convergence in distribution where ⇒ denotes convergence in distribution, D is the usual function space of right-continuous real-valued functions on [0, ∞) with left limits and the J 1 topology as in [50], and D ρ ≡ {D ρ (t) : t ≥ 0} is a stationary point process.Indeed, this stationary point process D ρ is the topic of the present paper.
In general, the stationary departure process may be well defined even if the limit in (1.1) does not exist, e.g., because the interarrival-time distribution is deterministic, creating periodic structure.
In this paper we consider the usual HT scaling of space and time applied to the stationary process; i.e., we consider the scaled stationary departure process defined by and we let ρ ↑ 1. Observe that convergence D * ρ ⇒ D * in D, which we establish in §5 is equivalent to the iterated limit i.e., we first let s ↑ ∞ and then afterwards we let ρ ↑ 1 with the scaling.
1.3.The Variance Function and the BRAVO Effect.Because of our interest in the IDC, we are especially interested in the associated variance function of the stationary departure process, defined in the obvious way by V ρ (t) ≡ V ar(D ρ (t)), t ≥ 0. As a consequence of the HT stochastic-process limit D * ρ ⇒ D * and appropriate uniform integrability, we will obtain the associated limit for the appropriately scaled variance functions, i.e., (1.4 where V * d (t) ≡ V ar(D * (t)), t ≥ 0. There is a recent related result for the M/GI/1 queue in Proposition 6 of [30], as part of an investigation into the BRAVO (balancing reduces asymptotic variance of outputs) effect; see [38], [26] and the earlier Theorem 4.1 of [9] and [54].For the M/GI/1 model, a two-term asymptotic expansion is developed for the variance function V ρ (t) ≡ V ar(D ρ (t)) as t → ∞ for fixed ρ < 1.There is no direct heavy-traffic scaling, but the scaling emerges in the parameters.We discuss the connections between our result and this earlier result in Remark 5.3.

Organization.
Here is how this paper is organized: We start in §2 by providing a brief review of stationary point processes, focusing especially on the variance function.In §3 we use Laplace transforms (LT's) of the stationary departure process in the GI/M/1 queue derived by [18,19] to derive the HT limit of its variance function.In §4 we use the HT limit for the Palm version of the mean function derived by [45] to derive the HT limit of the stationary variance function.(In §2.2 we review the application of the Palm-Khintchine equation to express the stationary variance in terms of the Palm mean function.)In §5 we establish the HT limit for the stationary departure process in the GI/GI/1 queue (Theorem 5.2) and its variance function (Theorem 5.3).In §6 we provide a brief overview of the application of Theorem 5.3 to support our RQNA-IDC developed in [53], which is briefly outlined in §6 of [51].We discuss extensions in §7.Finally, we present postponed proofs in §8.

Review of Stationary Point Processes.
In this section we review basic properties of stationary point processes; see [20] and [43] for more background.In §2.1 we review renewal processes and their Laplace transforms.In §2.2 we review the Palm-Khintchine equation and use it to express the variance function of a stationary point process in terms of the mean function of the Palm version.
2.1.Renewal Processes and the Laplace Transform.We start with a rateλ renewal process N ≡ {N (t) : t ≥ 0}.Let F be the cumulative distribution function (cdf) of the interval U between points (the interarrival time in a GI arrival process), having mean E[U ] = λ −1 and finite second moment.As a regularity condition for our queueing application, we also assume that F has a probability density function (pdf) f , where F (t) = t 0 f (u) du, t ≥ 0. Throughout this paper, we assume that the interarrival-time distribution of our renewal arrival processes has a pdf.That pdf assumption ensures that the equilibrium renewal process arises as the time limit of the ordinary renewal process; e.g., see §3.4 and §3.5 of [42].
The stationary or equilibrium renewal process differs from the ordinary renewal process only by the distribution of the first interarrival times.Let F e be the cdf of the equilibrium distribution, which has pdf f e (t) = λ(1 − F (t)).Let E e [•] denote the expectation under the stationary distribution (with first interval distributed according to F e ) and let E 0 [•] denote the expectation under the Palm distribution (with first interval distributed as F ).
Conditioning on the first arrival, distributed as F under the Palm distribution or as F e under stationary distribution, the renewal equations for the mean and second moment of N (t), the number of points in an interval [0, t], are: Throughout the paper, we use the Laplace Transform (LT) instead of the Laplace-Stieltjes Transform (LST).The LT of f (t) and the LST of F , denoted by L(f )(s) ≡ f (s), are as must be true for any stationary point process.
Let V (t) ≡ Var e (N (t)) be the variance process of N (t) under timestationary distribution.(We omit the e superscript on V (t) because we will only discuss stationary variance functions.)Combining (2.5) and (2.6), we have The variance function then can be obtained from the numerical inversion of the Laplace transform, e.g., see §13 of [3], [4] and [52].Term by term inversion shows that we can express V (t) in terms of the renewal function m(t) In §2.2 we show that the Palm-Khintchine equation can be used to derive a generalization of (2.8) for general stationary and ergodic point processes.
2.2.The Palm-Khintchine Equation.We now consider a continuous-time stationary point process, i.e., having stationary increments.The main idea is the Palm transformation relating continuous-time stationary processes to the associated discrete-time stationary processes.An important manifestation of that relation is the Palm-Khintchine equation; see Theorem 3.4.II. of [20].It is important here because it can be applied to generalize the variance formula discussed in §2.1; see §2.4 of [17] and §3.4 of [20].
We focus on orderly stationary ergodic point processes with finite intensity.(Orderly means that the points occur one at a time.)Let N (s, t] denote the number of events in interval (s, t], and N (t) ≡ N (0, t]. where q k (t) is the probability of exactly k arrivals in (0, t] under the Palm distribution, i.e., (2.10) Under ergodicity, the Palm distribution is equivalent to the event stationary distribution, so that q k (t) = P 0 (N (t) = k).
We now apply Theorem 2.1 to generalize (2.8) and (2.7) to the case of orderly stationary ergodic point process.
Corollary 2.1.(Variance of a stationary ergodic point process) For a general stationary ergodic point process with rate λ and finite second moment, the variance function is where and its LT is where m(s) is the LT of m(t).
Proof.Let (2.14) p k (t) = P e (N (t) = k), for k = 0, 1, 2, . . .so that k i=1 p i (t) = P e (N (t) ≤ k).With Theorem 2.1, we can write where m(t) ≡ E 0 [N (t)] as in (2.12).Taking the Laplace Transform, we obtain 3. The Departure Variance in the GI/M/1 Queue.We now start considering the queueing models.In particular, we focus on the GI/GI/1 queue, which has unlimited waiting space, the first-come first-served service discipline and independent sequences of i.i.d.interarrival times and service times distributed as random variables U and V , respectively, where U has a pdf f (t).Let λ ≡ 1/E[U ] be the arrival rate; let f (s) ≡ E e −sU be the LT of the interarrival-time pdf f (t); let µ ≡ 1/E[V ] be the service rate; and let ρ ≡ λ/µ be the traffic intensity, assuming that ρ < 1.
Daley [18,19] derived the LST of the variance V d (t) of the stationary departure process in a GI/M/1 queue.The associated LT of where ξ(s) is the root with the smallest absolute value in z of the equation and δ = ξ(0) is the unique root in (0, 1) of the equation which appears in the distribution of the stationary queue length in a GI/M/1 queue.Useful properties of ξ(s) and δ = ξ(0) are contained in Lemma 8.1.We now establish a HT limit for the departure variance function in the GI/M/1 model.To do so, we consider a family of GI/M/1 models parameterized by ρ, where λ ≡ 1/E[U ] and where γ ρ are positive constants such that lim ρ↑1 γ ρ = γ > 0. Note that if γ ρ = 1/ρ, then we come to the usual case of λ/µ = ρ.We allow this general scaling so that we can gain insight into reflected Brownian motion (RBM) with non-unit drift.Let the HT-scaled variance function be just as in (1.4).Throughout the paper, we use the asterisk (*) superscript with ρ subscript to denote HT-scaled items in the queueing model, as in (3.4), and the asterisk without the ρ subscript to denote the associated HT limit.
As should be expected from established HT limits, e.g., as in §5.7 and Chapter 9 in [50], the HT limit of the variance function V * d,ρ (t) in (3.4) depends on properties of the normal distribution and RBM.Let φ(x) be the pdf and Φ(x) the cdf of the standard normal variable N (0, 1).Let Φ c (x) ≡ 1 − Φ(x) be the complementary cdf (ccdf).Let R(t) be canonical RBM (having drift −1, diffusion coefficient 1) and let R e (t) be the stationary version, which has the exponential marginal distribution for each t with mean 1/2.The correlation function c * (t) of R e is defined as where H * 2 (t) is the second-moment cdf of canonical RBM in [1], which has mean 1 and variance 2.5; see Corollaries 1.1.1 and 1.3.4 of [1] and Corollary 1 of [2].The correlation function c * (t) has LT [1].Equivalently, the Gaussian terms in (3.5) can be re-expressed as φ( , where erf is the error function.This LT can be explicitly inverted, yielding By Corollary 1.3.5 of [1], the correlation function has tail asymptotics according to From the correlation function, we define It then follows from the explicit expression of c * (t) that One may verify that w * (t) is a increasing function satisfying 0 ≤ w * (t) ≤ 1.
As we shall see in Theorem 3.1, this w * (t) serves as the weight function that appears in the limiting departure variance function.
We now present the main result for the departure variance in the GI/M/1 special case.The idea of the proof is to exploit the explicit form of the LT V * d,ρ (t) of the scaled stationary departure variance and derive its HT limit.
We then obtain the convergence of the HT-scaled variance function V * d,ρ (t) by applying continuity theorem of the LT, see Theorem 2(a) in Chapter XIII of [21].The proof of the theorem can be found in §8.1.
Theorem 3.1.(HT limit for the GI/M/1 departure variance) Consider the GI/M/1 model with where γ ρ are positive constants such that lim ρ↑1 γ ρ = γ > 0. Assume that E[U 3 ] < ∞ so that a two-term Taylor series expansion of the LT f (s) about the origin is valid with asymptotically negligible remainder.Then the HTscaled variance function with ξ * (s) being the unique root with non-negative real part of the quadratic equation In addition, We shall want to relate our HT limit for the departure variance function to associated HT limits for the variance functions of the arrival and service processes.For that step, it is significant that the functional central limit theorem (FCLT) for any stationary point process has the same form as the FCLT for the associated Palm process, as was shown by [39].Extra uniform integrability is required to get the associated limit for the variance function.To be relatively self-contained, we will directly derive the desired result from the transform of the equilibrium renewal process in (2.7).
For that purpose, let N a (t) denote the arrival renewal process and let V a (t) ≡ Var e (N a (t)) denote its variance process under the stationary distribution.Similarly, we define N s (t) and V s (t) for the service renewal process.The following lemma states that the terms c 2 a λt and c 2 s λt in (3.14) can be interpreted as the limiting variance function of the arrival and service renewal processes, respectively.This implies that the limiting departure variance function V * d is a convex combination of the arrival and service variance functions with a scaled version of the time-varying weight function w * (t).This convex combination result is consistent with the more elementary approximation used in QNA; see [47,48]; there the departure variability parameter is approximated by a convex combination of the arrival and service variability parameters.Lemma 3.1.(Limiting variance function of stationary renewal processes) Let N (t) be a renewal process with rate λ and let c 2 N be the scv of the interrenewal distribution.Consider the HT-scaled stationary variance function Proof.Let f denote the inter-renewal distribution.Recall the expression for the LT of a stationary renewal process in (2.7), we have The result follows from inverting the LT, i.e., V * s (t) = λc 2 N t.To derive a pre-limit approximation, define the weight function where V a (t) and V s (t) are the variance functions associated with the equilibrium arrival and service renewal processes, both with rate λ.Define the HT-scaled weight function Combining Theorem 3.This justifies the following approximation for the variance function of the stationary departure process from a GI/M/1 queue: We conclude this section with the tail asymptotic behavior of the variance function.To start, we re-write V * d in terms of c * and H * 2 , where c 2 s = 1.
Example 3.1.(numerical experiments) We now evaluate the approximation stemming from Theorem 3.1.First, we use numerical transform inversion as in [3,4] to calculate the limiting variance function.
Figure 1 (left) reports V * (t) − λc 2 a t for four sets of parameters such that the limiting constant (1 − c 4 a )/2γ 2 in Corollary 3.2 will be 1.5, 0.375, −1.5 and −1.5, respectively.Figure 1 (right) confirms Theorem 3.1 by comparing simulation estimates of the HT-scaled departure variance function a t for ρ = 0.8 and 0.9 from simulation with the theoretical limit V * (t)−λc 2 a t for the E 2 /M/1 model, where c 2 a = 0.5, showing that the theoretical limit in (3.12) serves as a good approximation of the HT-scaled variance function.On the left is V * (t) − λc 2 a t for four sets of parameters calculated from numerically inverting (3.12).On the right is the the HT-scaled variance a t for ρ = 0.8 and 0.9 in the E2/M/1 model, estimated by simulation, compared with the theoretical limit V * (t) − λc 2 a t.
4. The Departure Variance in the M/GI/1 Queue.We now prove that the HT limit for the stationary departure variance in (3.14) also holds true for the M/GI/1 model.Of course, here we restrict our attention to c 2 a = 1 instead of c 2 s = 1 before.Theorem 5.3 will show that the same formula is valid for GI/GI/1 with general c 2 a and c 2 s .Recall from (2.13) that the Laplace Transform of the variance function of a general stationary and ergordic point process is In the case of the M/GI/1 model, [45] (on p. 78) derived an expression for md (s).
, where ĝ(s) = E e −sV is the Laplace Transform of the service time pdf g(t), ν(s) is the root with the smallest absolute value in z of the equation and is the probability generating function of the distribution of the stationary queue length Q.
, is exactly the Laplace Transform of the mean process of the service renewal process.Now, we state the HT limit in terms of the HT-scaled variance function defined in (3.4) for the M/GI/1 special case.The result parallels that for the GI/M/1 case.The proof can be found in §8.2.Theorem 4.2.(HT limit for the M/GI/1 departure variance) Consider an M/GI/1 model with where γ ρ are positive constants such that lim ρ↑1 γ ρ = γ > 0. Assume that E[V 3 ] < ∞ so that a two-term Taylor series expansion of the LT ĝ(s) about the origin is valid with asymptotically negligible remainder.Then the HTscaled variance function V * d,ρ (t) defined in (3.4) converges as ρ ↑ 1, i.e., (4.4) where the limit V * d (t) is a continuous nonnegative function with LT with ν * (s) being the unique root with positive real part of the equation In addition, With the same technique as in Corollary 3.2, one can prove the following corollary, which yields exactly the same asymptotic behavior.Proposition 6 of [30] developed a two-term asymptotic expansion for the variance function V ρ (t) ≡ V ar(D ρ (t)) as t → ∞ for for the M/GI/1 queue and fixed ρ < 1.We discuss the connections between our result and this earlier result in Remark 5.3.

Heavy-Traffic Limit for the Stationary Departure Process.
In this section, we establish an HT limit for the stationary departure process and its variance function in a GI/GI/1 queue.To do so, we apply the recent HT results for the stationary queue length (number in system) in [25] and [12] together with the HT limits for the general single-server queue in §9.3 of [50] and the general reflection mapping with non-zero initial conditions in §13.5 of [50].As in [50], a major component of the proof is the continuous mapping theorem.
The corresponding limit starting out empty is contained in Theorem 2 of [33].There has since been a substantial literature on that case; see [26,34,50].As can be seen from §9.3 and §13.5 of [50], for the queue length, the key map is the reflection map ψ applied to a potential net-input function x, with a ∧ b ≡ min {a, b}, so that ζ(x) ≤ 0 and ψ(x)(t) ≥ x(t) for all t ≥ 0.
The key point is that we now allow x(0) = 0.

5.1.
A General Heavy-Traffic Limit for the G/G/1 Model.For the general G/G/1 single-server queue with unlimited waiting space and service provided in the order of arrival, we consider a family of processes indexed by the traffic intensity ρ, where ρ ↑ 1.Let Q ρ (t) be the number of customers in the system at time t; let A ρ (t) count the number of arrivals in the interval [0, t], which we assume to have rate λ; let S ρ (t) be a corresponding counting process for the successive service times, applied after time 0, to be applied to the initial Q ρ (0) customers and to all new arrivals; let B ρ (t) be the cumulative time that the server is busy in the interval [0, t].Then the queue-length process can be expressed as where the three components are typically dependent.(For simplicity, we assume that A ρ (0) = S ρ (0) = B ρ (0) = 0 w.p.1.) We have in mind that the system is starting in steady-state.Thus the triple (Q ρ (0), A ρ (•), S ρ (•)) is in general quite complicated for each ρ.Even in the relatively tractable GI/GI/1 cases, which we shall primarily treat, the residual interarrival time and service time at time 0 will be complicated, depending on ρ and Q ρ (0).We will need to make assumptions ensuring that these are uniformly asymptotically negligible in the HT limit.
By flow conservation, the departure (counting) process can be represented as Directly, or by combining (5.3) and (5.4), (5.5) be a net-input process, acting as if the server is busy all the time, and thus allowing X ρ (t) to assume negative values.The cumulative busy time B ρ (t) is then relate to X ρ (t) by As a consequence of the assumptions above, X ρ (0) = Q ρ (0).Roughly, for ψ in (5.1), but the exact relation breaks down because the service process shuts down when the system becomes idle, so that a new service time does not start until after the next arrival.While (5.8) does not hold exactly for each ρ, it holds in the HT limit, as shown in Theorem 9.3.4 of [50].It would hold exactly if we used the modified system in which we let the continuoustime service process run continuously, so that equation (5.8) holds as an equality, as done by [11] and then again in §2 of [32].Because the modified system has been shown to be asymptotically equivalent to the original system for these HT limits in [11] and [32], that is an alternative approach.
We now introduce HT-scaled versions of these processes, for that purpose, let Let D be the function space of all right-continuous real-valued functions on [0, ∞) with left limits, with the usual J 1 topology, which reduces to uniform convergence over all bounded intervals for continuous limit functions.Let D k be the k-fold product space, using the product topology on all product spaces.Let ⇒ denote convergence in distribution.Let e be the identity function in D, i.e., e(t) ≡ t, t ≥ 0.
Proof.First, note that (5.13) ρ and S * ρ have different translation terms in (5.9), ensuring that the potential rate out is λ/ρ, which exceeds the rate in of λ, consistent with a stable model for each ρ, 0 < ρ < 1.Hence, under the assumption, ) is obtained by exploiting the relationship in (5.7).The limits for Q * ρ and D * ρ then follow from the continuous mapping theorem after carefully accounting for the busy and idle time of the server; see the proof of Theorem 9.3.4 and preceding material in [50].

5.2.
A Heavy-Traffic Limit for the Stationary Departure Process.Theorem 5.1 is not easy to apply to establish HT limits for stationary processes because condition (5.10) is not easy to check and the limit in (5.11) and (5.12) is not easy to evaluate.
In order to establish a tractable HT limit for the stationary departure process, we apply the recent HT limits for the stationary queue length in [25] and [12].Their HT limits are for generalized open Jackson networks of queues, which for the single queue we consider reduce to the GI/GI/1 model.Following [12], we assume that the interarrival times and service times come from independent sequences of i.i.d.random variables with uniformly bounded third moments (2 + ǫ would do).
Theorem 5.2.For the GI/GI/1 model indexed by ρ, assume that (i) the interarrival-time cdf has a pdf as in §2.1 and (ii) the interarrival times and service times have means λ and λ/ρ, scv's c 2 a and c 2 s , without both being 0, and uniformly bounded third moments.Then: (a) For each ρ, 0 < ρ < 1, the process Q * ρ can be regarded as a stationary process, while the process D * ρ can be regarded as a stationary point process (with stationary increments).
(b) Condition (5.10) in Theorem 5.1 holds with where B a and B s are independent standard (mean 0, variance 1) Brownian motions (BM's) that are independent of Q * (0), which is distributed as R e (0) with R e being a stationary RBM with drift −λ and variance λc 2 x ≡ λc 2 a +λc 2 s , and so an exponential marginal distribution, i.e., (5.15) (c) The limits in Theorem 5.1 hold, where (5.16) , B a and B s being mutually independent.(d) We have  Proof.First, recall that the HT limit as ρ → 1 starting empty is the RBM which converges as t → ∞ to the exponential distribution in (5.15).We will be applying [25] and [12] to show that the two iterated limits involving ρ → 1 and t → ∞ are equal.Toward that end, we observe that, by § §X.3-X.4 of [7], the queue-length process has a proper steady-state distribution for each ρ.As on p. 63 of [25], we add the residual interarrival times and service times to the state description for Q ρ (t) to make it a Markov process that has a unique steady-state distribution for each ρ.These residual interarrival and service times are asymptotically negligible in the HT limit.The associated departure process D ρ (t) then necessarily is a stationary point process for each ρ.We then can apply Theorem 8 of [25] to have a limit for the scaled stationary distributions, so that condition (5.10) holds with (5.14).Since strong moment-generating-function-condition is imposed in (1) and (2) on p. 62 of [25], we apply the extension in Theorems 3.1 and 3.2 of [12] to cover our moment condition.Hence, we can apply Theorem 5.1 with these special initial distributions to get the associated process limits in the space D.
We now establish an HT limit for the variance of the stationary departure process.The form of that limit is already given in Theorem 3.1.The proof can be found in §8.3.
Theorem 5.3.(Limiting variance) Under the conditions of Theorem 5.2 plus the usual uniform integrability conditions, for which it suffices for the interarrival times and service times to have uniformly bounded fourth moments, As a by-product, the covariance formulas in (5.21) can be generalized to describe the covariance of between a stationary RBM and a BM, where the underlying BMs are correlated.
Remark 5.1.(the quasireversible case) The limit process where as in Theorem 5.2, can be called the Brownian queue; see [27,28,29,40].The Brownian queue is known to be quasireversible if and only if c 2 a = c 2 s .In that case, the stationary departure process is a BM and the departures in the past are independent of the steady-state content.Consistent with that theory, V * d (t) = c 2 a λt, t ≥ 0 in (5.20) if and only if c 2 a = c 2 s .
Remark 5.2.We considered only the case where µ ρ = λ/ρ in this section.Now, we list the results for a slightly more general case as in §3, where we have µ ρ = λ(1 + (1 − ρ)γ ρ ) and lim ρ↑1 γ ρ = γ.One can easily check that Theorem 5.1 holds with Theorem 5.2 holds with and Theorem 5.3 holds with We conclude this section with the tail asymptotics of the departure variance function.Just as in Corollaries 3.
Remark 5.3.Hautphenne et al. [30] developed explicit expressions for the y-intercept bθ of the linear asymptote for the variance of the stationary departure from M/GI/1 see Proposition 6 there.Their result (i) holds for M/GI/1 case; (ii) depends on the third moment of the service distribution; (iii) holds for general traffic intensity.Even though there is no direct heavy-traffic scaling in their result, the scaling parameter emerges in their expression, see the definition of bθ there.In specific, the scaling constant ρ/(1 − ρ) 2 coincides (up to a constant ρ) with the one we use in (1.4).
On the other hand, our result here (i) coincides with their y-intercept (after scaling) in the HT limit in the M/GI/1 case, i.e., let ρ = 1, γ = 1 and c a = 1; (ii) holds for general GI/GI/1 cases but only under HT limit; (iii) has explicit characterization of the remainder term, again only under HT limit.
6. Application to a Robust Queueing Network Analyzer.We conclude by explaining the important role that Theorem 5.3 plays in our Robust Queueing Network Analyzer (RQNA) based on the index of dispersion for counts (IDC), which we refer to as RQNA-IDC.In §6.1 we briefly review the robust queueing (RQ) approximation for the mean steady-state workload of a G/G/1 queue developed in [51], which requires the IDC of the arrival process as model data in the G/GI/1 model.Then, in §6.2 we review the approximation for the IDC of the departure process that we propose in [53], which is supported by this paper.
6.1.The Mean Steady-State Workload at a G/G/1 Queue.In this section we review the RQ approximation for the mean steady-state workload at a G/G/1 queue developed in [51].We start with a rate-λ stationary arrival process A ≡ {A(t) : t ≥ 0}, having stationary increments.Thus, for a renewal arrival process, we work with the associated equilibrium renewal process.We also start with a stationary and ergodic sequence of rate-µ service times {V k : k ≥ 1} with finite scv c 2 s and possibly correlated with the arrival process.For model stability, we assume that ρ = λ/µ < 1.Let Y (t) be the associated total input process, defined by (6.1) which also has stationary increments.Let Z ρ (t) be the workload at time t in model ρ with arrival rate ρ and service rate 1, then As in §3.1 of [51], we can apply a reverse-time construction to write the steady-state workload Z ρ as a simple supremum.
For the RQ approximation of the mean steady-state workload, we use the associated index of dispersion for work (IDW) I w ≡ {I w (x) : x ≥ 0} introduced by [22], which is defined by (6.2) , t ≥ 0; see §4.3 of [51] for key properties.Let Z ρ denote the steady-state workload in the G/G/1 model when the traffic intensity is ρ, then the RQ approximation (based on this partial model characterization) is for the IDW in (6.5).The approximation (6.3) comes from ( 28) of [51], assuming that we set the parameter b f = √ 2, which makes the approximation asymptotically correct for the GI/GI/1 model in both the heavy-traffic and light-traffic limits; see Theorem 5 of [51].Notice that the approximation in (6.3) is directly a supremum of a real-valued function, and so can be computed quite easily for any given tuple (ρ, I w ).
In what follows, we primarily focus on the G/GI/1 special case, where the arrival process is general stationary and ergodic, and the service times are i.i.d.We assume that the arrival process A is partially characterized by its index of dispersion for counts (IDC) I a (t), where (6.4) [16].(To explain, just as the distribution of a random variable is not fully characterized by its mean and variance, so the distribution of the stationary stochastic process {A(t) : t ≥ 0} is not fully characterized by its mean function E[A(t)] and its variance function V ar(A(t)), t ≥ 0.) In this context, (6.5) as noted in §4.3.1 of [51].The corresponding RQ approximation is then Hence, to apply the RQ algorithm in the G/GI/1 model, we need only develop an approximation for the arrival IDC I a (t).
6.2.The IDC of a Stationary Departure Process.The main challenge in developing a full RQNA-IDC involving a decomposition approximation is calculating or approximating the required IDC for the arrival process at each queue.For a renewal arrival process, the IDC I a can be computed by inverting the LT V (s) in (2.7) using the LT m(s) of the renewal function m(t) in (2.2), which only requires the LT f (s) of the interarrival-time pdf in (2.1).Numerical algorithms for calculating and simulations algorithms for estimating the IDC are discussed in [52].
As discussed in §6 of [51], the main challenge is approximating the IDC of a departure process from a G/GI/1 queue, partially specified by the tuple (λ, ρ, I a , c 2 s ).As in §6 of [51], we propose approximating the IDC of the departure process, I d,ρ (t) by the weighted average of the IDC's of the arrival and service processes, i.e., (6.7) where I a (t) and I s (t) are the IDCs associated with the equilibrium arrival and service renewal processes, both with rate λ.We thus require (λ, ρ, I a , I s ) as model data.Based on initial study, a candidate weight function w ρ (t) was suggested in (43) of [51].
We now show how the further study in Theorem 5.3 suggests a new weight function (6.8) w ρ (t) ≡ w * ((1 − ρ) 2 λt/(ρc 2 x )), where c 2 x ≡ c 2 a + c 2 s and w * is given in (3.9) with c * in (3.5).To start, let I d,ρ denote the departure IDC and define the weight function (6.9) where V a (t) and V s (t) are the variance functions associated with the equilibrium arrival and service renewal processes, both with rate λ.Note that this is exactly the same weight function we defined in (3.9), thus we have the same HT-scaled weight function w * ρ as in (3.16).We then apply Theorem 5.3 to obtain a paralleling corollary.Corollary 6.1.(Limiting weight function) Under the assumptions in Theorem 5.3, we have w * ρ (t) ⇒ w * (λt/c 2 x ) for w * defined in (3.9).To obtain the pre-limit approximation, we re-arrange terms in (6.9), and assume (6.10) . One remaining issue is that the approximation (6.10) does not automatically yield a correct light traffic limit, in which case we must have I d,0 (t) = I a (t) since the service time is negligible.As a remedy, we propose to add a constant ρ −1 correction in the weight function, so that we have (6.8) as the final weight function.
We conclude this section with two simulation examples.
Example 6.1.(simulation experiments) We now consider two GI/GI/1 queues, where neither the interarrival time nor service time has an exponential distribution.Let H 2 (c 2 ) be the H 2 (hyperexponential) distribution with scv c 2 and balanced means, as in (3.7) on p. 137 of [46].Consider the H 2 (4)/E 2 /1 model with λ = 2, Figure 2 (left) reports the simulated departure IDCs for three different traffic intensities ρ = 0.95, 0.8, 0.5, as well as the approximation (6.7) with (6.8).The simulation estimation of the departure IDC is obtained from a single run of length 10 9 time units, with the first 10 6 time units are discarded in order for the system to approach steady-state.The reference IDCs I a and I s is calculated by numerically inverting the LT in (2.7). Figure 2 (right) is the corresponding plot for the E 2 /H 2 (4)/1 model with λ = 2.

7.
Extensions.The approximation for the departure IDC I d (t) in (6.7) and (6.8) should be good for much more general models than GI/GI/1, with the independence conditions relaxed and more than 1 server.We also conjecture that the HT limit of the variance function in Theorem 5.3 extends to a larger class of models as well.Indeed, we conjecture that the limits established for GI/GI/1 extend in that way.First, Theorem 5.1 extends quite directly by exploiting [32,33].For the extension of Theorem 5.2, there is a large class of models for which the HT-scaled arrival and service processes have the limits where B a and B s are independent standard (mean 0, variance 1) Brownian motions (BM's) that are independent of the initial queue length.What is needed is the extension of [25] and [12] to more general models.We conjecture that can be done for GI/GI/s and other models with regenerative structure in the arrival and service processes.For GI/GI/s the queue-length process again becomes a Markov process if we append the s elapsed service times as well as the elapsed interarrival time, but it remains to do the hard technical analysis leading to an appropriate Lyapunov function.
It is also of interest to establish related results for departure processes in models with non-renewal arrival processes, as in [23] and references therein.It also remains to establish new HT limits for stationary departure processes from a queue within a network, obeying the HT FCLT in [41].
The relevant approximation for the stationary departure process from a many-server GI/GI/s queue evidently is quite different, being more like the service process than the arrival process.We conjecture that the relevant many-server heavy-traffic limit for the stationary departure process is a Gaussian process with the covariance function of the stationary renewal processes associated with the service times, as in the CLT for renewal processes in Theorems 7.2.1 and 7.2.4 of [50].Partial support comes from [5], Appendix F of [6] and [24].

Proofs.
8.1.Proof of Theorem 3.1.We now review a useful lemma on the properties of ξ(s) and δ; see p. 113 of [45] or Appendix 6 of [15].(The notation here is slightly different.)that has the smallest absolute value is This root ξ(s) is a continuous function of s for Re(s) ≥ 0. Furthermore, z = ξ(s) is the only root in the unit circle |z| ≤ 1 if at least one of two conditions is satisfied (i) Re(s) > 0, or (ii) Re(s) ≥ 0 and λ/µ < 1.Specifically, δ = ξ(0) is the smallest positive real root of the equation If λ/µ < 1, then δ < 1 and if λ/µ ≥ 1 then δ = 1.
As a consequence, we pick the root of (8.6) with non-negative real part.In particular, for real s, we have For complex s, the square root in (8.7) corresponds to two complex roots, which are also the roots of γ − γ 2 + 2(c where ξ * is defined in (8.5).Moreover, we have Combining everything into the Laplace Transform of ( Plugging in (8.7), we obtain where we pick the root such that 1 + 2(c 2 a + 1)s/(λγ 2 ) − 1 /(s/(λγ 2 )) has non-negative real part.We used the fact that Re(z) ≥ 0 if and only if Re(1/z) ≥ 0 for z = 0.
For the explicit inversion, one can exploit the LT of the correclation function in (3.6) and note that for any constant a = 0 and any function f with LT f .For our case here, we use a = λγ 2 /(c 2 a + 1). where One can easily show that F (1) ρ (s) converges to 0 as ρ ↑ 1.Note also that ĝ(0) = 1 and ĝ′ (0 Furthermore, a Taylor series expansion around s = 0 yields where we used the fact that ĝ′′ (0 With essentially the same argument as in the proof of Theorem 3.1, one can also show that ν * (s) is the only root of (8.12) with positive real part, furthermore It remains to show that F (2) ρ (s) converges (pointwise) to a proper limit.
To this end, we write where we apply (8.12) to obtain the simplified expression in (8.15).
To obtain the explicit inversion, we write Then, one exploit the LT of the correclation function in (3.6) and note that L(f (at))(s) = f (t/a)/a, for any constant a = 0 and any function f with LT f .For our case here, we use a = µγ 2 /(1 + c 2 s ).Proof.By combining Theorems 2.1 and 4.2 in Chapter X of [7], we deduce that the k th moment of the steady-state queue length is finite if the (k + 1) st moments of the interarrival time and service time are finite.We add the extra uniformly bounded fourth moment to provide the uniform integrability needed to get convergence of the moments in the HT limit.We use (5.4) to obtain the corresponding result for the departure process.
To get (5.20),combine (5.19) and (5.17 see §2 of [2] or Theorem 5.7.11 of [50].Inserting (8.17) into (8.16)yields the first line in (5.20) above.To establish the second limit, we do a space-time transformation of the limit, so that the limit is the same as one of the models analyzed directly.Let us re-scale space and time so that the general result is in terms on B a instead of c a B a (assuming that c a > 0), so that we can apply the established result for the M/GI/1 model.(Essentially the same argument works for GI/M/1.)The first step is to observe that the HT limit for the departure process {D * (t) : t ≥ 0} can be written as a function Ψ : R × D 3 → D of the vector process {(Q * (0), c a B a • λe, c s B s • λe, −λe)}; i.e., by (5.17If we replace the basic vector process (Q * (0), c a B a • λe, c s B s • λe, −λe) by another that has the same distribution as a process, then the distribution of D * will be unchanged.
By the basic time and space scaling of BM, for c a > 0, the stochastic processes have equivalent distributions as follows {Q * (0), c a B a (λt), c s B s (λt), −λt} where u = λt/c 2 a .After this transformation, to describe the system at time u, the associated RBM has drift −1 and variance coefficient 1+(c 2 s /c 2 a ) = c 2 x /c 2 a .Note that the mean of the steady-state distribution associated with the new RBM is the diffusion coefficient divided by twice the absolute value of the drift, which is c 2 x /(2c 2 a ).As a result, Q * (0)/c 2 a is exactly the steady-state distribution needed for the new RBM.From above, we see that x ) = (λ, 1, c 2 x /c 2 a ), so that w * (t) = w * (c 2 a λt/c 2 x ).
We now turn to the variance.By applying (4.7), we obtain
.18) Combining (3.8) and (3.18), we obtain the asymptotic behavior of the departure variance function.Corollary 3.2.(Asymptotic behavior of the departure variance function) Under the assumptions in Theorem 3.1, Fig 1.On the left is V * (t) − λc 2 a t for four sets of parameters calculated from numerically inverting (3.12).On the right is the the HT-scaled variance(1 − ρ) 2 V ((1 − ρ) −2 t) − λc 2a t for ρ = 0.8 and 0.9 in the E2/M/1 model, estimated by simulation, compared with the theoretical limit V * (t) − λc 2 a t.

Theorem 4 . 1 .
(Laplace transform of the Palm mean function) For the departure process from a M/GI/1 queue,

Corollary 4 . 1 .
(Asymptotic behavior of the departure variance curve) Under the assumptions in Theorem 4.2, we have the limit in(3.19),except now c 2 a = 1 and c 2 s is general.
2 and 4.1, we have Corollary 5.2.(Asymptotic behavior of the departure variance function) Under the assumptions in Theorem 5.3 and Remark 5.2, (5.22)

5 Fig 2 .
Fig 2.On the left is the IDCs of the departure from H2(4)/E2/1 model with λ = 2 and three different traffic intensities, together the two reference IDCs, one for H2(4) and one for E2, displayed in broken black lines.The approcximation (6.7) and (6.8) are displayed in broken colored lines.On the right is a similar plot for the E2/H2(4)/1 model.

(7. 1 )
A * ≡ c a B a and S * ≡ c s B s ,
is the correlation function in (3.5) and w * (t) is the weight function in (3.9); i.e., V * d (t) is given in (3.14) with γ = 1, but allowing general c 2 s .Moreover, we have the covariance formulas