Weather impact on containership routing in closed seas: A chance-constraint optimization approach

Transportation Research Part C 55 (2015) 139–155

Contents lists available at ScienceDirect

Transportation Research Part C journal homepage: www.elsevier.com/locate/trc

Weather impact on containership routing in closed seas: A chance-constraint optimization approach Konstantinos Kepaptsoglou a,⇑, Grigorios Fountas b, Matthew G. Karlaftis b,1 a b

School of Rural and Surveying Engineering, National Technical University of Athens, 9, Iroon Polytechniou str, 15770 Athens, Greece School of Civil Engineering, National Technical University of Athens, 5, Iroon Polytechniou str, 15773 Athens, Greece

a r t i c l e

i n f o

Article history: Received 29 October 2014 Received in revised form 25 January 2015 Accepted 28 January 2015 Available online 11 March 2015 Keywords: Ship routing Containerships VRP Stochastic travel times Chance-constrained model Pick-ups and deliveries Time deadlines

a b s t r a c t Weather conditions have a strong effect on the operation of vessels and unavoidably influence total time at sea and associated transportation costs. The velocity and direction of the wind in particular may considerably affect travel speed of vessels and therefore the reliability of scheduled maritime services. This paper considers weather effects in containership routing; a stochastic model is developed for determining optimal routes for a homogeneous fleet performing pick-ups and deliveries of containers between a hub and several spoke ports, while incorporating travel time uncertainties attributed to the weather. The problem is originally formulated as a chance-constrained variant of the vehicle routing problem with simultaneous pick-ups and deliveries and time constraints and solved using a genetic algorithm. The model is implemented to a network of island ports of the Aegean Sea. Results on the application of algorithm reveal that a small fleet is sufficient enough to serve network’s islands, under the influence of minor delays. A sensitivity analysis based on alternative scenarios in the problem’s parameters, leads to encouraging conclusions with respect to the efficiency and robustness of the algorithm. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Reliability is an important challenge in shipping operations, since it guides shippers/consignees into arranging accordingly their production, transportation and distribution plans (Wang and Meng, 2012a). According to Vernimmen et al. (2007), vessels delays can be caused by factors beyond the control of shipping firms, such as weather conditions, port congestion and labor strikes. The adverse effect of weather conditions on the operation of commercial vessels has a significant influence on the total time at sea, affecting schedule reliability and increasing operational costs. In particular, the wind affects the performance of vessels by increasing both travel and service time at port and consequently the overall cost of transportation. The Aegean Sea is a case where weather and especially, wind and waves play a critical role for maritime activities (Mazarakis et al., 2012). In terms of freight shipping operations, maritime transportation between the islands of the Aegean Sea and the Greek mainland is classified as short sea shipping, where liner companies seek to minimize transportation time of cargo.

⇑ Corresponding author. Tel.: +30 2107723167; fax: +30 210 7722404. 1

E-mail addresses: [email protected] (K. Kepaptsoglou), [email protected] (G. Fountas), [email protected] (M.G. Karlaftis). Tel.: +30 2107721280; fax: +30 210 7722404.

http://dx.doi.org/10.1016/j.trc.2015.01.027 0968-090X/Ó 2015 Elsevier Ltd. All rights reserved.

140

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

In order to establish a stable freight transport network, disengaged from the problems of Ro–Ro passenger liners, Sambracos et al. (2004) introduced the potential use of small containers for carrying freight to a network of Aegean islands and considered the associated vehicle routing problem (VRP) for planning a small containership service. Later, Karlaftis et al. (2009) investigated routing for a fleet of small containerships between hub-and-spoke ports of the Aegean Sea by additionally considering simultaneous pick-ups and deliveries and time deadlines. In this paper, we extend previous work by Karlaftis et al. (2009) on containership routing by introducing weather effects on travel times between ports. Indeed, uncertainty in weather conditions contributes to the stochastic nature of vessel’s travel time and therefore incorporation of its impact in the routing process allows more realistic outcome. Nhe consideration of the stochastic nature of travel times variable (a) offers a better representation of the real problem conditions, (b) allows a better evaluation of whether vessels arrive within preset time deadlines or not (6) and (c) aids in a better estimation of possible delays and additional transportation cost for maritime companies. For this purpose, we adopt a chance-constrained programming model, in order to formulate the problem and to convert its stochastic nature to a deterministic form (Zhang et al., 2012). We also develop and implement a genetic algorithm, adopted to the requirements of the stochastic travel time vehicle routing problem with simultaneous pick-ups and deliveries and time deadlines (STT-VRPSPDTD). The proposed approach focuses on important aspects on the strategic design of a dedicated containership network in the Aegean Sea, and particularly the introduction of weather impacts on establishing service efficient vessel routes. Nhe paper is organized as follows. Section 2 discusses the background which deals with STT-VRPSPDTD and wind distributions in Aegean Sea. In Section 3, the problem is described in detail and the mathematical model is presented, while in Section 4 the solution process is described in detail, including the genetic algorithm details and the method used for the estimation of stochastic travel times. The computational results of the application to the Aegean Sea are reported in Section 5. Section 6 concludes our paper and discusses suggestion for future research.

2. Related work In the maritime industry, three modes of shipping operations are mentioned: industrial, tramp and liner shipping (Lawrence, 1972). According to Sigurd et al. (2005), the transportation of containers via sea ways is considered as a part of liner shipping. Liner vessels follow a fixed route, according to a published itinerary (Christiansen et al., 2013). A number of reviews concerning research in maritime transportation have been published in the last 30 years (Ronen, 1983, 1993; Christiansen et al., 2004; Christiansen et al., 2007, 2013). In addition, Kjeldsen (2011) has investigated routing and scheduling of ships and cargo in liner shipping. Meng et al. (2013) prepared a review paper on containership routing and scheduling problems and Tran and Haasis (2013) surveyed work on container shipping and indicated three major planning problems on the topic: container routing, fleet management and network design. A wide range of papers on liner shipping has been published in the last decade, mainly due to the increase in cargo containerization (Christiansen et al., 2013). According to Christiansen et al. (2013) the research on liner shipping operations can be categorized in four groups, depending on the type of the network design model exploited (see Table 1). The first group focuses on models with a single route or set of routes without transshipment. Among recent studies in this category, Sambracos et al. (2004) formulated a capacitated vehicle routing problem considering a network of 13 island of the Aegean Sea, in order to eliminate the dependence of freight transport from the passenger shipping. Shintani et al. (2007) investigated the design of containership routes taking into account the repositioning of empty containers and developed a genetic algorithm for this purpose. Karlaftis et al. (2009) extended the work of Sambracos et al. (2004) and proposed a genetic algorithm approach for routing of a containership service considering simultaneous pickups and deliveries at ports and time deadlines. Meng and Wang (2011a) developed a dynamic programming model in order to solve efficiently a multiperiod liner ship fleet programming problem. They created a scenario decision tree for determining the optimal fleet deployment plan. Chuang et al. (2010) constructed a fuzzy genetic algorithm in order to determine optimal containerships routes taking into account demand, travel and berthing time. The network models, in which all feeder ports are connected to a hub port, belong to the second group. Meng and Wang (2011b) considered an intermodal hub and spoke network and formulated a mathematical program with equilibrium constraints (MPEC), proposing a hybrid genetic algorithm for its solution. The third category concentrates on models where some ports of the network are considered as hub ports, without any restriction on the number of hub and non-hub ports a vessel can visit. In this context, Reinhardt and Pisinger (2012) investigated a network design problem in combination with a fleet assignment problem. The formulation of the problem allowed simple and butterfly routes with transshipment in the central port and a branch-and-cut method was adopted for the solution of the model. The last group consists of networks where the ports are neither classified as hub and non-hub. Imai et al. (2009) formulated one multi-port calling (MPC) network and one hub and spoke (H&S) network and examined the container management cost including empty container repositioning. For the first network they developed a genetic algorithm and for the H&S network a brute-force heuristic method was used. The joint optimization of vessel routing and deployment of fleet was investigated by Álvarez (2009). He formulated a mixed integer programming model, which considers multiple types of vessels reflecting differences in terms

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

