Distributed Robust Formation Tracking Control for Quadrotor UAVs with Unknown Parameters and Uncertain Disturbances

1. Introduction

In recent years, cooperative control of quadrotor unmanned aerial vehicles (UAVs) has received considerable attention due to its broad applications in areas such as wireless communication, nuclear radiation detection, and agricultural mapping. Formation control is one of the most important research areas in the field of cooperative control for quadrotor UAVs. A formation composed of multiple low-cost UAVs can replace an expensive multi-functional UAV to complete complex tasks. In addition, UAV formations can provide system redundancy and reconfiguration ability [1]. Formation control of quadrotor UAVs has garnered significant research attention owing to its prospective applications in both military and civilian domains [2,3,4]. From the perspective of control mechanisms, the existing methods of quadrotor UAV formation control include leader-follower method [5], artificial potential method [6], behavior-based method [7], etc. Recently, the work in [8] studied dynamic formation collision avoidance control for quadrotor UAVs using the virtual structure method. In [9], a consensus-based method was used to design a time-varying formation tracking control scheme for quadrotor UAVs. However, the quadrotor UAV model considered in the above literature is simplified, and the designed formation control schemes rely on the model parameters of the quadrotor UAV. In many practical applications of quadrotor UAVs, it can be difficult to accurately obtain model parameters. As a result, it is important to design a formation control scheme for quadrotor UAVs that does not rely on the use of model parameters. In addition, quadrotor UAVs are very sensitive to uncertain disturbances, and it is necessary to design effective disturbance rejection flight control schemes. The disturbance rejection control of individual UAVs has been extensively studied in the existing literature [10,11,12,13,14,15]. For quadrotor UAV swarms, uncertain disturbances acting on each UAV will affect neighboring UAVs through the communication network. Hence, designing disturbance rejection control schemes for quadrotor UAV swarms is a more complex task. The existing literature has not extensively investigated the problem of disturbance rejection control for quadrotor UAV swarms, which underscores the significance of the research presented in this paper.

In this paper, a distributed robust formation tracking control method is proposed for quadrotor UAVs with unknown parameters and uncertain disturbances. The proposed method has the following novelties. First, a more practical formation tracking control method is proposed in this paper, which does not need to use the model parameters of the quadrotor UAV. Second, an adaptive disturbance rejection control scheme for quadrotor UAV swarms is developed. In the presence of uncertain disturbances, this scheme can still achieve formation tracking control for quadrotor UAV swarms, and the tracking error can eventually converge to zero.

The structure of this paper is arranged in the following manner. In Section II and Section III, a distributed formation tracking control scheme and an adaptive disturbance rejection attitude control method are designed for quadrotor UAVs. The efficacy of the proposed control method is validated in Section IV. Finally, Section V concludes the paper.

2. Distributed Robust Formation Tracking Control for Quadrotor UAVs

In this section, a distributed formation flight control method is developed for quadrotor UAVs to achieve the following three control objectives: 1) form the desired formation; 2) track the desired flight trajectory; 3) reduce the influence of uncertain disturbances.

2.1. Graph Theory

The communication topology among a group of N quadrotor UAVs is considered as an undirected graph $\mathrm{G}=(\mathbb{W},\mathbb{S})$, where $\mathbb{W}\triangleq \{1,2,\cdots ,N\}$ denotes the vertex set and $\mathbb{S}\triangleq \{(i,j):i\in \mathbb{W},j\in {\mathcal{N}}_i\}$ denotes the edge set. The neighbor set of the ith UAV is ${\mathcal{N}}_i\triangleq \{j\in \mathbb{W}$: there is a communication link between UAV i and UAV j, $j\ne i\}$. Define a weight ${\mathbf{a}}_{ij}$ for each edge $(i,j)\in \mathbb{S}$, ${\mathbf{a}}_{ij}=1$ if $j\in {\mathcal{N}}_i$, and ${\mathbf{a}}_{ij}=0$ otherwise. The Laplacian matrix is $\mathcal{L}=\left[{\mathbf{w}}_{ij}\right]\in {\mathbb{R}}^{N\times N}$, where ${\mathbf{w}}_{ii}={\sum }_{j=1,j\ne i}^N{\mathbf{a}}_{ij}$ and ${\mathbf{w}}_{ij}=-{\mathbf{a}}_{ij}$$(j\ne i)$. The leader adjacency matrix is $\mathcal{D}=diag\{{\mathbf{d}}_1,\cdots ,{\mathbf{d}}_N\}$, where ${\mathbf{d}}_i>0$ if UAV i can obtain the desired flight trajectory and ${\mathbf{d}}_i=0$ otherwise. An undirected graph is considered connected if there is a path between every pair of distinct vertices.

Next, two useful lemmas are introduced.

 Lemma 1.

[16]. If the undirected graph $\mathrm{G}$ is connected, and at least one UAV can obtain the desired flight trajectory, then the symmetric matrix $\mathcal{L}+\mathcal{D}$ is positive definite.

 Lemma 2.

[17]. For any $\kappa >0$ and $\varsigma \in \mathbb{R}$, the inequality $0\le |\varsigma |-\frac{{\varsigma }^2}{\sqrt{{\varsigma }^2+{\kappa }^2}}\le \kappa$ holds.

2.2. Quadrotor UAV Position Dynamic Model

In this paper, define $\mathcal{E}={[\varphi ,\theta ,\phi ]}^T$ as the attitude of the quadrotor UAV, where $\varphi ,$$\theta$ and $\phi$ denote the angles of roll, pitch and yaw, respectively. The rotation matrix that describes the transformation from the body-fixed frame to the earth-fixed frame is denoted as

