A Ship Route Planning Method under the Sailing Time Constraint

Abstract

This paper realizes the simultaneous optimization of a vessel’s course and speed for a whole voyage within the estimated time of arrival (ETA), which can ensure the voyage is safe and energy-saving through proper planning of the route and speed. Firstly, a dynamic sea area model with meteorological and oceanographic data sets is established to delineate the navigable and prohibited areas; secondly, some data are extracted from the records of previous voyages, to train two artificial neural network models to predict fuel consumption rate and revolutions per minute (RPM), which are the keys to route optimization. After that, speed configuration is introduced to the optimization process, and a simultaneous optimization model for the ship’s course and speed is proposed. Then, based on a customized version of the A* algorithm, the optimization is solved in simulation. Two simulations of a ship crossing the North Pacific show that the proposed methods can make navigation decisions in advance that ensure the voyage’s safety, and compared with a naive route, the optimized navigation program can reduce fuel consumption while retaining an approximately constant time to destination and adapting to variations in oceanic conditions.

1. Introduction

With improvements in ship intelligence, big data technology, and navigation-related sensors, Maritime Autonomous Surface Ships (MASS) have attracted significant attention. At the same time, safety and energy-saving issues related to maritime navigation have gradually become the focus of attention in this new low-carbon era. In 2018, the International Maritime Organization adopted the IMOInitial Strategy for Greenhouse Gas Emission Reduction from Ships, sending a strong signal to the international community that the shipping industry is changing to a low-carbon industry [1]. Reducing fuel consumption through proper route planning is an important means to respond to the requirement for a low-carbon strategy. However, due to the length of ships’ transoceanic voyages, which can range from two weeks to more than a month, weather and sea state forecasts cannot be accurately made, and the accuracy of these forecasts decreases as the time increases. Therefore, it is required that the ship route optimization should fully consider the dynamically meteorological marine environment, and the route needs to be re-evaluated and updated when forecast data are updated [2].

Weather routing can be defined as finding an optimum route in consideration of the ETA, sailing waypoints, sailing speed, and fuel consumption based on the weather forecast data and ship performance [3]. According to this definition, some researchers mainly extract and analyze sailing records to find the association of vessels to routes [4] or summarize a data-driven optimal route in a specific navigation region [5]. Others focus on solving a path-finding problem at the operation level. The main goal of this operational problem is to plan a navigation route from one port to another based on a predefined purpose [6], in which conditions (weather, waves, and other environment variables) are along the route.

Since the ship’s navigational task often has complex influences, planning a ship route is often regarded as a pathfinding problem, which assumes that the ship’s speed or main engine power is maintained at a fixed value and only the ship’s geographical path is considered as a variable. In this way, route optimization has been researched based on various methods or algorithms, such as the modified isochrones method [7], dynamic programming [8,9,10], Dijkstra’s algorithm [11], and the A* algorithm [12], and some evolutionary algorithms such as the genetic algorithm [13,14] and the swarm algorithm [15] have also been used. However, for ocean-going ships, the meteorological and marine environment changes rapidly, and it is not easy to maintain a stable and uniform speed for a long time on the route. The above research ignores many constraints of the meteorological environment and the ship itself. Thus, it is difficult to ensure the result is the globally optimal solution [2].

In order to solve such problems, some 3D algorithms have been developed, in which the time variable is considered as the third dimension. Taking minimal fuel consumption as the goal, some scholars have used genetic algorithms to optimize meteorological routes under the premise of ensuring the ships’ safety in adverse sea conditions [16,17,18]. Evolutionary algorithms, such as genetic algorithms, have the potential to balance local searches and global searches. However, when the ocean voyage is so long that it has many feasible routes, the limited search capability in the optimization process may fail to obtain the optimal global solution. Another research interest is to add a time variable to an existing two-dimensional route searching deterministic algorithm so that it becomes a three-dimensional path searching deterministic algorithm, such as the modified three-dimensional isochrones method with weighting factors [19] and three-dimensional dynamic programming algorithms that include meteorological factors [20,21,22,23,24]. These algorithms mainly use the idea of staged optimization in dynamic programming. In each stage of the optimization process, only a few nodes with the best results can enter the next stage, while the other sub-path nodes are eliminated to reduce the computational effort [6]. Meanwhile, in each stage of the recursion process, the discarded nodes may be associated with the optimal global solution with deterministic constraints. In contrast, the A* algorithm has the advantage of retaining all nodes in each recursion and finding the optimal nodes in stages; thus, it increases the possibility of finding the optimal global solution. However, it also results in an exponential increase in the number of search nodes with each stage of the recursion process, which makes it difficult to add a time-like dimension to the route optimization problem.

With the continuous development of artificial intelligence technology, its use to achieve ship weather route optimization has become a new research interest. The idea is to generate recommended routes after summarizing and analyzing a large amount of ship navigation data. Some scholars process AIS data to recognize key turning regions and connect these turning regions via cluster similarity measuring to generate reasonable routes for different types of ships [25,26]. Others try to use artificial neural networks and machine learning algorithms to predict the fuel consumption of a ship under different sailing conditions to achieve the goal of determining the expected duration or saving energy [27]. For example, artificial neural networks have been used to optimize a ship’s speed based on large amounts of ship operation data to reduce fuel consumption [28,29]. Three different statistical models have been used to forecast and optimize the speed of container ship routes with the goal of improving navigational safety and determining expected duration by integrating meteorological information such as wave height, wave period, wind speed, and other information [30]. From the above research, it can be shown that different weather route optimization methods often require different applicable conditions and generally need to assume some ideal conditions, such as a more stable ship navigation state, fewer optimization nodes or dimensions, etc.

