NOMENCLATURE

$GMT$

Ground Moving Target

$UAV$

Unmanned Aerial Vehicle

$MPC$

Model Predictive Control

$FOV$

Field of View

$LQR$

Linear Quadratic Programing

$UGV$

Unmanned Ground Vehicle

$\phi, \theta, \psi$

Rotational Variables of Quadrotor

$a_\phi, a_\theta, a_\psi$

Rotational Acceleration of Quadrotor

$w_\phi, w_\theta, w_\psi$

Rotational Disturbance

$x, y, z$

Actual Translational Variables of Quadrotor

$\hat x,\hat y,\hat z$

Measured Translational Variables of Quadrotor

$w_x, w_y, w_z$

Translational Disturbance

$a_x, a_y, a_z$

Translational Acceleration of Quadrotor

$x_{mt}, y_{mt}, z_{mt}$

Actual Translational Variables of Moving Target

$\hat x_{mt},\hat y_{mt},\hat z_{mt}$

Measured Translational Variables of Moving Target

$a^{mt}_x, a^{mt}_y, a^{mt}_z$

Translational Acceleration of Moving Target

$e_{x}, e_{y}, e_{z}$

Target Tracking Errors

$\eta$

Rotational state vector

$T_\eta$

Sampling Time of Control Loop

$\xi$

Translational state vector

$T_\xi$

Sampling Time of Guidance Loop

$B_x,B_y,B_z$

Body-fixed Coordinate

$E_x,E_y,E_z$

Earth-fixed Coordinate

$T_i, i=1,\ldots,4$

Thrust Forces

$I_{xx},I_{yy},I_{zz}$

Moment of Inertia of the Quadrotor

$m$

Quadrotor's Mass

$l$

Length of Quadrotor's Arm

$g$

Gravitational Acceleration

$\textit{v}_i, i=1,\ldots,4$

Normalised Input Voltage of each Motor

$L^o_{FOV}$

Actual Length of Field of View

$L_{FOV}$

Length of Truncated Field of View

$L_C$

Length of Critical Region

$L_{TR}$

Length of Trusted Region

$z_r$

Altitude's Reference of Quadrotor

$\delta^{uav}_{x},\delta^{uav}_{y},\delta^{uav}_{z}$

Measuring Fault of Quadrotor's Position

$\delta^{mt}_{x},\delta^{mt}_{y},\delta^{mt}_{z}$

Measuring Fault of Target's Position

$\Delta_{x},\Delta_{y},\Delta_{z}$

Maximum Measuring Fault of Position

$\epsilon$

Constant Adjustable Parameter

$P_\mu$

Prediction Horizon

$N_\mu$

Control Horizon

$J_\mu$

Cost Function

$E_{a_\mu}$

Input Vector of Cost Function

$Y_{e_\mu}$

Output Vector of Cost Function

$R_\mu$

Penalty Factor of Inputs

$Q_\mu$

Penalty Factor of Outputs

2.0 QUADROTOR NONLINEAR DYNAMICS

It is assumed that the quadrotor has a rigid body and symmetrical structure. Also, the propellers of the motors are considered rigid which means that the flapping effects of propellers are negligible. Furthermore, the origin of the body-fixed frame coincides with the center of the gravity of the quadrotor and the thrust forces are proportional to the input voltage of the actuators. Two coordinate systems have been utilised including the body-fixed frame $ B=[B_{1},B_{2},B_{3}] $ and the Earth-fixed frame as illustrated in Fig. 1. The dynamic model of the quadrotor, which is used in the current paper, is presented in detail in Ref. (Reference Waslander, Hoffmann, Jang and Tomlin23). This dynamic model is explained briefly in the following.

By employing the Euler-Lagrange equation, the nonlinear dynamic model is derived in the general form of $ \dot{X}=F(X,U)+W$ , in which $ X \in R^{12},U \in R^{4} $ and $ W\in R^{12}$ are states, inputs, and aerodynamic disturbances, respectively^{(Reference Alexis, Nikolakopoulos and Tzes24)}. The state vector consists of the quadrotor mass center position and the corresponding linear velocities of the quadrotor, i.e. $\xi=[z,\dot{z},x,\dot{x},y,\dot{y}]$ , which is expressed in the Earth-fixed frame and the system attitude state vector, which is shown as $\eta=[\phi,\dot{\phi},\theta,\dot{\theta},\psi,\dot{\psi}]$ .

According to the nature of the quadrotor dynamic model, it is usual to consider two separate subsystems, named rotational and translational subsystems containing $\xi$ and $\eta$ variables as it is mentioned by Bouabdallah et al.^{(Reference Bouabdallah and Siegwart25)}. It should be noted that this consideration remains admissible for the relatively slow translational velocities.

Based on these descriptions, the state space representation of the rotational subsystem dynamic model is as follows^{(Reference Alexis, Nikolakopoulos and Tzes24)}.

(1) \begin{equation}\begin{bmatrix}\dot{\eta}_1 \\[3pt]\dot{\eta}_2 \\[3pt]\dot{\eta}_3 \\[3pt]\dot{\eta}_4 \\[3pt]\dot{\eta}_5 \\[3pt]\dot{\eta}_6 \\v\end{bmatrix}=\begin{bmatrix}\eta_2\\[3pt]a_1\eta_6\eta_4+b_1(T_3-T_1)\\[3pt]\eta_4\\[3pt] a_2\eta_6\eta_2+b_2(T_4-T_2)\\[3pt] \eta_6\\[3pt] a_3\eta_4\eta_2+b_3(T_4+T_2-T_3-T_1)\end{bmatrix}+\begin{bmatrix}0\\[3pt] w_{\phi}\\[3pt] 0\\[3pt] w_\theta\\[3pt] 0\\[3pt] w_\psi \end{bmatrix}\end{equation}

Table 1 Quadrotor model parameters^{(Reference Alexis, Nikolakopoulos and Tzes24)}

Figure 1. Body and earth coordinate systems used for quadrotor dynamic modeling.

Also, the state space representation of the translational subsystem dynamic model is shown in Equation(2).

(2) \begin{equation}\begin{bmatrix}\dot{\xi}_1 \\[4pt]\dot{\xi}_2 \\[4pt]\dot{\xi}_3 \\[4pt]\dot{\xi}_4 \\[4pt]\dot{\xi}_5 \\[4pt]\dot{\xi}_6 \\[4pt]\end{bmatrix}=\begin{bmatrix}\xi_2\\[3pt]-g+\cos\eta_1 \cos\eta_3 \sum_{i=1}^{4}T_i/{m}\\[3pt]\xi_4\\[3pt] (\cos\phi \sin\theta \cos\psi+ \sin\phi \sin\psi)\sum_{i=1}^{4}T_i/{m}\\[3pt] \xi_6\\[3pt] (\cos\phi \sin\theta \sin\psi - \sin\phi \cos\psi)\sum_{i=1}^{4}T_i/{m}\end{bmatrix}+\begin{bmatrix}0\\[3pt] w_{z}\\[3pt] 0\\[3pt] w_x\\[3pt] 0\\[3pt] w_y\end{bmatrix}\end{equation}

