Trajectory tracking of a quadrotor using model predictive and backstepping control with anti-disturbance compensation subjected to variations in the center of gravity

Abstract

The trajectory tracking control of a quadrotor unmanned aerial vehicle is done in this paper in the presence of changes in the center of gravity due to the variations of a connected load mass and its position, external disturbances, and parametric uncertainties. At first, dynamic equations are obtained using Newton–Euler relations, and then a two-loop control architecture is designed for the trajectory tracking task. Backstepping and model predictive control are used in the inner and outer loops, respectively. An extended disturbance observer is adopted to improve disturbance rejection capabilities in both loops. For the six-degrees-of-freedom unmanned aerial vehicle, the effects of time-varying mass, which leads to the change in the moment of inertia and center of gravity of the system, together with external disturbances and parametric uncertainties are taken into account. It is possible to mount an arbitrary number of different point masses at arbitrary positions on the unmanned aerial vehicle. The proposed control scheme is evaluated using simulation in the presence of uncertainties and external disturbances. According to the simulation results, the developed control scheme can achieve stable flight under different conditions.

Keywords

Quadrotor model predictive control backstepping trajectory tracking variable center of gravity extended disturbance observer

Introduction

In recent years, with progress of technology and science, the application of unmanned aerial vehicles (UAVs) has developed greatly (Zhao et al., 2021). Quadrotors, as one of the UAVs with vertical-takeoff-and-landing capability, have attracted significant interest in applications including search-and-rescue operations (Alonso-Mora et al., 2012), videography, evaluating agricultural fields, fire detection, mapping, delivery, and use in inaccessible environments. The control of these vehicles is complicated since they are coupled and underactuated (ul Amin et al., 2019).

In real-world scenarios, the weight of a UAV or its load may change over time. This can have negative impacts on the UAV. Variation in the system’s center of gravity (CG) is one of the challenges of aerial grasping. Hence, the control system must be able to cope with the dynamic changes in the CG of the system (Ariyanto et al., 2016). The deviation of the vehicle CG can be caused, for example, by releasing a load or by transporting equipment mounted outside the geometric center. In addition, the installation of photo/videography devices and robotic manipulators can also shift the CG with respect to the geometric center of the UAV.

Different control methods have been utilized for the quadrotors, from linear control methods, such as proportional-integral-derivative (PID) (Bouabdallah et al., 2004; Li and Li, 2011), linear quadratic regulator (Reyes-Valeria et al., 2013), to nonlinear control methods, such as sliding mode control (SMC) (Davoudi Dehkordi and Danesh, 2022), and its variant methods (Danesh et al., 2023; Liang and Yang, 2023), backstepping (BS) (Mohammadi et al., 2018), feedback linearization (Jebelli et al., 2018; Palunko and Fierro, 2011). In Li et al. (2020), a hierarchical approach was used for quadrotor trajectory tracking in the presence of disturbances and parametric uncertainties. A disturbance observer-based BS was employed to control the outer loop of an outer–inner loop structure, and an adaptive nonsingular terminal sliding mode controller was designed for the inner loop to deal with parametric uncertainties. The attitude control of a quadrotor was performed by Wang and Zhao (2022) by adopting adaptive control which uses a linear active disturbance rejection control (LADRC) in the presence of external disturbances. The adaptive control was used to tune the parameters of the LADRC.

In the research by Xuan-Mung and Hong (2019), robust BS controllers were employed for trajectory tracking using a quadrotor. The extended state observer (ESO) was utilized to estimate both disturbances and system state variables. A nonlinear disturbance observer was used in conjunction with a BS controller (Tripathi et al., 2016) for position, altitude, and attitude regulation controls. Simulation results show the effectiveness of the composite controller. The works by Chen et al. (2016), Dong et al. (2014), and Moeini et al. (2022) used disturbance observers along with BS controllers to enhance trajectory tracking performance against disturbances and uncertainty. In the work of Derakhshan et al. (2024), disturbance observers were designed for a coaxial octorotor in conjunction with the model predictive control (MPC) approach to attenuate the effects of disturbances and time-varying mass. Moreover, in the work of Thanh et al. (2022), a controller with integration of an extended state/disturbance observer with a nominal BS controller was formed to improve attitude and position control performance.

The disturbance problem caused by changes in CG has been addressed in different ways. The study accomplished by Zhang et al. (2014) to control an aircraft involves a situation in which the system’s mass changes abruptly. Using linear matrix inequality-based nonlinear adaptive robust control as a control method for quadrotors delivering packages with unknown masses in a windy environment was another work done by Kun and Hwang, (2016). A study was conducted by Sadeghzadeh et al. (2012) to investigate load release operations on a quadrotor. The UAV was controlled using gain-scheduled PID and MPC. In Wu et al. (2019), the attitude of a quadrotor was controlled in the presence of time-varying mass by SMC theory.

A quadrotor with variable mass was controlled by nonsingular fast terminal SMC based on an adaptive integral BS control strategy (Zhao et al., 2021). However, the effect of time-varying mass on the system moment of inertia was not considered. Ariyanto et al. (2016) analyzed the effects of a change in the UAV CG when a payload was connected to the UAV. PID controllers were utilized to control the altitude and attitude of the system.

MPC uses a cost function that can consider the constraints imposed on the system to produce the desired control inputs. It is possible to incorporate the future trajectory in the cost function for path tracking, which can improve path tracking performance and disturbance rejection. Additionally, it will lead to smoother control signals (Mayne et al., 2000). MPC is considered a relatively practical method for linear and slow systems such as industrial processes (Mayne, 2014). However, with advances in numerical capabilities, this method is also applicable to nonlinear processes such as aerospace and robotic systems (Zhou et al., 2019).

In the research by Ma et al., (2016), the use of active disturbance rejection and predictive control techniques was explored for trajectory tracking of a quadrotor UAV in the presence of disturbances. To estimate existing disturbances, the predictive controller utilizes an ESO. Active disturbance rejection control and MPC were used to control rotational and translational movements, respectively. Additionally, the study by Wang et al. (2020) used a predictive control method to stabilize quadrotor attitude dynamics with a payload on the UAV. In the study by Urbina-Brito et al. (2021), an MPC was proposed for a quadrotor with a hanging load, in which the controller considers the dynamics of the suspended load, and the model is expressed as an input-affine system to design the controller.

The review of literature shows that many works conducted on load transfer using a quadrotor have merely investigated its attitude control and the moment of inertia changes have not been taken into account. Furthermore, in the majority of the previous studies, the system was either linearized around certain operating points or small angles assumption was used. Therefore, controlling the overall dynamics of the system in the presence of external disturbances and uncertainties can be challenging.

Although MPC is known as a powerful control method, as shown in the simulations, this method is sensitive to uncertainties in the model and model mismatches. Moreover, BS control is not robust to disturbances. Hence, in this research, a robust approach is developed based on a combination of MPC, BS, and an extended disturbance observer (EDO) to improve the performance of the system in disturbance rejection. The main contributions of the present work can be summarized as follows:

An EDO is designed to estimate various disturbances across all degrees of freedom (DOF). These disturbances are compensated in advance to increase the robustness of the UAV. Accurate dynamic modeling is presented in six degrees of freedom (DOF). To bring the conditions closer to the real world, the effects of variation in mass, moment of inertia, and CG on the overall dynamics are considered along with disturbances and uncertainties. It is also possible to mount an arbitrary number of different point masses at arbitrary positions on the UAV. A two-loop control structure is proposed using MPC with final state cost and BS control methods enhanced by EDO without knowing the upper bound of disturbances. Small angle simplification is not used in the control design process and robust trajectory tracking is achieved in the presence of uncertainties, CG changes, and external disturbances.

The remainder of this paper is organized as follows: In Section “Dynamic model of quadrotor,” the dynamic model of the quadrotor system is derived considering mass and CG changes. In Section “Control design based on disturbance observer,” a disturbance observer is designed and the proposed control scheme including translational and rotational controllers is developed. Simulation results are given in Section “Simulation results” to evaluate the designed control scheme. Finally, conclusions are given in Section “Conclusion.”

Dynamic model of quadrotor

The structure of a quadrotor is shown in Figure 1, in which $Ω_{i}, i = 1, 2, \dots, 4$ illustrates the speed of the $i$ th UAV rotor. The mathematical model of the system is derived using the Newton–Euler approach. The inertial and body fixed frames which, respectively, denoted by I and B are considered. The inertial frame $(O_{I}, x_{I}, y_{I}, z_{I})$ is used to determine the linear position $P$ and the angular position $η$ . The body frame $(O_{B}, x_{B}, y_{B}, z_{B})$ is attached to the UAV and moves with it. The vector $P = {[x y z]}^{T}$ represents the position of the body frame measured in the inertial frame. The translational and rotational dynamics can be defined by

\frac{d}{d t} (m_{v} \dot{P}) = \frac{d}{d t} L = F_{B}

(1)

\frac{d}{dt} (I_{v} ω_{B}) = \frac{d}{dt} H = M_{B}

(2)

where $L$ and $F_{B}$ are linear momentum and external force, respectively. $H$ and $M_{B}$ are angular momentum and external torque, respectively.