141

Table 1 Recent papers on liner shipping. Paper

Model category

Scope

Objective function

Formulation and solution method

Sambracos et al. (2004) Shintani et al. (2007)

1st group: Single route or sets of routes without transshipment 1st group: Single route or sets of routes without transshipment 2nd group: Hub and feeder routes 1st group: Single route or sets of routes without transshipment 4th group: Ports are not classified as hub and non-hub

Fleet size

Minimum costs

Linear Programming (LP) + Mixed Integer Programming (MIP) & metaheuristic

Routing

Maximum profit

Heuristic

Routing

Minimum costs Minimum costs

Metaheuristic

Minimum costs

Heuristic + Brute-force method

Maximum profit

MIP + Heuristics + Benders Decomposition

Ship routing

Maximum profit

Fuzzy genetic algorithm

Routing and fleet deployment

Minimum costs Maximum profit

MIP + metaheuristics

Karlaftis et al. (2009) Sigurd et al. (2005) Imai et al. (2009) Agarwal and Ergun (2008) Chuang et al. (2010) Álvarez (2009) Meng and Wang (2011a) Meng and Wang (2011b) Reinhardt and Pisinger (2012)

4th group: Ports are not classified as hub and non-hub 1st group: Single route or sets of routes without transshipment 4th group: Ports are not classified as hub and non-hub 1st group: Single route or sets of routes without transshipment 2nd: Hub and feeder routes

3rd: Some ports of the network are hub ports

Fleet composition, visit-schedule and ship routing Comparison of network models taking into account container management Ship scheduling and cargo routing

Multi-period liner

Optimal hub-and-spoke network design

Minimum costs

Network design and fleet assignment

Minimum costs

DW decomposition + heuristic method (branch and price)

Dynamic programming (DP)

Mathematical program with equilibrium constraints (MPEC) + hybrid genetic algorithm MIP + branch and cut algorithm

of cost and operational characteristics and allows the transshipment of cargo. By fixing the values of some variants, the problem was transformed to a multi commodity flow problem, which was solved by an interior point algorithm. In existing research devoted to liner ship route design and scheduling, travel time of ships (time at sea) is usually considered as a deterministic variable. Wang and Meng (2012b) take into account uncertainty at sea, calculating sea contingency time for each voyage leg, and uncertainty at port, formulating the available voyage time as a random variable. In this work, we employ a network design model considering uncertainty at sea, by estimating the impact of weather conditions on sailing speed and consequently on sailing time. Therefore, travel times between the island ports are assumed to be stochastic variables, in order to maintain the reliability of the container service schedules.

3. Problem description and model formulation In the context of designing a freight transport system in the Aegean Sea, Sambracos et al. (2004) formulated a capacitated vehicle routing problem (CVRP) considering a homogeneous fleet of containerships which serves a network of 13 islands using small containerships. These island ports cannot be supplied and are characterized only by inbound cargo. In real-life conditions, however, there are periods during the year when the supply of the island ports is quantitatively comparable to demand, a fact demonstrates the need for a seasonal redesign of routes (Sambracos, 2000). In this direction, Karlaftis et al. (2009), extended the initial problem and considered additional parameters such as pick-ups and deliveries. Indeed, larger islands periodically have a significant export activity, fed by local production. Due to sensitivity of transported goods (products of primary and secondary sector of the economy), the paper by Karlaftis et al. (2009) included time deadlines in order to reduce the delays of the deliveries and to limit the possibility of product obsolescence. In this case, the islands are characterized in the network, not only by demand but also by supply at a time period. A fundamental assumption of these formulations is the deterministic nature of travel times, i.e. the times between the network nodes (islands and central port) are considered fixed and only depend on the corresponding distances and the predetermined speed of the containership. But, in real conditions there are several factors which can significantly affect the

142

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

travel or service time, making them uncertain variables (Wang and Meng, 2012b). This uncertainty may completely alter routing conditions, causing successive delays to arrivals and ultimately a significant increase of transportation costs. These factors can be unpredictable, (such as mechanical damage to the ship or a random event at the port facility), or expected (traffic jam at berthing positions, condition of the loading area). Weather conditions also contribute to that uncertainty and constitute a common reason for vessels delays (Vernimmen et al., 2007). More specifically, wind and waves are the major factors affecting ship performance and provoking speed loss (Kwon, 1981). Therefore, the consideration of stochastic travel times is more appropriate, on the basis of the above assumptions. The background of previous research and the proposed extensions of the problem result in the formulation of a stochastic vehicle routing problem of a homogeneous containership fleet with the following constraints satisfied:  Each vehicle route begins from and terminates to a central depot. In this case, the central mainland port plays the role of depot.  The containerships not only deliver goods to the island ports, but also pick-up some products.  The equilibrium supply-demand at each node of the network cannot exceed the vessel’s capacity.  Each node of the network is served only once, by containerships with predetermined capacity.  Time deadlines are adopted, since the arrival at each destination port cannot exceed a predetermined period after the departure.  Service times are modeled as deterministic variables. The objective of this STT-VRPSPDTD consists of scheduling a containership service, to minimize the required number of routes as well as to minimize the corresponding total transportation cost. The stochastic VRP can be cast within the framework of stochastic programming (Gendreau et al., 1996). Such a stochastic problem can be modeled either as a chance constrained programming model (CCP) or a stochastic programming model with recourse (SPR) (Gendreau et al., 1996). In this paper, we develop a chance-constrained programming model (CCP), by extending the formulation proposed by Karlaftis et al. for the vehicle routing problem with pick-ups and deliveries. In particular, the proposed model aims to the determination of routes which correspond to minimal transportation costs under stochastic chance constraints, with respect to weather affected travel times. The STT-VRPSPDTD can be defined as follows. Let G = (Vo, A) be a directed graph, where Vo = (0, 1, 2, . . ., n) is the set of nodes representing the island ports and A is the set of arcs. A fleet of containerships is based at a mainland port (depot), which is labeled as node 0. Each node corresponds to a known demand and supply and to a fixed service time. The travel time on each arc in A is a continuous random variable, which follows a probability distribution. Also, travel time and service time variables are assumed to be independent. Before providing the formulation of the problem, we first define the notations used. Let: C: Set of Nodes (Ports) W: Set of Containerships i, j, m: Nodes 2 C k: Containership 2 W di: Demand for node i pj: Supply for node j cij: Travel cost for arc (i, j) Qk: Load capacity for containership k tijk: Travel time between nodes i and j for containership k sik: Service time for node i by containership k Tk: Maximum route travel time for containership k ai, aj: Arrival times at ports i, j ltj: Latest arrival time at port j l0k: Load for ship k when leaving the mainland port ljk: Load for ship k after leaving port j A: Probability value M: Arbitrarily large number

