Distributed 3D Navigation of Swarms of Non-Holonomic UAVs for Coverage of Unsteady Environmental Boundaries

Abstract

A team of non-holonomic constant-speed under-actuated unmanned aerial vehicles (UAVs) with lower-limited turning radii travel in 3D. The space hosts an unknown and unpredictably varying scalar environmental field. A space direction is given; this direction and the coordinate along it are conditionally termed as the “vertical” and “altitude”, respectively. All UAVs should arrive at the moving and deforming isosurface where the field assumes a given value. They also should evenly distribute themselves over a pre-specified range of the “altitudes” and repeatedly encircle the entirety of the isosurface while remaining on it, each at its own altitude. Every UAV measures only the field intensity at the current location and both the Euclidean and altitudinal distances to the objects (including the top and bottom of the altitudinal range) within a finite range of visibility and has access to its own speed and the vertical direction. The UAVs carry no communication facilities, are anonymous to one another, and cannot play distinct roles in the team. A distributed control law is presented that solves this mission under minimal and partly inevitable assumptions. This law is justified by a mathematically rigorous global convergence result; computer simulation tests confirm its performance.

1. Introduction

The need to explore various environmental boundaries has motivated extensive research on using mobile robotic platforms for such a purpose; see, e.g., Refs. [1,2,3,4,5,6,7,8,9,10,11,12,13] and the literature therein. A typical mission is to find and arrive at the level set where an unknown environmental field assumes a specific value and then sweep the entirety of this set, thus exhibiting and putting under control the border of the region with the greater field values. Examples include finding the flows of air pollutants or contaminant clouds [14] and tracking zones of turbulence or high radioactivity level, to name just a few. In such missions, typical challenges include a paucity of a priori information about the field, obsolescence of the data collected online due to the field changes, and the capacity of the available sensors to measure only the field value at the current location via immediate contact with the sensed entity, e.g., with a transparent gas.

Recently, much attention has been given to navigation algorithms that enable mobile robots to localize, approach, and cover the environmental boundary of interest. A large group of the algorithms relies on access to the field’s gradient [8,12,15,16,17,18]. This group is exemplified by, e.g., the methods based on multi-agent estimation of the gradient [17], cooperative contour estimation [2,16], and gradient-based artificial potentials [18]. However, the possibility to directly measure the field’s gradient is uncommon, whereas reliable gradient estimation from noisy measurements of the field value is still an intricate challenge in a practical setting [19,20]. Also, such estimation calls for measurements in neighboring locations distributed across all dimensions, whereas exploration of an environmental boundary motivates to place the sensors on this lower-dimensional structure. Finally, communication constraints may hinder transfers of field measurements to the gradient estimator, wherever it may be built.

The alternative, gradient-free methods do not attempt to assess the gradient and are well fitted to the situation of pointwise measurement of the field value only. Some such methods (exemplified by [4,21]) implement oscillations around a profitable path, thus enabling the robot to collect field data from a whole corridor hosting this path. This approach raises concerns about a waste of resources due to systematic and mutually nullifying shifts sideways. Common image segmentation techniques are used in [1] to monitor a forest fire-front by a team of UAVs. However, these findings are not rigorously and completely justified. A PID controller empowered by an extended Kalman filter and adaptive crossing angle correction scheme is justified in [22] for a holonomic planar mobile robot. To drive a Dubins car-like robot along an isoline of a planar field, a PD controller is presented in [23], and its local convergence in radial harmonic fields are proven, whereas [24,25] offer sliding mode controllers whose global convergence is rigorously justified for generic smooth steady [24] and time-varying [25] fields. The findings of [24] are extended to the case of multiple robots in [26], where the algorithm also ensures effective self-deployment of the robots over the environmental boundary.

The expansion of drone technology motivates the interest to navigate drones using all three dimensions. However, the literature on the sensor-based robotic tracking environmental boundaries has focused so far on the case of 2D or those reducible to 2D. Few exceptions [12,27] deal with a single tracking robot; it is supported by a sensing robot in the context of [12]. The controller from [12] assumes communicating robots modeled as simple integrators and also that the field evolves subject to an advection-diffusion equation with the fully known constant parameters; these assumptions about the field are generally challenged in practice. The control algorithm from [12] requires a computationally expensive online solution of a partial differential equation, and the completeness of the coverage of the level set is not addressed. In [27], a gradient-free control law is presented that drives a non-holonomic underactuated mobile robot to an unknown and unsteady environmental boundary in 3D and then ensures its exhaustively sweeping.

Contrary to the scenario of a single robot, the strength of drone technology greatly stems from the use of large teams of simple and low-cost devices. Reaping this benefit requires multi-agent control strategies that are robust, fault-tolerant, distributed, and homogeneous in the sense of identical roles of the teammates. Other requirements include low consumption of energy, computational, and communication resources, as well as rigorous guarantees for global convergence. Among the survey of papers on sensor-based robotic tracking environmental boundaries in 3D, the authors, however, failed to come across one addressing these issues.

This paper seeks to fill this gap while combining the above issues with that of constraints due to non-holonomy, under-actuation, and a limited control range of the robots. Whether the plenty of the identified factors and requirements allows solving the mission by a computationally inexpensive, low-level controller that directly converts the current observation into the current control and is, nevertheless, justified by a rigorous global convergence result? The paper answers in the affirmative and offers respective details, including the techniques and concepts of justification.

Specifically, we consider a swarm of UAVs whose kinematics are described by the Dubins vehicle model [28,29]. Every UAV moves in 3D with a constant speed in the longitudinal direction and is steered by the yawing and pitching rates, which are limited in magnitude. This model applies to, e.g., fixed-wing UAVs, torpedo-like underwater drones, surface vessels, and various rotorcraft [28,29].

The UAVs cannot distinguish among their peers and cannot play distinct roles in the team; they are unaware of the team’s size. Any UAV has access to the vertical direction, is aware of its own speed, can assess the altitudinal and Euclidean distances to the objects within a finite visibility range, and measures the value of an unsteady and unpredictable scalar field at the current location. All UAVs should find and reach the moving and deforming level set (isosurface) where the field assumes a given value. They should also distribute themselves into the densest net across a pre-specified altitudinal range (this may be, e.g., the range of particularly important altitudes or those at which the UAVs can operate). After this, every UAV should repeatedly circumnavigate the isosurface at its own altitude selected in a distributed fashion, thus forming an altitudinally densest and horizontally complete dynamic barrier around the isosurface for exposition, surveillance, processing, or protection. All these goals should be achieved via independent decisions of the UAVs according to a common rule and with no communication among them.

The paper presents a gradient-free navigation law that meets the above requirements. Moreover, we first disclose conditions necessary for the mission to be achievable. Then we show that the proposed law solves the mission under only a slight and somewhat inevitable enhancement of these conditions. This is done by means of a mathematically rigorous global convergence result. Basic theoretical findings are confirmed and complemented by computer simulation tests.

This paper develops some ideas reported in [25,27]. However, Ref. [25] deals with a planar workspace so that its findings are insufficient to cope with special challenges inherent in 3D environments. Meanwhile, [27] considers a single-robot scenario and so contributes nothing to the issue of inter-UAV cooperation, which is the major focus of the current paper. Moreover, despite some similarity in processing the field measurements, analysis of the entailed behavior with respect to isosurface finding and tracking aspects must be fully updated as compared with [27] due to the critical coupling of the concerned control loop with that of regulating the inter-UAV altitudinal gaps.

The body of the paper is organized as follows. Section 2 and Section 3 introduce the problem setup and the control law, respectively. Section 4 is devoted to necessary conditions for the mission feasibility and the assumptions of our theoretical analysis. The main results are stated in Section 6. Section 7 reports on computer simulations, and Section 8 offers brief conclusions. All proofs are placed in Appendix A, Appendix B, Appendix C and Appendix D.

The following general notations are used in the paper: $:=$ means “is defined to be”, “is used to denote”; the dot · in notations like $f(·)$ is the placeholder of the argument of the function f; the symbols $〈·;·〉$, $\Vert ·\Vert$, and × stand for the standard inner product, Euclidean norm, and cross product in ${\mathbb{R}}^3$, respectively.

2. Problem of Sweep Coverage of an Isosurface

