Attitude control precision analysis for hypersonic vehicles based on covariance analysis describing equation technique

Abstract

This article investigates a control precision analysis based on covariance analysis describing equation technique for the generic hypersonic vehicle attitude tracking system. For a nominal generic hypersonic vehicle attitude system, a global sliding mode control strategy is pointed out first and then the stabilization of the control strategy is analyzed based on Lyapunov stability theory. In the presence of external disturbance moment and internal parameter perturbation, the statistical attitude model with the proposed control scheme is obtained using the covariance analysis describing equation technique. Moreover, 500 Monte Carlo tests are simulated for the hypersonic vehicle with the same disturbance and those mean and standard deviation are same as the results of the one-time covariance analysis describing equation technique simulation. Finally, the comparison of the covariance analysis describing equation technique and the Monte Carlo tests results shows that the covariance analysis describing equation technique is more effective and timesaving.

Introduction

Recently, the hypersonic vehicle is a more reliable and cost-efficient way to access space. ¹ This kind of vehicle has complicated environmental characteristics, large flight envelope, high flying Mach number, and high maneuverability. Its special designed structure of the propulsion airframe integration causes heavy coupling between the propulsive and aerodynamic forces, which makes the control of the generic hypersonic vehicle (GHV) a difficulty. In reentry phase, due to the large flying altitude changes, the dynamic pressure and the Mach number of the aircraft are drastically changed. ² Thus, the aerodynamic characteristics of the aircraft are difficult to be measured and estimated. In addition, the rapid change of flying altitude also causes the unknown interference by the thermal effects. These characteristics make the traditional linear control methods no longer applicable. In order to improve the GHV attitude control performance, various kinds of nonlinear control methods have been proposed such as sliding mode control (SMC),^3–5 backstepping method,^6,7 adaptive control,^8,9 and disturbance-based control.^10–12

Backstepping control is an effective tool for the design of controllers of nonlinear system. The idea of backstepping is to design a virtual controller using some state variables and designing for them intermediate control laws. ¹ It shows unique advantages in dealing with nonlinear uncertain systems, whether the uncertainties are matched or mismatched. The stability of the backstepping method can be obtained by proving each virtual controllers’ stability with selected Lyapunov functions. However, a disadvantage of the backstepping method is that the derivative of the virtual controller is usually needed in the next step controller design. Thus, a command filter is usually used. ¹³ For the attitude model of the GHV, the backstepping method with different controllers such as SMC or adaptive controller is widely used in recent years. In Liu et al., ¹⁴ based on the backstepping method, a new, fast, and singularity-avoided terminal sliding mode controller was proposed. In Wang et al., ¹⁵ an adaptive robust attitude tracking control scheme for variable structure near space vehicle was presented. The adaptive fuzzy systems were employed for estimating unknown functions in the dynamic model, and those parameters were updated online. In Zhang et al., ¹⁶ based on the backstepping method, an adaptive sliding mode controller combined with radial base function neural network was used to solve attitude tracking problem.

For all the above references, most of their simulation verifications only simulated once in the case of external disturbance and internal parameters perturbation. More attention is paid to design a better robust and effective control method. Tian et al. ¹⁰ delivered an adaptive multivariable disturbance compensation control strategy based a finite-time convergence disturbance observer for the attitude tracking problem of the reentry vehicle where the bounds of the perturbations are not known. In Han, ¹⁷ the active disturbance rejection control (ADRC) technique was introduced. An account of each component of ADRC as well as its structure and philosophy was given. Shao and Wang ¹¹ presented a control scheme combined trajectory linearization control (TLC) and ADRC for reentry vehicle attitude tracking system considering bounded uncertainties. They compared three controllers TLC, TLC-ADRC, and an adaptive robust controller with neutral networks in their simulation to demonstrate the superiority of the TLC-ADRC. Cheng et al. ¹² delivered a predictive sliding mode control (PSMC) method with nonlinear disturbance observer (NDO) for attitude control of hypersonic vehicle. Different control schemes like SMC, PSMC, and NDO-PSMC are compared in their simulation. Considering the variation of the aerodynamic coefficients, they also carried out 100 times Monte Carlo tests via SMC and PSMC. However, the number of their simulation times is not enough, and they did not take into account the disturbance moment during the Monte tests. Therefore, to improve the degree of confidence of the control methods, plenty of simulations in different situations are necessary. That is, hundreds of Monte Carlo tests are needed to be carried out. However, numerous simulation trials will cost a lot of time and the efficiency is not high. Therefore, a faster statistical analysis method is required.

As we mentioned above, the Monte Carlo method, as a traditional precision analysis method, needs a large number of computer trials with different disturbance and initial states to obtain the target system statistical properties. In other words, the Monte Carlo method is not suitable for situations where rapid analysis is required. The covariance analysis describing equation technique (CADET), developed at The Analytic Science Corporation in the 1970s, is a direct method to analyze the statistical properties of the nonlinear system with the random inputs and can save more computer proceeding time compared with the traditional approach, Monte Carlo method. Moreover, the CADET method has been proven to be a powerful mean for directly evaluating the statistical behavior of a guidance system. The pivotal early work of the CADET is from Gelb and Warren. ¹⁸ Taylor ¹⁹ delivered a guiding book that presented detailed instructions for the application of the CADET in a tactical missile guidance system. This book is the basis of the CADET as well as this article. In the recent years, the CADET has been applied in different problems. Wang and Xie ²⁰ presented the application of the CADET in the assessment of missile hit probability in combination with a nonlinear model of missile. Jiang et al. ²¹ developed a modified CADET method to solve the problem of analyzing a nonlinear system with parametric perturbation. In Zhang et al., ²² considering the line of sight (LOS) rate accuracy, acceleration saturation capabilities of interceptor, and sampling rate, Zhang and Li used CADET analyze the endgame performance of an interceptor.

Inspired by the above research, this article investigates a control precision analysis for the hypersonic vehicle attitude tracking system based on CADET. To the best of authors’ knowledge, the CADET is usually used in analyzing the precision of the guidance system. There is no application in attitude control accuracy analysis in article, especially for strong nonlinear, coupling systems such as hypersonic vehicles. The main contribution of this article is that in the presence of external disturbance and parameter perturbation the statistical characteristics of the attitude system are obtained using the CADET. To present the effectiveness of CADET, 500 Monte Carlo tests are simulated under the different kinds of noise. The disturbance and parameter perturbation are regarded as Gaussian noises in following simulations. A sliding mode backstepping control method is proposed in this article, and a low-pass filter is used to obtain the derivative of the first-loop virtual control. The contributions of this article are listed in the following:

1. In this article, the CADET is used to analyze the statistical characteristics for the GHV attitude tracking control system. In order to get a good control accuracy and robustness, a global sliding mode backstepping controller is used to solve the problem of attitude control of GHV.

2. Considering the external disturbance and parameter perturbation, CADET is used to linearize the controlled attitude system, and the statistical model of the system is obtained.

3. Compared with the 500 Monte Carlo simulation results, the results of CADET has the same variation trend and similar accuracy with less simulation time, under the action of external disturbances and parameter perturbation.

The rest of this article is organized as follows. Section “The attitude model of GHV” describes the dynamic attitude model of the GHV. The control strategy is presented in section “Sliding mode controller design of the GHV.” In section “Application of CADET in attitude control of the GHV,” the methodology of CADET and linearized model of the controlled system are derived. Simulation based on MATLAB/Simulink platform is included in section “Simulation.” Finally, this article is completed with conclusion.

The attitude model of GHV

In this article, the attitude model of the GHV under consideration is taken from. ¹⁵ The kinematics and dynamics equations are shown as follows

$$\begin{array}{l}\begin{array}{l}\stackrel{.}{\mathbf{\Omega }}={\mathit{f}}_s+{\mathit{g}}_s\omega +{\mathit{d}}_s\hfill \\ \stackrel{.}{\mathit{\omega }}={\mathit{f}}_f+{\mathit{g}}_{f}M+{\mathit{d}}_f\hfill \end{array}\hfill \end{array}$$

(1)

$\mathbf{\Omega }=[\alpha ,\beta ,\mu ]$ represents the attitude vector which contains the angle of attack $\alpha$ , the sideslip angle $\beta$ , and the bank angle $\mu$ . The $\mathbf{\omega }=[p,q,r]$ is the angular rate vector; and the $p$ , $q$ , and $r$ represent the roll rate, pitch rate, and yaw rate, respectively. The ${\mathit{d}}_s$ and ${\mathit{d}}_f$ are the external disturbances. And

$${\mathit{f}}_s={[f_{\alpha }\mathrm{,}{f}_{\beta }\mathrm{,}{f}_{\mu }]}^T$$

