Voyage optimization combining genetic algorithm and dynamic programming for fuel/emissions reduction

Deterministic optimization algorithms generate optimal routes/paths and speeds along ship voyages. However, a ship can rarely follow pre-defined speeds because dynamic sea environments lead to continuous speed variation. In this paper, a voyage optimization method is proposed to optimize ship engine power to reduce fuel and air emissions. It is a combination of dynamic programming and genetic algorithm to solve voyage planning in three-dimensions. In this method, the engine power is discretized into several levels. The potential benefit of using this algorithm is investigated by a medium-size chemical tanker. A ship ’ s actual sailing is used to demonstrate benefits of the proposed method. On average 3.4% of fuel-saving and emission reduction can be achieved than state-of-the-art deterministic methods. If compared with the actual full-scale measurements, on average 5.6% reduction of fuel consumption and GHG emissions (about 275 tons) can be expected by the proposed method for the six case study voyages.


Introduction
Maritime transport emitted more than 1000 Mt CO 2 every year accounting for 2.89% of global greenhouse gas (GHG) emissions.The GHG emissions from shipping could increase 90-130% of 2008 emissions by 2050 due to the growth of world maritime trade (IMO 2020).To help decarbonize shipping industry, IMO introduced several regulations for both ship design and ship operations, such as EEDI and SEEMP both entered into force in 2013.Furthermore, IMO MEPC72 (2018) set up an initial strategy to reduce GHG emissions from the shipping industry by at least 50% by 2050 compared to 2008.To achieve the goal of emission reduction from shipping, various energy efficiency measures should be developed and tested in the maritime industry.
A comprehensive review of various energy efficiency solutions ranging from ship design maintenance, operation and logistics, etc., as well as their potential CO 2 reduction is presented in Bouman et al. (2017).For example, Kontovas (2014) proposed a conceptual approach to solve ship routing and scheduling problem by identifying important parameters that can be used to accurately estimate a ship's air emissions and fuel consumption.As noted in Kontovas (2014), in most of maritime transport related literature, a ship's fuel consumption is assumed to be described to ship speed to the power between 3 and 4.But from ocean engineering perspectives, a physical model to describe a ship's fuel consumption and air emissions is much more complex and includes at least the modelling of a ship's mechanics system, propulsion system and thermodynamic system (Tillig et al., 2017).To avoid complex physical modelling, ship data-based regression and machine learning methods become more popular to estimate the fuel and emission problems, such as Mao et al. (2016aMao et al. ( , 2016b)), Merien-Paul et al. (2018) and Wang et al. (2018aWang et al. ( , 2018b)).It is well-known that a ship's speed can significantly influence the fuel consumption and air emissions along a voyage.Finding optimal sailing/service speed for ships of similar sizes can be extremely beneficial to reduce the fuel and air emissions from the entire maritime transport sector (Wang and Meng 2012, Wen et al., 2017, Venturini et al., 2017, etc.).
Furthermore, for individual ships and their navigation, a market study DNVGL (2015) to more than 80 shipping companies shows that one of the most promising measures is the computer aided voyage optimization that can help reduce fuel consumption, air emissions, and comply with IMO environmental regulations (SEEMP).A voyage optimization system is used to choose the "best" sailing course/route and corresponding speeds for each individual voyage based on forecasted sea environments (Simonsen et al. 2015, Wang et al. 2019, etc.).This paper focuses on the voyage optimizations with objectives of reducing fuel consumption (air emissions).Most voyage optimization systems use the speed as a control variable for the voyage planning, such as the conventional Isochrone methods (Hagiwara 1989, Lin et al. 2013), 2D dynamic programming-based methods (Dijkstra 1959, Chen 1978, Calvert, 1990, etc.), and 3D dynamic programming methods (Shao and Zhou, 2012, Zaccone et al., 2018, Wang et al. 2019 etc.).The 2D methods are mainly used for the course optimization, while the 3D method can allow for speed optimization along optimized ship routes but may lead to locally optimal results.Large research efforts have been devoted to developing evolutionary algorithms for multi-objective global optimization (Hinnenthal 2008, Maki, 2011, Vettor and Guedes Soares 2016, Lee et al. 2018, etc.).
However, these genetic algorithms may be not efficient enough to solve a long voyage planning problem.Furthermore, the outputs (control variables) of almost all the above voyage optimization algorithms are ship routes/paths and a series of speeds along the paths.Due to dynamic sea environments encountered during actual sailing, a ship's speed can vary significantly even under one same stationary sea state with the same engine power/RPM setting (Mao et al., 2016a(Mao et al., , 2016b)).Furthermore, it is not energy efficiency and impossible for seafarers to continuously adjust engine settings to follow a series of given speeds, because more than 90% of a ship's navigation at open sea is on autopilot mode that can actively adjust ship course to follow planned path (ABB 2020) while engine settings are manually adjusted normally once a day to change ship speeds.The difference between a ship's planned and actual sailing speeds from voyage optimizations may lead to serious discrepancy in arrival punctuality (Skoglund, 2016).Moreover, a ship is often navigated by setting her engine power or RPM, but voyage optimization using engine settings as decision variables is very seldom researched because of many iterations to solve the complex relationship between engine power and ship speed under dynamic wave environments.In this study, a combined optimization algorithm is proposed to consider engine power as a decision variable for voyage planning to lower fuel cost (air emissions) and increase arrival punctuality.The proposed method employs the deterministic optimization algorithm into a genetic algorithm framework for fast convergence and reliable optimal results.In the following, the basic concept and mathematical formulation of the proposed optimization method are presented in Section 2. The scheme and major components of the methods are presented in Section 3. Its implementation is introduced in Section 4. The effectiveness of the method for fuel saving and GHG emissions reduction is investigated in Section 5 and 6 through comparing with full-scale measurement of 6 typical case study voyages.The paper is concluded in Section 6.

