Ship routing optimisation based on forecasted weather data and considering safety criteria

Abstract

A weather routing tool is presented based on forecasted weather data along the route and considering safety aspects. The tool aims to determine the optimal path for the minimisation of the fuel oil consumption, ensuring a safe passage. It is developed in MATLAB and considers detailed ship characteristics. Specifically, ship's motions and fuel oil consumption of the main engine during a potential path are estimated. For the latter, a physics-based model for a specific vessel is developed where tools of different level of detail are utilised to calculate the various resistance components. A speed management strategy along the route is specified as well as safety criteria representing acceptable limits of ship's responses. When the set criteria and constraints have been set, a genetic algorithm is used to find the optimal route by means of ship's heading or by considering both heading and ship's power settings as variables to minimise the fuel oil consumption. The search space of the algorithm lies within a predefined envelop, but still the evolutionary optimisation approach used has no pre-assigned values to any possible candidate waypoint.

1. Introduction

Global climate change has already had observable effects on the environment and poses serious threat to our future. In the past years and for the maritime sector, a regulatory framework was introduced concerning the energy efficiency of new ships and recently for existing ships. Furthermore, International Maritime Organization (IMO) has already set a strategy to reduce greenhouse gases related to ships’ operation based on short-, medium- and long-term measures (IMO, 2018). As many studies show, fuel cost constitutes by far the greatest part (up to 70%) of a ship's operational cost and, as a result, the reduction of fuel oil consumption per voyage is a priority for all shipping companies. In recent years, many measures have been proposed for the reduction of greenhouse gases. Two main categories of measures exist. The first one is related with ship design characteristics aiming at the improvement of propulsion characteristics, the reduction of the total resistance (hull and propeller), the improvement of the machinery power performance and the usage of alternative energy sources. The second one regards measures concerning operational efficiency. An overview of the reduction potential achieved by the candidate measures has been presented in Brynolf et al. (2016) and Bouman et al. (2017). According to Bouman et al. (2017), measures that control operational efficiency can lead to great reduction of CO₂ emissions estimated up to 48–60%. The most efficient short-term measures identified were route planning, which is employed to avoid severe weather conditions, and the ship's speed optimisation to reduce the main engine's fuel oil consumption. In addition, optimising a ship's route and speed based on weather conditions and seakeeping behaviour can improve scheduled trade and, of course, the safety of the ship, her crew and the cargo. However, the decision-making process is not always based on data and/or models that reflect the ship's performance under the prevailing weather conditions, resulting to low prediction potential.

In this paper, a weather routing system is presented. The main aim is the identification of an optimised route, in terms of fuel oil consumption, from a set of generated feasible routes. A physics-based model is developed for the calculation of the fuel consumption of the main engine considering the wave, wind and sea currents conditions along the route, while seakeeping criteria are applied to ensure a safe passage. A genetic algorithm is employed to solve the optimisation problem in terms of optimal heading selection and optimal set of speeds between every two waypoints. Section 2 provides a brief review of works covering several aspects of this study. Section 3 presents in detail the developed methodology, and Section 4 presents case studies, demonstrating the optimised results in comparison with the orthodromic and the loxodromic paths. Finally, Section 5 highlights the main findings of the study, and future developments are discussed.

2. Literature review

While weather routing has been a well discussed topic in past decades and industry applications have been developed, there are several aspects of the problem examined by later studies. There are many methods for generating the shortest path between two points. The utilisation of the Dijkstra algorithm (e.g. Sen and Padhy, 2015) and the A-star algorithm (e.g. Bentin et al., 2016) are well-known examples, while Joo et al. (2012) used cell-based methods by firstly considering a constant ship's speed along the route, and then by adjusting the speed properly on the generated route. There also studies proposing the simultaneous optimisation of route and speed, such as Zaccone et al. (2018), where a dynamic programming approach was followed using available forecast maps. Wang et al. (2020) presented a new read-coded genetic algorithm which mainly focused on the determination of the minimum travelling time path while considering involuntary speed loss due to wind and waves factors.

The isochrone method was used by James (1957) to estimate time-efficiently an optimal route corresponding to specific weather conditions. This method had been quite popular and was the basis for several studies that proposed improved versions, such as Szlapczynska and Smierzchalski (2007), where the revolutionary approach was used to generate routes bereft of land crossings. Moreover, Lin et al. (2013) employed an algorithm using a 3D modified isochrone method (3DMI), while Lee et al. (2018) proposed a similar method where speed optimisation is performed and fuel oil consumption is estimated, considering the weather conditions in comparison with theoretical methods that ignore weather impact. In addition to well-known methods and their improvements, several new methods have been presented. According to Szlapczynski (2006), a new route-finding technique was introduced based on raster planes. This algorithm seeks the optimal route while trying to minimise the number of course changes each time, giving advantage to the method against the original one previously presented by Lee (1961) and its improved method introduced by Chang et al. (2003). Turning penalties are introduced, meaning that the shortest path may not be the optimal for both the original and the improved method. However, in many studies, the focus is not only on the proposed optimisation method but also on a ship's performance and safety. Lu et al. (2015) developed a ship's resistance model to predict ship's performance at various loadings and sea states. The model is based on Kwon's method in conjunction with noon reports and sea trial data for a specific ship type. Vettor and Soares (2016) presented a weather routing code that can optimise a ship's route in terms of minimum fuel oil consumption, while ship's exposure to risks due to severe weather conditions can be avoided (e.g. slamming event, high vertical accelerations at the bridge, etc.). The hydrodynamic calculations, concerning response amplitude operators (RAOs) of the motions and added resistance, are carried out by using a strip theory-based code, while calm water resistance is estimated according to Holtrop and Mennen's method.