A team of N UAVs travels in 3D. Every UAV moves with a constant speed in the longitudinal direction over paths of bounded curvatures, driven by pitching and yawing rates limited in absolute value. There is an unknown and unsteady environmental field described by a scalar function $F(t,\mathit{r})\in \mathbb{R}$ of time t and space location $\mathit{r}\in {\mathbb{R}}^3$. All UAVs should arrive at the locus of points (called isosurface) $S_t\left(f_{\star }\right):=\{\mathit{r}:F(t,\mathit{r})=f_{\star }\}$ with the pre-specified field value $f_{\star }$. Then they should sweep this isosurface while uniformly distributing themselves over it. The UAVs are not equipped with communication facilities. So coordination of their motions should be based on only individual sensory data and be achieved in a fully distributed fashion.

Further specification of the targeted distribution assumes that a certain space direction is given by a unit vector $\mathit{h}\in {\mathbb{R}}^3$; the associated coordinate $h\left(\mathit{r}\right)$ of point $\mathit{r}\in {\mathbb{R}}^3$ is loosely referred to as altitude. The mission is confined to a certain altitudinal range ${\mathcal{H}}_{al}=[{ɦ}_{-},{ɦ}_{+}],{ɦ}_{-}<{ɦ}_{+}$; for example, this may be the range of particularly important altitudes or altitudes at which the UAVs are able to operate. Self-distribution of the UAVs into the densest net across the range ${\mathcal{H}}_{al}$ should be achieved, whereas each of them should fully circumnavigate the isosurface at its own altitude.

The ith UAV has access only to the field value $f_i\left(t\right):=F[t,{\mathit{r}}_i\left(t\right)]$ at its own current location ${\mathit{r}}_i={\mathit{r}}_i\left(t\right)$, and has no idea about the distance to or bearing of the targeted isosurface $S_t\left(f_{\star }\right)$. The ith UAV also has access to $\mathit{h}$ in its frame of reference, is aware of its own speed $v_i$, and can assess both the altitudinal and Euclidean distances to the objects, including the top/bottom ${ɦ}_{\pm }$ of the altitudinal range ${\mathcal{H}}_{al}$, that lie within a given “visibility” range $d^{\mathrm{vis}}>0$. The UAVs do not know their total number and cannot distinguish between their peers or play different roles in the team. The last circumstance implies that all UAVs should be driven by a common control rule and cannot be assigned individual serial numbers that influence the control input.

To unveil the structure of the densest net across ${\mathcal{H}}_{al}$, we denote by

$$\mathbf{mesh}\left[\mathcal{H}\right]:=\underset{h\in {\mathcal{H}}_{al}}{max}\underset{h^{\prime }\in \mathcal{H}}{min}|h-h^{\prime }|$$

the worst-case distance from a point of ${\mathcal{H}}_{al}$ to a finite set $\mathcal{H}\subset {\mathcal{H}}_{al}$. It is easy to see that the minimum of $\mathbf{mesh}\left[\mathcal{H}\right]$ over all sets $\mathcal{H}\subset {\mathcal{H}}_{al}$ with N elements is attained at the set

$$\mathcal{H}:=\{{ɦ}_i:={ɦ}_{-}+\Delta ɦ+2i\Delta ɦ,i\in [0:N-1]\},\mathrm{where}\Delta ɦ:=[{ɦ}_{+}-{ɦ}_{-}]/\left(2N\right).$$

(1)

Now we flesh out the targeted collective behavior of the considered team of UAVs.

Definition 1.

The team is said to form the densest horizontally sweeping net on the isosurface $S_t\left(f_{\star }\right)$ in the range from ${ɦ}_{-}$ to ${ɦ}_{+}$ if the UAVs can be labeled with i from 0 to $N-1$ so that $h_i\left(t\right)\to {ɦ}_i,f_i\left(t\right)\to f_{\star }$ as $t\to \infty$, where $h_i\left(t\right):=h\left[{\mathit{r}}_i\left(t\right)\right]$ is the altitude of the ith UAV.

The imposed information constraints mean that any UAV is unaware of the altitude ${ɦ}_i$ assigned to this UAV in Definition 1, may have no access to the end-points of the altitudinal range since they are beyond its range of visibility, and cannot receive these data from more informed teammates (if they exist) due to the lack of communication facilities.

It is required to design a common control rule by executing which every UAV individually builds its own control based on the available data, whereas the entire team acquires the property described in Definition 1. Given an ever-growing use of mass-produced, cheap, relatively small-sized, and, as a result, energy and computationally constrained drones, this rule is welcome to be computationally inexpensive and exhibit a regular energy-efficient behavior. Whether the entire range of the above rather diverse and partly contradictory wishes and goals may be compromised and attained?

For theoretical analysis, we employ a truncated model of the kinematics of a rigid-body robot moving with a constant speed in the longitudinal direction [30,31,32,33]. This model disregards the roll motion and is used in vector form borrowed from [31,32]:

$${\dot{\mathit{r}}}_i=v_i{\mathit{e}}_i,{\dot{\mathit{e}}}_i={\mathit{u}}_i\in {\mathbb{R}}^3,〈{\mathit{u}}_i;{\mathit{e}}_i〉=0,\Vert {\mathit{u}}_i\Vert \le {\widehat{u}}_i.$$

(2)

Here ${\mathit{u}}_i$ is the control, the upper bound ${\widehat{u}}_i>0$ on its magnitude is given, $v_i>0$ is the constant speed of the UAV, ${\mathit{e}}_i$ is the unit vector along the centerline of the ith UAV (see Figure 1), and the third equation in (2) “keeps” the length of ${\mathit{e}}_i$ constant. The model (2) captures the robot’s capacity to travel over space paths whose curvature radii $\ge v_i/{\widehat{u}}_i$. The scope of applicability of this model is discussed at length in Remark 2.1 in [31] and includes fixed-wing UAVs, various rotorcraft, and torpedo-like underwater drones. Due to a one-to-one correspondence ${\mathit{u}}_i\leftrightarrow (q_i,r_i)$ in Remark 2.1 in [31], the vector control input ${\mathit{u}}_i$ can be replaced by the pitching $q_i$ and yawing $r_i$ rates. In fact, (2) is a 3D extension of the standard Dubins-vehicle model of an aircraft or boat in a plane [34,35,36,37].

3. Proposed Hybrid Controller

It uses the following tunable free parameters

$$d_{-}^{\mathrm{vis}}\in (0,d^{\mathrm{vis}}),T_{\mathrm{tr}}>0,{\Delta }_h>0,{\delta }_f>0,\gamma >0,{\widehat{u}}^h>0,{\widehat{u}}_i^f>0$$

(3)

and two functions $\chi$ and ð, where $\chi$ maps $\mathbb{R}$ to $\mathbb{R}$ and ð maps $[0,\infty )$ to $[0,1]$. The choice of these parameters and functions is discussed in Section 6.

For any i, the ith UAV builds and updates the set $E_i=E_i\left(t\right)$ of its essential neighbors by executing the following instructions:

$$\begin{array}{c}{E}_i(0-):=\varnothing ,\hfill \\ \mathrm{whenever}\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert \le d_{-}^{\mathrm{vis}}\mathrm{for}\mathrm{some}j\notin E_i,\mathrm{add}\mathrm{this}j\mathrm{to}\mathrm{the}\mathrm{set}{E}_i,\\ \hfill \mathrm{whenever}\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert \ge d^{\mathrm{vis}}\mathrm{for}\mathrm{some}j\in E_i,\mathrm{remove}\mathrm{this}j\mathrm{from}{E}_i.\end{array}$$

(4)

The set $E_i(0-)$ can be immediately replenished due to the second line in (4). The use of not equal but different parameters $d_{-}^{\mathrm{vis}}<d^{\mathrm{vis}}$ in the second and third lines, respectively, is aimed at suppressing excessive sliding-mode phenomena via preventing the situations where a just enrolled peer j should be immediately excluded and vice-versa. The sets $E_i^{+}$ and $E_i^{-}$ of essential higher and lower neighbors, respectively, are defined as

$$E_i^{+}:=\left\{j\in E_i:h_j-h_i>0\right\},E_i^{-}:=\left\{j\in E_i:h_j-h_i<0\right\}.$$

(5)