Concept of the proposed voyage optimization method
In this study, a voyage optimization method is proposed to help plan the "best" route/path and engine settings along the path to achieve minimum fuel consumption (GHG emissions) and arrival punctuality, etc.The proposed method is composed of four elements as in Fig. 1, i.e., basic voyage information, ship models (cost functions), metocean data, and voyage optimization algorithms.A ship's voyage optimization problem can be mathematically formulated by the three components (Ng 2019), i.e., 1) decision variables, 2) Fig. 1.Overview of the proposed voyage optimization method.
H. Wang et al. objective function and 3) constraints.

The decision variables/vectors
First, voyage information is used to discretize potential sailing area into an optimization graph as in Fig. 2.Then, optimization algorithms can formulate, evaluate, and select candidate voyage planning solutions based on the estimation of their costs.In the proposed optimization method, the candidate voyages are referred to the following decision variables/vectors, N i,ji = (x i,ji ,y i,ji ): waypoints with longitude x i,ji and latitude y i,ji at the j i -th state of the i-th stage in the discretized graph system (as in Fig. 2) representing the sailing/searching area of a voyage.Note that the waypoint N i,ji and all the waypoints in its preceding stage N i− 1,ji− 1 are connected by rhumb lines to form sub-paths/edges E i,ji (j i− 1 ), j i− 1 = 1, 2, ⋯, J i− 1 , where J i− 1 is the number of discretized waypoints in the (i-1)-th stage.Let assume that the sailing area is discretized into M stages in the graph system.

The objective function
In a voyage optimization system, the most optimal voyage is selected based on pre-defined objectives.For a specific voyage, let T ETA denote the sailing time corresponding to a ship's expected time of arrival (ETA), and T denote the sailing times of all candidate voyages (represented by decision vectors) in the optimization process.Let f fuel ) , i = 2, 3, ⋯, M denote the cost functions to estimate a ship's fuel consumption and travel time along the sub-path/edge E i,ji (j i− 1 ).Models to estimate the cost functions are presented in Section 2.4.The task of a voyage optimization algorithm is to find a series of waypoints, sub-paths E i,ji (j i− 1 ) and corresponding power settings P(E i,ji (j i− 1 )) to fulfil the following objectives, i.e., to minimize fuel consumption (GHG emissions) and increase arrival punctuality in this study, ) ) Then the sub-paths at different stages E i,j i (j i− 1 ) are connected to form a complete ship path.However, the ETA and fuel consumption (air emissions) of an optimal planning (decision vector) always conflict with each other.It means that a decision vector with a shorter ETA often leads to higher fuel consumption (air emissions), and vice versa.In the proposed method, the Pareto front with multiobjective outputs is used in the iteration process to generate optimal decision vectors as in Section 4. Fig. 2.An example of the grid/graph system and its labeling in the optimization method.
H. Wang et al.

Constraints for the voyage optimization
In the proposed optimization method as in Fig. 1, various optimization algorithms need to search huge amount of decision vectors (combination of sub-paths and engine settings).But too many decision vectors as inputs to the genetic algorithm can make the estimation of voyage optimization extremely time-consuming.Furthermore, concerning a ship's maneuverability and actual navigation limitations, the following constraints are considered in the optimization method: where G represents the discretized waypoint grid/graph system for the admissible geographical sailing area of a voyage, where the land-crossing constraints, obstacles avoidance along sub-paths and traffic separation zones should be considered.While α notes course changes at the waypoint N i,ji between two adjacent sub-paths, i.e., α(N i,ji ) :=< (E i,ji (j i− 1 ), E i+1,ji+1 (j i )) , and P = [p 1 , p 2 , ⋯ , p C ] represents the pre-defined power levels from a ship's engine configuration.In addition, the ECA (emission control area) is considered in a simple way as a constraint in this method.If a ship leaves the ECA in an ocean crossing voyage, she will not enter the same ECA in the same voyage because too much effort is needed to change fuels for sailing in the ECA.There are also some safety related constraints in an actual weather routing system, such as maximum allowed motions at sea, which are not considered in this study.

