IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

353

An Algorithm of Multi-Subpopulation Parameters
With Hybrid Estimation of Distribution for
Semiconductor Scheduling With
Constrained Waiting Time
Hung-Kai Wang, Chen-Fu Chien, Member, IEEE, and Mitsuo Gen

Abstract—Scheduling for wafer fabrication of advanced technology nodes entails complicated constraints such as limited
waiting times. Focusing on real settings, this paper aims todevelop a novel genetic algorithm of multi-subpopulation parameters with hybrid estimation of distribution (MSPHEDA) tosolve the present problem effectively and efficiently. To estimate the validity of this approach, ten scenarios were simulatedon the basis of empirical data as the basis to compare theperformance of MSPHEDA and other heuristic methods for minimizing makespan and reducing the total exceeded limited waitingtime. The results have shown practical viability of the proposed
approach.
Index Terms—Hybrid estimation of distribution, limited waiting time constraint, multi-subpopulation, semiconductor manufacturing scheduling.

I. I NTRODUCTIONEMICONDUCTOR critical dimension shrinks rapidly
via advanced equipment driven by Moore’s Law [1].
Thus, semiconductor manufacturing continuously migrates
to advanced technology nodes that are increasingly
complicated [2]–[4]. Effective job scheduling to enhance
tool productivity and reduce production cycle time iscritical for maintaining competitive advantage [5], [6].
Mönch et al. [7] classified scheduling problems into six
types: batching problems, problems with auxiliary resources,
multiple orders per job problems, scheduling of cluster
tools, scheduling approaches for individual work areas, andS

Manuscript received October 11, 2014; revised April 12, 2015; accepted
May 25, 2015. Date of publication June 1, 2015; date of current version July 31, 2015. This work was supported in part by the Ministryof Science and Technology under Grant NSC 100-2628-E-007-017-MY3
and Grant MOST 103-2218-E-007-023, in part by the Semiconductor
Technologies
Empowerment
Partners
Consortium
under
Grant NSC 102-2622-E-007-013, in part by the Advanced Manufacturing
and Service Management Research Center under Grant 104N2074E1,
in part by the Taiwan Semiconductor Manufacturing Company, Hsinchu
Science Park, Taiwan under Grant 102A26, and in part by the Japan
Society of Promotion of Science: Grant-in-Aid for Scientific Research underGrant 24510219.0001.
H.-K. Wang and C.-F. Chien are with the Department of Industrial
Engineering and Engineering Management, National Tsing Hua University,
Hsinchu 30013, Taiwan (e-mail: cfchien@mx.nthu.edu.tw).
M. Gen is with both Fuzzy Logic Systems Institute and Research Institutefor Science and Technology, Tokyo University of Science, Japan.
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TSM.2015.2439054

job-shop problems. Indeed, the semiconductor manufacturingscheduling problem (SMSP) is an expanded model of theflexible job-shop scheduling problem (FJSP) that is more
complicated owing to additional constraints such as limited
waiting times [6]–[8]. However, few studies have addressed
the present problem for advanced wafer fabrication [7], [8].
As manufacturing contexts change rapidly, the schedulingsolution must be able to obtain a near-optimal solution withina short time. This study aims to propose an efficient algorithm that can obtain an effective solution with shorter limited
time than previous studies [6]–[8]. Conventional mixed-integer
linear programming (MILP) model cannot obtain a solution
quickly in practice [9]. Metaheuristics have been developed
to solve the scheduling problems. However, as the constraints increase, the efficiency of metaheuristics may beaffected [10].
Considering limited waiting time constraints in real settings,
this study aims to develop a novel genetic algorithm ofmulti-subpopulation parameters with hybrid estimation of distribution to solve the SMSP effectively. The present SMSP
is a single stage scheduling problem for unrelated parallelmachines with dedications. The proposed approach is validated
via simulations of ten scenarios on the basis of realistic ion
implantation process data and practical manufacturing situations including those related to process defects, new products,
ramping up, urgent jobs from customers, and abnormal eventsin advanced fabs. The results have shown practical viabilityof the proposed approach.
The rest of this paper is organized as follows: Section II
reviews related studies on the SMSP and evolutionary algorithms. Section III formulates an MINLP model of thepresent problem. Section IV presents the proposed algorithm.
Section V compares the proposed method with existing metaheuristics via simulations of different scenarios for validation.
Section VI concludes with a discussion of contributions andfuture research directions.
II. L ITERATURE R EVIEW
A. SMSP
Mathematical programming models have been developed
to solve scheduling problems. For instance, Jia et al. [11]proposed a slack-based mixed-integer programming modelc 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
0894-6507 
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

354

IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

for the scheduling of reentrant batch-processing machinewith incompatible job families in wafer fab. However,their study considered only up to eight job families.
Klemmt and Mönch [12] addressed the SMSP as a flow-shop
scheduling problem (FSP) with waiting time constraints andproposed a mixed-integer programming model with a scheduling and decomposition approach. They considered the time
constraints across several stages that can be nested. However,
it is hard to estimate the waiting time between different
stages and thus accurately handle the SMSP between several stages for modern Giga-fabs. For a medium-size problem,
Aguirre et al. [13] developed an MILP model for wet-etch
operation scheduling that can achieve near-optimal solutions
within reasonable time. Yan et al. [14] used a branch-and-cut
method with convex hull analyses for two-phase lithographymachine scheduling that can obtain near-optimal solutions
within reasonable time. Jung et al. [15] integrated an MILP
model into a real-time dispatcher for diffusion process scheduling with time window constraint. Since it will take several
hours for MILP model to obtain an optimal solution for thefull size problem, they firstly assigned a small number of runs
to the tools to reduce the problem scope and then extended thesolutions during several iterations. However, this study only
compared the proposed MILP with historical dispatching rule,
without comparing with other heuristics that can schedule all
the jobs at one time.
Simulation has been employed to determine schedulingrules. Qiao et al. [16] simulated a minifab and employeda hierarchical colored timed Petri net and an extended genetic
algorithm for scheduling. However, simulation approaches
have limitations in real settings of modern Giga-fabs.
Evolutionary algorithms [17] have been developed toaddress the large-size SMSP. For example, Wu and Chien [18]developed a GA for the scheduling of semiconductor final testing and implement it in practice. Song et al. [19] developed
an ant colony optimization to minimize the total unsupported demand and machine conversion time of a bottleneck station in semiconductor assembly and testing facilities.
Driessel and Mönch [20] considered the constraints of parallelmachines, sequence-dependent setup times, and the precedence and ready time of jobs, and proposed a variableneighborhood search approach to minimize the total weighted
tardiness. Chien et al. [21] developed an evolutionary approachfor the rehabilitation patient scheduling that is a hybrid shop.
Dauzere-Peres and Mönch [22] used a random key genetic
algorithm to solve a single-machine batch scheduling for
diffusion.
Attar et al. [23] solved a multi-objective hybrid flexible flow shop (HFFS) problem considering limited waitingtimes and machine setup time constraints. However, theirstudy assumed no machine breakdowns and maintenance andnegligible transportation time between stages and unlimited
intermediate storages that are not realistic in a practice.
Moreover, unlike HFFS, the SMSP is an FJSP in which eachjob has different routes with reentrant processes. Indeed, few
studies have addressed the SMSP with practical constraints
such as limited waiting time constraints. Chien and Chen [8]developed a batch sequencing genetic algorithm to solve thescheduling of oxide-nitride-oxide (ONO) stacked film fabrication processes, with consideration of waiting time constraints,
frequency-based setups, and capacity preoccupation.
B. EDA
Unlike the crossover and mutation in conventional GA, the
EDA [24] produces a population over multiple generations
to obtain improved solutions, via acquiring the probability
distribution from current elitist solutions and sampling new
solutions based on the probability distribution. After sampling,
the probability distribution evolves with previous probability
distribution and current elitist solutions. According to variable
dependency, the EDA can be classified as univariate, bivariate,
and multivariate distribution models. The population-based
incremental learning [25] and univariate marginal distribution
algorithm [26] belong to the univariate distribution model.
Mutual information maximization for input clustering [27]and bivariate marginal distribution algorithm [28] belong
to the bivariate distribution model. Bayesian optimization algorithm [29] belongs to the multivariate distribution
model.
Since EDA adapts and evolves from improved solutions
and probability distribution, EDA is more statistically valid
than conventional GA to obtain superior solutions [30], [31].
However, the EDA has limited exploitation ability [32]. Thus,
a number of studies have integrated EDA with local searchmethods to enhance its exploitation ability. For example,
Wang et al. [32] divided one population into two subpopulations to perform a local search independently, and combined
two subpopulations when no improvement occurred for specific generations. Gao et al. [33] used two types of bottleneck
shifting methods for neighborhood search. Gao et al. [34]proposed the variable neighborhood descent algorithm for
the FJSP. The EDA has been applied to various problems such as FJSP [30], FSP [35], single-machine schedulingproblem [36], traveling salesman problem [37], and multiobjective resource-constrained project scheduling [38]. However,few studies have employed the EDA to address the SMSP.
III. M IXED -I NTEGER N ONLINEAR P ROGRAMMING M ODEL
A. Characteristics of the SMSPAn MINLP model is formulated for the present problem
with the following constraints:
1) Limited Waiting Time Constraints: The advanced wafer
fabrication consists of processing stages with limited waiting time constraints. For example, the nitride layer should
be deposited within a strictly defined waiting time after thepolysilicon oxide layer is deposited since the surfaces ofpolysilicon oxide and nitride layers are unstable [7].
2) Machine Status and Recovery Time Constraints: The
machine status changes rapidly due to scheduled preventive maintenance or abnormal shutdown, while equipment
engineers will repair and recover it for rescheduling.
3) Different Processing Times for Jobs in Different
Machines: One job contains at most 25 wafers in one lot witha specific recipe by a machine containing the conditions suchas temperature, pressure, humidity, and the chemical elements.

