Incremental Backstepping Sliding-Mode Trajectory Control for Tailless Aircraft with Stability Enhancer

Abstract

This paper presents an incremental backstepping sliding-mode (IBS) controller for trajectory control of a tailless aircraft with unknown disturbances and model uncertainties. The proposed controller is based on a nonlinear dynamic model of the tailless aircraft. A stability enhancer (SE) that limits both the rate and amplitude of the virtual control input is proposed. The stability enhancer consists of two layers. When the virtual control input approaches the edge, the first layer SE would be activated to modify the trajectory tracking error; when the virtual control input exceeds the edge, the second layer SE would reduce the control gains to make sure the virtual control input drops within the edge as soon as possible. With the help of SE, the incremental control method could be extended to outer-loop control without considering the dynamics of the inner-loop system. In addition, an adaptive estimator for state derivatives is proposed, together with IBS, allowing the controller to show excellent robustness. Finally, two simulations are presented. The first simulation shows that the system is insensitive to external disturbances and model uncertainties, and the effectiveness of SE is proved in the second simulation.

1. Introduction

In past decades, tailless aircraft have attracted widespread attention. Due to excellent stealth performance, tailless aircraft have already been successfully used in the military. The future battlefield environment needs the next-generation combat aircraft to be low observable as well as super maneuverable. This goal prompts next-generation combat aircraft to adopt a tailless design with a small aspect ratio. However, this kind of aerodynamic layout brings stronger coupling, nonlinear, and non-affine effects [1,2,3], challenging its controller design.

The most widely used framework in flight control is the gain-scheduled linear feedback controller with Jacobian linearization around specific operational points. Traditionally, this framework is an effective way to deal with nonlinearity and is still used now because well-established linear control theories can support its design. Such approaches, however, have failed to address the significant nonlinearities throughout the flight envelope, and, because of the linearized model it uses, the stability between the operational points is always questionable.

Nonlinear control methods do not suffer from the above shortcomings. The essence of these nonlinear control methods is to cancel the nonlinearity through feedback, and then the linear control technique could be applied. Backstepping (BS) is a typical nonlinear control method. To exactly cancel the nonlinear dynamic, BS highly depends on accurate model knowledge. However, for tailless aircraft, it is not easy to extract useful model knowledge, such as control effectiveness, from aerodynamic data. Specifically, due to the interaction between the state variables and control inputs, the structure of the aircraft may not be very apparent, and the control input may appear non-affine, making the BS very hard to deal with.

Of course there have been many studies concerning non-affine systems using BS, but it should be noted that a strict decoupling condition $0⩽f_l⩽\frac{\partial f(x,u)}{\partial u}⩽f_u$ is used in many of them [4,5,6,7,8,9]. However, it is hard to prove that such a condition is satisfied in practice, and the application of the condition will cause the many accessible details of actuator dynamics to be underused, making it a very difficult work to choose the controller’s parameters. In addition, this condition is very difficult to extend to multiple control input systems, especially if there are coupling effects between control inputs.

To reduce the dependence on the model information, incremental backstepping (IB) was proposed by [10,11,12,13] and others. As an effective control method, IB has received considerable scholarly attention in the flight control field ever since it was proposed, and has achieved quite fruitful results [14,15,16,17,18]. Instead of relying on model information, IB cancels the nonlinearity of the system through the feedback of state derivatives. Because state derivatives contain all the information of the model and disturbances, IB no longer needs accurate model information. The only model information that IB requires is control effectiveness. Recent studies [19,20,21] have found a condition under which the stability of the IB-controlled system will not be degraded by the uncertainties in control effectiveness; by increasing the sample rate, the effect of uncertainty in control effectiveness can be largely attenuated. Generally, IB obtains control effectiveness through numerical differentiation of aerodynamic data. With this condition, the designer of the IB controller would know how accurate the control effectiveness needs to be to keep the system stable, so the IB will not waste time in obtaining over-accurate control effectiveness, making IB response more faster and more applicable.

In addition, IB is also a powerful tool to cope with a non-affine system, as it transforms a non-affine system to an incremental affine form by Taylor expansion. Therefore, there is no need to consider complex model structures in IB control law design. In recent years, tailless aircraft tend to take a multiple control surfaces design to improve maneuverability. For example, the Innovative Control Effector (ICE) aircraft is equipped with 11 control surfaces with overlapping functionality [2]. Because of the compact layout, there are strong coupling effects between control surfaces and state variables, making the aerodynamic model very complex, and the control inputs appear non-affine. Therefore, for this kind of non-affine system with multiple control inputs, the incremental control methods seems to be the only option.

However, the derivation of IB is based on the time scale separation principle, which means that the control action should change far faster than the state [22]. Thus, the actuator dynamics are very important to IB. When the command cannot be implemented fast enough, the stability of the system will degrade [23]. Because the dynamic of the inner-loop state is usually non-negligible, it would be hard to tell whether the command can be implemented fast enough. Therefore, most studies only use IB in the last step of controller design, and the application of an incremental control method in the outer loop of the cascade system is hardly seen. Because of the time delay caused by inner-loop dynamics, the time scale separation principle is longer applicable. Specifically, due to the strong nonlinearity of a tailless vehicle, incremental control methods should be an effective tool that can be used in trajectory tracking, but the time scale separation principle limits most recent studies from using an incremental control method in attitude control. Therefore, it would be a very meaningful work to find a effective way that could alleviate the time scale separation requirement and extend the incremental control to outer-loop control.

From another perspective, however, if the rate and amplitude of the virtual control input produced by the outer-loop controller are limited, it would be more reasonable to assume the instantaneous reaction of the inner-loop dynamics, which would guarantee the stability of the system. This brings to mind that many anti-windup measures [24,25] could meet such needs. Thus, a two-layer stability enhancer (SE) is proposed in this paper. The first layer SE works the same as a traditional anti-windup measure, which modifies the feedback of the outer-loop system when the virtual control input is going to run out of edge. When the virtual control input runs out of edge, the second layer SE will reduce the control gain to make sure the virtual control input falls within the edge as soon as possible.

With the rapid development of hardware in recent years, the controller’s sample rate is getting higher. This allows the incremental control method to be applied more often, and some satisfactory results are seen in [26,27,28,29]. The advantages of incremental control noted above make it a promising application in the field of flight control. Up to now, there are only a few kinds of incremental control methods, such as incremental nonlinear dynamic inversion and incremental backstepping. It would be very meaningful work to develop an incremental version for other classic control methods.

Motivated by the aforementioned discussion, this paper proposes an incremental backstepping sliding-mode (IBS) trajectory controller with stability enhancer (SE) for tailless aircraft featuring strong nonlinearity and non-affine input. The SE considers both rate and amplitude limits. When the virtual control input approaches the virtual edge, the first layer SE will be activated to modify the tracking error; when the virtual control input exceeds the virtual edge, the second layer SE will be activated, and the control gain will be reduced to make sure the virtual control input is within the limits. With the help of the SE, the proposed IBS is extended to outer-loop control. Considering the condition proposed in [21], only a linear spline is used to obtain the control effectiveness, which reduces the computational load and maintains the performance. Together with an adaptive estimator for state derivation, the proposed control method shows robustness against external disturbances and model uncertainties.

The paper is arranged as follows: Section 2 presents the aircraft model. The description of the problem is given in Section 3, which is followed by the main result in Section 4, where the control law and stability analysis is presented. Section 5 gives the simulation results. With the proposed stability enhancer, the incremental control methods could be extended to outer-loop control for the cascade system. The proposed control method also show robustness against model uncertainties and disturbance.

2. Nonlinear Tailless Aircraft Model

In this section, the tailless aircraft model used for controller design is introduced. This model originated from Lockheed Martin’s Innovative Control Effector (ICE) project [30]. The unique shape of the ICE aircraft makes it a potential option for future air combat aircraft. In [2,3], the aerodynamic characteristics of the ICE aircraft were analyzed, and it was found that the strong nonlinearity of the ICE aircraft yields great challenge for its flight control. Based on the aerodynamic data given in [1], we could build a high-fidelity data-based model. The 6-DOF aircraft motion equation [31] is given as follows.

Define the position vector in the earth fixed coordinate ${\mathit{x}}_{\mathbf{0}}={\left[XYZ\right]}^T$, and the derivatives of ${\mathit{x}}_{\mathbf{0}}$ are:

$${\dot{\mathit{x}}}_0=\left[\begin{array}{c}\dot{X}\\ \dot{Y}\\ \dot{Z}\end{array}\right]=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right]=\left[\begin{array}{c}V{c}_{\gamma }{c}_{\chi }\\ V{c}_{\gamma }{s}_{\chi }\\ -V{s}_{\gamma }\end{array}\right]$$

(1)

where $v_x$, $v_y$, $v_z$ are velocity components in the x, y, z directions; $V=\sqrt{v_x^2+v_y^2+v_z^2}$ is the total velocity, $\chi =arctan\left(\frac{v_y}{v_x}\right)$ is the ground tracking angle, $\gamma =arcsin\left(\frac{v_z}{V}\right)$ is the flight path angle, $c_{★}$ represent $cos★$, and $s_{★}$ represent $sin★$. Define ${\mathit{x}}_{\mathbf{1}}=\left[v_x{v}_y{v}_z\right]$, and the dynamics of ${\mathit{x}}_{\mathbf{1}}$ are [32]:

$${\dot{\mathit{x}}}_1=\left[\begin{array}{cccc}& c_{\chi }{c}_{\gamma }& -V{s}_{\chi }{c}_{\gamma }& -V{c}_{\gamma }{s}_{\gamma }\\ & s_{\chi }{c}_{\gamma }& V{c}_{\chi }{c}_{\gamma }& -V{s}_{\gamma }{s}_{\gamma }\\ & -s_{\gamma }& 0& -V{c}_{\gamma }\end{array}\right]·\left[\begin{array}{c}\dot{V}\\ \dot{\chi }\\ \dot{\gamma }\end{array}\right]$$

(2)

and the derivative of ${\left[V\chi \gamma \right]}^T$ is:

