# A Novel Fuzzy-Neural Slack-Diversifying Rule Based on Soft Computing Applications for Job Dispatching in a Wafer Fabrication Factory

Authors: Toly Chen; Richard Romanowski
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This study proposes a slack-diversifying fuzzy-neural rule to improve job dispatching in a wafer fabrication factory. Several soft computing techniques, including fuzzy classification and artificial neural network prediction, have been applied in the proposed methodology. A highly effective fuzzy-neural approach is applied to estimate the remaining cycle time of a job. This research presents empirical evidence of the relationship between the estimation accuracy and the scheduling performance. Because dynamic maximization of the standard deviation of schedule slack has been shown to improve performance, this work applies such maximization to a slack-diversifying fuzzy-neural rule derived from a two-factor tailored nonlinear fluctuation smoothing rule for mean cycle time (2f-TNFSMCT). The effectiveness of the proposed rule was checked with a simulated case, which provided evidence of the rule’s effectiveness. The findings in this research point to several directions that can be exploited in the future.


---

## Body

## 1. Introduction

Wafer fabrication is one of the most important steps in semiconductor manufacturing [1]. Wafer fabrication factories exist all over the world; many are in Taiwan. The production operations of a wafer fabrication factory are very expensive, and the factory must be fully utilized to stay in business. If the factory is to remain in operation, its capacity must not substantially exceed the demand. Factory managers must plan the use of the existing capacity to shorten the cycle time. To maximize the product turnover rate is an important goal. In this regard, scheduling is undoubtedly a very useful tool. Kim et al. [2] simultaneously considered three issues: release control, mask scheduling, and batch scheduling. However, Chen et al. [3–5] noted that job dispatching is very difficult task in a wafer fabrication factory. Traditionally, a scheduling problem is formulated as a mathematical programming problem. The optimal solution of the mathematical programming problem gives the optimal scheduling of the manufacturing system. However, the mathematical programming problem of scheduling a wafer fabrication factory is large and effectively intractable. In practice, many wafer fabrication factories suffer from lengthy cycle times and cannot speed up their deliveries to their customers.Semiconductor manufacturing can be divided into four stages: wafer fabrication, wafer probing, packaging, and final testing. The most important and most time-consuming stage is wafer fabrication. This study investigates job dispatching for this stage.An effective scheduling and dispatching algorithm is an urgent necessity for wafer fabrication. This field includes many different methods, including dispatching rules, heuristics, data mining-based approaches [6, 7], agent technologies [6, 8–10], and simulation. The prevalent methods for practical applications include dispatching rules (e.g., first-in first out (FIFO), earliest due date (EDD), least slack (LS), shortest processing time (SPT), shortest remaining processing time (SRPT), critical ratio (CR), the fluctuation smoothing rule for the mean cycle time (FSMCT), the fluctuation smoothing rule for cycle time variation (FSVCT), FIFO+, SRPT+, and SRPT++), all of which have received a lot of attention over the last few years [6–8]. For details on traditional dispatching rules, please refer to Lu et al. [11].Recently, Chen proposed the one-factor tailored nonlinear fluctuation smoothing rule for mean cycle time (1f-TNFSMCT). This rule contains an adjustable parameter that allows it to be customized for a target wafer fabrication factory. Chen et al. [12] proposed the two-factor tailored nonlinear fluctuation smoothing rule for mean cycle time (2f-TNFSMCT), which outperformed four existing rules in scheduling a wafer fabrication factory.Magnifying the difference in the slack seems to improve scheduling performance, especially with respect to the average cycle time. To exploit this advantage, Wang et al. [13] derived the slack-diversifying nonlinear fluctuation smoothing rule by diversifying the slack in the 1f-TNFSVCT rule. To extend this advantage and to enhance the scheduling of wafer fabrication factories, a slack-diversifying fuzzy-neural rule is proposed in this study. The objective function includes the average and standard deviation of the cycle time, and therefore the problem can be denoted by J/rj, dj/CTj¯, and σCTj.The proposed methodology applies several soft computing techniques, including fuzzy classification and artificial neural network prediction. First, the remaining cycle time of a job needs to be estimated by the slack-diversifying fuzzy-neural rule. In this research, an innovative and highly effective fuzzy-neural method estimates the remaining cycle time of a job. The fuzzy-neural approach is based on the fuzzy c-means and back propagation network (FCM-BPN) approach [14]. According to Chen and Wang [4], improvements to the accuracy of remaining cycle time estimation can significantly improve the scheduling performance of a fluctuation smoothing rule. In the original study, Chen and Wang used a time-consuming and not very accurate gradient search algorithm to train the BPN. In this study, we use the Levenberg-Marquardt algorithm [15] to achieve the same purpose; it is more efficient and more accurate than the algorithm in Chen and Wang’s study. In addition, we also found some empirical evidence regarding the relationship between the estimation accuracy and the scheduling performance.The slack-diversifying nonlinear fluctuation smoothing rule is modified from 2f-TNFSMCT [12]; the new rule maximizes the difference in the slack, as measured by the standard deviation of the slack. Slack is a fuzzy concept, and in this study it is defined in a way that is conducive to scheduling performance. The factor values for achieving this must be determined, and that calculation turns out to be a complex optimization problem. We applied a polynomial fitting technique to convert it into a more tractable form for which several optimal solutions can be found. After screening some values out of the specified range, we used the remaining values to construct an optimized 2f-TNFSMCT rule. It is possible that some jobs have very large or small slacks, which will distort the optimization results. For this reason, such jobs are excluded. Further, the values of parameters influence the range of the slack. For a fair comparison, the range of the slack should be considered to determine the optimal parameter values. The differences between the proposed methodology and Wang et al.’s method are summarized in Table 1.Table 1
The differences between the proposed methodology and Wang et al.’s method.
Method
Rule base
Degree of freedom
Exclude extreme values
Considering the range ofSKju
Wang et al.’s method
1f-TNFSVCT
1
No
No
The proposed methodology
2f-TNFSMCT
2
Yes
YesThis paper is arranged as follows. Section2 reviews the existing approaches to scheduling a wafer fabrication factory. Section 3 provides the details of the proposed methodology. In Section 4, a simulated case is used to validate the effectiveness of the slack-diversifying fuzzy-neural rule. The performance levels of some of the existing rules in this field are also tested with the simulated data. Finally, we draw our conclusions in Section 5 and discuss some worthwhile topics for future work.


## 2. Related Work

Some earlier work in this field is relevant. Mönch et al. [16] classified the scheduling problems in a semiconductor manufacturing factory into six categories: batch scheduling problems, parallel machine scheduling problems, job shop scheduling problems, scheduling problems with auxiliary resources, multiple orders per job scheduling problems, and scheduling problems related to cluster tools. In Yao et al. [17], a decentralised multiobjective scheduling methodology was presented for semiconductor manufacturing, in which global objectives were decentralised into local ones of work stations. Lee et al. [18] adopted the Petri nets to accurately model the semiconductor manufacturing activities. Through representing the token movements in a Petri net with the well-established scheduling model for batch chemical processes, the optimal schedule of the given semiconductor process could be determined accordingly. Yugma et al. [19] proposed an efficient heuristic algorithm based on iterative sampling and simulated annealing for solving a complex batching and scheduling problem in a diffusion area of a semiconductor plant. Altendorfer et al. [20] proposed the work in parallel queue (WIPQ) rule to maximize throughput with a low level of work in process (WIP). Zhang et al. [21] proposed the dynamic bottleneck detection (DBD) approach that classifies workstations into several categories and then applies different dispatching rules to these categories. They used three dispatching rules including FIFO, the shortest processing time until the next bottleneck (SPNB), and CR. In view of the uncertainty in the classification of workstations, Chen [22] proposed the fuzzy DBD approach.Considering the current conditions in a wafer fabrication factory, Hsieh et al. [7] chose one approach from FSMCT, FSVCT, largest deviation first (LDF), one step ahead (OSA), or FIFO. Chen [23] modified FSMCT and proposed the nonlinear FSMCT (NFSMCT) rule, in which he smoothed the fluctuation in the estimated remaining cycle time and balanced it with that of the release time or the mean release rate. To diversify the slack, he applied the “division” operator. This was followed by Chen [24], in which he proposed the one-factor tailored NFSMCT (1f-TNFSMCT) rule and the one-factor tailored nonlinear FSVCT (1f-TNFSVCT) rule. Both rules contain adjustable parameters that allow them to be customized for a target wafer fabrication factory. Dabbas and Fowler [25] and Dabbas et al. [26] combined some dispatching rules into a single rule by forming their linear combination with relative weights. However, that research lacked a systematic procedure to determine the weights of those rules. In a multiple-objective study, Chen and Wang [27] proposed a biobjective nonlinear fluctuation smoothing rule with an adjustable factor (1f-biNFS) to optimize both the average cycle time and the cycle time variation at the same time. More degrees of freedom seem to enhance the performance of customizable rules. For this reason, Chen et al. [14] extended 1f-biNFS to a bi-objective fluctuation smoothing rule with four adjustable factors (4f-biNFS). One drawback of these rules is that only static factors are used, and they must be determined in advance. To this end, most studies (e.g., [14, 23, 24, 27]) performed extensive simulations. This is not only time-consuming, but it also fails to consider a sufficient number of factor combinations.Chen [28] established a mechanism that was able to adjust the values of the factor in 1f-biNFS dynamically (dynamic 1f-biNFS). However, even though satisfactory results were obtained in his experiment, there was no theoretical basis supporting the proposed mechanism. Chen [29] attempted to relate the scheduling performance to the factor values using a back propagation network (BPN). If that had worked, then the factor values contributing to the optimal scheduling performance could have been found. However, the explanatory ability of the BPN was not sufficient.At the same time, Chen [28] stated that a nonlinear fluctuation smoothing rule uses the divisor operator instead of the subtraction operator, which diversifies the slack and makes the nonlinear fluctuation smoothing rule more responsive to changes in the parameters. Chen and Wang [27] proved that the effects of the parameters are balanced better in a nonlinear fluctuation smoothing rule than in a traditional one if the variation in the parameters is large.