(2)

where

$$\begin{array}{c}{f}_{\alpha }=\frac{1}{MV\cos \beta }(-qS{C}_{L\mathrm{,}\alpha }+Mg\cos \gamma \cos \mu -T\sin \alpha )\\ f_{\beta }=\frac{1}{MV}(qS{C}_{Y\mathrm{,}\beta }\beta \cos \beta +Mg\cos \gamma \sin \mu -T\sin \beta \cos \alpha )\\ f_{\mu }=-\frac{g}{V}\cos \gamma \cos \mu \tan \beta +\frac{1}{MV}qS{C}_{Y\mathrm{,}\beta }\beta \tan \gamma \cos \mu \cos \beta \\ +\frac{1}{MV}qS{C}_{L\mathrm{,}\alpha }(\tan \gamma \sin \mu +\tan \beta )\\ +\frac{T}{MV}[\sin \alpha (\tan \gamma \sin \mu +\tan \beta )-\cos \alpha \tan \gamma \cos \mu \sin \beta ]\\ {\mathit{f}}_f={[f_p\mathrm{,}{f}_q\mathrm{,}{f}_r]}^T\end{array}$$

(3)

where

$$\begin{array}{l}\begin{array}{l}{f}_p=I_{xx}^{-1}\left[qr(I_{yy}-I_{zz})-{\stackrel{.}{I}}_{xx}p\right]\hfill \\ f_q=I_{yy}^{-1}\left[pr(I_{zz}-I_{xx})-{\stackrel{.}{I}}_{yy}q\right]\hfill \\ f_r=I_{zz}^{-1}\left[pq(I_{yy}-I_{zz})-{\stackrel{.}{I}}_{zz}r\right]\hfill \end{array}\\ {\mathit{g}}_s=\left[\begin{array}{ccc}-\tan \beta \cos \alpha & 1& -\tan \beta \sin \alpha \\ \sin \alpha & 0& -\cos \alpha \\ \sec \beta \cos \alpha & 0& \sec \beta \sin \alpha \end{array}\right]\end{array}$$

(4)

$${\mathit{g}}_f=diag(I_{xx}^{-1}\mathrm{,}{I}_{yy}^{-1}\mathrm{,}{I}_{zz}^{-1})$$

(5)

where $M$ and $V$ are the mass and velocity of the GHV, respectively; $q=(1/2)\rho V^2$ is the dynamic pressure and $g$ is the gravitation acceleration. $S$ is the wing area and $\gamma$ is the flight-path angle. $\mathit{M}$ is the control moment vector which contains the moments generated by the reaction control system (RCS) and the aerodynamic surfaces. $C_{L,\alpha }$ and $C_{Y,\beta }$ are the lift coefficients and side force coefficients, respectively. In this article, the coefficients we used are from S Keshmiri et al.^23,24 Based on the wind tunnel data and computational fluid dynamics (CFD) studies in Shaughnessy et al., ²⁵ they developed an aerodynamic database for the blunt body of the GHV. In their work, they provided both force and moment coefficients in different speeds of sound. In this article, the system satisfies the following assumption.

Assumption 1

For simplifying the attitude model of the GHV, we assume that there is no thrust and the moment of inertia do not change, which means the thrust $\mathit{T}=0$ and ${\stackrel{·}{I}}_{\mathrm{xx}}={\stackrel{·}{I}}_{\mathrm{yy}}={\stackrel{·}{I}}_{\mathrm{zz}}=0$ .

Sliding mode controller design of the GHV

In this section, a sliding mode controller is designed for the attitude model of the GHV. Here, the backstepping method is used. Based on the multi-time scale attitude model that the dynamic response rate of the inner loop is faster than the outer loop, the controller system can be divided into two parts. ¹⁴ Thereby, two SMC controllers are designed for the slowloop and the fastloop, respectively. Figure 1 shows the structure of the SMC control for the attitude model of the GHV. We use $\Vert ·\Vert$ for the standard Euclidean norm of the vector. The following lemma will be used.

Figure 1.

The SMC control structure of the GHV.

Lemma 1

If $V\ge 0$ , $\forall t\ge t_0\ge 0$ , the solution of the inequality $\stackrel{·}{V}\le -\alpha V+f$ is ²⁶

$$V(t)\le e^{-\alpha (t-t_0)}V(t_0)+{\displaystyle \underset{t_0}{\overset{t}{\int }}{e}^{-\alpha (t-\tau )}f(\tau )d\tau }$$

where $a$ is an arbitrary constant.

Slowloop controller

In this part, the desired angular rate vector ${\mathbf{\omega }}_c$ is obtained using a global sliding mode controller. The sliding surface is selected as

$${\mathit{s}}_1={\mathbf{\Omega }}_e-{\mathbf{\Omega }}_e(0)e^{-k_{e1}t}$$

(6)

where the ${\mathbf{\Omega }}_e={\mathbf{\Omega }}_c-\mathbf{\Omega }$ is the error of the attitude angle. A global sliding mode controller is utilized by adding a attenuation term $-{\mathbf{\Omega }}_e(0)e^{-k_{e1}t}$ . What’s more, the exponential reaching law is used

$${\stackrel{·}{\mathit{s}}}_1=-k_1{\mathit{s}}_1-{ϵ}_1\mathrm{sig}{n}_1({\mathit{s}}_1)$$

(7)

Substituting the derivative of equation (6) into equation (7), thus

$${\stackrel{·}{\mathbf{\Omega }}}_e+k_{e1}{\mathbf{\Omega }}_e(0)e^{-k_{e1}t}=-k_1{\mathit{s}}_1-{ϵ}_1\mathrm{sig}{n}_1({\mathit{s}}_1)$$

(8)

which yields to

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{w}}_c=}& {\mathit{g}}_s^{-1}\{k_1{\mathbf{\Omega }}_e+(k_{e1}-k_1)e^{-k_{e1}t}{\mathbf{\Omega }}_e(0)\hfill \\ \multicolumn{1}{c}{}& +{ϵ}_1\mathrm{sign}({\mathit{s}}_1)+{\stackrel{·}{\mathbf{\Omega }}}_c-{\mathit{f}}_s\hfill \end{array}\}$$

(9)

Fastloop controller

In this part, the control input ${\mathit{M}}_c$ is obtained. Similarly to the slowloop part, the same sliding mode controller is used. The sliding surface is selected as follows

$${\mathit{s}}_2={\mathit{w}}_e-{\mathit{w}}_e(0)e^{-k_{e2}t}$$

(10)

where the ${\mathit{w}}_e={\mathit{w}}_c-\mathit{w}$ . Thus

$${\stackrel{·}{\mathit{w}}}_e+k_{e2}{\mathit{w}}_e(0)e^{-k_{e2}t}=-k_2{\mathit{s}}_2-{ϵ}_2\mathrm{sig}{n}_2({\mathit{s}}_2)$$

(11)

Similarly to equation (8), the control input ${\mathit{M}}_c$ is

$${\mathit{M}}_c={\mathit{g}}_f^{-1}\left[k_2{\mathit{w}}_e+{ϵ}_{2}sign({\mathit{s}}_2)+(k_{e2}-k_2)e^{-k_{e2}t}{\mathit{w}}_e(0)+{\stackrel{.}{\mathit{w}}}_c-{\mathit{f}}_f\right]$$

(12)

Here, to obtain the ${\stackrel{·}{\mathit{w}}}_c$ , a low-pass filter is used ²⁷

$${\mathit{w}}_m=\frac{1}{\mathrm{Ts}+1}{\mathit{w}}_c$$

Then

$${\stackrel{·}{\mathit{w}}}_m=\frac{1}{T}{\mathit{w}}_c-\frac{1}{T}{\mathit{w}}_m$$

(13)

where the $T$ is the time constant.

Stability analysis

The stability of the above sliding mode controllers is analyzed in this part. The following theorem holds.

Theorem 1

Consider the system (1) with the proposed SMC controller (9) and (12), and Assumption 1 holds. If $k_i>0$ , ${ϵ}_i>0$ , and $k_{e_i}>0$ , $i=1,2$ , the system is stable and the tracking error ${\Omega }_e$ and ${\omega }_e$ can asymptotically converge to zero.

Proof

The Lyapunov function is selected as

$$\begin{array}{cc}\multicolumn{1}{c}{V}& =V_1+V_2+V_3\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{2}{\mathit{s}}^T\mathit{s}+\frac{1}{2}{\mathbf{\Omega }}_e^T{\mathbf{\Omega }}_e+\frac{1}{2}{\mathbf{\omega }}_e^T{\mathbf{\omega }}_e+\frac{1}{2}{\mathit{e}}^T\mathit{e}\hfill \end{array}$$

