2.1.1. Basic Concepts

A\* algorithm is a heuristic algorithm widely used for finding an optimal path in static road network presented by [18], which is derived from the Dijkstra algorithm [19] and the Greedy algorithm [20]. The Dijkstra algorithm can find the shortest path, but has to traverse the entire network with low efficiency, and the Greedy algorithm has fast search speed but cannot guarantee to find the best path. The A\* algorithm can balance both search speed and global optimality by using the specific utility function *f*(*n*), which consists of a kind of cost function *g*(*n*) and a kind of cost function *q*(*n*):

$$f(n) = \gcd(n) + q(n) \tag{1}$$

where *g*(*n*) represents the actual cost from initial node to the current node, and *q*(*n*) is the estimated cost from the current node to the end node. When *q*(*n*) = 0, only *q*(*n*) works, then the A\* algorithm degenerates to the Dijkstra algorithm, which can only guarantee finding the optimal route. When *h*(*n*) ≤ *q*(*n*), then the A\* algorithm can maintain the search speed and the global optimality, and the search speed will be slower when the value of *q*(*n*) becomes smaller. When *h*(*n*) *q*(*n*), then the A\* algorithm degenerates to the Greedy algorithm, which can run faster but may fall into local optimum.

## 2.1.2. Work Flow

The flow of the algorithm can be seen in Figure 1.


**Figure 1.** Work flow of the traditional A\* algorithm.

#### *2.2. Basic Concepts of Hesitant Fuzzy Set*

#### 2.2.1. Hesitant Fuzzy Set

Hesitant fuzzy set, proposed by [21], is a more general fuzzy set. An HFS is defined in terms of a function that returns a set of membership values for each element in the domain [21].

**Definition 1** ([21]). *A hesitant fuzzy set A on X is a function h<sup>A</sup> that when applied to X returns a finite subset of [0,1], which can be represented as the following mathematical symbol:*

$$A = \{ \langle \mathbf{x}, h\_A(\mathbf{x}) \rangle | \mathbf{x} \in X \}, \tag{2}$$

*where hA*(*x*) *is a set of some values in [0,1], denoting the possible membership degrees of the element x* ∈ *X to the set A. For convenience, hA*(*x*) *is named a hesitant fuzzy element (HFE)* [22]*.*

**Definition 2** ([21]). *For a hesitant fuzzy set represented by its membership function h* = *hA*(*x*), *we define its complement as follows:*

$$h^c = \cup\_{\mathcal{I}^c \in h} \{1 - \gamma\}.\tag{3}$$

**Definition 3** ([22]). *For an HFE h, Sc*(*h*) = 1*lh* ∑*γ*∈*<sup>h</sup> γ, is called the score function of h, where lh is the number of elements in h and Sc*(*h*) ∈ [0, 1]*. For two HFEs h*1 *and h*2*, if Sc*(*h*1) > *Sc*(*h*2)*, then h*2 ≺ *h*1*, if Sc*(*h*1) = *Sc*(*h*2)*, h*2 ≈ *h*1.

Some operations on the HFEs (*h*, *h*1 and *h*2) and the scalar number *λ* are defined by [22]:

$$h\_1 \oplus h\_2 = \cup\_{\gamma\_1 \in h\_1, \gamma\_2 \in h\_2} \{\gamma\_1 + \gamma\_2 \gamma\_1 \gamma\_2\},\tag{4}$$

$$h\_1 \odot h\_2 = \cup\_{\gamma\_1 \in h\_1 \cap \gamma\_2 \in h\_2} \{\gamma\_1 \gamma\_2\},\tag{5}$$

$$h^{\lambda} = \cup\_{\gamma \in h} \left\{ \gamma^{\lambda} \right\}\_{\prime} \tag{6}$$

$$h^{\lambda} = \cup\_{\gamma \in h} \{\gamma^{\lambda}\}.\tag{7}$$

#### 2.2.2. The Aggregation Operators for Hesitant Fuzzy Information

Reference [23] proposed an aggregation principle for HFEs:

**Definition 4** ([23]). *Let A* = {*h*1, *h*2, ..., *hn*} *be a set of n HFEs,* Θ *be a function on A,* Θ : [0, 1]*<sup>N</sup>* → [0, 1] *, then*