WANG et al.: MSPHEDA FOR SEMICONDUCTOR SCHEDULING WITH CONSTRAINED WAITING TIMETABLE IDATA F ORM OF THE SMSP

Parameters:ntotal number of jobsmtotal number of machinespriority of job i
ωiprocessing time of job i in machine jpijsetup time between job i and job ksikarrival time of job iairecovery time of machine j
tjr
tiw
P
w
M
Ej

Different jobs fabricated with the corresponding recipes indifferent machines have different processing times.
4) Different Setup Time Between Jobs: The machine
requires a setup time between the jobs with different recipes.
5) Different Arrival Time: The jobs being processed in
the previous stage that will arrive in the near future and thejobs queued in this stage are considered in the forthcoming
schedule.
Table I lists the data of an example in fab. One job represents a number of wafers in a lot in the front opening
unified pod (FOUP). A FOUP may contain different numbersof wafers up to the full capacity of 25 wafers. The processing time of each job is affected by the number of wafers, therecipes, and the machines. The more wafers in a job, the longer
processing time required. Each job is conducted by a recipe.
The weighting of a job is assigned based on its importance for
production control. The arrival time indicates the job arriving
time at this stage.
B. MINLP Formulation
Keha et al. [39] modeled a single machine scheduling problem and Unlu and Mason [40] modeled a parallel machinescheduling problems as both sequence position based and timeindexed based. The sequence position based formulation withbinary decision variables xijk that are 1 if job j is on position
i on machine k, to avoid the nonlinearity with two multiplied
decision variable. This study considers real constraints suchas job arrival time, machine recovery time, and limited waiting time, in which the time indexed based formulation withtwo multiplied decision variable, job weight and job completion time can handle these constraints directly. To estimate
the validity of the proposed approach, an MINLP model wasformulated considering the constraints of operation allocation,
process time, machine recovery time, machine capacity, jobarrival time, limited waiting time, and setup time to optimize
the allocation and sequencing of jobs. The indices, parameters,
and decision variables are listed as follows:
Indices:i, kindex of jobs, i, k = 1, 2, . . . , njindex of machines, j = 1, 2, . . . , m355

limited waiting time of job i
penalty function when exceeding the limited waitingtime
coefficient of penalty function
a large positive integer number
the set of jobs that can be processed in machine j,i.e., the available jobs for machine j.

Decision variables:xij
xij = 1 if job i is processed in machine j; otherwise,
xij = 0process beginning time of job i in machine j
tijb
tijc
process completion time of job i in machine jzikj = 1 if job i precedes job k, which are both
zikj
processed in machine j; otherwise, zikj = 0
total time of job i exceeding tiw
Tiw
The objective function (1) is to minimize the weighted
sum of completion time for all the jobs and the penalty of
total exceeded limited waiting time for all jobs. Since different jobs have different priority, the domain experts set
the objective as job weight multiplied by job completiontime to ensure early processing of urgent jobs. A numberof multiobjective approaches such as Non-dominated Sorting
GA-II (NSGA-II) [41] and NSGA-III [42] can be employed.
Based on realistic needs, this study combined two objectives
into a single objective with a penalty P when limited waiting time constraints are violated [43]. In practice, minimizing
the total exceeded time of the limited waiting time is more
important than minimizing the weighted sum of completion
time. However, the fabs need to follow a schedule, even whensome jobs violate limited waiting times. This study models
the penalty as soft constraint to avoid infeasible solutions.
Objective function:
MinimizeZ =

n 
m


ωi tijc xij + P(1)

i=1 j=1

subject to1) Operation Allocation Constraints:m

xij = 1,

∀i

(2)

∀i, j

(3)

j=1

xij ∈ {0, 1},

where xij is a binary variable that assigns each job i to
machine j, where each job i will be processed in only one
machine j.
2) Process Time Constraints:
tijb + pij xij = tijc ,tijb , tijc ≥ 0,∀i ∈ Ej

(4)

∀i, j

(5)

356

IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

Equation (4) is the process time constraint that determines
the process beginning time of each job i, where its process
beginning time plus processing time should equal to its completion time in each machine. Equation (5) restricts the process
beginning time and completion time of each job i to greaterthan or equal to zero (i.e., nonnegative restriction).
3) Machine Recovery Time Constraints:
tijb + M 1 − xij ≥ tjr ,∀i ∈ Ej(6)

Equation (6) restricts the process beginning time of eachjob i in each machine j to be greater than or equal to therecovery time of machine j. In other words, one job can beprocessed in the machines that are recovered.
4) Machine Capacity Constraints:btijb + pij xij − M 1 − zikj ≤ tkj, ∀i, k ∈ Ej , i = kb
+ pkj xkj − M 1 − zkij ≤ tijb , ∀i, k ∈ Ej , i = ktkj
zikj + zkij = xij ∗ xkj , ∀i, k ∈ Ej , i = k

(8)
(9)

zikj ∈ {0, 1},(10)∀i, k, j(7)Equations (7) to (10) define all the preceding constraints forjob i and job k. In other words, one machine can process only
one job at a time.
5) Job Arrival Time Constraints:
tijb + M 1 − xij ≥ ai ,∀i ∈ Ej

P=w

Tiw

(13)

i=1

