- Email: [email protected]

Contents lists available at ScienceDirect

Transportation Research Part E journal homepage: www.elsevier.com/locate/tre

A robust optimization approach to the integrated berth allocation and quay crane assignment problem Xiao Ting Shang a, Jin Xin Cao b,⇑, Jie Ren a a b

School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, PR China Institute of Transportation Engineering, Inner Mongolia University, Hohhot 010070, PR China

a r t i c l e

i n f o

Article history: Received 18 September 2015 Received in revised form 20 June 2016 Accepted 28 June 2016

Keywords: Berth allocation Quay crane assignment Robust optimization Price constraints Genetic algorithm Insertion heuristic algorithm

a b s t r a c t This paper investigates the integrated berth allocation and quay crane assignment problem in container terminals. A deterministic model is formulated by considering the setup time of quay cranes. However, data uncertainties widely exist, and it may cause the deterministic solution to be infeasible. To handle the uncertainties, a robust optimization model is established. Furthermore, to control the level of conservativeness, another robust optimization model with the price constraints is proposed. A genetic algorithm and an insertion heuristic algorithm are suggested to obtain near optimal solutions. Computational experiments indicate that the presented models and algorithms are effective to solve the problems. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction In the background of economic globalization, transportation logistics is rapidly developing. Container transportation industry has become to a mainstay of the transportation and logistics, and container terminals play an indispensable role in the multi-modal transportation network. Major operational problems faced in seaside operations are comprehensively described by Vis and Koster (2003), Steenken et al. (2004), and Stahlbock and Voß (2008), which can be typically identified as the berth allocation problem (BAP), the quay crane assignment problem (QCAP), and the quay crane scheduling problem. The BAP and the QCAP are two decision problems that are highly dependent on each other, the quay crane (QC) capacity demand is influenced by the berthing position, and the handling time varies with the QC assignment. Berths and QCs are two types of the scarce and expensive resources, whereas they are also the key resources in the BAP and the QCAP. Therefore, it is an urgent problem for container terminals to reasonably configure and to effectively utilize these resources. In some practical situations, the QC setup time cannot be negligible in the decision-making process due to the QC speed limitations of shifting along the quay. These situations include emergency situations, priority level, long shoreline, complex layout, etc. Besides, the QC setup time influences the starting and ending time of service, so it is influential on the QC assignment and may lead to the infeasible schedule. Hence, taking the QC setup time of shifting along the quay into account, a deterministic model is proposed in this paper for the integrated berth allocation and quay crane assignment problem (BACAP).

⇑ Corresponding author. E-mail address: [email protected] (J.X. Cao). http://dx.doi.org/10.1016/j.tre.2016.06.011 1366-5545/Ó 2016 Elsevier Ltd. All rights reserved.

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

45

Fig. 1. An illustrative example for the BACAP.

In addition, the operational environment of QC is in the open air, so the QC productivity is affected by a variety of external factors, including ambient temperature, humidity, degree of salt spray, icing, wind power, etc. These external factors occur frequently, which often make the planned schedule infeasible. Temporary adjustments of scheduling will decrease the efficiency and increase the cost. Therefore, it is necessary to develop approaches capable of handling these data uncertainties to improve the robustness of plan. This paper proposes two robust optimization models to solve the BACAP under data uncertainties in QC productivities. The BACAP can be depicted as a space-time diagram, as is shown in Fig. 1. In our paper, it divides the time in periods of 1-h and the quay in segments of 10-m length. The height of each rectangle corresponds to the length of a vessel and the width corresponds to the required handling time. The lower-left vertex of a rectangle shows the vessel’s berthing position and the starting handling time. A gray rectangle represents a QC assignment to the vessel, and the number of gray rectangles corresponds to the number of QC assignment. From the computational complexity point of view, the BACAP is quite complex, since the BAP is an NP-hard problem (Park and Kim, 2003). It is more difficult to deal with the BACAP under data uncertainties. To effectively solve the large scale BACAP, it is necessary to develop effective solution algorithms. In this paper, considering the QC setup time of shifting along the quay, we propose a deterministic model for the BACAP. Recognizing that the common uncertainties may cause the deterministic solution to be infeasible, two robust optimization models are proposed to deal with data uncertainties. Due to the computational complexity, we develop a genetic algorithm (GA) and an insertion heuristic algorithm to obtain near optimal solutions.Computational experiments show that the solution quality and the computational time of the GA and the insertion heuristic algorithm are all comparable, even for the large-scale instances. Comparing these two algorithms, the GA is more competitive in term of solution quality, while the insertion heuristic algorithm is time-efficient in term of computational time. In addition, the proposed robust optimization method is risk-averse so that the solution will have robustness to uncertainties. We can adjust the price constraint to control the robustness of the model, thus the failure probability of plan can be decreased. The following of this paper is organized as follows. In Section 2, related literatures are reviewed. Problem description and model formulation are given in Section 3. A GA and an insertion heuristic algorithm are developed in Sections 4 and 5, respectively. To test the efficiency of the proposed method, computational experiments are conducted in Section 6. Finally, Section 7 concludes this paper.

2. Literature review The BAP plays an important role in container terminal operations. According to the layout of berths, the BAP can be classified into the discrete BAP and the continuous BAP. Lai and Shih (1992) proposed a discrete BAP inspired by the usage of more efficient terminal in Hong Kong International Terminal, and the problem was solved by a heuristic algorithm. Imai et al. (1997) studied the statistic discrete BAP to minimize the sum of the waiting time and the handling time taking the priority order into consideration. Then Imai et al. (2001) proposed a dynamic BAP with the same objective function as Imai et al. (1997), and the problem was solved by a Lagrange relaxation algorithm. Kim and Moon (2003) proposed a mixed-integer-linear programming model for the continuous BAP to minimize the penalty

46

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

cost due to delays in the departures of vessels and the additional handling costs due to non-optimal locations of vessels. A simulated annealing algorithm was developed to obtain the optimal solution of the problem. Cordeau et al. (2005) studied a dynamic discrete BAP, and the problem was formulated as a multi-depot vehicle routing problem with time window. Imai et al. (2007) addressed the BAP in a multi-user container terminal with indented berths, and solved the problem by GA. Lee et al. (2010) studied the dynamic continuous BAP to minimize the sum of weighted turnaround time for each incoming vessel. Two versions of greedy randomized adaptive search procedures were developed to search for near optimal solutions. The QCAP determines the handling time or cost for given berthing position. Dagnazo (1989) examined the QCAP to minimize the aggregated costs of delay with the assumptions that vessels were divided into holds and that only one QC can work on a hold at a time. Kim and Park (2004) proposed a mixed-integer programming model for the QCAP and developed a branch and bound method to obtain the optimal solution and a heuristic search algorithm to overcome the computational difficulties. Chen et al. (2011) addressed a mixed integer programming model for the group-based QCAP at indented berth to minimize the makespan of tasks. The problem took two unique features into account: one is that QCs are able to serve a vessel simultaneously from both sides of the vessel at indented berth and the other is that the handling time of a task depends on the assigned QCs. To solve the problem, a decomposition heuristic framework was developed and was enhanced by Tabu search. The BAP and the QCAP are two highly dependent problems in the operational process. It was Park and Kim (2003) that firstly studied the BACAP. They proposed an integer programming model for the BACAP based on the assumption that the berthing duration of a vessel was inversely proportional to the number of QCs assigned to the vessel. The model was solved by a two-phase solution procedure combined by the subgradient optimization technique and the dynamic programming technique. Imai et al. (2008) studied the BACAP with the minimization of the total service time. In their study, the number of QCs in the course of handling was allowed to vary if QC transfer took place. Due to the complexity of the problem, a GA-based heuristic was proposed to find a near-optimal solution. Meisel and Bierwirth (2009) investigated the BACAP considering several important real-word influential factors, such as the decrease of marginal productivity of QCs assigned to a vessel, and the increase of handling time because of the deviation between the desired and the actually chosen berthing position. A construction heuristic, a local refinement procedure, and two meta-heuristics were provided to solve the problem. Lee and Wang (2010) provided a mixed integer programming model for the BACAP and designed a GA to obtain near optimal solutions. Rodriguez-Molins et al. (2012) proposed a mixed integer linear programming model to minimize the total weighted service time of incoming vessels, and developed a GA to obtain near-optimal solutions. Lalla-Ruiz et al. (2014) proposed an alternative to solve the tactical BAP involved process (berth and quay crane allocation) jointly to minimize the housekeeping costs and to maximize the total value of the QC profiles assigned to the ships. To obtain good quality solutions with considerably short computational effort, a biased random key GA was developed for solving this problem. Turkogullari et al. (2016) studied the berth allocation, time-variant QC assignment and scheduling problem with QC setup time, taking the costs of deviation from the desired berthing position into the objective function. To solve the model, the cutting plane procedure was proposed. Iris et al. (2015a) proposed novel generalized set partitioning formulations solving the problem introduced by Meisel and Bierwirth (2009), which considered both time-variant and time-invariant QC assignment policies. Li et al. (2015) developed a nonlinear mixed-integer programming model for the continuous berth allocation and specific QC assignment problem with QC coverage range, and applied a novel heuristic method based on spatiotemporal conflict analysis to derive a high-quality berth-QCs plan. Extensive reviews on BACAP and its variants can be found in Bierwirth and Meisel (2015) and Iris et al. (2015a). The above studies did not take the uncertainties of the real-world into consideration. Recent years, some researchers have noticed the importance of uncertainties. Zhou and Kang (2008) proposed a berth and QC allocation model under stochastic environments, so as to minimize the average waiting time of containership in terminals. To obtain a good solution with considerably small computational effort, a GA is developed with a reduced search space. Karafa et al. (2013) formulated the BAP with stochastic vessel handling times as a bi-objective function. To solve the problem, an evolutionary algorithm-based heuristic and a simulation-based Pareto front pruning algorithm was proposed. Golias et al. (2014) studied the discrete BAP where vessel arrival and handling time were not known with certainty. A bi-objective optimization problem was formulated and a heuristic algorithm was proposed to solve the robust BAP. Rodriguez-Molins et al. (2014) proposed a multi-objective model to solve the BACAP under uncertainty which had two conflicting objectives: minimizing the total service time and maximizing the robustness or buffer time. A GA was developed for solving this multi-objective optimization problem. Zhen (2015) studied the tactical-level berth allocation under uncertain vessel’s operation time, and proposed two types of stochastic models to handle the uncertainty. To solve the problem, some meta-heuristic algorithms were proposed. Iris et al. (2015b) presented a stochastic programming approach for the quayside operation problem under uncertainties, which focused on the BAP, QC assignment and scheduling. To solve the problem, a decomposition algorithm was designed. To the best of our knowledge, the number of studies that consider the QC setup time of shifting along the quay for the BACAP and the BACAP under data uncertainties is limited. However, we believe that it is necessary to take these two aspects into the BACAP. Thus, this paper proposes a deterministic model considering the QC setup time of shifting along the quay for the BACAP, and provides two robust optimization models to handle data uncertainties in QC productivity. To obtain solutions of this problem, a GA and an insertion heuristic algorithms are proposed.