## 3. Methodology

The variables and parameters that will be used in the proposed methodology are defined as follows:(1)
Rj: the release time of jobj, j=1~n;(2)
RCTEju: the estimated remaining cycle time of job j from step u;(3)
SKju: the slack of job j at step u;(4)
λ: the mean release rate;(5)
xjp: inputs to the three-layer BPN of job j, p=1~P;(6)
hl: the output from hidden-layer node l, l=1~L;(7)
wlo: the connection weight between hidden-layer node l and the output node;(8)
wplh: the connection weight between input node p and hidden-layer node l, p=1~P; l=1~L;(9)
θlh: the threshold on hidden-layer node l;(10)
θo: the threshold on the output node.The proposed methodology includes the following six steps.Step 1.
Normalize the collected data [30].Step 2.
Use FCM to classify jobs. The required inputs for this step are job attributes. To determine the optimal number of categories, we use theS test [31]. The output of this step is the category of each job.Step 3.
Use the BPN approach to estimate the remaining cycle time of each job. Jobs of different categories will be sent to different three-layer BPNs. The inputs to the three-layer BPN include the attributes of a job, while the output is the estimated remaining cycle time of the job.Step 4.
Derive the 2f-TNFSMCT rule.Step 5.
Diversify the slack in the 2f-TNFSMCT rule.Step 6.
Incorporate the estimated remaining cycle time into the new rule.The flowchart of the proposed methodology is shown in Figure1.Figure 1
The flowchart of the proposed methodology.Table2 is used to compare the proposed methodology with the existing methods.Table 2
Comparison of the proposed methodology and some existing methods.
Method
Static or dynamic
Number of objectives
Number of adjustable parameters
Optimized
1f-TNFSMCT
Static
1
1
No
1f-TNFSVCT
Static
1
1
No
1f-biNFS
Static
2
1
No
2f-TNFSMCT
Static
1
2
No
4f-biNFS
Static
2
4
No
Dynamic 1f-biNFS
Dynamic
2
1
No
The proposed methodology
Dynamic
2
2
YesThe remaining cycle time of a job being produced in a wafer fabrication factory is the time still needed to complete the job (see Figure2). If the job has just been released into the wafer fabrication factory, then the remaining cycle time of the job is its cycle time [32–38]. Tai et al. [32] provided a statistical approach to calculate the cycle time for multilayer semiconductor final testing involving the sum of multiple Weibull-distributed waiting times. The remaining cycle time is an important attribute (or performance measure) for the WIP in the wafer fabrication factory. We need to estimate the remaining cycle time for each job because the remaining cycle time is an important input to the scheduling rule. Past studies (e.g., [12]) have shown that the accuracy of remaining cycle time estimation can be improved by job classification. Soft computing methods (e.g., [4, 12, 37, 38]) have received much attention in this regard.Figure 2
The concept of the remaining cycle time.


### 3.1. Step 1: Normalize the Collected Data

First, in order to facilitate the subsequent calculations and problem solving, all collected dataxj=[xjp] are normalized into [0.10.9] [38] as follows:
(1)zjp=NL+xjp-min(xjp)max(xjp)-min(xjp)·(NU-NL),
where zjp is the normalized value of xjp and NL and NU indicate the lower and upper bounds of the range of the normalized value, respectively. min(xjp) and max(xjp) are the minimum and maximum of xj, respectively. The formula can be written as
(2)xjp=zjp-NLNU-NL·(max(xjp)-min(xjp))+min(xjp),if the unnormalized value is to be obtained. Then, we place the (normalized) attributes of a job in vectorzj=[zjp]. To illustrate the proposed methodology, a real case containing the data of 35 jobs was used. For each job, twelve attributes were collected from the reports and databases of a production management information system (PROMIS). After a backward elimination by regression analysis, six attributes that were the most influential to the job cycle time were chosen (see Table 3).Table 3
An illustrative example.
Jobnumber
u
x
j
1
x
j
2
x
j
3
x
j
4
x
j
5
x
j
6
RCTE
j
u
1
197
96%
22
9
35
63
98
1181
2
192
97%
26
9
35
63
257
1194
3
188
96%
14
8
36
65
99
1260
4
202
97%
11
9
32
61
131
1240
5
187
97%
26
8
37
66
247
1180
6
184
97%
19
9
39
68
191
1227
7
197
96%
20
9
38
63
219
1236
8
184
96%
25
9
40
66
219
1215
9
178
97%
33
9
41
67
219
1228
10
212
96%
31
9
38
66
54
1266
11
177
96%
36
9
39
67
54
1285
12
186
96%
20
9
37
65
54
1272
13
199
96%
19
9
36
66
54
1310
14
195
96%
11
9
37
66
49
1265
15
223
97%
12
8
34
63
201
1308
16
223
96%
32
8
36
62
103
1331
17
206
95%
30
9
33
59
53
1294
18
212
95%
18
9
35
64
53
1314
19
168
97%
24
9
34
60
248
1321
20
200
97%
20
9
34
60
248
1353
21
140
95%
34
9
35
60
82
1226
22
156
96%
10
8
36
60
98
1301
23
206
96%
19
9
38
67
67
1280
24
152
96%
23
9
40
69
67
1286
25
187
95%
24
8
39
68
67
1252
26
185
97%
32
9
39
64
223
1214
27
137
97%
28
9
39
66
176
1251
28
151
98%
19
9
40
67
462
1222
29
148
97%
9
9
39
66
168
1187
30
173
96%
20
9
39
68
141
1205
31
161
96%
16
9
39
69
95
1120
32
108
97%
19
9
40
68
179
1133
33
115
96%
35
9
41
71
209
1130
34
100
96%
22
9
38
67
237
1113
35
97
99%
19
9
41
67
684
1107The results of partial normalization are shown in Table4.Table 4
The results of partial normalization.
Jobnumber
z
j
1
z
j
2
z
j
3
z
j
4
z
j
5
z
j
6
1
0.300
0.485
0.900
0.367
0.367
0.162
2
0.500
0.604
0.900
0.367
0.367
0.362
3
0.300
0.248
0.100
0.456
0.500
0.163
4
0.500
0.159
0.900
0.100
0.233
0.203
5
0.500
0.604
0.100
0.544
0.567
0.349
6
0.500
0.396
0.900
0.722
0.700
0.279
7
0.300
0.426
0.900
0.633
0.367
0.314
8
0.300
0.574
0.900
0.811
0.567
0.314
9
0.500
0.811
0.900
0.900
0.633
0.314
10
0.300
0.752
0.900
0.633
0.567
0.106
11
0.300
0.900
0.900
0.722
0.633
0.106
12
0.300
0.426
0.900
0.544
0.500
0.106
13
0.300
0.396
0.900
0.456
0.567
0.106
14
0.300
0.159
0.900
0.544
0.567
0.100
15
0.500
0.189
0.100
0.278
0.367
0.291
16
0.300
0.781
0.100
0.456
0.300
0.168
17
0.100
0.722
0.900
0.189
0.100
0.105
18
0.100
0.367
0.900
0.367
0.433
0.105
19
0.500
0.544
0.900
0.278
0.167
0.351
20
0.500
0.426
0.900
0.278
0.167
0.351
21
0.100
0.841
0.900
0.367
0.167
0.142
22
0.300
0.130
0.100
0.456
0.167
0.162
23
0.300
0.396
0.900
0.633
0.633
0.123
24
0.300
0.515
0.900
0.811
0.767
0.123
25
0.100
0.544
0.100
0.722
0.700
0.123
26
0.500
0.781
0.900
0.722
0.433
0.319
27
0.500
0.663
0.900
0.722
0.567
0.260
28
0.700
0.396
0.900
0.811
0.633
0.620
29
0.500
0.100
0.900
0.722
0.567
0.250
30
0.300
0.426
0.900
0.722
0.700
0.216
31
0.300
0.307
0.900
0.722
0.767
0.158
32
0.500
0.396
0.900
0.811
0.700
0.264
33
0.300
0.870
0.900
0.900
0.900
0.302
34
0.300
0.485
0.900
0.633
0.633
0.337
35
0.900
0.396
0.900
0.900
0.633
0.900


