Dynamic modeling and faster finite-time attitude stabilization of receiver aircraft for aerial refueling

Abstract

This paper studies the faster attitude stabilization problem of receiver aircraft during the refueling phase. Dynamic models of receiver aircraft subject to practical uncertainties and disturbances are established, where the uncertain inertia, the wind disturbance, the time-varying inertia, and the shift of the center of mass are explicitly considered. A faster finite-time stable system is preliminarily proposed and proved rigorously in comparison with the conventional fast finite-time stable system. Then, a nonsingular terminal sliding mode surface is synthesized to ensure a fast convergence rate. An adaptive faster finite-time controller is finally developed with the aircraft attitude stabilized into a small neighborhood of the origin with a faster convergence rate. Simulations are conducted to validate that the proposed solution can accomplish attitude stabilization maneuvers in a faster way.

1 Introduction

Aerial refueling is a mission to transfer fuel from the tanker aircraft to the receiver aircraft in flight [1]. It increases the endurance range and weight of receiver aircraft payload [2]. Note that attitude stabilization control of the receiver aircraft is the main objective in the refueling phase. At that time, the probe has been engaged with the drogue, as shown in Fig. 1. If the receiver aircraft’s attitude is not stabilized, the drogue may be separated from the probe. Then, the refueling task will fail. Moreover, the towing drogue may cause a collision to the receiver’s windshield. Therefore, the attitude stabilization problem [3] should be investigated for the receiver aircraft during the aerial refueling mission.

In the past three decades, many control approaches have been developed for the receiver aircraft such as the particle swarm-optimized PID control [4], the linear quadratic regulator (LQR) control [5], the state feedback control [6], ${\mathcal {L}}_1$ adaptive control [7, 8], the nonlinear dynamic inversion control [9], and the ADRC-based backstepping control [10, 11]. Those results were presented based on the assumption that the mass of the receiver aircraft is constant. However, in practice, when fuel continuously enters the receiver aircraft, its time-varying mass inevitably changes the inertia and the center of mass (CM). If these changes are not appropriately handled by the controller, it will deteriorate the attitude stabilization performance. Moreover, the external disturbance induced by the turbulence and the trailing vortices generated by a tanker aircraft will further affect attitude stabilization performance. Hence, the time-varying mass and the external disturbance should be explicitly addressed in the receiver aircraft’s attitude controller design.

To achieve high control performance of the receiver aircraft’s attitude system with time-varying mass, a gain-scheduling linear controller was designed to follow the desired receiver trajectory [12]. In addition, a LQR controller and an ADRC controller were presented in [13] and [14], respectively, to achieve attitude control in the refueling phase. Although the time-varying mass was explicitly considered, the linear controllers in [12] and [13] were developed using the linearized aircraft model. They lacked anti-disturbance capability and were less robust in the complex flow field during refueling phase. Moreover, the ADRC law in [14] was developed by assuming that the disturbance was constant or slowly time varying [15]. This assumption was merely satisfied in the aerial refueling engineering process.

Fig. 1

Refueling phase in the aerial refueling mission

Note that almost all of the previous control schemes lead to infinite settling time. This is not an optimal solution to accomplish the refueling task. Compared with asymptotic stabilization, the finite-time control strategy could provide faster convergence rate [16,17,18,19,20]. Hence, it could become a proper candidate control solution for attitude stabilization problem. However, the basic finite-time stable system in [21, 22] is slower than exponential convergence if the initial state is far from the origin. This drawback was modified by incorporating a linear term in [23, 24] to accelerate the convergence rate. More recently, a faster finite-time stable system was developed in [25, 26] by adding a power integrator. However, rigorous proof of the faster convergence rate was not provided. Only a numerical example was illustrated in [25].

With the above drawbacks in mind, this paper investigates the faster attitude stabilization control problem of the receiver aircraft during the refueling phase even in the presence of time-varying inertia, internal uncertainty, complex wind disturbances, and the shift of CM. The main contributions of this paper are listed as

(1) The mathematical models of the receiver aircraft’s attitude system, its time-varying and uncertain inertia, and disturbance in the aerial refueling task are established. Compared with [12, 13], the uncertain inertia and the shift of CM are explicitly considered. Moreover, the derived model is user-friendly for controller design.

(2) An adaptive faster finite-time attitude controller is designed to guarantee the practical finite-time stability of the attitude even in the presence of external disturbance and internal uncertainties. To the best knowledge of the authors, this may be the first result to achieve such good control performance for aerial refueling task.

(3) Compared with [22, 23], faster convergence rate is achieved by the proposed controller with rigorous proof provided, which is the extension of the result in [25, 26].

The rest of this paper is organized as follows. In Sect. 2, the nonlinear model with time-varying inertia of the aircraft attitude stabilization system and problem formulation are given. Section 3 presents the main results to achieve faster finite-time attitude stabilization control. Numerical simulations are shown in Sect. 4, and this paper is concluded in Sect. 5.

2 System modeling and problem formulation

Throughout this paper, the real and the positive real number sets are denoted as ${\mathbb {R}}$ and ${\mathbb {R}}_+$, respectively. ${\mathbb {N}}_0$ is the set of all nonnegative integers. The set of $m \in {\mathbb {N}}_0$ by $n \in {\mathbb {N}}_0$ real matrices is denoted as ${\mathbb {R}}^{m\times n}$. The inverse of the matrix ${\varvec{A}} \in {\mathbb {R}}^{m \times m}$ is defined as ${\varvec{A}}^{-1}$. ${\varvec{0}}$ is the zero matrix or vector with appropriate dimension. The skew-symmetric matrix of any vector ${\varvec{a}}=[a_1\ a_2 \ a_3]^{\mathrm{T}}\in {\mathbb {R}}^3$ is defined as ${\varvec{a}}^{\times }=\left[ {\begin{array}{*{20}{c}} {0}&{}\quad { - {a_3}}&{}\quad {{a_2}}\\ {{a_3}}&{}\quad {0}&{}\quad { - {a_1}}\\ { - {a_2}}&{}\quad {{a_1}}&{}\quad {0} \end{array}} \right]$. $\text {sig}({\varvec{a}})^{\gamma } = [ \text {sgn}(a_1)|a_1|^{\gamma } \quad \text {sgn}(a_2)|a_2|^{\gamma }\quad \text {sgn}(a_3)|a_3|^{\gamma } ]^{\mathrm{T}}$ is defined for any $\gamma \in {\mathbb {R}}_+$ with $\text {sgn}(\cdot )$ denoting the sign function. For any vector ${\varvec{v}}=[v_1\ v_2 \ \cdots \ v_n]^{\mathrm{T}}\in {\mathbb {R}}^n$, the diagonal matrix $\text {diag}({\varvec{v}}) = \text {diag}(v_i)=\left[ {\begin{array}{*{20}{c}} {{v_1}}&{}\quad 0&{}\quad \cdots &{}\quad 0\\ 0&{}\quad {{v_2}}&{}\quad \cdots &{}\quad 0\\ \vdots &{}\quad \vdots &{}\quad \ddots &{}\quad \vdots \\ 0&{}\quad 0&{}\quad \cdots &{}\quad {{v_n}} \end{array}} \right]$, and $\text {sgn}({\varvec{v}})= \text {diag}([\text {sgn}(v_1)\;\text {sgn}(v_2)\;\cdots \text {sgn}(v_n)]^\text {T})$ are defined. $||\cdot ||$ is the Euclidean norm of vectors in ${\mathbb {R}}^n$ or matrices in ${\mathbb {R}}^{n\times n}$ .

2.1 Modeling of the attitude control system of receiver aircraft

Let ${\varvec{\sigma }}=[\phi \quad \theta \quad \psi ]^{\mathrm{T}}\in {\mathbb {R}}^3$ denote the attitude Euler angles of the receiver aircraft with respect to the inertial frame ${\mathcal {F}}_i$. ${\varvec{\omega }}_b=[p \quad q \quad r]^{\mathrm{T}}\in {\mathbb {R}}^3$ described in the receiver aircraft’s body-fixed frame ${\mathcal {F}}_b$ represents the absolute angular velocity of the aircraft with respect to ${\mathcal {F}}_i$. Then, the attitude control system of the receiver aircraft can be modeled as [27,28,29]

$$\begin{aligned} {\varvec{\omega }}_b&={\varvec{R}}({\varvec{\sigma }})\dot{{\varvec{\sigma }}} , \end{aligned}$$

(1)

$$\begin{aligned} {\varvec{J}}\dot{{\varvec{\omega }}_b}&=-{\varvec{\omega }}_b^{\times }{\varvec{J\omega }}_b-\dot{{\varvec{J}}}{\varvec{\omega }}_b+{\varvec{u}}+{\varvec{d}} , \end{aligned}$$

(2)

where ${\varvec{J}} \in {\mathbb {R}}^{3 \times 3}$ is the total inertia matrix of the receiver aircraft, ${\varvec{u}} \in {\mathbb {R}}^3$ is total control torque generated by actuators, and ${\varvec{d}} \in {\mathbb {R}}^3$ is the total external disturbance. ${\varvec{R}}({\varvec{\sigma }}) \in {\mathbb {R}}^{3\times 3}$ is defined as