47

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

3. Problem description and model formulation 3.1. Assumptions Assumptions for the problem under study are summarized as follows: 1. 2. 3. 4. 5. 6. 7.

The arrival times of vessels are known, and the starting times are later than the arrival times. For the convenience of calculation, the continuous quay and time are partitioned by a fixed length. Once the vessel is in the berth, it does not move until all tasks are completed. The berth is continuous and the length of the berth is larger than the length of any vessel to be served. There is no difference in productivity among QCs in the deterministic scenario. Each vessel has its maximum and minimum number of QCs simultaneously. Once QC operation starts on a vessel, it cannot stop until all tasks are finished.

3.2. Deterministic model The following notations are used in the model formulation. Input V Q L T T1 T0 li 0 bi mi r min i r max i Ri ai

data set of vessels to be served, V ¼ f1; 2; . . . ; Ng, where N is the number of vessels in the container terminal number of QCs that can work simultaneously number of 10-m berth segments (length of the quay) set of 1-h periods, T ¼ f0; 1; . . . ; H 1g, where H is the number of planning periods of 1-h set of f0; 1; . . . ; H 2g set of periods that the QC productivity is subject to data uncertainties length of vessel i 2 V, given as a number of 10-m segment desired berthing position of vessel i 2 V QC capacity demand of vessel i 2 V, given as a number of QC-hours minimum number of QCs agreed to serve vessel i 2 V simultaneously maximum number of QCs allowed to serve vessel i 2 V simultaneously feasible range of QCs assignable to vessel i 2 V, Ri ¼ r min ; r max i i arrival time of vessel i 2 V l weighted coefficient, l 2 ½0; 1 a interference exponent b berth deviation factor M a large positive number Decision variables bi integer, berthing position of vessel i 2 V si integer, time of starting the handling of vessel i 2 V ei integer, time of ending the handling of vessel i 2 V r it binary, equal to 1 if at least one QC is assigned operating to vessel i 2 V at time t 2 T and 0 otherwise r itq binary, equal to 1 if exactly q QCs are assigned operating to vessel i 2 V at time t 2 T; q 2 Ri and 0 otherwise Dbi integer, deviation between the desired and actually chosen berthing position for vessel i 2 V; Dbi ¼ b0i bi integer, the number of QCs being set up to serve vessel i 2 V at time t 2 T 1 wit uit binary, equal to 1 if there is at least one QC is under setting up to serve vessel i 2 V at time t 2 T 1 and 0 otherwise yij binary, equal to 1 if vessel i is berthed in front of vessel j, i.e., bi þ li 6 bj , i – j; 8i; j 2 V and 0 otherwise zij binary, equal to 1 if handling of vessel i ends not later than handling of vessel j starts, 8i; j 2 V; i – j and 0 otherwise

The handling time of a vessel is not only influenced by the berth deviation, but is also influenced by the interferences of the assigned QCs considerably. According to Schonfeld and Sharafeldien (1985), for a given interference exponent a, the productivity obtained from assigning q QCs to a vessel for one hour is qa QC-hours. Given the berth deviation factor b, the required QC-hours are ð1 þ b Dbi Þ mi , and the minimum handling time needed is given as min

di

¼

& ’ ð1 þ b Dbi Þ mi ; max a ri

8i 2 V:

ð1Þ

48

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

The deterministic model of the BACAP with QC setup times is formulated as ½P 1 :

½P 1

min

Z1 ¼

X ½lðei si Þ þ ð1 lÞðsi ai Þ

ð2Þ

i2V

s:t:

XX qa r itq P ð1 þ b Dbi Þ mi

8i2V

ð3Þ

t2T q2Ri

X XX q r itq þ wit 6 Q i2V q2Ri

8t 2 T 1

ð4Þ

i2V

XX q r iðH1Þq 6 Q

ð5Þ

i2V q2Ri

X r itq ¼ r it

8i 2 V; t 2 T

ð6Þ

q2Ri

" # X X q r iðtþ1Þq q ritq P 1 M ð1 uit Þ þ q2Ri

8i 2 V;

t 2 T1

ð7Þ

q2Ri

" # X X q r iðtþ1Þq q r itq P 0 M uit q2Ri

8i 2 V;

t 2 T1

ð8Þ

q2Ri

M ð1 uit Þ þ wit P

X X q r iðtþ1Þq q r itq q2Ri

X

M ð1 uit Þ þ wit 6

8i 2 V;

t 2 T1

ð9Þ

q2Ri

q r iðtþ1Þq

q2Ri

X q r itq

8i 2 V;

t 2 T1

ð10Þ

q2Ri

M uit þ wit 6 0 8i 2 V; X r it ¼ ei si 8i 2 V

t 2 T1

ð11Þ ð12Þ

t2T

8i 2 V;

ð1 þ tÞ r it 6 ei

t2T

ð13Þ

8i 2 V;

t rit þ H ð1 rit Þ P si

t2T

ð14Þ

0 bi

8i 2 V

ð15Þ

Dbi P bi bi

8i 2 V

ð16Þ

Dbi P bi 0

8i; j 2 V;

bj þ M ð1 yij Þ P bi þ li

i–j

ð17Þ

sj þ M ð1 zij Þ P ei

8i; j 2 V;

i–j

ð18Þ

yij þ yji þ zij þ zji P 1

8i; j 2 V;

i–j

ð19Þ

bi 2 f0; 1; . . . ; L li g 8i 2 V

ð20Þ

si 2 fai ; . . . ; H 1g; ei 2 fai ; . . . ; Hg 8i 2 V

ð21Þ

wit 2 f0; 1; . . . ; Q g uit 2 f0; 1g

8i 2 V;

8t 2 T 1 ;

r itq ; r it ; yij ; zij 2 f0; 1g

t 2 T1

ð22Þ

i2V

8t 2 T;

ð23Þ q 2 Ri ;

i; j 2 V;

i–j

ð24Þ

Eq. (2) is the objective function of this model which aims to minimize the total weighted handling time and the waiting time for all vessels within the planning horizon given weighted coefficient l. Constraints (3) make sure that every vessel receives the required QC capacity with respect to productivity losses by QC interference and the berth deviation. Constraints (4) and (5) ensure that at most Q QCs are utilized during a period. Constraints (6) ensure the consistent setting of the corresponding variables rit and ritq . The QC setup time of shifting along the quay is taken into consideration by Constraints (7)–(11). They ensure that when there is a QC available and assigned to the vessel at time t 2 T 1 , the QC is able to start handling at time t þ 1. Constraints (7) and (8) ensure the relationship between uit and ritq . Constraints (9) and (10) give the equivalent relationship of wit and ritq if uit ¼ 1 and Constraints (11) force wit ¼ 0, when uit ¼ 0. Constraints (12)–(14) determine the starting time and ending time for handling without priority. Constraints (15) and (16) are to obtain the deviations from the desired berthing position. Constraints (17)–(19) are to avoid vessels overlapping at berth and time space. Constraints (20)–(24) give the domain of all decision variables.

