1. INTRODUCTION
The shipping industry has always been at the forefront of the world trade market, which is the barometer of international economic trends. Recently, due to a multitude of pressures such as the impact of natural disasters, geopolitical conflict, pirate attacks, and so on, shipping companies have started to seek alternate routes, in an attempt to reduce costs and increase profitability. In the late Twentieth Century, Arctic sea ice began to melt more quickly as a result of global warming. The Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report showed that Arctic sea ice may completely melt during summer in the middle of the Twenty-First Century and in extreme cases during winter by the end of 2070 (Stocker, Reference Stocker2014), which makes Arctic navigation no longer just a dream. The Arctic sea routes considered in this paper are the Northern Sea Route (NSR), the Northwest Passage (NWP) and the Central Passage (see Figure 1). Compared with the traditional sea routes, for many destinations, the Arctic routes have a shorter shipping distance, wider routes and fewer pirates, which will lead to more common usage.
Figure 1. The Arctic marine area.
The question posed by the Arctic Council in the Arctic Marine Shipping Assessment 2009 Report aims to find the optimal shipping routes considering both the Arctic and non-Arctic routes. In terms of sailing time, the distance from Asia (China or Japan) to the east coast of North America through the NWP for general types of ships can save 25% to 44% of total distance compared to routes using the Panama Canal, while navigation through the NSR from Asia to Europe can save 25% to 55% of total distance compared to using the Suez Canal (Andersen et al., Reference Andersen, Heggeli and Wergeland1995; Ragner, Reference Ragner2000; Zhang et al., Reference Zhang, Tu, Guo, Sun and Ling2009). In view of economic benefits, Somanathan et al. (Reference Somanathan, Flynn and Szymanski2009) simulated the NWP and the Panama Canal by computer and the result showed the economic potential of the NWP would gradually rise with the decline in the extent of sea ice. The NSR will become an important alternative route to the Suez Canal (Verny and Grigentin, Reference Verny and Grigentin2009; Chernova and Volkov, Reference Chernova and Volkov2010; Srinath, Reference Srinath2010; Liu and Kronbak, Reference Liu and Kronbak2010). Way et al. (Reference Way, Khan and Veitch2015) presented a probability simulation method to compare the optimal sailing speed of container ships on the NSR or on the Suez Canal based on the model presented by Liu and Kronbak (Reference Liu and Kronbak2010). When considering both sailing time and distance, a comparison of shipping costs was made between the NSR and the Suez Canal and results showed that the rise of the NSR will promote the agility and adaptability of the global shipping supply chain (Schøyen and Bråthen, Reference Schøyen and Bråthen2011).
In order to further enhance the competitiveness of Arctic routes, route optimisation issues have been studied by some researchers. The safe navigation speed was determined by the Ice Number (IN) index (McCallum, Reference McCallum1996) which became the criteria to plan the optimal route (Smith et al., Reference Smith and Stephenson2013; Melia et al., Reference Melia, Haines and Hawkins2016). Kotovirta et al. (Reference Kotovirta, Jalonen, Axell, Riska and Berglund2009) took sea ice thickness, ice ridge thickness and sea ice concentration into consideration to analyse their impact on ship speed. Choi et al. (Reference Choi, Chung, Yamaguchi and Nagakawa2015) presented a route planning model with a heuristic algorithm in which both extreme conditions and model uncertainty were considered. Other environmental factors (for example, wind, waves, current and topography) and economic factors were also taken as variables in some integrated models (Motte et al., Reference Motte, Burns and Calvert1988; Mulherin et al., Reference Mulherin, Eppler, Proshutinsky, Farmer and Smith1996; Choi et al., Reference Choi, Nam, Park, Ha and Jeong2010; Hogben and Lumb, Reference Hogben and Lumb1967; Ha et al., Reference Ha, Choi, Nam and Park2011; Nam et al., Reference Nam, Park, Lee, Mi, Choi and Seo2013; Aksenov et al., Reference Aksenov, Karcher, Proshutinsky, Gerdes, De Cuevas, Golubeva, Kauker, Nguyen, Platov, Wadley, Watanabe, Coward and George Nurser2016).
In the literature mentioned above, both Arctic and non-Arctic routes were assumed to be independent of each other and no interaction between them is considered in route optimisation. However, an obvious competitive relationship actually exists between the routes relative to any Origin-Destination (OD) pair including the competition among shipping companies and the interaction between shipping companies and customers. Therefore, when the potential competitiveness of the Arctic routes is considered, the impact of the non-Arctic routes should also be considered.
For the purpose of finding a quantitative model to illustrate the interaction between Arctic routes and the traditional routes, game theory is a powerful tool to characterise the elements of dynamic competition. This paper proposes a bi-level model to demonstrate the effect of the rise of Arctic routes on the traditional routes and the response of traditional routes. The model illustrates the entire competition process among shipping companies over different sea routes and the customers during two equilibrium stages, which can be presented in a Stackelberg form (Lam and Zhou, Reference Lam and Zhou2000). The upper-level model illustrates the competition among shipping companies who cut their own cost by speed optimisation and maximise their own profits, which is presented as a Generalised Nash Equilibrium Problem (GNEP) and can be transformed into a type of Quasi-Variational Inequality (QVI) problem. Harker (Reference Harker1991) pointed out that solving this QVI problem is equivalent to finding a Variational Inequality (VI) solution, which can be solved by the generalised reduced-gradient method (Lasdon et al., Reference Lasdon, Waren, Jain and Ratner1978). The lower-level model illustrates the responses of customers who reassign their demands on each shipping route with an elastic demand, which is presented as a logit-type multi-path assignment problem and can be solved using the balancing iteration method.
In the case study, the proposed model used to simulate a practical problem illustrates the interaction between these two types of routes. Two stages (before and after the rise of the Arctic routes) of one Origin to Destination (OD) pair (no stop on the way) are compared to demonstrate the effect of the rise of the Arctic routes on the traditional routes and the response of the traditional routes. The harsh weather conditions in the Arctic region are considered to modify the results of the bi-level model so as to make the simulation more realistic. The structure of this paper is as follows: Section 2 introduces some basic concepts and notations. In Section 3 this improved bi-level game model is introduced and its mathematical characteristics are analysed. The algorithm is introduced in Section 4. Section 5 presents a case study using the proposed model and practical climatic and hydrological data to reveal the interaction between the Arctic and traditional routes. Finally, the conclusion is presented in Section 6.
2. BASIC CONCEPTS
2.1. Total travel cost
The total travel cost Q s in sailing route s generally includes four components: cost of ships, voyage cost, operating cost, and fuel-consumption cost, which can be expressed as:
$$Q_{\rm s}={Pcc}_{s}+{Pvc}_{s}+{Poc}_{s}+{Pfc}_{s}\quad \forall s\in \hbox{S}_{w}$$where Sw represents the set of all sailing routes, Pcc s, Pvc s, Poc s, Pfc s represent the cost of ships, the voyage cost (including the route toll, ice-breaking fee, and pilotage, etc), the operating cost and the fuel consumption cost in sailing route s, respectively, whereas the cost of ships in sailing route s can be expressed as:
$$Pcc_s = Pc_s \times \displaystyle{{T_s} \over {T_a}},\quad \forall s\in S_w$$where Pc s represents the price of one ship in sailing route s, T s represents the travel time of a voyage in sailing route s and T a represents the total travel time in sailing route s.
The operating cost (including manning, hull and machinery insurance, protection and indemnity insurance, repairs and maintenance, administration, and others) in sailing route s can be expressed as:
$${Poc}_{s}={Po}_{s}\times T_{s},\quad \forall s\in \hbox{S}_{w}$$where Po s represents the operating cost per day in sailing route s.
The fuel consumption is determined by the sailing speed. Dykstra (Reference Dykstra2001) defined the relationship as:
$$F_s = \displaystyle{F \over {\hbox{V}^2}}V_s^2 ,\quad \forall s\in \hbox{S}_w$$where F and V represent the standard fuel consumption and standard sailing speed while F s and V s represent the fuel consumption and sailing speed in sailing route s.
On the other hand, the fuel-consumption cost in sailing route s can be obtained as:
$$\hbox{Pf}\hbox{c}_{\rm s} = \displaystyle{{\hbox{F} \times \hbox{V}_{\rm s}^2 } \over {\hbox{V}^2}} \times \hbox{S}\hbox{D}_{\rm s} \times \hbox{Pf},\quad \forall s\in \hbox{S}_{\rm w}$$where SD s represents the distance of sailing route s and Pf represents the fuel price.
Meanwhile, the total travel time of a voyage can be expressed as
$${T}_{s}={Tw}_{s}+{Tv}_{s},\quad \forall s\in \hbox{S}_{w}$$
$$Tv_s = \displaystyle{{SD_s} \over {V_s}},\quad \forall s\in S_w$$where Tw s and Tv s represent the waiting time of a voyage and the travel time of a voyage in sailing route s.
2.2. Generalised Nash Equilibrium Problem GNEP
GNEP, which was originally presented by Debreu (Reference Debreu1952), is a generalisation of the traditional Nash Equilibrium Problem (NEP). Three factors, players, their strategies and their utility functions, are usually considered in a standard game model. In the NEP, only the utility function of one player can be affected by the other players' strategies; the strategies of one player are independent of the other players' strategies. In the GNEP, the utility function and the set of strategies of one player both depend on the other players' strategies.
Definition: Let us consider a non-cooperative game with n players involved. Each player s ∈ S w is represented by a set of strategies 
$\hbox{Y}_{s}\in \Xi_{s}\subseteq {\open R}^{n}$, a point to set mapping 
$\emptyset_{s}\colon \Xi \to \Xi_{s}$, and a utility function U s : Y → ℝ, where 
$\emptyset=\emptyset_{s}$. The generalised Nash equilibrium 
$\hbox{Y}^{\ast }\in \Xi $ of this game is defined as a point where no one can increase their utility by changing their strategies.
$$U_{s}(\hbox{Y}_{s}^{\ast },\hbox{Y}_{-s}^{\ast })\ge U_{s}(\hbox{Y}_{s},\hbox{Y}_{-s}^{\ast }),\quad \forall \hbox{Y}_{s}\in \Xi_{s}(\hbox{Y}^{\ast })$$Thus, Ξ is a constraint set for all players.
Proposition 1:
Let us suppose that each player s ∈ S w satisfies the conditions that:
$\Xi_{s}\subseteq {\open R}^{n}$ is a non-empty, compact, and convex set,- $\emptyset_{s}\colon \Xi \to \Xi_{s}$ is a non-empty and closed map, and
U s : Y →ℝ is continuously differentiable and pseudo-concave.
Then, 
$\hbox{Y}^{\ast }\in \Xi \lpar \hbox{Y}^{\ast }\rpar $ is an optimal solution to the GNEP if and only if 
$\hbox{Y}^{\ast }\in \Xi \lpar \hbox{Y}^{\ast }\rpar $ is an optimal solution to the following QVI problem:
$$\sum_{s\in S_{w}} {{F_{s}(Y^{\ast})}^{T}(\hbox{Y}_{s}-\hbox{Y}_{s}^{\ast })} \ge 0,\quad \forall \hbox{Y}_{s}\in \Xi_{s}(\hbox{Y}^{\ast})$$Thus, F s(Y) can be expressed as:
$$F_{s}\lpar Y \rpar =-\nabla_{{\rm Y}_{s}}U_{s}(Y)$$3. BI-LEVEL MODEL WITH SPEED OPTIMISATION
To better clarify the proposed model, it is illustrated from two aspects. At the upper-level, shipping companies would obviously aim at maximising their own profits, which are mainly determined by the freight volume and unit cargo profits linked to unit freight revenue and travel costs. In recent years, the shipping industry has continued to slump; it is currently facing the pressure of weak trade growth and excess capacity. Therefore, an increasing number of scholars have started to study shipping-cost optimisation, which includes route selection, fleet management, and network design, to reduce freight costs and enhance economic competitiveness (Bijlsma, Reference Bijlsma2008; Tran and Haasis, Reference Tran and Haasis2015; Zhou et al., Reference Zhou, Wang, Song and Peng2016). Among these factors, the network design, which is concerned with the choice of ports and related infrastructure, is not considered in this study. This study focuses on the optimisation of route selection and fleet management. In general, the fuel consumption cost of shipping companies accounts for more than 50% of the total travel cost (Golias et al., Reference Golias, Saharidis, Boile, Theofanis and Ierapetritou2009; Notteboom, Reference Notteboom2006), whereas fuel consumption and sailing speed have a significant relationship (Ronen, Reference Ronen1982; Wang and Meng, Reference Wang and Meng2012; Psaraftis and Kontovas, Reference Psaraftis and Kontovas2013). Therefore, to cut the total cost, the fuel consumption needs to be reduced by adjusting the sailing speed or choosing a shorter route (Corbett et al., Reference Corbett, Wang and Winebrake2009; Lang and Veenstra, Reference Lang and Veenstra2010; Norstad et al., Reference Norstad, Fagerholt and Laporte2011). At the lower-level, the volume of freight in each route depends on the choice of the customers, which is a response to the results of the competition between the shipping companies. In addition, if the external environment (such as world economy, politics, law and technology) is assumed to remain unchanged, the total customer demand is related not only to the revenue but also to the travel time (Wang et al., Reference Wang, Meng and Liu2013).
To clarify this issue, this paper proposes a bi-level model with speed optimisation in the Stackelberg form in which the shipping companies are presented as leaders while the customers are presented as followers. The upper-level model studies the competition among the shipping companies, whereas the lower-level model studies the change in the total customer demand with the change in their shipping cost and the reassignment of customer demand in each route.
3.1. Lower-level model
For the customers, the actual shipping cost in sailing route s is composed of the charges imposed by shipping companies and the time cost of navigation, which can be obtained as:
$$C_{s}=\lambda P+{T}_{s},\quad \forall s\in S_{w}$$where λ is the weight coefficient that reflects the relationship between the time and economic costs in the actual shipping cost in sailing route s.
It is assumed that the customer demand monotonically decreases with the increase in the actual shipping cost, which depends on the distance of the routes and sailing speed of the ships when other factors are fixed. For the shipping companies, they can reduce their time cost by speed optimisation in any route. Therefore, the total customer demand can be obtained as follows:
$$D=\hbox{G}\lpar \bar{\hbox{E}} \rpar =D^{0}\exp\,(-\beta \bar{\hbox{E}})$$ D 0 is the baseline customer demand before the rise of the Arctic routes, β is a dispersion coefficient, and 
$\bar{\hbox{E}}$ is the expected minimum shipping cost for the customers, which can be expressed as (see Ben-Akiva and Lerman, Reference Ben-Akiva and Lerman1985):
$${\bar{\hbox{E}}} = E\lpar {min\lcub {C_s} \rcub } \rpar = -\displaystyle{1 \over \theta }\ln \left( {\sum\limits_{s\in S_w} {\exp } \;(-\bar{\theta }C_s)} \right),\quad \forall s\in S_w.$$ 
$\bar{\theta }$ represents the sensitivity of the different types of customers relative to their actual shipping cost.
In fact, the lower-level model describes the change in the total customer demand with the change in the shipping cost in response to the speed adjustment made by the companies.
3.2. Upper-level model
It is assumed that the shipping companies that use the traditional routes have reached equilibrium during competition before the rise of the Arctic routes, which we designate as Stage 1. Due to the rise of the Arctic routes, some companies change their shipping routes from the traditional ones to the Arctic routes. A new equilibrium among shipping companies will be finally reached, which is designated as Stage 2. For convenience, the shipping companies on the same route are classified as a type of player. Then, the problem of studying the competition among n players is simplified, which is equal to the number of elements in Sw. Therefore, the utility function of the shipping companies in sailing route s can be obtained by:
$$U_s = (P-Q_s) \times \displaystyle{{D_s} \over {CC_s}},\quad \forall s\in \hbox{S}_w$$From Equations (1)–(7), the total cost function in sailing route s can be obtained as follows:
$$Q_s = \left( {Tw_s + \displaystyle{{SD_s} \over {V_s}}} \right)\left( {\displaystyle{{Pc_s} \over {T_a}} + Po_s} \right) + Pvc_s + \displaystyle{{F \times V_s^2 } \over {V^2}} \times SD_s \times Pf,\quad \forall s\in S_w$$ For any s ∈ S w, Pvc s > 0 and 
$0\le V_{s}\le V_{{s}_{max}}$ is obtained.
The number of shipping events in each route is determined by the total customer demand, which is related to the actual shipping cost of the customers in each route. Generally, the number of shipping events in a route is small if the actual shipping cost to the customers in this route is high, which can be expressed according to the following logit-type multi-path assignment model:
$$D_s = \displaystyle{{e^{-\bar{\theta }\lpar {\lambda P + T_s} \rpar }} \over {\sum\nolimits_{k\in S_w} {e^{-\bar{\theta }(\lambda P + T_{\rm s})}} }} \times D,\quad \forall \hbox{s}\in \hbox{S}_w$$Let:
$$\hbox{X}_{s}=e^{-\bar{\theta }\lpar \lambda P+{ T}_{s}\rpar },\quad \forall \hbox{s}\in \hbox{S}_{w}$$
$$\hbox{X}_{-s} = \sum\limits_{k\in S_w,{\kern 1pt} {\rm k}{\kern 1pt} \not{ = }{\kern 1pt} {\rm s}} {e^{-\bar{\theta }(\lambda P + T_s)},\quad \forall \hbox{s}\in \hbox{S}_w}$$Then, Equation (15) can be rewritten as:
$$V_s = \displaystyle{{\bar{\theta }SD_s} \over {\ln \lpar {\hbox{D}-\hbox{D}_s} \rpar -\ln \lpar {\hbox{D}_s} \rpar -\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P}},\quad \forall \hbox{s}\in \hbox{S}_w.$$From Equation (14):
$$\eqalign{Q_s = & a_s\lpar {\ln \lpar {\hbox{D}-\hbox{D}_s} \rpar -\ln \lpar {\hbox{D}_s} \rpar -\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P} \rpar \cr & + \displaystyle{{b_s} \over {{\lpar {\ln \lpar {\hbox{D}-\hbox{D}_s} \rpar -\ln \lpar {\hbox{D}_s} \rpar -\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P} \rpar }^2}} + e_s}$$ Thus, 
$a_s = \displaystyle{1 \over {\bar{\theta }}}\left({\displaystyle{{Pc_s} \over {T_a}} + Po_s} \right)$, 
$b_s = \bar{\theta }^2\displaystyle{F \over {V^2}}SD_s^3 Pf$, and 
$e_s = \left({\displaystyle{{Pc_s} \over {T_a}} + Po_s} \right)Tw_s + Pvc_s$.
Then, the utility function of the shipping companies in sailing route s can be rewritten as:
$$\eqalign{U_s &= \left( {P-Q_s} \right) \times \displaystyle{{D_s} \over {CC_s}} \cr & = \left( {{\rm P}-a_s(\left( {\ln \left( {{\rm D}-{\rm D}_s} \right)-\ln \left( {{\rm D}_s} \right)-\ln \left( {{\rm X}_{-s}} \right)-\bar{\theta }\lambda P)} \right)} \right. \cr & \, \,-\left. {\displaystyle{{b_s} \over {{\left( {\ln \left( {{\rm D}-{\rm D}_s} \right)-\ln \left( {{\rm D}_s} \right)-\ln \left( {{\rm X}_{-s}} \right)-\bar{\theta }\lambda P} \right)}^2}}-e_s} \right) \times \displaystyle{{D_s} \over {CC_s}}} $$It is assumed that each type of company can optimise its total travel by speed adjustment to maximise its profit. Then, the utility maximisation model in sailing route s can be expressed as:
$$max_{V_s,\ D_s}\lpar U_{s}\rpar =U_{s}(V_{s},D_{s},D_{-s})$$subject to:
$$D_{s}+D_{-s}=\hbox{D},\quad \forall \hbox{s}\in \hbox{S}_{w}$$
$$0\le D_{s}\le D,\quad \forall \hbox{s}\in \hbox{S}_{w}$$Let Ξs be a set that satisfies Equations (22b) and (22c). Then U s can be expressed as:
$$max_{Y_s\in \Xi_s}\lpar U_{s}(\hbox{Y}_{s},\hbox{Y}_{-s}) \rpar $$where U = (U s, U −s), Y s = (D s, V s), and Y −s = (D −s, V −s).
3.3. Concavity of the utility function
The upper-level model can be presented as a multi-player non-cooperative game model. The existence of the generalised Nash equilibrium solution is analysed using the mathematical properties of its utility function.
Proposition 2:
Utility function U s in Equation (22a) is strictly concave with respect to Ys.
Proof. From s ∈ S w:
$$\eqalign{ \displaystyle{{\partial \hbox{U}_s} \over {\partial \hbox{D}_s}} = &\left( {\displaystyle{{a_s} \over {CC_s}}-\displaystyle{{2b_s} \over {CC_s}}\displaystyle{1 \over {{(\ln \lpar {D-D_s} \rpar -lnD_s-\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P)}^3}}} \right)\left( {\displaystyle{D \over {D-D_s}}} \right) \cr & + \lpar {\hbox{P}-a_s\lpar {\lpar {\ln \lpar {\hbox{D}-\hbox{D}_s} \rpar -\ln \lpar {\hbox{D}_s} \rpar -\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P} \rpar } \rpar } \cr & -\left. {\displaystyle{{b_s} \over {{\lpar {\ln \lpar {\hbox{D}-\hbox{D}_s} \rpar -\ln \lpar {\hbox{D}_s} \rpar -\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P} \rpar }^2}}-e_s} \right)\displaystyle{1 \over {CC_s}}}$$
$$\eqalign{ \displaystyle{{\partial ^2\hbox{U}_s} \over {\partial \hbox{D}_s^2 }} =& \left( {\displaystyle{{a_s} \over {CC_s}}-\displaystyle{{2b_s} \over {CC_s}}\displaystyle{1 \over {{(\ln \lpar {D-D_s} \rpar -lnD_s-\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P)}^3}}} \right. \cr & -\left. {\displaystyle{{6b_s} \over {CC_s}}\displaystyle{1 \over {{(\ln \lpar {D-D_s} \rpar -lnD_s-\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P)}^6}}} \right) \times \left( {\displaystyle{{D^2} \over {{\lpar {D-D_s} \rpar }^2D_s}}} \right)}$$As mentioned, fuel consumption accounts for 50% of the total travel cost; thus:
$$\left( {\displaystyle{{Pc_s} \over {T_a}} + Po_s} \right)T_s < \displaystyle{{FV_s^2 } \over {V^2}}SD_sPf_S$$and:
$$\left( {\displaystyle{{a_s} \over {CC_s}}-\displaystyle{{2b_s} \over {CC_s}}\displaystyle{1 \over {{(\ln \lpar {D-D_s} \rpar -\ln D_s-\ln \lpar {\hbox{X}_{-s}} \rpar -\bar{\theta }\lambda P)}^3}}} \right) < 0$$From Equation (19):
$$\displaystyle{{\partial ^2\hbox{U}_s} \over {\partial \hbox{D}_s^2 }} < 0$$The utility function is proven to be strictly concave.
As utility function U s is strictly concave with respect to Ys, 
$\emptyset_{s}$ is composed of the linear equality in Equation (22b) and the inequality in Equation (22c). Then, according to the conclusion proven by Harker (Reference Harker1991), the QVI formulation of Equation (9) can be equivalent to the following VI problem: finding 
$\hbox{Y}^{\ast }\in \Xi $ such that:
$$F\lpar Y^{\ast}\rpar ^{T}(Y-Y^{\ast})\ge 0,\quad \forall Y\in \varXi_{s}(Y^{\ast })$$where Y = (Y 1, …, Y n) is the optimal solution of the GNEP in Equation (23).
4. SOLUTION ALGORITHM
To summarise, the whole bi-level model proposed in this paper is aimed at studying the game process between the customers and shipping companies whose profits are based on both total travel cost and total demand of customers, which is equivalent to solving the following optimisation problem:
$$max_{V_s,\ D_s}\lpar U_{s} \rpar =U_{s}(V_{s},D)$$ As mentioned earlier, solving Equation (27) is equivalent to finding 
$V^{\ast }=\lpar V_{s}^{\ast }\comma \; V_{-s}^{\ast }\rpar $ such that:
$${F\lpar D^{\ast },V^{\ast }\rpar }^{T}(V^{\ast }-V)\ge 0,\quad \forall Y\in \varXi_{s}(Y^{\ast })$$ With equilibrium point 
$V^{\ast }=\lpar V_{s}^{\ast }\comma \; V_{s}^{\ast }\rpar $, the total customer demand can be expressed as:
$$D^{\ast }=D(V^{\ast })$$In this paper the generalised reduced-gradient method has been used to solve the upper-level model, which is presented as a GNEP of Equation (31a). The iterative balancing method (Bell, Reference Bell1995) is used to solve the lower-level problem, which can be described as a logit-type multi-path assignment problem. The formulation of the algorithm to solve this bi-level model with speed optimisation is proposed as follows:
Step 1 Initialise the sailing speed on each route vj, j = 1
Step 2 Find vj [solve Equation (31a) using the generalised reduced-gradient method, see Figure 2 such that:
(32)Thus, the gradient of the utility function can be expressed as:$$\sum_k F(\hbox{v}_{s}^{j},\hbox{v}_{-s}^{j},\hbox{D}(\hbox{v}_{s}^{j},\hbox{v}_{-s}^{j})) \ge 0,\quad \forall \hbox{v}=(\hbox{v}_{s}^{j}\hbox{v}_{-s}^{j})\in \hbox{V}$$(33)Moving steps αj can be solved by:$$F = \displaystyle{{dU_s} \over {dV_s}} = \displaystyle{{\partial U_s} \over {\partial V_s}} + \displaystyle{{\partial U_s} \over {\partial D}}\displaystyle{{\partial D} \over {\partial V_s}}$$(34a)subject to:$$ {min}\, f(v^{j}+\alpha^{j}(v^{j}-z))$$(34b)The total customer demand in Equation (31b) can be solved using the iterative balancing method.$$0< \alpha ^{j}< 1,\quad \forall \hbox{j}\in J.$$Step 3 Output optimal speed v bes = ( v bes1, …, v besn) and optimal utility function U max = (U max 1, …, U max n).
Figure 2. Algorithm flow of the generalised reduced-gradient method.
5. CASE STUDY
5.1. Study area and work flow
In this paper, a case study is used to investigate the effectiveness of the bi-level model. The OD pair is from Dalian to Rotterdam, two routes are chosen to connect these two ports, one is the Northern Sea Route (NSR), while the other is the traditional route via Suez Canal (SUE) (Figure 3). For shipping companies, two kinds of companies can be selected for transport.
Company 1: Using the traditional route all the year round.
Company 2: Using the NSR when it is navigable and choosing the traditional route in the rest of the year
Figure 3. Study area.
Work flow of the case can be seen in Figure 4. Firstly, the bi-level model has been proposed to find the optimal speed and demand change of each company when the equilibrium between two companies is reached. After considering the harsh climatic and hydrological conditions in the Arctic, the vessel speed should be adjusted which leads to the adjustment of elastic demand. Finally, the economic potential of each company in one year can be analysed with the adjusted speed and demand.
Figure 4. Flow of case study.
5.2. Bi-level model for optimal speed and elastic demand
The case presented is based on a scenario originally used by Liu and Kronbak (Reference Liu and Kronbak2010) (with partial data from Way et al. (Reference Way, Khan and Veitch2015)) so as to demonstrate the application of the proposed model and to test the convergence of the solution algorithm. The simulated scenario is presented as follows:
The sailing distance on the NSR from Dalian to Rotterdam is 7,931 nautical miles, whereas that on the traditional route via the Suez Canal (SUE) is 10,907 nautical miles.
A 4,300 Twenty-foot Equivalent Unit (TEU) container ship has been chosen on the traditional route while a CSC3 ice class vessel on the NSR is assumed to have the same main dimensions except the lightship weight. The maximum speed of a ship on each route is no more than 25 knots.
Three fuel prices are assumed for comparison purposes: USD 350, USD 700, and USD 900 per ton.
The average pilotage and ice-breaking fee on the NSR is USD 446,000 per transit.
The voyage fee on the SUE is USD 240,800 per transit.
The price of a 4,300 TEU ship on the NSR is USD 5.28 M per year, whereas that of a vessel on the SUE route is USD 4.4 M per year.
The operating cost (including manning, hull and machinery insurance, protection and indemnity insurance, repairs and maintenance, administration, etc) on the NSR is USD 8,925 per day, whereas that on the SUE route is USD 6,100 per day.
The expected revenues of all shipping companies are assumed to be the same, that is USD 3 M per transit.
The total customer demand in the first equilibrium state is 4.3 M TEU per year.
The above scenario is introduced into the proposed bi-level model and the first equlibrium is considered, which is the competition result among shipping companies before the Arctic sea route rise, as Stage 1. Similarly, Stage 2 represents another equilibrium, which is the result of the competition among shipping companies with the rise of the Arctic sea route. The game begins at Stage 1. Each company type maximises its profits through speed adjustment, which can reduce its total travel cost in the upper-level and attract more customer demand in the lower-level. Finally, the new equilibrium in Stage 2 is reached. Throughout this game, three factors are used to compare Stage 1 with Stage 2: speed, maximum profits and total customer demand.
First, the selection of model parameters is analysed. The relationships between the sailing speed and dispersion parameter 
$\bar{\theta }$ are shown in Figure 5. It can be seen that under different fuel-price scenarios, the sailing speed increases, and the difference in the sailing speed between Stages 1 and 2 becomes more obvious as the value of 
$\bar{\theta }$ increases. This result can be interpreted to mean that a larger 
$\bar{\theta }$ implies a better understanding of the shipping companies that can more markedly change their sailing speed according to the variation in the customer demand. Table 1 lists the relationship between sensitivity parameter β and the total customer demand. It can be seen that the value of the total customer demand decreases as the value of β increases under different fuel-price scenarios both in Stages 1 and 2. Meanwhile, the difference in the total customer demand among different fuel-price scenarios becomes more definite with the increase in parameter β, which implies that the total customer demand is more sensitive to a larger value of β. In addition, if the value of β exceeds 0.05, the results in both Stages 1 and 2 greatly deviate from the original data (4.3 M TEU).
Figure 5. Dispersion parameter 
$\bar{\theta }$ and sailing speed.
Table 1. Sensitivity parameter and total customer demand. (The values in the brackets are the fuel-price scenarios: USD 350, USD 700, and USD 900 per ton).
In this numerical experiment, the dispersion parameter is given as 0.1, whereas the sensitivity parameter is given as 0.01. The convergence of the proposed solution algorithm is shown in Figures 6 and 7. It can be seen that stable solutions are obtained by the solution algorithm. Finally, the optimal speed and customer demand in each stage under different fuel-price scenarios are listed in Table 2.
Figure 6. Convergence of proposed model in Stage 1. (Red lines: optimal speed and demand of Company 1. Green lines: optimal speed and demand of Company 2. Black lines: total customer demand).
Figure 7. Convergence of the proposed model in Stage 2. (Red lines: optimal speed and demand of Company 1. Green lines: optimal speed and demand of Company 2. Black lines: total customer demand).
Table 2. Optimal speed and profits of each shipping company and total customer demand. (The values in the brackets are the fuel-price scenarios: USD 350, USD 700, and USD 900 per ton. Company 1 is defined as the shipping company on the SUE. Company 2 is defined as the shipping company on the NSR.)
By comparing the data (see Table 2) in Stages 1 and 2, some interesting phenomena are observed. First, shipping companies on different sea routes would take different measures to join the competition in the event the equilibrium in Stage 1 is broken with the rise of the NSR. The shipping companies on the traditional route would decide to minimise their delivery time cost by increasing speed, aiming at attracting more customer demand, which is their optimal strategy given the competition with other shipping companies on the NSR, although the profits themselves actually decrease in Stage 2. Meanwhile, for the NSR, shipping companies take the opposite measure, that is, lowering the shipping speed because the shorter shipping distance on the NSR guarantees that shipping companies can cut their travel cost by appropriately lowering their shipping speed as well as maintaining competitiveness, although the shipbuilding cost, daily operating cost, and even the voyage fee on the NSR are all higher than those for the SUE route. Therefore, the profits of the shipping companies in using the NSR increase in Stage 2. However, the equilibrium in Stage 2 is fragile. In fact, only in the case where no communication exists between these two types of shipping companies can this equilibrium be maintained. Otherwise, the SUE shipping companies will transform into NSR shipping companies when they realise that the profits in using the NSR are much higher than those in the SUE route, irrespective of the strategies they choose to maximise the profits on their current sea route. As a conclusion, a relative decline in the SUE may be a result of a rise of the use of the NSR.
Secondly, the total customer demand increases according to the competition among the companies that use the SUE and NSR routes. If these two company types are considered as a whole, it can be seen that the total profits rise in Stage 2 compared with those in Stage 1. Therefore, this Stackelberg game becomes a win–win game in which both customers and shipping companies optimise their operating patterns.
The conclusion that traditional routes may decline in face of the rise of the Arctic routes and the total volume of freight can be increased as the Arctic routes rise have been proved in previous works (Bekkers et al., Reference Bekkers, Francois and Rojas-Romagosa2016; Countryman et al., Reference Countryman, Francois and Rojas-Romagosa2016). Therefore, according to this numerical case, this model can well reflect the interaction mechanism between the Arctic and more traditional routes.
5.3. Economic potential based on climatic and hydrological conditions
The bi-level model is an ideal mathematic model, and it does not consider the impact of climatic and hydrological conditions. For the foreseeable future, vessel speeds will be affected by harsh Arctic weather (Nam et al., Reference Nam, Park, Lee, Mi, Choi and Seo2013). Data related to water depth comes from a product called ETOPO1 provided by the National Geophysical Data Center (NGDC), with a resolution of 1 arc-minute. Wave height data comes from the 1979–2017 daily global atmospheric reanalysis product ERA-Interim daily data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF), with an accuracy of 0.125° × 0.125°. Sea ice concentration and thickness data comes from the 2006-2100 daily data of the Community Climate System Model version 4 (CCSM4) model under the medium Representative Concentration Pathway (RCP45) provided by the National Oceanic and Atmospheric Administration/ National Centre for Atmospheric Research (NOAA/NCAR), with a lattice precision of 1° × 0.5°. Wind field data comes from the 2006–2010 monthly data from the GFDL-ESM2G model provided by NOAA under the RCP45, with an accuracy 1° × 0.43°. Current field data is derived from the 1993–2014 AVISO+ multi-satellite fusion daily data provided by the French National Space Research Center (CNRS), with an accuracy of 0.75° × 0.37°. Due to different spatial accuracy and time scales of the climatic and hydrological data, all the data should be interpolated into the same time and space scales. In this paper, time scale is assumed to be annual average, and the spatial scale is 1° × 0.5°.
5.3.1. Optimal speed adjustment
5.3.1.1. Navigable time on the NSR
The Arctic Ice Regime Shipping System (AIRSS) is used to evaluate the navigable days for a CSC3 ice class vessel on the NSR. In this system, the Ice Number (IN) index reflecting the navigable situation under specific ice conditions can be calculated as:
$$IN=(C_{a} IM_{a} )+(C_{b} IM_{b} )+\ldots+(C_{n} IM_{n})$$where C a represents the ice concentration of ice type a, IM a is the multiplier which is used to calculate the impact weight of ice type a on vessel navigation. When the IN is larger than zero, it means the CSC3 ice class vessel can safely sail in this area, otherwise, it cannot navigate there. The details of all kinds of ice and ice multiplier are presented in Tables 3 and 4. If the IN of all girds on the NSR are larger than zero, that day can be calculated as a navigable day.
Table 3. Ice type (Transport Canada, 1998).
Table 4. Ice multiplier (knots) (Transport Canada, 1998).
5.3.1.2. Ice condition to speed adjustment
Sea ice concentration and thickness are the key factors to Arctic vessel speed. Generally, the thicker and more concentrated the sea ice becomes, the greater the reduction of the ship's sailing speed. In some harsh ice condition areas, an icebreaker has to be hired to help navigate, which will undoubtedly increase the shipping cost for the whole route. The relationship between sea ice thickness and concentration can be seen in Table 5. DS is the designed vessel speed. In an area where sea ice concentration is less than 30%, without considering other meteorological and hydrological factors, the CSC3 vessel can sail at the DS speed. Additionally, the italic part refers to the vessel speed with the help of an icebreaker.
Table 5. Relationship between sea ice and vessel speed (Nam et al., Reference Nam, Park, Lee, Mi, Choi and Seo2013).
5.3.1.3. Wind, wave and current to speed adjustment
In open water areas, vessel speed is also affected by the wind, wave and current conditions. The impact of these factors on vessels depends on their value and direction and is presented in Table 6. For example, if the sailing direction is consistent with the direction of wind and flow, vessel speed will be increased, otherwise, it will hinder the navigation of vessels. In addition, when the sea ice concentration is larger than 30%, these factors can be neglected.
Table 6. Relationship between vessel speed and current, wave, wind (knots) (Nam et al., Reference Nam, Park, Lee, Mi, Choi and Seo2013).
Firstly, the navigable time has been calculated, and then the optimal speed derived from the bi-level model is incorporated into the above models and the results can be seen in Figure 8. The navigable days in the year 2050 can reach to 250 and the average adjusted speed 
$\bar{V}_{s}^{A}$ on the NSR is 10.54 (Oil price =350), 8.97 (Oil price =700), 8.41 (Oil price =900) while the speed on the traditional route is unchanged.
Figure 8. Adjustment of vessel speed on the NSR.
5.3.2. Elastic demand adjustment
According to Equations (11)–(13), the total demand can be changed according to the variation of time cost manipulated by the change of vessel speed, the equation can be rewritten as:
$$D^A = \hbox{D}^0\exp \left\{ {min\left[ {\displaystyle{\beta \over \theta }\ln \left( {\sum\limits_{s\in {\rm S}_w} {\exp } (-\bar{\theta }(\lambda P + T_s^A ))} \right)} \right]} \right\},\quad \forall s{\rm \epsilon }S_w$$where T A is the navigation time according to the adjusted vessel speed and can be presented as:
$$T_s^A = T_w + \sum\limits_{1:n} {\displaystyle{{SD_{sn}} \over {V_{sn}^A }}} ,\quad \forall s{\rm \epsilon }S_w$$When the total demand is decided, the elastic demand can be derived for each company according to Equation (16):
$$\hbox{D}_s^A = \displaystyle{{e^{-\bar{\theta }\lpar {\lambda P + T_s^A } \rpar }} \over {\sum\nolimits_{k\in S_w} {e^{-\bar{\theta }(\lambda P + T_s^A )}} }} \times D,\quad \forall \hbox{s}\in \hbox{S}_w.$$Based on Equations (1)–(7) and Equation (14), the annual profits of these two companies can be derived with adjusted speed and demand. Results are shown in Table 7. Ships on the NSR cannot navigate all year round even in the middle of this century and cannot reach the optimal speed due to the harsh weather. Using the bi-level model, the actual demand of Company 2 will be less, as will the annual profits, while the demand and annual profits of Company 1 will increase. Additionally, the actual total demand will have a tiny decrease. Overall, in this practical case, the NSR has potential economic competitiveness with the SUE route and will be a viable alternative route for container shipping.
Table 7. Demand and speed adjustment and annual profits. (The values in the brackets are the fuel-price scenarios: USD 350, USD 700, and USD 900 per ton. Company 1 is defined as the shipping company on SUE. Company 2 is defined as the shipping company on NSR).
6. CONCLUSION
Due to climate change, Arctic sea ice has been decreasing, which has started to open Arctic sea routes. Compared with traditional routes, Arctic routes are shorter, have fewer pirates, less congested waters and the potential for more sea room. However, the worse weather conditions, more fragile ecosystem, and lack of infrastructure, ports and emergency rescue make this comparison complex and multi-dimensional. The most direct method to compare the economic potential of these two routes is to calculate the average annual profits on each route by taking all the main factors into consideration and by assuming that all routes are independent of each other. In fact, an obvious interaction exists among the routes between any OD pair. Therefore, when the potential competitiveness of Arctic routes is considered, the impact of non-Arctic routes needs to be taken into account. Thus, a bi-level model from the perspective of game theory is proposed to solve this issue.
In this bi-level game model, two levels are presented to depict the dynamic variation of the world liner shipping market linking shipping companies (their strategies and profits) and customers (their cost and demand). The upper-level model demonstrates the competition among shipping companies who maximise their own profits through speed optimisation. The lower-level model illustrates the response of customers who reassign their demand with elastic total demand. A case study is presented to demonstrate the application of the proposed model and to test the convergence of the solution algorithm. The bi-level model that has been described in this paper is an ideal model that assumes the Arctic ice is completely melted and vessels can navigate on that route all the year round. To make the study more practical, the impact of harsh weather conditions have been examined in Arctic regions on vessel speed and elastic demand of each company in the case study section. With the adjusted vessel speed and demand, annual profits of each company are analysed. The results reveal that in 2050, the NSR has potential economic competitiveness compared to the traditional route and will be an alternative route for container shipping. Furthermore, this model can be applied to simulate the variation of the world shipping market demand and provide decision support for shipping companies to develop future development strategies (including route selection, new shipbuilding amount and type, sailing speed setting, etc).
ACKNOWLEDGMENT
This work has received funding from the National Natural Science Foundation of China under grant agreement number 51609254.
