Document [original]

Robust Pareto – Optimum Routing of Ships

utilizing

Deterministic and Ensemble Weather Forecasts

vorgelegt von

Diplom-Ingenieur

Jörn Hinnenthal

aus Berlin

von der Fakultät V – Verkehrs- und Maschinensysteme

der Technischen Universität Berlin

zur Erlangung des akademischen Grades

Doktor der Ingenieurwissenschaften

– Dr.-Ing. –

genehmigte Dissertation

Promotionsausschuss:

Vorsitzender: Prof. Dr.-Ing. Jürgen Siegmann

Berichter: Prof. Dr.-Ing. Günther F. Clauss

Berichter: Prof. Dr.-Ing. Apostolos D. Papanikolaou

Tag der wissenschaftlichen Aussprache: 18.12.2007

Berlin 2008

D 83

Pareto Optimum Ship Routing

Acknowledgement

Inspired by my participation as student assistant in the European research project

SEAROUTES, I decided to continue the research into ship route optimization that I already

started within my diploma thesis. In this regard I owe my deepest gratitude to my promoters

Dr.-Ing. Stefan Harries and Prof. Dr.-Ing. Lothar Birk, who encouraged me to become

scientific assistant and to continue with this work up to a doctoral thesis.

I want to express my gratitude to my promoter and supervisor Prof. Dr.-Ing. Günther Clauss

for his encouragement to publish scientific results, for the outstanding composition of

guidance and freedom, and for generous support whenever I asked for. I also want to express

my gratitude to my second supervisor Prof. Dr.-Ing. Apostolos Papanikolaou for his

instantaneous agreement to contribute in the doctoral committee, and I only regret that I

didn’t involve him much earlier into this work. Also many thanks to the chairman of the

doctoral committee Prof. Dr.-Ing. Jürgen Siegmann for supporting a fast and frictionless

procedure.

Special thanks I owe to my former colleagues Claus Abt, Dr.-Ing. Justus Heimann, and

Henning Winter for reliable support and many fruitful discussions. The same applies to my

present, esteemed, and embosomed colleagues Felix Fliege, Gonzalo Tampier Brockhaus and

Dr.-Ing. Uwe Boettner who built for a long time the core of the workgroup at Technische

Universität Berlin for me, the fundament for the efficiency and the pleasure within my work.

Further on I would like to thank all my colleagues that I not mentioned by name, because they

not directly contributed to this work, I got a lot of support in other projects and found

wonderful companions.

Sincere thanks are given to Dr. Øvind Saetra whom I met within the SEAROUTES project,

who provided me with exquisite weather forecasts, and who inspired me to apply ensemble

forecasts for the routing problem. Sincere thanks are also given to Thor Marquardt, chief

officer at Hapag Lloyd, for his valuable support and consultancy regarding navigation and

operation of ships and basic conditions of maritime transport, and to Dr. Masaru Tsujimoto

from the National Maritime Research Institute in Tokyo, who became an esteemed colleague

during his time in Berlin and who inspired me a lot within our lively discussions on routing.

Special thanks also to Gabriele Schmitz for her excellent high-speed proofreading.

Finally and most of all I would like to express my gratitude to Nicole Reimer, my beloved

wife, who encourages me and believes in me and our partnership, also in times that are

everything else than mellifluous.

Jörn Hinnenthal – Berlin, March 2008

Pareto Optimum Ship Routing

Abstract

Sophisticated routing of ships is increasingly recognized as an important contribution to safe,

reliable, and economic ship operation. The more reliable weather forecasts and performance

simulation of ships in a seaway become, the better they serve to identify the best possible

route in terms of criteria like: ETA (estimated time of arrival), fuel consumption, safety (of

ship, crew, passengers. and cargo), and comfort. This establishes a multi-objective, non-

linear, and constrained optimization problem in which a suitable compromise is to be found

between opposing targets.

For its solution a new optimization approach to select the most advantageous route on the

basis of hydrodynamic simulation and sophisticated weather forecast is posed. Transfer

functions are employed to assess the operating behavior of a ship in waves. Probabilistic

ensemble forecasts, provided by the European Centre for Medium-Range Weather Forecasts

ECMWF, are applied to account for the stochastic behavior of weather. The routes in adverse

weather conditions are established as perturbations of a parent route in calm weather, which is

assumed to be the concatenated great circles between waypoints. Utilizing a B-spline

technique, the number of free variables for describing both the course and the velocity profile

is kept low. For solving the multi-objective, non-linear, and constrained optimization problem

the commercial package modeFRONTIER is successfully applied. In order to balance

opposing criteria and to value the performance of optimized routes an approach suggested by

the Italian economist Vilfredo PARETO is adopted. A multi-objective genetic algorithm turns

out to be a suitable optimization method to identify PARETO optimum routes for a

sustainable support of a conscious decision-making process.

An elaborated example is given for an intercontinental container service, employing the

Panmax container vessel CMS HANNOVER EXPRESS, between Europe and North

America. Different weather situations for the North Atlantic are taken into account.

Sensitivity studies are applied to validate the set-up and to conduct plausibility tests. The

robustness of optimized routes against weather changes, time loss, fuel consumption,

accelerations, slamming, and parametric rolling are taken into account.

Pareto Optimum Ship Routing

Kurzdarstellung

Zur Gewährleistung eines sicheren, zuverlässigen und wirtschaftlichen Schiffsbetriebes

kommen zunehmend Routing-Systeme zur Anwendung. Um sichere Aussagen über die

Reisedauer, den Brennstoffverbrauch, Sicherheit oder auch den Komfort an Bord treffen zu

können, spielen die Zuverlässigkeit von Wettervorhersagen und die der Simulation der

Seegangseigenschaften des Schiffes eine große Rolle. Dieses führt letztlich zu einem nicht-

linearen Optimierungsproblem mit mehreren, oftmals antagonistischen Gütekriterien. Der

Lösungsraum ist zudem durch Nebenbedingungen beschränkt.

Zur Lösung dieser Aufgabe wird ein neu entwickelter Ansatz vorgestellt. Auf Basis von

Bewegungssimulationen des Schiffes im Seegang und detaillierten Seegangsprognosen

werden optimale Kurse und zugehörige Geschwindigkeitsprofile eines Schiffes identifiziert.

Das Betriebsverhalten des Schiffes im Seegang wird dabei mit Hilfe von

Übertragungsfunktionen und Seegangsspektren bestimmt. Zunächst kommen deterministische

Seegangsvorhersagen zum Einsatz; später werden diese durch Ensemble-Vorhersagen

erweitert, um das stochastische Verhalten möglicher Wetterentwicklungen abzubilden. Die

Seegangsprognosen werden vom European Centre for Medium-Range Weather Forecasts

ECMWF zur Verfügung gestellt.

Zur Beschreibung von Routenvarianten wird eine Perturbationsmethode verwendet. Sie

basiert auf der Darstellung des Kurses und eines zugehörigen Geschwindigkeitsprofils durch

B-Splines. Für die Lösung der Optimierungsaufgabe kommt die generische Optimierungs-

Software modeFRONTIER zum Einsatz. Zur Bewertung optimierter Routen unter dem

Aspekt sich widersprechender Gütekriterien, wird ein von dem italienischen Ökonomen

Vilfredo PARETO vorgeschlagenes Konzept verwendet. Ein heuristisches Suchverfahren, der

Genetische Algorithmus für mehrere Gütekriterien, zeigt gute Ergebnisse bei der Suche nach

PARETO optimalen Routen, die zur Unterstützung der Routenplanung herangezogen werden

können.

Die neu entwickelte Optimierungsmethode wird am Beispiel eines Containerservices mit der

CMS HANNOVER EXPRESS im Nord Atlantik erläutert. Die Auswirkungen verschiedener

Wetterszenarien auf ein Optimierungsergebnis werden dargestellt. Empfindlichkeits- und

Variationsstudien dienen der Validierung und zur Plausibilitätskontrolle.

Pareto Optimum Ship Routing

Content

Acknowledgement ........................................................................................................i

Abstract .......................................................................................................................ii

Kurzdarstellung........................................................................................................... iii

Content.......................................................................................................................iv

List of Figures.............................................................................................................vi

List of Tables ...............................................................................................................x

Nomenclature .............................................................................................................xi

1 Introduction to weather routing of ships............................................................... 1

1.1 State of the art in applied research ............................................................................................2

1.2 Routing services......................................................................................................................... 8

1.3 A new approach to route optimization...................................................................................... 10

2 Basic principles and terms of optimization......................................................... 12

2.1 SIMPLEX algorithm .................................................................................................................. 14

2.2 Genetic algorithms, GA and MOGA ......................................................................................... 15

3 Modeling of ship route and environmental conditions........................................ 19

3.1 Route description and perturbation .......................................................................................... 19

3.2 Deterministic weather forecast and parametric wave model.................................................... 22

3.3 Ensemble weather forecast...................................................................................................... 24

4 Assessment of ships responses in waves ......................................................... 26

4.1 Ship motion and transfer function............................................................................................. 26

4.2 Statistical evaluation of the ship motion in irregular waves...................................................... 30

4.3 Acceleration on the bridge........................................................................................................ 32

4.4 Slamming.................................................................................................................................. 33

4.5 Fuel consumption and load to the main engine ....................................................................... 35

4.5.1 Determination of the over-all resistance ............................................................................................ 35

4.5.2 Propeller characteristics and operation point ..................................................................................... 37

4.5.3 Feasibility of the operation point of the main engine......................................................................... 39

4.5.4 Specific fuel consumption.................................................................................................................. 40

4.6 Practical calculation of ship responses .................................................................................... 41

4.7 Avoiding irregular frequencies..................................................................................................45

4.8 Add-on for motion sickness incidence...................................................................................... 46

4.9 Add-on for parametric rolling .................................................................................................... 49

4.9.1 Introduction of criteria from the Shin approach.................................................................................49

4.9.2 Additional criteria provided by the Krueger approach....................................................................... 51

4.9.3 Combination of the approaches.......................................................................................................... 54

4.10 Considered ship responses...................................................................................................... 57

Pareto Optimum Ship Routing

5 A multi objective, stochastic approach for the route optimization ...................... 58

5.1 Optimization setup.................................................................................................................... 58

5.2 Simplex versus Genetic Algorithm ........................................................................................... 62

5.3 Potential fuel savings by route optimization ............................................................................. 64

6 Validations and extensions................................................................................ 67

6.1 Various weather conditions ......................................................................................................71

6.2 Wave spectra variation............................................................................................................. 74

6.3 Modified hull shape................................................................................................................... 78

6.4 Constraint and threshold variations.......................................................................................... 82

6.4.1 Variation of the threshold for the vertical acceleration on the bridge................................................ 82

6.4.2 Variation of the threshold for the slamming probability.................................................................... 84

6.4.3 Comparison to thresholds posed by NORDFORSK .......................................................................... 86

6.4.4 Comparison of the motion sickness incidence MSI and significant values for vertical acceleration . 87

6.4.5 Optimizations with modified constraints ........................................................................................... 88

6.5 Parametric rolling...................................................................................................................... 91

6.5.1 Initial investigations........................................................................................................................... 91

6.5.2 Assessment of the sensitivity for the parametric rolling parameter ................................................... 93

6.5.3 Increasing of the threshold value for the forward perpendicular ....................................................... 97

6.6 Robust optimization.................................................................................................................. 98

6.6.1 Deterministic, mean- and ensemble forecast......................................................................................99

6.6.2 Robustness as constraint .................................................................................................................. 104

6.6.3 Robustness as objective ................................................................................................................... 107

6.6.4 Rectification for using analyzed weather as deterministic forecast ................................................. 109

7 Summary and outlook...................................................................................... 112

7.1 Initiation of a new approach ................................................................................................... 112

7.2 Discussion on the performance.............................................................................................. 116

7.3 Conclusions ............................................................................................................................ 119

8 References...................................................................................................... 121

Appendix 1 Stochastic evaluation..................................................................... 126

Appendix 2 Irregular frequencies...................................................................... 129

Appendix 3 Root mean square values.............................................................. 132

Pareto Optimum Ship Routing

List of Figures

Fig. 1: APL China on arrival in the Port of Seattle (by courtesy of R. Ahern, Los Angeles, CA)..................1

Fig. 2: Breakdown of operational costs (MEIJERS, 1980) ...............................................................3

Fig. 3: Isochrone method (HAGIWARA et al. 1999) ............................................................................ 3

Fig. 4: Grid for BELLMANs dynamic programming,

North Atlantic route (by courtesy of ECMWF)...............................................................................5

Fig. 5: Ensemble ship routing (by courtesy of ECMWF) ...................................................................... 6

Fig. 6: CMS HANNOVER EXPRESS (by courtesy of Hapag Lloyd).................................................. 11

Fig. 7: PARETO frontier and area of feasible solutions ........................................................13

Fig. 8: Test function, perspective and contour plot................................................................ 16

Fig. 9: SIMPLEX and GA search pattern ..............................................................................17

Fig. 10: Optimization characteristics .....................................................................................17

Fig. 11: Definition of angles and directions...........................................................................19

Fig. 12: Atlantic Express / ATX (North Atlantic) (source: www.hapag-lloyd.com) ..............................20

Fig. 13: Parent route and maximum perturbations.................................................................21

Fig. 14: Shift spline and perturbed parent route.....................................................................21

Fig. 15: Ship route in different forecasts................................................................................ 25

Fig. 16: Transfer functions for the roll motion

= 90°, Vs = 23kn.......................................28

Fig. 17: Sections CMS HANNOVER EXPRESS..................................................................28

Fig. 18: Representative transfer functions and response functions

for the added resistance in waves..................................................................................... 29

Fig. 19: Encounter frequency versus wave frequency ...........................................................31

Fig. 20: Calm water resistance...............................................................................................36

Fig. 21: Response functions for the added resistance due to waves ...................................... 36

Fig. 22: Added resistance due to waves.................................................................................37

Fig. 23: Full-scale propeller characteristic............................................................................. 38

Fig. 24: Determination of the operating point........................................................................ 39

Fig. 25: Main engine characteristic........................................................................................ 40

Fig. 26: Specific fuel oil consumption, sfoc........................................................................... 40

Fig. 27: PIERSON-MOSKOWITZ spectra for H1/3 = 1m ..................................................... 41

Fig. 28: Significant amplitudes of vertical acceleration on the bridge

for different encounter angles .......................................................................................... 42

Pareto Optimum Ship Routing

vii

Fig. 29: Significant amplitudes of vertical acceleration on the bridge

for different ship speeds...................................................................................................43

Fig. 30: Slamming probability for different encounter angles............................................... 43

Fig. 31: Added resistance due to waves.................................................................................44

Fig. 32: Brake power at reduced and design speed................................................................ 44

Fig. 33: Specific fuel oil consumption at reduced and design speed ..................................... 45

Fig. 34: Specific fuel consumption at reduced and design speed...........................................45

Fig. 35: Transfer function before and after smoothing .......................................................... 46

Fig. 36: Ship responses, disturbed by irregular frequencies (left) and corrected(right)......... 46

Fig. 37: 10% MSI isolines ...................................................................................................... 48

Fig. 38: MSI and significant amplitudes of vertical acceleration on the bridge..................... 48

Fig. 39: GMT and roll period ..................................................................................................50

Fig. 40: Lever curves for calm water, trough and crest condition .........................................51

Fig. 41: Lever curves and capsize probability........................................................................52

Fig. 42: Variation of the wave peak period............................................................................54

Fig. 43: Absolute and relative motion at FP, COG, and AP.................................................. 54

Fig. 44: Ship lines at FP and AP............................................................................................55

Fig. 45: Threshold exceeding generally and for parametric rolling conditions .....................56

Fig. 46: modeFRONTIER process flow chart........................................................................59

Fig. 47: Feasible route designs and PARETO frontier...........................................................60

Fig. 48: Optimized westbound North Atlantic crossing......................................................... 61

Fig. 49: MOGA and SIMPLEX optimization results.............................................................63

Fig. 50: Variation of free variables ........................................................................................ 63

Fig. 51: PARETO frontiers for optimizations with fixed and with unbound course............ 65

Fig. 52: Courses, velocity profiles, and specific fuel consumption

for optimizations with fixed and with unbound course....................................................65

Fig. 53: Optimization results at different wave conditions....................................................72

Fig. 54: PARETO frontiers for calm, medium and rough sea................................................73

Fig. 55: Fastest routes for different wave conditions............................................................. 73

Fig. 56: Typical storm spectra for the North Atlantic,

PIERSON-MOSKOWITZ, JONSWAP and ECMWF-1D spectrum...............................74

Fig. 57: Optimization result applying a PIERSON-MOSKOWITZ

and JONSWAP spectrum................................................................................................. 75

Pareto Optimum Ship Routing

viii

Fig. 58: Representative response functions for different ship speeds

(enlarged diagram of Fig. 18)...........................................................................................76

Fig. 59: Representative transfer functions for the vertical motion on the bridge and

for the relative motion between water surface and bow

(enlarged diagrams of Fig. 18)......................................................................................... 76

Fig. 60: Time minimum routes for the PIERSON-MOSKOWITZ

and the JONSWAP spectrum........................................................................................... 77

Fig. 61: Lines plans of VERSLUIS.051, CMS HANNOVER EXPRESS

and JOURNÉE.044 (S175) ..............................................................................................79

Fig. 62: Added resistances due to waves for different hull shapes ........................................79

Fig. 63: Slamming probabilities for different hull shapes......................................................79

Fig. 64: Significant amplitudes of vertical accelerations on the bridge

for different hull shapes....................................................................................................80

Fig. 65: PARETO frontiers and ETAmin routes for HANNOVER EXPRESS,

scaled VERSLUIS.051 and JOURNÈE.044 (S175) ........................................................ 80

Fig. 66: Velocity profiles for ETAmin routes of HANNOVER EXPRESS,

scaled VERSLUIS.051 and S175..................................................................................... 81

Fig. 67: Variation of the threshold for vertical acceleration on the bridge C_zacc ..............83

Fig. 68: Time minimum routes for varied C_zacc ................................................................. 83

Fig. 69: Optimization result for deactivated and enabled C_zacc..........................................84

Fig. 70: Parameter variations influencing the slamming probability..................................... 85

Fig. 71: PARETO frontiers applying different thresholds for the vertical acceleration ........88

Fig. 72: Optimization results for modified constraints .......................................................... 89

Fig. 73: Distance over ETA for modified constraints.............................................................89

Fig. 74: Fastest routes for modified constraints..................................................................... 90

Fig. 75: PARETO frontiers for optimizations considering parametric rolling ...................... 92

Fig. 76: Routes of minimum ETA ..........................................................................................93

Fig. 77: Histogram plots for parametric rolling criteria.........................................................94

Fig. 78: Feasible and infeasible route..................................................................................... 95

Fig. 79: Parametric rolling, route comparison........................................................................96

Fig. 80: Polar plots, left side: limiting significant wave heights depending on ship speed

and encounter angle, right side: screen shot from OCTOPUS RESONANCE

