Predicting performance measures for Markovian type of manufacturing systems with product failures

European Journal of Operational Research 184 (2008) 725–744 www.elsevier.com/locate/ejor

O.R. Applications

Predicting performance measures for Markovian type of manufacturing systems with product failures Salil Pradhan a, Purushothaman Damodaran

b,*

, Krishnaswami Srihari

c

a

b

FormFactor, Livermore, CA 94551, United States Department of Industrial and Systems Engineering, Florida International University, Miami, FL 33174, United States c Electronics Manufacturing Research and Services, Department of Systems Science and Industrial Engineering, State University of New York, Binghamton, NY 13902-6000, United States Received 4 April 2006; accepted 15 November 2006 Available online 10 January 2007

Abstract Optoelectronic products are typically assembled and tested in a flow shop environment with multiple processors at each stage. The first few stages are dedicated for assembly and the later stages are dedicated for calibration and testing. Whenever a product (or job) fails at a stage, it is routed back to one of the downstream stages or to the same stage (depending upon the nature of the failure). Consequently, the product could circulate several times between the current stage and the preceding stage(s) before moving to the next stage. Estimating the performance measures (such as WIP and flow time) of such manufacturing systems is not trivial. This paper presents analytical approximations to estimate the performance measures of a manufacturing system with multiple product classes, job circulations due to failures, and some resources being shared among different product classes. The analytical approximations were verified using simulation on several problem instances. The experimental study indicates that these approximations can be used by operations managers to estimate the performance measures of a manufacturing system with product failures.  2006 Elsevier B.V. All rights reserved. Keywords: Manufacturing; Queuing; Optoelectronics assembly

1. Introduction Electronics manufacturing service (EMS) providers are traditionally responsible for the assembly of surface mount components (SMCs) on printed circuit boards (PCBs). The EMS providers assemble/manufacture different types of electronic products ranging from miniature sized PCBs used in medical products to highly complex electronic assemblies for high-end networking products used in harsh environments. The industry has undergone several technological changes and challenges with respect to stabilizing the processes. Around the end of the 20th century, optical assembly was one another area embraced by the EMS providers. This *

Corresponding author. Tel.: +1 305 345 6912. E-mail address: damodarp@fiu.edu (P. Damodaran).

0377-2217/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2006.11.016

726

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

Fig. 1. Process flow in optical assembly operations.

expanded the scope of the manufacturing services that an EMS provider could offer to an original equipment manufacturer (OEM). Optoelectronics assembly and test operations are labor intensive and heavily dependent on operator skills. The assembly operations include splicing optical fibers, which are partially automated, and fiber management operations, which are labor intensive. These assembly operations are preferably performed in a single stage due to the constraints imposed by in-line insertion loss measurement processes (Pradhan et al., 2001). On completion of the assembly operations, the products are passed to the subsequent stages for calibration and testing. Some jobs may fail at the test stations for several reasons. The failed jobs are typically routed back to one of the preceding stations or to the same station (depending on the nature of the failure) for troubleshooting and rework. Consequently, jobs can circulate between the current stage and the preceding stage(s) several times before moving to the next stage. Fig. 1 presents the schematics of a typical optoelectronics assembly. The first stage is the splicing stage where the product is assembled. Stages subsequent to the first stage are the testing stages such as calibration (stage 2), qualification (stage 3), and burn-in (stage 4). Yield fallouts introduce some level of uncertainty in predicting the throughput. Some jobs might have to circulate multiple times between the stations until all the parameters meet the required specifications. For example, if a job fails at the qualification stage, it could be routed to the calibration stage for retesting, or to a splicing station for a rework. Irrespective of whether the job is routed to the first or second stage, it has to follow the regular sequence until it passes the qualification test. 2. Problem description Job circulation is not limited to optoelectronic manufacturing and test operations. Job circulations are commonly observed in semiconductor integrated circuit (IC) manufacturing (or re-entrant lines). In re-entrant lines, a job may visit one or more stages/machines on multiple occasions, depending upon the process flow. Although the jobs make multiple visits to a machine, the route is deterministic. Estimating the performance measure of a re-entrant line is not trivial. Consequently, estimating the performance measures of the manufacturing system under study is also not trivial as it considers backward flows due to probabilistic job failures. Since most low mix high volume production has moved to low cost regions, US manufacturers are left to deal with high variety and low volume manufacturing. Consequently, this paper considers multiple product classes. The equipment sets required for assembling and testing the optoelectronic products are customer specific.

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

727

