Data Mining for Cycle Time Key Factor Identification and Prediction in Semiconductor Manufacturing

Proceedings of the 13th IFAC Symposium on Information Control Problems in Manufacturing Moscow, Russia, June 3-5, 2009

Data Mining for Cycle Time Key Factor Identification and Prediction in Semiconductor Manufacturing

Y. Meidan, B. Lerner, M. Hassoun, G. Rabinowitz Department of Industrial Engineering and Management, Ben-Gurion University of the Negev, Beer Sheva, Israel (yairme, boaz ,[email protected])

Abstract: We suggest a data-driven methodology to identify key factors of the cycle time (CT) in a semiconductor manufacturing plant and to predict its value. We first extract a data set from a simulated fab and describe each operation in the set using 182 features (factors). Then, we apply conditional mutual information maximization for feature selection and the selective naïve Bayesian classifier for further selection and CT prediction. Prediction accuracy of 72.6% is achieved by employing no more than 20 features. Similar results are obtained by neural networks and the C5.0 decision tree. Keywords: industrial plant control, probabilistic models, machine learning, naïve Bayesian classifier.

1. INTRODUCTION Market conditions continually force the need to increase operational performance, e.g., through increased equipment utilization and productivity. To maintain a competitive edge in semiconductor manufacturing (SM), the cycle time (CT) – the cumulative time for completion of a work cycle (either a single step or a sequence of steps) – should be predicted and reduced (Chien et al., 2005; Backus et al., 2006). By identifying its key factors, the CT can be accurately predicted and minimized. Previous approaches to the optimization of production planning and control have mostly relied on rules of thumbs and heuristics (Chien et al., 2005) or on simulation, statistical analysis methods, analytical methods and hybrids of the previous three (Chung and Huang, 2002). Recently, another approach has been introduced into this field, known as machine learning and data mining (MLDM) (Chien et al., 2005; Backus et al, 2006). MLDM can be used to analyze complex problems represented by a large number of variables and high amounts of data (Bishop, 2007). In this study, MLDM is used to map the main CT factors and construct process models and estimations of CT. A first step towards CT reduction is the construction of a CT model which is descriptive to the role of the key factors affecting CT and quantitative with respect to CT prediction. For this purpose, we analyzed simulated SM data and used an MLDM approach – the naïve Bayesian classifier (NBC) – to transform raw data, typically collected and stored during the manufacturing process, into knowledge for supporting timely decisions on CT. The NBC was applied to CT explanation and prediction, providing also knowledge on CT key factors.

978-3-902661-43-2/09/$20.00 © 2009 IFAC

A broadly accepted reference model to designing MLDM, known as CRISP-DM (Cross-Industry Standard Process for Data Mining), has six primary phases: 1. Business Understanding: Getting familiar with business objectives and demands to convert them to MLDM problems. 2. Data Understanding: Describing data origins, background, nature and quality to gain familiarity with the data. 3. Data Preparation: Describing activities of converting raw data into a data set ready to be used by MLDM algorithms. 4. Modeling: Applying MLDM algorithms while calibrating their parameters to get optimal prediction models. 5. Evaluation: Assessing models' quality, appropriateness and validity to assure proper realization of business objectives. 6. Deployment: Organizing and presenting gained knowledge in a way that the customer can use and benefit from. The paper summarizes efforts taken to model and estimate CT according to the above phases for the benefit of manufacturing engineering in semiconductor fabs. 2. BUSINESS UNDERSTANDING Cycle time in SM is the time interval from bare silicon wafer start to final test (Akcali et al., 2001). It consists of queuing time for equipment, waiting time due to preventive or breakdown maintenance or engineering holds, and processing, inspection and transportation times. Within the harsh environment of this industry, characterized by high production, aggressive competition and fast changes, the CT of production lots is one of the most important performance measures. Reducing CT mean and variance is strategically important as it can meaningfully contribute to cost reduction, e.g., by decreasing the amount of work in progress (WIP). It also provides a competitive advantage by a faster time-tomarket, assisting in adapting to demand changes and in

217

10.3182/20090603-3-RU-2001.0466

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

