Precise and Accurate Job Cycle Time Forecasting in a Wafer Fabrication Factory with a Fuzzy Data Mining Approach

Abstract

Many data mining methods have been proposed to improve the precision and accuracy of job cycle time forecasts for wafer fabrication factories. This study presents a fuzzy data mining approach based on an innovative fuzzy backpropagation network (FBPN) that determines the lower and upper bounds of the job cycle time. Forecasting accuracy is also significantly improved by a combination of principal component analysis (PCA), fuzzy c-means (FCM), and FBPN. An applied case that uses data collected from a wafer fabrication factory illustrates this fuzzy data mining approach. For this applied case, the proposed methodology performs better than six existing data mining approaches.

1. Introduction

To forecast the job cycle time in a wafer fabrication factory is difficult because each wafer fabrication factory is a complicated production system with idiosyncratic features such as changing demand, a variety of product types and priorities, equipment unreliability, unbalanced capacity, job reentry into machines, alternative machines, sequence-dependent setup times, and shifting bottlenecks. An average job cycle time is several months with hundreds of hours of standard deviation. Many studies have shown that accurately predicting the cycle/completion times for such large systems is very difficult [2–4].

There are two approaches to forecasting the cycle time of a job. The first approach, the input-output relationship approach, is to determine certain factors (e.g., average waiting time, queuing length, utilization, future release plans, bottlenecks, etc.) that influence the job cycle time, and then to apply different approaches (e.g., multiple linear regression (MLR), artificial neural networks (ANN), etc.) to model the relationship between the job cycle time and these factors, in order to forecast the cycle time of a new job. The second approach, the time-series approach, is to treat fluctuations in the job cycle time as a type of time series. Many approaches, for example, moving average (MA), weighted moving average (WMA), exponential smoothing (ES), MLR, ANN, autoregressive integrated moving average (ARIMA), and others, can be applied to forecast the cycle time of a new job. The cycle time of a job is usually predicted when the job is released into the factory, but several months are required for the job to complete all operations. In such a case, the production controller may need to forecast the cycle time of job number 𝑘+10000 according to the historical data of jobs 1∼𝑘. This is a difficult problem for time-series methods; the present study contributes a useful approach to this problem.

Existing approaches to job cycle time estimation for a wafer fabrication factory can be classified into six categories: statistical analysis, production simulation (PS), ANN, case-based reasoning (CBR), fuzzy modeling methods, and hybrid approaches [5]. Many approaches are in fact special types of data mining methods. Table 1 summarizes references related to these categories. Data mining methods that have been extensively applied to job cycle time forecasting include clustering, genetic algorithms (GA), ANNs, decision trees, and hybrid approaches. The relevant references are reviewed below.

Table 1 

References relevant to the six categories.

Backus et al. [1] compared the performances of clustering, 𝐾-nearest neighbors (i.e., CBR), classification and regression tree (CART), and ANN; CART gave the best performance of that experiment. Meidan et al. [6] focused on the choice of key factors in cycle time forecasting. They compared multiple linear regression (i.e., statistical analysis), ANN, and a selective naive Bayesian classifier (SNBC); they found that their ANN had the best performance.

Wang and Mendel [7] proposed a fuzzy modeling method called the WM method; the first step of that method was partitioning the range of each input variable into several fuzzy intervals. Chang et al. [8] modified the first step of the WM method with a simple GA and proposed an evolving fuzzy rule (EFR) approach to predict the cycle time of a job in a wafer fabrication factory. Their EFR approach outperformed CBR and ANN in terms of forecasting accuracy.

A BPN is an effective tool for modeling complex physical systems described by sets of different equations. BPNs are useful for prediction, control, and design purposes; they can offer both effectiveness (high forecasting accuracy) and efficiency (short execution time). Chang and Hsieh [9], Chang et al. [8], and Sha and Hsu [10] all predicted the cycle times of jobs in wafer fabrication factories with backpropagation networks (BPNs) that each had a single hidden layer. These BPNs delivered better average forecasting accuracy than statistical analysis approaches, as measured by root mean squared error (RMSE). For example, an improvement of about 40% in RMSE was achieved in Chang et al. [8]. Chen [5] constructed a fuzzy backpropagation network (FBPN) that incorporated expert opinions to modify the inputs of the FBPN. Chen’s FBPN surpassed a crisp BPN especially with respect to efficiency. Chen [11] incorporated the job release plan of a wafer fabrication factory into a BPN, and constructed a “look-ahead” BPN for the same purpose, which led to an average reduction of 12% in RMSE. Beeg [12] and Chen et al. [13, 14] predicted the cycle time of a job in a ramping-up wafer fabrication factory. Chen et al. used a BPN-based method; Beeg tried to figure out the impact of utilization for the cycle time. Tirkel [15] applied two data mining approaches—decision trees and ANN—and concluded that ANN was more accurate.