$$\begin{aligned} {\varvec{R}}({\varvec{\sigma }}) = \left[ {\begin{array}{*{20}{c}} 1&{}\quad 0&{}\quad { - \sin \theta }\\ 0&{}\quad {\cos \phi }&{}\quad {\cos \theta \sin \phi }\\ 0&{}\quad { - \sin \phi }&{}\quad {\cos \theta \cos \phi } \end{array}} \right] . \end{aligned}$$

(3)

Remark 1

In this paper, the receiver aircraft’s attitude is represented by Euler angles but not the unit quaternion [30]. Although singularity may occur in the case of $\theta = \frac{\pi }{2}$, the attitude of the receiver aircraft during the refueling phase will not encounter this case. Hence, applying Euler angles to describe the receiver aircraft’s attitude is reasonable and more intuitive.

In practice, ${\varvec{J}}$ is unknown due to limited modeling technique. Moreover, fuel transferred from the tanker aircraft will lead ${\varvec{J}}$ to be significantly time varying. In this work, let ${\varvec{J}}_0 \in {\mathbb {R}}^{3\times 3}$, ${\varvec{J}}_{\Delta }\in {\mathbb {R}}^{3\times 3}$, and ${\varvec{J}}_{\text {var}} \in {\mathbb {R}}^{3\times 3}$ denote the nominal constant part, the unknown constant part, and the time-varying uncertain part of the aircraft’s inertia, respectively. Then, it follows that:

$$\begin{aligned} {\varvec{J}}&={\varvec{J}}_0+{\varvec{J}}_{\Delta }+{\varvec{J}}_{\text {var}}, \end{aligned}$$

(4)

$$\begin{aligned} \dot{{\varvec{J}}}&=\dot{{\varvec{J}}}_{\text {var}}. \end{aligned}$$

(5)

The probe-and-drogue aerial refueling mission is considered in this work. Hence, the external disturbance acting on the receiver aircraft mainly consists of the wind disturbance torque ${\varvec{d}}_{\text {wind}} \in {\mathbb {R}}^3$ [1, 31] and the additional aerodynamic torque ${\varvec{d}}_{\text {aero}} \in {\mathbb {R}}^3$. Then, the external disturbance ${\varvec{d}}$ in (2) can be modeled as

$$\begin{aligned} {\varvec{d}}={\varvec{d}}_{\text {wind}}+{\varvec{d}}_{\text {aero}}. \end{aligned}$$

(6)

Assumption 1

The unknown inertia ${\varvec{J}}_{\Delta }$ is bounded. There is an unknown scalar $\kappa _1\in {\mathbb {R}}_+$ such that $||{\varvec{J}}_{\Delta }|| \le {\kappa _1}$.

2.2 Modeling of the time-varying uncertain inertia

The time-varying uncertain inertia ${\varvec{J}}_{\text {var}}$ is mathematically modeled in this subsection. The explicit modeling method is to theoretically calculate ${\varvec{J}}_{\text {var}}$ according to the definition of inertia using the mass distribution and the position of fuel tanks. However, this method is complicated and time-consuming. It is not applicable for real-time attitude stabilization of the receiver aircraft. In contrast, a novel modeling method is presented in this paper to model ${\varvec{J}}_{\text {var}}$ via the varying fuel mass [32]. It is achieved by modeling ${\varvec{J}}_{\text {var}}$ as

$$\begin{aligned} {\varvec{J}}_{\text {var}}=\text {diag}({\varvec{\rho }})m_{\text {var}}, \end{aligned}$$

(7)

where ${\varvec{\rho }}=[\rho _1 \quad \rho _2 \quad \rho _3]^{\mathrm{T}}\in {\mathbb {R}}^3$ is the correlation coefficient vector between the adding fuel $m_{\text {var}}\in {\mathbb {R}}_+$ and the time-varying inertia.

Note that the time-varying mass mainly changes the principal moment of inertia in three axes, i.e., $I_x\in {\mathbb {R}}_+$, $I_y\in {\mathbb {R}}_+$, and $I_z \in {\mathbb {R}}_+$. The change of the product of inertia, i.e., $I_{xz} \in {\mathbb {R}}$, $I_{yz} \in {\mathbb {R}}$, and $I_{xy} \in {\mathbb {R}}$ can be omitted. Figure 2 illustrates a typical curve between $I_x$ and $m_{\text {var}}$. With the increase in the fuel transferred from the tanker aircraft, $I_x$ rises almost linearly with different speed with respect to the $m_{\text {var}}$.

Fig. 2

Typical relationship of $I_x$ versus $m_{\text {var}}$

Remark 2

It should be stressed that the fuel transferred from the tanker aircraft is sophistically distributed by the fuel system to guarantee that the moment of inertia along x, y, and z axes and CM vary smoothly. There are multiple fuel tanks distributed in the receiver aircraft’s fuselage and the left/right wings. Fuel is added to those fuel tanks with the prescribed sequence to keep balance. Therefore, the correlation coefficient curve in Fig. 2 is a broken line with different slopes.

Assumption 2

Although ${\varvec{J}}_{\text {var}}$ is uncertain, its varying rate is bounded in practice. There is a scalar $\kappa _2\in {\mathbb {R}}_+$ satisfying:

$$\begin{aligned} ||\dot{{\varvec{J}}}_{\text {var}}||=||\frac{d{\varvec{J}}_{\text {var}}}{\mathrm{d}m_{\text {var}}} \frac{\mathrm{d}m_{\text {var}}}{\mathrm{d}t}|| \le \sqrt{3}\rho _{\text {max}}{\dot{m}}_{\text {var}}=\kappa _2, \end{aligned}$$

(8)

where $\rho _{\text {max}}=\text {max}\left\{ \rho _1 , \rho _2 , \rho _3\right\} >0$.

2.3 Modeling of wind disturbance torque

The wind disturbance torque ${\varvec{d}}_{\text {wind}}$ is mainly induced by the atmospheric turbulence and the tanker aircraft’s vortex effect, which are to be modeled in this part.

2.3.1 Atmospheric turbulence

It is known that the turbulence represents the stochastic part of the wind. As recommended in the MIL-F-8785C standard, the Dryden model can be used to describe the turbulence effect. The translational and the rotational power spectral density (PSD) functions $\varPhi ({\omega })$ in terms of the angular frequency $\omega \in {\mathbb {R}}$ can be computed by

$$\begin{aligned} \varPhi _u(\omega )&=\frac{2\sigma _u^2L_u}{\pi V_T}\frac{1}{1+(L_u\frac{\omega }{V_T})^2}, \end{aligned}$$

(9)

$$\begin{aligned} \varPhi _v(\omega )&=\frac{2\sigma _v^2L_v}{\pi V_T}\frac{1+12(L_v\frac{\omega }{V_T})^2}{1+4(L_v\frac{\omega }{V_T})^2},\end{aligned}$$

(10)

$$\begin{aligned} \varPhi _w(\omega )&=\frac{2\sigma _w^2L_w}{\pi V_T}\frac{1+12(L_w\frac{\omega }{V_T})^2}{1+4(L_w\frac{\omega }{V_T})^2},\end{aligned}$$

(11)

$$\begin{aligned} \varPhi _p(\omega )&=\frac{\sigma _w^2}{L_wV_T}\frac{0.4(\frac{\pi L_w}{2b})^\frac{1}{3}}{1+(\frac{4b\omega }{\pi V_T})^2},\end{aligned}$$

(12)

$$\begin{aligned} \varPhi _q(\omega )&=\left( \frac{\omega }{V_T}\right) ^2\frac{\phi _v(\omega )}{1+(\frac{3b\omega }{\pi V_T})^2}\end{aligned}$$

(13)

$$\begin{aligned} \varPhi _r(\omega )&=\left( \frac{\omega }{V_T}\right) ^2\frac{\phi _w(\omega )}{1+(\frac{4b\omega }{\pi V_T})^2}, \end{aligned}$$

(14)

where $V_T \in {\mathbb {R}}_+$ is the airspeed and $b \in {\mathbb {R}}_+$ is the span length. $\sigma _u \in {\mathbb {R}}_+$, $\sigma _v \in {\mathbb {R}}_+$, and $\sigma _w \in {\mathbb {R}}_+$ are the root-mean-square turbulence intensities. $L_u \in {\mathbb {R}}_+$, $L_v \in {\mathbb {R}}_+$, and $L_w \in {\mathbb {R}}_+$ are the scale lengths.

Applying the spectrum decomposition to the PSD functions (9)–(14), the corresponding shaping filters of PSD functions are determined by

$$\begin{aligned} G_u(s)&=\sigma _u\sqrt{\frac{2L_u}{\pi V}}\frac{1}{{\frac{L_u}{V}s+1}}, \end{aligned}$$

(15)

$$\begin{aligned} G_v(s)&=\sigma _v\sqrt{\frac{2L_v}{\pi V}}\frac{1}{\frac{2L_v}{\sqrt{3}V}s+1},\end{aligned}$$

(16)

$$\begin{aligned} G_w(s)&=\sigma _w\sqrt{\frac{2L_w}{\pi V}}\frac{1}{\frac{2L_w}{\sqrt{3}V}s+1},\end{aligned}$$

(17)

$$\begin{aligned} G_p(s)&=\frac{K_p}{T_ps+1},\end{aligned}$$

(18)

$$\begin{aligned} G_q(s)&=\frac{K_{q}s G_w(s)}{T_qs+1},\end{aligned}$$

(19)

$$\begin{aligned} G_r(s)&=\frac{1}{T_rs+1}, \end{aligned}$$