$$R_t=\left[\begin{array}{ccc}{C}_{\phi }{C}_{\theta }& C_{\phi }{S}_{\varphi }{S}_{\theta }-S_{\phi }{C}_{\varphi }& C_{\phi }{C}_{\varphi }{S}_{\theta }+S_{\phi }{S}_{\varphi }\\ S_{\phi }{C}_{\theta }& S_{\phi }{S}_{\varphi }{S}_{\theta }+C_{\phi }{C}_{\varphi }& S_{\phi }{C}_{\varphi }{S}_{\theta }-C_{\phi }{S}_{\varphi }\\ -S_{\theta }& S_{\varphi }{C}_{\theta }& C_{\varphi }{C}_{\theta }\end{array}\right]$$

(1)

where $S_{(·)}$ and $C_{(·)}$ denote sin$(·)$ and cos$(·)$, respectively.

Define $h={[x,y,z]}^T$ as the position of the quadrotor UAV. As described in [18], the translational dynamic equations are given as

$$\begin{array}{c}\hfill \left[\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{z}\end{array}\right]=R_t\left[\begin{array}{c}0\\ 0\\ U_s\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ -g\end{array}\right]-\frac{1}{m}\left[\begin{array}{c}{d}_x\dot{x}\\ d_y\dot{y}\\ d_z\dot{z}\end{array}\right]\end{array}$$

(2)

where m is the quadrotor mass; $d_x$, $d_y$, $d_z$ are the air drag coefficients; $U_s=b({\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)/m$, b is the lift coefficient and ${\Omega }_k$$(k=1,2,3,4)$ are the rotor speed; g is the acceleration of gravity.

In this paper, the formation tracking control problem of quadrotor UAVs is studied. From (2), the position dynamic system of the ith $(i=1,\cdots ,N)$ UAV can be described as

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{h}}_{i,l}=v_{i,l}+w_{h_{i,l}}\hfill \\ {\dot{v}}_{i,l}=u_{i,l}+{\mathsf{\Theta }}_{i,l}{v}_{i,l}+w_{v_{i,l}},l=1,2,3\hfill \end{array}\right.\end{array}$$

(3)

where $[h_{i,1},h_{i,2},h_{i,3}]=[x_i,y_i,z_i]$ and $[v_{i,1},v_{i,2},v_{i,3}]=[{\dot{x}}_i,{\dot{y}}_i,{\dot{z}}_i]$ are the position and velocity of UAV i, respectively; $[{\mathsf{\Theta }}_{i,1},{\mathsf{\Theta }}_{i,2},{\mathsf{\Theta }}_{i,3}]=[-\frac{d_x}{m},-\frac{d_y}{m},-\frac{d_z}{m}]$ are the unknown system parameters; $u_{i,1},u_{i,2},u_{i,3}$ are the control inputs, and $u_{i,1}=(C_{{\varphi }_i}{S}_{{\theta }_i}{C}_{{\phi }_i}+S_{{\varphi }_i}{S}_{{\phi }_i})U_{s_i},$ $u_{i,2}=(C_{{\varphi }_i}{S}_{{\theta }_i}{S}_{{\phi }_i}-S_{{\varphi }_i}{C}_{{\phi }_i})U_{s_i},$$u_{i,3}=C_{{\varphi }_i}{C}_{{\theta }_i}{U}_{s_i}-g.$ In addition, $w_{h_{i,l}}$ and $w_{v_{i,l}}$ represent uncertain disturbances.

 Assumption 1.

The uncertain disturbances satisfy

$$\begin{array}{c}\hfill |w_{h_{i,l}}|\le {\overline{w}}_h,|w_{v_{i,l}}|\le {\overline{w}}_v,i=1,\cdots ,N\end{array}$$

(4)

where ${\overline{w}}_h>0$ and ${\overline{w}}_v>0$ are unknown constants.

Definition: A time-varying formation formed by a group of N UAVs is specified by $\mathbf{F}(t)={[{\mathbf{F}}_1^T(t),\cdots ,{\mathbf{F}}_N^T(t)]}^T\in {\mathbb{R}}^{3N}$, where ${\mathbf{F}}_i(t)={[{\mathbf{F}}_{i,1},{\mathbf{F}}_{i,2},{\mathbf{F}}_{i,3}]}^T\in {\mathbb{R}}^3$ is the piecewise continuously formation vector. Formation tracking control of quadrotor UAVs can be achieved if

$$\begin{array}{c}\hfill \underset{t\to +\infty }{lim}[h_i(t)-{\mathbf{F}}_i(t)-\mathbf{n}(t)]=0,i=1,\cdots ,N\end{array}$$

(5)

where $h_i(t)={[x_i,y_i,z_i]}^T$; and $\mathbf{n}(t)={[{\mathbf{n}}_1,{\mathbf{n}}_2,{\mathbf{n}}_3]}^T$ represents the desired flight trajectory.

 Assumption 2.

The desired flight trajectory satisfies

$$\begin{array}{c}\hfill |{\mathbf{n}}_l|\le \overline{\mathbf{n}},l=1,2,3\end{array}$$

(6)

where $\overline{\mathbf{n}}>0$ is an unknown constant.

2.3. Distributed Formation Tracking Controller Design

For the ith UAV, define two error variables

$$\begin{array}{c}\hfill \begin{array}{c}{\chi }_{i,l}=\sum _{j=1}^N{\mathbf{a}}_{ij}(h_{i,l}-{\mathbf{F}}_{i,l}-h_{j,l}+{\mathbf{F}}_{j,l})\hfill \\ +{\mathbf{d}}_i(h_{i,l}-{\mathbf{F}}_{i,l}-{\mathbf{n}}_l)\hfill \end{array}\end{array}$$

(7)

$$\begin{array}{c}\hfill {\eta }_{i,l}=v_{i,l}-{\alpha }_{i,l},l=1,2,3\end{array}$$

(8)

where ${\alpha }_{i,l}$ are the virtual control functions. The detailed design procedure is given as follows:

Step 1: By defining ${\chi }_l={[{\chi }_{1,l},\cdots ,{\chi }_{N,l}]}^T$, one can obtain ${\chi }_l=(\mathcal{L}+\mathcal{D})e_l$, where $e_l={\underline{h}}_l-{\underline{\mathbf{F}}}_l-{\underline{\mathbf{n}}}_l$ with ${\underline{h}}_l={[h_{1,l},\cdots ,h_{N,l}]}^T$, ${\underline{\mathbf{F}}}_l={[{\mathbf{F}}_{1,l},\cdots ,{\mathbf{F}}_{N,l}]}^T$, and ${\underline{\mathbf{n}}}_l={[{\mathbf{n}}_l,\cdots ,{\mathbf{n}}_l]}^T$. Then, the derivative of $e_l$ satisfies

$$\begin{array}{c}\hfill {\dot{e}}_l=\left[\begin{array}{c}{\alpha }_{1,l}+{\eta }_{1,l}+w_{h_{1,l}}-{\dot{\mathbf{F}}}_{1,l}-{\dot{\mathbf{n}}}_l\\ ⋮\\ {\alpha }_{N,l}+{\eta }_{N,l}+w_{h_{N,l}}-{\dot{\mathbf{F}}}_{N,l}-{\dot{\mathbf{n}}}_l\end{array}\right].\end{array}$$

(9)

The virtual control function ${\alpha }_{i,l}$ is chosen as

$$\begin{array}{c}\hfill \begin{array}{c}{\alpha }_{i,l}=-k_{h_{i,l}}{\chi }_{i,l}+{\dot{\mathbf{F}}}_{i,l}+{\mathbf{d}}_i{\dot{\mathbf{n}}}_l-\frac{{\chi }_{i,l}}{\sqrt{{\chi }_{i,l}^2+{\delta }_{i,l}^2}}{\widehat{\mu }}_{h_{i,l}}\hfill \end{array}\end{array}$$

(10)

where $k_{h_{i,l}}>0$ is a design constant; ${\widehat{\mu }}_{h_{i,l}}$ is the estimate of ${\mu }_{h_{i,l}}={\overline{w}}_h+(1-{\mathbf{d}}_i)\overline{\mathbf{n}}$; and ${\delta }_{i,l}(t)$ is a positive continuous function satisfying ${lim}_{t\to \infty }{\int }_{t_0}^t{\delta }_{i,l}(\tau )d\tau \le {\overline{\delta }}_{i,l}<+\infty$, and ${\overline{\delta }}_{i,l}$ is a positive constant. Consider the Lyapunov function

$$\begin{array}{c}\hfill \begin{array}{c}{V}_{h_l}=\frac{1}{2}{e}_l^T(\mathcal{L}+\mathcal{D})e_l+\sum _{i=1}^N\frac{1}{2{\lambda }_{h_{i,l}}}{\tilde{\mu }}_{h_{i,l}}^2\hfill \end{array}\end{array}$$

(11)

where the estimation error ${\tilde{\mu }}_{h_{i,l}}={\mu }_{h_{i,l}}-{\widehat{\mu }}_{h_{i,l}}$; and ${\lambda }_{h_{i,l}}>0$ is a design parameter.

From (9)-(11), the derivative of $V_{h_l}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{h_l}\le \sum _{i=1}^N[-k_{h_{i,l}}{\chi }_{i,l}^2+(|{\chi }_{i,l}|-\frac{{\chi }_{i,l}^2}{\sqrt{{\chi }_{i,l}^2+{\delta }_{i,l}^2}}){\mu }_{h_{i,l}}\hfill \\ +{\chi }_{i,l}{\eta }_{i,l}+\frac{1}{{\lambda }_{h_{i,l}}}{\tilde{\mu }}_{h_{i,l}}(\frac{{\lambda }_{h_{i,l}}{\chi }_{i,l}^2}{\sqrt{{\chi }_{i,l}^2+{\delta }_{i,l}^2}}-{\dot{\widehat{\mu }}}_{h_{i,l}})].\hfill \end{array}\end{array}$$

(12)

The parameter update law is chosen as

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\widehat{\mu }}}_{h_{i,l}}=\frac{{\lambda }_{h_{i,l}}{\chi }_{i,l}^2}{\sqrt{{\chi }_{i,l}^2+{\delta }_{i,l}^2}}.\hfill \end{array}\end{array}$$

(13)

Then, by applying Lemma 2, we have

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{h_l}\le \sum _{i=1}^N(-k_{h_{i,l}}{\chi }_{i,l}^2+{\chi }_{i,l}{\eta }_{i,l}+{\delta }_{i,l}{\mu }_{h_{i,l}}).\hfill \end{array}\end{array}$$

(14)

Step 2: From (3) and (10), the derivative of ${\eta }_{i,l}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\eta }}_{i,l}=u_{i,l}+{\mathsf{\Theta }}_{i,l}{v}_{i,l}+w_{v_{i,l}}-\frac{\partial {\alpha }_{i,l}}{\partial h_{i,l}}(v_{i,l}+w_{h_{i,l}})\hfill \\ -\frac{\partial {\alpha }_{i,l}}{\partial {\delta }_{i,l}}{\dot{\delta }}_{i,l}-\frac{\partial {\alpha }_{i,l}}{\partial {\widehat{\mu }}_{h_{i,l}}}{\dot{\widehat{\mu }}}_{h_{i,l}}-{\mathbf{d}}_i\frac{\partial {\alpha }_{i,l}}{\partial {\mathbf{n}}_l}{\dot{\mathbf{n}}}_l\hfill \\ -{\mathbf{d}}_i\frac{\partial {\alpha }_{i,l}}{\partial {\dot{\mathbf{n}}}_l}{\ddot{\mathbf{n}}}_l-\sum _{j=1}^N{\mathbf{a}}_{ij}\frac{\partial {\alpha }_{i,l}}{\partial h_{j,l}}(v_{j,l}+w_{h_{j,l}}).\hfill \end{array}\end{array}$$

