Discrete Event Dynamic Systems, 4, 129-169 (1994) 9 1994 Kluwer Academic Publishers, Boston. Manufactured in The Netherlands.

On Using Continuous Flow Lines to Model Discrete Production Lines RAJAN SURI AND BOR-RUEY FU

University of Wisconsin-Madison, Department of lndustrial Engineering, 1513 University Avenue. Madison, W153706 Received October 5, 1992; Revised November 18, 1993.

Abstract. We extend previous research on the use of models of continuous tandem (CT) lines for performance analysis of discrete tandem (DT) production lines. We formalize the translation of input parameters from the DT line to the CT model, as well as the translation of performance measures (PMs) obtained from the CT model back to the DT line. We show that although the CT model conceptually represents a line with continuous fluid, it can be represented as a generalized semi-Markov process (GSMP). This representation leads to a simple and concise simulation algorithm for a CT model. We investigate the accuracy of the CT model for prediction of PMs in the DT line, and show that, with proper translation of parameters and PMs, the CT model provides reasonable estimates for the DT line PMs. We provide preliminary results on gradient estimation for CT models via infinitesimal perturbation analysis. The aim of the paper is to provide a firm foundation for the future exploration of CT models as a means to parameter optimization for DT lines. Keywords: Transfer lines, tandem queues, performance analysis, simulation

1.

1.1.

Introduction

Tandem Production Lines

In the manufacturing world there exist two distinct types of manufacturing systems: discrete and continuous. The two terms indicate whether material moves through the processes as discrete entities (e.g. components and assemblies in an automobile factory) or as continuous fluid (e.g. liquid constituents in a chemicals plant). In this paper we are concerned with a particular configuration of production line, namely a tandem line, either discrete or continuous. Such a configuration consists of a sequence of processing machines in series (Figure 1). The product, whether discrete or continuous, arrives from an external source and starts its processing at the first machine. After being processed, it goes to the second machine, and so on, in order, until it is processed by the last machine, after which it departs from the system. Machines may have different processing rates, and between each pair of machines is a buffer of a specified size. Also, while a machine is processing it may fail, and once failed, take some time to be repaired; both these occurrences are characterized by specified random variables. If a buffer gets full (due to failure, or just slow speed, of a downstream machine) then upstream machines cannot process product at their normal rate. We will discuss these dynamics more precisely later. However, it should be noted that all these characteristics result in complex dynamics of a tandem production line, which makes it difficult to analyze these lines and predict their performance.

130

R. SURI AND B.-R. FU

Machine Buffer Machine

S1

B1

Buffer Machine

82

U n i n t e r r u p t e ~ ~ ~ l ~ Supply of Products

~ ** *

BM-1

SM ~ l

Unlimited Storage for Products

Figure1. Tandemproductionline.

This study is one step in our long term research aim to enhance the set of tools available for modeling and parameter optimization of discrete tandem (DT) lines. These form an important subset of modern production systems: two instances are transfer lines, which are high speed automated lines for manufacturing high volume components such as crankshafts or engine blocks, andfinal assembly lines in a variety of industries. Hence, the analysis of DT lines has received a lot of attention in the literature (see section 1.3).

1.2.

Use of Continuous Tandem Line Models

Our approach in this paper is to focus on the modeling of DT lines using models of continuous tandem (CT) lines. Other researchers have also suggested this avenue (see section 1.3 for a detailed review). The reasons why it is worth exploring this avenue are manifold, and we preview them here. They will be discussed further in the body of the paper. First, there is the possibility of considerable increase in computational efficiency. The analysis and optimization of real world DT lines continues to require large amounts of computational effort, with the modeling tool of choice being discrete event simulation (e.g. see Wei et al. 1989). (The reasons why simulation continues to be the tool of choice are discussed in section 1.4.) As will be seen, simulations of CT line models can often be done much more efficiently. Second, all the parameters in the CT model are continuous (e.g. buffer sizes). From the point of view of optimization, techniques for continuous parameter optimization are much more efficient than those for discrete parameter optimization. In fact, this is even more critical for optimization of stochastic systems, where continuous parameter optimization algorithms have been studied, while little is known about discrete optimization (Glynn 1986). Thus one would hope that the CT model could be optimized with less computational effort than that needed for the DT model. Of course, this still requires understanding how the optimum of the CT model relates to that of the DT line, but this is precisely the type of research result that we are interested in. Third, whenever one has continuous parameters, there is the possibility of using gradient information to speed up the optimization algorithm. Recently there have been many advances in gradient estimation for stochastic systems (Ho 1987, Glynn 1987, Suri 1989), and there is empirical as well as theoretical evidence for the fact that such gradient estimates considerably improve the convergence rate of stochastic optimization algorithms.

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTIONLINES

] 31

Input (Parameter)

\

\ tAJ 9

I%-

~ %

_

_&

_

_

_

F" . . . . . - 1 I~ e c i s i o n su;Dort. t.. _ ~_ _ i , fas~ t~tttler ~lze : I Optimization I

Input

a

onl.

Figure 2. Overviewof researchconcept.

Finally, some previous papers have also considered the use of CT models to analyze DT lines. However, in our review of past work (see section 1.3) we found there did not exist a generally accepted and sufficiently formal basis for this approach--each set of researchers had used different assumptions, and often these were not clearly stated, making it difficult for other researchers to replicate the studies. We feel that a more formal framework, with clearly and explicitly stated assumptions and translations schemes, would assist in better understanding the relationship between the two types of models, and would be beneficial to the research field. All the above ideas are summarized in Figure 2. The figure indicates that production lines are currently analyzed using simulations of DT lines (A). The first step is to study the relation between DT and CT models, and then to obtain methods for optimizing CT models (B). With such tools available, we could then hope to use CT models themselves as decision support tools for production lines (C). Figure 3 outlines the main research steps needed to obtain this goal. As seen then, the first step in our research plan is to thoroughly understand the representation of DT lines by CT models. Only after this step, and after convincing ourselves of its merit, can we proceed with issues of performance analysis and parameter optimization of DT lines using CT models. This issue of representation is the focus of the current paper.

132

R. SURI AND B.-R. FU

Investigate Accuracy of CT Representation

II Gradient E s t i m a t i o n for CT Model

II Optimization Algorithms for CT Model

II T r a n s l a t i o n to DT Line Parameters

i Figure 3. Outline of ongoing research steps.

The contributions of this paper are as follows. (i) We formalize the correspondence between the two types of models--this includes the translation of input parameters from the actual DT line to the CT model, and the translation of performance measures (PMs) obtained from the CT model back to the original DT line. It will be seen that this translation is not trivial, and a few subtle issues need to be addressed carefully. (ii) We clearly state all the assumptions used in both the DT model and the CT model. This is particularly important for the CT model since, as we stated earlier, previous papers have not specified all the details of their CT model simulations (see the discussion in sections 1.3 and 1.4). (iii) Based on the two previous items, we present a simple and efficient algorithm for simulating a CT model. With our approach, the algorithm will be seen to be quite concise and easy to program. We hope that these first three contributions will make it possible for other researchers to completely replicate our empirical work. (iv) We show the possibly surprising result that

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETEPRODUCTIONLINES

133

(with appropriate choice of state and events) a CT model can in fact be represented as a discrete event dynamic system (DEDS)--the surprise arises from the fact that the CT model represents a line with continuous fluid, while the term DEDS conjures up the notion of only modeling discrete entities (e.g. Ho 1987). (v) We investigate the accuracy of the CT model for prediction of PMs for the DT line. The proper translation of parameters and PMs, as described above, is key to the satisfactory results of this step. In fact, earlier studies on this topic may have underestimated the accuracy and potential of CT models due to inadequate translations of parameters or PMs. Finally, overall, we hope that this paper builds on and formalizes the previous '~ork on this topic to establish firmly that the use of CT models can be an efficient and sufficiently accurate alternative for the analysis of DT lines. In doing so, we hope to have completed the first step in our own research plan, and provided the foundation for additional work by ourselves and other researchers working on this topic.

1.3.

Previous Results

Before reviewing the literature it is important to clarify one aspect of our terminology. Throughout this paper we will use the terms "discrete" (or "DT") and "continuous" (or "CT") to describe whether the material in the production line consists of"discrete entities" or "continuous fluid". The terms "discrete" and "continuous" can also be used to describe the way that time is being represented. In the context of production line analysis, discretetime models usually apply to rather special cases. For example, a line with machines having deterministic cycle times and geometric failures/repairs, can be represented by a discrete-time model provided that cycle times at all machines are equal. Similarly, other special cases can be found where discrete-time models would be appropriate. Since we are interested in a more general class of models (which usually require simulation, see below) discrete-time representations will not be considered and so our use of the terms "discrete" and "continuous" will, unless otherwise stated, refer to the material in the line." There is a substantial body of literature on the analysis of tandem production lines. For pre-1978 work, a good bibliography can be found in Buzacott and Hanifin (1978), while Dallery and Gershwin (1992), and Suri et al. (1993), provide more recent references. Even so, analytic results are available only for 2-machine and 3-machine lines. The analysis of longer lines has involved the use of heuristic approximations or simulation models. We begin by reviewing work on analysis of DTlines. Analytic solutions for 2-machine DT lines were given in Buzacott (1967), Buzacott and Hanifin (1978), Gershwin and Berman (1981), Gershwin and Schick (1983), Gershwin (1987). The results there are limited to lines where either the machines have the same cycle times ("homogeneous lines") and geometric failures and repairs (Buzacott 1967, Buzacott and Hanifin 1978, Gershwin and Schick 1983), or exponential service times and exponential failures or repairs--see Gershwin and Berman (1981) for homogeneous lines, and Gershwin (1987) for non-homogeneous lines. These models assume "communication blocking" which may be inappropriate for production lines, where "manufacturing blocking" is the common practice. (For a definition of the types of blocking see Altiok and Stidham 1982.)

134

R. SURI AND B.-R. FU

For longer DT lines, approximate solution algorithms were proposed in Gershwin (1987), Choong and Gershwin (1987), Dallery et al. (1988), and Gershwin (1989). These algorithms solve a set of nonlinear balance equations and boundary equations iteratively, using the 2machine line solution as a "building block", until convergence is achieved. The first of these papers, Gershwin (1987) proposed a decomposition algorithm to approximate long discretetime DT lines. This algorithm sometimes had poor performance, so Dallery et al. (1988) provided an improved algorithm based on the same decomposition. A similar decomposition algorithm for continuous-time DT lines Was presented in Cboong and Gershwin (1987). It had been reported that the algorithm in Choong and Gershwin (1987) fails to converge for some cases. Gershwin (1989) improved their algorithm based on the approach used in Dallery et al. 1988. For all these approximation methods the convergence properties and error bounds remain open questions. In the "real-world" of manufacturing, a large class of DT lines is that where machines have deterministic, but unequal, cycle times, along with failure and repair characteristics. For example, "transfer lines" which are high-volume lines used in the automobile industry, have such characteristics. Since good analytical models are not available for this common and important case, an alternative approach to efficient DT line analysis has been to use simulation, along with gradient estimation and optimization techniques, to obtain optimal parameter values in a more efficient manner than repeated simulations. The first paper on gradient estimation of throughput with respect to buffer sizes in DT lines was by Ho et al. (1979) who also introduced the technique of perturbation analysis (PA) for gradient estimation in DEDS. Ho et al. (1983) presented results on the gradients of throughput with respect to buffer sizes, cycle times, and repair times in DT lines, again using PA methods. Cao and Ho (1987) studied PA gradient algorithms for cyclic DT lines with finite buffers and random service times, and provided some theoretical results on the consistency of the PA estimates. They also used the PA estimates along with stochastic approximation techniques to optimize the throughput of the cyclic DT line. Another set of papers focused on more efficient simulations of DT lines. Kouikoglou and Phillis (1989, 1991) proposed an exact simulation algorithm for DT lines which is more efficient than conventional piece-by-piece simulators. They then applied PA to obtain the gradient of throughput with respect to repair rates, in order to find the optimal distribution of repair rates among machines which maximized the throughput. Note however, that they erroneously cite the results in Cao and Ho (1987) as justification for their approach--the proofs in Cao and Ho (1987) are specifically for cyclic DT lines. (On the other hand, there have subsequently been several papers, e.g., Hu (1992), Glynn (1992), and Robinson (1993)) which indeed provide results on consistency of PA gradients for DT lines.) Finally, in this series of papers, Ramsden et al. (1986) demonstrated that fast simulation rates can be obtained using a hardware simulator for DT lines. We now review work on analysis of CT lines. For CT lines with time-dependent failures (as defined by Buzacott and Hanifin 1978), exact solutions for 2-machine models can be found in Sevast'yanov (1962), Wijngaard (1979), and Yeralan et al. (1986). 2-machine models with operation-dependent failures (see A4 in section 2.2 for a definition) have been solved by Gershwin and Schick (1980) for both homogeneous and non-homogeneous lines, and independently by Dubois and Forestier (1982) for homogeneous lines. An approxi-

