Enhanced trajectory tracking for quadrotors: disturbance observer state feedback control

TL;DR

This paper introduces a new control method for quadcopters that improves their ability to follow complex flight paths by reducing the impact of disturbances.

Contribution

A novel disturbance observer-based control algorithm is proposed to enhance quadcopter trajectory tracking in the presence of uncertainties.

Findings

The DOBC algorithm effectively compensates for disturbances in quadcopter dynamics.

Simulations show accurate tracking of both simple and complex 3D trajectories.

The method improves quadcopter performance under uncertainties and external disturbances.

Abstract

This study investigates the dynamics and trajectory tracking of quadcopters by utilizing the Disturbance Observer-based Control (DOBC) algorithm. The quadcopter’s dynamic model, which consists of six degrees of freedom, includes both disturbances and uncertainties in the model. The DOBC approach incorporates the disturbance model into the system by introducing it at the input variables. It then compensates for disturbances to achieve accurate tracking of different reference trajectories. The simulated trajectories span a range of motion, varying from simple straight paths to complex spiral paths. In order to verify and evaluate the efficacy of the suggested control technique, simulations are performed using MATLAB. The simulations conclusively show that the disturbance observer-based method effectively achieves the tracking of specified reference trajectories in three-dimensional space. The study highlights the effectiveness of the DOBC algorithm in reducing the effects of disturbances and uncertainties, thereby improving the quadcopter’s capacity to accurately track various trajectories.

Figures

Reference signal (brown) and following results (blue) with a steep slope track in the x-, y-, and z-axes.

Reference signal (brown) and following results (blue) with a cylindrical spiral track (A) in the x-, y-, and z-axes and (B) in 3D space.

Reference signal (brown) and following results (blue) with a conic track (A) in the x-, y-, and z-axes and (B) in 3D space.

Reference signal (brown) and following results (blue) based on the integral type D with conic track in the y-axis (A) and in 3D space (B).

Reference signal (brown) and following result (blue) with conic y-axis track presentation for the mass (A) and inertial (B) change conditions.

Full text

Introduction

The widespread use of quadcopters in various industries highlights their ability to be used in different applications, including important roles in military and law enforcement operations as well as essential functions in fire control departments (Jha, 2017; Sebbane, 2018; Castiglioni et al., 2017). Quadcopters are strategically deployed in various sectors and have a crucial role in revolutionizing agriculture by offering valuable monitoring capabilities. Furthermore, their proficiency in quickly capturing incidents or events has established them as essential instruments in the field of electronic media coverage (Sriram Reddy & Nippun Kumaar, 2021).

Manufacturers are producing a large number of quadcopters in response to the increasing demand for these aerial vehicles (Deepak & Singh, 2016). As production goes up, it becomes clearer how important it is to develop new algorithms that can accurately track references and make it easier for quadcopters to ignore disturbances (Mendoza-Soto, Corona-Sánchez & Rodríguez-Cortés, 2018; Song, Zhao & Theil, 2023; Garlick & Bradley, 2021). The need to increase the effectiveness of quadcopters in various fields is what drives the primary emphasis on improving algorithms in this constantly changing environment. The pursuit of resilient algorithms is in line with the primary objective of fully exploiting the capabilities of quadcopter technology and guaranteeing its smooth incorporation into various operational settings.

The distinctive features of rotorcrafts and helicopters, including their proficiency in vertical takeoff and landing combined with exceptional maneuverability, make them essential in both manned and unmanned aerial vehicle applications (García Carrillo et al., 2013). Significantly, in the field of unmanned aerial vehicles (UAVs), their forceful movements caused by a small moment of inertia increase the difficulty of determining their current state and require accurate implementation of controllers (Ji, Ma & Ge, 2020; Dalwadi, Deb & Ozana, 2023; Mahony, Kumar & Corke, 2012; Derafa, Madani & Benallegue, 2006).

The quadcopter, a sophisticated device that combines electronics, mechanics, and the fundamental principles of aviation, represents this intricate nature (Chen et al., 2016). Therefore, the development of control algorithms requires a comprehensive approach that takes into account multiple factors. It is important to prioritize the resolution of input and output disturbances (Ma et al., 2023) while ensuring the implementation of effective and efficient control actions (Cole & Wickenheiser, 2019; Wang & Guan, 2022; Xuan-Mung et al., 2022). This comprehensive strategy considers the complex interaction among the device’s mechanical structure, electronic systems, and the dynamic flight environment it operates in.

The design of an efficient control algorithm is crucial to successfully navigating the complex details of quadcopters, with the ultimate objective of maximizing their performance in various operational situations (Li, Ge & Lee, 2021; Zhang, Zhou & Li, 2023; Zhang et al., 2023). This approach guarantees both stability and precision during flight and highlights the ongoing development of control methods to meet the changing requirements of unmanned aerial vehicle technology (Shi et al., 2023; Chen et al., 2022; Lu et al., 2022).

Moreover, within the wide range of research studies focused on quadcopters, a crucial area of investigation is the analysis of how disturbances affect the state inputs of these aerial vehicles (Islam, Liu & El Saddik, 2015). This factor has a substantial impact on the simulation results in terms of their flight paths. This study investigates the utilization of a disturbance observer-based control method (Li et al., 2014; Chen et al., 2016; Wang & Guan, 2022; Liu, Chen & Shi, 2022) to effectively address this crucial factor. The quadcopter’s behavior is mathematically described by a set of nonlinear equations (Joos et al., 2017; Huang, 2022), which are sometimes approximated as linear equations when certain assumptions are made.

