Robust schedule design for liner shipping services

Transportation Research Part E 48 (2012) 1093–1106

Contents lists available at SciVerse ScienceDirect

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

Robust schedule design for liner shipping services Shuaian Wang, Qiang Meng ⇑ Department of Civil and Environmental Engineering, National University of Singapore, Singapore 117576, Singapore

a r t i c l e

i n f o

Article history: Received 8 December 2011 Received in revised form 4 April 2012 Accepted 17 April 2012

Keywords: Scheduling Liner shipping Schedule robustness Port time Sample average approximation

a b s t r a c t This paper examines the design of liner ship route schedules that can hedge against the uncertainties in port operations, which include the uncertain wait time due to port congestion and uncertain container handling time. The designed schedule is robust in that uncertainties in port operations and schedule recovery by fast steaming are captured endogenously. This problem is formulated as a mixed-integer nonlinear stochastic programming model. A solution algorithm which incorporates a sample average approximation method, linearization techniques, and a decomposition scheme, is proposed. Extensive numerical experiments demonstrate that the algorithm obtains near-optimal solutions with the stochastic optimality gap less 1.5% within reasonable time. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Liner shipping companies try to provide regular shipping services with fixed schedules. They publish the information in advance up to a period of 3–6 months on port rotations, arrival and departure dates at each port of call, and sometimes also the ships that are deployed. For example, the ship routes provided by global liner shipping companies Maersk (2011) and OOCL (2011) are available in their websites. A regular and reliable shipping schedule is important for shippers/consignees to arrange their production, transportation and distribution plans. Nevertheless, liner shipping services are not as reliable as one might imagine. According to a survey by Drewry Shipping Consultants (2006) based on no less than 5410 ship arrivals, 21% of ships arrive 1 day behind schedule and 22% arrive two or more days behind schedule. The low reliability of liner schedules can be explained by a number of factors, such as congestion at ports, fluctuation of container handling time, and cascading effects from previous ports of call. The schedule unreliability of liner shipping services would pose substantial losses for shippers/consignees. Consider an automotive manufacturing company which has components shipped from overseas on a weekly basis. If the delivery of the components is delayed and the manufacturing company does not keep enough inventories, the production line may have to stop. This will further have cascading effects on the delivery of the end products (cars) to distributors and consumers. As is acknowledged by Psaraftis (2004), former CEO of the Piraeus Port Authority, ‘‘the name of the game of all major container lines is their ability to meet their schedules, as they incur enormous costs, both real and intangible, in case they do not’’. According to a survey conducted by Notteboom (2006), 93.6% of delayed schedules are attributable to port access and terminal operations. To improve the reliability of liner schedules, port operators can enhance their capacity by expanding the berthing area and deploying more quay cranes, or by employing more sophisticated berth and yard planning systems (see, e.g., Vis and de Koster, 2003; Steenken et al., 2004; Stahlbock and Voß, 2008; Bierwirth and Meisel, 2010; Zhen et al., 2011a,b). However, the rising container volumes and capacity constraints in many ports around the world mean that berth availability on arrival at a port is not always guaranteed. Therefore, in this paper we examine how to design robust ⇑ Corresponding author. Tel.: +65 6516 5494; fax: +65 6779 1635. E-mail address: [email protected] (Q. Meng). 1366-5545/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tre.2012.04.007

1094

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

schedules for liner ship routes to hedge against the uncertainties arising in port operations. Here we define robust schedule as schedule that is designed while considering the adverse effects of late arrival as a result of uncertainties in the shipping operations. In fact, it is widely accepted that including some buffer time (contingencies) in the schedule can improve the schedule integrity. However, liner shipping companies are reluctant to do so because buffer time is very expensive (Notteboom, 2006): the daily cost of a ship amounts to tens of thousands of dollars. The objective of this paper is to design robust schedules for liner ship routes by determining the optimal trade-off between buffer time and schedule robustness (reliability, integrity, stability). In the literature on optimization of liner shipping operations, ship route schedules and uncertainties at ports are seldom considered (see Ronen, 1983, 1993; Christiansen et al., 2004, 2007; for reviews). Existing research efforts are mainly devoted to liner ship fleet planning (e.g., Meng and Wang, 2011a), liner ship route design (e.g., Shintani et al., 2007; Reinhardt and Pisinger, 2012), liner shipping network design (e.g., Fagerholt, 1999, 2004; Alvarez, 2009; Jepsen et al., 2011; Meng and Wang, 2011b), liner ship fleet deployment (e.g., Gelareh and Pisinger, 2011), container ship speed optimization and emission (e.g., Ronen, 2011; Song and Xu, 2012), container routing (e.g., Brouer et al., 2011), container fleet sizing (e.g., Dong and Song, 2009, 2012), empty container repositioning (e.g., Song and Zhang, 2010; Song and Dong, 2011a,b), and liner shipping network analysis (e.g., Bell et al., 2011). These studies usually assume that the port time is fixed or is a linear function of the number of containers handled. Port congestion and port time variability are not taken into account. There are a few research efforts on the analysis of liner schedules. Notteboom (2006) assessed the causes of schedule unreliability: 93.6% are attributable to port operations, 5.3% are due to weather conditions or mechanical problems at sea, 0.9% are because of delay when transiting canals, and 0.2% results from unexpected wait time at bunkering sites. Vernimmen et al. (2007) examined the adverse effects of schedule unreliability on shippers/consignees. Bell and Bichou (2008) investigated how to avoid ship bunching through improving the handling efficiency at ports whereas they did not consider the schedule recovery by increasing sailing speed. Yan et al. (2009) investigated the ship scheduling problem at the operational level without taking into account port time uncertainty. Wang and Meng (2011) analysed the design of ship route schedules by assuming fixed sailing speeds at sea and fixed time at ports. Wang and Meng (2012a) further considered variable sailing speed and stochastic port time. However, late arrival at a port is not allowed in their model. Qi and Song (2012) examined a schedule design problem considering the probability of on-time arrivals. In this seminal work, they obtained the structural properties of the optimal schedule under certain conditions with useful managerial insights regarding the impact of port uncertainties. According to the literature review, the robust schedule design problem which takes into account the penalty for late arrival/late container handling resulted from uncertain port time has practical significance and is worth research efforts. The objective of this paper is to design a schedule while capturing the uncertainties in port operations, which include the uncertain wait time due to port congestion and the uncertain container handling time. The designed schedule is robust in that the uncertainties in port operations and schedule recovery by fast steaming are captured endogenously. To address the problem, we formulate a mixed-integer nonlinear stochastic programming model. To solve the model, we propose a solution algorithm that incorporates a sample average approximation method, linearization techniques, and a decomposition scheme. Extensive numerical experiments demonstrate that the algorithm obtains near-optimal solutions with the stochastic optimality gap less 1.5% within reasonable time. The remainder of this paper is organized as follows. In addition to describing notation and assumptions, Section 2 introduces the concept of ship route schedule, schedule recovery, and the trade-off between schedule robustness and total costs. Section 3 develops a mixed-integer nonlinear stochastic programming model for the robust ship route schedule design problem. Section 4 proposes a solution methodology for the robust schedule design problem. Section 5 carries out numerical experiments based on an Asia–America–Europe ship route to assess the computational performance of the solution method. Conclusions are presented in Section 6. 2. Problem description A ship route has fixed port rotations where the ports of call are visited in a predetermined sequence. For example, an Asia–America–Europe (AAE) ship route provided by a global liner shipping company, as shown in Fig. 1, can be represented

