Multi-UAV Coverage Path Planning Based on Hexagonal Grid Decomposition in Maritime Search and Rescue

Abstract

In the event of a maritime accident, surveying the maximum area efficiently in the least amount of time is crucial for rescuing survivors. Increasingly, unmanned aerial vehicles (UAVs) are being used in search and rescue operations. This study proposes a method to generate a search path that covers all generated nodes in the shortest amount of time with multiple heterogeneous UAVs. The proposed model, which is a mixed-integer linear programming (MILP) model based on a hexagonal grid-based decomposition method, was verified through a simulation analysis based on the performance of an actual UAV. This study presents both the optimization technique’s calculation time as a function of the search area size and the various UAV routes derived as the search area grows. The results of this study can have wide-ranging applications for emergency search and rescue operations.

1. Introduction

In Korea, an average of 2631 marine accidents have occurred per year over the last six years (Figure 1). Out of 553 fatal accidents in 2020, 126 people were killed or are missing [1]. As the environment changes, including weather and ocean currents, marine accidents involving fishing boats or passenger ships are becoming more common [2]. Because marine accidents usually occur in environments where direct human access is difficult, the risk of death also increases, creating danger for the victim and rescuer alike [3]. To reduce the number of fatalities, a system that can respond to marine accidents under normal circumstances should be established, and first responders should be prepared to use it in an emergency [4,5].

Unlike land, the ocean is a difficult environment for humans to directly access in the event of an accident [6]. The time required for search and rescue (SAR) operations is determined by the amount and performance of available resources, and various tools and resources have been developed to address this issue [7,8,9,10,11,12]. Unmanned aerial vehicles (UAVs) are being introduced in several countries to improve maritime SAR operations [13]. Humans lose efficiency over time, whereas robots are able to maintain their efficiency at all times, thereby reducing the uncertainty that can arise in maritime SAR [14].

This study focuses on search operations that require rapid and efficient coverage of a target area in maritime SAR. A search operation, which can also be defined as a coverage path planning (CPP) problem, requires navigating a large area in an abbreviated period of time to rescue survivors. In general, robot path planning is classified into global path planning and local path planning [15,16]. Global path planning creates a complete path from the initial starting point to the target point based on deterministic or stochastic information before the mobile robot starts operating [17,18]. In contrast, local path planning generates an operation plan in real-time, responding to changes in operating environments and data collected by local sensors during navigation [19,20,21,22]. In this study, we propose a method for global CPP with deterministic information of UAVs to search for survivors at sea.

A CPP problem first divides the target area into subareas that need to be covered and then determines the order in which the divided area will be covered. As shown in Figure 2, the two methods for segmenting the region are cellular decomposition and grid-based decomposition [23,24,25]. In Figure 2a, the cellular decomposition method divides the area based on the characteristics of the search area and the UAV’s speed, weight, battery endurance, etc. [26,27,28]. Although the cellular decomposition method is simple to determine the order of cells to cover as it considers the characteristics of the region and the UAV when dividing, owing to the various characteristics, it is difficult to define the decomposition criterion. As shown in Figure 2b, the grid-based decomposition method divides an area into small pieces using a predetermined unit. The grid-based decomposition method divides the region without considering the region’s characteristics or the UAV, making it simple to divide the region but increasing the complexity of planning the coverage path.

For the CPP problem, the method of dividing a region has a direct impact on the method of coverage path planning [24,28,29]. For studies that use the cellular decomposition method [30,31,32,33,34,35], the coverage path is simply obtained by assuming that the order (i.e., path) in which the UAV moves to cover each cell is back-and-forth. The sweep direction is determined by the amount of energy the robot consumes moving in reverse when moving back and forth [26,27,28,30,33,34,35,36,37,38,39,40]. The cellular decomposition method has the advantage of quickly establishing a path at each cell using an algorithm with low computational complexity. However, when multiple UAV types are used or an irregularly shaped region is searched, this method has difficulty finding an optimal path, and a local optimal solution is derived [41,42]. When the area is divided using the grid-based decomposition method, a large number of grid cells (i.e., the number of nodes) are created. Thus, it takes a relatively long time to generate a path. However, the grid-based decomposition method does generate an optimal path [43,44,45,46,47]. This study focuses on generating the optimal search path for multiple UAVs by segmenting the target area using the grid-based decomposition method.

