Fault estimation and compensation for hypersonic flight vehicle via type-2 fuzzy technique and cuckoo search algorithm

Abstract

This study proposes a fault estimation and compensation scheme for hypersonic flight vehicle (HFV) attitude system with body flap fault. To achieve the fault compensation capability, an interval type-2 fuzzy estimator is designed to approximate the unknown nonlinear function caused by modeling uncertainties and the fault. The fault-tolerant controller via the terminal sliding mode and dynamic surface techniques is developed to ensure the tracking accuracy of attitude angles for HFV and the desired control torque can be achieved. Then, a novel control allocation scheme with a cuckoo search algorithm and linear programming method is proposed to distribute the desired torque to aerodynamic surfaces and reaction control system jets. The stability of the proposed scheme is analyzed using the Lyapunov stability theory. The validity of the method is verified by a series of comparisons on numerical simulation results.

Introduction

Hypersonic flight vehicles (HFVs) refer to vehicles that travel at a velocity greater than 5 Mach. ¹ Such technique provides a reliable and cost-effective means for access to near space and has attracted great attention owing to their potential civilian and military applications such as prompt global response. ^2,3 Since HFV works in a harsh environment where flight conditions and parameters often change radically, faults may occur in actuators such as elevators and rudders. The high cost of research and development of HFV also means that the failure of flight test will cause immeasurable loss of property. ⁴ Therefore, the fault-tolerant control (FTC) of HFV under actuator faults is highly significant.

In recent years, a variety of control methods have been proposed to deal with the adverse effects due to the nonlinearity and actuator faults in HFV. ^5,6 Adaptive FTC is a common approach, which can make the flight state track the given command ^7,8 A fault compensation scheme was investigated through a finite-time observer which was used to estimate actuator faults and uncertainties for hypersonic gliding vehicles. ⁹ Sliding mode control-based FTC was developed to ensure the tracking performance via suitable sliding mode surface. ^10,11 Back-stepping fault-tolerant output tracking control was utilized for HFV subject to actuator constraints, disturbances, and uncertainties. ¹² To handle the problem of parameter uncertainties or unknown parameters, neural network was introduced to enhance the robustness of the controller. ^13,14 An adaptive output feedback FTC was developed for the longitudinal dynamics of HFV in the presence of parameter uncertainties, actuator faults, and external disturbances. ¹⁵ In addition to actuator faults, the modeling uncertainties are inevitably rooted in the control-oriented models due to the lack of wind tunnel test data at hypersonic speeds and the difficulty to obtain sufficient on-flight test data. ^16,17 Therefore, the fuzzy estimation algorithm that has high approximation capability for identifying nonlinear systems is considered to calculate the unknown effects of modeling error and fault. ¹⁸ However, the membership functions are usually defined according to the experience of designers and difficult to be calculated precisely. Compared with the conventional fuzzy systems, interval type-2 fuzzy system uses an interval rather than a number in their membership functions to define the membership degree. Thus, it provides more degrees of freedom that could reduce the effects of uncertainty in the fuzzy rules. ^{19 –21} Given such advantages, an interval type-2 fuzzy estimation algorithm is designed to represent the unknown adverse effects of modeling error and faults in this article.

The estimation results are further used in the design of the FTC to obtain the needed control torque. Such control torque will be allocated to eight aerodynamic surfaces. However, these aerodynamic surfaces are usually too weak to satisfy control commands. ²² Thus, reaction control system (RCS) will be introduced to provide the residual control torque. ²³ In the literature, ¹⁰ the derived control torque is distributed by pseudoinverse method directly. With this allocator, the utilization ratio of aerosurfaces is definitely low because of the fixed calculation. This has resulted in fast energy consumption of the RCS system. A new feature of this article is the utilization of cuckoo search algorithm (CSA) in torque allocation part. The CSA is a relatively novel and promising metaheuristic algorithm. ^24,25 It was inspired by the obligate brood parasitism of some cuckoo species by laying their eggs in the nests of other host birds (of other species). Cuckoo search idealized its breeding behavior and thus can be applied for various optimization problems such as face recognition, ²⁶ structural optimization, ²⁷ optimal placement and sizing of static volt-ampere reactive (VAR) compensator, ²⁸ and irrigation station optimization. ²⁹ Some studies have demonstrated that its search capability is better than many other metaheuristic algorithms. ^30,31 For HFVs, the control allocation problem is described as a linear programming problem and solved by CSA in our work.

This article is organized as follows. The attitude model of HFV and FTC scheme is first described in the second section. The third section presents an interval type-2 estimation algorithm and then develops an FTC via terminal sliding mode and dynamic surface control techniques. The fourth section proposes a control allocation algorithm which combines the linear programming problem and CSA. Several contrast simulations are conducted to demonstrate the validity of the proposed fuzzy fault estimation and allocation schemes in the fifth section. Finally, conclusions are presented at the end of this article.

Problem formulation

Reentry attitude dynamic model

The most common attitude dynamic model of HFV is shown as follows ³²