where $V_1=(1/2){\mathit{s}}^T\mathit{s}$ , $V_2=(1/2){\mathbf{\Omega }}_e^T{\mathbf{\Omega }}_e$ , $V_3=(1/2){\mathbf{\omega }}_e^T{\mathbf{\omega }}_e+(1/2){\mathit{e}}^T\mathit{e}$ , and $\mathit{e}={\mathbf{\omega }}_m-{\mathbf{\omega }}_c$ is the error caused using the low-pass filter.

For $V_1$ , on the derivative

$${\stackrel{·}{V}}_1={\mathit{s}}^T\stackrel{·}{\mathit{s}}$$

Substituting equation (7) or equation (11), thus

$${\stackrel{·}{V}}_1={\mathit{s}}^T[-k\mathit{s}-ϵ\mathrm{sign}(s)]=-k\Vert \mathit{s}{\Vert }^2-ϵ\Vert \mathit{s}\Vert <0$$

(14)

For $V_2$

$${\stackrel{·}{V}}_2={\mathbf{\Omega }}_e^T{\stackrel{·}{\mathbf{\Omega }}}_e$$

Considering equation (8), we have

$${\stackrel{·}{V}}_2={\mathbf{\Omega }}_e^T(-k_{e_1}{\mathbf{\Omega }}_e(0)e^{-k_{e_1}t}-k_1{\mathit{s}}_1-{ϵ}_1\mathrm{sign}({\mathit{s}}_1))$$

If $t\to \infty$ , $e^{-k_{e_1}t}\to 0$ , thus

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{V}}_2}& ={\mathbf{\Omega }}_e^T(-k_1{\mathbf{\Omega }}_e-{ϵ}_1\mathrm{sign}({\mathbf{\Omega }}_e))\hfill \\ \multicolumn{1}{c}{}& =-k_1\Vert {\mathbf{\Omega }}_e{\Vert }^2-{ϵ}_1\Vert {\mathbf{\Omega }}_e\Vert <0\hfill \end{array}$$

(15)

As for the $V_3$ , if $V_3(0)\le p$ and $p\ge 0$ , $\Vert {\stackrel{·}{\mathbf{\omega }}}_c{\Vert }^2$ is bounded and $\Vert {\stackrel{·}{\mathbf{\omega }}}_c{\Vert }^2\le M^2$ . ²⁸ Thus, $(\Vert {\stackrel{·}{\mathbf{\omega }}}_c{\Vert }^2/M^2)-1\le 0$ .

Considering the derivative of the $V_3$

$${\stackrel{·}{V}}_3={\mathbf{\omega }}_e^T{\stackrel{·}{\mathbf{\omega }}}_e+{\mathit{e}}^T\stackrel{·}{\mathit{e}}$$

With $t\to \infty$ , which yields to

$$\begin{array}{c}\stackrel{.}{V}={\mathit{\omega }}_e^T(-k_2{\mathit{\omega }}_e-{ϵ}_{2}sign({\mathit{\omega }}_e))+{\mathit{e}}^T{\left(-\frac{\mathit{e}}{T}-{\stackrel{.}{\mathit{\omega }}}_c\right)}_3\\ =-k_2{\Vert {\mathit{\omega }}_e\Vert }^2-{ϵ}_2\Vert {\mathit{\omega }}_e\Vert -\frac{1}{T}{\Vert \mathit{e}\Vert }^2-{\mathit{e}}^T{\stackrel{.}{\mathit{\omega }}}_c\\ \le -k_2{\Vert {\mathit{\omega }}_e\Vert }^2-{ϵ}_2\Vert {\mathit{\omega }}_e\Vert -\frac{1}{T}{\Vert \mathit{e}\Vert }^2+\left|{\mathit{e}}^T\right|\left|{\stackrel{.}{\mathit{\omega }}}_c\right|\\ \le -k_2{\Vert {\mathit{\omega }}_e\Vert }^2+\left(-\frac{1}{T}+{\Vert {\stackrel{.}{\mathit{\omega }}}_c\Vert }^2\right){\Vert \mathit{e}\Vert }^2+\frac{1}{4}\\ \le -k_2{\Vert {\mathit{\omega }}_e\Vert }^2+\left(-\frac{1}{T}+{\mathit{M}}^2\right){\Vert \mathit{e}\Vert }^2+\frac{1}{4}\end{array}$$

Choose $k_2\ge r$ and $(1/T)\ge M^2+r$ , where $r>0$ . Then

$$\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{V}}_3}& \le -r{\Vert {\mathbf{\omega }}_e\Vert }^2+({\Vert {\stackrel{·}{\mathbf{\omega }}}_c\Vert }^2-M^2-r){\Vert \mathit{e}\Vert }^2+\frac{1}{4}\hfill \\ \multicolumn{1}{c}{}& =-r{V}_3+M^2(\frac{{\Vert {\stackrel{·}{\mathbf{\omega }}}_c\Vert }^2}{M^2}-1){\Vert \mathit{e}\Vert }^2+\frac{1}{4}\hfill \\ \multicolumn{1}{c}{}& \le -r{V}_3+\frac{1}{4}\hfill \end{array}$$

(16)

Considering $V_3=p$ , thus ${\stackrel{·}{V}}_3\le -\mathrm{rp}+(1/4)$ . If we choose $r\ge (1/4p)$ , ${\stackrel{·}{V}}_3\le 0$ . Therefore, $V_3(t)\le p$ , if $r\ge (1/4p)$ . Meanwhile, with using Lemma 1, we obtain the solution of equation (16)

$$V_3(t)\le e^{-r(t-t_0)}{V}_3(t_0)+\frac{1}{4r}(1-e^{-r(t-t_0)})$$

With $t\to \infty$

$$\underset{t\to \infty }{lim}{V}_3(t)\le \frac{1}{4r}$$

Therefore, $V_3$ converges asymptotically, and the convergence latitude depends on $r$ . What’s more, considering $(1/T)\ge M^2+r$ , if $T\to 0$ , $r$ can be chosen to $\infty$ and the limit of the $V_3$ tends to 0. From equations (14), (15), and (16), when $t\to \infty$ , it is clear that

$$\stackrel{·}{V}={\stackrel{·}{V}}_1+{\stackrel{·}{V}}_2+{\stackrel{·}{V}}_3<0$$

(17)

So far, the stability of the system including the slowloop and the fastloop is proved using the Lyapunov function theory, and the tracking errors are able to converge to zero.

Application of CADET in attitude control of the GHV

Methodology of the CADET

The covariance analysis describing function technique, developed by The Analytic Science Corporation in 1970s, is a method for analyzing the statistical characters of a nonlinear system with random inputs. ¹⁹ In this part, the basic methodology is presented according to Taylor. ¹⁹ Normally, a nonlinear system can be presented as follows

$$\stackrel{·}{\mathit{x}}(t)=\mathit{f}(\mathit{x},t)+\mathit{G}(t)\mathit{w}(t)+\mathit{D}(t)$$

(18)

where $\mathit{x}(t)$ and $\mathit{w}(t)$ is the system state vector and system input vector as well as the external disturbance. The $\mathit{f}(\mathit{x},t)$ is the nonlinear system function and we assume that $\mathit{f}(\mathit{x},t)$ and $\mathit{G}(t)$ are all continuous with respect to t.

The inputs $\mathit{w}(t)$ is composed of a deterministic value $\mathit{b}(t)$ and a random component $\mathit{u}(t)$ , thus

$$\mathit{w}(t)=\mathit{b}(t)+\mathit{u}(t)$$

The $\mathit{w}(t)$ is assumed to be a white noise process, and its spectral density matrix is $\mathit{Q}$ , thereby statistical characters of the $\mathit{w}(t)$ are specified by

$$\begin{array}{c}E[\mathit{w}(t)]=\mathit{b}(t)\hfill \\ E[\mathit{w}(t)\mathit{w}{(t)}^T]=\mathit{Q}\delta (t)(t-\tau )\hfill \end{array}$$

Similarly, the state vector $\mathit{x}(t)$ has a deterministic value $\mathit{m}(t)$ and a random part $\mathit{r}(t)$ . Then, $\mathit{x}(t)$ is described as follows

$$\begin{array}{c}\mathit{x}(t)=\mathit{m}(t)+\mathit{r}(t)\hfill \\ \mathit{m}(t)=E[\mathit{x}(t)]\hfill \\ \mathit{P}(t)=E[\mathit{r}(t)\mathit{r}{(t)}^T]\hfill \end{array}$$