(20)

where $K_p=\sigma _w \sqrt{\frac{0.4}{L_w V}(\frac{\pi L_w}{2b})^\frac{1}{3}}$, $T_p=T_q=\frac{4b}{\pi V}$, $K_q=K_r=\frac{1}{V}$, and $T_r=\frac{3b}{\pi V}$.

Based on the shaping filters (15)–(20), the induced translational velocities $u_{w1} \in {\mathbb {R}}$, $v_{w1} \in {\mathbb {R}}$, and $w_{w1} \in {\mathbb {R}}$ as well as the induced angular velocities $p_{w1} \in {\mathbb {R}}$, $q_{w1} \in {\mathbb {R}}$, and $r_{w1} \in {\mathbb {R}}$ introduced by the wind turbulence are derived in Fig. 3, where white noise is the input of those transfer functions.

Fig. 3

Structure of the induced velocities generator

2.3.2 Vortex effect

During the aerial refueling mission, the receiver aircraft will encounter the tanker aircraft’s trailing vortex. This will affect the receiver aircraft’s attitude dynamics. Figure 4 illustrates the static velocity contours of a typical trailing vortex effect.

Fig. 4

Static velocity contours in CFD

The tanker aircraft produces two straight and semi-infinite trailing vortex filaments in rear regions. This vortex-induced wind field acting on the receiver aircraft is highly non-uniform. It leads to changes in the aerodynamic forces and the moments. However, direct computation of those forces and moments is not effective, since the standard aerodynamics derived using the airspeed, the angles of attack, and the sideslip are not suitable in the non-uniform situation. To solve this problem, a feasible approach is to compute the induced wind velocities and the wind gradients as a function of the relative distance and the orientation by invoking a modified horseshoe vortex model [33]. Then, the averaging technique can be applied to transform the non-uniform wind components and gradients into the uniform ones. Consequently, the equivalent wind components $u_{w2}\in {\mathbb {R}}$, $v_{w2}\in {\mathbb {R}}$, and $w_{w2}\in {\mathbb {R}}$ as well as the equivalent gradients $p_{w2}\in {\mathbb {R}}$, $q_{w2} \in {\mathbb {R}}$, and $r_{w2} \in {\mathbb {R}}$ can be introduced into equations of motion and the aerodynamics. The detailed modeling procedures can be referred to [31, 34].

2.3.3 Calculation of wind disturbance torque

Based on the preceding results, the wind disturbance torque ${\varvec{d}}_{\text {wind}}$ can be finally modelled as

$$\begin{aligned} {\varvec{d}}_{\text {wind}}=QS_w\text {diag}({\varvec{g}}) \left[ {\begin{array}{*{20}{c}} {{c_{l*}(\beta ,Ma,p,r)} - {c_{l}}}\\ {{c_{m*}(\alpha ,Ma,q)} - {c_{m}}}\\ {{c_{n*}(\beta ,Ma,p,r)} - {c_{n}}} \end{array}} \right] , \end{aligned}$$

(21)

where ${\varvec{g}}=[b\quad c\quad b]^{\mathrm{T}}$, $Q \in \mathbb {R_+}$ is the dynamic pressure, $S_w \in {\mathbb {R}}_+$ is reference area. $b \in {\mathbb {R}}_+$ and $c \in {\mathbb {R}}_+$ are the span and the chord length, respectively. $c_l \in {\mathbb {R}}$, $c_m\in {\mathbb {R}}$, and $c_n\in {\mathbb {R}}_+$ are the nominal moment coefficients for the roll, the pitch, and the yaw axis, respectively. $c_{l*} \in {\mathbb {R}}$, $c_{m*}\in {\mathbb {R}}$, and $c_{n*}\in {\mathbb {R}}$ are moment coefficients induced by the wind disturbance for the roll, the pitch, and the yaw axis, respectively. $\alpha \in {\mathbb {R}}$ is the angle of attack, $\beta \in {\mathbb {R}}$ is the sideslip angle, and $Ma \in {\mathbb {R}}_+$ is the Mach number; moreover, they are calculated as

$$\begin{aligned} \alpha&=\arctan \left( \frac{w_*}{u_*}\right) , \end{aligned}$$

(22)

$$\begin{aligned} \beta&=\arcsin \left( \frac{v_*}{\sqrt{u_*^2+v_*^2+w_*^2}}\right) , \end{aligned}$$

(23)

$$\begin{aligned} Ma&=\frac{1}{a}\sqrt{u_*^2+v_*^2+w_*^2}, \end{aligned}$$

(24)

with $u_*=u+u_{w1}+u_{w2}$, $v_*=v+v_{w1}+v_{w2}$, and $w_*=w+w_{w1}+w_{w2}$. $u \in {\mathbb {R}}$, $v \in {\mathbb {R}}$, and $w \in {\mathbb {R}}$ represent the translational velocities in ${\mathcal {F}}_b$. a denotes the speed of sound.

Assumption 3

The wind disturbance torque ${\varvec{d}}_{\text {wind}}$ is bounded by $||{\varvec{d}}_{\text {wind}}||\le \kappa _3$, where $\kappa _3 \in {\mathbb {R}}_+$ is a scalar.

2.4 Modeling of additional aerodynamic torque

Let ${\varvec{F}}_{\text {aero}} \in {\mathbb {R}}^3$ be the aerodynamic force acting on the center of pressure (CP), which is usually defined at the position of 0.25 chord. When CM does not coincide with CP, a rotation torque called the additional aerodynamic torque will be generated by ${\varvec{F}}_{\text {aero}}$. Let ${\varvec{d}}_{\text {aero}}\in {\mathbb {R}}^3$ represent that aerodynamic torque, then it can be mathematically modeled as

$$\begin{aligned} {\varvec{d}}_{\text {aero}}=-({\varvec{r}}_{\text {var}})^\times {\varvec{F}}_{\text {aero}}, \end{aligned}$$

(25)

where ${\varvec{r}}_{\text {var}}={\varvec{x}}_{\text {CM}}-{\varvec{x}}_{\text {CP}}$ is the position error between CM and CP.

Assumption 4

During the whole attitude stabilization process in the aerial refueling, ${\varvec{d}}_{\text {aero}}$ is bounded. An unknown constant $\kappa _4 \in {\mathbb {R}}_+$ exists such that $||{\varvec{d}}_{\text {aero}} || \le \kappa _4$.

Remark 3

In the refueling process, the airspeed and the attitude of the receiver aircraft almost remain stable. In this case, the changes in the angle of attack, the sideslip, and the Mach number are relatively small even in the presence of the wind disturbance. Hence, the aerodynamic force ${\varvec{F}}_{\text {aero}}$ and the moments are bounded. Moreover, CM only varies in a small range under the prescribed refueling sequence. Assumption 3 and Assumption 4 are therefore reasonable [10, 11, 30].

2.5 Problem formulation

Suppose that there are attitude sensors and gyros fixed in the receiver aircraft to measure its attitude ${\varvec{\sigma }}$ and its angular velocity ${\varvec{\omega }}_b$. In this paper, we aim at the faster finite-time attitude stabilization of receiver aircraft in the presence of external disturbance ${\varvec{d}}$, the uncertain inertia ${\varvec{J}}_{\Delta }$, and the time-varying uncertain inertia ${\varvec{J}}_{\text {var}}$. The control objective is to design a feedback control law to ensure that the closed-loop system (1) - (2) is finite-time stable, that is, $\lim _{x \rightarrow t_f} {\varvec{\sigma }}={\varvec{0}}$ and $\lim _{x \rightarrow t_f} {\varvec{\omega }}_b={\varvec{0}}$. Moreover, the settling time $t_f$ should be smaller than the conventional fast finite-time stable system [24].

2.6 Useful lemmas

Lemma 1

[24] Consider a system described as

$$\begin{aligned} {\dot{z}}=-m_1z-m_2z^{q_2}, \end{aligned}$$

(26)

where z is the system state, $m_1 \in {\mathbb {R}}_+$, $m_2 \in {\mathbb {R}}_+$, $0<q_2<1$, and $z_0=z(0)$ is the initial value. Then, system (26) is globally fast finite-time stable with the settling time $T\in {\mathbb {R}}_+$ bounded by

$$\begin{aligned} T \le \frac{1}{m_1(1-q_2)}\text {ln}\left( \frac{m_1z_0^{1-q_2}+m_2}{m_2}\right) . \end{aligned}$$

(27)

Lemma 2

([30, 35]) For ${\varvec{z}}=[z_1 \quad z_2 \quad \cdots \quad z_N]^{\text {T}} \in {\mathbb {R}}^N$ and constants $q_1 \in {\mathbb {R}}_+$, $q_2 \in {\mathbb {R}}_+$ satisfying $q_1 > 1$ and $0<q_2<1$, it follows that:

$$\begin{aligned} N^{q_1-1} \sum _{i=1}^{N} |z_i|^{q_1}&\ge \left( \sum _{i=1}^{N} |z_i| \right) ^{q_1},\end{aligned}$$

(28)

$$\begin{aligned} \sum _{i=1}^{N} |z_i|^{q_2}&\le N^{1-q_2} \left( \sum _{i=1}^{N} |z_i| \right) ^{q_2},\end{aligned}$$

(29)

$$\begin{aligned} \sum _{i=1}^{N} |z_i|^{q_2}&\ge \left( \sum _{i=1}^{N} |z_i| \right) ^{q_2}. \end{aligned}$$