ON USING CONTINUOUS FLOW LINES TO MODELDISCRETEPRODUCTIONLINES

135

mation algorithm for longer CT lines based on the 2-machine solution was presented in DaUery et al. (1989). While this algorithm looks promising, it works by transforming any line to a homogeneous line, which can introduce significant errors under some parameter settings. Mitra (1988) analyzed the case ofm machines feeding a single buffer followed by n machines. D'Angelo et al. (1985, 1988) presented a simulation algorithm for CT lines and showed that such simulations could be much more efficient than simulating DT lines. However, they did not explicitly study issues of how the (input) model parameters and the (output) PMs of the two types of models are related. Caramanis (1987) combined simulations of CT lines, the zeroth order PA algorithm for the gradients of the throughput with respect to buffer sizes in DT lines (Ho et al. 1979), and nonlinear programming techniques to solve CT line simulation optimization problems. However, there is no derivation or justification of how the PA algorithm from DT lines (Ho et al. 1979) can be applied to a CT line. (Our work shows that, at the very least, certain algorithmic extensions are needed.) It is also likely that the zeroth order PA algorithm, even if correctly extended for the CT case, would still provide biased gradient estimates. Further, it is well known in the stochastic approximation literature that standard nonlinear programming algorithms will fail in such a situation, unless they are substantially modified for the stochastic nature of the simulations. As a generalization of the above simulation-based research, Kouikoglou and Phillis (1990), and Phillis and Kouikoglou (1991) used a continuous flow simulation model to approximate networks of discrete-material queues where products Can be merged/split before/after service. Kouikoglou and Phillis (1992) extended the work by D'Angelo et al. (1988) and demonstrated that DT lines with machines having exponential service times, geometric failures, and general repair distributions can be approximated by a similar continuous flow simulation model. The exponential service times are approximated in the flow model by adjusting (with a pre-specified frequency q) the maximum production rate for each of the machines during simulation runs. They showed that larger q gives more accurate throughput estimates but the computer run times also increase. Therefore, a trade-offbetween accuracy and efficiency has to be made. Some control-theoretic approaches also use CT models for production networks--however they have a different analysis aim (see section 1.4). David et al. (1990) showed the interesting result that the throughput of 2-machine DT lines with a buffer of capacity B (not including the spaces on the machines) can be bounded by certain CT models with buffer capacities B and B + 1. A comparison of 2-machine DT and CT models can be found in De Koster and Wijngaard (1989) where communication blocking and time-dependent failures were assumed. The relative errors in throughput and average buffer level were shown. Their results indicated that the use of a CT model to approximate a DT line is often justified. Simulation results of CT and DT lines with manufacturing blocking and operation dependent failures can be found in Alvarez et al. (1991). They compared the throughputs of a number of 5-machine lines, with all the machines in a given line having identical reliability parameters. Several empirically-based conclusions were drawn from their study; for example, one of their conclusions is: "the approximation is good (error about 1% or less) as long as the buffer capacity is larger than 10, whatever the other parameters".

136

1.4.

R. SURI AND B.-R. FU

Perspective on Prior Results

Since the prior work on both DT and CT lines is so extensive, it is important for us to clarify our contribution, so we devote a few paragraphs to providing a perspective on the previous research. We begin by raising the issue of the need for clear representation in order to ensure replicability of research efforts. In our search of the literature, we have not found prior work that provides a generally accepted and sufficiently formal basis for the use of CT models to analyze DT lines. Three main subjects need to be addressed for such a basis to be obtained: first, one needs an unambiguous translation scheme for parameters; second one needs a complete description of assumptions needed to simulate the two models; and third, one needs a clear translation of PMs between the two models (not just for throughput, but including more detailed PMs such as work-in-process and percentage of time a machine is blocked). We now elaborate on each of these three subjects. Regarding the translation of model parameters between DT and CT models, there are three issues to be addressed. First, how does one translate cycle times? Although Gershwin and Schick (1980) presented a 8-transformation which provides a link between homogeneous DT and CT models, it cannot be applied to non-homogeneous D T~and CT lines. Second, how does one translate machine reliability parameters? As we will see in section 2, the machines in DT lines always work at "full speed" (when they work) while those in CT lines sometimes work at "reduced speeds (flow rates)". This needs to be accommodated in any reasonable translation. Third, how does one translate buffer sizes? In a DT line there is space for one work piece at each of the machines, while no such capacity is present in a CT line. These points will be discussed later in section 2.2.3 where a suitable translation of the model parameters will be presented. We now discuss reasons for use of simulation. First, accurate analytic methods are not available for the case of lines with deterministic (but different) cycle times, i.e., nonhomogeneous deterministic lines with failures. Yet this is the most common case in realworld transfer lines. So the analytic methods do not work for the most important real-world situations! A second reason is that there are many variations in the real-world configurations, such as, one buffer feeding multiple machines, or re-sequencing of parts at certain stations (this happens in GM Engine assembly lines and has been studied by many researchers). In such situations, given the current state of analytic methods, we have to resort to simulation. So continued investigation of simulation is important. Our study is but one small step in this continuing investigation of the various modeling approaches. Next is the issue of the detailed assumptions needed for simulation. Many modeling issues, such as the blocking mechanisms and failure dynamics, need to be specified explicitly in order that other researchers can duplicate the simulations. These and other assumptions can be found in section 2.2, and we are not aware of previous papers that provide such a comprehensive statement for CT models. Last on the subject of replicability is the translation of PMs. Two main questions arise when we consider this issue. First, how does one translate work-in-process (WIP) in the line? In a DT line, the WIP includes products at the buffers and at the machines, while in a CT line, there is no WIP at a machine. Second, how does one translate machine utilizations?

ON USING CONTINUOUS FLOWLINES TO MODEL DISCRETEPRODUCTIONLINES

137

The interpretation of machine utilization in a CT line is more subtle since a machine can work at a reduced flow rate. We will explain this concept and also answer these questions in detail in sections 2 and 4. Based on these questions we believe that the seemingly obvious relations between the two types of lines are deceptive--there are some subtle issues to be resolved if we are to perform sound research on the use of CT lines to accurately model DT lines. Other researchers could also benefit from having a clearly stated framework. We feel this justifies our further study and detailed clarification of this topic despite the previous work on it. Next is the subject of accuracy of the CT model for DT line analysis. Several papers simply mention this point with a couple of illustrative examples. A few more comprehensive studies are available for short lines. Most of the studies only compare a few main PMs. Specifically, we did not find any work which thoroughly compares many detailed PMs for the cases of longer lines. On the topic of gradient estimation, there are algorithms available for DT lines but none have been presented that are specifically derived for CT lines (see the detailed discussion in section 1.3). Recent control-theoretic approaches to production line optimization also use continuous flow models, but they focus on characterizing the structure of the resulting system dynamics (e.g., Malham6 and Boukas 1989, Sharifnia et al. 1989) or the structure of the optimal control policies (e.g., Akella and Kumar 1986, Sharifnia 1988). None of these methods relate directly to the design optimization problem in DT lines. However, the longer term impact of our (and other researchers) work may indeed provide a tie-in to these approaches, since we are building a foundation that justifies the continuous fluid approximation that is used in those control-theoretic models. In this paper, we will focus on DT lines with constant (non-homogeneous) cycle times, operation-dependent failures, and manufacturing blocking (details are in section 2.2)--as mentioned in section 1.3, this is an important situation in manufacturing and good analytical models are not available, so simulation is the only means for performance analysis. We will formally define the translation of the input parameters as well as output PMs between such DT lines and CT models. We will clearly state all the assumptions used in simulating both types of lines. The accuracy of several PMs that are not compared in De Koster and Wijngaard (1989), and Alvarez et al. (1991) will be investigated. We also extend the comparison to longer lines. Lastly, we provide a detailed gradient estimation algorithm tailored to the CT line, along with some empirical results on its performance.

1.5. Overview of Paper The rest of this paper is organized as follows. Section 2 discusses the assumptions for tandem line models and the translation of parameters between discrete and continuous models. A formal model for CT lines based on a GSMP formulation is presented in section 3. Section 4 describes the translation of the PMs between the two models. The accuracy of the CT model representation for a DT line is evaluated via numerical results in section 5. Section 6 presents preliminary results on gradient estimation for CT lines. Section 7 concludes this paper and discusses the directions for future research.

13 8

2.

R. SUR| AND B.-R. FU

Model Description

A tandem line, whether discrete or continuous, consists of M processing machines ($1, $2, . . . . SM) in series, connected by M - 1 buffers (Bl, B2. . . . . BM-l) of finite sizes (see Figure 1). A tandem line in which the cycle times at all machines are the same is called a homogeneous line. In this paper we will allow the cycle times to be different (also known as a non-homogeneous line). In this section we discuss the modeling assumptions and parameters used to characterize both the DT and CT lines.

2.1. A Motivating Example We will begin by discussing a short sample path for a 2-machine system, for each of the two models. The dynamics of the two lines are quite different, and especially for readers used to dealing with DT models, this example will help in bringing out many of the concepts and assumptions for CT models, which we will then present formally in the next subsection. The two sample paths are shown in Figure 4. In order to illustrate both the similarities and the differences between the two types of models, we have made the sample paths evolve in a "similar" manner. First we discuss the sample path for the DT line, i.e. Figure 4(a). Let the cycle time for S1 be 1 and for $2 be 2, and let the capacity of B1 be 1. The products being made are labelled as P 1, P2, and so on. Assume that, at time 0, both $1 and $2 are working (not failed), with P2 starting a cycle at $1 and P 1 starting a cycle at $2, and that B1 is empty. At time 1, P2 departs from Sl but it must then wait in B1 since $2 is busy, and P3 immediately starts a cycle at $1. At time 2, S1 and $2 both complete their cycles. P 1 leaves $2, P2 moves from the buffer to $2, P3 moves from Sl to the buffer, and P4 starts a cycle at S1. Now, at time 3, a new type of event occurs. P4 completes its cycle at Sl but since B1 is full, P4 cannot leave, and $1 cannot continue to work. In this case, S~ is forcibly idled, and is said to be blocked (we assume "manufacturing blocking" see section 2.2). At time 4, $2 completes its cycle, allowing P2 to depart and P3 to move from the buffer to $2. This allows P4 to move out of $1 into the buffer, and $1 is able to start working on the next product--it is now said to be unblocked. Similarly, the reader can follow the sample path through time 6, when another blocking/unbloCking cycle is completed at $1. Now suppose the first failure of S1 is destined to occur after it has operated for 4.5 time units since the beginning of this sample path. This is called the operating time to failure. Then we see in Figure 4(a) that $1 fails at time 6.5 (since it was not working for 2 time units while it was blocked), while still processing P6. In this condition, $1 is said to be down. Meanwhile, $2 continues to operate, but at time 10 we see that it has no more products to work on, and goes idle. It is said to be starved. Suppose $1 is repaired 4.5 time units after it failed. This is called the time to repair. Then it starts working again at time 11, finishes its work on P6 at time 11.5, at which point P6 goes immediately to $2 which can now start working again (i.e. its starvation is terminated). This sample path has illustrated most of the basic concepts and states ofa DT line. We will now do the same for a CT line, which will introduce some additional concepts. Figure 4(b) shows a sample path for a CT line. As mentioned, we have chosen parameters and operating

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

S1 Occupied I

I

I I// A

VIA

1

9.

.q

,t

