# Management of a fleet of autonomous underwater gliders for area coverage: From simulation to real-life experimentation

Aurélien Merci, Cédric Anthierens, Nadège Thirion-Moreau, Yann Le Page

## ► To cite this version:

Aurélien Merci, Cédric Anthierens, Nadège Thirion-Moreau, Yann Le Page. Management of a fleet of autonomous underwater gliders for area coverage: From simulation to real-life experimentation. Robotics and Autonomous Systems, 2024, 183, pp.104825. 10.1016/j.robot.2024.104825. hal-04943841

HAL Id: hal-04943841

https://hal.science/hal-04943841v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Management of a fleet of autonomous underwater gliders for area coverage: From simulation to real-life experimentation ☆

Aurélien Merci<sup>a</sup>, Cédric Anthierens<sup>c,*</sup>, Nadège Thirion-Moreau<sup>d</sup>, Yann Le Page<sup>b</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Unmanned Underwater Vehicles (UUV)
"SeaExplorer" Underwater Gliders (UG)

Coverage path planning

Automatic management of a UUV fleet

Graph methods

Mission plan

## ABSTRACT

This article deals with underwater gliders whether there are operated in a fleet or individually. They constitute the most affordable and energy-saving autonomous observation/data acquisition platform, making long-duration ocean exploration missions possible. In this article, theoretical researches are led to solve the path planning problem of multi-point exploration missions of this type of vehicle. We focus on the area coverage type missions i.e. all points of a given area must be visited only once. We suggest a new path planning method for area coverage i.e. the fleet of glider is sized and the optimized glider trajectories are calculated according to selected criteria (mission duration, energy consumption or traveled distance). Our proposed approach combines weighted graph theory with our underwater glider simulator whose main interest is to be capable of integrating time-varying 3D environmental data (4D). Our method is tested in simulation and then in a dynamic real-life context (Mediterranean Sea) on Alseamar's SeaExplorer autonomous underwater gliders. Finally, a comparison with the expertise of a glider pilot and a more conventional approach, exploiting only the distance between the waypoints in the operation area, confirms the relevance and effectiveness of the suggested method. The experimental mission demonstrates the interest and benefits of the approach and the ease of operational implementation in an industrial context.

## 1. Introduction

This article deals with the automatic trajectory planning of underwater vehicles, which, in our case, are underwater gliders (UG). One objective is to automatically build their mission plan in order to make the mission safer and more efficient to achieve its goal of mapping. Basically the mission plan is edited by pilots that can rely on 3D map of currents and bathymetry to make their decisions. So the mission plans are often based on empirical results whereas they can be optimized several ways depending on the considered criterion. This optimization problem is part of "Routing Problems" (see [1] for an overview of this type of problems). The VRP (Vehicle Routing Problem) is a combinatorial optimization and operational research problem belonging to the category of transport problems. It is a subset of mathematical optimization that focuses on discrete problems. It consists in finding the best solution from a finite set of possibilities [2] i.e. the optimal routes for several vehicles visiting a set of locations. When a single vehicle is considered, the problem is reduced to a Traveling Salesman Problem (TSP) [3].

Widely used in logistics, planning, networks design and routing, production management, information technology and other fields where decisions have to be made, combinatorial optimization generally involves the maximization or minimization of a well-chosen cost function under specific constraints linked to the considered problem. To obtain an optimal solution in a discrete, finite set, one method is to examine all the possible solutions and compare them to identify the best one. However, the number of possible solutions may be too large to be compared in polynomial time.

Some combinatorial optimization problems can be solved exactly in polynomial time, using techniques such as dynamic programming [4] or reformulation as a linear optimization problem however most of these problems cannot be solved in these ways. Yet, they can be reformulated using different approaches among which are integer linear optimization, graph theory [5] and so on. Here, the approach that we suggest is based on weighted graph theory. Moreover, since in most cases these problems turn out to be NP-hard, they require the use of dedicated solvers. In our case, OR-tools [6] will be used.

<footnote><sup>a</sup> DGA TN, Av. de la Tour Royale, Toulon, France</footnote>

<footnote><sup>b</sup> Alseamar, 60 Av. O. Perroy, Rousset, 13790, France</footnote>

<footnote><sup>c</sup> COSMER lab, UTLN, Av. de l'Université, La Garde, 83130, France</footnote>

<footnote><sup>d</sup> LIS Toulon, UTLN, AMU, CNRS UMR, 7020, France</footnote>

<footnote>☆ This work has been carried out in the framework of a CIFRE-DGA grant (n°2020858) allowed by the French Ministry of Army - Defense Innovation Agency.</footnote>

<footnote>* Corresponding author.
E-mail address: cedric.anthierens@univ-tln.fr (C. Anthierens).</footnote>
Area coverage also known as “Coverage Path Planning” (CPP) [7] consists in finding a path that passes through all the points in an area or volume of interest only once, while avoiding possible obstacles. CPP is quite common in UG missions like for the detection of hydrocarbon seepage [8]. In such missions, the underwater gliders have to acquire measurements at depth in each of the cells of a predefined zone. Although CPP is a fairly classic UG mission, very few solutions have been suggested to solve it especially for an UG fleet [9,10]. Moreover such missions still rely on glider pilots experience instead of being based on an automatically generated mission plan.

In this article, our aim is to solve a VRP in the context of area coverage considering one or multiple underwater gliders. The objective is thus to generate an optimized trajectory for a single glider as well as for a set of several gliders (or fleet of gliders) in order to make the whole mission safer and more efficient. Our aim is thus to develop a method that can automate the construction of the mission plan for a fleet of gliders, i.e. to calculate the optimized trajectories in terms of a certain criterion for all the gliders in the fleet. To achieve this objective, it is suggested to solve the CPP problem by combining weighted graph theory with the simulator originally developed for SeaExplorer gliders [11], which takes into account a 4D variable environment of the area to be covered. This area may contain exclusion zones. This approach involves the construction of a cost function to minimize. This cost function can depend on one of the following quantities: the mission duration, the energy consumption of each UG or the distance traveled by one or the fleet of UG.

The real novelty of the approach that we suggest is the use of our UG simulator as a direct internal model to have an easy access to these quantities of interest and also to better take into account the operation area environmental data (ocean currents, bathymetry levels, exclusion areas) and energy consumption constraints related to the SeaExplorer gliders (a deep discharge below 15% must be avoided). Our approach provides an optimal deployment and recovery point and computes an optimized trajectory between these two points for each glider while achieving the coverage of the whole area. When these deployment and recovery points are already provided, it generates the most optimized trajectory for each glider between these points. It also allows the fleet sizing by avoiding the drop to critical levels of energy for each glider.

Finally, it should be noted that in this article, we will focus on the mission preparation phase (offline mode) which consists in automatically editing the optimal mission plan for the fleet. The supervision of the mission itself (online mode) whose aim is to deal with possible unexpected events such as the lost, failure or alarm of a glider or the variability of the environmental conditions (large ocean current variations within the whole water column) is another issue that will be addressed in an upcoming article.

In the remainder of this article, we will use the following abbreviations:

### Abbreviations

* AUV: Autonomous (or semi-autonomous) Underwater Vehicles

* UUV: Unmanned Underwater Vehicles

* UG: Underwater Gliders

* CPP: Coverage Path Planning

* TSP: Traveling Salesman Problem

* (C or D)VRP(TW): (Capacity constraints or Dynamic) Vehicle

* Routing Problem (with Time Windows)

* N4, N8, N16: Four, Eight or Sixteen Neighbors

* GLIMPSE: GLIder Mission Piloting SystEm

Finally, this article is organized as follows. After an introduction of the problem under consideration and its context, the second section is dedicated to recalls on underwater gliders and on the state of the art of coverage path planning algorithms. It also reviews the objectives

pursued through the management of a fleet of underwater gliders, as well as the constraints related to this type of vehicle. Then, in the third section, some necessary notions of graph theory and modeling of vehicle routing problems are recalled. The fourth section details the suggested approach based on weighted graph theory and the construction of adjacency matrices specific to the problem under consideration. Three different ways to build these adjacency matrices are introduced and explained. We will also explain how this new approach can be used to find optimal deployment and recovery points if desired, in addition to the fleet sizing. In the fifth section, numerical simulations are performed to confirm the relevance and effectiveness of the proposed method. The results are discussed and also compared with the mission plan designed by a glider pilot. The sixth section is devoted to real-life implementation of the method on a fleet of four underwater gliders operating in the Mediterranean Sea. Finally a conclusion is drawn and perspectives are given as upcoming work to extend this method to other use cases.

## 2. Recalls on gliders, typical missions and coverage path planning

### 2.1. Basics on gliders

It was H. Stommel in 1989 in [12] who first came up with the concept of underwater gliders. These devices are, nowadays, classified as Autonomous (or semi-autonomous) Underwater Vehicles (AUVs). Their purpose was to acquire data profiles (conductivity, temperature, pressure, etc.) by moving through the water column on large scales covering thousands of kilometres during weeks or months of mission. To that end, they have to be easy to use, customizable with different sensors and more “affordable” than most other underwater vehicles. Instead of thrusters, they rely on a ballast system to manage their depth allowing to alternate between up and down movements, hence the name “yo” often used to characterize their flight profile. Various technological breakthroughs have led to the development of the first four underwater gliders on the market: “Slocum” [13], “Spray” [14], “Seaglider” [15] and “SeaExplorer” [11,16] developed by Alseamar which is the subject of this study.

Among the most recent research works on UG, some of them are devoted to the dynamic modeling process to improve the dynamic simulation accuracy [11,17]. Other aim at improving their structure and the materials that are used [18,19] or the control of these vehicles [20,21]. [22–24] presented optimization methods of the UG control parameter values for improving the UG motion accuracy and energy utilization rate simultaneously. Other works are dedicated to the development of path planning algorithms capable of taking into account a dynamic environment [9,25–27]. For an overview of the recent advances in path planning for underwater gliders see [28]. Some works are also dedicated to the management in a fleet whether it involves formation control, coordination or cooperation of multiple UG [29,30]. Finally, another important industrial issue is also to make them dive deeper.

