*Article* **Development of a New Ship Adaptive Weather Routing Model Based on Seakeeping Analysis and Optimization**

#### **Silvia Pennino \*, Salvatore Gaglione, Anna Innac, Vincenzo Piscopo and Antonio Scamardella**

Department of Science and Technology, University of Naples "Parthenope", Centro Direzionale Isola C4, 80143 Naples, Italy; salvatore.gaglione@uniparthenope.it (S.G.); anna.innac@uniparthenope.it (A.I.); vincenzo.piscopo@uniparthenope.it (V.P.); antonio.scamardella@uniparthenope.it (A.S.) **\*** Correspondence: silvia.pennino@uniparthenope.it; Tel.: +39-081-547-6686

Received: 17 March 2020; Accepted: 8 April 2020; Published: 10 April 2020

**Abstract:** This paper provides a new adaptive weather routing model, based on the Dijkstra shortest path algorithm, aiming to select the optimal route that maximizes the ship performances in a seaway. The model is based on a set of ship motion-limiting criteria and on the weather forecast maps, providing the sea state conditions the ship is expected to encounter along the scheduled route. The new adaptive weather routing model is applied to optimize the scheduled route in the Northern Atlantic Ocean of the S175 containership, assumed as a reference vessel, based on the weather forecast data provided by the Global WAve Model (GWAM). In the analysis, both wave and combined wind/swell wave conditions are embodied to investigate the incidence on the optimum route assessment. Furthermore, the effect of the vessel speed on the optimum route detection is also investigated. Current results clearly show that it is possible to achieve appreciable improvements, up to 50% of the ship seakeeping performances, without excessively increasing the route length and the voyage duration.

**Keywords:** adaptive weather routing; Seakeeping Performance Index; route optimization; Dijkstra algorithm

#### **1. Introduction**

Maritime trades are strictly dependent on the weather conditions the ship is expected to encounter along her route, provided that excessive ship motions may cause damage or cargo loss, as well as the decrease of the onboard comfort level, with a negative impact on the safety of navigation for both cargo and passenger ships. In this respect, even if the route selection is entrusted to the ship master, the adaptive weather routing algorithms can provide a significant support to the decision-making process, to ensure the proper balancing between the safety of navigation and the related economic impact, in terms of voyage duration and fuel savings.

Nevertheless, as further discussed in Section 2, the impact of weather routing on fuel savings is expected to be marginal, provided that obtained results are not particularly encouraging, as no substantial fuel savings can be gathered if the ship route is varied, as regards the minimum distance route, without reducing the vessel speed. Hence, in the following research, a new adaptive weather routing model is developed to select the best route that maximizes the ship performances in a seaway, based on a set of proper seakeeping criteria. In fact, even if the route selection has a low impact in terms of fuel savings, an appreciable increase of the ship seakeeping performances can be gathered without excessively changing the route, as regards the minimum distance route. The new route optimization procedure, based on the shortest path algorithm, is developed in this paper to detect the best route, based on the weather forecasts the ship is expected to encounter up to the arrival in port. The algorithm

is applied to the S715 containership, assumed as a reference vessel, to investigate the effect of the selected seakeeping parameters and vessel speed on the best route selection, as well as on the expected voyage duration and path length. Hence, the main aim of the current research is to develop a new adaptive weather routing model aiming to: (i) maximize the seakeeping performances in a seaway, (ii) investigate the incidence of unimodal/bimodal spectra and vessel speed on the improvement of the seakeeping performances, and (iii) carry out a sensitivity analysis to investigate the effect of each seakeeping parameter on the optimum route selection.

#### **2. Literature Review**

The weather routing problem was addressed by many authors and different approaches were embodied in the past. Following the pioneering works by James [1], Zoppoli [2], and Papadakis and Perakis [3], the first generation of weather routing criteria aimed to minimize the voyage duration, and consequently, the fuel consumption, mainly neglecting the impact of the best route selection on the ship performances in a seaway. Nevertheless, in the last years, the voyage optimization problem was approached in a wide sense, taking into account the ship seakeeping performances, apart from the fuel consumption. In this respect, the second-generation algorithms, developed as general optimization tools and applied to the ship voyage optimization problem, are generally categorized into: dynamic programming, genetic algorithms, and pathfinding methods.

