Backstepping-Based Inverse Optimal Attitude Control of Quadrotor

Abstract

Abstract Input saturation must be taken into account for applying rapid reorientation in the large angle manoeuvre of a quadrotor. In this paper, a backstepping-based inverse optimal attitude controller (BIOAC) is derived which has the property of a maximum convergence rate in the sense of a control Lyapunov function (CLF) under input torque limitation. In the controller, a backstepping technique is used for handling the complexity introducing by the unit quaternion representation of the attitude of a quadrotor with four parameters. Moreover, the inverse optimal approach is employed to circumvent the difficulty of solving the Hamilton-Jacobi-Bellman (HJB) equation. The performance of BIOAC is compared with a PD controller in which the input torque limitation is not considered under the same unit quaternion representation using numerical simulation while the results show that BIOAC gains faster convergence with less control effort. Next, BIOAC is realized on a test bed and the effectiveness of the control law is verified by experimental studies.

1. Introduction

Compared with a fixed-wing aircraft, a rotorcraft has special advantages in city and indoor environments due to a capability for VTOL (vertical take-off and landing) and hovering flight. As a new kind of rotorcraft, a quadrotor has certain unique advantages over other traditional rotorcrafts. Firstly, it avoids the complexity of swash plate and mechanical linkages. Secondly, it is much safer for indoor flying because its propellers are smaller and can be wrapped up. In addition, it has a greater thrust-weight ratio and better manoeuvrability performance. Therefore, a great deal of work has been done to research the control of the quadrotor. The projects in [1] and [2] use commercially-available toys to derive attitude controllers and position controllers respectively. Some nonlinear control methods, such as backstepping [3], sliding mode [4] and inverse dynamics [5] are used for better control performance of the quadrotor. The OS4 project [6] and X-4 [7] develop their own prototypes for generating more thrust and better stability. Most of the works presented in the literature use simplified models, where nonlinear effects and the performance of the actuators are ignored. However, when a rapid reorientation problem is considered for a large angle manoeuvre of a quadrotor, input saturation will cause control failure and must be taken into account. This is one of the unavoidable problems in the actual control system [8].

In this paper, the maximum rate attitude control problem under the input saturation is addressed. Moreover, a backstepping-based inverse optimal attitude controller (BIOAC) is derived which has the property of a maximum convergence rate within the meaning of a control Lyapunov function under input torque limitation. For designing BIOAC, the representation of the attitude of the quadrotor must first be determined. In all representations, the unit quaternion is a pervasive one. Distinct from minimal attitude representations, such as Euler angles - which have the disadvantages of singularities and complicated trigonometric expressions - the unit quaternion provides a globally non-singular representation of attitude. Compared with a classical direction cosine matrix, it has only 4 numbers and the multiplication and propagation rule can be calculated more effectively [9–11]. However, unit quaternion representation for the attitude of a quadrotor introduces some difficulties in controller design. For simplifying the design procedure of the controller, a backstepping technique [12,13] is employed in this work.

The optimal control problem for nonlinear systems is usually attributed to the solvability of the Hamilton-Jacobi-Bellman(HJB) equation. Moreover, it is particularly difficult to solve the HJB equation in the optimal control problem while taking into account the attitude kinematics and dynamics at the same time. An alternative method is the so-called inverse optimal control method in designing an optimal feedback controller [14–17]. Based on a control Lyapunov function (CLF), the inverse optimal approach circumvents solving the HJB equation. In addition, a feedback controller for a set of reasonable cost functions can be achieved.

This article is organized as follows. In section two, a mathematical model of a quadrotor on the attitude platform is given based upon the assumptions of a single rigid body and a symmetrical structure. Next, a maximum rate attitude controller is designed under input saturation conditions in section three. The performance of the proposed BIOAC is compared with a PD controller using simulations in section four, which is followed by experimental studies in section five. Finally, the conclusion is outlined in the last section.

2. Attitude Dynamics of the Quadrotor

2.1 Unit quaternion kinematics

A unit quaternion $\stackrel{\rightharpoonup }{q}$ is a complex number with four parameters. It has the norm of constant one, and can be seen as a point of the unit sphere S³

$\begin{array}{l}\stackrel{\rightharpoonup }{q}={\mathrm{q}}_0+{\mathrm{q}}_1\mathrm{i}+{\mathrm{q}}_2\mathrm{j}+{\mathrm{q}}_3\mathrm{k}\\ {\mathrm{q}}_0^2+{\mathrm{q}}_1^2+{\mathrm{q}}_2^2+{\mathrm{q}}_3^2=1\end{array}$

(1)

Let q represent its vector part, i.e. q˜ = q₁i + q₂j + q₃k, then the unit quaternion can be written in the following form:

$\stackrel{\rightharpoonup }{q}=\left[{\mathrm{q}}_0,\tilde{q}\right]$

(2)

Given a unit quaternion, there is a correspondence between the unit quaternion and rotation transformation:

$\mathrm{R}\left({\mathrm{q}}_0,\tilde{q}\right)=\left({\mathrm{q}}_0^2-{\tilde{q}}^{\mathrm{T}}\tilde{q}\right)\mathrm{I}+2{\mathrm{q}}_0\mathrm{S}\left(\tilde{q}\right)+2\tilde{q}{\tilde{q}}^{\mathrm{T}}$