(30)

3 Faster finite-time controller design

In this section, a faster finite-time attitude stabilization control approach is presented based on the faster finite-time stable system and the nonsingular terminal sliding mode control technique.

To facilitate the controller design, the receiver aircraft’s attitude control system (1)–(2) with the time-varying inertia (7) and the external disturbance (6) are transformed into

$$\begin{aligned} \left\{ \begin{array}{ll} \begin{aligned} \dot{{\varvec{x}}}_1&{}={\varvec{x}}_2,\\ \dot{{\varvec{x}}}_2&{}=-{\varvec{M}}^{-1}({\varvec{x}}_1){\varvec{C}}({\varvec{x}}_1,{\varvec{x}}_2){\varvec{x}}_2 + {\varvec{\tau }} + \bar{{\varvec{d}}}, \end{aligned} \end{array}\right. \end{aligned}$$

(31)

where ${\varvec{x}}_1={\varvec{\sigma }}=[x_{11} \quad x_{12} \quad x_{13}]^{\mathrm{T}}$, ${\varvec{M}}({\varvec{x}}_1)={\varvec{J}}_0{\varvec{R}}({\varvec{x}}_1)$, ${\varvec{x}}_2=\dot{{\varvec{\sigma }}}=[x_{21} \quad x_{22} \quad x_{23}]^{\mathrm{T}}$, ${\varvec{C}}({\varvec{x}}_1,{\varvec{x}}_2) ={\varvec{J}}_0\dot{{\varvec{R}}}({\varvec{x}}_1)+{\varvec{\omega }}^{\times }_b{\varvec{J}}_{0}{\varvec{R}}({\varvec{x}}_1)$, ${\varvec{\tau }}={\varvec{M}}^{-1}({\varvec{x}}_1){\varvec{u}}$, and

$$\begin{aligned} \begin{aligned} \bar{{\varvec{d}}}=&{\varvec{M}}^{-1}({\varvec{x}}_1)( {\varvec{d}}_{\text {aero}}+{\varvec{d}}_{\text {wind}} - {\varvec{J}}_{\Delta }\dot{{\varvec{\omega }}}_b \\&- {\varvec{J}}_{\text {var}}\dot{{\varvec{\omega }}}_b - {\varvec{\omega }}^{\times }_b{\varvec{J}}_{\Delta }{\varvec{\omega }}_b - {\varvec{\omega }}^{\times }_b{\varvec{J}}_{\text {var}}{\varvec{\omega }}_b - \dot{{\varvec{J}}}_{\text {var}}{\varvec{\omega }}_b), \end{aligned} \end{aligned}$$

(32)

which denotes the lumped disturbance.

Remark 4

Because the attitude and angular velocity of the receiver aircraft can be measured by the fixed sensors, ${\varvec{x}}_1$ , ${\varvec{x}}_2$, ${\varvec{M}}({\varvec{x}}_1)$, and ${\varvec{C}}({\varvec{x}}_1,{\varvec{x}}_2)$ are measurable. They can be directly applied in the attitude stabilization controller design.

3.1 Development of the faster finite-time stable system

In this subsection, a faster finite-time stable system is presented as the backstone of the subsequent design of controllers. It is developed as

$$\begin{aligned} {\dot{z}}=-m_1z^{q_1}-m_2z^{q_2}, \end{aligned}$$

(33)

where $z\in {\mathbb {R}}_+$ is the system state, $m_1 \in {\mathbb {R}}_+$, $m_2$ $\in {\mathbb {R}}_+$, $q_1 > 1$ and $0<q_2<1$. $z_0=z(0)$ is the initial state.

Theorem 1

System (33) is finite-time stable for any $z_0$ with the settling time $T_1 \in {\mathbb {R}}_+$, i.e., $z(t)\equiv 0$, $\forall \ t \ge T_1$, where

$$\begin{aligned} T_1=\frac{1}{m_2(1-q_2)}+\frac{1-z_0^{1-q_1}}{m_1(q_1-1)}. \end{aligned}$$

(34)

Moreover, the convergence rate is faster than the finite-time system presented in [24].

Proof

Follow the method in [25] and [26], solving (33) can directly prove $z(t)\equiv 0$ for $\forall \ t \ge T_1$.

Introducing $y=z^{1-q_2}$ and applying (33) results in

$$\begin{aligned} {\dot{y}}=(1-q_2)z^{-q_2}{\dot{z}}=-(1-q_2)(m_2+m_1y^{\frac{q_1-q_2}{1-q_2}}) . \end{aligned}$$

(35)

Defining $\eta =\frac{q_1-q_2}{1-q_2}$, one can rewrite (35) as

$$\begin{aligned} -\frac{\mathrm{d}y}{(1-q_2)(m_2+m_1y^{\eta })}=\mathrm{d}t. \end{aligned}$$

(36)

Then, the settling time of the system (33) satisfies $T_1=\int _0^{y_0} \frac{\mathrm{d}y}{(1-q_2)(m_2+m_1y^{\eta })}$, where $y_0=y(0)=z(0)^{1-q_2}$. For the finite-time system ${\dot{z}}=-m_1z-m_2z^{q_2}$ given in [24], its settling time is denoted as $T_2$ and it follows that $T_2=\int _0^{y_0} \frac{\mathrm{d}y}{(1-q_2)(m_2+m_1y)}$.

Define ${\widetilde{T}}=T_2-T_1$, then it leads to

$$\begin{aligned} {\widetilde{T}}=&\int _0^{y_0} \frac{1}{(1-q_2)}\Big (\frac{1}{m_2+m_1y}-\frac{1}{m_2+m_1y^{\eta }}\Big )\mathrm{d}y \nonumber \\ =&\int _1^{y_0} \frac{m_1(y^{\eta }-y)\mathrm{d}y}{(1-q_2)(m_2+m_1y)(m_2+m_1y^{\eta })}\nonumber \\&-\int _0^{1} \frac{m_1(y-y^{\eta })\mathrm{d}y}{(1-q_2)(m_2+m_1y)(m_2+m_1y^{\eta })}. \end{aligned}$$

(37)

Because $\frac{1}{(m_2+m_1y)(m_2+m_1y^{\eta })}$ is a decreasing function and $\eta >1$, invoking the second mean value theorem for integrals, we can prove that:

$$\begin{aligned} \int _0^{1} \frac{(y-y^{\eta })\mathrm{d}y}{(m_2+m_1y)(m_2+m_1y^{\eta })}&=\frac{1}{m_2^2}\int _0^{\xi _1} (y-y^{\eta })\mathrm{d}y\nonumber \\&<\frac{\frac{1}{2}-\frac{1}{\eta +1}}{m_2^2}, \end{aligned}$$

(38)

where $\xi _1$ is a constant satisfying $0< \xi _1 < 1$. On the other hand, it yields

$$\begin{aligned}&\int _1^{y_0} \frac{(y^{\eta }-y)\mathrm{d}y}{(m_2+m_1y)(m_2+m_1y^{\eta })}=\frac{\int _1^{\xi _2} (y^{\eta }-y)\mathrm{d}y}{(m_1+m_2)^2}\nonumber \\&=\frac{(\frac{1}{2}-\frac{1}{\eta +1})+(\frac{\xi _2^{\eta +1}}{\eta +1}-\frac{\xi _2^2}{2})}{(m_1+m_2)^2} . \end{aligned}$$

(39)

where $\xi _2$ is a scalar such that $1< \xi _2 < y_0$.

Denote $l(\xi _2)= \frac{(\frac{1}{2}-\frac{1}{\eta +1})+(\frac{\xi _2^{\eta +1}}{\eta +1}-\frac{\xi _2^2}{2})}{(m_1+m_2)^2}-\frac{\frac{1}{2}-\frac{1}{\eta +1}}{m_2^2}$, $l(\xi _2)>0$ can be guaranteed by appropriately choosing $m_1$, $m_2$, and $\eta$ with $y_0=z_0^{1-q_2}$. Then, it is proved from ${\widetilde{T}} >0$ that the convergence rate of system (33) is faster than the finite-time system given in [24]. $\square$

Theorem 2

Consider the system described by

$$\begin{aligned} {\dot{z}}=-m_1z^{q_1}-m_2z^{q_2}+\epsilon , \end{aligned}$$

(40)

where $z\in {\mathbb {R}}_+$, $m_1\in {\mathbb {R}}_+$, $m_2\in {\mathbb {R}}_+$, $q_1\in {\mathbb {R}}_+$, $q_2\in {\mathbb {R}}_+$ are defined as same as (33), and $0< \epsilon < \infty$; then, the trajectory of this system is practical finite-time stable. Moreover, the residual set of the solution of the system is given by

$$\begin{aligned} \left\{ \mathop {\lim }\limits _{t \rightarrow {T_R}}\Big |z(t) \le z_{\text {min}}\right\} , \end{aligned}$$

(41)

where $z_{\text {min}} = \text {min} \left\{ (\frac{\epsilon }{m_1(1-\mu )})^{\frac{1}{q_1}}, (\frac{\epsilon }{m_2(1-\mu )})^{\frac{1}{q_2}} \right\}$, $0<\mu <1$ is a constant. The settling time $T_R$ is bounded as

$$\begin{aligned} T_R \le \text {max} \Big \{ \frac{1}{m_2(1-q_2)}+\frac{z(0)^{1-q_1}}{m_1\mu (q_1-1)} \;, \nonumber \\ \;\frac{1}{m_2\mu (1-q_2)}+\frac{z(0)^{1-q_1}}{m_1(q_1-1)} \Big \}. \end{aligned}$$

(42)

Proof

For any scalar $\mu$ satisfying $0<\mu <1$, it leaves (40) as

$$\begin{aligned} {\dot{z}}=-m_1 \mu z^{q_1}-m_2z^{q_2}+\epsilon -m_1(1-\mu )z^{q_1} \end{aligned}$$

(43)

or

$$\begin{aligned} {\dot{z}}=-m_1z^{q_1}-m_2 \mu z^{q_2}+\epsilon -m_2(1-\mu )z^{q_2}. \end{aligned}$$

(44)

It follows from (43) that ${\dot{z}} \le -m_1 \mu z^{q_1}-m_2z^{q_2}$, if $z^{q_1} > \frac{\epsilon }{m_1(1-\mu )}$. Invoking the comparison principle, Theorem 1, and reference [3], we can prove that after finite-time the system state z will fall into the set:

$$\begin{aligned} \left\{ z(t)^{q_1} \le \frac{\epsilon }{m_1(1-\mu )}\right\} , \end{aligned}$$

(45)

and the reaching time is given by

$$\begin{aligned} T_{r1} \le \frac{1}{m_2(1-q_2)}+\frac{z(0)^{1-q_1}}{m_1\mu (q_1-1)}. \end{aligned}$$

(46)

Based on (44), it also yields ${\dot{z}} \le -m_1 z^{q_1}-m_2 \mu z^{q_2}$, if $z^{q_2} > \frac{\epsilon }{m_2(1-\mu )}$. Similarly, the system state z will fall into the following set:

$$\begin{aligned} \left\{ z(t)^{q_2} \le \frac{\epsilon }{m_2(1-\mu )}\right\} , \end{aligned}$$

(47)

with the reaching time:

$$\begin{aligned} T_{r2} \le \frac{1}{m_2 \mu (1-q_2)}+\frac{z(0)^{1-q_1}}{m_1(q_1-1)}. \end{aligned}$$

(48)

Hence, the conclusion in Theorem 2 is proved based on (45) and (47). $\square$

Remark 5

Theorem 1 is the extension of the results presented in [25, 26]. In comparison with [25, 26], the faster convergence rate achieved by system (33) is rigorously proved in this paper, if parameters are chosen to satisfy $l(\xi _2)>0$. There must exist an initial condition z(0) with $1<\xi _2<z(0)$ ensuring that the convergence time $T_1$ is shorter than $T_2$. However, Remark 1 in [25] only qualitatively demonstrated the faster convergence performance using a numerical example. It also proves that the convergence time of the system (33) is not always shorter than the finite-time system [25, 26].

Remark 6

The fixed-time stable system developed in [36] almost has the similar form of the proposed faster finite-time system (33). However, the convergence time achieved in this paper is shorter than [36]. This can be verified by $T_1 - T_{\text {fixd}}= - \frac{z_0^{1-q_1}}{m_1(q_1-1)} < 0$ due to $z_0=z(0)>0$. Thus, $T_1 < T_{\text {fixed}}$ is valid.

3.2 Development of the fast finite-time sliding mode surface

Based on Lemma 1, a nonsingular sliding mode surface ${\varvec{s}}=[s_1 \quad s_2 \quad s_3]^\text {T}$ is synthesized to overcome the singular problem [37]:

$$\begin{aligned} s_i=h(x_{1i})x_{1i}+\text {sgn}(x_{2i})|x_{2i}|^{\lambda }, i=1,2,3, \end{aligned}$$

(49)

where $h(x_{1i}) = (m_1|x_1|^{q_{s1}-\frac{1}{\lambda }}+m_2|x_1|^{q_{s2}-\frac{1}{\lambda }})^{\lambda }$ with $m_1 \in {\mathbb {R}}_+$, $m_2 \in {\mathbb {R}}_+$, $1< \lambda <2$, $q_{s1}=1$, and $0<q_{s2}<1$.

Theorem 3

If a controller can govern the system states in (31) to reach the sliding mode surface ${\varvec{s}}={\varvec{0}}$ and stay on it thereafter, then the system states can converge to ${\varvec{0}}$ after a finite time $T_{s} \in {\mathbb {R}}_+$, which is bounded as

$$\begin{aligned} \begin{aligned} T_{s}=\max _{i=1,2,3}\left\{ \frac{\ln \left( \frac{m_1x_{1i}(0)^{1-q_{s2}}+m_2}{m_2}\right) }{m_1(1-q_{s2})}\right\} . \end{aligned} \end{aligned}$$

(50)

Proof

Since $h(x_{1i})>0$ always holds, once the states reach the sliding mode surface ${\varvec{s}}=0$, it follows that:

$$\begin{aligned} h(x_{1i})x_{1i}=-\text {sgn}(x_{2i})|x_{2i}|^{\lambda }=-\text {sgn}({\dot{x}}_{1i})|{\dot{x}}_{1i}|^{\lambda }, \end{aligned}$$

(51)

and

$$\begin{aligned} {\dot{x}}_{1i}&=-h(x_{1i})^{\frac{1}{\lambda }}|x_{1i}|^{\frac{1}{\lambda }}\text {sgn}(x_{1i})\nonumber \\&=-(m_1|x_{1i}|^{1-\frac{1}{\lambda }}+m_2|x_{1i}|^{q_{s2}-\frac{1}{\lambda }})|x_{1i}|^{\frac{1}{\lambda }}\text {sgn}(x_{1i}) \nonumber \\&=-m_1|x_{1i}|\text {sgn}(x_{1i})-m_2|x_{1i}|^{q_{s2}}\text {sgn}(x_{1i}). \end{aligned}$$

(52)

Invoking Lemma 1, it can prove that the system states will converge to zero in a finite time $T_{s_i}, i=1,2,3$, i.e.,

$$\begin{aligned} T_{s_i} =\frac{1}{m_1(1-q_{s2})}\ln \frac{m_1x_{1i}(0)^{1-q_{s2}}+m_2}{m_2}. \end{aligned}$$

(53)

It leaves the convergence time as $T_s=\max \{T_{s_i}\}$. $\square$

3.3 Adaptive attitude stabilization controller design

According to Assumptions 1–4, one has:

$$\begin{aligned} \Vert {\bar{{\varvec{d}}}} \Vert =&|| {\varvec{M}}^{-1}({\varvec{x}}_1)( {\varvec{d}}_{\text {aero}}+{\varvec{d}}_{\text {wind}} - {\varvec{J}}_{\Delta }\dot{{\varvec{\omega }}}_b - {\varvec{J}}_{\text {var}}\dot{{\varvec{\omega }}}_b \nonumber \\&- {\varvec{\omega }}^{\times }_b{\varvec{J}}_{\Delta }{\varvec{\omega }}_b - {\varvec{\omega }}^{\times }_b{\varvec{J}}_{\text {var}}{\varvec{\omega }}_b - \dot{{\varvec{J}}}_{\text {var}}{\varvec{\omega }}_b ) || \nonumber \\ \le&|| {\varvec{M}}^{-1}({\varvec{x}}_1) || \big ( \kappa + \kappa _1 \left\| \dot{{\varvec{\omega }}}_b\right\| + \kappa _2 \left\| {\varvec{\omega }}_b\right\| \nonumber \\&+ \kappa _1\left\| {\varvec{\omega }}_b\right\| ^2 + \sqrt{3}\rho _{\text {max}}m_{\text {var}}(\left\| \dot{{\varvec{\omega }}}_b\right\| +\left\| {\varvec{\omega }}_b\right\| ^2 )\big )\nonumber \\ \le&\zeta _1+ \zeta _2 D_1+ \zeta _3 D_2, \end{aligned}$$

(54)

where $\kappa =\kappa _3+\kappa _4$, $\zeta _1 \in {\mathbb {R}}_+$ and $\zeta _2 \in {\mathbb {R}}_+$ are two unknown scalars, $\zeta _3 \in {\mathbb {R}}_+$ is a known bound related to $||{\varvec{M}}^{-1}({\varvec{x}}_1)|| \sqrt{3}\rho _{\text {max}}m_{\text {var}}$. $D_1 = \max \{||\dot{{\varvec{\omega }}}_b||,||{\varvec{\omega }} _b||,||{\varvec{\omega }} _b||^2\}$ and $D_2= \left\| {{\dot{{\varvec{\omega }}} _b}} \right\| + {{\left\| {{{\varvec{\omega }} _b}} \right\| }^2}$ are measurable.

Design the attitude stabilization controller as

$$\begin{aligned} {\varvec{\tau }}=&-\sum _{j=1}^{2}\gamma _j \text {sig}({\varvec{s}})^{q_j}+{{\varvec{M}}^{ - 1}}( {{{\varvec{x}}_1}}){\varvec{C}}( {{{\varvec{x}}_1},{{\varvec{x}}_2}}){{\varvec{x}}_2}\nonumber \\&- \frac{1}{\lambda } \left( {\varvec{H}}({\varvec{x}}_1)+\bar{{\varvec{H}}}({\varvec{x}}_1) \right) \text {sig}({\varvec{x}}_2)^{2-\lambda }+{\varvec{u}}_\text {adp}, \end{aligned}$$

(55)

where $\gamma _1 \in {\mathbb {R}}_{+}$, $\gamma _2 \in {\mathbb {R}}_{+}$, $q_1>1$, $0<q_2<1$, $1<\lambda <2$, ${\varvec{H}}({\varvec{x}}_1)=\text {diag}([h(x_{11}) \quad h(x_{12}) \quad h(x_{13})]^{\mathrm{T}})$, and $\bar{{\varvec{H}}}({\varvec{x}}_1) = \text {diag}([{\bar{h}}(x_{11}) \quad {\bar{h}}(x_{12}) \quad {\bar{h}}(x_{13})]^\text {T})$ with

$$\begin{aligned} {\bar{h}}(x_{1i})= \,&\lambda \big (\sum _{j=1}^{2}m_j|x_{1i}|^{q_{sj} - \frac{1}{\lambda }} \big )^{\lambda - 1}\big ( {m_1}( {{q_{s1}} - \frac{1}{\lambda }})|{x_{1i}} |^{ q_{s1}- \frac{1}{\lambda }} \nonumber \\&+ {m_2} ( {q_{s2}} - \frac{1}{\lambda } ) |{x_{1i}}|^{{{q_{s2}} - \frac{1}{\lambda }}} \big ), i=1,2,3. \end{aligned}$$

(56)

Moreover, ${\varvec{u}}_{\text {adp}}$ is determined as

$$\begin{aligned} {\varvec{u}}_{\text {adp}}=-\text {sgn}({\varvec{s}})( {\hat{\zeta }}_1+ {\hat{\zeta }}_2D_1 + \zeta _3D_2), \end{aligned}$$

(57)

and updated by

$$\begin{aligned} \dot{{\hat{\zeta }}}_1&= \frac{1}{\vartheta _0c_1}(\lambda \sqrt{3}^{2-\lambda } ||{\varvec{s}}|| ||{\varvec{x}}_2||^{\lambda -1} -g_1 {\hat{\zeta }}_1), \end{aligned}$$

(58)

$$\begin{aligned} \dot{{\hat{\zeta }}}_2&= \frac{1}{\vartheta _0c_2}(\lambda \sqrt{3}^{2-\lambda } ||{\varvec{s}}|| ||{\varvec{x}}_2||^{\lambda -1}D_1 -g_2 {\hat{\zeta }}_2 ), \end{aligned}$$

(59)

where $0<\vartheta _0<1$ is a constant satisfying $\vartheta _0^{(1+q_2)/2}+\vartheta _0-1=0$, $c_1= \frac{2g_1 (2 \varepsilon _1-1)}{2 \varepsilon _1}$, $c_2=\frac{2g_2 (2 \varepsilon _2-1)}{2 \varepsilon _2}$, $\varepsilon _1 > 0.5$, $\varepsilon _2 > 0.5$, $g_1 \in {\mathbb {R}}_+$, and $g_2 \in {\mathbb {R}}_+$.

Theorem 4

For the transformed receiver aircraft’s dynamics (31) with the lumped disturbance $\bar{{\varvec{d}}}$ induced by the wind turbulence, the vortex effect, the time-varying inertia, and the shift of CM, applying the adaptive controller (57) , then the closed-loop attitude system is practically finite-time stable. The convergence time is faster than the conventional fast finite-time stable system [24].

Proof

Applying (31) and (49), we can calculate $\dot{{\varvec{s}}}$ as

$$\begin{aligned} \dot{{\varvec{s}}}=\,&{\varvec{H}}({\varvec{x}}_1){\varvec{x}}_1+{\varvec{H}}({\varvec{x}}_1){\varvec{x}}_2+\lambda \text {diag} \left( |x_{2i}|^{\lambda -1} \right) \nonumber \\&\left( - {{\varvec{M}}^{ - 1}}\left( {{{\varvec{x}}_1}} \right) {\varvec{C}}\left( {{{\varvec{x}}_1},{{\varvec{x}}_2}} \right) {{\varvec{x}}_2} + {\varvec{\tau }} + \bar{{\varvec{d}}}\, \right) \nonumber \\ =\,&\text {diag} \Big ( \lambda \big (m_1|x_{1i}|^{q_{s1}-\frac{1}{\lambda }} + m_2|x_{1i}|^{q_{s2}-\frac{1}{\lambda }} \big )^{\lambda -1} \nonumber \\&\text {sgn}(x_{1i}) \big ( m_1(q_{s1}-\frac{1}{\lambda })|x_{1i}|^{q_{s1}-\frac{1+\lambda }{\lambda }} \nonumber \\&+ m_2(q_{s2}-\frac{1}{\lambda } ) |x_{1i}|^{q_{s2}-\frac{1+\lambda }{\lambda }} \big ) {\dot{x}}_{1i} \Big ){\varvec{x}}_1 \nonumber \\&+{\varvec{H}}({\varvec{x}}_1){\varvec{x}}_2+\lambda \text {diag} \left( |x_{2i}|^{\lambda -1} \right) \big ( {\varvec{\tau }}+\bar{{\varvec{d}}} \nonumber \\&- {\varvec{M}}^{-1}({\varvec{x}}_1){\varvec{C}}({\varvec{x}}_1,{\varvec{x}}_2){\varvec{x}}_2 \big ) \nonumber \\ =\,&\text {diag} \Big ( \lambda \big (m_1|x_{1i}|^{q_{s1}-\frac{1}{\lambda }} + m_2|x_{1i}|^{q_{s2}-\frac{1}{\lambda }} \big )^{\lambda -1} \nonumber \\&\sum _{j=1}^{2} \big ( m_j(q_{sj}-\frac{1}{\lambda })|x_{1i}|^{q_{sj}-\frac{1}{\lambda }}\big ) \Big ) {\varvec{x}}_2 + {\varvec{H}}({\varvec{x}}_1){\varvec{x}}_2 \nonumber \\&+ \lambda \text {diag} \left( |x_{2i}|^{\lambda -1} \right) \big ( - {\varvec{M}}^{-1}({\varvec{x}}_1){\varvec{C}}({\varvec{x}}_1,{\varvec{x}}_2) \nonumber \\&+ {\varvec{x}}_2 +{\varvec{\tau }}+\bar{{\varvec{d}}} \, \big ). \end{aligned}$$

(60)

Define ${\tilde{\zeta }}_1=\zeta _1-{\hat{\zeta }}_1$, ${\tilde{\zeta }}_2=\zeta _2-{\hat{\zeta }}_2$, and ${\varvec{\Lambda }}({\varvec{x}}_2)^{\lambda -1} = \text {diag} ([|x_{21}|^{\lambda -1} \; |x_{22}|^{\lambda -1} \;|x_{23}|^{\lambda -1}]^{\text {T}})$, and let a Lyapunov candidate be chosen as

$$\begin{aligned} V={\varvec{s}}^{\text {T}}{\varvec{s}}+\vartheta _0(c_1{\tilde{\zeta }}_1^2 + c_2{\tilde{\zeta }}_2^2). \end{aligned}$$

(61)

Then, using (55) and (57) yields

$$\begin{aligned} {\dot{V}}=&2{\varvec{s}}^{\text {T}}\dot{{\varvec{s}}}+2\vartheta _0c_1{\tilde{\zeta }}_1\dot{{\tilde{\zeta }}}_1+2\vartheta _0c_2{\tilde{\zeta }}_2\dot{{\tilde{\zeta }}}_2 \nonumber \\ =&-2\lambda {{\varvec{s}}^T}{{\varvec{\Lambda }}}( {{{\varvec{x}}_2}})^{\lambda - 1}(\gamma _1 \text {sig}({\varvec{s}})^{q_1} + \gamma _2 \text {sig}({\varvec{s}})^{q_2}) \nonumber \\&- 2{\vartheta _0}{c_1}{{{{\tilde{\zeta }}} }_1}{{\dot{{{\hat{\zeta }}}} }_1} - 2{\vartheta _0}{c_2}{{{{\tilde{\zeta }}} }_2}{{\dot{{{\hat{\zeta }}}} }_2} + 2\lambda \sum \limits _{i = 1}^3 { {{s_i}{{\left| {{x_{2i}}} \right| }^{\lambda - 1}}{{{\bar{d}}}_i}} } \nonumber \\&- 2\lambda ( {{{{{\hat{\zeta }}} }_1} + {{{{\hat{\zeta }}} }_2}{D_1} + {\zeta _3}{D_2}} )\sum \limits _{i = 1}^3 { {{s_i}{{\left| {{x_{2i}}} \right| }^{\lambda - 1}}{\mathop {\hbox {sgn}}} \left( {{s_i}} \right) } } \nonumber \\ \le&-2\lambda {{\varvec{s}}^T}{{\varvec{\Lambda }}}\left( {{{\varvec{x}}_2}} \right) ^{\lambda - 1}\left( \gamma _1 \text {sig}({\varvec{s}})^{q_1} + \gamma _2 \text {sig}({\varvec{s}})^{q_2} \right) \nonumber \\&- 2{\vartheta _0} ( {c_1}{{{{\tilde{\zeta }}} }_1}{{\dot{{{\hat{\zeta }}}} }_1} + {c_2}{{{{\tilde{\zeta }}} }_2}{{\dot{{{\hat{\zeta }}}} }_2} ) - 2\lambda ( {{{{{\hat{\zeta }}} }_1} + {{{{\hat{\zeta }}} }_2}{D_1} + {\zeta _3}{D_2}} ) \nonumber \\&\sum \limits _{i = 1}^3 { {{|s_i|}{{\left| {{x_{2i}}} \right| }^{\lambda - 1}}}} + 2\lambda |{\bar{d}}_{\text {max}}| \sum \limits _{i = 1}^3 { {|{s_i}|{{\left| {{x_{2i}}} \right| }^{\lambda - 1}}} } \nonumber \\ \le&-2\lambda {{\varvec{s}}^T}{{\varvec{\Lambda }}}\left( {{{\varvec{x}}_2}} \right) ^{\lambda - 1}\left( \gamma _1 \text {sig}({\varvec{s}})^{q_1} + \gamma _2 \text {sig}({\varvec{s}})^{q_2} \right) \nonumber \\&- 2{\vartheta _0} {c_1}{{{{\tilde{\zeta }}} }_1}{{\dot{{{\hat{\zeta }}}} }_1} - 2{\vartheta _0}{c_2}{{{{\tilde{\zeta }}} }_2}{{\dot{{{\hat{\zeta }}}} }_2} \nonumber \\&- 2\lambda ( {{{{{\hat{\zeta }}} }_1} + {{{{\hat{\zeta }}} }_2}{D_1} + {\zeta _3}{D_2}} )\sum \limits _{i = 1}^3 { {{|s_i|}{{\left| {{x_{2i}}} \right| }^{\lambda - 1}}} } \nonumber \\&+ 2\lambda ({{\zeta _1} + {\zeta _2}{D_1} + {\zeta _3}{D_2}}) \sum \limits _{i = 1}^3 { {|{s_i}|{{\left| {{x_{2i}}} \right| }^{\lambda - 1}}} }, \end{aligned}$$

(62)

where $|{\bar{d}}_{\text {max}}| = \text {max} \{ |\bar{d_i}|\}$ for $i=1,2,3$. Substituting (58), (59) in (62) and invoking Lemma 2 results in:

$$\begin{aligned} {\dot{V}} \le&-2\lambda {{\varvec{s}}^T}{{\varvec{\Lambda }}}( {{{\varvec{x}}_2}})^{\lambda - 1}(\gamma _1 \text {sig}({\varvec{s}})^{q_1} + \gamma _2 \text {sig}({\varvec{s}})^{q_2}) \nonumber \\&+2g_1 \tilde{\zeta _1}{\hat{\zeta }}_1 +2g_2 \tilde{\zeta _2}{\hat{\zeta }}_2 \nonumber \\ =&-2\lambda {{\varvec{s}}^T}{{\varvec{\Lambda }}}( {{{\varvec{x}}_2}})^{\lambda - 1}(\gamma _1 \text {sig}({\varvec{s}})^{q_1} + \gamma _2 \text {sig}({\varvec{s}})^{q_2}) \nonumber \\&-2g_1 {\tilde{\zeta }}_1^2-2g_2 {\tilde{\zeta }}_2^2 + 2g_1 \zeta _1 \tilde{\zeta _1} + 2g_2 \zeta _2 \tilde{\zeta _2}. \end{aligned}$$

(63)

Using Young’s inequalities, (63) can be simplified as

$$\begin{aligned} {\dot{V}} \le&-2\lambda {{\varvec{s}}^T}{{\varvec{\Lambda }}}( {{{\varvec{x}}_2}} )^{\lambda - 1}(\gamma _1 \text {sig}({\varvec{s}})^{q_1} + \gamma _2 \text {sig}({\varvec{s}})^{q_2}) \nonumber \\&- c_1 {\tilde{\zeta }}_1^2 + g_1 \varepsilon _1 \zeta _1^2 - c_2 {\tilde{\zeta }}_2^2 + g_2 \varepsilon _2 \zeta _2^2 \nonumber \\ \le&- 2\lambda {\gamma _1} \Lambda _{\text {min}} {3^{\frac{{1 - {q_1}}}{2}}}{\Big ( {\sum \limits _{i = 1}^3 {{{\left| {{s_i}} \right| }^2}} } \Big )^{\frac{{1 + {q_1}}}{2}}} - (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_1}{2}} \nonumber \\&- (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_2}{2}} - 2\lambda {\gamma _2} \Lambda _{\text {min}} {\Big ( {\sum \limits _{i = 1}^3 {{{\left| {{s_i}} \right| }^2}} } \Big )^{\frac{{1 + {q_2}}}{2}}} \nonumber \\&- (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_1}{2}} - (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_2}{2}} + (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_2}{2}} \nonumber \\&- \vartheta _0^{\frac{1+q_2}{2}} c_1 {\tilde{\zeta }}_1^2 - c_1 {\tilde{\zeta }}_1^2 + g_1 \varepsilon _1 \zeta _1^2 + \vartheta _0^{\frac{1+q_2}{2}} c_1 {\tilde{\zeta }}_1^2 \nonumber \\&+ (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_1}{2}} + g_2 \varepsilon _2 \zeta _2^2 + \vartheta _0^{\frac{1+q_2}{2}} c_2 {\tilde{\zeta }}_2^2 - c_2 {\tilde{\zeta }}_2^2 \nonumber \\&+ (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_1}{2}} + (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_2}{2}} - \vartheta _0^{\frac{1+q_2}{2}} c_2 {\tilde{\zeta }}_2^2, \end{aligned}$$

