Control law design of hypersonic vehicles using the elastic surrogate model

Abstract

This paper presents a novel control design strategy for the hypersonic vehicle using the elastic surrogate model. First, the parametric model is established for the rigid mode of the hypersonic vehicle based on the engineering estimation and panel methods. Then, the beam surrogate model is applied to identify the elastic mode of the hypersonic vehicle, and the complete parametric model including the rigid and elastic modes is obtained accordingly. Afterward, the control-relevant model is acquired based on the Morris sensitivity analysis method. Furthermore, the control system using the surrogate model is proposed for the hypersonic vehicle to suppress elastic disturbances and maintain system stability. Finally, an illustrative example of the hypersonic vehicle is provided to verify the effectiveness of the presented methods.

Introduction

The air-breathing hypersonic vehicle is considered to be the critical stage for achieving reliable affordable access to space. Compared with traditional vehicles, hypersonic vehicles are featured by significant aero-themo-elastic-propulsion interactions, which bring the challenges for the control system design. In particular, the air-breathing hypersonic vehicle tends to use the lightweight material to decrease the takeoff weight; however, this will result in the strong coupling between the elastic mode and rigid mode. Not only that, the mutual actions between the elasticity and propulsion make the forebody position changes, leading to the change in the shock wave of the forebody to affect the inlet airflow and propulsive efficiency. Also, the structural flexibility has a significant impact on the control action, thus changing the elevator deflection angle and effective control area. As a result, the elastic dynamics need to be fully considered for the hypersonic vehicle to discern the flexible effect on system stability and overall performances.

The hypersonic vehicle is a very complex system when considering the elastic dynamics, and some modeling methods based on the physical principles and experimental data were used to construct the rigid-elastic coupling model of the hypersonic vehicle. In this aspect, the mean axis method was presented in Waszak and Schmidt ¹ to establish the fuselage/structure-coupling model of the elastic aircraft without reflecting the inherent coupling characteristics between the normal acceleration and elastic deformation of the fuselage. Afterward, an integrated modeling method with the airframe/propulsion/elasticity couplings was proposed in Chavez and Schmidt, ² and the classical two-dimensional Newton theory was used to estimate aerodynamic forces and moments, whereas the one-dimensional Rayleigh flow was applied to construct the simplified engine model. Also, the lumped mass method was adopted to describe the features of structural dynamics, and the coupling between rigid bodies and elastic modes was emerged into aerodynamic and propulsive forces. Furthermore, two cantilever beams were employed to describe the structural deformations with respect to the forebody and afterbody bodies in Bolender et al., ³ and, accordingly, an elastic model was established by considering the coupling relationship between structural elasticity and propulsion system to embody the coupling characteristics among subsystems. Following that, some other factors are considered in the modeling process, such as viscous effects, ⁴ unsteady flow, ⁵ and aerodynamic heating, ⁶ thereby improving the modeling accuracy. Beyond this, the modeling strategies using computational fluid dynamics (CFD) tool and experimental data were more and more widely used to acquire the elastic results which can ameliorate the mechanism model. ⁷

On the other hand, the aerodynamic heating will result in the change in thermal stress and material modulus, thereby affecting the system stability. ⁸ To this end, the effects of aerodynamic thermoelasticity were considered in Mackunis et al., ⁹ and the elastic model was described as a linear parameter varying (LPV) model of the state equation matrix regarding the temperature, and further a continuous robust tracking controller was designed using the Lyapunov method. In addition, the coupling influence between the rigid body and low frequency elastic modes was analyzed for the control system design in Buschek and Calise, ¹⁰ and the integrated controller was designed based on the robust control method when considering the elasticity as the uncertain disturbances. Also, the linear quadratic regulator (LQR) controller in combination with the adaptive structure filter was proposed to restrain the elastic mode in Levin et al., ¹¹ and the results showed that the presented controller was satisfactory for guaranteeing system stability. However, the elastic control issues of the hypersonic vehicle are difficult to deal with, and this is because the elastic mode is related with the flight states and external conditions. For instance, the adaptive structure filter depends on the elastic frequency and damping ratio, which are not easy to acquire, and also the elastic mode of the airframe is subjected by the surface temperature and pressure in relation to the aerodynamic forces, thus affecting the control efficiency and accuracy. As a result, this paper presents a novel control design strategy for the hypersonic vehicle using the elastic surrogate model, and the developed surrogate model can adaptively match with the elastic dynamics, such that the performance of the elastic control is dramatically improved accordingly. Not only that, the control law can remove the coupling influence between the rigid and elastic modes by constructing the rigid surrogate model and elastic surrogate model, respectively. In addition, the notch filter and robust control law are designed to coordinate the mutual action of the rigid and elastic mode, and the robust adaptability of the control system is improved for the hypersonic vehicle.

The outline of this paper is provided as follows: “Rigid modeling for hypersonic vehicles” section deals with building the parametric model of the hypersonic vehicle based on the mechanism modeling method. “Elastic mechanism model for hypersonic vehicles” section relates to the surrogate model of the elastic mode using the vibration beam theory. “Elastic surrogate model for hypersonic vehicles” section involves the performance evaluation and control law design associated with the obtained surrogate model. “Adaptive control law design for hypersonic vehicles” section provides an illustrative example to verify the effectiveness of the proposed elastic controller. Some concluding remarks are given in the final section.