(3)

where:

$\mathrm{S}\left(\tilde{q}\right)=\left[\begin{array}{ccc}0& -{\mathrm{q}}_3& {\mathrm{q}}_2\\ {\mathrm{q}}_3& 0& -{\mathrm{q}}_1\\ -{\mathrm{q}}_2& {\mathrm{q}}_1& 0\end{array}\right]$

(4)

Notice R(q₀,q˜) = R(–q₀,–q˜). Thus, the unit sphere S³ represented by the unit quaternion is a double cover of the attitude space SO(3). The phenomenon is clearer using the angle form of the unit quaternion. Euler's theorem states that any rotation in space can be achieved by rotating an angle θ around an eigenaxis $\stackrel{\rightharpoonup }{n}$ ; then, an unit quaternion can be written as:

$\stackrel{\rightharpoonup }{q}=\left(\cos (\theta /2),\sin (\theta /2)\stackrel{\rightharpoonup }{n}\right)\begin{array}{cc}& \theta \in (-2\pi ,2\pi ]\end{array}$

(5)

It is obvious that rotating the angle 0 < θ₀ ≤ 2π around eigenaxis $\stackrel{\rightharpoonup }{n}$ is equivalent to a rotation of the 2π – θ₀ angle around an eigenaxis of $-\stackrel{\rightharpoonup }{n}$ . In Fig. 1, the X point and O point are different by quaternion representation, but have the same physical meaning. Actually, $\stackrel{\rightharpoonup }{q}=\left(\cos ({\theta }_0/2),\sin ({\theta }_0/2)\stackrel{\rightharpoonup }{n}\right)$ and $-\stackrel{\rightharpoonup }{q}=\left(\cos (\pi +{\theta }_0/2),\sin (\pi +{\theta }_0/2)\stackrel{\rightharpoonup }{n}\right)$ have the same physical meaning.

Figure 1.

Quaternion representation along some fixed axis

The conjugate of a unit quaternion is:

${\stackrel{\rightharpoonup }{q}}^{\ast }=\left({\mathrm{q}}_0,-\tilde{q}\right)$

(6)

Given two unit quaternion ${\stackrel{\rightharpoonup }{q}}_1=\left({\mathrm{q}}_{01},{\tilde{q}}_1\right)$ and ${\stackrel{\rightharpoonup }{q}}_2=\left({\mathrm{q}}_{02},{\tilde{q}}_2\right)$ , the multiplication ⊗ is defined as:

${\stackrel{\rightharpoonup }{q}}_1\otimes {\stackrel{\rightharpoonup }{q}}_2=\left[\begin{array}{c}{\mathrm{q}}_{01}{\mathrm{q}}_{02}-{\tilde{q}}_1^{\mathrm{T}}{\tilde{q}}_2\\ {\mathrm{q}}_{01}{\tilde{q}}_2+{\mathrm{q}}_{02}{\tilde{q}}_1+\mathrm{S}\left({\tilde{q}}_1\right){\tilde{q}}_2\end{array}\right]$

(7)

The unit quaternion multiplication is equivalent to the multiplication of the rotation matrix. If ${\stackrel{\rightharpoonup }{q}}_1$ and ${\stackrel{\rightharpoonup }{q}}_2$ correspond to the rotation matrices R₁ and R₂ in SO(3), then there are the following correspondences:

$\begin{array}{l}{\mathrm{R}}_1{\mathrm{R}}_2\iff \pm {\stackrel{\rightharpoonup }{q}}_1\otimes {\stackrel{\rightharpoonup }{q}}_2\\ {\mathrm{R}}_1\stackrel{\rightharpoonup }{v}\iff {\stackrel{\rightharpoonup }{q}}_1\otimes \bar{v}\otimes {\stackrel{\rightharpoonup }{q}}_1^{\ast }\end{array}$

(8)

where $\bar{v}={\left(0,{\stackrel{\rightharpoonup }{v}}^{\mathrm{T}}\right)}^{\mathrm{T}}$ .

When a rigid body is rotating with an angular velocity vector $\stackrel{\rightharpoonup }{w}$ , the time derivative of the unit quaternion can be given by the unit quaternion propagation rule:

$\stackrel{.}{\stackrel{\rightharpoonup }{q}}=\frac{1}{2}\stackrel{\rightharpoonup }{q}\otimes \bar{w}$

(9)

2.2 System model

In order to establish the dynamic equations of the quadrotor attitude system, coordinate frames must be establish first of all. This article uses two coordinate frames: an inertial coordinate frame and a body coordinate frame. As seen in Fig. 2, the origin of the body coordinate frame O_b – X_bY_bZ_b is located in the centre of the cross beam of the quadrotor: axis X pointing towards propeller 1, axis Z is perpendicular to the cross beam plane, and the three axis is to form a right-handed coordinate system.

Figure 2.

Coordinate frames of the quadrotor

The quadrotor dynamics are based upon two assumptions:

The quadrotor is a rigid body.

Any small error in the manufacturing and installation of the vehicle body is not considered, so the quadrotor is assumed to be symmetric.