$$\left[\begin{array}{c}\dot{V}\\ \dot{\chi }\\ \dot{\gamma }\end{array}\right]=\left[\begin{array}{cccc}& \frac{1}{M}& 0& 0\\ & 0& \frac{1}{MV{c}_{\gamma }}& 0\\ & 0& 0& \frac{-1}{MV}\end{array}\right]·\left[{\mathit{T}}_{\mathit{ve}}{\mathit{G}}_{\mathit{f}}+{\mathit{T}}_{\mathit{vb}}{\mathit{T}}_{\mathit{f}}+\left[\begin{array}{cccc}& 1& 0& 0\\ & 0& s_{\mu }& c_{\mu }\\ & 0& c_{\mu }& -s_{\mu }\end{array}\right]·{\mathit{F}}_{\mathit{a}}\right]$$

(3)

where ${\mathit{G}}_{\mathit{f}}={\left[00Mg\right]}^T$ represents gravitational forces, ${\mathit{T}}_{\mathit{f}}={\left[T00\right]}^T$ represents engine thrust, and ${\mathit{F}}_{\mathit{a}}={\left[DLY\right]}^T$ represents aerodynamic force, which is defined by:

$${\mathit{F}}_{\mathit{a}}=\left[\begin{array}{c}D\\ L\\ Y\end{array}\right]=\overline{q}S\left[\begin{array}{c}{C}_D\\ C_L\\ C_Y\end{array}\right]=\overline{q}S\left[\begin{array}{c}{C}_{D,base}(\alpha ,\beta ,V)+\sum _{i=1}^j{C}_{D,i}(\alpha ,\beta ,\delta )\\ C_{L,base}(\alpha ,\beta ,V)+\sum _{i=1}^j{C}_{L,i}(\alpha ,\beta ,\delta )\\ C_{Y,base}(\alpha ,\beta ,V)+\sum _{i=1}^j{C}_{Y,i}(\alpha ,\beta ,\delta )\end{array}\right]$$

(4)

Because the ICE aircraft has 11 independent control surfaces, $\mathbf{\delta }={[{\delta }_1,{\delta }_2,\cdots ,{\delta }_{11}]}^T$ represents the control surface deflection vector, and $C_{·,base}$ and ${\sum }_{i=1}^j{C}_{·,i}(\alpha ,\beta ,\delta )$ are the aerodynamic coefficients generated by the body and control surfaces.

Considering the aircraft capacity and flight safety, the amplitude of $T_f,\alpha ,\mu$ as well as their rates should be limited. Define ${\mathit{x}}_{\mathbf{2}}^{\prime }={\left[T_f\alpha \mu \right]}^T\doteq {\left[x_{2,1}^{\prime }{x}_{2,2}^{\prime }{x}_{2,3}^{\prime }\right]}^T$ and denote its upper and lower limits of amplitude as $\overline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }}$ and $\underline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }}$ and the upper and lower limits of rate as $\overline{{\dot{\mathit{x}}}_{2s}^{\prime }}$ and $\underline{{\dot{\mathit{x}}}_{2s}^{\prime }}$.

The expression ${\mathit{T}}_{\mathit{vb}}$ is the transformation matrix from the body frame to the velocity frame, and ${\mathit{T}}_{\mathit{ve}}$ is the transformation matrix from the earth frame to the velocity frame. These matrices are given as [33]:

$${\mathit{T}}_{vb}=\left[\begin{array}{cccc}& c_{\alpha }{c}_{\beta }& s_{\beta }& s_{\alpha }{c}_{\beta }\\ & -c_{\alpha }{s}_{\beta }{c}_{\mu }+s_{\alpha }{s}_{\mu }& c_{\beta }{c}_{\mu }& -s_{\alpha }{s}_{\beta }{c}_{\mu }-c_{\alpha }{s}_{\mu }\\ & -c_{\alpha }{s}_{\beta }{s}_{\mu }-s_{\alpha }{c}_{\mu }& c_{\beta }{s}_{\mu }& -s_{\alpha }{s}_{\beta }{s}_{\mu }+c_{\alpha }{c}_{\mu }\end{array}\right]$$

(5)

$${\mathit{T}}_{ve}=\left[\begin{array}{cccc}& c_{\chi }{c}_{\gamma }& s_{\chi }{c}_{\gamma }& s_{\gamma }\\ & -s_{\chi }& c_{\chi }& 0\\ & c_{\chi }{s}_{\gamma }& s_{\chi }{s}_{\gamma }& c_{\gamma }\end{array}\right]$$

(6)

where $\mu$ is the bank angle, $\alpha$ is the angle of attack, and $\beta$ is the sideslip angle.

Define ${\mathit{x}}_2={\left[\mu \alpha \beta \right]}^T$. The derivatives of ${\mathit{x}}_2$ are given as:

$${\dot{\mathit{x}}}_2={\left[\begin{array}{cccc}& c_{\alpha }{c}_{\beta }& 0& s_{\alpha }\\ & s_{\beta }& 1& 0\\ & s_{\alpha }{c}_{\beta }& 0& -c_{\alpha }\end{array}\right]}^{-1}\left({\mathit{T}}_{bv}\left[\begin{array}{c}-\dot{\chi }{s}_{\gamma }\\ \dot{\gamma }\\ \dot{\chi }{c}_{\gamma }\end{array}\right]+\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\right)$$

(7)

where p, q, and r are the body-axis roll, pitch, and yaw rates.

Further define the angular rate vector as ${\mathit{x}}_3={\left[pqr\right]}^T$. Its kinematics are described as:

$${\dot{\mathit{x}}}_3={\mathit{J}}^{-1}({\mathit{M}}_{\mathit{a}}-{\mathit{x}}_3\times \mathit{J}{\mathit{x}}_3)$$

(8)

where $\mathit{J}$ is moment of inertia matrix, defined as:

$$\mathit{J}=\left[\begin{array}{cccc}& I_{xx}& -I_{xy}& -I_{xz}\\ & -I_{xy}& I_{yy}& -I_{yz}\\ & -I_{xz}& -I_{yz}& I_{zz}\end{array}\right]$$

(9)

and ${\mathit{M}}_{\mathit{a}}=\left[lmn\right]$ is the aerodynamic moment, defined by:

$${\mathit{M}}_{\mathit{a}}=\left[\begin{array}{c}l\\ m\\ n\end{array}\right]=\overline{q}S\left[\begin{array}{c}b·C_l\\ \overline{c}·C_m\\ b·C_n\end{array}\right]=\overline{q}S\left[\begin{array}{c}b·(C_{l,base}(\alpha ,\beta ,V)+\sum _{i=1}^j{C}_{l,i}(\alpha ,\beta ,\delta ))\\ \overline{c}·(C_{m,base}(\alpha ,\beta ,V)+\sum _{i=1}^j{C}_{m,i}(\alpha ,\beta ,\delta ))\\ b·(C_{n,base}(\alpha ,\beta ,V)+\sum _{i=1}^j{C}_{n,i}(\alpha ,\beta ,\delta ))\end{array}\right]$$

(10)

Remark 1.

The aircraft model noted before is constructed using the spline technique based on the discrete aerodynamic data provided in [1]. Considering the controller’s efficiency, only linear interpolation is used in the controller design. Of course, the aerodynamic coefficients obtained through linear interpolation cannot be very accurate. In the following sections, we will prove that these errors will not affect the performance of the controller we propose.

3. Description of Problem

Based on the model in Section 2, the trajectory tracking problem will be discussed in the following context.

