Scheduling manufacturing systems with work-in-process inventory control: single-part-type...

1. Introduction

Manufacturing systems are complex. Many types of machines, workers, and parts are often involved. The variety of random events makes the scheduling of manufacturing systems difficult. For example, in a semiconductor fabrication factory, dozens of part types are produced simultaneously

by hundreds of workers on dozens of machines. Each part type follows a predefined process that consists of hundreds of operations. Machines are subject to random failures, and need setup changes for different part types. Maintenance and rework must be considered. Workers are absent at random. These factors can result in long throughput time (also called cycle time or lead time), large work-in-process (WIP) inventory, and significant tardiness.

The WIP inventory in a manufacturing system affects significantly the throughput time and the tardiness. To improve the efficiency of production, we would like to reduce inventory, throughput time, and tardiness simultaneously. In this paper, we study how WIP can be used with a sophisticated production control policy. One key question is: what is the minimal necessary WIP and how should it be allocated in a manufacturing system to make the production effective?

1.1. Previous research

There is a large body of literature in production scheduling. Much of it is surveyed in Graves (1981). Many of the articles before the early 1980's are based on combinatorial optimization/integer programming or mixed-integer methods (Afentakis et al., 1984; Lageweg et al., 1977, 1978; Newson, 1975a,b; Papadimitriou and Kannelakis, 1980; Wagner and Whitin, 1958). Some other works are based on queuing network models (Chen et al., 1988; Glassey and Resende, 1988; Harrison and Wein, 1990; Wein, 1988, 1990).

Because many kinds of machines, workers, part types, and operations are involved in real production systems, hierarchical structures have been proposed for production control in order to reduce the problem size and complexity (Bitran et al., 1981; Dempster et al., 1981; Gershwin, 1989; Graves, 1982; Hax and Meal, 1975). The goal is to replace one large problem by a set of many small ones, because each of the latter is invariably easier to solve. Tactical planning models are also developed to analyse job shops (Graves, 1986; Fine and Graves, 1990). In those models, work flow is aggregated during each period so that the complexity is reduced. Even so, the enormous size of production systems and the variety of random events associated with manufacturing procedures make traditional optimization methods, in many cases, inadequate or inappropriate for production scheduling, especially in real time.

Since the early 1980's, production flow models have been developed to reduce further the complexity of scheduling problems. In those formulations, part movement in a production system is treated as a continuous flow so that the dimension of the model is reduced dramatically. Furthermore, the system dynamics of the production flow models are in a form that is appropriate for control theory and techniques.

Using Rishel's (1975) methodology, Kimemia and Gershwin (1983) investigated the optimal flow controller's structure and determined that it is a hedging-point, feedback control policy. Tsitsiklis (1982) proved the convexity of the value function that satisfies the Hamilton-Jacobi-Bellman equations and determines the optimal controller. Gershwin et al. (1985) proposed a heuristic approximation of the value function. Akella and Kumar (1986) solved the Hamilton-Jacobi-Bellman equation analytically to obtain the optimal value function for a simple single-part-type, one-machine system. Van Ryzin (1987) studied the delay of the production flow in a buffer and obtained a numerical solution for a single-part-type, two-machine system. In short, much effort has been directed to the development of the production flow control models, for both analytical solutions and approximation methods (see also Akella et al., 1984; Maimon and Choong, 1987; Maimon and Gershwin, 1988; Sharifnia, 1988; Eleftheriu, 1989).

Work-in-process (WIP) inventory plays a very important role in production scheduling. It has drawn a great deal of attention from researchers. Conway et al. (1988) studied the effects of WIP in serial production lines. Burman et al. (1986) investigated the relation between the WIP level and the system performance of integrated circuit manufacturing lines. Zeghmi (1985) studied inventory buffers in a production line. Because of the complex way in which WIP interacts with all the random events, WIP control in a dynamic environment, such as production systems scheduled in real time, is still not well solved and understood.

1.2. Outline of the paper

In this paper a real-time feedback control algorithm is developed for scheduling manufacturing systems. We consider three classes of activities: operations, failures and repairs, and starvation and blockage. The scheduling objectives are to keep the actual production close to the demand and to keep the WIP inventory level low and cycle time short. In the algorithm we combine strategic decision-making, tactical planning, and operational scheduling into a hierarchy to respond to events with different frequencies. The major contribution of this paper is that we explicitly introduce buffer sizes, average buffer levels, and starvation and blockage fractions as control parameters for long-term capacity planning as well as for short-term scheduling. For long-term capacity planning, the buffer sizes are minimized such that the system has enough capacity to achieve the demand. A desirable system state, called the hedging point, is also determined at this stage, which will be used to regulate the production rates. For short-term scheduling, the controller recalculates the production rates in real time whenever a machine fails or is starved or blocked. The system always attempts to move toward the hedging point and tries to stay there if possible. Then the production rates are used to determine the loading times for individual parts, according to a simple staircase strategy. The production scheduling algorithm is evaluated by using computer simulations for a variety of cases. Compared with a transfer line operation policy (i.e., performing an operation whenever possible), a significant improvement in system performance is observed.

In the next section we study the two-machine, single-part-type systems. We investigate the buffer effects and the relation between WIP and starvation or blockage. The following steps are carried out in the development of the real-time production control algorithm.

Step 1. We formulate the scheduling problem as a dynamic program. We observe that the optimal solution of the DP formulation is characterized by the partition of the state space.

Step 2. We develop a heuristic partition of the state space such that the system behaves in a desirable manner, which will be explained in Section 2.2. We demonstrate on intuitive grounds that the boundaries in the state space should be orthogonal to the axes.

Step 3. The location of the boundaries and of the buffer sizes are determined by using a nonlinear program.

The results of the simple case are extended to N-machine, single-part-type systems in Section 3. Each part travels in the same fixed sequence: machine 1, buffer 1, . . ., buffer N - 1, machine N. A machine in the middle of the production line can be either starved or blocked. In Section 4, the system performance estimated by the real-time algorithm is compared with that of a tandem queueing model. Concluding remarks are in Section 5.

2. Two-machine, single-part-type systems

In this section we study two-machine, one-buffer, single-part-type systems (see Fig. 1). Both machines are subject to random failures. For machine i (i = 1,2), the time to fail and time to repair are modeled by exponentially distributed random variables with means 1/[[Rho].sub.i] and 1/[r.sub.i], respectively. One part type is produced. Each part needs an operation with processing time [[Gamma].sub.1] on machine 1 and then an operation with time [[Gamma].sub.2] on machine 2. The processing times are assumed to be deterministic and small compared with the time to fail and the time to repair. Therefore the movement of parts in the system is modeled as a continuous flow. A buffer is located between the two machines. We assume that machine 1 is never starved and machine 2 is never blocked. In addition, we assume that a machine can be starved or blocked only when it is operational.

[Figure 1 ILLUSTRATION OMITTED]

Denote by [x.sub.i] the production surplus. The rate of change of [x.sub.i] is given by

(1) [x.sub.i] = [u.sub.i] - d (i = 1,2),

where [u.sub.i] is the production rate of the ith operation and d is the demand rate.

Denote by b the buffer level, i.e., the number of parts in the buffer. Because we represent part movement as continuous flow, the buffer level b is a real number, not restricted to integers. The rate of change of the buffer level is governed by

(2) b = [u.sub.i] - [u.sub.2].

The buffer level, b, is a function of the surplus x, which can be determined by

(3) b(t) = [x.sub.i](t) - [x.sub.2](t),

if the initial buffer level b([t.sub.0]) is

(4) b([t.sub.0]) = [x.sub.1]([t.sub.0]) - [x.sub.2]([t.sub.0]).

The buffer level is also subject to the constraint

(5) B [is greater than or equal to] b [is greater than or equal to] 0,

where B is the buffer size.

To maximize the efficiency of the production system, we wish to minimize the WIP level and cycle time, and to meet demand on time. Suppose that g(x,b) is a convex function that penalizes x(t) and b(t) for not being close to zero. In practice, different systems may require different penalty functions. Furthermore, in the development of the scheduling algorithm below we shall largely ignore exact form of the penalty function. Therefore in the following we formulate a dynamic program for production control without specifying precisely the g function. The solution of this program has a partition structure regardless of g (as long as g is convex, has a minimum at x = 0, and goes to [infinity] as |x| [right arrow] [infinity]).

