The Transient Response of the M/Ek/2 Queue and Steady-State Simulation Joseph R. Murray W. David Kelton Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, Michigan 48109-2117 Technical Report 85-29

The Transient Response of the M,/Ek/2 Queue and Steady-State Simulation* Joseph R. Murray W. David Kelton Department of Industrial and Operations Engineering The University of Michigan Ann Arbor, Michigan 48109-2117 The probabilistic structure for the transient MjEk/j2 queue is derived in discrete time, where Ek denotes a k-stage Erlang distribution. This queue has a two-dimensional state space. Expressions in terms of transition probabilities are formulated for the expected delay in queue. Results are numerically evaluated for one case. The convergence behavior is similar to that seen in previous work on queues with one-dimensional state spaces. The implications for bivariate initialization of steady-state simulations are discussed. Key Words: Queueing systems, Transient results, Simulation, Initialization of simulation models. * This research was supported by contract fromi Electroili Data Systems Decision Technologies Division t< thet IlUniversity of Michigan.

1. INTRODUCTION A review of the literature in queueing theory reveals an abundance of results for steadystate conditions and relatively few results for the transient phase of a queueing system. One reason for this is the complexity and intractability of the mathematics involved in solving the transient problem, and it is not uncommon to see results stated in terms of transforms which are very difficult, if not impossible, to invert. Results, when left as transforms, seem somewhat less satisfactory in terms of ease of interpretation. Analytical solutions for the transient characteristics of queueing models are useful for studying the finite-time properties of systems that are accurately represented by such models. Even if exact analytical transient results are not known, it would be useful to know, in some fashion, the rate and manner (e.g., monotonic or oscillatory) of convergence of the system to steady-state. Analytical transient results are also valuable in the evaluation of alternative start-up strategies for simulations aimed at estimating steady-state parameters. Kelton and Law [6] and Kelton [5] present transient results for M/M/s, M/Ek/l, and EkI/MI1 queues and use them as benchmarks to evaluate alternative initialization methods for simulation of similar systems. Enlarging the range of benchmark models to include systems with multivariate state spaces and multiple servers with non-exponential service times served as the main motivation for this paper. Another simulation-related area where transient results can be applied is the external control variates technique for variance reduction. When examining the transient behavior of a complex, analytically intractable system in a terminating simulation, a simpler system with known transient behavior is simulated alongside the system of interest. The results from the two systems, assuming the use of common random numbers, would be expected to be correlated, leading to a variance reduction. The larger the class of systems for which transient results are known analytically, the greater the similarity possible, leading to stronger variance reductions. Known transient results can be classified according to whether or not the time measure is continuous (real time) or discrete (indexing by customer number). Continuous time results for various queues can be found in Saaty [15], Kleinrock [9], Odoni and Roth [12], Grassman [3], Pegden and Rosenshine [13], etc. Continuous time results describe the behavior of the parameters of the system, e.g., the number of customers in the system, at every point in time. Discrete time results, on the other hand, focus on the state of the 1

system at certain transition points, e.g., at the point of arrival of the nt1 customer, or at the point of departure of the nth customer. The treatment of queues in discrete time is especially relevant from the standpoint of simulation. Standard measures of performance of general GI/G/s queueing systems include the expected delay in queue (excluding service time), denoted by d, the expected wait in system (delay in queue plus service time), w, the expected number of customers in the queue, Q, and the number of customers in the system, L. Estimates for all of these quantities can be made directly during a simulation, but Carson and Law [2] have shown that it is preferable (in terms of achieving a reduction in variance) to estimate w, Q, and L indirectly from a direct estimate of d using the conservation equations: w = d + E(S), Q = Ad and L = A(d + E(S)) where E(S) is the expected service time and A is the arrival rate. For this reason, most simulation application and methodological research has focused on the delay-in-queue process, which clearly is in discrete time. Hence, analytical results for queueing systems in discrete time can be more easily related to simulation methodology. Morisaku [11], Kelton and Law [6], Kelton [5], Moore [10], Heathcote and Wiener [4], Stanford et.al. [16] and Bhat and Sahin [1] present discrete-time results for various queueing systems. In this paper, we present discrete-time transient results for an M/Ek/2 queue, where Ek denotes the k-Erlang distribution. The state variable for this system has to be expressed as a tuple. The method of analysis used here admits arbitrary (deterministic) initial states for the queue; this allows a numerical evaluation of the effect of alternative initial states on the nature of convergence to steady-state, a general question of interest in simulation aimed at estimating steady-state parameters. In Section 2 the analytical results for the M/Ek/2 system are derived. Section 3 reports some of the results based on the numerical evaluation of these results. Section 4 contains some conclusions, and the Appendix contains proofs of the results presented in Section 2. 2