The total exceeded time of job i equals its completion timein machine j minus the limited waiting time of job i; otherwise,
the total exceeded time of job i equals zero, as in (12), which
is a soft constraint. If (12) is violated, the objective value will
have a penalty to minimize the total exceeded time, as in (13).
The coefficient of penalty function w can be determined basedon realistic needs.
7) Setup Time Constraints:b,tijb + pij xij + sik − M 1 − zikj ≤ tkjbb
tkj + pkj xkj + ski − M 1 − zkij ≤ tij ,Procedure of the MSPHEDA.

Fig. 2.

Chromosome representation.

(11)

A job can be processed in a machine only when the job
has arrived. Equation (11) restricts the beginning time of eachjob i in machine j to be greater or equal to the arrival time ofjob i.
6) Limited Waiting Time Constraints:(12)
Tiw = max tijc − tiw , 0 , ∀i ∈ Ejn
Fig. 1.

∀i, k ∈ Ej , i = k

(14)∀i, k ∈ Ej , i = k

(15)

The setup time is different between different jobs. A job
can be processed only after the machine is set up.
Equations (14) and (15) show that the process completion timeof job i plus the setup time between job i and subsequent job kshould be less than or equal to the beginning time of job i ineach machine j.

IV. M ULTI -S UBPOPULATION PARAMETERS W ITH H YBRIDE STIMATION OF D ISTRIBUTION A LGORITHM
To solve the SMSP effectively and efficiently, the proposed
algorithm of multi-subpopulation parameters with hybrid
estimation of distribution (MSPHEDA) designed chromosome representation, a decoding procedure, and a univariate
EDA model combining cooperative subpopulation with parameter settings and a local search mechanism as shown in Fig. 1.
A. Chromosome Representation
The designed chromosome consists of two parts, in whicha random key sequencing [44] is employed in both parts. Part I
represents the machine assignment, where the length is thenumber of jobs, and each gene represents each job on thesame site. For each gene, an integer number from one tototal number of machines m is used to assign each job to
its available machines. The relationship between the chromosome and gene is shown in Fig. 2. However, each job is notavailable for every machine. To avoid infeasible solutions ofmachine assignment, the integer number does not represent
the number of machines to be assigned, and is thus normalized between zero and one for decoding. The reason for using

WANG et al.: MSPHEDA FOR SEMICONDUCTOR SCHEDULING WITH CONSTRAINED WAITING TIME

Fig. 4.

Fig. 3.

Example of machine assignment.

an integer number instead of a random number between zeroand one is because the appearance time of every integer in the
EDA probability model is calculable. Part II represents the job
sequence, where the length of the chromosome is the numberof jobs, and each gene represents the sequence processing inthe same machine. Each gene is produced by a random number
from zero to one that determines the job sequence according
to the genes in Part II.
B. Initialize the Probability Matrix and Initial Population
All genes in Part I are independent and thus univariate
EDA model can be used to save computational time. An
element Pij (t) of the probability matrix P(t) represents theprobability of job i assigned to machine j. Equation (16)initializes the probability matrix P(t), representing the uniform distribution of job assignments. The initial population is
generated as follows:
Pij (1) = 1/m,

∀i, j

(16)

C. Decoding Procedure and Evaluation
Each gene is decoded by assigning each job to its available
machines. For example, there are five jobs and ten machines inTable I. For Job 1, there are six available machines; the probability of assigning Job 1 to each machine is 1/6 = 0.1667.
Equation (17) normalizes the gene number gi in part I of thechromosome to a normalized number Nori between zero andone to assign its available machines. In this example, the total
number of machines m is ten. If the gene number of Job 1 isthree, it can be normalized to 0.272 by (17) to assign Job 1
to Machine 4 through roulette wheel selection in Fig. 3.
Nori = gi /(m + 1),

∀i

(17)

When the machine assignment is complete, the genes inpart II are used to prioritize the processing sequence of the jobs
in the same machine. After the machine assignment and job
sequencing, the jobs are arranged into a timetable to determine
the process beginning time for all the jobs. The following
constraints are considered to ensure that the processing timeof each job is feasible.

357

Example of inserting jobs into a previous idle time interval.

1) Machine Recovery Time: All jobs can be arranged only
after the machine available time, as shown in (6).
2) Arrival Time of Each Job: All jobs are arranged after
their arrival time; if no jobs arrive at a specific time, themachines will be idle to wait for jobs, as shown in (11).
3) Setup Time Between Jobs: All jobs can be processed
only after the machines are set up. Therefore, one job will bearranged after the completion time of its previous job plus the
setup time between them, as shown in (14) and (15).
4) Inserting Jobs Into the Idle Time Interval Between
Previous Jobs: As illustrated in Fig. 4, assume that a longsetup time exists between Jobs 1 and 2. If the jobs after Job 2
can be inserted in this interval, the machine idle time can be
reduced and the machine utilization can be enhanced. All theinserting conditions of all the jobs in all the machines are considered to complete the decoding procedure. Job 3 cannot be
inserted because its process time is longer than this interval,
while Job 4 can be inserted since the summation of the setup
time between Jobs 1 and 4, its process time, and the setup
time are shorter than this interval.
D. Parameter Settings in Each Subpopulation
The multi-subpopulation employs several subpopulations toprevent a single population from jumping into local optima.
These subpopulations evolve separately and are cooperated
with each other after certain generations [45].
For the present problem, efficiency is critical for implementing the algorithm in real settings. Thus, the proposed
MSPHEDA employs different parameters in each subpopulation to enhance efficiency. Numerous studies have shown thata high crossover rate and high elitist selection lead to rapid
convergence [10]. However, a high mutation rate and highrank selection cause the solution to diversify and converge
slowly [25]. The proposed approach can balance the global
search ability and local search ability in each subpopulation.
For Subpopulation 1, a high crossover rate and low mutation
rate with high elitist selection and low rank selection are used.
For Subpopulation 2, a middle crossover rate and middle mutation rate with middle elitist selection and middle rank selection
are used. For Subpopulation 3, a low crossover rate and highmutation rate with low elitist selection and high rank selection
are used.
E. Select the Superior Subpopulation
The superior subpopulation is firstly selected to calculate
the superior appearance matrix A(t) as in (18) to calculate the358IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

each part of the chromosome to assist the current solution
jumping out of local optima to enhance solution quality.
I. Produce Superior Subpopulation for the Next Generation

Fig. 5.

Two-cut points crossover in part II of the chromosomes.

appearance time of each gene from the superior population.
SP denotes the total number of superior populations. δijs is anindicator vector of superior population s.
1  s
δij ,
SP
SP

Aij =

∀i, j

s=1

where
1,
δijs =
0,

if job i assigns to position j of total machines
otherwise(18)

F. Update the Probability Matrix
Equation (19) is used to update the probability matrix forthe next generation P(t + 1) that is the weighted sum of theprobability matrix P(t) and superior appearance matrix A(t),where α is the weight of Aij (t) for the learning rate of thesuperior population.
Pij (t + 1) = (1 − α)Pij (t) + αAij (t),
where
0≤α≤1

∀i, j,(19)

G. Create New Part I and Part II Offspring
The EDA creates new offspring by sampling from the probability matrix. For each gene in part I of the chromosome, theroulette wheel is used to sample an integer number from probability matrix P(t + 1). However, the univariate EDA modelcannot be used in part II of the chromosome since its genes
are not integer numbers. Two-cut points crossover to generate
new offspring that randomly select two chromosomes and two
cutting points to cut each chromosome into three segments,
and then exchange every gene in the middle segment between
the two chromosomes, as shown in Fig. 5.
H. Create New Offspring by Mutation
A probability model and updating mechanism perform
effectively in early generations. However, when the solution
falls into local optima in later generations, all the populations
may become identical. It is difficult to leave local optima only
by sampling from a probability matrix or by crossover. Thus,
two-cut point mutation is used by randomizing one gene in