In most grid-based decomposition studies [43,44,45,46,47,48,49,50], the shape of the grid cell is set to a square, which is the motivation for this study. The search path may derive differently as the location and number of adjacent nodes that the UAV can move to from each node vary depending on the shape of the grid cell. Additionally, the more sides the grid cell has, the more directions it allows for movement when planning the path. This allows paths to be explored efficiently but increases the computational complexity of the problem. Thus, this study proposes a hexagonal grid-based decomposition method and a mixed-integer linear programming (MILP) model to generate a search path that covers all generated nodes in the shortest amount of time with multiple heterogeneous UAVs.

The contributions of our study are twofold as follows. First, to the best of our knowledge, this is the first study to propose a MILP model for CPP based on hexagonal grid decomposition. The model proposed by [44] is a comprehensive model that does not consider the angle between adjacent nodes as a constraint, but it is time-consuming to derive the solution. Therefore, we develop a model that lowers the computational volume by utilizing the angles between adjacent grid cells. Second, we present a coverage path by considering multiple UAVs that operate simultaneously. Most studies on the multi-UAV CPP problem decomposed the target area into subregions and hierarchically interpreted the problem as a multiple CPP problem of a single UAV for each subregion. This method can generate an inefficient path because it cannot plan a path that entirely considers the performance of heterogeneous UAVs.

The remainder of this paper is organized as follows. Section 2 provides a detailed description of the problems and describes the proposed model. The computational experiments are explained and interpreted in Section 3. Finally, conclusions are drawn in Section 4.

2. Problem Description

The objective of this study is to enhance the efficiency and completeness of search operations in maritime SAR systems. Figure 3 shows the overall architecture of the SAR system, which consists of a monitoring system (subsystem 1), a global planning system (subsystem 2), and a local planning system (subsystem 3).

Subsystem 1 identifies the accident and gathers information, including the accident area. Then, subsystem 2 is triggered to create a global SAR plan. Subsystem 2 generates a map for the given search boundary; on this map, a coverage path plan and mission plan are generated for the initial search operation and rescue operation, respectively. These global plans provide a reference for subsystem 3 to create a local SAR plan (i.e., navigation). Subsystem 3 continuously modifies the existing global plans as sensors update environment information.

This study focuses on map generation and coverage path planning for the search operation of subsystem 2. At the level of subsystems 1 and 2, the details regarding the accident and environment are typically unknown or uncertain; thus, the coverage plans of the global strategies for initial search operations must ensure an exhaustive search within a boundary around the accident area. Because there is no room for delay, the completion time of the search is the most crucial consideration in a successful SAR. Therefore, this study aims to produce coverage path plans that minimize the search completion time, which is achieved by the proposed method for map generation.

Through a two-phase process, this study creates a path for multiple UAVs to search a given area. In phase 1, path planning begins by dividing the search area using the hexagonal grid decomposition method. Figure 4 shows the camera footprint in the UAV expressed as a hexagonal grid cell.

Phase 2 transforms the generated hexagonal grid-based map into a graph to obtain input values for path generation. The graph is comprised of nodes centered at each hexagonal grid cell and edges connecting each grid and adjacent grids. The grid-decomposed area is converted to graph form, as shown in Figure 5, and all grids corresponding to the search area are covered. In this study, we set up a MILP model that creates a search path for multiple heterogeneous UAVs that covers all the nodes in the shortest time.

2.1. Phase 1: Hexagonal Grid-Based Decomposition

The search area in this study is divided into hexagonal grid cells, assuming each UAV moves through the grid’s center. Hexagonal grid-based decomposition is used accordingly from the point where the UAV departs. The size of the hexagon is the image size (a) as determined by the camera quality of the UAV (see Figure 6). Covered and uncovered grids are defined in the divided search area and used as input values to create the search path in phase 2.

2.2. Phase 2: Multi-UAV Coverage Path Planning

To formulate a MILP model for CPP based on hexagonal grid decomposition, several assumptions are made. All the nodes in the search area that belong to set $N^{-}$ must be covered. All nodes are covered only once by a UAV, and thus a collision-free path between UAVs is always guaranteed. Because the search operation is terminated once all the hexagonal nodes are covered, a dummy node $i^{-1}$ is introduced as an ending node to indicate the completion of the search. $N^{+}$ is a set of searching nodes including the dummy ending node $i^{-1}$. The starting node $i^0\left(u\right)$ of the UAVs is outside the search area, and it is not in $N^{-}$ or $N^{+}$. The UAVs must move from the current node to adjacent nodes, which indicate the physically neighboring nodes along the sides $e=1,2,\cdots ,6$ of the hexagon shown in Figure 7. The adjacent nodes and dummy ending node belong to $N^{ADJ}\left(i\right)$ for node i. The dummy ending node only has itself as an adjacent node, i.e. $N^{ADJ}\left(i^{-1}\right)=\left\{i^{-1}\right\}$. To move to an adjacent node, cruising and rotational movements are taken into account. In detail, if a UAV passes the nodes $i\to j\to k$ in sequence, the UAV straightly cruises for the paths $i\to j$ and $j\to k$, and rotates as the angle between the vectors $〈i,j〉$ and $〈j,k〉$ at node j. Thus, the travel times of UAVs can be calculated as the sum of cruising time with speed $v_u$ and rotational time with speed $w_u$. Exceptionally, rotational time at the starting node is not considered.

