A decision model for berth allocation under uncertainty

European Journal of Operational Research 212 (2011) 54–68

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Production, Manufacturing and Logistics

A decision model for berth allocation under uncertainty Lu Zhen a,⇑, Loo Hay Lee b, Ek Peng Chew b a b

School of Management, Shanghai University, Shanghai 200444, PR China Department of Industrial & System Engineering, National University of Singapore, Singapore 117576, Singapore

a r t i c l e

i n f o

Article history: Received 24 December 2009 Accepted 15 January 2011 Available online 21 January 2011 Keywords: Scheduling Port operation Berth allocation Container terminals Meta-heuristic

a b s t r a c t This paper studies the berth allocation problem (BAP) under uncertain arrival time or operation time of vessels. It does not only concern the proactive strategy to develop an initial schedule that incorporates a degree of anticipation of uncertainty during the schedule’s execution, but also studies the reactive recovery strategy which adjusts the initial schedule to handle realistic scenarios with minimum penalty cost of deviating from the initial schedule. A two-stage decision model is developed for the BAP under uncertainties. Moreover, a meta-heuristic approach is proposed for solving the above problem in large-scale realistic environments. Numerical experiments are conducted to validate the effectiveness and efficiency of the proposed method. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction The berth allocation problem (BAP) is an important issue for the operation management in container terminals (Murty et al., 2005). It is mainly referred to the assignment of quay space and service time to vessels that have to be unloaded and loaded at a terminal. For a comprehensive overview on BAP, we refer to the review work given by Bierwirth and Meisel (2010). Most of researches on BAP concern how to obtain an initial schedule (baseline schedule, or berth allocation template) in a static and deterministic environment with complete information. In realistic situations, however, these assumptions will hardly be satisfied, there are a lot of uncertain factors and unexpected events, e.g., deviation of arrival time, operation time, inserting another vessel into the already arranged berthing plan. These unexpected events will result in last minute change of plans. Therefore, port planners have to consider the adverse effects of possible disruptions in their planning of the initial schedule. This paper studies how to obtain an initial schedule (baseline schedule) for berth allocation under uncertain arrival time or operation time of vessels. Here the uncertainty is represented by a finite set of discrete scenarios. There are usually two strategies for coping with uncertainties: proactive and reactive scheduling. This paper covers both two strategies for berth allocation scheduling. It does not only concern the proactive strategy that incorporates a degree of anticipation of uncertainty and variability during the schedule’s execution, but also studies the reactive recovery strategy which adjusts the baseline schedule to handle realistic scenarios with minimum penalty cost of deviating from the initial schedule. A two-stage decision model is developed to balance the cost associated with the initial schedule and the expected cost of deviation from the initial schedule. Moreover, a meta-heuristic approach is proposed for solving the above problem in large-scale realistic environments. Numerical experiments are conducted to validate the effectiveness and efficiency of the proposed method. 2. Literature review There are numerous studies on the berth allocation problem (BAP) in the past two decades. Imai et al. addressed the static BAP (SBAP) in commercial ports in Asia (Imai et al., 1997), and then they extended the SBAP to DBAP (dynamic BAP) (Imai et al., 2001). It should be mentioned that the stochastic BAP studied in this paper is totally different from the DBAP. The DBAP’s input is essentially deterministic, yet this study concerns the strategies for scheduling berth allocation under uncertainties. BAP could be classified into two types: discrete and continuous problems (Imai et al., 2005). As to the discrete problems, the quay is partitioned into a number of sections (berths), where ⇑ Corresponding author. Tel.: +86 13816890752. E-mail addresses: [email protected] (L. Zhen), [email protected] (L.H. Lee), [email protected] (E.P. Chew). 0377-2217/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2011.01.021

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

55

one ship could be handled at a time. Whereas in the continuous problems, ships could be served wherever empty spaces are available, which resembles a two-dimensional un-rotated bin packing problem (Lim, 1998). As to the solution methodology for the continuous BAP, Kim and Moon proposed a simulated annealing method (Kim and Moon, 2003); Park and Kim employed a sub-gradient optimization method (Park and Kim, 2002). Guan et al. developed a heuristic with minimizing total weighted completion time (Guan and Cheung, 2004). Wang and Lim developed a novel beam search method for BAP (Wang and Lim, 2007). Nishimura et al. proposed a genetic algorithm for BAP to obtain a good solution with small computational effort (Nishimura et al., 2001). Imai et al. investigated BAP for the indented berths, where mega-containerships could be served from two sides (Imai et al., 2007). Cordeau et al. studied BAP in some specific quay which consists of variable berths (Cordeau et al., 2005). Hansen et al. proposed a variable neighbor search method for BAP, in which vessels’ handing time and cost depend on berthing positions (Hansen et al., 2008). The integration of BAP and quay crane (QC) assignment has been widely studied in recent years. Park and Kim developed a two-phase solution procedure, and a dynamic programming technique was applied to solve the problem (Park and Kim, 2003). Kim and Park addressed a crane assignment problem for minimizing the turnaround time of containerships (Kim and Park, 2004). Meisel and Bierwirth treated the BAP-QC assignment as a multi-mode resource constrained project scheduling problem. The squeaky wheel optimization (SWO) and tabu search meta-heuristics were employed to solve the problem (Meisel and Bierwirth, 2009). Imai et al. studied the simultaneous berth and QC allocation problem, which took into account QCs could not move freely along berth track. When one vessel is handling, QCs could not pass or bypass from its one side to the other (Imai et al., 2008). In addition, discrete event simulation methods were utilized to solve berth planning problems (Legato and Mazza, 2001). For berth template design, Moorthy and Teo developed a model from two aspects: the percentage of vessels berthed on arrival and the connectivity cost (Moorthy and Teo, 2006). A novel sequence pair based approach was proposed. In all, BAP is a well researched domain. There are plenty of articles in this area. However, to our best knowledge, few scholars have specially studied the berth allocation scheduling under uncertainties. This study could be regarded as a supplement to the ordinary deterministic BAP.

3. Problem background For the traditional deterministic BAP, port planners face two interrelated decisions: where and when the vessels should be moored. The continuous BAP is represented in two-dimensional space, shown in Fig. 1. The objective function is usually to minimize the waiting time (or tardiness time) cost in the time axis and position cost in berth axis. Here the position cost is used to minimize deviations between the chosen berthing positions and the desired positions, which aim to reduce the travel distances for the transport vehicles in loading and unloading activities. As to BAP, there are a lot of different mathematical models. In this study, we refer to the model proposed by Kim and Moon (2003). In their model, the objective function is to minimize tardiness cost and position cost, c1i(xi + ti  di)+ + c2ijyi  bij, here c1i, c2i are cost coefficients, xi, yi are berthing time and position for vessel i, ti is the operation time, di is the requested departure time, and bi is the berthing position with the lowest cost. In this study, we minimize the waiting cost instead of the tardiness cost. Thus the cost function for vessel i is expressed as c1i(xi  ai) + c2ijyi  bij. There are plenty of uncertain factors in BAP. In this study, we mainly consider two uncertain factors: deviation of vessels’ arrival time and operation time. The first factor means that the vessels’ actual arrival time deviates from their estimated arrival time. Because of delay in previous ports or unforeseen events in voyage courses, vessels may arrive later than the estimated time. The second one means that the vessels’ actual operation time deviates from their estimated operation time. The operation time is usually estimated according to the number of loading and unloading containers. Sometimes, the actual number deviates from the one that is reported from the vessels to the port due to some unforeseen changes. In addition, the stack position of containers in vessels may cause unproductive reshuffles. So vessels’ actual operation time may deviate from their estimated one.

4. Model formulation V L ai

the set of all vessels, V = {V1, V2, . . . , VN}, N is the number of vessels the length of the wharf the estimated arrival time of vessel i

Fig. 1. Berth-time space for BAP.

56

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

ti bi li c1i

the the the the

c2i

the cost rate, resulting from deviation from its lowest cost position

þ c1i ; c1 i

the cost rate (or reward rate), resulting from more (or less) delay in berthing in recovery process, c1 < 0; here i

þ c2i ; c2 i

X p(xs) S ai(xs) ti(xs) U(xs)

s(xs)

estimated operation time of vessel i berthing position with the lowest cost for vessel i length of vessel i cost rate, resulting from delay in berthing beyond its arrival time

þ

(it will be explained later) c1i P c1i P c1 i the cost rate (or reward rate), resulting from more (or less) deviation from its lowest cost position in recovery process, c2 < 0; here c2þ P c2i P c2 i i i {x1, . . . , xs, . . . , xS}, the set of discrete future scenarios the probability of scenario xs the number of scenarios the actual arrival time of vessel i in scenario xs the actual operation time of vessel i in scenario xs the set of vessels that have varied arrival time or operation time in scenario xs. If ai(xs) – ai, or ti(xs) – ti, then i 2 U(xs); for i 2 Uðxs Þ; ai ðxs Þ ¼ ai , and ti(xs) = ti the frozen period point in scenario xs ; sðxs Þ ¼ minj2Uðxs Þ aj . If a vessel arrives before this point, its schedule is frozen; otherwise the schedule may be adjusted

