Quaternion-based state-dependent differential Riccati equation for quadrotor drones: Regulation control problem in aerobatic flight

Abstract

The quaternion is a powerful and common tool to avoid singularity in rotational dynamics in three-dimensional (3D) space. Here it has been particularly used as an alternative to Euler angles and rotation matrix. The application of the quaternion is exercised in quadrotor modeling and control. It changes the dynamics and represents a singularity-free attitude model. Here for the first time (for the best knowledge of authors), the state-dependent differential Riccati equation (SDDRE) control has been implemented on the quaternion-based model of a quadcopter. The proposed control structure is capable of aerobatic flight, and the Pugachev’s Cobra maneuver is chosen to assess the capability of the quaternion-based SDDRE approach. The introduced control simulator is validated by comparison with conventional dynamics based on Euler angles, controlled using a proportional-derivative (PD) controller on a normal regulation flight. The simulator successfully performed the Cobra maneuver and also validated the proposed structure. The more precision in regulation along with lower energy consumption demonstrated the superiority of the introduced approach.

1. Introduction

The Euler angles and rotation matrix in three-dimensional (3D) space are vulnerable to singularities, considering dynamics, especially in an aerobatic flight. Here we focus, particularly on multirotor drones and quadcopters. The dynamics of the multirotor drones are usually subjected to hovering assumptions to guarantee a stable flight [1]. To perform agile maneuvers with sudden changes in attitude, geometric control was introduced that avoids singularities and computes the rotation matrix in another manifold [2]. Krishna et al. used the geometric control for helicopter trajectory tracking in agile flight regarding the attitude of the system [3]. The sliding mode was employed as the controller and the error function in the geometric domain constructed the sliding condition to prove the stability. The geometric control was also applied on quadrotor drones subject to wind disturbance and sudden changes in attitude dynamics [4]. A flip is representative of a challenging maneuver that includes passing through a singularity that was addressed by the geometric approach [5]. The $ \pi \left(\textrm{r}\textrm{a}\textrm{d}\right)$ flip in roll angle was considered, sudden motion between two stable roll angles. The change in the manifold and solving singularity cost the user more complicated integration methods [6]. Variational integration is an effective method [7]. Improper treatment of the integration may cause numerical drift in the results, especially in big orientation changes.

The quaternion is another representation of the complex numbers in mathematics, with a wide range of usage in theory and application. The focus of this work is to implement the alternative representation of the Euler angles and rotation matrix in 3D space. The application of quaternion in control was reported on aircraft [8, 9], orientation and translation control of manipulators [10, 11], control of autonomous underwater vehicles [12–15], helicopter attitude control [16], etc. Terze et al. used quaternion representation of the rotational dynamics for aircraft simulators and introduced shifting update process to ensure precise integration in long flight simulations [8]. Kinematics control of a robotic arm was presented using dual quaternion, aimed to present a robust controller without decoupling of the translational and orientation dynamics [10]. Cooke et al. presented a thorough implementation of the quaternion dynamics for flight simulation [17]. The performance of the quaternion-based control was validated by an omnidirectional intelligent navigator for an underwater platform [12]. Suzuki et al. presented the simulation and experimental studies on Lepton-Ex unmanned helicopter using quaternion feedback in control [16]. Quadrotor and multicopters were also employed quaternions [18–23]. A study compared three controllers, linear quadratic regulator (LQR), proportional-derivative (PD), and model predictive controller (MPC) for trajectory tracking that revealed the PD and the LQR gained better results in an ideal condition, and the MPC was more accurate in presence of disturbance [18]. Quaternion variables in the attitude controller provided an advantage, employment of low-cost sensors due to the high-speed capability of the singularity-free controller [20]. Palomo et al. presented an observer-based method based on position-yaw for the ellipsoid method [21]. Linear matrix inequalities were used to optimize the feedback of the observer, and experimental results showed successful tracking in high speed maneuver. Euler angles are popular in UAV modeling though they suffer from gimbal lock when two orthogonal axes align and lock together [24]. Another disadvantage of Euler angles is the computational cost by computing so many trigonometric functions [25], on the contrary, quaternion mathematics only involves algebraic computations [26]. Sanchez et al. used quaternion dynamics for quadrotor control based on receiving gesture commands [27], where it was important to cope with the unpredictability of the gesture and quaternion made it more reliable due to insensitivity to singularity.

Performing experiment on aggressive maneuvers requires more safety concerns in addition to solving singularities. Gillula et al. used the Hamilton–Jacobi differential game approach for finding a reachable set for aerobatic maneuvers to guarantee safety [28]. Considering the constraint of the actuators in the flight experiment is an important issue since, during the flip, it is highly probable to put the rotors in saturation [29]. Here in this work, the simulation of the Cobra maneuver is done and for experimentation, actuator limits and speed of the maneuver must be considered to avoid undesired saturations in the middle of trajectory that could be dangerous.

The state-dependent Riccati equation (SDRE) is a closed-loop optimal nonlinear controller introduced in the 1960s [30]. The utilization of the quaternion in SDRE was explored in several platforms such as satellite [31, 32], spacecraft [33–36], spacecraft in proximity operations [37], attitude control of a rigid body [38], and remote sensing CubeSats [39]. Here we present the state-dependent differential Riccati equation (SDDRE) control based on quaternion for quadcopter dynamic systems. For the best knowledge of authors, a quaternion-based SDDRE controller has not been used for quadrotors in the literature so far which makes the first contribution of this work. The SDDRE is a finite time controller that penalizes the states (error variables) by a final boundary condition [40].

The singularity-free attitude control is the principal advantage of the quaternion. The Cobra is a challenging aerobatic flight maneuver, performed by a jet aircraft [41, 42]. The motion turns the aircraft vertical (even for pitch angle $ \theta \gt \pi /2$ ) to perform the maneuver along with sudden deceleration; the thrust of the jet engine helps to avoid the system from falling.

A tail-sitter drone is a good choice, with a vertical thrust option, to perform the Cobra maneuver [43]. Xu et al. presented iterative learning control for a tail sitter unmanned aerial vertical-take-off-and-landing system for Pugachev’s Cobra maneuver [43]. They used the acceleration model that resulted in simple dynamics without system identification. The Cobra was exercised by a quadrotor using adaptive control [44]. The quadrotor possessed $ 28(\textrm{g})$ , a small lightweight platform. An adaptive backstepping controller with a modified recursive least-square was employed to control the system.

Contributions: (1) presenting a quaternion-based SDDRE control peculiar to a quadrotor. (2) Implementing Cobra maneuver in pitch angle in a forward flight. The performed Cobra in ref. [44] was done in vertical ascending flight. In a forward flight, conducting a Cobra maneuver, the system will lose its thrust ( $ \theta \approx \pi /2$ ) and is subjected to fall. Here in this work, using SDDRE, the Cobra in forward flight is done.

Section 2 describes the preliminaries in quaternion mathematics. Section 3 states the system dynamics details. Section 4 presents the SDDRE control structure. Section 5 expresses the approximate closed-form solution to the SDDRE. Section 6 presents the implementation and method and cascade design. Simulations are reported in Section 7 which includes validation and aerobatic flight, and concluding remarks are summarized in Section 8.

Notations: $ {\mathbb{R}}^{n}$ denotes the $ n$ -dimensional Euclidean space, $ {\mathbb{H}}^{n}$ denotes the $ n$ -dimensional Hamilton space, and ${\left(\cdot \right)}^{*}$ performs conjugate transpose. $ {\mathbb{R}}^{n\times m}$ is the set of $ n\times m$ real matrices; $ {\left(\cdot \right)}^{T}$ is the transpose of a matrix or a vector; $\bigotimes$ denotes Kronecker product; $ \textrm{d}\textrm{i}\textrm{a}\textrm{g}\left(\cdot \right)$ means a diagonal matrix; $ {\textbf{I}}_{n\times n}$ and ${\textbf{0}}_{n\times n}$ denote $ n\times n$ identity and zero matrices.