For a voyage’s safety and economy, ship course and speed are the two key variables in saving energy and controlling arrival time. In an actual navigation situation, distinct course and speed choices may cause a vessel to encounter radically different sea states, so a ship’s course and speed should be planned according to the dynamic sea environment, and most commercial software do so, such as the Bon voyage system, Seaware enroute, SPOS, and so on [31]. Normally, some environmental factors, like the resistance of calm water [32], waves [33], and winds [34], can significantly affect a ship’s speed. Calculating the values of ship dynamics is extremely complex because of the consideration of hull shape, seakeeping characteristics of the ship, the sea spectrum, and other parameters. Thus, it is too laborious and time-consuming for ship weather routing. In addition, under severe sea conditions, fuel consumption can be reduced exponentially by properly reducing a ship’s speed, thus the ship can adjust its speeds in different sea states to reduce its fuel consumption [27].

Therefore, this paper proposes a ship route planning method under a sailing time constraint, which introduces speed as a variable in the ship route optimization process to realize the co-optimization of ship course and speed. In the next section, a three-dimensional sea environment model is constructed, which can load weather data in real time to inform the route optimization process. Section 3 builds an ANN-based model that relates ship parameters to fuel consumption rate and RPM for use in the optimization process. Section 4 introduces a speed variable to the route optimization process to jointly optimize ship course and speed, to design the voyage route, and to plan for deceleration or acceleration to effectively reduce fuel consumption. As a result of adding a variable, it can be understood from the above description of the 3D deterministic algorithms that the number of sub-paths increases exponentially as the algorithm iterates [35,36]. For ship route planning oriented towards transoceanic voyages, this will increase the searching time significantly. Thus, to compute a navigation strategy under the sailing time constraint in an acceptable computational time, Section 5 proposes a customized A* algorithm, which uses ETA to constrain the search space and a suitable heuristic function to guide the direction of the search. Two simulations are performed in Section 6 to test and highlight the main features of the method, followed by conclusions in Section 7.

2. Dynamic Sea Environment Model

2.1. Three-Dimensional Sea Environment Model

In current navigation situations, route optimization should comply with the following restrictions:

• The ship should navigate in an area with sufficient water depth.
• The ship should not sail in dangerous wind and wave conditions.
• The route plan should consider dynamic meteorological and sea conditions for long-distance and long-term navigation.

To satisfy the above constraints, a three-dimensional sea environment model is established for ship route optimization. Most of the forecast data and reanalysis data of the current meteorological and sea state are in “GRIB” format, and this type of data is used to divide the ocean area into a sequence of dense grids. Each grid cell has uniform meteorological and oceanographic properties; hence, grid cells can be taken as the route nodes and discretize ship routes into sequences of waypoints.

The environmental grid data in this paper include meteorological and oceanographic data from the European Centre for Medium-Range Weather Forecasts [37] and ocean bathymetry data from the General Bathymetric Chart of the Oceans [38]. Each two-dimensional grid datum at a single time is read from the environmental data and recorded in the form of a matrix. These multiple time matrices are stored in a 3D matrix, as shown in Figure 1, where $t_1,t_2,\cdots t_n$ represent the meteorological and oceanographic information in the grid cells at discrete times t, and the time interval between them is generally fixed. Moreover, the environmental data in row x and column y at time t are defined by $E(x,y,t)$. The environmental effect on a ship can be expressed as follows:

$$E(x,y,t)=f(\phi ,\lambda ,v,t_0)$$

(1)

Here, $(\phi ,\lambda )$ is the ship’s location, v is the ship’s speed, and $t_0$ is the ship’s departure time. After the calculation of the waypoints sequence and the corresponding speed configuration from the departure time $t_0$, the goal grid $(x,y)$ will be found to arrive at time t. In addition, the index $(x,y)$ is also related to the range and accuracy of the environmental data.

2.2. Binarization of Navigation Area

After reading and storing the environmental data, the data representation is shown in Figure 2. The white contour lines in the figure are contours of wave height, the direction of the arrows indicates the wind direction, their length indicates the value of the wind speed, and the elevation is represented by the isosurface.

As can be seen in Figure 2, there are obstructive conditions in the map, including static obstacles such as islands and shoals in the marine environment, and dynamic restricted navigation zones corresponding to areas of strong winds and high waves (e.g., the contour-dense area in Figure 2 has wave heights above 6 m). Grid cells with obstructive conditions should be marked as unnavigable cells through which the route is forbidden to pass, while others marked as navigable, in which the ship can sail through unconstrained. Significantly, the classification and marking of the cells should be pre-processed considering the ship’s characteristics and performance factors.

A binary method is used in grid pre-processing to divide the navigation area, and the navigable cells are set to 1; unnavigable cells are set to 0, which are also called restricted grids. The grid pre-processing is shown in the Formula (2). Combining the ship draft and elevation information, when the water depth of the grid $(x,y)$ is not sufficient to ensure the safety of shipping, it is set as unnavigable and expressed as $Navi(x,y)=0$. Also, when the wind and wave of the grid $(x,y)$ exceeds the ship’s anti-wave ability at time t, it is set as unnavigable and expressed as $Navi(x,y,t)=0$.