Hinnenthal and Harries (2004) selected the most advantageous route on the basis of hydrodynamic simulations, where a strip theory code was utilised to calculate in advance the added wave resistance and the ship motion response in the frequency domain. Pachero and Soares (2007) presented the use of geographic information system (GIS) for the optimal ship's route selection. The best sailing route is selected by making small deviations from the shortest path, which is the orthodrome. A wave forecast model is implemented, and the system can predict ship's responses (roll, pitch, heave) at various sea conditions. Bentin et al. (2016) examined the effect of a wind-assisted ship propulsion system (WASP) within a weather-routing optimisation tool and investigated the fuel oil reduction potential when using such systems in contrast to a conventional diesel propulsion system.

Lamprinidis and Belibassakis (2017) presented a direct method for modelling ship weather routing. The downhill simplex method was selected to solve the optimisation problem, while for the calculation of the added resistance and the vertical response at the stern, the radiated energy method in conjunction with strip theory was used. Kobayashi et al. (2015) proposed a new advanced weather routing optimisation method . It targeted the minimisation of fuel oil consumption and the energy efficiency operational indicator (EEOI). An entire manoeuvring simulation was performed along transoceanic trajectories by solving differential equations concerning ship's motions.

In addition, data-driven models have also been implemented for the improvement of ship's operational efficiency, where they take advantage of the vast amount of operational data available. Karagiannidis and Themelis (2021) utilised high-frequency operational data and test several pre-processing techniques to examine their effect on the accuracy of data-driven models for the prediction of fuel oil consumption using a feed-forward neural network. Gkerekos and Lazakis (2020) presented a heuristic framework for ship weather routing, where the main engine's FOC is predicted by a data-driven model. Data-driven energy efficiency analysis in shipping was also presented by Zhang et al. (2019), focusing on ships sailing the Arctic waters. Zis et al. (2020) provided a taxonomy by reviewing several ship weather routing studies, while Wen et al. (2017) examined the problem by considering also other parameters, such as the fuel price and freight rate, but this type of analysis is beyond the scope of the current study.

3. Methodology

The developed methodology is presented in Figure 1. As mentioned in the introduction, the basic idea lies on the evaluation of a cost function that represents the main engine's fuel oil consumption for several candidate routes considering the encountered weather conditions and safety restrictions. Thus, as a first step, the ports of departure and arrival are defined as the starting and ending points, respectively, as well as a number of waypoints. Next, a predefined number of random routes are generated based on functions developed in the MATLAB environment. Then, the acquisition of all the necessary seakeeping and weather data and their pre-processing to save time during optimisation is analysed. For the estimation of the fuel oil consumption of the main engine, the total resistance is calculated, while the main engine's and the propeller's characteristics are utilised. The voluntary and involuntary speed loss due to weather conditions is then discussed and the handling of speed deviations is presented. A detailed analysis of the previously mentioned steps follows.

Figure 1. Weather routing application flowchart

3.1 Routes generation

The presented weather-routing tool has been developed in the MATLAB environment employing functions and toolboxes, such as the mapping toolbox which provides functions for analysing geographic data and creating map displays. In Zaraphonitis et al. (2021), a brief analysis of the developed tool has been presented. The vector format is used for the loaded geographic data, which form sequentially ordered pairs of latitude and longitude coordinates. The vector map can be displayed using a suitable map projection. For navigation purposes, the Mercator projection is the most suitable one, since on this display the loxodromic paths are presented by straight lines on the created map.

By having defined the map environment, an algorithm has been developed generating random routes and displaying them on the map. First, the starting and the ending point of a route are provided, which correspond to the coordinates of the departure and the arrival ports, respectively. Thereafter, there are two options for the generation of the routes. In the first one and based on the n waypoints, a developed algorithm randomly locates them within a predefined envelope. In the second one, the location of all or some of the waypoints is specified. In both options, the group of all these points creates n + 1 legs, where each one can be separately analysed. For instance, it is possible to divide each leg of the whole route into several equidistant points with respective coordinates. These coordinates are unique and based on the available functions, information for each point, such as whether it lies in the sea or at land, the sea depth at that location and its distance from the nearest coast, can be obtained.

