Robust Hierarchical Formation Control of Unmanned Aerial Vehicles via Neural-Based Observers

Herein, we investigate the robust formation control problem for a group of unmanned aerial vehicles (UAVs) with system uncertainty. A hierarchical formation control strategy is introduced to ensure the uniform ultimate boundedness of each UAV’s reference tracking error. First, a group of saturated high-level virtual agents are defined to act as the trajectory planners that offer feasible position references to the actual UAVs. A sliding mode neural-based observer is then constructed to estimate the nonlinear uncertainty in the UAV model. Furthermore, sliding mode controllers are designed for both the position loop and the attitude loop of the UAV. To attenuate the chattering phenomenon in the control input, a saturated and smoothed differentiator is proposed along with an observation introduction function. The effectiveness of the proposed control scheme is validated by both the Lyapunov stability theory and numerical simulations based on a multiple-UAV system.

In [16], the time-varying formation tracking problem for a multiple-UAV system subjected to switching communication topologies is investigated by solving an algebraic Riccati equation. Similar results were found for the time-varying formation-containment control problem in [21]. Although the results presented in [16,21] can illustrate the uniform ultimate boundedness of formation tracking errors, the corresponding designs were based on the ideal homogeneous second-order dynamics, which are only applicable in the highlevel path planning perspective, not low-level motion controller design.
Based on the nonlinear cascade model presented in [22], a finite-time leader-follower formation tracking scheme was proposed in [23] for a group of UAVs with undirected topologies. Furthermore, a heterogeneous multi-robot system that includes both UAVs and UGVs was analysed in [24], and a neural-based fault tolerant formation controller was introduced to perform adaptive formation tracking. However, neither of the above two works included the path planning module for UAVs, leading to rapid and excessive changes in the angular status, which is also referred to as aggressive flight [25]. Hence, it is necessary to develop a path planning scheme that offers smooth and suitable reference states to the formation controller to reduce the frequency of aggressive motions.
Model uncertainty and external disturbance are also essential factors to consider for practical applications [26][27][28]. Regarding nonlinear systems with model switching phenomena because of the linearisation operation, the methods that are based on the linear parameter varying modelling theory are helpful for maintaining the robustness of the controller design [29][30][31]. Alternatively, building up nonlinear adaptive observers is necessary to compensate for the system uncertainty if we choose to analyse a system without linearisation. An adaptive terminal sliding mode observer was introduced in [32] to estimate the system uncertainties within finite time. Although the observer given in [32] was found to be effective, this was based on the assumption of correlation with the uncertainty's Lipschitz constant.
Compared to the conventional adaptive observers [33], the neural-based observer presented in [34] has less constraints for the system uncertainty and is more adaptive. The structure in [34] was further modified in [35] to estimate the unknown velocity and the system uncertainty simultaneously for second-order systems. However, the designs in [34] and [35] both use the proportional error reduction term, which makes it hard to guarantee fast convergence of observer states. Therefore, whether we can integrate the neural-based design in [34] with the sliding mode technique in the conventional observer design is a topic worthy of discussion.
Due to the special cascading relationship between a UAV's position loop and attitude loop, the motion controller of a UAV is usually designed in a two-loop structure [15] that requires using the virtual control input in the position loop to get the desired angular status through flight control algorithms. With a strong coupling relationship between the two loops, the control input is sensitive to the sudden changes in system states that can potentially lead to oscillation and chattering [32]. The work in [35] also pointed out that the neural network (NN) output is filled with chattering if the neural weight adjusting law is equipped with static parameters. Hence, it is necessary to find a method to attenuate the potential chattering phenomenon in the control input.
Based on our previous discussions, the contributions of this paper are listed as follows: 1.
Inspired by the work in [36], a hierarchical formation control scheme is employed in this paper. To make the reference provided by the high-level systems feasible for low-level UAVs, a saturated high-level formation controller is proposed to reduce the frequency of aggressive motion.

2.
To improve the error converging speed of the observer structure in [34], the sliding mode technique is integrated with an artificial NN to construct an adaptive observer to estimate the unknown nonlinearities in UAV dynamics, and a fully error-related update law is employed.

