1. Introduction

The study of production systems was motivated by the rapid development of manufacturing automation in the past decades. Obtaining accurate estimates of a system's performance, such as the average throughput rates and buffer levels, is a costly but necessary task for the design and optimal operation of any modern factory.

Simulation is a class of solution techniques that, although accurate, suffer from high computation times. In this paper we develop a hybrid simulation/analytic model for production systems. The model is a flow approximation of discrete-part production networks but it turns out to be quite accurate and faster than conventional simulators.

2. Review of the literature

Queueing network analysis and simulation have been used to describe complex production systems. Exact Markovian modeling of such systems is rarely efficient because the dimensions of the model become enormous even for simple systems. Decomposition is another analytic tool for approximating the performance of large systems. Here the network is broken into subsystems that are analyzed separately. The accuracy of the method depends crucially on the complexity of the network.

Exact solutions for the analysis of queueing networks with exponential interarrival and service time distributions were first reported by Jackson (1957). The steady-state probabilities have a sum-of-products form. Later these results were further extended to include more general networks. Decomposition methods (Kuehn, 1979; Takahashi et al., 1980; Bocharov, 1987), based on studying the first, second or higher moments of the subsystems in isolation, often give good approximations of performance measures.

Recently, simulation has been used to estimate performance sensitivities with respect to various design parameters. Perturbation analysis, which first appeared in Ho et al. (1979), is used in optimization problems with the aid of a single simulation run. However, simulation algorithms are very slow owing to their piece-driven logic. An attempt towards a new discrete-event logic (D'Angelo et al., 1988) for production lines led to a very fast simulator. The model owes part of its speed to the fact that it ignores various transients occurring during the operation of the line. The study of such transient phenomena was incorporated into another model, which is exact (Kouikoglou and Phillis, 1991), quite fast, and appropriate for the study of transient performance of long production lines.

3. System description

To illustrate the kind of production systems to be studied here, consider Fig. 1. The system is represented by a network of five nodes and directed arcs. The nodes represent workstations [W.sub.1], [W.sub.2],..., [W.sub.5]. To each arc we assign one buffer. Thus buffer [B.sub.1,4] is between workstations [W.sub.1] and [W.sub.4], which are called adjacent, the first one being the upstream, and the second one the downstream workstation.

Raw parts enter the network from an infinite source [B.sub.I] and after being processed they eventually leave the system on one of the directed arcs to the output buffer [B.sub.O]. A workstation [W.sub.j] consists of a number of parallel machines labelled [M.sub.i,j], which may produce at different rates, and fail and are repaired randomly. Each workstation performs a specific type of operation.

There are several part-types a, b, c, d, and e that require a fixed set of operations, each following a certain route in the production system. All types of raw parts are available from the source buffer [B.sub.I], and the output buffer [B.sub.O] can accept all the types of finished products.

Workstations [W.sub.4] and [W.sub.5] accept several types of semi-finished parts and combine them to produce composite part-types. A single part of the composite product needs a constant number of parts from each part-type, which will be called an assembly factor. For example, one item of type e in [W.sub.5] is composed of [[delta].sub.2,5] b-parts supplied by [W.sub.2], and [[delta].sub.3,5] c-parts from [W.sub.3]. Similarly composite parts are produced by the workstation [W.sub.4].

The flow of a part-type may also split into two or more streams. This is done by connecting a workstation with several buffers downstream. Workstation [W.sub.2] is the only splitting node of the production system. To each one of [B.sub.2,4] and [B.sub.2,5] we assign distribution factors (routing probabilities) [d.sub.2,4] and [d.sub.2,5], representing the long-run percentage of items completing service at [W.sub.2] that are routed to [B.sub.2,4] and [B.sub.2,5] respectively. However, if one of the two buffers happens to be full for some time, then the production of [W.sub.2] will be sent to the non-full buffer.

4. Assumptions

1. The production network is an acyclic graph, and has several part-types. Merging or splitting of part flows may take place at each node.

2. There can be any number of parallel machines in each workstation producing at deterministic rates. Although they perform the same kind of operation, their production rates need not be the same. The allocation of parts to the machines is made according to a first completed -- first supplied protocol. That is, a workpart arriving at a starved workstation will be sent to the machine that was starved first. Similarly if a workstation is blocked and a unit space is available downstream, a workpart will be released by the machine that had been blocked first.

3. A machine may fail during its operation on a workpart. A failed machine is repaired after an elapsed time. The time-to-failure and time-to-repair for each machine are exponential random variables with known means.

4. Between adjacent workstations there is a buffer of finite capacity, which accepts one class of parts. There is an infinite source at the beginning and an infinite sink of workparts at the end of the network.

Assumption 3 can be relaxed by considering more general distributions for the uptimes and the downtimes (gamma, Weibull, etc.). We assume that failures are operation dependent, that is, they occur only after a certain number of parts have been processed. Time-dependent failures, such as power failures, caused by exogenous factors can be incorporated into the model with minimum effort.

5. The discrete-event method

Piece-by-piece simulation models create an ordered list of every arrival or departure that occurs at a machine during a sample trajectory of the network, and carry out the actions appropriate for the next event. These events occur frequently and thus the cost of running piece-by-piece simulators may be high for most applications.

In our discrete-event model, we consider the following events:

a, machine_is_repaired;

b, machine_fails;

c, buffer_becomes_full;

d, buffer_becomes_empty;

e, buffer_becomes_not full;

f, buffer_becomes_not empty.

Each event corresponds to a change of the system state. The machines have two states (up and down), corresponding to the events a and b respectively, and the buffers have three (full, empty, and partly full), corresponding to events c, d, and the pair (e, f) respectively.

In discrete-part systems, machines may become starved or blocked, because they have different processing times. A starved or blocked machine operates on a workpart for a production cycle, then it becomes idle for some time waiting for a new part or an empty space, and then it becomes active again operating on the next workpart, and so on. The model is a continuous-flow approximation of a discrete-part network. As a result, blockage or starvation phenomena are treated by appropriately reducing the machine rate. We shall see later that this approximation produces excellent results in most practical situations. A machine that is up can be in any one of the following substates:

unforced, producing at full rate;

slowed down, producing at a reduced rate, because either it has limited supply of workparts, or all the downstream buffers alternate between full and partly full states for some time;

idle, if the reduced rate due to blockage or starvation happens to be zero. This phenomenon takes place when the workstation that causes the rate reduction is under repair.

In the following, idle states will be treated as special cases of slowed down states.

The main idea of the method is to use analysis to schedule the next event of a machine or buffer, rather than scheduling arrival and departures of every workpart separately. As will be made clear in the next section, the assumption of deterministic production rates is central in our analysis.

The above idea was first introduced in D'Angelo et al. (1988) for serial production systems with single machines between adjacent buffers. In that model it is possible for a sequence of simultaneous events to be examined as an aggregate event and execute it by one procedure. For example when a machine is repaired, and has a full buffer upstream, this buffer will become not full immediately, and the same not full events will be generated in the upstream buffers until a partly full buffer is encountered and the propagation is terminated. The dual not empty events propagate downstream until a partly full buffer is encountered. Also, when a buffer becomes full (empty) additional buffer_full (buffer_empty) events are generated upstream (downstream) until a partly full buffer is encountered.

However, in the more general systems examined here, the events propagate through several branches of the network instead of the single path of a line. To avoid the use of combinatorial search techniques for the allocation of these simultaneous events we suggest a simpler sequential procedure. After determining the next event that will occur in the network, the corresponding routine is executed and a number of instantaneous events are scheduled. These events are then executed successively in arbitrary order. Each event may also generate new instantaneous events that are treated likewise.

6. System parameters and state variables

Consider a production network with K buffers and N workstations. We use the following notation:

[n.sub.j]:         total number of machines in workstation
                   [W.sub.j], j= 1,2,...,N.
[NR.sub.i,j]:      nominal production rate of machine
                   [M.sub.i,j], i = 1,...,[n.sub.j].
1/[p.sub.i,j]:     mean time-to-failure of [M.sub.i,j] operating at
its
                   nominal rate.
1/[r.sub.i,j]:     mean time-to-repair of [M.sub.i,j].
[BC.sub.j,q]:      capacity of buffer [B.sub.j,q] that connects
workstation
[W.sub.j] to
                   [W.sub.q], j, q = 1,...,N.
[[delta].sub.j,q]  assembly factor of a part-type that is
transported
                   from [W.sub.j] and assembled into a composite
                   new part-type at [W.sub.q].
[d.sub.j,q]:    distribution factor of buffer [B.sub.j,q].

The production network is a discrete-event dynamic system (DEDS). The state vector consists of the following components:

(1) machine and buffer states and production rates, which remain unchanged between events, called the long-term variables;

(2) total production, buffer levels, remaining time-to-failure and time-to-repair which change continuously, called the short-term variables.

Thus long-term variables are event-driven whereas short-term variables are time-driven. We use the following notation for the system state at instant t:

M(i,j,t) - {1, if machine [M.sub.i,j] is up,

{0, if [M.sub.i,j] is under repair.

B(j, q, t,) = {1, if buffer [B.sub.j,q] is partly full,

= {2, if [B.sub.j,q] is full,

= {0, if [B.sub.j,q] is empty.

R(i, j, t): production rate of machine [M.sub.i,j]. RT(j, t): total production rate of workstation [W.sub.j]. RI(j, q, t): input flow rate of buffer [B.sub.j,q]. RO(j, q, t): output flow rate of [B.sub.j,q]. BL(i, q, t): level of buffer [B.sub.j,q]. TP(i, j, t): total production of machine [M.sub.i,j] by time t. PTF(i, j, t): remaining parts-to-failure of [M.sub.i,j]. TTR(i, j, t): remaining time-to-repair of [M.sub.i,j].

At any time t, the remaining time-to-failure TTF(i, j, t) of an operating machine [M.sub.i,j] is provided by

TTF(i, j, t) = PTF(i, j, t) / R(i, j, t),

by the assumption of operation-dependent failures.

7. Generation of events

When an event occurs, it causes the workstations to interact and triggers new events to occur in the system. In discrete-part systems, the time needed for the generation of new events, caused by buffers becoming full or empty and machine failures or repairs, is not negligible. For example, consider two machines and an empty intermediate buffer in a serial connection. Suppose that the upstream machine breaks on a workpart, whereas the last one is still working on another workpart. Then the buffer_empty event will occur after a transient time when the downstream machine finishes its workpart and is ready to receive one more. Such transient phenomena have been extensively studied in Kouikoglou and Phillis (1991) for production lines, and their consideration in modeling the system is of great importance, especially in cases of long lines with small buffer capacities. The resulting model is quite complex and its extension to production networks becomes prohibitive owing to its combinatorial nature.

This problem is resolved here by developing the continuous-flow model, which is the limit of the discrete-part model as the part size approaches zero. As a result of the continuous flow, the transient times needed for the generation of such events is infinitesimal. Next we describe the event generation rules:

1. Suppose that in the subsystem of Fig. 2. buffer [B.sub.j,q] becomes full. Then, its input rate is reduced instantly to its output rate. Workstation [W.sub.j] will attempt to send the overflowing production to other downstream buffers that are partly full. These buffers will tend to fill. However, if all the downstream buffers happen to be full, the workstation will be slowed down. As a result, the output rates of its upstream buffers [B.sub.i,j] will be reduced, and these buffers will tend to fill. Thus the buffer_full event will propagate upstream.

2. When buffer [B.sub.i,j] becomes empty, its output rate is reduced instantly. The downstream workstation [W.sub.j] will be slowed down and therefore the inventories of the downstream buffers [B.sub.j,p], [B.sub.j,q], etc., will tend to be exhausted. Thus the buffer_empty event propagates downstream and is dual to the buffer_full event. However, if different parts are assembled in workstation [W.sub.j], their flow rates will be reduced by the assembly protocol and the corresponding upstream buffers [B.sub.m,j] will tend to fill. Yet, if any of them happens to be full, a new buffer_becomes_full event will take place instantly.

3. By the assumption of operation-dependent failures, machines fail only during busy periods and do not deteriorate during idle periods. This means that a failure takes place right after a certain amount of production. As a consequence, when a machine is slowed down, its remaining time-to-failure is prolonged.

4. If a workstation is slowed down and one of its machines breaks down, then the production rates of the other machines will be increased. If, however, the workstation is neither starved nor blocked, its total production rate will be reduced instantly. The upstream buffers will tend to fill and the downstream ones will tend to empty.

5. In a dual fashion, if a machine is repaired and the workstation is not slowed down, then its upstream buffers will tend to empty and its downstream ones will tend to fill.

6. When [W.sub.j] is slowed down owing to blockage or starvation, the rates of its operating machines are reduced simultaneously. By the continuous flow principle and assumption 2 it turns out that the new machine rates are proportional to the nominal rates. For example, suppose that [W.sub.j] consists of two machines [M.sub.1,j] and [M.sub.2,j], which produce at maximum rates [NR.sub.1,j] and [NR.sub.2,j] respectively. If at time t the workstation happens to slow down to the rate RT(j, t), then the new machine rates are computed as follows:

(2) R(i,j,t) = RT(j, t) [NR.sub.i,j] / ([NR.sub.1,j] + [NR.sub.2,j]), i = 1, 2.

8. The DEDS model

Our system consists of K buffers and [Mathematical Expression Omitted] machines. Let [t.sub.k] be the time of the kth event occurrence in the system, and [X.sub.L] (t) and [X.sub.S](t) the vectors of long- and short-term states of the DEDS herein at time t. Let also E(t) denote the vector of next events that are scheduled to occur at every machine and buffer after time t, and TE(t) the vector of the associated event times. In the following we shall drop the argument t for convenience. The vector TE has dimension M + K and the form

TE = [[tau][M.sub.1,1] . . . [tau][M.sub.i,j] . . . [tau][B.sub.j,q]],

where [tau][M.sub.i,j] is the time of next event at machine [M.sub.i,j], and [tau][B.sub.j,q] is the time of next event at buffer [B.sub.j,q].

The proposed simulation model works as follows:

1. At time [t.sub.k], right after the execution of the kth event, the next event is determined. This is the one with the smallest event time. Thus

(1) [Mathematical Expression Omitted]

2. Right before the execution of the (k + 1)st event, the short-term states of machines and buffers related to that event, are updated. The corresponding update equations have the form

(2) [Mathematical Expression Omitted]

where [gamma] is a vector function whose explicit form will be derived in the next paragraph.

3. In the interval [t.sub.k], [t.sub.k+1]), the long-term variables remained unchanged. However, at time [t.sub.k+1] both [X.sub.S] and [X.sub.L] are partly adjusted. That is, the occurrence of the (k + 1)st event at a component causes changes to the states of adjacent components only. The corresponding adjusting equations are:

(3a) [X.sub.L]([t.sub.k+1]) = [symbol omitted]1([X.sub.L]([t.sub.k]), E([t.sub.k]), TE([t.sub.k])),

and

(3b) [Mathematical Expression Omitted]

[symbol omitted]1 and [symbol omitted]2 will be derived in subsequent paragraphs, and [xi] is a stochastic process. We can think of [xi] as representing all the random variables in the DEDS. In simulation it is a sequence of independent variables from the uniform distribution on [0,1] used to generate a possible realization of the uptimes and downtimes of machines, which belong to the short-term states of the system.

4. Having computed the operational characteristics of the components that are affected by the (k + 1)st event, it is then possible to find the next event that will occur at each one of them. As will be shown in following sections, by the event-generation rules and elementary algebra we can obtain the event-scheduling equations. These have the form

(4) E([t.sub.k+1]) = [phi])(Xs([t.sub.k+1]), [X.sub.L]([t.sub.k+1])),

and

TE([t.sub.k+1]) = 1 [t.sub.k+1] + [pi]([X.sub.S]([t.sub.k+1]), [X.sub.L]([t.sub.k+1])),

where 1 denotes the (M + K)-unit vector whose elements are all equal to one, and [phi], [pi] are functions to be derived in subsequent paragraphs.

It becomes clear that the above model represents a decomposable system, and we shall use this property in building the event-driven algorithm.

9. Event scheduling

In this section we develop explicit expressions for the event scheduling Equations 4 and 5. Consider the production network shown in Fig. 2.

9.1. Scheduling events in buffers

We examine a buffer [B.sub.j,q] at time t between successive events. Let [tau][B.sub.j,q] be the time of next event of the buffer.

Case a: If the input rate of buffer [B.sub.j,q] is greater than the ouput rate, then the buffer will fill in time

(6) [tau][B.sub.j,q] = t + ([BC.sub.j,q] - BL(j, q, t)) / (RI(j, q, t) - RO(j, q, t)),

for RI(j, q, t) > RO(j, q, t).

Case b: If the output rate of the buffer is greater than its input rate, then the buffer will empty in time

(7) [tau][B.sub.j,q] = t + BL(j, q, t) / (RO(j, q, t) - RI(j, q, t)),

for RO(j, q, t) > RI(j, q, t).

Case c: If the buffer rates are equal, then the buffer will stay at the same state, i.e.

(8) [tau][B.sub.j,q] = [infinity], for RO(j, q, t) = RI(j, q, t).

Case d: We consider buffer [B.sub.j,q], which is full because the downstream workstation [W.sub.q] is slowed down owing to a machine failure or blockage. The input rate of the buffer has been reduced and set equal to the output rate. If at time t the machine is repaired, or the workstation is not blocked any more, then buffer [B.sub.j,q] becomes not full immediately. We set

(9) [tau][B.sub.j,q] = t.

In a dual fashion, if buffer Bl,j is empty and the production rate of [W.sub.l] is increased at time t, then [B.sub.l,j] becomes not empty at time

(10) [tau][B.sub.l,j] = t.

A not full (not empty) event occurs only after a machine is repaired, or after an upstream buffer becomes not full (a downstream buffer becomes not empty).

9.2. Machine event scheduling

The time-to-failure of an operating machine [M.sub.i,j] is computed by the parts-to-failure, which is a variable with known distribution. The time [tau][M.sub.i,j] at which a failure occurs is given by

(11) [tau][M.sub.i,j] = t + PTF(i,j,t) / R(i,j,t)

When a machine is under repair, the time-of-next-event is computed by

(12) [tau][M.sub.i,j] = t + TTR(i, j, t),

where TTR(i, j, t) is the time-to-repair of [M.sub.i,j]. This time is a random variable of known distribution.

10. Updating and adjusting the state

The short-term variables of the system need to be updated just before the execution of a scheduled event. This is necessary because Equations 3b-5 depend on the short-term state [X.sub.S]. In addition, buffer levels and total production of the workstations determine the system performance in the intervals between successive events and must be computed. After updating the short-term variables relevant to a scheduled event, we execute it by adjusting the state variables. First we proceed with the update procedure.

10.1. Update equations

Let [t.sub.k] be the time of the previous event that caused a change to the production rate of [M.sub.i,j], and t > [t.sub.k]. Then, Equation 2 for the short-term variables takes the form:

(13) TP(i, j, t) = R(i, j, [t.sub.k])(t - [t.sub.k]) + TP(i, j, [t.sub.k]), PTE(i, j, t) = - R(i, j, [t.sub.k])(t- [t.sub.k]) + PTE(i, j, [t.sub.k]), if [M.sub.i,j] is up, TTR(i. j, t) = TTR(i,j, [t.sub.k]) - (t - [t.sub.k]), if [M.sub.i,j] is down.

Also the new level of buffer [B.sub.l,j] is

(14) BL(l, j, t) = BL(l, j, [t.sub.k]) + [RI(l, j, [t.sub.k]) - RO(l, j, [t.sub.k])](t - [t.sub.k]).

10.2. Effects of changes in machine states

The total production rate of a workstation when its operating machines produce at their nominal production rates will be called availability of the workstation. Let A(j, t) be the availability of workstation [W.sub.j] at time t. If [W.sub.j] is forced down owing to an empty buffer upstream, or because the buffers downstream are all full, then its production rate RT(j, t) is less than A(j, t). First we examine the adjustments in the long-term variables when a machine is repaired.

Case a: Let [t.sub.k - 1] be the time of the previous (k - 1)st event that caused a change in the rate of workstation [W.sub.j]. Suppose that at time [t.sub.k] machine [M.sub.i,j] is repaired. Then we mark the new machine state and adjust the availability of [W.sub.j],

(15) M(i, j, [t.sub.k]) = 1,

A(j, [t.sub.k]) = A(j, [t.sub.k - 1]) + [NR.sub.i,j].

Then, the number of parts-to-failure is found from the exponential variate generator

(16) PTE(i, j, [t.sub.k]) = - ([NR.sub.i,j] / [p.sub.i,j]) ln [xi]

where [xi] is a uniform random number on [0,1].

During the interval [[t.sub.k - 1], [t.sub.k]) workstation [W.sub.j] could be slowed down owing to blockage or starvation. This situation may be encountered even when [M.sub.i,j] is down, if the availability of [W.sub.j] is still large enough and either all downstream buffers become full, or at least one upstream buffer empties. Then, after the repair of [M.sub.i,j], the total production rate RT(j, [t.sub.k]) of [W.sub.j] remains unchanged and we set the new machine rates separately:

[lambda] = RT(j, [t.sub.k]) / A(j, [t.sub.k]) < 1,

(17) R(n,j, [t.sub.k]) =[lambda][NR.sub.nj] if [M.sub.n,j] is up, n = 1,.. .,[n.sub.j].

In this case the flow rates of the adjacent buffers are not modified.

If, however, [lambda] = 1, then the workstation is not slowed down and we set

(18a) R(i, j, [t.sub.k]) = [NR.sub.i,j],

and the total rate of [W.sub.j]

(18b) RT(i, j, [t.sub.k]) = RT(i, j, [t.sub.k - 1]) + [NR.sub.i,j].

The new output rates of the upstream buffers [B.sub.m,j ]are proportionally adjusted:

(19a) RO(m, j, [t.sub.k]) = [[delta].sub.m,j] RT(j, [t.sub.k]).

The input rates of the full downstream buffers are not modified; the increase in the production rate of [W.sub.j] is allocated to the non-full downstream buffers according to their distribution factors, that is,

(19b) [Mathematical Expression Omitted] if B(i, p, [t.sub.k]) [not equal to] 2.

Case b: Suppose now that at time [t.sub.k] machine [M.sub.i,j] breaks down. Then the new state variables become

(20) M(i, j, [t.sub.k]) = 0,}

R(i, j, [t.sub.k]) = 0,}

A(j, [t.sub.k]) = A(j, [t.sub.k - 1]) - [NR.sub.i,j],}

and the time-to-repair is

(21) TTR(i, j, [t.sub.k]) = - (1 /ri,j) ln [xi],

where [xi] is an independently sampled uniform random number on [0,1].

If the workstation is slowed down in the interval [[t.sub.k - 1], [t.sub.k]), then [lambda] < 1 and the new machine rates are computed from Equation 17 and there is no further effect on the adjacent buffers.

If the workstation is not slowed down, it is possible that some of its downstream buffers are full. Allocation of the input rates to the downstream buffers is made iteratively. Let [F.sub.j] be the set of full buffers downstream of [W.sub.j], and [G.sub.j] the set of partly full and empty buffers. The input rate of a full buffer must not exceed its output rate, whereas a non-full buffer of [G.sub.j] can accept parts at any rate by time [t.sub.k]. If after the time of failure the input rates to the full buffers become less than the output ones, then their states are set equal to l instantly and there is no further propagation upstream. Such phenomena do not occur when a machine is repaired (see case a). The rate allocation procedure (RAP) for the downstream buffers is as follows:

Step 1. Initialize. Let D = [F.sub.j][UG.sub.j] and L be the set and the total number of downstream buffers to allocate respectively. The total undisposed rate is TUR = A(j, [t.sub.k]).

Step 2. Set input rates for all buffers [B.sub.j,p] in D equal to

[Mathematical Expression Omitted]

Step 3. Trace the buffers in D:

If a buffer is full and its output rate is smaller

than its input rate, then set the input rate equal to

its output rate, reduce the total undisposed rate

by this amount, exclude this buffer from the set

D, and reduce L by one. Repeat for all elements of D.

If L has been reduced and L > 0, go to step 2. Else go to step 4.

Step 4. If L = 0 but TUR > 0, then all downstream buffers are full and the workstation is slowed down immediately.

The new output rates of the upstream buffers are computed according to Equation 12a.

10.3. Effects of changes in buffer states

We consider the four events: buffer_full, buffer_not full, buffer_empty, and buffer_not empty.

Case a: the buffer_full algorithm. If buffer [B.sub.j,q] fills at time [t.sub.k], its input rate is set equal to its output rate:

(22) B(j, q, [t.sub.k]) = 2,}

