Optimal ice routing of a ship with icebreaker assistance

Applied Ocean Research 86 (2019) 177–187

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Optimal ice routing of a ship with icebreaker assistance A.G. Topaj , O.V. Tarovik, A.A. Bakharev, A.A. Kondratenko ⁎

T

Krylov State Research Centre, 196158, 44, Moskovskoe shosse, Saint-Petersburg, Russia

ARTICLE INFO

ABSTRACT

Keywords: Ship routing Ice navigation Path finding Graph-based algorithms Wave-based methods Icebreaker assistance

Offshore development and growing prospects of commercial shipping in the Arctic pose the challenge of optimal ship routing in ice. Route selection in spatially distributed ice conditions significantly affects the voyage time and determines the efficiency of shipping. Most of the applied methods of ice routing solve the problem of a single vessel route selection without considering icebreaker support. At the same time, the real practice of ice navigation is closely connected with icebreaker assistance. It allows reducing the voyage time and fuel consumption, while having additional costs for icebreaker services. Such opposite trends set an optimization task that has not been studied in detail before. In this article, we presented the formulation of a Single Vessel and Icebreaker Assisted Ice Routing optimization problem in non-stationary ice conditions. We considered the icebreaker assistance as an integral part of the overall route optimization problem, and used the economic criterion to optimize both ship route and amount of icebreaker involvement. The article contains the adapted mathematical formulations of classical graph-based and wave-based routing problems in order to consider icebreaker assistance. To prove the practical applicability of these formulations, we developed special subject-oriented research software and implemented there both graph-based (Dijkstra, A*) and the wave-based ice routing methods. Using this developments, we conducted several case studies and made the analysis of strengthens and weaknesses of the alternative routing methods in case of ice operation. The results of the study may serve an additional step to the practical implementation of ice routing technologies and planning of icebreaker resources.

1. Introduction Navigation in ice-covered waters has a number of specific features in contrast to open water operation. It is especially relevant for the problem of route selection, since the ship path in ice significantly influences voyage time. Active development of ice routing methods began in the 2000s and was governed by the growth of Arctic navigation and new possibilities of radar satellite imagery [1]. Unlike weather-based routing in open water, ice routing is based on other principles and needs more data sources and factors to be considered for reliable decision-making [2]. The principal peculiarities of ice navigation are the following: a) instability of the environment in which ships move, both in temporal and spatial scales, b) critical influence of ice conditions on ship performance. The essence of the operational and tactical ice routing is the effective use of local areas with weak ice or ice lanes, even in case if it leads to a significant increase in sailing distance. Indeed, the real tracks of ice-going ships and icebreakers in ice are usually the complex ⁎

labyrinthine curves; especially for the case of high ice class ships (see www.marinetraffic.com). Recently, the problem of ship management in ice and, particularly, ice routing has attracted increased attention of researchers. The first attempts to develop ice routing tools were done in the late 1990s in the frame of ARCDEV project [3]. However, the majority of studies are dated by the middle 2000s. Kotovirta et al. [2] presented an algorithm of step-by-step “greedy” optimization with the use of Powell’s method. They set the task of ice routing in three-dimensional space (X–Y coordinates and time) taking into account a number of spatial restrictions, such as land area, pathways, and areas of shallow water. The algorithm interpolates the values of ice parameters between the nodes and represents vessel movement in a non-discrete manner with smooth routes. In the algorithm, travel time straightly determines the cost function due to the assumption that the ship moves at maximum speed. Choi et al. [4] contributed to the problem of ice routing by introducing the genetic optimization algorithm instead of the “greedy” one. They also analyzed the advantages of continuous routing in comparison with vessel movement through the fixed nodes, and presented a new formulation of the objective function as the weighted sum of distance, smoothness, time and clearness. The latter parameter represents ship safety in ice.

Corresponding author. E-mail addresses: [email protected], [email protected] (A.G. Topaj).

https://doi.org/10.1016/j.apor.2019.02.021 Received 15 November 2018; Received in revised form 24 January 2019; Accepted 26 February 2019 Available online 08 March 2019 0141-1187/ © 2019 Elsevier Ltd. All rights reserved.

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Nam et al. [5] set the objective function as the total cost defined as the sum of fuel price, operational and capital costs, port charges and icebreaker fee. They applied Dijkstra algorithm to choose the best way on a set of segments along the Northern sea route (NSR). In that study, ship engage an icebreaker for a fixed fee in any segment of the route in case if the predefined condition is met, i.e. optimization of the amount of icebreaker assistance is out of the investigation. Ice routing methods based on the A* heuristic algorithm have got much interest in recent time [6,7] due to the favorable calculation time. Some studies contain modifications of the basic A* algorithm for better applicability to the considered problem [7,8]. Implementation of the ensemble ice forecasts caused a series of studies intended to solve ice routing problem under uncertain or stochastic conditions [8–10]. Ensemble forecast differs from the traditional one by the fact that there is not a single predicted scenario, but a range of possible scenarios of future ice dynamics. In addition to the use of traditional optimization methods, recently there appeared several research studies devoted to the development of original ice routing approaches, such as application of the Voronoi diagram [11], the use of isochrones method [12] and the Finite Element Method-based potential theory [13]. Besides the pure scientific studies, nowadays there are several industrial projects in the field of ice routing. For example, IceTrails is the Internet service developed by Navidium company (Finland) in collaboration with Finnish Meteorological Institute, European Space Agency and other institutions. The main declared objective is the planning of ship navigation and optimization of ship routes in ice in the Baltic region (www.navidium.com). KARS is an Arctic safe routing system, developed by Korea Institute of Ships and Ocean Engineering. KARS main functions are the monitoring and prediction of weather and ice conditions, ship performance prediction using ship transit model and riskbased optimization of ship route [14]. The main area of KARS usage is the Europe-Asia transits along the Northern Sea Route (www.kriso.re. kr). The ice routing problem involves several interrelated issues, therefore any onboard or coastal system of tactical ice routing of ships should consist of several key components:

the time of convoy operations. The reason is that in many cases vessels can operate in ice either independently or under icebreaker assistance. In the latter case, ships increase the speed and save fuel, but also have extra costs for icebreaker services. Obviously, there is some economically proved level of icebreaker assistance with the minimum total costs of passage through ice-covered waters. Thus, the general mathematical formulation of routing problem could be done as follows: T

S, B

opt

= arg min S, B

C (S (t ), B (t )) dt

(1)

0

where S(t ) = X (t ); Y (t ) is the vector of ship (or convoy) positions on the X–Y plane of geographical coordinates at times t, i.e. temporal dynamics of S(t) determines the route of the ship; B(t) – boolean process, indicating that the ship is assisted by icebreaker at time t; C – total cost per time unit. Optimization problem (1) has a set of constraints. Boundary conditions may be written as: (2)

S (0) = S0 , S (T ) = ST

where S0 and ST are the fixed initial and destination points of the route. The right time boundary T can be either free or fixed value. In the latter case, the one will need to put an additional constraint on the given date of voyage end. Operational limitations on the ship speed in specific ice conditions may be formulated as: (3)

|| S (t ) || = V (K (S, t ), B (t ))

where V is the vessel or convoy speed, calculated from the ship transit model for a given vector K of environmental parameters (ice, weather, sea depth, etc.) at times t and in geographical locations S. Speed V may be either the attainable (maximum) or operational speed of the ship or a convoy. We suppose to know both the ice parameters at any time t in any geographical point S and the ship transit model, i.e. the internal logic of V (K, B) algorithm. In formulation (3) we made a simplifying assumption that ice environment is isotropic and speed V does not depend on the direction of movement, while in the most general case of ice routing it is not so because of ice lanes and cracks in ice cover. Integrand in (1) is the law of total unit cost calculation; it forms the optimization criterion. In the simplest case C 1 and the problem (1) reduces to finding the fastest way and minimizing the total voyage time. The solution obtained in this case will be similar to the one that corresponds to the use of the economic criterion of minimum freight costs for the independent movement of the transport ship, i.e. C cS , where cS is a freight rate of the ship. In cases that are more sophisticated, the value C should take into account additional components of operational costs, such as fuel consumption, icebreaker assistance, etc. In general, C can be expressed as:

1) geographical component includes the tools for monitoring of the current ice state and the models for ice dynamics forecasting at different time horizons, 2) technological component relates to the ship transit modeling in case of independent and icebreaker-assisted operation, 3) mathematical component is the methods and algorithms for optimal path finding by chosen criteria and under a number of constraints. In this article, we address the issues related to the third aspect, i.e optimal ice routing of ships and convoys in arbitrary non-stationary ice environment using mathematical programming methods. 2. Statement of ice routing problem under icebreaker assistance

j csupp (S , t )

C (S (t )) = cS + cIB B (t ) + cfuel f (K (S , t )) +

Three basic elements define any optimization problem; they are the control variables, the optimization criterion, and the constraints. The ice route optimization implies the possibility of varying three variables: ship speed, the trajectory of movement (route) and the use of icebreaker. Traditionally, speed and route are considered as the control variables in ship weather routing in open water conditions [15]. Somewhere only the trajectory varies for a predetermined speed [12], somewhere ‒ only the speed for a given route [16], and somewhere ‒ both characteristics together [17]. At the same time, ice navigation is characterized by the increased risks of ship hull and propeller damage or unplanned delays during the voyage. Therefore, vessels always tend to pass ice areas minimizing sailing time and having the highest level of safety. Since icebreakers serve these purposes, the new object for optimization appears in this task. Apart from the optimum route in ice, it is also required to find the optimum amount of icebreaker assistance or

j

(4)

where cS is a freight rate of the ship, cIB – freight rate for the selected icebreaker, cfuel – fuel price, f – total fuel consumption per time unit of the ship or convoy in the specified environmental conditions, j csupp – any supplementary costs, which can be both spatially and temporary dependent. For example, additional costs may include a one-time fee for involving an icebreaker. Thus, any formal structure of the optimization criterion that contains time, fuel and involved resources can be reduced to a versatile monetary expression (4). Another important feature of the practical ice navigation is the spatial environment in which ships move, that is, the type of K(S) 178

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

function. In general, this two-dimensional relationship cannot be considered as spatially continuous. This is not only because of the presence of a coastline, but also due to the need to take into account a large number of fairways or predefined routes – the narrow segments in which the possibility of movement and the speed differ significantly from the surrounding area. Good examples are the navigation channels in shallow waters or renewable ice channels in fast ice. The task stated by Eqs. (1)–(4) can be formally classified as a calculus of variations problem. Its analytical or even numerical solution in a complete formulation is impossible in real practice. However, there are various approaches that allow making spatial and/or temporal discretization of an originally continuous problem as well as the different methods to solve the obtained discrete task. This article presents the author’s experience in the use and comparative analysis of two such classical approaches to the optimal ship and convoy routing in non-stationary ice conditions. They are the graph-based and wave-based methods. The principal feature that distinguish our approaches from the similar studies published earlier by Kotovirta et al. [2], Choi et al. [10], Nam et al. [5] and other researchers are the fairway accounting and the use of a universal economic criterion to optimize both ship route and the amount of icebreaker assistance. At the same time, in our formulation we do not consider the stochastic factors and related risks of ice navigation, since these issues turned out to be not of the primary interest for navigators and ship owners.