49

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

3.3. Robust optimization model Model ½P 1 is a deterministic model taking the QC setup time of shifting along the quay into account. However, the QC productivity variations owing to data uncertainties may lead to the planned schedule infeasible. To withstand data uncertainties in QC productivities, a robust optimization model is proposed. It is considered that the mean productivity obtained from assigning q QCs to a vessel for one hour is qa QC-hours. The maximum productivity is q QC-hours, and the minimum productivity is 0 QC-hours. So it is reasonable that we assume that the productivity shows the symmetric distribution with ^ ¼ minfqa ; q qa g. So let q ~ denote the actual productivity of a group of q QCs under uncerthe mean value qa and the error q ~ takes value in ½qa q ^; qa þ q ^. We introduce the auxiliary variable xitq ; 8t 2 T; q 2 Ri ; i 2 V, which is greater tainties, where q than or equal to the decision variable r itq and is used to ensure the highest robustness. Constraints (3) are replaced to the robust optimization Constraints (25) and (26) which are capable of absorbing data uncertainties.

XX XX ^ xitq P ð1 þ b Dbi Þ mi ; q qa r itq

8i 2 V

ð25Þ

t2T 0 q2Ri

t2T q2Ri

xitq P ritq ;

8t 2 T;

q 2 Ri ;

i2V

ð26Þ

Then we have the robust optimization model ½P 2 :

½P 2

min Z 2 ¼

X

½lðei si Þ þ ð1 lÞðsi ai Þ

ð27Þ

i2V

s:t: ð4Þ—ð26Þ: 3.4. Robust optimization model with the price constraints Robust optimization addresses the BACAP under data uncertainties by guaranteeing the feasibility of the solution for the worst case of the parameters. However, feasibility often comes at the cost of the performance, the robust optimization often has much conservativeness. To address the over-conservativeness, another robust optimization model with the price constraint was firstly presented by Bertsimas and Sim (2004), which is linear without excessively affecting the objective function value. The price constraint Ci is introduced to adjust the level of conservativeness of the solutions, which is not necessarily integer and takes value in the interval ½0; jT 0 j, where jT 0 j is the cardinality of set T 0 . Speaking intuitively, up to ^. Let gi denote the maximum bCi c of these coefficients are allowed to change, and one coefficient changes by ðCi bCi cÞ q data uncertainties determined by the price constraint Ci ; Si denote a set whose cardinality is equal to bCi c, and vi denote a time belongs to set T 0 n Si ; 8i 2 V. Based on a general case presented by Bertsimas and Sim (2004), the maximum data uncertainties gi for vessel i determined by the price constraint Ci can be expressed as Eq. (28):

( ) XX ^ ^ q r itq þ ðCi bCi cÞq r ivi q : gi ðritq ; Ci Þ ¼ S max fSi fvi gjSi # T 0 ;jSi j¼bCi c;vi 2T 0 nSi g t2Si q2Ri

ð28Þ

Constraints (25) can be transformed to Constraints (29) to obtain a robust optimization with the price constraints:

XX qa r itq t2T q2Ri

( ) X XX ^ xitq þ ^ xivi q P ð1 þ b Dbi Þ mi ; q max ðCi bCi cÞq S fSi fvi gjSi # T 0 ;jSi j¼bCi c;vi 2T 0 nSi g t2Si q2Ri q2Ri

8i 2 V:

ð29Þ

3.5. Reformulation n

In order to obtain a linear optimization model, we need to reformulate nonlinear Constraints (29). Given a solution o 8t 2 T; q 2 Ri ; i 2 V , we introduce the auxiliary variable g itq ; 8t 2 T; q 2 Ri ; i 2 V, so Eq. (30), that is the protection

ritq j

function of the ith element of Constraints (29)

gi ritq ; Ci ¼

fSi

S

max

fvi gjSi # T 0 ;jSi j¼bCi c;vi 2T 0 nSi g

( ) XX ^ ritq þ ðCi bCi cÞq ^ riv q q i

ð30Þ

t2Si q2Ri

equals to the following linear optimization problem:

gi ritq ; Ci ¼ max

( ) XX ^ q r itq g itq 0

ð31Þ

q2R

t2T X X i s:t: g itq 6 Ci

ð32Þ

t2T 0 q2Ri

0 6 g itq 1;

8t 2 T 0 ;

q 2 Ri :

ð33Þ

50

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

We introduce the dual variables pitq ; v iq ; 8t 2 T 0 ; q 2 Ri ; i 2 V. By the strong duality, the dual of Eqs. (31)–(33) can be expressed by Eqs. (34)–(36), and the objective values of these two forms coincide.

XX

min

pitq þ Ci

t2T 0 q2Ri

v iq

ð34Þ

q2Ri

X

X X ^ r itq ; pitq P q

q2Ri

q2Ri

v iq þ

s:t:

X

X

8t 2 T 0 ;

i2V

ð35Þ

q2Ri

pitq P 0;

q2Ri

X

v iq P 0; 8t 2 T 0 ;

i 2 V:

ð36Þ

q2Ri

Therefore, the equivalent linear formulations of Constraints (29) are as Constraints (35)–(37).

X XX XX qa r itq Ci v iq pitq P ð1 þ b Dbi Þ mi ; t2T q2Ri

8i 2 V

ð37Þ

t2T 0 q2Ri

q2Ri

Then we have the linear robust optimization model with price constraints ½P 3 :

½P 3

min Z 3 ¼

X ½lðei si Þ þ ð1 lÞðsi ai Þ

ð38Þ

i2V

s:t: ð4Þ—ð24Þ; ð26Þ; ð36Þ—ð37Þ X X X v iq þ pitq P q^ xitq ; q2Ri

q2Ri

8t 2 T 0 ;

i2V

ð39Þ

q2Ri

3.6. Probability bounds of constraint violation As mentioned in Section 3.4, if up to bCi c of the coefficients change within their bounds, and one coefficient qa changes by ^, the solutions of model ½P 3 will remain feasible. In this subsection, we show that even if more than bCi c change, ðCi bCi cÞ q the robust solutions of model ½P 3 will remain feasible with high probability. ^Þ=q ^ 2 ½1; 1; 8i 2 V; t 2 T 0 , which follow an unknown symmetric distribuWe define the random variables nit ¼ ðqa q tion. The following proposition, based on a general proposition proposed by Bertsimas and Sim (2004), gives the probability bounds of constraint violation, and proves the robustness of model ½P 3 . To present a more concise conclusion, the quantities cit ; 8i 2 V; t 2 T 0 are introduced in the following proposition. n o Proposition. Let ritq j8t 2 T; q 2 Ri ; i 2 V be an optimal solution of model ½P 3 . Let Si and vi be the set and the time respectively that achieve the maximum for Eq. (30). The probability that the ith element of Constraints (25) is violated satisfies:

P

XX ~ ritq < ð1 þ b Dbi Þ mi q

cit ¼

6P

XX t2T 0 q2Ri

t2T q2Ri

where

!

8 > < 1; P

!

cit nitq P Ci 6 exp

C2i

2jT 0 j

! ð40Þ

if t 2 Si ^r q

itq q2Ri 0 > : P q^r ; if t 2 T n Si q2Ri

ik q

and

k ¼ arg

^ r ikq : min S q fvi g

k2Si

Proof. See Appendix A for proof. 4. Genetic algorithm GA is a fairly complete method to handle combinational problems. It adopts the model of nature evolution, including selection, crossover and mutation. In recent years, many researchers have utilized GA to solve the BACAP, for example, Rodriguez-Molins et al. (2012) designed a GA for the deterministic BACAP and Rodriguez-Molins et al. (2014) developed a GA for the robust BACAP under uncertainties.

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

51

Algorithm 1 gives the structure of the GA. At the beginning, population initialization produces a series of random solutions and calculates the fitness of each individual, called the first generation. If the first generation cannot satisfy the termination criterion, a new generation is created by the processes of selection, crossover and mutation. When the termination criterion is satisfied, the searching process of GA ends. Algorithm 1. Genetic algorithm. 1. 2. 3. 4. 5. 6. 7. 8. 9.

Randomly generated the initial population; Evaluate the initial population; while the termination criterion is not satisfied do Selection; Crossover; Mutation; Evaluate the chromosomes generated by these reproduction operators; end while return the schedule from the best chromosome evaluated so far;