RI(j, q, [t.sub.k]) = RO(j, q, [t.sub.k]).}

This event is propagated upstream according to rule 1. If workstation [W.sub.j] can redistribute the excess of production to the non-full downstream buffers, then the input rates of these buffers are computed by using Equation 19b and there is no further effect to the workstation. If, however, all the downstream buffers are full, then the production rates are reduced:

(23) RT(j, [t.sub.k]) = RT(j, [t.sub.k - 1]) + RI(j, q, [t.sub.k]) - RI(j, q, [t.sub.k - 1]),}

[lambda] = RT(j, [t.sub.k]) / RT(J, [t.sub.k - 1]),}

R(i, j, [t.sub.k]) = [lambda] R(i, j, [t.sub.k - 1]), i = 1,...,[n.sub.j].}

This change propagates upstream. The output rates of the upstream buffers are reduced:

(24) RO(l, j, [t.sub.k]) = [[delta].sub.l,j], RT (j, [t.sub.k]).

If one buffer is full, a new buffer_full event is scheduled immediately.

Case b: the buffer_not full algorithm.

We give an example here. Suppose that the previous event (k - 1) occurred after a machine in workstation [W.sub.q] broke down and buffer [B.sub.j,q] became full. When the machine is repaired, the output flow from [B.sub.j,q] increases, and a buffer_not full event is scheduled immediately (see Equation 9). If workstation [W.sub.j] is not starved (it is possible for a workstation to be starved only when at least one of its downstream buffers is not full), its production rate is set equal to its available one, and the rates of its upstream buffers are adjusted as in Equation 24. The new input flow rates of the buffers downstream to [W.sub.j] are found using the RAP procedure described in case b of the previous paragraph.