$$\Theta\_A = \cup\_{\gamma \in \{h\_1 \times h\_2 \times \dots \times h\_n\}} \{\Theta(\gamma)\}. \tag{8}$$

Based on Definition 4, some new aggregation operators for HFEs were given by [22]:

**Definition 5** ([22]). *Let hi (i* = 1, 2, ... , *n) be a collection of HFEs, w* = (*<sup>w</sup>*1, *w*2,..., *wn*)*<sup>T</sup> be the weight vector of them, such that wi* ∈ [0, 1] *and* ∑*ni*=<sup>1</sup> *wi* = 1*. A generalized hesitant fuzzy weighted averaging (GHFWA) operator is a mapping Hn* → *H , and*

$$\text{GHFWA}\_{\lambda}(h\_1, h\_2, \dots, h\_n) = \left(\oplus\_{i=1}^n \left(w\_i h\_i^{\lambda}\right)\right)^{\frac{1}{\lambda}} = \cup\_{\gamma\_1 \in b\_1, \gamma\_2 \in b\_2, \dots, \gamma\_n \in b\_n} \left\{ \left(1 - \prod\_{l=1}^n \left(1 - \gamma\_l^{\lambda}\right)^{w\_l} \right)^{\frac{1}{\lambda}} \right\}.\tag{9}$$

**Definition 6** ([22]). *Let hi (i* = 1, 2, ... , *n) be a collection of HFEs, w* = (*<sup>w</sup>*1, *w*2,..., *wn*)*<sup>T</sup> be the weight vector of them, such that wi* ∈ [0, 1] *and* ∑*ni*=<sup>1</sup> *wi* = 1*. A generalized hesitant fuzzy weighted geometric (GHFWG) operator is a mapping Hn* → *H , such that*

$$\begin{aligned} \text{GHFNG}\_{\lambda}(h\_1, h\_2, \dots, h\_n) &= \frac{1}{\lambda} \Big( \otimes\_{i=1}^n (\lambda h\_i)^{w\_i} \Big) \\ &= \cup\_{\gamma\_1 \in h\_1, \gamma\_2 \in h\_2, \dots, \gamma\_n \in h\_n} \Big\{ 1 - \left( 1 - \prod\_{i=1}^n \left( 1 - (1 - \gamma\_i)^{\lambda} \right)^{w\_i} \right)^{\frac{1}{\lambda}} \Big\}. \end{aligned} \tag{10}$$

#### 2.2.3. Decision Making Based on Hesitant Fuzzy Information

The decision method based on the above definitions can be derived as follows:

**Step 1** The possible alternative *Xi* of the attribute *Ai* provided by decision makers or other sources are denoted by the hesitant fuzzy elements *hi*(*i* = 1, 2, . . . , *<sup>n</sup>*).


**Step 4** Choose the optimal alternative *X*∗ by the comparison of *Sc*(*hi*).

**Example 1.** *The vessel, which includes four directions Xi*(*i* = 1, 2, 3, 4) *to navigate, is planed to determine the optimal Arctic route for the following year. Suppose there are three factors Ai*(*i* = 1, 2, 3) *that affect the decision making—A*1*: navigation time; A*2*: economic cost; A*3*: navigation risk. It should be noted that all of them are of the minimization type. The weight vector of the attributes is w* = (0.3, 0.4, 0.3)*<sup>T</sup>*.

Then, the optimal route can be determined by using the mentioned method.


$$\begin{aligned} h\_1^A &= G HFWA\_2(h\_{11}, h\_{12}, h\_{13}) = \left( \oplus\_{j=1}^3 \left( w\_j h\_{1j}^2 \right) \right)^{\frac{1}{2}} \\ &= \cup\_{\gamma\_{11} \in h\_{11}, \gamma\_{12} \in h\_{12}, \gamma\_{13} \in h\_{13}} \left\{ \left( 1 - \prod\_{j=1}^3 \left( 1 - \gamma\_{1j}^2 \right)^{w\_j} \right)^{\frac{1}{2}} \right\} \end{aligned}$$

= {0.2157, 0.2652, 0.3503, 0.4351, 0.4577, 0.5040, 0.5966, 0.6099, 0.6383, 0.2630, 0.3039, 0.3789, 0.4566, 0.4788, 0.5213, 0.6093, 0.6220, 0.6492, 0.3166, 0.3503, 0.4149, 0.4847, 0.5040, 0.5442, 0.6263, 0.6383, 0.6639}.