The most used techniques include multi-objective genetic algorithms [4–7] that stochastically solve a discretized nonlinear optimization problem [8], based on the ship course and her velocity profile that, in turn, is represented by a set of parametric curves. In this respect, Szłapczynska and Smierzchalski [9] proposed a multi-criteria evolutionary weather-routing algorithm, based on an iterative process of population development, resulting in a Pareto-optimal set of solutions. Shao and Zhou [10] presented a three-dimensional (3D) dynamic-programming based on the original work by De Wit [11], devoted to optimizing the ship speed and the heading angle in the route planning. Zaccone et al. [12] developed a 3D dynamic-programming-based ship voyage optimization method, that attempts to select the optimal path and speed profile, based on the expected weather forecasting along the vessel route.

Another technique, devoted to managing the weather routing optimization procedure, is the path-finding method. In this respect, Padhy et al. [13] embodied the Dijkstra [14] algorithm to solve the ship routing problem, while Veneti [15] provided an improved solution based on the exact time-dependent bi-objective shortest path algorithm, to optimize two different conflicting objectives, namely, the fuel consumption and the expected risk along the route. In Reference [16], an overview of weather routing and safe ship handling approaches is presented, addressing the importance of the complementarity between both approaches in such situations and their extensive implementation to achieve optimal and safe navigation conditions in shipping. In Reference [17], some modifications are proposed to find the optimal path in ice-covered waters in order to take into account the possible involvement of an icebreaker. In this view, the icebreaker assistance is considered as an integral part of the overall formulation of the route optimization problem in the frame of graph-based and wave-based numerical routing algorithms.

A solution method for a real-life problem of routing and scheduling ships, motivated by the production and logistic activities of an oil company operating in Brazil, is described in Reference [18], where the authors proposed a multi-start heuristic that makes use of intensification and diversification strategies to minimize both transportation and docking costs, as well as avoid consecutive dockings of the same ship in different platforms. In Reference [19], the authors presented a study about a data-driven framework for decision support in the selection of optimal ship routes based on available weather forecasts and the calculation of fuel oil consumption over alternative routes. In detail, a modified version of the Dijkstra algorithm was recursively applied until the optimal route was detected.

#### **3. Adaptive Weather Routing Model**

#### *3.1. Shortest Path Detection*

The adaptive weather routing model, developed in Section 3, is based on the Dijkstra algorithm, that allows for the detection of the shortest path between two points of known coordinates, by maximizing a reference non-negative cost-function which is, in the current analysis, the Seakeeping Performance Index (*SPI*), provided by Equation (1). The Dijkstra algorithm is based on two problems, namely, the construction of the tree of minimum total length between a set of nodes, and the research of the best path between the start and the endpoint by minimizing an assigned cost function. The points, connected by a set of branches, form the required tree, from which a uniquely defined path between two nodes is obtained. In the current analysis, the cost function is the complement to 1 of the *SPI* provided by Equation (1). Besides, the algorithm is not applied statically but it is iteratively updated by changing the starting point, based on the elapsed time and the path already travelled by the ship. Before assessing the ship performances in a seaway, a network of nodes, between the start and the endpoints of the scheduled route, needs to be preliminarily provided to define a region that constrains the ship route. The network, generated as detailed by Lee et al. [20], is a special form of a Dijkstra graph consisting of nodes and arcs, associated to a set of parameters such as the path length, the cost, or the transit time. The network is generated according to the following steps: (i) the center point C of the great circle, black line in Figure 1a, between S and E, is determined, (ii) around C, a set of center points (red dots), with a fixed step distance, are considered, belonging to the rhumb-line perpendicular to the great circle, and (iii) a network of nodes, blue dots in Figure 1a, are generated, considering, in the first half of the great circle, a set of perpendicular rhumb-lines, as well as in the second half part. On each rhumb-line, a proper step distance between nodes is adopted. Each node is characterized by static information, such as the geographic coordinates, a dynamic grid with the weather data, namely, the significant wave height, the mean wave period, and the prevailing wave direction, which are systematically updated once a new dataset is available. After providing the network of nodes, a set of arcs, connecting two adjacent grid points, is created, and a cost-function, depending on the *SPI* values at the arc extremities, is assigned to detect the best route, as schematically shown in the flow chart reported in Figure 1b. Each segment of the optimal route corresponds to the maximum of the Seakeeping Performance Index among the set of possible values referenced to the calculation point.

