TRANSIENT EXPONENTIAL S TEADY-'STATE - ERLANG QUEUES AND S IMULATION by W. David Kelton Department of Industrial and Operations Engineering Technical Report 84-7 April 1984

Transient Exponential-Erlang Queues and Steady-State Simulation W. David Kelton Department of Industrial and Operations Engineering The University of Michigan ABSTRACT The transient probabilistic structure of arbitrarily initialized M/E /1 and m~~~~~~~~~~~~~~~m Em/M/l queues is derived in discrete time. Computational algorithms for obtaining the required probabilities are provided, and their application in calculating a variety of system performance measures is illustrated. The results are used to investigate the question of initializing simulations of systems such as these in order to promote rapid convergence to steady state, if that is the object of the simulation. These results are consistent with earlier studies for transient queueing systems, such as the M/M/s, but allow greater flexibility in specification of interarrival or service-time models inherent in the Erlang distributions.

1 1. INTRODUCTION Analytical results for transient characteristics of queueing models are not as widely available as are steady-state results, but are directly useful for studying the finite-time properties of systems accurately represented by such models. There are several additional indirect reasons for having exact transient results in simulation applications and methodological research: * To serve as a controlling system in the external control variates technique for variance reduction (see, for example, Gaver and Shedler [4]). If a terminating simulation is to be done for a system resembling a simpler system with known transient behavior, this second system's output from a simulation would be expected to be strongly correlated with that of the system of interest (assuming the use of common random numbers for both simulations), leading to large variance reductions. The wider the class of analytically tractable models from which to choose, the greater the similarity possible, leading to better variance reductions. * To serve as benchmark models on which to test techniques for controlling startup bias in steady-state simulations (see Gafarian, Ancker and Morisaku [3], Schruben [18], or Kelton and Law [9]). * To serve as benchmark models for studying alternative methods for initializing simulations aimed at estimating steady-state parameters (see Wilson and Pritsker [21] or Kelton and Law [10]). This final reason served as the main motivation for the present paper, and is treated in more detail later. Available transient results for queueing models may be classified according to whether the time measure is continuous (real time) or discrete

2 (indexing by customer number). Continuous-time results for M/M/1 and M/M/s queues appear in Morse [13], Saaty [17], Rothkopf and Oren [16], Clark [2], van Doorn [19], Whitt [20], Halfin and Whitt [6], Pegden and Rosenshine [15], Grassmann [5], and Odoni and Roth [14]; see also references in these papers. Whereas continuous-time analysis is useful for studying questions such as the queue length at a particular time or the experience of a customer who might arrive at a certain point, discrete-time results are more relevant if we are interested in the experience of, say, the nth arrival to a system or the state just after the nth arrival. Such is the case in many simulations, where one typically focuses on estimating properties of customers' delays in queue, other continuous-time parameters (e.g., mean queue length) being estimable indirectly from delay statistics via conservation equations (see Carson and Law [1]). Papers dealing with discrete-time transients of queueing systems include Heathcote and Winer [7], Morisaku [12], and Kelton and Law [10]. In this paper we extend the body of discrete-time transient results to include M/Em/1 and E /M/l queues, where E denotes an m-Erlang distribution. m m m Also, our results permit arbitrary initial states of the system in terms of the number of Erlang stages present; this allows a numerical evaluation of the effect of alternative initial conditions on the nature of convergence to steady state, a general question of interest in simulation aimed at estimating steady-state parameters. In Sections 2 and 3 the analytical results for the two classes of models are derived, with specific algorithms for computational implementation and application. Section 4 reports on numerical evaluation of these results to address questions of initialization to promote rapid convergence to steady state in simulation experiments. Some conclusions are drawn in Section 5, and the Appendix contains proofs of the results in Sections 2 and 3.