### 3.2. Step 2: Classify Jobs Using FCM

In the proposed methodology, jobs are classified intoK categories using FCM. If a crisp clustering methods were applied, then it would be possible that some clusters would have very few examples. By contrast, in a fuzzy clustering method, an example belongs to multiple clusters to different degrees, which provides a solution to this problem. Similarly, in probability theory the naïve Bayes method provides the probability that a given item belongs to each class. However, the application of FCM can consider subjective issues in job classification.FCM classifies jobs by minimizing the following objective function:(3)Min∑k=1K∑j=1nμj(k)mej(k)2,
where K is the required number of categories; n is the number of jobs; μj(k) indicates the degree of membership with which job j belongs to category k; ej(k) measures the distance from job j to the centroid of category k; m∈[1,∞) is a parameter to adjust the fuzziness and is usually set to 2. The procedure of FCM is as follows:(1)
produce a preliminary clustering result;(2)
(for some number of iterations) calculate the centroid of each category as(4)z-(k)={z-(k)p},k=1~K,z-(k)p=∑j=1nμj(k)mzjp∑j=1nμj(k)m,k=1~K,p=1~P,μj(k)=1∑g=1K(ej(k)/ej(g))2/(m-1),j=1~n,k=1~K,ej(k)=∑p=1P(zjp-z-(k)p)2,j=1~n,k=1~K,
where z-(k) is the centroid of category k. μj(k)(t) is the degree of membership that job j holds in category k after the tth iteration;(3)
remeasure the distance from each job to the centroid of each category and then recalculate the corresponding membership;(4)
stop if the following condition is met. Otherwise, return to step(3) as follows:
(5)maxkmaxj|μj(k)(t)-μj(k)(t-1)|<d,
where d is a real number representing the threshold for the convergence of membership.Finally, the separate distance test (S test) proposed by Xie and Beni [31] can be applied to determine the optimal number of categories K as follows:
(6)minSsubjecttoJm=∑k=1K∑j=1nμj(k)mej(k)2,emin2=mink1≠k2(∑p=1P(z-(k1)p-z-(k2)p)2),S=Jmn×emin2,K∈Z+.The Fuzzy Logic Toolbox of MATLAB is used to implement the FCM approach. The FCM program code for the illustrative example is shown in Algorithm1. The S value can also be obtained using this program.Algorithm 1:The FCM program code in MATLAB.
A = [0.300  0.485  0.900  0.367  0.367  0.162; …;  0.900  0.396  0.900  0.900 0.633  0.900]
K = 4
[
center,  U,  obj_fun] = fcm(A,  K);
Jm  =  min(obj_fun)
e2_min  =  9999;
for  i  =  1:  K
for  j  =  i  +  1 : K
e2_sum  =  0;
for  p  =  1:  3
e2_sum  =  e2_sum + (center(i,  p)   −  center(j,  p))∧2;
end
if  e2_sum  <  e2_min
e2_min  =  e2_sum;
end
end
end
e2_min
S  =  min(Jm)/(35*e2_min)The results of theS test are summarized in Table 5. In this case, the optimal number of job categories was 3. The threshold for the convergence of membership (d) was set to 0.01. A common practice is to set a threshold of membership μL to determine whether a job belongs to each category. For example, the classification results for the assumption that μL=0.5 are shown in Table 6. With each decrease in the threshold, each category will contain more jobs. For example, the classification results for the assumption that μL=0.1 are shown in Table 7. This can solve the problem of an insufficient number of examples. In addition, a job may be classified under more than one category, which makes FCM different from crisp classification methods. The disadvantage of this approach is that some jobs do not belong to any category, that is, the outliers. Another way is to assign a job to the category for which it has the highest degree of membership. The disadvantage of this approach is that some categories will have very few jobs.Table 5
The results of theS test.
Number of categories (K)
J
m
e
min
2
S
2
0.0739
0.0027
0.7946
3
0.0444
0.0028
0.4586
4
0.0316
0.0014
0.6298
5
0.0247
0.0014
0.4882
6
0.0205
0.0009
0.6336
7
0.0175
0.0005
1.1102
8
0.0154
0.0003
1.2577Table 6
The classifying results (μL=0.5).
Category
Jobs
1
6, 8, 9, 10, 11, 23, 24, 26, 27, 28, 29, 30, 31, 32, 33, 34
2
3, 5, 15, 16, 22, 25
3
1, 2, 4, 7, 12, 13, 17, 18, 19, 20, 21Table 7
The classifying results (μL=0.1).
Category
Jobs
1
2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35
2
3, 4, 5, 11, 14, 15, 16, 17, 20, 21, 22, 25, 28, 29, 33, 35
3
1, 2, 4, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 31, 33, 34, 35The application of FCM has the following problems:(1)
how to determine the membership thresholdμL. In theory, setting μL to 1/K ensures that each job can be classified under some category with absolute certainty. Conversely, if μL is set to a large value, it is possible that some jobs cannot be classified under any category with absolute certainty;(2)
whether to consider the degree of membership in the remaining cycle time estimation. If two jobs belong to the same category but have different degrees of membership, how can one guarantee that the BPN of this category has the same predictive power for the two jobs? To address this issue, the degree of membership can be considered when estimating the remaining cycle time with the BPN.


### 3.3. Step 3: Estimate the Remaining Cycle Time with a BPN

After clustering, for each category, a three-layer BPN is built to estimate the remaining cycle times of jobs in the category. Some of the jobs in the category are input as “training examples” to the BPN to determine the parameter values. The configuration of the three-layer BPN is as follows. First, the inputs are theP parameters associated with the jth job. These parameters have to be normalized before they are fed into the three-layer BPN. Subsequently, there is only a single hidden layer with neurons that are twice as many as the input layer. Finally, the output from the three-layer BPN is the (normalized) remaining cycle time estimate (RCTEj) of the example. The activation function used in each layer is the Sigmoid function,
(7)f(x)=11+e-x.The procedure for determining the parameter values is now described. Two phases are involved at the training stage. At first, in the forward phase, inputs are multiplied with weights, summed, and transferred to the hidden layer. Then activated signals are output from the hidden layer as(8)hl=11+e-nlh,
where
(9)nlh=Ilh-θlh,Ilh=∑p=1Pwplhxp,
where hl values are also transferred to the output layer with the same procedure. Finally, the output of the BPN is generated as
(10)o=11+e-no,
where
(11)no=Io-θo,Io=∑l=1Lwlohl.
Some algorithms are applicable for training a BPN in the backward phase, such as gradient descent algorithms, 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 in the following.The network parameters are placed in vectorβ=[w11h,…,wPLh,θ1h,…,θLh,w1o,…,wLo,θo]. The network output oj can be represented with f(xj,β). The objective function of the BPN is to minimize RMSE or equivalently the sum of squared error (SSE) as follows:
(12)SSE(β)=∑j=1n(N(RCTEju)-f(xj,β))2.The Levenberg-Marquardt algorithm is an iterative procedure. At the beginning, the user should specify the initial values of the network parametersβ. It is common practice to set βT=(1,1,…,1). At each step, the parameter vector β is replaced by a new estimate β+δ, where δ=[Δw11h,…,ΔwPLh,Δθ1h,…,ΔθLh,Δw1o,…,ΔwLo,Δθo]. The network output becomes f(xj,β+δ); it is approximated by its linearization as
(13)f(xj,β+δ)≈f(xj,β)+Jjδ,
where
(14)Jj=∂f(xj,β)∂β
is the gradient vector of f with respect to β. Substituting (13) into (12),
(15)SSE(β+δ)≈∑j=1n(N(RCTEju)-f(xj,β)-Jjδ)2.
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 give
(16)(JTJ)δ=JT(N(RCTEju)-f(xj,β)),
where J is the Jacobian matrix containing the first derivative of network error with respect to the weights and biases. Equation (16) includes a set of linear equations that can be solved for δ.Finally, the three-layer BPN can be applied to estimate the remaining cycle time of a job. When a new job is released into the factory, the parameters associated with the new job are recorded. Then the new job is classified into a category, and the three-layer BPN of the category can be applied to estimate the remaining cycle time of the new job.Consider category 1 of the previous example. WhenμL=0.5, there are 16 jobs in this category. These jobs are split into two parts: the training data (the first 12 jobs) and the testing data (the remaining jobs). A three-layer BPN estimates the remaining cycle time of jobs in this category according to their six attributes with the following settings: single hidden layer,the number of neurons in the hidden layer: 12,convergence criterion: mean squared error (MSE)<10-4.The Neural Network Toolbox of MATLAB is used to implement the BPN approach. The BPN program code for the illustrative example is shown in Algorithm2. The estimation results are shown in Figure 3.Algorithm 2:The BPN code in MATLAB (category 1).
t
n_input  =    [0.500  0.300  0.500…;  0.396  0.574  0.811…;  0.900  0.900  0.900…;  0.722
0.811  0.900…;  0.700  0.567  0.633…;  0.279  0.314  0.314…]
tn_target  =  [0.341  0.383  0.598⋯0.100]
net    =    newff([0  1;  0  1;  0  1;  0  1;  0  1;  0  1],[12,  1],{“logsig”,  “logsig”},  “trainlm”);
net    =    init(net);
net.trainParam.show  =  10;
net.trainParam.lr  =  0.1;
net.trainParam.epochs  =  1000;
net.trainParam.goal  =  1e  −  4;
[
net,  tr]  =  train(net,  tn_input,  tn_target);
tn_output  =  sim(net,  tn_input)
te_input  =  [0.300  0.500  0.300…;  0.307  0.396  0.870…;  0.900  0.900  0.900…;
0.722  0.811  0.900…;  0.767  0.700  0.900…;  0.158  0.264  0.302…]
te_output  =  sim(net,  te_input)Figure 3
The estimation results by FCM-BPN.The following indexes evaluate the estimation accuracy of both training and testing data:mean absolute error (MAE) = 43 (hrs),mean absolute percentage error (MAPE) = 3.7%,RMSE = 62 (hrs).In contrast, if the remaining cycle times of these jobs are predicted in association with those of the other unclassified jobs, then the estimation performance will be mean absolute error (MAE) = 47 (hrs),mean absolute percentage error (MAPE) = 4.1%,root mean squared error (RMSE) = 70 (hrs)which is much poorer. Table 8 compares the performance levels of the gradient descent algorithm and the Levenberg-Marquardt algorithm.Table 8
Training algorithm performance (convergence criterion: MSE < 10−4).
Algorithm
Number of epochs
MAE
MAPE
RMSE
Gradient descent
10000
79
6.7%
98
Levenberg-Marquardt
<100
47
4.1%
70How does the accuracy of remaining cycle time estimation affect the performance of the dispatching rule? To answer this question, we add a noise factor to the remaining cycle time before feeding it to the dispatching rule as follows:(17)RCTEju+εju,εju∈N(0~200,0~50),
where εju follows the normal distribution. After a simulation study using FSMCT, the estimated remaining cycle time was compared with the actual value to evaluate the estimation accuracy. The relationship between the scheduling performance and the estimation accuracy is shown in Figure 4. Forty (40) replications of a simulation experiment were run with FSMCT. The scheduling performance and the estimation accuracy were evaluated with the average cycle time and RMSE, respectively. Obviously, the more accurately the remaining cycle time can be estimated, the better the schedule will perform. Therefore, the efforts made in this paper to improve the accuracy of remaining cycle time estimation make sense.Figure 4
The relationship between the scheduling performance and the estimation accuracy.


