Ascent guidance for air-breathing hypersonic vehicles based on dynamic inversion

Abstract

This article presents an adaptive ascent guidance strategy based on dynamic inversion for air-breathing hypersonic vehicles. Since dynamic inversion is effective to deal with nonlinear problems, this method is employed here to design the ascent guidance command. Compared with conventional dynamic inversion, an adaptive scheme is added into the guidance strategy, which significantly improves the robust performance of this method. And the stability analysis is given to prove the convergence of the nonlinear system. Besides, considering the real-world dynamic motions of the vehicle should be aware to calculate the guidance command instead of referring to the nominal model, the least squares algorithm with forgetting factor is adopted in this article to estimate the real thrust acceleration. And the reference trajectory generated online by polynomial fitting shows better adaptation to the change of flight environment and less computation time than other trajectory generation methods. Results of simulation are provided to demonstrate the feasibility and effectiveness of this approach.

Introduction

Researches on hypersonic vehicles using air-breathing propulsion have attracted significant interests over the last decades. ^1,2 This kind of vehicles can extremely reduce the flight time to achieve the desired destination due to the high speed compared with other traditional aircrafts. With the development of hypersonic technologies, researchers have achieved great progress to ensure the feasibility and efficiency of flight vehicles. ^{3 –5} Some of them are potential to be implemented in orbital transportation to access the space ⁶ and missile design to attack or intercept the target more precisely. ⁷

The ascent flight of the air-breathing hypersonic vehicle studied in this article is mainly in the dense atmosphere since the altitude of the aircraft is relatively low. Therefore, the resultant of lift, drag force and the engine thrust impacts the ascent longitudinal trajectory, which makes the dynamic models highly nonlinear. Besides, some uncertain parameters and fast changes of the state variables in real flight environment bring more challenges to obtain the online guidance command that needs to achieve the destination and satisfy all the path constraints. So, it becomes practically difficult to solve the ascent guidance problem of hypersonic vehicles and more attempts are necessary to overcome those difficulties.

Researchers have developed many feasible approaches about the ascent guidance problem. The strategies can be generally categorized into three manners according to the information needed to design the guidance command. The first manner is the open-loop guidance way, and this approach is simply following the reference trajectory, which satisfies all the constraints and is generated offline. Apparently, this kind of methods provide less capability of adapting to fast change of flight conditions and are not suitable to be taken to solve the ascent guidance problem. The closed-loop guidance technique that utilized offline optimized trajectory is regarded as the second manner here. Since this closed-loop way owns certain adaptability to the real flight environment, this guidance strategy has been researched and implemented very widely. An adaptive predictive guidance method based on nonlinear programming was presented in the literature, ⁸ and the optimized parameters are performed before launch utilized for online optimization flight. Brinda et al. ⁹ developed an adaptive guidance law by using a feedback loop based on a second-order rate controller and the followed reference optimal trajectory is generated offline. Zhang et al. ¹⁰ proposed a strategy that the guidance command is achieved by online optimal trajectory generation based on pseudo-spectral method, and the offline optimized trajectory is used as an initial value of the pseudo method. The closed-loop guidance scheme that only needs the information measured by navigation systems is the last manner. The advantage of this way is that the preparation time before flight is reduced extremely and the guidance command calculated by the real-time states variables can easily meet all the demands. Lu et al. ¹¹ translated the constrained optimal atmospheric guidance problem to a two-point boundary-value problem by optimal control theory and used the classical finite difference method to solve it. Besides, another research ¹² on ascent guidance was presented, which utilized the primer vector theory to determine the optimal trajectory.

By comparing these three guidance manners, it is noted that the third way qualifies the most adaptability under non-nominal mission conditions for hypersonic vehicles. Besides, in order to further reduce the computation amount and avoid calculating the costate variables based on the optimal control theory, some control methods are combined with the guidance technique to generate a robust guidance command, which has been a trend for guidance research. Vachon et al. ¹³ designed a trajectory tracking guidance law by the application of discrete multi-model predictive control method. Christian and Yuri ¹⁴ developed a new approach that the second-order sliding mode control technique was adopted to get the commanded velocity. Qiao and Chen ¹⁵ presented an adaptive dynamic inversion guidance law to fly along the optimal trajectory for a hypersonic vehicle. Compared with other control methods, the procedure of dynamic inversion is less complicated, and this technique can be easily implemented. Therefore, in this work, the advantage of dynamic inversion is taken to design the ascent guidance command. And the basic definition of this method has been introduced in Lakshmikanth et al. ¹⁶ So, this article concentrates on the implementation of dynamic inversion on ascent guidance problem.