|E| – edge length, VS (E, t) – the attainable speed of ship independent navigation on the edge E at time t. The fact that the edge price depends on the time causes serious difficulties when using classical methods for path optimization on weighted oriented graphs, such as dynamic programming [19] or Dijkstra algorithm [20]. These approaches are based on a successive marking of vertices by tentative distance values. For so-called backward algorithms (dynamic programming), this operation is done in reverse order "from the end to the beginning", while the forward algorithms, like Dijkstra’s one, make it in a direct way. In both cases, the edge weight must be a constant value, which is not time-dependent and known in advance. The latter is impossible for unsteady conditions, when the environmental parameters, ship speed, transit times and associated costs depend on when the ship reaches a particular point, i.e. from the entire prehistory of its voyage. The widely used approach to overcome this restriction is the use of 3D-graph instead of 2D representation [21], when the third dimension is the discrete time. This common approach allows taking into account any type of environmental nonstationarity in the frame of conservative methods of path optimization on graphs, but at the same time, it causes an increase in dimension and computational complexity of the problem. Another way to cope with time-dependence of the edge price is the use of various heuristic methods of search by the first best match on the graph, e.g. the widely used A* algorithm. Although this algorithm is approximate and does not guarantee to find the best route, it can be easily adapted to solve the problem of path optimization in deterministic time-dependent dynamic networks [22]. The icebreaker-oriented modification of the path optimization problem on an oriented graph is illustrated in Fig. 2 as a two-layer network [23]. The upper layer describes an independent ship operation and corresponding edge costs C determined by Eq. (5). The lower layer represents the alternative edge price for the case of icebreaker-assisted operation and considers the convoy speed, as well as the convoy operational costs. Fig. 2 shows the simplest case of calculation of operational costs, which are estimated based on the ship and icebreaker freight rates only, with no fuel consumption or additional costs. In this case, the cost of the lower layer edge could be expressed as:

3. Adaptation of ice routing methods to account for icebreaker assistance 3.1. Graph-based approach The essence of a classical graph-based (or cell-based) method of path finding is the replacement of a continuous space with a predefined set of discrete points that form a network with a number of edges. There are many ways to build such a network [18]. The simplest method of spatial discretization is the use of a regular orthogonal two-dimensional grid, which is oriented along the orthodromic line and has a constant cell size in the latitudinal and longitudinal directions. In general, the less the cell size, the higher the solution accuracy (due to the higher spatial resolution of the network) and, on the other hand, the higher the computational complexity of the problem. To transform a grid to a graph, the network vertices should be connected by the edges. Each point can be connected with its nearest vertical and horizontal neighbors only (the rectangular graph), although the setting of diagonal edges between vertices of several distant layers is also possible (the graph with a higher coherence). In order to take into account the fairways, some modifications of any chosen method of graph construction should be undertaken. In addition to the basic points in a regular grid, the special points for every fairway need to be constructed. These points lie exactly on the fairway and connected by the edges (possibly curved) that exactly follow the path. The distance between the graph vertices on a fairway should be chosen in a way to correspond to the resolution of the main grid. Internal fairway points are connected with each other only, while the input and output points are stitched with the nearest adjacent points of the main grid. Fig. 1 demonstrates the technique of such graph construction. The weight of every edge (i.e. the total unit cost C) may be calculated by means of direct simulation of vessel movement along the corresponding route segment. Ship speed and all operational costs are determined by the ice- and weather parameters at the corresponding geographic point and at the particular time. Generally, the total unit cost CE of some edge E can be expressed as: j csupp (E ,

CE = C (S (t )) tE = cS + cfuel f (E , t ) + j

|E| t) VS (E , t )

CE = (cS + cIB )

|E| VSIB (E , t )

(6)

where VSIB(E, t) – the attainable speed of the convoy (ship + icebreaker) on the edge E at time t. The transition from the upper layer to the lower one is possible at any point and transition cost CE is a product of two parameters:

CE = cIB TIB

(7)

where TIB – the conditional time that icebreaker needs to reach the point where it meets the ship. The cost of involving an icebreaker cIB · TIB approximately consider the fact that icebreaking resources are limited and ship cannot get an icebreaker instantly and free of charge at the required point, as well as the ship is not able to involve an icebreaker as many times as necessary. The price of the reverse transition to the upper layer from the lower one is assumed to be zero, i.e. the release of an icebreaker is free of charge. By means of changing TIB value, the one can manage the degree of segmentation of route parts, assisted by an icebreaker. The more TIB is, the less the number of icebreaker-assisted sections on the route. Since each case of icebreaker request is rather expensive, its maximum use for each attraction will be cost-effective. 3.2. Wave-based approach The family of the cell-free numerical methods for path finding combines different approaches that do not require prior explicit spatial discretization. In particular, below we consider the wave-based

(5)

where ΔtE is the time of passing the edge E, 179

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Fig. 1. An example of graph grid construction. The path (bold line) connects to the fairway (dashed line).

methods, i.e. the group of algorithms where the solution of routing problem is based on the consistent construction of equal level lines in geographical space. A transparent physical analogy to this approach is the wave propagation in the inhomogeneous medium. We may note here the well-known teleonomic interpretation of the law of geometric optics – the Fermat’s principle. It states that the path taken between two points by a ray of light is the path that can be traveled in the shortest time. An analogy with waves on the liquid surface may be even more evident. Let us imagine that we initiated a surface wave in the nonviscous fluid at the initial point of the route S0 (for example, by dropping a heavy stone). We consider the speed of wave propagation at a specific spatial location as an analogy of the attainable speed of the model ship. The same as ship speed is determined by ice severity, the wave speed depends on fluid density at the current point. Eventually, the wavefront in such nonviscous fluid will sooner or later reach the destination point. The path, along which the first wave reached the end point, can be interpreted as the optimal route. Perhaps, the most logical way of mathematical modeling the wave

