- Email: [email protected]
Minimizing the fuel consumption and the risk in maritime transportation: A bi-objective weather routing approach
Accepted Manuscript
Minimizing the fuel consumption and the risk in maritime transportation: a bi-objective weather routing approach Aphrodite Veneti, Angelos Makrygiorgos, Charalampos Konstantopoulos, Grammati Pantziou, Ioannis Vetsikas PII: DOI: Reference:
S0305-0548(17)30193-4 10.1016/j.cor.2017.07.010 CAOR 4291
To appear in:
Computers and Operations Research
Received date: Revised date: Accepted date:
27 September 2016 6 July 2017 13 July 2017
Please cite this article as: Aphrodite Veneti, Angelos Makrygiorgos, Charalampos Konstantopoulos, Grammati Pantziou, Ioannis Vetsikas, Minimizing the fuel consumption and the risk in maritime transportation: a bi-objective weather routing approach, Computers and Operations Research (2017), doi: 10.1016/j.cor.2017.07.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • A time-dependent algorithm is presented for the ship weather routing
CR IP T
problem. • Given a maximum travel duration constraint, fuel consumption and risk are minimized.
• Improvements of the approach are proposed for faster execution but similar solutions.
AN US
• Techniques are also presented that improve the practicality of the derived
AC
CE
PT
ED
M
solutions.
1
ACCEPTED MANUSCRIPT
Minimizing the fuel consumption and the risk in maritime transportation: a bi-objective weather routing approach
CR IP T
Aphrodite Venetia,∗, Angelos Makrygiorgosc , Charalampos Konstantopoulosa , Grammati Pantzioub , Ioannis Vetsikasc a Department
of Informatics, University of Piraeus, Greece of Informatics, Technological Educational Institute of Athens, Greece c Institute of Informatics & Telecommunications, NCSR Demokritos, Agia Paraskevi, Greece
AN US
b Department
Abstract
The paper presents an improved solution to the ship weather routing problem based on an exact time-dependent bi-objective shortest path algorithm. The two objectives of the problem are the minimization of the fuel consumption and the total risk of the ship route while taking into account the time-varying sea and
M
weather conditions and an upper bound on the total passage time of the route. Safety is also considered by applying the guidelines of the International Maritime
ED
Organization (IMO). As a case study, the proposed algorithm is applied for finding ship routes in the area of the Aegean Sea, Greece. Enhancements of the proposed algorithm are also presented which improve the efficiency of our
PT
approach.
Keywords: Ship weather routing, Multi-criteria optimization, Time
CE
dependent networks, Resource-constrained shortest path.
AC
1. Introduction Maritime transportation is crucial for the economic growth of many coun-
tries as well as for the global economy, in general. During the past two decades, ∗ Corresponding
author Email address: [email protected] (Aphrodite Veneti)
Preprint submitted to Computers & Operations Research
July 19, 2017
ACCEPTED MANUSCRIPT
there has been a growing interest in methods for calculating optimal ship tra5
jectories. The traditional economic objective was to minimize the total traveled distance while later, additional energy efficiency objectives were added such as
CR IP T
the minimization of the fuel consumption and gas emissions. Besides energy efficiency, navigation safety is also considered since sea transportation poses a high environmental risk due to the possibility of maritime accidents, especially if haz10
ardous materials (i.e. chemicals or oil) are transferred. Moreover, ship routing methods have to adhere to the safety guidelines specified by the International Maritime Organization [1].
AN US
The problem of optimal route planning in maritime transportation (or ship
routing problem) takes into consideration the different objectives and constraints 15
set by the ship owners, the national regulations, the international organizations, etc. Regardless of the specific objectives/constraints, the main factor that makes the ship routing problem difficult to solve, is the time-varying weather conditions [2]. The parameters and cost factors of the problem as well as the optimization
20
M
criteria are strongly affected by the weather and sea conditions and therefore, finding an optimal ship route should take into consideration the weather fore-
ED
casts and sea conditions during the course of the route. To accomplish this, a reliable model that predicts the ship’s response to the weather and sea conditions is required, while an effective solution approach for the ship routing problem
25
PT
should determine the heading control and ideally the optimal ship power settings according to the weather and sea conditions. Thus, broadly speaking, the
CE
ship routing problem is a dynamic, time-dependent problem which is usually referred to as the ship weather routing problem. The optimization objectives in the ship routing problem are usually the
AC
minimization of the voyage time, fuel consumption and voyage risk. The ship
30
routing algorithms appeared in the literature are classified into two categories, namely the exact and the heuristic approaches. The exact algorithms [3, 4, 5, 6, 7, 8, 9, 10, 2, 11, 12, 13, 14] determine the optimal solution at the expense of the execution time while heuristic approaches [15, 16, 17, 18, 19, 20] run faster but search only inside a subspace of the search space for the best solution. 3
ACCEPTED MANUSCRIPT
35
Also, for coastal navigation and trans-oceanic seafaring, the methodology of the calculus of variations [21] has also been employed and is based on long-term weather forecasts. Finally, a classic method for optimizing ship routes is still
CR IP T
the isochrone method [22]. In this paper, we consider the bi-objective time-constrained ship weather 40
routing problem as a time-dependent bi-objective point-to-point shortest path
problem with a constraint on the total voyage time. The two objectives of the problem are the minimization of the fuel consumption and the total risk. In
practice, frequent change of ship power or equivalently of nominal speed during
45
AN US
the voyage is avoided [3] especially in the case of coastal navigation. Thus, in
this work, the navigation speed is not a control variable; it is assumed to be constant during the whole journey and should be given as an input parameter. Notably, when the objective is the fastest route between two ports, the networks modelling the sea transportation commonly have the FIFO property1 . However, in our problem setting, the objectives are different and leaving from a node immediately after the arrival may not be the best option, for instance, when
M
50
the costs (fuel consumption and risk) along the outgoing link are decreasing
ED
in the next time interval [23]. Although, computing shortest paths in FIFO networks is a polynomially solvable problem [24, 25], this is not the case when optimal waiting policy should be determined for other cost objectives [23]. On the other hand, frequent stops in the middle of the sea or alternating ship speed
PT
55
frequently is not a common practice in ship routing [3]. Thus, although it is
CE
generally better to wait at a node, in our maritime setting, we assume that there is neither the option of waiting nor the option of speed decrease; by decreasing the speed, we could simulate the forbidden waiting and shift the arrival at the node later exactly at the optimal departure time.
AC
60
To sum up, in this work, we address the time-dependent bi-objective shortest 1A
network has the FIFO property if all of its arcs have that property i.e., for each arc
(i, j) of the network, earlier departure from i always leads to earlier arrival at j. Therefore, the arrival events at j are in the same chronological order as the departure events at i.
4
ACCEPTED MANUSCRIPT
path problem on a network with objectives other than the travel time, with fixed departure time, no waiting at nodes, constant nominal speed and a constraint on the total travel time. Also, since the problem is time-dependent, it is assumed that all arc costs are deterministic functions of the time. However, considering
CR IP T
65
the low execution time of the proposed techniques, dynamic situations where the
weather changes unpredictably during the journey can be handled by running the algorithms afresh.
We first give a non linear integer programming formulation of the problem, 70
and in the sequel we present a new exact algorithm for solving it. To the best
AN US
of our knowledge, only a few algorithmic approaches for the multi-objective
time-dependent shortest path problem have been proposed in the literature [26],[27],[28] are not efficient when waiting at nodes is forbidden. We employ a forward label-setting approach to efficiently tackle the aforementioned time75
dependent bi-objective shortest path problem. The number of the labels being processed during the algorithm execution is kept at a reasonable level, reducing
M
the computational overhead. The algorithms is compared against the backward algorithm of Hamacher et al’s [28], the best so far algorithm in the literature
80
ED
that solves the problem at hand, and the experimental results confirm the faster execution time of the proposed algorithms. A preliminary version of the algorithm appears in [11].
PT
Furthermore, we present critical modelling issues that directly affect the performance of our problem solution approach. We introduce a novel grid structure
CE
that is used as a base layer for solving the ship weather routing problem and 85
significantly improves the efficiency of the proposed algorithm. We also proposed a number of methods to further improve the efficiency of our approach,
AC
namely: (i) a dynamically partitioned grid for reduced graph size as well as an alternative reduction technique where grid is pruned along a given voyage plan, (ii) a heuristic function used to transform the initial algorithm into a bi-
90
objective A* algorithm thereby reducing the search space and (iii) a technique for deriving loopless ship routes, thus making the solutions conformable to the usual practice in this industry. In order to test and validate our approach in a 5
ACCEPTED MANUSCRIPT
realistic setting, we applied the proposed algorithm as well as the speeding up techniques and improvements on the algorithm for the optimal ship route plan95
ing in the Aegean Sea, Greece. The experimental results confirm the efficiency
for some of the afrementioned techniques appear in [29].
CR IP T
of the approach and its applicability in a realistic scenario. Preliminary results
The rest of the paper is organized as follows. Section 2 overviews the related work. Section 3 presents a non-linear integer programming formulation 100
of the ship weather routing problem which is considered as a time-dependent bi-objective point-to-point shortest path problem, and gives an efficient exact
AN US
algorithm for solving the problem. Section 4 presents critical modeling issues while Section 5 includes the speeding up techniques and improvements on the
proposed algorithm. Section 6 presents the empirical evaluation of the proposed 105
algorithm and its enhancements. Finally, Section 7 concludes our work.
M
2. Related Work
The ship routing problem is a multi-objective, non-linear optimization problem with constraints which, due to the importance of its applications, has been
110
ED
widely studied in the literature. In [30] the common objectives of the optimization problems in the transportation industry are classified into three main
PT
categories: classical/economic objectives, environmental (ecosystem)/regional fairness/health objectives, and climate and sustainability objectives. Until recently, in ship navigation the traditional economic objective was to minimize
CE
the total distance travelled [31]. Other relevant objectives were the minimiza-
115
tion of the fuel consumption and gas emissions. However, both these metrics are very complex to estimate. Nowadays, beyond energy efficiency, navigation
AC
safety is also a main concern, since sea transport is the main carrier of crude oil and other hazardous materials. Furthermore, the proposed methodologies have to adhere to the safety guidelines specified by the International Maritime
120
Organization. Exact as well as heuristic approaches have been proposed in the literature for
6
ACCEPTED MANUSCRIPT
ship route planning. Exact algorithms derive optimal solutions of the problem, however, they usually have heavy computation demands. In contrast, a heuristic approach searches for a feasible solution only within a subspace of the search space hopefully, close to the optimal solution. Hence the execution time of a
CR IP T
125
heuristic algorithm is generally much lower than that of an exact algorithm. Since long-term weather forecasts are available, calculus of variations methods
have also been proved suitable for both coastal navigation and trans-oceanic
seafaring [32, 21, 33]. Another old but popular method for optimizing ship 130
route is the isochrone method [34, 22].
AN US
Dynamic Programming and Label Setting methods are commonly met in exact approaches for solving the single objective ship routing problem with or without constraints. In [8] the Optimum Track Ship Routing (OTSR) system is presented, an automation of the weather ship routing service provided by the 135
US navy. A binary heap version of Dijkstra’s algorithm determines the optimum route taking into consideration the weather and sea conditions. Safety of navi-
M
gation is ensured by eliminating grid points where wind speed and wave height exceed ship’s predefined limits. Based on dynamic programming, the authors
140
ED
in [3] presented a modified version of the maze routing method first presented in [35]. The main objective is the minimization of the travel time while turn penalties and collision avoidance techniques are considered for enhancing the
PT
navigation safety. The proposed approach uses relatively few control variables and is space and time efficient. However, it does not take into account the exist-
CE
ing weather and sea conditions. In [9] a method based on Dijkstra’s algorithm is 145
proposed for finding routes with minimum fuel consumption. The input grid is built starting from a given voyage plan and then additional nodes are added on
AC
lines perpendicular to the input route. The authors in [10] incorporate Dreyfus approach [36] so as to solve the ship routing problem in a dynamic environment using as input, measurements of the ocean wave-field surrounding the vessel.
150
In [5], a forward 3-dimensional dynamic programming method was employed where both ship speed and ship course are control variables. The objective is to minimize the fuel consumption while taking into account geographical, 7
ACCEPTED MANUSCRIPT
weather, operational and safety constraints. It was shown that optimizing only the ship course is not as effective as optimizing both the course and the speed. 155
In [4] a Decision Support System is presented which uses meteo-marine and
CR IP T
oceanographic operational information data for all relevant environmental field variables (wave height, direction, period and currents). It uses a modification
of Dijkstra’s algorithm in order to find the optimal route with respect to the
total navigation time in the presence of time-varying fields. Safety restrictions 160
are also taken into account.
Another exact method with fast convergence for solving the ship routing
AN US
problem is the A∗ search algorithm. In [37] a two-tiered A∗ approach is employed
where firstly, the A∗ algorithm is applied to a macro data map to find a nondetailed path from the origin to the destination, and then the A∗ algorithm 165
is applied to obtain an optimal path for every edge of the suboptimal path computed in the first phase, using micro map data. The authors in [14] run the A∗ algorithm on a grid with eight neighbors per node for finding paths
M
optimized in terms of the navigation time subject to safety distance and turnradius constraints and in the presence of obstacles. The heuristic information of A∗ is the shortest possible travel time between a node and the destination. As
ED
170
the path given by the A∗ algorithm is piece-wise linear, the algorithm applies a smoothing technique on this path so that the final route is a realistic route that
PT
can be followed by ships in practice. In [13] the A∗ algorithm is employed to tackle the route planning problem in the presence of obstacles. The objective is varying depending on the particular scenario at hand (e.g. the minimization
CE
175
of the route length, the number of manoeuvres, etc.). The authors adapt the A∗ algorithm to find the least cost route using the exact distance between the
AC
nodes as the heuristic function. The risk of the path is also taken into account.
180
The main heuristic methods that have been employed in the literature for
solving the ship routing problem, are the evolutionary computation and the simulated annealing method. The fast convergence of these methods to a near optimal solution makes them suitable for solving the ship routing problem in real-time. Despite the potential of these heuristic methods in solving multi8
ACCEPTED MANUSCRIPT
objective problems, they were mainly used for solving the single objective ship 185
routing problem and for minimizing a weighted sum of multiple different objectives [15]. In the sequel we briefly refer to the research works which employ
CR IP T
heuristic methods for solving the multi-objective ship weather routing problem. In [19] a heuristic method is presented for multi-criteria ship weather routing optimization with respect to the voyage time, the fuel consumption, and the 190
voyage risk. The approach is based on the Strength Pareto Evolutionary Algo-
rithm (SPEA), which was designed for combinatorial and continuous function
multiple objective optimization problem instances. The solution routes are de-
AN US
fined as a sequence of waypoints from the origin to the destination. The initial population consists of a set of basic routes: an orthodrome, a loxodrome, a 195
time-optimized isochrone route and a route given by applying fuel-optimization to the time-optimized isochrone route. In [17] a Genetic Algorithm is employed for retrieving a Pareto set of optimal ship routes with respect to the voyage time and the fuel consumption. The initial population is chosen as a random sequence
200
M
of grid nodes. The authors take into account the time dependency of the sea state and wind field applying linear interpolation in time when measurements
ED
are missing. In [20] an evolutionary algorithm is proposed for obtaining a number of low cost routes which are also of high safety and avoid narrow waters, bad weather, foggy areas, congested and fishing areas. The initial candidate
205
PT
route population is automatically generated by GIS. In [38] a genetic algorithm for solving the bi-objective time-constrained ship weather routing problem is
CE
presented. The derived routes are optimized with respect to two conflicting objectives, the total fuel consumption and the mean risk of the route while taking
AC
into account the weather conditions.
210
Route planning in the maritime setting is a complex problem due to a number
of dynamic phenomena and navigation constraints. Time-dependent algorithms can efficiently handle the inherently time-varying environment of ship routing. Despite the importance of the time-dependent shortest path problem and its applications in maritime navigation, there are relatively few works in the literature on these issues ([39, 17, 19, 5, 10, 11]). In contrast, most literature on 9
ACCEPTED MANUSCRIPT
215
ship routing and scheduling focuses on the time-invariant shortest path problem [40, 21, 8, 3, 2, 9, 37, 14, 13]. However, there are also works which consider the inherent uncertainty in the behaviour of ocean currents and weather conditions
to minimize the expected fuel consumption [6, 7].
CR IP T
and treat ship routing as a stochastic optimization problem whose objective is
In this paper, we propose a new efficient exact algorithm for the bi-objective
220
time-constrained ship weather routing problem. To the best of our knowledge, it is the first time that an exact multi-objective time-dependent algorithm is employed for solving the ship weather routing problem. Apart from finding op-
225
AN US
timal routes with respect to economical objectives, we are interested in reducing the possibility of an accident taking also care of marine resources conservation. The algorithm returns the whole Pareto set and the derived routes are opti-
mized with respect to two conflicting objectives, the total fuel consumption and the mean risk of the route while taking into account the prevailing weather
230
M
conditions.
3. The bi-objective time-constrained ship weather routing problem
ED
As already mentioned, we consider the bi-objective time-constrained ship weather routing problem as a time-dependent bi-objective time-constrained short-
PT
est path problem. We deal with the most realistic point-to-point instance of the problem, where the ship speed is not a control variable and waiting at nodes is 235
not permitted. In this section, we first give a non linear integer programming
CE
formulation of the problem, and then we present a new forward label-setting time-dependent bi-objective shortest path algorithm to solve the problem. Let G = (V, A) be a directed graph, where V is the set of nodes (|V | = n)
AC
and A is the set of edges (|A| = m). In the bi-objective time-constrained ship
240
weather routing problem, there are three time-dependent positive costs c1ij (t), c2ij (t) and c3ij (t) associated with each edge (i, j) ∈ A, namely, the travel time, the fuel consumption and the risk respectively, incurred along this edge when leaving from node i at time t. We follow the frozen arc model of [23], where
10
ACCEPTED MANUSCRIPT
the above edge costs are assumed to be constant during the traversal of the 245
edge. Let also tstart denote the departure time from source node s ∈ V , and T rT (P ), F C(P ), R(P ) denote the total travel time, fuel consumption and risk
CR IP T
respectively, accumulated along a path P . The problem is to find a shortest path P between a source s ∈ V and a
destination d ∈ V that minimizes the two objectives F C(P ) and R(P ), leaving 250
from the source node at time tstart , while not violating the constraint T rT (P ) ≤
T , where T is the maximum acceptable duration of a route. In the following
formulation, we assume an upper bound L on the number of edges a solution
AN US
path may contain. This constraint is necessary since in the case of forbidden waiting at nodes, the optimal solution may contain infinite number of edges 255
at least in principle, as a result of infinite loops created along the route in an attempt to simulate waiting at nodes which is not permitted [23]. Thus, actually, the derived solution may be a walk and not necessarily a path. Since, this potentially infinite number of edges cannot be expressed by the finite number
260
M
of variables in the above network flow problem, an additional constraint on the number of edges in a solution path should be placed.
ED
We have used a number of variables in the formulation below. Specifically, xi ∈ {0, 1} with xi = 1 when the solution path has exactly i edges. We also
define the binary variable erij with erij = 1 when the edge eij is the r-th edge
265
PT
starting from s along the route. We also define trj as the arrival time at node j through the edge eij when erij = 1. If such an edge does not exist, trj is set to a
AC
CE
large value. We also set t0s = tstart .
11
ACCEPTED MANUSCRIPT
r=1...L (i,j)∈A
X
R(P ) =
X
c3ij (tr−1 )erij i
r=1...L (i,j)∈A
X
T rT (P ) =
X
r=1...L (i,j)∈A
X
c1ij (tir−1 )erij ≤ T
e1sj = 1
AN US
(i,j)∈A
erij ≤ 1, r = 2 . . . L
X
erdj = 0, r = 2 . . . L
j∈V −{d}
X
xi = 1
i=1..L
erjd = xr , r = 1 . . . L
M
X
j∈V −{d}
X
erji −
X
j∈V
er+1 = 0, r = 1 . . . L, i ∈ V − {d} ij
ED
j∈V
(2)
trj =
X
PT CE AC
X
(4)
(6) (7) (8) (9) (10)
erij (tr−1 + c1ij (tr−1 )) i i
i∈{k|(k,j)∈A}
+ (1 −
(3)
(5)
j∈V −{s}
X
(1)
CR IP T
min z = (F C(P ), R(P )) X X F C(P ) = c2ij (tr−1 )erij i
i∈{k|(k,j)∈A}
(11) erij )M, j ∈ V, r = 1 . . . L
t0s = tstart
(12)
cqij (M )
(13)
= M0
erij , xr ∈ {0, 1}, tri ∈ R+ , i, j ∈ V, r = 1 . . . L
Constraint (4) requires that the total travel time of the solution path must
not exceed the maximum acceptable duration T . Constraint (5) states that the first edge of the solution path is an outgoing edge of the start node s, while
270
constraint (6) states that at most one of the edges of the graph can be the r-th edge of the solution path. Constraint 7 ensures that the solution path cannot 12
ACCEPTED MANUSCRIPT
extend beyond the destination node d regardless of its actual length. Constraint (8) determines the length of the solution path. Constraint 9 ensures that the solution path will end up at the destination node d after traversing the right number of edges. Constraint (10) is the ordinary flow conservation constraint
CR IP T
275
ensuring that there is only one incoming and one outgoing edge at each node
along the solution path and also that these edges are successive in the path. Equation (11) gives the arrival time at node j when this node is the head of the
r-th edge of the solution path. In this equation, we assume that M is a large 280
number and also that cqij (M ) = M 0 (q = 1 · · · 3) where M 0 is a large number
AN US
again.
Each solution z contains the values of the two objective functions for the corresponding path. Also, although the proposed Integer Programming model is non linear, when the time dependent edge costs functions c1ij (t), c2ij (t) and 285
c3ij (t) are linear, the above formulation can be converted into a linear model [41]. Since the objective functions may be conflicting, a single solution that simul-
M
taneously optimizes each objective may not exist and therefore, the problem is actually to find the set of non dominated solutions (Pareto optimal solutions).
290
ED
A non dominated solution is the solution where none of the objectives can be improved without worsening the others. Formally, a path p is said to dominate another path p0 (p ≺ p0 ) if (F C(p) < F C(p0 ) and R(p) ≤ R(p0 )) or
PT
(F C(p) ≤ F C(p0 ) and R(p) < R(p0 )).
In the sequel, we present a new forward label-setting algorithm (Algorithm
CE
1) for solving the time-dependent bi-objective shortest path problem with no 295
waiting at the nodes and finding the set of Pareto optimal paths. Noteworthily, although we handle the bi-objective case, the proposed algorithm turns out to
AC
be applicable in the case of multiple objectives, as well. The Ship Routing Algorithm 1 keeps a set of labels li (t) where i is a node and t is a time instance. Each label li (t) = (F C, R, j, prev ptr) corresponds to a path from the source
300
node s to the node i arriving at i at time t and comprises four components. Specifically, the total fuel consumption F C and risk R represents the total fuel consumption and risk of a so-far constructed path from s to i, the integer j is 13
ACCEPTED MANUSCRIPT
the predecessor of i on the path from s to i and the pointer prev ptr points to the Pareto optimal label of the predecessor j whose extension gave this label. As has already been mentioned, leaving as earlier as possible from a node
305
CR IP T
does not necessarily reduce the fuel consumption or risk of the outgoing edge. Thus, there is no easy way to know in advance the ideal time of leaving a node for obtaining the minimum cost path toward the destination node. Consequently,
for each node and for each realized arrival (and hence departure) time instance 310
at this node, we should keep all the non-dominated labels referring to that particular arrival time. In contrast, in the corresponding static bi-criteria shortest
AN US
path problem there is no need of partitioning the labels at each node according
to their arrival time and only a single set of non-dominated labels needs to be kept at each node.
The Ship Routing Algorithm 1 iterates over all integer time instances, with
315
step size of 1 minute, in the interval [tstart , tstart + T ], where tstart is the departure time from the source node s and T is the maximum permissible voy-
M
age duration. At each node, the algorithm maintains two groups of label lists namely, the permanent lists and temporary lists. Each of these lists corresponds to a certain realized arrival time and either contains the final Pareto-optimal
ED
320
paths “arriving” this node at that particular time instance or contains the sofar Pareto-optimal paths reaching this node at that instance. When the outer
PT
loop scans a time instance, all the labels of the temporary list of a node for that instance are transferred one by one to the corresponding permanent list. However, just before the transfer, these labels are extended thereby creating
CE
325
new labels. Each of these new labels is kept only if it does not violate the total travel time constraint (line 19) and is not dominated by another label belonging
AC
to the same node with the same arrival time (lines 22-38). Also, labels already in the temporary list for that arrival instance are compared against this newly
330
generated label and discarded when found to be dominated by the new label.
14
ACCEPTED MANUSCRIPT
Algorithm 1 The time-dependent bi-objective shortest path algorithm Require: G = (V, A), and C, the cost matrix for all edges (i, j) ∈ A Ensure: All Pareto optimal paths from the node s to the node d s: the start node d: the destination node tstart : the departure time T: the maximum allowed total travel time
CR IP T
1:
trt(i, j, t): the travel time for traversing the edge (i, j) departing from the node i at time instance t
AN US
tjar : an arrival time instance at the node j
li (t): a label of the node i corresponding to the path from s to i, arriving at i at time instance t.
Ltempi (t): the list of temporary labels of node i at time instance t Lpermi (t): the list of permanent labels of node i at time instance t card(li (t), Lpermi (t)): the position of li (t) in the list of permanent labels
M
of node i
th lp component of a label li (t) i (t): the p
ED
/* Initialization of temporary and permanent labels of every node and for all t ∈ {tstart , tstart + 1, ..., T + tstart } Ltempi (t), Lpermi (t) ← ∅,
3:
Ltemps (tstart ) ← {[(0, 0); null; null]}
4:
∀i ∈ V, ∀t ∈ {tstart , tstart + 1, ..., T + tstart }
for tar = tstart to T + tstart do
CE
5:
PT
2:
6:
AC
7:
8: 9:
10: 11: 12: 13:
*/
while (∪i∈V Ltempi (tar ) 6= ∅) do Select a pivot label li∗ (tar ) from ∪i∈V Ltempi (tar )
/* Remove li∗ (tar ) from Ltempi (tar ) and add it to Lpermi∗ (tar )
*/
Ltempi∗ (tar ) ← Ltempi∗ (tar )\{li∗ (tar )} Lpermi∗ (tar ) ← Lpermi∗ (tar ) ∪ {li∗ (tar )} /* Store the position of label li∗ (tar ) in the list Lpermi∗ (tar ) pos ← card(li∗ (tar ), Lpermi∗ (tar ))
if i∗ 6= d then
15
*/
ACCEPTED MANUSCRIPT
14:
/* Label all the successors of i∗
*/
∗
for all do(i , j) ∈ E
15:
tjar ← tar + trt(i∗ , j, tar )
16:
/* Check the total duration constraint
18:
g←
tjar
19:
h ← duration of the shortest path from node j to node d
20:
f ←g+h
− tstart
if tjar − tstart + h ≤ T then lj (tjar )
22:
r(i∗ , j, tar )); i∗ ; pos)
if 6
24:
∃lj0 (tjar )
∈
Ltempj (tjar )
:
˜ j (tjar ) lj0 (tjar )≺l
Ltempj (tjar )
←
Ltempj (tjar )
∪
*/
then
/* Store the label lj (tjar ) of node j at time instance tjar as temporary
26:
*/
{lj (tjar )}
/* Delete all temporary labels of node j at time instance tjar dominated by
M
27:
li∗ 2 (tar ) +
/* Check that there is no label lj0 (tjar ) of node j at time instance tjar domi-
nating label lj (tjar )
25:
((li∗ 1 (tar ) + f c(i∗ , j, tar ),
←
AN US
21:
23:
*/
CR IP T
17:
lj (tjar )
Ltempj (tjar )
ED
28:
Ltempj (tjar )
←
and
end if
29:
*/ Ltempj (tjar )\{lj0 (tjar )
lj (tjar )
≺
lj0 (tjar )}
end if
PT
30:
end for
31:
end if
CE
32:
end while
33:
end for
35:
/* Output:Lpermd , all Pareto optimal paths from the node s to the node d
AC
34:
*/
36:
return Lpermd
The order by which the labels of the current time instance are extended is not actually important, since all labels of that instance should be extended first 16
∈
ACCEPTED MANUSCRIPT
before proceeding to the next iteration. Therefore, the label to be extended (pivot label) can be selected randomly (line 8). By the end of the execution, 335
all Pareto-optimal solution paths will have been stored as permanent labels of
CR IP T
the destination node d. Notice also that for the node d, no label extension takes place (line 14) since the algorithm does not need to go further down the destination node due to the fact that the cost functions are non-negative-valued and hence going further only adds to the total cost of the path. Finally, Ship 340
Routing Algorithm 1 assumes a single departure time instance from the source
node. Apparently, the algorithm can be easily extended to handle the possibility
AN US
of more than one departure time instances.
Also, for improving the performance of the algorithm, a heuristic function, denoted by h, is used for estimating a lower bound of the travel time from the 345
current position to the node d. Summing the current cost of the path, denoted by g, and the h value, a lower bound, f , for the travel time from s to d is obtained. This estimation is used to check if the extension of the current, partial, path
M
will violate the constraint of the total travel time. Using this heuristic, we can prune paths from the very early stages of the algorithm, since it is certain that these paths will not reach d in time. The h function is computed by the single-
ED
350
objective one-to-all Dijkstra, i.e. all optimal paths in terms of travel time are determined from d to any other node of the grid. The edge weights used in this
PT
computation are the least travel time of edges occurring within the time window [tstart , tstart + T ].
In quest of an efficient algorithm which will optimally solve the ship routing
CE
355
problem, we considered several enhancements of the algorithms in [11]. First, we converted the Continuous time Algorithm of [11] to a time-dependent multi-
AC
objective A∗ algorithm namely, to a time-dependent N AM OA∗ -like algorithm [42], by employing heuristic functions for the objective values of risk and fuel
360
consumption [29]. However, this contributed to a slightly faster algorithm only. Then, we applied the pruning technique with the h function described above to
the Continuous time Algorithm and the new algorithm was considerably faster than the original one [29] in the maritime setting. Since, the Discrete time 17
ACCEPTED MANUSCRIPT
Algorithm in [11] is generally faster than the Continuous time algorithm, we 365
finally decided to use the heuristic function h in the Discrete time algorithm. Notice also that techniques based on the heuristic functions for fuel consumption
CR IP T
and risk cannot be easily incorporated in the Discrete-Time algorithm because its logic differs a lot from that of the a time-dependent N AM OA∗ . The efficiency of the new version of the discrete time Ship Routing Algorithm (Algorithm 1) is 370
evaluated through a number of experiments. Also, both versions of the Discrete time algorithm (the initial and the new) were compared against the backward
algorithm of Hamacher et al’s [28] which was slightly adjusted to the problem at
AN US
hand. Hamacher et al’s algorithm tackles the same problem but with arbitrary arrival time at the destination node, while not taking into account possible 375
constraints on the arrival time at the destination.
Specifically, the algorithms were tested on the same set of randomly generated graphs. The size of the random graphs in the experiments was set to 500, 750, 1000 and 1500 nodes while the number of ingoing and outgoing edges
380
M
was 12, 14, 16, 18 and 20 in different experiments. All the nodes in the graph had the same number of incoming and outgoing edges. The single departure
ED
time instance from the source node tstart was set to 0. The time-dependent edge costs were randomly selected from the set {0, 1, ..., 10} for each time instance t ∈ {0, 1, ...T }. The maximum permissible total travel time T was set equal 385
PT
to 100 and 500. We computed the average execution time over 50 different experiments. The computational results are presented in Table 1.
CE
The experimental results demonstrate that the solution pruning based on the heuristic estimation of the remaining travel time reduces the execution time. In all the experiments, Ship Routing Algorithm 1 outperforms both Discrete time
AC
Algorithm of [11] and Hamacher et al’s algorithm. Specifically, the Discrete
390
time Algorithm of [11] is slower than Ship Routing Algorithm 1, but in some cases namely, in small graphs, appears to be faster. A possible reason for this is that the pruning technique causes more delay than speeding-up when applied on small graphs. Also, the performance of all algorithms are affected by the size of the random graph and its topology, but Ship Routing Algorithm 1 seems 18
ACCEPTED MANUSCRIPT
Table 1: Comperative computational results - departure time tstart = 0, maximum travel time T = 100
1500
12
1
1.251
14
1.24
1.26
16
2.6
2.5
18
5
4.74
20
6.2
5.22
12
1.54
1.6
14
2.24
2.18
16
3.82
18
6.94
20
7.98
12
2.96
14
4.12
16
5.74
18
7.68
20
9.2
12
8.78
8.34
98.56
14
10.46
9.82
123.22
13.22
12.4
141.4
18.60
17.74
179.34
20.38
18.66
200.86
16 18
395
23.42 36.4
50.86 61.36 83.2 42.6
63.44
3.36
78.9
6.12
93.76
7.08
110
2.9
53.24
3.74
69.98
4.98
90.5
6.88
104.76
8.6
135.7
PT
20
Hamacher’s Algorithm
AN US
1000
Algorithm 1
M
750
Discrete time Algorithm of [11]
ED
500
m/n
CR IP T
Average CPU Time (in seconds) n
to be less affected relative to the other two algorithms. Finally, Ship Routing
CE
Algorithm 1 is much faster than Hamacher et al’s algorithm which implies that a forward algorithm can solve this particular problem more efficiently than a
AC
backward one.
400
Admittedly, the use of the heuristic function does not improve the execution
time significantly when the algorithms run on random graphs. However, in the maritime setting, the situation is completely reversed and the use of the heuristic function greatly speeds up the computation (Table 2). Notice that in random graphs the edge costs are accordingly random but in graphs with maritime data 19
ACCEPTED MANUSCRIPT
Table 2: Comparison of Discrete time Algorithm of [11] and Algorithm 1 in maritime setting. The experiment setup is the same as that of Table 8 in Section 6. Discrete time Algorithm of [11]
Algorithm 1
(in hours)
CPU time (in seconds)
CPU time (in seconds)
Piraeus
Paros
3
4
Milos
Nauplio
4.5
20
Piraeus
Serifos
3
12
Chania
Paros
5
29
Piraeus
Syros
3.5
12
Piraeus
Syros
4
35
3 8 5
18 7
13
AN US
Destination
CR IP T
Max Travel Time
Start
the edge costs, especially the time-dependent travel time cost, are correlated 405
with the weather and sea conditions. Thus, in this context, the range of values of passage time cost is relatively short, depending on the weather changes. As a result, the heuristic estimate of the total travel time is commonly a good ap-
M
proximation. In contrast, in random graphs there is no similar guarantee about the heuristic estimate and thus there is a high probability that this estimate will 410
not be close to the real total travel time. Finally it should be mentioned that
ED
when the algorithms run on graphs with maritime data, heavier computation is usually required for processing each label than that in random graphs. Thus, the decrease in the number of labels considered due to the use of the heuristic
when working on random graphs.
CE
415
PT
function has a greater impact on the execution time in maritime setting than
4. Modelization Issues
AC
There are many critical modelling issues that directly affect the performance
of a ship weather routing method. The input grid and the kind of information stored in this structure affect the performance of the algorithm that will be used
420
eventually. It is also important to have an accurate model that predicts the ship response to the weather conditions as well as the impact of these conditions on the actual speed of the ship. A realistic, yet easily implemented model for 20
AN US
Figure 1: Adjacent nodes
CR IP T
ACCEPTED MANUSCRIPT
the fuel consumption of a ship is also needed. Finally, since our main goal is the design of a green ship routing method, the guidelines of the International 425
Maritime Organization (IMO) about voyage safety should also be taken into account.
M
4.1. The Grid
As has been mentioned above, the underlying grid structure is important for
430
ED
the performance of the ship routing algorithm. The grid used in our algorithm consists actually of two different grids namely, the Static and the Dynamic Grid, whose edges correspond to moving directions allowed between grid points. A
PT
lot of different data are spatially stored in the nodes and the edges of these two grids, such as the risk of traversing an edge, its length etc. Furthermore,
CE
the two grids together offer a fast and efficient way to retrieve and summarize 435
important information regarding maritime, geographic and weather data from the area of the Aegean Sea in Greece. In the following, more details about the
AC
two grids are given. 4.1.1. Static Grid
440
The Static Grid is the initial graph structure used for storing the static data.
These data include not only geographic and bathymetric information, but also information about the statutory protected areas of the Aegean Sea. The graph 21
ACCEPTED MANUSCRIPT
is composed of a two-dimensional array. Each array cell corresponds to a specific point on the map. However, not all points are valid; specifically, the following
• the point should be located on the sea,
445
CR IP T
conditions should hold:
• the distance of the point from the shore should be greater than a particular threshold, which is given as an input parameter,
• the point should be located in sufficiently deep waters; the minimum depth is also given as an input parameter,
AN US
• the point should not be located past the country’s maritime borders.
450
If a point satisfies all the conditions above, a new node is created and assigned to corresponding array cell. The so created nodes are connected with each other via bidirectional edges. A node can be connected to at most 16 adjacent nodes with different directions (Figure 1). The coordinates of the grid nodes are such that the distances between successive grid nodes are approximately the same
M
455
and still allowing a ship to move along the edges with constant bearing (rhumb
ED
line navigation). Again, some edges may not be valid and should be removed. Specifically, for an edge to be valid, the conditions mentioned above for points must hold for every point along the edge. It is now clear that for constructing the Static Grid the following parameters must be given as input:
PT
460
• the longitude and latitude of top-left corner point,
CE
• the edge length (for horizontal and vertical edges only),
AC
• the number of rows and columns in the array,
465
• the aforementioned parameters regarding minimum depth and distance from land.
4.1.2. Dynamic Grid The Dynamic Grid enriches the static information from the Static Grid with dynamic data. These data include weather forecast, such as wind direction and 22
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 2: Dynamic Grid
speed, significant wave height, temperature, as well as risk estimations resulting from these meteorological data [43]. For the weather information, both current
ED
470
and future forecasts are retrieved from the National Weather Service. Since, these predictions have a certain time granularity, for time instances lacking
PT
forecast, the weather conditions at each node and edge are determined by interpolating the nearest in time predictions. Also, for each edge, the risk a vessel 475
incurs along the edge given the predicted weather conditions is estimated.
CE
After finishing this processing, the dynamic and the static grid is combined
into a single grid. Figure 2 depicts that grid for some region of the Aegean
AC
Sea, namely the Cyclades Archipelagos. The total number of grid nodes for the whole Aegean Sea is 10, 000, the length of horizontal and vertical edges is 3
480
miles while the maximum distance (long diagonal lines) between two adjacent nodes is equal to 6.9 miles. At this point, it is worth mentioning that while the Static Grid is only produced once for the entire area of the Aegean sea, the Dynamic Grid should 23
ACCEPTED MANUSCRIPT
be rebuild regularly whenever newly updated data are available. For instance, 485
in our application scenario, dynamic grid should be rebuild every three hours
CR IP T
since this is the update rate of weather forecasts for the Aegean sea. 4.2. Control variables
The course of the ship is the only control variable used in our approach.
Specifically, the ship power or equivalently the nominal speed does not change 490
during the voyage, although there are routing methodologies [44], [19], [17], in
which the power and the nominal speed are control variables. The frequent
AN US
change of the ship speed is not recommended, especially in case of coastal navigation. Admittedly, ship power/speed may need to change in order to delay or speed up the passage through a region that currently has or will have ad495
verse weather conditions respectively. However, this can be handled by simply restarting the execution of the proposed algorithm giving as an input new operational parameters. Moreover, changing speed in order to avoid adverse weather
ED
4.3. Ship response
M
and sea conditions is more relevant in transoceanic seafaring.
In order to calculate the navigation time between two points, we have to
500
consider the actual ship speed which is usually lower than the nominal one
PT
due to the added resistances induced by irregular waves and wind during ship navigation. To this end, the ship model described in [4] is used. This model is general, independent of specific ship features. According to this model, the speed reduction depends only on the significant wave height H and the ship-
CE
505
wave relative direction Θ. The actual ship speed VR is calculated using the
AC
following equation: VR = V (H, Θ) = V0 − f (Θ) · H 2
The values of the coefficient f are reported in table 3.
510
If more specific ship characteristics are available, a more realistic model for the vessel response to weather and sea conditions could be used. Specifically,
24
ACCEPTED MANUSCRIPT
o
45 < Θ < 135
o
135o ≤ Θ ≤ 180o
f [kn/ft2 ]
following seas
0.0083
beam seas
0.0165
head seas
0.0248
Table 3: Values of the coefficient f
CR IP T
0o ≤ Θ ≤ 45o
Configuration name
we could use the approximate method established by Kwon [45] for predicting
AN US
the loss of speed caused by the added resistance in irregular waves and wind. 4.4. Fuel consumption
The fuel consumption F measured in kg per hour is calculated using the
515
following equation [9], [46]:
F = sf c · P
M
where sf c and P is the specific fuel consumption and the engine power in BHP (kW) of the ship respectively.
Thus, the total fuel consumption along a route is the product of its travelling
ED
520
time T rT and the fuel consumption per hour F :
PT
F Ctotal = T rT · F
4.5. Voyage safety/risk
CE
Voyage safety is ensured by taking into account the total risk of a route during route calculation as well as the IMO recommendations [1] for avoiding dangerous situations. There may be different ways of defining the risk, but the
AC
exact meaning is hard to be captured. According to the FSA Guidelines (MSC 83/INF.2) [47] risk “is the combination of the frequency and the severity of the consequence”. In line with the FSA guidelines, risk is defined as the probability of an event multiplied by its consequences (14): R=P ×S 25
(14)
ACCEPTED MANUSCRIPT
The proposed model combines the results of a marine accident probability (P ) 525
with the severity of its consequences (S). Concerning hazard identification, marine accidents are considered as poten-
CR IP T
tial hazards. These accidents can be divided into nine categories: foundered, missing vessel, fire/explosion, collision, contact, grounding/wrecked/stranded, war loss/hostilities, hull/machinery damage and miscellaneous. Fault and/or 530
Event tree Analysis, Fuzzy logic and Bayesian networks are the most frequently
used methodologies [48]. In our approach, we employ the model proposed in [43] which takes into account dynamic information about a significant navigation
using a Bayesian network method.
AN US
area as well as the prevailing weather conditions in that area and is developed
Besides the diversity of methods in risk assessment of maritime transporta-
535
tion, there are also many different factors that should be considered when evaluating the probability of an accident. Complex methodologies that take into account weather and environmental conditions, information about the vessel
540
M
and the human factor have been introduced to calculate the accidental risk of shipping in a sea area and identify its causes. In particular, in our approach risk
ED
assessment considers information such as the vessel type, size, age and flag, navigation area characteristics, traffic data, historic data about accidents, weather and sea conditions [48].
545
PT
With respect to the risk of a particular edge, this depends on the weather conditions at the location of the edge and is estimated using the model in [43],
CE
as well. Finally, the risk of a particular route is the total risk encountered along each edge of the route. Regarding IMO recommendations, surf-riding and broaching-to situations
AC
should be avoided when navigating in severe weather conditions. Surf-riding
550
and broaching-to occur when the conditions below are satisfied: √ 1.8 L 135o < Θ < 225o and VR > cos(180 − Θ) where Θ, VR and L are the relative ship-wave angle, the ship speed and the ship length, respectively. 26
ACCEPTED MANUSCRIPT
Since ship speed is constant, the only way to avoid surf-riding and broaching555
to situations is by rejecting the edges where the above conditions are arising. Parametric rolling motions is another precarious situation that should be
CR IP T
avoided. It will occur when the encounter wave period TE is nearly equal to the natural rolling period of ship TR or the encounter period TE is close to one half of the ship roll period TR . The period of encounter TE is calculated by the 560
following formula: TE =
3TW
2 3TW (sec) + VR cos(Θ)
ship speed (in knots). 4.6. Manoeuvring characteristics
AN US
where TW is the wave period, Θ is the relative ship-wave angle and VR is the
In the proposed algorithm, the turning ability of a ship was not considered
565
during path planning. Thus, when a solution path does not comply with the ship
M
turning ability, that path is smoothed by using Cubic Spline interpolation [49]. In particular, smoothing takes place locally between the two consecutive route
570
ED
edges that form the sharp turn and thus the obtained smoothed route will be within the three grid cells that the original route crosses. Therefore, it is reasonable to assume that the smoothed route will not violate any of the conditions
PT
listed in Section 4.1.1. Also, the total risk and the total fuel consumption of the smoothed route are expected to be roughly the same as those in the original route.
CE
If the maximum design rudder angle has to be considered during the cal-
575
culation of Pareto optimal routing paths, the approaches in [14] and [50] could
AC
be easily incorporated in our technique. In [14], a modified network graph is constructed by creating copies of the original grid nodes and each of these copies models a allowable turn at the corresponding original node. On the other hand,
580
the authors in [50] propose a modified Dijkstra algorithm which uses arc labelling for taking into account turn restrictions. Our bi-objective shortest path algorithm can be easily adjusted so that it can work with either of the two 27
ACCEPTED MANUSCRIPT
approaches. However, the execution time would be considerably longer since the modified graph would be much larger than the original one [14] or in the 585
case of [50] because working on arcs would be computationally much heavier
CR IP T
than working on nodes due to their much higher number. Moreover, since sharp turns were seldom observed in our case study, navigation in the Aegean sea, a simple solution with local curve smoothing was preferred eventually.
5. Speeding Up and Improvement Techniques
So far, the ship routing problem has been handled as a bi-criteria, time-
590
AN US
dependent and resource-constrained path problem whose parameters have been
properly set considering the particular application scenario. By closer examination of the common practices and assumptions in this application area, the basic algorithm can be improved so that it can run faster and return more prac595
tical solutions, that is, solutions with higher probability of being adopted by the
M
end-users. In this section, we present speeding up techniques and techniques which improve the practicality of the derived solutions. More specifically, we propose a technique for bounding the solution space around a voyage plan given
600
ED
in advance. We also present a grid partitioning method which can lower the execution time of the algorithm with the caveat that it can be applied only in
PT
the case that the cost functions in the input grid are time-invariant. Moreover, a technique that ensures that there will be no loops in the solution routes is presented. Finally, a selection tool is developed for presenting the most repre-
CE
sentative solutions to the end user when the number of Pareto-optimal solutions
605
is too large for clear illustration.
AC
5.1. Grid Pruning It is a fact that out of experience, shipping companies usually know the
best route between two ports and the voyage plan given to the captain of a ship reflects this knowledge. Since, most companies are rather reluctant to
610
accept large changes in the customary route they follow between two ports, the
28
AN US
Figure 3: Grid pruning
CR IP T
ACCEPTED MANUSCRIPT
proposed routing algorithm should derive solutions that do not much deviate from the original voyage plan of the ship. For this purpose, the original grid is pruned leaving only a “corridor” around the voyage plan (Figure 3). This
615
M
drastically reduces the solution space for finding the optimal solution. Also, the pruning can be easily performed online, during the execution of the routing algorithm. In comparison to the off-line execution, the on-line pruning examines
ED
fewer grid points, saving computation time, as a result. For evaluating the benefits of the grid pruning, we ran the routing algorithm with and without grid pruning. The thickness of the corridor around the voyage plan, distV P , was set to two different distinct values, 42 miles and 21 miles
PT
620
respectively (third line of Table 4). These distance values from the voyage plan
CE
were selected empirically. The tests were carried out on a server with a 4-core 3.6 GHz processor and 8 GB RAM and are listed in Tables 4 and 5. The experiment setup is the same as that of Table 8 in Section 6. In the case of 42 miles distance, the number of the derived solutions is the same as that returned
AC
625
in the original grid and thus we can still find the whole Pareto set. However, if distV P is set to a value lower than 42 miles, some Pareto optimal solutions are not returned (Table 4). Concerning execution time, clearly grid pruning along the voyage plan can reduce the execution time dramatically (Table 5).
29
ACCEPTED MANUSCRIPT
Table 4: Algorithm performance: original and pruned grid
Destination
Number of Routes
distV P
distV P
Max Travel Time (in hours) ∞
42 miles
21 miles
∞
42 miles
21 miles
33
33
29
47
47
47
26
26
26
CR IP T
Departure
CPU time (in seconds)
Paros
3
3
3
3
Milos
Nauplio
4.5
11
8
7
Piraeus
Serifos
3
5
5
5
Chania
Paros
5
36
18
17
50
50
50
Piraeus
Syros
3.5
7
7
7
38
38
36
Piraeus
Syros
4
14
13
13
38
38
36
Spetses
Amorgos
4
53
49
AN US
Piraeus
15
13
13
53
Table 5: Algorithm performance: original and pruned grid (distV P = 42 miles)
Execution time
Destination
M
Start
Original Grid
Pruned Grid
Elafonisi
01:49:42
00:01:09
Kymi
Kalamata
00:01:00
00:00:27
Chios
00:09:43
00:00:59
Cythera
00:04:43
00:00:50
Hydra
Samos
00:00:27
00:00:16
Rhodes
Western Crete
01:17:27
00:05:14
SW Crete
CE
PT
Piperi
ED
Kos
5.2. Dynamically Partitioned Grid
AC
630
In this section, we propose another speed-up technique which can be used
together with the basic ship routing algorithm. However, the main limitation is that it can be applied only for static grids, i.e., with constant costs at the edges and nodes. For short trips, where the weather conditions in the travel area are
635
not likely to change due to short travel time, static grids can be used for solving
30
ACCEPTED MANUSCRIPT
the routing problem. The main objective is the removal of nodes and edges from the input grid in such a way that most of the solutions from the original Pareto set can also be
640
CR IP T
obtained when working on the reduced grid. In addition, the routing algorithm should run faster on the new grid.
This size reduction is feasible since in the Aegean sea, there are broad areas
with few or even no islands where the weather conditions are practically the same. Thus, the objective is to locate convex polygons with same weather
within. Apparently also, a straight line between any two nodes inside such a polygon lies also in that polygon and thus all the points of the line will have the
AN US
645
same weather. Therefore, we can reduce the size of the grid, by maintaining only the boundary nodes of each polygon, termed also as gateways, and connecting them pairwise with straight line edges. The weather and risk condition along each of the new edges can be approximated by averaging the values of these 650
metrics at the end vertices.
M
A previous version of this speeding up technique can be found in [29]. In that approach, we extracted areas with uniform weather and risk conditions
ED
and then applied a convex partitioning algorithm for finding a more accurate partition. However, the number of partitions obtained were relatively high and 655
hence the size of the grid was not reduced a lot. In this work, much smaller grid
PT
is obtained by extracting only the largest convex polygons from each uniform area, maintaining the remaining nodes and edges of the original grid.
CE
The algorithm of creating the Dynamically Partitioned Grid consists of 3 phases.
5.2.1. Phase 1 - Partitioning the Grid into Sectors with Uniform Weather and
AC
660
Risk Conditions The main steps of this phase are given in Algorithm 2. For each node,
depending on the similarity of that node with the neighboring sectors, either it is assigned to an existing adjacent sector or a new sector is created, containing that
665
node. A node and a sector are considered similar if the risk, wind direction and 31
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 4: Overlapping sectors
wave height at the node and within the sector are also similar. For a sector, the
M
values of these parameters are estimated as the mean of the corresponding values at nodes of this sector and all these averages are recalculated upon insertion or
670
ED
removal of a node from the sector. Also, sectors connected diagonally with a node, are not considered unless no sectors exist at the north and west of the node. Otherwise, overlapping areas would be possible (Figure 4).
PT
After the assignment of nodes to sectors has been completed, the similarity of each node with the home sector, as well as with its neighbouring sectors,
CE
is calculated and the node is moved to the highest-scored sector provided that 675
the previous sector still remains connected after this movement. This node
AC
re-affiliation continues until no further change is possible. 5.2.2. Phase 2- Extraction of the largest convex polygon(s) from each sector Fischer’s algorithm [51] is used for extracting the maximum convex polygon
from each sector. This is a O(n4 logn) dynamic programming algorithm taking
680
as input a set of positive and negative points and returns the maximum convex polygon which encloses only positive points. For applying this algorithm, first,
32
ACCEPTED MANUSCRIPT
Algorithm 2 Partitioning the Grid in Sectors with Uniform Conditions Require: The 2-dimensional grid array containing all nodes
CR IP T
Ensure: All nodes with similar risk and weather conditions belong to the same sector 1:
Ni , Wi , NWi , SWi : the northern, western, north western and south western neighbour of node i sci : the sector node i belongs to
sim(arg1, arg2): compares either two sectors or a sector with a node and
AN US
returns a measure of the similarity of the two entities with regard to weather
and risk conditions; the more similar they are the higher the value returned same(arg1, arg2): compares either two sectors or a sector and a node and returns true if the similarity between the two entities is above a certain threshold.
M
merge(sc1, sc2): merges the two sectors given as arguments /* Initialization of sectors
3:
for each node i do/* the grid is parsed from left to right and top to bottom */
5: 6:
add i to scNi
end if
if ∃ Ni , Wi and same(scWi , i) and sim(scWi , i) > sim(scNi , i) then
CE
7:
if ∃ Ni and same(scNi , i) then
PT
4:
add i to scWi
8:
9:
AC
10:
11: 12:
*/
ED
2:
end if if ∃ Ni , Wi and same(scNi , i) and same(scWi , i) and same(scNi , scWi )
then merge(scNi , scWi ) end if
33
13: 14:
if 6 ∃ Ni , Wi then if ∃ N Wi and same(scN Wi , i) then add i to scN Wi
15:
CR IP T
ACCEPTED MANUSCRIPT
end if
17:
if ∃ SWi and same(scSWi , i) and sim(scSWi , i) > sim(scN Wi , i) then add i to scSWi
18: 19:
end if
AN US
16:
20:
end if
21:
if i not in any sector yet then
22:
Create a new sector for i end if
M
23:
end for
25:
/* Associate each node with a more similar local sector if exists
26:
Change = F alse
27:
repeat
29: 30:
if ∃ adjacent sector more similar than the current one then move i to the more similar adjacent sector Change = T rue
CE
31:
for each node i do
PT
28:
ED
24:
32: 33:
AC
34:
end if
end for
until Change = F alse
34
*/
ACCEPTED MANUSCRIPT
the bounding box of each sector sci is determined. Then, any node belonging to sci is considered as a positive point while any other node in the bounding box is treated as a negative point. Also, each sector should be connected and contain 300 nodes at most. Although, these requirements are not necessary for
CR IP T
685
the correctness of the algorithm, they both help in faster execution. Thus, in
case of not connected sector, the sector is split into its connected components. If a sector contains more than 300 nodes, the sector is split along the horizontal
or vertical line that partitions the nodes as evenly as possible and generates the 690
fewest possible gateways.
AN US
However, even for sectors with fewer than 300 nodes, the computation of
maximum convex polygon may not be fast. Thus, for each sector, a sparse sampling of representative positive points of the sector is given as an input to Fischer’s algorithm. Specifically, if the number of nodes in the sector is 695
100 to 200, nodes with odd row and column number are omitted, while for a number between 200 and 300, nodes with either odd row and column number,
M
or even row and column number are ignored. In both cases, nodes lying on the boundaries of the sector are never removed.
700
ED
After this preprocessing, Fischer’s algorithm is executed. If the algorithm is unsuccessful or the resulting polygon has too many gateways compared to the total number of nodes contained in the polygon, the polygon is rejected.
PT
Otherwise, the polygon is stored along with the nodes it contains. A node in a sector is a gateway if at least one of its 8 incident shorter edges is also incident
CE
to a node outside the sector. The 8 long diagonal edges are ignored at this 705
phase, because the gateways should be near the sector perimeter. In addition, we ignore nodes having only one short diagonal edge which is incident with an
AC
external node and also passes between two other gateways. The reason is that these nodes are located inside the sector (Figure 5) and are not part of the polygon perimeter.
710
Now, the remaining nodes of the bounding box not belonging to the polygon
just formed create a new sector and the previous steps are repeated again. This iteration continues until it is not possible to find a maximum convex polygon 35
CR IP T
ACCEPTED MANUSCRIPT
with low number of gateways.
AN US
Figure 5: A case where a node is not considered as gateway.
Finally, for further reducing the execution time, parallel processing is em715
ployed since the computation above is amenable to parallel processing as each sector can be independently processed by a separate thread in a multi-core sys-
M
tem.
5.2.3. Phase 3 - Building the new Grid
720
ED
In this phase, the final grid is obtained by adding or removing some edges from the original grid. Specifically, if a node is not a gateway, all its incident edges are removed. Nevertheless, these isolated nodes are kept, since it might
PT
be chosen as the starting or destination point of a route. In this case, new edges from or to this node are created dynamically, connecting it with the gateways
CE
of its sector. Essentially, the new grid contains only the gateways of the sectors 725
and the nodes that do not belong to any sector. Each pair of gateways belonging to the same sector is connected with a new bidirectional edge while nodes in
AC
different sectors are connected with edges of the original grid, including also the long ones. The costs (risk, fuel, time) of each new edge are calculated as the average of the corresponding costs of the nodes of the original grid met along
730
the edge, gateways inclusive. The derived partitioned grid (Figure 6), contains about 56% less nodes and 24% less edges compared to the original grid. The routing algorithm running on 36
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 6: The partitioned grid.
this grid finds a close approximation of the original Pareto set in considerably
735
M
less time. However, there is a better alternative. After running the routing algorithm on the partitioned grid, we extract the minimum risk and the mini-
ED
mum fuel path as well as any other Pareto-optimal path which is significantly different from the two paths above. Then, these paths are used as a sort of voyage plan and any nodes not near these paths are removed, following again
740
PT
the technique described in Section 5.1. Next, the ship routing algorithm is run on the pruned grid with much shorter execution time. Again, most of the initial
CE
Pareto optimal solutions are extracted, if not the entire set. 5.3. Avoiding Loops It is a well known fact that if waiting at nodes is not an option, the optimal
AC
solutions for the shortest path problem may contain loops [24]. In our appli-
745
cation scenario, this is rather a rare case and it mainly happens when there is an abrupt change in the weather conditions in a local region. Since, the spatial as well as the temporal transition between different weather conditions is rather smooth most of the time, routing loops are not frequently observed in
37
ACCEPTED MANUSCRIPT
the derived routes. Indeed, only in two out of 100 execution runs of the ship 750
routing algorithm on random instances2 , few Pareto-optimal solutions are found to contain loops.
CR IP T
However, a route containing a loop is viewed undesired in ship navigation and it is usually dismissed. A possible approach for deriving loop-less ship routes is to directly modify the Algorithm 1 so that it will maintain only paths 755
that do not form loops at each execution step. However, as has been pointed
out in [23, 24] where special cases of the problem at hand were treated, finding
routes without loops requires an enumeration of all paths between source and
AN US
destination, which increases the state space that should be explored dramatically and hence the execution time. The main reason is that the algorithm should 760
“remember” the intermediate nodes of each path so that it will be easy to check then if an extension of the path is actually feasible, that is, it does not contain nodes previously met in the path. In particular, in our scenario, it is not certain that keeping only Pareto-optimal solutions at each node for each arrival time
765
M
will give the Pareto-optimal routes between the source and the destination, eventually. Thus, non-Pareto-optimal solutions should be kept as well, and this
ED
leads to the complete enumeration of the routing paths mentioned above. Even if we disregard these technical difficulties, the modified Algorithm 1 may derive routes that will not be sought-after again. Specifically, the algorithm
770
PT
may return routes that do not contain “exact” loops but they may pass again through regions very close to locations already traversed. Apparently, this sort
CE
of routes can be also be viewed as routes with routing loops. However, detecting these cases would add additional complexity to the Algorithm 1. For all these
AC
reasons, the option of modifying the basic algorithm was not followed.
775
A more practical approach is to modify the input grid around the voyage
plan so that it will be made acyclic and hence routing loops will not be possible in the returned solutions. This can be done on-line, during the algorithm execution; specifically, all edges of the grid whose angle with the current di2 randomly
chosen source and destination ports and weather conditions
38
ACCEPTED MANUSCRIPT
rection of voyage plan is greater than 90o are discarded. It can be easily seen that the resulting grid after these edge removals becomes acyclic. Applying this 780
technique, optimal solutions may be discarded since grid edges are removed,
CR IP T
however, no noticeable difference in the solution quality was observed in our experiments. More precisely, the proposed technique was tested in 100 different experiments, with random departure and destination nodes. In 94% of the test
cases, the algorithm running on the modified grid retrieved the entire Pareto 785
set which is returned when the algorithm is running without the loop avoidance
technique. In all other test cases, where the whole Pareto set was not returned,
AN US
the quality of the solution routes was only slightly affected, since 92% of the original Pareto set was found again, on average. For instance, in the test cases of Table 4 the algorithm with the loop-avoidance technique found the whole 790
Pareto set in all itineraries, except in the case where the departure and destination ports are the cities of Chania and Paros respectively (line 4 of Table 4); specifically, the algorithm returned 48 solution routes and retrieved 47 out of
M
50 optimal routes that the basic algorithm found. In Figure 7, the solutions of the basic algorithm and those returned by our approach for the voyage Chania to Paros are plotted with blue and red color respectively. Clearly, the difference
ED
795
of the two Pareto-fronts is hardly discernible. Moreover, if the first approach of the direct modification of Algorithm 1 had been followed, the Pareto-front that
PT
would be returned in that case would closely match the Pareto front derived by our approach. This is because the Pareto front derived when the basic algorithm is running without a loop avoidance technique always “dominates” the Pareto-
CE
800
fronts of any loop-avoidance implementation. Since also, the Pareto-front of the directly modified Algorithm 1 would normally be in between the two fronts of
AC
Figure 7, the close matching between the fronts depicted in that figure certainly entails a corresponding matching in the results of our approach and those that
805
would be returned with a direct modification of Algorithm 1. Finally, it is worth mentioning that if the geographical morphology of a
region requires that the route should form an exact or “approximate” loop, the proposed approach can derive this sort of loop provided that the voyage 39
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 7: Chania to Paros: Pareto fronts of the basic algorithm and loop avoidance technique.
plan contain the required cyclic segments. This is possible, since the grid lines 810
following the current direction of the voyage plan are always included in the
5.4. User Selection
M
working grid.
ED
After the execution of the ship routing algorithm, the set of Pareto optimal routes is illustrated on the map (see Figure 12, for instance). However, if the 815
number of these solutions is very large, it is more convenient to illustrate a
PT
smaller number of representative solutions. Specifically, three routes are selected from the Pareto optimal solutions in this case as follows:
CE
1. The first route minimizes the risk and is selected among these solutions where the fuel consumption is no higher than c times the minimum fuel
AC
820
consumption attained over all the Pareto optimal routes where c is an input parameter with a default value between 1.2 and 1.25.
2. The second route minimizes the fuel consumption conditional on the risk of this route being no higher than c times the minimum risk achieved among all the routes of the Pareto set.
40
ACCEPTED MANUSCRIPT
3. The third route minimizes a weighted average of both the risk and fuel
825
consumption among the remaining routes of the Pareto set. Apparently, by changing the value of parameter c, the whole Pareto front
CR IP T
can be gradually returned. Another possibility is to use a measure of route similarity for comparing routes and eliminating routes similar to previously 830
selected routes. Thus, only one route will be kept out of each group of routes following broadly the same trajectory in the final output.
The selected routes are returned to the user, as a sequence of grid points
plotted on a map along with the coordinates of each point. Moreover, an expla-
835
AN US
nation about each route is included in the output. An explanation of this kind
could read as “This route minimizes the risk and has only 17% higher fuel consumption than the most fuel economical route (this route would be 40% more risky)”. Finally, this output is illustrated on a graphical user interface and then the user can choose the most preferable of these routes. For instance, Figure 8
840
M
illustrates the output of the algorithm with Spetses and Amorgos as the departure and destination port respectively. The input voyage plan is illustrated in
ED
Figure 3 and the corridor around the voyage plan was set to 42 miles.
6. Experimental Results
PT
In this section, the experimental evaluation of the baseline algorithm as well as its variants is presented. Firstly, the algorithm based on the dynamically 845
partitioned grid is tested and compared with the case where the algorithm uses
CE
only the original grid as input. Then, the performance of the algorithm enhanced with the loop avoidance and the grid pruning technique is tested in a number
AC
of real-life test instances.
850
In order to test the Dynamically Partitioned Grid, we generated 100 random
test cases. For each test case, a random set of starting and destination ports were chosen and the maximum travel time for each voyage was set proportional to the euclidean distance from the source to the destination port. Since, the technique of Dynamically Partitioned Grid cannot be applied in time-dependent 41
ED
M
AN US
(a) The route minimizing the risk
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
(b) The route minimizing the fuel consumption
(c) The route minimizing the weighted average of risk and fuel cost Figure 8: The three solutions returned to the user (distV P = 42 miles).
42
ACCEPTED MANUSCRIPT
grid networks, the weather and sea conditions do not vary in time but they can 855
be spatially different in these test scenarios. The other test parameters were fixed as follows: speed = 30 kn, s.f.c. = 200 g/kWh, break power = 4000 kWh.
CR IP T
Also, the tests were carried out on a server equipped with a 4-core 3.6 GHz processor and 8 GB RAM.
Some indicative experimental results of these tests are listed in Table 6. It 860
should be mentioned that the baseline algorithm runs on the original grid without grid pruning and loop avoidance technique. Also, in the estimation of the execution time, we did not take into account the time for creating the Dynam-
AN US
ically Partitioned Grid. This is because in our test setting there is no temporal
change in weather conditions and thus the dynamically partitioned grid is cre865
ated only once during initialization. The experimental results suggest that the dynamically partitioned grid technique accelerates the ship routing algorithm, without conspicuously decreasing the quality of the derived solutions. Specifically, the technique of using solution paths from the dynamically partitioned
870
M
grid for creating the voyage plan on the original grid is 2.218 times faster on average than the baseline algorithm; at the same time it is able to generate the
ED
88.32% of the Pareto set approximately. Noteworthily, when the execution time of the baseline algorithm is 5 minutes or more, the algorithm on the dynamically partitioned grid is 4.166 times faster with the greatest improvement being 70-
875
PT
fold. Regarding the solution quality, when using the Dynamically Partitioned Grid, there is an average increase of 0.037% in the risk and 0.051% in the fuel
CE
consumption of the safest and the most fuel economical route respectively in comparison to the original grid. Next, we evaluated the ship routing algorithm when using the grid pruning
AC
(Sect. 5.1) and loop avoidance technique (Sect. 5.3). Besides the standard grid,
880
in these experiments we generated two other grids for the Aegean Sea with coarser and finer resolution respectively (Table 7) for studying the impact of grid resolution on the execution time and the output of the algorithm. For each grid, we generated 100 random test cases, that is, for each case, a random departure and destination port was chosen and the maximum travel time of 43
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
ED
Figure 9: Source-Destination ports of Table 6.
Table 6: Comparison of the algorithm performance on the Original and on the Dynamically
CE
Partitioned Grid
Relative Increase in the Dynam-
Execution Time
ically Partitioned Grid over the Original Grid
Max Travel Time
Percent of Common Routes
Best Risk
Best Fuel Consumption
Original Grid
Dynamically Partitioned Grid
Kos - Elafonisi
10 hrs
82.76%
0%
0%
01:49:42
00:01:33
Kymi - Kalamata
10 hrs
100%
0%
0%
00:01:00
00:00:29
SW Crete - Chios
11 hrs
68.33%
1.64%
0.003%
00:09:43
00:01:14
Piperi - Cythera
10 hrs
100%
0%
0%
00:04:43
00:00:50
Hydra - Samos
8 hrs
83.33%
0%
1.04%
00:00:27
00:00:22
Rhodes - W Crete
13 hrs
99.11%
0%
0%
01:17:27
00:13:10
AC
Route
44
ACCEPTED MANUSCRIPT
Table 7: Grid Information
Coarse
Fine
Grid
Grid
Grid
(SG)
(CG)
(FG)
10,000
6,400
40,000
Distance of Nodes (in miles)
2
2.5
1
Number of Neighbor Nodes
16
16
16
ED
M
AN US
Grid Nodes
CR IP T
Standard
Figure 10: Source-Destination ports of Table 8.
voyage was set again proportional to euclidean distance of the departure and
PT
885
destination port. In all experiments, we set the maximum distance distV P from
CE
the voyage plan equal to 42 miles. Also, in contrast to the previous tests, now the input grids are time-dependent. However, the same weather conditions were assumed for all grids and thus the same time-dependent cost functions were used in all tests concerning each particular grid. Table 8 shows the results of these
AC 890
experiments. Clearly, the proposed algorithm enhanced with the grid pruning and loop
avoidance technique can efficiently solve the point-to-point ship routing problem. The increase in the maximum travel time affects the execution time (lines 5-6,
45
ACCEPTED MANUSCRIPT
895
Table 8) but it does not change the number of the derived solutions. As a side note, since the route planning algorithm is fast when working with a voyage plan, the algorithm can run from scratch and return new Pareto-optimal routes
CR IP T
in case of an unforeseeable situation during the journey. Concerning the impact of grid resolution on the execution time of the algo900
rithm, we can see from Table 8) that the algorithm runs faster on the Coarse
Grid than on Standard Grid except in two cases where the execution time is the same. On other hand, the execution time of the algorithm on the Fine Grid is longer in most tests; it is the same only in one case. The largest increase
905
AN US
in run-time is 16% on the Fine Grid, while the largest decrease is 12% on the
Coarse Grid. This should be expected since the number of edges and vertices of the grid directly affects the number of iterations of the algorithm. However, it should be mentioned again that the algorithm works on a pruned grid around the voyage plan. Otherwise, the impact of the grid resolution on the execution time of the algorithm would be more noticeable since in that case the algorithm search space would be the whole grid and the three grids differ a lot in the
M
910
number of nodes and edges that they contain.
ED
With regard to the solutions returned, the number of the Pareto-optimal solutions is only slightly affected by the resolution of the grid (Table 8). Figure 11 also illustrates the Pareto-fronts obtained for a particular journey, namely the journey from Piraeus Chania. We can see that solutions derived on the Fine
PT
915
Grid are slightly superior and those of the Coarse Grid a little inferior to the
CE
solutions of the Standard Grid. This invariance of the solution quality can be attributed to the use of the pruned grid and the fact that the weather conditions do not abruptly change on a local scale, most of times. Thus, an increase of decrease of detail by using a finer or coarser grid respectively may not be
AC
920
essential due to the same weather over large parts of the sea. Finally, we ran the algorithm with the two enhancements above for finding
the Pareto-optimal routes between two specific ports, namely Piraeus and Chania of the Crete island. In Figure 12, we can see the Pareto front derived by
925
the algorithm, with the maximum distance from the voyage plan set equal to 21 46
ACCEPTED MANUSCRIPT
Departure
Destination
CR IP T
Table 8: Computational results: Execution time
CPU time (in seconds)
Number of Routes
SG
CG
FG
SG
CG
FG
4
33
30
34
9
47
44
49
5
26
26
26
Max Travel Time (in hours)
Paros
3
3
3
Milos
Nauplio
4.5
8
7
Piraeus
Serifos
3
5
4
Chania
Paros
5
Piraeus
Syros
3.5
Piraeus
Syros
4
AN US
Piraeus
16
21
50
45
51
7
6
8
38
36
38
13
13
15
38
36
38
AC
CE
PT
ED
M
18
(a) Standard Grid
(b) Coarse Grid
(c) Fine Grid
Figure 11: Pareto Front from Chania to Paros
47
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 12: Grid pruning and loop avoidance technique: Pareto front for the voyage between
M
Piraeus and Chania.
miles and the maximum duration of the voyage up to 3 hours. The ship routing
ED
algorithm found the whole Pareto set which comprises 48 routes.
PT
7. Conclusion
In this paper, we developed an exact algorithm for solving the time dependent, multi-objective and time constrained ship routing problem with fixed
CE
930
departure time. The final routes are optimal with respect to two different conflicting objectives namely, the fuel consumption and the risk along the route.
AC
In order to enhance the safety of the voyage, we have also considered the IMO regulations. Due to heavy computational demands, we focused on reducing the
935
execution time by reforming the grid structure. We have also enhanced the quality of the solution routes by using techniques which ensure that no loops will exist in the derived solutions and also increase the similarity of the solution routes to a given voyage plan. The experimental evaluation of the proposed ap48
ACCEPTED MANUSCRIPT
proach confirms that our techniques efficiently solve the ship-routing problem 940
and they are suitable for real-time applications.
CR IP T
References [1] M. IMO, 1/circ. 1228, Revised guidance to the master for avoiding dangerous situations in adverse weather and sea conditions, adopted 11th January.
[2] C. P. Padhy, D. Sen, P. K. Bhaskaran, Application of wave model for weather routing of ships in the north indian ocean, Natural Hazards 44 (3)
945
AN US
(2008) 373–385.
[3] R. Szlapczynski, A new method of ship routing on raster grids, with turn penalties and collision avoidance, Journal of Navigation 59 (01) (2006) 27– 42. 950
[4] G. Mannarini, G. Coppini, P. Oddo, N. Pinardi, A prototype of ship
M
routing decision support system for an operational oceanographic service, TransNav, the International Journal on Marine Navigation and Safety of
ED
Sea Transportation 7 (1) (2013) 53–59. [5] W. Shao, P. Zhou, S. K. Thong, Development of a novel forward dynamic programming method for weather routing, Journal of marine science and
955
PT
technology 17 (2) (2012) 239–251. [6] A. Azaron, F. Kianfar, Dynamic shortest path in stochastic dynamic net-
CE
works: Ship routing problem, European Journal of Operational Research 144 (1) (2003) 138–156.
[7] H. K. Lo, M. R. McCord, Adaptive ship routing through stochastic ocean
AC
960
currents: General formulations and empirical results, Transportation Research Part A: Policy and Practice 32 (7) (1998) 547–561.
[8] A. A. Montes, Network shortest path application for optimum track ship routing, Tech. rep., DTIC Document (2005).
49
ACCEPTED MANUSCRIPT
965
[9] K. Takashima, B. Mezaoui, R. Shoji, On the fuel saving operation for coastal merchant ships using weather routing, in: Proceedings of Int. Symp. TransNav, Vol. 9, 2009, pp. 431–436.
CR IP T
[10] I. S. Dolinskaya, Optimal path finding in direction, location, and time dependent environments, Naval Research Logistics (NRL) 59 (5) (2012) 325–339.
970
[11] A. Veneti, C. Konstantopoulos, G. Pantziou, Continuous and discrete time
label setting algorithms for the time dependent bi-criteria shortest path
AN US
problem, in: Operations Research and Computing: Algorithms and Software for Analytics, 2015, pp. 62–73. 975
´ Saux, C. Claramunt, A bidirectional path-finding algorithm [12] D. Tsatcha, E. and data structure for maritime routing, International Journal of Geographical Information Science (ahead-of-print) (2014) 1–23.
M
[13] L. Babel, T. Zimmermann, Planning safe navigation routes through mined waters, European Journal of Operational Research. [14] I. Ari, V. Aksakalli, V. Aydogdu, S. Kum, Optimal ship navigation with
ED
980
safety distance and realistic turn constraints, European Journal of Operational Research 229 (3) (2013) 707–717.
PT
[15] O. Kosmas, D. Vlachos, Simulated annealing for optimal ship routing, Computers & Operations Research 39 (3) (2012) 576–581. [16] S. Harries, J. HEIMANN, J. Hinnenthal, Pareto optimal routing of ships,
CE
985
transformation 772 (2003) 1.
AC
[17] S. Marie, E. Courteille, et al., Multi-objective optimization of motor vessel
990
route, in: Proceedings of the Int. Symp. TransNav, Vol. 9, 2009, pp. 411– 418.
[18] A. Dec` o, D. M. Frangopol, Real-time risk of ship structures integrating structural health monitoring data: Application to multi-objective optimal ship routing, Ocean Engineering 96 (2015) 312–329. 50
ACCEPTED MANUSCRIPT
[19] J. Szlapczynska, R. Smierzchalski, Multicriteria optimisation in weather routing, Marine Navigation and Safety of Sea Transportation (2009) 423. 995
[20] M.-C. Tsou, Integration of a geographic information system and evolution-
Navigation 63 (02) (2010) 323–341.
CR IP T
ary computation for automatic routing in coastal navigation, Journal of
[21] S. Bijlsma, A computational method in ship routing using the concept of limited manoeuvrability, Journal of Navigation 57 (03) (2004) 357–369. 1000
[22] M.-C. Fang, Y.-H. Lin, The optimization of ship weather-routing algorithm
AN US
based on the composite influence of multi-dynamic elements (ii): Optimized routings, Applied Ocean Research 50 (2015) 130–140.
[23] A. Orda, R. Rom, Minimum weight paths in time-dependent networks, Networks 21 (3) (1991) 295–319. 1005
[24] A. Orda, R. Rom, Shortest-path and minimum-delay algorithms in net-
ED
37 (3) (1990) 607–625.
M
works with time-dependent edge-length, Journal of the ACM (JACM)
[25] D. E. Kaufman, R. L. Smith, Fastest paths in time-dependent networks for intelligent vehicle-highway systems application, Journal of Intelligent Transportation Systems 1 (1) (1993) 1–11.
PT
1010
[26] M. M. Kostreva, M. M. Wiecek, Time dependency in multiple objective
CE
dynamic programming, Journal of Mathematical Analysis and Applications 173 (1) (1993) 289–307.
AC
[27] T. Getachew, M. Kostreva, L. Lancaster, A generalization of dynamic
1015
programming for pareto optimization in dynamic networks, RAIROOperations Research 34 (01) (2000) 27–47.
[28] H. W. Hamacher, S. Ruzika, S. A. Tjandra, Algorithms for time-dependent bicriteria shortest path problems, Discrete optimization 3 (3) (2006) 238– 254. 51
ACCEPTED MANUSCRIPT
1020
[29] A. Makrygiorgos, I. A. Vetsikas, S. Perantonis, Accelerating multi-objective ship routing using a novel grid structure and a simple heuristic, in: Proceedings of the 6th International Conference on Information, Intelligence,
CR IP T
Systems and Applications, IEEE, 2015. [30] N. Touati, V. Jost, On green routing and scheduling problem, arXiv preprint 1203.1604.
1025
[31] M. Christiansen, K. Fagerholt, D. Ronen, Ship routing and scheduling: Status and perspectives, Transportation science 38 (1) (2004) 1–18.
AN US
[32] S. Bijlsma, A computational method for the solution of optimal control problems in ship routing, Navigation 48 (3) (2001) 145–154. 1030
[33] O. Bolza, Vorlesungen u ¨ber Variationsrechnung, BoD–Books on Demand, 2013.
M
[34] R. W. James, Application of wave forecasts to marine navigation. [35] K.-Y. Chang, G. E. Jan, I. Parberry, A method for searching optimal routes
ED
with collision avoidance on raster charts, The Journal of Navigation 56 (03) (2003) 371–384.
1035
[36] S. E. Dreyfus, An appraisal of some shortest-path algorithms, Operations
PT
research 17 (3) (1969) 395–412. [37] A. Warner-Grieve, Optimal route finding with collision avoidance in marine
CE
navigation, Tech. rep. (2006).
[38] A. Veneti, C. Konstantopoulos, G. Pantziou, An evolutionary approach to
AC
1040
multi-objective ship weather routing, in: Information, Intelligence, Systems and Applications (IISA), 2015 6th International Conference on, IEEE, 2015, pp. 1–6.
[39] A. N. Perakis, N. A. Papadakis, Minimal time vessel routing in a time-
1045
dependent environment, Transportation Science 23 (4) (1989) 266–276.
52
ACCEPTED MANUSCRIPT
[40] H. Lee, G. Kong, S. Kim, C. Kim, J. Lee, Optimum ship routing and its implementation on the web, in: Advanced Internet Services and Applications, Springer, 2002, pp. 125–136.
CR IP T
[41] J. Bisschop, AIMMS-Optimization modeling, Paragon Decision Technology, 2011.
1050
[42] L. Mandow, J. L. P. De La Cruz, Multiobjective a* search with consistent heuristics, Journal of the ACM (JACM) 57 (5) (2010) 27.
[43] I. Koromila, Z. Nivolianitou, T. Giannakopoulos, Bayesian network to pre-
AN US
dict environmental risk of a possible ship accident, in: Proceedings of the
7th International Conference on PErvasive Technologies Related to Assis-
1055
tive Environments, ACM, 2014, p. 44.
[44] S. Wei, P. Zhou, 23. development of a 3d dynamic programming method for weather routing, Methods and Algorithms in Navigation: Marine Nav-
1060
M
igation and Safety of Sea Transportation (2011) 181.
[45] Y. Kwon, Speed loss due to added resistance in wind and waves, Naval
ED
Architect (2008) 14–16.
[46] K. Avgouleas, Optimal ship routing, Ph.D. thesis, Massachusetts Institute
PT
of Technology (2008).
[47] I. MSC, 83/inf. 2: Formal safety assessmentconsolidated text of the guidelines for formal safety assessment (fsa) for use in the imo rule-making pro-
CE
1065
cess (msc/circ. 1023- mepc/circ. 392), International Maritime Organisation:
AC
London, UK.
[48] Z. S. Nivolianitou, I. A. Koromila, T. Giannakopoulos, Bayesian network to
1070
predict environmental risk of a possible ship accident, International Journal of Risk Assessment and Management 19 (3) (2016) 228–239.
[49] C. De Boor, A practical guide to splines, Vol. 27, Springer-Verlag New York, 1978. 53
ACCEPTED MANUSCRIPT
[50] E. Guti´errez, A. L. Medaglia, Labeling algorithm for the shortest path problem with turn prohibitions with application to large-scale road networks, Annals of Operations Research 157 (1) (2008) 169–182.
1075
AC
CE
PT
ED
M
AN US
putation Theory, Springer, 1993, pp. 234–243.
CR IP T
[51] P. Fischer, Finding maximum convex polygons, in: Fundamentals of Com-
54
Minimizing the fuel consumption and the risk in maritime transportation: a bi-objective weather routing approach Aphrodite Veneti, Angelos Makrygiorgos, Charalampos Konstantopoulos, Grammati Pantziou, Ioannis Vetsikas PII: DOI: Reference:
S0305-0548(17)30193-4 10.1016/j.cor.2017.07.010 CAOR 4291
To appear in:
Computers and Operations Research
Received date: Revised date: Accepted date:
27 September 2016 6 July 2017 13 July 2017
Please cite this article as: Aphrodite Veneti, Angelos Makrygiorgos, Charalampos Konstantopoulos, Grammati Pantziou, Ioannis Vetsikas, Minimizing the fuel consumption and the risk in maritime transportation: a bi-objective weather routing approach, Computers and Operations Research (2017), doi: 10.1016/j.cor.2017.07.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • A time-dependent algorithm is presented for the ship weather routing
CR IP T
problem. • Given a maximum travel duration constraint, fuel consumption and risk are minimized.
• Improvements of the approach are proposed for faster execution but similar solutions.
AN US
• Techniques are also presented that improve the practicality of the derived
AC
CE
PT
ED
M
solutions.
1
ACCEPTED MANUSCRIPT
Minimizing the fuel consumption and the risk in maritime transportation: a bi-objective weather routing approach
CR IP T
Aphrodite Venetia,∗, Angelos Makrygiorgosc , Charalampos Konstantopoulosa , Grammati Pantzioub , Ioannis Vetsikasc a Department
of Informatics, University of Piraeus, Greece of Informatics, Technological Educational Institute of Athens, Greece c Institute of Informatics & Telecommunications, NCSR Demokritos, Agia Paraskevi, Greece
AN US
b Department
Abstract
The paper presents an improved solution to the ship weather routing problem based on an exact time-dependent bi-objective shortest path algorithm. The two objectives of the problem are the minimization of the fuel consumption and the total risk of the ship route while taking into account the time-varying sea and
M
weather conditions and an upper bound on the total passage time of the route. Safety is also considered by applying the guidelines of the International Maritime
ED
Organization (IMO). As a case study, the proposed algorithm is applied for finding ship routes in the area of the Aegean Sea, Greece. Enhancements of the proposed algorithm are also presented which improve the efficiency of our
PT
approach.
Keywords: Ship weather routing, Multi-criteria optimization, Time
CE
dependent networks, Resource-constrained shortest path.
AC
1. Introduction Maritime transportation is crucial for the economic growth of many coun-
tries as well as for the global economy, in general. During the past two decades, ∗ Corresponding
author Email address: [email protected] (Aphrodite Veneti)
Preprint submitted to Computers & Operations Research
July 19, 2017
ACCEPTED MANUSCRIPT
there has been a growing interest in methods for calculating optimal ship tra5
jectories. The traditional economic objective was to minimize the total traveled distance while later, additional energy efficiency objectives were added such as
CR IP T
the minimization of the fuel consumption and gas emissions. Besides energy efficiency, navigation safety is also considered since sea transportation poses a high environmental risk due to the possibility of maritime accidents, especially if haz10
ardous materials (i.e. chemicals or oil) are transferred. Moreover, ship routing methods have to adhere to the safety guidelines specified by the International Maritime Organization [1].
AN US
The problem of optimal route planning in maritime transportation (or ship
routing problem) takes into consideration the different objectives and constraints 15
set by the ship owners, the national regulations, the international organizations, etc. Regardless of the specific objectives/constraints, the main factor that makes the ship routing problem difficult to solve, is the time-varying weather conditions [2]. The parameters and cost factors of the problem as well as the optimization
20
M
criteria are strongly affected by the weather and sea conditions and therefore, finding an optimal ship route should take into consideration the weather fore-
ED
casts and sea conditions during the course of the route. To accomplish this, a reliable model that predicts the ship’s response to the weather and sea conditions is required, while an effective solution approach for the ship routing problem
25
PT
should determine the heading control and ideally the optimal ship power settings according to the weather and sea conditions. Thus, broadly speaking, the
CE
ship routing problem is a dynamic, time-dependent problem which is usually referred to as the ship weather routing problem. The optimization objectives in the ship routing problem are usually the
AC
minimization of the voyage time, fuel consumption and voyage risk. The ship
30
routing algorithms appeared in the literature are classified into two categories, namely the exact and the heuristic approaches. The exact algorithms [3, 4, 5, 6, 7, 8, 9, 10, 2, 11, 12, 13, 14] determine the optimal solution at the expense of the execution time while heuristic approaches [15, 16, 17, 18, 19, 20] run faster but search only inside a subspace of the search space for the best solution. 3
ACCEPTED MANUSCRIPT
35
Also, for coastal navigation and trans-oceanic seafaring, the methodology of the calculus of variations [21] has also been employed and is based on long-term weather forecasts. Finally, a classic method for optimizing ship routes is still
CR IP T
the isochrone method [22]. In this paper, we consider the bi-objective time-constrained ship weather 40
routing problem as a time-dependent bi-objective point-to-point shortest path
problem with a constraint on the total voyage time. The two objectives of the problem are the minimization of the fuel consumption and the total risk. In
practice, frequent change of ship power or equivalently of nominal speed during
45
AN US
the voyage is avoided [3] especially in the case of coastal navigation. Thus, in
this work, the navigation speed is not a control variable; it is assumed to be constant during the whole journey and should be given as an input parameter. Notably, when the objective is the fastest route between two ports, the networks modelling the sea transportation commonly have the FIFO property1 . However, in our problem setting, the objectives are different and leaving from a node immediately after the arrival may not be the best option, for instance, when
M
50
the costs (fuel consumption and risk) along the outgoing link are decreasing
ED
in the next time interval [23]. Although, computing shortest paths in FIFO networks is a polynomially solvable problem [24, 25], this is not the case when optimal waiting policy should be determined for other cost objectives [23]. On the other hand, frequent stops in the middle of the sea or alternating ship speed
PT
55
frequently is not a common practice in ship routing [3]. Thus, although it is
CE
generally better to wait at a node, in our maritime setting, we assume that there is neither the option of waiting nor the option of speed decrease; by decreasing the speed, we could simulate the forbidden waiting and shift the arrival at the node later exactly at the optimal departure time.
AC
60
To sum up, in this work, we address the time-dependent bi-objective shortest 1A
network has the FIFO property if all of its arcs have that property i.e., for each arc
(i, j) of the network, earlier departure from i always leads to earlier arrival at j. Therefore, the arrival events at j are in the same chronological order as the departure events at i.
4
ACCEPTED MANUSCRIPT
path problem on a network with objectives other than the travel time, with fixed departure time, no waiting at nodes, constant nominal speed and a constraint on the total travel time. Also, since the problem is time-dependent, it is assumed that all arc costs are deterministic functions of the time. However, considering
CR IP T
65
the low execution time of the proposed techniques, dynamic situations where the
weather changes unpredictably during the journey can be handled by running the algorithms afresh.
We first give a non linear integer programming formulation of the problem, 70
and in the sequel we present a new exact algorithm for solving it. To the best
AN US
of our knowledge, only a few algorithmic approaches for the multi-objective
time-dependent shortest path problem have been proposed in the literature [26],[27],[28] are not efficient when waiting at nodes is forbidden. We employ a forward label-setting approach to efficiently tackle the aforementioned time75
dependent bi-objective shortest path problem. The number of the labels being processed during the algorithm execution is kept at a reasonable level, reducing
M
the computational overhead. The algorithms is compared against the backward algorithm of Hamacher et al’s [28], the best so far algorithm in the literature
80
ED
that solves the problem at hand, and the experimental results confirm the faster execution time of the proposed algorithms. A preliminary version of the algorithm appears in [11].
PT
Furthermore, we present critical modelling issues that directly affect the performance of our problem solution approach. We introduce a novel grid structure
CE
that is used as a base layer for solving the ship weather routing problem and 85
significantly improves the efficiency of the proposed algorithm. We also proposed a number of methods to further improve the efficiency of our approach,
AC
namely: (i) a dynamically partitioned grid for reduced graph size as well as an alternative reduction technique where grid is pruned along a given voyage plan, (ii) a heuristic function used to transform the initial algorithm into a bi-
90
objective A* algorithm thereby reducing the search space and (iii) a technique for deriving loopless ship routes, thus making the solutions conformable to the usual practice in this industry. In order to test and validate our approach in a 5
ACCEPTED MANUSCRIPT
realistic setting, we applied the proposed algorithm as well as the speeding up techniques and improvements on the algorithm for the optimal ship route plan95
ing in the Aegean Sea, Greece. The experimental results confirm the efficiency
for some of the afrementioned techniques appear in [29].
CR IP T
of the approach and its applicability in a realistic scenario. Preliminary results
The rest of the paper is organized as follows. Section 2 overviews the related work. Section 3 presents a non-linear integer programming formulation 100
of the ship weather routing problem which is considered as a time-dependent bi-objective point-to-point shortest path problem, and gives an efficient exact
AN US
algorithm for solving the problem. Section 4 presents critical modeling issues while Section 5 includes the speeding up techniques and improvements on the
proposed algorithm. Section 6 presents the empirical evaluation of the proposed 105
algorithm and its enhancements. Finally, Section 7 concludes our work.
M
2. Related Work
The ship routing problem is a multi-objective, non-linear optimization problem with constraints which, due to the importance of its applications, has been
110
ED
widely studied in the literature. In [30] the common objectives of the optimization problems in the transportation industry are classified into three main
PT
categories: classical/economic objectives, environmental (ecosystem)/regional fairness/health objectives, and climate and sustainability objectives. Until recently, in ship navigation the traditional economic objective was to minimize
CE
the total distance travelled [31]. Other relevant objectives were the minimiza-
115
tion of the fuel consumption and gas emissions. However, both these metrics are very complex to estimate. Nowadays, beyond energy efficiency, navigation
AC
safety is also a main concern, since sea transport is the main carrier of crude oil and other hazardous materials. Furthermore, the proposed methodologies have to adhere to the safety guidelines specified by the International Maritime
120
Organization. Exact as well as heuristic approaches have been proposed in the literature for
6
ACCEPTED MANUSCRIPT
ship route planning. Exact algorithms derive optimal solutions of the problem, however, they usually have heavy computation demands. In contrast, a heuristic approach searches for a feasible solution only within a subspace of the search space hopefully, close to the optimal solution. Hence the execution time of a
CR IP T
125
heuristic algorithm is generally much lower than that of an exact algorithm. Since long-term weather forecasts are available, calculus of variations methods
have also been proved suitable for both coastal navigation and trans-oceanic
seafaring [32, 21, 33]. Another old but popular method for optimizing ship 130
route is the isochrone method [34, 22].
AN US
Dynamic Programming and Label Setting methods are commonly met in exact approaches for solving the single objective ship routing problem with or without constraints. In [8] the Optimum Track Ship Routing (OTSR) system is presented, an automation of the weather ship routing service provided by the 135
US navy. A binary heap version of Dijkstra’s algorithm determines the optimum route taking into consideration the weather and sea conditions. Safety of navi-
M
gation is ensured by eliminating grid points where wind speed and wave height exceed ship’s predefined limits. Based on dynamic programming, the authors
140
ED
in [3] presented a modified version of the maze routing method first presented in [35]. The main objective is the minimization of the travel time while turn penalties and collision avoidance techniques are considered for enhancing the
PT
navigation safety. The proposed approach uses relatively few control variables and is space and time efficient. However, it does not take into account the exist-
CE
ing weather and sea conditions. In [9] a method based on Dijkstra’s algorithm is 145
proposed for finding routes with minimum fuel consumption. The input grid is built starting from a given voyage plan and then additional nodes are added on
AC
lines perpendicular to the input route. The authors in [10] incorporate Dreyfus approach [36] so as to solve the ship routing problem in a dynamic environment using as input, measurements of the ocean wave-field surrounding the vessel.
150
In [5], a forward 3-dimensional dynamic programming method was employed where both ship speed and ship course are control variables. The objective is to minimize the fuel consumption while taking into account geographical, 7
ACCEPTED MANUSCRIPT
weather, operational and safety constraints. It was shown that optimizing only the ship course is not as effective as optimizing both the course and the speed. 155
In [4] a Decision Support System is presented which uses meteo-marine and
CR IP T
oceanographic operational information data for all relevant environmental field variables (wave height, direction, period and currents). It uses a modification
of Dijkstra’s algorithm in order to find the optimal route with respect to the
total navigation time in the presence of time-varying fields. Safety restrictions 160
are also taken into account.
Another exact method with fast convergence for solving the ship routing
AN US
problem is the A∗ search algorithm. In [37] a two-tiered A∗ approach is employed
where firstly, the A∗ algorithm is applied to a macro data map to find a nondetailed path from the origin to the destination, and then the A∗ algorithm 165
is applied to obtain an optimal path for every edge of the suboptimal path computed in the first phase, using micro map data. The authors in [14] run the A∗ algorithm on a grid with eight neighbors per node for finding paths
M
optimized in terms of the navigation time subject to safety distance and turnradius constraints and in the presence of obstacles. The heuristic information of A∗ is the shortest possible travel time between a node and the destination. As
ED
170
the path given by the A∗ algorithm is piece-wise linear, the algorithm applies a smoothing technique on this path so that the final route is a realistic route that
PT
can be followed by ships in practice. In [13] the A∗ algorithm is employed to tackle the route planning problem in the presence of obstacles. The objective is varying depending on the particular scenario at hand (e.g. the minimization
CE
175
of the route length, the number of manoeuvres, etc.). The authors adapt the A∗ algorithm to find the least cost route using the exact distance between the
AC
nodes as the heuristic function. The risk of the path is also taken into account.
180
The main heuristic methods that have been employed in the literature for
solving the ship routing problem, are the evolutionary computation and the simulated annealing method. The fast convergence of these methods to a near optimal solution makes them suitable for solving the ship routing problem in real-time. Despite the potential of these heuristic methods in solving multi8
ACCEPTED MANUSCRIPT
objective problems, they were mainly used for solving the single objective ship 185
routing problem and for minimizing a weighted sum of multiple different objectives [15]. In the sequel we briefly refer to the research works which employ
CR IP T
heuristic methods for solving the multi-objective ship weather routing problem. In [19] a heuristic method is presented for multi-criteria ship weather routing optimization with respect to the voyage time, the fuel consumption, and the 190
voyage risk. The approach is based on the Strength Pareto Evolutionary Algo-
rithm (SPEA), which was designed for combinatorial and continuous function
multiple objective optimization problem instances. The solution routes are de-
AN US
fined as a sequence of waypoints from the origin to the destination. The initial population consists of a set of basic routes: an orthodrome, a loxodrome, a 195
time-optimized isochrone route and a route given by applying fuel-optimization to the time-optimized isochrone route. In [17] a Genetic Algorithm is employed for retrieving a Pareto set of optimal ship routes with respect to the voyage time and the fuel consumption. The initial population is chosen as a random sequence
200
M
of grid nodes. The authors take into account the time dependency of the sea state and wind field applying linear interpolation in time when measurements
ED
are missing. In [20] an evolutionary algorithm is proposed for obtaining a number of low cost routes which are also of high safety and avoid narrow waters, bad weather, foggy areas, congested and fishing areas. The initial candidate
205
PT
route population is automatically generated by GIS. In [38] a genetic algorithm for solving the bi-objective time-constrained ship weather routing problem is
CE
presented. The derived routes are optimized with respect to two conflicting objectives, the total fuel consumption and the mean risk of the route while taking
AC
into account the weather conditions.
210
Route planning in the maritime setting is a complex problem due to a number
of dynamic phenomena and navigation constraints. Time-dependent algorithms can efficiently handle the inherently time-varying environment of ship routing. Despite the importance of the time-dependent shortest path problem and its applications in maritime navigation, there are relatively few works in the literature on these issues ([39, 17, 19, 5, 10, 11]). In contrast, most literature on 9
ACCEPTED MANUSCRIPT
215
ship routing and scheduling focuses on the time-invariant shortest path problem [40, 21, 8, 3, 2, 9, 37, 14, 13]. However, there are also works which consider the inherent uncertainty in the behaviour of ocean currents and weather conditions
to minimize the expected fuel consumption [6, 7].
CR IP T
and treat ship routing as a stochastic optimization problem whose objective is
In this paper, we propose a new efficient exact algorithm for the bi-objective
220
time-constrained ship weather routing problem. To the best of our knowledge, it is the first time that an exact multi-objective time-dependent algorithm is employed for solving the ship weather routing problem. Apart from finding op-
225
AN US
timal routes with respect to economical objectives, we are interested in reducing the possibility of an accident taking also care of marine resources conservation. The algorithm returns the whole Pareto set and the derived routes are opti-
mized with respect to two conflicting objectives, the total fuel consumption and the mean risk of the route while taking into account the prevailing weather
230
M
conditions.
3. The bi-objective time-constrained ship weather routing problem
ED
As already mentioned, we consider the bi-objective time-constrained ship weather routing problem as a time-dependent bi-objective time-constrained short-
PT
est path problem. We deal with the most realistic point-to-point instance of the problem, where the ship speed is not a control variable and waiting at nodes is 235
not permitted. In this section, we first give a non linear integer programming
CE
formulation of the problem, and then we present a new forward label-setting time-dependent bi-objective shortest path algorithm to solve the problem. Let G = (V, A) be a directed graph, where V is the set of nodes (|V | = n)
AC
and A is the set of edges (|A| = m). In the bi-objective time-constrained ship
240
weather routing problem, there are three time-dependent positive costs c1ij (t), c2ij (t) and c3ij (t) associated with each edge (i, j) ∈ A, namely, the travel time, the fuel consumption and the risk respectively, incurred along this edge when leaving from node i at time t. We follow the frozen arc model of [23], where
10
ACCEPTED MANUSCRIPT
the above edge costs are assumed to be constant during the traversal of the 245
edge. Let also tstart denote the departure time from source node s ∈ V , and T rT (P ), F C(P ), R(P ) denote the total travel time, fuel consumption and risk
CR IP T
respectively, accumulated along a path P . The problem is to find a shortest path P between a source s ∈ V and a
destination d ∈ V that minimizes the two objectives F C(P ) and R(P ), leaving 250
from the source node at time tstart , while not violating the constraint T rT (P ) ≤
T , where T is the maximum acceptable duration of a route. In the following
formulation, we assume an upper bound L on the number of edges a solution
AN US
path may contain. This constraint is necessary since in the case of forbidden waiting at nodes, the optimal solution may contain infinite number of edges 255
at least in principle, as a result of infinite loops created along the route in an attempt to simulate waiting at nodes which is not permitted [23]. Thus, actually, the derived solution may be a walk and not necessarily a path. Since, this potentially infinite number of edges cannot be expressed by the finite number
260
M
of variables in the above network flow problem, an additional constraint on the number of edges in a solution path should be placed.
ED
We have used a number of variables in the formulation below. Specifically, xi ∈ {0, 1} with xi = 1 when the solution path has exactly i edges. We also
define the binary variable erij with erij = 1 when the edge eij is the r-th edge
265
PT
starting from s along the route. We also define trj as the arrival time at node j through the edge eij when erij = 1. If such an edge does not exist, trj is set to a
AC
CE
large value. We also set t0s = tstart .
11
ACCEPTED MANUSCRIPT
r=1...L (i,j)∈A
X
R(P ) =
X
c3ij (tr−1 )erij i
r=1...L (i,j)∈A
X
T rT (P ) =
X
r=1...L (i,j)∈A
X
c1ij (tir−1 )erij ≤ T
e1sj = 1
AN US
(i,j)∈A
erij ≤ 1, r = 2 . . . L
X
erdj = 0, r = 2 . . . L
j∈V −{d}
X
xi = 1
i=1..L
erjd = xr , r = 1 . . . L
M
X
j∈V −{d}
X
erji −
X
j∈V
er+1 = 0, r = 1 . . . L, i ∈ V − {d} ij
ED
j∈V
(2)
trj =
X
PT CE AC
X
(4)
(6) (7) (8) (9) (10)
erij (tr−1 + c1ij (tr−1 )) i i
i∈{k|(k,j)∈A}
+ (1 −
(3)
(5)
j∈V −{s}
X
(1)
CR IP T
min z = (F C(P ), R(P )) X X F C(P ) = c2ij (tr−1 )erij i
i∈{k|(k,j)∈A}
(11) erij )M, j ∈ V, r = 1 . . . L
t0s = tstart
(12)
cqij (M )
(13)
= M0
erij , xr ∈ {0, 1}, tri ∈ R+ , i, j ∈ V, r = 1 . . . L
Constraint (4) requires that the total travel time of the solution path must
not exceed the maximum acceptable duration T . Constraint (5) states that the first edge of the solution path is an outgoing edge of the start node s, while
270
constraint (6) states that at most one of the edges of the graph can be the r-th edge of the solution path. Constraint 7 ensures that the solution path cannot 12
ACCEPTED MANUSCRIPT
extend beyond the destination node d regardless of its actual length. Constraint (8) determines the length of the solution path. Constraint 9 ensures that the solution path will end up at the destination node d after traversing the right number of edges. Constraint (10) is the ordinary flow conservation constraint
CR IP T
275
ensuring that there is only one incoming and one outgoing edge at each node
along the solution path and also that these edges are successive in the path. Equation (11) gives the arrival time at node j when this node is the head of the
r-th edge of the solution path. In this equation, we assume that M is a large 280
number and also that cqij (M ) = M 0 (q = 1 · · · 3) where M 0 is a large number
AN US
again.
Each solution z contains the values of the two objective functions for the corresponding path. Also, although the proposed Integer Programming model is non linear, when the time dependent edge costs functions c1ij (t), c2ij (t) and 285
c3ij (t) are linear, the above formulation can be converted into a linear model [41]. Since the objective functions may be conflicting, a single solution that simul-
M
taneously optimizes each objective may not exist and therefore, the problem is actually to find the set of non dominated solutions (Pareto optimal solutions).
290
ED
A non dominated solution is the solution where none of the objectives can be improved without worsening the others. Formally, a path p is said to dominate another path p0 (p ≺ p0 ) if (F C(p) < F C(p0 ) and R(p) ≤ R(p0 )) or
PT
(F C(p) ≤ F C(p0 ) and R(p) < R(p0 )).
In the sequel, we present a new forward label-setting algorithm (Algorithm
CE
1) for solving the time-dependent bi-objective shortest path problem with no 295
waiting at the nodes and finding the set of Pareto optimal paths. Noteworthily, although we handle the bi-objective case, the proposed algorithm turns out to
AC
be applicable in the case of multiple objectives, as well. The Ship Routing Algorithm 1 keeps a set of labels li (t) where i is a node and t is a time instance. Each label li (t) = (F C, R, j, prev ptr) corresponds to a path from the source
300
node s to the node i arriving at i at time t and comprises four components. Specifically, the total fuel consumption F C and risk R represents the total fuel consumption and risk of a so-far constructed path from s to i, the integer j is 13
ACCEPTED MANUSCRIPT
the predecessor of i on the path from s to i and the pointer prev ptr points to the Pareto optimal label of the predecessor j whose extension gave this label. As has already been mentioned, leaving as earlier as possible from a node
305
CR IP T
does not necessarily reduce the fuel consumption or risk of the outgoing edge. Thus, there is no easy way to know in advance the ideal time of leaving a node for obtaining the minimum cost path toward the destination node. Consequently,
for each node and for each realized arrival (and hence departure) time instance 310
at this node, we should keep all the non-dominated labels referring to that particular arrival time. In contrast, in the corresponding static bi-criteria shortest
AN US
path problem there is no need of partitioning the labels at each node according
to their arrival time and only a single set of non-dominated labels needs to be kept at each node.
The Ship Routing Algorithm 1 iterates over all integer time instances, with
315
step size of 1 minute, in the interval [tstart , tstart + T ], where tstart is the departure time from the source node s and T is the maximum permissible voy-
M
age duration. At each node, the algorithm maintains two groups of label lists namely, the permanent lists and temporary lists. Each of these lists corresponds to a certain realized arrival time and either contains the final Pareto-optimal
ED
320
paths “arriving” this node at that particular time instance or contains the sofar Pareto-optimal paths reaching this node at that instance. When the outer
PT
loop scans a time instance, all the labels of the temporary list of a node for that instance are transferred one by one to the corresponding permanent list. However, just before the transfer, these labels are extended thereby creating
CE
325
new labels. Each of these new labels is kept only if it does not violate the total travel time constraint (line 19) and is not dominated by another label belonging
AC
to the same node with the same arrival time (lines 22-38). Also, labels already in the temporary list for that arrival instance are compared against this newly
330
generated label and discarded when found to be dominated by the new label.
14
ACCEPTED MANUSCRIPT
Algorithm 1 The time-dependent bi-objective shortest path algorithm Require: G = (V, A), and C, the cost matrix for all edges (i, j) ∈ A Ensure: All Pareto optimal paths from the node s to the node d s: the start node d: the destination node tstart : the departure time T: the maximum allowed total travel time
CR IP T
1:
trt(i, j, t): the travel time for traversing the edge (i, j) departing from the node i at time instance t
AN US
tjar : an arrival time instance at the node j
li (t): a label of the node i corresponding to the path from s to i, arriving at i at time instance t.
Ltempi (t): the list of temporary labels of node i at time instance t Lpermi (t): the list of permanent labels of node i at time instance t card(li (t), Lpermi (t)): the position of li (t) in the list of permanent labels
M
of node i
th lp component of a label li (t) i (t): the p
ED
/* Initialization of temporary and permanent labels of every node and for all t ∈ {tstart , tstart + 1, ..., T + tstart } Ltempi (t), Lpermi (t) ← ∅,
3:
Ltemps (tstart ) ← {[(0, 0); null; null]}
4:
∀i ∈ V, ∀t ∈ {tstart , tstart + 1, ..., T + tstart }
for tar = tstart to T + tstart do
CE
5:
PT
2:
6:
AC
7:
8: 9:
10: 11: 12: 13:
*/
while (∪i∈V Ltempi (tar ) 6= ∅) do Select a pivot label li∗ (tar ) from ∪i∈V Ltempi (tar )
/* Remove li∗ (tar ) from Ltempi (tar ) and add it to Lpermi∗ (tar )
*/
Ltempi∗ (tar ) ← Ltempi∗ (tar )\{li∗ (tar )} Lpermi∗ (tar ) ← Lpermi∗ (tar ) ∪ {li∗ (tar )} /* Store the position of label li∗ (tar ) in the list Lpermi∗ (tar ) pos ← card(li∗ (tar ), Lpermi∗ (tar ))
if i∗ 6= d then
15
*/
ACCEPTED MANUSCRIPT
14:
/* Label all the successors of i∗
*/
∗
for all do(i , j) ∈ E
15:
tjar ← tar + trt(i∗ , j, tar )
16:
/* Check the total duration constraint
18:
g←
tjar
19:
h ← duration of the shortest path from node j to node d
20:
f ←g+h
− tstart
if tjar − tstart + h ≤ T then lj (tjar )
22:
r(i∗ , j, tar )); i∗ ; pos)
if 6
24:
∃lj0 (tjar )
∈
Ltempj (tjar )
:
˜ j (tjar ) lj0 (tjar )≺l
Ltempj (tjar )
←
Ltempj (tjar )
∪
*/
then
/* Store the label lj (tjar ) of node j at time instance tjar as temporary
26:
*/
{lj (tjar )}
/* Delete all temporary labels of node j at time instance tjar dominated by
M
27:
li∗ 2 (tar ) +
/* Check that there is no label lj0 (tjar ) of node j at time instance tjar domi-
nating label lj (tjar )
25:
((li∗ 1 (tar ) + f c(i∗ , j, tar ),
←
AN US
21:
23:
*/
CR IP T
17:
lj (tjar )
Ltempj (tjar )
ED
28:
Ltempj (tjar )
←
and
end if
29:
*/ Ltempj (tjar )\{lj0 (tjar )
lj (tjar )
≺
lj0 (tjar )}
end if
PT
30:
end for
31:
end if
CE
32:
end while
33:
end for
35:
/* Output:Lpermd , all Pareto optimal paths from the node s to the node d
AC
34:
*/
36:
return Lpermd
The order by which the labels of the current time instance are extended is not actually important, since all labels of that instance should be extended first 16
∈
ACCEPTED MANUSCRIPT
before proceeding to the next iteration. Therefore, the label to be extended (pivot label) can be selected randomly (line 8). By the end of the execution, 335
all Pareto-optimal solution paths will have been stored as permanent labels of
CR IP T
the destination node d. Notice also that for the node d, no label extension takes place (line 14) since the algorithm does not need to go further down the destination node due to the fact that the cost functions are non-negative-valued and hence going further only adds to the total cost of the path. Finally, Ship 340
Routing Algorithm 1 assumes a single departure time instance from the source
node. Apparently, the algorithm can be easily extended to handle the possibility
AN US
of more than one departure time instances.
Also, for improving the performance of the algorithm, a heuristic function, denoted by h, is used for estimating a lower bound of the travel time from the 345
current position to the node d. Summing the current cost of the path, denoted by g, and the h value, a lower bound, f , for the travel time from s to d is obtained. This estimation is used to check if the extension of the current, partial, path
M
will violate the constraint of the total travel time. Using this heuristic, we can prune paths from the very early stages of the algorithm, since it is certain that these paths will not reach d in time. The h function is computed by the single-
ED
350
objective one-to-all Dijkstra, i.e. all optimal paths in terms of travel time are determined from d to any other node of the grid. The edge weights used in this
PT
computation are the least travel time of edges occurring within the time window [tstart , tstart + T ].
In quest of an efficient algorithm which will optimally solve the ship routing
CE
355
problem, we considered several enhancements of the algorithms in [11]. First, we converted the Continuous time Algorithm of [11] to a time-dependent multi-
AC
objective A∗ algorithm namely, to a time-dependent N AM OA∗ -like algorithm [42], by employing heuristic functions for the objective values of risk and fuel
360
consumption [29]. However, this contributed to a slightly faster algorithm only. Then, we applied the pruning technique with the h function described above to
the Continuous time Algorithm and the new algorithm was considerably faster than the original one [29] in the maritime setting. Since, the Discrete time 17
ACCEPTED MANUSCRIPT
Algorithm in [11] is generally faster than the Continuous time algorithm, we 365
finally decided to use the heuristic function h in the Discrete time algorithm. Notice also that techniques based on the heuristic functions for fuel consumption
CR IP T
and risk cannot be easily incorporated in the Discrete-Time algorithm because its logic differs a lot from that of the a time-dependent N AM OA∗ . The efficiency of the new version of the discrete time Ship Routing Algorithm (Algorithm 1) is 370
evaluated through a number of experiments. Also, both versions of the Discrete time algorithm (the initial and the new) were compared against the backward
algorithm of Hamacher et al’s [28] which was slightly adjusted to the problem at
AN US
hand. Hamacher et al’s algorithm tackles the same problem but with arbitrary arrival time at the destination node, while not taking into account possible 375
constraints on the arrival time at the destination.
Specifically, the algorithms were tested on the same set of randomly generated graphs. The size of the random graphs in the experiments was set to 500, 750, 1000 and 1500 nodes while the number of ingoing and outgoing edges
380
M
was 12, 14, 16, 18 and 20 in different experiments. All the nodes in the graph had the same number of incoming and outgoing edges. The single departure
ED
time instance from the source node tstart was set to 0. The time-dependent edge costs were randomly selected from the set {0, 1, ..., 10} for each time instance t ∈ {0, 1, ...T }. The maximum permissible total travel time T was set equal 385
PT
to 100 and 500. We computed the average execution time over 50 different experiments. The computational results are presented in Table 1.
CE
The experimental results demonstrate that the solution pruning based on the heuristic estimation of the remaining travel time reduces the execution time. In all the experiments, Ship Routing Algorithm 1 outperforms both Discrete time
AC
Algorithm of [11] and Hamacher et al’s algorithm. Specifically, the Discrete
390
time Algorithm of [11] is slower than Ship Routing Algorithm 1, but in some cases namely, in small graphs, appears to be faster. A possible reason for this is that the pruning technique causes more delay than speeding-up when applied on small graphs. Also, the performance of all algorithms are affected by the size of the random graph and its topology, but Ship Routing Algorithm 1 seems 18
ACCEPTED MANUSCRIPT
Table 1: Comperative computational results - departure time tstart = 0, maximum travel time T = 100
1500
12
1
1.251
14
1.24
1.26
16
2.6
2.5
18
5
4.74
20
6.2
5.22
12
1.54
1.6
14
2.24
2.18
16
3.82
18
6.94
20
7.98
12
2.96
14
4.12
16
5.74
18
7.68
20
9.2
12
8.78
8.34
98.56
14
10.46
9.82
123.22
13.22
12.4
141.4
18.60
17.74
179.34
20.38
18.66
200.86
16 18
395
23.42 36.4
50.86 61.36 83.2 42.6
63.44
3.36
78.9
6.12
93.76
7.08
110
2.9
53.24
3.74
69.98
4.98
90.5
6.88
104.76
8.6
135.7
PT
20
Hamacher’s Algorithm
AN US
1000
Algorithm 1
M
750
Discrete time Algorithm of [11]
ED
500
m/n
CR IP T
Average CPU Time (in seconds) n
to be less affected relative to the other two algorithms. Finally, Ship Routing
CE
Algorithm 1 is much faster than Hamacher et al’s algorithm which implies that a forward algorithm can solve this particular problem more efficiently than a
AC
backward one.
400
Admittedly, the use of the heuristic function does not improve the execution
time significantly when the algorithms run on random graphs. However, in the maritime setting, the situation is completely reversed and the use of the heuristic function greatly speeds up the computation (Table 2). Notice that in random graphs the edge costs are accordingly random but in graphs with maritime data 19
ACCEPTED MANUSCRIPT
Table 2: Comparison of Discrete time Algorithm of [11] and Algorithm 1 in maritime setting. The experiment setup is the same as that of Table 8 in Section 6. Discrete time Algorithm of [11]
Algorithm 1
(in hours)
CPU time (in seconds)
CPU time (in seconds)
Piraeus
Paros
3
4
Milos
Nauplio
4.5
20
Piraeus
Serifos
3
12
Chania
Paros
5
29
Piraeus
Syros
3.5
12
Piraeus
Syros
4
35
3 8 5
18 7
13
AN US
Destination
CR IP T
Max Travel Time
Start
the edge costs, especially the time-dependent travel time cost, are correlated 405
with the weather and sea conditions. Thus, in this context, the range of values of passage time cost is relatively short, depending on the weather changes. As a result, the heuristic estimate of the total travel time is commonly a good ap-
M
proximation. In contrast, in random graphs there is no similar guarantee about the heuristic estimate and thus there is a high probability that this estimate will 410
not be close to the real total travel time. Finally it should be mentioned that
ED
when the algorithms run on graphs with maritime data, heavier computation is usually required for processing each label than that in random graphs. Thus, the decrease in the number of labels considered due to the use of the heuristic
when working on random graphs.
CE
415
PT
function has a greater impact on the execution time in maritime setting than
4. Modelization Issues
AC
There are many critical modelling issues that directly affect the performance
of a ship weather routing method. The input grid and the kind of information stored in this structure affect the performance of the algorithm that will be used
420
eventually. It is also important to have an accurate model that predicts the ship response to the weather conditions as well as the impact of these conditions on the actual speed of the ship. A realistic, yet easily implemented model for 20
AN US
Figure 1: Adjacent nodes
CR IP T
ACCEPTED MANUSCRIPT
the fuel consumption of a ship is also needed. Finally, since our main goal is the design of a green ship routing method, the guidelines of the International 425
Maritime Organization (IMO) about voyage safety should also be taken into account.
M
4.1. The Grid
As has been mentioned above, the underlying grid structure is important for
430
ED
the performance of the ship routing algorithm. The grid used in our algorithm consists actually of two different grids namely, the Static and the Dynamic Grid, whose edges correspond to moving directions allowed between grid points. A
PT
lot of different data are spatially stored in the nodes and the edges of these two grids, such as the risk of traversing an edge, its length etc. Furthermore,
CE
the two grids together offer a fast and efficient way to retrieve and summarize 435
important information regarding maritime, geographic and weather data from the area of the Aegean Sea in Greece. In the following, more details about the
AC
two grids are given. 4.1.1. Static Grid
440
The Static Grid is the initial graph structure used for storing the static data.
These data include not only geographic and bathymetric information, but also information about the statutory protected areas of the Aegean Sea. The graph 21
ACCEPTED MANUSCRIPT
is composed of a two-dimensional array. Each array cell corresponds to a specific point on the map. However, not all points are valid; specifically, the following
• the point should be located on the sea,
445
CR IP T
conditions should hold:
• the distance of the point from the shore should be greater than a particular threshold, which is given as an input parameter,
• the point should be located in sufficiently deep waters; the minimum depth is also given as an input parameter,
AN US
• the point should not be located past the country’s maritime borders.
450
If a point satisfies all the conditions above, a new node is created and assigned to corresponding array cell. The so created nodes are connected with each other via bidirectional edges. A node can be connected to at most 16 adjacent nodes with different directions (Figure 1). The coordinates of the grid nodes are such that the distances between successive grid nodes are approximately the same
M
455
and still allowing a ship to move along the edges with constant bearing (rhumb
ED
line navigation). Again, some edges may not be valid and should be removed. Specifically, for an edge to be valid, the conditions mentioned above for points must hold for every point along the edge. It is now clear that for constructing the Static Grid the following parameters must be given as input:
PT
460
• the longitude and latitude of top-left corner point,
CE
• the edge length (for horizontal and vertical edges only),
AC
• the number of rows and columns in the array,
465
• the aforementioned parameters regarding minimum depth and distance from land.
4.1.2. Dynamic Grid The Dynamic Grid enriches the static information from the Static Grid with dynamic data. These data include weather forecast, such as wind direction and 22
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 2: Dynamic Grid
speed, significant wave height, temperature, as well as risk estimations resulting from these meteorological data [43]. For the weather information, both current
ED
470
and future forecasts are retrieved from the National Weather Service. Since, these predictions have a certain time granularity, for time instances lacking
PT
forecast, the weather conditions at each node and edge are determined by interpolating the nearest in time predictions. Also, for each edge, the risk a vessel 475
incurs along the edge given the predicted weather conditions is estimated.
CE
After finishing this processing, the dynamic and the static grid is combined
into a single grid. Figure 2 depicts that grid for some region of the Aegean
AC
Sea, namely the Cyclades Archipelagos. The total number of grid nodes for the whole Aegean Sea is 10, 000, the length of horizontal and vertical edges is 3
480
miles while the maximum distance (long diagonal lines) between two adjacent nodes is equal to 6.9 miles. At this point, it is worth mentioning that while the Static Grid is only produced once for the entire area of the Aegean sea, the Dynamic Grid should 23
ACCEPTED MANUSCRIPT
be rebuild regularly whenever newly updated data are available. For instance, 485
in our application scenario, dynamic grid should be rebuild every three hours
CR IP T
since this is the update rate of weather forecasts for the Aegean sea. 4.2. Control variables
The course of the ship is the only control variable used in our approach.
Specifically, the ship power or equivalently the nominal speed does not change 490
during the voyage, although there are routing methodologies [44], [19], [17], in
which the power and the nominal speed are control variables. The frequent
AN US
change of the ship speed is not recommended, especially in case of coastal navigation. Admittedly, ship power/speed may need to change in order to delay or speed up the passage through a region that currently has or will have ad495
verse weather conditions respectively. However, this can be handled by simply restarting the execution of the proposed algorithm giving as an input new operational parameters. Moreover, changing speed in order to avoid adverse weather
ED
4.3. Ship response
M
and sea conditions is more relevant in transoceanic seafaring.
In order to calculate the navigation time between two points, we have to
500
consider the actual ship speed which is usually lower than the nominal one
PT
due to the added resistances induced by irregular waves and wind during ship navigation. To this end, the ship model described in [4] is used. This model is general, independent of specific ship features. According to this model, the speed reduction depends only on the significant wave height H and the ship-
CE
505
wave relative direction Θ. The actual ship speed VR is calculated using the
AC
following equation: VR = V (H, Θ) = V0 − f (Θ) · H 2
The values of the coefficient f are reported in table 3.
510
If more specific ship characteristics are available, a more realistic model for the vessel response to weather and sea conditions could be used. Specifically,
24
ACCEPTED MANUSCRIPT
o
45 < Θ < 135
o
135o ≤ Θ ≤ 180o
f [kn/ft2 ]
following seas
0.0083
beam seas
0.0165
head seas
0.0248
Table 3: Values of the coefficient f
CR IP T
0o ≤ Θ ≤ 45o
Configuration name
we could use the approximate method established by Kwon [45] for predicting
AN US
the loss of speed caused by the added resistance in irregular waves and wind. 4.4. Fuel consumption
The fuel consumption F measured in kg per hour is calculated using the
515
following equation [9], [46]:
F = sf c · P
M
where sf c and P is the specific fuel consumption and the engine power in BHP (kW) of the ship respectively.
Thus, the total fuel consumption along a route is the product of its travelling
ED
520
time T rT and the fuel consumption per hour F :
PT
F Ctotal = T rT · F
4.5. Voyage safety/risk
CE
Voyage safety is ensured by taking into account the total risk of a route during route calculation as well as the IMO recommendations [1] for avoiding dangerous situations. There may be different ways of defining the risk, but the
AC
exact meaning is hard to be captured. According to the FSA Guidelines (MSC 83/INF.2) [47] risk “is the combination of the frequency and the severity of the consequence”. In line with the FSA guidelines, risk is defined as the probability of an event multiplied by its consequences (14): R=P ×S 25
(14)
ACCEPTED MANUSCRIPT
The proposed model combines the results of a marine accident probability (P ) 525
with the severity of its consequences (S). Concerning hazard identification, marine accidents are considered as poten-
CR IP T
tial hazards. These accidents can be divided into nine categories: foundered, missing vessel, fire/explosion, collision, contact, grounding/wrecked/stranded, war loss/hostilities, hull/machinery damage and miscellaneous. Fault and/or 530
Event tree Analysis, Fuzzy logic and Bayesian networks are the most frequently
used methodologies [48]. In our approach, we employ the model proposed in [43] which takes into account dynamic information about a significant navigation
using a Bayesian network method.
AN US
area as well as the prevailing weather conditions in that area and is developed
Besides the diversity of methods in risk assessment of maritime transporta-
535
tion, there are also many different factors that should be considered when evaluating the probability of an accident. Complex methodologies that take into account weather and environmental conditions, information about the vessel
540
M
and the human factor have been introduced to calculate the accidental risk of shipping in a sea area and identify its causes. In particular, in our approach risk
ED
assessment considers information such as the vessel type, size, age and flag, navigation area characteristics, traffic data, historic data about accidents, weather and sea conditions [48].
545
PT
With respect to the risk of a particular edge, this depends on the weather conditions at the location of the edge and is estimated using the model in [43],
CE
as well. Finally, the risk of a particular route is the total risk encountered along each edge of the route. Regarding IMO recommendations, surf-riding and broaching-to situations
AC
should be avoided when navigating in severe weather conditions. Surf-riding
550
and broaching-to occur when the conditions below are satisfied: √ 1.8 L 135o < Θ < 225o and VR > cos(180 − Θ) where Θ, VR and L are the relative ship-wave angle, the ship speed and the ship length, respectively. 26
ACCEPTED MANUSCRIPT
Since ship speed is constant, the only way to avoid surf-riding and broaching555
to situations is by rejecting the edges where the above conditions are arising. Parametric rolling motions is another precarious situation that should be
CR IP T
avoided. It will occur when the encounter wave period TE is nearly equal to the natural rolling period of ship TR or the encounter period TE is close to one half of the ship roll period TR . The period of encounter TE is calculated by the 560
following formula: TE =
3TW
2 3TW (sec) + VR cos(Θ)
ship speed (in knots). 4.6. Manoeuvring characteristics
AN US
where TW is the wave period, Θ is the relative ship-wave angle and VR is the
In the proposed algorithm, the turning ability of a ship was not considered
565
during path planning. Thus, when a solution path does not comply with the ship
M
turning ability, that path is smoothed by using Cubic Spline interpolation [49]. In particular, smoothing takes place locally between the two consecutive route
570
ED
edges that form the sharp turn and thus the obtained smoothed route will be within the three grid cells that the original route crosses. Therefore, it is reasonable to assume that the smoothed route will not violate any of the conditions
PT
listed in Section 4.1.1. Also, the total risk and the total fuel consumption of the smoothed route are expected to be roughly the same as those in the original route.
CE
If the maximum design rudder angle has to be considered during the cal-
575
culation of Pareto optimal routing paths, the approaches in [14] and [50] could
AC
be easily incorporated in our technique. In [14], a modified network graph is constructed by creating copies of the original grid nodes and each of these copies models a allowable turn at the corresponding original node. On the other hand,
580
the authors in [50] propose a modified Dijkstra algorithm which uses arc labelling for taking into account turn restrictions. Our bi-objective shortest path algorithm can be easily adjusted so that it can work with either of the two 27
ACCEPTED MANUSCRIPT
approaches. However, the execution time would be considerably longer since the modified graph would be much larger than the original one [14] or in the 585
case of [50] because working on arcs would be computationally much heavier
CR IP T
than working on nodes due to their much higher number. Moreover, since sharp turns were seldom observed in our case study, navigation in the Aegean sea, a simple solution with local curve smoothing was preferred eventually.
5. Speeding Up and Improvement Techniques
So far, the ship routing problem has been handled as a bi-criteria, time-
590
AN US
dependent and resource-constrained path problem whose parameters have been
properly set considering the particular application scenario. By closer examination of the common practices and assumptions in this application area, the basic algorithm can be improved so that it can run faster and return more prac595
tical solutions, that is, solutions with higher probability of being adopted by the
M
end-users. In this section, we present speeding up techniques and techniques which improve the practicality of the derived solutions. More specifically, we propose a technique for bounding the solution space around a voyage plan given
600
ED
in advance. We also present a grid partitioning method which can lower the execution time of the algorithm with the caveat that it can be applied only in
PT
the case that the cost functions in the input grid are time-invariant. Moreover, a technique that ensures that there will be no loops in the solution routes is presented. Finally, a selection tool is developed for presenting the most repre-
CE
sentative solutions to the end user when the number of Pareto-optimal solutions
605
is too large for clear illustration.
AC
5.1. Grid Pruning It is a fact that out of experience, shipping companies usually know the
best route between two ports and the voyage plan given to the captain of a ship reflects this knowledge. Since, most companies are rather reluctant to
610
accept large changes in the customary route they follow between two ports, the
28
AN US
Figure 3: Grid pruning
CR IP T
ACCEPTED MANUSCRIPT
proposed routing algorithm should derive solutions that do not much deviate from the original voyage plan of the ship. For this purpose, the original grid is pruned leaving only a “corridor” around the voyage plan (Figure 3). This
615
M
drastically reduces the solution space for finding the optimal solution. Also, the pruning can be easily performed online, during the execution of the routing algorithm. In comparison to the off-line execution, the on-line pruning examines
ED
fewer grid points, saving computation time, as a result. For evaluating the benefits of the grid pruning, we ran the routing algorithm with and without grid pruning. The thickness of the corridor around the voyage plan, distV P , was set to two different distinct values, 42 miles and 21 miles
PT
620
respectively (third line of Table 4). These distance values from the voyage plan
CE
were selected empirically. The tests were carried out on a server with a 4-core 3.6 GHz processor and 8 GB RAM and are listed in Tables 4 and 5. The experiment setup is the same as that of Table 8 in Section 6. In the case of 42 miles distance, the number of the derived solutions is the same as that returned
AC
625
in the original grid and thus we can still find the whole Pareto set. However, if distV P is set to a value lower than 42 miles, some Pareto optimal solutions are not returned (Table 4). Concerning execution time, clearly grid pruning along the voyage plan can reduce the execution time dramatically (Table 5).
29
ACCEPTED MANUSCRIPT
Table 4: Algorithm performance: original and pruned grid
Destination
Number of Routes
distV P
distV P
Max Travel Time (in hours) ∞
42 miles
21 miles
∞
42 miles
21 miles
33
33
29
47
47
47
26
26
26
CR IP T
Departure
CPU time (in seconds)
Paros
3
3
3
3
Milos
Nauplio
4.5
11
8
7
Piraeus
Serifos
3
5
5
5
Chania
Paros
5
36
18
17
50
50
50
Piraeus
Syros
3.5
7
7
7
38
38
36
Piraeus
Syros
4
14
13
13
38
38
36
Spetses
Amorgos
4
53
49
AN US
Piraeus
15
13
13
53
Table 5: Algorithm performance: original and pruned grid (distV P = 42 miles)
Execution time
Destination
M
Start
Original Grid
Pruned Grid
Elafonisi
01:49:42
00:01:09
Kymi
Kalamata
00:01:00
00:00:27
Chios
00:09:43
00:00:59
Cythera
00:04:43
00:00:50
Hydra
Samos
00:00:27
00:00:16
Rhodes
Western Crete
01:17:27
00:05:14
SW Crete
CE
PT
Piperi
ED
Kos
5.2. Dynamically Partitioned Grid
AC
630
In this section, we propose another speed-up technique which can be used
together with the basic ship routing algorithm. However, the main limitation is that it can be applied only for static grids, i.e., with constant costs at the edges and nodes. For short trips, where the weather conditions in the travel area are
635
not likely to change due to short travel time, static grids can be used for solving
30
ACCEPTED MANUSCRIPT
the routing problem. The main objective is the removal of nodes and edges from the input grid in such a way that most of the solutions from the original Pareto set can also be
640
CR IP T
obtained when working on the reduced grid. In addition, the routing algorithm should run faster on the new grid.
This size reduction is feasible since in the Aegean sea, there are broad areas
with few or even no islands where the weather conditions are practically the same. Thus, the objective is to locate convex polygons with same weather
within. Apparently also, a straight line between any two nodes inside such a polygon lies also in that polygon and thus all the points of the line will have the
AN US
645
same weather. Therefore, we can reduce the size of the grid, by maintaining only the boundary nodes of each polygon, termed also as gateways, and connecting them pairwise with straight line edges. The weather and risk condition along each of the new edges can be approximated by averaging the values of these 650
metrics at the end vertices.
M
A previous version of this speeding up technique can be found in [29]. In that approach, we extracted areas with uniform weather and risk conditions
ED
and then applied a convex partitioning algorithm for finding a more accurate partition. However, the number of partitions obtained were relatively high and 655
hence the size of the grid was not reduced a lot. In this work, much smaller grid
PT
is obtained by extracting only the largest convex polygons from each uniform area, maintaining the remaining nodes and edges of the original grid.
CE
The algorithm of creating the Dynamically Partitioned Grid consists of 3 phases.
5.2.1. Phase 1 - Partitioning the Grid into Sectors with Uniform Weather and
AC
660
Risk Conditions The main steps of this phase are given in Algorithm 2. For each node,
depending on the similarity of that node with the neighboring sectors, either it is assigned to an existing adjacent sector or a new sector is created, containing that
665
node. A node and a sector are considered similar if the risk, wind direction and 31
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 4: Overlapping sectors
wave height at the node and within the sector are also similar. For a sector, the
M
values of these parameters are estimated as the mean of the corresponding values at nodes of this sector and all these averages are recalculated upon insertion or
670
ED
removal of a node from the sector. Also, sectors connected diagonally with a node, are not considered unless no sectors exist at the north and west of the node. Otherwise, overlapping areas would be possible (Figure 4).
PT
After the assignment of nodes to sectors has been completed, the similarity of each node with the home sector, as well as with its neighbouring sectors,
CE
is calculated and the node is moved to the highest-scored sector provided that 675
the previous sector still remains connected after this movement. This node
AC
re-affiliation continues until no further change is possible. 5.2.2. Phase 2- Extraction of the largest convex polygon(s) from each sector Fischer’s algorithm [51] is used for extracting the maximum convex polygon
from each sector. This is a O(n4 logn) dynamic programming algorithm taking
680
as input a set of positive and negative points and returns the maximum convex polygon which encloses only positive points. For applying this algorithm, first,
32
ACCEPTED MANUSCRIPT
Algorithm 2 Partitioning the Grid in Sectors with Uniform Conditions Require: The 2-dimensional grid array containing all nodes
CR IP T
Ensure: All nodes with similar risk and weather conditions belong to the same sector 1:
Ni , Wi , NWi , SWi : the northern, western, north western and south western neighbour of node i sci : the sector node i belongs to
sim(arg1, arg2): compares either two sectors or a sector with a node and
AN US
returns a measure of the similarity of the two entities with regard to weather
and risk conditions; the more similar they are the higher the value returned same(arg1, arg2): compares either two sectors or a sector and a node and returns true if the similarity between the two entities is above a certain threshold.
M
merge(sc1, sc2): merges the two sectors given as arguments /* Initialization of sectors
3:
for each node i do/* the grid is parsed from left to right and top to bottom */
5: 6:
add i to scNi
end if
if ∃ Ni , Wi and same(scWi , i) and sim(scWi , i) > sim(scNi , i) then
CE
7:
if ∃ Ni and same(scNi , i) then
PT
4:
add i to scWi
8:
9:
AC
10:
11: 12:
*/
ED
2:
end if if ∃ Ni , Wi and same(scNi , i) and same(scWi , i) and same(scNi , scWi )
then merge(scNi , scWi ) end if
33
13: 14:
if 6 ∃ Ni , Wi then if ∃ N Wi and same(scN Wi , i) then add i to scN Wi
15:
CR IP T
ACCEPTED MANUSCRIPT
end if
17:
if ∃ SWi and same(scSWi , i) and sim(scSWi , i) > sim(scN Wi , i) then add i to scSWi
18: 19:
end if
AN US
16:
20:
end if
21:
if i not in any sector yet then
22:
Create a new sector for i end if
M
23:
end for
25:
/* Associate each node with a more similar local sector if exists
26:
Change = F alse
27:
repeat
29: 30:
if ∃ adjacent sector more similar than the current one then move i to the more similar adjacent sector Change = T rue
CE
31:
for each node i do
PT
28:
ED
24:
32: 33:
AC
34:
end if
end for
until Change = F alse
34
*/
ACCEPTED MANUSCRIPT
the bounding box of each sector sci is determined. Then, any node belonging to sci is considered as a positive point while any other node in the bounding box is treated as a negative point. Also, each sector should be connected and contain 300 nodes at most. Although, these requirements are not necessary for
CR IP T
685
the correctness of the algorithm, they both help in faster execution. Thus, in
case of not connected sector, the sector is split into its connected components. If a sector contains more than 300 nodes, the sector is split along the horizontal
or vertical line that partitions the nodes as evenly as possible and generates the 690
fewest possible gateways.
AN US
However, even for sectors with fewer than 300 nodes, the computation of
maximum convex polygon may not be fast. Thus, for each sector, a sparse sampling of representative positive points of the sector is given as an input to Fischer’s algorithm. Specifically, if the number of nodes in the sector is 695
100 to 200, nodes with odd row and column number are omitted, while for a number between 200 and 300, nodes with either odd row and column number,
M
or even row and column number are ignored. In both cases, nodes lying on the boundaries of the sector are never removed.
700
ED
After this preprocessing, Fischer’s algorithm is executed. If the algorithm is unsuccessful or the resulting polygon has too many gateways compared to the total number of nodes contained in the polygon, the polygon is rejected.
PT
Otherwise, the polygon is stored along with the nodes it contains. A node in a sector is a gateway if at least one of its 8 incident shorter edges is also incident
CE
to a node outside the sector. The 8 long diagonal edges are ignored at this 705
phase, because the gateways should be near the sector perimeter. In addition, we ignore nodes having only one short diagonal edge which is incident with an
AC
external node and also passes between two other gateways. The reason is that these nodes are located inside the sector (Figure 5) and are not part of the polygon perimeter.
710
Now, the remaining nodes of the bounding box not belonging to the polygon
just formed create a new sector and the previous steps are repeated again. This iteration continues until it is not possible to find a maximum convex polygon 35
CR IP T
ACCEPTED MANUSCRIPT
with low number of gateways.
AN US
Figure 5: A case where a node is not considered as gateway.
Finally, for further reducing the execution time, parallel processing is em715
ployed since the computation above is amenable to parallel processing as each sector can be independently processed by a separate thread in a multi-core sys-
M
tem.
5.2.3. Phase 3 - Building the new Grid
720
ED
In this phase, the final grid is obtained by adding or removing some edges from the original grid. Specifically, if a node is not a gateway, all its incident edges are removed. Nevertheless, these isolated nodes are kept, since it might
PT
be chosen as the starting or destination point of a route. In this case, new edges from or to this node are created dynamically, connecting it with the gateways
CE
of its sector. Essentially, the new grid contains only the gateways of the sectors 725
and the nodes that do not belong to any sector. Each pair of gateways belonging to the same sector is connected with a new bidirectional edge while nodes in
AC
different sectors are connected with edges of the original grid, including also the long ones. The costs (risk, fuel, time) of each new edge are calculated as the average of the corresponding costs of the nodes of the original grid met along
730
the edge, gateways inclusive. The derived partitioned grid (Figure 6), contains about 56% less nodes and 24% less edges compared to the original grid. The routing algorithm running on 36
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 6: The partitioned grid.
this grid finds a close approximation of the original Pareto set in considerably
735
M
less time. However, there is a better alternative. After running the routing algorithm on the partitioned grid, we extract the minimum risk and the mini-
ED
mum fuel path as well as any other Pareto-optimal path which is significantly different from the two paths above. Then, these paths are used as a sort of voyage plan and any nodes not near these paths are removed, following again
740
PT
the technique described in Section 5.1. Next, the ship routing algorithm is run on the pruned grid with much shorter execution time. Again, most of the initial
CE
Pareto optimal solutions are extracted, if not the entire set. 5.3. Avoiding Loops It is a well known fact that if waiting at nodes is not an option, the optimal
AC
solutions for the shortest path problem may contain loops [24]. In our appli-
745
cation scenario, this is rather a rare case and it mainly happens when there is an abrupt change in the weather conditions in a local region. Since, the spatial as well as the temporal transition between different weather conditions is rather smooth most of the time, routing loops are not frequently observed in
37
ACCEPTED MANUSCRIPT
the derived routes. Indeed, only in two out of 100 execution runs of the ship 750
routing algorithm on random instances2 , few Pareto-optimal solutions are found to contain loops.
CR IP T
However, a route containing a loop is viewed undesired in ship navigation and it is usually dismissed. A possible approach for deriving loop-less ship routes is to directly modify the Algorithm 1 so that it will maintain only paths 755
that do not form loops at each execution step. However, as has been pointed
out in [23, 24] where special cases of the problem at hand were treated, finding
routes without loops requires an enumeration of all paths between source and
AN US
destination, which increases the state space that should be explored dramatically and hence the execution time. The main reason is that the algorithm should 760
“remember” the intermediate nodes of each path so that it will be easy to check then if an extension of the path is actually feasible, that is, it does not contain nodes previously met in the path. In particular, in our scenario, it is not certain that keeping only Pareto-optimal solutions at each node for each arrival time
765
M
will give the Pareto-optimal routes between the source and the destination, eventually. Thus, non-Pareto-optimal solutions should be kept as well, and this
ED
leads to the complete enumeration of the routing paths mentioned above. Even if we disregard these technical difficulties, the modified Algorithm 1 may derive routes that will not be sought-after again. Specifically, the algorithm
770
PT
may return routes that do not contain “exact” loops but they may pass again through regions very close to locations already traversed. Apparently, this sort
CE
of routes can be also be viewed as routes with routing loops. However, detecting these cases would add additional complexity to the Algorithm 1. For all these
AC
reasons, the option of modifying the basic algorithm was not followed.
775
A more practical approach is to modify the input grid around the voyage
plan so that it will be made acyclic and hence routing loops will not be possible in the returned solutions. This can be done on-line, during the algorithm execution; specifically, all edges of the grid whose angle with the current di2 randomly
chosen source and destination ports and weather conditions
38
ACCEPTED MANUSCRIPT
rection of voyage plan is greater than 90o are discarded. It can be easily seen that the resulting grid after these edge removals becomes acyclic. Applying this 780
technique, optimal solutions may be discarded since grid edges are removed,
CR IP T
however, no noticeable difference in the solution quality was observed in our experiments. More precisely, the proposed technique was tested in 100 different experiments, with random departure and destination nodes. In 94% of the test
cases, the algorithm running on the modified grid retrieved the entire Pareto 785
set which is returned when the algorithm is running without the loop avoidance
technique. In all other test cases, where the whole Pareto set was not returned,
AN US
the quality of the solution routes was only slightly affected, since 92% of the original Pareto set was found again, on average. For instance, in the test cases of Table 4 the algorithm with the loop-avoidance technique found the whole 790
Pareto set in all itineraries, except in the case where the departure and destination ports are the cities of Chania and Paros respectively (line 4 of Table 4); specifically, the algorithm returned 48 solution routes and retrieved 47 out of
M
50 optimal routes that the basic algorithm found. In Figure 7, the solutions of the basic algorithm and those returned by our approach for the voyage Chania to Paros are plotted with blue and red color respectively. Clearly, the difference
ED
795
of the two Pareto-fronts is hardly discernible. Moreover, if the first approach of the direct modification of Algorithm 1 had been followed, the Pareto-front that
PT
would be returned in that case would closely match the Pareto front derived by our approach. This is because the Pareto front derived when the basic algorithm is running without a loop avoidance technique always “dominates” the Pareto-
CE
800
fronts of any loop-avoidance implementation. Since also, the Pareto-front of the directly modified Algorithm 1 would normally be in between the two fronts of
AC
Figure 7, the close matching between the fronts depicted in that figure certainly entails a corresponding matching in the results of our approach and those that
805
would be returned with a direct modification of Algorithm 1. Finally, it is worth mentioning that if the geographical morphology of a
region requires that the route should form an exact or “approximate” loop, the proposed approach can derive this sort of loop provided that the voyage 39
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 7: Chania to Paros: Pareto fronts of the basic algorithm and loop avoidance technique.
plan contain the required cyclic segments. This is possible, since the grid lines 810
following the current direction of the voyage plan are always included in the
5.4. User Selection
M
working grid.
ED
After the execution of the ship routing algorithm, the set of Pareto optimal routes is illustrated on the map (see Figure 12, for instance). However, if the 815
number of these solutions is very large, it is more convenient to illustrate a
PT
smaller number of representative solutions. Specifically, three routes are selected from the Pareto optimal solutions in this case as follows:
CE
1. The first route minimizes the risk and is selected among these solutions where the fuel consumption is no higher than c times the minimum fuel
AC
820
consumption attained over all the Pareto optimal routes where c is an input parameter with a default value between 1.2 and 1.25.
2. The second route minimizes the fuel consumption conditional on the risk of this route being no higher than c times the minimum risk achieved among all the routes of the Pareto set.
40
ACCEPTED MANUSCRIPT
3. The third route minimizes a weighted average of both the risk and fuel
825
consumption among the remaining routes of the Pareto set. Apparently, by changing the value of parameter c, the whole Pareto front
CR IP T
can be gradually returned. Another possibility is to use a measure of route similarity for comparing routes and eliminating routes similar to previously 830
selected routes. Thus, only one route will be kept out of each group of routes following broadly the same trajectory in the final output.
The selected routes are returned to the user, as a sequence of grid points
plotted on a map along with the coordinates of each point. Moreover, an expla-
835
AN US
nation about each route is included in the output. An explanation of this kind
could read as “This route minimizes the risk and has only 17% higher fuel consumption than the most fuel economical route (this route would be 40% more risky)”. Finally, this output is illustrated on a graphical user interface and then the user can choose the most preferable of these routes. For instance, Figure 8
840
M
illustrates the output of the algorithm with Spetses and Amorgos as the departure and destination port respectively. The input voyage plan is illustrated in
ED
Figure 3 and the corridor around the voyage plan was set to 42 miles.
6. Experimental Results
PT
In this section, the experimental evaluation of the baseline algorithm as well as its variants is presented. Firstly, the algorithm based on the dynamically 845
partitioned grid is tested and compared with the case where the algorithm uses
CE
only the original grid as input. Then, the performance of the algorithm enhanced with the loop avoidance and the grid pruning technique is tested in a number
AC
of real-life test instances.
850
In order to test the Dynamically Partitioned Grid, we generated 100 random
test cases. For each test case, a random set of starting and destination ports were chosen and the maximum travel time for each voyage was set proportional to the euclidean distance from the source to the destination port. Since, the technique of Dynamically Partitioned Grid cannot be applied in time-dependent 41
ED
M
AN US
(a) The route minimizing the risk
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
(b) The route minimizing the fuel consumption
(c) The route minimizing the weighted average of risk and fuel cost Figure 8: The three solutions returned to the user (distV P = 42 miles).
42
ACCEPTED MANUSCRIPT
grid networks, the weather and sea conditions do not vary in time but they can 855
be spatially different in these test scenarios. The other test parameters were fixed as follows: speed = 30 kn, s.f.c. = 200 g/kWh, break power = 4000 kWh.
CR IP T
Also, the tests were carried out on a server equipped with a 4-core 3.6 GHz processor and 8 GB RAM.
Some indicative experimental results of these tests are listed in Table 6. It 860
should be mentioned that the baseline algorithm runs on the original grid without grid pruning and loop avoidance technique. Also, in the estimation of the execution time, we did not take into account the time for creating the Dynam-
AN US
ically Partitioned Grid. This is because in our test setting there is no temporal
change in weather conditions and thus the dynamically partitioned grid is cre865
ated only once during initialization. The experimental results suggest that the dynamically partitioned grid technique accelerates the ship routing algorithm, without conspicuously decreasing the quality of the derived solutions. Specifically, the technique of using solution paths from the dynamically partitioned
870
M
grid for creating the voyage plan on the original grid is 2.218 times faster on average than the baseline algorithm; at the same time it is able to generate the
ED
88.32% of the Pareto set approximately. Noteworthily, when the execution time of the baseline algorithm is 5 minutes or more, the algorithm on the dynamically partitioned grid is 4.166 times faster with the greatest improvement being 70-
875
PT
fold. Regarding the solution quality, when using the Dynamically Partitioned Grid, there is an average increase of 0.037% in the risk and 0.051% in the fuel
CE
consumption of the safest and the most fuel economical route respectively in comparison to the original grid. Next, we evaluated the ship routing algorithm when using the grid pruning
AC
(Sect. 5.1) and loop avoidance technique (Sect. 5.3). Besides the standard grid,
880
in these experiments we generated two other grids for the Aegean Sea with coarser and finer resolution respectively (Table 7) for studying the impact of grid resolution on the execution time and the output of the algorithm. For each grid, we generated 100 random test cases, that is, for each case, a random departure and destination port was chosen and the maximum travel time of 43
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
PT
ED
Figure 9: Source-Destination ports of Table 6.
Table 6: Comparison of the algorithm performance on the Original and on the Dynamically
CE
Partitioned Grid
Relative Increase in the Dynam-
Execution Time
ically Partitioned Grid over the Original Grid
Max Travel Time
Percent of Common Routes
Best Risk
Best Fuel Consumption
Original Grid
Dynamically Partitioned Grid
Kos - Elafonisi
10 hrs
82.76%
0%
0%
01:49:42
00:01:33
Kymi - Kalamata
10 hrs
100%
0%
0%
00:01:00
00:00:29
SW Crete - Chios
11 hrs
68.33%
1.64%
0.003%
00:09:43
00:01:14
Piperi - Cythera
10 hrs
100%
0%
0%
00:04:43
00:00:50
Hydra - Samos
8 hrs
83.33%
0%
1.04%
00:00:27
00:00:22
Rhodes - W Crete
13 hrs
99.11%
0%
0%
01:17:27
00:13:10
AC
Route
44
ACCEPTED MANUSCRIPT
Table 7: Grid Information
Coarse
Fine
Grid
Grid
Grid
(SG)
(CG)
(FG)
10,000
6,400
40,000
Distance of Nodes (in miles)
2
2.5
1
Number of Neighbor Nodes
16
16
16
ED
M
AN US
Grid Nodes
CR IP T
Standard
Figure 10: Source-Destination ports of Table 8.
voyage was set again proportional to euclidean distance of the departure and
PT
885
destination port. In all experiments, we set the maximum distance distV P from
CE
the voyage plan equal to 42 miles. Also, in contrast to the previous tests, now the input grids are time-dependent. However, the same weather conditions were assumed for all grids and thus the same time-dependent cost functions were used in all tests concerning each particular grid. Table 8 shows the results of these
AC 890
experiments. Clearly, the proposed algorithm enhanced with the grid pruning and loop
avoidance technique can efficiently solve the point-to-point ship routing problem. The increase in the maximum travel time affects the execution time (lines 5-6,
45
ACCEPTED MANUSCRIPT
895
Table 8) but it does not change the number of the derived solutions. As a side note, since the route planning algorithm is fast when working with a voyage plan, the algorithm can run from scratch and return new Pareto-optimal routes
CR IP T
in case of an unforeseeable situation during the journey. Concerning the impact of grid resolution on the execution time of the algo900
rithm, we can see from Table 8) that the algorithm runs faster on the Coarse
Grid than on Standard Grid except in two cases where the execution time is the same. On other hand, the execution time of the algorithm on the Fine Grid is longer in most tests; it is the same only in one case. The largest increase
905
AN US
in run-time is 16% on the Fine Grid, while the largest decrease is 12% on the
Coarse Grid. This should be expected since the number of edges and vertices of the grid directly affects the number of iterations of the algorithm. However, it should be mentioned again that the algorithm works on a pruned grid around the voyage plan. Otherwise, the impact of the grid resolution on the execution time of the algorithm would be more noticeable since in that case the algorithm search space would be the whole grid and the three grids differ a lot in the
M
910
number of nodes and edges that they contain.
ED
With regard to the solutions returned, the number of the Pareto-optimal solutions is only slightly affected by the resolution of the grid (Table 8). Figure 11 also illustrates the Pareto-fronts obtained for a particular journey, namely the journey from Piraeus Chania. We can see that solutions derived on the Fine
PT
915
Grid are slightly superior and those of the Coarse Grid a little inferior to the
CE
solutions of the Standard Grid. This invariance of the solution quality can be attributed to the use of the pruned grid and the fact that the weather conditions do not abruptly change on a local scale, most of times. Thus, an increase of decrease of detail by using a finer or coarser grid respectively may not be
AC
920
essential due to the same weather over large parts of the sea. Finally, we ran the algorithm with the two enhancements above for finding
the Pareto-optimal routes between two specific ports, namely Piraeus and Chania of the Crete island. In Figure 12, we can see the Pareto front derived by
925
the algorithm, with the maximum distance from the voyage plan set equal to 21 46
ACCEPTED MANUSCRIPT
Departure
Destination
CR IP T
Table 8: Computational results: Execution time
CPU time (in seconds)
Number of Routes
SG
CG
FG
SG
CG
FG
4
33
30
34
9
47
44
49
5
26
26
26
Max Travel Time (in hours)
Paros
3
3
3
Milos
Nauplio
4.5
8
7
Piraeus
Serifos
3
5
4
Chania
Paros
5
Piraeus
Syros
3.5
Piraeus
Syros
4
AN US
Piraeus
16
21
50
45
51
7
6
8
38
36
38
13
13
15
38
36
38
AC
CE
PT
ED
M
18
(a) Standard Grid
(b) Coarse Grid
(c) Fine Grid
Figure 11: Pareto Front from Chania to Paros
47
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 12: Grid pruning and loop avoidance technique: Pareto front for the voyage between
M
Piraeus and Chania.
miles and the maximum duration of the voyage up to 3 hours. The ship routing
ED
algorithm found the whole Pareto set which comprises 48 routes.
PT
7. Conclusion
In this paper, we developed an exact algorithm for solving the time dependent, multi-objective and time constrained ship routing problem with fixed
CE
930
departure time. The final routes are optimal with respect to two different conflicting objectives namely, the fuel consumption and the risk along the route.
AC
In order to enhance the safety of the voyage, we have also considered the IMO regulations. Due to heavy computational demands, we focused on reducing the
935
execution time by reforming the grid structure. We have also enhanced the quality of the solution routes by using techniques which ensure that no loops will exist in the derived solutions and also increase the similarity of the solution routes to a given voyage plan. The experimental evaluation of the proposed ap48
ACCEPTED MANUSCRIPT
proach confirms that our techniques efficiently solve the ship-routing problem 940
and they are suitable for real-time applications.
CR IP T
References [1] M. IMO, 1/circ. 1228, Revised guidance to the master for avoiding dangerous situations in adverse weather and sea conditions, adopted 11th January.
[2] C. P. Padhy, D. Sen, P. K. Bhaskaran, Application of wave model for weather routing of ships in the north indian ocean, Natural Hazards 44 (3)
945
AN US
(2008) 373–385.
[3] R. Szlapczynski, A new method of ship routing on raster grids, with turn penalties and collision avoidance, Journal of Navigation 59 (01) (2006) 27– 42. 950
[4] G. Mannarini, G. Coppini, P. Oddo, N. Pinardi, A prototype of ship
M
routing decision support system for an operational oceanographic service, TransNav, the International Journal on Marine Navigation and Safety of
ED
Sea Transportation 7 (1) (2013) 53–59. [5] W. Shao, P. Zhou, S. K. Thong, Development of a novel forward dynamic programming method for weather routing, Journal of marine science and
955
PT
technology 17 (2) (2012) 239–251. [6] A. Azaron, F. Kianfar, Dynamic shortest path in stochastic dynamic net-
CE
works: Ship routing problem, European Journal of Operational Research 144 (1) (2003) 138–156.
[7] H. K. Lo, M. R. McCord, Adaptive ship routing through stochastic ocean
AC
960
currents: General formulations and empirical results, Transportation Research Part A: Policy and Practice 32 (7) (1998) 547–561.
[8] A. A. Montes, Network shortest path application for optimum track ship routing, Tech. rep., DTIC Document (2005).
49
ACCEPTED MANUSCRIPT
965
[9] K. Takashima, B. Mezaoui, R. Shoji, On the fuel saving operation for coastal merchant ships using weather routing, in: Proceedings of Int. Symp. TransNav, Vol. 9, 2009, pp. 431–436.
CR IP T
[10] I. S. Dolinskaya, Optimal path finding in direction, location, and time dependent environments, Naval Research Logistics (NRL) 59 (5) (2012) 325–339.
970
[11] A. Veneti, C. Konstantopoulos, G. Pantziou, Continuous and discrete time
label setting algorithms for the time dependent bi-criteria shortest path
AN US
problem, in: Operations Research and Computing: Algorithms and Software for Analytics, 2015, pp. 62–73. 975
´ Saux, C. Claramunt, A bidirectional path-finding algorithm [12] D. Tsatcha, E. and data structure for maritime routing, International Journal of Geographical Information Science (ahead-of-print) (2014) 1–23.
M
[13] L. Babel, T. Zimmermann, Planning safe navigation routes through mined waters, European Journal of Operational Research. [14] I. Ari, V. Aksakalli, V. Aydogdu, S. Kum, Optimal ship navigation with
ED
980
safety distance and realistic turn constraints, European Journal of Operational Research 229 (3) (2013) 707–717.
PT
[15] O. Kosmas, D. Vlachos, Simulated annealing for optimal ship routing, Computers & Operations Research 39 (3) (2012) 576–581. [16] S. Harries, J. HEIMANN, J. Hinnenthal, Pareto optimal routing of ships,
CE
985
transformation 772 (2003) 1.
AC
[17] S. Marie, E. Courteille, et al., Multi-objective optimization of motor vessel
990
route, in: Proceedings of the Int. Symp. TransNav, Vol. 9, 2009, pp. 411– 418.
[18] A. Dec` o, D. M. Frangopol, Real-time risk of ship structures integrating structural health monitoring data: Application to multi-objective optimal ship routing, Ocean Engineering 96 (2015) 312–329. 50
ACCEPTED MANUSCRIPT
[19] J. Szlapczynska, R. Smierzchalski, Multicriteria optimisation in weather routing, Marine Navigation and Safety of Sea Transportation (2009) 423. 995
[20] M.-C. Tsou, Integration of a geographic information system and evolution-
Navigation 63 (02) (2010) 323–341.
CR IP T
ary computation for automatic routing in coastal navigation, Journal of
[21] S. Bijlsma, A computational method in ship routing using the concept of limited manoeuvrability, Journal of Navigation 57 (03) (2004) 357–369. 1000
[22] M.-C. Fang, Y.-H. Lin, The optimization of ship weather-routing algorithm
AN US
based on the composite influence of multi-dynamic elements (ii): Optimized routings, Applied Ocean Research 50 (2015) 130–140.
[23] A. Orda, R. Rom, Minimum weight paths in time-dependent networks, Networks 21 (3) (1991) 295–319. 1005
[24] A. Orda, R. Rom, Shortest-path and minimum-delay algorithms in net-
ED
37 (3) (1990) 607–625.
M
works with time-dependent edge-length, Journal of the ACM (JACM)
[25] D. E. Kaufman, R. L. Smith, Fastest paths in time-dependent networks for intelligent vehicle-highway systems application, Journal of Intelligent Transportation Systems 1 (1) (1993) 1–11.
PT
1010
[26] M. M. Kostreva, M. M. Wiecek, Time dependency in multiple objective
CE
dynamic programming, Journal of Mathematical Analysis and Applications 173 (1) (1993) 289–307.
AC
[27] T. Getachew, M. Kostreva, L. Lancaster, A generalization of dynamic
1015
programming for pareto optimization in dynamic networks, RAIROOperations Research 34 (01) (2000) 27–47.
[28] H. W. Hamacher, S. Ruzika, S. A. Tjandra, Algorithms for time-dependent bicriteria shortest path problems, Discrete optimization 3 (3) (2006) 238– 254. 51
ACCEPTED MANUSCRIPT
1020
[29] A. Makrygiorgos, I. A. Vetsikas, S. Perantonis, Accelerating multi-objective ship routing using a novel grid structure and a simple heuristic, in: Proceedings of the 6th International Conference on Information, Intelligence,
CR IP T
Systems and Applications, IEEE, 2015. [30] N. Touati, V. Jost, On green routing and scheduling problem, arXiv preprint 1203.1604.
1025
[31] M. Christiansen, K. Fagerholt, D. Ronen, Ship routing and scheduling: Status and perspectives, Transportation science 38 (1) (2004) 1–18.
AN US
[32] S. Bijlsma, A computational method for the solution of optimal control problems in ship routing, Navigation 48 (3) (2001) 145–154. 1030
[33] O. Bolza, Vorlesungen u ¨ber Variationsrechnung, BoD–Books on Demand, 2013.
M
[34] R. W. James, Application of wave forecasts to marine navigation. [35] K.-Y. Chang, G. E. Jan, I. Parberry, A method for searching optimal routes
ED
with collision avoidance on raster charts, The Journal of Navigation 56 (03) (2003) 371–384.
1035
[36] S. E. Dreyfus, An appraisal of some shortest-path algorithms, Operations
PT
research 17 (3) (1969) 395–412. [37] A. Warner-Grieve, Optimal route finding with collision avoidance in marine
CE
navigation, Tech. rep. (2006).
[38] A. Veneti, C. Konstantopoulos, G. Pantziou, An evolutionary approach to
AC
1040
multi-objective ship weather routing, in: Information, Intelligence, Systems and Applications (IISA), 2015 6th International Conference on, IEEE, 2015, pp. 1–6.
[39] A. N. Perakis, N. A. Papadakis, Minimal time vessel routing in a time-
1045
dependent environment, Transportation Science 23 (4) (1989) 266–276.
52
ACCEPTED MANUSCRIPT
[40] H. Lee, G. Kong, S. Kim, C. Kim, J. Lee, Optimum ship routing and its implementation on the web, in: Advanced Internet Services and Applications, Springer, 2002, pp. 125–136.
CR IP T
[41] J. Bisschop, AIMMS-Optimization modeling, Paragon Decision Technology, 2011.
1050
[42] L. Mandow, J. L. P. De La Cruz, Multiobjective a* search with consistent heuristics, Journal of the ACM (JACM) 57 (5) (2010) 27.
[43] I. Koromila, Z. Nivolianitou, T. Giannakopoulos, Bayesian network to pre-
AN US
dict environmental risk of a possible ship accident, in: Proceedings of the
7th International Conference on PErvasive Technologies Related to Assis-
1055
tive Environments, ACM, 2014, p. 44.
[44] S. Wei, P. Zhou, 23. development of a 3d dynamic programming method for weather routing, Methods and Algorithms in Navigation: Marine Nav-
1060
M
igation and Safety of Sea Transportation (2011) 181.
[45] Y. Kwon, Speed loss due to added resistance in wind and waves, Naval
ED
Architect (2008) 14–16.
[46] K. Avgouleas, Optimal ship routing, Ph.D. thesis, Massachusetts Institute
PT
of Technology (2008).
[47] I. MSC, 83/inf. 2: Formal safety assessmentconsolidated text of the guidelines for formal safety assessment (fsa) for use in the imo rule-making pro-
CE
1065
cess (msc/circ. 1023- mepc/circ. 392), International Maritime Organisation:
AC
London, UK.
[48] Z. S. Nivolianitou, I. A. Koromila, T. Giannakopoulos, Bayesian network to
1070
predict environmental risk of a possible ship accident, International Journal of Risk Assessment and Management 19 (3) (2016) 228–239.
[49] C. De Boor, A practical guide to splines, Vol. 27, Springer-Verlag New York, 1978. 53
ACCEPTED MANUSCRIPT
[50] E. Guti´errez, A. L. Medaglia, Labeling algorithm for the shortest path problem with turn prohibitions with application to large-scale road networks, Annals of Operations Research 157 (1) (2008) 169–182.
1075
AC
CE
PT
ED
M
AN US
putation Theory, Springer, 1993, pp. 234–243.
CR IP T
[51] P. Fischer, Finding maximum convex polygons, in: Fundamentals of Com-
54