A series consisting of the decision variable $x_{uijes}$ can represent the path of the UAVs. Figure 7 shows an illustrative example to aid in understanding of the decision variable $x_{uijes}$. As can be seen in Figure 7, the exit direction at node i is toward side $e=3$ which connects to node j while the initial heading direction is toward side $e=2$. In this situation, the decision variable $x_{uijes}$ is 1 at $e=2$, and the rotational time at node i is proportional to an angle between the vectors $l\left(e=2\right)$ and $〈i,j〉$. Furthermore, for the travel between the starting node and first visited node, the decision variable $x_{uijes}$ is 1 at $e=7$. Additionally, for the virtual travel between the dummy ending nodes, the decision variable $x_{uijes}$ is 1 at $e=8$.

Mixed-Integer Linear Programming Model

In this section, we describe the MILP model for CPP based on hexagonal grid decomposition.

$$min{t}^{max}$$

(1)

$$s.t.\sum _{i\in N^{+}}\sum _{j\in N^{+}}\sum _{e\in E}{x}_{uijes}=1,foru\in U,s\in S$$

(2)

$$\sum _{u\in U}\sum _{j\in N^{+}}\sum _{e\in E}\sum _{s\in S}{x}_{uijes}=1,fori\in N^{-}$$

(3)

$$x_{uijes}=0,foru\in U,i\in N^{+},j\in N^{+}\backslash \left\{N^{ADJ}\left(i\right)\right\},e\in E,s\in S$$

(4)

$$\sum _{i\in N^{+}}\sum _{e\in E}{x}_{uijes}-\sum _{i\in N^{+}}\sum _{e\in E}{x}_{uije\left(s-1\right)}=0,foru\in U,j\in N^{+},s\in S\backslash \left\{1\right\}$$

(5)

$$\begin{array}{c}\hfill \sum _{i\in N^{ADJ}\left(i^0\left(u\right)\right)}{y}_{ui}=1,foru\in U\end{array}$$

(6)

$$\begin{array}{c}\hfill \sum _{j\in N^{+}}{x}_{uijes}-y_{ui}=0,foru\in U,i\in N^{+},e=7,s=1\end{array}$$

(7)

$$\begin{array}{c}\hfill x_{uijes}=0,foru\in U,i=i^{-1},j=i^{-1},e\in E\backslash \left\{8\right\},s\in S\end{array}$$

(8)

$$\begin{array}{c}\hfill \sum _{\begin{array}{c}k\in N^{ADJ}\left(j\right)\\ e=e\left(〈i,j〉\right)\end{array}}{x}_{uijes}-\sum _{e\in E}{x}_{uije\left(s-1\right)}=0,foru\in U,i\in N^{-},j\in N^{ADJ}\left(i\right)\backslash \left\{i^{-1}\right\},s\in S\backslash \left\{1\right\}\end{array}$$

(9)

$$\begin{array}{c}\hfill t_u^0=\sum _{\begin{array}{c}i=i^0\left(u\right)\\ j\in N^{ADJ}\left(i^0\left(u\right)\right)\end{array}}\frac{d_{ij}}{v_u}\times y_{ui}+\sum _{i\in N^{ADJ}\left(i^0\left(u\right)\right)}\sum _{\begin{array}{c}j\in N^{ADJ}\left(i\right)\\ e=7\\ s=1\end{array}}\frac{\left|r\left(\overrightarrow{i^0\left(u\right)i},\overrightarrow{ij}\right)\right|}{w_u}\times x_{uijes},foru\in U\end{array}$$

(10)

$$\begin{array}{c}\hfill t_u=t_u^0+\sum _{i\in N^{+}}\sum _{j\in N^{ADJ}\left(i\right)}\sum _{s\in S}\left[\sum _{e\in E}\frac{d_{ij}}{v_u}+\sum _{e\in E\backslash \left\{7,8\right\}}\frac{\left|r\left(l\left(e\right),\overrightarrow{ij}\right)\right|}{w_u}\right]\times x_{uijes},foru\in U\end{array}$$