The system model in quadcopter position control accounts for disturbances that affect the state inputs, utilizing closed-loop feedback control. From the perspective of control system design, it is particularly effective to develop a state feedback control system that includes input disturbance estimation (Wang & Guan, 2022). A comparison of simulation results using various disturbance models demonstrates the superiority of simulation outcomes when the disturbance model is in line with the scenario. This methodology exhibits impressive tracking abilities that are well-suited for different trajectory needs, indicating significant potential in the field of quadcopter aviation.

Quadcopter Model

Quadcopter dynamics

In quadrotor control, the body frame allows for the definition of the total thrust and rotational torques asthe total thrust T and rotational torques in the body frame can be defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} & {\mathrm{T}}_{\mathrm{total}}=\sum _{i=1}^{4}{T}_{i}={b}_{t}\sum _{i=1}^{4}{\omega }_{i}^{2}, & {\tau }_{\phi }={d}_{mm}{b}_{t} \left( {\omega }_{4}^{2}-{\omega }_{2}^{2} \right) , & {\tau }_{\theta }={d}_{mm}{b}_{t} \left( {\omega }_{3}^{2}-{\omega }_{1}^{2} \right) , & {\tau }_{\psi }={k}_{d} \left( {\omega }_{1}^{2}+{\omega }_{3}^{2}-{\omega }_{2}^{2}-{\omega }_{4}^{2} \right) ,\end{eqnarray*}\end{document}

where kd is the drag constant, ωi is the angular speed of the ith rotor, dmm denotes the distance extending from the rotor to the center of mass of the quadcopter, bt represents the thrust constant, which is dependent on variables including air density, as well as the dimensions of the quadcopter propeller, specifically its length and radius.

The rotor torque τMotori__ around its axis of rotation is opposite to the aerodynamic drag, which leads to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{\tau }_{{M}_{i}}=b{\omega }_{i}^{2}+{I}_{M}{\dot {\omega }}_{i},\end{eqnarray*}\end{document}

where b is also the drag constant and IM is the inertia moment associated with the quadcopter’s rotor. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {\omega }\end{document} can be neglected under quasi-static state, so that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{\tau }_{{M}_{i}}=b{\omega }_{i}^{2}.\end{eqnarray*}\end{document}

For each rotor, the rotation axis also moves in the fixed frame of the fuselage, thus generating the gyroscopic torque

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{\tau }_{gyro}=\sum _{i=1}^{4}{I}_{M} \left( \Omega \times {U}_{z} \right) {\omega }_{i},\end{eqnarray*}\end{document}

where Uz is the unit vector along inertial z-axis and Ω is the angular velocity.

When transitioning from the quadcopter’s body coordinate system to the inertial coordinate system, a rotation matrix is employed below,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{\mathrm{R}}_{\mathrm{o}}= \left[ \begin{array}{@{}ccc@{}} \displaystyle {C}_{\theta }{C}_{\psi }&\displaystyle {S}_{\theta }{S}_{\phi }{C}_{\psi }-{C}_{\phi }{S}_{\psi }&\displaystyle {S}_{\theta }{C}_{\psi }{C}_{\phi }+{S}_{\phi }{S}_{\psi }\\ \displaystyle {C}_{\theta }{S}_{\psi }&\displaystyle {S}_{\theta }{S}_{\phi }{S}_{\psi }+{C}_{\phi }{C}_{\psi }&\displaystyle {S}_{\theta }{S}_{\psi }{S}_{\phi }-{S}_{\phi }{C}_{\psi }\\ \displaystyle -{S}_{\theta }&\displaystyle {S}_{\phi }{C}_{\theta }&\displaystyle {C}_{\phi }{C}_{\theta } \end{array} \right] , \end{eqnarray*}\end{document}

where Sx = sin(x) and Cx = cos(x), with x = θ or ψ. In addition, from the angular velocities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {\eta }=[\dot {\phi },\dot {\theta },\dot {\psi }]^{T}\end{document} of inertial frame to the rates of angular motion Ω, where Ω = [p, q, r]^T^, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\Omega =\mathbi{F}\dot {\eta },\end{eqnarray*}\end{document}

where the transformation matrix is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\mathbi{F}= \left[ \begin{array}{@{}ccc@{}} \displaystyle 1&\displaystyle 0&\displaystyle -\sin \nolimits \theta \\ \displaystyle 0&\displaystyle \cos \nolimits \phi &\displaystyle \sin \nolimits \phi \cos \nolimits \theta \\ \displaystyle 0&\displaystyle -\sin \nolimits \phi &\displaystyle \cos \nolimits \phi \cos \nolimits \theta \end{array} \right] .\end{eqnarray*}\end{document}

For the translational dynamics of the quadcopter establishing. According to Newton’s second law,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}m \left[ \begin{array}{@{}c@{}} \displaystyle \ddot{x}\\ \displaystyle \ddot{y}\\ \displaystyle \ddot{z} \end{array} \right] = \left[ \begin{array}{@{}c@{}} \displaystyle 0\\ \displaystyle 0\\ \displaystyle -mg \end{array} \right] +{\mathrm{R}}_{\mathrm{o}} \left[ \begin{array}{@{}c@{}} \displaystyle 0\\ \displaystyle 0\\ \displaystyle {\mathrm{T}}_{\mathrm{total}} \end{array} \right] .\end{eqnarray*}\end{document}