quickly detecting yield decreases and their causes. Beyond reducing CT, predicting it is of tactical importance (Backus et al., 2006) as it has a practical viability for supporting production planning decisions (capacity requirement planning, resource scheduling, etc.), assessing customer due dates, controlling and monitoring CT, and improving factory decision factors such as inventory, unit cost and flexibility (Raddon and Grigsby, 1997). Cunningham et al. (1995) suggested grouping factors affecting CT using three categories: Process/product (wafer size, number of mask layers, etc.), facility (clean room class and size, facility age and region) and practice (SPC, CAM, yield modeling practices, etc.). Additional common factors previously analyzed are WIP levels (Chien et al., 2005), bottlenecks (Ramamurthi et al., 2005), MTBF and MTTR of production tools (Akcali et al., 2001), setups, wafer starts, lot size, priority & scheduling policies, and dispatch rules. Altogether, twenty nine distinct variables were reviewed from the literature as potential CT factors, a great deal of which is covered by the current research while introducing several additional ones. The literature survey revealed that most research of CTrelated performance measures have focused on the whole CT rather than any of its components. However, due to the assumption of rigid processing times on which the industrial engineer usually has no effect, the focus of the current research is set on the reduction of waiting time (WT) rather than CT. Identifying the top WT factors is a step towards eliminating waste from the manufacturing line and a means to narrow the range of possible actions for CT reduction. The set of key factors to be revealed from amongst all available factors throughout this research is ought to meet several basic requirements. First, they should be able to profoundly assist in distinguishing between different levels of operation WTs. Second, they ought to be as limited in number as possible, for both the benefits of model simplicity and generality. Third, they are desired to be feasibly under the influence of the manufacturing engineers, promoting model's usefulness. The basic business objectives of identifying and quantifying primary WT factors, designated to enable WT estimation and reduction, can be translated into MLDM objectives of feature selection and classification. Classification is a process which assigns a class label H(i) to an instance i described by a set of features. In our case, each instance is a feature (factor) representation of a fab operation and the classifier assigns one of a finite number of WT classes – each corresponds to a different range of WT – to the instance. The notion is that factors proving to accurately classify data into their correct WT level are the relevant WT key factors. Feature selection is the effort to extract the smallest set of key factors which contribute the most to classification accuracy.

3. DATA UNDERSTANDING AutoSched AP (Brooks Automation Inc.), which is a finite capacity planning and scheduling tool, was selected to model the working fab. SEMATECH Data set One, proposed by SEMATECH (a global, non-profit consortium for research in SM technology) was chosen for implementation. It represents a reduced scale fab of flash non-volatile memory in which tools, routing, products and flow characteristics are taken from actual fabs. Both the simulation model and software are well-known benchmarks in the semiconductor industry and the academic research (Johal, 1998; Ramamurthi et al., 2005). The SEMATECH One model is characterized by two high volume routes (technologies, products) varying from one another mainly by the path they follow in the manufacturing line. The number of operations (production steps) needed to complete routes one and two is 208 and 243 respectively, described with numerous characteristics of a true fab-level complexity. Over 120 simulation experiments were conducted, differentiated by three "tuning variables": Lot starts (time intervals between releasing new lots into the process) and distributions of MTBF & MTTR. Following the simulation runs, few dozens of features were collected or calculated, forming a vector representing each operation. A total of 451 vectors were extracted from each experiment, and summing to a total of 55,473 vectors for all experiments. To test the assumption that an operation's WT depends also on the characteristics of prior operations, every operation instance holds some information regarding its preceding operations. That is, an instance is composed of features of the operation (in time t) and the two previous operations (in times t-1 and t-2), together forming a segment. In this way, a total of 182 features were gathered for each instance as the raw data. These features may be divided into: simulationspecific identifiers (operation ID, tool ID, etc.), infrastructural variables available former to the simulation run (operation position in the process, tool type running the operation, batch restrictions, etc.), performance measures (indices of the line performance under certain scenarios, such as CT, WT, X-factor, WIP, etc.) and calculated variables (industry standards such as utilization over availability ratio, expected workload, etc.). Some of these features characterise an operation while others describe the tool on which the operation took place. Due to the data set's synthesized nature, the vast majority of feature instances were 100% complete. This data set represents a large variety of production scenarios, operations, routes, tool types, release rates, MTBFs and MTTRs, and it was validated with fab experts. Hence, it is considered to be representative. No errors were found in data; however some features (having alphanumeric or continuous values) needed pre-processing and transformation into nominal values, appropriate for analysis by the NBC. 4. DATA PREPARATION Preliminary data preparation included filtering experimentspecific identifiers and WT-related unimportant performance measures of current operation and segment (thus reducing the

218

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

number of candidate factors from 182 to 152), filling in missing values and reclassifying alphanumeric values. A special attention was devoted to a couple of pre-processing additional actions: Discretization of continuous features and an extensive feature selection. The drive for the former is to meet data type requirements induced by the NBC (which is more easily modelled using discrete, rather than continuous, data). The main motivation for feature selection is to improve accuracy (refraining from over-fitting the data) and reduce the high computational cost of the classifier. Discretization of continuous-valued features can be performed by employing numerous methods (Dougherty et al., 1995; Kohavi and Sahami, 1996). These can be categorized as supervised (utilizing class labels) vs. unsupervised, global (performed prior to induction) vs. local and static (determining the maximum number of intervals for each feature independently of the others) vs. dynamic. In this paper, we considered two of these methods to discretize WT factors and three for discretizing the WT variable. Concerning the predictors, the first method is simple equal frequency binning, which is unsupervised, global and static. It divides each continuous feature into k bins where each bin contains m/k instances of the total m instances in the data set. Since this method does not utilize class labels in setting bin boundaries, it is likely that classification information will be lost as a result of combining in the same bin instances that should be associated with different classes.