propagation is the use of Huygens-Fresnel principle, which states that every point on a wavefront is itself a source of a spherical wave; the sum of these secondary waves determines the form of the wave at any subsequent time (see Fig. 3). The first, classical and more natural interpretation of the wavebased approach assumes that the objective function is the total voyage time. In this case, isoline is the envelop curve of distances that the ship can reach for the same voyage time, so the approach is called the isochron method [24]. Formalization of the problem in its general meaning allows extending the mentioned interpretation to any optimization criteria. Thus, in the case of minimizing the fuel consumption, isochrons are converted into isopones [25], and in the case of selecting the purely economic criterion – into the isocosts, i.e. the lines that the ship can reach with equal cost. When doing spatial discretization, the wavefront lines are replaced by the sets of points that form the discretized equal-level lines of the optimization criterion C. In practical implementation, each step of numerical wave propagation consists of two stages: a) generation of all

Fig. 2. Graph-based formulation of ice routing problem considering icebreaker assistance. 180

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

areas of prohibited navigation (shallow water, coast), can be surely excluded from further consideration. Then, the dimension of a wavefront should be reduced by eliminating all internal points, since they obviously less promising. There are several techniques proposed to solve the task of point reduction in wave-based numerical methods. The first one is a sub-channel-based method introduced in [26] and later extended and improved by several researchers [27]. We will describe the principal idea of this method according to [12]. N parallel lines of equal spacing ΔD are arranged on both sides of the global direction vector (e.g., orthodromic line) to form 2 N stripes or subchannels (see Fig. 5). The approach stipulates that ship can sail only within the area with the width of 2 N × ΔD. At each discrete time step (i.e. for the certain iteration of wave propagation), only one “best” point remains on each stripe and forms the next isocost. The best point is the one that is closest in distance to the end point. Such a method of point number reduction leads to the so-called "greediness" of the algorithm. One of the promising alternative methods to reduce the number of points in a wavefront is a Concave Hull computation technology [28]. The idea of the algorithm is to preserve only the external points that form the Pareto subset on the current wavefront after propagation stage. It could be done by constructing a convex or non-convex curve, which envelopes a set of points and covers the area occupied by them. For the case of non-convex shape, the problem of constructing an envelope curve does not have a single solution and the one can obtain several polygon boundaries by varying the internal parameter of the method. We also need to note here the original approach to reduce the number of points and model wave propagation by means of geo-algorithms [29]. According to this approach, secondary wavefront generated by each point is not a set of points, but a continuous circle of a given radius. Using standard geo-techniques, the algorithm merges all newly plotted circles and forms a new polygonal object, which boundary (polyline) represents a new isochrone. Then the spatial polynomial filter, based on the well-known Savitzky-Golay algorithm, smoothes the boundary polyline of the obtained object. After that, the algorithm determines a set of new reference points on a new polyline. The distance between reference points is a variable that is equal to the shortest distance between the current point and previous isochrone. The latter leads to a higher spatial resolution of the isochrone shape in the areas with slow movement (e.g., in severe ice conditions). The main disadvantage of this method is a big amount of complicated and longtime calculations with the use of geo-algorithms. The final stage of the algorithm in Fig. 4 is a check whether the propagating wave reached the end point. This could be done in a simple way. If the distance from progenitor point to some heritor point, generated on the next step, is less than the distance to the end point, then all heritor points are substituted by the destination point. After that, the resulting optimal route is built step by step using the reverse motion principle, i.e. starting from the destination and consistently passing from each point to its progenitor. The improvement of the isocost method in order to take into account icebreaker services requires its significant modification. Firstly, let us define several variables in addition to the ones defined in Eqs. (4)–(6): ΔQ – step of isocost propagation in monetary terms; CS – unit cost of vessel operation at the current time and location, determined by Eq. (4) presuming cIB = 0; CIB – unit cost of icebreaker operation at the current time and location, determined by Eq. (4) presuming cS = 0 and B (t ) 1 (i.e. CIB does not consider the cost of involving an icebreaker, which is cIB · TIB) Each point on the isocost or wavefront is described by an extended set of parameters, including geo-coordinates, progenitor point, time to reach it from the start point and the status. The latter parameter is a specific integer value that characterizes the type and genesis of the point. The Semantics of possible status values is as follow:

Fig. 3. The layout of wave propagation (subsequent wavefronts) in the problem of path optimization.

points in a new wavefront and b) reduction of number of points by excluding “unpromising” ones. Fig. 4 describes the wave-based algorithm flowchart for the case of independent vessel movement. The initial set of wavefront points contains the start point only. Then, the wavefront expands from step to step according to the following rules. Each point of the previous set generates the “fan” of points of the next level set (see Fig. 3). The angle between neighboring rays determines the desired spatial resolution. The distance Ri of a single step i, that is the radius of propagation of the secondary wave from each point S, is calculated as:

Ri = Vi (K (S, t )) ti (S , t ) = Vi (K (S, t ))

C Ci (S (t ))

(8)