.~

a

7

~

9

10 11 12 13 Time

I

1

2

3

4

5

6

7

8

10 11 12 13 Time

Occu edl

pl

I I I I

9

Empty

1

2

Empty B1

)

:" .:~i:~ ' ~E"~%~~i~% I ~ ~ i ~ i 4

139

. t

I

I

r~

1

S2

P2 3

P3

4

5

P4

6

7

8

P5 9

P6 10

I

P7

11 12 13 Tim~" ~

(a) D i s c r e t e T a n d e m L i n e

Vl

~

~

1

B 2"

~ v2 0.5

3

4

5

" ' ~"~ 6

7

8

9

'~ ~

~''"

10 11 12 13 Time

/ j ~k

1

2

/

P1

!

1

2

Blocked

t

V

3

4

5

6

7

8

9

10 11 12 13 Time

,

P2 t, P3 ,t P4 I, P5 . . . . . P5 i, P6 t t 3 4 5 6 7 8 9 10 11 12 13 Time (b) C o n t i n u o u s T a n d e m Line

Down

Starved

Cycle

Reduced Cycle

Figure 4. Sample paths of 2-machinelines.

conditions in a way such that the sample path "resembles" that of the DT line. (This will be formalized later). The processing capacity of $1 is 1 unit of volume per unit time, and that of $2 is 0.5. (Note that this means that $2 produces 1 unit of product in 2 time units, which corresponds to the cycle time of $2 in the DT line.) The capacity of Bl is 1 as before. Let vl and v2 be the actual processing rate of S1 and $2 respectively, at a given point in time. Figure 4(b) shows plots of vl and 02 over time. Also shown is the level of product in the buffer. We will now explain these plots in detail. We start at time 0 with the buffer empty and S1 working at its full rate of 1. As soon as Sl begins producing, due to the nature of the CT line, the product flowing out of S] is immediately available to Sz which begins working

140

R. SURI AND B.-R. FU

at its full rate of 0.5. Since the processing rate of $2 is less than that of St, the buffer level starts to rise at the rate of 0.5 units per unit of time (the difference in the processing rates). At time 2, $1 has produced 2 units of volume, and $2 has produced I. Although the product is continuous in this case, we use dashed vertical lines, and labels P 1, P 2 . . . . . to denote the volumes equivalent to unit product completions (so we can compare the sample path with the DT line). Now we have a new phenomenon occurring. At this point the buffer is full, so $1 is constricted by $2 and the processing rate at $1 is reduced to 0.5 (see the solid line for Vl). Note, however, that unlike in the discrete case where S1 got blocked and completely stopped processing, in this case product still keeps flowing through Sl, Bl, and $2. We wilt say Sl is in a reduced flow condition. While such a phenomenon was not seen in the DT line, there is still a strong correspondence between the two sample paths. Consider time units 2 through 6 for $I in both the lines. We see that the average production rate in both cases is 0.5. The reduced flow in the CT line thus mimics, in a "smoothed out" manner, the blocking in the DT line. We will explain the shaded area labelled "B" later. Next, for the CT line we suppose that the failure of Sl occurs after it produces 4.5 units of volume. We will call this the operating volume to failure. Because of the periods of reduced flow rate of S1, the reader can verify that this failure will occur at time 7. $2 can still keep processing since there is product in the buffer, and the buffer level decreases at a rate of 0.5. At time 9 the buffer becomes empty and $2 is starved. As before, we suppose the repair of $I occurs 4.5 time units after the failure, so in this case SI is repaired at time 11.5. Note that $2 also resumes its operation instantaneously. We have now illustrated enough of the sample path to explain all our basic concepts below. Although the sample path for the CT line only introduced reduced flow rate due to downstream constriction, it is also possible for a machine to have a reduced flow rate due to upstream throttling. For instance, let us reverse the capacities, so that of $1 is 0.5 and $2 is 1. Then, at the very start of the sample path, SI will start processing at its maximum rate of 0.5, and since the buffer is empty, $2 will be forced to process at the same rate' of 0.5, rather than its full rate of 1. $2 is now in a reduced flow rate condition, due to upstream throttling. In this case, the reduced flow rate mimics, in an average sense, the idle time that would occur at $2 due to starvation in the corresponding DT line. These observations will form the core of our method, described later, for translating performance measures from the CT model to the DT line.

2.2. Definitionsand Assumptions Now we formalize our models for the two types of lines. The following assumptions apply to both types of models. AI. There is an unlimited supply of material available to $1, that is, Sl is never starved. A2. There is an unlimited storage area following SM, that is, SM is never blocked. A3. There is no transfer delay from machines to buffers, within the buffers, or from buffers to machines.

ON USING CONTINUOUSFLOW LINES TO MODEL DISCRETEPRODUCTIONLINES

141

A4. A machine can fail only when it is operational. This is called the assumption of operation-dependent failure and it is considered the most common type of failure in production lines (Buzacott and Hanifin 1978). A5. The repair time for each machine is exponentially distributed. (For the purpose of our simulation model, a general distribution can be used here. However, we choose the exponential distribution in order to compare our results with analytic ones, where available.)

2.2.1.

The DT Model

Now we describe the model characteristics that are specific for the DT line. We have the following additional assumptions. A6D. The cycle times of machines are deterministic. A7D. The operating time to failure for each machine is exponentially distributed. (The remark in A5, about using a general distribution, applies here too.) A8D. Machines operate and are blocked via the "manufacturing blocking" procedure (AItiok and Stidham 1982, Suri and Diehl 1986). This states that if a buffer is full, the preceding machine may begin work on the next piece, but if it finishes its cycle and the buffer is still full, it will then be blocked. With the above assumptions, our DT line model can be completely characterized by the following parameters: T/ fi

ri B(j)

Cycle time of S i Mean operating time between failures for Si Repair rate for Si Buffer size of Bj (non-negative integer)

w h e r e i = 1,2 . . . . . M, a n d j = l , 2 . . . . . M - I . In the literature on DT lines, one usually finds an alternative parameter for characterizing machine failures, namely Pi, the probability of failure during a cycle at Si. In this case, failures are assumed to be geometric. The relation between such a model and ours is easily seen. In the geometric failure model, the mean number of cycles between failures at Si is 1/pi. Thus the mean operating time between failures for Si is Ti/Pi. Hence, fi = Ti/Pi, or Pi = Ti/fi. Thus we can convert parameters from other models in the literature, into our framework, for the purposes of comparison of results.

2.2.2.

The CT Model

For the CT model we have the following additional assumptions. A6C. The processing rate at a machine can range from 0 to its maximum processing capacity.

142

R. SURI AND B.-R. FU

A7C. The operating volume to failure for each machine is exponentially distributed. (The remark in A5, about using a general distribution, applies here too.) It is useful to review the reason for using the operating volume to failure for machines in the CT line, when we used the operating time to failure for the DT line. The motivating example we gave in section 2.1 is important to understanding this. The basic idea is that we assume the machine fails after producing a given quantity of product. In the continuous line, machines can work at rates below their capacity, so it requires us to measure the volume produced, instead of just their operating time, to correctly determine the next failure. Note also that for the CT model there is no need for an analog of assumption A8D which specifies the particular blocking definition used in the DT model. The different blocking mechanisms in DT models arise from the space available for a work piece at a machine. Specifically, the communication blocking mechanism prevents products from moving into a machine if the succeeding buffer is full, while the manufacturing blocking mechanism prevents the products from moving out of a machine if the buffer is full. Since the difference arises from the way of using the space at the machines, and there is no such space in a CT line, the communication and manufacturing blocking mechanisms reduce to the same behavior in CT lines. With the above assumptions, our CT line model can be completely characterized by the following parameters:

Ci wi ri B(j)

Maximum flow rate (processing capacity) of Si Mean operating volume to failure for Si Repair rate for Si Buffer size of By (non-negative real number)

w h e r e / = 1,2 . . . . . M, and j = 1,2 . . . . . M -

2.2.3.

l.

Translation of Parameters Between Discrete and Continuous Models

We now discuss our proposal for constructing an "equivalent" CT model from DT line parameters. We use one unit of volume of product in the CT line to represent one item produced in the DT line. In this case, the obvious translation of machine capacities is Ci = 1/T/. The repair rates translate directly. In converting buffer sizes, one faces the problem that in the DT line, there is space for one workpiece in each machine. No such capacity is present in the CT line. Should this be accounted for by adding to the buffer sizes for the CT line? For now, we suggest simply using the given DT line buffer sizes. However, we will argue this point later in the paper. For the translation of the failure parameters, note that in the DT line the quantity produced when Si works for fi time units is fi/Ti. Thus we suggest wi = fi / Ti. In summary, we have: 1

c, = ~ ,

(1)

w, = ~.,

(2)

ri and B(j) remain the same.

(3)

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

3.

143

A Formal Model for a CT Line

Now we provide a formal model for the dynamics of a CT line. For DT lines, such a model is commonly available, for example in standard textbooks on discrete event simulation (Law and Kelton 1982, and Banks and Carson 1984), or in publications on DEDS (Suri 1987, 1989, and Cao and Ho 1990). Thus we will not present a DT line model here. However, for CT lines we feel it serves an important purpose to present a clearly stated model. We elaborate on this point in section 3.1. We emphasize that we are interested in defining a CT model that gives reasonable performance estimates for DT lines, rather than attempting to make the CT model behave in detail like a DT line. While elaborate schemes for doing the latter are available (for example, see Kouikoglou and Phillis 1989, 1991, and Phillis and Kouikoglou 1991) they make the simulation algorithm complex, the translation of PMs more complicated, and the implementation of PA algorithms even more difficult. Our aim is to keep everything as simple as possible while preserving a reasonable level of accuracy. We will return to this point again in section 7.

3.1.

Motivation for GSMP Representation

A generalized semi-Markov process (GSMP) is a mathematical framework for studying DEDS. While it will be useful if the reader is familiar with the basics of GSMPs, that will not be necessary in order to understand the main points of our model and the simulation algorithm that follows later. Our explanation is self-contained. Readers interested in learning more about GSMPs are referred to Glynn (1989) for an introduction and bibliography. In this section we show that it is possible to construct a GSMP representation for the dynamics of a CT line. This result may be surprising since one associates GSMP models with DEDS which typically have discrete entities (e.g. customers or parts), and as stated by Glynn (1989), "DEDS are frequently used as models of systems having piecewise constant trajectories". Thus, one does not expect a model with continuous fluid in it, and with a timevarying trajectory as in Figure 4(b), to be representable as a GSMP. To clarify this further, we are not just showing a loose relation between CT lines and D E D S - - w e are showing that a suitable DEDS can be constructed which provides an exact representation for a CT line. Once we show how this is done, it will be seen to be quite obvious. Making this formal connection, even if not very deep, does have some benefits. First, it provides a formal model for a CT line that might serve as a standard reference for future researchers. Such a standard seems to be lacking in the CT line literature (sections 1.3 and 1.4). Second, the resulting simulation algorithm will be found to be both simple and concise. Third, many results are available in the context of GSMPs and these could be brought to bear on CT models where needed. The papers by Glynn (1989) and Whitt (1980) provide instances of such results. A concrete example of such results that might prove useful in our context is the work by Glasserman (1991). This provides several sets of results on conditions that ensure the consistency of infinitesimal perturbation analysis (IPA) for gradient estimation in GSMPs. Similarly, Fu and Hu (1992) present a general method for deriving smoothed perturbation analysis (SPA) gradient estimators when the system being studied is represented as a GSMP.

144

R. SURI AND B.-R. FU

Table I. Summary of physical state variables. Variable

Description

cei(t)

State of Si at time t Flow rate of Si at time t Level of Bj at time t

vi (t) xj(t)

Table 2. Summary of events. Symbol

Event Represented

5r/ 7~i

Failure of Si Repair of Si By becomes Empty By becomes Full Termination of the simulation

t3Cj ~.~j "Tf

It is possible that these results, or potential extensions of them, could be brought to bear on CT models after our formalism is established.

3.2.

A GSMP Representation for a CT Line

Now we will describe the dynamics of a CT line in a form that enables GSMP representation. We begin by defining the variables that correspond to the physical state of the line. Let t denote the time along the sample path, and o~i (t) be the state of Si where

