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

: 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 ﬁeld. 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 ﬁeld assumes a given value. They also should evenly distribute themselves over a pre-speciﬁed 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 ﬁeld 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 ﬁnite 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 justiﬁed by a mathematically rigorous global convergence result; computer simulation tests conﬁrm its performance.


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. Sections 2 and 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 Appendices A-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 ·; · , · , and × stand for the standard inner product, Euclidean norm, and cross product in R 3 , respectively.

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, r) ∈ R of time t and space location r ∈ R 3 . All UAVs should arrive at the locus of points (called isosurface) S t ( f ) := {r : F(t, r) = f } with the pre-specified field value f . 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 h ∈ R 3 ; the associated coordinate h(r) of point r ∈ R 3 is loosely referred to as altitude. The mission is confined to a certain altitudinal range H al = [h − , h + ], h − < h + ; 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 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 (t) := F[t, r i (t)] at its own current location r i = r i (t), and has no idea about the distance to or bearing of the targeted isosurface S t ( f ). The ith UAV also has access to 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 h ± of the altitudinal range H al , that lie within a given "visibility" range d 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 H al , we denote by 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 is the altitude of the ith UAV.
The imposed information constraints mean that any UAV is unaware of the altitude h 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]: Here u i is the control, the upper bound u i > 0 on its magnitude is given, v i > 0 is the constant speed of the UAV, 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 e i constant. The model (2) captures the robot's capacity to travel over space paths whose curvature radii ≥ v i / 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 u i ↔ (q i , r i ) in Remark 2.1 in [31], the vector control input 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].

Proposed Hybrid Controller
It uses the following tunable free parameters and two functions χ and ð, where χ maps R to R and ð maps [0, ∞) 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 (t) of its essential neighbors by executing the following instructions: 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 vis − < d 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 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 higherh + i and lower h − i altitudinal gaps near the ith UAV are defined to be Here either + should be put in place of ± everywhere, or the same should be performed with − instead of +. If E ± i = ∅ and d vis > |h i − h ± |, there are no essential neighbors between the UAV at hands and the top/bottom of the altitudinal range H al ; thenh ± i is twice the altitudinal distance to the top/bottom. Given a UAV, the gapsh ± i 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 u i must lie in the plane e ⊥ i normal to the centerline vector e i (the pitch-yaw plane). We define u i in a special orthonormal basis h where h py i is the orthogonal projection of the vertical vector h onto the pitch-yaw plane e ⊥ i , which projection is then normalized to the unit length: Here θ i is the angle between the centerline e i and the vertical h. As will be shown in (i) of Theorem 1, our controller ensures that any UAV is always non-vertical, so both the vector h py i and the considered basis are well-defined. Our control law is hybrid with two modes: A (approaching the isosurface) and M (main mode). Any UAV i starts in A and then switches to M. The role of mode 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 by not more than a pre-specified and nominally small δ f from (3). Vertical distribution of the team is the task of the next mode M so that the global search of the isosurface during mode A is not disturbed by control signals aimed for other purposes. The switching rule is as follows: The control rule invokes (like (8) does) parameters from (3) and is as follows: where v i := 0 in mode A, Here χ and ð are functions that are to be chosen by the designer of the controller (see (20)- (29) for details), T sw i is the time of the transition A → M in (8), and is a linear function of L ≥ 0 with saturation at the threshold ∆ h > 0 mentioned in (3).
The ith UAV can compute the derivativeḣ i = v i h; e i since it is aware of its own speed v i and has access to the vertical vector h in its local frame illustrated in Figure 1. Numerical differentiation can be used to assess the time-derivativeḟ i of the sensor readings f i (t). 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 h in the local frame of reference of the UAV at hand. It thus provides access to the angle θ i between h and the centerline 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 M and A, respectively.  Except for the coefficient u f i , the control rule (9) is common for all robots. Remark 2 discusses when complete uniformity can be achieved by picking u f i common for all i.