Fig. 1. An Asia–America–Europe (AAE) ship route.

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

1095

by its port calling sequence: Tokyo (1) ? Kobe (2) ? Chiwan (3) ? Hong Kong (4) ? Kaohsiung (5) ? Busan (6) ? Kobe (7) ? Tokyo (8) ? Balboa (9) ? Manzanillo (10) ? Miami (11) ? Jacksonville (12) ? Savannah (13) ? Charleston (14) ? New York (15) ? Antwerp (16) ? Felixstowe (17) ? Bremerhaven (18) ? Rotterdam (19) ? Le Havre (20) ? New York (21) ? Norfolk (22) ? Charleston (23) ? Manzanillo (24) ? Balboa (25) ? San Pedro (26) ? Oakland (27) ? Tokyo (1). In fact, we can choose any port of call as the first one. Let I denote the total number of ports of call and define set I :¼ f1; 2; . . . ; Ig. The voyage between two consecutive ports of call on a ship route is referred to as a leg. The ith leg is defined as the voyage from the ith port of call to the (i + 1) th port of call when i = 1, 2, . . . , I  1 and the Ith leg is from the Ith port of call to the 1st port of call. Nearly all the ship routes operated by global liner shipping companies have weekly services, that is, each port of call is visited on the same day every week. For example, if the round-trip journey time is 63 days, then nine ships are required to maintain the weekly service. 2.1. Ship route schedule and time components in a round-trip journey Liner shipping companies need to publish the arrival date at each port of call. Since shippers/consignees are more concerned about the actual start handling (discharging and loading of containers) time, it makes no difference for shippers/consignees whether a ship has arrived at a port but has to wait for berth, or the ship arrives late. To be more precise, we let ti (hours) be the published start handling time at the ith port of call (the time point 0 can be arbitrarily chosen). We assume that ships will not start handling containers before the published time ti. If the actual start handling time yi is later than ti, the liner shipping company may be subject to real and intangible costs. We assume that the cost function Fi(yi  ti) has the following expression:

F i ðyi  ti Þ ¼



C i =Ci  ðyi  t i Þ; yi  t i 6 Ci 8i 2 I C i ; y i  t i > Ci

ð1Þ

A schematic illustration of Fi(yi  ti) is shown in Fig. 2. The cost increases linearly when the delay does not exceed Ci. Ci can be e.g. 1 day. When the delay is longer than Ci, the cost remains at a constant value Ci. The rationale behind using such a piecewise-linear cost function is that shippers/consignees might just wait for the delivery of their cargo when the delay is not significant. Therefore the cost increases linearly with the delay. When the delay exceeds a threshold Ci, shippers/consignees must take measures such as emergent purchasing and production rescheduling, and hence the cost does not change any more. In general, Ci is larger at the ports of call with more discharged containers, such as major US ports. It should be noted that cost functions other than that in Eq. (1) can easily be accommodated in the subsequent formulations. Even if a ship arrives at the ith port of call at time ti, it may not immediately start berthing for container handling due to ~ i (hours) the random wait time at the ith port of call. w ~ i is related to the berthing strategy port congestion. We represent by w of the port. h For example, i the wait time is shorter if the company has dedicated berths or terminals. Let wi be a realization of ~ i ; wi 2 W min . To hedge against port congestion, ships try to arrive earlier than ti. Let ui be the company’s target w ; W max i i arrival time at the ith port of call. Whereas ti is published, ui is solely for the company’s internal use. It is reasonable to design ui and ti satisfying

W min 6 t i  ui 6 W max i i

8i 2 I

ð2Þ

Let xi be the actual arrival time at the ith port of call. As we assume that ships will not start handling containers before the published time ti, the actual start handling time yi can be calculated by

yi ¼ maxfxi þ wi ; ti g 8i 2 I

ð3Þ

~ (hours), is also random due to the variability associated The container handling time at the ith port of call, denoted by h i with the number of quay cranes deployed, the efficiencies of quay crane operators, the number of containers handled, and h i ~ ; h 2 Hmin ; Hmax . After discharging and loading containthe storage plan of containers on ships. Let h be a realization of h i

i

i

i

i

ers, a ship departs from the port. Let zi be the actual departure time from the ith port of call. It follows that

Fig. 2. Cost function for delayed container handling.

1096

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

zi ¼ yi þ hi

8i 2 I

ð4Þ

Fig. 3 plots a few examples of the time components and time points in a round trip. The circle represents the port rotations in a round trip: a ship returns to the 1st port of call after visiting the last port of call. Fig. 3 shows that at the 1st port of call (Tokyo), the ship arrives exactly at the target arrival time, namely, x1 = u1. The wait time w1 is shorter than the time interval between the published start handling time and the target arrival time, that is, w1 < t1  u1. Hence, the ship still waits for a period of t1  u1  w1 hours even when the berth is available. The ship starts to discharge and load containers at time y1 = t1 and the handling time is h1. After that, the ship leaves the port at time z1 = t1 + h1. At the 2nd port of call (Kobe), the ship arrives at the target arrival time u2. However, the wait time is too long and the actual start handling time y2 is later than the published start handling time t2. At the third port of call (Chiwan), the ship arrives later than the target arrival time, that is, x3 > u3, and the start handling time is also delayed, namely, y3 > t3. Let uI+1 be the target arrival time at the 1st port of call after a round trip journey. tI+1, xI+1, yI+1, and zI+1 are defined in a similar ^ , satisfies manner. To ensure the weekly service frequency, the number of ships deployed on the ship route, denoted by n

^ tIþ1  t1 ¼ uIþ1  u1 ¼ 168 h=week  n

ð5Þ

We further define vector

^ ; t i ; ui Þi2I k ¼ ðn

ð6Þ