~i(t) =

D, O, S, B,

if if if if

Si Si Si Si

is is is is

down (failed), operational (at its full capacity of Ci), starved (reduced flow or zero flow), blocked (reduced flow or zero flow).

Let vi(t) be the flow rate of Si at time t. Let xj(t) be the level of Bj at time t. For ease of reference, the set of physical state variables is summarized in Table 1. Next we consider the events that occur in a CT line. (We use the term events in the formal GSMP sense of the term; however, readers not familiar with GSMPs can assume the usual DEDS interpretation of this term.)-The set of event types is E = {~., ~i, 13Cj, 13.T'j, Ty: i = I, 2 . . . . . M, j = 1,2 . . . . . M - 1}, and each of these events is defined in Table 2. In our model we will assume that the sample path is terminated when the CT line achieves a given production target, Q, (i.e. when the last machine has output Q units of product). In a GSMP, each event is associated with a clock representing the residual lifetime of that event, and each clock has a rate at which it runs down (called the clock speed). When the clock associated with an event runs down to 0, that event occurs. Upon the occurrence of an event, changes may occur to the physical state, clock settings, and clock speeds, according to defined rules. Now we will describe the clocks suitable for use in the CT line model.

ON USING CONTINUOUSFLOW LINES TO MODEL DISCRETEPRODUCTIONLINES

145

In the following, we will use the fact that whenever a clock speed is set to 0, the clock value is irrelevant. Hence whenever we set a clock speed to 0, we will not need to specify the setting of the clock value. First we discuss the clocks for the failure and repair events. As discussed in A4 and A7C, after a machine is repaired, its next failure will occur after it produces a given volume of product. This operating volume to failure is determined by sampling from a probability distribution, described in more detail later. With this model in mind, whenever Si is operational, let Wi (t) be the remaining volume to be produced by Si until its next failure. Then Wi (t) decreases at the processing rate of Si, namely vi(t), and the failure of Si occurs when Wi (t) is 0. Therefore, when Si is operational, Wi (t) is a suitable clock for the event f / and its speed is vi (t). When Si is failed, the speed of this clock is set to 0. Next, whenever Si is failed, let Ui (t) be the remaining time to the repair of Si. Then Ui (t) decreases at the same rate as actual time, that is, a rate of 1. The event 7~i occurs when Ui (t) is 0. Hence, when Si is failed, Ui (t) is an appropriate clock for 7~i, with speed l. When Si is operational, the speed of this clock is set to 0. Now we consider the filling and emptying of buffers. The level of buffer Bj decreases if Sj+ 1 works faster than Sj, and the rate at which it decreases is vj+l (t) - vj (t). The buffer will become empty when this level reaches 0. Thus, whenever Vj+l (t) > Vy(t), a suitable clock for event 13,fj is the buffer level, xj (t), with speed vj+t (t) - vj (t), and its speed is 0 otherwise. By a similar argument, whenever vj (t) > Vj+l (t) a suitable clock for 1337j is B(j) - xj (t) (i.e., the excess capacity of By at time t) with speed Vy(t) - Vj+l (t), and with speed 0 otherwise. (Note that if vj+l (t) = vj (t) then both clocks will have speed 0 which is appropriate, since the buffer level remains constant in this case, and so neither the buffer full nor buffer empty event can be scheduled.) Finally, let q(t) be the remaining amount of product to be output by SM in order to achieve the production target. The rate at which q(t) decreases is VM(t). Clearly q(t) is a suitable clock for the termination event Tf and its speed is vM(t). This defines all the clocks needed for our GSMP representation of the CT model. In order to provide a quick reference for the reader, we summarize the clock definitions and speeds in Table 3. Note that the units for the clock readings are different; for example, Wi (t) is in volume units and Ui(t) is in time units. The choice of units for the clock readings makes our representation and resulting simulation algorithm easy to comprehend because all the clocks have phYSical interpretations. (The speed parameter essentially converts the various clock units to standard time units.) Now that we have defined the variables that characterize the physical system states and the clocks for the model, we proceed to describe how the system evolves. Let

k(e, t) r(e, t) E(t)

be the reading of the clock for event e at time t, be the speed of the clock for event e at time t (k(e, t)), be the set of events for which r (e, t) # 0.

E(t) can be thought of as the current set of possible events. At time t, the additional time to the possible occurrence of event e is defined as At(e, t) --

k(e, t) r(e, t)

for e 6 E(t).

(4)

146

R. SURI A N D B.-R. FU

Table 3. Clock definitions and clock speeds. Symbol

Wi (t)

u~(t) xj(t) B(j) --xj(t) q(t)

Clock Description Remaining operating volume to failure for Si at time t Remaining repair time for Si at time t Level of By at time t Excess capacity of By at time t Remaining volume to be produced by SM at time t

Clock Speed t

vl (t)

l{ul (t)=D} [vj+l(t) - vj(t)] + vj(t) - vj+l (t)] + VM(t)

tNote: 1(.) is the indicator function of set (.), [x] + = max(0, x).

The reason that we say "possible occurrence" is that if another event, say e', occurs before e, it is possible that the clock speed or clock value would be reset for e, and the value in (4) would no longer be relevant. With the above definition, the next event to occur is given by

e*(t) = argmin{At(e, t) : e 6 E(t)}. e

If more than one e satisfies the "argmin" on the RHS then, in principle, any scheme can be used to resolve ambiguity. In this case, in order to provide a canonical model, we suggest using lexicographic ordering by event name (i.e. the names in Table 2). In the GSMP model, all clocks continue to run down at their stated speeds until an event occurs. Now let us consider the effect of an event, say e, on the system. Suppose this event occurs at time t. Let the argument t - used with any quantity denote its value the instant before e occurs, and let the argument t denote values after the event has affected the physical states and clocks of the system. Before we consider the effect on the whole system, we start by understanding the "local" dynamics of a CT line, i.e. how buffers interact with their immediate upstream and downstream machines, and how machines interact with their immediate upstream and downstream machines via the buffers. First we will consider the buffer Bi and the machines that are immediately upstream and downstream of it (i.e. the configuration of Si, Bi, and Si+l). We will discuss how the states at Si and Si+l, are affected by the "local" events Bori, and BEi, which are the events related to B i. Then we will consider the machine Si and the buffers that are immediately upstream and downstream of it (i.e. the configuration of Bi-1, S/, and Bi), and the effect o n S i of the "local" events 0rl, and 7~i, which are the events related to Si. (The effect of the above events on the rest of the line, and the effect on Si of events from the rest of the line, is referred to as "global" dynamics and will be discussed later.) We now discuss the buffer events. If Bi becomes full at time t, then v i ( t - ) must have been greater than Vi+l ( t - ) . So Si is now blocked by Si+l. In this case vi+l (t) remains the same as vi+l ( t - ) and we set vi (t) = vi+l (t). Similarly, if Bi becomes empty at time t, then Si + j is now starved by Si . So vi ( t ) = vi ( t - ) and we set vi + l ( t ) = vi ( t ). We can summarize this as follows:

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

147

Table 4. Flow rate update table for machine Si.

If Xi l(t) = = > >

0 0 0 0

then Xi(t) = < = <

B(i) B(i) B(i) B(i)

Vi(t) min(vi-l(t), Q , vi+l(t)) min(vi_l(t), Ci) rain(G, Vi+l(t)) Ci