Remark 1. $m_{v}$ contains the total mass of the quadrotor and the load or loads connected to the drone where $m_{v} = m_{c} + m_{a}$ . The fixed mass of the quadrotor is $m_{c}$ , and $m_{a}$ indicates the overall mass of the connected load or loads. $I_{v}$ is the total moment of inertia of the UAV. This moment of inertia consists of two fixed and variable parts as $I_{c}$ and $I_{a}$ , respectively, that is, $I_{v} = I_{c} + I_{a}$ .

Figure 1.

Body and inertial coordinate systems of the quadrotor with load.

The relation between linear velocities in the body frame $v_{B} = {[u v w]}^{T}$ and the inertial frame $\overset{\cdot}{P} = {[\overset{\cdot}{x} \overset{\cdot}{y} \overset{\cdot}{z}]}^{T}$ is

\overset{\cdot}{P} = R (η) v_{B}

(3)

Three consecutive rotations around the main axes can be performed to obtain the angular position vector $η = [ϕ θ ψ]^{T}$ , which contains the roll ( $ϕ$ ), pitch ( $θ$ ), and yaw ( $ψ$ ) Euler angles. The rotation matrix $R (η)$ is defined as

R (η) = R_{ψ} R_{θ} R_{ϕ} = [\begin{matrix} c_{ψ} c_{θ} & - s_{ψ} c_{ϕ} + c_{ψ} s_{θ} s_{ϕ} & s_{ψ} s_{ϕ} + c_{ψ} s_{θ} c_{ϕ} \\ s_{ψ} c_{θ} & c_{ψ} c_{ϕ} + s_{ψ} s_{θ} s_{ϕ} & - c_{ψ} s_{ϕ} + s_{ψ} s_{θ} c_{ϕ} \\ - s_{θ} & c_{θ} s_{ϕ} & c_{θ} c_{ϕ} \end{matrix}]

(4)

where $c_{k} = \cos (k), s_{k} = \sin (k), k = ϕ, θ, ψ$ . The relation between angular velocities in the body frame $ω_{B} = {[p q r]}^{T}$ and the inertial frame $\dot{η} = {[\dot{ϕ} \dot{θ} \dot{ψ}]}^{T}$ is expressed as

\overset{\cdot}{η} = W (η) ω_{B}

(5)

W (η) = [\begin{matrix} 1 & s_{ϕ} t_{θ} & c_{ϕ} t_{θ} \\ 0 & c_{ϕ} & - s_{ϕ} \\ 0 & s_{ϕ} / c_{θ} & c_{ϕ} / c_{θ} \end{matrix}]

(6)

where $W (η)$ is the translation matrix and $t_{θ} = \tan (θ)$ . The moments of inertia of the connected loads are considered as point masses on the UAV with respect to the body frame, and they are defined as follows:

I_{a} = [\begin{matrix} I_{xx} & - I_{xy} & - I_{xz} \\ - I_{yx} & I_{yy} & - I_{yz} \\ - I_{zx} & - I_{zy} & I_{zz} \end{matrix}]

(7)

I_{xx} = \sum_{i = 1}^{N} m_{li} (y_{li}^{2} + z_{li}^{2}), I_{yy} = \sum_{i = 1}^{N} m_{li} (x_{li}^{2} + z_{li}^{2})

I_{zz} = \sum_{i = 1}^{N} m_{li} (x_{li}^{2} + y_{li}^{2}), I_{xy} = I_{yx} = \sum_{i = 1}^{N} m_{li} x_{li} y_{li}

I_{xz} = I_{zx} = \sum_{i = 1}^{N} m_{li} x_{li} z_{li}, I_{yz} = I_{zy} = \sum_{i = 1}^{N} m_{li} y_{li} z_{li}

(8)

where $x_{li}$ , $y_{li}$ , and $z_{li}$ are the position coordinates of the $i$ th load, $N$ represents the number of the point masses rigidly mounted on the UAV and $m_{li} (kg)$ is the $i$ th load mass. Invoking (1), we have

m_{v} \overset{\cdot\cdot}{P} = - {\overset{\cdot}{m}}_{v} \overset{\cdot}{P} + F_{B}

(9)

The resulting external force in the inertial frame is $F_{B} = f_{G} + f_{c} + f_{l}$ with $f_{G} = [0 0 - m_{c} g]^{T}$ and

f_{c} = [\begin{matrix} (c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ}) u_{1} \\ (c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ}) u_{1} \\ (c_{θ} c_{ϕ}) u_{1} \end{matrix}]

(10)

where $f_{G}$ and $f_{c}$ are the gravity force vector of the quadrotor and force vector including the forces generated by the rotors, respectively. $u_{1}$ is the control input that denotes the total thrust. The force vector created by connecting the masses to the UAV is $f_{l} = \sum_{i = 1}^{N} f_{li}$ where $f_{li} = [0 0 - m_{li} g]^{T}$ . The relation (2) can be rewritten as

I_{v} {\overset{\cdot}{ω}}_{B} = - {\overset{\cdot}{I}}_{v} ω_{B} - ω_{B} \times I_{v} ω_{B} + M_{B}

(11)

The total external torque is $M_{B} = M_{u} + M_{h}$ , where $M_{u} = {[\begin{matrix} u_{2} & u_{3} & u_{4} \end{matrix}]}^{T}$ is the control torque vector and $M_{h} = \sum_{i = 1}^{N} M_{hi}$ is the torque vector created due to the addition of the point masses with

M_{h i} = r_{l i} \times f_{d i} = r_{l i} \times R {(η)}^{- 1} f_{l i} = [\begin{matrix} m_{l i} g (z_{l i} c_{θ} s_{ϕ} - y_{l i} c_{ϕ} c_{θ}) \\ m_{l i} g (x_{l i} c_{ϕ} c_{θ} + z_{l i} s_{θ}) \\ - m_{l i} g (x_{l i} c_{θ} s_{ϕ} + y_{l i} s_{θ}) \end{matrix}]

(12)

where $u_{i}, i = 2, 3, 4$ are control inputs of the rotational dynamics. $r_{li} = {[x_{li} y_{li} z_{li}]}^{T}$ is the coordinate vector of the $i$ th connected mass with respect to the body frame, and $f_{d}$ is the gravity force vector in the body frame. According to Figure 1, control inputs are obtained as

{\begin{matrix} u_{1} = f_{1} + f_{2} + f_{3} + f_{4} \\ u_{2} = b (f_{3} - f_{1}) \\ u_{3} = b (f_{4} - f_{2}) \\ u_{4} = Q_{2} + Q_{4} - Q_{1} - Q_{3} \end{matrix}

(13)

where $f_{i}$ and $Q_{i}$ for $i = 1, \dots, 4$ are the lift force and reaction torque produced by the rotors, respectively, and are given as (14). The distance from the geometric center of the UAV to each rotor is denoted by $b$ .

\begin{matrix} f_{i} = C_{l} Ω_{i}^{2} \\ Q_{i} = C_{d} Ω_{i}^{2} \end{matrix}

(14)

where $C_{l}$ and $C_{d}$ are thrust and drag constants, respectively. Substituting (14) into (13) yields

[\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \\ u_{4} \end{matrix}] = [\begin{matrix} C_{l} & C_{l} & C_{l} & C_{l} \\ - b C_{l} & 0 & b C_{l} & 0 \\ 0 & - b C_{l} & 0 & b C_{l} \\ - C_{d} & C_{d} & - C_{d} & C_{d} \end{matrix}] [\begin{matrix} Ω_{1}^{2} \\ Ω_{2}^{2} \\ Ω_{3}^{2} \\ Ω_{4}^{2} \end{matrix}]

(15)

The translational and rotational dynamics can be rewritten in the following form:

\overset{\cdot\cdot}{P} = (f_{G} + f_{c} + f_{l}) / m_{v} - ({\overset{\cdot}{m}}_{v} \overset{\cdot}{P}) / m_{v}

(16)

{\overset{\cdot}{ω}}_{B} = {I_{v}}^{- 1} (- ω_{B} \times I_{v} ω_{B} + M_{u} + M_{h}) - {I_{v}}^{- 1} ({\overset{\cdot}{I}}_{v} ω_{B})

(17)

In Figure 1, a general connection of a load to the UAV is shown.

Remark 2. The connected load can create torque around the body frame axes of the UAV.

Rewriting (16) as a set of equations, (18) is obtained that considers the effects of time-varying mass and moment of inertia. Finally, the dynamic equations of the UAV can be given as

\begin{matrix} \overset{\cdot\cdot}{x} = (c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ}) \frac{u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{x}}{m_{v}} \\ \overset{\cdot\cdot}{y} = (c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ}) \frac{u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{y}}{m_{v}} \\ \overset{\cdot\cdot}{z} = (c_{ϕ} c_{θ}) \frac{u_{1}}{m_{v}} - g - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{z}}{m_{v}} \end{matrix}

(18)

\begin{matrix} {\overset{\cdot}{ω}}_{B} = {I_{v}}^{- 1} (- ω_{B} \times I_{v} ω_{B} + M_{u} + M_{h}) - {I_{v}}^{- 1} ({\overset{\cdot}{I}}_{v} ω_{B}) \end{matrix}

(19)

Taking (18) into account, we have

\begin{matrix} \overset{\cdot\cdot}{x} = (c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ}) \frac{u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{x}}{m_{v}} + d_{1} \\ \overset{\cdot\cdot}{y} = (c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ}) \frac{u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{y}}{m_{v}} + d_{2} \\ \overset{\cdot\cdot}{z} = (c_{ϕ} c_{θ}) \frac{u_{1}}{m_{v}} - g - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{z}}{m_{v}} + d_{3} \end{matrix}