Given an initial surplus state x([t.sub.0]) and machine state [Alpha]([t.sub.0]), the production flow rate control problem can be formulated as a dynamic optimization problem

(6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(7) [u.sub.1] [is less than or equal to] [U.sub.1] [[Alpha].sub.1]; [u.sub.2] [is less than or equal to] [U.sub.2] [[Alpha].sub.2]; [u.sub.1] [is greater than or equal to] 0, [u.sub.2] [is greater than or equal to] 0;

where [[Alpha].sub.i] is the machine state, which is a binary random variable at any given time instant t. It is 0 if machine i is down and 1 otherwise. [U.sub.i] is the maximum rate for operation i (i = 1, 2). In this case, [U.sub.i] = 1/[[Tau].sub.i].

The dynamic program (6) is to control a stochastic process with jump Markov disturbances. This class of problem has been studied in previous work (see Rishel, 1975; Tsitsiklis, 1982; Kimemia and Gershwin, 1983; Akella and Kumar, 1986; Van Ryzin, 1987). Although there is no general technique available to solve the problem analytically, previous results have shown the characteristics of the optimal solution. The most important finding is that the optimal control is a switching policy. That is, the state space is partitioned into mutually exclusive regions and production rates are constant within each region. Control actions take place when the system moves from one region to another. This result will serve as a guideline for the design of a production control scheme in this paper.

The dynamic programming formulation implies that the optimal production policy is based on a partition of the state space. However, it is difficult to determine the optimal partition from dynamic programming considerations. Consequently, we develop a heuristic partition in this section; first we demonstrate on intuitive grounds that the boundaries should be orthogonal to the axes, and then we locate them by using a nonlinear program.

If J were easy to calculate, we would derive the partition from it. Because it is not, we use only dynamic programming to tell us that the control policy is based on a partition of the space. We derive a partition on heuristic grounds, and do not try to construct J.

2.1. Partition of the state space

If the optimal value function J(x, [Alpha], t) is differentiable with respect to x and t, the optimal production rate u can be determined by solving the linear programming problem (see Rishel 1975; Kimemia and Gershwin, 1983; Bai, 1991b)

(8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to (7), and the system dynamics are governed by (1), (2), and (5). By inspecting the linear program above, we can make the following observations.

Observation 1. The right-hand side of the inequalities in the constraint set (7) are random parameters. Therefore the shape and size of the capacity constraint set change randomly with the machine status [Alpha](t).

Observation 2. The linear program (8) represents a real-time feedback control law because the LP is determined by x and [Alpha]. When the production surplus x(t) and machine state [Alpha](t) are fed back from the shop floor at time t, a new production rate, u(t), is generated by solving the linear program.

Observation 3. The objective function is linear in the production rate u. For each instance of [Alpha], the constraint set is a convex polyhedron in u-space. Therefore the optimal production policy u(t) takes on values at the extreme points of the constraint set. The coefficients of u in the objective are functions of the production surplus x. Therefore for each machine state [Alpha], an optimal production policy divides the x-space into a set of regions in which the production rate is constant (see Fig. 2). Each region corresponds to an extreme point of the constraint set. These regions cover the whole x-space except where the gradient of J is zero or where the subgradient of J contains 0, if J has discontinuous derivatives. (In the following, to simplify the discussion, we assume that J has continuous derivatives.) In that situation, the controller (8) becomes singular and a unique optimal solution may not exist. A set of conditional constraints is introduced below to override the singularity.

[Figure 2 ILLUSTRATION OMITTED]

In the case that both machines are operational ([[Alpha].sub.1] = 1 and [[Alpha].sub.2] = 1) Fig. 2 illustrates the regions in x-space. The two straight lines are the zero-buffer and full-buffer boundaries because the buffer level is a linear function of x. The feasible region of x-space lies between the zero-buffer boundary (b = 0) and the full-buffer boundary (b = B). The other two curves are the coefficient boundaries, which are the sets of points in which one of the coefficients of the objective function in (8) is zero ([Delta]J/[Delta][x.sub.i] = 0). The feasible area of the x-space is divided into four mutually exclusive regions that correspond to the four extreme points of the constraint set. In each region, the production rate u is constant. The intersection of the coefficient boundaries is called the hedging point (HP), which is the desirable operating state of the system. The feedback controller (8) always attempts to drive the system to the hedging point, and to keep it there.

Observation 4. The coefficient boundaries are attractive. (See Gershwin et al. (1985) for a similar situation.) That is, when the system state reaches a boundary, it would cross the boundary back and forth while it moves along on a boundary towards the hedging point if u were determined only by (8). This phenomena is referred to as chattering on a boundary, and would severely affect the efficiency of the production control algorithm. To avoid chattering, a set of conditional constraints is used to guide the system to move to the hedging point when production surplus x reaches a coefficient boundary. The conditional constraints are stated in Section 2.4.

Observation 5. Because the production rates are constant in every region in x-space, the production rates do not have to be calculated at every time instant. They need only to be computed when machine state [Alpha] changes or when production surplus x(t) reaches a boundary. In fact, once the boundaries are known, it is easy to calculate when x(t) will reach a boundary (if a change in [Alpha]does not occur first) and which region x(t) will enter next.

We have completed the first step for the development of the production controller listed at the end of Section 1. We have discussed some properties of the production controller. The linear program (8) always attempts to drive the system toward hedging point along the steepest descent direction of the cost function, over a feasible region in x-space. We have observed that the controller is characterized by the shape and position of the coefficient boundaries in x-space and the buffer size. In other words the production policy is completely determined by the partition of x-space. In the sequel of this section, we construct a heuristic partition of the x-space such that the system behavior is satisfactory.

2.2. System behavior specification

In the previous section we observed that the optimal control is a switching policy that is characterized by the partition of the x-space. The optimal value function J is difficult to find. However, we care only about the shapes and positions of the boundaries and need not determine J. Therefore we work with the boundaries in x-space directly and do not construct the J function. In the following, we specify a list of requirements regarding system behavior, which serves as a guide for shaping the boundary in x-space in the subsequent subsections. The desirable system behavior specifications are:

(a) when machine 1 fails, keep machine 2 producing without changing its production plan until the buffer is empty;

(b) when machine 2 fails, keep machine 1 producing without changing its production plan until the buffer is full.

These behavior requirements derive from the considerations of spatial decomposition. For example, for a large factory, when a single machine fails, we do not wish to change the production plan of the entire factory unless we have to. Consequently we would like to separate the machines as much as possible to reduce the effects of machine failures. This consideration is essential for dividing a system into several subsystems in a hierarchical structure. For instance, machines in a factory can be grouped into cells and separated by storage areas. The behavior requirements ensure that the machine status in a cell would not affect adjacent cells, provided the storage buffers are not empty or full. Therefore short-term disruptions could be absorbed locally and the production would have to be rescheduled only for a single cell. In this case, the state of machine 1 does not affect machine 2 if the buffer is not empty, and the state of machine 2 does not affect machine 1 if the buffer is not full. In later sections we shall see more precisely how these requirements reduce the complexity of the production control problem.

2.3. The desirable boundary shape in x-space

From the discussion in the previous sections we see that the production controller (8) is characterized by the shape and position of the coefficient and the buffer-size boundaries in x-space. In this section: (1) we determine the proper boundary shape such that the system behaves as specified in Section 2.2; and (2) we find a linear program of the form

min{[[Alpha].sub.1] (x, [Alpha])[u.sub.1] + [a.sub.2](x, [Alpha])[u.sub.2]} subject to (7)

that produces the same control actions as these boundaries. It is important to note that the boundaries are heuristic, and not intended as close approximations of the optimal boundaries. Therefore this linear program is not intended to be a close approximation to (8), and the coefficients [a.sub.i] are therefore not intended to be close approximations of [Delta]J/[Delta][x.sub.i]. Rather, we select them solely to produce the desired behavior outlined above.

The behavior specification of the previous section says that the choice of [u.sub.1] should be independent of [[Alpha].sub.2], and that the choice of [u.sub.2] should be independent of or [[Alpha].sub.1]. Furthermore, the constraints of the linear program require that if [[Alpha].sub.1] = 0, then [u.sub.1] = 0. Therefore we can assume that [a.sub.i] is independent of [Alpha].

Suppose that the system is not at the hedging point. In particular, assume that x is in regions (1) or (4) of Fig. 3. In that case, [u.sub.2] = 0, so the coefficient [a.sub.2](x) must be positive. Similarly, if x is in regions (2) or (3), [u.sub.2] is at its maximum value, so [a.sub.2](x) must be negative. If we assume that [a.sub.2](x) is continuous, then [a.sub.2](x) = 0 on the boundary between regions (1) and (2) and on the boundary between regions (3) and (4). By the same reasoning, [a.sub.1](x) = 0 on the boundary between legions (1) and (4) and on the boundary between regions (2) and (3). In particular, [a.sub.1](z) = [a.sub.2](z) = 0, where z is the hedging point. (Note that we have not yet specified u on the boundaries.)

[Figure 3 ILLUSTRATION OMITTED]

Suppose the system has reached the hedging point and both machines are operational. The production rate decision is then [u.sub.1] = [u.sub.2] = d, so the system stays at the hedging point indefinitely. Then, suppose that at time t after the system reaches the hedging point, machine 1 fails ([[Alpha].sub.1] = 0) and machine 2 is still operational ([[Alpha].sub.2] = 1). According to the behavior requirement (a) of Section 2.2 and the capacity constraints of (8), the production decision should be [u.sub.1] = 0 and [u.sub.2] = d until the buffer is empty. Before the buffer is empty, [x.sub.1] decreases at rate d and [x.sub.2] is constant. Therefore the system state ([x.sub.1], [x.sub.2]) moves along a horizontal line towards the zero-buffer boundary ([x.sub.1] = [x.sub.2]). This line must be the boundary between regions (3) and (4). If it were not, [u.sub.2] would be 0 or [U.sub.2], and [x.sub.2] would not be constant.

Until the system reaches the zero-buffer boundary, since [u.sub.2] = d and [x.sub.2] does not change, the coefficient of [u.sub.2] remains zero ([a.sub.2](x) = 0). Therefore the coefficient boundary defined by [a.sub.2](x) = 0 must be straight and horizontal, and must go through the hedging point. Similarly, the other boundary must be straight and vertical, and must go through the hedging point z. Thus, regions (1) and (2) are where [x.sub.1] [is greater than] [z.sub.1], and regions (3) and (4) are where [x.sub.1] [is less than] [z.sub.1]. Similarly, regions (1) and (4) are where [x.sub.2] [is greater than] [z.sub.2], and regions (2) and (3) are where [x.sub.2] [is less than] [z.sub.2].

We have established that [a.sub.i][is greater than] 0 when [x.sub.i][is greater than] [z.sub.i] and [a.sub.i] [is less than] 0 when [x.sub.i] [is less than] [z.sub.i]. The simplest expression that meets this condition is

[a.sub.i] = [x.sub.i] - [z.sub.i].

The desirable boundary shape in x-space is illustrated in Fig. 3. Regions (1), (2), (4) and the dashed boundaries are transient states. This is because after the system reaches the hedging point, it will never go above the horizontal boundary or to the right of the vertical boundary. On the horizontal boundary, [x.sub.2] is constant ([u.sub.2] = d). On the vertical boundary, [x.sub.1] is constant ([u.sub.1] = d). When the buffer is empty, machine 2 cannot produce faster than machine 1. When the buffer is full, machine 1 cannot produce faster than machine 2.

2.4. The conditional constraints

To avoid chattering on the coefficient boundaries and to respect the buffer constraints, we need to add a set of conditional constraints to the linear program (8).

Denote by ([z.sub.1], [z.sub.2]) the hedging point, which is the desirable value of ([x.sub.1], [x.sub.2]). The components of the hedging point are unknown and will be determined in Section 2.6.2. The conditional constraint set is

(9) if [x.sub.1] = [z.sub.1] and B [is greater than] b [is greater than] 0, then [u.sub.1] = d[[Alpha].sub.1]; if [x.sub.2] = [z.sub.2] and B [is greater than] b [is greater than] 0, [u.sub.2] = d[[Alpha].sub.2]; if b = 0, [u.sub.1] [is greater than or equal to] [u.sub.2]; if b = B, [u.sub.1] [is less than or equal to] [u.sub.2];

which says that when the system reaches the vertical boundary, the production rate of machine 1, [u.sub.1], should be equal to the demand d. When the system reaches the horizontal boundary, the production rate of machine 2, [u.sub.2], should be equal to the demand. When the buffer is empty, machine 2 cannot produce faster than machine 1. When the buffer is full, machine 1 cannot produce faster than machine 2.

2.5. The linear program for real-time feedback control

The following linear program ensures that the coefficient boundaries in x-space are horizontal and vertical and go through the hedging point:

(10) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to (7) and (9). The system dynamics are governed by (1), (2), and (5). Note that the coefficient of [u.sub.1] is not a function of [x.sub.2]. Therefore, if the buffer is neither empty nor full, (10) indicates that the production flow rate, [u.sub.1], is independent of the flow rate [u.sub.2], the machine state [[Alpha].sub.2], and the surplus state [x.sub.2]. The same observation can be made for [u.sub.2]. Coupling occurs only when the buffer is either empty or full as specified in the conditional constraints. Because the production rates are constant in every region as well as on every boundary in x-space, (10) yields a piecewise constant solution u(t).

When the system reaches the hedging point, the scheduler generates the same behavior as in a Kanban system. That is, the buffer level stays constant. Moreover, if the system drifts away from the hedging point owing to random events such as machine failures or starvation or blockage, new policies are generated such that the system recovers as soon as possible.

We have now completed step 2 in the list at the end of Section 1. We have constructed the partition of the state space and ended with the feedback controller in the form of linear program. In the linear program there are three parameters that we need to determine: the components of the hedging point ([z.sub.1], [z.sub.2]) and the buffer size B; this will be done in the following subsections.

2.6. Control parameter estimation

In this subsection we estimate the unknown parameters of the feedback control linear program (10). To do so we develop a method to formulate an approximate relationship among starvation, blockage, buffer hedging level, buffer hedging space, machine parameters, and demand. Consequently, a nonlinear optimization problem is set up to determine the buffer size. The hedging point is determined such that the average final product inventory is minimized. In the following we present the formulations without derivations. Interested readers are referred to Appendix A or Bai (1991b).

2.6.1. Starvation and blockage

Define [z.sup.b] to be the buffer hedging level (see Fig. 3). It is the buffer level when the system reaches the hedging point, which satisfies

(11) [z.sup.b] = [z.sub.1] - [z.sub.2]

Define [z.sup.s] to be the buffer hedging space. It is the room left for more parts in the buffer when the system reaches the hedging point, which satisfies

(12) [z.sup.s] = B - [z.sup.b]

Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be the fractions of time that machine i (i = 1, 2) is starved or blocked. A machine can be starved or blocked only when it is operational. The starvation and blockage fractions and buffer hedging level and space can be determined as follows (see Appendix A).

(13) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [D.sub.i] is the isolated capacity of machine i, given by [D.sub.i] = [[r.sub.i]/ ([r.sub.i] + [p.sub.i])] [U.sub.i] (i = 1, 2)

2.6.2. The hedging point

The components of the hedging point can be determined recursively:

(14) [z.sub.2] = [[Delta].sub.2]; [z.sub.1] = [z.sup.b] + [z.sub.2];

where [[Delta].sub.i] (i = 1, 2) is the average surplus loss at machine i, which is approximately given (see Appendix A) by

(15) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

2.6.3. The buffer size and average buffer level

The buffer size is defined in (12) as

B = [z.sup.b] + [z.sup.s],

The average buffer level is different from the buffer hedging level. Let [bar]b be the average buffer level, which can be obtained by taking the time average of (3):

(16) [bar]b = [[bar]x.sub.1] - [[bar]x.sub.2]

The relation between the hedging buffer level and the average buffer level is given by

[bar]b = [z.sup.b] + ([[Delta].sub.2] - [[Delta].sub.1])

In the formulation above, the average buffer level should be truncated to ensure the satisfaction of (5). If the estimates of the surplus losses are used,

(17) [bar]b = min{max (0, [z.sup.b] + [[Delta].sub.2] - [[Delta].sub.1]), B}.

2.6.4. The average cycle time and WIP level

The cycle time includes the waiting times in buffers and processing times on machines. Each part travels through the system with an average rate d. The average waiting time in the buffer is [bar]b/d. The total processing time is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Denote by [bar]CT the average cycle time that a part spends in the system:

(18) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Denote by W the average WIP level in the system. By Little's law, we have

(19) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We have now constructed the real-time scheduler for the two-machine, single-part-type systems. In Section 2.5 the feedback controller was established as a linear programming problem with three unknown parameters, namely, the buffer size B and the components of the hedging point, [z.sub.1] and [z.sub.2] The components of the hedging point were determined in Section 2.6.2. The buffer size was determined in Section 2.6.3.

2.7. The production scheduling algorithm and the hierarchical policy

In this section, we summarize the steps of the algorithm of the production scheduler and describe the hierarchical structure for implementation.

Step 1. Collect the input data set, which consists of the failure rate [p.sub.i], the repair rate [r.sub.i], and the processing time [[Tau.sub.i], for machine i (i = 1, 2), and the demand [d.sub.i].

Step 2. Calculate the buffer hedging level [z.sup.b] and hedging space [z.sup.s], the starvation fraction [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], and the blockage fraction [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] according to (13). Then calculate the buffer size B by summing the buffer hedging level and hedging space.

Step 3. Calculate the components [z.sub.1] and [z.sub.2] of the hedging point according to (14).

Step 4. Using the feedback information of production surplus [x.sub.i](t) and machine state [[Alpha].sub.1](t) (i = 1, 2), calculate the production rates, [u.sub.i](t) (i = 1, 2), in real time by solving the linear program (10).

Step 5. The loading times for each machine are determined by a heuristic policy such as the staircase strategy (Akella et al., 1984). That is, whenever the actual cumulative production is less than the integral of the production rate, load a part into the machine.

Step 6. If the demand or any one of the machine parameters changes, go to step 2. If the machine state changes the production surplus reaches a boundary, go to step 4.

This production scheduling algorithm can be divided into a three-level hierarchy (Kimemia and Gershwin, 1983). At the top level of the hierarchy is a control parameter estimator. It calculates the buffer size and the hedging point given the demand and the machine parameters (as we state in steps 1, 2, and 3 of the algorithm). A real-time production rate generator fills the middle level of the hierarchy. It calculates the production rates in real time, while the machine states and the production surplus are fed back from the shop floor (step 4 of the algorithm). The bottom level of the hierarchy is a loading instruction generator. It determines the loading times for each machine by using the staircase strategy (step 5 of the algorithm).

2.8. Simulation examples

To verify the algorithm we have run computer simulations of production control using the real-time scheduler for a variety of cases. In this section we demonstrate several simulation examples of the production control algorithm.

We consider three different cases for the system described at the beginning of this section. In the first case, the system is balanced in the sense that both machines have the same failure and repair rate and processing times. In the second case, however, the second machine is the bottleneck, which is slower than the first on the average. Conversely, the first machine is the bottleneck in case 3. The machine parameters and processing times are chosen as follows:

case 1: [r.sub.1] = 0.5,   [p.sub.1] = 0.1,     [[Tau].sub.1] = 0.5;
        [r.sub.2] = 0.5,   [p.sub.2] = 0.1,     [[Tau].sub.2] = 0.5;

case 2: [r.sub.1] = 0.05,  [p.sub.1] = 0.01,    [[Tau].sub.1] = 0.4;
        [r.sub.2] = 0.5,   [p.sub.2] = O.1,     [[Tau].sub.2] = 0.8;

case 3: [r.sub.1] = 0.5,   [p.sub.1] = 0.1,     [[Tau].sub.1] = 0.8;
        [r.sub.2] = 0.05,  [p.sub.2] = 0.01,    [[Tau].sub.2] = 0.4;

where the unit of r and p is 1/day. The unit of [Tau] is a day. The unit of parts is a lot.

For each case, the buffer size and hedging point have been calculated for different demand values and are listed in Tables 1-3, where the buffer size is rounded off to an integer.

[TABULAR DATA 1-3 NOT REPRODUCIBLE IN ASCII]

The simulation program that we use is called HIERCSIM; it was developed by Darakananda (1989). Fig. 4 illustrates the cumulative production results of a single simulation run for d = 1.6 in case 1. The straight line is the cumulative demand. The upper curve is the cumulative input of the raw parts at machine 1. The lower curve is the cumulative output of the final products at machine 2. The dashed lines are the middle-level results, which are the integrals of the flow rates. The staircase-like graphs are the bottom-level results, which are the actual counts of cumulative production. It is almost impossible to tell the difference between the middle (level 2) and bottom-level (level 3) results. Fig. 5 illustrates buffer level against time. The buffer level consists of the parts in the buffer and in the work space of machine 1. (For this reason, the buffer level is sometimes 6 while the buffer size B is 5.) The dashed line is the middle-level (level 2) result. That is, it is the buffer level b(t), which is governed by (3). The solid line is the bottom-level (level 3) result, which is the actual count of the parts in the buffer. The actual count and b(t) rarely differ by more than 1.

[Figures 4 and 5 ILLUSTRATION OMITTED]

To verify the algorithm we ran simulations for the three cases above. For a given demand in each case, five simulation runs were made, using different seeds for the random number generator. The starvation and blockage fraction of time, the average buffer level, and the average WIP inventory were averaged over the five runs and are listed in Tables 1-3. The in-process inventory W was estimated by using (13), (17), and (19), and the results are also presented in Tables 1-3. We see that the simulation results agree well with the estimated values in general. However, the following patterns are recognizable. First, the real-time algorithm provides better prediction on system performance (starvation, blockage, buffer level, WIP) for balanced production lines (case 1) than for unbalanced ones (cases 2 and 3). Secondly, simulations agree with the estimates better when demand is high, relative to the system capacity, than when demand is very low.

3. N-machine, single-part-type systems

In this section we study single-part-type systems consisting of N machines and N - 1 buffers. For machine i (i = 1, 2, ..., N), the time to fail and time to repair are assumed to be exponentially distributed random variables with means 1/[p.sub.i], and 1/[r.sub.i] respectively. One part type is produced. Each part needs an operation with processing time [[Tau].sub.i] on machine i. The parts travel in a fixed sequence: machine 1, buffer 1, machine 2, buffer 2, ..., machine N. The buffers are located between machines. We assume that machine 1 is never starved and machine N is never blocked. A machine can be starved or blocked only when it is operational.

In this case, a machine in the middle of the production line can be either starved or blocked. The relations among machines are more complex than the previous case, because a machine failure can starve or block more than one machine. The technique developed in the previous section is extended to deal with the serial production line.

Given the system as described above, the production control problem is formulated as follows:

(20) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to

(21) [u.sub.i] [is less than or equal to] [U.sub.i][[Alpha].sub.i] (i = 1, 2,..., N);

[u.sub.i] [is greater than or equal to] 0 (i = 1, 2,..., N);

where the system dynamics and buffer constraints are

(22) [x.sub.i] = [u.sub.i] - d (i = 1, 2,..., N); [b.sub.i] = [u.sub.i] - [u.sub.i + 1] (i= 1, 2,..., N - 1); [B.sub.i] [is greater than or equal to] [b.sub.i] [is greater than or equal to] 0 (i = 1, 2,..., N - 1);

where [B.sub.i], is the buffer size, which is to be determined; [[Alpha].sub.i] is the state of machine i, which is a binary random variable; and [U.sub.i] is the maximum service rate of machine i. In this case, [U.sub.i] = 1/[[Tau].sub.i].

As we observed in Section 2.1, the coefficient boundaries defined by [Delta]J/[Delta][x.sub.i] = 0 divide the x-space into mutually exclusive regions in which the production rate u(t) is constant. The scheduler drives the system towards the hedging point, which is the intersection of the coefficient boundaries in x-space. For different positions of the hedging point and boundary shape, the system behaves differently. In the following sections we design a heuristic partition of x-space.

3.1. System behavior specifications

In a way similar to the classical control system design procedure, we specify the system behavior requirements in this section. Then the x-space is partitioned into regions according to the specifications. The following are specifications on system behavior:

(1) when machine 1 fails, keep machine 2 producing without changing the production plan until buffer 1 is empty;

(2) when machine i (i = 2, ..., N - 1) fails, keep machine i - 1 producing without changing the production plan until buffer i - 1 is full, and keep machine i + 1 producing without changing the production plan until buffer N - i is empty;

(3) When machine N fails, keep machine N - 1 producing without changing the production plan until buffer N - 1 is full.

The behavior requirements impose maximum independence among machines for a given WIP distribution. This consideration is essential for system decomposition as well as for the implementation of the production control policy. Consequently we are to a large extent able to make decisions locally.

3.2. The boundary shape in x-space

To ensure that the system behaves as specified, we need to determine the boundaries in x-space. Denote by z = ([z.sub.1], [Z.sub.2], ..., [Z.sub.N]) the hedging point. Consider the situation when the system has reached the hedging point, [x.sub.i] = [z.sub.i] (i = 1, 2,..., N). The production rate decision is [u.sub.i] = d (i = 1, 2, ..., N), so the system stays at the hedging point indefinitely. Then, suppose that at time t after the system has reached the hedging point, all machines but machine k are down ([[Alpha].sub.i] = 1, [[Alpha].sub.i] = 0 for i [is not equal to] k). According to the behavior requirements in the preceding subsection and the capacity constraints in (21), the production decision should be

(23) [u.sub.k] = d end [u.sub.i] = 0 (i [is not equal to] k),

until machine k is either starved or blocked or fails, or one of the other machines is repaired. Before the starvation or blockage or failure or repair occurs, [x.sub.k] is equal to a constant [z.sub.k] and the other production surpluses, [x.sub.i] (i [is not equal to] k), decrease at rate d, according to the system dynamics. This suggests that the coefficient boundary defined by [Delta]J/ [Delta][x.sub.k] = 0 must be perpendicular to the [x.sub.k] axis in x-space. Otherwise, on the boundary, we have [x.sub.k] [is not equal to] constant, which contradicts (23). From the discussion above, to satisfy the system behavior requirements of Section 3.1, the coefficient boundaries in x-space must have the following properties:

(a) the coefficient boundary must be perpendicular to the [x.sub.i] axis in x-space;

(b) all coefficient boundaries intersect each other at the hedging point.

Fig. 6 illustrates the boundary shape in the ([x.sub.i - 1], [x.sub.i]) subspace. Note that the behavior requirements restrict the hedging point to be independent of machine state [Alpha].

[Figure 6 ILLUSTRATION OMITTED]

3.3. The conditional constraints

To prevent chattering, we need a set of conditional constraints to guide the system moving along the boundaries in x-space. The conditional constraints are:

(24) if [x.sub.i] = [z.sub.i] and [B.sub.i] [is greater than] [b.sub.i] [is greater than] 0, then u[[Alpha].sub.i] = d[[Alpha].sub.i] (i = 1, 2,..., N);}

if [b.sub.i] = 0, [u.sub.i] [is greater than or equal to] [u.sub.i] + 1 (i = 1, 2,..., N - 1);}

if [b.sub.i] = [B.sub.i], [u.sub.i] [is less than or equal to] [u.sub.i] + 1 (i = 1, 2,..., N - 1);}