xijk

8 > < 1 if ship k uses arcði; jÞ ¼ 0 otherwise > :

The mathematical formulation of the problem is the following:

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

minimize

XXX

ðcij  xijk Þ

143

ð1Þ

i2C j2C k2W

Subject to XX

xijk ¼ 1 8j 2 C

ð2Þ

xijk ¼ 1 8i 2 C

ð3Þ

i2C k2W

XX j2C k2W

X X ximk ¼ xmjk i2C

j2C

X X di xijk i2C

8m 2 C; k 2 W

ð4Þ

8k 2 W

ð5Þ

! 6 Qk

j2C

X X XX sik xijk þ tijk  xijk 6 T k i2C

P X

j2C

i2C j2C

X X XX sik xijk þ tijk  xijk i2C

j2C

8k 2 W

!

ð6Þ

! 6 Tk

P A 8k 2 W

ð7Þ

i2C j2C

x0jk 6 1 8k 2 W

ð8Þ

j2C=f0g

X

xi0k 6 1 8k 2 W

ð9Þ

i2C=f0g

xijk 2 S 8i 2 C; j 2 C; k 2 W ( ) XX xijk 6 jBj  1 for B # C=f0g; jBj P 2 S ¼ xijk :

ð10Þ ð11Þ

i2B j2B

aj P ai þ sik þ tijk  ð1  xijk Þ  T k aj P ai þ sik þ tijk þ ð1  xijk Þ  T k

8i; j 2 C; k 2 W 8i; j 2 C; k 2 W

ð12Þ ð13Þ

a0 ¼ 0

ð14Þ

aj 6 lt j 8j 2 C XX dj  xijk l0k ¼

ð15Þ

8k 2 W

ð16Þ

i2C j2C

ljk P l0k  dj þ pj  M  ð1  x0jk Þ 8j 2 C; k 2 W ! X ljk P lik  dj þ pj  M  1  xijk 8i; j 2 C; i – j

ð17Þ ð18Þ

k2W

l0k 6 Q k 8k 2 W ljk 6 Q k 8j 2 C; k 2 W

ð19Þ ð20Þ

Objective function (1) minimizes the total travel cost of all routes. The set of constraints (2)–(11), except constraint (7), are part of the typical formulation for a capacitated VRP (Karlaftis et al., 2009). Constraints (2) and (3) guarantee that each port is visited only once by ship, while constraint (4) indicates the conservation of flow at each demand node. According to constraint (5), cargo cannot exceed the ship capacity, while in constraint (6) we bound the maximum route duration. The constraints (8) and (9) ensure that only the available number of ships can be used, while constraints in Eqs. (10) and (11) eliminate sub-tours. Constraint (7) is the chance-constrained condition, which assures that total travel time of a ship cannot exceed the maximum duration of route Tk, for a probability A. Probability A is the confidence level provided as an appropriate safety margin by the decision-maker (Li et al., 2010). The definition of this level is associated with the probability of excess of time deadline. If we consider a low value for A, this means that a route can exceed the maximum duration and the appearance of significant delays is probable. This chance-constrained condition permits the transformation of dynamic problem to a deterministic problem (Zhang et al., 2012). In this context, a feasible solution can be found in a deterministic equivalent. The travel time on each arc (i, j) follows a probability distribution, which is not known a priori. It is assumed that the probability distribution of each arc is independent and is not affected by other links’ travel time distributions. This fact constitutes a fundamental assumption of the proposed chance-constrained model. Constraints (12)–(15) are derived from the consideration of time deadlines (Wiley, 2000). The first two constraints ensure the successiveness of arrivals in ports, while constraint (15) indicates that arrival times are limited by latest possible times. Finally, constraints (16)–(20) are additional conditions for simultaneous pick-ups and deliveries. These constraints are proposed by Dethloff (2001) and are adopted by Karlaftis et al. (2009). Constraint (16) defines the ship load when leaving the central port (depot), while Eqs. (17) and (18) determine the ship load when arriving at the first and the next ports of the route. The last two constraints are equivalent with condition (5) and are related to the ship capacity.

144

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