2. Preliminaries: Quaternion mathematics

Here we consider the quaternion definition with real-scalar part $ {q}_{0}$ and vector-imaginary part $ {\textbf{q}}_{\textrm{v}}={q}_{1}\textbf{i}+{q}_{2}\textbf{j}+{q}_{3}\textbf{k}$ , set all together

(1) \begin{align}\textbf{q}=\left[\begin{array}{c}{q}_{0}\\ \\[-7pt] {\textbf{q}}_{\textrm{v}}\end{array}\right]=\left[\begin{array}{c}\textrm{cos}\dfrac{\vartheta }{2}\\ \\[-7pt] \textbf{r}\textrm{sin}\dfrac{\vartheta }{2}\end{array}\right]\!,\end{align}

where $ \textbf{r}\in {\mathbb{R}}^{3}$ is a unit rotation vector and $ \vartheta $ is the corresponding rotation angle about that ref. [20]. To solve the ambiguity of the direction of quaternions, $ 0\le {q}_{0}\le 1$ is chosen. A conjugate transpose of the quaternion (1) is presented as

(2) \begin{align}{\textbf{q}}^{\textbf{*}}=\left[\begin{array}{c}{q}_{0}\\ \\[-7pt] -{\textbf{q}}_{\textrm{v}}\end{array}\right].\end{align}

The multiplication product of two arbitrary quaternions $ \textbf{q}$ and $ \textbf{p}$ is defined through Kronecker product

(3) \begin{align} \nonumber\\[-25pt]\textbf{q}\otimes \textbf{p}=\textbf{Q}(\textbf{q})\textbf{p}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{q}_{0} & -{q}_{1} & -{q}_{2} &-{q}_{3}\\ \\[-7pt] {q}_{1} &{q}_{0} &-{q}_{3} & {q}_{2}\\ \\[-7pt] {q}_{2} & {q}_{3} & {q}_{0} & -{q}_{1}\\ \\[-7pt] {q}_{3} & -{q}_{2} & {q}_{1} &{q}_{0}\end{array}\right]\left[\begin{array}{c}{p}_{0}\\ \\[-7pt] {p}_{1}\\ \\[-7pt] {p}_{2}\\ \\[-7pt] {p}_{3}\end{array}\right].\end{align}

A unit quaternion can build a rotation transformation by two multiplications by Kronecker product that could rotate an arbitrary vector $ \textbf{v}$ from the global coordinate to a moving coordinate $ \textbf{q}$ as in the form of $ \textbf{q}\otimes {\left[0,{\textbf{v}}^{T}\right]}^{T}\otimes {\textbf{q}}^{\textrm{*}}$ . So, using definitions (2) and (3), the rotation matrix is built by replacing $ x,y,z$ in turn into $ \textbf{v}$ :

\begin{align*}{\textbf{R}}_{x}(\textbf{q})=\textbf{q}\otimes {\left[\textrm{0,1},\textrm{0,0}\right]}^{T}\otimes {\textbf{q}}^{\textbf{*}},\\[3pt]{\textbf{R}}_{y}(\textbf{q})=\textbf{q}\otimes {\left[\textrm{0,0,1,0}\right]}^{T}\otimes {\textbf{q}}^{\textbf{*}},\\[3pt]{\textbf{R}}_{z}(\textbf{q})=\textbf{q}\otimes {\left[\textrm{0,0,0,1}\right]}^{T}\otimes {\textbf{q}}^{\textbf{*}},\end{align*}

which form [25]:

\begin{align*}\textbf{R}(\textbf{q})=\left[{\textbf{R}}_{x}(2:\ 4),{\textbf{R}}_{y}(2:\ 4),{\textbf{R}}_{z}(2:\ 4)\right]=\left[\begin{array}{c@{\quad}c@{\quad}c}{q}_{0}^{2}+{q}_{1}^{2}-{q}_{2}^{2}-{q}_{3}^{2} & 2{q}_{1}{q}_{2}-2{q}_{0}{q}_{3} & 2{q}_{0}{q}_{2}+2{q}_{1}{q}_{3}\\ \\[-7pt] 2{q}_{0}{q}_{3}+2{q}_{1}{q}_{2} &{q}_{0}^{2}-{q}_{1}^{2}+{q}_{2}^{2}-{q}_{3}^{2} & 2{q}_{2}{q}_{3}-2{q}_{0}{q}_{1}\\ \\[-7pt] 2{q}_{1}{q}_{3}-2{q}_{0}{q}_{2} &2{q}_{0}{q}_{1}+2{q}_{2}{q}_{3} & {q}_{0}^{2}-{q}_{1}^{2}-{q}_{2}^{2}+{q}_{3}^{2}\end{array}\right].\end{align*}

where that is $ {\textbf{R}}_{x}(2:\ 4)$ selects three last components of $ {\textbf{R}}_{x}$ .

The angular velocity vector of the quaternion is accessible (supposedly) in local moving coordinate; it is presented by the vector $ {\boldsymbol\omega }\left(\textbf{q},\dot{\textbf{q}}\right)\!:\mathbb{H}\times {\mathbb{R}}^{4}\to {\mathbb{R}}^{3}$ where $ {\boldsymbol\omega }={\left[{\omega }_{1},{\omega }_{2},{\omega }_{3}\right]}^{T}\left(\dfrac{\textrm{r}\textrm{a}\textrm{d}}{\textrm{s}}\right)$ . Then the derivative of the quaternion is found [45]:

(4) \begin{align}{\dot{\textbf{q}}}_{{\boldsymbol\omega }}=\dfrac{1}{2}\textbf{q}\otimes \left[\begin{array}{c}0\\ \\[-7pt] \boldsymbol\omega \end{array}\right]=\frac{1}{2}\textbf{Q}(\textbf{q})\left[\begin{array}{c}0\\ \\[-7pt] \boldsymbol\omega \end{array}\right]\!,\end{align}

in which $ {\dot{\textbf{q}}}_{{\boldsymbol\omega }}(\textbf{q},{\boldsymbol\omega }):\ \mathbb{H}\times {\mathbb{R}}^{3}\to {\mathbb{R}}^{4}$ and $ \textbf{Q}(\textbf{q})$ has been introduced in Eq. (3).

The relation between the Euler angles $ (\varphi ,\theta ,\psi )$ and quaternions is also needed for the control design [25]:

(5) \begin{align}\textbf{q}=\left[\begin{array}{c}\textrm{cos}\dfrac{\varphi }{2}\textrm{cos}\dfrac{\theta }{2}\textrm{cos}\dfrac{\psi }{2}+\textrm{sin}\dfrac{\varphi }{2}\textrm{sin}\dfrac{\theta }{2}\textrm{sin}\dfrac{\psi }{2}\\ \\[-7pt] \textrm{sin}\dfrac{\varphi }{2}\textrm{cos}\dfrac{\theta }{2}\textrm{cos}\dfrac{\psi }{2}-\textrm{cos}\dfrac{\varphi }{2}\textrm{sin}\dfrac{\theta }{2}\textrm{sin}\dfrac{\psi }{2}\\ \\[-7pt] \textrm{cos}\dfrac{\varphi }{2}\textrm{sin}\dfrac{\theta }{2}\textrm{cos}\dfrac{\psi }{2}+\textrm{sin}\dfrac{\varphi }{2}\textrm{cos}\dfrac{\theta }{2}\textrm{sin}\dfrac{\psi }{2}\\ \\[-7pt] \textrm{cos}\dfrac{\varphi }{2}\textrm{cos}\dfrac{\theta }{2}\textrm{sin}\dfrac{\psi }{2}-\textrm{sin}\dfrac{\varphi }{2}\textrm{sin}\dfrac{\theta }{2}\textrm{cos}\dfrac{\psi }{2}\end{array}\right]\!,\end{align}

