Optimal 3D trajectory planning for AUVs using ocean general circulation models

Ocean Engineering 188 (2019) 106266

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Optimal 3D trajectory planning for AUVs using ocean general circulation models Sultan Albarakati a ,∗, Ricardo M. Lima a , Loïc Giraldi a ,1 , Ibrahim Hoteit b , Omar Knio a a b

Computer, Electrical, and Mathematical Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia Physical Science and Engineering Division, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia

ARTICLE

INFO

Keywords: AUV Trajectory planning Optimization Ocean general circulation model Red Sea

ABSTRACT In this paper, we consider the autonomous underwater vehicle (AUV) trajectory planning problem under the influence of a realistic 3D current as simulated by an ocean general circulation model (OGCM). Attention is focused on the case of a deterministic steady OGCM field, which is used to specify data for both the ocean current and for ocean bathymetry. A general framework for optimal trajectory planning is developed for this setting, accounting for the 3D ocean current and for static obstacle avoidance constraints. A nonlinear programming approach is used for this purpose, which leads to a low complexity discrete-time model that can be efficiently solved. To demonstrate the efficiency of the model, we consider the optimal time trajectory planning of an AUV operating in the Red Sea and Gulf of Aden, with velocity, and bathymetric data provided by an eddy-resolving MITgcm. Different optimal-time trajectory planning scenarios are implemented to demonstrate the capabilities of the model to identify trajectories that adapt to favorable and adverse currents and to avoid obstacles corresponding to a complex bathymetry environment. The simulations are also used to evaluate the performance of the proposed approach, and to illustrate the application of advanced visualization tools to interpret the model predictions.

1. Introduction Autonomous underwater vehicles (AUVs) are self controlled vehicles that can navigate in under water environments. They are currently utilized for the purpose of collecting data on coastal ecosystems, search and rescue operations, oil and gas exploration, inspection of underwater cables and pipelines, operation of intake and discharge of thermal plants, management of shipping operations, as well as characterization of the ocean state (Blidberg et al., 1991; Aghababa, 2012; Bayat et al., 2017). AUVs are also being used in military applications, including surveillance of installations and sea vessels (Soulignac et al., 2008; Lolla et al., 2012; Yuh, 2000). With recent technological advances, AUVs are becoming increasingly economical to build and operate, and their range is being expanded to deeper seas (Wynn et al., 2014). Whereas AUVs aim to operate in a self-contained and intelligent fashion, it is still essential to provide trajectory planning information so that the specified mission can be safely and effectively achieved. Besides, optimal path planning solutions play an important role in increasing the range of operation of the AUV for a fixed amount of available energy (Zeng et al., 2015). A case-in-point concerns the fact

that the strength of the ocean current can be equal or larger than the vehicle relative velocity (Schmidt et al., 1996; Elisseeff et al., 1999), so it is essential to consider trajectory planning when designing specific missions. Various objectives have been considered, for instance based on minimizing the risk of encountering hazardous areas (Warren, 1990; Dijkstra, 1959; Ferguson and Stentz , 2006; Garau et al., 2005; Lolla et al., 2012; Petres et al., 2005; Soulignac et al., 2008), energy consumption (Subramani et al., 2016; Alvarez et al., 2004; Garau et al., 2009; Koay and Chitre, 2013; Kruger et al., 2007), total travel time (Petres et al., 2007; Thompson et al., 2010; Lolla et al., 2014), or maximizing the quality of data collected by the AUV (Yilmaz et al., 2008; Hollinger et al., 2012). See Yang et al. (2016) for an overview. Most of the planning methods for AUVs consider 2D environments (Wang et al., 2016; Petres et al., 2007; Sun et al., 2010; Yilmaz et al., 2008; Wang et al., 2019; Warren, 1990; Soulignac, 2011; Garau et al., 2005; Lolla et al., 2012), which limit the trajectory of the vehicle to a constant depth. Therefore, 2D methods are not adequate for complex 3D environments, for example, for navigation along irregular seabeds, nor for complex flow fields having substantial variation with depth.

∗ Corresponding author. E-mail addresses: [email protected] (S. Albarakati), [email protected] (R.M. Lima), [email protected] (L. Giraldi), [email protected] (I. Hoteit), [email protected] (O. Knio). 1 Present address: Technocentre Renault, 1 avenue du Golf, 78084 Guyancourt, France.

https://doi.org/10.1016/j.oceaneng.2019.106266 Received 28 January 2019; Received in revised form 4 June 2019; Accepted 4 August 2019 Available online 26 August 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

point of view, most works in the literature have relied on one or more of the following idealizations: (1) the models and case studies consider only static synthetic obstacles (Alvarez et al., 2004; Kruger et al., 2007; Ataei and Yousefi-Koma, 2015; Liu et al., 2015; Ma et al., 2019) or dynamic synthetic obstacles (Witt and Dunbabin, 2008; Aghababa, 2012; Zeng et al., 2014; Li et al., 2019); (2) synthetic 3D current fields (Alvarez et al., 2004; Kruger et al., 2007; Witt and Dunbabin, 2008; Aghababa, 2012; Zeng et al., 2014; Li et al., 2019) or 2D current fields (Smith et al., 2010; Zeng et al., 2015) are involved in the models; or (3) 3D currents are not considered, just a 3D space with obstacles (Ataei and Yousefi-Koma, 2015; Liu et al., 2015; Ma et al., 2019). This reveals a gap in addressing 3D trajectory planning problems involving realistic obstacles and ocean currents.

Nomenclature Sets, Indices  

indices for the time steps indices for time horizon

Parameters 𝜒 𝛥𝑡 𝐱(𝑡𝑓 ) 𝐱(𝑡0 ) 𝐱(𝑡𝑤 ) 𝐱𝑓 𝐱0 𝑎max𝑧 𝑎max 𝑁 𝑡0 𝑡𝑖 𝑣max𝑧 𝑣max 𝑥𝑓max 𝑥𝑓min 𝑥max 𝑥min 𝑦𝑓max 𝑦𝑓min 𝑦max 𝑦min 𝑧𝑓max 𝑧𝑓min 𝑧max 𝑧min

the spatial domain time step (constant) final position vector initial position vector coordinate of the way point (𝑥𝑤 , 𝑦𝑤 , 𝑧𝑤 ) destination starting point upper bound on the acceleration of the vehicle relative to the current in the 𝑧 direction upper bound on the acceleration of the vehicle relative to the current in 𝑥𝑦 plane number of grid points in the time discretization initial time grid point time grid point upper bound on the velocity of the vehicle relative to the current in the 𝑧 direction upper bound on the velocity of the vehicle relative to the current in 𝑥𝑦 plane maximum 𝑥 coordinate of the destination target minimum 𝑥 coordinate of the destination target maximum 𝑥 coordinate of 𝜒 minimum 𝑥 coordinate of 𝜒 maximum 𝑦 coordinate of the destination target minimum 𝑦 coordinate of the destination target maximum 𝑦 coordinate of 𝜒 minimum 𝑦 coordinate of 𝜒 maximum 𝑧 coordinate of the destination target minimum 𝑧 coordinate of the destination target maximum 𝑧 coordinate of 𝜒 minimum 𝑧 coordinate of 𝜒

From the solution methodology perspective, trajectory planning methods can be divided into four categories (Wang et al., 2019): (1) sampling-based methods; (2) grid-based methods; (3) random searchbased methods; (4) mathematical programming-based methods. Rapidly-exploring Random Trees (RRT) (Kuffner and LaValle, 2000) is a popular sampling-based method that randomly builds a spacefilling tree to represent the collision-free areas. The method is easy to implement, fast to get feasible solutions and to sample high dimensional spaces, but it may lead to suboptimal solutions. Fast marching methods (Petres et al., 2007), level set methods (Sethian, 1999; Lolla et al., 2014), and the A* method (Hart et al., 1968) are extensively used as trajectory planning methods. These methods work by finding collisionfree paths over a discrete version of the space and can be adjusted to account for current fields. One of the drawbacks of these methods resides in the lack of scalability for high-dimensionality problems. Yilmaz et al. (2008) used a mixed integer linear programming model to plan the trajectory of an AUV with the objective of maximizing a learning index from the environment. In their model, the vehicle is limited to a discrete set of directions from one period to the next one. Recently, Wang et al. (2019) proposed a detailed mixed integer nonlinear programming model and solution approach for 2D trajectory planning for problems with nonlinear flow fields. There is also a considerable amount of work on testing and comparing a plethora of metaheuristics based on various types of phenomena (evolution, jumps of frogs, refraction of light, flow of water, an orchestra playing, flower pollination, behavior of swarms, bats, bees, ants, birds, flies, wolves, etc.), based on simulated annealing, and tabu search (e.g. Witt and Dunbabin, 2008; Besada-Portas et al., 2013; Alvarez et al., 2004). See Li et al. (2019) for a discussion on some of these algorithms applied to trajectory planning for AUVs, Sörensen (2015) for a general discussion on metaheuristics, and Hooker (1995) on considerations on testing and comparing algorithms. In this work, we use an optimization framework that relies on solving a nonlinear programming (NLP) problem with an interior point NLP solver: IPOPT (Wächter and Biegler, 2006). Selection of this approach is motivated by the mathematical foundations of NLP solvers, the improved performance of these solvers, and their potential to unlock capabilities afforded by the underlying optimization framework for trajectory planning purposes. NLP solvers are gradient- and Hessianbased solvers that seek solutions that respect the KKT conditions, thus guaranteeing local optimality for the corresponding solution. The main advantage of state-of-the-art NLP solvers is that they can handle high-dimensional search spaces in contrast to stochastic based search algorithms. Note that despite the fact that IPOPT does not guarantee convergence to global optima, previous studies in trajectory planning have found that it leads to near optimal solutions (Borrelli et al., 2006). We thus adopt an optimization framework built on the top of a stateof-the-art NLP solver, IPOPT, and implement strategies to overcome some of the limitations of this solver. Specifically, we leverage insights