which imply that when the production surplus [x.sub.i] reaches its component of the hedging point, [z.sub.i], at machine i, the flow rate should equal the demand. This is so that chattering on the attractive boundary can never occur (Gershwin et al., 1985). When buffer i is empty, machine i + 1 cannot be faster than the upstream machine. When buffer i is full, machine i cannot be faster than the downstream machine.

3.4. The linear program for real-time control

To ensure that the coefficient boundaries in x-space are perpendicular to axes, and go through the hedging point, we construct the production controller,

(25) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to (21) and (24). The system dynamics are governed by (22). It is possible to show that the linear program above generates the desirable system behavior.

In the linear program, the hedging point ([z.sub.1], ..., [z.sub.N]) and, the buffer sizes, [B.sub.i] (i = 1, ..., N - 1), are still unknown. We show how they may be found in the next subsection.

3.5. Control parameter estimation

In the preceding subsection we formulated a real-time production scheduler as a linear program. In this subsection, we determine the unknown parameters of the controller, namely the hedging point and buffer sizes, by extending the formulations in Section 2.6 to the N-machine, one-part-type case.

3.5.1. Starvation and blockage

Define [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] to be the hedging level of buffer i (see Fig. 6). It is the number of parts in buffer i when the system reaches the hedging point, which satisfies