**Figure 1.** Adaptive weather routing algorithm: (**a**) Grid point network, (**b**) flow chart. *SPI*: Seakeeping Performance Index.

#### *3.2. Assessment of Ship Seakeeping Performances*

In the current analysis, the ship seakeeping performances are determined on the basis of five reference criteria, namely: (i) the amplitude of the pitch motion, (ii) the relative vertical acceleration at the ship forward perpendicular, (iii) the probability of slamming occurrence, (iv) the probability of green water on deck, and (v) the Motion Sickness Incidence (*MSI*), according to the following equation:

$$SPI = \max\left\{0; \left(1 - \frac{rms\_p}{rms\_{p,l}}\right) \cdot \left(1 - \frac{rms\_a}{rms\_{a,l}}\right) \cdot \left(1 - \frac{p\_{sl}}{p\_{sl,l}}\right) \cdot \left(1 - \frac{p\_{wd}}{p\_{wd,l}}\right) \cdot \left(1 - \frac{MSI}{MSI\_l}\right)\right\} \tag{1}$$

where *rmsp* (*rmsa*) is the RMS, root mean square, of the pitch motion amplitude (relative vertical acceleration), *psl* (*pwd*) is the probability of occurrence of slamming (water on deck), *MSI* is the Motion Sickness Incidence, while the denominators represent the relevant limit values. It is noticed that when any seakeeping index is greater than the limit value, the *SPI* is null, so satisfying the non-negative condition is required to apply the Dijkstra method. The limit values of the pitch amplitude and the *MSI* are based on the NATO STANAG, Standardization Agreement, 4154 criteria, while the remaining ones comply with the NORDFORSK 1987 criteria [21,22], as detailed in Table 1.

**Table 1.** General operability limiting criteria for ships.


\* at intermediate values, linear interpolation is applied.

The RMS of the pitch motion amplitude is determined according to the following equation, based on the ship Response Amplitude Operator [23]:

$$rms\_p = \sqrt{\int\_0^\infty \left| H\_5(\omega\_\varepsilon) \right|^2 S\_\varepsilon(\omega\_\varepsilon)}\tag{2}$$