Case c: the buffer_empty algorithm.

If buffer [B.sub.l,j] empties at time [t.sub.k] there is no effect upstream because the buffer can still accept parts, but its output rate is reduced as follows:

(25) B(l, j, [t.sub.k]) = 0,}

RO(l, j, [t.sub.k]) = RI(l, j, [t.sub.k]),}

and [lambda] = RO(l, j, [t.sub.k]) / RO(l, j, [t.sub.k - 1]) is the adjustment ratio.

The rates of workstation [W.sub.j] and its machines [M.sub.i,j] are reduced. By the assembly protocol this reduction influences the output rates of the upstream buffers [B.sub.m,j] and we have

(26) RT(j, [t.sub.k]) = [lambda] RT(j, [t.sub.k - 1]),}

R(i, j, [t.sub.k]) = [lambda] R(i, j, [t.sub.k - 1]),}

RO(m, j, [t.sub.k]) = [lambda] RO(m, j, [t.sub.k - 1]).}

If any of the other upstream buffers [B.sub.m,j] had been empty by time [t.sub.k], its state would become partially full immediately owing to the instant reduction of output rates. Again we apply the rate allocation procedure to adjust the input flows of buffers downstream to [W.sub.j].

Case d: the buffer_not empty algorithm.

In a dual fashion to case c, a buffer [B.sub.l,j] becomes not empty right after a machine of workstation [W.sub.l] is repaired (see Equation 10). Now workstation [W.sub.j] can produce faster and its rate is determined by the new input availability of [B.sub.l,j] and its assembly factor

