IIE Transactions (2002) 34, 891–903

Optimal design of stochastic production lines: a dynamic programming approach KAREN L. DONOHUE1 , WALLACE J. HOPP2 and MARK L. SPEARMAN3 1

Department of Operations and Management Science, The Carlson School of the University of Minnesota, Minneapolis, MN 55455, USA E-mail: [email protected] 2 Department of Industrial Engineering and Management Sciences, Northwestern University, Evanston, IL 60208-3119, USA E-mail: [email protected] 3 School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0205, USA E-mail: [email protected] Received September 2000 and accepted December 2001

We consider the problem of choosing the number and type of machines for each station in a new production line where the sequence of processes (i.e., manufacturing recipe) has already been established. We formulate a model to minimize cost (investment plus operating) subject to constraints on throughput and cycle time. Using queueing network approximations within a dynamic programming framework, we develop a line design algorithm that works in station-wise fashion. For computational tractability, we must discretize a continuous state space. However, we are able to compute bounds on the error in the cost function as a guide to the appropriate choice of grid size. We conclude by applying our algorithm to an industrial problem that motivated this work.

1. Introduction Investment in process technology is a key determinant of a plant’s performance and a potentially significant source of competitive leverage for the firm. Therefore, we were surprised to find in discussions with engineers and managers in several manufacturing facilities (e.g., printed circuit board, flat panel screens) that the procedures used to select process equipment are often very crude. Most approaches are variants of the following: determine a target throughput rate for the production line and select the cheapest combination of equipment at each stage whose rated capacity (including detractors for downtime, setups, etc.) is greater than or equal to this target. Some designers do not even account for possible detractors. Not surprisingly, these same engineers and managers complained that the throughput of their installed lines fell well-short of design targets with average cycle times (flow times) considerably longer than management would like. The reason, of course, was that in most cases no consideration was given to variability and the congestion it causes between process stages. The plants we encountered who did consider variability did so via simulation and, for the most part, only performed simulation studies after almost all process equipment selections had been made. The long run times required for statistical validity makes 0740-817X


2002 ‘‘IIE’’

simulation better suited to fine tuning a line (e.g., establishing buffer sizes, staffing levels, incremental capacity changes) than searching through the large number of equipment options available in the initial line design phase. In response to these shortcomings, we introduce a procedure for designing a new production line where the sequence of operations (i.e., manufacturing recipe) has already been determined and the basic problem is to choose the type and number of machines to install at each station. In keeping with the concerns of our industrial clients, we frame this problem as one of minimizing cost (purchase plus operating) subject to constraints on throughput and cycle time. Because, in most plants, cycle time consists largely of queue time, it is essential for our approach to explicitly consider queueing. Since choosing equipment options from a very large set requires us to evaluate many configurations, we do this using analytic queueing network approximations instead of simulation. Finally, to facilitate the process of searching through configurations, we develop a Dynamic Programming (DP) procedure. The result is an efficient procedure for designing tandem production lines that provides better consideration of variability effects than existing methods that consider only capacity and is much faster than simulation-based

892 methods. This procedure gives the designer of a line a tool with which to sift through the many equipment options quickly and identify a small number of configurations that can be further refined using more detailed simulation analysis. It also allows explicit examination of the tradeoff between variability and capacity (i.e., a less variable line requires less capacity to meet given throughput and cycle time constraints) which was lacking in previous design models. The remainder of the paper is organized as follows. Section 2 gives an overview of previous work, while Section 3 defines the problem and associated notation. Section 4 describes the DP solution procedure, discretization approach, and associated error bound. Section 5 presents a numerical example based on a printed circuit board facility design problem, and Section 6 concludes the paper.

2. Previous work The problem of optimizing the design of manufacturing systems is not new. Over the past 10 years various researchers have worked to incorporate advances in queueing theory into a larger optimization framework for selecting optimal process rates, work loads, or number of servers. Walrand (1993) offers an excellent survey of queueing network models; Kouvelis and Tirupati (1991) give a thorough review of capacity decision models for open and closed manufacturing systems. How queueing models can be used to guide design decisions depends on the decision variables involved and whether the system is open or closed. Motivated by the printed circuit board facilities we studied, we are interested primarily in open manufacturing lines where production releases are controlled by an aggregate production rate. Our goal is to solve a ‘‘targeting problem’’ which minimizes the cost of resource investment while meeting targets for average cycle time and other possible performance measures. In the short-term, manufacturing managers can try to achieve a given cycle time target by adjusting the processing rate of the resources at their disposal. Viewed in this way, capacity choices are modeled as continuous but subject to feasibility bounds which vary by resource. Bitran and Tirupati (1989a) examine this problem, choosing process rates to minimize capacity cost subject to a desired performance target, in the context of an open network of single server queues. They simplify the problem by assuming the coefficients of variation of the process and interarrival times at each station are fixed (i.e., independent of capacity choices). The problem is then stated as a convex program which can be solved using conventional convex programming methods or a greedy heuristic. This simplifying approximation is most reasonable for large networks where the addition of capacity at one station is unlikely to have a large affect on other

Donohue et al. stations. Bitran and Sarkar (1994) propose an exact iterative algorithm which relaxes this assumption. However, this causes the objective function to be non-convex, so their search is only guaranteed to yield a locally optimal solution. Other work that considers process rate selection includes Schweitzer and Seidmann (1991), who offer an economic analysis of capacity utilization rates in a Flexible Manufacturing System (FMS) producing a single part type, and Lawrence and Buss (1992), who develop indices for the optimal capacity and associated utilization rates for stations in a queueing network. Both papers concentrate on the economic aspects of capacity decisions and assume exponentially distributed service and interarrival times within the manufacturing system. Another line of research in manufacturing system design focuses on the longer term decision of selecting the number of resources to purchase of a pre-specified type to, again, minimize cost subject to performance targets. Boxma et al. (1990) develop a greedy heuristic for guiding resource quantity selection in a manufacturing system consisting of a network of M=M=m queues. Van Vliet and Rinnooy Kan (1990) extend the results of Boxman et al. (1990) to systems with general service rates. Like Bitran and Tirupati (1989a), they assume the coefficients of variation of the interarrival times at each station are not affected by the number of resources. Connors et al. (1993) consider the problem of designing a new wafer fabrication line using queueing network approximations. They present a greedy heuristic to determine the minimum number of tools needed to achieve a mean cycle time target for a general system (which includes multiple routings, rework, recycle). Although they do not give numerical results to support the effectiveness of this approach nor establish whether their objective function is convex, they do report that their algorithm has been used successfully at IBM to arrive ‘‘at a tool set that had a significantly lower cost than the set that had been proposed.’’ Finally, Frenk et al. (1994) offer an improved algorithm for the exponential system considered by Boxma et al. (1990). Their algorithm takes into account a major limitation of the greedy approach in solving discrete choice problems. Namely, it recognizes when the solution is in close proximity of feasibility, and therefore when a small investment with little ‘‘bang’’ may be most desirable. Our research focuses on the initial greenfield stage of manufacturing system design where designers start from a clean slate in terms of their resource choices (e.g., there is no carry over of machinery from a previously configured production line). At this stage, choices must be made for both the types and quantities of resources to purchase. The targeting problem is more difficult in this context since resource choices are discrete and the associated cost function is not convex. Bitran and Tirupati (1989b) consider this problem when capacity choices vary by service time and number of resources. Their solution

893

