APPROXIMATION ANALYSIS FOR OPEN TANDEM QUEUES WITH BLOCKING: EXPONENTIAL AND GENERAL SERVICE DISTRIBUTION Stephen M. Pollock John R. Birge Jeffrey M. Alden Department of Industrial & Operations Engineering University of Michigan Ann Arbor, MI 48109 October, 1985 Technical Report 85-30

Approximation Analysis for Open Tandem Queues with Blocking: Exponential and General Service Distributions Stephen M. Pollock, John R. Birge, and Jeffrey M. Alden University of Michigan, Ann Arbor, Michigan Calculating occupancy probabilities for open finite tandem queueing systems using Markov Process analysis is impractical for all but small systems due to the large number of states. We offer an iterative procedure for approximating the marginal occupancy probabilities for each queue of the system. The procedure is easy to implement, requires little memory, and is computationally fast. Two implementations are presented: one for exponential service distributions and one for general service distributions. Both implementations give better accuracy than other approximation methods. Introduction Systems of servers in tandem (series), where the output of one server is the input to the next one in line, can represent many processes of interest in manufacturing, computer systems, telecommunications, etc. This is particularly true when the service times are random variables. When the storage ("buffer") space between servers is finite, so that an upstream server can be "blocked" due to unavailability of space in the next buffer, the resulting system becomes notoriously hard to analyze. (For an excellent review, see Perros [1984].) In particular, exact solutions (i.e., ones that provide the usual measures of queueing performance) are available only for either infinite intermediate buffers (Gordon and Newell [1967], Hordijk and Van Dijk [1981])-in which case solutions are obtained by using product-form representations of state probabilities-or small numbers of servers and buffer sizes (Asare [1978], Caseau and Pujolle [1979], Foster and Perros [1980], Gershwin [1983], Gordon and Newell [1967], Hillier and Boling [1967], Konheim and Reiser [1976], and Labetoulle and Pujolle [1977]). In these studies, service times at each server are assumed to be negative exponentially distributed random variables or, in the case of Gershwin [1983], a probability mixture of a discrete geometric random variable and a constant. 1

One reason for the lack of success in studying this important class of problems is that a comprehensive representation of the state-space of the system requires at least an enumeration (if not a computation!) of state variables that grows combinatorially with the number of servers. This unfortunate fact has led to a number of analyses based upon approximation methods. Most of these (Altiok [1982], Boxma and Konheim [1981], Hillier and Boling [1967], Suri and Diehl [1983], Takahashi et al. [1980]) collapse or aggregate system states into "super-states" or "marginal states" that in turn serve to characterize (to some degree of accuracy) the desired system performance measures. We present here an approximation method that is similar to those presented in Perros and Altiok (forthcoming) and Takahashi et al. (1980), but with simplifications (each server only considers the behavior of-or information from-its immediate neighbors), generalizability (it can be extended to non-exponential service times), and increased accuracy. 1.0 The M-server Tandem System: Exponential Service The fundamental model we will study is identical to the model in Takahashi et al. [1980]. There are M servers in tandem (see Figure 1). Each server has a buffer (waiting room) for holding arrivals when it is busy. Each server has an exponentially distributed service time. Units are assumed to arrive only at the first queue as a Poisson process with a rate A. |- N, __ 1- N2 -e NM e- B1 e B2 H BM - H lit queue 2nd queue H HMth queues Figure 1. Tandem queueing system. Let pi and Bi be the service rate and buffer capacity for the ith server, where i = 1,2, *, M. The ith queue is the combination of the ith server and its buffer. Thus, if Ni is the capacity of the ish queue (assumed finite for all but possibly the first queue), then 2