The following of this section gives the concrete procedures of the GA. 4.1. Chromosome representation In this study, a multi-level encoding scheme with non-negative integers is applied to represent possible solutions. The first row of the chromosome represents the permutation of vessels, the second row represents the berthing positions, and the third row represents the minimum number restriction of QC assignment r i which takes value in the range r min ; r max . i i Note that, a vessel may have different number of QCs during the handling time. So the number of QCs each period takes value randomly. An example of chromosome structure is shown in Fig. 2, the permutation of vessels is in the range r i ; rmax i ð5; 4; 3; 6; 1; 2Þ, the berthing positions are ð4; 9; 1; 15; 23; 12Þ, and the minimum number restriction of QC assignment is ð3; 2; 2; 1; 2; 3Þ, respectively. 4.2. Initialization The GA begins with a group of initial individuals as the first generation. To simplify the generation of initial individuals, a series of randomly generated chromosomes are used as the first generation. The permutation of vessels is a random sequence. The berthing position is randomly selected in the range ½0; L li . The minimum number restriction of QC assign ment is randomly generated in the range rmin ; rmax . All the chromosomes of the first generation are generated in accordance i i with this rule, and the initialization ends when the number of chromosomes reaches the population size. 4.3. Decoding and fitness evaluation Considering the structure of the chromosome, the permutation of vessels is used as dispatching rule. We use the following decoding algorithm: the genes are visited from left to right in the chromosome sequence. Due to the setup time of QCs, the judgment is performed whether there are ri available QCs at time interval ½si 1; ei 1Þ. For each gene ði; bi ; ri Þ, if the berthing position bi during ½si ; ei Þ is idle and ri QCs during ½si 1; ei 1Þ are available, vessel i is scheduled at the earliest starting time si . Otherwise, the starting time would experience a delay of one time period. In this way, none of the constraints will be violated. In case that QCs are scheduled to vessel i, the number of idle QCs decreases to ensure that no other vessel uses these QCs. It is also ensured that no other vessel moors berthing positions ½bi ; bi þ li Þ during the handling period of vessel i. Once a valid starting time si has been assigned to each vessel i, the waiting time and the handling time are obtained, accordingly. Let ci be a chromosome, and Z : ðci Þ be the objective value represented by the chromosome. The objective function

Fig. 2. An example of the chromosome structure.

52

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

of this study is to minimize the weighted sum of service time of all vessels, so we take the reciprocal of the objective function as the fitness evaluation function: fitðci Þ ¼ 1=Z : ðci Þ. 4.4. Selection According to certain selection probability, a new population is selected from the parent population, which is better than the parent population. The roulette wheel approach is adopted as the selection procedure in this study, which selects a new population with regard to the probability distribution based on fitness values (Gen and Cheng, 1996). A single chromosome is selected as a member of new population by the roulette wheel approach in the following way: Step 1: Generate a random number rd from the range ½0; 1. P Step 2: Calculate the fitness value for each chromosome fitðci Þ and the total fitness for the population i fitðci Þ. P Step 3: Calculate selection probability pi for each chromosome: pi ¼ fitðci Þ i fitðc i Þ. Pi Step 4: Calculate the cumulative probability hi for each chromosome: hi ¼ j¼1 pj . Step 5: If rd < h1 , then select the first chromosome; otherwise, select the ith chromosome such that hi1 < rd < hi .

4.5. Crossover Crossover plays an important role in GA. An individual of the new population is selected when the randomly generated number is smaller than the presupposed crossover rate pc , and then it exchanges the genes with another selected individual according to certain crossover rules. The common crossover methods are one-cut-point, two-cut-point, multi-cut-point, uniform, ordered crossover (Gen and Cheng, 1996). To avoid the occurrence of infeasible offspring, this paper adopts an ordered crossover, and the major procedures are the following: Step 1: Randomly select two crossover columns, k1 and k2 ð1 6 k1 6 k2 6 NÞ. Step 2: Each gene in Column k of chromosome Parent 1 and Parent 2, k1 6 k 6 k2 is copied to the same position in chromosome Offspring 1 and Offspring 2, respectively. Step 3: Delete the columns in Parent 2 and Parent 1 which present the genes of vessel V k of Parent 1 and Parent 2, k1 6 k 6 k2 . Step 4: Place the columns of Parent 2 and Parent 1 into the unfixed columns of Offspring 1 and Offspring 2 from left to right according to the order of the chromosome to produce offspring, respectively. An example of chromosome crossover is illustrated in Fig. 3. After the crossover, the fitness of the two parents and the two offspring need to be compared. If the fitness value of the offspring is worse than its parent, the offspring is replaced with its parent. 4.6. Mutation Mutation forces GA searching by altering one or more genes with a probability equal to the mutation rate pm . It is beneficial for GA to avoid having premature convergence and facilitating the discovery of global optimal solution. Owing to the non-negative integers and multi-level encoding approach, we take different levels mutation procedure for selected mutation columns. For the first row, swap mutation is used, which swap the genes of the selected columns. For the second row, half of chromosomes take the desired berthing position mutation, that is to say the original berthing position is replaced with the desired berthing position, and the other half take random mutation in accordance with the initialization, correspondingly. For the third row, half of chromosomes take the maximum number of QCs mutation, that is to say the original number of QCs is replaced with the maximum number of QCs, and the other half take random mutation in accordance with the initialization, correspondingly. An example of chromosome mutation is illustrated in Fig. 4, which takes random mutation for the second and third rows.

5. Insertion heuristic algorithm In this section, an insertion heuristic algorithm, inspired by the construction heuristic algorithm (Meisel and Bierwirth, 2009), is designed to solve the BACAP. The essence of this BACAP problem is to obtain the berthing positions bi ; 8i 2 V, the starting times si ; 8i 2 V, and the numbers of assigned QCs that can minimize the objective function. The algorithm is

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

53

Fig. 3. An illustration of the crossover operator.

based on First Come First Served and Longest Service Time regulation, that is to say the priority sequence ðV 1 ; V 2 ; . . . ; V N Þ is determined by the arrival time ai and the QC capacity demand mi . So the vessel with earlier arrival time is prior to be served. If two vessels have the same arrival time, the one with larger QC capacity demand is prior to be served. However, if the one with larger QC capacity demand fails to be inserted and the other one has feasible insertion, the priority sequence changes and the feasible insertion is the schedule for the vessel with less QC capacity demand. The pseudo code of the insertion heuristic algorithm is presented by Algorithm 2.

54

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

Fig. 4. An illustration of the mutation operator.

Algorithm 2. Insertion heuristic algorithm. n o min 1. Initialization: set zV i 1; j i; sV j aV j ; K 1 maxi2V di 2. while sV j < aV j þ K 1 do P P 3. Calculate the number of idle QCs at time sV j 1 : idlesV j 1 ¼ Q j1 q2Ri q r kðsV 1 Þq ; k¼1 j

4. if idlesv j 1 P rmin and there is available berth then j 5. Select the berth bv j with the minimum berth deviation; 6. QC assignment; 7. Calculate the costzV j ; 8. zV i min zV j ; zV i ; 9. else 10. while j < N 1 and ajþ1 ¼ aj and zV i ¼ 1 do 11. j j þ 1; 12. if idlesV j 1 P rmin and there is available berth then j 13. Select the berth bV j with the minimum berth deviation; 14. QC assignment; 15. Calculate the cost zV j ; 16. zV i zV j ; 17. The priority sequence ð ; V i ; . . . ; V j Þ is replaced with 18. ð ; V j ; V i ; Þ; 19. j N; t 2 T; 20. End if 21. end while 22. if there is no feasible schedule for vessel having the same arrival time then 23. j i; 24. end if 25. end if 26. sV j sV j þ 1; 27. end while 28. return the best the solution for vessel V i ;

In the initialization, the total service time is initially set to infinity, and the starting time is set to the arrival time, n o min . To select a berthing position, the berth bV i for inserting vessel V i . The constant K 1 is determined by K 1 ¼ maxi2V di is picked with the minimum berth deviation such that the overlapping conflict is resolved. In QC assignment, if the max QCs assigned to the vessel at time t; otherwise, there number of idle QCs idlet is greater than rmax V i , there would be r V i would be idlet QCs assigned to the vessel at time t. This continues until QC capacity demand is satisfied. Once the QC assignment is provided, an ending time for its handling is fixed, so the total service time of vessel V i is determined by Eq. (2). If the termination condition is satisfied, the procedure of inserting vessel V i ends, and the best solution is determined.

55

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65 Table 1 Example vessel data. i

li

mi

bi

0

r min i

r max i

ai

1 2 3 4

5 5 4 7

6 6 5 8

5 12 8 1

1 1 1 2

2 2 2 3

3 6 4 6