### 3.4. The Slack-Diversifying Fuzzy-Neural Rule

In traditional fluctuation smoothing (FS) there are two different formulation methods with two separate strengths [11]. One method is aimed at minimizing the mean cycle time with FSMCT as follows:
(18)SKju(FSMCT)=jλ-RCTEju.
The other method is aimed at minimizing the variance of cycle time with FSVCT as follows:
(19)SKju(FSVCT)=Rj-RCTEju.
Jobs with the smallest slack values are given high priority. These two rules and their variants have been proven to be very effective in shortening the cycle times of wafer fabrication factories [11, 14, 24, 27]. In the traditional FSMCT rule, RCTEju might be much greater than j/λ. As a result, the slack of a job is determined solely by RCTEju. To tackle this problem, both terms in the FSMCT rule are normalized as follows:
(20)Nor(RCTEju)=(RCTEju-min(RCTEju))(max(RCTEju)-min(RCTEju)),Nor(jλ)=(j/λ-min(j/λ))(max(j/λ)-min(j/λ))=(j/λ-1/λ)(n/λ-1/λ)=(j-1)(n-1).
After normalization, both terms now range from 0 to 1. Subsequently, to improve the responsiveness of the FSMCT rule, the division operator is applied instead of the traditional subtraction operator as follows:
(21)nonlinearFSMCTrule:SKju=Nor(j/λ)Nor(RCTEju).
However, this rule cannot be tailored to the wafer fabrication factory that is to be scheduled. To address this problem, the transition from a traditional FSMCT rule to its nonlinear form is analyzed as follows. The nonlinear form can be rewritten as(22)SKju=Nor(j/λ)Nor(RCTEju)=(j-1)/(n-1)(RCTEju-min(RCTEju))/(max(RCTEju)-min(RCTEju))=βγ·j-1RCTEju-min(RCTEju)=βλγ·(j/λ)-(1/λ)RCTEju-min(RCTEju)=βλγ(RCTEju-min(RCTEju))·(jλ-RCTEju+RCTEju-1λ)=(γ(RCTEju-min(RCTEju))βλ)-1·(jλ-RCTEju+(RCTEju-min(jλ))·1),where β=max(RCTEju)-min(RCTEju) and γ=n-1. Conversely, the linear form can also be rewritten as
(23)SKju=jλ-RCTEju=(γ(RCTEju-min(RCTEju))βλ)-0·(jλ-RCTEju+(RCTEju-min(jλ))·0).
These two formulas can be generalized into the following form:
(24)SKju=(γ(RCTEju-min(RCTEju))βλ)-ξ·(jλ-RCTEju+(RCTEju-min(jλ))·ζ).Wang et al. derived [13] the slack-diversifying nonlinear fluctuation smoothing rule by diversifying the slack in the 1f-TNFSVCT rule. In this study, we diversify the slack in the 2f-TNFSMCT rule as follows:
(25)SKju=(βλγ(RCTEju-min(RCTEju)))ξ·(jλ-RCTEju+ζ(RCTEju-1λ))=ajuξ(bju+cjuζ)=bjuajuξ+cjuζajuξ,
where
(26)aju=βλγ(RCTEju-min(RCTEju)),bju=jλ-RCTEju,cju=RCTEju-1λ.
There are many possible models that can form the combination of ξ and ζ. For example,
(27)(Linearmodel)ξ=ζ,(28)(Nonlinearmodel)ξ=ζk,k≥0,(29)(Logarithmicmodel)ξ=ln(1+ζ)ln2.
However, (25) is difficult to deal with. For this reason, the following polynomial fitting techniques are used to convert it into a more tractable form:
(30)xξ≅(1.15-0.11ξ)+(-0.07+0.57ξ)x.
The mean absolute percentage error (RMSE) of (30) is less than 9% when x≤10. RMSE will not be a serious problem, since it is the ξ value associated with the maximum σSKju that is to be found, not the SKju value. These polynomial fitting techniques are especially effective when x exceeds 1.Applying (30) to (25) gives
(31)SKju≅bju(1.15-0.11ξ+(-0.07+0.57ξ)aju)+cjuζ(1.15-0.11ξ+(-0.07+0.57ξ)aju)=(1.15bju-0.07ajubju)+(-0.11bju+0.57ajubju)ξ+(1.15cju-0.07ajucju)ζ+(-0.11cju+0.57ajucju)ξζ=dju+ejuξ+fjuζ+gjuξζ,
where
(32)dju=1.15bju-0.07ajubju,eju=-0.11bju+0.57ajubju,fju=1.15cju-0.07ajucju,gju=-0.11cju+0.57ajucju.
To diversify the slack, the standard deviation of the slack is to be maximized as follows:
(33)σSKju=∑j=1n(SKju-SKju¯)2n-1=1(n-1)∑j=1nSKju2-1n(∑j=1nSKju)2
which is equivalent to maximizing the following term:
(34)∑j=1nSKij2-1n(∑j=1nSKij)2=∑j=1n(dju+ejuξ+fjuζ+gjuξζ)2-1n(∑i=1n(dju+ejuξ+fjuζ+gjuξζ))2=w1+w2ζ+w3ξ+w4ζ2+w5ξ2+w6ζξ+w7ζξ2+w8ζ2ξ.
By taking the derivatives of (34) with respect to ζ and ξ and setting them equal to zero, we can obtain
(35)w2+2w4ζ+w6ξ+w7ξ2+2w8ζξ=0,(36)w3+2w5ξ+w6ζ+2w7ζξ+w8ζ2=0.
In addition, ζ and ξ must follow some model, for example, (27), (28), or (29). For example, if the nonlinear model in (28) is satisfied (k=2), then according to (35),
(37)w2+2w4ζ+w6ζ2+2w8ζ3+w7ζ4=0;ξ=ζ2,
while according to (36)
(38)w3+w6ζ+(2w5+w8)ζ2+2w7ζ3=0;ξ=ζ2.
The general solutions are too long to be presented here. Further, to guarantee a maximum, the second-order derivative must be negative as follows:
(39)2w4+2w8ξ≤0,2w5+2w7ζ≤0.
In addition, it is possible that some jobs have very large or very small slacks, either of which will distort the calculation of σSKju. For this reason, such jobs are excluded. Further, the values of ξ and ζ influence the range of SKju. For a fair comparison, σSKju should be divided by the range of SKju to determine the optimal values of ξ and ζ.An example is given in Table9 to illustrate the procedure mentioned above, in which α=244, β=195, and λ=1.18. Assume that the nonlinear model with k=2 is used for the relationship between ξ and ζ. The optimal solution is ξ* = 0.53 and ζ*=0.28 with the maximum σSKju equal to 9.21. In Figure 5 the results of the proposed methodology are compared with those of Wang et al.’s method. Obviously, the proposed methodology performed better at diversifying the slacks of jobs, which reduces the risk of misscheduling and promotes scheduling performance.Table 9
A demonstrative example.
j
R
j
RCTE
j
u
1
158
54
2
154
109
3
161
118
4
131
52
5
162
198
6
248
25
7
220
11
8
87
106
9
22
31
10
178
118
11
229
95
12
148
105
13
8
8
14
228
55
15
242
132
16
99
51
17
4
3
18
246
154
19
191
157
20
131
148Comparison of the results by the two methods.
(a)
(b)