Define ${\mathit{x}}_{\mathbf{2}}^{\prime }={\left[T\mu \alpha \right]}^T\doteq {\left[x_{2,1}^{\prime }{x}_{2,2}^{\prime }{x}_{2,3}^{\prime }\right]}^T$. Accroding to Section 2, we have ${\mathit{x}}_0={\left[XYZ\right]}^T$, ${\mathit{x}}_1={\left[v_x{v}_y{v}_z\right]}^T$, ${\mathit{x}}_2={\left[\mu \alpha \beta \right]}^T$, ${\mathit{x}}_3={\left[pqr\right]}^T$, and ${\mathit{x}}_1$ and ${\mathit{x}}_3$ appear affine in ${\dot{\mathit{x}}}_0$ and ${\dot{\mathit{x}}}_2$, ${\mathit{x}}_2^{\prime }$ and $\delta$ appear non-affine in ${\dot{\mathit{x}}}_1$ and ${\dot{\mathit{x}}}_3$. We rewrite the system in the following compact form [34]:

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_0& ={\mathit{x}}_1\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_1& ={\mathit{f}}_1\left({\mathit{x}}_1\right)+{\mathit{g}}_1({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_2& ={\mathit{f}}_2({\mathit{x}}_1,{\mathit{x}}_2,\delta )+{\mathit{g}}_2·{\mathit{x}}_3\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_3& ={\mathit{f}}_3\left({\mathit{x}}_3\right)+{\mathit{g}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\delta )\hfill \end{array}$$

(14)

Remark 2.

As to ${\mathit{f}}_1$, ${\mathit{f}}_2$, ${\mathit{g}}_2$, and ${\mathit{f}}_3$, their corresponding definitions are obvious in Section 2. Because ${\mathit{x}}_2^{\prime }$ and δ appear nonlinearly in $C_{·,base}$ and $C_{·,i}$, ${\mathit{g}}_1$ and ${\mathit{g}}_3$ are seen as non-affine functions here.

Assumption 1.

The aircraft system is sufficiently time scale separated. This means that, in a small time increment, changes of control input cause the state derivative that is directly affected to change much faster than the controlled one.

According to Assumption 1, it is obvious that we see ${\mathit{x}}_1$, ${\mathit{x}}_3$, and $\delta$ as the virtual control input of Equations (11), (13), and (14). Note that from Equations (2)–(4), we know ${\mathit{g}}_1$ should be the function of ${\mathit{x}}_1$, ${\mathit{x}}_2$, $\delta$, and T, we choose ${\mathit{x}}_2^{\prime }={\left[T_f\alpha \mu \right]}^T$ as the virtual control input of Equation (12). Because the aerodynamic forces contributed by $\delta$ are relatively small as they are designed to produce moments, and for the flight safety, the command signal of $\beta$ is usually set as ${\beta }^c=0$. However, due to the unique control surface setting of the ICE aircraft, the aerodynamic force $\delta$ generated by the control surface cannot be neglected in ${\mathit{g}}_1$.

During flight, disturbances are ubiquitous, which will affect the stability and performance of the flight control system. Neither disturbances from modeling error nor external disturbances are easy to measure in practice, making the disturbance rejection a challenge of controller design. In this paper, the modeling error comes from the acquisition of aerodynamic force and moment, which are involved in ${\mathit{g}}_1$, ${\mathit{f}}_2$, and ${\mathit{g}}_3$; thus we define:

$$\begin{array}{cc}\hfill & {\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })={\mathit{g}}_1({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })-{\mathit{g}}_{1,e}({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & {\overline{\mathit{f}}}_2({\mathit{x}}_1,{\mathit{x}}_2,\delta )={\mathit{f}}_2({\mathit{x}}_1,{\mathit{x}}_2,\delta )-{\mathit{f}}_{2,e}({\mathit{x}}_1,{\mathit{x}}_2,\delta )\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & {\overline{\mathit{g}}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\delta )={\mathit{g}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\delta )-{\mathit{g}}_{3,e}({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\delta )\hfill \end{array}$$

(17)

where ${\mathit{g}}_{1,e}$, ${\mathit{f}}_{2,e}$, and ${\mathit{g}}_{3,e}$ represent the model uncertainties, ${\overline{\mathit{g}}}_1$, ${\overline{\mathit{f}}}_2$ and ${\overline{\mathit{g}}}_3$ are konwn dynamics used for control law design.

The external disturbance is assumed to be time varying and exists in the form of aerodynamic force and moment, which are defined as ${\mathit{d}}_F={[d_{F1},d_{F2},d_{F3}]}^T$ and ${\mathit{d}}_M={[d_{M1},d_{M2},d_{M3}]}^T$, respectively. Equations (12) and (14) can be modified as:

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_1& ={\mathit{f}}_1\left({\mathit{x}}_1\right)+{\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })+{\mathit{d}}_1\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_2& ={\overline{\mathit{f}}}_2({\mathit{x}}_1,{\mathit{x}}_2,\delta )+{\mathit{g}}_2·{\mathit{x}}_3+{\mathit{d}}_2\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_3& ={\mathit{f}}_3\left({\mathit{x}}_3\right)+{\overline{\mathit{g}}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\delta )+{\mathit{d}}_3\hfill \end{array}$$

(20)

where ${\mathit{d}}_1=M^{-1}{\mathit{d}}_F+{\mathit{g}}_{1,e}({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })$, ${\mathit{d}}_2={\mathit{f}}_{2,e}$, and ${\mathit{d}}_3={\mathit{J}}^{-1}{\mathit{d}}_M+{\mathit{g}}_{3,e}({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })$.

Considering the non-affine input in ${\dot{\mathit{x}}}_1$ and ${\dot{\mathit{x}}}_3$, the incremental control law is only used in velocity and angular rate control loops. Thus, using a one-order Taylor expansion, we rewrite $\dot{{\mathit{x}}_1}$ and $\dot{{\mathit{x}}_3}$ in incremental form:

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_1& ={\dot{\mathit{x}}}_{1,0}+{\mathit{A}}_1\Delta {\mathit{x}}_1+{\mathit{B}}_1\Delta \beta +{\mathit{C}}_1\Delta \delta +(\overline{{\mathit{G}}_1}+{\mathit{G}}_{1,e})\Delta {\mathit{x}}_2^{\prime }+\Delta {\mathit{d}}_1\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_3& ={\dot{\mathit{x}}}_{3,0}+{\mathit{A}}_3\Delta {\mathit{x}}_1+{\mathit{B}}_3\Delta {\mathit{x}}_2+{\mathit{C}}_3\Delta {\mathit{x}}_3+(\overline{{\mathit{G}}_3}+{\mathit{G}}_{3,e})\Delta \delta +\Delta {\mathit{d}}_3\hfill \end{array}$$

(22)

$$\begin{array}{c}\hfill {\mathit{A}}_1=\frac{\partial {\mathit{f}}_1}{\partial {\mathit{x}}_1}+\frac{\partial {\mathit{g}}_1}{\partial {\mathit{x}}_1}{\mathit{B}}_1=\frac{\partial {\mathit{g}}_1}{\partial \beta }{\mathit{C}}_1=\frac{\partial {\mathit{g}}_1}{\partial \delta }{\overline{\mathit{G}}}_1=\frac{\partial {\overline{\mathit{g}}}_1}{\partial {\mathit{x}}_2^{\prime }}{\mathit{G}}_{1,e}=\frac{\partial {\mathit{g}}_{1,e}}{\partial {\mathit{x}}_2^{\prime }}\end{array}$$

$$\begin{array}{c}\hfill {\mathit{A}}_3=\frac{\partial {\mathit{f}}_3}{\partial {\mathit{x}}_3}+\frac{\partial {\mathit{g}}_3}{\partial {\mathit{x}}_3}{\mathit{B}}_3=\frac{\partial {\mathit{g}}_3}{\partial {\mathit{x}}_1}{\mathit{C}}_3=\frac{\partial {\mathit{g}}_3}{\partial {\mathit{x}}_2}{\overline{\mathit{G}}}_3=\frac{\partial {\overline{\mathit{g}}}_3}{\partial \delta }{\mathit{G}}_{3,e}=\frac{\partial {\mathit{g}}_{3,e}}{\partial \delta }\end{array}$$

where ${\mathit{A}}_1,{\mathit{A}}_3,{\mathit{B}}_3,{\mathit{C}}_3,\overline{{\mathit{G}}_1},{\mathit{G}}_{1,e}\in {\mathfrak{R}}^{3\times 3}$; ${\mathit{B}}_1\in {\mathfrak{R}}^{3\times 1}$; and ${\mathit{C}}_3,{\overline{\mathit{G}}}_3,{\mathit{G}}_{3,e}\in {\mathfrak{R}}^{3\times 13}$.

According to Assumption 3, $\Delta {\mathit{x}}_2^{\prime }$ changes much faster than $\Delta {\mathit{x}}_1$, so the latter can be neglected here. Because we set ${\beta }^c=0$, with sufficiently high sample rate, both $\Delta \beta$ and $\Delta {\mathit{d}}_1$ can be neglected as well. It should be noted that because the control surface obviously change faster than the attitude angle, $\Delta \delta$ cannot be neglected here. Therefore, we modify Equation (21) as:

$$\begin{array}{cc}\hfill {\dot{\mathit{x}}}_1& ={\dot{\mathit{x}}}_{1,0}+{\mathit{C}}_1\Delta \delta +(\overline{{\mathit{G}}_1}+{\mathit{G}}_{1,e})\Delta {\mathit{x}}_2^{\prime }\hfill \end{array}$$

(23)

For the same reason, $\Delta {\mathit{x}}_1$, $\Delta {\mathit{x}}_2$, and $\Delta {\mathit{x}}_3$ can also be neglected in Equation (22), so it can be rewritten as:

$${\dot{\mathit{x}}}_3={\dot{\mathit{x}}}_{3,0}+(\overline{{\mathit{G}}_3}+{\mathit{G}}_{3,e})\Delta \delta$$

(24)

The objective of control is to design a control law of thrust and $\delta$ such that the output ${\mathit{x}}_0$ can track ${\mathit{x}}_0^c$, where ${\mathit{x}}_0^c$, ${\dot{\mathit{x}}}_0^c$, and ${\ddot{\mathit{x}}}_0^c$ are bounded.

The structure of the control system is shown in the Figure 1. For convenience, we denote ${·}^c$ as the command signal that was given to controllers to follow, and ${·}^d$ as the virtual control input produced by the controller.

The IBS is applied in the outer-loop control, and it follows the trajectory command $x_0^c$ and produces thrust command T and virtual control input ${\alpha }^d$, ${\mu }^d$ as the desired attitude signal, while the SE is used to limit ${\alpha }^d$, ${\mu }^d$, T, and ${\dot{\alpha }}^d$, ${\dot{\mu }}^d$, $\dot{T}$. The inner-loop control consists of attitude control and angular rate control, where BS and IB are applied, respectively. The attitude control loop is to follow the attitude signals given in the velocity control loop and provides the desired angular rate signal $x_3^d$. Finally, the objective of the angular rate control loop follows $x_3^d$ and generates the effector deflection command $\delta$.

Prior to controller design, the following assumptions are made:

Assumption 2.

Accurate knowledge of the system states ${\mathit{x}}_0$, ${\mathit{x}}_1$, ${\mathit{x}}_2$, and ${\mathit{x}}_3$ is accessible.

Assumption 3.

The model uncertainties ${\mathit{d}}_1$, ${\mathit{d}}_2$, and ${\mathit{d}}_3$ are time-varying unknown yet bounded variables.

Assumption 4.

In the design of the outer control loops, it is believed that, with proper design of the inner control loops, and within the amplitude and rate limits of the virtual control input, instantaneous response to ${\mathit{x}}_2^{\prime }$ is assumed. This means the difference between ${\mathit{x}}_2^{\prime }$ and ${\mathit{x}}_2^{\prime c}$ is negligible [22].

Assumption 5.

The aerodynamic data we used are precise enough, so that $\parallel {\mathit{G}}_1{\overline{\mathit{G}}}_1^{-1}-\mathit{I}\parallel ⩽b_1⩽1$ and $\parallel {\mathit{G}}_3{\overline{\mathit{G}}}_3^{-1}-\mathit{I}\parallel ⩽b_3⩽1$, where ${\overline{\mathit{G}}}_3^{-1}$ is the generalized inversion of ${\overline{\mathit{G}}}_3$ [35].

Remark 3.

In practice, the control effectiveness is obtained through numerical differentiation of aerodynamic data. For example:

$${\overline{\mathit{G}}}_1=\frac{\partial {\overline{\mathit{g}}}_1}{\partial \beta }=\frac{{\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta +ϵ,\delta ,{\mathit{x}}_2^{\prime })-{\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta -ϵ,\delta ,{\mathit{x}}_2^{\prime })}{2ϵ}$$

(25)

where ${\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta +ϵ,\delta ,{\mathit{x}}_2^{\prime })$ and ${\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta -ϵ,\delta ,{\mathit{x}}_2^{\prime })$ is obtain through interpolating of discreet aerodynamic data, and ϵ is small positive number. Note that the aerodynamic data are essentially static and are well tailored in advance, so the numerical differentiation used here would not be affected by disturbances or local chatter.

The accuracy of the control effectiveness is affected by two factors, namely the quality of the aerodynamic data and the interpolation method. More detailed wind tunnel tests and more advanced interpolation algorithms could help improve both the accuracy of the control effectiveness. In theory, if the wind tunnel test is accurate enough, ${\overline{\mathit{G}}}_1$ can infinitely approach ${\mathit{G}}_1$. However, as this paper is controller design oriented, the wind tunnel test is not taken into account, so Assumption 5 is made.

According to [36], zero order of continuity would be enough for the controller design, so only linear interpolation is used in this paper. In fact, a certain degree of error in controller effectiveness will not affect the stability of the system, which will be discussed subsequently.

4. Main Result

In this section, the control strategies of the outer control loops and inner control loops are introduced. The detail of the design procedure are shown as follows.

4.1. Outer Control Loops

Based on a slide-mode incremental backstepping technique and a novel two-stage SE scheme, an adaptive controller is constructed for outer control loops. The design procedure is divided into two steps. In step one, we introduce the nominal adaptive slide-mode incremental backstepping controller, and then a two-stage SE scheme design is presented in step two.

4.1.1. Nominal Sliding-Mode Incremental Backstepping Controller Design

First, we choose the sliding surface $\mathit{s}$, where

$$\mathit{s}={\mathit{c}}_0{\mathit{e}}_0+{\dot{\mathit{e}}}_0$$

(26)

where ${\mathit{e}}_0={\mathit{x}}_0-{\mathit{x}}_0^c$ is the position error, and ${\mathit{c}}_0$ is a positive definite diagonal matrix.

Differentiating $\mathit{s}$ and using Equation (11), we get:

$$\begin{array}{c}\hfill \dot{\mathit{s}}={\mathit{c}}_0({\dot{\mathit{x}}}_0-{\dot{\mathit{x}}}_0^c)+({\ddot{\mathit{x}}}_0-{\ddot{\mathit{x}}}_0^c)={\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c)+({\dot{\mathit{x}}}_1-{\ddot{\mathit{x}}}_0^c)\end{array}$$

(27)

Substituting Equation (23) into Equation (27), we get:

$$\dot{\mathit{s}}={\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c)-{\ddot{\mathit{x}}}_0^c+{\mathit{f}}_1\left({\mathit{x}}_1\right)+\overline{{\mathit{g}}_1}({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })+{\xi }_1+{\overline{\mathit{G}}}_1\Delta {\mathit{x}}_2^{\prime }$$

(28)

where ${\xi }_1={\mathit{d}}_1+{\mathit{C}}_1\Delta \delta +{\mathit{G}}_{1,e}\Delta {\mathit{x}}_2^{\prime }$.

Remark 4.

Note that the control surface is designed to produce aerodynamic moments. Compared with the attitude angles, the contribution of control surfaces to aerodynamic forces is relatively small, which means ${\mathit{C}}_1$ is smaller than ${\mathit{G}}_1$. In addition, in small time increments, $\Delta \delta$ is also bounded. Thus, it is believed that ${\mathit{C}}_1\Delta \delta$ is bounded.

According to Propositions 1 and 2 in [8,35,37], if the sample rate is sufficiently high, and $\parallel {\mathit{G}}_1{\overline{\mathit{G}}}_1^{-1}-\mathit{I}\parallel ⩽b⩽1$, ${\mathit{G}}_{1,e}\Delta {\mathit{x}}_2^{\prime }$ would be globally bounded. Considering the condition is satisfied by Assumption 5, we know ${\mathit{G}}_{1,e}\Delta {\mathit{x}}_2^{\prime }$ is bounded. Together with Assumption 3, it is concluded that $\xi ={\mathit{d}}_1+{\mathit{C}}_1\Delta \delta +{\mathit{G}}_{1,e}\Delta {\mathit{x}}_2^{\prime }$ is bounded, that is, $|{\xi }_1|<{\overline{\xi }}_1$.

The desired incremental virtual control input is designed as:

$$\Delta {\mathit{x}}_2^{\prime d}={\overline{\mathit{G}}}^{-1}(-{\mathit{k}}_{\mathit{s},\mathbf{1}}\mathit{s}-{\mathit{k}}_{\mathit{s},\mathbf{2}}\mathit{sgn}\left(\mathit{s}\right)-\mathit{msgn}\left(\mathit{s}\right){\widehat{\xi }}_1-{\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c)+{\ddot{\mathit{x}}}_0^c-{\mathit{f}}_1\left({\mathit{x}}_1\right)-{\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta ,{\mathit{x}}_2^{\prime }))$$

(29)

where $\overline{\mathit{G}}=\partial {\mathit{g}}_1/\partial {\mathit{x}}_2^{\prime }$, $\mathit{sgn}\left(\mathit{s}\right)={[sgn\left(s_1\right),sgn\left(s_2\right),sgn\left(s_3\right)]}^T$ for $\mathit{s}={\left[s_1{s}_2{s}_3\right]}^T$, $\mathit{msgn}\left(\mathit{s}\right)=diag\left(\mathit{sgn}\right(s\left)\right)$, and ${\widehat{\xi }}_1$ is the adaptive compensation ${\xi }_1$. The adaptive law is designed as follows:

$${\dot{\widehat{\xi }}}_1=-{\mathit{c}}_{\xi 1}({\widehat{\xi }}_1+\left|\mathit{s}\right|)$$

(30)

where ${\mathit{c}}_d$ is a positive defined symmetric matrix; with sufficiently high sample rate, the derivative of virtual control input can be approximated as:

$${\dot{\mathit{x}}}_2^{\prime d}\approx \frac{\Delta {\mathit{x}}_2^{\prime c}}{\Delta t}$$

(31)

Integrating ${\dot{\mathit{x}}}_2^{\prime d}$, ${\mathit{x}}_2^{\prime d}$ obtains the desired signal of ${\mathit{x}}_2^{\prime }$.

4.1.2. Stability Enhancer

For the inner control loops, the instantaneous action of the actuator is assumed. For the outer control loops, however, such a simple assumption is inappropriate, because the dynamic of the attitude angle is obviously non-negligible. To meet the condition of Assumption 4, both amplitude and rate of attitude angle should be limited. Therefore, a two-layer stability enhancer is introduced.

As noted in Section 2, ${\mathit{x}}_2^{\prime c}$ and ${\dot{\mathit{x}}}_2^{\prime c}$ should be strictly limited in the ranges $[\underline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }},\overline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }}]$ and $[\underline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{s}}^{\prime }},\overline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{s}}^{\prime }}]$. For prevention, we set a virtual edge $[\underline{{\mathit{x}}_{\mathbf{2}\mathit{v}}^{\prime }},\overline{{\mathit{x}}_{\mathbf{2}\mathit{v}}^{\prime }}]$ and $[\underline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{v}}^{\prime }},\overline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{v}}^{\prime }}]$, so some action can be taken before the virtual control input reaches the strict edge.