(64)

where $\Lambda _{\text {min}}$ is the minimum eigenvalue of ${\varvec{\Lambda }}({\varvec{x}}_2)^{\lambda -1}$.

Introduce $v_0= {\bar{q}}_2 ^ {\frac{{\bar{q}}_2}{1-{\bar{q}}_2} }- {\bar{q}}_2^{\frac{1}{1-{\bar{q}}_2}}$ with ${\bar{q}}_2=\frac{1+q_2}{2}$, then:

$$\begin{aligned} {\dot{V}} \le&- 2\lambda {\gamma _1} \Lambda _{\text {min}} {3^{\frac{{1 - {q_1}}}{2}}}{( {\varvec{s}}^{\text {T}}{\varvec{s}} )^{\frac{{1 + {q_1}}}{2}}} - (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_1}{2}} \nonumber \\&- 2\lambda {\gamma _2} \Lambda _{\text {min}} {({\varvec{s}}^{\text {T}}{\varvec{s}} )^{\frac{{1 + {q_2}}}{2}}} - (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_2}{2}} \nonumber \\&- (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_1}{2}} - (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_2}{2}} + \vartheta _1 \nonumber \\&+ (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_1}{2}} + (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_1}{2}} \nonumber \\&- \vartheta _0c_1 {\tilde{\zeta }}_1^2 - \vartheta _0c_2 {\tilde{\zeta }}_2^2 \nonumber \\ \le&-\eta _1V^{\frac{1+q_1}{2}} - \eta _2 V^ {\frac{1+q_2}{2}} + \vartheta _1 + \mathrm {\Delta }_{\zeta }, \end{aligned}$$