The parameters used in the above equations are described in Table 1. Other parameters used for simplification in equations presentation are defined as follows.

\begin{equation*}\begin{array}{l}a_1=(I_{yy}-I_{zz})/I_{xx}, a_2=(I_{zz}-I_{xx})/I_{yy}\\[3pt] a_3=(I_{xx}-I_{yy})/I_{zz}, b_1=l_a/I_{xx}, b_2=l_a/I_{yy}\\[3pt] b_3=1/I_{zz}.\end{array}\end{equation*}

To apply the practical consideration, thrust forces are modeled as the following form.

(3) \begin{equation}T_i=k_1v_i ;\,\,\,\,\,\,\, {\rm with}\,\, 0\le v_i\le 1\,\, {\rm and}\,\, i=1,2,3,4\end{equation}

Also, the duty cycle of the actuators are considered using (4).

(4) \begin{equation}D_i=(h_0v_i+h_1)/T_{pwm} ;\,\,\,\,\,\,\, {\rm for} \,\, i=1,2,3,4\end{equation}

In Equation (3), $v_i$ is normalised input voltage that corresponds to ith motor, which is typically in the range of zero to one, $T_{pwm}$ is the duty cycle period, $h_0$ and $h_1$ are parameters relate to the dead zone and the saturation limit of the actuators, and $k_1$ is motor parameter obtained by a simple linear identification.

Substituting Equation (3) in Equations (2) and (1) yields the dynamic model of the system in the following form.

(5) \begin{equation}\begin{bmatrix}\dot{\eta}_1 \\[4pt]\dot{\eta}_2 \\[4pt]\dot{\eta}_3 \\[4pt]\dot{\eta}_4 \\[4pt]\dot{\eta}_5 \\[4pt]\dot{\eta}_6 \\[4pt]\end{bmatrix}=\begin{bmatrix}\eta_2\\[4pt]a_1\eta_6\eta_4+b_1k_1u_{\phi}\\[4pt]\eta_4\\[4pt] a_2\eta_6\eta_2+b_2k_2u_{\theta}\\[4pt] \eta_6\\[4pt] a_3\eta_4\eta_2+b_3k_3u_{\psi}\end{bmatrix}+\begin{bmatrix}0\\[4pt] w_{\phi}\\[4pt] 0\\[4pt] w_\theta\\[4pt] 0\\[4pt] w_\psi\end{bmatrix}\\[4pt] \end{equation}

(6) \begin{equation}\begin{bmatrix}\dot{\xi}_1 \\[4pt]\dot{\xi}_2 \\[4pt]\dot{\xi}_3 \\[4pt]\dot{\xi}_4 \\[4pt]\dot{\xi}_5 \\[4pt]\dot{\xi}_6 \\[4pt]\end{bmatrix}=\begin{bmatrix}\xi_2\\[4pt]-g+k_1\cos\eta_1 \cos\eta_3 u_z/{m}\\[4pt]\xi_4\\[4pt]k_1 u_xu_z/{m}\\[4pt]\xi_6\\[4pt]k_1 u_yu_z/{m}\end{bmatrix}+\begin{bmatrix}0\\[4pt]w_{z}\\[4pt]0\\[4pt] w_x\\[4pt] 0\\[4pt] w_y \end{bmatrix}\end{equation}

in which

(7) \begin{equation}\begin{array}{ll}u_x= \cos\phi \sin\theta \cos\psi+ \sin\phi \sin\psi \\[3pt]u_y= \cos\phi \sin\theta \sin\psi - \sin\phi \cos\psi\\[3pt]u_z=\sum_{i=1}^{4}v_i/{m}\end{array}\end{equation}

$u_x$ and $u_y$ are considered as virtual inputs while $u_{z}$ and $u_{\phi},$ $u_{\theta},$ $u_{\psi}$ are real inputs that correspond to the total thrust and the attitude of the quadrotor. The real inputs of the system are defined by the following equations.

\begin{equation*}\begin{bmatrix}v_1\\[4pt]v_2\\[4pt]v_3\\[4pt]v_4\end{bmatrix}=T^{-1}\begin{bmatrix}u_{\phi}\\[4pt]u_{\theta}\\[4pt]u_{\psi}\\[4pt] u_z\end{bmatrix}; \,\,\,\,\,\,T=\begin{bmatrix}-1&\quad 0&\quad 1&\quad 0\\[4pt] 0&\quad -1&\quad 0&\quad 1\\[4pt] -1&\quad 1&\quad -1 &\quad 1\\[4pt] 1&\quad 1&\quad 1&\quad 1\end{bmatrix}\end{equation*}

It should be mentioned that the orientations of the quadrotor are restricted as follows.

\begin{equation*}\begin{bmatrix}-\pi/4\\[4pt]-\pi/4\\[4pt]-\pi\end{bmatrix}<\begin{bmatrix}\phi\\[4pt]\theta\\[4pt]\psi\end{bmatrix}<\begin{bmatrix}\pi/4\\[4pt]\pi/4\\[4pt]\pi\end{bmatrix}\end{equation*}

The dynamic model of the ground moving target is considered as follows, which corresponds to the accelerated three-dimensional motion of the target on the ground

(8) \begin{equation}\begin{bmatrix}\dot{\xi}^{mt}_1 \\[5pt]\dot{\xi}^{mt}_2 \\[5pt]\dot{\xi}^{mt}_3 \\[5pt]\dot{\xi}^{mt}_4 \\[5pt]\dot{\xi}^{mt}_5 \\[5pt]\dot{\xi}^{mt}_6 \\[5pt]\end{bmatrix}=\begin{bmatrix}\xi^{mt}_2\\[4pt]0\\[4pt]\xi^{mt}_4\\[4pt]0\\[4pt]\xi^{mt}_6\\[4pt]0\end{bmatrix}+\begin{bmatrix}0\\[4pt]a^{mt}_z\\[4pt]0\\[4pt]a^{mt}_x\\[4pt]0\\[4pt]a^{mt}_y\end{bmatrix}\end{equation}

in which, $\xi^{mt}=[z_{mt},\dot{z}_{mt},x_{mt},\dot{x}_{mt},y_{mt},\dot{y}_{mt}] $ .

3.0 PROBLEM FORMULATION

