1. INTRODUCTION
The air-breathing hypersonic launch vehicle has many inherent features to reduce the cost of payloads for future space transportation (Henry and McLellan, Reference Henry and McLellan1971). Optimal steering for an air-breathing hypersonic launch vehicle exhibits completely different behaviour from conventional rocket flight. Both thrust and lift depend on the atmospheric dynamic pressure (Yamamoto and Kawaguchi, Reference Yamamoto and Kawaguchi2007), and an optimal ascent guidance algorithm should be adopted during the flight to achieve a fuel-optimal control trajectory.
Ascent phase trajectory optimisation and guidance of an air-breathing hypersonic vehicle is a difficult non-linear problem with various constraints (Vinh, Reference Vinh1981; García, Reference García2005). Aerodynamics, propulsion, vehicle and attitude are strongly coupled during the ascent phase. The traditional approach for trajectory guidance is online tracking of the optimal reference profile which was designed offline (Li et al., Reference Li, Jiang and Liu2014; Saraf et al., Reference Saraf, Leavitt, Chen and Mease2004), but this method does not have enough capability to deal with external disturbances such as aerodynamic and thrust errors.
When solving optimal guidance trajectories in endo-atmospheric flight, it is difficult to solve the Two Point Boundary Value Problem (TPBVP) in time. Both indirect and direct methods have been mostly used to solve the offline trajectory optimisation problem in previous literature (Conway, Reference Conway2012). The numerical algorithms of non-linear trajectory optimisation for vehicles were summarised and systematically analysed by Betts (Reference Betts1998) and Huang et al. (Reference Huang, Lu and Nan2012). Dukeman (Reference Dukeman2002) and Dukeman and Calise (Reference Dukeman and Calise2003) developed a closed-form ascent guidance algorithm, to cyclically solve the TBBVP during endo-atmospheric flight. The work of Murillo and Lu (Reference Murillo and Lu2010) adopted an indirect scheme to achieve a proposed closed-loop ascent guidance algorithm by solving TPBVP; they first derived the necessary conditions of the optimal ascent trajectory. However, indirect methods suffer from difficulties in finding an appropriate initial guess and obtaining a convergent solution for the TPBVP.
In recent years, with improvements in on board computer performance, a direct optimal method named Gauss Pseudo-spectral Method (GPM) has been proved to be an effective way to solve the real-time trajectory optimisation problem (Bollino et al., Reference Bollino, Ross and Doman2006). GPM not only has high precision and fast convergence in solving the optimal control problems, but also proved the equivalence between the Karush-Kuhn-Tucker (KKT) conditions and the Hamiltonian Boundary Value Problem (HBVP) first-order optimality conditions (Benson, Reference Benson2005; Huntington, Reference Huntington2007). This removes the need for computing analytical gradients of aerodynamic coefficients. This approach to optimal feedback control based on the Carathodory-π trajectory concept is capable of solving industrial-strength problems and it offers the ability to solve the Non-linear Programming (NLP) problem online (Ross et al., Reference Ross, Sekhavat, Fleming and Gong2008). Two algorithms (free sampling frequency and fixed sampling frequency) are proposed in a real-time optimal feedback control method based on a domain transformation technique and a Radau-based Pseudo-spectral method. The techniques are used for a flexible robot arm and a benchmark inverted pendulum problem (Ross et al., Reference Ross, Gong, Fahroo and Kang2006) as well. A rapid trajectory optimisation method via GPM was recently studied for air-breathing hypersonic vehicles (Zhang et al., Reference Zhang, Lu, Liu and Wang2012). However, a significant barrier to this optimisation-based control strategy is computational delay, especially when the solution time of the NLP problem exceeds the sampling time.
Recently proposed non-linear programming sensitivity methods have reduced online computational costs and can lead to significantly improved performance. Sensitivity update strategy has been successfully used for Non-linear Model Predictive Control (NMPC) problems (Zavala and Biegler, Reference Zavala and Biegler2009; Yang and Biegle Reference Yang and Biegler2013). NMPC is an approach to feedback design that is based on the solution of an Optimal Control Problem (OCP) at each controller update step (Grüne and Pannek, Reference Grüne and Pannek2011). The basic sensitivity strategy for NLP problem solvers is derived through application of the Implicit Function Theorem (IFT) to the KKT conditions of a parametric NLP problem (Fiacco, Reference Fiacco1976). However, the air-breathing hypersonic launch vehicle exhibits severe non-linearity, the sample-and-hold technique used in traditional NMPC problems cannot be adopted.
In this paper, an online guidance algorithm for an air-breathing hypersonic vehicles with process constraints and terminal states constraints is proposed based on the Pseudo-spectral method and sensitivity analysis. An approximated closed-loop optimal control problem solution is computed in the following way: in each sampling interval, the nominal non-linear model was used to predict the next sampling states and an optimal control problem is solved based on the predicted states. When the system dynamics are simulated to the next sampling instant, based on the difference between predicted states and measured states, the resulting optimal control sequence which is applied as input for the next sampling interval(s), is corrected via sensitivities. This procedure is then repeated iteratively. The proposed ascent optimal guidance method can successively generate online open-loop suboptimal controls without the design procedure of an inner-loop feedback controller. This new proposed guidance method has the ability of online trajectory reconstruction for the optimal performance index and has high accuracy and robustness via sensitivity updates.
This paper is organised as follows. In Section 2 the nominal differential equations of an air-breathing hypersonic launch vehicle are presented. In Section 3 the online optimal guidance algorithm based on the Pseudo-spectral method and sensitivity updates is proposed. In Section 4, considering model parameters uncertainties and external disturbance, a sampling theorem is proposed which indicates the effect of the Lipschitz constant of the dynamics on sampling frequency. In Section 5, simulation results indicate that the method proposed above offers improved performance and has promising ability for online computation. In Section 6, some conclusions and indications of future research are given.
2. PROBLEM FORMULATION
During the ascent phase for a point-mass vehicle model, the hypersonic vehicle dynamics are modelled using Equations (1) to (4) for a non-rotating Earth assumption (Prasanna et al., Reference Prasanna, Ghose, Bhat, Bhattacharyya and Umakant2005).
$$\displaystyle{{dh} \over {dt}} = v\sin\gamma $$
$$\displaystyle{{dv} \over {dt}} = \displaystyle{{T\cos\alpha - D} \over m} - \displaystyle{{\mu \sin\gamma} \over {(R_e + h)^2}} $$
$$\displaystyle{{d\gamma} \over {dt}} = \displaystyle{{T\sin\alpha + L} \over {mv}} + \cos\gamma \left[ {\displaystyle{v \over {(R_e + h)}} - \displaystyle{\mu \over {v(R_e + h)^2}}} \right]$$
$$\displaystyle{{dm} \over {dt}} = - \displaystyle{T \over {I_{sp} g}}$$where, h is the altitude, v is the velocity, α is the angle of attack, γ is the flight path angle, m is the mass, T is the aircraft thrust, I sp is the engine specific, D and L are drag and lift, respectively. μ is the gravitational constant, g is the gravitational acceleration and R e is the radius of the earth.
Lift force L and drag force D are given as follows:
$$L = \displaystyle{1 \over 2}\rho v^2 S_{ref} C_L (M_a, \alpha )$$
$$D = \displaystyle{1 \over 2}\rho v^2 S_{ref} C_D (M_a, \alpha )$$where ρ, S ref are the air density of the current altitude and the reference area respectively. C L and C D are the lift and drag coefficients respectively, which are the non-linear functions of the attack angle α and the mach M a.
Considering the thrust T and aerodynamic forces as a composition of forces, we define the longitudinal, lateral accelerations and dynamic pressure as a long, alat and Q respectively.
$$a_{long} = \displaystyle{{T\cos\alpha - D} \over m} - \displaystyle{{\mu \sin\gamma} \over {(R_e + h)^2}} $$
$$a_{lat} = \displaystyle{{T\sin\alpha + L} \over m} + \cos\gamma \left[ {\displaystyle{{v^2} \over {(R_e + h)}} - \displaystyle{\mu \over {(R_e + h)^2}}} \right]$$
$$Q = \displaystyle{1 \over 2}\rho v^2 $$The 1976 US Standard Atmosphere is used as the model for the density ρ which is applied to compute lift and drag force. The speed of sound c used to calculate the Mach number M a is also provided in the model.
Reduction in fuel will make the vehicle climb to a predetermined height; if we put minimum fuel consumption as a objective of ascent flight phase, the cruise segment of flight range will be greatly increased. For these reasons, the optimisation objective is to minimise fuel consumption (or maximise mass at the end) during the ascent phase, that is equivalent to the following formula:
$$\min J = - m_f $$3. OPTIMAL GUIDANCE ALGORITHM BASED ON PSEUDO-SPECTRAL METHOD AND SENSITIVITY UPDATES
The idea of a recursive open-loop solution is compelling, because control is available at any point that has been defined as optimal. Once the state is perturbed from the expected optimal path, correcting back to that trajectory is likely not optimal from the disturbed position, and a new optimal path originating from the current state should be introduced. However, a significant barrier to this optimisation-based control strategy is computational delay, especially when the solution time of the optimal control problem exceeds the sampling time. Sensitivity analysis can be adopted to reduce online computational costs.
3.1. Optimal control problem
The optimal control problem, which is often turned into a Non-linear Programming (NLP) problem, can be written in the following general form (the 1976 US atmosphere model which as a lookup table was implemented in the optimisation problem (Murillo and Lu, Reference Murillo and Lu2010)). The state, x(t) = [h v γ m], control, u(t) = α, and final time, t f, that minimise the cost function
$$J = M[x(t_f ),t_f ] + \int_{t_0} ^{t_f} g(x(t),u(t),t)dt$$subject to the constraints, differential dynamics, terminal constraints and path constraints:
$$\dot x(t) = f\left[ {x(t),u(t),t} \right]$$
$$\varphi \left[ {x(t_0 ),t_0, x(t_f ),t_f} \right] = 0$$
$$c[x(t),u(t),t] \le 0$$3.2. The basic principle of the Gauss Pseudo-spectral Method
The trajectory optimisation problem is formulated over the time interval [t 0, t f], and the Legendre-Gauss (LG) points lie in the interval [−1, +1]. Consider now the following transformation of the independent variable, t, to the variable τ ϵ [−1, +1] (Benson, Reference Benson2005) :
$$t = \displaystyle{{t_f - t_0} \over 2}\tau + \displaystyle{{t_f + t_0} \over 2}$$Furthermore, we choose the collocation points in the phase to be the set of LG points, (τ 1, …, τ N), which are the roots of an Nth-degree Legendre polynomial, P N(τ), given as
$$P_N = \displaystyle{1 \over {2^N N!}}\displaystyle{{d^N} \over {d\tau ^N}} [(\tau ^2 - 1)^N ]$$The state and control variables are approximated by using a basis of N + 1 Lagrange interpolating polynomials, L i(τ)(i = 0, 1, …N)
$$x\left( {\tau (t)} \right) = \sum\limits_{i = 0}^N x\left( {\tau _i} \right)L_i \left( \tau \right)$$
$$u\left( {\tau (t)} \right) = \sum\limits_{i = 1}^N u\left( {\tau _i} \right)L_i \left( \tau \right)$$The derivative of Equation (17) at the LG points τ k(k = 1, 2…, N) can be represented as follows
$$\dot x(\tau _k ) = \displaystyle{2 \over {t_f - t_0}} \sum\limits_l^N D_{kl} x_l $$where D = [D kl] is a (N) × (N + 1) differential approximation matrix,
$$D_{kl} = \dot L_l (\tau _k ) = \sum\limits_{i = 0}^N \displaystyle{{\prod\limits_{\,j = 0,\,j \ne l,i}^N (\tau _k - \tau _j )} \over {\prod\limits_{\,j = 0,\,j \ne l}^N (\tau _l - \tau _j )}}$$where k = 1, …, N and l = 0, 1, …N. So the state dynamic constraint Equation (12) is transcribed into algebraic constraints via the differential approximation matrix as follows
$$\sum\limits_{l = 0}^N D_{kl} x_l - \displaystyle{{t_f - t_0} \over 2}f\,(x_k, u_k ) = 0$$Also terminal constraints and path constraints are transcribed as follows
$$e(x_0, x_N, \tau _0, \tau _f ) = 0$$
$$c(x_k, u_k ) \le 0$$The cost function can be approximated with a Gauss quadrature, resulting in
$$J = E(x_0, x_N, \tau _0, \tau _f ) + \displaystyle{{t_f - t_0} \over 2}\sum\limits_{k = 1}^N F(x_k, u_k )w_k $$Here, w k is the Gauss weights defined as
$$w_k = \displaystyle{2 \over {1 - \tau _l^2}} [\dot L_N (\tau _l )]^2 $$Finally, the optimal control problem is transformed into a NLP problem, it can be summarised as: finding the state parameters XN = [x 0, x 1, …, x N], control parameters UN = [u 1, …, u N] to make the cost function of Equation (24) minimised under the condition that all the constraints of Equations (21) to (23) are satisfied.
Although the Pseudo-spectral Method can be used to calculate an open-loop optimal trajectory in a very short time, for the online vehicle ascent phase optimal guidance problem under model parameters uncertainties and external disturbance Equation (26), a significant barrier to this optimisation-based control strategy is computational delay, especially when the solution time of the NLPP exceeds the sampling time.
$$\dot x = f\,(x,u,t;p) + d(t)$$where p represents the model parameters and d(t) represents the external disturbance.
3.3. NLP Sensitivity Analysis
In recent decades, NLP sensitivity analysis has been applied to numerous applications in process engineering with parameters that represent uncertain data or unknown inputs (Zavala and Biegler, Reference Zavala and Biegler2009; Yang and Biegler, Reference Yang and Biegler2013). In the MPC feedback context, Zavala and Biegler (Reference Zavala and Biegler2009) analysed the impact of sensitivity updates on the stability of the closed loop in the presence of stabilising terminal constraints and Lyapunov-type terminal costs.
After an optimal control problem has been transformed into a NLP problem, we can write the NLP as P(p) which is parametric in the initial state p: = x 0 via Equation (27)
$$\matrix{ {\min J_N (z,p)} \hfill \cr {subject \; to \; C_i (z,p) = 0,i = 1, \ldots ,n_e} \hfill \cr {C_i (z,p) \ge 0,i = n_e + 1, \ldots ,n_c} \hfill \cr} $$where z denotes the non-linear program variable.
The Lagrangian function can be denoted as Equation (28)
$$L(z,\mu, p) = J_N (z,p) + \mu ^T C(z,p)$$The index set of active constraints is shown as Equation (29)
$$\eqalign{& A(z,p) = \{ 1, \ldots ,{n_e}\} \cup \cr & \{ i \vert {C_i}(z,p) = 0,i = {n_e} + 1, \ldots ,{n_c}\} }$$we denote C A(z,p) and μ A(z,p) as the active constraints and the corresponding multipliers respectively.
The NLP Sensitivity Theorem (Fiacco, Reference Fiacco1976) states sufficient conditions for the relationship between an optimal solution z(p) and the parameter p. The theorem states:
Lemma 1: Considing Equation (28), suppose that J N and C are at least twice differentiable in a neighborhood of the nominal solution z *(p 0). If the Linear Independence Constraint Qualification (LICQ), the Sufficient Second order Optimality Conditions (SSOC) and the Strict Complementarity Condition (SCC) are satisfied in this neighborhood, then we have that:
• z *(p 0) is an isolated local minimizer and the respective Lagrange multipliers are unique.
• for where p is in a neighborhood of p 0 there exists a unique local minimizer z *(p) which satisfies LICQ, SSOC and SCC and is differentiable with respect to p.
• there exists a Lipschitz constant L J such that such that
(30)
$$\left\vert {J_N (\,p) - J_N (\,p_0 )} \right\vert \le L_{Jp} \left\Vert {\,p - p_0} \right\Vert$$• there exists a Lipschitz constant L z such that for the solution z
(31)$$\left\vert {z(\,p) - z(\,p_0 )} \right\vert \le L_{zp} \left\Vert {\,p - p_0} \right\Vert$$
Based on the implicit function theorem (Fiacco, Reference Fiacco1976), we can obtain Equation (32)
$$\eqalign{& \cr & \left[ {\matrix{ {\nabla _{zz}^2 L(z^{\ast}, \mu ^{\ast}, p_0 )} & {\nabla _z C_{A(z^{\ast},\, p_0 )} (z^{\ast}, p_0 )^T} \cr {\nabla _z C_{A(z^{\ast},\, p_0 )} (z^{\ast}, p_0 )} & 0 \cr}} \right] \bullet \left[ {\matrix{ {\displaystyle{{\partial z} \over {\partial p}}(\,p_0 )} \cr {\displaystyle{{\partial \mu _{A(z^{\ast}, \;p_0 )}} \over {\partial p}}(\,p_0 )} \cr}} \right] \cr & = - \left[ {\matrix{ {\nabla _{zp}^2 L(z^{\ast}, \mu ^{\ast}, p_0 )} \cr {\nabla _p C_{A(z^{\ast},\, p_0 )} (z^{\ast}, p_0 )} \cr}} \right]} $$
where the matrix 
$\left[ {\displaystyle{{\partial z} \over {\partial p}}(p_0 ),\displaystyle{{\partial \mu _{A(z^{\ast},\, p_0 )}} \over {\partial p}}(p_0 )} \right]^T $ is called the sensitivity matrix. It is simple to obtain a first-order approximation of the optimal solution for a perturbed parameter via Equation (33)
$$z(\,p) = z^{\ast} (\,p_0 ) + \displaystyle{{\partial z} \over {\partial p}}(\,p_0 ) \cdot (\,p - p_0 )$$
where 
$\displaystyle{{\partial z} \over {\partial p}}(p_0 )$ is the sensitivity. During the computation of the optimal control problem, sensitivity information can be obtained very easily with improvements in fast Newton-based barrier methods. The software Interior Point Optimiser (IPOPT) and Sensitivity IPOPT (sIPOPT) (Pirnay et al., Reference Pirnay, López-Negrete and Biegler2012) are used to compute the sensitivity information, where IPOPT solves a non-linear program in the background and sIPOPT provides an estimate of a neighboring solution. Moreover, in Equation (33), 
$\displaystyle{{\partial z} \over {\partial p}}(p_0 )$ is directly available in factored form from the optimal solution by IPOPT, so the sensitivity can be calculated through a simple back solution by sIPOPT.
3.4. Real-time optimal guidance algorithm
If the sampling frequency of the system is fixed, an approximated closed-loop solution of an optimal control problem is computed in the following way: in each sampling interval, we use the nominal non-linear model to predict the next sampling states, an optimal control problem is solved based on the predicted states. When the system dynamics are simulated to the next sampling instant, based on the difference between predicted states and measured states, the resulting optimal control sequence is corrected via sensitivities and is applied as input for the next sampling interval(s). This procedure is then repeated iteratively.
We denote t 0 as the initial time, t i(i = 1, 2, …) represents the sampling instant, ΔT is the fixed sampling period; x(t 0) is the initial state, x(t i) and 
$\hat x({t_i})(i = 1,2, \ldots )$ represent the measured state and the predicted state respectively; 
${u^{\ast}}(\hat x({t_i}),t)$ is the optimal control via GPM and 
$u(\hat x({t_i}),t)$ is the corrected control via sensitivity updates; we denote t pi(i = 1, 2, …) as the predicted time from t i−1 to t i and t gi(i = 1, 2, …) as the computational time of the optimal control problems from t i−1 to t f with the initial state 
$\hat x({t_i})$. The real-time optimal guidance algorithm design is summarised in the following algorithm and is shown in Figure 1.
Figure 1. Real-time optimal guidance algorithm.
Step 1: Choose the discrete parameter N , the fixed sampling period ΔT and the first sampling time t 1; Collect the initial state x 0 = x(t 0) and set i = 1. Use the Pseudo-spectral method to calculate the optimal controller u 0* = u *(x 0,t) and apply it to the system until t 1 = t 0 + ΔT and get the the predicted state 
$\hat x({t_1})$;
Step 2: Collect the real measurement x(t 1), correct the control via sensitivity updates Equation (34)
$$u(x(t_1 ),t) = u_0^{\ast} + \displaystyle{{\partial u_0^{\ast}} \over {\partial \hat x(t_1 )}}(x(t_1 ) - \hat x(t_1 ))$$ Step 3: Propagate the system dynamics with the control u(x(t i−1),t) with x(t i−1) as the initial condition to time t i and get the the predicted state 
$\hat x({t_i})$;
Step 4: Use the Pseudo-spectral method to calculate the optimal controller 
${u^{\ast}}(\hat x({t_i}),t)$ with 
$\hat x({t_i})$ as the initial condition;
Step 5: If t pi + t gi ≤ ΔT, collect the measurement x(t i), correct the control via sensitivity updates to get the control u(x(t i), t) as Equation (35) and go to Step 3; If t pi + t gi > ΔT, during time period [t i, ti +1], apply the control u(x(t i−1), t) and get the predicted state 
$\hat x({t_{i + 1}})$, go to Step 3;
$$u(x(t_i ),t) = u_i^{\ast} + \displaystyle{{\partial u_i^{\ast}} \over {\partial \hat x(t_i )}}(x(t_i ) - \hat x(t_i ))$$We denote that if the control values violate their corresponding constraints, then the control command u(t) is limited to satisfy the active constraint (Jiang et al., Reference Jiang, Teo and Duan2012) as Equation (36).
$$u(t) = \{ v \in U\left\vert {c\left[ {x(t),u(t),t} \right]} \right. = 0\} $$4. FEEDBACK-BASED REQUIREMENTS ON SAMPLING FREQUENCY
Consider two trajectories x(t) and r i(t). x(t) is the system trajectory under the real control u i(t) = u(x(t i),t), r i(t) is the is the optimal trajectory starting from the same initial condition x(t i) and driven by the optimal control u i*(t) = u *(x(t i), t). For each sampling interval [t i, ti +1], where x(t) satisfies
$$\matrix{ {\dot x(t) = f\,(x(t),u_i (t),p) + d(t)} \hfill \cr {x(t_i ) = x(t_i )} \hfill \cr} $$and r i(t) satisfies
$$\matrix{ {\dot r_i (t) = f\,(r_i (t),u_i^{\ast} (t),p_0 )} \hfill \cr {r_i (t_i ) = x(t_i )} \hfill \cr} $$Assumption 1: For the non-linear system Equation (26), model parameters uncertainties and external disturbance are bounded, there exist ε p > 0, ε d > 0, such that
$$\left\Vert {\,p - p_0} \right\Vert \le \varepsilon _p, \left\Vert d \right\Vert_{L_\infty} \le \varepsilon _d $$Assumption 2: The non-linear vector field f (x, u, t;p) is Lipschitz continuous, i.e., there are constants L fx > 0, L fu > 0, L fp > 0, such that:
$$\left\Vert {\,f\,(x,u,p) - f\,(y,v,q)} \right\Vert \le L_{\,fx} \left\Vert {x - y} \right\Vert + L_{\,fu} \left\Vert {u - v} \right\Vert + L_{\,fp} \left\Vert {\,p - q} \right\Vert$$Assumption 3: Given any two initial conditions x 1 and x 2, we denote u *(x 1, t) and u *(x 2, t) as the optimal controls with the corresponding initial conditions. It is assumed that there are constants L ux > 0, such that for all t ϵ [t 0, ∞)
$$\left\Vert {u^{\ast} (x_1, t) - u^{\ast} (x_2, t)} \right\Vert \le L_{ux} \left\Vert {x_1 - x_2} \right\Vert$$This assumption requires that the difference in the optimal controls be linearly bounded by the difference in the initial conditions.
Assumption 4: Errors between the predicted state and the real state caused by model parameters uncertainties and outside disturbance are bounded. There are positive constants δ, such that for all i and t ϵ [t 0, ∞):
$$\left\Vert {x(t_i ) - \hat x(t_i )} \right\Vert \le \delta $$Assumption 5: The relationship between the corrected control via sensitivity updates and the optimal control is as follows, there are positive constants ε u, such that for all i and t ϵ [t i, ti +1]:
$$\left\Vert {u(x(t_i ),t) - u^{\ast} (\hat x(t_i ),t)} \right\Vert = \displaystyle{{\partial u_i^{\ast}} \over {\partial \hat x(t_1 )}}\left\Vert {x(t_i ) - \hat x(t_i )} \right\Vert \le \displaystyle{{\partial u_i^{\ast}} \over {\partial \hat x(t_1 )}}\delta \le \varepsilon _u $$
where 
$\displaystyle{{\partial u_i^{\ast}} \over {\partial \hat x(t_1 )}}$ is the sensitivity. ε u is used to represent the sensitivity and the computational error in the calculation of the optimal controller.
Definition 1: The multivalued function W(x), x ϵ R, given implicitly by
$$W(z)e^{W(z)} = z$$is called the Lambert W function (Ross et al., Reference Ross, Sekhavat, Fleming and Gong2008). For x ≥ 0, W(x) is single-valued.
Lemma 2: Let 
$\left[ {{t_0},{t_f}} \right] \to y(t) \in {R_ + }$ be an integrable function that satisfies Gronwalls inequality (Ross et al., Reference Ross, Sekhavat, Fleming and Gong2008):
$$y(t) \le a(t) + \int_{t_0} ^{t_f} b(s)y(s)ds$$
where a(t), b(t) are continuous, nonnegative, bounded functions, with 
$t: \to a(t)$ nondecreasing over the interval 
$\left[ {{t_0},{t_f}} \right]$; then
$$y(t) \le a(t)e^{B(t)} $$
$$B(t) = \int_{t_0} ^{t_f} b(s)ds$$The process of the proof is similar to the notes (Ross et al., Reference Ross, Sekhavat, Fleming and Gong2008), but additionally, we consider the effect of model parameters, uncertainties, external disturbance and sensitivities. We note that all of the statements below are under the assumption that t pi + t gi ≤ ΔT; if t pi + t gi > ΔT, the ones below should be progressed as indicated in Step 5, Section 3.4.
Theorem 1: Let the controller be designed according to the proposed algorithm. Under Assumption 1–5, for any given ε > 0 if
$$\Delta T \le \Delta T_m $$where ΔT m satisfies
$$\Delta T_m = \displaystyle{{W(z_m )} \over {L_{\,fx}}} $$
$$z_m = \displaystyle{{\varepsilon L_{\,fx}} \over A}$$
$$A = L_{\,fu} (\varepsilon _u + L_{ux} \delta ) + L_{\,fp} \varepsilon _p + \varepsilon _d $$then,
$$\left\Vert {x(t_{i + 1} ) - r_i (t_{i + 1} )} \right\Vert \le \varepsilon $$ Proof: By the Principle of Optimality, 
${u^{\ast}}(\hat x({t_i}),t)$ is the optimal control with the initial condition 
$\hat x({t_i})$. Therefore, from Assumption 3 and Assumption 5, 
$\forall t \in \left[ {{t_i},{t_{i + 1}}} \right]$,
$$\eqalign{\left\Vert {u_i (t) - u_i^{\ast} (t)} \right\Vert & \le \left\Vert {u_i (t) - u^{\ast} (\hat x(t_i ),t)} \right\Vert + \left\Vert {u^{\ast} (\hat x(t_i ),t) - u_i^{\ast} (t)} \right\Vert \cr & \le \varepsilon _u + L_{ux} \left\Vert {\hat x(t_i ) - x(t_i )} \right\Vert \le \varepsilon _u + L_{ux} \delta} $$By Assumption 1 and Assumption 2, we can obtain
$$\eqalign{\left\Vert {x(t) - r_i (t)} \right\Vert & \le L_{\,fx} \int_{t_i} ^t \left\Vert {x(\tau ) - r_i (\tau )} \right\Vert d\tau + L_{\,fu} \Delta T\left\Vert {u_i (t) - u_i^{\ast} (t)} \right\Vert \cr &\quad + L_{\,fp} \Delta T\left\Vert {\,p - p_0} \right\Vert + \Delta T\left\Vert d \right\Vert_\infty \cr &\quad \le L_{\,fx} \int_{t_i} ^t \left\Vert {x(\tau ) - r_i (\tau )} \right\Vert d\tau + A\Delta T}$$where A = L fu(ε u + L uxδ) + L fpεp + ε d. Moreover, by Gronwalls inequality (Ross et al., Reference Ross, Sekhavat, Fleming and Gong2008), from Equation (46) we can get:
$$\left\Vert {x(t) - r_i (t)} \right\Vert \le A\Delta Te^{L_{\,fx} \Delta T} $$When t = t i+1,
$$\left\Vert {x(t_{i + 1} ) - r_i (t_{i + 1} )} \right\Vert \le (A/L_{\,fx} )L_{\,fx} \Delta Te^{L_{\,fx} \Delta T} $$Based on Definition 1 of the Lambert W function, it can be easily verified that for z m, c ϵ R +, if c ≤ W(z m), then cec ≤ z m. Hence letting c = ΔTL fx, from Equations (42) to (44), we have
$$\left\Vert {x(t_{i + 1} ) - r(t_{i + 1} )} \right\Vert \le (A/L_{\,fx} )z_m = \varepsilon $$According to the property of the W function, there exists a maximal ΔT m satisfying
$$\Delta T_m = \displaystyle{{W(\varepsilon L_{\,fx} /A)} \over {L_{\,fx}}} $$Remark 1: From Equations (47) and (48) we can see that errors between the real state and the optimal state of the nominal model are bounded and are affected by the sampling period, the sensitivity, the computational error, model parameters uncertainties and external disturbance.
Remark 2: The sample period is mainly determined by the Lipschitz constants of the system. The W function is an increasing function. The maximal ΔT m can be obtained by Equation (42).
Remark 3: For the infinite-horizon feedback control problem, optimal control has the advantage of guaranteeing stability through the Bellman value function serving as a natural Lyapunov function (Clarke, Reference Clarke2004). However, the online vehicle guidance problem is a finite-horizon feedback control problem. Thus, if open-loop optimal controls are generated starting at t i ϵ π for each x(t i), it is possible to guide a vehicle to its target by generating a Carathodory-π feedback solution (Ross et al., Reference Ross, Sekhavat, Fleming and Gong2008). The proposed method can generate a corrected Carathodory-π feedback solution online.
5. SIMULATION RESULTS AND ANALYSIS
5.1. Simulation parameters
In this section, we present the trajectory optimisation problem proposed in the previous sections. Gauss Pseudospectral Optimisation Software (GPOPS) (Rao et al., Reference Rao, Benson, Darby, Patterson, Francolin, Sanders and Huntington2010; Patterson and Rao, Reference Patterson and Rao2014) and sIPOPT (Pirnay et al., Reference Pirnay, López-Negrete and Biegler2012) are used as the simulation software. During the ascent phase of the trajectory, the parameters used in solving the problem are as follows (Prasanna et al., Reference Prasanna, Ghose, Bhat, Bhattacharyya and Umakant2005): reference area S = 1 m 2, the empty mass of the aircraft m 1 = 2000 kg, the quality of fuel m 2 = 1600 kg, the initial mass of the aircraft m 0 = m 1 + m 2 = 3600 kg, the gravitational constant μ = 3·9853 × 1014m 3/s 2, g = 9·8 m/s 2 and R e = 6378000 m is the radius of the earth. The number of LG point is 20. The other parameters used in simulation are shown in Tables 1 to 6.
Table 1. Boundary Conditions.
Table 2. Constraint Conditions.
Table 3. Lift Coefficients.
Table 4. Drag Coefficients.
Table 5. Specific impulse (I sp in seconds).
Table 6. Mass flow rate of air (-kg/s) corresponds to altitude 32·5 km.
At 32·5 km, let the density be ρ 32·5, speed of sound be C 32·5, and the speed corresponding to a given mach number M is v 32·5 = MC 32·5. 
${\dot m_{air(32\cdot5)}}$ is shown in Table 6.
The hypersonic propulsion system has been modelled using data on air mass flow rate (
${\dot m_{air}}$) and specific impulse (I sp) as functions of Mach number (M), angle of attack (α) and altitude (h), respectively. From (
${\dot m_{air}}$) the thrust (T) can be found using Equations (50) and (51), where 
${\dot m_f} = {\dot m_{air}}\phi /15$, φ is the equivalence ratio and during ascent phase its value is 1. A is the nozzle area.
$$\underbrace {\dot m_{air}} _{at \; given \; altitude} = \displaystyle{{\dot m_{air(32\cdot5)}} \over {\rho _{32\cdot5} v_{32\cdot5} A}}\underbrace {\rho vA}_{at \; given \; altitude}$$
$$T = I_{sp} \dot m_f g$$5.2. Analysis of the optimised results
Solving optimal control problems by means of the Pseudo-spectral method implies that the differential equation is approximated by a discrete time counterpart. In order to assess the accuracy of the discretization, we may simulate the system using the optimal control profile as input. With this approach, the state profiles are computed with high accuracy and the result may then be compared with the profiles resulting from optimisation. Notice that this procedure does not verify the optimality of the resulting optimal control profiles, but only the accuracy of the discretization of the dynamics. The comparison of the two states is shown in Figures 2 and 3 and the number of discretization points used in the optimal control solution is 30 in Figures 2 and 3. The magnitude of relative errors between them are 10−3 ~ 10−4, so the optimised results are credible.
Figure 2. Feasibility analysis of the altitude curve.
Figure 3. Feasibility analysis of the speed curve.
5.3. Simulation of online trajectory reconstruction
The external disturbance will lead the trajectory to deviate from the standard trajectory during the flight. Based on the current measured states, the trajectory should be re-planned to meet the terminal and path constraints from the proposed method. The sampling period for starting the online trajectory reconstruction algorithm and updating the guidance commands is defined as 5 s. To assess the capability of the proposed method to deal with disturbance, the aerodynamic coefficients variations in the simulations are given by following cases:
• Closed-loop simulation based on Pseudo-spectral method optimal control with sensitivity updates: set drag coefficient error is +10% and lift coefficient error is −10%
• Closed-loop simulation base on Pseudo-spectral method optimal control without sensitivity updates: set drag coefficient error is +10% and lift coefficient error is −10%
• Open-loop simulation with off-line optimal control command: set drag coefficient error is +10% and lift coefficient error is −10%
• Simulation with off-line optimal control command in the normal case without error
The trajectory without reconstruction is defined as the open-loop trajectory. As a contrast, Table 7 lists the deviation of the re-planning trajectory with sensitivity updates and the re-planning trajectory without sensitivity updates for the terminal constraints and the max trajectory re-planning time during flight. The results are shown from Figure 4 to Figure 8. Figure 4 illustrates the altitude-versus-time profiles in the four cases. Figure 5 illustrates the velocity-versus-time profiles in the four cases. As expected, the flight trajectory with the online trajectory reconstruction algorithm can arrive at the terminal altitude. If there is no trajectory reconstruction, the open-loop flight trajectory cannot meet the terminal altitude constraint. Figure 6 illustrates the control-versus-time profiles in the four cases. The oscillatory behavior of the angle of attack near the end of the trajectory is affected by the sampling period of the online trajectory reconstruction algorithm. In our paper, the sampling period for starting the online trajectory reconstruction algorithm and updating the guidance commands is defined as 5 s. This is because the simulation environment is Matlab; the sampling period should be much shorter in a C or C++ language environment and the oscillatory behavior of the angle of attack near the end of the trajectory should be improved. Figure 7 illustrates the dynamic presssure-versus-time profiles in the four cases. If there is no trajectory reconstruction, the open-loop flight trajectory will exceed the path constraints. Figure 8 illustrates the mass-versus-time profiles in the four cases. The online trajectory reconstruction algorithm with sensitivity updates can consume less mass. The time to reconstruct the trajectory is shorter than 5 s, it should be much faster in a C or C++ language environment. The re-planning trajectory with sensitivity updates has indeed improved the performance indexes from the data in Table 7.
Figure 4. The altitude curve.
Figure 5. The velocity curve.
Figure 6. The control curve.
Figure 7. The dynamic pressure curve.
Figure 8. The mass curve.
Table 7. Reconstruction precision and time on disturbance on-line.
6. CONCLUSIONS AND FUTURE WORK
For hypersonic ascent phase optimal guidance problems, an online optimal guidance algorithm based on Pseudo-spectral method and sensitivity analysis is proposed. The proposed closed-loop feedback method can successively generate online open-loop suboptimal controls without the design procedure of an inner-loop feedback controller. Considering model parameter uncertainties and external disturbance, a sampling theorem is proposed which indicates the effect of the Lipschitz constant on the dynamics of sampling frequency. The simulation results indicate that the method proposed above offers improved performance and has ability for online calculation. Future work will propose a guidance law for tracking the re-planning ascent phase trajectory and estimating aerodynamic parameters and density online to get higher predicted accuracy.
ACKNOWLEDGMENTS
This work was supported in part by the National Nature Science Foundation of China nos. 61473124, 61203081 and 61174079, Doctoral Fund of Ministry of Education of China no. 20120142120091, and Precision manufacturing technology and equipment for metal parts no. 2012DFG70640.