The robust schedule design problem aims to determine the vector k to balance the trade-off between ship cost, bunker cost, and late start handling cost. 2.2. Schedule recovery After a ship leaves the ith port of call at time zi, it tries to sail at a speed such that the actual arrival time xi+1 matches the target arrival time ui+1 at the next port of call. Let Li (n mile) be the oceanic distance of the ith leg and hence the sailing speed vi (knot) should be Li/(ui+1  zi). In case of delayed departure, either because of port congestion, or long container handling time, or cascading effects from previous ports of call, the ship has to sail at a very high speed to catch up with the target arrival time ui+1. However, the sailing speed vi cannot exceed the designed maximum speed Vmax (knot) as a result of engine speed and power. If Vmax < Li/(ui+1  zi), even if the ship sails at its highest speed, it still cannot catch up with the target arrival time ui+1. Therefore, we require that if Vmax P Li/(ui+1  zi), the sailing speed vi should equal Li/(ui+1  zi). Otherwise the sailing speed vi = Vmax. It is reasonable to assume that in the best situation, that is, when a ship starts handling containers on time and the container handling time takes its minimum value Hmin , the ship should be able to arrive at the target arrival time at the next i port of call. Therefore we have

Fig. 3. Time components in a round-trip journey.

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

  Li = uiþ1  ti  Hmin 6 V max i

8i 2 I

1097

ð7Þ

It is evident that a longer leg with larger Li is more capable of absorbing the delay from previous ports of call. As a result, when the company designs the schedule, it requires that ships must be able to arrive at the target arrival time at the port of call following a long leg. In the example of AAE, legs 27, 8, 15, and 20 are such long legs. Hence, the target arrival times at the 1st (Tokyo), 9th (Balboa), 16th (Antwerp), and 21st (New York) ports of call must be maintained. In fact, these ports of call are all the 1st calling ports in their respective regions. Let K be the number of the 1st calling ports in regions and define K :¼ f1; 2; . . . ; Kg. We further define Ik to be the calling sequences corresponding to these ports of call. That is, I1 = 1, I2 = 9, I3 = 16, I4 = 21 for AAE. The voyage between the ports of call Ik and Ik+1 is referred to as segment k, and the port of call Ik is defined as the start port of segment k. To simplify the subsequent mathematical representation, we choose one of the start port of a segment as the 1st port of call in the itinerary, for example, Tokyo in AAE. Since it is the time intervals ui+1  ui rather than each time point ui that actually determine the schedule, without loss of generality, we can stipulate

u1 ¼ 0

ð8Þ xmax i

zmax i

Define IK+1 :¼ I + 1. Let and represent the latest arrival and departure time at the ith port of call, respectively. We assume that ships will not arrive at the (i + 1)st port of call earlier than the planned time ui+1. It is evident that the following equations hold

xmax ¼0 1 max þ W max þ Hmax zIk ¼ xmax Ik Ik Ik

8k 2 K

¼ xmax þ W max þ Hmax zmax i i i i

8k 2 K;

ð9Þ ð10Þ i ¼ Ik þ 1; . . . ; Ikþ1  1

max xmax þ Li =V max g 8i 2 I iþ1 ¼ maxfuiþ1 ; zi

ð11Þ ð12Þ

To ensure that ships arrive at each start port of a segment at the target arrival time, that is, xIk ¼ uIk for each k 2 K, we have

xmax Ikþ1 6 uIkþ1

8k 2 K

ð13Þ

2.3. Trade-offs in robust schedule design In the robust schedule design problem, the vector k is the here-and-now decision vector and xi, yi, zi and vi are wait-andsee decision variables. That is, contrary to k, the variables xi, yi, zi and vi are determined according to the realization of the ~ i ; i 2 I . To highlight the dependency of x , y , z and v on the random parameters w ~i , we ~ i and h ~ i and h uncertain parameters w i i i i ~ ¼ ðh ~ Þ , and use x ðw; ~ y ðw; ~ z ðw; ~ and v ðw; ~ in place of x , y , z and v , respec~ ¼ ðw ~ ~ ~ ~ ~ i Þi2I and h define vectors w hÞ; hÞ; hÞ, hÞ i i2I i i i i i i i i ~j Þ, where j 6 i, and both leg j ~ j; h tively. It should be noted that the values of xi, yi, zi and vi only depend on the sample values ðw ~ only to simplify the notation. We further define vector ~ hÞ and leg i are on the same segment. We use the whole vector ðw;

~ ¼ ðx ðw; ~ ~ z ðw; ~ ~ ~ hÞ ~ hÞ; ~ hÞÞ hðw; i ~ hÞ; yi ðw; i ~ hÞ; v i ðw; i2I

ð14Þ

The trade-offs in schedule design is as follows. If the schedule is tight, namely, ti  ui is small and the sailing speeds on each leg are high, then the time required to finish a round-trip journey is short. Consequently, fewer ships are required to maintain a weekly service according to Eq. (5). At the same time, the bunker cost is high because the daily bunker consumption is approximately proportional to the third power of sailing speed (Notteboom, 2006). In case of delay at one port, schedules at the subsequent ports of call are prone to be unreliable because ships already sail at very high speeds and thus the sailing time cannot be significantly shortened. The situation is opposite when the schedule is slack. Therefore, the robust ~ for each realization of the ~ hÞ schedule design problem aims to determine the vector k while considering the vector hðw; ~ so as to minimize the total expected costs, including ship costs, bunker costs, and late start han~ hÞ uncertain parameter ðw; dling costs. 3. Mixed-integer nonlinear stochastic programming model Let copr (USD/week) represent the fixed operating cost of a ship deployed on the ship route, a (USD/ton) be the bunker price, and fi(vi) be the bunker consumption rate (ton/n mile) at the speed vi on the ith leg. Since the daily bunker consumption is proportional to the third power of sailing speed, we can write fi(vi) = ai(vi)2 as fi(vi) is expressed in terms of tons per nautical mile. The auxiliary decision variables are combined into a vector





p ¼ tIþ1 ; uIþ1 ; xmax ; zmax i i

 i2I

; xmax Iþ1



ð15Þ

The robust schedule design (RSD) problem can be formulated as a mixed-integer nonlinear stochastic programming model:

1098

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

X ~ þ F i ðy ðw; ~  t i ÞÞ ~ hÞÞ ~ hÞ ^ þ E ða  Li  fi ðv i ðw; ½RSD c ¼ mincopr  n i k;p

ð16Þ

i2I

s.t.:

~ ¼ y ðw; ~ þh ~ 8i 2 I ~ hÞ ~ hÞ zi ðw; i i ~ ~ þw ~ hÞ ¼ maxfxi ðw; ~ hÞ ~ i ; ti g 8i 2 I yi ðw; ~ ~ V max g 8i 2 I ~ ~ hÞÞ; v i ðw; hÞ ¼ minfLi =ðuiþ1  zi ðw; ~ ¼0 ~ hÞ x ðw;

