- Email: [email protected]
Single server stochastic recirculation systems
Cmnpur.Opns Res. Vol. 14, No. 1, pp. 75-84, 1987 Printed in Great B&in. AI1rights reserved
SINGLE
SERVER
0305-0548/87 $3.00+ 0.00 Copyright@ 1987Pergmon Journals Ltd
STOCHASTIC
RECIRCULATION
SYSTEMS
DAVID SONDERMAN~*** and BEHNAM POURBABAI~*~ ‘Department of Industrial Engineering and Operations Research, University of Massachusetts, Amherst, MA 01003 and %zpartment of Statistics and Operations Research, New York University, New York, NY 10006, U.S.A. (Receiued November 1984; in revised
furmJune 1985)
Scope and Parpase-Material handling systems play an important role in the performance of production systems. In this paper a queueing model is developed for the performance analysis of a production system consisting of a single processor, a loading tid an unloading stations linked by a closed loop material handling system. Our primary objective is to study the flows of units along the material handling system and those which leave the production system. We call such a system a single server stochastic recirculation system and develop a computationally e&ient methodology for studying its performance. Abstract-A recursive methodology is suggested for approximating the asymptotic performance of a single server recirculation system. A single server recirculation system is modeled by a G/G/l/O queueing loss system with generally distributed interarrival times, generally distributed service times, a single server, no waiting room, first come first served queueing discipline, and retrials. In this queueing system, the units which initially are not prooessed by the server are not lost. These units retry to receive service by retrying. Furthermore, numerical results are provided and the approbation outcomes are compared against those from a simulation study.
1. INTRODUCTION
AND OBJECTIVES
The broad objective of this paper is to examine the stochastic behavior of a single server, finite waiting room queueing system whose overflows are not lost to the system but rather return to try for service after a specific fixed period of time. We define such a system to be a single server stochastic recircdution system. Through the use of two parameter approximation techniques, the important stochastic processes in the system are examined for a variety of values of input and service parameters. When compared with the corr~~nding simulation runs, these appro~mation techniques are shown to be very accurate for predicting the- behavior of the system. However, they employ vastly less computational time, and are therefore useful in situations when simulations are not practical. This analysis can be used to provide answers to questions about design and operation of the system, either as a separate entity or as a node in a larger network. The results of this paper are extensions of a model proposed in Sonderman [ 1J for a recirculating conveyor in a materials handling situation. In that model it was assumed that both external interarrival times and service times were Gquences of independent, identically distributed exponential random variables. In this case both arrivals and service times are allowed more general distributions. The paper is organized in the following manner. In Section 2 we introduce the model and the necessary notation. In Section 3 we review the approximation techniques utilized to analyze the stochastic processes encounters in our model, and present an algo~t~ to calculate the steady-state behavior of the system. In Section 4 we examine the results of our model for a variety of parameter values, and compare its performance in predicting system parameters against a simulation of the system. We close in Section 5 with an example of a network of stochastic recirculation systems in tandem. *David Sonderman was an Assistant Professor in the Department of Industrial Engineering and Operations Research at the University of Massachusetts at Amherst. We earned a B.S. in IE/OR from Cornell University and a Ph.D. in Operations Research from Yale University. He has published in Advancesin AppliedProbability,Mathematicsof OperationsResearchand other journals. Dr Sonderman unexpectedly passed away in 1985. tBehnam Pourbabai is an Assistant Professor in the Department of Statistics and Operations Research at New York University. He earned a B.S. in Civil Engineering and an M.S. in Industrial and Operations Engineering from the University of Michigan at Ann Arbor, and a Ph.D. in Operations Research from the University of Massachusetts at Amherst. He has published in the Journal of AppliedPr~~i~~y, European J~r~al of O~r~ions Research, l~ter~t~~l J~r~l of Systems Science, ~~~i~ns Re~~ch Letter, M~~sicai M~elli~, Jnter~i~al Journal of S~l~i~n and M~elli~ and Computersand Mohair. 75
DAVID SONDERMAN and BEHNAM POURBABAI
Superpositiw
arrival
process
(Aa, c,) Primary
arrival
(A, Overflow
Fig. 1. A stochastic
2.
STATEMENT
I Cd process
recirculation
OF
THE
system
MODEL
In this section we introduce the notation for our model, and present a model of the general stochastic recirculation system. In the following sections we examine both an analytical approximation and a simulation model of the system. Consider Fig. 1. Units entering the system (also called the primary input to the system) do so according to a stochastic process {Aj, j 2 l}, where Aj = interarrival time between (j - 1)st and jth unit, j 2 1. These units are transported for a fixed period of time to the service station. At the service subsystem, units finding the server free are served by a mechanism, the service time being a random variable. Let Sj = time required to service the jth unit available for service, j 2 1. The server will of course not necessarily be free, and we will assume there are no waiting spaces available. If a unit arrives at the server to find it busy, the unit overflows, being transported for a fixed period of time back to the source of primary input. At this point the stream of recirculating units superimposes itself with the primary input process. At this point there are no assumptions regarding the physical size of the system Fig. 1 is only a schematic representation. We need some assumptions regarding the stochastic nature of the arrival and service processes. It is assumed that the arrival process is a general renewal process whose interevent times are characterized by their first two moments, and that the service times are independent, identically distributed random variables, also characterized by their first two moments. Later we will see that it is possible to relax the renewal assumption. Because, the input process can also be approximated by a stationary counting process, as will be demonstrated in the following section. Let (i,,, c,) = the rate of the primary input process and the coefficient of variation of the associated interarrival times, respectively, and (pb, cb) = the rate of the service process and the coefficient of variation of the service times, respectively. The nature of the two different stochastic processes along the top and bottom legs of the system will be of primary interest to us in the next three sections. In our analytical model we will be approximating these general counting processes by renewal processes, which in turn are characterized by their first two moments. We start with an empty system, and define one cycle to
Single server stochastic recirculation systems
77
consist of the total time that would be required for a unit to go from the source of primary input to the service subsystem and back to the source of primary input, assuming it is refused service. To describe the approximating renewal processes we need to define the following distributions: G,(t) = the cumulative distribution function of the interevent times of the approximating renewal process along the top (i = 1) and bottom (i = 2) legs, as the kth cycle, k=l,2,.... Throughout this paper we will generically use G*(s) to denote the Laplace Stieltjes transform of any c.d.f. G(t). Let
s a,
G*(s) =
exp( -st) dG(t).
0
From the Laplace Stieltjes transform we can derive the rate of the process and the coefficient of variation of the interevent times as follows. i,, = (-G%(O))-’ cik =
&(Gr(O)
(2)
(G$'(0))2)"2.
-
The details of the approximation techniques utilized in our analysis are discussed in the following section. At this junction it will be sufftcient to define two distributions that will be of use in the approximations. Although these distributions may be used to approximate the service times and/or the times between events on either leg on any cycle, for convenience we suppress the i, k notation. The first distribution, used to approximate distributions with a coefficient of variation greater than or equal to one, is the two-parameter hyperexponential distribution (H2) with balanced means. The form of the density is
dH,(p, dt
a;t) =Pexp(--&t)+
(1 -p)8, exp(-8,t),
tao
where PA
=
(1
-
(4)
P)P**
The parameters p and e1 are referred to as the shape and intensity parameters, respectively. The second distribution, used to approximate distributions with a coefficient of variation less than one, is the delayed exponential distribution M’(d, /?; t). The form of the density is
d”‘fi” t,=j3 exp( - fi(t -
d)),
t>dBO.
The parameters d and /? are referred to as the scale and intensity parameters, respectively.
3. APPROXIMATION
TECHNIQUES
Our primary task is to approximate the behavior of the stochastic processes of the system. We extend a method developed in Sonderman [l] to approximate the rate and coefficient of variation of the stochastic processes along the two legs of the system. This two parameter approximation procedure worked well in Sonderman [ 11, and our simulation results indicate it also performs well in this case with more general assumptions. Consider the queueing subsystem formed by the single server. This subsystem is treated as a GZ/H,/l/l queueing system. The main interests are to investigate the nature of the departure and
DAVID SONDERMANand BEHNAMPOURBABAI
78
overflow processes from the queueing subsystem, and then the superposition of the overflow process with the primary input. Naturally these latter two processes of interest interact with one another; this interaction is considered in our model. Whitt [2,3] introduced the following procedure to approximate a stationary counting process by a renewal process, given that certain regularity conditions are satisfied. Here we assume the processes along both legs can be represented by stationary counting processes. The methodology is to find several moments of the times between the renewal events and then fit a convenient distribution to those moments. To know what distribution is appropriate, the value of the coefficient of variation of these times must be known. With a coefficient of variation greater than one it is recommmended by Whitt [2,3] to approximate the interevent distribution by a hyperexponential distribution with balanced means, while with a coefficient of variation less than one it is recommended to approximate the interevent distribution by a delayed exponential distribution. This approximation technique affects both the behavior of the queueing subsystem as well as the superposition of processes at the source of primary input. We first consider the queueing subsystem. Since we will consider interarrival distributions with coefficients of variation either greater or less than one and service distributions with coefficients of variation greater than one, in general we need a method to handle the overflows from a GI/H,/l/l queueing system. The methodology employed is based on the results of Haltin [4], which extends the method first developed by Palm [5,6]. We will need expressions for the first two moments of the interoverflow processes from two special versions of the GZ/H,/l/l subsystem. These expressions are summarized in the following proposition.
Proposition
1
Consider the following GI/H,/l/l
queueing systems.
(4 H21HzIlIl
@I M’IH,IlIl. Arrivals form a sequence of @dependent, identically distributed random variables with c.d.f. A(t) = H2(p, 8; t) (or M’(d, fi; t)) and service times form a sequence of i.i.d. random variables with c.d.f. G(t) E H,(q, c$; t). Then the first two moments of the interoverflow distribution with Laplace Stieltjes transform G*(s) are:
- G*‘(O)= -A*‘(O) G*“(O)
= A*“(O)
1
- +
(6)
(7)
where x = [2N*(O)l-*‘(O)A*‘(O)+ N*(O)r*(O)A*“(O) - r*(O)A*“(O) -2r*yo)A*yo)]p x
[N*(o)r*(o)A*qo)
Y = [l -
- w(o)12 + 2~*yo)[i - r*(o)A*yo)]
N*(o)]4
- N*(O)]
(8) (9)
and for (a) HZ/Hz/l/l
(10)
~- 01
N*(s)=pq [ (e, +s)
0, (e, +f#Q +s) 1
(11)
Single server stochastic recirculation systems
82
-- 02
+(l-p)q [ (e, +s) ~- 4
+
(12)
(e, +r#Jt +s) 1
4 (e,++,+s)
+p(l-q) [ (et+4
[
79
1
e2 ~- e2 + S) (e, + 4z +
(’- p)U - q,(e,
s) 1
where
r*w=(1 +~)~*(s)-(~)[~*(s)+r~(s)l qA*(s+d4
r:(S)=
A*(($)
1_
+
u--qM*(s+4,) 1 - A*(44
(13)
(l-4) P=l-J*(f$r)+l-A*($,)-l and for (b) M’/H,/l/l N*(s)
=
B exp( - sd) @+s)
-9(8+Sl+~~)exPl-d(s+d,)l
-(1-4)(8+sB,m,)exP[-d(.+m,)l
(14)
Proof
From Haltin [4] the Laplace Stieltjes transform G*(s) of the interoverflow distribution from a GZ/G/l/l system with interarrival (service) distribution A(t)@(t)) is
G*(s) = ,4*(s) -
I-*(s)[l - A*(s)] 1 -N*(s) (
(15)
where
(16)
’B(x) d/t(x)
(17)
[l-B(x)] dU,(x)
(19)
N(t) =
s0
P=
s ;;
g = sup@: B(t) < 1).
(20)
The desired moments are achieved by substitution and differentiation. The results of Proposition 1 make it possible to approximate the behavior of the stochastic process along the bottom leg of the system, given the behavior of the stochastic process along the top leg.
DAVIDSONDERMAN and BEHNAMPOURBABAI
80
However, this latter process is determined by the superposition of the primary input process with the process along the bottom leg. We need a method to approximate this superposition of stochastic processes. The approximation procedure which is used in this paper is the stationary interval method of Whitt [2,3]. In this method, the dependence among events is considered. This method is selected since the dependence among events in steady state is of primary interest to us here. As was shown in Proposition 1, we can obtain an approximation for the first two moments of the overflow process from the queueing subsystem along the bottom leg. To obtain the arrival rate of the superposition arrival process at any cycle we add the arrival rate of the primary input process to the rate of the overflow process from the queueing subsystem at the previous cycle:
(21) To find the coefficient of variation of the superposition
arrival process we use
(22) The first term represents the contribution
to variability from the primary input process, and the
The primary arrival process with (A., , Co b
Write
h,
Yes
,qk I 1
for
i = 1 and 2
1
i. A$woximata
Approximate the process along the top
the process along the top
leg by a renewal prcces with
leg by a renewal process with indspendsnt and delayed exponentially
independent and hypeexponentially
dktributd
distributi interevent times with balanced mean*
interwent times 2A
1 FindlXzk prcc&
, Czk1from
26
k I Find (x2 k,
the overflow
C2k )
from the overflow
process from an Hz / Hz /l
from a M’ /Hz / 1 I 1
queueing svst~m
qususing system
0 k-k+
1
the superposition of the primary input with overflows from the quweing system
Fig. 2. Flowchart of the approximation algorithm.
I 1
Single server
stochasticrecirculationsystems
81
second term represents the contribution from the process representing leg 2 at the previous cycle. We refer the interested reader to Whitt [2,3] for details. An overview of the composite approximation technique is provided in Fig. 2. We start in box 1 by approximating the stationary counting process along the top leg by a renewal process whose interevent distribution has the same first two moments as the counting process. Branching to either box 2A (cp < 1) or box 2B (ca, 2 l), we analyze the behavior at the queueing subsystem using Proposition 1. We should note that given the results of box 1, there is no further approximation in boxes 2A or 2B. We finish the cycle by approximating the superposition of the overflow process with the primary input in box 3, after which we return to box 1. In our tests of the method we started the system with no units in the system, and went through a number of cycles until a steady state situation was achieved. The decision node before box 2 is necessary because during any one application of the algorithm it is possible to use both boxes a number of different times on the way to steady state. An algorithm based on Fig. 2 was coded, and run for various service and primary input rates and coeffkients of variation. A comparison of the algorithm with a simulation is discussed in Section 5. 4. SIMULATION
RESULTS
In this section we compare the results of our approximation algorithm with a simulation of the stochastic recirculation system for various values of input and service parameters. A summary of both the simulation and algorithm results for these various parameter values is provided in Fig. 3. The algorithm results were obtained by starting with an empty system, and applying the algorithm until a steady-state was reached. The simulation results were obtained by averaging ten different runs, each. run lasting until the steady state was achieved. The numbers in parentheses represent the standard deviations of the ten different runs. Consider Fig. 3. Comparisons can be drawn between the rates and coefficients of variation along the top and bottom legs of the system for low (p = 0.3), medium (p = 0.5) and high (p = 0.7) traffic intensities. Two types of service discipline were considered: exponential (cb = 1) and general hyperexponential (cb = 1.5). Comparing the simulation with the approximation shows an excellent match for the rates i,,, &. Results for the coeffkients of variation cl, c2, although perhaps less important to the design and operation of the system, are also quite satisfactory. All the steady-state results display equilibrium between the rate of overflows from the server and the rate of superposition of processes at the source of primary input. Consider the parameter values I, = 0.7, c, = La = cb = 1. A rate of R, = 3.19 along the top leg, along with the approximated cl, yields an ovefflow stream with a
0.5 0.5 0.5
1 1 1
1 1 1
0.33 0.82 2.55
0.34(0.00)
0.57 0.77 1.13
0.58(0.00)
0.86(0.02) 2.50(0.08)
0.77(0.01) 1.07(0.02)
0.03 0.32 1.85
0.05(0.00) 0.36(0.02) l.El(O.08)
1.00 1.06 1.29
l.OO(O.08) 1.17(0.03) 1.38(0.00)
1 1
1 1
1 1
0.44
1
11
1.11 3.19
0.45(0.01) 1.15(0.03) 3.32(0.15)
1.08 1.20 1.44
1.07(0.10) l.lE(O.04) 1.40(0.11)
0.14 0.61 2.49
0.15(0.01) 0.65(0.03) 2.63(0.15)
1.22 1.35 1.55
1.28(0.02) 1.46(0.04) 1.69(0.05)
fX 0:7
1.5 1.5 1.5
1 1 1
1 1 1
0.49 1.31 4.11
0.48(0.01) 1.25(0.03) 3.41(0.16)
1.47 1.56 1.82
1.35(0.02) 1.39(0.03) 1.54(0.03)
0.19 0.81 3.41
0.19(0.01) 0.76(0.02) 2.72(0.15)
1.41 1.60 1.88
1.46(0.02) 1.64(0.05) l.El(O.03)
0.3 0.5 0.7
0.5 0.5 0.5
1 1 1
1.5 1.5 1.5
0.35 0.91 2.61
0.36(0.00) 0.90(0.00) 2.61(0.11)
0.65 0.89 1.17
0.61(0.01) O.EOfO.01) 1 .i 9io.27j
0.05 0.41 1.91
O.OS(O.00) 0.4310.091 i.92io.i 1j
1.20 1.21 1.34
1.22(0.05) 1.36(0.03) 1.50(0.04)
E.35 0:7
1 1 1
1 1 1
1.5 1.5 1.5
0.44 1.11 3.11
0.44(0.01) 1.16(0.04) 3.26(0.18)
1.10 1.24 1.45
1.01(0.01) 1.22(0.02) 1.43(0.04)
0.14 0.61 2.41
0.15(0.01) 0.67(0.03) 2.58(0.20)
1.29 1.40 1.55
1.35(0.03) 1.55(0.03) 1.74(0.07)
0.3 0.5 0.7
1.5 1.5 1.5
1 1 1
1.5 1.5 1.5
0.48 1.28 3.85
0.46(0.01) 1.26(0.06) 3.51(0.26)
1.49 1.58 1.79
1.35(0.03) 1.43(0.02) 1.57(0.04)
0.18 0.78 3.15
O.lE(O.01) 0.77(0.06) 2.81(0.25)
1.47 1.63 1.85
1.51(0.03) 1.71(0.03) 1.86(0.04)
0.3 0.3 0.7 8:: 0.7
*An approximation outcome. tA simulationoutcome.
Fig. 3. Comparison of approximation and simulation outcomes.
DAVID SONDERMANand BEHNAMPOURBABAI
82
I2 22,1.441
f2.49,1.55)
(2.19,1431
Fig. 4. Three recirculating conveyors.
rate of i,, = 2.49, and an approximated c2. In turn, these values of (L,, cJ superimpose with the primary input stream to form a process with the values (n,, cl). This demonstrates how a steady-state can be reached even when the tragic intensity at the server is greater than one. This is similar to an ordinary queueing system with a finite waiting room, but in this model we have explicitly dealt with the overflows. 5. APPLICATION:
MATERIALS
HANDLING
SYSTEMS
IN SERIES
In this section, a tandem recirculation system is considered. The system consists of three single server recirculating systems. Each recirculating system could represent a production system consisting of a single processor, a loading and an unloading stations linked by a material handling system. The material handling system could be a closed loop conveyor system. The overall behavior of the system is of interest to us in this section. Consider Fig. 4, where we display a schematic diagram of a sequence of 3 stochastic recirculation systems in tandem, as might be found in an automated, flexible materials handling setting. The departure process from the first system becomes the arrival process to the second system, and so forth. It is of interest here to characterize the overall stochastic behavior of such a system, including the distribution of the number of items on the various individual conveyors. We present results for the case of equal exponential service rates at each service subsystem. We first need a result concerning the departure process from a single server queueing system with exponential service times and no waiting room. Proposition
2
Consider an Hz/M/l/1 queueing system. Arrivals form a sequence of i.i.d. random variables with c.d.f. H,(p, fi; t) and the one server is exponential at rate p. Then the first two moments of the interdeparture distribution with laplace stieltjes transform D*(s) are: (23)
. ..(.,,;li+,,,k,(~-~)+2k~k~(~-~)
(24)
k, =
(25)
where --r!F[l + m*(P)] ( P +e, )
k2=H k 3=
-
( eI
2 [l + m*(P)] p )
(26)
(27)
Single server stochastic recirculation systems
k4=
83
(-l-- > 62
P-
m*(P)=
WP, & P) 1 - fm, 6;PO’
(29)
Proof
The proposition is a special case of equation (2.4) of Daley [7]. The results of this proposition are combined with our algorithm to approximate the behavior of the stochastic processes representing units on both legs of each conveyor. In Fig. 4 we have displayed the results of our algorithm for the input parameters i,, = 0.7, c, = 1. Starting with Poisson input, we note that the departure process from the first conveyor has less variability than the input process (c = 0.82). Using this departure process as input to the second conveyor, the same phenomenon occurs to a lesser degree with the second departure process (c = 0.80). At this point we reach an equilibrium state where an input of (0.7, 0.80) yields a departure process with the same two parameters. In other words, the conveyors in series create a situation where the stochastic behavior of the departure processes from the individual conveyors becomes less variable as the units proceed “downstream”, and at some point an equilibrium state is reached. It will also be of interest to examine the overall congestion pattern in the series of conveyors. To do so we need to assign some physical characteristics to the individual systems. We assume that the amount of time required for a unit to transit each leg of each conveyor is 100 time units, although in general these distances may all be different. We then apply Smith’s [8] central limit theorem for renewal processes to calculate the approximate distribution of the number of items to be found on each conveyor. Let N,,(t) = number of items found on leg i of conveyor j, during time period t, i = 1,2, j = 1,2, . . . , t 2 0. Then for sufficiently large t, N,,(t) N Normal&
= i&t, a$ = lijtC2,)
where (n,, cij) are the parameters of the processes along leg i of conveyor j. For example, consider conveyor 1. It can be seen that N,,(lOO)-Normal(~,,
=319,0:,
=25.7’)
N,, (100) - Normal(p, 1 = 249,o: 1 = 24.52). Defining Nj(t)
=
Nij(t)
+
N2j(t)
we obtain N,(lOO) - Normal(pi = 568,~: = 35.5’). In Fig. 4, we have labeled the parameters of the various distributions of the number of items in the systems along with the appropriate 95 ‘A Confidence Intervals (CIs). The lessening of variability as the units move down the series of conveyors is also evident in the distributions of the various N,(t). The complete characterization of the system in this manner enables one to make changes in the design or operation of the system.
Acknowledgement-The
authors would like to acknowledge the financial support of NSF Grant ECS-8105954.
84
DAVIDSONDERMAN and BEHNAMPOURLIABAI REFERENCES
1. D. Sonderman, An analytical model for recirculating conveyors with stochastic input and outputs. I.J.P.R. 20, 591-605 (1982). 2. W. Whitt, Approximating a point process by a renewal process, I: a generate framework. Bell Laboratories, Holmdel, N.J. (1980). 3. W. Whitt, Approximating a point process by a renewal process, I: two basic methods. Opns Res. 30, 125-147 (1982). 4. S. Hallin, Distribution of the interoverflow time for the GI/G/l loss system. Math. Opns Res. 6, 563-569 (1981). 5. C. Palm, Intensity fluctuations in telephone traffic. Ericsson 7&h. 44, l-189 (1943). 6. C. Palm, Research on telephone traffic carried by full availability groups. Tele. No. 1, I-107 (1957). 7. D. J. Daley, Characterizing pure loss GI/G/l queues with renewal output. Proc. Camb. Phil. Sot. 75, 103-107 (1974). 8. W. L. Smith, Renewal theory and its ramifications. JI R. Statist. Sot. Ser. B 20, 243-302 (1958).
SINGLE
SERVER
0305-0548/87 $3.00+ 0.00 Copyright@ 1987Pergmon Journals Ltd
STOCHASTIC
RECIRCULATION
SYSTEMS
DAVID SONDERMAN~*** and BEHNAM POURBABAI~*~ ‘Department of Industrial Engineering and Operations Research, University of Massachusetts, Amherst, MA 01003 and %zpartment of Statistics and Operations Research, New York University, New York, NY 10006, U.S.A. (Receiued November 1984; in revised
furmJune 1985)
Scope and Parpase-Material handling systems play an important role in the performance of production systems. In this paper a queueing model is developed for the performance analysis of a production system consisting of a single processor, a loading tid an unloading stations linked by a closed loop material handling system. Our primary objective is to study the flows of units along the material handling system and those which leave the production system. We call such a system a single server stochastic recirculation system and develop a computationally e&ient methodology for studying its performance. Abstract-A recursive methodology is suggested for approximating the asymptotic performance of a single server recirculation system. A single server recirculation system is modeled by a G/G/l/O queueing loss system with generally distributed interarrival times, generally distributed service times, a single server, no waiting room, first come first served queueing discipline, and retrials. In this queueing system, the units which initially are not prooessed by the server are not lost. These units retry to receive service by retrying. Furthermore, numerical results are provided and the approbation outcomes are compared against those from a simulation study.
1. INTRODUCTION
AND OBJECTIVES
The broad objective of this paper is to examine the stochastic behavior of a single server, finite waiting room queueing system whose overflows are not lost to the system but rather return to try for service after a specific fixed period of time. We define such a system to be a single server stochastic recircdution system. Through the use of two parameter approximation techniques, the important stochastic processes in the system are examined for a variety of values of input and service parameters. When compared with the corr~~nding simulation runs, these appro~mation techniques are shown to be very accurate for predicting the- behavior of the system. However, they employ vastly less computational time, and are therefore useful in situations when simulations are not practical. This analysis can be used to provide answers to questions about design and operation of the system, either as a separate entity or as a node in a larger network. The results of this paper are extensions of a model proposed in Sonderman [ 1J for a recirculating conveyor in a materials handling situation. In that model it was assumed that both external interarrival times and service times were Gquences of independent, identically distributed exponential random variables. In this case both arrivals and service times are allowed more general distributions. The paper is organized in the following manner. In Section 2 we introduce the model and the necessary notation. In Section 3 we review the approximation techniques utilized to analyze the stochastic processes encounters in our model, and present an algo~t~ to calculate the steady-state behavior of the system. In Section 4 we examine the results of our model for a variety of parameter values, and compare its performance in predicting system parameters against a simulation of the system. We close in Section 5 with an example of a network of stochastic recirculation systems in tandem. *David Sonderman was an Assistant Professor in the Department of Industrial Engineering and Operations Research at the University of Massachusetts at Amherst. We earned a B.S. in IE/OR from Cornell University and a Ph.D. in Operations Research from Yale University. He has published in Advancesin AppliedProbability,Mathematicsof OperationsResearchand other journals. Dr Sonderman unexpectedly passed away in 1985. tBehnam Pourbabai is an Assistant Professor in the Department of Statistics and Operations Research at New York University. He earned a B.S. in Civil Engineering and an M.S. in Industrial and Operations Engineering from the University of Michigan at Ann Arbor, and a Ph.D. in Operations Research from the University of Massachusetts at Amherst. He has published in the Journal of AppliedPr~~i~~y, European J~r~al of O~r~ions Research, l~ter~t~~l J~r~l of Systems Science, ~~~i~ns Re~~ch Letter, M~~sicai M~elli~, Jnter~i~al Journal of S~l~i~n and M~elli~ and Computersand Mohair. 75
DAVID SONDERMAN and BEHNAM POURBABAI
Superpositiw
arrival
process
(Aa, c,) Primary
arrival
(A, Overflow
Fig. 1. A stochastic
2.
STATEMENT
I Cd process
recirculation
OF
THE
system
MODEL
In this section we introduce the notation for our model, and present a model of the general stochastic recirculation system. In the following sections we examine both an analytical approximation and a simulation model of the system. Consider Fig. 1. Units entering the system (also called the primary input to the system) do so according to a stochastic process {Aj, j 2 l}, where Aj = interarrival time between (j - 1)st and jth unit, j 2 1. These units are transported for a fixed period of time to the service station. At the service subsystem, units finding the server free are served by a mechanism, the service time being a random variable. Let Sj = time required to service the jth unit available for service, j 2 1. The server will of course not necessarily be free, and we will assume there are no waiting spaces available. If a unit arrives at the server to find it busy, the unit overflows, being transported for a fixed period of time back to the source of primary input. At this point the stream of recirculating units superimposes itself with the primary input process. At this point there are no assumptions regarding the physical size of the system Fig. 1 is only a schematic representation. We need some assumptions regarding the stochastic nature of the arrival and service processes. It is assumed that the arrival process is a general renewal process whose interevent times are characterized by their first two moments, and that the service times are independent, identically distributed random variables, also characterized by their first two moments. Later we will see that it is possible to relax the renewal assumption. Because, the input process can also be approximated by a stationary counting process, as will be demonstrated in the following section. Let (i,,, c,) = the rate of the primary input process and the coefficient of variation of the associated interarrival times, respectively, and (pb, cb) = the rate of the service process and the coefficient of variation of the service times, respectively. The nature of the two different stochastic processes along the top and bottom legs of the system will be of primary interest to us in the next three sections. In our analytical model we will be approximating these general counting processes by renewal processes, which in turn are characterized by their first two moments. We start with an empty system, and define one cycle to
Single server stochastic recirculation systems
77
consist of the total time that would be required for a unit to go from the source of primary input to the service subsystem and back to the source of primary input, assuming it is refused service. To describe the approximating renewal processes we need to define the following distributions: G,(t) = the cumulative distribution function of the interevent times of the approximating renewal process along the top (i = 1) and bottom (i = 2) legs, as the kth cycle, k=l,2,.... Throughout this paper we will generically use G*(s) to denote the Laplace Stieltjes transform of any c.d.f. G(t). Let
s a,
G*(s) =
exp( -st) dG(t).
0
From the Laplace Stieltjes transform we can derive the rate of the process and the coefficient of variation of the interevent times as follows. i,, = (-G%(O))-’ cik =
&(Gr(O)
(2)
(G$'(0))2)"2.
-
The details of the approximation techniques utilized in our analysis are discussed in the following section. At this junction it will be sufftcient to define two distributions that will be of use in the approximations. Although these distributions may be used to approximate the service times and/or the times between events on either leg on any cycle, for convenience we suppress the i, k notation. The first distribution, used to approximate distributions with a coefficient of variation greater than or equal to one, is the two-parameter hyperexponential distribution (H2) with balanced means. The form of the density is
dH,(p, dt
a;t) =Pexp(--&t)+
(1 -p)8, exp(-8,t),
tao
where PA
=
(1
-
(4)
P)P**
The parameters p and e1 are referred to as the shape and intensity parameters, respectively. The second distribution, used to approximate distributions with a coefficient of variation less than one, is the delayed exponential distribution M’(d, /?; t). The form of the density is
d”‘fi” t,=j3 exp( - fi(t -
d)),
t>dBO.
The parameters d and /? are referred to as the scale and intensity parameters, respectively.
3. APPROXIMATION
TECHNIQUES
Our primary task is to approximate the behavior of the stochastic processes of the system. We extend a method developed in Sonderman [l] to approximate the rate and coefficient of variation of the stochastic processes along the two legs of the system. This two parameter approximation procedure worked well in Sonderman [ 11, and our simulation results indicate it also performs well in this case with more general assumptions. Consider the queueing subsystem formed by the single server. This subsystem is treated as a GZ/H,/l/l queueing system. The main interests are to investigate the nature of the departure and
DAVID SONDERMANand BEHNAMPOURBABAI
78
overflow processes from the queueing subsystem, and then the superposition of the overflow process with the primary input. Naturally these latter two processes of interest interact with one another; this interaction is considered in our model. Whitt [2,3] introduced the following procedure to approximate a stationary counting process by a renewal process, given that certain regularity conditions are satisfied. Here we assume the processes along both legs can be represented by stationary counting processes. The methodology is to find several moments of the times between the renewal events and then fit a convenient distribution to those moments. To know what distribution is appropriate, the value of the coefficient of variation of these times must be known. With a coefficient of variation greater than one it is recommmended by Whitt [2,3] to approximate the interevent distribution by a hyperexponential distribution with balanced means, while with a coefficient of variation less than one it is recommended to approximate the interevent distribution by a delayed exponential distribution. This approximation technique affects both the behavior of the queueing subsystem as well as the superposition of processes at the source of primary input. We first consider the queueing subsystem. Since we will consider interarrival distributions with coefficients of variation either greater or less than one and service distributions with coefficients of variation greater than one, in general we need a method to handle the overflows from a GI/H,/l/l queueing system. The methodology employed is based on the results of Haltin [4], which extends the method first developed by Palm [5,6]. We will need expressions for the first two moments of the interoverflow processes from two special versions of the GZ/H,/l/l subsystem. These expressions are summarized in the following proposition.
Proposition
1
Consider the following GI/H,/l/l
queueing systems.
(4 H21HzIlIl
@I M’IH,IlIl. Arrivals form a sequence of @dependent, identically distributed random variables with c.d.f. A(t) = H2(p, 8; t) (or M’(d, fi; t)) and service times form a sequence of i.i.d. random variables with c.d.f. G(t) E H,(q, c$; t). Then the first two moments of the interoverflow distribution with Laplace Stieltjes transform G*(s) are:
- G*‘(O)= -A*‘(O) G*“(O)
= A*“(O)
1
- +
(6)
(7)
where x = [2N*(O)l-*‘(O)A*‘(O)+ N*(O)r*(O)A*“(O) - r*(O)A*“(O) -2r*yo)A*yo)]p x
[N*(o)r*(o)A*qo)
Y = [l -
- w(o)12 + 2~*yo)[i - r*(o)A*yo)]
N*(o)]4
- N*(O)]
(8) (9)
and for (a) HZ/Hz/l/l
(10)
~- 01
N*(s)=pq [ (e, +s)
0, (e, +f#Q +s) 1
(11)
Single server stochastic recirculation systems
82
-- 02
+(l-p)q [ (e, +s) ~- 4
+
(12)
(e, +r#Jt +s) 1
4 (e,++,+s)
+p(l-q) [ (et+4
[
79
1
e2 ~- e2 + S) (e, + 4z +
(’- p)U - q,(e,
s) 1
where
r*w=(1 +~)~*(s)-(~)[~*(s)+r~(s)l qA*(s+d4
r:(S)=
A*(($)
1_
+
u--qM*(s+4,) 1 - A*(44
(13)
(l-4) P=l-J*(f$r)+l-A*($,)-l and for (b) M’/H,/l/l N*(s)
=
B exp( - sd) @+s)
-9(8+Sl+~~)exPl-d(s+d,)l
-(1-4)(8+sB,m,)exP[-d(.+m,)l
(14)
Proof
From Haltin [4] the Laplace Stieltjes transform G*(s) of the interoverflow distribution from a GZ/G/l/l system with interarrival (service) distribution A(t)@(t)) is
G*(s) = ,4*(s) -
I-*(s)[l - A*(s)] 1 -N*(s) (
(15)
where
(16)
’B(x) d/t(x)
(17)
[l-B(x)] dU,(x)
(19)
N(t) =
s0
P=
s ;;
g = sup@: B(t) < 1).
(20)
The desired moments are achieved by substitution and differentiation. The results of Proposition 1 make it possible to approximate the behavior of the stochastic process along the bottom leg of the system, given the behavior of the stochastic process along the top leg.
DAVIDSONDERMAN and BEHNAMPOURBABAI
80
However, this latter process is determined by the superposition of the primary input process with the process along the bottom leg. We need a method to approximate this superposition of stochastic processes. The approximation procedure which is used in this paper is the stationary interval method of Whitt [2,3]. In this method, the dependence among events is considered. This method is selected since the dependence among events in steady state is of primary interest to us here. As was shown in Proposition 1, we can obtain an approximation for the first two moments of the overflow process from the queueing subsystem along the bottom leg. To obtain the arrival rate of the superposition arrival process at any cycle we add the arrival rate of the primary input process to the rate of the overflow process from the queueing subsystem at the previous cycle:
(21) To find the coefficient of variation of the superposition
arrival process we use
(22) The first term represents the contribution
to variability from the primary input process, and the
The primary arrival process with (A., , Co b
Write
h,
Yes
,qk I 1
for
i = 1 and 2
1
i. A$woximata
Approximate the process along the top
the process along the top
leg by a renewal prcces with
leg by a renewal process with indspendsnt and delayed exponentially
independent and hypeexponentially
dktributd
distributi interevent times with balanced mean*
interwent times 2A
1 FindlXzk prcc&
, Czk1from
26
k I Find (x2 k,
the overflow
C2k )
from the overflow
process from an Hz / Hz /l
from a M’ /Hz / 1 I 1
queueing svst~m
qususing system
0 k-k+
1
the superposition of the primary input with overflows from the quweing system
Fig. 2. Flowchart of the approximation algorithm.
I 1
Single server
stochasticrecirculationsystems
81
second term represents the contribution from the process representing leg 2 at the previous cycle. We refer the interested reader to Whitt [2,3] for details. An overview of the composite approximation technique is provided in Fig. 2. We start in box 1 by approximating the stationary counting process along the top leg by a renewal process whose interevent distribution has the same first two moments as the counting process. Branching to either box 2A (cp < 1) or box 2B (ca, 2 l), we analyze the behavior at the queueing subsystem using Proposition 1. We should note that given the results of box 1, there is no further approximation in boxes 2A or 2B. We finish the cycle by approximating the superposition of the overflow process with the primary input in box 3, after which we return to box 1. In our tests of the method we started the system with no units in the system, and went through a number of cycles until a steady state situation was achieved. The decision node before box 2 is necessary because during any one application of the algorithm it is possible to use both boxes a number of different times on the way to steady state. An algorithm based on Fig. 2 was coded, and run for various service and primary input rates and coeffkients of variation. A comparison of the algorithm with a simulation is discussed in Section 5. 4. SIMULATION
RESULTS
In this section we compare the results of our approximation algorithm with a simulation of the stochastic recirculation system for various values of input and service parameters. A summary of both the simulation and algorithm results for these various parameter values is provided in Fig. 3. The algorithm results were obtained by starting with an empty system, and applying the algorithm until a steady-state was reached. The simulation results were obtained by averaging ten different runs, each. run lasting until the steady state was achieved. The numbers in parentheses represent the standard deviations of the ten different runs. Consider Fig. 3. Comparisons can be drawn between the rates and coefficients of variation along the top and bottom legs of the system for low (p = 0.3), medium (p = 0.5) and high (p = 0.7) traffic intensities. Two types of service discipline were considered: exponential (cb = 1) and general hyperexponential (cb = 1.5). Comparing the simulation with the approximation shows an excellent match for the rates i,,, &. Results for the coeffkients of variation cl, c2, although perhaps less important to the design and operation of the system, are also quite satisfactory. All the steady-state results display equilibrium between the rate of overflows from the server and the rate of superposition of processes at the source of primary input. Consider the parameter values I, = 0.7, c, = La = cb = 1. A rate of R, = 3.19 along the top leg, along with the approximated cl, yields an ovefflow stream with a
0.5 0.5 0.5
1 1 1
1 1 1
0.33 0.82 2.55
0.34(0.00)
0.57 0.77 1.13
0.58(0.00)
0.86(0.02) 2.50(0.08)
0.77(0.01) 1.07(0.02)
0.03 0.32 1.85
0.05(0.00) 0.36(0.02) l.El(O.08)
1.00 1.06 1.29
l.OO(O.08) 1.17(0.03) 1.38(0.00)
1 1
1 1
1 1
0.44
1
11
1.11 3.19
0.45(0.01) 1.15(0.03) 3.32(0.15)
1.08 1.20 1.44
1.07(0.10) l.lE(O.04) 1.40(0.11)
0.14 0.61 2.49
0.15(0.01) 0.65(0.03) 2.63(0.15)
1.22 1.35 1.55
1.28(0.02) 1.46(0.04) 1.69(0.05)
fX 0:7
1.5 1.5 1.5
1 1 1
1 1 1
0.49 1.31 4.11
0.48(0.01) 1.25(0.03) 3.41(0.16)
1.47 1.56 1.82
1.35(0.02) 1.39(0.03) 1.54(0.03)
0.19 0.81 3.41
0.19(0.01) 0.76(0.02) 2.72(0.15)
1.41 1.60 1.88
1.46(0.02) 1.64(0.05) l.El(O.03)
0.3 0.5 0.7
0.5 0.5 0.5
1 1 1
1.5 1.5 1.5
0.35 0.91 2.61
0.36(0.00) 0.90(0.00) 2.61(0.11)
0.65 0.89 1.17
0.61(0.01) O.EOfO.01) 1 .i 9io.27j
0.05 0.41 1.91
O.OS(O.00) 0.4310.091 i.92io.i 1j
1.20 1.21 1.34
1.22(0.05) 1.36(0.03) 1.50(0.04)
E.35 0:7
1 1 1
1 1 1
1.5 1.5 1.5
0.44 1.11 3.11
0.44(0.01) 1.16(0.04) 3.26(0.18)
1.10 1.24 1.45
1.01(0.01) 1.22(0.02) 1.43(0.04)
0.14 0.61 2.41
0.15(0.01) 0.67(0.03) 2.58(0.20)
1.29 1.40 1.55
1.35(0.03) 1.55(0.03) 1.74(0.07)
0.3 0.5 0.7
1.5 1.5 1.5
1 1 1
1.5 1.5 1.5
0.48 1.28 3.85
0.46(0.01) 1.26(0.06) 3.51(0.26)
1.49 1.58 1.79
1.35(0.03) 1.43(0.02) 1.57(0.04)
0.18 0.78 3.15
O.lE(O.01) 0.77(0.06) 2.81(0.25)
1.47 1.63 1.85
1.51(0.03) 1.71(0.03) 1.86(0.04)
0.3 0.3 0.7 8:: 0.7
*An approximation outcome. tA simulationoutcome.
Fig. 3. Comparison of approximation and simulation outcomes.
DAVID SONDERMANand BEHNAMPOURBABAI
82
I2 22,1.441
f2.49,1.55)
(2.19,1431
Fig. 4. Three recirculating conveyors.
rate of i,, = 2.49, and an approximated c2. In turn, these values of (L,, cJ superimpose with the primary input stream to form a process with the values (n,, cl). This demonstrates how a steady-state can be reached even when the tragic intensity at the server is greater than one. This is similar to an ordinary queueing system with a finite waiting room, but in this model we have explicitly dealt with the overflows. 5. APPLICATION:
MATERIALS
HANDLING
SYSTEMS
IN SERIES
In this section, a tandem recirculation system is considered. The system consists of three single server recirculating systems. Each recirculating system could represent a production system consisting of a single processor, a loading and an unloading stations linked by a material handling system. The material handling system could be a closed loop conveyor system. The overall behavior of the system is of interest to us in this section. Consider Fig. 4, where we display a schematic diagram of a sequence of 3 stochastic recirculation systems in tandem, as might be found in an automated, flexible materials handling setting. The departure process from the first system becomes the arrival process to the second system, and so forth. It is of interest here to characterize the overall stochastic behavior of such a system, including the distribution of the number of items on the various individual conveyors. We present results for the case of equal exponential service rates at each service subsystem. We first need a result concerning the departure process from a single server queueing system with exponential service times and no waiting room. Proposition
2
Consider an Hz/M/l/1 queueing system. Arrivals form a sequence of i.i.d. random variables with c.d.f. H,(p, fi; t) and the one server is exponential at rate p. Then the first two moments of the interdeparture distribution with laplace stieltjes transform D*(s) are: (23)
. ..(.,,;li+,,,k,(~-~)+2k~k~(~-~)
(24)
k, =
(25)
where --r!F[l + m*(P)] ( P +e, )
k2=H k 3=
-
( eI
2 [l + m*(P)] p )
(26)
(27)
Single server stochastic recirculation systems
k4=
83
(-l-- > 62
P-
m*(P)=
WP, & P) 1 - fm, 6;PO’
(29)
Proof
The proposition is a special case of equation (2.4) of Daley [7]. The results of this proposition are combined with our algorithm to approximate the behavior of the stochastic processes representing units on both legs of each conveyor. In Fig. 4 we have displayed the results of our algorithm for the input parameters i,, = 0.7, c, = 1. Starting with Poisson input, we note that the departure process from the first conveyor has less variability than the input process (c = 0.82). Using this departure process as input to the second conveyor, the same phenomenon occurs to a lesser degree with the second departure process (c = 0.80). At this point we reach an equilibrium state where an input of (0.7, 0.80) yields a departure process with the same two parameters. In other words, the conveyors in series create a situation where the stochastic behavior of the departure processes from the individual conveyors becomes less variable as the units proceed “downstream”, and at some point an equilibrium state is reached. It will also be of interest to examine the overall congestion pattern in the series of conveyors. To do so we need to assign some physical characteristics to the individual systems. We assume that the amount of time required for a unit to transit each leg of each conveyor is 100 time units, although in general these distances may all be different. We then apply Smith’s [8] central limit theorem for renewal processes to calculate the approximate distribution of the number of items to be found on each conveyor. Let N,,(t) = number of items found on leg i of conveyor j, during time period t, i = 1,2, j = 1,2, . . . , t 2 0. Then for sufficiently large t, N,,(t) N Normal&
= i&t, a$ = lijtC2,)
where (n,, cij) are the parameters of the processes along leg i of conveyor j. For example, consider conveyor 1. It can be seen that N,,(lOO)-Normal(~,,
=319,0:,
=25.7’)
N,, (100) - Normal(p, 1 = 249,o: 1 = 24.52). Defining Nj(t)
=
Nij(t)
+
N2j(t)
we obtain N,(lOO) - Normal(pi = 568,~: = 35.5’). In Fig. 4, we have labeled the parameters of the various distributions of the number of items in the systems along with the appropriate 95 ‘A Confidence Intervals (CIs). The lessening of variability as the units move down the series of conveyors is also evident in the distributions of the various N,(t). The complete characterization of the system in this manner enables one to make changes in the design or operation of the system.
Acknowledgement-The
authors would like to acknowledge the financial support of NSF Grant ECS-8105954.
84
DAVIDSONDERMAN and BEHNAMPOURLIABAI REFERENCES
1. D. Sonderman, An analytical model for recirculating conveyors with stochastic input and outputs. I.J.P.R. 20, 591-605 (1982). 2. W. Whitt, Approximating a point process by a renewal process, I: a generate framework. Bell Laboratories, Holmdel, N.J. (1980). 3. W. Whitt, Approximating a point process by a renewal process, I: two basic methods. Opns Res. 30, 125-147 (1982). 4. S. Hallin, Distribution of the interoverflow time for the GI/G/l loss system. Math. Opns Res. 6, 563-569 (1981). 5. C. Palm, Intensity fluctuations in telephone traffic. Ericsson 7&h. 44, l-189 (1943). 6. C. Palm, Research on telephone traffic carried by full availability groups. Tele. No. 1, I-107 (1957). 7. D. J. Daley, Characterizing pure loss GI/G/l queues with renewal output. Proc. Camb. Phil. Sot. 75, 103-107 (1974). 8. W. L. Smith, Renewal theory and its ramifications. JI R. Statist. Sot. Ser. B 20, 243-302 (1958).