(15)

The formation flight controller is designed as

$$\begin{array}{c}\hfill \begin{array}{c}{u}_{i,l}=-k_{v_{i,l}}{\eta }_{i,l}-{\widehat{\mathsf{\Theta }}}_{i,l}{v}_{i,l}+\frac{\partial {\alpha }_{i,l}}{\partial h_{i,l}}{v}_{i,l}+\frac{\partial {\alpha }_{i,l}}{\partial {\delta }_{i,l}}{\dot{\delta }}_{i,l}\hfill \\ +\frac{\partial {\alpha }_{i,l}}{\partial {\widehat{\mu }}_{h_{i,l}}}{\dot{\widehat{\mu }}}_{h_{i,l}}+{\mathbf{d}}_i\frac{\partial {\alpha }_{i,l}}{\partial {\mathbf{n}}_l}{\dot{\mathbf{n}}}_l+{\mathbf{d}}_i\frac{\partial {\alpha }_{i,l}}{\partial {\dot{\mathbf{n}}}_l}{\ddot{\mathbf{n}}}_l\hfill \\ +\sum _{j=1}^N{\mathbf{a}}_{ij}\frac{\partial {\alpha }_{i,l}}{\partial h_{j,l}}{v}_{j,l}-\frac{{\eta }_{i,l}{\varpi }_{i,l}^2}{\sqrt{{\eta }_{i,l}^2{\varpi }_{i,l}^2+{\delta }_{i,l}^2}}{\widehat{\mu }}_{v_{i,l}}-{\chi }_{i,l}\hfill \end{array}\end{array}$$

(16)

where $k_{v_{i,l}}>0$ is a design constant; ${\varpi }_{i,l}=\sqrt{1+{(\frac{\partial {\alpha }_{i,l}}{\partial h_{i,l}})}^2+{\sum }_{j=1}^N{\mathbf{a}}_{ij}{(\frac{\partial {\alpha }_{i,l}}{\partial h_{j,l}})}^2}$; and ${\widehat{\mu }}_{v_{i,l}}$ is the estimate of ${\mu }_{v_{i,l}}=max\{{\overline{w}}_h,{\overline{w}}_v,\overline{\mathbf{n}}\}$.

Construct the following Lyapunov function

$$\begin{array}{c}\hfill \begin{array}{c}{V}_{v_l}=V_{h_l}+\frac{1}{2}\sum _{i=1}^N({\eta }_{i,l}^2+\frac{1}{{\lambda }_{h_{i,l}}}{\tilde{\mu }}_{v_{i,l}}^2+\frac{1}{{\lambda }_{{\mathsf{\Theta }}_{i,l}}}{\tilde{\mathsf{\Theta }}}_{i,l}^2)\hfill \end{array}\end{array}$$

(17)

where the estimation errors ${\tilde{\mu }}_{v_{i,l}}={\mu }_{v_{i,l}}-{\widehat{\mu }}_{v_{i,l}}$ and ${\tilde{\mathsf{\Theta }}}_{i,l}={\mathsf{\Theta }}_{i,l}-{\widehat{\mathsf{\Theta }}}_{i,l}$; ${\lambda }_{v_{i,l}}>0$ and ${\lambda }_{{\mathsf{\Theta }}_{i,l}}>0$ are design parameters.

From (14)-(17), the derivative of $V_{v_l}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{v_l}\le \sum _{i=1}^N[-k_{h_{i,l}}{\chi }_{i,l}^2+(|{\eta }_{i,l}|{\varpi }_{i,l}-\frac{{\eta }_{i,l}^2{\varpi }_{i,l}^2}{\sqrt{{\eta }_{i,l}^2{\varpi }_{i,l}^2+{\delta }_{i,l}^2}}){\mu }_{v_{i,l}}\hfill \\ -k_{v_{i,l}}{\eta }_{i,l}^2+\frac{1}{{\lambda }_{v_{i,l}}}{\tilde{\mu }}_{v_{i,l}}(\frac{{\lambda }_{v_{i,l}}{\eta }_{i,l}^2{\varpi }_{i,l}^2}{\sqrt{{\eta }_{i,l}^2{\varpi }_{i,l}^2+{\delta }_{i,l}^2}}-{\dot{\widehat{\mu }}}_{v_{i,l}})\hfill \\ +\frac{1}{{\lambda }_{{\mathsf{\Theta }}_{i,l}}}{\tilde{\mathsf{\Theta }}}_{i,l}({\lambda }_{{\mathsf{\Theta }}_{i,l}}{v}_{i,l}{\eta }_{i,l}-{\dot{\widehat{\mathsf{\Theta }}}}_{i,l})+{\delta }_{i,l}{\mu }_{h_{i,l}}].\hfill \end{array}\end{array}$$

(18)

The adaptive update laws are chosen as

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\widehat{\mu }}}_{v_{i,l}}=\frac{{\lambda }_{v_{i,l}}{\eta }_{i,l}^2{\varpi }_{i,l}^2}{\sqrt{{\eta }_{i,l}^2{\varpi }_{i,l}^2+{\delta }_{i,l}^2}},{\dot{\widehat{\mathsf{\Theta }}}}_{i,l}={\lambda }_{{\mathsf{\Theta }}_{i,l}}{v}_{i,l}{\eta }_{i,l}.\hfill \end{array}\end{array}$$

(19)

Then, by applying Lemma 2, we have

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{v_l}\le \sum _{i=1}^N(-k_{h_{i,l}}{\chi }_{i,l}^2-k_{v_{i,l}}{\eta }_{i,l}^2+{\delta }_{i,l}{\mu }_{h_{i,l}}+{\delta }_{i,l}{\mu }_{v_{i,l}}).\hfill \end{array}\end{array}$$