(27) RO(l, j, [t.sub.k]) = min{RI(l, j, [t.sub.k]), A(j, [t.sub.k]) [[delta].sub.l,j]},}

RT(j, [t.sub.k]) = RO(l, j, [t.sub.k]) / [[delta].sub.l,j],}

R(i, j, [t.sub.k]) = (RT(j, [t.sub.k]) / A(j, [t.sub.k]))[NR.sub.i,j], if [M.sub.i,j] is up.}

The output rates of the buffers [B.sub.m,j] upstream to [W.sub.j] are increased:

(28) RO(m, j, [t.sub.k]) = [[delta].sub.m,j] TR)(j, [t.sub.k]).

As in previous cases we adjust the input rates of the downstream buffers with the RAP procedure.

11. The algorithm

In the two previous sections we obtained the state equations of the event-driven model. These equations were derived by considering the event-generation rules and assuming infinitesimal generation times. Both the event-driven and the conventional piece-by-piece simulation algorithms work as follows:

1. Initialize the network. Set total production period, number of buffers, their capacities, and initial levels. Enter network topology, and assembly and distribution factors for each buffer. For each machine and buffer schedule the next events.

2. Determine the next event. Find the event that occurs first. If the time of the event does not exceed the total simulation time, execute the appropriate event routine, and return to step 2, otherwise terminate the simulation.