2. THE M/Ek/2 SYSTEM The arrival rate to the queue is denoted by A, and the service times have a k-Erlang distribution with mean 1/f,. The service time distribution at both servers is assumed to be identical. It is also assumed that a customer professes no preference for one or the other of the servers, but enters service at the first available server. An arriving customer who finds both servers idle chooses one of the servers at random. For purposes of analysis, each complete service period is modeled as k consecutive independent exponential stages, each at rate k1c.. The traffic intensity is p = A/(2p). Let T n = 1, 2,3... be a random variable that represents the time of arrival of the nth customer to the system. Let At be the number of service stages yet to be completed at server 1 at each point of continuous time t, t > 0, and let Bt be the total number of service stages present in the system at time t. (Server 1 is idle at time t if and only if At = 0.) The pair X(t) = (At, Bt) is sufficient to describe the state of the system at time t, since other quantities, such as the number of customers in queue at time t or the number of service stages yet to be completed at server 2 at time t, are functions of X(t). Note that 0 < At < k and At < Bt. The process X(t) renews at each point of continuous time t (i.e., the evolution of (A,, B,) for s > t is a function only of (At, B,), and is independent of all that happened in [0, t) ), because both the interarrival time and the service stage distributions are exponential. X(t) is, in fact, a continuous time Markov chain. It is easily seen that the T 's are stopping times for the process X(t), and X(t) therefore renews at each time Tn. In other words, Xn= X(T,) = (AT,, BTn) is a Markov chain. Similarly, the process X(t) also renews at each (random) epoch in time when a service stage is completed. In the next Section, the transition probabilities for X, are presented. In Section 2.2 various quantities of interest are derived in terms of these transition probabilities. 3