(20)

where $d_{1}$ , $d_{2}$ , and $d_{3}$ represent external disturbances in x, y, and z channels, respectively. To obtain the desired form, new terms are added to (20), hence, it can be rewritten as

\begin{matrix} \overset{\cdot\cdot}{x} = (c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ}) \frac{u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{x}}{m_{v}} - (c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ}) \frac{u_{1}}{m_{c}} \\ + (c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ}) \frac{u_{1}}{m_{c}} + d_{1} \end{matrix}

\begin{matrix} \overset{\cdot\cdot}{y} = (c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ}) \frac{u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{y}}{m_{v}} - (c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ}) \frac{u_{1}}{m_{c}} \\ + (c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ}) \frac{u_{1}}{m_{c}} + d_{2} \end{matrix}

\overset{\cdot\cdot}{z} = (c_{ϕ} c_{θ}) \frac{u_{1}}{m_{v}} - g - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{z}}{m_{v}} - (c_{ϕ} c_{θ}) \frac{u_{1}}{m_{c}} + (c_{ϕ} c_{θ}) \frac{u_{1}}{m_{c}} + d_{3}

(21)

Define

[\begin{matrix} u_{x} \\ u_{y} \\ u_{z} \end{matrix}] = [\begin{matrix} c_{ϕ} s_{θ} c_{ψ} + s_{ϕ} s_{ψ} \\ c_{ϕ} s_{θ} s_{ψ} - s_{ϕ} c_{ψ} \\ c_{ϕ} c_{θ} \end{matrix}]

(22)

The terms that cause uncertainty can be summed up in the lumped disturbance term. Therefore, the translational dynamics can be rewritten as

\overset{\cdot\cdot}{x} = \frac{u_{x} u_{1}}{m_{c}} + d_{u 1}

(23)

\overset{\cdot\cdot}{y} = \frac{u_{y} u_{1}}{m_{c}} + d_{u 2}

(24)

\overset{\cdot\cdot}{z} = \frac{u_{z} u_{1}}{m_{c}} - g + d_{u 3}

(25)

where

\begin{matrix} d_{u 1} = \frac{u_{x} u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{x}}{m_{v}} - \frac{u_{x} u_{1}}{m_{c}} + d_{1} \\ d_{u 2} = \frac{u_{y} u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{y}}{m_{v}} - \frac{u_{y} u_{1}}{m_{c}} + d_{2} \\ d_{u 3} = \frac{u_{z} u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{z}}{m_{v}} - \frac{u_{z} u_{1}}{m_{c}} + d_{3} \end{matrix}

(26)

New disturbance term $p_{t} = {[p_{1} p_{2} p_{3}]}^{T}$ is introduced to describe the remaining disturbances such as parametric uncertainties and the unmodeled dynamics of the system. Finally, the disturbances in the translational dynamics can be summarized as

\begin{matrix} d_{ζ 1} = \frac{u_{x} u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{x}}{m_{v}} - \frac{u_{x} u_{1}}{m_{c}} + p_{1} + d_{1} \\ d_{ζ 2} = \frac{u_{y} u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{y}}{m_{v}} - \frac{u_{y} u_{1}}{m_{c}} + p_{2} + d_{2} \\ d_{ζ 3} = \frac{u_{z} u_{1}}{m_{v}} - \frac{{\overset{\cdot}{m}}_{v} \overset{\cdot}{z}}{m_{v}} - \frac{u_{z} u_{1}}{m_{c}} + p_{3} + d_{3} \end{matrix}

(27)

These lumped disturbances are expressed in a compact form as $d_{t, l} = d_{u} + p_{t}$ . From (19), the rotational dynamics are rewritten as

\begin{matrix} {\overset{\cdot}{ω}}_{B} = {(I_{v})}^{- 1} (- ω_{B} \times I_{v} ω_{B} + M_{h}) - {(I_{v})}^{- 1} ({\overset{\cdot}{I}}_{v} ω_{B}) \\ + {(I_{v})}^{- 1} M_{u} + d_{λ} + p_{r} \end{matrix}

(28)

where $d_{λ}$ contains the external disturbances $d_{4}$ , $d_{5}$ , and $d_{6}$ in the attitude channels. The remaining uncertainties such as the neglected gyro effects, uncertainties created by the aerodynamic coefficients, and unmodeled dynamics are expressed by $p_{r} = {[p_{4} p_{5} p_{6}]}^{T}$ . The above relation can be summarized as

{\overset{\cdot}{ω}}_{B} = {[\overset{\cdot}{p} \overset{\cdot}{q} \overset{\cdot}{r}]}^{T} = {[g_{s 1} g_{s 2} g_{s 3}]}^{T}

(29)

where $g_{si} (.)$ for $i = 1, 2, 3$ includes $i$ th component of all defined terms of (28). The angular accelerations in (28) have been defined in the body frame; however, the control inputs are determined using the angular accelerations in the inertial frame. Therefore, the angular accelerations are calculated in the inertial frame. The first time derivative of (5) is

\overset{\cdot\cdot}{η} = \overset{\cdot}{W} (η) ω_{B} + W (η) {\overset{\cdot}{ω}}_{B}

(30)

where it can be expressed as

\begin{matrix} \overset{\cdot\cdot}{ϕ} = \overset{\cdot}{p} + \overset{\cdot}{q} s_{ϕ} t_{θ} + \overset{\cdot}{ϕ} q c_{ϕ} t_{θ} + \overset{\cdot}{θ} q (1 + t_{θ}^{2}) s_{ϕ} + \overset{\cdot}{r} c_{ϕ} t_{θ} - r \overset{\cdot}{ϕ} s_{ϕ} t_{θ} \\ + \overset{\cdot}{θ} r (1 + t_{θ}^{2}) c_{ϕ} \\ \overset{\cdot\cdot}{θ} = \overset{\cdot}{q} c_{ϕ} - \overset{\cdot}{ϕ} q s_{ϕ} - \overset{\cdot}{r} s_{ϕ} - \overset{\cdot}{ϕ} r c_{ϕ} \\ \overset{\cdot\cdot}{ψ} = \overset{\cdot}{q} s_{ϕ} s q_{θ} + \overset{\cdot}{ϕ} q c_{ϕ} s q_{θ} + \overset{\cdot}{θ} q s_{ϕ} t_{θ} s q_{θ} + \overset{\cdot}{r} c_{ϕ} s q_{θ} \\ - r \overset{\cdot}{ϕ} s_{ϕ} s q_{θ} + \overset{\cdot}{θ} r c_{ϕ} s q_{θ} t_{θ} \end{matrix}

(31)

where $s q_{θ}$ stands for $\sec (θ)$ . Substituting (29) into (31), the rotational relations can be written as

\begin{matrix} \overset{\cdot\cdot}{ϕ} = u_{2} / I_{x} + d_{ζ 4} \\ \overset{\cdot\cdot}{θ} = u_{3} / I_{y} + d_{ζ 5} \\ \overset{\cdot\cdot}{ψ} = u_{4} / I_{z} + d_{ζ 6} \end{matrix}

(32)

where

\begin{matrix} d_{ζ 4} = g_{s 1} + (r c_{ϕ} + q s_{ϕ}) \overset{\cdot}{θ} (1 + t_{θ}^{2}) \\ + [g_{s 2} s_{ϕ} + g_{s 3} c_{ϕ} + \overset{\cdot}{ϕ} (q c_{ϕ} - r s_{ϕ})] t_{θ} - u_{2} / I_{x} \\ d_{ζ 5} = g_{s 2} c_{ϕ} - \overset{\cdot}{ϕ} (q s_{ϕ} + r c_{ϕ}) - g_{s 3} s_{ϕ} - u_{3} / I_{y} \\ d_{ζ 6} = [g_{s 2} s_{ϕ} + g_{s 3} c_{ϕ} + \overset{\cdot}{ϕ} (q c_{ϕ} - r s_{ϕ}) + \overset{\cdot}{θ} t_{θ} (q s_{ϕ} + r c_{ϕ})] \\ s q_{θ} - u_{4} / I_{z} \end{matrix}

(33)

Control design based on disturbance observer

A disturbance estimation approach is employed to improve the performance of the designed controllers. After estimating the overall system disturbances, the control scheme uses them to make the system robust against lumped disturbances.

Disturbance observer design

Consider a second-order system defined as

\begin{matrix} {\overset{\cdot}{x}}_{1} = x_{2} \\ {\overset{\cdot}{x}}_{2} = f (x) + g (x) u + d_{w} (x, u, t) \end{matrix}

(34)

where $x_{1}, x_{2} \in R$ are the state variables and $u \in R$ is the control input, $d_{w} (x, u, t)$ is the lumped disturbance. It is assumed that $d_{w}$ is continuous and

| \frac{d^{j} d_{w} (t)}{d t^{j}} | \leq β for j = 0, 1, 2, \dots, r

(35)

where $β$ is a positive number and does not necessary to be known, and $r$ is the order of the disturbance observer.

In this study, a second-order EDO (Ginoya et al., 2013) is used. Using (23) and (34), the disturbance observer relations for position $x$ of the system can be written as

{\hat{d}}_{ζ 1} = p_{1 a} + M_{1 a} \overset{\cdot}{x}