3. Terminate the simulation. For each machine and buffer update all the short-term variables.

The event routines of the proposed model are:

(i) Machine [M.sub.i,j] is repaired:

a. Update all the machines of workstation [W.sub.j], following Equations 13 and 14.

b. Adjust the system state, as suggested by Equations 15-19.

c. Schedule the next event of the machines. If [lambda] = 1 the buffer rates will change; using Equations 14 and 6-10. update and schedule the upstream and downstream buffers.

(ii) Machine [M.sub.i,j] fails:

a. Update workstation [W.sub.j].

b. Adjust the state applying Equations 20 and 21, and the RAP.

c. Schedule the next events of [W.sub.j]. If [lambda] = 1, update and schedule the upstream and downstream buffers.

(iii) Buffer [B.sub.j,q] fills:

a. Update workstation [W.sub.j] and its downstream buffers [B.sub.j,p].

b. Adjust the state as in Equations 22-24.

c. If [lambda] = 1, i.e. at least one buffer [B.sub.j,p] is not full, use Equation 6 to compute the times of next events in buffers. Otherwise if [lambda] < 1, update the upstream buffers of [W.sub.j]. Adjust the buffer rates using Equation 24 and schedule the next events of [W.sub.j] and its upstream buffers.

(iv) Buffer [B.sub.j,q] becomes not full:

a. Update the downstream buffers of workstation [W.sub.j].

b. Using the RAP, adjust the rates of the downstream buffers.

c. If there is a further effect upstream of [W.sub.j], update its machines and buffers [B.sub.m,j] and adjust the buffer output rates as in Equation 24.

d. Schedule next events.

(v) Buffer [B.sub.l,j] empties:

a. Update workstation [W.sub.j] and all upstream and downstream buffers.

b. Following Equations 25 and 26, and the RAP, adjust the states of [W.sub.j] and its adjacent buffers.

c. Schedule next events.

(vi) Buffer [B.sub.l,j] becomes not empty.

a. Update [W.sub.j] and all upstream and downstream buffers.

b. Following Equations 27 and 28, and the RAP, adjust the system states.

c. Schedule next events.

The piece-by-piece model observes all part arrivals and departures at the machines. The corresponding event routines are:

(i) Machine [M.sub.i,j] is ready to receive parts:

a. Remove items from all the upstream buffers according to their assembly factors.

a1. If any of these buffers happens to be full, look for blocked upstream machines; if there are blocked machines, then schedule a departure event instantly at the machine with the highest priority.

a2. Also, if any of these buffers happens to be exhausted, then [M.sub.i,j] becomes the starved machine with the lowest priority, and its next event time becomes [tau][M.sub.i,j] = [infinity].