Variables 𝐮(𝐱) 𝐚(𝑡) 𝐯(𝑡) 𝑔(𝐱)

𝑡𝑓 𝑡𝑤

non-homogeneous steady ocean current field (𝑢𝑥 (𝐱), 𝑢𝑦 (𝐱), 𝑢𝑧 (𝐱)) acceleration of the vehicle with respect to the current (𝑎𝑥 (𝑡), 𝑎𝑦 (𝑡), 𝑎𝑧 (𝑡)) instantaneous velocity of the AUV with respect to the current (𝑣𝑥 (𝑡), 𝑣𝑦 (𝑡), 𝑣𝑧 (𝑡)) smoothed indicator function, which defines whether the vehicle is inside or outside the obstacle time to get to destination time to get to the way point

Three-dimensional trajectory planning methods are more appropriate in these situations, at the expense of more elaborate models, and concomitantly additional algorithmic difficulties. From the modeling 2

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

from the specific trajectory planning problem studied to improve the performance of NLP solvers. Furthermore, as we discuss below, we make a relevant contribution to the incorporation of OGCM data into an NLP problem formulation. To our knowledge this is the first work that addresses establishment of such capability. The goal of this study is to develop a 3D trajectory planning framework that can account for multiple static obstacles and the ocean state, as provided by an Ocean General Circulation Model (OGCM) simulation or forecast. Attention is focused on the minimal travel time problem, and a general formulation is developed that accounts for vehicle kinematic constraints. Following the 2D approach of Wang et al. (2019), we adopt a nonlinear programming (NLP) approach that incorporates a diversified initialization/parallel solution strategy that aims to avoid local minima and efficiently provide solutions close to the global optima. Here, we develop alternative methodologies that enable us to (1) effectively accommodate general obstacle avoidance constraints, and (2) specify both the ocean current and the obstacle constraints directly based on OGCM data. To illustrate the application of the resulting framework, we consider a case study of minimal time trajectory planning problems for AUVs operating in the Red Sea and Gulf of Aden, and rely on MITgcm model simulations to describe the ocean environment. This paper is organized as follows. Section 2, defines the minimum time trajectory planning problem for an AUV operating in a steady, deterministic ocean current. The NLP approach to this problem is discussed in Section 3. Solution techniques developed to solve the NLP problem are then outlined in Section 4. Section 5 presents the case study results, which are also used to assess the capabilities of the proposed framework. Main conclusions are discussed in Section 6.

1 𝑑𝑡 = 𝑡𝑓 − 𝑡0 .

In this section, we describe the discretization scheme for the continuous trajectory planning problem defined above. We also outline the discrete version of the objective function and the target region, and the methodology accounting for OGCM data in the resulting NLP model. 3.1. Time discretization The problem formulation in Section 2 provides the continuous objective function and constraints for solving the AUV trajectory planning problem. This section formulates a discrete time version of the problem and uses it to solve the AUV trajectory planning problem. Let  = {1, 2, 3, … , 𝑁} be a set of indices, and  = {𝑡0 , 𝑡1 , 𝑡2 , … , 𝑡𝑁 } be a collection of discrete times. We use 𝑡0 = 0 as initial time, and assume that the travel time 𝑡𝑓 ∈  . A standard approach in trajectory planning models (Schouwenaars et al., 2001; Subramani and Lermusiaux, 2016; Wang et al., 2016, 2019) is to use a fixed time step, in which the time grid is divided into steps of equal size, 𝛥𝑡. We thus have

(1)

subject to: 𝑑𝐱(𝑡) 𝑑𝑡

(2)

= 𝐯(𝑡) + 𝐮(𝐱(𝑡)),

(12)

3. NLP formulation

𝑡𝑓

∫𝑡0

𝐯(𝑡), 𝐚(𝑡) ∈ R3 ,

In order for the trajectory to avoid all obstacles, 𝑔(𝐱(𝑡)) is required to be positive at all times, 𝑡0 ≤ 𝑡 ≤ 𝑡𝑓 . The constraint 𝐱(𝑡) ∈ 𝜒 forces the vehicle to navigate within the problem domain.

We tackle the specific problem of 3D underwater trajectory planning. The objective is to minimize the total travel time of a single AUV from starting point to destination, in the presence of current, and subject to obstacle avoidance constraints as well as kinematic constraints that limit the magnitude of the acceleration and velocity vectors. We focus on the case of steady deterministic current field in a fixed environment, assumed known throughout the AUV mission, and that the position and size of the obstacles are stationary and given. Accordingly, the trajectory planning problem is expressed as

𝐱(𝑡),𝐯(𝑡),𝐚(𝑡),𝑡𝑓

(11)