## 3.1. Step 1: Normalize the Collected Data

First, in order to facilitate the subsequent calculations and problem solving, all collected dataxj=[xjp] are normalized into [0.10.9] [38] as follows:
(1)zjp=NL+xjp-min(xjp)max(xjp)-min(xjp)·(NU-NL),
where zjp is the normalized value of xjp and NL and NU indicate the lower and upper bounds of the range of the normalized value, respectively. min(xjp) and max(xjp) are the minimum and maximum of xj, respectively. The formula can be written as
(2)xjp=zjp-NLNU-NL·(max(xjp)-min(xjp))+min(xjp),if the unnormalized value is to be obtained. Then, we place the (normalized) attributes of a job in vectorzj=[zjp]. To illustrate the proposed methodology, a real case containing the data of 35 jobs was used. For each job, twelve attributes were collected from the reports and databases of a production management information system (PROMIS). After a backward elimination by regression analysis, six attributes that were the most influential to the job cycle time were chosen (see Table 3).Table 3
An illustrative example.
Jobnumber
u
x
j
1
x
j
2
x
j
3
x
j
4
x
j
5
x
j
6
RCTE
j
u
1
197
96%
22
9
35
63
98
1181
2
192
97%
26
9
35
63
257
1194
3
188
96%
14
8
36
65
99
1260
4
202
97%
11
9
32
61
131
1240
5
187
97%
26
8
37
66
247
1180
6
184
97%
19
9
39
68
191
1227
7
197
96%
20
9
38
63
219
1236
8
184
96%
25
9
40
66
219
1215
9
178
97%
33
9
41
67
219
1228
10
212
96%
31
9
38
66
54
1266
11
177
96%
36
9
39
67
54
1285
12
186
96%
20
9
37
65
54
1272
13
199
96%
19
9
36
66
54
1310
14
195
96%
11
9
37
66
49
1265
15
223
97%
12
8
34
63
201
1308
16
223
96%
32
8
36
62
103
1331
17
206
95%
30
9
33
59
53
1294
18
212
95%
18
9
35
64
53
1314
19
168
97%
24
9
34
60
248
1321
20
200
97%
20
9
34
60
248
1353
21
140
95%
34
9
35
60
82
1226
22
156
96%
10
8
36
60
98
1301
23
206
96%
19
9
38
67
67
1280
24
152
96%
23
9
40
69
67
1286
25
187
95%
24
8
39
68
67
1252
26
185
97%
32
9
39
64
223
1214
27
137
97%
28
9
39
66
176
1251
28
151
98%
19
9
40
67
462
1222
29
148
97%
9
9
39
66
168
1187
30
173
96%
20
9
39
68
141
1205
31
161
96%
16
9
39
69
95
1120
32
108
97%
19
9
40
68
179
1133
33
115
96%
35
9
41
71
209
1130
34
100
96%
22
9
38
67
237
1113
35
97
99%
19
9
41
67
684
1107The results of partial normalization are shown in Table4.Table 4
The results of partial normalization.
Jobnumber
z
j
1
z
j
2
z
j
3
z
j
4
z
j
5
z
j
6
1
0.300
0.485
0.900
0.367
0.367
0.162
2
0.500
0.604
0.900
0.367
0.367
0.362
3
0.300
0.248
0.100
0.456
0.500
0.163
4
0.500
0.159
0.900
0.100
0.233
0.203
5
0.500
0.604
0.100
0.544
0.567
0.349
6
0.500
0.396
0.900
0.722
0.700
0.279
7
0.300
0.426
0.900
0.633
0.367
0.314
8
0.300
0.574
0.900
0.811
0.567
0.314
9
0.500
0.811
0.900
0.900
0.633
0.314
10
0.300
0.752
0.900
0.633
0.567
0.106
11
0.300
0.900
0.900
0.722
0.633
0.106
12
0.300
0.426
0.900
0.544
0.500
0.106
13
0.300
0.396
0.900
0.456
0.567
0.106
14
0.300
0.159
0.900
0.544
0.567
0.100
15
0.500
0.189
0.100
0.278
0.367
0.291
16
0.300
0.781
0.100
0.456
0.300
0.168
17
0.100
0.722
0.900
0.189
0.100
0.105
18
0.100
0.367
0.900
0.367
0.433
0.105
19
0.500
0.544
0.900
0.278
0.167
0.351
20
0.500
0.426
0.900
0.278
0.167
0.351
21
0.100
0.841
0.900
0.367
0.167
0.142
22
0.300
0.130
0.100
0.456
0.167
0.162
23
0.300
0.396
0.900
0.633
0.633
0.123
24
0.300
0.515
0.900
0.811
0.767
0.123
25
0.100
0.544
0.100
0.722
0.700
0.123
26
0.500
0.781
0.900
0.722
0.433
0.319
27
0.500
0.663
0.900
0.722
0.567
0.260
28
0.700
0.396
0.900
0.811
0.633
0.620
29
0.500
0.100
0.900
0.722
0.567
0.250
30
0.300
0.426
0.900
0.722
0.700
0.216
31
0.300
0.307
0.900
0.722
0.767
0.158
32
0.500
0.396
0.900
0.811
0.700
0.264
33
0.300
0.870
0.900
0.900
0.900
0.302
34
0.300
0.485
0.900
0.633
0.633
0.337
35
0.900
0.396
0.900
0.900
0.633
0.900


## 3.2. Step 2: Classify Jobs Using FCM