3.
To attenuate the chattering phenomenon in the control input [32], a saturated and smoothed differentiator is proposed along with an observation introduction function to reduce the oscillations caused by the differentiating process and the neural-based observer.
The article is structured as follows. The modelling of the multi-UAV system, the basics of graph theory and some background on artificial NN estimation are given in Section 2. The hierarchical formation control scheme and the sliding mode neural-based observer are presented in Section 3. The set-up, results and discussion regarding the numerical simulations of a group of UAVs are illustrated in Section 4. The conclusions are drawn in Section 5.
Notation: In this paper, the term ⊗ represents Kronecker production, I n denotes an identity matrix with the dimension of n. For a matrix M, M F represents its Frobenius norm. If M is a square matrix, we have σ(M) and σ(M) as the maximum and minimal eigenvalues of M, respectively.

Dynamic Model of Uavs
Consider a distributed heterogeneous multi-drone system consists of N(N > 1) UAVs, where the dynamics of the ith UAV is expressed as: [24] where p xi , p yi and p zi represent the global coordinates of the ith UAV; φ i , θ i and ψ i denote the roll angle, pitch angle and yaw angle, respectively; u i,j (j ∈ [1,4]) represents the combined thrust or force provided by the jth motor; R i is the distance between the centre of the drone and the centre of the rotor; and the rest of the parameters are defined in Table 1.

Term(s) Definition
Moments of inertia around the axis m i The mass of the UAV g The gravity constant υ i The drag force coefficient For the sake of simplicity, we make the following definitions to divide the control input of the system into T i , τ i,1 , τ i,2 and τ i,3 : Due to the strong coupling between channels, the position channel and the attitude channel should not be combined into an overall second-order model [23]. Instead, design and analysis based on individual loops are required. Define x i,p = [p xi , p yi , p zi ] T and v i,p = [ṗ xi ,ṗ yi ,ṗ zi ] T . We then have the dynamics of the position loop as for which we have the following equations: If we have w i,p = f i,p +w i,p as the overall system uncertainties in the position loop, Equation (3) can be simplified to the following version: Similarly, if we define x i,a = [φ i , θ i , ψ i ] T and v i,a = [φ i ,θ i ,ψ i ] T , then the dynamics in the attitude loop has the following expression: where For the sake of convenience while analysing the cluster formation tracking behaviour in later parts, it is still necessary to have a unified cluster dynamics expression. For the ith UAV, define T ; then, we have the following simplified version: where we have Definition 1 ([37]). Consider a vector X. We have a correlated continuous Lyapunov function V(X). Then the vector X is said to be semi-globally uniformly ultimately bounded (SGUUB) if V(X) satisfies V(X) = 0 only when X = 0, and there exists a positive boundary b X and a time t X (X(t 0 ), b X ) such that V(X) ≤ b X for all t ≥ t X and X(t 0 ) ∈ Ω V X , where t 0 is the initial time, X(t 0 ) is the initial value of X and Ω V X is a compact set of X.
There are two sets of reference for each UAV to follow, the position reference and the attitude reference. Throughout this paper, we use x d i,p ∈ R 3 and ψ d i ∈ R 1 to represent the position reference and the desired yaw angle for the ith UAV, respectively. The goal of this article is to achieve semi-global uniform ultimate boundedness for each UAV's position tracking error and attitude tracking error, which is illustrated as where both µ p and µ a are small positive constants, and Ω x is a compact set of x i . The following assumptions are made for the system of Equation (7): The trajectory reference x d i,p and its derivativesẋ d i,p andẍ d i,p are all bounded and accessible to the ith UAV. The yaw angle reference ψ d i is bounded and known to the ith UAV.

Graph Theory
In this paper, the distributed communication topology of the multi-UAV system is illustrated by a directed graph G = {R, E, A}, where R = {r 1 , r 2 , . . . , r N } represents the set of nodes, E ⊆ R × R stands for the set of edges and A = [a ij ] ∈ R N×N is the adjacency matrix with nonnegative entries. An edge of the graph G is illustrated as e ij = (r i , r j ), which stands for edge points from node r i to node r j . Self-loops are not considered in this paper, and a ji = 1 if and only if e ij ∈ E. We define deg in (r i ) = ∑ N j=1 a ij to be the in-degree of the ith node, and the degree matrix of the graph is illustrated as D = diag{deg in (r 1 ), deg in (r 2 ), . . . , deg in (r N )}. The Laplacian matrix of the graph is defined as L = D − A. If there always exists a directed path between a pair of distinguished nodes, then the directed graph G is said to be strongly connected. The following lemma is useful for the stability analysis of the proposed control scheme.