(26) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Define [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] to be the hedging space of buffer i. It is the room left for more parts in buffer i when the system reaches the hedging point, which satisfies

(27) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be the fractions of time that machine i (i = 1, ..., N) is starved or blocked. The parameters defined above can be determined by solving the following non-linear program (see Appendix A):

(28) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII];

where

[D.sub.i] = [[r.sub.i]/([r.sub.i] + [p.sub.i])] [U.sub.i] (i = 1, 2,..., N).

3.5.2. The buffer sizes and average buffer levels

By definition, the buffer sizes are given by

(29) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

As we discussed in Section 2.6.3, the average buffer levels are

(30) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where [[Delta].sub.i] is the surplus loss at machine i. The average buffer should be truncated to

(31) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

if the following estimate of surplus loss is used (see Appendix A):

(32) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

3.5.3. The hedging point

The hedging point ([z.sub.1], ..., [z.sub.N]) is given (see Appendix A) by

(33) [z.sub.i] = [[bar]x.sub.i] + [[Delta].sub.i] (i = 1, 2,..., N).

To reduce the final product inventory, we would like to minimize the absolute value of surplus [x.sub.N]. The criterion we use to select the hedging point ([z.sub.1], [z.sub.2], ..., [z.sub.N]) is such that

(34) [[bar]x.sub.N] = 0.