5.1. An example for the insertion heuristic algorithm To illustrate the algorithm, we take an example that there are 4 vessels at a container terminal with L ¼ 20; T ¼ 12; Q ¼ 6; b ¼ 0:1; a ¼ 0:9. The weighted coefficient is set to l ¼ 0:6 and the constant is set to K 1 ¼ 3. The sailing data of every vessel is shown in Table 1. Firstly, we consider the deterministic model ½P 1 . The schedule for vessel 1 ðV 1 Þ and vessel 3 ðV 2 Þ is shown in Fig. 5 (a), next vessel 4 ðV 3 Þ and vessel 2 ðV 4 Þ have to be inserted. According to the priority sequence, vessel 4 is considered prior to P P vessel 2. The number of idle QCs at time t ¼ sV 3 1 ¼ aV 3 1 is idlesV 3 1 ¼ Q j1 q2Ri q r kðsv 3 1Þq ¼ 2, which is equal to k¼1 r min V 3 . However, there is no available berth for vessel 4 as is shown in Fig. 5(b). Due to the same arrival time as vessel 4, vessel 2 is considered. The berthing position b2 ¼ 14 is selected with the minimum berth deviation Db2 ¼ 2, so the QC capacity demand increases from m2 ¼ 6 to ð1 þ b Db2 Þ ¼ 7:2. As depicted in Fig. 5(c), the QC assignment delivers r 262 ¼ r 272 ¼ r282 ¼ r 292 ¼ 1, and the QC assignment is sufficient owing to 20:9 4 ¼ 7:46. So the original priority sequence ð1; 3; 4; 2Þ is replaced with the new priority sequence ð1; 3; 2; 4Þ. The schedule for vessel 4 (now, V 4 ) continues. Because P P 0 idle6 ¼ Q 3k¼1 q rv k 6q ¼ 1 and idle7 ¼ Q 3k¼1 q r V k 7q ¼ 3; sV 4 ¼ 8; bV 4 ¼ bV 4 and the QC assignment is shown in Fig. 5 (d). So zV 4 ¼ zV 4 ¼ 0:6 3 þ 0:4 2 ¼ 2:6. Then if sV 4 ¼ 9; zV 4 ¼ 0:6 3 þ 0:4 3 ¼ 3:0. So zV 4 ¼ 2:6 and the algorithm for

the deterministic model terminates with Z 1 ¼ 9:8. Secondly, we consider the robust optimization model ½P 2 . Assuming T 0 ¼ f8; 9; 10; 11g; 8i 2 V. Vessel 1 and vessel 3 are without the influence of uncertainties, so the schedule is same with the deterministic model ½P 1 shown in Fig. 6(a). For vessel 2, the solution in the deterministic model is infeasible as a result of 20:9 2 þ ½20:9 ð2 20:9 Þ 2 ¼ 7:19 < 7:2. To satisfy the QC capacity demand, the feasible QC assignment with zV 3 ¼ 0:6 5 ¼ 3:0 is shown in Fig. 6(a). For vessel 4, the solution

Fig. 5. Example positioning of vessels for model ½P 1 .

Fig. 6. Example positioning of vessels for model ½P 2 .

56

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

Fig. 7. Example positioning of vessels for model ½P 3 .

of the deterministic model can not satisfy the QC capacity demand as a result of ½30:9 ð3 30:9 Þ 3 ¼ 7:13 < 8. The feasible schedule is shown in Fig. 6(b), and it is the returned optimal solution with zV 4 ¼ 3:2. The algorithm for the robust optimiza-

tion model ½P 2 terminates with Z 2 ¼ 11:0, which is larger than that of the deterministic model. Thirdly, we solve the robust optimization model with the price constraint ½P 3 assuming Ci ¼ 1; 8i 2 V. For vessel 2, nP o a maxt2T 0 ¼ 2 20:9 ¼ 0:134, and 20:9 4 0:134 ¼ 7:326 > 7:2. So the solution for vessel 2 is q2R2 ðq q Þ r itq nP o a r262 ¼ r 272 ¼ r282 ¼ r 292 ¼ 1 as is shown in Fig. 7(a), and the cost is zV 3 ¼ 2:4. For vessel 4, maxt2T 0 ¼ q2R4 ðq q Þ r itq 3 30:9 ¼ 0:312, and 30:9 3 0:312 ¼ 7:75 < 8. So the feasible schedule for vessel 4 is r383 ¼ r 393 ¼ r 3ð10Þ3 ¼ r3ð11Þ1 ¼ 1 as is shown in Fig. 7(b), and the cost is zV 4 ¼ 0:6 4 þ 0:4 2 ¼ 3:2. The algorithm for the robust optimization model with

the price constraints terminates with Z 3 ¼ 10:4, which lies in between the objective values of model ½P 1 and model ½P 2 , as a result of controlling of the conservativeness. 5.2. QC resource leveling and spatial shift

The insertion heuristic algorithm generates a feasible berth plan with respect to a given priority sequence of vessels. Vessels which are inserted early have good prospects to place at their desired position. To alleviate the preferential treatment of earlier inserted vessels, QC resource leveling and spatial shifts of vessels toward the berth’s border are performed (Meisel and Bierwirth, 2009). Moreover, the QC resource leveling and the spatial shift are applied to all problems. The first refinement procedure is QC resource leveling. Given the priority sequence ðV 1 ; V 2 ; . . . ; V N Þ, vessel V 1 is inserted once for every resource level within the range RV 1 . Each of these incomplete berth plans is completed by subsequently inserting the remaining vessels V 2 to V N using the insertion procedure without a restricting resource level. Due to the resource restriction for vessel V 1 , other vessels receive higher priority regarding the QC assignment in the completed berth plans. Possibly, the saved QC capacity has not been completely exhausted by these vessels and can therefore be reassigned to vessel V 1 . The process is continued for every vessel up to vessel V N1 . A further refinement is spatial shift. Given the priority sequence ðV 1 ; V 2 ; . . . ; V N Þ, vessel V i is inserted in the new berth one n o after the other from the set bV i 1; . . . ; bV i K 2i ; bV i þ 1; . . . ; bV i þ K 3i for given positive constant K 2i < bV i ; K 3i < L bV i lV i . Once vessel V i is inserted, the schedules of vessels subsequently are generated in accordance with the insertion heuristic algorithm. The procedure continues until vessel V N1 is inserted. Fig. 8(a) shows the solution of insertion heuristic algorithm

Fig. 8. Schedule of after spatial shift.

57

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

and vessel 1 will be shifted to illustrate the procedure for given constants K 21 ¼ 3; K 31 ¼ 0. So the new berth set is f4; 3; 2g. If bV 1 ¼ 4, the QC capacity demand increase to mV 1 ð1 þ 1 bÞ ¼ 6:6, and the QC assignment is r 132 ¼ 1; r 142 ¼ 1; r 152 ¼ 1; r 162 ¼ 1. Then vessels subsequently are inserted shown in Fig. 8(b), and the cost of all vessels is 9.2. If bV 1 ¼ 3, the berth deviation of vessel 1 increase to 2 while the berth deviation of vessel 3 and vessel 2 decrease to 0, and the QC assignment is shown in Fig. 8(c). If bV 1 ¼ 2, vessel 3 and vessel 2 keep unchanged, because the insertion heuristic algorithm makes them search for the berthing position with the minimum deviation, and QC assignment is shown in Fig. 8(d). So the optimal berthing position of vessel 1 is bV 1 ¼ 3. 6. Computational experiments In this section, a comprehensive set of computational experiments are conducted to test the performance of proposed method. All models are solved by using the CPLEX 12.6 solver. The GA and the heuristic algorithm are implemented in C++. All tests are run on a 64 core AMD Opteron at 2.10 GHz and 2.0 GB of RAM. All running times are measured in seconds. 6.1. Benchmark instances and parameter setting This paper uses the benchmark provided by Meisel and Bierwirth (2009) to demonstrate the power of our solution approach. To obtain the data of Meisel and Bierwirth (2009), the link ‘‘http://www.scm.bwl.uni-kiel.de/de/forschung/research-data” is used. These instances divided vessels into three categories: Feeder, Medium, Jumbo. In each instance, there are 60% of the vessels belong to Feeder, 30% of the vessels belong to Medium, and 10% of vessels belong to Jumbo. The benchmark consists of 30 instances, and contains ten instances of 20, 30, and 40 vessels, respectively. The planning horizon T is set to one week (168 h). The terminal data is L ¼ 100 (1000 m) and Q ¼ 10. The values of parameters used are set as follows:

The The The The

weighted coefficient l ¼ 0:6. interference exponent of QC productivity a ¼ 0:9. berth deviation factor b ¼ 0:01. set that the productivity is subject to data uncertainties T 0 ¼ f0; 1; 2; . . . ; 47g.

Because the property of GA relies on the genetic manipulation, the convergence performance is influenced by the parameters of GA. Fig. 9 shows the convergence performance of the GA, where Fig. 9(a) shows the GA convergence of a representative instance, and Fig. 9(b) shows the generations for convergence of the 10 instances. From the figure, it can be seen that the convergence performance varies with the crossover rate and mutation rate, and the case pc ¼ 0:9; pm ¼ 0:1 is better. In this paper, the experiments are implemented with crossover rate pc ¼ 0:9 and mutation rate pm ¼ 0:1.

a

70

b The generation for convergence

68

Objective function value

66 64 62 60 58 56 54

70 60 50 40 30 20

52

50

80

0

20

40

60

80

100