Chang and Liao [26] combined a self-organizing map and the Wang and Mendel method (SOM-WM); a job was classified using a SOM before the job cycle time was predicted by the WM method. Chen [11] constructed a look-ahead k-means-fuzzy backpropagation network (kM-FBPN) for the same purpose and discussed the effects of using different look-ahead functions. More recently, Chen et al. [23] proposed a look-ahead SOM-FBPN approach for job cycle time forecasting in a semiconductor factory; a set of fuzzy inference rules was developed to evaluate the achievability of a cycle time forecast. Subsequently, Chen et al. [3] added a selective time allowance to the look-ahead SOM-FBPN approach to determine internal due dates. Further, Chen [24] showed that a combined SOM and FBPN could be improved by feeding back the forecasting error by the FBPN to adjust the classification results of the SOM. Chen et al. [13, 14] proposed a postclassification FBPN-BPN approach in which a job was not preclassified; the job was classified after the cycle time had been predicted. Experimental results showed that the post-classification approach was better than some preclassification approaches in some cases. In order to combine the advantages of preclassifying and post-classifying approaches, Chen [25] proposed a bidirectional classifying approach that combined FCM, FBPN, and a radial basis function network (RBF), in which jobs are not only preclassified but also postclassified. Instead of bi-directional classification, Chien et al. [4] used nonlinear regression equations and then related the forecasting error to some factory conditions and job attributes with a backpropagation network (BPN) to improve the forecasting accuracy.

Table 2 summarizes the data mining methods in this field. The advantages and disadvantages of these methods are compared in Table 3.

Table 2 

Data mining methods for job cycle time forecasting.

Table 3 

Advantages and disadvantages of various data mining approaches for job cycle time forecasting.

To improve the precision and accuracy of job cycle time forecasting in a wafer fabrication factory, this study presents a fuzzy data mining approach that emphasizes both accuracy and precision.(1) Accuracy: the forecast value should be as close as possible to the actual value.(2) Precision: a narrow interval containing the actual value should be established.Chen [30] defined several methods to measure the precision and accuracy of a fuzzy forecasting method and noted that precision has seldom been considered. However, past studies have confirmed that precision and accuracy are closely related. Song and Chissom [31] proposed a first-order time-variant model to forecast enrollment and considered the effect of the model basis on the forecasting precision. Specifically, if the space of feasible solutions is narrowed, it becomes more likely that the actual value can be found [30, 32].

The fuzzy data mining approach has the following innovative characteristics(1) Some factors used to estimate the cycle time are dependent on each other, which may cause problems in classifying jobs and in fitting the relationship between the job cycle time and these factors. To solve this problem, PCA, FCM, and FBPN are combined to predict the job cycle time, which is an innovation in this field. (2) The job cycle time is predicted by an effective FBPN approach. For the job cycle time, there are upper and lower bounds that constrain the forecast to be within a possible range. The FBPN approach from Chen and Wang’s FBPN research [33] has been modified for this purpose. In Chen and Wang’s FBPN approach two nonlinear programming (NP) models are solved to determine the upper and lower bounds of the job cycle time. However, the NP models involve complicated constraints and therefore are difficult to solve. In addition, the NP models will become too huge if many jobs are to be considered. To solve these problems, in this study, the upper and lower bounds are determined in a simpler way.The differences between the proposed methodology and various previous methods are summarized in Table 4.

Table 4 

The differences between the proposed methodology and the previous methods.

The remainder of this paper is organized as follows. Section 2 introduces the proposed methodology, which is composed of six steps. Section 3 demonstrates the proposed method with an illustrative example, then presents another case with data collected from a real wafer fabrication factory, and makes comparisons with some existing approaches. Finally, concluding remarks and directions for future research are given in Section 4.

2. Methodology

The fuzzy data mining procedure consists of several steps.

Step 1. Replacing parameters using PCA.

Step 2. Classifying jobs using fuzzy c-means (FCM).

Step 3. Forecasting the cycle times of jobs in each category using a FBPN.

Step 4. Updating the upper and lower bounds.

Step 5. Evaluating and comparing the forecasting precision and accuracy of the proposed methodology.

Step 6. If the improvement in the forecasting precision or accuracy becomes negligible, stop; otherwise, return to Step 3.

A flow chart of the proposed methodology is shown in Figure 1. These steps can be compared to the basic steps of data mining, as shown in Table 5.

Table 5 

The proposed methodology compared to data mining.

Figure 1 

The flowchart of the proposed methodology.

2.1. Step 1: Variable Replacement Using PCA