ð17Þ

~ ¼ z ðw; ~ ~ 8i 2 I n fIg ~ hÞ ~ hÞ xiþ1 ðw; i ~ hÞ þ Li =v i ðw;

ð21Þ

W min i

6 t i  ui 6 8i 2 I   min Li = uiþ1  ti  Hi 6 V max 8i 2 I

ð22Þ

u1 ¼ 0

ð24Þ

1

W max i

xmax ¼0 1

ð18Þ ð19Þ ð20Þ

ð23Þ ð25Þ

zmax ¼ uIk þ W max þ Hmax Ik Ik Ik

8k 2 K

ð26Þ

¼ xmax þ W max þ Hmax 8k 2 K; i ¼ Ik þ 1; . . . ; Ikþ1  1 zmax i i  i max i  max xiþ1 ¼ max uiþ1 ; zi þ Li =V max 8i 2 I

ð27Þ

xmax Ikþ1 6 uIkþ1

ð29Þ

8k 2 K

^ tIþ1  t1 ¼ uIþ1  u1 ¼ 168n ^ positive integers n

ð28Þ ð30Þ ð31Þ

The objective function (16) minimizes the expected weekly total costs. The first term is the ship costs and the second term is ~ y ðw; ~ v ðw; ~ ~ hÞ; ~ hÞ; the expected bunker costs plus late start handling costs. Eqs. (17)–(21) define the variables zi ðw; i ~ hÞ, and i ~ respectively. Eq. (22) states the relationship between the target arrival time and the published start handling time. ~ hÞ, xi ðw; Eq. (23) imposes that in the best situation, that is, when a ship starts handling containers on time and the container handling time takes its minimum value, the ship should be able to arrive at the target arrival time at the next port of call. Eq. (24) states that the target arrival time at the 1st port of call is time 0. Eqs. (25)–(29) enforce that ships are able to arrive at the target arrival time at the start port of each segment. Eq. (30) imposes the weekly service frequency constraint. Eq. (31) defines that the number of ships takes a positive integer value. 4. Solution methodology The mixed-integer nonlinear stochastic programming model RSD can hardly be addressed by conventional optimization methods or solvers. In this section, we first apply the sample average approximation method to approximate the expected value function in the objective (16). We subsequently adopt several linearization techniques to linearize the nonlinear constraints and the nonlinear objective function. Finally, we propose a decomposition scheme to improve the computational efficiency. 4.1. Sample average approximation ~ has only a small number of realizations (scenarios), it is usually impossible to solve ~ hÞ Unless the uncertain parameter ðw; model RSD exactly. A possible approach to address the difficulty is the sample average approximation (SAA) method. The SAA method solves stochastic optimization problems using Monte Carlo simulation. In this technique the objective function of the stochastic program is approximated by a sample average estimate derived from a random sample. The resulting sample average approximating problem is then solved by deterministic optimization approaches. This process is repeated with different samples to obtain a good candidate solution along with the statistical estimate of its optimality gap (see e.g. Mak et al., 1999; Verweij et al., 2003). ~ according to its probability distribution. Let S be the set of the S ~ hÞ To solve model RSD, we first generate S observations of ðw; ~ in scenario s 2 S is denoted by ws ; hs  . The decision variables x ðw; ~ ~ z ðw; ~ ~ hÞ ~ hÞ; scenarios. The realization of ðw; i ~ hÞ; yi ðw; i ~ hÞ, i i2I i s s s s ~ ~ and v i ðw; hÞ corresponding to scenario s 2 S is denoted by xi ; yi ; zi and v i , respectively. Define vector

  hS :¼ xsi ; ysi ; zsi ; v si i2I ;s2S

ð32Þ

Then the expected value objective function (16) is approximated with the sample average function. The approximating deterministic model is:

^þ ½SAA cS ¼ mincopr  n k;p

    1 XX a  Li  fi v si þ F i ysi  ti S s2S i2I

ð33Þ

1099

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

s.t.: s

zsi ¼ ysi þ hi 8s 2 S; 8i 2 I   ysi ¼ max xsi þ wsi ; ti 8s 2 S; 8i 2 I    max  s s v i ¼ min Li = uiþ1  zi ; V 8s 2 S;

ð34Þ ð35Þ

8i 2 I

ð36Þ

xs1 ¼ 0 8s 2 S

ð37Þ

xsiþ1 ¼ zsi þ Li =v si

8s 2 S;

8i 2 I n fIg

ð38Þ

and constraints (22)–(31). The optimal value cS to model SAA is actually a random variable depending on the set S. The expected value of cS is no ~ each ~ hÞ, greater than c⁄, namely, E½cS  6 c (Mak et al., 1999). Consequently, we can generate N independent samples of ðw; of size S, and obtain N optimal objective values of model SAA, denoted by cnS ; n ¼ 1; 2; . . . ; N. A statistical lower bound for c⁄   P can be estimated by cS :¼ Nn¼1 cnS =N. Note that cS is also a random variable. Let Var cnS be the sample variance of     n cS ; n ¼ 1; 2; . . . ; N. When N is large (e.g., 20), cS can be considered as normally distributed Normal E½cS ; Var cnS =N (Strictly qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  n ffi speaking, ðcS  E½cS Þ= Var cS =N has a t-distribution with N  1 degrees of freedom). Therefore, in practice we can consider qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   LB :¼ cS  3  Var cnS =N as a lower bound for the optimal value c⁄ of model RSD according to the 3-sigma rule. A total of N schedule variables k are obtained after solving the N SAA models with different samples. We can choose the one with the lowest objective value cnS , denoted by k⁄, for deriving an upper bound as follows. First, a new sample S0 , whose size denoted by S0 is much larger than S, is generated. Then we compute the cost, denoted by cs0 , with fixed schedule variables k⁄ for eachpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi scenario s0 2 S 0 . Let Varðcs0 Þ be the sample variance of cs0 ; s0 2 S 0 . Similar to the lower bound, P UB :¼ cs0 þ 3 Varðcs0 Þ=S0 can be considered as an upper bound for c⁄ where cs0 :¼ s0 2S0 cs0 =S0 . 4.2. Linearization Model SAA has nonlinear objective function (33) and nonlinear constraints 28, 35, 36, 38 which are difficult to tackle. We linearize them by taking advantage of the problem structure, introducing indicating binary variables, reformulating decision variables and approximating nonlinear functions. It is easy to see that the ‘‘=’’ in constraints (28) and (35) can be replaced by ‘‘P’’. Hence, Eq. (28) can be replaced by