where 𝐱(𝑡) = (𝑥(𝑡), 𝑦(𝑡), 𝑧(𝑡)) denotes the instantaneous position of the AUV at time 𝑡. The vehicle is assumed to travel in a non-homogeneous, steady current field 𝐮(𝐱(𝑡)) = (𝑢𝑥 (𝐱(𝑡)), 𝑢𝑦 (𝐱(𝑡)), 𝑢𝑧 (𝐱(𝑡))). The variable 𝐯(𝑡) = (𝑣𝑥 (𝑡), 𝑣𝑦 (𝑡), 𝑣𝑧 (𝑡)) is the instantaneous velocity of the AUV with respect to the current, and the variable 𝐚(𝑡) = (𝑎𝑥 (𝑡), 𝑎𝑦 (𝑡), 𝑎𝑧 (𝑡)) is the acceleration of the vehicle with respect to the current. (Note that the total, ground-referenced, AUV velocity is 𝐮+𝐯.) The starting and ending time of the AUV trajectory are denoted by 𝑡0 and 𝑡𝑓 , respectively. The corresponding starting point and ending point are 𝐱(𝑡0 ) = (𝑥0 , 𝑦0 , 𝑧0 ), 𝐱(𝑡𝑓 ) = (𝑥𝑓 , 𝑦𝑓 , 𝑧𝑓 ) respectively. The vehicle’s velocity and acceleration relative to the current in the 𝑥, 𝑦 plane are respectively bounded by 𝑣max and 𝑎max ; the maximum velocity and acceleration of the AUV relative to the current in the 𝑧 direction are respectively bounded by 𝑣max𝑧 and 𝑎max𝑧 . The smoothed indicator function 𝑔 is defined on the domain 𝜒, so that: { 𝑔(𝐱) > 0 𝐱 is in the fluid, (13) 𝑔(𝐱) ≤ 0 otherwise.

2. Problem statement

min

𝐱(𝑡) ∈ 𝜒 ⊂ R3 ,

𝛥𝑡 = 𝑡1 − 𝑡0 = 𝑡2 − 𝑡1 = ⋯ = 𝑡𝑁 − 𝑡𝑁−1 , 𝑑 2 𝐱(𝑡) 𝑑𝑡2

= 𝐚(𝑡) + ∇𝐮(𝐱(𝑡)) ⋅

𝐱(𝑡0 ) = 𝐱0 ,

𝑑𝐱(𝑡) , 𝑑𝑡

(3)

𝑡𝑖 = 𝑖 ⋅ 𝛥𝑡,

(4)

The discrete version of the kinematic equations is derived based on a second-order Taylor expansion at 𝑡𝑖 . Specifically, we use:

𝑖 ∈ .

𝐱(𝑡𝑖+1 ) ≈ 𝐱(𝑡𝑖 ) + (𝐯(𝑡𝑖 ) + 𝐮(𝐱(𝑡𝑖 ))) ⋅ 𝛥𝑡 + 𝐱(𝑡𝑓 ) = 𝐱𝑓 , √ √

(14)

1 𝐚(𝑡 ) ⋅ 𝛥𝑡2 , 2 𝑖

𝑖 ∈ ,

(15)

(5)

𝑣2𝑥 (𝑡) + 𝑣2𝑦 (𝑡) ≤ 𝑣max ,

(6)

𝑎2𝑥 (𝑡) + 𝑎2𝑦 (𝑡) ≤ 𝑎max ,

(7)

𝐯(𝑡𝑖+1 ) ≈ 𝐯(𝑡𝑖 ) + 𝐚(𝑡𝑖 ) ⋅ 𝛥𝑡,

𝑖 ∈ ,

(16)

‖𝑣𝑧 (𝑡)‖ ≤ 𝑣max𝑧

(8)

We consider two approaches for specifying the AUV destination. In the first approach, the target destination is specified as a single point, (𝑥𝑑 , 𝑦𝑑 , 𝑧𝑑 ). In this case, the constraint of reaching the target point is imposed:

‖𝑎𝑧 (𝑡)‖ ≤ 𝑎max𝑧 ,

(9)

𝐱(𝑡𝑓 ) = 𝐱𝑑 .

𝑔(𝐱(𝑡)) > 0,

(17)

In the second approach, a relaxed version of the target constraint is considered. For this purpose, we define a region 𝛺(𝑥𝑑 , 𝑦𝑑 , 𝑧𝑑 ) around

(10) 3

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 1. Schematic illustration of data flow, showing the use of current and obstacle functions to communicate OGCM data to the NLP model.

the target point (𝑥𝑑 , 𝑦𝑑 , 𝑧𝑑 ). Whereas the region 𝛺(.) can have a general shape, in this work we restrict our attention to cubical regions defined by the two corners (𝑥𝑓𝑚𝑖𝑛 , 𝑦𝑓𝑚𝑖𝑛 , 𝑧𝑓𝑚𝑖𝑛 ) and (𝑥𝑓𝑚𝑎𝑥 , 𝑦𝑓𝑚𝑎𝑥 , 𝑧𝑓𝑚𝑎𝑥 ) (Wang et al., 2019). To ensure that the vehicle reaches the target region at the final time, 𝑡𝑓 , we impose the following inequalities: 𝑥𝑓𝑚𝑖𝑛 ≤ 𝑥(𝑡𝑓 ) ≤ 𝑥𝑓𝑚𝑎𝑥 𝑦𝑓𝑚𝑖𝑛 ≤ 𝑦(𝑡𝑓 ) ≤ 𝑦𝑓𝑚𝑎𝑥 𝑧𝑓𝑚𝑖𝑛

𝑓

≤ 𝑧(𝑡 ) ≤

(18)

𝑧𝑓𝑚𝑎𝑥 .

3.2. Evaluation of steady current field and obstacle avoidance constraints We demonstrate the method within a steady-state framework using a single snapshot extracted from an assimilative ocean general circulation model (OGCM) simulation of the Red Sea. A detailed description of the modeling framework used to generate the underlying data set is given in Toye et al. (2017). A brief outline of this framework, largely adapted from Toye et al. (2017) and El Mohtar et al. (2018) is provided in the case study, see Section 5.1. The trajectory planning model used is based on two variables extracted from the OGCM current and obstacles data present in the domain topography. We extract these variables from the OGCM output files as schematically illustrated in Fig. 1. Because the OGCM data is defined on a relatively coarse grid, we interpolate the OGCM nodal values to obtain a smooth representation of the continuous problem, using a tricubic interpolation in three dimensions (Lekien and Marsden, 2005). Specifically, the interpolation is applied in each cell of the grid using the nodal values associated with the corresponding corner points, see Fig. 2; this can be expressed as 𝑓 (𝑥, 𝑦, 𝑧) =

3 ∑

𝑎𝑖,𝑗,𝑘 𝑥𝑖 𝑦𝑗 𝑧𝑘 ,

Fig. 2. Schematic for tricubic interpolation. 𝑝1 , 𝑝2 , . . . , 𝑝8 are the eight corners of the cube, and 𝑝 is the point inside the cube bounded by points 𝑝1 , 𝑝2 , … , 𝑝8 .

first associating +1 value for nodes in the fluid, and -1 for the remaining nodes. Tricubic interpolation is then applied to the resulting discrete indicator field, leading to a sufficiently smooth indicator function 𝑔 that satisfies the formulation in Eq. (13). A sample of discrete indicator field is depicted in Fig. 3, and contrasted with the corresponding interpolated field. 3.3. The NLP problem To sum up, the trajectory planning problem we solve is modeled as the following NLP problem with the objective function: min

𝐱(𝑡𝑖 ),𝐯(𝑡𝑖 ),𝐚(𝑡𝑖 ),𝑡𝑓

(19)

𝑖,𝑗,𝑘=0

𝑡𝑓 ,

(20)

subject to (a) the discrete kinematic equations

where 𝑓 (𝑥, 𝑦, 𝑧) is the value of the interpolated function inside the cubic cell with corners denoted as 𝑝1 , 𝑝2 , . . . , 𝑝8 , and the values of 𝑓 are specified together with its first and second derivatives. See Appendices A and B for full details. The information extracted from the OGCM data file consist of nodal values of velocity, as well as discrete (binary) values that indicate whether or not a node lies in the sea. The velocity nodal values are interpolated directly as outlined above. This interpolation leads to a representation of the ocean current field that is sufficiently smooth, and that can be readily exploited to determine continuous gradient and Hessian fields. This is a key feature of the approach because the applied NLP solvers require the gradient and Hessian of the objective function and of the constraints. For the obstacle constraints, we proceed in a similar fashion as for the ocean current field. Specifically, starting from the information extracted from the topography data, we define a discrete indicator field by

1 𝐱(𝑡𝑖+1 ) = 𝐱(𝑡𝑖 ) + (𝐯(𝑡𝑖 ) + 𝐮(𝐱(𝑡𝑖 ))) ⋅ 𝛥𝑡 + 𝐚(𝑡𝑖 ) ⋅ 𝛥𝑡2 , 2 𝐯(𝑡𝑖+1 ) = 𝐯(𝑡𝑖 ) + 𝐚(𝑡𝑖 ) ⋅ 𝛥𝑡,

𝑖 ∈ ,

𝑖 ∈ ,

∀𝑡𝑖 ∈  , (21)

∀𝑡𝑖 ∈  , (22)

(b) the discrete kinematic constraints √ 𝑣2𝑥 (𝑡𝑖 ) + 𝑣2𝑦 (𝑡𝑖 ) ≤ 𝑣max , ∀𝑡𝑖 ∈  , √ 𝑎2𝑥 (𝑡𝑖 ) + 𝑎2𝑦 (𝑡𝑖 ) ≤ 𝑎max , ∀𝑡𝑖 ∈  ,

(23)

‖𝑣𝑧 (𝑡𝑖 )‖ ≤ 𝑣max𝑧 ,

∀𝑡𝑖 ∈  ,

(25)

‖𝑎𝑧 (𝑡𝑖 )‖ ≤ 𝑎max𝑧 ,

∀𝑡𝑖 ∈  ,

(26)

(24)

for 𝑖 ∈ , (c) the obstacle avoidance constraint, 𝑔(𝐱(𝑡𝑖 )) > 0, 4

∀𝑡𝑖 ∈  ,

(27)

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 3. Level set contours showing the obstacle locations in a plane of fixed depth; left: contours based on raw OGCM data, right: contours based on interpolated level set function. The raw OGCM data is mapped so that a level of +1 corresponds to a location in the fluid, whereas a level of −1 corresponds to a location within an obstacle. The smooth level set field identifies fluid and solid regions by level sets larger or equal to 1 and smaller or equal to −1, respectively. Note the smooth but sharp transition of the level set function near the boundaries in the left plot.

(d) the starting point at time 𝑡0 ,

4.1. Procedure to solve the full and reduced NLP problems (28)

𝐱(𝑡0 ) = 𝐱0 ,

This sub-section describes a procedure to handle the reduced and full NLP problems with unknown final time. By construction 𝑡𝑓 ∈  , which means that the final time is one of the discrete time points postulated in the time discretization. Therefore, we propose a solution procedure where a set of final times are enumerated, and for each final time an NLP subproblem is solved. This solution procedure is systematized by starting the enumeration at a lower bound of the final time defined as

and, finally (e) either of the target constraints 𝐱(𝑡𝑓 ) = (𝑥𝑑 , 𝑦𝑑 , 𝑧𝑑 ), OR 𝐱(𝑡𝑓 ) ∈ 𝛺(𝑥𝑑 , 𝑦𝑑 , 𝑧𝑑 )

(29)

𝐱(𝑡𝑖 ) ∈ 𝜒 ⊂ R3 ,

(30)

3

(31)

𝐯(𝑡𝑖 ), 𝐚(𝑡𝑖 ) ∈ R ,

𝑡𝑓0 = 𝐷𝐱0 →𝐱𝑓 ∕(‖𝑉 ‖ + 𝑢max ),

respectively applied according to whether a strict or relaxed target formulation is considered. We refer to the above problem as the full NLP problem, and as reduced NLP problem the problem with 𝐮(𝐱(𝑡𝑖 )) = 0.

(32)

where 𝐷𝐱(𝑡0 )→𝐱(𝑡𝑓 ) is the distance between the starting point and the √ 𝑣2𝑚𝑎𝑥 + 𝑣2𝑚𝑎𝑥 is the maximum magnitude of the destination, ‖𝑉 ‖ = 𝑧 vehicle velocity, and 𝑢max is the maximum magnitude of the ocean current velocity. The solution procedure enumerates a finite set of final times, and for each final time, it solves an NLP subproblem with fixed final time. If a solution is not found, a new finite set of final times is then postulated. This procedure is implemented iteratively, where in each iteration 𝑀 values of the final time are postulated and 𝑀 NLP subproblems are solved in parallel. For the first iteration, we consider the travel times {𝑡𝑓0 , 𝑡𝑓0 + 𝛥𝑡, … , 𝑡𝑓0 + (𝑀 − 1)𝛥𝑡}. If multiple solutions are found from distinct final times, the solution corresponding to the minimum final time is stored; otherwise 𝑀 additional NLP subproblems are solved, each with a final time from {𝑡𝑓0 + 𝑘𝑀, 𝑡𝑓0 + (𝑘𝑀 + 1)𝛥𝑡, … , 𝑡𝑓0 + ((𝑘 + 1)𝑀 − 1)𝛥𝑡}, 𝑘 = 1, 2, …. Note that this iterative procedure is repeated until a solution is found, or until it reaches the last element 𝑡𝑁 in the set  . This procedure is applied to solve the reduced NLP problem in the initialization step and the full NLP problem in the optimization step.

4. Solution approach The NLP formulation presented in Section 3 is a non-convex optimization problem with unknown final time and final conditions. Besides, the formulation includes an indicator function 𝑔 and a current field function 𝑢 that are not available algebraically to the NLP solver. In preliminary experiments, the full NLP problem posed difficulties to different NLP solvers to find an initial feasible solution, and consequently converge to locally optimal solutions. Therefore, we developed robust solution strategy compromising the following three procedures: • a solution procedure to solve the full and reduced NLP problems, whereby a set of NLP problems with fixed final time, denoted by NLP subproblems, are solved to find the minimum time, see Section 4.1; • an initialization step based on the utilization of a reduced NLP problem and on the enforcement of one waypoint, which is used only in this step; • a diversification strategy based on the definition of alternative waypoints that induces alternative starting paths, see Section 4.3.

4.2. Initialization procedure The full NLP problem includes a set of variables defined over the time horizon that are difficult to initialize without using a systematic approach. However, it is well known that NLP solvers require good initialization schemes to find feasible solutions in order to reduce the computational time required to converge to a local solution. We propose an initialization step that is based on (1) the solution of a reduced NLP problem, which does not consider the current field; and (2) the enforcement of one waypoint in the reduced NLP problem. Note that the reduced NLP problem is solved through the solution of reduced

The solution approach is then divided into two steps: (1) initialization; and 2) optimization. In the first step, the reduced NLP problem is used to estimate trajectory plans, whereas in the second step the full NLP problem uses the estimated trajectory plans as the starting point; see Fig. 4 for a diagram of the approach. The following sub-sections describe the three procedures above and the two steps within the proposed solution strategy. 5

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 4. Flowchart of the solution procedures, illustrating the solution of the NLP problems, initialization step, and diversification strategies. The solution of NLP problems is based on the parallel solution of multiple NLP problems in parallel for a fixed final time. The initialization uses a reduced NLP problem without considering the current field. The diversification is implemented through the use of waypoints. For each waypoint a sequence of initialization step and optimization step is performed in parallel.