3 2. THE M/E /1 QUEUE m Let X be the exponential arrival rate, u be the m-Erlang service rate (where a complete service time is composed of m consecutive exponential stages each at rate mn), and let p = X/g; it is not necessary to assume that p < 1 for any of the results of this paper, so that the rate of explosion of these queues could be studied in the case that no steady state exists. For n > 1, let tn be the time of arrival of the nth customer to the system. The "method of stages" analysis of this system proceeds by using as a state variable the number of exponential service stages (rather than customers) present in the system, i.e., if there are c customers present in the system (including the one in service, if any, so c > 0) and the customer in service is in the dth of his m service stages, the system state would be cm - d + 1. 2.1 Mass Functions The embedded discrete-time process used is defined (following Morisaku [12]) as X = the number of service stages present in the system at time t, including the m stages arriving at time t, for n > 1. Letting k be the number of stages in the system at time 0 (k > 0) and noting that the range on X is then {m, m+l,..., k+nm}, the probability mass function of Xn, conditional on k, is Pk(n,i) = P(Xn = i I X0 = k), where X0 is the initial number of stages. The first arrival occurs at time tl, which is exponentially distributed at rate X; thus, t1 > 0 and the first arrival finds at most k service stages already present. The following three propositions (proved in the Appendix) are sufficient for computation of the Pk(n,i)'s.