The torpedo-shaped SeaExplorer glider is hydro-dynamically profiled for long-duration missions. Its front cone carries the payload, while the rear is dedicated to navigation control. Finally, the tail is a foldable mast that allows the glider to transmit its position and data to a ground base and receive orders each time it is surfacing.

### 2.2. Gliders typical missions & objectives of mission planning

Underwater Gliders (UGs) are an efficient solution for a wide range of military or civil missions, involving either a single vehicle or a fleet of vehicles. Two kind of gliders networks can be considered: coordinated and cooperative-unaware networking [31]. In the case of coordinated networks, each glider influences the behavior of the other gliders (e.g. master-slave relationship or geometric shape formation). In the cooperative-unaware case, each glider executes its own part of a global mission without influencing the mission of the other gliders.
The aim of most single-vehicle and multi-vehicles glider missions is to update oceanographic data provided by physical models or satellite images. Yet, the considered areas are often so vast or the duration of the mission so long that the energy autonomy of a single glider is insufficient to cover them and, thus, require the use of a fleet of UG. An important constraint for UG is related to the autonomy of the vehicles: a battery level below 15% for each glider must be avoided to guaranty that the communication and localization functions remain effective.

To achieve UG missions, path planning algorithms taking into account a time-varying environment are often considered [25-27,32-36]. In this case, the aim is to find the optimal path between two points: a starting one and an arrival one. Yet, the area coverage [9,10,37,38] also turns out to be a very important type of mission particularly for mapping (coverage of bathymetric surfaces [39]), deep ocean exploration or search purposes (for example in the purpose of oil and gas [8] for searching hydrocarbon seepages). In such a problem, a cooperative-unaware networking is usually used, especially in the case of underwater gliders, since they do not all surface at the same time and thus they cannot communicate with each other (their communication is asynchronous via a land base).

In the single vehicle case, the UG has to perform a measurement acquisition in depth by validating once all the cells of the grid corresponding to the operation area. In the multi-vehicle case, the cells in the area to be covered are distributed among the different UG. Underwater gliders typical missions often imply the investigation of several hundred kilometres wide areas with a fixed spatial sampling that can basically vary from 5 km to 12 km. It may require several dives between two cells in the zone. Classical dives are up to 1000 m deep. As their missions can last for more than one month, the mission plan often involves a significant number of waypoints. Most of the time, the mission plan is built by glider pilots but this preliminary stage can become an extremely complicated task since gliders trajectories, due to their low speed, may be very sensitive to water currents. For this reason, the spatio-temporal variability of the 3D environment in which the UG operate [40,41] has to be taken into account in this mission preparation process.

Usually, while the mission is running, the pilot can adjust the flight profile of the glider depending on the water current of its surrounding area: a slow mode for energy saving in current free area, a normal mode which corresponds to the nominal speed and a fast mode to escape strong water currents. Depending on the piloting mode, the vertical speeds $V_z$ and the pitch angles vary:

* slow flight profile: $V_z = 0.10 \text{ m/s}, pitch = 15^\circ$,

* normal flight profile: $V_z = 0.15 \text{ m/s}, pitch = 20^\circ$,

* fast flight profile: $V_z = 0.25 \text{ m/s}, pitch = 30^\circ$.

These flight profiles have been selected empirically by Alseamar, based on their piloting experience of underwater gliders. An analogy could be made with a car gearbox but with only 3 possible modes, which is sufficient for our needs. The parameter values of these flight profiles could be further optimized as in [22-24].

Obviously, the preparation mission task becomes even much more complex and difficult as the number of gliders increases. So far, very few methods and tools could be found in the literature [41] to help the pilots in this preliminary mission preparation stage, hence the interest in methods to automate this task. Yet, in [41], even if the impact of ocean currents is taken into account, the lack of a glider simulator does not make it possible to automate the generation of a mission plan. Moreover, the CPP is not addressed.

Here, we focus on the problem of the automatic mission plan generation considering area coverage type missions. It means:

* to determine the size of the UG fleet required for the mission,

* to find the optimal deployment and recovery points unless they are already provided (chosen by gliders pilots, customers or shipowners),

* and finally to calculate the optimal or slightly sub-optimal single or multi-UG trajectories between these two points

under the constraints that all the points of the area, discretized according to a given spatial sampling, have to be visited one single time while taking into consideration bathymetry and possible obstacles such as coastline borders or area prohibited to navigation. To that aim, the environmental conditions have to be taken into account and more specifically currents evolving over time and over the entire water column. In our case, the data of the Marine Copernicus database [42] are used. Moreover, energy consumption constraints related to the SeaExplorer gliders also have to be taken into account: a deep discharge below 15% must be avoided.

To tackle that problem, a simulator of underwater gliders mission has been developed first. We briefly recall in the next section what it provides as output data (for more details see [11]).

## 2.3. On the interest of an UG simulator

Based on the glider's behavioral model and the model of the environment in which it is evolving (3D bathymetry and 4D sea current), the simulator computes the main variables of interest: travel time, energy consumption, distance traveled for a given trajectory for several flight profiles between any given waypoint and its neighbors. As shown in [11], the results of a mission tested using the developed simulator are very close to those obtained during the real mission, which confirms the relevance and good calibration of the simulator, making it a particularly interesting tool for trajectory planning. Not all UG functionalities have yet been implemented, resulting in a slight time lag in the results (the simulated glider is faster than the real one).

In a weighted graph theory based approach, we will show that it facilitates the construction of the weighted adjacency matrices, which are the representation used to minimize the cost of the route in relation to the selected variable: distance, time or energy.

As well as obtaining the optimized trajectory (planning level), the simulator also helps to improve the navigation task by evaluating the effectiveness of different flight profiles. Indeed, to avoid a strong water current at shallow depths, it is sometimes more profitable to fly at greater depths. Finally, the optimized trajectory can then be transferred to Alseamar "GLIMPSE" mission management system with a flight profile adapted throughout the trajectory.

## 2.4. State of the art of coverage path planning algorithms

The existing state-of-the-art Coverage Path Planning (CPP) algorithms dedicated to robotics can be classified into different categories: deterministic algorithms versus stochastic algorithms, classical algorithms versus heuristic-based algorithms, etc. In [43], this second choice of classification is made, with algorithms that can be either deterministic or stochastic. Classical algorithms aim to find the optimal solution, i.e. the best possible solution according to given criteria. They are appropriate for problems in which it is possible to exhaustively explore all possible solutions and often guarantee the optimal solution, but can be computationally time-consuming, particularly for complex problems. These algorithms include the Random Walk (RW) [44-46], the spanning-tree coverage (STC) [47,48], Dynamic Programming (DP) [49] and many more.

Heuristic algorithms, on the other hand, are more flexible. Their objective is to find a solution close to the optimum within a reasonable time, but with no real guarantee of finding the optimum solution. They are commonly used for complex problems where it is not possible to examine all the solutions in a reasonable time. They often rely on rules of thumb, approximations, or intelligent search strategies to quickly converge towards an acceptable solution. These algorithms include the greedy algorithm [50], graph search algorithms such as the depth-first search (DFS), the Dijkstra's algorithm, A* [51], D* [52], Theta* [53].
The notable differences between these algorithms are their cost in terms of computing time and their application cases.

If this CPP problem has been extensively studied in the aerial and ground robotics fields, it is less the case in the underwater robotics field. Concerning this last field, more specifically for AUVs, most CPP algorithms focus on the operation area definition itself since complex area are considered. Cell decomposition coverage methods [54, 55], boustrophedon cellular decomposition [56] or Glasius Bio-inspired Neural Network (GBNN) [57] have been used, then, a boustrophedon path is often chosen to cover the different regions obtained with the previous methods that make up the complex area under investigation. In our case, the focus is not on the definition of the operation area itself but on the generation of trajectories for each UG. Moreover, regarding the operation area, a grid map is built by discretizing the underwater environment into cells of the same size and shape on a single area. Each cell can have two states: obstacle (exclusion zone) or target. Concerning CPP in the specific case of UGs, relatively few solutions have been suggested. A comparison of different generic solvers on an asymmetric TSP with the same pre-fixed start and end points is performed on the mission considered in [38]. Yet, the authors only consider one single glider within an area limited to a maximum of 30 waypoints. In [9,10], an ant colony bio-inspired optimization solution is proposed, yet, the impact of ocean currents is not taken into account and cells can be revisited.

Before presenting the details of the suggested method to build automatically UG mission plan, first, some necessary notions of graph theory and modeling of vehicle routing problems are presented in the next section.

# 3. Graph theory and routing problems notions for coverage path planning

## 3.1. Graph theory

In the case of area covering, the area is represented by a grid of $N \times M$ waypoints to be visited by the fleet of UG. This grid can be represented by a graph whose nodes are the waypoints considered.

As indicated in [58], a graph $G = (V, E)$ is a mathematical structure composed of two sets $V$ and $E$, respectively called vertices and edges. Each edge is associated to vertices, which correspond to the waypoints in the center of the cells of the grid. In graph theory, an adjacency matrix is used to represent which vertices of a graph are adjacent to which other vertices [59]. The adjacency matrix is a $NM \times NM$ square matrix, denoted by $A = (a_{ij})$ whose elements are defined by:

$$a_{ij} = \begin{cases} 1 & \text{if there is an edge from node } j \text{ to node } i \\ 0 & \text{otherwise} \end{cases} \quad (1)$$

To create the edges between the vertices of our representation in order to fill the adjacency matrix, different types of “neighborhood” ($Nk$ with $k = 4, 8$ or $16$) can be considered on the grid of waypoints with a Cartesian mesh, which is usually used to discretize a geographical area:

1. N4: connections in row and in column, the four nearest neighbors are considered

2. N8: connections in row, column and diagonal, the eight nearest neighbors are considered

3. N16: connections in row, column, diagonal as well as N4 connections of each diagonal

4. Complete graph: all the connections between all the nodes of the graph are admissible

This notion of connectivity [60] is widely used in image processing. Fig. 1 gives an overview of the connections involved in our method.