(65)

where $\vartheta _1 = g_1 \varepsilon _1 \zeta _1^2 + g_2 \varepsilon _2 \zeta _2^2 + 2 \vartheta _0 ^ {{\bar{q}}_2} v_0$, $\Delta _{\zeta } = (\vartheta _0 c_1 {\tilde{\zeta }}_1^2)^{\frac{1+q_1}{2}} + (\vartheta _0 c_2 {\tilde{\zeta }}_2^2)^{\frac{1+q_1}{2}} - \vartheta _0c_1 {\tilde{\zeta }}_1^2 - \vartheta _0c_2 {\tilde{\zeta }}_2^2$, $\eta _1 = \text {min} \{ 2\lambda {\gamma _1} \Lambda _{\text {min}}, 1 \}$, and $\eta _2 = \text {min} \{ 2\lambda {\gamma _2} \Lambda _{\text {min}}, 1 \}$.

Assuming that an unknown constant $\varpi$ and a compact set ${\mathbb {U}}$ exist and satisfy ${\mathbb {U}}= \{ ({\tilde{\zeta }}_1,{\tilde{\zeta }}_2) \, | \, {\tilde{\zeta }}_1 \le \varpi , {\tilde{\zeta }}_2 \le \varpi \}$, it leads to

$$\begin{aligned} {\dot{V}} \le -\eta _1V^{\frac{1+q_1}{2}} - \eta _2 V^ {\frac{1+q_2}{2}} + \vartheta _2, \end{aligned}$$