(11)

$$\begin{array}{c}\hfill t_u\le p_u,foru\in U\end{array}$$

(12)

$$\begin{array}{c}\hfill t^{max}\ge t_u,foru\in U\end{array}$$

(13)

$$\begin{array}{c}\hfill t_u^0\ge 0,t_u\ge 0,t^{max}\ge 0\end{array}$$

(14)

$$\begin{array}{c}\hfill x_{uijes},y_{ui}\in \left\{0,1\right\}\end{array}$$

(15)

The objective function (1) is the completion time that is the longest among the travel times of all the UAVs, and the function is minimized. Constraint (2) ensures that each movement step for each UAV has only one path and state. Constraint (3) forces each node in the defined search area to be covered only once, preventing duplicate paths from being reused by different UAVs. Constraint (4) limits movement to nodes that are adjacent to one another. Constraint (5) states that if the path of each UAV is $i\to j$, $j\to k$ becomes the next path in a sequential manner. Constraint (6) defines a movement path from the starting node to one of the nearby nodes. Constraints (7)–(9) are relevant to the definition of the index e in the decision variable $x_{uijes}$. Constraint (7) links the path-related decision variable $y_{ui}$ from the starting node to the search area and the path-related decision variable $x_{uijes}$ within the search area. If the path from the starting node is $i^0\to i$, the first movement step ($s=1$) in the search area sequentially leads to $i\to j$, satisfying the definition of ($e=7$) as the first movement of index e. Likewise, Constraint (8) satisfies the definition of ($e=8$) as the ending movement of index e. Constraint (9) allows the index e to be the index $e\left(〈i,j〉\right)$ of the side facing from the previous path, as defined by the decision variable $x_{uijes}$. Constraint (10) specifies the time required for each UAV to travel from the starting node to the search area, while Constraint (11) defines the time required for each UAV to complete the search. The required time is the sum of the travel time between nodes in a linear distance and the rotation time of the UAV according to the angle formed by the continuous path. Constraint (12) ensures that each search completion time of a UAV is always less than or equal to the maximum duration time. Constraint (13) defines the time needed to complete a search after all nodes have been visited. Finally, Constraints (14)–(15) define the non-negative real and binary variable types of the decision variables.

3. Numerical Experiments

3.1. Experiment Design

We conducted numerical experiments to verify and analyze the methodology presented in this study. The hexagonal grid map generated in phase 1 is defined as an input value of the mathematical model presented in Section 2.2. Figure 8 is one of the target search areas used in the experiments, and the grid cells marked in green are used as the input value.

As mentioned in Section 1, the hexagonal grid decomposition method has a direct impact on the computational complexity of the mathematical model for path generation as well as the target value of the optimal solution. Hence, the time required to calculate the optimal solution in response to changes in the size of the search area is analyzed. In addition, the optimal UAV search path generated in response to changes in the shape of the search area is presented.

All UAVs maintain a constant altitude during a mission, and thus the camera footprint of all UAVs in this experiment is the same as a circle with a radius of 20 m. The performance of each UAV is shown in

Table 1. UAV performance.

UAV	${\mathit{v}}_{\mathit{u}}$ (m/s)	${\mathit{w}}_{\mathit{u}}$ (radian/s)	${\mathit{p}}_{\mathit{u}}$ (s)
UAV 1 (Quad-rotor)	4	0.5	1800
UAV 2 (Quad-rotor)	4	1	1200
UAV 3 (Quad-rotor)	4	0.5	1500