First, PCA is used to replace the inputs to the FBPN. A series of linear combinations of the original variables forms a new variable so that these new variables are as unrelated to each other as possible. Although there are more advanced applications of PCA, in this study PCA is used to enhance the efficiency of FBPN training. PCA consists of the four following steps.(1) Raw data standardization: to eliminate dimensional differences and large numerical differences the original variables are standardized [13, 14]: 𝑥^∗_𝑗𝑖=𝑥_𝑗𝑖−𝑥_𝑖𝜎_𝑖,𝑥_𝑖=∑^𝑛_𝑗=1𝑥_𝑗𝑖𝑛,𝜎_𝑖=∑^𝑛_𝑗=1𝑥_𝑗𝑖−𝑥_𝑖²,𝑛−1 (1) where 𝑥_𝑗𝑖 is the 𝑖th attribute of job 𝑗, 𝑗=1∼𝑛; 𝑥_𝑖 and 𝜎_𝑖 indicate the mean and standard deviation of variable 𝑖, respectively.(2) Establishment of the correlation matrix 𝑅 is as follows: 1𝑅=𝑋𝑛−1^∗𝑇𝑋^∗, (2) where 𝑋^∗ is the standardized data matrix. The eigenvalues and eigenvectors of 𝑅 are calculated and represented as 𝜆₁∼𝜆_𝑚 and 𝑢₁∼𝑢_𝑚, respectively; 𝜆₁≥𝜆₂≥⋯≥𝜆_𝑚.(3) Determination of the number of principal components: the variance contribution rate is calculated as 𝜂_𝑞=𝜆_𝑞∑^𝑚_𝑟=1𝜆_𝑟⋅100%, (3) and the accumulated variance contribution rate is 𝜂_Σ𝑝=^𝑝_𝑞=1𝜂_𝑞. (4) Choose the smallest 𝑃 value such that 𝜂_Σ(𝑝) rangesfrom 85% to 90%. A Pareto analysis chart can be used to show the percentage of variability that can be attributed to each principal component.(4) Formation of the following matrices: 𝑈_𝑚×𝑝=𝑢₁,𝑢₂,…,𝑢_𝑝,𝑍_𝑛×𝑝=𝑋^∗_𝑛×𝑚𝑈_𝑚×𝑝. (5)𝑍_𝑛×𝑝=[𝑧_𝑗𝑞]  (𝑗=1∼𝑛;𝑞=1∼𝑝) are the component scores. These scores contain the coordinates of the original data in the new coordinate system defined by the principal components. They will be used as the new inputs to the FBPN.

2.2. Step 2: Classifying Jobs Using FCM

After PCA, examples are classified by FCM. If a crisp clustering method is applied instead, then it is possible that some clusters will have very few examples. FCM avoids the problems of crisp clustering, because in FCM, an example belongs to multiple clusters to different degrees. Similarly, in probability theory the naïve Bayes method provides a probability that each item belongs to a class. However, the classification of jobs in FCM can include subjective judgments.

FCM classifies jobs by minimizing the following objective function: Min^𝐾𝑛_𝑘=1𝑗=1𝜇^𝑚_𝑗(𝑘)𝑒²_𝑗(𝑘), (6) where 𝐾 is the required number of categories; 𝑛 is the number of jobs; 𝜇_𝑗(𝑘) indicates the degree of membership that job 𝑗 has in category 𝑘; 𝑒_𝑗(𝑘) measures the distance from job 𝑗 to the centroid of category 𝑘;  𝑚∈[1,∞) is a parameter to adjust the fuzziness and is usually set to 2. The procedure of FCM is as follows.(1) Normalize the input data.(2) Produce a preliminary clustering result.(3) Iterations calculate the centroid of each category as 𝑧_(𝑘)=𝑧_(𝑘)𝑞,𝑘=1∼𝐾,𝑧_(𝑘)𝑞=∑^𝑛_𝑗=1𝜇^𝑚_𝑗(𝑘)𝑧_𝑗𝑞∑^𝑛_𝑗=1𝜇^𝑚_𝑗(𝑘)𝜇,𝑘=1∼𝐾,𝑞=1∼𝑝,_𝑗(𝑘)=1∑^𝐾_𝑔=1𝑒_𝑗(𝑘)/𝑒_𝑗(𝑔)^2/(𝑚−1)𝑒,𝑗=1∼𝑛,𝑘=1∼𝐾,_𝑗(𝑘)=^𝑝_𝑞=1𝑧_𝑗𝑝−𝑧_(𝑘)𝑝²,𝑗=1∼𝑛,𝑘=1∼𝐾, (7) where 𝑧_(𝑘) is the centroid of category 𝑘. 𝜇^(𝑡)_𝑗(𝑘) is the degree of membership that job 𝑗 has in category 𝑘 after the 𝑡th iteration.(4) Remeasure the distance from each job to the centroid of each category and then recalculate the corresponding membership.(5) Stop if the following condition is met. Otherwise, return to Step (3): max_𝑘max_𝑗𝜇_𝑗(𝑘)^(𝑡)−𝜇_𝑗(𝑘)^(𝑡−1)<𝑑, (8) where 𝑑 is a real number representing the membership convergence threshold.The performance of FCM is affected by the settings of the initial values; initialization can be repeated multiple times in order to find the optimal solution. Finally, the separate distance test (𝑆 test) proposed by Xie and Beni [34] can be applied to determine the optimal number of categories 𝐾: Min𝑆subjectto𝐽_𝑚=^𝐾𝑛_𝑘=1𝑗=1𝜇^𝑚_𝑗(𝑘)𝑒²_𝑗(𝑘),𝑒²_min=min_{𝑘1≠𝑘2}^𝑝_𝑞=1𝑧_(𝑘1)𝑞−𝑧_(𝑘2)𝑞²,𝐽𝑆=^𝑚𝑛×𝑒²_min,𝐾∈𝑍⁺. (9) The 𝐾 value that minimizes 𝑆 determines the optimal number of categories.