Models and metocean inputs for fuel consumption (GHG emissions)
A ship's cost functions to describe her speed-power relationship is essential for any voyage optimization methods.For the cost analysis, a ship's voyage is often assumed to be composed of a series of stationary sea states lasting for from 30 min to several hours.In this study, it is chosen as the same resolution as metocean data inputs.Each stationary sea state is represented by the significant wave height H s , wave period T p , and wind speed U w , etc.Under a stationary sea state, a ship's energy system is assumed to be quasi-static.Therefore, one-to-one relationship between engine power and (mean) speed under a stationary sea state should be constructed in terms of various sea parameters.For the power-based voyage optimization proposed in this study, engine power is the main input variable.A B-spline interpolation method is used to inversely estimate ship speed from a given power setting.To perform the inverse analysis, the speed-power relationship is modelled as follows for voyage optimization.Firstly, under a stationary sea state, the required power P b to push a ship forward at (mean) speed V to overcome the total resistance R TOTAL is evaluated as: where η prop , η s , η h , η o ,η r represent a ship's propulsive efficiency, the shaft efficiency, the hull efficiency, open water efficiency, and relative rotate efficiency, respectively.Under such a stationary sea state, the total resistance R TOTAL is estimated by, where R CALM ,R AA and R AW denote the calm water resistance, wind resistance, and added resistance due to waves, respectively.A chemical tanker is used as case study ship to investigate the effectiveness of the proposed method.For this ship, the calm water resistance and all the propulsive efficiency are available from the towing tank tests, and the wind-induced added resistance is calculated as ISO (2015).The estimation of added resistance due to waves R AW is one of the important tasks for the cost analysis and voyage optimization.Under a stationary sea state, R AW can be estimated by, where S(ω|H s , T P ) represent the wave spectrum of the sea state.The wave spectrum is modeled as a standard JONSWAP spectrum (Lewis 1988).The transfer function of added wave resistance R aw (ω) depends on different ship parameters and loading conditions.In this study, the semi-empirical model for R aw (ω) developed by Lang and Mao (2020) is used.It should be noted that this model has been validated by the ship's full-scale measurements.
Finally, the fuel consumption rate r fuel (kg/hour) and GHG emission rate r CO2 (kg/hour) for the required engine power P b due to R TOTAL is estimated by: where SFOC represents the specific fuel oil consumption, which is obtained from the ship's energy performance measurement.Since the case study voyages in the following analysis are measured in the North Atlantic and only the heavy fuel was used, it is reasonable to assume that various air emissions, such as GHG emissions, NO x , SO x , PMs, are proportional to the ship's fuel consumption.Therefore, H. Wang et al. the cost of fuel consumption is used to represent both fuel cost and air emissions.In addition, only GHG emissions are chosen as a representative to demonstrate the benefits of the proposed method to reduce air emissions.The emissions factor is set to be 3.1144 kg CO 2 /kg fuel in this study (Kontovas 2014).Finally, for a sub-path E i,ji with engine power P λ , the corresponding sailing time and fuel consumption should be estimated iteratively by the following simplified formula, where ‖E i,ji (j i− 1 )‖ represents the sailing distance along the sub-path E i,ji (j i− 1 ).Every sub-path E i,ji (j i− 1 ) may be associated with various wind and wave conditions due to different arrival time.In the voyage planning process, one of the most necessary inputs is the surrounding metocean conditions, such as wind, wave, and current, etc., which are used to analyze a ship's safety, fuel consumption and sailing time as above.For actual voyage optimization, metocean forecast is updated every 24 h.In this study, the metocean parameters of wind (wind speed and direction) and waves (wave height H s and period T p ) are extracted from ERA5 dataset (ECMWF 2019), and the parameters of current (speed and direction) from the Copernicus server (Copernicus 2019).The temporal resolution of the data is 3 h, and the spatial resolution of the data is 0.5 • × 0.5 • , i.e., about 55 km × 55 km.An example of various metocean parameters is shown in Fig. 3.

Scheme and major components of the proposed method
In this study, the proposed voyage optimization method searches for optimal sailing paths associated with the best engine power configurations by an evolutionary algorithm, which has the potential to provide global and multi-objective voyage optimization solutions.However, the shortcomings of using evolutionary algorithms in the voyage optimization problem are also quite noticeable.For example, the initialization inputs of those algorithms can significantly influence the convergence and robustness of voyage optimizations.The randomly generated samples (refer to the decision vectors in this study) evolved in the voyage optimization process usually contain many variables.It may cause optimization convergent to a locally optimal solution or even without a proper solution.
To overcome the shortcomings, in this paper several techniques from deterministic voyage optimization methods are integrated into the genetic algorithm, which is supposed to generate a ship's optimal sailing path associated with proper engine power settings.The workflow of the proposed optimization method is shown in Fig. 4. The main input parameters and control variables are the so-called decision vectors, which contain a potential sailing path and its corresponding engine power settings.Here, three different decision vector generators are included: 1) Heuristic decision vector generator (HDVG), 2) Deterministic decision vector generator (DDVG), 3) Stochastic decision vector generator (SDVG).
In the optimization process as in Fig. 4, HDVG is only used in the initial population generation (Step I as the right part of Fig. 4), while DDVG and SDVG are used in both Steps I and IV to generate new populations (decision vectors).The following subsections will present the main ideas of the three population generators, and the decision vectors as the population samples for the voyage optimization.

Decision vector
In a ship's voyage optimization process, the potential sailing area between the departure and destination should be first discretized into a grid/graph system that composes a series of waypoints/nodes.In the grid/graph system, sub-paths/edges connected between adjacent waypoints are usually directed toward the destination.Those sub-paths/edges are assigned by both attributes and weights.The attributes refer to the ship's operational control variables, i.e., ship speed, engine power/RPM.The weights represent the costs of the sub-paths/edges related to different objectives, i.e., fuel consumption, air emissions, sailing time.The target of the voyage optimization algorithms is to find a combination of sub-paths/edges forming a feasible route from the departure to destination with optimal attributes for each sub-path.Fig. 3.One example of the metocean data in the North Atlantic.
H. Wang et al.Firstly, the potential sailing area is discretized into a series of time stages (say M stages).The waypoints/nodes N in all the time stages and sub-paths/edges E connecting waypoints at adjacent time stages can form a grid/graph system G = (N, E).An example of the grid/graph system deployed in this study is shown in Fig. 2. The upper plot demonstrates the graph system and the bottom plot shows the corresponding labeling system.For example, the waypoint/node at the j i -th state in the i-th stage of the graph system is denoted by N i,ji with longitude x i,j and latitude y i,j .The sub-path/edges E i,ji associated with the waypoint N i,ji is formed by connecting it with all waypoints in its preceding stage, where the symbol "→" means connecting two waypoints to form a sub-path, and J i-1 is the number of waypoints in the preceding (i-1)th stage.A complete route/path-vector w from the departure to the destination is formulated by selecting one adjacent sub-path at each individual stage as, Secondly, a ship's engine power is configured such that only C discrete power levels P can be set, For one complete route vector w, the number of possible path-vectors is 1 By assigning different power levels to each sub-path of w can generate C M− 1 possible voyage planning candidates.Then, the total number of the candidates becomes huge for this 3D voyage optimization problem.Therefore, it is necessary to eliminate some unrealistic, bad performed sub-paths by introducing sailing constraints and optimization methods such as dynamic programming or Dijkstra algorithms.Finally, the voyage optimization is to choose the two best vector variables, i.e., the path vector w * and the power vector P f (w * ) along the path vector w * .In the proposed method, one sample of the route population is formulated by the two vector variables and is named as a decision vector D v , where the subscript v represents the sample index within the population of decision vectors, P(λ i,ji,ji− 1 ) is the power setting applied for the sub-pathE i,ji (j i− 1 ).In the following, to ease the description, the above decision vector D v is simplified as, where E i = 2, 3,…, M represents all sub-paths at the i-th stage, e i is the sub-path index among all sub-paths at the i-th stage, and λ i is the power index λ i ∈ [1, 2, ⋯, C] such that P(λ i ) ∈ P as in Eq.( 10).