xmax 8i 2 I iþ1 P uiþ1 max þ Li =V max xmax iþ1 P zi

8i 2 I

ð39Þ ð40Þ

ysi P xsi þ wsi 8s 2 S; 8i 2 I ysi P ti 8s 2 S; 8i 2 I

ð41Þ ð42Þ

Eq. (35) can be replaced by

To linearize Eqs. (36) and (38), we first take the reciprocal of sailing speed v si , denoted by usi , as the new decision variable, namely,

usi ¼ 1=v si

8s 2 S;

8i 2 I

ð43Þ

Now Eq. (38) becomes

xsiþ1 ¼ zsi þ Li usi

8s 2 S;

8i 2 I n fIg

ð44Þ

and Eq. (36) becomes

usi ¼ max





uiþ1  zsi =Li ; 1=V max

Introducing binary decision variable

usi usi usi



P uiþ1 



=Li

8s 2 S;

s bi

8s 2 S;

8i 2 I

ð45Þ

2 f0; 1g; s 2 S; i 2 I and using the big-M method, Eq. (45) can be linearized as

8i 2 I

ð46Þ

max

P 1=V 8s 2 S; 8i 2 I    s s 6 uiþ1  zi =Li þ M1 1  bi 8s 2 S;

usi 6 1=V s bi

zsi



max

þ

s M 1 bi

2 f0; 1g 8s 2 S;

8s 2 S;

ð47Þ

8i 2 I

ð48Þ

8i 2 I

ð49Þ

8i 2 I

ð50Þ

s where M1 is a big number. Since bi is binary, Eqs. (46) and (47) are    max s s (48) and (49) are equivalent to ui 6 max . iþ1  zi =Li ; 1=V

u

equivalent to

usi

P max



uiþ1 

zsi

  =Li ; 1=V max and Eqs.

1100

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

The cost function Fi(yi  ti) is a piecewise linear concave function. We introduce new binary variables   s di 2 f0; 1g; s 2 S; i 2 I . The term F i ysi  ti is replaced with a new decision variable F si and the following constraints are added using the big-M method: s

F si P C i  M 2 di 8s 2 S; 8i 2 I    s 8s 2 S; F si P C i =Ci  ysi  ti  M2 1  di s di

2 f0; 1g 8s 2 S;

ð51Þ

8i 2 I

ð52Þ

8i 2 I

ð53Þ

    where M2 is a big number. Constraints (51) and (52) are equivalent to P min C i =Ci  ysi  t i ; C i . To linearize the bunker consumption function fi(vi), we first define gi(ui) :¼ fi(1/ui) = ai(ui)2. Hence, gi(ui) is a nonlinear, monotonically decreasing, nonnegative, and convex univariate function, as shown in Fig. 4, where gi represents the bunker consumption on the ith leg. Given an approximation error eLA (the subscript LA means linear approximation), we can generate a piecewise-linear function for approximating the nonlinear function gi(ui) as shown in Fig. 4. It should be highlighted that the approximation approach in Fig. 4 outperforms the outer approximation approach proposed by Wang and Meng (2012b) because the former needs fewer approximation lines. The piecewise-linear function can be efficiently obtained by the bisection search method. As a consequence of the convexity of gi(ui), the approximating piecewise-linear function is also convex. Suppose that there are a total of Ki approximating lines in the piecewise-linear function for gi(ui): F si

j

j

g i ¼ slopei ui þ intercepti

j ¼ 1; 2; . . . ; Ki