Hence, most of them are provided by the OEM and are limited. It is imperative to utilize these resources effectively to meet the deadlines and to survive competition. It is critical for the planners and managers to plan for resources at various stations to improve the throughput in the presence of product failures. In summary, the operations department of such manufacturing systems would need a tool to predict performance measures such as WIP, flow time, throughput, and system bottlenecks in order to provide better service to the end customers. This research effort deals with a manufacturing system in which the jobs from multiple product classes could stochastically flow between different stations in the system due to product failures. The principal objective of the research is to predict the primary performance measures such as WIP, flow time, and throughput of such a manufacturing system using analytical methods, and verify the results using discrete event simulation. The analytical models can eventually be used in conjunction with a capacity planning model to determine the number of servers required at each stage in the presence of product failures and shared resources. The rest of the paper is organized as follows. A brief overview of the literature is presented in the next section. Then the assumptions made in this research, a brief discussion on Jackson networks and the analytical approximations for the problem under study are presented. Later, the performance measures for several problem instances are estimated using proposed analytical approximations and simulation. The error in performance measures between the two methods is analyzed and the possible sources of error between the analytical and simulation results are explained. 3. Literature review The flow of jobs between different stations in a semiconductor IC manufacturing industry can be represented deterministically and the jobs could make multiple visits to different stations. However, the number of visits would be fixed by the original process sequence. Such manufacturing systems are referred to as reentrant lines in the literature. Reentrant lines have been extensively studied in the past by Kumar (1993), Narahari and Khan (1995), Morrison and Kumar (1998), Yao (1994) and many other researchers with the objective of obtaining the performance measures of the manufacturing system. Most of these efforts have utilized queueing network theory as the principal tool for performance estimation. The problem under consideration deals with the stochastic flow of jobs between different stations due to job failures. Narahari and Khan (1995) studied the performance of reentrant lines using mean value analysis (MVA). This technique has been used for both open and closed queueing networks. In an open queueing network, jobs arrive externally at a specific rate into the system and the WIP in the system need not be constant. In a closed queueing network, a fixed number of jobs circulate in the system, as in a manufacturing system with conveyors (Kobza et al., 2002). Narahari and Khan (1996) studied the effect of inspection on reentrant lines using MVA. Park et al. (2000) used the MVA technique to analyze closed reentrant flow shops with single product class and batch machines. The analysis of such a system could render both exact and approximate solutions. Exact solutions are obtained when the job arrivals are Poisson and their processing times at different stations are exponentially distributed. Narahari and Khan (1998) presented an approximate but efficient method to compute the asymptotic loss of buffer priority scheduling policies in closed reentrant lines. The asymptotic loss refers to a situation where the population of the customers/jobs reaches infinity and the throughput reaches the maximum possible value i.e. the throughput of the bottleneck station. MVA was used to compute the asymptotic loss for closed queueing networks. Exact solutions were pioneered by Jackson (1957, 1963). Jackson (1963) proved that under the assumptions of Poisson arrivals and exponential service times, the performance of a class of networks could be predicted using the standard formulations for M/M/1 queues, and the results could be extended for M/M/c queues as well, where c is the number of servers at that stage. The primary result of Jackson network models is the equilibrium joint probability distribution of the component M/M/c queues. There exists a rich body of literature that applies Jackson networks to model and analyze a manufacturing system under various configurations such as flow shop, job shop, and flexible manufacturing system. Based on Jackson networks, many other queueing networks such as BCMP networks (Baskett et al., 1975) and Kelly networks (Kelly, 1979) have been developed with their principal application in the telecommunication industry.

728

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

Queuing networks for deterministic flow of jobs with multiple product classes have been extensively studied by Whitt (1983a), Segal and Whitt (1989), Bitran and Sarkar (1994) and Whitt (1992). Whitt (1983a) proposed the first version of the queueing network analyzer (QNA) with multi-server nodes, first-in-first-out (FIFO) service discipline, and infinite buffer. The primary idea in the paper was to compute the variability parameters for arrival and service rates and derive the performance measures based on these values. Segal and Whitt (1989) modified the original QNA to analyze manufacturing lines. The objective was to obtain the performance measures such as capacity, WIP, and production intervals to design or change a manufacturing line. The original QNA was later modified to account for machine breakdowns, batch service, lot sizes, and product testing with rework and partial yields. Bitran and Sarkar (1994) developed mathematical models to quantify the tradeoffs between capacity requirements and process improvements through variance reduction and throughput. Multiple product classes with single server queues were considered in the problem. Whitt (1992) extended the approximations developed by Bitran and Tirupati (1988) to develop approximations for the departure process of multi-class open queueing networks. In summary, the MVA technique has been extensively utilized for manufacturing systems with constant WIP and Markovian arrival and service distributions. To our knowledge, performance analysis of complex manufacturing systems (multiple classes, shared nodes, probabilistic routes) has not been dealt systematically in a way which would serve as a decision making tool for planners and managers in a high mix-low volume type of manufacturing environment. Consequently, this research deals with the manufacturing systems that represent the current state of the manufacturing industry in the US. However, the analytical formulations and the design of the mathematical model are generic enough to be applied to high volume-low mix systems as well. 4. Flow shops with job failures As mentioned earlier, optoelectronic assembly and test operations form a flow shop and job circulations occur due to failures. This section presents the analytical models developed for a multi-product system with shared nodes and job failures. The assumptions made are: • • • • • • • •

Poisson arrivals and exponential processing times. There exists an infinite buffer at each station in order to accommodate the WIP. Jobs follow a FIFO queuing discipline. There exists a set of dedicated servers at each station. In other words, the resources are not shared between different stations. Jobs could be routed back for rework, troubleshooting, and retest at the same station. Jobs are not preempted. Jobs which fail at any step during the process do not form a separate queue. The system is of an open type.

The probability of the job moving from one station to another is dependent upon the yield at that station. The yields are represented by a Markovian matrix (P) where the elements of the matrix represent the probabilities. The matrix is also called as the transition probability matrix (TPM). Several performance measures co-exist in a manufacturing system. In this research, throughput, WIP, and flow time were the performance measures of interest. Some of the secondary performance measures, such as scrap and yields, depend upon the performance of the product, which is a function of the manufacturing processes and product design. 4.1. Jackson networks Consider a system with two stations and one server at each station (see Fig. 2). Such a system can be denoted as M/M/1 ! ./M/1. Let k0 be the external arrival rate to the system and l1 and l2 be the processing rates at stations 1 and 2, respectively. The processing times at both the stations follow an exponential distribution. Jobs that successfully leave the first stage join the queue for service at the second stage. Jobs that fail at a stage are usually routed to the same stage or to one of the preceding stages. A node represents each stage and

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