$$\begin{aligned} h\_1^G &= GHFWG\_2(h\_{11\prime}, h\_{12\prime}, h\_{13}) = \frac{1}{2} \left( \otimes\_{j=1}^3 \left( 2h\_{1j} \right)^{w\_j} \right) \\ &= \cup\_{\gamma 11 \in h\_{11\prime}\gamma 12 \in h\_{12\prime}\gamma 13 \in h\_{13}} \left\{ 1 - \left( 1 - \prod\_{j=1}^3 \left( 1 - \left( 1 - \gamma\_{1j} \right)^2 \right)^{w\_j} \right)^{\frac{1}{2}} \right\} \end{aligned}$$

= {0.1814, 0.2541, 0.2915, 0.2672, 0.3857, 0.4514, 0.2848, 0.4142, 0.4875, 0.1958, 0.2754, 0.3168, 0.2898, 0.4226, 0.4983, 0.3092, 0.4551, 0.5411, 0.2065, 0.2915, 0.3361, 0.3071, 0.4514, 0.5361, 0.3279, 0.4875, 0.5851}.

**Step 3** The score values Sc(*hi*), *i* = 1, 2, 3, 4 are calculated by **Definition 3,** which can be seen in Table 2. **Step 4** From Table 2, *X*4 will be chosen as the optimal direction based on both the GHFWA operator and GHFWG operator where *λ* is set as 2.

**Table 1.** Hesitant fuzzy decision matrix.


**Table 2.** Score values obtained by the GHFWA operator and GHFWG operator and the rankings of alternatives (*λ* = 2).


#### **3. An Improved A\* Algorithm (HFS-A\*)**

In light of the harsh weather conditions in the Arctic region, the primary task for route planning is to identify the obstacles (e.g., sea ice). The Coupled Model Intercomparison Project, phase 5 (CMIP5) provided 39 Global Climate Models (GCMs) to predict sea ice data, from history to 21st century, under different representative concentration pathways (RCPs) [23,24]. Unlike the route planning in other regions, for the current sea ice forecasts in the Arctic region, there exists large uncertainty among these GCMs [25,26], which leads to the uncertainty of the length of the navigation season and the economic risk of exploiting the Arctic routes [27]. Therefore, only treating the shortest distance as the optimal route in the Arctic region is not reasonable; more factors, including the navigation risk, the navigation time, and the economic cost during navigation should be considered. Compared to the traditional A\* algorithm, the HFS-A\* algorithm is used to tackle the multi-criteria decision-making (MCDM) problem with large uncertainty derived from multi-model outputs and expert knowledge. The improved parts mainly focus on t, the identification of obstacles, and the construction of utility function.

#### *3.1. Navigability of the Arctic Routes*

With the impact of global warming, the extent of Arctic sea ice continues to decline [27]. Human's enthusiasm to explore and develop the Arctic routes are aroused by shorter sailing distance, longer navigation season and increased access to natural resources. There are three criteria related to sea ice conditions for evaluating the navigability in the Arctic area.

**Criterion 1** (navigation uncertainty)**.** *Sea ice concentration is considered only for no ice-breaking or ice-strengthening ships, and it is navigable when sea ice concentration is less than 15%* [28–31]*.*

**Criterion 2** (navigation time)**.** *Sea ice thickness derived from an empirical regression model is considered and for no ice-breaking or ice-strengthening ships, it is navigable when sea ice thickness is no more than 1.2 m* [32,33]*.*

**Criterion 3** (navigation economic cost)**.** *Both sea ice concentration and thickness are considered by computing the Ice Numeral (IN) index from the Arctic Ice Regime Shipping System (AIRSS) provided by the Canada Transport* [34–37]. *The Ice Numeral is given by*

$$IN = \mathbb{C}\_a IM\_a + \mathbb{C}\_b IM\_b + \dots + \mathbb{C}\_n IM\_n \tag{11}$$