According to the introduction of this dynamic inversion, the principle of this method is to transform the originally nonlinear problem to a linear system. Thus, the system deviation and disturbance during the translation should be taken into account when designing the guidance command. There are many techniques to compensate the uncertainties such as disturbance observer ¹⁷ designed by backstepping method, ¹⁸ neural network ¹⁹ and other methods. ^20,21 Due to the high speed of hypersonic vehicles, a feasible and less complex strategy to estimate the uncertainties is necessary to generate guidance command immediately. So, in this article, a simple adaptive scheme is developed in which an adaptive variable is defined to describe all the uncertainties. And proof of convergence is given in this article to demonstrate the robustness and feasibility of the proposed adaptive scheme. Besides, this article adopts the least square algorithm with forgetting factor to estimate the acceleration of the engine thrust in a real flight condition in order to get accurate guidance command. And the proposed algorithm has been used extensively in filtering, identifying and predicting its good convergence property and well adapting to the measurement noise. ^22,23 Furthermore, since the reference altitude trajectory is required during the employment of the dynamic inversion, an online reference trajectory generation strategy is utilized in which the polynomial with respect to time represents the reference altitude trajectory according to the ascent trend of altitude. This manner not only performs good approximation property but also takes little computation time. Therefore, the designed scheme as the third guidance manner mentioned above only takes advantage of the measurable information to get the guidance command, which exhibits high adaptability.

With the description of the guidance strategy, this article is organized as follows: the problem formulation of the ascent guidance is presented in the second section. The procedure of the designed guidance strategy is exhibited in the third section, which included three parts about application of dynamic inversion, the estimating the acceleration of thrust and online generation of reference trajectory. Finally, some simulation results in the fourth section are provided to explain the feasibility of this method and fifth section is the conclusion.

Formulation of ascent guidance problem

This section will mainly discuss the dynamic model adopted to describe the motion of an air-breathing hypersonic vehicle and the flight constraints that should be satisfied during the design of guidance command.

By assuming that there is no vehicle side-slip or roll, the dynamic equations of a point mass with three degrees of freedom for non-rotating earth are employed to describe the motion of a hypersonic vehicle. And the point mass equations of motion are