It is possible to associate a weight $w_{ij}$ with each edge of a graph in order to represent travel costs between nodes ($i$ and $j$) when passing through these edges (these costs can be linked to a simple distance calculation or to other quantities like time for instance). It is therefore possible to build a weighted undirected graph (the edges have no direction) by calculating a weight on the edges between the nodes.

The weighted adjacency matrix remains a $NM \times NM$ square matrix $A = (a_{ij})$ whose elements are defined by:

$$a_{ij} = \begin{cases} 0 & \text{if there is an edge between node } j \text{ and node } i \\ 0 & \text{if } i = j \\ z & \text{otherwise} \end{cases} \quad (2)$$

with $z$ a fixed number with a high value ($10^{13}$ for example, to penalize or prohibit certain paths).

## 3.2. Modeling of routing problems

Routing problems are related to logistics and involve the need for one or more vehicles to navigate in a transport network in order to deliver goods to customers or to cover network routes. The VRP (Vehicle Routing Problem) has been widely studied in the past [61] due to its complexity of solution and the multiple practical applications in logistics. A second dimension concerns the economic and environmental benefits associated with reducing transport costs. Solving this type of problem involves determining a set of routes, each carried out by a single vehicle that starts and ends at its own depot, in such a way that all customer requirements are met, all operational constraints are respected and the overall cost of transport is minimized. The VRP was initially defined by G. B. Dantzig and J. H. Ramser in [62] as “The Truck Dispatching Problem” which is a generalization of the traveling Salesman Problem. The latter is a simplified variant in a single vehicle case.

According to [63], the classic VRP can be defined as follows:

Let $G = (V, A)$ be a graph where $V = 1...n$ is a set of vertices representing the $n$ cities whose starting point (also called “depot”) is located at vertex 1, and $A$ is the set of arcs. For all arcs $(i, j)$ with $i \neq j$, it is possible to define a non-negative weighted adjacency matrix $C = (c_{ij})$. In some contexts, $c_{ij}$ can be interpreted as a travel cost or travel time between nodes in the graph. The VRP involves designing a set of least-cost vehicle routes such that:

1. Each town in $V$ (not including the depot) is visited exactly once by a single vehicle

2. All vehicle routes start and end at the depot

3. Certain constraints are satisfied

Existing constraints include restrictions on vehicle capacity [63]: a request for a certain volume to be transported is attached to each city (not including the depot) and the sum of the weights of any vehicle route must not exceed the vehicle capacity. VRPs subject to capacity constraints are called “CVRPs”. There are also variants of the VRP in which the constraint consists in a total journey time for each vehicle that must not be exceeded (“DVRP” stands for Dynamic VRP).

Time windows can also be defined for certain cities to be visited (“VRPTW” stands for Vehicle Routing Problem with Time Windows). A precedence relationships between pairs of cities can also be defined: city $i$ may have to be visited before city $j$ (to avoid having to unload to pick up something from the bottom of the transport vehicle, for example). These different constraints can be used to define variants of the VRP. Several mathematical formulations of these variants are given by [64], as an example a model of VRP known as “vehicle flow formulations” is given here which is a model frequently used for basic versions of VRP. According to [1], this model uses integer variables, associated with each arc or edge of the graph, which count the number of times the arc or edge is crossed by a vehicle. They are particularly suited to cases where the cost of the solution can be expressed as the sum of the costs associated with the arcs.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Neighborhoods 4, 8, 16, complete graph.</figcaption></figure>

This model can be described as an integer program. In the latter, $c_{ij}$ represents the cost of moving from node $i$ to node $j$, $x_{ij}$ is a binary variable that takes the value 1 if the arc from $i$ to $j$ is considered part of the solution and 0 otherwise, $K$ is the number of vehicles available and $r(S)$ is the minimum number of vehicles needed to serve the set $S$. It is also assumed that 0 is the depot node. This formulation is given by [1]:

$$ \text{Minimization: min} \quad \sum_{i \in V} \sum_{j \in V} c_{ij} x_{ij} \qquad (3a) $$

$$ \text{Under constraints:} \quad \sum_{i \in V} x_{ij} = 1, \quad \forall j \in V \setminus \{0\} \qquad (3b) $$

$$ \sum_{j \in V} x_{ij} = 1, \quad \forall i \in V \setminus \{0\} \qquad (3c) $$

$$ \sum_{i \in V} x_{i0} = K, \quad \forall i \in V \qquad (3d) $$

$$ \sum_{j \in V} x_{0j} = K, \quad \forall j \in V \qquad (3e) $$

$$ \sum_{i \notin S} \sum_{j \in S} x_{ij} \geq r(S), \quad \forall S \subseteq V \setminus \{0\}, S \neq \emptyset \qquad (3f) $$

$$ x_{ij} \in \{0, 1\}, \quad \forall i, j \in V \qquad (3g) $$

Constraints specify the conditions that the solution must satisfy. Constraints (3b) and (3c) stipulate that an arc enters and exits each vertex associated with a customer. The constraints (3d) and (3e) stipulate that the number of vehicles that leave the depot is the same as the number of vehicles that enter it. The constraint (3f) stands for a capacity cut constraint requiring that the routes are connected and that the demand on each route does not exceed the capacity of the vehicles. Finally, the last constraint (3g) specifies that the variables $x_{ij}$ are binary: if a connection between the cities $i$ and $j$ is established, the variable is worth 1 and 0 otherwise.

A wide variety of solvers have been developed to obtain approximate solutions of variants of the VRP (e.g. [65]). We have shown that our problem can be modeled as a VRP and we have chosen to use OR-tools [6], the solver developed by Google, to solve it. It will be presented in the Section 4.3.

## 4. Details of the suggested approach

As part of our developments, a general framework was designed in Python Object language. It takes as input model an area to cover represented by a grid of waypoints to be visited by UG. As explained previously in the Section 3.1, this grid can be represented by a graph [58] whose nodes are the considered waypoints. An adjacency matrix [59] (cf. Section 3.1) is built from this modeling using one of the different types of neighborhood (cf. Section 3.1) between the nodes. Based on the connections created between the nodes and depending on the variable of interest chosen (duration, energy or distance), it is possible to establish a cost of these connections. It is represented by the weighting adjacency matrix (cf. Section 3.1). This cost can be estimated by different methods. Here, we introduce three different ones. Either our UG missions simulator is used as a direct internal model or a distance or a time cost estimation between the nodes connected by edges is

performed. These three methods will be presented in Section 4.2. Then a generic solver named OR-Tools [6] is used to solve the considered problem. The result is a feasible and optimized path going through all the nodes of the defined graph. Finally, this path can be tested on our UG missions simulator in order to compare the obtained variables of interest and to confirm the relevance of the built paths. Fig. 2 provides a synoptic diagram of the suggested coverage path planning framework. Its main blocks/functions are detailed in the next subsections.

### 4.1. Preparation of the input data

During the mission preparation phase, SeaExplorer's pilots define each waypoint to target by its latitude and longitude coordinates and by a validation distance [11]. The underwater gliders are very sensitive to water currents. To properly operate the SeaExplorer gliders, the pilot should take into consideration a 4D varying environment (latitude (N), longitude (M), immersion (z) and time t) which can be composed of exclusion zones such as coasts or areas with low bathymetry or prohibited to navigation. The developed method takes as input the geographical positions of the waypoints (vectors of N latitudes and M longitudes) to be reached by the glider, a 4D water current dataset and a 3D bathymetry dataset in the NetCDF format of the zone of interest. Concerning water current and bathymetry data, they can be retrieved from the Marine Copernicus database [42], where the water current data are available from 01/11/2020 to the present day, with 9 days of forecasts.

The method generates a $N \times M$ grid of the area to cover with the $NM$ waypoints. Their validation distance is the center of the cells of the grid. The 4D water current dataset is then prepared in order to obtain an averaged 3D-tensor of water current over the area of operation: at each point in the water “cube” on the dates available in the 4D dataset (see Fig. 3), the East and North components of the currents at each point are averaged over time.

During the grid generation phase, the method uses the water current dataset and the bathymetry level dataset to check if all the waypoints are reachable or whether some of them are part of the exclusion areas. If a waypoint has to be excluded from the zone, it is removed from the grid. Once the grid of the area of operation has been created, a graph can be generated from which it will be possible to construct different adjacency matrices depending on the considered neighborhood, see Section 3.1 (resp. weighted adjacency matrices depending on the considered neighborhood and the chosen cost (see next Section 4.2)).

### 4.2. Construction of different weighted adjacency matrices

The presence of strong water currents in the operation area remains the main constraint in underwater gliders missions: it is cheaper to go as the water current than against it. This is why the binary weight in the adjacency matrix has to be replaced by a weight with a real physical meaning. Thus, we suggest to take into account both the behavioral model of the glider and the model of the environment (3D bathymetry and 4D water current datasets) to simulate a glider trajectory under the influence of the water current and as close to reality as possible.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Overview of the suggested coverage path planning framework.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Partitioning the 4D ocean current data (latitude, longitude, immersion, time) into n “cubes” according to the n corresponding dates.</figcaption></figure>

Three different approaches to calculate the weights $w_{ij}$ (see Section 3.1) have been developed. They are presented in the next subsections. The first method simply consists in using the Haversine formula [66] to calculate the distances between all the nodes of the generated complete graph. It does not take into account physical parameters (either linked to the vehicle or to the environment). On the contrary, the two other methods exploit a water current model averaged over the duration of the planned mission (see Section 4.1). It is necessary to use an average current information because it is not possible to know in advance what the optimum sequence of waypoints will be, and, therefore, when the simulated glider will actually pass through a given waypoint in the area to be covered. The first of these two methods is based on the simulator presented in [11] and the second one relies on a travel time based cost estimation. Finally, it can be seen that the weighted adjacency matrices built with these last two methods are asymmetrical since the cost between the node i to j is not the same as the cost between the node j to i due to the influence of the current, whereas the first one is symmetrical.

### 4.2.1. Distance based cost estimation

From the area to cover used as input of the model, the method extracts the geographical positions of all the waypoints composing the grid. From these positions, a Haversine formula is used to calculate the distance between nodes of the complete graph generated from the area to cover. The distance between two waypoints i and j is given by Eq. (4):