NLP subproblems according to the procedure described in Section 4.1. Therefore, within the initialization there are two steps:

4.3. Diversification search Given the non-convexity of the full NLP problem, we implement a space diversification over the sequence initialization, optimization. The space diversification is based on using alternative waypoints, one waypoint per sequence, and solve the sequence initialization, optimization in parallel. Note that one of these sequences does not consider a fixed waypoint. In Fig. 4, we represent each sequence initialization, optimization for an alternative waypoint by a column. The best solution obtained corresponds to the minimum final time obtained across the full NLP subproblems solved at the end of each sequence.

1. Solve each reduced NLP subproblem, with a fixed waypoint at the final time: 𝐱(𝑡𝑓 ) = (𝑥𝑤 , 𝑦𝑤 , 𝑧𝑤 ), instead of the final target; 2. Solve each reduced NLP subproblem, with 𝐱(𝑡𝑓 ) = 𝐱𝑓 and 𝐱(𝑡𝑤 ) = (𝑥𝑤 , 𝑦𝑤 , 𝑧𝑤 ), where 𝑡𝑤 is the final time obtained from the previous reduced NLP subproblem solution. The main output of the above procedure is a starting trajectory plan to initialize the variables 𝐱(𝑡𝑖 ), 𝐯(𝑡𝑖 ) and 𝐚(𝑡𝑖 ) in the full NLP problem considering the ocean current. 6

Ocean Engineering 188 (2019) 106266

S. Albarakati et al. Table 1 Domain size and parameters for Case 1. Parameters

Values

Units

Domain 𝜒

(𝑥min , 𝑦min , 𝑧min ) (𝑥max , 𝑦max , 𝑧max )

(1620000, 160000, −2234) (1720000, 280000, −2)

m m

Instance I

𝐱0 𝐱𝑓

(1688658, 162436, −1900) (1700000, 258958, −1800)

m m

Instance II

𝐱0 𝐱𝑓

(1700000, 258958, −1800) (1688658, 162436, −1900)

m m

𝐷𝐱0 →𝐱𝑓 𝛥𝑡 𝑁 𝑣max 𝑣max𝑧 𝑎max 𝑎max𝑧 𝑢max

97213 610 180 1 0.5 1 0.5 0.47

m s

Size of the target box

(𝑥𝑓 − (𝑥𝑓 +

Maximum time seta a

Initialization step Optimization step

m/s m/s m/s2 m/s2 m/s

𝛥𝑡⋅𝑣max 2 𝛥𝑡⋅𝑣max 2

, 𝑦𝑓 − , 𝑦𝑓 +

𝛥𝑡⋅𝑣max 2 𝛥𝑡⋅𝑣max 2

, 𝑧𝑓 − , 𝑧𝑓 +

𝛥𝑡⋅𝑣max 10 𝛥𝑡⋅𝑣max 10

) )

100 1000

s s

Maximum time set per iteration in the solution procedure described in Section 4.1.

Table 2 Lower bound of the final time for the initialization (𝑡̂𝑓0 ) and NLP (𝑡𝑓0 ) problems in Case I; top: Instance I, bottom: Instance II. Lower bound

No wp

wp1 (1705383, 219177, −1634)

wp2 (1697560, 211740, −2)

wp3 (1684301, 211410, −1423)

𝑡̂𝑓0 (h) 𝑡𝑓0 (h)

24.14 16

24.67 16

24.21 16

24.65 16

Lower bound

No wp

wp1 (1687334, 201337, −221)

wp2 (1697560, 211740, −2)

wp3 (1704523, 238212, −1509)

𝑡̂𝑓0 (h) 𝑡𝑓0 (h)

24.14 16

24.34 16

24.68 16

24.51 16

wp: waypoint.

Table 3 Case1: size of the models and computed results. Instance I Step

Initialization

#EQ

#VAR

No wp

wp1

wp2

wp3

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

2168

1812

27.11

80

27.96

81

27.11

10

28.30

84

2306 2711

2004 2355

26.60 26.43

1067 1235

26.60 26.60

1113 1211

26.77 26.60

1369 1431

26.94 26.77

1205 1164

#EQ

#VAR

No wp 𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

Optimization Point target Region target Instance II Step

Initialization Optimization Point target Region target

wp1

wp2

wp3

2168

1812

27.11

18

27.45

42

32.53

3

28.30

84

2306 2711

2004 2355

22.71 22.54

1113 1235

22.71 22.37

1154 1264

22.77 22.54

1369 1401

22.71 22.54

1197 1258

wp: waypoint; #EQ: number of equations; #VAR: number of variables; T: wall clock time (𝑠); bold: best solution.

this purpose. In all cases, the simulations are performed on a workstation with a 24-core, 2.5 GHz processor and 128 GB RAM. The stopping criteria for each iteration of the procedure described in Section 4.1 is 100 s, and 1000 s for each iteration of the same procedure in the optimization step. The simulations explored both the strict and relaxed target constraints. The solution of the latter constraint has a travel time that is smaller or equal to that obtained with the former.

4.4. Optimization procedure Th full NLP problem considers the specified ocean current field and the trajectory planning obtained in the initialization step. The solution of the full NLP problem relies on the solution of a set of subproblems, as described in Section 4.1. Note that the waypoints are only part of the initialization step and are dropped in the full NLP problem. 4.5. Implementation details

4.6. Remarks

In the application below, we have utilized the local solver IPOPT 3.12 (Wächter and Biegler, 2006) to solve each reduced and full NLP subproblem, and relied on the GAMS 24.8.2 modeling framework for

1. The procedure outlined in Section 4.1 to find the minimum time uses a systematic enumeration approach starting from a 7

Ocean Engineering 188 (2019) 106266

S. Albarakati et al. Table 4 Domain size and parameters for Case 2. Parameters

Values

Units

Domain 𝜒

(𝑥min , 𝑦min , 𝑧min ) (𝑥max , 𝑦max , 𝑧max )

(635500, 1708032, −2000) (697772, 1770400, −2)

m m

Instance I

𝐱0 𝐱𝑓 𝐷𝐱0 →𝐱𝑓

(638683, 1767854, −700) (688182, 1721856, −2) 67575.54

m m m

Instance II

𝐱0 𝐱𝑓 𝐷𝐱0 →𝐱𝑓

(696814, 1735314, −2) (639176, 1710196, −1740) 62897.3

𝛥𝑡 𝑁 𝑣max 𝑣max𝑧 𝑎max 𝑎max𝑧 𝑢max

500 180 1 0.5 1 0.5 0.57

Size of the target box

(𝑥𝑓 − (𝑥𝑓 +

Maximum time seta a Maximum

Initialization step Optimization step

m s m/s m/s m/s2 m/s2 m/s

𝛥𝑡⋅𝑣max 2 𝛥𝑡⋅𝑣max 2

, 𝑦𝑓 − , 𝑦𝑓 +

𝛥𝑡⋅𝑣max 2 𝛥𝑡⋅𝑣max 2

, 𝑧𝑓 − , 𝑧𝑓 +

𝛥𝑡⋅𝑣max 10 𝛥𝑡⋅𝑣max 10

) )

100 1000

s s

time set per iteration in the solution procedure described in Section 4.1.

Table 5 Lower bound of the final time for the initialization (𝑡̂𝑓0 ) and NLP (𝑡𝑓0 ) problems in Case II; top: Instance I, bottom: Instance II. Lower bound

No wp

wp1 (642065, 1745911, −100)

wp2 (647694, 1748955, −200)

wp3 (646646, 1736229, −200)

𝑡̂𝑓0 (h) 𝑡𝑓0 (h)

18.56 11.87

20.39 11.87

19.14 11.87

21.03 11.87

Lower bound

No wp

wp1 (688266, 1721577, −2)

wp2 (656424, 1740288, −230)

wp3 (686798, 1712080, −130)

𝑡̂𝑓0 (h) 𝑡𝑓0 (h)

17.28 11.05

18.30 11.05

20.72 11.05

20.05 11.05

Table 6 Case 2: size of the models and computed results. Instance I Step

#EQ

Initialization

#VAR

No wp

wp1

wp2

wp3

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

1928

1612

19.58

267

20.97

233

19.86

252

21.39

280

2306 2411

2004 2095

19.17 18.89

1482 1360

20.14 20.00

1131 1225

19.44 19.31

1904 1350

20.28 19.72

1175 1456

#EQ

#VAR

No wp

Optimization Point target Region target Instance II Step

Initialization Optimization Point target Region target

wp1

wp2

wp3

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

𝑡𝑓 (h)

T (s)

1928

1612

18.06

300

18.61

261

20.97

300

20.42

288

2306 2411

2004 2095

17.64 17.64

2113 2146