ð54Þ

    Due to the convexity of the approximating piecewise-linear function, we can replace g i usi (which is fi v si Þ in the objective function by a new decision variable g si and add the following constraints: j

j

g si P slopei usi þ intercepti

8s 2 S;

8i 2 I ; j ¼ 1; 2; . . . ; Ki

ð55Þ

The linearized sample average approximation model is as follows:

½L-SAA capprox ¼ S

min

k;hS ;p;usi ;bsi ;dsi ;g si ;F si

^þ copr  n

 1 XX a  Li  g si þ F si S s2S i2I

ð56Þ

subject to constraints (22)–(27), (29)–(31), (34), (38)–(42), (44), (46)–(53), and (55). 4.3. Decomposition L-SAA is a mixed-integer linear programming model with a large number of integer variables and constraints due to the sample average approximation and linearization. We note that the time components on segment k, namely, between the ~ on other segments. The only complicating constraint ~ hÞ ports of call Ik and Ik+1, are not affected by the realizations of ðw; combining all segments is Eq. (30). Therefore, for any given time interval T k :¼ uIkþ1  uIk , we can evaluate the bunker cost (the second term in Eq. (56)) plus the delayed container handling cost (the third term in Eq. (56)) for segment k, denoted by b k ðT k Þ. After that we combine all the segments and obtain the optimal ship number n ^ . This decomposition scheme is elabC orated below. The theoretical minimum value T min for Tk can be derived from Eqs. (25)–(29). Whereas Tk can be infinitely large in theory, k b k ðT k Þ is a nonnegative decreasing it has an upper limit T max in practice for the following reasons. It is easy to see that C k

Fig. 4. Linear approximation of bunker consumption function.

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

1101

function. Given Tk, the number of ships required on the segment is between bTk/168c and dTk/168e. Hence, for any Tk no less than T min k , its total costs (ship costs, bunker costs, and late start handling costs) is no greater than

b k ðT k Þ þ dT k =168ecopr C

8k 2 K

ð57Þ

  b k T  þ dT  =168ecopr provides an upper bound on the total costs for segment k. Hence, given an which is no less than C k k On the other hand, the total costs with given Tk is no less than T min k ,

T k

b k ðT k Þ þ bT k =168ccopr C

8k 2 K

ð58Þ

b k ðT k Þ, the total cost is also no less than Due to the nonnegativity of C

bT k =168ccopr

8k 2 K

ð59Þ

Therefore, all the Tk fulfilling the following conditions can be excluded from consideration:

 b k ðT  Þ þ T  =168 copr bT k =168ccopr P C k k

8k 2 K

ð60Þ

As a result, we can simply set the upper limit

T max ¼ 168  k

! b k T   þ T  =168copr C k k þ 1 8k 2 K copr

ð61Þ

  b k ðT k Þ at many discrete points T min ¼ T 0 6    6 T 1 6    6 T Hk ¼ T max . Let C b 1 :¼ C b k T 1 ; 1 ¼ 0; 1; 2 . . . Hk . We We evaluate C k k k k k k k b use a piecewise linear function to approximate C k ðT k Þ, as shown in Fig. 5. Let Tk be a decision variable representing the time interval uIkþ1  uIk for segment k. We further introduce auxiliary continuous decision variables DT 1k and binary indication decision variables e1k . To obtain the optimal Tk for each segment and the ^ on the ship route, we can formulate L-SAA as a decomposed model: optimal ship number n

½DEC-L-SAA

capprox S

¼

min c

^ ;T k ;DT 1 ;e1 n k

k

opr

^þ n

X k2K

b0 þ C k

Hk b 1 b 11 X C C k

1¼1

T 1k 

k T k11

! 1

 DT k

ð62Þ

s.t.:

T k ¼ T min þ k

Hk X

DT 1k

8k 2 K

ð63Þ

1¼1

DT 1k 6 T 1k  T 0k

8k 2 K

Hk

DT k P 0 8 k 2 K   DT 1k P T 1k  T k11 e1k 8k 2 K; 1 ¼ 1; 2; . . . ; Hk  1   DT k1þ1 6 T k1þ1  T 1k e1k 8k 2 K; 1 ¼ 1; 2; . . . ; Hk  1 1

ek 2 f0; 1g 8k 2 K; X ^ T k ¼ 168n

1 ¼ 1; 2; . . . ; Hk  1

ð64Þ ð65Þ ð66Þ ð67Þ ð68Þ ð69Þ

k2K

^ positive integers n

ð70Þ

b 1 . Nevertheless, the number of decision variables After the decomposition, we have to solve many problems to obtain C k and constraints of each problem is about 1/K of the original L-SAA model. The decomposed model is preferable when K is large. 5. Numerical experiments We carry out numerical experiments to validate the effectiveness of the proposed solution approach. The mixed-integer linear programming models are solved by CPLEX-12.1 running on a 3 GHz Dual Core PC with 3 GB of RAM. 5.1. Dataset description The test instances in the numerical experiments are based on the AAE ship route. We choose this ship route because it is the most comprehensive ship routes (OOCL, 2011; Maersk, 2011). It covers two trade lanes: trans-Pacific and trans-Atlantic, and three regions: Asia, America, and Europe. In contrast, most ship routes either serve ports within a region, or serve one trade lane (two regions). The sequence of ports of call and the oceanic distance of each voyage leg Li (n mile) is as follows: Tokyo (1: 348) ? Kobe (2: 1353) ? Chiwan (3: 46) ? Hong Kong (4: 354) ? Kaohsiung (5: 898) ? Busan (6: 341) ? Kobe (7: 348) ? Tokyo (8: 7696) ? Balboa (9: 72) ? Manzanillo (10: 1165) ? Miami (11: 300) ? Jacksonville (12: 109)

1102

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

Fig. 5. Discretization and linear approximation of bunker and late handling cost function.

? Savannah (13: 52) ? Charleston (14: 611) ? New York (15: 3269) ? Antwerp (16: 136) ? Felixstowe (17: 327) ? Bremerhaven (18: 246) ? Rotterdam (19: 438) ? Le Havre (20: 3082) ? New York (21: 273) ? Norfolk (22: 402) ? Charleston (23: 1537) ? Manzanillo (24: 72) ? Balboa (25: 2898) ? San Pedro (26: 360) ? Oakland (27: 4547) ? Tokyo (1). We use the following parameters based on data provided by a global liner shipping company. 8000-TEU ships with Vmax = 26.5 knots are deployed on the ship route. The operating cost of a ship copr = 700,000 USD/week. The bunker consumption coefficient ai is randomly generated from [0.02/24, 0.03/24], that is, the daily bunker consumption at the speed of 20 knots is between 160 and 240 tons. Bunker price a = 800 USD/ton. The late delivery penalty parameters Ci (million USD) ~ i satisfies and Ci (hours) are randomly generated from uniform distribution [0.5, 1] and [12, 24], respectively. The wait time w ~ satisfies the truncated normal the truncated normal distribution Normal (6, 3) on the interval [0, 12]. The handling time h i distribution Normal (36, 3) on the interval [24, 48]. In the SAA method we choose S = 15, N = 20, S0 = 200. 5.2. Computational performance First, we examine the effect of the linear approximation (55) on the objective value (total costs). Given an approximation P P error eLA, the largest absolute error in the total costs is i2I aLi eLA and the relative error is i2I aLi eLA divided by the total costs. We randomly generate eight test instances with different values of ai, Ci and Ci, and solve each instance with different values of eLA 2 {103, 5  104, 104}. Table 1 reports the average number of linear pieces for approximating the bunker consumption on each leg and the relative error in total costs. We observe that the linear approximation approach in Fig. 4 requires only a small number of linear pieces to approximate the nonlinear bunker consumption functions. This means that the additional computational burden posed by Eq. (55) is marginal. At the same time, the accuracy of the solution is high. When eLA = 104, the relative error is less than 0.02%. In the sequel, we set eLA = 104 for all computations. The main optimality gap comes from the SAA method and the relative mixed-integer programming gap tolerance (MIP gap) associated with CPLEX. Based on a few trial experiments, we set the MIP gap at 0.5% so that it is roughly the same  as the  optimality gap of SAA. Table 2 reports for each of the above eight test instances the average value cS and variance Var cnS of the N = 20 optimal objective values of model SAA, the probabilistic lower bound (mean value minus three times the standard deviation), the average value cs0 and variance Varðcs0 Þ of the S0 = 200 realizations with a given schedule, the probabilistic upper bound (mean value plus three times the standard deviation), the relative optimality gap of the SAA method, the MIP gap, the sum of the two gaps, and the CPU time (min) for obtaining both the lower bound and the upper bound. According to the results, the relative optimality gap of the SAA method is very tight (less than 1%). The total relative optimality gap of the overall solution algorithm is less than 1.5%. Moreover, all instances can be solved within 1 h of CPU time. Therefore, the proposed solution approach is effective and efficient for solving practical problems. 5.3. Impacts of the number of deployed ships on schedule robustness The optimal number of ships deployed for the first test instance is 20. To investigate the impacts of the number of ships on ^ at different values and solve the schedule design problem. Table 3 the schedule robustness, we fix the number of ships n reports the percentage of ports of call where the actual start handling time yi is later than the published ti and the total costs ^ . According to the results, the percentage of late start handling decreases with the number of for different numbers of ships n ships deployed. This means that liner shipping companies can improve their schedule robustness by adding more ships to

1103

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106 Table 1 Results of different linear approximation errors for bunker consumption functions.

eLA = 103

ID

1 2 3 4 5 6 7 8

eLA = 5  104

eLA = 104

# Pieces per leg

Relative error (%)

# Pieces per leg

Relative error (%)

# Pieces per leg

Relative error (%)

17.6 17.0 17.3 17.3 17.1 17.3 17.4 17.1

0.1318 0.1333 0.1330 0.1311 0.1327 0.1308 0.1318 0.1333

24.6 24.0 24.1 24.3 24.1 24.3 24.4 24.0

0.0659 0.0666 0.0665 0.0656 0.0663 0.0654 0.0659 0.0667

54.0 52.5 53.1 53.3 52.9 53.4 53.6 52.6

0.0132 0.0133 0.0133 0.0131 0.0133 0.0131 0.0132 0.0133

Table 2 Relative optimality gap and CPU time (min) of the SAA method. ID

1 2 3 4 5 6 7 8

Lower bound (106)   cS Var cnS

LB

cs0

19.03 18.81 18.83 19.12 18.94 19.12 19.03 18.80

19.02 18.80 18.81 19.10 18.92 19.11 19.01 18.79

19.11 18.90 18.90 19.24 19.02 19.21 19.15 18.89

759 327 434 969 479 448 464 484

Upper bound (106) Varðcs0 Þ 16497 11899 11395 20015 12636 11661 18373 9926

(UB–LB)/LB (%)

MIP gap (%)

Total gap (%)

CPU time

0.65 0.66 0.60 0.86 0.65 0.67 0.88 0.65

0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

1.15 1.16 1.10 1.36 1.15 1.17 1.38 1.15

36 29 31 30 41 35 31 30

UB 19.14 18.92 18.93 19.27 19.05 19.23 19.18 18.91

their ship routes. It should be noted that 20 ships corresponds to the lowest total costs. Therefore adding more ships would result in more operating costs, which explains why most liner shipping companies do not include in their weekly schedules sufficient buffer time as a hedge against port delay. The schedule robustness is a differentiating factor of liner shipping services. For instance, MSC keeps time buffers relatively low while Maersk keeps sufficiently high time buffers to cope with unexpected disruptions (Notteboom, 2006). As a result, Maersk is reputed for its high schedule reliability and high freight rate. 5.4. Analytical analysis of the effects of container handling time variability To obtain analytical insights into the problem, we focus on a single voyage leg. Assume that the wait time at the port ~ i ¼ 0 and delayed container handling penalty is infinity. We focus on the handling operations at port i and the voyage w of leg i. To simplify the notation, we remove the subscript i. The voyage distance of leg i is L. We first assume that the handling time at port i is deterministic, denoted by h, and hence represent by v the planned sailing speed on leg i. The formulation for minimizing the total cost is

L minaLav 2 þ ðh þ Þcopr =168

v

v

ð71Þ

The first term in Eq. (71) is the bunker cost and the second term is the ship operating cost where copr/168 is the hourly operating cost of a ship. Taking the first order derivative, we obtain the optimal speed

v ¼



1=3 copr 332aa

ð72Þ

Eq. (72) indicates that in the above setting, the optimal sailing speed is not related to the voyage distance L or the handling time h. v⁄ increases with ship operating cost copr and decreases with bunker price a and bunker consumption rate coefficient a. The increase of v⁄ with (copr)1/3 is linear as a consequence of the cubic relationship. When the handling time is random, we minimize the total expected cost. Represent by t :¼ ti+1  ti the planned time inter~ depending on h. ~ The val between the arrival at port i and port i + 1. The sailing speed would be a random variable L=ðt  hÞ problem can be formulated as

" #

2 L 1 ~ þ t  hÞc ~ opr =168 ¼ minaL3 aE minE aLa þ ðh þ tcopr =168 ~ ~ 2 t t th ðt  hÞ ~ is uniformly distributed over [Hmin, Hmax]. Evidently, t > Hmax. Therefore To obtain a tractable result, we assume that h

ð73Þ

1104

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

Table 3 Impacts of ship number on schedule robustness. Number of deployed ships

19

20

21

22

23

24

% Late start handling Total costs (106 USD)

12.8% 19.20

10.1% 19.07

8.6% 19.11

6.4% 19.27

5.6% 19.52

3.5% 19.84

Table 4 The optimal planned time interval and the average speed.

r

L = 200

t⁄ v 

E

L = 2000

0

10

20

30

0

10

20

30

50.5 13.8

54.3 12.3

62.1 10.1

71.0 8.5

181.3 13.8

181.8 13.7

183.1 13.7

185.4 13.6

1

¼

~ 2 ðt  hÞ

1 ðt  Hmax Þðt  Hmin Þ

ð74Þ

 :¼ ðHmax  Hmin Þ=2. Since t > Hmax, t  H > r  . Eq. (74) can be reformulated as Define H :¼ ðHmax þ Hmin Þ=2 and r

E

1

¼

~ 2 ðt  hÞ

1 2 ðt  HÞ2  r

ð75Þ

Eq. (75) indicates that compared with deterministic handling time which takes the average value H, the variability of han would lead to a higher fuel consumption for the same planned time interval t. This means that the variability of dling time r handling time would increase the cost of liner shipping companies. From the viewpoint of port operators, providing a consistent level of service would attract more liner shipping companies. The optimal planned time t⁄ can be obtained by letting  . For example, when the first-order derivative of Eq. (73) equal 0. The optimal planned time t⁄ increases with r H ¼ 36; a ¼ 0:001 and all other parameters are the same as in Section 5.1, Table 4 shows the results with two different voy . Evidently, t⁄ increases with r  and the increase is more significant when L is smaller. This image distances and different r  is large and L is small. plies that liner shipping companies would design a large t⁄ (slacker schedules) when r   can be calculated by After having obtained the optimal planned time t⁄, t⁄ > Hmax, the average speed denoted by v

v  ¼ E

L ~ t h 

¼

L t   Hmin ln  t   Hmax 2r

ð76Þ

Eq. (76) indicates that a larger t⁄results in a slower average sailing speed, as is demonstrated in Table 4. This result is consistent with the practice: the computed sailing speed according to the published schedules of liner shipping companies is usually very small for short voyage legs (see, e.g., Maersk, 2011; OOCL, 2011). It should be noted that the above analysis is based on a simplified setting. In practice, we need to consider the maximum sailing speed Vmax and the weekly service frequency. Nevertheless, the analytical analysis provides useful managerial insights into the problem. 5.5. Analytical analysis of the consequence of wait time variability ~ at To obtain insights into the consequence of wait time variability, we focus on a single port. Assume that the wait time w port i is uniformly distributed over [Wmin, Wmax]. Represent by t :¼ ti  ui the buffer time at port i. Denote by C the penalty per unit time of delay. We minimize the sum of the expected idle cost at port plus late handling penalty:

minE½tcopr =168jw 6 t þ E½wcopr =168jw > t þ E½Cðw  tÞjw > t t

ð77Þ

where the first item is the expected ship cost when the ship is in idle at port given w 6 t; the second item is the expected ship cost when the ship is in idle at port given w > t, and the third item is the expected penalty cost for late handling given w > t. We can analytically derive that the optimal buffer time t⁄ satisfies:

W max  t  

t W

min

¼

copr =168 C

ð78Þ

Eq. (78) can be intuitively interpreted as follows: the optimal buffer time t⁄ is designed such that the ship idle cost as a result ~ takes the minimum value Wmin equals the late handling penalty when w ~ takes the maximum value of early arrival when w ^ :¼ ðW max þ W min Þ=2 and r ^ whether t⁄ will increase or decrease with r ^ :¼ ðW max  W min Þ=2. Given W, ^ deWmax. Define W ^ and t⁄ decreases with r ^ . If copr/168 < C, t⁄ increases pends on the relative value of copr/168 and C. If copr/168 > C, then t < W ^. with r

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

1105

6. Conclusions This paper has examined the design of ship route schedules that hedge against the uncertainties in port operations, which include the uncertain wait time due to port congestion and the uncertain container handling time. The designed schedule is robust in that the uncertainties in port operations and schedule recovery by fast steaming are captured endogenously. A mixed-integer nonlinear stochastic programming model is developed. To solve this model, a solution algorithm is proposed, which incorporates a sample average approximation method, linearization techniques, and a decomposition scheme. Extensive numerical experiments based on an Asia–America–Europe ship route demonstrate that the algorithm obtains near-optimal solutions with the optimality gap less 1.5% within 1 h of CPU time. Considering that the robust schedule design problem is a tactical-level planning issue, the proposed solution algorithm is effective and efficient for solving practical problems. Numerical experiments also demonstrate that liner shipping companies can improve their schedule robustness by adding more ships to their ship routes while incurring more costs. Therefore, the proposed model and algorithm in this paper provide a quantitative analysis tool for liner shipping companies to obtain the optimal trade-off between schedule robustness and total costs.

Acknowledgements We would like to thank an anonymous reviewer for her/his thorough and helpful suggestions. This study is supported by the research Grants – WBS No. R-264-000-244-720 and WBS No. R-302-000-014-720 – from the NOL Fellowship Programme of Singapore.

References Alvarez, J.F., 2009. Joint routing and deployment of a fleet of container vessels. Maritime Economics & Logistics 11, 186–208. Bell, M.G.H., Bichou, K., 2008. An analysis into speed and schedule stability for container liner services. Maritime Economics and Logistics 10, 175–184. Bell, M.G.H., Liu, X., Angeloudis, P., Fonzone, A., Hosseinloo, S.H., 2011. A frequency-based maritime container assignment model. Transportation Research 45B, 1152–1161. Bierwirth, C., Meisel, F., 2010. A survey of berth allocation and quay crane scheduling problems in container terminals. European Journal of Operational Research 202, 615–627. Brouer, B.D., Pisinger, D., Spoorendonk, S., 2011. Liner shipping cargo allocation with repositioning of empty containers. INFOR 49 (2), 109–124. Christiansen, M., Fagerholt, K., Ronen, D., 2004. Ship routing and scheduling: status and perspectives. Transportation Science 38, 1–18. Christiansen, M., Fagerholt, K., Nygreen, B., Ronen, D., 2007. Maritime transportation. In: Barnhart, C., Laporte, G. (Eds.), Handbook in OR & MS, vol. 14. Elsevier, pp. 189–284. Dong, J.X., Song, D.P., 2009. Container fleet sizing and empty repositioning in liner shipping systems. Transportation Research 45E, 860–877. Dong, J.X., Song, D.P., 2012. Quantifying the impact of inland transport time on container fleet sizing in a liner shipping service with uncertainties. OR Spectrum 34, 155–180. Drewry Shipping Consultants, 2006. The Drewry Container Shipper Insight-Further Quarter 2006. Drewry Shipping Consultants, London. Fagerholt, K., 1999. Optimal fleet design in a ship routing problem. International Transactions in Operational Research 6, 453–464. Fagerholt, K., 2004. Designing optimal routes in a liner shipping problem. Maritime Policy and Management 31 (4), 259–268. Gelareh, S., Pisinger, D., 2011. Fleet deployment, network design and hub location of liner shipping companies. Transportation Research 47E, 947–964. Jepsen, M.K., Løfstedt, B., Plum, C.E.M., Pisinger, D., Sigurd, M.M., 2011. A path based model for a green liner shipping network design problem. In: Proceedings of the International MultiConference of Engineers and Computer Scientists, 2011. Maersk, 2011. Maersk Line Sailing Schedules. (accessed 18.04.11). Mak, W.K., Morton, D.P., Wood, R.K., 1999. Monte Carlo bounding techniques for determining solution quality in stochastic programs. Operations Research Letters 24, 47–56. Meng, Q., Wang, T., 2011a. A scenario-based dynamic programming model for multi-period liner ship fleet planning. Transportation Research 47E, 401–413. Meng, Q., Wang, S., 2011b. Liner shipping service network design with empty container repositioning. Transportation Research 47E, 695–708. Notteboom, T.E., 2006. The time factor in liner shipping services. Maritime Economics and Logistics 8, 19–39. OOCL, 2011. Service Routes.
1106

S. Wang, Q. Meng / Transportation Research Part E 48 (2012) 1093–1106

Verweij, B., Ahmed, S., Kleywegt, A.J., Nemhauser, G., Shapiro, A., 2003. The sample average approximation method applied to stochastic routing problems: a computational study. Computational Optimization and Applications 24, 289–333. Vis, I.F.A., de Koster, R., 2003. Transshipment of containers at a container terminal: an overview. European Journal of Operational Research 147, 1–16. Wang, S., Meng, Q., 2011. Schedule design and container routing in liner shipping. Transportation Research Record 2222, 25–33. Wang, S., Meng, Q., 2012a. Liner ship route schedule design with sea contingency time and port time uncertainty. Transportation Research 46B, 615–633. Wang, S., Meng, Q., 2012b. Sailing speed optimization for container ships in a liner shipping network. Transportation Research 48E, 701–714. Yan, S., Chen, C.-Y., Lin, S.-C., 2009. Ship scheduling and container shipment planning for liners in short-term operations. Journal of Marine Science and Technology 14 (4), 417–435. Zhen, L., Lee, L.H., Chew, E.P., 2011a. A decision model for berth allocation under uncertainty. European Journal of Operational Research 212, 54–68. Zhen, L., Chew, E.P., Lee, L.H., 2011b. An integrated model for berth template and yard template planning in transshipment hubs. Transportation Science 45 (4), 483–504.