(left side CLAUSS (2008), right side http://www.amarcon.com)......................................................................97

Fig. 81: Histogram plots for a threshold for hrel of 7m and 8m .............................................97

Pareto Optimum Ship Routing

Fig. 82: PARETO optimum routes for deterministic and mean forecast............................... 99

Fig. 83: Time minimum routes for mean- and deterministic forecast..................................100

Fig. 84: Significant wave heights for minimum ETA routes

in mean- and deterministic forecast ...............................................................................100

Fig. 85: Assessment of constraints in different forecasts..................................................... 102

Fig. 86: PARETO frontiers and optimum routes in different forecasts ............................... 103

Fig. 87: Fuel consumption and mean fuel consumption ......................................................105

Fig. 88: PARETO frontiers for optimizations in a deterministic

and an ensemble forecast................................................................................................106

Fig. 89: Courses and velocity profiles of routes optimized in a deterministic

and an ensemble forecast................................................................................................107

Fig. 90: PARETO frontiers of the ensemble forecast optimization..................................... 108

Fig. 91: Courses and velocity profiles of routes with different robustness..........................108

Fig. 92: Ensemble results in analyzed weather ...................................................................110

Fig. 93: PARETO frontier and time minimum route at rough seas ..................................... 112

Fig. 94: Significant wave heights of a severe winter storm in North Atlantic,

22. January, 2002 ........................................................................................................... 114

Fig. 95: Rough weather route optimization for deterministic and ensemble forecasts.......115

Fig. 96: Symmetric autospectrum and asymmetric spectrum

of a random stationary process....................................................................................... 126

Fig. 97: 2D values for mass- M33 and damping N33 coefficients, S175 ............................... 129

Fig. 98: 2D values for mass- M33 and damping N33 coefficients, section 5 .........................130

Pareto Optimum Ship Routing

List of Tables

Tab. 1: Exemplary compilation of routing service or decision support systems .....................9

Tab. 2: Free variables and constants of the route perturbation ..............................................22

Tab. 3: Maximum wave heights of considered weather scenarios.........................................23

Tab. 4: Propeller data ............................................................................................................. 37

Tab. 5, Joint capsize probability for TP = 16.5s, Vs = 23kn, hcrit = 10m,

head- and following seas.................................................................................................. 53

Tab. 6: Considered ship attributes and summary of ranges of extreme responses................. 57

Tab. 7: Parameter of the optimization setup .......................................................................... 58

Tab. 8: Route comparison to assess the potential fuel saving................................................66

Tab. 9: Applied forecasts and optimization results................................................................ 71

Tab. 10: Comparison of JONSWAP and PIERSON-MOSKOWITZ spectra........................ 75

Tab. 11: Recalculation comparing JONSWAP and PIERSON-MOSKOWITZ spectra ....... 76

Tab. 12: Hull properties of compared ships........................................................................... 78

Tab. 13: NORDFORSK, general operability criteria (JOURNÉE, 2001)........................................ 86

Tab. 14: NORDFORSK, operability criteria for various types of work (JOURNÉE, 2001)...........87

Tab. 15: Applied constraints for the investigation of parametric rolling...............................92

Tab. 16: Optimization applying all parametric rolling constraints ........................................ 94

Tab. 17: Parametric rolling parameter for compared routes ..................................................96

Tab. 18: Optimized routes evaluated in ensemble forecasts................................................ 103

Tab. 19: Constraint for robust route optimization................................................................105

Pareto Optimum Ship Routing

Nomenclature

Units which are dependent on the particular assignment of a symbol are referred to as [ dep. ].

latin symbols

and

abbreviation

name unit or value

a index [ - ]

a, i, j, m, n control variable [ - ]

a parameter of a harmonic or random process [ dep. ]

ã amplitude of a harmonic function [ dep. ]

a1/3 significant wave amplitude [ m ]

Ae/A0 propeller blade area ratio [ - ]

Am midship cross-sectional area [ m2 ]

Am midship cross-sectional area coefficient [ - ]

AP, APP aft perpendicular -

arms root mean square of the wave amplitude [ m ]

B breadth of a ship [ m ]

C peak enhancement function of a power spectrum [ - ]

CB block coefficient [ - ]

COG center of gravity -

CP pressure coefficient [ - ]

C_parameter parameter used as constraint [ dep. ]

CWP water plane area coefficient [ - ]

D propeller diameter [ m ]

dep. unit depends on the appropriate symbol or abbreviation -

df, DF deterministic forecast -

DIST sailed distance in nautical miles [ nm ]

DOE design of experiment -

DP draft at point P [ m ]

E threshold exceedings [ times/h ]

ef, EF ensemble forecast -

ETA estimated time of arrival,

here a synonym for the duration of a journey [ h ]

ETAmin-DF time minimum route based on a deterministic forecast -

ETAmin-EF time minimum route based on an ensemble forecast -

ETAmin-EF# time minimum route based on ensemble forecast no. # -

ETAmin-MF time minimum route based on a mean ensemble forecast -

f frequency [ 1/s ]

FP, FPP forward perpendicular -

g gravity constant 9.807 m/s2

g inequality constraint [ dep. ]

GA genetic algorithm -

GMT transverse metacentric height [ m ]

Pareto Optimum Ship Routing

xii

GT gross tonnage [ m3 ]

h equality constraint [ dep. ]

H1/3, h1/3 significant wave height / double wave amplitude [ m ]

Hcrit, hcrit critical wave height [ m ]

Hrms root mean square of the wave height [ m ]

J propeller advance ratio [ - ]

J term within the peak enhancement function

of a power spectrum [ - ]

k form factor [ - ]

k wave number, 2π/λ [ rad/m ]

KQS torque coefficient [ - ]

KTS thrust coefficient [ - ]

kxx radius of inertia for the roll motion [ m ]

lat latitude [ deg ]

lon longitude [ deg ]

LCB longitudinal center of buoyancy -

LCF longitudinal center of flotation -

LPP, LPP length between perpendiculars [ m ]

Lwave wave length [ m ]

LWL length of water line [ m ]

M33 potential mass coefficient of the heave motion,

non-dimensionalized with ρ Am [ - ]

mf, MF mean ensemble forecast -

mSn nth order moment of the motion power spectrum [ (m/sn/2)2, (rad/sn/2)2 ]

n nth order moment of the wave power spectrum [ (m/sn/2)2, (rad/sn/2)2 ]

0 0th order moment of velocity power spectrum [ (m/s)2, (rad/s)2 ]

0 0th order moment of acceleration power spectrum [ (m/s2)2, (rad/s2)2 ]

MCR maximum continuous rating of a main engine [ MW at rpm ]

MOGA multi objective genetic algorithm -

MSI motion sickness incidence [ % ]

n number of revolutions of a propeller [ 1/s ]

N33 potential damping coefficient of the heave motion,

non-dimensionalized with ρ Am sqrt(2g/B) [ - ]

nef number of ensemble forecasts [ - ]

nph number of times per hour [ 1/h ]

nph@ap number of events per hour suspicious for

parametric rolling at AP [ 1/h ]

nph@cog number of events per hour suspicious for

parametric rolling at COG [ 1/h ]

nph@fp number of events per hour suspicious for

parametric rolling at FP [ 1/h ]

nph4hcrit number of wave height exceedings

suspicious for parametric rolling [ 1/h ]

Obj objective function [ dep. ]

O_parameter parameter used as objective [ dep. ]

p probability density distribution [ - ]

P parameter vector [ dep. ]

Pareto Optimum Ship Routing

xiii

P probability [ -, % ]

pcr critical pressure [ kN/m2 ]

P/D pitch-diameter ratio of a propeller [ - ]

PR parametric rolling -

PS shaft power [ kW ]

r lever vector, the vertical distance from center of gravity

to a considered point [ m ]

rAW /

response function of the added resistance in waves,

normalized mean value of added resistance in regular

waves

[ kN / m2]

RA air resistance [ kN ]

RADD additional resistance [ kN ]

RAO response amplitude operator, squared transfer function [ - ]

RAPP resistance of appendages [ kN ]

RAW, RAW added resistance due to waves [ kN ]

RB additional resistance due to the bulbous bow [ kN ]

RF response function for the added resistance

due to waves [ kN/m2 ]

RF frictional resistance [ kN ]

RMS root mean square value [ dep. ]

rpm rounds per minute, propeller revolutions [ 1/min ]

RT calm water resistance [ kN ]

RTD retarded time of departure [ h ]

RV viscous resistance [ kN ]

ζζ

(

) autocorrelation function regarding the stochastic

parameter ζ [ m2 ]

sa ship response amplitude, a = 1 – 6 [ m, rad ]

sa1/3 significant amplitude of the ship response [ m, rad ]

sa /

a transfer function [ m/m, rad/m ]

scale maximum orthogonal perturbation in percent of the

parent route length [ - ]

sfoc specific fuel oil consumption [ g/kWh ]

SMCR specified maximum continuous rating of a main engine [ MW at rpm]

slprob slamming probability [ % ]

sP vector of translatory motion at point P [ m ]

sP3rel relative vertical motion between a point P

and water surface [ m ]

spef specific fuel consumption [ t/h ]

sR rotation vector of the center of gravity [ rad ]

sT translation vector of the center of gravity [ m ]

S vector of spatial shift parameter [ - ]

SIMPLEX downhill Simplex-algorithm according Nelder and Mead -

SS (

) power spectrum of the ship response [ dep. ]

SSP (

) power spectrum of the ship response at point P [ dep. ]

(

) power spectrum of swell [ m2/(rad/s) ]

ζζ

(

) autospectrum regarding the stochastic parameter ζ [ m2/(rad/s) ]

t thrust coefficient [ - ]

t time [ h ]

Pareto Optimum Ship Routing

xiv

T draft of a ship [ m ]

T time period,

e.g. sample period of recorded wave elevation [ s ]

T0 mean zero-up-crossing wave period [ s ]

T0S mean zero-up-crossing period of the ship motion [ s ]

T1 mean wave period [ s ]

Te natural period of the roll motion [ s ]

Tenc wave encounter period [ s ]

Tp peak period of the power spectrum [ s ]

Tpitch pitch period of the ship motion [ s ]

U wind speed at 19m above the sea level [ m/s ]

vcr critical velocity [ m/s ]

vmax maximum for the velocity perturbation [ kn ]

vmin minimum for the velocity perturbation [ kn ]

V vector of velocity shift parameter [ - ]

VS ship speed [ m/s ]

w wake fraction [ - ]

x, y, z variable [ - ]

xacc significant amplitude of acceleration, x-direction [ m/s2 ]

X vector of decision variables [ dep. ]

yacc significant amplitude of acceleration, y-direction [ m/s2 ]

zacc significant amplitude of acceleration, z-direction [ m/s2 ]

Z number of propeller blades [ - ]

greek symbols name unit or value

α course angle [ deg ]

Phillips constant, used within power spectra 0.0081

decay parameter, used within power spectra 0.74

incident wave angle of swell [ deg ]

phase shift [ rad ]

shape parameter, used within power spectra [ - ]

wave encounter angle [ deg ]

open water efficiency [ - ]

relative rotative efficiency [ - ]

shaft efficiency 0.99

wave length [ m ]

density of sea water 1025 kg/m3

variance of a stochastic process [ dep. ]

standard deviation of the wave amplitude [ m ]

standard deviation of the wave surface displacement [ m ]

time lag, variable of the autocorrelation function [ s ]

frequency [ rad/s ]

Pareto Optimum Ship Routing

zero-up-crossing frequency [ rad/s ]

wave encounter frequency [ rad/s ]

peak- or modal frequency [ rad/s ]

wave elevation, wave surface displacement [ m ]

wave amplitude [ m ]

rel crit

critical relative wave amplitude [ m ]

rms

root mean square of the wave surface displacement [ m ]

Pareto Optimum Ship Routing, Chapter 1

1 Introduction to weather routing of ships

During the last years ship monitoring, routing assistance, and decision support became

important topics for safe and reliable sea transport. Decision support systems for a route

optimization regarding weather and sea state are important tools to improve the reliability of

sea trade within a transport chain or to increase comfort and safety of crew and passengers.

Generally an optimum route complies with the desired time of arrival at minimum fuel

consumption and maximum safety. Permissible and reasonable loads on the ship, the cargo,

and the crew are not allowed to exceed. This establishes a multi-objective, non-linear, and

constrained optimization problem in which a suitable compromise is to be found between

opposing targets.

Following the casual statistics of IMO (2000) for ships of 100GT and above, there were 533

registered serious and very serious casualties at sea, 64 thereof in heavy weather. These

comprise 145 total losses, 38 thereof in heavy weather. Even 3 years later in 2003, the number

of severe casualties amount to 150, all in all 109 total losses of ships greater than 100GT were

registered, thereof 40 in heavy weather. The most affected ships are general cargo ships,

accounting for nearly 20% of the world merchant fleet, but suffering over 40% of the total

losses. However, total loss is only one aspect in the consequences of heavy weather, structural

damages and the loss of cargo are some other. E.g. Container ships: in comparison to general

cargo ships they seem to be much more safe regarding the danger of total loss. In a way

container ships possess of an unmeant safety device, their stability increases in consequence

of loosing on-deck container when roll angles or transverse accelerations exceed a critical

limit. In 1998, the M/V APL CHINA gained notoriety and became the most famous victim of

the Pacific typhoon BABS. Sole the loss of cargo amounted for more than 100M US$. 406

containers were lost and more than 1000 damaged. The ship was able to reach the harbor on

its own.

Fig. 1: APL China on arrival in the Port of Seattle (by courtesy of R. Ahern, Los Angeles, CA)

Besides loss and damage of cargo and ship, each casualty bears the risk and sometimes results

in damage to persons and loss of life. Therefore each accident poses the question if it could

have been avoided. It is a concern of ship design, maintenance, operation, and prudent

seamanship. Notwithstanding the commercial pressures imposed by shipping schedules, it is

the master’s sole discretion to take whatever action he/she sees appropriate to maintain the

safety of the crew, vessel, and cargo. Driving the vessel too hard may result in damage to the

cargo or to the vessel itself. Being overcautious annoys the charterer. In this context, situation

awareness regarding oncoming weather condition and the expected impact on the ship is of

greatest importance for a route planning. It is a ship’s master who finally has to face this

challenge and who solely is responsible for route decisions and for the safety of crew, ship,

and cargo.

Pareto Optimum Ship Routing, Chapter 1

As rough weather condition are often not avoidable, it is of greatest importance for a master

to have reliable weather data available on board to reduce the risk of accidents excited by the

sea condition. Additionally, as the average temperature currently rises, extreme weather

conditions may become more probable.

To support the decision making when planning a journey various services and systems are

available. They reach from guidance by meteorologists to on-board decision support systems.

Nevertheless, an improvement of currently available decision support systems is desirable,

e.g. to increase the reliability of weather forecasts and the prediction of consequences or to

enable optimization procedures to handle more than one objective target simultaneously. The

latter one would support a master to easily decide, if he/she wants to put the focus of a route

optimization on an arrival on schedule, on fuel savings, on the avoidance of weather

conditions that may become hazardous, or on any combination of these criteria.

The aim of this work, therefore, is to suggest a novel approach to the optimization of ship

routes in adverse weather and to investigate its merits, without claiming completeness of the

model. It will be shown, however, that an extension of the optimization scheme is rather

straightforward, allowing the quick incorporation of further analyses. Finally, posing and

solving a route optimization problem as a multi objective and constraint optimization task

supports conscious decision making, also in complex situations, where various objective

targets have to be balanced.

The following overview on scientific work and routing services serve to give an impression of

the current state of the art. It is not intended to mention all available publications and services

but to illustrate the diversity of routing services and current development in research.

1.1 State of the art in applied research

Design and operation of sea trade frequently uses modern optimization and simulation

methods. They are employed in tasks like fleet planning, fleet sizing and scheduling, the

integration of the sea transport into a supply chain, acceleration of loading- and unloading

operations, and improvements of the production line. Further they support to solve particular

design tasks, e.g. for the optimization of a hull shape or the reduction of steel weight. Here

optimization applied within decision support systems for weather routing purposes is

considered. The term routing itself is used with different meanings. CHRISTIANSEN et al.

(2004) give a comprehensive overview on optimization in sea trade. They define routing as

the assignment of sequences of ports to be visited by the ships. Ship weather routing as

addressed here is called environmental routing.

With regard to the considered period of time environmental or weather routing can be:

• short term routing, regarding the current sea condition and its consequences for the

next minutes up to hours, to give advice of suited countermeasures,

• medium term routing, regarding the sea condition predicted by forecasts, to provide

recommendation of favorable routes for a particular voyage,

• long term routing, regarding long-term weather statistics, e.g. applied for fleet

planning.

In the following simply routing will be addressed, this implies weather- or environmental

routing applying medium range weather forecasts, i.e. medium term routing is considered

here.

Pareto Optimum Ship Routing, Chapter 1

The awareness of current and oncoming weather and sea conditions, at all times, are of

greatest importance for the route planning and navigation of a ship’s master. Formerly this

knowledge was derived solely by direct observations on-board a ship and by the individual

knowledge of a master about characteristic weather conditions in a particular sea area. Later

on, this was supported by means of statistical evaluation of long-term observations. Ever

since the invention of radio technology, it is possible to supply a ship during the whole

journey with weather forecasts and supplementary weather observations from other ships.

Since the sixties of the last century, ship officers are able to use routing advises from weather

routing departments of meteorological institutes, where the prediction of the ship reaction in

wind and waves is based on the officers experience with the considered ship or similar ships.

JOURNÉE and MEIJERS (1980) discuss the influence of natural speed loss, defined by rough

sea resistance and maximum load of the main engine, and voluntary speed loss, to prevent

severe motions, shipping green water, slamming, or propeller racing. In the second part of

their publication, the influence of routing on the fractions of operational costs is shown.

The biggest influence can be expected on the fuel and lubrication oil consumption. The

breakdown of operation costs by MEIJERS is illustrated in Fig. 2. A program to assess the

speed loss in a seaway is applied to the evaluation of different routes. Therein the ship

performance is assessed by means of a seakeeping method based on strip theory. For a given

weather condition advantages of a route decision are discussed. Costs for damage of cargo

and ship are neglected, as they are hardly known.

Fig. 2: Breakdown of operational costs

(MEIJERS, 1980)

Fig. 3: Isochrone method

(HAGIWARA et al. 1999)

HAGIWARA et al. (1999) present route optimizations on a 50kn high-speed catamaran for

the Sea of Japan. They use an optimization algorithm, the isochrone method, to find the time

minimum route for sea conditions described by currents and wave forecasts. This method,

illustrated by Fig. 3, repeatedly computes an isochrone, i.e. the attainable time front that

describes the outer boundary reachable from the departure point after a certain time. The

assessment of the ship performance at sea is achieved by sea trials with a smaller ship.

Minimum fuel routes for a given time of arrival are accomplished by adjusting the speed of

the time minimum route.

Pareto Optimum Ship Routing, Chapter 1

In the nineties, a group of Dutch researchers and companies created an expert system, the

Ship Performance Optimization System SPOS (SPAANS and STOTER, 2000). Already in

2000 they could look back to their experience with 75 installations of this system. SPOS

consciously is an on-board routing system. It should combine the forecast of meteorologists

about the oncoming weather with the experience of the master about the ship, who finally is

responsible for navigational planning. Initially the system used the isopone method by

HAGIWARA and SPAANS (1987), an extension of the isochrone method. An isopone is the

plane of equal fuel consumption that defines the outer boundary of the attainable region in a

three-dimensional space, i.e. position and time. But finally the isochrone method is applied,

because this method is less complex and easier to comprehend. The main reason for this

decision is the higher acceptance of the isochrone method by the navigational staff. SPOS

supplies a graphical user interface to display geographic- and weather data. Up to 5 routes can

be generated and balanced against each other, the great circle route, rhumb line, time

minimum route, and two user defined routes. The ship performance is supplied by a numerical

model, which has to be calibrated according to a particular ship by the crew. Comparisons

between ship operation with and without SPOS show that the time of operating in bad

weather can be reduced to a quarter.

To control the stochastic behavior of weather, i.e. the uncertainty of the weather forecast, in

route optimization TREBY (2002) suggests the application of ensemble forecasts instead of

their deterministic counterparts. An ensemble is a set of forecasts with an equal probability

for the point in time where the forecast starts. At later time steps the probability is adjusted

according to the concordance of observed weather and the prediction of each member of the

ensemble forecast. TREBY applies a dynamic programming algorithm, originally developed

by BELLMAN1, to detect the fastest course of a yacht in a racing competition. (Further

information on the ensemble forecast system is given in section 3.3)

The EU 5th framework research project SEAROUTES (2003) bundles the whole knowledge

on ship routing. To describe the project at best is to look for the full title:

“SEAROUTES – Advanced Decision Support for Shiprouting

based on Full-scale Shipspecific Responses as well as

Improved Sea and Weather Forecasts including Synoptic,

High Precision and Realtime Satellite Data”.

One of the most important topics in routing is to obtain reliable forecast data. The best code is

useless if the upcoming weather condition is not known sufficiently. Therefore the project

addresses possible and necessary improvements of forecasts extensively. SEAROUTES

establishes a decision support system SEAROUTES DSS for the optimization of ship routes

based on medium range weather forecasts. The assessment of the seakeeping behavior of a

particular ship is performed by strip theory- or panel methods. Comparisons of both methods

and different codes also with model trials are conducted. State-of-the-art technology for

analysis and forecasting of sea states is applied to get most accurate and reliable wave and

weather forecasts to simulate the ship performance on a particular route. For the route

optimization BELLMANs dynamic programming is used. It is a graph algorithm starting from

an initial guess, e.g. the great circle route. From this initial route an orthogonal grid is build.

Fig. 4 serves to illustrate this method.

1 Richard Ernest Bellman, applied mathematician and theoretical physicist, 1920-1984, USA

Pareto Optimum Ship Routing, Chapter 1

Fig. 4: Grid for BELLMANs dynamic programming,

North Atlantic route (by courtesy of ECMWF)

The Figure represents the first guess route and the adjunctive search grid. This grid is

evaluated step by step from the origin towards the destination of the route. For example,

starting at the most eastern point of Fig. 4, all routes to the next column of grid points are

evaluated. For this first step a cost function is calculated. According to this cost function, a

preset number of best routes to this column are stored and taken as initial courses to the next

column of grid points. Again the cost function for all the combinations of courses is

calculated, values are integrated for the whole routes from the origin to the considered column

of the search grid. Again a number of best routes so far are stored and the step to the next

column is executed. This procedure continues until the destination point is reached. If finally

a better route than the initial guess is found, the same procedure starts with this improved

route as first guess. If no further improvement is possible, the first iteration loop is done and a

second iteration on a refined grid starts. The initial mesh size is 1-2°, it will be refined during

the optimization. The most important advantage of the BELLMAN method, compared to

gradient based methods, is its ability to transcend relative optimum solutions. The method is

known to be very fast, but it is single objective and the cost function is often not conceivable.

A more comprehensive description of the method, applied to routing, is given by

HOFFSCHILDT et al. (1999).

“The benefits of any ship routing system, however good the model is, will always depend on

the quality of the forecasted weather parameters that are used to force the system”,

SAETRA (2004). Therefore, a measure for the certainty of a forecast is desirable. SAETRA

applies ensemble forecasts, produced by the ensemble prediction system EPS of the European

Centre for Medium-Range Weather Forecasts ECMWF, to give an estimate of probable

forecast errors. He shows that the spread of the significant wave height in the ensemble

forecast can be related to the spread of routes optimized within its forecast members. Fig. 5

represents a typical optimization result. The solid blue line represents the optimum route

within the high-resolution deterministic forecast. It is overlaid with error bars. The boxes

represent the 10 and 90 percentiles for the spread of ship routes, i.e. 80% of all routes can be

found inside these boxes. The whiskers on the error bars show the position of the outliers. In

that way probable deviations from the initially identified optimum route are visualized. A big

error bar represents a high probability of weather changes and vice versa. A high probability

of weather changes implies a high probability of necessary route adaptations during the

journey and low reliability for initially predicted fuel consumption and time of arrival.

Pareto Optimum Ship Routing, Chapter 1

Fig. 5: Ensemble ship routing (by courtesy of ECMWF)

Parallel to the developments in Europe, the Japanese National Maritime Research Institute

NMRI developed a navigation system called WAN, Weather Adaptive Navigation,

TSUJIMOTO and TANIZAWA (2006). The system uses the augmented LAGRANGE2

multiplier method in order to handle objectives and constraints during the optimization. In

that way optimum routes possess minimum fuel consumption for an arrival on schedule and

e.g. restricted vertical acceleration at the forward perpendicular. TSUJIMOTO and

TANIZAWA show that dependent on the objective of the optimization, i.e. whether the

objective is the minimization of traveling time or the minimization of fuel consumption, a

ship encounters harsh or calmer sea conditions. The common preconception that ships with

routing systems always encounter higher waves cannot hold. For a Panmax container vessel

on a transpacific route they fond fuel savings up to 26% compared to the great circle route at

maximum attainable speed.

Beside these developments on ship route optimization, other, related projects address on

monitoring and short-term prediction for hazardous situations.

PAPANIKOLAOU et al. (2000) put their focus in a direction that may be described as an

initial en-route monitoring or short term decision support. They establish a Seakeeping

Information Booklet SIB, to be understood as an upgraded stability booklet that describes the

“safe operational envelope” due to speed, heading, and wave condition for different types of

bulk carriers. Motivated by the loss of MV DERBYSHIRE in 1980 the considered ship

responses include bending moments. Apparently, hull strength becomes a parameter for

navigational decisions.

As a supplement to the criteria for intact stability of ships specified in IMO (1988), resolution

A.749(18), the Maritime Safety Committee approved a “guidance to the master to avoid

dangerous situations in following and quartering seas”, IMO (1995). This became necessary

as basic stability criteria turned out to be insufficient for a safe ship operation in adverse sea

conditions.

2 Joseph Louis Lagrange, comte de l’Empire , 1736-1813, mathematician and astronomer, Italy

Pareto Optimum Ship Routing, Chapter 1

The guidelines give recommendations to prevent capsize and heavy roll motion due to:

• surf-riding and broaching-to,

• reduction of intact stability caused by riding on the wave crest at midship,

• synchronous and parametric roll motion,

• and combinations of further various dangerous phenomena.

Meanwhile these guidelines are extended and supersede by the “revised guidance to the

master for avoiding dangerous situations in adverse weather and sea conditions”, IMO (2007).

As before, it is designed to accommodate all types of merchant ships. Therefore the Maritime

Safety Committee explicitly recommends tailor-made software, which takes into account

characteristics of an individual ship.

The guidelines are restricted to preventive measures for extreme roll motion and capsize,

hazards and risks in adverse weather like e.g. slamming, longitudinal and torsional stresses, or

collision and stranding are not considered. In this regard the decision-making that takes into

account all possible risks becomes a quite complex task. For this reason the demand for a

supporting tool that identifies suited countermeasures in any hazardous situation further

increases.

In October 2002, the classification society GERMANISCHER LLOYD (2002) installed a

Shipboard Routing Assistance System SRAS on-board of a post-Panmax container vessel of

the Greek shipping company COSTAMARE. The OCTOPUS software of AMARCON B.V.

builds the framework of SRAS. It serves as control unit for a wave- and a hull response

monitoring system, contains tools for route planning and evaluation regarding to weather data

and to resulting ship responses, and provides the graphical user interface for SRAS. The

system assesses the surrounding wave field and the probability of particular ship responses in

half-hourly intervals. It reports to the master if the probability of these responses exceeds

predefined threshold values. In the initial setup warnings are given for the slamming

probability and for accelerations at the bow, meanwhile the assessment of parametric rolling

is included, RATHJE and BEIERSDORF (2004).

During the last years Technische Universität Berlin contributed in several scientific projects

addressing ship stability in irregular seas, CLAUSS et al. (2003) and CRAMER et al. (2004).

Even if the focus is put on the improvement of ship design, an important result from the

operational point of view are tailor-made data bases for particular ships, e.g. represented by

polar plots that enable a master to identify dangerous wave conditions and give advice, how to

avoid them at the best. Current projects, like LaSSe – Lasten auf Schiffe im Seegang, directly

address this topic. They serve to establish an on-board system that provides a short-term

prediction of the pressure field at the outer hull regarding the approaching waves measured by

a wave monitoring system. In this way, it is possible to predict the load on the ship and ship

motions a few minutes before they occur; hazardous events that could not be predicted by any

forecast will be identified, instantaneous countermeasures can be executed, q.v. CLAUSS

(2008).

Pareto Optimum Ship Routing, Chapter 1

The EU 6th framework research project ADOPT (2005) currently develops a risk-based

decision support system for ship operation in rough seas. The ADOPT DSS – Advanced

Decision Support System For Ship Design, Operation, And Training accounts for

environmental, hydrodynamic, structural and operational factors that influence the decision

making on-board a ship. Based on measurements and acquisition of all relevant data the

system displays the current and also short, and medium term predictions of operating

conditions. Further on the system conducts an evaluation by means of simulations and risk-

based safety assessment techniques, accounts for the uncertainty of related information and

for consequences of navigational measures. By this means guidance for a conscious decision

making in terms of safe and economic navigation is given.

From different points of view it is useful to combine a forecast based route decision support

system with a short-range wave monitoring system, if this is already installed on-board. As

regards content, both systems can use similar methods and databases for the assessment of the

ship motion in a seaway. Integrated bridge design demands for the reduction of devices. Last

but not least, it is favorable for the route decision support to access data of the wave

monitoring system. In that way the reliability of a forecast can be checked on-site.

1.2 Routing services

A characteristic difference of services for routing is whether the system or support comes

from ashore or if there is an installation on-board a ship. In the first case it can be simply a

forecast for weather and waves, or the service comprises advises from meteorologists and

experienced masters based on their personal knowledge or acquired by automated

optimization methods. Sometimes the service includes vessel tracing to enable the shipping

company to watch the current positions and operating conditions of their ships.

The latter case consists of any on-board installed routing system, in most cases, a software

running on a PC to display defined routes and resulting ship responses according to a

weather- and wave forecast. It often comprises optimization functionality for the automated

detection of favorable routes. In some cases the on-board system is integrated into a

monitoring environment, i.e. an integrated system of a processing framework and

measurement devices, e.g. to measure the surrounding wave field or sensors for loads and

accelerations. A fully integrated system therefore includes medium term and short-term

response prediction, i.e. response prediction based on forecasts and measured wave fields. In

that way it provides support for navigation and maneuvering based on predicted and measured

ship responses. Generally these systems also perform voyage data recording. Tab. 1 gives an

overview on some available routing systems and services. It serves to illustrate their diversity

and is an arbitrary selection. Therefore it does not aim to give a complete overview on routing

system or service providers.

Pareto Optimum Ship Routing, Chapter 1

Tab. 1: Exemplary compilation of routing service or decision support systems

service provider service / system

weather forecast

route planning

route optimization

warning system

ship monitoring

data recording

vessel tracking

Aerospace and Marine

International

(USA)

Weather 3000, internet service, maps

displaying fleet and weather information x x x

Applied Weather Technology

(USA)

BonVoyage System, on-board system x x

Deutscher Wetterdienst

(Germany)

MetMaster, MetFerry, ashore routing

systems, advice on demand x x

Euronav

(UK)

seaPro, on-board system, software or

fully integrated bridge system x x x

Finish Meteorological Institute

(Finland)

weather and routing advice from ashore

for the Baltic sea x x

Fleetweather

(USA)

meteorological consultancy from ashore x x x

Force Technology

(Denmark)

SeaSense, real-time on-board decision

support system managing wave-induced

structural loads and ship motions x x x

Germanischer Lloyd,

Amarcon B.V.

(Germany, Netherlands)

Shipboard Routing Assistance System

SRAS x x x x x

Meteo Consult

(Netherlands)

SPOS, on-board routing system x x x x

Metworks Ltd.

(UK)

meteorological consultancy from ashore x x

Norwegian met office,

C-Map

(Norway, Italy)

C-Star, on-board system x x

Oceanweather INC.,

Ocean Systems INC.

(USA)

Vessel Optimization and Safety System,

VOSS, on-board system x x x x

Swedish met and hydrology

institute (Sweden)

Seaware Routing TM, Seaware Routing

Plus TM and Seaware EnRoute Live TM,

on-board systems, and support ashore x x x x x

US Navy

(USA)

STARS, on-board system x x x x

Weather News International,

Oceanwaves (USA, Japan)

Voyage planning system VPS and

ORION, combination of ashore and on-

board routing and optimization software x x x

Weather Routing Inc.

(USA)

routing advice from ashore and Dolphin

navigation program combined with a

web-based interactive site x x

Transas

(UK)

Ship Guard SSAS, on-board system,

software or integrated to bridge system x x x x x

If at all automated route optimization methods are applied, all these systems and services

operate single objective. That means, the optimization follows the ambition to improve a

route in one matter, e.g. to find the route with the earliest estimated time of arrival ETA or to

find the route with the lowest fuel consumption for a given ETA. However, in the latter

example the objective function already becomes unhandy. To handle arrival on schedule and

Pareto Optimum Ship Routing, Chapter 1

lowest fuel consumption simultaneously, a delay has to be penalized and weighted against the

fuel consumption. In that way, it is a combination of at least two objectives in one objective

function. This is necessary for the applied optimization algorithms but unfavorable for a

conscious decision support of a master. For this purpose a separate handling and assessment

of the applied objectives would be better, i.e. compared to a single objective method, a multi-

objective approach would be favorable.

1.3 A new approach to route optimization

So far, it appears that a route optimization procedure should be able to operate multi

objective. Furthermore it has to consider a large amount of parameters that have an influence

on a ship route decision:

• There are environmental parameters like restrictions in the navigable water, swell,

wind waves, wind, currents, and drift ice. Most of these parameters are brought into

the decision making by forecasts. Obviously an assessment of the reliability of the

forecasts is desirable.

• Ship-sided parameters like ship size, hull shape, and the load condition, have a direct

influence the ship characteristics in waves. The available main engine power forces a

reduction of speed in higher waves, structural conditions may be regarded, too.

• Operational parameters cover the whole range of situations that might cause harm to

the ship, the cargo, to passengers, or crew. For a reliable sea transport, situations that

bear the risk of capsizing by parametric rolling, surf riding, or stability loss in

following waves have to be prevented. Just as well, operational conditions that enable

overcritical slamming, propeller racing, large amounts of green water on deck, and

extreme accelerations should be avoided. Normally on-site countermeasures like

course changes or speed reduction are employed to avoid damage of the ship or cargo

or to guarantee the well-being of passengers and crew.

Dependent on a particular service the main aspect of a route decision is variable. Cruise liners

may put focus on the well-being of passengers, a container service looks for safety of cargo

and for being on schedule. Both are interested in reducing operating costs, mainly fuel- and

lube oil consumption. Furthermore the global load on the ship structure becomes more and

more deciding for a route decision, too.

This study establishes a new approach to route optimization. It is multi-objective and

therefore optimization results are easy to overview. The method is able to produce meaningful

results, e.g. to identify fuel minimum routes even if an arrival on schedule is impossible due

to severe weather conditions. For the set-up two objectives are considered. They are fuel

consumption representing operational costs and estimated time of arrival ETA. Later on the

approach will be extended, considering the sensitivity of optimized routes to probable weather

changes. Regarding the environmental conditions, it is decided to put the focus on swell as

being the crucial factor for the ship motion in waves. For the time being, wind, wind waves,

and currents are left out of consideration. The load on the main engine is observed to

guarantee that the operational conditions agree to the main engine characteristic. Initially the

operational parameters, as posed in the SEAROUTES project, are adopted, i.e. thresholds for

slamming probability and accelerations on the bridge. These parameters, defining the “safe

operational envelope”, are extended at a later point of this study. The presented optimizations

belong to a westbound Atlantic crossing of the HAPAG LLOYD Panmax container vessel

CMS HANNOVER EXPRESS. This ship is used for the set-up of the route optimization

process since she is well known from the SEAROUTES project.

Pareto Optimum Ship Routing, Chapter 1

Fig. 6: CMS HANNOVER EXPRESS (by courtesy of Hapag Lloyd)

It is the purpose of this work to document the set-up and to make a feasibility study for this

new approach to ship route optimization, not to develop turnkey ready software. As a start,

chapter 2 introduces some basic principles and terms of optimization that are important for the

understanding of the applied optimization procedures. Chapter 3 illustrates navigational

aspects, i.e. route planning: determination of course and ship speed and its numerical

representation. Furthermore the applied forecasts and the numerical modeling of swell is

shown. Following, the fundamentals and the set-up of the mathematical model of the ship in

waves are presented in chapter 4. The modeling of two extensions, namely for motion

sickness incidence and parametric rolling, and some numerical recipes are given here, too.

Chapter 5 continues with the set-up of the optimization procedure and shows the

implementation of the numerical ship into a generic optimization environment. In chapter 6,

the capability of the route optimization approach is assessed by means of different sensitivity

studies. Finally a summary of this work and some aspects for further research are given in

chapter 7.

Pareto Optimum Ship Routing, Chapter 2

2 Basic principles and terms of optimization

Before transforming the routing problem into an optimization task, it is necessary to introduce

some basic terms and methods of optimization. The principle of optimization itself is very

old. It is a matter of opinion, if optimization exists since the beginning of life on earth or since

the beginning of humankind. However, since humans are aware about their capability to

decide, optimization is a conscious behavior.

In a mathematical definition, the term optimization refers to the study of problems in which a

real function is to be maximized or minimized by systematically choosing the values of

variables from within an allowed set and by means of acceptable computational resources.

The first optimization technique, which is known as steepest descent, goes back to GAUSS3.

Since that time, reams of optimization algorithms have been developed, generic ones as well

as ones that are tailor-made for a particular problem.

The problem statement in modern engineering technology bases on optimization tasks

initially posed within the field of operations research. These became solvable due to the

progress in computer technology. In this respect the modern age of numerical optimization

began in the middle of the last century. Since that time the progress in optimization is closely

related to the progress in computational capability and additional applied methods, e.g.

simulation tools for computational fluid dynamics. Parallel to the developments in other fields

of engineering, optimization in maritime design and operation stepwise improved to solve

assignments of increasing complexity. Based on the research studies at the University of

Michigan, see e.g. BENFORD (1965), NOWACKI et al. (1970) optimized main dimensions

of a tanker to achieve minimum required freight rates. Numerous projects of applied research

at Technische Universität Berlin represent the progress in optimization with regard to ship

and offshore design since the seventies of the last century; e.g. NOWACKI and LESSENICH

(1976) optimized main dimensions of a tanker, bulk carrier and a general cargo ship with

regard to the required freight rates, HARRIES (1998) established a comprehensive method for

parametric design and hull form optimization regarding ship resistance coupling computer

aided design, methods for computational fluid dynamics, and optimization, BIRK and

CLAUSS (2001) applied optimization techniques to a parametric hull form variation of

offshore structures to minimize downtime.

For a continuative introduction to optimization and a comprehensive overview on

optimization methods the study of BIRK and HARRIES (2003) is recommended, in particular

because it is devoted to the optimization in marine design.

The least common denominator of al assignments in optimization is to minimize or maximize

an objective function. Therefore the minimum requirements to state a problem as an

optimization task, is to transform it into one objective function:

()

PXfO ,=, (2.1)

i.e. a closed mathematical description of the problem in that way that the minimum

(maximum) of this function represents the desired optimum. Then, optimization is the process

of finding the set of decision variables

for a set of given parameters

that represent the

minimum (maximum) of the objective function. Objective functions can be linear, uni- or

multi-modal. In the first case the optimum (minimum or maximum) is mostly achieved at a

boundary of the decision variables. Unimodal functions contain one single optimum, whereas

multi-modal functions can afford suboptimal minima or maxima or more than one absolute

3 Carl Friedrich Gauß, mathematician, astronomer, geodesist and physicist, 1777 – 1855, Germany

Pareto Optimum Ship Routing, Chapter 2

extremum of the same magnitude, named local optima. If there is one optimum that exceeds

all the others, it is called global optimum. Certainly, it is the challenge of an optimization to

prove that there is a global optimum and to identify it.

If the solution of the problem is subject to one or more boundary conditions it is called

constrained. In principle two different types of constraints, equality- and inequality

constraints, exist:

()

NiPXh i,...,1,0, == , (2.2)

()

MjPXg j,...,1,0, =≥ , (2.3)

A set of decision variables that violates these constraints, i.e. one or more constraints obtain

0≠

h or 0<

g, is called infeasible. The constraints that cause infeasibility are also

referred to as active constraints. If no violation of constraints exists, these constraints are

inactive. Sometimes it is favorable to have the possibility to switch constraints on or off. In

this case they are referred to as enabled or deactivated. An enabled constraint can be active

or inactive, a deactivated constraint has no influence on the optimization.

In general, the decision variables represent a particular design of the subject matter that is to

be optimized. The term decision variables, therefore, is the equivalent to design variables.

The set of all designs that can be built by the set of permissible permutations of the design

variables is called design space. In case of active constraints it is subdivided into areas of

feasible designs and areas of infeasible designs. Furthermore constraints can generate local

and global optima although the objective function is linear or unimodal.

If more than one objective function is addressed, the optimization is called multi-objective.

In doing so, it is possible that two or more objectives represent a certain contradiction and the

design space is bordered by a set of solutions for which an improvement in one objective can

only be realized by impairing another objective. Such a border is called PARETO4 frontier

and the designs of this set are PARETO optimal designs, i.e. there is no design that is the

best with regard to all objectives.

Fig. 7: PARETO frontier and area of feasible solutions

Fig. 7 serves to illustrate an optimization aiming a minimization of two opposing objective

functions. The solid line marks the PARETO frontier.

4 Vilfredo Pareto, engineer ,economist and sociologist ,1848-1923, Italy

Pareto Optimum Ship Routing, Chapter 2

All designs that are represented by function values on the left side and below this frontier are

infeasible. Designs represented by function values above and right from the frontier are not

necessarily feasible, but feasible designs can only be found in this area. The designs that mark

the PARETO frontier, or designs that are at least close to these, represent optimum solutions.

Dependent on the weighting of the objectives, it is now possible to consciously select a

particular solution. Methods of multi-criteria decision-making can follow up but are

consciously ignored within this study.

Related to the optimization task, it is necessary to find a suited optimization algorithm.

Generally they can be divided into two types: deterministic and stochastic algorithms. To

step forward towards the optimum, deterministic methods use e.g. function values, gradients,

and higher derivations to define a new set of design variable. These algorithms are fast but

tend to stick to local optima. Stochastic methods always produce a portion of their designs by

a random process. Thus they are able to "jump" trough the solution space. This enables them

to avoid premature convergence at local optima. In return, an optimization mostly takes more

time.

In the following, two algorithms will be discussed more detailed as they are applied within

this study on route optimization. They are a SIMPLEX algorithm representing a deterministic

method and a genetic algorithm GA, a stochastic method. Both methods require a number of

initial designs that depends on the dimension of the solution space. n

Xℜ∈ yields a n-

dimensional solution space.

2.1 SIMPLEX algorithm

The SIMPLEX method by NELDER and MEAD (1965) is a commonly used, nonlinear

optimization algorithm for optimizing an unconstraint objective function in an n-dimensional

space. The algorithm works in principle as follows:

1. It requires 1+n initial designs or starting points, i.e. for a 2-dimensional problem

they build a triangle, a tetrahedron for 3-dimensional, and a polytope (simplex) with

1+nvertices for n-dimensional problems.

2. The objective function is evaluated for all 1+n points. If one of the designs is good

enough the algorithm stops.

3. Otherwise the worst point is deleted and replaced by a new one. In this way a new

polytope is built.

4. Continue the loop with step (2).

The core of the algorithm is the strategy in step (3) to build the new simplex:

a) Reflect the worst point about the centroid of all others.

b) If this point is now better than all the others, expand the step in the same direction.

c) If it is simply better than before, continue with step (2).

d) However, if the point gets worse in either case, shrink the simplex e.g. to half size

around the best point and continue at step (2).

The algorithm terminates if the attainable improvement at successive optimization loops falls

below a preset convergence limit. Meanwhile the SIMPLEX method is widely extended, e.g.

to preserve an equable shape of the simplex or even some random capabilities are adopted to

overcome the convergence at local optima. For the discussion on the convergence properties

Pareto Optimum Ship Routing, Chapter 2

see LAGARIAS et al. (1988). Although it is not formally proved for more than 2=n, the

SIMPLEX algorithm is widely used in high dimensional optimization problems.

2.2 Genetic algorithms, GA and MOGA

“The systematic utilization of stochastic events is one of the recipes of success of evolution”,

SCHÖNEBURG et al. (1994). It is a process of searching mainly based on three principles:

mutation, recombination, and selection. By inspiration of the biological archetype, two

types of evolutionary algorithms have been built independently, one is called evolutionary

strategy ES by RECHENBERG (1973) and the other one genetic algorithms GA by

HOLLAND (1975). In the beginning ES and GA were rather different, e.g. ES used a real-

valued coding and a variation of the genotype was solely done by recombination whereas GA

used binary-valued coding and a lot of effort was put to various mutation schemes. But even

more, it was a philosophical dispute; ES supporters aimed to develop a strong optimization

tool for engineering problems; GA supporters wanted to improve the understanding of natural

adaptation processes and to design artificial systems having properties similar to natural

systems. The genetic algorithms nowadays often are a merging of both strategies (as is

known, some who are deeply involved in the matter would strictly disagree). Nevertheless,

the ability to overcome local optima by applying stochastic schemes unifies both algorithms.

For the initial set-up to solve the route optimization problem, it is decided to apply a

multi-objective genetic algorithm MOGA of the generic optimization environment

modeFRONTIER by ES.TEC.O (1999). In the following, the focus is put on the genetic

algorithm. Surely any other stochastic method may work, too.

The basic idea of a genetic algorithm is as follows: the genetic pool of a given population

potentially contains the solution, or a better solution, to a given adaptive problem. This

solution is not "active" because the genetic combination on which it relies is split between

several subjects. Only the association of different genomes can lead to the solution.

HOLLAND’s method is especially effective because it does not only consider mutation

(mutations improve very seldom the algorithms), it also utilizes genetic recombination

(crossover). This recombination, the crossover of partial solutions, greatly improves the

capability of the algorithm to approach and eventually find the optimum, EMMECHE (1994).

In fact, the desired solution may happen not to be present inside a given genetic pool, even a

large one. If so, mutations allow the generation of new genetic configurations that widen the

gene pool and improve the chances to find the optimal solution. Mutation diversifies the

genetic pool and prevents premature convergence, whereas selection affects a harmonization.

To obtain a satisfying performance of the algorithm it is necessary to align the associated

mechanisms for mutation, selection, and recombination accurately. As genetic algorithms are

based on heuristics, no general convergence is proved.

The principle steps of a genetic algorithm are as follows:

1. Encode the problem in a binary string, i.e. transform the free variables that describe a

possible solution design (phenotype) into a suited binary representation (genotype).

2. Randomly generate a set of initial designs, the first generation or start population.

This one includes the genetic pool representing a group of possible designs.

3. Evaluate the

fitness value for each design, i.e. a ranking of the designs. Mostly the

value of the objective function relative to the one of the other designs is used. In case

of a multi-objective optimization it is deciding if a design is a member of the

PARETO frontier or if it is closed to it.

4. Select the designs that will mate according to their share in the population global

fitness. The probability to be taken into consideration to create descendants is

Pareto Optimum Ship Routing, Chapter 2

dependant on the fitness. The consideration of the best designs within the next

generation is only ensured if elitism is applied.

5. The designs that fail the selection are deleted. All others are used to fill up the

population again. As a general rule, the number of designs in each generation remains

the same. The reproduction to build new designs is performed by crossover and

mutation mechanisms, i.e. a permutation and random variation of the binary strings

that represent the free variable.

6. Start again from point (3) if no convergence limit or the number of maximum

generations is reached.

The potential to produce design variants directly depends on the number of designs per

generation. But as the mechanism that brings forward the process towards better solutions

(step 4 and 5) is accessed only once per generation, a large number of designs decelerates the

propagation towards a better solution. In a sense, the setup of a genetic algorithm is an

optimization process itself.

To illustrate the mode of operation and the performance of a SIMPLEX an a GA, both

methods are applied to maximize the following bivariate test function, available as the peaks

function within MATLAB (2002):

() ( )

(

)

()()

()

(

)

2253

2y-1x-exp31-y-x-expy-x-5x-10-1y-x-expx-13z +⋅+⋅= . (2.4)

Fig. 8 represents a perspective and a contour plot of equation (2.4), a function with three

humps and two hollows.

Fig. 8: Test function, perspective and contour plot

Both algorithms used the same start designs represented by square dots at the lower left edge

of the contour plots in Fig. 9. Because the SIMPLEX requires only three starting points

whereas the GA needs at least four, the first optimization step of the SIMPLEX is used as the

fourth starting point for the GA. Both methods are supposed to find the maximum value of the

test function at x = 1.6 and y = 0.0.

Pareto Optimum Ship Routing, Chapter 2

Fig. 9: SIMPLEX and GA search pattern

Obviously the SIMPLEX converges directly towards the nearest local maximum. The

maximum is crossed slightly but not far enough to detect the humps behind. The optimization

stops due to an achieved convergence limit. The GA indeed is able to detect the global

optimum. Admittedly it needs much more designs to reach it. It mostly preserves the

orthogonal characteristic given by the start designs. Fig. 10 illustrates the convergence

process of both optimization methods. The SIMPLEX, represented by red asterisks, quickly

propagates towards the nearest local optimum. It needs only 80 designs to converge.

Normally nothing is known about the function to be optimized, therefore it is only proved that

a local optimum is found. However, the GA needs some time "to come out of the corner.“ For

a GA the start population is chosen unfavorable. It has to increase the gene pool, to spread out

and become more moveable by mutating and recombining the initial designs. At design

number 415 it reaches the global optimum.

Fig. 10: Optimization characteristics

As there is no guarantee that really the global optimum is found, it is recommended to

continue the optimization to support at least the assumption that the optimization succeeded.

In this case, the price for reaching the global optimum is the tenfold number of required

designs. Nevertheless, it shows the ability of a GA to prevent a premature convergence and

not to stick to local optima. It is interesting to see that both, local maxima are apparently

Pareto Optimum Ship Routing, Chapter 2

ignored by the GA. If there has been a very slim maximum, it probably would have been

overlooked. Finally, there is no guarantee, it is a trade-off between computational effort, and

strategies to improve the global convergence, e.g. test various start designs or make variations

of the parameters that control the algorithm.

In spite of all difficulties connected to optimization, without these methods a route

optimization problem as attended here would become intractable. Therefore it is very

impressive to realize the number of possible designs that are considered. In section 5.1 it is

shown that at least 17 free variables are required for the applied routing example. Assuming

that each variable can take 20 discrete values (essentially it is much more), the number of

possible design variants would be:

designspossibleofnumberstepsdiscrete iablesfreeofnumber =

var , (2.5)

which results to 2217 103.120 ⋅= designs. If one design evaluation requires 2.5s (it is currently

the time required for the assessment of the ship behavior at 60 points of a single route) the

evaluation of all possible designs would take 1012 millennia. The present version of the route

optimization runs 20 000 designs and requires at most 14h. Of course this is too much, but is

has to be borne in mind: the calculations are conducted on a simple 2.8GHz personal

computer. There have been no efforts with regard to the programming language to optimize

the run-time, and the GA is not adapted to the route optimization problem. In most cases it

requires 3000 to 5000 designs to identify the PARETO frontier, all further designs are used to

prove the result. It is therefore accepted to be a reliable and efficient method to optimize ship

routes.

Pareto Optimum Ship Routing, Chapter 3

3 Modeling of ship route and environmental conditions

This chapter introduces the applied methods to generate a ship route, section 3.1, and to

describe the corresponding weather conditions, sections 3.2 and 3.3. The used definitions of

angles and directions are illustrated in Fig. 11:

• The course angle

describes the forward direction of the ship. It is 0° for a ship

sailing in northern direction, increasing towards East, i.e. 90° means sailing towards

East.

• The incident wave angle