If event at time t is J~{,~i then {Vi(t) = Vi+l(t);

~i(t) = B;}

If event at time t is BEi then { V i + l ( t ) = V i ( t ) ; oti+l(t) = S; }

For notational convenience, in the logic below, we define pseudo-buffers B0 and BM with xo(t) = o0 (for all t) and B(M) = o~ (in words, we create artificial buffers so that the first machine is never starved and the last machine is never blocked). Now consider the machine events. If Si fails, we set its flow rate to zero, and start its repair clock. (If the failure of Si causes blocking or starvation of its neighboring machines, this will be considered in the global dynamics later.) Thus: If event at time t is .T/then

{v i(t) = O; ~i(t) = D; k ( ~ i , t) = sample(1/ri); } where sample(#) denotes a sample taken from an exponential distribution with mean/x. If Si is repaired, we start its failure clock. In updating the flow rate of Si however, the situation is more complicated, since it depends on the states of its neighboring buffers. The reader can verify that the following logic holds: If event at time t is Ri then {

k(~i, t) = sample(wi); set V i ( t ) according to Table 4; S, ifxi_l(t) = Oand vi(t) = Vi-l(t), oti(t) = B, if xi(t) = B(i) and vi(t) = O i + l ( t ) , O, otherwise;

(5)

Finally we discuss the "global" dynamics, that is, how an event at Si or Bt affects buffers and machines in the rest of the line. We define a "downstream chain" of empty buffers as a set of consecutive buffers following Si such that each is empty, and the buffer following

148

R. SURI AND B.-R. FU

this set is non-empty. Similarly, an upstream chain of full buffers is a set of consecutive buffers preceding Si such that each is full, and the buffer preceding this set is not full. If, as a result of the local dynamics described above, the flow rate at any machine (say Si) is changed, then we check to see if there is a downstream chain of empty buffers. If there is such a chain (say, Bi . . . . . Bi+n) then the change of vi will affect Si+l . . . . . Si+n+l as shown below. We also check to see if there is an upstream chain of full buffers, and if there is such a chain (say, Bi_ 1. . . . . Bi_n) then the change of vi will affect S i - l . . . . . Si-n as shown below. If oi(f) is changed at time t then do { m=i;

while (Xm(t) = 0) and (m < M - 1) do { vm+l (t) = min(vm (t), C,~+1); if vm+t (t) = Vm(t) then Otm+I(t) = S; else Otm+l(t) = O; m=m+l;

} re=i-l; while (Xm(t) = B ( m ) ) and (m > 1) do { v,,(t) = min(C,,, vm+, (t)); if Vm (t) = V,,+I (t) then ~,~+l (t) = B; else otto(t) = O; m=m-1;

} }

/* if */

Any clock, clock speed, or other physical variable not updated in the above equations simply retains its value prior to the event. After this event, the system proceeds along its sample path until the next event occurs, and the process repeats itself until the termination event. This completes the formal description of the dynamics of our CT model. A complete simulation algorithm implementing the above logic is given in Appendix 1. It can be seen that with the approach taken here the algorithm is simple and concise.

4.

Translation of Performance Measures

As mentioned in the introduction, our primary motivation for this work is to use CT models to analyze DT lines. In order to do that effectively, we need to know how to use the performance measures (PMs) obtained from CT models to estimate the PMs for DT lines. In this paper, we will assume that the PMs to be estimated for the DT line are steady state values of: line throughput, average buffer levels, machine utilizations, average product flow time, and average work-in-process (WlP). We will also assume that these steady state values will be estimated via finite horizon experiments, and standard experimental techniques will be used to minimize the effects of initial conditions (e.g. see Law and Kelton 1982, Banks and Carson 1984). Thus below we will define PMs for a finite observation interval, with

ON USING CONTINUOUSFLOW LINES TO MODEL DISCRETEPRODUCTIONLINES

149

the understanding that these observations will be used with a suitable estimation scheme. For instance, a standard scheme would be to run the simulation for a "warm up" period, then observe the PMs for a specified additional "run length". There is extensive literature on determining suitable lengths for both these periods, so we will not discuss that here. The definitions of the above mentioned PMs for DT lines are well accepted, and will be used below. On the other hand, for CT lines it will be seen that the definition of certain PMs (e.g. those relating to blocking and starvation) are not straightforward. While several definitions could be justified, in each case we propose one that (hopefully) will most closely track the DT line PM. The merit of our definitions is eventually dependent on the numerical results on our resulting estimates. Let to, q . . . . be the times of the 0th, 1st . . . . occurrences of events in a sample path. Let the observation period for the PMs be the interval from tm to tn (where tm < tn). Define Ati

=

ti-I

ti - -

and

T = t~ - tin.

Note that the machine states do not change in [tn-l, t,~) and hence the flow rate of machines is piecewise constant. We will use this fact below. We now define the line PMs. In the following, we will use the argument "DT" with a PM to denote the definition for a DT line. PMs without this argument denote our suggested definitions of PMs derived from the CT line model. Throughput: For a DT line, the throughput is defined by TPDT =

number of pieces produced by S M during [tm, tn] T

For the CT line, we define throughput in the natural way, namely, TP = ~ .

vM(t)dt.

Since VM is piecewise constant, we can write this in a form easier for simulation computation,

TP = -~ 9

VM(ti-I) " Ati. i=m+l

Average buffer level: In a DT line, if xj (t) represents the number of parts in Bj, then the average buffer level is defined by Xj.or = "~

xj(t) dt.

(6)

For the CT line, we also define the average buffer level of Bj in the natural way, namely by the same equation above. However, some further observations here make it simple to calculate this integral in a GSMP-based simulation. Observe that if we start with an empty

150

R. SURI AND B.-R. FU

system, then the level of Bj at time t is the difference in the total volume produced by Sj and Sj+l up to time t. So we can write

f0 t [vj(s)

xj(t) =

- vj+l(s)]ds + xj(O),

(7)

where xj (0) is the initial buffer level of Bj at time 0. Since vj (t) is piecewise constant, from (7) we see that xj (t) must be piecewise linear and continuous. Note that the area under a linear segment in a graph can be obtained by averaging the heights of the endpoints. Thus we can replace the integral in (6) by the summation

l

~

[Xj(ti-l) -1-xj(ti)] 9Ati.

t=m+l

This relation will be useful in calculating the value of Sj during our simulation. Machine utilizations: In a DT line, the utilization of a machine is broken out into four categories: cycle (c), starved (s), blocked (b), and down (d). Each of these measures represent the proportion of time the machine spends in the corresponding state. Thus for the DT line, the blocked utilization would be defined as

Ub.or(i) = -~ Jtm lt~'~t)=BIdt. For the CT line, however, the situation is more complicated, since machines can be in a reduced flow condition (see the sample path in Figure 4 of section 2). Although we could define a new set of PMs for this situation, our aim is to estimate the PMs for the DT line. Therefore we propose certain definitions here. We begin with the cycle utilization. For a machine in the DT line, the cycle utilization represents the proportion of time the machine is busy. This is not the quantity we should observe for the CT line, since a machine can be busy but working at a highly reduced rate. In that case, we would intuitively feel that the machine is "not utilized well". To get to a suitable definition for the CT line, consider again the DT line definition which is:

Uc.or(i) -

time in cycle total time

Dividing numerator and denominator by the cycle time we get

Uc.or(i) =

(time in cycle/T/) (total time/T/)

Now the numerator is the number of pieces actually produced in the observation period, while the denominator represents the potential number of pieces that could have been produced if the machine could be in cycle during the entire observation period. Thus we propose the same definition for the CT line, namely, the ratio of the actual volume produced

ON USING CONTINUOUSFLOW LINES TO MODEL DISCRETE PRODUCTIONLINES

151

during T to the total volume thatcould be produced if it worked at full capacity for T. Thus,

Uc(i)

--

vi(t) dt

CiT cir

vi(tk-~)" Ark, k=m+l

where we used the piecewise constant property to derive the last equality. In a DT line, the starved utilization, Us.or(i), is the proportion of time the machine is starved. A similar argument to the one just used for cycle utilization shows that this equals the proportion of the potential capacity of Si not being used for production because of starvation. This motivates the following definition for a CT line:

Us(i)

= -

's

CiT CJ

[Ci - vi(t)l . l/~,(t)=sldt

~

[Ci - vi(tk-1)] 9 l~,(t~_l)=Sj 9 Ark.

k=m+l

To understand this, note that Ci T is Si's potential capacity, and [ G - vi (t)] is its unused capacity at time t. Recall that Si might also work at a reduced flow rate when it is blocked, so the indicator function ensures the integral takes into account only the times when Si is starved. Analogously, we propose the definition of the block utilization, Ub(i), as

Ub(i)

--

cir

[Ci - vi(t)l . l~,~t)=Bldt

CiT1 k='-~,+'~l[Ci- vi(tk-1)l 9 ll~,(t~_,)=Bj 9 Ark.

The only difference between this and the starved utilization is the indicator function which accounts for reduced flow during blocked times only. For DT lines, the down utilization, Ud.or(i), is the fraction of time that Si is under repair. This definition also makes sense for a CT line, so that

1 ['~ 1 ~ Ud(i) = T Jt~ l{c~(t)=D}dt = -~ y ~

l{~,(t~_,)=D} 9 Ark.

k=m+l

Average WIP: This is the average quantity of products in the whole line. In a DT line, this includes products at the buffers and at the machines. In our CT model, however, there is no WlP at a machine. Thus we need to account for this difference. Several possibilities exist, including modifying the buffer size in the CT model (e.g. make all the CT model buffer sizes larger by one unit compared to the DT line). We will discuss this further in section 7. For now we propose keeping the same buffer sizes, but adding an extra term to the CT line WIP as explained next. Observe that in a DT line the average WlP in a machine Si is given by 1 - Us(i). Since we already have an estimate for the DT line value of U~(i) from our CT model, we can use

152

R. SURI AND B.-R. FU

this to get an estimate of the average WIP for the DT line by adding a term to the CT line WlP, giving M-I

M

WIP = ~ Yj + )--~[1 - Us(i)]. j=l i=l Again, the merit of this proposed definition will be tested via numerical results. Average Flow Time: Since the material in a CT model is continuous fluid the definition and calculation of product flow time would require complex integration to account for the time spent in the system by each infinitesimal piece of material. Rather than compute this, we suggest using an indirect estimate of the average flow time via Little's law. Such indirect estimates have been proposed even for discrete systems (Glynn and Whitt 1989). Therefore, we define a pseudo flow time as Flow Time --

WIP TP

This completes our set of proposed definitions of PMs for the CT line.

5.

Numerical Results

Now we present numerical results for simulations of CT lines. The objectives of this section are twofold: (1) to attempt to verify the correctness of the simulation algorithm, and (2) to investigate the accuracy of the CT model for prediction of PMs for a DT line. Section 5.1 verifies the simulation algorithm for 2-machine and 6-machine CT lines, while section 5.2 compares the PMs estimated from simulations of DT lines and CT lines. An extensive study for 15-machine lines is given in section 5.3. Our simulation program for CT lines is an implementation of the algorithm of Appendix 1 in the C language. The PMs shown in the tables below are obtained from the average of 10 independent replications. For each replication, the simulation is run for a production quantity of 10,000 (warm up) to eliminate the effects of the initial transient, and then statistics are collected for a production quantity of 60,000.

5.1.

Verification of Simulations for CT Models

Here we compare the output of our simulation algorithm for CT lines with previously published numerical results. Four cases of CT lines are considered. Cases 1 and 2 are taken from Gershwin and Schick (1980). Case 1 is a 2-machine homogeneous line and Case 2 is a 2-machine non-homogeneous line. Case 3 and 4 are 6-machine non-homogeneous lines from Dallery et al. (1989). The input parameters for the CT lines are listed in Table 5. The last column of Table 5 lists the isolated throughput (TPi) of Si. This is the throughput that this machine could achieve if it were working in a stand alone mode with an unlimited

ON USING C O N T I N U O U S F L O W LINES TO M O D E L DISCRETE PRODUCTION LINES

153

Table 5. Input parameters for cases 1-4.

Case

Machine i

Flow Mean Oper~ing Rep~r R~e Volume to R~e Ci Fmlure wi ri

Buffer Capacity B(j)

Isolated Throughput TP i

1

I 2

1 1

10.000 1.429

0.7 0.1

4

0.875 0.125

2

1 2

3 5

2.500 2.976

0.40 0.32

5

0.750 0.800

I 2 3 4 5 6

2.809 3.571 3.571 3.571 3.571 2.882

25.955 160.714 160.714 160,714 160.714 26.628

1.316 0.200 0.200 0.200 0.200 1.316

4 2 2 2 4

2.596 3214 3.214 3214 3.214 2.663

1 2 3 4 5 6

2.857 4.000 3.333 3.125 3.333 3.226

85.714 56.000 116.667 125.000 120.000 45.161

0.250 0.154 0.100 0.118 0.083 0.286

100 100 150 250 250

2.521 2.733 2.592 2.579 2.498 2.581

supply of material and an unlimited buffer for finished parts. It is given by (Gershwin and Berman 1981) toi

TPi = [OillC i q- l/l"ri"

(8)

The purpose of listing the isolated throughput is to give an indication of how well the average capacities of the machines are matched, i.e., how "balanced" the line is. Case 1 is a fairly unbalanced line since the isolated throughput of $1 is much greater than that of $2. In Cases 2, 3 and 4 the lines are more closely balanced. Table 6 shows the PMs for the 2-machine models (Cases 1 and 2). Specifically, the line throughput, average buffer level, WIP, probability of blocking at Sl, and probability of starvation of $2 are displayed. The column labelled "CTL simulation" shows the PMs obtained from our simulation, along with 95% confidence intervals. The column labelled "CTL formula" shows the PMs calculated by the exact analysis given in Gershwin and Schick (1980). The table shows that for all the PMs compared, the confidence intervals contain the analytical values. The simulation outputs for Case 3 and 4 are shown in Table 7, and compared with simulation results from Dallery et al. (1989). Here we only compare estimates of throughput and average buffer levels since these are the results provided in Dallery et al. (1989). It can be seen that the two sets of results are in close agreement. Tables 6 and 7 together give us reasonable confidence that our CT line simulation algorithm is correctly written.

154

R. SURI A N D B.-R. FU

Table 6. Comparison of PMs for 2-machine CT lines simulation and analytical results (Case 1 and 2). CTL Simulation

CTL Formula

Throughput Case 1 Case 2

0.125 4- 0.00l 0.561 4- 0.003

0.125 0.561

Average Buffer Case 1 Case 2

3.970 4- 0.001 2.340 4- 0.01 !

3.971 2.334

Case 1 Case 2

5.970 4- 0.001 4.042 4- 0.014

5.970 4.036

Prob. of Blocking Case 1 Case 2

0.857 4- 0.001 0.253 4- 0.002

0.857 0.252

Prob. of Starvation Case 1 Case 2

0 0.298 4- 0.002

0 0.298

WIP

Table 7. Comparison of PMs for 6-machine CT line simulation and approximation results (Case 3 and 4). CTL Simulation

Dallery et al. Simulation

Throughput Case3 Case 4

2.1104-0.008 2.331 4- 0.008

2.116 2.3

Average Buffer Case 3 B1 B2 B3 B4 B5

1.41 0.89 0.87 0.86 2.20

1.4 0.9 0.88 0.88 2.26

Case 4 B1 B2 B3 B4 B5

5.2.

q- 0.02 + 0.02 q- 0.02 -4- 0.02 • 0.03

52 4- 2 54 4- 2 63 4- 6 95 4- 12 59 4- 4

55 54 61 89 52

Comparison of D T Line and CT Line Representations

I n o r d e r to c o m p a r e t h e C T l i n e r e p r e s e n t a t i o n w i t h a D T line, w e a l s o w r o t e a s i m u l a t i o n o f a D T l i n e . S i n c e t h i s is j u s t a c o n v e n t i o n a l d i s c r e t e e v e n t s i m u l a t i o n o f a t a n d e m p r o d u c t i o n

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTIONLINES

155

line, discussed in many textbooks on simulation, we do not describe it in detail. The simulation program for DT lines is written in the C language. While special discrete event simulation languages are available, we chose to use C in order to have a fair comparison of the speed of execution of the DT line simulation with the CT line simulation (also written in C). Using a simulation language would have made the comparison unfair to the DT line simulation. For example, using the SIMAN language (Pegden et al. 1990) the computer run time for the DT line in Case 3 is more than 13 hours on a PSI2 Model 80 (with a 80387 math co-processor), while our C program for the same line takes only 41 minutes. We will compare the execution times of our DT and CT line simulations later. First we compare the PMs for the two representations. Again, statistics are based on l0 independent replications with warm up and production quantities the same as for the CT line simulation. For Cases 1-4, we translate the parameters to their DT line equivalents (see section 2) and then simulate the resulting DT lines. Table 8 compares the main PMs from simulations of DT lines and CT lines for all four cases, along with 95% confidence intervals. The last column shows the percentage deviation of the average CT line simulation estimate from the average DT line simulation estimate, i.e.,

PMcT - PMDT • 100%. PMDT We see that the throughput, WIP, and the average buffer level estimates are all within 3%. Tables 9, 10, and 11 show the details of machine utilizations for Cases 1-4. As defined in section 4, the machine utilizations are classified into blocked, cycle, down, and starved. For Case 1 (Table 9), all the machine utilizations of the CT line agree with those of the DT line to within 10 -3. The estimates of starved utilization of $2 are not significantly different from zero (< 2 • 10 -4) in both DT line and CT line, so an error term is not calculated here. From Tables 9, 10, and 11 we see that for Cases 2, 3, and 4 the relative errors of the detailed machine utilizations for blocked, cycle and down are within 11%. The errors for starved utilizations are higher. We will now discuss this further. Of course, there are fundamental differences in the operation of CT and DT lines. Beyond the simple fact that entities are discrete in one case and continuous in another, is the fact that this affects the detailed dynamics of the line. For instance, suppose the machine Si is starved because of the failure of Si-1. In a DT line, after Si-ll is repaired, Si still has to wait for a workpiece to complete its processing at Si-1 before Si can start working again. On the other hand, in a CT line, as soon as Si-i is repaired it starts producing, and Si is also immediately operational. These (and other) differences in the dynamics mean that we can certainly expect some discrepancies between the detailed performance measures of the lines. Also, as we said in section 3, we deliberately did not try to modify the CT model to mimic some of these details, since we wanted to keep the CT model simple. What is interesting then, is that while the detailed utilizations of the machines are not as accurate, we find the main PMs (throughput, and WIP) are still surprisingly accurate, given the fundamental difference in the physical behavior of the lines. In fact, for practical system design, throughput and WIP are the main PMs of interest. These points are discussed further in section 7. We also remark on the relation between the time scale and the accuracy of the CT model for prediction of PMs in the DT lines. As discussed in Gershwin (1987), De Koster and

156

R. SURI AND B.-R. FU

Table 8. PMs of DT line and CT line simulations.

Throughput Case 1 Case 2 Case3 Case 4

CTL Simulation

DTL Simulation

0.125 4- 0.001 0.561 4- 0.003 2.1104-0.008 2.331 4- 0.008

0.125 4- 0.001 0.574 4- 0.003 2.1254-0.008 2.332 4- 0.008

Error

0.0% -2.3% -0.7% -0.1%

WIP Case 1 Case 2 Case3 Case 4 Avg. Buffer Case 1 - B 1 Case2-B1 Case 3 B1 B2 B3 B4 B5 Case 4 BI B2 B3 B4 B5

5.97044.042 411.2124329.940 4-

0.001 0.014 0.112 17.847

5.978 44.054411.3344331.004 4-

0.001 0.014 0.109 17.463

-0.1% -0.3% -1.1% -0.3%

3.970 4- 0.001 2.3404-0.011

3.978 4- 0.001 2.3354-0.011

-0.2% 0.2%

1.406 0.894 0.871 0.857 2.204

1.383 0.892 0.878 0.875 2.262

4- 0.016 4- 0.017 4- 0.021 4- 0.024 4- 0.027

1.7% 0.2% -0.8% -2.1% -2.6%

44444-

0.0% -0.4% -0.5% -0.3% -0.5%

4- 0.016 -4- 0.017 4- 0.022 4- 0.023 4- 0.026

52.3 454.4 463.4495.5 458.7 4-

1.9 1.9 5.8 12.2 3.5

52.3 54.6 63.7 95.8 59.0

1.9 2.0 5.8 12.0 3.7

Table 9. Machine utilizations of DT line and CT line simulations (Case 1 and 2). Case

1

Utilization Machine 1 Block Cycle Down Machine 2 Cycle Down Starve Machine 1 Block Cycle Down Machine 2 Cycle Down Starve

CTL Simulation

DTL Simulation

Error

0.857 4- 0.001 0.125 4- 0.001 0.018 4- 0.001

0.857 4- 0.001 0.125 4- 0.001 0.018 4- 0.001

0.0% 0.0% 0.0%

0.125 4- 0.001 0.875 4- 0.001 0

0.125 4- 0.001 0.875 4- 0.001 0

0.0% 0.0% --

0.253 + 0.002 0.187 4- 0.001 0.560 4- 0.002

0.236 4- 0.002 0.191 4- 0.001 0.573 4- 0.002

7.2% -2.1% -2.3%

0.1124-t:0.001 0.590 4-4-0.002 0.298 :t: 0.002

0.115+0.001 0.604 -t- 0.002 0.281 4- 0.003

-2.6% -2.3% 6.0%

ON USING CONTINUOUS FLOW LINES TO MODELDISCRETEPRODUCTIONLINES

157

Wijngaard (1989), and Dallery and Gershwin (1992), there are two time scales in tandem lines: the time scale of the operations, and the time scale of the failures and repairs. It is believed by some researchers that if the time scale of the failures and repairs is much larger than the time scale of the operations (at least one order of magnitude) then the CT model will be agood approximation of the DTline(e.g. see Dallery and Gershwin 1992). However, this statement seems overly simplistic; the actual issue is more subtle and additional research is needed to make such a statement accurate. For example, De Koster and Wijngaard (1989) present a 2-machine line for which the error in throughput follows a bell-shaped curve as the time scale of the second machine decreases. Suri and Fu (1993) study a 15-machine line with 192 combinations of parameter settings, and again it is shown that the error in throughput follows a more complex pattern than being a monotonic function of the time scales. Clearly, any a priori statement about accuracy would be useful for practitioners, but equally clearly, we see this to be a topic that needs additional work. Table 12 lists the computer run times for the DT and CT line simulations on a Sun SPARC station. The last column shows the ratio of the computer run times. It can be seen that, for the longer lines, CT line simulations run much faster than DT line simulations. This is of particular interest for designers of production lines such as transfer lines, which can have tens and even hundreds of stations (Ho et al. 1983, and Wei et al. 1989). Therefore we next explore further this issue of run times.

5.3. An Extensive Study of Run Times for One Configuration We present here an extensive comparison of the run times for DT and CT line simulations. We will do this for a large number of configurations of 15-machine lines. We choose this length because most practical production lines are much longer than this, and if our study shows significant advantages of the CT line approach for this length, then the advantages will be more pronounced for longer lengths. We intend to fully explore other line lengths in the future. However, as will be seen, even studying one line length extensively requires a large amount of computation time, so we do this as a first step. It should be pointed out that the study by Alvarez et al. (1991) also focused on only one line length. In order to study the run times, we used an experimental design which resulted in 192 configurations of lines with randomly generated parameters. The complete details for the experimental design can be found in Suri and Fu (1993). Here we simply mention that the design covered all combinations of balanced, unbalanced, homogeneous and nonhomogeneous lines. Further, all the CT line parameters (flow rates, mean operating volumes to failure, repair rates, buffer capacities) are varied randomly according to the specified design, in order to generate the 192 configurations. Note that Alvarez et al. (1991) only considered lines with machines having identical reliability parameters. Our experimental design also produces cases where machines have non-identical reliability parameters. We studied the run times for 10 independent replications with warm up and production quantities as before. For each line configuration, we calculated the ratio of the DT line simulation run time to that of the corresponding CT line. Figure 5 is a histogram of these ratios. For this set of experiments we see that, except for three cases where the ratios are slightly less than 1, CT line simulations run much faster than DT line simulations. Two-

158

R. SURf AND B.-R. FU

Table 10. Machine utilizations of DT line and CT line simulations (Case 3),

Utilization

CTL Simulation

DTL Simulation

Error

Machine I Block Cycle Down

0.187 4- 0.003 0.751 -4- 0,003 0.062 4- 0.002

0.182 • 0,003 0.756 + 0,003 0.062 + 0,002

2.7% -0.7% 0.0%

Machine 2 Block Cycle Down Starve

0.180 0.591 0.065 0.164

4- 0.003 4- 0,002 4- 0.003 4- 0.002

0.194 0.595 0.066 0.145

4- 0.003 4- 0.002 4- 0.003 4- 0.002

-7.2% -0.7% 1.5% 13,1%

Machine 3 Block Cycle Down Starve

0.158 0,591 0.067 0.185

4- 0.004 4- 0,002 4- 0.003 4- 0.004

0.165 0.595 0.067 0.173

4- 0.004 4- 0.002 4- 0.003 4- 0.004

-4.2% -0.7% 0.0% 6.9%

Machine 4 Block Cycle Down Starve

0.127 0.591 0.065 0.217

4- 0.005 4- 0.002 4- 0.002 4- 0,005

0.135 0.595 0.065 0.205

4- 0.005 4- 0.002 4- 0,002 4- 0,004

-5.9% -0.7% 0.0% 5,9%

Machine 5 Block Cycle Down Starve

0.097 0.591 0.066 0.246

4- 0.002 4- 0.002 4- 0.003 4- 0.005

0.108 0.595 0,066 0,231

4- 0,002 + 0,002 • 0,003 4- 0.005

- 10.2% -0.7% 0,0% 6.5%

0.737 + 0.003 0,060 4- 0.001 0,202 -t- 0,003

-0.7% 0.0% 3.0%

Machine 6 Cycle

Down Starve

0,732 4- ~,003 0.060 4- 0.001 0.208 d: 0.003

thirds of the CT lines run at least 10 times faster than the DT lines. In ten of these cases, the ratios even go above 80. Thus we can expect that, for long production lines, in almost all cases the CT line simulation approach is the preferred one, particularly in light of the demonstrated accuracy of our translation scheme for PMs. We can give an intuitive explanation for the difference in computer run times. The reason for the large ratios comes from the fact that, with certain parameter settings, the number of events simulated in CT line models is much smaller than that in DT line models. To illustrate this, consider the previous Figure 4, where each of the vertical solid lines in the machine sample paths represents the occurrence of an event. We can see that the number of events in the sample path of the CT line is less than its counterpart in the sample path of the DT line. It is also worth speculating on why the ratio of run times is less than one in a few situations. In those cases, the number of events is about the same, but the CT line simulation requires more divisions/multiplications for event scheduling compared with the DT line simulation which requires mostly additions/subtractions. This difference stems

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

159

Table 11. Machine utilizations of DT line and CT line simulations (Case 4).

Utilization

CTL Simulation

DTL Simulation

Error

Machine 1 Block Cycle Down

0.076 4- 0.004 0.816 4- 0.002 0.108 4- 0.004

0.076 4- 0.004 0.816 4- 0.002 0.108 4- 0.004

0.0% 0.0% 0.0%

Machine 2 Block Cycle Down Starve

0.098 4- 0.009 0.583 4- 0.002 0.269 • 0.008 0,050 4- 0.006

0.100 4- 0.009 0.583 4- 0.002 0.269 4- 0.008 0.047 4- 0.006

-2.0% 0.0% 0.0% 6.4%

Machine 3 Block Cycle Down Starve

0.038 4- 0.008 0.699 4- 0.002 0.204 4- 0.007 0.059 4- 0.006

0.038 4- 0.008 0.699 4- 0.002 0.204 4- 0.007 0.059 4- 0.007

0.0% 0.0% 0.0% 0.0%

Machine 4 Block Cycle Down Starve

0.023 4- 0.008 0.745 4- 0.002 0.159 4- 0.008 0.073 4- 0.009

0.023 4- 0.005 0.745 4- 0.002 0.159 4- 0.002 0.073 4- 0.004

0.0% 0.0% 0.0% 0.0%

Machine 5 Block Cycle Down Starve

0.003 4- 0.002 0.699 4- 0.002 0.232 4- 0.011 0.067 4- 0.012

0.003 4- 0.002 0.700 4- 0.002 0.232 4- 0.011 0.066 4- 0.012

0.0% -0.1% 0.0% 1.5%

Machine 6 Cycle Down Starve

0.723 4- 0.002 0.181 4- 0.005 0.096 4- 0.005

0.723 • 0.002 0.181 4- 0.005 0.096 4- 0.005

0.0% 0.0% 0.0%

Table 12. Computer run times for simulations of Cases 1 ~ on a Sun SPARC station (in CPU seconds).

Case

DT Simul~ion (1)

CT Simulmion (2)

Rmio (= (1)/(2))

1 2 3 4

196 181 564 553

230 218 107 55

0.85 0.83 5.29 9.99

f r o m the G S M P f o r m u l a t i o n u s e d as the basis for the C T s i m u l a t i o n m o d e l . H o w e v e r , w e also w i s h to e m p h a s i z e that these are heuristic e x p l a n a t i o n s , and our o n g o i n g r e s e a r c h a i m s to c h a r a c t e r i z e m o r e clearly t h o s e c o n f i g u r a t i o n s o f lines w h e r e the C T line a p p r o a c h will b e m o r e efficient.

160

R. SURI AND B.-R. FU

70 60 5O o

40 30 20

X.\\\NI

o 0.9

,

1

lO

,

,

20

3o

40

50

60

70

80

90

100

Ratio of Computer R u n Times Figure 5. Ratio of computer run times of DT and CT lines.

6.

Preliminary Results on Gradient Estimation for CT Lines

In this section we preview some of our results on gradient estimation for CT lines. The aim here is to show that, although the IPA technique for gradient estimation was originally developed for systems with discrete entities, it can be applied to CT lines as well. We consider estimating the gradient of steady state throughput with respect to the flow rates of the M machines (i.e., the M-vector of values O TP/OCi, for i = 1. . . . . M). The IPA algorithm for calculating this gradient vector is concise and easily inserted in the simulation code, as seen in Appendix 2. The derivation of the algorithm will not be discussed here, since we are only interested in illustrating some of the potential of CT line models. For a complete discussion on the IPA algorithm for CT lines see Suri and Fu (1992). We have implemented this IPA algorithm in our simulation algorithm for CT lines. While we are still working on an analytical proof of consistency of this gradient estimate, preliminary numerical results show that it seems accurate. As an example, we compare the IPA estimate with a central finite difference (CFD) estimate obtained by multiple simulation runs for a 6-machine non-homogeneous line (denoted "CFD+Simulation") taken from Dallery et al. (1989). The input parameters for this line (Case 4) are listed in Table 5. (We chose the finite difference ~ = 4-0.01Ci to obtain CFD estimates.) The numerical results are shown in Table 13. The IPA gradient vector estimate is obtained from observing 10 independent replications. For each replication, the simulation is run for a quantity of 100,000 products (warm up) to eliminate the effect of the initial transient, and then statistics are collected for a production quantity of 5,000,000. On the other hand, to get n gradient estimates using CFD we need 2n sets of replications, each set having the run lengths just stated. The gradient estimates are shown in the table along with their 95% confident intervals. Note that we used common random numbers to obtain the

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

161

Table 13. Comparison of IPA and CFD+simulationgradient estimates. Case 4 0 TP/OCj 0 TP/OC2 TP/OC3 0 TP/OC4 TP/OC5 TP/OC6

IPA (1.427 4- 0.028) • (4.513 4- 0.052) x (1.029 4- 0.009) • (9.143 4- 0.145) x (8.148 4- 0.234) x (%849 4- 1.143) x

CFD+Simulation 10-l 10-2 10-1 10-2 10-2 10-3

(1.434 • 0.024) x (4.518 4- 0.048) • (1.029 4- 0.010) x (9.176 4- 0.123) x (8.193 4- 0.198) • (7.883 4- 1.230) x

10-1 10-~10-1 I0 -e 10-2 10-3

Simulation time for computing entire gradient vector: IPA: 5090 seconds CFD: 33822seconds CFD/IPA = 6.64

C F D gradient estimates, as suggested by L'Ecuyer and Perron (1992), in order to reduce variance. As shown in the table the widths of the 95% confident intervals for both CFD and IPA estimates are about the same. It can be seen that the IPA estimate compares well with the "CFD+simulation" estimate. However, the CFD+simulation estimate takes 6.6 times the computational effort of the IPA estimate. (The CPU times are for a Sun SPARC workstation.) Thus we see that the IPA algorithm already proven for DT lines (Ho et al. 1979, Ho et al. 1983, Cao and Ho 1987) offers substantial computational savings for gradient computation in CT lines as well.

7.

Conclusion

We have formalized and extended previous work to show that DT lines can be reasonably represented by CT models. From a practical point of view our results show that the CT model should be an acceptable analysis tool for an actual production line. The reason is that one can expect much larger errors in the estimates of input parameters (especially failure and repair rates). The inaccuracies introduced by such errors would overshadow those observed in Tables 8-11. We have, in fact, experimented with modifying the behavior of the CT model in attempts to make it behave "more like a DT line". For instance, in the example at the end of section 5, we can implement a restriction in the model that whenever a machine (such as Si above) is starved, it has to wait for a whole unit of product to arrive (i.e. the buffer level must reach unity) before it can start working. Other similar ideas have been explored too. Some of these implementations do appear to give estimates that are somewhat closer to the PMs of the DT line, but at the cost of several complications. First, we realize that all such ideas introduce more types of events as well as the consideration of a discrete quantity (a whole product) into CT models. The resulting model will be more like a simulation model of a DT line which implies that more events would have to be simulated. This is contrary to one o f our goals in using the CT model. Second, the translation o f PMs back to the DT line becomes more complicated as well. Third, we have successfully implemented IPA

162

R. SURI AND B.-R. FU

algorithms for gradient estimation on the basic CT line. With the more complex simulations, implementation of any type of PA algorithms will be much less obvious. Fourth, the CT model that we analyze corresponds to the way an actual CT line would physically operate. So our present and ongoing work (e.g. gradient estimation) can be used by people designing such lines as well. This would not be the case if we modified the CT model with special rules. Finally, the CT model that we described in this paper provides a simple and intuitive way to translate the input parameters and PMs from/to a DT line. We feel the ease and simplicity of translation/implementation, and speed of the simpler simulation, outweigh the added accuracy that might be achieved. Our study also suggests the possibility of applying optimization algorithms to CT models as an approximation to the optimization of DT lines (e.g. the optimal design of buffer sizes and cycle times). It is conceivable that one could use a gradient estimation algorithm, along with the "single run" optimization methods as in Suri and Leung (1987, 1989) on a CT model, and then through a final translation algorithm obtain fast estimates of optimal parameter settings for DT lines. Our ongoing research is exploring these ideas (Suri and Fu 1993). We hope our study will provide the foundation for such research, and also stimulate additional work on this topic.

Acknowledgment This work was partly supported by NSF Grant DDM-9201813, and by a grant from Ford Motor Company. We thank several anonymous reviewers for pointing out missing references, and for their constructive comments on our original manuscript.

Appendix 1: Notation and Simulation Algorithm The definitions for the following variables are described in the text: M, t, e*, q, Q, Ci, vi, oti, B(j), Wi, Ui, and xj, for i = 1. . . . . M, and j = 1. . . . . M - 1, where the argument t is dropped for notational convenience. The following notation is used in the CT line simulation algorithm.

INF: At: EM[i]: EB[j]: TM[i]: TB[i]: KM, KB: K*: L(j):

A very large constant. Time to the occurrence of triggering event. Next possible event at Si. Next possible event at Bj, and EB[M] denotes the possible event T/ if vM # 0. Time to the occurrence of EM[i]. Time to the occurrence of EB[i], and TB[M] is the time to the occurrence of EB[M]. Constants. Index of triggering event. Accumulator for average level of Bj.

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

163

or, C, v, W, U, x, and L are the corresponding column vectors. The simulation algorithm is as follows: O. INITIALIZATION t = 0; At = 0; ~_ = O; x = 0; L = 0; q = Q; x0 = INF; B ( M ) = INF; 1)1 = for i = 2 to M do vi = min(vi_l, Ci); Generate W; EB[M] = 7-f ; 1.

NEXT EVENT AT MACHINES AND BUFFERS (local events) for j =

ltoM-ldo

case

oj > vj+l: {TB[j] = s(j)-x/ uj-v~+, ; E B [ j ] = B.Tj; } vj < vj+l: {TB[j] = vj+l-vj xj '. EB[j] = 13Cj; } vy = vy+l: TB[j] = I N F ; end case for/ = 1 t o M d o { if (oti ~= D) then if(vi = 0) then TM[i] = INF; else {EMIl] = ~,.; TM[i] = W i / v i ; } else {EM[i] = T~i; TM[i] = Ui; }

} if (VM ~ O) then TB[M] = -~-" UM

else TB[M] = INF; 2.

NEXT EVENT IN CTL (global event) KM = argmin(TM[i]); i = 1,2 . . . . . M. K s = argmin(TB[i]);

i =

l, 2 . . . . . M.

if ( T M [ K M ] > T B [ K s ] ) then {K* = K s ; A t = TB[K]; e* = EB[K]; } else {K* = KM; A t = TM[K]; e* = EM[K]; } 3.

UPDATE t=t+At; q = q - u M 9 At;

Cl" ,

164

R. SURI AND B.-R. FU

W=W--v.

At;

u=u--At-1; for j =

ItoM-ldo

{temp = xj; xj = xj + (vj - vj+l) 9 At; L ( j ) = L ( j ) +

(temp+x)). At. 2 ' }

case e* = ~'K* :

otK, = D; vm = 0; Ui = s a m p l e ( 1 / r i ) ; Update_Downmacs (K*); Update_Up_mcs (K* - 1);

e* = ~K*:

Set_Flow_Rate (K*); Wi = sample(wi); Update_Down_mcs (K*); Update_Up_mcs (K* - 1);

e* = B ~ x * :

Update_Up_mcs (K*);

e* = BEK*:

Update_Down_mcs (K*);

end case. 4.

TEST if (e* = Tf) then go to OUTPUT; else go to NEXT EVENT AT MACHINES A N D BUFFERS;

5.

OUTPUT Throughput = ~ ; Average buffer levels

L. t'

p r o c e d u r e Set_Flow_Rate (n) case x n - l = 0 and Xn = B(n): on = min(vn_l, Cn, V~+l);

xn-1 = 0 and xn < B(n): On = min(vn_l, C,); Xn-1 > 0 and xn = B(n): on = min(C,,, vn+l); otherwise:

{vn = Cn; o t i = O; }

end case; if (Vn = vn-1) then an = S; i f ( v , = On+l)then otn = B; p r o c e d u r e Update_Down_mcs (n) while ((xnis empty) and (n < M - 1)) do {

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

165

n=n+l; S e t , l o w _ R a t e (n);

} procedure Update_Up_mcs (n) while ((Xn is full) and (n > 1)) do { n=n-l; Set_Flow_Rate (n);

} Appendix 2: Algorithm for Computation of the Gradient Vector of a CT Line with non-homogenous flow rates The notation is the same as defined in Appendix 1. TSum[k], QSum[i][k] and USum[i][k] are PA accumulators where i is the index of the machine and k is the index of the machine with perturbed flow rate. Gt is the gradient of throughput with respect to Ck. dvdc(i, j), i, j = 1 . . . . M are indicator functions where dvdc(i, j ) = 1, if v[i] = Cj, and 0, otherwise. Insert this after Step 0 in the Simulation Algorithm. 0A. I N I T I A L I Z A T I O N f o r k = 1 to M do {

TSum[k] = 0; f o r i = 1 to M do {

QSum[i][k] = 0; USum[i][k] = 0;

} } Insert 3A after Step 3 in the Simulation Algorithm. 3A. PERTURBATION GENERATION If (e* = St'K,) then fork=ltoMdo{

temp[k] = - ~-Wri l 9 Iv 9dvdc( K*, k) + QSum[ K*][k ]] ;

166

R. SUR1 AND B.-R. FU

USum[ K* ] [k] = O;

} If (e* = 7EK,) then for k = 1 to M do temp[k] = USum[K*][k]; If (e* =/35rK,) then fork=

ItoMdo

temp[ k ] =

l v[K*I-v[K*+I] 9 {r 9 [ d v d c ( K * , k) - d v d c ( K * + 1, k)] +QSum[K*][k] - Q S u m [ K * + i][k]};

If (e* = BEK,) then fork=

ltoMdo

temp[ k ] =

t v[K*WI]-v[K*] 9 {r 9 [ d v d c ( K * + 1, k) - d v d c ( K * , k)] + Q S u m [ K * + l][k] - QSum[K*][k];

If (e* = ~ ) then for k = I to M do temp[k] = - viM] I 9 {3 - d v d c ( M , k) + QSum[M][k]}"

fork=

ltoMdo{

TSum[k] = TSum[k] + temp[k];

for/=

l to M d o {

if (Si is failed) then USum[i][k] = USum[i][k] - temp[k];

else QSum[i][k] = QSum[i][k] + r 9 d v d c ( i , k) + v[i] . temp[k];

I l Insert P1 in procedure Update_Down_mcs(n) before Set_Flow_Rate(n). PI. PERTURBATION PROPAGATION (downstream) fori=

ltoMdo

QSum[n][i] = QSum[n - 1][i];

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

-] 67

Insert P2 in p r o c e d u r e U p d a t e _ U p _ m c s ( n ) b e f o r e Set_Flow_Rate(n). P2. P E R T U R B A T I O N P R O P A G A T I O N ( u p s t r e a m ) f o r i = 1 to M d o

QSum[n][i] = Q S u m [ n + 1][i]; Insert this after Step 5 in the S i m u l a t i o n A l g o r i t h m . 5A. O U T P U T for k = 1 to M d o

Gk = - t ~ 9 TSum[k];

References R. Akella and E R. Kumar. Optimal control of production rate in a failure prone manufacturing system. IEEE Trans. Aurora. Control, 31(2):116-126, 1986. T. Altiok and S. Stidham, Jr. A note on transfer lines with unreliable machines, random processing times, and finite buffers, l i e Trans., 14(2): 125-127, 1982. R. Alvarez, Y. Daltery, and R. David. An experimental study of the continuous flow model of transfer lines with unreliable machines and finite buffers. In IMACS-IFAC Symposimn, Modelling and Control of Technological Systems, Lille, France, 1991. J. Banks and J. S. Carson, II. Discrete-EventSystem Simulation. Prentice-Hall, 1984. J. A. Buzacott. Automatic transfer lines with buffer stocks. Int. J. Prod. Res., 5(3): 183-200, 1967. J. A. Buzacott and L. E. Hanifin. Models of automatic transfer lines with inventory banks: A review and comparison. AIIE Trans.. 10(2):197-207, 1978. X. R. Cap and Y. C. Ho. Sensitivity analysis and optimization of throughput in a production line with blocking. IEEE Trans. Aumm. Control., AC-32(I 1):959-967, 1987. X. R. Cap and Y. C. Ho. Models of discrete event dynamic systems. IEEE ControISystems Magazine, 10(4):69-76, 1990. M. Caramanis. Production system design: A discrete event dynamic system and generalized Benders' decomposition approach. Int. J. Prod. Res., 25(8):1223-1234, 1987. Y. E Choong and S. B. Gershwin. A decomposition method for the approximate evaluation of capacitated transfer lines with unreliable machines and random processing times, liE Trans., 19(2):150-159, 1987. Y. Dallery, R. David, and X.-L. Xie. An efficient algorithm for analysis of transfer lines with unreliable machines and finite buffers. IIE Trans., 20(3):280-283, 1988. Y. Dallery, R. David, and X.-L. Xie. Approximate analysis of transfer lines with unreliable machines and finite buffers. IEEE Trans. Aurora. Control., 34(9):943-953, 1989. Y. DaUery and S. B. Gershwin. Manufacturing flow line systems: A review of models and analytical results. In Qtteueing Systems: Theory and Applications, Special Issue on Queueing Models of Manufacturing Systems, Vol. 12, pp. 3-94, 1992. H. D'Angelo, M. Caramanis, S. Finger, A. Mavretic, Y. A. Phillis, and E. Ramsden. Event-driven model of an unreliable production line with storage. In Proc. 24th IEEE Conf. on Decision and Control, Ft. Lauderdale, FL, December, pp. 1694-1698, 1985. H. D'Angelo, M. Caramanis, S. Finger, A. Mavretic, Y. A. Phillis, and E. Ramsden. Event-driven model of unreliable production lines with storage. Int. J. Prod. Res., 26(7): 1173-1182, 1988. R. David, X. L. Xie, and Y. Dallery. Properties of continuous models of transfer lines with unreliable machines and finite buffers. IMA J. of Math. Applied hz Business and Indnstt3,, 6:281-308, 1990. R. De Koster and J. Wijngaard. Continuous vs. discrete models for production lines with blocking. In H. G. Perros and T. Altiok, editors, Queueing Networks with Blocking. North-Holland, 1989.

168

R. SURI AND B.-R. FU

D. Dubois and J. E Forestier. Productivit6 et en-cours moyens d'un ensemble de deux machines s6par6es par une z6ne de stockage. RAIRO Automatique, 16(2): 105-132, 1982. M. C. Fu and J.-Q. Hu. Extensions and generalizations of smoothed perturbation analysis in a generalized semiMarkov process framework. IEEE Trans. Autom. Control., AC-37(2): 1483-1500, 1992. S. B. Gershwin. An efficient decomposition method for the approximate evaluation of tandem queues with finite storage and blocking. Oper. Res., 35(2):291-305, 1987. S. B. Gershwin. An efficient decomposition algorithm for unreliable tandem queueing systems with finite buffers. In H. G. Perros and T. Altiok, editors, Queueing Networks with Blocking. North-Holland, 1989. S. B. Gershwin and O. Berman. Analysis of transfer lines consisting of two unreliable machines with random processing times and finite storage buffers. AIIE Trans., 13(1):2-11, 1981. S. B. Gershwin and I. C. Schick. Continuous model of an unreliable two-stage material flow system with a finite interstage buffer. MIT Technical Report LIDS-R-1039, 1980. S. B. Gershwin and I. C. Schick. Modeling and analysis of three-stage transfer lines with unreliable machines and finite buffers. Oper Res., 31(2):354-380, 1983. P. Glasserman. Gradient Estimation Via Perturbation Analysis. Kluwer Academic Publishers, 1991. P. W. Glynn. Optimization of stochastic systems. In Proc. 1986 Winter Simulation Co~, pp. 52-59, 1986. P. W. Glynn. Likelihood ratio gradient estimation: An overview. Proc. 1987 Winter Simulation Conf., pp. 366-374, 1987. P. W. Glynn. A GSMP formalism for discrete event systems. Proc. IEEE., 77(1):14-23, 1989. P. W. Glynn. Pathwise convexity and its relation to convergence of time-average derivatives. Management Sci., 38(9):1360-1366, 1992. P. W. Glynn and W. Whitt. Indirect estimation via L = ~.W. Oper. Res., 37(1):82-103, 1989. Y. C. Ho. Performance evaluation and perturbation analysis of discrete event dynamics systems. IEEE Trans. Autom. Control., AC-32(7):563-572, 1987. Y. C. Ho, M. A. Eyler, and T. T. Chien. A gradient technique for general buffer storage design in a production line. Int. J. Prod. Res., 17(6):557-580, 1979. Y. C. Ho, M. A. Eyler, and T. T. Chien. A new approach to determine parameter sensitivities of transfer lines. Management Sci., 29(6):700-714, 1983. J.-Q. Hu. Convexity of sample path performance and strong consistency of infinitesimal perturbation analysis estimates. IEEE Trans. Autom. Control., 37(2):258-262, 1992. V. S. Kouikoglou and Y. A. Phillis. An exact efficient discrete-event model for production lines with buffers. Proc. 28th IEEE Conf. on Decision and Control, Tampa, FL, pp. 65-70, December 1989. V. S. Kouikoglou and Y. A. Phillis. An efficient discrete-event model for production networks of general geometry. Proc. 29th IEEE Conf. on Decision and Control, Honolulu, HI, pp. 3446-3451, December 1990. V. S. Kouikoglou and Y. A. Phillis. An exact discrete-event model and control policies for production lines with buffers. IEEE Trans. Autom. Control., 36(5):515-527, 1991. V. S. Kouikoglou and "1(.A. Phillis. A serial finite unreliable queue model for production lines. Proc. 31stlEEE Conf. on Decision and Control, Tucson, AZ, pp. 1659-1664, December 1992. A. M. Law and W. D. Kelton. Simulation Modeling and Analysis. McGraw-Hill, 1982. P. L'Ecuyer and G. Perron. On the convergence rates of IPA and FDC derivative estimators. Submitted. for publication, 1992. R. Malham6 and E. K. Boukas. Transient and steady-states of the statistical flowbalance equations in manufacturing systems. In Proc. Third ORSA/TIMS Conf. on Flexible Manufacturing Systems: Operations Research Models and Applications. K. E. Stecke and R. Suri eds. Elsevier, 339-345, 1989. D. Mitra. Stochastic theory of a fluid model of producers and consumers coupled by a buffer. Adv. Appl. Prob., 20:646-676, 1988. C. D. Pegden, R. E. Shannon, and R. D. Sadowski. Introduction to Simulation Using SIMAN. McGraw-Hill, 1990. Y. A. Phillis and V. S. Kouikoglou. Techniques in modeling and control policies for production networks. In C. T. Le.ondes, Control and Dynamic Systems, Vol. 47, Part 3, pp. 1-50, 1991. E. Ramsden, M. Ruane, A. Mavretic, M. Caramanis, and H. D'Angelo. The use of hardware simulators in modeling production networks. Large Scale Systems, Vol. 11, pp. 149-164, 1986. S. M. Robinson. Convergence of snbdifferentials under strong stochastic convexity. Submitted to Management Science, 1993. B. A. Sevast'yanov. Influence of storage bin capacity on the average standstill time of a production line. Theory of Probability and Its Applications, 7(4):429--438, 1962.

ON USING CONTINUOUS FLOW LINES TO MODEL DISCRETE PRODUCTION LINES

169

A. Sharifnia. Production control of a manufacturing system with multiple machine states. IEEE Trans. Autom. Control., 33(7):620-625, 1988. A. Sharifnia, M. Caramanis, and S. B. Gershwin. Dynamic set-up scheduling and flow control in flexible manufacturing systems. In K. E. Stecke and R. Suri, editors, Proc. Third ORSA/TIMS Conference on Flexible Manufacturing Systems: Operations Research Models and Applications, pp. 327-332. Elsevier, 1989. R. Suri. Infinitesimal perturbation analysis for general discrete event systems. Z ACM, 34(3):686-717, 1987. R. Suri. Perturbation analysis: The state of the art and research issues explained via the GI/G/I queue. Proc. IEEE, 77(1):114-137, 1989. R. Suri and G. W. Diehl. A variable buffer-size model and its use in analyzing closed queueing networks with blocking. Management Sci., 32(2):206-224. R. Suri and B.-R. Fu. Flow rate gradient estimation for continuous flow production lines. University of WisconsinMadison. Working paper, 1992. R. Suri and B.-R. Fu. Using continuous flow models to enable rapid analysis and optimization of discrete production lines--A progress report. In Proc. 19th Annual NSF Grantees Conf. on Design and Manufacturing Systems Research, pp. 1229-1238, 1993. R. Suri and Y. T. Leung. Single run optimization of a S1MAN model for closed loop flexible assembly systems. In Proc. 1987 Winter Simulation Conf., pp. 738-748, 1987. R. Suri and Y. T. Leung. Single run optimization of discrete event simulations-An empirical study using the M/M/I queue, l i e Trans., 21(1):35--49, 1989. R. Suri, J. Sanders, and M. Kamath. Performance evaluation of production networks. In S. C. Graves, A. H. G. Rinnooy Kan, and P. H. Zipkin, editors, Handbooks in Operations Research, Vol. 4. Elsevier, 1993. K. C. Wei, Q. Q. Tsao, and N. C. Otto. Determining buffer size requirements using stochastic approximation methods. Technical Report SR-89-73. Ford Research, 1989. W. Whitt. Continuity of generalized semi-Markov processes. Math. of Oper. Res., 5(4):494-50 I, 1980. J. Wijngaard. The effect of interstage buffer storage on the output of two unreliable production units in series, with different production rates. AIIE Trans., 11(t):4247, 1979. S. Yeralan, W. E. Franck Jr., and M. A. Quasem. A continuous materials flow production line model with station breakdown. European J. of Operational Res., 27:289-300, 1986.