Heuristic decision vector generator
To increase the convergence efficiency of the genetic algorithm, the proposed method introduces a heuristic decision vector generator to construct the initial population.For deterministic approaches such as dynamic programming, it is difficult to optimize engine power settings along a sailing route (sub-path vectors).Furthermore, adding extra dimensions to the Euclidean space will exponentially increase the number of searching candidates in the deterministic approaches.It could lead to local optimal results or may even cause the optimization process crashed.However, those deterministic approaches can be implemented to find the optimal path vector for a fixed engine power (speed) as a conventional 2-dimensional optimization problem.In this study, Dijkstra's algorithm is employed to find such an optimal path vector w * h,c with minimum fuel consumption for each individual power level inP = [p 1 , p 2 , ⋯, p C ].A heuristic decision vector D * h,c for the fixed power P(λ) = p λ is denoted by, where λ = 1, 2, …, C corresponds to various power levels, and sub-path index e h* i (i = 2, 3, …, M) refers to the optimized path vector that can be identified by a dynamic programming method.Then, optimal routes for all power settings within P form an optimal route vector W * H as, Combining all these optimal routes/paths can form a blue region shown in Fig. 3, which further confines the search space from the original graph system G.The heuristic decision vectors can decrease the search space W and eliminate the possibility of generating decision vectors with low qualities.Such an example is presented in Fig. 5.The confined search space G c is generated between the maximum and minimum label numbers of the nodes in each stage through the heuristics.It shows that the size of the search space has been greatly reduced by the heuristic method.

Deterministic decision vector generator
The deterministic decision vector generator (DDVG) generates optimized decision vectors D * h,λ by a dynamic programming method.In this study, Dijkstra's method like that proposed in Wang et al. (2019) is used for the DDVG, which optimizes the power settings along a fixed path.Firstly, the path vectors of optimized routings (decision vectors) from the last iteration of the genetic algorithm are used as input s to the DDVG to generate new population samples.For example, for the initial population generation (Step 1 in Fig. 4), the path vectors w * h,λ of D * h,λ , λ = 1, 2,…, C generated from the Heuristic method are used as inputs for the DDVG method to generate new population samples.While for the evolution of genetic algorithm (Step 4 in Fig. 4), the path vectors of optimized decision vectors selected from the preceding generation are used for the DDVG method.Secondly, based on the inherited path vectors, ship engine power settings along with each path vector W* are optimized with minimum fuel consumption by a dynamic programming approach, which begins with an input path vector w * i from a decision vector.In the population initialization step, the decision vectors are coming from D * h,λ (λ = 1, 2, …, C).Let take a typical path vector for example as in Fig. 6 to demonstrate the DDVG method.For the chosen path vector, i.e., the blue path in Fig. 6, the DDVG will optimize a ship's power settings P f (w * ) along the path vector w * .
Let simplify the notation of the node at the i-th stage by N i , which is assigned with various power settings P for sailing to the next stage.Since different power settings lead to different arrival times and fuel consumption (air emissions), a ship's arrival time at each stage is divided into K intervals, corresponding to K different fuel consumption (GHG emissions).The arrival times within the same time interval are assumed to have the same arrival time (ship status), which can be denoted by It contains various arrival times T i,k and different fuel consumptions F i,k when reaching at N i .The variation of the status S i is caused by different engine power settings during the optimization process.The initial status at the departure point, i.e., the first stage, is defined as S 1 = (0, 0).The procedure to implement a dynamic programming optimization method is described in the following steps: 1. Assign every the engine power levels p λ (p λ ∈ P) from the departure to the second stage along the selected "optimized" sub-path . All power levels are kept fixed along the sub-path E * 2 When a ship advances to the third stage with C 2 possible arrival times (C 2 > K), the arrival times will be re-discretized into K new time intervals as in Eq. ( 15).For each time interval, only keep the power setting leading to the lowest fuel consumption.
2. Let the ship advances to the i-th stage at the node N i .By assigning all powers in P to N i , the ship's status of the next stage, i.e., the i + 1 stage, can be calculated as, Fig. 6.An example of an input path vector for DDVG with simplified notations of nodes and status.
H. Wang et al.
where ΔT(S i,2 , p λ ) is the sailing time from the status of S i,2 at the node N i sailing to N i+1 using power setting p λ , and ΔF(S i,2 , p λ ) is the corresponding fuel consumption for such a sailing.
3. There are C × K elements of arrival times in the i + 1-th stage.These arrival times are again divided into K time intervals/sectors.In each sector, select the status with the minimum fuel consumption (air emissions), save it as one status in S i+1 and reserve corresponding power settings into P * i+1 , where the indexes λ * ,j are obtained by the dynamic programming optimization algorithm.4. Let i = i + 1, repeat the above process from Step 2 until i = M. Finally, the dynamic programming method generate J groups (series) of engine power settings for the chosen "optimized" path vector w * .Combining the path vector with every series of engine power setting from the above steps, J more "optimized" decision vectors are obtained.

Stochastic decision vector generator
In order to increase the coverage of the path vector samples rather than the C fixed path vectors W * H from the Heuristic method, a stochastic decision vector generator (SDVG) is proposed to generate more diverse/fresh path vectors and decision vectors in both the population initiation and evolution processes of the genetic algorithm.
In the SDVG method, one sub-path at every stage over the confined graph G c is randomly selected to generate a complete route (a path vector w).The random sub-path selection is based on the algorithm proposed by Yen (1971), which finds M shortest sub-paths from the confined graph G c .Then, each generated path vector w will be assigned with a series of power settings along the path vector.The engine power setting vector is generated by a random walk approach.In this study, the step length in the random walk process follows a discrete uniform or Gaussian distribution.The limit of a maximum of 2 step lengths is set to control the change of power switch in a reasonable range.In the random walk process, if the generated power setting exceeds the boundary of the power setting P, the power setting at this step will be assigned as the boundary value of P. Finally, by combining the randomly generated power series with path vectors, various decision vectors D v are generated in this stochastic way for the iteration of the voyage optimization.