specify an interval out of all possible intervals. As oppose to the equal-frequency binning, here each bin is prone to be of a different size, depending on its level of entropy. The target mean WT originally ranges from 0 (no waiting time) to 176 hours. Weibull distribution was found to fit the data most accurately (Chi square test p<0.005, square error=0.000787), suggesting:

Mean _ WT = −0.001 + Weibull (4.84, 0.659).

(3)

Three discretization methods were considered for WT to set three WT levels, i.e., low ("1"), medium ("2") and high ("3"): a. Equal-frequency binning that maintains equal prior probabilities for the three levels. b. Manual selection of thresholds according to major visible changes in the WT distribution. c. L-Level Lloyd-Max quantization (Lloyd, 1982) that minimizes the mean square quantization error. Figure 1 depicts the distribution of the WT lower values (0 to 20 hours) and the thresholds induced by the above three methods, transforming WT from continuous to nominal.

- Equal Frequency - Manual Selection - Lloyd's

The second method applied to the continuous predictors, known as recursive minimal entropy partitioning (Fayyad and Irani, 1993), is global, static, and supervised. It uses the class information entropy

(1)

Fig. 1. Distribution and discretization thresholds of WT

to

select bins' boundaries for discretization, where P (Ci , S ) is the proportion of instances in data set S that

belong to class Ci , i = 1, k , and the logarithm expression measures the amount of information needed to specify Ci in S. Similarly to binary decision trees, this discretization method is based on the notion that any continuous feature A can be separated into two intervals – having S1 and S2 instances, respectively – by selecting a threshold TA. The optimal threshold is selected as the one which minimizes the weighted average class entropy

Feature selection. Preliminary experimentation proved an exponential growth of computation times for modeling. Redundancy and irrelevance of features were also apparent, showing that after selecting an optimal set of features to classify WT, not only that adding features did not contribute to the hit rate, but it actually offended it. This phenomenon indicated that when too many features were selected as predictors the classifier over fitted the training data and its generalization ability was negatively affected. Results of these runs had revealed that the number of features needed to reach the maximal hit rate was 47 at most. Therefore, it was decided to reduce the number of features from 152 down to 50 to decrease over fitting and to increase classifier induction speed.

(2)

of the resulting data sets S1 and S2 , having sizes |S1| and |S2|, respectively, among all candidate TA's. The method is then applied recursively to create multiple intervals on A. The stopping criterion is the minimal description length (MDL), defined as the minimum number of bits required to uniquely

Feature selection methods can be categorized as filters or wrappers (Kohavi and John, 1997). Filters are not dedicated to any specific type of classification method. Standard filters rank features, or sets of features, according to their predictive power, which can be estimated by means such as Fisher score, Kolmogorov-Smirnov test and Pearson correlation.

219

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