With the developed algorithm, a targeted number of different routes are generated, as shown in Figure 2, having as input the starting and the ending points and the random position of the n waypoints. In addition, constraints can be defined such as not allowing the operation within a specific distance from the coast or specifying a minimum depth along the route. Any route that does not satisfy the constraints is characterised as unfeasible and it is not considered further in the process.

Figure 2. Generation of 10 random routes

3.2 Weather and seakeeping data

3.2.1 Weather data

Weather data, and specifically wave, wind and current data, can be acquired at different temporal and spatial resolution. Temporal resolution refers to the time needed to update the information for the same location. During this time step, the weather conditions are assumed to remain constant, which means that having a low temporal resolution can result in poor estimations considering that weather conditions are time sensitive. On the other hand, a high temporal resolution can lead to a time-consuming process with not a considerable change in the quality of the results. Similarly, spatial resolution refers to the distance between every two locations of the grid where weather data are provided.

In case the wave, wind and current data are not available at the same spatial or temporal resolution, the tool sets as the main grid for the optimisation problem the wave's spatial resolution. The intersection points of the main grid as the ship crosses each ‘grid square’ are identified for the examined route, as shown in Figure 3. For these points, shown as red stars, the wave values are calculated through linear interpolation based on the coordinates of these points and their distance from the four picks forming each ‘grid square’ (black dots). Similarly, wind and current values at the intersection points are calculated through linear interpolation based on the nearest grid square formed by their spatial resolution. Finally, between every two intersection points weather conditions are assumed to remain constant, which is an adequate assumption when a dense grid is used.

Figure 3. Main grid formed from wave data resolution equals to 0 ⋅ 5° × 0 ⋅ 5° superposed to 1° × 1° wind resolution. The yellow line represents a random ship route intersecting the example grid section and the red stars represent the intersection points of the random route with the main grid where the desirable values are calculated through linear interpolation

3.2.2 Seakeeping data

During the assessment of each route, constraints are set based on ship-motion criteria. Therefore, ship motions shall be calculated during the passage of the ship in each route-leg considering the prevailed weather conditions. Safety aspects could range from motion, manoeuvring and structural response, as mentioned in the review. At this stage, ship-motion criteria are introduced aiming at identifying the subroute sections with unsafe ship response or excessive accelerations in critical locations onboard. To implement such criteria, the estimation of the ship's responses in real sea conditions is necessary using seakeeping analysis, based on time or frequency domain models. For the presented application, a frequency-domain model has been used and the respective calculation procedure is presented in the next section. The calculations concerning the ship's motions are performed in advance and are utilised in the main application through appropriate linear interpolation.

3.3 Fuel oil consumption calculation

3.3.1 Calm water and added wave resistance

Calm water resistance is a resistance component corresponding to the absence of the effect of wind and waves. It depends on the vessel's submerged volume, form and speed. Cleanness of the hull is also a main factor that affects the calm water resistance, since its frictional part increases due to algae, microorganisms and other materials that can stick to the hull. Still (or calm) water resistance is a major part of the total resistance and can be calculated by using experimental resistance curves extracted by towing-tank experiments, semi-empirical formulas (e.g. Holtrop and Mennen, 1982) or by computational hydrodynamic calculations (CFD). The operation in waves imposes an added resistance and in severe seas the attainable speed for the same engine power will be decreased. For the application presented in this paper, the added wave resistance is pre-calculated based on a detailed numerical code (NEWDRIFT+), presented next.

3.3.2 Wind resistance

The wind acts as an external force on the vessel, resulting in the wind resistance component. There are several calculation methods, which are based either on semi-empirical formulas, experimental results of models in tunnel tests or CFD results. This study considers only the longitudinal component of the wind resistance. Nevertheless, it is worth mentioning that the transverse component could impose a drift motion that can also increase the resistance or affect ship's manoeuvrability, which, however, are not considered in this study. Figure 4 describes the calculation process.

Figure 4. Wind resistance flowchart

3.3.3 Impact of sea currents

Sea currents play a major role in route selection by significantly affecting ship's speed and thus the fuel oil consumption. The effect of sea currents is considered by correcting the vessel's speed, while calm water resistance is also corrected by adjusting the speed over ground using the current speed.

3.3.4 Propeller and main engine modelling

A ship's total resistance $({{R_{\textrm{tot}}}} )$ can be derived by the several components according to the following:

(1)\begin{equation}{R_{\textrm{tot}}} = {R_{\textrm{calm}}} + {R_{\textrm{wave}}} + {R_{\textrm{wind}}}\end{equation}

where ${R_{\textrm{calm}}}$, ${R_{\textrm{wave}}}$ and ${R_{\textrm{wind}}}\;$ refer to calm water, added wave and wind components of the resistance, respectively. Having determined the ship's total resistance, the brake power can be calculated (e.g. Wang et al., 2020):