Fig. 2 illustrates an example to explain the above defined U(xs) and s(xs). For scenario xs, as a3(xs) – a3 and a4(xs) – a4, U(xs) = {Vessel 3, 4}. Thus sðxs Þ ¼ minj2Uðxs Þ aj ¼ a3 , which means: as to the vessels that arrive before this point, i.e., Vessel 1 and 2, their schedules are frozen; whereas the schedules of Vessel 3, 4 and 5 may be adjusted in the recovery plan for scenario xs. The reason for setting the frozen period lies in: the original schedule will be executed until the frozen period point when the first unexpected event is observed, i.e., Vessel 3 is late in the example of Fig. 2. Thus the plan from that time point will be rescheduled. Port planners could generate the scenarios and probabilities based on the historical data of vessels’ past berthing records, which include the actual and estimated arrival times (or operation times) when they moored at the port. In practice, p(xs) could be derived from the probabilities of unpunctuality in the past records; ai(xs) (or ti(xs)) could also be obtained by adding the estimated time and the average deviation of actual arrival time (or operation time) in the past records. Decision variables: xi the yi the  the xþ ð x Þ; x ð x Þ s s i i the yþ ðxs Þ; y i ðxs Þ i xDþ ðxs Þ; xDi  ðxs Þ the i

yDi þ ðxs Þ; yDi  ðxs Þ dxij

start berthing time of vessel i in the baseline schedule berthing position of vessel i in the baseline schedule increment (or decrement) with respect to xi in scenario xs increment (or decrement) with respect to yi in scenario xs increment (or decrement) with respect to ‘xi  ai’ in scenario xs

the increment (or decrement) with respect to jyi  bij in scenario xs dxij ¼ 1, if vessel i is located in the left of vessel j in the 2-dimensional time-berth plane of baseline schedule (as shown in Fig. 1); dxij ¼ 0, otherwise

dyij

dyij ¼ 1, if vessel i is located below vessel j in the 2-dimensional time-berth plane of baseline schedule; dyij ¼ 0,

dxij ðxs Þ

otherwise dxij ðxs Þ ¼ 1, if vessel i is located in the left of vessel j in the 2-dimensional time-berth plane of the scenario xs’s schedule; dxij ðxs Þ ¼ 0, otherwise

Fig. 2. An example for illustrating the frozen period in scenario xs.

57

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

dyij ðxs Þ

dyij ðxs Þ ¼ 1, if vessel i is located below vessel j in the 2-dimensional time-berth plane of the scenario xs’s schedule; dyij ðxs Þ ¼ 0, otherwise

Mathematical model for BAP under uncertainties: N X 

ðM 0Þ Minimize

( ) S N X  X  1þ Dþ  c1i  ðxi  ai Þ þ c2i  jyi  bi j þ pðxs Þ  ci  xi ðxs Þ þ c1  xiD ðxs Þ þ c2þ  yiDþ ðxs Þ þ c2  yiD ðxs Þ i i i

i¼1



s¼1



8i; j 2 V; i – j; Subject to : xi þ ti 6 xj þ M 1  dxij   y yi þ li 6 yj þ M 1  dij 8i; j 2 V; i – j; dxij

þ

dyij

þ

dxji

þ

dyji

ð4:2Þ ð4:3Þ

P 1 8i; j 2 V; i – j;

ð4:4Þ

yi þ li 6 L 8i 2 V; xi P a i

ð4:5Þ

8i 2 V;

ð4:6Þ

xi þ xþi ðxs Þ  xi ðxs Þ þ t i ðxs Þ 6 xj þ xþj ðxs Þ  xj ðxs Þ þ Mð1  dxij ðxs ÞÞ 8i; j 2 V; i – j; 8xs 2 X; yi þ

yþi ð

xs Þ 

yi ð

ð4:1Þ

i¼1

xs Þ þ li 6 yj þ

yþj ð

xs Þ 

yj ð

xs Þ þ Mð1 

dyij ð

xs ÞÞ 8i; j 2 V; i – j; 8xs 2 X;

dxij ðxs Þ þ dyij ðxs Þ þ dxji ðxs Þ þ dyji ðxs Þ P 1 8i; j 2 V; i – j; 8xs 2 X;

ð4:7Þ ð4:8Þ ð4:9Þ

li 6 yi þ yþi ðxs Þ  yi ðxs Þ þ li 6 L 8i 2 V; 8xs 2 X;

ð4:10Þ

xi þ xþi ðxs Þ  xi ðxs Þ P ai ðxs Þ 8i 2 V; 8xs 2 X;

ð4:11Þ

xi  ai þ xiDþ ðxs Þ  xDi  ðxs Þ ¼ xi þ xþi ðxs Þ  xi ðxs Þ  ai ðxs Þ 8i 2 V; 8xs 2 X;

ð4:12Þ

jyi  bj þ yDi þ ðxs Þ  yDi  ðxs Þ ¼ jyi þ yþi ðxs Þ  yi ðxs Þ  bi j 8i 2 V; 8xs 2 X;

ð4:13Þ

xþi ðxs Þ; xi ðxs Þ; yþi ðxs Þ; yi ðxs Þ 6 Mðai ðxs Þ  sðxs ÞÞþ

ð4:14Þ

8i 2 V; 8xs 2 X:

Decision variables : xi ; yi ; xþi ðxs Þ; xi ðxs Þ; yþi ðxs Þ; yi ðxs Þ; xDi þ ðxs Þ; xiD ðxs Þ; yiDþ ðxs Þ; yiD ðxs Þ P 0 8i 2 V; 8xs 2 X; dxij ; dyij ; dxij ; ðxs Þ; dyij ; ðxs Þ 2 f0; 1g 8i; j 2 V; i – j:

ð4:15Þ ð4:16Þ

Due to (4.1) and (4.13) contain absolute value forms, we linearize them and rewrite M0 as:

Minimize

N X  1   ci  ðxi  ai Þ þ c2i  hþi þ hi i¼1

þ

S X

(

N n o X þ D 2þ Dþ 2 D pðxs Þ  c1i  xiDþ ðxs Þ þ c1 i  xi ðxs Þ þ c i  yi ðxs Þ þ c i  yi ðxs Þ

s¼1

) ð4:17Þ

i¼1

Subject to : Constraints ð4:2Þ—ð4:12Þ andð4:14Þ—ð4:16Þ; hþi þ hi þ yiDþ ðxs Þ  yiD ðxs Þ ¼ nþi þ ni ;

ð4:18Þ

yi  bi ¼ hþi  hi ;

ð4:19Þ

yi þ yþi ð s Þ  yi ð s Þ hþi ; hi ; nþi ; ni P 0:

x

x  bi ¼

nþi



ni ;

ð4:20Þ ð4:21Þ

The objective is to minimize the cost of baseline schedule and the expected value of the recovery costs. The baseline schedule, also named by tactical berth allocation plan or home berth template, is the basis for making other tactical/operational plans in port operations, e.g., yard template, storage allocation plans, and some resources plans in operational level (Moorthy and Teo, 2006). Thus, if the realistic schedule deviates from the baseline schedule, this adjustment will incur extra cost that is named by ‘recovery cost’ in this study. Thus the objective is to minimize the baseline schedule cost and recovery cost simultaneously.  1þ Dþ  1þ As mentioned in the notation, c1þ P c1i P c1 and c2þ P c2i P c2  xi ðxs Þ i i . It means: the coefficients ‘c i ’ of the ‘penalty’ cost c i i i   and ‘c1 the ‘reward’ c1  xiD ðxs Þ in recovery process are different from the counterpart ‘c1i ’ in baseline schedule’s cost i ’ of i  1  ci  ðxi  ai Þ . This difference is caused by the extra cost of adjusting other tactical/operational plans in the terminal. For example, we assume there are two different plans. In Plan_a, the waiting time for vessel i is 2 hours in the baseline schedule, then turns to 5 hours due to adjustment in realistic scenario (i.e., xi  ai ¼ 2; xDi þ ðxs Þ ¼ 3; xiD ðxs Þ ¼ 0). In Plan_b, the waiting time for vessel i is 5 hours in the baseline schedule, and has no adjustment in realistic scenario (i.e., xi  ai ¼ 5; xiDþ ðxs Þ ¼ 0; xDi  ðxs Þ ¼ 0). If the cost coefficients of the baseline and recovery adjustments are the same, i.e., c1þ ¼ c1i ¼ c1 i , the cost of Plan_a and Plan_b would have no difference. It disobeys the common i sense (or the usual practice) that the Plan_b should be better than Plan_a, because Plan_a has adjustment in recovery process that will cause a lot of changes to other resources plan, e.g., prime movers, quay cranes, human resources, etc. Hence, Plan_a should incur a higher cost than Plan_b in realistic environments. Therefore, the coefficients of the above ‘penalty’ cost and ‘reward’ should be different from the counterparts in baseline schedule’s cost. More specifically, the coefficient c1þ should be larger than c1i ; and the absolute value of the coefficient  2 i Dþ    1 1 ci should be smaller than ci . The analysis on the ‘penalty’ cost ci  yi ðxs Þ and ‘reward’ c2i  yiD ðxs Þ is also similar. Constraints (4.2)–(4.4) ensure there is no overlap among ‘rectangles’ in the 2-dimensional time-berth plane; here these rectangles represent the schedules for vessels, as shown in Fig. 1. Constraint (4.5) implies the positions of vessels are restricted by the length of the wharf. Constraint (4.6) enforces the berthing time cannot be earlier than the arrival time. Constraints (4.7)–(4.9) ensure there is no overlap among ‘rectangles’ in the time-berth plane after adjusting vessels’ schedules. Constraint (4.10) implies the newly planned positions of vessels after adjusting vessels’ schedules are restricted by the length of the wharf. Constraint (4.11) ensures the newly planned berthing time should not    be earlier than the actual arrival time. (4.12) builds the relationship between the adjustment of berthing time xþ i ðxs Þ; xi ðxs Þ  Dþ Constraint  D and the change of waiting time xi ðxs Þ; xi ðxs Þ . Constraint (4.13) is the relationship between the adjustment of berthing position