After evaluating the new offspring, the next subpopulation
is produced from the offspring and current subpopulation. The
improved solutions are updated to the next generation. Percent
of elitist selection PES and percent of rank selection PRS areused. Elitist selection reserves the current optimal solutions
from the superior subpopulations, thereby ensuring the continuous improvement of the solution and also maintaining theconvergent speed. Rank selection selects the solutions fromnonelitists to maintain the diversity of the subpopulations.
That is, the MSPHEDA maintains the exploitation ability via elitist selection and the exploration ability via rank
selection.
J. Local Search
The proposed approach designed local search mechanism toenhance the EDA exploitation ability [32] as follows.
1) Beginning Time: When the solution does not improvefor consecutive TLS generations, all the chromosomes become
identical and thus cannot generate different new offspring.
Thus, a local search mechanism will be initiated to diversify
the chromosomes to possibly enhance the solutions.
2) Identifying the Critical Machine: The critical machineis the largest objective value among all the machines withthe largest loading. When switching one job from the criticalmachine to another feasible machine, the objective value Z
in (1) is likely reduced.
3) Switching the Machine Assignment and Job Sequence:
Select one job from the critical machine, assign it to anotherfeasible machine, and replace its job sequence with a random
number between zero and one.
4) Stopping Criteria: Retain the obtained superior solution
and select the next job in the critical machine until all the jobs
have been selected.
K. Exhaustive Local Search
In later generations, all the chromosomes may converge,causing the solutions into local optima. When the superior
solution is not improved for consecutive TELS generations,where TELS is greater than TLS indicating a local searchwithout improvement, an exhaustive local search should
be used.
The term “exhaustive” means to select one job at one timefrom the first job to the last job, assign it from Machine 1 tomachine m, and randomly generate its job sequence to assess
whether it obtains an enhanced solution. When an enhanced
solution is obtained, retain it and continue to assess othermachine assignments and job sequences until all the jobs havebeen searched. Although an exhaustive local search is time
consuming, it is an effective mechanism to find an improved
solution at a later stage. When all the solutions have been
searched, the current optimal solution will be replaced into
other subpopulations to ensure that all the subpopulations aresearched on the basis of the current optimal solution.

WANG et al.: MSPHEDA FOR SEMICONDUCTOR SCHEDULING WITH CONSTRAINED WAITING TIME

359

TABLE IIR ESULT OF LINGO AND THE MSPHEDA

TABLE IIIE XAMPLE DATA OF T EN J OBS AND T HREE M ACHINES

Fig. 6.
Example of coordinating two chromosomes between threesubpopulations.

L. Coordinate Chromosomes Between Subpopulations
In conventional metaheuristics, using only one population may converge into local optima. The proposed
MSPHEDA employs a multi-subpopulation Npop in which each
subpopulation evolves independently to increase diversity.
During specific generations TEX , the best top Nbest chromosomes in each subpopulation are selected. After sorting thesolutions of Nbest ∗ Npop chromosomes, the Nbest chromosomes
replace the chromosomes of each subpopulation. Fig. 6 illustrates an example with Npop equal to three and Nbest equal
to two. After incorporating the superior chromosomes, each
subpopulation evolves continuously for TEX generations.
V. VALIDATION
To validate the proposed MSPHEDA, several small-scale
problems and a large problem in a real setting were employed
for comparison. For the small problems, the LINGO software was used to solve the formulated MINLP models to
obtain the optimal solutions for comparison with those of
the proposed MSPHEDA. For the large problem in realsetting, the MSPHEDA was compared with other metaheuristics such as the GA, particle swarm optimization (PSO),
artificial bee colony (ABC), and EDA. Each simulation was
implemented using C++ program on Windows 7 X64 withIntel(R) Xeon(R) CPU E3-1230-V2 at 3.30GHz and an 8 GB
RAM processor.
A. Experimental Result of Small Problems
The following five small problems were designed for comparison: Problem 1 has seven jobs and three machines;
Problem 2 has ten jobs and three machines; Problem 3 has
ten jobs and five machines; Problem 4 has 15 jobs and threemachines; and Problem 5 has 15 jobs and five machines.
Equation (20) defines the optimality gap to denote the deviation percentage of the solution of the MSPHEDA from theoptimal solution of the MINLP obtained using LINGO:MSHEDA − MINLP
× 100%(20)
Optimality gap(%) =MINLP

Fig. 7.

Gantt chart of the example of ten jobs and three machines.

Domain experts determine that the setting of w equals to1000 in (13) for both MINLP and MSPHEGA. The remaining parameters of MSPHEDA are identical to those listed inTable IX. Table II lists the objective value, optimality gap, andCPU time of LINGO and the MSPHEDA. The CPU time of
LINGO is more than 8 hours in Problems 1, 2 and 3. The
MSPHEDA obtains the same solutions within 2 seconds. In
Problems 4 and 5, LINGO cannot obtain the optimal solution
within 1 day.
Table III shows the example data of ten jobs and threemachines. Without loss of generality, all ten jobs are 25 pieces,arriving at time 0, with the same limited waiting time of
120 minutes. M denotes a large processing time to preventthe jobs being processed in unavailable machines. The setup
time between two jobs with different recipes is 10 minutes.
For the example in Table III, both the MINLP and
MSPHEDA obtain the same optimal solution shown in
Table II. Fig. 7 illustrates the Gantt chart of the optimalscheduling result. Identical backgrounds of the jobs indicate identical recipes. Jobs 10, 9, 1, and 2 are processed in
Machine 1 because they all use Recipe A. Jobs 6, 4, 5, and 7
are processed in Machine 2 since they all use Recipe B.

360

IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