The multilayer PHC scheme works as follows: when both ${\mathit{x}}_2^{\prime c}$ and ${\dot{\mathit{x}}}_2^{\prime c}$ are within the virtual edge, no stability enhancer is activated, so the nominal controller works; when ${\mathit{x}}_2^{\prime c}$ or ${\dot{\mathit{x}}}_2^{\prime c}$ exceeds the virtual edge, the first layer SE is activated to modify $\mathit{s}$; when ${\mathit{x}}_2^{\prime c}$ or ${\dot{\mathit{x}}}_2^{\prime c}$ exceeds the strict edge, the first and second layer SE work simultaneously, where the second layer adjusts the gain to force the control input to reduce effectively.

Before elaborating on the two-layer SE, the definitions of the saturation functions are given. First, we define the scalar rate saturation function as:

$$sa{t}_{r,i}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\left(i\right)\right)=\left\{\begin{array}{c}\begin{array}{cc}\hfill \overline{{\dot{\mathit{x}}}_{\mathbf{2}★}^{\prime }}\left(i\right),& {\dot{\mathit{x}}}_2^d\left(i\right)>\overline{{\dot{\mathit{x}}}_{\mathbf{2}★}^{\prime }}\left(i\right)\hfill \\ \hfill {\dot{\mathit{x}}}_2^d\left(i\right),& \overline{{\dot{\mathit{x}}}_{\mathbf{2}★}^{\prime }}\left(i\right)>{\dot{\mathit{x}}}_2^d\left(i\right)>\underline{{\dot{\mathit{x}}}_{\mathbf{2}★}^{\prime }}\left(i\right)\hfill \\ \hfill \underline{{\dot{\mathit{x}}}_{\mathbf{2}★}^{\prime }}\left(i\right),& \underline{{\dot{\mathit{x}}}_{\mathbf{2}★}^{\prime }}\left(i\right)>{\dot{\mathit{x}}}_2^d\left(i\right)\hfill \end{array}\hfill \end{array}\right.(i=1,2,3)$$

(32)

where ★ could be v or s and represent the virtual or strict edge.

Based on the above definition, the control input vector rate saturation function is given as:

$${\mathit{sat}}_{\mathit{v}}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\right)={\left[sa{t}_{v,1}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\left(1\right)\right)sa{t}_{v,2}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\left(2\right)\right)sa{t}_{v,3}^{★}\left({\dot{x}}_2^{\prime d}\left(3\right)\right)\right]}^T$$

(33)

In the same way, the control input vector amplitude saturation function is defined as:

$${\mathit{sat}}_{\mathit{a}}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\right)={\left[sa{t}_{a,1}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\left(1\right)\right)sa{t}_{a,2}^{★}\left({\dot{\mathit{x}}}_2^{\prime d}\left(2\right)\right)sa{t}_{a,3}^{★}\left({\dot{x}}_2^{\prime d}\left(3\right)\right)\right]}^T$$

(34)

In the following, the details of the first and second layers of SE are introduced.

First layer SE