Expending the equation, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\ddot{x} & =- \frac{{k}_{x}}{m} \dot {x}+ \frac{{\mathrm{T}}_{\mathrm{total}}}{m} (\sin \nolimits \theta \cos \nolimits \psi \cos \nolimits \phi +\sin \nolimits \phi \sin \nolimits \psi ), & \ddot{y} & =- \frac{{k}_{y}}{m} \dot {y}+ \frac{{\mathrm{T}}_{\mathrm{total}}}{m} (\sin \nolimits \theta \sin \nolimits \psi \cos \nolimits \phi -\sin \nolimits \phi \cos \nolimits \psi ),\ddot{z} & =- \frac{{k}_{z}}{m} \dot {z}+ \frac{{\mathrm{T}}_{\mathrm{total}}}{m} \cos \nolimits \phi \cos \nolimits \theta -g,\end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \ddot{x}\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \ddot{y}\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \ddot{z}\end{document} represent the acceleration in their respective axes. kx, ky and kz denote the coefficients of air friction, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {x},\dot {y},and\dot {z}\end{document} are the velocities in the x −, y −, and z − axes, respectively.

In accordance with Euler’s equation of motion, the rotational dynamics of quadcopter is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \left[ \begin{array}{@{}c@{}} \displaystyle \dot {p}\\ \displaystyle \dot {q}\\ \displaystyle \dot {r} \end{array} \right] = \left[ \begin{array}{@{}ccc@{}} \displaystyle 1/{I}_{x}&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 1/{I}_{y}&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 1/{I}_{z} \end{array} \right] \left[ \begin{array}{@{}c@{}} \displaystyle {\tau }_{\phi }\\ \displaystyle {\tau }_{\theta }\\ \displaystyle {\tau }_{\psi } \end{array} \right] + \left[ \begin{array}{@{}l@{}} \displaystyle qr({I}_{y}-{I}_{z})/{I}_{x}\\ \displaystyle pr({I}_{z}-{I}_{x})/{I}_{y}\\ \displaystyle pq({I}_{x}-{I}_{y})/{I}_{z} \end{array} \right] . \end{eqnarray*}\end{document}

Simplify the above equation, which leads to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\dot {p}= \frac{{\tau }_{\phi }}{{I}_{x}} + \frac{rq({I}_{y}-{I}_{z})}{{I}_{x}} ,\dot {q}= \frac{{\tau }_{\theta }}{{I}_{y}} + \frac{rp({I}_{z}-{I}_{x})}{{I}_{y}} ,\dot {r}= \frac{{\tau }_{\psi }}{{I}_{z}} + \frac{qp({I}_{x}-{I}_{y})}{{I}_{z}} ,\end{eqnarray*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {p}\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {q}\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {r}\end{document} stand for the angular acceleration parameters within the quadcopter body-centric coordinate system. Ix, Iy and Iz represent the inertia tensors associated with the x −, y −, and z − axes, respectively.

Moreover, the relationship between the Euler angular velocities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {\phi }\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {\theta },\dot {\psi }\end{document} and the angular velocities (p, q, r) is expressed as follows

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} & \left[ \begin{array}{@{}c@{}} \displaystyle \dot {\phi }\\ \displaystyle \dot {\theta }\\ \displaystyle \dot {\psi } \end{array} \right] = \left[ \begin{array}{@{}ccc@{}} \displaystyle 1&\displaystyle {T}_{\theta }{S}_{\phi }&\displaystyle {T}_{\theta }{C}_{\phi }\\ \displaystyle 0&\displaystyle {C}_{\phi }&\displaystyle -{S}_{\phi }\\ \displaystyle 0&\displaystyle \frac{{S}_{\phi }}{{C}_{\theta }} &\displaystyle \frac{{C}_{\phi }}{{C}_{\theta }} \end{array} \right] \left[ \begin{array}{@{}l@{}} \displaystyle p\\ \displaystyle q\\ \displaystyle r \end{array} \right] .\end{eqnarray*}\end{document}

Expending the equation, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\dot {\phi }=p+q{T}_{\theta }{S}_{\phi }+r{T}_{\theta }{C}_{\phi }, & \dot {\theta }=q{C}_{\phi }-r{S}_{\phi },\dot {\psi }=q \frac{{S}_{\phi }}{{C}_{\theta }} +r \frac{{C}_{\phi }}{{C}_{\theta }} ,\end{eqnarray*}\end{document}

where Tx = tan(x), Sx = sin(x), and Cx = cos(x), with x = ϕ or θ; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {\phi },\dot {\theta }\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {\psi }\end{document} are the angular velocities within the inertial coordinate system along the x −, y −, and z − axes, respectively.

State space model