A 10 minute setup time exists between Job 8 and 3 inMachine 3 with different recipes.
The job with a higher weight divided by a shorter processing
time has a higher priority on the same machine. For example,
in machine 1, Job 10 (10/22 = 0.833) has a higher priority
than Job 9 (9/31 = 0.29); Job 9 has a higher priority thanJob 1 (1/6 = 0.167); and Job 1 has a higher priority than Job 2(2/21 = 0.095). The result is the same for Machines 2 and 3.
B. Experimental Design of Practical Problems
The proposed MSPHEDA was validated in a realistic projectsponsored by a leading semiconductor manufacturing company
in Taiwan. The data was collected from an ion implantation process. For confidentiality, all the continuous empirical
data such as processing times, setup times, and weights are
fitted and transformed with normal distribution N(μ, σ 2 ), andall the discrete empirical data such as recipes, number ofpieces, and available machines are fitted and transformed
with Bernoulli distribution. Without loss of generality, all the
specific terms regarding the machines, jobs, and recipes are
replaced by general terms. The following practical scenarios
were investigated.
1) Normal Situation: Scenario 1 denotes a normal case infab. Normal machine breakdown and maintenance rates are
10%, and these machines are repaired by equipment engineers
and returned to production within a specific time.
2) Process Defect Effect: Scenarios 2, 3, and 4 exhibit
a severe process defect in the previous stage, causing all thejobs to stack up and queue for a long time. The process defecteffect reduces the remaining waiting time and increases theconcentrated job arrival time, causing it to exceed that of
the normal situation, thereby increasing constraints for thescheduling problem.
Scenario 2 assumes that the jobs that arrive at this stage
have a shorter remaining waiting time than that of the normal
situation. In particular, the remaining waiting time of each job
is set as 60 minutes. Thus, each job should be processed within
60 minutes or it will become a defect.
Scenario 3 assumes that all the machines are prohibited producing any jobs until the defect is removed. The arrival time
has half the standard deviation of that in the normal situation,
causing the jobs to arrive in an increased concentrated time.
Scenario 4 combines Scenarios 2 and 3, while the remaining
waiting time of each job is 60 minutes when it arrives, andeach job has half the standard deviation of the arrival timecompared with the normal situation.
3) Machine Process Time Effect: Scenario 5 assumes thatthe variation of the machine process time is double that of
the normal situation. That is, advanced machines are released
for production, creating a diverse process time between theexisting and newly released machines. In other words, newmachine selection and installation affect the scheduling andproduction.
4) Urgent Jobs Effect: Scenario 6 enhances the mean priority of 30% of the jobs to be five times greater, reflecting
the surge in demand of critical customers. With the shortening product life cycle of consumer electronics, the proposedscheduling algorithm should address the needs of super-hot
lot to reduce the time-to-market and satisfy customer needs toenlarge the market share.
5) New Product Effect: Scenario 7 extends double recipes
for the ramping up of new technologies and new machines,
causing the setup times to increase between jobs. The proposedscheduling algorithm should reduce the setup time to enhance
the machine utilization and reduce the production cycle time.
6) New Fab Ramping Up Effect: Scenario 8 considers new
fab ramping up with the mean priority of 30% of the jobs tobe five times greater with double recipes for new technologies
and new machines compared with the normal situation. For
a new ramping fab, there is more super-hot lot to test newprocesses and new machines to enhance the yields of new
products. Therefore, more recipes exist for these new products
and new machines.
7) Abnormal Events Effect: Scenario 9 considers abnormal
events such as earthquakes and power and gas supply interruptions; this causes several machines to become unavailable,
and the unavailable rate increases to 20%.
8) Least Favorable Situation: Scenario 10 is the worst casethat has all the considerations of Scenarios 2 to 9 to address
all the possible constraints in fab to assess the performance ofthe proposed algorithm.
C. Experimental Result of Practical Problems
The ten aforementioned scenarios consider a practical problem of 500 jobs and 71 machines, in which more than
17.8 million decision variables are considered in the formulated MINLP model, subjected to more than 35.6 million
constraints. Thus, the scenarios cannot be solved to obtain the
optimal solution within a reasonable time. Since the heuristics
can obtain only the near optimal solutions of the large problem within a reasonable time, each scenario was repeatedly
performed ten times with the same testing data to compare the average performance of all the heuristic algorithms.
Since online scheduling requires a satisfactory solution within
3 minutes, each run should stop at 3 minutes.
To estimate the validity, three conventional metaheuristics
were compared: the GA [10], PSO [46], and ABC [47]. In thisstudy, GA uses the same chromosome representation, two-cut
points crossover, and two-cut points mutation as MSPHEDA in
Section IV. Two-cut points mutation randomizes one gene ineach part of the chromosome. PSO uses one particle to represent one solution in optimization problem. Many particles
exist simultaneously to simulate the swarm path by record thecurrent position, local best position, global best position andcurrent velocity to determine the velocities and the positions
for next iteration. ABC used the food source position to represent one solution. The onlooker bee will choose the new
food by the probability value calculated from fitness value.
The bees will search for new candidate food position basedon current position and a random number direction.
Indeed, the parameter setting of the GA is the same as
GA_H in Table VI, where PES is the percent of elitist selection,
PRS is the percent of rank selection, pc is the crossover rate
and pm is the mutation rate. The chosen numbers of theseWANG et al.: MSPHEDA FOR SEMICONDUCTOR SCHEDULING WITH CONSTRAINED WAITING TIMETABLE IV
M EAN OF O BJECTIVE FOR THE T EN S CENARIOSOF THE GA, PSO, AND ABC

TABLE VM EAN OF THE T OTAL E XCEEDED L IMITED WAITINGT IME OF THE GA, PSO, AND ABC

parameters are determined based on trial and error on a small
set of problem instances. The parameter setting of PSO is thatthe population size equals to 100, global acceleration constantequals to 2 and local acceleration constant equals to 2 [46].
The parameter setting of ABC is that initial population equals
to 50, population size equals to 100 and random number
between [−1, 1] [47]. Table IV shows that the GA improved
the PSO by 41% and the ABC by 43% for the average objectives of the ten scenarios. In addition, the GA reduced 46% of
the total exceeded limited waiting time in the PSO and 49% in
the ABC, as shown in Table V. Therefore, this study attempted
to improve the proposed algorithm from the GA to EDA as
follows.

361

TABLE VI
PARAMETER S ETTINGS OF GA S AND THE MSPGAThe GA fairly outperforms PSO and ABC. Thus,
this study improved GA with the proposed MSPGA(Multi-subpopulation Parameters with Genetic Algorithm)
to determine the different parameters. Subsequently, GA_H
was designed with a high crossover rate and low mutation rate. GA_M was designed with a medium crossover
rate and medium mutation rate. GA_L was designed witha low crossover rate and high mutation rate. In the proposed MSPGA, three subpopulations exist with the same
parameter settings as GA_H, GA_M, and GA_L. In GA_H,
GA_M, GA_L and MSPGA, the same chromosome representation, two-cut points crossover and two-cut points mutation
as MSPHEDA were used. Table VI shows the parameter settings of GAs and the MSPGA. Table VII shows the meanand standard deviation of the objective values of GAs and the
MSPGA, with the improvement rates of the MSPGA compared with those of GAs. In particular, the MSPGA has morethan 8% improvement in the objective mean than three GAs.
The MSPGA has 71% improvement in the standard deviation than three GAs. The results showed that the proposed
MSPGA can effectively improve GA with the same parameters. Table VIII shows the total exceeded limited waitingtime to analyze the constraint of limited waiting time. The
MSPGA obtains the least total exceeded limited waiting timethan GAs, in Scenarios 2, 3, 4 and 10. The MSPGA reduces
more than 16% average total exceeded limited waiting timeand more than 9% standard deviation of the total exceeded
waiting time.
To enhance the conventional EDA, the authors firstly proposed Hybrid Estimation of Distribution Algorithm with
Multiple Subpopulations (HEDA-MS) [48] that can outperform the EDA and GA. HEDA-MS used the same chromosome representation, initializing the probability matrix,
initial population, decoding procedure, evaluation, creating
new offspring and local search mechanism as MSPHEDA.
HEDA-MS used fuzzy logic controller (FLC) [49] to determine the crossover rate and mutation rate, in which FLC
increases the crossover rate and decreases mutation rate
when the fitness from crossover offspring performs well.
On the other hand, FLC decreases the crossover rate andincreases mutation rate when the fitness from mutation offspring performs well. This study improves the HEDA-MS
as MSPHEDA by using different parameter settings ineach subpopulation to replace FLC. Also, the MSPGA was362IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

TABLE VII
M EAN OF O BJECTIVE FOR T EN S CENARIOSB ETWEEN THE GA S AND MSPGATABLE VIIIM EAN OF THE T OTAL E XCEEDED L IMITED WAITING T IME FOR
T EN S CENARIOS B ETWEEN GA S AND THE MSPGA