Lemma 1 ([38]
). Let the graph G be strongly connected and B ∈ R n×n be a non-negative diagonal matrix with at least one positive element. The matrix (L + B) is considered as an irreducible nonsingular M-matrix. Define We then obtain that P = diag{p i } = diag{1/q i } is a positive definite matrix. Then the matrix Q defined as Q = P(L + B) + (L + B) T P is symmetric and positive definite.

Artificial Nn Approximation
In this study two-layer artificial NNs were employed to estimate the overall system uncertainty w i . According to the universal approximation theorem, a single-layer artificial NN can be used to approximate the unknown function whose variables are restricted to a compact set [39]. Hence, the uncertainty w i can be given in the following form: where ϕ(·) is the activation function vector of the NN, W i is the optimal weight matrix, i is the bounded network approximation bias that satisfies i ≤ M , M is a bounded positive number and Ω v is a compact set of v i . To reduce the complexity of the NN design, we chose the activation function as the identity function, which further leads to and i ∈ R 6 . For the ith UAV, the NN estimation of w i is given as: where W i denotes the estimated weight matrix.
Then the estimation error w i is given as The following assumption is made to ensure the boundedness of the NN output: Remark 1. In a practical task, the external disturbance may include functions that do not use the system states x i and v i as variables. For example, the external disturbancew i can be a function that only relates to the task time t (such asw i = sin(t)). However, since the system states x i and v i are correlated with external variables such as t, and therefore can be expressed as a function whose variables include the external variables, the external variables can also be seen as a function that uses the system states x i and v i as its variables. Hence, we are still able to employ the NN to perform a unified estimation ofw i , which validates the implementation of Equation (9).

Main Results
In this paper, the robust formation control problem of multi-UAV systems is considered, along with the collision avoidance issue when each UAV contains model uncertainty. To reduce the complexity of the controller design, we propose to use a hierarchical two-level formation controller design to separate the concerns.
In high-level designs, virtual agents are generated according to the second-order nominal model of UAVs to act as the reference generators that provide feasible commends to the low-level design, and the low-level controllers are responsible for the actual motion control of physical UAVs. A new neural-based observer design is also proposed in this section for the estimation and compensation of the nonlinear model uncertainties.

High-Level Formation Controller Design
To ensure that the formation controller of each UAV is offered with a sufficient number of state references, each high-level virtual agent is defined to have the homogeneous second-order dynamics as follows: wherex i ∈ R n andv i ∈ R n are the position and velocity information of the virtual agent, respectively;ū i ∈ R n is the control input of the virtual agent; and S(ū i ,Ū Mi ) ∈ R n is the actuator saturation phenomenon. Define S(ū i (j),Ū Mi ) to be the jth element of S(ū i ,Ū Mi ). Then we have whereŪ Mi is a positive constant that represents the saturation limitation. Accordingly, we have the following cluster expression: T Regarding the virtual system Equation (10), we define the high-level tracking errors The tracking error dynamics for the ith virtual agent is given as: Then we have the virtual local formation tracking errorsē xi andē vi as follows (respectively): where b i is the ith diagonal element of B. Defineλ i to be a positive constant. Then the virtual sliding surface is designed as whereτ e is a positive constant,ψ e is a very small positive constant andS (ē xi ,τ e ,ψ e ) ∈ R n is a bounded smooth projection function whose jth element is expressed as: Then we have the time derivative of the virtual sliding surface as:s where the jth element inS d (ē xi ,τ e ,ψ e ) has the following expression: Based on the discussions about the tracking error (Equation (13)) and the sliding surface design (Equation (15)) of the virtual system (Equation (10)), we have the nominal high-level controller design as follows: wherec i andk i are both positive constants.
To ensure that the amplitudes of the control input stay within the saturation limitation, we have the following saturated high-level formation controller: whereτ u andψ u are both positive constants. Now we are ready to present our result within high-level controller design: (10)-where Assumption 1 holds due to the sliding surface design (Equation (15)) and the sliding mode controller (Equation (18)); the variablesS,ē x andδ x are all uniformly ultimately bounded (UUB).

