International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 3
Number 2
June 2009
133
1 INTRODUCTION
Ship weather routing develops an optimum track for
ocean voyages based on forecasts of weather, sea
conditions, and ship's characteristics. Within
specified limits of weather and sea conditions, the
term optimum is used to mean maximum safety and
crew confort, minimum fuel and oil consumption,
minimum time or distance underway, or any desired
combination of these factors.
There exists many works dealing with the
problem of optimal ship weather routing, in which
motor or sailing vessels are considered. The
approaches in the literature may be divided into four
categories : isochrone construction, application of
the calculus of variation, dynamic programming and
evolutionary algorithms.
James (1957) proposed a scheme for non-variable
weather conditions using lines of equal-time, or
isochrones to achieve minimal-time objective. This
method is very efficient for solving the problem of
determinist minimal-time weather routing but is not
adapted for minimization of the vessel's fuel
consumption. Numerous variation of this type of
method have been presented. For exemple Hagiwara
& Spaans (1987) presented an improvement of
isochrone definition for sail-assisted motor vessel
routing. In addition to time objective, the author
tried to minimize fuel consumption but the method
is not very efficient because the propeller speed is
kept constant on the track.
Haltiner, Hamilton and Árnason (1962) solved
the same problem by using the calculus of
variations. This method is based on parametric
curves obtained by solving the associated Euler
differential equation by relaxation methods. Bleick
& Faulkner (1965) extended this approach to the
case of deterministic variations of sea state. This
definition is mathematicaly elegant but impractical
for ship weather routing because of convergence
problems.
The third approach was proposed by Zappoli
(1972), who treated the minimal time problem as a
decision process solvable using dynamic
programming. The sailing domain is discretized
using grid refinement techniques. Allsopp & al.
(2000), modified this method integrating branching
scenario structure to model the manner the weather
will evolve in time. The main advantage of that kind
of method is that the problem is divided into a set of
linked stages and the optimal decision depends on
decisions made in the previous stages. But, for fine
grid, the calculation time may be very high and the
amount of data very large.
The most recent approaches use B-spline technics
and evolutionary algorithms in order to minimize
fuel consumption and maximize safety. Harries,
Heimann and Hinnenthal (2003) proposed such a
method for large motor vessels using Multi-
Objective Genetic Algorithms. Hinnethal and
Saertra (2005) improved this method considering the
stochastic nature of weather along the route. Böttner
(2007) recently used this work combined with
Multi-Objective Optimization of Motor Vessel
Route
S. Marie & E. Courteille
Department of Mechanical Engineering and Automation, Institut National des
Sciences Appliquées,
Rennes, France
ABSTRACT: This paper presents an original method that allows computation of the optimal route of a motor
vessel by minimizing its fuel consumption. The proposed method is based on a new and efficient meshing
procedure that is used to define a set of possible routes. A consumption prediction tool has been developed in
order to estimate the fuel consumption along a given trajectory. The consumption model involves the effects
of the meteorological conditions, the shape of the hull and the power train characteristics. Pareto-optimization
with a Multi-Objective Genetic Algorithm (MOGA) is taken as a framework for the definition and the
solution of the multi-objective optimization problem addressed. The final goal of this study is to provide a
decision helping tool giving the route that minimizes the fuel consumption in a limited or optimum time.
134
dynamic programming for costal approach in order
to propose the best possible track from harbor to
harbor. These methods are characterized by a low
number of free variables to describe both the course
and the speed of the boat. Moreover, a high number
of route variants can be considered from which
Pareto optimal solutions may be identified.
In the four different approaches, the scheme for
optimisation is almost the same :
− Mathematical modeling of the ship to compute
the objectives,
− 3-D interpolation (time and space) of weather and
sea state data,
− Parametric route definition,
− Optimization of the route using an algorithm.
Our work is based on a new discretization of the
research space based on few physical parameters.
This parametric definition of the griding makes it
understandable and easy to tune. Moreover this kind
of meshing may be applicable for all type of vessels,
journeys and all weather conditions are easily taken
into account. The gridding of the sail area is
systematic and uniform since it is based on spherical
geometry and accepts constrains like bathymetric
data. Our work is related to the method proposed by
Harries & al. because we kept their definition of
both the geography and speed since it allows the
complete location of the boat in space and time. The
modeling of motor vessels is not the purpose of this
work and will just be briefly presented. Future works
concern the identification of the vessel model using
fuzzy logic technique in order to obtain an accurate
consumption model. Since this identification is not
achieved yet, a model from the literature will be
used to present this new gridding method.
The meshing of the explored area is presented
first. Then, the way of constructing routes is shown
and a sensitivity of the meshing is presented. Next,
the way of modeling a motor vessel is introduced.
Finally, results of numerical optimizations are
presented in order to evaluate the limitations and
benefits of the proposed meshing method.
2 MESHING OF THE EXPLORED AREA
The new automatic meshing method that we propose
is based on spherical rhombus where two of the
opposite vertexes are the departure and the arrival
points. The main advantages of this discretization of
the sailing area are :
− The genericity of its construction taking into
account the sea-beds geography, the time
dependant meteorological data and the
characteristics of the vessel.
− The systematic gridding of the explored area with
few physical parameters.
− The automation of its calculation leading to
optimizable routes.
− The possible reactualization of the rhombus to
change the routing policy during the sailing.
2.1 Rhombus definition
In this part denotes the departure point, the
arrival one and is the center of Earth considered
spherical. (
) is the plane containing the lines
() and () carrying respectively vectors and
. is the angle between these two vectors and is
their bisectrix. is one vector normal both to and
. We also define the two related unit vectors
and
. Let be (
) the plane containing , , the two
remain vertexes of the rhombus and and point
. The lines () and () are directed by vectors
and . We defined
(
,
)
= 2. This notation is
recalled in the figure 1.
Figure 1. Main planes to define the meshing
Point denotes the intersection between the great
circles
and
. (
) is the plane containing
and , we have
(
,
)
= . This angle is the image of
the orthodromic distance between and .
Knowing the maximal speed of the vessel
,
imposed by the design of the hull and the power
train characteristics, and the desired time of sailing
, the maximal distance that can cross the vessel
during the time window is :
=
.
By this mean we can set the maximal distance on
the great circle route :
=
, (1)
135
where is an dimensionless factor lower than the
unit used to get some margin. From this equation
(1), angle is calculable as follow:
=
, (2)
with the Earth mean radius. To compute the value
of angle, the following vectorial equations are
used:
=
, (3)
=
+
, (4)
. = .
(
)
. (5)
Developing this relation (5), one can write:
= cos
, (6)
= tan
= sin
. (7)
Using the inverse spherical transformation (7)
and the definition of and vectors (3,4), it is
possible to compute the spherical coordinates of
and :
= tan
, (8)
= sin
cos
sin
. (9)
Where
is the longitude of point and
denotes its latitude. The same kind of relations can
also be written for point .
2.2 Meshing's levels calculation
The previous construction is extended in order to
define the
levels of the meshing. For that
purpose, planes
are defined,
[
1,
2
]
.
These planes contain the vectors and
(Fig. 2).
The components of vector
are calculated
according to:
.
= cos
.
= 0
. = cos
. (10)
By solving the system (10) the components of
are computed. To complete the spherical pyramid
(Fig. 2), three more plans comparable to (
) are
defined : (
) contains and , (
) contains and
, (
) contains and .
Figure 2. Definition of planes
For each of the above mention planes, their
normal
,
,
,
are defined,
the normal to
(
) is also calculated. The intersections between
(
), (
) and (
) (=
{
3,5
}
) are the two lines
called (
) and (
) oriented respectively by
vectors
and
. The components of these
vectors are calculated according to :
− If
(
,
)
,
then
=
=
− If
(
,
)
>
,
then
=
=
.
The related unit vectors
and
are
defined. In each level of the meshing i.e. plane
(Fig. 3), nodes of the level are defined. The distance
between 2 nodes is constant and defined by the
parameter
. As a result, in a level,
nodes
,
are defined with :
=
+ 1. (11)
Figure 3. Vectors and nodes in plane
The angle
between vectors
and
is
defined by :
= cos
.
. (12)
136
For each node in the level, directing vector
,
of
the line (
,
) is defined by the following system :
,
,
= cos
,
,
= cos
,
,
= 0
. (13)
So the cartesian coordinates of the
node of the
level are :
,
=
[
1,
2
]
1,
. (14)
Using the inverse spherical transformation (7),
the coordinates of
,
(14) are expressed within the
spherical coordinate system.
2.3 Limitation of the possible nodes
From vectorial cartographies, the matrix giving the
depth of water according to longitude and latitude of
meshing's points is compared with the draught of the
boat. The nodes at which the depth is insufficient are
removed from the grid (14).
A criterion of maximum course between two
successive nodes is also defined. This criterion
makes it possible to exclude the nodes which move
the ship too away from its final destination .
3 ROUTE DEFINITION
3.1 Geographical definition of routes
3.1.1 Definition of routes
Each possible route is defined using a navigable
node of each level of the meshing. The nodes
(,
, . . .
, ) define the control points of the
associated Bézier curve. We choose Bézier curve
since it begins at and end at in addition the
curve is always in the convex hull of the control
polygon. For recall the Bézier curve is defined by :
(
)
(
)
,
[
0,1
]
.
(
)
(15)
Where
(
)
(
)
stand for the Bernstein basis
polynomials.
3.1.2 Discretization of Bézier curve
Moreover because of its parametric definition the
Bézier curve is discreditable. This discretization is
done relatively to a parameter
corresponding to
the maximal number of course changes per hour.
The number of segments of a route
is defined
from the minimal distance i.e. the orthodromic
distance
sailed at the maximum speed of the
vessel
:
=
. (16)
As a result, the points
defining the route are
computed as presented below :
(
)
(
)
(
)
, =
, = {0, . . . ,
}. (17)
The discretization of Bézier curves is done
because on each facet of the route, we consider that
both the weather and sea state remain constant. It
also allows to define loxodromic courses between
and
which leads to a route defined by
waypoints and courses.
must be set with great
care because it corresponds to a compromise
between the number of course changes that the
captain has to perform and the approximation of the
weather field along the course.
3.2 Velocity along a route
In order to locate the vessel both in time and space,
the velocity on the course is set. As a result along
each facet of the route, the time dependant weather
data are known. Moreover the crossing time
is easily calculable.
The target speeds of the vessel
are included
between two boundaries :
and
.
has
to be tuned by the captain. The number of target
speeds is
. Each target speed is valid on several
segments of the discretized route. The number of
facets
on which
is valid is defined by :
=
+ 1. (18)
3.3 Meteorological conditions along the route
The weather and the sea state at the current position
of the vessel are extracted from GRIB files defined
with a regular 1.25
× 1
grid downloaded from
NOAA ftp
1
. These files gathered the meteorological
data for a 180 hours time window with a 3 hours
step.
The decoding of these files allows the
construction of true wind, current and waves fields
for the sailing time window. The parameters of these
fields are presented in Table 1.
1
://. . . ///_/
137
Table 1: Parameters of weather and sea state
___________________________________________________
True wind speed .
True wind direction
Current speed .
Current direction
Swell height
Swell period
Swell direction
___________________________________________________
As the meteorological grid weather does not
correspond to the route's points
, a space
interpolation of the data is done. For that purpose,
we used 2D linear interpolation techniques to
estimate the encountered conditions during the
whole time window for each
. In addition, the
time dependency of the fields can be easily taken
into account by interpolation in time.
4 MESHING SENSITIVITY TO ITS
PARAMETERS
The main advantage of the method compared to the
previous approaches is that the rhombus definition is
based on physical parameters (Tab. 2) easily
adjusted and tuned by the captain. The journey is
defined by
and
. The design of the boat
imposes
and the four other meshing's
parameters define its geometry. On Fig. 4, two
different values of each adjustable parameters are
shown. The dots corresponds to the navigable nodes
and the cross to the non navigable ones These
parameters are left in user's hand but a great
attention must be brought to their value since they
define the manner the research space is discretized
and then the fineness of the solution.
Table 2: Physical parameters defining the meshing
___________________________________________________
Departure point °, °
Arrival point °, °
Maximum speed of the vessel .
Objective time
Number of levels
Distance between nodes
Course changes per hour
___________________________________________________
Moreover the discretization of the research space
is uniform because of the spherical definition of the
meshing. This property is interesting for long
voyages since the routes will be equi-distributed on
Earth's surface. Of course for short travels, the
Earth's roundness is too small to have an influence
on the meshing definition so a planar definition of
the rhombus is more convenient.
Figure 4. Sensitivity to meshing's parameters
5 CONSUMPTION MODELING
At this point, the set of possible routes has been
defined and the physical magnitudes have been
calculated. In the perspective of optimizing the route
of the vessel, a cost function has to be defined in
order to select the best track.
In the literature, various approaches are proposed.
For example Zappoli (1972) chose the sailing time
as cost function, Hinnenthal & al. (2003) chose both
the consumption and the estimate time of arrival.
But to compute these functions, the ship
performances in a seaway must be accurately known
Journée & Mejiers (1980). Most of the works use
parametric models in order to estimate the boat's
behavior. The scheme of calculation in the recent
approaches is often the same one :
− Estimation of the resistances acting on the hull :
138
− Still water resistance based on regression
analysis of model tests and full-scale data
(Haltrop, Van Oomertsen, ...),
− Wind resistance of the emmerged part of the
vessel (Isherwood),
− Added resitance due to waves (Gerristma,
Boese,...),
− Operating point of the propeller with the ITTC
power prediction method,
− Operating point of the main engine to estimate
the consumption.
These parametric models are well known and
controlled but their identification requires specific
equipments such has wind tunnels or towing tanks,
realization of instrumental models, instrumentation
of the ship for full-scale measurements. The
calculation time may be rather important and the
measurement cost very high to get an efficient
estimator of performances so applications of these
methods are limited. Moreover the parameters are
time dependent if we take into account the fouling of
the hull leading to an increase of resistances for
example. In addition, these models are not generic,
the calculation method must be adapted to each
hull's shape.
To analyze the efficiency of the proposed
meshing method, we introduce a parametric model
that will react to the encounter sea conditions. We
are aware that such a raw model could not describe
finely the resistances that encounter a boat in a real
seaway but it as to be implanted to test the general
approach we propose. The consumption model will
further be identified on board using fuzzy logic
technics. As the meshing method is uncoupled from
the modeling, any function cost can be used to
estimate the efficiency of a particular route. We used
the ship presented in Hagiwara (1987) since many
information are available and the time calculation
cost is low (3). For our routing scheme, we chose
the consumption and the sailing time as estimators
of the vessel efficiency.
5.1 Scheme of calculation
The scheme of consumption calculation is presented
on the figure 5. As presented earlier, both the
geography and the velocity are defined for a route,
as a result, the inputs of the consumption calculation
are the GPS position of the boat and the desired
seabed velocities
.
Figure 5. Scheme calculation of the consumption
Classically, the composition between true wind
and seabed velocity of the vessel leads to the
apparent wind. In order to get the desired seabed
velocity, the current effects must be compensated if
the resulting relative speed does not exceed the
limitation due to the vessel design i.e. the maximal
speed of the vessel
. The drifting of the vessel is
neglected here.
5.2 Resistances computation
In this part the relations used to compute the
resistances acting on the hull are presented rapidly.
They will allow to estimate the power that has to be
delivered by the vessel's engine.
− Still water resistance
is the first resistance
acting on the vessel. It corresponds to the energy
used to overcome the frictional resistance of the
hull plus the one used to create the bow wave.
− Aerodynamic resistance
is the action of the
wind onto the emerged part of the vessel.
− Added resistance due to waves
is the
increase of resistance due to the encountered
waves. It depends of their average period,
significative height and mean direction.
The total resistance
is the global resistance
acting on the hull. Its value is calculated according
to:
=
+
+
. (19)
The numerical values of the model parameters
and the vessel design, are presented in Hagiwara &
Spaans (1987).
139
5.3 Propulsion characteristics
As the resistances acting on the hull are known, the
propulsion system has to be modeled in order to
evaluate the torque and power that must be delivered
by the engine at the target velocity
. For the
routing exemple we chose a 8 diameter fixed pitch
Wageningen B5-75 screw series propeller, Carlton
(2007) and a 10000 Wärtsilä engine (2007).
Propeller thrust, torque, rate of revolution The
propeller thrust
, torque
and power
are
calculated using the ITTC scheme of calculation
(1978). Their calculation are well known and will
not be discussed here. The propeller's rate of
revolution is adjusted to obtain the proper thrust
.
Engine power and consumption For a known
propeller and engine, a fixed propeller revolution
rate and a given ship speed, the engine power can be
calculated as follows:
=
, (20)
with
=
the efficiency of the hull and the
thrust deduction fraction due to suction of the water
in front of the propeller, the wake fraction,
the
open water efficiency,
the relative rotative
efficiency and
the mechanical efficiency of shaft
bearing.
The consumption of the engine is calculated
knowing the delivered power. For that purpose, we
use the specific consumption law
given by the
engine manufacturer. For a given rate of revolution
and power
of the engine, the hourly
consumption is given by :
=
(
). (21)
The unit of
is .
.
6 OPTIMIZATION SCHEME
6.1 Search method
6.1.1 Performances indices
The goal of the optimization is to minimize the
antagonist objectives : the consumption and the
sailing time . A numerical optimization of a route
off Corea is presented hereafter (Fig. 6). The journey
is defined between
= (133
40
) and
= (122. 5
32. 5
). The number of level is
= 7, the distance between nodes of a level is
= 25 and the number of course changes per
hour is
= 1
. The number of target speeds is
= 8. For this application we used the
meteorological data of the 23
April 2008 at 0: 00
GMT which will be the departure time.
6.1.2 Pareto-optimal solutions
Solving this optimization problem with
conflicting objectives across a high-dimensional
research space is a difficult goal. Instead of a single
optimum, there is rather a set of alternative trade-
offs, generally known as Pareto-optimal solutions.
Various evolutionary approaches to multi-objective
optimization have been proposed since 1985,
capable of searching for multiple Pareto-optimal
solutions concurrently in a single simulation run ,
Valdhuizen & Lamont (2000). The optimization
program FRONTIER
®2
and the technical computing
software MATLAB
®
are used to set up the
framework of the multi-objective design
optimization study of weather routing. The Multi-
objective Genetic Algorithm (MOGA), implemented
first by Fonseca & Fleming (1998), is used to
perform the optimization problem.
6.1.3 Design parameters
The number of parameters necessary to define a
route is
(
2 +
)
. For each level of the
meshing the associate parameter is the index of the
node
,
. Concerning the target speeds, the
parameters are within the boundary previously
presented and their step is 0.1.
6.1.4 Global optimization process
The algorithm will attempt a number of
evaluations equal to the size of the initial population
for the MOGA multiplied by the number of
generation. The initial population is generated by a
random sequence of 60 designs. The major
disadvantage of the MOGA is mainly related to the
number of evaluations necessary to obtain
satisfactory solutions. The search for the optimal
solutions extends in all the directions from design
space and produces a rich data base and there is not
a true stop criterion. The numerical evaluation of the
performances calls upon MATLAB codes is not so
expensive in terms of computing time (about 2 ). In
an attempt to solve the optimization problem in an
acceptable timeframe, the number of generations
evaluated is almost 70, i.e. 4000 designs in all. The
required computation time for the global
optimization process is about 2 hours (2.4 GHz / 3.0
Gb RAM). Integrating a Response Surface
Methodology to reduce the computation time could
be an interesting extension of our work especially if
one wants to achieve on board routing.
2
http://www.esteco.com/
140
6.2 Numerical optimizations
The time dependency of the sea state and wind field
is taken into account with linear interpolation in
time. Figure 6 presents an exemple of weather
routing. On these maps, two positions of the vessel
are shown. For the analysis of the optimization, we
compute the shortest path and we tune the constant
velocity on this path so that the sailing time is the
same as the optimal one. The comparison between
these two routes is done in Table 3.
Figure 6. Comparison between shortest and optimal routes in
wave field
Table 3: Comparison of the optimum route to the shortest one
___________________________________________________
Route Distance () Time () Fuel ()
________________________________________________
Shortest 1337.7 56.8 198.4
Optimal 1342.1 56.8 188.4
________________________________________________
This example puts forward the interest of the
routing weather for motor vessel. We have shown
that a longer distance of journey does not imply
automatically an increase of consumption. The
example above shows that a saving of 5% is
realizable on the chosen journey. The economy are
not very large because the distance to be traversed is
too short. Moreover, the sea is not rough, the
maximum height of waves during the journey is only
2. Another problem that might occur is the spatial
interpolation of the conditions. In the NOAA's
GRIB, on earth, the fields are not known so an
arbitrary value is set (for instance 0). This
unknowing of the fields leads to inconsistent value
of the advancing resistance along the coast that
might distorted the results.
7 CONCLUSION & FUTURE WORK
This paper presents a brief outlook to motor vessel
routing using deterministic weather forecasts. A
method for spatial and temporal generation of route
variants based on a generic and automatic meshing
method has been presented. The major advantage of
this technic is the physic based definition of the
rhombus and the low number of free variables used
to define a route. The ship route was optimized
using multi-objective algorithm for a deterministic
weather case. The reasonable computation time will
allow the use of this method for on board routing
applications. Pareto-optimization may be considered
as a tool providing a set of efficient solutions among
different and conflicting objectives, under different
constraints. The final choice remains always
subjective and is let to the user's hands.
Future works concern the improvement of the
consumption prediction model using on board
measurements and fuzzy logic identification
technics. We will also investigate the routing of
wail-assisted motor vessels. The use of these
technics will allow the quick establishment of
reliable models of different types of vessels keeping
the measurement cost very low since no tawing tank,
wind-tunnel tests or models will be necessary. An
other aspect that has not been discussed in this paper
is the stochastic nature of weather and sea conditions
that must been taken into account. This aspect will
also be investigated using robust optimization
technics. The reactualization of the routing policy
during the journey will also be considered for long
travels. Thus we will be able to integrate the most
recent weather data.
ACKNOWLEDGEMENTS
The research in this paper is supported by the Région
Bretagne and some local community. This work
takes place in the Grand Largue project.
REFERENCES
Allsopp, T. 2000. Optimising Yacht Routes Under Uncertainty.
Proceedings of the 2000 Fall National Conference of the
Operations Reaseach Society of Japan : 176-183.
141
Bleick, W. & Faulkner F. 1965. Minimum-Time Ship Routing.
Journal of Applied Meteorology, 4 : 217-221.
Braddock, R.D. 1970. On Meteorological Navigation. Journal
of Applied Meteorology 9(1) : 149-153.
Böttner, C.U. 2007. Weather Routing for Ships in Degraded
Condition. International Symposium on Maritime Safety,
Security and Environmental Protection, Athens, Greece.
Carlton, J. 2007. Marine Propellers and Propulsion. Butter-
worth-Heinemann.
Fonseca, C.M. & Fleming, P.J. 1998. Multiobjective Optimiza-
tion and Multiple Constraint Handling with Evolutionary
Algorithms. IEEE Trans. On Systems, Man and Cybernetics
28 : 26-37.
Hagiwara, H. & Spaans, J. 1987. Practical Weather Routing of
Sail-Assisted Motor Vessels. Journal of Navigation 40(1) :
96-119.
Haltiner, G.J., Hamilton, H.D. & Árnason, G. 1962. Minimal-
Time Ship Routing. Journal of Applied Meterology 1(1) : 1-
7.
Harries, S., Heinmann, J. & Hinnenthal J. 2003. Pareto-
Optimal Routing of Ships. International Conference on
Ship and Shipping Research.
Hinnentham, J. & Saerta, Ø. 2005. Robust Pareto-Optimal
Routing of Ships Utilizing Ensemble Weather Forecasts.
Maritime Transportation and Exploitation of Ocean and
Coastal Resources : 1045-1050.
ITTC. 1978. 1978 ITTC Performance Prediction Method for
Single Screw Ships. Proceedings of the 15
th
International
Towing Tank Conference, 1978, The Hague, Netherland.
James, R.W. 1957. Application of Wave Forecasts to Marine
Navigation. U.S. Navy Hydrographic Office.
Journée, J.M.J. & Mejijers, J.H.C. 1980. Ship Routeing for Op-
timum Performance. IME Transactions : 1-17.
Wärtsilä. 2007. Project Guide Wärtsilä 46.
Valdhizen, D.A.V. & Lamont, G.B. 2000. Multiobjective Evo-
lutionary Algorithms : Analyzing the State-of-the-art. Evo-
lutionary Computation 8 : 125-147.
Zappoli, R. 1972. Minimum-Time Routing as an N-Stage De-
cision Process. Journal of Applied Meteorology 11(3) : 429-
435.