$\{\begin{array}{l}\dot{\omega }=-J^{-1}\Omega J\omega +J^{-1}T+g(\omega (t))\hfill \\ \dot{\theta }=R\omega +d\hfill \end{array}$

1

where $\theta \text{=[}\varphi ,\alpha ,\beta {\text{]}}^T$ is the attitude angle vector, with $\varphi$ , α, and $\beta$ as the bank angle, angle of attack, and sideslip angle, respectively; $\omega \text{=[}p,q,r{\text{]}}^T$ is the angular rate vector, with p, q, and r as roll rate, pitch rate, and yaw rate, respectively; $J\in R^{3\times 3}$ is the symmetric positive definite moment of inertia tensor; $T\in R^{3\times 1}$ is the control torque vector; $g(\omega (t))$ is an unknown nonlinear function caused by modeling error and external disturbance; and d denotes the external disturbance. The other matrices $\Omega$ and R are given as follows

$\Omega \text{=}\left[\begin{array}{ccc}0& -r& q\\ r& 0& -p\\ -q& p& 0\end{array}\right], R=\left[\begin{array}{ccc}\text{cos}\alpha & 0& \text{sin}\alpha \\ \text{sin}\alpha & 0& -\text{cos}\alpha \\ 0& 1& 0\end{array}\right]$

2

The control torque T is generated by 4 types of control surfaces and 10 RCS jets, which can be described by

$T=T_{\delta }+T_{\text{RCS}}$

3

The torque $T_{\delta }$ produced by aerodynamic surfaces is expressed as follows

$T_{\delta }\text{=}B{\delta }_e$

4

where ${\delta }_e={[{\delta }_{\text{rei}}{\delta }_{\text{lei}}{\delta }_{\text{rft}}{\delta }_{\text{lft}}{\delta }_{\text{lvr}}{\delta }_{\text{lvr}}{\delta }_{\text{reo}}{\delta }_{\text{leo}}]}^T$ with ${\delta }_{\text{rei}},\cdots , {\delta }_{\text{leo}}$ representing right and left inboard elevons, right and left body flaps, right and left rudders, and right and left outboard elevons, respectively; and $B\in R^{3\times 8}$ is a sensitivity matrix. ³⁰

The torque $T_{\text{RCS}}$ produced by RCS jets is described as follows

$T_{\text{RCS}}=\Psi u_{\text{rcs}}$

5

where $\Psi \in R^{3\times 10}$ is a weight matrix, and its elements ${\Psi }_{ij}$ represents the control torque generated by the $j\text{th}$ RCS on the $i\text{th}$ axis. $u_{\text{rcs}}\text{=}\left[u_1,\cdots ,u_{10}\right]$ is the command vector and $u_k\text{=}0/1$ ( $k=1,2,\cdots ,10$ ) represents the off/on state of the $k\text{th}$ RCS jet.

Actuator fault model

Aerodynamic surface faults may occur since HFV faces harsh environment in reentry phase, resulting in that the actual control torque cannot reach the desired one. Body flap is one of the actuators, which has great impact on the control torque. Therefore, HFV model with a right body flap fault is considered in this study, which is described as follows

$\{\begin{array}{l}\dot{\omega }=-J^{-1}\Omega J\omega +J^{-1}T+g(\omega (t))+J^{-1}B{f}_{\delta }\hfill \\ \dot{\theta }=R\omega +d\hfill \end{array}$

6

where $f_{\delta }\in R^{8\times 1}$ represents the additive fault of the actuator and its third element is not always zero and other elements always zero.

Assumption 1

The external disturbance d is bounded and its derivative is also bounded. That is, there exist two positive scalars d ₁ and d ₂ such that $\Vert d\Vert \le d_1$ and $\Vert \dot{d}\Vert \le d_2$ .

The control objective is to design an FTC law with respect to the control torque generated by aerodynamic surfaces and RCS jets such that the attitude angles can track the given commands asymptotically even if the body flap fault and disturbances exist. Normally, the whole attitude control system is divided into two subsystems: the inner angular rate loop and the outer attitude angle loop. The overall control scheme is shown in Figure 1.

Figure 1.

Block diagram of the FTC scheme. FTC: fault-tolerant control.

Design of FTC

Interval type-2 fuzzy estimation

As long as there are enough fuzzy conditions, any nonlinear function can be approximated with arbitrary precision. Therefore, a type-2 fuzzy estimation algorithm is designed to approximate the lumped effect of the unknown nonlinear function $g(\omega (t))$ and the actuator fault $f_{\delta }$ in this subsection. With reference to Rahmani and Shahriari-Kahkeshi, ²⁰ the fuzzy rules of the estimation algorithm can be written as follows

$\begin{array}{l}{R}^{\left(j\right)}: \text{If} x_1 \text{is} {\stackrel{\sim }{F}}_1^j \text{and}\cdots \text{and} x_n is {\stackrel{\sim }{F}}_n^j\\ \text{then} h \text{is} G^j, j=1,2,\cdots ,M\end{array}$

7

where $x_1,x_2,\cdots ,x_n$ are input variables, h is the output variable, and M represents the total number of fuzzy rules. ${\stackrel{\sim }{F}}_i^j\left(i=1,2,\cdots ,n\right)$ and G^j represent interval type-2 antecedent and the consequent fuzzy sets, respectively. Through combining the firing strength of each fuzzy rule, the type-2 inference engine performs a nonlinear map from input interval type-2 fuzzy set to the output interval type-2 fuzzy set. The upper and lower bounds of the firing strength of the $P\text{th}$ rule can be calculated as follows

$\begin{array}{l}{\overline{f}}^P={\overline{\mu }}_{F_1^P}(x_1)\times {\overline{\mu }}_{F_2^P}(x_2)\times \mathrm{...}\times {\overline{\mu }}_{F_n^P}(x_n)\\ {\underline{f}}^P={\underline{\mu }}_{F_1^P}(x_1)\times {\underline{\mu }}_{F_2^P}(x_2)\times \mathrm{...}\times {\underline{\mu }}_{F_n^P}(x_n)\end{array}$

8

where ${\overline{\mu }}_{F_i^P}(x_i)$ and ${\underline{\mu }}_{F_j^P}(x_j)$ represent the upper and lower membership functions corresponding to the $P\text{th}$ fuzzy rule, respectively. The defuzzified crisp point can be calculated via the fuzzy weights as follows

$h_l=\frac{{\displaystyle {\sum }_{P=1}^M{\underline{f}}^P{h}_l^P}}{{\displaystyle {\sum }_{P=1}^M{\underline{f}}^P}},h_u=\frac{{\displaystyle {\sum }_{p=1}^M{\overline{f}}^P{h}_u^P}}{{\displaystyle {\sum }_{P=1}^M{\overline{f}}^P}}$

9

For convenience of analysis, h_l and h_u can be written as the linear regression form as follows

$h_l={\phi }_l^T{\xi }_l(x),h_u={\phi }_u^T{\xi }_u(x)$

10

where ${\phi }_l={[h_l^1{h}_l^2\cdots h_l^M]}^T,{\xi }_l(x)={[{\xi }_l^1(x){\xi }_l^2(x)\cdots {\xi }_l^M(x)]}^T$ , ${\phi }_u={[h_u^1{h}_u^2\cdots h_u^M]}^T,\hspace{0.17em}{\xi }_u(x)={[{\xi }_u^1(x){\xi }_u^2(x)\cdots {\xi }_u^M(x)]}^T$ , and ${\xi }_u(x)={[{\xi }_u^1(x){\xi }_u^2(x)\cdots {\xi }_u^M(x)]}^T$ are the fuzzy basis functions derived by the following methods

${\xi }_l^j(x)=\frac{{\underline{f}}^j}{{\displaystyle {\sum }_{P=1}^M{\underline{f}}^P}}, {\xi }_u^j(x)=\frac{{\overline{f}}^j}{{\displaystyle {\sum }_{P=1}^M{\overline{f}}^P}}$

11

where $j=1,2,\cdots ,M$ . Then, the final defuzzified crisp point is given by

$\begin{array}{l}h(x)=[{\phi }_l^T{\xi }_l(x)+{\phi }_u^T{\xi }_u(x)]/2\\ ={\phi }_l^T{\xi }_l^{*}(x)+{\phi }_u^T{\xi }_u^{*}(x)\end{array}$

12

where ${\xi }_l^{\text{*}}(x)={\xi }_l(x)/2$ and ${\xi }_u^{*}(x)={\xi }_u(x)/2$ .

Let

$G_f\left(x\right)=g(\omega (t))+J^{-1}B(\cdot )f_{\delta }$

13

For any constant $\delta >0$ and the unknown nonlinear term $G_f(x)$ in the angular rate loop, there exists the following fuzzy logic system

${\widehat{G}}_f(x)={\phi }_{lf}^{*T}{\xi }_{lf}^{*}(x)+{\phi }_{uf}^{*T}{\xi }_{uf}^{*}(x)$

14

such that

$\underset{x\in U}{\text{sup}}|{\widehat{G}}_f(x)-G_f(x)|\le \delta$

15

where ${\phi }_{lf}^{*}$ and ${\phi }_{uf}^{*}$ are optimal parameter matrices, and ${\xi }_{lf}^{*}(x)$ and ${\xi }_{uf}^{*}(x)$ are fuzzy basis functions. Then, $G_f(x)$ can be written as follows

$G_f(x)={\phi }_{lf}^{*T}{\xi }_{lf}^{*}(x)+{\phi }_{uf}^{*T}{\xi }_{uf}^{*}(x)+{\epsilon }_f$

16

where ${\epsilon }_f$ is the error between the original unknown function and the estimation function obtained by using the optimal parameter matrices. However, since those optimal parameter matrices are difficult to obtain, they are replaced by their estimation values ${\widehat{\phi }}_{lf}$ and ${\widehat{\phi }}_{uf}$ . Then, the fuzzy estimation result is given by

${\widehat{G}}_f(x)={{\widehat{\phi }}_{lf}}^T{\xi }_{lf}^{*}(x)+{{\widehat{\phi }}_{uf}}^T{\xi }_{uf}^{*}(x)$

17

where the adaptive updating laws of ${\widehat{\phi }}_{lf}$ and ${\widehat{\phi }}_{uf}$ will be derived from the Lyapunov function in the next section.

Outer loop terminal sliding mode controller

For the attitude control of HFV, system outputs are three attitude angles which are actually desired to track some predefined reference commands. Three angular rates, which can affect the system outputs, are treated as interval variables. And normally, the control torque is the system input. To obtain the needed control torque, it is necessary to calculate the corresponding angular rates via the system outputs. Thus, we can choose the following terminal sliding mode surface

$S=\stackrel{\sim }{\theta }+{\displaystyle {\int }_0^t[K_1{|\stackrel{\sim }{\theta }|}^{{\eta }_1}\text{sgn}(\stackrel{\sim }{\theta })+K_2{|\stackrel{\sim }{\theta }|}^{{\eta }_2}\text{sgn}(\stackrel{\sim }{\theta })]}d\tau$

18

$\stackrel{\sim }{\theta }=\theta -{\theta }_d$

19

where ${\eta }_1\ge 1,0<{\eta }_2<1$ ; ${\theta }_d$ denotes the desired attitude angle vector; $K_i=\text{diag}\left\{k_{i1},k_{i2},k_{i3}\right\}(i=1,2)$ is the diagonal outer loop gain matrix. And let $|\stackrel{\sim }{\theta }|=\text{diag}\left\{|{\stackrel{\sim }{\theta }}_1|,|{\stackrel{\sim }{\theta }}_2|,|{\stackrel{\sim }{\theta }}_3|\right\}$ .

The reaching law of this terminal sliding mode surface is designed with a parameter $0<{\eta }_3<1$ as follows

$\dot{S}=-{\epsilon }_{1}S-{\epsilon }_2{|S|}^{{\eta }_3}\text{sgn}\left(S\right)-\stackrel{\sim }{d}$

20

where ${\epsilon }_1$ and ${\epsilon }_2$ are two positive scalars. The virtual control input can be derived as follows

$\begin{array}{l}{\omega }_c=R^{-1}[{\dot{\theta }}_d-K_1{|\stackrel{\sim }{\theta }|}^{{\eta }_1}\text{sgn}(\stackrel{\sim }{\theta })-K_2{|\stackrel{\sim }{\theta }|}^{{\eta }_2}\text{sgn}(\stackrel{\sim }{\theta })\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\epsilon }_{1}S-{\epsilon }_2{|S|}^{{\eta }_3}\text{sgn}(S)-\widehat{d}]\end{array}$

21

where $\widehat{d}$ represents the estimation of the external disturbance d.

The adaptive law of $\widehat{d}$ is shown as follows

$\dot{\widehat{d}}={\Lambda }_{1}S$

22

where ${\Lambda }_1=\text{diag}\left\{{\lambda }_{11},{\lambda }_{12},{\lambda }_{13}\right\}$ is the adaptive learning rate.

Theorem 1

The terminal sliding mode surface S can converge to an expected region and then $\stackrel{\sim }{\theta }$ is uniformly ultimately bounded under the adaptive law expressed in equation (22).

Proof 1

Define the disturbance estimation error $\stackrel{\sim }{d}=\widehat{d}-d$ and the following Lyapunov function is selected

$V_S=\frac{1}{2}{S}^{T}S+\frac{1}{2}{\stackrel{\sim }{d}}^T{\Lambda }_1^{\text{-}1}\stackrel{\sim }{d}$

23

The derivative of V_S is finally shown as follows

$\begin{array}{l}{\dot{V}}_s=-S^T[{\epsilon }_{1}S\text{+}{\epsilon }_2|S{|}^{{\eta }_3}\text{sgn}(S)+\stackrel{\sim }{d}]+{\stackrel{\sim }{d}}^T{{\Lambda }_1}^{-\; 1}\dot{\stackrel{\sim }{d}}\\ =-S^T[{\epsilon }_{1}S\text{+}{\epsilon }_2|S{|}^{{\eta }_3}\text{sgn}(S)+\stackrel{\sim }{d}]+{\stackrel{\sim }{d}}^T{{\Lambda }_1}^{-\; 1}(\dot{\widehat{d}}-\dot{d})\\ =-{\epsilon }_1||S|{|}^2-{\epsilon }_2||S|{|}^{{\eta }_3+1}-{\stackrel{\sim }{d}}^T{{\Lambda }_1}^{-\; 1}\dot{d}\\ \le -{\epsilon }_1||S|{|}^2-{\epsilon }_2||S|{|}^{{\eta }_3+1}+\frac{1}{{\lambda }_{\text{min}}({\Lambda }_1)}||{\stackrel{\sim }{d}}^T\dot{d}||\\ \le -{\epsilon }_1||S|{|}^2-{\epsilon }_2||S|{|}^{{\eta }_3+1}+ \frac{\kappa }{{\lambda }_{\text{min}}({\Lambda }_1)}\end{array}$

24

where $\kappa =||{\stackrel{\sim }{d}}^T\dot{d}||$ and ${\lambda }_{\text{min}}(\cdot )$ denote the minimum eigenvalue of the corresponding matrix. Then, the system can converge to the following interval

$||S||\le \text{min}\left\{\begin{array}{cc}{\left(\frac{\kappa }{{\epsilon }_1{\lambda }_{\text{min}}(\Lambda )}\right)}^{1/2},& {\left(\frac{\kappa }{{\epsilon }_2{\lambda }_{\text{min}}(\Lambda )}\right)}^{1/({\eta }_3+1)}\end{array}\right\}$

25

Therefore, when appropriate parameters ( ${\epsilon }_1$ , ${\epsilon }_2$ , and $\Lambda$ ) are adopted, the system can converge to an expected region. That is, $\stackrel{\sim }{\theta }$ is uniformly ultimately bounded with the virtual control (21).

Here, we only give simple information. Refer to the work of Wang et al. ³³ for the detailed information of the design and certification process of the outer loop controller.

Remark 1

The sign function in equation (21) is a discontinuous function with a significantly high-frequency switching, which may result in the chattering of the angular rate command signal. Therefore, the sign function can be approximately replaced by the following continuous function

$\text{sat}\left(\frac{S_i}{{\vartheta }_i}\right)=\{\begin{array}{ll}1& \text{if\hspace{0.17em}}{S}_i>{\vartheta }_i \\ \frac{S_i}{{\vartheta }_i}& \text{if\hspace{0.17em}}|S_i|\le {\vartheta }_i\\ -1& \text{if\hspace{0.17em}}{S}_i<{\vartheta }_i\end{array}\text{, }i=1,2,3$

26

where ${\vartheta }_i\left(i=1,2,3\right)$ are adjustable parameters which can be chosen in the form shown as follows

${\vartheta }_i=\{\begin{array}{l}1/\left({\omega }_{c_i}+\upsilon \right),{\omega }_{c_i}\ge \text{0}\\ -1/\left({\omega }_{c_i}-\upsilon \right), {\omega }_{c_i}<\text{0 }\end{array}$

27

A negative scalar $\upsilon$ which is extremely small needs to be added here to avoid the appearance of the singular point since the value of the virtual control may be equal to zero in actual applications.

Inner loop dynamic surface controller

When the virtual control input is determined, the desired angular rate command can be derived by introducing a first-order filter. Then, the dynamic surface control can be utilized to calculate the needed control torque and stabilize the inner loop angular rate system.

Let the virtual control input (21) pass through the following first-order filter

$\{\begin{array}{l}\tau {\dot{\omega }}_d+{\omega }_d={\omega }_c\\ {\omega }_c(0)={\omega }_d(0)\end{array}$

28

where $\tau >0$ denotes the time constant. Then

$\begin{array}{l}\dot{\stackrel{\sim }{\theta }}=\dot{\theta }-{\dot{\theta }}_d\\ =R\left(\omega -{\omega }_d+{\omega }_d-{\omega }_c+{\omega }_c\right)+d-{\dot{\theta }}_d\\ =R\left(\stackrel{\sim }{\omega }+\overline{\mu }+{\omega }_c\right)+d-{\dot{\theta }}_d\\ =-K_1{|\stackrel{\sim }{\theta }|}^{{\eta }_1}\text{sgn}(\stackrel{\sim }{\theta })-K_2{|\stackrel{\sim }{\theta }|}^{{\eta }_2}\text{sgn}(\stackrel{\sim }{\theta })-{\epsilon }_1{S}_1\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\epsilon }_2{|S_1|}^{{\eta }_3}\text{sgn}(S_1)+\stackrel{\sim }{d}+R\stackrel{\sim }{\omega }+R\overline{\mu }\end{array}$

29

where $\overline{\mu }={\omega }_d-{\omega }_c$ represents the filter error.

Define the angular rate tracking error vector as follows

$\stackrel{\sim }{\omega }\hspace{0.17em}=\hspace{0.17em}\omega \hspace{0.17em}-\hspace{0.17em}{\omega }_d$

30

$\begin{array}{l}\dot{\stackrel{\sim }{\omega }}\hspace{0.17em}=\hspace{0.17em}\dot{\omega }\hspace{0.17em}-\hspace{0.17em}{\dot{\omega }}_d\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}-J^{-1}\hspace{0.17em}\Omega J\omega \hspace{0.17em}+\hspace{0.17em}{J}^{-1}T\hspace{0.17em}+\hspace{0.17em}{G}_f\hspace{0.17em}-\hspace{0.17em}{\dot{\omega }}_d\end{array}$

31

Considering equations (17) and (31), the FTC for the angular rate loop can be designed as follows

$\begin{array}{l}T=-J[{\Lambda }^{*}\stackrel{\sim }{\omega }-J^{-1}\Omega J\omega +R^T\stackrel{\sim }{\theta }\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\dot{\omega }}_d+({{\widehat{\phi }}_{lf}}^T{\xi }_{lf}^{*}+{{\widehat{\phi }}_{uf}}^T{\xi }_{uf}^{*})]\end{array}$

32

where ${\Lambda }^{*}=\text{diag}\{{\lambda }_1^{*},{\lambda }_2^{*},{\lambda }_3^{*}\}({\lambda }_i^{*}>0, i=1,2,3)$ is the diagonal gain matrix to be designed.

Substituting equation (32) into equation (31), we can obtain

$\dot{\stackrel{\sim }{\omega }}=-{\Lambda }^{*}\stackrel{\sim }{\omega }-R^T\theta +{{\stackrel{\sim }{\phi }}_{lf}}^T{\xi }_{lf}^{*}+{{\stackrel{\sim }{\phi }}_{uf}}^T{\xi }_{uf}^{*}+{\epsilon }_f$

33

where ${\stackrel{\sim }{\phi }}_{lf}={\phi }_{lf}^{*}-{\widehat{\phi }}_{lf}$ and ${\stackrel{\sim }{\phi }}_{uf}={\phi }_{uf}^{*}-{\widehat{\phi }}_{uf}$ .

For the whole HFV attitude system, the Lyapunov function candidate is chosen as follows

$\begin{array}{l}V=\underset{V_1}{\underbrace{\frac{1}{2}{\stackrel{\sim }{\theta }}^T\stackrel{\sim }{\theta }+\frac{1}{2}{\stackrel{\sim }{d}}^T{\Lambda }^{-1}\stackrel{\sim }{d}+\frac{1}{2}{\overline{\mu }}^T\overline{\mu }}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{+}\underset{V_2}{\underbrace{\frac{1}{2}{\stackrel{\sim }{\omega }}^T\stackrel{\sim }{\omega }+\frac{1}{2}tr\left\{{{\stackrel{\sim }{\phi }}_{lf}}^T{{\Upsilon }_{lf}}^{-1}{\stackrel{\sim }{\phi }}_{lf}\right\}+\frac{1}{2}tr\left\{{{\stackrel{\sim }{\phi }}_{uf}}^T{{\Upsilon }_{uf}}^{-1}{\stackrel{\sim }{\phi }}_{uf}\right\}}}\end{array}$

34

where $tr\{\cdot \}$ represents the trace of the corresponding matrix, and ${\Upsilon }_{lf}$ and ${\Upsilon }_{uf}$ are diagonal gain matrices to be designed. Let ${\eta }_1=1$ and take the time derivative of V ₁, then

$\begin{array}{l}{\stackrel{\sim }{\theta }}^T\dot{\stackrel{\sim }{\theta }}={\stackrel{\sim }{\theta }}^T[-K_1|\stackrel{\sim }{\theta }|\text{sgn}(\stackrel{\sim }{\theta })-K_2|\stackrel{\sim }{\theta }{|}^{{\eta }_2}\text{sgn}(\stackrel{\sim }{\theta })\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\epsilon }_{1}S-{\epsilon }_2|S{|}^{{\eta }_3}\text{sgn}(S)+\stackrel{\sim }{d}+R\stackrel{\sim }{\omega }+R\overline{\mu }]\\ \hspace{0.17em}\le -{\stackrel{\sim }{\theta }}^T{K}_1\stackrel{\sim }{\theta }-K_2||{\stackrel{\sim }{\theta }||}^{\frac{{\eta }_2+1}{\text{2}}}-{\stackrel{\sim }{\theta }}^T{\epsilon }_{1}S\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\stackrel{\sim }{\theta }}^T{\epsilon }_2|S{|}^{{\eta }_3}\text{sgn}(S)+{\stackrel{\sim }{\theta }}^T\stackrel{\sim }{d}+{\stackrel{\sim }{\theta }}^{T}R\stackrel{\sim }{\omega }+{\stackrel{\sim }{\theta }}^{T}R\overline{\mu }\end{array}$