$$Navi(x,y,t)=\left\{\begin{array}{cc}0,\hfill & insufficientwaterdepth,\hfill \\ 0,\hfill & windandwavesoverthethreshold,\hfill \\ 1,\hfill & meettheshi{p}^{\prime }snavigationconditions.\hfill \end{array}\right.$$

(2)

In this way, the marine environment in Figure 2 can be binarized, and the restricted grids are filled with black while the free grids are filled with white. Thus, the restricted and free grids can be distinguished. And the accuracy of grids is often the same as the acquired environment data, which is set to $0.5°\times 0.5°$ in this paper, as shown in Figure 3. In this figure, a sea area with wave heights over 6 m is set as a unnavigable area, as framed in the red box.

3. Predictive Model for Ship Parameters

Since the sea status during the voyage is not always static, the relevant parameters of the ship will constantly change (e.g., fuel consumption rate, speed, course, etc.), which affects the corresponding sailing decisions. Therefore, establishing a reliable and effective model to predict these effects is necessary for accurate route optimization. Artificial neural networks can be used for modeling complex systems. They are highly adaptable, robust, fault-tolerant, and surface-fitting and are well-suited to predictive analysis. Since the prediction model needs to be continuously invoked in route planning, a prediction model with a relatively simple structure and high prediction accuracy is required, and an artificial neural network model with a single hidden layer can meet such needs [39].

3.1. Pre-Processing of Data

To accurately predict the ship’s sailing state under various external factors, a large amount of sailing data and an artificial neural network-based ship parameter prediction model are needed. In this work, the training and testing data are from actual navigation records of one container ship’s voyages from 1 January to 28 February 2021. The details of that container ship are shown in

Table 1. Ship parameters.

Parameter	Value
Ship type	Container ship
Displacement (t)	169,700
Length (m)	348
Beam (m)	51.2
Draft (m)	13.5
Pitch (mm)	8668

.

A total of 36,440 consecutive records were collected during the voyages. Each record has required data types, such as date, latitude, longitude, speed over ground (SOG), course, draft, fuel consumption, and RPM, etc. Due to errors in equipment and instruments, network signals, and other factors, there are noisy and erroneous data in the records. The predictions of the model may have large errors compared with the actual values if all records are directly used for model training, making it difficult to reflect the actual fuel consumption under many influencing factors. So, a data cleaning method proposed by ISO based on the idea of removing track anomalies was designed to filter the data and eliminate useless data [40]; the major steps are shown below:

• Remove anomalous records in which data are incomplete, such as missing latitude and longitude data.
• To reduce the complexity of the artificial neural network mentioned later, the absolute directions of weather are converted to the directions to the ship reference frame:

  $$di{r}_{wr}=|180-|(di{r}_c-di{r}_{wa})mod360||$$

  (3)

  where $di{r}_{wr}$ represents directions of weather components (wind and wave) to the ship’s reference frame, $di{r}_c$ is the ship’s course, and $di{r}_{wa}$ is the absolute directions of weather components.
• According to this method [40], extract sequential records of approximately 10 min and combine them into one segment.
• Calculate the mean value $\mu$ of each type of data $d_i$ in a segment.
  For data that are not measured as angles, the mean $\mu$ for the N data points in a segment with values $d_i$ is computed by

  $$\mu =\frac{1}{N}{\sum }_{i=1}^N{d}_i$$

  (4)
  For data that are measured as angles, the mean $\mu$ is computed by

  $$\mu =arctan\left(\frac{{\sum }_{i=1}^{N}sin{d}_i}{N}/\frac{{\sum }_{i=1}^{N}cos{d}_i}{N}\right)$$

  (5)
• Calculate the standard error of the mean $\sigma$ of each type of data in a segment. The standard error of the mean $\sigma$ is computed by

  $$\sigma =\sqrt{\frac{1}{N}{\sum }_i^N{\Delta }_i^2}$$

  (6)
  For data not measured as angles, the difference ${\Delta }_i$ is computed by

  $${\Delta }_i=|d_i-\mu |$$

  (7)
  For data measured as angles, the difference ${\Delta }_i$ is computed by

  $$\left\{\begin{array}{cc}{\Delta }_i={360}^{\circ }-r_i,\hfill & r_i=mod\left(\right|d_i-\mu |,{360}^{\circ })>{180}^{\circ }\hfill \\ {\Delta }_i=r_i,\hfill & r_i=mod\left(\right|d_i-\mu |,{360}^{\circ })\le {180}^{\circ }\hfill \end{array}\right.$$

  (8)
• Two thresholds are used to delete segments in which the standard error of the mean of certain data is too large. In particular, if the $\sigma$ of RPM is greater than 3 min${}^{-1}$ or the $\sigma$ of SOG is greater than 0.5 kt, then all records in that segment are deleted. After this operation, several uniform sailing segments are obtained.
• We use Chauvenet’s criterion to determine if data are anomalous in each uniform sailing segment. Once a data point is identified as an outlier based on the criterion, it is removed. A new mean and standard error of the mean can be calculated based on the remaining values and the new sample size.
  The probability for the occurrence of any $d_i$ is computed by

  $$P\left(d_i\right)=\mathrm{erfc}\left(\frac{{\Delta }_i}{\sigma ·\sqrt{2}}\right)$$

  (9)

  where $P\left(d_i\right)$ is the probability of occurrence of $d_i$ and $\mathrm{erfc}$ is the complementary error function.
  A datum is considered an outlier if Formula (10) is fulfilled.

  $$P\left(d_i\right)\times N<0.5$$

  (10)
  This method cleans the original voyages’ records to obtain several navigation segments with a duration of about 10 min and a uniform speed. The mean values of these segments’ data can be used for subsequent training and prediction of the model.
• The different types of input data have different units. It would reduce the performance and convergence of the predictive model if the filtered data of the segments were input into the neural network directly. Therefore, the z-score standardization method is adopted to standardize the data:

  $$d_{ty}^{*}={\displaystyle \frac{d_{ty}-{\mu }_{ty}}{{\sigma }_{ty}}}$$

  (11)

  where $d_{ty}^{*}$ is one type of data in a segment after the standardization, $d_{ty}$ is this type of data unprocessed in the segment, and ${\mu }_{ty}$ and ${\sigma }_{ty}$ are the mean and standard error value of this type of data in all segments.

3.2. Predictive Model for Ship Parameters

An artificial neural network (ANN) has good nonlinear mapping capability, adaptive learning capability, and parallel information processing capability. Compared with traditional system identification methods, ANNs do not require knowledge of the physical causal relationships between the observed system variables, providing a way to model the system without information about the internal state of the system. In this paper, an ANN with a single hidden layer is used as a prediction tool for fuel consumption rate and RPM because of its simplicity, robustness, and ideal prediction effect [39,41]. In addition, the ship’s sailing performance also changes as time passes, and the ANN can be reconstructed quickly based on more recent sailing data when necessary.

3.2.1. Structure of the Predictive Model

An artificial neural network consists of several connection layers, which can be divided into an input layer, hidden layer, and output layer according to their position and function. All neurons in the network are connected to adjacent neurons located in different hierarchies. The structure of an ANN with a single hidden layer is shown in Figure 4.

Here, $w_{ij}$ is the connection weight from the ith neuron in the input layer to the jth neuron in the hidden layer. $v_{jl}$ is the connection weight from the jth neuron in the hidden layer to the lth neuron in the output layer. ${\theta }_j$ is the bias of the jth neuron in the hidden layer. ${\gamma }_l$ is the bias of the lth neuron in the output layer. f is the activation function of the hidden layer, and g is the activation function of the output layer.

An ANN establishes a mapping from N-dimensional to Q-dimensional data through the above process, and the functional relationships can be expressed by the following formulas:

$$b_j=f\left({\sum }_{i=1}^n{w}_{ij}{a}_i+{\theta }_j\right)$$

(12)

$$c_l=g\left({\sum }_{j=1}^p{v}_{jl}{b}_j+{\gamma }_l\right)$$

(13)

The activation function of the hidden layer $f\left(x\right)$ must be a bonded nonlinear function that is continuous, smooth, and monotonically increasing. So, the Relu activation function is used for its advantages of simple computation, simple derivatives, fast convergence, one-sided inhibition, and broader boundaries. The activation function of the output layer $g\left(x\right)$ does not necessarily need to be nonlinear, so a linear function is used. The relevant activation functions are shown as follows:

$$f\left(x\right)=max\left(0,x\right)$$

(14)

$$g\left(x\right)=x$$

(15)

3.2.2. Error Back-Propagation Algorithm

The error back-propagation algorithm [42] is a method that uses the difference between the actual output and the desired output to correct the weight parameters layer by layer from back to front, which applies to supervised learning. When learning samples are provided to the network, if there is a gap between the desired output and the output based on current weights and biases, the error function of the current output and the desired output is calculated, and then the error signal is back-propagated along the original connection path to update the neuron weight parameters of each layer.

Assuming that M learning samples are provided to the network, the error function of the network for each learning sample is constructed as follows.

The goal of learning is to minimize a given error function, and network training is generally stopped when the error function value terminates a downward trend and starts to rise:

$$E^k=\frac{1}{2}{\sum }_{l=1}^P{\left(y_l^k-c_l^k\right)}^2$$

(16)

where k is the learning sample number ($k=1\sim M$), l is the dimension of the output layer ($l=1\sim Q$), $E^k$ is the error function of the kth learning sample, $y_l^k$ is the expected output in the lth dimension of the kth learning sample, and $c_l^k$ is the actual output of the lth dimension of the kth learning sample.

The network calculates the global error after obtaining all M samples and updates the parameters uniformly based on the global error E. The global error E of the network is the sum of all the learning sample errors:

$$E={\sum }_{k=1}^M{E}^k$$

(17)

This work uses the Levenberg–Marquardt (LM) algorithm to update the parameters. The LM algorithm combines the advantages of the gradient descent method and the Newton method, as its error function decreases faster at the beginning, and provides the ideal search direction when the error function is near the optimal value [43]. The parameter updating formula of the LM algorithm is as follows:

$$\omega \left(i+1\right)=\omega \left(i\right)-{\left(H+\lambda I\right)}^{-1}\frac{\partial E}{\partial \omega }$$

(18)

Here, $\omega$ represents the neural network parameters, including connection weights between different layers and the offsets of each layer; i is the number of learning times; H is the Hessian Matrix of the error function; and I is the unit matrix.

The parameter $\lambda$ is updated by taking larger values in the early stages of learning so that the parameters are updated along the inverse direction of the gradient, while the direction of iteration of the step is shifted to the direction of the Newton method as $\lambda$ gradually decreases to zero in the later stages [43]. Compared with other iterative methods, the LM algorithm has a relatively good convergence rate for training networks of medium size [44].

3.3. Selection of Model Variables

The two essential criteria for evaluating the utility of a sailing route are sailing time and fuel consumption. Except for voyage distance, the variables most closely related to these two criteria are the speed over ground (SOG) and fuel consumption rate. SOG is affected by the navigation environment in a natural sea state, and differs from the calm water speed due to added resistance from waves, winds, and currents, as well as the reduction in propulsive efficiency caused by waves and increased resistance, which is often known as ship speed loss. Ship speed loss is a critical parameter in route planning for voyage time and fuel consumption evaluation, and it is calculated based on the calm water speed and the actual sea state in the general route optimization process; then, the ship’s actual SOG is obtained. But in voyage records, it can be noted that the expected calm water speed is not recorded, so the traditional way to calculate speed loss using calm water speed and the SOG in natural sea state is no longer practical. Fortunately, there is the real-time RPM data in the voyage record, which is a measurement closely related to the expected calm water speed.

So, the actual SOG is taken as one optimization variable in this work, and once the SOG is assigned, the sailing time can easily be calculated, and the recommended value of RPM can be predicted through the ANN model, as shown in Figure 4, which can make the ship sail at the assigned SOG. In this way, the output variables of the ANN model are RPM and fuel consumption rate, and they are used to calculate the total sailing time and fuel consumption in the ship route optimization procedure.

To meet the route optimization task, the calculation and prediction of ship-related parameters, which generally need to be updated frequently to serve the route planning algorithm, should not occupy too much running time. Thus, the input variables of the ANN should have the following characteristics:

• They can represent the current navigation state or navigation environment of the ship.
• They can be obtained or predicted in some way in the route planning process.
• How the data are obtained or predicted cannot be too complex or take up too much running time.

Thus, eight input variables which meet the above characteristics have been selected to establish an ANN-based model for fuel consumption rate and RPM prediction, including ship data such as the SOG, ship course, fore draft, aft draft, and sea environment data such as relative wind direction, relative wave direction, wind speed, and wave height. The relevant variable mapping relationship is shown in Figure 5.

4. Co-Optimization of Ship Course and Speed

The route between the starting node and end node consists of several rhumb lines, and each rhumb line connects two adjacent nodes, as shown in Figure 6. Each node includes three variables: position, time, and speed. Then, a sequence of n nodes can represent the whole voyage, including the speed along each leg.

The course is represented by the orientation of the rhumb line connecting two adjacent nodes; by adjusting the position of the nodes, the course is changed. The time when the ship arrives at the next node is determined by the course and speed, and the meteorological and sea conditions that the ship encounters vary as time goes by. So, the route and its speed are the keys to ship route decision, and by adjusting the position of nodes and the speed between them, the whole route can be optimized.

Figure 7 shows the process of route generation. The centers of grids marked with 83, 67, 45 indicate nodes where the ship is located; the surrounding arrows refer to alternative courses from one node to the next; and $t_1$, $t_2$, $t_3$ are the times when the ship reaches these nodes in sequence. When the ship chooses different courses and sails with different speeds, the nodes and their arrival times are different too, so the ship will encounter different wind and wave states, and the sailing time and fuel consumption may differ accordingly. The wind and wave state the ship encounters can be obtained using the 3D sea environment model in Section 2.1, while the fuel consumption can be calculated using the predictive model in Section 3.

Thus, node selection and speed configuration are the two key methods in ship voyage optimization, and optimizing ship course as well as speed is an efficient way to plan a voyage, so a method for the co-optimization of a ship’s course and speed is proposed in this work.

4.1. Speed Configuration

Usually, the speed during the voyage will be maintained in a range according to the ship’s performance. To ensure safety and economic efficiency, the ship will sail at a suitable speed in a certain sea area, and the set of permissible speeds is presented as follows:

$$v=\left[v_{min},v_{min}+v_{cd},v_{min}+2v_{cd},\dots ,v_{max}\right]$$

(19)

Here, v is a row matrix vector including n permissible speeds, its interval is $v_{cd}$, $v_{min}$ is the minimum value in v, and $v_{max}$ is the maximum value in v.

Accordingly, since there are multiple permissible speeds in each node, a leg from one node to an adjacent node is expanded from 1 to n. Alternatively, one node i is extended to n “ext-nodes” $i_w$, described as follows:

$$i_w=\left\{P\left(i,v_w\right)|v_w\in v\right\}$$

(20)

Here, $i_w$ means the ship sails from node i to the next adjacent node with fixed speed $v_w$, and P contains all navigation states at node i. In addition, one sailing node is expanded into a row matrix. Thus, this method can be used to achieve parallel computation and speed up the algorithm operation.

Therefore, a speed configuration can be a part of an optimized route, as shown in Figure 8. The expression form of the optimized route from the starting node S to the end node T evolves from the original $S-A-B-C-D-E-T$ to $S_{13}-A_{12}-B_{13}-C_{14}-D_{12}-E_{15}-T$. That is, the ship sails from the origin node S to node A at 13 kt, from node A to node B at 12 kt, and continues sailing until from node E to the end node T at 15 kt. With this method, the optimized route not only contains the optimized route (Figure 9) but also includes the optimal speed in each segment (Figure 10) to realize the synergistic optimization of the ship’s course and speed.

5. Weather Route Optimization for Ships with Sailing Time Constraints

5.1. A* Algorithm

The A* algorithm is based on Dijkstra’s algorithm, with an added heuristic function for predicting the future cost of movement to influence the search direction and narrow the search space, so that it can find the shortest path in a certain sense in a short time.

The main process of the A* algorithm is to build an “$Openlist$” node set to store nodes to be checked and a “$Closedlist$” node set to store checked nodes, and iterate continuously to select the node with the lowest movement cost to obtain the shortest path. The basic form of the function to evaluate the movement cost is shown in the following formula:

$$F\left(i\right)=G\left(i\right)+H\left(i\right)$$

(21)

Here, i is the current node. F is the estimated total movement cost from the starting node to the end node. G is the minimum movement cost recorded from the starting node to the current node, and the heuristic function H is the estimated movement cost from the current node to the end node.

The A* algorithm calculates and iterates the node selection according to this evaluation function. Its basic process is as follows: the total movement cost of surrounding nodes is calculated from the starting node. The node with the minimum movement cost is selected in each iteration until the end node is found. The process is shown in Figure 11.

The A* algorithm has characteristics of fast calculation speed to an optimal solution, but there are several problems that need to be solved for route optimization tasks, as follows:

• The two-dimensional A* algorithm takes distance as the evaluation criterion. At the same time, ship fuel consumption and navigation time are the main factors to evaluate the navigation cost of a route, and taking distance as the evaluation criterion cannot accurately quantify the navigation cost.
• The two-dimensional A* algorithm does not consider speed. In contrast, the ship does not sail at a fixed speed over the whole voyage; dynamically adjusting speed according to the sea state is desirable.
• The two-dimensional A* algorithm does not consider dynamic unnavigable areas with heavy winds and waves, which will affect the safety and efficiency of the voyage.

Thus, the A* algorithm needs to be modified to meet the actual needs of ship navigation.

5.2. Evaluation Functions Concerning Time and Fuel Consumption

The Estimated Time of Arrival (ETA) is the time when a ship or vessel is expected to arrive at a specific destination. An accurate ETA can make entire supply chains more efficient and reliable, and it also helps to determine the expected duration of a vessel’s route.

Generally, as the two most important factors in evaluating a route, voyage time and fuel consumption need to be fully considered in the evaluation function. In fact, the voyage schedule is the crucial constraint for merchant vessels, and in the scope of reducing gas emissions, the effect of an increased sailing speed could be enhanced by the waiting time for a free berth, therefore accounting for even more emissions. The ship must arrive at the specific destination at a time close to the ETA, and thus the expected duration is taken as a constraint in the route design process, and minimum fuel consumption is the optimization objective.

In order to meet the above description, two evaluation functions for sailing time and fuel consumption are proposed:

$$F{T}_{iw}=G{T}_{iw}+H{T}_{iw}$$

(22)

$$F{F}_{iw}=G{F}_{iw}+H{F}_{iw}$$

(23)

Here, $G{T}_{iw}$ and $G{F}_{iw}$ are the actual time and fuel consumption of the route from the starting node to the current ext-node $i_w$, described as follows:

$$G{T}_{iw}=\sum _{i=1}^n\frac{d_i}{v_w}$$

(24)

$$G{F}_{iw}=\sum _{i=1}^n\frac{d_i}{v_w}\times q_{iw}$$

(25)

Let $h_w$ and $i_w$ describe two adjacent ext-nodes and $h_w$ be the previous node; then, $d_i$ is the distance between $h_w$ and $i_w$, and $\frac{d_i}{v_w}$ refers to the sailing time from $h_w$ to $i_w$ at the speed $v_w$. $q_{iw}$ refers to the fuel consumption rate when the ship sails from $h_w$ to $i_w$. The value of $q_{iw}$ can be calculated using the fuel consumption rate model in Section 4.1.

In addition, it is well-known that since the heuristic function H affects the search efficiency of the A* algorithm, a suitable H function can guarantee a faster search speed and the quality of the route [45]. To ensure the solution quality with a short running time, heuristic functions with suitable search efficiency, named $HT$ and $HF$, are proposed as follows:

$$H{T}_{iw}=d{E}_i/v_w$$

(26)

$$H{F}_{iw}=\frac{d{E}_i}{v_{eco}}·f_q\left(v_{eco}\right)$$

(27)

where $d{E}_i$ is the distance from node i to the end node; $H{T}_{iw}$ is the estimated time when the ship reaches the end node from the current ext-node $i_w$ at the speed $v_w$; $H{F}_{iw}$ is the estimated fuel consumption from the current ext-node $i_w$ to the end node at speed $v_{eco}$; $v_{eco}$ is the economic speed of the ship, usually artificially set; and $f_q$ is a function of the fuel consumption rate and speed.

The relationship of fuel consumption rate with speed is generally a cubic relationship, as shown in the Formula (28):

$$f_q\left(v\right)=a{v}^3+b{v}^2+cv+d$$

(28)

where a, b, c, and d are fitting coefficients that vary with the ship’s characteristics.

5.3. Customized A* Algorithm with Time Constraint

After reconstructing the evaluation function, the process of the A* algorithm is also customized accordingly. The process is as follows, which is shown in Figure 12:

• Initialize $Openlist$ (the set of ext-nodes to be checked) and $Closedlist$ (the set of ext-nodes checked).
• All ext-nodes belonging to the start node S are added to the $Openlist$, and their $GT$ and $GF$ are set to 0, while the $HT$ and $HF$ are calculated and recorded.
• If it is judged that the $Openlist$ is not empty, the operation goes to Step 4; otherwise, there is no solution to this problem, which means no available path exists between S and T, and the procedure is terminated.
• The total navigation time $FT$ of each ext-node in the $Openlist$ is calculated. If the total navigation time $FT$ of one ext-node is greater than the pre-set ETA, the node is retained in the $Openlist$ but will not be executed in the rest of this step. Then, the ext-node $i_w$ with the minimum fuel consumption of the whole voyage is found in these nodes which meet the time limit. And the ext-node $i_w$ will be removed from the $Openlist$ and added into the $Closedlist$, setting the SOG record to $v_w$.
• Traverse all the adjacent nodes of the navigation node i. There are two situations for the ext-nodes of these adjacent nodes:
  • If there is no ext-node of the adjacent node j in $Openlist$ and $Closedlist$, all ext-nodes of node j are added to $Openlist$, the G and H value of these ext-nodes are calculated and recorded, and the parent node of these ext-nodes is set as node $i_w$.
  • If there is an ext-node $j_w$ of the adjacent node j in the $Openlist$, and the $GT$ and $GF$ values of $j_w$ calculated in this step are both smaller than those stored in the $Openlist$, the G and H values of $j_w$ in the $Openlist$ will be updated correspondingly in this step. Also, the ext-node $j_w$’s parent node is set as node $i_w$. Conversely, all information of the ext-node $j_w$ stored in $Openlist$ will remain unchanged.
• Check whether there is any ext-node in $Openlist$ that is subordinate to the end node T. If not, repeat Step 3; if it already exists, a path composed of ext-nodes is found by querying the parent node of the current ext-node step by step until the starting node. The sailing route and speed configuration are obtained after mapping like in Section 4.1.
• The weather forecast data and navigation data at each ext-node on the route are extracted based on the waypoint sequence and speed configuration. Then, the RPM prediction model is called, and the above data are modified as the input of the prediction model to calculate the recommended RPM of each ext-node on the route. Then, the RPM recommendation scheme is output for the whole voyage.

5.4. Overall Optimization Flow

The weather route optimization method proposed in this paper consists of four parts: ship historical data processing, the ship parameter prediction model, ship navigation data processing, and ship route optimization. Briefly, the work flow is as follows, which is shown in the Figure 13:

• The historical data processing part uses historical ship voyage records and corresponding meteorological historical data, etc. These data are standardized by pre-processing data methods, and the pre-processed data are used as input for the ship parameter prediction models.
• The ship parameter prediction model part is used to build two artificial neural network models with the same frame. The input contains eight variables: SOG, fore draft, aft draft, course, wind speed, wind direction, wave height, and wave direction; the output is fuel consumption rate or RPM. The trained models are used for calculating fuel consumption rate and RPM in the voyage afterward.
• The ship navigation data processing part obtains and divides the ship’s characteristics and the weather forecast data into eight input variables which meet the requirements of the ship parameter prediction model and are used in the route planning process afterward.
• The ship route optimization part is based on the customized A* algorithm, using the ship fuel consumption rate prediction model and the current sailing data to calculate the fuel consumption for all nodes to be checked in order to find the node with the minimum fuel consumption that satisfies the ETA constraint and thus obtain the optimal ship route. Finally, the RPM recommendation scheme for the whole range is obtained based on the sailing data and the generated route.

Notably, the ship route can be frequently updated during the voyage. The frequency of calculation during the sailing voyage would depend on the fleet’s will and the updating frequency of the forecast weather data. If the voyage is too long to acquire the forecast data at the later part of the voyage, this method would find the optimum route based on the current weather information. In addition, the updated route may have great difference with the prior one. But both routes are planned based on the assurance of sailing safety and economy, which can be identified as the optimal sailing strategy at their respective states.

6. Simulation Experiment and Results Analysis

6.1. Performance Analysis of Ship Parameter Prediction Models

In order to evaluate the prediction effect of the artificial neural network model designed in this paper, the data set is divided into a training set and a test set. The training set is used to train the model, and the test set is used to test the generalization ability of the model. Each ten-minute segment (extracted in Section 3.1) is taken as a data point by calculating its average value. A total of 1976 data points are obtained from the data set, of which 1580 data points are randomly picked and used as the training set, and the remaining 396 data points are used as the test set. The parameters of fuel consumption rate model and RPM model are shown in

Table 2. ANN parameters.

Parameters	Value
Number of neurons in hidden layer of fuel consumption rate model	45
Number of neurons in hidden layer of RPM model	140
Activation function of the hidden layer	Relu
Optimizer	Levenberg–Marquardt

. The prediction effect of the two models’ test sets are shown in Figure 14 and Figure 15, and the blue dots are the actual value of the test set.

A perfect predictive model should have a predicted value equal to the true value, which means the points would lie on the black diagonal line. The vertical distance from the line to any point is the prediction error for that point. A good model has minor errors, as Figure 14 and Figure 15 show; all the blue dots surround the black lines, and the average relative error of the fuel consumption rate model and RPM model are 0.0594 and 0.0101, respectively. Thus, the above models can meet the navigation requirements in practice and can be used in the ship route optimization process to calculate fuel consumption and RPM.

6.2. Simulation Experiment of Ship Weather Route Optimization System

The proposed ship route optimization method under the sailing time constraint, described in detail in Section 5, has been validated and verified in a series of computer simulations based on actual voyage records. Two scenarios have been selected with different passages (voyage departure and destination points) and times (different weather conditions, having a direct impact on the optimization process). The simulations are programmed using MATLAB, and the accuracy of the relevant meteorological data is ${0.5}^{\circ }\times {0.5}^{\circ }$, as mentioned in Section 2.

6.2.1. Scenario 1: Voyage from Osaka to Los Angeles on 3 February 2021

In this scenario, a voyage from Osaka with a departure on 3 February 2021, at 18:00 UTC with a destination of Los Angeles has been simulated; its expected duration is 370 h, and its economic speed is set to 13.5 kt. The primary purpose of this simulation is to verify the ability to avoid wind and wave areas.

The experimental parameters are listed in

Table 3. Experimental parameters.

Parameters	Value
Start	35.5${}^{\circ }$ N, 141${}^{\circ }$ E
End	34${}^{\circ }$ N, 120${}^{\circ }$ W
SOG range (knot)	[8, 19]
Maximum allowable wind speed (m/s)	20
Maximum allowable wave height (m)	6
Minimum draught (m)	14

, in which the threshold values are set according to the actual demand of the ship’s voyage.

The ship’s course and speed were optimized using the methods proposed in this work, and an optimal route was generated. As shown in Figure 16, the optimal route (red curve) is compared with a great circle (GC, green curve) with SOG of 13.5 kt for the whole voyage. For Figure 16, the wave heights are represented by the isosurface and the wind speeds are represented by contour lines. The four charts describe the wind and wave conditions and the routes at different times (80 h, 165 h, 230 h, and 367 h, referring to the time after the ship departs).

From

Table 4. Comparison of routes.

	Time (h)	Danger Time (h)
Optimized route	367.15	0
Great circle route	361.38	26.26

, it can be concluded that when the ship sails along the great circle route, there are 26.26 h of sailing time in heavy wind and wave areas, accounting for about 7.27% of the total sailing time, which is very unfavorable to the ship’s navigational safety, while the optimized route generated by the proposed method in this work avoids the dangerous sea area with heavy wind and waves.

The distribution of and variation in the wind and waves during the voyage are demonstrated in Figure 16, where “80 h” refers to the situation that the optimized route chooses to sail to higher latitudes in advance based on the analysis of meteorological information; at about 165 h and 230 h, there are heavy winds and waves on the great circle route, and the optimized route avoids this area in advance; “367 h” is the comparison of the two routes after reaching the destination point.

The speed distribution is shown in Figure 17, and the wave heights on the two routes are shown in Figure 18. Define the route between each two adjacent nodes as a “segment”, and the segment numbers from the starting node to the end node increase in sequence.

It can be seen from Figure 18 that the optimized route encounters a significantly lower incidence of high significant wave heights compared with the great circle route, and the maximum wave height is also smaller. Therefore, the algorithm can handle the upcoming navigational risks in advance. The speed configuration is also more reasonable than the case of uniform speed for the whole voyage.

6.2.2. Scenario 2: Voyage from Tainan to Los Angeles on 17 January 2021

In this scenario, a voyage from Tainan (23${}^{\circ }$ N, 120${}^{\circ }$ E) with a departure on 17 January 2021 at 07:00 UTC with a destination of Los Angeles (34${}^{\circ }$ N, 120${}^{\circ }$ W) was simulated, and its expected duration is 410 h. The primary purpose is to compare the optimized route with the historical route to evaluate the effectiveness of the algorithm. Some experimental parameters of the ship are the same as in Scenario 1.

As shown in Figure 19, the optimized route (red curve) is compared with the historical route (green curve). The economic speed of the optimal route is set to 13.5 kt. It is known that the ship’s safety was not compromised on the historical route, so the comparison focuses on the sailing time and fuel consumption (FOC).

As seen from

Table 5. Comparison of routes.

	Time (h)	FOC (t)
Optimized route	408.69	1674.27
Historical route	404.34	1805.20

, compared with the historical route, the optimized route saves about 130.93 tons of fuel consumption (7.25%) within the specified sailing time. And the increase in voyage time is by 4.35 h (1.08%), which is tolerable for the route planning. Thus, it achieves the purpose of saving energy by designing the optimal route and reasonably configuring speeds. In addition, the speed configuration (Figure 20) and required RPM (Figure 21) for the whole voyage are also generated. Moreover, the output RPMs of the whole voyage will also guide the ship to reach the planned SOGs, which provides more assistance for real-time ship navigation.

6.2.3. Computational Complexity of This Method

This paper adds one dimension about the SOG to the route planning method, which may cause a exponential increase in the computational cost. But with the rational design of the customized A* algorithm’s evaluation function and operation process, the computational cost is significantly reduced, as

Table 6. Computational cost of two scenarios.

	Computational Time (s)	Data Size (MB)
Scenario 1	991.917	1.10
Scenario 2	731.639	1.04

shows. When the optimal route is output, the computation times for the two scenarios are 991.917 s and 731.639 s, which is tolerable in the practical application. In addition, the storage sizes of the data (includes $Openlist$ and $Closedlist$) are 1.10 MB and 1.04 MB. Given all the above, the computational complexity of this method is acceptable.

7. Conclusions

As the main influencing factor of navigation efficiency in the process of ship route optimization, speed control gives ample space for optimization and should be fully considered in ship route design. Based on this, this paper proposes a ship route planning method under a sailing time constraint that designs the ship course and speed simultaneously to realize an optimal navigation strategy, including route and speed configurations for the whole voyage. This navigation strategy can ensure that a ship arrives at its destination close to its ETA. Such an optimized navigation program can not only avoid the navigational risks brought by a heavy sea state but can also generate an optimal route with minimized fuel consumption while arriving at the required ETA.

In this paper, a rasterized dynamic sea area model was constructed by processing meteorological information, and speed configuration was added to the route optimization; then, single nodes were extended to collections of “ext-nodes” for optimization. Based on the actual voyages’ records, an artificial network model was modeled to predict the fuel consumption rate and RPM to assist the route optimization. Then, the evaluation function and operation process of an A* algorithm were designed to achieve the simultaneous optimization of vessel course and speed based on the premise of arriving at ETA. The results of two optimization simulations of the winter route in the North Pacific Ocean showed that the algorithm proposed in this paper could not only guarantee navigational safety but could also reduce fuel consumption by 130.93 tons or about 7.25% compared with the historical route. Therefore, the methods in this paper can help ships reduce their fuel consumption while avoiding windy and stormy areas and still arrive at their ETA, which has particular practical significance in accomplishing complex transoceanic tasks.

In the future, we will consider more specific complexities of the sailing environment (such as currents, ice-covered waters, etc.) by obtaining and analyzing data from more voyages or different types of ships. In addition, while an ANN with a single hidden layer has great performance, there is space to optimize the prediction of the fuel consumption rate, which can be further investigated. We have up to now only used the three-dimensional A* algorithm to plan a low-fuel-consumption route under a time constraint, while the needs of ship navigation may be very diverse, and future attempts will be made to use such three-dimensional algorithms to solve multi-objective ship navigation tasks and further improve the computational speed.