where ΔC is the selected constant resolution of contour lines in terms of the criterion. In case of finding the fastest route, Eq. (8) degenerates to Vi (K (S, ti )) t that is the classical problem of isochrone method when t is a constant predetermined value. Generation of multiple descendants takes place only for the parent points located in free water areas. In order to consider fairways and features of ice navigation, some modifications of the common fan-based numerical algorithm should be done. The first one is that point lying on the fairway produces only two new points on the fairway curve (linear propagation forward and backward). In turn, if the fan of rays, emitted by a parent point, crosses the fairway line, then the corresponding point should be assigned to the intersection. The second modification is a permission to a virtual ship to stand in a place during the limited time in case if the ship has no ability to move at the current moment due to the local ice obstructions (e.g., ice compressions). When such an event happens, the current point is not removed from the next isocost, because a rapid change in ice environment may soon resuscitate this path for further wave propagation. Finally, it is important to note that each point of the new wavefront contains information about its progenitor. The necessity of the reduction stage is due to the fundamental problem of any wave-based method. It is the overcoming the "curse of dimension", i.e. the choice of a suitable algorithm to reduce the number of points on the next level line and construct the envelope curve. Indeed, if we leave all the generated points on the next isocost, their total number will grow exponentially and the corresponding algorithm will be computationally impracticable. The reasonable idea here is to leave only the most promising points on the current wavefront and exclude all others. Firstly, all new points, which turned out to be in 181

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Fig. 4. The flowchart of wave-based routing algorithm for the case of ship independent navigation.

3.2.1. The preliminary calculation of several supplementary characteristics NIB is the “initial delay of an icebreaker” – the integer value, which determines the number of isocost steps ΔQ to compensate the cost of icebreaker involvement, i.e.

NIB = 1 +

cIB TIB Q

(9)

QIB is the “remained charge” – the float value that defines the amount of conditional money remaining to ship movement after paying the last portion of icebreaker involvement cost, i.e.

QIB = NIB Q

Fig. 5. The sub-channel-based method to reduce the point number in a wavefront.

cIB TIB = Q 1

cIB TIB Q

(10)

In Eqs. (9) and (10), the square and curly brackets are the operators of taking the integer and fractional parts of the expression respectively.

Status -1 means that vessel moves from the point independently, without icebreaker assistance; Status 0 means that vessel moves from the point accompanied by an available icebreaker, i.e. without the need to pay for icebreaker involvement; Status greater than or equal to 1 means that the current point is fictitious and represents the rest of the payment to involve an icebreaker. The positive value of such status determines the number of elementary trenches in ΔQ portions, which is left to pay at the current point to get an icebreaker. The proposed algorithm of wave-based isocost propagation with accounting for the ability of icebreaker assistance consists of four main steps:

3.2.2. Generation of the first isocost At the start location of the route the algorithm generates two points that form the wavefront of initial isocost. Status of the first point is set to -1 (independent navigation), the status of the second one is NIB (such point corresponds to the case, when an icebreaker is requested right at the start of the route). 3.2.3. Run of the recursive algorithm of wave propagation The algorithm consequently generates the next isocost from progenitor points, until fulfilling the condition of reaching the end point or 182

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Fig. 6. Rules of isocost propagation from the points with different status values.

facing the situation when there are no points on the next isocost (complete stop due to heavy ice). The technique of point generation at each step is closely connected with the status of progenitor point.

algorithm for the reduction of a number of points in a wavefront is applied separately for the sets of points with the equal status values. So, the resulting set of points that form the next isocost is the assembly of reduced sets of points with status values varied from -1 to NIB.

a For the point with status -1 (independent navigation) a fan of isocost propagation paths is generated according to the common principles of the wave-based algorithm (see Fig. 6A). The distance of propagation R-1 is calculated by:

R

1

=

VS Q CS

3.2.4. Checking the end point of the route and building the resulting optimal route This is done in the same way as for the basic algorithm. The only difference is that statuses of the points of a received route additionally determine the segments where icebreaker escort is profitable.

(11)

At each destination point of the obtained fan, two sub-points are generated: one with the status -1 (continue of movement without icebreaker) and the other with status NIB (request to involve an icebreaker). The current point becomes the progenitor for both newly generated points (the same is true for cases in Fig. 6B and D). b For the point with status 0 (convoy movement) a fan of isocost propagation is also created (see Fig. 6B), but the distance R0 considers the convoy speed:

R0 =

VSIB Q CS + CIB

4. Test cases In order to demonstrate and compare different numerical methods of ship and convoy routing in ice, we developed special subject-oriented research software. It allows us to preprocess and perform calculations as well as visualize the results of running various path-finding algorithms. The application encloses the built-in GIS system and spatially referenced database of weather and ice conditions. The software has a modular structure and contains several components from the earlier developed integrated software solution for simulation and multidisciplinary analysis of marine transportation systems [23]. They are:

(12)

The same as the in a previous case, at each point of the fan two subpoints are generated having status 0 (continue of convoy movement with an icebreaker) and -1 (icebreaker release and the start of independent navigation). We do not consider the case of releasing an icebreaker and ordering it again immediately (i.e. generation the point with status NIB from the point with status 0), because such strategy looks evidently unprofitable. c The point with status K, where 1 < K ≤ NIB, represents the internal stages of sequential payment of the “credit” for icebreaker involvement. It generates only one point in the next isocost with a status decremented by one (see Fig. 6C). The newly generated point has the same coordinates as the previous one, i.e. next isocost contains the same geographical points, but their statuses are decreased by one. d When the point status is 1 (the last stage of credit payment for involving an icebreaker), the credit is fully paid and virtual vessel gets an icebreaker, which was ordered a number of steps ago. Thus, the rest of the cost according to (9) can be spent on the movement of the new-formed convoy. Therefore, a standard fan of isocost propagation paths could be created (see Fig. 6D); the distance of propagation R1 is calculated by:

R1 =

VSIB QIB CS + CIB

• Ship transit model “Mechanic” – the external library to calculate the

(13)

•

The same as before, two points are generated at each destination point of the fan: the one with status 0 (continue of convoy movement) and the other with status -1 (icebreaker release). The last strategy can take place, for example, when a vessel needs an icebreaker to pass only through a narrow strip of hummocks. e When point generation is done, the algorithm eliminates the internal points and generates the front line of the next isocost. The chosen