From (26), (33), and (34), the hedging point must satisfy

(35) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

3.5.4. The average cycle time and WIP level

The formulation developed in Section 2.6.4 can be easily extended to N-machine production lines. The cycle time includes the waiting times in buffers and processing times on machines. Each part travels through the system with n average rate d. The average waiting time in buffer i is [[bar]b.sub.i]/d. The total waiting time on the average is given by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The total processing time is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Denote by [bar]CT the average cycle time that a part spends in the system:

(36) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Denote by W the average WIP level in the system. By Little's law, we have

(37) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3.6. The algorithm

We have extended the real-time feedback control algorithm to N-machine, single-part-type production lines. The steps of the algorithm are summarized in the following.

Step 1. Collect the input data set, which consists of the failure rates [p.sub.i], the repair rates [r.sub.i], and the processing time [[Tau].sub.i] for machine i(i = 1, 2, ..., N), and the demand d.

Step 2. Calculate the buffer hedging levels, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (i = 1, 2, ..., N), and hedging spaces, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (i = 1, 2, ..., N), and the starvation and blockage fractions for each machine by solving the nonlinear program (28). Then calculate the buffer size for each machine by summing the buffer hedging level and hedging space.

Step 3. Calculate the components of the hedging point, ([z.sub.1], [z.sub.2], [z.sub.N]) according to (35).

Step 4. Using the feedback information of surplus [x.sub.i] and machine state [[Alpha].sub.i] (i = 1, 2, ..., N) calculate the production rates, [u.sub.i] (i = 1, 2, ..., N), in real time by solving the linear program (25).

Step 5. The loading times for each machine are determined by the heuristic staircase strategy. That is, whenever the actual cumulative production is less than the integral of the production rate, load a part into the machine.

Step 6. If the demand or any of the machine parameters changes, go to step 2. Otherwise if the machine state changes or the production surplus reaches a boundary, go to step 4.

The hierarchical structure for production control described in Section 2.7 is applicable here.

3.7. Simulation examples