35

We can denote that ${\text{lim}}_{t\to \infty }{\Vert S\Vert }^{\text{2}}\le \varpi$ with $\varpi =\text{max}\left\{\begin{array}{cc}\left(\frac{\kappa }{{\epsilon }_1{\lambda }_{\text{min}}(\Lambda )}\right),& {\left(\frac{\kappa }{{\epsilon }_2{\lambda }_{\text{min}}(\Lambda )}\right)}^{\text{2}/({\eta }_3+1)}\end{array}\right\}$ because it is known that S can converge to a bounded range according to Theorem 1. Then, equation (35) can be rewritten as follows

$\begin{array}{l}{\stackrel{\sim }{\theta }}^T\dot{\stackrel{\sim }{\theta }}\le -{\stackrel{\sim }{\theta }}^T(K_1-\frac{3}{2}I-\frac{1}{2}R{R}^T)\stackrel{\sim }{\theta }-K_2||\stackrel{\sim }{\theta }|{|}^{\frac{{\eta }_2+1}{\text{2}}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{2}{\stackrel{\sim }{d}}^T\stackrel{\sim }{d}+\frac{1}{2}{\overline{\mu }}^T\overline{\mu }+\frac{1}{2}{{\epsilon }_1}^2\varpi +\frac{3}{2}{{\epsilon }_2}^2{\varpi }^{{\eta }_3}+{\stackrel{\sim }{\theta }}^{T}R\stackrel{\sim }{\omega }\end{array}$

36

With reference to the work of Qian et al., ³⁴ the following inequality holds for the filter error item

${\overline{\mu }}^T\dot{\overline{\mu }}\le -\frac{1}{\tau }{\overline{\mu }}^T\overline{\mu }+\frac{1}{2}{\overline{\mu }}^T\overline{\mu }+\frac{1}{2}{\nu }^T\nu$

37

where $\nu (\stackrel{\sim }{\theta },\dot{\stackrel{\sim }{\theta }},\stackrel{\sim }{\omega },{\theta }_d,{\dot{\theta }}_d,{\ddot{\theta }}_d)$ is a continuous and bounded function. By substituting equations (36) and (37) into the derivative of V ₁, we can get

$\begin{array}{l}{\dot{V}}_1\le -{\stackrel{\sim }{\theta }}^T(K_1-\frac{3}{2}I-\frac{1}{2}R{R}^T)\stackrel{\sim }{\theta }-K_2||\stackrel{\sim }{\theta }|{|}^{\frac{{\eta }_2+1}{\text{2}}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-(\frac{1}{\tau }-1){\overline{\mu }}^T\mu +\frac{1}{2}({{\epsilon }_1}^2+1)\varpi +\frac{3}{2}{{\epsilon }_2}^2{\varpi }^{{\eta }_3}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\stackrel{\sim }{\theta }}^{T}R\stackrel{\sim }{\omega }+\frac{1}{2}{\nu }^T\nu +{\stackrel{\sim }{d}}^T\stackrel{\sim }{d}+\frac{\kappa }{{\lambda }_{\text{min}}({\Lambda }_1)}\end{array}$

38

Similar to the derivative of V ₁, ${\dot{V}}_2$ is calculated as follows

$\begin{array}{l}{\dot{V}}_2={\stackrel{\sim }{\omega }}^T\dot{\stackrel{\sim }{\omega }}+tr\{{{\stackrel{\sim }{\phi }}_{lf}}^T{{\Upsilon }_{lf}}^{-1}{\dot{\stackrel{\sim }{\phi }}}_{lf}\}+tr\{{{\stackrel{\sim }{\phi }}_{uf}}^T{{\Upsilon }_{uf}}^{-1}{\dot{\stackrel{\sim }{\phi }}}_{uf}\}\\ ={\stackrel{\sim }{\omega }}^T(-{\Lambda }^{*}\stackrel{\sim }{\omega }-R^T\theta +{{\stackrel{\sim }{\phi }}_{lf}}^T{\xi }_{lf}^{*}+{{\stackrel{\sim }{\phi }}_{uf}}^T{\xi }_{uf}^{*}+{\epsilon }_f)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+tr\{{{\stackrel{\sim }{\phi }}_{lf}}^T{{\Upsilon }_{lf}}^{-1}({\dot{\phi }}_{lf}^{*}-{\dot{\widehat{\phi }}}_{lf})\}+tr\{{{\stackrel{\sim }{\phi }}_{uf}}^T{{\Upsilon }_{uf}}^{-1}({\dot{\phi }}_{uf}^{*}-{\dot{\widehat{\phi }}}_{uf})\}\\ =-{\stackrel{\sim }{\omega }}^T{\Lambda }^{\text{*}}\stackrel{\sim }{\omega }-{\stackrel{\sim }{\omega }}^T{R}^T\stackrel{\sim }{\theta }+{\stackrel{\sim }{\omega }}^T{{\stackrel{\sim }{\phi }}_{lf}}^T{\xi }_{lf}^{*}+{\stackrel{\sim }{\omega }}^T{{\stackrel{\sim }{\phi }}_{uf}}^T{\xi }_{uf}^{*}+{\stackrel{\sim }{\omega }}^T{\epsilon }_f\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-tr\left\{{{\stackrel{\sim }{\phi }}_{lf}}^T{{\Upsilon }_{lf}}^{-1}{\dot{\widehat{\phi }}}_{lf}\right\}-tr\left\{{{\stackrel{\sim }{\phi }}_{uf}}^T{{\Upsilon }_{uf}}^{-1}{\dot{\widehat{\phi }}}_{uf}\right\}\end{array}$

39

To eliminate the terms ${\stackrel{\sim }{\omega }}^T{{\stackrel{\sim }{\phi }}_{lf}}^T{\xi }_{lf}^{*}$ and ${\stackrel{\sim }{\omega }}^T{{\stackrel{\sim }{\phi }}_{uf}}^T{\xi }_{uf}^{*}$ in equation (39), the adaptive updating laws of the parameter matrices in fuzzy estimation part can be designed as follows

$\begin{array}{l}{\dot{\widehat{\phi }}}_{lf}={\Upsilon }_{lf}{\xi }_{lf}^{*}{\stackrel{\sim }{\omega }}^T\\ {\dot{\widehat{\phi }}}_{uf}={\Upsilon }_{uf}{\xi }_{uf}^{*}{\stackrel{\sim }{\omega }}^T\end{array}$

40

and ${\Upsilon }_{lf}=\text{diag}\{{\gamma }_{lf1},{\gamma }_{lf2},{\gamma }_{lf3}\}({\gamma }_{lfi}>0, i\text{=1,2,3})$ and ${\Upsilon }_{uf}=\text{diag{}{\gamma }_{uf1},{\gamma }_{uf2},{\gamma }_{uf3}\}({\gamma }_{ufi}>0, i\text{=1,2,3})$ .

Then, ${\dot{V}}_2$ can be transformed into the following form

$\begin{array}{l}{\dot{V}}_2=-{\stackrel{\sim }{\omega }}^T{\Lambda }^{\text{*}}\stackrel{\sim }{\omega }-{\stackrel{\sim }{\omega }}^T{R}^T\theta +{\stackrel{\sim }{\omega }}^T{\epsilon }_f+{\stackrel{\sim }{\omega }}^T{{\stackrel{\sim }{\phi }}_{lf}}^T{\xi }_{lf}^{*}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\stackrel{\sim }{\omega }}^T{{\stackrel{\sim }{\phi }}_{uf}}^T{\xi }_{uf}^{*}-tr\left\{{{\stackrel{\sim }{\phi }}_{lf}}^T{\xi }_{lf}^{*}{\stackrel{\sim }{\omega }}^T\right\}-tr\left\{{{\stackrel{\sim }{\phi }}_{uf}}^T{\xi }_{uf}^{*}{\stackrel{\sim }{\omega }}^T\right\}\\ =-{\stackrel{\sim }{\omega }}^T{\Lambda }^{\text{*}}\stackrel{\sim }{\omega }-{\stackrel{\sim }{\omega }}^T{R}^T\theta +{\stackrel{\sim }{\omega }}^T{\epsilon }_f\\ \le -{\stackrel{\sim }{\omega }}^T{\Lambda }^{\text{*}}\stackrel{\sim }{\omega }-{\stackrel{\sim }{\omega }}^T{R}^T\theta +\frac{1}{2}{\stackrel{\sim }{\omega }}^T\stackrel{\sim }{\omega }+{\Vert {\epsilon }_f\Vert }^2\\ =-{\stackrel{\sim }{\omega }}^T\left({\Lambda }^{\text{*}}-\frac{1}{2}I\right)\stackrel{\sim }{\omega }-{\stackrel{\sim }{\omega }}^T{R}^T\theta +{\Vert {\epsilon }_f\Vert }^2\end{array}$