729

Fig. 2. A two node-single server queuing system.

an arc represents the movement of the job between stages. Therefore, a manufacturing system can be represented in the form of a network consisting of a set of nodes and arcs. An arc (i, j) is associated with the transitional probabilities (qij), i.e. qij is the probability of a job to move from station i to station j. The state of the system at any time instance t is given by (n1, n2), where n1 and n2 are the number of jobs at stages 1 and 2, respectively. The general assumption made in many of the queuing models is that at steady state the service efficiency or traffic intensity (q = k/l) is less than 1. Therefore, assuming that the service efficiency at stage 1 (q1 = k/l1) and stage 2 (q2 = k/l2) is less than 1, then the steady state is reached for some large t and limt!1p(n1, n2, t) exists. It has been proven that (Medhi, 1991) pðn1 ; n2 Þ ¼ qn1 qn2 pð0; 0Þ and pð0; 0Þ ¼ ð1  q1 Þð1  q2 Þ; where p(0, 0) represent a system with zero jobs at stages 1 and 2 p(n1, n2, t) represent a system with n1 and n2 jobs at stages 1 and 2, respectively at any large t (t ! 1) Thus, pðn1 ; n2 Þ ¼ ½ð1  q1 Þðqn11 Þ½ð1  q2 Þðqn22 Þ:

ð1Þ

The two stages in the system are independent of each other because the arrival and service rates at each stage are exponentially distributed. Consequently, according to Burke’s theorem, the arrival rate at stage 2 is k (Burke, 1956). This result can be extended to a finite number of M/M/1 queues in series or in tandem as shown in Eq. (2): k Y ð1  qi Þqni i ; ð2Þ pðn1 ; n2 ; . . . ; nk Þ ¼ i¼1

where qi ¼ ki =li : The results can also be generalized for a series of M/M/c type of systems in a product form and can be extended to a more generic network in which jobs can enter the system at one or multiple nodal points, join the queue, receive service, and exit the system. Additional results for Markovian traffic flows between nodes have also been presented in the literature for the following: • Open and closed systems (Morrison and Kumar, 1998). • Systems with overflow discipline wherein the jobs which do not receive immediate service at the next node (j + 1) overflow the node and go to server (j + 2) (Chan, 2004). • Systems with loss discipline in which the jobs which do not receive service at node (j + 1) leave the system (Ku and Yi, 2002). • Systems with blocking discipline in which jobs that do not gain immediate access to server (j + 1) stay at server j until a space at (j + 1) is available (Koizumi, 2002). The generic model of network of Markovian queues is the Jackson network model. The input to the ith node consists of the input from the other nodes and/or the external input k0. The net arrival rate kj at node j is computed using the traffic equation (Jackson, 1954)

730

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

kj ¼ k0 q0j þ

X

qij ki

8j 2 I:

ð3Þ

i2I

Eq. (3) is also termed as flow balance equations or conservation equations (Medhi, 1991). For the problem under study, a job could visit a station more than once depending upon the yields. In such a case, some of the nodes could have more than one queue with different arrival rates and the jobs ultimately line up and form a single queue in front of each node. In order to compute the performance measures, the total number of visits made by a job to each node (vj) needs to be calculated. This can be calculated using an equation that is similar to the traffic equation for the net arrival rate (Jackson, 1954): X qij vi 8j 2 I: ð4Þ vj ¼ q0j þ i2I

Eqs. (3) and (4) for the net arrival rate and the number of visits respectively are valid for each internal node. By simultaneously solving these two sets of equations, the net arrival rate and the number of visits to each stage can be computed. 4.2. System with multiple product classes Eqs. (3) and (4) are applicable to a manufacturing system with a single product class with probabilistic flows. In a typical manufacturing system, it is common to produce multiple products simultaneously on different lines. Consequently, analytical results for a multiple product class scenario will be discussed. In such a system, multiple products could enter the system at different nodes (external) that are unique to each class. The jobs follow the predefined processing sequence until they exit the system from their respective exit nodes. Additional complexity arises if some of the internal nodes in the network are shared by two or more products. This situation often arises in the manufacturing facility when different product variants are being manufactured for the same customer. Such a system (for two product classes) is represented in Fig. 3. In this paper, the performance analysis of a system with exponential interarrival and service times (M/M/c) is described. The analysis of an M/M/c system is based on the modified version of the queueing network analyzer (QNA) developed by Whitt (1983a). QNA was originally developed by Bell Labs to calculate the performance measures for network of queues in telecommunications (with deterministic routes). The original version of the QNA has a much wider scope by including the G/G/c type of problems. The arrival and service rates are represented by a variability parameter, which is termed as the squared coefficient of variation (SCV). The SCV is given by the variance divided by the square of its mean. The performance measures of the system largely depend on these parameters. The parameters for the internal flows in the queue depend upon three basic network operations:

Fig. 3. Multiple product classes with one shared node.

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

731

