Performance evaluation for general queueing networks in ...

Author's personal copy

European Journal of Operational Research 221 (2012) 328–339

Contents lists available at SciVerse ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Production, Manufacturing and Logistics

Performance evaluation for general queueing networks in manufacturing systems: Characterizing the trade-off between queue time and utilization Kan Wu a,⇑, Leon McGinnis b a b

School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Ave., Singapore 639798, Singapore School of Industrial and Systems Engineering, Georgia Institute of Technology, USA

a r t i c l e

i n f o

Article history: Received 8 November 2011 Accepted 8 March 2012 Available online 17 March 2012 Keywords: Manufacturing Production Productivity and competitiveness Queueing

a b s t r a c t Performance evaluation plays a key role in manufacturing system design and productivity improvement. Characterizing performance objectively is the first step. Inspired by the underlying structure of tandem queues, we have derived an approximate model to characterize the system performance. The model decomposes system queue time and variability into bottleneck and non-bottleneck parts while capturing the dependence among workstations. Compared the new model with prior approaches, the new model not only is more accurate but also requires less information. The property of manufacturing system performance is given based on the insight from the model. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction A practical manufacturing system usually faces batches, assembly, interruptions, reentry, product mixes, and complicated dispatching disciplines. Due to the complexity of the system, it is difficult to describe its behavior exactly. Instead of finding out how a manufacturing system behaves by considering every possible detail, it might be more practical to model a manufacturing system from a macroscopic view while capturing its main underlying structure. An objective evaluation of the system performance is an essential input for manufacturing system design and productivity improvement. To quantify system performance, performance curves (and their associated variabilities) play a key role since they characterize the trade-off between queue time and utilization. Because the curve characterizes the performance of a manufacturing system, it is also called a characteristic curve, trade-off curve, operation curve or queueing curve. Bitran and Tirupati (1989) used performance curves to describe the relationship between work-inprocess (WIP), cycle time and capacity. Sattler (1996) used performance curves to determine productivity improvements of a semiconductor fab. She assumed variability is independent of utilization and approximated the curve by using a constant k to replace the variability term in Kingman’s G/G/1 approximation (Kingman, 1965). Boebel and Ruelle (1996), Collins et al. (1997), Fowler et al. (1997), Ruelle (1997), and Rose (2001) used this curve to quantify the productivity improvement of a fab. ⇑ Corresponding author. Tel.: +65 6790 5584; fax: +65 6792 4062. E-mail address: [email protected] (K. Wu). 0377-2217/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ejor.2012.03.019

Since performance curves are commonly used to quantify the trade-off between queue time and utilization, it would be nice if we can generate the curves analytically. However, due to the non-renewal departure process, queueing networks can be analyzed exactly only for some special Markovian cases. Landmark papers in this regard are Jackson (1957) and Baskett et al. (1975). Hence, approximations are commonly used for the analysis of general queueing networks. Most of the approximation methods can be classified as either aggregation or decomposition approaches (Li et al., 2009). Kerbache and Smith (1987) developed the generalized expansion method (GEM) to model a general open finite queueing network, and applied the GEM to examine the optimal routing in layout and location problems (Kerbache and Smith, 2000). Due to the complexity of practical manufacturing systems, it is difficult to take all the details into account. Hence, people sometimes model a manufacturing system by the aggregation approaches with a macroscopic view. Chandy et al. (1975) applied parametric analysis to queueing networks. They constructed an equivalent network in which all queues except those in a subsystem were replaced by a single composite queue. The behavior of the subsystem was then studied through the equivalent network. The subsystem in the equivalent network was also called the flow equivalent server. The authors showed that the behavior of the subsystem in their constructed equivalent network is the same as in the given network for a limited class of systems. Rather than creating the curve analytically, Rose (2001), Nazzal and Mollaghasemi (2001), and Park et al. (2001) developed the fab performance curve through simulations. Because simulation studies are time consuming, Sattler (1996) and Collins et al. (1997)

Author's personal copy

329

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

attempted to describe the performance curve of a fab based on Kingman’s G/G/1 approximation and demonstrated how to enhance productivity through variability reduction. However, rather than a single server, a fab is composed of a sequence of operations executed by a series of workstations. Predicting the performance curve simply based on the G/G/1 queue is not adequate. When the buffer sizes are limited, Cruz et al. (2010) solved a buffer allocation and throughput trade-off problem by means of multi-objective genetic algorithms and GEM. Under the assumption of Poisson arrivals and exponential service times, Hulett and Damodaran (2011) derived analytical approximations to predict the performance curves of manufacturing systems with parallel processing. Motivated by Kingman’s approximation, Yang et al. (2007) proposed empirical algorithms to generate performance curves through simulation and metamodeling. Two unknown scalars and one unknown vector have to be determined through complex procedures. The algorithm performs well for M/M/1 systems with various dispatching rules, and is extended to predict the performance curve of manufacturing systems by ignoring the dependence among workstations. In practice, congestion at a workstation usually implies congestion at its downstream stations in a later period. Dependence among workstations plays an important role in understanding the behavior of manufacturing systems. Furthermore, there can be reentry, rework, batches, shift schedules, and multiple products with different priorities. All these factors make modeling the dependence among workstation even harder. A workstation may also consist of heterogeneous servers. Harrison (1998) studied a specific type N-Model, which consists of two customer classes and two heterogeneous servers working in parallel. Gurvich and Whitt (2009) proposed a routing rule called queue-and-idlenessratio (QIR). Similar to the study of Harrison and Lopez (1999), Gurvich and Whitt started their study with solving a linear program to insure that the parallel server systems achieve heavy traffic under QIR. However, not all servers in a heterogeneous server system can achieve heavy traffic simultaneously in general. When some servers are not in the heavy traffic regime, queue time will also depend on dispatching rules and as a result, the queue time approximation becomes very difficult. Hence, we assume all workstations have homogeneous servers in the derivation of our approximate models. It is difficult to compute the performance curve of a practical manufacturing system exactly. Our goal is to have a good approximation of the performance curve by capturing the main underlying structure of a manufacturing system. Rather than assuming stochastic independence, we can capture the dependence among workstations by the intrinsic ratio discovered by Wu (2009). Based on the intrinsic ratio, an approximate model is derived, and it performs very well in the examined cases with complex routing, scheduling, and batching. Our model takes a ‘‘black box’’ approach which treats the whole system as a black box, where only system throughput rate and mean cycle time are known. The objective of the black box approach is to find a model which can describe factory behavior through regression analysis so that we can use only a few known throughput rates and mean queue times to fit the model parameters. After the parameters are determined, the model can be used to predict mean cycle time at any other utilization. The derivation of the black box approach is given in Section 2. Model verification is given in Section 3. Conclusions are given in Section 4.