41

By combining equations (38) and (41), the derivative of V can be expressed as follows

$\begin{array}{l}\dot{V}={\dot{V}}_1+{\dot{V}}_2\\ \le -{\stackrel{\sim }{\theta }}^T(K_1-\frac{3}{2}I-\frac{1}{2}R{R}^T)\stackrel{\sim }{\theta }-K_2||\stackrel{\sim }{\theta }|{|}^{\frac{{\eta }_2+1}{\text{2}}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-(\frac{1}{\tau }-1){\overline{\mu }}^T\overline{\mu }+\frac{1}{2}({\epsilon }_1{ }^2+1)\varpi +\frac{3}{2}{\epsilon }_2{ }^2{\varpi }^{{\eta }_3}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\stackrel{\sim }{\theta }}^{T}R\stackrel{\sim }{\omega }+\frac{1}{2}{\nu }^T\nu +{\stackrel{\sim }{d}}^T\stackrel{\sim }{d}+\frac{\kappa }{{\lambda }_{\text{min}}({\Lambda }_1)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\stackrel{\sim }{\omega }}^T\left({\Lambda }^{\text{*}}-\frac{1}{2}I\right)\stackrel{\sim }{\omega }-{\stackrel{\sim }{\omega }}^T{R}^T\theta +{\Vert {\epsilon }_f\Vert }^2\\ =-{\stackrel{\sim }{\theta }}^T(K_1-\frac{3}{2}I-\frac{1}{2}R{R}^T)\stackrel{\sim }{\theta }\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-(\frac{1}{\tau }-1){\overline{\mu }}^T\overline{\mu }-{\stackrel{\sim }{\omega }}^T\left({\Lambda }^{\text{*}}-\frac{1}{2}I\right)\stackrel{\sim }{\omega }+{\Vert {\epsilon }_f\Vert }^2\end{array}$