For the simulation, we still use the three-level hierarchical policy described in Section 2.7. At the top level, the buffer sizes and hedging point are calculated by using a commercially available software package GAMS (Brook et al., 1988). HIERCSIM (Darakananda, 1989) is used for the next two level simulations. The system consists of five machines and four buffers. Each part follows the route: machine 1, buffer 1, machine 2, ..., buffer 4, machine 5. The parameters are chosen as follows:

case 4. [r.sub.1] = 0.5,   [p.sub.1] = 0.3,   [[Tau].sub.1] = 0.5;
        [r.sub.2] = 0.2,   [p.sub.2] = 0.05,  [[Tau].sub.2] = 0.3;
        [r.sub.3] = 0.3,   [p.sub.3] = 0.2,   [[Tau].sub.3] = 0.6;
        [r.sub.4] = 1.2,   [p.sub.4] = 0.1,   [[Tau].sub.4] = 0.4;
        [r.sub.5] = 0.3,   [p.sub.5] = 0.1,   [[Tau].sub.5 = 0.7;

where the unit of r and p is 1/day. The unit of [Tau] is a day. The unit of parts is a lot.

Given that the demand is d = 0.7 lots/day, the buffer sizes and hedging point are calculated at the top level by solving (28) and (35), and listed in Table 4.

[TABULAR DATA 4 NOT REPRODUCIBLE IN ASCII]

Fig. 7 illustrates the results of a single-run simulation of the cumulative production. The straight line is the cumulative demand. The upper curve is the input of the raw parts at machine 1. The lower curve is the output of the final products at machine 5. The dashed lines are the results at the middle level by solving (25).

[Figure 7 ILLUSTRATION OMITTED]

Fig. 8 shows the history of the level of buffer 3, which lies between machine 3 and machine 4. The dashed lines are the second-level result, which is determined by [b.sub.3](t) = [x.sub.3](t) - [x.sub.4](t). The solid lines are the actual count of the parts in the buffer.

[Figure 8 ILLUSTRATION OMITTED]

To see the effects of buffer levels and sizes, we increase the demand to 0.85 lots/day without changing the buffer sizes and the hedging point (Table 4). The simulation result in Fig. 9 shows that the production fell behind the demand. That is, with the buffer levels and sizes in Table 4, the system is starved or blocked too much to achieve the demand, 0.85.

[Figure 9 ILLUSTRATION OMITTED]

Given that the demand is 0.85 lots/day, the desirable buffer sizes and hedging point are calculated and listed in Table 6. With the appropriate buffer sizes and hedging point (Table 6), the actual production follows the demand closely (see Fig. 10).

[Figure 10 ILLUSTRATION OMITTED]

[TABULAR DATA 6 NOT REPRODUCIBLE IN ASCII]

In order to evaluate further the performance of the algorithm, we perform simulations for the following three cases of a five-machine production line. In case 5, all machines have the same failure and repair rate and the same processing time. The system is unbalanced in cases 6 and 7. Each machine is faster than its downstream neighbor in case 6, and the configuration is the other way around in case 7. The system parameters are chosen as in Table 5.

Table 5. The system parameters for cases 5-7

         [r.sub.1]   [r.sub.2]   [r.sub.3]   [r.sub.4]   [r.sub.5]
Case 5      0.5         0.5         0.5         0.5         0.5
Case 6      0.5         0.5         0.5         0.5         0.5
Case 7      0.5         0.5         0.5         0.5         0.5

         [p.sub.1]   [p.sub.2]   [p.sub.3]   [p.sub.4]   [p.sub.5]
Case 5      0.1         0.1         0.1         0.1         0.1
Case 6      0.05        0.1         0.2         0.3         0.4
Case 7      0.4         0.3         0.2         0.1         0.05

         [[Tau].sub.1]  [[Tau].sub.2]   [[Tau].sub.3]  [[Tau].sub.4]
Case 5       0.5             0.5             0.5             0.5
Case 6       0.5             0.5             0.5             0.5
Case 7       0.5             0.5             0.5             0.5

         [[Tau].sub.5]
Case 5      0.5
Case 6      0.5
Case 7      0.5

For case 5, the buffer sizes and the hedging point are calculated and listed in Table 7 for different demand values. The buffer sizes are rounded off to integers. For a given demand, five simulation runs are made with different random number seeds. The starvation and blockage fractions, buffer levels, and WIP inventory are averaged over the five runs. Those quantities are also calculated by using (28), (31), and (37). The results for case 5 are presented in Tables 8-10. Readers are referred to Bai and Gershwin (1993), for the simulation results for the other cases. The observations made in Section 2.8 can also be made for these cases. That is, simulations agree with the estimates better when demand is relatively high than when demand is low, and better for balanced systems than unbalanced systems.

[TABULAR DATA 8-9 NOT REPRODUCIBLE IN ASCII]

Table 7. The hedging point and buffer sizes for case 5

 d    [z.sub.1]   [z.sub.2]   [z.sub.3]   [z.sub.4]   [z.sub.5]
1.6     13.47       10.61       7.76         3.92       1.38
1.4      9.38        8.03       6.44         1.4        1.27
1.2      5.74        5.74       1.61         1.48       1.48
1.0      2.2         2.2        1.88         1.67       1.67
0.6      0.76        0.76       0.76         0.76       0.76

 d    [B.sub.1]   [B.sub.2]   [B.sub.3]   [B.sub.4]
1.6       6           6           8           5
1.4       3           3           10          1
1.2       1           4           5           1
1.0       1           1           1           1
0.6       1           1           1           1
Table 10. The average buffer levels and WIP inventory for case 5

 d     [[bar]b.sub.1]            [[bar]b.sub.2]
            Est.         Sim.       Est.           Sim.
1.6         2.8          4.1        2.9            3.2
1.4         1.1          2.0        1.8            1.3
1.2         0.0          0.7        3.6            2.1
1.0         0.0          0.7        0.2            0.5
0.6         0.0          0.3        0.0            0.3

 d    [[bar]b.sub.3]          [[bar]b.sub.4]            W
          Est.         Sim.      Est.           Sim.   Est.   Sim.
1.6       3.8          3.8       2.6            1.6    16.1   16.5
1.4       5.0          4.6       0.1            0.4    11.6   10.7
1.2       0.6          2.3       0.5            0.4     7.8    7.5
1.0       0.3          0.4       0.6            0.2     3.6    5.0
0.6       0.1          0.2       0.3            0.2     1.8    3.9

4. Comparison with a transfer line model

As a comparison, we consider an alternative production policy. Gershwin (1987) modeled the system described at the beginning of the last section as a tandem queue. The control policy is to produce at the maximum speed of each machine whenever it is possible. The failure and repair rates, processing times, and buffer sizes are given. Then the average production rate of the system and average buffer levels are evaluated. Dallery et al. (1988) developed an efficient algorithm for the transfer line model above.

As an example, we demonstrate the simulation comparison for the following production line.

case 8: [r.sub.1] = 0.3,   [p.sub.1] = 0.1,   [[Tau].sub.1] = 1.0;
        [r.sub.2] = 0.3,   [p.sub.2] = 0.1,   [[Tau].sub.2] = 1.0;
        [r.sub.3] = 0.3,   [p.sub.3] = 0.1,   [[Tau].sub.3] = 1.0;
        [r.sub.4] = 0.3,   [p.sub.4] = 0.1,   [[Tau].sub.4] = 1.0;
        [r.sub.5] = 0.3,   [p.sub.5] = 0.1,   [[Tau].sub.5] = 1.0.

For the transfer line model, we determine the optimal values of the average buffer levels and the average WIP with the use of the algorithm of Dallery et al. (1988) as well as the simulation package XCELL + (Conway et al., 1990). We manually change the buffer sizes until the minimal total buffer space is found such that a specified production rate is achieved. The production line is also run under the control policy developed in the previous section. The quantities of interest are estimated using the formulations developed in the present paper (B&G) as well as the simulation package HIERCSIM. The results are listed in Table 11. From the simulation results, we see that to achieve a production rate d = 0.6, the total buffer space is at least 73 lots and the average total amount of material in buffers is 32 lots in the queueing system. Those quantities are reduced to 20 and 9 lots respectively if the real-time algorithm is applied. That is a 72.6% reduction in total buffer space and a 60.3% reduction in average amount of material in buffers. More simulation cases can be found in Bai and Gershwin (1993). Reductions in total buffer space and average buffer levels can be observed when the real-time control policy is adopted for all the cases we have considered. The reductions range from 20.8% to 72.6% in total buffer space and from 38.1% to 60.3% in total buffer level.

Table 11. The buffer sizes and average buffer levels for
case 8 (d = 0.6)

Case 8, d = 0.6   [B.sub.1]     [B.sub.2]     [B.sub.3]   [B.sub.4]
DDX (Est.)           6              8             7            5
XCELL+ (Sim.)       12             18            21           22
B&G (Est.)           2              6             6            2
HIERCSIM (Sim.)      3              8             7            2

Case 8, d = 0.6   [[bar]b.sub.1]  [[bar]b.sub.2]    [[bar]b.sub.3]
DDX (Est.)              3.8           4.4                3.5
XCELL+ (Sim.)           8.9           11                 8.4
B&G (Est.)              0.5           3.0                2.2
HIERCSIM (Sim.)         1.1           3.3                3.9

Case 8, d = 0.6   [[bar]b.sub.4]  [Sigma][B.sub.i]
DDX (Est.)             2.0               26
XCELL+ (Sim.)          3.7               73
B&G (Est.)             0.5               14
HIERCSIM (Sim.)        0.7               20

Case 8, d = 0.6    [Sigma][[bar]b.sub.i]
DDX (Est.)                 13.7
XCELL+ (Sim.)              32.0
B&G (Est.)                  6.2
HIERCSIM (Sim.)             9.0

5. Summary

A real-time algorithm has been developed for scheduling single-part-type production lines with failure prone machines. The objectives are to keep the actual production close to demand and to keep the WIP inventory level low. In the algorithm we combine strategic decision making, tactical planning, and operational scheduling into a hierarchy to respond to events with different frequencies. For long-term capacity planning, the buffer size and the average buffer level are determined by a control parameter estimator such that the system has enough capacity to achieve a given demand. Production rates are adjusted in real time by the scheduler whenever a machine fails or is starved or blocked. The loading times for each machine are determined by using a staircase dispatching rule in accordance with the production rates. The production flow model and the bottom-up control system design approach provide good tractability and predictability of system performance as well as system behavior. The simulation results verify that the algorithm works well.

The method developed in this paper is extended to multiple-part-type and reentrant systems in Bai (1991b) and Bai and Gershwin (1990, 1991). The algorithm is also modified in Bai (1991a) such that the formulation for buffer hedging levels and spaces becomes linear. In many cases, the linear formulation gives us equally good results without the difficulty of solving the nonlinear programming problem. Further extensions to general manufacturing networks and application to semiconductor fabrication processes are future research topics.

Acknowledgments

The authors thank an anonymous reviewer whose comments were very helpful in improving the paper. We thank B. Darakananda and J.H. Burhanpurwala for their help on simulations. We also thank Professor D.J. Elzinga and Professor R.L. Francis for reading the manuscript, and for their thoughtful comments. This research was supported by the Defense Advanced Research Project Agency under contracts N00014-85-K-0213 and MDA-972-88-K-008.

References

Afentakis, P., Gavish, B. and Karmarkar, U. (1984) Computationally efficient optimal solutions to the lot-sizing problem in multi-stage assembly systems. Management Science, 30(2), 222-239.

Akella, R. and Kumar, P.R. (1986) Optimal control of production rate in a failure prone manufacturing system. IEEE Transactions on Automatic Control, AC-31(2), 116-126.

Akella, R., Choong, Y.F. and Gershwin, S.B. (1984) Performance of hierarchical production scheduling policy. IEEE Transactions on Components, Hybrids, and Manufacturing Technology, CHMT-7(3)

Bai, S.X. (1991a) Linear formulation of control parameter estimation for real-time production scheduling. Master's thesis, MIT.

Bai, S.X. (1991b) Scheduling manufacturing systems with work-in-process inventory control. Ph.D. thesis, MIT.

Bai, S.X. and Gershwin, S.B. (1990) Scheduling manufacturing systems with work-in-process inventory control: multiple-part-type systems. MIT technical report OR-230-90, Operations Research Center.

Bai, S.X. and Gershwin, S.B. (1991) Scheduling manufacturing systems with work-in-process inventory control: reentrant systems. MIT technical report OR-240-91, Operations Research Center.

Bai, S.X. and Gershwin, S.B. (1993) Scheduling manufacturing systems with work-in-process inventory control: single-part-type systems. University of Florida technical report 93-16, Department of Industrial and System Engineering.

Bitran, G.R., Haas, E.A. and Hax, A.C. (1981) Hierarchical production planning: a single-stage system. Operations Research, 29(4), 717-743.

Brook, A., Kendrick, D. and Meeraus, A. (1988) GAMS: A User's Guide. The Scientific Press.

Burman, D.Y., Gurrola-Gal, F.J., Nozari, A., Sathaye, S. and Sitarik, J.P. (1986) Performance analysis techniques for IC manufacturing lines. AT&T Technical Journal, 65(4).

Chen, H., Harrison, M., Mandelbaum, A., Van Ackere, A. and Wein, L. (1988) Empirical evaluation of a queueing network model for semiconductor wafer fabrication. Operations Research, 36(2).

Conway, R., Maxwell, W., McClain, J.O. and Thomas, L.J. (1988) The role of work-in-process inventory in serial production lines. Operations Research, 36(2).

Conway, R., Maxwell, W.L., McClain, J.O. and Worona, S.L. (1990) User's Guide to XCELL+: Factory Modeling System (Release 4.0). The Scientific Press.

Dallery, Y., David, R. and Xie, X.-L. (1988) An efficient algorithm for analysis of transfer lines with unreliable machines and finite buffers. IIE Transactions, 20(3).

Darakananda, B. (1989) Simulation of manufacturing process under a hierarchical control structure. Master's thesis, MIT.

Dempster, M.A.H., Jansen, L., Fisher, M.L., Lageweg, B.J., Lenstra, J.K. and Rinnooy Kan, A.H.G. (1981) Analytical evaluation of hierarchical planning systems. Operations Research, 29(4), 707-716.

Eleftheriu, M.N. (1989) On the analysis of hedging point policies of multi-stage production manufacturing systems. Ph.D. thesis, Rensselaer Polytechnic Institute.

Fine, C.H. and Graves, S.C. (1990) A tactical planning model for manufacturing subcomponents of mainframe computers. Journal of Manufacturing and Operations Research, 2, 4-34.

Gershwin, S.B. (1987) An efficient decomposition method for the approximate evaluation of tandem queues with finite storage space and blocking. Operations Research, 35(2), 291-305.

Gershwin, S.B. (1989) Hierarchical flow control: a framework for scheduling and planning discrete events in manufacturing systems. Proceedings of the IEEE, Special Issue on Dynamics of Discrete Event Systems, 77(1), 195-209.

Gershwin, S.B., Akella, R. and Choong, Y.F. (1985) Short-term production scheduling of an automated manufacturing facility. IBM Journal of Research and Development, 29(4), 392-400.

Glassey, C.R. and Resende, M.G.C. (1988) Close-loop job release for VLSI circuit manufacturing. IEEE Transactions on Semiconductor Manufacturing, 1(1).

Graves, S.C. (1981) A review of production scheduling. Operations Research, 29(4), 646-675.

Graves, S.C. (1982) Using lagrangean relaxation techniques to solve hierarchical production planning problems. Management Science, 28(3), 260-275.

Graves, S.C. (1986) A tactical planning model for a job shop. Operations Research, 34(4), 522-533.

Harrison, J.M. and Wein, L.M. (1990) Scheduling networks of queues: heavy traffic analysis of a two-station closed network. Operations Research, 38(6).

Hax, A.C. and Meal, H.C. (1975) Hierarchical integration of production planning and scheduling. North Holland/TIMS, Studies in Management Sciences, vol. 1 (Logistics).

Kimemia, J. and Gershwin, S.B. (1983) An algorithm for the computer control of a flexible manufacturing system. IIE Transactions, 15(4), 353-362.

Lageweg, B.J., Lenstra, J.K. and Rinnooy Kan, A.H.G. (1977) Job-shop scheduling by implicit enumeration. Management Science, 24(4), 441-450.

Lageweg, B.J., Lenstra, J.K. and Rinnooy Kan, A.H.G. (1978) A general bounding scheme for the permutation flow-shop problem. Operations Research, 26(1), 53-67.

Maimon, O.Z. and Choong, Y.F. (1987) Dynamic routing in reentrant flexible manufacturing system. Robotic and Computer-Aided Manufacturing, 3, 295-300.

Maimon, O.Z. and Gershwin, S.B. (1988) Dynamic scheduling and routing for flexible manufacturing systems that have unreliable machines. Operations Research, 36(2), 279-292.

Newson, E.F.P. (1975a) Multi-item lot size scheduling by heuristic, part I: with fixed resources. Management Science, 21(10), 1186-1193.

Newson, E.F.P. (1975b) Multi-item lot size scheduling by heuristic, part II: with variable resources. Management Science, 21(10), 1194-1203.

Papadimitriou, C.H. and Kannelakis, P.C. (1980) Flowshop scheduling with limited temporary storage. Journal of the ACM, 27(3).

Rishel, R. (1975) Dynamic programming and minimum principles for systems with jump Markov disturburces. SIAM Journal on Control, 13(2).

Sharifnia, A. (1988) Production control of a manufacturing system with multiple machine states. IEEE Transactions on Automatic Control, AC-33(7), 620-625.

Tsitsiklis, J. (1982) Convexity and characterization of optimal policies in a dynamic routing problem. MIT technical report LIDS-R-1178.

Van Ryzin, G. (1987) Control of manufacturing systems with delay. Master's thesis, MIT.

Wagner, H.M. and Whitin, T.M. (1958) Dynamic version of the economic lot size model. Management Science, 5(1), 89-96.

Wein, L.M. (1988) Scheduling semiconductor wafer fabrication. IEEE Transactions on Semiconductor Manufacturing, 1(3).

Wein, L.M. (1990) Scheduling networks of queues: heavy traffic analysis of a two-station closed network with controllable inputs. Operations Research, 38(6).

Zeghmi, A.H. (1985) Inventory buffers for a production line with controllable production rates. Ph.D. thesis, MIT.

Appendix A

Here we briefly outline the derivation of the formulations for the control parameters discussed in Section 3. Readers are referred to Bai (1991b) for more details.

Starvation and blockage

The starvation fraction of machine 1. Because we assume that machine 1 is never starved, we have

(38) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The starvation fraction of machine i (i = 2, ..., N). Consider the average up-down cycle of machine i - 1, which has a length of [L.sub.i - 1] = 1/[r.sub.i-1] + 1/[p.sub.i-1]. The average portion of [L.sub.i - 1] in which machine i - 1 is down or starved is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. We divide this into three parts, [[Beta].sub.1], [[Beta].sub.2], and [[Beta].sub.3]. Let [[Beta].sub.1] be the portion of [L.sub.i-1] during which machine i is down or blocked when, machine i - 1 is down or starved. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first factor is the length of the average up down cycle of machine i - 1. The second is the fraction of the cycle in which machine i - 1 is down or starved. The third factor is the fraction of the cycle in which machine i is down or blocked.

Let [[bar]u.sub.1] be the average production rate of machine i when it is operating. It is determined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Let [[Beta].sub.2] be the average amount of time during [L.sub.i - 1] in which machine i operates while machine i - 1 is down or starved. We assume that the level of buffer i - 1 is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] when machine i - 1 fails. Machine i operates until it removes this much material from the buffer, and [[Beta].sub.2] is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Let [[Beta].sub3] be the portion of [L.sub.i -1] in which machine i is starved while machine i - 1 is down or starved. Similarly to [[Beta].sub.1] it is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first factor is again the length of the average up-down cycle of machine i- 1. The second is 1 because machine i can be starved only if machine i - 1 is down or starved. The third factor is the fraction of the cycle in which machine i is starved. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The [Beta] satisfy