10

0

2

Generations Fig. 9. Convergence performance of the GA.

4

6

Instances

8

10

58

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

Table 2 Performance of the GA and the insertion heuristic algorithm on problem ½P.

UBbest /LBbest : the best upper/lower bounds are obtained from Meisel and Bierwirth (2009) and Iris et al. (2015a); IH: the upper bound given by the insertion heuristic algorithm. Z MB : the best solution of an instance for three heuristics of Meisel and Bierwirth (2009). TGA /TMB : the computational time of the GA/Meisel and Bierwirth (2009). Gap0 = (UBbest -LBbest )/LBbest , Gap1 = (GA-LBbest )/LBbest , Gap2 = (IH-LBbest )/LBbest , Gap3 = (CH-LBbest )/LBbest , Gap4 = (ZMB -LBbest )/LBbest .

6.2. Computational results and analysis 6.2.1. Comparison with the Meisel and Bierwirth (2009) and Iris et al. (2015a) The deterministic model ½P 1 and the insertion heuristic algorithm are developed on top of Meisel and Bierwirth (2009). However, ½P 1 uses different objective function, considers QC setup time of shifting along the quay and assumes that the starting time is known, which becomes to a new problem. In order to compare ½P 1 with Meisel and Bierwirth (2009) and Iris et al. (2015a), we set setup time to zero, replace arrival time with the earliest starting time and use the same objective function as Meisel and Bierwirth (2009). In this way, the transformed model ½P is the same as Meisel and Bierwirth (2009), which is shown in Appendix B. We now present the results of our GA and insertion heuristic algorithm for the problem ½P. The insertion heuristic algorithm is kept the same to solve the problem ½P, however GA requires small modifications to cope with ½P. The current decoding of our GA leads to that almost every vessel has a speedup cost, which results in a poor GA performance. Hence, the earliest starting time is replaced with the expected arrival time in the decoding of GA. This adjustment is only for the problem ½P, the presented GA mechanism is used for the remaining problems. Table 2 reports the results provided by the GA, the insertion heuristic algorithm, Meisel and Bierwirth (2009) and Iris et al. (2015a). Columns 3, 4, 6, 7, 9, 11, 13 demonstrate the total service cost, and Columns 5, 7, 10, 12, 14 demonstrate the computational gap against the best lower bound. Since the computational times of the insertion heuristic algorithm are all less than a second, computational times are reported (in seconds) for the GA only.

59

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

To compare the aggregate level of these heuristics, the observed computational gaps and computational times are averaged over every ten test instances. It can be seen that the GA and the insertion heuristic algorithm are better than the construction heuristic with local refinements procedures (CH), but are worse than the other two heuristics of Meisel and Bierwirth (2009), in term of solution quality. It is noted that almost all best bounds are obtained from Iris et al. (2015a). In term of computational time, the insertion heuristic algorithm is comparable with the CH, since their computational time is less than a second for every instance. Besides, the computational times of the GA and the insertion heuristic algorithm are all shorter than those of the other two heuristics of Meisel and Bierwirth (2009). The GA and the insertion heuristic algorithm are designed for three models in this paper, whereas the heuristics of Meisel and Bierwirth (2009) and Iris et al. (2015a) are designed to solve model ½P, but not applied to the models in this paper. 6.2.2. Computational results Tables 3–5 show the objective function values and computational times of 20, 30 and 40 vessels for models ½P 1 ½P 3 , respectively. For the small-scale instances (N ¼ 20), merely four instances (e.g., n ¼ 1; n ¼ 2; n ¼ 5; n ¼ 6) for model ½P 1 , three instances (e.g., n ¼ 2; n ¼ 5; n ¼ 6) for models ½P 2 ½P 3 are proven to be solved to optimality by CPLEX within the given runtime limit. Most of the medium-scale instances (N ¼ 30) remain unsolved using CPLEX, while not a single integer feasible solution has been found for the large scale instances (N ¼ 40). Obviously, the computational time of CPLEX rapidly increases, and the performance clearly turns bad with the increase of instance scale. However, the GA and the insertion heuristic algorithm always obtain feasible solutions for all instances. In term of solution quality, the GA is prior to the insertion heuristic algorithm, while in term of computational time, the

Table 3 Results for model ½P 1 . Example

CPLEX

GA

Heuristic

N

n

UB1

LB

UB1

Time (s)

Gap (%)

UB1

Time (s)

Gap (%)

20

1 2 3 4 5 6 7 8 9 10

109.8 88.0 110.8 101.1 89.0 99.0 103.0 90.4 101.2 106.6

109.8 88.0 99.7 96.1 89.0 99.0 101.2 88.6 95.1 100.2

110.0 88.0 106.2 102.6 89.0 99.0 103.4 89.6 97.2 104

7.84 7.82 9.98 8.91 8.92 8.44 8.52 7.51 7.76 8.05

0.2 0.0 6.5 6.8 0.0 0.0 2.2 1.1 2.0 3.8

113.4 88.4 111.4 104.2 89.0 99.0 105.6 92.4 103.8 111.8

0.16 0.14 0.14 0.16 0.15 0.13 0.15 0.17 0.13 0.14

3.3 0.5 11.7 8.4 0.0 0.0 4.4 4.3 9.2 11.6

8.38

2.3

0.15

5.3

21.58 15.07 26.23 23.37 27.93 25.24 22.44 24.95 21.75 28.16

10.2 4.2 2.8 2.7 7.9 6.2 2.6 5.8 8.5 2.9

0.25 0.24 0.31 0.26 0.23 0.24 0.26 0.24 0.23 0.23

15.3 6.6 6.7 10.7 11.5 8.6 7.3 6.1 16.2 10.6

23.67

5.4

0.25

10.0

56.36 82.72 73.16 69.00 61.17 70.37 65.15 61.34 53.91 53.44

20.9 14.7 34.7 30.9 15.8 21.7 22.0 30.3 12.0 12.9

0.29 0.32 0.36 0.38 0.34 0.43 0.35 0.37 0.40 0.34

32.5 20.9 41.8 41.1 21.2 27.1 28.4 29.5 20.8 28.2

64.66

21.6

0.36

29.2

Average 30

11 12 13 14 15 16 17 18 19 20

– 134.4 – 144.6 – – – – – –

146.8 129.7 144.4 137.9 154.1 139.4 143.9 141.4 167.2 151.0

161.8 135.2 148.4 141.6 166.2 148.0 147.6 149.6 181.4 155.4

Average 40

21 22 23 24 25 26 27 28 29 30

– – – – – – – – – –

184.4 185.0 177.1 193.3 176.1 194.8 194.2 204.6 182.5 199.0

223.0 212.2 238.6 252.6 204.0 237.0 237.0 266.6 204.4 224.6

Average UB1 /LB1 : the upper/lower bound of model ½P 1 within 10 h limit. Gap = ðUB1 LB1 Þ=LB1 . –: not a single integer feasible solution has been found within 10 h limit.

169.2 138.2 154.0 152.6 171.8 151.4 154.4 150.0 194.2 167.0

244.4 223.6 251.2 272.8 213.4 247.6 249.4 265.0 220.4 255.2

60

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

Table 4 Results for model ½P 2 . Example

CPLEX

GA

Heuristic

N

n

UB2

LB2

UB2

Time (s)

Gap (%)

UB2

Time (s)

Gap (%)

20

1 2 3 4 5 6 7 8 9 10

113.0 89.2 116.0 – 89.6 100.8 104.2 94.6 102.4 –

110.9 89.2 100.8 96.7 89.6 100.8 102.5 91.8 96.3 105.1

113.0 89.4 109.2 106.6 89.6 100.8 105.8 97.2 102.0 112.8

8.23 8.39 7.92 7.52 8.65 8.17 9.10 7.89 8.24 8.73

1.9 0.2 8.3 10.2 0.0 0.0 3.2 5.9 5.9 7.3

119.6 92.6 113.6 114.6 89.6 100.8 107.8 99.4 108.6 119.6

0.17 0.16 0.15 0.15 0.16 0.19 0.14 0.16 0.14 0.15

7.8 3.8 12.7 18.5 0.0 0.0 5.2 8.3 12.8 13.8

8.28

4.3

0.16

8.3

25.62 16.57 27.99 22.98 22.95 28.77 25.65 26.95 28.77 21.44

13.2 6.3 3.9 5.6 13.5 11.6 6.5 9.9 14.4 4.6

0.24 0.25 0.28 0.26 0.25 0.24 0.24 0.26 0.23 0.20

15.7 8.6 8.7 12.9 15.6 12.6 10.5 10.0 19.2 12.6

24.77

9.0

0.25

12.6

76.28 97.92 87.97 131.20 95.60 72.68 61.63 73.01 62.09 60.20

37.5 17.2 38.5 35.4 18.9 26.8 24.0 32.4 19.1 18.9

0.36 0.36 0.38 0.37 0.39 0.35 0.37 0.36 0.35 0.34