where $\mathit{P}$ is the covariance matrix. According to Taylor ¹⁹

$$\begin{array}{cc}\multicolumn{1}{c}{\stackrel{·}{\mathit{m}}}& =E[\stackrel{·}{\mathit{x}}(t)]=E[\mathit{f}(\mathit{x},t)+\mathit{G}(t)\mathit{w}(t)+\mathit{D}(t)]\hfill \\ \multicolumn{1}{c}{}& =E[\mathit{f}(\mathit{x},t)]+\mathit{G}(t)E[\mathit{w}(t)]+E[\mathit{D}(t)]\hfill \\ \multicolumn{1}{c}{}& =\hat{\mathit{f}}(\mathit{x},t)+\mathit{G}(t)\mathit{b}(t)+\mathit{D}(t)\hfill \end{array}$$

(19)

and

$$\begin{array}{cc}\multicolumn{1}{c}{\stackrel{·}{\mathit{P}}}& =\frac{d}{\mathrm{dt}}E[\mathit{r}(t){\mathit{r}}^T(t)]\hfill \\ \multicolumn{1}{c}{}& =E[\stackrel{·}{\mathit{r}}(t){\mathit{r}}^T(t)]+E[\mathit{r}(t){\stackrel{·}{\mathit{r}}}^T(t)]\hfill \\ \multicolumn{1}{c}{}& =E[\mathit{f}{\mathit{r}}^T(t)]+E[\mathit{r}(t){\mathit{f}}^T]+\mathit{G}(t)\mathit{Q}{\mathit{G}}^T(t)\hfill \\ \multicolumn{1}{c}{}& =\mathit{NP}+\mathit{P}{\mathit{N}}^T+\mathit{G}(t)\mathit{Q}{\mathit{G}}^T(t)\hfill \end{array}$$

(20)

where $\mathit{NP}=E[\mathit{f}{\mathit{r}}^T(t)]$ , $\mathit{N}$ is defined as auxiliary matrix and

$$\mathit{N}=E[\mathit{f}{\mathit{r}}^T(t)]{\mathit{P}}^{-1}$$

(21)

because $\mathit{P}$ is usually positive definite, thus a unique ${\mathit{P}}^{-1}$ exists. Moreover, the nonlinear system function can be approximated by the following equation ¹⁹

$$\mathit{f}(\mathit{x},t)\approx \hat{\mathit{f}}+\mathit{N}(\mathit{x}-\mathit{m})$$

which satisfies the condition for any positive semi-definite matrix $\mathit{S}$

$$E[{\mathit{e}}^T\mathit{se}]=\mathrm{minimum}$$

where $\mathit{e}=\mathit{f}(\mathit{x},t)-\stackrel{·}{\mathit{f}}-\mathit{N}(\mathit{x}-\mathit{m})$ . Thus, the transmit equations of the mean vector and covariance matrix are

$$\begin{array}{c}\stackrel{·}{\mathit{m}}(t)=\hat{\mathit{f}}+\mathit{G}(t)\mathit{b}(t)\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mathit{P}}(t)=\mathit{N}(t)\mathit{P}(\mathit{t})+\mathit{P}(\mathit{t}){\mathit{N}}^T(t)+\mathit{G}(t)\mathit{Q}{\mathit{G}}^T(t)}\end{array}$$

(22)

With the initial conditions $\mathit{m}(0)$ and $\mathit{P}(0)$ , the mean state vector and covariance matrix can be obtained by integrating the above equations.

GHV attitude model statistical linearization

In this part, a statistical linearized model is obtained using the CADET. It is essential first to mention the basic assumption of the CADET. According to Taylor, ¹⁹ all of CADET developments are based on the assumption that the state variables are jointly normal.

Assumption 2

If the state variables $x(t)$ are jointly normal, thus any nonlinearity $f(x(t),t)$ has ¹⁹

$$\begin{array}{cc}\multicolumn{1}{c}{\hat{f}}& =E[f(x(t),t)]=\underset{-\infty }{\overset{+\infty }{\int }}f(x(t),t)P(x)\mathrm{dx}\hfill \\ \multicolumn{1}{c}{N}& =\frac{d\hat{f}}{d{m}_x}\hfill \end{array}$$

where

$$P(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{1}{2}{\left(\frac{x-m_x}{\sigma }\right)}^2}$$

Based on the attitude model presented in section “The attitude model of GHV” and the control scheme proposed in section “Sliding mode controller design of the GHV,” considering the disturbance moment and parameter perturbation, equation (1) became

$$\begin{array}{c}\stackrel{·}{\mathbf{\Omega }}={\mathit{f}}_s(\mathbf{\Omega },\mathit{a})+{\mathit{g}}_s\mathbf{\omega }\hfill \\ \stackrel{·}{\mathbf{\omega }}={\mathit{f}}_f+{\mathit{g}}_f({\mathit{M}}_c+{\mathit{M}}_d)+{\mathit{d}}_f\hfill \end{array}$$

(23)

where ${\mathit{M}}_c$ equals to the control moment in equation (12) and ${\mathit{M}}_d$ is the disturbance moment. The vector $\mathit{a}$ is the parameter vector in ${\mathit{f}}_s$ . Here, we regard the parameter as a combination of the mass, velocity, dynamic pressure, and aero-coefficients, thus $\mathit{a}=[a_1,a_2,a_3{]}^T=[\mathrm{qS}{C}_{L,\alpha }/\mathrm{MV},g/V,\mathrm{qS}{C}_{Y,\beta }/\mathrm{MV}{]}^T$ . To statistically linearize the above system is to linearize ${\mathit{f}}_s(\mathbf{\Omega },\mathit{a})$ , ${\mathit{g}}_s\mathbf{\omega }$ , ${\mathit{f}}_f$ , and the control moment ${\mathit{M}}_c$ . The linearization of the several complex sine terms in ${\mathit{f}}_s$ and ${\mathit{g}}_s\mathbf{\omega }$ can be realized using Taylor expansion.

The method of handling the internal parameters perturbation comes from Jiang et al. ²¹ Considering a normal nonlinear system disturbed by external noises and internal parameter perturbation

$$\stackrel{·}{\mathit{x}}=\mathit{f}(\mathit{x},a)+\mathit{GW}$$

(24)

Remark 1

In this article, the parameter $a$ is equal to

$$a=\overline{a}(1+{\Delta }_a{p}_a)$$

(25)

where $\overline{a}$ is the nominal value. $p_a$ is the deviation percentage, and ${\Delta }_a$ is an indeterminate number from $-1$ to $1$ .

The nonlinear in equation (24) can be linearized as follow

$$\mathit{f}=\hat{\mathit{f}}+\mathit{N}(\mathit{x}-\mathit{m})+{\mathit{N}}_a(a-\overline{a})$$

(26)

Substituting equations (25) and (26) into equation (24), we have

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathit{x}}=\hat{\mathit{f}}+\mathit{N}(\mathit{x}-\mathit{m})+{\mathit{N}}_a\overline{a}{\Delta }_a{p}_a+\mathit{GW}}\\ \multicolumn{1}{c}{=\hat{\mathit{f}}+\mathit{N}(\mathit{x}-\mathit{m})+\mathit{G}\prime W\prime }\end{array}$$

(27)

where $\mathit{G}\prime =[{\mathit{N}}_a\overline{a}{p}_a\mathit{G}]$ and $\mathit{W}\prime =[{\Delta }_a\mathit{W}{]}^T$ . Thus, the new transmit equations of the mean vector and covariance matrix are

$$\begin{array}{c}\stackrel{·}{\mathit{m}}=\hat{\mathit{f}}\hfill \\ \stackrel{·}{\mathit{P}}=\mathit{NP}+\mathit{P}{\mathit{N}}^T+\mathit{G}\prime Q{\mathit{G}\prime }^T\hfill \end{array}$$

(28)

where

$$\mathit{Q}=\left[\begin{array}{cc}{Q}_{{\Delta }_a}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathit{Q}}_W\hfill \end{array}\right]$$

The standard deviation ${\sigma }_{{\Delta }_a}$ is obtained using “Three $\sigma$ Principle.” In this way

$${\sigma }_{{\Delta }_a}=\frac{1}{3}$$

Thus, applying above method on equation (23), we have

$${\mathit{m}}_{f_s}={\mathit{m}}_{\stackrel{·}{\mathbf{\Omega }}}+{\mathit{N}}_{\mathit{a}}(\mathit{a}-\overline{\mathit{a}})$$

Remark 2

In this article, we define the disturbance moment ${\mathit{M}}_d$

