Computationally-Efficient Distributed Algorithms of Navigation of Teams of Autonomous UAVs for 3D Coverage and Flocking

: This paper proposes novel distributed control methods to address coverage and ﬂocking problems in three-dimensional (3D) environments using multiple unmanned aerial vehicles (UAVs). Two classes of coverage problems are considered in this work, namely barrier and sweep problems. Additionally, the approach is also applied to general 3D ﬂocking problems for advanced swarm behavior. The proposed control strategies adopt a region-based control approach based on Voronoi partitions to ensure collision-free self-deployment and coordinated movement of all vehicles within a 3D region. It provides robustness for the multi-vehicle system against vehicles’ failure. It is also computationally-efﬁcient to ensure scalability, and it handles obstacle avoidance on a higher level to avoid conﬂicts in control with the inter-vehicle collision avoidance objective. The problem formulation is rather general considering mobile robots navigating in 3D spaces, which makes the proposed approach applicable to different UAV types and autonomous underwater vehicles (AUVs). However, implementation details have also been shown considering quadrotor-type UAVs for an example application in precision agriculture. Validation of the proposed methods have been performed using several simulations considering different simulation platforms such as MATLAB and Gazebo. Software-in-the-loop simulations were carried out to asses the real-time computational performance of the methods showing the actual implementation with quadrotors using C++ and the Robot Operating System (ROS) framework. Good results were obtained validating the performance of the suggested methods for coverage and ﬂocking scenarios in 3D using systems with different sizes up to 100 vehicles. Some scenarios considering obstacle avoidance and robustness against vehicles’ failure were also used.