Define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {x}=a\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {y}=b\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \dot {z}=c\end{document} , the state space model contains 12 state variables for the highly coupled quadcopter system:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\dot {x}=a, & \dot {y}=b,\dot {z}=c,\dot {a}=- \frac{{k}_{x}}{m} a+ \frac{{\mathrm{T}}_{\mathrm{total}}}{m} (\sin \nolimits \theta \cos \nolimits \psi \cos \nolimits \phi +\sin \nolimits \phi \sin \nolimits \psi ),\dot {b}=- \frac{{k}_{y}}{m} b+ \frac{{\mathrm{T}}_{\mathrm{total}}}{m} (\sin \nolimits \theta \sin \nolimits \psi \cos \nolimits \phi -\sin \nolimits \phi \cos \nolimits \psi ),\dot {c}=- \frac{{k}_{z}}{m} c+ \frac{{\mathrm{T}}_{\mathrm{total}}}{m} \cos \nolimits \phi \cos \nolimits \theta -g,\dot {p}= \frac{{\tau }_{\phi }}{{I}_{x}} + \frac{rq({I}_{y}-{I}_{z})}{{I}_{x}} ,\dot {q}= \frac{{\tau }_{\theta }}{{I}_{y}} + \frac{rp({I}_{z}-{I}_{x})}{{I}_{y}} ,\dot {r}= \frac{{\tau }_{\psi }}{{I}_{z}} + \frac{qp({I}_{x}-{I}_{y})}{{I}_{z}} ,\dot {\phi }=p+q{T}_{\theta }{S}_{\phi }+r{T}_{\theta }{C}_{\phi },\dot {\theta }=q{C}_{\phi }-r{S}_{\phi },\dot {\psi }=q \frac{{S}_{\phi }}{{C}_{\theta }} +r \frac{{C}_{\phi }}{{C}_{\theta }} . \end{eqnarray*}\end{document}

From the above equations, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\dot {X}=f(X,U), \end{eqnarray*}\end{document}

where X = (x, y, z, a, b, c, p, q, r, ϕ, θ, ψ)^T^ represents the state vertor, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} U= \left( {\mathrm{T}}{\mathrm{total}},{\tau }{\phi },{\tau }{\theta },{\tau }{\psi } \right) \end{document} denotes the thrust vector acting as control input, and Y = (x, y, z, θ) is the output vector. To ensure precise trajectory tracking, the initial focus lies on confirming the alignment of the first three vectors concerning position control. Given the interconnected nature of the quadcopter, an additional check involves validating an angular velocity along the y −axis. This acts as an indirect measure to ascertain the accuracy of the quadcopter’s interlinked system.

Disturbance observer design

The disturbance observer (DOB) is an advanced control approach used in engineering to accurately estimate and effectively prevent unexpected interruptions that occur within control systems. The main purpose of the system is to measure disruptions that arise from several sources, including external environmental factors, oscillations within the system, or errors in the modeling process. These disruptions provide a potential danger to the performance of the control system and require proactive measures to be taken through the observer’s assessments.

The fundamental principle of DOB is typically based on the mathematical depiction of the system, encompassing its state variables and outputs. By conducting a comparison between the projections of the model and the actual behavior of the system, DOB is able to detect inconsistencies, which are frequently caused by disturbance signals. By evaluating these disruptions, DOB produces a corrective signal that is incorporated into the input of the control system. This effectively reduces the negative effects of the disturbances and maintains the desired performance and stability of the system.

Quadcopter input disturbance observer design

Initially, it presupposes an input disturbance denoted as ν(i) in the system, which subsequently results in a transformation of the state-space model.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{x}_{n}(i+1)={A}_{n}{x}_{n}(i)+{B}_{n}(u(i)+\nu (i)), & y(i)={C}_{n}{x}_{n}(i),\end{eqnarray*}\end{document}

where u(i) and ν(i) stand for the control inputs, y(i) represents the system responses, and xn(i) represents the state variable set.

The input disturbance ν(i) is conceptualized as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\nu (i)= \frac{{j}^{-1}\epsilon (i)}{D({j}^{-1})} ,\end{eqnarray*}\end{document}

where ϵ(i) is consistent with the dimension of vector ν(i).

Let the disturbance model D(j^−1^) be

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}D({j}^{-1})=1+{d}_{1}{j}^{-1}+{d}_{2}{j}^{-2}+{d}_{3}{j}^{-3}+\ldots +{d}_{\gamma }{j}^{-\gamma }.\end{eqnarray*}\end{document}

Next, use the form of difference equation to express the disturbance vector ν(i + 1):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\nu (i+1)=-{d}_{1}\nu (i)-{d}_{2}\nu (i-1)-\ldots -{d}_{\gamma }\nu (i-\gamma +1)+\epsilon (i),\end{eqnarray*}\end{document}

where ϵ(i) represents white noise with a mean value of zero.

Following our previous study (Wang & Guan, 2022), by introducing matrix h(i) of size nγ × 1 to represent the disturbance group and presenting it with a state space mode, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}h(i)={ \left[ \begin{array}{@{}llll@{}} \displaystyle {\nu }^{T}(i)&\displaystyle {\nu }^{T}(i-1)&\displaystyle \cdots &\displaystyle {\nu }^{T}(i-\gamma +1) \end{array} \right] }^{T}.\end{eqnarray*}\end{document}

Then, using the state space model below to obtain the disturbance

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}h(i+1)={A}_{dis}h(i)+{B}_{\epsilon }\epsilon (i),\nu (i)={C}_{\epsilon }h(i).\end{eqnarray*}\end{document}

For matrix Adis, it shows the form as below:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{A}_{dis}= \left[ \begin{array}{@{}ccccc@{}} \displaystyle -{d}_{1}I&\displaystyle -{d}_{2}I&\displaystyle \ldots &\displaystyle -{d}_{\gamma -1}I&\displaystyle -{d}_{\gamma }I\\ \displaystyle I&\displaystyle O&\displaystyle \ldots &\displaystyle \ldots &\displaystyle O\\ \displaystyle O&\displaystyle I&\displaystyle O&\displaystyle \ldots &\displaystyle \vdots \\ \displaystyle \vdots &\displaystyle \ddots &\displaystyle \ddots &\displaystyle \ddots &\displaystyle \vdots \\ \displaystyle O&\displaystyle \ldots &\displaystyle O&\displaystyle I&\displaystyle O \end{array} \right] , \end{eqnarray*}\end{document}