(39) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Substituting the [Beta] into (39), and manipulating, leads to

(40) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The blockage fraction of machine i (i = 1, 2, . . ., N-1). Assume that the average spare space in buffer i is z,. at the instant that machine i + 1 goes down. By similar reasoning as for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the blockage fraction of machine i is governed by

(41) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The blockage fraction of machine N. Because we assumed that machine N is never blocked, the blockage fraction of machine N is

(42) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

To ensure that the system has enough capacity to achieve the demand, the starvation and blockage fractions must satisfy

(43) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

By putting (38), (40), (41), (42), and (43) together, we form an optimization problem to minimize the sum of the buffer spaces:

(44) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

subject to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that the nonlinear program above is separable when N = 2. In that case, it has a closed-form solution as presented in (13).

The hedging point and surplus loss

Because machines are unreliable and can be starved or blocked, there is a difference between the hedging point ([z.sub.1],...,[z.sub.N]) and the average surplus ([[bar]x.sub.1], . . .,[[bar]x.sub.N]), which is the time average over the planning horizon of ([x.sub.1], . . .,[x.sub.N]). The relation can be written as

(45) [z.sub.i] + [[bar]x.sub.1] + [[Delta].sub.i] (1 = 1,...,N),

where [[Delta].sub.i] (i = 1, ...,N) is the average surplus loss at machine i, which is the average amount that [x.sub.i], deviates from [z.sub.i].