Optimal design of stochastic production lines approach assumes once again that the coefficient of variation of process and interarrival times remain constant with changes in capacity. Even with this assumption, the underlying integer program is difficult to solve. They propose an efficient heuristic to solve the approximation based on a linear relaxation of the IP. We are interested in a more general version of the greenfield design problem where the trade-offs between process capacity and variability are explicitly measured. This restricts us from using the fixed output coefficient of variation assumption used by Bitran and Tirupati (1989b) and others. Our discussions with manufacturing managers lead us to believe that resources must be compared on a series of dimensions including reliability, natural process variance, and changeover time, as well as average process time and cost. Our work appears to be the first to consider this wide range of discrete resource characteristics in the design decision. These added complexities make the mixed integer program difficult to solve using linear relaxation or greedy heuristics. Donohue (1994) shows, for example, that a greedy heuristic performs poorly as a tool to guide greenfield design decisions in this general setting since the productive benefits of earlier decisions are likely to be negated by later ones. The dynamic programming procedure described in this paper solves the general greenfield design problem for open manufacturing lines by discretizing a continuous state space defined by our cycle time target. This approach allows us to capture information for sensitivity analysis as part of the algorithm’s output (much like the greedy approach proposed by Bitran and Tirupati (1989a) and others for previous problems), as well as control the numerical complexity of the algorithm itself.

3. Problem statement and DP formulation The basic facility design problem is described as follows. Let i denote a proposed workstation (i.e., activity or operation) where i ¼ 1; . . . ; N . We are given a set of possible resources Qi (e.g., machine or worker options) to select from for each workstation i, where jQi j denotes the number of resource types in the set. We assume each resource available at station i is completely characterized by the first and second moment of its effective process time. In other words, its effective mean ðsei Þ and squared coefficient of variation ðc2ei Þ. Since we are speaking of effective processing times, sei and c2ei must reflect a wide range of capacity detractors and variability sources including machine breakdowns and repairs, setup times, operator-induced outages, in addition to non-uniformity in the process itself. These additional sources of capacity loss and variability inflation can be divided into two classes: preemptive outages, which interrupt job processing, and non-preemptive outages, which occur between jobs. For each type of pre-

emptive outage (e.g., equipment failures), we need to estimate the Mean Time to Failure (MTTF) and Mean Time To Repair (MTTR). For each type of non-preemptive outage (e.g., setups) we need to estimate the mean and standard deviation of outage time, as well as the average number of jobs between outages. With these, the sei and c2ei parameters can be computed using the methods given in Hopp and Spearman (2000). Note, however, that since we may be speaking of new equipment which has not been put into use yet, we may need to estimate some of these parameters on the basis of historical performance (e.g., by using values seen in the past for similar types of equipment). We assume the cost of a resource, ki , reflects both its annual operating and amortized investment costs. Let xi denote a particular resource description (ki ; sei ; c2ei ) that could be chosen to support operations at station i, i.e., xi 2 Qi . For simplicity of presentation, we assume equipment cost is a linear function of the quantity purchased, so that Cðxi ; mi Þ ¼ ki mi ; represents the cost of purchasing and operating mi resources of type xi . However, any reasonable (i.e., nondecreasing in mi ) cost function could be substituted without affecting the validity of our model framework. For practical purposes, we assume each resource set Qi , includes only non-dominated resources. We define a resource option xai as dominated if there exists xbi 2 Qi such that sbei  saei , cbei  caei , and kib  kia . That is, if option xbi is cheaper, faster, and less variable than option xa1 . Because of the linear structure of our objective function, we know that any dominated resource is sub-optimal. (Of course, if we use another objective function, then a different definition of dominance is required, but the general concept still applies.) Let sðxi ; mi Þ denote the average sojourn time at station i given mi resources of type xi . An exact expression for sojourn time in a general multi-server queue ðG=G=mi Þ is not known. However, several good approximations exist (see Buzacott and Shanthikumar (1993) for an overview). One class of approximations that is particularly wellsuited to our problem makes use of a parametric decomposition approach based on the simplifying assumption that the departure process of each station is renewal. Under this assumption, the arrival process to a station is simply the departure process from an upstream station. This allows each station to be analyzed separately, then connected to its proceeding station by the resulting departure process. For our analysis, we make use of the approximations developed by Whitt (1993). For example, we compute the squared coefficient of variation (scv) of the departure process at station i by q2ei 2 c2di ¼ 1 þ ð1  q2ei Þðc2di1  1Þ þ pffiffiffiffiffi ðc  1Þ: mi ei

ð1Þ

894

Donohue et al.

Of course, other approximations based on the parametric decomposition principle could also be used. Note that here the scv of the department process at station i is a modified weighted average of the squared coefficients of variation of the station i’s arrival and production processes (c2di1 and c2ei , respectively). The weighting depends on the station’s utilization rate, with c2di more closely approximating c2ei as utilization increases. As noted previously, our objective is to determine the optimal resource type and quantity for each station to minimize total resource cost while meeting or exceeding an average cycle time target, CTT , and throughput target HT for the system. We also consider a maximum threshold for the scv of the departure process, c2T , as a control on service. The general design problem (IP-G) is an integer programming with nonlinear constraints: N X (IP-G) min Cðxi ; mi Þ; i¼1

Subject to

N X

sðxi ; mi Þ  CTT ;

i¼1

ðhsei Þ=mi  1; 8 i; h ¼ HT ; c2dN  c2T ; xi 2 Qi ; i ¼ 1; . . . ; N ; mi 2 Iþ ; where h denotes the release rate of products into the system and c2dN denotes the squared coefficient of variation of the departure process from the system (i.e., from the last station), and sei is defined by resource choice xi . This problem can be transformed into a forward dynamic program, with the following additional notation. Define fi ðRi ; c2di Þ as the minimum possible cost to achieve an average remaining cycle time Ri and squared departure process coefficient of variation c2di through station i. We can express the coefficient of variation of the departure process at station i  1 as a function of c2di , xi ¼ ðsei ; c2ei ; ki Þ, and mi by simple manipulation of Equation (1),   q2ei 2 1 2 2 2 cdi1 fxi ; mi ; cdi g ¼ cdi  1  pffiffiffiffiffi ðcei  1Þ þ 1: mi ð1  q2ei Þ This expression tells us what the departure scv needs to be at station i  1 to support a particular c2di at station i under resource choice ðxi ; mi Þ. Similarly, Ri1 ¼ Ri þ

sðxi ; mi ; c2di1 fxi ; mi ; c2di gÞ;

defines what the average remaining cycle time needs to be at station i  1 to support a particular Ri under resource choice ðxi ; mi Þ. These expressions can be used together to define the feasible search space for each state ðRi ; c2di Þ at stage i. The general problem is expressed by the following DP recursion,

ðDP-G1Þ f0 ðCTT ; c2a Þ ¼ 0;

n fi ðR; c2di Þ ¼ minðxi 2Qi ;mi 2Iþ Þ Cðxi ; mi Þ  þfi1 ½R þ sðxi ; mi ; c2di1 fxi ; mi ; c2di gÞ; o ½c2di1 fxi ; mi ; c2di g ;

where fN ð0; c2T Þ denotes the final solution. However, this formulation is not tractable since both R and c2di are continuous variables and thus the problem has an infinite state space. In the next section, we develop a method for discretizing the state space to make the problem tractable and develop a means to track the error introduce by this discretization.