• 183

attainable ship speed and fuel consumption in various environmental conditions. General functional of this module covers resistance calculations, analysis of propeller regime and shaft power predictions for any feasible combinations of the following cases: calm water, wind and waves of arbitrary direction (considering mandatory speed decrease due to slamming); restricted sea depth; independent operation in ice (including hull ice strength constraints); operation under icebreaker assistance and in freezing channels. Open water calculations are done using the existing semiempirical formulas to predict ship resistance, while the model of ship operation in ice is based on the empirical-statistical model described in [30] and [31]. In this model, the attainable ship speed in ice is calculated by the regression equations derived from the fullscale on-board observations in the period of 1970s. We made several improvements of this model in the following areas: accounting for double-acting ships (stern- and bow-forward motion regimes); reconstructing the values of ice resistance using the predicted ship speed and net trust curve; accounting for high tonnage ice-going ships. However, the detailed description of the ship transit model is out of the scope of this article, since it does not influence the mathematical formulation of vessel and convoy routing problem. Expandable and editable repository of informational models of different transport ships and icebreakers. Each model contains a structured description of parameters that define propulsion performance of a particular ship: main dimensions, loading cases, characteristics of propeller groups, power plant parameters and others. Databases of ice conditions in the Russian Arctic for the light, medium, heavy and extreme scenarios of ice severity. Each scenario

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

• •

represents a benchmark year-round spatial picture of the temporal dynamics of main ice parameters, such as concentration, age, shape, ridging, snow cover, decay rate, compression, etc. Ice parameters are presented in a grid form; time step is a month; spatial resolution is 0.25°. All ice data was provided by the Russian Arctic Antarctic Research Institute (AARI) based on the archive ice charts from the 1960–2014 period. Geo-layers that represent in a simplified form the navigational restrictions due to ship draught. To consider such restrictions we use smoothed 20 m isobatic line. Plugins of both graph- and wave-based computational methods for searching the optimal ship routes in ice. Each method has its own setting parameters that define spatial resolution, form or optimization criterion and details of the calculation process.

because for an arbitrarily small cell size the total “Manhattan distance” remains constant. There are many approaches intended to overcome such restriction of cell-based methods. The basic idea of the method of “smooth routes” between discrete nodes [2] was successfully implemented by Zhang et al. [7] for improved A* algorithm. According to this approach, not the center of the current grid cell serves as a graph node, but the whole edge of each cell is considered as an infinite set of nodes. Unfortunately, all such modifications are very computationally expensive. In turn, the wave-based method has its own limitations and drawbacks (see Fig. 8B). Many problems are caused by the chosen algorithm to reduce the number of the wavefront points. For example, the evident bottleneck of the sub-channel-based method is the “greediness” of this algorithm. The typical example is represented in Fig. 8B, where the modified isochrone algorithm fails to find an optimal route in specially designed non-uniform conditions. We set difficult ice conditions on the great circle route and form a trap, i.e. the dead-end channel of open water area. The greedy approach for selecting the best points in current isochrone leads to the route through the dead-end open water corridor and then through the heavy ice. Such route takes much more time in comparison with the bypass route in open water provided by Dijkstra algorithm. More advanced approaches to remove surplus points from the wavefront, e.g. Convex Hull algorithms, may overcome this drawback. At the same time, their use requires serious computational resources that make it difficult to apply such algorithms in practical operational planning of vessel navigation with high spatial and temporal resolution. Finally, we would note here that mentioned weak points of various routing methods are due to the basic logic of corresponding numerical algorithms. At the same time, the proposed modifications to consider the possibility of icebreaker support in the cell-based and cell-free approaches do not introduce any additional computational or methodological troubles.

We made several calculations to demonstrate the results of various routing approaches. Fig. 7 and Table 1 represent optimal ship routes for the selected pair of source and destination points at the Northern Sea Route for the month of January. In order to visualize the multi-parameter ice conditions in some graphical form, we used the attainable speed of a powerful nuclear icebreaker LK-60 as a single ice severity index. This index is displayed in Fig. 7 by a color: the lighter the color of a grid cell the less the speed of icebreaker and, therefore, the more severe the ice. All presented results are obtained for the ice-going tanker “Kirill Lavrov” (project P-70046) and diesel-electric icebreaker LK-25 (project 22,600) under the assumption that ship and convoy move with the maximum (attainable) speed. In addition, we assume that the freight rates of tanker and icebreaker are equal (cS = cIB) and the conditional overpayment time for the icebreaker (TIB) is 2 days. Fig. 7A and 7B show the optimal routes for independent (when icebreaker is totally unavailable) and icebreaker-assisted voyages obtained by A* graph-based and isocost wave-based methods for the light type of ice conditions. One can see that the obtained paths demonstrate the efficiency of using the local areas of weak ice (large polynya to the west of Novaya Zemlya), even despite the increase in sailing distance. In case if the icebreaker is available, the shape of optimal routes varies insignificantly, but there appear several icebreaker-assisted route segments (bold line) in the areas where icebreaking support is cost-effective. As it could be seen, these segments correspond to the heavy ice regions. The routes in Fig. 7A and 7B differ in their trajectories, but in a given pattern of spatial ice distribution, all of them pass through the Cape Zhelaniya. Calculations for another environmental pattern that describes heavy ice conditions (see Fig. 7C) lead to fundamentally different geography of optimal routes. In this case, the optimal way passes through the Kara gate. Moreover, icebreaker escorting turns out to be cost effective in two sections of the route (close to the end point and when going around Yamal), despite the one-time costs for icebreaker involvement.