18.47 18.47

2178 2283

19.58 18.75

2418 2639

18.75 20.00

2516 2116

wp: waypoint; #EQ: number of equations; #VAR: number of variables; T: wall clock time (𝑠); bold: best solution.

lower bound on the minimum time. Of course, alternative approaches can be applied that offer to yield the minimum time more efficiently. For example, starting from a large value of the minimum time and then solving the NLP problems backwards by setting decreasing minimum times. That approach requires an educated guess on the minimum time to avoid distant estimates or infeasible minimum times. Another approach may involve sampling a set of minimum times, solving the corresponding problems and based on the (in)feasible solutions decide to move up or down on the values of the minimum times. Given that such enhancements are case and sample dependent, we focused in the present demonstrations on the systematic enumeration approach starting from a lower bound on the minimum time.

2. The computational performance of the overall solution approach also depends on the length of the time step 𝛥𝑡 and 𝑀. Generally, selection of an excessively 𝛥𝑡 would lead to an increase in the problem size and in the computational solution effort. Note, however, that the selection of very small time steps can be naturally avoided in the present context, as OGCM-specified current fields reflect coarse space and time discretizations that can be effectively exploited to set an appropriate range of time steps. 3. The parallel tasks involved in the solution approach are independent and do not share data between tasks. The implementation does not have a structure with master and subproblems, and

8

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

consequently does not require task synchronization or communication. Though in some steps of the solution approach there may be a waiting time for individual tasks to terminate, the parallel implementation is quite efficient, as waiting times would also be present in a sequential approach.

results of Fig. 7, which suggest that in Instance II the AUV seeks higher elevations where favorable currents are present, whereas in Instance II the AUV essentially navigates around the ocean floor where the (adverse) current magnitude is weak. To gain additional insights into these predictions, we plot in Fig. 7 profiles of the AUV velocity. It can be observed that in both instances, the velocity of the AUV with respect to the current is constant, and coincides with the specified peak velocity. This is not surprising because the formulation seeks to minimize the total travel time and because the peak AUV velocity is substantially larger than that of the (steady) ocean current prevailing in the domain. However, for most of the mission in Instance II the AUV is operating in a condition of favorable current, whose magnitude can reach about 40% of the peak AUV velocity. On the other hand, in Instance I the AUV is navigating in regions where the ocean current is weak, and can be either favorable or adverse. It is also interesting to note that the travel time for the best case in Instance II is actually lower than 𝑡̂𝑓0 . In other words, in this instance the AUV takes less time to reach its destination than would be required in the absence of current and obstacles. Thus, the observed changes in the vertical elevation of the AUV towards regions of a favorable current result in reduction in the overall travel time. This illustrates the ability of the trajectory planning algorithm to select 3D trajectories that effectively exploit the prevailing structure of the current field. Another key feature of the trajectory planning algorithm concerns its ability to select trajectories that enable the AUV to navigate in a complex environment. This occurs in Instance I, and is further illustrated in the close-up view provided in Fig. 8. As noted earlier, in this instance the AUV travels around the seabed so as to avoid (otherwise strong) adverse currents, and as seen in Fig. 8, the trajectory selected by the solver threads a passage bounded by several neighboring obstacles.

5. Application In this section, we present a case study of 3D AUV trajectory planning for different areas in the Red Sea and Gulf of Aden. The regions are selected based on complex terrain changes and strongly spatially current fields, which enables to test the algorithm in multiple challenging settings. 5.1. OGCM field The data set used for the demonstrations below originate from an assimilative OGCM of the Red Sea (Toye et al., 2017). The underlying ensemble assimilation system operates sequentially as cycles of forecast-analysis steps, using the MIT ocean general circulation model (MITgcm) for forecasting the Red Sea circulation, and an ensemble Kalman filter (EnKF) based on the Data Assimilation Research Testbed (DART) software for updating the forecast every time new observations become available (Hoteit et al., 2013). The system domain extends from 30◦ E to 50◦ E and from 10◦ N to 30◦ N, covering the whole Red Sea, the Gulf of Suez, the Gulf of Aqaba, and the Gulf of Aden. The system was configured on a spherical grid with a resolution of 0.04◦ ×0.04◦ , resulting in 500 × 500 grid points, and 50 vertical layers ranging from 4 m at the surface to 300 m near the bottom. The model bathymetry was obtained from the gridded General Bathymetric Chart of the Ocean (GEBCO), and the model was forced with 6-hourly atmospheric fields from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis. After validating the system with available in situ and satellite remote sensing observations (Yao et al., 2014b,a), this assimilation system was integrated to generate an ensemble of 3-daily flow-fields conditioned by satellite sea surface temperature and sea surface height for the month of January 2006. From this data set, we selected a single snapshot, corresponding to the ensemble mean (or filter estimate) for 28-01-2006. The 3D velocity field corresponding to the ocean state is used for the purpose of the demonstrations below. To measure Euclidean distance in the problem domain, we consider 30◦ longitude, 10◦ latitude, and the sea level as the origin, and measure all the distances from this point.

5.3. Case 2 In this case, we tackle the problem of trajectory planning in a region of the Red Sea. As illustrated in Fig. 5, this region is on average deeper than the Gulf of Aden, and has more complex subsea terrain than the region considered in Case 1. Similar to Case 1, we consider two instances in Case 2 that correspond to distinct starting points and destinations. The input data for these instances are given in Table 4. In the initialization step, time optimal paths are generated without current, with and without waypoints; see Fig. 9. The waypoints were selected randomly to diversify different initialization and avoid local minima. The coordinates of waypoints, and the lower bounds for the initialization phase (𝑡̂𝑓0 ) and full NLP solution (𝑡𝑓0 ) are given in Table 5. Fig. 10 shows the paths for the best trajectory calculated by the trajectory planning algorithm for both instances of Case 2. The figure illustrates the 3D view of the AUV with the starting point and end point of the journey, along with the local topography and the ocean current. Quantitative results of the travel time of the trajectory planning problem and the computational cost of the algorithms are summarized in Table 6. The best solution for Instance I is 18.89 h, and the best solution for Instance II is 17.64 h. The interesting point to note here is that in order to take advantage of the strong favorable current the AUV spends most of its journey near the surface. In Instance I, the AUV rises to the surface at the start of its mission, whereas in Instance II the AUV dives towards its destination at the end of the journey.

5.2. Case 1 We first consider a trajectory planning problem in the Gulf of Aden region that is illustrated in Fig. 5. The selected region, identified in Table 1, includes several subsea obstacles, and the current field exhibits substantial spatial variability. We consider two instances of the planning problem, which are also specified in Table 1. Note that the starting location and destination location are swapped between Instance I and Instance II, whereas the remaining parameters are identical. Provided in Table 2 are the coordinates of the waypoints used in the initialization scheme, as well as the lower bounds of the travel time for both the initialization phase (𝑡̂𝑓0 ) and the full NLP solution (𝑡𝑓0 ). The paths corresponding to the best trajectories determined by the trajectory planning algorithm are illustrated in Fig. 6 for both Instances I and II; animations for these cases are provided in the Supplementary Material. The figure provides 3D views of the AUV path, the starting point and destination, as well as illustrations of the local topography and of the ocean current. Quantitative estimates of travel times and solver performance metrics for the initialization and solution phases are summarized in Table 3. The results indicate that in the best case of Instance I the AUV reaches its destination in about 26.4 h, whereas in the case trajectory of Instance II the AUV completes its mission in 22.4 h. These differences appear to be well correlated with the

5.4. Remarks 1. Even though we are using multiple initializations, the optimal solution is achieved efficiently, particularly because individual solutions are determined in parallel. Specifically, as reflected in Tables 3 and 6, the solution is reached in less than 20 min of clock time, i.e. at a small fraction of a typical mission time. 9

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 5. Bathymetry contours in the Red Sea and in the Gulf of Aden. Insets show regions selected for the case studies.