(20)

Now, we present the analysis results.

 Theorem 1.

Consider the quadrotor UAV swarm system (3), the formation tracking controller (16), and the adaptive laws (13) and (19). All the signals in the closed-loop system are globally bounded, and the quadrotor UAV swarm can achieve time-varying formation flying and track the virtual leader.

$Proof:$ Integrating both sizes of (20), it follows that

$$\begin{array}{c}\hfill \begin{array}{c}{V}_{v_l}(t)+k_{h_{i,l}}{\int }_0^t{\chi }_{i,l}^2(\tau )d\tau +k_{v_{i,l}}{\int }_0^t{\eta }_{i,l}^2(\tau )d\tau \hfill \\ \le V_{v_l}(0)+({\mu }_{h_{i,l}}+{\mu }_{v_{i,l}}){\overline{\delta }}_{i,l}.\hfill \end{array}\end{array}$$

(21)

From the definition of $V_{v_l}$ in (17), one can get that ${\chi }_{i,l}$, ${\eta }_{i,l}$, ${\widehat{\mu }}_{h_{i,l}}$, ${\widehat{\mu }}_{v_{i,l}}$, and ${\widehat{\mathsf{\Theta }}}_{i,l}$$(l=1,2,3)$ are bounded. From (10), (16), and Lemma 1, ${\alpha }_{i,l}$, $u_{i,l}$, and $h_{i,l}$ are bounded. Therefore, the boundedness of all the signals is guaranteed, and ${\dot{\chi }}_{i,l}$ is bounded. By applying Barbalat’s lemma, one has ${lim}_{t\to \infty }{\chi }_{i,l}(t)=0$. From the definition of ${\chi }_{i,l}$ and Lemma 1, it follows that formation tracking control of quadrotor UAVs can be achieved, i.e. ${lim}_{t\to +\infty }[h_i(t)-{\mathbf{F}}_i(t)-\mathbf{n}(t)]=0$. This completes the proof.

 Remark 1.

When the distributed formation tracking controller $u_{i,1},u_{i,2},u_{i,3}$ is designed, and the desired yaw angle ${\phi }_{i,0}$ is treated as an additional reference signal, then the desired roll angle ${\varphi }_{i,0}$, the desired pitch angle ${\theta }_{i,0}$, and the control input $U_{s_i}$ can be obtained in the following way

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\theta }_{i,0}=arctan(\frac{C_{{\phi }_{i,0}}{u}_{i,1}+S_{{\phi }_{i,0}}{u}_{i,2}}{u_{i,3}+g})\hfill \\ {\varphi }_{i,0}=arctan(\frac{(S_{{\phi }_{i,0}}{u}_{i,1}-C_{{\phi }_{i,0}}{u}_{i,2})C_{{\theta }_{i,0}}}{u_{i,3}+g})\hfill \\ U_{s_i}=\sqrt{u_{i,1}^2+u_{i,2}^2+{(u_{i,3}+g)}^2}.\hfill \end{array}\right.\end{array}$$

(22)

Since $u_{i,1},u_{i,2},u_{i,3}$, and ${\phi }_{i,0}$ are continuous and bounded, it is known that ${\theta }_{i,0},{\varphi }_{i,0}$ and $U_{s_i}$ are bounded.

3. Disturbance Rejection Control of Quadrotor UAV Attitude

In this section, an adaptive disturbance rejection attitude control method will be designed for the quadrotor UAV.

3.1. Quadrotor UAV Attitude Dynamic Model

The angular velocity with respect to the attitude is given as $\mathcal{W}={[p,q,r]}^T$. As described in [19], The correlation between the attitude angle and angular velocity can be denoted by

$$\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\phi }\end{array}\right)=\left[\begin{array}{ccc}1& T_{\theta }{S}_{\varphi }& T_{\theta }{C}_{\varphi }\\ 0& C_{\varphi }& -S_{\varphi }\\ 0& S_{\varphi }/C_{\theta }& C_{\varphi }/C_{\theta }\end{array}\right]\left(\begin{array}{c}p\\ q\\ r\end{array}\right)$$

(23)

where $T_{(·)}$ denotes tan$(·)$.

By employing the Newton-Euler formulation, the rotational dynamic equations can be derived as

$$\begin{array}{c}\hfill \begin{array}{c}{J}_b\dot{\mathcal{W}}=-\mathcal{W}\times J_b\mathcal{W}-M_g-M_d+M_b\hfill \end{array}\end{array}$$

(24)

$$\begin{array}{c}\hfill \begin{array}{c}{M}_g=\sum _{i=1}^4{J}_r(\mathcal{W}\times e_3){(-1)}^{i+1}{\Omega }_i\hfill \end{array}\end{array}$$

(25)

$$\begin{array}{c}\hfill \begin{array}{c}{M}_d=diag\{d_{\varphi },d_{\theta },d_{\phi }\}\dot{\mathcal{E}}\hfill \end{array}\end{array}$$

(26)

$$\begin{array}{c}\hfill M_b=\left[\begin{array}{c}lb({\Omega }_4^2-{\Omega }_2^2)\\ lb({\Omega }_3^2-{\Omega }_1^2)\\ d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)\end{array}\right]\end{array}$$

(27)

where $J_b=diag\{J_x,J_y,J_z\}$; $M_g$ is the resultant torque; $M_d$ is the aerodynamic frictions torque; $M_b$ is the rotor torque; l is the distance between rotor and center of mass; d denotes the reverse moment coefficient; $J_r$ is the rotational inertia of each rotor; $J_x$,$J_y$,$J_z$ are the rotary inertia; and $d_{\varphi },d_{\theta },d_{\phi }$ are the drag coefficients.