In the proposed methodology, jobs are classified intoK categories using FCM. If a crisp clustering methods were applied, then it would be possible that some clusters would have very few examples. By contrast, in a fuzzy clustering method, an example belongs to multiple clusters to different degrees, which provides a solution to this problem. Similarly, in probability theory the naïve Bayes method provides the probability that a given item belongs to each class. However, the application of FCM can consider subjective issues in job classification.FCM classifies jobs by minimizing the following objective function:(3)Min∑k=1K∑j=1nμj(k)mej(k)2,
where K is the required number of categories; n is the number of jobs; μj(k) indicates the degree of membership with which job j belongs to category k; ej(k) measures the distance from job j to the centroid of category k; m∈[1,∞) is a parameter to adjust the fuzziness and is usually set to 2. The procedure of FCM is as follows:(1)
produce a preliminary clustering result;(2)
(for some number of iterations) calculate the centroid of each category as(4)z-(k)={z-(k)p},k=1~K,z-(k)p=∑j=1nμj(k)mzjp∑j=1nμj(k)m,k=1~K,p=1~P,μj(k)=1∑g=1K(ej(k)/ej(g))2/(m-1),j=1~n,k=1~K,ej(k)=∑p=1P(zjp-z-(k)p)2,j=1~n,k=1~K,
where z-(k) is the centroid of category k. μj(k)(t) is the degree of membership that job j holds in category k after the tth iteration;(3)
remeasure the distance from each job to the centroid of each category and then recalculate the corresponding membership;(4)
stop if the following condition is met. Otherwise, return to step(3) as follows:
(5)maxkmaxj|μj(k)(t)-μj(k)(t-1)|<d,
where d is a real number representing the threshold for the convergence of membership.Finally, the separate distance test (S test) proposed by Xie and Beni [31] can be applied to determine the optimal number of categories K as follows:
(6)minSsubjecttoJm=∑k=1K∑j=1nμj(k)mej(k)2,emin2=mink1≠k2(∑p=1P(z-(k1)p-z-(k2)p)2),S=Jmn×emin2,K∈Z+.The Fuzzy Logic Toolbox of MATLAB is used to implement the FCM approach. The FCM program code for the illustrative example is shown in Algorithm1. The S value can also be obtained using this program.Algorithm 1:The FCM program code in MATLAB.
A = [0.300  0.485  0.900  0.367  0.367  0.162; …;  0.900  0.396  0.900  0.900 0.633  0.900]
K = 4
[
center,  U,  obj_fun] = fcm(A,  K);
Jm  =  min(obj_fun)
e2_min  =  9999;
for  i  =  1:  K
for  j  =  i  +  1 : K
e2_sum  =  0;
for  p  =  1:  3
e2_sum  =  e2_sum + (center(i,  p)   −  center(j,  p))∧2;
end
if  e2_sum  <  e2_min
e2_min  =  e2_sum;
end
end
end
e2_min
S  =  min(Jm)/(35*e2_min)The results of theS test are summarized in Table 5. In this case, the optimal number of job categories was 3. The threshold for the convergence of membership (d) was set to 0.01. A common practice is to set a threshold of membership μL to determine whether a job belongs to each category. For example, the classification results for the assumption that μL=0.5 are shown in Table 6. With each decrease in the threshold, each category will contain more jobs. For example, the classification results for the assumption that μL=0.1 are shown in Table 7. This can solve the problem of an insufficient number of examples. In addition, a job may be classified under more than one category, which makes FCM different from crisp classification methods. The disadvantage of this approach is that some jobs do not belong to any category, that is, the outliers. Another way is to assign a job to the category for which it has the highest degree of membership. The disadvantage of this approach is that some categories will have very few jobs.Table 5
The results of theS test.
Number of categories (K)
J
m
e
min
2
S
2
0.0739
0.0027
0.7946
3
0.0444
0.0028
0.4586
4
0.0316
0.0014
0.6298
5
0.0247
0.0014
0.4882
6
0.0205
0.0009
0.6336
7
0.0175
0.0005
1.1102
8
0.0154
0.0003
1.2577Table 6
The classifying results (μL=0.5).
Category
Jobs
1
6, 8, 9, 10, 11, 23, 24, 26, 27, 28, 29, 30, 31, 32, 33, 34
2
3, 5, 15, 16, 22, 25
3
1, 2, 4, 7, 12, 13, 17, 18, 19, 20, 21Table 7
The classifying results (μL=0.1).
Category
Jobs
1
2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35
2
3, 4, 5, 11, 14, 15, 16, 17, 20, 21, 22, 25, 28, 29, 33, 35
3
1, 2, 4, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 31, 33, 34, 35The application of FCM has the following problems:(1)
how to determine the membership thresholdμL. In theory, setting μL to 1/K ensures that each job can be classified under some category with absolute certainty. Conversely, if μL is set to a large value, it is possible that some jobs cannot be classified under any category with absolute certainty;(2)
whether to consider the degree of membership in the remaining cycle time estimation. If two jobs belong to the same category but have different degrees of membership, how can one guarantee that the BPN of this category has the same predictive power for the two jobs? To address this issue, the degree of membership can be considered when estimating the remaining cycle time with the BPN.


## 3.3. Step 3: Estimate the Remaining Cycle Time with a BPN

After clustering, for each category, a three-layer BPN is built to estimate the remaining cycle times of jobs in the category. Some of the jobs in the category are input as “training examples” to the BPN to determine the parameter values. The configuration of the three-layer BPN is as follows. First, the inputs are theP parameters associated with the jth job. These parameters have to be normalized before they are fed into the three-layer BPN. Subsequently, there is only a single hidden layer with neurons that are twice as many as the input layer. Finally, the output from the three-layer BPN is the (normalized) remaining cycle time estimate (RCTEj) of the example. The activation function used in each layer is the Sigmoid function,
(7)f(x)=11+e-x.The procedure for determining the parameter values is now described. Two phases are involved at the training stage. At first, in the forward phase, inputs are multiplied with weights, summed, and transferred to the hidden layer. Then activated signals are output from the hidden layer as(8)hl=11+e-nlh,
where
(9)nlh=Ilh-θlh,Ilh=∑p=1Pwplhxp,
where hl values are also transferred to the output layer with the same procedure. Finally, the output of the BPN is generated as
(10)o=11+e-no,
where
(11)no=Io-θo,Io=∑l=1Lwlohl.
Some algorithms are applicable for training a BPN in the backward phase, such as gradient descent algorithms, 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 in the following.The network parameters are placed in vectorβ=[w11h,…,wPLh,θ1h,…,θLh,w1o,…,wLo,θo]. The network output oj can be represented with f(xj,β). The objective function of the BPN is to minimize RMSE or equivalently the sum of squared error (SSE) as follows:
(12)SSE(β)=∑j=1n(N(RCTEju)-f(xj,β))2.The Levenberg-Marquardt algorithm is an iterative procedure. At the beginning, the user should specify the initial values of the network parametersβ. It is common practice to set βT=(1,1,…,1). At each step, the parameter vector β is replaced by a new estimate β+δ, where δ=[Δw11h,…,ΔwPLh,Δθ1h,…,ΔθLh,Δw1o,…,ΔwLo,Δθo]. The network output becomes f(xj,β+δ); it is approximated by its linearization as
(13)f(xj,β+δ)≈f(xj,β)+Jjδ,
where
(14)Jj=∂f(xj,β)∂β
is the gradient vector of f with respect to β. Substituting (13) into (12),
(15)SSE(β+δ)≈∑j=1n(N(RCTEju)-f(xj,β)-Jjδ)2.
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 give
(16)(JTJ)δ=JT(N(RCTEju)-f(xj,β)),
where J is the Jacobian matrix containing the first derivative of network error with respect to the weights and biases. Equation (16) includes a set of linear equations that can be solved for δ.Finally, the three-layer BPN can be applied to estimate the remaining cycle time of a job. When a new job is released into the factory, the parameters associated with the new job are recorded. Then the new job is classified into a category, and the three-layer BPN of the category can be applied to estimate the remaining cycle time of the new job.Consider category 1 of the previous example. WhenμL=0.5, there are 16 jobs in this category. These jobs are split into two parts: the training data (the first 12 jobs) and the testing data (the remaining jobs). A three-layer BPN estimates the remaining cycle time of jobs in this category according to their six attributes with the following settings: single hidden layer,the number of neurons in the hidden layer: 12,convergence criterion: mean squared error (MSE)<10-4.The Neural Network Toolbox of MATLAB is used to implement the BPN approach. The BPN program code for the illustrative example is shown in Algorithm2. The estimation results are shown in Figure 3.Algorithm 2:The BPN code in MATLAB (category 1).
t
n_input  =    [0.500  0.300  0.500…;  0.396  0.574  0.811…;  0.900  0.900  0.900…;  0.722
0.811  0.900…;  0.700  0.567  0.633…;  0.279  0.314  0.314…]
tn_target  =  [0.341  0.383  0.598⋯0.100]
net    =    newff([0  1;  0  1;  0  1;  0  1;  0  1;  0  1],[12,  1],{“logsig”,  “logsig”},  “trainlm”);
net    =    init(net);
net.trainParam.show  =  10;
net.trainParam.lr  =  0.1;
net.trainParam.epochs  =  1000;
net.trainParam.goal  =  1e  −  4;
[
net,  tr]  =  train(net,  tn_input,  tn_target);
tn_output  =  sim(net,  tn_input)
te_input  =  [0.300  0.500  0.300…;  0.307  0.396  0.870…;  0.900  0.900  0.900…;
0.722  0.811  0.900…;  0.767  0.700  0.900…;  0.158  0.264  0.302…]
te_output  =  sim(net,  te_input)Figure 3
The estimation results by FCM-BPN.The following indexes evaluate the estimation accuracy of both training and testing data:mean absolute error (MAE) = 43 (hrs),mean absolute percentage error (MAPE) = 3.7%,RMSE = 62 (hrs).In contrast, if the remaining cycle times of these jobs are predicted in association with those of the other unclassified jobs, then the estimation performance will be mean absolute error (MAE) = 47 (hrs),mean absolute percentage error (MAPE) = 4.1%,root mean squared error (RMSE) = 70 (hrs)which is much poorer. Table 8 compares the performance levels of the gradient descent algorithm and the Levenberg-Marquardt algorithm.Table 8
Training algorithm performance (convergence criterion: MSE < 10−4).
Algorithm
Number of epochs
MAE
MAPE
RMSE
Gradient descent
10000
79
6.7%
98
Levenberg-Marquardt
<100
47
4.1%
70How does the accuracy of remaining cycle time estimation affect the performance of the dispatching rule? To answer this question, we add a noise factor to the remaining cycle time before feeding it to the dispatching rule as follows:(17)RCTEju+εju,εju∈N(0~200,0~50),
where εju follows the normal distribution. After a simulation study using FSMCT, the estimated remaining cycle time was compared with the actual value to evaluate the estimation accuracy. The relationship between the scheduling performance and the estimation accuracy is shown in Figure 4. Forty (40) replications of a simulation experiment were run with FSMCT. The scheduling performance and the estimation accuracy were evaluated with the average cycle time and RMSE, respectively. Obviously, the more accurately the remaining cycle time can be estimated, the better the schedule will perform. Therefore, the efforts made in this paper to improve the accuracy of remaining cycle time estimation make sense.Figure 4
The relationship between the scheduling performance and the estimation accuracy.