$$w_{ij} = D_{ij} = E_r \cdot \left(2 \cdot \arctan2(\sqrt{a}, \sqrt{1-a})\right) \quad (4)$$

with $a = \sin^2(\frac{\text{Lat}_j - \text{Lat}_i}{2}) + \cos(\text{Lat}_i) \cdot \cos(\text{Lat}_j) \cdot \sin^2(\frac{\text{Lon}_j - \text{Lon}_i}{2})$, $E_r$ stands for the Earth's radius (mean radius = 6371 km), $\text{Lat}_i$ and $\text{Lon}_i$ (resp.

$\text{Lat}_j$ and $\text{Lon}_j$) are the geographical positions of the waypoint i (resp. waypoint j).

This method is used to calculate the distance on each connection of the adjacency matrix. In this case, the optimization is performed on the distance covered by the SeaExplorer gliders. The aim is to minimize this distance.

### 4.2.2. Simulator as a direct internal model

It has been shown in [11] that our UG mission simulator is relevant and reliable to estimate the most important variables of interest such as the battery level, the mission duration and the distance traveled by the SeaExplorer glider based on the behavioral model of the glider and the models of the environment. It is possible to simulate any path with different flight profiles between any waypoints. Our coverage path planning method has been designed to be able to call this simulator in order to calculate the weight of each edge by taking as inputs the operation area and the adjacency matrix of the associated graph. For each edge, the method uses the geographical coordinates of the two nodes connected by the considered edge to start a SeaExplorer simulation of a straight trajectory between those two nodes with a given flight profile. Once the simulation is completed, the result obtained on the metric to be minimized (duration, energy or distance) is saved. This value is used as weight on the considered edge. This method is the most flexible since the optimization can be performed on one of the three different calculated metrics, yet, it remains the most expensive in terms of computing time. With this method, we go through all the non-zero elements of the adjacency matrix to establish the value of all the weights $w_{ij}$ of the weighted adjacency matrix.

### 4.2.3. Travel time based cost estimation

To keep calculation times under 1 h, the use of the previous method is limited to neighborhoods 8 and 16. To overgo this limit, a faster
<figure><img src="image_7.png" /><figcaption>Fig. 4. Vector projection.</figcaption></figure>

approximations based method has been developed. It consists in calculating an estimate of the time required to complete the path between all the pairs of nodes connected by the edges. Instead of relying on the simulator as a direct internal model, it uses a vector projection between the displacement vector $V_d$, the mean water current vector $V_c$ (with $\alpha$ the angle between $V_c$ and $V_d$) and the glider speed vector $V_g$ (see Fig. 4).

Let consider the geographical coordinates of the two waypoints $Wp_i$ and $Wp_j$, the distance $D_{ij}$ and the bearing (clockwise from North) between them. The path is spatially discretized with a chosen spatial sampling $d$ to provide $N$ estimations of the mean current vector on the water column between both considered waypoints (see Fig. 5) with $d = \frac{D_{ij}}{N}$. An average of the water currents calculated along the route between the two waypoints is generated. It corresponds to the norm of the mean water current vector $V_c$. The norm of the vector $V_g$ corresponds to the horizontal speed of the SeaExplorer. It is computed from the hydrodynamic model which involves the target vertical speed and the target pitch as inputs. Using the vector projection depicted in Fig. 4, the following second order polynomial $V_d^2+V_d \times (-2V_c \cos(\alpha))+(V_c^2-V_g^2) = 0$ has to be solved to get the norm of the displacement vector, $V_d$. Then, the time-cost of each edge (complete graph) is obtained with: $w_{ij} = \frac{D_{ij}}{V_g}$. This method is cheaper in terms of calculation time than the previous method but it is limited to mission duration optimization.

It is also possible to use this method to complete the weighted adjacency matrix obtained with the simulator as a direct internal model. Instead of using a fixed high value for $z$ in Eq. (2), it can be calculated using this method to obtain a complete graph.

Once the weighted adjacency matrix has been built, the OR-Tools solver presented in the next subsection is used to solve the considered VRP. It can handle problems involving a large number of waypoints.

### 4.3. OR-tools resolution

OR-Tools is an open-source solver developed by Google [6]. It is a software suite for optimization in constraint programming, linear programming, vehicle routing and so on.

Our method uses the vehicle routing solver of OR-Tools which generates an optimal route, a route being a sequence of nodes to visit. First, it creates a data model which stores the data for the problem to solve. In our case, the problem at hand is to find an optimized path that goes through all the nodes of the operation area once, either a single glider or a fleet of gliders is considered. The data model is composed of three elements: the number of vehicles to use, a “depot” parameter which corresponds to the start and end location of the route and a cost matrix which, in our case, is the weighted adjacency matrix generated for a given metric to be minimized by the solver.

Once the data model instantiated, the method creates a routing model thanks to the "pywrapcp" module of the API. The vehicle routing solver is using a "transit callback" function to take any pair of nodes and to return the cost of the edges between those nodes by retrieving this information directly in our weighted adjacency matrix.

In the case of multi-vehicles path planning, a “dimension” parameter is added in the solver parameters which computes the cumulative cost of the route traveled by each vehicle. This parameter allows the solver to find a route with approximately the same global cost for each vehicle of the fleet.

The search process of the solver is based on a “first solution strategy” approach. By this way the solver finds an initial solution which will be optimized during the search. Several heuristic methods are available in the solver such as “PATH_CHEAPEST_ARC” or “PATH_MOST_CONSTRAINED_ARC”, whose descriptions are given in [67]. It is also possible to use as parameter a local search such as the “GUIDED_LOCAL_SEARCH” metaheuristic in order to escape local minima during the search process. A time limit for the solver has to be set when this metaheuristic is used.

Initially, the vehicle routing part of the solver was designed to solve the traveling salesman problem which corresponds to the search for the best (shortest) route going through all the nodes of the graph once and returning to the starting point (cycle). However, to set different start and end locations for each vehicle involved in the problem, it remains possible to modify the data model provided as input to the solver. A path is then built since there is no return to the starting node. To that aim, two vectors containing the indices of the start and end locations need to be added in the data model.

Instead of fixing the deployment and recovery points to generate such a path, another approach can be considered. We can still search for a cycle, but the idea is to add a dummy start and end node in the original graph used to create the adjacency matrix. This dummy node will be removed later. To be included in the weighted adjacency matrix, it has to have edges connected to all the other nodes of the graph but these edges must have a null weight. Finally, it implies to add a row and a column of zeros in the weighted adjacency matrix of the graph. By setting the depot at this dummy node, the solver will find a cycle starting at the dummy node and ending in the same node. By deleting the first step and the final step of the generated route, the dummy node is eliminated and it becomes possible to find the best deployment and recovery nodes.

### 4.4. Outputs generation

The obtained route for each glider of the fleet is saved by the method and is converted to a “Mission” object for each glider as input of the UG mission simulator [11]. Once the global mission plan has been optimized, our simulator uses multi-threading to generate the trajectories simulations for the fleet of gliders.

From the grid of the operation area to cover (Section 4.1), it has access to all the geographical positions of the waypoints composing this grid as well as their validation distance ($V_d$). In this way, the method generates a list of target waypoints in the format: $[Lat_{wp}, Lon_{wp}, V_d]$, where $Lat_{wp}$ (resp. $Lon_{wp}$) stands for the latitude (resp. the longitude). After creating the waypoints list, it generates the “segment” attribute of the Mission object as input of the simulator. This attribute consists in setting a list of initial waypoints, waypoints to target with a given flight profile number and a given payload number considered for each segment to be carried out by the simulated glider on the optimized route generated. The same flight profile (vertical speed, pitch, immersions, number of yos underwater, time at the surface to communicate) and payload profile (activation/deactivation of sensors) given as input of the method for the creation of the adjacency matrix are used for the outputs generation. Therefore, the bathymetry level dataset and the water current dataset are the same as the datasets used for the generation of the weighted adjacency matrices with the simulator. But, instead of using the 3D water current model averaged over the whole duration of the planned mission, the 4D model with a starting date is given as input of the simulator. Thus it is possible to simulate the influence of the real spatio-temporal variation of the water current influencing the trajectory of the simulated glider. The variables of
<figure><img src="image_8.png" /><figcaption>Fig. 5. Calculation of the weighted average current.</figcaption></figure>

interest, such as battery level, mission duration and distance covered by the simulated glider can be predicted as closely as possible to the ground truth, as shown in [11].

Now, to illustrate the relevance and effectiveness of the proposed method, numerical simulations are performed in the next section. In the last section, the method is then tested under real operating conditions in the Mediterranean Sea.

# 5. Numerical simulations

In this case study, we consider an operation area located off the coast of Iceland. This zone is very disturbed in terms of water currents with areas of low bathymetric level. The considered mission is an area coverage on a grid of 211 cells. It started the 25th of October 2021. The dimensions of the operation area are 201 km × 333 km with a 12.57 km of spacing between the nodes along the longitudinal axis and 27.8 km of spacing between the nodes along the latitudinal axis. An averaged mean water current on the duration of the mission is represented by a field of vectors in gray on the figures. The length of the vectors is proportional to the norm of the water current. The available energy of the simulated glider is also represented on the figures using a linear variation of the color between green (100% of battery) and red (15% of battery) then the black color is used below 15% of battery.

## 5.1. The single glider case

### 5.1.1. Case of deployment and recovery points chosen by a glider pilot

First, given the operation area to cover, a pilot has been asked to prepare this mission. He has to set a time optimized route for one

glider from a starting waypoint to a final waypoint arbitrarily chosen by accounting the 4D water current dataset of the area [68] and the bathymetry levels. The pilot chooses the starting node (depicted by a blue star) on the North-West of the area and the ending node (depicted by a red star) in the South-East. As observed on Fig. 6, the pilot draws a lawn-mower trajectory for the glider, choosing a normal flight profile.

This path has been played on our UG mission simulator with the starting date of the simulation and the 4D water current dataset already used by the pilot in the mission preparation phase. The output results have been generated in 483 s and can be found in Table 1.

From now on, all the resolutions involved in our method, are performed using the same OR-tools parameters: "PATH_MOST_CONSTRAINED_ARC" heuristic, "GUIDED_LOCAL_SEARCH" metaheuristic and running time of 900 sec.