Based upon assumption one, the following mathematical model can be derived from the Newton-Euler equations and quaternion kinematics:

$J\stackrel{.}{\stackrel{\rightharpoonup }{w}}=\stackrel{\rightharpoonup }{u}-\stackrel{\rightharpoonup }{w}\times (J\stackrel{\rightharpoonup }{w})$

(10a)

$\{\begin{array}{c}\stackrel{.}{\tilde{q}}=\frac{1}{2}\left({\mathrm{q}}_0\mathrm{I}+\mathrm{S}\left(\tilde{q}\right)\right)\overrightarrow{w}\\ {\stackrel{.}{q}}_0=-\frac{1}{2}{\tilde{q}}^{\mathrm{T}}\overrightarrow{w}\end{array}$

(10b)

where $\stackrel{\rightharpoonup }{q}={\left[{\mathrm{q}}_0 {\tilde{q}}^{\mathrm{T}}\right]}^{\mathrm{T}}$ represents the attitude quaternion of a body coordinate frame with respect to an inertial coordinate frame, $\stackrel{\rightharpoonup }{w}$ is the angular velocity in body coordinate frame, and $\stackrel{\rightharpoonup }{u}$ is the torque acting on the quadrotor.

According to the second assumption, the correspondence between four propeller thrusts f_i (i = 1,2,3,4) and the torque u(i = 1,2,3) acting on the quadrotor is given as follows:

$\left[\begin{array}{c}\mathrm{f}\\ {\mathrm{u}}_1\\ {\mathrm{u}}_2\\ {\mathrm{u}}_3\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& \mathrm{d}& 0& -\mathrm{d}\\ -\mathrm{d}& 0& \mathrm{d}& 0\\ -{\mathrm{c}}_{\tau \mathrm{f}}& {\mathrm{c}}_{\tau \mathrm{f}}& -{\mathrm{c}}_{\tau \mathrm{f}}& {\mathrm{c}}_{\tau \mathrm{f}}\end{array}\right]\left[\begin{array}{c}{\mathrm{f}}_1\\ {\mathrm{f}}_2\\ {\mathrm{f}}_3\\ {\mathrm{f}}_4\end{array}\right]$

(11)

where f is the total thrust and d is the distance between the propeller and the centre of the cross beam and c_τf is the torque thrust ratio of the propeller.

Our model does not include the propeller gyroscopic torque and the flapping dynamic term; thus, $\stackrel{\rightharpoonup }{u}$ is consists of just three axis control torques u₁,u₂,u₃:

$\stackrel{\rightharpoonup }{u}=\left[\begin{array}{c}{\mathrm{u}}_1\\ {\mathrm{u}}_2\\ {\mathrm{u}}_3\end{array}\right]$

(12)

Considering the low-speed characteristics of the quadrotor, the following simplified equation is used to derive the propeller thrust [1]:

${\mathrm{f}}_{\mathrm{i}}={\mathrm{c}}_{\mathrm{f}}{\mathrm{\Omega }}_{\mathrm{i}}^2$

(13)

where c_f is a constant coefficient.

3. Controller Design

The objective is to design a time optimal attitude controller under the limitations of input torque. In order to circumvent the difficulty of solving the Hamilton-Jacobi-Bellman (HJB) equation, the inverse optimal control approach is used. However, since the unit quaternion has four elements, it is difficult to implement the inverse optimal theorem directly. In order to simplify the design complexity, the backstepping technique is employed in the controller design procedure.

3.1 Inverse Optimal Theorem

The inverse optimal control approach is originated by Kalman to establish certain gain and phase margins of linear quadratic regulators. However, until Freeman uses the approach in [15] to develop a methodology for designing a robust nonlinear controller, it has been long dormant. In this section, only the main results of the approach are given, more detailed material and proofs can be found in [14–16].

Consider a control affine nonlinear system:

$\stackrel{.}{x}=\mathrm{f}\left(\mathrm{x}\right)+\mathrm{g}\left(\mathrm{x}\right)\mathrm{u}$

(14)

where u ∈ ℝⁿ is a control input, and f(x) and g(x) are continuous function matrices.

Definition 1: control Lyapunov function(CLF)

Let V be a C¹, proper, positive definite function V: ℝⁿ → ℝ₊, such that:

$\underset{\mathrm{u}}{\inf }\left[{\mathrm{L}}_{\mathrm{f}}\mathrm{V}\left(\mathrm{x}\right)+{\mathrm{L}}_{\mathrm{g}}\mathrm{V}\left(\mathrm{x}\right)\mathrm{u}\right]<0$

(15)

for all x ≠ 0. The existence of a CLF for the system (14) is equivalent to the existence of a globally asymptotically stabilizing control law u = k(x).

Theorem 1: inverse optimal theorem

Assume that the static state feedback control law:

$\mathrm{u}=\kappa \left(\mathrm{x}\right):=-{\mathrm{N}}^{-1}\left(\mathrm{x}\right){\left({\mathrm{L}}_{\mathrm{g}}\mathrm{V}\right)}^{\mathrm{T}}$

(16)

stabilizes the system in (14) with respect to a positive definite radially unbounded Lyapunov function V (x), where N: ℝⁿ → ℝⁿ is a positive definite matrix-valued function. Then the control law:

$\mathrm{u}={\kappa }^{\ast }\left(\mathrm{x}\right):=\beta \kappa \left(\mathrm{x}\right)\begin{array}{cc}& \beta \ge 2\end{array}$

(17)

is optimal with respect to the cost:

$\mathrm{J}={\displaystyle {\int }_0^{\infty }\left\{\mathrm{l}\left(\mathrm{x}\right)+{\mathrm{u}}^{\mathrm{T}}\mathrm{R}\left(\mathrm{x}\right)\mathrm{u}\right\}dt}$

(18)

where:

$\begin{array}{l}\mathrm{l}\left(\mathrm{x}\right)=-2\beta {\mathrm{L}}_{\mathrm{f}+\mathrm{g}\kappa }\mathrm{V}+\mathrm{\dots }\\ \mathrm{\dots }\beta \left(\beta -2\right){\mathrm{L}}_{\mathrm{g}}\mathrm{V}\cdot {\mathrm{R}}^{-1}\left(\mathrm{x}\right){\left({\mathrm{L}}_{\mathrm{g}}\mathrm{V}\right)}^{\mathrm{T}}>0\end{array}$

(19)

3.2 Backstepping-based Inverse Optimal Attitude Controller

The controller design procedures of backstepping can be divided into two phases:

Design a virtual control law to stabilize virtual quadrotor kinematics.

Design the true control law which would be implemented on the system based upon the quadrotor dynamics and the virtual control law designed in phase 1.

3.2.1 Control of the Kinematics Subsystem

The kinematics subsystem in (10b) is controlled only indirectly through the angular velocity vector w, so the virtual system is considered:

$\{\begin{array}{c}\stackrel{.}{\tilde{q}}=\frac{1}{2}\left({\mathrm{q}}_0\mathrm{I}+\mathrm{S}\left(\tilde{q}\right)\right)\stackrel{\rightharpoonup }{v}\\ {\stackrel{.}{q}}_0=-\frac{1}{2}{\tilde{q}}^{\mathrm{T}}\stackrel{\rightharpoonup }{v}\end{array}$

(20)

where $\stackrel{\rightharpoonup }{v}$ is the virtual control. Let the attitude error be:

$\begin{array}{c}{\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}={\stackrel{\rightharpoonup }{q}}_{\mathrm{d}}^{\ast }\otimes \stackrel{\rightharpoonup }{q}\hfill \\ =\left[\begin{array}{c}{\mathrm{q}}_{\mathrm{d}0}{\mathrm{q}}_0+{\tilde{q}}_{\mathrm{d}}^{\mathrm{T}}\tilde{q}\\ {\mathrm{q}}_{\mathrm{d}0}\tilde{q}-{\mathrm{q}}_0{\tilde{q}}_{\mathrm{d}}-\mathrm{S}\left({\tilde{q}}_{\mathrm{d}}\right)\tilde{q}\end{array}\right]\end{array}$

(21)

Then the derivative of attitude error is:

$\begin{array}{l}{\stackrel{.}{\stackrel{\rightharpoonup }{q}}}_{\mathrm{e}}=\left[\begin{array}{c}{\mathrm{q}}_{\mathrm{d}0}{\stackrel{.}{\tilde{q}}}_0+{\tilde{q}}_{\mathrm{d}}^{\mathrm{T}}\stackrel{.}{\tilde{q}}\\ {\mathrm{q}}_{\mathrm{d}0}\stackrel{.}{\tilde{q}}-{\stackrel{.}{q}}_0{\tilde{q}}_{\mathrm{d}}-\mathrm{S}\left({\tilde{q}}_{\mathrm{d}}\right)\stackrel{.}{\tilde{q}}\end{array}\right]\\ =\frac{1}{2}\left[\begin{array}{c}-{\mathrm{q}}_{\mathrm{d}0}{\tilde{q}}^{\mathrm{T}}+{\mathrm{q}}_0{\tilde{q}}_{\mathrm{d}}^{\mathrm{T}}+{\left(\mathrm{S}\left({\tilde{q}}_{\mathrm{d}}\right)\tilde{q}\right)}^{\mathrm{T}}\\ \left[{\mathrm{q}}_{\mathrm{d}0}\mathrm{I}-\mathrm{S}\left({\tilde{q}}_{\mathrm{d}}\right)\right]\left[{\mathrm{q}}_0\mathrm{I}+\mathrm{S}\left(\tilde{q}\right)\right]+{\tilde{q}}_{\mathrm{d}}{\tilde{q}}^{\mathrm{T}}\end{array}\right]\stackrel{\rightharpoonup }{v}\\ =\frac{1}{2}\left[\begin{array}{c}-{\tilde{q}}_{\mathrm{e}}^{\mathrm{T}}\\ {\mathrm{q}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left({\tilde{q}}_{\mathrm{e}}\right)\end{array}\right]\stackrel{\rightharpoonup }{v}\end{array}$

(22)

Equation (22) can be rewritten as:

$\{\begin{array}{c}{\stackrel{.}{\tilde{q}}}_{\mathrm{e}}=\frac{1}{2}\left({\mathrm{q}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left({\tilde{q}}_{\mathrm{e}}\right)\right)\stackrel{\rightharpoonup }{v}\\ {\stackrel{.}{q}}_{\mathrm{e}0}=-\frac{1}{2}{\tilde{q}}_{\mathrm{e}}^{\mathrm{T}}\stackrel{\rightharpoonup }{v}\end{array}$

(23)

Consider a Lyapunov function candidate for the virtual system:

${\mathrm{V}}_1={\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]$

(24)

V₁ ≥ 0 and V₁ = 0 if and only if ${\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}={\left[\pm 1 0 0 0\right]}^{\mathrm{T}}$ .

The derivative of V₁ is:

${\stackrel{.}{V}}_1={\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\tilde{q}}_{\mathrm{e}}^{\mathrm{T}}\\ {\mathrm{q}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left({\tilde{q}}_{\mathrm{e}}\right)\end{array}\right]\mathrm{v}$

(25)

Here, sgn(·) is not a traditional sign function, and is defined as:

$\mathrm{sgn}\left(\lambda \right)=\{\begin{array}{c}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\lambda \ge 0\\ -1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\lambda <0\end{array}$

Let:

${\mathrm{G}}^{\mathrm{T}}\left({\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}\right)=\left[\begin{array}{c}\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\tilde{q}}_{\mathrm{e}}^{\mathrm{T}}\\ {\mathrm{q}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left({\tilde{q}}_{\mathrm{e}}\right)\end{array}\right]$

(26)

Then, we have:

$\begin{array}{c}\mathrm{G}\left({\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}\right)\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]={\left[\begin{array}{c}\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\tilde{q}}_{\mathrm{e}}^{\mathrm{T}}\\ {\mathrm{q}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left({\tilde{q}}_{\mathrm{e}}\right)\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]\\ =\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\tilde{q}}_{\mathrm{e}}\end{array}$

(27)

Design the virtue control input as:

$\stackrel{\rightharpoonup }{v}=-{\mathrm{k}}_1\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\tilde{q}}_{\mathrm{e}}$

(28)

with any k₁ > 0. Then:

$\begin{array}{c}{\stackrel{.}{V}}_1={\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]}^{\mathrm{T}}{\mathrm{G}}^{\mathrm{T}}\left({\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}\right)\left(-{\mathrm{k}}_1\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\tilde{q}}_{\mathrm{e}}\right)\\ =-{\mathrm{k}}_1{\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]}^{\mathrm{T}}{\mathrm{G}}^{\mathrm{T}}\left({\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}\right)\mathrm{G}\left({\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}\right)\left[\begin{array}{c}1-\left|{\mathrm{q}}_{\mathrm{e}0}\right|\\ {\tilde{q}}_{\mathrm{e}}\end{array}\right]\end{array}$

(29)

It can be seen that V̇₁ < 0 and V̇₁ = 0 if and only if ${\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}={\left[\pm 1 0 0 0\right]}^{\mathrm{T}}$ . According to the Lyapunov Theorem, the virtual system is asymptotically stable.

After introducing the virtual control $\stackrel{\rightharpoonup }{v}$ , the virtual system (20) becomes a variable structure system - it switches at θ = ±π along any rotating axis. As shown in Fig.3, ${\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}$ converges on [±1 0 0 0]^T with a minimal angle θ_min, where:

Figure 3.

${\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}$ convergence direction

${\theta }_{\min }=\{\begin{array}{c}\begin{array}{cc}\min \left\{\theta ,2\pi -\theta \right\}& 0\le \theta \le 2\pi \end{array}\\ \begin{array}{cc}\min \left\{-\theta ,2\pi +\theta \right\}& -2\pi \le \theta <0\end{array}\end{array}$

It is already known that the unit quaternion $\stackrel{\rightharpoonup }{q}=\left(\cos (\theta /2),\sin (\theta /2)\stackrel{\rightharpoonup }{n}\right)$ and $-\stackrel{\rightharpoonup }{q}=\left(\cos (\pi +\theta /2),\sin (\pi +\theta /2)\stackrel{\rightharpoonup }{n}\right)$ represent the same attitude physically; as such, the differential equation of the attitude error can be formulated as (a similar derivation can be seen in [11]):

$\{\begin{array}{c}\stackrel{.}{\tilde{q}}{\text{'}}_{\mathrm{e}}=\frac{1}{2}\left(\mathrm{q}{\text{'}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left(\tilde{q}{\text{'}}_{\mathrm{e}}\right)\right)\stackrel{\rightharpoonup }{v}\\ \stackrel{.}{q}{\text{'}}_{\mathrm{e}0}=-\frac{1}{2}{\left(\tilde{q}{\text{'}}_{\mathrm{e}}^{}\right)}^{\mathrm{T}}\stackrel{\rightharpoonup }{v}\end{array}$

(30)

where $\stackrel{\rightharpoonup }{q}{\prime }_{\mathrm{e}}=\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\stackrel{\rightharpoonup }{q}}_{\mathrm{e}}$ , with the virtual control:

$\stackrel{\rightharpoonup }{v}=-\mathrm{sgn}\left({\mathrm{q}}_{\mathrm{e}0}\right){\mathrm{k}}_1{\tilde{q}}_{\mathrm{e}}=-{\mathrm{k}}_1\tilde{q}{\text{'}}_{\mathrm{e}}$

(31)

3.2.2 Control of the Full Rigid Body Model

Consider now the angular velocity error:

$\stackrel{\rightharpoonup }{z}=\stackrel{\rightharpoonup }{v}-\stackrel{\rightharpoonup }{w}=-\left(\stackrel{\rightharpoonup }{w}+{\mathrm{k}}_1\tilde{q}{\text{'}}_{\mathrm{e}}\right)$

(32)

Then, the differential Equation (30) can be rewritten as:

$\{\begin{array}{c}\stackrel{.}{\tilde{q}}{\text{'}}_{\mathrm{e}}=\frac{1}{2}\left(\mathrm{q}{\text{'}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left(\tilde{q}{\text{'}}_{\mathrm{e}}\right)\right)\left(\stackrel{\rightharpoonup }{v}-\stackrel{\rightharpoonup }{z}\right)\\ \stackrel{.}{q}{\text{'}}_{\mathrm{e}0}=-\frac{1}{2}{\left(\tilde{q}{\text{'}}_{\mathrm{e}}^{}\right)}^{\mathrm{T}}\left(\stackrel{\rightharpoonup }{v}-\stackrel{\rightharpoonup }{z}\right)\end{array}$

(33)

Moreover, as shown above, it is asymptotically stable for $\stackrel{\rightharpoonup }{z}=0$ . The differential equation for $\stackrel{\rightharpoonup }{z}$ is:

$\begin{array}{l}\stackrel{.}{\stackrel{\rightharpoonup }{z}}=-\stackrel{.}{\stackrel{\rightharpoonup }{w}}-{\mathrm{k}}_1\stackrel{.}{\tilde{q}}{\text{'}}_{\mathrm{e}}\\ ={\mathrm{J}}^{-1}\left(\mathrm{S}\left(\stackrel{\rightharpoonup }{w}\right)J\stackrel{\rightharpoonup }{w}-\mathrm{u}\right)+\frac{1}{2}{\mathrm{k}}_1^2\left(\mathrm{q}{\text{'}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left(\tilde{q}{\text{'}}_{\mathrm{e}}\right)\right)\tilde{q}{\text{'}}_{\mathrm{e}}\\ =\mathrm{f}\left(\tilde{q}{\text{'}}_{\mathrm{e}},\stackrel{\rightharpoonup }{z}\right)-{\mathrm{J}}^{-1}\stackrel{\rightharpoonup }{u}\end{array}$

(34)

where:

$\begin{array}{l}\mathrm{f}\left(\tilde{q}{\text{'}}_{\mathrm{e}},\stackrel{\rightharpoonup }{z}\right)=-{\mathrm{k}}_1{\mathrm{J}}^{-1}\mathrm{S}\left(-{\mathrm{k}}_1\stackrel{.}{\tilde{q}}{\text{'}}_{\mathrm{e}}-\stackrel{\rightharpoonup }{z}\right)J\tilde{q}{\text{'}}_{\mathrm{e}}-\\ \mathrm{\dots }{\mathrm{J}}^{-1}\mathrm{S}\left(-{\mathrm{k}}_1\stackrel{.}{\tilde{q}}{\text{'}}_{\mathrm{e}}-\stackrel{\rightharpoonup }{z}\right)J\stackrel{\rightharpoonup }{z}-\frac{1}{2}{\mathrm{k}}_1^2\left(\mathrm{q}{\text{'}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left(\tilde{q}{\text{'}}_{\mathrm{e}}\right)\right)\tilde{q}{\text{'}}_{\mathrm{e}}\end{array}$

(35)

We want to find $\mathrm{u}=\mathrm{u}\left(\tilde{q}{\prime }_{\mathrm{e}},\stackrel{\rightharpoonup }{z}\right)$ , such that the system of (33) and (34) is globally asymptotically stable. Let $\mathrm{x}={\left[\begin{array}{cc}\tilde{q}{\prime }_{\mathrm{e}}& \stackrel{\rightharpoonup }{z}\end{array}\right]}^{\mathrm{T}}$ , the system model becomes:

$\stackrel{.}{x}=\left[\begin{array}{c}-\frac{1}{2}{\mathrm{k}}_1\left(\mathrm{q}{\text{'}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left(\tilde{q}{\text{'}}_{\mathrm{e}}\right)\right)\tilde{q}{\text{'}}_{\mathrm{e}}\\ \mathrm{f}\left(\tilde{q}{\text{'}}_{\mathrm{e}},\stackrel{\rightharpoonup }{z}\right)\end{array}\right]-\left(\begin{array}{cc}0& 0\\ 0& {\mathrm{J}}^{-1}\end{array}\right)\left(\begin{array}{c}0\\ \mathrm{u}\end{array}\right)$

(36)

Let:

$\mathrm{F}\left(\mathrm{x}\right)=\left[\begin{array}{c}-\frac{1}{2}{\mathrm{k}}_1\left(\mathrm{q}{\text{'}}_{\mathrm{e}0}\mathrm{I}+\mathrm{S}\left(\tilde{q}{\text{'}}_{\mathrm{e}}\right)\right)\tilde{q}{\text{'}}_{\mathrm{e}}\\ \mathrm{f}\left(\tilde{q}{\text{'}}_{\mathrm{e}},\stackrel{\rightharpoonup }{z}\right)\end{array}\right]$

(37)

and:

$\mathrm{G}\left(\mathrm{x}\right)=\left(\begin{array}{cc}0& 0\\ 0& -{\mathrm{J}}^{-1}\end{array}\right)$

(38)

The system model can be rewritten as:

$\stackrel{.}{x}=\mathrm{F}\left(\mathrm{x}\right)+\mathrm{G}\left(\mathrm{x}\right)\left(\begin{array}{c}0\\ \stackrel{\rightharpoonup }{u}\end{array}\right)$

(39)

Since the derivative of the control Lyapunov function is:

$\stackrel{.}{V}\left(\mathrm{x}\right)={\mathrm{L}}_{\mathrm{F}}\mathrm{V}\left(\mathrm{x}\right)+{\mathrm{L}}_{\mathrm{G}}\mathrm{V}\left(\mathrm{x}\right)\stackrel{\rightharpoonup }{u}$

(40)

Provided that applying the inverse optimal control theorem:

$\stackrel{\rightharpoonup }{u}=-\beta \cdot {\mathrm{N}}^{-1}\left(\mathrm{x}\right){\left({\mathrm{L}}_{\mathrm{G}}\mathrm{V}\right)}^{\mathrm{T}}$

(41)

could make V̇(x) negative with a sufficiently high gain, then the maximum rate control law is given by:

$\stackrel{\rightharpoonup }{u}=-{\mathrm{u}}_{\max }\cdot SGN\left({\left({\mathrm{L}}_{\mathrm{G}}\mathrm{V}\right)}^{\mathrm{T}}\right)$

(42)

where u_max represents the maximum value of input saturation.

Notice that the SGN function is an extended sign function for a multidimensional situation. Since the sign function can be seen as the direction function of one dimension, i.e.:

$\mathrm{sgn}\left(\lambda \right)=\{\begin{array}{c}\lambda /\left|\lambda \right|\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\lambda \ne 0\\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\lambda =0\end{array}$

(43)

we extend it to a three dimensional situation, such as:

$SGN\left(\mathrm{X}\right)=\{\begin{array}{l}\mathrm{X}/\Vert \mathrm{X}\Vert \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{X}\ne [0\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}0{]}^{\mathrm{T}}\\ [0\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}0{]}^{\mathrm{T}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{X}=[0\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}0{]}^{\mathrm{T}}\end{array}$

(44)

Consider the new Lyapunov function:

${\mathrm{V}}_2={\mathrm{x}}^{\mathrm{T}}\mathrm{x}$

(45)

Then, the maximum rate control law (with respect to the designed Lyapunov function) is:

$\begin{array}{c}\stackrel{\rightharpoonup }{u}=-{\mathrm{u}}_{\max }\cdot SGN\left(-\mathrm{J}\left(\tilde{q}{\text{'}}_{\mathrm{e}}+\mathrm{z}\right)\right)\\ =-{\mathrm{u}}_{\max }\cdot SGN\left(\mathrm{J}\left(\left({\mathrm{k}}_1-1\right)\tilde{q}{\text{'}}_{\mathrm{e}}+\mathrm{w}\right)\right)\end{array}$

(46)

Notice that the backstepping-based inverse optimal controller (BIOAC) that we designed is a sliding mode controller with the sliding surface:

$\mathrm{s}=\mathrm{J}\left(\left({\mathrm{k}}_1-1\right)\tilde{q}{\text{'}}_{\mathrm{e}}+\stackrel{\rightharpoonup }{w}\right)$

(47)

It has the inherent property of robustness as to matched model uncertainty and explicit disturbances. However, since the sign function SGN(·)will cause a high frequency oscillation in practice, the sliding mode control law will suffer the chattering problem. For removing chattering in sliding mode control, there are a few methods. A more common approach is to replace the sign function by a smooth function with a boundary layer of the form:

$\eta \left(\mathrm{s}\right)=\frac{\mathrm{s}}{\epsilon +\Vert \mathrm{s}\Vert }$

(48)

where ε is a constant vector to represent the boundary layer. However, with a boundary layer the BIOAC sacrifices the robustness in relation model uncertainty and explicit disturbances.

4. Numerical Simulation

In this section, the performance of the designed BIOAC is compared with a PD controller using numerical simulation in the MATLAB environment. The system parameters used in the simulation are:

$\begin{array}{l}{\mathrm{I}}_1=0.0079kg\cdot {\mathrm{m}}^2,{\mathrm{I}}_2=0.0079kg\cdot {\mathrm{m}}^2,{\mathrm{I}}_3=0.0093kg\cdot {\mathrm{m}}^2\\ \mathrm{d}=0.2\mathrm{m},{\mathrm{c}}_{\mathrm{f}}=0.00003,{\mathrm{c}}_{\tau \mathrm{f}}=0.01\end{array}$

The initial conditions are:

$\begin{array}{l}\stackrel{\rightharpoonup }{q}\left(0\right)=(0.159,0.57,0.57,0.57{)}^{\mathrm{T}}\\ \stackrel{\rightharpoonup }{w}\left(0\right)={\left(0,0,0\right)}^{\mathrm{T}}\end{array}$

And the desired condition are:

$\begin{array}{l}{\stackrel{\rightharpoonup }{q}}_{\mathrm{d}}={\left(1,0,0,0\right)}^{\mathrm{T}}\\ {\stackrel{\rightharpoonup }{w}}_{\mathrm{d}}={\left(0,0,0\right)}^{\mathrm{T}}\end{array}$

This corresponds to the initial attitude errors of 161.7 degrees around axis (1,1,1)^T.

For a more obvious illustration of control performance, the PD controller which provides a near-eigenaxis rotation [18] is used as a comparison. The PD algorithm is described as:

${\stackrel{\rightharpoonup }{u}}_{pd}=\mathrm{S}\left(\stackrel{\rightharpoonup }{w}\right)J\stackrel{\rightharpoonup }{w}-kJ{\tilde{q}}_{\mathrm{e}}-dJ\stackrel{\rightharpoonup }{w}$

(49)

In the simulation process of the PD controller, the control effort limitation is never settled, and in a maximum rate control simulation the maximum control effort is (0.05N · m, 0.05N · m, 0.05N · m).

Appearing in the attitude quaternion curves in Fig. 4 and the output control torque curve in Fig. 5, the BIOAC is faster than the PD with a control effort under the maximum limitations. Moreover, the PD controller seeks a greater control torque when the error of orientation is large, and will suffer the problem of control failure in practice.

Figure 4.

Time histories of the quaternion

Figure 5.

Time histories of the control torques

5. Experiments

Given a good performance in the MATLAB simulation, BIOAC is realized on a test bench to show its effectiveness in practice. The physical device being used in the work discussed here is shown in Fig. 6. The experimental system consists of a posture platform, an aircraft body, an onboard controller, an attitude and heading reference system (AHRS), a wireless serial communication interface, a remote controller and receiver, and a computer.

Figure 6.

Physical devices for the experimental system

The aircraft body selected is the XAircraft X-650 quadrotor, which is made of carbon fibre, equipped with four brushless DC motors and a motor speed control device. It has an empty weight of about 800g and the maximum single propeller thrust is 4N. The Attitude and Heading Reference System is the Xsens company MTi inertial measurement unit (IMU), which can provide 3 axes of driftless orientation and 3 axes of calibrated angular velocity. Its maximum update rate is 100Hz, with a 0.5 degree attitude measurement accuracy, a 1 degree heading measure accuracy, and a 300 degree per second maximum measurable angular velocity. The remote controller and receiver is a WFLY FT06-A with a 5-way continuous pulse channel and a switching channel. TI's TMS320F2809 chip with a 100MHz core frequency is the core of the onboard controller.

The connection diagram of experimental system signals is shown in Fig 7. In every control cycle (10ms), the F2809 chip is interrupted to receive orientation and angular velocity information from the IMU and then calculate the control pulse width and transfer it to the motor speed control devices. Outside the control loop, the F2809 is interrupted to turn the parameters and sends a set of state data to the computer in every 100ms. The work of receiving the remote control instruction is done in the same interrupter program.

Figure 7.

Signal connection diagram

In the experiment, the attitude of the quadrotor switch is between:

$\begin{array}{l}{\stackrel{\rightharpoonup }{q}}_1=(0.9848,0.1003,0.1003,0.1003{)}^{\mathrm{T}}\\ {\stackrel{\rightharpoonup }{w}}_1={\left(0,0,0\right)}^{\mathrm{T}}\end{array}$

and:

$\begin{array}{l}{\stackrel{\rightharpoonup }{q}}_2={\left(1,0,0,0\right)}^{\mathrm{T}}\\ {\stackrel{\rightharpoonup }{w}}_2={\left(0,0,0\right)}^{\mathrm{T}}\end{array}$

i.e., the quadrotor rotates 20 degrees around the axis of (1,1,1)^T

In the experiment, J and c_τf is unknown. However, since BIOAC is a sliding mode controller, it can work effectively with the settings of J = I and c_τf = 1.

The attitude curves in Fig. 8 show that the controlled system is stable and, after the switch instruction, the attitude of the quadrotor converges to the desired attitude fast. There are small static errors in Fig. 8 and Fig. 9. The reasons for this are the incomplete compensation and the existence of the boundary layer. The results would be improved with a more accurate estimation of J and c_τf.

Figure 8.

Time histories of the quaternion

Figure 9.

Time histories of the input instructions

6. Conclusions

The maximum rate attitude control problem of the quadrotor is studied under the constraint of input saturation. Due to the difficulty in obtaining closed-form solutions to the HJB equation, the inverse optimal approach is used. Following the backstepping phases, a global stable controller is derived. With respect to a classic PD controller, the numerical simulation results show that BIOAC could achieve faster convergence with control effort limitation. The realization of the BIOAC is gained in an experimental system and shows the practical effectiveness of being free from computational issues, which is highly desirable given the limited resources onboard.

7. Acknowledgements

The authors would like to thank Dr. HE Yunze in NUDT for providing advice and Dr. SHEN Tao and SONG Baoquan in NUDT for their help in constructing the test bench.