4 Proposition 1. For n > 1, [p/(p + m)]n if k > 1 Pk(n,k+nm) = n[p/(p + m)] if k = O Proposition 1 represents a boundary contition i.e., that X is at its maximum. n The following proposition establishes another boundary condition concerning the mass function of the system state just after the first arrival. Proposition 2. For k > 1 and m + 1 < i < k + m - 1, Pk(l,i) = [(m/(p + m)] [p/(p + m)]. Note that k = 0 is excluded in Proposition 2, but in this case X1 = m almost surely. The following result is the main recursion, and is true regardless of whether k is zero or positive. Proposition 3. For n > 2 and m + 1 < i < k + nm - 1, k+(n-l)m Pk(n,i) = [p/(p + m)] Z [m/(p + m)] -imp (n-lj). j=max(i-m,m) Note that the case i = k + nm, omitted in Proposition 3, is covered by Proposition 1, and that for the case i = m, we obtain k+nm Pk(n,m) = 1 - Z Pk(n,i). i=m+l 2.2 Computational Algorithm The formulas in Propositions 1 - 3 above can be used directly to compute the mass function of X, but it is possible to manipulate them into the following more efficient algorithm. The procedure must be invoked in the order n = 1, 2, 3,..., and takes as input the values of n, m, p, k, and Pk(n-l,i) for m < k + (n-l)m (unless n = 1); the output is Pk(n,i) for m < i < k + nm. k - -

5 procedure M/E /1 [n, m, p, k, Pk(n-l,i); Pk(n,i)] if n = 1 then if k = O then Pk(l,m) - 1; return else h 4- m/(p + m); Pk(l,k+m) 4 p/(p + m); s 4 Pk(l,k+m) for i + k + m - 1 to m.+ 1 by -1 do Pk(l,i) 4 hPk(l,i+l); s - s + Pk(l,i) end do Pk(l,m) - 1 - s; return end if else a - p/(p + m); b 4 1 - a; Pk(n,k+nm) 4 an1 s 54 if k > 0 then Pk(n,k+nm) - aPk(n,k+nm) for i - k + nm - 1 to 2m by -1 do Pk(n,i) - aPk(n-l,i-m) + bPk(n,i+l); s 4 s + Pk(n,i) end do for i 4 2m -1 to m + 1 by -1 do Pk(n'i) 4 bPk(n,i+l); s - s + Pk(n,i) end do Pk(n,m) - 1 - s; return end if end procedure M/E /1 The principal storage requirements of this algorithm are two vectors of maximal length k + (n* - l)m + 1, where n* is the largest value of n for which the mass function of X is desired; the first vector holds Pk(n-l,i) and the second holds Pk(n,i). After the nth invocation, the previous mass function

6 Pk(n-l,i) is replaced by the newly computed Pk(n,i) for input into the (n+l)st invocation. Translation of the algorithm into any structured language should be immediate; the computations in Section 4 were carried out using VS Fortran, a subset of Fortran 77. It was found that double precision (64-bit) was necessary to avoid buildup of roundoff in the recursive computations by the time large values of n (e.g., 500) were reached. 2.3 Applications Given the mass function Pk(n,i) of Xn, it is possible to develop simple formulas for several standard measures of queueing performance. Immediately, the expectation and cumulative distribution function of Xn are given as k+nm Ek(Xn) = iPk(n,i) n i=m and LxJ Pk(X < x) = 2 Pk(n,i), n - i=m where LJ denotes the greatest integer function, Ek and Pk respectively denote the conditional expectation and probability measure conditioned on the event X0 = k, x is any real number, and an empty sum is defined as zero. More easily interpreted than X (the units of which is service stages) is Yn, the number of customers present in the system just after time tn. The range on Yn is thus all integers between 1 and p'(k,m,n) = Lk/mJ + n + 1 inclusively, and the relation between Y and X is given by ~~~~~n n Y = p if and only if pm < X < (p + l)m - 1, if 1 < p < p'(k,m,n) - 1 and Y = p'(k,m,n) if and only if p'(k,m,n)m < X < k + nm. Thus, letting Qk(n,p) denote the mass function of Y, we get k~ n

7 (p+l)m-l z P,(n,i) if 1 < p < p'(k,m,n) - 1 i=pm Qk(n,p) = k+nm | z Pk(n,i) if p = p'(k,m,n) i=pm The expectation, for example, of the number of customers (not service stages) in system just after tn is thus p'(k,m,n) Ek(Yn) = 2 pQk(n,p) p=l and the cumulative distribution, variance, etc. of Y could be found n similarly. Further, if Z denotes the number of customers in queue just after n t, then n Y - 1 if Y > 1 from which the mass function, expectation, etc. of Zn can be found using the probabilities Q (n,p). As a final application that is of most interest to simulation, let D be the delay in queue (excluding service time) of the nth arriving customer. If X = m, then customer n arrives to find the system empty, so D =0. However, n n if X = i > m, then at least one service stage remains at the time of the nth arrival, so customer n is delayed in queue for the remainder of the inprogress service stage, plus i - m - 1 additional complete service stages. By exponentiality of service stages, the remainder of the in-progress stage also is exponential, so that Dn is the sum of i - m independent exponential service stages, each at rate mg, i.e., Dn is an (i-m)-Erlang random variable with mean (i-m)/(mu). The cumulative distribution function of Dn is then obtained by conditionitiing on Xn,

8 k+nm Pk(Dn< x) = Pk(D < I X = i) Pkn,i) i=m k+nm = Pk(n,m) + Z G.i(x;m) Pk(n,i), i=m+l where Gq(x;n) is the q-Erlang cumulative distribution function with mean q/n, q-1 G (x;t) = 1 - e X (7x)3/j! q - j=0 for any x > O. Similarly, the expected delay in queue of the nth customer is k+nm Ek(Dn) = [l/(m)] Z (i-m) Pk(n,i). i=m+l Section 4 discusses results of evaluating Ek(Dn) over a range of system parameters. Finally, the distribution of the total system wait of the nth customer, W = D + S, where S is an independent m-Erlang service time, can n 3. THE E /M/1 QUEUE m As before, let X and g denote the arrival and service rates, and define p = X/1. Here, however, service times are simply exponential, and we think of an arrival as occurring in m consecutive independent exponential stages, each at rate mX, with exactly one customer in some stage of arrival at all times; see, for example, Kleinrock [11]. As soon as an arriving customer finishes the mth stage of arrival (and thus physically arrives to the system), another customer begins the first stage of his arrival. The state of the system is the number of exponential arrival stages present, counting m for each customer physically present. Thus, if c customers are physically present (including the customer in service, if any) and the customer currently in the arrival process is in the dth stage of arrival (1 < d < m), the system state would be cm + d - 1.

9 3.1 Mass Functions Let t. be the time of the jth arrival stage completion, for j > 1, and let X be the number of arrival stages present just after time t.. If there are k (k > 0) arrival stages present at time 0, then the time of the nth (n > 1) customer arrival (physical) to the system is tmk+mk/m. Note that at each tj the system state rises by 1, and at the time of each service completion the state falls by m. For j > 1 and k > 0, let Pk(ji) = P(Xj = i | X0 = k), (1) where X0 is the number of arrival stages present at time 0. Again, tl is not zero, but is exponential with mean l/(mX). To determine the range of Xj, first note that it is maximally j + k, an attainable bound in the event that no departures occur in [0, tj]. At the other extreme, the minimal possible value for X. occurs if there are no customers in the system just after time t. other than the one completing an arrival stage at that time. If j + k is divisible by m, then a customer physically arrives at time tj, so the minimal X. is m; otherwise, the minimal X. is in {1, 2,..., m-l}. In either case, the minimal X. is attained if the maximal number of departures in [0, tj] occurs, which is Lj + k - l)/mj. Finally, since each customer departure drops the system state by m, it is not possible for Xj to take on all integral values between its minimum and j + k. Thus, the general range on X. is {j + k - fm: 0 < f < L(j + k - 1)/mj}, (2) the range of i over which the mass function in (1) must be computed. Here, f represents the number of service completions occurring in [0, tj]. The following three results, proved in the Appendix, are sufficient for calculating the mass function in (1) over the range in (2):

