Bottleneck Identification and Prediction of Wafer Fabrication Systems in Transient States

ZHOU Bing－h(huán)ai(周炳海)，YIN Meng(殷萌)，SUN Chao(孫超)

School of Mechanical Engineering，Tongji University，Shanghai201804，China

GU Ping(古平)*，LUO Xin(羅辛)，YANG Rui－long(楊瑞龍)，ZHANG Cheng(張程)

Department of Computer Science and Technology，Chongqing University，Chongqing 400044，China

LU Zhi－qiang(陸志強)*，WUWen(吳文)，HAN Xiao－le(韓笑樂)，NIRAVONG Juliane

School of Mechanical Engineering，Tongji University，Shanghai201804，China

Bottleneck Identification and Prediction of Wafer Fabrication Systems in Transient States

ZHOU Bing－h(huán)ai(周炳海)*，YIN Meng(殷萌)，SUN Chao(孫超)

School of Mechanical Engineering，Tongji University，Shanghai201804，China

According to theory of constraints(TOCs)，the performance of a com plex manufacturing system，such as a wafer fabrication system，ismainly determ ined by its bottleneck machine．A method of the identification and prediction of the bottleneck machine was proposed in transient states of a system．Firstly，the bottleneck index was formulated based on the workloads and the variability in wafer fabrication system s．Second ly，main factors causing the variability and their influenceson the bottleneck machine in transient states of the system were analyzed and discussed．An effective bottleneck identification and prediction model was presented，which incorporated the variability and queuing theory，and took machine breakdowns and setups into considerations．Finally，the proposed bottleneck prediction method was verified by simulation experiments．Results indicate that the proposed bottleneck prediction method is feasible and effective．

wafer fabrications;transient states;bottleneck;prediction

Introduction

Due to increases of custom ization demands in the semiconductor industry，product types and changing over frequencies also increase.Thesemay lead to the increase of the changing over time in transient states.Themachine breakdowns and setups are the main causes in transient states of wafer fabrication systems.Machine repair time and setup time affect themean and the standard deviation of effective processing time of wafers.They will cause random fluctuations in the system，and result in the bottleneck shifting.Bottleneck identifications are vital，and will do help for the efficiency improvement of wafer fabrication systems.

The identification of a bottleneckmachine inmanufacturing systems has been paid much attention to in both academ ia and industry［1－2］.Law and Kelton［3］defined amachine which had the longest waiting time in queue as the system bottleneck. Zhou and Rose［4］identified the bottleneck resources by considering the utilization rate of machines in a system.Chang et al.［5］introduced a bottleneck identification method based on machines'sensitive coefficients.Chiang et al.［6］proposed a machine－h(huán)ungry/blocked－based bottleneck identification method.But thismethod was only applicable to the flow shop systems.Li et al.［7］put forward a concept of the instantaneous bottleneck，and analyzed the efficiency of potential bottlenecks of the system by mathematical and simulation methods. However，the exact throughput analysis only existed for simple two stations and one buffer system with finite buffer capacities or with no buffer.Wang etal.［8］combined a scheduling policy with the bottleneck identification，and used interval forms to describe these uncertain attributes of machines.Sengupta et al.［9］analyzed inter－departure time from differentmachines to identify and rank the bottlenecks.For bottleneck shifting problems，Roser et al.［10］believed that the bottleneck shifting was caused by the cumulative fluctuations of manufacturing systems.Thus，they proposed a detection method based on the equipment activity duration.Wang et al．［11］introduced an economic factor for the bottleneck identification，and proposed a new bottleneck indicator.Liu et al.［12］aimed to solve the prediction problem of production logistics bottleneck shifting. But some critical factors such as machine failures and setups，weren't considered in above studies.

Only for serial production lines，Betterton and Silver［13］suggested a practical approach to detect bottlenecks in open，asynchronous serial production lineswith finite buffers.

Although lots of studies have reported about the bottleneck identification and prediction methods，as detailed in the above literature，few researches have been done on evaluating the influence of performance on the perspective of the shifting bottleneck considering machine failures and setups in wafer fabrication systems.Therefore，this paper aimed to address the bottleneck shifting in transient states of the systems caused by the system random variability，and presented a bottleneck predictionmodel.

1 Problem Descriptions

This paper describes a method to detect and monitor the bottleneck machine in non－steady－state production systems subject to random variations.To describe the problem domain effectively，the follow ing assumptions are made.(1)The release time of wafers follows the Poisson distribution.The processing time is normally distributed，and wafers are dispatched according to the first in first out(FIFO)rule initially.(2)Machines must be repaired immediately when they break down，and the failure rates of machines follow the Weibull distribution.(3)The setup times follow the negative exponential distribution.(4)Each machine can process only one wafer at a time，and each wafer can be processed by only onemachine at a time.

The notationswere defined as follows.

i is the index of wafer type，i=1，2，…，Nr．

j is the operation index，j=1，2，…，pi．

mkis the index of amachine，k=1，2，…，M．is the j th operation of wafer i conducted on machineis the former operation ofof wafer i．is the latter operation ofof wafer i．is the latter operation ofon machine mk．is the former operation ofonmachine mk．

Niis the total number of wafer i．

Nris the number of differentwafer types.

M is the total number ofmachines.

IBN(mk)is the bottleneck index of mk．

IBN1(mk)is the workload index of mk．

IBN2(mk)is the system variability index of mk．

γhis theweight index of each factor in transientstates，h= 1，2，…，u．

λi(t)is the failure rate function of the i th periodΔt．

ωl(t)is the failure rate function after the l thmaintenance.

a is the improvement factor in age of the system reduced by maintenance(0＜a＜1)．

b is the improvement factor in failure rate of the system increased by maintenance(b＞1)．

Tl－1is the average time between two breakdowns when λi(t)=ωl－1(t)．

fi+1(t)is the failure probability density function of the (i+1)thΔt．is the averagemaintenance time of mk．

2 Bottleneck Prediction Mode l

The bottleneck index is defined based on theworkloads and the variability in a wafer fabrication system.The larger the bottleneck index is，the higher the probability of being bottleneck becomes.γ1andγ2are the weight coefficients of the two parts in transient states.

The throughput of a bottleneck machine determ ines the total throughput of a system.A machine that has significant effect on throughput should be regarded as a bottleneck machine.Thus，workloads play a significant role in the bottleneck identification.The expression of workloads of machine mkis depicted as:

To identify a bottleneck machine on the perspective of the variability，the bottleneck machine is caused by the follow ing factors，which occur unexpectedly during processing:machine breakdowns，setups，new arrivals，and other external influences.These factors serve to inflate both themean and the standard deviation of effective processing times，and the effect of the variability can be quantified as a random variable via them.

A variability index in awafer fabrication system in transient states can be depicted as:

where these external factors are related with time t and u is an integer value.Here，the twomain factors，machine failures and setups，are considered.The detailswillbe presented as follows.

2．1 Variability from machine failures

Machine failures will disrupt steady states of a wafer fabrication system and result in transient operations［14］. Unscheduled downtimes can greatly inflate both the mean and the standard deviation of effective processing times.These are themain causes of the variability.

Supposemachine mkmay have three states at time point T in a wafer fabrication system.(1)The machine mkis functioning well at time point T.(2)The machine mkbreaks down，and the repair is undergoing at thatmoment.(3)The machine mkis broken，and has just finished the repair or maintenance activity.

Supposing that the machine mkhas made maintenance activities l times at the time periodΔt(which starts from time point T)，the failure rate functionλi(t)is:

where Tl－1is the average failure timewhen themachine's failure probability isωl－1(t)．

State 1When the machine mkis functioning well，its failure rate will not change.Thus，the failure rateλi+1(t)of the machine mkin the next period can be depicted as:

Convertingλi+1(t)to the failure probability density function fi+1(t)of themachine mkin that period，the follow ing equation can be given:

where e is a constant.

The average time between two failures of themachine mkin that period can be described as follows:

where tlis the time point after the machine mkfinishes the l th maintenance activity.

The availability of themachine mkin the(i+1)th period is:

State 2Themachine mkbreaks down，and is just right during repairs，then the failure rateλi+1(t)of the machine mkin the next period can be depicted as:

Similarly，converteλi+1(t)to the failure probability density function fi+1(t)of themachine mkin that period:where e is a constant.

The average time between two failures of themachine mkin that period is:

State 3The machine mkbreaks down，and just finishes the repair ormaintenance activity，that is:tm，i+1=0，

The availability of themachine mkin the(i+1)th period is:

According to the variability theory［15］，the mean of effective processing time of themachine mkconsidering machine failures in the(i+1)th period is:

The variance of effective processing time of themachine mkin that period is:

The coefficient of the variation of effective processing time of themachine mkis:

2．2 Variability from setups

Setups can be regarded as non－preemptive outages when they occur due to the changes in the production process.On the basis of machine failures，this section analyzes the variability from setups.The FIFO rule was introduced as the initial scheduling policy.The earliest start processing time and the latest finish processing time can be calculated as follows:

It is known that in the time interval［T，T+Δt］，the first waferand the latestwaferon themachine m should satisfykthe follow ing constraints:

The total setup time of themachine mkin that period is:

where

The total number of setups of themachine mkis:

Themean of the effective processing time with setups is:

The variance of the effective processing time of the machine mkwith setups is:

The variability index of the machine mkin transient states is:

IBN2(mk)where tm，i+1is the repair time of machine mkin the(i+1)th period.

The availability of themachine mkin the(i+1)th period can be described as follows:

According to Eq.(1)，there is:

A bottleneckmachine can be identified as themachinewith the largest bottleneck index Skin that period:

The core of the bottleneck identification and prediction method can be summarized as follows.Firstly，the workload index IBN1(mk)of eachmachine can be calculated based on Eq. (2).Secondly，according to the machine failure conditions，there are three kinds of circumstances.The availability of the machine mkin the(i+1)th period can be calculated according to the current state of the system.What'smore，the variance and coefficient of the standard deviation of an effective processing time for themachine mkcould be obtained.Thirdly，calculate the total time and the number of setups of themachine mkin that period，and then compute the availability index IBN2(mk)considering themachine failures and setups.Fourthly，calculate the total bottleneck index on the basis of Eq.(1).The machine mkwhich has the largest bottleneck index Skis the bottleneck machine.Finally，take the current period as the initial state，and continue to identify the bottleneck machine in the next period.

3 Experim ents and Analysis

3．1 Bottleneck shifting analysis

The simulationmodel adopted in this paper is based on the model HP－24 of Wein［16］.The experiments are designed as follows.(1)There are 24 differentmachines and three types of wafers Pa，Pb，and Pcin the wafer fabrication system.Wafers are processed according to their processing routes，and the processing time follows normal distribution.(2)The failure function of the machine mkfollows the Weibull distribution U(2，4)and the scale parameterθiis 10 times of the value of uniform distribution U(60，130).(3)The improvement factors，a is1.05 and b is 0.02.(4)The setup times follow negative exponential distribution，whereλ=7．It is assumed thatγ1∶γ2=1∶1．

The experiments have been done on the HP personal computer with 160 GB hard disk，2GB DDR3 memory，and 1.58 GHz Intel CoreTM 2 processor.The simulation lasts for eight consecutive time intervalΔt，Δt=10 h．The results showed in Table 1 are the bottleneck indexes of 24 machines at eight different time periods.All the data are the average of 10 experiments.The machine which has the largest bottleneck index is deemed to be the bottleneck machine.

Table 1 Bottleneck index changes over time

It is showed in the above table that the bottleneck of the system has shifted over time.In the first four periods T1－T4，the largest bottleneck index is the bottleneck index of machine m6．So the bottleneck ismachine m6in this period.But in period T5，the bottleneck shifts to machine m12，and the machine m6becomes the sub－bottleneck.In the last period，the bottleneck shifts back to machine m6．This is caused by the random variation generated in the system，and it inflates themean of the effective processing time，resulting in the bottleneck shifting.

3．2 Sensitivity analysis

Themain causes of variability in wafer fabrication system aremachine failures and setups.This experiment analyzes their influence on the bottleneck shifting respectively.

Change the parameters of the machine failures while the setup times are constant，and observe the system bottleneck shifting situation.The results are presented in Table 2.The machine failures are 0.01λ，0.1λ，1λ，10λ，and 100λ．

Table 2 Bottleneck shifting with the failure rate

As shown in Table2，when themachine failure rate is0.01 λ，the system bottleneck is always m6．When the failure rate ofmachine increases to 0.1λ，the bottleneck changes during T5and T6periods，and machine m12becomes the bottleneck of the system.Butwhen themachine failure rate continues increasing，the bottleneck shiftsmore frequently:the bottleneck shifts from machine m6tomachine m18and then tomachine m12，and at last it shifts back to m6again.Hence，the failures ofmachines have a great influence on the system bottleneck，and the fluctuation accumulated by machine failures is the reason leading to bottleneck shifts，too.With the increase of machine failure rate，the frequency of machine breakdown increases，so the stability of the system reduces gradually.

The follow ing experimentexplains the relationship between the target value and differentsetup timeswhen the failure ratesλ is a constant.As shown in Table 3，similar to the influence of themachine failure rates，with the increase of setup times，the bottleneck shifting trend is remarkable.In addition，when the machine failure rates and the setup times are both very small，the system bottleneck remains the same.Itmeans thatwhen the system is in a stable state，the utilization of machine m6is the largest.

Table 3 Bottleneck shifting with the setup times

3．3 Comparison with the benchmark

To further analyze the effectiveness of the model，a simulation experiment is conducted as a benchmark on the Arena simulation platform.Take the utilization of machines as the index of performance to identify the system bottlenecks.The simulation time is one year，and the pre－simulation time is 3 months.The results of the prediction model are compared with the simulation results through the experiment to test the accuracy of the proposed algorithm.Take the rootmean square difference (RMSE)R to calculate the deviations of the prediction and simulation results.

where diis the deviation of the prediction and the simulation results of each time and n is the time epochs of the experiments.

The statistical results are shown in Fig.1.The error curve is falling with the increase of experiment times.Itmeans that the accuracy of the proposed prediction method increases gradually，and keeps balance at last.The error is within the acceptable range.

Fig.1 Error curve of themodel

Analyze the target valuewhen the failure rates ofmachines become 0.2λ，0.4λ，0.6λ，0.8λ，λ，1.2λ，1.4λ，1.6λ，1.8λ，2λ，and the setup times are 0.1 t，1 t，10 t．Here，t is the basic setup time for processing a wafer.The number of experimental groups is 50.

Fig.2 Prediction error curve with coefficient

As Fig.2 shows，the greater the probability of the system failure is，the longer the setup times are，and themore unstable the system is，the better the proposed predictionmethod is.

4 Conclusions

(1)In this paper，a bottleneck shifting predictionmodel is established considering the system variability in transient states caused by machine breakdowns and setups.

(2)The proposed bottleneck prediction model is verified by simulation experiments and the comparison analysis.Results indicate that the proposed bottleneck prediction method is feasible and effective.

(3)The sensitivity analysis shows that the frequency of bottleneckmachine shifting increaseswith the failure rates or the setup times.

References

［1］Zhang R，Wu C.Bottleneck Machine Identification Method Based on Constraint Transformation for Job Shop Scheduling with Genetic Algorithm［J］.Information Sciences，2012，188:236－252.