b. If [M.sub.i,j] is not starved, then schedule the production of a part after a processing cycle; this cycle may include downtimes due to the occurrence of one or more failures.

(ii) Machine [M.sub.i,j] is ready to produce one part:

a. Find the non-full downstream buffer that will accept the part.

a1. If all the downstream buffers are full, then [M.sub.i,j] becomes the blocked machine with the lowest priority, and we set [tau][M.sub.i,j] = [infinity].

a2. If the destination buffer happens to be empty, then look for starved downstream machines. If there are starved machines, then schedule a part arrival at the machine with the highest priority instantly.

b. If [M.sub.i,j] is not blocked, then the machine becomes ready to receive parts instantly.

12. Model validation

In this section several results are reported, comparing the accuracy of our event-driven (e-d) algorithm and its speed to those of a piece-by-piece (p-p) simulator described in the previous paragraph. Both algorithms were written in FORTRAN 77 code for a HP 9000 Workstation. About 1000 lines of source code were written for the e-d model and 500 for the pep model. To make the comparisons under the same experimental conditions, identical random numbers were used in both algorithms. We define the computing index by the ratio:

CI = CPU time of the p/p model/

CPU time of the e-d model

First we study a production line with n workstations and n - 1 intermediate buffers with capacity 20 workparts. Each workstation consists of three identical machines, each having production, failure and repair rates respectively 10, 0.1, and 1. The initial buffer level is 10 workparts and the simulation horizon 3000 time units.

In Fig. 3 we demonstrate the increase of the relative speed of our model with increasing length n. The computing index is 13 for a line with two workstations and an intermediate buffer, and 22.7 for a line with n = 50 workstations. The computing ratio increases with system size. The model provides very accurate estimates for the average throughput rate with errors never exceeding 0.08%. The maximum errors in average buffer levels were less than 2%.