$${\mathit{M}}_d=M_d\left(\begin{array}{c}{w}_1\hfill \\ \multicolumn{1}{c}{w_2}\\ \multicolumn{1}{c}{w_3}\end{array}\right)$$

where $M_d$ is the amplitude of the moment disturbance and $w_1,w_2,\mathrm{and}{w}_3$ are consistent with the normal distribution of white noise.

Then, for linearizing the ${\mathit{M}}_c$ , the CADET is used. Here, to reduce chattering phenomenon caused by the SMC controller, the saturation function is used to replace the sign function in equations (9) and (12)

$$\begin{array}{c}\mathrm{sat}(v)=\{\begin{array}{cc}\mathrm{sign}(v)\hfill & \mathrm{if}\left|v\right|>\psi \hfill \\ \multicolumn{1}{c}{\mathrm{kv}}& \mathrm{if}\left|v\right|\le \psi \hfill \end{array}\hfill \end{array}$$

(29)

where $k=1/\psi$ .

Thus, equation (2) and ${\mathit{g}}_s^{-1}$ into equation (9), we have

$$\begin{array}{cc}\multicolumn{1}{c}{p_c=}& \sin (\alpha )\{k_1{\beta }_e+(k_e-k_1)e^{-k_{e}t}{\beta }_e(0)\hfill \\ \multicolumn{1}{c}{}& +{ϵ}_1\mathrm{sa}{t}_1({\mathit{s}}_{1-\beta })-f_{\beta }\}+\cos (\alpha )\cos (\beta )\{k_1{\mu }_e\hfill \\ \multicolumn{1}{c}{}& +(k_e-k_1)e^{-k_{e}t}{\mu }_e(0){ϵ}_1\mathrm{sa}{t}_1({\mathit{s}}_{1-\mu })-f_{\mu }\}\hfill \\ \multicolumn{1}{c}{q_c=}& \{k_1{\alpha }_e+(k_e-k_1)e^{-k_{e}t}{\alpha }_e(0)+{ϵ}_1\mathrm{sa}{t}_1({\mathit{s}}_{1-\alpha })-f_{\alpha }\}\hfill \\ \multicolumn{1}{c}{}& +\sin (\beta )\{k_1{\mu }_e+(k_e-k_1)e^{-k_{e}t}{\mu }_e(0){ϵ}_1\mathrm{sa}{t}_1({\mathit{s}}_{1-\mu })-f_{\mu }\}\hfill \\ \multicolumn{1}{c}{r_c=}& -\cos (\alpha )\{k_1{\beta }_e+(k_e-k_1)e^{-k_{e}t}{\beta }_e(0)\hfill \\ \multicolumn{1}{c}{}& +{ϵ}_1\mathrm{sa}{t}_1({\mathit{s}}_{1-\beta })-f_{\beta }\}+\sin (\alpha )\cos (\beta )\{k_1{\mu }_e\hfill \\ \multicolumn{1}{c}{}& +(k_e-k_1)e^{-k_{e}t}{\mu }_e(0)+{ϵ}_1\mathrm{sa}{t}_1({\mathit{s}}_{1-\mu })-f_{\mu }\}\hfill \end{array}$$

(30)

Similarly, substituting equations (3), (12), (13), and (30) into the $\stackrel{·}{\mathit{w}}$ in equation (1), thus

$$\begin{array}{cc}\multicolumn{1}{c}{\stackrel{·}{p}=}& (k_2+1/T)p_c-k_{2}p-\frac{1}{T}{p}_m+{ϵ}_2\mathrm{sa}{t}_2({\mathit{s}}_{2-p})\hfill \\ \multicolumn{1}{c}{}& +(k_{e2}-k_2)e^{-k_{e2}t}{p}_e(0)\hfill \\ \multicolumn{1}{c}{\stackrel{·}{q}=}& (k_2+1/T)q_c-k_{2}q-\frac{1}{T}{q}_m+{ϵ}_2\mathrm{sa}{t}_2({\mathit{s}}_{2-q})\hfill \\ \multicolumn{1}{c}{}& +(k_{e2}-k_2)e^{-k_{e2}t}{q}_e(0)\hfill \\ \multicolumn{1}{c}{\stackrel{·}{r}=}& (k_2+1/T)r_c-k_{2}r-\frac{1}{T}{r}_m+{ϵ}_2\mathrm{sa}{t}_2({\mathit{s}}_{2-r})\hfill \\ \multicolumn{1}{c}{}& +(k_{e2}-k_2)e^{-k_{e2}t}{r}_e(0)\hfill \end{array}$$

(31)

So far, statistically linearizing the above equations, main objectives are the saturation terms, $\mathrm{sa}{t}_1({\mathit{s}}_1)$ and $\mathrm{sa}{t}_2({\mathit{s}}_2)$ . Here, we regard error variables as state variables, so that differential equations of the error of the angles and angular velocities are necessary to obtain. With using the sliding mode controllers we proposed, the sliding mode variables ${\mathit{s}}_1$ and ${\mathit{s}}_2$ are the function of the error ${\mathbf{\Omega }}_e$ and ${\mathit{w}}_e$ , which make the linearization much easier to be realized.

The differential equations of the error of the angles are

$${\stackrel{·}{\mathbf{\Omega }}}_e={\stackrel{·}{\mathbf{\Omega }}}_c-\stackrel{·}{\mathbf{\Omega }}=\left(\begin{array}{c}{\stackrel{·}{\alpha }}_e\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\beta }}_e}\\ \multicolumn{1}{c}{{\stackrel{·}{\mu }}_e}\end{array}\right)$$

(32)

And, the differential equations of the error of the angular velocities are

$$\begin{array}{c}{\stackrel{·}{\mathit{w}}}_e={\stackrel{·}{\mathit{w}}}_c-\stackrel{·}{\mathit{w}}=\left(\begin{array}{c}{\stackrel{·}{p}}_e\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{q}}_e}\\ \multicolumn{1}{c}{{\stackrel{·}{r}}_e}\end{array}\right)\hfill \end{array}$$

(33)

Moreover, for linearizing the saturation term $\mathrm{sa}{t}_i$ , where $i=1,2$ , some auxiliary functions based on the normal density function in Taylor ¹⁹ are required

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{PF}(w)}& =\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{w}^2}\hfill \\ \multicolumn{1}{c}{\mathrm{PI}(w)}& =\underset{-\infty }{\overset{w}{\int }}\mathrm{PF}(w)\mathrm{dw}\hfill \\ \multicolumn{1}{c}{G(w)}& =\underset{-\infty }{\overset{w}{\int }}\mathrm{PF}(w)\mathrm{dw}\hfill \\ \multicolumn{1}{c}{}& =\mathrm{wPI}(w)+\mathrm{PF}(w)\hfill \end{array}$$

What’s more, the equations of ${\mathit{s}}_1$ and ${\mathit{s}}_2$ are similar, so that we only calculate the linearization of $\mathrm{sa}{t}_1({\mathit{s}}_1)$ as an example. First, we linearize the normal saturation function $\mathrm{sat}(v)$ in equation (29)

$$\hat{f}=\frac{1}{\sqrt{2\pi }\sigma }\underset{-\infty }{\overset{+\infty }{\int }}\mathrm{sat}(v)e^{-\frac{1}{2}{\left(\frac{v-m_v}{\sigma }\right)}^2}\mathrm{dv}$$

(34)

Using the change of variable

$$u=\frac{v-m_v}{\sigma }$$

yields

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{f}}_{\mathrm{sat}(v)}=\frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{+\infty }{\int }}\mathrm{sat}(\sigma u+m_v)e^{-\frac{1}{2}{u}^2}\mathrm{du}}\\ \multicolumn{1}{c}{=-\frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{a_1}{\int }}{e}^{-\frac{1}{2}{u}^2}\mathrm{du}+\frac{k}{\sqrt{2\pi }}}\\ \underset{a_2}{\overset{a_2}{\int }}{e}^{-\frac{1}{2}{u}^2}\mathrm{du}+\frac{1}{\sqrt{2\pi }}\underset{a_2}{\overset{+\infty }{\int }}{e}^{-\frac{1}{2}{u}^2}\mathrm{du}\hfill \\ \multicolumn{1}{c}{=\mathrm{PI}\left(\frac{\psi +m_v}{\sigma }\right)-\mathrm{PI}\left(\frac{\psi -m_v}{\sigma }\right)}\\ +k{m}_v[\mathrm{PI}\left(\frac{\psi +m_v}{\sigma }\right)+\mathrm{PI}\left(\frac{\psi -m_v}{\sigma }\right)]\hfill \\ +k\sigma [\mathrm{PF}\left(\frac{\psi +m_v}{\sigma }\right)-\mathrm{PI}\left(\frac{\psi -m_v}{\sigma }\right)]-k{m}_v\hfill \\ \multicolumn{1}{c}{n_{\mathrm{sat}(v)}=k[\mathrm{PI}\left(\frac{\psi +m_v}{\sigma }\right)+\mathrm{PI}\left(\frac{\psi -m_v}{\sigma }\right)]-k}\end{array}$$