and also with inverse mapping, one could find:

\begin{align*}\left[\begin{array}{c}\varphi \\ \\[-7pt] \theta \\ \\[-7pt] \psi \end{array}\right]=\left[\begin{array}{c}\textrm{a}\textrm{t}\textrm{a}\textrm{n}2\left(2\left({q}_{0}{q}_{1}+{q}_{2}{q}_{3}\right),{q}_{0}^{2}-{q}_{1}^{2}-{q}_{2}^{2}+{q}_{3}^{2}\right)\\ \\[-7pt] \textrm{a}\textrm{s}\textrm{i}\textrm{n}\left(2\left({q}_{0}{q}_{2}-{q}_{1}{q}_{3}\right)\right)\\ \\[-7pt] \textrm{a}\textrm{t}\textrm{a}\textrm{n}2\left(2\left({q}_{0}{q}_{3}+{q}_{1}{q}_{2}\right),{q}_{0}^{2}+{q}_{1}^{2}-{q}_{2}^{2}-{q}_{3}^{2}\right)\end{array}\right].\end{align*}

3. Dynamics of the system

Consider a “plus-shaped” quadrotor drone with two moving and fixed reference coordinates, body frame denoted by $ \mathcal{B}$ , and global frame marked with $ \mathcal{M}=\left\{X,Y,Z\right\}$ , with respect, see Fig. 1. The position of the moving coordinate is defined through the vector $ {{\boldsymbol\xi }}_{1}(t)={\left[{x}_{\textrm{c}}(t),{y}_{\textrm{c}}(t),{z}_{\textrm{c}}(t)\right]}^{T}\left(\textrm{m}\right).$ The kinematics equation is

(6) \begin{align}{\dot{{\boldsymbol\xi }}}_{1}(t)=\textbf{R}(\textbf{q}(t)){{\boldsymbol\upsilon }}_{1}(t),\end{align}

where $ {{\boldsymbol\upsilon }}_{1}(t)={\left[u(t),v(t),w(t)\right]}^{T}\left(\frac{\textrm{m}}{\textrm{s}}\right)$ , and $ \textbf{R}(\textbf{q}):\mathbb{H}\to {\mathbb{R}}^{3\times 3}$ is the quaternion-based rotation matrix; and $ \left\{\varphi ,\theta ,\psi \right\}\left(\textrm{r}\textrm{a}\textrm{d}\right)$ are Euler angles set in global coordinate. The local angular velocity vector set on the body frame is also named $ {{\boldsymbol\upsilon }}_{2}(t)={\left[p(t),q(t),r(t)\right]}^{T}\left(\frac{\textrm{r}\textrm{a}\textrm{d}}{\textrm{s}}\right)$ .

Fig. 1. The definition of the reference coordinates for a sample quadrotor drone.

To find the dynamics of the system, Newton–Euler equation could be used that results in two sets of dynamic equations, the first set is translational:

(7) \begin{align}m{\ddot{{\boldsymbol\xi }}}_{1}(t)={\textbf{R}}_{3}(\textbf{q}(t)){T}_{\textrm{B}}(t)-mg{\textbf{e}}_{3},\end{align}

where $ {\textbf{R}}_{3}(\textbf{q})$ is the last column of $ \textbf{R}(\textbf{q})$ , $ g=9.81\left(\frac{\textrm{m}}{{\textrm{s}}^{2}}\right)$ is the gravity constant, $ m\left(\textrm{k}\textrm{g}\right)$ is the total mass of the drone, and $ {\textbf{e}}_{3}={\left[\textrm{0,0,1}\right]}^{T}$ . The second set is rotational dynamic:

(8) \begin{align}\textbf{J}{\dot{{\boldsymbol\upsilon }}}_{2}(t)={{\boldsymbol\tau }}_{\textrm{B}}(t)-{{\boldsymbol\upsilon }}_{2}(t)\times \textbf{J}{{\boldsymbol\upsilon }}_{2}(t),\end{align}

where $ {{\boldsymbol\tau }}_{\textrm{B}}(t)={\left[{\tau }_{x}(t),{\tau }_{y}(t),{\tau }_{z}(t)\right]}^{T}\left(\textrm{N}\textrm{m}\right)$ is the input torque vector and $ \textbf{J}=\textrm{d}\textrm{i}\textrm{a}\textrm{g}\left({I}_{xx},{I}_{yy},{I}_{zz}\right)\textrm{}\left(\textrm{k}\textrm{g}{\textrm{m}}^{2}\right)$ is the inertia matrix assigned in the body frame.

The state-vector is

\begin{align*}{\left[\textbf{x}(t)\right]}_{13\times 1}={\left[{{\boldsymbol\xi }}_{1}^{T}(t),{\textbf{q}}^{T}(t),{\dot{{\boldsymbol\xi }}}_{1}^{T}(t),{\dot{{\boldsymbol\upsilon }}}_{2}^{T}(t)\right]}^{T},\end{align*}

which provides the state-space representation of the multi-copter using Eqs. (4)–(8) and setting $ {\boldsymbol\omega }(t)={{\boldsymbol\upsilon }}_{2}(t)={\left[p,q,r\right]}^{T}$ in Eq. (4):

(9) \begin{align}\dot{\textbf{x}}(t)=\left[\begin{array}{c}\textbf{R}(\textbf{q}(t)){{\boldsymbol\upsilon }}_{1}(t)\\ \\[-7pt] \dfrac{1}{2}\textbf{Q}(\textbf{q}(t))\left[\begin{array}{c}0\\ \\[-7pt] {{\boldsymbol\upsilon }}_{2}(t)\end{array}\right]\\ \\[-7pt] \dfrac{1}{m}\left({\textbf{R}}_{3}(\textbf{q}(t)){T}_{\textrm{B}}(t)-mg{\textbf{e}}_{3}-\textbf{D}{\dot{{\boldsymbol\xi }}}_{1}(t)\right)\\ \\[-7pt] {\textbf{J}}^{-1}({{\boldsymbol\tau }}_{\textrm{B}}(t)-{{\boldsymbol\upsilon }}_{2}(t)\times \textbf{J}{{\boldsymbol\upsilon }}_{2}(t))\end{array}\right].\end{align}

In state-space Eq. (9), $ \textbf{D}\in {\mathbb{R}}^{3\times 3}$ collects drag and aerodynamics parameters of the quadcopter model and has been added to complete the model.

4. The state-dependent differential Riccati equation controller design

Consider the time-invariant affine-in-control nonlinear system

(10) \begin{align}\dot{\textbf{x}}(t)=\textbf{f}(\textbf{x}(t))+\textbf{g}\left(\textbf{x}(t),\textbf{u}(t)\right),\end{align}

where $ \textbf{x}(t)\in {\mathbb{R}}^{n}$ is the state vector, $ \textbf{u}(t)\in {\mathbb{R}}^{m}$ is the input vector, $ \textbf{f}(\textbf{x}(t)):\ {\mathbb{R}}^{n}\to {\mathbb{R}}^{n}$ and $ \textbf{g}(\textbf{x}(t),\textbf{u}(t)):\ {\mathbb{R}}^{n}\times {\mathbb{R}}^{m}\to {\mathbb{R}}^{n}$ represents the dynamics and they satisfy local Lipschitz condition. $ n$ is the dimension of the state vector and $ m$ is the total number of actuators. The system Eq. (10) is transformed into state-dependent coefficient (SDC) parameterization [46]:

(11) \begin{align}\textbf{f}(\textbf{x}(t))=\textbf{A}(\textbf{x}(t))\textbf{x}(t),\end{align}