In the following, we introduce a heuristic to estimate the average surplus loss. The average surplus loss [[Delta].sub.i] (i = 1, . . ., N) consists of three components caused by failure, starvation, and blockage. In the heuristic, the three components are calculated separately. Fig. A1 illustrates the typical surplus loss due to failures. The shaded area is the total surplus loss due to failure during an average machine i up down interval. The area of the shaded region is equal to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [t.sub.ci] is the average catch-up time needed for machine i to move back to the hedging point after repair. It is given by

[t.sub.ci] = d/[r.sub.i] ([U.sub.i] - d),

where d/[r.sub.i] is the average surplus loss during a machine i breakdown and ([U.sub.i] - d) is the average speed for machine i to recover after repair.

[Figure A1 ILLUSTRATION OMITTED]

Let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be the average surplus loss due to failures for machine i. Dividing the shaded area of Fig. A1 by the average time that machine i is up once and down once, (1/[r.sub.i] + 1/[p.sub.i]), leads to

(46) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Let and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] be the average surplus losses due to starvation and blockage. With similar reasoning as for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], they are approximated by

(47) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(48) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The average surplus loss then given by

(49) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In addition to minimizing WIP inventory, we wish to optimize the system performance by delivering on time. We would like to minimize both final product inventory and backlog. A plausible choice of the hedging point ([z.sub.1], . ., [z.sub.N]) is such that

(50) [[bar]x.sub.N] = 0.

From (11), (45), and (50), the hedging point should satisfy

[z.sub.N] = [[Delta].sub.N],

(51) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that we have calculated these parameters without reference to the penalty function g of (20). Consequently, this is a crude approximation of the optimal parameters. However, the important performance measures are in process inventory and cycle times, and the simulation results indicate that we get good performance with parameters computed this way. Further research should be devoted to improved approximations.

Biographies

Sherman X. Bai received his Ph. D. in Operations Research from the Sloan School of Management of MIT, in September 1991. Since September 1991, he has joined the faculty of Industrial and Systems Engineering Department at the University of Florida, Gainesville. His research interests are in the areas of Manufacturing Systems Modeling and Analysis, Digital Simulation, Optimal Control, and Dynamic Game Theory.

Stanley B. Gershwin received the B. S. degree in Engineering Mathematics from Columbia University in 1966, and the M. A. and Ph.D. degrees in Applied Mathematics from Harvard University in 1967 and 1971 respectively. In 1970-1971, he was employed by the Bell Labs and from 1971 to-1975 he worked for the Drape Laboratory for Information and Decision Sciences of MIT. In 1986-1987, he was a Professor in the Manufacturing Engineering Department of Boston University . He is currently a Senior Research Scientist in the Department of Mechanical Engineering of MIT, and is affiliated with the Laboratory for manufacturing and Productivity. His major research interest is on systems issues in manufacturing.

Received January 1992 and accepted June 1993