2. Performance evaluation for systems with single-server bottlenecks Due to the variant machine cost, there is usually a distinct bottleneck in a manufacturing system, which has relatively low capacity (due to the higher cost). The bottleneck plays a key role in

system performance especially in heavy traffic. To clearly present the derivation and insight of our models, we begin with a simple system where the system bottleneck is a single server. The model will be extended to a multiple server bottleneck afterwards. 2.1. Models for multiple single-server stations in series Wu and McGinnis (2012) identified the nice property of intrinsic ratio in tandem queues, and then extended the results to n single-server stations in series with infinite buffers. Their model captures the dependence among tandem queues through so called intrinsic ratios and considerably outperforms the previous approximate models based on parametric-decomposition or diffusion approximations. Based on the property of intrinsic ratio, Procedure 1 explains how to approximate the mean queue time of each station in n single-server stations in series. It consists of two stages: decomposition and computation. Stage 1 identifies the main system bottleneck first, and then identifies the next bottleneck within each subsystem. A subsystem is composed of the stations from the first station to the newest identified bottleneck (not included). At the beginning, when no bottleneck has been identified, the subsystem is the same as the original system. The subsystem then gradually becomes smaller until it is composed of solely one single station, which is the first station of the tandem queue. Procedure 1. (Queue Time for Each Station in Serial Queues) Stage I: Decomposition by bottlenecks 1. Identify the index of the system bottleneck station (BN1), where lBN1 ¼ min li , for i = 1 to n. Let k = 1. – If more than one station has the minimum service rate, BN1 = min i, where li ¼ lBN1 . 2. Identify the index of the next bottleneck station (BNk+1) in front of the previous one (BNk), where lBNkþ1 ¼ min li, for i = 1 to BNk  1. – If more than one station has the minimum service rate, BNk+1 = min i, where li ¼ lBNkþ1 . 3. If BNk+1 = 1, stop. Otherwise, let k = k + 1, go to 2. Stage II: Computing the mean queue time for each station 4. Let QT 1 ¼ a1



q1

1q1



1

l1 ;

and i = 2, where a1 ¼



c2a1 þc2s1 2

 .

5. If station i is marked as a bottleneck, then QT i ¼   P ai 1qiqi l1i  ð1  yi Þ i1 k¼1 QT k ;   c2 þc2  Otherwise, QT i ¼ xi ai 1qiq l1 ; where ai ¼ a12 si . i

i

6. If i = n, stop. Otherwise, let i = i + 1, go to 5. Fig. 1 shows the bottleneck decomposition result. In Procedure 1, yi is the intrinsic ratio of a bottleneck in a subsystem, and xi is the intrinsic ratio of a non-bottleneck. Wu and McGinnis (2012) discovered that in tandem queues, the intrinsic ratio behaves approximately linearly. The term, ai ðqi =ð1  qi ÞÞ=li , in step 5 is called ASIA (i.e., all see the initial arrival process) system queue time of station i, since the arrival process squared coefficient of variation (SCV) in is always the SCV of the initial external arrival process. Wu (2009) gave more insights on ASIA system queue time. Summing up the mean queue times of all stations, system mean queue time can be approximated by n P i¼1

QT i ¼ f1 a1













q1 1 q2 1 qn 1 þ f2 a2 þ   þ fn an ; 1  q1 l1 1  q 2 l2 1  qn ln ð1Þ

Author's personal copy

330

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

(λ, Ca1)

(μ1, Cs1) (λ, Cd1) … 1 BNp