(12) \begin{align}\textbf{g}(\textbf{x}(t),\textbf{u}(t))=\textbf{B}(\textbf{x}(t))\textbf{u}(t),\\[7pt] \nonumber \end{align}

in which $ \textbf{B}(\textbf{x}(t)):\ {\mathbb{R}}^{n}\to {\mathbb{R}}^{n\times m}$ and $ \textbf{A}(\textbf{x}(t)):\ {\mathbb{R}}^{n}\to {\mathbb{R}}^{n\times n}$ are held. The pair of $ \left\{\textbf{A}(\textbf{x}(t)),\textbf{B}(\textbf{x}(t))\right\}$ is a controllable parameterization of system (10) [47].

The cost function of the SDDRE is structured as

(13) \begin{align}J\left(\cdot \right)=\frac{1}{2}{\textbf{x}}^{T}(t)\textbf{F}\textbf{x}(t)+\frac{1}{2}\underset{0}{\overset{{t}_{\textrm{f}}}{\int }}\left[{\textbf{x}}^{T}(t)\textbf{Q}(\textbf{x}(t))\textbf{x}(t)+{\textbf{u}}^{T}(t)\textbf{R}(\textbf{x}(t))\textbf{u}(t)\right]\textrm{d}t,\end{align}

where ${t}_{\textrm{f}}\in {\mathbb{R}}^{+}$ is the final time of the control task, $ \textbf{Q}(\textbf{x}(t)):\ {\mathbb{R}}^{n}\to {\mathbb{R}}^{n\times n}$ , $ \textbf{R}(\textbf{x}(t)){\mathbb{R}}^{n}\to {\mathbb{R}}^{m\times m}$ , and $ \textbf{F}\in {\mathbb{R}}^{n\times n}$ are weighting matrices, for states and inputs in $ t\in [0,{t}_{\textrm{f}})$ and states at the final time $ {t}_{\textrm{f}}$ with respect. The pair of $ \left\{\textbf{A}(\textbf{x}(t)),{\textbf{Q}}^{\frac{1}{2}}(\textbf{x}(t))\right\}$ is an observable parameterization of system (10) where $ {\textbf{Q}}^{\frac{1}{2}}(\textbf{x}(t))$ is the Cholesky decomposition of weighting matrix in (13).

The control law is defined by applying $\frac{\partial \mathcal{H}\left(\textbf{x}(t),\textbf{u}(t),{\boldsymbol\lambda }(t)\right)}{\partial \textbf{u}(t)}=\textbf{0}$ , stationary condition, on Hamiltonian function $\mathcal{H}\left(\textbf{x}(t),\textbf{u}(t),{\boldsymbol\lambda }(t)\right)=J\left(\cdot \right)+{{\boldsymbol\lambda }}^{T}(t)\dot{\textbf{x}}(t)$ :

(14) \begin{align}\textbf{u}(t)=-{\textbf{R}}^{-1}(\textbf{x}(t)){\textbf{B}}^{T}(\textbf{x}(t))\textbf{K}(\textbf{x}(t))\textbf{x}(t),\end{align}

where $ {\boldsymbol\lambda }(t)=\textbf{K}(\textbf{x}(t))\textbf{x}(t)$ is the co-state vector and $\textbf{K}(\textbf{x}(t)):\ {\mathbb{R}}^{n}\to {\mathbb{R}}^{n \times n}$ is the control gain of the control equation, SDDRE.

Using the necessary condition for optimality, $ \frac{\partial \mathcal{H}\left(\textbf{x}(t),\textbf{u}(t),{\boldsymbol\lambda }(t)\right)}{\partial \textbf{x}(t)}=-\dot{{\boldsymbol\lambda }}(t)$ and the derivative of the co-state vector $ \dot{{\boldsymbol\lambda }}=\dot{\textbf{K}}\textbf{x}+\textbf{K}\dot{\textbf{x}}$ , one could find:

\begin{align*}\dot{\textbf{K}}(\textbf{x})+\textbf{K}(\textbf{x})\textbf{A}(\textbf{x})-\textbf{K}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})\textbf{K}(\textbf{x})+{\textbf{A}}^{T}(\textbf{x})\textbf{K}(\textbf{x})+\textbf{Q}(\textbf{x})+{\boldsymbol\Pi }(\textbf{x})=\textbf{0},\end{align*}

where $ {\boldsymbol\Pi }(\textbf{x})$ collects the derivative terms (see the complete equation in ref. [46]), then the SDDRE is represented as [48]:

(15) \begin{align}-\dot{\textbf{K}}(\textbf{x})=\textbf{Q}(\textbf{x})-\textbf{K}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})\textbf{K}(\textbf{x})+\textbf{K}(\textbf{x})\textbf{A}(\textbf{x})+{\textbf{A}}^{T}(\textbf{x})\textbf{K}(\textbf{x}),\end{align}

with $ \textbf{K}\left(\textbf{x}\left({t}_{\textrm{f}}\right)\right)=\textbf{F}$ , a final boundary condition.

5. An approximate closed-form solution to the SDDRE

The quadrotor dynamics in Eq. (9) must be controlled by the SDDRE approach. The extended linearization model of (9) is so-called state-dependent coefficient parameterization (or apparent linearization) [49]. To solve the SDDRE (15) and to find the control gain, several methods could be used such as backward integration (BI) [50], state-transition matrix [51], and Lyapunov-based approach [52]. The BI imposes a two-round solution that might not be proper for systems that need frequent changes in initial and final conditions, nonrepetitive systems. The STM approach is not computationally robust when the final time is long [53], then the Lyapunov-based approach has been selected for this work. It works with both positive and negative solutions to the Riccati equation and delivers an approximate closed-form answer, in just one round. This work uses Lyapunov-based method via negative root to the related SDRE, $ {\textbf{K}}_{\textrm{s}\textrm{s}}^{-}(\textbf{x})$ to find the symmetric positive-definite solution to the SDDRE, $ \textbf{K}(\textbf{x}(t))$ [50]. Subtracting (15) from

(16) \begin{align}{\textbf{A}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})+{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\textbf{A}(\textbf{x})+\textbf{Q}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})=\textbf{0},\end{align}

generates

(17) \begin{align}-\dot{\textbf{K}}(\textbf{x})&=\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\textbf{A}(\textbf{x})+{\textbf{A}}^{T}(\textbf{x})\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right] \nonumber\\[3pt]&\quad -\textbf{K}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})\textbf{K}(\textbf{x})\nonumber\\[3pt]&\quad +{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x}).\end{align}

Adding and subtracting $ \textbf{K}\textbf{B}{\textbf{R}}^{-1}{\textbf{B}}^{T}{\textbf{K}}_{\textrm{s}\textrm{s}}$ , $ {\textbf{K}}_{\textrm{s}\textrm{s}}\textbf{B}{\textbf{R}}^{-1}{\textbf{B}}^{T}\textbf{K}$ and $ {\textbf{K}}_{\textrm{s}\textrm{s}}\textbf{B}{\textbf{R}}^{-1}{\textbf{B}}^{T}{\textbf{K}}_{\textrm{s}\textrm{s}}$ to (17) result in:

(18) \begin{align}-\dot{\textbf{K}}(\textbf{x})&={\textbf{A}}^{T}(\textbf{x})\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]+\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\textbf{A}(\textbf{x})\nonumber\\&\quad -\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\nonumber\\&\quad -{\left[\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]}^{T}\left[\textbf{K}\right(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x}\left)\right]\nonumber\\&\quad -\textbf{K}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})\textbf{K}(\textbf{x})\nonumber\\&\quad +\textbf{K}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\nonumber\\&\quad +{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})\textbf{K}(\textbf{x})\nonumber\\&\quad -{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x}).\end{align}

Rewriting (18):