extended to the MSPHGA (Multi-subpopulation Parameterswith Hybrid Genetic Algorithm) with two proposed local
search mechanisms. MSPHGA used the same chromosome
representation, two-cut point crossover, two-cut point mutationand three subpopulations as MSPGA. Thus, MSPHGA andMSPHEDA were compared with the same parameter settings.
Table IX shows different PES , PRS , pc , and pm settings
among HEDA-MS, MSPHGA, and MSPHEDA. In MSPHGA
and MSPHEDA, a high crossover rate with low elitist selection was set for Subpopulation 1, a medium crossover ratewith medium elitist selection was set for Subpopulation 2,
and a low crossover rate with high elitist selection was set for
Subpopulation 3. The proposed approach can maintain both
exploration and exploitation abilities in each subpopulation.
Table X shows the mean objectives of the HEDA-MS,
MSPHGA, and MSPHEDA among the ten scenarios. The
proposed MSPHEDA improves 4% of the mean and 51%of the standard deviation comparing with the HEDA-MS, on
average. The MSPHEDA also improves 31% of the mean and

TABLE IXPARAMETER S ETTINGS OF THE HEDA-MS, MSPHGA, AND MSPHEDA

TABLE XM EAN OF O BJECTIVE AND E XECUTED S OLUTIONS FOR T EN S CENARIOSA MONG THE MSPHGA, HEDA-MS, AND MSPHEDA

92% of the standard deviation compared with the MSPHGA,
on average. The proposed MSPHEDA can generate superior
solutions compared with the HEDA-MS and MSPHGA.
Table XI shows that the MSPHEDA can obtain 6% and 50%
less total exceeded limited waiting time than the HEDA-MS
and MSPHGA can, respectively. The proposed MSPHEDA has
robust performance than the HEDA-MS and MSPHGA.
Fig. 8 shows the box plots of the mean objective values ofthe HEDA-MS and MSPHEDA for Scenarios 1 to 10, indicating that the proposed MSPHEDA has smaller mean values andthe significantly smaller standard deviation, indicating crucial
stability in practice.
D. Experimental Result of Academic Instances
Five instances are designed for medium size problems.
Each instance contains 50 jobs assigning to 15 machines.

WANG et al.: MSPHEDA FOR SEMICONDUCTOR SCHEDULING WITH CONSTRAINED WAITING TIME

363

TABLE XIM EAN OF THE T OTAL E XCEEDED L IMITED WAITING T IME FOR T EN
S CENARIOS A MONG THE MSPHGA, HEDA-MS, AND MSPHEDA

Each data is generated from the same distribution of practical
data. Instance 1 is the normal case. Instance 2 generates themean of the machine process time double that of the normal
situation. Instance 3 assumes that the mean of the limited waiting time is 80% that of the normal situation. Instance 4 doubles
the number of recipes from the normal situation. Instance 5
combines all the considerations of Instance 1 to 4.
The parameter settings of all the algorithms are identical to the setting for practical problems in Section V-C.
Tables XII and XIII show the experimental results for the ninealgorithms. MSPHEDA performs the best among the others
for both the mean of the objectives and the mean of the total
exceeded limited waiting time. MSPHEDA improves 3% of
the mean and 23% of the standard deviation compared with
the HEDA-MS.
E. Discussion
The proposed approach designed the parameter settings canoutperform three GAs. The effectiveness of the proposed local
search mechanisms was validated. Tables VII and X showed
20% improvement from MSPGA to MSPHGA, on average.
Moreover, the MSPHGA can be enhanced to the
MSPHEDA that can obtain 31% improvement and 45% more
solutions comparing to the MSPHGA with the same parametersetting in Table X. Since the EDA samples new solutions on
the basis of a probability matrix, the CPU execution speed isfaster than the crossover in GA. Thus, the EDA embedded in
MSPHGA can obtain better solutions in a global search thanconventional GA.
Finally, the performance between the HEDA-MS andMSPHEDA was compared. The FLC in the HEDA-MS iseffective in obtaining improved solutions since it can automatically set pc and pm . However, it is time consuming toobtain suitable pc and pm . The proposed MSPHEDA canset different parameters of pc , pm , PES , and PRS in eachFig. 8.

Box plot of ten scenarios.

subpopulation to save the calculation time and thus effectively
improve the executed solutions by 3% on average. Table X and
Fig. 8 showed MSPHEDA obtains a 51% smaller standard
deviation of the objective values than HEDA-MS, indicating
critical robustness for online scheduling in real settings.
VI. C ONCLUSION
This study presented a novel algorithm called the
MSPHEDA for semiconductor manufacturing scheduling withconstrained waiting times. To estimate its validity, an experiment was designed to compare it with alternative approaches innumerous scenarios for the SMSP in real settings. An MINLP
model was formulated for the small-scale SMSP problems by
LINGO to obtain optimal solutions for validation. The results
showed that the proposed approach can obtain the optimal
solution for small-scale problems. Furthermore, ten common
scenarios in real settings were considered. The results showedthat the average mean and standard deviation of the proposed
algorithm were smaller than those of other heuristics suchas the PSO, ABC, GA, MSPGA, MSPHGA, and HEDA-MS.

364

IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015

TABLE XIIM EAN OF O BJECTIVE FOR F IVE I NSTANCESA MONG A LL THE N INE A LGORITHMSTABLE XIV
ACADEMIC I NSTANCES I

TABLE XIIIM EAN OF THE T OTAL E XCEEDED L IMITED WAITING T IME FOR F IVEI NSTANCES A MONG A LL THE N INE A LGORITHMS

The proposed algorithm has the shortest total exceeded limited
waiting time.
Future research can apply the proposed algorithm to othercomplicated batch production problems such as those related
to etching and furnaces for semiconductor manufacturing.
Because the waiting time of jobs between different process
stages is difficult to estimate in current Giga-fabs, most existing semiconductor manufacturing studies have focused on
individual work areas such as oxidation, deposition, diffusion,
photolithography, etching, ion implantation, and planarization, and have not covered the scheduling process loopthat has increasingly complicated constraints such as waiting times, frequency-based setups, and capacity preoccupation. For instance, Mason et al. [50] developed a modified
shifting bottleneck heuristic to minimize the total weighted
tardiness for a complex job shop scheduling problem ina fab. Furthermore, when the waiting time between two
consecutive stations can be estimated more precisely, the proposed approach can be extended to address multiple-station
scheduling problems for wafer fabrication such as ONO
stacked film processes.

A PPENDIX
Five instances are designed for the medium size problems
in Section V-D. Tables XIV and XV illustrate two of theinstances. All the 15 machines are available at time zero. The
setup time are set as 10 minutes between different recipes andno setup time between the same recipes.

WANG et al.: MSPHEDA FOR SEMICONDUCTOR SCHEDULING WITH CONSTRAINED WAITING TIMETABLE XVACADEMIC I NSTANCES IIR EFERENCES[1] G. E. Moore, “Cramming more components onto integrated circuits,”
Electronics, vol. 38, no. 8, pp. 114–117, 1965.
[2] R. C. Leachman, S. Ding, and C.-F. Chien, “Economic efficiency analysis of wafer fabrication,” IEEE Trans. Autom. Sci. Eng., vol. 4, no. 4,
pp. 501–512, Oct. 2007.