10 Proposition 4. For j > 1, 1 if k < m and j < m - k Pk(j,j+k) = [mp/(mp + )]+km if k < m and j > m - k [mp/(mp + 1)]j if k > m Proposition 5. For k-> m and 0 < f < lk/mj, Pk(l,k+l-fm) = mp/(mp + 1)f+. Proposition 6. For j > 2 and 1 < f < Lj + k - 1)/mj - 1, ffPk(j,j+k-fm) = [mp/(mp + 1)] z [l/(mp + 1)] gPk(j-l,j-l+k-gm). g=0 Note that the analogue to Proposition 5 for the case k < m is covered by the first two branches of Proposition 4, since in this case the range of f in the set (2) for j = 1 is simply f = 0; thus Pk(l,k+l-fm) = Pk(1,1+k), which is in Proposition 4 for j = 1. Similarly, if f were 0 in Proposition 5, then Proposition 4 would instead apply. Finally, as a result of Proposition 5 we obtain Lk/mJ -1 Pk(l,k+l-Lk/mj m) = 1 - 2 Pk(l,k+l-fm), f=0 completing calculation of the mass function of X1. In Proposition 6, the cases j = 1 or f = 0 are covered by Proposition 4, and we obtain L j+k-l)/mj -1 Pk(jj+k- Lj+k-l)/m m) = 1 - Pk(j,j+k-fm). f=0 3.2 Computational Algorithm As for the M/Em/1 model, we can derive from Propositions 4 - 6 a more efficient recursive algorithm for computing the values for the Pk(j,i) mass function for i in the set (2). The following procedure is entered in the

11 order j = 1, 2,..., takes as input j, m, p, k, and the values in Pk(j-l,i) (unless j = 1) and returns the values Pk(j,i) for i in (2). procedure E /M/1 [j, m, p, k, Pk(j-l,i); Pk(ji)] a 4 mp/(mp + 1); b 4- 1 - a if j = 1 then if k < m then Pk(l,k+l) - 1; return else Pk(l,k+l) 4 a; s + a for f - 1 to k/mi - 1 do Pk(l,k+l-fm) 4 bPk(l,k+l-(f-l)m); s s + P(l,k+l-fm) end do Pk(lk+l-Lk/mJ m) 4 1 - s; return end if else if k < m then if j < m - k then Pk(jk+j) 4 1 else Pk(jk+j) a -m end if else Pk(j,k+j) - aj end if s 4- Pk(j,k+j) if (j + k - l)/mj = 1 then go to LAST if k > m or j > m - k + 1 then