58

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

 þ   Dþ  D yi ðxs Þ; y i ðxs Þ and the change of deviation from the lowest-cost position yi ðxs Þ; yi ðxs Þ . Constraint (4.14) ensures that: if the vessels’ arrival time ai(xs) 6 s(xs), their schedules have no change, and will be executed according to the initial baseline schedule. Here, s(xs) is the frozen period point in scenario xs. If a vessel’s arrival time is before this point, its schedule is frozen; otherwise it may be adjusted. 5. Meta-heuristic for large-scale problems The large-scale deterministic BAP is hard to solve by commercial software, not to mention M0, especially with a lot of scenarios. Thus a meta-heuristic is suggested to solve large-scale M0. 5.1. General framework   þ   The solution for M0 is to obtain a baseline schedule (xi, yi) and recovery schedules xþ i ðxs Þ; xi ðxs Þ; yi ðxs Þ; yi ðxs Þ . Since solving the problem is not easy, we will solve it in stages. Firstly, we try to find a baseline schedule. Then based on the baseline schedule, we find the recovery schedules for different scenarios in order to evaluate the performance of the schedule. For obtaining a baseline schedule in large-scale problems, it is time-consuming to solve a deterministic BAP model with P  ‘ Ni¼1 c1i  ðxi  ai Þ þ c2i  jyi  bi j ’ as the objective function and (4.2)–(4.6) as constraints. In this study, a sequence-based method is used to obtain a baseline schedule. It means: when a sequence (or priority sequence) of vessels is given, we insert these vessels into the 2dimensional berth-time plane sequentially, and determine a baseline schedule (denoted by BS). Then based on BS, we determine a series of recovery schedules RSs for all scenarios xs 2 X. The BS and RSs are obtained in two stages, so it may cause some loss of optimality of the original problem ðM0Þ which optimizes BS and RSs simultaneously. Therefore, given the BS and RSs, we enhance the solution by solving M0   with some integer variables fixed, and obtain BS⁄ and RSs . Then we calculate the objective value of the solution BS ; RSs according to (4.1). From the above, when a sequence of vessels is given, we can obtain a solution (including a baseline schedule and a series of recovery schedules) and its objective value. We can improve the solution by consider other sequences of vessels. In this study, the proposed metaheuristic approach is built on a framework of simulated annealing (SA), which is a generic probabilistic meta-heuristic for the global optimization problem. The typical simulated annealing approach involves a pair of nested loops and some additional parameters, e.g., cooling rate, 0 < r < 1, and temperature length, R, which represents the size of neighborhood of a solution. The following describes the procedure: Step 1: Obtain an initial solution, let T = T0 (the initial temperature in SA). Step 1.1: Generate an initial sequence seq. Step 1.2: Based on seq, obtain a baseline schedule BS. Step 1.3: Based on BS, obtain a series of recovery schedules RSs for all scenarios xs 2 X.    Step 1.4: Enhance the solution (BS, RSs), and obtain a new BS ; RSs .  solution   Step 1.5: Calculate the objective value of the solution BS ; RSs , denoted by F(seq). Step 2: Repeat the following steps until one of stopping conditions becomes true. Step 2.1: Generate R neighbor sequences of seq, i.e., seqn, n 2 {1, 2, . . . , R}. Step 2.2: For n = 1 to R, repeat the following steps: Step 2.2.1: Based on seqn, obtain a baseline schedule BS. Step 2.2.2: Based on BS, obtain a series of RSs for all scenarios xs 2 X.  Step 2.2.3: Enhance the solution (BS, RSs), and obtain a new solution BS ; RSs .    Step 2.2.4: Calculate the objective value of BS ; RSs , denoted by F(seqn). Step 2.2.5: Let d = F(seqn)  F(seq). Step 2.2.6: If d < 0, set seq = seqn. Step 2.2.7: If d P 0, generate a random number x 2 (0, 1); If x < ed/T,seq = seqn. Step 2.3: Set T = r  T. The stopping criterions: (a) the fitness value of a solution is zero; (b) temperature becomes less than a given threshold value; (c) the best value of the fitness has not been improved during a given number of external loops. There are some key issues for implementing the above SA procedure: (1) (2) (3) (4)

How How How How

to to to to

obtain a baseline schedule with a sequence of vessels given. obtain recovery schedules for scenarios with a baseline schedule given. enhance a solution which includes a baseline schedule and recovery schedules. define the neighborhood set of a sequence.

The following Sections 5.2–5.5 will address the above four issues respectively in detail. 5.2. Obtain an baseline schedule for a given sequence of vessels This section introduces the procedure on obtaining a good baseline schedule if an inserting sequence of vessels is given. This procedure tries to insert all the N vessels into the 2-dimensional time-berth plane one by one according to the sequence. The initial sequence could be set according to the increasing order of vessels’ estimated arrival time. Based on this sequence, the procedure could obtain an initial baseline schedule. This procedure consists of N iterations for inserting the N vessels. Let the inserting sequence of vessels be {V1, . . . , Vn, . . . , VN}. Each iteration contains three steps: (1) inserting one vessel, (2) greedy selection, and (3) look-ahead selection.

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

59

Fig. 3. Insert a vessel in its searching space.

In the first iteration, vessel V1 is inserted into the plane within its searching space. The searching space for vessel V1 is defined as: {(x, y)ja1 6 x 6 a1 + s  t1; 0 6 y 6 L  l1}. Here s is a pre-defined number. Fig. 3 shows an example of the searching space for vessel V3; the searching space for the first vessel V1 is formulated similarly. The searching space is a discretized area. The lattice size is pre-defined according to the scale of the problem and the computer’s capacity. When inserting V1, we enumerate all the feasible points (x1, y1) in its searching space, and select the top NG best points (x1, y1) according to the cost of baseline schedule, i.e., c11  ðx1  a1 Þ þ c21  jy1  b1 j. Here, NG is a pre-defined parameter, e.g., NG = 100 in this study. Then, among these NG points, we select the top NS best points (x1, y1) in the stage of ‘look-ahead selection’ by considering some more vessels inserted into the plane. Here, NS is also a pre-defined parameter, e.g., NS = 30 in this study. For how to evaluate the goodness of the above NG points so as to select the top NS best ones, it will be explained later (model (5.1)–(5.8) with n = 1). The set of these NS obtained points (x1, y1), denoted by S1, will be the basis for the second iteration. The nth iteration is similar to the first iteration, except that in the nth iteration, we need to repeat NS times, as from the previous n1th iteration, we have NS solutions in Sn1. We elaborate the detailed process of the nth iteration in the procedure, which is shown in Fig. 4. (1) Inserting vessel Vn: When inserting Vn, vessels V1 to Vn1 are already in the plane. The locations of these n1 vessels in the plane are denoted by an array of n1 points, {(x1, y1), . . . , (xn1, yn1)}. There are NS arrays retained in each iteration, which means there are NS possible sets of locations for V1 to Vn1. In Fig. 4, Sn1 is the set of NS arrays (or solutions) retained in the previous iteration. As shown in Fig. 4, each solution in Sn1 is an array of n1 vessels’ locations {(x1, y1), . . . , (xn1, yn1)}. Based on each solution in Sn1, vessel Vn is inserted into the plane within its searching space. And we repeat the process for all the NS arrays (or solutions) in the set Sn1. In more details, given each array, i.e., n1 vessels’ locations {(x1, y1), . . . , (xn1, yn1)}, we insert Vn and enumerate all the feasible points (xn, yn) in its discretized searching space. Here a feasible point means Vn is placed at this point without overlaps with vessels V1 to Vn1. We repeat this inserting process for all the NS arrays (or solutions) in the set Sn1. (2) Greedy selection: As shown in Fig. 4, the above NS enumerations leads to a solution pool containing many arrays of n vessels’ locations {(x1, y1), . . . , (xn, yn)}, which are generated by enumerating all the feasible search spaces for all the NS arrays of n1 vessels’ locations. Among the solution pool, we select the top NG best ones according to the cost of baseline schedule, i.e.,  Pn  1 2 i¼1 ci  ðxi  ai Þ þ ci  jyi  bi j . (3) Look-ahead selection: Based on each of the obtained NG arrays of the n vessels’ locations {(x1, y1), . . . , (xn, yn)}, vessels Vn+1, Vn+2, . . . , Vn+DS are inserted into the plane. Here DS is a parameter which means ‘searching depth’, e.g., DS = 3 in this study. We   define hnk ; k 2 f1; 2; . . . ; NGg, as one of the NG arrays of the n vessels’ locations {(x1, y1), . . . , (xn, yn)}. The goodness of hnk ; G hnk , is evaluated by:

Fig. 4. The nth iteration in the procedure.

60

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68 nþDS X    G hnk ¼ Min c1i  ðxi  ai Þ þ c2i  jyi  bi j i¼1



ð5:1Þ



Subject to : xi þ t i 6 xj þ M 1  dxij 8i; j 2 fV 1 ; . . . ; V nþDS g; i – j;   y 8i; j 2 fV 1 ; . . . ; V nþDS g; i – j; yi þ li 6 yj þ M 1  dij dxij

þ

dyij

þ

dxji

þ

dyji

ð5:2Þ ð5:3Þ

P 1 8i; j 2 fV 1 ; . . . ; V nþDS g; i – j;

ð5:4Þ

yi þ li 6 L 8i 2 fV nþ1 ; . . . ; V nþDS g;

ð5:5Þ

xi P ai 8i 2 fV nþ1 ; . . . ; V nþDS g; Decision variables : xi ; yi ; P 0 8i 2 fV nþ1 ; . . . ; V nþDS g;

ð5:6Þ ð5:7Þ

dxij ; dyij 2 f0; 1g 8i; j 2 fV 1 ; . . . ; V nþDS g; i – j:

ð5:8Þ

Since (5.1) contains the absolute value form, we linearize it and rewrite the model as: nþDS X    G hnk ¼ Min c1i  ðxi  ai Þ þ c2i  ðhþi þ hi Þ

ð5:9Þ

i¼1

Subject to : Constraints ð5:2—5:8Þ; yi  bi ¼ hþi  hi ;

ð5:10Þ

hþi ; hi P 0:

ð5:11Þ

Since (x1, y1), . . . , (xn, yn) are the known data in the above model, the number of decision variables is not large if DS is reasonable. Thus the solving process is fast.   According to G hnk , we select the top NS best hnk from the set of NG hnk ; k 2 f1; 2; . . . ; NGg. Here, hnk is an array of the n vessels’ locations. These NS obtained arrays of n vessels’ locations {(x1, y1), . . . , (xn, yn)}, denoted by Sn in Fig. 4, will be the basis for the next iteration. The core idea behind the look-ahead selection lies in: when inserting a vessel, we consider the influence on its following DS vessels instead of just simply selecting the best locations for the vessel within its searching space. This procedure inserts all the N vessels into the time-berth plane one by one according to a given sequence. After inserting all the N vessels, choose the best solution from SN as the baseline schedule. In this study, we set NG = 100, NS = 30, and DS = 3. These parameters could be adjusted manually. The increase of these parameters could enlarge possibility of finding better near-optimal solutions, but it makes the solving process become more time-consuming. The proper parameter setting depends on the computer’s capacity. If we set NG = 1, NS = 1 and DS = 0, the proposed procedure will be degenerated to the normal greedy search, i.e., inserting vessels at their best locations within the remaining space one by one. This way is fast but incurs much loss of optimality. 5.3. Obtain recovery schedules under scenarios   þ   This section addresses how to obtain recovery schedule xþ i ðxs Þ; xi ðxs Þ; yi ðxs Þ; yi ðxs Þ for the scenario xs with a baseline schedule (xi, yi) given. The model is as follows:

Minimize

N X 

 Dþ 1 D 2þ Dþ 2 D c1þ i  xi ðxs Þ þ c i  xi ðxs Þ þ c i  yi ðxs Þ þ c i  yi ðxs Þ

ð5:12Þ

i¼1

s:t:

xi þ xþi ðxs Þ  xi ðxs Þ þ t i ðxs Þ 6 xj þ xþj ðxs Þ  xj ðxs Þ þ Mð1  dxij ðxs ÞÞ 8i; j 2 V; i – j; yþi ð

yi ð

yþj ð

yj ð

dyij ð

xs Þ  xs Þ þ li 6 yj þ xs Þ  xs Þ þ Mð1  xs ÞÞ 8i; j 2 V; i – j; xs Þ þ dyij ðxs Þ þ dxji ðxs Þ þ dyji ðxs Þ P 1 8i; j 2 V; i – j; li 6 yi þ yþi ðxs Þ  yi ðxs Þ þ li 6 L 8i 2 V; xi þ xþi ðxs Þ  xi ðxs Þ P ai ðxs Þ 8i 2 V; xi  ai þ xiDþ ðxs Þ  xDi  ðxs Þ ¼ xi þ xþi ðxs Þ  xi ðxs Þ  ai ðxs Þ 8i 2 V; jyi  bi ðxs Þj þ yiDþ ðxs Þ  yiD ðxs Þ ¼ jyi þ yþi ðxs Þ  yi ðxs Þ  bi ðxs Þj 8i 2 V; xþi ðxs Þ; xi ðxs Þ; yþi ðxs Þ; yi ðxs Þ 6 Mðai ðxs Þ  sðxs ÞÞþ 8i 2 V; xþi ðxs Þ; xi ðxs Þ; yþi ðxs Þ; yi ðxs Þ; xiDþ ðxs Þ; xiD ðxs Þ; yiDþ ðxs Þ; yDi  ðxs Þ P 0 8i 2 V: yi þ

dxij ð

ð5:13Þ ð5:14Þ ð5:15Þ ð5:16Þ ð5:17Þ ð5:18Þ ð5:19Þ ð5:20Þ ð5:21Þ