Given the start and end nodes selected by the pilot, we use our coverage path planning method in order to obtain a distance optimized path with a normal flight profile between these two nodes. As expected, we obtain the same lawn-mower trajectory drawn by the pilot. This means that, unconsciously, the pilot optimizes the trajectory in distance rather than in time.

Then, a time optimized path with a normal flight profile between these two nodes is calculated with our method. This path is obtained with an adjacency matrix generated with the N8 neighborhood and using the simulator as a direct internal model (Section 4.2.2 with the duration metric) to generate the weighted adjacency matrix used as input of the resolution. The obtained path is shown in Fig. 7.

To be compared, the obtained paths have been tested on the simulator with the same inputs as the simulation of the pilot plan previously described. The output results can be found in Table 1.
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Lawn-mower trajectory generated by a glider pilot between a starting point (represented by a blue star) and an ending point (red star). The path color indicates the available energy from green (100% of battery) to red (15% of battery). Black color means below 15% of battery. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

As observed in this table, the comparison of our coverage path planning method with the expertise of a real glider pilot emphasizes the strength of the method for the mission preparation even when managing a single glider. The time-optimized trajectory generated by the suggested method saves approximately 43 out of 144 days of mission time initially planned with the same deployment and recovery points. This is due to the complexity for the pilot to define a trajectory that takes into consideration the strong water currents to which the glider might be confronted along its whole trajectory. One advantage of our method is to take the water current and the hydrodynamic model of the glider into consideration directly through the weighted adjacency matrix. As observed in Fig. 7, the trajectory generated by the method is often following the water current vectors. If it is not the case, the path is more expensive locally but this ultimately results in a globally optimal path. Finally, by comparing the path generated by the pilot and the time-optimized trajectory obtained thanks to the method, one can notice that on a total of 210 segments, about 72% of them are in the same direction, 22% are the same but in the opposite direction and only 6% of them are in a totally different direction.

### 5.1.2. Optimized deployment and recovery points

Previously the starting and ending node were chosen by the pilot. Yet, as explained in Section 4.3, the suggested method can also determine the optimal deployment and recovery points on the operation area. This feature of our method has been used on the study case with a normal flight profile. The result is shown in Fig. 8.

This result was obtained with an adjacency matrix generated with a N8 neighborhood and by using the simulator as a direct internal model (with the duration metric) as well as the time based cost estimation method in order to extend the adjacency matrix to a complete graph for the resolution. The output obtained through the test of the obtained route on the simulator are presented in the Table 1. Taking into account the sea current in a smart way allows the method to find a more optimal trajectory than a pilot would choose. The trajectory obtained with an optimal deployment and recovery points generation presented in Fig. 8 leads to a saving of approximately 46 days of mission time compared to the pilot's trajectory. Yet, although deploying the simulated glider towards the center of the area and by recovering it in the North-East of the zone is more optimal, this is not necessarily a very convenient or economical solution in most operational cases, since in most glider missions, deployments close to the coast are preferred.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Time-optimal trajectory generated with N8 and simulator as direct internal model (duration metric). Starting and ending points (blue and red stars) defined by the glider pilot. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 1**
Single-vehicle outputs comparison.</caption><thead><tr><td>Methods</td><td>Duration of the simulated mission</td><td>Gain</td><td>Simulated glider energy level</td></tr></thead><tbody><tr><td>Lawn-mower trajectory defined by the pilot</td><td>144 days, 12 h, 23 min</td><td>N/A</td><td>-214.72%</td></tr><tr><td>Distance-optimized trajectory with deployment and recovery points defined by the pilot</td><td>144 days, 12 h, 23 min</td><td>N/A</td><td>-214.72%</td></tr><tr><td>Time-optimized trajectory with deployment and recovery points defined by the pilot</td><td>101 days, 6 h, 35 min</td><td>-43 days</td><td>-130.63%</td></tr><tr><td>Time-optimized trajectory with optimal deployment and recovery points</td><td><b>98 days, 16 h, 17 min</b></td><td>-46 days</td><td>-130.69%</td></tr></tbody></table>

We notice that in a single vehicle case, the best trajectory in terms of energy level was obtained with the generation of the optimal deployment and recovery points with $E = -130.69\%$ of available energy at the end of the simulated mission. So we can assume that at least $i > (\frac{E}{85} + 1)$, i.e. 3 gliders are required to cover the whole area.

To conclude this first study, all the results obtained with a single vehicle clearly highlight the interest of our method to generate optimized trajectories. It also emphasizes the relevance of the fleet management when the operation area is too wide to be covered by one single vehicle, since all the time-optimized trajectories generated with our method, lead to a negative energy level of the simulated glider. When a fleet of

gliders is considered, the developed method makes it possible to test a large variety of configurations by playing upon the number of vehicles, the deployment and recovery points, the flight profiles and so on. It is illustrated in the next subsection.

## 5.2. The fleet of gliders case

Thanks to our path planning method, a distance optimized and a time optimized trajectory have been generated for each vehicle composing a fleet of 3 and 4 vehicles considering a normal flight profile.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Time-optimal trajectory obtained with the simulator as internal model on N8 neighborhood + time based cost estimation to extend to a complete graph. Estimation of the optimal starting and ending points.</figcaption></figure>

The optimal deployment and recovery points obtained in the single vehicle case were taken as fixed starting and ending points of the fleet of 3 gliders whose trajectories are shown in Figs. 9 & 10. The problem at hand is now too difficult to be solved by a glider pilot.

For obvious economic reasons and to account operational constraints, the starting and ending points are placed near to the coast as depicted on Fig. 11. Nevertheless the fleet has to be increased to 4 gliders. The output of the time-optimized trajectories presented in the Table 2 were obtained with an adjacency matrix generated with a neighborhood 8 and by using the simulator as a direct internal model as well as the time based cost estimation method in order to get an adjacency matrix of a complete graph. The distance-optimized output have been obtained with an adjacency matrix generated from a complete graph and using the distance based cost estimation. These variables of interest are the outputs of the test of the obtained routes for each glider composing the fleets on the simulator. The routes were generated with a running time of 2400 s using the “PATH_CHEAPEST_ARC” heuristic and the “GUIDED_LOCAL_SEARCH” metaheuristic (see Fig. 12).

The route were tested on the simulator and the output of these tests are presented in the Table 2.

As the problem is too difficult to be solved by a glider pilot, the distance between the waypoints have been used to fill the adjacency matrices in order to obtain distance-optimized trajectories for a fleet of three gliders as depicted in Fig. 9 and for a fleet of four gliders as shown in Fig. 11. The trajectories obtained are more affected by the water current. The outputs clearly show that not taking the water current into account when generating the matrices can lead to more

time-consuming trajectories. Without taking the current into account, it is also possible that one or several vehicles complete their mission in a critical state, as shown by the energy levels below 15% of the vehicles 1, 2 and 3 (for the fleet of three gliders), of the vehicles 4 for the fleet of four gliders. With the time-optimized mission plan generated by the method, each vehicle of the fleet of 4 gliders arrives at the recovery point near the coast at approximately the same time and takes from 31 to 33 days to complete its trajectory. This mission plan greatly reduces the duration of the coverage mission on the operation area compared to the single vehicle case. Finally, all these results emphasize the ability of our method to define precisely the size of the fleet required to cover the operation area and to automate the construction of feasible mission plan. In this example, the simulations confirm that 4 gliders have to be used to complete the mission. Playing upon flight profiles, it would be possible to further improve the remaining energy levels. The method also offers the possibility to test various sets of parameters (deployment and recovery points, fleet size, optimization criteria).

### 5.3. Impact of current averaging on results

Finally, the time-optimized trajectories of the four gliders with a deployment and recovery points near the coast have been played on the simulator but using a water current model averaged over 32 days instead of using a 4D water current dataset. The aim of this test is to evaluate the impact of the spatio-temporal variability of the marine environment. The outputs are presented in the Table 3.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Distance optimal trajectories for a fleet of three gliders calculated with a complete graph + distance based cost estimation. Starting and ending points are fixed. Blue for the first glider, green for the second glider and red for the third glider. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 2**
Multi-vehicles outputs comparison between fleets of *i* = 3 and *i* + 1 = 4 gliders.</caption><thead><tr><th>Methods</th><th>UG</th><th>Duration of the simulated mission</th><th>Simulated glider energy level</th></tr></thead><tbody><tr><td>Distance-optimized trajectories of the 3 gliders with the deployment and recovery points obtained in the single vehicle case</td><td>UG 1<br>UG 2<br>UG 3</td><td>64 days, 8 h, 36 min<br>48 days, 4 h, 19 min<br>41 days, 17 h, 15 min</td><td>-37.51%<br>14.16%<br>-19.65%</td></tr><tr><td>Time-optimized trajectories of the 3 gliders with the deployment and recovery points obtained in the single vehicle case</td><td>UG 1<br>UG 2<br>UG 3</td><td>42 days, 17 h, 18 min<br>41 days, 5 h, 1 min<br>40 days, 23 h, 8 min</td><td>-5.38%<br>18.59%<br>0.65%</td></tr><tr><td>Distance-optimized trajectories of the 4 gliders with deployment and recovery points near the coast</td><td>UG 1<br>UG 2<br>UG 3<br>UG 4</td><td>31 days, 22 h, 44 min<br>33 days, 19 h, 2 min<br>32 days, 10 h, 34 min<br>40 days, 0 h, 14 min</td><td>22.23%<br>18.79%<br>25.10%<br>0.40%</td></tr><tr><td>Time-optimized trajectories of the 4 gliders with deployment and recovery points near the coast</td><td>UG 1<br>UG 2<br>UG 3<br>UG 4</td><td>31 days, 1 h, 54 min<br>31 days, 4 h, 12 min<br>32 days, 10 h, 11 min<br>33 days, 21 h, 19 min</td><td>18.28%<br>27.84%<br>18.52%<br>26.48%</td></tr></tbody></table>