2.3. Step 3: Forecasting the Cycle Times of Jobs in Each Category Using a FBPN

All categories in the present study are considered by the FBPN. There are different types of FBPNs, for example, Nomura et al. [35], Chen [5], Lin [36], and Chen and Wang [33]. Among them, Chen and Wang’s FBPN is unique since it is aimed at optimizing precision. For this reason, it is chosen for this study. This FBPN is a proven technology that has incorporated many features of prior artificial neural networks, which have solved a wide variety of problems characterized by sets of different equations. Some special problems can benefit from specialized tools, such as compositional pattern-producing networks, cascading neural networks, and dynamic neural networks, but for the present job cycle time problem, a well-trained FBPN with an optimized structure can still produce very good results. All parameters in the FBPN are expressed in triangular fuzzy numbers (TFNs). There are various types of fuzzy numbers with different shapes; among them, the triangular type is easily implemented and has been universally applied to numerous applications (e.g., [37, 38]). Although the proposed methodology uses TFNs, it can be easily modified to use other types of fuzzy numbers, such as trapezoidal fuzzy numbers or bounded LR-type fuzzy numbers.

The FBPN is configured as follows.(1) Inputs: there are 𝑗 sets of inputs; each set consists of the new factors determined by PCA from the 𝑗th job. These factors have to be partially normalized so that their values fall within [0.1, 0.9] [13, 14]. (2) Single hidden layer: generally one hidden layer gives the best convergence results for the FBPN.(3) For simplicity, the number of neurons in the hidden layer is twice the number in the input layer. An increase in the number of hidden-layer nodes lessens the output errors for the training examples, but increases the errors for novel examples. This phenomenon is often called “overfitting.” Some research has considered the connections between the complexity of a FBPN, the performance for the training data, and the number of examples; noteworthy research includes Akaike’s information criterion (AIC) [39] and Rissanen’s minimum description length (MDL) [40]. (4) Output: the (normalized) cycle time forecast of the example.

The procedure for determining the parameter values is now described. After pre-classification, a portion of the adopted examples in each category is fed as “training examples” into the FBPN to determine the parameter values. Two phases are involved at the training stage. First, in the forward phase, inputs are multiplied with weights, summed, and transferred to the hidden layer. Then, activated signals are outputted from the hidden layer as ℎ_𝑙=ℎ_𝑙1,ℎ_𝑙2,ℎ_𝑙3=11+𝑒⁻𝑛^ℎ_𝑙=11+𝑒^{−𝑛ℎ𝑙1},11+𝑒^{−𝑛ℎ𝑙2},11+𝑒^{−𝑛ℎ𝑙3}, (10) where 𝑛^ℎ_𝑙=𝑛^ℎ_𝑙1,𝑛^ℎ_𝑙2,𝑛^ℎ_𝑙3=𝐼^ℎ_𝑙̃(−)𝜃^ℎ_𝑙=𝐼^ℎ_𝑙1−𝜃^ℎ_𝑙3,𝐼^ℎ_𝑙2−𝜃^ℎ_𝑙2,𝐼^ℎ_𝑙3−𝜃^ℎ_𝑙1,𝐼^ℎ_𝑙̃=𝐼^ℎ_𝑙1,𝐼^ℎ_𝑙2,𝐼^ℎ_𝑙3=_all𝑘𝑤^ℎ_𝑘𝑙⋅𝑥_𝑘=_all𝑘min𝑤^ℎ_𝑘𝑙1𝑥_𝑘,𝑤^ℎ_𝑘𝑙3𝑥_𝑘,_all𝑘𝑤^ℎ_𝑥𝑘2𝑥_𝑘,_all𝑘max𝑤^ℎ_𝑘𝑙1𝑥_𝑘,𝑤^ℎ_𝑘𝑙3𝑥_𝑘 (11)ℎ_𝑙 is the output from hidden-layer node 𝑙,𝑙=1∼𝐿;  𝜃^ℎ_𝑙 is the threshold for screening out weak signals by hidden-layer node 𝑙; 𝑤^ℎ_𝑘𝑙 is the weight of the connection between input node 𝑘 and hidden-layer node 𝑙,𝑘=1∼𝑝;𝑙=1∼𝐿;  𝑥_𝑖𝑗 is the 𝑘th input, 𝑘=1∼𝑝. The remaining parameters are transition variables. (−) and (×) denote fuzzy subtraction and multiplication, respectively.