42

Therefore, according to Qian et al., ³⁴ appropriate parameters can be selected to make the attitude system uniformly bounded.

Control allocation

In the above section, we have obtained the desired control torque to make the HFV system stable even if an actuator fault occurs through the FTC. In this section, the desired torque needs to be distributed into 8 aerodynamic surfaces and 10 RCS jets. The main goal is to minimize the error between the desired torque and the torque provided by aerodynamic surfaces to reduce the use of RCS jets which consume more fuel in the actual flight.

The desired control torque T is calculated by equation (32). The control allocation problem is addressed to minimize the index J_T by applying CSA

$\text{min }{J}_T=\Vert Q\left(T-B{\delta }_e\right)\Vert$

43

where J_T denotes the error index of the control torque and Q is the weight matrix. The admissible deflection angles ${\delta }_e$ are given in Table 1.

Table 1.

Admissible ranges of aerosurfaces.

Aerosurfaces	Minimum (rad)	Maximum (rad)
Inboard elevons	−0.44	0.44
Body flaps	−0.26	0.45
Rudders	−0.52	1.04
Outboard elevons	−0.52	0.52

The nest location of the new generation is updated via the following two approaches:

Global random walk via Lévy flight

$r_j^{d+1}=r_j^d+{\alpha }_j^{*}\otimes L(\lambda )$