4. Discretizing the state space To facilitate our discussion, we focus first on a simple case where the system’s interarrival process and all candidate machines process times are exponentially distributed. 4.1. Exponential systems In the exponential case, the DP formulation reduces to ðDP-E1Þ f0 ðT Þ ¼ 0; fi ðRÞ ¼ minðxi 2Qi ;mi 2Iþ Þ fCðxi ; mi Þ þfi1 ðR þ sðxi ; mi ÞÞg: This problem is much easier, since c2di ¼ 1 8 i and thus the state space is defined only in terms of R. To examine the impact of discretizing R, suppose an upper bound exists on the amount of money we can invest at any station i. This is an assumption most companies make, either explicitly or implicitly, when designing a factory. It then follows that there are a finite number of ‘‘breaks’’ in the sojourn time range ð0; T Þ at station i, where each break corresponds to a change in the minimum cost resource choice. Note that if a resource/ quantity combination can achieve a sojourn time of S, it can also maintain any sojourn time S 0  S simply by waiting longer to move finished products. Let gi ðSÞ represent the solution to the subproblem ðSP-EÞ min Cðxi ; mi Þ; Subject to sðxi ; mi Þ  S; ðhsei Þ=mi  1; h ¼ HT ; xi 2 Qi ; mi 2 Iþ ; where gi ðSÞ signifies the cheapest resource/quantity combination that supports an average sojourn time of S or less at station i. This subproblem exhibits two convenient properties.

895

Optimal design of stochastic production lines

Fig. 1. Error tracking at stage i  1 for exponential case.

Proposition 1. The optimal resource cost at a given station i is non-increasing in S. Proof. The proof of Proposition 1 is straight forward. Note that gi ðSÞ represents the optimal resource cost given S. An increase in S corresponds to a loosening of the constraint in (SP-E) which implies the optimal resource cost will either decrease or remain unchanged. j This leads to the following Corollary.

between Rak and Rk1 represents R values which could be overlooked by the algorithm. These gaps represent potential error sources at stage i  1. Tracking a bound on the error passed from stage to stage in (DP-E2) therefore yields a lower bound on the optimal solution to problem (DP-E1). Letting Li ðRj Þ denote a lower bound on fi ðRj Þ for a given state j and stage i, we have Theorem 1. The error introduced by discretizing R is bounded above by fNd ð0Þ  LN ð0Þ, where LN ð0Þ is given by the following recursion.

Corollary 1. For a given station i, if fi ðR2 Þ  fi ðR1 Þ and R2 > R1 then fi ðR1 Þ is suboptimal. Corollary 1 suggests that if we discretize R, then the solution for a particular state, say fi ðR0 Þ, will also be exact for some limited range of states. In other words, there exists  > 0 such that fi ðR0 þ dÞ ¼ fi ðR0 Þ for 0  d  . This structure is useful for tracking a bound on the error introduced by discretizing R. Suppose we choose a grid Gi for R, where jGi j denotes the number of points in the grid and Gi ¼ ð0; T =jGi j; 2T =jGi j; . . . ; T Þ: We assume the grid size for R is identical at each stage. (Stage dependent grids could also be supported within our framework in cases where resource capabilities vary severely from station to station. However, the notation needed to describe stage dependent grids is sufficiently more complex so for simplicity we focus here on the single grid G.) Let Rj denote a specific remaining cycle time in the grid. Assuming grid points are indexed in increasing order, the discrete-version of the dynamic program is expressed as ðDP-E2Þ

f0d ðT Þ fid ðRj Þ

¼ 0;

n o d ¼ mink  j Cðgi ðRak;i1  Rj ÞÞ þ fi1 ðRk Þ ; j ¼ 1; . . . ; jGj; i ¼ 1; . . . ; N :

where Rak;i1 is the actual remaining cycle time corresponding to grid point k at stage i  1 and fNd ð0Þ is the solution to the discretized problem. Note that Rak;i1  Rk . Typically, we expect Rak;i1 to be strictly greater than Rk because of the discrete nature of the function gi . The solution to the discrete problem (DP-E2) provides an upper bound on the solution to the continuous problem (DP-E1). Corollary 1 shows that no error is introduced in the interval Rk  R  Rak for stage i  1 (see Fig. 1). The gap

L1 ðRj Þ ¼ f1d ðRj Þ; 8j ¼ 1; . . . ; jGj:

Li ðRj Þ ¼ min Cðgi ðRkþ1  Rj ÞÞ þ Li1 ðRk Þ ;

ð2Þ

j ¼ 1; . . . ; jGj; i ¼ 2; . . . ; N :

ð3Þ

kj

Proof. From stage 0 to stage 1 no error is passed so the functional equation is exact for each Rj . Thus, Equation (2) follows immediately. Equation (3) is shown by induction on i. Case i ¼ 2: Consider the calculation of f2d ðRj Þ for any j ¼ 1; . . . ; jGj. This value may differ from f2 ðRj Þ because values of R in ranges ðRa1 ; R2 Þ; . . . ; ðRaj1 ; Rj Þ are not considered. Let 1k ¼ ðRak ; Rkþ1 Þ denote a specific gap in R at stage i. Note that f1 ðRÞ  f1 ðRak Þ ¼ f1 ðRk Þ and Cðgi ðR  Rj ÞÞ  CðgðRkþ1  Rj ÞÞ for R 2 1k and any j. Let h2j be the minimum cost to support Rj at stage 2, given that R 2 1k is achieved at stage 1. Then h2k  Cðg2 ðRkþ1  Rj ÞÞ þ f1d ðRk Þ; and f2 ðRj Þ  min Cðg2 ðRkþ1  Rj ÞÞ þ f1d ðRk Þ; kj

¼ min Cðg2 ðRkþ1  Rj ÞÞ þ L1 ðRk Þ; kj

as desired. The result holds for all j ¼ 1; :::; jGj: Case i ¼ n þ 1: Assume Equation (3) holds for i ¼ n. Consider the calculation of fnþ1 ðRj Þ for any j ¼ 1; :::; jGj. Again, error ranges can be defined and denoted by nþ1;k where k ¼ j þ 1; :::; jGj. Let hnþ1;k be the min cost to support Rj at stage n þ 1 given that R 2 nk at stage n: Note that a lower bound on the cost to support R 2 nk at stage n is

896

Donohue et al. Ln ðRk Þ by our induction assumption. Also, Cðgnþ1 ðR  Rj ÞÞ  Cðgnþ1 ðRkþ1  Rj ÞÞ holds for all R 2 nk : Therefore, hnk  Cðgnþ1 ðRkþ1  Rj ÞÞ þ Ln ðRk Þ;

and fnþ1 ðRj Þ  min Cðgnþ1 ðRkþ1  Rj ÞÞ þ Ln ðRk Þ; kj

j

follows immediately.

Li ðRj Þ is easily computed along with fid ðRj Þ as the recursion in (DP-E2) evolves. The error bound fNd ð0Þ  LN ð0Þ then provides a measure of the solution’s accuracy. If the bound is unacceptably large, the problem can be rerun with a denser grid (i.e., a larger number of grid points, jGj). A denser grid will reduce the gap ranges that are overlooked in the recursion. Obviously a tradeoff exists between the accuracy and speed of the solution Table 1. Design parameters for exponential machines, by station Cost per machine, ($) Parts/hr

ST 1

ST 2

ST 3 ST 4

30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0