Rigid modeling for hypersonic vehicles

The geometry parameters of the hypersonic vehicle mainly adopt the class and shape function transformation method, ¹² and the advantage of this method is that the application of less parameters can approximately describe the smoother curves which construct the complicated configuration of the waverider. In particular, the class function is regarded as a two-dimensional curve function with two parameters N₁ and N₂

$$C_{N_2}^{N_1}\left(\overline{x}\right)=\frac{1}{2^{N_1+N_2}}{\overline{x}}^{N_1}{\left(1-\overline{x}\right)}^{N_2},x\in \left[0,1\right]$$

(1)

The transverse section of the hypersonic vehicle is obtained using the single parameter N₁ because the geometry of the hypersonic vehicle is symmetrical. After that, the three-dimensional surface is generated by stretching the transverse section. Furthermore, the geometric shape of the hypersonic vehicle is composed of different parts, including the fuselage, wing, and engine which need to be handled, respectively, using the class and shape function transformation method. Accordingly, the typical parameterized shape of the hypersonic vehicle is shown in Figure 1. ¹³

Figure 1.

Typical parameterized shape of hypersonic waverider.

After identifying the geometric shape of the hypersonic vehicle, the following work is to estimate the aerodynamic forces using the engineering methods. In particular, the panel method is used to divide the contour surfaces into the triangular surfaces including the three vertices Q₁, Q₂, and Q₃. Accordingly, the normal vector n, centroid coordinates Q_cg, and area of the surface element S_pan are calculated as follows