(36)

{\overset{\cdot}{p}}_{1 a} = - M_{1 a} (u_{x} u_{1} / m_{c} + {\hat{d}}_{ζ 1}) + {\hat{\overset{\cdot}{d}}}_{ζ 1}

(37)

{\hat{\overset{\cdot}{d}}}_{ζ 1} = p_{1 b} + M_{1 b} \overset{\cdot}{x}

(38)

{\overset{\cdot}{p}}_{1 b} = - M_{1 b} (u_{x} u_{1} / m_{c} + {\hat{d}}_{ζ 1})

(39)

where ${\hat{d}}_{ζ 1}$ and ${\hat{\overset{\cdot}{d}}}_{ζ 1}$ are the estimates of $d_{ζ 1}$ and ${\overset{\cdot}{d}}_{ζ 1}$ , respectively, $p_{1 a}$ and $p_{1 b}$ are auxiliary variables. The constants $M_{1 a}$ , $M_{1 b}$ are observer gains. The disturbance observer for position $y$ is similar to position $x$ , so that in (39), $x$ is replaced by $y$ as well as the disturbance, auxiliary variables, and adjustment coefficients notation is changed to ${\hat{d}}_{ζ 2}$ , $p_{2 a}$ , $p_{2 b}$ , $M_{2 a}$ , and $M_{2 b}$ , respectively.

Similarly, the lumped disturbance estimation for position $z$ using (25) and (34) is given as

\begin{matrix} {\hat{d}}_{ζ 3} = p_{3 a} + M_{3 a} \overset{\cdot}{z} \\ {\overset{\cdot}{p}}_{3 a} = - M_{3 a} (u_{z} u_{1} / m_{c} - g + {\hat{d}}_{ζ 3}) + {\hat{\overset{\cdot}{d}}}_{ζ 3} \\ {\hat{\overset{\cdot}{d}}}_{ζ 3} = p_{3 b} + M_{3 b} \overset{\cdot}{z} \\ {\overset{\cdot}{p}}_{3 b} = - M_{3 b} (u_{z} u_{1} / m_{c} - g + {\hat{d}}_{ζ 3}) \end{matrix}

(40)

To estimate the lumped disturbances of the attitude motions, the subscripts $i = 4, 5, 6$ are used, which indicates the design of the disturbance observer for $ϕ$ , $θ$ , and $ψ$ channels, respectively. Hence, using (32) and (34) we have

\begin{matrix} {\hat{d}}_{ζ i} = p_{ia} + M_{ia} \overset{\cdot}{q} \\ {\overset{\cdot}{p}}_{ia} = - M_{ia} (u_{m} / τ + {\hat{d}}_{ζ i}) + {\hat{\overset{\cdot}{d}}}_{ζ i} \\ {\hat{\overset{\cdot}{d}}}_{ζ i} = p_{ib} + M_{ib} \overset{\cdot}{q} \\ {\overset{\cdot}{p}}_{ib} = - M_{ib} (u_{m} / τ + {\hat{d}}_{ζ i}) \end{matrix}

(41)

where for $i = 4$ the defined parameters in the observer are $u_{m} = u_{2}$ , $τ = I_{x}$ , and $q = ϕ$ . Similarly, for $i = 5$ the related parameters are $u_{m} = u_{3}$ , $τ = I_{y}$ , and $q = θ$ . Finally, for $i = 6$ the parameters are $u_{m} = u_{4}$ , $τ = I_{z}$ , and $q = ψ$ .

Stability analysis of the disturbance observer

Defining disturbance and disturbance derivative estimation errors as ${\tilde{d}}_{ζ 1} = d_{ζ 1} - {\hat{d}}_{ζ 1}$ and ${\tilde{\overset{\cdot}{d}}}_{ζ 1} = {\overset{\cdot}{d}}_{ζ 1} - {\hat{\overset{\cdot}{d}}}_{ζ 1}$ , respectively, the error vector is expressed as

\tilde{e} = [{\tilde{d}}_{ζ 1} {\tilde{\overset{\cdot}{d}}}_{ζ 1}]^{T}

(42)

From (23), (36), and (37), one can get

{\overset{\cdot}{\hat{d}}}_{ζ 1} = M_{1 a} {\tilde{d}}_{ζ 1} + {\hat{\overset{\cdot}{d}}}_{ζ 1}

(43)

Subtracting both sides of the recent equation from ${\overset{\cdot}{d}}_{ζ 1}$ yields

{\overset{\cdot}{\tilde{d}}}_{ζ 1} = - M_{1 a} {\tilde{d}}_{ζ 1} + {\overset{\cdot}{d}}_{ζ 1} - {\hat{\overset{\cdot}{d}}}_{ζ 1} = - M_{1 a} {\tilde{d}}_{ζ 1} + {\tilde{\overset{\cdot}{d}}}_{ζ 1}

(44)

Similarly, from (23), (38), and (39), the following relation can be obtained:

{\dot{\tilde{\dot{d}}}}_{ζ 1} = - M_{1 b} {\tilde{d}}_{ζ 1} + {\overset{\cdot\cdot}{d}}_{ζ 1}

(45)

Differentiating (44) and using (45) leads to

{\overset{\cdot\cdot}{\tilde{d}}}_{ζ 1} = - M_{1 a} {\overset{\cdot}{\tilde{d}}}_{ζ 1} - M_{1 b} {\tilde{d}}_{ζ 1} + {\overset{\cdot\cdot}{d}}_{ζ 1}

(46)

By recalling (44) and (45), finally, the dynamics of the estimation error can be written as

\overset{\cdot}{\tilde{e}} = D_{1} \tilde{e} + S_{1} {\overset{\cdot\cdot}{d}}_{ζ 1}

(47)

where

D_{1} = [\begin{matrix} - M_{1 a} & 1 \\ - M_{1 b} & 0 \end{matrix}], S_{1} = [\begin{matrix} 0 \\ 1 \end{matrix}]

(48)

According to (47) and (48), by choosing $M_{1 a}$ and $M_{1 b}$ such that the eigenvalues of $D_{1}$ are located on the left-hand side of the imaginary axis, a positive definite symmetric matrix $P_{1}$ (Ginoya et al., 2013) can be selected as

D_{1}^{T} P_{1} + P_{1} D_{1} = - M_{1}

(49)

For any positive definite symmetric matrix $M_{1}$ , assuming $λ_{\min}$ and $λ_{\max}$ as the smallest and greatest eigenvalues of $M_{1}$ , respectively, then

λ_{\min} ‖ \tilde{e} ‖^{2} \leq {\tilde{e}}^{T} M_{1} \tilde{e} \leq λ_{\max} ‖ \tilde{e} ‖^{2}

(50)

Defining the Lyapunov function candidate as

V (\tilde{e}) = {\tilde{e}}^{T} P_{1} \tilde{e}

(51)

Differentiating $V$ and using (47) and (35) yields

\begin{matrix} \overset{\cdot}{V} (\tilde{e}) = {\tilde{e}}^{T} (D_{1}^{T} P_{1} + P_{1} D_{1}) \tilde{e} + 2 {\tilde{e}}^{T} P_{1} S_{1} {\overset{\cdot\cdot}{d}}_{ζ 1} \\ \leq - {\tilde{e}}^{T} M_{1} \tilde{e} + 2 ‖ P_{1} S_{1} ‖ ‖ \tilde{e} ‖ β \\ \leq - λ_{\min} {‖ \tilde{e} ‖}^{2} + 2 ‖ P_{1} S_{1} ‖ ‖ \tilde{e} ‖ β \\ \leq - ‖ \tilde{e} ‖ (λ_{\min} ‖ \tilde{e} ‖ - 2 ‖ P_{1} S_{1} ‖ β) \end{matrix}

(52)

Finally, the estimation error norm is bounded as

‖ \tilde{e} ‖ \leq \frac{2 ‖ P_{1} S_{1} ‖ β}{λ_{\min}}

(53)

Remark 3. It should be noted that based on (35), the disturbance and its two first derivatives are required to be bounded; however, it is not necessary to know the bounds. Moreover, as seen, no special condition is assumed for the presented disturbance observer regarding the form of disturbances.

Controller design

To stabilize the UAV, a two-loop control structure is presented in which the system attitude dynamics are controlled in the inner loop, and the outer loop controls the system position variables. Figure 2 illustrates the trajectory tracking control structure.

Figure 2.

Block diagram of the proposed control structure for quadrotor.

Translational dynamics control

As a first step in designing the MPC of the translational subsystem, an internal model for prediction is derived. Defining

x_{t} (t) = {[\begin{matrix} x & \overset{\cdot}{x} & y & \overset{\cdot}{y} & z & \overset{\cdot}{z} \end{matrix}]}^{T} \in R^{6}

u_{c} (t) = {[\begin{matrix} u_{x} u_{1} & u_{y} u_{1} & u_{z} u_{1} - m_{c} g \end{matrix}]}^{T} = [v_{1} v_{2} v_{3}]^{T} \in R^{3}

d_{t, l} (t) = {[\begin{matrix} d_{ζ 1} & d_{ζ 2} & d_{ζ 3} \end{matrix}]}^{T} \in R^{3}

(54)

where $x_{t} (t)$ is the state vector, $u_{c} (t)$ is the control input vector, and $d_{t, l} (t)$ is the lumped disturbance vector. The following state space representation of relations (23)–(25) is given for the translational subsystem in the presence of all existing disturbances.