1494 2116 2743 3348 4110 4372 4500 4590 6000 7000

2000 6750 3000 7000 4400 7500 4625 8000 6375 8625 7000 8750 7500 9000 7750 9625 8500 9750 10 000 10 000

300 1020 1380 1740 2100 2280 2460 2640 2820 3000

ST 5

ST 6

1120 1360 1620 1900 2200 2520 2860 3220 3600 4000

200 250 260 270 300 310 320 330 370 400

algorithm, which should be considered when choosing jGj. To demonstrate the effect of jGj on the error of the DP, consider the following six station production line design problem, with requirements HT ¼ 70 units/hour and CTT ¼ 30 minutes. Table 1 lists the possible machine options for each station, while Table 2 lists the DP algorithm’s minimum cost solution for different grid sizes. Note that because state information is passed from stage to stage (i.e., station to station) the algorithm provides a minimum cost solution for every break in the cycle time grid, basically for free. Table 2 lists the actual average flow times (T a ), minimum machine cost (Cost), and error bound (EB) for grid points between 18 and 30 minutes. Table 2 shows that jGj ¼ 3 ði:e:; G ¼ ð0; 10; 20; 30ÞÞ yields a minimum cost of fNd ð0Þ ¼ $29 696 and lower bound of LN ð0Þ ¼ $30 496, which represents an error of at most 2:7%. The arrows in columns 2–4 of Table 2 mark the feasible range for this solution, T ¼ ð25:85–30Þ. As the number of points in the grid increases, information on solutions for other ranges of T begin to emerge and the solution error rapidly decreases. The lower bound and objective function finally converge when jGj ¼ 15. The results of this example are identical to most problems tested. The algorithm locates the optimal solution with a small number of grid points, in this case 10, but takes longer to equate the objective function and lower bound. 4.2. General system A few subtle differences make the general problem significantly more difficult than the exponential case. First, while the cost to decrease R by s is independent of R, it does vary with c2di . Recall in the exponential model that sojourn time could be viewed as a range with breaks

Table 2. Performance of algorithm with various grids T

30 29 28 27 26 25 24 23 22 21 20 19 18

Grid = 3

Grid = 5

Ta

Cost ($)

EB ($)

Ta

25.85 + + + + – – – – – – – –

29 696 # # # # – – – – – – – –

800 # # # # – – – – – – – –

27.37 + + 23.15 + + + 17.83 + + + + +

Cost EB ($) ($) 29 096 200 # # # # 29 326 430 # # # # # # 30 216 1290 # # # # # # # # # #

Grid = 10 Ta 27.07 + + 25.85 + 23.15 + 20.15 + + 17.78 + +

Grid = 15

Cost EB ($) ($) 28 926 # # 29 126 # 29 326 # 29 576 # # 29 916 # #

30 # # 200 # 200 # 250 # # 340 # #

Grid = 30

Ta

Cost ($)

EB ($)

Ta

Cost ($)

EB ($)

27.07 + + 25.85 + 23.15 + 21.37 + 19.63 + 17.78 +

28 926 # # 29 126 # 29 326 # 29 376 # 29 596 # 29 916 #

0 # # 0 # 200 # 50 # 220 # 340 #

27.07 + + 25.85 + 24.37 23.15 22.63 21.37 20.15 19.63 18.99 17.78

28 926 # # 29 126 # 29 126 29 326 29 346 29 376 29 567 29 596 29 776 29 916

0 # # 0 # 0 0 20 10 0 20 160 40

897

Optimal design of stochastic production lines corresponding to changes in gi ðsÞ for 0  s  CTT . In the general model, such breaks are only determined once c2di is set. Second, the cost of moving along the second dimension of our state, c2di , depends on the absolute values of c2di1 and c2di as well as their relative difference. This means that cost breaks corresponding to changes in c2di must be examined separately for each value of c2di1 . Lastly, machine types cannot simply be ranked from fastest to slowest; we must also consider differences in process variance. This increases the numerical complexity of the problem. We attack the problem at stage i by assuming a particular arrival process is specified. Suppose that q is an index into a set V of possible values for c2di1 : A subproblem can then be developed in terms of each Vq . Let gi ðS; D; Vq Þ represent the cheapest resource combination that supports an average sojourn time of S or less and a squared coefficient of variation for the departure process of D or less. More specifically, the function gi ðS; D; Vq Þ represents the solution to the subproblem ðSP-GÞ min Cðxi ; mi Þ; Subject to sðxi ; mi ; Vq Þ  S; cdi ðxi ; mi ; Vq Þ  D; ðhsei Þ=mi  1; 8 i; h ¼ HT ; xi 2 Qi ; mi 2 Iþ ; where c2di ðxi ; Vq Þ is given by (1). By studying the structure of this subproblem, it is easy to show that fi ðRi ; c2di Þ is non-decreasing in Ri (i.e., Proposition 1 and Corollary 1 still hold in the general case). Similar results hold in terms of D, as stated below. Proposition 2. The optimal resource cost at a given station i, given a fixed arrival process, is non-increasing in D. Corollary 2. For a given station i, value R, and fixed arrival process, if fi ðR; Da Þ  fi ðR; Db Þ where Da < Db then fi ðR; Db Þ is suboptimal. Proposition 2 suggests that if we want to decrease the coefficient of variation of the departure process significantly we need to invest in better (faster and/or more precise) resources. Corollary 2 further refines the structure to suggest that for sufficiently small decreases, our choice of resource (and thus the associated cost) will not change. In other words, there exists  > 0 such that fi ðR; D  dÞ ¼ fi ðR; DÞ for 0  d  , R > 0, and D > 0. This is similar to the structure that facilitated tracking a bound on the error introduced by discretizing R in the exponential case. Additional notation is needed to formulate the discreteversion of the general dynamic program. Let V denote the grid for c2di where V ¼ ðV1 ; V2 ; . . . ; VjV j Þ. As before, we

Fig. 2. Error tracking at stage i  1 for general case.

assume the grid is identical at each stage. G will continue to denote the grid for R. Assuming the points in both grids are indexed in increasing order, problem (DP-G1) can now be redefined in discrete terms as: ðDP-G2Þf0d ðT ; c2a Þ ¼ 0 fid ðRj ; V‘ Þ ¼ min

h

a ÞÞ: minfCðgi ðRak;i1  Rj ; V‘ ; Vq;i1 kj i d þfi1 ðRk ; Vq Þg ; j ¼ 1; . . . ; jGj;

q¼1;:::;jV j

‘ ¼ 1; . . . ; jV j; i ¼ 1; . . . ; N :

where Rak;i1 again represents the actual remaining cycle a time corresponding to grid point k and Vq;i1 is the actual coefficient of variation of the arrival process corresponding to grid point q, both at stage i  1. Note that a Rak;i1  Rk while Vq;i1  Vq since cost increases with remaining flow time but decreases with the variation of the output process. By keeping track of all gaps in R and V created by the discrete approximation, we can calculate a lower bound on the original problem at each stage. Figure 2 provides an illustration. The small black boxes mark the least cost a solutions ðRak;i1 ; Vq;i1 Þ for a particular grid point. For example, the least cost solution that achieves at most VjV j and at least R0 is located in the cell bounded by VjV j , VjV j1 , R2 and R3 . Since cost increases as we move downward or to the left, this is also the least cost solution for grid points ðVjV j ; R1 Þ and ðVjV j ; R2 Þ. The non-shaded regions represent gaps where error may be introduced by the grid. The costs of points in each gap region are bounded above by the nearest solution (i.e., black square) above or to the left. We use this structure to develop a lower bound for the problem which is similar in spirit to the lower bound created in the exponential case. Let Li ðRj ; V‘ Þ denote a lower bound on fi ðRj ; V‘ Þ for a given state ðj; ‘Þ and stage i. The error bound is then computed as follows.