• Superposition: Jobs from different product classes merge along different arcs to a node. • Splitting: A job from a particular class splits itself along different arcs to multiple nodes. • Departure: Jobs from different product classes exit the network. Splitting is not applicable for the problem under study since a job does not split itself into multiple jobs. The applicability of the splitting operation pertains to the communication industry where data packets are split into multiple packets at different nodes. Superposition and departure are the key operations under consideration since they are directly related to the queue behavior wherein the jobs merge from several nodes to a common node and later flow through the queue (departure process). The most important inputs to the QNA are the SCV of the arrival (Ca) and service (Cs) rates for each node that is visited by a product class. In order to compute the performance measures using the parametric decomposition approach, these parameters need to be aggregated at each node. The parametric decomposition method is an approach to approximating the steady state performance measures of open queuing networks with non-Poisson arrival process and non-exponential service time distributions (Whitt, 1983a,b). The idea is to analyze the individual queues separately after approximately characterizing the arrival processes. The arrival process is typically characterized by the rate and its variability. The arrival rates are obtained by solving the traffic equations (see Eq. (3)) and are exact. The approximating variability parameters (for example, SCV of interarrival times) are obtained by solving the traffic variability equations, which are approximations. The total network performance is computed by aggregating the individual queue performances. 4.2.1. Formulation of the analytical model for the multiple product class problem This section presents our analytical formulations to predict the performance measures of multi-class manufacturing systems with product failures and shared nodes. The formulations are based on the work done by Whitt (1983a) and are modified to incorporate the stochastic flow of jobs between different nodes in the network. The notation used is first presented. Notations k product class index, k 2 K Ik set of nodes visited by product class k total number of nodes visited by product class k nk number of servers at node i ci external arrival rate at node j for class k k0jk service rate at node j for class k ljk service time at node j for class k (1/ljk) sjk traffic intensity at node j qj number of visits made to node j by class k vjk transition probability to node j of class k from an external node (a virtual node from where the jobs q0jk enter a product class for the first time in the network) transition probabilities from node i to node j for class k qijk Ca0jk SCV for the interarrival rate at node j for class k SCV for the service rate at node j for class k Csjk The TPMs for all the product classes need to be specified along with the sequence of operations for each class. The net arrival rate of product k at node j is defined as the summation of its external arrival rate at node j (k0jk) and the rate at which it arrives from other nodes (which depends upon the transition probability). The net arrival rate to node j for class k can be computed using the set of simultaneous Eq. (5), which are similar to Eq. (3) X kjk ¼ k0jk q0jk þ qijk kik 8k 2 K: ð5Þ i2I k

Similarly, the number of visits made to node j by class k can be computed by solving the following set of simultaneous equations:

732

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

vjk ¼ q0jk þ

X

qijk vik

8k 2 K:

ð6Þ

i2I k

The actual arrival rate for class k (kijk) from node i to node j can be calculated using the net arrival rates calculated using Eq. (5). It is the proportion of the job arrivals (from product class k) which visit node j immediately after visiting node i. The jobs from product class k can visit node j from any other node depending upon the process route and the nature of the failure. The proportion of arrivals for product class k (pijk) from node i to node j can be calculated using Eq. (8): kijk ¼ kik qijk ; kijk pijk ¼ : kjk

ð7Þ ð8Þ

The traffic intensity for a node in a queuing network with one product class and no product failure considerations is given by the ratio of arrival rate to service rate. However, for the problem under study, depending upon the sequence of operations and the TPMs products from different classes could visit a node multiple times. Therefore, the service rates at each node need to be aggregated in order to compute the net traffic intensity at each node. The net traffic intensity for node j is the ratio of the arrival rate to the service rate considering all the jobs from different product classes which flow through the node. The aggregate service time at node j can be computed using Eq. (9): P P P kijk sjk vjk k2K k0jk sjk þ Pk2K Pi2I k sj ¼ P 8j 2 J: ð9Þ k2K k0jk þ k2K i2I k kijk vjk Eq. (9) essentially calculates the weighted average of the service time based on the external arrival rate and the actual arrival rate at node j. Eq. (9) is the modified version of its deterministic equivalent (Whitt, 1983a) shown in Eq. (10): P P k2K i2I K kk sjk F ij P sj ¼ P 8j 2 J ; ð10Þ k2K i2I k kk F ij where Fij is an indicator variable; Fij = 1, if product class k visits node j and Fij = 0, otherwise. Eq. (9) differs from Eq. (10) with respect to the stochastic flow of jobs that are represented by the TPMs of different products. The first term in the numerator and the denominator of Eq. (9) is not multiplied by vjk since it represents the external arrival process with vjk = 1. In Eq. (10), job circulations due to product failures are not considered and hence, the number of visits made by a job to a node is not taken into account (or vjk = 1). The node variability parameters (i.e. the aggregate SCV at node j) are calculated using the property that the ‘second moment of a mixture of distributions is the mixture of the second moments’ (Whitt, 1983a). Therefore, the aggregate value of Csj (due to superposition) is calculated using P P P 2 kijk s2jk ðCsjk þ 1Þvjk k2K k0jk sjk ðCsjk þ 1Þ þ 2 P Pk2K Pi2I k sj ðCsj þ 1Þ ¼ 8j 2 J : ð11Þ k2K k0jk þ k2K i2I k kijk vjk Eq. (11) calculates the aggregate SCV of the service rates by weighing the arrival rates with the square of the service times. Similar to Eq. (10), the deterministic version of Eq. (11) is given by Whitt (1983a) P P kk s2 ðCsj þ 1ÞF ij k2K 2 Pi2I k P jk sj ðCsj þ 1Þ ¼ 8j 2 J : ð12Þ k2K i2I k kk F ij The net traffic intensity at node j is given by P kjk qj ¼ k2K : sj c j

ð13Þ

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

The aggregate SCV of the external arrival process to node j is given by (  ) X k0jk Ca0jk P Ca0j ¼ ð1  wj Þ þ wj ; k2K k0jk k2K

733

ð14Þ