(35)

where

$$a_1=-\frac{\psi +m_v}{\sigma }\mathrm{and}{a}_2=\frac{\psi -m_v}{\sigma }$$

Then, to linearize the saturation term $\mathrm{sa}{t}_1({\mathit{s}}_{1-i})$ , where $i=\alpha ,\beta ,\mu$ . We take the $\mathrm{sa}{t}_1({\mathit{s}}_{1-\alpha })$ as an example

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{f}}_{\mathrm{sa}{t}_1(s_{1-\alpha })}=\frac{1}{\sqrt{2\pi }\sigma }\underset{-\infty }{\overset{+\infty }{\int }}\mathrm{sa}{t}_1(V)e^{-\frac{1}{2}{\left(\frac{V-m_V}{\sigma }\right)}^2}}\\ \multicolumn{1}{c}{V=s_{1-\alpha }={\alpha }_e-{\alpha }_e(0)e^{-k_{e1}t}}\end{array}$$

(36)

Thus

$$\begin{array}{cc}\multicolumn{1}{c}{{\hat{f}}_{\mathrm{sa}{t}_1(s_{1-\alpha })}=}& \frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{+\infty }{\int }}\mathrm{sa}{t}_1(\sigma u+m_{{\alpha }_e})e^{-\frac{1}{2}{u}^2}\mathrm{du}\hfill \\ \multicolumn{1}{c}{}& =-\frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{a{\prime }_1}{\int }}{e}^{-\frac{1}{2}{u}^2}\mathrm{du}+\frac{k}{\sqrt{2\pi }}\hfill \\ \multicolumn{1}{c}{}& \underset{a{\prime }_1}{\overset{a{\prime }_2}{\int }}{e}^{-\frac{1}{2}{u}^2}\mathrm{du}+\frac{1}{\sqrt{2\pi }}\underset{a{\prime }_2}{\overset{+\infty }{\int }}{e}^{-\frac{1}{2}{u}^2}\mathrm{du}\hfill \\ \multicolumn{1}{c}{}& =\mathrm{PI}\left(\frac{-a{\prime }_1}{\sigma }\right)-\mathrm{PI}\left(\frac{a{\prime }_2}{\sigma }\right)\hfill \\ \multicolumn{1}{c}{}& +k{m}_{{\alpha }_e}[\mathrm{PI}\left(\frac{-a{\prime }_1}{\sigma }\right)+\mathrm{PI}\left(\frac{a{\prime }_2}{\sigma }\right)]\hfill \\ \multicolumn{1}{c}{}& +k\sigma [\mathrm{PF}\left(\frac{-a{\prime }_1}{\sigma }\right)-\mathrm{PI}\left(\frac{a{\prime }_2}{\sigma }\right)]-k{m}_{{\alpha }_e}\hfill \\ \multicolumn{1}{c}{n_{\mathrm{sa}{t}_1(s_{1-\alpha })}=}& k[\mathrm{PI}\left(\frac{-a{\prime }_1}{\sigma }\right)+\mathrm{PI}\left(\frac{a{\prime }_2}{\sigma }\right)]-k\hfill \end{array}$$

(37)

where

$$\begin{array}{c}a{\prime }_1=-\frac{\psi +m_{{\alpha }_e}-{\alpha }_e(0)e^{-k_{e1}t}}{\sigma }\hfill \\ a{\prime }_2=\frac{\psi -m_{{\alpha }_e}+{\alpha }_e(0)e^{-k_{e1}t}}{\sigma }\hfill \end{array}$$

Next, for the other complex nonlinear terms in equation (30) like ${ϵ}_1\sin \alpha \mathrm{sa}{t}_1({\mathit{s}}_{1-\beta })$ and ${ϵ}_1\cos \alpha \cos \beta \mathrm{sa}{t}_1({\mathit{s}}_{1-\beta })$ . They are difficult to be linearized statistically by integration with Gaussian probability distribution function or using Taylor expansion method. In order to simplify the linearization, the following conclusion in Taylor ¹⁹ and assumption are used.

Lemma 2

For a nonlinearity that has a linear factor in one variable ¹⁹

$$f(v_1,v_2)=v_{1}g(v_2)$$

Thus, its random input describing functions (RIDFs) is

$$\hat{f}=m_1\hat{g}+p_{12}{n}_g$$

and that the two random component RIDFs are

$$\begin{array}{c}{n}_1=\frac{\partial f}{\partial m_1}=\hat{g}\hfill \\ \multicolumn{1}{c}{n_2=m_1\frac{\partial \hat{g}}{\partial m_2}+p_{12}\frac{\partial n_g}{\partial m_2}}\\ \multicolumn{1}{c}{=m_1{n}_g+p_{12}\frac{\partial n_g}{\partial m_2}}\end{array}$$

Assumption 3

If a angle $\phi$ is small enough, thus $\sin \phi =\phi$ and $\cos \phi =1$ .

Thereby, the RIDFs of the nonlinear terms above can be linearized as follows

$$\begin{array}{l}E[{ϵ}_1\sin \alpha sa{t}_1({\mathit{s}}_{1-\beta })]\approx {ϵ}_{1}E[\alpha sa{t}_1({\mathit{s}}_{1-\beta })]\hfill \\ ={ϵ}_1{m}_{\alpha }{\widehat{f}}_{sa{t}_1(s_{1-\beta })}+P(m_{\alpha }\mathrm{,}{m}_{{\alpha }_e})n_{{\mathit{s}}_{1-\beta }}\hfill \\ E[{ϵ}_1\cos \alpha \cos \beta sa{t}_1({\mathit{s}}_{1-\beta })]\approx {ϵ}_{1}E[sa{t}_1({\mathit{s}}_{1-\beta })]\hfill \\ ={ϵ}_1{\widehat{f}}_{sa{t}_1(s_{1-\beta })}\hfill \end{array}$$

We take the $\Omega$ , ${\Omega }_e$ , ${\mathit{w}}_m$ , ${\mathit{w}}_c$ , and ${\mathbf{\omega }}_e$ as the state vector $\mathit{X}$ . By statistically linearizing the above equations, the covariance describing functions system of the GHV attitude model is obtained

$$\begin{array}{c}{\stackrel{·}{\mathit{m}}}_{\mathit{X}}=\mathit{F}{\mathit{m}}_{\mathit{X}}+{\hat{\mathit{f}}}_{\mathit{X}}+\mathit{b}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mathit{P}}=\mathit{NP}+\mathit{P}{\mathit{N}}^T+\mathit{G}\prime Q{\mathit{G}}^{\prime T}}\end{array}$$

(38)

where

$$\mathit{N}=\mathit{F}+{\mathit{N}}_f$$

and the $\mathit{F}$ stands for linear terms of the state ${\mathit{m}}_{\mathit{X}}$ in equation (23). The vector $\mathit{b}$ is the constant terms in equations (13), (31), and (33). The matrices $\mathit{F}$ , ${\mathit{N}}_a$ , and ${\mathit{N}}_f$ are shown as following, where

$$\begin{array}{c}{b}_1=(k_{e_1}-k_1){\alpha }_e(0)e^{-k_{e1}t}\hfill \\ \multicolumn{1}{c}{b_2=(k_{e_2}-k_2)p_e(0)e^{-k_{e2}t}}\\ \multicolumn{1}{c}{b_3=(k_{e_2}-k_2)q_e(0)e^{-k_{e2}t}}\\ \multicolumn{1}{c}{b_4=(k_{e_2}-k_2)r_e(0)e^{-k_{e2}t}}\end{array}$$