{\overset{\cdot}{x}}_{t} (t) = A_{c} x_{t} (t) + B_{c} u_{c} (t) + M_{c} d_{t, l} (t)

y_{t} (t) = C_{c} x_{t} (t)

(55)

where the matrices $A_{c} \in R^{6 \times 6}$ , $B_{c} \in R^{6 \times 3}$ , $C_{c} \in R^{6 \times 6}$ , and $M_{c} \in R^{6 \times 3}$ become

\begin{matrix} A_{c} = [\begin{matrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}], B_{c} = [\begin{matrix} 0 & 0 & 0 \\ 1 / m_{c} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 / m_{c} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 / m_{c} \end{matrix}], \\ M_{c} = [\begin{matrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{matrix}], C_{c} = I_{6 \times 6} \end{matrix}

(56)

The lumped disturbances are not initially considered in the controller design process of the translational movements; however, later feedforward compensation is applied to counteract them. Therefore, (55) yields

{\overset{\cdot}{x}}_{t} (t) = A_{c} x_{t} (t) + B_{c} u_{t} (t)

y_{t} (t) = C_{c} x_{t} (t)

(57)

Using the feedforward Euler method to discretize these equations produces

x_{t} (k + 1) = A_{d} x_{t} (k) + B_{d} u_{t} (k)

y_{t} (k) = C_{d} x_{t} (k)

(58)

with

A_{d} = I_{6 \times 6} + Δ t_{s} A_{c}, B_{d} = Δ t_{s} B_{c}, C_{d} = C_{c}

(59)

which results in

\begin{matrix} A_{d} = [\begin{matrix} 1 & Δ t_{s} & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & Δ t_{s} & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & Δ t_{s} \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}], \\ B_{d} = [\begin{matrix} 0 & 0 & 0 \\ Δ t_{s} / m_{c} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & Δ t_{s} / m_{c} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & Δ t_{s} / m_{c} \end{matrix}] \end{matrix}

(60)

where $Δ t_{s}$ is discretization sampling time. Using the recursive model (58), it is possible to obtain the system states along the prediction horizon by employing the following relations:

\begin{matrix} x_{t} (k | k) = x_{t} (k) \\ x_{t} (k + 1 | k) = A_{d} x_{t} (k) + B_{d} u_{t} (k | k) \\ x_{t} (k + 2 | k) = A_{d} x_{t} (k + 1 | k) + B_{d} u_{t} (k + 1 | k) \\ = A_{d}^{2} x_{t} (k) + A_{d} B_{d} u_{t} (k | k) + B_{d} u_{t} (k + 1 | k) \end{matrix}

⋮

\begin{matrix} x_{t} (k + N_{p} | k) = A_{d}^{N_{p}} x_{t} (k) + A_{d}^{N_{p} - 1} B_{d} u_{t} (k | k) \\ + A_{d}^{N_{p} - 2} B_{d} u_{t} (k + 1 | k) + \dots + \\ A_{d}^{N_{p} - N_{u}} B_{d} u_{t} (k + N_{u} - 1 | k) \end{matrix}

(61)

The system outputs along the prediction horizon can be obtained as

y_{t} (k + 1 | k) = C_{d} x_{t} (k + 1 | k)

y_{t} (k + 2 | k) = C_{d} x_{t} (k + 2 | k)

⋮

y_{t} (k + N_{p} | k) = C_{d} x_{t} (k + N_{p} | k)

(62)

where $N_{p}$ and $N_{u}$ represent the prediction and control horizon, respectively, with $N_{u} \leq N_{p}$ .

Remark 4. A prediction horizon denotes the number of future samples for which MPC can make a prediction.

The predicted state vector ${\tilde{x}}_{t} \in ℛ^{n_{m} N_{p}}$ , the predicted output vector ${\tilde{y}}_{t} \in R^{n_{m} N_{p}},$ and the future control input vector ${\tilde{u}}_{t} (k) \in R^{n_{c} N_{u}}$ with $n_{m} = 6$ and $n_{c} = 3$ are considered as follows:

\begin{matrix} {\tilde{x}}_{t} (k) = [\begin{matrix} x_{t} (k + 1 | k) \\ x_{t} (k + 2 | k) \\ ⋮ \\ x_{t} (k + N_{p} | k) \end{matrix}], {\tilde{u}}_{t} (k) = [\begin{matrix} u_{t} (k | k) \\ u_{t} (k + 1 | k) \\ ⋮ \\ u_{t} (k + N_{u} - 1 | k) \end{matrix}], \\ {\tilde{y}}_{t} (k) = [\begin{matrix} y_{t} (k + 1 | k) \\ y_{t} (k + 2 | k) \\ ⋮ \\ y_{t} (k + N_{p} | k) \end{matrix}] \end{matrix}

(63)

where $u_{t} (k + i | k)$ , $x_{t} (k + i | k),$ and $y_{t} (k + i | k)$ are the inputs, system state variables and outputs obtained at time $k + i$ with the information up to time $k$ . $x_{t} (k | k)$ is the initial condition at the beginning of each prediction horizon. Using (61) and (62), the predicted states and outputs can be written as

{\tilde{x}}_{t} (k) = Γ x_{t} (k) + Λ {\tilde{u}}_{t} (k)

{\tilde{y}}_{t} (k) = S {\tilde{x}}_{t} (k)

(64)

The matrices $Γ \in R^{n_{m} N_{p} \times n_{m}}$ , $Λ \in R^{n_{m} N_{p} \times n_{c} N_{u}},$ and $S \in R^{n_{m} N_{p} \times n_{m} N_{p}}$ are expressed as

Γ = {[{(A_{d})}^{T} {(A_{d}^{2})}^{T} \dots {(A_{d}^{N_{u}})}^{T} \dots {(A_{d}^{N_{p}})}^{T}]}^{T}

Λ = [\begin{matrix} B_{d} & 0 & 0 & \dots & 0 \\ A_{d} B_{d} & B_{d} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ A_{d}^{N_{u} - 1} B_{d} & A_{d}^{N_{u} - 2} B_{d} & A_{d}^{N_{u} - 3} B_{d} & \dots & B_{d} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ A_{d}^{N_{p} - 1} B_{d} & A_{d}^{N_{p} - 2} B_{d} & A_{d}^{N_{p} - 3} B_{d} & \dots & A_{d}^{N_{p} - N_{u}} B_{d} \end{matrix}]

(65)

S = diag (C_{d} C_{d} \dots C_{d})

(66)

As all state variables are controlled, hence, $y_{t} (k) = x_{t} (k)$ . It is assumed that ${\tilde{u}}_{t} (k + ρ | k) = 0$ for $ρ \geq N_{u}$ .

The next step is to present a cost function that produces the control inputs at each time step using optimization. The reference trajectory vector at time step $k$ is selected as

r_{w}^{T} (k) = [x_{d} (k) u_{x, d} (k) y_{d} (k) u_{y, d} (k) z_{d} (k) u_{z, d} (k)]

where $u_{x, d} (k), u_{y, d} (k),$ and $u_{z, d} (k)$ are the desired velocities in $x$ , $y$ , and $z$ channels, respectively. The optimization problem with the quadratic objective function is defined as

\begin{matrix} \min \\ {\tilde{u}}_{t} \end{matrix} J (k) = \sum_{i = 1}^{N_{p} - 1} ‖ r_{w} (k + i | k) - x_{t} (k + i | k) ‖_{Q}^{2} + \sum_{i = 0}^{N_{u}} ‖ u_{t} (k + i | k) ‖_{R}^{2} + (r_{w} (k + N_{p} | k) - x_{t} (k + N_{p} | k)) Q_{f} (r_{w} (k + N_{p} | k) - x_{t} (k + N_{p} | k))

(67)

The first part of the above objective function is a penalty on the controlled variables, the second part penalizes the increase in the control inputs, and the last part describes the tracking final cost of tracking variables. $Q \in R^{n_{m} \times n_{m}}$ and $R \in R^{n_{c} \times n_{c}}$ with $Q \geq 0$ , $R > 0$ are weighting matrices to impose design control objectives and control behavior regulation, respectively. $Q_{f} \in R^{n_{m} \times n_{m}}$ is a weighting matrix which indicates emphasis on the final cost and $Q_{f} \geq 0$ . Using (63), it can rewrite (67) as

\begin{matrix} \min \\ {\tilde{u}}_{t} \end{matrix} J (k) = (R_{s} (k) - {\tilde{x}}_{t} (k))^{T} \bar{Q} (R_{s} (k) - {\tilde{x}}_{t} (k)) + {\tilde{u}}_{t}^{T} \bar{R} {\tilde{u}}_{t}

(68)

where

\begin{matrix} \bar{Q} = diag (Q Q \dots Q Q_{f}) \in R^{n_{m} N_{p} \times n_{m} N_{p}} \\ \bar{R} = diag (R R \dots R) \in R^{n_{c} N_{u} \times n_{c} N_{u}} \end{matrix}

(69)

The reference signals are provided as follows:

R_{s} (k) = {[r_{w}^{T} (k + 1 | k) r_{w}^{T} (k + 2 | k) \dots r_{w}^{T} (k + N_{p} | k)]}^{T}

(70)

Substituting (64) into (68) yields