Remark 1.
To facilitate understanding the procedure for determiningh ± i , its first step was pictured as building the sets E i , E ± i 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.

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: • ∇F, spatial gradient of the field; • N(t, r) = ∇F(t,r) ∇F(t,r) , unit vector normal to the associated isosurface (AI) that passes through location r at time t; • α h = arcsin N; h , angle from this vector N to the horizontal planes; κ, maximal (in absolute value) eigenvalue of this quadratic form; • λ(t, r), front velocity of the associated isosurface; • α(t, r), front acceleration of the associated isosurface; • ω(t, r), angular velocity of rotation of the associated isosurface; • ρ(t, 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 ρ (t, r), proportional growth rate of the density ρ with time; • ∇ ∇ρ(t, r), proportional tangential gradient of the density ρ, • n ρ (t, r), proportional growth rate of the density ρ 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 hor t ( f |h) = {r ∈ S t ( f ) : h(r) = h} at a fixed altitude h = h i ∈ H al . Since the size N of the team and so the altitudes h i 's are unknown to the team members and no UAV is assigned its own altitude h i a priori, it is fair to require that any UAV be able to trace the whole of S hor t ( f |h) within the operational zone OZ at any altitude h ∈ H al in the given altitudinal range 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 OZ in terms of the extreme values f − < f + ( f ∈ ( f − , f + )) achievable by the field F in this zone: 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 hor t ( f |h), h ∈ 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 ∇F = 0, and the horizontal section S hor t ( f |h) of the f -isosurface at any altitude h ∈ H al is not empty.
Conditions necessary for the mission feasibility are as follows. Lemma 1. Let the isosurface S t ( f ) be trackable by the ith UAV. Then at any temporal-spatial point of S t ( f ) ∩ OZ, the following inequality holds: v i cos α h ≥ |λ|. (13) If, in addition, the normal N to the associated isosurface is not vertical cos α h = 0, then where 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 → ∞ or (if applicable) r → ∞.

Assumption 2.
In the operational zone OZ, Assumption 1 and inequalities (13) and (14) hold in an enhanced form: there exist ∆ λ , ∆ u , b ρ > 0 such that ∇F(t, r) ≥ b −1 ρ , and for all i, Here (16) guarantees that cos α 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: The control objective pursued in this paper tacitly assumes that the isosurface S t ( f ) 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 hor ∈ (0, ∞) such that for any t, the distance between the horizontal projections of any two points from S t ( f ) ∩ OZ does not exceed d 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 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 u i ; h = 0 ∀t, then the ith UAV does not leave its initial horizontal plane. Let, in addition, the control input u i continuously depends on time and where u i is the upper bound from (2). Then the ith UAV moves over the boundary of one of two horizontal discsD ± i , making a full turn for 2π/ u i units of time. We assume that these discs lie in 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 normal to the associated isosurface rotates through an angle that does not exceed 2π(n i − 1); (ii) The initial discs lie inside the operational zone (12) during the time interval I i : By (16), the normal N is not vertical. Hence its horizontal projection is nonzero, so the vector N hor [t, r i (0)] and its rotation angle are well-defined. If the field is steady F(r, t) = F(r), this angle is zero, and (i) does hold with n i := 1.

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)ẋ ∈ R(t, x) on a Riemannian manifold M, where the RHS is a convex, compact, and nonempty subset of the tangential (to M at x) space and the map 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) with r i ∈ R 3 and e i from the unit sphere in R 3 centered at the origin, and R(t, x) is obtained via the Fillipov's covexification procedure [40]).
Let T be a (maybe infinite) time interval. A solution x(t), t ∈ T of the DI is said to be fully chimerical if for any finite subinterval T * = [τ − , τ + ] ⊂ T, τ − < τ + , there is δ > 0 such that almost all (with respect to the Lebesgue measure or, equivalently, with respect to any continuous probability distribution) initial states x in from a sufficiently small ball centered at x(τ − ) give rise to a solution x(t), t ∈ T * whose maximal deviation from x(·) on T * is no less than δ, irrespective of how close x in is to x(τ − ). The solution is said to be chimerical if it is fully chimerical on some subinterval of 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 x(·) is firmly corporeal if, for any closed finite subinterval T * ⊂ T, that maximum deviation goes to zero whenever x(τ − ) → x(τ − ), and corporeal if there exists a finite set F such that the solution is firmly corporeal on every connected component of T \ F. If T = [0, ∞), the latter means existence of times τ 1 < τ 2 < . . . < τ s ∈ T such that the solution is firmly corporeal on [0, τ 1 ), (τ 1 , τ 2 ), . . . , (τ s−1 , τ s ), and (τ s , ∞). 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.

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: 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, ∞), 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 θ i (t) = 0 ∀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 al and the operational zone (12); (v) The team forms the densest horizontally sweeping net on the targeted isosurface S t ( f ) in the range from h − to h + , as is specified by Definition 1.
Moreover, let a compact set Q 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 in .
The proofs of all theoretical results are given in Appendices B-D. Theorem 1 means that our control law ensures the attainment of the posed objective. By Assumption 4, the requirement d vis > h + −h − N + d hor means that if the UAVs are close to the targeted isosurface S t ( f ) and the even distribution over the altitudinal range 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 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 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 Appendices B-D (see Lemmas A8, A18 and A19), such trajectories possess at least one of the following two features. (1) On some time interval [0, τ], τ > 0, the state remains on a two-side repelling (like in Figure 3a (2) There exist T ≥ 0 and i = j such that the UAVs i and j constantly remain at a common and constant altitude for t ≥ 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 in . (The case of a single initial state is that of a singleton set Q in .) Preparatory choice of u f i from (9). In (19), we put some u , with a view to possibly push it closer to u i afterward.
General requirements to the functions χ and ð ð ð from (9) and (10). We use continuous and piecewise smooth functions that map R to R and [0, ∞) to R, respectively, and are such that Examples include ð(t) = 1 − e −kt (k > 0) and Switching δ f and saturation ∆ h thresholds from (8) and (10), respectively, and the parameter d vis − from (4) are chosen so that: Here b ρ and ∆ λ are taken from Assumption 2, and d hor from Assumption 4. In (23), min{ f + − f ; f − f − } assesses the remoteness of the targeted isosurface S t ( f ) 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 respectively. An auxiliary parameter η ∈ (0, 1) is recommended to be picked at this tuning stage to simplify the subsequent choice of the basic parameters.
The slope γ of Ξ(·) in (10), the parameters u f i , u h in (9), the upper bounds χ, χ in (21), and the upper bound ð in (22) are subjected to the following constraints: Here ∆ λ , ∆ u , b ρ are taken from Assumption 2, b κ , b ω , b n , b g , b ∇ from Assumption 3, and v i , u i from (2). At least, the requirement (26) to γ and χ ρ means that its left-hand side is less than ∆ 2 λ . If this is satisfied, (26) gives an upper bound on the choice of the auxiliary parameter η. Putting this bound to (27) in place of η results in an "η-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 η.
These considerations give guidelines for experimentally tuning the control law.  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 γ > 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).