The comparison between the trajectories of the four gliders simulated on a 4D and a 3D water current dataset highlights that the approximation we made (use of a 3D water current to generate the weighted adjacency matrices, see Fig. 3) makes sense for mission preparation since the outputs presented in the Table 3 are relatively close to each other. In this example, sea current conditions certainly remained sufficiently stable (or stationary) over the entire duration of

the mission on the operation area to avoid impacting too much the results.

Yet, under real-life conditions, and in the event of major fluctuations in sea currents compared with the predictions used to calculate the optimized trajectories of the gliders' fleet, these trajectories may no longer be as optimal. This potential variability of the water current that may occur during the mission (online mode) should be taken in
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Time-optimal trajectories for a fleet of three gliders calculated with N8 + simulator + time based cost estimation. Starting and ending points are fixed. Blue for the first glider, green for the second glider and red for the third glider. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>**Table 3**

Outputs comparison between the trajectories simulated on a 4D and a 3D water current dataset.</caption><thead><tr><th>Methods</th><th>UG</th><th>Duration of the simulated mission</th><th>Simulated glider energy level</th></tr></thead><tbody><tr><td rowspan="4">Time-optimized trajectories of the 4 gliders played on a 4D water current dataset</td><td>UG 1</td><td>31 days, 1 h, 54 min</td><td>18.28%</td></tr><tr><td>UG 2</td><td>31 days, 4 h, 12 min</td><td>27.84%</td></tr><tr><td>UG 3</td><td>32 days, 10 h, 11 min</td><td>18.52%</td></tr><tr><td>UG 4</td><td>33 days, 21 h, 19 min</td><td>26.48%</td></tr><tr><td rowspan="4">Time-optimized trajectories of the 4 gliders played on a 3D water current dataset</td><td>UG 1</td><td>30 days, 19 h, 39 min</td><td>20.91%</td></tr><tr><td>UG 2</td><td>31 days, 20 h, 53 min</td><td>24.11%</td></tr><tr><td>UG 3</td><td>32 days, 7 h, 9 min</td><td>17.02%</td></tr><tr><td>UG 4</td><td>32 days, 18 h, 18 min</td><td>27.90%</td></tr></tbody></table>

consideration by a supervision method. One possibility is to update the water current model used every day as soon as the forecast data are available in order to generate an optimal path for the fleet of gliders removing the already visited waypoints. This part is beyond the scope of this article and will be presented in an upcoming article.

Finally, on an Intel(R) Core(TM) i7-9750H CPU, 16 GB of RAM, it took 5 min to simulate the whole trajectory for one glider, and about 10 min for the mission plan of 4 gliders. The global calculation time (running time of the resolution plus the simulation of the generated mission) amounts to a maximum of 50 min which is much shorter than the duration between two surfacing operations of a SeaExplorer glider (around 3 h).

## 6. Experimental validation in the mediterranean sea

In order to show the interests and advances provided by this automatic path planning method for a UG fleet, a real mission was carried out in July 2023 in Mediterranean sea. Its goal was to upgrade the acoustic propagation model within a square area of 55 km long side including a square exclusion area of 15 km long side arbitrary chosen for the test.

For the needs of this experimental campaign, the deployment and recovery points were chosen as identical and not very far from the coast to be reached easily with a light semi-rigid boat (as depicted in Fig. 13).
<figure><img src="image_14.png" /><figcaption>Fig. 11. Fleet of four gliders trajectories — Distance optimized.</figcaption></figure>

The oceanographic requirements gave a pitch of 5 km between longitude and latitude waypoints that implies 12 waypoints on each side with validation radius of 2,5 km.

The main currents in this area are a downstream for north-east parallel to the coast (liguro-provençal current) and a counterclockwise rotational one in the bottom left of the considered area. After importing the bathymetry map and the forecast of the water current, the mission has been simulated with a single glider to assess the need in energy that presizes the fleet.

It turned out that a 3-glider fleet could meet the energy criterion, but the duration of the mission would not fit the availability constrain of all the gliders. So, to make it shorter, it was decided to operate a 4-glider fleet for this mission. Our method was used to generate time-optimized trajectories for each gliders, enabling them to cover the desired area by starting from and returning to the same meeting point as depicted on Fig. 14.

Following the joint launch of the four gliders, this mission was carried out without human assistance and proceeded as planned, without any errors (see Fig. 15). All the waypoints were reached only once, and all four gliders were recovered after 6 days of mission at the desired point in a time window shorter than 8 h.

Table 4 compares the mission duration and the remaining energy of the planned mission based on the forecast current, the actual values measured after the mission is over, and the mission that would be time-optimized by relying on the actual water currents.

We note that the current forecast are accurate enough because both simulation results with actual and forecast currents are almost the

<table><caption>**Table 4**
Comparison of parameters of interest regarding the actual mission led in sea and the simulated mission.</caption><thead><tr><th>Method</th><th>Glider</th><th>Mission duration</th><th>Remaining energy</th></tr></thead><tbody><tr><td rowspan="4">Time optimized trajectories for 4 gliders with forecast current</td><td>SEA006</td><td>6 days, 12 h, 7 min</td><td>86.55%</td></tr><tr><td>SEA030</td><td>6 days, 0 h, 41 min</td><td>87.53%</td></tr><tr><td>SEA037</td><td>6 days, 8 h, 18 min</td><td>86.87%</td></tr><tr><td>SEA038</td><td>6 days, 8 h, 18 min</td><td>86.87%</td></tr><tr><td rowspan="4">Actual paths of 4 gliders carried out while the real mission</td><td>SEA006</td><td>7 days, 7 h, 17 min</td><td>86%</td></tr><tr><td>SEA030</td><td>6 days, 14 h, 42 min</td><td>87%</td></tr><tr><td>SEA037</td><td>7 days, 5 h</td><td>86%</td></tr><tr><td>SEA038</td><td>7 days, 5 h, 3 min</td><td>86%</td></tr><tr><td rowspan="4">Time optimized trajectories for 4 gliders with measured current</td><td>SEA006</td><td>6 days, 8 h, 18 min</td><td>86.87%</td></tr><tr><td>SEA030</td><td>6 days, 0 h, 41 min</td><td>87.53%</td></tr><tr><td>SEA037</td><td>6 days, 8 h, 18 min</td><td>86.87%</td></tr><tr><td>SEA038</td><td>6 days, 8 h, 18 min</td><td>86.87%</td></tr></tbody></table>

same. The gliders' order of arrival and remaining energy are the same as those simulated. Finally, the real glider is slightly slower than its optimistic simulated hydrodynamic model, this is the reason why the simulated mission durations are shorter than the actual ones.

More details and an analysis of the experimental campaign will be addressed in an up-coming article.
<figure><img src="image_15.png" /><figcaption>Fig. 12. Fleet of four gliders trajectories — Time optimized.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. Deployment of 4 UG with a semi-rigid boat.</figcaption></figure>

## 7. Conclusion & perspectives

In this article, we have presented a new method dedicated to the coverage path planning for underwater gliders fleet. This method is based on the graph theory and on the simulator of underwater gliders

missions developed for the Alseamar SeaExplorer gliders. The latter takes into consideration a 4D varying environment, the bathymetry levels of the area to cover and the hydrodynamic model of the gliders which makes it the genuine tool to calculate the weights of the edges in the weighted adjacency matrices. Thanks to its modularity, this
<figure><img src="image_17.png" /><figcaption>Fig. 14. Mission plan for 4-glider fleet.</figcaption></figure>

method is suitable for all types of underwater vehicle provided their
hydrodynamic model is available. This tool can also be used as output
of our method in order to validate the generated mission plan.

By combining the input data of the operation area, the adjacency
matrices definition, their weighted version and the resolution based
on a generic solver, the suggested method automatically produces an
optimized trajectory for each vehicle of the fleet in the area to be
covered. The latter may contain exclusion areas such as coasts or zones
prohibited for navigation. The generated trajectories prove to be much
more efficient than those provided by glider pilots. Moreover, in fleet
navigation, the trajectories generated are far too complex to be defined
by a human operator (they happen to be very different from classic
boustrophedon trajectories).

It is also possible to calculate optimal deployment and recovery
points on this area. Many scenarios can be tested on the method allow-
ing to define precisely the size of the glider fleet, and other parameters
of interest, for instance, in this article, mission duration was optimized.
A real experimental validation in the context of a multi-vehicles optimal
area coverage mission has also been performed. It was the first time that
Alseamar's SeaExplorer UG were operated completely autonomously
during a 8 day mission without any glider pilot intervention. Given the
evolving performance of current solvers, this method is even efficient
online (calculation time shorter than time between two surfacings).

A perspective of evolution would be to exploit the developed
method for other types of mission than area coverage. It would be
interesting to test the approach on other shapes of map than grids too.
The method is used, here, for the mission preparation phase (offline)
but it can also be used in the supervision phase (online). It consists in

the generation of optimized trajectories during the mission depending
on the waypoints already validated. To visit the remaining waypoints,
an updated water current model of the operation area or the water
current measured by the glider during the mission could be used. In
this supervision phase, the problem of losing a vehicle during a mission
could also be addressed.

CRediT authorship contribution statement

**Aurélien Merci:** Writing – original draft, Software, Methodology, Conceptualization. **Cédric Anthierens:** Writing – review & editing, Validation, Supervision, Resources, Methodology, Conceptualization. **Nadège Thirion-Moreau:** Writing – review & editing, Validation, Supervision, Resources, Methodology, Conceptualization. **Yann Le Page:** Supervision, Resources, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Acknowledgments

This study was conducted as part of a Cifre-Défense PhD thesis with
the Alseamar company, designer and manufacturer of the SeaExplorer
gliders. We would like to thank the French Ministry of the Armed Forces
and the Defense Innovation Agency for the technical and financial sup-
port they provided to this project. We also want to thank F. Imperadori
<figure><img src="image_18.png" /><figcaption>Fig. 15. Actual paths of the 4-glider fleet planned mission.</figcaption></figure>

and C. Vuilmet, DGA tutors of the thesis, for their support. This study has been conducted using E.U. Copernicus Marine Service Information: https://doi.org/10.48670/moi-00016.

## Data availability

The authors do not have permission to share data.

## References

[1] P. Toth, D. Vigo, An overview of vehicle routing problems, Veh. Routing Probl. (2002) 1-26.