4. Solution method 4.1. Overview As noted earlier, we treat the stochastic vehicle routing problem as a chance-constrained formulation; in that sense the probabilistic nature of travel times is transformed into a deterministic constraint depicted in Eq. (7). The problem is then tackled as follows:  A pre-processing step described later on is developed for estimating expected (deterministic) travel times per segment for a fixed confidence level a; this step is used for tackling constraint (7).  An optimization algorithm is then used for solving the associated deterministic VRP for these travel times. It should be noted that due to the assumed deterministic nature of service times, the confidence level affects only travel times when transforming constraint (7) into a deterministic equation. In this case, we select an appropriate confidence level and use it to estimate travel times accordingly. These are then fed as deterministic travel times in the optimization process. 4.2. Solution techniques for the VRP The literature exhibits a variety of heuristic and metaheuristics algorithms focusing on efficiently solving the VRP and its variants. Traditional heuristics were mostly developed between 1960 and late 1990s (Mole and Jameson, 1976; Wark and Holt, 1994; Laporte et al., 2000), while metaheuristics, which combine advances in computing power and scientific evolution of algorithms, have been developed and implemented over the last 20 years. Gendreau et al. (1999) refer to a variety of metaheurisics techniques for the VRP and its extensions: Ant colony optimization, Genetic algorithms, Greedy Randomized Adaptive Search Procedure, Simulated Annealing, Tabu Search, Variable neighborhood search. Tabu Search is often used for the VRP and many researchers have proposed algorithms based on this method (Gendreau et al., 2002, 1994; Rego, 2001; Taillard, 1993; Cordeau et al., 2001). Li et al. (2010) proposed a tabu search based heuristic for solving a stochastic vehicle routing problem with time windows. Application of genetic algorithms have been reported not only for the classic VRP but also for some of its variants (Blanton and Wainwright, 1993; Thangiah, 1995; Bräysy and Gendreau, 2001; Baker and Ayechew, 2003; Jih and Hsu, 2004; Prins, 2004; Karlaftis et al., 2009). Also, scatter search approaches have been reported in the last few years for some types of VRP (Russell and Chiang, 2006; Zhang et al., 2012). Zhang et al. (2012) in particular implemented a scatter search method to address the combinatorial nature of a stochastic VRP with simultaneous pick-ups and deliveries. In this paper, we use a genetic algorithm for solving the STT-VRPSPDTD. The wide application of this technique for solving routing problems with time windows (Bräysy et al., 2004), its stochastic nature (Li et al., 2010) and the need for a solution of good quality in an acceptable computational time (Jih and Hsu, 2004; Karlaftis et al., 2009) are major drivers toward selecting Gasfor the problem at hand. 4.3. Genetic algorithm Genetic algorithms (GA), first introduced by Holland (1975) are a widespread metaheuristic, which can be easily applied to complex optimization problems. In this paper, we adopt a steady state genetic algorithm (Luke, 2013). Capacity constraints are considered as hard constraints, while time deadlines can be violated (soft constraints) and delays are penalized. Each chromosome is represented by the sequence of nodes a containership must visit. The length of the string coincides with the size of the network. For example, if the network consists of 9 nodes, a possible chromosome could be: Considering the depot as node 0, we can identify a number of possible trips within this chromosome, given the constraints of the algorithm. For instance, for the case of Fig. 1, a set of three candidate trips would be 0-3-5-28-0, 0-1-6-4-0, 0-7-9-0. The fitness measure is assumed to be the total travel time (travel and service time) for each possible set of routes. A penalty is also added for time delays and is considered as a percentage of time delay. The fitness function is given by:

FF ¼

X X r2R n2V r # V

t r;ðn1Þn þ

X X sm þ pen  drm m2V

ð21Þ

r2R;m2V

where pen is the penalty percentage and drm is the time delay in reaching port m for route r. The selection of the initial population is derived from a random procedure. The individuals are selected as parents according to their fitness ranking (Eiben and Smith, 2003). A steady state approach is used for parent selection and replacement. Steady state GAs exhibit improved computational performance over generational ones: at each generation only two parents are selected for recombination based on their ranking. The best two individuals of these two parents and their offspring, are used for replacing parent chromosomes and thus the population size remains constant (Luke, 2013). During the process of crossover (Davis, 1991; Eiben and Smith, 2003) two crossover points are generated randomly at the first parent and the genes between the points are copied to the first offspring. The other genes of the first parent are copied in the first offspring with the appearance sequence of the second parent, starting from the second crossover point. The same process is applied for

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

145

Sequence=voyage order of a containership if only one could serve the network

3

5

2

8

1

6

4

7

9

Fig. 1. Example of chromosome representation.

the second offspring. Furthermore, the algorithm proceeds to the procedure of mutation, where two random genes in each chromosome are reversed in order to maintain the diversity of individuals (see Fig. 2). Finally, a route extraction algorithm is used, which extracts possible routes for each sequence based on vessel capacity and calculates total travel and service cost at each port of the route; the interested reader is referred to Karlaftis et al. (2009) for more details on the extraction procedure. The overall approach (including the pre-process step) is shown in Fig. 3. 4.4. Estimation of stochastic travel times 4.4.1. Overview While in the deterministic problem, the vessel’s speed is assumed constant and can be defined a priori, the same cannot be considered in stochastic routing, where the speed of the ship is dependent on weather conditions. In such cases, the vessel’s speed follows a probability density distribution and the travel time between two nodes varies accordingly. Travel time distribution is originally not known since it depends on the complex mechanism of interaction between the wind and wave action and the vessel performance and particularly on the spatially varying wind characteristics (Kwon, 1981). In this work, we apply an approximate process for deriving travel time distributions, which is summarized in the following Fig. 4. 4.4.2. Estimation of wind speed probability distribution For the purpose of assessing the impact of weather conditions on the vessel’s performance and thus on travel time, we have to identify a set of wind parameters, such as average wind velocity and wind direction. The combined action of wind and waves is a major factor affecting operations and performance of a vessel (Kwon, 1981). The wind resistance is a function of the square of wind speed, so an increase in wind velocity makes cruising significantly more difficult (Kwon, 1981). Waves are also generated with the emergence of wind, which are developed in line with increasing ship pressure and affect ship performance. When taking into consideration wind conditions in routing, simple knowledge of the average velocity of the wind at the region is not sufficient. Long-term detailed wind speed data is necessary to extract the probability distribution of wind velocity at a time period. For the case at hand, the extent of data and the complexity of subsequent statistical processing, in combination with the lack of accurate measurements in many regions of the Aegean Sea, lead to the adoption of generalized functions which describe the wind profile of a region, considering a small number of parameters. The probability distributions most broadly used are those of ‘‘Weibull’’ and ‘‘Rayleigh’’ (Kaldelis, 1999). The analytical expression of ‘‘Weibull’’ distribution is given by the equation:

(  )  k1 k k V V f ðVÞ ¼  exp  C C C

ð22Þ

where f(V) is the probability density function of wind speed, V is wind speed and k, C are Weibull parameters. The Weibull distribution is considered as rather sufficient for an altitude of up to 100 m, using two main parameters, namely the scale ‘‘C’’ and the shape ‘‘k’’ factors so as to determine the probability density (Zafirakis et al., 2012). In the case of the Aegean Sea, researchers have created maps, which describe the spatial distribution of Weibull parameters (k, C). Figs. 5 and 6 present the curves with the values of parameters (Kaldelis, 1999). The Wind and Waves Atlas of the North-Eastern Mediterranean Sea (1992) provides a map with 33 wave homogeneous regions covering the main part of the eastern Mediterranean Sea (Fig. 7). These regions have been derived by statistically testing wave climate homogeneity and each region is characterized by stability with respect to its main wind and waves parameters. Wind statistics of these zones are published in the corresponding Atlas part.

(3,5,2,8,1,6,4,7,9) (4,9,1,5,2,6,8,3,7)

(4,9,5,8,1,6,2,3,7) (3,8,1,5,2,6,4,7,9)

3

5

2

8

1

6

4

9

7

3

5

2

6

1

8

4

9

7

Fig. 2. Examples of crossover and mutation operators.

146

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

Fig. 3. Process flowchart.