Ni = Bi + 1. Note that the queues are numbered from 1 to M starting from the input queue. Arrivals to any queue are served in a FIFO manner. An arrival enters the first queue if it arrives at a time when the queue is not full; otherwise, the arrival is lost. When a unit completes service at the ith queue, it proceeds to the (i + 1)et queue if space is available. However, if the (i + 1).t queue is full at that time, the unit must wait in the ith server until a departure occurs from the (i + 1)t queue. During this time the ith server cannot serve other units that might be waiting in its buffer. In this case, the ith server is blocked, and the (i + 1)"t queue is blocking. Thus, the first queue cannot be blocking, and the Mth server cannot be blocked. A complete (steady-state) analysis of this system would produce all state occupancy probabilities, as well as related measures such as the expected number of units in each queue (or buffer), the waiting time distribution at each server, and the throughput of the system (since some arrivals can be turned away, the rate of units through the system is less than A). These measures are functions of the 2M+ 1 parameters A, uAL and Bi (i = 1,2, * *, M). The system state is the vector S = (S1i,5 2,... SM), where si = 0, 1, 2,, Ni, is the number of units in the ith queue and si = Ni + 1 indicates that the ith queue is blocking the (i - 1)" server. The total number of possible states is then TI 1 (Ni +2). For even a reasonably sized system, e.g., M = 8, N, = 4 (i = 1, 2,. -, 8), this number is over 1.7 x 106, and represents a considerable challenge for the computations associated with a straightforward Markov Process analysis. 1.1 The Approximation Method: General Approach Our general approach (labelled SIMP for "simple iterative myopic procedure") is to focus on the measures that seem to be important (for example, individual queue state descriptions) and to seek an approximation that will produce these measures fairly accurately. We sacrifice the possibility of having good information about less important events (such as a specific subset of servers being simultaneously idle). 3

To do this, we analyze each queue separately, using information from only its nearest neighbors in order to approximate, by a simple model, what is in fact a very complicated and dependent state of affairs. In particular, we assume that the model for the ith queue is a simple M/M/1/N queue with the following conditions: al) there are Poisson arrivals, with rate A*, as long as the ith queue is not blocking (that is, as long as the number of units in the ith queue is less than or equal to Ni). When the ith queue is blocking (si = Ni + 1), there are no arrivals. a2) The (random variable) time to complete service has two components: the actual service time (exponentially distributed with rate pi) plus a term due to the occasional and probabilistic delay caused by blocking downstream. a3) the service completion time (having two components, see a2) is exponentially distributed, with rate p. These approximations are, of course, a distortion of what is actually happening in the system. Indeed, al) is by itself a heroic assumption-the input to each queue but the first, being the (buffered) output of the preceding queue, is anything but Poisson; and a3) actually contradicts a2). Nonetheless, these assumptions allow the computation of approximate values for the ith queue (steady-state) occupancy probabilities. We can then use these approximate values to evaluate measures of interest. 1.2 Definitions and Underlying Relationships Let: Si - {0,1,2, * *, NNi + 1} be the state space for the ith queue, si state of ith queue (si E S. is the number of units in the ith queue, except when si = Ni + 1 in which case the ith queue is blocking the (i - l)t server), A* = arrival rate to queue i, Ii* — service rate of server i, f, Pr{ih queue is full} = Pr{si = Ni}+ Pr{si = N, + 1}, b, - Pr{sth queue is blocking} = Pr{s, = Ni + 1}, 4

Ei - expected time between acceptances into queue i, Ti _ (random variable) time for server i to complete service and pass unit on to next queue, and E[Ti] — the expected value of Ti. These definitions, and the structure of the system, produce the following relationships: 1. The rate at which arrivals join the system is A times Pr{queue 1 is not full} = A(1 - f). Thus, El= (1 E1 = A(1- ) (1) 2. The expected time to complete service at server i is the unblocked expected service time 1/pgi plus the expected delay due to blocking: Pr{a served unit at server i will see queue (i+ 1) full} times the expected time for server (i+l1) to complete service. The terms required for the computation are obtained from two additional (approximation) assumptions: a4) A unit at server i at the instance service is finished sees queue (i+ 1) in "steady-state". a5) The expected sojourn time in the blocking state for a queue is equal to its expected service completion time. Given these assumptions, the expected delay due to blocking is equal to server i's unconditional expected time to complete service. Thus, E[Ti]= + fl- b E[Ti+ ], i= 1,2,...,M, (2) lis 1 - hi+ where E[TM+] = 0. Since the viewpoint is that of server i at completion of unblocked service, it is necessary to condition the (i + 1).t queue probability on not being blocking, which produces the denominator of the second term. For the following discussion, we define the steady-state probability that queue (i + 1) blocks server i at the end of service as +1 - bi+l 5