$$\{\begin{array}{l}\mathbf{n}=\frac{{\mathbf{T}}_1\times {\mathbf{T}}_2}{\Vert {\mathbf{T}}_1\times {\mathbf{T}}_2\Vert }\hfill \\ Q_{\text{cg}}=\frac{Q_1+Q_2+Q_3}{3}\hfill \\ S_{\text{pan}}=\frac{1}{4}\sqrt{(a+b+c)(a+b-c)}\sqrt{(a-c+b)(b+c-a)}\hfill \end{array}$$

(2)

where $\mathbf{n}$ is the normal vector of the surface element, ${\mathbf{T}}_1=Q_2-Q_1,{\mathbf{T}}_2=Q_3-Q_2,a=\left|Q_2-Q_1\right|,b=\left|Q_3-Q_2\right|$ , and $c=\left|Q_1-Q_3\right|$ . Furthermore, the pressure of the surface element is expressed as

$$P_i=C_{\mathit{pi}}{q}_{c,i}+P_{\infty }$$

(3)

where $P_{\infty }$ and $q_{c,i}$ denote the local static pressure and dynamic pressure, respectively, whereas C_pi indicates the aerodynamic coefficients of each element using the estimation method. Afterward, the total aerodynamic force and moment are obtained as

$$\{\begin{array}{l}{F}_x=-{\displaystyle \sum P_i}{n}_{\mathit{xi}}{S}_{\text{pan,}i}\\ F_z=-{\displaystyle \sum P_i}{n}_{\mathit{zi}}{S}_{\text{pan,}i}\\ M_y={\displaystyle \sum P_i}\left(d_{\mathit{xi}}{n}_{\mathit{zi}}-d_{\mathit{zi}}{n}_{\mathit{xi}}\right)S_{\text{pan,}i}\end{array}$$

(4)

where ${\mathbf{d}}_i=d_{\mathit{xi}}{\mathbf{i}}_{\text{b}}+d_{\mathit{yi}}{\mathbf{j}}_{\text{b}}+d_{\mathit{zi}}{\mathbf{k}}_{\text{b}}$ indicates the distance vector from the vehicle centroid to the surface centroid, and n_xi, n_yi, and n_zi are the components of $\mathbf{n}$ . Furthermore, the aerodynamic forces are decomposed along the body axis to obtain the following relationships

$$\{\begin{array}{l}{F}_x=L\sin \hspace{0.17em}\alpha -D\cos \hspace{0.17em}\alpha \hfill \\ F_z=-L\cos \hspace{0.17em}\alpha -D\sin \hspace{0.17em}\alpha \hfill \end{array}$$

(5)

where α is the angle of attack and L and D denote the lift and drag, respectively. Beyond this, the engine model is simplified to a two-dimensional structure consisting of a front body compression section, an inlet, an isolation chamber, a combustion chamber, an inner nozzle, and an outer nozzle coupled with the fuselage. Correspondingly, the simplified propulsion system model is provided as ¹⁴

$${\mathbf{F}}_T={\dot{m}}_{\text{a}}{\mathbf{V}}_1-{\dot{m}}_{\text{a}}\left(1+{\delta }_T{f}_{\text{st}}\right){\mathbf{V}}_{\text{5}}+P_5{\mathbf{S}}_{\text{th5}}+P_1{\mathbf{S}}_{\text{th1}}$$

(6)

where ${\mathbf{S}}_{\text{th5}}$ is the area vector of the nozzle exit, ${\mathbf{S}}_{\text{th1}}$ denotes the entrance area vector of the inlet, P₁ and P₅ indicate the pressures of the air flow at the inlet and nozzle, respectively, ${\mathbf{V}}_1$ and ${\mathbf{V}}_{\text{5}}$ denote the velocities of the air flow at the inlet and nozzle, respectively, f_st indicates the fuel-air ratio, δ_T is fuel equivalence ratio, and ${\dot{m}}_{\text{a}}$ denotes the air mass flow rate. Also, the pitching moment produced by the thrust is provided as

$${\mathbf{M}}_T={\mathbf{r}}_{\mathit{TG}}\times {\mathbf{F}}_T$$

(7)

where ${\mathbf{r}}_{\mathit{TG}}$ denotes the distance vector from the center of gravity to the thrust center. Based on equations (4) to (6), the aerodynamic forces and moments can be obtained using the complex parametric model of the hypersonic vehicle, and the model database can be constructed as the functions in relation to the flight condition and control inputs. They are expressed by

$$\{\begin{array}{l}L=q_s{S}_r{C}_L\left(\mathit{Ma},\alpha ,h,{\delta }_e,{\delta }_T\right)\\ D=q_s{S}_r{C}_D\left(\mathit{Ma},\alpha ,h,{\delta }_e,{\delta }_T\right)\\ T=q_s{S}_c{C}_T\left(\mathit{Ma},\alpha ,h,{\delta }_e,{\delta }_T\right)\\ M_y=M_T+q_s{S}_r{L}_r{C}_m\left(\mathit{Ma},\alpha ,h,{\delta }_e,{\delta }_T\right)\end{array}$$

(8)

where δ_e denotes the deflection angle of the elevon, Ma indicates the flight Mach, h expresses the flight altitude, S_r is the reference area, S_c denotes the capture area of the inlet, L_r indicates the reference length, q_s is the dynamic pressure, and $C_L,C_D,C_T,C_m$ denote the model coefficients regarding the lift L, drag D, thrust T, and pitch moment M_y.

However, the forms in equation (8) need to be specified as the polynomial surrogate model depending on state and control variables. Furthermore, the goodness of fit between the original model and polynomial surrogate model is calculated, and the forms and cross forms, which have an insignificant effect on the goodness of fit, will be eliminated. As a result, the feasible expression forms are acquired by keeping the important forms. In this work, the forms of state and control variables are analyzed, and the polynomial surrogate model is provided as

$$\{\begin{array}{l}{C}_L=C_{L0}+C_{L\alpha }\alpha +C_{L\delta }{\delta }_e+C_{\mathit{LM}\alpha }Ma·\alpha +C_{\mathit{LM}}Ma+C_{L2\delta }{\delta }_{e}2\hfill \\ C_D=C_{D0}+C_{D2\alpha }{\alpha }^2+C_{D2\delta }{\delta }_{e}2+C_{D\alpha \delta }\alpha ·{\delta }_e+C_{D\alpha }\alpha +C_{D\delta }{\delta }_e+C_{\mathit{DM}}Ma\hfill \\ \begin{array}{c}{C}_m=C_{m0}+C_{m\delta e}{\delta }_e+C_{m\alpha }\alpha +C_{\mathit{mM}}Ma+C_{m\delta t}{\delta }_T+C_{\mathit{mM}\delta e}Ma{\delta }_e\\ \hspace{1em}\hspace{1em}\hspace{1em}+C_{m\alpha \delta e}\alpha {\delta }_e+C_{\mathit{mM}\delta t}Ma{\delta }_T+C_{m2\alpha }{\alpha }^2+C_{m2\delta }{\delta }_{e}2+C_{m2M}M{a}^2\end{array}\hfill \end{array}$$

(9)

where $C_{L0},C_{D0}$ , and $C_{m0}$ are coefficients concerning the lift, drag, and pitch moment. After that, the rigid model is established as

$$\{\begin{array}{l}\dot{V}=\frac{T\cos \hspace{0.17em}\alpha -D}{m}-g\sin \hspace{0.17em}(\theta -\alpha )\hfill \\ \dot{\alpha }=q-\frac{T\sin \hspace{0.17em}\alpha +L}{\mathit{mV}}+\frac{g}{V}\cos \hspace{0.17em}(\theta -\alpha )\hfill \\ \dot{h}=V\sin \hspace{0.17em}(\theta -\alpha )\hfill \\ \dot{q}=\frac{M_y}{I_y}\hfill \\ \dot{\theta }=q\hfill \end{array}$$

(10)

where V, θ, and q denote the flight speed, pitch angle, and pitch angle rate, respectively, and m, I_y, and g are the mass, moment of inertia, and gravitational constant, respectively.

Elastic mechanism model for hypersonic vehicles

The slender fuselage of the hypersonic vehicle is considered to be two cantilever beams which represent the forebody and afterbody and further assume that the elastic vibration occurs in the central axis of each section of the beam without the change in the moment of inertia and the shear deformation. Furthermore, the micro section of the beam dx driven along the z direction is expressed by

$$\overline{m}\left(x\right)\frac{{\partial }^{2}z}{\partial t^2}+\frac{\partial Q}{\partial x}=f\left(x,t\right)$$

(11)

where $\overline{m}\left(x\right)$ is the mass density, Q denotes the shear force, M indicates the bending moment, and f(x, t) expresses the force acting on the surface of the beam structure. Furthermore, considering $Q=\frac{\partial M}{\partial x}$ , equation (11) is rewritten as

$$\overline{m}\left(x\right)\frac{{\partial }^{2}z}{\partial t^2}+\frac{{\partial }^{2}M}{\partial x^2}=f\left(x,t\right)$$

(12)

Also, the relation between the bending moment and deflection can be expressed as

$$M\left(x,t\right)=\mathit{EI}\left(x\right)\frac{{\partial }^{2}z\left(x,t\right)}{\partial x^2}$$

(13)

where E is Young modulus and I denotes the moment of inertia of the section along with the x axis. Substituting equation (13) into equation (12) yields

$$\overline{m}\left(x\right)\frac{{\partial }^{2}z\left(x,t\right)}{\partial t^2}+\frac{{\partial }^2}{\partial x^2}\left[\mathit{EI}\left(x\right)\frac{{\partial }^{2}z\left(x,t\right)}{\partial x^2}\right]=f\left(x,t\right)$$

(14)

Furthermore, z(x, t) is reshaped as

$$z\left(x,t\right)=\phi \left(x\right)\eta \left(t\right)$$

(15)

where $\phi \left(x\right),\eta \left(t\right)$ indicate the mode shape function and generalized coordinate, respectively. When considering free vibration of the beam with $f(x,t)=0$ and substituting equation (15) into equation (14), we get

$$\frac{d^2\eta \left(t\right)}{d{t}^2}+{\omega }^2\eta \left(t\right)=0$$

(16)

$${\omega }^2=\frac{1}{\overline{m}\left(x\right)\phi \left(x\right)}\frac{d^2}{d{x}^2}\left[\mathit{EI}\left(x\right)\frac{d^2\phi \left(x\right)}{d{x}^2}\right]$$

(17)

The general solution for equation (16) is expressed by

$$\eta \left(t\right)=A\sin \hspace{0.17em}\omega t+B\cos \hspace{0.17em}\omega t=C\sin \hspace{0.17em}\left(\omega t+a\right)$$

(18)

where A and B are the integral constants determined by two initial values.

Furthermore, considering the forced vibration of the beam of equation (14), the two-order differential equations of η are obtained as ¹⁵

$$\begin{array}{cc}{\ddot{\eta }}_k\left(t\right)+2\zeta {\omega }_k{\dot{\eta }}_k\left(t\right)+{\omega }_k^2{\eta }_k\left(t\right)=N_k\left(t\right)& k=1,2,3\dots \end{array}$$

(19)

where ${\eta }_k,{\omega }_k,\zeta$ are, respectively, the generalized coordinate, natural frequency, and damping coefficient with regard to the Kth beam. In turn, the generalized force regarding the Kth beam $N_k(t)$ is defined by

$$N_k\left(t\right)=\frac{{\displaystyle {\int }_0^l{\phi }_k\left(x\right)f\left(x,t\right)\text{d}x}}{{\displaystyle {\int }_0^l\overline{m}\left(x\right){\phi }_k^2\left(x\right)\text{d}x}}$$

(20)

where ${\phi }_k(x)$ denotes the mode shape function regarding the Kth beam which is decided by equation (17). However, ${\phi }_k(x)$ is related to $\overline{m}\left(x\right)$ and $\mathit{EI}\left(x\right)$ , such that it is difficult to identify the specific form. As a result, the finite element method needs to be applied to estimate the elastic modes of the heterogeneous beams. Based on the differential equation of the beam deflection curve $z{\left(x\right)}^{\left(4\right)}=0$ , the vibration mode function of the unit beam is expressed by

$$z\left(x\right)=c_1+c_2\frac{x}{l}+c_3{\left(\frac{x}{l}\right)}^2+c_4{\left(\frac{x}{l}\right)}^3$$

(21)

where l denotes the length of the unit beam, and $c_1,c_2,c_3,c_4$ are the according constants determined by

$$\{\begin{array}{l}{c}_1=z_i\\ c_2=l{\vartheta }_i\\ c_1+c_2+c_3+c_4=z_j\\ c_2+2c_3+3c_4=l{\vartheta }_j\end{array}$$

(22)

where $z\left(0\right)=z_i,\dot{z}\left(0\right)={\vartheta }_i,z\left(l\right)=z_j,\dot{z}\left(l\right)={\vartheta }_j$ , and the vibration mode function of the unit beam is identified by the linear combination of nodal displacements $z_i,z_j,{\vartheta }_i,{\vartheta }_j$ , and it is shown as

$$z\left(x,t\right)=\mathbf{N}\left(x\right){\mathbf{q}}_e\left(t\right)$$

(23)

where ${\mathbf{q}}_{\mathit{ei}}\left(t\right)={\left[z_i\left(t\right),{\vartheta }_i\left(t\right),z_j\left(t\right),{\vartheta }_j\left(t\right)\right]}^T$ , and the shape function vector $\mathbf{N}\left(x\right)=\left[N_{\mathit{iz}}\left(x\right),N_{i\vartheta }\left(x\right),N_{\mathit{jz}}\left(x\right),N_{j\vartheta }\left(x\right)\right]$ is computed as

$$\{\begin{array}{c}{N}_{\mathit{iz}}\left(x\right)=1-3{\left(\frac{x}{l}\right)}^2+2{\left(\frac{x}{l}\right)}^3\\ N_{i\vartheta }\left(x\right)=x-2l{\left(\frac{x}{l}\right)}^2+l{\left(\frac{x}{l}\right)}^3\\ N_{\mathit{jz}}\left(x\right)=3{\left(\frac{x}{l}\right)}^2-2{\left(\frac{x}{l}\right)}^3\\ N_{j\vartheta }\left(x\right)=-l{\left(\frac{x}{l}\right)}^2+l{\left(\frac{x}{l}\right)}^3\end{array}$$

(24)

Furthermore, the kinetic energy of the unit beam can be expressed as

$$E_{\mathit{ke}}=\frac{1}{2}{\displaystyle {\int }_0^l\rho {\left(\frac{\partial z\left(x,t\right)}{\partial t}\right)}^2}\text{d}x=\frac{1}{2}{\displaystyle {\int }_0^l\rho {\mathbf{N}}^T\mathbf{N}\text{d}x{\dot{\mathbf{q}}}_e=\frac{1}{2}{\dot{\mathbf{q}}}_e^T\mathbf{m}{\dot{\mathbf{q}}}_e}$$

(25)

$$\mathbf{m}={\displaystyle {\int }_0^l\rho {\mathbf{N}}^T\mathbf{N}\text{d}x}=\frac{\rho l}{420}\left[\begin{array}{l}15622l54-13l\hfill \\ 22l4l^{2}13l-3l^2\hfill \\ 5413l156-22l\hfill \\ -13l-3l^2-22l4l^2\hfill \end{array}\right]$$

(26)

Similarly, the potential energy of the unit beam is demonstrated as

$$E_{\mathit{pe}}=\frac{1}{2}{\displaystyle {\int }_0^l\mathit{EI}{\left(\frac{{\partial }^{2}z\left(x,t\right)}{\partial t^2}\right)}^2\text{d}x}=\frac{1}{2}{\mathbf{q}}_{\mathit{e}}^T{\displaystyle {\int }_0^l\mathit{EI}{\left({\mathbf{N}}^{\prime \prime }\right)}^T{\mathbf{N}}^{\prime \prime }\text{d}x{\mathbf{q}}_e=\frac{1}{2}{\mathbf{q}}_{\mathit{e}}^T\mathbf{k}{\mathbf{q}}_e}$$

(27)

$${\mathbf{k}}_{\mathbf{a}}={\displaystyle {\int }_0^l\mathit{EI}{\left({\mathbf{N}}^{\prime \prime }\right)}^T{\mathbf{N}}^{\prime \prime }\text{d}x}=\frac{\mathit{EI}}{l^3}\left[\begin{array}{l}126l-126l\hfill \\ 6l4l^2-6l2l^2\hfill \\ -12-6l12-6l\hfill \\ 6l2l^2-6l4l^2\hfill \end{array}\right]$$

(28)

Correspondingly, the total mass matrix and total stiffness matrix of the beam are expressed as follows

$$\{\begin{array}{l}\mathbf{M}={\displaystyle \sum _{i=1}^n{\mathbf{R}}_{i}T}\mathbf{m}{\mathbf{R}}_i\hfill \\ \mathbf{K}={\displaystyle \sum _{i=1}^n{\mathbf{R}}_{i}T}{\mathbf{k}}_{\mathbf{aR}}i\hfill \end{array}$$

(29)

where i is the number of unit beams, n represents the amount of element beams, $\mathbf{q}={\left(z_1,{\vartheta }_1,z_2,{\vartheta }_2\cdots z_n,{\vartheta }_n\right)}^T$ , and ${\mathbf{q}}_{\mathit{ei}}={\mathbf{R}}_i\mathbf{q},i=1,2\cdots n$ . In case of free vibration, the airframe beam can be obtained using the Lagrange equation, and it is expressed by

$$\left({\omega }^2\mathbf{I}-{\mathbf{M}}^{-1}\mathbf{K}\right)\mathbf{q}=0$$

(30)

where ω is the natural frequency. In particular, the node displacements corresponding to the forebody and afterbody cantilever beams are estimated with the constraints $z_i={\vartheta }_i=0$ .

Elastic surrogate model for hypersonic vehicles

The computational results using the finite element methods are very complicated and difficult to use for the control design, and the elastic surrogate model needs to be established to adapt the control-oriented design and analysis. To this end, the surrogate model of the forebody and afterbody beam is considered to provide the analytical expressions for elastic frequencies and vibration modes. For simplicity, assume that the cross section moment of inertia of the airframe beam is fixed with the change in the mass and body temperature. Accordingly, the surrogate model of the elastic frequency and vibration mode is expressed by

$$\{\begin{array}{l}\omega =\omega \left(m,T\right)\hfill \\ \phi =\phi \left(m,T\right)\hfill \end{array}$$

(31)

Specially, the elastic mode of the vibration beam of the forebody is determined by $({\omega }_f,{\phi }_f)$ , whereas the elastic mode of the afterbody beam is decided by $({\omega }_a,{\phi }_a)$ .

Also, the elastic surrogate model is deduced as ⁹

$$\{\begin{array}{c}{\ddot{\eta }}_f=-2\zeta \omega f\hspace{0.17em}{\dot{\eta }}_f-{\omega }_f^2{\eta }_f+N_f\\ {\ddot{\eta }}_a=-2\zeta {\omega }_a{\dot{\eta }}_a-{\omega }_a^2{\eta }_a+N_a\end{array}$$

(32)

where ${\omega }_f,{\omega }_a,N_f,N_a$ represent the natural frequencies and generalized forces with respect to the forebody and afterbody and ζ is the elastic damping coefficient.

Furthermore, the generalized forces can be considered as a function of the flight state and control input, and they are shown as

$$\{\begin{array}{l}{N}_f=q_S{S}_r{C}_{N_f}\hfill \\ N_a=q_S{S}_r{C}_{N_a}\hfill \\ C_{N_f}=q_S{S}_r\left(b_{\mathit{Nf}0}+b_{\mathit{Nf}\alpha }\alpha +b_{\mathit{NfM}}Ma+b_{\mathit{Nf}\eta f\hspace{0.17em}}{\eta }_f+b_{\mathit{Nf}\eta a}{\eta }_{\text{a}}\right)\hfill \\ C_{N_a}=q_S{S}_r\left(b_{\mathit{Na}0}+b_{\mathit{Na}\alpha }\alpha +b_{\mathit{NaM}}Ma+b_{\mathit{Na}\eta f\hspace{0.17em}}{\eta }_f+b_{\mathit{Na}\eta a}{\eta }_{\text{a}}+b_{\mathit{Na}\delta }{\delta }_e\right)\hfill \end{array}$$

(33)

where $C_{N_f}$ and $C_{N_a}$ are the generalized force coefficient with regard to the forebody and afterbody, respectively. Beyond this, two indexes including the variance ratio and mean square error are used to evaluate the matching properties between the mechanism model and surrogate model, and they are shown as

$$\mathit{VAF}=1-\frac{{\Vert y_{\mathit{obs}}-y_{\mathit{fit}}\Vert }_2^2}{{\Vert y_{\mathit{obs}}\Vert }_2^2}$$

(34)

$$\mathit{RMSE}=\frac{{\Vert y_{\mathit{obs}}-y_{\mathit{fit}}\Vert }_2}{N}$$

(35)

where y_obs, y_fit denote the model data corresponding to the elastic mechanism model and surrogate model, respectively.

Based on the aerodynamic, structural, and propulsion model of the hypersonic vehicle, the data regarding the force and moment of the hypersonic vehicle in the whole flight envelope are obtained to determine the specific form of the aerodynamic coefficients in equation (33).

Furthermore, the establishment of the polynomial surrogate model of the force and model consists of four steps: the sample point design, model structure determination, model parameter identification, and model verification. After obtaining the sample space, the structure of the surrogate model is required to determine based on the sample space data and to identify the critical items included in the surrogate model. Specially, the Morris sensitivity analysis method is used to analyze the impact of each item on the results, thereby determining the polynomial form according to the influence level. Beyond this, the fitting results can be verified based on two indexes of goodness of fit in equations (34) and (35).

Adaptive control law design for hypersonic vehicles

Considering the speed and altitude response of the hypersonic vehicle, the steady-state tracking errors are expressed by ¹⁶

$$\{\begin{array}{c}{x}_V={\displaystyle {\int }_0^t\left(\Delta V_{\mathit{ref}}-\Delta V\right)\text{d}t}\\ x_h={\displaystyle {\int }_0^t\left(\Delta h_{\mathit{ref}}-\Delta h\right)\text{d}t}\end{array}$$

(36)

where $\Delta V_{\mathit{ref}}$ and $\Delta h_{\mathit{ref}}$ denote the reference commands regarding V and h. The augmented system based on the linearized model is represented as

$$\{\begin{array}{l}\Delta {\dot{\mathbf{X}}}_{\mathit{aug}}={\mathbf{A}}_{\mathit{aug}}\Delta {\mathbf{X}}_{\mathit{aug}}+{\mathbf{B}}_{1,\mathit{aug}}\mathbf{W}+{\mathbf{B}}_{2,\mathit{aug}}\Delta \mathbf{U}\hfill \\ \Delta \mathbf{Z}={\mathbf{C}}_{1,\mathit{aug}}\Delta {\mathbf{X}}_{\mathit{aug}}\hfill \\ \Delta {\mathbf{Y}}_{\mathit{aug}}={\mathbf{C}}_{2,\mathit{aug}}\Delta {\mathbf{X}}_{\mathit{aug}}\hfill \end{array}$$

(37)

where $\Delta {\mathbf{X}}_{\mathit{aug}}={[\begin{array}{ccc}\Delta \mathbf{X}& x_V& x_h\end{array}]}^T$ and $\mathbf{W}={[\begin{array}{cc}\Delta V_{\mathit{ref}}& \Delta h_{\mathit{ref}}\end{array}]}^T$ , whereas ${\mathbf{A}}_{\mathit{aug}},{\mathbf{B}}_{1,\mathit{aug}},{\mathbf{B}}_{2,\mathit{aug}},{\mathbf{C}}_{1,\mathit{aug}},{\mathbf{C}}_{2,\mathit{aug}}$ indicate the model matrices of the hypersonic vehicle. For equation (37), the control law is designed as ¹⁷

$$\Delta \mathbf{U}={\mathbf{K}}_{\mathbf{c}}\Delta {\mathbf{X}}_{\mathit{aug}}={\mathbf{K}}_H{\mathbf{Q}}^{-1}\Delta {\mathbf{X}}_{\mathit{aug}}$$

(38)

where

$$\left[\begin{array}{l}{\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H+\mathbf{Q}{\mathbf{A}}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T{\mathbf{B}}_{1,\mathit{aug}}\mathbf{Q}{\mathbf{C}}_{1,\mathit{aug}}^T\hfill \\ {\mathbf{B}}_{1,\mathit{aug}}T-\gamma I\mathbf{0}\hfill \\ {\mathbf{C}}_{1,\mathit{aug}}\mathbf{Q}\mathbf{0}-\gamma I\hfill \end{array}\right]<0$$

(39)

$${\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H+\mathbf{Q}{\mathbf{A}}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T-2{\sigma }_s\mathbf{Q}<0$$

(40)

$$\left[\begin{array}{l}-{\omega }_n\mathbf{Q}{\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H\hfill \\ \mathbf{Q}{\mathbf{A}}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T-{\omega }_n\mathbf{Q}\hfill \end{array}\right]<0$$

(41)

$$\left[\begin{array}{l}\sin \hspace{0.17em}\vartheta \left(\begin{array}{l}{\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H+\hfill \\ \mathbf{Q}{A}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T\hfill \end{array}\right)\cos \hspace{0.17em}\vartheta \left(\begin{array}{l}{\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H-\hfill \\ \left(\mathbf{Q}{\mathbf{A}}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T\right)\hfill \end{array}\right)\hfill \\ \cos \hspace{0.17em}\vartheta \left(\begin{array}{l}\mathbf{Q}{\mathbf{A}}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T-\hfill \\ \left({\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H\right)\hfill \end{array}\right)\sin \hspace{0.17em}\vartheta \left(\begin{array}{l}{\mathbf{A}}_{\mathit{aug}}\mathbf{Q}+{\mathbf{B}}_{2,\mathit{aug}}{K}_H+\hfill \\ \mathbf{Q}{\mathbf{A}}_{\mathit{aug}}^T+K_{\mathit{H}}^{\mathit{T}}{\mathbf{B}}_{2,\mathit{aug}}^T\hfill \end{array}\right)\hfill \end{array}\right]<0$$

(42)

where γ, ${\omega }_n$ , and σ_s are the design parameters and $\mathbf{Q}={\mathbf{Q}}^T$ denotes the selected matrix. Based on Yuan, ¹⁷ the closed-system stability can be guaranteed if the control law in equation (38) is used for the system in equation (37). In addition, the designed controller can suppress the disturbances of the uncertainties, but the elastic mode will change the model structure as a result of leading to the presence of the unexpected flight states and control qualities. To this end, the adaptive control system is designed in combination with the control law in equation (38) and the adaptive notch filter based on the surrogate model in order to improve the stability robustness and performance robustness. In particular, the frequency of the notch filter should be matched with the elastic frequency to restrain the elastic response, but the elastic frequency tends to be difficult to measure and to change dynamically with the different flight conditions. Thus, the notch filter is established as

$$\begin{array}{l}H\left(s\right)=\frac{s^2+\varpi s+1}{s^2+\varpi \tau s+{\tau }^2}0<\tau <1\hfill \end{array}$$

(43)

where s indicates the Laplace’s unit, τ is the design parameter, ϖ is the frequency the notch filter and $\varpi =-2\cos \hspace{0.17em}\left(2\pi f_{0}T\right)$ , where f₀ denotes the notch frequency. In fact, ϖ is an important design parameter to filter the adverse effect out of the elastic response, and this parameter can be estimated by the on-line estimator associated with the surrogate elastic model. For instance, the linearized height model is built as

$$h\left(t\right)=\frac{N_{e,k}\left(s\right)}{D_{e,k}\left(s\right)}{\delta }_e\left(t\right)$$

(44)

where $N_{e,k}\left(s\right),D_{e,k}\left(s\right)$ are the model parameters. In this case, $D_{e,k}\left(s\right)$ is expressed by

$$D_{e,k}\left(s\right)=\left[s^2+2{\xi }_r{\omega }_{r}s+{\omega }_r^2\right]\left[s^2+2{\xi }_D{\omega }_{D}s+{\omega }_D^2\right]$$

(45)

where ω and ξ are the frequency and damping, and the subscripts r and D represent the rigid and elastic modes, respectively. Accordingly, the design purpose of the online estimator is to obtain adaptively the frequency ω_D as the design frequency ϖ of the notch filter in equation (43), and ω_D is estimated based on the surrogate elastic model in equations (31) and (32). The design process is provided in Figure 2.

Figure 2.

Design process of control law.

Figure 2 shows that the proposed control system depends on the surrogate models of the rigid and elastic mode, and these surrogate models are identified based on the parametric modeling methods. Specifically, the notch filter is designed to suppress the elastic effect in line with the frequency of the elastic surrogate model, whereas the control-oriented model is used to construct the robust control law to eliminate the uncertainties induced by the modeling errors, external disturbances, and elastic deformation. Different from the traditional controller design process, the designed controller considers fully the rigid-elastic coupling relations by analyzing the model characteristics and building the overall surrogate model, and the notch filter and control law are completely based on the analysis results. As a result, the presented controller can effectively restrain the elastic modal effect and guarantee robust stability of the hypersonic vehicle.

Illustrative example

This paper applies the typical structure of the hypersonic waverider to verify the effectiveness of the control strategy, and the parameters in Figure 1 are selected in Table 1. ¹³

Table 1.

Design parameters of typical structure of hypersonic waverider.

Total length	L (m)	30
Forebody length	L1 (m)	14
Mass	m (kg)	42094
Elevator area	Se (m2)	1.58
Moment of inertia	Iyy (km m2)	979106.44
Engine length	Le (m)	6
Engine height	hi (m)	1
Forebody angle	τ1 (°)	3

According to the calculating data, the coefficients in equation (9) are identified using the least square method, and the polynomial surrogate model of the force and moment of the hypersonic vehicle are obtained accordingly. These are expressed by

$$\{\begin{array}{l}{C}_L=0.0191+5.9994\alpha +0.5140{\delta }_e-0.2382Ma·\alpha +0.0008Ma+0.1930{\delta }_e^2\hfill \\ C_D=0.0622+0.6861{\alpha }^2+5.2791{\delta }_e^2+1.4939\alpha ·{\delta }_e+0.0968\alpha -0.0008{\delta }_e-0.0031Ma\hfill \\ \begin{array}{l}{C}_T=-0.9265-1.5997{\delta }_T+0.0690Ma-7.3694\alpha +49.8126\alpha M{a}^{-1}\\ \hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+11.5948\alpha {\delta }_T+26.9843{\delta }_{T}M{a}^{-1}+2.5405M{a}^{-1}-5.0748{\alpha }^2\end{array}\hfill \\ \begin{array}{l}{C}_m=0.7600-5.0096{\delta }_e+1.7932\alpha -0.1160Ma-0.2502{\delta }_T+0.1595Ma\delta e\\ \hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+0.0102\alpha {\delta }_e-2.3851Ma{\delta }_T+14.6598{\alpha }^2-1.2175{\delta }_e^2+0.0070M{a}^2\end{array}\hfill \end{array}$$

(46)

Furthermore, the elastic surrogate model in relation the cantilever beam is deduced on the basis of the sample point design, model structure determination, model parameter identification, and model verification. Accordingly, the coefficients in equation (35) are determined, and the specific surrogate model of the generalized forces is acquired as

$$\{\begin{array}{c}\begin{array}{l}{C}_{N_f}=1.7320\times {10}^{-4}+\text{0}.\text{0039}\alpha -7.4198\times {10}^{-6}Ma+5.2459\times {10}^{-6}{\eta }_f+9.9208\times {10}^{-6}{\eta }_f\hfill \end{array}\\ \begin{array}{l}{C}_{N_a}=1.4058\times {10}^{-4}+\text{0}.\text{0021}\alpha -1.6499\times {10}^{-5}Ma+9.5780\times {10}^{-4}{\delta }_e-5.7786\times {10}^{-6}{\eta }_a+1.2829\times {10}^{-5}{\eta }_a\hfill \end{array}\end{array}$$

(47)

Furthermore, the comparative results between the surrogate model and original model are shown for the pitch moment coefficient, lift coefficient, drag coefficient, and thrust coefficient are obtained, and they are shown in Figure 3.

Figure 3.

Comparative results between the surrogate model and actual model.

Figure 3 demonstrates that the surrogate model is matched with the original model; thus, the control system is designed based on the surrogate model reflecting the dynamic characteristics of the original. To this end, the adaptive control system is applied with the following control gain

$$\begin{array}{c}{\mathbf{K}}_c=[\begin{array}{l}-0.0812919-2.83-13.22-29890.06480.27\hfill \\ 2.214-7451872.22339.1176344-1.46-6.53\hfill \end{array}\\ \begin{array}{l}\hspace{1em}\hspace{1em}7.73444.25-0.0380.016-0.1682.1870.1393.06\hfill \\ \hspace{1em}-184.8-89940.94-0.373.968-69.52-4.07-77.97\hfill \end{array}]\end{array}$$

(48)

Thus, the response curves regarding the velocity and height commands are provided in Figure 4, and the changes in the control inputs are shown in Figures 5 and 6.

Figure 4.

Response curves regarding the velocity and height commands.

Figure 5.

Changes in the elevon deflection angle.

Figure 6.

Changes in the engine input.

From Figures 4 to 6, we observe that the control results are satisfactory even in case of the presence of the elastic disturbances, and the velocity and altitude tracking errors are very small as the system enters to the steady state. Also, the control inputs keep stable in the control process with the coordinated change with the track errors. This shows that the control action can not only guarantee system stability, but also relieve the effect of the elastic mode.

Conclusion

This paper proposes the control strategy using the elastic surrogate model for the hypersonic vehicle, and the innovative points lie in that the elastic surrogate model is constructed based on the computational data of the finite element method. In addition, the control law is designed with the notch filter which is built with the elastic surrogate model. Furthermore, the frequency of the notch filter can be adaptively acquired using the results of the elastic surrogate model as a result that the elastic mode is filtered out to ensure system stability for the hypersonic vehicle, and the results with regard to the illustrative example show that the presented control law is valid to suppress the disturbances of the elastic mode and keep system stability.