The ith UAV can compute these sets from the available sensory data since the definitions of these sets use only the relative altitudes and the distances to the teammates within the visibility range, which data are accessible to the ith UAV. The scaled higher ${\mathring{h}}_i^{+}$ and lower ${\mathring{h}}_i^{-}$ altitudinal gaps near the ith UAV are defined to be

$${\mathring{h}}_i^{\pm }:=\left\{\begin{array}{cc}{\displaystyle \underset{j\in E_i^{\pm }}{min}\pm (h_j-h_i)}\hfill & \mathrm{if}{E}_i^{\pm }\ne \varnothing ,\mathrm{otherwise}\hfill \\ 2|{ɦ}_{\pm }-h_i|\hfill & \mathrm{if}{d}^{\mathrm{vis}}>|h_i-h_{\pm }|,\hfill \\ {\Delta }_h\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(6)

Here either + should be put in place of ± everywhere, or the same should be performed with − instead of +. If $E_i^{\pm }=\varnothing$ and $d^{\mathrm{vis}}>|h_i-{ɦ}_{\pm }|$, there are no essential neighbors between the UAV at hands and the top/bottom of the altitudinal range ${\mathcal{H}}_{al}$; then ${\mathring{h}}_i^{\pm }$ is twice the altitudinal distance to the top/bottom. Given a UAV, the gaps ${\mathring{h}}_i^{\pm }$ do not depend on the enumeration of the UAVs and so are computable in the current situation where the teammates are anonymous to one another; see Remark 1 for more details.

By the third equation in (2), the control input ${\mathit{u}}_i$ must lie in the plane ${\mathit{e}}_i^{\perp }$ normal to the centerline vector ${\mathit{e}}_i$ (the pitch-yaw plane). We define ${\mathit{u}}_i$ in a special orthonormal basis ${\mathit{h}}_i^{\mathrm{py}},{\mathit{h}}_i^{\mathrm{py}}\times {\mathit{e}}_i$ of ${\mathit{e}}_i^{\perp }$, where ${\mathit{h}}_i^{\mathrm{py}}$ is the orthogonal projection of the vertical vector $\mathit{h}$ onto the pitch-yaw plane ${\mathit{e}}_i^{\perp }$, which projection is then normalized to the unit length:

$${\mathit{h}}_i^{\mathrm{py}}:=(\mathit{h}-{\mathit{e}}_{i}cos{\theta }_i)/sin{\theta }_i.$$

(7)

Here ${\theta }_i$ is the angle between the centerline ${\mathit{e}}_i$ and the vertical $\mathit{h}$. As will be shown in (i) of Theorem 1, our controller ensures that any UAV is always non-vertical, so both the vector ${\mathit{h}}_i^{\mathrm{py}}$ and the considered basis are well-defined.

Our control law is hybrid with two modes: $\mathfrak{A}$ (approaching the isosurface) and $\mathfrak{M}$ (main mode). Any UAV i starts in $\mathfrak{A}$ and then switches to $\mathfrak{M}$. The role of mode $\mathfrak{A}$ is to find the targeted isosurface and to arrive at its close vicinity defined as the locus of points where the field value differs from $f_{\star }$ by not more than a pre-specified and nominally small ${\delta }_f$ from (3). Vertical distribution of the team is the task of the next mode $\mathfrak{M}$ so that the global search of the isosurface during mode $\mathfrak{A}$ is not disturbed by control signals aimed for other purposes. The switching rule is as follows:

$$\mathfrak{A}\mapsto \mathfrak{M}\mathrm{whenever}|f_i-f_{\star }|\le {\delta }_f\mathrm{and}t\ge T_{\mathrm{tr}}.$$

(8)

The control rule invokes (like (8) does) parameters from (3) and is as follows:

$${\mathit{u}}_i={\widehat{u}}_i^f·\mathbf{sgn}\left[{\dot{f}}_i+\chi (f_i-f_{\star })\right]{\mathit{h}}_i^{\mathrm{py}}\times {\mathit{e}}_i-\frac{{\widehat{u}}^h}{v_{i}sin{\theta }_i}·\mathbf{sgn}\left[{\dot{h}}_i-v_i\right]{\mathit{h}}_i^{\mathrm{py}},$$

(9)

where

$$v_i:=\left\{\begin{array}{cc}0\hfill & \mathrm{in}\mathrm{mode}\mathfrak{A},\hfill \\ ð(t-T_i^{\mathrm{sw}})·\left[\mathsf{\Xi }\left({\mathring{h}}_i^{+}\right)-\mathsf{\Xi }\left({\mathring{h}}_i^{-}\right)\right]\hfill & \mathrm{in}\mathrm{mode}\mathfrak{M}.\hfill \end{array}\right.$$

(10)

Here $\chi$ and ð are functions that are to be chosen by the designer of the controller (see (20)–(29) for details), $T_i^{\mathrm{sw}}$ is the time of the transition $\mathfrak{A}\mapsto \mathfrak{M}$ in (8), and

$$\mathsf{\Xi }\left(L\right)=\gamma min\{{\Delta }_h;L\}$$

(11)

is a linear function of $L\ge 0$ with saturation at the threshold ${\Delta }_h>0$ mentioned in (3).

The ith UAV can compute the derivative ${\dot{h}}_i=v_i〈\mathit{h};{\mathit{e}}_i〉$ since it is aware of its own speed $v_i$ and has access to the vertical vector $\mathit{h}$ in its local frame illustrated in Figure 1. Numerical differentiation can be used to assess the time-derivative ${\dot{f}}_i$ of the sensor readings $f_i\left(t\right)$. Estimating the derivatives from noisy data is a well-researched discipline offering many methods; any of them is acceptable and welcome to implement the controller. Among these methods are, e.g., optimal schemes based on stochastic models, observers with sliding modes, difference methods; see [19,38,39] for a survey.

The proposed design of the control system is illustrated in Figure 2. The block conditionally entitled “inclinometer” is responsible for access to the vertical vector $\mathit{h}$ in the local frame of reference of the UAV at hand. It thus provides access to the angle ${\theta }_i$ between $\mathit{h}$ and the centerline ${\mathit{e}}_i$ of the UAV and the vector (7). The block in the lower right corner of the diagram illustrates the procedure (4) of forming the set of essential neighbors. The lower and upper positions of the switch in the diagram correspond to modes $\mathfrak{M}$ and $\mathfrak{A}$, respectively.

Except for the coefficient ${\widehat{u}}_i^f$, the control rule (9) is common for all robots. Remark 2 discusses when complete uniformity can be achieved by picking ${\widehat{u}}_i^f$ common for all i.

Remark 1.

To facilitate understanding the procedure for determining ${\mathring{h}}_i^{\pm }$, its first step was pictured as building the sets $E_i,E_i^{\pm }$ of labels j. However, the robots cannot figure out these labels. So actually, these “absolute” labels are not involved: on its own choice, robot i labels the available relative distances $h_j-h_i$, uses time-invariant labels for continuously changing distances, and processes exactly these labels. All of that is performable based on the available data.

Since the control law (9) and (10) is discontinuous, the solution of the closed-loop system is meant as that of the differential inclusion obtained via Filippov’s convexification method [40]. Given an initial state, a solution exists and does not blow up in a finite time due to the boundedness of the controls.

4. Mission Feasibility and Assumptions

To avoid overly restrictive assumptions in our theoretical analysis, we first disclose conditions necessary for the mission feasibility. They display the necessary balance between the level of maneuverability of the UAVs and the challenges from the contortions and motion of the targeted isosurface. Our assumptions will be only slight and partly inevitable enhancements of these necessary conditions. To disclose them, we need the following characteristics of the unsteady environmental field:

• $\nabla F$, spatial gradient of the field;
• $N(t,\mathit{r})=\frac{\nabla F(t,\mathit{r})}{\Vert \nabla F(t,\mathit{r})\Vert }$, unit vector normal to the associated isosurface (AI) that passes through location $\mathit{r}$ at time t;
• ${\alpha }_h=arcsin〈N;\mathit{h}〉$, angle from this vector N to the horizontal planes;
• $N_{\mathrm{hor}}=(N-\mathit{h}sin{\alpha }_h)/cos{\alpha }_h$, projection of this vector N onto the horizontal plane normalized to the unit length;
• ${\mathit{h}}_{\tan }=(\mathit{h}-Nsin{\alpha }_h)/cos{\alpha }_h$, normalized projection of the unit vertical vector $\mathit{h}$ onto the plane tangent to the associated isosurface;
• $S_t^{\mathrm{hor}}\left(f_{\star }\right|h):=\{\mathit{r}\in S_t\left(f_{\star }\right):h\left(\mathit{r}\right)=ɦ\}$, horizontal section (of the $f_{\star }$-isosurface) at the altitude h;
• $\overrightarrow{\tau }=\mathit{h}\times N/cos{\alpha }_h$, unit vector tangential to the horizontal section that passes through location $\mathit{r}$ at time t;
• $\mathit{II}\left[V\right]$, second fundamental (quadratic) form of AI, i.e., the signed curvature of the intersection of AI with the span of a unit tangent vector V and N; see Section 4 in [41];
• $\varkappa$, maximal (in absolute value) eigenvalue of this quadratic form;
• $\lambda (t,\mathit{r})$, front velocity of the associated isosurface;
• $\alpha (t,\mathit{r})$, front acceleration of the associated isosurface;
• $\overrightarrow{\omega }(t,\mathit{r})$, angular velocity of rotation of the associated isosurface;
• $\rho (t,\mathit{r})$, density of isosurfaces, which evaluates their number (assessed by the range of the field values) within the unit distance from the associated isosurface;
• $g_{\rho }(t,\mathit{r})$, proportional growth rate of the density $\rho$ with time;
• $\nabla \nabla \rho (t,\mathit{r})$, proportional tangential gradient of the density $\rho$,
• $n_{\rho }(t,\mathit{r})$, proportional growth rate of the density $\rho$ under the normal shift.

The last seven quantities are rigorously defined in Appendix A in [27].

By Definition 1, ideally carrying out the mission includes moving over the horizontal section $S_t^{\mathrm{hor}}\left(f_{\star }\right|ɦ)=\{\mathit{r}\in S_t\left(f_{\star }\right):h\left(\mathit{r}\right)=ɦ\}$ at a fixed altitude $ɦ={ɦ}_i\in {\mathcal{H}}_{al}$. Since the size N of the team and so the altitudes ${ɦ}_i$’s are unknown to the team members and no UAV is assigned its own altitude ${ɦ}_i$ a priori, it is fair to require that any UAV be able to trace the whole of $S_t^{\mathrm{hor}}\left(f_{\star }\right|ɦ)$ within the operational zone$\mathcal{OZ}$ at any altitude $ɦ\in {\mathcal{H}}_{al}$ in the given altitudinal range ${\mathcal{H}}_{al}$. If this requirement is met, the isosurface is said to be trackable by this UAV. For the sake of convenience, we define the zone $\mathcal{OZ}$ in terms of the extreme values $f_{-}<f_{+}$ ($f_{\star }\in (f_{-},f_{+})$) achievable by the field F in this zone:

$$\mathcal{OZ}:=\{(t,\mathit{r}):f_{-}\le F(t,\mathit{r})\le f_{+},h\left(\mathit{r}\right)\in {\mathcal{H}}_{al}\}.$$

(12)

Since any UAV (2) can trace only regular (i.e., differentiable with a nonzero derivative) curves, the above trackability may hold only if any curve $S_t^{\mathrm{hor}}\left(f_{\star }\right|ɦ),ɦ\in {\mathcal{H}}_{al}$ is regular. This may be violated not only because of the non-smoothness of the field but also due to the zero spatial gradient. Hence the following is compelled by necessity.

Assumption 1.

In an open vicinity of the operational zone (12), the field F is twice continuously differentiable, is not singular $\nabla F\ne 0$, and the horizontal section $S_t^{\mathrm{hor}}\left(f_{\star }\right|ɦ)$ of the $f_{\star }$-isosurface at any altitude $ɦ\in {\mathcal{H}}_{al}$ is not empty.

Conditions necessary for the mission feasibility are as follows.

Lemma 1.

Let the isosurface $S_t\left(f_{\star }\right)$ be trackable by the ith UAV. Then at any temporal-spatial point of $S_t\left(f_{\star }\right)\cap \mathcal{OZ}$, the following inequality holds:

$$v_{i}cos{\alpha }_h\ge \left|\lambda \right|.$$

(13)

If, in addition, the normal N to the associated isosurface is not vertical $cos{\alpha }_h\ne 0$, then

$$\begin{array}{c}{\widehat{u}}_i\sqrt{v_i^2{cos}^2{\alpha }_h-{\lambda }^2}\ge |\mathit{II}\left[V_i^{\pm }\right]+\alpha +2〈\overrightarrow{\omega };V_i^{\pm }〉|,\end{array}$$

(14)

$$\begin{array}{c}\mathrm{where}{V}_i^{\pm }:=\lambda tan{\alpha }_h{\mathit{h}}_{\tan }\pm \overrightarrow{\tau }\sqrt{v_i^2-{\lambda }^2{cos}^{-2}{\alpha }_h}\end{array}$$

(15)

and the inequality holds with any sign in ±.

This lemma is immediate from Proposition 4.1 in [27].

Inequality (13) means that projected onto the normal to AI, the speed of the ith UAV is enough to compensate for the normal displacement of AI in order to remain on the moving AI. Meanwhile, (14) means that while keeping its altitude unchanged, the UAV can remain on AI by meeting the challenges from the translational acceleration of AI and the Coriolis and centrifugal accelerations caused by the motion of the UAV over AI. We slightly enhance inequalities (13) and (14) by assuming that they hold with > put in place of ≥ and do not regress as $t\to \infty$ or (if applicable) $\Vert \mathit{r}\Vert \to \infty$.

Assumption 2.

In the operational zone $\mathcal{OZ}$, Assumption 1 and inequalities (13) and (14) hold in an enhanced form: there exist ${\Delta }_{\lambda }$, ${\Delta }_u$, $b_{\rho }>0$ such that $\Vert \nabla F(t,\mathit{r})\Vert \ge b_{\rho }^{-1}$, and for all i,

$$\begin{array}{c}{v}_{i}cos{\alpha }_h\ge \left|\lambda \right|+{\Delta }_{\lambda },\end{array}$$

(16)

$$\begin{array}{c}{\widehat{u}}_i\sqrt{v_i^2{cos}^2{\alpha }_h-{\lambda }^2}\ge |\mathit{II}\left[V_i^{\pm }\right]+\alpha +2〈\overrightarrow{\omega };V_i^{\pm }〉|+{\Delta }_u.\end{array}$$

(17)

Here (16) guarantees that $cos{\alpha }_h>0$ and so the normal N is not vertical, i.e., the prerequisite for (14) is met.

In the real world, the next assumption is typically satisfied.

Assumption 3.

In the operational zone, the basic characteristics of the field stay bounded:

$$\begin{array}{c}|g_{\rho }|\le b_g<\infty ,|n_{\rho }|\le b_n<\infty ,\Vert \nabla \nabla \rho \Vert \le b_{\nabla }<\infty ,\left|\lambda \right|\le b_{\lambda }<\infty ,\hfill \\ \hfill \parallel \overrightarrow{\omega }\Vert \le b_{\omega }<\infty ,\left|\varkappa \right|\le b_{\varkappa }<\infty \forall (t,\mathit{r})\in \mathcal{OZ}.\end{array}$$

(18)

The control objective pursued in this paper tacitly assumes that the isosurface $S_t\left(f_{\star }\right)$ can be horizontally circumnavigated and so is “horizontally” bounded. However, Assumptions 1–3 do not guarantee this. So we need another assumption.

Assumption 4.

There exists a constant $d_{\mathrm{hor}}\in (0,\infty )$ such that for any t, the distance between the horizontal projections of any two points from $S_t\left(f_{\star }\right)\cap \mathcal{OZ}$ does not exceed $d_{\mathrm{hor}}$.

Before applying the control law (9) and (10), the UAVs are to be driven to or put in special postures. Since this is trivially performable, we do not come into implementation details and merely describe those postures.

Assumption 5.

Initially, all UAVs are oriented horizontally and are in the interior of the operational zone $\mathcal{OZ}$ at distinct altitudes.

The requirement to the altitudes can be met with probability 1 if, for example, every teammate is instructed to preliminarily reach its own altitude that is independently drawn for any of them from a common continuous probability distribution.

If Assumption 5 holds and $〈{\mathit{u}}_i;\mathit{h}〉=0\forall t$, then the ith UAV does not leave its initial horizontal plane. Let, in addition, the control input ${\mathit{u}}_i$ continuously depends on time and

$$\Vert {\mathit{u}}_i\left(t\right)\Vert ={\widehat{u}}_i\forall t,$$

(19)

where ${\widehat{u}}_i$ is the upper bound from (2). Then the ith UAV moves over the boundary of one of two horizontal discs ${\mathring{\mathcal{D}}}_i^{\pm }$, making a full turn for $2\pi /{\widehat{u}}_i$ units of time. We assume that these discs lie in $\mathcal{OZ}$ and that the UAV’s turning rate exceeds the mean rate (over some initial period of time) at which the isosurface rotates about the vertical axis.

Assumption 6.

For any UAV i, there exists a natural number $n_i$ such that the following statements hold for the time interval ${\mathfrak{I}}_i=\left[0,2\pi n_i/{\widehat{u}}_i\right]$:

(i) 

During this interval, the horizontal projection $N_{\mathrm{hor}}[t,{\mathit{r}}_i\left(0\right)]$ of the unit vector $N[t,{\mathit{r}}_i\left(0\right)]$ normal to the associated isosurface rotates through an angle that does not exceed $2\pi (n_i-1)$;

(ii) 

The initial discs lie inside the operational zone (12) during the time interval${\mathfrak{I}}_i$:

$$f_{-}<F(t,\mathit{r})<f_{+}\forall t\in {\mathfrak{I}}_i,\mathit{r}\in {\mathring{\mathcal{D}}}_i^{\pm }.$$

By (16), the normal N is not vertical. Hence its horizontal projection is nonzero, so the vector $N_{\mathrm{hor}}[t,{\mathit{r}}_i\left(0\right)]$ and its rotation angle are well-defined. If the field is steady $F(\mathit{r},t)=F\left(\mathit{r}\right)$, this angle is zero, and (i) does hold with $n_i:=1$.

5. Chimerical Solutions

Under the control law (9) and (10), the closed-loop system is described by an ordinary differential equation (ODE) with a discontinuous right-hand side (RHS). In the theory of such ODE’s, studies on the phenomenon of sliding have been primarily confined to the case of attractive discontinuity surfaces up to now. Only the slightest attention has been paid to non-attractive ones. For them, the discussion has been typically brief and limited to a reference to the very possibility of sliding, with a two-side repelling surface S exemplified in Figure 3a being the most popular subject of focus.

However, there is much more diversity in sliding surfaces than mentioned above. Figure 3b,c presents two examples, where the surface of interest S contains a single (green) point and has the zero dimension. This surface hosts a sliding solution $ss$ (“staying still at S”), whereas some other solutions reach S in a finite time and can be continued by $ss$. These are those starting in the pink domain, which has a nonzero area in Figure 3c. In Figure 3a,b, the sliding solution is non-viable and can be treated as nonexistent in reality since it is catastrophically sensitive to arbitrarily small disturbances. Specifically, almost all (for a continuous probability) of them bring the state in the white domain, after which the state essentially deviates from the sliding solution on any finite time interval. This deviation is not small for small disturbances, and a nonzero lower bound on the deviation is determined by the interval. Moreover, disturbance causes an immediate repulsion from the sliding solution in Figure 3a, which also holds in Figure 3b if the state is brought to the angles A or C. If it is brought to the angles B or D, repulsion still occurs with probability 1. Still, it commences after a transient, whose duration is proportional to the disturbance magnitude. Meanwhile, Figure 3c shows that the overall diversity of repulsive behaviors is richer than those just discussed. For example, suppose all directions of disturbance are equiprobable. In that case, the disturbance brings the state to the pink domain in Figure 3c with a nonzero probability, and then the solution returns to $ss$ in a finite time. Simultaneously, the disturbance brings the state in the white domain with a nonzero probability, and then the solution diverges far away from $ss$.

The apathy of the classic theory to the detailed classification of the entire range of behaviors possible near sliding solutions partly stems from being aimed at building a controller that imparts a useful feature to the system, e.g., reduced dimension, robustness against disturbances, etc., which objective calls for an attractive sliding surface. With no intention to fill the identified gap, we offer a general concept of sliding solutions that “do not exist in reality” in the fashion similar to $ss$ in Figure 3a,b. The rationale for this selection is that such solutions can be formally found in our closed-loop system, whereas they can be ignored in a practical setting and so in the results.

With these in mind, we consider a differential inclusion (DI) $\dot{\mathit{x}}\in \mathcal{R}(t,\mathit{x})$ on a Riemannian manifold $\mathcal{M}$, where the RHS is a convex, compact, and nonempty subset of the tangential (to $\mathcal{M}$ at $\mathit{x}$) space and the map $\mathcal{R}(·)$ is upper-semicontinuous; these properties imply local solvability of the initial value problem by Theorem 6.2 in [42]. (For the swarm of the UAVs (2) driven by (9) and (10), the points of $\mathcal{M}$ are the team’s states $\mathit{x}={\{{\mathit{r}}_i,{\mathit{e}}_i\}}_{i=0}^{N-1}$ with ${\mathit{r}}_i\in {\mathbb{R}}^3$ and ${\mathit{e}}_i$ from the unit sphere in ${\mathbb{R}}^3$ centered at the origin, and $\mathcal{R}(t,\mathit{x})$ is obtained via the Fillipov’s covexification procedure [40]).

Let $\mathcal{T}$ be a (maybe infinite) time interval. A solution $\overline{\mathit{x}}\left(t\right),t\in \mathcal{T}$ of the DI is said to be fully chimerical if for any finite subinterval ${\mathcal{T}}_{*}=[{\tau }_{-},{\tau }_{+}]\subset \mathcal{T},{\tau }_{-}<{\tau }_{+}$, there is $\delta >0$ such that almost all (with respect to the Lebesgue measure or, equivalently, with respect to any continuous probability distribution) initial states ${\mathit{x}}_{\mathrm{in}}$ from a sufficiently small ball centered at $\overline{\mathit{x}}\left({\tau }_{-}\right)$ give rise to a solution $\mathit{x}\left(t\right),t\in {\mathcal{T}}_{*}$ whose maximal deviation from $\overline{\mathit{x}}(·)$ on ${\mathcal{T}}_{*}$ is no less than $\delta$, irrespective of how close ${\mathit{x}}_{\mathrm{in}}$ is to $\overline{\mathit{x}}\left({\tau }_{-}\right)$. The solution is said to be chimerical if it is fully chimerical on some subinterval of $\mathcal{T}$. In Figure 3b, the sliding solution $ss$ is fully chimerical, whereas examples of chimerical solutions are given by those that start on the pink line, then arrive at S, then stay at S, and then (possibly) depart from $ss$ either to the left or to the right. Fully chimerical and so chimerical solutions are nonexistent in practice since they are not stable against inevitable disturbances and errors in sensors, computational units, and actuators, including quantization errors.

With no intention to fully categorize the converse case, we say that the solution $\overline{\mathit{x}}(·)$ is firmly corporeal if, for any closed finite subinterval ${\mathcal{T}}_{*}\subset \mathcal{T}$, that maximum deviation goes to zero whenever $\mathit{x}\left({\tau }_{-}\right)\to \overline{\mathit{x}}\left({\tau }_{-}\right)$, and corporeal if there exists a finite set F such that the solution is firmly corporeal on every connected component of $\mathcal{T}\setminus F$. If $\mathcal{T}=[0,\infty )$, the latter means existence of times ${\tau }_1<{\tau }_2<\dots <{\tau }_s\in \mathcal{T}$ such that the solution is firmly corporeal on $[0,{\tau }_1),({\tau }_1,{\tau }_2),\dots ,({\tau }_{s-1},{\tau }_s)$, and $({\tau }_s,\infty )$. Figure 3b shows that the solutions starting in the white domain are firmly corporeal. Any initial state from the pink line gives rise to exactly two corporeal solutions: they go over this line to point S and then immediately leave it either to the left or right. The set F consists of a single time when the trajectory passes through S; at this time, the solution branches in two directions.

6. Main Results

We first skip tedious tuning details and show that the proposed navigation scheme is enough to solve the mission under minimal and partly inevitable assumptions.

Theorem 1.

Let Assumptions 1–6 hold, and the visibility range be not overly small: $d^{\mathrm{vis}}>\frac{{ɦ}_{+}-{ɦ}_{-}}{N}+d_{\mathrm{hor}}$. Then the parameters (3) of the control law and the maps χ and ð in (9) and (10) can be chosen so that the closed-loop system has a corporeal solution defined on $[0,\infty )$, and for any such solution and moreover, for any non-chimerical one the following claims are true:

(i) 

Any UAV is never vertically oriented: $sin{\theta }_i\left(t\right)\ne 0\forall t,i$;

(ii) 

The output of the control law (9) is feasible, i.e., the third and fourth relations in (2) do hold;

(iii) 

The team members do not collide with one another;

(iv) 

They are always in the pre-specified altitudinal range $H_{\mathrm{al}}$ and the operational zone (12);

(v) 

The team forms the densest horizontally sweeping net on the targeted isosurface $S_t\left(f_{\star }\right)$ in the range from ${ɦ}_{-}$ to ${ɦ}_{+}$, as is specified by Definition 1.

Moreover, let a compact set $Q_{\mathrm{in}}$ of initial states be given such that any its element satisfies Assumptions 5 and 6. Then common values of the controller parameters (including the functions χ and ð) can be chosen so that (i)–(v) hold whenever the initial state of the team is in $Q_{\mathrm{in}}$.

The proofs of all theoretical results are given in Appendix B, Appendix C and Appendix D.

Theorem 1 means that our control law ensures the attainment of the posed objective.

By Assumption 4, the requirement $d^{\mathrm{vis}}>\frac{{ɦ}_{+}-{ɦ}_{-}}{N}+d_{\mathrm{hor}}$ means that if the UAVs are close to the targeted isosurface $S_t\left(f_{\star }\right)$ and the even distribution over the altitudinal range ${\mathcal{H}}_{al}$ is nearly attained, the “altitudinally adjacent” robots “see” each other. Meanwhile, at the initial time the UAVs are permitted to be arbitrarily distributed over the range ${\mathcal{H}}_{al}$ (modulo that different UAVs should be at distinct altitudes by Assumption 5) so that some “altitudinally adjacent” robots may not “see” each other since the distance between them exceeds $d^{\mathrm{vis}}$. However, (v) in Theorem 1 and the first sentence in the current paragraph imply that this unwanted situation is eventually eliminated under the action of the proposed algorithm.

Theorem 1 neglects chimerical trajectories. According to Appendix B, Appendix C and Appendix D (see Lemmas A8, A18 and A19), such trajectories possess at least one of the following two features. (1) On some time interval $[0,\tau ],\tau >0$, the state remains on a two-side repelling (like in Figure 3a) surface described by ${\dot{f}}_i=-\chi (f_i-f_{\star })$ with some i. (2) There exist $T\ge 0$ and $i\ne j$ such that the UAVs i and j constantly remain at a common and constant altitude for $t\ge T$. Chimericality means that both features are unrealizable in the closed-loop due to instability against arbitrarily small perturbations, errors, and noises.

When discussing controller tuning, we have in mind the situation of multiple initial states from $Q_{\mathrm{in}}$. (The case of a single initial state is that of a singleton set $Q_{\mathrm{in}}$.)

Preparatory choice of ${\widehat{\mathit{u}}}_{\mathit{i}}^{\mathit{f}}$ from (9). In (19), we put some ${\widehat{u}}_i^f<{\widehat{u}}_i$ in place of ${\widehat{u}}_i$ and thus acquire larger disks ${\mathring{\mathcal{D}}}_i^{\pm }\left({\widehat{u}}_i^f\right)\supset {\mathring{\mathcal{D}}}_i^{\pm }$. If ${\widehat{u}}_i^f\approx {\widehat{u}}_i$, they are close to ${\mathring{\mathcal{D}}}_i^{\pm }$ (uniformly over all initial states from $Q_{\mathrm{in}}$) and so lie inside the operational zone for $t\in [0,2\pi n_i/{\widehat{u}}_i]$ by (ii) of Assumption 6, which thus remains true if ${\widehat{u}}_i$ is replaced by ${\widehat{u}}_i^f<{\widehat{u}}_i,{\widehat{u}}_i^f\approx {\widehat{u}}_i$. We pick such ${\widehat{u}}_i^f\in (0,{\widehat{u}}_i)$, with a view to possibly push it closer to ${\widehat{u}}_i$ afterward.

General requirements to the functions $\chi$ and ð from (9) and (10). We use continuous and piecewise smooth functions that map $\mathbb{R}$ to $\mathbb{R}$ and $[0,\infty )$ to $\mathbb{R}$, respectively, and are such that

$$\begin{array}{c}\chi \left(0\right)=0,\chi \left(z\right)<0\forall z<0,\chi \left(z\right)>0\forall z>0,\end{array}$$

(20)

$$\begin{array}{c}\overline{\chi }:=\underset{z\in \mathbb{R}}{sup}\left|\chi \left(z\right)\right|<\infty ,{\overline{\chi }}^{\prime }:=\underset{z\in \mathbb{R}}{sup}\left|{\chi }^{\prime }(z\pm )\right|<\infty ,\end{array}$$

(21)

$$\begin{array}{c}ð\left(0\right)=0,{ð}^{\prime }(t\pm )>0\forall t,ð\left(t\right)\to 1\mathrm{as}t\to \infty ,\overline{{ð}^{\prime }}:=\underset{t\ge 0}{sup}{ð}^{\prime }(t\pm )<\infty .\end{array}$$

(22)

Examples include $ð\left(t\right)=1-e^{-kt}(k>0)$ and

$$\chi \left(z\right)=aarctan\left(bz\right),\frac{az}{1+b\left|z\right|},max\{-a;min\{a;bz\}\}(a,b>0).$$

Switching ${\mathbf{\delta }}_{\mathit{f}}$ and saturation ${\Delta }_{\mathit{h}}$ thresholds from (8) and (10), respectively, and the parameter ${\mathit{d}}_{-}^{\mathrm{vis}}$ from (4) are chosen so that:

$$\begin{array}{c}0<{\delta }_f<min\{f_{+}-f_{\star };f_{\star }-f_{-}\},{\Delta }_h>\frac{{ɦ}_{+}-{ɦ}_{-}}{N},\end{array}$$

(23)

$$\begin{array}{c}{\left[d_{\mathrm{hor}}+2\frac{{\delta }_f{b}_{\rho }{v}_i^2}{{\Delta }_{\lambda }^2}\right]}^2+{\Delta }_h^2<{\left(d_{-}^{\mathrm{vis}}\right)}^2,d_{-}^{\mathrm{vis}}<d^{\mathrm{vis}}.\end{array}$$

(24)

Here $b_{\rho }$ and ${\Delta }_{\lambda }$ are taken from Assumption 2, and $d_{\mathrm{hor}}$ from Assumption 4. In (23), $min\{f_{+}-f_{\star };f_{\star }-f_{-}\}$ assesses the remoteness of the targeted isosurface $S_t\left(f_{\star }\right)$ from the borders of the operational zone (12). Thanks to the assumption introduced in the body of Theorem 1, the conditions (23) and (24) are feasible: they can be satisfied by picking ${\Delta }_h>\frac{{ɦ}_{+}-{ɦ}_{-}}{N},d_{-}^{\mathrm{vis}}$ and ${\delta }_f>0$ close enough to $\frac{{ɦ}_{+}-{ɦ}_{-}}{N},d^{\mathrm{vis}}$ and 0, respectively.

An auxiliary parameter$\eta \in (0,1)$ is recommended to be picked at this tuning stage to simplify the subsequent choice of the basic parameters.

The slope $\gamma$ of $\mathsf{\Xi }(·)$ in (10), the parameters ${\widehat{\mathit{u}}}_{\mathit{i}}^{\mathit{f}},{\widehat{\mathit{u}}}^{\mathit{h}}$ in (9), the upper bounds $\overline{\chi },{\overline{\chi }}^{\prime }$ in (21), and the upper bound $\overline{{ð}^{\prime }}$ in (22) are subjected to the following constraints:

$$\begin{array}{c}\gamma {\Delta }_h<v_i\forall i,{\chi }_{\rho }:=b_{\rho }\overline{\chi }<{\Delta }_{\lambda },\end{array}$$

(25)

$$\begin{array}{c}\xi :={(\gamma {\Delta }_h+{\chi }_{\rho }+b_{\lambda })}^2-b_{\lambda }^2<(1-{\eta }^2){\Delta }_{\lambda }^2;\\ \frac{v_i{\widehat{u}}^h}{\sqrt{v_i^2-{\left(\gamma {\Delta }_h\right)}^2}}+v_i[{\widehat{u}}_i-{\widehat{u}}_i^f]+2\left[v_i{b}_{\varkappa }+b_{\omega }\right]\frac{v_i}{{\Delta }_{\lambda }}\sqrt{\frac{{\xi }^2}{4{\eta }^2{\Delta }_{\lambda }^2}+{({\chi }_{\rho }+\gamma {\Delta }_h)}^2}\end{array}$$

(26)

$$\begin{array}{c}+{\chi }_{\rho }({\chi }_{\rho }{b}_n+2v_i{b}_{\nabla }+2b_g+{\overline{\chi }}^{\prime }/b_{\rho })+\frac{\xi {\widehat{u}}_i}{2\eta {\Delta }_{\lambda }}<{\Delta }_u,\end{array}$$

(27)

$$\begin{array}{c}{\left({\widehat{u}}_i^f\right)}^2+\frac{{\left({\widehat{u}}^h\right)}^2}{v_i^2-{\gamma }^2{\Delta }_h^2}\le {\left({\widehat{u}}_i\right)}^2\forall i,\end{array}$$

(28)

$$\begin{array}{c}4{\gamma }^2{\Delta }_h+2{\overline{ð}}^{\prime }\gamma {\Delta }_h<{\widehat{u}}^h\forall i.\end{array}$$

(29)

Here ${\Delta }_{\lambda },{\Delta }_u,b_{\rho }$ are taken from Assumption 2, $b_{\varkappa },b_{\omega },b_n,b_g,b_{\nabla }$ from Assumption 3, and $v_i,{\widehat{u}}_i$ from (2). At least, the requirement (26) to $\gamma$ and ${\chi }_{\rho }$ means that its left-hand side is less than ${\Delta }_{\lambda }^2$. If this is satisfied, (26) gives an upper bound on the choice of the auxiliary parameter $\eta$. Putting this bound to (27) in place of $\eta$ results in an “$\eta$-free” form of the conditions on the controller parameters, whose format is, however, rather cumbersome and so is not user-friendly. This is the rationale for using $\eta$.

Inequalities (25)–(29) are feasible. Indeed, it suffices to note that the left-hand sides of (25)–(27) go to 0 as $\gamma \to 0+,\overline{\chi }\to 0+,{\overline{\chi }}^{\prime }\to 0+,{\widehat{u}}_i^f\to {\widehat{u}}_i-,{\widehat{u}}^h\to 0+$, whereas the RHS’s are positive constants. So to meet (25)–(27), it suffices to pick $\gamma ,\overline{\chi },{\overline{\chi }}^{\prime },{\widehat{u}}^h$ small enough and ${\widehat{u}}_i^f$ close enough to ${\widehat{u}}_i$. Then (28) can be ensured by further decreasing $u^h$. After this, (29) can be satisfied (while not violating (25)–(28)) by further decreasing $\gamma$.

These considerations give guidelines for experimentally tuning the control law.

The final choice of the functions $\chi$ and ð from (9) and (10): They are chosen subject to the above general requirements and the already chosen upper bounds $\overline{\chi },{\overline{\chi }}^{\prime },{\overline{ð}}^{\prime }$.

The time ${\mathit{T}}_{\mathrm{tr}}$ from (8) is chosen so that $T_{\mathrm{tr}}>2\pi n_i/{\widehat{u}}_i^f\forall i$, where $n_i$ is taken from Assumption 6.

Theorem 2.

The statement of Theorem 1 is true if the parameters (3) of the control law and the maps $\chi (·)$ and $ð(·)$ in (9) and (10) are chosen based on the above recommendations.

So these recommendations can be used for analytically tuning the controller if the involved bounds on the parameters of the field or their estimates are available.

Remark 2.

Implementing the above recommendations results in controller parameters common for all robots, with the only exception of ${\widehat{u}}_i^f$. For homogeneous robotic teams (but not only for them), ${\widehat{u}}_i^f$’s can also be chosen common since $v_i$ and ${\widehat{u}}_i$ do not depend on i.

Remark 3.

Under the assumptions of Theorem 2, the statement (i) of Theorem 1 can be specified: the pitch angle of the ith UAV never exceeds $arcsin\frac{\gamma {\Delta }_h}{v_i}$ in absolute value. (Here $arcsin\frac{\gamma {\Delta }_h}{v_i}$ is well-defined due to the first inequality in (25).) Meanwhile, the above recommendations on the choice of the controller parameters are not violated by decreasing the coefficient $\gamma >0$. By using this, the controller can be tuned so that not only the statement of Theorem 1 is true but also the pitch angle of every UAV is always within a given bound, which can be chosen as small as desired. This observation is of interest whenever large pitch angles are unwelcome, unacceptable, or challenge the employed model (2).

7. Computer Simulation Tests

The numerical values of the basic parameters used in the tests are as follows:

$N=10$	$v_i=10.0\frac{\mathrm{m}}{\mathrm{s}}$	${\widehat{u}}_i=1.0\frac{\mathrm{rad}}{{\mathrm{s}}^2}$	${\overline{u}}^h=0.01\frac{\mathrm{rad}}{{\mathrm{s}}^2}$	${\overline{u}}_i^f=\sqrt{{\widehat{u}}_i^2-{\left({\widehat{u}}^h\right)}^2}$,
${\delta }_f=0.1$ ❢	$T_{\mathrm{tr}}=10.0$ s	$\gamma =0.05\frac{1}{\mathrm{s}}$	${\Delta }_h=40.0$ m	$f_{\star }=10$ ❢
$d^{\mathrm{vis}}=70$ m	$d_{-}^{\mathrm{vis}}=69$ m	$\chi \left(f\right)=\left\{\begin{array}{cc}\mu f/\delta \hfill & \mathrm{if}\left|f\right|\le \delta ,\hfill \\ \mathbf{sgn}\left(f\right)\mu \hfill & \mathrm{otherwise}\hfill \end{array}\right.$
$\delta =5.0$ ❢	$\mu =0.4\frac{❢}{\mathrm{s}}$	$ð\left(t\right)=1-e^{-10.0t}$

Here ❢ is the unit of measurement of the field value f. Zero-mean Gaussian additive noises corrupt the measurements of this value and the altitude h with the standard deviation of $0.1$ ❢ and $0.1$ m, respectively. No noise reduction techniques were applied, and the simplest two-point Newton’s quotient $\left[f\right(t)-f(t-\tau \left)\right]/\tau$ was employed to assess the time derivative $\dot{f}\left(t\right)$ in (9). The simulations were performed in MATLAB. In Figure 4, Figure 5, Figure 6 and Figure 7, the robots, their paths, and the targeted isosurface are depicted in green, blue, and gray, respectively; obsolete parts of the paths are erased; the targeted isosurface is treated as opaque. In fact, what is depicted is not the targeted isosurface but a close one; otherwise, the UAV’s path would seem overly dashed due to invisible portions that appear whenever the UAV dives, even slightly, behind the isosurface. Multimedia of all tests are available at https://cutt.ly/fTZsmjC, (accessed on 20 October 2021).

Figure 4 illustrates an experiment where an environmental field slowly translates along the x-axis so that its isosurfaces retain their form and size. Among the purposes of this experiment, there is complementing Theorems 1 and 2 via testing the capacity of the control law to cope with the non-smoothness of the field and isosurface. Specifically, the cross-like isosurface from Figure 4a is not smooth (and so Assumption 1 violates) on the red curves, where (except for the yellow point) the field gradient abruptly changes when crossing the curve. Initially, ten UAVs are organized into two groups of five; each group is aligned vertically and evenly distributed. Meanwhile, their vertical spacing is far from the desired one (which is associated with the even deployment from ${ɦ}_{-}=10$ m to ${ɦ}_{-}=190$ m), and the groups are out of “eye contact” with each other.

By Figure 4a, the moment of $t=20$ s can be viewed as when all UAVs localize the isosurface, though with a small degree of approximation, as can be seen in Figure 4h. Figure 4i shows that at this moment, the UAVs over-populate the altitudinal range from 40 m to 60 m, which condition is far from the targeted even distribution from 10 m to 190 m. Approximately at this time, all UAVs pass to mode $\mathfrak{M}$ and so regulation of the altitudes towards the even distribution is commenced according to (10). By Figure 4d,i, this goal is attained from $t=150$ s. The outbursts of the field value errors for two UAVs at ≈$40$ s occur when these UAVs have to pass from encircling the horizontal “beam” of the cross to dealing with the vertical one, as can be seen in Figure 4c. The height regulation module initially drives them upwards so intensively that they find themselves far enough from the targeted isosurface represented by its vertical beam, which is fairly distant from the just traced tip of the horizontal beam. These outbursts are promptly fixed and never repeated, as shown by Figure 4h. So the discussed episode can be related to the fact that by $t=40$ s, an even distribution of the altitudes has not yet been achieved, as shown in Figure 4i. Overall, the control law ensures the attainment of the control objective despite sensor errors and the non-smoothness of the field.

The experiment in Figure 5 complements Theorems 1 and 2 from another standpoint: the field gradient is not vertical at the red points in Figure 5a and their antipodes with respect to the centers of the holes. This means violation of (16) in Assumption 2 and implies that the number K of the connected components of the horizontal section varies as the cutting horizontal plane runs over the altitudinal range of interest (from 30 m to 170 m); four sample sections A, B, C, and D with different K’s are depicted in light blue in Figure 5a. The field and the targeted isosurface rotate about the pink axis. As a result, the number K varies over time for some fixed altitudes, as illustrated in Figure 5h,i. For example, the red section in Figure 5i has three components, unlike Figure 5a, with no more than two components. When keeping both the field value and the altitude constant, any UAV can trace only a single connected component and so has to “select” it from a variety of those (if applicable), e.g., has to select one of the two loops constituting B in Figure 5a. From time to time, the UAVs are forced to “reselect” because of changes in circumstances, e.g., alterations in K. Another trouble is highlighted by considering the horizontal cutting plane that gives rise to B in Figure 5a. As this plane approaches the upper red point from above, the curvature of both B-loops near that point increases without limits. It so exceeds, sooner or later, the maximal turning capacity of the UAVs. Hence, there are horizontal sections that can be traced by no means due to the limited turning radius of the UAVs. Though the proposed algorithm is not intended to handle the described troublesome issues, the experiment in Figure 5 is aimed to form an initial pre-judgment on the intrinsic potential of this algorithm for reasonably treating them.

Figure 5j,k show that the above extra challenges do not visibly worsen the performance of the control law with respect to the primary goals of finding and sweeping the isosurface and even self-distribution of the team over the altitudinal range. Moreover, Figure 5b–g provides evidence that the algorithm manages to attend all “simple” components of the topologically complex surface: two UAVs circumnavigate the “half-donut” B1 from Figure 5f, one UAV encircles B2, another UAV encircles C1, and two more UAVs circumnavigate the “half-donut” C2, whereas the remaining four UAVs go around the central part of the eight-shaped surface. This trait of “attending all parts” may also be identified at the other stages of the experiment, albeit to various degrees. In Figure 5j, the splash of the field error for a UAV at ≈$550$ s is due to being too close to a point with an excessive “curvature demand”, as is described in the previous paragraph. However, this splash is promptly fixed and never repeats within the duration of the experiment. Overall, this experiment shows that the algorithm more or less satisfactorily copes with the above extra challenges.

Figure 6 is concerned with an experiment whose purpose is to test the algorithm’s robustness against failures of the team members and its performance when dealing with a deforming isosurface. This isosurface has the form of a curved tube. The tube performs oscillatory displacements along the x-axis, alternates increasing and decreasing in size, and reshapes, e.g., changes the number of the “waves on the surface”, becomes a right cylinder or a “bottle” at some times, etc. The number of the UAVs is increased up to 20; the targeted altitudinal range is from 0 m to 200 m. Starting from the initial deployment shown in Figure 6a, all UAVs individually reach the targeted isosurface as early as at ≈$20$ s, according to Figure 6b,i. By Figure 6c,j, an even self-distribution over the altitudes is achieved later at ≈$220$ s. Then at $t=250$ s and $t=500$ s, a group of five UAVs is withdrawn from the mission (and their color is changed from green to red), whereas the remaining peers continue to run the algorithm “as usual” with ignoring the missed members. As can be inferred from Figure 6d,e,j, these peers need only ≈$50$ s to rebuild the even distribution with lesser team size. For the missing episode at $t=500$ s, this entails a slight temporal impairment in the field tracking performance, which is fixed for ≈$40$ s, as shown in Figure 6i.

Overall, the algorithm exhibits robustness to failures of the team members in a sophisticated scenario with a deforming and moving isosurface.

The last experiment tests the capacity of the algorithm to automatically manage the admission of new team members (newcomers). The deforming isosurface from the previous test is handled, though without displacement along the x-axis. The team consists of 10 members initially. Five extra members appear on stage at $t=250$ s and then another five at $t=500$ s, as shown in Figure 7. Meanwhile, the algorithm is run “as usual” at both newcomers and “oldies”, taking into account the UAVs currently present at the stage.

As follows from Figure 7f, the UAVs autonomously rebuild an even distribution over the altitudes for ≈$45$ s in both events of admission. By Figure 7e, the second admission implies detrimental effects in terms of the field value. However, they are minor in value and are overcome for ≈$150$ s. This demonstrates the algorithm’s capacity to incorporate extra UAVs in the team on the fly.

8. Conclusions and Future Work

This study aimed to design and analyze a distributed navigation and collision avoidance strategy for a team of UAVs traveling in a 3D environment. The strategy enables the team to first find the isosurface where an unknown and unpredictably varying scalar field assumes a given value and then form the vertically densest net-like barrier horizontally sweeping the isosurface. Among the complicating factors was the lack of access to the field gradient, absence of communication facilities, non-holonomy, under-actuation, and a finite control range of the UAVs. It was shown that even in such circumstances, the mission could be solved by a computationally inexpensive strategy justified by a mathematically rigorous global convergence theorem. Computer simulation tests confirmed the convergence and performance of the algorithm.

The algorithm is individually executed by each UAV and consists of two stages (operating modes). The main objective of the first and second stages is to find and arrive at the isosurface and, respectively, to track and circumnavigate it while distributing the team into the vertically densest net. The proposed regulation rule conforms to the sliding mode control paradigm at any stage. This paradigm has attracted significant interest from industry and academia thanks to well-known benefits such as high insensitivity to disturbances and noises, robustness against uncertainties, good dynamic response, and simple implementation (we refer the reader to [43,44,45,46,47,48,49] for a survey). The major problem with the practical implementation of sliding-mode controllers is identified as the possibility of a chattering phenomenon. The ever-increasing popularity of the sliding-mode approach to motion control is partly due to the development of rather effective general techniques of chattering elimination and suppression; see, e.g., [45,50,51] for a survey. Among them, there is a smooth approximation of the discontinuous controller using low-pass filters, adaptive controllers, and higher-order sliding modes. Whether the harmful chattering is encountered when implementing the proposed controller, it can be subjected to treatment via these methods. Their practical effectiveness has been widely reported, whereas the phenomenon does not necessarily occur in experiments with real mobile robots. Some examples are reported in, e.g., [46,52], where control laws that are similar in some respects to the law proposed in this paper are considered.

A fairly common approach to the design of control systems implements the idea of a two-level hierarchical structure, where a kinematic-level controller generates a reference signal to be tracked by low-level controllers. The findings of this paper are concerned with the first stage and are based on a model of the UAV’s kinematics. Implementation issues concerned with the second stage and controllers somewhat similar to that from this paper are addressed in, e.g., [46,53].

Future work includes an extension of the findings of this paper to the case where along with all the previous control goals, the UAVs should be “horizontally synchronized” in some sense, e.g., should all be ultimately contained by a common vertical moving plane. Consideration of more sophisticated isosurfaces and models of UAVs is also on the agenda.