Computer Simulation Tests
The numerical values of the basic parameters used in the tests are as follows: Here f 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 f and 0.1 m, respectively. No noise reduction techniques were applied, and the simplest two-point Newton's quotient [ f (t) − f (t − τ)]/τ was employed to assess the time derivativeḟ (t) in (9). The simulations were performed in MATLAB. In Figures 4-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 6 January 2022). 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 h − = 10 m to h − = 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 M and so regulation of the altitudes towards the even distribution is commenced according to (10). By Figures 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.  Figures 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 Figures 6b,i. By Figures 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 Figures 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.

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.

Acknowledgments:
We acknowledge the useful assessments and corrections from the anonymous reviewers, as well as the Journal Editors.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

UAV Unmanned Aerial Vehicle AI Associated Isosurface ODE Ordinary Differential Equation RHS Right-Hand Side DI
Differential Inclusion a.a.
"for almost all", i.e., for all but a set of the zero Lebesgue measure

Appendix A. Characteristics of the Field
This section offers rigorous definitions of field characteristics used in the theoretical results and technical facts underlying the proofs of these results. We start with the former.
• r + (∆t|t, r) and r +, f (∆ f |t, r), nearest (to r) point where the axis drawn from r in the direction of the normal N to the AI intersects the time-and space-displaced S t+∆t [ f † ] and S t ( f † + ∆ f ) isosurface, respectively, where f † := F(t, r); • p(∆t|t, r) and q(∆ f |t, r), coordinates of these respective points along that axis; • λ(t, r), front velocity of the isosurface, i.e., lim ∆t→0 ∆t −1 p(∆t|t, r); • α(t, r), front acceleration of the isosurface: • ω(t, r), angular velocity of rotation of the isosurface: , r), density of the isosurfaces: ρ(t, r) := lim ∆ f →0 ∆ f /q(∆ f |t, r); • g ρ (t, r), proportional growth rate of this density with time: • n ρ (t, r), normal proportional growth rate of the density: n ρ (t, r) := 1 • ∇ ∇ρ(t, r), tangential proportional gradient of the density, i.e., the tangential (to the AI) vector such that for any tangential vector V, shape operator, where D V N is the derivative of the vector-field N in the direction V tangential to the isosurface in Section 4 in [41].
Informal comments on these definitions are available in [54]. In this section, we adopt Assumptions 1 and 2. Then the above quantities are well-defined in the operational zone [54].
From now on, the notation A C = B means that the equation holds by the fact stated or referenced above =; here = can be replaced by the symbol of any binary relation, e.g., >, ≤, etc. The symbol ⇒ means "implies that"; |S| stands for the number of elements in the set S. The arguments of the form (t), (t, r), and [t, r i (t)] are omitted whenever this does not confuse.
The first two lemmas offer technical formulas concerned with the mechanics of environmental fields and the motion of UAVs in them.