12 Pk(j,k+j-m) - bPk(j,j+k) + aPk(j-l,j-l+k-m) else Pk(j,k+j-m) - a[b + Pk(j-l,j-l+k-m)] end if s 4 s + Pk(j,k+j-m) for f - 2 to L(j + k - l)/m - 1 do Pk(j,j+k-fm) - bPk(j,j+k-(f-l)m) + aPk(j-l,j-l+k-fm) s - s + Pk(j,j+k-fm) end do LAST: Pk(j,j+k- (j+k-l)/m m) 4 1 - s end if return end procedure E /M/l m As before, two vectors of storage are required, one for the previous j and one for the current j. Also, the length of the vector for storing Pk(j,i) is not j + k (the maximal i), but is Lj + k - l)/mJ + 1, the number of values in (2). 3.3 Applications Given Pk(j,i), several measures of performance of the transient E /M/l queue m are possible which parallel those for the M/E /1 described in Section 2.3. Again letting D be the delay in queue of the nth (physically) arriving customer, recall that tj(n) is the time the nth customer physically joins the system, where j(n) = nm - k + mLk/ml; thus, we are concerned only with the probabilities Pk(j(n),k+j(n)-fm) = Pk(j(n),(n+Lk/mJ -f)m), for 0 < f < n + Lk/mJ - 1. If f = n + Lk/mJ - 1, then the arrival at tj(n) finds the system empty, so D = 0. If f < n + lk/mJ - 1, then this arrival finds n