Implementation procedures of the proposed optimization method
The genetic algorithm (Goldberg 1989) is inspired by Darwin's theory of evolution (Darwin 1859).It uses an evolutionary process to solve problems.The evolutionary process aims at obtaining better generations compared to previous ones.In the proposed method, the number of individuals/decision vectors in the population has a great influence on the optimization results.Pareto front points (with optimal fuel cost and ETA) are considered as the qualified results and selected as the individuals for the next generation during the optimization iterations (Step 4 in Fig. 4).To ensure the searching quantity in each iteration, the Dijkstra's algorithm to get global optimal route for fixed power settings, and the dynamic programming method for fixed route/path, are introduced to initially confine potential decision vectors (composed of a series of waypoints) for the genetic algorithm.Furthermore, the number of individuals/ decision vectors should set as large as possible.In this study, 500 decision vectors are set in the population.
H. Wang et al.

Initial population generation
The genetic algorithm starts with generating a set of decision vectors (initial population).To improve the quality of the initial population, three decision vector generators are used as presented in Fig. 4: 1. Execute HDVG, which generates C decision vectors.Each decision vector is corresponding to one optimized ship routing for one engine power input within the set P in Eq. ( 10).In addition, the HDVG will also generate a confined graph (search space) as in Fig. 5.The routes/path vectors from the HDVG method is used as the basis for the DDVG and the confine graph is used as the basis for the SDVG to generate more samples of decision vectors for the initial population in Step 1. 2. Execute DDVG to generates 20% of the initial population.Based on the C optimal routes (only the path vectors are used) generated from the HDVG, this DDVG method assigns all discretized powers within P to all the sub-paths along the C optimal routes.For one optimal path, the dynamic programming method proposed in Wang et al. ( 2019) is used to generate 20 optimal routings of minimum fuel cost for various ETAs.The optimal routings generated by the DDVG based on the C optimal routes as inputs are then putting together to form a Pareto front with respect to fuel consumption and ETA.Finally, about 100 optimal routes (corresponding to 20% of the initial population) along the Pareto front with their ETAs closed to the required ETA are selected as the initial population in Step 1. 3. The rest of the population is generated by SDVG, to increase the diversity of the initial population and prevents early convergence.
In this method, sub-paths from the confined graph and power levels are picked randomly for each stage to form many decision vectors.Then 80% of the initial population of decision vectors are selected from the randomly generated decision vectors with their estimated ETA close to the required ETA.

Cost evaluation and Pareto front of decision vectors
After obtaining the next generation of population as in Fig. 4, the ETAs, and fuel consumption (GHG emissions) of all the decision vectors within the population are evaluated to form Pareto front (minimum fuel consumption for various ETAs), which is prepared for the evolution of the next generation in the proposed genetic algorithm.The ETA and fuel consumption (GHG emissions) of each decision vector in the population are estimated by the methods presented in Section 2.4.For each decision vector, its ETA is plotted against fuel consumption as one point in a two-dimensional design space as in Fig. 7, where a contradiction trend between the ETA and fuel consumption can be observed.
As illustrated in Fig. 7, for each ETA corresponding to many decision vectors, one decision vector with the lowest fuel consumption (GHG emissions) could be selected as an optimal route solution.The decision vectors with the lowest fuel consumption and emissions in terms of a series of ETAs will form a boundary in the two-dimensional design space.Such a boundary is called Pareto front and the set of solutions is called Pareto optimal solutions.The evaluation and selection of the Pareto solutions (decision vectors) in each iteration are based on the following criterion.A sample decision vector with ETA and fuel consumption, denoted by [ T pareto , F pareto ] , can be chosen as one of the pareto solutions if and only if the ETA of every other sample decision vector, whose fuel consumption is less than or equal to F pareto , is greater than T pareto .Normally, for multi-objective voyage optimizations with varying ETAs, all the generated Pareto front solutions for the current evolution should be used as parent decision vectors in the population of next evolution.
However, in the following analysis, case study voyages with full-scale measurements are used to demonstrate the effectiveness of the proposed method.The objective related to a ship's travel times in a voyage is re-defined such that the ETA should be close to the ship's measured travel time.Therefore, in this study, only 80% of the optimal solutions in the Pareto front, whose ETAs are the closest to the measured ETA, are chosen and remained.The selected decision vectors are illustrated as the blue dots in the magnified plot of Fig. 7.In Step IV, these selected decision vectors are used for the next population generation by the DDVG and SDVG methods.

Population evolution (transfusion and crossover) and results
In Step IV, the population evolution is coming from three sources.The first source is the 80% Pareto optimal solutions inherited from the last iteration.The second source is obtained from the transfusion of new decision vectors by the DDVG and SDVG methods.The two methods work similarly as in the initial population generation of Step 1.The only difference is that in Step 1, the optimized routes from the HDVG method are used as inputs for the DDVG and SDVG, while in Step 4, it is the 80% Pareto optimal decision vectors from last iteration are used for the inputs.Finally, the third source is from the crossover operator based on the same Pareto optimal solutions as the DDVG and SDVG methods.Different from conventional crossover methods, the proposed method implements the crossover for the sailing path and corresponding power settings, separately.The flowchart of the crossover operator is shown as Fig. 8.
The Pareto optimal decision vectors are referred to as parents in the crossover operator.Firstly, two parents are randomly selected from the remained Pareto optimal solutions and named "parent 1" and "parent 2".Each parent is composed of path and power configurations.Then, a new path vector is generated by the following steps: 11.Randomly select Node A from path 1 and Node B from path 2; 21.Check if there is a feasible sub-path to connect Node A and Node B; 31.If the connection sub-path exists, the new path is formulated by combining the sub-path from the departure to Node A in the "path 1", the connection sub-path and the sub-path from Node B to the destination in the "path 2"; 41.If feasible sub-path connection A and B does not exist, go to the above step P1 to repeat the process.
Next, the above process is used to generate new power configuration to be assigned in the new path to form a new decision vector in the new population.
Since parts of the results from the deterministic method and the heuristic method are highly optimized results, the evolution of the population is expected to converge within a small number of iterations.In this study, no explicit target (refer to fuel consumption) is available for the optimization and it is difficult to evaluate the quality of evolution.Thus, a maximum/stop number of iterations is used as the stop criterion and it is set to be 5 iterations.