Introduction
In recent years, there has been an increasing interest in mobile wireless sensor networks (MWSNs), where a number of networked autonomous vehicles can be deployed in different environments to achieve sensing tasks. Advances in communication made MWSNs more appealing where vehicles (sensors) can share information to perform cooperative monitoring, sensing, detection and exploration.
The problem of controlling multi-vehicle systems to achieve a global objective is referred to as cooperative control. There exist different subset problems related to this area depending on the common goal to be achieved by the vehicles. Coverage control deals with problems where a certain region of interest needs to be surveyed whether using a single vehicle or multiple (i.e., using MWSNs). Tackling these problems with multi-vehicle systems has given rise to new challenges to traditional cooperative control in the field of coverage control [1][2][3][4][5][6][7][8]. Unmanned aerial vehicles (UAVs) have become a popular choice in this area to form MWSNs especially in places inaccessible by ground vehicles. In general, multi-UAV systems have been emerging in various applications such as precision agriculture [9][10][11][12][13], aerial manipulation and transportation [14][15][16][17][18], search and rescue [14,19,20], firefighting [21], surveillance and monitoring [22,23], mapping and exploration [24][25][26], etc.
Coverage control problems can be classified as either static or dynamic. Another classification is based on [27] where coverage problems are categorized into Blanket coverage, Barrier coverage and Sweeping coverage, which are defined as follows: • Blanket coverage is forming a static arrangement to maximize the detection rate of events through an area of interest. • Barrier coverage is a static formation over some region (i.e., a barrier) to minimize intrusions or maximizing detection of objects going through it. • Sweeping coverage is the formation of dynamic arrangements moving across a region of interest for maximal detection/exploration along the whole region.
Clearly, blanket and barrier coverage problems belong to the static class while sweeping is a dynamic coverage problem.
Static coverage can be applied in different applications such as surveillance, real-time monitoring of crops, live stock and pollution, intruders detection in the security domain, etc. In these scenarios, the combined sensing field of view of all sensors is sufficient to monitor the regions of interest with a static configuration that the MWSN converges to. On the other hand, dynamic coverage is based on applications where the sensors need to remain mobile in order to survey the regions of interest. Examples of these applications include [28], but not limited to, autonomous search & rescue, monitoring of dynamically-changing environments, real-time monitoring of large areas with limited number of vehicles, visual inspection using multiple vehicles with limited combined field-of-view (FOV), firefighting in forests using a group of vehicles with limited spray area, autonomous irrigation with limited spray area, etc.
Flocking is another class of problems where cooperative control is applied, which is a subset of formation control, to deal with movement of large number of vehicles (swarms). In general, controlling a large number of vehicles to a achieve a common motion objective is motivated by collective behavior of animals such as bird flocks and fish schools. In this case, vehicles can make motion decisions by following a set of rules contributing towards the collective behavior (i.e., flocking). These rules are: flock centering (cohesion), collision avoidance (separation) and velocity matching (alignment) according to Reynold's model [29,30].
The main contribution of this work is to develop novel distributed control strategies addressing cooperative control problems to achieve 3D coverage and flocking behavior. Specifically, the proposed methods address the 3D barrier and sweep coverage problems motivated by some of the ideas in [1,6,31]. In a 3D environment, a barrier can be defined as a static arrangement of sensors with overlapping sensing zones [32] forming a surface or a 3D region. The suggested control strategies rely on estimated centroidal Voronoi configurations over virtual regions generated by the sensors locations in a distributed manner depending only on shared information from neighbor vehicles. Furthermore, the suggested methods are also extended to address flocking problems in 3D spaces where the movement of the 3D virtual region is decided collectively by the group to achieve the common motion objective.
The designed control laws require relative distances with neighbor vehicles to be shared over communication channels which makes the overall problem related to Networked Control Systems (NCSs) [33][34][35][36][37][38][39]. The current work assumes that such information is available to the control system. However, several challenges related to communication channels imperfections needs to be considered when evaluating the performance of the overall networked control system. Thus, it is important to design a reliable communication system to allow for advanced swarming behaviors of robotic systems [40].
Overall, the vehicles' collective motion becomes constrained within a specific dynamic region under the application of the suggested control methods similar to region-based shape control methods [41]. The dynamic region can be selected as a collection of non-overlapping polygons in 3D. Furthermore, one can control the dynamics of the dynamic region to generate 3D sweeping behavior which is the key idea used in the developed 3D sweeping coverage strategy. This is also considered to handle obstacle avoidance where vehicles can collaboratively control the dynamics of the virtual barrier and even apply deformations to its shape in real-time, which is then communicated through the networked multi-vehicle system. Furthermore, bounded control laws are proposed, which is important in practice to satisfy limits on the vehicles' velocities and accelerations. The main advantages of the suggested approaches can be highlighted as follows: • collision avoidance among vehicles and connectivity is ensured by the adopted Voronoi-based approach • the approach is highly scalable and robust against vehicles' failure • obstacle avoidance can be managed in a decomposed and distributed manner A general 3D kinematic model is adopted in the design which is applicable to different UAV types and autonomous underwater vehicles (AUVs). A 6 degrees-of-freedom (DOF) dynamical model for quadrotors is further considered to show a possible way of implementation with low-level control design. Several simulations were carried out to validate the performance of the suggested methods in addition to showing its scalability and robustness. Moreover, software-in-the-loop (SITL) simulations were also performed using the Gazebo robotic simulator based on the quadrotor full dynamical model to evaluate the computational complexity of the implemented algorithms with particular interest in applications related to precision agriculture.
The organization of this paper is as follows. A summary of related works is given in Section 2 highlighting some of the drawbacks that motivated the proposed design. Section 3 introduces some essential concepts related to graph theory, locational optimization and Voronoi Partitions which are used in our control strategy, and the tackled 3D coverage problems are defined in Section 4. After that, distributed barrier and sweeping coverage control strategies are proposed in Section 5 considering a general 3D kinematic model. These approaches are validated through several simulation cases in Section 6. Furthermore, a generalization of the suggested methods is presented in Section 7 to address 3D flocking problems along with validation. Further implementation details considering quadrotors dynamics with low-level control design are presented in Section 8, which is evaluated using software-in-the-loop simulations. Finally, this work is concluded in Section 9 with a suggestion for a potential direction of future work.
An example algorithm of static deployment of cameras (i.e., non-mobile sensors) to achieve coverage with optimal resource allocation was proposed in [44]. In [43], modified Lloyd-like algorithms were developed to address the 2D static deployment of mobile sensors while considering optimizing the overall power consumption. The tradeoff between coverage performance and energy consumption of mobile heteregoneous wireless sensor networks was also investigated in [45].
Optimal placement of fixed and mobile networked sensors have been tackled in some works using optimization-based techniques. For example, the work in [46] adopted a particle swarm optimization-based approach to address the static deployment of cameras over 2D and 3D monitoring spaces for optimal visual monitoring. On the other hand, 2D static deployment of mobile wireless sensors was addressed in [47] using the Multiobjective Immune Algorithm (MIA). The problem formulation targeted the maximization of the coverage area while minimizing energy consumption due to movement. The authors have also considered the connectivity perseverance through limited movement within specified communication ranges. However, collision avoidance between mobile sensors and obstacles avoidance were not considered.
Another class of approaches addressing coverage problems of MWSNs is based on artifical potential field, which makes it easier to include collision avoidance rules as potential functions. In this case, control laws are designed based on gradient descent algorithms. For example, a decentralized approach was developed in [48] addressing the problem of deploying mobile sensors with limited sensing capeabilities to achieve the best possible coverage over an n-dimensional Euclideon space. The mobility of sensors was modeled as single integrator. The devleoped control law relies on forcing the inter-agent distances to reach a desired distance while avoiding collisions with other vehicles. Obstacle avoidance was not considered in this work. On the other hand, another method was presented in [49] using artificial potential fields considering both forces between the mobile sensors and obstacles in the environment to achieve both static deployment and collision avoidance. The method was validated considering only 2D environments.
A similar class of methods existing in the literature is based on Voronoi diagrams. The formulation of of such methods relies also on gradient descent algorithms adopting potential functions that can encode optimal coverage in a distributed manner, which was well developed in [1,31,59] with rigorous mathematical proofs. In these works, Lloyd's algorithm was adopted to develop distributed control methods for motion coordination of MWSNs to achieve optimal coverage. Only 2D coverage problems were considered, and obstacle avoidance was not considered. The approaches were distributed in nature where each vehicle can compute its own Voronoi cell based on information shared from its Voronoi neighbours. Similarly, an adaptive 2D approach was proposed in [5] based on Voronoi diagrams. However, the authors of this work suggested a decentralized control law design, where each vehicle can approximate the centroid of its Voronoi cell based only on its sensory measurements rather than relying on information received from its Voronoi neighbors. This can be more practical in cases where Voronoi neighbors are at distances larger than the vehicle's communication range. This approach was also considered in [50]; however, a different weighting function, which defines the sensing performance over a planar region, was considered.
Many of the above approaches consider planar circular sensing models. Some works have considered different sensing models such as [51,52]. In these works, non-uniform unequally scaled 2D sensing footprints were considered. Additionally, Voronoi partitioning was performed only with respect to the sensed space, rather than the whole region of interest as in the previous methods.
Some research works have also considered additional requirements in the coverage problem formulation such as graph connectivity perseverance and information routing. For example, the 2D static coverage problem tackled in [53] considered that the mobile sensors are required to send the collected sensory information to a set of destinations over the network while performing the coverage task. Their solution was based on a constrained optimization problem considering Voronoi partitioning as part of the solution.
The above methods are mainly developed to address static coverage problems of MWSNs, where the sensors need to converge to a static arrangement with optimal coverage. There also exist methods addressing dynamic coverage problems such as [28,[60][61][62][63][64][65][66]. In [28], a multi-agent system was considered with a 2D nonholonomic agent model. The approach relied on a leader-follower paradigm, where the vehicles try to maintain desired distances from a position targeted by the leader. This can be more challenging in 3D and for systems with large number of vehicles. Another approaches considering 2D nonhonolomic agents were developed in [61,62] proposing control laws based on a special form of Lyapunov-like functions to achieve satisfactory coverage level and collision avoidance.
Motivated by the Voronoi-based static coverage methods, the work in [63] suggested a planar approach based on modified Voronoi partitions with time-varying density functions, which describe the coverage performance, to address the planar dynamic coverage problem. However, obstacle avoidance was not considered in this work. A similar approach was also proposed in [63] where the authors suggested a Bayesian prediction method to estimate the unknown time-varying density function.
Overall, many of the existing static and dynamic coverage control approaches consider only two-dimensional sensing fields, and the literature lacks a proper analysis of sensor networks deployed in three-dimensional (3D) sensing fields [67]. Even those methods proposed for multi-UAV and multi-AUV systems assume that the vehicles are moving at a fixed altitude/depth with a planar sensing footprint without utilizing the full capabilities of such vehicles. It is hence more motivating to work towards addressing 3D coverage problems exploiting the rich geometric properties of 3D MWSNs [67]. Some efforts have been made in that area such as [32,[67][68][69][70][71][72][73]. Thus, one of the main motivations of this work is to contribute towards the state-of-the-art 3D static and coverage problems using MWSNs.
On the other hand, the flocking problem has also attracted great interest. Many researchers have tackled the flocking control problem where the global group objective is to move the whole group towards some goal region or to follow a certain motion pattern (for example, see [29,[74][75][76][77][78][79][80][81][82][83]). Most of state-of-the-art flocking control methods rely completely on artificial potential field to handle local interactions between vehicles, which have some disadvantages such as being prone to getting stuck at a local minimum especially due to conflicts between forces required to achieve moving towards goal, collision avoidance and obstacle avoidance control objectives. To address these limitations, the suggested control strategy separates the control components for these objectives on different layers as will be shown later.

Preliminaries
The proposed methods in paper relies on concepts from graph theory, locational optimization and Voronoi partitions. A summary of these concepts is provided in this section based on [1,29,31,59]. Note that when considering a multi-UAV system as a mobile wireless sensor network, UAVs may interchangeably referred to throughout the paper as sensors, nodes, agents or vehicles.

Graph Theory
A multi-UAV sensor network consisting of n UAVs can generally be characterized using a set of nodes/vertices U = {1, 2, · · · , n} and a set of edges (paired vertices) E ⊆ {(i, j) : i, j ∈ U , j = i}. Note that i and j will be used in the rest of the paper to refer to the index of a node/vertex within a set (i.e., index of a UAV within the group). Each vertex corresponds to a single UAV/sensor, and edges represent interaction between UAVs which are within communication or detection range from each others. The overall network topology is then described using a graph G = (U , E ) which can be directed or undirected. In an undirected graph, an edge exists from vertex i to vertex j if and only if an edge exists from j to i (i.e., (i, j) ∈ E ↔ (j, i) ∈ E ). Otherwise, the graph is called directed. Generally, homogeneous multi-UAV systems can be described using undirected graphs since all UAVs have same communication and sensing capabilities. Moreover, a path between two vertices i 0 and i k is defined as a sequence of vertices {i 0 , i 1 , · · · , i k } ⊂ U where an edge exists between each subsequent vertices in the sequence such that (i l , i l+1 ) ∈ E , ∀l ∈ {0, 1, · · · , k − 1}. If every pair of vertices in U is connected by a path, the graph G(U , E ) is then called connected. Clearly, a crucial part for MWSNs is to maintain network connectivity all the time.
Furthermore, define a neighborhood N i around a vertex i as the set of all vertices which have edges with vertex i such that: For a homogeneous system, let r > 0 denote the communication range for all UAVs. Hence, all UAVs within a spherical region of radius r around UAV i belong to its neighborhood such that where p i ∈ IR 3 is the position of UAV i, and || · || is the Euclidean norm in IR 3 . Similarly, all norms in the rest of the paper are Euclidean norms in IR 3 .

Locational Optimization
Deployment of mobile sensors in an environment to achieve optimal sensor coverage is regarded as a multicenter problem from locational optimization (i.e., spatial resourceallocation problem). A brief description about some of the facts related to this class of problems is summarized next based on [1,59].
Consider a bounded region of interest Q, including its interior, defined in a space IR d with dimension d. A partition of Q consists of a group of n non-overlapping polytopes W = {W 1 , · · · , W n } such that W 1 ∪ · · · ∪ W n = Q. Furthermore, let φ : Q → IR + , (IR + = {a ∈ IR : a > 0}) be defined as a distribution density function representing a measure of information or the likelihood of an event to take place over Q. The sensing performance of a sensor located at some position p i as seen from any point q ∈ Q depends mostly on the distance ||q − p i ||. Clearly, as this distance increases, the sensing performance degrades. Hence, one can describe the sensing performance at location q of the sensor p i using a non-increasing piecewise continuously differentiable function f (||q − p i ||) : IR + → IR. Thus, the larger the value of f , the better the sensing performance at q is.
Using the above definitions, one can define a multicenter cost function characterizing the average coverage provided by a set of n sensors at p 1 , · · · , p n over an point in Q as follows: The above function provides a measure of the sensing performance expected value provided by all sensors at any point q ∈ Q [31]. Now, in order to find the optimal placement for all sensors, an optimization problem needs to be solved to maximize the value of H(p 1 , · · · , p n ). Remark 1. Note that there are slightly different definitions for f in the references [1,31,59] where it can be either considered as a representation of sensing degradation or sensing performance over Q (as considered here). This does not affect the overall analysis done here except that the considered optimization problem will either be minimization (of sensing degradation) or maximization (of sensing performance).

Voronoi Partitions (Tessellation)
This subsection highlights some key points about Voronoi partitions needed for our problem formulation. A Voronoi partition/tessellation is the subdivision of a space into a number of regions generated by a set of points (see Figure 1 for a 2D example). Consider that we have n sensors located at fixed locations P = {p 1 , · · · , p n } ∈ Q n . A voronoi partition of Q consists of a set of disjoint Voronoi regions/cells V (P) = {V 1 , · · · , V n } generated by these sensors where and It has been established that this Voronoi partition is the optimal partition of Q among all other partitions [59]. For any sensor located at a position p i , its Voronoi neighbors N V,i ⊂ P are defined as the sensors corresponding to adjacent Voronoi cells such that: Considering the above definition, one can rewrite (3) as: By taking the partial derivative of (6) with respect to p i , the following is obtained: where it is assumed that f does not have any discontinuities. Furthermore, considering f (x) = −x 2 , the multicenter cost function in (6) becomes: where J V i ,p i is the polar moment of inertia of V i about p i . Consequently, (8) reduces to: where M V i ∈ IR + and C V i ∈ IR 3 are the mass and center of mass (centroid) of the corresponding Voronoi partition V i with respect to the density function φ(q). It is clear from (9) that the critical points of H V (P, V (P)) are the configurations P ∈ Q n where p i = C V i ∀i, which are referred to as centroidal Voronoi configurations.

3D Coverage Problems
Consider a 3D bounded region of interest S ⊂ IR 3 . A multi-vehicle system can perform coverage tasks over S where coverage objective may vary according to the problem in hand. Definitions of the considered barrier and sweep coverage problems in 3D are defined next. Problem 1. (3D Barrier Coverage) Deploy a network of vehicles/sensors to form a static arrangement over some region B ⊂ S (i.e., a barrier) maximizing the sensing performance of the overall network to detect any intruder going through the barrier.
A special case of the above problem is when deploying the sensors over a planar region within S defined by B = {(x, y, z) : σ(x, y, z) = 0} where σ(x, y, z) is the plane's equation.
Note that for this problem to be solvable, the number of vehicles/sensors needed depends on the size of S and the sensing range of all vehicles (assuming a homogeneous system).

Problem 2.
(3D Sweep Coverage) Consider a group of vehicles whose overall sensing range is not large enough to achieve complete coverage over S. It is required to survey the region S by moving the whole group across S as a dynamic formation over some region F (t) ⊂ S.
The coverage task can be considered completed once the whole region S is scanned/surveyed. Alternatively, the coverage task can be done repetitively over a specified period of time (i.e., continuously survey S for several times).
Note that F (t) can be of any 3D shape in which its motion along S following a certain pattern can achieve complete coverage. A special case is when F (t) is a plane which is considered in this work. In this case, F (t) will referred to as the sweeping plane. The dynamics of the sweeping planeḞ (t) can be determined in a way to achieve complete coverage over S. It is also assumed that F (t) can change size and shape over time which can be utilized for other motion objectives such as obstacle avoidance as will be shown later.

Problem Formulation
The aim of this work is to develop distributed control laws for multi-UAV systems to address Problems 1 and 2. We consider a system of n homogeneous vehicles (UAVs/AUVs) with a single integrator motion model given by: where p i ∈ IR 3 is the i-th vehicle position defined in some inertial frame {I}, and u i ∈ IR 3 is its control input (velocity) where All vehicles can sense events in the environment within a sensing range r s > 0. Furthermore, any vehicle can exchange information with nearby vehicles within some communication range r c > 0. Obviously, it is assumed that r c > r s so that it is possible to design control laws which can maintain the connectivity of the network with minimal sensors overlapping. It is assumed that the vehicles are homogeous in terms of sensing and communication ranges. However, they can be of different types, sizes and/or weights.

Distributed Coverage Control Strategies
The proposed control schemes to address Problems 1 and 2 are based on a selfdeployment method for the multi-vehicle system over a planar region B ⊂ S which is static for barrier coverage problems and dynamic for sweep coverage problems. We consider a set of l vertices E to describe the boundary of B as a polygon such that E = {e 1 , e 2 , · · · , e l } ⊂ S; clearly, l > 2. Lloyd's algorithm is adopted in the designed controllers to guide the vehicles to reach the instantaneous centroids of their associated Voronoi regions over B (i.e., reaching the centroidal Voronoi configuration). Once this is reached, an optimal coverage over B is achieved. Furthermore, for sweeping problems, the designed dynamics of B will achieve coverage over S.

Online Computation of Centroidal Voronoi Configurations
The developed control law requires each vehicle to be able to compute the centroid of its Voronoi region in a distributed fashion based only on information exchanged with vehicles within its neighbourhood. We extend the approach proposed in [1] to compute Voronoi cells for planar regions in 3D in a distributed fashion.
To simplify the mathematical development, a new 3D coordinate frame {B} attached to B is needed. The origin of {B} can be selected to be one of the barrier vertices defined as O B = e 1 . Furthermore, the axes of {B} are defined using the orthonormal basis { a 1 , a 2 , a 3 } where: where h( w 1 , w 2 ) = w 2 − ( w 1 · w 2 ) w 1 is a mapping function which gives an orthogonal vector to w 1 ∈ IR 3 directed towards w 2 ∈ IR 3 .
All computations needed to find Voronoi centroids are carried out in the {B} frame through a transformation between {B} and the inertial frame {I}. Let v I ∈ IR 3 be a vector defined in the {I} frame. This vector can be transformed to the {B} coordinate frame using a transformation matrix T B I as follows: where T B I is a 4 × 4 affine transformation matrix given by: Note that we represent the transformation using an augmented matrix to consider both rotation and translation in a single matrix multiplication. Furthermore, a 1 , a 2 , a 3 , and O B are column vectors.
From now on, vectors represented in the inertial frame will be represented without the I superscript for simplicity. Given a UAV at position p i , it is required to compute instantaneous Voronoi centroid C V i of its projection onto B. First, the following assumption is made.
The proposed approach can now be described in these steps: S1: Transform the position p i into the {B} frame to obtain i using the following [1]: where These equations are obtained considering Assumption 1 and a constant distribution density function φ(q) = 1. Voronoi cell vertices {v B 1 , v B 2 , · · · , v B m } can be determined based on locations of Voronoi neighbours in a distributed fashion (see Remark 2). (13).

Remark 2.
The vertices of a Voronoi cell V i can be found as the circumcenters of triangles formed byp B i and any two of its Voronoi neighbours. A triangle made by three points p 1 , p 2 and p 3 with an area of A has a circumcenter at [1]: where α ls = p l − p s .

Barrier Coverage Control Design
In order to present our control design, some technical assumptions need to be made as follows.
Assumption 2. The communication graph G(U , E ) remains connected for t ≥ 0.

Assumption 3.
Each UAV is capable of estimating the Voronoi cell centroid of its projected position onto B at any time using only information from UAVs within its communication range (i.e., UAVs in its neighborhood N i ).

Assumption 4. The Voronoi neighbours ofp B
i correspond to UAVs within the neighbourhood N i .

Assumption 5.
The initial configuration of the multi-UAV system satisfies the following condition: , · · · , n}, i = j where a 3 is the barriers normal as defined in (12).
Assumption 2 is made to make sure that updated information about the barrier B is available to all vehicles at any time during the motion. This is essential in cases where B is dynamic. For example, a decision could be made by one of the UAVs to apply changes to the shape of B based on some detected obstacles. Such information needs to be shared among all vehicles so that they can compute their Voronoi regions accordingly. It is possible to ensure the connectivity of G(U , E ) by enforcing constraints on B based on the number of vehicles and communication ranges to ensure that Voronoi neighbors are within the communication range of each others. Assumptions 3 and 4 ensures that UAVs can determine Voronoi centroids in a distributed fashion. Finally, Assumption 5 ensures that all UAVs are initially located at positions with unique projections onto B. Now, the main results of this section can be presented. Consider the following control law u i : IR 3 → IR 3 based on Lloyd's algorithm: where the parameterK i is a diagonal positive definite gain matrix. We also propose a more practical bounded control law u i : IR 3 → Γ which can satisfy bounds u max on the control input such that Γ = {u ∈ IR 3 : ||u|| ≤ u max }. It is given as follows: where K i = diag{k i,x , k i,y , k i,z } is a positive definite diagonal matrix, γ i > 0, and tanh(v) is the hyperbolic tangent function defined element-wise for any vector v ∈ IR 3 . Clear, the bound of this control law depends on the gain matrix K i as follows: Theorem 1. Consider a multi-UAV system with n vehicles whose motion are modelled as (10). Under Assumptions 2-5, the distributed control law (16) (or (17) for bounded inputs) along with the algorithm in S1-S4 solves the 3D barrier coverage problem defined in Problem 1.

Proof.
Let P e be a centroidal Voronoi configuration. We define a Lyapunov canidate function as where H V is defined in (8) with φ(q) = 1, andH = H V (P e , V (P e )) which is constant. Hence, V 1 (P e ) = 0. Furthermore,H ≥ H V (P(t), V (P(t))) ∀P ∈ S n \{P e } since a centroidal Voronoi configuration is optimal for H V among all other configurations (see Proposition 2.14 in [59]). This indicates that V 1 (P(t)) > 0 which makes it a valid Lyapunov function.
The time derivative of (19) can be obtained using (9) and (10) and the control law (16) as: . Hence, the set of centroidal Voronoi configurations P e is locally asymptotically stable, and lim Similarly, for the bounded control law, the time derivative of V 1 (P(t)) can be obtained by substituting (17) into (20) as follows: It is evident from (24) thatV 1 < 0 ∀P ∈ S n \{P e } since the hyperbolic tangent function is an odd function and γ i > 0. This implies that P e is also locally asymptotically stable under the application of the bounded control law (17).
Therefore, the control law (16) (or (17)) guarantees that the vehicles will converge to centroidal Voronoi configurations which maximizes the sensing performance over the barrier B according to Proposition 2.13 in [59]. Furthermore, Assumptions 2-5 ensures that all vehicles can compute the centroids of their Voronoi cells in a distributed fashion at all times with no overlapping following the algorithm in S1-S4. This completes the proof.
Note that the proposed control law ensures that the vehicles will converge to centroidal Voronoi configurations generated by their projections onto B. Thus, all the vehicles will eventually reach B such that lim t→∞ p i =C V i ∈ B even if they are initially deployed at some positions p i (0) / ∈ B. Additionally, the trajectory of each vehicle remains within its Voronoi region which does not intersect with any other regions by definition. This guarantees that vehicles motions are collision-free using the proposed control laws.

Sweep Coverage Control Design
Theorem 1 shows that the proposed control laws can force the vehicles to reach a specified region within the 3D space (i.e., the barrier) and constrain their motion within that region. This is extended in this section to address the sweep coverage control problems. It can be achieved by enforcing vehicles to deploy over some dynamical "virtual" region whose motion is determined by the group.
In general, motion coordination control laws for coverage problems needs to satisfy the following objectives: O1 Avoid collisions with other vehicles while maintaining a certain formation as a group O2 Avoid collisions with obstacles within the environment O3 Achieve optimal coverage of the targeted environment S collaboratively The proposed sweeping algorithm targets these objective on two levels. At a lower level, the vehicles motions are constrained within a dynamic "sweeping "region F (t) ⊂ S, and they maintain an optimal formation over that region for maximal sensing. This achieves objective O1. At a higher level, decisions can be made in real-time collaboratively by the vehicles to decide the dynamics of F (t) (i.e.,Ḟ (t)). Note that it is also possible to apply deformations to F (t) (will be shown in simulations) as long as the deformed region is large enough for the vehicles to distribute over with safe spacing. This provides a good way in addressing objectives O2 and O3. In particular, the trajectory of F (t) will result in sweeping the whole environment S providing optimal coverage. Moreover, obstacle avoidance can be achieved by only changing the dynamics of F (t) rather than having each vehicle reacting independently to obstacles.
For r ∈ F (t), the motion dynamics of the sweeping region can be described as follows: r = g(t, r), g(t, r) ≤ḡ, ġ(t, r) ≤ḡ, where g(t, r) ∈ IR 3 is a desired velocity profile,ḡ,ḡ > 0 are upper bound design parameters to account for dynamic limitations where the following condition must hold: In other words, the speed of F (t) should not be larger than the maximum physical speed that can be achieved by the vehicles. Note that g(t, r) ∈ IR 3 can have any direction. However, for simplicity in achieving sweeping coverage, it assumed that F (t) is a planar region, and it is defined similar to B (see Section 5.1) at different time instants. A simple example is moving F (t) in the direction of its normal (i.e., a 3 ) with a constant sweeping speed g 0 > 0 such that: g(t, r) = g 0 a 3 .
More complex movements can be achieved depending on the considered environment shape and nearby obstacles. One can also adopt a 3D holonomic or non-holonomic model for (25) utilizing the available literature in obstacle for these models. Generally, obstacle avoidance can be achieved either by varying the dynamics in (25) or by dynamically deforming F (t). Simulation cases showing both approaches will be shown later. Note that the distributed behavior of the proposed algorithm is maintained in all these cases since information about g(t, r) are exchanged over the connected network following Assumption 2. At this point, we will leave out the design of (25) to a high-level controller shared among the vehicles while assuming the following: Assumption 6. The sweeping plane F (t) remains all the time within the sensing environment (i.e., F (t) ⊂ S, ∀t > 0), and its movement governed by (25) will completely span the volume of the sensing environment S.
The main results of this section can now be presented.

Theorem 2.
Consider a multi-UAV system of size n where each vehicle's motion model is represented by (10). The control law (16) along with the algorithm in S1-S4 defined for a sweeping region F (t) whose dynamics is governed by (25) solves the 3D sweeping coverage problem defined in Problem 2 under Assumptions 2-6 and the condition (26).

Proof. Let
where u i (t) is given by (16). Now, define a Lyapunov candidate function as follows: , · · · , e C n (t),ė C 1 (t), · · · ,ė C n (t)) where V 1 is defined in (19), and k 1 and k 2 are positive kinematic constants with units of m 2 and m 2 · s 2 to match the physical units of V 1 which is m 4 . The choice (29) guarantees that V 2 (t) > 0 for (e C 1 (t), · · · , e C n (t),ė C 1 (t), · · · ,ė C n (t)) = (0, · · · , 0), and V 2 (0, · · · , 0) = 0 is also true. The time derivative of (29) is obtained using (10) and (16) as follows: Recall thatV 1 < 0 which was established in (24). Therefore, the time derivative of V 2 is negative outside the compact set B Γ = B Γ 1 ∪ · · · ∪ B Γ n where B Γ i = {e C i ∈ IR 3 ,ė C i ∈ IR 3 : e C i > Γ i,1 , ė C i > Γ i,2 } and Γ i,1 and Γ i,2 are given by: which represent closed balls with radius Γ i,1 and Γ i,2 respectively. Thus, starting from any initial condition outside the set B Γ , the errors will converge to the closed set B Γ in finite time and stay there forever. That is, which means that the errors will be uniformly ultimately bounded with respect to B Γ i . Moreover, the radius of B Γ i can be made arbitrary small by increasingK i such that λ min (K i ) ḡ and λ min (K i ) ḡ . Furthermore, sinceC V i ∈ F (t), its derivative follows (25) (i.e., lim sup t→∞ g(t,C V i ) − u i ≤ Γ i,2 ). Hence, all the vehicles will converge to their centroidal Voronoi configurations over F (t) and follow its trajectory associated with (25). According to Assumption 6, the movement of the multi-UAV along the trajectory of F (t) solves the sweeping coverage problem. This completes the proof.
Theorem 3. Consider a multi-UAV system of size n where each vehicle's motion model is represented by (10). Furthermore, consider using the algorithm in S1-S4 defined for a sweeping region F (t) whose dynamics is governed by (25) so that the vehicles can compute their centroidal Voronoi configurations. The bounded control law (17) solves the 3D sweeping coverage problem defined in Problem 2 under Assumptions 2-6 and the condition (26).
Proof. Again, consider the errors definitions in (28). We now define a Lyapunov candidate function as follows: , · · · , e C n (t),ė C 1 (t), · · · ,ė C n (t)) where log(v) ∈ IR 3 and cosh(v) ∈ IR 3 are defined element-wise for any vector v ∈ IR 3 , k 3 and k 4 are positive kinematic constants, and u i is defined using (17). Furthermore, let Sech 2 (v) : IR 3 → IR 3×3 be a mapping function, based on the hyperbolic secant function, which maps the vector v = [v 1 , v 2 , v 3 ] T into a diagonal matrix as follows: The time derivative of V 3 can then be obtained as follows: where Similar to the previous analysis,V 3 is negative outside the compact set 1 and Υ i,2 are given by: Thus, the system trajectories will converge to the set D Υ starting from any initial condition, and stay there forever leading to the following: Therefore, the tracking errors are uniformly ultimately bounded with respect to D Υ i . By choosing, λ min (K i ) ḡ and γ i λ min (K i ) ḡ , the errors can be made arbitrarily small. (25) and S1-S4. By Assumption 6, the multi-UAV system solves the sweeping coverage problem which completes the proof.