Fig. 6. Case 1; top: Instance I; bottom: Instance II. Composite plot showing a 3D view of the optimal path (green) and the obstacle locations (black). The structure of the ocean ) and the current on a fixed vertical plane is illustrated using arrows, and a color scale for the velocity magnitude as indicated. The starting point is indicated by a red circle ( target destination is specified in terms a rectangular box ( ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2. We have introduced a relaxation scheme for the AUV destination, in which the latter can be specified as a cubical target

centered at the end point. The use of a relaxed target formulation always results in a travel time that is smaller or equal than that

10

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 7. Case 1; top: Instance I, bottom: Instance II. Profiles of the AUV and ocean current velocity magnitudes. The profile of the relative AUV velocity magnitude is also plotted.

Fig. 8. Case1, Instance I. Composite closeup showing the optimal path (green line), the obstacle locations (black), and the structure of the ocean current, namely using arrows and color scale for the velocity magnitude as indicated. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

6. Conclusions

obtained for the point target case. This is intuitively obvious as relaxing constraints generally lead to improvement of the objective function. We did not, however, see significant changes in either the CPU time or the predicted solutions between the two formulations.

We have developed a framework for the 3D trajectory planning problem based on an NLP formalism that accounts for realistic ocean currents and obstacles, as well as kinematic constraints. A new general 11

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 9. Case 2; top: Instance I, bottom: Instance II. Optimal paths with and without waypoints. The waypoints are represented by yellow spheres ( a red sphere ( article.)

) , the starting point by

), and the destination by a green box ( ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this

approach that exploits OGCM outputs to specify a 3D ocean current function and to formulate obstacle avoidance constraints is proposed. This approach builds a continuous and differentiable ocean current function and a smooth indicator function applied to nodal values extracted directly from OGCM data. Overall, the utilization of the processed OGCM data within the NLP formulation enables us to address 3D trajectory planning problems under realistic settings. Furthermore, the optimization framework relies on a solution approach that involves initialization and diversification strategies to enhance the performance of the solver.

visualization tools to simultaneously illustrate predictions, the local current structure and topography.

In the present work, we focused on the minimal-travel-time trajectory planning problem, and restricted our attention to steady, deterministic ocean currents. However, the framework developed may be naturally extended to accommodate unsteady ocean currents, and work is underway to demonstrate this capability. Another avenue that we plan to pursue is to accommodate for variability of ocean forecasts, as characterized by OGCM simulation ensembles. We also plan to consider other objective functions, and consequently assess the merits of the approach for a wider class of applications. From the framework implementation perspective, we plan to explore (1) improving the diversification strategy by using stochastic search-based methods to provide feasible initializations to NLP solution approach developed, and (2) extending the proposed enumeration strategy to determine minimum times to more elaborate search strategies.

We applied the proposed framework to the solution of trajectory planning problems in the Red Sea and Gulf of Aden. The computational results highlight the ability of the solution approach to determine trajectories that (1) effectively exploit the 3D ocean current structure to minimize travel time, and (2) as needed, navigate in or around complex topography. Key features of the 3D trajectory planning framework were demonstrated through (i) quantitative analysis of predictions, (ii) assessment of solver performance, (iii) application of advanced 12

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Fig. 10. Case 2; top: Instance I, bottom: Instance II. Composite plot showing a 3D view of the optimal path (green) and the obstacle locations (black). The structure of the ocean current is illustrated using arrows, and a color scale for the velocity magnitude as indicated. The starting point is indicated by a red circle ( ) and the target destination is depicted using a green box ( ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Acknowledgments

A.1. 𝑥 and 𝑦 directions

Research reported in this publication was supported by research funding from King Abdullah University of Science and Technology (KAUST), Saudi Arabia, and used resources of the KAUST Core Labs. The authors are also grateful to Mr. Thomas Theußl and Ms. Dina Garatly for their help with visualization of the results.

Given the regular discretization in these directions, the derivatives are estimated using a second-ordered centered scheme 𝑓𝑥,𝑖𝑗𝑘 ≈ 𝑓̃𝑥,𝑖𝑗𝑘 = 𝛿𝑥− 𝑓(𝑖−1)𝑗𝑘 + 𝛿𝑥+ 𝑓(𝑖−1)𝑗𝑘 , 𝑓𝑦,𝑖𝑗𝑘 ≈ 𝑓̃𝑦,𝑖𝑗𝑘 = 𝛿𝑦− 𝑓𝑖(𝑗−1)𝑘 + 𝛿𝑦+ 𝑓𝑖(𝑗+1)𝑘 , where 𝛿𝑥+ = 1∕(2ℎ𝑥 ) = −𝛿𝑥− and 𝛿𝑦+ = 1∕(2ℎ𝑦 ) = −𝛿𝑦− .

Appendix A. Numerical estimation of the first-order derivatives A.2. 𝑧 direction Consider a tensorized grid of points (𝑥𝑖 , 𝑦𝑗 , 𝑧𝑘 )𝑖,𝑗,𝑘 and a function 𝑓 specified at each point of the grid. We denote for short 𝑓 (𝑥𝑖 , 𝑦𝑗 , 𝑧𝑘 ) = 𝑓𝑖𝑗𝑘 . The goal of this section is to estimate the first-order derivatives 𝑓𝑥 , 𝑓𝑦 , 𝑓𝑧 , 𝑓𝑥𝑦 , 𝑓𝑥𝑧 , 𝑓𝑦𝑧 , and 𝑓𝑥𝑦𝑧 at every point of the grid using second-order finite difference schemes. The grid of interest is such that 𝑥1 ≤ … ≤ 𝑥𝑛𝑥 , 𝑦1 ≤ … ≤ 𝑦𝑛𝑦 , and 𝑧1 ≤ … ≤ 𝑧𝑛𝑧 , and the spacing in the 𝑥 and 𝑦 directions is constant. We denote ℎ𝑥 = 𝑥𝑖+1 − 𝑥𝑖 , ℎ𝑦 = 𝑦𝑗+1 − 𝑦𝑗 , and ℎ𝑧𝑘 = 𝑧𝑘+1 − 𝑧𝑘 .

The finite difference schemes are more elaborate when considering irregular grid. For the sake of clarity, we omit the indices 𝑖𝑗, and denote ℎ = ℎ𝑧 . We have 𝑓𝑘+1 = 𝑓𝑘 + ℎ𝑘 𝑓𝑧,𝑘 +

ℎ2𝑘 2

𝑓𝑘−1 = 𝑓𝑘 − ℎ𝑘−1 𝑓𝑧,𝑘 + 13

𝑓𝑧𝑧,𝑘 +

ℎ2𝑘−1 2

ℎ3𝑘 6

𝑓𝑧𝑧𝑧,𝑘 + ◦(ℎ3𝑘 ),

𝑓𝑧𝑧,𝑘 −

ℎ3𝑘−1 6

𝑓𝑧𝑧𝑧,𝑘 + ◦(ℎ3𝑘−1 ).

(33) (34)

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Dividing Eq. (33) (resp. (34)) by ℎ2𝑘 (resp. ℎ2𝑘−1 ) yields at the first order ) ℎ 1 1 1 ( 𝑓𝑘+1 − 𝑓𝑘 ≈ 𝑓 + 𝑓 + 𝑘 𝑓𝑧𝑧𝑧,𝑘 , ℎ𝑘 𝑧,𝑘 2 𝑧𝑧,𝑘 6 ℎ2𝑘 ) ℎ 1 1 1 ( 𝑓𝑘−1 − 𝑓𝑘 ≈ − 𝑓 + 𝑓 − 𝑘−1 𝑓𝑧𝑧𝑧,𝑘 . ℎ𝑘−1 𝑧,𝑘 2 𝑧𝑧,𝑘 6 ℎ2

The set of equations are deduced from the following derivatives of the polynomial 𝑝(𝑢, 𝑣, 𝑤) with respect to 𝑢, 𝑣, and 𝑤:

(35) 𝑝𝑢 (𝑢, 𝑣, 𝑤) =

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 (𝛼𝑢𝛼−1 )𝑣𝛽 𝑤𝛾 ,

𝛼=1 𝛽,𝛾=0

(36)

𝑘−1

𝑝𝑣 (𝑢, 𝑣, 𝑤) =

Subtracting Eq. (36) from Eq. (35) gives ) ) 1 ( 1 ( 𝑓𝑘+1 − 𝑓𝑘 − 𝑓𝑘−1 − 𝑓𝑘 ℎ2𝑘 ℎ2𝑘−1 ℎ + ℎ𝑘 1 ≈ 𝑘−1 𝑓 + (ℎ + ℎ𝑘 )𝑓𝑧𝑧𝑧,𝑘 , ℎ𝑘−1 ℎ𝑘 𝑧,𝑘 6 𝑘−1

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 𝑢𝛼 (𝛽𝑣𝛽−1 )𝑤𝛾 ,

𝛽=1 𝛼,𝛾=0

𝑝𝑤 (𝑢, 𝑣, 𝑤) =

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 𝑢𝛼 𝑣𝛽 (𝛾𝑤𝛾−1 ),

𝛾=1 𝛼,𝛽=0

𝑝𝑢𝑣 (𝑢, 𝑣, 𝑤) =

and finally

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 (𝛼𝑢𝛼−1 )(𝛽𝑣𝛽−1 )𝑤𝛾 ,

𝛼,𝛽=1 𝛾=0

ℎ − ℎ𝑘−1 ℎ𝑘 ℎ𝑘−1 𝑓 + 𝑘 𝑓 − 𝑓 ℎ𝑘 (ℎ𝑘−1 + ℎ𝑘 ) 𝑘+1 ℎ𝑘−1 ℎ𝑘 𝑘 ℎ𝑘−1 (ℎ𝑘−1 + ℎ𝑘 ) 𝑘−1 ℎ ℎ ≈ 𝑓𝑧,𝑘 + 𝑘−1 𝑘 𝑓𝑧𝑧𝑧,𝑘 . 6 We deduce the second-order approximation of 𝑓𝑧,𝑖𝑗𝑘

𝑝𝑢𝑤 (𝑢, 𝑣, 𝑤) =

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 (𝛼𝑢𝛼−1 )𝑣𝛽 (𝛾𝑤𝛾−1 ),

𝛼,𝛾=1 𝛽=0

𝑝𝑣𝑤 (𝑢, 𝑣, 𝑤) =

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 𝑢(𝛽𝑣𝛽−1 )(𝛾𝑤𝛾−1 ),

𝛽,𝛾=1 𝛼=0

− + + 𝑓𝑧,𝑖𝑗𝑘 ≈ 𝑓̃𝑧,𝑖𝑗𝑘 = 𝛿𝑧𝑘 𝑓𝑖𝑗(𝑘−1) + 𝛿𝑧𝑘 𝑓𝑖𝑗𝑘 + 𝛿𝑧𝑘 𝑓𝑖𝑗(𝑘+1) ,

𝑝𝑢𝑣𝑤 (𝑢, 𝑣, 𝑤) =

3 ∑

𝑎𝛼𝛽𝛾 𝛼𝛽𝛾(𝛼𝑢𝛼−1 )(𝛽𝑣𝛽−1 )(𝛾𝑤𝛾−1 ).

𝛼,𝛽,𝛾=1

where − 𝛿𝑧𝑘 =−

For any 𝑢, 𝑣, and 𝑤, we deduce the second derivatives

ℎ − ℎ𝑘−1 ℎ𝑘−1 ℎ𝑘 + , 𝛿 = 𝑘 , and 𝛿𝑧𝑘 = . ℎ𝑘−1 (ℎ𝑘−1 + ℎ𝑘 ) 𝑧𝑘 ℎ𝑘−1 ℎ𝑘 ℎ𝑘 (ℎ𝑘−1 + ℎ𝑘 )

𝑝𝑢𝑢 (𝑢, 𝑣, 𝑤) =

𝑝𝑣𝑣 (𝑢, 𝑣, 𝑤) =

The cross derivative in the 𝑥𝑦 directions is approximated as: 𝜕 𝑓 , 𝜕𝑥 𝑦,𝑖𝑗𝑘 − ≈ 𝛿𝑥 𝑓𝑦,(𝑖−1)𝑗𝑘 + 𝛿𝑥+ 𝑓𝑦,(𝑖+1)𝑗𝑘 ,

𝑝𝑤𝑤 (𝑢, 𝑣, 𝑤) =

𝑎𝛼𝛽𝛾 𝑢𝛼 ((𝛽 − 1)𝛽𝑣𝛽−2 )𝑤𝛾 ,

3 ∑ 3 ∑

𝑎𝛼𝛽𝛾 𝑢𝛼 𝑣𝛽 ((𝛾 − 1)𝛾𝑤𝛾−2 ).

𝛾=2 𝛼,𝛽=0

≈ 𝛿𝑥− 𝑓̃𝑦,(𝑖−1)𝑗𝑘 + 𝛿𝑥+ 𝑓̃𝑦,(𝑖+1)𝑗𝑘 , 𝛿𝑥− 𝛿𝑦− 𝑓(𝑖−1)(𝑗−1)𝑘

3 3 ∑ ∑ 𝛽=2 𝛼,𝛾=0

𝑓𝑥𝑦,𝑖𝑗𝑘 =

Using the fact that 𝑢 (resp. 𝑣, 𝑤) depends only on 𝑥 (resp. 𝑦, 𝑧), we deduce the gradient and Hessian of the polynomial interpolant 𝑞(𝑥, 𝑦, 𝑧) ≈ 𝑓 (𝑥, 𝑦, 𝑧),

+ 𝛿𝑥− 𝛿𝑦+ 𝑓(𝑖−1)(𝑗+1)𝑘 ,

= 𝑓̃𝑥𝑦,𝑖𝑗𝑘 .

𝑞𝑥 (𝑥, 𝑦, 𝑧) = 𝑢𝑥 𝑝𝑢 (𝑢, 𝑣, 𝑤),

Similarly, we deduce approximations for the cross derivatives in the 𝑥𝑧, 𝑦𝑧, 𝑥𝑦𝑧 directions, as:

𝑞𝑦 (𝑥, 𝑦, 𝑧) = 𝑣𝑦 𝑝𝑣 (𝑢, 𝑣, 𝑤), 𝑞𝑧 (𝑥, 𝑦, 𝑧) = 𝑤𝑧 𝑝𝑤 (𝑢, 𝑣, 𝑤),

𝑓𝑥𝑧,𝑖𝑗𝑘 ≈ 𝛿𝑥− 𝑓̃𝑧,(𝑖−1)𝑗𝑘 + 𝛿𝑥+ 𝑓̃𝑧,(𝑖+1)𝑗𝑘 ,

𝑞𝑥𝑥 (𝑥, 𝑦, 𝑧) = 𝑢2𝑥 𝑝𝑢𝑢 (𝑢, 𝑣, 𝑤),

𝑓𝑦𝑧,𝑖𝑗𝑘 ≈ 𝛿𝑦− 𝑓̃𝑧,𝑖(𝑗−1)𝑘 + 𝛿𝑦+ 𝑓̃𝑧,𝑖(𝑗+1)𝑘 , 𝑓𝑥𝑦𝑧,𝑖𝑗𝑘 ≈

𝑎𝛼𝛽𝛾 ((𝛼 − 1)𝛼𝑢𝛼−2 )𝑣𝛽 𝑤𝛾 ,

𝛼=2 𝛽,𝛾=0

A.3. Cross derivatives

=

3 ∑ 3 ∑

− ̃ 𝛿𝑧𝑘 𝑓𝑥𝑦,𝑖𝑗(𝑘−1)

𝑞𝑦𝑦 (𝑥, 𝑦, 𝑧) = 𝑣2𝑦 𝑝𝑣𝑣 (𝑢, 𝑣, 𝑤),

+ ̃ + 𝛿𝑧𝑘 𝑓̃𝑥𝑦,𝑖𝑗𝑘 + 𝛿𝑧𝑘 𝑓𝑥𝑦,𝑖𝑗(𝑘+1) .

𝑞𝑧𝑧 (𝑥, 𝑦, 𝑧) = 𝑤2𝑧 𝑝𝑤𝑤 (𝑢, 𝑣, 𝑤), 𝑞𝑥𝑦 (𝑥, 𝑦, 𝑧) = 𝑢𝑥 𝑣𝑦 𝑝𝑢𝑣 (𝑢, 𝑣, 𝑤),

Appendix B. Tricubic interpolation

𝑞𝑥𝑧 (𝑥, 𝑦, 𝑧) = 𝑢𝑥 𝑤𝑧 𝑝𝑢𝑤 (𝑢, 𝑣, 𝑤), 𝑞𝑦𝑧 (𝑥, 𝑦, 𝑧) = 𝑣𝑦 𝑤𝑧 𝑝𝑦𝑧 (𝑢, 𝑣, 𝑤).

Using the technique of Lekien and Marsden (2005), we can now interpolate the function 𝑓 on the Cartesian product of intervals 𝐼 = [𝑥𝑖 , 𝑥𝑖+1 ] × [𝑦𝑗 , 𝑦𝑗+1 ] × [𝑧𝑘 , 𝑧𝑘+1 ] ⊂ R3 given the values of 𝑓 , 𝑓𝑥 , 𝑓𝑦 , 𝑓𝑧 , 𝑓𝑥𝑦 , 𝑓𝑥𝑧 , 𝑓𝑦𝑧 , and 𝑓𝑥𝑦𝑧 at the eight corners of the box. The interpolation function, 𝑝, is a third-order polynomial of the form 𝑓 (𝑥, 𝑦, 𝑧) ≈ 𝑞(𝑥, 𝑦, 𝑧) = 𝑝(𝑢, 𝑣, 𝑤) =

3 ∑

Appendix C. Supplementary data Supplementary material related to this article can be found online at https://doi.org/10.1016/j.oceaneng.2019.106266.

𝑎𝑗𝑘𝓁 𝑢𝑗 𝑣𝑘 𝑤𝓁 ,

References

𝑗,𝑘,𝓁=0

with the affine change of coordinates, 𝑦 − 𝑦𝑗 𝑧 − 𝑧𝑘 𝑥 − 𝑥𝑖 , 𝑣= , 𝑤= , 𝑢= 𝑥𝑖+1 − 𝑥𝑖 𝑦𝑗+1 − 𝑦𝑗 𝑧𝑘+1 − 𝑧𝑘

Aghababa, M.P., 2012. 3D Path planning for underwater vehicles using five evolutionary optimization algorithms avoiding static and energetic obstacles. Appl. Ocean Res. 38, 48–62. Alvarez, A., Caiti, A., Onken, R., 2004. Evolutionary path planning for autonomous underwater vehicles in a variable ocean. IEEE J. Ocean. Eng. 29 (2), 418–429. Ataei, M., Yousefi-Koma, A., 2015. Three-dimensional optimal path planning for waypoint guidance of an autonomous underwater vehicle. Robot. Auton. Syst. 67, 23–32. Bayat, B., Crasta, N., Crespi, A., Pascoal, A.M., Ijspeert, A., 2017. Environmental monitoring using autonomous vehicles: a survey of recent searching techniques. Curr. Opin. Biotechnol. 45, 76–84. Besada-Portas, E., de la Torre, L., Moreno, A., Risco-Martín, J.L., 2013. On the performance comparison of multi-objective evolutionary UAV path planners. Inform. Sci. 238, 111–125.

having constant derivatives, 𝑢𝑥 =

1 , 𝑥𝑖+1 − 𝑥𝑖

𝑣𝑦 =

1 , 𝑦𝑗+1 − 𝑦𝑗

𝑤𝑧 =

1 . 𝑧𝑘+1 − 𝑧𝑘

The change of coordinates proposed above maps I to [0, 1]3 . In order to perform the interpolation, the 64 coefficients of the tensor 𝑎 must be identified, which corresponds to the number of values that are available (𝑓 and the 7 derivatives at 8 corner points). 14

Ocean Engineering 188 (2019) 106266

S. Albarakati et al.

Schouwenaars, T., Moor, B.D., Feron, E., How, J., 2001. Mixed integer programming for multi-vehicle path planning. In: 2001 European Control Conference. ECC, pp. 2603–2608. Sethian, J.A., 1999. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials, vol. 3. Cambridge University Press. Smith, R.N., Pereira, A., Chao, Y., Li, P.P., Caron, D.A., Jones, B.H., Sukhatme, G.S., 2010. Autonomous underwater vehicle trajectory design coupled with predictive ocean models: A case study. In: Robotics and Automation (ICRA), 2010 IEEE International Conference on, pp. 4770–4777. Sörensen, K., 2015. Metaheuristics-the metaphor exposed. Int. Trans. Oper. Res. 22 (1, SI), 3–18. Soulignac, M., 2011. Feasible and optimal path planning in strong current fields. IEEE Trans. Robot. 27 (1), 89–98. Soulignac, M., Taillibert, P., Rueher, M., 2008. Adapting the wavefront expansion in presence of strong currents. In: Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on, pp. 1352–1358. Subramani, D.N., Haley, P.J., Lermusiaux, P.F.J., 2016. Energy-optimal path planning in the coastal ocean. J. Geophys. Res. Oceans 122 (5), 3981–4003. Subramani, D.N., Lermusiaux, P.F., 2016. Energy-optimal path planning by stochastic dynamically orthogonal level-set optimization. Ocean Model. 100, 57–77. Sun, F., Xu, W., Jin, L., Li, J., 2010. Path planning of autonomous underwater vehicles for optimal environmental sampling. In: OCEANS 2010 IEEE - Sydney. pp. 1–4. Thompson, D.R., Chien, S., Chao, Y., Li, P., Cahill, B., Levin, J., Schofield, O., Balasuriya, A., Petillo, S., Arrott, M., Meisinger, M., 2010. Spatiotemporal path planning in strong, dynamic, uncertain currents. In: Robotics and Automation (ICRA), 2010 IEEE International Conference on, pp. 4778–4783. Toye, H., Zhan, P., Gopalakrishnan, G., Kartadikaria, A.R., Huang, H., Knio, O., Hoteit, I., 2017. Ensemble data assimilation in the Red Sea: Sensitivity to ensemble selection and atmospheric forcing. Ocean Dyn. 67 (7), 915–933. Wächter, A., Biegler, L.T., 2006. On the implementation of an interior-point filter linesearch algorithm for large-scale nonlinear programming. Math. Program. 106 (1), 25–57. Wang, T., Le Maître, O.P., Hoteit, I., Knio, O.M., 2016. Path planning in uncertain flow fields using ensemble method. Ocean Dyn. 66 (10), 1231–1251. Wang, T., Lima, R.M., Giraldi, L., Knio, O.M., 2019. Trajectory planning for autonomous underwater vehicles in the presence of obstacles and a nonlinear flow field using mixed integer nonlinear programming. Comput. Oper. Res. 101, 55–75. Warren, C.W., 1990. A technique for autonomous underwater vehicle route planning. In: Autonomous Underwater Vehicle Technology, 1990. AUV ’90. Proceedings of the (1990) Symposium on, pp. 201–205. Witt, J., Dunbabin, M., 2008. Go with the flow: Optimal AUV path planning in coastal environments. In: Proceedings of the 2008 Australasian Conference on Robotics and Automation, ACRA 2008, pp. 1–9. Wynn, R.B., Huvenne, V.A., Bas, T.P.L., Murton, B.J., Connelly, D.P., Bett, B.J., Ruhl, H.A., Morris, K.J., Peakall, J., Parsons, D.R., Sumner, E.J., Darby, S.E., Dorrell, R.M., Hunt, J.E., 2014. Autonomous underwater vehicles (AUVs): Their past, present and future contributions to the advancement of marine geoscience. Mar. Geol. 352, 451–468, 50th Anniversary Special Issue. Yang, L., Qi, J., Song, D., Xiao, J., Han, J., Xia, Y., 2016. Survey of robot 3D path planning algorithms. J. Control Sci. Eng. 2016, 7426913. Yao, F., Hoteit, I., Pratt, L., Bower, A., Köhl, A., Gopalakrishnan, G., Rivas, D., 2014a. Seasonal overturning circulation in the Red Sea: Part 2. Winter circulation. J. Geophys. Res. Oceans 119, 2263–2289. Yao, F., Hoteit, I., Pratt, L., Bower, A., Zhai, P., Köhl, A., Gopalakrishnan, G., 2014b. Seasonal overturning circulation in the Red Sea: Part 1. Model validation and summer circulation. J. Geophys. Res. Oceans 119, 2238–2262. Yilmaz, N.K., Evangelinos, C., Lermusiaux, P.F.J., Patrikalakis, N.M., 2008. Path planning of autonomous underwater vehicles for adaptive sampling using mixed integer linear programming. IEEE J. Ocean. Eng. 33 (4), 522–537. Yuh, J., 2000. Design and control of autonomous underwater robots: A survey. Auton. Robots 8 (1), 7–24. Zeng, Z., Lammas, A., Sammut, K., He, F., Tang, Y., 2014. Shell space decomposition based path planning for AUVs operating in a variable environment. Ocean Eng. 91, 181–195. Zeng, Z., Lian, L., Sammut, K., He, F., Tang, Y., Lammas, A., 2015. A survey on path planning for persistent autonomy of autonomous underwater vehicles. Ocean Eng. 110 (A), 303–313.

Blidberg, D.R., Turner, R.M., Chappell, S.G., 1991. Autonomous underwater vehicles: Current activities and research opportunities. Robot. Auton. Syst. 7 (2), 139–150. Borrelli, F., Subramanian, D., Raghunathan, A.U., Biegler, L.T., 2006. MILP and NLP techniques for centralized trajectory planning of multiple unmanned air vehicles. In: 2006 American Control Conference, VOLS 1-12. In: Proceedings of the American Control Conference, vol. 1-12, p. 644+. American Control Conference 2006, Minneapolis, MN, JUN 14–16, 2006. Dijkstra, E.W., 1959. A note on two problems in connexion with graphs. Numer. Math. 1 (1), 269–271. El Mohtar, S., Hoteit, I., Knio, O., Issa, L., Lakkis, I., 2018. Lagrangian Tracking in stochastic fields with application to an ensemble of velocity fields in the Red Sea. Ocean Model. 131, 1–14. Elisseeff, P., Schmidt, H., Johnson, M., Herold, D., Chapman, N.R., McDonald, M.M., 1999. Acoustic tomography of a coastal front in Haro Strait, British Columbia. J. Acoust. Soc. Am. 106 (1), 169–184. Ferguson , D., Stentz , A.T., 2006. Anytime RRTs. In: Proceedings of the 2006 IEEE/RSJ International Conference on Intelligent Robots and Systems. IROS ’06, pp. 5369–5375. Garau, B., Alvarez, A., Oliver, G., 2005. Path planning of autonomous underwater vehicles in current fields with complex spatial variability: an A* approach. In: Proceedings of the 2005 IEEE International Conference on Robotics and Automation, pp. 194–198. Garau, B., Bonet, M., Alvarez, A., Ruiz, S., Pascual, A., 2009. Path planning for autonomous underwater vehicles in realistic oceanic current fields: Application to gliders in the Western Mediterranean Sea. J. Maritime Res. 6 (2), 5–22. Hart, P.E., Nilsson, N.J., Raphael, B., 1968. A formal basis for the heuristic determination of minimum cost paths. IEEE Trans. Syst. Sci. Cybern. 4 (2), 100–107. Hollinger, G.A., Choudhary, S., Qarabaqi, P., Murphy, C., Mitra, U., Sukhatme, G.S., Stojanovic, M., Singh, H., Hover, F., 2012. Underwater data collection using robotic sensor networks. IEEE J. Sel. Areas Commun. 30 (5), 899–911. Hooker, J.N., 1995. Testing heuristics: We have it all wrong. J. Heuristics 1 (1), 33–42. Hoteit, I., Hoar, T., Gopalakrishnan, G., Collins, N., Anderson, J., Cornuelle, B., Köhl, A., Heimbach, P., 2013. A MITgcm/DART ensemble analysis and prediction system with application to the Gulf of Mexico. Dyn. Atmos. Oceans 63, 1–23. Koay, T.B., Chitre, M., 2013. Energy-efficient path planning for fully propelled AUVs in congested coastal waters. In: OCEANS - Bergen, 2013 MTS/IEEE. pp. 1–9. Kruger, D., Stolkin, R., Blum, A., Brìganti, J., 2007. Optimal AUV path planning for extended missions in complex, fast-flowing estuarine environments. In: Proceedings of the IEEE Int. Conf. on Robotics and Automation, pp. 4265–4270. Kuffner, J.J., LaValle, S.M., 2000. RRT-connect: An efficient approach to singlequery path planning. In: Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065), vol. 2, pp. 995–1001. Lekien, F., Marsden, J., 2005. Tricubic interpolation in three dimensions. Internat. J. Numer. Methods Engrg. 63 (3), 455–471. Li, D., Wang, P., Du, L., 2019. Path planning technologies for autonomous underwater vehicles-A review. IEEE Access 7, 9745–9768. Liu, C., Zhao, Y., Gao, F., Liu, L., 2015. Three-dimensional path planning method for autonomous underwater vehicle based on modified firefly algorithm. Math. Probl. Eng. 2015, 561394. Lolla, T., Lermusiaux, P.F.J., Ueckermann, M.P., Haley, P.J., 2014. Time-optimal path planning in dynamic flows using level set equations: theory and schemes. Ocean Dyn. 64 (10), 1373–1397. Lolla, T., Ueckermann, M.P., Yigit, K., Haley, P.J., Lermusiaux, P.F.J., 2012. Path planning in time dependent flow fields using level set methods. In: Robotics and Automation (ICRA), 2012 IEEE International Conference on, pp. 166–173. Ma, Y.-N., Gong, Y.-J., Xiao, C.-F., Gao, Y., Zhang, J., 2019. Path planning for autonomous underwater vehicles: An ant colony algorithm incorporating alarm pheromone. IEEE Trans. Veh. Technol. 68 (1), 141–154. Petres, C., Pailhas, Y., Patron, P., Petillot, Y., Evans, J., Lane, D., 2007. Path planning for autonomous underwater vehicles. IEEE Trans. Robot. 23 (2), 331–341. Petres, C., Pailhas, Y., Petillot, Y., Lane, D., 2005. Underwater path planing using fast marching algorithms. In: Europe Oceans 2005, vol. 2. pp. 814–819. Schmidt, H., Bellingham, J.G., Johnson, M., Herold, D., Farmer, D.M., Pawlowicz, R., 1996. Real-time frontal mapping with AUVs in a coastal environment. In: OCEANS 96 MTS/IEEE Conference Proceedings. the Coastal Ocean - Prospects for the 21st Century, vol. 3, pp. 1094–1098.

15