Case studies and overall results from the power-based voyage optimization
One way to verify if the proposed method can give optimal solutions is to compare with full enumeration for all possible solutions, Fig. 8. Flowchart of the crossover operator.
H. Wang et al. but current computability does not allow us for the full enumeration.In this study, full-scale measurements from a medium-size chemical tanker are used to investigate the capability of the proposed power-based optimization method for fuel-saving and emission reduction.Then, the proposed method is compared with the actual sailing (also aided by a market weather routing system) and the heuristic planning results.The effectiveness and benefits of this heuristic method is validated by two recent studies (Wang et al. 2019;2020).

Measured voyages involved in the case study
The case study ship is navigated by setting engine power configured with 8 ECO power levels, i.e., 8490KW, 7880KW, 7270KW, 6660KW, etc.The ship's main characteristics are listed in Table 1.A conventional voyage optimization system was installed for her voyage planning.Therefore, the measured routes are combinations of the ship captain's experience and advice from the onboard voyage planning system.The measured performances along the actual sailing routes are used to investigate potential benefits brought by the proposed method.Since the North Atlantic is known as one of the most challenging areas for ship operation, a good voyage optimization system would be greatly beneficial to safety and energy efficiency for ships sailing in this area.The full-scale measurement of the studied ship when sailing in the North Atlantic from the year 2015 to 2016 is used to verify the effectiveness of the proposed optimization method for voyage planning.Six case study voyages are selected as in Fig. 9, while three of them are eastbound voyages and the others are westbound voyages.The red triangle-up and the blue triangle-down represent the departure and destination of the voyages, respectively.The ship's full-scale measurement data contains various dynamic sailing and performance information, e. g., waypoints and their arrival times, heading, shaft power, ship speeds, etc.These measured parameters are considered as a reference to be compared with the optimized results by the proposed method.

Overall results from different route planning methods
In order to study the capabilities of the proposed method for voyage optimization in terms of minimum fuel consumption and GHG emissions with pre-defined sailing time (ETA), two conventional route planning methods are taken for the verification and comparison purpose, i.e., the actual sailing routes from the full-scale measurements, and the fixed power 2D voyage optimization method.The actual sailing routes were planned by a weather routing system from the shipping market.The fixed power optimization method uses Dijkstra's algorithm.It is implemented based on our recent development from Wang et al. (2019), and the optimized voyage solution is one of the optimized routes generated by the Heuristic method as in Fig. 4. In this study, the Heuristic method will generate 8 optimum routes, corresponding to 8 different ECO power levels of the case study chemical tanker.The optimum route with sailing time (ETA) closest to the actual route is chosen for the comparison.
It should be noted that to allow for a fair comparison, fuel consumption and GHG emissions in all the three voyage planning methods are calculated by the cost models presented in Section 2.4.The actual routes were optimized using the forecast sea environments, while both methods that generate "Heuristic route" and "Optimal route" are using the hindcast data.Since during the actual sailing, the onboard conventional voyage optimization system updated their forecast data every 24 h, and weather model can predict well of the sea environments 3-5 days ahead in particular in open sea areas, it is reasonable to assume that the difference in weather inputs for optimization between the actual sailing and the proposed voyage optimization will not influence the comparison results.
In this study, only parts of voyages measured in the open sea (longitude between − 70 • W to − 10 • W) are selected as case study voyages for the optimization and comparison.Table 2 lists the results of ETA, fuel consumption, GHG CO 2 emissions and sailing distances from the three voyage planning methods.It shows that the proposed method has a great potential of saving fuel and GHG emissions up to approximately 13% compared with the actual sailing route for the voyage 20160317.In general, the fuel consumption along the optimized routes by the proposed method is on average 5.6% and 3.4% less than those from the actual sailings and the heuristic routes, respectively.For these voyages, it is corresponding to 275 tons GHG emission reduction.
Normally, more improvement of fuel-saving and GHG emission reduction should be expected by a voyage optimization method for the westbound voyages than the eastbound voyages, since the storms are usually moving towards the west in the North Atlantic with harsh encountered sea environments.While the large fuel-saving and emission reduction for the Eastbound voyage 20,160,307 is because encountering a very large storm by the actual sailing.It will be further investigated in the following.For the westbound voyages, the proposed method reduces a bit less percentages of fuel and air emissions than the eastbound voyages.This might be contributed from careful voyage planning by captains and the onboard weather routing systems.However, the proposed method can still help to reduce 4.6% fuel consumption than the actual routes.In addition, the sailing distance along the optimized routes are almost the same as the proposed method and the Heuristic routes, but the proposed method can still help to further save more than 4% by adjusting power setting along the Heuristic routes.Moreover, the reduction of the total amount of fuel and emissions in tons along westbound voyages is more than the eastbound voyages.