Estimation of wind speed probability distribution

Estimation of vessel speed probability distribution

Extraction of travel time probability distribution

Fig. 4. Process of determining travel times distributions.

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

147

Fig. 5. Geographical distribution of parameter ‘‘k’’ (Kaldelis, 1999).

Fig. 6. Geographical distribution of parameter ‘‘C’’ (Kaldelis, 1999).

Based on the above assumptions and using data from relevant research efforts (Bagiorgas et al., 2012; Zafirakis et al., 2012; Kaldelis, 1999), we apply a grid on the Aegean sea, with geographically defined ‘‘homogeneous wind areas’’, based on associated values of Weibull parameters for the wind distribution. Thus, the Aegean Sea is divided in 45 equally sized squares (52 nautical miles each part side), and for each part Weibull parameters are estimated (Fig. 8). These values are fixed and reflect the probability distributions of wind velocity per square of the grid. 4.4.3. Estimation of vessel speed probability distribution During a sea trip, a vessel’s speed may be reduced. This reduction may be voluntary (voluntary speed loss) or unintentional (involuntary speed loss). The weather causes of these speed losses may broadly be described as the effects of wind and waves (Kwon, 1981). Involuntary speed reduction results from increased resistance in a seaway (wave added resistance, wind resistance) while voluntary speed reduction is the deliberate reduction in speed by the ship’s captain in order to ensure that the ship’s wave induced responses are within acceptable safety limits (Chen and Padhy, 2010). The determination of

148

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

Fig. 7. Wind and wave homogeneous regions (Wind and Waves Atlas of the North-eastern Mediterranean Sea, 1992).

Fig. 8. Grid of wind homogeneous zones in Aegean Sea.

ship speed loss, due to weather conditions is not an easy procedure. In light winds (less than 20-knots), ships lose speed in head winds and gain speed slightly in tail winds. For higher wind speeds, ship speed is reduced in both head and tail winds. This is due to the increased wave action, which even under tail winds, results in increased drag from steering corrections, and indicates the importance of sea conditions in determining ship performance (Bowditch, 2002). In addition, the loading status of the vessel (type, surface and volume of the cargo) is an important factor for its dynamic behavior. For example, strong winds will have a greater negative effect on a large, fully loaded containership than a fully loaded tanker of similar length (Bowditch, 2002). In order to estimate speed loss due to combined wind and wave impact, we adopt an approximate formula, which was developed for practical use by Townsin and Kwon (Kwon, 1981). The formula was developed particularly for

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

149

closed-seas (as the case of the Aegean), where waves created due to wind in distant areas do have practically any effect to routes (as islands intervene between route segments). A percentage of speed loss for the case of containerships is given by:

DV BN 6:5 ¼ 0:7  BN þ 2 V 22  r3

ð23Þ

where DV/V is the speed loss due to head wind and waves, BN the wind force given by Beaufort Scale, r the volume of displacement (m3). This formula leads to speed loss estimation due to weather impact, considering only wind data and without introducing wave parameters, on the basis that the reference vessel operates under conditions of short sea waves (such as in the case of the Aegean sea). Given this approximate formula, we can derive the vessel’s performance curve, which is used to estimate the ship speed, depending on prevailing weather conditions. For a typical, small containership, the speed curve has the following form (see Fig. 9). It is noted that the above formula is used for the determination of involuntary speed loss, so for higher Beaufort numbers (BN > 7) the speed loss cannot computed by the approximate relation (Kwon, 1981). Above this wind level, involuntary speed loss is less significant than voluntary (Kwon, 1981). 4.4.4. Estimation of travel times probability distributions According to the aforementioned assumptions, an analytical procedure is followed for extracting travel time probability distributions between ports. The steps of this procedure are described below: Step 1: Geographical survey of routes For each pair of successive ports, we identify the path followed in real conditions by passenger or cargo vessels (paths are readily available in nautical charts). We then determine grid squares shown in Fig. 8 and estimate a contribution um of each square along a path. This contribution of each square is assumed to be the proportion of the length of the path’s segment within that square over the total path length. Step 2: Determination of wind speed distribution for each grid part of the route Each grid square along the path is related to some predetermined values of Weibull parameters. This means that for each segment of the path contained in a grid square, the wind velocity probability distribution is known. Therefore, wind velocity distributions can be derived along all segments of a path: for a path segment m included in a grid square, its wind velocity probability distribution fm(V) can be determined. Given that probability distribution is derived from differentiated Weibull parameters, which were spatially estimated for each grid part, the distribution of each path segment is independent and unaffected from the wind conditions of the other parts. The similarities and the interactions of the adjacent sea areas concerning prevailing wind conditions have been incorporated in the process of parameters determination. Step 3: Determination of total wind speed distribution for the route between ports The probability ftot(V) of a specific wind velocity to appear is estimated as the weighted sum of individual probabilities fm(V) of each path segment m, using the corresponding contribution values um as weights (estimated in Step 1): 0

f tot ðVÞ ¼

X

um  f m

ð24Þ

m

Step 4: Estimation of speed loss and determination of ship speed probability distribution

SHIP SPEED (KNOTS)

14.00 12.00 10.00 8.00 6.00 4.00 2.00 0.00 0.00

2.00

4.00

6.00

8.00

10.00

WIND FORCE (BEAUFORT) Fig. 9. Speed curve for containership similar to reference vessel.

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

PROBABILITY DENSITY f (V)

150

0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0.00

5.00

10.00

15.00

20.00

25.00

30.00

35.00

TRAVEL TIME (HOURS) Fig. 10. Travel time probability distribution.