\begin{align*}-\dot{\textbf{K}}(\textbf{x})&=\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\textbf{A}(\textbf{x})+{\textbf{A}}^{T}(\textbf{x})\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right] \\&\quad -\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\\&\quad -{\left[\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]}^{T}\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\\&\quad -\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})\left[\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x})\right]\!,\end{align*}

and introducing a new variable

\begin{align*}{\textbf{P}}^{-1}(\textbf{x})=\textbf{K}(\textbf{x})-{\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x}),\end{align*}

and expressing the closed-loop matrix of the system

\begin{align*}{\textbf{A}}_{\text{cl}}(\textbf{x})=\textbf{A}(\textbf{x})-\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{K}}_{\textrm{s}\textrm{s}}(\textbf{x}),\end{align*}

one could present

(19) \begin{align}-\dot{\textbf{K}}(\textbf{x})={\textbf{P}}^{-1}(\textbf{x}){\textbf{A}}_{\text{cl}}(\textbf{x})+{\textbf{A}}_{\text{cl}}^{T}(\textbf{x}){\textbf{P}}^{-1}(\textbf{x})-{\textbf{P}}^{-1}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}){\textbf{P}}^{-1}(\textbf{x})\end{align}

Regarding that $ \frac{\textrm{d}}{\textrm{d}t}\left({\textbf{P}}^{-1}(\textbf{x})\right)=-{\textbf{P}}^{-1}(\textbf{x})\dot{\textbf{P}}(\textbf{x}){\textbf{P}}^{-1}(\textbf{x})$ , Eq. (19) should be changed to

\begin{align*}{\textbf{P}}^{-1}(\textbf{x})\dot{\textbf{P}}(\textbf{x}){\textbf{P}}^{-1}(\textbf{x})={\textbf{P}}^{-1}(\textbf{x}){\textbf{A}}_{\text{cl}}(\textbf{x})+{\textbf{A}}_{\text{cl}}^{T}(\textbf{x}){\textbf{P}}^{-1}(\textbf{x})-{\textbf{P}}^{-1}(\textbf{x})\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{}(\textbf{x}){\textbf{P}}^{-1}(\textbf{x}),\end{align*}

and consequently, results in a state-dependent differential Lyapunov Eq. [50]:

(20) \begin{align}\dot{\textbf{P}}(\textbf{x})={\textbf{A}}_{\text{cl}}(\textbf{x})\textbf{P}(\textbf{x})+\textbf{P}(\textbf{x}){\textbf{A}}_{\text{cl}}^{T}(\textbf{x})-\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x}),\end{align}

with a final boundary condition $ \textbf{P}\left({t}_{\text{f}}\right)={\left[\textbf{F}-{\textbf{K}}_{\text{ss}}(\textbf{x}(t))\right]}^{-1}$ . A solution to (20) is

(21) \begin{align}\textbf{P}(\textbf{x})=\textbf{E}(\textbf{x})+\textrm{e}\textrm{x}\textrm{p}\left\{{\textbf{A}}_{\text{cl}}(\textbf{x})\left(t-{t}_{\text{f}}\right)\right\}\left[\textbf{P}\left({t}_{\text{f}}\right)-\textbf{E}(\textbf{x})\right]\textrm{e}\textrm{x}\textrm{p}\left\{{\textbf{A}}_{\text{cl}}^{T}(\textbf{x})\left(t-{t}_{\text{f}}\right)\right\},\end{align}

in which $ \textbf{E}(\textbf{x}(t))$ is an answer to a state-dependent algebraic Lyapunov equation:

(22) \begin{align}\textbf{E}(\textbf{x}){\textbf{A}}_{\text{cl}}^{T}(\textbf{x})+{\textbf{A}}_{\text{cl}}(\textbf{x})\textbf{E}(\textbf{x})-\textbf{B}(\textbf{x}){\textbf{R}}^{-1}(\textbf{x}){\textbf{B}}^{T}(\textbf{x})=\textbf{0}.\end{align}

Proof of (21) could be checked by the substitution of Eq. (21) into (20):

(23) \begin{align}&{\dot{\textbf{E}}}({\textbf{x}}) + {{\textrm{d}\left( {{{\textbf{A}}_{\textrm{cl}}}({\textbf{x}})\left( {t - {t_{\textrm{f}}}} \right)} \right)} \over {\textrm{d}t}}{\textrm{exp}}\{ {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}})(t - {t_\textrm{f}})\} [{\textbf{P}}({t_\textrm{f}}) - {\textbf{E}}({\textbf{x}})]{\textrm{exp}}\{ {\textbf{A}}_{\textrm{cl}}^T({\textbf{x}})(t - {t_\textrm{f}})\}\nonumber\\[3pt]&\quad - {\textrm{exp}}\{ {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}})(t - {t_\textrm{f}})\} {\dot{\textbf{E}}}({\textbf{x}}){\textrm{exp}}\{ {\textbf{A}}_{\textrm{cl}}^T({\textbf{x}})(t - {t_\textrm{f}})\}\nonumber\\[3pt] &\quad + {\textrm{exp}}\{ {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}})(t - {t_\textrm{f}})\} [{\textbf{P}}({t_\textrm{f}}) - {\textbf{E}}({\textbf{x}})]{\textrm{exp}}\{ {\textbf{A}}_{\textrm{cl}}^T({\textbf{x}})(t - {t_\textrm{f}})\} {{\textrm{d}\left( {{\textbf{A}}_{\textrm{cl}}^T({\textbf{x}})\left( {t - {t_\textrm{f}}} \right)} \right)} \over {\textrm{d}t}}\nonumber\\[3pt] &= {\textbf{E}}({\textbf{x}}){\textbf{A}}_{\textrm{cl}}^T({\textbf{x}}) + {\textrm{exp}}\{ {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}})(t - {t_\textrm{f}})\} [{\textbf{P}}({t_\textrm{f}}) - {\textbf{E}}({\textbf{x}})]{\textrm{exp}}\{ {\textbf{A}}_{\textrm{cl}}^T({\textbf{x}})(t - {t_\textrm{f}})\} {\textbf{A}}_{\textrm{cl}}^T({\textbf{x}}) + {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}}){\textbf{E}}({\textbf{x}})\nonumber\\[3pt] &\quad + {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}}){\textrm{exp}}\{ {{\textbf{A}}_{\textrm{cl}}}({\textbf{x}})(t - {t_\textrm{f}})\} [{\textbf{P}}({t_\textrm{f}}) - {\textbf{E}}({\textbf{x}})]{\textrm{exp}}\{ {\textbf{A}}_{\textrm{cl}}^T({\textbf{x}})(t - {t_\textrm{f}})\} - {\textbf{B}}({\textbf{x}}){{\textbf{R}}^{ - 1}}({\textbf{x}}){{\textbf{B}}^T}({\textbf{x}}).\end{align}

From (22) we have

(24) \begin{align}{\textbf{B}}({\textbf{x}}){{\textbf{R}}^{ - 1}}({\textbf{x}}){{\textbf{B}}^T}({\textbf{x}}) = {{\textbf{A}}_{{\textrm{cl}}}}({\textbf{x}}){\textbf{E}}({\textbf{x}}) + {\textbf{E}}({\textbf{x}}){\textbf{A}}_{{\textrm{cl}}}^T({\textbf{x}}).\end{align}