where Bϵ is a nγ × n matrix and Cϵ is a n × nγ matrix, they both have the n × n identity matrix and the rest vertical and horizontal elements of the matrix are zeros. The leftover rows and columns in both matrices consist of zero matrices. O is a n × n zero matrix and I is the identity matrices with same dimensions.

For the disturbance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} D \left( {j}^{-1} \right) \end{document} , we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}det \left( {j}^{-1}I-{A}_{d} \right) =D({j}^{-1})^{n},\end{eqnarray*}\end{document}

where D(j^−1^) is the disturbance model in backward shift form. Matrix Adis includes n groups of eigenvalues which related to the zeros of the D(j^−1^).

The incorporation of the disturbance vector ν(i) into the state space model serves to mitigate the impact of external interference on the state estimate matrix xn(i). The introduced augmented state matrix is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}x(i)={ \left[ \begin{array}{@{}ll@{}} \displaystyle {x}_{n}^{T}(i)&\displaystyle h(i)^{T} \end{array} \right] }^{T}. \end{eqnarray*}\end{document}

We rewrite the previous state space model as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{x}_{n}(i+1)={A}_{n}{x}_{n}(i)+{B}_{n}{C}_{\epsilon }h(i)+{B}_{n}u(i). \end{eqnarray*}\end{document}

Under this condition, an augmented system space model is derived:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \left[ \begin{array}{@{}c@{}} \displaystyle {x}_{n}(i+1)\\ \displaystyle h(i+1) \end{array} \right] & =\overbrace{ \left[ \begin{array}{@{}c@{}} \displaystyle {A}_{n} & {B}_{n}{C}_{\epsilon }\\ \displaystyle {O}_{1} & {A}_{dis} \end{array} \right] {}}^{\mathcal{A}} \left[ \begin{array}{@{}c@{}} \displaystyle {x}_{n}(i)\\ \displaystyle h(i) \end{array} \right] & & +\overbrace{ \left[ \begin{array}{@{}c@{}} \displaystyle {B}_{n}\\ \displaystyle {O}_{2} \end{array} \right] {}}^{B}u(i)+ \left[ \begin{array}{@{}l@{}} \displaystyle {O}_{3}\\ \displaystyle {B}_{\epsilon } \end{array} \right] \epsilon (i), \end{eqnarray*}\end{document}

where O1,2,3 are zero matrices to adapt to corresponding dimensions; the integrated A and B matrices are regarded as the augmented matrices of the system.

The output is computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}y(i)=\overbrace{ \left[ \begin{array}{@{}l@{}} \displaystyle {C}_{n} & {O}_{4} \end{array} \right] {}}^{C} \left[ \begin{array}{@{}c@{}} \displaystyle {x}_{n}(i)\\ \displaystyle h(i) \end{array} \right] , \end{eqnarray*}\end{document}

where C is the augmented output matrix, and O4 is a zero matrix like O1,2,3.

After establishing the augmented model, we will design the observer on this basis, A and C matrices will be used in the design of this step. First define an observer gain Kob, it is a matrix with the same dimension as A and C. To get the value of Kob, we use the discrete-time LQR algorithm. The selected observer gain is to stabilize the state matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \left( A-{K}_{ob}C \right) \end{document} in the closed-loop observer error system deviation dynamics. As a result, the augmented state matrix x(i) becomes amenable to estimation.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\hat {x}(i+1)=A\hat {x}(i)+Bu(i)+{K}_{ob}(y(i)-C\hat {x}(i)).\end{eqnarray*}\end{document}

Noting that in order to ensure that A and C are observable, it is first necessary to ensure that the An and Cn matrices are observable, and the system has no zeroes which related to D(j^−1^). Under this condition, the observer can operate effectively.

Disturbance observer-based discrete-time state feedback control

Before appling augmented system state representation, converting continuous model to discrete model. The state space model in continuous time can be formulated as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{\dot {x}}_{n}(k)={A}_{c}{x}_{n}(k)+{B}_{c}u(k),\end{eqnarray*}\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}y(k)={C}_{c}{x}_{n}(k).\end{eqnarray*}\end{document}

The discrete-time formulation of the state space model can be succinctly articulated as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{x}_{n}(i+1)={A}_{n}{x}_{n}(i)+{B}_{n}u(i),\end{eqnarray*}\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}y(i)={C}_{n}{x}_{n}(i).\end{eqnarray*}\end{document}

Based on the design of adding disturbance to the input and its compensation in the state estimation, the design of the state feedback control is undertaken through the subsequent stepwise procedure. It continues to adopt the state-space formulation detailed in the preceding section, following the design of model:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{x}_{n}(i+1)={A}_{n}{x}_{n}(i)+{B}_{n}(u(i)+\nu (i)), & y(i)={C}_{n}{x}_{n}(i). \end{eqnarray*}\end{document}

Noting that the input vectors are u(i) and ν(i). It assumes that the observational properties are satisfied for matrices An and Cn and that the controllability conditions are met for the pair (An, Bn). Moreover, the system avoid the zeros which related to the disturbance occurring at the D(j^−1^).

We define the intermediary control signal at this stage:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\tilde {u}(k)=u(k)+\mu (k). \end{eqnarray*}\end{document}