(2)\begin{equation}{P_\textrm{B}} = \frac{{{R_{\textrm{tot}}}{V_\textrm{S}}}}{{k{\eta _\textrm{s}}{\eta _0}{\eta _\textrm{H}}{\eta _\textrm{R}}}}\end{equation}

where ${R_{\textrm{tot}}}$ is the total resistance, ${V_S}$ is the ship's speed, k corresponds to the number of the propellers, ${\eta _s}$ corresponds to the shaft efficiency, ${\eta _0}$ refers to the propeller's open water efficiency which is equal to (k _TJ)/(k _Q2π), ${\eta _H}$ refers to the hull efficiency equals to (1–t)/(1–w), J is the advance coefficient, and t and w are the trust deduction and wake coefficient, respectively, while ${\eta _R}$ corresponds to the propeller's total rotational coefficient.

When a propeller's characteristics (diameter, pitch ratio, number of blades and the ratio of the expanded blade area) are known, and thrust and torque coefficient curves are provided, then the quantities ${k_\textrm{T}}$, k _Q, $J$, $n$ are possible to be determined from a propeller's open water diagram, which corresponds to calm water conditions. Firstly, the following quantity is calculated:

(3)\begin{equation}\frac{{{k_T}}}{{{J^2}}} = \frac{{T/({\rho {n^2}{D^4}} )}}{{{{({V/nD} )}^2}}} = \frac{T}{{\rho {V^2}{D^2}}} = CC\end{equation}

where $\rho$ is the water density, n is the propeller's revolutions, D is the propeller's diameter and T is the thrust. From the intersection of the curve ${k_\textrm{T}} = CC{J^2}$ with the curve (${k_\textrm{T}} - J$) of the open water diagram of the propeller, the values of $J$, ${k_\textrm{T}}$, k _Q and ${\eta _0}$ can be determined. Then,

(4)\begin{equation}n = \frac{{{\textrm{V}_{ad}}}}{{\textrm{JD}}}\end{equation}

where ${V_{\textrm{ad}}} = {V_s}({1 - w} )$ is the advance speed of the propeller when operating in the ship's wake.

The fuel oil consumption will be calculated based on the main engine's model by using the operating point (required brake power from Equation (2) and main's engine revolution from Equation (4) when no gear box is used) and the corresponding specific fuel oil consumption (SFOC). Using the manufacturer's manual,¹ the engine's SFOC map is created following the instructions presented in pages 2 ⋅ 09–2 ⋅ 11. Specifically, the SFOC for an arbitrary load within the loading diagram is needed. Knowing the SFOC value for the nominal MCR L1 rating, reduction rates are provided in the manual for the propeller curve and the constant speed curve for a range of loads. Then, for the operating points that lie between these curves, interpolation is carried out, while for the others, extrapolation. Having calculated SFOC values for a grid of (P,n) pairs, the SFOC contours shown in Figure 5 along with the loading diagram can be derived.

Figure 5. Main engine's SFOC map

In our case, when the point of operation exceeds the main engine's diagram,² the route is defined as unfeasible. The process is repeated and when, finally, the ship's speed is attainable, the SFOC value for the specific operating conditions is calculated by the engine's SFOC diagram. Note that when a shaft generator or a shaft power in system are used, respective corrections in the fuel oil consumption are needed. Furthermore, this analysis assumes no engine degradation, while propulsion characteristics are not changed under the effect of weather conditions (e.g. due to ship motions).

3.4 Speed and route management

As mentioned, wave-added resistance can lead to speed loss for the same engine load. Apart from the adverse weather conditions, the captain can also decide to reduce the speed to avoid unsafe effects, such as slamming or propeller immersion. On the other hand, the estimated time of arrival (ETA) is a major factor leading captains to speed up or slow down during the voyage to reach the destination port at the arranged time and avoid port traffic. Moreover, storm avoidance through speeding up or slowing down is also a well-known method for weather routing applications.

In the current framework, the speed management is carried out in the following ways:

• • The ship maintains the commanded speed for the entire route or for a route segment by examining whether it can be achieved under the encountered weather conditions.

• • For a specific ETA, the algorithm can optimise not only the selected route but also the speed between waypoints or any other selected points, aiming at avoiding storms or other dangerous phenomena and reaching the destination port on time.

• • The virtual arrival concept is also considered. Information about port delays can be sent on board during the voyage and an updated speed optimisation can be derived by taking advantage of the voluntary speed loss concept.

It is obvious from these scenarios that route and speed optimisation are interwoven issues that must be handled as one to succeed significant fuels savings.

4. Case study

This section presents a route optimisation study of a containership whose main characteristics are presented in Table 1. The selected route is from Portugal (Port of Lisbon) to Canada (Port of Halifax, Nova Scotia) and the following cases are examined:

• • Case A, where a constant speed of 16 kn is assumed, while no safety criteria are considered.