As mentioned earlier, the possible altitude variation of the quadrotor can be used to prevent GMT loss.This section presents a guidance law and a sufficient condition to adjust the minimum quadrotor altitude for preserving the target in the FOV in long-term tracking. For this purpose, the FOV of the quadrotor is divided into two critical and trusted regions. The dimensions of the trusted region depend on the controller performance that is evaluated by the distance of the quadrotor from the target. The critical region is a predefined margin defined by an operator to prevent target loss in situations in which the target speed is temporally faster than the quadrotor speed. Also, it depends on the upper bounds of the quadrotor navigation errors and the target localisation errors. The proposed guidance law generates the minimum desired altitude of the quadrotor to adjust the FOV of the quadrotor to the union of the trusted and the critical regions.

Figure 2. Illustration of the parameters used in the problem formulation.

Let P and $\hat{P}$ be considered as the real and the estimated position vector, respectively. mt and uav subscriptions are used through the paper to indicate the moving target and the unmanned aerial vehicle.

\begin{equation*}P_{uav}=\begin{bmatrix} x& y& z \end{bmatrix}, P_{mt}=\begin{bmatrix}x_{mt} & y_{mt} & z_{mt} \end{bmatrix},\end{equation*}

\begin{equation*}\hat{P}_{uav}=\begin{bmatrix}\hat{x}& \hat{y}& \hat{z}\end{bmatrix}, \hat{P}_{mt}=\begin{bmatrix}\hat{x}_{mt} & \hat{y}_{mt} &\hat{z}_{mt}\end{bmatrix}\end{equation*}

Definition 1. Definition 1. By considering $z_r$ as the desired altitude reference; $e_z$ , $e_x$ and $e_y$ , i.e. the estimated tracking errors, are defined by the following equations.

(9) \begin{equation}\begin{array}{l}e_z=\hat{z}-z_r\\[3pt]e_x=\hat{x}-\hat{x}_{mt}\\[3pt]e_y=\hat{y}-\hat{y}_{mt}.\end{array}\end{equation}

The uncertainty of the position measurements, which are denoted by $\delta^{uav}_i$ and $\delta^{mt}_i$ for the quadrotor and moving target positions, respectively, are shown in Equation (10). Without loss of generality, the upper bounds of these uncertainties are considered constant.

(10) \begin{equation}\begin{array}{l}\hat{x}_{mt}=x_{mt}\pm \delta^{mt}_x\\[4pt]\hat{y}_{mt}=y_{mt}\pm \delta^{mt}_y\\[4pt]\hat{z}_{mt}=z_{mt}\pm \delta^{mt}_z\\[4pt]\hat{x}=x\pm \delta^{uav}_x\\[4pt]\hat{y}=y\pm \delta^{uav}_y\\[4pt]\hat{z}=z\pm \delta^{uav}_z.\\[4pt]\end{array}\end{equation}

As illustrated in Fig. 2, $e_{x-y}$ is the target tracking error in the horizontal plan and $L_{FOV}^o$ is the actual size of the smaller dimension of the field of view (FOV). To simplify the derivation and improvement of the proposed approach, a square FOV with $L_{FOV}$ dimensions is defined as a subset of the actual FOV which results the following relation;

(11) \begin{equation}L_{FOV}\le L^o_{FOV}\end{equation}

The considered square FOV is divided into critical and trusted regions. The dimensions of these regions are shown with $L_C$ and $L_{TR}$ . Based on this description, it can be written;

(12) \begin{equation}L_{FOV}(\hat{z})=L_{TR} (\hat{z})+L_C\end{equation}

As denoted in Equation (12), the dimensions of the critical region remain constant through an operation, which are set by the operator. However, $L_{FOV}$ and $L_{TR}$ are functions of the altitude of the quadrotor. In fact, $L_{FOV}$ will be equal to the $L_C$ when the tracking error ( $e_{x-y}$ in Fig. 2) becomes zero (13).

(13) \begin{equation}L_{FOV}(z_0)=L_C.\end{equation}

In Equation (13), $z_0$ is the required altitude of the quadrotor to provide a FOV with $L_C$ dimension.

The following assumptions are considered to develop the proposed strategy:

3.1 Critical region design

$L_C$ is a constant value composed of $\epsilon$ and $\Delta_{x-y-z}$ as follows.

(14) \begin{equation}L_C=\epsilon+\Delta_{x-y-z};\quad \epsilon>0,\quad \Delta_{x-y-z}>0\end{equation}

In the Equation (14), $\epsilon$ is used to prevent the target loss when the target speed is temporarily faster than the quadrotor speed. Its value is designed by the operator based on the maximum estimated duration and the amount of difference between the quadrotor and the target speed. Also, $\Delta_{x-y-z}$ is an upper bound for the effect of measurement errors on the horizontal plane, which is defined by the following equation

(15) \begin{equation}\Delta_{x-y-z}=2*(\Delta_z+\Delta_{x-y}),\end{equation}

where $\Delta_{z}$ is an equivalent value for the quadrotor and target al.titude measurement errors ( $\bar{\Delta}_{z}$ ) on the horizontal plane containing the target.

\begin{equation*}\begin{array}{l l}\Delta_{x-y}=\max({|\bar{\delta}^{uav}_x|+|\bar{\delta}^{mt}_x|,|\bar{\delta}^{uav}_y|+|\bar{\delta}^{mt}_y|}).\\[5pt]\bar{\Delta}_{z}=|\bar{\delta}^{mt}_z|+|\bar{\delta}^{uav}_z|\end{array}\end{equation*}

As previously mentioned about Equation (10), $\bar{\delta}^{uav}_i$ and $\bar{\delta}^{mt}_i$ are the upper bounds of the quadrotor and the moving target localisation errors, respectively. i in this parameters denote the x,y,z directions.

It should be mentioned that the average speed of the quadrotor should be bigger than the average speed of the target in a long-term tracking scenario. The proposed algorithm is designed to be robust to the temporary increase in the target velocity in relation to the quadrotor velocity.

3.2 Trusted region design

$L_{TR}$ is an ascending function of the quadrotor altitude that means increasing in the quadrotor altitude causes the larger dimentions of $L_{TR}$ . Considering the relation between $L_{FOV}(z)$ , and the critical and the trusted regions dimensions mentioned in Equation (12) and the Assumption 2, a suitable function between $L_{TR}$ and the quadrotor altitude that satisfies the marginal constraint in Equation (13) can be as follows.