898

Donohue et al.

Theorem 2. The error introduced by discretizing R is bounded above by fNd ð0; c2T Þ  LN ð0; c2T Þ, where LN ð0; c2T Þ is given by the following recursion. L1 ðRj ; V‘ Þ ¼ f1d ðRj ; V‘ Þ; j ¼ 1; . . . ; jGj; l ¼ 1; . . . ; jV j; h Li ðRj ; V‘ Þ ¼ min min Cðgi ðRkþ1  Rj ; V‘ ; Vq1 ÞÞ q kj

i þ Li1 ðRk ; Vq Þ ; j ¼ 1; . . . ; jGj; l ¼ 1; . . . ; jV j; i ¼ 2; . . . ; N :

ð4Þ

b2kq ðj; lÞ ¼ Cðg2 ðRkþ1  Rj ; Vl ; Vq ÞÞ þ f1d ðVq ; Rk Þ: Finally, considering scenario (iii), let b3kq ðj; lÞ denote a lower bound on the minimum cost to support Rj and Vl at stage 2 given v 2 V1q and r 2 R1k ðvÞ. Then we have, b3kq ðj; lÞ ¼ Cðg2 ðRkþ1  Rj ; Vl ; Vq1 ÞÞ þ f1d ðVq ; Rk Þ; ¼ Cðg2 ðRkþ1  Rj ; Vl ; Vq1 ÞÞ þ L1 ðRk ; Vq Þ:

ð5Þ ð6Þ

Proof. The proof involves extending the logic used in the proof of Theorem 1 to two-dimensions. Equation (4) follows immediately, as before, since no error is passed from stage 0 to stage 1. Equation (6) is, once again, shown by induction on i. Case i = 2: Consider the calculation of f2d ðRj ; Vl Þ for any j ¼ 1; . . . ; jGj and l ¼ 1; . . . ; jV j: This value may differ from f2 ðRj ; Vl Þ because values in ranges ð0; V1a Þ; . . . ; ðVjV j1 ; VjVa j Þ are not considered and for each grid point Vs , ranges ðRa1 ðVs Þ; R2 ðVs Þ; . . . ; Raj1 ðVs Þ; Rj ðVs ÞÞ are not considered. Let V1q denote a specific gap ðVq1 ; Vqa Þ in V at stage 1 and R1k ðVqa Þ denote a specific gap ðRak ; Rkþ1 Þ in R at stage 1 given Vqa . Three different error scenarios must be considered: (i) a solution is overlooked involving v 2 V1q where q ¼ 1; . . . ; jV j; (ii) given some CV partition value Vqa , a solution is overlooked involving r 2 R1k ðVq Þ where k ¼ 1; . . . ; jQj; and (iii) the two previous scenarios combined, namely a solution is overlooked involving v 2 V1q where q ¼ 1; . . . ; jV j and r 2 R1k ðvÞ. First, consider scenario (i). Note that

It is easy to show that b3kq ðj; lÞ provides a lower bound on both b1kq ðj; lÞ and b2kq ðj; lÞ. Taking the minimum of b3kq ðj; lÞ over all possible q and k then gives  

f2 ðRj ; Vl Þ  min min b3kq ðj; lÞ ; q kj  ¼ min min Cðg2 ðRkþ1  Rj ; Vl ; Vq1 ÞÞ q kj 

þL1 ðRk ; Vq Þ ; as desired. Case i = n+1: The proof is straightforward, combining the logic of the previous case with case i ¼ n þ 1 of Theorem 1, and is therefore omitted. j Li ðRj ; V‘ Þ is easily computed along with fid ðRj ; V‘ Þ as the recursion for the discrete problem evolves. The complexity of the general dynamic program is OðN jQjjGj2 jV j2 Þ. In practice, the algorithm usually runs much faster for two reasons. First, only a fraction of the grid may be feasible at a given stage because of resource and cost restrictions. Second, the solution for a given feasible grid point often dominates higher grid points.

f1 ðv; Rk Þ  f1 ðVqa ; Rk Þ; where v 2 V1q since v  Vqa and the cost function f increases as the coefficient of variation improves (i.e., as V decreases). Furthermore, Cðg2 ðRk  Rj ; Vl ; vÞÞ  Cðg2 ðRk  Rj ; Vl ; Vq1 ÞÞ; since Vq1  v and g2 and C are non-decreasing in V and g2 , respectively. Therefore, if we let b1kq ðj; lÞ denote a lower bound on the minimum cost to support Rj and Vl at stage 2 given v 2 V1q then b1kq ðj; lÞ  Cðg2 ðRk  Rj ; Vl ; Vq1 ÞÞ þ f1d ðVq ; Rk Þ: Now consider scenario (ii) which is similar to errors arising in the exponential case. Note again that f1 ðVq ; rÞ  f1d ðVq ; Rk Þ; and Cðg2 ðr  Rj ; Vl ; Vq ÞÞ  Cðg2 ðRkþ1  Rj ; Vl ; Vq ÞÞ: Therefore, if we now let b2kq ðj; lÞ denote a lower bound on the minimum cost to support Rj and Vl at stage 2 given r 2 R1k then

5. Numeric example To illustrate how the general algorithm could be used in practice, consider the problem of selecting machines for the Protective Coating process of a Printed Circuit Board (PCB) plant. This process takes place near the end of the PCB line, before panels are sized and packaged. The function of a Protective Coating is to apply a protective epoxy resin coating to the surface of each panel to guard against damage to the circuitry. Figure 3 illustrates the steps of the Protective Coating Process while Table 3 offers a list of possible resource options, for each station, based on industry data. 5.1. Base case and sensitivity analysis We wish to find the lowest cost design that supports an average throughput rate of HT ¼ 3000 panels/day, CTT ¼ 8 hours, and c2dN  c2a ¼ 1. The process runs three shifts per day with roughly six and a half productive

899

Optimal design of stochastic production lines

Fig. 3. The protective coating process (Cost in 1000’s).

Table 3. Design parameters for sample machines and resulting solutions Design parameters Tool Coat

C1 C2 C3 C4 Expose E1 E2 E3 E4 E5 E6 E7 E8 Develop D1 D2 D3 D4 D5 Inspect I1 I2 Bake B1 B2 B3 B4 B5 Final inspection F1 F2 F3 Total cost (·$1000) a

Cost ($1000) 1500 1600 850 1100 100.0 120.0 133.0 140.0 102.0 122.0 135.0 142.0 150.0 160.0 165.0 90.0 190.0 40.5 40.0 400 410 460 450 500 40 32 37

Time per job (56 panels/job);

b

sa

c2

MTBF

1.0 1.0 2.0 2.0 100.0 120.0 90.0 100.0 100.0 120.0 90.0 100.0 1.0 1.0 0.75 2.0 1.0 15.0 20.0 2.0 2.0 2.0 1.0 1.0 150.0 180.0 165.0

0.25 0.25 0.25 0.25 0.5 0.25 1.0 0.25 0.5 0.25 1.0 0.25 0.25 0.25 0.25 0.25 0.25 1.0 0.5 0.25 0.25 0.25 0.25 0.25 1.0 1.0 0.5