The next production system studied is sketched in Fig. 4. The numbers in boxes indicate the number of parallel machines of each workstation and the numbers in circles correspond to the buffer capacities. The assembly factors are read along the arcs in front of the assembly workstations. Also, the distribution factors at the flow splitting workstations are equal. That is, if buffers [B.sub.7,10] and [B.sub.7,11] are not full, then each one of these will accept half of the production of [W.sub.7]. The production and repair rates are 10 and I respectively, and failure rates take values on the interval [0.08,1]. The computing index is plotted in Fig. 5. The relative speed increases almost exponentially with the machine reliability. When [p.sub.i,j] = I the computing index of our model is 1.88. In this case, although still more efficient, our algorithm updates and schedules the state variables almost as frequently as the piece-by-piece simulator, because machines fail and are repaired very often. However, actual production systems are relatively reliable and can be efficiently analyzed by our model. The discrepancy in the estimation of average buffer levels was less than 0.8 workparts, although, when their levels were low, this gave a maximum error of 20% which was observed in one of the boundary buffers, i.e. those connected with the external nodes [B.sub.I] or [B.sub.O]. The errors for the other buffers were less than 5%. The errors of throughput rates for the various products of the system were less than 1.77%.

Next we considered the network of Fig. 6. It consists of single-machine workstations with production and repair rates 10 and 1 respectively. The assembly factors are all 1, and the distribution factors at the flow-splitting nodes are equal. The buffer capacities were all set equal to 10. As shown in Fig. 7, CI degrades very quickly with increasing failure rate. For the buffer B of Fig. 6 the discrepancy of the average level estimate was about 0.5 workparts (16%); for the other buffers the errors were less than 5%. The estimation error of the throughput rate was less than 0.15% in all cases.

In general our algorithm is much faster than piece-by-piece simulators and very accurate for a wide variety of production system topologies and parameters.

13. Conclusion

In this paper we have developed a hybrid simulation/ analytic model for multiproduct production systems with assembly operations. The mathematical representation of the system by a set of continuous-time algebraic equations and a set of difference equations revealed the decomposability property of the model, and directly led to an event-driven algorithm with reduced computation cost. Average throughput rate accuracy is quite remarkable. The model can be incorporated into design algorithms of production systems to compute their performance sensitivities efficiently with respect to various decision parameters. The model can also simulate on-line-controlled production networks. At a time when a specific production policy is decided, the routing protocol and the throughput rate of any product can be modified by appropriately adjusting the production rates of machines.

[Figures 1 to 7 ILLUSTRATION OMITTED]

Acknowledgement

This research was supported in part by two grants from the Ministry of Education in Athens, Greece.

References

Bocharov, P. P. (1987) Approximate method of design of open nonexponential queueing networks of finite capacity with losses or blocking. Automation and Remote Control, 44-53.

D'Angelo, H.. Caramanis, M., Finger, S., Mavretic, A., Phillis, Y. and Ramsden, E. (1988) Event-driven model of unreliable production lines with storage. International Journal of Production Research 26, 1173-1182.

Ho, Y. C., Eyler, M. A. and Chien, T. T. (1979) A gradient technique for general buffer storage design in a production line. International Journal of Production Research 17, 557-580.

Jackson, J. R. (1957) Networks of waiting lines. Operations Research 5, 518-521.

Kouikoglou, V. S. and Phillis, Y. A. (1991) An exact discrete-event model and control policies for production lines with buffers. IEEE Transactions on Automatic Control AC-36, 515-527.

Kuehn, P. J. (1979) Approximate analysis of general queueing networks by decomposition. IEEE Transactions on Communications COM-27, 113-125.

Takahashi, Y., Miyahara, H. and Hasegawa, T. (1980) An approximation method for open restricted queueing networks. Operations Research 28, 594-602.

Biographies

Yannis A. Phillis received a diploma in electrical and mechanical engineering from the National Technical University of Athens, Greece, in 1973, and M.S. and Ph.D. degrees in engineering systems from the University of California, Los Angeles, in 1978 and 198O, respectively. From 1980 to 1986 he was with Boston University, Boston, Massachusetts. Since 1986, he has been with the Department of Production Engineering and Management, Technical University of Crete, Chania, Greece, where he is a Professor and President. In 1986 he received a `Professor of the Year Award' at Boston University. He is also a known writer in Greece. His research interests are in schochastic control, discrete-event systems, and applications in manufacturing networks.

Vassilis S. Kouikoglou was born in Nafplion, Greece, in 1961. He received a diploma in electrical engineering from the National Technical University of Athens, Greece, in 1985 and a Ph.D. degree in production engineering and management from the Technical University of Crete, Chania, Greece, in 1989. Current]y, he is an Assistant Professor in the Department of Production Engineering and Management, Technical University of Crete. His research interests are in modeling discrete-event systems, optimization, and queueing systems.