44

where $r_j^d$ represents the $j\text{th}$ cuckoo in the population of $d\text{th}$ generation; ${\alpha }_j^{*}$ is the amplification factor; $\otimes$ is the Hadamard multiplication; and $L(\lambda )$ denotes the Lévy distribution, which represents a kind of random walk of animals.

Local random walk via the detection probability p_a

$r_j^{d+1}=r_j^d+{\beta }_j^{*}\cdot \left(r_m^d-r_n^d\right)\otimes H\left(p_a-h\right)$

45

where ${\beta }_j^{*}$ is the amplification factor; $r_m^d$ and $r_n^d$ are two different randomly chosen solutions (not equal to $r_j^{d+1}$ ). $h\sim U(0,1)$ is a random number; and $H(\cdot )$ denotes the Heaviside function.

From the above description, the detailed steps of control allocation based on CSA can be shown as follows:

Step 1. Initialize the search parameters including the number of initial allocation schemes (i.e. the number of the nests), the maximum number of iterations $d_{\text{max}}$ , the given accuracy of the control torque error ${\epsilon }_J$ (i.e. the fitness of nest position), and the detection probability p_a .

Step 2. Randomly select N sets of initial deflection angles of aerosurfaces $r_j^0(j=1,2,\cdots ,N)$ , which represents the nest position of the first generation. It is worth noting that $r_j^0(j=1,2,\cdots ,N)$ should be within the given ranges as presented in Table 1 as much as possible. Then calculate the fitness of each set of deflection angles according to equation (43).