describes the direction of swell. It is 0° for wave crests

propagating from North to South, increasing towards East, i.e. 90° means that wave

crests propagate from East to West.

• The wave encounter angle

describes the direction of approaching waves relative to

the sailing direction of the ship. It is 0° for following waves and 180° for heading

waves. Encounter angle of 90° represent beam seas from either starboard or portside.

• Locations on-board the ship are given in a ship-bound coordinate system with positive

x-direction ahead, y-direction to portside, and z-direction upward. The origin of this

coordinate system is at the center of gravity of the ship.

• All routes and maps are represented on a Cartesian coordinate grid, ordinate and

abscissa depict degrees of latitude and longitude, respectively.

Fig. 11: Definition of angles and directions

3.1 Route description and perturbation

The schedule of the example container service, used for the set-up of the route optimization,

is given in Fig. 12. It depicts the stations of the round trip and the durations of each step of the

journey, e.g. the trip from Southampton to New York takes 8 days.

Pareto Optimum Ship Routing, Chapter 3

Fig. 12: Atlantic Express / ATX (North Atlantic)

(source: www.hapag-lloyd.com)

The route optimization focuses on the open water part of the journey. Pilot time and estuary

traveling are deliberately left out of the investigation as they admit only small variations. For

that reason the route optimization starts at the western end of the English Channel

(5°W, 49°N) and ends in the estuary of New York (70°W, 41°N). It spans a distance of

roughly 3000nm. Sailed at 23kn it would take 130h, approximately 5.5 days. Compared to the

schedule of Fig. 12, the journey from Southampton to New York comprises 2.5 days pilot

time and 5.5 days at open water. Only for the open water part, the ship route is varied in time

and space. Thus, the route optimizations presented here necessitate a time period of less than

130h for the optimized part of the journey to arrive on schedule.

Necessary parameters to describe the route are longitude, latitude, course angle, and speed as

time-dependent values. For a discrete route description, as applied here, these parameters have

to be given in a sufficient resolution that ensures that all necessary environmental properties

are included. Further parameters that are regarded for the decision support are derived

according to this description, e.g. the significant wave height or the fuel consumption.

To design both, the spatial and the temporal pattern of a route in terms of free variable, the

course and the velocity profile are expressed as B-splines. An initial course and an

appropriate velocity profile, referred to as parent route, are given in advance. Perturbations of

the parent route in time and space are realized by superposing GREVILLE5-spaced shift

splines to the parent splines. The solid line in Fig. 13 shows the course spline of the parent

route for a westbound North Atlantic crossing. The route follows the great circle from

Bishops Rock to the region south of Newfoundland. Afterwards it continues rhumb line to

New York. The dashed lines represent the maximum northern and southern perturbations. For

the investigations presented here, the maximum allowed perturbation is set to 10% of the arc

length of the parent route to either starboard or portside direction. The northern perturbation is

constricted to avoid land collision at Newfoundland and at Southern Ireland. Alternatively a

validation of a route in terms of geographical feasibility is performed during the optimization

process. A variable representing a later departure is introduced to model a waiting of the ship

for better sea conditions.

5 Thomas Nall Eden Greville, mathematician, 1910 – 1998, University of Wisconsin, USA

Pareto Optimum Ship Routing, Chapter 3

Fig. 13: Parent route and maximum perturbations

Fig. 14 presents shift splines and example routes as realized with seven parameters for the

spatial shift. The horizontal positions of the shift spline vertices are computed according to

the method of GREVILLE-abscissa. The vertical offsets are taken as the free variables of the

optimization task.

The variation of the velocity profile is accomplished in the same way. Here the parent spline

straightens up and represents a ship speed of 24kn. The perturbation is performed by

superposing a shift spline in order to define a reduction of the design speed.

Fig. 14: Shift spline and perturbed parent route

Applying the GREVILLE-abscissa provides a harmonization of the variability of the splines

that are used for the course- and the velocity profile perturbation. The method couples the step

size of the control variables to the arc length of the splines. E.g. all splines are parameterized

from 0 to 1. A control variable of 0.5 represents a point in the middle of the splines. This

means, the number of passed vertices and of vertices lying ahead is the same. That is why the

variability of the spline segments also is the same. This is maintained for each value of the

control variable. Thus the variability of passed and remaining spline segments for the course

and the velocity profile is the same. A detailed description of the applied spline technique is

given in FARIN (2001).

Pareto Optimum Ship Routing, Chapter 3

Tab. 2 summarizes free variables and constants as used for the route perturbation.

Tab. 2: Free variables and constants

of the route perturbation

type name content

free variables S = [s1 – sn ] vector of spatial shift parameter

V = [v1 – vn ] vector of velocity shift parameter

RTD retarded time of departure

constants scale maximum orthogonal perturbation in

percent of the parent route length

vmin minimum for the velocity perturbation

vmax level of the parent velocity spline and

maximum for the velocity perturbation

Former investigations concerning the number of free variables showed that at least nine

variables for the description of the velocity profile and seven variables to describe the course

are necessary to obtain a sufficient temporal and spatial resolution of the route that correlates

to the wave pattern as well as to the geographical conditions in the North Atlantic,

cp. HINNENTHAL and SAETRA (2005).

3.2 Deterministic weather forecast and parametric wave model

Meteorological institutes like the European Centre for Medium-Range Weather Forecasts

operate global and associated local climatologic models. These are initialized and tuned with

atmospheric and environmental data, measured by satellites, meteorological stations ashore

and at sea, or special devices for particular measurements like wave buoys allocated over the

whole globe. Dependent on required accuracy and computational effort, these models make

use of a more or less dense numerical grid for which they calculate a wide range of

climatologic data. These data are used to produce e.g. weather- or wave forecasts.

In principle, the ECMWF provides three different types of deterministic sea state forecasts.

Therein wave data for a particular position and time are available as:

• Directional spectra, representing the energy distribution of swell and wind sea

depending on frequency and directional distribution,

• Integrated directional spectra, representing the energy distribution depending on the

frequency for a prevailing direction (referred to in the following as ECMWF-1D), or

• Wave parameter derived from ECMWF-1D; significant wave height H1/3, peak period

Tp and prevailing direction

. For practical purposes, e.g. when evaluation ship

responses, these parameters are used in combination with standard spectra like the

JONSWAP spectrum and the modified PIERSON-MOSKOWITZ spectrum,

respectively.

For the work presented here, three different weather scenario illustrated in Tab. 3 are applied.

They describe a calm weather condition, a rough sea state, and an intermediate condition. The

data cover the North Atlantic from 0.0°- 82.5°W and 18.0°- 70.5°N. They are given on a grid

of 1.5° an in time steps of 12h. Actually, these data are re-analysis data from the data archive

of ECMWF. Anyhow, they are used as deterministic forecasts for the set-up and within most

investigations of this study.

Pareto Optimum Ship Routing, Chapter 3

Tab. 3: Maximum wave heights of considered weather scenarios

sea condition date maximum H1/3 [m] corresponding Tp [s]

calm 6.6. – 16.6.2001 5.4 10.2

medium 1.1. – 11.1.2003 12.0 15.4

rough 20.1. – 30.1.2002 15.0 14.9

For practical purposes only the parametric description, i.e. significant wave height H1/3, peak

period Tp and wave angle

, is used. This is necessary to keep the computation time and the

amount of data in a reasonable magnitude. Anyhow, for the computation of the ship motion in

waves as applied here, a representation of the sea condition by means of a wave energy

spectrum is needed. Various parametric representations are available. They all refer to a form

given by PIERSON and MOSKOWITZ (1964), see e.g. CLAUSS et al. (1988 / 1994) and

HASSELMANN et al. (1975):

()































−= −

52 exp

βωαω

gS for 0>

. (3.1)

With:

0081.0=

, and 74.0=

Ug=

, and =U wind speed 19.5m above the sea level,

the original form of the PIERSON-MOSKOWITZ spectrum is obtained.

The parameter

and

are fixed for closest fit to measured data of developed wind sea for the

North Atlantic. The transformation of (3.1) to the commonly used two-parameter form is

shown by CLAUSS et al. (1994).

()











⋅−⋅= −− 4

316exp4 3

ωπω

S (3.2)

represents the PIERSON-MOSKOWITZ spectrum for the parameters significant wave height

H1/3 and zero up crossing period T0. Following the relation:

99771.0

1TgH ⋅⋅⋅=

, (3.3)

equation (3.2) is identical to (3.1). Meanwhile, this form is applied for independent values of

H1/3 and T0 resulting in a modification of the spectra. For that reason, it is called modified

PIERSON-MOSKOWITZ spectrum.

()











⋅−⋅= −− 4

328.22exp57.5 3

ωπω

S, (3.4)

represents the PIERSON-MOSKOWITZ spectrum for the parameters significant wave height

H1/3 and mean period T1, also known as BRETSCHNEIDER spectrum, BRETSCHNEIDER

(1957), or ITTC two-parameter spectrum. To apply equation (3.2) or (3.4) to the wave data

given by ECMWF, the relation p

TT ⋅= 711.0

0 or p

TT ⋅= 772.0

1 can be used, JOURNÉE

(2000).

Within the Joint North Sea Wave Project in 1968 - 1969, extensive wave measurements in the

North Sea were carried out. It was found out that multiplying (3.4) with a peak enhancement

function is most suitable to describe the measured spectra. Equation (3.5) represents the

JONSWAP formulation given by LLOYD (1998):

Pareto Optimum Ship Routing, Chapter 3

()

S⋅⋅











⋅−⋅= −− 658.028.22exp57.5 4

ωπω

, (3.5)

where C⋅658.0 is the peak enhancement function,

with J

C3.3= and 





























−

−=

0.1

exp

ωω

pp T

2= is the frequency of the spectral peak or modal frequency,

and a parameter

with 07.0=

for p

<, and 09.0=

for p

>, to scale the peak width

at frequencies above and below the modal frequency. According to JOURNÉE (2000) the

relation p

TT ⋅= 834.0

1 can be used to provide the mean period used within equation (3.5).

In 1978 the JONSWAP spectrum was advised, by the 15th ITTC, for coastal waters with

limited fetch. For a given significant wave height and wave period both spectra, JONSWAP

and PIERSON-MOSKOWITZ, have the same area under the spectral curve and therefore they

represent the same amount of wave energy. Thus, for underdeveloped wind sea both

formulations are admissible. It will be illustrated in section 6.2 (cp. Fig. 56 on page 74) that

the ECMWF-1D often lies, so to speak, “in between” the PIERSON-MOSKOWITZ and the

JONSWAP spectrum. Standard practice for the North Atlantic is the application of a

PIERSON-MOSKOWITZ spectrum.

3.3 Ensemble weather forecast

The deterministic forecast represents one likely development of the sea state but gives no

information about the probability of its occurrence. In contrast, the operational ensemble

prediction system EPS generates a set of 50 forecasts. The data format, i.e. grid and time

resolution, are the same like in deterministic forecasts. These ensemble members have an

almost equal probability of occurrence. They are generated by a superposition of small

perturbations to the operational analysis before launching the forecast calculation. The

method is based on the assumption that the forecast calculation is able to properly predict the

weather development, if it is based on a proper data set of the current weather situation. In

this view the error of the forecast is caused by errors in the measurement of the current

weather situation, and the perturbations of the start data are related to this error in the

measurements.

The ensemble spread measures the “differences” between the ensemble members, or the

differences between particular parameters given by the forecasts. It is related to the forecast

uncertainty. Small spread indicates low forecast uncertainty, and vice versa. Generally the

uncertainty increases while going further into the forecast, BIDLOT et al. (2002). Ensemble

forecasts are used to assess the amount of possible weather changes for the operation of ships

and offshore devices. Here they will be applied to evaluate the robustness of optimized ship

routes against weather changes.

Fig. 15 shows significant wave heights predicted for one particular ship route. The Figure

serves to illustrate the differences between deterministic- (re-analysis) and ensemble forecast.

Furthermore a mean ensemble forecast is shown. This one can be used instead of a

deterministic forecast to get a more likely estimate of the weather development but it is to

exercise with caution. The building of a mean value attenuates the time correlation of the data

and therefore a mean forecast undervalues the sea condition when maximum values show a

distinct time shift.

Pareto Optimum Ship Routing, Chapter 3

Fig. 15: Ship route in different forecasts

The route in Fig. 15 takes 143 hours. During this time three areas of adverse weather will be

crossed. It can be seen that the spread in the predicted wave height increases with an

increasing forecast horizon. In other words, the accuracy of the forecast decreases while it is

progressing in time. But also, and this is much more harmful for ship routing, the spread

increases in adverse weather situations. Commonly used as a rule of thumb, it is said that

deterministic forecasts are reliable within a three days horizon. This obviously is not the case

in the situation depicted in Fig. 15. Here the ensemble forecast demands on caution when

approaching the first extreme wave field approximately 30 hours after departure.

Pareto Optimum Ship Routing, Chapter 4

4 Assessment of ships responses in waves

It is expected that route optimization at harsh sea conditions is strongly affected by a

restricted operability of the ship. For this reason all factors that can cause harm to the ship, the

cargo, to passengers, or to the crew should be considered to support safe ship operation:

• regarding the absolute motion of the ship there are extreme accelerations in

longitudinal, transverse, and vertical direction and different mechanisms of extreme

roll motion that cause large roll angle or accelerations, e.g. resonant rolling,

parametric rolling, stability loss in following waves, surf riding and broaching, as well

as capsize,

• regarding the relative motion between ship and water surface there are propeller

racing, slamming, or large amounts of green water on deck,

• regarding operational parameters there are structural loads and main engine load due

to the increasing ship resistance in waves.

Here only a choice of parameters is considered just to investigate the functioning of the

suggested new routing approach, they are slamming, vertical accelerations on the bridge and

main engine load. Because an optimization requires numerous calculations of the ship

responses in an actual sea state, a fast method is needed to describe the ship in waves. In this

case the statistical evaluation method, using transfer functions provided by strip-theory and

wave spectra, is applied. This serves to estimate the ship responses for a particular route and

corresponding weather conditions.

Initially this chapter introduces the application of transfer functions and summarizes the basic

equations of this method, as far as they are applied for this route optimization approach.

Following, the conversion into a “response database” for the ship motion is shown. This step

is favorable to accelerate the calculation of ship responses during the route evaluation. The

chapter continues with a critical view on the applied methods and closes by introducing two

approaches that are not common in classic strip-theory. These address motion sickness

incidence and parametric rolling.

4.1 Ship motion and transfer function

In general, the response or output signal of a linear system is calculated from the product of

an input signal and a response function. In this case the input signal is a spectrum of an

irregular sea state, the response function is the squared transfer function of the ship motion,

and the output signal represents a ship motion spectrum.

() ()

ωζ

s⋅=

(4.1)

The input signal is the energy density spectrum of a particular sea state S

(

), in short the

wave spectrum. It can be derived by means of a FOURIER6 analysis of a measured sea state.

In this way, irregular long crested waves are represented as a superposition of harmonic

waves of frequency

, amplitude

a, and with a random phase shift. The spectrum contains

information about the amount of energy of a sea state and about the distribution of energy

over the frequency range of the spectrum. It can be used to calculate e.g. the significant wave

height H1/3, peak period TP, and other characteristic values of a sea state.

6 Jean Baptiste Joseph Fourier, mathematician and physicist, 1768 – 1830, France

Pareto Optimum Ship Routing, Chapter 4

Ss (

) is the energy density spectrum of the ship motion, in short the response spectrum. Like

the wave spectrum, it does not contain any information about the phase relation of the

harmonic components, i.e. there is no information about the motion in time. Also the response

spectrum will be evaluated statistically. For example the significant amplitudes of the motions

in the six degrees of freedom (the interest is mainly on heave-, roll- and pitch motion) and

characteristic periods (e.g. mean zero-up-crossing periods) can be determined.

The transfer function |sa/

a| represents the relation of amplitudes of the ship motion sa and the

exciting harmonic wave component

a. It is dependent on the frequency of excitation. For

each degree of freedom of the ship motion, a transfer function can be calculated. Three

translatory degrees of freedom are considered: surge, sway, and heave (index a = 1,2,3) and

three rotational degrees of freedom: roll, pitch, and yaw (index a = 4,5,6). |s3a/

a| e.g. is the

transfer function of the heave motion. In ship hydrodynamics a squared transfer function is

commonly denoted as response amplitude operator RAO.

In the strict sense a transfer function consists of magnitude and phase, i.e. for each frequency

of the input signal a magnification factor and a phase shift of the response is given. Generally

transfer functions are computed for the center of gravity motion. As long as the ship is

regarded as a rigid body they are valid to describe rotary motions at all positions on-board a

ship. However, the transfer functions for the translatory motion of the centre of gravity are

only valid for this particular point. For all other positions on-board a ship the transfer

functions of rotary- and translatory motions have to be superimposed taking into account the

frequency dependent phase angle. This e.g. is done for the calculation of transfer functions for

the vertical motion on the bridge. The statistical evaluation of ship motions regards only the

magnitude of the transfer function, i.e. the amplitude characteristics is employed whereas the

frequency characteristic is discarded because the phase relation of the wave spectrum is not

known.

In principle, the response characteristic is assumed to be linear, i.e. regarding a harmonic

excitation:

• a doubled amplitude of the exciting wave results in a doubled amplitude of the ship

response,

• and the frequencies of the input and the output signal are the same.

The first condition is only partly fulfilled. As long as restoring- and mass forces are dominant

the system behaves linear. When viscous effects arise the system appears increasingly non

linear. This applies to motions in the region of natural frequencies and in particular to roll

motions. Fig. 16 serves to illustrate the effect of non-linear damping coefficients; it shows

transfer functions for the roll motion at beam seas. In the natural frequency region the

magnitude of the ship response is governed by viscous damping. Because damping

coefficients depend not only on the height and frequency of exciting waves but also on the

ship motion itself, they have to be calculated iteratively. Higher waves cause larger ship

motion in connection with increasing damping coefficients, which finally reduces the

magnitude of the transfer function. During the calculation of a transfer function, linearized

terms are used to consider viscous effects. The linearization is done for a mean expected wave

height (variation from 1 to 7m in Fig. 16). In that way motions excited by higher waves are

overestimated because damping effects are underestimated. This can be considered as a little

safety factor in the determination of resonant roll motions. Motions like parametric rolling,

surf riding and broaching are even more difficult to assess, as they are highly non-linear.

Furthermore, within generic seakeeping codes based on strip theory these motions are

decoupled from the exciting forces.

Pareto Optimum Ship Routing, Chapter 4

Fig. 16: Transfer functions for the roll motion

= 90°, Vs = 23kn

For the calculation of transfer functions the desktop-program SEAWAY is utilized,

JOURNÉE (2001). SEAWAY is a strip-theory code, bases on potential theory, and operates

in the frequency domain. It solves a three-dimensional problem by integrating two-

dimensional solutions (strips) over the ship length. In principle the ship is assumed to be a

slender body at zero-speed, i.e. fluid flow velocities in the transverse direction are much

greater than in the longitudinal direction. SEAWAY provides two different methods to

account for forward speed effects referred to as ordinary- and modified strip theory.

Furthermore the program offers different semi-empirical methods to account for the influence

of viscous damping on the roll motion. A comprehensive overview on the program and all

implemented features is given by JOURNÉE and ADEGEEST (2003). Fig. 17 shows the

frame representation of the CMS HANNOVER EXPRESS, as it is used by SEAWAY.

Fig. 17: Sections CMS HANNOVER EXPRESS

Besides transfer functions for the center of gravity motions for all six degrees of freedom and

transfer functions for the motions of selected points, absolute and relative to the free surface,

Pareto Optimum Ship Routing, Chapter 4

two different methods are available to determine response functions for the added resistance

due to waves. The methods are the radiated energy method according to GERRITSMA and

BEUKELMANN (1972) and the integrated pressure method according to BOESE (1970). The

numerical results initially used within the route optimization procedure are:

• The transfer functions of the translatory motions on the bridge. They are used to

calculate the significant amplitudes of lateral and vertical accelerations on the

bridge (loading on the crew by means of accelerations).

• The transfer functions of the motion at a point "10% of LPP behind FP", relative to

the free surface are used to calculate the slamming probability.

• The response functions for the added resistance due to waves are used within the

calculation of the fuel consumption (preferably according to the integrated pressure

method of BOESE which displays a smother distribution over the wave frequency

than its radiated energy counterpart).

Fig. 18 shows some transfer function and response functions for the added resistance due to

waves. From left to right they are: transfer functions for the vertical motion on the bridge,

transfer functions for the relative vertical motion at the bow, and response functions for the

added resistance due to waves. They are given for a wave encounter angle of 180° (head sea)

and for ship speeds of 13, 18, and 23kn. The figure serves to illustrate the principle shape of

transfer functions and response functions.

Fig. 18: Representative transfer functions and response functions for the added

resistance in waves

It is decided to employ a basic but representative model of the ship in waves, as the focus is

put on the capability of the route optimization approach. Because roll motions are strongly

depend on viscous effects it is more complex to integrate them to the model compared to

other modes, e.g. vertical acceleration. A first step can be to apply the appropriate transfer

function with reference to the exciting wave height, cp. Fig. 16. But for the time being it is

decided to take the prediction of natural roll motions out of consideration. However, a method

to account for situations that are suspicious for parametric rolling is suggested in section 4.9.

Regarding the added resistance in waves, there exist meanwhile more sophisticated methods

based on direct pressure integration over the instantaneous wetted surface, taking into account

second order wave drift forces, see e.g. JOURNÉE and MASSIE (2001) These are not applied

here, simply because they are not offered by the employed strip theory code. Furthermore the

employed methods provide satisfactory results as wind forces and current are left out of

consideration yet.

It will be a task of future work to identify suited and fast methods to predict all ship

characteristics that are relevant for a routing decision support.

The employed mathematical model and procedures to describe ship responses in waves are

presented in the following sections.

Pareto Optimum Ship Routing, Chapter 4

4.2 Statistical evaluation of the ship motion in irregular waves

For the description of the sea state a wave spectrum according to PIERSON and

MOSKOWITZ (1964) is used. It describes the distribution of wave energy over harmonic

frequency components and is given as a function of the significant wave height H1/3 and the

average zero crossing period T0. Directional parameters are neglected. Therefore the spectrum

describes long crested waves, swell, encountering from one prevailing direction:

()











⋅

−

⋅⋅

⋅

=−− 4

21.496

exp

03.124 3

ωωω

S. (4.2)

The evaluation of the ship motion uses the 0th and higher order moments of the response

spectrum of the ship. It is the product of the squared transfer function with the wave spectrum.

However, the encounter frequencies of the harmonic components of a wave spectrum are

modified due to the ship speed and the wave encounter angle:

cos

se Vk−= . (4.3)

Since the statistical evaluation uses the response spectrum for the encounter frequency, the

wave spectrum has to be transformed from the wave frequency- to the encounter frequency

representation:

() ()

es S

ωζ

⋅=

. (4.4)

Claiming for an equivalent amount of energy in the frequency bands

∆ω

and

∆ω

e results in

the equation:

(

)

()

ζζ

dSdS ee ⋅=⋅ . (4.5)

From equation (4.5) follows equation (4.6) that can be used for the transformation of the wave

spectrum to its encounter frequency representation:

()

ωζ

e=. (4.6)

Applying the response spectrum obtained by equation (4.4), the nth order moment of the

response spectrum is calculated by:

()

∫

∞

⋅=

eessn dSm

ωωω

. (4.7)

The transformation of the wave spectrum is unproblematic in head- and quartering seas. In

contrast, the denominator of equation (4.6) shows a zero point for encounter angle of less than

90°. This results in singularities in the integrand of equation (4.7). Fig. 19 illustrates the zero-

gradient for following waves. At a wave frequency of approximately

= 0.4 the gradient

dd e becomes zero.

Pareto Optimum Ship Routing, Chapter 4

Fig. 19: Encounter frequency versus wave frequency

But equations (4.6) and (4.7) also show the validity of the following equation to determine the

nth order moments of the response spectrum:

()

∫∫ ∞∞

⋅=⋅=

ωωωωωω

dSdSm n

ese

eessn , (4.8)

with equations (4.1) and (4.4) for the calculation of response spectra.

In this way the integration over the wave frequency avoids singularities in the integrand.

The linearity of the applied method is used within the determination of significant amplitudes

of velocities and accelerations. Wave- and response spectra are considered as a superposition

of harmonic components. For that reason, the calculation of a velocity- and acceleration

response spectra is not necessary. The 0th moment of the velocity response spectrum equals

the 2nd moment of the motion response spectrum:

Smm 2

&. (4.9)

And the 0th moment of the acceleration response spectrum is equals to the 4th moment of the

motion response spectrum:

Smm 4

&& . (4.10)

A detailed derivation of equations (4.9) and (4.10) is given in Appendix 1.

Finally the calculation of the ship motion is achieved by means of a statistical evaluation of

these moments. The significant double amplitude of the ship motion is:

ams 0

42 3

1⋅= . (4.11)

And the mean zero-up-crossing period of the motion is:

022

ππ

. (4.12)

Pareto Optimum Ship Routing, Chapter 4

4.3 Acceleration on the bridge

Generally transfer functions are related to the center of gravity motion. To investigate the ship

motion at another point, the transfer function for this point has to be built. Since the ship is

regarded as a rigid body, the transfer functions for the rotatory motions are unaffected. The

translatory motions of a particular point on-board the ship are described by a superposition of

rotatory- and translatory center of gravity motions, considering the lever to this point and the

angular phase shift of the considered motions.

As motions are regarded to be small, this can be done in the following way:

The vector of the translatory motion at a point P is built by adding the vector of the

translatory center of gravity motion and the cross product of the vector of the rotatory center

of gravity motion and the appropriate lever-vector:

rsss RTP ×+= , (4.13)

with:

()

PPP

Pssss 321 ,,= the vector of the translatory motion of point P,

()

Tssss 321 ,,= the vector of the translatory center of gravity motion,

()

Rssss 654 ,,= the vector of the rotatory center of gravity motion and

()

PPP zyxr ,,= the lever-vector from COG to P.

Each element of the motion vectors contains a phase term j

that has to be considered

within the superposition, cp. JOURNÉE (2000):

ajj e

jeess

⋅⋅= with: j = (1,...,6). (4.14)

In this way the strip-theory code SEAWAY calculates transfer functions for arbitrary points

on-board a ship. Here a point on the bridge, 61m in front of the aft perpendicular, and 44m

above the keel line, is considered. These transfer functions are used to derive the response

spectra for the motions in x-, y-, and z-direction:

() ()

ωζ

Ps ⋅=

. (4.15)

Actually the 4th order moments according to equation (4.10) are used to calculate the

significant amplitudes of acceleration for these motions:

aP ms 4

1⋅=

&& . (4.16)

When optimizing a ship route, the occurring accelerations on the bridge have to stay within a

defined limit. For the time being, thresholds imposed within the European project

SEAROUTES (2003) are used. Thus, the maximum allowed significant amplitude of

acceleration on the bridge in the lateral directions is 0.2g, in vertical direction 0.15g

(g = 9.807m/s2). These constraints are imposed to ensure the well-being of the crew.

Pareto Optimum Ship Routing, Chapter 4

4.4 Slamming

Slamming is defined as the coincidence of two events: according to OCHI it is the emergence

of the bow, 10% of the ship length LPP behind the forward perpendicular FP and its

subsequent dunking above a critical velocity. In JOURNÉE (2000) the critical velocity is

given by:

ppcr Lgs ⋅⋅= 0928.0

&. (4.17)

Both, the emergence of the bow and the velocity of its dunking, depend on the relative motion

between ship and free water surface. That means, the transfer function of the relative motion

for the considered point 10% of LPP behind FP at the keel line has to be calculated.

The relative vertical motion between ship and wave at a point P (xP, yP) and with a wave

encounter angle

is:

()()

ββ

sincos

33 PP yxki

aPrelP ess +

⋅−= . (4.18)

This equation is used within SEAWAY to deduce the transfer function of the relative motion:

()()

ββ

ζζ

sincos

33 PP yxki

relP e

ss +

−= . (4.19)

In addition the amplitudes of the relative motion relP

s3 are adjusted by allowances for static-

and dynamic swell up, i.e. allowances that account for the deflection of the water surface

caused by the forward speed and the oscillating motion of the ship. Within SEAWAY these

allowances are calculated from empirical formulae based on model experiments given by

TASAKI, cp. JOURNÉE (2000).

Transfer functions, according to equation (4.19), and wave spectrum are used to calculate the

response spectrum of the relative motion:

() ()

relP

relPs ⋅=

3. (4.20)

Actually the 0th and 2nd order moments of the response spectrum are needed.

Since the ship motion represents a narrow banded GAUSSian process (narrow banded: the

number of maxima is equal to the number of zero-up crossing events), the response spectra of

the relative motion and of the relative velocity can be represented by a RAYLEIGH7

distribution. Generally, the RAYLEIGH probability density distribution is:



















−

⋅= 2

exp

σσ

relPrelP

Rayleigh

p. (4.21)

Here the variance of the process

2 is the 0th moment of the response spectrum considering

the relative motion and the 2nd moment for the velocity, respectively. The probability of an

emergence of the bow is the probability for a relative motion that exceeds the draft DP at 10%

LPP behind FP:

{}



















−

relPS

PrelaP m

DsP

exp . (4.22)

7 John William Strutt, 3rd Baron Rayleigh, physicist, 1842 – 1919, England

Pareto Optimum Ship Routing, Chapter 4

The probability for a re-entrance velocity that is higher than the critical velocity is calculated

by:

{}



















−

relPS

crrelaP m

ssP

exp &

&& . (4.23)

Presuming statistical independence of these events, the probability of slamming is calculated

from:

{}



















−

−

relPS

slamP

exp &. (4.24)

With a zero-up-crossing period of the relative motion:

relPS

relPS m

=, (4.25)

the number of slams per hour is calculated from:

[]

{}

slamP

hourslamsN

relPS

⋅=

3600

/. (4.26)

The maximum allowed slamming probability within the route optimization is set to 3%. This

threshold as well is taken from the SEAROUTES project.

Pareto Optimum Ship Routing, Chapter 4

4.5 Fuel consumption and load to the main engine

In the following the method for the calculation of the fuel consumption during a journey is

described. It is adapted to the equipment and the hydrodynamic characteristics of CMS

HANNOVER EXPRESS, who is used for the set-up of the route optimization procedure, and

considers the design load condition given in SAMSUNG (1991). Since no shaft generator is

installed, the considered power demand is restricted to the demand for propulsion. Fuel

demand of auxiliary engines is neglected yet.

Thus the operation point of the main engine, i.e. the machine power PS at a particular rate of

revolution rpm, depends on the over-all resistance and the ship speed. For a ship speed

required within the optimization process, it is necessary to control if the operation point is

permitted. Finally the specific fuel consumption for this operation point is determined.

4.5.1 Determination of the over-all resistance

The over-all resistance Rtotal is regarded as the sum of calm water resistance RT and mean

added resistance due to waves RAW:

AWTtotal RRR += . (4.27)

According to the hypothesis of FROUDE8, the calm water resistance consists of:

ADDWVT RRRR ++= , (4.28)

or rather, following the method of HOLTROPP and MENNEN (1984), that is applied here:

()

BAAPPWFT RRRRRkR +++++= 1 . (4.29)

The components are:

RV viscous resistance = (1+k) RF,

k form factor,

RF frictional resistance,

RW wave resistance,

RADD additional resistances = RAPP+RA+RB,

RAPP resistance of appendages (bilge keel, rudder),

RA roughness allowance and still air resistance, and

RB resistance due to the bulbous bow.

Additional the applied method provides the thrust deduction fraction t, the wake fraction w

and the relative rotative efficiency

r, cp. also HOLTROP (1978). The following figure

represents the composite of the calm water resistance.

8 William Froude, engineer, hydrodynamicist, and naval architect, 1818 – 1879, England

Pareto Optimum Ship Routing, Chapter 4

Fig. 20: Calm water resistance

Response functions for the added resistance are applied to determine the added resistance due

to waves. SEAWAY provides two different methods to calculate these response functions.

The method of BOESE (1970) integrates the oscillating pressure on the wetted surface.

GERRITSMA and BEUKELMANN (1972) determine the added resistance by means of the

radiated wave energy. Fig. 21 represents response functions for both methods. The left side

depicts response functions for the method of BOESE, the right side those of GERRITSMA

and BEUKELMANN, respectively.

Fig. 21: Response functions for the added resistance due to waves

The shape of the response functions in the area of

= 1.2 to 1.4 results from a numerical

imperfection of the method used within strip-theory to calculate potential mass- and damping

coefficients, commonly known as irregular frequencies9. Thus, for the determination of the

added resistance, the functions according to BOESE are preferred. The quality of the

numerical results is better, i.e. less affected by numerical irregularities.

9 Irregular frequencies are caused by singularities in the GREEN’s function. They occur during the determination

of hydrodynamic coefficients. Because they are strongly dependent on the frame shape and the frequency of a

motion, they are difficult to avoid, cp. JOURNÉE and ADEGEEST (2003). For further information see

Appendix 2.

George Green, mathematician and physicist, 1793 – 1841, England

Pareto Optimum Ship Routing, Chapter 4

The mean added resistance in waves is determined from:

()

∫

∞

⋅⋅=

ωω

ζζ

AW , (4.30)

with: 2

the response function for the added resistance due to waves.

Fig. 22 serves to illustrate the mean added resistance as a function of the ship speed for a

significant wave height of 7m and a peak period of 14s in head sea. At a ship speed of 15kn

the mean added resistance amounts to the same magnitude like the calm water resistance.

Even if it decreases at a higher ship speed, the total resistance at this sea state assumedly

increases the capacity of the main engine. This will be checked within the determination of

the engine operating point.

Fig. 22: Added resistance due to waves

4.5.2 Propeller characteristics and operation point

The “full scale HANNOVER EXPRESS” propeller characteristic is derived from the

propeller characteristic of a similar model scale propeller given by YASAKI (1962).

Tab. 4: Propeller data

propeller data HANNOVER EXPRESS model propeller

diameter D 8.25m 0.25m

number of blades Z 6 6

slope ratio P/D 1.036 1

area ratio Ae/A0 0.795 0.7

The full-scale propeller characteristics are provided by applying corrections for the

REYNOLDS10-number and roughness to the model propeller data, according to ITTC (1978).

Fig. 23 shows the propeller characteristic that is used for the CMS HANNOVER EXPRESS.

10 Osborne Reynolds, physicist, engineer and mathematician, 1842 – 1912, England

Pareto Optimum Ship Routing, Chapter 4

Fig. 23: Full-scale propeller characteristic

Thrust, torque, and propeller advance velocity are represented dimensionless:

⋅

=, (4.31)

42 Dn

KTS

=, (4.32)

52 Dn

KQS

=, (4.33)

with:

J advance ratio,

KTS thrust coefficient,

KQS torque coefficient,

VA advance speed,

n revolution frequency of the propeller,

D propeller diameter,

density of sea water, and

0S open water efficiency, QSTS KJK ⋅⋅

The determination of the operating point of the propeller is accomplished according to thrust

identity. Therefore the ratio KTS /J 2 has to be built:

()( )

211 DVwt

totalTS

−−

=. (4.34)

Following, equation (4.34) is multiplied by J 2 and depicted in the propeller characteristic. The

intersection with the curve for KTS represents the operating point of the propeller. Fig. 24

illustrates the determination for two different ship speeds, 10kn and 23kn, at a sea state of 6m

significant wave height, 12s peak period, and in head waves.

Pareto Optimum Ship Routing, Chapter 4