Speed loss is calculated using Eq. (23). Assuming a normal value of vessel speed (when cruising under minimal wind effects) and taking into account the wind velocity probability distribution, we can straightforwardly determine the vessel’s speed probability distribution. Step 5: Travel time probability distribution As the distance between ports is known and fixed, travel times and their corresponding probability distributions can be calculated for different vessel speeds and consequently for their corresponding wind velocities. An indicative travel time probability distribution for a pair of ports is illustrated at the following graph (Fig. 10). The value-starting point of the distribution coincides with the deterministic value of travel time and as it shown, the majority of probable times is observed near this value. 5. Application 5.1. Overview The proposed algorithm is applied for routing a homogeneous fleet of small containerships from the mainland port of Piraeus to a number of islands in the Aegean Sea. The investigated network consists of 26 island ports, where each island port corresponding to a single Aegean island, with the exception of the larger island of Crete, where two ports are considered. The reference vessel used is a containership with a capacity of 150 small containers, a volume displacement of 20,100 m3 and a normal cruising speed of 12 knots under calm wind conditions. Supply and demand data for each port are presented in Table 2, along with the estimation of service times (berthing and handling time) per port. Service times are outputs of a function of port characteristics, demand and supply; these are derived empirically assuming that a number of small containers must be loaded or unloaded and transferred within the port (Karlaftis et al., 2009). As for time deadlines, a limit of 35–40 h after departure for supplying an island is considered as acceptable by local authorities. In this application, a 40 h time deadline is selected. Travel time is assumed to be a stochastic variable; a confidence level a = 90% is initially assumed for the chance constrained formulation (expected travel time per route with a confidence of 90%). It should be noted that such assumption provokes increases of travel times from 4.5% to 12% compared to corresponding deterministic ones. 5.2. GA parameters Different combinations of population size, crossover rate and mutation rate are selected for GA experiments. Two population sizes (26 and 50), crossover rates of 0.2, 0.4, 0.5 and 0.6 and mutation rates 0.05, 0.1 of 0.15 are assumed and 24 experiments are formed and used for GA calibration. The GA is terminated when there is no significant improvement to the fitness function value (<1% for a number of 20,000 generations). Due to time deadlines, a penalty parameter is used, if a vessel arrives later than time deadlines. In this case, the fitness function is penalized by pen = 50% for delays in arrival time. A relatively high percentage is selected, in order to avoid large excess of time deadlines. 6. Results Table 3 shows best GA results (fitness function) for each experiment, the algorithm’s running time and the number of routes derived. It is noted that each experiment was executed for 10 times and the best value per set of 10 executions as also the average value over the multiple runs were kept.

151

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155 Table 2 Demand, supply and service time per island port. Destination

Demand (small containers)

Supply (small containers)

Service time (h)

Rethimno Iraklio Kithira Milos Kea Thira Paros Naxos Syros Mykonos Tinos Andros Chios Ikaria Samos Astypalea Karpathos Rhodes Kos Kalymnos Skyros Mytilini Limnos Amorgos Sifnos Leros

12 70 5 10 7 33 15 16 12 37 9 3 22 5 21 9 8 58 23 11 2 46 21 5 4 9

7 57 1 43 0 2 19 2 8 1 1 1 9 1 35 1 1 15 3 0 0 21 6 1 0 1

1.8 6.1 1.3 3.1 1.3 2.4 2.4 1.7 1.8 2.5 1.4 1.2 2.2 1.2 3.2 1.4 1.4 3.9 2.0 1.4 1.1 3.7 2.1 1.2 1.2 1.4

Table 3 Results for various GA parameters trials. Number of trial

Population

Crossover rate

Mutation rate

Best fitness function value

Average fitness function value

Maximum running time (min)

Number of routesvessels

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

26 26 26 26 26 26 26 26 26 26 26 26 50 50 50 50 50 50 50 50 50 50 50 50

0.20 0.20 0.20 0.40 0.40 0.40 0.50 0.50 0.50 0.60 0.60 0.60 0.20 0.20 0.20 0.40 0.40 0.40 0.50 0.50 0.50 0.60 0.60 0.60

0.05 0.10 0.15 0.05 0.10 0.15 0.05 0.10 0.15 0.05 0.10 0.15 0.05 0.10 0.15 0.05 0.10 0.15 0.05 0.10 0.15 0.05 0.10 0.15

257.94 258.87 272.07 265.54 258.90 269.13 252.51 268.85 257.73 263.99 265.35 273.69 254.99 251.62 267.51 254.97 263.95 288.26 254.97 260.68 283.59 257.39 254.97 274.94

263.80 267.15 275.98 268.75 266.52 283.97 265.35 271.21 286.08 276.83 278.73 279.71 265.46 268.69 277.00 263.02 267.24 291.87 266.91 263.84 289.83 261.28 267.38 288.29

9.20 7.00 15.20 16.10 19.80 18.60 7.40 6.00 18.90 9.40 7.20 17.30 18.30 20.7 18.70 19.50 18.50 19.20 13.40 13.70 12.40 12.10 11.20 9.30

4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

Results of the various experiments do not reflect any significant variability (coefficient of variation = 4% for best values per parameters combination) and therefore can be characterized as adequately robust. The best result was obtained for experiment #14, whose outcome is illustrated in Fig. 11. Due to the selected vessel capacity, four (4) routes are required in order to cover the needs of the ports. As for the spatial dispersion of routes, each route of the best set serves in general ports from the same Aegean region. Furthermore, due to the algorithm’s structure and soft time windows, the number of island ports included in each route does not differ significantly (6–8 ports per route). In addition, the high capacity of a vessel has an effect not only on the emergence of extended routes but also on delays. An average delay of 5.8 h is observed for the best set of routes (maximum delay for route to Dodecanese islands).

152

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

Fig. 11. Optimal ship routes.

It is also noted that the best solution shown in Fig. 11 yields crossing route segments for two of the VRP routes. A property of the Euclidian TSP optimal tour is that it should not include crossing edges, since, when a closed polygonal line (the TSP tour) includes edges that cross, there exists always another polygonal line with the same number of vertices and smaller length. However, in our case, additional constraints such as capacity constraints combined with pick-ups and deliveries and time deadlines, as well as asymmetric travel costs along route segments (due to weather impacts), affect service priority and edge costs among candidate ports to be visited and therefore this property does not necessarily hold. 6.1. Sensitivity analysis Sensitivity analysis is performed for investigating impacts of vessel capacity, time deadlines, penalty parameter and the chance constraint formulation confidence level. As shown in Fig. 12, when vessel capacity increases, the number of routes decreases (a containership can serve more island ports), while average total delay per route shows an increasing trend due to extended length of routes. Besides, the capacity constraint is the vital factor during the production of routes, in opposite with time deadlines which can modestly be violated. In addition, while time deadlines become larger (Fig. 12), the average delays decrease because the vessel has more time available to serve the islands of the network and violations of time constraints are getting shorter. Penalizing for delays does not guarantee their elimination but the level of expected delays is highly affected by the penalty parameter (Fig. 12). For example, when a time deadline of 35 h and vessel capacity of 150 small containers are adopted, the algorithm provides results as follows: smaller penalty percentages yield higher average delays, while for the strictest value of parameter (a = 1), average delay is the half compared to value of 0.05. Last, we examine the performance of the algorithm for a stricter confidence level (Fig. 12). A level corresponding to 95% of the distribution of travel time results leads, as expected, to an increase of the fitness measure. In contrast, average delay decrease since the new level leads to a redesign of routes. In this case, extended routes include smaller numbers of nodes in comparison to the initial route design. 6.2. Comparison with past work We further compare model outputs with those of the deterministic approach proposed by Karlaftis et al. (2009), for the same network of ports and problem parameters (vessel capacities, penalties and so on, time deadlines). Table 4 compares the results (quality-exactness) for the same problem under the deterministic and stochastic cases. Within the scope of the

153

360.00

8.00

340.00

7.00

Average delay per route (hours)

Fitness Function Value

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