[2] B.H. Korte, J. Vygen, B. Korte, J. Vygen, Combinatorial Optimization, Vol. 1, Springer, 2011.

[3] M. Jünger, G. Reinelt, G. Rinaldi, Chapter 4 the traveling salesman problem, in: Network Models, in: Handbooks in Operations Research and Management Science, vol. 7, Elsevier, 1995, pp. 225-330, http://dx.doi.org/10.1016/S0927-0507(05)80121-5, URL https://www.sciencedirect.com/science/article/pii/S0927050705801215.

[4] R. Bellman, The theory of dynamic programming, Bull. Amer. Math. Soc. (60) (1954) 503-516.

[5] J.A. Bondy, U.S.R. Murty, Graph Theory with Applications, American Elsevier, New York, 1976, Ch. 1.

[6] V. Furnon, L. Perron, Google OR-tools routing library, 2023, v9.7. https://developers.google.com/optimization/routing/.

[7] E. Galceran, M. Carreras, A survey on coverage path planning for robotics, Robot. Auton. Syst. 61 (12) (2013) 1258-1276.

[8] W.P. Meurer, J. Blum, G. Shipman, Volumetric mapping of methane concentrations at the bush hill hydrocarbon seep, gulf of Mexico, Front. Earth Sci. 9 (2021) 747, http://dx.doi.org/10.3389/feart.2021.604930.

[9] H. Ji, H. Hu, X. Peng, Multi-underwater gliders coverage path planning based on ant colony optimization, Electronics 11 (19) (2022) http://dx.doi.org/10.3390/electronics11193021, URL https://www.mdpi.com/2079-9292/11/19/3021.

[10] G. Han, Z. Zhou, T. Zhang, H. Wang, L. Liu, Y. Peng, M. Guizani, Ant-colony-based complete-coverage path-planning algorithm for underwater gliders in ocean areas with thermoclines, IEEE Trans. Veh. Technol. 69 (8) (2020) 8959-8971, http://dx.doi.org/10.1109/TVT.2020.2998137.

[11] A. Merci, C. Anthierens, N. Thirion-Moreau, Y. Le Page, A simulator of underwater glider missions for path planning, Ocean Eng. 269 (2023) 113514, http://dx.doi.org/10.1016/j.oceaneng.2022.113514, URL https://www.sciencedirect.com/science/article/pii/S0029801822027974.

[12] H. Stommel, The slocum mission, Oceanography (1989) URL https://doi.org/10.5670/oceanog.1989.26.

[13] D. Webb, P. Simonetti, C. Jones, SLOCUM: An underwater glider propelled by environmental energy, IEEE J. Ocean. Eng. 26 (2001) 447-452, http://dx.doi.org/10.1109/48.972077.

[14] J.T. Sherman, R.E. Davis, W.B. Owens, J.R. Valdes, The autonomous underwater glider "spray", IEEE J. Ocean. Eng. 26 (2001) 437-446.

[15] C. Eriksen, T. Osse, R. Light, T. Wen, T. Lehman, P. Sabin, J. Ballard, A. Chiodi, Seaglider: A long-range autonomous underwater vehicle for oceanographic research, IEEE J. Ocean. Eng. 26 (2001) 424-436, http://dx.doi.org/10.1109/48.972073.

[16] O. Fommervault, F. Besson, P. Lattes, SeaExplorer underwater glider: A new tool to measure water velocity, Mar. Technol. (2018).

[17] G. Jing, L. Lei, Y. Gang, Dynamic modeling and experimental analysis of an underwater glider in the ocean, Appl. Math. Model. 108 (2022) 392-407, http://dx.doi.org/10.1016/j.apm.2022.03.034, URL https://www.sciencedirect.com/science/article/pii/S0307904X22001561.

[18] Y. He, B. Song, Y. Cao, Structural optimization of an underwater glider with blended wing body, Adv. Mech. Eng. 9 (2017) 168781401772327, http://dx.doi.org/10.1177/1687814017723279.

[19] M. Yang, Y. Wang, S. Yang, L. Zhang, J. Deng, Shape optimization of underwater glider based on approximate model technology, Appl. Ocean Res. 110 (2021) 102580, http://dx.doi.org/10.1016/j.apor.2021.102580, URL https://www.sciencedirect.com/science/article/pii/S0141118721000572.

[20] R. da Silva Tchilian, E. Rafikova, S.A. Gafurov, M. Rafikov, Optimal control of an underwater glider vehicle, Procedia Eng. 176 (2017) 732-740, http://dx.doi.org/10.1016/j.proeng.2017.02.322, URL https://www.sciencedirect.com/
science/article/pii/S1877705817308317. Proceedings of the 3rd International Conference on Dynamics and Vibroacoustics of Machines (DVM2016) June 29-July 01, 2016 Samara, Russia.

[21] R. Juan, T. Wang, S. Liu, Y. Zhou, W. Ma, W. Niu, Z. Gao, High-precision motion control of underwater gliders based on reinforcement learning, Ocean Eng. 310 (2024) 118603, http://dx.doi.org/10.1016/j.oceaneng.2024.118603, URL https://www.sciencedirect.com/science/article/pii/S0029801824019413.

[22] H. Wu, W. Niu, S. Wang, S. Yan, An optimization method for control parameters of underwater gliders considering energy consumption and motion accuracy, Appl. Math. Model. 90 (2021) 1099-1119, http://dx.doi.org/10.1016/j.apm.2020.10.015, URL https://www.sciencedirect.com/science/article/pii/S0307904X20305989.

[23] H. Wu, W. Niu, Y. Zhang, S. Wang, S. Yan, Multidisciplinary optimization-based path planning for underwater gliders executing multi-point exploration missions, Ocean Eng. 266 (2022) 113022, http://dx.doi.org/10.1016/j.oceaneng.2022.113022, URL https://www.sciencedirect.com/science/article/pii/S0029801822023058.

[24] H. Wu, W. Niu, Y. Zhang, S. Wang, S. Yan, A hybrid polynomial-based optimization method for underwater gliders with parameter uncertainty, Appl. Ocean Res. 133 (2023) 103486, http://dx.doi.org/10.1016/j.apor.2023.103486, URL https://www.sciencedirect.com/science/article/pii/S014118723000305.

[25] J. Cai, F. Zhang, S. Sun, T. Li, A meta-heuristic assisted underwater glider path planning method, Ocean Eng. 242 (2021) 110121, http://dx.doi.org/10.1016/j.oceaneng.2021.110121, URL https://www.sciencedirect.com/science/article/pii/S002980182101444X.

[26] X. Liu, D. Ma, M. Yang, X. Xia, P. Guo, Modified block A\* path-planning method for hybrid-driven underwater gliders, IEEE J. Ocean. Eng. 47 (1) (2022) 20-31, http://dx.doi.org/10.1109/JOE.2021.3097414.

[27] W. Lan, X. Jin, T. Wang, H. Zhou, Improved RRT algorithms to solve path planning of multi-glider in time-varying ocean currents, IEEE Access 9 (2021) 158098-158115, http://dx.doi.org/10.1109/ACCESS.2021.3130367.

[28] R. Zhang, B. He, Y. Wang, W. Ma, S. Yang, Recent advances in path planning for underwater gliders: A comprehensive review, Ocean Eng. 299 (2024) 117166, http://dx.doi.org/10.1016/j.oceaneng.2024.117166, URL https://www.sciencedirect.com/science/article/pii/S0029801824005031.

[29] H. Wen, H. Zhou, J. Fu, X. Zhang, B. Yao, L. Lian, Consensus protocol based attitudes coordination control for underwater glider formation, Ocean Eng. 262 (2022) 112307, http://dx.doi.org/10.1016/j.oceaneng.2022.112307, URL https://www.sciencedirect.com/science/article/pii/S0029801822016055.

[30] X. Ma, Y. Wang, G. Zhang, S. Yang, S. Li, Discrete-time formation control of multiple heterogeneous underwater gliders, Ocean Eng. 258 (2022) 111728, http://dx.doi.org/10.1016/j.oceaneng.2022.111728, URL https://www.sciencedirect.com/science/article/pii/S00298018220160800.

[31] A. Alvarez, B. Mourre, Cooperation or coordination of underwater glider networks? An assessment from observing system simulation experiments in the ligurian sea, J. Atmos. Ocean. Technol. 31 (10) (2014) 2268-2277, http://dx.doi.org/10.1175/JTECH-D-13-00214.1, URL https://journals.ametsoc.org/view/journals/atot/31/10/jtech-d-13-00214.1.xml.

[32] W. Lan, X. Jin, X. Chang, T. Wang, H. Zhou, W. Tian, L. Zhou, Path planning for underwater gliders in time-varying ocean current using deep reinforcement learning, Ocean Eng. 262 (2022) 112226, http://dx.doi.org/10.1016/j.oceaneng.2022.112226, URL https://www.sciencedirect.com/science/article/pii/S0029801822015372.

[33] D. Kularatne, S. Bhattacharya, M.A. Hsieh, Going with the flow: A graph based approach to optimal path planning in general flows, Auton. Robots 42 (7) (2018) 1369-1387, http://dx.doi.org/10.1007/s10514-018-9741-6.

[34] B. Garau, A. Alvarez, G. Oliver, 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, 2005, pp. 194-198, http://dx.doi.org/10.1109/ROBOT.2005.1570118.

[35] D.N. Subramani, P.J. Haley Jr., P.F.J. Lermusiaux, Energy-optimal path planning in the coastal ocean, J. Geophys. Res.: Oceans 122 (5) (2017) 3981-4003, http://dx.doi.org/10.1002/2016JC012231, arXiv:https://agupubs.onlinelibrary.wiley.com/pdf/10.1002/2016JC012231. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016JC012231.

[36] B. Hadi, A. Khosravi, P. Sarhadi, Deep reinforcement learning for adaptive path planning and control of an autonomous underwater vehicle, Appl. Ocean Res. 129 (2022) 103326, http://dx.doi.org/10.1016/j.apor.2022.103326, URL https://www.sciencedirect.com/science/article/pii/S014118722002589.