365[3] C.-F. Chien et al., “Modeling and analysis of semiconductor manufacturing in a shrinking world: Challenges and successes,” Eur. J. Ind. Eng.,
vol. 5, no. 3, pp. 254–271, 2011.
[4] C.-F. Chien, Y.-J. Chen, C.-Y. Hsu, and H.-K. Wang, “Overlay error
compensation using advanced process control with dynamically adjusted
proportional-integral R2R controller,” IEEE Trans. Autom. Sci. Eng.,
vol. 11, no. 2, pp. 473–484, Apr. 2014.
[5] C. Kuo, C.-F. Chien, and J. Chen, “Manufacturing intelligence to exploit
the value of production and tool data to reduce cycle time,” IEEE Trans.
Autom. Sci. Eng., vol. 8, no. 1, pp. 103–111, Jan. 2011.
[6] C.-F. Chien and C. Chen, “Using GA and CTPN for modeling the
optimization-based schedule generator of a generic production scheduling system,” Int. J. Prod. Res., vol. 45, no. 8, pp. 1763–1789, 2007.
[7] L. Mönch, J. W. Fowler, S. Dauzere-Peres, S. J. Mason, and O. Rose,“A survey of problems, solution techniques, and future challenges
in scheduling semiconductor manufacturing operations,” J. Schedul.,
vol. 14, no. 6, pp. 583–599, 2011.
[8] C.-F. Chien and C.-H. Chen, “A novel timetabling algorithm for a furnace process for semiconductor fabrication with constrained waiting andfrequency-based setups,” OR Spectr., vol. 29, no. 3, pp. 391–419, 2007.
[9] Y. Demir and S. K. Isleyen, “Evaluation of mathematical models forflexible job-shop scheduling problems,” Appl. Math. Model., vol. 37,
no. 3, pp. 977–988, 2013.
[10] M. Gen, R. Cheng, and L. Lin, Network Models and Optimization:
Multiobjective Genetic Algorithms. London, U.K.: Springer, 2008.
[11] W. Jia, Z. Jiang, and Y. Li, “Combined scheduling algorithm forre-entrant batch-processing machines in semiconductor wafer manufacturing,” Int. J. Prod. Res., vol. 53, no. 6, pp. 1–14, 2015.
[12] A. Klemmt and L. Mönch, “Scheduling jobs with time constraints
between consecutive process steps in semiconductor manufacturing,” inProc. Winter Simulat. Conf., Berlin, Germany, 2012, pp. 1–10.
[13] A. M. Aguirre, C. A. Mendez, and P. M. Castro, “A novel optimization method to automated wet-etch station scheduling in semiconductor manufacturing systems,” Comput. Chem. Eng., vol. 35, no. 12,
pp. 2960–2972, 2011.
[14] B. Yan, H. Y. Chen, P. B. Luh, S. Wang, and J. Chang, “Litho machinescheduling with convex hull analyses,” IEEE Trans. Autom. Sci. Eng.,
vol. 10, no. 4, pp. 928–937, Oct. 2013.
[15] C. Jung, D. Pabst, M. Ham, M. Stehli, and M. Rothe, “An effective
problem decomposition method for scheduling of diffusion processes
based on mixed integer linear programming,” IEEE Trans. Semicond.
Manuf., vol. 27, no. 3, pp. 357–363, Aug. 2014.
[16] F. Qiao, Y. M. Ma, L. Li, and H. X. Yu, “A Petri net and extended
genetic algorithm combined scheduling method for wafer fabrication,”
IEEE Trans. Autom. Sci. Eng., vol. 10, no. 1, pp. 197–204, Jan. 2013.
[17] X. Yu and M. Gen, Introduction to Evolutionary Algorithm. Berlin,
Germany: Springer, 2010, p. 430.
[18] J.-Z. Wu and C.-F. Chien, “Modeling semiconductor testing job scheduling and dynamic testing machine configuration,” Expert Syst. Appl.,
vol. 35, nos. 1–2, pp. 485–496, 2008.
[19] Y. Song, M. T. Zhang, J. Yi, L. Zhang, and L. Zheng, “Bottleneckstation scheduling in semiconductor assembly and test manufacturing
using ant colony optimization,” IEEE Trans. Autom. Sci. Eng., vol. 4,
no. 4, pp. 569–578, Oct. 2007.
[20] R. Driessel and L. Mönch, “Variable neighborhood search approaches
for scheduling jobs on parallel machines with sequence-dependent setup
times, precedence constraints, and ready times,” Comput. Ind. Eng.,
vol. 61, no. 2, pp. 336–345, 2011.
[21] C.-F. Chien, F. Tseng, and C. Chen, “An evolutionary approach to rehabilitation patient scheduling: A case study,” Eur. J. Oper. Res., vol. 189,
no. 3, pp. 1234–1253, 2008.
[22] S. Dauzere-Peres and L. Mönch, “Scheduling jobs on a single batch processing machine with incompatible job families and weighted number oftardy jobs objective,” Comput. Oper. Res., vol. 40, no. 5, pp. 1224–1233,
2013.
[23] S. F. Attar, M. Mohammadi, R. Tavakkoli-Moghaddam, and S. Yaghoubi,“Solving a new multi-objective hybrid flexible flowshop problem withlimited waiting times and machine-sequence-dependent set-up time constraints,” Int. J. Comput. Integr. Manuf., vol. 27, no. 5, pp. 450–469,
2014.
[24] P. Larrañaga and J. A. Lozano, Estimation of Distribution Algorithms:
A New Tool for Evolutionary Computation. New York, NY, USA:
Springer, 2002.
[25] S. Baluja, “Population-based incremental learning: A method for integrating genetic search based function optimization and competitive
learning,” School Comput. Sci., Carnegie Mellon Univ., Pittsburgh, PA,
USA, Tech. Rep. CMU-CS-94-163, Jun. 1994.

366