. The experiment was conducted using an Intel Core i7-7700HQ 2.80 GHz, 16.0 GB RAM, Windows 10 operating system environment, and a Gurobi 8.0.1 solver (https://www.gurobi.com (accessed on 17 September 2021)) was used to solve the MILP model. The Python programming language was used in all experiments.

3.2. Time Needed to Calculate the Optimal Solution in Response to Changes in the Size of the Search Area

The multi-UAV CPP model proposed in this study has high computational complexity. The complexity in finding an optimal path increases in direct proportion to the size of the search area. This section discusses the time needed for the mathematical model presented in Section 2.2 to calculate the optimal solution using the commercial solver (Gurobi 8.0.1). Twenty-seven instances with different search area sizes were used in this experiment while retaining similar shapes. Figure 9 shows the instances with the smallest and largest number of nodes, respectively.

Figure 10 shows the relationship between the time required to calculate the optimal solution using the commercial solver and the number of search nodes. When the number of search nodes was less than 27, an optimal solution can be found within a minute based on the instance used in this experiment. When the number of search nodes reaches 40, the commercial solver takes more than an hour to find the optimal solution. Figure 10 shows that as the size of the search area grows, the computational time grows exponentially, indicating an exponential growth in the computational complexity.

3.3. Search Path According to Changes in Search Area Shape

The optimal search path derived from the mathematical model presented in Section 2 is discussed in this section. As shown in Figure 11, we used various search areas to test the proposed model in this experiment. Figure 11a–c are instances used in Section 3.2, and the shapes of the search areas are similar but the sizes are different. Figure 11d–f are instances that have different irregularly shaped search areas.

Despite the similar shapes of the search areas, the search paths in Figure 12a–c appear in a completely different manner. As shown in Figure 12d–f, the mathematical model presented in this study can generate optimal routes for several types of search areas, even if they are not standardized.

3.4. Search Path According to the Change in the Shape of Grid Cells

In this section, the effectiveness achieved by the hexagonal grid-based area decomposition is validated. For the same instances in Figure 11, we tested the proposed MILP model on a map with square grid cells and on a map with hexagonal grid cells and compared the results. There is no structural difference in the MILP model except for considering four adjacent sides instead of six adjacent sides in set E. Figure 13 shows the optimal search paths on the square grid map derived by the MILP model.

Compared to the results in Figure 13, on the hexagonal grid map in Figure 12, the search times of the served UAVs are reduced for all test instances. The search completion times on the hexagonal grid map are on average 10.59% lower than those on the square grid map. For instance (a), the worst case, the reduction reaches 23.39%. Furthermore, on the hexagonal grid map, the average difference between the search times of the served UAVs is 1.30 seconds, compared with an average difference of 4.02 seconds on the square grid map. Thus, the reduction is 67.7%. All these results are intuitively reasonable because the hexagonal grid decomposition allows more movable nodes for the served UAVs and thus creates relatively flexible paths.

3.5. Search Path According to the Change in the Number of Served UAVs

In this section, the performance of the proposed model with three UAVs is verified. Similar to the previous experiment, for the same instances in Figure 11, we tested the proposed MILP model with three UAVs. Figure 14 shows the optimal three-UAV search paths derived by the MILP model.

Compared to the results for two UAVs in Figure 12, the reduction in search completion times when using three UAVs is on average 28.4%. The difference between the maximum and minimum search times of the UAVs is 2.75 seconds on average. These results are encouraging in that the proposed method efficiently distributes UAV operations without any collisions despite the increasing number of served UAVs.

4. Conclusions

Marine accidents in Korea result in casualties and ongoing environmental and property damage. Because survivors can be exposed to increasingly dangerous conditions and uncertainty as time passes, it is critical for first responders to respond quickly. If an accident occurs, the search for survivors is fraught with difficulties resulting from various uncontrollable factors in the marine environment. Hence, public institutions such as the Ministry of Oceans and Fisheries, the Korea Coast Guard, and the Navy are making efforts to address this issue. In this study, multiple heterogeneous UAVs were utilized to search for casualties in a short period of time to protect lives from catastrophic marine accidents and to prevent the spread of damage. We proposed a hexagonal grid decomposition to represent the search area as a hexagonal grid map, and developed the MILP model that reduces the time required for all UAVs to complete a search using the graph for this grid map as an input. The proposed model in the numerical experiment was verified by a simulation analysis based on the performance of an actual UAV. This study presents both the optimization technique’s calculation time as a function of the search area size and the multiple UAV paths for each of the various search areas. The experimental results of this study provided the following two insights. First, the proposed model can be comprehensively operated even with changes such as the shape or size of the search area and the diversity or number of UAVs. Second, because the hexagonal grid cell has more adjacent grids than the square grid cell, the proposed model has a larger solution space and can derive a more efficient path.

This study confirmed the relationship between the number of search nodes and computational complexity based on a small sample size; subsequent experiments on more diverse environments could provide more comprehensive insight. In future works, related algorithms (dynamic programming, reinforcement learning, etc.) can be studied to solve larger-scale problems with reasonable computational time. In addition, this study which focuses on deterministic global planning can be extended in consideration of the factors that change over time such as weather conditions and ocean currents. Furthermore, CPP for search operations affects mission planning for rescue operations and local planning systems in maritime SAR systems, and future research can be conducted in an integrated manner. Lastly, not only can this study be used to aid in the decision-making of a quick and accurate initial response to marine accidents, but we also expect the results from this study to lead to a wider application in land-based emergencies.