3. The expected time between acceptances at queue i is 1/A, plus E[Ti] times Pr{an arrival sees queue i full}, since in this latter case the arrival must wait for server i to complete its current service before being accepted into the queue. We again use assumptions a4) and a5) to give 1 fi - bi Ei = A + EiE[Ti] i=2,3,..,M. (3) A~ 1 - bi As before, the probability that queue i is full is conditioned on queue i not blocking queue (i - 1), since otherwise there could not be an arrival. Figure 2 shows a schematic representation of this process. Delay in queue i [Queue (i + 1) blocking] Service begins Queue (i + 1) full Service completion in queue i ({i} in queue (i + 1) Service time in queuei J [Queue (i + 1) not blocking] Done Queue (i + 1) not full {1-a ci Done Time Figure 2. Service completion time in queue i. 4. By conservation of unit flow, Ei = E, i = 1,2,.,M. (4) 5. The service completion rate is the reciprocal of the expected service completion time. Thus, 1= E= 2, (5) =~ E[T]' i =1,2,...,M, (5) 6

1.3 Finding the Occupancy Distributions: The M/M/1/N Model Equations 1 through 4 together with standard M/M/1/N analysis suggest the following iterative procedure to find the arrival rate (A*), expected completion time (E[Ti]), and occupancy probability distribution (Pr{si = n}) for each queue. 0. (Setup) Set the values for A, pi, and Ni for i = 1,2, * *, M. 1. (Initialization) Set E[Ti] = 1/pi, bi = 0, A* = A, and A* = A(1- f) for i = 2,3,.., M where fl is given by Equation 6a. 2. (Find full and blocking probabilities) With pi = A E[Ti], set ( (1- pl)+pN f1 I _ - )P, (6a) 1- Nj+l (1- pi)p,+ fi 1-p +2 i=2,3, -,M, and (6b) bi = pif, i = 2,3,-..,M. (6c) 3. (Calculate completion times) Use Equation 2 to solve for E(Ti) in the order i - M,M - 1, —.,1. 4. (Update arrival rates) With Ei = E1 = 1/A(1- fi) (from Equations 1 and 4), use Equation 3 to solve for A* in the order i = 1,2, *, M. 5. (Convergence check) If updated values of As show little change (i.e., convergence), then go to Step 6, else go to Step 2. 6. (Calculate occupancy probabilities) Set Pr{s1 = n} ( n = 0,1,,N1, and I-N+ n=,,..,, andl 1 P Pr{si = n} = ( P i = 2,3,...,M and n = 0, 1,.,Ni + 1, 7. Stop. Queue 1 is given a slightly different treatment in Steps 2 and 6 because the queue never blocks a previous queue. To achieve convergence in the iterative procedure, it was sometimes necessary to introduce a damping factor to reduce excessive oscillation in the values of consecutive iterates. 7

1.5 Computational Results: M/M/1/N Model The approximation method was tested on three sets of problems in the literature with exponential servers. In these problems, analytic solutions are often not available. In such cases, the literature has used the thresults of large simulations to replace analytical values. We use these results in our comparisions as previous authors have used them. Table 1 shows the throughput for a 3-queue system each with a buffer of size one as reported in Takahashi et al. [1980]. As can be seen, SIMP (using Equation 2) performs better (relative to the simulation results) than Takahashi's method, and almost as well as Hillier and Boling's (Hillier and Boling [1967]) method which is more complicated than SIMP, which requires only O(M) elementary operations per iteration. Table 2 shows a comparison with the method of Altiok and Perros (forthcomning) for 3-queue systems with different buffer sizes. Again, SIMP (using Equation 2) performs favorably as compared to the more complicated method. Finally, Table 3 shows similarly good results for highly unbalanced systems of 3 and 5 queues. 2.0 Residual Service Time Considerations The above approximation, although it yields suggestively good results, can be improved by relaxing one assumption. Consider relationship (2) and Assumption a5). This represents a correction to the expected service completion time in the event that queue i "sees" server (i + 1) blocking it. Assumption a4) allows us to compute the probability of this event as aci. Strictly speaking, the (random variable) service completion time for server i is the sum of the unblocked service time (with pdf uie-"^it ) and the (random variable) length of time that server (i + 1) will take to become unblocking, given server i completes service and finds server (i+ 1) blocking. This second component is the (random variable) residual time for server (i + 1) to become unblocking. We now assume a5') the unblocked service of the unit in queue i ends at a random point in server (i + l)'s service completion period. 8

This assumption is less restrictive than a5) and it will allow the consideration of general service time distributions (Section 3). The following relevant result is well known (see, for example Section 5.2 of Cooper [1981]): Let T be a recurrence time with cdf FT(t) and expectation E[T], and let X be the (random variable) time from a random point in the recurrence interval to the recurrence event, then X (called the residual recurrence time) has the pdf x ( - FT(x) (7a) 9x(2)=~}* (7a) E[T] As a consequence, it can be shown that the moments of X are given by E[Xn] E[Tn+l] (7b) (n + 1)E[T]' (7b) and in particular, E[X] - E[T2] (8) 2E[T]' To use this result, we note that, if Ti - (random variable) time for server i to complete service, and Xi - (random variable) residual time for server i to complete service, then Ti is the sum of the unblocked service time (with pdf uieIzit) plus the residual completion time Xi+l, conditional on server (i + 1) blocking. Thus, given a5'), E[Ti] =-+ E[X+ ], (9) Pi which, by Equation 8, gives 1 LE[T1+ ] E[Ti] = - + i, i=M,M-1,..,1. (10) Pji 2E[Ti+i]' Under a5'), Equation 10 should be used, in place of Equation 2, to augment the (unblocked) average service time by a term that accounts for blocking. Similarly, E[Ti] in Equation 3 should be replaced by E[T}2]/2E[Ti]. Note that Equation 10 requires that E[T, J1] be available. Indeed, the recursive use of Equation 10, starting with i = M, requires 9

the terms E[Tin] for n = 1,2, * * *, i. The Appendix shows, by a simple transform argument, that these moments are given by: E[Ti] =n + A, (n) (n- i)E[T,+'1+] (11) E[T] -- "'.j +, j +i-jE[T,+ I]' or written as a recursion in the index n, n E[Trn1] cii E[T\1+ E[Ti] }- E -] + (nn 1)E[T+ n1 2,=.. (12) Thus, accounting for randomly interrupted residual service time requires replacing Equation 2 of Step 3 in the algorithm with Equation 12 (with E[T + 1] = 0 for all j > 1) and replacing E[Ti] with E[T,2]/2E[Ti] in Equation 3 (used in Step 4). By using Equation 12, the complexity of SIMP now becomes O(M2) per iteration. The results, using this modification are also given in Tables 1,2 and 3 under the heading "SIMP Eq. (12)". This modification slightly improves the accuracy of the model. Since, in these cases, we have exponential servers that are frequently unblocked the expected residual service time is appoximately equal to the expected service time and hence we do not expect much change due to this modification. However, the results derived in this section are used in Section 3. 3.0 General Service Time Distributions The above discussion allows us to adapt the approximation method for general service time distributions and eliminate Assumption a3). We are aware of only one analysis in the literature not requiring exponential servers (Gershwin [1983]). Gershwin's approach, however, is restricted to a particular form of non-exponential service. Our method replaces the M/M/1/N analysis in Step 2 (finding bi and fi) and Step 6 (calculating occupancy distributions) with an M/G/1/N analysis. This generalization provides wider applicability at the expense of more computation and greater data requirements-as will be seen, it requires the Laplace transform of the service time distribution (and its derivatives) for each queue. 10

Traditional M/G/1/N analysis (see Section 5.9 of Cooper [1981]) examines the queue at departure epochs. Following this approach, define for each queue i the following (the subscript i is suppressed here for clarity): S (random variable) unblocked service time, g(t) _ service completion time probability distribution function, A arrival rate, _= service rate (inverse of the first moment of g(.)), Pn = steady-state probability that a departure leaves n units behind, an _ probability of n arrivals during a service time, and IIn steady-state probability of n units in the queue. Since the number of units left by a departure is given by the number of arrivals during service and the number left by the previous departure, n Pn+ = aol(Pn-aPo -E an-j+lPj) n =0,1, * **, N - 1, (13) j=l where ~ 0 0. By defining 8n = Pn/Po, S6, can be calculated recursively from (13) after dividing by Po. Using the 8n, Po can obtained from Po N -- (14) zn=O En, and then Pn = SPo. (15) Pn is also the steady-state probability that an arrival sees n units in the queue (see Section 5.3 of Cooper [1981]), and is simply the steady-state probability of n units in the queue conditioned on the queue's not being full. This, combined with conservation of flow, gives Pn ILn = pn = 0- 1= O l, -,N-1, and (16a) Po+p' IIN = - (16b) Po+ where p = A//i. 11

The above equations provide a quick and easy method for obtaining the steady-state probabilities In, given the values of an. To find the values of an, let fi(t), hi(t), and gi(t) be the probability density functions of Ti, Si, and Xi and fi*() (s), h (n) (), and g9( (a) be the nth derivatives of their Laplace transforms. By definition, for the ith queue an = | (,,!) e-((t)dt, or equivalently, an= (! f()()- (17) In the Appendix, we show that the nth derivative of the Laplace transform of the service completion time can be recursively obtained from fi" (s) = (1 - )h(n) (s) + a, () hi (f -i)( *( (s). (18) j=1( The functions h*() (s) are given as data, but the g*i) (s) must be evaluated. Taking the Laplace transform and then the jth derivative of both sides of Equation 7a (with X = Xi+1) gives g W(s) 1 ((1)! - f ( ( j-"(j - ic)! i +x (s) -si E[T+ Sj-k+l f.+ 1 (S) which is equivalent to the following recursion in the index j, *(j+1) 1 *( ( +) + (s) I (_. Zi+ (s)+ E[Ti+ ] where.(o) 1__ _-__()_ g, () - E [T,+ () (19b) Since, in the last queue, fji)(s) = h*()(s) for all j > 0, we can recursively solve for g(j+1) (Ai) and f7*(n) ) using Equations 19 and 18, in the order I = M- 1, M- 2,, i. We then use the results in Equation 17. We also need the first moment E[Ti] to determine the utilization pi = AiE[Ti] used in Equation 16. The generalization of Equation 11 for general service distributions is (see the Appendix) E[Tl ] = E[SnL] + at ( j) +S E[Ti] (20) 12 12

By setting E[T +l]- = 0 for all n > 0 and noting that E[Sr] = (-l)nh(") (0), Equation 20 recursively gives, in the order I = M, M - 1,...,i, all the moments necessary to calculate Pi 3.1 Finding the Occupancy Distributions: The M/G/1/M Model The iterative procedure in the case of general service distributions follows. 0. (Setup) Set the values for A, Ni, and provide a routine to calculate hi*() (s) (i 1,2, *.,M). 1. (Initialization) Set b1 = 0, A; = A, AT = A(l - fi) (i = 2,3,-..,M) where fi is given by Equation 6a, fi = N1 - 1, and ni = Ni (i = 2,3,..,M). (ii is the maximum number of arrivals during a service at queue i.) 2. (Calculate moments) With E[T ] = (-l)hn hn) (0) (n = 1, 2, * * *, M) and ai = (fi+1 - bi+i)/(l- bi+1) use Equation 20 to calculate E[Ti] (n = 1,2, *.,i) recursively in the order i= M- 1,M - 2, -, 1. 3. (Find full and blocking probabilities) For each queue i = 1,2,- *, M let pi = A* E[Ti] and perform Steps 3a to 3d. 3a. Set f(n)(A) = h(n) (A;) for n = 0, 1, *, ni. If i = M then go to Step 3d. 3b. In the order I = M- 1, M- 2,, i use Equations 19 and 18 to calculate g 1) (A*) and f (A)(A*) for n = 0,1,., i. 3c. Calculate a, (n = 0,1, *, ni,) from Equation 17. 3d. Set N = ni + 1 and use Equations 13 to 16 to solve for nN and IN- 1- If i = 1 then set fi = IIN, else set bi = IIN and fi = 1IN + 1N -1. 4. (Update arrival rates) With Ei = El = 1/A(1 - fi) (from Equations 1 and 4), use Equation 3 with E[Ti] replaced by E[T?2]/2E[Ti] to solve for A, in the order i = 1, 2,., M. 5. (Convergence check) If the updated values of A* show little change (i.e., convergence), then go to Step 6, else go to Step 2. 6. (Calculate occupancy probabilities) For each queue i = 1,2, * *, M, repeat Steps 3a to 13

3d (giving IIn) and set Pr{si =n}=Il,, n=0, 1, —,ni +1. 7. Stop. 3.2 Computational Results: The M/G/1/N Model The SIMP M/M/1/N and M/G/1/N Model results are compared with the only available non-exponential service results in the literature: Gershwin [1983]. The compamiibn is summarized in Table 4. The service time at queue i in Gershwin's model is one period with probability 1 - pi and n periods with probability pi ri, n = 1,2,... The server of the first queue is never idle. We achieved this requirement by giving the first queue a buffer of 4 and making the arrival rate large compared to the average service time. In Table 4, B is the vector of buffer sizes, and p and r are the service time parameters (pi,r1) for each queue. Again, simulation values replace exact values when analytical results are not available. We note that the Gershwin model requires, amoung other things, the solution of a nonlinear equation in one unknown in each iteration whereas we have only 0(M3 + ZM=1 N2) computations per iteration (from equation 20 and M/G/1/N analysis). Note also that the form of Gershwin's service distribution has a high variance, which leads to occasional, but very long blocking delays, a phenomenon the M/M/1/N model would miss. The lower throughput values for the M/G/1/N versus the M/M/1/N analysis indicate an improvement in this area. The results in Table 4 are encouraging given the generality of the SIMP model relative to the specific model of Gershwin and given the ease of SIMP's computations. 4.0 Conclusions We have presented an extremely simple-minded, yet fast and apparently accurate approximation method for analyzing open tandem queues with blocking. The ability, to incorporate non-exponential service times presents an opportunity to study the behavior of more realistic systems than those considered previously, without the need to develop special procedures for idiosyncratic service distributions. 14

The conclusions about accuracy, of course, rest upon the evidence obtained to date of the approach's ability to produce, within acceptable percent differences, performance measures of interest for previously analyzed systems. It remains to be seen whether we can guarantee error bounds, or provide explicit advice on parameter values for which the approximation is unequivocally recommended. However, due to the simplicity of the approximation, its straight-forward structure, and its ready use for general service distributions, it holds promise to be a useful tool in the study of systems of queues. In particular, we are currently engaged in extending the general approach to closed tandem queues, and more general open networks. Appendix. Computation of Completion Time Moments Ti, the completion time for the ith server, is the weighted sum of two random variables: Si - (unblocked) service time for server i, and Xi - residual completion time for a randomly encountered (busy) server i, the latter weighted by ai = Pr{server i will encounter a "block" from server i + 1}. Let fi(t), hi (t), and gi (t) be the probability density functions, and fi*(s), hi (s), and g*(s) be the Laplace transforms, for the random variables Ti, Si, and Xi, respectively. Then f;*(s) = (1 - ci)h (s) + C a,h (s)g (s). (Al) From the moment-generating property of the Laplace transform, the nth moment of Ti is found by taking the nth derivative of fi (s) and evaluating it at s = O. Taking the nh derivative of both sides of (Al) (where f(n)(s) = df(a)) gives fi( ()() = (1- ai)h(n)" ) + ai ( h(ni)(- ) () (A2) j=1 Since Xi is the residual time for the random variable Ti, we have from Equation 7b E[T'n+] i+ ) E[X,+,]- (, + 1)E[Ti+,]' s8=0 (n+l i Thus, evaluating both sides of (A2) at s = 0 gives nrpr^ n \[Sin-j]E \ \ (j+ 1 E[T (1-E[ =(1- )[Sin] + a, E ( S E[T+1 (A3) = - ]E[15T+ 15

Since the j = 0 term in the sum is simply E[Si,], Equation A3 can also be written as Equation 20. In the special case where the Si are exponentially distributed with rate 1l, then E[S$ ] = " and A3 reduces to Equation 11. ~2T 16

TABLE 1 Throughput Comparisons (from Takahashi et al. 1980) JL 2 A3 Simul- Hillier- Takahashi SIMP SIMP SIMP ation Boling et al. Eq. (2) Eq. (12) M/G/1 1.1 1.2 1.3.709.712.645.689.691.694 1.2 1.4 1.6.762.767.706.746.748.751 1.3 1.6 1.9.798.807.755.790.790.793 1.4 1.8 2.2.828.837.793.823.824.827 1.5 2.0 2.5.855.861.824.850.850.853 1.6 2.2 2.8.877.880.849.871.872.874 1.7 2.4 3.1.889.896.870.889.889.891 1.8 2.6 3.4.902.909.887.904.904.906 1.9 2.8 3.7.915.920.901.915.915.917 2.0 3.0 4.0.929.929.913.925.926.927 Max. Abs. Deviation.000.009.064.020.018.015 Ave. Abs. Deviation.000.005.032.007.006.005 Ave. Abs. % Deviation.000 0.6 4.0 0.8 0.8 0.6 Arrivals to queue 1 constitute a Poisson process with rate 1; B, = 1 except B1 = 2. 17

TABLE 2 Comparisons with the Exact Solution and the Approximation of Altiok and Perros Case Measure Exact Altiok SIMP SIMP SIMP Solution and Perros Eq. (2) Eq. (12) M/G/1/N (B1,B2,B3)=(oo,2,2) P (0).281.299.298.299.299 (i1,t2,p3)=(2,2,2) P2(0).326.331.331.331.331 A=1.2 P2(3).262.237.239.239.237 P3 (0).400.400.400.400.400 P3(3).182.176.176.176.176 (B1,B2,B3)=(oo,l,1) P1 (0).166.173.167.175.175 (P,/12,P3)=(2,2,2) P2 (0).267.268.268.268.268 A=1.2 P2(2).483.474.477.476.472 P3 (0).400.400.400.400.400 P3 (2).323.323.323.323.323 (B1,B2,B3)=(oo,1,1) P1 (0).386.402.399.402.402 (/1,U2,/s3)=(3,2.5,3.5) P2 (0).396.396.396.396.396 A=1.4 P2(2).338.325.328.327.325 P3 (0).600.600.600.600.600 P3 (2).149.149.149.149.149 (B1,B2,B3)=(1,,1) Pi (0).118.127.131.131.128 (/i1,/2,P3)=(2,2,2) P1 (2).574.588.591.589.588 A=3 P2 (0).213.239.245.240.239 P2(2).524.512.507.513.512 P3 (0).361.382.386.383.382 P3 (2).356.342.338.341.343 (B1,B2,B3)=(1,,1) P (0).209.223.224.222.222 (/,1,2,/3 )=(3,4,2) P1 (2).451.459.456.459.459 A=3 P2 (0).243.274.268.276.274 P2(2).503.473.476.467.474 P3 (0).176.188.185.189.189 P3(2).608.589.593.586.587 Max. Abs. Deviation.000.031.032.036.031 Ave. Abs. Deviation.000.011.011.012.012 Ave. Abs. % Deviation.000 4.0 3.9 4.2 4.0 18

TABLE 3 Throughput Comparisons with Simulation and the Approximation of Altiok and Perros A B Simulation Altiok SIMP SIMP SIMP p_, (Exact) and Perros Eq. (2) Eq. (12) M/G/1/N 1.5 (1,1,2).916.851.939.930.928 (2,3,1) 2 (3,2,1) 1.780 1.553 1.669 1.672 1.668 (4,3,2) 2 (1,1,2).954.951.965.959.959 (2,3,1) 3 (4,3,2,1,1) 1.277 1.196 1.248 1.262 1.265 (22,2,2,22) 3 (1,1,2,2,2) 1.438 1.357 1.400 1.405 1.407 (3,3,2,2,2) 2.5 (1,2,3,2,1) 1.295 1.113 1.259 1.263 1.263 (2,2,2,2,2) 3 (3,3,2,2,1) 1.647 1.399 1.682 1.681 1.682 (2,3,4,3,2) 2 (1,2,3,2,1) (1.294) 1.344 1.207 1.208 1.208 (2,2,2,2,2) 2 (1,1,1,1,1) (1.131) 1.086 1.079 1.088 1.086 (2,2,2,2,2) Max. Abs. Deviation.000.248.111.108.112 Ave. Abs. Deviation.000.109.046.041.041 Ave. Abs. % Deviation.000 7.7 3.4 3.0 2.9 19

TABLE 4 Throughput Comparisons with the Approximation of Gershwin (1983) B '4 'I. I *4VI lrl i1 VI' IF wV s r4 nl Ir v,% (4,5,5) (4,5,5) All.002 All.1 (4,5,5) All.002 All.05 (4,5,5) All.1 All.1 (4,10,10) All.1 All.1 (4,10,10) (.01,.013,.007) (.07,.1,.05) (4,5,5,5) All.05 All.45 (4,10,10,10) All.1 All.1 (4,10,10,10,10) All.1 All.1 bimulatlon Gershwin IMP bIMP f IMa (Exact) (.4545) (.8993) (.3109) (.3556) (.7539).8153.3364.3204.4542.8993.3105.3563.7546.8716.3352.3228 Eq. (2).5519.7980.4139.4469.7864.7232.4393.4351 Eq. (12).5525.7969.4144.4466.7866.7243.4395.4354 M/G/1/N.4881.7591.3531.3939.7512.7716.3786.3697 I B1 = 4, Bi = 10(i = 2,3,...,12) (.01,.013,.007,001,.13,.007,.01,.013,.007,.01,.013,.007).5916.6098.7494.7495.6892 (.07,.1,.05,.07,.1,.05,.07,.1,.05,.07,.1,.05) B1 = 4, B = 10(i = 2,3,...,20) All.1.2767.2924.4243.4243.3498 All.1 (4,40,40,40) All.005.8362.8337.8754.8754.8484 All.05 Max. Abs. Deviation.0000.0182.1578.1579.1402 Ave. Abs. Deviation.0000.0040.0986.0983.0523 Ave. Abs. % Deviation.0000.52 23.52 23.51 11.40 The service time at queue i is probability pi ri. 1 period with probability 1 - pi and n periods with 20

References 1. Altiok, T., (1982) "Approximate analysis of exponential tandem queues with blocking", European Journal of Operations Research 11, 390-398. 2. Altiok, T., and Stidham, S., (1982) "A note on transfer lines with unreliable machines, random processing times and finite buffers", AIIE Transactions 14, 125-127. 3. Asare, B. K., (1978) "Queue networks with blocking", Ph.D. dissertation, Trinity College, Dublin, Ireland. 4. Boxma, 0. J., and Konheim, A. G., (1981) "Approximate analysis of exponential queueing systems with blocking", Acta Informatica 15, 19-66. 5. Buzacott, J. A., (1967) "Markov chain analysis of automatic transfer lines with buffer stock", Ph.D. dissertation, Department of Engineering Production, University of Birmingham. 6. Caseau, P., and Pujolle, G., (1979) "Throughput capacity of a sequence of queues with blocking due to finite waiting room", IEEE Transactions on Software Engineering SE5, 631-642. 7. Cooper, R. B., (1981) Introduction to Queueing Theory. 2nd ed. North-Holland, New York. 8. de Souza e Silva, E. Lavenberg, S. S., and Munty, R.R., (1983), "A perspective on iterative methods for the approximate analysis of closed queueing networks", IBM Report RC 10141. 9. Foster, F. G., and Perros, H. G., (1980) "On the blocking process in queue networks", European Journal of Operations Research 5, 276-283. 10. Gershwin, S. B., (1983) "An efficient decomposition method for the approximate evaluation of tandem queues with finite storage and blocking", manuscript, Laboratory for Information and Decision Sciences, MIT. 11. Gordon, W. J., and Newell, G. F., (1967) "Cyclic queueing systems with restricted length queues", Operations Research 15, 266-278. 21

12. Hillier, F. S., and Boling, R., (1967) "Finite queues in series with exponential or Erlang service times-A numerical approach", Operations Research 15, 286-303. 13. Hordijk, A., and Van Dijk, N., (1981) "Networks of queues with blocking", Performance '81, F. J. Kylstra (Ed.,), North-Holland, New York. 14. Konheim, A. G., and Reiser, M., (1976) "A queueing model with finite waiting room and blocking", Association for Computing Machinery 23, 328-341. 15. Labetoulle, J., and Pujolle, G., (1977) "Modelling of packet switching communication networks with finite buffer size at each node", Computer Performance, Chandy and Reiser (Eds.), North-Holland 515-535. 16. Latouche, G., and Neuts, M. F., (1980) "Efficient algorithmic solutions to exponential tandem queues with blocking", SIAM Journal of Algebraic Discrete Methods 1, 93-106. 17. Marchal, W. G., (1984) "Numerical Performance of Approximate Queueing Formulas with Application to Flexible Manufacturing Systems", University of Toledo. 18. Perros, H. G., (1981) "A symmetrical exponential open queue network with blocking and feedback", IEEE Transactions on Software Engineering Vol. SE-7, 395-402. 19. Perros, H. G., (1984) "Queueing networks with blocking: A bibliography". Performance Evaluation Review, A CM SIGMETRICS 12, 8-12. 20. Perros, H. G., and Altiok, T., "Approximate analysis of open networks of queues with blocking: Tandem configurations", IEEE Transactions on Software Engineering (forthcoming). 21. Pinedo, M., and Wolf, R. W., (1982) "A Comparison between tandem queues with dependent and independent service times", Operations Research 30, 464-479. 23. Suri, R., and Diehl, G. W., (1983) "A variable buffer-size model and its use in analytic closed queueing networks with blocking", manuscript, Division of Applied Science, Harvard University. 23. Takahashi, Y., Miyahara, H., and Hasegawa, T., (1980) "An approximation method for open restricted queueing networks", Operations Research 28, 594-602. 22

24. Towsley, D., (1979) "Queueing network models with state-dependent routing", Association for Computing Machinery 27, 323-337. 23