IEEE TRANSACTIONS ON SEMICONDUCTOR MANUFACTURING, VOL. 28, NO. 3, AUGUST 2015[26] H. Mühlenbein and G. Paaß, “From recombination of genes to theestimation of distributions I: Binary parameters,” in Parallel Problem
Solving From Nature—PPSN IV (LNCS 1141). Berlin, Germany:
Springer, 1996.
[27] J. S. De Bonet, C. L. Isbell, and P. Viola, “MIMIC: Finding optima by
estimating probability densities,” in Proc. Conf. Adv. Neural Inf. Process.
Syst, vol. 9. Denver, CO, USA, 1997, pp. 424–430.
[28] M. Pelikan and H. Mühlenbein, Advances in Soft Computing:Engineering Design and Manufacturing. London, U.K.: Springer, 1999.
[29] M. Pelikan, D. E. Goldberg, and E. Cantu-Paz, “BOA: The Bayesian
optimization algorithm,” in Proc. Genet. Evol. Comput., San Francisco,
CA, USA, 1999, pp. 525–532.
[30] J. Ceberio, E. Irurozki, A. Mendiburu, and J. A. Lozano, “A review on
estimation of distribution algorithms in permutation-based combinatorial
optimization problems,” Progr. Artif. Intell., vol. 1, no. 1, pp. 103–117,
2012.
[31] X. C. Hao, J.-Z. Wu, C.-F. Chien, and M. Gen, “The cooperative estimation of distribution algorithm: A novel approach for semiconductor final
test scheduling problems,” J. Intell. Manuf., vol. 25, no. 5, pp. 867–879,
2014.
[32] L. Wang, S. Y. Wang, Y. Xu, G. Zhou, and M. Liu, “A bi-population
based estimation of distribution algorithm for the flexible job-shop
scheduling problem,” Comput. Ind. Eng., vol. 62, no. 4, pp. 917–926,
2012.
[33] J. Gao, M. Gen, L. Y. Sun, and X. H. Zhao, “A hybrid of genetic
algorithm and bottleneck shifting for multiobjective flexible job shop
scheduling problems,” Comput. Ind. Eng., vol. 53, no. 1, pp. 149–162,
2007.
[34] J. Gao, M. Gen, L. Y. Sun, and X. H. Zhao, “A hybrid genetic andvariable neighborhood descent algorithm for flexible job shop schedulingproblems,” Comput. Oper. Res., vol. 35, no. 9, pp. 2892–2907, 2008.
[35] B. Jarboui, M. Eddaly, and P. Siarry, “An estimation of distribution algorithm for minimizing the total flowtime in permutation flowshop scheduling problems,” Comput. Oper. Res., vol. 36, no. 9, pp. 2638–2646,
2009.
[36] S. H. Chen, M. C. Chen, P. C. Chang, Q. F. Zhang, and Y. M. Chen,“Guidelines for developing effective estimation of distribution algorithms in solving single machine scheduling problems,” Expert Syst.
Appl., vol. 37, no. 9, pp. 6441–6451, 2010.
[37] V. Robles, P. D. Miguel, and P. Larrañaga, “Solving the traveling salesman
problem with EDAs,” Estimation of Distribution Algorithms: A New Tool
for Evolutionary Computation. New York, NY, USA: Springer, 2002.
[38] L. Wang, C. Fang, C.-D. Mu, and M. Liu, “A Pareto-archived estimationof-distribution algorithm for multiobjective resource-constrained projectscheduling problem,” IEEE Trans. Eng. Manage., vol. 60, no. 3,
pp. 617–626, Aug. 2013.
[39] A. B. Keha, K. Khowala, and J. W. Fowler, “Mixed integer programming
formulations for single machine scheduling problems,” Comput. Ind.
Eng., vol. 56, no. 1, pp. 357–365, 2009.
[40] Y. Unlu and S. J. Mason, “Evaluation of mixed integer programming
formulations for non-preemptive parallel machine scheduling problems,”
Comput. Ind. Eng., vol. 58, no. 4, pp. 785–800, 2010.
[41] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist
multiobjective genetic algorithm: NSGA-II,” IEEE Trans. Evol. Comput.,vol. 6, no. 2, pp. 182–197, Apr. 2002.
[42] K. Deb and H. Jain, “An evolutionary many-objective optimization
algorithm using reference-point-based nondominated sorting approach,part I: Solving problems with box constraints,” IEEE Trans. Evol.
Comput., vol. 18, no. 4, pp. 577–601, Aug. 2014.
[43] Ö. Yeniay, “Penalty function methods for constrained optimization withgenetic algorithms,” Math. Comput. Appl., vol. 10, no. 1, pp. 45–56,
2005.
[44] J. C. Bean, “Genetic algorithms and random keys for sequencing andoptimization,” ORSA J. Comput., vol. 6, no. 2, pp. 154–160, 1994.
[45] J. K. Cochran, S. M. Horng, and J. W. Fowler, “A multi-population
genetic algorithm to solve multi-objective scheduling problems for parallel machines,” Comput. Oper. Res., vol. 30, no. 7, pp. 1087–1102,
2003.
[46] J. Kennedy and R. Eberhart, “Particle swarm optimization,” in Proc.
IEEE Int. Conf. Neural Netw., vol. 4. Piscataway, NJ, USA, 1995,
pp. 1942–1948.
[47] D. Karaboga and B. Basturk, “A powerful and efficient algorithm
for numerical function optimization: Artificial bee colony (ABC)algorithm,” J. Glob. Optim., vol. 39, no. 3, pp. 459–471, 2007.

[48] H.-K. Wang, C.-F. Chien, and M. Gen, “Hybrid estimation of distribution
algorithm with multiple subpopulations for semiconductor manufacturing scheduling problem with limited waiting-time constraint,” in Proc.
IEEE Int. Conf. Autom. Sci. Eng., Taipei, Taiwan, 2014, pp. 101–106.
[49] Y. Yun and M. Gen, “Performance analysis of adaptive genetic algorithms with fuzzy logic and heuristics,” Fuzzy Optim. Decis. Making,vol. 2, no. 2, pp. 161–175, 2003.
[50] S. J. Mason, J. W. Fowler, and W. M. Carlyle, “A modified shiftingbottleneck heuristic for minimizing total weighted tardiness in complex
job shops,” J. Schedul., vol. 5, no. 3, pp. 247–262, 2002.
Hung-Kai Wang received the M.S. degree inindustrial engineering and engineering management from National Tsing Hua University (NTHU),
Taiwan, in 2010, where he is currently pursuing the Ph.D. degree with the IEEM, NTHU. His
research interests include advanced process control, genetic algorithms, modeling and analysis forsemiconductor manufacturing, and intelligent manufacturing. His research works appear in the IEEET RANSACTIONS ON AUTOMATION S CIENCE AND
E NGINEERING. He was a recipient of the Best
Thesis Award.
Chen-Fu Chien (M’03) received the B.S. degree(with Phi Tao Phi Hons.) in industrial engineering and electrical engineering from National
Tsing Hua University (NTHU), Taiwan, in 1990,
the M.S. and Ph.D. degrees from the Universityof Wisconsin-Madison, USA, in 1994 and 1996,
respectively. He was a Fulbright Scholar with the
University of California, Berkeley, from 2002 to 2003
and received PCMPCL Training with Harvard
Business School in 2007. He is a Tsing Hua Chair
Professor and the Director of the NTHU-Taiwan
Semiconductor Manufacturing Company (TSMC) Center for Manufacturing
Excellence. He is the Principal Investigator for the Semiconductor Technologies
Empowerment Partners Consortium funded by the Ministry of Science andTechnology, Taiwan. From 2005 to 2008, he was on-leave as the Deputy Directorof Industrial Engineering Division in TSMC. He holds eight USA invention
patents and published three books, over 120 journal papers, and a numberof case studies in Harvard Business School. His research efforts center on
decision analysis, big data analytics, modeling and analysis for semiconductor manufacturing, and advanced intelligent manufacturing systems. He
was a recipient of the National Quality Award, the Distinguished ResearchAward and the Tier 1 Principal Investigator Award (Top 3%) from NSC, theDistinguished University–Industry Collaborative Research Award from the
Ministry of Education, the University Industrial Contribution Award from
the Ministry of Economic Affairs, the Distinguished University–Industry
Collaborative Research Award and the Distinguished Young Faculty ResearchAward from NTHU, the Distinguished Young Industrial Engineer Award, theBest IE Paper Award, and the IE Award from CIIE, the Best Engineering
Paper Award and the Distinguished Engineering Professor by the Chinese
Institute of Engineers in Taiwan, and the 2012 Best Paper Award from theIEEE T RANSACTIONS ON AUTOMATION S CIENCE AND E NGINEERING.
Mitsuo Gen received the B.E., M.E., and Ph.D.
degrees in electronic engineering from Kogakuin
University in 1969, 1971, and 1975, respectively, the Ph.D. degree in informatics from
Kyoto University, Japan, in 2006. He is a SeniorResearch Scientist with Fuzzy Logic Systems
Institute and a Visiting Professor at Tokyo Universityof Science, Japan. He was a Professor with
the Graduate School of Information, Productionand Systems, Waseda University, Japan. He wasa Visiting Professor with the IEEM, National Tsing
Hua University, Taiwan, from 2012 to 2013, a Hanyang Chair Professor with
Hanyang University, Korea, from 2011 to 2012, and a Visiting Professor with
the University of California, Berkeley, from 1999 to 2000 and Texas A&M
University in 2000. His research interest includes evolutionary algorithms,
manufacturing scheduling, logistics network, and decision making. He has coauthored five books such as Introduction to Evolutionary Algorithms in 2010
and Network Models and Optimization: Multiobjective Genetic Algorithm
Approach (Springer) in 2008. He is an Area Editor of Computers & Industrial
Engineering and an Editorial Board Member of several international journals.