Then, the following dynamic equations can be derived

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\dot{p}={\tau }_{1}qr-{\tau }_2\mathcal{O}q-{\tau }_{3}p+U_p\hfill \\ \dot{q}={\tau }_{4}pr+{\tau }_5\mathcal{O}p-{\tau }_{6}q+U_q\hfill \\ \dot{r}={\tau }_{7}pq-{\tau }_{8}r+U_r\hfill \end{array}\right.\end{array}$$

(28)

where

• ${\tau }_1=\frac{J_y-J_z}{J_x},{\tau }_2=\frac{J_r}{J_x},{\tau }_3=\frac{d_{\varphi }}{J_x},{\tau }_4=\frac{J_z-J_x}{J_y},$

• ${\tau }_5=\frac{J_r}{J_y},{\tau }_6=\frac{d_{\theta }}{J_y},{\tau }_7=\frac{J_x-J_y}{J_z},{\tau }_8=\frac{d_{\phi }}{J_z},$

• $\mathcal{O}={\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4,U_p=lb({\Omega }_4^2-{\Omega }_2^2)/J_x,$

• $U_q=lb({\Omega }_3^2-{\Omega }_1^2)/J_y,U_r=d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)/J_z.$

Consider a group of N quadrotor UAVs, define $m_{i,1}={\varphi }_i$, $m_{i,2}={\theta }_i$, $m_{i,3}={\phi }_i$, $n_{i,1}=p_i$, $n_{i,2}=q_i$, and $n_{i,3}=r_i$, then the following unified attitude system can be obtained

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{m}}_{i,l}=n_{i,l}+f_{i,l}+w_{m_{i,l}}\hfill \\ {\dot{n}}_{i,l}=U_{n_{i,l}}+{\mathsf{\Phi }}_{i,l}^T{\xi }_{i,l}(n_i)+w_{n_{i,l}},l=1,2,3\hfill \end{array}\right.\end{array}$$

(29)

where $w_{m_{i,l}}$ and $w_{n_{i,l}}$ represent uncertain disturbances, and

• $f_{i,1}=T_{m_{i,2}}{S}_{m_{i,1}}{n}_{i,2}+T_{m_{i,2}}{C}_{m_{i,1}}{n}_{i,3},$

• $f_{i,2}=C_{m_{i,1}}{n}_{i,2}-S_{m_{i,1}}{n}_{i,3}-n_{i,2},$

• $f_{i,3}=S_{m_{i,1}}/C_{m_{i,2}}{n}_{i,2}+C_{m_{i,1}}/C_{m_{i,2}}{n}_{i,3}-n_{i,3},$

• ${\xi }_{i,1}(n_i)={[n_{i,2}{n}_{i,3},-{\mathcal{O}}_i{n}_{i,2},-n_{i,1}]}^T,$

• ${\xi }_{i,2}(n_i)={[n_{i,1}{n}_{i,3},{\mathcal{O}}_i{n}_{i,1},-n_{i,2}]}^T,$

• ${\xi }_{i,3}(n_i)={[n_{i,1}{n}_{i,2},-n_{i,3}]}^T,{\mathsf{\Phi }}_{i,1}={[{\tau }_1,{\tau }_2,{\tau }_3]}^T,$

• ${\mathsf{\Phi }}_{i,2}={[{\tau }_4,{\tau }_5,{\tau }_6]}^T,{\mathsf{\Phi }}_{i,3}={[{\tau }_7,{\tau }_8]}^T.$

 Assumption 3.

The uncertain disturbances satisfy

$$\begin{array}{c}\hfill |w_{m_{i,l}}|\le {\overline{w}}_m,|w_{n_{i,l}}|\le {\overline{w}}_n\end{array}$$

(30)

where ${\overline{w}}_m>0$ and ${\overline{w}}_n>0$ are positive constants.

 Remark 2.

Note that each UAV has to estimate the desired yaw angle ${\phi }_0$ by the information obtained from its neighbors. Inspired by [20], design the following distributed estimator

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\phi }}_{i,0}=-{ϵ}_1[\sum _{j=1}^N{\mathbf{a}}_{ij}({\phi }_{i,0}-{\phi }_{j,0})+{\mathbf{d}}_i({\phi }_{i,0}-{\phi }_0)]\hfill \\ -{ϵ}_2\mathrm{sgn}[\sum _{j=1}^N{\mathbf{a}}_{ij}({\phi }_{i,0}-{\phi }_{j,0})+{\mathbf{d}}_i({\phi }_{i,0}-{\phi }_0)]\hfill \end{array}\end{array}$$

(31)

where ${\phi }_{i,0}$ is an estimate of ${\phi }_0$; ${ϵ}_1>0$ and ${ϵ}_2>0$ are design parameters; $\mathrm{sgn}$ is the signum function. From Theorem 3.1 in [19], one can get that ${\phi }_{i,0}\to {\phi }_0$ in finite time.

3.2. Disturbance Rejection Attitude Controller Design

For the ith UAV, define two error variables

$$\begin{array}{c}\hfill \begin{array}{c}{\epsilon }_{i,l}=m_{i,l}-{\vartheta }_{i,l},{\rho }_{i,l}=n_{i,l}-{\beta }_{i,l},l=1,2,3\hfill \end{array}\end{array}$$

(32)

where ${\vartheta }_{i,1}={\varphi }_{i,0}$, ${\vartheta }_{i,2}={\theta }_{i,0}$, and ${\vartheta }_{i,3}={\phi }_{i,0}$; and ${\beta }_{i,l}$ are the virtual control functions. From (22) and (31), there exists an unknown constant $\overline{\vartheta }>0$ such that $|{\dot{\vartheta }}_{i,l}|\le \overline{\vartheta }$. The detailed design procedure is given as follows:

Step 1: The derivative of ${\epsilon }_{i,l}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\epsilon }}_{i,l}=n_{i,l}+f_{i,l}+w_{m_{i,l}}-{\dot{\vartheta }}_{i,l}.\hfill \end{array}\end{array}$$

(33)

The virtual control function ${\beta }_{i,l}$ is chosen as

$$\begin{array}{c}\hfill \begin{array}{c}{\beta }_{i,l}=-k_{m_{i,l}}{\epsilon }_{i,l}-f_{i,l}-\frac{{\epsilon }_{i,l}}{\sqrt{{\epsilon }_{i,l}^2+{\delta }_{i,l}^2}}{\widehat{\mu }}_{m_{i,l}}\hfill \end{array}\end{array}$$

(34)

where $k_{m_{i,l}}>0$ is a design constant; ${\widehat{\mu }}_{m_{i,l}}$ is the estimate of ${\mu }_{m_{i,l}}={\overline{w}}_m+\overline{\vartheta }$. Consider the Lyapunov function

$$\begin{array}{c}\hfill \begin{array}{c}{V}_{m_{i,l}}=\frac{1}{2}{\epsilon }_{i,l}^2+\frac{1}{2{\lambda }_{m_{i,l}}}{\tilde{\mu }}_{m_{i,l}}^2\hfill \end{array}\end{array}$$

(35)

where the estimation errors ${\tilde{\mu }}_{m_{i,l}}={\mu }_{m_{i,l}}-{\widehat{\mu }}_{m_{i,l}}$; and ${\lambda }_{m_{i,l}}>0$ is a design parameter.

From (33)-(35), the derivative of $V_{m_{i,l}}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{m_{i,l}}\le -k_{m_{i,l}}{\epsilon }_{i,l}^2+(|{\epsilon }_{i,l}|-\frac{{\epsilon }_{i,l}^2}{\sqrt{{\epsilon }_{i,l}^2+{\delta }_{i,l}^2}}){\mu }_{m_{i,l}}\hfill \\ +{\epsilon }_{i,l}{\rho }_{i,l}+\frac{1}{{\lambda }_{m_{i,l}}}{\tilde{\mu }}_{m_{i,l}}(\frac{{\lambda }_{m_{i,l}}{\epsilon }_{i,l}^2}{\sqrt{{\epsilon }_{i,l}^2+{\delta }_{i,l}^2}}-{\dot{\widehat{\mu }}}_{m_{i,l}}).\hfill \end{array}\end{array}$$

(36)

The parameter update law is chosen as

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\widehat{\mu }}}_{m_{i,l}}=\frac{{\lambda }_{m_{i,l}}{\epsilon }_{i,l}^2}{\sqrt{{\epsilon }_{i,l}^2+{\delta }_{i,l}^2}}.\hfill \end{array}\end{array}$$

(37)

Then, by applying Lemma 2, we have

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{h_l}\le \sum _{i=1}^N(-k_{m_{i,l}}{\epsilon }_{i,l}^2+{\epsilon }_{i,l}{\rho }_{i,l}+{\delta }_{i,l}{\mu }_{m_{i,l}}).\hfill \end{array}\end{array}$$

(38)

Step 2: From (29) and (32), the derivative of ${\rho }_{i,l}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\rho }}_{i,l}=U_{n_{i,l}}+{\mathsf{\Phi }}_{i,l}^T{\xi }_{i,l}(n_i)+w_{n_{i,l}}-\frac{\partial {\beta }_{i,l}}{\partial {\delta }_{i,l}}{\dot{\delta }}_{i,l}\hfill \\ -\frac{\partial {\beta }_{i,l}}{\partial {\widehat{\mu }}_{m_{i,l}}}{\dot{\widehat{\mu }}}_{m_{i,l}}-\frac{\partial {\beta }_{i,l}}{\partial m_{i,l}}(n_{i,l}+f_{i,l}+w_{m_{i,l}}).\hfill \end{array}\end{array}$$

(39)

The attitude controller is designed as

$$\begin{array}{c}\hfill \begin{array}{c}{U}_{n_{i,l}}=-k_{n_{i,l}}{\rho }_{i,l}-{\epsilon }_{i,l}-{\widehat{\mathsf{\Phi }}}_{i,l}^T{\xi }_{i,l}(n_i)+\frac{\partial {\beta }_{i,l}}{\partial {\widehat{\mu }}_{m_{i,l}}}{\dot{\widehat{\mu }}}_{m_{i,l}}\hfill \\ +\frac{\partial {\beta }_{i,l}}{\partial {\delta }_{i,l}}{\dot{\delta }}_{i,l}+\frac{\partial {\beta }_{i,l}}{\partial m_{i,l}}(n_{i,l}+f_{i,l})-\frac{{\rho }_{i,l}{\psi }_{i,l}^2}{\sqrt{{\rho }_{i,l}^2{\psi }_{i,l}^2+{\delta }_{i,l}^2}}{\widehat{\mu }}_{n_{i,l}}.\hfill \end{array}\end{array}$$

(40)

where $k_{n_{i,l}}>0$ is a design constant; ${\widehat{\mu }}_{n_{i,l}}$ is the estimate of ${\mu }_{n_{i,l}}=max\{{\overline{w}}_m,{\overline{w}}_n,\overline{\vartheta }\}$; and ${\psi }_{i,l}=\sqrt{1+{(\frac{\partial {\beta }_{i,l}}{\partial m_{i,l}})}^2}.$

Construct the following Lyapunov function

$$\begin{array}{c}\hfill \begin{array}{c}{V}_{n_{i,l}}=V_{m_{i,l}}+\frac{1}{2}({\rho }_{i,l}^2+\frac{1}{{\lambda }_{n_{i,l}}}{\tilde{\mu }}_{n_{i,l}}^2+\frac{1}{{\lambda }_{{\mathsf{\Phi }}_{i,l}}}{\tilde{\mathsf{\Phi }}}_{i,l}^T{\tilde{\mathsf{\Phi }}}_{i,l})\hfill \end{array}\end{array}$$