320.00 300.00 280.00 260.00 240.00 220.00 200.00 60.00

110.00

160.00

6.00 5.00 4.00 3.00 2.00 1.00 0.00 60.00

210.00

110.00

Ship capacity

Average delay per route (hours)

Number of routes

7.00 6.00 5.00 4.00 3.00 2.00 1.00 100.00

150.00

10.00 8.00 6.00 4.00 2.00 0.00 27.00

200.00

Ship capacity

32.00

37.00

42.00

47.00

Time deadline (h)

7.00

270.00

6.00

Fitness Function Value

Average delay per route (hours)

210.00

12.00

8.00

0.00 50.00

160.00

Ship capacity

5.00 4.00 3.00 2.00 1.00 0.00 0.00

265.00 260.00 255.00 250.00 245.00 240.00

0.50

1.00

Penalty parameter

0.90

0.95

Confidence level

Fig. 12. Sensitivity analysis graphs for various scenarios (time deadlines, penalty, confidence level).

Table 4 Comparison with Karlaftis et al. (2009). Results

Deterministic model (Karlaftis et al., 2009)

Stochastic model

Fitness function value

260.22

271.07

Number of routes

5.00

5.00

Average delay per route (h)

0.00

2.05

Comparison

154

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

stochastic model, the computation of the average delay is more exact and closer to real routing conditions. In particular, results for a confidence level of 90% are summarized in Table 4. Fitness function value shows an increase of 4.2% in the case of the stochastic model. Also, the average delay is significantly increased, as expected due to tighter consideration of sailing times. Furthermore, it is noted that the spatial distribution of routing slightly differs, since routes to larger and more distant islands include less service nodes, while routes to islands near Piraeus serve more nodes. Overall, the stochastic approach provides more conservative outputs but manages to incorporate realistic parameters that the corresponding deterministic models cannot consider. 7. Conclusions In this paper, we investigated a containership problem considering pick-ups, deliveries and time deadlines and taking into account the impact of weather conditions to the stochastic nature of travel times. We formulated the problem as a chanceconstrained model and developed a genetic algorithm approach for solving the problem. The algorithm was implemented for a homogeneous fleet of containerships transporting cargo in the Aegean Sea, in an effort to establish a set of optimal routes. Consideration of stochastic travel times leads to increased expected values ranging from 4.5% to 12% compared to their deterministic counterparts. GA computational experiment showed that a small fleet of containerships with capacities of 150 containers can provide efficient services, with an expected average excess of time deadline of about 5.4 h per route. Such delays are tolerable given the nature and sensitivity of carried goods with respect to time. Sensitivity analysis showed that results are consistent in changes in major parameters, such as vessel capacity and time deadline, reflecting practical significance in real-life conditions. The adoption of a tighter confidence level respecting travel times consequently provides routes characterized by bigger duration and smaller delays. All instances were solved between 7 and 20 min, so the proposed approach is adequate for solving practical problems of similar scale. With respect to future extensions, it should be noted that there could be cases that time deadlines may need to be violated, in favor of destination ports with increased importance (for example, in cases that higher than expected demand for serving a port is encountered for a particular reason or season-related factors). In such cases, the solution approach (and particularly the penalty function) could be redesigned for relaxing time deadlines for some ports. This could be combined with promoting routes which access high demand ports by priority. Future work could also focus on how to improve determination of stochastic travel times, in combination with the introduction of new stochastic variables, such as service times. Thus, the problem’s complexity is increased, but real-life conditions are better represented in the context of the liner shipping industry. Acknowledgements The paper is devoted to the memory of our advisor and co-author, Prof. Matthew Karlaftis, who met an untimely death in June 2014 at the age of 45. References Agarwal, R., Ergun, Ö., 2008. Ship scheduling and network design for cargo routing in liner shipping. Transport. Sci. 42 (2), 175–196. Álvarez, J.F., 2009. Joint routing and deployment of a fleet of container vessels. Marit. Econ. Logist. 11 (2), 186–208. Bagiorgas, h.S., Mihalakakou, g., Rehman, s., Al-Hadhrami, l.M., 2012. Offshore wind speed and wind power characteristics for ten locations in Aegean and Ionian Seas. J. Earth Syst. Sci. 121 (4), 975–987. Baker, B.M., Ayechew, M.A., 2003. A genetic algorithm for the vehicle routing problem. Comput. Oper. Res. 30 (5), 787–800. Blanton Jr., J.L., Wainwright, R.L., 1993. Multiple vehicle routing with time and capacity constraints using genetic algorithms. In: Proceedings of the 5th International Conference on Genetic Algorithms. Morgan Kaufmann Publishers Inc., pp. 452–459. Bowditch, N., The national imagery and mapping agency, 2002. The American Practical Navigator: An Epitome of Navigation. Bräysy, O., Gendreau, M., 2001. Genetic algorithms for the vehicle routing problem with time windows. Arpakannus 1, 33–38. Bräysy, O., Dullaert, W., Gendreau, M., 2004. Evolutionary algorithms for the vehicle routing problem with time windows. J. Heurist. 10 (6), 587–611. Chen, D., Padhy, C.P., 2010. Development of a ship weather-routing algorithm for specific application in North Indian Ocean region. In: The International Conference on Marine Technology BUET, Dhaka, Bangladesh, pp. 21–27. Christiansen, M., Fagerholt, K., Ronen, D., 2004. Ship routing and scheduling: status and perspectives. Transport. Sci. 38 (1), 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. Christiansen, M., Fagerholt, K., Nygreen, B., Ronen, D., 2013. Ship routing and scheduling in the new millennium. Eur. J. Oper. Res. 228 (3), 467–483. Chuang, T.N., Lin, C.T., Kung, J.Y., Lin, M.D., 2010. Planning the route of container ships: a fuzzy genetic approach. Expert Syst. Appl. 37 (4), 2948–2956. Cordeau, J.F., Laporte, G., Mercier, A., 2001. A unified tabu search heuristic for vehicle routing problems with time windows. J. Oper. Res. Soc. 52 (8), 928– 936. Davis, L. (Ed.), 1991. Handbook of Genetic Algorithms, vol. 115. Van Nostrand Reinhold, New York. Dethloff, J., 2001. Vehicle routing and reverse logistics: the vehicle routing problem with simultaneous delivery and pick-up. OR Spektrum 23, 79–96. Eiben, A.E., Smith, J.E., 2003. Introduction to Evolutionary Computing. Springer, Berlin. Gendreau, M., Hertz, A., Laporte, G., 1994. A tabu search heuristic for the vehicle routing problem. Manage. Sci. 40 (10), 1276–1290. Gendreau, M., Laporte, G., Séguin, R., 1996. Stochastic vehicle routing. Eur. J. Oper. Res. 88 (1), 3–12. Gendreau, M., Laporte, G., Potvin, J-Y., 1999. Metaheuristics for the Vehicle Routing Problem. Technical Report G-98-52, GERAD. Gendreau, M., Laporte, G., Potvin, J.Y., 2002. Metaheuristics for the capacitated VRP. In: The Vehicle Routing Problem, vol. 9. pp. 129–154. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. University of Michigan Press, Ann Arbor. Imai, A., Shintani, K., Papadimitriou, S., 2009. Multi-port vs. hub-and-spoke port calls by containerships. Transport. Res. Part E 45, 740–757. Jih, W.-R., Hsu, J.Y.-J., 2004. A family competition genetic algorithm for the pickup and delivery problems with time window. Bull. College Eng. Natl. Taiwan Univ. 90 (89–98), 2004.