Then the original xn(i + 1) becomes

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{x}_{n}(i+1)={A}_{n}{x}_{n}(i)+{B}_{n}\tilde {u}(i).\end{eqnarray*}\end{document}

Following our previous study (Wang & Guan, 2022), the principle behind state feedback control using a disturbance observer involves designing the state feedback controller K specifically for the intermediate control signal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tilde {u}(k)\end{document} . This controller aims to handle the estimated variable to effectively compensate for the input disturbance ν(k).

It integrates (xn(i + 1), h(i + 1)) and y(i) as a complete augmented model

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \left[ \begin{array}{@{}c@{}} \displaystyle {x}_{n}(i+1)\\ \displaystyle h(i+1) \end{array} \right] & =\overbrace{ \left[ \begin{array}{@{}c@{}} \displaystyle {A}_{n} & {B}_{n}{C}_{\epsilon }\\ \displaystyle {O}_{1} & {A}_{dis} \end{array} \right] {}}^{A} \left[ \begin{array}{@{}c@{}} \displaystyle {x}_{n}(i)\\ \displaystyle h(i) \end{array} \right] & & +\overbrace{ \left[ \begin{array}{@{}c@{}} \displaystyle {B}_{n}\\ \displaystyle {O}_{2} \end{array} \right] {}}^{B}u(i)+ \left[ \begin{array}{@{}c@{}} \displaystyle {O}_{3}\\ \displaystyle {B}_{\epsilon } \end{array} \right] \epsilon (i)y(i) & =\overbrace{ \left[ \begin{array}{@{}l@{}} \displaystyle {C}_{n} & {O}_{4} \end{array} \right] {}}^{C} \left[ \begin{array}{@{}c@{}} \displaystyle {x}_{n}(i)\\ \displaystyle h(i) \end{array} \right] . \end{eqnarray*}\end{document}

Based on this model, we can use the disturbance compensation strategy to get \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\hat {x}}_{n}(i)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {\nu }(i)\end{document} from previous disturbance observer design section.

Deriving the controller K utilizing matrices An and Bn to ensure the stability of the closed-loop control system matrix An − BnK, the interim control variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tilde {u}(i)\end{document} is computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\tilde {u}(i)=-K{\hat {x}}_{n}(i).\end{eqnarray*}\end{document}

Secondly, get the observer gain Kob using matrices A and C allowing for the estimation of both the input disturbance ν(i) and the state variable xn(k) within the augmented state space model, resulting in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {\nu }(i)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\hat {x}}_{n}(i)\end{document} .

After that, deducting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {\nu }(i)\end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tilde {u}(i)\end{document} results in the control variable u(i):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}u(i)=\tilde {u}(i)-\hat {\nu }(i).\end{eqnarray*}\end{document}

Compare a fixed signal r(i) with our output signal, considering the reference signal r(i), where the observer formulation takes the following form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\hat {x}(i+1)=A\hat {x}(i)+Bu(i)+{K}_{ob}(y(i)-r(i)-C\hat {x}(i)). \end{eqnarray*}\end{document}

The (A, B, C) matrices are augmented matrices which contains the input disturbance compensation. Furthermore, it is imperative to validate that the entirety of the eigenvalues pertaining to the observer error system matrix A − KobC remains strictly bounded within the confines of the unit circle, thereby guaranteeing system stability and observability.

Determined the closed-loop eigenvalue of the state feedback control system based on disturbance observer to ensure the stability of the system. For both xn(i) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {p}(i)\end{document} , we adopt the observation error system

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \left[ \begin{array}{@{}c@{}} \displaystyle {E}_{xp} \left\{ {\tilde {x}}_{n}(i+1) \right\} \\ \displaystyle {E}_{xp}\{ \tilde {h}(i+1)\} \end{array} \right] = \left( A-{K}_{ob}C \right) \left[ \begin{array}{@{}c@{}} \displaystyle {E}_{xp} \left\{ {\tilde {x}}_{n}(i) \right\} \\ \displaystyle {E}_{xp}\{ \tilde {h}(i)\} \end{array} \right] , \end{eqnarray*}\end{document}

where Exp represents the error variables xn(i) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {h}(i)\end{document} , with their expectations denoted as respective values.

In addition, the proposed disturbance observer state feedback control strategy is compared to the traditional Proportional-Integral-Derivative (PID) control strategy based on important performance indicators, including trajectory tracking accuracy, response time, resilience to external disturbances, and algorithm complexity. At first, the proposed technique demonstrates superior precision in accurately following a desired path compared to the classic PID control, especially when dealing with external disturbances. Furthermore, it demonstrates greater proficiency in maneuvering intricate and ever-changing surroundings. In terms of response time, this technique efficiently adjusts control parameters to accommodate variations in external variables, surpassing the performance of conventional PID control. Our system, which includes the use of a disturbance observer, is highly effective in resisting external disturbances. It is particularly adept at withstanding environmental changes like strong winds, ensuring excellent flight stability and trajectory accuracy. This sets it apart from conventional methods. Furthermore, although the complexity of our method’s control algorithm exceeds that of PID control, its benefits in terms of other performance metrics justify this complexity.

Experiments

For simulation, the mass of the UAV is set to 0.5 kilograms, The scale factor of gravity is g = 9.80122 N/kg. The system model is linearized using the Taylor series expansion approach, specifically targeting equilibrium points and utilizing the Jacobian linearization methodology for simulation purposes. The sampling instant is chosen to 900.