(41)

where the estimation errors ${\tilde{\mu }}_{n_{i,l}}={\mu }_{n_{i,l}}-{\widehat{\mu }}_{n_{i,l}}$ and ${\tilde{\mathsf{\Phi }}}_{i,l}={\mathsf{\Phi }}_{i,l}-{\widehat{\mathsf{\Phi }}}_{i,l}$; ${\lambda }_{n_{i,l}}>0$ and ${\lambda }_{{\mathsf{\Phi }}_{i,l}}>0$ are design parameters.

From (38)-(41), the derivative of $V_{n_{i,l}}$ satisfies

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{V}}_{n_{i,l}}\le -k_{m_{i,l}}{\epsilon }_{i,l}^2+(|{\rho }_{i,l}|{\psi }_{i,l}-\frac{{\rho }_{i,l}^2{\psi }_{i,l}^2}{\sqrt{{\rho }_{i,l}^2{\psi }_{i,l}^2+{\delta }_{i,l}^2}}){\mu }_{n_{i,l}}\hfill \\ -k_{n_{i,l}}{\rho }_{i,l}^2+\frac{1}{{\lambda }_{n_{i,l}}}{\tilde{\mu }}_{n_{i,l}}(\frac{{\lambda }_{n_{i,l}}{\rho }_{i,l}^2{\psi }_{i,l}^2}{\sqrt{{\rho }_{i,l}^2{\psi }_{i,l}^2+{\delta }_{i,l}^2}}-{\dot{\widehat{\mu }}}_{n_{i,l}})\hfill \\ +\frac{1}{{\lambda }_{{\mathsf{\Phi }}_{i,l}}}{\tilde{\mathsf{\Phi }}}_{i,l}^T({\lambda }_{{\mathsf{\Phi }}_{i,l}}{\rho }_{i,l}{\xi }_{i,l}(n_i)-{\dot{\widehat{\mathsf{\Phi }}}}_{i,l})+{\delta }_{i,l}{\mu }_{m_{i,l}}.\hfill \end{array}\end{array}$$

(42)

The parameter update laws are chosen as

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{\widehat{\mu }}}_{n_{i,l}}=\frac{{\lambda }_{n_{i,l}}{\rho }_{i,l}^2{\psi }_{i,l}^2}{\sqrt{{\rho }_{i,l}^2{\psi }_{i,l}^2+{\delta }_{i,l}^2}},{\dot{\widehat{\mathsf{\Phi }}}}_{i,l}={\lambda }_{{\mathsf{\Phi }}_{i,l}}{\rho }_{i,l}{\xi }_{i,l}(n_i).\hfill \end{array}\end{array}$$

(43)

Now, we present the analysis results.

 Theorem 2.

Consider the quadrotor UAV attitude system (29), the attitude controller (40), and the adaptive laws (37) and (43). All the signals in the closed-loop system are globally bounded, and the tracking error of the attitude angle system can converge to zero.

$Proof:$ The proof is similar to Theorem 1.

 Remark 3.

The proposed distributed formation tracking control scheme does not require the use of the quadrotor model parameters. Therefore, the proposed scheme is significant for achieving distributed formation tracking control of heterogeneous quadrotor UAV swarms.

4. An Illustrative Example

In this section, we consider a swarm system consisting of five quadrotor UAVs, and the model parameters of quadrotor UAVs are borrowed from literature [18]. The communication topology among UAVs is shown in Figure 1. The desired flight trajectory are chosen as $\mathbf{n}(t)={[0.1t,0.001t^2,0.1t]}^T$, and the desired yaw angle ${\phi }_0=t$. The reference formation shape vectors are given by ${\mathbf{F}}_i(t)={[cos(\frac{2i\pi }{5}+\frac{\pi }{50}t),sin(\frac{2i\pi }{5}+\frac{\pi }{50}t),0]}^T$ $(i=1,\cdots ,5).$ The controller parameters are chosen as $k_{h_{i,l}}=1$, $k_{v_{i,l}}=1$$(l=1,2,3)$, ${\lambda }_{h_{i,l}}=0.01$, ${\lambda }_{v_{i,l}}=0.01$, ${\delta }_{i,l}(t)=0.1e^{-t}$, ${ϵ}_1=2$, ${ϵ}_2=2$, $k_{m_{i,l}}=1,$$k_{n_{i,l}}=1$, ${\lambda }_{m_{i,l}}=0.01$, and ${\lambda }_{n_{i,l}}=0.01$. In addition, some uncertain disturbances are considered as $w_{h_{i,l}}=w_{v_{i,l}}=0.25sin(t)cos(t)$ and $w_{m_{i,l}}=w_{n_{i,l}}=0.15co{s}^2(t)+0.25sin(t)$.

The quadrotor UAVs’ flight trajectories in the 3-D space are displayed in Figure 2. As can be seen from Figure 2, by applying the proposed control scheme, the five UAVs form a desired formation shape and track the desired flight trajectory. Figure 3 shows the reference formation shapes and the actual flight formation of quadrotor UAVs. The response curves of formation tracking errors are also shown in Figure 3. It can be seen that the formation tracking error of each UAV converges to zero.

Figure 3. Formation shapes and formation tracking errors.

Figure 3. Formation shapes and formation tracking errors.

5. Conclusion

In this paper, a distributed formation tracking control method has been proposed for quadrotor UAVs. The proposed control scheme does not need to use the model parameters of quadrotor UAVs, which has wider practicality. The effectiveness of the proposed method has been verified by a numerical example. Our future work includes time-varying formation tracking control of heterogeneous quadrotor UAVs under switched communication topologies.