$$\begin{array}{l}\mathit{F}=\left[\begin{array}{ccccccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\frac{1}{T}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{1}{T}*k_1& 0& 0& 0& -\frac{1}{T}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -\frac{1}{T}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\frac{1}{T}& 0& 0& -k_2& 0& 0& 0& 0& 0\\ 0& 0& 0& \left(k_2+\frac{1}{T}\right)*k_1& 0& 0& 0& -\frac{1}{T}& 0& 0& -k_2& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& -\frac{1}{T}& 0& 0& -k_2& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& k_2& 0& 0& 0& 0& 0\\ 0& 0& 0& -k_2*k_1& 0& 0& 0& 0& 0& 0& k_2& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& k_2& 0& 0& 0\end{array}\right]\\ {\mathit{N}}_f=\left[\begin{array}{ccccccccccccccc}{N}_{{\stackrel{.}{\alpha }}_{\alpha }}& N_{{\stackrel{.}{\alpha }}_{\beta }}& N_{{\stackrel{.}{\alpha }}_{\mu }}& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{\alpha }}_p}& 0& N_{{\stackrel{.}{\alpha }}_r}& 0& 0& 0\\ N_{{\stackrel{.}{\beta }}_{\alpha }}& N_{{\stackrel{.}{\beta }}_{\alpha }}& N_{{\stackrel{.}{\beta }}_{\alpha }}& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{\beta }}_p}& 0& N_{{\stackrel{.}{\beta }}_r}& 0& 0& 0\\ N_{{\stackrel{.}{\mu }}_{\alpha }}& N_{{\stackrel{.}{\mu }}_{\alpha }}& N_{{\stackrel{.}{\mu }}_{\alpha }}& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{\mu }}_p}& 0& N_{{\stackrel{.}{\mu }}_r}& 0& 0& 0\\ -N_{{\stackrel{.}{\alpha }}_{\alpha }}& -N_{{\stackrel{.}{\alpha }}_{\beta }}& -N_{{\stackrel{.}{\alpha }}_{\mu }}& 0& 0& 0& 0& 0& 0& -N_{{\stackrel{.}{\alpha }}_p}& 0& -N_{{\stackrel{.}{\alpha }}_r}& 0& 0& 0\\ -N_{{\stackrel{.}{\beta }}_{\alpha }}& -N_{{\stackrel{.}{\beta }}_{\alpha }}& -N_{{\stackrel{.}{\beta }}_{\alpha }}& 0& 0& 0& 0& 0& 0& -N_{{\stackrel{.}{\beta }}_p}& 0& -N_{{\stackrel{.}{\beta }}_r}& 0& 0& 0\\ -N_{{\stackrel{.}{\mu }}_{\alpha }}& -N_{{\stackrel{.}{\mu }}_{\alpha }}& -N_{{\stackrel{.}{\mu }}_{\alpha }}& 0& 0& 0& 0& 0& 0& -N_{{\stackrel{.}{\mu }}_p}& 0& -N_{{\stackrel{.}{\mu }}_r}& 0& 0& 0\\ N_{{\stackrel{.}{p}}_{m_{\alpha }}}& N_{{\stackrel{.}{p}}_{m_{\beta }}}& N_{{\stackrel{.}{p}}_{m_{\mu }}}& 0& N_{{\stackrel{.}{p}}_{m_{{\beta }_e}}}& N_{{\stackrel{.}{p}}_{m_{{\mu }_e}}}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& N_{{\stackrel{.}{q}}_{m_{\beta }}}& N_{{\stackrel{.}{q}}_{m_{\mu }}}& N_{{\stackrel{.}{p}}_{m_{{\alpha }_e}}}& 0& N_{{\stackrel{.}{p}}_{m_{{\mu }_e}}}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ N_{{\stackrel{.}{r}}_{m_{\alpha }}}& N_{{\stackrel{.}{r}}_{m_{\beta }}}& N_{{\stackrel{.}{r}}_{m_{\mu }}}& 0& N_{{\stackrel{.}{r}}_{m_{{\beta }_e}}}& N_{{\stackrel{.}{r}}_{m_{{\mu }_e}}}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ N_{{\stackrel{.}{p}}_{\alpha }}& N_{{\stackrel{.}{p}}_{\beta }}& N_{{\stackrel{.}{p}}_{\mu }}& 0& N_{{\stackrel{.}{p}}_{{\beta }_e}}& N_{{\stackrel{.}{p}}_{{\mu }_e}}& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{p}}_{p_e}}& 0& 0\\ 0& N_{{\stackrel{.}{q}}_{\beta }}& N_{{\stackrel{.}{q}}_{\mu }}& N_{{\stackrel{.}{p}}_{{\alpha }_e}}& 0& N_{{\stackrel{.}{p}}_{{\mu }_e}}& 0& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{p}}_{q_e}}& 0\\ N_{{\stackrel{.}{r}}_{\alpha }}& N_{{\stackrel{.}{r}}_{\beta }}& N_{{\stackrel{.}{r}}_{\mu }}& 0& N_{{\stackrel{.}{r}}_{{\beta }_e}}& N_{{\stackrel{.}{r}}_{{\mu }_e}}& 0& 0& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{r}}_{r_e}}\\ N_{{\stackrel{.}{p}}_{e_{\alpha }}}& N_{{\stackrel{.}{p}}_{e_{\beta }}}& N_{{\stackrel{.}{p}}_{e_{mu}}}& 0& N_{{\stackrel{.}{p}}_{e_{{\beta }_e}}}& N_{{\stackrel{.}{p}}_{e_{{\mu }_e}}}& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{p}}_{e_{p_e}}}& 0& 0\\ 0& N_{{\stackrel{.}{q}}_{e_{\beta }}}& N_{{\stackrel{.}{q}}_{e_{\mu }}}& N_{{\stackrel{.}{q}}_{e_{{\alpha }_e}}}& 0& N_{{\stackrel{.}{q}}_{e_{{\mu }_e}}}& 0& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{q}}_{e_{q_e}}}& 0\\ N_{{\stackrel{.}{r}}_{e_{\alpha }}}& N_{{\stackrel{.}{r}}_{e_{\beta }}}& N_{{\stackrel{.}{r}}_{e_{\mu }}}& 0& N_{{\stackrel{.}{r}}_{e_{{\beta }_e}}}& N_{{\stackrel{.}{r}}_{e_{{\mu }_e}}}& 0& 0& 0& 0& 0& 0& 0& 0& N_{{\stackrel{.}{r}}_{e_{r_e}}}\end{array}\right]\\ {\mathit{N}}_a={\left[\begin{array}{ccccccccccccccc}{N}_{{\stackrel{.}{\alpha }}_{a_1}}& 0& N_{{\stackrel{.}{\mu }}_{a_1}}& -N_{{\stackrel{.}{\alpha }}_{a_1}}& 0& -N_{{\stackrel{.}{\mu }}_{a_1}}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ N_{{\stackrel{.}{\alpha }}_{a_2}}& N_{{\stackrel{.}{\beta }}_{a_2}}& N_{{\stackrel{.}{\mu }}_{a_2}}& -N_{{\stackrel{.}{\alpha }}_{a_2}}& -N_{{\stackrel{.}{\beta }}_{a_2}}& -N_{{\stackrel{.}{\mu }}_{a_2}}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& N_{{\stackrel{.}{\beta }}_{a_3}}& 0& 0& -N_{{\stackrel{.}{\beta }}_{a_3}}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]}^T\\ \begin{array}{l}{\mathbf{m}}_{\mathbf{X}}={\left(\begin{array}{ccccccccccccccc}{m}_{\alpha }& m_{\beta }& m_{\mu }& m_{{\alpha }_e}& m_{{\beta }_e}& m_{{\mu }_e}& m_{p_m}& m_{q_m}& m_{r_m}& m_p& m_q& m_r& m_{p_e}& m_{q_e}& m_{r_e}\end{array}\right)}^T\hfill \end{array}\\ \begin{array}{l}{\widehat{\mathit{f}}}_X={\left(\begin{array}{ccccccccccccccc}{\widehat{f}}_{\alpha }& {\widehat{f}}_{\beta }& {\widehat{f}}_{\mu }& {\widehat{f}}_{{\alpha }_e}& {\widehat{f}}_{{\beta }_e}& {\widehat{f}}_{{\mu }_e}& {\widehat{f}}_{p_m}& {\widehat{f}}_{q_m}& {\widehat{f}}_{r_m}& {\widehat{f}}_p& {\widehat{f}}_q& {\widehat{f}}_r& {\widehat{f}}_{p_e}& {\widehat{f}}_{q_e}& {\widehat{f}}_{r_e}\end{array}\right)}^T\hfill \end{array}\\ \begin{array}{l}\mathit{b}={\left(\begin{array}{ccccccccccccccc}0& 0& 0& 0& 0& 0& 0& \frac{b_1}{T}& 0& b_2& b_3+\left(k_2+\frac{1}{T}\right)*b_1& b_4& -b_2& -b_3-k_2*b_1& -b_4\end{array}\right)}^T\hfill \end{array}\end{array}$$

Extended statistical model with colored noise

Considering the methodology of the CADET in section “Methodology of the CADET,” equation (22) is obtained with the white noise $\mathit{w}(t)$ . However, in engineering practice, the noise encountered is often not white noise. So that a extended statistical model is presented if the colored noise is generated by white noise. ²⁹ For a colored noise $\xi (t)$