In the context of state feedback control utilizing disturbance observers, we choose a specific polynomial D(j^−1^) as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}D({j}^{-1})=1-2\cos \nolimits ~{\omega }_{d}{j}^{-1}+{j}^{-2}, \end{eqnarray*}\end{document}

where ωd = 0.5. When the aerial reference trajectory follows a linear path, the term cosωd will be removed.

We assume the quadcopter system has the input disturbance on each axis, for x − and y −axes, the input disturbance μ(k) can be expressed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\nu (i)= \frac{\epsilon (i)}{ \left( 1-2\cos \nolimits ~{\omega }_{d}{j}^{-1}+{j}^{-2} \right) } .\end{eqnarray*}\end{document}

This is the main disturbance signal that we test in this paper. In addition, if the controller is an integrator type, it chooses step or random walk type disturbance signal 1 − q^−1^ as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} D \left( {q}^{-1} \right) \end{document} , if the controller is an integrator plus sinusoidal mode, D(j^−1^) will also be changed to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}(1-{j}^{-1})(1-2\cos \nolimits ~{\omega }_{d}{j}^{-1}+{j}^{-2}). \end{eqnarray*}\end{document}

According to the augmented state space equation, it can accomplish the state feedback control system with an estimation of the input disturbance.

Initially, it obtains the augmented state space model (A, B, C) with (An, Bn, Cn). Since we have 12 state variables, four input variables and four outputs for quadcopter modeling at the begining. (An, Bn, Cn) are 12 × 12, 12 × 4 and 4 × 12 matrices respectively. The dimensions of A, B, and C are 20 × 20, 20 × 4 and 4 × 20, respectively.

Utilizing the MATLAB function dlqr.m with Q as the identity matrix and R set to one, the state feedback controller is systematically designed based on the pair (An, Bn), yielding

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}K={ \left[ \begin{array}{@{}rrrr@{}} \displaystyle 0&\displaystyle 0&\displaystyle 0.3172&\displaystyle 0\\ \displaystyle 0&\displaystyle -0.3172&\displaystyle 0&\displaystyle 0\\ \displaystyle 0.9766&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0.3841&\displaystyle 0\\ \displaystyle 0&\displaystyle -0.3841&\displaystyle 0&\displaystyle 0\\ \displaystyle 1.1623&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 1.6566&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 1.6566&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 0.4875\\ \displaystyle 0&\displaystyle 0.3536&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0.3536&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 0.5004 \end{array} \right] }^{T}. \end{eqnarray*}\end{document}

This state feedback regulator K is used for manipulating the interim control parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tilde {u}(i)\end{document} , thereby offsetting the external disturbance ν(i) using the estimated parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {\nu }(i)\end{document} . In addition, ensure that the matrix An − BnK lies within the unit circle, thereby maintaining the stability of the closed-loop control system. Computing the intermediary control variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tilde {u}(i)\end{document} based on the previous equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\tilde {u}(i)=-K{\hat {x}}_{n}(i).\end{eqnarray*}\end{document}

With the augmented state space model strategy using the disturbance model, derive the observer gain Kob utilizing the augmented matrices A and C, followed by estimating both the input disturbance ν(i) and the state variable xn(i), resulting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {\nu }(i)\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} {\hat {x}}{n}(i)\end{document} , where the Kob_ equals to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*} \left[ \begin{array}{@{}rrrr@{}} \displaystyle 1.2350&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 1.2350&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0.8012&\displaystyle 0\\ \displaystyle 25.0089&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 25.0089&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 6.2022&\displaystyle 0\\ \displaystyle 0&\displaystyle -33.2883&\displaystyle 0&\displaystyle 0\\ \displaystyle 33.2883&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 1.3037\\ \displaystyle 0&\displaystyle -275.4516&\displaystyle 0&\displaystyle 0\\ \displaystyle 275.4516&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 28.1125\\ \displaystyle 0&\displaystyle 0&\displaystyle 14.1530&\displaystyle 0\\ \displaystyle 0&\displaystyle -10.4817&\displaystyle 0&\displaystyle 0\\ \displaystyle 10.4817&\displaystyle -0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 4.9595\\ \displaystyle 0&\displaystyle 0&\displaystyle 13.7672&\displaystyle 0\\ \displaystyle 0&\displaystyle -10.0984&\displaystyle 0&\displaystyle 0\\ \displaystyle 10.0984&\displaystyle 0&\displaystyle 0&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 4.3870 \end{array} \right] . \end{eqnarray*}\end{document}

Once got the intermediary control signal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \tilde {u}(i)\end{document} and estimated input \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \hat {\nu }(i)\end{document} , the previous difference equation is used to obtain u(i)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}u(i)=\tilde {u}(i)-\hat {\nu }(i).\end{eqnarray*}\end{document}

Comparing the output vector y(i) with the target signal r(i), the observer equation can be written as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}\hat {x}(i+1)=A\hat {x}(i)+Bu(i)+{K}_{ob}(y(i)-r(i)-C\hat {x}(i)). \end{eqnarray*}\end{document}

A sensitivity analysis is performed to evaluate the effectiveness of the control algorithm under different parameter settings. This is achieved by modifying various parameters of the quadcopter unmanned aerial vehicle, such as mass, moment of inertia, and aerodynamic properties, within a simulated environment. This study includes a wide range of parameter adjustments, including the mass moment of inertia.