The algebraic Lyapunov Eq. (22) results in ${\dot{\textbf{E}}}({\textbf{x}}) = \textbf{0}$ . Regarding frozen computation at each simulation time-step, we neglect the time derivative of ${{\textbf{A}}_{{\textrm{cl}}}}({\textbf{x}})$ and as a result, ${{{\textrm{d}}\left( {{{\textbf{A}}_{{\textrm{cl}}}}({\textbf{x}})\left( {t - {t_{\textrm{f}}}} \right)} \right)} \over {{\textrm{d}}t}} = {{\textbf{A}}_{{\textrm{cl}}}}({\textbf{x}}) + \underbrace {{\textrm{ }}{{{\dot{\textbf{A}}}}_{{\textrm{cl}}}}({\textbf{x}})\left( {t - {t_{\textrm{f}}}} \right)}_{\textbf{0}}$ . Substituting (24) into (23), mathematical manipulation cancels all terms and shows that the solution (21) holds for Eq. (20). Since we neglected ${{\dot{\textbf{A}}}_{{\textrm{cl}}}}({\textbf{x}})\left( {t - {t_{\textrm{f}}}} \right) \approx \textbf{0}$ in the derivative $\frac{{{\textrm{d}}\left( {{{\textbf{A}}_{{\textrm{cl}}}}({\textbf{x}})\left( {t - {t_{\textrm{f}}}} \right)} \right)}}{{{\textrm{d}}t}}$ and considered ${\dot{\textbf{E}}}({\textbf{x}}) = \textbf{0}$ , this approach is so-called a closed-form approximate solution.

The positive gain of the SDDRE (15) is regarded as ${\textbf{K}}\left( {{\textbf{x}}(t)} \right) = {{\textbf{K}}_{{\textrm{ss}}}}\left( {{\textbf{x}}(t)} \right) + {{\textbf{P}}^{ - 1}}\left( {{\textbf{x}}(t)} \right),$ in which ${{\textbf{K}}_{{\textrm{ss}}}}\left( {{\textbf{x}}(t)} \right)$ could be a negative definite ${\textbf{K}}_{{\textrm{ss}}}^ - \left( {{\textbf{x}}(t)} \right)$ or positive definite ${\textbf{K}}_{{\textrm{ss}}}^ + \left( {{\textbf{x}}(t)} \right)$ solution to the SDRE (16). It works with both of them.

The details of the positive and negative roots of the SDRE are reported in Sections 3.1.1 and 3.3.2 of Ref. [50]. The negative root is computationally more robust than the positive one, and it has been used here in this work. Finally, notice that the negative definite solution to the SDRE (16) is the positive definite answer to

\begin{align*}- {\textbf{K}}_{\textrm{n}}^ + ({\textbf{x}}){\textbf{A}}({\textbf{x}}) - {{\textbf{A}}^T}({\textbf{x}}){\textbf{K}}_{\textrm{n}}^ + ({\textbf{x}}) - {\textbf{K}}_{\textrm{n}}^ + ({\textbf{x}}){\textbf{B}}({\textbf{x}}){{\textbf{R}}^{ - 1}}({\textbf{x}}){{\textbf{B}}^T}({\textbf{x}}){\textbf{K}}_{\textrm{n}}^ + ({\textbf{x}}) + {\textbf{Q}}({\textbf{x}}) = \textbf{0},\end{align*}

with consideration of ${\textbf{K}}_{{\textrm{ss}}}^ - \left( {{\textbf{x}}(t)} \right) = - {\textbf{K}}_{\textrm{n}}^ + \left( {{\textbf{x}}(t)} \right)$ .

6. Implementation of the SDDRE on quaternion-based dynamics

The state-space equation of the system (9) must be represented by SDC forms (11) and (12). On the one hand, the dimension of the state vector (13 states) is not compatible with the cascade approach, commonly used in quadrotor control (12 states) [5]. On the other hand, the separation of rotational and translational dynamics was reported helpful in the control design due to different speeds of them, slow and fast dynamic [54]. So, here we propose two subcontrollers for the translation and rotational dynamics, connected through cascade design. For the translation dynamic, the set of SDC parameterization is

\begin{align*}{{\textbf{A}}_\textrm{t}}( {{\textbf{x}}(t)}) = \left[ {\begin{array}{c@{\quad}c}{{\textbf{0}_{3 \times 3}}} {}& {{\textbf{R}}({{\textbf{q}}(t)})}\\ \\[-7pt] {{\textbf{0}_{3 \times 3}}}& {}{ - \dfrac{{\textbf{D}}}{m}}\end{array}} \right]\!,{{\textbf{B}}_\textrm{t}} = \left[ {\begin{array}{*{20}{c}}{{\textbf{0}_{3 \times 3}}}\\ \\[-7pt] {\dfrac{{{{\textbf{I}}_{3 \times 3}}}}{m}}\end{array}} \right]\!,\end{align*}

in which hovering condition has been assumed to find ${\dot {\boldsymbol\xi} _1} = \underbrace {{\textbf{R}}({\textbf{q}})}_{\textbf{I}}{\boldsymbol\upsilon _1} \approx {\boldsymbol\upsilon _1}$ , to make the factorization possible, and “t” stands for translation. For the orientation section, we neglect the scalar part of the quaternion, the first column, and the row of ${\textbf{Q}}({{\textbf{q}}(t)})$ in (9), then the SDC parameterization is

\begin{align*}{{\textbf{A}}_{\textrm{o}}}( {{\textbf{x}}(t)}) = \left[ {\begin{array}{c@{\quad}c}{{\textbf{0}_{3 \times 3}}}& {{{\left[ {{\textbf{Q}}\left( {2:\ 4,2:\ 4} \right)} \right]}_{3 \times 3}}}\\ \\[-7pt] {{\textbf{0}_{3 \times 3}}} &{ - {{\textbf{J}}^{ - 1}}{{\left[ {{{{\boldsymbol\upsilon }}_2} \times {\textbf{I}}} \right]}_{3 \times 3}}}\end{array}} \right]\!,{{\textbf{B}}_{\textrm{o}}} = \left[ \begin{array}{c}{{\textbf{0}_{3 \times 3}}}\\ \\[-7pt] {{{\textbf{J}}^{ - 1}}}\end{array} \right]\!,\end{align*}

where ${\textbf{Q}}\left( {2:\ 4,2:\ 4} \right)$ collects the second to fourth components of ${\textbf{Q}}({{\textbf{q}}(t)})$ in columns and rows, and “o” stands for orientation. Based on (14), the translation control law for regulation to the desired condition rather than zero is

(25) \begin{align}{{\textbf{u}}_\textrm{t}}(t) = - {\textbf{R}}_\textrm{t}^{ - 1}\left( {{\textbf{x}}(t)} \right){\textbf{B}}_\textrm{t}^T{{\textbf{K}}_\textrm{t}}\left( {{\textbf{x}}(t)} \right)\left[ {\begin{array}{*{20}{c}}{{{{\boldsymbol\xi }}_1}(t) - {{{\boldsymbol\xi }}_{1,{\textrm{des}}}}}\\ \\[-7pt] {{{{\dot{\boldsymbol{\xi}}}}_1}(t) - {{{\dot{\boldsymbol{\xi}}}}_{1,{\textrm{des}}}}}\end{array}} \right]\!,\end{align}

where ${{\textbf{R}}_{}}$ , ${{\textbf{Q}}_\textrm{t}}$ , ${{\textbf{F}}_\textrm{t}}$ , and ${{\textbf{K}}_\textrm{t}}$ are the corresponding matrices (with proper dimension) for the translation and Riccati equation and “des” defines the desired condition.

The computation of input law (25) is based on a fully actuated system, which is not possible in reality. So, to transform the ${\left[ {{{\textbf{u}}_{\textrm{t}}}(t)} \right]_{3 \times 1}}$ to the scalar thrust input considering the gravity as well, we use cascade design [55]:

\begin{align*}{T_{\textrm{B}}}(t) = m\left( {{{\left[ {{{\textbf{R}}_3}({{\textbf{q}}(t)})} \right]}_1}{u_{\textrm{t},1}}(t) + {{\left[ {{{\textbf{R}}_3}({{\textbf{q}}(t)})} \right]}_2}{u_{\textrm{t},2}}(t) + {{\left[ {{{\textbf{R}}_3}({{\textbf{q}}(t)})} \right]}_3}\left( {{u_{\textrm{t},3}}(t) + g} \right)} \right),\end{align*}