The above model is drawn out from the second stage of M0 . The notations for the variables and the explanation for the constraints are the same as Section 4. In this study, the above model will be solved by CPLEX. However, it may be time consuming in large-scale problems. Thus we propose some tactic to quicken the solving process by reducing the number of decision variables. For the vessel Vi 2 U(xs), which is defined in Section 4, we define its Adjacent Area AðV i Þ ¼ fðx; yÞjx 2 ½ai ðxs Þ; ai ðxs Þ þ a  t i ; y 2 ½yi  b  li ; yi þ b  li g. a and b are parameters to determine the impacting range in the time and position axes. For vessels þ   V k 2 fV n jðxn ; yn Þ 2 [V i 2Uðxs Þ AðV i Þg, their xþ k ðxs Þ; xk ðxs Þ; yk ðxs Þ; yk ðxs Þ are decision variables in the model; and for other vessels, these variables are fixed at zero.

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

61

5.4. Solution enhancement According to the proposed method in the above two subsections, we obtain a solution which includes: a baseline schedule (xi, yi) and a þ   series of recovery schedules xþ i ðxs Þ; xi ðxs Þ; yi ðxs Þ; yi ðxs Þ for all scenarios xs 2 X. However, it is a two-stage solving strategy, which may cause some loss of optimality of the original problem ðM0Þ which optimizes the baseline schedule and recovery schedules simultaneously. Therefore, we need to enhance the obtained solutions. When the numbers of vessels and scenarios are large, the original problem is too complex to be solved by CPLEX directly. Here we could solve a large scale problem of M0 with d parameters fixed, i.e., dxij ; dyij ; dxij ðxs Þ; dyij ðxs Þ. The computing complexity of M0 of mainly lies in the d-related constraints. If these d parameters are fixed, the d-related constraints will be simplified into normal constraints without the large number M, and the number of constrains will also be evidently reduced. The fixing of these d parameters is actually to fix the position relationships among vessels in the time-berth plane. For the obtained baseline schedule BS, and the recovery schedules RSs for scenarios xs 2 X, we calculate d the parameters: dxij ; dyij for BS, and dxij ðxs Þ; dyij ðxs Þ for all the RSs. According to these parameters, we simplify the d-related constraints in M0 as follows:

If dxij ¼ 1; If If If If

xi þ t i 6 xj

8i; j 2 V; i – j;

ð5:22Þ

dyij ¼ 1; yi þ li 6 yj i; j 2 V; i – j; dxij ð s Þ ¼ 1; xi þ xþi ð s Þ  xi ð s Þ þ t i ð s Þ 6 xj þ xþj ð s Þ  xj ð s Þ dyij ð s Þ ¼ 1; yi þ yþi ð s Þ  yi ð s Þ þ li 6 yj þ yþj ð s Þ  yj ð s Þ i; j y y x x dij ; dij ; dij ð s Þ; dij ð s Þ ¼ 0; these related constraints are omitted:

8

x x

x x

x

x x

x

x

x

x

x

8

8i; j 2 V; i – j; 8xs 2 X;

ð5:23Þ ð5:24Þ

2 V; i – j; 8xs 2 X:

ð5:25Þ

x

Moreover, the constraints dxij þ dyij þ dxji þ dyji P 1, and dxij ðxs Þ þ dyij ðxs Þ þ dxji ðxs Þ þ dyji ðxs Þ P 1 could also be ignored. Then M0 with the above simplification could be solved by CPLEX. In this way, both the baseline schedule and the recovery schedules are improved simultaneously. It should be mentioned that the above d-fixed enhancement method may not improve the solution substantially, because it does not enable changing the relative positions of vessels in the time-berth plane, which may need to be changed by the means of adjusting the sequence of vessels. However this approach is practical as the decision made during the recovery phase should not deviate too much from the original plan or else it will be too disruptive. 5.5. Neighborhood set of a sequence Section 5.2 addresses a heuristic to obtain a near-optimal baseline schedule with a sequence of vessels given. Pairwise exchange (i.e., swapping two randomly selected vessels in the sequence) is often used to obtain neighbor solutions in BAP (Kim and Moon, 2003). If the number of vessels is N, the size of neighborhood is N (N  1)/2. If N is large, the whole procedure will become extremely time-consuming due to the large number of neighbors. Different from the strategy of random exchange, we refine the selection process of neighbors to reduce the size of neighborhood. In the heuristic mentioned in Section 5.2, the vessels in the front of sequence will have more opportunities to occupy better spaces than the ones in the rear. Here, occupying a better space means the vessel would have lower operation cost. In order to strengthen the overall performance, weak performing elements are assigned with higher priority in the solution process by moving them toward the front of the inserting sequence, which borrows some core idea from the method of SWO (Meisel and Bierwirth, 2009). For a baseline schedule BS, i.e., {(x1, y1), . . . , (xN, yN)}, it is obtained based on a inserting sequence seq0. We assume the size of neighborhood is R, and it is also the temperature length R in simulated annealing. The R neighbors of seq0, i.e., seqn, n 2 {1, 2, . . . , R}, are obtained by following steps: (1) calculate the cost of each vessel, i.e., costðV i Þ ¼ c1i  ðxi  ai Þ þ c2i  jyi  bi j; (2) calculate cost gap between adjacent vessels in the sequence seq0 : Gap(i) = cost(Vi+1)  cost(Vi), here Vi means the ith vessel in the sequence; (3) select the vessels Vi with the largest Gap(i), and exchange Vi and Vi+1 so as to obtain a new inserting sequence seq1. In the same way, we obtain seq2, . . . , seqR with the second largest, . . ., the Rth largest Gap(i). Among these R newly obtained sequences, if one or more have been generated in a previous iteration, we will ignore these repeated ones and generate more sequences with the (R+1)th largest Gap(i) and so forth. 6. Computational experiments The computational experiments in this study include two parts: the performance analysis on the proposed meta-heuristic, and the numerical investigation on the proposed BAP model M0. 6.1. Performance analysis on the proposed meta-heuristic This section is the performance analysis on the proposed meta-heuristic. It contains two subsections: Section 6.1.1 concerns the comparison with CPLEX; and Section 6.1.2 concerns the analysis on convergence, parameters tuning, and CPU time of the proposed metaheuristic. 6.1.1. Comparison between CPLEX and the meta-heuristic approach We conduct some experiments to compare the results obtained by the proposed meta-heuristic and the optimal objective value solved by CPLEX directly. The mathematic model is implemented by CPLEX11.0 with concert technology of C# (VS2008) on a PC (Intel Core 2 Duo, 2.67 GHz; Memory, 4G). Table 1 shows the results on 20 small-scale instances, in which 8–20 vessels are considered. In Table 1, N and S are the number of vessels and scenarios. For the proposed meta-heuristic approach, we set T = 40, r = 0.6. From Table 1, it is noted that the proposed meta-heuristic approach could obtain optimal or near-optimal solutions with comparing the results calculated by CPLEX. The gap between the two results is not evident, and the average gap value is only 1.3%. It implies that the pro-

62

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

Table 1 Performance comparison between CPLEX and the proposed meta-heuristic. Problems

CPLEX

Meta-heuristic

GAP (%)

Case_ID

N

S

OBJC

CPU time

OBJH

CPU time

OBJH OBJ C OBJ G

c8_1 c8_2 c8_3 c8_4 c10_1 c10_2 c10_3 c10_4 c12_1 c12_2 c12_3 c12_4 c15_1 c15_2 c15_3 c15_4 c15_5 c15_6 c20_1 c20_2

8 8 8 8 10 10 10 10 12 12 12 12 15 15 15 15 15 15 20 20

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 8 8 5 5

206 256 452 371 172 291 242 276 264 250 283 302 225 125 124 394 LB = 265 LB = 309 LB = 653 LB = 518

17 46 1703 1071 61 1147 324 601 42 97 1132 359 473 114 82 2081 Out of Out of Out of Out of

206 263 452 378 178 292 242 282 264 254 289 304 225 130 124 403 295 318 713 552

42 80 45 47 47 52 49 51 135 139 112 127 164 146 157 151 185 191 248 226

0 2.7 0 1.9 3.5 0.3 0 2.2 0 1.6 2.1 0.7 0 4.0 0 2.3 – – – –

memory memory memory memory

Average gap (%)

1.3

Fig. 5. The evolution process and convergence of the meta-heuristic.

Fig. 6. The impacts of parameters (T & r) on objective values.

posed meta-heuristic approach is an effective method for solving the proposed BAP model M0. As to the computing time of two methods, the CPU time of CPLEX method varies very obviously even for the instances with the same problem scale (e.g., number of vessels, number of scenarios). Whereas the CPU time of the meta-heuristic method does not vary evidently, and it mainly depends on the problem scale. From Table 1, we could observe that CPLEX could not solve some instances when the number of vessels or the number of scenarios exceeds some limits, but the meta-heuristic approach could solve them in the reasonable computing time. 6.1.2. Performances of the meta-heuristic approach 6.1.2.1. The evolution process and convergence. We use a problem instance as the example to illustrate the convergence of the proposed meta-heuristic. This instance includes 40 vessels and 10 scenarios. We force the program to run 30 generations. The size of neighborhood

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

63

Fig. 7. The CPU time of the proposed meta-heuristic approach.

for solutions is 10, thus there are 10 solutions obtained in each generation. We plot the minimum, average and maximum objective values of each generation in Fig. 5. The convergence of the minimum objective value is evident. 6.1.2.2. The impacts of parameters on the objective value. The initial temperature T and the cooling rate r are two important parameters in the proposed approach. It needs to find a suitable combination of T and r for the proposed method. Five different problem instances with 40 vessels are tested. As to each problem instance, we solve it with five different initial temperatures T (20, 30, 40, 50, and 60) and four different cooling rates r (0.5, 0.6, 0.7, and 0.8). Thus these five instances are solved 20 times respectively. We calculate the average values of five instances as to these 20 combinations (T and r) respectively. Based on the 20 average values, we plot the minimum, average, and maximum values as to five different initial temperatures T and four different cooling rates r in Fig. 6, respectively. It is observed that the best performance is found when T = 40, r = 0.6. Thus these values are used in the following experiments. 6.1.2.3. The CPU time of the proposed meta-heuristic approach. It is necessary to investigate the computing speed (CPU time) of the proposed meta-heuristic approach. The CPU time is influenced by many aspects, among which the number of vessels and the number of scenarios may be the main factors. Here four different numbers of vessels (N = 20, 30, 40, and 50) and four different numbers of scenarios (S = 5, 10, 15, and 20), i.e., 16 combinations (N and S), are analyzed. For each combination, five problem instances are generated and solved by the proposed meta-heuristic approach. We calculate the average values of five instances as to these 16 combinations respectively. Based on the 16 average values, we plot the minimum, average, and maximum values as to four different numbers of vessels N and four different numbers of scenarios S in Fig. 7, respectively. From Fig. 7, it is observed that the average CPU time increases with theN and S grows. The increase trends are different. The impact of the number of scenarios on CPU time is not as significant as the number of vessels. This experimental result shows that the number of vessels may be an important factor of influencing the proposed meta-heuristic’s computation time. 6.2. Numerical investigation on the proposed BAP model The second part of experiments is the numerical investigation on the proposed BAP model. It has two sections: Section 6.2.1 concerns the numerical investigation on small-scale problems by using CPLEX; and Section 6.2.2 concerns the large-scale problems by using the proposed meta-heuristic. 6.2.1. Numerical investigation on small-scale problems Experiments are conducted under varied arrival time and operation time respectively. In these experiments, three different methods are compared. Method 1: It is the proposed model M0. The optimal objective value is denoted by Z1, which is the sum of the baseline schedule cost P and the expected recovery cost in realistic scenarios. Z 1 ¼ CðbslÞ þ x2X fpðxÞC rec ðx; bslÞg, here C(bsl) is the cost of baseline schedule bsl, p(x) is the probability of scenario x, and Crec(x, bsl) is the recovery cost in scenario x based on bsl. Method 2: According to the estimated arrival time ai and operation time ti, we solve the deterministic BAP and obtain a baseline schedule bsl0 whose cost measure is C(bsl0 ). Based on bsl0 , we calculate the expected value of recovery costs under scenarios, P P 0 0 0 x2X fpðxÞC rec ðx; bsl Þg. The total cost of Method 2 is denoted by Z 2 ¼ Cðbsl Þ þ x2X fpðxÞC rec ðx; bsl Þg. The major difference between Method 1 and Method 2 lies in: Method 1 solves a baseline schedule and recovery plans simultaneously; while Method 2 solves them in two steps. The solution obtained by Method 2 is also a feasible solution for the model M0 whose optimal objective is Z1. So Z1 is no greater than Z2. The gap between them evaluates the value of stochastic solution, which measures the cost of ignoring uncertainties while making a decision (Avriel and Williams, 1970).

Val Stochas ¼ Z 2  Z 1 :

ð6:1Þ

Method 3: According to the actual arrival time ai(x) and operation time ti(x), we solve a series of deterministic BAPs, each of which is  related with a particular scenario. For example, as to scenario x, the optimal schedule    is schx . We use the expected value of P these optimal schedules’ costs as the final cost of Method 3: Z 3 ¼ x2X pðxÞC schx . We can see Method 3 does not apply (or exist) in reality, because its assumption lies in the port operator can predict the actual arrival time and operation time in advance. Thus there is no recovery process in Method 3. However, Method 3 can supply us a lower bound for Method 1. P P For Method 1, Z 1 ¼ CðbslÞ þ x2X fpðxÞC rec ðx; bslÞg. As x2X pðxÞ ¼ 1, we have:

64

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

Z1 ¼

X

fpðxÞ  ½CðbslÞ þ C rec ðx; bslÞg:

x2X

In scenario x, the baseline schedule bsl is adjusted to be another schedule schx after recovery process. Due to c1þ P c1i P c1 and i i 2 2 P ci P ci , the sum of bsl’s cost and the adjustment cost in recovery process is no less than the cost of schedule schx,  C(bsl) + Crec(x, bsl) P C(schx). In addition, the cost of the schedule schx is no less than the optimal solution schx by solving the deterministic  BAP that is related with the scenario x. So CðbslÞ þ C rec ðx; bslÞ P Cðschx Þ. Hence: c2þ i

Z1 ¼

X

fpðxÞ  ½CðbslÞ þ C rec ðx; bslÞg P

x2X

X   pðxÞ  Cðschx Þ ¼ Z 3 :

x2X

The gap between the Z1 and Z3 evaluates the value of perfect information:

Val Info ¼ Z 3  Z 1 :

ð6:2Þ

As aforementioned, the key assumption of Method 3 is that the port operator can predict the actual arrival time and operation time in advance. The gap that Z3 is lower than Z1, just reflects the value of perfect information (here the ‘perfect’ means to exactly predict vessels’ arrival time and operation time in advance). It measures the maximum amount that port planners are willing to pay for knowledge value of random variables (e.g., vessels’ arrival time, operation time) before making their decision (Avriel and Williams, 1970). As to BAP under stochastic variables, jVal_Infoj may mean the cost of estimating exact operation time or the cost of influencing the vessels to commit themselves to their estimated arrival time. In this case, port planners can schedule the plan optimally without adjustments (or recovery process) during plan’s execution. A large jVal_Infoj means the randomness plays an important role in BAP, thus stochastic factors should not be ignored.

Table 2 Results of small-scale problems with varied arrival time. Problems

Method 3

Val_Stochas

N

S

Vvar

p(xs)

Mvar

Z1

Method 1 C M1 bsl

T

Z2

C M2 bsl

Z3

Z2  Z1

Gap Z1

Case 1 8 vessels

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2

0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15

2 2 2 3 3 3 4 4 4 2 2 2 3 3 3 4 4 4

103.0 110.0 118.0 110.0 117.8 129.5 115.0 123.8 135.1 116.0 125.6 136.5 148.0 160.0 165.8 166.2 179.2 198.8 136.6

68 68 103 71 71 71 68 71 79 68 68 84 88 88 88 98 101 101

2 4 2 3 3 3 2 2 7 19 33 56 557 884 3600 243 1420 1799

103.0 110.0 120.5 110.0 118.4 131.0 115.0 124.4 138.5 116.0 125.6 140.0 148.0 164.0 188.0 169.0 189.2 219.5 140.6

68 68 68 68 68 68 68 68 68 68 68 68 68 68 68 68 68 68

71.2 71.8 72.8 74.9 76.3 78.4 75.0 76.4 78.5 77.9 79.9 82.9 84.0 87.2 92.0 83.8 87.0 91.7 80.1

0 0 2.5 0 0.6 1.5 0 0.6 3.4 0 0 3.5 0 4.0 22.2 2.8 10.0 20.7

5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15

1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3

58 58 58 72 72 72 67 67 67 58 58 58 72 72 72 67 67 67

4 3 2 12 9 22 53 9 49 6 5 8 6 11 13 22 7 7

71.2 73.8 77.8 97.4 105.3 117.1 104.2 113.4 127.3 71.2 73.8 77.8 97.4 108.6 115.4 100.4 115.7 134.2 99.0

58 58 58 58 58 58 58 58 58 58 58 58 58 58 58 58 58 58

58.4 58.5 58.6 65.4 66.9 69.1 68.2 70.2 73.3 57.0 56.8 56.5 63.3 64.4 66.0 66.1 67.7 70.2 64.3

0 0 0 8.4 12.9 19.6 4.2 7.8 12 0 0 0 8.4 12.2 20.9 1.2 7.9 10.8

Average of case 2

71.2 73.8 77.8 89.0 92.4 97.5 100.0 105.6 115.3 71.2 73.8 77.8 89.0 96.4 94.5 99.2 107.8 123.4 92.0

Average of two cases

114.3

Average of case 1 Case 2 10 vessels

Method 2

119.8

72.2

Val_Info Z3  Z1

Gap Z1

0 0 2 0 1 1 0 0 3 0 0 3 0 3 13 2 6 10 2

31.8 38.2 45.2 35.1 41.5 51.1 40.0 47.4 56.6 38.1 45.7 53.6 64.0 72.8 73.8 82.4 92.2 107.1

31 35 38 32 35 39 35 38 42 33 36 39 43 46 45 50 51 54 40

0 0 0 9 14 20 4 7 10 0 0 0 9 13 22 1 7 9 7

12.8 15.3 19.2 23.6 25.5 28.4 31.8 35.4 42.0 14.2 17.0 21.3 25.7 32.0 28.5 33.1 40.1 53.2

18 21 25 27 28 29 32 34 36 20 23 27 29 33 30 33 37 43 29

5

(%)

(%)

34

Notes: (1) S: the number of scenarios in which some vessel’s arrival time deviates from the estimated one; it does not include a scenario, x0, in which all vessels arrive punctually; (2) Vvar: the number of varied vessels in each scenario; (3) p(xs): the probability of each scenario in which some vessel’s arrival time deviates from estimated one; so in the above table, you may notice that p(xs)  S < 1; this gap, i.e., 1  p(xs)  S is actually the p(x0); (4) Mvar: the magnitude of variation, i.e., the average of deviations of the M2 actual arrival times to the estimated times in all scenarios; (5) C M1 bsl : the cost of baseline schedule in Method 1; (6) T: CPU computing time; (9) C bsl : the cost of baseline schedule in Method 2; (7) Val_Stochas: the value of stochastic solution; (8) Val_Info: the value of perfect information.

65

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68 Table 3 Results of small-scale problems with varied operation time. Problems

Method 1

Method 2

Method 3

Val_Stochas

N

S

Vvar

p(xs)

Mvar

Z1

C M1 bsl

Case 1 8 vessels

5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2

0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15

2 2 2 3 3 3 4 4 4 2 2 2 3 3 3 4 4 4

137.2 151.0 164.8 176.2 194.6 215.3 206.0 233.2 263.0 198.2 207.8 215.8 271.0 279.0 289.5 326.6 342.7 353.4 234.7

68 68 100 84 84 116 70 70 137 138 156 190 162 237 237 168 300 300

4 11 9 14 41 26 17 46 47 1578 590 313 1912 407 615 3600 584 161

137.2 151.0 171.8 177.2 199.0 231.8 209.2 237.4 279.8 258.2 296.2 353.3 287.2 331.0 396.8 352.2 409.0 443.3 273.4

68 68 68 68 68 68 68 68 68 68 68 68 68 68 68 68 68 68

77.8 79.8 82.7 79.4 81.7 85.1 80.1 82.5 86.2 84.5 87.8 92.8 94.2 99.4 107.3 102.8 109.8 120.2 90.8

0 0 7.0 1.0 4.4 16.5 3.2 4.2 16.8 60.0 88.4 137.5 16.2 52.0 107.3 25.6 66.3 89.9

5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15 0.10 0.12 0.15

1 1 1 2 2 2 3 3 3 1 1 1 2 2 2 3 3 3

58 58 58 72 72 72 72 72 91 63 63 63 72 72 72 72 72 93

3 3 8 23 59 40 61 168 325 3 8 5 128 46 104 210 321 140

91.2 97.8 107.8 147.4 165.3 192.1 179.2 203.4 239.8 96.2 103.8 115.3 152.4 171.3 199.6 189.2 215.4 254.8 162.3

58 58 58 58 58 58 58 58 58 58 58 58 58 58 58 58 58 58

60.5 61.0 61.8 66.8 68.6 71.2 69.6 71.9 75.4 60.5 61.0 61.8 66.8 68.6 71.2 69.6 71.9 75.4 67.4

0 0 0 13.4 18.9 27.1 7.2 11.4 19.8 1.2 2.4 4.3 8.4 12.9 19.6 2.2 5.4 13

Average of case 2

91.2 97.8 107.8 134.0 146.4 165.0 172.0 192.0 220.0 95.0 101.4 111.0 144.0 158.4 180.0 187.0 210.0 241.8 153.0

Average of two cases

193.9

Average of case 1 Case 2 10 vessels

T

Val_Info

Z2

C M1 bsl

Z3

Z2  Z1

Gap Z1

217.9

79.1

Z3  Z1

Gap Z1

0 0 4 1 2 8 2 2 6 30 43 64 6 19 37 8 19 25 15

59.4 71.2 82.1 96.8 112.9 130.2 125.9 150.7 176.8 113.7 120.0 123.0 176.8 179.6 182.2 223.8 232.9 233.2

43 47 50 55 58 60 61 65 67 57 58 57 65 64 63 69 68 66 60

0 0 0 10 13 16 4 6 9 1 2 4 6 8 11 1 3 5 6

30.7 36.8 46.0 67.2 77.8 93.8 102.4 120.1 144.6 34.5 40.4 49.2 77.2 89.8 108.8 117.4 138.1 166.4

34 38 43 50 53 57 60 63 66 36 40 44 54 57 60 63 66 69 53

10

(%)

(%)

56

Notes: (1) the problem cases are the same as the ones in Table 1; (2) S: the number of scenarios in which some vessel’s operation time deviates from the estimated one; (3) Mvar: the magnitude of variation, i.e., the average of deviations of the actual operation times to the estimated times in all scenarios.

In contrast to Method 3, Method 1 is primarily based on a baseline schedule and secondarily on a fluctuation of scenarios. In addition, Method 1 assumes primarily one specific scenario for the baseline schedule and considers other ones with implicitly much less probabilities in experiments. Tables 2 and 3 illustrate the experiment results with varied arrival time and operation time, respectively. Z1 is no greater than Z2 in Tables 2 and 3. It validates the proposed model can generate a berth allocation schedule that is better at handling the uncertainties than the one obtained by the traditional deterministic BAP. In Tables 2 and 3, C M1 bsl denotes the cost of the baseline schedule obtained by the pro  posed BAP model, and it is no less than the cost C M2 of schedules obtained by the deterministic BAP model. However, with considering bsl the adjustment costs in recovery process, the total cost Z1 of the proposed model is no greater than Z2. This result demonstrates: (1) the best schedule, which is obtained by deterministic model according to the estimated arrival time and operation time, may not be the best one when facing uncertainties; (2) the proposed BAP model could generate the best schedule that minimize the sum of the baseline schedule cost and recovery cost in scenarios simultaneously. Method 2 is actually to minimize these two parts (baseline schedule cost & recovery cost) independently. Besides the above implications, from Tables 2 and 3, we could also observe some phenomena: (1) The impact of probability of scenarios in which some vessels’ arrival or operation time varies. In this experiment, three different probabilities p(xs) for each scenario are used, 0.1, 0.12, and 0.15. Z1, Z2, and Z3 increase with the probability growing. The value of stochastic solution and perfect information is also increasing. (2) The impact of the magnitude of variation for arrival or operation time. The variation for arrival (or operation) time represents how much time a vessel is later (or longer) than its estimated arrival (or operation) time. In this experiment, the magnitude of variation Mvar is the average value for all vessels in all scenarios. Here, three different magnitude of variation Mvar is used: 2, 3, 4 for the case 1 (8 vessels) and 1, 2, 3 for the case 2 (10 vessels). Z1, Z2, and Z3 increase with this magnitude of variation growing. The value of perfect information is also mounting up, but the value of stochastic solution has no positive or negative correlation with Mvar.

66

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

Table 4 Results of large-scale problems with varied arrival time or varied operation time. Problems

Experiments with varied arrival time Costs of methods S

Mvar

10 10 10 10 10 10 10 10 10 10 6 7 8 9 10 11 12 13 14 15

3 4 5 6 7 8 9 10 11 12 8 8 8 8 8 8 8 8 8 8

684 743 773 768 807 756 789 753 744 766 776 789 765 759 756 731 700 692 692 687 747

10 10 10 10 10 10 10 10 10 10 6 7 8 9 10 11 12 13 14 15

3 4 5 6 7 8 9 10 11 12 8 8 8 8 8 8 8 8 8 8

Costs of methods

Improvement

(Z2  Z1)/Z1 (%)

(Z3  Z1)/Z1 (%)

Z1

Z2

Z3

(Z2  Z1)/Z1 (%)

(Z3  Z1)/Z1 (%)

684 781 808 838 832 814 816 774 752 779 830 873 830 839 814 795 714 718 745 753 789

827 906 1053 1087 1065 1158 1154 1216 1235 1281 1172 1216 1194 1120 1158 1084 1119 1021 997 1013 1104

0 5 5 9 3 8 3 3 1 2 7 11 8 11 8 9 2 4 8 10 6

21 22 36 42 32 53 46 61 66 67 51 54 56 48 53 48 60 48 44 47 48

748 765 776 827 909 918 936 1031 1076 1110 956 980 960 964 918 939 786 816 905 811 907

748 799 852 929 991 1008 1021 1039 1076 1158 956 980 995 996 1008 1030 850 875 940 840 955

998 1086 1106 1257 1324 1476 1431 1508 1620 1577 1469 1534 1487 1483 1476 1362 1376 1298 1247 1278 1370

0 4 10 12 9 10 9 1 0 4 0 0 4 3 10 10 8 7 4 4 5

33 42 43 52 46 61 53 46 51 42 54 57 55 54 61 45 75 59 38 58 51

Average of case 2

1341 1387 1407 1419 1465 1498 1447 1466 1544 1525 1487 1506 1449 1448 1498 1388 1375 1326 1304 1286 1407

1436 1412 1567 1586 1618 1638 1591 1576 1634 1587 1642 1671 1616 1679 1638 1574 1546 1416 1438 1377 1560

1843 1824 1969 2064 2106 2177 2162 2458 2502 2528 2418 2456 2234 2269 2177 2064 2106 1964 1978 1887 2155

7 2 11 12 10 9 10 8 6 4 10 11 12 16 9 13 12 7 10 7 9

37 32 40 45 44 45 49 68 62 66 63 63 54 57 45 49 53 48 52 47 51

1479 1506 1534 1586 1638 1675 1697 1745 1736 1769 1603 1648 1683 1669 1675 1642 1599 1606 1568 1553 1631

1503 1568 1614 1761 1841 1849 1864 1837 1896 1934 1656 1739 1766 1803 1849 1855 1774 1721 1675 1634 1757

2065 2128 2320 2234 2271 2426 2356 2613 2646 2768 2837 2889 2754 2698 2426 2386 2351 2247 2278 2191 2444

2 4 5 11 12 10 10 5 9 9 3 6 5 8 10 13 11 7 7 5 8

40 41 51 41 39 45 39 50 52 56 77 75 64 62 45 45 47 40 45 41 50

Average of two cases

1087

1175

1630

8

49

1269

1356

1907

7

50

Average of case 1 Case 2 40 vessels

Z2

Improvement Z3

Case 1 40 vessels

Z1

Experiments with varied operation time

(3) The impact of number of scenarios. In Tables 2 and 3, S is the number of scenarios in which some vessels vary their arrival time or operation time. In this experiment, we set the number of scenarios S with 5 or 6 for the case 2. It does not show that S obviously impacts the total cost value. As to the value of perfect information and the value of stochastic solution, none direct correlation is observed with S. (4) The comparison between the impacts of varied arrival time and operation time. The problem cases and parameters setting are the same for experiments with varied arrival time and operation time. So we could make a comparison between the impacts of these two varied factors. From the average data in Tables 2 and 3, we could observe that the uncertain factor of operation time may have more positive influence on Z1, Z2, and Z3 than the uncertain factor of arrival time. Moreover, as to the value of stochastic solution and perfect information, the average of experiments with varied operation time is evidently larger than the ones with varied arrival time. 6.2.2. Numerical investigation on large-scale problems As to large-scale problems, the numerical investigation is conducted by using the proposed meta-heuristic. The number of vessels is 40. The results with varied arrival time or operation time are shown in the left and right parts in Table 4, respectively. Here, three methods are compared. Method 1: It is the proposed model, which is a two-stage optimization problem to minimize the total cost of baseline schedule and the expected value of recovery costs for scenarios. The best objective value is denoted by Z1 in Table 4. Method 2: It is the same as the ‘Method 2’ mentioned in Section 6.2.1. Firstly, the best baseline schedule is obtained according to the estimated time. Then the recovery adjustments are obtained according to the actual time in scenarios. The sum of the two costs is denoted by Z2 in Table 4.

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

67

Fig. 8. The influence of parameters on the performances.

Method 3: It is to perform the recovery process by manual adjustments. For example, if vessel i is late in arrival time or its operation time is prolonged, planners will simply postpone the vessels that succeed vessel i and occupy some berth space of vessel i. This is the ‘right-shift’ strategy which is commonly used in job scheduling. Firstly, the best baseline schedule is obtained according to the estimated arrival time and operation time. Then, the above manual adjustments are performed according to the actual arrival time and operation time in scenarios. The sum of the baseline schedule cost and recovery cost is denoted by Z3 in Table 4. In Table 4, two series of experiments are performed with the number of scenarios (S) changing, or the magnitude of variation (Mvar) for arrival time (or operation time) changing. From the results, the proposed Method 1 obviously outperforms the two other ones. The percentages of the improvement are about 5%–9% by comparing Method 1 to Method 2, and about 48%–51% by comparing Method 1 to Method 3. Experiments in Table 4 use the same problem cases and parameters setting (S and Mvar). So we could compare the impacts of these two varied factors: arrival time and operation time. From the results, we could find the costs of Method 1 (Z1) in the left part of Table 4 are smaller than the right part. It demonstrates the impact of varied operation time on the final performance is evidently larger than varied arrival time. More specifically, it means: the event (denoted by E1) in which the operation time of vessel i is prolonged for 10 more hours will incur more additional costs than the event (denoted by E2) in which vessel i is late in arrival for 10 hours. In the case of E2, some empty ‘time-berth’ space spared out by vessel i’s lateness in arrival will be used for other vessels which may moor at better berthing positions or earlier time. That is the reason for why the varied operation time causes worse performances than the varied arrival time under the same magnitude of variation. The reason for why Z2 in the left part of Table 4 is smaller than the right part is similar as the above. We plot the average value of results of two problem cases in Fig. 8 to investigate the influence of different parameters on the final performances. In Fig. 8(1), the variation of arrival time (Mvar) is from 3 to 12, and the number of scenarios is 10. Here the variation of arrival time represents how much time a vessel is later than its estimated arrival time; hence Mvar is the average value for all vessels in all scenarios. The costs of Method 1 and Method 2 are not changing evidently with the variation of arrival time growing. The cost of Method 3 is obviously growing with the variation increasing. It shows the increase of lateness for the arrival time does not incur the recovery costs in Method 1 and Method 2. The reason for this phenomenon is similar as the above analysis. Some ‘berth-time’ space of the late vessel will be spared for other vessels that may berth at better positions or earlier time, which will decrease the recovery costs. In Fig. 8(3), the variation of operation time (Mvar) is from 3 to 12, and the number of scenarios is 10. Different from Fig. 8(1), the costs of Method 1 and Method 2 are increas-

68

L. Zhen et al. / European Journal of Operational Research 212 (2011) 54–68

ing evidently with the variation of operation time growing. It demonstrates the increase of prolongation for the operation time will cause growing recovery costs in Method 1 and Method 2. Fig. 8(2) and (4) is the data with the changing number of scenarios (S) from 6 to 15, and the variation of arrival time or operation time (Mvar) is set at 8. From Fig. 8(2) and (4), we could see costs of Method 1 and Method 2 are not changing obviously with S varying. Though S grows, the probability for each scenario is decreasing; the total probability for scenarios that contain vessels with varied arrival time or operation time always remains at 0.5. Therefore, S has no evident influence on the total cost. 7. Conclusions This paper develops a two-stage decision model for berth allocation problem (BAP) under uncertainties. The model is designed to balance the cost associated with the initial (baseline) schedule and the expected cost of deviation from the initial schedule. By comparing with other scholars’ work in this area, the major contribution mainly includes two aspects: (1) Most of researches in BAP concern how to obtain an optimal baseline schedule in a static and deterministic environment. This paper makes an exploratory study on BAP under uncertainties. A two-stage decision model is developed. Some computational investigations are conducted to analyze the performance of berth allocation process under uncertainties. (2) The deterministic BAP is NP-hard, not to mention the BAP under uncertainties. We propose a meta-heuristic approach for solving the above problem in large-scale realistic environments. Numerical experiments are conducted to validate the effectiveness of the proposed method. However, there are some limitations for the current method. In this study, we use the re-schedule strategy to handle the conflicts between the initial plan and the realistic environments. The deterministic BAP is already NP-hard, thus the re-schedule problem on BAP is even more complex. This study is the integration of these two problems, so the proposed model seems a bit complex. In fact, port operators could also use compression strategy when these conflicts are not serious. Here the compression strategy means reducing the vessels’ operation time by allocating more resources, e.g., quay cranes, prime movers. In this way, the stochastic berth allocation problem could become easier to handle. Then we could conduct studies on BAP with stochastic arrival time that could follow arbitrary probability distributions. The scenario-based method used in this study could be circumvented. That will be our research direction in the future. Acknowledgements The authors thank three anonymous referees for their suggestions which have helped to improve this paper. This research is supported by the MPA (Maritime Port Authority, Singapore) and NOL (Neptune Orient Lines) fellowships. This work was done when the first author was at the Department of Industrial & Systems Engineering, National University of Singapore. References Avriel, M., Williams, A.C., 1970. The value of information and stochastic programming. Operations Research 18, 947–954. Bierwirth, C., Meisel, F., 2010. A survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research 202 (3), 615– 627. Cordeau, J.F., Laporte, G., Legato, P., Moccia, L., 2005. Models and tabu search heuristics for the berth allocation problem. Transportation Science 39, 526–538. Guan, Y., Cheung, R.K., 2004. The berth allocation problem: Models and solution methods. OR Spectrum 26, 75–92. Hansen, P., Oguz, C., Mladenovic, N., 2008. Variable neighborhood search for minimum cost berth allocation. European Journal of Operational Research 191, 636–649. Imai, A., Nagaiwa, K., Chan, W., 1997. Efficient planning of berthing allocation for container terminals in Asia. Journal of Advanced Transportation 31, 75–94. Imai, A., Nishimura, E., Papadimitriou, S., 2001. The dynamic berth allocation problem for a container port. Transportation Research Part B 35, 401–417. Imai, A., Sun, X., Nishimura, E., Papadimitriou, S., 2005. Berth allocation in a container port: Using a continuous location space approach. Transportation Research Part B 39, 199–221. Imai, A., Nishimura, E., Hattori, M., Papadimitriou, S., 2007. Berth allocation at intended berths for mega-containerships. European Journal of Operational Research 179, 579– 593. Imai, A., Chen, H., Nishimura, E., Papadimitriou, S., 2008. The simultaneous berth and quay crane allocation problem. Transportation Research Part E 44, 900–920. Kim, K.H., Moon, K.C., 2003. Berth scheduling by simulated annealing. Transportation Research Part B 37, 541–560. Kim, K.H., Park, Y.M., 2004. A crane scheduling method for port container terminals. European Journal of Operational Research 156, 752–768. Legato, P., Mazza, R., 2001. Berth planning and resources optimization at a container terminal via discrete event simulation. European Journal of Operational Research 133, 537–547. Lim, A., 1998. The berth planning problem. Operations Research Letters 22, 105–110. Meisel, F., Bierwirth, C., 2009. Heuristics for the integration of crane productivity in berth allocation problem. Transportation Research Part E 45, 196–209. Moorthy, R., Teo, C.P., 2006. Berth management in container terminal: The template design problem. OR Spectrum 28, 495–518. Murty, K.G., Liu, J., Wan, Y.W., Linn, R., 2005. A decision support system for operations in container terminal. Decision Support Systems 39, 309–332. Nishimura, E., Imai, A., Papadimitriou, S., 2001. Berth allocation planning in the public berth system by genetic algorithms. European Journal of Operational Research 131, 282–292. Park, K.T., Kim, K., 2002. Berth scheduling for container terminals by using a sub-gradient optimization technique. Journal of the Operational Research Society 53, 1054–1062. Park, Y.M., Kim, K.H., 2003. A scheduling method for berth and quay cranes. OR Spectrum 25, 1–23. Wang, F., Lim, A., 2007. A stochastic beam search for the berth allocation problem. Decision Support Systems 42, 2186–2196.