4800 4800 15 360 21 120 16 770 16 770 16 770 16 770 16 770 16 770 16 770 16 770 18 828 18 828 18 828 18 828 27 028 – – 18 828 18 828 24 828 4800 4800 – – –

Mean Time to Failure;

c

b

Solutions (mi , sojourn time) MTTR

c

240 120 360 420 498 498 498 498 258 258 258 258 372 192 372 192 192 – – 492 372 192 252 252 – – –

Mean Time to Repair;

conv

d

45 45 45 45 – – – – – – – – 4.0 4.0 4.0 8.0 6.0 – – 100.0 100.0 95.0 100.0 95.0 – – –

d

Setup – – – – 15 15 15 15 6 6 6 6 – – – – – – – – – – – – – – –

e

Prob1

Prob2

Prob3

1, 47.60

1, 47.60

1, 47.65

6, 117.53

6, 125.39 7, 100.56

1, 10.21

1, 10.20

1, 10.23

1, 27.85

1, 27.85

1, 102.67

1, 102.67

2, 21.8 1, 102.73

8, 173.17

8, 173.17

9, 164.67

2312.5

2645.5

2392.5

Conveyor traversal; and e Setup time between jobs.

900

Donohue et al.

Table 4. Optimal resource costs ($1000) for various values on CTT and c2T Bound on squared coefficient of variation of output process (c2T )

CTT 0.7 500 495 490 485 480 475 470 465 460 455 450 445 440 435 430 425 420

2312.5 (24) + # # \\\\\\\\\\\\\\ 2352.0 (41) # 2352.5 (30) # 2392.0 (41) 2453.0 (61) 2513.0 (120) 2553.0 (101) 2633.0 (165) 2695.0 (187) 2791.0 (223) 2926.0 (220)

0.6 2322.5 (32) # 2325.5 (35.5) 2352.5 (30) + # # # # ! ! ! ) ) ) ) )

0.5 2352.5 + # 2362.5 + # 2392.5 # # 2432.0 2433.0 2493.0 ! ! ! ! !

0.4 (51.5)

(32.5)

(31.5)

(41) (80.5) (100)

hours per shift. To determine the minimum cost solution, we considered average flow times of up to 500 minutes (8.3 hours) using grids of jGj ¼ 500 and V ¼ ð0:0; 0:1; 0:2; . . . ; 1:1; 5:0; 10:0Þ. The minimum cost machine configuration identified by our algorithm is listed in Table 3 (Prob1). The associated cost is $2312:5K with an error bounded above by $22:5K. Table 4 lists the minimum costs identified by the algorithm for each grid range. The shaded area denotes the solution to Prob 1 (i.e., satisfying CTT  8 hours and c2T  1). The numbers in parenthesis denote bounds on the error introduced by the two-dimensional grid. Error bounds range from 0:04 to 9:4%, which is typical of most problems we tested. One can perform sensitivity analysis on capability targets CTT and c2T by moving in various directions within the grid. For example, loosening the average flow time bound (i.e., moving upward from the shaded box) by 20 minutes does not affect total cost. However, tightening the average flow time target (i.e., moving downward) by 10 or 20 minutes results in a cost increase of $39:5K and $79:5K, respectively. Tightening the output coefficient of variation bound (i.e., moving to the right) from 0.7 to 0.6 increases cost by $39:5K. One can use this information to weigh the costs and benefits of changing CTT and c2T either separately or simultaneously. Performing sensitivity analysis on HT requires re-running the algorithm with a new input rate. Table 5 displays the effect of increasing throughput by varying degrees. This table is discussed in more detail in the next subsection. One can easily modify the general algorithm to include station-specific constraints to impose more control on the final solution. For example, consider the unique problems the Expose work center experiences in its clean room environment. Process Control (PC) managers often

2411 (37) 2412.5 (38.5) # 2436.5 (26.5) + # 2472.5 (37.5) + # 2512.0 (41) # 2513.0 (40.5) 2573.0 (61) 2633.0 (120) 2735.0 (162) 2831.0 (198) 2966.0 (278)

0.3 2491 (43) 2502.5 (53.5) # 2510.5 (20.5) + # 2550.0 (39.5) # 2551.0 (40.5) 2592.0 (42) # 2593.0 (42) 2653.0 (61) 2713.0 (120) 2815.0 (147) 2911.0 (183) 3046.0 (180)

0.2 2720.5 2732.5 + # # # 2772.0 # 2773.0 2832.0 # 2833.0 2893.0 2953.0 3055.0 3151.0 3286.0

(50.5) (2.5)

(39.5) (40.5) (60) (40.5) (61) (120) (147) (183) (180)

monitor the number of jobs maintained at such operations to avoid the possibility of contamination. As we saw in Table 3, our current solution (Prob1) maintains an average of six jobs at the Expose work center, with an average queueing delay of 117.53 minutes. If the PC manager feels this inventory level is too high, he can restrict it by adding a constraint on the average number of jobs maintained at the Expose work center. Table 3 shows how the solution changes when the average number of jobs at Expose is limited to 4.5 (Prob2). This additional constraint increases resource costs by 14.4%. 5.2. Adding a service level constraint The capacity solutions discussed thus far support production of 3000 panels per day on average (i.e., HT ¼ 30 000). This means the line may fall short of a 3000 panel daily quota 50% of the time. A service level constraint offers a way to improve the chances of meeting this production quota. In order to do this, we need to characterize the distribution of output over some time period P (e.g., a day or week). Recall that the queueing approximations used throughout this paper are based on the decomposition approach which assumes the departure process at each station is approximately renewal. Continuing with this assumption, if the output process is a renewal process with rate hN , and squared coefficient of variation c2dN , then the mean and variance of the number of outputs MP over time period ½0; P  is expressed as E½MP  ¼ P hN ; ð7Þ

2  Var½MP  ¼ P cdN hN   3 5 4 2 3 2 þ 2cdN þ þ cdN  l3 hN þ oð1Þ ; ð8Þ 4 4 3

901

Optimal design of stochastic production lines Table 5. Optimal resource cost for various input rates and c2dN ’s, assuming CTT ¼ 8 hours 2 CdN

Average daily throughput (panels/day) 3000

1.0

0.9

0.8

0.7

0.6

0.5

0.4

0.3

(3000,50%) (2728,95%) $2312.5 (24) (3000,50%) (2742,95%) + (3000,50%) (2757,95%) # (3000,50%) (2773,95%) $2312.5 (24) (3000,50%) (2789,95%) $2352.5 (30) (3000,50%) (2808,95%) $2392.5 (32.5) (3000,50%) (2828,95%) $2436.5 (26.5) (3000,50%) (2851,95%) $2510.5(20.5)

3050 (3000,61.9%) (2778,95.0%) $2352.0 (39.5) (3000,62.5%) (2792,95.0%) + (3000,63.2%) (2807,95.0%) # (3000,64.1%) (2822,95.0%) $2352.0 (39.5) (3000,65.2%) (2839,95.0%) $2352.5 (39) (3000,66.5%) (2858,95.0%) $2362.5 (22) (3000,68.4%) (2878,95.0%) $2436.5 (36) (3000,70.9%) (2900,95.0%) $2543.0 (40.5)