2.1 Transition probabilities The process being analyzed is defined as X,, = system state at T,, including the k stages of the nth arriving customer, for n > 1. Let x(, = (ao, b,) be the system state at time 0. The probabilities of interest are: Px, (x; n) = P(X, = x Xo = xo), n = 1,2,3,... Note that the range of values x can take is determined by x( and n. The first arrival occurs at T1, which is exponentially distributed with mean 1/A. The following Propostions, the proofs for which are given in the Appendix, are sufficient for the computation of the P,, (x; n)'s. For convenience, let aC = A/(A+ —iky}) and /P = 1-ai, for i = 1,2. a l is the probability that, given only one server is busy, the next customer arrives before the next service stage completion. Similarly, ua2 is the probability that, given both servers are busy, the next customer arrives before the next service stage completion. /1 and /2 are the probabilities of the complementary events. (Note: Xo = (a,, bo) and x = (a, b) below.) Proposition 1: P(.() ((k k); 1) = 1 Proposition 2: (a) for 1 < b, < k P(o.,,) ((k, k);1) -=Al f (b) for 1 < bo < k and k+ 1 < b < b(, + k P(o.,, ((k, bl); 1) = Cl ~ - 1l+k This propostition represents the initial condition where server 1 is idle, server 2 is busy, and there is no queue. Proposition 3: (a) for 1 < bo < k P(0.10) ((k, k); 1) = 1b~ (b) for <j < bo < k P(bo.bo) ((jj, + k); 1) = a1lj31-3 This proposition represents the initial condition where server 1 is busy, server 2 is idle and there is no queue. 4

Proposition 4: For 1 < ao < b(, < k+ 1, (a) P(,,)) ((k, k); 1) = 31'<' + f''-"(' - 1 n1l =() tM)-=-) + S 2 S &2fi/2 ( ll 2 /2~) 3 (b) for 1 <j < ao, P(f).,,o) ((, k +j); l) =a,", -j -C 2 (2/2)- (2/2)nl + l ao -) n1 =O (c) for 1 < j < (bo - ao), P(LO,.,()) ((k, k + j); 1) = l/~"~-' -?.1=0 ni (d) for 1 < a < a, and a + k + 1 < bo + k, P(o.0) ((a, b);)) = C2( 32/2)ob —(k (bo- b + k) This Proposition holds for the initial conditions where server 1 and server 2 are both busy and there is no queue waiting for service. The Proposition also holds for the case when bo > k + 1 as long as bo - a( < k. Proposition 5: For 1 < ao < k and bo > 2k, (a) for b> 2k+ 1, P(ao.b) ((a, b); 1) = 22(2/2)botk n (o- b + bk n=0 c / where: c = ao- aifao > a and c = a- a+ kifao < a n' = (bo - b+ k- c)/kJ ([zj)denotes the greatest integer < x.) and, if b < bo, then 1 < a< k, else, if b > bo, then a E { aj,j = 0,1, 2,...,j'}, where j'= bo- b+k 5

aj = o, - j if a,) > j aj = ao - j + k if ao < j. (b) for b < 2k +, P(. l) ) (({a, b); l) = (f2/2)bo k-1 f (b- k- 1 P(^.k+l) ((a, b); ) j = 1 n=O where: P(,L.k+l) ((a, b); 1) is found from Proposition 4. c = a - a, if ao > al and c= ao - a, + k if ao < a,, and, n'= (bo - k- 1- c)/kJ. This proposition represents the initial conditions where both servers are initially busy, and there is at least one customer in the queue. The proposition also holds for the case with bo < 2k as long as bo - ao > k. Propositions 1 through 5 have covered the cases where b, < k + 1 or bo > 2k. There are two subcases when k + 2 < bo < 2k: (a) if b6 - a,, < k, then both servers are busy and the queue is empty initially, hence Proposition 4 applies, and (b) if b0 - a0 > k, then both servers are busy and there is one customer in queue initially, hence Proposition 5 applies. Proposition 6: bo)+(n-l)k k P(lao.10) ((a, b); n) = E P (o.o ((ij); n - I)P(ij, j((a, b); 1) j=k i=l This follows directly from the fact that the process X(t) renews at every T. The quantities P(i.j) ((a, b); 1) can be found using Propositions 1 through 5. 6

2.2 Applications If the probability mass function P() (x,; n) of X, is known, then formaulae for several standard measures of queueing performance can be developed easily. Some examples of these measures are: the expected total number of stages present in the system just after T., the expected number of customers present in the system just after T,, the expected number of customers in queue just after Tn and the expected delay in queue for the nth customer. Only the last performance measure is considered in detail below. Let Dn denote the delay in queue of the nth arriving customer. Then, Ex((Dn) = E E(Dn | X, = Xn)Pxo(xn; n) xn The quantity of interest is therefore E(Dn Xn = Xn). Obviously, if b6 < k + 1, then E(D, I (a,, bn)) = 0. Also, if (b, - an) < k, then E(Dn I (a,, b)) = 0. Suppose x, is such that the nth customer has to wait in queue before being served. The earliest time when this customer can enter service is after bn - 2k service stage completions, and the latest time when this customer can enter service is after bn - (k + 1) service stage completions. Let Q(a, b) be the probability that there will be exactly a service stages remaining at server 1 and (b- a) service stages remaining at server 2 when the nt" customer enters service. It is clear that either a or (b - a) or both will be equal to k. Let EZ, = Expected time for j service stage completions, when both servers remain busy throughout the period required for these j stage completions. The rate of service stage completions when both servers are busy is 2ku/, so EZy = j/(2kus). Conditioning on the total number of remaining service stages at server 1 and server 2 when the nth customer just enters service, it can be seen that 2k-1 E(Dn | (an,, b)) = EZb,,b(Q(k, b) + Q(b - k, b)) + EZb,-2kQ(k, 2k) b=k+l (2k —1 O -b (b- 2k) = ( ^-2k (Q(k, b)+ Q(b-, b)) + k-^ Q(k, 2k) b= i i i i ill The formulae for Q(a, b) are given and derived in the appendix. Finally, 2k-1 Ex (Dn)= Po(xn; n){ i (Q(k, b) + Q(b - k, b)) xn:(bn,-an )>k b=k+l + Q(k, 2k) b - 2k} 7

The cumulative distribution function of Dn can be calculated in a similar manner, using the Q(a, b)'s. We empirically confirmed our results by simulating an M/E2/2 queue, initially empty and idle, with p = 0.7 and A = 1. The delay in queue of the first 25 customers was computed by averaging over 100,000 replications and the difference between the observed and theoretical values was less than 0.7%. 3. IMPLICATIONS FOR INITIALIZATION OF SIMULATIONS In general, the goal of steady-state simulation is to estimate properties of the steadystate distribution. A judicious choice of the initial state can result in a reduction in the time required to reach steady-state. Clearly, the best choice for the initial state would be the steady-state value(s), but lack of knowledge precludes this choice. Kelton [6] and Kelton and Law [7] present results for systems with a one-dimensional state-space (i.e., X, is a scalar), and show that it is better to choose an initial state other than the popular empty-and-idle state to promote convergence to steady-state. The results in this paper extend the work to include two-dimensional state space systems, which require a bivariate initial state. The assumption of Erlang service time distributions lends realism to the model, since it appears that for many processes an Erlang-shaped histogram arises from the service time data. Figure 1 shows a plot of Ex, (D,) as a function of n, for an M/E2/2 system with traffic intensity p = 0.7 and A = 1., and x(, = {(0,0),(2,4), (2,6), (2,8), (2,10), (2,14), (2,20) }. Kelton and Law [7], Kelton [6] and Bhat and Sahin [1] took advantage of the special structure of the transition matrix for systems with a one-dimensional state-space and efficiently computed results for large values of n. We could not find a similar efficient computational algorithm for our results. Obtaining the numerical results, therefore, was computationally very expensive. * The value for the expected steady-state delay in queue for this system (dashed line) was found from tables in Hillier and Yu [5]. The convergence of E,,(Dn) to steady-state is highly dependent on x0 and is nonmonotonic in some cases. Similar behavior was observed using discrete time analysis by Kelton [6], Kelton and Law [7] and Stanford et.al., [16], and using continuous time analysis by Grassman [3] and Odoni and Roth [12]. Odoni and Roth identified four types of * We were limited to maximum values of n sH 30. A program written in Fortran 77 took h 200 minss. of CPU time to evaluate the delays of 30 customers on a Harris 800 computer and the VOS operating system. 8

behaviour: i) monotonic convergence from below, the function being concave in time, ii) initial decrease in the function, followed by a monotonic increase to the steady state value, iii) monotonic convergence from above, the function being convex in time, and, iv) monotonic convergence from above, the function being convex in time, but with a linear decrease initially. No claim was made that these types of behavior were exhaustive, and the basis for the characterization was empirical observations. The four types of behavior were observed for: i) an empty and idle or near-empty initial state, ii) initial state near steady-state value, iii) initial state > steady-state value, and iv) initial state >> steady-state value, respectively, as illustrated in Fig. 1. It is clear, as in Kelton [5] and Kelton and Law [6] that the time for the expected delays to fall within a specified tolerance zone around the steady-state value is greatly influenced by the initial state. It is therefore advisable to investigate alternative initializations for simulation of systems such as these, in order to reduce bias and to shorten the length of the non-productive warm-up periods. A method similar to the one discussed in Kelton [5] that uses a series of preliminary runs, can be employed to find reasonable values for initialization of actual production runs. 4. CONCLUSIONS Results for the transition probabilities for the transient, discrete time M/Ek/2 system have been presented in this paper. Further results using these transition probabilities have also been derived. It is shown that the results are in close agreement with previously published results, which were solely for single dimensional state-space systems. A significant conclusion is that the method of analysis used here, though efficient for a scalar state-space appears computationally inefficient for higher dimensional state-spaces. Simulation of models like the one investigated is greatly influenced by the choice of the initial state, thus warranting some experimentation to identify good starting conditions. Good starting conditions would result in a quicker approach to steady state and consequent reduction in computer run-time for the simulations. The same method of analysis can be used for the M/Ek/s system for s > 2, but the benefit of such an analysis is questionable, unless some method for reducing the computation time is found. Other multi-dimensional state-space models, e.g., M/M/1/M/1 queue can 9