Fig. 24: Determination of the operating point

By means of the operating point of the propeller, the required power and rate of revolution of

the main engine are determined from:

()

⋅

−⋅

=1, (4.35)

KDn

ηη

=, (4.36)

with:

PS shaft power,

S shaft efficiency (0.99),

R relative rotative efficiency.

4.5.3 Feasibility of the operation point of the main engine

Since the ship velocity is a given value within the optimization, it has to be controlled if the

operating point of the main engine is in accordance with the main engine characteristic,

MAN B&W (2000). Fig. 25 represents the normalized main engine characteristic. It shows

the permitted operation area for the main engine of CMS HANNOVER EXPRESS, a

9K90MC. This area is bordered by 4 lines:

• a boundary line for engine overload caused by the limited ability of the turbo chargers

to compress air (~ n 2),

• a boundary line for the maximum allowable mean effective pressure (~ n),

• a boundary line for the nominal power (MCR, ~constant), and

• a boundary line for the maximum rate of revolution at MCR. At most 5% above the

nominal rate of revolution.

Former two lines specify the available power for an operation point below the nominal rate of

revolution. They are used as the characteristic of the control unit of the main engine.

Pareto Optimum Ship Routing, Chapter 4

Additional the curve for the power demand of the propeller under nominal operating

conditions (~ n3) and the specified MCR point are shown (SMCR, maximum rating for

continuous operation required by the yard or the owner).

Fig. 25: Main engine characteristic

The lower limit of the main engine operability is not explicitly modeled. Since the ship uses a

fixed pitch propeller and the ship speed is bounded below by the perturbation of the velocity

profile at 10kn, it is not possible to derate the main engine. Furthermore the masters

individual preferences in operating the main engine, e.g. to avoid soot, can be modeled by

increasing the minimum allowed ship speed. For the time being these aspects are neglected, as

for the purpose of this study a quiet simple main engine model is absolutely sufficient.

4.5.4 Specific fuel consumption

Fig. 26 shows the specific fuel oil consumption sfoc according to MAN B&W (2000).

Fig. 26: Specific fuel oil consumption, sfoc

Pareto Optimum Ship Routing, Chapter 4

It is dependent on the rate of revolution of the main engine. The specific fuel consumption

spef follows from:

[]

kWP

hkW

sfoc

spef S0000001

⋅























. (4.37)

The over-all fuel consumption for a whole journey is calculated from:

[] []

htspeftnconsumptiofueloverall

arrival

departure

∑∆⋅= , (4.38)

with: ∆t: required time for a route segment.

The over-all fuel consumption of a journey is one objective of the route optimization.

Lubrication oil and fuel consumption of the auxiliary engines are not considered. To keep

things simple, the consumption of lubrication oil is regarded to be proportional to the fuel oil

consumption of the main engine, therefore savings by route optimization are proportional to

fuel savings as well. The fuel consumption of the auxiliary engines is time dependent. Faster

routes will be the more fuel saving ones. Actually this type of fuel consumption is not known

exactly, e.g. it would be necessary to have the number of reefer container on-board. An

extension of the model is straightforward but nonetheless omitted as it is not necessary for

this study.

4.6 Practical calculation of ship responses

To accelerate the optimization process, the determination of ship responses due to waves on a

particular route is pre-processed as far as possible. Since the equations to determine the

moments of the response spectra (4.8 – 4.11) are linear in terms of the squared significant

wave height,

ii mHm ⋅= 2

1, (4.39)

it is reasonable to calculate a ship response due to a significant wave height H1/3 of 1m and to

scale it afterwards with the actual wave height. Fig. 27 shows PIERSON-MOSKOWITZ

spectra for H1/3 = 1m.

Fig. 27: PIERSON-MOSKOWITZ spectra for H1/3 = 1m

Pareto Optimum Ship Routing, Chapter 4

All of the spectra are equal in terms of m0, the area under the curve, i.e. they contain the same

amount of wave energy. Whereas the peak is more pronounced when peak periods Tp increase

or simply “the longer the waves get”. These spectra are used to calculate a “ship response data

base”. Therein the ship responses for H1/3 = 1m and peak periods TP in steps of 1s are stored.

The actual ship responses are calculated by interpolations for TP, encounter angle

, and ship

velocity VS and by scaling the appropriate values with the appearing wave height.

Applying the moments of the response spectra according to a significant wave height H1/3 of

1m, acceleration, slamming probability, and added resistance due to waves are accounted by:

[]

ap 807.9

⋅⋅

&& , (4.40)

{}

[]

%100

exp

⋅































+⋅−=

relps

slamP &, (4.41)

[]

RAW

AW 1000

1⋅

=. (4.42)

Additional they are transformed to the units of the threshold values used within the route

optimization.

Furthermore, the ship responses can be depicted as functions of the significant wave height

H1/3 and the peak period Tp. Fig. 28 – Fig. 34 serve to illustrate the dependency of the ship

response on these two crucial factors for selected operating conditions. The figures cover a

range of Tp = 0 – 30s and H1/3 = 1 – 15m. Furthermore a borderline, the wave-breaking limit,

for a maximum significant wave height H1/3 max is given. Here the physically maximum wave

height Hmax is assumed to be 1/10 of the wave length and to be approximately the double of

the significant wave height H1/3. The area of interest is on the left side of the borderline.

Fig. 28: Significant amplitudes of vertical

acceleration on the bridge for different

encounter angles

Fig. 28 shows the vertical acceleration on the

bridge for wave encounter angles of 180°,

140° and 90°. Expectedly the maximum of

the response shifts from longer periods in

head seas towards shorter ones in quartering

seas. In head seas a maximum response can

be observed at a peak period about TP = 16s,

corresponding to a wave length of 400m.

Pareto Optimum Ship Routing, Chapter 4

Having in mind the ship length LPP = 281m and the natural period of the pitch motion of

8.2s, this seems quite plausible. For quartering seas the maximum shifts towards a peak

period TP = 8.5s, the natural frequency of the heave motion. The influence of the encounter

angle is relatively small at longer waves (TP > 14s), whereas in shorter waves (TP < 14s) a

significant increase of the ship response can be observed when taking course towards

quartering seas.

Fig. 29 represents the vertical acceleration on the bridge in head seas for a ship speed of 16

and 23kn, the design speed. Obviously only a distinct reduction of the ship speed is suited to

clearly reduce accelerations. The little shift of the maximum responses from Tp = 16s to 15s is

caused by the reduction of the wave encounter frequency due to the reduced ship speed.

Fig. 29: Significant amplitudes of vertical acceleration on the bridge

for different ship speeds

Fig. 30: Slamming probability for

different encounter angles

Fig. 30 illustrates numerical results for the

slamming probability at design speed for

wave encounter angle of 180°, 140°, and

100°. Course variations around 180° do not

show any notable influence on the slamming

probability. Finally at a wave encounter

angle of 140° a clear effect can be observed.

As expected, quartering seas are not able to

excite any slamming, even not at an

Pareto Optimum Ship Routing, Chapter 4

encounter angle of 100° where we have still a ship response according to the transfer function,

but physically no exciting waves are possible. Nevertheless, in contrast to the accelerations on

the bridge, the slamming probability can be reduced by changing course to more quartering

seas.

Similar to the responses for the acceleration on the bridge presented in Fig. 29 there is a

distinct reduction of the ship speed necessary to clearly reduce the slamming probability.

The mean added resistance due to waves RAW for head sea conditions at design speed is

presented in Fig. 31. It shows maximum values for a peak period Tp = 15s. This corresponds

to a wave length of 350m, about 70m more than the ship length. Supported by Fig. 21 it can

be seen that the added resistance is less affected by changes of the wave encounter angle and

increases when speed is reduced (compare Fig. 22). Surely the decrease in the calm water

resistance overcompensates the increase of the added resistance due to waves at lower speed.

Fig. 31: Added resistance due to waves

Further on Fig. 32 presents the required brake power of the main engine at 22kn and at design

speed of 23kn. The line for 41MW marks the maximum continuous rating of the main engine

(MCR) and cannot be exceeded.

Fig. 32: Brake power at reduced and design speed

Pareto Optimum Ship Routing, Chapter 4

Fig. 33 illustrates the specific fuel oil consumption for these two ship speeds.

Fig. 33: Specific fuel oil consumption at reduced and design speed

Finally Fig. 34 shows the specific fuel consumption according to these operating conditions.

The trembling shape of the lines is caused by the discretization of the propeller model,

established as a good compromise between accuracy and computational effort.

Fig. 34: Specific fuel consumption at reduced and design speed

4.7 Avoiding irregular frequencies

Irregular frequencies are a well-known numerical disturbance during the determination of

hydrodynamic coefficients e.g. used within strip-theory methods. If no access to the source

code of the program is given, they are hardly to avoid, cp. Appendix 2. In that case they can

be identified by a plausibility check, and anyhow the numerical result can be emended. Fig.

35 shows a response amplitude operator for the vertical motion on the bridge. The left part of

the figure shows the transfer function like calculated by the strip-theory program. Physically

there is no reason for the peak at a wave frequency 1.3, just as the shape of the transfer

function between

= 1.0 – 1.3 can be doubted. Assuming that the ship response at

frequencies above

= 1 (corresponding wave length 60m) cannot be higher than the one at

= 0.8 (corresponding wave length 100m, the natural frequencies of the heave motion lies at

= 0.74 and for the pitch motion at

= 0.76) it is reasonable to smooth this transfer function

by setting the ship responses to zero for

> 1.

Pareto Optimum Ship Routing, Chapter 4

Fig. 35: Transfer function before and after smoothing

Fig. 36 shows resulting ship responses like presented in section 4.6. Again the left part

represents the calculation using the disturbed transfer function, the right part the smoothed

transfer function respectively. The extreme amplification of responses, in particular in shorter

waves, is caused by 4th order moments used when calculating the significant amplitudes of

acceleration.

Fig. 36: Ship responses, disturbed by irregular frequencies (left)

and corrected(right)

A similar post processing is applied to all transfer functions before calculating the “response

data base” whenever it is necessary.

4.8 Add-on for motion sickness incidence

Beside the recommendation given within the SEAROUTES project, the Institute of Naval

Medicine gives a threshold value for a maximum vertical acceleration, PINGREE (1988). If

no data for the frequency of excitation and its duration are available, a value of

1.75m/s2 RMS (root mean square) is suggested as the approximate level at which disturbances

may occur to tasks such as writing, control operation, etc.

The RMS of a harmonic function a with the amplitude a

is defined as:

Pareto Optimum Ship Routing, Chapter 4

()

adatda

RMS