Theorem 1. Consider the virtual cluster equation-Equation
Proof. Consider a Lyapunov candidate as follows: The time derivative of V 1 is given as: First, we rule out the saturation phenomenon (Equation (11)) and haveū i =ū nom i instead to test if the nominal controller is able to ensure the uniform ultimate boundedness of boths i andē xi . Then we have the modified version ofV 1 as: By Lemma 1 and the inequality thatS (ē x ,τ e ,ψ e ) ≤ē x , we have the following norm form: Hence, we get thatV 1 will remain negative until S = S (ē x ,τ e ,ψ e ) = 0 is achieved. By the characteristics of the smooth projection functionS (ē x ,τ e ,ψ e ), ē x will also converge to the value of 0. According to Equation (14), we have that δ x = 0 is ultimately achieved as well. By the definition of uniformly ultimately bounded (UUB) in [38], we have that S , ē x and δ x are all UUB, which indicates that the design ofū nom i is able to achieve convergence of the virtual tracking error.
If we recall the saturation phenomenon (Equation (11)) and the saturated controller design (Equation (18)), similar results are also expected and the states S , ē x and δ x are all UUB, which completes the proof. (10) is constructed with a saturation phenomenon Equation (11) to ensure that the statesx i ,v i andū i are a set of suitable and feasible reference vectors to prevent the low-level UAVs from causing aggressive motion (creating large pitch angles or rotational angles [25]).

Neural-Based Observer Design
Differently from the high-level design, it is vital to consider the system uncertainties w i for the low-level design. To maintain the robustness of the formation tracking process, one popular way is to employ the neural-based observer designs [34,35] to estimate the unknown terms and then perform compensation in the controller design. On the basis of the work in [35], the sliding mode technique is integrated with an artificial NN to approximate the unknown factor w i in system Equation (6).
Although the structures in [34,35] are effective for an arbitrary kind of uncertainty, the design of only using the NN output as the estimation value relies too much on the accuracy of the NN. In other words, if we were to use the designs in [34,35], then the difference between w i and w i would be no less than i , which illustrates their limitation. To overcome this weakness, we propose to have an alternative way of analysing the problem. First, we build up an imaginary second-order observation system according to the actual system (Equation (6)): where u i is the imaginary control input, and vectors x i and v i represent our estimations of states x i and v i , respectively. As u i is the to be designed controller and g i is known in advance, the term g i u i is treated as the known dynamics for the imaginary system. By comparing the difference between the imaginary system Equation (19) and the actual UAV dynamics Equation (6), we have the following tracking error dynamics for the imaginary system: where In theory, we have w i = u i when both x i = 0 and v i = 0 are satisfied. Hence, our goal of building up an adaptive observer to estimate w i is also equivalent to designing a tracking controller u i that reduces the value of x i and v i as much as possible. To achieve the uniformly ultimate boundedness of x i and v i , we define the observation sliding surface as where λ i is a positive constant.
We have the derivative of the observation sliding surface aṡ Based on our previous discussion about the artificial NN estimation of w i Equation (9), we have the following neural adaptive sliding mode controller design for the imaginary system: where c i and k i are both positive constants, and the update law of the NN iṡ where η 1 and η 2 are both positive constants. Now we are ready to present our result of the neural-based sliding mode observer design. (19), where Assumption 2 is satisfied by the observation sliding variable (Equation (21)), the NN estimation (Equation (21)), the adaptive neural weight tuning law (Equation (23)) and the imaginary control input (Equation (22)). We have that the error states s i , x i and W i are all SGUUB if the compact set conditions of the NN hold such that we have x i (t) ∈ Ω x and u i (t) ∈ Ω u when t ≥ t 0 for the ith UAV.