(66)

where $\vartheta _2=\vartheta _1$, if $\varpi \le \min \left\{ {\frac{1}{{\sqrt{{\vartheta _0}{c_1}} }},\frac{1}{{\sqrt{{\vartheta _0}{c_2}} }}} \right\}$; $\vartheta _2={\vartheta _1} + 2\big ( {{{( {{\vartheta _0}\max \left\{ {{c_i}} \right\} {\varpi ^2}} )}^{\frac{{1 + {q_1}}}{2}}} - {\vartheta _0}\max \left\{ {{c_i}} \right\} {\varpi ^2}}\big )$, if $\varpi > \min \left\{ {\frac{1}{{\sqrt{{\vartheta _0}{c_1}} }},\frac{1}{{\sqrt{{\vartheta _0}{c_2}} }}} \right\}$.

According to Theorem 2, it can be proved from (66) that the closed-loop attitude stabilization system is practically stable. V will converge to the residual set ${\mathbb {D}}_1 = \{ \mathop {\lim }\limits _{t \rightarrow {T_V}} \,|\, V \le \text {min} \big \{ (\frac{\vartheta _2}{\eta _1(1-\mu )})^{\frac{2}{1+q_1}} \;,\; (\frac{\vartheta _2}{\eta _2(1-\mu )})^{\frac{2}{1+q_2}} \} \big \}$ and the settling time is given by