Wrappers rely on the performance of a classifier evaluating the quality of the feature or feature set. The main drawback of wrapper methods is their computation cost (time and memory). On the other hand, filters do not necessarily provide the required information to accurately discriminate between the classes or independence between features. In the current research, a recently-developed approach, which applies the criterion of conditional mutual information maximization (CMIM) (Fleuret, 2004) to feature selection, was implemented in MATLAB, with our own extension from two-class to multiclass classification problem. This filter method was followed by a wrapper method based on the selective NBC (SNBC) that is described in the Modeling section, as illustrated in Fig. 2.

other given the class variable (Langley et al., 1992). It consists of a network structure (Fig. 3) and a set of parameters that quantify it. Finding the structure is trivial due to the independence assumption. Variables X1, X2, …, Xm in the network structure are represented by nodes onto which all arcs from the class variable C are directed. In predicting WT, the observable and class variables represent the WT factors and WT variable, respectively.

Fig. 3. Network structure of the naive Bayesian classifier The NBC parameters are estimated by the class-conditional probability distributions P( x j | Ci ) for each variable xj given each state of the class variable Ci. All parameters are then set in a structure known as a conditional probability table, upon which any new instance I is classified using Bayes' theorem

Fig. 2. Design of feature selection stages The CMIM filter method iteratively picks features that maximize their mutual information with the class to predict conditionally to any feature already picked. Thereby, in contrast to filters proposed earlier, it guarantees a good tradeoff between independence and discrimination. Within a forward selection search for the top-N features, it makes certain that a feature similar to those already picked will not be selected, even if it is individually powerful. The motivation is that this feature does not carry extra information about the class to predict, other than the information already captured by features previously chosen. Fleuret (2004) reports of experiments in which CMIM and NBC (employed here) outperformed other combinations of feature selection and classification methods, while proving robustness when challenged by noisy data sets and achieving error rates similar or lower than state-of-the-art classifiers. 5. MODELING Previous CT studies in the semiconductor industry which employed MLDM techniques have relied mainly on decisiontrees and neural networks (Chien et al., 2005; Backus et al., 2006). We suggest here CT modeling and prediction using a version of NBC. NBC is probably the simplest approach to probabilistic classifier induction (Langley et al., 1992). Although not as accurate as a neural net or as intuitive as a decision tree, it has advantages in terms of simplicity, learning and classification speed, storage space and problem scalability (Domingos and Pazzani, 1997). Moreover, experimental studies of NBC revealed performance competitive with state-of-the-art classifiers (Domingos and Pazzani, 1997; Friedman et al., 1997). NBC is a Bayesian network that is restricted to variables independent of each

(4)

Due to the independence assumption for each class, the computation of the conditional probability P( ∧ x j | Ci ) can be simplified by Π j P( x j | Ci ) . Hence, the posterior probability P (Ci | I ) of each class is estimated based on the training set (5) The predicted class Ci of WT is selected as the class having the highest posterior probability P (Ci | I ) . The naive assumption of feature independence for each class almost never holds for natural data sets and can harm the classification process when violated. Therefore the selective NBC (SNBC) (Langley and Sage, 1994) was chosen, dealing with correlated features by incorporating only those features that when added to the classifier enlarge its accuracy. This forward selection induction scheme embeds the NBC in an algorithm that carries out a greedy search through the space of features, thus performing feature selection through the process of classifier induction, extracting the features most relevant to WT classification. Practically, selecting these features (i.e., best SNBC) is performed using a validation set that is independent of the training set. Each trained model (corresponding to a different feature sub-set) is evaluated using this validation set and the model achieving the highest classification accuracy is selected as the SNBC. The final

220

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

performance is then reported on a third set, independent of the previous two sets – the test set. Figure 4 demonstrates the SNBC classification accuracy using a validation set for increasing numbers of selected features. It can be seen that on average, 16 features enabled the SNBC to reach the maximal hit rate, which is about 72.6%. Not only that a model having too few features underfits the data but also a model having too many features overfits the data. In both cases, performance is inferior to that of the optimal model.

Fig. 4. Hit rate of SNBC (mean+/-std.) for increasing numbers of features