Theorem 2. Consider the system of Equation
Proof. Consider the following Lyapunov candidate: Then we have the derivative of V 2 as follows: We can further modify Equation (24) into the following norm form: Hence,V 2 is said to be negative when the following condition is met: By Definition 1, we have that the vector χ 1 is SGUUB within the following neighbourhood: Hence, the error states s i and x i are both SGUUB. According to the Lyapunov stability theory extension presented in [40], the correlated state W i is also SGUUB, which completes the proof.
As a result, we are confident to have w i ≤ w M , where w M is a positive constant, to support our result in the low-level formation controller design.

Low-Level Formation Controller Design
Regarding the low-level design, we need to first focus on the position loop to provide an essential reference for the attitude control loop [23]. Accordingly, the position loop dynamics (Equation (4)) and the attitude loop dynamics (Equation (5)) of the ith UAV are investigated for the low-level designs.
By Theorem 1, we have that the statesx i ,v i andū i will converge to x d i,p ,ẋ d i,p andẍ d i,p , respectively. Hence, it is reasonable to use statesx i ,v i andū i to act as the references for the ith low-level system. Define δ xi,p and δ vi,p to be the low-level reference tracking errors as follows: Then we have the tracking error dynamics as: Differently from the nominal system designs [35,41], we did not design u i directly from the perspective of Equation (26). Instead, we first defined u i,p = g i,p u i to denote the nominal control input of the position loop. Then we could rewrite Equation (26) as: The position loop sliding surface is constructed as follows: where λ p i is a positive constant. The time derivative of the sliding surface is given as: By the neural-based observer design, we have u i = w i , where w i = [ w T i,p , w T i,a ] T represents our estimation of the system uncertainty w i .
Based on the discussion about the neural-based observer Equation (22) and the potential-based position loop sliding variable (Equation (28)), we have the following nominal controller design: where k p i ∈ R + and c p i ∈ R + . Now we are ready to present our results in the low-level nominal controller design for the position loop Equation (4) of the ith UAV.  (22)) and the nominal control law (Equation (29)). The states δ xi,p and s i,p are both SGUUB if the compact set conditions of the NN hold such that we have x i (t) ∈ Ω x and u i (t) ∈ Ω u when t ≥ t 0 for the ith UAV.
Proof. Consider the following Lyapunov candidate for the ith UAV: The time derivative of the Lyapunov candidate (Equation (30)) is: Alternatively, we have the following matrix form: Hence,V i,p is guaranteed to be negative within the following region: By Definition 1, we have that the vector χ 2 is SGUUB within the following region: As in the case of χ 2 , the values of s i,p and δ xi,p are both SGUUB, which completes the proof.
After obtaining the nominal control u i,p , the next step is to use the flight control techniques to calculate the reference for the row angle and the yaw angle. Inspired by [24], with u i,p = [u i,x , u i,y , u i,z ] T , we have a new guidance law for T d i : After both φ d i and θ d i are calculated, we can use the conventional approach of calculating the slope between the current value and the previous value to get the value ofφ (32) where ξ i,1 is the estimation ofẋ d i,a , t represents the current time and t step is the control step size. Similarly, we have the following structure to get the second-order derivative as: Accordingly, we have ξ i,1 =ẋ d i,a and ξ i,2 =ẍ d i,a when t > 2t step . Similarly to the position loop, define δ xi,a and δ vi,a to be the reference tracking errors in the attitude loop, which further leads to If we define u i,a = g i,a u i to be the nominal control input for the attitude loop, we have the error dynamics as follows: Differently from the position loop, the vector δ vi,a cannot be directly obtained becausė x d i,a is unknown to the ith UAV. Instead, the UAV can only gain access to the estimated error vector δ vi,a as Accordingly, we have the actual sliding surface for the attitude loop as s i,a = δ vi,a + λ a i δ xi,a , and the estimated sliding variable is given as where λ a i is a positive constant. The time derivative of the actual sliding surface s i,a iṡ Based on our discussion about the neural-based observer Equation (22) and the estimated attitude loop sliding variable Equation (36), we have the following nominal controller design: where k a i and c a i are both positive constants. Now we are ready to present our result in the low-level nominal controller design for the attitude loop of the ith UAV:  (22)), the first-order differentiator (Equation (32)), the second-order differentiator (Equation (33) ) and the nominal control law (Equation (37)). The states δ xi,a and s i,a are both SGUUB if the compact set conditions of the NN hold such that we have x i (t) ∈ Ω x and u i (t) ∈ Ω u when t ≥ t 0 for the ith UAV.
Proof. Consider the following Lyapunov candidate.
Then we have the time derivative of V i,a aṡ where w i,p = w i,p − w i,p . We have ξ i,1 =ẋ d i,a , ξ i,2 =ẍ d i,a , δ vi,a = δ vi,a and s i,a = s i,a when t > 2t step . Therefore, Equation (38) should be rewritten as: Similarly, we also have the following matrix form: Hence,V i,a is negative within the following region: By Definition 1, we have that the vector χ 3 is SGUUB within the following region: As in the case of χ 3 , the values of s i,a and δ xi,a are both SGUUB, which completes the proof.
As the position loop and the attitude loop are strongly coupled, a subtle fluctuation in the output of any part of the controller design can lead to oscillation in the system states. Regarding the designs covered by Theorems 1-4, there are two factors that can lead to potential oscillations in the system's control input:

1.
The values of θ d i and ψ d i are not guaranteed to be smooth because they are generated by the reverse calculation (Equation (31)). Hence, the output values of the direct derivative structures (Equations (32) and (33)) can be filled with chattering because of the discontinuity of their input.

2.
Before W i − W i F is settled within a neighbourhood around 0, the output of the NN is usually filled with oscillations [35].
Before we introduce the saturated and smoothed differentiator, we need to first make the following assumption: Consequently, we have the saturated and smoothed differentiator designs, as shown in Algorithm 1, where l ξ is initially set as 0. Accordingly, we have where t diff is the differentiate step size.
To reduce the negative effect brought about by the fluctuations in the NN's output, the following observation introduction function is proposed to smoothly introduce the uncertainty estimations w i,p and w i,a into the low-level controller designs.
where γ 1 and γ 2 are both positive constants. Hence, instead of using w i,p and w i,a directly, we have the following smoothed position controller and attitude controller: Algorithm 1: Saturated and smoothed differentiator D(t delay , t diff , ξ + in , t, ξ M ).
Ultimately, we have the following motion controller design: In all, the hierarchical formation controller is illustrated as Figure 1. Now, we are ready to present our final result in the overall system design.

To the Lower Level
From the Higher Level ∈ Ω x and u i (t) ∈ Ω u when t ≥ t 0 for the ith UAV.
Proof. By Theorem 1, we have Since the value of the introduction function Equation (40) will converge to 1 ultimately, the result and analysis related to Theorem 3 remain the same.
Although the smoothed design in Algorithm 1 will lead to differences between the set {ξ i,1 , ξ i,2 } and {ẋ d i,a ,ẍ d i,a }, the difference should be small and bounded if the value of t diff is chosen properly. Suppose we have ξ i,1 −ẋ d i,a ≤ ξ 1 M and ξ i,2 −ẍ d i,a ≤ ξ 2 M , similar to the proof of Theorem 4. We have that where ξ 1 M and ξ 2 M are both small positive constants. With the actual motion controller chosen as Equation (43), we have that: which matches the goal of this paper (Equation (8)) and concludes the proof.