where wj ¼

1

ð15Þ

1 þ 4ð1  qj Þ2 ðmj  1Þ

and "

X  k0jk 2 P mj ¼ k2K k0jk k2K

#1 ð16Þ

:

Eqs. (14)–(16) are based on the hybrid approximations for the superposition process (Whitt, 1982), which are based on the following methods: • Asymptotic method. • Stationary-interval method. The external arrival process to node j is the superposition of the external arrival process to node j from different classes. The SCV of the arrival process at node j for class k is calculated using the set of simultaneous equations described in Whitt (1983a) X Cajk ¼ ajk þ Caik bijk 8k 2 K; ð17Þ i2I k

where ( ajk ¼ 1 þ w1jk ðp0jk Ca0j  1Þ þ

X

) pijk ½ð1  qijk Þ þ qijk q2i xi  ;

ð18Þ

i2I k

xi ¼ 1 þ

maxðCsi ; 0:2Þ  1 ; pffiffiffiffi ci

ð19Þ

bijk ¼ w1jk pijk qijk ð1  q2i Þ; 2

ð20Þ 1

w1jk ¼ ½1 þ 4ð1  qj Þ ðmjk  1Þ ; " #1 Ik X 2 mjk ¼ pijk ;

ð21Þ ð22Þ

i¼0

where xi, mjk, and w1jk depend upon the data determined previously i.e. qj, cj, and Csj. The variable xi and mjk are used to specify the departure operations whereas w1j is used to specify the superposition operation (Whitt, 1983a). 4.2.2. Computation of performance measures According to the Jackson theorem, a product form solution exists for the number of jobs at different nodes in the system as long as the system is Markovian. The performance measures can be computed by simply adding the WIP and queue lengths multiplied by the number of visits made by the job to that particular node. Therefore, for a series of M/M/c systems with feedback due to product failures, the expected number of jobs in the system for class k in steady state is given by Vohra (1997) X Lsk ¼ Lsik 8k 2 K; ð23Þ i2I k

734

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

where kik Lsik ¼ Lqik þ lik 2  ci kik lik

8k 2 K; qi

3

5P 0ik 8k 2 K; Lqik ¼ 4 2 ci !ð1  qi Þ " ! # D c 1 ci 1 X ðkik =lik Þ ðkik =lik Þ i P 0ik ¼ þ D! ci !ð1  qi Þ D¼0 Lq Ls P0ik

ð24Þ

ð25Þ

8k 2 K:

ð26Þ

expected number of jobs in the queue expected number of jobs in the system probability of zero jobs at node i in class k

The aggregate arrival rates and net traffic intensity used in Eqs. (24)–(26) are obtained from Eqs. (5) and (13), respectively. The expected waiting time in the system is given by X Wsk ¼ Wsik vik ; ð27Þ i2I k

where Wsik ¼ Wqik þ

1 : lik

ð28Þ

In order to calculate the expected number of jobs in the system (node), Eq. (28) could be directly used. However, it should be noted that the expected number in the queue at node i is affected by multiple product classes, since a node could serve different classes at a given time (depending upon the routes). Therefore, the average number in the queue at node j needs to be considered. The average waiting time in the queue at node i can be calculated using P Wq Wqi ¼ Pk2K ik 8i 2 J ; ð29Þ k2K Rik where

  Lqi Caik þ Csi Wqi ¼ kik 2

ð30Þ

and Rik is 1 if node i is visited by class k and 0 otherwise. Therefore, the expected waiting time at node i is given by   1 ð31Þ Wsik ¼ Wqi þ vik : lik The multiplicative factor of the average of the SCV of the arrival and service rate equals one for a perfect M/M/c type of distribution. 4.3. Experiments on multi-class M/M/c systems Multiple problem instances were generated in order to predict the performance measures of a manufacturing system with multiple product classes, product failures, and shared nodes. The purpose of this experimental study was to obtain the performance measures using the analytical formulations discussed in the aforementioned sections. Eqs. (9) and (11) were presented as the modified version of their deterministic equivalents (Eqs. (10) and (12), respectively) developed by Whitt (1983a) to capture the effect of job circulations due to product failures. Many of the original formulations developed earlier were largely based on approximations

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

735

Table 1 Overview of problem instances for Jackson networks Classes

Min nodes

Max nodes

Total runs

1 2 3 4 5

2 2 3 3 3

5 5 5 5 5

4 6 4 5 6

and simulation experiments that were refined over time. In addition, the exact nature of the system for which the approximations were developed was not described by the previous researchers. Many of these formulations were developed for a single node system unlike the complex manufacturing systems presented in this research. Consequently, it is critical to verify the accuracy of the revised formulations using simulation for the manufacturing scenario under consideration. Developing a simulation model to analyze a large manufacturing system consumes a significant amount of resources in terms of computational time and model development. Furthermore, validating and verifying these models require a significant effort. The analytical formulations discussed in earlier sections are relatively easier to solve. If the error between the analytical formulations and simulation is shown to be negligible or insignificant, these formulations can be used to estimate the performance measures in a short time. Table 2 Problem instances for M/M/c systems RUNS

NODES_1

One class 1 2 3 4

NODES_2

NODES_3

NODES_4

2 3 4 5

Two classes 5 6 7 8 9 10

2 2 2 3 3 5

2 3 5 3 5 5

Three classes 11 12 13 14

3 3 3 5

3 3 5 5

3 5 5 5

Four classes 15 16 17 18 19

3 3 3 3 5

3 3 3 5 5

3 3 5 5 5

3 5 5 5 5