## 3.4. The Slack-Diversifying Fuzzy-Neural Rule

In traditional fluctuation smoothing (FS) there are two different formulation methods with two separate strengths [11]. One method is aimed at minimizing the mean cycle time with FSMCT as follows:
(18)SKju(FSMCT)=jλ-RCTEju.
The other method is aimed at minimizing the variance of cycle time with FSVCT as follows:
(19)SKju(FSVCT)=Rj-RCTEju.
Jobs with the smallest slack values are given high priority. These two rules and their variants have been proven to be very effective in shortening the cycle times of wafer fabrication factories [11, 14, 24, 27]. In the traditional FSMCT rule, RCTEju might be much greater than j/λ. As a result, the slack of a job is determined solely by RCTEju. To tackle this problem, both terms in the FSMCT rule are normalized as follows:
(20)Nor(RCTEju)=(RCTEju-min(RCTEju))(max(RCTEju)-min(RCTEju)),Nor(jλ)=(j/λ-min(j/λ))(max(j/λ)-min(j/λ))=(j/λ-1/λ)(n/λ-1/λ)=(j-1)(n-1).
After normalization, both terms now range from 0 to 1. Subsequently, to improve the responsiveness of the FSMCT rule, the division operator is applied instead of the traditional subtraction operator as follows:
(21)nonlinearFSMCTrule:SKju=Nor(j/λ)Nor(RCTEju).
However, this rule cannot be tailored to the wafer fabrication factory that is to be scheduled. To address this problem, the transition from a traditional FSMCT rule to its nonlinear form is analyzed as follows. The nonlinear form can be rewritten as(22)SKju=Nor(j/λ)Nor(RCTEju)=(j-1)/(n-1)(RCTEju-min(RCTEju))/(max(RCTEju)-min(RCTEju))=βγ·j-1RCTEju-min(RCTEju)=βλγ·(j/λ)-(1/λ)RCTEju-min(RCTEju)=βλγ(RCTEju-min(RCTEju))·(jλ-RCTEju+RCTEju-1λ)=(γ(RCTEju-min(RCTEju))βλ)-1·(jλ-RCTEju+(RCTEju-min(jλ))·1),where β=max(RCTEju)-min(RCTEju) and γ=n-1. Conversely, the linear form can also be rewritten as
(23)SKju=jλ-RCTEju=(γ(RCTEju-min(RCTEju))βλ)-0·(jλ-RCTEju+(RCTEju-min(jλ))·0).
These two formulas can be generalized into the following form:
(24)SKju=(γ(RCTEju-min(RCTEju))βλ)-ξ·(jλ-RCTEju+(RCTEju-min(jλ))·ζ).Wang et al. derived [13] the slack-diversifying nonlinear fluctuation smoothing rule by diversifying the slack in the 1f-TNFSVCT rule. In this study, we diversify the slack in the 2f-TNFSMCT rule as follows:
(25)SKju=(βλγ(RCTEju-min(RCTEju)))ξ·(jλ-RCTEju+ζ(RCTEju-1λ))=ajuξ(bju+cjuζ)=bjuajuξ+cjuζajuξ,
where
(26)aju=βλγ(RCTEju-min(RCTEju)),bju=jλ-RCTEju,cju=RCTEju-1λ.
There are many possible models that can form the combination of ξ and ζ. For example,
(27)(Linearmodel)ξ=ζ,(28)(Nonlinearmodel)ξ=ζk,k≥0,(29)(Logarithmicmodel)ξ=ln(1+ζ)ln2.
However, (25) is difficult to deal with. For this reason, the following polynomial fitting techniques are used to convert it into a more tractable form:
(30)xξ≅(1.15-0.11ξ)+(-0.07+0.57ξ)x.
The mean absolute percentage error (RMSE) of (30) is less than 9% when x≤10. RMSE will not be a serious problem, since it is the ξ value associated with the maximum σSKju that is to be found, not the SKju value. These polynomial fitting techniques are especially effective when x exceeds 1.Applying (30) to (25) gives
(31)SKju≅bju(1.15-0.11ξ+(-0.07+0.57ξ)aju)+cjuζ(1.15-0.11ξ+(-0.07+0.57ξ)aju)=(1.15bju-0.07ajubju)+(-0.11bju+0.57ajubju)ξ+(1.15cju-0.07ajucju)ζ+(-0.11cju+0.57ajucju)ξζ=dju+ejuξ+fjuζ+gjuξζ,
where
(32)dju=1.15bju-0.07ajubju,eju=-0.11bju+0.57ajubju,fju=1.15cju-0.07ajucju,gju=-0.11cju+0.57ajucju.
To diversify the slack, the standard deviation of the slack is to be maximized as follows:
(33)σSKju=∑j=1n(SKju-SKju¯)2n-1=1(n-1)∑j=1nSKju2-1n(∑j=1nSKju)2
which is equivalent to maximizing the following term:
(34)∑j=1nSKij2-1n(∑j=1nSKij)2=∑j=1n(dju+ejuξ+fjuζ+gjuξζ)2-1n(∑i=1n(dju+ejuξ+fjuζ+gjuξζ))2=w1+w2ζ+w3ξ+w4ζ2+w5ξ2+w6ζξ+w7ζξ2+w8ζ2ξ.
By taking the derivatives of (34) with respect to ζ and ξ and setting them equal to zero, we can obtain
(35)w2+2w4ζ+w6ξ+w7ξ2+2w8ζξ=0,(36)w3+2w5ξ+w6ζ+2w7ζξ+w8ζ2=0.
In addition, ζ and ξ must follow some model, for example, (27), (28), or (29). For example, if the nonlinear model in (28) is satisfied (k=2), then according to (35),
(37)w2+2w4ζ+w6ζ2+2w8ζ3+w7ζ4=0;ξ=ζ2,
while according to (36)
(38)w3+w6ζ+(2w5+w8)ζ2+2w7ζ3=0;ξ=ζ2.
The general solutions are too long to be presented here. Further, to guarantee a maximum, the second-order derivative must be negative as follows:
(39)2w4+2w8ξ≤0,2w5+2w7ζ≤0.
In addition, it is possible that some jobs have very large or very small slacks, either of which will distort the calculation of σSKju. For this reason, such jobs are excluded. Further, the values of ξ and ζ influence the range of SKju. For a fair comparison, σSKju should be divided by the range of SKju to determine the optimal values of ξ and ζ.An example is given in Table9 to illustrate the procedure mentioned above, in which α=244, β=195, and λ=1.18. Assume that the nonlinear model with k=2 is used for the relationship between ξ and ζ. The optimal solution is ξ* = 0.53 and ζ*=0.28 with the maximum σSKju equal to 9.21. In Figure 5 the results of the proposed methodology are compared with those of Wang et al.’s method. Obviously, the proposed methodology performed better at diversifying the slacks of jobs, which reduces the risk of misscheduling and promotes scheduling performance.Table 9
A demonstrative example.
j
R
j
RCTE
j
u
1
158
54
2
154
109
3
161
118
4
131
52
5
162
198
6
248
25
7
220
11
8
87
106
9
22
31
10
178
118
11
229
95
12
148
105
13
8
8
14
228
55
15
242
132
16
99
51
17
4
3
18
246
154
19
191
157
20
131
148Comparison of the results by the two methods.
(a)
(b)


## 4. Numerical Simulation and Results