13 n + Lk/ml - f - 1 other customers in the system, each of which still requires an independent exponential service with mean 1/1, so D is in this case an (n-[k/mJ -f-l)-Erlang random variable with mean (n + Lk/n! - f - 1)/1. Thus, given x > 0, the cumulative distribution function of Dn is n+ k/mj -2 Pk(Dn < x) f=Z Gn+Lk/J -f-l(X; ) Pk(j(n),k+j(n)-fm) + Pk(j(n),nm) and the expected delay in queue is n+ Lk/mJ -2 Ek(Dn) = (1/O) (n+Lk/mj -f-l) Pk(j(n),(n+Lk/mJ -f)m), f=0 where an empty sum is taken to be zero. Again, properties of the total system wait W of the nth customer may be derived from these expressions. As a second application, let Y be the number of customers present in the n system just after tj(n). Thus, Y = Xj/m, and is a discrete random variable on {1, 2,..., n+[k/m } with mass function Pk(Yn=i) = Pk(j(n),im), enabling computation of its cumulative distribution, moments, etc. 4. EVALUATION OF ALTERNATIVE INITIALIZATIONS FOR SIMULATION The principal motivation for this study was to examine the effect of alternative simulation initialization policies on the nature of convergence of simulation output to steady state, for models similar to those analyzed above. In a steady-state simulation, the goal is to estimate some property of the steady-state distribution (assuming it exists) of the output stochastic process, often its expectation. One of the most difficult problems facing an analyst in this case is choosing initial conditions which are in some sense representative of steady state. Of course, it is impossible to do this exactly since knowledge of the steady state distribution would be required. Results such as those obtained above can be used to address the question of

14 prudent choice of initial condidtions. In [10] we carried out such an investigation for M/M/s queues, and found that for that system, it may be wise to initialize in something other than the popular empty-and-idle state to reduce the length of time required for nearsteady-state conditions to be attained. One limitation of that study was that all interarrival time and service time distributions were assumed to be exponential. The results of this paper, while limited to single-server systems, allow a much richer class of distributions to be used, the Erlang, providing more realism; an Erlang assumption for service-time distributions is especially attractive, since it appears that for many processes an Erlangshaped histogram arises from service time data. For [10] as well as this paper, the critical point is that analytical transient results were obtained allowing an arbitrary initial state specification, so we can evaluate the resulting functions for various choices of initialization and observe convergence behavior. Values for Ek(Dn) were computed for both the M/Em/1 and Em/M/l systems k n m m across a range of parameter values: p = 0.5, 0.8, and 0.9 (always setting X = 1), and m = 2, 4, 8, and 9; k was chosen as described below. Figures 1 and 2 show plots of Ek(Dn), as functions of n, for various initialization schemes, for m = 4 and p = 0.9. The number of customers (not stages) physically present in the system initially is c, shown next to the corresponding curves. For simplicity, we assumed for the M/Em/1 model in Figure 1 that the customer in service (if any) was just beginning his first service stage at time 0; similarly, for the E /M/l model of Figure 2, the customer in the arrival m process at time zero was assumed to be at the start of his first stage of arrival. In both figures, the dashed line is at the expected steady-state delay in queue, found from Kleinrock [11] for Figure 1 and from tables in

15 (%4 C' 0 V) 0 w) a3 9, 100 200 300 400 n 50o Figure 1. Ek(Dn) for the M/E4/1 Queue with p = 0.9 and k = 4c. j~k n 4*

16 o 0 VI O 4) Cl LJ mo "b 200 300 400 500 n Figure 2. Ek(Dn) for the E4/M/1 Queue with p = 0.9 and k = 4c.

17 Hillier and Yu [8] for Figure 2. In both cases, Ek(Dn) converges to the steady-state expected delay monotonically from below if c = 0, but for other choices of c the approach may be nonmonotonic (e.g. c = 5 in both figures). Such behavior has been noticed for the discrete-time M/M/s queues in [10], as well as in continuous-time transient results by Grassmann [5] and Odoni and Roth [14]. As pointed out in these last two papers, the initial decrease in the curves (even though they begin below the steady-state mean) is attributable to the fact that c > 0 implies that the server is initially busy (which is not necessarily the case in steady state), increasing the probability of a downward state transition with respect to its steady-state value. Plots of other cases exhibited similar behavior; sometimes Ek(D ) would begin above the dashed line, decrease through it, then turn and converge from below. As in [10], it is clear that nonempty initialization can greatly reduce the time for the expected delays to fall within a specified band about the steady-state value. Thus, in simulating systems such as these, it would be advisable to investigate alternative initializations, especially if a large number of replications are to be made for purposes of statistical analysis, in order to reduce bias and shorten the length of nonproductive warmup periods. Although optimality studies such as those in [10] could be carried out, it is anticipated that the results and recommendations would be similar. In particular, empty and idle (c = 0) initialization may not be a good idea unless p is quite small. 5. CONCLUSIONS In this paper we have derived new discrete-time transient results for two classes of single-server queueing systems that admit a rich family of

18 interarrival or service time distributions; explicit algorithms for computation of the required probability mass functions have also been provided. Given these probabilities, several system performance measures were derived and one of these, the expected delay in queue, was numerically evaluated. The results of this evaluation were used to investigate the choice of initial conditions in simulations of models such as these aimed at estimating steady-state characteristics, The warmup period may be greatly affected by this choice, indicating that some initial experimentation with a given simulation model may prove fruitful in identifying good starting conditions. APPENDIX To prove Propositions 1 - 3, let A denote an exponential interarrival time random variable and let S be an m-Erlang service time; thus, E(A) = 1/X and S = S + + Sm where the S. s are independent exponential with common mean l/(mu). Then P(A < Si) = X/(X + mg) = p/(p + m). Proof of Proposition 1: Pk(n,k+nm) = Pk(X = k+nm), and since X < k + nm, ~ k n n - this is the probability of no service stage completions in [O,t ], a period of n n consecutive interarrivals. If k > 1 then a service stage is in progress at time 0, at which time it renews (due to exponential memorylessness), and will again renew at each arrival. Thus, Pk(n,k+nm) is the probability of n interarrivals during a single service stage, so is P(A < S1' A2 < S,.. An < S1) = [P(Ai < S1)]n = [p/(p + m)]1, where A. is the ith interarrival. On the other hand, if k = 0 then no departure in [0,tl] is possible, so X = k + nm is equivalent to interarrivals A2,..., A occurring during a single service stage, having probability n

19 [p/(p + m)]n 1. Proof of Proposition 2: The event here is that exactly i service stages are present just after the first customer arrival. This is equivalent to there being k - i + m service stage completions in [O,t1] and the arrival must occur before the next service stage completion. Thus, the event is {s < Al S > A} {1 Sk-i+m A Sk-i+m+l> ' since the interarrival time A renews at the time of each service stage completion. Thus, the required probability is [m/(p + m)] [p/(p + m)]. Proof of Proposition 3: The event is that X = i, and we will condition on n X = j, in which case j must be at least i - m if X is to equal i. Further, since m < X < k + (n - l)m, the range of j is - n-1 - max{i-m,m} < j < k + (n - l)m. Thus, k+(n-l)m Pk(n i)= Pk(X = i I Xn = j Pk(n-j). j=i-m To compute the conditional probability, note that if Xn1 = j, then X = i exactly if j - i + m service stage completions occur in the space of the single interarrival tn t1 and the nth arrival occurs before an additional service stage completion. Thus, Pk(X = i Xn- = j) = [m/(p + m)]-i+m[p/(p + m)], the desired result after simplification. To prove Propositions 4 - 6, it is convenient to let A. and S denote an interarrival stage completion and a service time, respectively; an interarrival time is thus A = A + + A. With this notation, P(Ai < S) = l( m mX/(mX + u) = mp/(mp + 1).

20 Proof of Proposition 4: Case 1: k < m. Since k < m, no customers are initially present (physically) in the system, so the initial departure rate is 0. The required probability is that of the event of no departures in [0,tj], a period of j consecutive exponential interarrival stages. Since the first arrival occurs with the (m-k)th interarrival stage completion, j < m - k implies that at time t. no customer has arrived, so X. must be equal to j + k. If j > m - k, then the required probability is that of no departures during j interarrival stages, the final j - (m - k) of which are during a period when the server is busy, which is [P(Ai < S)]J-(m-k) = [mp/(mp + 1)]J+k-m Case 2: k > m. Here, the server is initially busy and, as in Case 1, we want the probability of no departures in the span of j interarrival stages; this is [P(Ai < S)]J = [mp/(mp + 1)]j Proof of Proposition 5: X = k + 1 - fm exactly if there were f departures in [0,tl], occurring if f service completions take place during one interarrival stage but the next event is the interarrival stage completion. Thus, Pk(l,k+l-fm) = [P(Ai > S)] P(A. < S) = [l/(mp + 1) f[mp/(mp + 1)], as desired. Proof of Proposition 6: Conditioning on Xjl, L( j+k-2)/mj Pk(j,j+k-fm) = Z Pk(X=j+k-fm I Xj =j-l+k-gm) Pk(j-l,j-l-k-gm). (Al) kg=0 Case 1: j + k - 1 is divisible by m. Then Xj > m, so departures in [tjl,tj] are possible regardless. If the number of such departures is f - g, then

21 0 < f - g < Lj - 1 + k - gm)/m, implying that g < f as well as g < L(j + k - 2)/m on the range of the summation in (Al). However, since f < (j + k - l)/mj - 1 by assumption, it is easy to see that g < f is the binding upper bound on g, so the range on g in (Al) can be reduced to {0, 1,..., f}. Finally, the probability of exactly f - g departures in [tjltj] is [l/(mp + 1)f-g[mp/(mp + 1)], (A2) being f - g departures followed by an interarrival stage completion. Case 2: j + k - 1 is not divisible by m. Thus, m must be at least 2. Then the minimal Xj_1 is at most m - 1, creating the possibility of a zero exit rate; this is true for the g = Ej + k - 2)/mj term in the sum in (Al). In this case, X. must be j + k - gm, so 1 if f = L(j+k-2)/mJ Pk(Xj = j+k-fm I Xj_ = j-l+k- L(j+k-2)/mm) = L(j O otherwise However, if f = L(j + k - 2)/m, then due to the upper bound on f in the proposition statement, we must have f = L(j + k - 2)/mJ < l(j + k - l)/mj - 1. (A3) To show that (A3) cannot hold, note that since j + k - 1 is not divisible by m, there are integers h and r with 1 < r < m such that j + k - 1 = hm + r. Thus, Lj + k - 2)/mJ = h + (r - 1)/mj = h, since 1 < r < m, and Lj + k - l)/ml - 1 = h + L(r - m)/mj = h + Lr/mJ - 1 = h - 1, since r < m. Thus (A3) is contradicted, so that we must have f <

22 Lj + k -2)/mJ, for f in the range considered by the proposition. Therefore, the g = ij + k - 2)/mJ term in the sum in (Al) drops out, and the range on g can be restricted to 0 < g < Lj + k - 2)/mj - 1. For such g, Xj1 > m, i.e., there is at least one customer physically in the system after t.j. In this case, X. = j + k - fm if and only if f - g departures occur in [tjl,tj]. Also, g < min{f, 1j+k-2)/mj-1}, and this minimum is easily seen to be f. Thus, the range on g in the sum in (Al) is {0,1,...,f}, and Pk(Xj = j+k-fm | Xj_1 = j-l+k-gm) is given by (A2), completing the proof. REFERENCES 1. Carson, J.S. and Law, A.M. Conservation Equations and Variance Reduction in Queueing Simulations. Operations Research 28 (1980), 535-546. 2. Clark, G.M. Use of Polya Distributions in Approximate Solutions to Nonstationary M/M/s Queues. Commun. ACM 24 (1981), 206-217. 3. Gafarian, A.V., Ancker, C.J., Jr., and Morisaku, T. Evaluation of Commonly Used Rules for Detecting 'Steady State' in Computer Simulation. Naval Res. Logist. Quart. 25 (1978), 511-529. 4. Gaver, D.P. and Shedler, G.S. Control Variable Methods in the Simulation of a Multiprogrammed Computer System. Naval Res. Logist. Quart. 18 (1971), 435-450. 5. Grassman, W.K. Transient and Steady State Results for Two Parallel Queues. Omega 8 (1980), 105-112. 6. Halfin, S. and Whitt, W. Heavy-Traffic Limits for Queues with Many Exponential Servers. Operations Res. 29 (1981), 567-588. 7. Heathcote, C.R. and Winer, P. An Approximation for the Moments of Waiting Times. Operations Res. 17 (1969), 175-186. 8. Hillier, F.S. and Yu, O.S. Queueing Tables and Graphs. North Holland, New York, 1981. 9. Kelton, W.D. and Law, A.M. A New Approach for Dealing with the Startup Problem in Discrete Event Simulation. Naval Res. Logist. Quart. 30 (1983), 641-658.