\begin{matrix} \min \\ {\tilde{u}}_{t} \end{matrix} J (k) = (R_{s} - Γ x_{t} (k) - Λ {\tilde{u}}_{t})^{T} \bar{Q} (R_{s} - Γ x_{t} (k) - Λ {\tilde{u}}_{t}) + {\tilde{u}}_{t}^{T} \bar{R} {\tilde{u}}_{t} = (R_{s} - Γ x_{t} (k))^{T} \bar{Q} (R_{s} - Γ x_{t} (k)) - 2 {\tilde{u}}_{t}^{T} Λ^{T} \bar{Q} (R_{s} - Γ x_{t} (k)) + {\tilde{u}}_{t}^{T} (Λ^{T} \bar{Q} Λ + \bar{R}) {\tilde{u}}_{t} = {\tilde{u}}_{t}^{T} (Λ^{T} \bar{Q} Λ + \bar{R}) {\tilde{u}}_{t} - 2 {\tilde{u}}_{t}^{T} Λ^{T} \bar{Q} (R_{s} - Γ x_{t} (k)) + Q_{a}

(71)

where $Q_{a} = (R_{s} - Γ x_{t} (k))^{T} \bar{Q} (R_{s} - Γ x_{t} (k))$ . Since the optimization is with respect to the future system inputs sequence, $Q_{a}$ has no effect on the optimal response, thus, it can be removed from the objective function. Applying the necessary condition of minimizing $J (k)$ , that is, $\partial J (k) / \partial {\tilde{u}}_{t} = 0$ , gives

\begin{matrix} \frac{\partial J (k)}{\partial {\tilde{u}}_{t}} = \frac{\partial}{\partial {\tilde{u}}_{t}} (- 2 {\tilde{u}}_{t}^{T} Λ^{T} \bar{Q} (R_{s} - Γ x_{t} (k)) \\ + {\tilde{u}}_{t}^{T} (Λ^{T} \bar{Q} Λ + \bar{R}) {\tilde{u}}_{t}) = 0 \end{matrix}

(72)

Finally, the optimal response is obtained as

{\tilde{u}}_{t} = (Λ^{T} \bar{Q} Λ + \bar{R})^{- 1} Λ^{T} \bar{Q} (R_{s} - Γ x_{t} (k))

(73)

Based on the receding horizon control law, only the first part of the optimal control sequence is used, while the rest is ignored. Therefore, the control law is rewritten as

u_{t} = V_{p} (Λ^{T} \bar{Q} Λ + \bar{R})^{- 1} Λ^{T} \bar{Q} (R_{s} - Γ x_{t} (k))

(74)

where

V_{p} = [I_{n_{c} \times n_{c}} 0 \dots 0]

In the following, the rotational-dynamic controller is designed, and the overall control strategy is developed.

Rotational-dynamics control

In this section, to control the rotational motion of the quadrotor, BS controllers are designed in the presence of uncertainties and disturbances. The rotational subsystem relationships in the presence of lumped disturbances can be written in the following state space form:

x_{r} (t) = {[\begin{matrix} ϕ & \overset{\cdot}{ϕ} & θ & \overset{\cdot}{θ} & ψ & \overset{\cdot}{ψ} \end{matrix}]}^{T} = {[\begin{matrix} x_{1} & x_{2} & x_{3} & x_{4} & x_{5} & x_{6} \end{matrix}]}^{T}

(75)

{\overset{\cdot}{x}}_{r} (t) = (\begin{matrix} x_{2} \\ a_{1} u_{2} + d_{ζ 4} \\ x_{4} \\ a_{2} u_{3} + d_{ζ 5} \\ x_{6} \\ a_{3} u_{4} + d_{ζ 6} \end{matrix})

(76)

where $a_{1} = 1 / I_{x}$ , $a_{2} = 1 / I_{y}$ , and $a_{3} = 1 / I_{z}$ . Considering the state variables $x_{1}$ and $x_{2}$ related to the roll movement as follows:

{\overset{\cdot}{x}}_{1} = x_{2}

{\overset{\cdot}{x}}_{2} = a_{1} u_{2} + d_{ζ 4}

(77)

Defining the tracking error as $e_{1} = x_{1} - x_{1 d}$ , with $x_{1 d}$ as the reference roll angle, the derivative of the error is written as

{\overset{\cdot}{e}}_{1} = {\overset{\cdot}{x}}_{1} - {\overset{\cdot}{x}}_{1 d} = x_{2} - {\overset{\cdot}{x}}_{1 d}

(78)

The Lyapunov function candidate is introduced as follows:

V_{1} = \frac{1}{2} e_{1}^{2}

(79)

The derivative of this function is

{\overset{\cdot}{V}}_{1} = e_{1} {\overset{\cdot}{e}}_{1} = e_{1} ({\overset{\cdot}{x}}_{1} - {\overset{\cdot}{x}}_{1 d}) = e_{1} (x_{2} - {\overset{\cdot}{x}}_{1 d})

(80)

Taking $x_{2}$ as a virtual control law in (80), its desired value is considered as $x_{2 c} = {\overset{\cdot}{x}}_{1 d} - c_{1} e_{1}$ , where $c_{1} > 0$ is a positive scalar that determines the speed of convergence. The second error variable is considered as follows:

e_{2} = x_{2} - {\overset{\cdot}{x}}_{1 d} + c_{1} e_{1}

(81)

Therefore, ${\overset{\cdot}{V}}_{1}$ can be rewritten as follows:

{\overset{\cdot}{V}}_{1} = e_{1} (e_{2} + {\overset{\cdot}{x}}_{1 d} - c_{1} e_{1} - {\overset{\cdot}{x}}_{1 d}) = e_{1} e_{2} - c_{1} e_{1}^{2}

(82)

To obtain the stability requirements, the next step is to consider the following augmented Lyapunov function:

V_{2} = V_{1} + \frac{1}{2} e_{2}^{2}

(83)

Taking the time derivative of $V_{2}$ yields

{\overset{\cdot}{V}}_{2} = {\overset{\cdot}{V}}_{1} + e_{2} {\overset{\cdot}{e}}_{2} = e_{1} e_{2} - c_{1} e_{1}^{2} + e_{2} ({\overset{\cdot}{x}}_{2} - {\overset{\cdot\cdot}{x}}_{1 d} + c_{1} {\overset{\cdot}{e}}_{1})

(84)

Substituting ${\overset{\cdot}{x}}_{2}$ model into (84) leads to

{\overset{\cdot}{V}}_{2} = e_{1} e_{2} - c_{1} e_{1}^{2} + e_{2} (a_{1} u_{2} + d_{ζ 4} - {\overset{\cdot\cdot}{x}}_{1 d} + c_{1} {\overset{\cdot}{e}}_{1})

(85)

The stability condition can be considered as follows:

{\overset{\cdot}{V}}_{2} = - c_{1} e_{1}^{2} - c_{2} e_{2}^{2} \leq 0

(86)

where $c_{2} > 0$ . Comparing (85) with (86), and using (78) and (81), finally, the control input to control the roll motion is determined as follows:

u_{2} = \frac{1}{a_{1}} ({\overset{\cdot\cdot}{x}}_{1 d} - e_{1} - d_{ζ 4} - c_{2} e_{2} - c_{1} {\overset{\cdot}{e}}_{1})

(87)

The pitch motion equations are given by

\begin{matrix} {\overset{\cdot}{x}}_{3} = x_{4} \\ {\overset{\cdot}{x}}_{4} = a_{2} u_{3} + d_{ζ 5} \end{matrix}

(88)

Similarly, the control input becomes

u_{3} = \frac{1}{a_{2}} ({\overset{\cdot\cdot}{x}}_{3 d} - e_{3} - d_{ζ 5} - c_{4} e_{4} - c_{3} {\overset{\cdot}{e}}_{3})

(89)

with

e_{3} = x_{3} - x_{3 d}, e_{4} = x_{4} - {\overset{\cdot}{x}}_{3 d} + c_{3} e_{3}

(90)

The yaw motion dynamics are given by

\begin{matrix} {\overset{\cdot}{x}}_{5} = x_{6} \\ {\overset{\cdot}{x}}_{6} = a_{3} u_{4} + d_{ζ 6} \end{matrix}

(91)

The control input is obtained as follows:

u_{4} = \frac{1}{a_{3}} ({\overset{\cdot\cdot}{x}}_{5 d} - e_{5} - d_{ζ 6} - c_{6} e_{6} - c_{5} {\overset{\cdot}{e}}_{5})

(92)

with

e_{5} = x_{5} - x_{5 d}, e_{6} = x_{6} - {\overset{\cdot}{x}}_{5 d} + c_{5} e_{5}

(93)

The values of $d_{ζ 4}$ , $d_{ζ 5},$ and $d_{ζ 6}$ in (87), (89), and (92) are not known, and the disturbance observer is used to estimate them. Regarding (54), the control input $u_{1}$ , the reference roll and pitch angles for the inner loop can be calculated as

\begin{matrix} u_{1} = \sqrt{v_{1}^{2} + v_{2}^{2} + {(v_{3} + m_{c} g)}^{2}} \\ x_{3 d} = θ_{d} = ta n^{- 1} (\frac{v_{1} \cos ψ_{d} + v_{2} \sin ψ_{d}}{v_{3} + m_{c} g}) \\ x_{1 d} = ϕ_{d} = ta n^{- 1} (\cos θ_{d} (\frac{v_{1} \sin ψ_{d} - v_{2} \cos ψ_{d}}{v_{3} + m_{c} g})) \end{matrix}