These results show that although more complex (i.e., having more features), WT models for a single product perform better than those for two products in terms of classification accuracy. This matches expectations as modeling the behaviour of a single product is prone to be less complicated than finding similarities between two non-identical products. For each of the two product combinations, the accuracy using the test set was always significantly (α=0.05) lower than using the validation set, especially for the model based on Lloyd's discretization of WT. The examination of the effect of non-WT feature discretization reveals that the models do not differ significantly from one another. This matches results of Dougherty et. al. (1995) and Hsu et al. (2003), whose conclusions, supported by empirical studies of a variety of discretization methods, suggested that entropy-based methods perform slightly better in some data sets, but on average, all discretization methods perform roughly the same. Out of the top-20 resultant WT factors, the one which was found most influential was the mean number of times the tool was in an idle state. The proportion of the three WT levels for each of the seven levels of this factor is portrayed in Fig. 5, where visually comparing levels 4&5 clearly reveals differences in the probability of the mean WT conditional to the factor's level. This was followed by the tool's mean availability, the number of times that the product had gone through the previous tool (number of loops), the number of operations shared by the tool upon which the operation was carried out, and the standard deviation of the length of queue.

Results for eight settings of classifiers are presented for the test set in Table 1. The first column of the table specifies the experimentation setting, structured as -. The second column reflects model complexity by showing the number of selected features. The third column denotes the maximal hit rate achieved on the validation sets and the forth shows the mean hit rate of models using the test sets. The figures specify the averages and standard deviations (in brackets) of 20 models from a 10 hold-out 2-fold cross-validation experiment design. Table 1. Complexity and accuracy of resultant models Experiment

Number of Selected Features

Maximal HitRate using Validation Set

Hit-Rate using Test Set

1-Eq.-Eq.

25.35 (13.83)

75.8 (0.7)

75.3 (0.4)

2-Eq.-Eq.

36.75 (7.14)

73.7 (0.4)

73.3 (0.2)

1,2-Eq.-Eq.

18.55 (9.06)

72.9 (0.5)

72.2 (0.4)

1,2-Eq.-Ent.

20.3 (7.82)

73.2 (0.55)

72.6 (0.5)

1,2-Man.-Eq.

29.7 (9.17)

71.6 (0.5)

70.7 (0.4)

1,2-Man.-Ent.

29 (6.97)

71.6 (0.6)

70.6 (0.4)

1,2-Lloyd's-Eq.

21.8 (6.01)

69.0 (0.3)

64.7 (0.6)

1,2-Lloyd's-Ent.

15.5 (7.71)

69.2 (0.5)

64.6 (0.8)

Fig. 5. WT (‘1’-low, ‘2’-medium, ‘3’-high) vs. number of idle states (‘1’-low, ..., ‘7’-high) SNBC's classification power was then assessed by comparison with the most common MLDM classifiers. For similar experimentation and performed on the same test set, a C5.0 decision tree achieved a mean hit rate of 74.0% (0.3) and a neural net accomplished 73.2% (0.74) using a similar number of features. This similarity in accuracies and numbers of selected features may imply to the problem complexity and the robustness of the resolved WT factors.

221