Five classes 20 21 22 23 24 25

3 3 3 3 3 5

3 3 3 3 5 5

3 3 3 5 5 5

3 3 5 5 5 5

NODES_5

SN 0 0 0 0

1 3 4 2 2 4

5 4 6 4

3 5 7 8 9

3 5 5 5 5 5

9 6 8 8 7 8

736

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

Twenty five problem instances were generated to measure the accuracy of the analytical estimates. The product classes were varied from one to five. Table 1 presents an overview of these problem instances along with the minimum and maximum number of nodes in each class. The number of nodes indicates the number of process steps (internal nodes) considered in the problem. Multiple experiments were designed within each category by varying the number of nodes within the class. The number of nodes that could be shared within a network is generated from a Discrete Uniform distribution with parameters A and B, where A and B are 25% and 50% of the total number of nodes in the network for all the classes. Arrival and service rates were selected such that the service rates at all the stations (before solving the problem) were greater than the external arrival rate. Later, the aggregate service rates would take precedence over the individual service rates while computing the traffic efficiency at each node. Table 2 summarizes the characteristics of the problems that were generated. The size (number of nodes in the queuing network) of each problem instance can be determined by using the information provided in this table. For each problem instance, the number of product classes, the number of process steps for each product class, and the total number of shared nodes is listed in Table 2. The notation Nodes_X refers to the number of nodes in the Xth class in a problem instance. Experiments 1–4 are for one product class. Therefore, there are no shared nodes for these experiments. The TPMs were randomly generated for different product classes for each problem instance. Some sample TPMs for different node configurations are shown in Table 3. For each matrix, the first (and the last) row and column represent the entry (and the exit) nodes. A sample problem for three classes is shown in Fig. 4 and the corresponding inputs (for the model) in Table 4. The backward flows resulting due to product failures are not shown in Fig. 4. Note that the SCV for arrival and service rates are unity for a Markovian system. The performance measures are primarily based on the internal nodes. The simultaneous equations for the net arrival rate (Eq. (5)) and the number of visits (Eq. (6)) were independently solved and the performance measures (using Eqs. (23) and (31)) were calculated using Matlab 6.5. Simulation models were Table 3 Sample transition probability matrices (yield configurations) Two nodes 0 0 0 0

1 0.3 0.05 0

0 0.7 0.43 0

0 0 0.5 1

Three nodes 0 0 0 0 0

1 0.3 0.2 0.4 0

0 0.7 0.3 0.05 0

0 0 0.5 0.05 0

0 0 0 0.5 1

Four nodes 0 0 0 0 0 0

1 0.1 0.15 0.04 0 0

0 0.9 0.47 0.23 0 0

0 0 0.38 0.17 0 0

0 0 0 0.56 0 0

0 0 0 0 1 1

Five nodes 0 0 0 0 0 0 0

1 0.4 0.15 0.13 0.02 0 0

0 0.6 0.36 0.04 0.39 0 0

0 0 0.49 0 0.07 0 0

0 0 0 0.83 0.09 0.3 0

0 0 0 0 0.43 0.3 0

0 0 0 0 0 0.4 1

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

737

Fig. 4. Sample problem with three classes and shared nodes.

developed in Arena 5.0 for all the problem instances. Ten replications were conducted for each experiment. The SEEDS element in Arena defines the seed values for random number streams and determines how the stream should be reinitialized between replications. By default, when a new replication begins the final seed value from the previous replication is used as the new initial seed. The mean of the arrival and service times Table 4 Example of inputs to the program for analytical model (three classes)

738

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

were set equal to the values used in the analytical models. Several preliminary experiments were conducted to determine the warm-up period. These experiments were conducted to determine the run length during which the simulation model was unsteady. Our preliminary experiments indicated that a warm-up period of 100 h would discount the unsteady state for all the experimental runs. Consequently, the run length of the simulation model was set as 100 h of warm-up plus 1000 h. The product flow was represented using the SEQUENCE module. 5. Results and discussion The primary performance measures, such as WIP and flow time, were obtained from the analytical and the simulation models. The percentage error of the absolute difference between the analytical and the simulation results was calculated. The error between simulation and analytical results was analyzed by observing the trends between different factors such as total nodes in the system (TN), number of shared nodes (SN), sharing ratio (SR), and the maximum traffic intensity (q maxk) for class k. The sharing ratio is defined as the ratio of the number of shared nodes to the total number of nodes in the system Sharing ratio ¼

Shared nodes : Total nodes in the network

ð32Þ