When ${\mathit{x}}_2^{\prime d}$ or ${\dot{\mathit{x}}}_2^{\prime d}$ exceeds the virtual edge but is within the strict edge, the first layer SE is activated.

The first layer SE modifies the position error as:

$${\mathit{e}}_0={\mathit{x}}_0-{\mathit{x}}_0^c-{\mu }_1$$

(35)

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{\mu }}_1& =-{\mathit{c}}_{p1}{\mu }_1+{\mu }_{\mathbf{2}}\hfill \\ \hfill {\dot{\mu }}_2& =-{\mathit{c}}_{p2}{\mu }_2+{\overline{\mathit{G}}}_1({\tilde{\mathit{x}}}_{2,v}^{\prime d}+{\widetilde{\dot{\mathit{x}}}}_{2,v}^{\prime d})\hfill \end{array}\hfill \end{array}\right.$$

(36)

where ${\tilde{\mathit{x}}}_{2,v}^{\prime d}={\mathit{x}}_2^{\prime d}-{\mathit{sat}}_{\mathit{a}}^{\mathit{v}}\left({\mathit{x}}_{2,i}^{\prime d}\right)$, ${\widetilde{\dot{\mathit{x}}}}_{2,v}^{\prime d}={\dot{\mathit{x}}}_2^{\prime d}-{\mathit{sat}}_{\mathit{r}}^{\mathit{v}}\left({\dot{\mathit{x}}}_{2,i}^{\prime d}\right)$, ${\mathit{c}}_1$ and ${\mathit{c}}_2$ are the designed positive define matrix.

Differentiating $\mathit{s}$, we get:

$$\dot{\mathit{s}}={\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c-{\dot{\mu }}_1)-{\ddot{\mathit{x}}}_0^c+{\mathit{f}}_1\left({\mathit{x}}_1\right)+\overline{{\mathit{g}}_1}({\mathit{x}}_1,\beta ,\delta ,{\mathit{x}}_2^{\prime })+{\xi }_1-{\ddot{\mu }}_1+{\overline{\mathit{G}}}_1\Delta {\mathit{x}}_2^{\prime }$$

(37)

The modified desired incremental virtual control input with one layer SE is designed as:

$$\Delta {\overline{\mathit{x}}}_2^{\prime d}={\overline{\mathit{G}}}^{-1}(-{\mathit{k}}_{\mathit{s},\mathbf{1}}\mathit{s}-{\mathit{k}}_{\mathit{s},\mathbf{2}}\mathit{sgn}\left(\mathit{s}\right)-\mathit{msgn}\left(\mathit{s}\right){\widehat{\xi }}_1-{\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c-{\dot{\mu }}_1)+{\ddot{\mathit{x}}}_0^c-{\mathit{f}}_1\left({\mathit{x}}_1\right)-{\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta ,{\mathit{x}}_2^{\prime })+{\ddot{\mu }}_1)$$

(38)

Then, the derivative of virtual control input is modified as:

$${\dot{\overline{\mathit{x}}}}_2^{\prime d}\approx \frac{\Delta {\overline{\mathit{x}}}_2^{\prime d}}{\Delta t}$$

(39)

Integrating ${\dot{\overline{\mathit{x}}}}_2^{\prime d}$, ${\overline{\mathit{x}}}_2^{\prime d}$ obtains the desired signal of ${\mathit{x}}_2^{\prime }$.

Second layer SE

When ${\mathit{x}}_2^{\prime d}$ or ${\dot{\mathit{x}}}_2^{\prime d}$ exceeds the strict edge, the second layer SE is activated. Before introducing the second layer SE, a vector function is defined as:

$$\mathbf{tanh}\left(\mathit{x}\right)=[tanh\left(x_1\right){[tanh\left(x_2\right)\cdots tanh\left(x_n\right)]}^T,\mathit{x}={[x_1{x}_2\cdots x_n]}^T$$

(40)

where $x_i\in R$.

In order to make the control input rapidly decrease when it exceeds the strict edge, the second layer SE is designed as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \dot{\mathbf{\upsilon }}& ={\mathit{c}}_{p3}\mathbf{tanh}({\mathit{c}}_{p4}{\tilde{\mathit{x}}}_{2,s}^{\prime d}+{\mathit{c}}_{p5}{\widetilde{\dot{\mathit{x}}}}_{2,s}^{\prime d}-\mathbf{\upsilon })\hfill \\ \hfill \mathbf{\nu }& =\frac{1}{2}diag(\mathbf{tanh}\left(\mathbf{\upsilon }\right)+{\left[111\right]}^T)\hfill \end{array}\hfill \end{array}\right.$$

(41)

where ${\tilde{\mathit{x}}}_{2,s}^{\prime d}={\mathit{x}}_2^{\prime d}-{\mathit{sat}}_{\mathit{a}}^{\mathit{s}}\left({\mathit{x}}_{2,i}^{\prime d}\right)$, ${\widetilde{\dot{\mathit{x}}}}_{2,s}^{\prime d}={\dot{\mathit{x}}}_2^{\prime d}-{\mathit{sat}}_{\mathit{r}}^{\mathit{s}}\left({\dot{\mathit{x}}}_{2,i}^{\prime d}\right)$, ${\mathit{c}}_3$, ${\mathit{c}}_4$, and ${\mathit{c}}_5$ are the designed matrix.

From Equation (41) we could know that, when ${\mathit{x}}_2^{\prime d}$ or ${\dot{\mathit{x}}}_2^{\prime d}$ exceeds the strict edge, $\mu$ will rise form 0 to 1 rapidly, and when ${\mathit{x}}_2^{\prime d}$ or ${\dot{\mathit{x}}}_2^{\prime d}$ drops within the strict edge, $\mu$ goes back to 0. Thus, the incremental virtual control input with two-layer SE is modified as:

$$\Delta {\overline{\overline{\mathit{x}}}}_2^{\prime d}={\overline{\mathit{G}}}^{-1}(-({\mathit{k}}_{\mathit{s},\mathbf{1}}-{\mathit{k}}_{\mathit{s},\mathbf{3}}\mathbf{\nu })\mathit{s}-({\mathit{k}}_{\mathit{s},\mathbf{2}}-{\mathit{k}}_{\mathit{s},\mathbf{4}}\mathbf{\nu })\mathit{sgn}\left(\mathit{s}\right)-\mathit{sgn}\left(\mathit{s}\right){\widehat{\xi }}_1-{\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c-{\dot{\mathbf{\mu }}}_1)+{\ddot{\mathit{x}}}_0^c-{\mathit{f}}_1\left({\mathit{x}}_1\right)-{\overline{\mathit{g}}}_1({\mathit{x}}_1,\beta ,{\mathit{x}}_2^{\prime })+{\ddot{\mathbf{\mu }}}_1)$$

(42)

Then, with sufficiently high sample rate, the derivative of virtual control input is modified as:

$${\dot{\overline{\overline{\mathit{x}}}}}_2^{\prime d}\approx \frac{\Delta {\overline{\overline{\mathit{x}}}}_2^{\prime d}}{\Delta t}$$

(43)

Integrating ${\dot{\overline{\overline{\mathit{x}}}}}_2^{\prime d}$, ${\overline{\overline{\mathit{x}}}}_2^{\prime d}$ obtains the desired signal of ${\mathit{x}}_2^{\prime }$.

Remark 5.

Considering the controller in the form of ${\dot{\overline{\overline{\mathit{x}}}}}_2^{\prime d}$, it is concluded that, if the strict edge is not exceeded, ν would stabilize at 0, so the ${\dot{\overline{\overline{\mathit{x}}}}}_2^{\prime d}$ turns to ${\dot{\overline{\mathit{x}}}}_2^{\prime d}$. For the same reason, ${\dot{\overline{\mathit{x}}}}_2^{\prime d}$ would turn to ${\dot{\mathit{x}}}_2^{\prime d}$ if the virtual edge is not exceeded, as ${\mathbf{\lambda }}_{\mathbf{1}}$ and ${\mathbf{\lambda }}_{\mathbf{2}}$ would also stabilize at 0 in this case. Therefore, the controller ${\overline{\overline{\mathit{x}}}}_2^{\prime d}$ will be equivalent to ${\overline{\mathit{x}}}_2^{\prime d}$ and ${\mathit{x}}_2^{\prime d}$. In conclusion, the controller ${\overline{\overline{\mathit{x}}}}_2^{\prime d}$ is capable of dealing with all situations.

Remark 6.

Because ${\overline{\overline{\mathit{x}}}}_2^{\prime d}={\left[{\overline{\overline{T}}}_f^d{\overline{\overline{\alpha }}}^d{\overline{\overline{\mu }}}^d\right]}^T$, ${\overline{\overline{T}}}_f^d$ is the desired signal of thrust that would be directly sent to engine, whereas ${\overline{\overline{\alpha }}}^d$ and ${\overline{\overline{\mu }}}^d$ are the virtual control inputs of the inner control loops, together with ${\beta }^c=0$, define ${\mathit{x}}_2^d=\left[{\overline{\overline{\mu }}}^d{\overline{\overline{\alpha }}}^d{\beta }^c\right]$ as the virtual control input vector of inner control loops.

4.1.3. Stability Analysis of Outer Control Loops

Theorem 1.

Consider the outer control loops, satisfying Assumptions 1–4, if the conditions

$$\begin{array}{c}{\lambda }_{min}({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3})>\frac{1}{2}\end{array}$$

(44)

$$\begin{array}{c}{\lambda }_{min}\left({\mathit{c}}_{\xi 1}\right)>1\end{array}$$

(45)

$$\begin{array}{c}{\lambda }_{min}({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4})>0\end{array}$$

(46)

are satisfied, the proposed adaptive controller together with a two-layer SE scheme guarantees the asymptotic stability of system.

Proof. 

Select the Lyapunov function candidate for the outer control loops as:

$$V_s=\frac{1}{2}{\mathit{s}}^T\mathit{s}+\frac{1}{2}{\tilde{\mathbf{\xi }}}_1^T{\mathit{c}}_{\xi 1}^{-1}{\tilde{\mathbf{\xi }}}_1$$

(47)

where and ${\tilde{\mathbf{\xi }}}_1={\widehat{\mathbf{\xi }}}_1-{\overline{\mathbf{\xi }}}_1$.

Considering the first layer SE and adaptive law, according to Equations (30) and (37) and Remark 4, and taking the derivation of $V_s$ we obtain:

$$\begin{array}{cc}\hfill {\dot{V}}_s& ={\mathit{s}}^T({\mathit{c}}_0({\mathit{x}}_1-{\dot{\mathit{x}}}_0^c-{\dot{\mathbf{\mu }}}_1)-{\ddot{\mathit{x}}}_0^c+{\mathit{f}}_1\left({\mathit{x}}_1\right)+\overline{{\mathit{g}}_1}({\mathit{x}}_1,\beta ,\mathbf{\delta },{\mathit{x}}_2^{\prime })+{\mathbf{\xi }}_1-{\ddot{\mathbf{\mu }}}_1+{\overline{\mathit{G}}}_1\Delta {\mathit{x}}_2^{\prime })+{\tilde{\mathbf{\xi }}}_1^T(-{\widehat{\mathbf{\xi }}}_1+\left|\mathit{s}\right|)\hfill \end{array}$$

(48)

Considering Assumption 4, we have $\Delta {\mathit{x}}_2^{\prime }=\Delta {\mathit{x}}_2^{\prime c}$. Noticing Equations (23) and (42) and Remark 4, we get:

$$\begin{array}{c}\begin{array}{cc}\hfill {\dot{V}}_s& ={\mathit{s}}^T(-({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3}\mathsf{\nu })\mathit{s}-({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4}\mathsf{\nu })\mathit{sgn}\left(\mathit{s}\right)+{\mathsf{\xi }}_1-\mathit{msgn}\left(\mathit{s}\right){\widehat{\mathsf{\xi }}}_1)+{\tilde{\mathsf{\xi }}}_1^T(-{\widehat{\mathsf{\xi }}}_1+\left|\mathit{s}\right|)\hfill \\ \hfill & ⩽{\mathit{s}}^T(-({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3}\mathsf{\nu })\mathit{s}-({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4}\mathsf{\nu })\mathit{sgn}\left(\mathit{s}\right)+\mathit{msgn}\left(\mathit{s}\right){\overline{\mathsf{\xi }}}_1-\mathit{msgn}\left(\mathit{s}\right){\widehat{\mathsf{\xi }}}_1)+{\tilde{\mathsf{\xi }}}_1^T(-{\widehat{\mathsf{\xi }}}_1+\left|\mathit{s}\right|)\hfill \\ \hfill & ⩽{\mathit{s}}^T(-({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3}\mathsf{\nu })\mathit{s}-({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4}\mathsf{\nu })\mathit{sgn}\left(\mathit{s}\right))+\mathit{msgn}\left(\mathit{s}\right){\mathit{s}}^T(\overline{\mathsf{\xi }}-{\widehat{\mathsf{\xi }}}_1)+{\tilde{\mathsf{\xi }}}_1^T(-{\widehat{\mathsf{\xi }}}_1+\left|\mathit{s}\right|)\hfill \\ \hfill & ={\mathit{s}}^T(-({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3}\mathsf{\nu })\mathit{s}-({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4}\mathsf{\nu })\mathit{sgn}\left(\mathit{s}\right))-{\tilde{\mathsf{\xi }}}_1^T{\widehat{\mathsf{\xi }}}_1\hfill \\ \hfill & ⩽{\mathit{s}}^T(-({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3}\mathsf{\nu })\mathit{s}-({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4}\mathsf{\nu })\mathit{sgn}\left(\mathit{s}\right))-\frac{1}{2}{\widehat{\mathsf{\xi }}}_1^T{\widehat{\mathsf{\xi }}}_1+\frac{1}{2}{\overline{\mathsf{\xi }}}_1^T{\overline{\mathsf{\xi }}}_1\hfill \end{array}\hfill \end{array}$$

(49)

According to Equation (41), we have $0<|{\lambda }_{max}\left(\mathsf{\nu }\right)|<1$, and since ${\lambda }_{min}\left({\mathit{c}}_{\mathsf{\xi }1}\right)>1$, we have ${\lambda }_{max}\left({\mathit{c}}_{\mathsf{\xi }1}^{-1}\right)>1$, therefore:

$$\begin{array}{c}\hfill {\dot{V}}_s⩽-{\mathit{s}}^T({\mathit{k}}_{s,1}-{\mathit{k}}_{s,3})\mathit{s}-{\mathit{s}}^T({\mathit{k}}_{s,2}-{\mathit{k}}_{s,4})\mathit{sgn}\left(\mathit{s}\right)-\frac{1}{2}{\tilde{\mathsf{\xi }}}_1^T{\mathit{c}}_{\mathsf{\xi }1}{\tilde{\mathsf{\xi }}}_1+\frac{1}{2}{\overline{\mathsf{\xi }}}_1^T{\mathit{c}}_{\mathsf{\xi }1}{\overline{\mathsf{\xi }}}_1\end{array}$$

(50)

Considering Equations (44)–(46), we have:

$$\begin{array}{cc}\hfill {\dot{V}}_s& ⩽-\frac{1}{2}{\mathit{s}}^T\mathit{s}-\frac{1}{2}{\tilde{\mathsf{\xi }}}_1^T{\tilde{\mathsf{\xi }}}_1+{\sigma }_s\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill & =-V_s+{\sigma }_s\hfill \end{array}$$

(52)

where ${\sigma }_s=\frac{1}{2}{\overline{\mathsf{\xi }}}_1^T{\overline{\mathsf{\xi }}}_1$, then we obtain

$$0\le V_s\left(t\right)\le {\sigma }_s+[V_s\left(0\right)-{\sigma }_s]e^{-t}$$

(53)

Therefore, $V_s$ is bounded, and the stability of the outer-loop control is proved. □

4.2. Inner Control Loops

The objective of the inner control loops is to follow the virtual control input ${\mathit{x}}_2^d$. The design of the inner control loops is divided into two steps: the attitude control and the angular rate control. The attitude control follows ${\mathit{x}}_2^d$ and generates the desired angular rate ${\mathit{x}}_3^d$. Because there are no uncertainties in Equation (13), and control input is affine, BS is used to control the attitude angle. The angular rate control follows ${\mathit{x}}_3^d$ and generates the desired control surface deflection. Due to strong nonlinearity of the control surface and outside disturbances, IBS is used for angular rate control, with a adaptive model-based estimator to obtain the state derivation. The details of inner control loops are shown as follows:

4.2.1. Controller Design

Attitude control

Define the attitude error ${\mathit{e}}_2$ as:

$${\mathit{e}}_2={\mathit{x}}_2-{\mathit{x}}_2^d$$

(54)

Take the time derivation of ${\mathit{e}}_{\mathbf{2}}$ as:

$${\dot{\mathit{e}}}_2={\overline{\mathit{f}}}_2({\mathit{x}}_1,{\mathit{x}}_2,\mathsf{\delta })+{\mathsf{\xi }}_2+{\mathit{g}}_2·{\mathit{x}}_3$$

(55)

where ${\mathsf{\xi }}_2={\mathit{d}}_2-{\dot{\mathit{x}}}_2^d$.

According to Assumption 3, ${\mathit{f}}_{2,e}({\mathit{x}}_1,{\mathit{x}}_2,\mathsf{\delta })$ is bounded; because ${\dot{\mathit{x}}}_2^d$ is a continuous function of $\mathit{s}$, $\tilde{\mathsf{\xi }}$, ${\mathsf{\mu }}_1$, and $\mathsf{\nu }$, it is also bounded. Thus, it is concluded that ${\mathsf{\xi }}_2$ is bounded, that is, ${\mathsf{\xi }}_2<{\overline{\mathsf{\xi }}}_2$.

The virtual control input is designed as:

$${\mathit{x}}_3^d={\mathit{g}}_2^{-1}\left(-{\mathit{k}}_{\mathbf{2}}{\mathit{e}}_{\mathbf{2}}-\mathit{msgn}\left({\mathit{e}}_2\right){\widehat{\mathsf{\xi }}}_2-{\overline{\mathit{f}}}_2({\mathit{x}}_1,{\mathit{x}}_2,\mathsf{\delta })\right)$$

(56)

where ${\mathit{k}}_{\mathbf{2}}$ is a designed matrix, and ${\widehat{\mathsf{\xi }}}_2$ is the adaptive compensation for ${\mathsf{\xi }}_2$, which is designed as:

$${\dot{\widehat{\mathsf{\xi }}}}_2=-{\mathit{c}}_{\mathsf{\xi }2}({\widehat{\mathsf{\xi }}}_2+|{\mathit{e}}_2\left|\right)$$

(57)

where ${\tilde{\mathsf{\xi }}}_2={\widehat{\mathsf{\xi }}}_2-{\overline{\mathsf{\xi }}}_2$.

Angular rate control

Define the angular rate error as:

$${\mathit{e}}_3={\mathit{x}}_3-{\mathit{x}}_3^d$$

(58)

Taking the derivation of ${\mathit{e}}_3$ and writing it in incremental form according to Equations (20) and (24), we get:

$${\dot{\mathit{e}}}_3={\mathit{f}}_3\left({\mathit{x}}_3\right)+{\overline{\mathit{g}}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\mathsf{\delta })+{\mathsf{\xi }}_3+{\overline{\mathit{G}}}_3\Delta \mathsf{\delta }$$

(59)

where ${\mathsf{\xi }}_3={\mathit{d}}_3+{\mathit{G}}_{1,e}\Delta \mathsf{\delta }-{\dot{\mathit{x}}}_3^d$, for the same reason that ${\mathsf{\xi }}_1$ and ${\dot{\mathit{x}}}_2^d$ are bounded, ${\mathsf{\xi }}_3$ and ${\dot{\mathit{x}}}_3^d$ are bounded also, so we can conclude that there exists a 1*3 vector ${\overline{\mathsf{\xi }}}_3={\left[{\xi }_{3,1}{\xi }_{3,2}{\xi }_{3,3}\right]}^T$ where ${\xi }_{3,1},{\xi }_{3,2},{\xi }_{3,3}\in {\mathfrak{R}}^{+}$, let $|{\mathsf{\xi }}_3|<{\overline{\mathsf{\xi }}}_3$.

The desired moment increment $\Delta {\mathit{M}}_d$ is designed as:

$$\Delta {\mathit{M}}_d=-{\mathit{k}}_{\mathbf{3}}{\mathit{e}}_{\mathbf{3}}-{\mathit{f}}_3\left({\mathit{x}}_3\right)-{\overline{\mathit{g}}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\mathsf{\delta })-\mathit{msgn}\left({\mathit{e}}_3\right){\widehat{\mathsf{\xi }}}_3-\frac{1}{{\tau }_3}({\mathit{x}}_3^d-{\mathbf{\alpha }}_3)$$

(60)

where ${\widehat{\mathsf{\xi }}}_3$ is the adaptive compensation of ${\mathsf{\xi }}_3$, which is designed as:

$${\dot{\widehat{\mathsf{\xi }}}}_3=-{\mathit{c}}_{\mathsf{\xi }3}({\widehat{\mathsf{\xi }}}_3+|{\mathit{e}}_3\left|\right)$$

(61)

where ${\tilde{\mathsf{\xi }}}_3={\mathsf{\xi }}_3-{\widehat{\mathsf{\xi }}}_3$.

Thus, the desired control surface deflection could be obtained as:

$$\Delta {\mathsf{\delta }}_d={\overline{\mathit{G}}}_3^{+}\Delta {\mathit{M}}_d$$

(62)

where ${\overline{\mathit{G}}}_3^{+}={\overline{\mathit{G}}}_3^T{\left({\overline{\mathit{G}}}_3{\overline{\mathit{G}}}_3^T\right)}^{-1}$ is a Moore–Penrose pseudoinversion of ${\overline{\mathit{G}}}_3$.

Remark 7.

The 11 independent control surfaces of the ICE aircraft makes it an over-actuated system. Mathematically, as ${\overline{\mathit{G}}}_3$ is not square and generally has a non-trivial null space, there are an infinite number of rudder surface deflection combinations that could meet specific moment commands. Therefore, some approaches must be taken to select the ’best’ among the infinite combinations. These approaches are known as the ’control allocation’ method. The most intuitive way to conduct control allocation is to introduce a secondary goal. For example, control allocation could be described as an optimization problem as follows:

$$\begin{array}{c}min\frac{1}{2}\Delta {\mathsf{\delta }}^T\Delta \mathsf{\delta }\\ s.t.{\overline{\mathit{G}}}_3\Delta \mathsf{\delta }={\mathit{M}}_d\end{array}$$

The unique solution of this problem is $\Delta {\mathsf{\delta }}_d={\overline{\mathit{G}}}_3^{+}{\mathit{M}}_d$, and ${\overline{\mathit{G}}}_3^{+}$ is known as the Moore–Penrose pseudoinversion [29] of ${\overline{\mathit{G}}}_3$. In the control allocation field, this method is known as generalized inverse, which not only ensures the uniqueness of the solution, but also brings additional benefits, that is, minimizing the increment of the control input. At present, is in-depth research on control allocation. Many advance control allocation methods, such as reinforcement learning control allocation and various kind of multi-objective optimization, are proposed. However, control allocation is not the focus of this paper. We consider the simplicity of the mathematical definition of the generalized inverse method and its successful application in practice. It is believed that the generalized inverse method is capable of dealing with all situations involved in this paper.

4.2.2. Stability Analysis of Inner Control Loops

Theorem 2.

Considering the inner control loops with Assumptions 1–5. If the conditions

$$\begin{array}{c}-{\lambda }_{min}\left({\mathit{k}}_2\right)+\frac{1}{2}⩽-\frac{1}{2}\end{array}$$

(63)

$$\begin{array}{c}-{\lambda }_{min}\left({\mathit{k}}_3\right)+\frac{1}{2}{\lambda }_{max}\left({\mathit{g}}_2^T{\mathit{g}}_2\right)⩽-\frac{1}{2}\end{array}$$

(64)

$$\begin{array}{c}{\lambda }_{min}\left({\mathit{c}}_{\mathsf{\xi }2}\right)⩾1\end{array}$$

(65)

$$\begin{array}{c}{\lambda }_{min}\left({\mathit{c}}_{\mathsf{\xi }3}\right)⩾1\end{array}$$

(66)

are satisfied, the proposed controller can guarantee the boundedness of all error signals.

Proof. 

Considering the attitude angle tracking error and filter error, define the candidate Lyapunov function as:

$$V=\frac{1}{2}{\mathit{e}}_2^T{\mathit{e}}_2+\frac{1}{2}{\tilde{\mathsf{\xi }}}_2^T{\mathit{c}}_{\mathsf{\xi }2}^{-1}{\tilde{\mathsf{\xi }}}_2+\frac{1}{2}{\mathit{e}}_3^T{\mathit{e}}_3+\frac{1}{2}{\tilde{\mathsf{\xi }}}_3^T{\mathit{c}}_{\mathsf{\xi }3}^{-1}{\tilde{\mathsf{\xi }}}_3$$

(67)

Take the time derivation of V:

$$\begin{array}{c}\begin{array}{cc}\hfill \dot{V}& ={\mathit{e}}_2^T{\dot{\mathit{e}}}_2+{\tilde{\mathsf{\xi }}}_2^T{\mathit{c}}_{\mathsf{\xi }2}^{-1}{\dot{\widehat{\mathsf{\xi }}}}_2+{\mathit{e}}_3^T{\dot{\mathit{e}}}_3+{\tilde{\mathsf{\xi }}}_3^T{\mathit{c}}_{\mathsf{\xi }3}^{-1}{\dot{\widehat{\mathsf{\xi }}}}_3\hfill \end{array}\hfill \end{array}$$

(68)

According to Equations (55)–(58), we have:

$$\begin{array}{c}\begin{array}{cc}\hfill {\mathit{e}}_2^T{\dot{\mathit{e}}}_2+{\tilde{\mathsf{\xi }}}_2^T{\dot{\widehat{\mathsf{\xi }}}}_2& ={\mathit{e}}_2^T({\overline{\mathit{f}}}_2({\mathit{x}}_1,{\mathit{x}}_2,\mathsf{\delta })+{\mathsf{\xi }}_2+{\mathit{g}}_2·({\mathit{x}}_3^d+{\mathit{e}}_3))+{\tilde{\mathsf{\xi }}}_2^T(-{\widehat{\mathsf{\xi }}}_2+|{\mathit{e}}_2\left|\right)\hfill \\ \hfill & ={\mathit{e}}_2^T(-{\mathit{k}}_2{\mathit{e}}_2+{\mathsf{\xi }}_2-\mathit{msgn}\left({\mathit{e}}_2\right){\widehat{\mathsf{\xi }}}_2+{\mathit{g}}_2{\mathit{e}}_3)-{\tilde{\mathsf{\xi }}}_2^T{\widehat{\mathsf{\xi }}}_2+{\tilde{\mathsf{\xi }}}_2^T\left|{\mathit{e}}_2\right|\hfill \\ \hfill & ⩽-{\mathit{e}}_2^T{\mathit{k}}_2{\mathit{e}}_2+|{\mathit{e}}_2^T|{\mathsf{\xi }}_2-|{\mathit{e}}_2^T|{\widehat{\mathsf{\xi }}}_2+{\mathit{e}}_2^T{\mathit{g}}_2{\mathit{e}}_3-{\tilde{\mathsf{\xi }}}_2^T{\widehat{\mathsf{\xi }}}_2+{\tilde{\mathsf{\xi }}}_2^T\left|{\mathit{e}}_2\right|\hfill \\ \hfill & ⩽-{\lambda }_{min}\left({\mathit{k}}_2\right){\mathit{e}}_2^T{\mathit{e}}_2+{\mathit{e}}_2^T{\mathit{g}}_2{\mathit{e}}_3-{\tilde{\mathsf{\xi }}}_2^T{\widehat{\mathsf{\xi }}}_2\hfill \\ \hfill & ⩽(-{\lambda }_{min}\left({\mathit{k}}_2\right)+\frac{1}{2}){\mathit{e}}_2^T{\mathit{e}}_2+\frac{1}{2}{\lambda }_{max}\left({\mathit{g}}_2^T{\mathit{g}}_2\right){\mathit{e}}_3^T{\mathit{e}}_3-\frac{1}{2}{\widehat{\mathsf{\xi }}}_2^T{\widehat{\mathsf{\xi }}}_2+\frac{1}{2}{\overline{\mathsf{\xi }}}_2^T{\overline{\mathsf{\xi }}}_2\hfill \end{array}\hfill \end{array}$$

(69)

According to Equations (59)–(61), we have:

$$\begin{array}{c}\begin{array}{cc}\hfill {\mathit{e}}_3^T{\dot{\mathit{e}}}_3+{\tilde{\mathsf{\xi }}}_3^T{\dot{\widehat{\mathsf{\xi }}}}_3& ={\mathit{e}}_3^T({\mathit{f}}_3\left({\mathit{x}}_3\right)+{\overline{\mathit{g}}}_3({\mathit{x}}_1,{\mathit{x}}_2,{\mathit{x}}_3,\mathsf{\delta })+{\mathsf{\xi }}_3+{\overline{\mathit{G}}}_3\Delta \mathsf{\delta })+{\tilde{\mathsf{\xi }}}_3^T(-{\widehat{\mathsf{\xi }}}_3+|{\mathit{e}}_3\left|\right)\hfill \\ \hfill & ={\mathit{e}}_3^T(-{\mathit{k}}_3{\mathit{e}}_3+{\mathsf{\xi }}_3-\mathit{msgn}\left({\mathit{e}}_3\right){\widehat{\mathsf{\xi }}}_3)-{\tilde{\mathsf{\xi }}}_3^T{\widehat{\mathsf{\xi }}}_3+{\tilde{\mathsf{\xi }}}_3^T\left|{\mathit{e}}_3\right|\hfill \\ \hfill & ⩽-{\lambda }_{min}\left({\mathit{k}}_3\right){\mathit{e}}_3^T{\mathit{e}}_3-\frac{1}{2}{\tilde{\mathsf{\xi }}}_3^T{\tilde{\mathsf{\xi }}}_3+\frac{1}{2}{\overline{\mathsf{\xi }}}_3^T{\overline{\mathsf{\xi }}}_3\hfill \end{array}\hfill \end{array}$$

(70)

Substituting Equations (69) and (70) into Equation (68), we get:

$$\dot{V}⩽\left(-{\lambda }_{min}\left({\mathit{k}}_2\right)+\frac{1}{2}\right){\mathit{e}}_2^T{\mathit{e}}_2-\frac{1}{2}{\tilde{\mathsf{\xi }}}_2^T{\tilde{\mathsf{\xi }}}_2+\left(-{\lambda }_{min}\left({\mathit{k}}_3\right)+\frac{1}{2}{\lambda }_{max}\left({\mathit{g}}_2^T{\mathit{g}}_2\right)\right){\mathit{e}}_3^T{\mathit{e}}_3-\frac{1}{2}{\tilde{\mathsf{\xi }}}_3^T{\tilde{\mathsf{\xi }}}_3+\frac{1}{2}{\overline{\mathsf{\xi }}}_2^T{\overline{\mathsf{\xi }}}_2+\frac{1}{2}{\overline{\mathsf{\xi }}}_3^T{\overline{\mathsf{\xi }}}_3$$

(71)

According to Equations (63)–(66), we have:

$$\begin{array}{cc}\hfill \dot{V}& ⩽-\frac{1}{2}{\mathit{e}}_2^T{\mathit{e}}_2-\frac{1}{2}{\tilde{\mathsf{\xi }}}_2^T{\tilde{\mathsf{\xi }}}_2-\frac{1}{2}{\mathit{e}}_3^T{\mathit{e}}_3-\frac{1}{2}{\tilde{\mathsf{\xi }}}_3^T{\tilde{\mathsf{\xi }}}_3+\sigma \hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill & =-V+\sigma \hfill \end{array}$$

(73)

where $\sigma =\frac{1}{2}{\overline{\mathsf{\xi }}}_2^T{\overline{\mathsf{\xi }}}_2+\frac{1}{2}{\overline{\mathsf{\xi }}}_3^T{\overline{\mathsf{\xi }}}_3$, then we obtain

$$0\le V\left(t\right)\le \sigma +[V\left(0\right)-\sigma ]e^{-t}$$

(74)

Therefore, V is bounded, and the stability of inner-loop control is proved. □

5. Simulation

The proposed incremental sliding-mode controller with two-layer SE has two main advantages: its robustness and the two-layer SE’s contribution to stability and state saturation. Thus, two simulations are presented in this section. The first simulation aims to verify the robustness of the proposed controller, and the second simulation aims to verify the effectiveness of two-layer SE.

The trajectory command signal is generated by the following system:

$${\dot{\mathit{x}}}_0^c=\left[\begin{array}{c}{V}^{c}cos{\gamma }^{c}cos{\chi }^c\\ V^{c}cos{\gamma }^{c}sin{\chi }^c\\ V^{c}sin{\gamma }^c\end{array}\right]$$

(75)

The trajectory command signal in 3-dimensional space is given in Figure 2; ${\gamma }^c$ and ${\chi }^c$ are given in Figure 3, and $X^c\left(0\right)=0$ ft, $Y^c\left(0\right)=0$, $Z^c\left(0\right)=5000$ ft. In addition, $V^c=627$ ft/s, ${\beta }^c=0$. This means that the aircraft should fly at a constant speed and perform coordinated turns.

The control parameters are chosen as: ${\mathit{c}}_0=diag(0.21,0.25,0.13)$, ${\mathit{k}}_{s,1}=diag(0.13,0.13,0.16)$, ${\mathit{k}}_{s,2}=diag(3.50,3.50,3.50)$, ${\mathit{k}}_{s,3}=diag(0.76,0.76,0.85)$, ${\mathit{k}}_{s,4}=diag(3.00,3.00,3.00)$, ${\mathit{c}}_{\mathsf{\xi }1}=diag(1.6,1.6,1.6)$, ${\mathit{c}}_{\mathsf{\xi }2}=diag(1.5,1.5,1.5)$, ${\mathit{c}}_{\mathsf{\xi }3}=diag(1.5,1.5,1.5)$, ${\mathit{c}}_{p1}=diag$$(0.50,0.50,0.50)$, ${\mathit{c}}_{p2}=diag(0.50,0.50,0.50)$, ${\mathit{c}}_{p3}=diag(1.50,1.50,1.50)$, ${\mathit{c}}_{p4}=diag(0.10,0.30,0.20)$, ${\mathit{c}}_{p5}=diag(0.10,0.10,0.10)$, ${\mathit{k}}_2=diag(1.5,2,1.5)$, and ${\mathit{k}}_3=diag(10,10,10)$. The strict edge is set to be $\underline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }}={[0-5-60]}^T$, $\overline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }}={\left[400004060\right]}^T$, $\underline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{s}}^{\prime }}={[-1000-0.1-0.1]}^T$, $\overline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{s}}^{\prime }}={\left[20000.10.1\right]}^T$, and the virtual edge is set to be $\underline{{\mathit{x}}_{\mathbf{2}\mathit{v}}^{\prime }}={[-2000-4-50]}^T$, $\overline{{\mathit{x}}_{\mathbf{2}\mathit{s}}^{\prime }}={\left[390003550\right]}^T$, $\underline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{s}}^{\prime }}={[-1800-0.05-0.05]}^T$, $\overline{{\dot{\mathit{x}}}_{\mathbf{2}\mathit{s}}^{\prime }}={\left[18000.050.05\right]}^T$.