where that is ${\left[ {{{\textbf{R}}_3}({{\textbf{q}}(t)})} \right]_1}$ is the first component of ${{\textbf{R}}_3}({{\textbf{q}}(t)})$ . The cascade design delivers the necessary desired Euler angles

(26) \begin{align}{\theta _{{\textrm{des}}}}(t) = {\textrm{ta}}{{\textrm{n}}^{ - 1}}\left( {\frac{{{u_{\textrm{t},1}}{\textrm{cos}}{\psi _{{\textrm{des}}}} + {u_{\textrm{t},2}}{\textrm{sin}}{\psi _{{\textrm{des}}}}}}{{{u_{\textrm{t},3}} + g}}} \right),\end{align}

(27) \begin{align}{\phi _{{\textrm{des}}}}(t) = {\textrm{si}}{{\textrm{n}}^{ - 1}}\left( {\frac{{{u_{\textrm{t},1}}{\textrm{sin}}{\psi _{{\textrm{des}}}} - {u_{\textrm{t},2}}{\textrm{cos}}{\psi _{{\textrm{des}}}}}}{{\sqrt {u_{\textrm{t},1}^2 + u_{\textrm{t},2}^2 + {{\left( {{u_{\textrm{t},3}} + g} \right)}^2}} }}} \right),\\ \nonumber\end{align}

in which ${\psi _{{\textrm{des}}}}$ could possess an arbitrary value.

The orientation control law is also presented as

(28) \begin{align}{{\textbf{u}}_{\textrm{o}}}(t) = - {\textbf{R}}_{\textrm{o}}^{ - 1}\left( {{\textbf{x}}(t)} \right){\textbf{B}}_{\textrm{o}}^T{{\textbf{K}}_{\textrm{o}}}\left( {{\textbf{x}}(t)} \right)\left[ {\begin{array}{*{20}{c}}{\mathcal{F}{{\textbf{e}}_{\textbf{q}}}(t)}\\ \\[-7pt] {{{{\boldsymbol\upsilon }}_2}(t) - {{{\boldsymbol\upsilon }}_{2,{\textrm{des}}}}}\end{array}} \right]\!,\end{align}

where ${{\textbf{e}}_{\textbf{q}}}(t)$ is quaternion error and

\begin{align*}\mathcal{F} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}0 &0& 0& 0\\ \\[-7pt]0& 1 &0& 0\\ \\[-7pt] 0& 0& 1 &0\\ \\[-7pt] 0 &0 &0& 1\end{array}} \right].\end{align*}

The quaternion error in (28) is defined as

\begin{align*}{{\textbf{e}}_{\textbf{q}}}(t) = - {{\textbf{q}}_{{\textrm{des}}}}(t) \otimes {{\textbf{q}}^*}(t) = - {\textbf{Q}}\left( {{{\textbf{q}}_{{\textrm{des}}}}(t)} \right){{\textbf{q}}^*}(t),\end{align*}

where ${{\textbf{q}}_{{\textrm{des}}}}(t)$ is defined by substituting (26), (27), and ${\psi _{{\textrm{des}}}}$ into Eq. (5).

7. Simulations

7.1. Validation

To validate the proposed controller design and the simulator, a comparison has been done with the PD controller and a conventional sliding mode control (SMC). We consider a system based on the Euler angles in rotational dynamics, controlled by a simple PD to have it as a reference. Then the SDDRE controller is implemented on the quaternion-based dynamics to check the performance and also to validate the correctness of the implementation. The SMC has been frequently used in UAV control [56], in the context of robustness or combined with other techniques [57]. Based on that, the SMC is selected for comparison to add the more detailed result.

The mass of the system is $m = 0.468\left( {{\textrm{kg}}} \right)$ , the drag coefficient matrix is ${\textbf{D}} = {\textrm{diag}}\left( {0.25,0.25,0.25} \right)\left( {\frac{{{\textrm{kg}}}}{{\textrm{s}}}} \right)$ , the components of the inertia matrix are ${I_{xx}} = 4.856 \times {10^{ - 3}}\left( {{\textrm{kg}}{{\textrm{m}}^2}} \right)$ , ${I_{yy}} = 4.856 \times {10^{ - 3}}\left( {{\textrm{kg}}{{\textrm{m}}^2}} \right)$ , and ${I_{zz}} = 8.801 \times {10^{ - 3}}\left( {{\textrm{kg}}{{\textrm{m}}^2}} \right)$ . The time of the simulation has been set ${t_\textrm{f}} = 10\left( {\textrm{s}} \right)$ , and the drone flies from zero coordinate to the desired position in Cartesian coordinate $\left( { - 2,3,1.5} \right)\left( {\textrm{m}} \right)$ . All of the initial conditions are set zero including position, velocity, Euler angles, and angular velocity. The initial condition of the quaternions is found by substituting $\left( {\phi \left( 0 \right),\theta \left( 0 \right),\psi \left( 0 \right)} \right)$ into (5) to reach ${\textbf{q}}\left( 0 \right) = {\left[ {1,0,0,0} \right]^T}$ .

The PD control gains are set as ${{\textbf{K}}_{\textrm{P,t}}} = {{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{K}}_{\textrm{D,t}}} = 2 \times {{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{K}}_{\textrm{P,o}}} = {{\textbf{I}}_{3 \times 3}}$ , and ${{\textbf{K}}_{\textrm{D,o}}} = 0.5 \times {{\textbf{I}}_{3 \times 3}}$ . The SDDRE controller gains are selected as follows: ${{\textbf{R}}_\textrm{t}} = {{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{Q}}_\textrm{t}} = {\textrm{diag}}\left( {0.1,0.1,0.1,0.2,0.2,0.2} \right)$ , ${{\textbf{F}}_\textrm{t}} = 100 \times {{\textbf{Q}}_\textrm{t}}$ , ${{\textbf{R}}_{\textrm{o}}} = {{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{Q}}_{\textrm{o}}} = {\textrm{diag}}\left( {2,2,2,0,0,0} \right)$ , ${{\textbf{F}}_\textrm{o}} = 10 \times {{\textbf{Q}}_\textrm{o}}$ .