The control approach extensively depends on the Dryden turbulence model to accurately measure the impact of wind disturbances. The Power Spectral Density atmospheric turbulence model represents wind disturbances as stochastic fluctuations in a continuous random process. White noise, which has a consistent power spectral density and equal strength across certain frequency ranges, is subjected to shaping filtration in order to create the intended turbulence effects.

Disturbances are intentionally produced in a predetermined direction for a quadcopter drone operating in space in our simulated scenario. The transfer functions of the filter are expressed as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{eqnarray*}{H}_{u}(s)={\sigma }_{u}\sqrt{ \frac{2{L}_{u}}{\pi V} } \frac{1}{1+ \frac{{L}_{u}}{V} s} ,\end{eqnarray*}\end{document}

where Lu and σu represent the scale length of turbulence and turbulence intensity, respectively; V represents the air velocity.

Once the white noise is filtered, it is included into the simulation model, and the resulting simulation results are thoroughly evaluated.

Results

Trajectory types

First of all, we have assumed three cases for the reference trajectory in 3D space. The first case is a single straight-line trajectory. The second case is the cylindrical spiral track, which contains sine and cosine curves and can represent the turning or hovering situation of the UAV in general. The third case is the second conic trajectory, which can be understood as a high-order extension in the second case. Because the conic curve contains elliptical, parabolic, and hyperbolic trajectories, it can represent more realistic UAV flight trajectories and has certain representativeness in the aviation field. The simulation results are presented below.

Reference following. Figure 1 shows the reference tracking performance with the integrator mode controller and random walk disturbances: D(j^−1^) = 1 − j^−1^.

Reference signal (brown) and following results (blue) with a steep slope track in the x-, y-, and z-axes.

Figure 2 illustrates the simulation results with a sinusoidal mode controller and the disturbance in the corresponding sinusoidal mode: D(j^−1^) = 1 − 2cos ωdj^−1^ + j^−2^.

Reference signal (brown) and following results (blue) with a cylindrical spiral track (A) in the x-, y-, and z-axes and (B) in 3D space.

Figure 3A shows that the reference signal is extended to a conic track with a variable turning radius. In addition, Fig. 3B shows an intuitive 3D space simulation result. Figure 1 to Fig. 3 show that for different reference trajectories and controllers, selecting appropriate corresponding disturbance models D(j^−1^) can complete the disturbance compensation issue and make the simulation results consistent with the reference trajectories.

Reference signal (brown) and following results (blue) with a conic track (A) in the x-, y-, and z-axes and (B) in 3D space.

To experiment with a contrast test on the premise of conic trajectory, if the D(j^−1^) model is changed to not correspond to the controller, the quadcopter will not be able to effectively complete the compensation of disturbance and trajectory tracking. The simulation results are shown in Fig. 4, taking the y −axis as an example.

Reference signal (brown) and following results (blue) based on the integral type D with conic track in the y-axis (A) and in 3D space (B).

Unmanned aerial vehicles may undergo alterations in their mass and moment of inertia during actual use, which can be caused by factors such as cargo variations, battery usage, or ambient conditions. Therefore, it is essential to evaluate the stability and efficiency of the control strategy in light of these variations in parameters. This study thoroughly examines the resilience of the suggested control approach when faced with different aircraft parameter fluctuations.

The analysis considers the fluctuation in the weight of the unmanned aerial vehicle, utilizing simulated scenarios that incorporate a gradual rise in weight up to 1.5 kg. Our simulation results show that the suggested control method maintains a high level of trajectory tracking performance even when there are changes in mass that affect how the system moves. This validates the efficacy of the technique in addressing fluctuations in mass. The study examines the impact of changes in the weight of the unmanned aerial vehicle, specifically focusing on cases where the weight decreases. Following a 30% decrease in mass, the trajectory tracking, as illustrated in Fig. 5A, continues to provide very good performance. The simulation effectively retains its tracking performance, similar to its performance before the drop in mass, using the y-axis as an illustration. It ultimately approaches the reference trajectory after approximately two sinusoidal oscillation cycles. The results of the simulation show that the suggested control method consistently achieves a high level of trajectory tracking performance, even when there are changes in mass that affect how the system moves.

Reference signal (brown) and following result (blue) with conic y-axis track presentation for the mass (A) and inertial (B) change conditions.

Furthermore, a thorough examination of the influence of changes in the moment of inertia on the control system is conducted. After a specific decrease of 15% in inertia, the trajectory tracking, as seen in Fig. 5B, exhibits a comparable pattern. Illustrating the y-axis as an instance, at a small scale, the simulation presently reaches a state of convergence with the reference trajectory after around three cycles of sinusoidal oscillation. At a large scale, the general ability to accurately follow a desired path continues to be exceptional, as the error values consistently stay within the range of 0.002 m.

Conclusions

In conclusion, this article presents a control method for tracking trajectories that utilizes disturbance-observer-based control. This method is showcased in the context of quadcopter trajectory control by simulating it across three different reference trajectories. The results demonstrate the method’s ability to precisely monitor sinusoidal trajectories in three-dimensional space, with particularly impressive performance in tracking conic curves. Significantly, its flexibility encompasses intricate flight paths for UAVs, even when faced with irregular turning radii. Moreover, its potential uses range from accurate position manipulation to broader effects in controlling UAVs in aviation. This approach offers favorable opportunities in the fields of Unmanned Aerial Vehicle technology and flight control.

Supplemental Information

10.7717/peerj-cs.1861/supp-1Supplemental Information 1MATLAB code and data used in the study