((μl, Csl)

(λ, Cdl)

l

(μm, Csm) (λ, Cdm) … m

…

BN2

(μn, Csn) n

BN1

Fig. 1. n single queues in series.

where QTi is the mean queue time of station i in steady states, qi is the utilization of station i, li is the service rate of station i, and ai is the variability of station i in the ASIA system. fi is called contribution factor, since it represents the percentage of station i’s ASIA system queue time contributing to the overall system queue time as shown in Eq. (1). fi can be determined based on Procedure 2. Procedure 2. (System Queue Time for Single Queues in Series): Stage I: Decomposition by bottlenecks 1. Identify the system bottleneck (BN1), where lBN1 ¼ min li , for i = 1 to n. Let k = 1. – If there is more than one BN1 (with the same l), mark the one with the smallest sequence number. 2. Identify the next bottleneck (i.e., BNk+1) in front of the previous one (i.e., BNk), where lBNkþ1 ¼ min li , for i = 1 to BNk  1. – If there is more than one BNk+1 (with the same l), mark the one with the smallest sequence number. If BNk+1 = 1, stop. Otherwise, let k = k + 1, go to 2. Determining the parameters Let k = n, fi = 1 for i = 1 to k. If station k is marked as a bottleneck, fi yk fi for i = 1 to k  1. Otherwise, fk xk fk. Stop if k = 2. 6. Otherwise, let k = k  1, go to 5.

3. Stage II: 4. 5.

Proposition 1. The system mean queue times from Procedure 1 and 2 are identical. Proof. See Appendix. In Eq. (1), the value of fi always equals 1 for the system bottleneck, which implies that a unit weight is always given to the system bottleneck. However, there will be a (non-unit) weight on all other stations’ ASIA system mean queue times. When both xk and yk are smaller than 1, the contribution factor will be smaller than 1 and behaves like a discount factor. In this situation, reducing bottleneck service time SCV brings greater improvement on queue time than reducing non-bottleneck service time SCV. 2.2. Models for manufacturing systems with single-server bottlenecks A practical factory is more complex than single-server queues in series. It may have long process flows with reentry and rework. Each workstation may have multiple servers with different capabilities and each server may have complicated configuration and suffer different types of interruptions. Dispatching rules other than workconserving policies may be used. Modeling the behavior of a factory through analyzing all activities in detail can be a formidable task. Since the intrinsic ratio captures the dependence among stations (Wu, 2009), our approach is to abstract a simple form from Eq. (1) to describe the behavior of practical manufacturing systems. We hope this simple form can capture the underlying struc-

ture of system performance. We assume system stability can be achieved when the utilization of each station is smaller than one. Because fi is 1 for the system bottleneck, system mean cycle time can be expressed as

  n P QT i þ PT f ¼ fi ai 1qiq l1 þ PT f i i i¼1 i¼1     P 1 ¼ aBN 1qBN fi ai 1qiq l1 þ PT f ; q l þ

CT ¼

n P

BN

BN

i–BN

i

ð2Þ

i

where CT is mean cycle time, BN stands for bottleneck. Hence, qBN is the bottleneck utilization (or system utilization). PTf is mean total processing time (of a job in a factory). Total processing time can be approximated by the cycle time of a job in light traffic when the station has no cascading, and there is no batch and assembly. The first term of Eq. (2) is the ASIA system mean queue time of the bottleneck, and the second term is the gross mean queue time of the non-bottlenecks. Queue time of tandem queues is contributed by all stations and so does its variability. From Eq. (2) system variability is mainly composed of two parts: the bottleneck and non-bottlenecks, with the bottleneck part playing a key role in heavy traffic. The coefficient of the bottleneck (aBN) is the same as the variability in its ASIA system. The coefficients of the non-bottlenecks (fiai) are the weighted (by their contribution factors) variabilities in their ASIA systems. When all service times are exponential and the initial external arrival process is Poisson, all intrinsic ratios are one. Hence, fi is 1 and Eq. (2) reduces to the cycle time of a Jackson network. When all service times are deterministic, the intrinsic ratios are zero for the system bottleneck and the stations behind it. Hence, fi is zero for all non-bottlenecks and Eq. (2) reduces to the cycle time of a Friedman’s tandem queue (Friedman, 1965). In Eq. (2), since queue time is dominated by the bottleneck in heavy traffic, we will group those non-bottleneck queue times, and replace the (n  1) non-bottleneck performance curves by (n – 1) identical composite non-bottleneck performance curves. Furthermore, within the non-bottleneck performance curves, system queue time is dominated by the performance curve of the second bottleneck (i.e., the one with the highest utilization among all non-bottlenecks) in heavy traffic. Therefore, a reasonable choice is to pick up the second bottleneck performance curve to represent the composite non-bottleneck performance curves. Based on the above observations, Eq. (2) can be simplified as

    P 1 CT ¼ aBN 1qBN fi ai 1qiq l1 þ PT f qBN lBN þ i i i–BN     0 qBN k=k ffi k1 1q l1 þ ðn  1Þk2 1k=k3 3 k13 þ PT f BN BN     1 k 1 ¼ k1 1qBN q l þ k2 k3 k k3 þ PT f ; BN

ð3Þ

BN

where k is the initial external arrival rate at the first station. In Eq. (3), the first term on the right-hand side corresponds to the bottleneck mean queue time where k1 is the bottleneck variability. The second term corresponds to the mean queue time of the non-bottleneck stations where k2 approximates the variability of the composite station (representing the (n  1) non-bottleneck stations), and k3 represents the capacity of the composite station. Note that k1 is the same as aBN, the bottleneck variability in the ASIA system if the system is a tandem queue. Since there are three parameters in Eq. (3),

Author's personal copy

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

we call it the 3-parameter model. When there is reentry or rework, capacity can be computed as follows,

1=li ¼

l P

wj  ST j ;

ð4Þ

j¼1

where l is the total reentry or rework frequency at station i, wj is the rework percentage (or rework rate) when STj is the length of a rework, and wj is 1 when STj is the length of a reentry. In addition to its simplicity, Eq. (3) considerably reduces the burden of data collection in Eq. (2). In order to apply Eq. (2), we assume the service time SCV is available. However, in a setting of complex machine configurations, robot scheduling, interruptions, resource contention, reentry, and product mix, finding out the variance of service time, where service time is the reciprocal of capacity rather than processing time (Wu et al., 2011), can be a formidable task. Newell (1979) stated, ‘‘In fact, in most applications, one is lucky if one has a good estimate of the service rates (to within 5% say); the variance rates are often known only to within a factor of 2, seldom to within an accuracy of 20%.’’ Although SCV of service time is well defined in theory, it may not be accessible in practice. A common way to estimate service time SCV is through analyzing historical data. Although it is difficult to have a reliable estimate of service time SCV, it is relatively easy if we just need to estimate the mean queue time. In practice, observable mean queue time (first moment estimator) is much more accessible than the intangible service times SCV (second moment estimator). In order to make queueing models more accessible in practice, it is important to have an approach which does not rely on the service time SCV explicitly. The historical mean queue time is a good alternative. Indeed, except when we construct a brand new factory, the historical performance is usually accessible in practical applications of queueing theory. Hence, it is practical to develop an approximate model which can predict future performance based on few reliable historical data estimates. Compared with the approximate model by Wu (2005), Eq. (3) gauges the variability of a manufacturing system from the viewpoint of the bottleneck while adding a correction term to consider the impact from non-bottlenecks. 2.3. Model validation through regression analysis In this section, we validate the cycle time approximate model in Eq. (3) by the simulation model of a real manufacturing system from industry. Planning and managing major defense acquisition programs (DAP) requires balancing and synchronizing design, and manufacturing across a network of distributed activities performed by independent commercial entities. To achieve this goal, it is important to have the capability to describe the performance curve of each independent entity accurately. Hence, a simulation model is constructed by one of the manufacturers in a DAP using ARENAÒ. The model (illustrated in Fig. 2) describes the behavior of a manufacturing facility for a specific product. The output of this facility supplies critical parts to downstream assembly lines in the supply chain. There are 14 workstations (R1, R3, R4, R5, R6, R8, R11, R12, R13, R14, R15, R16, R17, R19) arranged in 11 main process groups (i.e., 9 process groups use 1 workstation, 1 group needs 2 workstations, and 1 group uses 3 workstations). Four workstations have multiple servers (R5 = 4, R11 = 2, R15 = 2, R16 = 6), and the remaining workstations have only a single server. In addition to workstations, each process group requires operators from one of the three operator types (R9, R10 and R18). While the operators have their own shift schedules (i.e., 8 h a day), machines work 24 h a day as long as a job has been loaded by operators.

331

In the simulation model, service times follow triangular distributions. Dispatching rules at critical workstations follow the shortest remaining processing time instead of First-Come-First-Serve (FCFS). Reentry and rework existing in the system are shown in Fig. 2. 12 different raw parts arrive every 12 days with random batch sizes following a Poisson distribution. System utilization is determined by the mean batch size rather than the arrival intervals. Before the process can be started, raw parts (Z1 and g2–g7) need to be assembled in front of the first process step. Subsequently, an incomplete job has to be assembled again with some other raw parts (g8–g12) in the middle of the process flow. The cycle time of a job is the duration between its process start and departure. Due to the many factors such as shift schedule, operator availability, batch arrivals and assembly lines, finding a reliable estimate of inter-arrival time and service time SCVs can be very difficult, not to mention computing the cycle time exactly. To get a reliable approximation of cycle time, we resort to simulation. Through experimentation, the system bottleneck has been identified to be R8, which is composed of a single server, and system capacity estimated as 13.7 jobs per 12 days (indeed, even finding the bottleneck and its true capacity is not trivial). In total, system performance at 14 different utilizations have been simulated. The performance curve with 99% confidence intervals is shown in Fig. 3. The cycle time at each utilization is the average of 30 batches from one long simulation run. In each simulation run, the first 1000–10,000 data points are discarded for warm-up based on the utilizations (e.g. 1000 in light traffic and 10,000 in heavy traffic), and each batch consists of 1500–500,000 data points. In the following paragraphs, we will fit the performance curve in Fig. 3 by three different methods. The first one is suggested by Sattler (1996), Collins et al. (1997), Wu (2005), who attempted to describe the performance curve based on Kingman’s approximation,

 CT ¼ k1



qBN 1 þ PT f : 1  qBN lBN

ð5Þ

Since there is only a single parameter in Eq. (5), we call it the singleparameter model. Motivated by Kingman’s approximation, Yang et al. (2007) proposed the following model to describe the cycle time of manufacturing systems:

Pt CTðx; c; pÞ ¼

k k¼0 ck x p

ð1  xÞ

;

ð6Þ

where x is system throughput rate. p, t and the vector c = (c0, c1, . . . , ct) are unknown parameters in the model. We call Eq. (6) YAN’s model (from the initials of the three authors). Since t can increase without limit, we limit the value of t to 2 to have a fair comparison with Eq. (3). Hence, Eq. (6) becomes

CTðx; c; pÞ ¼

c0 þ c1 x þ c2 x2 : ð1  xÞp

ð7Þ

There are three terms in Eq. (7), which is the same as Eq. (3) and four parameters (i.e., p, c0, c1, c2) need to be determined, which is one more than Eq. (3). In order to compare the performance of the three models, the simulated cycle time in Fig. 3 is used to fit the parameters in Eqs. (5), (7) and (3) by the statistical software, R. The value of k1 in Eq. (5) (i.e., G/G/1 curve in Fig. 4) is 0.642, and the largest error is 11.77% at 91% utilization. The values of c0, c1, c2 and p in Eq. (7) (i.e., YAN’s model curve in Fig. 4) are 639.452, 0, 0 and 0.200, respectively, and the largest error is 11.87% at 4% utilization. The values of k1, k2 and k3 in Eq. (3) (i.e., 3-Parameter model curve in Fig. 4) are 0.356, 422.951 and 14.636, and the largest error

Author's personal copy

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

Fi na lP

ro du c

t

332

Z1

B=

g2

P= 1, 1 B=1, P=0.05

g3 B~Pois(9), λ_C=1/12days

R10, R17 Δ(1,2,3)Hrs

B~Pois(9), λ_C=1/12days

A S S E M B L Y

B=1, P=0.95 D=Δ(4,5,6) days

1,

g5

R9, R12 Tria(3,4,5)M Δ(12,16,20)M Tria(3,4,5)M

B=1, P=0.95

B=

R10, R16(6) Δ(7.3,9.6,12.9)M Δ(7.3,9.6,12.9)Hrs Δ(2,3,4)+Δ(1,2,3)H Δ(1,2,3)H Δ(2,4,6)+Δ(1,1.5,2)H Δ(1,1.5,2)H 2*Δ(2,3,4)H

g4

P=

B=1 B=1 P=0.98 P=0.98

0. 95

B=1 P=0.05

g6 B=1, P=1

B=1, P=0.3//0.005

B=1, P=0.75

g7 B~Pois(9), λ_C=1/12days

g8

B=1, P=0.05

B=1, P=1 B=1, P=1

B=1, P=1 B=1, P=?? ?? B=1, P=1

99 B=1, P=

.5 1

R9, R11(2) Δ(1.5,2,3)M Δ(0.15,0.2,0.3)M Δ(0.15,0.2,0.3)M Δ(0.15,0.2,0.3)M Δ(1.5,2,3)M Δ(1.5,2,3)M Δ(0.15,0.2,0.3)M Δ(0.75,1,1.25)H

B=1,P=0.05 Loop-3

B=1, P=0.7

R9, R6 Δ(30,45,60)M Δ(30,45,60)M

B=1, P=1 B=1, P=1

R9, R19 Δ(1.5,2,3)H Δ(1.5,2,3)H

B=1, P=0.95 B=1, P=0.17

B=1, P=0.25

B=1, P=1

ASSEMBLY

2

P=0.0

P=

B=1, P=1

B=1,

1,

B=1, P=0.25

B=1, P=1

B=

R10, R15(2) Δ(0.1,0.2,0.4)H Δ(0.1,0.2,0.4)H Δ(0.15,0.2,0.4)H Δ(0.15,0.2,0.4)H Δ(0.5,0.75,1)H 2*Δ(0.5,0.75,1)H 2*Δ(0.5,0.75,1)H 2*Δ(0.5,0.75,1)H 2*Δ(0.5,0.75,1)H Δ(0.5,0.75,1)H

1, B=

1 P=

B=2

=0 ,P

B=2, P=0.98

B=1, P=1

1 B=

g10

.65

g9

B=1, P=1, if not 1st Time

R9, R1 Δ(5,6,7)M

R5(4)

g12 B~Pois(9), λ_C=1/12days

B=5, P=0.98

R9, R13, R14 Δ(35,40,45)M Δ(35,40,45)M

g11

B={1,2}

R9, R8 Δ(30,60,90)M Δ(60,90,120)M Δ(30,60,90)M Δ(30,60,90)M

B=1, P=0.02

B=1, P=0.35

(R9/R18) (R3/R4) B=1, P=??

4H

B=1, P=0.65

Notes: -- Priority of jobs is ignored for initial analysis -- whenreevr proceesing time is written as 2*( ) => that the same processing time distribution accurs twice In a row and not that the random no. is multiplied by 2!

B=1, P=0.02

Δ(1.5,2,3)M

Δ(2,3,5)M

Fig. 2. Process flows of a manufacturing facility in DAP.

2000 1800

Cycle Time (Hours)

1600 1400 1200 1000 800 600 400 200 0 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization

Fig. 5 shows the fitting errors of the above three models. The three-parameter model gives the best fit among the three. Furthermore, the three-parameter model performs very well for utilizations above 60%, which represents a realistic operating range in practice. YAN’s model gives large errors in light traffic since Eq. (6) is inspired by queue time instead of cycle time models. The ignored processing time is responsible for those errors. The simulation results exhibit the prediction capability of Eq. (3) in practical manufacturing systems. The question becomes how to approximate k1, k2 and k3 accurately without first knowing the performance curve since simulation models are not always available and we may not have so many historical data points spread out in a wide range of utilizations in practice. In the next section, we develop a procedure to approximate k1, k2 and k3 based on the insight introduced previously. For the purpose of illustration, we use a simpler model, where the capacity and bottlenecks are easier to identify.

Fig. 3. Performance curve of a manufacturing facility in DAP.

2.4. Implementation in manufacturing systems with single-server bottlenecks is 1.38% at 29% utilization. Note that the value of k3, 14.636, which represents the non-bottleneck capacity, is greater than the bottleneck capacity 13.7 as expected.

In practice, a complete performance curve is not available in general. Hence, we have to approximate system mean queue time

Author's personal copy

333

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

k3 ¼ min li : for i – BN

ð8Þ

After k1 and k3 are determined, k2 can be determined by regressing the historical factory cycle time against Eq. (3) with the values of k1 and k3 specified. In this approach, one point estimator for each of the historical factory and bottleneck queue time suffices. If multiple data points are available, regression analysis can be used to get a better estimate of k2. We call this method the k2 model. One nice feature of the k2 model is that we can approximate system mean queue times without knowing the service time SCV, which is difficult to obtain in general. However, the challenge comes from estimating of the bottleneck mean queue time in the ASIA system. Fortunately, this issue can be relieved by the heavytraffic bottleneck phenomenon (Iglehart and Whitt, 1970): the queue time distribution at the bottleneck is asymptotically the same as if the immediate arrival process were replaced by the external initial arrival process to the first queue. Since this condition is the same as our requirement for the ASIA system, we could approximate k1 by using the historical bottleneck queue time in heavy traffic. Based on this approach, we can predict factory performance based on the historical queue times without knowing service time variability and the entire performance curve. 2.4.1. Model validation through historical data approach In this section, the k2 model is validated by a flow shop simulator developed by DSN Innovations (or DSN). It can be quickly populated with data from a small manufacturer, which is called Doyle Center Model, and used to evaluate system performance. Although the kernel of the tool is ARENAÒ, it uses Excel as its input interface so working knowledge of ARENAÒ is not required. DSN provided a case based on a manufacturing system consisting of five workstations as shown in Fig. 6. While station 1 and 5 are visited only once, station 2, 3 and 4 are visited multiple times. Stations 2, 3 and 4 can execute multiple job functions: station 2 can do two different recipes, station 3 can do three and station 4 seven. There are total 14 process steps. There is a constant 20 minutes delay between the first and the second steps. Some steps may need to be reworked if the finished jobs are out of spec. All stations are composed of one single machine except for station 5, which con-

1900 G/G/1

1700

Cycle Time (Hours)

Simulation

1500

3-Parameter Model YAN's Model

1300 1100 900 700 500 0%

10%

20%

30%

40%

50%

60%

70%

Utilization Fig. 4. Fitting results of the three models.

80%

90%

100%

10% G/G/1 Errors 3-Parameter Model Errors

5%

YAN's Model Errors

Error %

by using the insight on Eq. (3), along with specific information about the first and second bottlenecks in the system. If the bottleneck queue time in its ASIA system is known, we can estimate k1 directly based on the first term in Eq. (3). If we let the second bottleneck represent all non-bottleneck stations, then k3 can be approximated by the capacity of the second bottleneck, i.e.,

0% 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

-5%

-10%

-15%

Utilization Fig. 5. Fitting errors of the three models.

tains 24 machines in parallel. The initial arrival process is Poisson. All dispatching rules are FCFS. The service time distribution is triangular and the parameters at each process step are shown in Table 1. Based on Table 1 and Fig. 6, the capacities of station 1–5 are 62.609, 51.429, 49.655, 49.021 and 54.857 jobs/day and service times are 23, 28, 29, 29.375 and 630 minutes, respectively. The forth station is the bottleneck. Because of the complex reentry and rework as well as the small service time variability, it is difficult to obtain reliable cycle time approximations for this small scale manufacturing system by conventional queuing theory approaches. Thus, a viable approach would be simulation. The simulated cycle times and their halfwidth 95% confidence interval are shown in Table 2. Each observation in Table 2 is the average of 100 replications. Each replication is the collection of the output data in 10,000 days after a warm-up period of 100 years. The total system mean processing time was estimated to be 761.05 minutes by simulating at a very low shop throughput rate. In the following illustrations, we validate the historical data approach by using the simulation results as the ‘‘historical’’ data. Because k1 can be approximated by the bottleneck variability in heavy traffic, we use the bottleneck mean queue time at 93.8% utilization to represent its performance in heavy traffic. Therefore, the bottleneck mean queue time is 23.046 minutes at 93.8% utilization. Through the first term in Eq. (3) (and assuming the other two terms are zero), k1 is 0.052. Since k3 is the capacity of the second bottleneck, k3 is 49.655 jobs/day (i.e., capacity of the third station). To estimate k2, we use the system mean queue time at 81.6% utilization to represent the historical performance. Hence, historical mean cycle time is 922.65 minutes at 81.6% utilization. Based on the values of k1, k3 and the system mean cycle time at 81.6% utilization, k2 is 1857.412. Through Eq. (3), mean cycle times at other utilizations can be approximated as shown in Table 2 (specified under the k2 model). The simulated and approximate cycle times are shown in Fig. 7. By using only few historical data points, we get a reliable approximation of system mean cycle time based on the insight on k1, k2, and k3 values. The biggest error among all utilizations is 0.97% as shown in Table 2. Since we compute k2 based on the system mean cycle time at 81.6% utilization, the error is zero. However, the result is attained by assuming k2 is a constant across all utilizations. It would be interesting to verify this assumption by computing the values of k2 across all utilizations. The results are shown in Fig. 8. The k2 values all fall between 1790 and 2525, which explains the good approximate results in Fig. 7. The almost-linear trend suggests that the approximation can be further improved by extrapolating the k2 values if more than two queue times are known at different utilizations.

Author's personal copy

334

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

Fig. 6. Process flow of the DNS model.

Table 1 Process times of each process step. Process number

1 2 3 4 5 6 7 8 9 10 11 12 13 14

Process description

Bond overlay and O-ring Integrate harness at Station 1 Integrate circuit boards at Station 2 Integrate harness at Station 3 Perform preliminary Func Test Install top cover Perform vacuum test Perform acceptance test Prepare for burn-in test Burn-in (environmental) test Perform final accept test Touch-up power supply Serialize power supply Package power supply

One possible explanation for the declining curve is that the true value of k3 is overestimated by the second bottleneck capacity. For a given arrival rate (k), the second term in the model is underestimated, and the effect is exaggerated in heavy traffic. It would be interesting to compare the k2 model with the regression results of the G/G/1 and YAN’s models. The simulated mean cycle time in Table 2 is used to fit the parameters in Eqs. (5) and (7). The results are shown in Fig. 9. The value of k1 in Eq. (5) (i.e., G/G/1 model) is 1.117, and the largest error is 1.81% at 73.4% utilization. The values of c0, c1, c2 and p in Eq. (7) (i.e., YAN’s model) are 718.164, 0, 0 and 0.178, respectively, and the largest error is 5.26% at 81.6% utilization.

Process resource

Station 2 Station 1 Station 2 Station 3 Station 3 Station 3 Station 4 Station 4 Station 4 Environ. chamber Station 4 Station 4 Station 4 Station 4

Processing times (minutes) Min

Mode

Max

6 22 20 18.5 2 4.5 6 2 4 600 2 1 7 3

7 23 21 19.5 3 5.5 7 3 5 600 3 2 7 3

8 24 22 20.5 4 6.5 8 4 6 600 4 3 7 3

The errors are summarized in Fig. 10. Although the k2 model uses fewer historical data points, it outperforms the fitted results from the G/G/1 and YAN’s models, which need much more historical information. Hence, the k2 model is a better model to predict the mean cycle time of manufacturing systems.

3. Performance evaluation for systems with multiple-server bottlenecks In practical manufacturing systems, a station may consist of multiple servers. An extension of Eq. (3) to a multiple server station

Author's personal copy

335

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339 Table 2 Cycle times from simulations and the historical data approach. Input rate (per day)

Utilization (%)

4 8 12 16 20 24 28 32 36 40 44 46

8.2 16.3 24.5 32.6 40.8 49.0 57.1 65.3 73.4 81.6 89.8 93.8

Simulation

k2 model

CT (minutes)

95% Cl

CT (minutes)

Error (%)

765.6 771.0 777.0 784.2 792.9 803.7 817.6 837.2 867.5 922.7 1054.9 1245.1

0.15 0.10 0.09 0.08 0.08 0.08 0.08 0.13 0.18 0.38 1.03 3.18

764.5 768.5 773.5 779.6 787.3 797.5 811.4 831.7 863.8 922.7 1065.2 1254.6

0.15 0.32 0.45 0.59 0.70 0.77 0.76 0.66 0.43 0.00 0.97 0.76

1300

Simulation 1200

k2 Model

Cycle Time

1100

1000

900

800

700 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization Fig. 7. Performance curve comparison for the DNS model.

is necessary. When a workstation is composed of multiple servers, Kingman’s approximation can be modified by the approximate model proposed by Sakasegawa (1977) as follows,

roughly the same among those non-bottlenecks with higher utilizations (e.g. the second and the third bottleneck stations), choosing the second bottleneck should suffice.

 2  pffiffiffiffiffiffiffiffiffiffiffiffi ca þ c2s q 2ðmþ1Þ1 1 QT ffi ; 2 ð1  qÞ ml

3.1. Model validation through historical data approach

ð9Þ

where m is the number of servers of the workstation and q = k/(ml). Whitt (1993) further improved this model by adding correction factors. Based on Eq. (9), Eq. (3) can be modified as follows,

0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2ðmBN þ1Þ1 BN A

q

CT ffi k1 @

ð1  qBN Þ

 1 þ k2 mBN lBN

k m 2 k3



ffiffiffiffiffiffiffiffiffiffiffiffiffi p2ðm 2 þ1Þ1

1

k m 2 k3



1 þ PT f m2 k3 ð10Þ

where mBN is the number of servers at the bottleneck and m2 is the number of servers at the non-bottleneck. Other parameters are the same as Eq. (3). When k3 is the second bottleneck capacity, m2 is the number of servers at the second bottleneck. A potential issue of Eq. (10) is that the non-bottleneck mean cycle time may not be dominated by the second bottleneck when its server numbers are considerably larger than the third (or forth, and so on) bottleneck workstations. In this situation, the selection of k3 may have to be adjusted. However, when server numbers are

In order to verify Eq. (10), the original DSN model has been modified. The servers of all workstations (as well as the input rate) become 5 times more than the original model. Therefore, all stations now are composed of 5 servers in parallel except for station 5, which contains 120 servers. The processing times and process flows are the same as the original case. From simulation results, the bottleneck mean queue time is 8.58 minutes at 94.8% utilization. Hence, the k1 value is 0.087. The k3 value can be approximated by the second bottleneck capacity 248.276 (i.e., capacity of the third workstation). The total system mean processing time is estimated to be 761.05 minutes by simulating at a very low shop throughput rate. We use the mean queue time at 81.6% utilization to represent the historical mean cycle time. Eq. (10) can be solved for k2 as 3899.9. The approximated mean cycle times and corresponding errors at other utilizations are summarized in Table 3 and Fig. 11. The largest error is 3.94% at 94.8% utilization. In the above approach, we assume k2 is a constant. It is important to examine how k2 behaves. When k1, k3 and the queue times

Author's personal copy

336

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339 5000 4500 4000

k2 Value

3500 3000 2500 2000 1500 1000 500 0 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization Fig. 8. k2 values at different utilizations.

1300

except in light traffic, where the queue time is negligible. This trend suggests that we may improve the previous results by using a two point approximation method. To illustrate, suppose the historical cycle times are known at 81.6% and 89.8% utilizations. Thus, the k2 values can be computed at those two utilizations. They are 3901.6 and 4592.1, respectively. The other k2 values can be extrapolated accordingly as shown in Table 4 (in the column of ‘‘k2_Exprapolation’’). The new performance curve is shown in Fig. 13, where the largest error is only 1.12% at 94.8% utilization. Since k2 should not be negative, we can further improve the results by replacing the negative extrapolated values with zeros in light traffic.

Simulation

1200

k2 Model

G/G/1 Model

Cycle Time (Min)

1100

YAN's Model

1000

900

800

700 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization

From Eqs. (3) and (10), we observe that system mean queue time can be decomposed into two parts: the bottleneck and the non-bottlenecks. System variability also corresponds to those two parts accordingly. The parameter k1 represents the variability of the bottleneck, and k2 represents the variability of the non-bottlenecks. Based on the derivation of Eq. (3), system variability is lower with less service time variability, initial arrival process variability, or the number of non-bottleneck stations. Manufacturing system performance holds the following properties:

Fig. 9. Fitting results of the three models.

10% 8%

G/G/1 Model

6%

YAN's Model

Error %

4%

4. Theory of system performance

k2 Model

2% 0% 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

-2%

Property 1. Relative Measure of System Performance. The performance of a manufacturing system is improved if either the value of k 1 or k2 becomes smaller across the operational rage of system utilizations.

-4%

Property 2. Absolute Measure of System Performance. Among a group of manufacturing systems, a system performs better if its k1 and k2 values are smaller than the values of other systems across the operational rage of system utilizations.

-6% -8% -10%

Utilization Fig. 10. Fitting errors of the three models.

at different utilizations are known, the corresponding k2 values can be computed as shown in Fig. 12 and Table 4 (in the column of ‘‘k2’’). The displayed values of k2 exhibit an almost-linear trend

Note that shorter cycle time does not necessarily imply better performance since total processing time, non-bottleneck service times, and the number of servers at a workstation also affect the values of k1 and k2. If system A has lower k1 and higher k2 than system B at a specific utilization, the bottleneck of system A performs better, but its non-bottleneck performs worse at that specific utilization.

Author's personal copy

337

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339 Table 3 Errors of the k2 model based on single historical data point. Input rate (per day)

Utilization (%)

20 40 60 80 100 120 140 160 180 200 220 232

8.2 16.3 24.5 32.6 40.8 49.0 57.1 65.3 73.4 81.6 89.8 94.8

Simulation

k2 model

CT (minutes)

95% Cl

CT (minutes)

Error (%)

761.1 761.2 761.2 761.6 762.2 763.3 765.3 768.4 774.0 785.1 818.2 892.8

0.07 0.04 0.03 0.03 0.03 0.03 0.03 0.03 0.02 0.04 0.12 0.43

761.1 761.2 761.6 762.1 763.0 764.4 766.4 769.6 774.9 785.1 810.2 857.6

0.00 0.01 0.04 0.07 0.11 0.14 0.15 0.15 0.12 0.00 0.98 3.94

950

Simulation 900

Cycle Time

k2 Model

850

800

750

700 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization Fig. 11. Performance curves based on single historical data point.

6000.0

5000.0

k2 Value

4000.0

3000.0

2000.0

1000.0

0.0 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization Fig. 12. k2 values at different utilizations.

5. Conclusion Through the nice property of intrinsic ratios, we developed analytical models to quantify the performance of manufacturing systems. From the simulation experiment results, we find the intrinsic ratios discovered in tandem queues can also be applied

to determine the performance of manufacturing systems. Comparing the proposed model with prior approaches, we find that Eq. (3) not only performs better but also permits insight into those parameters. Quantifying system performance exactly for practical manufacturing systems is difficult. On the other hand, approximating

Author's personal copy

338

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339

Table 4 Cycle times from simulations and the k2 model with two point approximation. Input rate (per day)

Utilization (%)

20 40 60 80 100 120 140 160 180 200 220 232

8.2 16.3 24.5 32.6 40.8 49.0 57.1 65.3 73.4 81.6 89.8 94.8

Simulation CT (minutes)

95% Cl

761.1 761.2 761.2 761.6 762.2 763.3 765.3 768.4 774.0 785.1 818.2 892.8

0.07 0.04 0.03 0.03 0.03 0.03 0.03 0.03 0.02 0.04 0.12 0.43

k2

k2_Extrapolation

k2 model CT (minutes)

Error (%)

2325.8 2202.2 1269.8 1765.9 2157.9 2615.4 3011.7 3330.0 3610.7 3901.6 4592.1 5459.6

2313.1 1622.6 932.0 241.5 449.0 1139.5 1830.0 2520.6 3211.1 3901.6 4592.1 5015.3

761.0 761.0 761.0 761.0 761.4 762.1 763.7 766.8 772.7 785.1 818.2 882.7

0.00 0.02 0.04 0.07 0.11 0.16 0.20 0.22 0.17 0.00 0.00 1.12

950

Simulation 900

Cycle Time

k2 Model

850

800

750

700 0%

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Utilization Fig. 13. Performance curves based on two point k2 method.

system performance by requiring too much information is not feasible. In this paper, we propose reliable methods to approximate system performance by only requiring few historical data without resorting to the second moment of service time. This important result constitutes a bridge between theory and practice. The almost-linear property of k2 value is used to improve the model accuracy as explained in Table 4. An intuitive explanation on this result is that the non-linear portion of mean queue time has been taken care by the utilization term, while the remaining portion of mean queue time dominated by variability behaves almost linearly. From our experiments, the almost-linear property does exist in both Figs. 8 and 12. However, similar to the linearity property of intrinsic ratios, this linearity property in our model should rely on the accuracy of Kingman’s approximation. In the examined cases with Poisson arrivals, Kingman’s approximation reduces to an M/G/1 queue and gives exact results. When the arrival process is not Poisson, although the k2 value would still behave more linearly than queue time (since a major non-linear portion from utilization term has been normalized), the pattern may not be as regular as the Poisson cases. A more detailed examination on this property in this situation is left for future research. Although our models perform well in the examined cases, the practical situation can be much more complicated. In real production lines, the values of k1, k2 and k3 can be affected by practical issues such as reentry, interruptions and dispatching rules, where the dispatching rules can be more complex than work-conserving or decentralized control policies. Furthermore, there can be multi-

ple products in a production system simultaneously. They may have different cycle times but share the same resources. Those situations have not been considered in our models and will be left for future research. Acknowledgement The authors are grateful to Dr. Hayriye Ayhan and Bert Zwart for their kindly help, and to Nikhil Kumar and Abhinav Dalal for their support in simulation experiments. The first author would like to thank Yawen Cheng for her generous help in improving the readability of this paper. Appendix A Proof of Proposition 1 (i) When n = 1 (i.e., the system has only one single server), based on the 4th and 5th step in Procedure 2, f1 is 1, which correctly describes this single server system queueing time. (ii) When n = 2 (i.e., adding the 2nd server to the single server system), the 2nd server can be either a bottleneck or a non-bottleneck.   q If it is a bottleneck, its queueing time is QT 2 ¼ a2 12q l1  2 2 P2 ð1  y2 ÞQT 1 . Total system queue time is i¼1 QT i ¼QT 1 þQT 2 ¼     P a2 1q2q l1 þy2 QT 1 ¼ an 1qnq l1 þyn n1 i¼1 QT i . 2

2

n

n

Author's personal copy

K. Wu, L. McGinnis / European Journal of Operational Research 221 (2012) 328–339



q



If it is a non-bottleneck, its queueing time is QT 2 ¼ x2 a2 12q 2 P2 i¼1 QT i ¼QT 1 þQT 2 ¼QT 1 þx2 a2 l2 . Total system queue time is     Pn1 q2 qn 1 1 i¼1 QT i þxn an 1q l . 1q l ¼ 1

2

n

2

n

(iii) When n = k, the kth server can be either a bottleneck or a non-bottleneck. If it is a bottleneck, based on Procedure 1, its queueing time is   P q QT k ¼ ak 1kq l1  ð1  yk Þ k1 i¼1 QT i ;. Total system queueing time k k   Pk Pk1 Pk1 qk 1 is i¼1 QT i ¼ i¼1 QT i þ QT k ¼ ak 1q l þ yk i¼1 QT i . A weight k

k

yk is given to all the servers in front of the kth server as described in Procedure 2.   If it is a non-bottleneck, its queue time is QT k ¼ xk ak 1qkq l1 . k k Pk P P Total system queue is QT i ¼ k1 QT i þ QT k ¼ k1 QT i þ i¼1 i¼1 i¼1   xk ak 1qkq l1 . A weight xk is given to the kth server as described k

k

in Procedure 2. (iv) When n = k + 1, the (k + 1)th server can be either a bottleneck or a non-bottleneck. If it is a bottleneck, based on Procedure 1, its queue time is   P qkþ1 1  ð1  ykþ1 Þ ki¼1 QT i . Total system queue 1qkþ1 lkþ1   Pk Pkþ1 qkþ1 1 time is i¼1 QT i ¼ i¼1 QT i þ QT kþ1 ¼ akþ1 1qkþ1 lkþ1 þ ykþ1 Pk i¼1 QT i . A weight yk+1 is given to all the servers in front of the kth server as described in Procedure 2. If it is a non-bottleneck, its queue time is QT kþ1 ¼ xkþ1   Pkþ1 Pk 1 akþ1 1qkþ1 i¼1 QT i ¼ i¼1 QT i þ qkþ1 lkþ1 . Total system queue time is   Pk qkþ1 1 QT kþ1 ¼ i¼1 QT i þ xkþ1 akþ1 1q l . A weight xk+1 is given to QT kþ1 ¼ akþ1

kþ1

kþ1

the (k + 1)th server as described in Procedure 2. References Baskett, F., Chandy, K.M., Muntz, R.R., Palacios, F.G., 1975. Open, closed, and mixed networks of queues with different classes of customers. Journal of the ACM (JACM) 22, 248–260. Bitran, G.R., Tirupati, D., 1989. Tradeoff curves, targeting and balancing in manufacturing queueing networks. Operations Research 37, 547–564. Boebel, F., Ruelle, O., 1996. Cycle time reduction program at ACL. In: Proceedings of the Advanced Semiconductor Manufacturing Conference, Cambridge, MA, USA, pp. 165–168. Chandy, K.M., Herzog, U., Woo, L., 1975. Approximate analysis of general queueing networks. IBM Journal of Research and Development 19, 43–49. Collins, D.W., Williams, K., Hoppensteadt, F.C., 1997. Implementation of minimum inventory variability scheduling 1-step ahead policy in a large semiconductor manufacturing facility. In: Proceedings of the 6th ETFA International Conference, pp. 497–504. Cruz, F.R.B., Van Woensel, T., Smith, J.M., 2010. Buffer and throughput trade-offs in M/G/1/K queueing networks: a bi-criteria approach. International Journal of Production Economics 125, 224–234.

339

Fowler, J.W., Brown, S., Gold, H., Schoeming, A., 1997. Measurable improvements in cycle-time-constrained capacity. In: Proceedings of the IEEE International Symposium on Semiconductor Manufacturing, pp. A21–A24. Friedman, H.D., 1965. Reduction methods for tandem queuing systems. Operations Research 13, 121–131. Gurvich, I., Whitt, W., 2009. Queue-and-idleness-ratio controls in many-server service systems. Mathematics of Operations Research 34, 363–396. Harrison, J.M., 1998. Heavy traffic analysis of a system with parallel servers: asymptotic optimality of discrete-review policies. Annals of Applied Probability, 822–848. Harrison, J.M., Lopez, M.J., 1999. Heavy traffic resource pooling in parallel server systems. Queueing Systems 33, 339–368. Hulett, M., Damodaran, P., 2011. Analytical approximations to predict performance measures of markovian type manufacturing systems with job failures and parallel processing. European Journal of Operational Research. Iglehart, D.L., Whitt, W., 1970. Multiple channel queues in heavy traffic. II: sequences, networks, and batches. Advances in Applied Probability 2, 355–369. Jackson, J.R., 1957. Networks of waiting lines. Operations Research 5, 518–521. Kerbache, L., Smith, J.M., 1987. The generalized expansion method for open finite queueing networks. European Journal of Operational Research 32, 448– 461. Kerbache, L., Smith, J.M., 2000. Multi-objective routing within large scale facilities using open finite queueing networks. European Journal of Operational Research 121, 105–123. Kingman, J.F.C., 1965. The heavy traffic approximation in the theory of queues. In: Smith, W.L., Wilkinson, W.E. (Eds.), Proceedings of the Symposium on Congestion Theory. University of North Caroline Press, University of North Carolina at Chapel Hill, pp. 137–159. Li, J., Blumenfeld, D.E., Huang, N., Alden, J.M., 2009. Throughput analysis of production systems: recent advances and future topics. International Journal of Production Research 47, 3823–3851. Nazzal, D., Mollaghasemi, M., 2001. Critical tools identification and characteristics curves construction in a wafer fabrication facility. In: Proceedings of the Winter Simulation Conference, pp. 1194–1199. Newell, G., 1979. Approximate Behavior of Tandem Queues. Springer. Park, S., Mackulak, G.T., Fowler, J.W., 2001. An overall framework for generating simulation-based cycle time-throughput curves. In: Proceedings of the Winter Simulation Conference, Vol. 2, pp. 1178–1187. Rose, O., 2001. The shortest processing time first (SPTF) dispatch rule and some variants in semiconductor manufacturing. In: Proceedings of the Winter Simulation Conference, Vol. 2, pp. 1220–1224. Ruelle, O., 1997. Continuous flow manufacturing: the ultimate theory of constraints. In: Proceedings of the Advanced Semiconductor Manufacturing Conference, pp. 216–221. Sakasegawa, H., 1977. An approximation formula Lq = a  qb/(1-q). Annals of the Institute of Statistical Mathematics 29, 67–75. Sattler, L., 1996. Using queueing curve approximations in a fab to determine productivity improvements. In: Proceedings of the Advanced Semiconductor Manufacturing Conference. IEEE, pp. 140–145. Whitt, W., 1993. Approximations for the GI/G/m queue. Production and Operations Management 2, 114–161. Wu, K., 2005. An examination of variability and its basic properties for a factory. IEEE Transactions on Semiconductor Manufacturing 18, 214–221. Wu, K., 2009. New results in factory physics – insights from the underlying structures of manufacturing systems. Ph.D Dissertation, Georgia Institute of Technology. Wu, K., McGinnis, L., 2012. Interpolation approximations for queues in series. IIE transactions. http://dx.doi.org/10.1080/0740817X.2012.682699. Wu, K., McGinnis, L., Zwart, B., 2011. Queueing models for a single machine subject to multiple types of interruptions. IIE Transactions 43, 753–759. Yang, F., Ankenman, B., Nelson, B., 2007. Efficient generation of cycle timethroughput curves through simulation and metamodeling. Naval Research Logistics 54, 78–93.