6. Conclusion Icebreaker support is a cornerstone of shipping in heavy ice conditions. The ability to seize an additional icebreaking resource in order to facilitate passage through the heavy ice is an important factor that needs to be taken into account when planning the routes of transport ships in the Arctic. In this paper, we proposed modifications of well-known methods of mathematical programming to find the optimal path in ice covered waters in order to take into account the possible involvement of an icebreaker. So far, the aspect of icebreaker assistance has been considered in ice routing problem just as an additional resource involved when the predetermined conditions are fulfilled. We considered the icebreaker assistance as an integral part of the overall formulation of the route optimization problem in the frame of graph-based (Dijkstra, A*) and wave-based numerical routing algorithms. The universal economic criterion forms an optimized objective function that includes both timedependent expenses and one-time payments. This allows us to take into account the fact of a limited number of the available icebreakers. At the same time, this approach has several uncertainties and the main one is the time TIB that icebreaker needs to reach the point of forming a caravan. Obviously, if there are several icebreakers and a number of cargo ships in the region, TIB value will differ for each spatial point and depend on the plan of icebreaker assistance of the whole fleet. Therefore, an accurate definition of TIB value is possible when considering a marine transport system as a whole. Based on these considerations, we assume that the algorithm of ice routing with icebreaker can serve as the component of complex computer-aided systems for planning and managing of icebreaker and transport fleets in the Arctic. Thus, there are several promising directions for the future development of the proposed approach:

5. Discussion The strengths and weakness of the Dijkstra and A* graph-based methods, as well as the isocost wave-based method, are demonstrated in Fig. 8 on the specially designed counterexamples. In all cases, the problem of finding the fastest path is under solution, i.e. the optimization criterion is the total voyage time. The possible limitations of graph-based methods are shown in Fig. 8A. The A* algorithm turned out to be unable to find any path, which allows minimizing the distance of passage through the fast ice. Apparently, this is due to the selected value of the “search width” heuristic parameter. In turn, Dijkstra algorithm can find a more profitable route than a great circle arc, but such a route is not optimal, because it consists of several almost straight sections, which can be smoothed. Such result is caused by the main methodological restriction of any graph-based method – an ability to move only along the edges of a predefined grid. This raises a known problem of the “Manhattan distance” that cannot be resolved by the traditional mesh grinding,

• application to the real problems of icebreaker support planning, taking into account the current location of icebreakers and the competition of cargo ships for limited icebreaking resources;