To evaluate the effectiveness of the slack-diversifying fuzzy-neural rule for job dispatching in a wafer fabrication factory, simulated data were used to avoid disturbing the regular operations of the wafer fabrication factory. Simulation is a widely used technology that can assess the effectiveness of a scheduling policy, especially when the proposed policy and the current practice are very different. The actual production environment is dedicated to make real products and is not available for algorithm testing. Therefore, a real wafer fabrication factory located in Taichung Scientific Park of Taiwan with a monthly capacity of about 25,000 wafers was simulated. That factory’s real time scheduling systems input information very rapidly into its production management information systems (PROMIS). The simulation program has been validated and verified by comparing the actual cycle times with the simulated values and by analyzing the trace report, respectively. The wafer fabrication produces more than 10 types of memory products and has more than 500 workstations for performing single wafer or batch operations using 58 nm ~ 110 nm technologies. Each job released into the fabrication factory is assigned one of three priorities, that is, “normal,” “hot,” and “super hot.” Jobs with the highest priorities will be processed first. The large scale of operations with reentrant process flows makes job dispatching in the wafer fabrication factory a very tough task. Currently, the longest average cycle time exceeds three months with a variation of more than 300 hours. The factory managers seek better dispatching rules to replace first-in first-out (FIFO) and EDD, in order to shorten the average cycle times and ensure on time deliveries to customers. One hundred replications of the simulation were successively run. The time required for each simulation replication was about 30 minutes using a PC with Intel Dual CPU E2200 2.2 GHz and 1.99G RAM. A horizon of twenty-four months was simulated.To assess the effectiveness of the proposed methodology and to make comparisons with some existing approaches—FIFO, EDD, shortest remaining processing time (SRPT), CR, FSVCT, FSMCT, and a nonlinear fluctuation smoothing rule (NFS) [4]—all of these methods were applied to schedule the simulated wafer fabrication factory with the data of 1000 jobs; the collected data were then separated by product types and priorities. That is about the amount of work that can be achieved with 100% of the monthly capacity. Some cases produced insufficient data, so the present paper does not discuss those cases.The due date of each job was determined as follows. The FCM-BPN approach was applied to estimate the cycle time; the Levenberg-Marquardt algorithm, rather than the gradient descent algorithm, was applied to speed up the network convergence. Then, we added a constant allowance of three days to the estimated cycle time, that is,κ=72, to determine the internal due date as follows:
(40)DDj=CTEj+72.Jobs with the highest priorities are usually processed first. In FIFO, jobs were sequenced on each machine first by their priorities, then by their arrival times at the machine. In EDD, jobs were sequenced first by their priorities, then by their due dates. With SRPT, the remaining processing time of each job was calculated. Then, jobs were sequenced first by their priorities, then by their remaining processing times. With CR, jobs were sequenced first by their priorities, then by their critical ratios as follows:(41)CRju=(DDj-t)RPTj.
FSVCT and FSMCT consisted of two stages. The first stage scheduled jobs by FIFO; the remaining cycle times of all jobs were recorded and averaged at each step. The second stage applied FSVCT or FSMCT policy to schedule the jobs based on the average remaining cycle times obtained earlier. In other words, jobs were sequenced on each machine first by their priorities, and then by their slack values, which were determined by (18) and (19). In NFS, the slack of a job was defined as follows:
(42)SKju=(Rj-minRjmaxRj-minRj)ξ·((j-1)/(n-1))1-ξ((RCTEju-minRCTEju)/(maxRCTEju-minRCTEju)),
where ξ∈[01]. After a preliminary experiment, ξ was set to 0.8. In the proposed methodology, the nonlinear model with k=2 was used.Subsequently, the average cycle time and cycle time standard deviation of all cases were calculated to assess the scheduling performance. Comparisons for the average cycle time used FSMCT policy; comparisons for cycle time standard deviation used FSVCT. The results are summarized in Tables10 and 11. Both products A and B are dynamic random access memory products. The main difference between the two products is their storage capacities.(1)
For the average cycle time, the proposed methodology outperformed the baseline approach, FSMCT policy. The average advantage was about 17%.(2)
In addition, the proposed methodology surpassed the FSVCT policy in reducing cycle time standard deviation. The most obvious advantage was 74%.(3)
As expected, SRPT performed well in reducing the average cycle times, especially for product types with short cycle times (e.g., product A) but sometimes gave an exceedingly bad performance with respect to cycle time standard deviation. If the cycle time is long, the remaining cycle time will be much longer than the remaining processing time; this makes SRPT ineffective. SRPT is similar to FSMCT. Both try to make all jobs equally early or late.(4)
The performance of EDD was satisfactory for product types with short cycle times. If the cycle time is long, it is more likely to deviate from the prescribed internal due date, which makes EDD ineffective. That becomes serious if the percentage of the product type is high in the product mix (e.g., product type A). CR has similar problems.(5)
The Wilcoxon signed-rank test [39], a commonly used nonparametric statistical hypothesis test for comparisons of two related samples or repeated measurements on a single sample, was used in this study to assess whether population means differed or not. The results are summarized in Table 12. The null hypothesis Ha was rejected at α=0.025, which showed that the slack-diversifying fluctuation-smoothing rule was statistically superior to four existing approaches to reduce the average cycle time. With regard to cycle time standard deviation, the advantage of the slack-diversifying fluctuation-smoothing rule over SRPT and FSVCT was significant.(6)
The slack-diversifying 2f-TNFSMCT rule was also compared with the traditional one without slack diversification. Taking product type A with normal priority as an example, the comparison results are shown in Figure6. Obviously, the slack-diversifying rules dominated most of the traditional rules without slack diversification. According to these results, the treatments carried out in this study did indeed improve on the performance levels of the traditional policies.(7)
The proposed methodology can be decomposed into two parts. In the first part, the remaining cycle time of a new job is estimated using FCM-BPN, based on the historical data retrieved from PROMIS. In the experiment, it took less than 2 minutes to run all MATLAB programs and generate the remaining cycle time estimate. In contrast, Chen and Wang’s method used the gradient search algorithm, and it took about 15 minutes to train the BPN. In the second part, the slacks of all jobs are calculated according to the slack-diversifying 2f-TNFSMCT rule. That can be done almost instantly.Table 10
The performance of various approaches in the average cycle time.
Avg. cycle time (hrs)
A (normal)
A (hot)
A (super hot)
B (normal)
B (hot)
FIFO
1256
401
320
1278
457
EDD
1087
346
306
1433
478
SRPT
966
350
309
1737
483
CR
1143
356
301
1497
470
FSMCT
1401
405
320
1408
430
FSVCT
1046
385
317
1745
519
NFS
1456
407
321
1452
421
The proposed methodology
1134
295
278
1208
399Table 11
The performance of various approaches in cycle time standard deviation.
Cycle time std. dev. (hrs)
A (normal)
A (hot)
A (super hot)
B (normal)
B (hot)
FIFO
56
24
23
87
40
EDD
130
25
23
50
39
SRPT
246
32
23
106
30
CR
68
30
19
58
37
FSMCT
42
44
23
35
28
FSVCT
319
35
28
222
55
biNFS
87
49
19
44
47
The proposed methodology
89
26
15
55
20Table 12
Results of the Wilcoxon signed-rank test.
H
a
0
H
b
0
(the average cycle time)
(cycle time standard deviation)
FIFO
2.02**
0.40
EDD
1.48
1.21
SRPT
0.94
2.02**
CR
2.02**
0.54
FSMCT
2.02**
-
0.54
FSVCT
1.48
2.02**
biNFS
2.02**
0.94
*
P
<
0.05.
**P<0.025.
***P<0.01.Figure 6
Comparing the slack-diversifying rule with traditional rules without slack diversification.The slack-diversifying 2f-TNFSMCT rule has been proven effective in real world situations because it has been applied to schedule a wafer fabrication factory. However, it has not been possible to experiment with a large assortment of scheduling policies in that wafer fabrication factory. In the past, FIFO and EDD have been applied to this wafer fabrication factory, and their performance statistics have been analyzed. For these reasons, the performance level of the slack-diversifying 2f-TNFSMCT rule was compared with those of FIFO and EDD. For example, consider the case in which the majority of the factory’s capacity is occupied. Figures7 and 8 compare the performance levels of the three approaches. With regard to the average cycle time, the slack-diversifying 2f-TNFSMCT rule surpassed the two old rules. The performance of the slack-diversifying 2f-TNFSMCT rule was also comparable to that of FIFO for reducing cycle time standard deviation. Both results supported the practicability of the proposed slack-diversifying 2f-TNFSMCT rule.Figure 7
The performances of various rules in the average cycle time.Figure 8
The performances of various rules in cycle time standard deviation.


## 5. Conclusions and Directions for Future Research

To optimize job dispatching in a wafer fabrication factory is a challenging but important task. This study attempts to innovate by a slack-diversifying fuzzy neural rule that optimizes the average cycle time and cycle time standard deviation for job dispatching in a wafer fabrication factory.The proposed methodology applies soft computing techniques like fuzzy classification and artificial neural network prediction. An effective fuzzy-neural approach is applied to estimate the remaining cycle time of a job, which is empirically shown to ameliorate scheduling performance. This slack-diversifying fuzzy-neural rule is a modification of the 2f-TNFSMCT rule.After a simulation study, we observed the following phenomena:(1)
through improving the accuracy of remaining cycle time estimation, the performance of a scheduling rule can indeed be strengthened;(2)
optimizing the adjustable factors in the 2f-TNFSMCT rule appears to be an appropriate tool to enhance the scheduling performance of the rule.However, a complete assessment of the effectiveness and efficiency of the proposed methodology requires actual application in a real-world wafer fabrication factory. Future studies can optimize other rules in the same way, taking into account the great uncertainty inherent in wafer fabrication systems [40] and possibly using a more effective data mining approach [41].



---

*Source: 980984-2013-05-20.xml*