ℎ_𝑙s are also transferred to the output layer with the same procedure. Finally, the output of the FBPN is generated as follows: 𝑜_𝑗̃=𝑜_𝑗1,𝑜_𝑗2,𝑜_𝑗3=11+𝑒⁻𝑛^𝑜=11+𝑒^−𝑛𝑜1,11+𝑒^−𝑛𝑜2,11+𝑒^−𝑛𝑜3, (12) where 𝑛^𝑜=𝑛^𝑜₁,𝑛^𝑜₂,𝑛^𝑜₃=𝐼^𝑜̃(−)𝜃^𝑜=𝐼^𝑜₁−𝜃^𝑜₃,𝐼^𝑜₂−𝜃^𝑜₂,𝐼^𝑜₃−𝜃^𝑜₁,𝐼^𝑜̃=𝐼^𝑜₁,𝐼^𝑜₂,𝐼^𝑜₃=_all𝑙𝑤^𝑜_𝑙(×)ℎ_𝑙≅_all𝑙min𝑤^𝑜_𝑙1ℎ_𝑙1,𝑤^𝑜_𝑙3ℎ_𝑙3,_all𝑙𝑤^𝑜_𝑙2ℎ_𝑙2,_all𝑙max𝑤^𝑜_𝑙1ℎ_𝑙1,𝑤^𝑜_𝑙3ℎ_𝑙3. (13)(14)𝑜_𝑗̃ is the network output, which is the normalized value of the cycle time forecast of job 𝑗; 𝜃^𝑜 is the threshold for screening out weak signals by the output node; 𝑤^𝑜_𝑙 is the weight of the connection between hidden-layer node 𝑙 and the output node;  𝑙=1∼𝐿.

Subsequently, in the backward phase, the training of the FBPN is decomposed into three subtasks: determining the center value and upper and lower bounds of the parameters.

First, to determine the center of each parameter (such as 𝑤^ℎ_𝑘𝑙2, 𝜃^ℎ_𝑙2, 𝑤^𝑜_𝑙2, and 𝜃^𝑜₂), the FBPN is treated as a crisp network. Algorithms applicable for this purpose include the gradient descent algorithms, the conjugate gradient algorithms, the Levenberg-Marquardt algorithm, and others. In this study, the Levenberg-Marquardt algorithm is applied. The Levenberg-Marquardt algorithm was designed for training with second-order speed without having to compute the Hessian matrix. It uses approximation and updates the network parameters in a Newton-like way, as described below.

The network parameters are placed in vector 𝜷=[𝑤^ℎ₁₁,…,𝑤^ℎ_𝑝𝐿,𝜃^ℎ₁,…,𝜃^ℎ_𝐿,𝑤^𝑜₁,…,𝑤^𝑜_𝐿,𝜃^𝑜]. The network output 𝑜_𝑗 can be represented with 𝑓(𝐱_𝑗,𝜷). The objective function of the FBPN is to minimize RMSE or equivalently the sum of squared error (SSE): SSE(𝜷)=^𝑛_𝑗=1(𝑁CT_𝑗−𝑓𝐱_𝑗,𝜷)². (15)

The Levenberg-Marquardt algorithm is an iterative procedure. In the beginning, the user should specify the initial values of the network parameter vector 𝜷. It is a common practice to set 𝜷^T=(1,1,…,1). In each step, the parameter vector 𝜷 is replaced by a new estimate 𝜷+𝜹, where 𝜹=[Δ𝑤^ℎ₁₁,…,Δ𝑤^ℎ_𝑝𝐿,Δ𝜃^ℎ₁,…,Δ𝜃^ℎ_𝐿,Δ𝑤^𝑜₁,…,Δ𝑤^𝑜_𝐿,Δ𝜃^𝑜]. The network output becomes 𝑓(𝐱_𝑗,𝜷+𝜹); it is approximated by its linearization as 𝑓𝐱_𝑗,𝜷+𝜹≈𝑓𝐱_𝑗,𝜷+𝐉_𝑗𝜹, (16) where 𝐉_𝑗=𝜕𝑓𝐱_𝑗,𝜷𝜕𝜷 (17) is the gradient vector of 𝑓 with respect to 𝜷. Substituting (16) into (15) gives SSE(𝜷+𝜹)≈^𝑛_𝑗=1𝑁CT_𝑗−𝑓𝐱_𝑗,𝜷−𝐉_𝑗𝜹². (18) When the network reaches the optimal solution, the gradient of SSE with respect to 𝜹 will be zero. Taking the derivative of SSE(𝜷+𝜹) with respect to 𝜹 and setting the result to zero gives 𝐉^T𝐉𝜹=𝐉^T𝑁CT_𝑗−𝑓𝐱_𝑗,𝜷, (19) where 𝐉 is the Jacobian matrix containing the first derivative of network error with respect to the weights and biases. Equation (19) includes a set of linear equations that can be solved for 𝜹.