possibly be studied in a similar manner, thereby providing more analytical test models for multivariate initialization heuristics for simulations. APPENDIX In the following proofs, when the initial state Xo is implicitly known, the subscript in Px0 (x; n) is dropped for notational convenience. Proposition 1 is trivially true. Proof of Proposition 2: (a) Let Ai denote the exponential arrival time of the first customer. Let Sn =1, 12,... be a random variable that denotes the time of the nth service stage completion. P((k, k); 1) is the probability that there will be bo service stage completions before the first arrival. Since be( > 1, there is a service stage in progress at time 0. Because this process is memoryless, it renews at time 0. The probability that the completion of the service stage in progress is before the t first arrival is kp/(A + ku), since the two competing processes have exponential distributions. At the instant of completion of the first service stage, i.e., at S1, the arrival process renews, and therefore P(S2 < A1) = ks/(A + kyu). Continuing in this fashion, P(S, < Al,S2 < A,...,S,(< A)= A ( = A - (b) P((k, bl); 1) is the probability that there will be exactly bo - (b1 - k) service stage completions before the first arrival. At the instant of the (b0 - (b61 - k))th service stage completion, the arrival process renews, and the probability that the arrival will be before the next service stage completion is A/(A + ky). As in (a), the probability that there will be bo - (bl - k) service stage completions before the first arrival is (ki /(A + ki)) ~,'-k The result follows immediately. Proof of Propositon 3: (a) Same as the proof of Proposition 2(a). (b) Note that if all of the b6 service stages in the system at time 0 are completed before the arrival of the first customer, then (a) applies. If, on the other hand, the number of service stage completions before the first arrival is i (O < i < bo - 1), then X1 = (bo - i, b0 - i+ k). Define j = bo - i. Since 0 < i < bo - 1, we have 1 < j < b0 and the result follows using Proposition 2(b). 10

Proof of Proposition 4: (a) P((k,k); ) is the probability that the first arriving customer finds the system empty. Define random variables Ci, C2 as follows: C, = no. of service stage completions from server i in [0, Ti]. Note that as long as server i is busy, Ci is a Poisson counting process with rate kmI. Further, C1 and C2 are independent. Initially, there are a, service stages remaining at server 1, and since bo < k + 1, there are bo - a(, service stages remaining at server 2. If the first arriving customer finds the system empty and idle, then C1 = ao and C2 = b6 - a(. Conditioning on T1, the arrival time of the first customer, we get: P(C = ao, C2 = b( - a,)) = ET1 (P(C = ao, C2 = bo - a I T1)) 00 = P(C = a,, Co = bo, - a(, T, = x) A e- A dx Since C1 and C2 are independent processes, P(C1 = ao, C2 = bo - ao T = x) = P(C, = a,, I T - x) P(C2 = b - aD T, = x). Since C1 cannot exceed ao, P(C T x) 1P(C = a(C < ao T1 = x),%-1 =1- E exp(-kx)(kll"x)/"n! nl=0 Similarly, P(C2 = bo - ao I T = x)= 1 - P(C2 < bo - ao I T1 = x) b0 - ao - 1 = 1 - E exp(- kzu)(kpXz)"2/n2! n2 =0 Finally, P(C1 = ao,C2 = bo- ao) = P 0-1 o, - ao - 1 /[ i - E k (k.Xz) /nlr! [1- e- b(klz) ]"/n A2!e- dx dz oo n= in an i iiin2= The result follows after integration and simplification. 11

(b) P((j,j + k); 1) is the probability that the first arriving customer finds j service stages remaining at server 1, and server 2 idle. This translates to C1 = a( - j and C2 = b6 - at. The rest of the proof follows the same lines as that of (a). (c) P((k, k+ j); 1) is the probability that the first arriving customer finds server 1 idle and j service stages remaining at server 2. This translates to Cl = a( and C2 = b(, - a0 - j. The rest of the proof follows the same lines as that of (a). (d) P((a, b); 1)) is the probability that the first arriving customer finds both servers busy. This translates to C1 = ao - a and C2 = (bo - a) - (b - a - k). The rest of the proof follows the same lines as that of (a). Proof of Proposition 5: The initial state for this proposition implies that at time 0 both the servers are busy and there is at least one customer in the queue. (a) The case b > 2k + 1 is examined here. Let C1 and C2 be as defined in the proof of Proposition 4(a). With c and n' as defined in the Proposition statement, it is easily seen that X1 = (a, b) if and only if C1 + C2 = bo - (b - k) and C1 E {c, c+ k, c+ 2k,..., c+ n'k}. n' can be physically interpreted as the maximum possible number of customer exits from server 1 during the first interarrival time. The rest of the proof follows the same lines as that of Propostion 4(a). It should be noted that if 2k + 2 < b < b0, any value of a E {1, 2,..., k} leads to a feasible X1 = (a, b). If, however, b > b6, then let j' = bo + k - b. Clearly, 0 < j' < k- 1. The feasible values for a are {a: j = 0,1,2,...,j' - 1} where a3 = ao -j if a, > j and aj = a - j + k if ao < j. (b) If b < 2k + 1, then the process has to have passed through the state (a,, k + 1) at some t E [0, T1), and a, E {1,2,...,k}. Designate this intermediate state as XI; XI is entered at some random time TI < T1. TI is a stopping time for X(t). Hence the process renews at time TI. To reach the state XI requires that (bo - (k + 1)) service stage completions occur before the first arrival. Defining c as in the proposition statement it can be seen that in order to reach the state (a,, k + 1) two conditions need to be satisfied: (a) bo - (k + 1) service stage completions occur before the first arrival and (b) the number of service stage completions from server 1 be c plus a multiple of k. Let Ctot = bo - (k + 1) and o = L(ctot - c)/kJ. n' can again be physically interpreted as the maximum possible number of customer departures from server 1 during the first interarrival time. Then the number of service stage completions at server 1 must be in the set {c, k + c,..., n'k + c}. Given that 12

there were c0t, service stage completions before the first arrival, the probability that any of these service stage completions is from server 1 is 1/2. Therefore, P(cl service stage completions at server 1 ctt)= (c tt)(1/2) tot which is the symmetric binomial probability mass function. Using this, we obtain P{bo - (k + 1) completions before first arrival, X = (a, k + 1)} A + 2 k tl) j ) (1,nk2+ c) ( k )totI Ctot n=0 Now, P(X1 = (a, b) | Xo = xo) = C P{bo - (k+1) comp. before first arr., XI = (al, k + 1)}. aI P{X1 = (a, b) I Xo = xo,X, = (a,,k+ 1)} Since TI < T1 and because X(t) renews at TI, P{X1 = (a, b) I Xo = xo,X = (a, k + 1)} = P{X1 = (a, b) IXo = (a,,k+ 1)}. The expression on the right hand side of the above equation can be evaluated using Proposition 4, and the result follows immediately. Derivation of Q(a, b)'s Q(a, b) is defined as the probability that, given the nth customer does not go into service immediately on arrival, there will be a and b - a service stages remaining at server 1 and server 2 respectively when the nth customer finally enters service. Though not explicit in the definition, Q(a, b) clearly depends on Xn. All arguments below are conditioned on the event that the nth customer has to wait in queue. Also, service stage completions, unless otherwise qualified, mean service stage completions since the arrival of the nth customer. Case 1: bn > 2k. At the instant when a total of b, - 2k service stage completions have occured, the nth customer is either the first in queue or has just entered service. Let Q(i, 2k) be the probability 13

that there will be i remaining sevice stages at server 1 when a total of b, - 2k service stage completions have occured. Both servers are busy during these b, - 2k service stage completions, and the rate of service stage completions is 2kyi. The probability that any of these completions is from a particular server is 1/2. If an > i, then let c = a, - i; else, if a, < i, then let c = a - i+ k. Let n' = L(bn - 2k- c)J. Then, using the same arguments as in the previous proofs, Q (i, 2k) =E ) (1/2) bn-2k n=O Clearly, Q(k,2k) = Q(k,2k). Let Q'(i,j) be the probability that after bn - j service stage completions there are i remaining service stages at server 1, conditioned on the event that the nth customer did not enter service after bn - j - 1 service stage completions. Because both servers are busy during the (b - j)th service stage completion, the following equations hold: For j = 2k - 1, '(1, 2k - 1) -(1/2)Q'(1,2k) + (1/2)Q'(2,2k) q(k- 2,2k- 1) =(1/2)Q'(k- 2,2k) + (1/2)Q'(k- 1,2k) Q'(k- 1,2k- 1) =(1/2)Q'(k- 1,2k) (k, 2k - 1) =(1/2)Q'(1, 2k) Similarly, for j = 2k - 2, Q'(1,2k - 2) =(1/2)Q'(1,2k - 1) + (1/2)Q'(2,2k - 2) q9(k - 3,2k - 2) =(1/2)Q'(k - 3,2k - 1) + (1/2)Q'(k - 2,2k - 1) Q(k - 2,2k- 2) =(1/2)Q'(k - 2,2k- 1) '(k, 2k - 2) =(1/2)Q'(1,2k - 1) And in general, for j such that k+ 1 < j <2k- 1, Q'(1,j) =(1/2)Q'(1, j + 1) + (1/2)Q'(2,j + 1) 14

Q(j - k - 1,j) =(1/2)Q'(j - k - 1,j + 1) + (1/2)Q'(j - k,j + 1) Q'(j - k,j) =(1/2)Q'(j - k,j + 1) Q(k,j) =(1/2)Q'(1,j + 1). It should be clear that for k + 1 < j < 2k - 1, Q(k,j) = Q(k,j) and Q(j - k,j) = Q (j - k,j). Since the Q2(i, 2k)'s are known, the Q(i,j)'s can be calculated. Case 2: b, < 2k. In this case, the nth customer is the first in queue upon arrival. Therefore, Q(i,j) = 0 for b, < j < 2k. Using the Q<'s as defined in Case 1, set Q(a,, bn) = 1, and Q(i, bn) = 0 for if an, i E {1, 2, **,bn - k, k}. Then, (i,j), for k+ 1 < j < b - 1, can be calculated recursively using the formulae in Case 1. Finally we have as in Case 1, for k+1 < j < b, - 1, Q(k,j) = Q(k,j), and Q(j- k,j) = Q(j - k,j). ACKNOWLEDGMENTS We would like to thank Electronic Data Systems Decision Technology Division for their support, and in particular express our appreciation to Ziv Barlach, Sam MacMillan and Michael Moore. 15

REFERENCES 1. Bhat, U.N., and Sahin, I. Transient behavior of queueing systems. M/D/1, M/Ek/1, D/M/1 and Ek/M/1. Tech. Memo. 135, Dept. of Operations Research, Case Western Reserve Univ., 1969. 2. Carson, J.S., and Law, A.M. Conservation equations and variance reduction in queueing simulations. Oper. Res. 28, 1980, 535 - 546. 3. Grassman, W.K. Transient and steady-state results for two parallel queues. Omega 8, 1980, 105 - 112. 4. Heathcote, C.R., and Winer, P. An approximation for the moments of waiting times. Oper. Res. 17, 1969, 175 - 186. 5. Hillier, F.S., and Yu, O.S. Queueing Tables and Graphs. North Holland. New York, 1981. 6. Kelton, W.D. Transient Exponential-Erlang Queues and Steady State Simulation. CACM 28, 1985, 741 - 749. 7. Kelton, W.D., and Law, A.M. The transient behavior of the M/M/s queue, with implications for steady-state simulation. Oper. Res. 33, 1985, 378 - 395. 8. Kelton, W.D., and Law, A.M. A new approach for dealing with the startup problem in discrete event simulation. Naval Res. Logist. Q. 30, 1983, 641 - 658. 9. Kleinrock, L. Queueing Systems, Vol. 1: Theory, Wiley, New York, 1975. 10. Moore, S.C. Approximating the behavior of nonstationary single-server queues. Oper. Res. 23, 1975, 1011 - 1032. 11. Morisaku, T. Techniques for data truncation in digital computer simulation. Ph.D. dissertation, Dept. of Industrial and Systems Engineering, Univ. of Southern Calif., Los Angeles, 1976. 12. Odoni, A.R., and Roth, E. An empirical investigation of the transient behavior of stationary queueing systems. Oper. Res. 31, 1983, 432 - 455. 13. Pegden, C.D., and Rosenshine, M. Some new results for the M/M/1 queue. Management Science 28, 1982, 821 - 828. 14. Rothkopf, M.H., and Oren, S.S. A closure approximation for the non-stationary M/M/s queue. Management Science, 25, 1979, 522 - 534. 16

15. Saaty, T.L. Elements of Queueing Theory with Applications, McGraw-Hill, New York, 1961. 16. Stanford, D.A., Pagurek, B., and Woodside, C.M. Optimal Prediction of Times and Queue Lengths in the GI,/M/1 Queue. Oper. Res. 31, 1983, 322 - 337. 17

(0 o0 LD 0 8 cM 5 10 15 20 n Figure 1. E,,((Dn) for the M/E2/2 Queue with p = 0.7