MEAN TIME FOR THE DEVELOPMENT OF LARGE WORKLOADS AND LARGE QUEUE LENGTHS IN THE GI/G/1 QUEUE by

We consider the GI/G/1 queue described by either the workload U(t) (unfinished work) or the number of customers N(t) in the system. We compute the mean time until U(t) reaches excess of the level K, and also the mean time until N(t) reaches N0. For the M/G/1 and GI/M/1 models, we obtain exact contour integral representations for these mean first passage times. We then compute the mean times asymptotically, as K and N0→∞, by evaluating these contour integrals. For the general GI/G/1 model, we obtain asymptotic results by a singular perturbation analysis of the appropriate backward Kolmogorov equation(s). Numerical comparisons show that the asymptotic formulas are very accurate even for moderate values of K and N0.


Introduction
Queueing models arise in a wide variety of applications such as computer systems and com- munications networks.The mathematical analysis of these models typically involves the compu- tation of certain performance measures, such as the steady-state queue length or workload distri- bution, the length of a busy period, etc.Of particular importance is the total number of custom- ers ("jobs") or the size of the workload in the queueing system.If this total size becomes very large, the system performance may deteriorate and jobs may be lost or suffer long delays.For example, in the design of high speed communications systems, the buffer size at a switch is of cru- cial importance.If the buffer capacity is too small, arriving jobs may be frequently lost.Even if the buffer has a large capacity, a large buffer size below capacity may develop, resulting in un- acceptably long delays.It is therefore natural to ask how long it will take before a large buffer size (as measured by either number of jobs or total workload) develops in a particular queueing system.
In this paper, we consider the GI/G/1 queue and compute the mean time needed for either the workload (total unfinished work) or the number of customers in the system, to reach or 1This research was supported in part by NSF Grant DMS-93-00136 and DOE Grant DE- FG02-93ER25168.
108 CHARLES KNESSL and CHARLES TIER exceed some specified level.We denote the workload at time t by U(t) and by N(t) the number of customers.These stochastic processes are not Markovian, but may be imbedded in higher di- mensional Markovian processes by using the supplementary variable technique.If we let Z(t) be the elapsed time since the last arrival and Y(t) be the elapsed service time of the customer pre- sently being served, then (U(t),Z(t)) and (N(t), r(t),Z(t)) are both Markov processes.
It is generally undesirable to have large queue lengths or large workloads.In a stable GI/G/1 queue the occurrence of these events is very rare.However, having accurate measures of the probabilities of such rare events may be an important measure of performance and reliability.
We therefore analyze the time for U(t) to reach or exceed some large level K, and also the time for N(t) to reach some large number N 0. For a model with a finite capacity of either workload or queue length, computing this time is the same as computing the time until the next customer is lost.
For the GI/G/1 model, we denote the interarrival time density by a(. and the service time density by b(. ).Throughout the paper we assume that the Laplace transforms are analytic for Re(0) < 0. moments by Thus, all the moments are finite, and we denote the respective m k J ykb(y)dy, M k-/ zka(z)dz, k>_ 1. 0 0 The traffic intensity is p-m 1/M 1 and we assume that p < 1, which guarantees that the queue is stable, and the processes N(t) and U(t) have steady-state distributions.For exponential arrivals we set M 1 1/A and for exponential service we set m I 1/#.Let p(w) be the steady-state workload density, and let p(n) be the steady-state probability that N(t) n.For these quantities there are exact integral representations.A good summary of known exact results for the GI/G/1 model can be found in the book of Cohen [4].From these integrals and our assumptions about the analyticity of the Laplace transforms, it is easy to obtain the asymptotic tail probabilities p(w) ale cw, w--,oo   (1.1) e CZa(z)dz (1.3) Thus, the tail exponent in the workload is precisely -c, and the tail exponent in the queue length distribution is log[3(c)].The constant c corresponds to a simple pole in the integral repre- sentations of the distribution functions.The constants c1, c 2 correspond to residues at this pole, and these are also readily computed.In sections 6 and 7 we explicitly identify al and a2" Now let N be the mean time until the workload reaches or exceeds the level K.For the GI/G/1 model, N will generally depend upon the initial workload U(0)= w, and also on the ini- tial value Z(0)= z.Thus, g N(w,z).We let T be the mean time until N(t)reaches N o This function depends on the initial values (g(0), Y(0), Z(0)) (n, y, z), i.e., T T(n, y, z).We will show, however, that asymptotically as K, N0---oc (and for initial values of U(0), N(0) not close to K, No) the mean first passage times are independent of w, n, y, z; and we have N(w,z) t1 ecK, K--,oo   (1.4) T(n,y,z) 2 eCUb(y)dy No-Oc.
(1.5) o Thus, N and T grow exponentially at precisely the same rates as the corresponding steady-state probabilities decay.While this seems to be well-known and can be argued, for example, by using (1.1)-(1.2) and renewal theory, the computation of the constants /1, 2 appears to be a much more difficult task.The purpose of this paper is to give explicit formulas for these constants.
The computation of 1 and /?2 is essential if one is to obtain accurate numerical approximations to N and T; this is further discussed in section 8.
Previous work on tail exponents and tail probabilities includes Cohen [2], Iglehart [8], Neuts and Takahashi [16] and Sadowsky and Szpankowski [17].These authors consider the m-server GI/G/m queue and/or various special cases of this model.Using the tail behavior as an exponen- tial approximation is also discussed by Fredericks [5], Gaver and Shedler [6], and in the book of Tijms [18].This type of approximation seems to be superior to the standard (exponential) heavy traffic approximation, and reduces to the latter in the heavy traffic limit (for the GI/G/1 model this is defined as pT1).The difficulty in using this approximation is the computation of the constants 1, O2" The asymptotic approach employed here makes use of singular perturbation methods such as boundary layer theory and asymptotic matching; general references for these techniques are Kevor- Man and Cole [9] and Bender and Orszag [1].Matkowsky and Schuss [15] developed a singular perturbation method for computing asymptotically first passage times for diffusion processes with small diffusion coefficients.We have extended this method to discrete random walks and Markov jump processes in [11][12][13].In particular, in [12, 13] we computed the first passage times N and T for a state-dependent M/G/1 model, where the arrival and service processes are allowed to depend on the present workload or queue length.Results for the standard M/G/1 model may be ob- tained simply by omitting the various state-dependence.Recently [10], we have computed the mean time for large queue lengths to develop in tandem Jackson networks.Numerical compari- sons in [10] show that the asymptotic results are in excellent agreement with exact (numerical) solutions.For this agreement to occur it is essential to compute the constants in the asymptotic expansions (i.e., the analogs of 2 in (1.5)).
For the M/M/1 model, it is trivial to compute exactly the first passage times N and T. For the more general M/G/1 and GI/M/1 models, we give exact integral representations for N and T (cf.Results 2-5).It is easy to then evaluate these integrals asymptotically and hence obtain (1.4) and (1.5).For these special models, we also obtain the asymptotics directly, by using perturba- tion techniques to analyze the appropriate backward Kolmogorov equation(s).Of course, the two methods yield identical results.For the general GI/G/1 model, we have not been able to obtain exact expressions for N and T, so that we use the perturbation method to obtain asymptotic for- mulas for K, N0--,oc.The main results for the GI/G/1 model are summarized in Results 6 and 7, and these appear in sections 6 and 7, respectively.
While we only compute mean first passage times, similar techniques may be used for higher moments.However, in the asymptotic limit considered here, the first passage times tend to be exponentially distributed, so that the mean is sufficient to characterize the entire distribution.This was shown explicitly for singularly perturbed diffusion processes by Williams [19], and this calculation can be easily adapted to discrete random walks and jump processes.
The assumption on the analyticity of the Laplace transforms if(-) and b (-) is essential to our analysis.If, say, the service time density had only an algebraic tail, then it is likely that the mean first passage times N and T would have only algebraic growth in K and No, and also be more sensitive to the initial conditions w and n.
2. Queue Length in the M/M/1 Queue We compute the mean time until N(t), the number of customers in the system, reaches some large number N o Thus we define Note that since customers arrive individually, the time to reach or exceed N o is the same as the time to reach N o The function T(n) satisfies the backward equation AT(n + I + #T(n-1) (A + #)T(n) -1, l <_ n <_ N o-1 (2.5) Here .kand # are the arrival and service rates, respectively.This difference equation is easily solved and then the result may be asymptotically evaluated in the limit N0-oo.We set p-A/# and give the results below.
Result 1: The mean time for N(t) to reach N o in the M/M/1 queue is given by Even this very simple model reveals some important insights into the structure of T(n).
First, we note that T(n) grows exponentially as N0oe.Also, T(n) is independent of the initial value N(0) n, as long as n is not close to the "exit boundary" N 0. In [10], we have shown how to obtain the asymptotic results in (a) and (b) directly from the difference equation, by using sin- gular perturbation techniques.This method was then extended to compute the time needed for large queue lengths to build up in tandem Jackson networks.For these problems exact results are not available, so that the direct asymptotic approach is needed.
We note that the last term in the exact expression for T(n) (the term linear in n) is uniform- ly smaller than the first term(s) in the limit N0--<x 0 <_ n < N 0. It is in fact, exponentially smaller.This last term is a particular solution to the difference equation (2.3), which has 1 as an inhomogeneous term.The parts of the exact expression for T(n) that grow exponentially as n, N0cx satisfy the homogeneous form of (2.3).These observations are useful for developing the perturbation method.
3. Workload in the M/G/1 Queue We consider an M/G/1 model with arrival rate A and service time density b(. ), with finite moments of all order.The workload U(t) is clearly a Markov process.We define the first passage time ', min[t:U(t) > K] (3.1) and its conditional mean N(w) E(r, U(O) w), 0_<w<K. (3.2) By definition, N(w) 0 for w _> K.Note that U(t) can jump across the level K without actually hitting it.In this respect, the jump process is different from a diffusion process or the discrete queue length process.
The backward equation satisfied by N(w) is The last equation may be replaced by the equivalent condition N'(0 +)-0, which is obtained by setting w-0 + in (a.3) and subtracting equation (3.4).(3.6) Here N(h') is understood to mean N(K-), as N(K + O.The solution to (3.5) and (3.6) is 1 e (" ')K[1 pe (" -,X)(K )] + #(w K) 1 (3.7) N(w) A(1-P) 2 #-A As was the case for T(n), we see that N(w)is exponentially large in K and nearly constant, except when K-w O(1), which corresponds to initial workloads close to the "capacity" K.
We observe that N(K-> 0 N(K + ), so that the first passage time has a discontinuity at w-K.This is because for initial workloads just below K, the system's unfinished work cannot exceed K until the next customer arrivers, and by the time this occurs the server has decreased the workload, possibly by a significant amount.
Asymptotically, as K---.oo,we have N(K-).. (1-p)N(O) so that N(K-) is in fact exponentially large in K, of the same order of magnitude as N(0), the mean first passage time starting with zero workload.The factor 1-p suggests that roughly a fraction 1-p of the sample paths that start at U(0) K-will have the workload decreased to some O(1) value, before finally undergoing the large deviation to cross U(t) K.The remaining fraction p of sample paths that start at U(0) K will cross U(t) K in a short time period, before the server has had the chance to empty the system of the large workload.Since the first set of paths weight the mean escape time by an exponentially large amount, N(K-) will also be asymptotically exponentially large.Now consider (3.3) for general service time densities.By setting u-K-w and using a Laplace transform over u, we obtain from (3.3) the contour integral representation where '(s) f e 0 N(w) 1 / e(K -s , _ : _ ( w ) s [ N ( K ) 1/s] ds Br -SZb(z)dz is the Laplace transform of service time density, N(K)-N(K-), and the integration contour is a Bromwich contour on which Re(s) > 0. Note that the integrand has a double pole at s = 0.The constant N(g) is determined from (3.4), or, equivalently, from the condition N'(0) 0; hence Ks 1---fBr cls , -+ ( UN(K) N(K) .f seKs ds DEN(K)" (3.9) 'B-+ () When b(z) ge-u, b(s) p/(g + s) and it is then easy to show that (3.8) and (3.9) reduces to (3.7).
We next examine the asymptotics of N(w) as K.We use two independent approaches.First, we expand he integrals in (3.8) and (3.9) as K. Then we obtain the identical results by asymptotically analyzing (3.3) and (3.4) by perturbation methods.
The asymptotic behavior of the integrals in (3.9) is determined by the singularities of the integrands with the largest real part.The integrand in DEN(K) is analytic at s 0, but has a pole at s -a, where a satisfies the (real) transcendental equation / a + a, a > 0. 0 The existence of the solution follows from our assumption that the moment generating function Mb(O --f eWb(w)dw is analytic for Re(0)> 0, and the stability condition ,rn < 1.The 0 numerator in (3.9) has a pole at s-0.Thus, computing the corresponding residues in (3.9) we obtain where 1 DEN(K) ,Ii(a)_ 1; K-cx  (3.11) NUM(K) 1 p' ae-Ka II(a)--/ weaWb(w)dw" o (3.12) It follows that N(K-) is again exponentially large as K---,o.Now consider the integral in (3.8).
Since N(K) is exponentially large, we can approximate N(K)-1Is N(K)in the integrand, with an error that is exponentially small.For K-w u O(1), we cannot simplify (3.8) any further.For K-w--c, the asymptotic behavior of N(w) is determined by the simple pole in (3.8) at s 0. Using (3.8) and (3.11) thus yields N(K-) AIl(a)-1 Ka; An alternative approach to the asymptotics is as follows.We introduce the scaled variables x W 1 N(w) n(X) ,--, and note that e0 with Kc.In terms of these new variables, (3.3) becomes n(X) + J o and we use the boundary condition n'(0)-0.For 1-X >, e, we expand n(X)as n(X) C()[no(X n t-e.nl(X +...] (3.16) and extend the upper limit on the integral in (3.15) to + o.Assuming further that C(e)-oo as -0 (which is a very mild assumption since we eventually find that C() grows exponentially as ---0), we obtain from (3.15) to leading order the problem (m 1 1)n(X) 0; n(0) 0. (3.17) It follows that no(X is constant, and without loss of generality we set no(X 1.The expansion (3.16) breaks down near the exit boundary, where X 1 (to be more precise, 1-X O()).
Setting u (1 X)/ K-w, with n(X) g(w) C(e) [0(u) + hCl(u +...], we obtain from (3.18) 0 The function No(U is referred to as a "boundary layer" approximation, since it is to be valid in the thin region 1-X-O(e), where N(w) varies appreciably.By asymptotically matching the two expansions, we obtain o(u)l as uoc. (3.19) Then (3.18) is easily solved using Laplace transforms, and we get the contour integral representa- tion jo(U) 127ri-P / eus ds. (3.20) s-A + Ab(s) Br It also follows that the approximation n(X)C(e).ffo((1-X)/e is valid uniformly, for X [0, 1] (i.e., to [0, K]).It remains only to determine the constant C(e).
We denote the steady-state unfinished work distribution by p(w)dw-lim Pr[U(t)E (w, w + dw)], w > 0 (3.21) It satisfies the Takcs integrodifferential equation The solution is well-known to be A 1-p and p(w) A(1 p) / 1 -b(s) .eWSds"Br (3.24) Now we multiply (3.3) by V(W), integrate from w-0 to w-K, and use (3.4), (3.22) and (3.23).After some simplification, we are led to p(I)N(h') A + J p(w)dw 1-J p(w)dw.This identity is exact for all K > 0. As K-oc, (3.25) can be replaced by the asymptotic relation N(K) N(K-lip(K). (3.26) Now, from the perturbation expansion, N(K),,,, C(e)N0(0 -C(e)(1-p).The asymptotics of p(K) as K--oo are easily obtained from (3.24), as the singularity with the largest real part is at s--a (cf.(3.10)).Computing the residue at this pole yields This result, of course, agrees with (3.13), which was obtained by asymptotically expanding the exact solution.We summarize our results below.
Result 2: The mean time for U(t) to reach or exceed K in the M/G/1 queue is given by N Here a > 0 is the solution to (3.10), Ii(a) is defined by (3.12), and b(. is the Laplace transform of the service time density. We again see that N(w) is asymptotically constant, except near the exit boundary w-K.
4. Queue Length in the M/G/1 Queue For the M/G/1 model, N(t)is no longer a Markov process.Thus we let Y(t) be the elapsed service time of the customer presently being served and consider the joint process (N(t),Y(t)), which is Markov.We denote the steady-state probabilities by p(n y)dy lim Pr[N(t) n, Y(t) e (y, y + dy)] n > 1 p(0) -lim Pr[N(t) 0] (4.1) p(n) J p(n,y)dy, n >_ l. 0 These satisfy the balance equations pu(n, y) ip(n 1, y) [i + #(y)]p(n, y), n >_ 2 p(1,O) Ap(O)+ / p(2, y)tt(y)dy 0 p(o) J" p(, )()d 0 (4.) p(n, O) / p(n + 1, y)#(y)dy, 0 + n--1 0 Subscripts are used to denote partial derivatives.Here #(y) is the conditional departure rate given Y(t) y, and is related to the service time density b(.via Assuming the stability condition m I ( 1, (4.2) is easily solved using generating functions to give p(0)-1-p and 27ri s ---]ds, n _> 1 (4.4)C where C is a small loop about s-0.The integrand is analytic for sl <1 and at s-1.As n--<, the asymptotic behavior of the integral is determined by the simple pole at s-1 + a/A (cf.(3.10)).Thus, A[AII(a 1] \a -t- Note that for the M/M/1 model, a-#-A and then 11 -#A -2 and (4.6) reduces to the well known (exact) geometric distribution for this model.
We analyze the time for N(t) to reach N 0. We again define 7 by (2.1) and set The function T satisfies the backward equation ,XT(n + 1, y) + #(y)T(n 1, O) [A + #(y)]T(n, y) + Tu(n, y) (4.10) Note that, by definition, T(No, y) 0 for all y.For exponential service, #(y)= tt =constant and then (4.8)-(4.10)admits a solution in- dependent of y, which leads to the expression in Result 1.
We solve (4.8)-(4.10)by using generating functions.The final result is fb(t)( )(t )dt We examine T(n, y) asymptotically as No-Cx.In this limit, the expansions of the integrals in (4.11) and (4.12) are determined by the singularities of the integrands whose distance from the origin is minimal.For DENI(N0) the dominant singularity is s 1 +a/A, and for NUMI (No)   the dominant singularity is s-1.
(4.25) Y In order for expansions (4.17) and (4.20) to be consistent, we must have 0(rn, y)--l as rn--oe for each y.As rnoc, the expansion of (4.25) is determined by the simple pole at s-1, which is the singularity of the integrand that is closest to the origin.Noting that the residue is independent of y, the matching condition implies that 0(1, 0) 1 p 1 Am 1. (4.26) It remains to determine the constant D(5).To this end, we use the steady-state balance equations (4.2).We multiply (4.8) by the steady-state probability p(n,y), integrate from y 0 to y ec, and sum from n 1 to n N o 1. Integrating some of the integrals by parts, shifting indices in the summations, and using the fact that p(n, y) satisfies the system of equations (4.2), we eventually obtain We evaluate this identity as Nooe.
an error that is exponentially small.
From the perturbation method, we found T(N o 1,0),-D( 5)ro(1,0 D(5)(1-p) and p(N o 1,0) can be evaluated using (4.5).solving (4.27) for 0(5) we obtain AIl(a)-i [ a)N0 -1 D(5) (1 p)2a , 1 +which agrees with (4.14).We summarize the main results.lesult 3: The mean time for N(t) to reach N O in the M/G/1 queue is given by fCb(t)e(s 1)(t Y)dt y The right side of (4.27) is asymptotically equal to 1, with that Thus, For No--oo asymptotic expansions for T(n, y) are (a) No--+oo When b(y) #e-u, b(A As) #/[# -t-A As] and then T is independent of y.By evaluat- ing the contour integrals, we can easily show that Result 3 reduces to the corresponding parts in Result 1.We again observe that T is exponentially large and approximately constant, except near the exit boundary.Near the exit boundary T is still exponentially large, but now depends upon y and m-N O -n.The expression for T(n, y) in Result 3 applies for all n >_ 0. Note that when n _> No, the contour integral vanishes for all y >_ 0, as now the integrand is analytic at 8--(} 5. Queue Length in the GI/M/1 Queue We consider the GI/M/1 model with service rate # and interarrival time density a(.).The process N(t)is not Markovian, but the joint process (N(t),Z(t)), where Z(t) denotes the elapsed time since the last arrival, is Markov.Assuming the stability condition p-(#M1)-1 < 1, we define the steady-state probabilities by p(n,z)dz lim Pr[N(t) n Z(t) e (z,z + dz)] n > 0 (5.1) p(n)--/ p(n,z)dz, 0 These satisfy the balance equations n>O. (5.2) Here A(z) is the conditional arrival rate given Z(t)z; it is related to the interarrival time densi- ty a(.via b,-J e(b*-l)za(z)dz, O < b, < l. 0 The marginal queue length probabilities are p(rt) )z/-,(b,)n-1, n )_ 1; p(0) 1 p 1 #M 1" set To compute the time until N(t) reaches No, we again define the stopping time r by (2.1) and T(n,z) E(7 IN(0) n, z(o) z).
We next evaluate T in the limit N0oe.The integrand in (5.13) now has a unique pole inside the unit circle Is < 1, at s-b. (eft (5.6)).For large No, this pole determines the asymptotic behavior of T(1,0), and we obtain For n, the contour integral in (5.12) also has a pole at s b,, and grows exponentially.For No-n this integral is smaller asymptotically than T(1,0) and we obtain T(n,z) T(1,0) for N0, N o-n.For N 0-n-m-O(1), the integral in (5.12) is of the same order of magnitude as T(1, 0).Evaluating the residue from s b,, we thus obtain -b(t-Z)a(t)dt An alternate approach to the asymptotics is to analyze (5.9)-(5.11)by the perturbation method.We set n x/5, 5 N o-1 with T(n,z) t(x,z) D(5)[to(x,z + 5tl(X,Z +...].Using these scaled variables in (5.9) and expanding for 5--0, we obtain to, z(X,Z + )(z)[to(x,O to(x,z)] 2.1t0 0 (5.16) "1tl Ox[to(X, z) (Z)to(X 0)].
(5.17) Solving these equations in a manner similar to that used to analyze (4.18) and (4.19), we even- tually find that to(x,z -to(x -1.Then the boundary condition (5.10)is also asymptotically satisfied.
Result 4: The mean time for N(t) to reach N o in the GI/M/1 queue is given by f For N0-oc asymptotic expansions for T(n, z) are (a) N0--,oo Here b-#(1-b,) where b, is given by (5.6), Jl(b)is defined in (5.15), and if(-)is the Laplace transform of the interarrival time density.
The overall asymptotic structure is similar to that for the M/M/1 and M/G/1 models.For the GI/M/1 model the structure of the solution in the boundary layer (result in part (b)) is simpler than the corresponding M/G/1 result.We also note that if N o 1, the exact result be- f a(t)dt z which is just the mean residual life on the renewal process that governs the arrivals.Now the first passage time measures the first arrival time to an empty system.
Workload in the GI/G/1 qnene We consider U(t) in the GI/G/1 model.We again let Z(t) be the elapsed time since the last arrival and analyze the Markov process (U(t),Z(t)).The steady-state distribution of this process is well-known to be equivalent to solving a Wiener-Hopf problem.Since our approach, which uses the supplementary variable technique, is different from that used in most books on queueing theory, we briefly give the details of the computation of the steady-state distribution.Our main goal is the computation of the time for U(t) to reach or exceed K.We can compute this exactly for the GI/M/1 model, but not for the general GI/G/1 case.For the latter case, we derive asymp- totic results by using the perturbation method that we outlined in the previous sections.This method requires that (i) we know the tail behavior of the steady-state distribution for the joint process (V(t),Z(t)) and (ii) we use the balance equations satisfied by this distribution.
Im( 0) A sketch of the integration contours in the complex w-plane.
A-f A(z)dz-1-p, as is well-known.0 From (6.16), it is easy to show that We now evaluate the tail behavior of p(w,z).Let c be the unique positive solution of the (real) transcendental equation e-CZa(z)dz eCyb(y)dy 1, c > 0.
(6.17) 0 0 Then the integrand in (6.15) has a pole at s--c, and this pole determines the asymptotic behavior p(w,z) as w-cx.We have [./'z 1 (1 p)cIo(c)ro(C )e-CW w--oc (6.18)   p(w,z -1]dw}. (6.20)The contour in (6.20) may be taken along the imaginary axis, with an indentation about w-0 in the right half-plane.Note that equation (6.17) may be compactly written as Jo(C)Io(c -1. the marginal density of U(t) has the tail V(W) J V(W, z)dz (1 p)(I o 1)roe ,, woo.
We compute the mean time for U(t) to reach or exceed K. Let r, be as in (3.1) and set N(, z) Z(, U(0) , Z(0) ).
The latter may be replaced by the equivalent "reflecting" condition Nw(O +, z) 0, z > 0. When a(z)-Ae -'xz we have A(z)-A-constant and then N(w,z)-N(w)so that (6.24)-(6.25) reduces to (3)3)-(3.4).Note that an arrival causes the clock on the renewal process Z(t) to be reset to zero and thus the second argument in N(.,. in (6.24) inside the integral is zero.Also, from the definition (3.1), we see that N(w,z) for w >_ K for all z.From now on, we will use N(K,z) to mean N(K-,z).
We can construct the (exact) solution to (6.24) for the GI/M/1 model with b(y)-#e-"Y.
It is then easy to evaluate this expression asymptotically as K-.cx.Below we give only the final results.
Result 5: The mean time for U(t) to reach or exceed K in the GI/M/1 queue is given by N(w,z) #M1/ # 3(s)esK Br For K---oc, asymptotic expansions for N(w, z) are (a) K--oo, g-w---,oo Here b, Jl(b) and if(.are as in Result 4. On the contour Br we have b < Re(s) < #.
Note that the first formula in Result 5 applies for 0 <_ w < K and all z.(w,z), we clearly have N(w,z)-0 by definition.

For other values of
We now consider the general GI/G/1 model.While we have not been able to solve (6.24)- (6.25) exactly in this case, the perturbation method in the previous sections can be readily extended to the general model.We begin by establishing an identity between the solution of (6.3)-(6.7)and N(w,z).As before, we multiply (6.24) by p(w,z)and integrate the resulting expression over the strip 0 < w < K, 0 < z < oc.To this we add A(z) times equation (6.25), integrated from z 0 to z oc.After some integration by parts and use of equations (6.3)-(6.6),we obtain 0 0 0 0 This identity generalizes (3.25).The right side is asymptotically equal to 1.To evaluate the left side as K--+oo, we use (6.18) with w K. Then we must know the asymptotic expansion of N(K-,z) which corresponds to computing the mean first passage time for an initial workload just below the exit boundary.We also note that for the GI/M/1 model, p(w,z) has the simple form _be-b p(w,z) -M-71 Zexp-( 7)dr e 0 bw while N(w,z)is quite complicated (cf.Result 5).By using (6.27) and Result 5 to evaluate the first integral in (6.26), we have verified that (6.26)is indeed satisfied.
Also, it is easy to show that as sl in the corresponding half-planes, we have exp[r.(s)](const.)/sand exp[ (s)] (const.)/s.It follows that the second integral in (6.40) vanishes, as can be seen by closing the Bromwich contour in the right half-plane, where the integrand is analytic.The first integral in (6.40) can be evaluated by closing the contour in the left half-plane.In this region, the only singularity is the simple pole at s -c.Hence, we have -czN(O,z)dz-(M1 ml)exp[I',( c)] (M 1 rrt I )r0(c).(6.41) Using (6.41) in (6.37) we solve for C(), and then the asymptotic expansion of N(w,z)is completely determined.We summarize the final results below.
Br + On the contour Br we have 0 < Re(s)< e0, and Br+ may be taken as the imaginary axis with an indentation about w-0 in the right half-plane.
The overall structure of N for the GI/G/1 model is similar to that we obtained for the special cases.The numerical evaluation of the constant C requires that we evaluate (numerically) the contour integral in r0(c).This can be done analytically, in closed form, if a(.) and b(.) have rational Laplace transforms.
For the GI/M/1 model, we have c b, b(y) #e-t*v, r,() ]/(1 p) and (IiJ0-IoJ1)/(bIoPo) (1-P)M1/b.Then Result 6 reduces to parts (a) and (b) in Result 5.The contour integral in part (b) of Result 6 may now be explicitly evaluated, as the integrand has but two poles in the region Re(s) < 0, at s 0 and s -b.The residues at these poles correspond to the two terms in Result 5, part (b).
7. Queue Length in the GI/G/1 Queue We compute asymptotically the mean time for N(t) to reach N o in the GI/G/1 model.We consider the Markov process (N(t),Y(t),Z(t)), where Y(t)is the elapsed service time and Z(t)is the age on the renewal process that governs the arrivals.The perturbation method requires that we know the tail behavior of the joint steady-state distribution of the 3-dimensional Markov process.The method also uses the balance equations satisfied by this distribution function.
We hence define (7.4)Note that we have decomposed the probability that N(t)-1 into two pieces (7.2) and (7.3).
The tail behavior as n--<x: is easily obtained by shifting the Br contours in (7.28) and (7.29) to the left, past the pole at 0--c.By computing the residues at this pole we obtain p(n, y, z) exp (t)dt #(') Here c, F0, I j, J j are as in Result 6.
Zfo(m,y,z)--l as rn--c, for all y and z.
summarize our results below.
We have thus determined D(6) and we Pesult 7: queue are given by (a) No--+oo N O n---+oo [/0 1 T(n, y, z) 1 l(c)Jo(c Io(c)J 1(c), eCUb(y Asymptotic expansions for the mean time for N(t) to reach N O in the GI/G/1 The various quantities are defined in Result 6.On the contour Br, 0 < Re(0) < 1, where fi'(-0) is analytic for Re(0) < 1" It is easy to show that Result 7 reduces to the asymptotic results in Result 3 for the M/G/1 model and to those in Result 4 for the GI/M/1 model.For the M/G/1 model, we have (-0) 1/(A-0) and then the contour integral in G(m,z) may be deformed into a small loop about 0 A. For the GI/M/1 model, the integrand in a(m,z) has only two poles (at 0 0 and 0 -c -b) in the half-plane Re(0) < el" 8. Discussion and Numerical Results We have computed the mean time needed for large workloads and large queue lengths to develop in the GI/G/1 queue.If either the arrival times or service times are exponentially distri- buted, we have derived exact contour integral representations for the mean first passage times.From these, the asymptotic results are easily obtained.For the general model, we have computed the mean time asymptotically by using singular perturbation methods to analyze the appropriate backward equation(s).
The asymptotic method requires that we know explicitly the tail behavior of the steady-state distribution of the queue length or workload.For the GI/G/1 model this may be easily obtained from the Wiener-Hopf theory.The method also requires that we carefully analyze the first passage time for initial conditions that are close to the exit boundary (i.e., N(0) _ No, U(0) _ K).
Using asymptotic methods, we have shown that this involves solving a second Wiener-Hopf problem.However, the second problem is very similar in structure to the computation of the steady-state distribution.
In principle, it should be possible to extend the asymptotic method to other models, such as the m-server GI/G/m queue.For this model, the tail exponent in the steady-state distribution(s) is easily computed, but it is very difficult to compute the constant in the asymptotic relation(s).
The latter would seem to require that we have some exact representation of the steady-state distribution(s).This is known for the GI/M/m model but not for the M/G/m model.Indeed, the analysis of even the M/G/2 queue is quite difficult (see Hokstadt [7]; Cohen [3]; and Knessl,  Matkowsky, Schuss and Tier [14]).Thus, while the asymptotic method reduces the first passage time problem to two simpler problems, the solution of the latter may itself be difficult.
We next discuss the numerical accuracy of the asymptotic results.We would like to get some idea as to how large N O and K must be before there is good agreement between the exact and asymptotic results.Let N and T be the exact values of the mean first passage times, and let N asy and T asy be the asymptotic formulas we give in Results 1-7, for initial conditions not close to the exit boundary.From Results 1-5, we can easily obtain the estimates N(w,z)/Nasy(W,Z 1 + O(EST), T(n, y, , z) 1 + O(EST),  where EST denotes terms that are exponentially small in the appropriate limits (and thus smaller than any power of K-1 or N 0-1).The estimate (S.1) is obviously true for the M/M/1 model (cf.Result 1 and (3.7)).It can easily be shown to hold for the M/G/1 and GI/M/1 models by estimating the various contour integrals in Results 2-5.We conjecture that (8.1) is also true for the general GI/G/1 case.This suggests that the asymptotic results should be highly accurate, even for moderate values of K and N 0.
In Tables 1-3 we compare the asymptotic and exact formulas in Result 2 for the M/E2/1 queue which has service time density b(y)-(2#)2ye 2uy with m I l/it.We set A-1 and in Tables 1, 2, and 3 we have #-4,2 and 4/3; respectively.The respective traffic intensities are thus p-0.25, 0.50, 0.75.We also set w-0 and compare the exact result to the asymptotic result N(0) C in part (a) in Result 2. The contour integrals in the exact answer are easily evaluated since the various integrands have only two or three poles.Also, the solution to (3.10) is a 1/2[4#--A--V/2+8#A]-[4-p-V/p2+ 8p].
In Table 1 we consider values of K in the range 5 _ K _ 10.The exact and asymptotic answers agree to 6 decimal places.Since the traffic intensity is quite small, the mean first passage time is very large (about 1010) even when K 5.
Tablel: A-l, #-4 In Table 2 we increase p to 0.5.There is still good agreement between exact and asymptotic answers, with the worst maximum relative error being less than 1% when K-5.When K-10, the two answers agree to 5 decimal places.In Table 3, we further increase p to 0.75.Now the error is unacceptable (about 40%) when K-5, but decreases to about 5% when K-10 and to under 1% when K-15.As p--.1, the asymptotics are no longer valid.
Table3: A-l, #-4/3 In Tables 4-6, we compare the asymptotic and exact formulas in Result 3 for the M/D/1 queue, which has service time density b(y)=5(y-1/#).We consider initial conditions (n,y) (0,0) and use the asymptotic result T(0,0),, D in part (a) in Result 3. We set A 1 and Tables 4, 5, and 6 have #=4(p=0.25),#--2(p--0.50),and #--4/3(p--0.75);respectively.The exact answer was obtained by using symbolic methods to compute the residues at s 0 for the various contour integrals.To get dependable numerical answers when p 0.25, it is necessary to do the calculation using 30 digits of precision.Now a >0 satisfies A + a Aexp(a/#), which we solved numerically.In Tables 4-6, we consider values of N O in the range 5 <_ N o < 15.When p 0.25, Table 4 shows that we get agreement within 1% even when N o 5.When N o 15, the asymptotic and exact answers agree to 6 decimal places.In Table 5, we increase p to 0.5.The error is about 3% when N O -5 but decreases to under 1% for N O >_ 7; when N o 15, we get agreement to 6 decimal places.6 has p 0.75.When N O 5, the error is an unacceptable 50%, but decreases to about 3% when N 0-10 and to about 0.3% when N 0-15.Throughout Tables 4-6, the error is always under 1% for N0(1p) >_ 4.