K. Kepaptsoglou et al. / Transportation Research Part C 55 (2015) 139–155

155

Kaldelis, K.I., 1999. Wind Energy Administration. Stamoulis Publications, Athens. Karlaftis, M.G., Kepaptsoglou, K., Sambracos, E., 2009. Containership routing with time deadlines and simultaneous deliveries and pick-ups. Transport. Res. Part E 45 (1), 210–221. Kjeldsen, K., 2011. Classification of ship routing and scheduling problems in liner shipping. INFOR 49 (2), 139–152. Kwon, Y.J., 1981. The Effect of Weather Particularly Short Sea Waves on the Ship Speed Performance. PhD Thesis, University of Newcastle, Newcastle Upon Tyne. Laporte, G., Gendreau, M., Potvin, J.-Y., Semet, F., 2000. Classical and modern heuristics for the vehicle routing problem. Int. Trans. Oper. Res. 7, 285–300. Lawrence, S.A., 1972. International Sea Transport: The Years Ahead. Lexington Books, Lexington, MA. Li, X., Tian, P., Leung, S.C., 2010. Vehicle routing problems with time windows and stochastic travel and service times: models and algorithm. Int. J. Prod. Econ. 125 (1), 137–145. Luke, S., 2013. Essentials of Metaheuristics. Lulu. Mazarakis, N., Kotroni, V., Lagouvardos, K., Bertotti, L., 2012. High-resolution wave model validation over the Greek maritime areas. Nat. Hazards Earth Syst. Sci. 12 (11), 3433–3440. Meng, Q., Wang, T., 2011a. A scenario-based dynamic programming model for multi-period liner ship fleet planning. Transport. Res. Part E 47, 401–413. Meng, Q., Wang, X., 2011b. Intermodal hub-and-spoke network design: incorporating multiple stakeholders and multi-type containers. Transport. Res. Part B 45 (4), 724–742. Meng, Q., Wang, S., Andersson, H., Thun, K., 2013. Containership routing and scheduling in liner shipping: overview and future research directions. Transport. Sci. 48 (2), 265–280. Mole, R.H., Jameson, S.R., 1976. A sequential route-building algorithm employing a generalised savings criterion. Operat. Res. Quart., 503–511 Prins, C., 2004. A simple and effective evolutionary algorithm for the vehicle routing problem. Comput. Oper. Res. 31, 1985–2002. Rego, C., 2001. Node-ejection chains for the vehicle routing problem: sequential and parallel algorithms. Parallel Comput. 27 (3), 201–222. Reinhardt, L.B., Pisinger, D., 2012. A branch and cut algorithm for the container shipping network design problem. Flex. Serv. Manuf. J. 24 (3), 349–374. Ronen, D., 1983. Cargo ships routing and scheduling: survey of models and problems. Eur. J. Oper. Res. 12, 119–126. Ronen, D., 1993. Ships scheduling: the last decade. Eur. J. Oper. Res. 71 (3), 325–333. Russell, R.A., Chiang, W.C., 2006. Scatter search for the vehicle routing problem with time windows. Eur. J. Oper. Res. 169 (2), 606–622. Sambracos, E., 2000a. Exploring operational problems of the goods supply chain in the Greek islands: towards a reengineering of the system, repositioning logistics. In: Proceedings of the 16th International Logistics Conference, Versailles, France. Sambracos, E., Paravantis, J.A., Tarantilis, C.D., Kiranoudis, C.T., 2004. Dispatching of small containers via coastal freight lines: the case of the Aegean Sea. Eur. J. Oper. Res. 152, 365–381. Shintani, K., Imai, A., Nishimura, E., Papadimitriou, E., 2007. The container shipping network design problem with empty container repositioning. Transport. Res. Part E 43, 39–59. Ship and Marine Hydrodynamics Laboratory, National Technical University of Athens, 1992. Wind and Wave Atlas of the North-eastern Mediterranean Sea. Authority of Hellenic Navy General Staff, Athens. Sigurd, M.M., Ulstein, N.L., Nygreen, B., Ryan, D.M., 2005. Ship scheduling with recurring visits and visit separation requirements. In: Desaulniers, G., Desrosiers, J., Solomon, M.M. (Eds.), Column Generation. Springer, US, pp. 225–245. Taillard, É., 1993. Parallel iterative search methods for vehicle routing problems. Networks 23 (8), 661–673. Thangiah, S., 1995. Vehicle routing with time windows using genetic algorithms. In: Chambers, L. (Ed.), Application Handbook of Genetic Algorithms: New Frontiers, vol. II. CRC Press, Boca Raton. Tran, N.K., Haasis, H.D., 2013. Literature survey of network optimization in container liner shipping. Flex. Serv. Manuf. J., 1–41 Vernimmen, B., Dullaert, W., Engelen, S., 2007. Schedule unreliability in liner shipping: origins and consequences for the hinterland supply chain. Marit. Econ. Logist. 9 (3), 193–213. Wang, S., Meng, Q., 2012a. Liner ship route schedule design with sea contingency time and port time uncertainty. Transport. Res. Part B: Methodol. 46 (5), 615–633. Wang, S., Meng, Q., 2012b. Robust schedule design for liner shipping services. Transport. Res. Part E: Logist. Transport. Rev. 48 (6), 1093–1106. Wark, P., Holt, J., 1994. A repeated matching heuristic for the vehicle routeing problem. J. Oper. Res. Soc., 1156–1167 Wiley, V.D., 2000. The Symmetric Group Class User’s Manual for Partitioning and Ordering Problems. The University of Texas at Austin, TX, USA. Zafirakis, D., Gavrilopoulou, E., Kavadias, K.A., Kaldellis, J.K., 2012. The need for the development of a new readjusted Weibull distribution for increased reliability of energy yield estimation. In: Online Conference Proceedings of EWEA. Zhang, T., Chaovalitwongse, W.A., Zhang, Y., 2012. Scatter search for the stochastic travel-time vehicle routing problem with simultaneous pick-ups and deliveries. Comput. Oper. Res. 39 (10), 2277–2290.