2.4. Step 4: Determining the Upper Bounds of Parameters in the FBPN

In Chen and Wang [33] and Chen and Lin [41], an NP model is constructed to adjust the connection weights and thresholds in the FBPN; such an NP model is not easy to solve. The proposed methodology only makes iterative adjustments to the threshold on the output node. This way is much simpler and can also achieve good results.

Substituting (13) into (12) gives 𝑜_𝑗2=11+𝑒^{−𝑛𝑜𝑗2}=11+𝑒^{−(𝐼𝑜𝑗2−𝜃𝑜2)}=11+𝑒^{𝜃𝑜2−𝐼𝑜𝑗2}. (20) Therefore, 1ln𝑜_𝑗2−1=𝜃^𝑜₂−𝐼^𝑜_𝑗2. (21) So 𝐼^𝑜_𝑗2=𝜃^𝑜₂1−ln𝑜_𝑗2−1. (22)

Assume that the adjustment made to the output node threshold is denoted as Δ𝜃^𝑜=𝜃^𝑜₃−𝜃^𝑜₂. After adjustment, the output from the new FBPN, 𝑜_𝑗3, determines the upper bound of the cycle time: 𝑜_𝑗3=11+𝑒^{−𝑛𝑜𝑗3}, (23) where 𝑛^𝑜_𝑗3=𝐼^𝑜_𝑗3−𝜃^𝑜₃=𝐼^𝑜_𝑗3−𝜃^𝑜₂+Δ𝜃^𝑜. (24) Substituting (24) into (23) gives 𝑜_𝑗3=11+𝑒^{−(𝐼𝑜𝑗3−𝜃𝑜2−Δ𝜃𝑜)}, (25) and substituting (22) into (25) gives 𝑜_𝑗3=11+𝑒^{−(𝜃𝑜2−ln(1/𝑜𝑗2−1)−𝜃𝑜2−Δ𝜃𝑜)}=11+𝑒^{ln(1/𝑜𝑗2−1)+Δ𝜃𝑜}=11+𝑒^Δ𝜃𝑜1/𝑜_𝑗2.−1 (26) Obviously, the maximum of Δ𝜃^𝑜 determines the lowest upper bound, since 𝑜_𝑗3 is the upper bound of the cycle time, 𝑜_𝑗3≥𝑁(CT_𝑗), 11+𝑒^{ln(1/𝑜𝑗2−1)+Δ𝜃𝑜}≥𝑁CT_𝑗,Δ𝜃^𝑜1≤ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1. (27)(28) Equation (28) holds for all jobs, so Δ𝜃^𝑜≤min_𝑗1ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1. (29) According to (26), the optimal value of Δ𝜃^𝑜 should be set to the maximum possible value: Δ𝜃^𝑜∗=min_𝑗1ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1. (30)

Then the FBPN’s optimization results depend on the initial conditions and therefore are different at every iteration. Assume that the optimal value of 𝑜_𝑗3 in the 𝑡th iteration is denoted by 𝑜_𝑗3(𝑡), then, after some iterations, 𝑜_𝑗3(alliterations)=min_𝑡𝑜_𝑗3(𝑡). (31) In this way, the upper bound of the cycle time is decreased gradually (see Figure 2). Another merit of this approach is that it does not rely on the parameters of the FBPN.

Figure 2 

Iterative reduction of the upper bound.

Theorem 1. The upper bound can be found by considering only jobs with 𝑁(CT_𝑗)≥𝑜_𝑗2.

Proof. Divide jobs into two categories 𝐶₁ (jobs with 𝑁(CT_𝑗)<𝑜_𝑗2) and 𝐶₂ (jobs with 𝑁(CT_𝑗)≥𝑜_𝑗2). Equation (30) becomes Δ𝜃^𝑜∗=minmin_𝑗∈𝐶11ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1,min_𝑗∈𝐶21ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1. (32) If 𝑁(CT_𝑗)<𝑜_𝑗2, then 1ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1>0. (33) Therefore, min_𝑗∈𝐶11ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1>0. (34) In contrast, if 𝑁(CT_𝑗)≥𝑜_𝑗2, then 1ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1≤0,min_𝑗∈𝐶21ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1≤0. (35) Finally, Δ𝜃^𝑜∗=min_𝑗∈𝐶21ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1. (36) Theorem 1 is proved.

2.5. Step 4: Determining the Lower Bounds of Parameters in the FBPN

In a similar way, the threshold on the output node can be modified to determine the lower bound of the cycle time forecast so that the actual value will be greater than the lower bound of the network output. The optimal value of Δ𝜃^𝑜 can be obtained as Δ𝜃^𝑜∗=max_𝑗1ln𝑁CT_𝑗1−1−ln𝑜_𝑗2−1. (37)