*where Cn is the concentration in tenths of ice type n, and IMn is the Ice Multiplier for ice type n*. *Ice type describes the specific stage of development of ice, which is closely related to the ice age. Ice Multipliers, determined by ship class and ice type, are a series of integers, which are used to reflect the impact of sea ice type to the specific vessel. A negative IM represents the obstacle effect of vessel sailing. Ice types are determined by* [34,38]*, which are presented in Appendix A. For no ice-breaking or ice-strengthening ships, it is navigable when the IN index is larger than zero. Details about vessel type and IM can also be seen in Appendix A.*

Additionally, geographical environment, including water depth and channel width, is also a key factor for ships to navigate, which is related to the vessel type and dimension.

Overall, various evaluation criteria will add to the uncertainty of the sea ice navigability projection. In this paper, we take all three criteria into consideration with geographical restriction to make sure the navigability of the Arctic region, in another words, for one region, can be defined as navigable if and only if the criteria mentioned above are all reached.

#### *3.2. Route Planning Criterion 1: Uncertainty of Sea Ice Condition*

In this paper, a series of GCMs have been chosen based on their reasonable projections for future sea ice conditions evaluated by literature studies [39–42] (see Appendix B). In order to obtain the uncertainty of each model, these model outputs were compared with the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) estimate data set, which is a reanalysis with good spatial and temporal consistency constrained by the quality of the assimilated observations [43–46].

Suppose we have *M* GCMs, let each model data-set be X*i*, *i* = 1, 2, ... , *M*, the PIOMAS data set be *X*, then the uncertainty of each model can be obtained as follows:

$$h\_{i1} = 0.5 - 0.5 \times \frac{\chi\_i \bar{X}}{\|X\_i\| \times \|\bar{X}\|} , \quad i = 1, 2, \dots, M \tag{12}$$

where, *Xi*·*X Xi*×*X* is called cosine similarity, which is a measure of similarity between two non-zero vectors introduced by [47]. *h*i1, *i* = 1, 2, ... , *M* represents the bias from the model to the "real state", which can also be regarded as model uncertainty. When *h*i1 is equal to zero, it represents that the model data sets can well reflect the "real state", while when *h*i1 is equal to one, the model uncertainty reaches its maximum.

#### *3.3. Route Planning Criterion 2: Time for Navigation on the Arctic Routes*

The navigation time for each grid can be depicted as follows:

$$T\_i = \frac{S\_i}{V\_i'} \text{ i } = 1, 2, \dots, M\_r \tag{13}$$

where, *Si* is the distance of each grid, and *vi* is the velocity of the vessel on each grid. *hi*2 is the normalization of *Ti*.

On the Northern Sea Route (NSR), the vessel speed is mainly impacted by sea ice conditions. A h-v curve presented by [48] can reflect, well, the relationship between sea ice thickness and vessel speed. In this h-v curve, the ice resistance and the net thrust of the engine to overcome the ice resistance should both be considered.

**Step 1** The ice resistance can be presented as follows:

$$\begin{split} \mathbf{R\_{ch}} &= 0.5 \mu\_{B} \rho\_{\Delta} \mathbf{g} \mathbf{H}\_{F}^{2} \mathbf{K\_{F}} \left( \frac{1}{2} + \frac{H\_{\rm M}}{2H\_{\rm F}} \right)^{2} \left[ \mathbf{B} + 2H\_{\rm F} (\cos \delta - 1/\tan \psi) \right] (\mu\_{h} \cos \phi + \sin \psi \sin \alpha) \\ &+ \mu\_{B} \rho\_{\Delta} \mathbf{g} \mathbf{K\_{0}} \mu\_{h} L\_{\rm par} H\_{\rm F}^{2} + \rho\_{\Delta} \mathbf{g} \left( \frac{17}{\mathcal{B}^{2}} \right)^{3} H\_{\rm M} A\_{\rm WF} F\_{\rm n}^{2} \end{split} \tag{14}$$

where *μB* = 0.5 represents the ice porosity factor, *ρ*Δ = 0.8 is the density difference between ice and water. In this model, these two factors are considered to be constant for simple situations, e.g., when the temperature changes. Variables include the waterline area of the foreship *AWF*, the Froude number *Fn*, the length *L*, the parallel midbody length at waterline *Lpar*, the width *B*, the friction coefficient *μh* and the vessel draft *T*. *KP* and *K*0 are mechanical factors of ice found by [49]. The thickness of the brash ice *HF* can be determined as follows:

$$H\_F = H\_M + \frac{B}{2}tan\gamma + (tan\gamma + tan\delta)\sqrt{B\left(\frac{H\_M + \frac{B}{4}tan\gamma}{tan\gamma + tan\delta}\right)}.\tag{15}$$

Both *γ* and *δ* represent the slope angles of the brash ice.

The formula can be simplified by an approximation when *B* > 10 m and *HM* > 4m[39]:

$$H\_F = 0.26 + (BH\_M)^{0.5}.\tag{16}$$

The flare angle *ψ* mentioned above can be obtained with the bow angles *φ* and *α*:

$$
\psi = \arctan\left(\frac{\tan\phi}{\sin\alpha}\right). \tag{17}
$$

**Step 2** The net thrust can be calculated by the following formula:

$$T\_{tot}(v)(1-t) = R\_{ow}(v) + R\_i(v),\tag{18}$$

$$R\_{tot}(v)(1-t) = R\_{av}(v) + R\_i(v),\tag{19}$$

where *Ttot* is the total thrust, (1 − *t*) represents the thrust deduction factor, and *Ri*(*v*) and *Row*(*v*) are the resistance in ice and in open water respectively.

The effect of vessel speed is approximated by a quadratic factor called bollard pull *TB* [50], with the maximum open water speed *vow*, and the net thrust can be rewritten as

$$T\_{nct}(\upsilon) = \left(1 - \frac{\frac{1}{3}\upsilon}{\upsilon\_{ow}} - \frac{2}{3}\left(\frac{\upsilon}{\upsilon\_{ow}}\right)^2\right)T\_{B\prime} \tag{20}$$

$$T\_B = K\_c (P\_D D\_P)^{\frac{2}{3}},\tag{21}$$

where *v* is the vessel speed in the ice, *Ke* is the bollard pull quality factor, *DP* is the propeller diameter and *PD* is the actual power delivered.

#### *3.4. Route Planning Criterion 3: Economic Cost for Navigation on the Arctic Routes*

The economic cost of navigation on the Arctic routes are consist of four parts, which are capital cost, fuel cost, operation cost, and transit cost:

$$\text{Cost}\_{\text{economic}} = \text{Cost}\_{\text{capital}} + \text{Cost}\_{fuel} + \text{Cost}\_{\text{operation}} + \text{Cost}\_{\text{transit}} \tag{22}$$

#### 3.4.1. Model for Capital Cost

Capital cost is related to the price of new ship building or the cost of a ship with loans and depreciations [51]. Generally, the annual capital cost can be computed as follows [52]:

$$\text{Cost}\_{capital} = \left( P \times (1 - \epsilon q) \times \frac{(1+r)^y \times r}{(1+r)^y - 1} \right) + (\epsilon q \times \mathbb{C}) / y \tag{23}$$

where *P* is a new ship price, *eq* is the equity, *r* is the interest rate, and *y* is the term of loan. The former item is called annual interests, while the latter item is called cash price.

#### 3.4.2. Model for Fuel Cost

Fuel cost is often depicted as the major cost on the marine transportation, which is the largest single cost factor to most related simulations, ranging from 36.7% to 61% [53].

Fuel cost is impacted by the rate of fuel consumption and the price of fuel, which can be descripted as follows:

$$Cost\_{fuel} = Fuel\_{consumption} \times Fuel\_{price} \tag{24}$$

The fuel consumption of a vessel can be influenced by ship dimensions (e.g., the ship size, hull design, engine profile, speed) and external factors (e.g., sea ice, wind, wave, current, foggy) [51]. The fuel consumption for a specific type of vessel is basically computed by multiplying SFOC (specific fuel oil consumption) (g/kWh), engine power (kw) and sailing hours (h) [54].

Most authors consider speed as the key factor impacting the fuel consumption when the type of a vessel is determined and a simple exponential law based on empirical data are derived: the fuel consumption per unit distance is proportional to the square of the speed [55–57], which can be presented as

$$\frac{V\_1^2}{F\_1} = \frac{V\_2^2}{F\_2} \, \, \, \, \tag{25}$$

where *F*1 and *F*2 are the fuel consumption rate under the velocity *V*1 and *V*2, respectively. Fuel price is affected by the fluctuation of the global economic market. Low fuel price indicates the depression of global economy, while high fuel price reflects the booming global market.

#### 3.4.3. Model for Operation Cost

Operation cost mainly includes crew cost, insurance cost and maintenance fee which is presented as follows:

$$\text{Cost}\_{operation} = \text{Cost}\_{crw} + \text{Cost}\_{maintenance} + \text{Cost}\_{insurance} \tag{26}$$

## 1. Crew cost

Crew cost is determined by the vessel type, automation level and numerous other factors [51]. Compared to the open water, the crews in the Arctic region require additional ice navigation experience and the ability to cope with hash weather conditions, which may increase the crew cost [53]. The increased crew cost may come from the higher wage [58] for each member or the larger size of members [59].

## 2. Maintenance cost

For the purpose of preventing the occurrence of breakdowns and following the scheduled maintenance program, the cost of regular maintenance is needed for vessels.

## 3. Insurance cost

In face of the risk of Arctic navigation (e.g., collision, engine damage, propeller damage, local hull damage, grounding, etc.) analyzed IN some studies [60,61], maritime insurance is a good tool to mitigate the associated risks, which can be approximately separated into three major components: protection and indemnity (P&I), hull and machinery (H&M), and cargo insurance. The third-party liabilities encountered during the commercial operation of a ship are charged by P&I. H&M covers the cost of damage done to the ship or its equipment. Cargo insurance provides the paymen<sup>t</sup> for the damage to the cargo itself [62].

#### 3.4.4. Model for Transit Cost

On the NSR, the transit fee based on vessel type and ice conditions mainly includes ice pilot fee and ice breaking fee and can be given as

$$\text{Cost}\_{transit} = \text{Cost}\_{iccbreaking} \times \text{Load} \times \text{Pay}\_{load} + \text{Cost}\_{icepilit} \tag{27}$$

These services are mainly provided by the Russian icebreaking service provider, Atomflot, and are compulsively charged subject to the law of the Russian Federation, which is dependent on the vessel type (e.g., the size and the ice-class of a vessel), and the navigation length and pilotage distance [62]. In general, higher ice-classed vessels are charged with lower icebreaking fees.

#### *3.5. Work Flow of HFS-A\* Algorithm*

In this HFS-A\* algorithm, the work flow can be seen in Figure 2. The method of obstacle identification has been discussed in Section 3.1 while the modified utility function can be described as follows:

$$f^\*(i) = \mathbf{g}^\*(i) + \mathbf{q}^\*(i) \tag{28}$$

When the current node *i* is determined, the actual cost *g*<sup>∗</sup>(*i*) is equal to the score value Sc(*hi*), which can be computed by Definition 3.

More specifically, we assume each selected node has three criteria *Cj*(*j* = 1, 2, 3) that affect the decision making— *C*1: uncertainty of sea ice condition; *C*2: navigation time; *C*3: navigation economic cost. It should be noted that all of them are of the minimization type. The weight vector of the attributes is *w* = (*<sup>w</sup>*1,..., *wj*)*<sup>T</sup>*.

**Figure 2.** Work flow of the HFS-A\* algorithm.

The heuristic estimated cost function can be approximately evaluated:

$$q^\*(i) = \mathcal{S}c(h\_i) \times D \tag{29}$$

where, *D* is the heuristic distance (Manhattan, Euclidean or Chebyshev) from the evaluated node to the END node [63].

**Step 1** Map initialization


**Step 2** The construction of utility function

• Each time, compare all the adjacent nodes *i* of the current node *n* by

$$f^\*(i) = \mathcal{Sc}(h\_i) \times (D+1) \tag{30}$$

where, *Sc*(*hi*) = *GHFWGλ*(*hi*1, *hi*2, *hi*3) or *Sc*(*hi*) = *GHFWAλ*(*hi*1, *hi*2, *hi*3), *hi*1 is the HFEs of navigation uncertainty (see Section 3.2), *hi*2 is the HFEs of navigation time (see Section 3.3), and *hi*3 is the HFEs of navigation economic cost (see Appendix **??**).

**Step 3–6** The same as the traditional A\* algorithm.

#### **4. Case Study and Conclusion**

#### *4.1. Study Area and Data Description*

This experiment is to find the optimal route on the Northern Sea Route (NSR) based on the proposed method from Shanghai to Bergen port for an IB-classed 3800TEU container vessel.

Data related to water depth is derived from a product called ETOPO1 provided by the National Geophysical Data Center (NGDC), with a resolution of 1 arc-minute [64]. Data related to AIRSS system can be seen in Appendix A. Data related to sea ice conditions (both sea ice thickness and sea ice concentration) can be seen in Appendix B. Data related to vessel information can be seen in Appendix C. Data related to economic cost can be seen in Appendix D. All the climate model outputs and data related to water depth are interpolated to the grid size of 360 × 120 for the comparison with PIOMAS estimate data set, which the spatial coverage is 45◦ N to 90◦ N and the temporal resolution is monthly.

#### *4.2. Route Planning by HFS-A\* Algorithm*

#### 4.2.1. Navigability of the NSR

The numerical simulation firstly examines the navigability of the IB-classed 3800 TEU container vessels on the NSR for each month in the year of 2050 (see Figure 3). According to the ensemble model predictions, the open time of the NSR for that vessel to access may last for 3 to 5 months in the year of 2050. Most model outputs show the navigable time starts from August to the October, while merely 2 to 4 models extend the navigable time (from July to November).

**Figure 3.** The navigability of IB-classed 3800 TEU container vessels on the NSR for each month in the year of 2050 based on multi-models. (The color of orange in the map reflects the geographic information, the white color represents the area that cannot access during that month, different blue colors reflect different amount of the models that give the navigable prediction for each grid during that month.).

#### 4.2.2. Selection of the Aggregation Operators and Route Optimization

Secondly, route planning criteria (uncertainty, time and economic cost) have been calculated based on the models mentioned in Sections 3.2–3.4, and the results can be seen in Appendix D. These factors can be normalized respectively between 0 and 1, which can be described as:

$$\frac{h\_i - \min(h\_i)}{\max(h\_i) - \min(h\_i)}, \text{ } i = 1, 2, \dots, M \tag{31}$$

Thirdly, the optimal routes for IB-classed 3800 TEU container vessels from Shanghai to Bergen port on the NSR in September of 2050 can be found by two kinds of the aggregation operators (GHFHA*<sup>λ</sup>*, *λ* = 1, 2, 6 and GHFHG*λ*, *λ* = 1, 2, 6), which are used to aggregate the normalized results mentioned above. The detailed results can be seen in Table 3. Compared these two kinds of aggregation operators, the performance of the GHFHG*λ*, *λ* = 1, 2, 6 operators in the route

optimization is better than the GHFHA*λ*, *λ* = 1, 2, 6 operators from the view of total sailing distance, sailing time, economic cost and average uncertainty. In the light of the comparison of different *λ* for each operator, the performance of the GHFHG*λ* becomes better with the *λ* increase, while the performance of the GHFHG*λ* becomes better when *λ* decreases. Therefore, the GHFHG1 operator has been examined as the best aggregation operator in this numerical study (The weights vectors for this experiment are assigned as 0.4, 0.3, and 0.3 for three criteria).

**Table 3.** Results of route optimization for IB-classed 3800 TEU container vessels on the NSR from Shanghai to Bergen port in September of 2050 with different aggregation operators based on HFS-A\* algorithm.


Finally, the optimal route determined by HFS-A\* algorithm with the GHFHG1 operator can be seen in Figure 4. Three other routes according to simple single criterion (uncertainty, time, or economic cost) are also drawn in Figure 4. The detailed results can be seen in Table 4, where it can be found that path planning based on a single factor shows a slight advantage in its related aspect, but shows significant disadvantage in any other aspect compared with the optimal route. In other words, the optimal route can better balance these three key factors and show more realistic performance of the proposed route planning algorithm than the other three single factor route planning.

**Table 4.** Results of route optimization for IB-classed 3800 TEU container vessels on the NSR from Shanghai to Bergen port in September of 2050 based on HFS-A\* algorithm and simple A\* algorithm (The percentages in brackets are compared with the values in GHFHG1).


**Figure 4.** Route optimization for IB-classed 3800 TEU container vessels on the NSR from Shanghai to Bergen port in September of 2050 by HFS-A\* algorithm. (The red line represents the route planning based on navigation uncertainty, the blue line based on the navigation time, the yellow line based on navigation economic cost, the dark line represents the optimal route integrated of these three criteria by the GHFHA1 operator.)