where *H*<sup>5</sup> is the speed-dependent pitch motion transfer function to be determined as detailed in Appendix A, *<sup>S</sup>*<sup>ς</sup> is the wave spectrum, and <sup>ω</sup>*<sup>e</sup>* = <sup>ω</sup> <sup>−</sup> <sup>ω</sup>2<sup>ψ</sup> is the encounter wave frequency that satisfies the Doppler shift equation, depending on the absolute wave frequency ω and the factor ψ = *Ucos*μ/*g*, where the vessel speed is denoted by *U* and μ denotes the heading angle. The RMS of the relative vertical acceleration is obtained by combining the heave and pitch motions at the ship forward perpendicular [23]:

$$
\sigma m \mathbf{s}\_{\mathbf{z}} = \sqrt{\int\_{0}^{\infty} \left| H\_{3}(\omega\_{\varepsilon}) - \overline{\mathbf{x}} H\_{5}(\omega\_{\varepsilon}) - e^{-i\overline{\mathbf{x}}\overline{\mathbf{x}}\cos\mu} \right|^{2} \omega\_{\varepsilon}^{4} S\_{\zeta}(\omega\_{\varepsilon})} \tag{3}
$$

where *H*<sup>3</sup> is the heave motion transfer function, *x* is the longitudinal distance (fwd+) of the ship forward perpendicular from its center of mass, and *k* = ω2/*g* is the wave number satisfying the deep-water condition. As concerns the slamming occurrence, it is determined according to the formula proposed by Faltinsen [24]:

$$p\_{sl} = \varepsilon^{-(\frac{v\_{cr}^2}{2\mathcal{C}\_s^2 m\_{2J}} + \frac{d^2}{2\mathcal{C}\_s^2 m\_{0,r}})}\tag{4}$$

where *vcr* = 0.093 % *gL* is the threshold velocity, *Cs* is the swell up coefficient equal to 1 up to *FN* = 0.30 [25], and *d* is the ship draught at the forward perpendicular. In Equation (4), *m*0,*<sup>r</sup>* and *m*2,*<sup>r</sup>* denote the 0 and second-order spectral moments of the ship relative motion as regards the sea surface. As concerns the probability of green water on forward deck structures, it is determined in a quite similar manner [24]:

$$p\_{rad} = e^{-\frac{f\_b^2}{2C\_s^2 m\_{0\varphi}}}\tag{5}$$

where the freeboard at the ship forward perpendicular is denoted by *fb*. Finally, the Motion Sickness Index is determined according to the formulation developed by O'Hanlon and McCauley [26] and modified by Colwell [27]:

$$MSI = 100\Phi(z\_d)\Phi(z'\_T) \tag{6}$$

where Φ is the standard normal cumulative distribution, while *za* and *z <sup>T</sup>* are determined by the following equations:

$$z\_a = 2.128 \log\_{10} \left( \frac{a\_{MSI}}{g} \right) - 9.277 \log\_{10} (f\_m) - 5.809 [\log\_{10} (f\_m)]^2 - 1.851 \tag{7}$$

$$z\_T' = 1.134z\_d + 1.989 \log\_{10} \left(\frac{T}{60}\right) - 2.904\tag{8}$$

In Equations (7) and (8), *T* is the exposure time in s, while *aMSI* and *fm* are the RMS and the mean frequency in Hz of the ship vertical acceleration [27].

#### *3.3. Assessment of Weather Forecasting Data*

The availability of daily operational weather forecast data is a key point to assess the reference scenario required in the adaptive weather routing algorithm. In this respect, several software packages, with different resolutions, domains (from regional to global), and quality [28] are available, as proven by the agreements that shipping companies stipulate with one or more meteorological institutes, to obtain updated and reliable weather forecast data. Anyway, nowadays, the community generally follows the standardization established by the World Meteorological Organization (WMO), delivering all information in a self-descripting GRIB (GRIdded Binary) format [29], making the reading of the input data very easy. Among the variety of weather forecast codes, the third-generation Global WAve Model (GWAM), initially developed in the mid-80s by an international group of wave modelers [30], is probably the most embodied one by all research institutions around the world. It explicitly solves the wave transport equation, without any assumptions about the shape of the wave spectrum, therefore properly representing the physics of the wave evolution in compliance with the knowledge about the full set of freedom degrees of two-dimensional wave spectra. The model can be applied to any given regional or global grid, based on a certain topographic dataset. Besides, the grid resolution can be arbitrarily set in both space and time, while the wave propagation can be performed on latitudinal-longitudinal or Cartesian grids. The model outputs, embodied in the route optimization procedure, are the significant wave height and the mean wave period and direction of both wind wave and swell components.

#### **4. Input Data**

#### *4.1. The S175 Containership*

The S175 containership [31–33] is assumed as a reference vessel in the performed case study. The ship's main dimensions are listed in Table 2, while the zero-speed added mass and radiation damping are determined by the open source code Nemoh [34] and corrected to account for the forward speed effect, as detailed in Appendix A [35].


**Table 2.** Main particulars of the S175 containership.

#### *4.2. Route Selection*

The reference route falls in the North Atlantic Sea, as this trade lane is one of the most important routes in the global shipping industry and most of the key port infrastructures in Europe are continually upgraded to reflect the demands of shipping across this corridor. Besides, the weather conditions in this area can be rough and vary with seasons, influencing the obtained results. The selected voyage has (34◦ N, 60◦ W) departure and (32◦ N, 20◦ W) arrival coordinates, as depicted in Figure 2, that also provides the great circle or orthodrome (black) and the rhumb-line or loxodrome (red) routes, equal to 2005 and 2018 nm respectively, provided that they are generally combined by the master during the navigation [36,37]. In fact, the great circle or orthodrome is the shortest path route, but it requires to constantly change the vessel heading. On the contrary, the rhumb-line or loxodrome is slightly longer, but it allows the master to keep a constant heading in the absence of wind and sea current, as they are generally combined for practical purposes. The voyage duration is equal to about 7 days, considering a reference vessel speed of 12 knots, corresponding to Froude number *FN* = 0.15. From Figure 2, it is gathered that the start and the endpoint are slightly far from the departure and arrival ports. This choice is mainly due to the marine traffic congestion the ship experiences when it approaches the coastline, that makes the application of adaptive routing models challenging, as additional restraints, mainly related to neighboring vessels, arise.

**Figure 2.** Great circle (black) and rhumb-line (red) between routes' departure and arrival points.

#### *4.3. Weather Forecasting Data*

In the current analysis, the GRIB files were downloaded for the period from 3 up to 10 January 2020, with 0.25◦ × 0.25◦ grid spacing, 3-hour forecast interval, and two observation periods equal to 1 and 7 days, respectively. The free software XyGrib, able to download and show the meteorological data, was embodied to obtain the GRIB files. An example of the images providing the changes of the environmental conditions in the selected area is reported in Figure 3, that provides a series of GRIB files showing the time trend of the significant wave height, *Hs* = ) *Hs*,*wind*<sup>2</sup> + *Hs*,*swell*<sup>2</sup> from 06:00 UTC (Coordinated Universal Time) on 3 January 2020 up to 03:00 UTC on 4 January 2020, where the wind wave (swell) component is denoted by *Hs*,*wind* (*Hs*,*swell*).

**Figure 3.** Time trend of the significant wave height from 06:00 UTC on 3 January 2020 up to 03:00 UTC on 4 January 2020.

#### **5. Case Study**

#### *5.1. Optimum Route Detection*

After defining the input data, two GRIB files are embodied in the route optimization procedure. The former, referenced as "Case 1", is a 7-day forecast GRIB file, ranging from 3 up to 10 January 2020. The latter, referenced as "Case 2", is a 1-day GRIB file, updated daily in the same reference period, according to the effective position of the ship during the voyage. The same GRIB files, referenced as "Case 3" and "Case 4", are re-analyzed, including the swell data, apart from the wind-generated waves. In the first two cases, the JONSWAP [38] spectrum is applied in the weather routing optimization procedure, while in the remaining ones, the two-peak Torsethaugen [39] spectrum is embodied, to model weather conditions obtained by combined wind sea and swell components, coming from different directions.

In this respect, Figure 4a–d provide the optimal routes with reference to the four selected conditions. In all cases, the great circle route is highlighted in red, the optimum route is depicted in black, while the grey circles represent the starting and the ending points. Besides, Table 3 provides the difference, Δ*c*, between the great circle and optimum route, that maximizes the ship performances, together with the percentage variation, Δ*SPI*, of the Seakeeping Performance Index (*SPI*), as regards the values corresponding to the great circle route. Based on the current results, the percentage increase of the route length ranges from 1.1% up to 2.7%, while the *SPI* increase is much higher and ranges from 10% up to 40%. These results clearly show that the ship seakeeping performances can be highly improved, without significantly increasing the voyage length and, consequently, the fuel consumptions.

After carrying out this preliminary analysis, Figure 5a–d and Table 4 provide the same calculations obtained by reversing the travel direction, therefore inverting the departure and arrival coordinates. Also, in this case, the percentage increase of the route length is low and ranges from 3.1% up to 4.9%, while the *SPI* increase is much higher and ranges from about 10% up to 68%. All calculations were performed based on a reference vessel speed equal to 12 kn, corresponding to *FN* = 0.15.

**Figure 4.** Minimum distance (great circle) and optimal route detection, *FN* = 0.15: (**a**) Case 1; (**b**) Case 2; (**c**) Case 3; (**d**) Case 4.


**Table 3.** Results of optimal route detection, *FN* = 0.15.

**Figure 5.** Minimum distance (great circle) and optimal route detection, *FN* = 0.15 (reversed travel direction): (**a**) Case 1; (**b**) Case 2; (**c**) Case 3; (**d**) Case 4.



#### *5.2. E*ff*ect of Vessel Speed*

The incidence of the vessel speed on the optimum route detection is investigated, increasing the Froude number from 0.18 up to 0.24, with 0.03 step. Obviously, the vessel speed increase leads to a corresponding decrease of the travelling time that diminishes up to 138, 119, and 104 h, as well as to a slight variation of the weather conditions that the ship is expected to encounter during the voyage. In this respect, Figure 6a–d, Figures 7a–d and 8a–d provide the optimal routes, while Table 5 summarizes the obtained results. In the first case, corresponding to *FN* = 0.18, the route length increases from 0.6% up to 2.0%, if compared with the great circle one, while the *SPI* increases from 3.5% up to 31.5%. In the second condition, corresponding to *FN* = 0.21, the route length rises up from 0.6% to 6.1%, while the *SPI* increase ranges from 3.4% up to 39.8%. In the last condition, instead, the route length increases from 0.3% up to 1.6% and the *SPI* rises from 6.7% up to 32.1%. Based on the current results, it is confirmed that it is possible to achieve considerable improvements of the ship seakeeping performances, without significantly affecting the voyage length. Furthermore, it is verified that the vessel speed plays a fundamental role in the assessment of the optimal route, as well as on the increase of the ship seakeeping performances.

**Figure 6.** Minimum distance (great circle) and optimal route detection, *FN* = 0.18: (**a**) Case 1; (**b**) Case 2; (**c**) Case 3; (**d**) Case 4.

**Figure 7.** Minimum distance (great circle) and optimal route detection, *FN* = 0.21: (**a**) Case 1; (**b**) Case 2; (**c**) Case 3; (**d**) Case 4.

**Figure 8.** Minimum distance (great circle) and optimal route detection, *FN* = 0.24: (**a**) Case 1; (**b**) Case 2; (**c**) Case 3; (**d**) Case 4.


**Table 5.** Incidence of vessel speed on optimal route detection.

#### *5.3. Sensitivity Analysis*

The incidence on the optimum route detection of each seakeeping index provided by Equation (1) is investigated. The optimum routes, corresponding to *FN* = 0.15, are reported in Figure 9a–j with reference to both JONSWAP (left) and Torsethaugen (right) wave spectra, separately considering the following criteria, namely: (i) the RMS of pitch amplitude (a,b), (ii) the RMS of vertical acceleration (c,d), (iii) the slamming probability (e,f), (iv) the water on deck probability (g,h), and (v) the *MSI* (i,j).

**Figure 9.** Sensitivity analysis at *FN* = 0.15 based on JONSWAP (left) and Torsethaugen (right) wave spectra: (**a**) RMS of pitch amplitude—JONSWAP; (**b**) RMS of pitch amplitude—Torsethaugen; (**c**) RMS of vertical acceleration—JONSWAP; (**d**) RMS of vertical acceleration—Torsethaugen; (**e**) Slamming probability—JONSWAP; (**f**) Slamming probability—Torsethaugen; (**g**) Water on deck probability—JONSWAP; (**h**) Water on deck probability—Torsethaugen; (**i**) *MSI*—JONSWAP; (**j**) *MSI*—Torsethaugen.

Table 6 provides a comparative analysis between the great circle and the optimum route, together with the percentage variation of the *SPI*. Based on the current results, the RMS of the vertical acceleration and the Motion Sickness Incidence seem to be the most sensitive parameters. This outcome is confirmed with reference to both JONSWAP and Torsethaugen spectra, even if in the latter case the improvements of the seakeeping performances are slightly lower.



#### **6. Discussion**

The results obtained in Section 5, obtained by a purposely developed code in Matlab, MathWorks, clearly highlight that there is wide space to improve the ship performances in a seaway by properly varying the vessel route, depending on the expected met-ocean conditions. In this respect, the following main outcomes have been achieved:


These outcomes clearly suggest that the adaptive weather routing model shall be specialized on a case-by-case basis, depending on the ship type and the key seakeeping performances that need to be monitored and improved.

#### **7. Conclusions**

The Dijkstra shortest path algorithm was applied to detect the optimum route, constrained to maximize a non-negative cost function, namely the *SPI* index provided by Equation (1), in order to maximize the ship performances in a seaway. Based on the current results, the proposed adaptive weather routing algorithm allows for consistently increasing the ship seakeeping performances up to 50%, without significantly affecting the voyage duration or the route length. Furthermore, it was verified that the vessel speed plays a fundamental role in the assessment of the optimum route, as well as on the improvement of the ship seakeeping performances. Besides, the existence of a swell component, coming from a different direction as regards the wind wave one, can affect the selection of the optimal route. Following these outcomes, the application of an adaptive weather routing criterion, based on ship seakeeping analysis and optimization, seems to be the most suitable way to detect the optimal route, provided that the currently embodied models have an almost negligible impact in terms of fuel savings. Nevertheless, additional fuel savings, also resulting in low gas emissions in the atmosphere, could be gathered by properly varying the ship speed, depending on the estimated time of arrival in port and the availability of the harbor infrastructures. In this respect, if the ship has to wait

at anchor, due to the unavailability of the harbor infrastructures for the scheduled loading/unloading operations, it is certainly preferable to reduce the vessel speed, resulting in fuel savings, and delay the time of arrival in port. All these suggestions could be embodied by ship-owners and port authorities, in order to establish a cooperating network between the harbor and the onboard vessel management, resulting in: (i) an increase of the ship performances in a seaway with a positive impact on the safety of navigation, (ii) a reduction of the fuel consumptions, resulting in savings for the ship-owner and in a possible reduction of the gas emissions in the atmosphere, and (iii) optimization of the scheduled loading/unloading operations in port.

Based on these outcomes, the current results are encouraging for further research activities, devoted to including additional seakeeping criteria in the adaptive weather routing algorithm. In this respect, the incidence of the selected seakeeping criteria on the optimum route assessment needs to be further investigated. In fact, the *SPI* index mainly depends on the ship type and mission, which implies that the proposed adaptive weather routing model needs to be properly specialized, on a case-by-case basis. Finally, the current results were obtained in a simulated environment, and therefore, they need to be further verified in real conditions, in order to carry out a comparative analysis between the optimum route, based on the current seakeeping adaptive weather routing model, and previous real ship routes. These topics will be the subject of future works.

**Author Contributions:** Conceptualization, S.G., V.P., and A.S.; Writing—original draft, S.P. and A.I.; Writing—review and editing, S.G., V.P., and A.S.; Methodology, S.P. and A.I.; Software and formal analysis, S.P.; Data curation and visualization, A.I.; Supervision, S.G., V.P., and A.S.; Validation S.G. and V.P.; Funding acquisition, A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The work was supported by "DORA—Deployable Optics for Remote Sensing Applications" (ARS01\_00653), a project funded by MIUR-PON "Research & Innovation"/PNR 2015–2020.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The heave/pitch motion equations in the encounter wave frequency-domain are provided by the following set of equations:

$$\begin{cases} \begin{bmatrix} -\omega\_{\varepsilon}^{2}(\Lambda + A\_{53}) + i\omega\_{\varepsilon}B\_{53} + \mathbb{C}\_{33} \\ -\omega\_{\varepsilon}^{2}A\_{53} + i\omega\_{\varepsilon}B\_{53} + \mathbb{C}\_{53} \end{bmatrix} \mathcal{H}\_{3} + \begin{bmatrix} -\omega\_{\varepsilon}^{2}A\_{55} + i\omega\_{\varepsilon}B\_{55} + \mathbb{C}\_{35} \\ -\omega\_{\varepsilon}^{2}(I\_{55} + A\_{55}) + i\omega\_{\varepsilon}B\_{55} + \mathbb{C}\_{55} \end{bmatrix} \mathcal{H}\_{5} = f\_{5}(\omega) \\ \end{cases} \tag{A1}$$

where Δ and *I*<sup>55</sup> are the ship displacement and pitch moment of inertia, *Aij*, *Bij*, and *Cij* are the added mass, hydrodynamic damping, and restoring coefficient for the i–jth mode, *H*<sup>3</sup> and *H*<sup>5</sup> are the heave and pitch motion transfer functions, while *f*<sup>3</sup> and *f*<sup>5</sup> are the heave force and pitch moment per unit wave amplitude. After solving the Equation System (A1), the modulus of the heave and pitch motion transfer functions is derived [23]:

$$\left| H\_{3}(\omega\_{\varepsilon}) \right| = \left| \frac{-\omega\_{\varepsilon}^{2}(I\_{\rm BS} + A\_{\rm BS}) + i\omega\_{\varepsilon}B\_{\rm S} + C\_{\rm S} \left| f\_{3}(\omega) - \left[ -\omega\_{\varepsilon}^{2}A\_{\rm S\rm S} + i\omega\_{\varepsilon}B\_{\rm S\rm S} + C\_{\rm S\rm S} \right] f\_{3}(\omega) \right|}{\left[ -\omega\_{\varepsilon}^{2}(\Lambda + A\_{\rm S\rm S}) + i\omega\_{\varepsilon}B\_{\rm S\rm S} + C\_{\rm S} \right] \left[ -\omega\_{\varepsilon}^{2}(I\_{\rm BS} + A\_{\rm S\rm S}) + i\omega\_{\varepsilon}B\_{\rm S\rm S} + C\_{\rm S} \right] \left[ -\omega\_{\varepsilon}^{2}A\_{\rm S\rm S} + i\omega\_{\varepsilon}B\_{\rm S\rm S} + C\_{\rm S} \right]} \right| \quad \text{(A2)}$$

$$\left| \mathrm{H\_3}(\omega\_{\mathbb{E}}) \right| \ = \left| \frac{-\omega\_x^2(\Lambda + A\_{33}) + i\omega\_x B\_{33} + C\_{33} \left[ f\_3(\omega) - \left[ -\omega\_x^2 A\_{33} + i\omega\_x B\_{33} + C\_{33} \right] f\_3(\omega) \right]}{\left[ -\omega\_x^2(\Lambda + A\_{33}) + i\omega\_x B\_{33} + C\_{33} \left[ -\omega\_x^2(I\_{33} + A\_{33}) + i\omega\_x B\_{33} + C\_{33} \right] - \left[ -\omega\_x^2 A\_{33} + i\omega\_x B\_{33} + C\_{33} \right] \left[ -\omega\_x^2 A\_{33} + i\omega\_x B\_{33} + C\_{33} \right] f\_3(\omega) \right] \right| \ (\Lambda^2 + A^2)^{-1}$$

In Equations (A2) and (A3), the frequency-dependent added masses and radiation dampings combine the zero-speed values with the speed-dependent corrective factors, with no transom correction [35]:

$$\begin{aligned} A\_{33} &= A\_{33}^{0}(\omega\_{\varepsilon}); \; B\_{33} = B\_{33}^{0}(\omega\_{\varepsilon}) \\ A\_{35} &= A\_{35}^{0}(\omega\_{\varepsilon}) - \frac{\underline{\iota}\underline{\iota}}{\omega\_{\varepsilon}^{2}} B\_{33}^{0}(\omega\_{\varepsilon}); \; B\_{35} = B\_{35}^{0}(\omega\_{\varepsilon}) + \underline{\iota}\underline{\iota} A\_{33}^{0}(\omega\_{\varepsilon}) \\ A\_{53} &= A\_{53}^{0}(\omega\_{\varepsilon}) + \frac{\underline{\iota}\underline{\iota}}{\omega\_{\varepsilon}^{2}} B\_{33}^{0}(\omega\_{\varepsilon}); \; B\_{53} = B\_{53}^{0}(\omega\_{\varepsilon}) - \underline{\iota}\underline{\iota} A\_{33}^{0}(\omega\_{\varepsilon}) \\ A\_{55} &= A\_{55}^{0}(\omega\_{\varepsilon}) + \frac{\underline{\iota}^{2}}{\omega\_{\varepsilon}^{2}} A\_{33}^{0}(\omega\_{\varepsilon}); \; B\_{55} = B\_{55}^{0}(\omega\_{\varepsilon}) + \frac{\underline{\iota}\underline{\iota}^{2}}{\omega\_{\varepsilon}^{2}} B\_{33}^{0}(\omega\_{\varepsilon}) \end{aligned} \tag{A4}$$

If the 1-to-3 multivalued problem between absolute and encounter wave frequencies occurs, the modulus of heave/pitch motion transfer functions is replaced by the following equation:

$$\left| H\_{\dot{f}}(\omega\_{\varepsilon}) \right| = \sqrt{\left| H\_{\dot{f}}^{1}(\omega\_{\varepsilon}) \right|^{2} + \left| H\_{\dot{f}}^{2}(\omega\_{\varepsilon}) \right|^{2} + \left| H\_{\dot{f}}^{3}(\omega\_{\varepsilon}) \right|^{2}} \tag{A5}$$

where *H<sup>i</sup> j* (ω*e*) is the j-th mode transfer function at the absolute wave frequency ω*i*, with *i* = 1, 2, or 3.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