Theorem 2. The lower bound can be found by considering only jobs with 𝑁(CT_𝑗)<𝑜_𝑗2.

The proof of Theorem 2 resembles that of Theorem 1.

Proof. Assume that the optimal value of 𝑜_𝑗1 in the 𝑡th replication is indicated with 𝑜_𝑗1(𝑡), then after some iterations, 𝑜_𝑗1(alliterations)=max_𝑡𝑜_𝑗1(𝑡). (38) In this way, the upper bound of the cycle time is increased gradually. Δ𝜃^𝑜∗does not rely on the parameters of the FBPN.

3. Application and Analyses

To demonstrate the application of the proposed methodology, an illustrative example containing the data of 40 jobs (see Table 6) was used.

Table 6 

The illustrative example.

We standardize the data and obtain the correlation matrix 𝑅. The eigenvalues and eigenvectors of 𝑅 are calculated as 𝜆₁=2.66,𝜆₂=1.15,𝜆₃𝜆=0.94,₄=0.83,𝜆₅=0.25,𝜆₆𝑢=0.02,₁=0.200.520.160.560.070.59,𝑢₂=,𝑢0.47−0.190.68−0.19−0.490.05₃=−0.270.31−0.320.03−0.85−0.02,𝑢₄=,𝑢0.81−0.02−0.570.06−0.04−0.15₅=−0.12−0.72−0.070.65−0.180.09,𝑢₆=,0.020.290.280.470.03−0.78 (39)(40) respectively. The variance contribution rates are 𝜂₁=46%;𝜂₂=20%;𝜂₃𝜂=16%;₄=14%;𝜂₅=4%;𝜂₆=0%. (41) Summing up 𝜂_𝑞’s, we obtain 𝜂_Σ(1)=46%;𝜂_Σ(2)=65%;𝜂_Σ(3)𝜂=81%;_Σ(4)=95%;𝜂_Σ(5)=100%;𝜂_Σ(6)=100%. (42) A Pareto analysis chart is used to compare the percentage of variability explained by each principal component (see Figure 3). There is a clear break in the amount of variance accounted for by each component between the first and second components. However, that component by itself explains less than 50% of the variance, so more components are considered. To meet the requirement 𝜂_Σ(𝑝)≥85%∼90%, 𝑝 is chosen as 3. We can see that the first three principal components explain roughly 80% of the total variability in the standardized data, so we reduce the issue to these three components in order to visualize the data.

Figure 3 

The Pareto analysis chart.

The component scores are recalculated as the new inputs to the FCM-FBPN (see Table 7). These scores contain the coordinates of the original data in the new coordinate system defined by the principal components.

Table 7 

New inputs to the FCM-FBPN.

Subsequently, jobs are classified using FCM based on the new variables. The results of the 𝑆 test are summarized in Table 8. In this case, the optimal number of job categories is 5. However, there will be some categories with very few jobs. For this reason, the second best solution is used, that is, 4 categories. A common practice is to set a threshold of membership 𝜇_𝐿 to determine whether a job belongs to each category. With the decrease in the threshold, each category contains more jobs. Because the overall forecasting performance is affected by the outliers, 𝜇_𝐿 should not be set too high. In this case, 𝜇_𝐿 is set to 0.3. The classifying results are shown in Table 9.

Table 8 

The results of the 𝑆 test.

Table 9 

The classifying results (𝜇_𝐿=0.3).

After job classification, 3/4 of the examples in each category are used as the training example. The remaining 1/4 is left for testing. A three-layer FBPN is then used to predict the cycle time of jobs in each category according to the new variables with the following settings. Single hidden layer The number of neurons in the hidden layer: 2∗3=6 Convergence criterion: SSE<10⁻⁶ or 10000 epochs have been run Iterations: 5

The forecasting results are shown in Figure 4. The performance level of the proposed methodology is compared with those of statistical analysis (i.e., multiple linear regression), CART, PCA-CART, BPN, FCM-BPN, and PCA-BPN in Table 10. The forecasting precision is evaluated with the average range: ∑averagerange=^𝑛𝑗=1𝑜_𝑗3−𝑜_𝑗1𝑛. (43) The average range of a nonbiased crisp approach can be measured with 6𝜎(or 6 RMSE): 𝜎²=RMSE²𝑛,𝑛−𝑘−1 (44) where RMSE=∑^𝑛_𝑗=1𝑎_𝑗−𝑓_𝑗²𝑛. (45)𝑘 is the number of independent variables; 𝑎_𝑗 and 𝑓_𝑗 denote the actual value and forecast of job 𝑗, respectively;  𝑛 is the total number of data. When 𝑛→∞, 𝜎²≈RMSE². Theoretically, the probability that such an interval contains the actual value is about 99.7%, under the assumption that residuals follow a normal distribution. The forecasting accuracy is measured with MAE, MAPE, and RMSE: ∑MAE=^𝑛𝑗=1𝑎_𝑗−𝑓_𝑗𝑛,∑MAPE=^𝑛𝑗=1𝑎_𝑗−𝑓_𝑗/𝑎_𝑗𝑛. (46) The CART algorithm is based on that of Breiman et al. [42]. The results of CART are shown in Figure 5. In PCA-CART, the new variables found by PCA are used as input variables to determine the job cycle time. The results of PCA-CART are shown in Figure 6. The adjustment of the lower and upper bounds for the proposed methodology is shown in Figure 7.