(16) \begin{equation}L_{TR} (\hat{z})=\left\{\begin{array}{l@{\quad}l}\hat{z}(k)-z_0 (k) & \hat{z}(k)>z_0 (k) \\[5pt] 0 & Else\end{array} \right.\end{equation}

A sufficient condition to preserve the target in the FOV will be proven through the three following lemmas. In the thired one, the desired altitude of the quadrotor is determined based on the designed trusted region. Also, a sufficient condition for the loss of target during the ususal target tracking without the FOV adjustment will be proven in the last lemma.

Lemma 1. Lemma 1. To satisfy the inequality (11), it is enough to determine $z_0$ as

\begin{equation*}z_0\ge L_C+\hat{z}_{mt}.\end{equation*}

Proof. Substituting Equation (16) in Equation (12) takes the following form.

(17) \begin{equation}L_{FOV}(\hat{z})=\hat{z}(k)-z_0+L_C=\hat{z}(k)-\hat{z}_{mt}-(z_0-\hat{z}_{mt}-L_C)\end{equation}

Considering the Assumption 2 and the $z_0$ determination, which means $z_0-\hat{z}_{mt}-L_C$ is not negative, (17) can be rewritten as the following inequality.

\begin{equation*}L_{FOV}(\hat{z})\le \hat{z}(k)-\hat{z}_{mt}\le L^o_{FOV}(k).\end{equation*}

Lemma 2. Lemma 2. The real and the estimated distances between the quadrotor and the target satisfy the following inequalities.

\begin{equation*}\begin{array}{l}|x-x_{mt}|<|e_x|+\Delta_{x-y} \\[4pt]|y-y_{mt}|<|e_y|+\Delta_{x-y}\\[4pt]|x-x_{mt}|>|e_x|-\Delta_{x-y} \\[4pt]|y-y_{mt}|>|e_y|-\Delta_{x-y}\\[4pt]\end{array}\end{equation*}

Proof. Substituting Equation (10) in Equation (9), the following equations are obtained.

(18) \begin{equation}\begin{array}{l}e_x=x-x_{mt}\pm \delta^{uav}_x \pm \delta^{mt}_x \\[5pt]e_y=y-y_{mt}\pm \delta^{uav}_y \pm \delta^{mt}_y,\end{array}\end{equation}

Considering the definition of $\Delta_{x-y}$ in the critical region design, it is straightforward to derive the Lemma 2 from Equation (18).

Lemma 3. Lemma 3. A sufficient condition to maintain the target in the FOV is defined as follows.

(19) \begin{equation}e_z(k)+\epsilon>0\quad ;e_z=\hat{z}-z_r\end{equation}

which is provided by

(20) \begin{equation}z_r(k)=2 \max (|e_x|,|e_y|)+z_0 (k)\end{equation}

Proof. To preserve the target in the field of view of the quadrotor, the following inequality should be satisfied.

(21) \begin{equation}L_{FOV}/2\gt\max(|x-x_{mt}|,|y-y_{mt}|), \quad for\,\, k\,\gt\,0\end{equation}

Applying the Lemma 2 to the inequality (21) results the following inequality.

(22) \begin{equation}L_{FOV}/2>\max(|e_x|+\Delta_{x-y},|e_y|+\Delta_{x-y})\end{equation}

Substituting Equations (15) and (12) in Equation (22) yields the following relation.

(23) \begin{equation}L_{TR}(\hat{z})/2+\epsilon/2>\max(|e_x|,|e_y|)\end{equation}

Defining $e_{L_{TR}}=L_{TR}(\hat{z})-L_{TR}(z_r)$ , Equation (23) takes the following form.

(24) \begin{equation}e_{L_{TR}}/2+\epsilon/2>\max(|e_x|,|e_y|)-L_{TR}(z_r)/2\end{equation}

To eliminate the right side of the inequality (24), $L_{TR}(z_r)$ is designed using the following form.

(25) \begin{equation}L_{TR}(z_r)=2\max(|e_x|,|e_y|).\end{equation}

According to Equations (16) and (25), the desired altitude of the quadrotor to preserve the target in the FOV is determined as follows.

(26) \begin{equation}z_r(k)=2\max(|e_x|,|e_y|)+z_0 (k)\end{equation}

Using Equations (16) and (26), $e_{L_{TR}}$ changes to the following form.

(27) \begin{equation}e_{L_{TR}}=\hat{z}(k)-z_0(k)-z_r(k)+z_0(k)=e_z(k)\end{equation}

Substituting Equations (25) and (27) in Equation (24), the Lemma 3 inequality is obtained.

Lemma 4. Lemma 4. Considering the assumption $ L^o_{FOV}= L_{FOV}$ , a sufficient condition for the loss of the target using the usual target tracking approach without the FOV adjustment by the altitude increasing is defined as follows.

(28) \begin{equation}e_z(k)+\epsilon-2\min (|e_x|,|e_y|) \lt 0;\quad e_z=\hat{z}-z_r\end{equation}

which is provided by

(29) \begin{equation}z_r(k)=z_0 (k)\end{equation}

Proof. To lose the target from the field of view of the quadrotor, the following inequality should be satisfied.

(30) \begin{equation} L_{FOV}/2<\min (|x-x_{mt}|,|y-y_{mt}|), \,\,\,\,\,\,\, for\,\, k\,>\,0\end{equation}

Applying the Lemma 2 to the inequality (30), the following inequality is obtained.

(31) \begin{equation} L_{FOV}/2<\min (|e_x|-\Delta_{x-y},|e_y|-\Delta_{x-y})\end{equation}

Substituting Equations (15) and (12) in Equation (31) yields the following relation.

(32) \begin{equation}L_{TR}(\hat{z})/2+\epsilon/2<\min (|e_x|,|e_y|);\quad \epsilon\pm\Delta_z\approx \epsilon\end{equation}

Defining $e_{L_{TR}}=L_{TR}(\hat{z})-L_{TR}(z_r)$ , Equation (32) takes the following form.

(33) \begin{equation}e_{L_{TR}}/2+\epsilon/2<\min (|e_x|,|e_y|)-L_{TR}(z_r)/2\end{equation}

In the usuall approaches without the proposed strategy in the inequality (33), $L_{TR}(z_r)$ is selected as follows.

(34) \begin{equation}L_{TR}(z_r)=0\end{equation}

According to Equations (16) and (34), the reference altitude of the quadrotor is determined as the follows.

(35) \begin{equation}z_r(k)=z_0 (k)\end{equation}

Using Equations (16) and (35), $e_{L_{TR}}$ changes to the following form.

(36) \begin{equation}e_{L_{TR}}=\hat{z}(k)-z_0(k)-z_r(k)+z_0(k)=e_z(k)\end{equation}

Substituting Equations (34) and (36) to Equation (33), yields the inequality of Lemma 4.

4.0 MPC-BASED GUIDANCE AND CONTROL

As shown in the previous section, to maintain the target in the FOV of the quadrotor, it is sufficient to satisfy the inequality of (19) when designing the guidance law. Given this fact, and without loss of generality, a feedback linearision (FL) is applied to the nonlinear model to achieve a decoupled linear description of the system. Furthermore, the position error vector integral is added to the model to achieve null steady errors in the presence of wind disturbances^{(Reference Raffo, Ortega and Rubio19)}. The proposed MPC-based guidance and control system consist of three main components, the quadrotor desired attitude generator, the reference $x-y$ planar motion generator, and the desired attitude generator in a cascade structure. It should be noted that such cascade architecture requires slower translational dynamics than rotational ones. The guidance and control diagram including the non-linear system dynamic and the actuator model is illustrated in Fig. 3.

Figure 3. Block diagram of the MPC-based guidance and control subsystem.

4.1 Linear decoupled dynamic of the quadrotor

Initially, a linear decoupled dynamics of the quadrotor must be derived. To meet this need, primitive nonlinear feedbacks for rotational and translational subsystems are considered in the following forms, respectively.

(37) \begin{align}u_{\phi} & =- (a_1\eta_4\eta_6-a_{\phi})/c_1\nonumber\\[3pt]u_{\theta} & =-(a_3\eta_2\eta_6-a_{\theta})/c_2\\[3pt]u_{\psi} & =-(a_5\eta_2\eta_4-a_{\psi})/c_3\nonumber\end{align}

(38) \begin{equation}\begin{array}{l}u_z={m}(g+a_z )/({k_1\cos\eta_1 \cos\eta_3})\\[4pt]u_x=m a_x/k_1u_z\\[4pt]u_y=m a_y/k_1u_z\\\end{array}\end{equation}

where, $V=[a_{\phi},a_{\theta},a_{\psi},a_z,a_x,a_y ] $ is the quadrotor movement acceleration vector.

Substituting Equation (37) in Equation (5) yields three decoupled identical linear dynamics for the rotational subsystem as shown in Equation (39).

(39) \begin{equation}\begin{array}{l}\dot\eta^s=A_m\eta^s+B_m a_{\eta}^s\quad \text{for}\ s=0,1,2\\[5pt]y_{\eta}^s=C_m \eta^s\end{array}\end{equation}

where

$\eta^0=[\phi,\dot{\phi} ]^T,\eta^1=[\theta,\dot{\theta} ]^T,\eta^2=[\psi,\dot{\psi}]^T$

$ A_m=\begin{bmatrix} 0&\quad 1 \\[5pt] 0&\quad 0 \end{bmatrix}$ $B_m=\begin{bmatrix} 0 \\[5pt] 1 \end{bmatrix}$ $C_m=\begin{bmatrix} 1 &0 \end{bmatrix}$

$a_{\eta}^0 = a_{\phi} $ $a_{\eta}^1 = a_{\theta} $ $a_{\eta}^2=a_{\psi} $

Also, substituting the Equation (38) into the Equation (6) yields three decoupled identical linear dynamics for the translational subsystem as follows.

(40) \begin{equation}\begin{array}{l}\dot\xi^s=A_m\xi^s+B_m a_{\xi}^s+B_m w_{\xi}^s\quad \text{for}\ s=0,1,2\\[5pt]y_{\xi}^s=C_m\xi^s\end{array}\end{equation}

in which

$\xi^0=[z,\dot{z} ]^T,\xi^1=[x,\dot{x}]^T,\xi^2=[y,\dot{y}]^T$

$\begin{bmatrix} a_{\xi}^0 \\[4pt]w_{\xi}^0 \end{bmatrix}=\begin{bmatrix} a_z \\[4pt]w_z \end{bmatrix}$ , $\begin{bmatrix} a_{\xi}^1 \\[4pt]w_{\xi}^1 \end{bmatrix}=\begin{bmatrix} a_x \\[4pt]w_x \end{bmatrix}$ , $\begin{bmatrix} a_{\xi}^2 \\[4pt]w_{\xi}^2 \end{bmatrix}=\begin{bmatrix} a_y \\[4pt]w_y \end{bmatrix}$

GMT dynamics in the Equation (8) can be written as

(41) \begin{equation}\dot\xi_{mt}^s=A_m\xi_{mt}^s+B_m a_{mt}^s \quad \text{for}\ s=0,1,2\end{equation}

in which

$\xi_{mt}^0=[z_{mt},\dot{z}_{mt} ]^T,\xi^1=[x_{mt},\dot{x}_{mt}]^T,\xi^2=[y_{mt},\dot{y}_{mt}]^T$

$a_{mt}^0=a_z^{mt}, a_{mt}^1=a_x^{mt}, a_{mt}^2=a_y^{mt}$

The virtual reference air vehicle defined in Ref. (Reference Raffo, Ortega and Rubio19) has a linear dynamic model similar to the quadrotor rotational subsystem.

(42) \begin{equation}\dot\eta_r^s=A_m\eta_r^s+B_ma_{\eta_r}^s \quad \text{for}\ s=0,1,2\end{equation}

where, $\eta_r^0=[\phi_r,\dot{\phi} _r]^T,\eta_r^1=[\theta_r,\dot{\theta}_r ]^T$ , and $ \eta_r^2=[\psi_r,\dot{\psi}_r]^T.$ Based on Equation (39), the linear acceleration reference values, i.e. $a_{\eta_r}^s$ can be defined as follows.

\begin{equation*}a_{\eta_r}^0=a_{\phi_r}=\ddot{\phi}_r, a_{\eta_r}^1=a_{\theta_r}=\ddot{\theta}_r, a_{\eta_r}^2=a_{\psi_r}=\ddot{\psi}_r\end{equation*}

To achieve zero steady state position errors, the augmented translational error vector is defined as follows.

(43) \begin{equation}\begin{array}{l}e_{\xi}^0=\begin{bmatrix}e_z \\[4pt] \dot{e}_z\\[4pt] \int{e_z}\end{bmatrix}=\begin{bmatrix}\hat{z}-z_r \\[4pt] \dot{z}-\dot{z}_{mt}-\dot{L}_c^{des}\\[4pt] \int{\hat{z}-z_r}\end{bmatrix}\\[30pt] e_{\xi}^1=\begin{bmatrix}e_x \\[4pt] \dot{e}_x\\[4pt] \int{e_x}\end{bmatrix}=\begin{bmatrix}\hat{x}-\hat{x}_{mt} \\[4pt] \dot{x}-\dot{x}_{mt}\\[4pt] \int{\hat{x}-\hat{x}_{mt}}\end{bmatrix}\\[30pt] e_{\xi}^2=\begin{bmatrix}e_y \\[4pt] \dot{e}_y\\[4pt] \int{e_y}\end{bmatrix}=\begin{bmatrix}\hat{y}-\hat{y}_{mt} \\[4pt] \dot{y}-\dot{y}_{mt}\\[4pt] \int{\hat{y}-\hat{y}_{mt}}\end{bmatrix},\end{array}\end{equation}

Substituting the Equations (40) and (41) in the Equation (43), the dynamics of the translational error is written as follows.

(44) \begin{equation}\begin{array}{l}\dot{e}_{\xi}^s(t)=A e_{\xi}^s+B e_{a_{\xi}}^s(t)+B w_{\xi}^s\ (t)\quad \text{for}\ s=0,1,2\\[3pt]e_{y_{\xi}}^s(t)=C e_{\xi}^s(t),\end{array}\end{equation}

in which

$e_{a_{\xi}}^0=\ddot{z}-\ddot{z}_r=a_z-a_z^{mt}-2\max(|e_{a_{\xi}}^1|,|e_{a_{\xi}}^2|)$ , $e_{a_{\xi}}^1=\ddot{x}-\ddot{x}_{mt}=a_x-a_x^{mt}$ , $$\eqalign{ & \quad e_{{a_\xi }}^2 = \ddot y - {{\ddot y}_{mt}} = \cr & {a_y} - a_y^{mt} \cr} $$

$A=\begin{bmatrix} 0&\quad 1&\quad 0\\[4pt]0&\quad 0&\quad 0\\[4pt]1&\quad 0&\quad 0 \end{bmatrix}, B=\begin{bmatrix} 0\\[4pt]1\\[4pt]0 \end{bmatrix}, C=\begin{bmatrix} 1&\quad 0&\quad 0 \end{bmatrix}$

By subtracting the Equation (39) from the Equation (42) and augmenting the attitude errors, the rotational error dynamics model is obtained as follows.

(45) \begin{align}\dot{e}_{\eta}^s(t) & =A e_{\eta}^s+B e_{a_{\eta}}^s(t)+B w_{\eta}^s\ (t)\quad \text{for}\ s=0,1,2\nonumber\\[5pt]e_{y_{\eta}}^s(t) & =C e_{\eta}^s(t),\end{align}

where, $e_{\eta}^s=\begin{bmatrix} \eta_1^s-{\eta_1}_r^s\\[5pt] \eta_2^s-{\eta_2}_r^s\\[5pt] \int{(\eta_1^s-{\eta_1}_r^s)} \end{bmatrix}$ , $e_{a_{\eta}}^s= a_{\eta}^s-a_{\eta_r}^s$

Considering the measurement rate of the sensors, $T_{\xi}$ and $T_{\eta}$ are defined as the sampling period for the guidance and control loops, respectively.

The discrete representation of the quadrotor error dynamics takes the following form^{(Reference Raffo, Ortega and Rubio19)} for $s=0,1,2$ .

(46) \begin{align} e_{\xi}^s(k+1)& =A_{\xi} e_{\xi}^s(k)+B_{\xi} e_{a_{\xi}}^s(k)+B_{\xi} w_{\xi}^s(k) \nonumber \\[5pt]e_{y_{\xi}}^s(k) & =C e_{\xi}^s(k),\end{align}

(47) \begin{align}e_{\eta}^s(k+1) & = A_{\eta} e_{\eta}^s(k)+B_{\eta} e_{a_{\eta}}^s(k)+B w_{\eta}^s(k) \nonumber \\[5pt]e_{y_{\eta}}^s(k) & = C e_{\eta}^s(k),\end{align}

in which

$A_{\xi}=\begin{bmatrix} 1&\quad T_{\xi}&\quad 0\\[5pt]0&\quad 1&\quad 0\\[5pt]T_{\xi}&\quad 0&\quad 1 \end{bmatrix}$ , $B_{\xi}=\begin{bmatrix} 0\\[5pt]T_{\xi}\\[5pt]0 \end{bmatrix}$ ,

$A_{\eta}=\begin{bmatrix} 1&\quad T_{\eta}&\quad 0\\[5pt]0&\quad 1&\quad 0\\[5pt]T_{\eta}&\quad 0&\quad 1 \end{bmatrix}$ , $B_{\eta}=\begin{bmatrix} 0\\[5pt]T_{\eta}\\[5pt]0 \end{bmatrix}$

Using $\mu$ to indicate both $\eta$ and $\xi$ subscriptions, the error dynamics, mentioned in Equations (46) and (47), can be written in a common form as follows.

(48) \begin{equation} \begin{array}{l}e_{\mu}^s(k+1)=A_{\mu} e_{\mu}^s(k)+B_{\mu} e_{a_{\mu}}^s(k) \\[4pt]e_{y_{\mu}}^s(k)=C e_{\mu}^s(k).\end{array}\quad {\rm for}\, \left\{ \begin{array}{l l} s=0,1,2 \\[4pt] \mu=\xi,\eta\end{array} \right.\end{equation}

4.2 Guidance and control synthesis

In the previous subsection, the guidance and control problem for a nonlinear high order system is simplified to a linear third-order model. In this section, the MPC formulation is discussed for the Equation (48), which is common for $s=0,1,2$ . It should be noted that unknown effects on the system such as unmeasurable disturbance are excluded from Equation (48).

$J_{\mu}$ is considered as a cost function that has the following form^{(Reference Sunan and Heng26)}.

(49) \begin{equation}J_{\mu}^s={Y_{e_{\mu}}^s}^TQ_{\mu}^sY_{e_{\mu}}^s+{E_{a_{\mu}}^s}^T R_{\mu}^s E_{a_{\mu}}^s\quad for\, \left\{ \begin{array}{l l} s=0,1,2 \\[4pt] \mu=\xi,\eta\end{array} \right.\end{equation}

in which

$Y_{e_{\mu}}^s=\begin{bmatrix}e_{y_{\mu}}^s(k+1|k) \\[4pt]e_{y_{\mu}}^s(k+2|k)\\[4pt] \vdots\\[4pt] e_{y_{\mu}}^s(k+P_{\mu}|k) \\[4pt] \end{bmatrix}$ $E_{a_{\mu}}^s=\begin{bmatrix}\hat{e}^s_{a_{\mu}}(k) \\[4pt] \hat{e}^s_{a_{\mu}}(k+1)\\[4pt] \vdots\\[4pt] \hat{e}^s_{a_{\mu}}(k+N_{\mu}-1) \\[4pt] \end{bmatrix},$

and the $P_{\mu}, N_{\mu} $ are the prediction and the control horizon parameters, respectively. Also, the $Q_{\mu}\geq 0 , R_{\mu}>0 $ are penalty factors related to the errors and the inputs, respectively. $\alpha_{\mu}$ is the pole of the first order filter on the reference signal to suppress its abrupt changes, as well as noise effects^{(Reference Camacho and Bordons27)}. The guidance law is obtained by solving the optimisation problem in the following form.

(50) \begin{equation}\begin{aligned}{E_{a_{\mu}}^s}^*=\ & \underset{E_{a_{\mu}}^s}{\text{minimize}} & & J_{\mu}^s\quad {\rm for} \left\{ \begin{array}{l l} s=0,1,2 \\[4pt] \mu=\xi,\eta\end{array} \right.\\[4pt] & \text{subject to} & & Y_{e_{\mu}}^s=F_{\mu}e_{\mu}^s(k)+G_{\mu}e_{a_{\mu}}^s(k)\\[4pt] & & & \underline{Y}_{e_{\mu}}^s \leq Y_{e_{\mu}}^s \leq \overline{Y}_{e_{\mu}}^s\\[4pt] & & & \underline{E}_{a_{\mu}}^s \leq E_{a_{\mu}}^s \leq \overline{E}_{a_{\mu}}^s \end{aligned}\end{equation}

where

\begin{equation*}\begin{array}{l}F_{\mu}=\begin{bmatrix}CA_{\mu}&\quad CA_{\mu}^2&\quad CA_{\mu}^3& \quad \ldots&\quad CA_{\mu}^P\end{bmatrix}^T\\[4pt]G_{\mu}=\begin{bmatrix}CB_{\mu}&\quad 0&\quad \ldots&\quad 0\\[4pt]CA_{\mu}B_{\mu}&\quad C_{\mu}B_{\mu}&\quad \ldots&\quad 0\\[4pt]CA_{\mu}^2B_{\mu}&\quad CA_{\mu}B_{\mu}&\quad \ldots&\quad 0\\[4pt] \vdots&&&\\[4pt] CA_{\mu}^{P-1}B_{\mu}&\quad CA_{\mu}^{P-2}B_{\mu}&\quad \ldots&\quad CA_{\mu}^{P-N}B_{\mu}\end{bmatrix}\end{array}\end{equation*}

Initial conditions of the dynamic model at instance $k{\rm th}$ is provided by the measurement of the system at the same time.

Optimisation methods such as quadratic programming can be used to solve the problem (50). In this approach, an optimisation should be performed at any time. The closed form is obtained by replacing of $ Y_{e_{\mu}}^s$ in the Equation (49) and derivation with respect to $e_{a_{\mu}}^s(k)$ , as the following equation^{(Reference Camacho and Bordons27)}:

(51) \begin{equation}{E_{a_{\mu}}^s}^*=-\left(G^T_{\mu}Q_{\mu}G_{\mu}+R_{\mu}\right)^{-1}G_{\mu}^TQ_{\mu}F_{\mu}e_{\mu}(k)\end{equation}

High computational burden and consumption time are prominent disadvantages of the optimisation methods. Given this, Equation (51) is employed as the solution of the optimisation problem (50) to derive the desired relative acceleration for the $x-y$ planar motion. This provides the admissible performance of the proposed algorithm. It should be noted that inequality constraints can be satisfied via adjusting $Q_{\mu}$ and ${R_{\mu}}$ parameters^{(Reference Camacho and Bordons27)}; however, in the considered case, it is assumed that inequality constraints remain passive using powerful actuators.

To obtain the desired relative acceleration for altitude, the inequality constraint (Reference Raffo, Ortega and Rubio19) should be satisfied. For this purpose, the following algorithm is proposed

• • Initially, the solution of Equation (51) is used to obtain $\hat{e}_{a_{\xi}}^0(k)$

• • If $e_z(k+1|k)+\epsilon>0 \to A_{\xi} e_{\xi}^0(k)+B_{\xi} \hat{e}_{a_{\xi}}^0(k)+\epsilon>0 $ then $(\hat{e}_{a_{\xi}}^0)^*=\hat{e}_{a_{\xi}}^0$

• • Else, numerical optimisation algorithms such as QP is employed for solving Equation (50) subject to $-\epsilon \begin{bmatrix}1&1&...&1 \end{bmatrix}< Y_{e_{\xi}}^0$ .

Subsequently, the appropriate translational acceleration, which reduces the distance between the quadrotor and the target, takes the following form.

(52) \begin{equation}\begin{array}{l}{a_{\xi}^{0}}^{des}=a_z^{des}(k)=\hat{a}_z^{mt}+\left(\hat{e}_{a_{\xi}}^0\right)^*+2\max\left(|e_{a_{\xi}}^1|,|e_{a_{\xi}}^2|\right)\\[8pt]{a_{\xi}^1}^{des}=a_x^{des}(k)=\hat{a}_x^{mt}+\left(\hat{e}_{a_{\xi}}^1\right)^*\\[8pt]{a_{\xi}^{2}}^{des}=a_y^{des}(k)=\hat{a}_y^{mt}+\left(\hat{e}_{a_{\xi}}^2\right)^*,\end{array}\end{equation}

where, $\hat{a}_z^{mt}$ , $\hat{a}_x^{mt}$ and $\hat{a}_y^{mt}$ are the measured acceleration of the moving target described by the following equations.

(53) \begin{equation}\begin{array}{l}\hat{a}_z^{mt}=a_z^{mt}\pm \delta_{a_z}^{mt}\\[7pt]\hat{a}_x^{mt}=a_x^{mt}\pm \delta_{a_x}^{mt}\\[7pt]\hat{a}_y^{mt}=a_y^{mt}\pm \delta_{a_y}^{mt},\end{array}\end{equation}

where $\delta_{a_z}^{mt}$ , $\delta_{a_z}^{mt}$ and $\delta_{a_z}^{mt}$ are estimation errors.

After calculating the desired translational acceleration, $u_z(k)$ and the reference of the roll and pitch angles, i.e. $\phi_r $ and $\theta_r $ , are derived using Equations (7) and (38). These references should be provided for the quadrotor rotational loop as follows.

(54) \begin{equation}\begin{array}{l}u_z(k)=m\dfrac{a_z^{des}(k)+g}{k_1\cos(\phi_k)\sin(\theta_k)}\\[14pt]\phi_r(k)=\arcsin\left(\dfrac{ma_x^{des}(k) }{k_1u_z(k)} \sin(\psi_k)-\dfrac{ma_y^{des}(k) }{k_1u_z(k)} \cos (\psi_k)\right)\\[14pt]\theta_r(k)=\arcsin\left(\dfrac {a_x^{des}(k) \cos(\psi_k)+a^{des}_y(k) \sin (\psi_k)} {\sqrt{\left(\dfrac{k_1u_z(k)}{m}\right)^2-\left(a^{des}_x(k) \sin(\psi_k)-a^{des}_y(k) \cos (\psi_k)\right)^2}}\right)\end{array}\end{equation}

Similarly, the desired rotational accelerations $a_\phi^{des}(k), a_\theta^{des}(k)$ and $a_\psi^{des}(k)$ have been computed to track attitude reference signals by the quadrotor.

(55) \begin{equation}{a_{\eta}^{s}}^{des}(k)=a_{\eta_r}^s(k)+\left(\hat{e}_{a_{\eta}}^s(k)\right)^* \,\,\,\,\,\,\, {\rm for}\,\, s=0,1,2\end{equation}

in which the $(\hat{e}_{a_{\eta}}^s(k))^*$ is the first component of ${E_{a_{\mu}}^s}^*$ in Equation (51).

General control inputs, i.e. $u_\phi(k), u_\theta(k)$ and $u_\psi(k)$ , can be generated by substituting Equation (55) in Equation (37).

The corresponding duty cycles of the quadrotor actuators are obtained using the Equation (4) presented in Section 4.

5.0 SIMULATION RESULTS

The MPC parameters presented in Table 2 are designed to achieve fast disturbance rejection and smooth output tracking. The sample time is selected to provide sufficient opportunity for the translational controller to track the predefined trajectory using the roll and pitch angles produced by the rotational controller.

Table 2 MPC-based guidance and control parameters

The initial condition is defined in such a way that satisfies Equation (19) at $k=0$ , i.e. $e_z(0)+\epsilon>0$ . The simulations are performed for the following initial condition of the translational and rotational subsystems, respectively.

$\begin{array}{l}\xi(0)=\begin{bmatrix} 200&\quad 0&\quad 50&\quad 0&\quad 50&\quad 0 \end{bmatrix}\\[5pt]\eta(0)=\begin{bmatrix} 0&\quad 0&\quad 0&\quad 0&\quad 0&\quad 0 \end{bmatrix}\end{array}$

Problem formulation parameters are listed in Table 3. The quadrotor parameters are selected based on a real model, which are as follows.

Table 3 Problem formulation parameters

Figure 4. The 3D trajectories of the ground moving target and the tracker UAV.

Figure 5. The 2D illustration of the ground moving target and the tracker UAV trajectories on the $\textit{x}\,-\,\textit{y}$ and $\textit{x}\,-\,\textit{z}$ planes.

Figure 6. The components of the position of the UAV and the target over time.

$I_{xx}=I_{yy}=0.0069{\rm kg.m}^2\quad I_{zz}=0.0137{\rm kg.m}^2\\$

$ \,J_r=0.0000551{\rm kg.m}^2 \, l_a=0.27{\rm m}\,$

Sinusoidal external disturbances on the six degrees of freedom are included. The magnitude of disturbances on the aerodynamic forces and moments have been set as in Ref. (Reference Alexis, Nikolakopoulos and Tzes24).

$\begin{array}{l}w_x=2u(t-70) \cos t ({\rm N/kg}) w_{\phi}=u(t-130) \sin t {\rm (N.m)} \\[5pt]w_y=2u(t-100) \cos t {\rm (N/kg)} w_{\theta}=u(t-170) \cos t {\rm (N.m)}\\[5pt]w_z=u(t-40) \sin t {\rm (N/kg)} w_{\psi}=u(t-200) \cos t \, {\rm (N.m)}\\[5pt]\end{array}$

Figures 4–11. show the tracking results in such a way that there is no error in estimating the target acceleration. This is not a general assumption and will be removed in the next simulation. This simulation evaluate the controller performance; however, the next simulation focuses on the guidance law, which should maintain the target in the FOV in the event of the controller errors. The 3D representation of GMT tracking is shown in Fig. 4, and its 2D representations are shown in Fig. 5. As can be seen, tracking and maintaining the moving target in the field of view has been done in various target manoeuver. Figure 6 shows the time history of the movement of the target and the UAV in x, y, and z directions. The tracking errors of the position and the attitude reference signals by the UAV are shown in Figs. 7 and 8, respectively. The jumps in the error signals are caused by sudden changes in the direction of the target trajectory. For example, at 25 and 35 seconds, the moving target changes its direction that causes jumps in the related error signals at the same time. Another reason for these error jumps are the wind turbulence. For example, in the $e_\psi$ error signal, perturbation occurs at the 200 seconds which causes error jumps at the same moment. It should be noted that despite the perturbation remaining, the controller has removed its effect well. These figures show that the proposed controller has good performance in tracking the reference generated trajectory by the guidance law.

Figure 7. The components of the position tracking error by the UAV.

Figure 8. The components of the attitude tracking errors by the UAV.

Figure 9. The control inputs of the quadrotor’s actuators.

Figure 10. The sufficient condition parameter values.

Figure 11. The parameters affecting the desired altitude calculation by the proposed guidance law.

Corresponding control inputs in this simulation are shown in Fig. 9. These signals indicate the PWMs of four actuators of the quadrotor. As can be seen, the amplitude of the inputs did not exceed the admissible interval $[0 \,1]$ , i.e. the actuators didn’t saturate. Figure 10 shows $\gamma=e_z+\epsilon $ is positive all over the simulation times. This means that Lemma 3 is satisfied during the tracking mission which guaranties the target preserving in the FOV. Figure 11 shows the parameters affecting the guidance law designed in accordance with Equation (20) in Lemma 3. It can be seen in this figure the relationship between the tracking errors on the $x-y$ plane, $z_{mt}$ , the proposed desired altitude and the results of the desired altitude tracking.

Due to some unpredictable reasons in the real world, such as target acceleration estimation error, quadrotor actuator error, and etc., the quadrotor position deviates from the target. Figure 12 shows the the x and y components of the relative distance between the quadrotor and the target caused by the presence of acceleration estimation errors as $\delta_{a_z}^{mt}=0.1{\rm m/s}^2$ , $\delta_{a_x}^{mt}=-0.5{\rm m/s}^2$ and $\delta_{a_y}^{mt}=-0.4 {\rm m/s}^2$ .

Figure 12. The x and y components of the relative distance between the quadrotor and the target.

In such circumstances, in order to examine the characteristics of the proposed strategy in relation to existing researches, a comparison is made between two quadrotors called $Quad_1$ and $Quad_2$ , in which the first quadrotor is equipped with the proposed guidance strategy and the second quadrotor is guided regardless of the altitude adjustment to hold the target.

For this purpose, $\gamma$ and $\lambda$ are defined as criterion parameters with the following form, which are defined according to the Lemmas 3 and 4, respectively;

\begin{equation*}\begin{array}{l} \gamma=e_z(k)+\epsilon\qquad {\rm with}\,\,z_r(k)=2 \max (|e_x|,|e_y|)+z_0 (k)\\[5pt] \lambda=e_z(k)+\epsilon-2\min (|e_x|,|e_y|) \qquad {\rm with}\,\, z_r(k)=z_0(k)\\ \end{array}\end{equation*}

According to Lemma 3, $\gamma(k)>0$ ensures the target retention in the FOV during the tracking mission by the $Quad_1$ and regarding Lemma 4, $\lambda(k)<0$ proves loss of target by the $Quad_2$ .

As it is illustrated in Fig. 13, $\gamma$ is positive during this simulation, which obtains the target retention in the FOV with $Qaud_1$ . On the other hand, $\lambda$ is negative, indicating the loss of target by $Quad_2$ .

Figure 13. The values of the criterion parameters to guarantee the target preserving by Quad1 and the target missing by the Quad2.

The desired altitude and their corresponding real values are illustrated in Fig. 14. As it can be seen, in contrast to $Quad_2$ , the altitude of $Quad_1$ increases to provide the FOV needed to ensure the target is maintained.

Figure 14. comparison between the reference altitudes of $\textit{Quad}_{\text{1}}$ and $\textit{Quad}_{\text{2}}$ and their tracking results.