42.6 24.1 46.9 46.2 24.0 33.1 30.0 33.4 29.5 31.2

81.86

26.9

0.36

34.1

Average 30

11 12 13 14 15 16 17 18 19 20

183.4 136.0 – 146.8 – – – – – –

148.9 131.2 146.1 139.2 156.8 141.9 144.4 142.9 171.2 152.8

168.6 139.4 151.8 147.0 178.0 158.4 153.8 157.0 195.8 159.8

Average 40

21 22 23 24 25 26 27 28 29 30

– – – – – – – – – –

188.6 187.9 185.0 195.1 176.8 197.4 196.6 207.8 185.5 200.9

259.4 220.2 256.2 264.2 210.2 250.2 243.8 275.2 221.0 238.8

Average

172.2 142.4 158.8 157.2 181.2 159.8 159.6 157.2 204.0 172.0

269.0 233.2 271.8 285.2 219.2 262.8 255.6 277.2 240.2 263.6

UB2 /LB2 : the upper/lower bound of model ½P 2 within 10 h limit. Gap = ðUB2 LB2 Þ=LB2 . –: not a single integer feasible solution has been found within 10 h limit.

insertion heuristic algorithm is clearly more competitive than the GA. Comparing these three models, Constraints (30) in model ½P 3 have higher computational complexity than Constraints (3) in model ½P 1 and Constraints (25) in model ½P 2 . Hence, model ½P 3 needs more computational time than models ½P 2 ½P 3 for the GA and the insertion heuristic algorithm, as is shown in these three tables. Besides, almost all the objective values of robust models are larger than the deterministic model, which means that the uncertainties would increase the objective values and alter the planned schedule. That is to say, the schedule of deterministic model ½P 1 may become infeasible once uncertainties occur. Comparing Table 4 with Table 5, the solutions of model ½P 2 are much worse than the solutions of model ½P 3 . Due to the conservativeness of robust model ½P 2 , model ½P 3 is more feasible to deal with data uncertainties.

6.3. Impact of Ci on experimental results As mentioned before, the price constraint Ci takes value in the interval ½0; jT 0 j, and it is introduced to adjust the conservativeness of model ½P 3 . Fig. 10 illustrates the impact of Ci on average results of ten instances and the number of infeasible deterministic schedules by CPLEX, the GA, and the insertion heuristic algorithm. It is easy to find that average objective values and the number of infeasible schedules increases with Ci . When Ci ¼ 8, the number of infeasible schedules by CPLEX, the GA, and the insertion heuristic algorithm are all up to 10. That is to say, the schedule of the deterministic model becomes infeasible under uncertainties.

61

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65 Table 5 Results for model ½P 3 . Example

CPLEX

GA

Heuristic

N

n

UB3

LB3

UB3

Time (s)

Gap (%)

UB3

Time (s)

Gap (%)

20

1 2 3 4 5 6 7 8 9 10

111.0 89.2 110.4 – 89.6 100.8 103.6 92.6 – –

109.8 89.2 100.2 96.6 89.6 100.8 101.8 89.8 95.7 102.1

111.0 89.4 107.4 103.4 89.6 100.8 104.6 92.6 99.6 109.4

164.18 117.80 111.86 129.80 147.29 151.88 143.92 102.98 126.29 114.02

1.1 0.2 7.2 7.0 0.0 0.0 2.8 3.1 4.1 7.2

116.4 91.0 112.2 114.6 89.6 100.8 106.8 94.0 104.6 115.8

4.91 5.25 4.70 4.08 4.52 5.48 4.71 4.99 4.26 4.10

6.0 2.0 12.0 18.6 0.0 0.0 4.9 4.7 9.3 13.4

131.00

3.3

4.70

7.1

319.70 196.48 261.62 205.81 245.14 213.81 223.40 289.18 286.28 319.95

11.9 3.5 3.3 3.5 10.6 8.1 3.5 7.5 9.9 4.7

6.75 7.33 8.14 7.52 8.92 6.93 7.43 7.53 6.38 6.57

16.6 6.0 7.2 10.5 13.7 8.7 7.5 6.9 15.0 13.1

256.14

6.7

7.35

10.5

464.67 408.79 593.81 443.91 453.45 468.17 463.10 644.32 480.20 498.40

34.2 19.4 35.6 35.1 17.0 25.0 22.4 33.8 18.4 17.5

9.47 9.26 9.32 8.47 7.63 8.94 8.65 9.17 7.32 8.37

38.4 18.0 42.2 43.3 22.4 32.6 28.8 30.0 29.7 30.4

491.90

25.8

8.66

31.6

Average 30

11 12 13 14 15 16 17 18 19 20

188.6 143.6 – 146.4 – – 150.6 – – –

150.2 131.2 145.0 138.6 154.1 141.5 144.6 142.7 169.8 149.4

168.0 135.8 149.8 143.4 170.4 153.0 149.6 153.4 186.6 156.4

Average 40

21 22 23 24 25 26 27 28 29 30

– – – – – – – – – –

187.5 186.6 184.4 194.3 176.2 195.2 196.2 205.2 182.7 200.0

251.6 222.8 250.0 262.4 206.2 244.0 240.2 274.6 216.4 235.0

Average

175.2 139.0 155.4 153.2 175.2 153.8 155.4 152.6 195.2 169.0

259.4 220.2 262.2 278.4 215.6 258.8 252.6 266.8 237.0 260.8

UB3 /LB3 : the upper/lower bound of model ½P 3 within 10 h limit. Gap = ðUB3 LB3 Þ=LB3 . –: not a single integer feasible solution has been found within 10 h limit.

6.4. Probability bounds of constraint violation The price constraints Ci control the tradeoff between the probability of violation. In Section 3.6, we present the probability bounds of constraint violation determined by the price constraints Ci , and in this subsection, the concrete analysis for different Ci under model ½P 3 is given. Table 6 illustrates the choice of Ci so that the probability of constraint violation is less than 1% for different jT 0 j. For jT 0 j ¼ 168, we need to use Ci ¼ 39:4, i.e., only 23.5% of the number of uncertain data, to guarantee violation probability of less than 1%. For constraints with fewer number of uncertain data, it is necessary to ensure full protection. Fig. 11 illustrates the probability bounds of constraint violation for jT 0 j ¼ 168. At the beginning, the probability bound of constraint violation rapidly reduces with the increase of Ci , and then it is stable in a very small positive number closed to 0. 7. Conclusions Since the BAP and the QCAP have close relationship, this paper studies the two problems in an integrated way. Considering the QC setup time of shifting along the quay, a deterministic model is proposed for the BACAP. However, as a result of the special operational environment of QC, the feasibility of schedule is affected by the uncertainties of the environment. To assure the validity of the planned schedule under data uncertainties, a robust optimization model based on the deterministic model is proposed to handle data uncertainties. Based on the robust optimization model, a new robust optimization model with the price constraints is formulated to control the level of conservativeness of the robust

62

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

Fig. 10. Impact of Ci on experimental results.

Table 6 Choice of Ci so that the probability of constraint violation is less than 1%. jT 0 j

5

10

15

24

48

72

96

120

144

168

Choice of Ci

5

9.6

11.8

14.9

21.1

25.8

29.8

33.3

36.5

39.4

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

63

Fig. 11. Probability bound of constraint violation for jT 0 j ¼ 168.

optimization. Due to the NP-hardness of the BACAP, we present a GA and an insertion heuristic algorithm to obtain near optimal solutions. Computational experiments show that the proposed models are reasonable and important, and the presented algorithms are effective. In term of computational time, the insertion heuristic algorithm is much shorter than the GA, while in term of solution quality, the GA is prior to the insertion heuristic algorithm. Further studies include the BACAP with the consideration of other uncertain factors and the development of exact solution algorithms to the BACAP under data uncertainties. Acknowledgements This work is supported by the National Nature Science Foundation of China (Grant No. 71262008), the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, and Young Talent on Science and Technology Program of University in Inner Mongolia (Grant No. NJYT-13-B02). Appendix A Proof. Let

n

o ritq j8t 2 T; q 2 Ri ; i 2 V ; Si and

vi be the solution of model ½P 3 . Firstly, the proof of the left-hand side of

Inequality (40) is as follows:

P

XX ~ r itq < ð1 þ b Dbi Þ mi q

!

t2T q2Ri

¼P

XX XX ^ r itq < ð1 þ b Dbi Þ mi qa r itq nit q 0

6 [email protected] 0 ¼ [email protected]

!

t2T 0 q2Ri

t2T q2Ri

1 XX X XX ^ r itq þ ^ r iv q < ^ ritq A q ðCi bCi cÞq nit q t2Si q2Ri

i

q2Ri

1 XX X XX ^ ritq ð1 nit Þ þ ^ r iv q < ^ r itq A q ðCi bCi cÞq nit q

0

t2Si q2Ri

i

q2Ri

t2T 0 nSi q2Ri