［2］Hu H T，Zhen L，Sun Z，et al.A Multi－stage Fluctuation Smoothing Method for Multiple Bottlenecks in Wafer Fabrication［J］.International Journal of Advanced Manufacturing Technology，2013，67(1/2/3/4):111－120.

［3］Law A M，Kelton W D.Simulation Modeling and Analysis［M］.3rd ed.New York:M cGraw－Hill，1991:35－49

［4］Zhou Z G，Rose O.A Bottleneck Detection and Dynam ic Dispatching Strategy for Sem iconductor Wafer Fabrication Facilities［C］.Proceedings of the 2009 W inter Simulation Conference，Austin，TX，USA，2009:1646－1656.

［5］Chang Q，Ni J，Bandyopadhyay P，et al.Supervisory Factory Control Based on Real－Time Production Feedback［J］.Journal of Manufacturing Science and Engineering，2007，129(3):653－660.

［6］Chiang SY，Kuo C T，Meerkov SM.DT－Bottlenecks in Serial Production Lines:Theory and Application［J］.Robotics and Automation，2000，16(5):567－580.

［7］Li L，Chang Q，Ni J.Data Driven Bottleneck Detection of Manufacturing Systems［J］.International Journal of Production Research，2009，47(18):5019－5036.

［8］Wang J Q，Chen J，Wang S，et al.Interval Multi－attribute Bottleneck Identification in Job Shop［J］.Computer Integrated Manufacturing System，2013，19(2):429－437.(in Chinese)

［9］Sengupta S，White T，Das K，et al.Analysis of a New Signal for Bottleneck Identification and Loss Allocation to Individual Machines［J］.International Journal of Industrial and Systems Engineering，2013，13(2):175－196.

( Continued on page 558)

Outlier Detection Algorithm Based on Iterative Clustering

GU Ping(古平)*，LUO Xin(羅辛)，YANG Rui－long(楊瑞龍)，ZHANG Cheng(張程)

Department of Computer Science and Technology，Chongqing University，Chongqing 400044，China

Abstract:A novel approach for outlier detection with iterative clustering(ICOD)in diverse subspaces is proposed．The proposed methodology comprises two phases，iterative clustering and outlier factor computation．During the clustering phase，multip le clusterings are detected alternatively based on an optim ization procedure that incorporates term s for cluster quality and novelty relative to existing solution．Once new clusters are detected，outlier factors can be estimated from a new definition for outliers(cluster based outlier)，which provides importance to the local data behavior．Experiment shows that the proposed algorithm can detect outliers which exist in different clusterings effectively even in high dimensional data sets．

Keywords:clustering;outlier detection;dimensional reduction

CLCnumber:TP301.6Documentcode:AArticleID:1672－5220(2015)04－0554－05

Introduction

Outlier detection is a fundamental task that is useful in applications of scientific，engineering or business databases.It deals with the problem of identifying rare points that widely diverge from the general behavior ormodel of the data.The task of outlier detection has been tackled by researchers in a variety of ways［1－4］，and these can be broadly classified into four approaches:distribution，distance，density，and clustering based.

Distribution based method［5］is mostly used in early studies，which assumes a distribution or constructs an explicit model of the data.Significant deviations from the assumed distribution or the constructed model are identified as outliers. The notion of distance－based outliers is proposed in Ref.［6］. An outlier is defined as“object O in dataset T is an outlier if at least a fraction q of the other objects in dataset T lies greater than distance D from O”.This definition can identify global outliers effectively，but cannot detect local outliers if the data set consists of clusters of diverse density.To overcome this problem，Breunig et al.［7］proposed a local outlier factor (LOF)for each object in the data set and declared the top n points in this ranking to be outliers.Papadimitriou et al.［8］proposed a variant of the LOF called multi－granularity deviation factor(MDEF).For a given record，its MDEF is calculated as the standard deviation among its local density and the local densities of its k nearest neighbors.In the same lines，Jin et al.［9］proposed another variant of LOF that improved efficiency by avoiding unnecessary calculations.They achieve this by calculating upper and lower bounds among the micro clusters detected.

In the last decades，several clustering algorithms［10－12］have been extended to identify objects significantly deviating from the others as outliers.In Ref.［10］，He et al．defined cluster－based local outlier factor(CBLOF)，a measure for identifying the significance of an outlier and an algorithm for discovering outliers with CBLOF.After that Duan et al.［11］proposed a cluster－based outlier detection algorithm which could detect both single point outliers and cluster－based outliers，and assign each outlier a degree of being an outlier.Shi and Zhang［12］proposed a cluster－outlier iterative detection(COID) algorithm.In thisalgorithm，clusterswere detected and adjusted according to the intra－relationship within clusters and the interrelationship between clusters and outliers iteratively.

However，in high dimensional spaces，traditional clustering methods suffer from the curse of dimensionality，and sim ilar and dissimilar objects cannot be discriminated well over the full attribute space.Though subspace clustering［13－15］can alleviate this problem，most algorithms work only with one clustering solution，while in many applications，complex data can be grouped and interpreted in many different ways，i．e．，the individual objects are described from several perspectives or several views，each of which highlights different characteristics of the outliers.

In this paper，we propose a novel approach for outlier detection with iterative clustering(ICOD).For each iteration，we try to find one clustering based on an optimization procedure that incorporates cluster quality and novelty into consideration，and then outliers are identified by exploratory outlier analysis. With multiple projected subspaces，we can discover outliers that are not detectable in previous clustering solutions.

1 Prelim inary

Inmany situations，the outlierness of an object can always be observed from different views.For example，consider a database of personswhich shares three combinationsof attributes X1，X2，and X3.In subspace corresponding to attribute X1and X2(Fig.1(a))，there are two clusters，and three persons are identified as an outliers(represented by“+”symbol) correctly，while from Fig.1(b)，it can be observed that addition of attribute X3has resulted in the concealment of the real clusters and real outliers and such a situation leads to incorrect outlier detection.This result suggests that different clustering from different subspaces can always lead to different outliers.

(a)Plot of the subspace corresponding to X1and X2

Received date:2014－03－28

Foundation items:Natural Science Foundation Project of CQ CSTC(Nos.cstc2012jjA40002，cstc2012jjA40016);Fundamental Research Funds for the Central Universities，China(No.0216005207016)

*Correspondence should be addressed to GU Ping，E－mail:guping2k@cqu.edu.cn

Fig.1 Plots of relationship between clustering and outliers in different subspaces

To adapt to this property，we improve traditional outlier detection methods by taking advantage of multiple clustering views.Consider a group of ndata points{ x1，x2，…，xn}，where xiis a data pointand d is the dimensionality of the feature space.Let current clustering bewhere c is the number of clusters，and cluster labelingmatrix iswith element lij=1 if xibelongs to cluster j in P1， otherw ise it is0.Our goal is to discover potential outliers from the new clustering solution U given previous labelingmatrix L=［L1，L2，…，Lt－1］in each iteration t.

2 The ICOD Algorithm

To detect outliers which are concealed in the multiple clustering solutions，the ICOD algorithm iteratively performs the follow ing two steps:(1)discover a novel clustering solution U in subspace W with spectral clustering［16］;(2)compute and rank objects based on the outlier score，and update outlier set incrementally.

2．1 Discovery of iterative clustering

During searching for the quality alternative clustering U，we combine the cluster quality and diversity criteria into a single optimization function，and compute the similarity between samples in a reduced dimensional subspace W．

2．1．1 Object function with cluster quality

