- Email: [email protected]
Using an empirical queueing approach to predict future flow times
Computers & Industrial Engineering 37 (1999) 809±821
www.elsevier.com/locate/dsw
Using an empirical queueing approach to predict future ¯ow times Yi-Feng Hung*, Ching-Bin Chang Department of Industrial Engineering and Engineering Management, National Tsing Hua University, Hsinchu, Taiwan
Abstract The traditional method used to describe a steady-state manufacturing system is a queueing model; whereas the common tool used to predict the future performance of a dynamic manufacturing system is a simulation model. This study proposes an empirical queueing approach to obtain the ¯ow time performance measures of a complex dynamic manufacturing system, such as semiconductor wafer fabrication. This approach is easier to implement than a simulation model and more straightforward than a queueing model. Initially, the empirical queueing curve of each work station, which de®nes the relationship between the utilization rate and the expected waiting time of a small time period, is obtained from the historical database. Then, an iterative scheme is used to predict the future system behavior. Several latest researches have reported that the prediction of future ¯ow times is important in the operation management of a complex manufacturing system. The approach proposed in this study can be easily implemented for such purposes, and the experimental results show that this approach can accurately predict the future ¯ow times. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Production management; Flow time prediction; Semiconductor manufacturing
1. Introduction and literature review In a steady-state manufacturing environment with the product mix being constant, the expected waiting times and the expected ¯ow times are time-independent. If the assumptions of a particular queueing model [3,4] are appropriate, it could be applied to obtain the long-term * Corresponding author. Tel.: +886-3-574-2939; fax: +886-3-572-2685. E-mail address: [email protected] (Y.-F. Hung). 0360-8352/99/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 3 6 0 - 8 3 5 2 ( 0 0 ) 0 0 0 1 3 - 9
810
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
average measures of the steady-state system. However, a dynamic manufacturing system, which involves time-varying product mix, time-varying expected waiting time and time-varying expected ¯ow times, would probably never enter steady state. In such an environment, queueing models are normally not appropriate for obtaining short-term measures of the system, since most queueing models assume long-run steady-state condition. An typical example of a dynamic manufacturing system is a multi-product semiconductor manufacturing facility, whose product mix shifts quickly from time to time, because of the short product life cycle and quick change in demands. A discrete-event simulation [10] is normally used to model this system and to predict the system performance. However, it takes engineers a lot of eort to build a detailed and satisfactory simulation model, which can accurately describe a complex manufacturing system. The diculties in using simulation in semiconductor manufacturing systems are discussed by Johri [9]. Our study, proposes an alternative method to predict the future ¯ow times of a system based on the idea from queueing theory. Flow time information is essential to production planning and scheduling for long ¯ow time systems, such as semiconductor manufacturing. The pioneering work on how to consider noninteger ¯ow times in production planning model is proposed by Leachman [11], whose modeling technique is now popular in today's semiconductor industry [12]. Since the planning problem deals with future production, the prediction of future ¯ow times is necessary. An iterative computation between the production planning model and simulation was proposed by Hung and Leachman [8]. They use simulation to predict the future time-dependent ¯ow times through which the production planning model can compute accurate production plan. In addition, Lu et al. [15] and Lin [13] use the predicted future ¯ow times to compute the slack, and then, apply the least-slack-rule [15,13] to dispatch the waiting jobs in queues. This method, has proved to be very eective in reducing the mean and standard deviation of ¯ow times, which is one of the top priorities of semiconductor manufacturing management. Lu et al. used simulation models to predict ¯ow times, while Lin used exponential smoothing to update waiting time forecasts. A lot is a group of the same product that travel together through the factory. A work station is a group of identical machines used to perform operations that add value to in-process material. An operation is the performance of a required processing activity on a lot by a particular work station. A route is a sequence of operations required to produce ®nished products. In semiconductor wafer fabrication, a route can consist of more than hundreds of operations. Flow time is the dierence between the time when a lot is released into shop-¯oor and the time when the lot is ®nished. A ¯ow time consists of the waiting time and processing time of each consecutive operation on the route. We assume here that transportation times are trivial, and processing times are known parameters, while waiting times are the only unknown parameters, which require estimation. Periods are small equally-spaced time intervals on time axis. Various common queueing models have curves describing the relationship between the expected waiting time and the utilization rate of the system. Tables and graphs of such can be found in [5]. However, these relationship curves are drawn based on the calculation of the mathematical equations of a particular queueing model with certain assumptions and simpli®cations. These analytical models are normally useful for system design dealing with the long-term steady-state capacity issue. However, as for the system operation (production planning and scheduling) mentioned previously, the analytical model because of its limitation
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
811
of steady-state assumption may not be appropriate for a dynamic production environment. The complexity of the semiconductor manufacturing are discussed by Johri [9]. Also, Connor et al. [2] classi®ed the complex manufacturing methods and proposed a detailed analytical model for the semiconductor manufacturing under steady-state and single-product assumptions. Connor's model is complex and requires many input parameters for computation. Thus, building the model for a dynamic and multiple-product manufacturing environment will be even more dicult, if not impossible. Normally, under such circumstances, a simulation model is used to predict the future factory status. This study proposes using an empirical utilization rate vs. waiting time curve for each work station. Previously, this curve is provided by a queueing model, which is derived under certain assumptions and simpli®cations of the system. The assumptions and simpli®cations of queueing models have been the major obstacles in their application. This empirical queueing approach bypasses this diculty of building an analytical model for a complicated actual system, and obtains the relationship curve between utilization rate and waiting time directly from the historical database. 2. Constructing the empirical curve There are two steps in constructing the empirical curve of each work station. Section 2.1 outlines the ®rst step, which collects the empirical data points; while Section 2.2 discusses the linear programming model used for constructing the empirical approximation curve, which describes the relationship between the utilization rate and expected waiting time of a work station. 2.1. Collecting the empirical data points The ®rst step involves collecting the empirical data points, which are the historical data points of utilization rate vs. waiting time. These data points can be computed from the manufacturing database, which is called WIP (Work-In-Process) Tracking system popularly implemented in semiconductor industry [17,18]. A WIP Tracking system records the arrival time at the work station, the starting time and the ending time of the processing, as well as the departure time from the work station for each operation of every lot produced in the past. A common queueing model provides the relationship between server utilization and expected waiting time. Most of the queueing models deal with a system under steady state. However, in a complex manufacturing system, owing to changing product mix, long processing route and unreliable machines, it may not be proper to assume steady state in the long run. However, it may be possible to assume steady state within a small time period, say within a half day. Even so, a queueing model is still not ideal because of three reasons. First, owing to the same reasons mentioned above, it is common that the bottleneck may shift from one work station to another. Further, the required machine utilization rate within a particular period may be more than 1.0; i.e., the workload requirement is more than the capacity available within a small time period. Most queueing models cannot deal with a utilization rate greater than 1.0 s, for most queueing models, such as M/M/1 or M/M/s, as the utilization rate approaches 1.0, the waiting time will approach in®nity. Third, Hou [6] pointed out that, based on his experiment, the
812
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
queueing model cannot predict waiting time accurately in a dynamic manufacturing environment. A waiting time is de®ned as the time between the time point of arrival and the time point of operation initiated. A utilization rate is de®ned as the workload required in a period divided by the number of machines and the length of a period. The empirical curve of a particular work station plays the same role as the relationship curve between utilization rate and expected waiting time provided by a queueing model. However, a queueing model is usually used in a long-term steady state; whereas, the empirical curve here is intended for use in short time periods. We also implicitly assume that the waiting time at the work station is independent on the operations to be processed and the product type of the lots. To describe how to obtain such empirical curve the following notations are used: t k b nk akt wkt ukt qkt pk ckt x kt ykt
time period index; t 1, 2, . . . ,T; work station index; k 1, 2, . . . ,K; the length (hours) of a small time period; the number of identical machines in work station k; the cumulative total workload (machine hours) arrival at work station k in period t; the cumulative waiting time of lots arriving at work station k in period t; the cumulative total number of lots arriving at work station k in period t; the amount of the queue (machine hours) of work station k at the end of period t; the up time ratio of work station k; the capacity (machine hours) of work station k in period t; the utilization rate of work station k in period t; the average waiting time of lots arriving at work station k in period t.
The following procedures are needed to obtain the historical data: 1. Initialize variables. set akt 0, wkt 0, and ukt 0, 8k 8t: 2. From a WIP Tracking database system, we can collect the arrival time at work stations of all lots and their processing times. Based on these data, we can compute the cumulative arrival workload akt s, and their corresponding cumulative total waiting times wkt s. While (there are lot arrival data left) remove and use the next lot arrival data of the lot e / compute the time period number / t d the arrival time b akt akt the processing time of the lot wkt wkt the waiting time of the lot ukt ukt 1 end while loop 3. We assume the capacity to be period-independent and estimate the average capacity of work station k by: ckt nk b pk ,
8t 8k:
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
813
4. We can estimate the queue amount at the end of each time period. The queue at the end of a period (the start of next period) is equal to the queue at the start of the period plus the arrival workload during the period minus the capacity during the period. However, it is impossible to have a negative queue. For k 1 to K For t 1 to T ÿ qkt max 0, qk, tÿ1 akt ÿ ckt : 5. For a particular work station k, the total workload of a particular period is the sum of initial queue and the arrival workload within the period; thus, the utilization rate can be computed as follows: x kt
qk, tÿ1 akt , nk b
8t 8k:
6. For all the lots arriving at work station k in period t, we can compute their average waiting time: ykt
wkt , 8t 8k: ukt
2.2. Fitting the empirical approximate curves Here, we have a data point x kt , ykt for each past period t of each work station k. For work station k, we can use the following linear programming model LPk and data points x kt , ykt , t 1, . . . ,T to ®t an empirical piecewise linear curve, as shown in Fig. 1. 1. Index: t: data point number (previously, period number); i: break point number on the piecewise linear approximation curve. 2. Parameters: x t = the utilization rate of data point t; it equals to x kt for LPk ; yt = the average waiting time of data point t; it equals to ykt for LPk ; max utilization rate max x t ; t
max waiting time max yt ; t
ri = the utilization rate of the ith break point, i 1, . . . ,H, where r0 0, r0 < r1 < r2 < < rH and rHÿ1 < max utilization rateRrH : 3. Decision variables: Ei = the expected waiting time of the ith break point, whose utilization rate is ri :
814
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
Fig. 1. The linear programming model used to ®t the piecewise linear curve.
4. Notations: yt = the y-value on the piecewise linear curve, when x-value is x t : 5. Objective function: T X jyt ÿ y^t j Minimize t1
We seek to minimize the sum of the absolute deviation between the waiting time projected by the piecewise linear curve y^t and the actual waiting time yt of each data point. 6. Constraints: (1) The y-value projection of x t on piecewise curve: y^t Ei
Ei1 ÿ Ei x t ÿ ri , in which i satisfies ri Rx t < ri1 , 8t: ri1 ÿ ri
(2) The slopes are non-decreasing as the utilization rate increases: Ei1 ÿ Ei Ei ÿ Eiÿ1 r , i 1, 2, ,H ÿ 1: ri1 ÿ ri ri ÿ riÿ1 (3) The maximum waiting time constraint: EH Rmax waiting time: (4) The slopes are non-negative: E1 rE0 : (5) The non-negativity constraints: Ei r0,
8i 0, 1, 2, ,H:
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
815
After solving the LPk , we have the break points ri , Ei ), i 1, 2, ,H: These points empirically depict the relationship between the utilization rate and waiting time for work station k. We believe that this curve is dependent only on the manufacturing method of the equipment, the number of machines in work station, and the machine failure pattern (the distribution of the time-between-failure and the distribution of the time-to-repair). Therefore, this curve need to be reconstructed only when the above dependent factors have changed. We do not need to resolve the LP every time when we do the iterative computations. In the above formulation, ri ÿ riÿ1 is the x-distance between two consecutive break points. The smaller the distance, the more break points, the more accuracy of the empirical curve. However, the longer it takes to compute in both obtaining and applying the curve. In our experiment, we use 0.04 as the x-distance of two consecutive break points. We believe 25 points is sucient to approximate the curve when the utilization rate is between 0.0 and 1.0. It took less than 2 min on a SUN UNIX computer work station to solve one LP problem, which covers the operation data of a work station for three years. 3. Iterative procedure to predict future ¯ow time Machine loading (utilization rate) will aect waiting time, which in turn determines the time of lot arrival at the following work station; that is, the machine loading of the following work station is changed. Then, the waiting time of the following work station is also changed. The determination of waiting times and utilization rates are dependent on each other; thus, these chain eects form a vicious circle. In order to resolve this problem, an iterative approach is proposed here. Initially, we assume the expected waiting time in each period on each work station to be zero. By using the processing times, we can then determine the arrival time at the work station of each operation for each lot, from which we can compute the arrival workload (machine hours) of each time period at each work station. Then, we can compute the estimated queue at the end points of periods, and the utilization rate in each period of each work station. Applying the piecewise curve obtained previously, we then compute a projected waiting time in each period of each work station from the utilization rate. At the end of an iteration, an exponential smoothing formula is used to update expected waiting times. The reason of using exponential smoothing is explained in Section 4. Once we have the new expected waiting time in each period of each work station, we can start another iteration. That is to recalculate the arrival time at the work station of each operation for each lot; thus, revising the arrival workload in each period of each work station. Then, the revised queue, the revised utilization rate, and revised expected waiting time can be computed. Through a series of iterations, we could converge to an accurate ¯ow time prediction for each lot. The following is the outline of the iterative procedure: Notations s variable for the iteration number S maximum number of iterations l variable for lot identi®cation
816
e r U
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
time point variable variable for operation identi®cation a set of lots
Begin wkt 0 8k, 8t for s 1 to S akt 0 8k, 8t x kt 0 8k, 8t put all the planned lots and WIP lots into set U while (set U is not empty) l = remove a lot from set U if (l is a WIP lot) e current time else e the planned release time endif while (there is operation left on l ) r = remove next operation of lot l k = the work station performing r t b be c / operation r start waiting on work station k in period t / e e wkt / delay for waiting time / e e + the processing time of r / delay for processing time / akt akt + the processing time of r / cumulate the workload / end while loop / operation / end while loop / lot / / Note: here we have the revised arrival workload akt array / for k 1 to K for t 1 to T qkt maxf0; qk; tÿ1 akt ÿckt g / revise queue qkt / x kt qk; tÿ1 akt nk b / revise utilization rate x kt / wtemp = the projected waiting time from utilization rate x kt using the empirical curve of work station k wkt wkt a wtemp ÿ wkt ; /use exponential smoothing to revise wkt / end for t end for k end for s End
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
817
4. Experiment A series of experiments have been done to validate our approach. We started with the determination of the parameters used; then, compared the ¯ow times predicted by our approach with those from a simulation model.
4.1. Determining the parameters This section focus on how to determine the values for a, b and S? We used the same data set as [7,8], in which there are 30 work stations and 10 product types, each of which follows a prede®ned distinct route. The average number of operations on a route is more than 130 operations. A particular machine type may be visited by the same lot several times. We varied the exponential smoothing factor a from 0.1, 0.2, 0.3,. . . through 1.0, and the length of small period b from 0.25, 0.5, 1.0, through 2.0 days, and we iterated 15 times. Let us assume that N is the total number of lots, and gsl the computed ¯ow time of lot l in iteration s. The closeness of two consecutive iterations is often used as a convergence criterion in numerical methods [16]. To measure this closeness, we de®ne ( ) N X gsl ÿ gsÿ1, l =N 100%: Ds gsÿ1, l l1 Fig. 2 shows the Ds values of various a values when b 0:5 days. Letting a 1:0 (no exponential smoothing) is equivalent to ®xed point iteration [16]. According to our experiment, when the a is larger than 0.7 (including 1.0), the Ds curves are too high (thus, not shown in the ®gure), and the iteration scheme could not converge. We conjectured that the function which governs the relationship between the machine loading of all machines and ¯ow times of all lots is not continuous, and is very complicated. The traditional ®xed point iteration would overreact to the changes between two iterations; thus, it could not converge. It can be seen that at iteration 9, the Ds value does not decrease greatly with increasing number of iterations; i.e., the iteration scheme converges after about nine iterations. Thus, in the following experiment, we set the number of iteration S 9: To measure the accuracy of the piecewise curve iteration approach, a simulation model that involved random machine failures was used. Thirty simulation runs with dierent initial seeds were executed to collect the mean ¯ow time of each individual lot. Let fl be the average simulated ¯ow time of lot l. These simulated ¯ow times were compared with those computed by the iteration approach. Two metrics were used to determine the accuracy of the iterative scheme with s iterations. They were the average ¯ow time dierence of iteration s, de®ned as ( ) N X gsl ÿ fl =N 100%, Ms fl l1 and the average absolute ¯ow time dierence of iteration s, de®ned as
818
( As
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
) N X jgsl ÿ fl j =N 100%: fl l1
Let p be the product type index, Fp the average ¯ow time of product type p from simulation, sp the standard deviation of ¯ow time of product type p from simulation, and rp the number of P lots that are product type p; thus, N p rp : A weighted average percentage of standard deviation of ¯ow time is computed to estimate the con®dence level for our ¯ow time estimation: X sp rp : SD 100% Fp N p Figs. 3 and 4 show the M9 and A9 values, respectively, of various a and b values. It can be seen that the performances are all satisfactory when a is between 0.2 and 0.7 and b is less than 2 days. However, in the following experiments, we set a 0:5 and b 0:5 days, whose M9 value is almost equal to zero. 4.2. Experimental validation This study used the production planning model proposed by Leachman [11] to obtain the monthly production plan within a planning horizon of two years. As long as there is demand, this model will fully utilize the resource capacity, which is the case in our experiment. Since the
Fig. 2. Ds values of various a values.
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
819
Fig. 3. M9 values of various a and b values.
demand data are time-varying, the planned production mix also varies from month to month. We called this production plan ``changing product mix''. If we sum up the planned production quantity of all months within the two-year planning horizon of a particular product type, and then, evenly release it within the planning horizon, we called this case the ``®xed product mix''. We can clearly see that ``®xed product mix'' production plan will drive the manufacturing system into steady state after a warm-up period; whereas, under the ``changing product mix'' production plan and the ¯ow time being as long as months, the system will not enter a steady
Fig. 4. A9 values of various a and b values.
820
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
Table 1 A comparison of the simulation ¯ow time and the iterative ¯ow time in the case of ®xed product mix (%) Release scale
80%
85%
90%
95%
100%
M9 A9 D9 SD
2.87% 3.87% 0.10% 4.36%
1.98% 3.23% 0.08% 4.98%
1.21% 2.79% 0.10% 4.87%
0.72% 2.83% 0.21% 5.08%
0.22% 3.06% 0.25% 5.04%
state, and thus, being a dynamic system. In order to experiment the various level of machine utilization rate (due to the various level of demand rates), we scaled the production quantity of each month by various percentages: 80%, 85%, 90%, 95%, and 100%. Owing to the unavailability of the actual factory waiting time data, simulation was used to generate the required data for constructing the relationship curve in our experiment. The ¯ow times obtained by the simulation of various production schedules were compared with those predicted by the iterative approach. Tables 1 and 2 show the case of ®xed product mix and changing product mix, respectively. From these tables, we can see that the dierences between the ¯ow times obtained by simulation and the ¯ow times predicted by the iterative approach was very small with mean absolute dierence less than 5%. They were all within one standard deviation of the simulated ¯ow time. If assuming normal distribution and using 95% con®dence interval, whose boundary are plus and minus two standard deviation away from the mean value, the ¯ow times estimated by the iterative approach were well within the con®dence interval. The risk of estimation error for our iterative approach is very small. The result shows that our approach performs well in both steady state environment (®xed product mix) and dynamic environment (changing product mix); and its performance is as good as that of the queueing models with a steady state assumption [1,2], which use complex mathematical models to derive measures. We used a SUN UNIX work station to run our experiments. The execution time of our proposed iterative approach is less than 10 s/iteration. Thus, it took about 1.5 min to converge in nine iterations. The simulation programs were written in C++ [14]. It took about 10 min/ simulation run and about 5 h to run 30 simulation runs to obtain the average ¯ow times. The time of a simulation run is about 60 times that of an iteration run in our approach.
Table 2 A comparison of the simulation ¯ow time and the iterative ¯ow time in the case of changing product mix (%) Release scale
80%
85%
90%
95%
100%
M9 A9 D9 SD
2.42% 3.66% 0.16% 6.49%
2.02% 3.59% 0.25% 6.65%
1.55% 4.62% 0.33% 7.14%
1.86% 4.13% 0.68% 8.58%
0.19% 5.10% 0.89% 15.65%
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
821
5. Conclusion and future direction Our approach bypasses the diculty of the mathematical modeling for various complex manufacturing methods. All analytical models aim to obtain the relationship between the machine utilization rate and the waiting time, which in this study are empirically obtained by direct construction from the historical data of machine utilization rate and waiting time. Based on the simulation experiment, our method works well in predicting the ¯ow times, which are important to the production planning and schedule for complex manufacturing system like semiconductor manufacturing. The industrial application of the iterative approach in realistic production plan and schedule problems requires further study. References [1] Chen H, Harrison JM, Mandelbaum A, Van Ackere A, Wein LM. Empirical evaluation of a queueing network model for semiconductor wafer fabrication. Operations Research 1988;36(2):202±15. [2] Connors DP, Feigin GE, Yao DD. A queueing newtwork model for semiconductor manufacturing. IEEE Transactions on Semiconductor Manufacturing 1996;9(3):412±27. [3] Gross D, Harris CM. Fundamentals of queueing theory, 2nd ed. USA: Wiley, 1985. [4] Hillier FS, Lieberman GJ. Introduction to operation research, 6th ed. New York: McGraw-Hill, 1995. [5] Hillier FS, Yu OS, Avis D, Fossett L, Lo F, Reiman M. Queueing tables and graphs. North-Holland, New York: Elsevier, 1981. [6] Hou M-C. Using queueing models to predict the cycle time of production planning model. Master Thesis, Department of Industrial Engineering, National Tsing Hua University, Hsinchu, Taiwan, 1995. [7] Hung Y-F. Corporate-level production planning with simulation feedback of parameters. Ph.D. Dissertation, College of Engineering, University of California, Berkeley, CA, 1991. [8] Hung Y-F, Leachman RC. A production planning methodology for semiconductor manufacturing based on iterative simulation and linear programming calculations. IEEE Transactions on Semiconductor Manufacturing 1996;9(2):257±69. [9] Johri PK. Practical issues in scheduling and dispatching in semiconductor wafer fabrication. Journal of Manufacturing Systems 1993;12(6):474±85. [10] Law AM, Kelton WD. Simulation modeling and analysis. New York: McGraw-Hill, 1991. [11] Leachman RC. Modeling technique for automated production planning in the semiconductor industry. In: Ciriani TA, Leachman RC, editors. Optimization in industry. England: Wiley, 1993. p. 1±30. [12] Leachman RC, Benson RF, Liu C, Raar DJ. IMPReSS: an automated production-planning and deliveryquotation system at Harris corporation Ð semiconductor sector. Interface 1996;26(1):6±37. [13] Lin C-Y, Shop ¯oor scheduling of semiconductor wafer fabrication using real-time feedback control and prediction. Ph.D. Dissertation, Engineering±Industrial Engineering and Operations Research, University of California at Berkeley, 1996. [14] Lippman SB. C++ primer, 2nd ed. MA, USA: Addison-Wesley, 1991. [15] Lu SCH, Ramaswamy D, Kumar PR. Ecient scheduling policies to reduce mean and variance of cycle-time in semiconductor manufacturing plants. IEEE Transactions on Semiconductor Manufacturing 1994;7(3):374±88. [16] Mathews JH. Numerical methods for mathematics, science, and engineering. Englewood Clis, NJ: PrenticeHall, 1992. [17] Uzsoy R, Lee CY, Martin-Vega LA. A review of production planning and scheduling models in the semiconductor industry. Part I: system characteristics, performance evaluation and production planning. IIE Transactions 1992;24(4):47±60. [18] Uzsoy R, Lee CY, Martin-Vega LA. A review of production planning and scheduling models in the semiconductor industry. Part II: shop-¯oor control. IIE Transactions 1994;26(5):44±55.
www.elsevier.com/locate/dsw
Using an empirical queueing approach to predict future ¯ow times Yi-Feng Hung*, Ching-Bin Chang Department of Industrial Engineering and Engineering Management, National Tsing Hua University, Hsinchu, Taiwan
Abstract The traditional method used to describe a steady-state manufacturing system is a queueing model; whereas the common tool used to predict the future performance of a dynamic manufacturing system is a simulation model. This study proposes an empirical queueing approach to obtain the ¯ow time performance measures of a complex dynamic manufacturing system, such as semiconductor wafer fabrication. This approach is easier to implement than a simulation model and more straightforward than a queueing model. Initially, the empirical queueing curve of each work station, which de®nes the relationship between the utilization rate and the expected waiting time of a small time period, is obtained from the historical database. Then, an iterative scheme is used to predict the future system behavior. Several latest researches have reported that the prediction of future ¯ow times is important in the operation management of a complex manufacturing system. The approach proposed in this study can be easily implemented for such purposes, and the experimental results show that this approach can accurately predict the future ¯ow times. 7 2000 Elsevier Science Ltd. All rights reserved. Keywords: Production management; Flow time prediction; Semiconductor manufacturing
1. Introduction and literature review In a steady-state manufacturing environment with the product mix being constant, the expected waiting times and the expected ¯ow times are time-independent. If the assumptions of a particular queueing model [3,4] are appropriate, it could be applied to obtain the long-term * Corresponding author. Tel.: +886-3-574-2939; fax: +886-3-572-2685. E-mail address: [email protected] (Y.-F. Hung). 0360-8352/99/$ - see front matter 7 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 3 6 0 - 8 3 5 2 ( 0 0 ) 0 0 0 1 3 - 9
810
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
average measures of the steady-state system. However, a dynamic manufacturing system, which involves time-varying product mix, time-varying expected waiting time and time-varying expected ¯ow times, would probably never enter steady state. In such an environment, queueing models are normally not appropriate for obtaining short-term measures of the system, since most queueing models assume long-run steady-state condition. An typical example of a dynamic manufacturing system is a multi-product semiconductor manufacturing facility, whose product mix shifts quickly from time to time, because of the short product life cycle and quick change in demands. A discrete-event simulation [10] is normally used to model this system and to predict the system performance. However, it takes engineers a lot of eort to build a detailed and satisfactory simulation model, which can accurately describe a complex manufacturing system. The diculties in using simulation in semiconductor manufacturing systems are discussed by Johri [9]. Our study, proposes an alternative method to predict the future ¯ow times of a system based on the idea from queueing theory. Flow time information is essential to production planning and scheduling for long ¯ow time systems, such as semiconductor manufacturing. The pioneering work on how to consider noninteger ¯ow times in production planning model is proposed by Leachman [11], whose modeling technique is now popular in today's semiconductor industry [12]. Since the planning problem deals with future production, the prediction of future ¯ow times is necessary. An iterative computation between the production planning model and simulation was proposed by Hung and Leachman [8]. They use simulation to predict the future time-dependent ¯ow times through which the production planning model can compute accurate production plan. In addition, Lu et al. [15] and Lin [13] use the predicted future ¯ow times to compute the slack, and then, apply the least-slack-rule [15,13] to dispatch the waiting jobs in queues. This method, has proved to be very eective in reducing the mean and standard deviation of ¯ow times, which is one of the top priorities of semiconductor manufacturing management. Lu et al. used simulation models to predict ¯ow times, while Lin used exponential smoothing to update waiting time forecasts. A lot is a group of the same product that travel together through the factory. A work station is a group of identical machines used to perform operations that add value to in-process material. An operation is the performance of a required processing activity on a lot by a particular work station. A route is a sequence of operations required to produce ®nished products. In semiconductor wafer fabrication, a route can consist of more than hundreds of operations. Flow time is the dierence between the time when a lot is released into shop-¯oor and the time when the lot is ®nished. A ¯ow time consists of the waiting time and processing time of each consecutive operation on the route. We assume here that transportation times are trivial, and processing times are known parameters, while waiting times are the only unknown parameters, which require estimation. Periods are small equally-spaced time intervals on time axis. Various common queueing models have curves describing the relationship between the expected waiting time and the utilization rate of the system. Tables and graphs of such can be found in [5]. However, these relationship curves are drawn based on the calculation of the mathematical equations of a particular queueing model with certain assumptions and simpli®cations. These analytical models are normally useful for system design dealing with the long-term steady-state capacity issue. However, as for the system operation (production planning and scheduling) mentioned previously, the analytical model because of its limitation
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
811
of steady-state assumption may not be appropriate for a dynamic production environment. The complexity of the semiconductor manufacturing are discussed by Johri [9]. Also, Connor et al. [2] classi®ed the complex manufacturing methods and proposed a detailed analytical model for the semiconductor manufacturing under steady-state and single-product assumptions. Connor's model is complex and requires many input parameters for computation. Thus, building the model for a dynamic and multiple-product manufacturing environment will be even more dicult, if not impossible. Normally, under such circumstances, a simulation model is used to predict the future factory status. This study proposes using an empirical utilization rate vs. waiting time curve for each work station. Previously, this curve is provided by a queueing model, which is derived under certain assumptions and simpli®cations of the system. The assumptions and simpli®cations of queueing models have been the major obstacles in their application. This empirical queueing approach bypasses this diculty of building an analytical model for a complicated actual system, and obtains the relationship curve between utilization rate and waiting time directly from the historical database. 2. Constructing the empirical curve There are two steps in constructing the empirical curve of each work station. Section 2.1 outlines the ®rst step, which collects the empirical data points; while Section 2.2 discusses the linear programming model used for constructing the empirical approximation curve, which describes the relationship between the utilization rate and expected waiting time of a work station. 2.1. Collecting the empirical data points The ®rst step involves collecting the empirical data points, which are the historical data points of utilization rate vs. waiting time. These data points can be computed from the manufacturing database, which is called WIP (Work-In-Process) Tracking system popularly implemented in semiconductor industry [17,18]. A WIP Tracking system records the arrival time at the work station, the starting time and the ending time of the processing, as well as the departure time from the work station for each operation of every lot produced in the past. A common queueing model provides the relationship between server utilization and expected waiting time. Most of the queueing models deal with a system under steady state. However, in a complex manufacturing system, owing to changing product mix, long processing route and unreliable machines, it may not be proper to assume steady state in the long run. However, it may be possible to assume steady state within a small time period, say within a half day. Even so, a queueing model is still not ideal because of three reasons. First, owing to the same reasons mentioned above, it is common that the bottleneck may shift from one work station to another. Further, the required machine utilization rate within a particular period may be more than 1.0; i.e., the workload requirement is more than the capacity available within a small time period. Most queueing models cannot deal with a utilization rate greater than 1.0 s, for most queueing models, such as M/M/1 or M/M/s, as the utilization rate approaches 1.0, the waiting time will approach in®nity. Third, Hou [6] pointed out that, based on his experiment, the
812
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
queueing model cannot predict waiting time accurately in a dynamic manufacturing environment. A waiting time is de®ned as the time between the time point of arrival and the time point of operation initiated. A utilization rate is de®ned as the workload required in a period divided by the number of machines and the length of a period. The empirical curve of a particular work station plays the same role as the relationship curve between utilization rate and expected waiting time provided by a queueing model. However, a queueing model is usually used in a long-term steady state; whereas, the empirical curve here is intended for use in short time periods. We also implicitly assume that the waiting time at the work station is independent on the operations to be processed and the product type of the lots. To describe how to obtain such empirical curve the following notations are used: t k b nk akt wkt ukt qkt pk ckt x kt ykt
time period index; t 1, 2, . . . ,T; work station index; k 1, 2, . . . ,K; the length (hours) of a small time period; the number of identical machines in work station k; the cumulative total workload (machine hours) arrival at work station k in period t; the cumulative waiting time of lots arriving at work station k in period t; the cumulative total number of lots arriving at work station k in period t; the amount of the queue (machine hours) of work station k at the end of period t; the up time ratio of work station k; the capacity (machine hours) of work station k in period t; the utilization rate of work station k in period t; the average waiting time of lots arriving at work station k in period t.
The following procedures are needed to obtain the historical data: 1. Initialize variables. set akt 0, wkt 0, and ukt 0, 8k 8t: 2. From a WIP Tracking database system, we can collect the arrival time at work stations of all lots and their processing times. Based on these data, we can compute the cumulative arrival workload akt s, and their corresponding cumulative total waiting times wkt s. While (there are lot arrival data left) remove and use the next lot arrival data of the lot e / compute the time period number / t d the arrival time b akt akt the processing time of the lot wkt wkt the waiting time of the lot ukt ukt 1 end while loop 3. We assume the capacity to be period-independent and estimate the average capacity of work station k by: ckt nk b pk ,
8t 8k:
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
813
4. We can estimate the queue amount at the end of each time period. The queue at the end of a period (the start of next period) is equal to the queue at the start of the period plus the arrival workload during the period minus the capacity during the period. However, it is impossible to have a negative queue. For k 1 to K For t 1 to T ÿ qkt max 0, qk, tÿ1 akt ÿ ckt : 5. For a particular work station k, the total workload of a particular period is the sum of initial queue and the arrival workload within the period; thus, the utilization rate can be computed as follows: x kt
qk, tÿ1 akt , nk b
8t 8k:
6. For all the lots arriving at work station k in period t, we can compute their average waiting time: ykt
wkt , 8t 8k: ukt
2.2. Fitting the empirical approximate curves Here, we have a data point x kt , ykt for each past period t of each work station k. For work station k, we can use the following linear programming model LPk and data points x kt , ykt , t 1, . . . ,T to ®t an empirical piecewise linear curve, as shown in Fig. 1. 1. Index: t: data point number (previously, period number); i: break point number on the piecewise linear approximation curve. 2. Parameters: x t = the utilization rate of data point t; it equals to x kt for LPk ; yt = the average waiting time of data point t; it equals to ykt for LPk ; max utilization rate max x t ; t
max waiting time max yt ; t
ri = the utilization rate of the ith break point, i 1, . . . ,H, where r0 0, r0 < r1 < r2 < < rH and rHÿ1 < max utilization rateRrH : 3. Decision variables: Ei = the expected waiting time of the ith break point, whose utilization rate is ri :
814
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
Fig. 1. The linear programming model used to ®t the piecewise linear curve.
4. Notations: yt = the y-value on the piecewise linear curve, when x-value is x t : 5. Objective function: T X jyt ÿ y^t j Minimize t1
We seek to minimize the sum of the absolute deviation between the waiting time projected by the piecewise linear curve y^t and the actual waiting time yt of each data point. 6. Constraints: (1) The y-value projection of x t on piecewise curve: y^t Ei
Ei1 ÿ Ei x t ÿ ri , in which i satisfies ri Rx t < ri1 , 8t: ri1 ÿ ri
(2) The slopes are non-decreasing as the utilization rate increases: Ei1 ÿ Ei Ei ÿ Eiÿ1 r , i 1, 2, ,H ÿ 1: ri1 ÿ ri ri ÿ riÿ1 (3) The maximum waiting time constraint: EH Rmax waiting time: (4) The slopes are non-negative: E1 rE0 : (5) The non-negativity constraints: Ei r0,
8i 0, 1, 2, ,H:
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
815
After solving the LPk , we have the break points ri , Ei ), i 1, 2, ,H: These points empirically depict the relationship between the utilization rate and waiting time for work station k. We believe that this curve is dependent only on the manufacturing method of the equipment, the number of machines in work station, and the machine failure pattern (the distribution of the time-between-failure and the distribution of the time-to-repair). Therefore, this curve need to be reconstructed only when the above dependent factors have changed. We do not need to resolve the LP every time when we do the iterative computations. In the above formulation, ri ÿ riÿ1 is the x-distance between two consecutive break points. The smaller the distance, the more break points, the more accuracy of the empirical curve. However, the longer it takes to compute in both obtaining and applying the curve. In our experiment, we use 0.04 as the x-distance of two consecutive break points. We believe 25 points is sucient to approximate the curve when the utilization rate is between 0.0 and 1.0. It took less than 2 min on a SUN UNIX computer work station to solve one LP problem, which covers the operation data of a work station for three years. 3. Iterative procedure to predict future ¯ow time Machine loading (utilization rate) will aect waiting time, which in turn determines the time of lot arrival at the following work station; that is, the machine loading of the following work station is changed. Then, the waiting time of the following work station is also changed. The determination of waiting times and utilization rates are dependent on each other; thus, these chain eects form a vicious circle. In order to resolve this problem, an iterative approach is proposed here. Initially, we assume the expected waiting time in each period on each work station to be zero. By using the processing times, we can then determine the arrival time at the work station of each operation for each lot, from which we can compute the arrival workload (machine hours) of each time period at each work station. Then, we can compute the estimated queue at the end points of periods, and the utilization rate in each period of each work station. Applying the piecewise curve obtained previously, we then compute a projected waiting time in each period of each work station from the utilization rate. At the end of an iteration, an exponential smoothing formula is used to update expected waiting times. The reason of using exponential smoothing is explained in Section 4. Once we have the new expected waiting time in each period of each work station, we can start another iteration. That is to recalculate the arrival time at the work station of each operation for each lot; thus, revising the arrival workload in each period of each work station. Then, the revised queue, the revised utilization rate, and revised expected waiting time can be computed. Through a series of iterations, we could converge to an accurate ¯ow time prediction for each lot. The following is the outline of the iterative procedure: Notations s variable for the iteration number S maximum number of iterations l variable for lot identi®cation
816
e r U
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
time point variable variable for operation identi®cation a set of lots
Begin wkt 0 8k, 8t for s 1 to S akt 0 8k, 8t x kt 0 8k, 8t put all the planned lots and WIP lots into set U while (set U is not empty) l = remove a lot from set U if (l is a WIP lot) e current time else e the planned release time endif while (there is operation left on l ) r = remove next operation of lot l k = the work station performing r t b be c / operation r start waiting on work station k in period t / e e wkt / delay for waiting time / e e + the processing time of r / delay for processing time / akt akt + the processing time of r / cumulate the workload / end while loop / operation / end while loop / lot / / Note: here we have the revised arrival workload akt array / for k 1 to K for t 1 to T qkt maxf0; qk; tÿ1 akt ÿckt g / revise queue qkt / x kt qk; tÿ1 akt nk b / revise utilization rate x kt / wtemp = the projected waiting time from utilization rate x kt using the empirical curve of work station k wkt wkt a wtemp ÿ wkt ; /use exponential smoothing to revise wkt / end for t end for k end for s End
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
817
4. Experiment A series of experiments have been done to validate our approach. We started with the determination of the parameters used; then, compared the ¯ow times predicted by our approach with those from a simulation model.
4.1. Determining the parameters This section focus on how to determine the values for a, b and S? We used the same data set as [7,8], in which there are 30 work stations and 10 product types, each of which follows a prede®ned distinct route. The average number of operations on a route is more than 130 operations. A particular machine type may be visited by the same lot several times. We varied the exponential smoothing factor a from 0.1, 0.2, 0.3,. . . through 1.0, and the length of small period b from 0.25, 0.5, 1.0, through 2.0 days, and we iterated 15 times. Let us assume that N is the total number of lots, and gsl the computed ¯ow time of lot l in iteration s. The closeness of two consecutive iterations is often used as a convergence criterion in numerical methods [16]. To measure this closeness, we de®ne ( ) N X gsl ÿ gsÿ1, l =N 100%: Ds gsÿ1, l l1 Fig. 2 shows the Ds values of various a values when b 0:5 days. Letting a 1:0 (no exponential smoothing) is equivalent to ®xed point iteration [16]. According to our experiment, when the a is larger than 0.7 (including 1.0), the Ds curves are too high (thus, not shown in the ®gure), and the iteration scheme could not converge. We conjectured that the function which governs the relationship between the machine loading of all machines and ¯ow times of all lots is not continuous, and is very complicated. The traditional ®xed point iteration would overreact to the changes between two iterations; thus, it could not converge. It can be seen that at iteration 9, the Ds value does not decrease greatly with increasing number of iterations; i.e., the iteration scheme converges after about nine iterations. Thus, in the following experiment, we set the number of iteration S 9: To measure the accuracy of the piecewise curve iteration approach, a simulation model that involved random machine failures was used. Thirty simulation runs with dierent initial seeds were executed to collect the mean ¯ow time of each individual lot. Let fl be the average simulated ¯ow time of lot l. These simulated ¯ow times were compared with those computed by the iteration approach. Two metrics were used to determine the accuracy of the iterative scheme with s iterations. They were the average ¯ow time dierence of iteration s, de®ned as ( ) N X gsl ÿ fl =N 100%, Ms fl l1 and the average absolute ¯ow time dierence of iteration s, de®ned as
818
( As
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
) N X jgsl ÿ fl j =N 100%: fl l1
Let p be the product type index, Fp the average ¯ow time of product type p from simulation, sp the standard deviation of ¯ow time of product type p from simulation, and rp the number of P lots that are product type p; thus, N p rp : A weighted average percentage of standard deviation of ¯ow time is computed to estimate the con®dence level for our ¯ow time estimation: X sp rp : SD 100% Fp N p Figs. 3 and 4 show the M9 and A9 values, respectively, of various a and b values. It can be seen that the performances are all satisfactory when a is between 0.2 and 0.7 and b is less than 2 days. However, in the following experiments, we set a 0:5 and b 0:5 days, whose M9 value is almost equal to zero. 4.2. Experimental validation This study used the production planning model proposed by Leachman [11] to obtain the monthly production plan within a planning horizon of two years. As long as there is demand, this model will fully utilize the resource capacity, which is the case in our experiment. Since the
Fig. 2. Ds values of various a values.
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
819
Fig. 3. M9 values of various a and b values.
demand data are time-varying, the planned production mix also varies from month to month. We called this production plan ``changing product mix''. If we sum up the planned production quantity of all months within the two-year planning horizon of a particular product type, and then, evenly release it within the planning horizon, we called this case the ``®xed product mix''. We can clearly see that ``®xed product mix'' production plan will drive the manufacturing system into steady state after a warm-up period; whereas, under the ``changing product mix'' production plan and the ¯ow time being as long as months, the system will not enter a steady
Fig. 4. A9 values of various a and b values.
820
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
Table 1 A comparison of the simulation ¯ow time and the iterative ¯ow time in the case of ®xed product mix (%) Release scale
80%
85%
90%
95%
100%
M9 A9 D9 SD
2.87% 3.87% 0.10% 4.36%
1.98% 3.23% 0.08% 4.98%
1.21% 2.79% 0.10% 4.87%
0.72% 2.83% 0.21% 5.08%
0.22% 3.06% 0.25% 5.04%
state, and thus, being a dynamic system. In order to experiment the various level of machine utilization rate (due to the various level of demand rates), we scaled the production quantity of each month by various percentages: 80%, 85%, 90%, 95%, and 100%. Owing to the unavailability of the actual factory waiting time data, simulation was used to generate the required data for constructing the relationship curve in our experiment. The ¯ow times obtained by the simulation of various production schedules were compared with those predicted by the iterative approach. Tables 1 and 2 show the case of ®xed product mix and changing product mix, respectively. From these tables, we can see that the dierences between the ¯ow times obtained by simulation and the ¯ow times predicted by the iterative approach was very small with mean absolute dierence less than 5%. They were all within one standard deviation of the simulated ¯ow time. If assuming normal distribution and using 95% con®dence interval, whose boundary are plus and minus two standard deviation away from the mean value, the ¯ow times estimated by the iterative approach were well within the con®dence interval. The risk of estimation error for our iterative approach is very small. The result shows that our approach performs well in both steady state environment (®xed product mix) and dynamic environment (changing product mix); and its performance is as good as that of the queueing models with a steady state assumption [1,2], which use complex mathematical models to derive measures. We used a SUN UNIX work station to run our experiments. The execution time of our proposed iterative approach is less than 10 s/iteration. Thus, it took about 1.5 min to converge in nine iterations. The simulation programs were written in C++ [14]. It took about 10 min/ simulation run and about 5 h to run 30 simulation runs to obtain the average ¯ow times. The time of a simulation run is about 60 times that of an iteration run in our approach.
Table 2 A comparison of the simulation ¯ow time and the iterative ¯ow time in the case of changing product mix (%) Release scale
80%
85%
90%
95%
100%
M9 A9 D9 SD
2.42% 3.66% 0.16% 6.49%
2.02% 3.59% 0.25% 6.65%
1.55% 4.62% 0.33% 7.14%
1.86% 4.13% 0.68% 8.58%
0.19% 5.10% 0.89% 15.65%
Y.-F. Hung, C.-B. Chang / Computers & Industrial Engineering 37 (1999) 809±821
821
5. Conclusion and future direction Our approach bypasses the diculty of the mathematical modeling for various complex manufacturing methods. All analytical models aim to obtain the relationship between the machine utilization rate and the waiting time, which in this study are empirically obtained by direct construction from the historical data of machine utilization rate and waiting time. Based on the simulation experiment, our method works well in predicting the ¯ow times, which are important to the production planning and schedule for complex manufacturing system like semiconductor manufacturing. The industrial application of the iterative approach in realistic production plan and schedule problems requires further study. References [1] Chen H, Harrison JM, Mandelbaum A, Van Ackere A, Wein LM. Empirical evaluation of a queueing network model for semiconductor wafer fabrication. Operations Research 1988;36(2):202±15. [2] Connors DP, Feigin GE, Yao DD. A queueing newtwork model for semiconductor manufacturing. IEEE Transactions on Semiconductor Manufacturing 1996;9(3):412±27. [3] Gross D, Harris CM. Fundamentals of queueing theory, 2nd ed. USA: Wiley, 1985. [4] Hillier FS, Lieberman GJ. Introduction to operation research, 6th ed. New York: McGraw-Hill, 1995. [5] Hillier FS, Yu OS, Avis D, Fossett L, Lo F, Reiman M. Queueing tables and graphs. North-Holland, New York: Elsevier, 1981. [6] Hou M-C. Using queueing models to predict the cycle time of production planning model. Master Thesis, Department of Industrial Engineering, National Tsing Hua University, Hsinchu, Taiwan, 1995. [7] Hung Y-F. Corporate-level production planning with simulation feedback of parameters. Ph.D. Dissertation, College of Engineering, University of California, Berkeley, CA, 1991. [8] Hung Y-F, Leachman RC. A production planning methodology for semiconductor manufacturing based on iterative simulation and linear programming calculations. IEEE Transactions on Semiconductor Manufacturing 1996;9(2):257±69. [9] Johri PK. Practical issues in scheduling and dispatching in semiconductor wafer fabrication. Journal of Manufacturing Systems 1993;12(6):474±85. [10] Law AM, Kelton WD. Simulation modeling and analysis. New York: McGraw-Hill, 1991. [11] Leachman RC. Modeling technique for automated production planning in the semiconductor industry. In: Ciriani TA, Leachman RC, editors. Optimization in industry. England: Wiley, 1993. p. 1±30. [12] Leachman RC, Benson RF, Liu C, Raar DJ. IMPReSS: an automated production-planning and deliveryquotation system at Harris corporation Ð semiconductor sector. Interface 1996;26(1):6±37. [13] Lin C-Y, Shop ¯oor scheduling of semiconductor wafer fabrication using real-time feedback control and prediction. Ph.D. Dissertation, Engineering±Industrial Engineering and Operations Research, University of California at Berkeley, 1996. [14] Lippman SB. C++ primer, 2nd ed. MA, USA: Addison-Wesley, 1991. [15] Lu SCH, Ramaswamy D, Kumar PR. Ecient scheduling policies to reduce mean and variance of cycle-time in semiconductor manufacturing plants. IEEE Transactions on Semiconductor Manufacturing 1994;7(3):374±88. [16] Mathews JH. Numerical methods for mathematics, science, and engineering. Englewood Clis, NJ: PrenticeHall, 1992. [17] Uzsoy R, Lee CY, Martin-Vega LA. A review of production planning and scheduling models in the semiconductor industry. Part I: system characteristics, performance evaluation and production planning. IIE Transactions 1992;24(4):47±60. [18] Uzsoy R, Lee CY, Martin-Vega LA. A review of production planning and scheduling models in the semiconductor industry. Part II: shop-¯oor control. IIE Transactions 1994;26(5):44±55.