Validation & Discussion
Simulations were carried out to validate the performance of the proposed 3D coverage control laws in (16) and (17) using the algorithm in S1-S4 to compute the centroidal Voronoi configurations. Additionally, more simulation cases were performed to demonstrate the robustness of the proposed method and how obstacle avoidance can be incorporated within the overall framework. The following subsections provide details of these simulations and the obtained results.

Simulation Cases 1-4: Performance Validation
In the first set of simulations, a multi-UAV system of size n = 20 was used. All vehicles have been initially deployed to random locations in some predefined region . The control design parameters were chosen to be K i = diag{2.5, 0.5, 0.5} andK i = diag{0.3, 0.3, 0.3} for all vehicles.
The 3D barrier coverage problem was considered in the first two simulation cases where the unbounded control law (16) and the bounded control law (17) are used. The goal was to achieve optimal coverage over a barrier region defined according to B = {(x, y, z) : x = 20, 0 ≤ (y, z) ≤ 10}. The obtained results for this case are presented in Figures 2-5. Figures 2 and 4 show the complete trajectories taken by all UAVs and their final locations which are optimally distributed over the barrier B (rectangular area highlighted in yellow). The multi-UAV system reaches the centroidal Voronoi configurations and all vehicles form a static arrangement over B. Note that the shape of B and the number of UAVs were chosen arbitrarily just as a proof of concept. However, in practical applications, the number of UAVs and the design of the barrier can be considered as a design problem which depends on the UAVs sensing and communications capabilities. Having a larger number of vehicles may result in some overlapping between sensors field-of-view (FOV). On the other hand, It may not be possible to completely cover B using lower number of UAVs than what is needed (i.e., the combined FOV of all sensors is less than the size (area/volume) of B. Overall, the sensing performance will be maximized using the developed strategy. It can also be seen that the resultant trajectories are collision-free.
The time evolution of control inputs for all vehicles are shown in Figures 3 and 5 for both the bounded and unbounded control laws, respectively. Using (16) would require choosing the controller gain matrixK i properly in order to satisfy the constraint (11). This depends mostly on how far the vehicle is from its centroidal Voronoi configuration initially which indicates that tuningK i could not be ideal in practice. Alternatively, using (17) provides an easier way of choosing K i to ensure that the physical limit on the vehicles velocity is respected (i.e., (11) is satisfied). This can be clearly seen in Figure 5 where the vehicles speed remains constant for the first 8 s until the barrier is reached at which the velocities drop down to zero, and the vehicles become statically distributed over B. The upper bounds on ||u i || in this case was u i,max = 2.6 m/s ∀i which is in accordance with (18).
In the next two simulation cases, the sweeping coverage problem was considered. The task was to completely scan a 3D sensing region S which was defined as S = S 1 ∪ S 2 ∪ S 2 where Based on the environment shape, an initial sweeping plane F (t = 0) was determined by the multi-UAV system as F (0) = {(x, y, z) ∈ S : x = 10}. The dynamics of F (t) was also considered to be the simplest case as in (27) where the sweeping plane is moving with a constant speed of g 0 = 1.5 m/s in a progressive direction that can result in a sweeping behavior, as will be shown. In more complex cases, some other patterns could be adopted for moving F (t) as a higher level decision making which can still be done in a distributed manner as the vehicles can exchange information over the connected network. For example, the literature on coverage path planning for single-vehicle systems can be utilized in this case. Figures 6 and 7 shows the results for the sweeping coverage problem when using (16), and the results obtained when applying (17) are shown in Figures 8 and 9. For both cases, the vehicles move from their initial positions to quickly deploy over F (t) reaching their centroidal Voronoi configurations. As the position of F (t) evolves over time according to (27), the vehicles corresponding centroidal Voronoi configurations evolve accordingly sinceC V i ∈ F (t). Hence, the vehicles will start to move with the same velocity as of F (t). This can be clearly seen from Figures 7 and 9. You can see that the vehicles starts with a higher speed to reach the moving sweeping plane and achieve optimal distribution. Once this is achieved (around t = 10 s), the vehicles are no longer moving within F (t) at which all velocities converge to g 0 = 1.5 m/s in the direction of the sweeping plane's movement. It is important to notice that the speed of F (t) (i.e., g 0 ) should be slower than the maximum velocity achievable by any vehicle (u max ≥ g 0 ). After scanning the desired region S, the sweeping plane F (t) becomes static which reflects on all the vehicles as can be seen from the results where all velocities converge to 0. The complete trajectories of the vehicles along with the scanned 3D volume (highlighted in yellow) are shown in Figures 6 and 8 which confirms that the sweeping coverage problem is achieved over S. These results clearly validate the performance of the proposed 3D coverage control strategy.

Simulation Case 5: Robustness
Another simulation case was considered to show how the proposed method perform very well in situations where the number of active UAVs within the multi-UAV system change over time during a sweeping coverage mission. For example, when a number of UAVs fail, the whole group should still be able to continue their coverage task as long as there is enough number of active UAVs to finish the mission. This goes the same way when adding new vehicles to the system during the mission; however, this case was not considered here for brevity.
We considered a similar environment S and sweeping plane choice F (t) as in simulations cases 3 and 4. Different time instants of the simulation are shown in Figure 10. The Voronoi regions generated by the UAVs over F (t) with their centroids are clearly highlighted to show how they change over time as the vehicles move. Once centroidal Voronoi configurations are reached, the vehicles move according to the plane's movement as was discussed earlier. However, some vehicles fail and become inactive at certain time instants as in Figure 10e,g,i,k. Whenever this occurs, the remaining active vehicles quickly adjust their distribution over F (t) while still moving in accordance withḞ (t). For example, at t =21s, one vehicle fail as shown in Figure 10e which directly indicates a change of the Voronoi partition of F (t). This results in a change of the centroidal Voronoi configurations which causes the vehicles to quickly adapt to the situation in a robust way by moving to the new centroidal Voronoi configurations. It can be noticed that once a vehicle fail, only the vehicles in the neighborhood of that vehicle will be affected (i.e., their Voronoi regions will be extended).
Initially, the multi-UAV system had a size of n = 9. The complete collision-free trajectories of all vehicles are shown in Figure 11 where 4 vehicles have failed during the mission (inactive UAVs), and the coverage task was completed efficiently by the remaining 5 vehicles (active UAVs). This clearly shows how robust and scalable our method can be. It is also worth mentioning that changes to F (t) can be applied in real-time if its size becomes larger/smaller than what the remaining vehicles can cover based on their combined sensing FOV. Such a decision can be autonomously made by the vehicles and shared among the connected network. The next simulation case shows how the proposed control laws work when such changes to F (t) are applied which is really important when considering obstacle avoidance.

Simulation Case 6: Obstacle Avoidance
An additional simulation case was performed to show one way of incorporating obstacle avoidance capabilities within the proposed method. The simulation results for this case are given in Figures 12-14. In this case, six vehicles were initially deployed, and a rectangular region was considered to be the sweeping plane F (t). This particular choice of F (t) can be used in scenarios where the vehicles are equipped with some downward-facing sensors (ex. cameras), and their sensing FOV are simply a footprint on the ground. Even though using six vehicles may be redundant in this case as there will be overlapping in the vehicles' sensing FOV, the aim here is only to show a particular way for the multi-UAV system to avoid obstacles. The adopted approach is simply to manipulate the shape, size, orientation and/or velocity of the sweeping plane F (t) to safely avoid detected obstacles. Note that in this case, only vehicles within a sensing range from obstacles can determine such required action which then sent to the remaining vehicles over the connected network. Obstacle avoidance is ensured in this case because vehicles are guaranteed to maintain their motion within F (t) once its reached using the proposed control laws (i.e., p i (t) ∈ F (t), ∀t > t * where t * is the time at which the vehicles reach F (t)).
It can be seen from Figure 12a,b that the vehicles are approaching F (t) since they were initially deployed away from it. Simultaneously, the vehicles move along F (t) to reach centroidal Voronoi configurations. At t = 15 s, two obstacles are detected along the way of the sweeping plane's movement by two nearby vehicles. The vehicles decide to dynamically change the size of F (t) to avoid both obstacles as shown in Figure 12g-j. A different decision is also made by the multi-UAV system when detecting another obstacle at t = 49.5 s (Figure 12k). At this time, the sweeping plane is tilted in the z-direction as can be seen from Figure 12l-n.
Overall, this approach has good potential for motion coordination of multi-vehicle systems especially when considering obstacle avoidance compared to artificial potential based approaches where sometimes there will be conflicts between achieving obstacle avoidance and avoiding collision with other vehicles.

Approach
The proposed strategy constraints the vehicles' movement within a desired plane. This can be sufficient for sweep coverage problems as the region's movement will ensure that the 3D environment is covered. However, when considering flocking problems and swarms with large number of vehicles, it becomes limiting to constrain the movement of all vehicles within a planar region. Thus, this section generalizes the proposed method by introducing the use of multi-regions.
Let F(t) = {F 1 (t), ·, F m (t)} be a list of m 3D polygonal regions defined by their vertices such that the following condition is satisfied: Each region has dynamics governed by: where r i and g i (t, r i ) are as defined in (25), and g i (t, r i ) is subject to the same constraint in (26). Furthermore, let Q max ∈ IR m be a vector representing the maximum capacity for a region (i.e., the maximum number of vehicles it can include without violating the safety condition). Additionally, define Q(T) ∈ IR m as a vector of the instantaneous number of vehicles allocated for each region. The proposed approach using Centroidal Voronoi tessellations along with the control law (16) (or (17)) can guarantee collision avoidance among vehicles moving towards the same region. However, the control law (16) (or (17)) need to be modified to ensure that vehicles can avoid collisions when moving towards centroidal Voronoi tessellations of different regions. Thus, the following control law is proposed (based on (17)): where K ij is a positive semi-definite diagonal matrix, and d c > 0 is some threshold distance at which the repelling force (the second component of (37)) is applied. Once all vehicles converge to their associated regions, this component vanishes given that the regions design is feasible. Note that the control law (37) allows us to tune the gains K i and K ij while ensuring that the condition (11) is satisfied.
To maintain the distributed nature of the overall strategy, it is assumed that F(t) is shared among the connected networked vehicles along with the regions' dynamical model, the maximum capacity vector Q max and the vector Q(t). Thus, whenever a vehicle decides to pick a region F i (t), it can chose the closest one where Q i (t) < Q max .

Simulations
Several simulations have been carried to verify the generalized approach presented in this section. The simulation cases show how the swarm converges to the desired structure and how the flocking objective can be achieved. Different swarm sizes, desired structures and movement patterns have been considered.
The obtained results for four cases are shown in Figures 15-18. In all cases, the vehicles are initially deployed to random locations. Figures 15-18a show the trajectories taken by all vehicles, which verify that the swarm can converge to and hold the desired structure while following the desired movement pattern. Figures 15-18b show the multi-regions considered and the final locations of the vehicles. Notice that there is no constraints on the shapes of the regions except that they should not intersect according to (35).

Implementation Using a Multi-Quadrotor System
The proposed control methods are developed based on the general kinematic model in (10) which is applicable to different UAV types and AUVs. An example way in implementing the sweep coverage control method is presented in this section for a particular application in precision agriculture using a group of quadrotor UAVs.

Quadrotor Dynamics
Now, we will extend the model in (10) to include the dynamics of quadrotor-type UAVs according to the following model: where the velocity vector u i is replaced by v i for convenience, g is the gravitational acceleration, e 3 = [1 1 1] T , m i is the UAV mass, and I i is the inertia matrix. Recall that the vehicle's position and velocities are expressed in the inertial coordinate frame {I}. Consider also another reference frame attached to the UAV body with its origin at the center of mass. This frame is referred to as the body-fixed frame {B i }. The orientation of the UAV is represented using a rotation matrix R i between {B i } and {I}. The rate of change in the vehicle's orientation is denoted by Ω i (i.e., angular velocity) which is expressed in the {B i } frame. The control inputs for this model are the collective thrust T i ∈ IR and the body-torques vector τ i ∈ IR 3 . For more details about quadrators modelling, the reader is referred to [84,85].

Tracking Control
There are different ways to implement the proposed sweeping control strategy. One possible approach is to directly couple the velocity commands in (16) (or (17)) into attitude and thrust control inputs. Another possible direction is to use the determined centroidal Voronoi configurations by the algorithm S1-S4 as goal positions with a position tracking controller. The latter approach is considered here for a control design based on the differential-flatness property of the quadrotor dynamics and the sliding mode control technique.
Let us redefine the position and velocity tracking errors as follows: consider sliding surfaces σ i such that: where µ 1 > 0, and L 1,i is a positive definite diagonal matrix. Note that this choice will ensure that the velocities can remain bounded by λ max (L 1,i ) to account for physical limits. A desired acceleration command can then be computed using: 1 m i a des,i =ge 3 + L 2,i tanh(µ 2 σ i ) where L 2,i and L ij are positive definite diagonal matrices. The additional term in (45) provides more safety as it represents a repelling force from vehicles closer than some critical distance d c > 0 in the neighbourhood of vehicle i. The differential-flatness property is then utilized to achieve the acceleration command in (45) while maintaining some desired yaw angle ψ re f ,i using the following equations: A low-level control is then used to compute τ i to ensure that the vehicle can achieve the desired attitude R des,i . Further details about the differential-flatness property of quadrotor dynamics can be found in [85,86].

Software-in-the-Loop Simulations
The performance of the suggested sweeping coverage control for multi-quadrotor systems was evaluated using Software-in-the-Loop (SITL) simulations using Gazebo and the Robot Operating System (ROS) framework. This allows us to test the computational performance of the production code which can be used directly on the real vehicle's onboard computer to implement our algorithms and control strategies in real-time. A simple scenario was considered where a group of three quadrotors were needed to survey a crops fieldS = {x, y, z ∈ IR : 0 ≤ x ≤ 32, 0 ≤ y ≤ 40, z = 0}. Each vehicle is assumed to have an onboard camera providing a sensing footprint A i ⊂ S dependent on its characteristics and the UAV's altitude. For simplicity, an obstacle-free environment was used as shown in Figure 19.
We used the open-source PX4 flight stack as an implementation for the low-level controller and an extended Kalman filter for states estimation considering noisy measurements. Our tracking control logic was implemented using C++ providing thrust and attitude inputs at a rate of 100Hz. Furthermore, centroidal Voronoi configurations computations were implemented using Python utilizing some of the available tools from the "scipy.spatial" Python module to determine Voronoi cells as a convex polygon in accordance with Assumption 1. Furthermore, the algorithm in S1-S4 was used to obtain the centroids. This turned out to be very computationally efficient since the computations rely on closed-form expressions.
In order to achieve the coverage task in hand, an initial rectangular sweeping plane was selected as: F (0) = {x, y, z ∈ IR : 0 ≤ x ≤ 8, 0 ≤ y ≤ 2, z = 2} where it was required for the UAVs to fly at a fixed altitude of 2m. The classical lawn-mower pattern was considered to generate the trajectory of F (t) to completely scanS while keeping the yaw angle at ψ re f = 0 without loss of generality. The following parameters were used in the simulation: L 1,i = diag{2, 2.5, 3}, L 2,i = diag{2.5, 3, 4}, L ij = diag{1, 1, 1}, d c =0.5 m, µ 1 = 2 and µ 2 = 1.5 (where diag{·} represents a diagonal matrix). During the simulation, all measurements and data were captured as a bag file using available ROS recording tools. The results are plotted using MATLAB which are shown in Figures 20-25. In Figure 20, coordinates of all vehicles are plotted with respect to time. The overall paths taken by the vehicles are shown in Figures 21 where the highlighted rectangular area represents the crops field to be surveyed (i.e.,S). Note that in this time, the vehicles were moving along some other region S whose projection on the ground floor isS. Thus, the trajectory of F (t) was designed to coverS given that the equipped sensors have a collective footprint which ensures the optimal coverage ofS. It is clear from these results that the overall motion is collision-free, and the lawn-mower pattern decided for F (t) was followed by the whole group. Velocity norms are shown in Figure 23. Overall the vehicles are moving with a constant speed g 0 = 1 m/s except for some changes in velocity at times where the vehicles are turning (t ≈ 60 s, 110 s,160 s) in order to avoid collisions with each others. Overall, velocities remain within the desired limits of v ≤ 3.5 m/s except at the start when vehicles were taking off which is based on a different control law. The high-level control commands, namely roll, pitch, yaw and thrust, are shown in Figures 24 and 25 where the thrust is normalized in the range [0, 1]. This simulation case shows that our control can be successfully implemented using multi-quadrotor systems.

Conclusions & Future Work
Distributed cooperative control methods were proposed in this paper for multi-UAV systems to address coverage problems in 3D spaces including barrier and sweeping problems. Additionally, a generalization of these approaches was given to address flocking problems. The development of control laws was based on general kinematic model which make them applicable to different UAV types and autonomous underwater vehicles as well. A region-based control mechanism was adopted to constrain the movement of the vehicles within some desired region providing a computationally efficient solution that is both robust and scalable. Mathematical analysis was performed to ensure the stability of the developed control laws. Several simulations have been performed confirming the performance of the developed methods for 3D static and dynamic coverage problems. Additional simulation cases were also presented to show how well the suggested approach can handle obstacle avoidance as well as being robust against vehicles' failure. To evaluate the computational performance, software-in-the-loop simulations, using Gazebo and ROS, were also carried out considering a special case of using a multi-quadrotor system in precision agriculture. To that end, specific implementation details were also proposed based on quadrotor dynamics. Future work can consider practical implementation and developing proper ways where the vehicles can collaboratively decide the shape of the dynamic region based on the multi-UAV system characteristics and the targeted sensing region size.