Optimized routes/paths from different methods
The optimized ship routes/paths by the three voyage planning methods are presented in Fig. 10.The voyages are categorized into eastbound voyages and westbound voyages.Most of the measured (actual sailing) routes went through either great circle routes or rhumb lines.It may indicate that the major part of these routes is associated with relatively calm sea environments.In this case, the shortest paths (the great circle routes) or the easiest paths for operation (the rhumb lines) are preferable by the shipmasters.However, the eastbound voyage 20,160,317 traveled differently from other voyages since the ship may encounter storms during the conventional routes.In this case, great involuntary speed reduction could happen if sailed through the great circle route.In Fig. 10, some sub-plots have not shown the red lines clearly, because the red lines coincide with the blue lines.It implies that the routes provided by the heuristic method are of high quality.The proposed method further optimized engine power configurations for the heuristic routes.

Efficiency of the proposed method
Both deterministic and stochastic methods are involved in the genetic algorithm evolution.In Fig. 11, the red vertical line represents the target ETA (measured total sailing time), while the red diamonds represent the fuel consumption and ETA of the actual routes.In each iteration of the evolution, DDVG generates new deterministic decision vectors and SDVG generates new random decision vectors.It is hard to define a fitness function that can evaluate the quality of the population evolution by the DDVG and SDVG.Fig. 11 demonstrates that the heuristic method can provide highly optimized results marked as black crosses.This might be the reason why the two-dimensional voyage optimization algorithms are still widely used in today's weather routing systems.Those results are considered as good resources for guiding the evolution of the proposed method towards fast convergence.This can also allow voluntary speed reduction during the voyage planning process.During the evolution iteration in the proposed method, the population evolves towards the target ETA.Fig. 11 shows that the results of the final generation are distributed around the target ETAs.
On the other hand, the voyage optimization problem defined in this study includes the optimization of two decision variable sets, i.  e., the sailing path-vector and corresponding power settings along the path-vector, to make a balance between the fuel consumption (GHG emissions) and the ETAs.Each of the decision variable set contains many variables, such as a path vector composed of many waypoints.Stochastic optimization algorithms such as the genetic algorithm for solving the problem with many variables could hardly obtain the same results every time when conducting the same optimization.However, by increasing the number of iterations, there will be more optimal or Pareto optimal points in the design space and more chances for obtaining identical solutions.Furthermore, the efficiency of a stochastic optimization-based method is often hard to measure due to uncertainty in searching convergency.In the proposed method, the genetic algorithm is combined with the deterministic dynamic programming methods.They can help to significantly reduce the searching space and thus reduce the runs of iterations to reach optimal results.In addition, a lot of information is stored in the computer memory (database) after running the deterministic optimization, such as the grid system and the decision vectors from HDVG and DDVG.For the evolution process in the genetic algorithm, only the corresponding passing times, metocean conditions and costs needed be updated.Furthermore, parallel computing is implemented for generating decision vectors to  H. Wang et al. reduce computational time.For the voyage optimization of 10 days sailing period with waypoint resolution of approximately 2 • , running 7 iterations of the proposed approach requires approximately 3 mins, using a laptop with Intel® Core™ i7-8550U processor.

Further analysis of optimizing three typical voyages
In order to further understand the details of various voyage optimization methods contributing to different fuel consumption and air emissions, three typical voyages, i.e., Voyage 20160317, Voyage 20160607, and Voyage20160806, are selected for the detailed analysis as follows.The Voyage 20,160,317 is chosen because of quite harsh sea environments encountered during this sailing period.The Voyage 20,160,607 is a very normal voyage with typical encountered sea conditions, while the Voyage 20,160,806 is a rather very (extremely) calm sailing voyage with the highest significant wave height of less than 2.5 m.The interim results from the voyage optimization by the proposed method are presented and compared with the actual sailing routes.

Voyage optimization of harsh sailing environment
For the eastbound voyage 20160317, the evolution of the optimized ship route is plotted against the ship's actual sailing in Fig. 12, which also presents the surrounding sea conditions at about one-quarter and half of the voyage.Along the ship's actual sailing route, two storms were encountered by the ship, while the optimal route by the proposed method has successfully helped the ship avoid those storms.Especially, in the middle of the voyage, the actual route was encountered very harsh sea environments with significant wave height up to 10 m.In order to take a deep look at how the proposed method searches the optimal route to avoid storms, the planned power settings along the optimized routes by the proposed method and the heuristic method are presented in the upper plot of Fig. 13, as well as the series of engine power along the actual sailing route.It should be noted that the used engine powers are calculated by the speed-power model presented in Section 2.4 using the measurement data along the sailing route, such as the sailing speed, encountered wave, wind, currents, and water depth, etc.Before the power calculation, the model is also validated with the measured power performance.To have a fair comparison, the engine powers optimized by the other two methods are also estimated by the same model.The bottom plot of Fig. 13 presents significant wave height encountered along routes from the three voyage planning methods.
As is shown that the marine engine power reached its maximum point when hitting the storms, which result in a large increase in fuel consumption and a dramatic drop in engine efficiency.On the other hand, if the ship sails along both the heuristic route and the route optimized by the proposed method, it would encounter much less severe weather conditions than the actual route.When foreseeing storms occur during the voyage, it is judicious to slow down the ship, let storms pass and speed up to catch the expected time of arrival, because the ship's sailing speed and the energy efficiency will dramatically decrease when sailing in harsh weather condition.The upper plot also shows that the shaft power settings optimized by the proposed method follow such a sailing strategy to increase energy efficiency.During the calm weather conditions, high shaft power is used helping to increase the ship's sailing speed and low traveling time.The saved time by speeding up in calm weather conditions can provide enough space to reduce engine power (sailing speed) when the ship sails in harsh weather conditions.In this way, it can significantly reduce a ship's total fuel/power consumption along the entire voyage, i.e., more than 10% for this Voyage 20060317.