(94)

After obtaining the lumped disturbance estimation using the EDO, the control laws for the translational subsystem can be rewritten as

u_{c} = u_{t} - B_{p}^{- 1} {\hat{d}}_{t, l}

(95)

where $u_{t}$ is the translational MPC output and ${\hat{d}}_{t, l} = {[{\hat{d}}_{ζ 1} {\hat{d}}_{ζ 2} {\hat{d}}_{ζ 3}]}^{T}$ describes the estimated lumped disturbance vector of the translational subsystem, and $B_{p} = diag (1 / m_{c}, 1 / m_{c}, 1 / m_{c})$ . For the rotational subsystem, the proposed control laws are as follows:

\begin{matrix} u_{2} = \frac{1}{a_{1}} ({\overset{\cdot\cdot}{x}}_{1 d} - e_{1} - {\hat{d}}_{ζ 4} - c_{2} e_{2} - c_{1} {\overset{\cdot}{e}}_{1}) \\ u_{3} = \frac{1}{a_{2}} ({\overset{\cdot\cdot}{x}}_{3 d} - e_{3} - {\hat{d}}_{ζ 5} - c_{4} e_{4} - c_{3} {\overset{\cdot}{e}}_{3}) \\ u_{4} = \frac{1}{a_{3}} ({\overset{\cdot\cdot}{x}}_{5 d} - e_{5} - {\hat{d}}_{ζ 6} - c_{6} e_{6} - c_{5} {\overset{\cdot}{e}}_{5}) \end{matrix}

(96)

Simulation results

This section defines various scenarios to evaluate the performance of the presented controller. The effectiveness of the controller is tested with and without the disturbance rejection strategy. The designed controller is evaluated under conditions such as CG deviation, uncertainties, and external disturbances. The quadrotor parameters are given in Table 1. The parameters of the controllers were tuned as

N_{p} = 10, N_{u} = 3, Δ t_{s} = 0.1 s

Q_{t} = diag (30 30 30 1 1 1), Q_{f} = 2 Q_{t}, R_{t} = 0.2 I_{3 \times 3}

c_{1} = 10, c_{2} = 5, c_{3} = 10, c_{4} = 5, c_{5} = 4, c_{6} = 2

Table 1.

Parameters of the quadrotor system.

Name	Parameter	Value	Unit
Mass	$m_{c}$	0.74	kg
Moment of inertia matrix	$I_{c}$	Diag(4e−3, 4e−3, 84e−4)	kg/m²
Thrust factor	$C_{l}$	2.9e−5	N/s²
Drag factor	$C_{d}$	1.1e−6	Nm/s²
Length of the quadrotor arm	$b$	0.21	m
Gravity acceleration	$g$	9.81	m/s²

The tuned parameters for the EDO are

M_{1 a} = M_{2 a} = M_{3 a} = 6, M_{1 b} = M_{2 b} = M_{3 b} = 1

M_{4 a} = M_{5 a} = M_{6 a} = 8, M_{4 b} = M_{5 b} = M_{6 b} = 2

Two scenarios are defined to investigate the performance of the proposed approach. In the first scenario, the proposed controller is evaluated under CG deviation caused by connecting a point mass with a time-varying position. A similar situation can occur when external loads such as robotic manipulators, cameras, and packages are rigidly attached to UAVs. All adjusting parameters of the reduced form controller are identical to the proposed controller. The second scenario compares the proposed method with integral backstepping (IBS) and PID controllers. It is worth noting that for the IBS, a tracking differentiator and first order command filter were used (Rashad et al., 2015). The parameters of the PID controllers were tuned using the MATLAB/Simulink optimization toolbox. The regulated parameters in both scenarios are identical.

Simulation example 1

The reference trajectory is defined as

x_{d} = 1.5 \sin (0.35 t), y_{d} = 1.5 \cos (0.35 t), z_{d} = 0.5 t, ψ_{d} = 0

The initial conditions are as $x_{t} (0) = {[0 0 1 0 0 0]}^{T},$ $x_{r} (0) = {[0 0 0 0 0 0]}^{T}$ . The amount of point mass connected to the system is 0.3 kg and the load position varies as follows:

\begin{matrix} {\begin{matrix} x_{l} = - 0.01 t + 0.05 t \in [0, 10) \\ x_{l} = - 0.05 t \in [10, 60] \end{matrix} \\ {\begin{matrix} y_{l} = - 0.01 t + 0.05 t \in [0, 10) \\ y_{l} = - 0.05 t \in [10, 60] \end{matrix} \end{matrix}

z_{l} = - 0.08

The applied disturbances in the first scenario are

\begin{matrix} d_{1} (t) = 0.5 \cos (0.4 t) & t = 5 - 50 s \\ d_{2} (t) = 0.5 \sin (0.5 t) & t = 5 - 50 s \\ d_{3} (t) = 0.5 \cos (0.5 t) & t = 5 - 50 s \end{matrix}

(97)

\begin{matrix} d_{4} (t) = 0.5 \sin (0.3 t) & t = 20 - 50 s \\ d_{5} (t) = 0.5 \cos (0.8 t) & t = 20 - 50 s \\ d_{6} (t) = 0.5 \sin (0.4 t) & t = 20 - 50 s \end{matrix}

(98)

The simulation results in the first scenario are shown in Figures 3 to 5. The c1 case represents the case of CG changes, and the c2 case represents the condition of imposing external disturbances along with CG changes. The position tracking is shown in Figure 4(a). Clearly, the proposed approach outperforms MPC-BS without EDO. Generally, both methods follow the reference path and are stable. In addition, the performance of the proposed approach without considering CG variation and external disturbances (n.c.) can be observed. The tracking performance of BS controller in the inner loop is presented in Figure 4(b). It is observed that the attitude angles track their references generated by the proposed MPC and EDO in the outer loop. The references in the roll and pitch attitude angles curves are related to the case c2.

Figure 3.

2D and 3D tracking view.

Figure 4.

Trajectory tracking in the presence of external disturbances and CG changes: (a) position responses and (b) attitude responses.

Figure 5.

Lumped disturbances estimation.

The proposed method can successfully track the reference trajectories. The reason for this can be explained by the accuracy of the disturbance observer in estimating lumped disturbances and rejecting them actively. The performance of the EDO is shown in Figure 5. As seen, the observer estimates disturbances well. In Figure 6, variation in the control inputs is shown. It is observed that control inputs change to deal with load effects and external disturbances. At the beginning of the simulation in case (n.c.), the control inputs vary to move the UAV toward the desired reference. Once it reaches the desired reference path, the control inputs remain consistent. In the case of CG changes, the MPC-BS tracked path deviates from the reference (see Figure 4(a), red dashed curve), even though the system is still stable. According to Figures 3 and 4(a), the tracking performance of the MPC-BS controller is significantly better when the disturbance observer is used. Figure 7 shows variations in the connected mass position and moment of inertia. As a result of the variation in mass position, the moment of inertia is changed. Despite the CG changes, trajectory tracking using the designed controller is accurate in all directions.

Figure 6.

Control inputs.

Figure 7.

Mass position changes in $xy$ plane and moment of inertia variation.

Simulation example 2

The reference trajectory is given in the following form:

x_{d} = \sin (0.1 t), y_{d} = \cos (0.1 t), z_{d} = 4 - 3 \cos (0.1 t), ψ_{d} = 0

The initial conditions of the translational and rotational subsystems are $x_{t} (0) = {[0 0 0.5 0 0 0]}^{T}$ and $x_{r} (0) =$ ${[0 0 0 0 0.3 0]}^{T}$ , respectively. In this scenario, system behavior is evaluated in the presence of uncertainties and external disturbances. Additionally, the proposed approach is compared with the IBS and PID methods. The applied disturbances and their start and end times are defined as

\begin{matrix} d_{1}, d_{2}, d_{3} = 1 t = 10 - 60 s \\ d_{4}, d_{5}, d_{6} = 0.5 t = 0 - 90 s \end{matrix}

Unmodeled dynamics and remaining uncertainties are assumed to always exist in the system and are as

$\begin{matrix} p_{1} (t) = 0.2 \cos (0.8 t) & p_{4} (t) = 0.4 \cos (0.4 t) \\ p_{2} (t) = 0.3 \cos (0.6 t) & p_{5} (t) = 0.5 \sin (0.7 t + 3) \\ p_{3} (t) = 0.2 \sin (0.4 t + 3) & p_{6} (t) = 0.4 \cos (0.5 t + 4) \end{matrix}$

The simulation results of this scenario are depicted in Figures 8 and 9. The 3D and 2D spaces of trajectory tracking are shown in Figure 8, where it is observed that the proposed strategy achieves superior performance under disturbances and uncertainties. Figure 9(a) shows tracking performance related to the position variables. The PID method fails to achieve accurate trajectory tracking. This is due to the effects of uncertainties and disturbances and the PID controller has no capability to reject them. The IBS and proposed methods can track the reference path accurately. According to the results, the proposed method outperforms the IBS method and deviates less from the reference path. Figure 9(b) illustrates attitude tracking performance. It can be seen that the attitude controller can follow the references generated by the presented approach. The IBS method has oscillations at the beginning of the motion and also while applying external disturbances at the $ϕ$ and $θ$ angles. However, this method shows a better response in stabilizing the ψ angle in terms of integral square error value, despite exhibiting larger oscillations.

Figure 8.

2D and 3D tracking view.

Figure 9.

Tracking comparison between the proposed approach and other methods in the presence of uncertainties and external disturbances: (a) position responses and (b) attitude responses.

The proposed controller is capable of ensuring high performance attitude and position tracking for the quadrotor in the presence of uncertainties and disturbances. The results reveal that the proposed approach exhibits more accurate tracking performance compared with PID and IBS controllers under sustained and time-varying disturbances and uncertainties. According to the simulation results, the proposed approach can present a fast and robust response. Consequently, accurate trajectory tracking and robustness have been achieved against uncertainties, external disturbances, and CG changes.

Conclusion

This study presented an anti-disturbance control method that utilizes MPC and BS control combined with EDO to control a quadrotor with variation in the CG. A double-loop control structure was presented in which in the outer loop MPC stabilizes the position variables of the system, and in the inner loop BS was designed to control the attitude angles. The disturbances due to parametric uncertainties, unmodulated dynamics, external disturbances, and time-varying CG were considered across all DOF channels and estimated using the EDO. To increase the robustness of the UAV, these disturbances were compensated in advance. Moreover, it is possible to mount an arbitrary number of different point masses at arbitrary positions on the UAV. The simulations were carried out under time-varying disturbances and uncertainties. By comparing the presented approach with PID and IBS approaches, it can be seen that the presented approach achieves superior and smoother responses. From the simulation results, it can be concluded that the proposed approach offers significant robustness against different types of disturbances that the UAV may encounter in the real world.

Footnotes

Author contributions

All authors contributed to the study conception and design. Authors’ contributions are as follows: Material preparation: Reza Ebrahimpour Derakhshan. Writing the first draft of the manuscript: Reza Ebrahimpour Derakhshan. Review and editing: Mohammad Danesh and Reza Ebrahimpour Derakhshan. Supervision: Mohammad Danesh. Analysis and interpretation of results: Reza Ebrahimpour Derakhshan and Mohammad Danesh. All authors read and approved the final manuscript.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Mohammad Danesh

Data availability statement

The datasets generated or analyzed during the current study are not available publicly.

References

Alonso-Mora

Schoch

Breitenmoser

, et al. (2012) Object and animation display with multiple aerial vehicles. In: 2012 IEEE/RSJ international conference on intelligent robots and systems, Vilamoura-Algarve, Portugal, pp. 1078–1083.

Ariyanto

Munadi Paryanto , et al. (2016) Mathematical and physical modelling with dynamic change in the center of gravity of quadrotor. In: 2016 8th international conference on information technology and electrical engineering (ICITEE), Yogyakarta, Indonesia, pp. 1–6.

Bouabdallah

Noth

Siegwart

(2004) PID vs LQ control techniques applied to an indoor micro quadrotor. In: 2004 IEEE/RSJ international conference on intelligent robots and systems (IROS) (IEEE Cat. No. 04CH37566), Sendai, Japan, pp. 2451–2456.

Chen

Lei

Zhang

, et al. (2016) A novel nonlinear resilient control for a quadrotor UAV via backstepping control and nonlinear disturbance observer. Nonlinear Dynamics 85(2): 1281–1295.

Danesh

Jalalaei

Derakhshan

(2023) Auto-landing algorithm for quadrotor UAV using super-twisting second-order sliding mode control in the presence of external disturbances. International Journal of Dynamics and Control 11: 2940–2957.

Davoudi Dehkordi

Danesh

(2022) Positionable rotor quadrotor: Dynamic modeling and adaptive finite-time sliding-mode controller design. Journal of Guidance, Control, and Dynamics 45(3): 424–433.

Derakhshan

Danesh

Moosavi

(2024) Disturbance observer-based model predictive control of a coaxial octorotor with variable centre of gravity. IET Control Theory & Applications 18: n/a.

Dong

Zhu

, et al. (2014) High-performance trajectory tracking control of a quadrotor with disturbance observer. Sensors and Actuators A: Physical 211: 67–77.

Ginoya

Shendge

Phadke

(2013) Sliding mode control for mismatched uncertain systems using an extended disturbance observer. IEEE Transactions on Industrial Electronics 61(4): 1983–1992.

10.

Jebelli

Yagoub

MCE

Dhillon

(2018) Feedback linearization approach to fault tolerance for a micro quadrotor. In: 2018 IEEE international conference on industrial technology (ICIT), Lyon, France, pp. 165–168.

11.

Kun

Hwang

(2016) Linear matrix inequality-based nonlinear adaptive robust control of quadrotor. Journal of Guidance, Control, and Dynamics 39(5): 996–1008.

12.

(2011) Dynamic analysis and PID control for a quadrotor. In: 2011 IEEE international conference on mechatronics and automation, Beijing, China, pp. 573–578.

13.

(2020) Robust trajectory tracking control for a quadrotor subject to disturbances and model uncertainties. International Journal of Systems Science 51(5): 839–851.

14.

Liang

Yang

(2023) Adaptive nonsingular terminal sliding mode tracking control of quadrotor systems with mass uncertainties and input delay. Transactions of the Institute of Measurement and Control 45: 01423312231161140.

15.

Xia

, et al. (2016) Active disturbance rejection and predictive control strategy for a quadrotor helicopter. IET Control Theory & Applications 10(17): 2213–2222.

16.

Mayne

Rawlings

Rao

, et al. (2000) Constrained model predictive control: Stability and optimality. Automatica 36(6): 789–814.

17.

Mayne

(2014) Model predictive control: Recent developments and future promise. Automatica 50(12): 2967–2986.

18.

Moeini

Lynch

Zhao

(2022) A backstepping disturbance observer control for multirotor UAVs: Theory and experiment. International Journal of Control 95(9): 2364–2378.

19.

Mohammadi

Abbasi

Ghayour

, et al. (2018) Composite control design for trajectory tracking of a quadrotor in the presence of unknown disturbances. In: 2018 6th RSI international conference on robotics and mechatronics (IcRoM), Tehran, Iran, pp. 267–272.

20.

Palunko

Fierro

(2011) Adaptive control of a quadrotor with dynamic changes in the center of gravity. IFAC Proceedings Volumes 44(1): 2626–2631.

21.

Rashad

Aboudonia

El-Badawy

(2015) Backstepping trajectory tracking control of a quadrotor with disturbance rejection. In: 2015 XXV international conference on information, communication and automation technologies (ICAT), Sarajevo, Bosnia and Herzegovina, pp. 1–7.

22.

Reyes-Valeria

Enriquez-Caldera

Camacho-Lara

, et al. (2013) LQR control for a quadrotor using unit quaternions: Modeling and simulation. In: CONIELECOMP 2013, 23rd international conference on electronics, communications and computing, Cholula, Puebla, Mexico, pp. 172–178.

23.

Sadeghzadeh

Abdolhosseini

Zhang

Y.M.

(2012) Payload drop application of unmanned quadrotor helicopter using gain-scheduled PID and model predictive control techniques. In: International conference on intelligent robotics and applications, Berlin, Heidelberg, pp. 386–395.

24.

Thanh

HLNN

Huynh

, et al. (2022) Quadcopter UAVs extended states/disturbance observer-based nonlinear robust backstepping control. Sensors 22(14): 5082.

25.

Tripathi

Behera

Verma

(2016) Disturbance observer based backstepping controller for a quadcopter. In: IECON 2016—42nd annual conference of the IEEE Industrial Electronics Society, Florence, Italy pp. 108–113.

26.

ul Amin

Inayat

Jun

(2019) Finite time position and heading tracking control of coaxial octorotor based on extended inverse multi-quadratic radial basis function network and external disturbance observer. Journal of the Franklin Institute 356(8): 4240–4269.

27.

Urbina-Brito

Guerrero-Sánchez

M-S

Valencia-Palomo

, et al. (2021) A predictive control strategy for aerial payload transportation with an unmanned aerial vehicle. Mathematics 9(15): 1822.

28.

Wang

Cai

Zhang

, et al. (2020) Disturbance attenuation predictive optimal control for quad-rotor transporting unknown varying payload. IEEE Access 8: 44671–44686.

29.

Wang

Zhao

(2022) Adaptive-based linear active disturbance rejection attitude control for quadrotor with external disturbances. Transactions of the Institute of Measurement and Control 44(2): 286–298.

30.

Xiao

(2019) Modeling and sliding mode-based attitude tracking control of a quadrotor UAV with time-varying mass. ISA Transactions 124: 436–443.

31.

Xuan-Mung

Hong

(2019) Robust backstepping trajectory tracking control of a quadrotor with input saturation via extended state observer. Applied Sciences 9(23): 5184.

32.

Zhang

, et al. (2014) Safety modeling and simulation of multi-factor coupling heavy-equipment airdrop. Chinese Journal of Aeronautics 27(5): 1062–1069.

33.

Zhao

Ding

Jiang

, et al. (2021) A novel control strategy for quadrotors with variable mass and external disturbance. International Journal of Robust and Nonlinear Control 31(17): 8605–8631.

34.

Zhou

Sun

, et al. (2019) Position control of a tail-sitter UAV using successive linearization based model predictive control. Control Engineering Practice 91: 104125.