Table 10 

Comparison of the forecasting performance levels.

Figure 4 

The forecasting results from the FBPN.

Figure 5 

Results of CART.

Figure 6 

Results of PCA-CART.

(a)

(b)

(c)

Figure 7 

The adjustment of the lower bounds.

According to the experimental results, we have the following.(1) The nonlinear nature of this problem is obvious, since the performance of statistical analysis is poor, and statistical analysis is a linear approach.(2) The accuracy of the job cycle time forecasting by the proposed fuzzy data mining approach is better than those of the compared approaches because the proposed approach has achieved 25%~85% reductions in MAPE. (3) Moreover, the precision level of the proposed methodology, measured in terms of the average range of forecasts, is also superior to those of the compared approaches. The proposed methodology can be used to specify a very narrow range of the job cycle time. Since the internal due date is based on the cycle time forecast, a narrow range promotes the reliability of the internal due date.(4) The proposed method is unusually slow in its modification of the lower bound (see Figure 7); a more effective way to solve this problem is needed.(5) PCA seems to be useful for methods with job classification like CART. The combination of PCA and CART has significantly better forecasting accuracy than either method in isolation.(6) The simple combination of PCA and BPN does not have much effect. The main effect of PCA is to improve the correctness of the job classification.

In order to test the performance of the proposed methodology on real world situations, it has been applied to schedule a problem of 457 jobs from a wafer fabrication factory located in Taichung City Scientific Park, Taiwan. There are more than ten products in the wafer fabrication factory. The wafer fabrication factory has a monthly capacity of 20,000 wafers. Jobs are regularly released into the wafer fabrication factory. Various types of priorities are assigned to these jobs. Each product has 150~200 steps and six to nine reentries to the most bottlenecked machine. The first-in first-out (FIFO) rule is used to dispatch the jobs. The production characteristic of “reentry,” which is highly relevant to the semiconductor industry, is clearly reflected in this problem. It also shows the difficulties facing production planners and schedulers who attempt to provide an accurate due date for a product with a very complicated routing. For each job, twelve parameters were collected or retrieved from production management information system (PROMIS) databases and reports. After a backward elimination by regression analysis, six parameters (the job size, factory utilization, the queue length on the route, the queue length before the bottleneck, the work-in-process (WIP), and the average waiting time) were chosen as the inputs to the FBPN. Here “the average waiting time” refers to the average of the waiting times of the three jobs that were completed most recently, not to the waiting time of the job under consideration. The proposed methodology and several existing methods have been applied to this case. The forecasting precision (measured in terms of the average range) and accuracy (measured in terms of RMSE) of these approaches are compared in Figures 8 and 9, respectively. The samples in this case exhibited notable variability, which limited the performance of the forecasting methods. Nevertheless, it is obvious that the proposed methodology still outperformed the six existing methods in both regards. It is also noteworthy that the methods based on job classification (such as CART, PCA-CART, FCM-BPN, and the proposed methodology) achieved better forecasting accuracy than methods without job classification. Such an advantage was further strengthened by variable replacement techniques like PCA.

Figure 8 

Comparison of the forecasting precision.

Figure 9 

Comparison of the forecasting precision.

4. Conclusions and Directions for Future Research

Job cycle time forecasting is an important task in any wafer fabrication factory. In this field, many data mining methods have been applied. A fuzzy data mining approach for precise and accurate job cycle time forecasting is presented in this study. The forecast of the job cycle time by the proposed methodology is a fuzzy value; to improve precision, an effective procedure tightens the upper and lower bounds of the forecast. The forecasting accuracy is enhanced by a combination of PCA, FCM, and FBPN.

This study presents an illustrative example and a real case with the data of 457 jobs from a wafer fabrication factory. According to the experimental results, we have the following.(1) the forecasting accuracy (measured with MAE, MAPE, and RMSE) of the PCA-FCM-FBPN was significantly better than the accuracy levels of the other approaches.(2) Methods based on job classification have better forecasting accuracy. The incorporation of PCA makes the relationships between variables clearer and strengthens the advantage of job classification.(3) It is possible to specify a very narrow range of the job cycle time using the proposed methodology.

The upper bound can be adjusted faster than the lower bound. An effective method is needed to tackle this issue. In addition, the FBPN approach focuses on the modification of the threshold on the output node. Similar concepts can be applied to modify other parameters in the FBPN.

Acknowledgment

This study was financially supported by the National Science Council of Taiwan.