$\{\begin{array}{c}\dot{h}=v\text{sin}\mathit{\gamma }\\ \dot{v}=\frac{T\cos \alpha -D}{m}-g\text{sin}\mathit{\gamma }\\ \dot{\mathit{\gamma }}=\frac{T\text{sin}\alpha }{mv}+\frac{L}{mv}-\frac{g}{v}\cos \mathit{\gamma }+\frac{v}{r}\cos \mathit{\gamma }\\ \dot{m}=f(\alpha ,Ma,\varphi ,h)\end{array}$

1

where T, L and D represent the engine thrust, aerodynamic lift and drag, respectively. h, v and γ denote the altitude, airspeed and flight path angle, respectively. $r=h+R_E$ , where R_E is the radius of the earth and g is the gravitational acceleration. Besides, α is the angle of attack and ϕ is the throttle of a specific air-breathing engine. m is defined as the mass of the vehicle. Based on the dynamic model of an air-breathing hypersonic vehicle. ²⁴ , the forces can be described as follows

$\begin{array}{l}T=0.5\rho v^2{s}_A{T}_C\hfill \\ L=0.5\rho v^2{s}_{ref}{C}_L\hfill \\ D=0.5\rho v^2{s}_{ref}{C}_D\hfill \end{array}$

2

where

$\begin{array}{l}{T}_C=T_C\left(Ma,\alpha ,\varphi \right)\hfill \\ C_L=C_L\left(\alpha ,Ma\right)\hfill \\ C_D=C_D\left(\alpha ,Ma\right)\hfill \end{array}$

3

According to the model, ²⁴ the lift coefficient C_L and drag coefficient C_D are functions of Mach number Ma and angle of attack α. The thrust coefficient T_C is influenced by Mach number Ma, angle of attack α and throttle ϕ. And the standard atmosphere model is used here to calculate the air density ρ. s_A is the effective area of the engine and s_ref is the reference area.

Based on the dynamic equations, it is observed that the angle of attack α and the throttle ϕ are control variables to determine the states. Here, in order to simplify this ascent guidance problem, the throttle ϕ is fixed as the maximum value ϕ_max within its bound. Therefore, the angle of attack α is the unique control variable, which is regarded as guidance command to be designed.

The aim to realize the ascent guidance is that the vehicle with designated initial states h_i, Ma_i and γ_i can achieve the required destination as well as satisfy all the path constraints by appropriate guidance command with off-nominal flight conditions. First, the terminal constraints including the altitude, Mach number and flight path angle are as follows

$h_f=h_f^{*},M{a}_f=M{a}_f^{*},{\mathit{\gamma }}_f={\mathit{\gamma }}_f^{*}$

4

Thus, the terminal states of vehicle h_f,Ma_f and γ_f should achieve the required states $h_f^{*},M{a}_f^{*}\text{ and }{\mathit{\gamma }}_f^{*}$ at the final time. Besides, the path constraints including restrictions on dynamic pressure q and angle of attack are given as

$q=0.5\rho v^2\le Q_{\text{max}},{\alpha }_{\text{min}}\le \alpha \le {\alpha }_{\text{max}},\left|\dot{\alpha }\right|\le {\dot{\alpha }}_{\text{max}}$

5

where Q_max is the allowable maximum dynamic pressure, α_max and α_min are the permitted maximum and minimum angle of attack, respectively, and ${\dot{\alpha }}_{\text{max}}$ is the allowable maximum change rate about angle of attack. After formulating the ascent guidance problem, the proposed guidance strategy is aimed to satisfy all the necessary conditions.

Guidance strategy based on dynamic inversion

Since the ascent guidance problem illustrated above is a highly nonlinear problem, as mentioned before, this article plans to adopt the dynamic inversion to deal with it. Considering that there exists disturbances during the application of dynamic inversion, an adaptive method to eliminate those uncertainties is proposed with stability analysis. Besides, since the force conditions on the vehicle should be aware to calculate the value of guidance command according to the expression, the least squares algorithm with forgetting factor is employed to estimate the actual acceleration of engine thrust. Additionally, the reference altitude trajectory to be tracked is better to be generated online, so the scheme of third order polynomial fitting with respect to time is designed to represent the reference altitude. Hence, three parts are included in this section to describe the procedure of proposed guidance method based on dynamic inversion.

Application of dynamic inversion

Based on the theory of basic dynamic inversion method, the nonlinear system with coupled input–output behaviour is fundamental to the control design. By differentiating the output state with respect to time until the input appears, the relation between the input and the output is set up. And the control law is formulated by inverting the relation. So, combined with the ascent guidance problem, the dynamic equations in equation (1) can be similarly viewed as a nonlinear system, where the altitude is defined as the output.

First, the dynamic equations in equation (1) involved with the altitude are given as

$\{\begin{array}{l}\dot{h}=v\text{sin}\mathit{\gamma }\hfill \\ \dot{v}=a_T\cos \alpha -a_D-g\text{sin}\mathit{\gamma }\hfill \\ \dot{\mathit{\gamma }}=\frac{\left(a_T\text{sin}\alpha +a_L-g\cos \mathit{\gamma }\right)}{v}+\left(\frac{v}{r}\right)\cos \mathit{\gamma }\hfill \end{array}$

6

where

$a_T=\frac{T}{m}\text{, }{a}_L=\frac{L}{m}\text{, }{a}_D=\frac{D}{m}$

7

$a_T\text{, }{a}_L\text{ and }{a}_D$ are the accelerations of engine thrust, lift and drag, respectively. According to the principle of dynamic inversion, the explicit relation is found by differentiating the altitude. And we obtain

$\begin{array}{l}\ddot{h}=v\dot{\mathit{\gamma }}\cos \mathit{\gamma }+\dot{v}\text{sin}\mathit{\gamma }\hfill \\ =\left(\frac{\left(a_T\text{sin}\alpha +a_L-g\cos \mathit{\gamma }\right)}{v}+\left(\frac{v}{r}\right)\cos \mathit{\gamma }\right)v\cos \mathit{\gamma }+\left(a_T\cos \alpha -a_D-g\text{sin}\mathit{\gamma }\right)\text{sin}\mathit{\gamma }\hfill \\ =a_T\text{sin}(\alpha +\mathit{\gamma })-g-a_D\text{sin}\mathit{\gamma }+a_L\cos \mathit{\gamma }+\left(\frac{v^2}{r}\right){{\displaystyle \cos }}^2\mathit{\gamma }\hfill \\ \hfill \end{array}$

Here, the control input of this nonlinear system can be defined as $u=\text{sin}(\alpha +\mathit{\gamma })$ . So, the relation is illustrated as follows

$\begin{array}{l}\ddot{h}=au+b\hfill \\ a=a_T\text{, }b=-g-a_D\text{sin}\mathit{\gamma }+a_L\cos \mathit{\gamma }+\left(\frac{v^2}{r}\right){{\displaystyle \text{cos}}}^2\mathit{\gamma }\hfill \\ u=\text{sin}(\alpha +\mathit{\gamma })\hfill \end{array}$

8

Therefore, the nonlinear system is translated into a linear system and the control input can be formulated. Furthermore, define δ = au + b as the pseudo input and the linear input–output feature is set up

$\ddot{h}=\mathit{\delta }$

9

Note that the accurate model is the key to get the ideal control input by the dynamic inversion method. And the model mismatch may cause disability to accomplish tracking or stabilizing control tasks. To eliminate the drawbacks of this method, some literatures introduce robust control strategies with dynamic inversion, which improves the performance significantly. ^25,26 Referring to the scheme of the adaptive robust control method, the adaptive variable Δ is defined in this article to cover all the uncertainties during the translation to a linear system. And we obtain

$\ddot{h}=\mathit{\delta }+\text{Δ}$

10

Then let h_c denote the reference altitude and define $\mathit{\tau }=\mathit{\delta }-{\ddot{h}}_c$ , the above expression is translated into

$\ddot{h}-{\ddot{h}}_c=\mathit{\tau }+\text{Δ}$

11

Defining the tracking error $e=h-h_c$ , the second-order error dynamics are given by

$\ddot{e}=\mathit{\tau }+\text{Δ}$

12

Since it is difficult to obtain the second-order error used to get the control law τ, by setting $\eta =\dot{e}+\mu e$ , this second order system is turned into a first-order system as follows

$\dot{\eta }=\mathit{\tau }+H\left(e,\dot{e}\right)+\text{Δ}$

13

where μ is a positive scalar and $H\left(e,\dot{e}\right)=\mu \dot{e}$ . And a simple linear controller can be used to deal with this first-order dynamic system. Besides, an adaptive strategy is adopted to eliminate the impact of the uncertain factor Δ. The control law τ can be designed as

$\mathit{\tau }=-\mathit{\lambda }\eta -H\left(e,\dot{e}\right)-{\text{Δ}}^{*}\text{sgn}\left(\eta \right)\left|\text{Δ}\right|\le {\text{Δ}}^{*}$

14

where Δ* is defined as the absolute maximum value of Δ. Define $\widehat{\text{Δ}}$ as the estimate of Δ* and an adaptive method is employed to get the value of $\widehat{\text{Δ}}$ according to η

$\dot{\widehat{\text{Δ}}}=\mathit{\kappa }\left|\eta \right|,\mathit{\kappa }>0$

15

λ and κ are positive scalars. So, the final control law τ is given as

$\mathit{\tau }=-\mathit{\lambda }\eta -H\left(e,\dot{e}\right)-\widehat{\text{Δ}}\text{sgn}\left(\eta \right)$

16

Substituting the above equation into equation (13), we obtain

$\dot{\eta }=-\mathit{\lambda }\eta -\widehat{\text{Δ}}\text{sgn}\left(\eta \right)+\text{Δ}$

17

Then the stability of the dynamics in equation (12) is proved by combining with proposed control law in equation (16) and the adaptive method. Based on Lyapunov method for stability analysis and combined with the relation $\eta =\dot{e}+\mu e$ , firstly define the positive definite Lyapunov function as

$V=\frac{1}{2}{\eta }^2+\frac{1}{2\mathit{\kappa }}{\widetilde{\text{Δ}}}^2$

18

where $\widetilde{\text{Δ}}={\text{Δ}}^{*}-\widehat{\text{Δ}}$ . Differentiating the function V with respect to time, we obtain

$\begin{array}{l}\dot{V}=\eta \cdot \dot{\eta }+\frac{1}{\mathit{\kappa }}\widetilde{\text{Δ}}\cdot \left({\dot{\text{Δ}}}^{*}-\dot{\widehat{\text{Δ}}}\right)\hfill \\ =\eta \cdot \left(-\mathit{\lambda }\eta -\widehat{\text{Δ}}\text{sgn}\left(\eta \right)+\text{Δ}\right)+\frac{1}{\mathit{\kappa }}\left({\text{Δ}}^{*}-\widehat{\text{Δ}}\right)\cdot \left(-\mathit{\kappa }\left|\eta \right|\right)\hfill \\ =-\mathit{\lambda }{\eta }^2-\widehat{\text{Δ}}\cdot \left|\eta \right|+\eta \cdot \text{Δ}-\left|\eta \right|\cdot {\text{Δ}}^{*}+\left|\eta \right|\cdot \widehat{\text{Δ}}\hfill \\ =-\mathit{\lambda }{\eta }^2+\eta \cdot \text{Δ}-\left|\eta \right|\cdot {\text{Δ}}^{*}\hfill \\ \le -\mathit{\lambda }{\eta }^2+\left|\eta \right|\cdot \text{Δ}-\left|\eta \right|\cdot {\text{Δ}}^{*}\hfill \\ \le -\mathit{\lambda }{\eta }^2\hfill \end{array}$

Based on equation (14), there results in $\dot{V}\le 0$ which proves that the system using the proposed control law and adaptive strategy scheme is asymptotically stable and the tracking error will converge to zero over finite time. Therefore, with the relation $\mathit{\tau }=\mathit{\delta }-{\ddot{h}}_c$ , δ can be formulated as

$\mathit{\delta }=-\mathit{\lambda }\eta -H\left(e,\dot{e}\right)-\widehat{\text{Δ}}\text{sgn}\left(\eta \right)+{\ddot{h}}_c$

19

And combined with δ = au + b and the expression of u, the angle of attack α as the guidance command in this article is given by

$\alpha =\arcsin \left(\frac{\mathit{\delta }-b}{a}\right)-\mathit{\gamma }$

20

Since the value of α + γ varies in the interval $\left(-\pi /2,\pi /2\right)$ during the ascent flight, this equation can be calculated to a unique solution of guidance command. Additionally, note that the sign function sgn(η) in equation (19) may cause chattering effect on the trajectory of angle of attack. Generally, the saturation function can be used to replace the sign function, which greatly improves the performance.

After deriving the specific expression of guidance command, it is worth noting that there are some unknown parameters of a and b in equation (8). The next part will propose an efficient method to identify the unknown information.

Estimation of thrust acceleration

Based on the above derivation of guidance command, it is crucial to get the values of accelerations $a_T\text{, }{a}_L\text{ and }{a}_D$ in real flight condition for computation. Due to fast change of the flight conditions for hypersonic vehicles, the dynamic model is significantly mismatched with the nominal one. So, it is necessary to get the actual values of the accelerations instead of referring to nominal dynamic model. Therefore, some identification methods are considered to estimate the accelerations. Since identifying all the values of $a_T\text{, }{a}_L\text{ and }{a}_D$ are difficult with available methods and considering that the engine thrust is the primary factor to determine the shape of ascent trajectory, this article suggests that the acceleration of thrust a_T is identified by the least squares algorithm with forgetting factor while the lift and drag accelerations $a_L\text{ and }{a}_D$ are obtained with nominal dynamic model.

As mentioned above, the least squares algorithm with forgetting factor has been applied extensively for its good convergence property and less computation time compared with other methods. Combining the principle of this method and the dynamic equations containing the thrust acceleration, we obtain

$\begin{array}{l}\dot{v}=a_T\cos \alpha -a_D-g\text{sin}\mathit{\gamma }\hfill \\ \dot{\mathit{\gamma }}=\frac{\left(a_T\text{sin}\alpha +a_L-g\cos \mathit{\gamma }\right)}{v}+\left(\frac{v}{r}\right)\cos \mathit{\gamma }\hfill \end{array}$

21

And the above equations can be obtained as

$\dot{Y}=F\left(Y\right)K+G\left(Y\right)$

22

where

$\begin{array}{l}Y={\left[v,\mathit{\gamma }\right]}^T\text{,}K={\left[a_T\cos \alpha ,a_T\text{sin}\alpha \right]}^T\\ G(Y)={\left[-a_D-g\text{sin}\mathit{\gamma },\frac{\left(a_L-g\cos \mathit{\gamma }\right)}{v}+\left(\frac{v}{r}\right)\cos \mathit{\gamma }\right]}^T\\ F(Y)=\left[\begin{array}{cc}1& 0\\ 0& \frac{1}{v}\end{array}\right]\end{array}$

Note that the derivative term in equation (22) is not measurable directly. So, a first-order low-pass filter is introduced here and the output of the filter Y_f is used to approximate the value Y. And the filter is designed as

$Y\cdot \frac{\mathit{\lambda }}{s+\mathit{\lambda }}=Y_f\Rightarrow {\dot{Y}}_f+{\mathit{\lambda }}_f{Y}_f={\mathit{\lambda }}_{f}Y$

23

where ${\mathit{\lambda }}_f=\text{diag}\left({\mathit{\lambda }}_v,{\mathit{\lambda }}_{\mathit{\gamma }}\right)$ is the coefficient matrix. By using Y_f to replace Y in equation (22), we obtain

${\mathit{\lambda }}_f\left(Y-Y_f\right)=F\left(Y_f\right)K+G\left(Y_f\right)$

24

So, the time-varying problem is translated into a static system and by identifying the vector K, the estimate of thrust acceleration is achieved. Furthermore, the above equation can be regarded as the following linear form

$W=X\mathit{\theta }+\epsilon$

25

Connected with the original least squares algorithm, W represents the desired output vector, X denotes the input vector, θ is the parameter vector to be estimated and ε is the measurement noise, which is assumed as a white Gaussian noise. Additionally, let ${\stackrel{⌢}{\mathit{\theta }}}_L$ as the estimate of θ. So, the least squares cost function is defined as

$J=(W-X{\stackrel{⌢}{\mathit{\theta }}}_L{)}^T{B}^{-1}(W-X{\stackrel{⌢}{\mathit{\theta }}}_L)$

26

where B is the variance matrix of the noise ε which is a diagonal matrix. Based on the necessary conditions for minimization of cost function, the formulation of ${\stackrel{⌢}{\mathit{\theta }}}_L$ is solved as

${\stackrel{⌢}{\mathit{\theta }}}_L={\left(X^T{B}^{-1}X\right)}^{-1}{X}^T{B}^{-1}W$

27

Then discrete the time domain and let ${\stackrel{⌢}{\mathit{\theta }}}_{L,k}$ as the estimate of θ for the kth time instant while ${\stackrel{⌢}{\mathit{\theta }}}_{L,k+1}$ as the estimate for the (k + 1) th time instant. So, $X_k\text{, }{X}_{k+1}\text{, }{W}_k\text{ and }{W}_{k+1}$ are defined in the similar way. And the expressions of ${\stackrel{⌢}{\mathit{\theta }}}_{L,k}$ and ${\stackrel{⌢}{\mathit{\theta }}}_{L,k+1}$ are

$\begin{array}{l}{\stackrel{⌢}{\mathit{\theta }}}_{L,k}={\left(X_k^T{B}^{-1}{X}_k\right)}^{-1}{X}_k^T{B}^{-1}{W}_k\hfill \\ {\stackrel{⌢}{\mathit{\theta }}}_{L,k+1}={\left(X_{k+1}^T{B}^{-1}{X}_{k+1}\right)}^{-1}{X}_{k+1}^T{B}^{-1}{W}_{k+1}\hfill \end{array}$

28

And ζ is denoted as the forgetting factor, which is a positive constant here. Referred to the least squares algorithm with forgetting factor, the presentations of X_k+1 and W_k+1 are

$X_{k+1}=\left[\begin{array}{c}\zeta X_k\\ x_{k+1}\end{array}\right]\text{, }{W}_{k+1}=\left[\begin{array}{c}\zeta W_k\\ w_{k+1}\end{array}\right]$

29

where x_k+1 and w_k+1 are the vectors representing additional information of vector X and vector W gained at the (k + 1) th time instant, respectively. Combining equations (28) and (29), the updated equation of ${\stackrel{⌢}{\mathit{\theta }}}_{L,k}$ is given as

$\begin{array}{l}{U}_k=P_k{x}_{k+1}^T{\left({\zeta }^{2}B+x_{k+1}{P}_k{x}_{k+1}^T\right)}^{-1}\hfill \\ P_{k+1}=\frac{1}{{\zeta }^2}{\left(P_k-U_k{x}_{k+1}P\right)}_k\hfill \\ {\stackrel{⌢}{\mathit{\theta }}}_{k+1}={\stackrel{⌢}{\mathit{\theta }}}_k+U_k\left(w_{k+1}-x_{k+1}{\stackrel{⌢}{\mathit{\theta }}}_k\right)\hfill \end{array}$

30

where $P_k={\left(X_k^T{B}^{-1}{X}_k\right)}^{-1}$ . With this updated scheme, the estimate ${\stackrel{⌢}{\mathit{\theta }}}_{L,k}\left(k=1\text{, }2\text{, }3\text{,}\dots \right)$ of θ at each time instant can be achieved. In order to estimate the parameter K in system (24), the specific forms of W, X and θ in the least squares algorithm are given as

$\begin{array}{l}W={\mathit{\lambda }}_f\left(Y-Y_f\right)-G\left(Y_f\right)\hfill \\ X=F\left(Y_f\right)\hfill \\ \mathit{\theta }=K\hfill \end{array}$

31

Besides, the expressions of Y, G, F and K are expressed by the dynamic states in equation (21) according to equation (22).

Consequently, the coefficient vector K can be estimated by the proposed least squares algorithm with forgetting factor, and then the thrust acceleration is able to be estimated. Here, let $\stackrel{⌢}{K}$ as the estimate of K, and we obtain

$\begin{array}{l}K\approx \stackrel{⌢}{K}={\left[{\stackrel{⌢}{K}}_1,{\stackrel{⌢}{K}}_2\right]}^T\hfill \\ \sqrt{{\stackrel{⌢}{K}}^2{}_1+{\stackrel{⌢}{K}}^2{}_2}\approx \sqrt{{\left(a_T\cdot \cos \alpha \right)}^2+{\left(a_T\cdot \text{sin}\alpha \right)}^2}=a_T\hfill \end{array}$

32

Finally, the estimate of the thrust acceleration a_T is obtained. As mentioned previously, the aerodynamic accelerations a_L and a_D depend on the real-time angle of attack α and state variables taking nominal dynamic model. Since the state variables are measurable, the accelerations can be viewed as the function of α. Furthermore, the estimate of thrust acceleration a_T associated with a_L, a_D can be also viewed as the function of α. Considering the equation δ = au + b and the expressions in equation (8), it is easy to find that both right terms au and b of this equation that rely on α are described as

$f\left(\alpha \right)=au+b-\mathit{\delta }=0$

33

In this article, the Newton iterative method is employed to get the command guidance α. And the next part will generate the online altitude reference trajectory to obtain the tracking error e in the expression δ.

Online reference trajectory generation

Due to the characters of hypersonic vehicles, the offline optimized trajectory as the reference trajectory to be tracked has little capability of adapting to the fast changing flight conditions. Besides, the first-order and second-order derivatives of reference altitude with respect to time should be aware in order to get the guidance command. Thus, an online reference altitude trajectory with respect to time is necessary. In this article, considering that altitude of the vehicle is almost increasing during the ascent phase, the scheme of polynomial fitting to generate the reference altitude trajectory is appropriately used here. Eventually, a third-order polynomial function with respect to time is used to generate the reference altitude trajectory

$h_c\left(t\right)=a_3{t}^3+a_2{t}^2+a_{1}t+a_0$

34

where $a_0\text{, }{a}_1\text{, }{a}_2\text{ and }{a}_3$ are the coefficients and h_c(t) represents the reference altitude. The state point where the reference trajectory is being generated is viewed as the starting point including the state variables $h_0\text{, }{v}_0\text{ and }{\mathit{\gamma }}_0$ and time t₀. And the desired terminal states $h_f^{*}\text{, }M{a}_f^{*}\text{ and }{\mathit{\gamma }}_f^{*}$ and final time t_f constitute the final point. Combined with the information of starting point and final point, the third order polynomial is determined as

$\begin{array}{l}{a}_3{t}_0^3+a_2{t}_0^2+a_1{t}_0+a_0=h_0\hfill \\ 3a_3{t}_0^2+2a_2{t}_0+a_1=v_0\text{sin}{\mathit{\gamma }}_0\hfill \\ a_3{t}_f^3+a_2{t}_f^2+a_1{t}_f+a_0=h_f^{*}\hfill \\ 3a_3{t}_f^2+2a_2{t}_f+a_1=v_f^{*}\text{sin}{\mathit{\gamma }}_f^{*}\hfill \end{array}$

35

where the desired terminal velocity $v_f^{*}$ is obtained by $h_f^{*}$ and $M{a}_f^{*}$ . And the coefficients are expressed as

$\begin{array}{l}{a}_0=h_0-a_3{t}_0^3-a_2{t}_0^2-a_1{t}_0\hfill \\ a_1=\frac{h_f^{*}-h_0}{t_f-t_0}-a_2\left(t_f+t_0\right)-a_3\left(t_f^2+t_0{t}_f+t_0^2\right)\hfill \\ a_2=\frac{v_f^{*}\text{sin}{\mathit{\gamma }}_f^{*}-v_0\text{sin}{\mathit{\gamma }}_0-3a_3\left(t_f^2-t_0^2\right)}{2\left(t_f-t_0\right)}\hfill \\ a_3=\frac{\left(t_f-t_0\right)\cdot \left(v_f^{*}\text{sin}{\mathit{\gamma }}_f^{*}+v_0\text{sin}{\mathit{\gamma }}_0\right)-2\left(h_f^{*}-h_0\right)}{{\left(t_f-t_0\right)}^3}\hfill \end{array}$

36

Thus, the first-order and second-order derivatives of reference altitude with respect to time are described as

$\begin{array}{l}{h}_c\left(t\right)=a_3{t}^3+a_2{t}^2+a_{1}t+a_0\hfill \\ {\dot{h}}_c\left(t\right)=3a_3{t}^2+2a_{2}t+a_1\hfill \\ {\ddot{h}}_c\left(t\right)=6a_{3}t+2a_2\hfill \end{array}$

37

Note that the final time t_f should be predetermined to calculate the values of third-order polynomial coefficients. However, the final time t_f is not fixed and given during the ascent phase. In order to achieve the required Mach number $M{a}_f^{*}$ , the similar scheme of polynomial fitting applied on the velocity trajectory is taken to estimate the final time t_f since the velocity is always increasing during the ascent flight. And a second order polynomial with respect to time is applied to approximate the velocity

$v\left(t\right)=b_2{t}^2+b_{1}t+b_0$

38

where b ₀, b ₁ and b ₂ are the coefficients. Here, the initial velocity v_i obtained by h_i and Ma_i as well as the present point where the final time is being estimated including the velocity v_e, flight path angle γ_e, angle of attack α_e, dynamic accelerations $a_{T_e}\text{ and }{a}_{D_e}$ and the time t_e are used to determine this second-order polynomial

$\begin{array}{l}{b}_0=v_i\hfill \\ 2b_2{t}_e+b_1=a_{T_e}\cos {\alpha }_e-a_{D_e}-g\text{sin}{\mathit{\gamma }}_e\hfill \\ b_2{t}_e^2+b_1{t}_e+b_0=v_e\hfill \end{array}$

39

So, the coefficients are

$\begin{array}{l}{b}_0=v_i\hfill \\ b_1=a_{T_e}\cos {\alpha }_e-a_{D_e}-g\text{sin}{\mathit{\gamma }}_e-2b_2{t}_e\hfill \\ b_2=\frac{t_e\left(a_{T_e}\cos {\alpha }_e-a_{D_e}-g\text{sin}{\mathit{\gamma }}_e\right)-\left(v_e-v_i\right)}{t_e^2}\hfill \end{array}$

40

Therefore, the terminal velocity can be expressed by the given polynomial with the final time

$v_f^{*}=b_2{t}_f^2+b_1{t}_f+b_0$

41

And the final time t_f is obtained as

$t_f=\frac{-b_1+\sqrt{b_1{}^2-4b_2\left(b_0-v_f^{*}\right)}}{2b_2}$

42

It is worth noting that at the beginning of the ascent flight, the final time should be estimated for reference trajectory generation. The above method may not be proper to get the estimation. Here, a rough way to estimate the final time is taken at the beginning. Combined with the dynamic motion equation of the altitude, we obtain

$\begin{array}{l}\frac{dh}{dt}=v\text{sin}\mathit{\gamma }\Rightarrow dh\approx \overline{v}\text{sin}\overline{\mathit{\gamma }}dt\hfill \\ dh\approx \left(h_f^{*}-h_i\right)\approx \overline{v}\text{sin}\overline{\mathit{\gamma }}\cdot \left(t_f-0\right)\hfill \end{array}$

43

where $\overline{v}=\left(v_i+v_f^{*}\right)/2$ and $\overline{\mathit{\gamma }}{=5}^{\circ }$ since the flight path angle γ varies mainly within the region [0,10°] during the ascent phase. So, the final time is obtained by the above equation at the beginning. With the flight time increasing, the proposed polynomial fitting method is adopted to estimate the final time more precisely.

Apparently, the estimation of the final time is along with the generation of reference altitude trajectory. And the mechanism to active the conduct to generate the trajectory is based on the difference Δ _h between the actual altitude h and the reference altitude h_c. If

$\text{Δ}h=\left|h-h_c\right|>100\text{\hspace{0.17em}m}$

44

the reference trajectory is started to be generated once again. Besides, considering the maximum dynamic pressure constraint, when the real dynamic pressure q is close to Q_max, the reference trajectory is refreshed by adding a small positive constant σ in equation (34)

$3a_3{t}_0^2+2a_2{t}_0+a_1=v_0\text{sin}\left({\mathit{\gamma }}_0+\sigma \right)$

45

This scheme is aimed to increase the altitude of vehicle quickly, which can further reduce the dynamic pressure q to satisfy the path constraint. Consequently, using the proposed guidance strategy based on dynamic inversion, the guidance command α is achieved by equation (33).

Numerical example

This section will present the results of simulation with five cases to test the guidance method implemented on the vehicle with air-breathing hypersonic propulsion during the ascent phase.

Conditions for simulation

The specified values of the constraints including the initial and terminal states as well as the path constraints are illustrated as follows. And it is worth mentioned that some flight parameters for simulation are different from the model recommended in ‘Formulation of ascent guidance problem’ section. Firstly, the initial and desirable terminal states are shown in Table 1.

Table 1.

The initial and final conditions.

	Altitude (km)	Mach number	Flight path angle (°)
Initial conditions	25	5	0
Final conditions	35	5.9	0

And the values of path constraints are

$-1°\le \alpha \le 5°,\left|\dot{\alpha }\right|\le 3°/s\text{, \hspace{0.17em}}{Q}_{\text{max}}=50\text{\hspace{0.17em}kPa}$

Besides, fives cases for simulations are as follows:

S0: nominal model;

S1: 10% error of the thrust coefficient;

S2: 10% error of the thrust coefficient, lift coefficient and drag coefficient, respectively;

S3: −10% error of the thrust coefficient;

S4: −10% error of the thrust coefficient lift coefficient and drag coefficient, respectively.

Performance of numerical results

Firstly, the comparison of the estimate of the thrust acceleration by the least squares algorithm with forgetting factor method with the actual value under nominal dynamic conditions is illustrated in Figure 1. As shown in Figure 1, the estimated value is quite different from the actual value at the beginning owing to less information aware of the dynamic nonlinear system. With the increasing knowledge about this nonlinear system, the estimate of thrust acceleration becomes more and more close to the real value. Hence, the results demonstrate that the least squares algorithm with forgetting factor method is efficient and feasible to identify the actual thrust acceleration.

Figure 1.

The performance of least squares algorithm with forgetting factor method.

Then the numerical results with five cases are presented to exhibit the performance of implemented guidance method. The trajectories of the altitude, velocity, flight path angle and dynamic pressure of the vehicle are shown in Figures 2 to 5.

Figure 2.

The trajectory of altitude.

Figure 3.

The dynamic response of velocity.

Figure 4.

The trajectory of flight path angle.

Figure 5.

The dynamic pressure.

From Figures 2 to 4, it is noted that the vehicle is gradually approaching to the ideal destination, which explains that the scheme about polynomial fitting used for the altitude and velocity shows great approximation property of the real dynamic response. And the results also indicate that the designed guidance command based on dynamic inversion provides significant capability of tracking the reference trajectory.

The allowable terminal states deviations from the desired values for altitude, Mach number and flight path angle are set as $\text{Δ}{h}_f=10\text{\hspace{0.17em}m, Δ}M{a}_f=0.1\text{ and Δ}\mathit{\gamma }=0.1°$ , respectively. The performance of the terminal states with the five cases is shown in Table 2.

Table 2.

The performance of the terminal states with the five cases.

Cases	Altitude (km)	Mach number	Flight path angle (°)
S0	35.000	5.90	0.00
S1	35.002	5.95	0.07
S2	35.000	5.97	0.00
S3	35.000	5.87	0.00
S4	35.000	5.88	0.00

It is observed that the actual deviations of the terminal states are all satisfactory. Since the final flight path angle is required to 0°, combined with the equation (34), the terminal required Mach number is achieved basically dependent on the scheme of the second-order polynomial fitting method applied for approximating the trajectory of velocity. Referring to Figure 3, the dynamic response of velocity shows that this estimation method is valid to predict the flight time according to the terminal constraint of Mach number. Besides, the restriction on the dynamic pressure is satisfied as shown in Figure 5.

The Monte Carlo simulation is implemented for running 500 times simulation which the errors of the thrust coefficient, lift coefficient and drag coefficient, respectively, and they randomly vary within the range of −10% to 10%. The results of the terminal states for each simulation are displayed in Figure 6.

Figure 6.

The distribution of terminal states for 500 times simulation.

As shown in Figure 6, combined with the allowable terminal states’ deviations, the cases satisfying the constraints are located in the described rectangle. And the statistics of the 500 times simulation show that the probability to get satisfactory results under off-nominal flight conditions is more than 93%. Consequently, all displayed simulation results demonstrate that the proposed ascent guidance method is significantly feasible and efficient for application.

Conclusion

A detailed solution of ascent guidance based on dynamic inversion for air-breathing hypersonic vehicles has been presented in this article. And the schemes about estimation of the actual thrust acceleration and the online reference trajectory generation are proposed to deal with the dynamic model mismatch caused by changing flight conditions of hypersonic vehicles. Additionally, the dynamic inversion with adaptive strategy shows good capability of tracking the reference trajectory.

It is worth noting that this article achieves several contributions for ascent guidance. First, the proposed dynamic inversion method is successfully implemented on hypersonic vehicles. Second, an adaptive algorithm to eliminate the uncertainties of the dynamic model is developed and the dynamic system is proved as asymptotically stable by Lyapunov stability theory. Besides, the least squares algorithm with forgetting factor method to identify the real acceleration of thrust is applied to get the accurate result with off-nominal conditions. Furthermore, the online reference trajectory generated by polynomial shows significant adaptation to the high speed of hypersonic vehicles. The procedure to get the guidance command is totally adaptive and automatic without manual adjustment, and the result of Monte Carlo simulation for 500 times simulation demonstrates the feasibility of the approach based on dynamic inversion. Thus, the proposed ascent guidance method offers great potential for practical applications.

Besides, lots of improvements about this method will be carried out in the future. For example, the difference between the real aerodynamic accelerations and the nominal values can be viewed as a specific uncertain factor distinguished with other disturbances such as measurement noises and similar adaptive strategies are available to indicate this uncertainty. Moreover, considering the time delay of real-world sensors, some steps should be taken to eliminate it.