• • Case B, where again a constant speed of 16 kn is assumed, while safety criteria are specified by setting the maximum of the significant heave values, derived from each route-leg, equals to 3 ⋅ 5 m and the maximum significant roll is equal to 5 deg.

• • Case C, where speed is optimised within a range of 14 to 18 kn and the maximum duration of the voyage is set to 152 ⋅ 5 h.

For all the optimisation cases, the objective function is the minimisation of fuel oil consumption of the main engine. Next are presented the tools and methods used for the required calculations described in the previous section.

Table 1. Main characteristics of the examined containership

4.1 Weather data

All weather data are obtained from Copernicus Marine Environment Monitoring Service (CMEMS). Wave data sets (including significant wave height and its direction as well as the wave peak period) are available at 0 ⋅ 2° × 0 ⋅ 2° spatial resolution and for 3-hour temporal resolution. In addition, wind data sets (including true wind speed and its direction) are available at 0 ⋅ 25° × 0 ⋅ 25° spatial resolution and for 6-hour temporal resolution, while wave current data sets (including current speed and its direction) at 0 ⋅ 25° × 0 ⋅ 25° spatial resolution and 12-hour temporal resolution.

4.2 Calm water resistance

Calm water resistance was determined by CFD calculations using FreSCo+ (Hafermann, 2007). Note that the CFD calculations were significantly close to the ship's model tests, while being utilised to cover a speed range lower than the available one of the model tests (see also Zaraphonitis et al., 2021).

4.3 Wind resistance

The wind speed experienced by a moving vessel at sea is called apparent wind speed and differs from the true wind speed that the weather forecasts provide. So, it is necessary to calculate the apparent wind speed ${V_\textrm{a}}$ and direction a, where according to Figure 6, V _s is the ship's speed, V_t is the true wind speed, and $\varTheta$ is the angle formed between the two vectors mentioned, Y = 90–$\varTheta$, and:

(5)\begin{gather}{V_a} = \sqrt {{{({{V_t}\cos Y} )}^2} + {{({{V_t}\sin Y + {V_s}} )}^2}} \end{gather}

(6)\begin{gather}a = {\tan ^{ - 1}}\left( {\frac{{{V_t}\cos Y}}{{{V_t}\sin Y + {V_s}}}} \right)\end{gather}

Figure 6. Apparent wind calculation

Apparent wind speed and its direction are two major parameters to estimate wind resistance. In the presented cases, the wind resistance was determined by using Blendermann's method and coefficients (Blendermann, 1994), which uses a semi-empirical loading function based on wind tunnel tests. Wind resistance coefficients were also compared (Figure 7), with results obtained using Fujiwara's regression formula (ITTC, 2017) as well as the reference values from the sample containership appeared in Annex F.3 of ITTC (2017).

Figure 7. Comparison of wind resistance coefficients

4.4 Added wave resistance

To calculate the added wave resistance in an irregular sea, the spectral analysis technique (Lloyd, 1998) was employed as presented next, while the JONSWAP wave spectrum was used. Given that the significant wave height ${H_\textrm{s}}$ and the modal wave period ${T_0}$ describe the sea state, the JONSWAP spectrum is defined as (e.g. Lloyd, 1998):

(7)\begin{equation}{S_{J\zeta }}(\omega )= 0.658\; C\; {S_{B\zeta }}(\omega )\end{equation}

where $C = {\gamma ^J}$, $J = \textrm{exp}\left[ {\frac{{ - 1}}{{2{\sigma^2}}}{{\left( {\frac{{\omega {T_0}}}{{2\pi }} - 1} \right)}^2}} \right]$, while: $\sigma = 0.07\,\textrm{for}\,\omega < \frac{{2\pi }}{{{{\rm T}_0}}}$ and $\sigma = 0.09\,\textrm{for}\,\omega > \frac{{2\pi }}{{{{\rm T}_0}}}$, ${S_B}\zeta (\omega )= \frac{A}{{{\omega ^5}}}exp \left( {\frac{{ - B}}{{{\omega^4}}}} \right)$, while: $A = 487\frac{{H_s^2}}{{T_0^4}}$ and $B = \frac{{1949}}{{T_0^4}}$.

For the current application the value of the peak-enhancement parameter γ is equal to 3.3, which is the typical value when the JONSWAP spectrum is used, however it shall be noticed that for fully developed and open seas a more refined estimation of this parameter could be considered (e.g. Ochi, 1998).

Added wave resistance depends on the wave frequency, and thus for a vessel advancing through sea waves, the encountered wave frequency is needed:

(8)\begin{equation}{\omega _e} = \omega - \frac{{{\omega ^2}{V_s}cos\mu }}{g}\end{equation}

where $\omega$ is the wave frequency, $\mu$ is the angle defined by the ship's heading and wave direction (which is defined as the relative wave heading), and $\textrm{g}$ = 9 ⋅ 81 m/sec².

A frequency interval $\delta \omega$ around the frequency $\omega$ of the wave energy spectrum is transformed as:

(9)\begin{equation}\delta {\omega _e} = \left( {1 - \frac{{2\omega {V_s}}}{g}cos\mu } \right)\delta \omega \; \; \end{equation}

The area below the spectrum within the interval $\mathrm{\delta \omega }$ is relevant to the energy carried by that frequency range. This energy is not changed by the spectrum transformation and, as a result, the area bounded by the frequency range $\delta \omega$ in the wave energy spectrum will be equal to the area bounded by the frequency range $\delta {\omega _\textrm{e}}$ in the encounter wave spectrum. If $\delta \omega$ and $\delta {\omega _\textrm{e}}$ are infinitesimal, the encounter wave spectrum ordinates are:

(10)\begin{equation}{S_\zeta }({{\omega_e}} )= {S_\zeta }(\omega )\frac{{d\omega }}{{d{\omega _e}}} = {S_\zeta }(\omega )\frac{g}{{g - 2\omega {V_s}\cos \mu }}\end{equation}

The total areas under the two spectra are equal as the total wave energy as well as the significant wave height are not affected by this conversion. Using spectral analysis, the amplitude of a sine wave representing a small range of frequencies $\delta {\omega _\textrm{e}}$ is

(11)\begin{equation}{\zeta _0} = \sqrt {2{S_\zeta }({{\omega_e}} )\delta {\omega _e}} \end{equation}

and the single sine wave results in an added resistance equal to

(12)\begin{equation}2{C_{aw}}{S_\zeta }({{\omega_e}} )\delta {\omega _e} = 2{C_{aw}}{S_\zeta }(\omega )\mathrm{\delta \omega }\end{equation}

where ${C_{\textrm{aw}}} = {R_{\textrm{aw}}}/{\mathrm{\zeta }_0}^2$ is the added resistance response function. Thus, when considering all the wave components, the total added resistance is (Lloyd, 1998)

(13)\begin{equation}{\bar{R}_{aw}} = 2\int_0^\infty {{C_{aw}}{S_\zeta }({{\omega_{}}} )d\omega } \end{equation}

The added resistance component that corresponds to Equation (12) is pre-calculated by NewDrift+ software, for a range of different wave frequencies, headings and speeds. NEWDRIFT+ is NTUA's in-house software for seakeeping analysis based on a 3D linear panel method and can calculate all ship's responses and added resistance, assuming that the ship is under the effect of regular waves (Liu et al., 2011).

For each one of the discrete frequencies of the generated wave spectrum, NewDrift+ data are interpolated first regarding the wave frequency, then the heading and finally regarding the ship's speed to derive the added resistance (Zaraphonitis et al., 2021) Finally, this methodology is used to calculate the total added resistance in a sea state (Figure 8).

Figure 8. Added wave resistance

4.5 Seakeeping data

The results presented herein are based on the use of the NEWDRIFT+ software for the evaluation of the ship's responses in waves. To represent the wetted surface, the hull is discretised by using trilateral or quadrilateral panels, as shown in Figure 9.

Figure 9. Containership's panelisation

Following Lloyd (1998) for the irregular sea state, similarly to the methodology for added wave resistance, the motion energy spectra must be calculated to estimate statistical quantities of the ship motion responses. The resulting motion amplitudes and phases depend on the vessel's speed, the waves’ heading and the encounter frequency. Within linear seakeeping theory and considering the ship as a rigid body, ship motion amplitudes are proportional to the wave amplitude and can be expressed in the form of nondimensional transfer functions. Transfer functions are defined as ship motion amplitudes divided by the wave amplitude (ζ₀) for linear motions and by the wave slope amplitude (k ζ₀) for angular motions, where k is the wave number.

Motion energy spectra are obtained by multiplying the encounter wave spectrum with the squared motion transfer functions of the corresponding encounter frequency. The motion energy spectrum ordinates for each encounter frequency ${S_\mathrm{\zeta }}({{\omega_\textrm{e}}} )$ is calculated by

(14)\begin{equation}{S_{{x_i}}}({{\omega_e}} )= {S_\zeta }({{\omega_e}} ){\left( {\frac{{{x_{i0}}}}{{{\zeta_0}}}} \right)^2}\end{equation}

where i = 1,2,3 correspond to surge, sway and heave motion, respectively. For roll, pitch and yaw motions (i = 4,5,6), the wave slope spectrum must be calculated for the nth sine wave component:

(15)\begin{equation}{S_a}({{\omega_n}} )= \frac{{\omega _n^4}}{{{g^2}}}{S_\zeta }({{\omega_n}} )\end{equation}

Then, the motion energy for roll, pitch and yaw is

(16)\begin{equation}{S_{{x_i}}}({{\omega_e}} )= {S_a}({{\omega_e}} ){\left( {\frac{{{x_{i0}}}}{{k{\zeta_0}}}} \right)^2}\end{equation}