3100 (3000,72.7%) (2828,95.0%) $2352.5 (39) (3000,73.8%) (2842,95.0%) + (3000,75.0%) (2857,95.0%) # (3000,76.5%) (2872,95.0%) # (3000,78.2%) (2889,95.0%) $2352.5 (39) (3000,80.4%) (2908,95.0%) $2362.5 (22) (3000,83.0%) (2928,95.0%) $2436.5 (22) (3000,86.5%) (2951,95.0%) $2543.0 (40.5)

3150 (3000,81.8%) (2878,95.0%) $2352.5 (40) (3000,83.0%) (2892,95.0%) + (3000,84.4%) (2907,95.0%) $2352.5 (40) (3000,86.1%) (2922,95.0%) $2352.5 (20) (3000,87.9%) (2939,95.0%) $2366.0 (40.5) (3000,90.0%) (2958,95.0%) $2399.5(17) (3000,92.2%) (2978,95.0%) $2436.5 (36) (3000,95.1%) (3000,95.0%) $2547.5 (27)

(Heyman and Sobel, 1982) where l3 is the third moment of the time between outputs. These distribution measurements are useful for calculating the probability, s, of meeting some specified target production quota, Q, over the time frame P . For example, if the distribution of MP is normal (a plausible assumption when P is large) then the relation between (7), (8), Q, and s is described by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi E½MP  þ zð1sÞ Var½MP   Q: ð9Þ In some cases, a minimum service level s may be specified a priori based on an analysis of the amount of overtime or late deliveries the system is willing to incur (Hopp et al., 1993). Equation (9) then determines the maximum production quota the production line can support at this service level. Now suppose we wish to support a production quota of 3000 panels per day with a service level of 95% while still maintaining an average cycle time of 8 hours. To improve the chances of meeting quota, we must increase the average throughput target, decrease the output variation target, or perform some combination of the two. Table 5 displays a mapping between process output characteristics and service level for the Protective Coating process, using Equation (9). The first two lines within each cell list feasible ðQ; sÞ combinations, first when Q is set at 3000 and second when s is set at 95%. The third line lists the optimal capacity cost determined by our algo-

3200 (3000,88.7%) (2928,95.0%) $2352.5 (52) (3000,89.9%) (2942,95.0%) + (3000,91.2%) (2957,95.0%) # (3000,92.6%) (2972,95.0%) $2352.5 (42) (3000,94.1%) (2989,95.0%) $2392.5 (52) (3000,95.6%) (3008,95.0%) $2432.5 (62) (3000,97.2%) (3028,95.0%) $2472.5 (48) (3000,98.6%) (3051,95.0%) $2547.5 (27)

3250 (3000,93.5%) (2978,95.0%) $2392.0 (52) (3000,94.4%) (2992,95.0%) + (3000,95.4%) (3007,95.0%) # (3000,96.4%) (3022,95.0%) $2392.0 (52) (3000,97.4%) (3039,95.0%) $2392.5 (52) (3000,98.4%) (3058,95.0%) $2432.5 (62) (3000,99.2%) (3078,95.0%) $2512.5(82) (3078,99.7%) (3101,95.0%) $2587.0 (66.5)

3300 (3000,96.5%) (3028,96.0%) $2392.0 (39.5) (3000,97.2%) (3042,95.0%) + (3000,97.9%) (3057,95.0%) # (3000,98.5%) (3072,95.0%) $2392.0 (39.5) (3000,99.0%) (3089,95.0%) $2392.5 (40) (3000,99.5%) (3108,95.0%) $2432.5 (50) (3000,99.8%) (3128,95.0%) $2512:5ð70Þ (3000,99.9%) (3151,95.0%) $2592.5 (60)

rithm. The thick zigzagged line running through the table separates the infeasible region on the left from the smaller feasible region on the right. We see that an average throughput rate of at least 3150 panels/day is needed to achieve the 95% service level. At this rate, output variability must be held closely in check (c2T ¼ 0:3). The total resource cost is $2547:5K, a 10.2% increase over the base case (i.e., (Prob1)). One can do better by increasing the average throughput rate to 3250 panels/day. This faster rate allows a higher level of variability in the output process (c2T ¼ 0:7). In turns out that controlling process variability is more expensive than increasing process speed given the resource choices available. Total resource costs reduce to $2392:0K under this option, a minor 3.4% increase over the base case. The machine configuration for this option is listed as (Prob3) in Table 3.

6. Conclusions and extensions This paper describes an efficient approach for finding the least cost resource configuration of a production line that attains specified throughput and cycle time constraints. Since our approach leads to an easily computed bound on the error that results from discretizing the state space for computational purposes, we can ensure that the optimization procedure does not introduce large amounts of

902 error. The procedure is essentially as accurate as the queueing approximations upon which it is based. Unless the squared coefficients of variation are excessively high (e.g., greater than five), these approximations should be at least as accurate as the input data. For such cases, this approach is a reasonable basis for a practical line design tool. Our analysis assumed that resource costs were independent and the sequence of operations was fixed. However, our approach could be extended to the case where resource costs are dependent (e.g., there are discounts for buying equipment from the same supplier). Doing this would require augmenting the action spaces to include ‘‘bulk buying’’ (i.e., purchasing sets of equipment at the discount price). Although this would increase the complexity of the problem, the basic structure would remain the same. Similarly, we could extend our results to the case where equipment choices affect the sequence (e.g., equipment exists that could accomplish multiple tasks). To do this, we would need to run our algorithm for different line lengths (e.g., to consider machine combinations for performing tasks 1–6 in six steps, five steps, four steps, etc.). Another, heuristic, approach is suggested by Askin and Zhou (1998) who study the deterministic version of this machine sequencing problem. After solving the problem to optimality using a shortest path algorithm, they propose a simple greedy heuristic which extends the state space at each stage to include the resource used at the previous stage, as well as the resource having the minimum incremental cost. This heuristic could easily be applied within our stochastic DP framework. Our work also focused exclusively on the single product case. Multiple product production lines could be handled, provided the products have the same, or almost the same, routing. Where this is the case, the processing time means and coefficients of variations can be combined to form a single representative job class (Buzacott and Shanthikumar, 1993). Then, the optimization procedure described in this paper could be applied directly. If the routings differ but are ‘‘unidirectional,’’ in the sense that no two work centers are visited in reverse order by different products, then our approach will still work with only minor modifications. However, if the multiple products have non-unidirectional routings then, even though interference models exist for predicting the cycle times of a given configuration (Bitran and Tirupati, 1988; Whitt, 1994), the sequential optimization approach required by the DP algorithm cannot be used. This is because we must solve the traffic equations and the variability traffic equations (i.e., the equations for the scv’s of the departure processes) for the underlying queueing network in stage-wise fashion. As we have noted, this is the case for tandem production lines, of which there are many in industry (including many of the printed circuit board lines we

Donohue et al. observed). But, for lines with feedback (e.g., rework or routings that pass through some resources more than once) or multiple product lines with non-unidirectional routings this does not hold and hence our method cannot be used directly. Further research is needed to devise an appropriate optimization method for such situations. In closing, it is worth mentioning that the solution procedure for any manufacturing design problem is dependent on the availability of accurate machine-specific information (e.g., expected process time, mean times to failure, etc.). Freit and Wu (2000) recently developed an interesting framework for dealing with situations where this information is uncertain or can only be gathered at significant cost. Further work is needed to extend these ideas to the constrained optimization problem studied here.