The initial trim conditions are selected as: $X_0=0$, $Y_0=0$, $Z_0=5000$ ft, $v_{x0}=627$ ft/s, $v_{y0}=0$, $v_{z0}=0$, ${\beta }_0=-0.0196$ deg, ${\alpha }_0=3.76$ deg, $\mathsf{\mu }=0.304$ deg, $p=0$, $q=0$, $r=0$.

5.1. Simulation 1

In the first simulation, the performance of the controller in the presence of disturbance and without disturbance is tested. The disturbances here consists of aerodynamic coefficient errors and external disturbances. Aerodynamic coefficient errors come from different interpolation methods. Note, the same aerodynamic data are used in the aircraft model build and controller design; however, to ensure high fidelity, an Akima spline is used to build the aircraft model, whereas linear point slope is used in the controller design. Different interpolation methods can cause errors up to 30 percent, as shown in Figure 4. The external disturbances are designed as:

$${\mathit{d}}_F=\left[\begin{array}{cccc}& 0.25sin\left(0.02t\right)& 0& 0\\ & 0& 0.25sin(0.02t+\pi /4)& 0\\ & 0& 0& 0.25sin(0.02t+\pi /2)\end{array}\right]·{\mathit{F}}_a·{\mathit{T}}_{\mathit{ve}}$$