The respective variance and route mean square (rms) values for each motion are calculated according to

(17)\begin{gather}{m_0} = \int_0^\infty {{S_{{x_i}}}({{\omega_e}} )d{\omega _e}} \end{gather}

(18)\begin{gather}{\sigma _0} = \sqrt {{m_0}} \end{gather}

The significant value of the motion energy spectrum is equal to four times the rms motion, since the wave spectra as well as the motion energy spectra are considered narrow-banded spectra. For some critical locations onboard, such as the bridge and the bow, the information concerning the motion displacement, velocity and acceleration could be crucial. The motions on any point on the ship can be calculated based on the motions of the centre of gravity as calculated previously. The angular motions (roll, pitch and yaw) are constant on every point on the vessel, whereas local linear motions (surge, sway and heave) depend on the location on the ship. For example, Figure 10 shows the lateral and vertical accelerations for two critical locations onboard for several wave periods and for a specific value of the significant wave height and relative wave heading $\mu$.

Figure 10. Lateral and vertical accelerations at bridge and at the bow for ship's speed 16 kn, significant wave height 4 m and relative wave heading equal to μ = 150 deg. (Head seas correspond to 180 deg.)

4.6 Optimisation cases

The optimisation is performed through the optimisation toolbox available in MATLAB using a genetic algorithm (see for example Matlab, 2022). Specifically, 100 generations with a population of 100 are chosen for the optimisation processes presented next. The search space is predefined to save computational time, but the optimisation procedure is not based on a fixed-grid technique, as happens in other optimisation approaches (e.g. path finding algorithms, conventional dynamic programming, etc.). Optimisation is also performed for optimal heading (Case A) and both for optimal heading and for nonconstant main engine's and propeller's revolutions (Case C), to approach the route planning problem under realistic conditions. Moreover, the speed optimisation in Case C is combined with a supposed required ETA.

In Case A, the speed is set constant to 16 kn and no safety criteria are considered. Figure 11(a) shows the genetic algorithm evolution during 100 generation of 100 population each.³. It is obvious from the figure that the algorithm converges to the optimal result during almost the 10th generation. On the other hand, in Case B, the speed is again kept constant to 16 kn, but the next safety criteria are added:

• • The threshold value for the significant heave is set to 3 ⋅ 5 m.

• • The respective value for the significant roll angle is set to 5°.

These values were chosen for the demonstration of the optimisation study, however operability criteria based on the ability for carrying out several activities onboard either for the crew or passengers can be implemented (e.g. Perera et al., 2012).

Figure 11. Algorithm's evolution for all cases, where the red circle indicates the optimum ones

The evolution of the algorithm is shown in Figure 11(b). Note also that the number of population and the number of generations are the same as in Case A. As expected, the number of unfeasible routes is higher, and the algorithm converges to the optimal result later in comparison with Case A. Finally, in Case C five new parameters that represent the speed for each route's leg are added to the optimisation process. A maximum travel time of 152 ⋅ 5 h is also added as a constraint. The evolution of the algorithm is shown in Figure 11(c).

A common approach is to compare the results of the optimal routes with the respective ones of the orthodromic and the loxodromic paths, while considering for both of them a constant speed of 16 kn. In the following figures (Figures 12–16) characteristics of these routes are shown in comparison with the optimal ones from each case. The points listed in x-axes represent the intersection points with the main grid, as described in Section 3.2.1.

Figure 12. Significant wave height encountered for the orthodrome, the loxodrome, the optimal path with no safety criteria (Case A) and the optimal path with safety criteria (Case B)

Figure 13. Relative wave heading for orthodrome, loxodrome and optimal routes (Case A & Case B)

Figure 14. Calm water resistance, accounting the effect of currents, for orthodrome, loxodrome and optimal routes (Case A & Case B)

Figure 15. Added resistance due to wind for orthodrome, loxodrome and optimal routes (Case A & Case B)

Figure 16. Wave added resistance for orthodrome, loxodrome and optimal routes (Case A & Case B)

As derived by the optimisation procedure, to save fuel, it is preferable to select a route that either passes from regions with low values of the significant wave height (Figure 12) or from regions where the ship encounters following seas (Figure 13). The optimal route in Case B in order to satisfy the set constraints completely avoids head seas, but the avoidance of the bad weather leads to an extended journey with higher fuel oil consumption. Similarly, Figure 14 presents the calm water resistance for the same four paths. Again, the optimal route for Case B follows the favourable path concerning the encountered sea currents, whereas the orthodrome and the loxodrome encounters unfavourable sea currents that increase the calm water resistance and, consequently, the total resistance. Figures 15 and 16 present the added wind and wave resistance, respectively. Obviously, even though the contribution of the wind resistance component to the total resistance is much lower compared with those of calm water and wave, the variations derived by the selections of the optimal route when considering safety criteria reveal that it is a significant element to be considered. Another issue to be discussed is that although the sailing sea state is characterised by heavy waves, the wind resistance is low, since the analysis considers the combined effect of wind and swell waves. As expected, both the optimal routes avoid the passage from regions where the wave-added resistance is high to reduce the required fuel. In addition, Figure 17 shows the resultant optimised speeds at each leg for the optimal path found by the genetic algorithm for Case C. The algorithm suggests reducing the speed in regions where the significant wave height is high or when encountering head seas. Finally, Figure 18 presents the weather evolution each day of travel.