Voyage optimization for normal sailing environment
For the Voyage 20160607, the ship encountered the harshest sea environment of 5 m significant wave height in the middle of the voyage.It is rather a normal sailing and encountered sea environments in the North Atlantic.The evolution of voyage planning by the actual sailing and the proposed optimization method is plotted against her surrounding sea environments (represented by significant wave height H s of a sea state) in Fig. 14.It shows that the optimal route starts to deviate from the actual route at the departure point, and the ship would encounter slightly calmer sea conditions if sailing along the optimal route by the proposed method.In the middle of the voyage, the route optimized by the proposed method advances faster than the actual route attempting to move away from the latitudinal-moving storm.Furthermore, the time series of shaft power setting along the routes optimized by the three voyage planning methods are presented in Fig. 15 (upper plot).The corresponding significant wave height H s of sea states encountered along the routes are plotted in the bottom plot of Fig. 15.The shaft power setting along the optimal route by the proposed method is recommended to be 6000 kW at the lowest level when the ship is encountering storms or locally highest H s of sea states, such as the beginning of the voyage and when sailing at the areas with longitude around − 25 • , − 35 • and − 45 • .After passing these areas of global/local highest H s , the engine power is set to be 7270 kW by the proposed method and speed up the ship to catch up with the design ETA.While the actual sailing is still following very typical ship navigation, i.e., speed up at the beginning of a voyage (especially with moderate sea environment) to allow some navigation margin at the end of the voyage in case of unexpected large storms.It shows that the ship's engine power is set quite high even encountering quite high waves (H s ).The lower plot of Fig. 18 shows that even though the optimal voyage from the proposed method encountered only a little calmer sea conditions than the other approaches, a proper shaft power setting/optimization by the proposed method can reduce around 5.2% fuel consumption and GHG emissions than the heuristic route.

Voyage optimization for calm sea sailing
The encountered sea environments (represented by H s ) along the Voyage 20,160,806 are rather calm, less than 2.5 m.As shown in Fig. 10 for this voyage, routes/paths from the three voyage planning methods are quite similar.Furthermore, results in Table 2 also  H. Wang et al. show that little improvement in fuel saving and emission reduction can be achieved by the proposed method than the actual sailing route.The small and invariant confined graph/grid along the heuristic route constructed in the proposed method as in Fig. 16 may also imply that the weather conditions are quite calm during the voyage (no need to diverge from the heuristic route to avoid storms).Thus, the proposed method focuses more on the shaft power optimization for engine efficiency to reduce fuel consumption and GHG emissions.Fig. 17 shows that the proposed method takes full advantage of the engine power optimization to reduce power levels at locally "high" sea states.Eventually, the proposed method can help to save about 2% of fuel and GHG emissions than the heuristic method (conventional 2D voyage optimization method), but the saving benefit is rather small in comparison with the actual sailing route.This might be due to the captain's experiences combined with their onboard weather routing system is good enough for navigation in calm weather conditions.

Conclusion
This paper proposed a voyage optimization method combining the benefits of genetic algorithms and the dynamic programming concept to allow for a three-dimensional voyage optimization in terms of minimum fuel consumption and air emissions.Different from other conventional 3D optimization method taking ship speeds as a control parameter for the voyage planning, the big advantage and contribution of the proposed method are that it can provide recommended sailing route and engine power settings along the voyage as direct outputs to guide practical ship navigation.A chemical tanker with full-scale measurements sailed in North Atlantic in the year 2015 and 2016 is taken as a reference to demonstrate the capability of the proposed method.The proposed method is compared with the conventional 2D voyage optimization method (the heuristic method) and the ship's actual sailing for the six case study voyages.By keeping the same ETAs as the measured voyages for the optimization, the proposed method shows great potential with respect to fuel saving and GHG emissions, i.e., on average 5.6% fuel saving and emission reduction than the actual sailing routes and 3.4% to the heuristic routes.It is corresponding to 275 tons GHG emission reduction than the actual six routes.For voyage optimization when sailing in the harsh sea environment, the proposed voyage optimization method can help to reduce about 12.5% fuel and air emissions than the actual sailing.While for voyage planning for ships sailing in a very calm sea environment (maximum significant wave height less than 2.5 m), the proposed method can help to save 2% fuel and GHG emissions than the 2D voyage optimization method by voluntarily adjusting the ship's engine power/speed.
In the meantime, through the evolution of the proposed method, more candidate optimal routes with pre-defined ETAs are generated.Those candidate optimal routes provide many alternative sailing schedules for the decision-makers to fulfil their sailing objectives.This will be very beneficial to combine with slow steaming when more stringent regulations are introduced to reduce air emissions.Furthermore, the case study has also implicitly shown that the proposed method has the capability of voluntary speed reduction under harsh sea conditions to reduce fuel consumption and air emissions along the entire voyage.

Fig. 5 .
Fig. 5.An example of the confined graph generated by the heuristic method.
H.Wang et al.

Fig. 7 .
Fig. 7. Illustration of Pareto front from the optimization and 80% of the Pareto samples are chosen as the basis for the evolution of next population generation.

Fig. 9 .
Fig. 9. Measured voyages used in the case study.

Fig. 10 .
Fig. 10.Comparison of routes/paths optimized by different route planning methods for the case study.

Fig. 11 .
Fig. 11.Results of fuel consumptions and ETAs of various decision vectors during the evolution of the genetic algorithm (note that the iteration 1 is the initial population).

Fig. 12 .
Fig. 12. Evolution of voyage planning along with the contour plot of significant wave height (H s ) of surrounding sea states along the Voyage 20160317.

Fig. 13 .
Fig.13.Shaft power configurations and the encountered H s along the Voyage 20,160,317.

Fig. 14 .
Fig.14.Evolution of voyage planning along with the contour plot of significant wave height (H s ) of surrounding sea states along the Voyage 20160607.

Fig. 15 .
Fig.15.Shaft power configurations and the encountered H s along the Voyage 20160607.

Fig. 16 .
Fig. 16.Original and confined graphs by the proposed optimization method for the Voyage 20160806.
H.Wang et al.

Fig. 17 .
Fig. 17.Shaft power configurations and the encountered H s along the Voyage 20160806.

Table 1
Main characteristics of the chemical tanker.

Table 2
Overall results of various optimization method to the case study voyages.