Step 3. Update the deflection angles of eight aerosurfaces via equation (44) and calculate the fitness of the new set of deflection angles. If the new fitness is better, keep the update; otherwise, give up the update.

Step 4. Randomly select a number $h\sim U(0,1)$ . If $h>p_a$ , jump to Step 5. Otherwise, update the deflection angles of eight aerosurfaces via equation (45) and calculate the fitness of the new set of deflection angles. If the new fitness is better, keep the update; otherwise, give up the update.

Step 5. Calculate the current optimal solution. If current error index J_T is less than the given accuracy ${\epsilon }_J$ or the number of iterations exceeds $d_{\text{max}}$ , end the search process; otherwise, return back to Step 2.

Remark 2

Similar to other metaheuristic algorithms, the process of the CSA is of randomicity. Compared with the allocation method proposed by Gao et al. ¹⁰ , it can improve the efficiency of the aerosurfaces and reduce the consumption of RCS fuel.

Numerical simulation

In this section, simulation results of the proposed methods are presented and discussed. Different cases, including interval type-2 fuzzy estimation, FTC, and control allocation, are considered.

The aerodynamic parameters of X-33 HFV can be obtained from the work of Zhao et al. ³⁵ The desired attitude angle vector ${\theta }_d$ is set as ${[\text{0}.8°,\text{ 2}°,\text{ 0}°]}^T$ , while the initial attitude angle vector is ${[\text{0}°,\text{ 0}°,\text{ 0}°]}^T$ . The moment of inertia tensor is given by

$J=\left[\begin{array}{ccc}554486& 0& -23002\\ 0& 1136949& 0\\ -23002& 0& 1376853\end{array}\right]\text{ kg}\cdot {\text{m}}^2$

46

Case A: interval type-2 fuzzy estimation

The right body flap is assumed to have an additive fault, which can be expressed as follows