To work on arbitrarily shaped clusters，we build our work on spectral clustering.In Ref.［17］，spectral clustering is considered as the k－way normalized cut problem.Given a graph G=(V，E，K)，where V is the set of vertices，E is the setof edges connecting vertices，and K is an edge affinity matrix. Suppose P，Q?V，and the normalized link ratio of P，Q is linkratio(P，Q)=link(P，Q)/link(P，V)，where link(P，

The k－way normalized cut problem is tominimize the links that escape a cluster relative to the totalweight of the cluster:

Let D be the diagonalmatrix whose(i，i)entry is the sum of the entries of row i inmatrix K．Relaxing the indicatormatrix to allow its entries to take on any real value，then the normalized cut criterion is equivalent to the follow ing trace maximization problem:

A well－known solution to this problem is obtained by setting thematrix U to be the top k eigenvectors of the matrix D－1/2KD－1/2．These eigenvectors are then used to compute a discrete partitioning of the points by applying k－means algorithm.

2．1．2 Object function with cluster novelty

This object is to ensure that the clustering found is novel compared with existing solutions，and we measure the dependence between them in terms of a kernel dependence measure，the Hilbert－Schm idt Independence Criterion (HSIC)［18］.Given two random variables x and y，let us define amappingφ(x)from x∈X to kernel space F，such that the inner product between vectors in that space is given by a kernel function k1(x，x')=〈φ(x)，φ(x')〉．Let G be a second kernel space on Y with kernel function k2(·，·)．With n independent observations，we can approximate the HSIC dependence by:

where K1and K2are the Gram matrices and Hij=δij－n－1．2．1．3 Object definition and optim ization

Given previous labeling matrix L=［L1，L2…，Lt－1］，in iteration t，we try to find an alternative clustering solution U in subspace W by optim izing the follow ing object function: where I is the identity matrix.Maximizing function(4)is achieved when the first term is large(corresponding to Eq.(2) which implies a good clustering)and the second term is small (small dependence between XW and L implies a diverse solution with L)．

The optimization of function(4)can be split into two parts:determ ining U under current subspace W and optim ization of W given solution U．Assume the feature projection W is fixed，and we can calculate the similarity and degreematrices K and D．Similar to spectral clustering，herewe relax the indicatormatrix U to take on real values.The solution for U is equal to the largest c eigenvectors of matrix D－1/2K D－1/2，where c is the number of clusters.

When U is fixed，optim izing the function(4)with respect to W is a nonlinear optim ization problem with orthonormality constraints.We utilize a dimension grow th algorithm to optim ize W．Initially one－dimensional vector w1is defined by random projection，and then we increase the dimensionality by one and get w2，which is initialized by random projection and then projected to the space orthogonal to w1．The process is continued for w1，w2，…，wpuntil the desired number of dimensions p is reached.

2．2 Cluster based outlier score

Based on the partition Ptfound in each clustering，we tried to assign to each objectan outlier score to demonstrate its degree of outlierness.

Definition 1(Large and small cluster)Let the current clustering isdenotes the number of objects inWe define b as the boundary of large and small clu ster if formu lasnβanddefine the setof large clusters as:and the set of small clusters as:

Ifβis set to 95%，we intend to regard clusters which contain 95%of data points as normal.However，this guideline is not always appropriate，in some cases abnormal cluster deviated from a large cluster might contain more points than a certain small normal cluster.We address this problem by

Definition 3(Cluster based outlier score)For each object q∈∈SC，let its nearest large cluster be CLN，d(q，CLN)be the distance between q and cluster CLN，and the cluster based outlier score can be defined as:

It is easy to see that the smaller of its cluster size，the farther the distance between the object and its nearest large cluster，the more dense its closest cluster is，the higher the Od(q)，and object q ismore likely to be an outlier.

2．3 Outline of algorithm ICOD

Input:existing cluster labeling L=［L1，L2…，Lt－1］，cluster number c，and the iteration number N'.

Output:top－m objects in outlier set O;

3 Com plexity of ICOD

Let the size of the data setbe n，and for each iteration，the complexity of clustering solution discovering is mainly determined by calculating of the kernel similarity matrix K and the eigen－decomposition of K．Since the eigenvalues of K drops fast，we can find low rank approximations of K with rank s (s＜＜n)，so the complexity is O(ns2)．The computation of subspace W has complexity of O(nsrp)，where p is the number of dimensions，and r is the steps in the dimension grow th procedure.In the second step of ICOD，the complexity of kmeans clustering is O(npc)，and the computation of outlier score needs O(n2)(with k－d tree，it can be reduced to O(n log n))．Hence，for t iterations，the total complexity of algorithm ICOD is O(ns2+nsrd+n log n+npc)t．

4 Experim ents

In this section，we conduct experiments to compare our ICOD algorithm with CBLOF［10］and COID［12］on three data sets(WBCDD，Iris，and CD27+).The WBCDD data set contains 569 records with 30 attributes in two classes(benign，malignant).We use all the 357 instances in benign class as normal objects，and random ly choose 20 instances from malignant class as outliers.The Iris data set contains 150 recordswith 4 attributes in 3 classes.We choose two classes (Setosa and Versicolour)as normal class，and 20 instances from class Virginica as outliers.CD27+is a gene expression data set，which contains9 209 recordswith 25 attributes(among them only 4 000 normal instances and 20 outliers are random ly selected).For all experiments，we apply the Gaussian kernel and set c four times the number of class，β=0．8．The evaluation ismainly based on twometrics:detection rate［20］and receiver operating characteristic(ROC)curve.

In the first experiment，the ICOD algorithm is evaluated with increasing number of clustering iterations for the top 32 objects.We can see from Fig.2，that detection rate of ICOD continuously increases along with the number of clustering iterations，and this confirms that iterative clustering can indeed find some new outliers which are ignored in the single clustering.For example，in data setof CD27+，we observe up to 15%improvementof detection rate in iteration 4.The reason behind is thatwe notonly find outliers in dimensions1 and 3 as obtained in iteration 3，but also obtain new outliers in dimensions 11，17，and 22.

Fig.2 Detection rate of ICOD as the clustering iterations is increased

As detection rate is highly dependenton the choice of m for top－m outliers，we conducted an experiment to compare the performance of all algorithms under different values of m.The result in Table 1 shows that，inmost cases，COID is superior to CBLOF.A closer look at the result shows that adjustment of clusters and outliers iteratively in COID can really contribute to the improved detection results，especially when outliers are close to the boundary of clusters.However，we cannot get better result through setting the value m higher，for example，even when m is set to 48 in WBCDD，there is still two outliers undiscovered.On the contrary，when m just reaches 40，algorithm ICOD finds all20 outliers(including two outliers exist in subspace spanned by dimensions6，14，20，and 22)，which reveals the importance of subspace analysis in uncovering complex deviations hidden in different subspace projections.On the other data set CD27+，similar results can be obtained: algorithm ICOD outperforms COID.Meanwhile，COID is superior to CBLOF inmost cases.

Figure 3 is the ROC curve of all outlier detection algorithms on data set WBCDD.Result shows that in most case，algorithm ICOD reaches a higher true positive rate earlier than other competitors，which proves that algorithm ICOD can obtain a better tradeoff between the coverage and detection rate.

Table 1 Outliers detection rate in data setsWBCDD and CD27+

Fig.3 ROC curves of ICOD，CBLOF，and COID on WBCDD

Table 2 shows the running time of algorithms on three datasets with varying sizes.As can be seen from the table，among all three algorithms，CBLOF performsmuch faster than the two others and ICOD achieves similar efficiency to COID. The result ismore obvious in particular in CD27+.The main difference lies in that CBLOF requires only single pass clustering in outlier detection，while COID and ICOD need to refine the set of clusters iteratively.Though inferior to CBLOF in efficiency，algorithm ICOD scales well with the increasing number of data，and this demonstrates the effectiveness of subspace analysis in iterative clustering even in faced with large dataset.

Table 2 Running time vs data set size

5 Conclusions

Clustering is an attractive approach for outlier detection. However，previous researches fail to utilize multiple clustering solutions to find complex outliers that are hidden in multiple subspaces.In this paper，a new outlier detection algorithm which relies onmultiple clustering results in different subspaces is proposed.For each iteration，we try to find one clustering based on an optim ization procedure that incorporates cluster quality and novelty into consideration，and then outliers are identified by exploratory outlier analysis.With these multiple projected subspaces，we can discover outliers that are not detectable in previous clustering solutions.Experimental results on several data sets also demonstrate the effectiveness of our method.For future work，we will integrate subspace clustering algorithm more tightly with outlier analysis to make the detecting processmore efficient.

References

［1］Hido S，Tsuboi Y，Kashima H，et al.Statistical Outlier Detection Using Direct Density Ratio Estimation［J］.Know ledgeInformation System，2011，26(2):309－336.

［2］Gao J，Hu W，Zhang Z F，et al.RKOF:Robust Kernel－Based Local Outlier Detection［J］.Advances in Know ledge Discovery and Data Mining，Lecture Notes in Computer Science，2011，6635:270－283.

［3］Wu M，Jermaine C.Outlier Detection by Sampling with Accuracy Guarantees［C］.Proceedings of the 12th ACM SIGKDD International Conference on Know ledge Discovery and Data M ining，Philadelphia，2006:767－772.

［4］Yang J，Zhong N，Yao Y Y，et al.Local Peculiarity Factor and Its Application in Outlier Detection［C］.Proceedings of the 14th ACM SIGKDD International Conference on Know ledge Discovery and Data M ining，Las Vegas，USA，2008:776－784.

［5］Hido S，Tsuboi Y，Kashima H，et al.Statistical Outlier Detection Using Direct Density Ratio Estimation［J］.Know ledge and Information Systems，2011，26(2):309－336.

［6］Sadik M S，Gruenwald L.DBOD－DS:Distance Based Outlier Detection for Data Streams［J］.Computer Technologies and Information Sciences，2011，6261:122－136.

［7］Breunig M M，Kriegel H P，Ng R T，et al.LOF:Identifying Density－Based Local Outliers［J］.Sigmod Record，2000，29 (2):93－104.

［8］Papadim itriou S，Kitagawa H，Gibbons P，et al.Faloutsos，Loci:FastOutlier Detection Using the Local Correlation Integral，Technical Report IRP－TR－02－09［R］.Pittsburgh，USA:Intel Research Laboratory，2002.

［9］Jin W，Tung A K H，Han J.M ining Top－n Local Outliers in Large Databases［C］.Proceedings of the 7th ACM SIGKDD International Conference on Know ledge Discovery and Data M ining，San Francisco，CA，USA，2001:293－298.

［10］He Z Y，Xu X F，Deng S C.Discovering Cluster－Based Local Outliers［J］.Pattern Recognition Letters，2003，24(9/10): 1641－1650.

［11］Duan L，Xu L D，Liu Y，et al.Cluster－Based Outlier Detection［J］.Annals of Operations Research，2009，168(1):151－168.

［12］Shi Y，Zhang L.COID:a Cluster－Outlier Iterative Detection Approach to Multi－dimensional Data Analysis［J］.Know ledge and Information Systems，2011，28(3):709－733.

［13］Muller E，Schiffer M，Seidl T.Statistical Selection of Relevant Subspace Projections for Outlier Ranking［C］.International Council for Open and Distance Education，Hannover，Germany，2011:434－445.

［14］Xue Z X，Shang Y L，F(xiàn)eng A F.Sem i－supervised Outlier Detection Based on Fuzzy Rough c－Means Clustering［J］. Mathematics and Computers in Simulation，2010，80(9):1911－1921.

［15］Seidl T，Muller E，Assent I，et al.Outlier Detection and Ranking Based on Subspace Clustering［M］.Dagstuhl，Germany:Uncertainty Management in Information Systems，2009.

［16］Niu D L，Dy JG，Jordan M I.Iterative Discovery of Multiple Alternative Clustering Views［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2014，36(7):1340－1353.

［17］Dhillon I S，Guan Y Q，Kulis B.Kernel k－Means:Spectral Clustering and Normalized Cuts［C］.Proceedings of the 10th ACM SIGKDD on Know ledge Discovery and Data M ining，New York，USA，2004:551－556.

［18］Gretton A，Bousquet O，Smola A，et al.Measuring statistical dependence with Hilbert－Schm idt Norms［C］.International Conference of Algorithm ic Learning Theory，Singapore，2005: 63－77.

［19］Duan L，Xu L，Guo F，et al.A Local－Density Based Spatial Clustering Algorithm with Noise［J］.Information Systems，2007，32(7):978－986.

［20］Cerioli A，F(xiàn)arcomeni A.Error Rates for Multivariate Outlier Detection［J］.Computational Statistics＆Data Analysis，2011，55(1):544－553.

［10］Roser C，Nakano M，Tanaka M.Shifting Bottleneck Detection［C］.Proceedings of the 2002 W inter Simulation Conference，Aichi，Japan，2002:1079－1086.

［11］Wang J Q，Chen J，Wang S，et al.Shifting Econom ic Bottleneck Identification［C］.IEEE International Conference on Industrial Engineering and Engineering Management(IEEM)，Singapore，2011:1760－1764.

［12］Liu Z，Tang J，F(xiàn)ei Z M.Close－Loop Prediction Method of Logistics Bottleneck Based on Bottleneck Polymorphism［J］. Computer Integrated Manufacturing System，2012，18(11): 2554－2561.(in Chinese)

［13］Betterton C E，Silver S J.Detecting Bottlenecks in Serial Production Lines—a Focus on Inter－Departure Time Variance［J］.International Journal of Production Research，2012，50 (15):4158－4174.

［14］Park K，Morrison JR.Control of Wafer Release in Multi Cluster Tools［C］.8th IEEE International Conference on Control and Automation，Xiamen，China，2010:1481－1487.

［15］Hopp W J，Spearman M L.Factory Physics［M］.3rd ed.New York:M cGraw Hill Higher Education，2002:213－287.

［16］Wein L M.Scheduling Sem iconductor Wafer Fabrication［J］. IEEE Transactions on Semiconductor Manufacturing，1988，1 (3):115－130.

Integration of Berth Allocation and Quay Crane Assignment in Tidal Container Ports

LU Zhi－qiang(陸志強)*，WUWen(吳文)，HAN Xiao－le(韓笑樂)，NIRAVONG Juliane

School of Mechanical Engineering，Tongji University，Shanghai201804，China

Abstract:Tide is a significant factor which interferes with the berthing and departing operations of vessels in tidal ports．It is a preferable way to incorporate this factor into the simultaneous berth allocation and quay crane(QC)assignment problem(BACAP)in order to facilitate the realistic decision-making process at container term inal．For this purpose，an integrated optim ization model is built with tidal time w indows as forbidden intervals for berthing or departing．A hind-and-fore adjustment heuristic is proposed and applied under an iterative optim ization framework．Numerical experiment shows the satisfying performance of the proposed algorithm．

Keywords:berth allocation;quay crane(QC)assignment;integrated scheduling;tidal container ports

CLCnumber:TP29;U691Documentcode:AArticleID:1672－5220(2015)04－0559－06

Introduction

Since seaside operations have significant impact on the overall production efficiency of container terminal，adequate distribution of seaside resources to vessels will benefit both the terminal and clients a lot.Moreover，in consideration of the lim itof term inal size and the pricey investment，berths and quay cranes(QCs)are typically the scarcest resources at container terminal.Thus，an outstanding integrated decision of berth allocation problem(BAP)and quay crane assignment problem (QCAP)willmake considerable sense.

Meanwhile，tide is an important element for the daily operation in tidal ports.For vessels in harbor，the scour of tides has negative impacton the efficiency of seaside operationswhen tides rise and fall.Fortunately，this problem can be easily addressed by taking some fastening measures in most cases. When it comes to the berthing and departing operations，rigorous water condition is also a must in addition to high coordination between the main vessel and tow ing vessels. Inappropriate timing for berthing or departing will make big trouble for the operation.In severe cases，both the vessels and the quay will be destroyed.So，in order to ensure the operational security and efficiency of berthing and departing operations in tidal port，optimal choice of time w indow is of vital importance.

The“double w indow departure”，i．e．，vessels being able to leave the portatboth the times of high and ebb tides，is a big challenge for all ports in the world，and has been broken through in Yangshan Deepwater Port(China)in 2013.This experimentation is conducive to cutting down the in－berth waiting time of vessels by 3 to 6 h.However，rigorous berth conditions and essential technique for such departing operations will not be achieved without huge investment.Regardless of large ports like Yangshan Deepwater Port，in hundreds of smallandmedium－sized ports affected severely by tides，vessels have to wait for hours to meet expedient berthing and departing conditions in order to guarantee the safety of operations，no matter when they want to berth or depart.Taicang Port(China) happens to be one of them.Besides，available depth is not always adequate for the movement of vessels，so they have to wait for sufficient depth at high tide in some other tidal ports. All these above show that tide is a crucial factor that cannot be ignored in the daily operation of the port.Consequently，it is necessary to give consideration to tidal effect in scheduling to provide a more realistic decision support for the operational management in a tidal port.

Few literatures about tide can be detected in researches of scheduling problem about container terminal.Xu et al．［1］considered a BAP of both static and dynamic cases in container terminals，in which the berth allocation to vesselswas lim ited by water depth and tidal condition.Then the time horizon is divided into two periods，named low－water period，where the vessel assignment is more restrictive，and high－water period，respectively.Barros et al．［2］considered the BAP in tidal bulk ports with stock level conditions，in which a berth was collocated discretely，the planning horizon was divided into several tidal time w indows(TTWs)and the draft conditions was decided by stock level additionally.

Literaturesmentioned above only considered tide in BAP. Ithad the same imperfection as that researchers studied BAPand QCAP separately at the very beginning.However，BAP and QCAP are basically interrelated.Number of the rest QCs will determine whether a vessel with restriction of certain QCs number can berth or not by solving the QCAP，while the allocation of berth position in BAP will conversely decide the available time and number of QCs.The strong interrelationship between the utilization of the quay space resource and the QCs rightly motivates an integrated solution approach to the two problems，namely berth allocation and quay crane assignment problem(BACAP).

BACAPwas firstly considered by Park and Kim［3］.This study dealt with BACAP as two independent problems:BAP and QCAP.Imai et al．［4］introduced a formulation for the simultaneous berth and crane allocation problem and found an approximate solution for the problem with genetic algorithm. BACAP with both variable in－time QC－to－vessel assignments and crane productivity losses by interference were modeled by Meisel and Bierw irth［5］.Moreover，the crane productivity depended on the berthing positionsof vessels.The objectivewas them inimization of service quality costs incurred by speedups of vessels and by tardy departures plus operational cost of the utilized QC－h(huán)ours.Raa etal．［6］presented amixed integer linear programming model for the BACAP taking into account vessel priorities，preferred berthing locations and handling time.Yang et al．［7］solved the BACAP in amulti－user container terminal. Interactions between berth allocation and QC assignment are important consideration.In Ref.［8］，a hybrid model was proposed to solve the collaborative scheduling of BAP and QCAPwith the objective to m inimize the sum of extra cost of berthing at non－optimal location and penalty cost of time delay. Recently，an integrated heuristics－based solution methodology was proposed that tackled the BACAP in Ref.［9］and the experimental results showed that the proposed approach could yield high quality solutions in a reasonable computational time.

Received date:2014－05－06

Foundation items:National Natural Science Foundations of China(Nos.70771065，71171130，61473211，71502129)

*Correspondence should be addressed to LU Zhi－qiang，E－mail:zhiqianglu@#edu.cn

Furthermore，somemore practical issues，which had significant effects on quay side scheduling problems，were taken into account.The impact of number of allocated QCs on port operational cost and a vessel's fuel consumption and em issions was analyzed in Ref.［10］Besides，in order to solve the problem of allocation of berth and QC simultaneously，Chang et al．［11］developed a dynamic allocation model using objective programm ing based on rolling－h(huán)orizon approach and Gui et al．［12］developed a two－level iterative algorithm by separating decision objections.Nonetheless，all the studies don't involve tide，which is a crucial practical issue in tidal.

Since tide makes a big influence on the berthing and departing allowance of vessels，itwill further cause variation of berth position and berth time，and QC numbers should be adjusted to get a more reasonable scheduling solution. Therefore，integrated scheduling with tidal time willmake a big difference.

The remainder of this paper is organized as follows.In section 1，based on several related assumptions，the mathematicalmodel for the integrated scheduling problem with special constraint is built.Then the mentioned iterative optimization algorithm framework and the corresponding two heuristics are proposed respectively in section 2.In section 3，numerical experiment is exam ined and the analysis of results is carried out.The conclusions are summarized in section 4.

1 Problem Description

There aremany typesof tides，such as diurnal，sem idiurnal and mixed tide.Ports in different places have different types of tides.In fact，the length of tidal time and the frequency of tide all change over time in the same port.In order to facilitate research，this paper proposes a tidal time w indow，in which berthing and departing are forbidden.The length of TTW is adjustable but the frequency is a constant(once every 12 h). TTW starts at time E(a parameter，used to indicate the start point of TTW)and ends at the end of each minimum time period(12 h).

As the whole planning horizon is divided into many sections，the complexity of the scheduling problem rises precipitously.The integrated scheduling of BACAP in a tidal port can be described as an integrated scheduling with TTW. Tide becomes a restriction factor for the berth allocation and release in this problem.Subsequently，the QC assignment will also be affected and the whole scheduling policy must be adjusted aswell.

1．1 Assum ptions

Some important assumptions are given as follows for modeling.

(1)We consider continuous berth layout for a high level utilization of berth resources and dynamic arrival of vessels to im itate amore realistic scheduling situation［13］.

(2)The number of QCs serving a vessel simultaneously is often restricted by a minimum number and a maximum number［5］.

(3)Handling time of a vessel depends on the number of QCs assigned to it and the handling tasks must be finished without interruption after starting.

(4)The cranes are supposed to be lined up alongside the quay in the same rail.

(5)The number of QCs assigned to a vessel does not change during its stay at the berth［14］.

(6)A reassignment of QCs is allowed among vessels in berth without changing the QCs number of assigned vessels when a new vessel arrivals.The setup time of QC is ignored in this research.

(7)Tide appears once every 12 h and the fluctuation of tidal time duration is negligible.

(8)TTW is set to 3 h for research in this paper.No berthing or departing is allowed in tidal time.

1．2 Notations and scenario exam p le

In order to develop this analysis and research easily，this paper divides the quay in segments of 10 m and the time in periods of 1 h.Notations are defined as follows.

Input data:

i is the index of vessels;

t is the index of time;

V is the set of vessels，V={1，2，…，i，…，V};

Q is the number of available QCs;

L is the length of quay(number of berth segments);

H is the planning horizon;

T is the set of time period，T={1，2，…，t，…，H};

E is the start point of tidal time;

Liis the length of vessel i;

Aiis the arrival time of vessel i;

Miis the crane capacity demand of vessel i given as a number of QC－h(huán)ours;

M is a large positive number;

and

are nonnegative integers.

Decision variables:

biis an integer，the berthing position of vessel i;

siis an integer，the start time of vessel i;

eiis an integer，the departure time of vessel i(not always finishing time);

ritis binary，set to 1 if at least one QC is assigned to vessel i at time t，otherw ise 0;

yijis binary，set to 1 if vessel i is berthed below vessel j，i．e．，bi+Lj≤bj，otherw ise 0;

zijis binary，set to 1 if vessel i departs not later than the startof vessel j，otherwise 0;

qiis the number of QCs assigned to vessel i;

To analyze this problem easily，four kinds of vessels' arrival and departure scenarios in tidal condition are shown in Fig.1.In this diagram，the abscissa axis shows time and ordinate axis represents the quay.The verticalbar in greymeans the length of timew indow.Priorities(i，j，k，p)，arrival times (Ai，Aj，Ak，Ap)，and QC number(qi，qj，qk，qp)of the four vessels are all given here.The priority list in this paper is an operational sequence of the given number vessels.As can be seen from Fig.1，the dotted rectangles represent the initial states of vesselswithout tidal factor.For instance，the initial state for vessel k is k0.According to the definition of the coordinate axes，the length and the width of a rectangle represent the service time and the berth amount allocated to the vessel，respectively.Since vessels generally berth along the length direction and the w idths of vessels have no effect on the quay side scheduling problem，vessels'w idths will not be considered in this paper.Obviously，berths are not enough for four vessels to berth incipiently.Solid rectangles are the adjusted position after tide is taken into consideration.The length of the oblique line filled rectangles represents the additional turnaround time after adjustment.The way shown in this chart is a w idely used one to deal with tide in tidal port now，which is named rightshifting manipulation.This manipulation means to get a complete scheduling solution without taking tide intoconsideration(TTW=0)firstly and then make an adjustment of the periodic solution to get a final resultwith the tidal factor (TTW=3 h).Besides，parameterswith superscript“'”mean the corresponding values have changed due to the adjustment. The given QC number for vessels i，j，p，and k are respectively set to 2，4，1，and 3 here.

Fig.1 Scenarios after right－shifting adjustment

1．3 M athematicalmodel

We formulate the integrated scheduling problem of BACAP with TTW as follows:

According to the assumption(1)，the dynam ic arrival of vesselsmeans that the arrival time of vessels is fixed，so the objective function is revised to:

subject to

Constraint(3)ensures that the QC capacity demand of all vessels can be satisfied.Resource constraint of QC is shown by constraint(4).Constraints(5)－(7)set the start times and the finish times for serving vessels without preemption.Constraints (8)－(9)determine that vessels can berth and depart safely with tide.Constraints(10)－(12)set variables yijand zijto avoid overlapping the handling of vessels in the space time diagram.The relationship among xit，rit，and qiis expressed by constraints(13)and(14).As shown in constraint(15)，the planning horizon H defines a lim it on the start time and the departure time of vessels.Constraint(16)ensures that each vessel is positioned within the berth boundaries.Further constraints(17)－(20)define the domains of the remaining variables and parameters.

2 The Proposed Solution Technique

2．1 Framework of the iterative algorithm

In the iterative optimization algorithm，the priority list and QC number(determ ined or suggested)assigned to each vessel are generated by a modified genetic algorithm.The follow ing procedure(in Fig.2)is the outline of the iterative optim ization framework.Parameters in Fig.2 are as follows:λis the number of iterations，λmaxis the stop criterion，φis the number of unimproved iterations，φmaxis the threshold of the regeneration， xbis the current global best solution，is the best solution of theλth iteration，C(xb)is the objective function value of xb，is the objective function value ofis thevalue obtained by hind adjustment heuristic(HAH) algorithm，andis thevalue obtained by fore adjustment heuristic(FAH)algorithm.

Fig.2 Procedure of the iterative framework

Details of HAH and FAH will be introduced in the section below.

The partial solution of this problem is represented as a chromosome，which is a string consisting of 2V integer items.The firstV items in this chromosome represent thepriority list and the lastV items are corresponding QC number assigned to each vessel(shown in Fig.3).

Fig.3 Example of a chromosome

Partial matched crossover and inversion mutation are adapted.The probability of crossover and mutation is calculated using the adaptive algorithm in Gui etal．［12］The conceptof the regeneration mechanism is also learned from that paper to avoid premature.For initialization，three additional chromosomes are included in the initial population.Their priority list is generated by the following three commonly used priority rules:FCFS，SPT，and LPT.Chromosomes are random ly selected and copied to the next generation using the Roulette Wheel selection procedure.The probability of a chromosome being copied is proportional to the fitness value related to the objective function value.The fitness value of individuals is calculated using Eq. (21)［15］，where F(ch)is the fitness value of chromosome ch andαis a problem scale related parameter.

2．2 HAH

In the proposed HAH algorithm，particular adjustment is made for some vessels after completing all the loading and unloading operations in order to meet the requirement of tidal constraints.Obviously，the additional turnaround time of vessels is cut down due to this adjustment in Fig.4.

Fig.4 Scenarios after HAH adjustment

Corresponding framework is presented in Fig.5 and details are introduced in the follow ing steps.

Fig.5 Flow chart of HAH

Step 1 Select the first unscheduled vessel i in the priority list produced by genetic algorithm in the upper layer.

Step 2 Mark the arrival time Aiof vessel i and set the start time siequal to Ai．

Step 3 If vessel i is not allowed to berth(siis in TTW)，postpone the start time of vessel i to the end of the current TTW and go to step 4;otherw ise，go to next step directly.

Step 4 Assign qiQCs to vessel i and calculate the departure time eiof it.

Step 5 Set the berthing position biof vessel i to zero.

Step 6 If overlapping exists，skip to step 8;otherw ise，go to next step.

Step 7 If QCs are sufficient，the operation planning of vessel i is finished，skip to step 12;otherw ise，go to the next step.

Step 8 Move the berthing place one segment every time along the quay.

Step 9 If the summation of the berthing place biand the vessel length Liis larger than the length of the quay L，go to next step;otherw ise skip to step 6.

Step 10 Postpone the start time sito the next1 h.

Step 11 If the start time siis larger than the planning horizon，fail to get the feasible solution;otherw ise go back to step 5.

Step 12 If the departure time eiis not in the TTW，then it is safe to depart from the berth;otherw ise，go to the next step.

Step 13 Wait until the end of this TTW to leave the port safely.

2．3 FAH

Since the operating time(total time of loading and unloading)of a vessel is contingent on the number of QCs assigned to it，the departure time will also change accordinglyunder the ratified start time.Thismeans that the departure time of vessels is possible in the TTW.In other words，the inalterable QC assignment decided in the upper layer will probably lead to additionalwaiting in the berth and the excessive consumption of QC resource.So，it is a good way to determine the number of the immediate available QCs timely asmany as possible and guarantee that the vessel can depart immediately after the accomplishment of all operations.Then several modifications have been carried out through trial and error. Finally，F(xiàn)AH has been proven effective enough in getting a better solution when combined with HAH.The result of the abovementioned scenarios is shown in Fig.6.The suggested QC numbers for vessel i，j，p，and k have been respectively changed to 3，4，2，and 3.The length of the vertical line filled rectangles represent the reduced turnaround time after FAH is implemented.The algorithm shows obvious superiority in this example.

The process of FAH is depicted in the flow chart shown in Fig.7.Compared with HAH，detailed steps of FAH won't be covered again here.Since the operation time of a vessel is mainly determ ined by the QC number assigned to it，the departure time adjustment of the vessel here is achieved by the change of its QC numbers.The number of QCs assigned to vessels can be changed to make sure that vessels can leave as soon as possible when they finish loading and unloading operations.In FAH，QC numbers searched in the upper layer are considered as suggested ones.

Fig.6 Scenarios after HAH adjustment

Fig.7 Flow chart of FAH

3 Com putational Experiments

In this section，the proposed method in section 2 is implemented through C#(VS2010)and the solutions are compared with the exact solution produced by Cplex 12.5.All computations are run on a PC with Intel T5750 and 4.0 GHz CPU.Moreover，the performance test and analysis of the algorithm for the integrated scheduling of BACAP with TTW problem are conducted in this section.

3．1 Test instance generation

The proposed algorithms(both HAH and FAH)are all applied to a suite of instances(Table 1)introduced by Meisel and Bierw irth in Ref.［5］.Three sets of test instances containing 20，30 and 40 vessels with 10 instances each are generated.

Table 1 Parameters of different vessel classes

In Table 1，U designates a uniform distribution of integer values in the specified interval.Within each instance，60%of the vessels belong to the feeder class，30%belong to the medium class，and 10%belong to the jumbo class.The population size of the three scale problem is 50，75，and 100，respectively.The maximum number of iteration is set to 200and the threshold of regenerating manipulation is 15，which follow the parameters setting in Ref.［12］.

Moreover，the length of quay is set to 100 m and the total number of available QC is 10.The planning horizon is one week，i．e．，168 h.

3．2 Experimental results

Tables 2－4 show the optimal solution(BC)by Cplex and the best solution of the proposed iterative algorithm(BI)with respect to tide(TTW=3)for different vessel scales.Besides，BRis the solution obtained by Cplex after right shifting manipulation.G1and G2can be calculated with formula G1= (BI－BC)/BCand G2=(BI－BR)/BR.

Table 2 Results for 20－vessel instances

Table 3 Results for 30－vessel instances

Table 4 Results for 40－vessel instances

Values with mark“LB”in Tables 3 and 4 is logged by Cplex after running 10 h and“Average”means the average value of 10 instances.Values(marked with“*”)are the average values of the available results.

3．3 Analysis of experimental results

The average values of G1in Tables2 and 3 are 0.00%and 1.40%，respectively.This indicates that solutions achieved by the proposed algorithm in this paper are very close to the exact solutions obtained by Cplex，when the vessel's scale is small and medium.Both reasonable gaps are effective enough to declare the validity of the proposed iterative algorithm.

Meanwhile，only 3 optimal solutions of the 10 instances can be gained by Cplex，when it comes to large－scaled problem (40 vessels).But the iterative algorithm developed in this paper can help to get near－optimal solutions of all instances in an acceptable period of time.Therefore，it is a powerfulmethod to solve the large－scaled problem.In addition，the gap 8.14% seems larger than ever achieved，but the real gap will be definitely smaller than this if the optimal solution can be achieved.

In actual operational management，operators are accustomed to do the scheduling firstly and then do some empirical adjustment when tide comes.In consideration of above，G2in Tables 1 and 2 signify the improvementmade by the proposed algorithm for the integrated scheduling problem with tidal factor comparing with traditional handling method. Therefore，taking the tidal factor into consideration ahead in the integrated scheduling problem is a preferable choice.

3．4 Com parison of the calculating time

The mean calculating time of 10 instances for different scales is shown in Table 5.

Table 5 Calculating time

In this table，tCand tIare values of the two solving methods，respectively.G3can be calculated with the formula G3=(tI－tC)/tC.

The value of G3shows that the proposed algorithm performs much better at the aspect of the computational efficiency for small－and medium－scaled problems.In terms of results'procurability for large－scaled problems，applying the iterative algorithm to the integrated scheduling problem can easily achieve a reasonable solution in an acceptable period time.

4 Conclusions

This paper builds an integrated scheduling model to represent BACAP in tidal port.To address this problem，two heuristics are developed to comprise an iterative algorithm based on genetic algorithm.Thiswell－designed algorithm is proved to be valid comparing with the results gained by Cplex. Furthermore，it has a great advantage over the commercial solver as the scale of vessels expands.

Surely，the solution of this integrated scheduling problem with tidal effectbenefits a lot from several idealized assumptions in this paper，such as the rigidly cyclical tide，deliberately ignored transferring time of QC and the selectively discarded uncertainty.All these above will be interesting research directions in the future.

［1］Xu D S，Li C L，Leung J Y T.Berth Allocation with Time－Dependent Physical Lim itations on Vessels［J］.EuropeanJournal of Operational Research，2012，216(1):47－56.

［2］Barros V H，Costa T S，Oliveira A C M，et al.Model and Heuristic for Berth Allocation in Tidal Bulk Ports with Stock Level Constraints［J］.Computers＆Industrial Engineering，2011，60(4):606－613.

［3］Park Y M，Kim K H.A Scheduling Method for Berth and Quay Cranes［J］.OR Spectrum，2003，25(1):1－23.

［4］Imai A，Chen H C，Nishimura E，etal.The Simultaneous Berth and Quay Crane Allocation Problem［J］.Transportation Research Part E:Logistics and Transportation Review，2008，44 (5):900－920.

［5］Meisel F，Bierw irth C.Heuristics for the Integration of Crane Productivity in the Berth Allocation Problem［J］.Transportation Research Part E:Logistics and Transportation Review，2009，45 (1):196－209.

［6］Raa B，DullaertW，Schaeren R V.An Enriched Model for the Integrated Berth Allocation and Quay Crane Assignment Problem［J］.Expert Systems with Applications，2011，38(11):14136－14147.

［7］Yang C X，Wang X J，Li Z F.An Optimization Approach for Coupling Problem of Berth Allocation and Quay Crane Assignment in Container Term inal［J］.Computers＆Industrial Engineering，2012，63(1):243－253.

［8］Liang X L，Li W F，Zhao W，et al.Multistage Collaborative Scheduling of Berth and Quay Crane Based on Heuristic Strategies and Particle Swarm Optim ization［C］.2012 IEEE 16th International Conference on Computer Supported Cooperative Work in Design(CSCWD)，Wuhan，China，2012:913－918.

［9］Elwany M H，Ali I，Abouelseoud Y.A Heuristics－Based Solution to the Continuous Berth Allocation and Crane Assignment Problem［J］.Alexandria Engineering Journal，2013，52(4):671－677.

［10］Hu Q M，Hu Z H，Du Y Q.Berth and Quay－Crane Allocation Problem Considering Fuel Consumption and Em issions from Vessels［J］.Computers＆Industrial Engineering，2014，70:1－10.

［11］Chang D F，Jiang Z，Yan W，et al.Integrating Berth Allocation and Quay Crane Assignments［J］.Transportation Research Part E:Logistics and Transportation Review，2010，46(6):975－990.

［12］Gui X Y，Lu Z Q，Han X L.Integrating Optim ization Method for Continuous Berth and Quay Crane Scheduling in Container Terminals［J］.Journal of Shanghai Jiaotong University: Science，2013，47(2):226－229.

［13］Bierw irth C，Meisel F.A Survey of Berth Allocation and Quay Crane Scheduling Problems in Container Term inals［J］. European Journal of Operational Research，2010，202(3):615－627.

［14］Türkoˇgullar?Y B，Ta?k?n Z C，A ras N，et al.Optimal Berth Allocation and Time－Invariant Quay Crane Assignment in Container Term inals［J］.European Journal of Operational Research，2014，235(1):88－101.

［15］Lu Z Q，Han X L，Xi L F.Simultaneous Berth and Quay Crane A llocation Problem in Container Term inal［J］.Advanced Science Letters，2011，4(6/7):2113－2118.

TP29

A

1672－5220(2015)04－0549－05

date:2014－05－16

s:National Natural Science Foundations of China(Nos.61273035，71471135)

*Correspondence should be addressed to ZHOU Bing－h(huán)ai，E－mail:bhzhou@#edu.cn