1 XX X XX ^ r ik q ð1 nit Þ þ ^ r ik q < ^ r itq A q 6 [email protected] ðCi bCi cÞq nit q 0

t2Si q2Ri

t2T 0 nSi q2Ri

q2Ri

1 0 1 X XX X ^ r ik q @ ð1 nit Þ þ ðCi bCi cÞA < ^ ritq A 6 [email protected] q nit q 0

ð41Þ

t2T 0 q2Ri

t2Si

q2Ri

1

X ^ r itq q P q2Ri nit þ nit A ^ q2Ri q r ik q t2Si t2T 0 nSi 0 1 X nit cit A: ¼ [email protected] < ¼ [email protected] <

X

P

t2T 0 nSi

t2T 0 nSi q2Ri

ð42Þ

64

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

n o The right-side of Inequality (40) provides a bound that is independent of the solution r itq j8t 2 T; q 2 Ri ; i 2 V and it is proven as follows: Let h > 0. Then

P

X nit cit P Ci t2T 0

¼ ¼

Q

P E exp h t2T 0 nit cit exp ðh Ci Þ

ð43Þ

nit cit Þ exp ðh Ci Þ i R 1 P1 h Q 2k k¼0 ðh cit nÞ =ð2kÞ! dF nit ðnÞ t2T 0 2 0 t2T 0

P1 h k¼0

ð44Þ ð45Þ

exp ðh Ci Þ i ðh cit Þ2k =ð2kÞ!

exp ðh Ci Þ 2 2 exp h cit =2 t2T exp ðh Ci Þ

Q 6

6

t2T 0 E½exp ðh

Q 6

!

0

6 exp jT 0 j

! h2 h Ci : 2

ð46Þ

Remark. 1. Inequality (41) follows from inequality (29), since r itq satisfies

8 9

i

i

8i 2 V:

i

^ r ikq . 2. Inequality (42) follows 1 nit P 0 and k ¼ arg mink2S Sfv g q i

i

3. Inequality (43) follows from Markov’s Inequality. Equalities (44) and (45) follow from the independence symmetric distribution assumption of the random variables nit . Inequality (46) follows from cit 1. Selecting h ¼ Ci =jT 0 j, we can obtain the right-side of Inequality (40). Appendix B Input data c1i ; c2i ; c3i c4 EFTi LFTi DETAi

service cost rate for vessel i given in units of 1000 USD per hour operation cost rate given in units of 1000 USD per QC-hour expected finishing time of vessel i latest finishing time of vessel i without penalty cost arising integer, required speed up of vessel i to reach its berthing time, DETAi ¼ ðETAi si Þþ integer, tardiness of vessel i; D EFTi ¼ ðei EFT i Þþ

DEFTi Decision variables u0i

½P

min

binary, set to 1 if the finishing time of vessel i exceed LFTi, 0 otherwise

Z¼

X

c1i

DETAi þ

c2i

DEFT i þ

c3i

u0i

XX þc q r itq

i2V

s:t:

ð3Þ; ð6Þ; ð12Þ—ð20Þ; ð24Þ XX q r itq 6 Q 8t 2 T i2V q2Ri

DETAi P ETAi si 8i 2 V DEFTi P ei EFTi 8i 2 V M u0i P ei LFTi 8i 2 V si 2 fEST i ; . . . ; H 1g; ei 2 fEST i ; . . . ; Hg:

4

t2T q2Ri

!

X.T. Shang et al. / Transportation Research Part E 94 (2016) 44–65

65

References Bertsimas, D., Sim, M., 2004. The price of robustness. Oper. Res. 52, 35–53. Bierwirth, C., Meisel, F., 2015. A follow-up survey of berth allocation and quay crane scheduling problems in container terminals. Eur. J. Oper. Res. 244, 675– 689. Chen, J.H., Lee, D.H., Cao, J.X., 2011. Heuristic for quay crane scheduling at indented berth. Transport. Res. Part E 47, 1005–1020. Cordeau, J.F., Laporte, G., Legato, P., Moccia, L., 2005. Models and tabu search heuristics for the berth-allocation problem. Transport. Sci. 39 (40), 526–538. Dagnazo, C.F., 1989. The cranes scheduling problem. Transport. Res. Part B 23 (3), 159–175. Gen, M., Cheng, R., 1996. Genetic Algorithms and Engineering Design. John Wiley, Singapore. Golias, M., Portal, I., Konur, D., Kaisar, E., Kolomvos, G., 2014. Robust berth scheduling at marine container terminals via hierarchical optimization. Comput. Oper. Res. 41, 412–422. Imai, A., Nagaiwa, K., Tat, C.W., 1997. Efficient planning of berth allocation for container terminals in Asia. J. Adv. Transport. 31 (1), 75–94. Imai, A., Nishimura, E., Papadimitriou, S., 2001. The dynamic berth allocation problem for a container port. Transport. Res. Part B 35 (4), 401–417. Imai, A., Nishimura, E., Hattori, M., Papadimitriou, S., 2007. Berth allocation at indented berths for mega-containerships. Eur. J. Oper. Res. 179 (2), 579–593. Imai, A., Nishimura, E., Papadimitriou, S., 2008. The simultaneous berth and quay crane allocation problem. Transport. Res. Part E 44 (5), 900–920. Iris, C., Pacino, D., Popke, S., Larsen, A., 2015a. Integrated berth allocation and quay crane assignment problem: set partitioning models and computational results. Transport. Res. Part E 81, 75–97. Iris, C., Jin, J.G., Lee, D.-H., 2015b. Quayside operations planning under uncertainty. In: European Conference for Operational Research 2015, Glasgow, United Kingdom. Karafa, J., Golias, M.M., Ivey, S., Saharidis, G.K.D., Leonardos, N., 2013. The berth allocation problem with stochastic vessel handling times. Int. J. Adv. Manuf. Technol. 65 (1–4), 473–484. Kim, K.H., Moon, K.C., 2003. Berth scheduling by simulated annealing. Transport. Res. Part B 57, 541–560. Kim, K.H., Park, Y.-M., 2004. A crane scheduling method for port container terminals. Eur. J. Oper. Res. 156, 752–768. Lai, K.K., Shih, K., 1992. A study of container berth allocation. J. Adv. Transport. 26, 45–60. Lalla-Ruiz, E., Gonzalez-Velarde, J., Meli¢n-Batistaa, B., Moreno-Vega, J., 2014. Biased random key genetic algorithm for the tactical berth allocation problem. Appl. Soft Comput. 22, 60–76. Lee, D.H., Chen, J.H., Cao, J.X., 2010. The continuous berth allocation problem: a greedy randomized adaptive search solution. Transport. Res. Part E 46, 1017– 1029. Lee, D.H., Wang, H.Q., 2010. Integrated discrete berth allocation and quay crane scheduling in port container terminals. Eng. Optimiz. 42 (8), 747–761. Li, F., Sheu, J.-B., Gao, Z.-Y., 2015. Solving the continuous berth allocation and specific quay crane assignment problems with quay crane coverage range. Transport. Sci. 49 (4), 968–989. Meisel, F., Bierwirth, C., 2009. Heuristic for the integration of crane productivity in the berth allocation problem. Transport. Res. Part E 45, 196–209. Park, Y.-M., Kim, K.-H., 2003. A scheduling method for berth and quay cranes. Oper. Res. Spectrum 25, 1–23. Rodriguez-Molins, M., Barber, F., Sierre, M., Puente, J., Salido, M., 2012. A genetic algorithm for berth allocation and quay crane assignment. In: Lecture notes in computer science, vol. 7637, pp. 601–610. Rodriguez-Molins, M., Ingolotti, L., Barber, F., Salido, M., Sierra, M., Puente, J., 2014. A genetic algorithm for robust berth allocation and quay crane assignment. Prog. Artif. Intell. 2 (4), 177–192. Schonfeld, P., Sharafeldien, O., 1985. Optimal berth and quay crane combinations in container ports. J. Waterw. Port Coast. Ocean Eng. 111 (6), 1060–1072. Stahlbock, R., Voß, S., 2008. Operations research at container terminals: a literature update. OR Spectrum 30, 1–52. Steenken, D., Voß, S., Stahlbock, R., 2004. Container terminal operation and operations research – a classification and literature review. OR Spectrum 26, 3– 49. Turkogullari, Y.B., Taskin, Z.C., Aras, N., Altinel, I.K., 2016. Optimal berth allocation, time-variant quay crane assignment and scheduling with crane setups in container terminals. Eur. J. Oper. Res. 254 (3), 985–1001. Vis, I.F.A., Koster, R., 2003. Transshipment of containers at container terminal: an overview. Eur. J. Oper. Res. 147, 1–16. Zhen, L., 2015. Tactical berth allocation under uncertainty. Eur. J. Oper. Res. 247, 928–944. Zhou, P.-F., Kang, H.-G., 2008. Study on berth and quay-crane allocation under stochastic environments in container terminal. Syst. Eng.-Theory Pract. 28 (1), 161–169.