Lemma A1 ([54]
). The following relations hold in the operational zone: ρ(t, r) = ∇F(t, r) , N = ∇F(t, r)/ ∇F(t, r) , λ(t, r) = −F t (t, r) ∇F(t, r) . For any vector V tangent to the isosurface, we have Lemma A2. Whenever the ith UAV moves in the operational zone, the following relations hold: Proof. The first two formulas in (A8) hold by (2); the third one is true sincė As is noted after Assumption 2, the vectors h and N are not co-linear. Hence h, N and τ = h × N/ cos α h form a basis in R 3 and so e i = xh + yN + z τ. By finding x, y, z from the first and third formulas in (A8), with using the equations e i = 1, h; N = sin α h , we arrive at (A9) since To prove (A12), we observe that by (2) and (A9), Since N = 1 everywhere, the derivatives (A2) and S(V i ) = −D V i N are perpendicular to N; so is (A10) by the definition of τ and h tan , and ∇ ∇ρ by its own definition. By combining this with (2), (A12), the third equation from (A8), and the formula II[V] = S(V); V in Section 4 in [41], we infer thatf The next lemma offers useful facts about the unit vectors (7), τ, and h tan .
Lemma A3. If the ith UAV moves in the operational zone and is not vertical, the following formulas hold: where the sign in ∓ is borrowed from (A10).
Proof. By (7), h py i × e i = h×e i sin θ i , so the first formula in (A15) is true, and also h py hence the second formula is valid as well. By using the triple product [A, B, C] = A; B × C ∀A, B, C ∈ R 3 , we see that Here (a) and (b) use the definition τ = h × N/ cos α h of τ, and (c) uses the definition h tan = (h − N sin α h )/ cos α h of h tan , by which τ; h tan = 0. Thus the third formula in (A15) is true.
The last lemma in this appendix displays a technical property of special systems of ODEs.

Appendix B. Proofs of the Results from Section: Convergence to the Isosurface
From now on, we assume that the UAVs are driven by the proposed control law, and the assumptions of Theorem 2 hold. By Assumption 5, the ith UAV is initially inside the operational zone OZ defined by (12). Lemma A5. While the ith UAV remains in the operational zone (12), the following is true: (i) While this UAV is in the initial mode A, it flies at a constant altitude from (h − , h + ), andḣ i = 0; (ii) If A → M at some time T sw i , then afterward |ḣ i | ≤ γ∆ h and the UAV is not vertically oriented; (iii) The ith UAV is never vertically oriented and | cos θ i | ≤ γ∆ h /v i ; (iv) The output of (9) meets the requirements from (2). (10) and soḧ iḣi < 0 near the discontinuity surface {(r i , e i ) :ḣ i = 0}, which is thereby sliding. Sinceḣ i (0) = 0 by Assumption 5 and (A8), sliding motion over this surface commences at t = 0 and continues until either A is switched off or robot i leaves OZ. So until this moment,ḣ i = 0. The constant altitude of the UAV belongs to (h − , h + ) by Assumption 5. (ii) By (10), (11) and (22)

Proof. Whenever the UAV is not vertical, we have by invoking the first two formulas in (A15):
i , τ , (A24) holds by the foregoing and soḣ i > γ∆ h ⇒ḧ i < 0 andḣ i < −γ∆ h ⇒ḧ i > 0. It follows that on this interval, {(r i , e i ) : −γ∆ h ≤ḣ i ≤ γ∆ h } is a trapping region, so the UAV remains there, and |ḣ i (t)| ≤ γ∆ h , where γ∆ h < v i by (25). Then the first equation in (A8) yields that | e i ; h | < 1, so the UAV is not vertical. If τ < T op i , then by letting t → τ− in |ḣ i (t)| ≤ γ∆ h < v i , we see that the UAV is not vertical at t = τ and, by the continuity argument, for all t ∈ T sw i , τ 1 , where τ 1 ∈ (τ, T op i ) is close enough to τ. However, this contradicts the definition of T sw i , τ as the connected component of T. Thus τ = T op i , which completes the proof of (ii).
(iii) summarizes some parts of (i) and (ii) modulo the first equation from (A8). (iv) Since in (9), the unit vectors (7) and h py i × e i are perpendicular to e i , the output u i of (9) meets the third requirement from (2). Since these vectors are mutually perpendicular, Thus we see that the last requirement from (2) is satisfied.
Since the ith UAV starts within the operational zone (12) by
is implied by (i) and (A10). (iii) We invoke S i from (A26) and will examine the limit points ofṠ i andf i when approaching a state s from S op i in such a way that S i = 0 and sgn S i is kept unchanged. We will retain the notationṡ S i andf i for these limit values, sgn S i for the constant value of the sign, and will compute the other quantities at the state s. Due to (9) and (A13), we see thaṫ Here ζ ∈ [−1, 1] is born by sgn ḣ i − v i in (9) via Filippov's convexification procedure [40].
By using (iii) in Lemma A5, we see that | sin | ≤ 1 since both N and h py i are unit vectors. Now we put V i := v 2 i cos 2 α h − λ 2 and Υ i := u f i V i ±Ṡ i sin θ i sgn S i , and note that the formula forṠ i can be rewritten in the form: To estimate A, B, and C, we first note that ≤ χ ρ (χ ρ b n + 2v i b ∇ + 2b g ).
By invoking (15), (A8), and (A10), we see that Now we introduce the function f (y) := v 2 cos 2 α h − λ 2 + y and note that V i = f (0), whereas V τ i = f (x) by (A28). Due to (A30), | f (y)| ≤ 1/(2η∆ λ ) for y between 0 and x. By applying the mean value theorem to the function f (·), we see that (A34) By (15) and the last equation from (A15), By retaining the notation II for the bilinear form II [V, W] associated with the quadratic form II [V] and noting that |II [V, W]| ≤ |κ| V W for the maximal (in absolute value) eigenvalue κ of the latter form, we see that By summarizing and invoking (A32), we see that So ∓Ṡ i S i > 0 by (A31), which completes the proof of (iii).
The next lemma employs the natural number n i from Assumption 6.
Lemma A8. (i) There exists a closed-loop trajectory such that the following claims hold for any UAV i:

1.
There exists time t sm i < T i := 2πn i / u f i such that during (0, t sm i ), the UAV is in mode A and moves at a constant altitude with S i = 0, and f − ≤ f i (t) ≤ f + ; 2.
For t ≥ t sm i , the UAV undergoes sliding motion over S We intend to show that there is t sm i for which 1 is true and the UAV is on S op i,− at t = t sm i . In the case A), it suffices to take t sm i := 0 (then (0, t sm i ) = ∅). In the case B), we focus on the situation where S i (0) > 0 ∀t > 0, t ≈ 0, the converse one S i (0) < 0 ∀t > 0, t ≈ 0 is handled likewise. Suppose that t sm i does not exist in the case B). Then (8) and the specification of T tr given just before Theorem 2, during [0, T i ] mode A is on andḣ i = 0, h i ∈ (h − , h + ). Also, S i (t) > 0 ∀t ∈ (0, τ] since otherwise S i arrives at 0 at a time τ ar ∈ (0, τ] and t sm i := τ ar the UAV is on S  (2) and (9), for t ∈ [0, τ], the UAV moves in a horizontal plane H with the angular velocity u f i and hence goes over the boundary of the discD + i ( u f i ) defined in Section 6 when preliminarily choosing u f i . By this choice, Hence τ = T i , and by the foregoing, We introduce a Cartesian frame concentric withD Since cos α h > 0 and so N hor = 0 in OZ by (16), and the set D ♦ := [0,  ϕ(t, r) is the polar angle of N hor (t, r) for any (t, r) ∈ D ♦ and ϕ[0, r i (0)] ∈ [0, 2π). By (i) in Assumption 6, As t ranges over [0, T i ], the vector e i (t) rotates counterclockwise with the angular rate u f i . So the continuous function ψ(t) := ϕ[t, r i (t)] runs over an interval whose length does not exceed 2πn i − 2π, whereas the polar angle of e i (t) continuously runs from 0 to 2πn i . So there exist t − , t + ∈ [0, T i ] such that the vector ±e i (t ± ) is aligned with N hor (t ± ) and so e i (t ± ) = ±N hor (t ± ). Hence Here (a) is due to Assumption 2 and Lemma A1. Thus at t = t ± , the continuous map S i (t) takes values of the opposite signs. Hence it vanishes between t − and t + , in violation of (A35). This contradiction proves that not only in the case (A) but also in (B), there exists t sm i for which 1 is true and the UAV is on S Now we "inflate" the targeted isosurface S t ( f ) into the surrounding layer: Lemmas A6, A8, and formula (8) yield the following.
Corollary A1. For any trajectory satisfying 1, 2 in Lemma A8, any UAV i is always inside the operational zone (12) and switches to M at some time T sw i ≥ T tr . In mode M, it undergoes sliding motion over S op i,− and is in the set (A36), whereas f i (t) monotonically converges to f as t → ∞ andḟ i = −χ( f i − f ).
The first statement in this corollary, (iii) in Lemma A5 and (ii) in Lemma A8 justify the pitch angle estimate from the first sentence in Remark 3.

Appendix C. Proofs of the Results from Section: Absence of Collisions among the UAVs
In Appendix C, we performed a purely mathematical study of the closed-loop trajectories. Unlike the real world, a coincidence of the locations of two UAVs at some time entails no implications for such study. Now we will show that this coincidence does not occur in effect. From now on, we focus on the non-chimerical trajectories of the team.
Lemma A9. Two UAVs may collide only if both of them are in mode M and are switched from A to M simultaneously and with the same field value.
Proof. In A, the UAVs are at distinct altitudes by Assumption 5 and (i) in Lemma A5 and so cannot collide. If one of them (say i) is in mode M at time t, then t ≥ T tr , | f i (t) − f | ≤ δ f by (8) and where f i (t) and f j (t) are the field values at the locations of the respective UAVs. Thus these locations are different, and so the UAVs do not collide.
Let both UAVs be in M. By Corollary A1, f i (t), t ≥ T sw i and f j (t), t ≥ T sw j solve a common ODE whose equilibrium f attracts all non-equilibrium solutions at a time-varying nonzero rate. Suppose that T sw i = T sw j , and let T sw where the last inequality excludes being at the same location. The cases where T sw are considered likewise (in the latter case, the arguments start just before the last ⇒ symbol).
Thus it remains to show that two UAVs cannot collide if both are in mode M and launch this mode simultaneously. This will be done by Lemma A13, which is prefaced by several intermediate technical facts.
Lemma A10. At any time t, the horizontal projections of any two points from the layer (A36) are separated by a distance not greater than d hor where (a) holds by Assumption 2. Here r ∈ OZ until F[t, r(s)] reaches f * at some s * , which does happen with Assumption 4 completes the proof. This lemma and Corollary A1 yield the following.
Corollary A2. Whenever two UAVs are both in mode M, their horizontal projections are separated by a distance not greater than d hor Let in M(t) and in A(t) stand for the set of UAVs that are in mode M and A, respectively.
Lemma A11. For any time t and UAV i ∈ in M(t), the following modification of the sets (5) does not alter the output u i in (9) and Ξ(h ± i ) in (10) if the sets (A37) are used and (6) is replaced bẙ Proof. It suffices to show that Ξ(h ± i ) are not altered. We focus on since their effect falls in the saturation zone of Ξ(·) due to (6), (11), (24). Let h i < h j < h i + ∆ h . If (5) considers j, so does (A37). Let j be considered in (A37). If robot j is in mode A, it is considered in (5). Let j be in M. Then the horizontal projections of r i and r j are separated by a distance ≤ d hor + 2δ f b ρ v 2 i ∆ −2 λ by Corollary A2, and the gap between their altitudes ≤ ∆ h . So r i − r j < d vis − by (24), j ∈ E i by (4), and so j is considered in (5). Hence Ξ(h + i ) is not altered.
Suppose that (A43) is untrue, i.e., Σ j|i (τ) ≤ 0 but the set {t ∈ T : Σ j|i (t) > 0} is not empty. Any its connected component is an interval T c whose left end τ − is a root of Σ j|i . However, Σ j|i (t) > 0 and soΣ j|i (t) < 0 on the interior of T c thanks to (A41). Then Σ j|i (τ − ) = 0 ⇒ Σ j|i (t) < 0 ∀t ∈ T c , in violation of the definition of T c . This contradiction completes the proof (A43).
If j ∈ in M(τ), (A44) is established likewise. Otherwise, Σ j (t) = 0 ∀t ≤ T sw j , and so the conclusion in (A44) trivially holds until T sw j . It remains to note that j ∈ in M[T sw j ].
Lemma A13. There are no collisions among the UAVs. If two UAVs start mode M simultaneously, they are always at different altitudes.

Appendix D. Proofs of the Results from Section: Distribution over the Altitudinal Range
In this appendix, the first lemma shows that whenever two UAVs go to different altitudes from a common one, they never return to a common altitude afterward. More precisely, this is true if both UAVs are in mode M and even in a more general situation.
Lemma A14. Suppose that on a time interval T = (τ − , τ + ), always h i (t) = h j (t), one of UAVs i, j is in mode M, whereas the other either is constantly in M or is constantly in A and in the set (4) associated with the other UAV. If h i (τ − ) = h j (τ − ), then h i (τ + ) = h j (τ + ).
This case is considered likewise.
The following lemma shows that for any trajectory of the team, since some time any, two UAVs either are constantly at a common altitude or are constantly at different altitudes.
Lemma A15. For any i = j,there is T such that either (i) h i (t) = h j (t) ∀t ≥ T or (ii) h i (t) = h j (t)∀t ≥ T.
The next lemma sheds light on the evolution of the function Σ m (·) from (A45). Lemma A16. Let on a time interval T, the ith UAV be in mode M, the set in A(t) ∩ E i does not alter, and for any j = i, either constantly h j = h i or constantly h j = h i . Then the function Σ i from (A39) is absolutely continuous on T andΣ i (t)sgn Σ i (t) ≤ −k a.a. t ∈ T, where k is defined by (A40).
This Lemma and Lemma A15 imply the following.
Corollary A3. There is time T such that for t ≥ T, every UAV i is in mode M, undergoes sliding motion with Σ i ≡ 0, and for any Let A(t) be a non-enlargeable group of UAVs from in M(t) that are at a common altitude at time t, and let A k (t), k = 1, . . . , K(t) be all such groups, enumerated in the order of increasing common altitude h c k (t). If for A k (t), there is UAV j ∈ in A(t) at the same altitude and in the set E i (t) of some i ∈ A k (t), the group A k (t) is augmented via adding this j. (By Assumption 5 and (i) in Lemma A5, there is no more than one such j.) A trajectory of the team is said to be altitudinally scattered at time t if the size |A k (t)| ≤ 1 ∀k, and weakly altitudinally scattered on a time interval T if |A k (t)| ≤ 1 ∀k anywhere on T, maybe, except for finitely many points t.
Proof. We observe the state of the team at t = τ and re-enumerate the UAVs with the index i = 0, . . . , N − 1 so that (1) the higher the altitude, the larger the index, (2) for the UAVs at a common altitude, the greater the speedḣ i (τ), the larger the index, and (3) if in the group of the UAVs at the altitude h c k (τ), there are UAVs from in A(τ) \ A k (τ), they are assigned indices lesser than any UAV from A k (τ). If in A k (τ), there are several UAVs with a common speed, including some j ∈ in M(τ) and r ∈ in A(τ+) ∩ E j (τ), then the rth UAV is assigned the least index among them and the jth UAV the next index. (If there are many such jth, an arbitrary j is used.) Let i k 1 be the least index in A k (τ). Then A k (τ) = {i k 1 , i k 1 + 1, . . . , i k n k := i k 1 + n k }, where n k := |A k (τ)|, and i k+1 1 ≥ i k n k + 1. For any k = 1, . . . , K(τ), we define the Lipschitz continuous function h − k of the variables h 0 , . . . , h N−1 ∈ R as the maximum of the following quantities h j (τ) if this max is over a nonempty set, if the first two conditions are not met.
We also define the Lipschitz continuous function h + k as the minimum of the following quantities h j (τ) if this min is over a nonempty set, if the first two conditions are not met. Now we continue the trajectory of the robotic team from its state at t = τ by applying the controller (8)-(10), where for i ∈ A k (τ) ∩ in M(τ) with any k, the second line in (10) is altered into To complete the proof, it suffices to show that the resultant trajectory has the following property: p) The greater the index k of the hosting (at t = τ) group A k (τ) or the index of the UAV within this group, the higher the altitude of the UAV for t > τ, t ≈ τ.
Indeed, for τ + > τ, τ + ≈ τ, (4) and (8) imply that in M(t), in A(t), E i (t) do not alter with t ∈ [τ, τ + ]. Hence any UAV i ∈ in A(τ) is driven by the original control law (8)- (10), and h i (t) = h i (τ) thanks to (i) in Lemma A5. If i ∈ in M(τ), then i ∈ A k (τ) ∩ in M(τ) for some k and so v i is given by (A49). It remains to note that this v i equals , whereh − i andh + i are computed from (A37) and (A38), and to invoke Lemma A11. Now we turn to justify p). Let A r k (τ), r = 1, . . . , R k (τ) be the partition of A k (τ) into groups with a common vertical speedḣ r k (τ) :=ḣ i (τ) ∀i ∈ A r k (τ). For two UAVs from the groups A r k (τ) with different (k, r)'s, the property p) holds due to the continuity argument, the definition of the velocity, (i) in Lemma A5, and Assumption 5. For the case of a common subgroup A r k (τ) with no less than two elements, we examine the behavior of h r k (τ) = 0 and the minimal and maximal i ∈ A r k (τ) ∩ in M(τ) is greater and lesser than i k 1 and i k n k , respectively: For i ∈ A r k (τ) ∩ in M(τ), we have Σ i (τ) = 0 by (A39) and so sliding motion with Σ i (·) ≡ 0 starts at t = τ, which is proven like in (A48). Hence for t > τ, t ≈ τ, the variables h i (t), i ∈ A r k (τ) ∩ in M(τ) obey a system of ODE's that meets the assumptions of Lemma A4 (up to shifts in t and i), where ϑ i := γð(t − T sw i ) due to (11). Then p) holds by Lemma A4. h r k (τ) = 0 but the other assumptions of the previous paragraph are not met: If i = i k n k ∈ A r k (τ) ∩ in M(τ) and T sw i = τ, then Σ i (τ) < 0 by (A39) andḧ i (t) = u h ∀t > τ, t ≈ τ by (A45). Similarly, if i = i k 0 ∈ A r k (τ) ∩ in M(τ) and T sw i = τ, thenḧ i (t) = −u h ∀t > τ, t ≈ τ. As before, we see that sliding motion with Σ i ≡ 0 ∀t > τ, t ≈ τ holds for all other indices i ∈ A r k (τ) ∩ in M(τ) (if they exist). The proof of p) is completed by applying Lemma A4 to them. h r k (τ) = 0: Letḣ r k (τ) > 0, the caseḣ r k (τ) < 0 is handled likewise. Then A r k (τ) ⊂ in M(τ) by (i) in Lemma. A5. For any i ∈ A r k (τ), i = i k n k and t ≈ τ, we have Σ i (t) > 0 by (A39) and soḧ i ≡ −u h by (A45). Since the trajectory is weakly altitudinally scattered on [0, τ], there is at most one index i ∈ A r k (τ) such that Σ i (t) > 0 ∀t ≈ τ. Hence set A r m (τ) consists of i and i + 1 = i k n k , insofar as |A r m (τ| ≥ 2, and Σ i+1 (τ) ≤ 0. If Σ i+1 (τ) = 0, we infer, like before, that sliding motion starts at t = τ in the system (A45) with m = i + 1. This implies thatḧ i+1 > −u h =ḧ i ∀t > τ, t ≈ τ and so p) does hold. If Σ i+1 (τ) < 0, thenḧ i+1 = u h > −u h =ḧ i ∀t > τ, t ≈ τ, and p) holds again.
The times when the set E i (t) is altered for some i do not accumulate since they are separated by periods of no less (d vis − d vis − )/v i due to (4). Since any trajectory is altitudunally scattered on [0, T tr ] thanks to (8), (i) in Lemma A5, and Assumption 5, Lemmas A14 and A17 imply the following.