13th IFAC INCOM (INCOM'09) Moscow, Russia, June 3-5, 2009

6. EVALUATION In response to objectives, the results of the model in which WT was equal-frequency discretized whereas the predictors were discretized by the entropy-based method demonstrate how the SNBC can assist in classifying WT. From a base line of 33.33% prior probability of any WT level (low, medium or high), the generated models have proven to be capable of reaching an average hit rate of 72.6% using the test set, reflecting an improvement of about 40% in WT estimation accuracy. Concerning the previous operations, it was found that out of the top-20, 12 factors describe the current operation, 6 describe its predecessor (4 of which in the top10) and 2 of the former, both in the top-10. Apparently, previous operations do have an effect on the average WT of the current operation whose extent is higher than expected. 7. DEPLOYMENT As this research assumes a truthful representation of a fab by the simulated data, its resultant models and conclusions should be validated on genuine data, collected along with routine manufacturing. Following that, the models could be integrated into manufacturing lines, automatically collecting fab-wide data and processing them into WT estimations. As NBC is computationally inexpensive, it can even serve as a basis for real-time decision support, identifying potential delays in production and their causes, as well as better assessing due dates and fundamentally enhancing production planning. 8. CONCLUSIONS This paper aimed at identifying and quantifying the key WT factors. It showed how the uncertainty in WT estimation can dramatically decrease by using MLDM methods while employing no more than 20 key factors which interestingly include descriptors of past operations as well. Altogether, the current research lays practical and methodological grounds for future research, which can be enriched and extended for a study of other target functions (including WT, CT, X-factor, yield, etc.) and sets of respective features for analysis. The modular structure of experimentation, offering feature discretization and selection, can be effortlessly adopted by semiconductor manufacturers. Acknowledgment: This work was supported in part by the Paul Ivanier Center for Robotics and Production Management, Ben-Gurion University, Beer-Sheva, Israel. REFERENCES Akcali, E., Nemoto, K. and Uzsoy, R. (2001). Cycle time improvements for photolythography process in semiconductor manufaturing. IEEE Transactions on Semiconductor Manufacturing, 14, 48-56. Backus, P., Janakiram. M., Mowzoon, S., Runger G. C. and Bhargava, A. (2006). Factory cycle time prediction with a

machine learning approach. IEEE Transactions on Semiconductor Manufacturing, 19, 252-258. Bishop, C. M. (2007). Pattern Recognition and Machine Learning, Springer, New York, NY. Chien, C. F., Hsiao, C. W., Meng. C., Hong, K. T. and Wang, S. T. (2005). Cycle time prediction and control based on production line status and manufacturing machine learning. IEEE International Symposium on Semiconductor Manufacturing, ISSM 2005. Chung, S.H. and Huang, H.W. (2002). Cycle time estimation for wafer fab with engineering lots. IIE Transactions, 34, 105-118. Cunningham, P.S., Spanos C.J. and Voros K. (1995). Semiconductor yield improvement: results and best practices. IEEE Transactions on Semiconductor Manufacturing, 8, 103-109. Domingos, P., and Pazzani, M. (1997). On the optimality of the simple Bayesian classifier under zero-one loss. Machine Learning, 29, 103-130. Dougherty, J., Kohavi, R. and Sahami, M. (1995). Supervised and unsupervised discretization of continuous features. Proc. of the 12th International Conference on Machine Learning. (Prieditis A. and Russell S. (Ed)), 194–202. Morgan Kaufmann, Lake Tahoe, Los Altos, CA. Fayyad, U. M. and Irani, K. B. (1993). Multi-interval discretization of continuous-valued attributes for classification learning. Proc. of the 13th International Joint Conference on Artificial Intelligence, 13, 1022-1027. Fleuret, F. (2004). Fast binary feature selection with conditional mutual information. The Journal of Machine Learning Research, 5, 1531-1555. Friedman, N., Geiger, D. and Goldszmidt. (1997). Bayesian network classifiers. Machine Learning, 29, 131-163. Hsu, C.N., Huang, H.J. and Wong, T.T. (2003). Implications of the Dirichlet assumption for discretization of continuous variables in naïve Bayesian classifiers. Machine Learning, 53, 235–263. Johal, S. (1998). Application Report: Simulation reduces product cycle time. Semiconductor Intl., 21, 101-102. Kohavi, R. and John, G. (1997). Wrappers for feature subset selection. Artificial Intelligence, 12, 273–324. Kohavi, R., and Sahami, M. (1996). Error-based and entropybased discretization of continuous features. Proc. of the Second International Conference on Knowledge Discovery and Machine learning, 114–119. Langley, P., Iba, W. and Thompson, K. (1992). An analysis of Bayesian classifiers. Proc. of the Tenth National Conference on Artificial Intelligence, 223–228. Langley, P. and Sage, S. (1994). Induction of selective Bayesian classifiers. Proc. of the 10th Conf. on Uncertainty in Artificial Intelligence. Seattle, WA: Morgan Kaufman. Lloyd, S.P. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28, 129-137. Raddon, A. and Grigsby, B. (1997). Throughput time forecasting model. Advanced Semiconductor Manufturing Conf. and Workshop, 1997. IEEE /SEMI, 430-433. Ramamurthi, V., Kuhl, M.E. and Hirschman, M.E. (2005). Analysis of production control methods for semiconductor research and development fabs using simulation. Proc. of the 2005 Winter Simulation Conference, 2177-2185.

222