The actual number of product classes that share a particular node was not taken into account in Eq. (32). For the sample problem shown in Table 4, nodes 2–4, and 7 are shared between the three product classes. Therefore, the sharing ratio is 4/10 or 25%. Node 4 is shared among all the classes. However, node 3 is shared only between product classes 1 and 2. Consequently, Eq. (32) does not differentiate between the nodes with different degrees of sharing. The effect is compensated by the traffic intensity at a node that would consider the effect of net arrival rates from different classes and product failures. Since it was not possible to measure the error based on the traffic intensity at each node, the maximum traffic intensity for a given class was considered for the purpose of analyzing the errors. In order to obtain the overall distribution of the errors for the two performance measures, histograms of the errors are plotted for all products classes as shown in Fig. 5. The error appears to be distributed between zero and 30%. It is evident from this figure that the error for most of the problem instances was smaller (less than 15%). Fig. 5 displays the results for individual product classes and hence, there are 35 data points although we had generated only 25 problem instances. If the accuracy of the analytical formulations were to be improved, it is required to identify the source of error. Hence, the error was analyzed by forming a logical group of problems based on the salient features of the problems such as shared nodes, total number of nodes, sharing ratio, and traffic intensity. In Fig. 6, the error is shown for all the problem classes in the increasing order of the number of shared nodes. The error appears to be increasing as the number of shared nodes increases with significant discontinuity in the trend. In Fig. 7, the average error in the performance measures is plotted by grouping the problem instances with same number of shared nodes. An increasing trend is apparent except for problems where the number of shared nodes equaled five. Further analysis of Fig. 7 was conducted to analyze the decrease in error for problems with five shared nodes. This was accomplished by analyzing the problems with 4–6 shared nodes. Besides the number of shared nodes, one of the critical factors that could affect the trends is the product yield, which is represented by the TPMs that are randomly generated during program execution. Since it is not possible to directly compare the yield of different products, the net arrival rates at all the nodes for the respective problem classes were compared. Figs. 8–10 display the trend in the total net arrival rate at each node in the network for all the product classes. Fig. 8 represents problems 10, 12, and 14 with four shared nodes; Fig. 9 represents problems 11 and 16 with five shared nodes; and Fig. 10 represents problems 13, 17, and 18 with six shared nodes. It is apparent from the graphs that problems 11 and 16 had an overall net arrival rate lower than their counterparts in Figs. 8 and 10. This implies that the yields for the product classes in problems 11 and 16 were relatively higher compared to those belonging to Figs. 8 and 10. Additionally, in Fig. 8, the net arrival rate for

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

Fig. 5. Overall distribution in the error for the performance measures.

Fig. 6. Overall trend in error for flow time and WIP against shared nodes.

739

740

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

Fig. 7. Trend in the average error in performance measures grouped by shared nodes.

Fig. 8. Total net arrival rate for problems with four shared nodes.

Fig. 9. Total net arrival rate for problems with five shared nodes.

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

741

Fig. 10. Total net arrival rate for problems with six shared nodes.

problem 10 is lesser than for problems 12 and 14 and so are the errors in the performance measures. Higher yields could be a contributing factor to the overall reduction of error in the performance measures. Fig. 11 presents the percentage error in flow time and WIP for each problem instance when the problem instances were arranged in ascending order of the total number of nodes. Fig. 12 displays the trends in error when the problem instances where grouped per the total number of nodes in the network. The total number of nodes will increase with the increase in the number of products. Fig. 11 shows that the overall error is less than 1% for problems with a single product class, and later continues to increase with significant discontinuity in the trend. The slope of the trend is lower compared to Fig. 7 for shared nodes. A significant increase in error was observed for problems with 16 nodes. Additional analysis was conducted by comparing the total net arrival rate for problems 13, 16, and 17. The design of the networks for these problems indicated a heavy concentration of product classes at fewer nodes in the network for problem instance 16 as compared to problems 13 and 17. The heavy concentration of product classes at fewer nodes, coupled with product yields, could contribute to a relatively higher error in the performance measures. The error in the performance measures appears to increase with the total number of nodes in the system. The error is less than 1% for problems with a single product class, and a significant variation in the trend is observed for higher product classes and total nodes (although the overall error appears to increase). The trend chart for sharing ratio (Fig. 13) displays a similar behavior. However, the amount of increase in error is not as

Fig. 11. Overall trend in error for flow time and WIP against total nodes.

742

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

Fig. 12. Trend in the average error in performance measures grouped by total nodes.

Fig. 13. Error trends in performance measures against sharing ratio.

significant as observed with shared nodes and total nodes. This is because the sharing ratio is normalized as the total nodes and shared nodes in the system increase. The effect of these factors on the error is evident from the trends. The trends indicate that the complexity in the problem could increase as the total number of nodes and the number of shared nodes increases in the system. The effect of maximum traffic intensity on the error in performance measures (for each run) is shown in Fig. 14. Since the traffic intensity for each node in the network could not be graphically represented for all the problem classes, the maximum traffic intensity for each product class was considered. Traffic intensity also has a significant effect on the error, which increases as the system approaches a heavy traffic scenario. Traffic intensity (by definition) captures the effects of most of the factors such as the net arrival rate (which considers the yields in the system), aggregate service rates, and the number of servers at the node with the highest traffic intensity. It should be noted that these results hold good for a system with the yields that were considered in the problem. The yield data for a given class was randomly generated while the program was executed. Depending upon the intensity of the yield and the backflows at different nodes, factors such as traffic intensity may not affect the error significantly. The claim was validated by comparing the error in performance measures

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

743

Fig. 14. Error trends in performance measures against maximum traffic intensity.

Table 5 Differences in the error in performance measures with and without product failures Problem no. 3 8 13

17

22

FT without failures (%)

FT with failures (%)

WIP without failures (%)

WIP with failures (%)

0.17 0.40 0.40 0.15 0.51 0.51 0.71 0.35 0.23 0.38 0.32 0.20 0.66 0.16 0.06

0.06 3.97 5.34 8.20 6.87 7.66 1.21 4.68 12.81 17.88 4.55 0.91 9.83 17.05 15.57

0.38 0.32 0.32 0.83 0.76 0.20 0.01 0.50 0.76 0.21 1.09 1.00 0.23 0.40 0.53

0.24 4.44 5.79 13.18 9.17 7.35 1.10 6.49 13.05 16.15 0.88 0.42 11.18 15.93 19.77