$\{\begin{array}{l}{f}_{{\delta }_3}=0,t<15s\hfill \\ f_{{\delta }_3}=20, t\ge 15s\hfill \end{array}$

47

The nonlinear function caused by the modeling error is assumed as the following form

$g(\omega (t))=J^{-1}\left[\begin{array}{c}1.8\times {10}^4\text{sin}(t/5)+2\times {10}^4(p^2+\text{exp}(q))+1.5\times {10}^{4}p\\ 5\times {10}^4+5\times {10}^{3}q\\ 2\times {10}^4+5\times {10}^4\text{cos}((t-20)/4)-{10}^{4}r\end{array}\right]$

48

The upper and lower membership functions of the fuzzy estimator are given by

$\{\begin{array}{c}\begin{array}{l}{\overline{\mu }}_{F_p}^1(p)=\text{exp}[-80\times {(p-0.11)}^2]\\ {\underline{\mu }}_{F_p}^1(p)=\text{exp}[-80\times {(p-0.09)}^2]\end{array}\\ \begin{array}{l}{\overline{\mu }}_{F_p}^2(p)=\text{exp}[-80\times {(p-0.03)}^2]\\ {\underline{\mu }}_{F_p}^2(p)=\text{exp}[-80\times {(p-0.02)}^2]\end{array}\\ \begin{array}{l}\hspace{0.17em}{\overline{\mu }}_{F_p}^3(p)=\text{exp}[-80\times {(p+0.09)}^2] \\ \hspace{0.17em}{\underline{\mu }}_{F_p}^3(p)=\text{exp}[-80\times {(p+0.11)}^2]\end{array}\end{array}$

49

$\{\hspace{0.17em}\begin{array}{c}\begin{array}{l}{\overline{\mu }}_{F_q}^1(q)=\text{exp}[-80\times {(q-0.11)}^2]\\ {\underline{\mu }}_{F_q}^1(q)=\text{exp}[-80\times {(q-0.09)}^2]\end{array}\\ \begin{array}{l}{\overline{\mu }}_{F_q}^2(q)=\text{exp}[-80\times {(q-0.03)}^2]\\ {\underline{\mu }}_{F_q}^2(q)=\text{exp}[-80\times {(q-0.02)}^2]\end{array}\\ \begin{array}{l}{\overline{\mu }}_{F_q}^3(q)=\text{exp}[-80\times {(q+0.09)}^2]\\ {\underline{\mu }}_{F_q}^3(q)=\text{exp}[-80\times {(q+0.11)}^2]\end{array}\end{array}$

50

$\{\begin{array}{c}\begin{array}{l}{\overline{\mu }}_{F_r}^1(r)=\text{exp}[-80\times {(r-0.11)}^2]\\ {\underline{\mu }}_{F_r}^1(r)=\text{exp}[-80\times {(r-0.09)}^2]\end{array}\\ \begin{array}{l}{\overline{\mu }}_{F_r}^2(r)=\text{exp}[-80\times {(r-0.03)}^2]\\ {\underline{\mu }}_{F_r}^2(r)=\text{exp}[-80\times {(r-0.02)}^2]\end{array}\\ \begin{array}{l}\hspace{0.17em}{\overline{\mu }}_{F_r}^3(r)=\text{exp}[-80\times {(r+0.09)}^2] \\ \hspace{0.17em}{\underline{\mu }}_{F_r}^3(r)=\text{exp}[-80\times {(r+0.11)}^2]\end{array}\end{array}$

51

The adaptive updating parameters in the interval type-2 fuzzy estimator are set as ${\gamma }_{lf}=\text{diag}\{200,200,200\}$ and ${\gamma }_{uf}=\text{diag}\{100,100,100\}$ . Then, the estimation result of the unknown nonlinear function $G_f(x)$ is shown in Figure 2.

Figure 2.

Estimation results of G_f (x).

From Figure 2, it can be seen that the estimation curve can track $G_f(x)$ after approximately 5 s with accuracy.

Thus, the type-2 fuzzy estimator can estimate the lumped effect of the nonlinear function caused by modeling error and the right body flap fault accurately and quickly.

Case B: FTC

The parameters of the terminal sliding mode controller are selected as $k_{1i}=0.2$ , $k_{2i}=0.5$ , ${\epsilon }_1=3, {\epsilon }_2=2$ , and ${\Lambda }_1=\text{diag}\left\{10,10,10\right\}$ . And as for the dynamic surface controller, we select ${\Lambda }^{*}=\text{diag}\{100,100,100\}$ .

In this case, the control torque T is assumed to be sufficient and control allocation is not considered. The response curves of attitude angles and angular rates under different control methods are presented in Figures 3 and 4, respectively.

Figure 3.

Angular rate response under different controllers.

Figure 4.

Attitude angle responses under different controllers.

In Figure 3, it is shown that three angular rates can remain stable under nominal control and FTC when the right body flap fault occurs. However, the actual attitude angles deviate from the desired ones under the nominal controller, whereas three attitude angles track the reference command accurately under self-healing controller, as is shown in Figure 4. Thus, the FTC proposed in this study can guarantee the stability and the tracking accuracy of HFV. Here, the nominal controller is designed based on the HFV system without faults. That is to say, the nominal controller means the control method does not consider the fault compensation signals.

Case C: CSA

In this case, the deflection angles of 8 aerosurfaces and the on/off states of 10 RCS jets under different allocation algorithms when the actuator fault exists are presented in Figures 5 and 6, respectively. Additionally, three attitude angle curves under different controllers considering control allocation are given in Figure 7.

Figure 5.

The deflection angle signals of aerosurfaces with different allocation algorithms.

Figure 6.

RCS on/off states with different allocation methods. RCS: reaction control system.

Figure 7.

Attitude angular responses with different control schemes.

Figure 5 shows the utilization efficiency of aerosurfaces under CSA is higher than under the pseudoinverse allocation algorithm proposed by Gao et al. ¹⁰ . When the right body flap fault occurs, the left and right rudders almost have no compensation effect under the algorithm in the work of Gao et al. ¹⁰ , whereas CSA distributes the control torque to more aerosurfaces. Figure 6 shows RCS jets do not have too much switching states under CSA, which is beneficial to the operation of the system. In addition, the running time of the RCS system under the pseudoinverse method is longer than under CSA.

In Figure 7, it can be seen that with the interval type-2 fuzzy estimator and CSA, the stability of the system can be ensured and three attitude angles can track the desired ones accurately and rapidly. The transient performance of the response curves with the self-healing controller proposed in this study is better than that with the controller in the work of Gao et al. ¹⁰

Conclusions

In this article, a fault estimation and FTC scheme is developed for HFV subject to body flap additive fault and modeling error. With an interval type-2 fuzzy estimator, the lumped effect of the body flap fault and modeling error can be effectively approximated at first. The virtual control input of angular rate loop is obtained by constructing a terminal sliding mode controller, which avoids the imaginary situation successfully. Then, the dynamic surface technique is utilized to calculate the desired control torque based on the fuzzy estimation result, thereby preventing trajectory derivation of attitude angles. Finally, the control allocation problem is solved by combining the CSA and linear programming methods. With this allocation scheme, the utilization efficiency of 10 aerosurfaces is effectively improved. Several comparative simulation results for different cases illustrate the effectiveness of the presented fault estimation and compensation scheme.