Acknowledgement This work was supported in part by The Center for Manufacturing and Logistics Research, Motorola, Inc., and the National Science Foundation under Grant DMI9322830.

References Askin, R.G. and Zhou, M. (1998) Formulation of independent flowline cells based on operation requirements and machine capabilities. IIE Transactions, 30, 319–329. Bitran, G. and Sarkar, D. (1994) Targeting problems in manufacturing queueing networks – an interative scheme and convergence. European Journal of Operational Research, 76, 501–510. Bitran, G. and Tirupati, D. (1988) Multiproduct queueing networks with deterministic routing: decomposition approach and the notion of interference. Management Science, 34(1), 75–100. Bitran, G. and Tirupati, D. (1989a) Trade-off curves, targeting and balancing in manufacturing queueing networks. Operations Research, 37(4), 547–564. Bitran, G. and Tirupati, D. (1989b) Capacity planning in manufacturing networks with discrete options. Annals of Operations Research, 17, 119–136. Boxma, O.J., Rinnooy Kan, A.H.G. and Van Vliet, M. (1990) Machine allocation problems in manufacturing networks. European Journal of Operational Research, 45, 47–54. Buzacott, J. and Shanthikumar, J.G. (1993) Stochastic Models of Manufacturing Systems, Prentice Hall, Englewood Cliffs, NJ. Connors, D., Feigin, G. and Yao, D. (1993) A queueing network model for semiconductor manufacturing. Working paper, IBM Research Division, T.J. Watson Research Center, Yorktown Heights, NY. Donohue, K. (1994) The economics of capacity and marketing measures in a simple manufacturing environment. Production and Operations Management, 3, 78–99. Feit, E.M. and Wu, S.D. (2000) Transfer line design with uncertain machine performance information. IEEE Transactions on Robotics and Automation, 16(5), 581–587. Frenk, H., Labbe, M., Van Vliet, M. and Zhang, S. (1994) Improved algorithms for machine allocation in manufacturing systems. Operations Research, 42, 523–530.

903

Optimal design of stochastic production lines Heyman, D. and Sobel, M. (1984) Stochastic Models in Operations Research, McGraw-Hill, New York, NY. Hopp, W. and Spearman, M. (2000) Factory Physics: Foundations of Manufacturing Management, 2nd edn, Irwin/McGraw-Hill, Burr Ridge, IL. Hopp, W., Duenyas, I. and Spearman, M. (1993) Economic production quotas for pull manufacturing systems. IIE Transactions, 25, 71–79. Kouvelis, P. and Tirupati, D. (1991) Approximate performance modeling and decision making for manufacturing systems: a queueing network optimization framework. Journal of Intelligent Manufacturing, 2, 107–134. Lawrence, S.R. and Buss, A.H. (1995) Economic analysis of production bottlenecks. Mathematical Problems in Engineering, 1, 341– 363. Schweitzer, P.J. and Seidmann, A. (1991) Optimizing processing rates for flexible manufacturing systems. Management Science, 37, 454– 466. Van Vliet, M. and Rinnooy Kan, A. (1991) Machine allocation algorithms for job shop manufacturing. Journal of Intelligent Manufacturing, 2, 83–94. Walrand, J. (1993) Queueing networks, in Handbooks in Operations Research and Management Science, Volume 2: Stochastic Models, Heyman, D.P. and Sobel, M.J. (eds), North Holland, New York, pp. 519–603. Whitt, W. (1993) Approximations for the GI/G/m queue. Production and Operations Management, 2, 114–161. Whitt, W. (1994) Towards better multi-class parametric-decomposition approximations for open queueing networks. Annals of Operations Research, 48, 221–248.

Biographies Karen Donohue is an Assistant Professor of Operations and Management Science at the University of Minnesota. She joined the Carlson School in 2000, after serving 6 years on the faculty of the Wharton School at the University of Pennsylvania. She earned a Ph.D. in Industrial Engineering from Northwestern University in 1993. She also holds an M.S. in Industrial Engineering and Management Science from Northwestern University and a B.A. in Mathematics and Economics from St. Olaf College. Professor Donohue’s research interests focus on the design and analysis of manufacturing systems and mechanisms for coordinating production, distribution, and sales decisions across a supply chain. Previous research projects include: (i) developing a new generation of facility design tools for the planning of semiconductor fabrication facilities; (ii) quantifying the effects of return contracts on production and inventory replenishment decisions; and (iii) investigating the possible operational benefits of dual distribution channels for short life cycle products. Her research has been supported by NSF, Sandia National Laboratory, the US Navy, and the Sloan Foundation Industry Study on Retailing. Professor Donohue serves on the Editorial Boards of the journals of Manufacturing and Service Operations

Management (M&SOM) and Production and Operations Management (POMS). Wallace J. Hopp is the Breed University Professor of Manufacturing Management in the Department of Industrial Engineering and Director of the Master of Management and Manufacturing Program, a joint program of the Kellogg Graduate School of Management and the McCormick School of Engineering, at Northwestern University. He has been at Northwestern since receiving his Ph.D. in Industrial and Operations Engineering from the University of Michigan in 1984. His research focuses primarily on the design, control and management of manufacturing systems. He has won a number of research and teaching awards, including the 1985 Nicholson Prize, the 1989 McCormick Teaching Award, the 1990 Scaife Award (with Mark Spearman), the Pentair-Nugent Professorship in Business Leadership in 1993, the Kellogg Top Five Professors Award in 1998, the 1998 IIE Joint Publishers Book-of-the-Year Award (for the book Factory Physics: Foundations of Manufacturing Management), and the 2001 Sargent Americanism Award from SME. He is Department Editor for Manufacturing, Distribution and Service Operations for the journal Management Science, a senior editor of M&SOM, an Associate Editor of IIE Transactions, and a member of IIE, INFORMS, POMS, and SME. He is an active industry consultant, whose clients have included Anixter, Bell & Howell, Black & Decker, Case, Dell, Ford, Eli Lilly, General Motors, John Deere, IBM, Intel, Motorola, Owens Corning, S&C Electric, SPX, Texas Instruments, Whirlpool, Zenith, and others. Mark L. Spearman is President and Chief Executive Officer of Factory Physics, Inc. and is on leave from Georgia Tech where he is Professor of Industrial and Systems Engineering. His research and teaching involve manufacturing enterprise integration, and supply chain management. Prior to joining Georgia Tech and founding Factory Physics, Inc. Spearman was a tenured professor in the Department of Industrial Engineering and Management Sciences at Northwestern University. Before he obtained his Ph.D., he worked over 5 years as an engineer for Monsanto Corporation, Superior Oil Company, and IBM Corporation. Spearman has been the principal investigator for $1:75 million and has acted as a contributor for $3 million in research support from both government and industrial sources and has published numerous articles in the area of production and operations management and industrial engineering. He is co-author, with Wallace J. Hopp, of the book, Factory Physics: The Foundations of Manufacturing Management, that has been named the IIE Book of the Year (1998). He is PastPresident of the Manufacturing and Services Operations Management Society and has served as Secretary and a Member of the Board of the Production and Operations Management Society. He is an Associate Editor for Management Science and has served as Area Editor for Production and Operations Management and as Senior Editor for Manufacturing and Service Operations Management. He holds a Ph.D. in Industrial Engineering, is a Senior Member of IIE, a full member of INFORMS, and is Certified in Production and Inventory Management by APICS. Contributed by the Manufacturing Systems Modelling Department