for problems without product failures (100% yields) with the problems in which the TPMs were randomly generated by the analytical model. One problem from each category (product class) was chosen for comparison. The results are shown in Table 5. The results indicate that the error in the performance measures was higher when product failures were considered in the system for most of the problem instances, which leads to the conclusion that the intensity of product yields could significantly affect the error in the performance measures. Additionally, the error could also be affected by the approximations for calculating the aggregate variability parameters at each node. Additional studies were conducted to validate the traffic equation (6). The number of visits made by different products to the different nodes was calculated using Jackson networks and compared with simulation results. A counter was used in simulation to determine how many times a product visits a node. The results were in close agreement thereby proving the validity of Eq. (6). 6. Conclusion The manufacturing scenario under study cannot be treated as a traditional flow shop. The stochastic nature of the product flow (due to product failures) and resource sharing necessitates new analytical formulations to

744

S. Pradhan et al. / European Journal of Operational Research 184 (2008) 725–744

determine the performance measures. Although analytical formulations exist for determining the flow time and WIP of an M/M/c system, the existing formulations may not be applicable to the problem under study. In this research, analytical models to compute the flow time and WIP for the complex manufacturing systems were proposed. Multiple problem instances were randomly generated and their performance measures were determined using the analytical formulations. Simulation models for these problems were also developed, and the results were compared to the analytical results. The error was less than 15% for many of the problem instances. The error between the analytical and simulation results was analyzed by logically grouping the problems with respect to some of the critical factors (such as shared nodes and total nodes in the system) that characterize the network. It was observed that error increases as the complexity of the system increases. System complexity was characterized by the total number of nodes and shared nodes in the network. Several discontinuities in these trends were observed, which were attributed to product yields and the concentration of product classes at fewer nodes in the network. The fact that the degree of error was dependent upon the intensity of the product yields was also validated. A system with high yields and lesser sharing of nodes between classes would yield accurate results when compared with simulation. References Baskett, F., Chandy, H., Muntz, M., Palacios, F., 1975. Open, closed, and mixed network of queues with different classes of customer. Journal of Association for Computer Machinery 22 (2), 248–260. Bitran, G., Sarkar, D., 1994. Throughput analysis of manufacturing networks. European Journal of Operational Research 74 (3), 448–465. Bitran, G., Tirupati, D., 1988. Multiproduct queueing networks with deterministic routing: Decomposition approach and the notion of interference. Management Science 34 (1), 75–100. Burke, P., 1956. The output of a queueing system. Operations Research 4, 699–704. Chan, R., 2004. Iterative methods in overflow queueing models. (accessed December 2004). Jackson, R., 1954. Queueing systems with phase type service. Operations Research Quarterly 5 (2), 109–120. Jackson, J., 1957. Network of waiting lines. Operations Research 5, 518–521. Jackson, J., 1963. Jobshop like queuing systems. Management Science 10, 131–142. Kelly, F., 1979. Reversibility and Stochastic Networks. Wiley Series in Probability and Mathematical Statistics, NY. Kobza, J., Ellis, K., Vittes, F., 2002. Improving throughput for an electronics assembly line using a constraints analysis methodology. Production Planning and Control 13 (3), 262–273. Koizumi, N., 2002. A queueing network model with blocking: Analysis of congested patients flows in mental health systems, Ph.D. Dissertation, University of Pennsylvania. Ku, C., Yi, C., 2002. Access control of parallel multi-server loss queues. Performance Evaluation 50 (4), 219–231. Kumar, P., 1993. Reentrant lines. Queueing Systems 13, 87–110. Medhi, J., 1991. Stochastic Models in Queueing Theory. Academic Press Inc, San Diego, CA. Morrison, J., Kumar, P., 1998. On guaranteed throughput and efficiency of closed reentrant lines. Queueing Systems 28, 33–54. Narahari, Y., Khan, L., 1995. Performance analysis of scheduling policies in reentrant manufacturing systems. Computers and Operations Research 23 (1), 37–51. Narahari, Y., Khan, L., 1996. Modeling reentrant manufacturing systems with inspection stations. Journal of Manufacturing Systems 15 (6), 367–378. Narahari, Y., Khan, L., 1998. Asymptotic loss of priority scheduling policies in closed re-entrant lines: A computational study. European Journal of Operational Research 110 (3), 586–596. Park, Y., Kim, S., Jun, C., 2000. Performance analysis of re-entrant flow ship with single job and batch machines using mean value analysis. Production Planning and Control 2 (6), 537–546. Pradhan, S., Balasubramanian, S., Arbulich, J., Srihari, K., 2001. Implementation of fiber optics assembly in a contract manufacturing facility. In: Proceedings of SMTA Optoelectronics and Telecom Revolution, Dallas, TX, pp. 28–34. Segal, M., Whitt, W., 1989. A queueing network analyzer for manufacturing. In: Proceedings of the 12th International Teletraffic Conference, vol. 2, pp. 1146–1152. Vohra, N., 1997. Quantitative Techniques in Management. Tata McGraw Hill, New Delhi, India. Whitt, W., 1982. Approximating a point process by a renewal process, I: Two basic method. Operations Research 30 (1), 125–147. Whitt, W., 1983a. Performance of a queuing network analyzer. Bell System Technical Journal 63, 1911–1979. Whitt, W., 1983b. The queueing network analyzer for manufacturing. The Bell System Technical Journal 62 (9), 2779–2815. Whitt, W., 1992. Towards better multiclass parametric decomposition approximations for open queueing networks. Annals of Operations Research 48, 221–248. Yao, D., 1994. Stochastic Modeling and Analysis of Manufacturing Systems. Springer-Verlag Publications, NY.