Figure 17. Optimized speed at each leg

Figure 18. Significant wave height evolution within the area of the optimal route (Case A) for each day of the voyage

4.7 Summary of optimisation results

Table 2 summarises the results for all the cases examined, where the maximum of the significant vertical acceleration values, as derived considering each route leg, at the bridge and at the bow have also been included. The optimal route with no safety criteria leads to 4 ⋅ 9% fuel savings compared with the orthodromic path, which is the shortest path and reaches the port of destination almost 2 h earlier than the optimal one. Moreover, when the safety constraints are considered, the optimal path found by the genetic algorithm requires 10% more fuel and 22 h more time than the orthodrome at the same constant speed (16 kn). When speed is optimised (Case C), while the total time duration is set to 152 ⋅ 5 h as a constraint, the optimal path ends in a slight decrease of fuel oil consumption compared with the optimal with constant speed and no safety criteria. All routes are shown in Figure 19.

Figure 19. The derived optimized routes, where the red line = loxodrome, red curve = orthodrome, green polyline = optimal with safety criteria, purple polyline = optimal with optimized speed, black polyline = optimal with no safety criteria

Table 2. Results from all optimisation cases

5. Concluding remarks

The paper presents a weather routing tool based on detailed ship's characteristics and forecasted weather data along the route by providing examples on the achieved fuel savings. The application presented consists of several tools connected to the main body of the application in a modular way. This is considered an advantage of the tool, as it provides the user the flexibility to utilise any data and tools available. Furthermore, the fact that optimisation is based on an evolutionary approach provides the benefit of exploring a variety of waypoint combinations lying in a predefined search space. The search space is fixed in terms of boundaries to save time during the optimisation process. Thus, the algorithm can test any possible position (longitude and latitude), compared with fixed-grid techniques where not only the search boundaries are predefined but also all the candidate grid points (waypoints). Furthermore, by looking the ship's route planning in a more realistic point of view, the presented tool can optimise the heading of the ship from one port to the destination, but also optimise simultaneously the speed between every two waypoints. Route and speed optimisation must be considered as parts of the same problem in order to achieve a significant amount of fuel savings. ETAs calculations are also performed by the tool in conjunction with speed management, a fact that leads not only to savings regarding fuel oil consumption but also by avoiding delays at ports.

Moreover, either historical or forecasted weather data can be used, depending on the needs of the examined study. The optimisation can be performed either at the beginning of the voyage or it can be repeated at any time during the sailing (for example, optimisation can be repeated when new forecasted data is available, which depends on the temporal resolution of the data). In addition, the tool entails safety aspects in terms of seakeeping criteria, by considering detailed ship characteristics and utilising an advanced numerical code. This addition results in treating the problem of weather routing not only from the energy efficiency viewpoint but also by imposing constraints due to safety, and this can be also considered a key asset of the developed tool.

Most of the time, the shortest path is neither the safest nor the most fuel efficient. A compromise must be achieved, and that means that the route planning is a crucial procedure that must take into consideration accurate weather forecasted data and the specific characteristics of each vessel, since every vessel is different and performs different even when talking about sister ships (e.g. they most likely would have different levels of fouling, different machinery set up, etc.). Seakeeping and safety criteria are bound to this procedure and when considered can result to significant changes in the selection of the optimised route.

The modular character of the presented tool allows for examining in the future several effects, by appropriately updating the models presented. Firstly, the current tool is based on a physics-based model for the calculation of the main's engine fuel oil consumption where, for the calculation of the resistance, tools of different level of detail can be employed either based on semi-empirical formulas or on detailed numerical codes. An alternative approach could be the utilisation of a data-driven model for the calculation of the fuel consumption based on a set of key features that would allow consideration of the current hull and propeller performance status. On the other hand, the impact of the hull and propeller fouling could be also implemented in a physics-based approach by modelling the impact of the fouling on ship's resistance and propeller's characteristics (see, for example, Liu et al., 2021). Moreover, the effect of ship motions on the propeller performance is an aspect worth investigating, as it could result in deviating from the uniform flow assumed when using the propeller's open-water characteristics, while the separation of wind and swell effect could also be explored. Sharp turns from one waypoint to the next can also be avoided, ensuring manoeuvrability even for large vessels by adding a turning penalty, which would also contribute to a smoother optimal path. In addition, the degradation of machinery performance can be also considered through updating the SFOC values. Finally, by employing different type of wave spectra, the encountered wave conditions could be appropriately captured.