184

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Fig. 7. Optimal routes determined by various methods for the cases of independent navigation (red line) and icebreaker assisted voyages (black line, icebreakerescorted segments are bold): A – A* graph-based method and ice conditions of a light type; B – isocost wave-based method and light ice conditions; C – isocost wavebased method and heavy ice conditions (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

185

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Table 1 The results of route optimization for test cases. No

Icebreaker type

Light ice conditions 1 – 2 22,600 3 – 4 22,600 Heavy ice conditions 5 – 6 22,600

Calc. method

Calc. time (sec)

Total cost (106 RUR)

Total time (hours)

Time in convoy (hours)

Overall distance (nm)

Distance in convoy (nm)

А* А* Isocost Isocost

1.42 1.71 4.11 13.49

32.7 30.9 31.9 29.2

157 126 153 100

0 19.8 0 23.3

1027 1004 935 937

0 140 0 173

Isocost Isocost

4.32 15.00

37.1 32.9

178 137

0 17.3

828 912

0 102

• operational •

management of caravan movement considering one icebreaker and several transport vessels, i.e. planning the points of forming a caravan and its release; estimation of the marginal cost of ice forecasts in terms of routing efficiency, i.e. estimation of how the accuracy and completeness of

•

ice forecasts will influence shipping and what is the reasonable price of the services for ice dynamics prediction; development and testing of the new methods to reduce the number of points in a wavefron of the wave-based methods.

Fig. 8. Counterexamples that demonstrate the advantages of wave-based (A) and graph-based methods (B). 186

Applied Ocean Research 86 (2019) 177–187

A.G. Topaj, et al.

Acknowledgements

[15] L. Walther, A. Rizvanolli, M. Wendebourg, C. Jahn, Modeling and optimization algorithms in ship weather routing, e-Navi 4 (2016) 031–045, https://doi.org/10. 1016/j.enavi.2016.06.004. [16] M. Sotnikova, E. Veremey, Algorithms for motion optimization on a given trajectory taking into account weather forecast and constraints, Proc. Of 17th IFAC Workshop on Control Applications of Optimization 17 (2018) 25b. [17] Y.-H. Lin, M.-C. Fang, R. Yeung, The optimization of ship weather-routing algorithm based on the composite influence of multi-dynamic elements, Appl. Ocean Res. 43 (2013) 184–194, https://doi.org/10.1016/j.apor.2014.12.005. [18] R.H. Motte, S. Calvert, On the selection of discrete grid systems for on-board microbased weather routing, J. Navigation 43 (1990) 104–117, https://doi.org/10.1017/ S0373463300013849. [19] C. de Wit, Proposal for low cost ocean weather routing, J. Navigation 03 (1990) 428–439, https://doi.org/10.1017/S0373463300014053. [20] X. Zhu, H. Wang, Z. Shen, H. Lv, Ship weather routing based on modified Dijkstra algorithm, Proc. of MMEBC 2016 6 (2016) 696–699, https://doi.org/10.2991/ mmebc-16.2016.147. [21] E.I. Veremei, M.V. Sotnikova, Optimal routing based on weather forecast, INJOIT 4 (2016) 55–61. [22] I. Chabini, S. Lan, Adaptations of the A* algorithm for the computation of fastest paths in deterministic discrete-time dynamic networks, IEEE Trans. Intell. Transp. Syst. 3 (2002) 60–74. [23] O.V. Tarovik, A. Topaj, A.A. Bakharev, A.V. Kosorotov, A.B. Krestyantsev, A.A. Kondratenko, Multidisciplinary approach to design and analysis of arctic marine transport systems, Proc. of OMAE-2017 8 (2017) 10, https://doi.org/10. 1115/OMAE2017-61951 OMAE2017-61951. [24] R.W. James, Application of wave forecasts to marine navigation, Navy Hydrographic Office, Washington DC, 1957. [25] M.B. Klompstra, Time Aspects in Games and in Optimal Control, Ph.D. Thesis, Delft University of technology, Delft, 1992. [26] H. Hagiwara, Weather Routing of (Sail-Assisted) motor Vessels, Ph.D. Thesis, Delft University of technology, Delft, 1989. [27] J. Szłapczyńska, R. Śmierzchalski, Adopted isochrones method improving ship safety in weather routing with evolutionary approach, Reliabil. Risk Anal. Theory Appl. 2 (2008) 139–145. [28] A. Moreira, M.Y. Santos, Concave hull: a k-nearest neighbours approach for the computation of the region occupied by a set of points, Proc. of the Second International Conference on Computer Graphics Theory and Applications Volume GM/R, 61-68 (2007) 8–11. [29] R. May, A. Rubchenia, V. Smolyanitsky, O. Tarovik, Typification of sea ice conditions in the arctic based on voyage times of ships on optimized routes, Proc. Of 18th SGEM-2018 (2018) 513–528, https://doi.org/10.5593/sgem2018/2.2/S08.065. [30] N.M. Adamovich, A.Ya. Buzuyev, V.E. Fedyakov, The empiric model of vessel movement in the ice and generalization of the experience of the model usage in hydrometeorological support of shipping in the arctic, Proc. of POAC-1995 2 (1995) 30–40. [31] A.I. Brovin, S.V. Klyachkin, S.U. Bhat, Application of an empirical-statistical model of ship motion in ice to new types of icebreakers and ships, Proc. of OMAE-1997 IV (1997) 43–49.

The research is carried out at the expense of the Russian Science Foundation (Project No. 17-79-20162). References [1] V.G. Smirnov, Ye.U. Mironov, Ice information support to marine transport system using modern technologies for delivery of hydrometeorological information to the end user, Proceedings of RAO/CIS Offshore-2009 2 (2009) 149–153 (In Russian). [2] V. Kotovirta, R. Jalonen, L. Axell, K. Riska, R. Berglund, A system for route optimization in ice-covered waters, Cold Reg. Sci. Technol. 55 (2009) 52–62, https:// doi.org/10.1016/j.coldregions.2008.07.003. [3] Final Public Report of the ARCDEV Project, (1998) https://trimis.ec.europa.eu/ sites/default/files/project/documents/arcdev.pdf. [4] M. Choi, H. Yamaguchi, L.W.A. De Silva, Application of genetic algorithm to ship route optimization in ice navigation, Proc. Of the 22nd POAC, (2013). [5] J.H. Nam, I. Park, H.J. Lee, M.O. Kwon, K. Choi, Y.K. Seo, Simulation of optimal arctic routes using a numerical sea ice model based on an ice-coupled ocean circulation method, Int. J. Naval Archit. Ocean Eng. 5 (2013) 210–226, https://doi. org/10.2478/IJNAOE-2013-0128. [6] R.E. Guinness, J. Saarimäki, L. Ruotsalainen, H. Kuusniemi, F. Goerlandt, J. Montewka, R. Berlund, V. Kotovirta, A method for ice-aware maritime route optimization, Proc. of PLANS 2014 (2014), https://doi.org/10.1109/plans.2014. 6851512. [7] M. Zhang, D. Zhang, S. Fu, X. Yan, C. Zhang, A method for planning arctic Sea routes under multi-constraint conditions, Proc. Of the 24 Nd POAC. POAC, (2017), pp. 17–042. [8] Y. Wang, R. Zhang, L. Qian, An improved A* algorithm based on hesitant fuzzy sets theory for multi-criteria arctic route planning, Symmetry 10 (2018) 765, https:// doi.org/10.3390/sym10120765. [9] P. Schutz, Dynamic routing through waters partially covered with Sea ice, Proc. Of the OTC. OTC-24658 (2014), https://doi.org/10.4043/24658-MS. [10] M. Choi, H. Chunga, H. Yamaguchi, K. Nagakawa, Arctic sea route path planning based on an uncertain ice prediction model, Cold Reg. Sci. Technol. 109 (2015) 61–69, https://doi.org/10.1016/j.coldregions.2014.10.001. [11] X. Liu, S. Sattar, S. Li, Towards an automatic ice navigation support system in the arctic Sea, ISPRS Int. J. Geoinf. 5 (2016) 36, https://doi.org/10.3390/ijgi5030036. [12] H. Wang, P. Li, Yu. Xue, M.V. Korovkin, Application of Improved Isochron Method in Ship’s Minimum Voyage Time Weather Routing 13 Vestnik Sankt-Peterburgskogo Universiteta, Prikladnaya Matematika, Informatika, Protsessy Upravleniya, 2017, pp. 286–299, https://doi.org/10.21638/11701/spbu10.2017.306 (In Russian). [13] H. Piehl, A.S. Milakovic, S. Ehlers, A finite element method-based potential theory approach for optimal ice routing, J. Offshore Mech. Arct. 139 (2017) 061502, https://doi.org/10.1115/1.4037141. [14] S.Y. Jeong, K.J. Kang, H.-S. Kim, J.-J. Kim, M.-I. Roh, A Study of Ship Voyage Planning in the Northern Sea Route, Proc. of the 13th Pacific-Asia Offshore Mechanics Symposium, (2019), pp. 579–584.

187