707.0)sin(

11 2

2⋅=== ∫∫

φφ

. (4.43)

The definition of RMS can also be applied to a stochastic process a, and in this case it equals

m0, the square root of the area under curve of the motion spectrum that describes this process,

i.e. the square root of the variance

2 of the stochastic process:

=== ∫0

1mtda

RMS

. (4.44)

The investigations on motion sickness are generally conducted with harmonic excitation.

Since stochastic motions of the same RMS as their harmonic counterpart show much higher

amplitudes, an equivalent harmonic excitation is built by using the significant amplitude of

the ship response. Being 0

1ma =, the significant amplitude of the ship response, the

equivalent RMS according to KRAPPINGER (1983) is calculated by:

707.0414.1 0amRMS ⋅=⋅= . (4.45)

Therefore the threshold value given within SEAROUTES of 0.15g for the significant vertical

acceleration equals 1.0m/s2 RMS. Recommendations applying a1/10, the mean value of the

one-tenth highest amplitudes, would result in 1.3m/s2 RMS.

Because of these discrepancies of the threshold values and the fact that thresholds represent

only a rough description of the effect of vertical accelerations on the crew, a refined method is

desirable. That is why the concept of motion sickness incidence MSI as given in ISO 2631 is

introduced to route optimization. The MSI gives the percentage of people that suffer from

severe discomfort dependent on the RMS of acceleration, the frequency, and the duration of

the excitation (in the following simply RMS is used for RMS of acceleration). The

investigations of KRAPPINGER (1983) provide a mean value applied as threshold where

kinetosis appears. It is dependent on the frequency of excitation and acts on the assumption

that log(RMS) is GAUSSian distributed:

)(log73.2)(log36.487.0 ff ++=

, (4.46)

with:

: mean value applied as threshold where kinetosis appears,

f: frequency of excitation.

E.g. 10µ is the RMS value at frequency f for a 2h exposure where 50% of the crew suffers

from motion sickness. The variance of log(RMS) is independent on the frequency, i.e.

221.0

. So that the distribution function becomes:

[]

()

()()

RMSd

RMS

RMSF

RMS

log

)(log

exp

)(log

∫

∞−











−−

⋅

σπ

. (4.47)

Expressed in percent this distribution function is known as

[]

)(log100 RMSFMSI ⋅= .

For durations of exposure different from 2h an equivalent RMS for a 2h exposure is calculated

according ISO 2631, assuming that .

2consttRMS =⋅ :

)(

)2(

2ttRMS

htRMS ⋅

== . (4.48)

Fig. 37 shows isolines for 10% MSI as a function of frequency and duration of exposure.

Pareto Optimum Ship Routing, Chapter 4

Fig. 37: 10% MSI isolines

The left diagram of Fig. 38 depicts the isolines of MSI for a 4h exposure at 180° wave

encounter angle and 23kn ship speed. The right diagram is the counterpart for the significant

vertical acceleration. Obviously the concept of applying a threshold value to the vertical

acceleration underestimates the load on the crew. Furthermore the application of MSI is more

meaningful to assess the consequence of vertical acceleration as it includes the frequency and

the duration of the exposure.

Fig. 38: MSI and significant amplitudes of vertical acceleration on the bridge

Pareto Optimum Ship Routing, Chapter 4

4.9 Add-on for parametric rolling

Auto-parametrically excited roll motion, or shortly parametric rolling, was initially classified

as a cock-and-bull story of seafarers to excuse ship accidents caused by human failure.

Nowadays the mechanism that generates these extreme roll motions is known and the ship-

wave parameter constellation that affords parametric rolling is identified. Actually strip-

theory and transfer functions, as applied here for a route optimization, are not capable to

predict parametric rolling or a probability of such an event. However, it should be possible to

identify situations that are suspicious for parametric rolling. Therefore this investigation puts

its focus on the impact on an optimization result, if such situations are avoided. To set up a

routine to identify situations that are susceptible to parametric rolling two contributions to the

discussion at the 13th SNAME ad hoc panel, SHIN et al. (2005) and KRUEGER et al. (2006),

have been combined. Both are presented below as far as they are used here. Actually both

approaches step much further, as they claim to be able to predict extreme roll motions or

capsizing in rough seas. In the following they are simply referred to as SHIN and KRUEGER

approach.

4.9.1 Introduction of criteria from the Shin approach

SHIN et al. (2005) point out that dependent on the relation between angle, length, and period

of the encountering waves on one hand and pitch and natural roll period on the other hand,

certain relations of these parameter are suspicious for parametric rolling.

Considering a ship in heading or following waves, there seem to be no external forces that

may induce roll motion. But this motion phenomenon depicts the difference between

mathematical model and real life. Even long crested waves are superimposed by wave

components with directions that are different from the prevailing direction. Wind load is a

further factor, so that at least a bit of roll motion is always induced and present. As long as

damping forces exceed exciting forces the amplitudes of this roll motion do not attract

attention. During parametric rolling the frequency of exciting forces matches the natural roll

frequency or a multiple of it. As damping forces are relatively low, the amplitudes of the roll

motion increase dramatically. A further amplification of the process is caused by means of

oscillating restoring forces due to a periodically changing metacentric height GMT. Changes

in the shape of the water line provide maximum uplifting forces during a ship sails in a wave

trough, whereas minima of GMT occur while a ship sails on a crest. In particular modern

container vessel and other ships with pronounced bow flare and extreme stern overhang are

affected.

A typical mechanism is as follows:

Assuming a ship is sailing in head sea, the wave length matches approximately LPP and the

ship is pitching down its bow into an approaching wave with a little heel angle to starboard. In

this situation the restoring forces for the roll motion are high, and if the wave is high enough

it takes the ship by its bow flare, raises the bow, and powerfully rolls the ship to portside

direction. The ship reaches the upright position during the wave crest passes the midship

region. Now as restoring forces are small the ship rolls unresisted onward to portside. Again it

pitches down the bow into the next wave crest, but with a much higher heel angle than on the

starboard side before. If this process has started and matches the natural roll period of the ship

only instantaneous course changes can avoid further increasing of roll amplitudes and finally,

the loss of cargo or a capsizing of the ship. Reducing speed is an adequate countermeasure if

sea states that may become critical are detected early enough. Because this course of motion

typically takes two pitch periods associated with one roll period it is commonly called the

Pareto Optimum Ship Routing, Chapter 4

2:1 type of parametric rolling. A more or less symmetric roll motion, i.e. large roll angles to

starboard and port side, characterizes this type. Another kind of parametric rolling is the

1:1 type. In this case one roll period is associated to one pitch period. An asymmetric roll

motion, i.e. large roll angles to starboard or port side, characterizes the 1:1 type. Dependent

on the ratio of oscillating restoring- and oscillating exciting forces the 1:1 or the 2:1 type of

parametric rolling is met. In this regard the ship size, ship speed, wave direction (i.e. head or

following seas), and the metacentric height GMT are deciding parameters. SHIN et al. (2005)

do not explicitly addressed the 1:1 type, even if their numerical model is able to predict it.

According to SHIN the following parameters are characteristic for parametric rolling:

• Heading or following sea.

• The wave length matches approximately the ship length.

• The time for one roll period takes two pitch periods

and matches approximately the natural roll period.

• A required wave height brings about exciting forces

that exceed the damping forces.

Regarding a route optimization, the first two parameters can directly be obtained from

weather forecast data and the course angle of the ship. Wave direction and course angle

provide the wave encounter angle and the wave length according LEWIS (1998, Vol. III) is:

2TT

LPwave ⋅=

. (4.49)

The third parameter, the natural roll period, is strongly affected by the wave height and the

wave period. Fig. 39 shows the metacentric height GMT for the CMS HANNOVER

EXPRESS in a STOKES-III11 wave of 7.5m wave height and 283m wave length,

approximately the ship length. Associated natural roll periods are assessed by the formula for

the natural period at calm water

eGMg

=, (4.50)

with the radius of inertia xx

k of 12.4m. This figure simply should illustrate the strong effect of

waves on the stability and the roll motion characteristics of the ship. The phase, given at the

abscissa, represents the shift of the wave crests. Phase = 0 or 1 means wave crests at FP and

AP. Phase = 0.5 represents a wave crest amidships.

Fig. 39: GMT and roll period

11 George Gabriel Stokes, mathematician and physicist, 1819 – 1903, Great Britain

Pareto Optimum Ship Routing, Chapter 4

In this position of the wave GMT becomes even negative which causes the gap in the curve of

the roll period on the right side, an indicator for vanishing upright stability.

At low wave heights the effects on the calm water natural roll frequency may be negligible.

While in a parametric rolling situation, where the wave length matches the ship length, the

modification of the natural roll frequency is dependent on the wave height. Are waves high

enough, the “natural frequency” becomes almost the wave encounter frequency.

SHIN detects situations that are susceptible for parametric rolling by transforming differential

roll equation to the MATHIEU12 equation and identifies bounded and unbounded solutions.

However, this method is too extensive in computing time to include it to a route optimization

where thousands of weather situations have to be evaluated. But latter two parameters, the

natural frequency and the required wave height, can be condensed into one, an overcritical

wave height. For this purpose, the approach suggested by KRUEGER is used.

4.9.2 Additional criteria provided by the Krueger approach

KRUEGER et al. (2006) carry out capsize simulations for various ship types and load

conditions. At this, roll angle of more than 50° are ranked as a capsize event. They conclude

that “the energy introduced into any specific hull form in a sea state may be expressed by the

alteration of the areas below the righting levers at trough and crest condition” and “that the

most important phenomena leading to large rolling angles can be directly accessed by linking

the righting lever changes between crest and trough condition to the minimum stillwater

stability requirements”. The method for the stability assessment proposed by KRUEGER does

not distinguish between capsizing due to e.g. stability loss in following seas or due to

parametric rolling (both 1:1 and 2:1 type). It is posed to assess a required safety level for ships

in a group of probable sea states. Here the method is used to assess a critical wave height.

Higher waves characterize a sea state that becomes suspicious for parametric rolling.

For this purpose the alteration of the lever curves for a design load condition at GMT = 0.73m

are studied. Fig. 40 shows two investigated cases. They are lever curves for calm water,

trough- and crest condition at a wave length of 283m and 5m or rather 10m wave height. The

figures nicely illustrate the increase in the alteration at higher waves. The shape of the calm

water lever curve is strongly dependent on the metacentric height GMT, whereas the

alterations of the lever curve are affected by the wave height and the hull shape.

Fig. 40: Lever curves for calm water, trough and crest condition

12 Émile Léonard Mathieu, mathematician, 1835 – 1890, France

Pareto Optimum Ship Routing, Chapter 4

Although the approach does not distinguish between capsizing due to stability loss in a wave

crest and due to parametric rolling conditions, and although capsize is the final result, i.e.

severe parametric rolling may also occur without capsizing, this approach can be used to

assess a threshold value for a critical wave height.

KRUEGER et al. use a capsize index ak to decide whether a ship can be operated safely in a

given sea state or not. Within this index they use the following capsize probability based on

their simulations and applicable to all ships:

max401 shsa ⋅−= , (4.51)

i.e. 1=a represents 100% capsize probability and 0=a is 0% respectively.

The factors under the root are as follow:











≤≤











≤≤

2maxmax1

0maxmax0

2maxmax0

maxmax

max

140401

040400

140400

4040

khhif

hhif

khhif

kaaif

aaif

kaaif

diff

With 40a being the area from 0 - 40deg under the lever curve at limiting stability conditions

and diff

a40 the difference of the areas of the lever curves at the current load condition. With

maxh and diff

hmax the same is put for the maximum uplifting levers. According to the

recommendation of KRUEGER the factors k1 and k2 are set to 0.7.

The left side of Fig. 41 shows the lever curves for GMT = 0.73m (design load condition, calm

water, crest and trough condition) and for GMT = 0.37m (minimum or limiting stability

condition according to IMO resolution A.749(18), code on intact stability). On the right side

the capsize probability dependent on the wave height is shown. The wave length is held

constant and matches approximately the ship length.

Fig. 41: Lever curves and capsize probability

Pareto Optimum Ship Routing, Chapter 4

For 10m wave height the capsize probability for a ship sailing at limiting stability conditions

is about 24%. For the design load condition considered here the capsize probability decreases

to 12%, i.e. the increased metacentric height of GMT = 0.73m already enhances safety as

regards the capsize probability.

The diagram on the right side of Fig. 41 is used to set a threshold value hcrit for the wave

height. For example: parametric rolling needs a group of consecutive waves showing an

overcritical wave height. That is why an event that causes a capsize probability of, reasonably

assumed, 12% should be tolerated only once per hour. This implies a threshold mhcrit 10=.

Assessing the wave height to be the double wave amplitude a

2= leads to the probability

that a critical amplitude is exceeded by:

()











−

crit

crita m

hP 2

exp

. (4.52)

Considering the relation of significant wave height and 0th moment of the wave spectrum:

1mh =, (4.53)

the probability of exceeding a critical wave height becomes:

()











−

=> 2

exp h

hhP crit

crit . (4.54)

The number of times per hour the wave height exceeds the critical wave height is calculated

by:

[]



















−

exp

3600

hourpernumberE crit

enc

. (4.55)

For the time being the value E is set to “1”, i.e. a threshold of maximum one exceeding per

hour is required. Therefore the tolerated probability P calculated according to equation (4.54)

also depends on the encounter period Tenc given in equation (4.55), i.e. the tolerated

probability is dependent on the ship speed and the wave encounter angle, cp. equation (4.3).

Tab. 5 gives an example for the HANNOVER EXPRESS at design load condition,

GMT = 0.73m. The critical wave height hcrit is set to 10m, related to a capsize probability of

12%, according to Fig. 41. A joint capsize probability can be calculated from the probability

calculated according to equation (4.54) multiplied with the probability of 12%.

Tab. 5, Joint capsize probability for TP = 16.5s, Vs = 23kn,

hcrit = 10m, head- and following seas

encounter angle 180°, head waves 0°, following waves

encounter period Tenc 10.0s 46.6s

probability for waves above hcrit 0.0028 0.013

capsize probability 0.12 0.12

joint capsize probability 3.3 * 10-4 1.6 * 10-3

The following section describes the assembly of the approaches and the set-up of the

procedure to assess parametric rolling within route optimization.

Pareto Optimum Ship Routing, Chapter 4

4.9.3 Combination of the approaches

Fig. 42 to Fig. 45 serve to illustrate the set-up of the parametric-rolling-module used within

the route optimization. They show calculations for CMS HANNOVER EXPRESS at head

seas and 23kn ship speed. The natural roll period in calm water is 29s, GMT is 0.73m, design

load condition. The wave peak period TP has been varied from 2s to 29s. The significant wave

height H1/3 has been set to a fortieth of the appropriate wave length calculated by equation

(4.49). The left diagram of Fig. 42 shows the wave length and significant wave height versus

the peak period. They are used as an input parameter for this exemplary variation study. The

right diagram shows the corresponding pitch period of the ship and the encounter period of

the waves. At a peak period of 16.5s the zero-up-crossing-pitch period matches the wave-

encounter period. The corresponding wave length is 290m, the significant wave height is

about 10m. Apparently the area around TP = 16s is suspicious for parametric rolling because

the parameter constellation mentioned in section 4.9.1 is fulfilled.

Fig. 42: Variation of the wave peak period

Fig. 43 depicts the ship response at FP, COG, and AP. Absolute motions and also vertical

motions relative to the wave surface are shown. As expected, motions increase with

increasing wave height and wave period. The higher amplitudes at the forward perpendicular

are caused by a center of the waterline area lying abaft, as typical for these ships (in calm

water the center of the waterline area is approximately 10 m behind of the center of gravity).

The magnitudes of the motions appear credible.

Fig. 43: Absolute and relative motion at FP, COG, and AP

Pareto Optimum Ship Routing, Chapter 4

The previous section describes how to set up a threshold value for the wave height. As a start

position for further investigations it is put to hcrit = 10m, representing a capsize probability of

12%. Such an event will be tolerated only once per hour.

Another attempt, aiming to increase the influence of the pitch motion, uses the relative

displacement between ship and free surface at FP, COG, and AP. Instead of a maximum

tolerable wave height hcrit a threshold for maximum relative motion

rel crit is posed. Similar to

equation (4.55), the number of times this threshold will be exceeded is calculated by:

[]



















−

rels

critrel

enc mT

hourpernumberE

exp

3600

. (4.56)

As aforementioned for equation (4.55) the value E is set to “1”, i.e. here too a threshold of

maximum one exceeding per hour is required.

If no ship motion is regarded a wave height of 10m causes approximately a relative

displacement of 5m at FP, COG or AP. The consideration of the hull shape (cp. Fig. 44)

suggests that a wave elevation at the stern may have a higher impact on the ship as the stern

overhang is reached more quickly. The bow shows steeper angles of entrance in the vertical

direction. For this reason a moderate adaptation of

rel crit seems reliable.

Fig. 44: Ship lines at FP and AP

For the time being the tolerable relative motion between ship and water surface

rel crit

is arbitrarily set to:

• 6m for FP

• 5m for COG

• 4m for AP

Pareto Optimum Ship Routing, Chapter 4

According to equation (4.55) or (4.56) the exceeding of a threshold value E is calculated if:

• the wave length is in the range of the ship length,

for a start 5.175.0 ⋅<<⋅ ppwavepp LLL is used,

• the wave encounter angle β approaches head or following sea,

for a start a range of °<°+ 1805

or °<°− 05

is used,

• the encounter period Tenc matches the pitch period Tpitch

for a start

(

)

1.0⋅<− pitchpitchenc TTT is used,

otherwise E is set to zero.

Maxima of exciting forces are detected by comparing encounter- and pitch periods (third

bullet point). The natural roll period is not considered, as it is strongly dependent on the wave

height. In this way the method is able to identify wave conditions that can excite both, the 1:1

or the 2:1 type of parametric rolling. The ranges, so far, are arbitrarily set and need further

investigation and validation by model tests or numerical simulation, e.g. following

IMO (2007) the range for the wave encounter angle

has surely to be increased. For the time

being they are used for the route optimization as aforementioned to apply and investigate the

proposed add-on.

In the way described above, sea states and ship responses are evaluated during the route

optimization. Parameter constellations that are suspicious for parametric rolling are

recognized and assessed by the frequency of exceeding thresholds for the absolute and the

relative wave height at FP, COG, or AP, respectively. The left part of Fig. 45 shows the

exceeding E for the variation presented in Fig. 42. In the right part, these functions are shaped

considering critical parameter constellations, i.e. E is evaluated only if wave length, pitch-

and wave period, and encounter angle indicate that parametric rolling is probable. Although

the threshold for the significant amplitude at FP is increased, the most events where the

threshold value is crossed happen at FP.

Fig. 45: Threshold exceeding generally and for parametric rolling conditions

Section 6.5 introduces first investigations upon the developed parametric-rolling-module in a

route optimization. Furthermore comparative studies on the utilization of hcrit and

rel crit for

the calculation of E are presented.

Pareto Optimum Ship Routing, Chapter 4

It is very important to line out that the above developed add-on should be a

first step to identify sea states and courses that are suspicious to parametric

rolling. It does not cover all the possible situations that are able to cause

capsize or extreme roll motions. For this purpose e.g. a much larger range of

wave encounter angle has to be considered.

4.10 Considered ship responses

Tab. 6 summarizes the characteristics of CMS HANNOVER EXPRESS considered in this

study. Besides these criteria, the table gives information about the focused locations on the

ship and the name of the parameter within the mathematical model. Additionally it is

specified if these parameters are sensitive to variations of the ship speed VS, the peak period

of the encountering waves TP, or the encounter angle

. Cells of Tab. 6 that are left blank

display that there is no pronounced sensitivity. Values or ranges of values display that there

are peaks or areas of extreme responses. All parameters show sensitivity within at least two of

the variations. As a start, all of them appear to be crucial for a route optimization.

Tab. 6: Considered ship attributes and summary of ranges of extreme responses

criteria at / of parameter name Vs [kn] Tp [s]

[deg]

accelerations bridge

xacc 18 / 24 9.5 150 / 170

bridge

yacc - 12.0 110

bridge

zacc - 10.0 100

motion sickness bridge MSI - 10.0 100

slamming probability 10% behind FP slprob 24 13.0 130 - 180

parametric rolling FP nph@fp

24 14.0 -16.0 175 - 180

naming convention is COG nph@cog

- 15.0 -16.0 175 - 180

given in section 6.5.1 AP nph@ap

10 12.0 –13.0 175 - 180

on page 91 wave heights nph4hcrit - - 175 - 180

added resistance ship Raw 14 - 16 13.0 - 15.0 130 - 180

Pareto Optimum Ship Routing, Chapter 5

5 A multi objective, stochastic approach for the route optimization

The first section of chapter 5 explains the setup of the route optimization procedure.

Following, section 5.2 addresses the optimization method. Here the application of

deterministic and stochastic optimization methods is discussed. The final section 5.3 focuses

on the amount of probable fuel saving for achievable estimated times of arrival.

5.1 Optimization setup

Tab. 7 illustrates parameters being used for an initial setup of the route optimization. The first

letter of a free variable s or v shows if it refers to a spatial or a velocity perturbation of the

route. Constraints are marked by a capital C followed by the threshold they refer to.

Objectives start with capital O followed by the parameter to be optimized. Constraints for

acceleration on the bridge and slamming use thresholds that are used within the

SEAROUTES project. Later on, during the course of the sensitivity studies presented in

chapter 6, further constraints for parametric rolling and motion sickness incidence will be

included. Finally a third objective to minimize the risk of unfavorable weather changes will

be established.

Three types of free variables are applied: shift variable for the spatial shape of the route, shift

variable for the shape of the velocity profile during the journey, and a variable that performs a

delayed departure. Four constraints are imposed for lateral and vertical accelerations and for

the slamming probability. Three additional constraints are needed to control the feasibility of

a route variant in terms of main engine operability, time range of available weather data, and

land collision. Finally, two objectives are applied, for the estimated time of arrival ETA13 and

for the fuel consumption, both of them are to be minimized, and furthermore ETA has to

match the schedule.

Tab. 7: Parameter of the optimization setup

parameter name value utilization

free variable s1 to s7 -1.0 to 1.0 spatial route perturbation

v1 to v9 0 to 1.0 velocity profile perturbation

RTD 0 to 12 later departure

constraint C_xacc 0.2 g upper limit for lateral acceleration

C_yacc 0.2 g upper limit for lateral acceleration

C_zacc 0.15 g upper limit for vertical acceleration

C_slprob 3% upper limit for slamming probability

C_spef dep. operability of main engine, etc.

C_fuel dep. over all fuel consumption

C_eta dep. time range of available weather data

objective O_eta --- arrival on schedule

O_fuel --- minimum fuel consumption

13 In the strict sense, the duration of a journey and not ETA is minimized. However, all route optimizations start

at the time t = 0h and aim to optimize the time that is needed to arrive at a destination point. For this reason ETA

is used as a synonym for the duration of a journey.

Pareto Optimum Ship Routing, Chapter 5

In principle it is possible to simply use the estimated time of arrival ETA and the fuel

consumption FUEL as objectives. But it turned out that tailored objective functions are able to

control, speed up, and improve the optimization result:

()

tETAvaluesetETAETAObjective ∆−−= , (5.1)













∆−

⋅−⋅= tETAvalueset

ETA

FUELFUELweightFUELObjective . (5.2)

Equations (5.1) and (5.2) show the objectives, applied in most of the presented optimizations.

Equation (5.1) uses the squared difference between ETA and a set value, representing an

arrival on schedule. The objective function considering the fuel consumption (5.2) includes a

term containing ETA as well. In this way slow routes are forced to save more fuel, whereas

the fuel saving of faster routes is not handled as strictly. The time deduction

∆

t is used to

obtain objective function values that still can be minimized, even if an arrival on schedule is

reached. The weight in equation (5.2) serves to provide values of equal magnitude for both

objective functions. In the terminology of evolutionism this modeling of objective functions

denotes an increasing selective pressure in favor of faster route designs. It is a common

technique to control the optimization process by modifying the fitness function instead of

applying the original objective. This method should not be mistaken for the application of

utility functions in multi-criteria decision-making. Utility functions are a combination of

objectives and assigned weights to emphasize single objectives. This is deliberately not

applied here, as it should be the conscious decision of a master which route he/she takes from

a set of feasible routes and which objective is considered to be more important than another.

Fig. 46: modeFRONTIER process flow chart

Pareto Optimum Ship Routing, Chapter 5

For the optimization task the generic optimization software modeFRONTIER of ES.TEC.O

s.r.l. (1999) is utilized. Fig. 46 shows the flow chart of modeFRONTIER as it is used for the

setup of the route optimization. Following the process flow, the diagram can be read from left

to right. Solid lines represent the data flow, dotted lines the process flow during the

optimization. With the nodes on the left side free variables are defined as input parameters.

These are transferred to a node that performs the generation and evaluation of a route variant,

called VirtualShip. The result is given to output parameter nodes and can now be passed over

to nodes that perform the evaluation in terms of constraints and objectives.

In the middle of the figure four other nodes can be seen. The DOE (design of experiment) is

needed to provide initial start designs for the optimization. The node below serves to select

the optimization method and settings like the number of loops or convergence parameter.

The two nodes below contain modules for the process control to decide whether a simulation

is successful or not.

Fig. 47: Feasible route designs and PARETO frontier

Fig. 47 shows a result of a multi-objective route optimization at rough sea conditions

(20. – 30.01.2002). Each dot represents a feasible route that can be taken into account for a

route decision considering ETA and FUEL. Infeasible solutions are not depicted.

The bold black circles represent the initially given designs, the set of feasible solutions to start

the optimization. This set has to be produced in advance of the optimization, for example by a

random search method.

A bold line, the PARETO frontier, borders the solution space. It is built by the set of solutions

for which a single objective cannot be further improved without deteriorating any other

objective. All feasible routes lie above and to the right of the PARETO frontier. Those

variants that are closest to the frontier display low passage time for reasonably low fuel

consumption. Apparently it is impossible to decrease ETA below certain limits without

impairing FUEL.

Pareto Optimum Ship Routing, Chapter 5

Fig. 48 shows the fastest variant of these routes (ETA 143h, FUEL 594t). Furthermore the

wave data of the deterministic forecast are depicted. The wave fields in the figure mostly

propagate from west to east. The contour lines represent the current significant wave height

during the crossing. The route and the actual position of the ship (a dot at the end of the bold

black line) are drawn, too.

Fig. 48: Optimized westbound North Atlantic crossing

Four stations of the journey are depicted:

• Upper left picture: 24h after departure a strong wave field passes the ship, it turns

towards North before reaching the ship.

• Upper right picture: A second wave field is surrounded northerly during the third day.

• Lower left picture: The protuberance of a third wave field towards Northwest at the

end of the fourth day causes extreme load to the ship and the crew, however,

constraints are not violated.

• Lower right picture: Finally the ship reaches the shallow and calm water in front of

New York.

Pareto Optimum Ship Routing, Chapter 5

5.2 Simplex versus Genetic Algorithm

For the choice of an appropriate optimization method various aspects are to be taken into

account:

• The optimizer should be able to handle more than one objective function. Although it

is commonly used to combine different objectives into one objective function, these

approaches generally produce intricate objective functions.

• The shape of the solution space and the one of the objective function is crucial for the

decision about applying a deterministic or a stochastic optimization method.

Deterministic methods are mostly sufficient for unimodal objective functions, whereas

multi-modal functions require methods that provide heuristic or stochastic attributes.

The introduction of constraints may change the characteristic of an objective function.

Within the presented ship routing, constraints distort a unimodal objective function to

one with local and global optima. Furthermore the spatial-temporal dependency,

between free variables describing the course and the velocity profile, causes a shift of

the range where free variables produce feasible solutions or even builds islands of

feasible combinations of free variables.

• Finally the optimizer should be able to produce the solution within reasonable time

expenditure.

The last aspect is not attended yet as it is the requirement of a final implementation. The

runtime improvements realized within the MATLAB routines that model the virtual ship are

mentioned in section 4.6.

In this section the focus is put on the qualitative aspects, in particular the one mentioned in the

second bullet point. For this purpose different optimization methods with enabled and

deactivated constraints are applied. The results for two representative investigations are given

below:

Fig. 49 shows four different optimization results. Two optimizations are conducted with a

SIMPLEX (representative for deterministic methods), and two apply a multi-objective genetic

algorithm MOGA (representative for stochastic methods). For SIMPLEX all investigated

route designs are depicted. The only requirement for SIMPLEX is to identify the minimum

ETA routes regardless of fuel consumption. For the MOGA optimizations only the PARETO

frontiers are shown. All feasible designs above and right from the PARETO frontier are not

depicted. The solid line is built by the PARETO optimum route designs where only the basic

constraints for main engine operability are enabled (the gray shaded ones in Tab. 7). The

dotted line marks the PARETO designs when also vertical accelerations are constrained. By

applying this constraint the minimum attainable ETA increases about 20h. Dots and triangles

mark the route designs that are investigated during corresponding SIMPLEX optimizations.

All optimizations start with relatively slow routes, i.e. at the right side of

Fig. 49. Identifying faster variants generally involves increasing fuel consumption.

Optimizations applying SIMPLEX perform 10 times faster than MOGA, but obviously they

are not able to find the assumed global optimum.

Pareto Optimum Ship Routing, Chapter 5

Fig. 49: MOGA and SIMPLEX optimization results

Fig. 50 is taken from investigations considering a separate variation of individual free

variables. In this case the two time minimum routes pointed out in Fig. 49 are considered.

These investigations serve to enlighten the differences in the optimization results of both

optimization methods:

Fig. 50: Variation of free variables

The figure depicts a variation of the variables s3 and v4, describing course and velocity profile

of a route segment at the middle of the journey. The left part of the figure is a variation of the

time minimum route found by MOGA, on the right hand side the one for the time minimum

route of the SIMPLEX, respectively. The s3-v4-combinations of the time minimum routes are

marked by a red circle for the MOGA and a blue square for the SIMPLEX optimization. The

contours represent ETA due to the variation of s3 and v4. Infeasible combinations of free

variables are shaded gray. Obviously both varied routes represent an optimum, i.e. the time

Pareto Optimum Ship Routing, Chapter 5

minimum route that is feasible regarding these two free variables. An improvement is only

possible by violating one or more constraints. Considering the SIMPLEX investigation at the

right side of Fig. 50, the combination of free variables for the time minimum route, found by

the MOGA optimization, produces an infeasible design. This means that variables are

dependent among each other, and being feasible or not depends on the combination of all

variables. Certainly Fig. 50 is a very malicious example for the characteristic of the objective

function and its limitation due to active constraints. But it nicely shows the extraordinary

dynamic behavior of the solution space that is after all responsible for the malfunction of the

SIMPLEX. It shows that the optimizer must be able to produce solutions that jump through

the solution space to overcome borders built by active constraints or to detect islands of

feasible combinations of variables. Furthermore this variation gives advice for necessary

investigations regarding the quality of an optimization result. An optimization result that

allows wider variations of free variable is more favorable than the one that immediately turns

infeasible. In this context also the severity of a constraint violation is of interest.

Investigations considering these aspects are given in the chapter 6.

Based on the presented considerations, it is decided to apply the MOGA as the primary

optimization method in the following. Certainly any other stochastic optimization method can

be adopted as well. Deterministic methods are applied only to control whether a local

refinement is possible.

5.3 Potential fuel savings by route optimization

Besides increasing the safety and reliability of ship operation and detecting the range of

attainable times of arrival ETA, fuel saving is a major target of route optimization. Two

different optimizations serve to assess the amount of possible fuel or time saving by applying

route optimization. The results are illustrated in Fig. 51 and Fig. 52:

• The first optimization uses unbound course and velocity profiles, i.e. the parameters

for the course and velocity perturbation si and vi are used as free variables. The

optimization serves to assess the maximum of possible fuel- or time savings. The

optimization result is represented by a PARETO frontier in Fig. 51 (solid line). For

this optimization achievable improvements become a trade-off between (i) fuel

savings by avoiding bad weather and (ii) additional fuel consumption because the

circumnavigation of bad weather areas causes longer routes. The achievable ETA and

the fuel consumption depend on both, the variation of the ship speed and the variation

of the sailed route.

• The second optimization is conducted with a fixed course, the unperturbed parent

route, i.e. only the parameters of the velocity perturbation vi are used as free variables.

This represents a strategy of ship operation where a fixed course is sailed and the

speed is reduced according to the weather conditions. The dashed line in Fig. 51

represents the resulting PARETO frontier. In this case ETA and the fuel consumption,

depend only on the perturbation of the velocity profile.

Pareto Optimum Ship Routing, Chapter 5

Fig. 52 depicts three routes, courses, velocity profiles, and the respective specific fuel oil

consumption, selected from the PARETO optimum routes pointed out in Fig. 51:

• route 1: fastest attainable route achieved by variation of both, course and ship speed,

• route 2: PARETO optimum route achieved by variation of both, course and

ship speed with the same ETA like route 3,

• route 3: fastest attainable route achieved only by variation the ship speed.

Fig. 51: PARETO frontiers for

optimizations with fixed and with unbound course

The left side of Fig. 52 shows the courses of routes 1, 2, and 3. The right side depicts the

respective ship speed and the specific fuel oil consumption during the journey for route

2 and 3, the routes with an ETA of 155h.

Fig. 52: Courses, velocity profiles, and specific fuel consumption for

optimizations with fixed and with unbound course

The numeric values for ETA, fuel consumption, and distance of the compared routes are given

in Tab. 8.

Pareto Optimum Ship Routing, Chapter 5

Tab. 8: Route comparison to assess the potential fuel saving

route ETA fuel consumption average ship speed distance free variables

no. [ h ] [ t ] [ kn ] [ nm ] [ - ]

1 143 594 19.7 2813 s and v

2 155 488 18.1 2810 s and v

3 155 567 17.6 2728 only v

It turns out that the spatial perturbation improves the optimization result. For the same time of

arrival (route 2 and 3) fuel savings of roughly 15% are possible. On the other hand it becomes

possible to arrive 12h earlier although the route takes a longer course (route 1).

Comparing the routes of Fig. 52 supports the opinion that an optimum route is a combination

of an optimum course and an optimum velocity profile. Due to the longer distance of route 2

the ship has to sail faster (see the averaged ship speed given in Tab. 8). During the first half of

the journey route 2 shows higher speed than route 3, however the specific fuel oil

consumption is comparable to the one of route 3. This is achieved by circumnavigating bad

weather and a resulting reduced added resistance due to waves. In the second half of the

journey route 3 has to speed up to meet ETA, consequently the fuel consumption increases.

For route 2 is admissible to maintain speed and the fuel consumption clearly falls below the

one of route 3.

Under rough weather conditions, here the optimization uses the forecast of a rough weather

condition in January 2002, the optimum course doesn’t necessarily follow the shortest track.

Furthermore, Fig. 48 and Fig. 52 illustrate the complexity of decision making in navigation. It

is evident that the benefit from route optimization increases in adverse weather situations. For

calm sea conditions the decision-making becomes easier. The additional resistance due to

waves diminishes and maximum fuel saving is achieved by sailing the shortest route as slow

as possible, i.e. below a certain sea state there is no need for route optimization. For example

at the medium weather conditions (January 2003) that are investigated within this study, the

attainable fuel saving reduces to 3%. In contrast to the optimization within the rough weather

scenario, this reduction of the fuel consumption is mainly achieved by reducing the sailed

distance. The influence of the weather condition on route optimization will be further

addressed in section 6.1.

Pareto Optimum Ship Routing, Chapter 6

6 Validations and extensions

The previous chapters present the setup of the route optimization procedure as a combination

of simulation and optimization techniques. Following, the reliability and capability of this

approach is analyzed. This is achieved by means of sensitivity analysis and

plausibility controls. General conditions, e.g. the applied weather forecast, are changed to

observe if the resulting optimization output is plausible and therefore reliable. On the other

hand, these studies serve to get an idea of the influence of various deciding parameters on a

ship route optimization result. Here parameter variations serve to assess the sensitivity of the

approach to changes of e.g. external conditions or parameters like threshold values and to

particular constraints. Thus the investigations also provide insight to the capability of the

optimization. In this way a first validation of the ship route optimization approach is achieved.

The following investigations will be presented in this chapter:

• Section 6.1 shows optimization results for a calm, a medium, and a rough weather

scenario. This investigation serves to get a first plausibility control of the optimization

results. As expected, the fastest route at calm sea conditions is the shortest course

sailed at maximum speed. Also at medium sea condition the fastest route takes the

shortest course, however, the ship speed has to be reduced temporarily. Within the

rough weather scenario the fastest route doesn’t follow the shortest course. Optimum

ship speed and optimum fuel consumption are achieved by circumnavigation of storm

areas. In doing so, a compromise between additional fuel consumption due to a longer

course and fuel saving by reducing the additional resistance due to waves is found. On

the one hand the ship speed is reduced due to operational requirements (posed

thresholds, permitted engine load), on the other hand it represents a tactical measure to

reduce the fuel consumption (await a storm area to pass by).

It turns out that the rough weather case is the most challenging case for a route

optimization. In this case an arrival on schedule seems impossible. During the

optimization numerous route variants are ranked as infeasible due to violated

constraints, and optimum routes show considerable deviations from the shortest

course. Therefore the following sections focus exclusively on the rough sea condition.

• Section 6.2 concentrates on the wave forecast and the modeling of swell by means of

wave spectra. The section focuses on the influence of different descriptions of long

crested waves on the optimization result, namely the spectra according to JONSWAP

and PIERSON-MOSKOWITZ.

Comparisons of 1D-spectra of the European Centre for Medium Range Weather

Forecast ECMWF and corresponding PIERSON-MOSKOWITZ spectra show

remarkable differences in the distribution of wave energy over the frequency. The area

under the spectral curves, i.e. the total amount of wave energy, is the same. During a

review of wave data it is observed that the JONSWAP and PIERSON-MOSKOWITZ

spectra represent so to speak edge cases of the 1D-ECMWF spectrum. In some cases

the energy distribution of the 1D-ECMWF spectrum resembles more the one of the

JONSWAP spectrum, sometimes it is more similar to the distribution of a PIERSON-

MOSKOWITZ spectrum. In the vicinity of storm areas the curve of the 1D-ECMWF

spectrum appears to lie in between the one of the JONSWAP and the one of the

PIERSON-MOSKOWITZ spectrum. Especially these areas are of interest for ship

routing as they are barely navigable.

The comparison of optimization results shows only small differences of the PARETO

frontiers for fast routes. At lower speed the frontiers are congruent. The number of

Pareto Optimum Ship Routing, Chapter 6

constraints for the slamming probability and engine overload that are violated

increases when a JONSWAP spectrum is applied. Comparisons of representative

transfer functions support this observation, both, spectra and transfer functions show

maxima in the same frequency range.

However, the optimization result is governed by another constraint, the one for

vertical acceleration on the bridge. This constraint shows no sensitivity regarding the

applied spectrum.

Based on these investigations it is decided to apply a PIERSON-MOSKOWITZ

spectrum in the following, as it is the standard spectrum for the North-Atlantic.

Certainly slamming and engine overload are not always of minor importance for a

route decision. For this reason the introduction of shape parameters, e.g. similar to the

parameters of the JONSWAP formulation (equation 3.5), can improve the quality of

the forecast data that are used within the route optimization. These parameters can be

provided as a supplement to the weather data, i.e. the forecast would exist of data for

the significant wave height, wave period, wave angle, and three to four shape

parameters.

• Section 6.3 focuses on the geometric description of the hull shape and its influence to

predicted ship responses and to optimization results. This is done to give an answer to

the open question on how accurate a particular hull shape has to be modeled for the

determination of ship responses within routing. For this purpose the seakeeping

characteristics of two further ships, scaled to Panmax-size and combined with the

calm water resistance and the engine characteristics of the HANNOVER EXPRESS,

are used for a route optimization in the rough weather scenario.

The optimizations show differences in the PARETO frontiers, in the courses, and in

particular in the velocity profiles of optimized routes. Advantages and disadvantages

of a particular hull shape result from the contributions of all considered seakeeping

criteria, i.e. optimality cannot be reduced to the contribution of a single criteria. Even

if the influence on the seakeeping behavior of the hull shape above the water line is

neglected by strip theory, it appears recommended to use a geometric description of

the ship hull that is as accurate as possible.

• Section 6.4 discusses the consequences of constraint- and threshold variations on the

optimization result. During the investigations so far, it is observed that the

optimization results for the rough sea scenario are strongly governed by constraints. In

particular the threshold for admissible vertical acceleration on the bridge shows a

strong influence on the shape of a PARETO frontier. The definition of reasonable

threshold values turns out to be one of the major tasks for the set-up of a route

optimization procedure. Here the variations serve to confirm or, if necessary, to adjust

the originally applied thresholds.

As expected, the variations of the threshold for the significant amplitude of vertical

acceleration on the bridge show a considerable influence on the optimization result. It

has to be decided if this threshold is set too tight. Comparisons with optimization

results applying a second evaluation method for vertical accelerations, namely the

motion sickness incidence MSI, advise against an increasing of the threshold for

significant amplitudes of vertical acceleration (currently 0.15g).

Regarding the threshold for the slamming probability, variations of threshold values

within a reasonable range have no influence on the PARETO frontier. The concept of

evaluating the slamming probability may need revision.

Pareto Optimum Ship Routing, Chapter 6

Following variations and combinations of thresholds for vertical acceleration and

slamming support the observations that:

- Admissible variations of the constraint regarding the slamming probability do

not affect the shape of the PARETO frontier, however, they influence the

optimization process.

- The constraint for vertical acceleration is dominant in any case.

- The fastest routes taken from optimizations with varied threshold values only

differ in 1h regarding the time of arrival and 3% regarding the fuel

consumption, although they show distinct differences regarding their courses,

i.e. the effect on the PARETO frontier is small. A cautious interpretation is:

Even if the routing-problem is characterized by high complexity, it seems to be

of “good nature”, i.e. deviations from an optimum course do not totally disturb

the merit of an optimized route.

• Section 6.5 discusses an extension of the optimization approach with a further

constraint, namely for parametric rolling.

Parametrically exited roll motion poses a serious threat to a safe operation of ships. In

particular ships with large bow flare and stern overhang, like container ships, are

affected. The motion sequence of parametric rolling is not comparable to the one of

resonant rolling. Therefore this phenomenon requires a separate treatment. Here the

model posed in section 4.9 is tested and discussed with regard to its influence on a

route optimization result.

The applied thresholds for the absolute wave height and the relative wave height at the

center of gravity, the fore-, and the aft perpendicular can be adjusted reasonably

among each other. The developed extension serves to avoid wave conditions that are

suspicious for parametric rolling, a prediction however is not possible. Regarding the

fact that differences in the weather development become even more probable the more

a forecast reaches into the future and that already relative small waves are able to

excite parametric rolling, it appears recommended to support the officer on watch by

means of a monitoring- and warning system.

• Section 6.6 addresses the robustness of an optimization result regarding the

uncertainty of a weather prediction. For this purpose a criterion representing the

robustness of an optimized route to probable weather changes is established. It is used

as constraint as well as objective function for the optimization.

A mean-ensemble forecast provides only information of restricted usability regarding

the robustness of a route. Whereas the ensemble forecast (the set of 50 forecasts with

an equal probability of occurrence) provides meaningful results concerning the

expected practicability of an optimized route. An optimization with simultaneous

consideration of three objectives, accounting for the time of arrival, the fuel

consumption, and for robustness, provides excellent results. Two extremely robust

PARETO-optimum routes that match the schedule are identified. They differ distinctly

from the time minimum route optimized within the deterministic forecast that is too

slow to match the schedule. A recalculation in analyzed weather points out that these

routes are only affected by minor and furthermore bearable violations of the constraint

for vertical acceleration. Therefore they appear to be feasible.

The results of this investigation emphasize the importance of reasonably appointed

thresholds. A constraint that is set too tight can exclude an optimization result that

probably turns out to be feasible and favorable later on, i.e. the range of permissible

solutions is narrowed unnecessarily.

Pareto Optimum Ship Routing, Chapter 6

From these optimizations follows a categorization of constraints:

- Hard constraints, e.g. engine overload, this constraint should be set tight to

keep a power reserve.

- Soft constraints, e.g. for vertical accelerations on the bridge; if minor

exceedings are permitted this constraint can be set less tight to avoid that

favorable route variants are discarded ab initio.

To compare single routes of different optimization results mostly time minimum routes are

used. Naturally all PARETO designs, and therefore also routes of different optimizations that

match the same ETA, can be used for comparisons. Actually a route that complies with the

schedule may be favorable to one that is faster but consumes more fuel. Nevertheless, for two

reasons it is decided to use time minimum routes for comparisons:

• In the applied rough weather case it turns out that an arrival on schedule is hardly

possible. In this case the time minimum route becomes interesting because it

represents the earliest possible arrival.

• From the optimization point of view, it is desirable to completely fill out the entire

solution space, i.e. to detect all feasible solutions. Generally optimizations start at slow

route designs and successively approach faster designs. In most cases it becomes more

difficult to build feasible routes the faster they get. Therefore, when addressing the

capability of an optimization set-up, it is meaningful to compare the time minimum

routes.

Pareto Optimum Ship Routing, Chapter 6

6.1 Various weather conditions

The first and simplest test case for a routing tool is a route optimization in a calm weather

scenario. In this case, when no constraints are active, the fastest ship route can only be the

shortest route at maximum speed. Furthermore it should be a PARETO design. Within the

following investigations three different weather conditions are used. Tab. 9 illustrates the

considered wave forecasts. A calm weather condition during summer 2001, medium and

rough sea conditions are taken from winter 2003 and 2002, respectively. Maximum occurring

significant wave heights and corresponding peak periods are mentioned to give an idea of the

weather conditions. These data refer to the first 5 days of the forecasts. The last four columns

of Tab. 9 depict the attained time minimum ETAs, the averaged ship speed during the journey,

the corresponding fuel consumption, and distances. These routes are found by a single

optimization applying a MOGA. At calm and medium wave conditions the fastest routes are

also the shortest ones. At medium sea condition the fastest route requires 3 hours more than at

calm sea condition. Here temporary speed reductions are necessary to avoid active constraints

and main engine overload. The fuel saving due to reduced speed and additional consumption

due to added resistance in waves almost balance. At rough sea condition the shortest route is

not the fastest anymore. Due to the wave conditions the ship speed has to be reduced

significantly and therefore an arrival on schedule is impossible (this requires an ETAmin of

approximately 130h). On the other hand, the fuel saving due to the reduced speed exceeds the

additional consumption due to the added resistance. Anyhow, this is no cause for delight, as

costs for delay and main engine maintenance arise.

Tab. 9: Applied forecasts and optimization results

sea

condition period H1/3 max

[m]

at Tp

[s]

ETA min

[h]

mean Vs

[kn]

FUEL

[t]

DIST

[nm]

calm 6.6. – 16.6.2001 5.4 10.2 114 23.8 694 2716

medium 1.1. – 11.1.2003 12.0 15.4 117 23.2 704 2716

rough 20.1. – 30.1.2002 15.0 14.9 143 19.7 585 2811

Fig. 53 depicts the optimization results for the three forecasts. The diagrams on top are fuel

consumption FUEL vs. estimated time of arrival ETA and distance DIST vs. ETA for the calm

weather condition. In the middle, the same diagrams are shown for the medium wave

condition. The diagrams for the rough sea condition are below.

Each blue dot represents a feasible route variant. The red dots mark feasible PARETO

optimum routes. Within the calm sea- and the medium sea case the optimization identifies the

shortest possible connection (2716nm) to be the fastest, and further more a PARETO

optimum route. Slower PARETO optimum routes do not necessarily take the shortest track.

Here the trade-off between (i) increasing fuel consumption due to additional resistance in

waves and a longer distance and (ii) fuel saving due to speed reduction is crucial. Taking

2728nm, the route at medium sea that matches schedule is still quite near to the shortest track.

Within the rough sea condition, no feasible routes are identified on the shortest track. This

does not mean that there aren’t any feasible routes, but they are not considered because they

are too time-consuming.

Pareto Optimum Ship Routing, Chapter 6

FUEL vs. ETA DIST vs. ETA

Fig. 53: Optimization results at different wave conditions

Fig. 54 compares the PARETO frontiers for the three weather conditions. It illustrates the

remarkable influence of the weather condition on the fuel consumption. At calm and medium

sea condition an arrival on schedule is possible, i.e. there are available routes of ETA = 130h.

But for the medium sea condition the fuel consumption increases about 15%.

Pareto Optimum Ship Routing, Chapter 6

Fig. 54: PARETO frontiers for calm, medium and rough sea

Fig. 55 shows the fastest routes for each weather condition. It clearly shows that in this case

the fastest route at rough sea condition is not the shortest one. Altogether the optimization

provided the expected result: as far as no constraints are active the shortest route is also the

fastest one. But the shortest route is not always necessarily a PARETO optimum route.

Fig. 55: Fastest routes for different wave conditions

Fig. 55 also illustrates a further characteristic of the optimization applying a multi objective

genetic algorithm MOGA. Between Newfoundland and New York the fastest rough sea route

doesn’t follow the shortest track, as both other routes do. But there are no extreme wave

conditions that prohibit taking the shortest way, neither in the rough weather case.

Furthermore it is no restriction of the applied B-spline modeling technique, like a southern

overshoot due to the northern perturbation in advance. This disfigurement is caused by the

applied optimization method, the MOGA. Taking the shortest track between Newfoundland

and New York yields an improvement of 0.5 hours in ETA, i.e. the attainable improvement is

relatively small. Applying elitism during the optimization ensures that the fastest variants are

taken over to the next generation. Nevertheless, the probability to be considered for

recombination, mutation, and crossover depends on the fitness relative to other members of

the generation. Therefore, when having flat shaped optima like here, an improvement of the

optimization result becomes less probable the more the optimum is reached.

Pareto Optimum Ship Routing, Chapter 6

It has to be kept in mind, that the applied MOGA is not at all adapted to the route

optimization problem. Fore sure this is an important task for further research.

The optimizations at various weather conditions show that the rough weather case is the most

challenging case for a route optimization. For this reason all following investigations focus

exclusively on the rough sea condition.

6.2 Wave spectra variation

Considering the applied wave data of ECMWF, peak periods of 13 – 15s (ω = 0.4 – 0.5 rad/s)

are typical in the vicinity of rough weather in the North Atlantic. Therefore the wave spectra

represented in Fig. 56 are typical for a sea state that can be found in rough weather areas. The

figure shows an ECMWF-1D spectrum compared to a JONSWAP and a PIERSON-

MOSKOWITZ spectrum built by the appropriate wave parameter H1/3 and Tp. Obviously the

spectral distribution of wave energy strongly depends on the applied spectra. With a peak

periods Tp of 13.5s, they show a mean wave length of about 300m that is approximately the

length of the considered ship. It is assumed that the application of different spectra causes

differences in the predicted ship motion. Consequently a route optimization result should be

affected as well, especially when it is governed by constraints.

Fig. 56: Typical storm spectra for the North Atlantic,

PIERSON-MOSKOWITZ, JONSWAP and ECMWF-1D spectrum

To analyze this phenomenon, two comparative optimizations are conducted. Both use the

deterministic rough-weather forecast. In one case the wave parameter are associated with a

JONSWAP spectrum, in the other case a PIERSON-MOSKOWITZ spectrum is applied. The

application of the rough weather forecast ensures that the optimization result is governed by

constraints. Fig. 57 depicts the PARETO designs for these optimizations. In the range of ETA

from 148 – 159 h no remarkable difference in the PARETO frontiers is found. However,

below 147h the frontiers separate and the application of the PIERSON-MOSKOWITZ

spectrum provides more feasible designs down to 143h. To prove that this is not just a defect

of the optimization, all designs evaluated during the MOGA optimization applying the

PIERSON-MOSKOWITZ spectrum are recalculated applying a JONSWAP spectrum.

Pareto Optimum Ship Routing, Chapter 6

The recalculation approves the former optimizations as shown in Fig. 5714.

Fig. 57: Optimization result applying a PIERSON-MOSKOWITZ

and JONSWAP spectrum

To investigate the reason for the differences in the optimization results, the focus is put on the

active constraints. Tab. 10 gives maximum, minimum, and mean values for significant

amplitudes of vertical accelerations on the bridge and for slamming probability regarding the

optimizations shown in Fig. 57. Below the number of routes, thereof infeasible routes and the

number of active constraints, absolute and as percentage of all routes are given. C_zacc for

example denotes the constraint for vertical acceleration on the bridge. Remarkable differences

are found for the constraints that control the maximum acceptable slamming probability

C_slprob and the specific fuel oil consumption C_spef; the latter becomes active in case of

main engine overload.

Tab. 10: Comparison of JONSWAP and PIERSON-MOSKOWITZ spectra

JONSWAP PIERSON - MOSKOWITZ

min mean max min mean max

zacc [m/s2] 0.09 0.16 0.27 0.09 0.16 0.25

slprob [%] 0 0.44 7.8 0 0.2 5.3

routes 18964 (100%) 18925 (100%)

infeasible 11483 (61%) 10682 (56%)

active C_zacc 10409 (55%) 10361 (55%)

active C_slprob 529 (3%) 229 (1%)

number

active C_spef 3725 (20%) 1544 (8%)

To verify this evaluation 33000 randomly produced designs presented in Tab. 11 are

investigated (for a sufficient coverage of the design space the variables RTD and v9 are set

constant. 15 free variables are remaining: 215 = 32768 Æ 33000 designs. This means that for

each free variable, there are at least two design variants that differ in only one variable).

14 Recalculated designs are called DOE, design of experiment. This terminology is borrowed from the generic

optimization tool modeFRONTIER. Here a simple evaluation of designs without using any optimization

technique is termed DOE.

Pareto Optimum Ship Routing, Chapter 6

Tab. 11: Recalculation comparing JONSWAP and PIERSON-MOSKOWITZ spectra

JONSWAP PIERSON - MOSKOWITZ

number of designs 33000 (100%)

thereof infeasible 30861 (94%) 28297 (86%)

active C_zacc 29090 (88%) 25877 (78%)

active C_slprob 7538 (23%) 1326 (4%)

active C_spef 20790 (63%) 14077 (43%)

It shows that JONSWAP spectrum produces 8% more infeasible routes. Here even an increase

of infeasible routes due to bridge acceleration is found. But while C_zacc is active on routes

of the whole solution space, the other two constraints affect mostly fast and fuel consuming

routes. It is obvious that fast routes are more vulnerable to main engine overload. Due to the

ship speed the calm water resistance is relatively high and only a limited sea margin is

possible. Further on, regarding active constraints, Fig. 58 and Fig. 59 show representative

transfer functions and response functions for the added resistance due to waves for a wave

encounter angle of 180°.

Fig. 58: Representative response functions for different ship speeds

(enlarged diagram of Fig. 18)

Fig. 59: Representative transfer functions for the vertical motion on the bridge and for

the relative motion between water surface and bow (enlarged diagrams of Fig. 18)

Pareto Optimum Ship Routing, Chapter 6

The response functions for the added resistance show a pronounced peak in the frequency

range where maximum wave energy can be found in a typical rough weather scenario, visible

in Fig. 56. The magnitude of the peak is speed dependent but shows nearly no sensitivity

within a variation of the encounter angle of up to 50° from head sea. Just the peak shifts

slightly due to the shift in the encounter frequency. Consequently the increase of the added

resistance due to waves at high speeds makes a main engine overload even more probable.

Regarding the transfer functions for the relative motion at the bow and accordingly the

slamming probability, here also a clear speed dependency can be found. For a wave encounter

angle of 180°, these transfer functions show a pronounced peak that matches the peak of the

rough weather spectra. This suggests that at rough seas in particular fast routes are affected.

The transfer functions for the vertical motion on the bridge rather show an inverse speed

dependency, i.e. the slower ship is affected more. Furthermore there are no pronounced peaks.

The humps at

= 0.6, corresponding to a wave length of approximately 170m, are caused by

the pitch motion characteristic of the vessel. This conclusion regarding the motions also holds

for vertical accelerations. Finally the shapes of the transfer functions support the observation

derived from the optimization results and from the recalculation.

Considering Fig. 57 the PARETO frontier applying the JONSWAP spectrum ends with a

steep increase in fuel consumption for the fastest routes. This could be validated by local

optimizations in this area. Thus the solution space offers the opportunity to take routes of

longer distances at higher speed and accordingly higher fuel consumption, whereas the

improvement in ETA is negligible (fastest design: ETA = 144.6h FUEL = 623t DIST =

2820nm, followed by: ETA = 144.9h FUEL = 574t DIST = 2812nm). Because of the

considerable fuel savings, a route decision would take the second PARETO design into

account.

To conclude, two time minimum routes are compared in Fig. 60. The one of the optimization

applying PIERSON-MOSKOWITZ spectrum takes 143h, the one for the JONSWAP

spectrum 145h respectively.

Fig. 60: Time minimum routes for the PIERSON-MOSKOWITZ and the JONSWAP

spectrum

Both optimizations found more or less the same geographical shape for the time minimum

route. The velocity profiles differ in particular at the beginning of the journey. Around

midway these differences are balanced and finally the difference in ETA becomes relatively

small. At least the basic trend of the velocity profiles is the same. Regarding the

Pareto Optimum Ship Routing, Chapter 6

investigations represented in the following sections, it can be said that taking a

PIERSON-MOSKOWITZ spectrum instead of the more accurate 1-D spectra causes a

negligible error. However with a view to reduce all involved sources of defect, a more

accurate representation of the sea condition is desirable. For the time being the application of

a PIERSON-MOSKOWITZ spectrum is absolutely sufficient to investigate this novel set-up

of a routing approach.

6.3 Modified hull shape

A frequent question in connection with route optimization regards the required knowledge of

the particular hull form and load condition. At least draft and trim are deciding for the

determination of transfer functions, whereas the metacentric height GMT and the radii of

inertia for pitch and roll can be changed in remarkable magnitude, SEAROUTES (2003).

Here the focus is put on the hull shape. To accomplish a hull form variation, two container

ships with a hull topology similar to the one of HANNOVER EXPRESS are taken from the

SEAWAY/OCTOPUS hull form database. They are called JOURNÈE.044, also known as

S175, and VERSLUIS.051. Both ships are scaled to the main dimensions of CMS

HANNOVER EXPRESS. To solely regard the influence of the hull form variation due to

waves, the calm water resistance and propulsion characteristics of the HANNOVER

EXPRESS are applied to both scaled ships. However, it is not the aim to conduct a systematic

variation on hull form parameters that are crucial for the seakeeping behavior of a ship, but:

• to give a rough assessment on the desirable accuracy of a hull description,

• and as the applied seakeeping method for the determination of transfer functions is

based on strip theory, the hull form variations also provide an insight into the

sensitivity of strip theory.

Tab. 12 shows the hull properties of CMS HANNOVER EXPRESS, scaled VERSLUIS.051,

and JOURNÈE.044. The lines plans available within the Seaway hull form editor are given in

Fig. 61. The differences in the frame spacing between HANNOVER EXPRESS and both

other ships visible in the bow and stern region in Fig. 61 do not affect the numerical result.

From the viewpoint of main dimensions they are similar ships. Nevertheless, there are distinct

differences in the displacement, the centers of buoyancy, and the centers of the water-plane

area. For VERSLUIS.051 and HANNOVER EXPRESS also block- and water-plane area

coefficient are similar or actually the same. In this way it is possible to compare three similar

ships.

Tab. 12: Hull properties of compared ships

hull parameter Versluis.051

scaled

CMS

Hannover Express

Journée.044 (S175)

scaled

LWL / B / T

(trim by stern) [m]

293.4 / 32.25 / 11.79

(0.14)

293.4 / 32.25 / 11.79

(0.14)

293.4 / 32.25 / 11.79

(0.14)

displacement [t] 77287 74474 62609

CB (LWL) 0.68 0.65 0.55

LCB to AP [m] 145.1 136.2 140.4

CWP (LWL) 0.80 0.80 0.67

LCF to AP [m] 140.8 129.1 134.4

Pareto Optimum Ship Routing, Chapter 6

VERSLUIS.051 HANNOVER EXPRESS JOURNÉE.044 (S175)

Fig. 61: Lines plans of VERSLUIS.051, CMS HANNOVER EXPRESS and

JOURNÉE.044 (S175)

Fig. 62 shows diagrams for the added resistance due to heading waves at a ship speed of

23kn. Compared from left to right, only a slightly increasing resistance can be observed.

Regarding the calm water resistance of 1900kN and a wave height of e.g. 6m at a peak period

of 15s the total resistance would amount about 2400kN. Therefore the deviation of ±50kN

would represent 2%, a small but remarkable amount.

Fig. 62: Added resistances due to waves for different hull shapes

Fig. 63 shows the slamming probability in head waves at 23kn. It is the combination of all

mentioned hull form parameter that impact on the heave and pitch motion and on their phase

shift and consequently on the slamming probability. Therefore it is not possible to point out

the deciding parameter here, but it can be seen that the probability of slamming is quite

different for these three ships.

Fig. 63: Slamming probabilities for different hull shapes

Pareto Optimum Ship Routing, Chapter 6

Fig. 64 depicts the vertical acceleration on the bridge in head waves for a ship speed of 23kn.

Here as well remarkable differences can be observed. S175 differs most distinctly from the

others. For VERSLUIS.051 and HANNOVER EXPRESS the ship responses are more similar

in smaller waves, but they become more different when wave heights increase.

VERSLUIS.051 HANNOVER EXPRESS JOURNÉE.044 (S175)

Fig. 64: Significant amplitudes of vertical accelerations on the bridge

for different hull shapes

Fig. 65 shows optimization results and courses for the time minimum routes at rough sea

conditions. On the left side, PARETO frontiers for ETA and fuel consumption, on the right

side, the courses for time minimum routes are given. The differences of PARETO frontiers

and courses are appreciable, compared to the influence of the weather situations, presented in

section 6.1, they are quite small.

Fig. 65: PARETO frontiers and ETAmin routes for HANNOVER EXPRESS,

scaled VERSLUIS.051 and JOURNÈE.044 (S175)

At least the velocity profiles of the ETAmin routes, presented in Fig. 66, show that there are

clear differences between these optimum routes. In particular at the beginning of the journey

CMS HANNOVER EXPRESS is able to sail faster than both other ships. In this period of the

journey the ships already have to withstand a first storm region. Mainly due to the constraint

posed for main engine load and vertical acceleration on the bridge, the ships have to reduce

velocity. The ship with the lowest acceleration response, i.e. HANNOVER EXPRESS, is able

to maintain the highest speed, and although the ship takes the longest route it arrives earlier

than the others.

Pareto Optimum Ship Routing, Chapter 6

Fig. 66: Velocity profiles for ETAmin routes of HANNOVER EXPRESS,

scaled VERSLUIS.051 and S175

Balancing the factors that are deciding for the accurateness of an optimization result suggests:

• There is no need to model the hull shape of the ship in waves with extreme accuracy

while underlying conditions do not serve to maintain this accuracy. For example,

balancing the effect of probable changes of the predicted weather development on the

PARETO frontier (this topic will be addressed in section 6.6) and the effects of hull

form variations on the seakeeping behavior (shown here) the ship behavior of a related

ship may be sufficient for the modeling within route optimization. Furthermore, by

applying strip theory, the hull shape above the calm water line is kept out of

consideration although it affects a lot of the seakeeping ability of a vessel. Therefore

neglecting the hull above the waterline may deteriorate improvements by regarding

the individual hull shape and load condition.

• On the other hand, even if the differences in the PARETO frontiers and courses are

relatively small, the differences of the proposed velocity profiles are remarkable. From

this point of view, a proper hull description is recommended and at least a closely

related hull should be used. The recommendation should be to minimize the over-all

error by minimizing the partial error and therefore to use an accurate hull description

whenever it is possible. In the case of rough weather conditions this becomes even

more important as optimization results are governed by constraints. In this regard

seakeeping prediction based on strip theory turns out to be a suited method. For the

considered hull form variations it is sensitive enough to take an influence of the hull

shape on the optimization result into account. It is decided to apply strip theory, as

there is a need for a fast method for the assessment of ship responses in waves. Known

disadvantages like in the prediction of roll motion or neglecting the upper hull have to

be taken into account by empirical and semi-empirical methods, as far as possible.

Pareto Optimum Ship Routing, Chapter 6

6.4 Constraint and threshold variations

Three further aspects are identified during the investigations focusing the optimization at

rough sea conditions (cp. Tab. 9 on page 71):

• The extreme number of infeasible routes due to active constraints for the vertical

acceleration on the bridge indicates that this constraint may be set too harsh.

• The constraint for the slamming probability is only active in cases of machine

overload and therefore it probably can be left out of consideration at least for the full

load or the design load case.

• The constraints for lateral acceleration C_xacc and C_yacc do not influence the

optimization at all. Independent of any admissible threshold variations, they only

became active when C_zacc distinctly exceeded its threshold. For this reason these

constraints are left out in the following.

To step forward in the matters addressed by the first two bullet points, the threshold values for

the maximum slamming probability and maximum vertical acceleration on the bridge are

varied separately. The influence on the optimization result is observed. If advisable,

modifications to the constraints are realized and investigated afterwards.

Regarding the last bullet point, roll motions substantially contribute to the transverse

acceleration on the bridge. Furthermore roll motion is surely an important parameter in the

operation of ships and in particular in the operation of container ships. As explained in

chapter 4 resonant rolling is not included to this first approach yet. The fact that in the

following neither roll motion nor transverse acceleration is further regarded does not imply

that these parameters are of subordinate importance. They can be the deciding parameter for

other routes, other load cases, or other weather conditions. As regards parametric rolling, it is

known that strip theory is not able to account for it. An approach to overcome this and at least

to identify situations that are suspicious for parametric rolling is presented in section 4.9, first

results are given in section 6.5

6.4.1 Variation of the threshold for the vertical acceleration on the bridge

Fig. 67 shows PARETO designs for an investigation considering the constraint for vertical

acceleration on the bridge, C_zacc. The threshold is varied from 0.10 to 0.20g. Expectedly a

tightening of this constraint increases the attainable ETA. For a particular time of arrival, it

produces higher fuel consumption mainly caused by longer distances. Whereas the PARETO

frontiers for 0.10, 0.13, and 0.15g (0.17g is very close to 0.15g and therefore left out in the

diagram) are relatively close to each other, the frontier for 0.20g clearly separates. Here ETA

can be reduced significantly but this necessitates more fuel due to higher ship speed. In view

of operational costs, for fuel consumption comparable to routes at tightened constraints, an

earlier arrival is possible. Actually there is only a small difference between the PARETO

frontier for a threshold value of 0.20g and the frontier when this constraint is deactivated. In

the second case the optimization result is completely governed by another constraint, namely

the one for machine overload.

Pareto Optimum Ship Routing, Chapter 6

Fig. 67: Variation of the threshold for

vertical acceleration on the bridge C_zacc

Fig. 68 illustrates the main reason for the reaction of the optimization result to the variation of

C_zacc. It is the reduced distance and the increase of the ship speed, made possible by

relaxing the constraint. The shift of the 0.20g-PARETO frontier towards faster route designs

at approximately ETA = 145h, visible in Fig. 67, results from finding a shorter, still feasible,

and faster route during the optimization.

Fig. 68: Time minimum routes for varied C_zacc

At the end of the journey depicted in Fig. 68 the same characteristic can be found as already

addressed in section 6.1. The routes do not follow the shortest track. This is the mentioned

disadvantage of the genetic algorithm.

However, another rather strange result of these investigations should be mentioned. It is

illustrated in Fig. 69.

Pareto Optimum Ship Routing, Chapter 6

Fig. 69: Optimization result for deactivated and enabled C_zacc

Two optimization results are shown. One has got a threshold for C_zacc of 0.2g and the other

one a deactivated constraint for vertical acceleration. The optimization with enabled

constraint shows PARETO designs at decreased FUEL and ETA (left part of Fig. 69).

Naturally these remain feasible when C_zacc is deactivated, i.e. enabling the constraint

improves the optimization result. On the right hand side of the figure the time minimum

routes of both optimizations are shown. The routes differ significantly. Within a range of ETA

from 123 to 145h all PARETO results are more or less variations of the respective time

minimum routes shown in Fig. 69. At ETA = 145h both optimizations identify PARETO

designs of at least the same distance. These observations suggest that the unconstrained

optimization identifies and finally sticks to a local optimum. Assuming that a stochastic

method is always able to find the global maximum (what is obviously wrong) the

optimization has to be ranked as incomplete. Furthermore, considering Fig. 50, it is

reasonable to ask for the quality of an optimization result. A local optimum with high

tolerance to variations of free variable and consequently a great potential to react on weather

changing, may be preferred to a global optimum that doesn’t allow any changes. On the other

hand, the shape of the ship response, as presented in Fig. 28, qualifies this position. The flat

characteristic of this response function denotes that changes in H1/3 or Tp cause relative small

changes in the ship response. For this reason the violation of a constraint may be small as well

and tolerable. This matter is certainly an important point within route optimization and should

be addressed in continuative research.

6.4.2 Variation of the threshold for the slamming probability

The second threshold variation considers the constraint for the slamming probability

C_slprob. During the investigations at rough sea conditions, it turned out that the threshold

for the tolerable slamming probability has to be reduced from 3% to less than 0.5% to get a

remarkable influence of this constraint on the optimization result. Otherwise the constraint

posed for the operability of the main engine C_spef becomes active long before C_slprob

approaches its threshold value. But without doubt, Panmax container vessels are able to slam

in such a way that this causes damage to the ship and to the cargo.

Pareto Optimum Ship Routing, Chapter 6

Various aspects may be considered to step forward in this matter:

• The assumption of a statistical independence of bow emergence and velocity of the re-

entry in the calculation of the slamming probability (equation 4.25) does not hold.

Therefore the probability is underestimated.

• The influence of horizontal velocity components is neglected, although they contribute

to the local pressure impact.

• The applied methods are not suited or may need a revision to adapt them to modern

ship hull forms and ship constructions.

The first two bullet points address themes that are beyond the scope of this study, however the

latter one should be addressed as far as applied methods are affected.

The applied definition of slamming follows OCHI (1964) and is conducted according to the

recommendations of JOURNÉE (2000). OCHI investigated ships with a TAYLOR15 bulb.

That is why it may be reasonable to shift the slamming point an amount towards the forward

perpendicular FP when ships with a pronounced bulbous bow are considered. Here the bulb

reaches out 2%LPP in front of FP. Therefore the slamming point may be put at 8% behind FP

instead of 10%. Furthermore, according to OCHI the critical velocity is 4.88m/s. JOURNÉE

represents a second method considering a critical velocity calculation based on a critical

pressure. This method is set up by CONOLLY (1974) and is based on experiments with cones

and wedges. With a deadrise angle of approximately 25° at 10% behind FP for the vessel

considered here, the pressure coefficient within this calculation amounts to Cp = 10.

Following JOURNÉE the critical pressure is given by:

15005.0 mkNLgp PPcr =⋅⋅⋅=

, (6.1)

this results to a critical velocity of:

cr 4.5

⋅

&. (6.2)

These considerations show that it is possible to vary the velocity threshold and the location of

the slamming point in a reasonable range. Fig. 70 illustrates SEAWAY results for a variation

of the slamming point and the threshold for the relative velocity, i.e. a variation of the two

characteristic parameter used by OCHI where one is varied, whereas the other one is held

constant.

Fig. 70: Parameter variations influencing the slamming probability

15 David Watson Taylor, naval architect and engineer, rear admiral of the US Navy, 1864-1940, USA

Pareto Optimum Ship Routing, Chapter 6

The slamming probability for three ship velocities is calculated at a wave encounter angle of

180° and for H1/3 = 10m and Tp = 15s, a very extreme sea condition where the wave length

approximately matches the ship length. Although slamming should be quite probable in such

extreme waves, the calculated slamming probabilities do not exceed the 3% level, neither at

reduced velocity threshold nor by shifting the slamming point towards the bow. However, a

clear dependency on the addressed parameters is obvious and the re-consideration of these

parameters and of the threshold for the slamming probability may be the most available way

to come to a contemporary and appropriate method to define slamming.

The considerations above regard the CMS HANNOVER EXPRESS at design draft of 12m.

Reducing the draft to ballast load condition of 9m the slamming probability would increase to

6.2%, i.e. 30 slams per hour (at H1/3 = 10m, Tp = 15s, 23kn ship velocity, 180° wave

encounter angle, and no trim).

In this context “no active constraint for slamming” is also a result of the optimization, this

constraint may become active for other load conditions. Nevertheless, a reduction of the

threshold value for the slamming probability seems to be recommended. This, inter alia, will

be addressed in the following section.

6.4.3 Comparison to thresholds posed by NORDFORSK

So far the threshold values posed within SEAROUTES are applied here. These threshold

values refer to maximum tolerable significant amplitudes of accelerations or probabilities of

occurrence in case of the slamming probability. Other and more classified thresholds are

established during the Nordic Cooperative Project, NORDFORSK (1987). Tab. 13 and

Tab. 14 show threshold values posed within this project. Herein the threshold values are given

as root mean square values, RMS.

Tab. 13: NORDFORSK, general operability criteria (JOURNÉE, 2001)

Following equation (4.44) the relation between RMS and the significant amplitude of an

irregular motion follows from:

RMStda

a⋅=⋅=⋅= ∫2

1. (6.3)

A detailed explanation regarding the conversion of significant and root mean square values is

given in Appendix 3.

Pareto Optimum Ship Routing, Chapter 6

With regard to the slamming probability, JOURNÉE (2001) recommends a linear

interpolation of the threshold values according to the ship length. This suggests a reduction of

the threshold for C_slprob to 1%. Even if this does not affect the PARETO frontier of the

rough weather optimization, it may affect the optimization procedure.

Tab. 14: NORDFORSK, operability criteria for various types of work (JOURNÉE, 2001)

Regarding the threshold for the vertical acceleration on the bridge, the threshold for

intellectual work as given in Tab. 14 should be applied to ensure situation awareness in

dangerous situations. This would result in a threshold of 0.2g for maximum allowable

significant amplitude of vertical acceleration, i.e. a considerable lowering of the constraint.

To summarize the above made considerations:

• The threshold value for the slamming probability may be tightened to 1%.

• Additional the slamming point can be shifted to 8% behind LPP.

• Following Tab. 14 the threshold for the significant amplitude of vertical acceleration

may be put up to 0.2g. According to Tab. 13 it is even 0.3g.

Results of the optimizations with adjusted constraints and thresholds are presented in section

6.4.5. But before, a further method considering vertical acceleration, the motion sickness

incidence MSI, should be discussed. This is to reassure that a lowering of the constraint for

acceleration on the bridge is admissible.

6.4.4 Comparison of the motion sickness incidence MSI and significant

values for vertical acceleration

According to Tab. 13 and Tab. 14, a RMS of 0.1g to 0.15g of vertical acceleration on the

bridge seems acceptable. This would result in 0.2g to 0.3g of significant amplitude of vertical

acceleration. Thus the constraint for vertical acceleration C_zacc applying a threshold of

0.15g seems to be set too harsh. On the other hand, the duration of the exposure to the motion

and the frequency of excitation, both important parameters to analyze kinetosis, are neglected.

Therefore the MSI producing an optimization result comparable to a threshold of 0.15g for

C_zacc is identified. By getting a further indicator for the load on the crew, it should be

possible to get a comprehensive view on this matter.

Pareto Optimum Ship Routing, Chapter 6

Fig. 71: PARETO frontiers applying different thresholds for the vertical acceleration

Fig. 71 shows both, the PARETO frontiers of Fig. 67 produced by optimizations applying

different thresholds for significant amplitudes of vertical acceleration and the results for

various thresholds applying MSI. Posing 0.2g for maximum significant vertical acceleration

would require an admissible MSI of approximately 50%, an extremely high value. But also

the MSI related to 0.15g lying at about 40% MSI is still enormous. It has to be said that the

applied concept to calculate MSI does not integrate the exposure to accelerations of

consecutive time steps of a ship route. That is why the load on the crew is rather

underestimated by the MSI given here. In this regard, a lowering of the constraint for the

vertical acceleration as suggested by the NORDFORSK data is refused. So far it seems

recommended to keep the threshold for vertical acceleration on the bridge at 0.15g for

significant amplitudes or at an MSI of 40% respectively.

6.4.5 Optimizations with modified constraints

Closing this matter, the following four optimizations with modified constraints will be

compared, i.e. moderate variations of the constraints for the slamming probability and the

acceleration on the bridge as applied so far. This investigation serves to gain insight to the

functioning and to the interaction of applied constraints and the optimization method. The

threshold value combinations and resulting PARETO frontiers, given in Fig. 72, are:

• slprob of 1% at 8% LPP behind FP and a MSI= 40% Æ Curve 1,

• slprob of 1% at 8% LPP behind FP and a MSI= 30% Æ Curve 2,

• slprob of 1% at 8% LPP behind FP and C_zacc = 0.15g Æ Curve 3,

• slprob of 3% at 10% LPP behind FP and C_zacc = 0.15g Æ Curve 4.

To enhance the overview, the PARETO designs of Fig. 72 are connected by lines and not

depicted as dots. Dashed lines mark the position of PARETO designs applying MSI, solid

ones represent the acceleration threshold counterpart.

Pareto Optimum Ship Routing, Chapter 6

Curve 4 depicts the PARETO frontier applying the previous constraints (zacc = 0.15g,

slprob = 3%). It is enveloped by the frontiers applying 30% and 40% MSI. Curve 3 represents

the PARETO frontier preserving the acceleration threshold of 0.15g but reducing the one for

the slamming probability to 1% at 8% behind LPP. Even if differences in the PARETO

frontiers are relatively small, less than 1% ETA and 3% FUEL for the fastest designs, there

are differences caused by the interaction of the applied constraints as well as by the applied

optimization method.

Fig. 72: Optimization results for modified constraints

Fig. 73 illustrates the achieved distances of the PARETO designs represented in Fig. 72. It

can be seen that referring to the distance, in the low ETA region the PARETO designs clearly

separate.

Fig. 73: Distance over ETA for modified constraints

Pareto Optimum Ship Routing, Chapter 6

Fig. 74 illustrates the courses of the fastest route from each of the optimizations, the ETAmin

routes of curve 1 – 4.

Fig. 74: Fastest routes for modified constraints

Four conclusions are supported by these optimizations:

• In no case, the PARETO frontier is directly affected by the constraint for the

slamming probability, i.e. filtering the designs following to the optimization by

enabling or disabling this constraint does not change the shape of the PARETO

frontier or changes the attribute of a PARETO member from feasible to infeasible. But

in contrast to the former threshold (3% at 10% LPP behind FP), the optimization

applying the tightened threshold (1% at 8% LPP behind FP) produces about 5%

infeasible designs that are infeasible due to this constraint exclusively. This shows that

a reasonable variation of this constraint as discussed above makes it affecting at least

the optimization procedure. Therefore it seems recommended to maintain this

constraint for a route optimization, anyhow.

• Comparing the optimization results for the threshold of 0.15g maximum significant

amplitude of vertical acceleration (curves 1 + 2), the PARETO frontiers do not clearly

separate. Nevertheless, even if PARETO frontiers are quite close to each other, the

shape of the fastest routes clearly differs, as can be seen in Fig. 74. Furthermore Fig.

73 shows that these routes take tracks of a bit more distance, compared to the routes

where MSI is applied. Cross calculations, i.e. evaluating the optimization result of one

optimization with the constraints of another optimization, suggest that this again may

be the weak point of the MOGA mentioned above. The designs of both PARETO

frontiers remain feasible when replacing the appropriate slamming constraint. Either

the MOGA sticks to a local optimum or it doesn’t proceed due to the insignificance of

the achievable improvement. However, a further adaptation of the MOGA or any other

suited multi objective optimization algorithm is necessary.

• In the range of ETA = 147h - 155h the PARETO frontiers for 40% MSI and 0.15g zacc

(cures 1 + 3) mostly overlap. At low speed, above ETA = 155h the frontiers separate.

This is caused by different start designs in the first generation of the optimizations. By

contrast, the separation of the frontiers below ETA = 147h is caused by the different

constraints referring to vertical acceleration. In this region the PARETO designs

Pareto Optimum Ship Routing, Chapter 6

achieving a threshold of 0.15g in significant amplitude of vertical acceleration show a

MSI above 40%. Finally both methods, applying thresholds for MSI or for significant

amplitudes of vertical acceleration, are useful to pose a constraint. It is a matter of

taste to put the focus on the duration of an exposure or on maximum values or on both.

• All routes in Fig. 74 differ only in the first half of the journey. Considering Fig. 48 on

page 61 shows that it is the circumnavigation of the first storm (upper right part of Fig.

48). Here the differences in the constraint for the vertical acceleration cause the

differences of these route designs. Besides differences in the velocity profiles, which

are not shown here, a remarkable lateral displacement is obvious. Although the routes

show a lateral shift of approximately 3° (180nm = 333km) their difference in ETA is

only 1h. Compared to the medium- and calm weather optimizations the difference in

fuel consumption of 3% is relatively small, too. So far it seems that sizeable variations

of the constraints influence, but not dramatically distort, an optimization result. In

other words, there may be a considerable range of admissible route variations without

deteriorating the qualities marked by the PARETO designs. This characteristic may

rapidly change when further constraints are applied and of course, these assumptions

deserve further study.

6.5 Parametric rolling

Parametric rolling is one of the main reasons for cargo losses in seaborne container

transportation. According to the importance of this topic it is worth to investigate how often

sea states that are suspicious for parametric rolling occur during the route optimization. Since

most mathematical models, e.g. strip theory, account only for a coupling of surge, pitch and

heave or sway, roll and yaw respectively, they are not able to predict this motion sequence.

Nevertheless, the module established in section 4.9 may serve to answer the question, if

parametric rolling should be included to route optimization or not. Even if it is not able to

predict, it may help to avoid this kind of extreme roll motions.

To investigate in this matter the constraints developed in section 4.9 are added to the

previously applied constraints. Thresholds considering the absolute and the relative wave

height are checked for plausibility. The effect of these constraints to optimization results is

evaluated.

6.5.1 Initial investigations

In the following the approach for parametric rolling and its two variants, applying a threshold

(i) for the encountering wave height and (ii) for the relative wave elevation at FP, COG, and

AP are investigated. Tab. 15 summarizes the constraints as posed for these investigations.

Besides general constraints controlling machine overload, geographic feasibility, and the

maximum acceptable duration of the journey also the constraints for slamming probability

and motion sickness incidence are used. These are set to thresholds identified in the previous

sections. All parameters used within the assessment of parametric rolling account for the

number of times per hour a situation occurs that may excite this motion. Therefore the

threshold names start with nph and the ending indicates if the criterions for FP, COG, AP, or

the critical wave height are faced (nph4hcrit, to be read as nph for hcrit, and nph@ap as nph at

AP).

Threshold values are set to <1>, i.e. maximal one exceeding per hour of the thresholds posed

in section 4.9 is tolerated.

Pareto Optimum Ship Routing, Chapter 6

Tab. 15: Applied constraints for the investigation of parametric rolling

parametric rolling

investigation

constrained

parameter content threshold

slprob slamming probability at 8% LPP

behind FP 1%

basic constraints

MSI motion sickness incidence 40%

nph@fp 1time x 6m

nph@cog 1 time x 5m

constraints

for the relative wave

elevation

nph@ap

times per hour a threshold value for

the relative motion between ship

and water surface is exceeded

1 time x 4m

constraints for the

absolute wave height nph4hcrit times per hour a threshold value for

the wave height is exceeded 1 time x 10m

Fig. 75 shows four PARETO frontiers. The right diagram of the figure is a cutout of the low

ETA region from the left diagram.

• The solid black line represents the optimization result where no constraints for

parametric rolling are posed.

• The blue-dashed line represents the result when constraints for the wave elevation

relative to the ship are set.

• Blue dots are related to constraints for the absolute wave height.

• The red-dash-dotted line represents the optimization result where both types of

constraints are used. It is congruent to the line representing the PARETO frontier for

the enabled constraint that observes the relative wave elevation.

Surprisingly the application of an additional constraint serves to improve the optimization

result. For the same ETA, routes of slightly lowered fuel consumption are found and ETA

itself could be reduced by 1h. In particular the constraints for the relative wave elevation seem

to serve in this regard. As a matter of course, the improved PARETO frontier stays feasible

when all parametric rolling constraints are deactivated. That means there is no obvious reason

that it was not found before, without applied parametric rolling constraints.

Fig. 75: PARETO frontiers for optimizations considering parametric rolling

In the right part of Fig. 75 three time minimum route designs are marked. Because of the

steep ascent of the PARETO frontier, two routes of the optimization without, and one time

Pareto Optimum Ship Routing, Chapter 6

minimum route for the run applying both types of parametric rolling constraints are chosen

for further comparison. Fig. 76 shows these three routes.

Those routes where no parametric rolling constraints are applied take nearly the same course.

The one where all parametric rolling constraints are applied takes a shorter course, which is

more northern in the beginning of the journey. Assumedly this again relates to the

characteristics of the genetic algorithm as observed in section 6.1, i.e. not completely

converging if achievable improvements are small.

A recalculation of fast routes taken from the optimization result with deactivated parametric

rolling constraints shows that several of these routes would be infeasible if these constraints

are enabled. Applying parametric rolling constraints serves to deplete the population in the

low ETA region. In this way selective pressure increases and the chances for alternative

solutions to become dominant improve. Although this is an unintended effect, it is beneficial

for the convergence of the optimization procedure.

Fig. 76: Routes of minimum ETA

In conclusion the attainable ETA and the fuel consumption do not vary greatly, and in this

regard the results of the optimizations remain nearly the same: An arrival on schedule is

impossible and the earliest possible arrival is around 143h. In contrast to the constraints for

MSI and accelerations that are very active on the shape of the PARETO frontier, the

constraints for parametric rolling are less dominant, i.e. there are only a few routes that are

ranked infeasible exclusively by these constraints and they are spread over the whole solution

space. Therefore the focus should be put on each of the applied parametric rolling parameter.

6.5.2 Assessment of the sensitivity for the parametric rolling parameter

It is obvious that the considered parametric rolling constraints influence the optimization

result, even if the effect on the PARETO frontier is relatively small. In the following step they

are focused in detail to assess their sensitivity and to get insight into their operating behavior.

For this purpose the result of the optimization applying all parametric rolling constraints is

evaluated more in detail. The number of feasible and infeasible routes for this optimization is

given in Tab. 16.

Pareto Optimum Ship Routing, Chapter 6

Tab. 16: Optimization applying all parametric rolling constraints

number of all designs 20.000

twofold generated designs 2.206

feasible designs 5.618

infeasible designs 12.176

thereof infeasible due to parametric rolling 3.276

All routes are evaluated at a number of discrete points. Regarding parametric rolling, for each

route point the number of threshold exceedings is calculated and the maximum value is given

back to the optimizer to decide whether a route is feasible or not. Except from about

30 routes, i.e. 1% of the infeasible routes due to parametric rolling, the most suspicious route

points according to the criterion that accounts for the absolute wave height are also the most

suspicious points that are identified by the constraints considering the relative wave elevation

at FP, COG, and AP. All the variants of the constraints posed for parametric rolling provide

almost the same ranking of route points regarding the constraint violations. However the

magnitude of the violation, i.e. the number of times the threshold is exceeded, is different.

Fig. 77 shows histogram plots for the designs that are infeasible due to parametric rolling

given in Tab. 16. The abscissa depicts the maximum number of exceedings for a particular

route. The height of the bars shows the number of routes with this maximum value. The

exceedings are grouped in bins from 1-5, 5-10, …, and 245-250. The most left bar represents

the route designs with 1-5 exceedings, those that are hardly infeasible. The ordinates are cut

off at number 400. Higher bars are described by giving the numeric value.

Fig. 77: Histogram plots for parametric rolling criteria

Pareto Optimum Ship Routing, Chapter 6

The histograms at the bottom (considering nph@ap and nph4hcrit) show a similar

characteristic. Most route designs are in the 1-5 bin and only a few show a higher number of

exceedings. The characteristic for COG is related too even if the number of routes in the 1-5

bin increases a little (nph@cog, upper right side). However, the histogram plot for the relative

wave elevation at FP (nph@fp, upper left side) is different from all others. The maximum bar

is at the 10-15 bin and also other bars nearby are relatively large. Further more there are some

routes that show up to 240 exceedings of the threshold per hour. It is interesting to observe

that the 1-5 bin is absolutely empty. Obviously the parameter that assesses the exceedings of

the relative wave height at the bow reacts much more sensitive. This could imply that the

threshold is set too harsh, on the other hand it is self-evident as in head seas like here, the

relative wave elevation at the bow is much larger than at other parts of the ship, because

encountering waves coincide with large vertical motions by pitching. Therefore high values

regarding the relative motion at the bow appear plausible. To step forward in this matter,

model tests or further simulations would be needed. Here the functioning of the parametric

rolling constraints will be illustrated by a further example and a resultant parameter variation

regarding the admissible relative wave height at the bow.

Fig. 78 compares the course of the feasible time minimum route from Fig. 76, the

optimization with all parametric rolling constraints, to an infeasible route due to at least one

active parametric rolling constraint close to this one, i.e. with a similar ETA and fuel

consumption.

Fig. 78: Feasible and infeasible route

Both routes take nearly the same course. Certainly, there is no information about the position

of the ship at a particular time, as no velocity profiles of the ships are given. For this reason

both routes are evaluated again step by step. The result is illustrated in Tab. 17. The route

point that caused the infeasibility of route id (7277) is identified. It is at 3d:3h:49min (75.8h)

after departure. Then, the closest route point in time of route id (16644), the feasible one, is

identified. It is at 3d:3h:21min (75.4h), half an hour earlier. At least both routes bring the ship

at the same time to the same sea area. As can be seen in Tab. 17 the wave conditions have not

changed except from the wave encounter angle. The latter one is the reason that the

parametric rolling parameter of the infeasible route have been evaluated according to

equations (4.55) and (4.56), whereas they are set to zero for the feasible one because the

encounter angle criterion is only active from 175 – 180°.

Pareto Optimum Ship Routing, Chapter 6

Tab. 17: Parametric rolling parameter for compared routes

route data feasible route

id(16644)

infeasible route

id(7277)

at time stamp [d:h:min] 3:3:21 3:3:49

H 1/3 [m] 4.3 4.3

TP [s] 14.8 14.9

enc. angle [deg.] 167 177

Vs [kn] 17 18

nph@fp [1/h] 0 11.35

nph@cog [1/h] 0 0

nph@ap [1/h] 0 0.01

nph4hcrit [1/h] 0 0.01

Fig. 79 shows both routes at the time given in Tab. 17, displayed in the wave field

interpolated for a time at 3d:4h. The figure serves to illustrate that both designs share the

same sea area at the same time and that both ships do not encounter extraordinary high waves.

Fig. 79: Parametric rolling, route comparison

Again the conclusion is that (i) either the constraint for the relative wave elevation at the

forward perpendicular is set too harsh or (ii) a harmless appearing situation may bear a risk of

parametric rolling. However, also both conclusions can turn out to be true. The first

conclusion has to be addressed in route optimization, whereas the second one relates to

operational aspects.

Regarding route optimization, section 6.5.3 presents results for variations of the threshold

value for the relative wave elevation at the forward perpendicular.

Regarding operational aspects, it is known that parametric rolling can cause abruptly

increasing roll angles and -accelerations. This motion behavior often occurs unexpected but

needs to be identified immediately to induce suited countermeasures. Analytical methods like

e.g. developed during the research project SinSee – evaluation of ship safety in severe seas,

cp. CLAUSS and HENNIG (2004), and finally implemented in several on-board monitoring

systems help to increase situation awareness and support in save navigation. In general polar

plots are employed to display hazardous operating conditions depending on wave height,

wave encounter angle, and ship speed. Fig. 80 gives two examples for these approaches.

Pareto Optimum Ship Routing, Chapter 6

Fig. 80: Polar plots, left side: limiting significant wave heights depending on ship speed

and encounter angle, right side: screen shot from OCTOPUS RESONANCE

(left side CLAUSS (2008), right side http://www.amarcon.com)

6.5.3 Increasing of the threshold value for the forward perpendicular

To take a step further in this matter, the threshold for the relative wave height at FP is

increased to 7m and 8m, respectively. For route id (7277) of Tab. 17 this would cause a

nph@fp of 3.2 respectively 0.74 exceedings per hour, i.e. for a threshold of 8m the route

becomes feasible. The histogram plots in Fig. 81 illustrate the effects on all routes that are

infeasible due to parametric rolling during the optimization. Already when increasing the

threshold to 7m, the 1-5 bin becomes the maximum one. Expectedly the number of routes in

the lower bins further increases, when the threshold is set to 8m. But in all cases, the

histograms for forward perpendicular FP show more routes with a high number of exceedings

than the histograms for COG and AP or those for the critical wave height.

Fig. 81: Histogram plots for a threshold for hrel of 7m and 8m

It has to be borne in mind that the used thresholds for the relative wave elevation are roughly

estimated. Anyhow, the constraint considering the relative motion at FP behaves different

from the other ones. Whereas a threshold of 6m produces a lot of infeasible designs a

threshold of 8m seems to be a quite risky and high value compared to Fig. 44.

Pareto Optimum Ship Routing, Chapter 6

On the other hand a threshold of 8m produces a similar histogram plot like the threshold for

the critical significant wave height and for this reason it seems permissible. Finally it would

be interesting to extend the focus on the relative motion within further investigations, i.e.

simulations and model tests.

Definitely, the proposed method is only a first step to establish a parameter that allows to

handle the risk of parametric rolling and to avoid hazardous situations. Considering the

assessment of threshold exceedings per hour it works and behaves suitable being used as

constraint. It serves to identify critical situations, and regarding the wave height it behaves

moderate, i.e. there are no jumps in the constraint functions. However, referring to encounter

angle, wave length, as well as encounter and pitch period, the proposed constraint functions

work like a switch, a behavior that is unfavorable to constraints. General in optimization it is

preferable to use constraint- and objective functions that show only minor changes of the

function value at minor changes in the design phenotype. To overcome this behavior, it is

possible to use e.g. a ramp instead of a jump function to pose a weight, related to the distance

of the current operating point to a range that is considered as critical. In this way a fuzzy

transition from feasible to infeasible can be realized.

In addition these constraints may become somewhat like the master’s choice. Based on the

numerical investigations the minimum safety requirements for a particular ship are posed.

Additional the master extends these requirements according to his/her experience.

6.6 Robust optimization

So far, analyzed weather data are applied as deterministic forecast. This will be

continued in the following, as no real deterministic forecast is available yet. The

comparisons of optimization results for deterministic-, mean- and ensemble forecasts

still use and also call the analyzed weather as deterministic forecast. This has no

influence on the argument and on the conclusions regarding the benefit of robust

optimization. Nevertheless this circumstance will be adjusted in section 6.6.4.

Robustness in optimization denotes the ability of an optimized design to stay optimal even if

essential parameters of the considered system are modified. In this consideration, a local

optimum with a great capacity to withstand parameter changes may be favorable compared to

a global optimum that loses its advantages if some system parameters deflect. Besides

numerical accuracy and the decision to consider some and neglect other system parameter, the

differences between the forecasted and real trend of weather are the most influencing

parameter for the robustness of a route optimization result. Therefore the focus in section 6.6

is put on the influence of probable weather changes on route optimization.

The assessment of robustness against weather changes uses ensemble forecasts, the set of 50

forecasts with an equal probability of occurrence. Here robustness will be used as:

• Constraint, i.e. a route has to stay feasible in a given number of ensemble forecasts.

• Objective, i.e. the number of ensemble forecasts where a route stays feasible is to be

maximized.

Besides, all other constraints and objectives stay active. Just like in the previous sections only

the rough sea condition at 20. – 30.1.2002 is considered.

Pareto Optimum Ship Routing, Chapter 6

To make the discussion of results more comprehensive and easier to read, the following

abbreviations are introduced:

• DF, MF, EF refers to

Deterministic- Mean-, and Ensemble Forecast, used

within the legends of figures.

• EF N refers to ensemble forecast number N, with N = (1…50).

• EF-id(M) refers to a route optimized by applying an ensemble forecast,

the route is identified by the ID number M, the number is

assigned by the optimizer.

• ETAmin-DF refers to the time minimum route identified by an optimization

applying the

Deterministic Forecast.

• ETAmin-MF refers to the time minimum route identified by an optimization

applying the Mean ensemble Forecast instead of the appropriate

deterministic forecast.

• ETAmin-EF N refers to the time minimum route identified by an optimization

applying ensemble forecast number N instead of the appropriate

deterministic forecast.

6.6.1 Deterministic, mean- and ensemble forecast

To step into this topic, the optimization result using the deterministic forecast is compared to

results produced by utilizing other forecasts types for the rough weather condition in January

2002. Fig. 82 shows the PARETO frontiers of the optimization applying the deterministic-

and simply the corresponding mean ensemble forecast. Nearly all constraints mentioned in

chapter 6 are posed, i.e. 1% slamming probability, maximum 1 event per hour that is

suspicious for parametric rolling and operability of the main engine. For the acceleration on

the bridge a threshold of 40% MSI is put. Here too, the constraints considering accelerations,

posed for MSI, and for main engine operability are the most active and deciding ones.

Fig. 82: PARETO optimum routes for deterministic and mean forecast

Pareto Optimum Ship Routing, Chapter 6

100

The optimization result firstly looks marvelous. A much faster journey than predicted by the

deterministic forecast seems possible, or with the same fuel consumption it is possible to

arrive 6 h earlier. Still ETAs like achieved in the medium or calm weather condition (114h,

117h) are not reached. But an ETA of 124h seems to be reachable and consequently an arrival

on schedule would be possible. As mentioned in section 3.3, caution is recommended because

the computation of a mean forecast implies the danger of damping extreme values. Therefore

it may roughly underestimate the coming weather situation. Fig. 83 shows the time minimum

routes of the optimizations applying the deterministic and the mean forecast. For comparison,

the shortest course is depicted as well. The right part of the figure shows the velocity profiles

for both routes. It can be seen that within the deterministic forecast a much lower speed level

is required. Fig. 84 depicts the occurring significant wave heights as predicted by the

deterministic-, mean-, and ensemble forecasts. On the left side the wave heights for the

ETAmin-MF route are shown, on the right the wave heights for ETAmin-DF respectively.

Fig. 83: Time minimum routes for mean- and deterministic forecast

Fig. 84: Significant wave heights for minimum ETA routes

in mean- and deterministic forecast

Both routes are only feasible within the corresponding forecast. As shown in Fig. 84 both

routes cross three areas of rough weather. ETAmin-MF gets infeasible in the third wave field of

the deterministic forecast (left diagram of Fig. 84). On the other hand, ETAmin -DF already

gets infeasible within the first rough weather area of the mean forecast (right diagram of Fig.

84).

Pareto Optimum Ship Routing, Chapter 6

101

The comparison of both routes, i.e. courses, velocity-profiles, and predicted waves, nicely

illustrates that a reduction of the ship speed does not only follow from current sea conditions

but serves also as a tactical maneuver to let strong wave fields pass. In this way it becomes

possible to save fuel, when faster variants would get infeasible due to wave conditions at a

later time of the journey. This at least gives an impression of the complexity of a routing

decision.

For example the route recommended within the deterministic forecast reduces speed at the

beginning of the journey to let the first storm field pass, Fig. 83 and Fig. 48. Later on it takes

a course more northerly and in addition it reduces speed again before crossing the tip of

Newfoundland. This is to reduce the impact of a third area of adverse weather.

Considered apart, both optimization results are plausible. To decide if the improvements

indicated by the mean ensemble forecast are trustable or not, the time minimum routes of both

optimizations are evaluated within the 50 corresponding ensemble forecasts.

The result is shown in Fig. 85. The diagrams of the left column belong to the ETAmin-MF

route, the right column depicts the results for the ETAmin-DF route, respectively. The three

most active constraints are considered. The diagrams in the upper line represent those for the

motion sickness incidence MSI, below those for the slamming probability, and at the bottom

line is the indicator for the main engine operability. Gray, horizontal lines represent

thresholds. Red-solid and blue-dashed, horizontal lines represent the values resulting from the

mean- and the deterministic forecast. Black dots represent the result for a particular ensemble

forecast. On the abscissa the ensemble number is given, on the ordinate the appropriate value

of the constraint.

Pareto Optimum Ship Routing, Chapter 6

102

ETAmin-MF evaluated in DF and EF ETAmin-DF evaluated in MF and EF

Fig. 85: Assessment of constraints in different forecasts

Regarding motion sickness incidence MSI and slamming probability splrob, even if

ETAmin-DF reaches the 40% MSI limit, both routes are feasible within both, the deterministic

and the mean forecast. In this case the main engine operability, represented by the specific

fuel oil consumption sfoc, is the crucial parameter. The deterministic forecast causes a

machine overload for the ETAmin-MF route and vice versa. Tab. 18 gives the numeric results

of this investigation.

Pareto Optimum Ship Routing, Chapter 6

103

Tab. 18: Optimized routes evaluated in ensemble forecasts

number of EF, rating

a route as infeasible ETAmin-MF ETAmin-DF

due to MSI 21 / 50 20 / 50

due to slamming 6 / 50 9 / 50

due to engine overload 36 / 50 24 / 50

due to all constraints 39 / 50 33 / 50

Finally 39 ensemble forecasts rate the ETAmin-MF route as infeasible and 33 the time

minimum route of the deterministic forecast. Admittedly the difference between both routes is

lower than expected and against all skepticism, the optimization utilizing the mean forecast

brings a remarkable result. At least, the mean forecast seems capable to assess if

improvements are possible, to be faster or less fast than the optimization result predicted

within the deterministic forecast. As long as the shift in time of maximum values of the

ensemble forecasts is small, this assumption may hold. Nevertheless, the spread in peak

periods and encounter angles is neglected from this point of view, although it is important for

the prediction of ship response.

Completing this first application of ensemble forecasts, two further optimizations are run. As

observed from Fig. 85, both time minimum routes of the optimizations above are feasible in

EF35 and infeasible in EF25. In a way these two forecasts represent edge cases of probable

weather development, both have the same probability of occurrence like the deterministic

forecast. Therefore these two ensemble forecasts are taken as forecast for a route

optimization. Fig. 86 shows the resulting PARETO frontiers and the appropriate time

minimum routes ETAmin-DF, ETAmin-MF, ETAmin-EF25 and ETAmin-EF35.

Fig. 86: PARETO frontiers and optimum routes in different forecasts

Both ensemble forecasts show an improvement of the PARETO frontier. Although

ETAmin-MF is feasible in ensemble EF35, the PARETO frontier of the optimization applying

EF35 does not match this design. A first hint for the reason is given by comparing the courses

in Fig. 86. The course of ETAmin-EF35 clearly differs from all other courses. Here the genetic

algorithm converges in a local optimum and it is not able to jump out of this suboptimal

region. A comprehensible matter as the probability to find a favorable design decreases the

more the phenotype has to change.

Pareto Optimum Ship Routing, Chapter 6

104

Further on, the course of ETAmin-EF25 is quite close to ETAmin-DF, nevertheless it requires 11

hours less. Accidentally ETAmin-DF is infeasible in EF25 although it lies in the region of

probable feasible solutions. This is no reason for concern, during an optimization most

infeasible designs are inside the area of feasible solutions, i.e. lying above and to the right

side of the PARETO frontier does not automatically include being feasible. On the other

hand, all designs below and left from a time minimum design are infeasible. Regarding the

PARETO frontier of the optimization within EF25, all PARETO designs of the optimizations

applying MF and EF35 that lie below or at the left are infeasible. However, as

aforementioned, there is no proof for a convergence.

Definitely the robustness of an optimization strongly depends on the quality of the weather

forecast. Fig. 86 clearly shows the diversity of possible solutions. In this case, the

deterministic forecast optimization seems to overvalue the achievable ETA. It is a common

strategy to repair this weak spot by re-optimizing from the current waypoint of a journey

using updates of the weather forecast. Preferably there are measurement devices on-board the

ship to recognize the difference between the present wave field and the forecasted one.

Thinking of an intelligent ship management system it would be possible:

• In the case of deteriorated weather, to demand for suited countermeasures, e.g. reduce

the ship speed, consequently falling behind the aimed route point, assess a reachable

route point, and optimize again with an updated forecast.

• In the case of more convenient weather, it is surely possible to reach the aimed

waypoint and furthermore it can be possible to step beyond. Admittedly, as long as

optimizations applying a genetic algorithm are time consumptive, a faster algorithm is

favorable as it is necessary to quickly optimize from the achieved waypoint applying

an updated forecast. Certainly the convergence characteristic should be regarded, too.

En-route optimizations with updated forecasts will remain necessary as long as the weather

pattern deviates from the forecasted weather. Here the benefit of the ensemble forecast system

should be focused more in detail in order to assess its capability to give the most

comprehensive insight for a route decision prior to the departure.

6.6.2 Robustness as constraint

Even by applying route optimization with updated forecasts, the actual achievable ETA is not

known until the end of a journey. Using ensemble instead of deterministic forecasts may serve

to get more comprehensive insight to the future and to assess achievable ETAs right at the

beginning. Therefore a further constraint is posed that covers all constraints considered within

the route optimization, called C_nef (Constraint observing the threshold posed for the number

of ensemble forecasts where a route has to stay feasible). This constraint counts in how many

of the 50 ensemble forecasts a particular route stays feasible. It becomes active when the

number of ensemble forecasts where the considered route stays feasible falls below a

threshold value. Tab. 19 shows the constraints posed for the following optimizations. The

number of ensemble forecasts where a route should stay feasible is varied.

Pareto Optimum Ship Routing, Chapter 6

105

Tab. 19: Constraint for robust route optimization

constraint name criteria threshold

C_slprob slamming probability 1% at 8% LPP behind FP

C_MSI motion sickness incidence 40%

C_spef main engine operability accordance with main engine

C_nef

number of ensemble forecasts

where a route stay feasible due

to C_slprob, C_MSI and C_spef

20, 25, 30, 35, 40 or 45

Differences in the weather conditions of the ensemble forecasts and resulting different

additional resistances due to waves for a given velocity and course angle cause a spread in the

fuel consumption of a considered route. For this reason the objective function uses a mean

fuel consumption instead of the fuel consumption. It is the mean value of fuel consumptions

taken from the ensemble forecast evaluations where a considered route stays feasible.

Infeasible results are left out of consideration. The objective function regarding the mean fuel

consumption is built by simply replacing the fuel consumption in equation (5.2) by the mean

fuel consumption.

Fig. 87 compares the PARETO frontier of ETA and fuel consumption of the optimization

applying the deterministic forecast with a frontier for ETA and mean fuel consumption for the

same routes.

Fig. 87: Fuel consumption and mean fuel consumption

A PARETO design is represented by a red square, a black square in vertical direction

represents the corresponding mean fuel consumption. Furthermore the figure contains

information about the robustness of the optimized routes against weather changes. The

number of blue crosses in the vertical direction relates to the number of ensemble forecasts

where this route stays feasible. If there are no markers vertically to a PARETO design, it is

infeasible in all ensemble forecasts. The spread in the fuel consumption may be utilized as

well to value robustness, however, this is not done here to keep things simple. All in all the

mean fuel consumption deviates not much from the fuel consumption, whereas the spread, i.e.

possible differences in the fuel consumption, is remarkable. Here 10% of the overall fuel

consumption simply depends on the weather development!

Pareto Optimum Ship Routing, Chapter 6

106

Fig. 88 shows PARETO frontiers for route optimizations in ensemble forecasts applying

different thresholds considering the number of ensemble forecasts where a route stays

feasible. In addition the PARETO designs of the deterministic forecast, evaluated for mean

fuel consumption, taken from Fig. 87, are shown. The right part of the figure shows an

enlarged sector of the left part. Minimum attainable ETAs are pointed out. Contrary to

expectations the PARETO frontier for routes feasible in 25 ensemble forecasts lies above the

frontier for routes that are feasible in 30 ensemble forecasts. This for sure is not plausible but

with regard to the PARETO frontiers depicted in the left part of Fig. 88 this defect is small. It

may be due to a limited capability of the optimization method, and therefore it should be

resolvable by further adjustments of the optimization method to the routing problem.

Fig. 88: PARETO frontiers for optimizations

in a deterministic and an ensemble forecast

As already discovered in section 6.6.1, the optimization within the deterministic forecast

seems to roughly overvalue the attainable ETA. Increasing the threshold of C_nef from 20 to

40 surprisingly increases ETAmin just from 124h to 129h. Only a further increase of the

threshold leads to remarkable increases in ETAmin. The PARETO frontiers for all thresholds

lie close together. Obviously, for a particular ETA the differences in the mean fuel

consumption are clearly lower than the total spread due to the ensemble forecasts, as given in

Fig. 87.

Compared to the optimization result of the deterministic forecast, even the time minimum

route that is feasible in 45 ensemble members arrives two hours earlier than the one of the

optimization within the deterministic one, feasible in only 17 ensemble members, cp. Tab. 18.

Furthermore, considering the fastest routes of the deterministic forecast optimization, the

predicted fuel consumption within the ensemble forecasts is much lower. Fig. 89 illustrates

the optimization results by plots of the courses and velocity profiles for the fastest routes from

the optimizations with different thresholds for C_nef (20, 30, and 40) and from the

deterministic forecast optimization.

Pareto Optimum Ship Routing, Chapter 6

107

Fig. 89: Courses and velocity profiles of routes optimized

in a deterministic and an ensemble forecast

Even if the courses seem to stay close together, there are deviations of about 100nm, a

remarkable distance considering that extreme wave fields are circumnavigated here.

Regarding the velocity profiles, in particular the difference between the deterministic and

ensemble forecast optimization results become obvious. Within the deterministic forecast the

ship speed has to be reduced much more to avoid heavy weather. Considering the ensemble

forecast, the differences of the velocity profiles appear plausible, e.g. maximum speed

reductions can be found at highest threshold values of the robustness constraint.

Finally all optimized routes clearly deviate from the shortest possible journey, the great circle

route. So maybe a master who utilizes the deterministic optimization result can sail faster due

to the current weather conditions. Finally he/she enters port with a comparable ETA as

predicted in the ensemble forecast optimization. But taking the shortest track would clearly

not succeed.

6.6.3 Robustness as objective

During the optimizations regarding robustness as a further constraint, all in all 126.000 routes

are evaluated (threshold of C_nef varied from 20 to 45, i.e. 6 optimizations and in each case

300 generations with 70 members). Applying robustness as an objective, this reduces from six

to one optimization. Here 500 generations each with 100 individuals are run, all together

50.000 routes are evaluated. Although this setup saves a lot of computational effort, the prime

reason for this investigation is to support the observations made above. Similar to the usage as

constraint, simply the number of ensemble forecast members where a route stays feasible is

applied as objective, called O_nef. To produce a value for this objective comparable to the

magnitude of the both others, the one for mean fuel consumption and ETA, it is set up as:

1













=∑

nNEFObjective , (6.4)

where ni = 1, if a route is feasible in the respective ensemble forecast and 0 otherwise.

Fig. 90 illustrates the result of an optimization utilizing robustness as objective of the

optimization. The left part of the figure shows the mean fuel consumption versus ETA, the

right part provides information on attainable ETAs relating to robustness against weather

changes.

Pareto Optimum Ship Routing, Chapter 6

108

Fig. 90: PARETO frontiers of the ensemble forecast optimization

Regarding the PARETO frontiers on the left side of the figure, the mean fuel consumption is

similar to the fuel consumption found when using robustness as constraint, as shown in Fig.

88. Here the PARETO frontiers for different threshold values are produced by filtering all

evaluated route designs following to the optimization. The right part of the figure depicts the

diagram for the objectives O_nef and ETA. Blue dots represent the time minimum routes for a

particular number of ensemble forecasts nef that rank these routes as feasible. PARETO

optimum designs regarding these two objectives are connected with a red line, i.e. for these

designs an increasing robustness deteriorates the time of arrival. Two route designs are

emphasized by green squares. Both are quite fast combined with reasonable fuel consumption,

as they are located closest to the PARETO frontier for ETA and the mean fuel consumption.

They would match the schedule, and furthermore they possess a high robustness against

weather changes. As can be seen in Fig. 90, decreasing ETA would rapidly decrease the

robustness, increasing robustness would strongly increase ETA. In a sense these two route

designs represent “over all optimum routes.“

Fig. 91 depicts courses and velocity profiles of these two routes. To compare, the ETAmin

routes of the deterministic and the one of the mean ensemble forecast are shown, too.

Fig. 91: Courses and velocity profiles of routes with different robustness

Except from ETAmin-DF, all of them are PARETO optimum regarding ETA and fuel

consumption. Even if the ETAmin-DF is quite slow and ETAmin-MF is the fastest one here, both

are less persuasive regarding the robustness against weather changes; the DF-optimized route

Pareto Optimum Ship Routing, Chapter 6

109

is feasible in 17, the MF-optimized in 11 ensemble forecasts, respectively. This clearly shows

the advantages of including robustness into a route optimization. Applying a mean forecast

does not automatically produce routes that are feasible in 50% of the ensemble forecasts;

slower routes are not automatically more robust.

It finally proves true what is indicated by the mean ensemble forecast optimization, compared

to the optimization in the deterministic forecast it is possible to decrease the duration of the

journey. Furthermore it is possible to simultaneously increase the robustness against weather

changes and reduce ETA. The required fuel consumption of the routes selected from the

ensemble forecast optimization (favorable designs in Fig. 90) is in between the fuel

consumption of ETAmin-DF and ETAmin-MF. Nevertheless, except from ETAmin-DF, they all

represent an optimum regarding the fuel consumption as they are close to the PARETO

frontier (Actually the PARETO frontier given in Fig. 90 is a PARETO surface now. It

represents a 3-dimensional border of the solutions space for three objectives, O_eta, O_fuel,

and O_nef where an improvement in one objective deteriorates or at most preserves the other

objectives).

Both optimum routes of the ensemble forecast, presented in Fig. 91, distinguish clearly in

course and velocity profile from their deterministic- and mean forecast counterpart. Regarding

further research, it would be interesting to see if a route optimization applying updates of a

deterministic forecast would lead to a comparable result.

6.6.4 Rectification for using analyzed weather as deterministic forecast

Unfortunately, there is no real deterministic forecast available yet and therefore the re-

analysis is used instead. This in fact is no mistake, as the re-analysis surely could have been

the real deterministic forecast. Consequently all conclusions made above are valid.

Fortunately, these weather data are much more than a forecast, they represent the best

available description of the real weather condition during the considered period.

Now, the optimization results for this forecast / re-analysis should be taken as what they are,

the re-calculation on routes that are possible at the best.

Regarding the PARETO frontiers of Fig. 86 and Fig. 90 it is obvious that the ETAmin-MF and

the two favorable routes selected from the ensemble forecast optimization are not feasible in

the weather condition that existed according to the re-analysis. For all of them one or more

constraints that disapprove these routes would become active. This implies that most likely,

suited countermeasures like course or speed changes would have become necessary.

Consequently the ship leaves the planned route and the aspired ETA cannot be kept. In this

case, a re-optimization from the current waypoint to the port of destination serves, e.g. to find

the fastest track at maximum fuel savings. In doing so it is not possible to get comprehensive

insight to an attainable ETA at the beginning of the journey. However, this is favorable e.g.

for a conscious disposition of cargo. To step forward in this matter, Fig. 92 serves to evaluate

the optimization results of the MF- and the EF-optimization. The figure shows the velocity

profile and constraints for ETAmin-MF and for both routes selected from the ensemble forecast

result EF-id(61) and EF-id(70) at re-analyzed weather conditions. Since the constraints for

parametric rolling did not become active, they are left out of consideration. The appropriate

courses can be seen in Fig. 91.

Pareto Optimum Ship Routing, Chapter 6

110

Since ETAmin-MF takes a shorter track at higher

ship speed, it is able to arrive 6h earlier than

EF-id(70). Consequently the fuel consumption is

higher but the route is anyway close to the

attainable PARETO frontier.

Regarding the specific fuel oil consumption spef,

both EF-routes stay feasible within the analyzed

weather. ETAmin-MF becomes infeasible due to

main engine overload. During the first and the

fourth day spef of ETAmin-MF exceeds the upper

level of the diagram. In these cases the load to the

main engine exceeds the nominal MCR. In this

regard ETAmin-MF becomes infeasible without

any doubt.

Concerning the motion sickness incidence MSI

and significant amplitudes of vertical

accelerations on the bridge, the ETAmin-MF route

does not violate any constraints. Here both routes

taken from the ensemble forecast optimization

produce active constraints. Three times the

threshold for MSI and two times the one for

significant amplitudes of vertical acceleration is

exceeded. However, these violations are small

and limited to a short period. Therefore they may

be tolerated. Admittedly, the graphs for zacc and

MSI look very similar. This becomes clear as

steps of the route evaluation are arranged more or

less equal in time by the routine for the route

perturbation. Regarding the acceleration around

t = 20h, a period of high accelerations spanning a

longer time can be observed. This certainly

produces a much higher MSI like reckoned here,

and an integration of time steps for the

assessment of MSI appears recommended.

Related to the slamming probability there seems

to be no risk for none of the routes to be affected

by slamming.

Fig. 92: Ensemble results

in analyzed weather

Pareto Optimum Ship Routing, Chapter 6

111

Altogether the optimization result produced by means of the ensemble forecast is the more

reliable one. It is possible to execute the journey according to the optimization result when

slightly higher accelerations are tolerated. Admittedly, this is just one example applying one

weather scenario and it is necessary to extend these investigations to more weather cases, to

see if the drawn conclusions hold. On the other hand, the computational effort for an

ensemble optimization is more than fifty times higher than an optimization applying a single

forecast. Notwithstanding all thinkable improvements to accelerate the optimization

procedure, an application of the current set-up for the ensemble route optimization on-board a

ship seems not possible, as it is too time-consuming. Nevertheless, regarding the explanatory

power, the approach to employ ensemble forecasts should be continued. Up to now it turned

out to be a reasonable attempt to look out for the future and to predict the unpredictable.

Pareto Optimum Ship Routing, Chapter 7

112

7 Summary and outlook

7.1 Initiation of a new approach

As shown by the operators of the Dutch Ship Performance Optimization System SPOS, it is

possible to achieve a significant reduction of time operating in bad weather by applying a

weather routing system. Weather routing and route optimization are meaningful methods to

improve the safety on-board a ship regarding threads by adverse weather. Dependent on the

set-up of objectives for a route optimization, it is furthermore possible to save operational

costs, mainly fuel- and lubrication oil, in order to find the fastest track or to save fuel for an

arrival on schedule. Inspired by the work of SEAROUTES the set-up of a new route

optimization approach is decided. This approach is able to handle multiple objectives

separately, as including opposing objectives into a single objective function becomes difficult

to comprehend. Already traveling time and operating costs, major objectives of any route

optimization and here represented by the estimated time of arrival ETA and the fuel

consumption of the main engine, depict opposing targets, as improving one normally

deteriorates the other. This suggests an optimization that is able to identify PARETO

optimum routes. As long as an arrival on schedule is possible the optimization has to identify

routes that comply with the schedule at minimum operational costs. When weather worsens

and an arrival on schedule becomes impossible, the optimization should depict routes of

minimum operational costs for achievable ETAs. That is why the knowledge about feasible

PARETO optimum routes, as depicted in Fig. 93, supports a conscious decision making in

terms of safer ship operation, reducing operational costs and the ability to rearrange the

sailing list if necessary.

Fig. 93: PARETO frontier and time minimum route at rough seas

Expectedly, at calm weather the shortest route is the optimum one. When weather worsens it

is mostly possible to maintain the optimum track by retaining the shortest route and reducing

speed if necessary. In a rough weather situation the routing decision becomes more complex,

as the shortest route is not necessarily the best one anymore. The optimum becomes a trade-

off between (i) additional fuel consumption due to a longer course and (ii) fuel saving by

reducing the additional resistance due to waves by means of circumnavigation of strong wave

fields. For the rough weather scenario, considered within this study, it turns out that applying

optimization serves to reduce the fuel consumption up to 15 %. Regarding the demands for

the optimization process, it is expected that objective functions show humps and hollows or at

Pareto Optimum Ship Routing, Chapter 7

113

least become quasi-multimodal when optimization results are governed by constraints.

Therefore the optimizer should be able to overcome local optima. Generally stochastic

methods serve for this purpose. Here a multi-objective genetic algorithm, available within the

generic optimization environment modeFRONTIER, is applied.

As a rule of thumb, it turned out that a stochastic method requires the tenfold number of

designs compared to their deterministic counterpart. Approximately half of the evaluated

designs are necessary to be sure that a global optimum is detected. Anyhow there is no proof.

A frequent requirement during the set-up of an optimization procedure is the reduction of

computational effort, generally to accelerate the procedure. In the case of genetic algorithms

the number of designs within one generation depends on the number of free variables and the

number of objectives. Furthermore, each route design has to be evaluated at a number of

points that are tight enough to recognize the weather pattern. Altogether, a fast method to

assess the ship behavior in waves is necessary. Therein the modeling of ship routes should

require as few free variables as possible. For this purpose the mathematical description of a

route design and its evaluation base on:

• a B-spline technique to model course and velocity profile of a route,

• standard spectra to describe the seaway,

• the theory of linear superposition to assess ship responses.

The applied B-spline modeling technique provides a small set of free variables that enables a

smart modeling of the course and the velocity profile. It is possible to extend the method to a

finer spatial and temporal resolution or to cut-off areas that do not allow passages, like land or

islands. In this way the perturbation for the maximum northern shift in Fig. 13 is controlled.

Fig. 13: Parent route and maximum perturbations

To describe predicted wave conditions, deterministic- and ensemble forecast data from the

European Centre for Medium-Range Weather Forecasts ECMWF are used. In both cases they

consist of significant data that are used in combination with standard spectra. A PIERSON-

MOSKOWITZ spectrum provides good results.

Pareto Optimum Ship Routing, Chapter 7

114

Fig. 94: Significant wave heights of a severe winter storm

in North Atlantic, 22. January, 2002

For the assessment of ship responses in waves, transfer functions for the ship motion and

response functions for the added resistance due to waves are applied. They are calculated by

means of the well-established strip theory code SEAWAY that is based on potential theory,

JOURNÉE (2001). Response functions for the added resistance due to waves are calculated

according to the integrated pressure method of BOESE (1970). Together with the calm water

resistance according to HOLTROP and MENNEN (1984), a characteristic of a similar

propeller provided by YASAKI (1962), and the main engine characteristic,

MAN B&W (2000), the fuel consumption for particular service conditions in a seaway is

assessed.

Transfer functions are used for the set-up of constraints within the optimization. They are

calculated for various points on-board the ship and serve to assess absolute motions like the

acceleration on the bridge and motions relative to the water surface, e.g. applied to a

constraint for the slamming probability.

Inspired by the current discussion on parametric rolling, an additional module is established to

detect situations that are suspicious for such events. The proposed model is not

straightforward strip theory, but a combination of strip theory and results from sophisticated

capsize simulations, provided by KRUEGER et al. (2006). In this way it becomes possible to

identify hazardous situations with reasonable computational effort.

It has to be pointed out that the applied constraints represent an arbitrary selection used to

demonstrate the feasibility of a new approach for route optimization. A sophisticated decision

support system should include all available features to assess the behavior of a ship in waves,

otherwise a route recommendation would become delusive by pretending a safety that does

not exist. In this regard it is a trade-off between numerical accuracy and computational time,

e.g. between the achievable accuracy of employed methods, the accuracy of underlying data,

and the computational effort to predict the ship behavior on a particular route.

Regarding the accuracy of underlying data, one of the most important topics in ship route

optimization is the forecast uncertainty, i.e. the increasing deviation of the upcoming from the

forecasted weather, the more the forecast reaches into the future. Deterministic medium range

forecasts are considered to be reliable for three days. The rough weather example applied in

this study implies that, under severe conditions, this assumption may not hold. Normally

en-route countermeasures are necessary to avoid danger to the crew or damage to the ship and

the cargo.

Pareto Optimum Ship Routing, Chapter 7

115

For conscious route planning and disposition of cargo it is favorable to get the most

comprehensive insight to the upcoming weather as soon as possible. To step forward in this

matter, the robustness of optimized routes to weather changes is assessed by means of

ensemble forecasts. Finally a definition of robustness against weather changes is established.

It is used to extend the optimization approach by applying robustness as further constraint or

objective. In this way it becomes possible to identify PARETO optimal routes that, in contrast

to the optimizations applying the deterministic forecast, comply with the schedule for a

reasonable violation of constraints, cp. Fig. 95.

Fig. 95: Rough weather route optimization

for deterministic and ensemble forecasts

Pareto Optimum Ship Routing, Chapter 7

116

7.2 Discussion on the performance

Various sensitivity analyses serve to assess the capability and the reliability of the route

optimization approach. Therein, the following aspects are addressed:

• Considering different weather conditions serves as a first plausibility check for the

optimization setup.

• Applying different spectra, namely a PIERSON-MOSKOWITZ and a JONSWAP

spectrum, is used to assess the error induced by the approximate description of swell.

• The required accurateness for the representation of the hull shape, within strip theory

based motion prediction, is addressed by hull form modifications.

• Threshold levels and constraints are evaluated by variations of thresholds and by

enabling and disabling constraints while the remaining parameters of the optimization

are held constant. In some cases reasonable modifications are installed,

e.g. for the assessment of the slamming probability.

• Benefits and performance of two add-ons, addressing parametric rolling and

robustness against weather changes, are evaluated by threshold variations.

• Besides all these investigations it is possible to assess the performance and usability of

strip theory and of the genetic algorithm within route optimization.

Regarding the variation of the weather condition, the optimizations show plausible results. As

sea conditions worsen the fuel consumption increases for a given ETA. At the same time the

minimum attainable ETA increases as well. Actually this is just an initial test for a qualitative

assessment of the optimization result. However, regarding the fuel consumption the results fit

very well to data measured during the SEAROUTES project. Therefore the resistance

prediction by means of the method provided by HOLTROP and MENNEN yields good

results for the design draft. Furthermore the usage of a similar propeller is admissible. At

minor draft load condition, e.g. ballast load condition, the calm water resistance is assessed

too low.

Regarding the modeling of the sea state, a method considering simply swell is sufficient as

long as the inaccuracies of other applied methods dominate the quality of the optimization

result, e.g. the uncertainty of the forecast is the much more influencing parameter. Surely,

applying 1-D or even 2-D spectra can improve the accuracy of the calculated ship responses.

But it has to be borne in mind that resulting little additional time requirements during the

evaluation would summarize to significant amounts when computed repeatedly. Therefore

simply significant data together with a PIERSON-MOSKOWITZ spectrum are preferred here.

A suited compromise, to improve the modeling of swell and to keep the amount of data in a

reasonable magnitude, could be the usage of parameters for scaling the peak width similar to

those applied in a JONSWAP spectrum. In that way significant wave data would consist of

data fields for wave heights, -periods, and -direction and additional three to four parameters

that serve to scale the spectrum, given on the same grid as the other ones. Finally, as long as

wind forces and currents are neglected, the application of a PIERSON-MOSKOWITZ

spectrum can be regarded as satisfactory.

Pareto Optimum Ship Routing, Chapter 7

117

Fig. 56: Typical storm spectra for the North Atlantic,

PIERSON-MOSKOWITZ, JONSWAP and ECMWF-1D spectrum

During an optimization about 20,000 routes at 60 route points are investigated. This amounts

to 1.2 million calculations of ship responses in a particular sea state. Obviously a fast method

for the assessment of ship responses is needed. This implies that the method of linear

superposition is recommended, applying wave spectra, transfer functions, and their statistical

evaluation. Here transfer functions for a particular load condition calculated by means of

strip-theory are used. It turned out that there is no need to consider the changes regarding

departure and arrival condition, i.e. to consider fuel consumption during the journey. The

usage of an intermediate load condition is absolutely sufficient.

Even if strip theory does not take the hull form above the water line into consideration, a

proper hull description is recommended, as is shown by hull form variation studies. However,

the upper hull shape has a clear impact on the seakeeping behavior of a ship. In some cases

the upper hull has to be included, e.g. by methods like the proposed parametric rolling module

that is a combination of strip theory with results from more sophisticated numerical

investigations. To include more of the ship behavior that is definitely not covered by the

numerical representation, it is still possible to employ the master’s knowledge and experience

to build constraints that represent so to speak master’s choice. This is performed for example

in the SPOS routing system.

The initial setup uses constraints for the load to the main engine, the vertical acceleration on

the bridge and for the slamming probability. For the latter two, thresholds comply with the

thresholds posed in the European project SEAROUTES (2003). It turns out during the

sensitivity analysis that the constraint posed for vertical acceleration on the bridge is the most

deciding one for the optimization and for the shape of the resulting PARETO frontier.

Comparisons to thresholds posed by NORDFORSK (1987) indicate that the applied threshold

might be too low. For this reason the model is extended by applying a constraint for the

acceleration on the bridge that is based on MSI. Including amplitude, frequency, and duration

of the exposure, the MSI provides a more comprehensive evaluation. In that way, a second

indicator for the load on the crew by means of accelerations is provided for a more

differentiated assessment. Following investigations show that the load on the crew in an

adverse weather situation can be enormous and with regard to maintain situation awareness,

the thresholds for vertical acceleration on the bridge are kept tight.

Pareto Optimum Ship Routing, Chapter 7

118

Considering the threshold for the slamming probability, no influence on the shape of the

PARETO frontier is recognized, even if this constraint is tightened significantly. Here the

study shows that commonly applied methods to assess slamming could need revision. It is

known that Panmax container vessel slam, but this is ignored by the initially applied

procedure, based on the method of OCHI. It turns out that the method could be improved by

permitted tuning of thresholds and parameters. But still the influence of this constraint to the

shape of the PARETO frontier is small, at least it affects the optimization process. Obviously

the concept of slamming probability deserves further investigation. It is the question whether

modern structures in ship design withstand the criteria that formerly were posed for slamming

or not.

The constraint for the load to the main engine is not changed, as the operating point has to

comply with the characteristic of the main engine. However, this constraint has a high impact

on the optimization result, similar to the one for vertical acceleration.

Regarding parametric rolling, first sensitivity analysis show that it is possible to establish and

adjust thresholds applying the absolute wave height and thresholds for relative wave height in

a reliable way. Furthermore wave length, encounter angle, and the period of encountering

waves, related to the pitch period of the ship, are applied to successfully identify situations

that are suspicious for parametric rolling. Similar to the constraint for the slamming

probability the constraints for parametric rolling show only minor influences to the shape of

the PARETO frontier. So far the operating behavior of the constraints posed for parametric

rolling are plausible, surely the reliability is to be proved comprehensively before the

application to a decision support system can be recommended.

Considering the robustness of optimized routes against probable weather changes, the

usability of a mean-ensemble- and an ensemble forecast is investigated. Achieved

optimization results are compared to results of optimizations utilizing the deterministic

forecast. For the weather case considered here, applying a mean-ensemble forecast indicates

that a considerable reduction of the duration of the journey is possible. Still these

optimizations yield no information whether optimized routes improve in terms of robustness.

For this purpose the ensemble forecast is applied and robustness is included as constraint as

well as objective to the optimization. In this way, routes that are PARETO optimum regarding

ETA, fuel consumption, and robustness against weather changes can be identified. Two

superior routes are taken for further investigations. Regarding the high number of ensemble

forecasts (84%) that rank these faster routes as feasible it seems reasonable to say that even if

there is a probability of 16% that constraints are violated, the severity of violation will be

small as its probability is reasonably low. Observations of constraints along optimized routes

in analyzed weather supported this assumption. Anyhow, this aspect needs to be proved by

further investigation. So far, employing ensemble forecasts appears to be a reliable method

for a conscious decision support. The major obstacle that has to be solved for a practical

application is reducing the required computational time.

For this first approach it is absolutely sufficient that constraints simply control the maximum

values of an investigated route. More sophisticated constraints that also account for the

duration of an incidence can be included quite simple. This is shown by the concept of MSI.

However, it is an important matter of the set-up to find reasonable threshold values. For

vertical acceleration on the bridge it turns out to be meaningful to account for the duration of

the exposure. Regarding e.g. the slamming probability, there is no need to account for

durations as the constraint already includes a time dependency.

Pareto Optimum Ship Routing, Chapter 7

119

The applied generic optimizer, here a genetic algorithm, produces reasonable results even if it

is not at all adapted to the route optimization problem. If merit functions are shaped relatively

flat, the genetic algorithm shows problems to identify an optimum. This is because the

probability of a design to be considered for producing the next generation is coupled to the

fitness relative to other designs and not to its absolute fitness. Up to certain limits, this can be

overcome by applying objective functions instead of the objectives themselves. Nevertheless,

an adaptation of the applied stochastic method is recommended to accelerate the code and to

yield convergence improvements. It turns out that even if constraints do not directly affect the

PARETO frontier, they are able to influence the optimization. In the case of MSI or

slamming, the enabled constraint helps to improve the optimization process, local optima are

overcome a lot easier, and consequently the procedure converges faster to the desired

PARETO designs. Even if the number of infeasible designs during a rough weather

optimization is already high, additional constraints can stabilize the optimization. In this

regard, improvements of the optimization result are observed when parametric rolling is

included. Applying this constraint, serves to deplete the population in the low ETA region.

Consequently selective pressure increases, i.e. making an achieved optimum less comfortable.

This serves to keep the genetic pool diversified and the chances for alternative solutions to

survive become more probable. As the probability to produce different route variants during

the recombination increases, it becomes easier for the optimization to overcome local optima.

Although this is an unmeant effect, it is beneficial for the optimization.

In the present, provisional setup the investigation of one route takes 2.5s, i.e. the whole

optimization would take nearly 14h. Fortunately, repeated designs are investigated only once,

this reduces the computational time to approximately 10h. This is still too time consuming for

an application on-board. However, there has been no runtime optimization in terms of

integrating the whole process into one program. Currently separate programs perform the

perturbation, the evaluation, and the optimization of a route and the data transfer takes a lot of

computation time. There is definitely a high potential to accelerate the procedure. However, if

the computational speed of workstations increases like it has done in last years, an application

of the proposed optimization method should be possible within the near future. On the other

hand, it is not necessary to apply a genetic algorithm. Other stochastic methods, e.g. particle

swarm optimization that is known for its high adaptability and speed, may serve better. This

for sure, is an important topic to be addressed in further research, in particular when an

integration of the routing system on-board a ship is desired.

7.3 Conclusions

Within this study the set-up of a new ship route optimization approach is documented.

Benefits and obstacles of the method are pointed out and discussed. The introduction of

further constraints and features, like e.g. wind resistance, current, and wind waves, is

straightforward. Regarding the method itself, enhancement of the numerical performance is

the most important topic for further improvements. To be competitive to other methods, a

considerable reduction of computation time is needed. As a matter of fact, a genetic algorithm

based method will never reach the computational speed of a deterministic or graph theory

based method. This disadvantage is compensated by some clear advantages. Instead of one

optimum route design, the result consists of a number of PARETO optimum routes.

As the optimization method is able to handle more than one objective, new objectives can be

easily introduced. For a merchant ship ETA and fuel consumption are major aspects for a

route decision. Cruise liner may prefer ETA and MSI as objectives to stay on schedule and

assure the well-being of passengers. By employing ensemble forecasts, also robustness

Pareto Optimum Ship Routing, Chapter 7

120

against weather changes can be used as further objective. By means of the PARETO ranking

it is possible to come to a more conscious decision. The master is provided with a complete

overview on what is feasible for which price or risk.

The topic of computational effort is closely connected to the discussion on where to situate a

decision support system. Requirements on computational power suggest a positioning ashore.

On the other hand, this is far away from the ship in operation, and a feedback of the master or

other related devices becomes difficult. Following the current development in ship operation,

integrated monitoring and decision support systems will clearly take place on-board. With

regard to uncertainties of a weather forecast, there are good reasons to integrate the route

optimization functionality into an integrated monitoring and decision support framework.

Forecasted and nowcasted weather provided by meteorologists and direct wave measurements

of a ship monitoring system give the officer a comprehensive view on current and upcoming

weather conditions. Route optimization supports to find the optimum course during the route

planning. At sea, deviations of the forecasted weather to on-site wave measurements indicate

if forecast updates and recalculations for optimum routes are necessary. Wave monitoring and

real-time seakeeping simulation serve to identify hazardous situations shortly before they

arise and enable instantaneous countermeasures. In this way monitoring and routing will

merge to support a prudent seamanship.

GODSPEED !

Pareto Optimum Ship Routing, References

121

8 References

ADOPT (2005), Advanced Decision Support System For Ship Design, Operation And

Training, EU 6th Framework Program, April 2005 – Sept. 2008, http://adopt.rtdproject.net

BENFORD, H. (1965), Fundamentals of Ship Design Economics, The University of

Michigan, Department of Naval Architecture an Marine Engineering, Ann Arbor, Michigan

BIDLOT, J.-R.; HERSBACH, H.; RICHARDSON, D.; SAETRA, Ø. (2002), Effects of

observation errors on the statistics for ensemble spread and reliability, Technical

Memorandum 393, Research Department of the ECMWF, Reading, UK

BIRK, L.; CLAUSS, G.F. (2001), Automated Hull Optimization of Offshore Structures based

on Rational Seakeeping Criteria, Proceedings of the 11th International Offshore and Polar

Engineering Conference, Stavanger, Norway, ISBN 1-880653-51-6

BIRK, L.; HARRIES, S. (2003), OPTIMISTIC – Optimization in Marine Design,

39th WEGEMT Summer School, Mensch & Buch Verlag, ISBN 3-89820-514-2

BOESE, P. (1970), Eine einfache Methode zur Berechnung der Widerstandserhöhung eines

Schiffes im Seegang, IFS-Bericht Nr. 258, Institut für Schiffbau der Universität Hamburg,

Germany

BRETSCHNEIDER, C.L. (1957), Review of Practical Methods for Observing and

Forecasting Ocean Waves by Means of Wave Spectra and Statistics, Trans. American

Geophysical Union, Vol. 38, No. 2, U.S. Navy Hydrographic Office, Publ. 603

CHRISTIANSEN, M.; FAGERHOLT, K.; RONEN, D. (2004), Ship Routing and Scheduling:

Status and Perspectives, Transportation Science, Vol. 38, pp. 1-18, ISSN 0041-1655

CLAUSS, G.; LEHMANN, E.; ÖSTERGARD, C. (1988), Meerestechnische Konstruktionen,

Springer-Verlag, ISBN 3-540-18964-5, ISBN 0-387-18964-5

CLAUSS, G.; LEHMANN, E.; ÖSTERGARD, C. (1994), Offshore Structures, Strength and

Safety for Structural Design, Vol.II, Springer-Verlag, ISBN 3-540-19770-2

CLAUSS, G.; RIEKERT, T. (1999), Stochastische Analyse meerestechnischer Systeme,

lecture notes, department for offshore technology, Technische Universität Berlin, Germany

CLAUSS, G.F.; HENNIG, J.; CRAMER, H.; KRUEGER, S. (2003), Development of Safer

Ships by Deterministic Analysis of Extreme Roll Motions in Harsh Seas, 22nd international

Conference on Offshore Mechanics and Arctic Engineering, Proceedings of OMAE’2003,

Mexico

CLAUSS, G.F.; HENNIG, J. (2004), Deterministic Analysis of Extreme Roll Motions and

Subsequent Evaluation of Capsizing Risk, International Shipbuilding Progress, Vol. 51, Issue

2/3, Stability

Pareto Optimum Ship Routing, References

122

CLAUSS, G.F. (2008), The Taming of the Shrew – Tailoring Freak Wave Sequences for Sea

keeping Tests, 29th Georg Weinblum Memorial Lecture, Journal of Ship Research, Society of

Naval Architects and Marine Engineers, to be published in September 2008

CONOLLY, J.E. (1974), Standards of Good Sea keeping for Destroyers and Frigates in Head

Seas, International Symposium on the Dynamics of Marine Vehicles and Structures in Waves,

No. 8, London, UK

CRAMER, H.; REICHERT, K.; HESSNER, K.; HENNIG, J.; CLAUSS, G.F. (2004),

Seakeeping Simulations and Seaway Models and Parameters Supporting Ship Design and

Operation, 9th International Symposium on Practical Design of Ships and other Floating

Structures, PRADS 2004, Germany

EMMECHE, C. (1994), Garden in the Machine. The Emerging Science of Artificial Life,

Princeton University Press

ES.TEC.O s.r.l. (1999), modeFRONTIER, Multi-Objective Design Environment, User

Manual, Engin Soft Tecnologie per l'Ottimizzazione, Italy

FARIN, G. (2001), Curves and surfaces for CAGD, a practical guide, Morgan Kaufmann

Publishers Inc., ISBN 1-55860-737-4, USA

FRANK, W. (1967), Oscillation of Cylinders in or below the Free Surface of Deep Fluids,

Naval Ship Research and Development Center, Technical Report 2375,Washington DC, USA

GERMANISCHER LLOYD (2002), annual report, Hamburg,

http://www.gl-group.com/magazines/ann_02/cover02.pdf

GERRITSMA, J.; BEUKELMAN, W. (1972), Analysis of the resistance increase in waves of

a fast cargo ship, International Shipbuilding Progress, Vol. 19, No. 217

HAGIWARA, H.; FUKUDA, H.; SUGAI, K.; KUSAKA, Y. (1999), A Study on Weather

Routing of High Speed Ships, FAST´99, 5th International Conference on Fast Sea

Transportation, Seattle

HAGIWARA, H.; SPAANS, J.A. (1987), Computation of Optimum Routes for Ship Weather

Routing, International Symposium on Weather Routing, Tokyo

HAPEL, K.-H. (1990), Festigkeitsanalyse dynamisch beanspruchter Offshore-

Konstruktionen, Vieweg & Sohn Verlagsgesellschaft mbH, Technische Universität Berlin,

Germany

HARRIES, S. (1998), Parametric Design and Hydrodynamic Optimization off Ship Hull

Forms, Dissertation at TUB, Institut für Schiffs- und Meerestechnik, ISBN 3-933346-24-X

HASSELMANN, K.; ROSS, D.B.; MÜLLER, P.; SELL, W. (1975), A Parametric Wave

Prediction Model, Journal of Physical Oceanography, Volume 6

Pareto Optimum Ship Routing, References

123

HINNENTHAL, J.; SAETRA, Ø. (2005), Robust Pareto – Optimal Routing of Ships utilizing

Ensemble Weather Forecasts, Proceedings of the 4th International Conference on Computer

and IT Applications in the Maritime Industries, COMPIT’05, Hamburg, Germany,

ISBN 3-00-014981-3

HOFFSCHILDT, M.; BIDLOT, J-R.; HANSEN, B.; JANSSEN, P.A.E.M. (1999), Potential

benefit of ensemble forecast for ship routing, ECMWF Techn. Memo. 287, Internal Report,

ECMWF, Reading/UK

HOLLAND, J. (1975), Adaptation in natural and artificial systems, University of Michigan

Press, reprinted in 1992 by MIT Press. Cambridge

HOLTROP, J. (1978), Statistical Data for the Extrapolation of Model Performance Tests,

International Shipbuilding Progress, Vol. 25, Rotterdam

HOLTROP, J.; MENNEN, G.G.J. (1984), An approximate power prediction method,

International Shipbuilding Progress, Vol. 31

IMO (1988), Code on Intact Stability for All Types of Ships, A.749(18), amended by

resolution MSC.75(69)

IMO (1995), Guidance to the Master for Avoiding Dangerous Situations in Following and

Quartering Seas, MSC/Circ.707

IMO (2000 and 2003), Casualty Statistics and Investigations, International Maritime

Organization, Ref. T1/2.02, FSI3/Circ.3 and 6, UK

IMO (2007), Revised Guidance to the Master for Avoiding Dangerous Situations in Adverse

Weather and Sea Conditions, MSC/Circ.1228

ITTC (1978), Performance prediction method for single screw ships, 15th Int. Towing Tank

Conf., The Hague

JOURNÉE, J.M.J. (2000), Theoretical Manual of SEAWAY, Rel. 4.18, Report 1216, Delft

University of Technology, The Netherlands

JOURNÉE, J.M.J. (2001), User Manual of SEAWAY, Rel. 4.19, Report 1212a, Delft

University of Technology, The Netherlands

JOURNÉE, J.M.J.; ADEGEEST, L.J.M. (2003), Theoretical Manual of SEAWAY, “SEAWAY

for Windows”, Report 1370, Delft University of Technology and AMARCON Advanced

Maritime Consulting, The Netherlands

JOURNÉE, J.M.J.; MASSIE, W.W. (2001), Offshore Hydromechanics, First Edition, Delft

University of Technology, http://www.shipmotions.nl

JOURNÉE, J.M.J.; MEIJERS, J.H.C. (1980), Ship Routing for Optimum Performance,

Conference on Operation of Ships in Rough Weather, Transactions IME, London

KRAPPINGER, O. (1983), Beurteilung des Einflusses von Schiffsbewegungen auf die

Besatzung, Schiffstechnik Bd. 30, Hamburgische Schiffbau Versuchsanstalt

Pareto Optimum Ship Routing, References

124

KRUEGER, S.; HINRICHS, R.; KLUVE, F.; BILLERBECK, H. (2006), Towards the

Development of Dynamic Stability Criteria, TUHH 3-14, Hamburg, also available as

Proposal of a probabilistic intact stability criterion, Revision of the Intact stability Code,

IMO, SLF 49/5/2

LAGARIAS, J. C.; REEDS, J. A.; WRIGHT, M. H.; WRIGHT, P. E. (1998), Convergence

Properties of the Nelder-Mead Simplex Method in low Dimensions, SIAM Journal on

Optimization, Vol. 9, available at wikipedia.org

LASSE (2006 –2008), scientific project funded by the German Ministry, http://www.naoe.tu-

berlin.de/MT/Forschung/forschung.html,

LEWIS, E.V. (Ed.) (1998), Principles of Naval Architecture - Motions in Waves and

Controllability, Vol. III, SNAME, Jersey City

LLOYD, A.R.J.M. (Ed.) (1998), Sea keeping, Ship Behavior in Rough Weather, RPM

Reprographics, Chichester, Sussex, United Kingdom, ISBN 0 9532634 0 1

MAN B&W (2000), K90MC MK6 Project Guide Two-stroke Engines, 5th Edition,

www.manbw.dk

MATLAB (2002), version 6.5.0, release 13, http://www.mathworks.com/

NELDER, J.A.; MEAD, R.A. (1965), A Simplex Method for Function Minimization.

Computer Journal Vol.7, Oxford University Press, UK

NORDFORSK (1987), Assessment of Ship Performance in a Seaway, Results of a Nordic

Co-Operative Project on Seakeeping Performance of Ships, Denmark, ISBN 87-982637-1-4

NOWACKI, H.; BRUSIS, F.; SWIFT, P. (1970), Tanker preliminary design – An

optimization problem with constraints, Transactions Vol. 78, SNAME

NOWACKI, H.; LESSENICH, J. (1976), Synthesis models for preliminary ship design,

Jahrbuch der Schiffbautechnischen Gesellschaft STG, Vol. 70

OCHI, M.K. (1964), Extreme Behavior of Ship in Rough Seas – Slamming and Shipping of

Green Water, Transactions Volume 72, The Society of Naval Architects and Marine

Engineers, New York

PAPANIKOLAOU, A. (1977), Zum nichtlinearen Problem eines vertikal schwingenden

zylindrischen Körpers beliebiger Form, Techische Universität Berlin, Doctoral Thesis, D83

PAPANIKOLAOU, A. (1982), Über Analytische und Numerische Methoden zur Behandlung

des Irregularitätenproblems bei der Anwendung von Integralgleichungsmethoden, Technische

Universität Berlin, TUB/ISM Bericht Nr.82/10

PAPANIKOLAOU, A.; GRATSOS, G.; BOULOUGOURIS, E.; ELIOPOULOU, E. (2000),

Operational Measures for Avoiding Dangerous Situations in Extreme Weather Conditions,

7th International Conference on the Stability of Ships and Ocean Vehicles STAB 2000,

Australia

Pareto Optimum Ship Routing, References

125

PIERSON, W.J.; MOSKOWITZ, L. (1964), A Proposed Spectral Form for Fully Developed

Wind Seas Based on the Similarity Theory of S.A. Kitaigorodskii, Journal of Geophysical

Research, Vol. 69

PINGREE, B.J.W. (1988), A review of human performance in a ship motion environment,

Warship Technology No. 4, Institute of Naval Medicine, Alverstoke Hampshire

RATHJE, H.; BEIERSDORF, C. (2004), Decision support for Container Ships in Heavy Seas

– Shipboard Routing Assistance, International Symposium Information on Ships ISIS 2004,

German Institute of Navigation DGON, Hamburg

RECHENBERG, I. (1973), Evolutionsstrategie, Optimierung technischer Systeme nach

Prinzipien der biologischen Evolution, Reihe Problemata Vol. 15, Friedrich Frommann

Verlag, Stuttgart

SAETRA, Ø. (2004), Ensemble Shiprouting, ECMWF Technical Memorandum 435,

Reading/UK

SAMSUNG SHIPBUILDING & HEAVY INDUSTRIES (1991), 4,400 TEU Container

Vessel, Trim and Stability Calc. (M/V Hannover Express), Koje Island, South Korea

SCHÖNEBURG, E.; HEINZMANN, F.; FEDDERSEN, S. (1994), Genetische Algorithmen

und Evolutionsstrategien – Eine Einführung in Theorie und Praxis der simulierten Evolution,

Addison-Wesley, ISBN 3-89319-493-2

SEAROUTES (2003), Routing Parameters and Priorities, Advanced Decision Support for

Shiprouting based on Full-scale Ship-specific Responses as well as Improved Sea and

Weather Forecasts Including Synoptic High Precision and Real-time Satellite Data, EU 5th

Framework Program, Contract G3RD-CT-2000-00309

SHIN, Y.S.; BELENKY, V.L.; PAULLING, J.R.; WEEMS, K.M.; LIN, W.M. (2005),

Criteria for Parametric Roll of Large Containerships in Longitudinal Seas, ABS, a

contribution to the SNAME 13th ad-hoc panel parametric rolling

SPAANS, J.A.; STOTER, P. (2000), Shipboard Weather Routing, International Symposium

Information on Ships ISIS 2000, German Institute of Navigation DGON, Wilhelmshaven

TREBY, A. (2002), Optimal Weather Routing Using Ensemble Weather Forecasts,

37th annual conference, Operational Research Society of New Zealand

TSUJIMOTO, M.; TANIZAWA, K. (2006), Development of a Weather Adaptive Navigation

System Considering Ship Performance In Actual Seas, 25th International Conference on

Offshore Mechanics and Arctic Engineering, Proceedings of OMAE 2006, Germany

YASAKI, A. (1962), Design diagrams for modern four, five, six and seven-bladed propellers

developed in Japan, 4th Naval Hydrodynamics Symposium, Office of Naval Research

Pareto Optimum Ship Routing, Appendix 1

126

Appendix 1 Stochastic evaluation

The stochastic evaluation of the ship motion uses equations (4.9) and (4.10) to assess the 0th

moments of the response spectra for velocity and acceleration:

Smm 2

& and S

Smm 4

&& . (4.9) and (4.10)

In the following a derivation of these equations will be given in detail. For the time being, the

random process of wave elevation of an irregular sea state

()

is focused. However, it is

shown that these considerations are valid for ship motions, too. The derivation bases on the

WIENER-CHINTSCHIN16-relations. They describe the relation of a spectrum (autospectrum)

and the appropriate autocorrelation function of a random stationary process:

the FOURIER transform of a the autocorrelation function gives the spectrum

() ()()

∫

+∞

∞−

−⋅=

τωττ

ωζζζζ

diRS exp

1, (A1.1)

the inverse FOURIER transform of a spectrum gives the autocorrelation function

() ( ) ( )

∫

+∞

∞−

⋅=

ωωτωτ ζζζζ

diSR exp , (A1.2)

cp. CLAUSS and RIEKERT (1999).

The relation between autospectrum and the wave spectrum that is used for the statistical

evaluation of ship motions is illustrated in Fig. 96 and given by equation (A1.3). The

autospectrum occurs because negative frequencies are permitted within the FOURIER

transform. In contrast to the mathematical evaluation, no negative frequencies exist from the

physical point of view. For this reason the asymmetric representation of the spectrum is

employed to assess e.g. wave elevation or ship motions. The transformation given by equation

(A1.3) follows the requirement that the area m0 under both curves has to be the same as it

represents a measure for the energy of the considered process.

Fig. 96: Symmetric autospectrum and asymmetric spectrum

of a random stationary process

16 Alexander Jakowlewitsch Chintschin, mathematician, 1894 – 1959, Russia

Norbert Wiener, mathematician, 1894 – 1964, North America

Pareto Optimum Ship Routing, Appendix 1

127

() ()

ζζζ

SS 2= for: ∞<<

0 (A1.3)

Another definition of the autocorrelation function is given by:

() () ( )

∫

−

∞→ +⋅=

lim

Tdttt

τζζτ

ζζ

. (A1.4)

Combining equations (A1.2) and (A1.4) result in equation (A1.5):

() () ( ) () ( )

∫∫ +∞

∞−

−

∞→ ⋅=+⋅=

ωωτωτζζτ ζζζζ

diSdttt

Texp

lim

. (A1.5)

Just as the surface elevation of irregular waves

, surface velocities

& and accelerations

represent random distributed GAUSSian processes. Thus the relation of the autocorrelation

function and the spectral density for the velocity of the process is given by:

() () ( ) () ( )

∫∫ +∞

∞−

−

∞→ ⋅=+⋅=

ωωτωτζζτ ζζζζ

diSdttt

Texp

lim

&&&&

&& . (A1.6)

Equations (A1.5) and (A1.6) are used by HAPEL (1990) to show the validity of equations

(4.9) and (4.10). Building the derivative of

()

ζζ

R with respect to

leads to:

() ( ) () ( )

∫∫ +∞

∞−

−

∞→ ⋅=+⋅=

ωωτωωτζζ

τζζ

ζζ

diSidttt

dR T

Texp

lim

&. (A1.7)

The substitution

−= tt yields:

()() () ( )

∫∫ +∞

∞−

−

∞→ ⋅=⋅−=

ωωτωωζτζ

τζζ

ζζ

diSidttt

dR T

Texp

lim

&. (A1.8)

A further derivative of

()

ζζ

R with respect to

and the substitution

+= tt leads to:

()() () ( )

∫∫ +∞

∞−

−

∞→ ⋅−=⋅−−=

ωωτωωζτζ

τζζ

ζζ

diSdttt

dR T

Texp

lim 2

2&& (A1.9)

and further to

() ( ) () ( )

∫∫ +∞

∞−

−

∞→ ⋅−=+⋅−=

ωωτωωτζζ

τζζ

ζζ

diSdttt

dR T

Texp

lim 2

2&& . (A1.10)

A comparison of equations (A1.6) and (A1.10) shows that following equations are valid:

()

τζζ

ζζ

R−=

&& , (A1.11)

() ()

ωωω ζζ

ζζ

SS 2

&& . (A1.12)

Repeating the process above yields the respective equation for the acceleration of the wave

elevation

&& :

() ()

ωωω ζζ

ζζ

SS 4

&&&& . (A1.13)

Equations (A1.12) and (A1.13) are also valid for asymmetric spectra like employed in the

study because the area under each spectrum, i.e. for motion, velocity, or acceleration, remains

the same and the transformation (A1.3) produces no shift of the spectral density with regard to

the absolute value of the frequency. For this reason the following equations are employed

within practical application:

() ()

ωωω ζ

SS 2

&, (A1.14)

() ()

ωωω ζ

SS 4

&& . (A1.15)

Consequently the 2nd and the 4th moment of the spectrum describing the surface

elevation

equal the 0th moments of the velocity- and acceleration spectrum:

Pareto Optimum Ship Routing, Appendix 1

128

() ()

ζζ

ωωωωω

&& 0

2mdSdSm === ∫∫ ∞∞

(A1.16)

() ()

ζζ

ωωωωω

&&&& 0

4mdSdSm === ∫∫ ∞∞

(A1.17)

CLAUSS et. al. (1994) show that linear functions of GAUSSian random variables like the

surface elevation of waves are themselves GAUSSian random variables. This applies to the

determination of ship responses by means of equation (4.1), i.e. the ship motion represents a

GAUSSian random process too.

() ()

s⋅=

(4.1)

Therefore the considerations made above are valid for the motion spectra, too.

The 0th moments of the velocity- and the acceleration spectrum of the ship equal the 2nd and

the 4th moment of the response spectrum for the ship motion:

() () ()

SS mdS

dSdSm 2

0∫∫∫ ∞∞∞

=⋅===

ωω

ωωωωωω ζ

&& (A1.18) or (4.9)

() () ()

SS mdS

dSdSm 4

0∫∫∫ ∞∞∞

=⋅===

ωω

ωωωωωω ζ

&&&& (A1.19) or (4.10)

Pareto Optimum Ship Routing, Appendix 2

129

Appendix 2 Irregular frequencies

In case of strip-theory methods applied to seakeeping problems, the determination of

hydrodynamic forces and moments requires the solution of a series of 2D potential theory

boundary value problems (BVP) for a number of cross-sections that are representative for the

three-dimensional hull shape of the ship. The solution of these BVPs leads to the

determination of the velocity potential and hence to hydrodynamic coefficients (added mass

and damping coefficients) which are dependent on the shape of the relevant cross-section, the

frequency of oscillation, and the direction of the considered motion, i.e. the sway-, heave-, or

roll motion. Hydrodynamic forces and moments, acting on the cross-section, the 2D-strip, are

derived by integration of hydrodynamic pressure deduced from the determined velocity

potentials by application of Bernoulli’s equation.

When solving the potential theory BVP by numerical solution for formulated integral

equations (i) with a distribution of pulsating sources along the wetted part of a cross-section

and (ii) for surface piercing cross-sections, numerical results are disturbed at particular

frequencies, the so called irregular frequencies, cp. JOURNÉE and MASSIE (2001). Fig. 97

shows non-dimensionalized mass- and damping coefficients for the heave motion of the S175,

presented in section 6.3. Mass- and damping coefficients are given with respect to the

frequency of the vertical motion and the position of the cross-section, where x/Lpp = 1

represents the bow, 0 the stern, respectively. The unrealistic peak at a frequency of about

2.5rad/s, in particular visible in the diagram for the damping coefficient, is caused by the

occurrence of an irregular frequency.

Fig. 97: 2D values for mass- M33 and damping N33 coefficients, S175

Fig. 98 depicts mass- and damping coefficients for the affected cross-section near the stern.

It is evident that both the mass- and damping coefficients are disturbed.

Pareto Optimum Ship Routing, Appendix 2

130

Fig. 98: 2D values for mass- M33 and damping N33 coefficients, section 5

As explained above, the determination of the hydrodynamic coefficients by means of

potential theory presumes the solution of a boundary value problem that can be solved by

different methods appropriate to the required accuracy. In case that a traditional conformal

mapping method is not able to satisfactorily resemble a cross-section shape, which is often the

case for realistic ship-like sections, it is useful to apply a 2D singularity distribution method

for the solution of the relevant 2D BVP, e.g. FRANK’s pulsating source method, see

FRANK (1967). In this case the calculation of the flow field around the cross-section is based

on the solution of a GREEN’s function integral equation. Numerically the solution is provided

by means of a discretized representation of the cross-section, i.e. the boundary conditions are

fulfilled at a number of collocation points on the cross-section contour. For the appropriate

discretized source distribution, given by a discretized representation of the GREEN’s function

integral, the unknown strengths of the pulsating sources have to be determined by solution of

a linear algebraic set of equations.

The solution of this problem consists of an internal and external solution describing the flow

field inside and outside the cross-section. The external solution is used for the determination

of the fluid pressure around and along the cross-section, whose integration leads to the

hydrodynamic coefficients. Dependent on the cross-sectional shape, the breadth/draught ratio

B/T and the cross-sectional area coefficient Am, the internal solution disposes eigenfrequencies

corresponding to sloshing effects (resonances) of the imaginary internal flow. In this case no

realistic solution for the unknown source strengths determining the external velocity potential

exists and therefore the external solution is disturbed (because of the so-called

duality problem of potential theory, see PAPANIKOLAOU (1977)). Definitely there is no

physical reason for the sudden change of mass- and damping coefficients at particular

frequencies, as visible in Fig. 98. Therefore the numerical results for the affected frequency

band around the irregular frequencies cannot be used.

From the numerical point of view it can be shown that this problem is caused by a zero (or

sudden decrease) of the determinant of the integral equation that has to be solved for the

unknown source strengths. This also simply explains why mass- and damping coefficients

determined by conformal mapping methods or for fully submerged cross-sections do not show

a distortion by irregular frequencies. The former do not use the GREEN’s function integral

equation method and the latter always disposes a unique solution, as the internal fluid

problem has no eigenfrequencies. Anyhow it is possible to avoid the occurrence of irregular

frequencies or at least to shift them out of the considered practical frequency range, also for

singularity panel methods.

Pareto Optimum Ship Routing, Appendix 2

131

PAPANIKOLAOU (1982) presents a series of different techniques to eliminate the

disturbance by irregular frequencies regarding the 2D problem. They are divided into three

classes:

• Analytical methods make use of a modified GREEN’s function. Additional source

terms in the GREEN’s function modify the FREDHOLM’s17 determinant that finally

serves to shift the eigenfrequencies, i.e. the frequencies where the determinant

becomes zero are shifted to higher frequencies that are out of the practical range for

the determination of ship motions.

• Analytic-numerical methods use a numerical lid-inside deck, i.e. a discretization of the

free surface inside the cross-section contour. Alternative boundary conditions for the

internal flow serve to control the internal solution and to eliminate eigenfrequencies in

the practical range of frequencies. In a way this method forces the internal solution to

behave like the one of a fully submerged cross-section. It finally results in a

modification of the FREDHOLM’s determinant, too.

• Numerical methods pre-assess the position of irregular frequencies by (i) identification

of typical constellations of the breadth/draught ratio B/T and the cross-sectional area

coefficient Am or (ii) they observe the value of the FREDHOLM’s determinant for this

purpose. For the identified irregular frequencies the values of the flow potentials are

interpolated, e.g. by a spline-interpolation using the surrounding values as

interpolation values.

All three sets of methods have in common that they require access to the source code of the

applied seakeeping program. In case of the code used within this thesis it was not possible to

implement the above techniques and a simplified method to deal with this problem was

adopted, as described in section 4.7.

17 Erik Ivar Fredholm, mathematician, 1866 – 1927, Sweden

Pareto Optimum Ship Routing, Appendix 3

132

Appendix 3 Root mean square values

The abbreviation RMS represents the root mean square value of a stochastic process. Within

the present study the root mean square values of two different but associated process factors

are used:

• the RMS of the wave surface displacement or a ship motion,

• the RMS of the wave amplitude or a ship motion amplitude.

The following nomenclature is used to describe the relation between these terms and to avoid

mistake:

symbol name unit

a parameter of a harmonic or a random process [ dep. ]

ã amplitude of a harmonic function [ dep. ]

a1/3 significant wave amplitude [ m ]

arms root mean square of the wave amplitude [ m ]

H1/3 significant wave height / double amplitude [ m ]

Hrms root mean square of the wave height [ m ]

m0 0th order moment of a wave power spectrum [ m2 ]

RMS root mean square value [ dep. ]

standard deviation of the wave amplitude [ m ]

standard deviation of the wave surface displacement [ m ]

T time, time period of a wave record [ s ]

ζ (

)

wave surface displacement [ m ]

rms

root mean square of the wave surface displacement [ m ]

The equations (A3.1) to (A3.3) are used as starting point of this explanation:

Equation (A3.1) represents the relation of root mean square of the wave surface displacement

of an irregular sea state, the standard deviation of the wave surface displacement, and the

0th moment of the wave spectrum.

()

1mtdt

rms === ∫

σζζ

(A3.1)

Equation (A3.2) denotes the relation of significant wave height (significant double amplitude)

and 0th moment of the power spectrum.

03/1 4mH ⋅= (A3.2)

Equation (A3.3) gives the relation of significant wave height and root mean square of the

wave heights, in other words the relation of significant double amplitudes and root mean

square of double amplitudes, respectively.

rms

HH ⋅= 2

3/1 (A3.3)

Pareto Optimum Ship Routing, Appendix 3

133

The relation of root mean square of the wave height and 0th moment of the power spectrum is

obtained from equating (A3.2) and (A3.3) with

22 mHrms ⋅= . (A3.4)

Equation (A3.1) shows the relation of standard deviation of the wave surface displacement

and 0th moment of the power spectrum 0

. In the same way equation (A3.4) can be

used to obtain the relation of standard deviation of the wave amplitude and power spectrum

, (A3.5)

with:

armsrms aH

⋅=⋅= 22,

the relation of root mean square of the wave height (double amplitude) and standard deviation

(root mean square) of the wave amplitude.

Equating the 0th moments of the power spectrum from the equations (A3.1) and (A3.5)

provides the relation of the standard deviations of wave amplitudes and wave surface

displacement

σσ

⋅= 2

a, (A3.6)

or rather the relation of root mean square of the wave amplitudes and of the wave surface

displacement

rmsrms

⋅= 2. (A3.7)

Setting up equation (A3.3) for amplitudes instead of double amplitudes provides

rms

aa ⋅= 2

3/1 . (A3.3a)

The relation of significant wave amplitude and root mean square of the wave surface

displacement is obtained from inserting equation (A3.7) in (A3.3a).

rmsrms

ζζ

⋅=⋅⋅= 222

3/1 (A3.8)

These considerations regarding the wave surface displacement and wave amplitudes are

transferable to the ship motions and amplitudes of motion obtained by linear theory.

Regarding equation (6.3) the terms of this equation have to be interpreted as follows:

RMStda

a⋅=⋅=⋅= ∫2

1, (6.3)

s significant amplitude of a motion,

2m⋅ 0

22 m⋅⋅ = ⋅2 (root mean square of the motion amplitude),

∫

⋅

tda

2 ⋅2 (root mean square of the motion),

RMS⋅2 ditto.

The thresholds posed within NORDFORSK (1987) given in Tab. 13 and Tab. 14 represent

root mean square values of vertical accelerations (not amplitudes of acceleration). Therefore

the conversion from RMS to significant amplitudes has to follow equation (6.3).