$$\begin{aligned} T_V \le \max&\left\{ \frac{2}{\eta _2(1-q_2)}+\frac{2V(0)^{\frac{1-q_1}{2}}}{\eta _1\mu (q_1-1)} \;, \right. \nonumber \\&\left. \frac{2}{\eta _2\mu (1-q_2)}+\frac{2V(0)^{\frac{1-q_1}{2}}}{\eta _1(q_1-1)} \right\} , \end{aligned}$$

(67)

where $0<\mu <1$ is a constant. Hence, the settling time of the attitude stabilization system is guaranteed by $T \le T_s+T_V$. Moreover, (66) satisfies the form of faster finite-time system; this completes the proof. $\square$

4 Numerical simulation

To verify the control performance of the proposed faster finite-time controller (55), numerical simulations are conducted where the nominal inertia of the receiver aircraft is ${{\varvec{J}}_0} = \left[ {\begin{array}{*{20}{c}} {500}&{}\quad {0}&{}\quad {-100}\\ {0}&{}\quad {500}&{}\quad {0}\\ {-100}&{}\quad {0}&{}\quad {500} \end{array}} \right] \text {kg} \cdot \text {m}^2$ and the uncertain inertia ${\varvec{J}}_{\Delta }=0.05{\varvec{J}}_0$. The correlation coefficient in the time-varying inertia (7) is selected as ${\varvec{\rho }}=[1\quad 0.8 \quad 0.9 ]^{\text {T}}$, and the fuel transferring speed is assumed as $m_{\text {var}}= 25$ kg/sec. The turbulence is generated through the white noise generator and the shaping filters shown in Fig. 3. The disturbance torque ${\varvec{d}}_{\text {wind}}$ is calculated as stated in Sect. 2.3. The additional aerodynamic moment ${\varvec{d}}_{\text {aero}}$ is obtained using (25) with the shift of CM limiting in $\pm \;4 \%$ chord. The initial states of the receiver aircraft are set as ${\varvec{\sigma }}_{0}= [10\quad 2\quad -5]^{\mathrm{{T}}}$ degrees and ${\varvec{\omega }}_{b0}= [2\quad -5\quad 6]^{\text {T}}$ deg/sec. Aerial refueling is conducted in the height of 12000 meters, and the velocity for aircraft is set as 100 m/sec.

Fortran is selected as the programming language to carry out the simulation. ODEPACK and BLAS packages are mainly utilized to solve the differential equations and perform vector and matrix operations. Since numerical solutions of the system is important [38, 39], we give Algorithm 1 to briefly introduce the steps of the numerical solutions of attitude stabilization system. In simulations, the proposed approach (named Faster FTC) is compared with the conventional fast finite-time control law [24] (named Fast FTC). The control gains of both controllers are selected as $m_1=6$, $m_2=6$, $q_1=2$, $q_{s1}=1$, $q_2=0.8$, $q_{s2}=0.8$, $\lambda =1.5$, $g_1=5$, $g_2=5$, $\varepsilon _1=0.51$, $\varepsilon _2=0.51$, $\gamma _1=5$, and $\gamma _2=10$.

4.1 Simulation results of the Faster FTC

During the refueling process, light turbulence is assumed because the aerial refueling task often occurs in good weather condition. Then, the receiver aircraft is exposed to the multiple flow perturbations (turbulence and vortex) shown in Fig. 5, which is to deteriorate the attitude stabilization control performance. After applying the Faster FTC, the resulting attitude stabilization performance is illustrated in Figs. 6 and 7. It is seen that the proposed controller (55) ensures a fast convergence rate even in the presence of wind perturbations. This verifies the conclusion in Theorem 4.

Fig. 5

The airflows that the receiver aircraft exposed

As shown in Figs. 6 and 7, the control accuracy of the attitude and the angular velocity are $\pm \, 0.2$ deg and $\pm \, 0.1 \,\text {rad}$ /sec, respectively, which is inferior to the results in [30, 40]. That is because the disturbance $\bar{{\varvec{d}}}$ in this paper is related to the high-frequency, non-smooth wind disturbance, while the continuous sine/cosine functions were used in [30, 40] to represent the external disturbance. Moreover, the magnitude of $|| \bar{{\varvec{d}}} ||$ ranges from $10^3$ to $10^4$ Nm. It is larger than the magnitude 1 Nm in [30] and $10^{-2}$ Nm in [40]. The control effort ${\varvec{\tau }}$ consumed to stabilize the attitude system is shown in Fig. 8. The required control torques are feasible through aerodynamic effects of aileron, elevator, and rudder. Hence, the control scheme is a potential basis for engineering practice.

Fig. 6

Time responses of the attitude

Fig. 7

Time responses of the angular velocity

Fig. 8

Overall control torque

4.2 Comparison with the Fast FTC

In this subsection, the faster convergence rate of Faster FTC compared with conventional Fast FTC is verified. Simulations are divided into two parts. Firstly, without loss of generality, we omit the lumped disturbance $\bar{{\varvec{d}}}$ and cancel the ${\varvec{u}}_\text {adp}$ for better illustrating the influence of parameter $q_1$. Then, practical faster characteristics using the adaptive law ${\varvec{u}}_{\text {adp}}$ is demonstrated as well.

4.2.1 Comparison of convergence rate in the absence of disturbance

Figure 9 shows the comparative results between Faster FTC (i.e., $q_1 > 1$) and conventional FTC (i.e., $q_1=1$) when $q_{s1}$ is fixed as 1 in the sliding mode surface (49). In both simulations, control gains and system parameters remain the same besides the value of $q_1$. We can observe that when Faster FTC is applied, the convergence rates in three channels increase and the settling time decrease, which offers faster attitude stabilization capability for receiver aircraft. This also verifies the motivation of the development of the faster finite time stable system in (33). The underlying reason of the faster rate is related to the decrease in approaching time of the sliding mode surface. Taking the roll angle $\phi$ as an example, Fig. 10 demonstrates the first component in the sliding mode surface ${\varvec{s}}$, which reveals that larger $q_1$ improves the convergence rate.

Fig. 9

Comparisons of attitude convergence speed when $q_{s1}=1$

However, it should be pointed out that Theorem 1 shows time difference ${\tilde{T}}$ between Faster FTC and Fast FTC is determined by $l(\xi _2)$. If we choose $q_{s1} > 1$ in the sliding manifold (49), the completely opposite tendency of convergence rate occurs, as shown in Fig. 11. Larger $q_1$ and $q_{s1}$ lead to slower convergence rate. This is not a counter example of the proposed faster finite-time stable system. It is a numerical verification of Theorem 1 from the other perspective. That is because in the approaching phase, larger $q_1$ will lead to shorter convergence time of ${\varvec{s}}$. However, when the sliding mode surface reaches the equilibrium point, the attitude ${\varvec{\sigma }}$ has already converged into the small set containing ${\varvec{0}}$, as we can see in Fig. 9. Since ${\varvec{x}}_{1}(0)$ is relatively small at that time, the faster convergence rate is not seen in the sliding mode phase according to Theorem 1.

Fig. 10

Comparisons of time response of sliding manifold when $q_{s1}=1$

Fig. 11

Comparisons of attitude convergence speed when $q_{s1}=q_1>1$

4.2.2 Comparison of convergence rate with disturbance

After analyzing the relationship between the parameter $q_1$ in Faster FTC and convergence rate under disturbance-free environment, the practical faster attitude stabilization effect with the lumped disturbance considered is shown in this subsection.

Due to the lumped disturbance $\bar{{\varvec{d}}}$, deterioration of the faster convergence rate of the Faster FTC occurs, as shown in Figs. 12 and 13. When $\bar{{\varvec{d}}}$ is considered, the adaptive control effort is compensated to handle the unknown disturbance. It affects the nominal convergence performance of the attitude stabilization system illustrated in Sect. 4.2.1 . Only a slight improvement in the convergence rate can be observed in the attitude and the sliding mode surface. Nevertheless, those results still authenticate the conclusion in Theorem 4 that the proposed controller could not only attenuate the unknown lumped disturbance adaptively but also decrease the settling time when compared with the conventional fast finite-time control schemes.

Fig. 12

Comparisons of attitude convergence speed with disturbance considered

Fig. 13

Comparisons of sliding manifold convergence speed with disturbance considered

5 Conclusion

Although there exist several approaches regarding the attitude control for receiver aircraft in the aerial refueling task, few of them could achieve finite-time convergence of the system states in the presence of multiple disturbances. In this paper, the attitude dynamics with uncertain and time-varying inertia, wind disturbances and the shift of CM was first established . On this basis, a nonsingular sliding mode-based adaptive control law was synthesized to guarantee the faster finite-time stability of the closed-loop system. The main feature of the proposed control strategy is smaller convergence time in comparison with the conventional fast finite-time control scheme; meanwhile, the rigorous proof was also provided. The unknown disturbance could be adaptively adopted without requiring its upper bound. As some of the future work, actuator constraints like input saturation and the transplantation of this control approach to commercial/open-source autopilots can be considered.