Two separate SMC controllers are considered for translation and orientation parts, consistent with Section 6. The sliding surface is $\textbf{s}_{i}(\textbf{X})=\dot{\tilde{\textbf{X}}}_{i}+\boldsymbol\Lambda_{i}\tilde{\textbf{X}}_{i}$ for $i = \left\{ {\textrm{t,o}} \right\}$ (translation and orientation) where $\tilde{\textbf{X}}_{i}=\boldsymbol\xi_{i}-\boldsymbol\xi_{i,\textrm{des}}$ and ${{{\boldsymbol\Lambda }}_i}$ is a strictly positive constant matrix. The control law is also in the form of ${{\textbf{u}}_i} = {\textbf{B}}_{i,{\textrm{SMC}}}^{ - 1}({\textbf{x}})\left( {{{{\ddot{\boldsymbol\xi}}}_{i,{\textrm{des}}}} - {{\textbf{f}}_{i,{\textrm{SMC}}}}({\textbf{x}}) - {{\textbf{K}}_{i,{\textrm{SMC}}}}\tanh \left( {\frac{{{{\textbf{s}}_i}({\textbf{x}})}}{\sigma }} \right)} \right)$ , where ${{\textbf{B}}_{\textrm{t},{\textrm{SMC}}}} = 1/m{{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{B}}_{\textrm{o},{\textrm{SMC}}}}({\textbf{x}}) = {\textbf{J}}_\textrm{c}^{ - 1}\left( {{{{\boldsymbol\xi }}_2}} \right)$ , ${{\textbf{f}}_{\textrm{t},{\textrm{SMC}}}}({\textbf{x}}) = - \frac{1}{m}{\textbf{D}}{{\dot{\boldsymbol{\xi}}}_1}(t)$ , ${{\textbf{f}}_{\textrm{o},{\textrm{SMC}}}}({\textbf{x}}) = - {{\textbf{J}}^{ - 1}}\left( {{{{\boldsymbol\xi }}_2}} \right){\textbf{C}}\left( {{{{\boldsymbol\xi }}_2},{{{\dot{\boldsymbol{\xi}}}}_2}} \right){{\dot{\boldsymbol{\xi}}}_2}$ , and ${{\textbf{K}}_{i,{\textrm{SMC}}}}$ is the correction gain of the SMC. More details on conventional dynamics, such as ${{\textbf{J}}_\textrm{c}}\left( {{{{\boldsymbol\xi }}_2}} \right) = {{\textbf{W}}^T}\left( {{{{\boldsymbol\xi }}_2}} \right){\textbf{IW}}\left( {{{{\boldsymbol\xi }}_2}} \right)$ , can be found in ref. [54]. To avoid chattering in SMC, $\tanh \left( {\frac{{{{\textbf{s}}_i}({\textbf{x}})}}{\sigma }} \right)$ is used where $\sigma = 0.2$ in this simulation. The SMC control parameters are selected as ${{{\boldsymbol\Lambda }}_\textrm{t}} = {\textrm{diag}}\left( {0.5,0.5,1} \right)$ , ${{{\boldsymbol\Lambda }}_\textrm{o}} = 1.5\times{{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{K}}_{\textrm{t},{\textrm{SMC}}}} = {\textrm{diag}}\left( {5,5,1} \right)$ , ${{\textbf{K}}_{\textrm{o},{\textrm{SMC}}}} = 2.5\times{{\textbf{I}}_{3 \times 3}}$ , and desired accelerations are set zero, ${{\ddot{\boldsymbol\xi}}_{i,{\textrm{des}}}} = \textbf{0}$ .

Simulating the system, the results are found in the following. The position variables of the drone are illustrated in Fig. 2 to Fig. 4 with respect. The roll and pitch angles of the multi-copter are plotted in Fig. 5 and Fig. 6 with respect. The input norm of the quadrotor is illustrated in Fig. 7 and the configuration and trajectories of the drones are shown in Fig. 8.

Fig. 2. x-axis regulation of the system, comparison with PD and SMC.

Fig. 3. y-axis regulation of the system, comparison with PD and SMC.

Fig. 4. z-axis regulation of the system, comparison with PD and SMC.

Fig. 5. The roll angle of the system, comparison with PD and SMC.

Fig. 6. The pitch angle of the system, comparison with PD and SMC.

Fig. 7. The input norm of the inputs of the system, comparison with PD and SMC.

Fig. 8. The configuration and trajectories of the quadrotor drones with PD, SMC, and SDDRE controllers.

The error of the regulation with PD controller was gained higher than the other two and the error of the SDDRE was obtained the least, see Table I. Since the norm of the inputs (representative of the energy consumption) of the SDDRE is less than the PD and SMC, the performance of the proposed system is satisfactory. The results also confirm the validity of the quaternion-based dynamics and also the control implementation.

Table I. Comparison of PD, SMC, and SDDRE controller.

Validation with previous work: To confirm the correctness of the quaternion-based dynamics and the SDDRE controller, an existing model will be employed for comparison. Xiong and Zhang used a global fast terminal sliding mode controller (TSMC) for quadrotor regulation and also compared the results with conventional SMC [58]. Here the parameters of the system are substituted into the quaternion model and the SDDRE controls the model. The results are similar to the ones in Fig. 2 of ref. [58] presented in this section, in Fig. 9. The regulation of translation control is quite like the TSMC and regulated to desired values around 2s, without overshoot. The controller parameters are similar to Sections 7.2.

Fig. 9. The validation results of the quaternion-based SDDRE with previous work in ref. [58].

It should be noted that since the loop is closed, it cannot be said that the dynamics are validated; however, observing the behavior of the system in comparison with the SMC, the closed-loop system is validated.

7.2. Cobra maneuver

The motivation of using quaternions and avoiding Euler angles in rotational dynamics are to gain a singularity-free controller, and as a result, obtain agile flight and aerobatic maneuvers. One of the hardest positions in quadrotor control in $\pi /2\left( {{\textrm{rad}}} \right)$ rotation either in roll or pitch angles. The Cobra maneuver is famous for a jet aircraft to perform aerobatic shows or in combat for sudden brake, etc. For the quadrotors, it is the first time that this motion is simulated in a forward flight; the Cobra in ascending motion was reported [44]. That is a challenge since, in $\phi = \pi /2$ or $\theta = \pi /2$ , there is no thrust to compensate the gravity; for an aircraft, a jet engine supports the gravity. This might cause a crash or fall for the multirotor. To perform the Cobra maneuver, ${\theta _{{\textrm{des}}}} = \pi /2$ is imposed in $t \in \left[ {1,1.35} \right]\left( {\textrm{s}} \right)$ instead of Eq. (26). After that, the multirotor drone tries to recover the stability and regulate to final condition. The simulation time is ${t_\textbf{f}} = 10\left( {\textrm{s}} \right),$ and the parameters of the control are ${{\textbf{R}}_\textbf{t}} = {{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{Q}}_\textbf{t}} = {\textrm{diag}}\left( {1,1,1,0.5,0.5,0.5} \right)$ , ${{\textbf{F}}_\textbf{t}} = 10 \times {{\textbf{Q}}_{}}$ , ${{\textbf{R}}_\textbf{o}} = {{\textbf{I}}_{3 \times 3}}$ , ${{\textbf{Q}}_\textbf{o}} = {\textrm{diag}}\left( {2,2,2,0,0,0} \right)$ , ${{\textbf{F}}_\textbf{o}} = 10 \times {{\textbf{Q}}_\textbf{o}}$ . The start point of the regulation is set at zero along with other initial conditions; the endpoint is chosen $\left( {5,1,1.25} \right)\left( {\textrm{m}} \right)$ . Simulating the drone, the error is found $7.75\left( {{\textrm{mm}}} \right)$ and the system successfully performed the maneuver. The position variables and attitude ones are shown in Fig. 10 and Fig. 11 with respect. The trajectory and configuration of the system are demonstrated in Fig. 12. The quaternions are plotted in Fig. 13. The input thrust and input torque signals are presented in Fig. 14 and Fig. 15 with respect.

Fig. 10. The position variables of the system in aerobatic maneuver.

Fig. 11. The Euler angles of the drone in aerobatic maneuver.

Fig. 12. The Cobra maneuver via the SDDRE and quaternion dynamics.

Fig. 13. The quaternions.

Fig. 14. The input thrust of the quadrotor.

Fig. 15. The input torque signals of the drone.

8. Conclusions

This work investigated the quaternion-based control design using the SDDRE to control a quadrotor in aerobatic flight. The Euler angles are vulnerable to big changes in attitude and rotational dynamics, specifically, the rotation matrix. The diagonal components of ${\textbf{R}}\left( {\phi ,\theta ,\psi } \right)$ become zero for either of $\phi ,\theta ,\psi $ at $\pi /2$ . This means the omission of thrust in flight and unstable conditions. Specifically, the controllability pair will be unsatisfied. To solve this problem, quaternion representation has been used to have a singular-free attitude control and rotation matrix. The introduced model has been validated through a comparison with the conventional Euler dynamics controlled by a PD input law. To show the application of the proposed method, a challenging maneuver, Cobra, has been simulated in the forward flight, successfully controlled. The Cobra maneuver has put the drone in a position without thrust to compensate for the gravity; however, this approach generated a stable motion to reduce the fall in that condition.