(76)

$${\mathit{d}}_M=\left[\begin{array}{cccc}& 0.3sin(0.2t++\pi /5)& 0& 0\\ & 0& 0.3sin(0.2t+\pi /4)& 0\\ & 0& 0& 0.3sin(0.2t+\pi /2)\end{array}\right]·{\mathit{M}}_a$$

(77)

The results of simulation 1 are shown in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14. It can be observed from Figure 5 and Figure 6 that the presence of disturbances do not significantly affect the trajectory tracking performance of the controller. From Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12, we know that the overshoot and regulation time of the rest of the signals are not affected very much, either. The robustness of the proposed controller is satisfactory.

5.2. Simulation 2

In the second simulation, the performance of the controller without two-layer SE is tested, and the results are shown in Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22.

From Figure 15 and Figure 16, we know that, in the absence of SE, the trajectory tracking error does not converge, and the aircraft ends up losing control.

From Figure 18 and Figure 20, we know the inner-loop control works well, so the problems come from the outer-loop control. From Figure 18, Figure 19, and Figure 21, we see that, compared with the controller with SE, the absence of SE causes saturation in $\alpha$, $\mathsf{\mu }$, and T, respectively. In Figure 21, T reaches the amplitude limits around 100 s, but the aircraft is still in control. It appears that the amplitude saturation endures to some degree, and the direct reason the aircraft loses control is that the virtual control input reaches the rate limit around 175 s, as shown in Figure 19. To conclude, the system is more sensitive to rate saturation.

From a theoretical point of view, if the amplitude of the virtual control input is saturated, at least the control direction is still correct, and the system will still remain stable while the saturation is not too severe. The reason why the system is more sensitive to rate saturation is that the stability of the system is based on Assumption 4. If the rate of the virtual control input is saturated, it means that the virtual control input cannot be responded to quickly enough, so Assumption 4 will no longer hold. Thereby, the virtual control input may severely chatter, and the system would end up losing control.

In this simulation, we can found that the proposed SE can compensate for the saturation and make the system more stable. Certainly, if saturation is not too severe, the system may remain stable without SE. However, the absence of SE will make the time delay of the inner-loop state become unbearable if the aircraft is at the edge of the flight envelope. The results of this simulation demonstrate the effectiveness of the SE.

6. Conclusions

This paper proposes a trajectory controller for a tailless aircraft. The trajectory controller is divided into three parts: outer-loop control, attitude control, and angular rate control. The incremental backstepping sliding-mode control and incremental backstepping control are applied in the outer-loop control and angular rate control, and the attitude control loops are treated by the backstepping approach. A two-layer stability enhancer that considers both amplitude and rate limits of the virtual control input is proposed. With the help of SE, the incremental control method could be extended to outer-loop control. An adaptive estimator for the state derivative is proposed, and, together with the incremental control method, the controller shows excellent robustness. Using Lyapunov theory, the system is proven to be globally asymptotic stable. In the simulations, the robustness of the controller and the effectiveness of SE are proven. Due to disturbances, sensor errors, and real-world bias, verifying the result in an experiment would be slightly more difficult than in simulation. It also should be noted that, because the aerodynamic data in the simulation are high-fidelity, we believe the result would still be well supported. Further investigation could be devoted to studying if there is noise in state signals and how to deal with it. In addition, for future research focusing on the robustness of incremental control methods, it is expected to see how the incremental control method performs under more realistic turbulence models, such as the Dryden model.