$$\stackrel{·}{\xi }(t)={\mathit{f}}_{\xi }(\mathit{X},w,t)+{\mathit{G}}_{\xi }w(t)$$

(39)

where $\mathit{w}(t)$ is a white noise satisfied $\mathit{w}(t)=b+r(t)$ . Thus, a extended state ${\mathit{X}}_{\mathit{w}}$ can be defined as follows

$${\mathit{X}}_w=\left(\begin{array}{c}\mathit{X}\hfill \\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{\xi (t)}\end{array}\right)$$

Then

$$\begin{array}{c}{\stackrel{·}{\mathit{X}}}_w=\left[\begin{array}{c}\mathit{f}({\mathit{X}}_w,w,t)\hfill \\ \multicolumn{1}{c}{\cdots \cdots \cdots \cdots }\\ \multicolumn{1}{c}{{\mathit{f}}_{\xi }({\mathit{X}}_w,w,t)}\end{array}\right]+\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{{\mathit{G}}_{\xi }w(t)}\end{array}\right]\hfill \\ \multicolumn{1}{c}{={\mathit{f}}_w({\mathit{X}}_w,w,t)+{\mathit{G}}_{w}w(t)}\end{array}$$

Therefore, based on the knowledge of section “Methodology of the CADET,” the transmit equations of the extended model is

$$\begin{array}{c}{\stackrel{·}{\mathit{m}}}_w(t)={\hat{\mathit{f}}}_w+{\mathit{G}}_w(t)\mathit{b}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{P}}}_w(t)={\mathit{N}}_w{\mathit{P}}_w+{\mathit{P}}_w{\mathit{N}}_w^T+{\mathit{G}}_w{\mathit{QG}}_w^T}\end{array}$$

(40)

Simulation

To evaluate the effectiveness and performance of CADET for the attitude model of the GHV using the proposed SMC controller, numerical simulations have been carried out based on the model described by equations (38) and (40). For the purpose of comparison, 500 Monte Carlo tests have also been carried out with same white noises. For all simulations presented in this section, the vehicle physical parameters including the mass and aerodynamic coefficients are referred to Keshmiri et al., ²⁴ which are listed in Table 1.

Table 1.

The physical parameters of the GHV.

Parameters	Value
Mass	140,000 lbs
Mach number	6
Height	30 km
Wing reference area	3603 ft2

GHV: generic hypersonic vehicle.

We assume that the constant desired aerodynamic angle command is ${\mathit{\Phi }}_c=[{\alpha }_c,{\beta }_c,{\mu }_c{]}^T=[3\circ ,0\circ ,2\circ {]}^T$ and the initial state is ${\mathit{\Phi }}_0=[0\circ ,2\circ ,0\circ {]}^T$ . There the same controller parameters are used in the three attitude angle channels and angular velocity channels. The parameters of the proposed SMC are listed in Table 2. The time constant $T$ in equation (13) is $0.01\mathrm{s}$ , and the flight-path angle is assumed as $\gamma =0\mathrm{degree}$ .

Table 2.

The parameters of the proposed SMC.

Slowloop parameters	Value	Fastloop parameters	Value
$k_1$	$15$	$k_2$	$15$
${ϵ}_1$	$0.1$	${ϵ}_2$	$0.1$
${\psi }_1$	$0.1$	${\psi }_2$	$0.1$
$k_{e1}$	$2.4$	$k_{e2}$	$2.4$

SMC: sliding mode control.

According to the different kinds of disturbance moment, simulations are divided into the following three cases:

• Case 1: normal Gaussian random disturbance moment and parameter perturbation are analyzed.

• Case 2: parameter perturbation and Gaussian noise with mean 1 and variance 2 disturbance moment are analyzed.

• Case 3: parameter perturbation and colored noise disturbance moment generated by the same noise in Case 2 passing through a transfer function $R/(\alpha s+1)$ are analyzed.

where $\alpha =0.1,R=2.5$ . The amplitude value of disturbance moment in each case is listed in Table 3. In this article, we assume the parameter perturbation is caused by a normal Gaussian noise in all simulation cases. The parameter deviation is $p_a=50\%$ .

Table 3.

The moment disturbances in different simulation cases.

Case	Mean	Variation	Bounds	Amplitude
$1$	$0$	$1$	$\pm \infty$	$2\times {10}^5\mathrm{N}\mathrm{m}$
$2$	$1$	$2$	$\pm 2.5$	$2\times {10}^5\mathrm{N}\mathrm{m}$
$3$	$1$	$2$	$\pm 2.5$	$2\times {10}^5\mathrm{N}\mathrm{m}$

The simulation results of tracking performances of the nominal attitude system with the proposed controller are shown in Figure 2. It is noted that the system is stable and its outputs are able to follow the attitude command. In this article, the control inputs are described by the control moment, which is shown in Figure 3.

Figure 2.

Tracking performances of attitude angles via SMC: (a) tracking performance of the angle of attack via SMC, (b) tracking performance of the sideslip angle via SMC, and (c) tracking performance of the bank angle via SMC.

Figure 3.

The responses of control torques via SMC.

In Figure 4, we show the mean of the attitude tracking error of 500 Monte Carlo tests and the results of the statistical system with external moment disturbances in Case 1. In Figure 5, the standard deviations of CADET and Monte Carlo tests are presented. From those figures, we can see that magnitudes of mean error and standard deviation error are about ${10}^{-3}$ degree. It means the CADET has a similar accuracy compared with Monte Carlo test. For Case 2, from Figures 6 and 7, for any Gaussian white noise the CADET also has a good performance. For colored noise which can be described by equation (39), the comparison of the CADET and Monte Carlo tests in Figures 8 and 9 shows the CADET is able to solve the colored noise and has a good engineering practical significance. From the above simulation diagrams, the CADET method results in one simulation has the same accuracy with 500 Monte Carlo if the mean and the variance of the disturbance are known.

Figure 4.

Case 1: comparison of tracking errors via CADET and Monte Carlo with external moment disturbances and parameter perturbation—(a) mean of the angle of attack error, (b) mean of the sideslip angle error, and (c) mean of the bank angle error.

Figure 5.

Case 1: comparison of standard deviations via CADET and Monte Carlo with external moment disturbances and parameter perturbation—(a) standard deviation of the angle of attack error, (b) standard deviation of the sideslip angle error, and (c) standard deviation of the bank angle error.

Figure 6.

Case 2: comparison of tracking errors via CADET and Monte Carlo with external moment disturbances and parameter perturbation—(a) mean of the angle of attack error, (b) mean of the sideslip angle error, and (c) mean of the bank angle error.

Figure 7.

Case 2: comparison of standard deviations via CADET and Monte Carlo with external moment disturbances and parameter perturbation—(a) standard deviation of the angle of attack error, (b) standard deviation of the sideslip angle error, and (c) standard deviation of the bank angle error.

Figure 8.

Case 3: comparison of tracking errors via CADET and Monte Carlo with external moment disturbances and parameter perturbation—(a) mean of the angle of attack error, (b) mean of the sideslip angle error, and (c) mean of the bank angle error.

Figure 9.

Case 3: comparison of standard deviations via CADET and Monte Carlo with external moment disturbances and parameter perturbation—(a) standard deviation of the angle of attack error, (b) standard deviation of the sideslip angle error, and (c) standard deviation of the bank angle error.

The differences in above figures may be caused by the following two points. The first reason is that the statistical linearization of the nonlinear term of the GHV attitude system is not accurate. During the process of linearization, we ignore the high-order items in Taylor expansion which will produce some errors. The second reason may be due to the poor accuracy of the noise mean and variance. In CADET, we only provide the mean and variance of noise to describe the effects of disturbances on the system. Therefore, inaccurate statistical characteristics of noise will cause differences between CADET and Monte Carlo tests. However, we can also directly change the mean and variance of noises in statistical model, to quickly analyze the system’s response under a new situation.

Conclusion

Considering the extra moment disturbances and internal parameters perturbation, the CADET is applied for the GHV attitude model with the proposed SMC strategy in this article. Compared with the result of the 500 Monte Carlo tests, the CADET’s result shows the same variation trend under the different kinds of noise, but the CADET just uses one-time simulation instead of 500 Monte Carlo tests. Therefore, the CADET can save more computer running time than the conventional Monte Carlo method. In engineering practice, it is better to use CADET and Monte Carlo together in the controller design phase. Moreover, CADET does not require the system to be continuously differentiable, so that many nonlinear controllers can be analyzed. The analysis of other complex guidance or control methods with noises would be an interesting possible future work.