Corollary A4.
There is a trajectory that starts with the given initial state and has the following property: p) The trajectory is defined on R + := [0, ∞), is weakly altitudinally scattered on R + , and meets 1 and 2 from Lemma A8 for any i. Lemma A18. Any trajectory for which p) from Corollary A4 holds is corporeal.
Proof. It suffices to show that the trajectory is corporeal on any interval T = [τ − , τ + ] that does not contain 0 and times when two UAVs are at a common altitude, or the set E i (t) alters or A → M for some i. Due to (A24) and 1, 2 in Lemma A8, there is only one trajectory of the closed-loop system that starts at t = τ − with the same state as the considered trajectory. Then by (9) and (10), and Corollary 1 in Section 8 from Chapter 2 in [40], the trajectory continuously depends on small perturbations of the team's state at t = τ − , which completes the proof. Lemma A19. Any trajectory for which (ii) from Lemma A15 holds with some i = j is chimerical, and the common altitude of UAVs i and j is constant for t ≥ T, where T is taken from Lemma A15.
Proof. Let T be the time from Corollary A3. By using (22) and increasing T if necessary, we ensure that ð(t − T sw i ) ≥ 1/2 ∀t ≥ T, i. There is a group G of UAVs such that |G| ≥ 2 and By Lemma A13, UAVs i = i from G do not start mode M simultaneously: T sw i = T sw i . For t ≥ T, we see thatḣ i (t) =ḣ i (t) andh ± i (t) =h ± i (t) by Lemma A11, whereas sliding motion with Σ r ≡ 0 occurs for any r since T is borrowed from Corollary A3. So (A39) implies that ð(t − (22). Hence Ξ(h + i ) − Ξ(h − i ) = 0 andḣ i (t) = 0 ∀i ∈ G, t ≥ T. Now we focus on an interval T = [τ − , τ + ], where τ + > τ − (≥ T) will be specified later on. It suffices to show that the trajectory is fully chimerical on T. Let this fail to be true. Almost all perturbations (no matter how small) of the team's state at time τ − bring the UAVs to different altitudes. So by the definition of the full chimericality, there is a sequence of states {x k } at t = τ − with the "different altitudes" property such that any trajectory emitted at t = τ − from x k converges to the original trajectory uniformly on T as k → ∞. For the emitted trajectory, the UAVs remain at different altitudes for t > τ − , t ≈ τ − . By Corollary A4, this trajectory can be extended on [τ − , ∞) to be weakly altitudinally scattered on [τ − , ∞). We shall consider this trajectory, assuming that k is large enough and dropping k in the notations whenever this does not confuse.
We consider a trajectory for which (i) in Lemma A15 is valid whenever i = j. Then there is T such that all UAVs are at different altitudes for t ≥ T. We re-enumerate them in the order of increasing altitude and put d −1 (t) : . We recall that ω-limit point of the bounded function D(·) := [d −1 (·), . . . , d N−1 (·)] is the limit lim k→∞ D(t k ) associated with any sequence {t k } ∞ k=1 such that t k → ∞ as k → ∞. The set Ω of all such points is nonempty and compact. Lemma A20. If (i) from Lemma A15 holds whenever i = j, the following inequality is true: Proof. Suppose to the contrary that the left-hand side of (A53) is less than σ h . We consider D furnishing min Ω in (A53) and the sequence {t k } ∞ k=1 associated with the ω-limit point D, and denote by J the set of all minimizers i for d i with this D. The set {−1, . . . , N − 1} \ J is not empty since So there is i ∈ J such that either i ≥ 0 and d i < d i−1 or i < N − 1 and d i < d i+1 . We focus on the first case; the second one is handled likewise. Due to (i) and (ii) in Lemma A5, |ḋ j (t)| ≤ 2γ∆ h and so |d j (t ) − d j (t )| ≤ 2γ∆ h |t − t | ∀j = −1, . . . , N − 1, t, t , t ≥ 0. (A55) We note that ∆ h > σ h by (23) and (A53). This permits us to pick ξ > 0 and then ε > 0 so that ξ < min{d i−1 ; ∆ h } − d i , ε < min{d i−1 ; ∆ h } − d i − ξ, and ε < ξ/6. We put δ := ε/(4γ∆ h ), t − k := t k − δ, note that d j (t k ) → d j ∀j as k → ∞, and consider so large k's that t − k > T, where T is taken from Due to (11) t) from Lemma A11 and t ∈ [t − k , t k ], we haveh + j (t) = d j (t),h − j (t) = d j−1 (t) and so