23 10. Kelton, W.D. and Law, A.M. The Transient Behavior of the M/M/s Queue, with Implications for Steady-State Simulation. Operations Res., forthcoming. 11. Kleinrock, L. Queueing Systems, Vol. 1:.Theory. Wiley, New York, 1975. 12. Morisaku, T. Techniques for Data Truncation in Digital Computer Simulation. Ph.D. Dissertation, Department of Industrial and Systems Engineering, University of Southern California, 1976. 13. Morse, P.M. Stochastic Properties of Waiting Lines. J. Operations Res. Soc. Amer. 3 (1955), 255-261. 14. Odoni, A.R. and Roth, E. An Empirical Investigation of the Transient Behavior of Stationary Queueing Systems. Operations Res. 31 (1983), 432 -455. 15. Pegden, C.D. and Rosenshine, M. Some New Results for the M/M/1 Queue. Management Sci. 28 (1982), 821-828. 16. Rothkopf, M.H. and Oren, S.S. A Closure Approximation for the Nonstationary M/M/s Queue. Management Sci. 25 (1979), 522-534. 17. Saaty, T.L. Time-Dependent Solution of the Many-Server Poisson Queue. Operations Res. 8 (1960), 755-772. 18. Schruben, L.W. Detecting Initialization Bias in Simulation Output. Operations Res. 30 (1982), 569-590. 19. van Doom, E. Stochastic Monotonicity and Queueing Applications in BirthDeath Processes. Springer-Verlag, New York, 1981. 20. Whitt, W. Comparing Counting Processes and Queues. Adv. Appl. Prob. 13 (1981), 207-220. 21. Wilson, J.R. and Pritsker, A.A.B. Evaluation of Startup Policies in Simulation Experiments. Simulation 31 (1978), 79-89.