Simulation Results and Discussion
To justify the performance of the proposed hierarchical formation controller design, comparative simulations based on a multi-UAV system were conducted.
Consider a multi-UAV system that contains six heterogeneous UAVs whose dynamics are expressed as Equation (1). The system's parameter values are given in Table 2 The gravity constant g was set to g = 9.81. The communication topology of the system was chosen as shown in Figure 2, and we used b i = 1 for i ∈ [1, 6]. The position reference for the ith UAV was chosen as The reference of the yaw angle was chosen as ψ d i = 0. The system uncertainties were chosen as follows: The parameters of the high level controller were chosen asŪ Mi = 0.2,τ e = 0.5, ψ e = 0.05,c i = 1,λ i = 2,k i = 1.5,τ u = 0.19 andψ u = 0.01. To illustrate the effectiveness of the high-level formation controller (Equation (18)), the norms of the high-level error-related vectorsδ xi ands i are given in Figures 3 and 4, respectively.  Although the performance of the sliding mode controller (Equation (18)) was affected by the saturation phenomenon (Equation (11)) and obvious state overshoots can be noted, the uniform ultimate boundedness of bothδ xi ands i was still achieved simultaneously, which validates the results presented in Theorem 1.
The parameters of the neural-based observer were selected as λ i = 2, c i = 4, k i = 4, η 1 = 20 and η 2 = 0.5. Define x , s and d as follows to represent the overall state observation error norm, overall observation sliding variable norm and uncertainty estimation error norm, respectively.  Figure 5, where we found that all three norms are SGUUB for 2 × 10 −3 . Hence, Theorem 2 was verified. However, both overshooting and chattering were observed in the output of the neural-based observer within the first 2 s, which validates our motivation of employing the observation introduction function Equation (40). The parameters of the low-level position controller were chosen as λ p i = 2, c p i = 2 and k p i = 2; and the parameters of the low-level attitude controller were selected as λ a i = 2, c a i = 4 and k a i = 4. To illustrate the stability of our low-level controller design, we define error norm vectors δ x,a , δ x,p , s a and s p as follows: The trends of δ x,a , δ x,p , s a and s p are recorded in Figure 6, where they are SGUUB for 4 × 10 −4 , 7 × 10 −4 , 3 × 10 −4 and 1.5 × 10 −3 . Hence, the both Theorem 3 and 4 are valid. To prove that the observation introduction function (Equation (40)) and the saturated differentiator (Algorithm 1) help attenuate the chattering in the control input, the following two comparative simulations were conducted: 1.
The original controller (TOC) where both Equations (32) and (33) are employed to provide the derivatives of the attitude reference x d i,a . The derivatives are saturated with ξ i,1 = S(ξ i,1 , 0.2) and ξ i,2 = S(ξ i,2 , 0.1).

2.
The smoothed controller (TSC) where the observation introduction function Equation (40) is employed with γ 1 = 4 and γ 2 = 2, and the attitude reference derivatives are obtained as   As the introduction function (Equation (40)) is not used in the TOC design, more chattering was observed in the control input for the first 4 s. Besides, the TSC design was found to have smoother control input ultimately (see UAV 2 and 6) due to the implementation of the smoothed differentiator (Algorithm 1). Therefore, the effectiveness of the TSC design and its superiority over the TOC design are both illustrated.
Finally, define x d p = [x d 1,p , x d 2,p , . . . , x d N,p ] T and x p = [x 1,p , x 2,p , . . . , x N,p ] T to represent the overall formation reference and the actual position of the multi-UAV system, respectively. Then we have the norm of the overall reference tracking error of the system with TSC design, as shown in Figure 9, where the semi-global uniform ultimate boundedness of x p − x d p proves the validity of Theorem 5. The trajectories of all UAVs are given in Figure 10 to illustrate the movement and formation status of the entire system. According to Equation (46), the desired formation is a rotating circle that keeps on elevating (see the dotted grey circle). In Figure 10, each UAV successfully reached its position reference within bounded error to form the desired formation, which illustrates the effectiveness of the proposed design structure in Figure 1. The angular status of each UAV is also given in Figure 11 to indicate that the saturation phenomenon employed in the higher level helps restrict the pitch and roll angle within [−30 • , 30 • ] to further reduce the aggressive motion of each UAV.

Conclusions
In this paper, the robust formation control problem for a group of UAVs with system uncertainties was investigated. A neural-based hierarchical formation controller was proposed to ensure the semi-global uniform ultimate boundedness of the system's formation tracking error. A cluster of saturated high-level agents were defined to act as the smooth reference generator, and a saturated sliding mode controller was proposed to ensure the convergence of the high-level tracking errors. The sliding mode technique was further integrated with the neural-based observer structure to estimate system uncertainties in both the position loop and the angular loop. Adaptive sliding mode observers were then introduced for both position and angular loops to perform adaptive estimation of the system uncertainties. To reduce the chattering and oscillation in the control input, a combination of the smoothed differentiator and the observation introduction function is employed for each UAV. The effectiveness of the robust and adaptive hierarchical formation control scheme was validated by both theoretical analysis and comparative simulations.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.

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