[37] M. Panda, B. Das, B. Subudhi, B. Pati, A comprehensive review of path planning algorithms for autonomous underwater vehicles, Int. J. Autom. Comput. 17 (2020) http://dx.doi.org/10.1007/s11633-019-1204-9.

[38] T.K. Avi, Glider Mission Planning Using Generic Solvers (Ph.D. thesis), Memorial University of Newfoundland, 2014.

[39] E. Galceran, M. Carreras, Planning coverage paths on bathymetric maps for in-detail inspection of the ocean floor, in: 2013 IEEE International Conference on Robotics and Automation, 2013, pp. 4159-4164, http://dx.doi.org/10.1109/ICRA.2013.6631164.

[40] M. Troesch, S. Chien, Y. Chao, J. Farrara, J. Girton, J. Dunlap, Autonomous control of marine floats in the presence of dynamic, uncertain ocean currents, Robot. Auton. Syst. (108) (2018) 100-114.

[41] M. Eichhorn, Optimal routing strategies for autonomous underwater vehicles in time-varying environment, Robot. Auton. Syst. (67) (2015) 33-43.

[42] Marine Copernicus database, 2015, https://resources.marine.copernicus.eu/products. (Accessed 02 August 2022).

[43] C.S. Tan, R. Mohd-Mokhtar, M.R. Arshad, A comprehensive review of coverage path planning in robotics using classical and heuristic algorithms, IEEE Access 9 (2021) 119310-119342, http://dx.doi.org/10.1109/ACCESS.2021.3108177.

[44] F. Bartumeus, M.G.E. da Luz, G.M. Viswanathan, J. Catalán, Animal search strategies: a quantitative random-walk analysis, Ecology 86 (2005) 3078-3087.

[45] C. Dimidov, G. Oriolo, V. Trianni, Random walks in swarm robotics: An experiment with kilobots, in: ANTS Conference, 2016.

[46] M. Kegeleirs, D. Garzón Ramos, M. Birattari, Random walk exploration for swarm mapping, in: K. Althoefer, J. Konstantinova, K. Zhang (Eds.), Towards Autonomous Robotic Systems, Springer International Publishing, Cham, 2019, pp. 211-222.

[47] Y. Gabriely, E.D. Rimon, Spanning-tree based coverage of continuous areas by a mobile robot, Ann. Math. Artif. Intell. 31 (2001) 77-98.

[48] N. Hazon, G.A. Kaminka, Redundancy, efficiency, and robustness in multi-robot coverage, in: ICRA-05, 2005.

[49] R. Wu, The application of dynamic programming in production planning, AIP Conf. Proc. 1839 (1) (2017) 020155, http://dx.doi.org/10.1063/1.4982520.

[50] Y. Wu, J. Wang, Algorithm Design Practice for Collegiate Programming Contests and Education, CRC Press, 2018, http://dx.doi.org/10.1201/9780429401855.

[51] X. Liu, D. Gong, A comparative study of A-star algorithms for search and rescue in perfect maze, in: 2011 International Conference on Electric Information and Control Engineering, 2011, pp. 24-27.

[52] A. Stentz, Optimal and efficient path planning for partially-known environments, in: Proceedings of the 1994 IEEE International Conference on Robotics and Automation, 1994, pp. 3310-3317, http://dx.doi.org/10.1109/ROBOT.1994.351061, vol.4.

[53] A. Nash, K. Daniel, S. Koenig, A. Felnner, Theta\*: Any-angle path planning on grids, in: AAAI Conference on Artificial Intelligence, 2007.

[54] E. Galceran, M. Carreras, Efficient seabed coverage path planning for ASVs and AUVs, in: 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, 2012, pp. 88-93, http://dx.doi.org/10.1109/IROS.2012.6385553.

[55] E. Galceran, Coverage Path Planning for Autonomous Underwater Vehicles (Ph.D. thesis), University of Girona, 2014.

[56] H. Choset, P. Pignon, Coverage path planning: The boustrophedon cellular decomposition, in: A. Zelinsky (Ed.), Field and Service Robotics, Springer London, London, 1998, pp. 203-209.

[57] D. Zhu, C. Tian, B. Sun, C. Luo, Complete coverage path planning of autonomous underwater vehicle based on GBNN algorithm, J. Intell. Robot. Syst. 94 (2019) http://dx.doi.org/10.1007/s10846-018-0787-7.

[58] G. Chartrand, P. Zhang, Introduction to Graph Theory, in: Walter Rudin student series in advanced mathematics, McGraw-Hill Higher Education, 2005, URL https://books.google.fr/books?id=kea8QgAACAAJ.

[59] H. Singh, R. Sharma, Role of adjacency matrix & adjacency list in graph theory, Int. J. Comput. Technol. 3 (2012) 179-183, http://dx.doi.org/10.24297/ijct.v3i1c.2775.

[60] R.C. Gonzalez, R.E. Woods, Digital Image Processing, Prentice Hall, Upper Saddle River, N.J., 2008, URL http://www.amazon.com/Digital-Image-Processing-3rd-Edition/dp/013168728X.

[61] W. Cao, W. Yang, A survey of vehicle routing problem, MATEC Web Conf. 100 (2017) 01006, http://dx.doi.org/10.1051/mateccconf/201710001006.

[62] G.B. Dantzig, J.H. Ramser, The truck dispatching problem, Manag. Sci. 6 (1) (1959) 80-91.

[63] G. Laporte, The vehicle routing problem: An overview of exact and approximate algorithms, European J. Oper. Res. 59 (3) (1992) 345-358, http://dx.doi.org/10.1016/0377-2217(92)90192-C, URL https://www.sciencedirect.com/science/article/pii/03772219290192C.

[64] G. Laporte, Y. Nobert, Exact Algorithms for the Vehicle Routing Problem, in: North-Holland Mathematics Studies, vol. 132, 1987, pp. 147-184, URL https://apisemanticscholar.org/CorpusID:117747488.

[65] A. Pessoa, R. Sadykov, E. Uchoa, F. Vanderbeck, A generic exact solver for vehicle routing and related problems, Math. Program. 183 (2020) 483-523, http://dx.doi.org/10.1007/s10107-020-01523-z.

[66] Haversine formulas, 2002, https://www.movable-type.co.uk/scripts/latlong.html. (Accessed 02 August 2022).

[67] OR-tools routing options, 2010, https://developers.google.com/optimization/routing/routing_options#firstSolution_strategy. (Accessed 17 March 2023).

[68] Global water current model - Global_analysis_forecast phys_001_024, 2015, https://doi.org/10.48670/moi-00016. (Accessed 02 August 2022).
<figure><img src="image_19.png" /></figure>

**Aurélien MERCI** was born in 1997 in Nogent-Sur-Marne, France. He graduated from the engineering school of SeaTech and he obtained a Master's degree in robotics and connected objects, both in 2020. In 2023, he obtained his Ph.D. in underwater robotics as part of a CifreDefense grant with the Alseamar company. Most of his works focused on the development of tools enabling the cooperation between several underwater gliders of a fleet in order to fulfill a data acquisition requirement over areas with very large spatiotemporal scales. Since December 2023, he has joined the French defense procurement agency (DGA) and he is in charge of expertise and architecture for underwater drones.

<figure><img src="image_20.png" /></figure>

**Cédric ANTHIERENS** is from the north of France where he started studies on mechanical engineering. He graduated in 1996 in the engineering school of INSA of Lyon before achieving a Ph.D. thesis in 1999 on microrobotics. He has been working in Seatech Engineering school of the University of Toulon since 2014 where he is head of the mechatronics and robotics school. Most of his research are led in the COSMER lab and address issues about marine and submarine robots, umbilical and the cooperation between autonomous surface vehicle and autonomous underwater vehicle.

So, his recent works deal with path planning for area coverage with underwater gliders or the automatic management of the umbilical that links a surface robot to an underwater robot.

<figure><img src="image_21.png" /></figure>

**Nadège Thirion-Moreau** was born in 1969 in Montbéliard, France. She obtained her diploma in physics engineering from the École nationale supérieure de physique (now PHELMA), and her Master's degree in signal processing, both in 1992, from the Institut National Polytechnique de Grenoble (INPG), France. In 1995, she obtained her Ph.D. in signal processing from INPG. From 1996 to 1998, she was assistant professor at the École supérieure des procédés électroniques et optiques (now Politech), Orléans, France. Since 1998, she has joined the University of Toulon, first as Assistant Professor in the Telecom. & Computer Sciences Dept (1998–2011), then as Professor in the Signal Process-

ing - Machine Learning Dept. Since 2017, she has been a full professor at the SeaTech engineering school and head of the Signal Image Modeling team at the LIS laboratory, UMR CNRS 7020 (Computer Science and Systems).

Her main areas of research are signal analysis and processing (deterministic and statistical), machine learning and computer vision. She is particularly interested in antenna processing, blind source separation/equalization, matrix and tensor decompositions (CP, Tucker), latent variable estimation, constrained optimization, non-stationary processes, (spatial quadratic) time-frequency representations, decision/classification, object recognition, perception, weakly supervised learning for applications in maritime surveillance, defense and autonomous vehicles/drones.

<figure><img src="image_22.png" /></figure>

**Yann Le Page** Trained as a robotics engineer (UTC, France), Yann Le Page specialized in underwater robotics by obtaining a Master of Science in Ocean Engineering in the United States (FAU, Florida).

His initial experience as a research engineer at the SeaTech laboratory attached to Florida Atlantic University confirmed his passion for this field. On the strength of this experience, in 2001 he joined Cybernetix (Marseille, France), where he was in charge of a number of projects, developing a strong oil seismic activity. His greatest achievement in this field was the invention and development of the Nautilus product, now marketed by Sercel.

He then joined ACSA in 2007 to develop a subsea robotics business. Among other projects, he was responsible for developing the non-American SeaExplorer glider-submarine, while also setting up the mechanical-robotics department. The machine is now making the headlines after breaking two world records. In 2012, he took over the technical management of the company and all the design offices, then became interim co-director of the company in 2013. In 2015, ACSA, BMTI and ATOB merged to become ALSEAMAR company of which he is now Technical Director of the Products Department.