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. [Reference Kang, Sadegh and Urschel20]. 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 [Reference Fresk and Nikolakopoulos25]:

\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 [Reference Diebel45]:

(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 [Reference Fresk and Nikolakopoulos25]:

(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 [Reference Nekoo46]:

(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) [Reference Nekoo47].

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. [Reference Nekoo46]), then the SDDRE is represented as [Reference Nekoo48]:

(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) [Reference Jayaram and Tadi49]. To solve the SDDRE (15) and to find the control gain, several methods could be used such as backward integration (BI) [Reference Korayem and Nekoo50], state-transition matrix [Reference Mathavaraj and Padhi51], and Lyapunov-based approach [Reference Heydari and Balakrishnan52]. 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 [Reference Korayem and Nekoo53], 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))$ [Reference Korayem and Nekoo50]. 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. [Reference Korayem and Nekoo50]:

(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. [Reference Korayem and Nekoo50]. 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) [Reference Nekoo, Acosta and Ollero5]. 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 [Reference Nekoo, Acosta, Gomez-Tamm and Ollero54]. 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 [Reference Nekoo, Acosta and Ollero55]:

\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 [Reference Xu and Ozguner56], in the context of robustness or combined with other techniques [Reference Abro, Zulkifli, Asirvadam and Ali57]. 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. [Reference Nekoo, Acosta, Gomez-Tamm and Ollero54]. 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 [Reference Xiong and Zhang58]. 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. [Reference Xiong and Zhang58] 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. [Reference Xiong and Zhang58].

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 [Reference Chen and Pérez-Arancibia44]. 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.