# Optimal Ascent Guidance for Air-Breathing Launch Vehicle Based on Optimal Trajectory Correction

Authors: Xuefang Lu; Yongji Wang; Lei Liu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

An optimal guidance algorithm for air-breathing launch vehicle is proposed based on optimal trajectory correction. The optimal trajectory correction problem is a nonlinear optimal feedback control problem with state inequality constraints which results in a nonlinear and nondifferentiable two-point boundary value problem (TPBVP). It is difficult to solve TPBVP on-board. To reduce the on-board calculation cost, the proposed guidance algorithm corrects the reference trajectory in every guidance cycle to satisfy the optimality condition of the optimal feedback control problem. By linearizing the optimality condition, the linear TPBVP is obtained for the optimal trajectory correction. The solution of the linear TPBVP is obtained by solving linear equations through the Simpson rule. Considering the solution of the linear TPBVP as the searching direction for the correction values, the updating step size is generated by linear search. Smooth approximation is applied to the inequality constraints for the nondifferentiable Hamiltonian. The sufficient condition for the global convergence of the algorithm is given in this paper. Finally, simulation results show the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

The development of space technology has given rise to the expectation that launchers will become low cost and fully reusable. The launch vehicle with hypersonic air-breathing propulsion is considered to reduce the cost of payloads taken to the Earth’s orbit. The air-breathing launch has inherent features that make it a candidate for future space transportation [1]. The impulse of air-breathing propulsion, which is approximately 3000 s [2], is significantly higher than that of a rocket (360 s). Air-breathing propulsion brings high impulse as well as strong nonlinear, coupling aerodynamic and thrust.The traditional design category for ascent guidance is to drive an optimal nominal trajectory off-board. The guidance problem is then transformed into a tracking problem for the designed optimal nominal trajectory. The design of a trajectory could be formulated as a global optimization task [3]. The methods for numerical optimization of continuous dynamic systems could be termed “Hamiltonian” (indirect method) and “Transcription” (direct method) [4]. For the optimization problem with inequality constraints, smoothing approximation was considered in [5, 6]. A filled function approach for nonsmooth constrained global optimization was presented in [7]. Linear-quadratic optimization was implemented to optimal control in [8]. Direct and indirect optimization methods were implemented to the off-board trajectory optimization problem in previous literature. The numerical algorithms of trajectory optimization for vehicles were summarized and systematically analyzed in [9, 10]. A new concept of pseudocontrol sets to solve optimal control problems was proposed in [11]. This approach reduces the calculation cost by combining large-scale linear programming algorithms with discretization of the continuous system dynamics on small segments. An algorithm for multiobjective optimization was presented in [12]. Intelligent algorithms can also be used for trajectory optimization problems. The particle swarm optimization (PSO) method was implemented to the space trajectory optimization in [13]. The simulation results showed the effective of PSO in finding the optimal solution to the space trajectory optimization problems, with great numerical accuracy. Approximate numerical methods of optimization were presented for multiorbit noncoplanar orbit transfers of low-thrust spacecraft in [14].Many control methods were implemented to trajectory tracking guidance and control problem [15]. In the previous literature, many researchers focused on the robust control method [16, 17]. It has obtained successful application in the industry [18, 19]. For the networked control system, a H-infinity step tracking control method was presented in [20]. An adaptive fuzzy robust control for a class of nonlinear systems was proposed in [21]. An adaptive guidance law and off-board trajectory optimization for air-breathing launch vehicle were addressed in [22]. In that paper, the optimal control problem was solved using SQP method. And the adaptive guidance law was developed using a feedback loop based on a second-order rate controller for angle of attack. A robust state feedback guidance law was generated in real time using the indirect Legendre pseudospectral feedback method in [23]. In that paper, the guidance problem was converted into a trajectory state regulation problem which is a linear time varying system.However, the accuracy of the trajectory tracking method is low with the disturbance and the modeling error. It lacks the autonomy and adaptability to cope with the nonnominal vehicle and mission conditions needed for future reusable launch vehicles [24]. To improve the performance and the accuracy of the guidance, the optimal guidance method generates the prospective trajectory by fast trajectory optimization on-board based on current flight state. This guidance method has become a research focus with the advance in on-board computation capability. Ping Lu used the indirect method to pose the trajectory optimization problem as a nonlinear two-point boundary value problem (TPBVP) in [25, 26]. From this model, the optimal thrust vector that satisfied the optimality condition was derived. A finite difference method was employed to solve the nonlinear TPBVP through numerical calculations. Similarly, a fast trajectory optimization for hypersonic air-breathing vehicles was presented in [27]. The indirect method was implemented for ascent trajectory optimization on-board considering the features of the air-breathing vehicle. Ping Lu’s work is significant for the optimal ascent guidance. However, the TPBVP is nonlinear and nondifferentiable for the air breathing launch vehicle. It is difficult to solve the nonlinear and nondifferentiable TPBVP on-board in every guidance cycle. Using the direct method, a new guidance concept based on nonlinear programming (NLP) method was proposed in [24]. NLP-based guidance concepts appear advantageous over conventional methods because the on-board guidance algorithm allows a single algorithm to be implemented for different vehicles and missions. The optimal control problem was parameterized into a nonlinear programming problem that was solved by the gradient projection algorithm. In [28], the reference trajectory was updated for disturbance by an on-board algorithm that satisfied the real-time requirement. A new real-time guidance method derived from the optimality condition was proposed in [29]. In that paper, the simple guidance parameters were updated in real time. In [30], a guidance method was presented for online launcher ascent trajectory updating based on neural networks. In the paper, the utilization of a neural network approximation was used online during the ascent flight, with a training process performed off-line.In this paper, we present an optimal guidance algorithm for air breathing launch vehicle based on optimal trajectory correction to reduce the on-board calculation cost. Considering the current vehicle state as the initial condition, the optimal trajectory correction problem is referred to as a nonlinear optimal control problem with state inequality constraints. For the real-time requirement of the on-board algorithm, the linear TPBVP is obtained for optimal trajectory correction by linearizing the optimality condition in this paper. The Simpson rule is applied to transform the linear TPBVP into linear equations. Considering the solution of the linear TPBVP as the searching direction for the correction values, the updating step size is generated by linear search. Smooth approximation is applied to the inequality constraints for the nondifferentiable Hamiltonian. The sufficient condition for the global convergence of the algorithm is given in this paper. Finally, simulation results for different cases of the modeling error show the effectiveness of the proposed algorithm. In summation, the main contributions of this paper are given as follows.(1)
We reduce the on-board calculation cost of the guidance method a lot. Comparing with the methods in previous literatures which solves nonlinear equations in every guidance cycle, the proposed guidance method in this paper solves linear equations once in every guidance cycle only.(2)
We obtain the relationship between the global convergence and the guidance cycle of the online algorithm. The sufficient condition of the global convergence of this algorithm is given.The remainder of this paper is organized as follows. Section2 presents the state normalized-energy differential equations of motion of the vehicle and the optimal control problem for optimal ascent guidance. Section 3 provides details on the optimal guidance algorithm. Section 4 presents an analysis of the global convergence of the proposed algorithm. Section 5 discusses the proposed differentiable approximation for Hamiltonian. Section 6 presents the simulations for the Generic Hypersonic Vehicle (GHV) model and scramjet engine. Section 7 presents the conclusion.


## 2. Problem Formulation

### 2.1. Ascent Dynamic

The formulation of motion of the air-breathing vehicle in the longitudinal plane is presented in this section. As shown in Figure1, the thrust T, gravitation g, and the aerodynamic lift force L and drag force D are acting on the vehicle. The angle of attack α is the angle between the body axis and the vector of the velocity v. The flight path angle γ is the angle between the vector of the velocity v and the ground plane.Figure 1
The forces acting on the vehicle.The motion of the vehicle in the longitudinal plane can be described as follows:(1)drdt=vsinγ,dvdt=Tcosα-Dm-gsinγ,dγdt=1v[Tsinα+Lmcosσ-gcosγ+v2rcosγ],dmdt=-TIsp,
where r, m, Isp, and g are the altitude, weight of the vehicle, impulse, and gravity acceleration, respectively.The lift forceL and the drag force D are given as follows:
(2)L=12ρv2SrefCL(Ma,α),D=12ρv2SrefCD(Ma,α),
where ρ and Sref are the air density of the current altitude and the reference area respectively, CL and CD are the lift coefficient and the drag coefficient, respectively, which are the nonlinear functions of the angle of attack α and Mach number Ma. For the air-breathing engine, the thrust is given by
(3)T=T(ρ,Ma,α,ϕ),
where T is the nonlinear function of ρ, α, Ma, and the throttle command ϕ.In [25], the equations were normalized to reduce the numerical calculation error. Variable substitutions are applied as follows:
(4)R=rr0,V=vg0r0,g=g0r02r2,
where r0 and g0 are Earth radius and the ground gravity acceleration, respectively. Considering the thrust T and aerodynamic forces as a composition of forces, we define the normalized accelerations aT and aL as follows:
(5)aT=Tcosα-Dmg0,aL=Tsinα+Lmg0.
We define
(6)τ=tr0/g0.
The normalized equations of motion are obtained from (1), (4), (5), and (6) as follows:
(7)dRdτ=Vsinγ,dVdτ=aT-sinγR2,dγdτ=1V[aL+(V2-1R)(cosγR)],dmdτ=-Tr0Ispg0.The mission of ascent is to send the vehicle to the required final state which is denoted as(Rf,Vf,γf). Define normalized-energy E as follows:
(8)E=V22-1R.
From (8) we obtain
(9)dEdτ=VaT,dRdE=sinγaT,dVdE=1aTV[aT-sinγR2],dγdE=1aTV2[aL+(V2-1R)(cosγR)],dmdE=-Tr0aTVIspg0.
The integration interval [E0,Ef] is fixed.


### 2.2. Optimal Feedback Control Problem for Ascent Guidance

In this section, the optimal feedback control problem for optimal ascent guidance is addressed. This problem differs from trajectory optimization off-board in that the initial state of the optimal control problem is obtained from the navigation system. We denote the state and the guidance command of the trajectory as(10)x(E)=[R(E)V(E)γ(E)m(E)]T,u(E)=[α(E)ϕ(E)]T.
With the feedback vehicle state xck obtained from navigation system, the optimal ascent guidance algorithm is used to generate new reference trajectory [x(E)u(E)]T and output u(E) from the following optimal control problem:
(11)minJ=ϕ(x(Ef))+∫Eck+1Ef‍g(x,u)dE,.s.t.00x˙=f(x,u),.s.t.00x(Eck)=xck,.s.t.00W(x(Ef))=0,.s.t.00C(x,u)≤0,
where Eck is the normalized energy for current vehicle state xck.


## 2.1. Ascent Dynamic

The formulation of motion of the air-breathing vehicle in the longitudinal plane is presented in this section. As shown in Figure1, the thrust T, gravitation g, and the aerodynamic lift force L and drag force D are acting on the vehicle. The angle of attack α is the angle between the body axis and the vector of the velocity v. The flight path angle γ is the angle between the vector of the velocity v and the ground plane.Figure 1
The forces acting on the vehicle.The motion of the vehicle in the longitudinal plane can be described as follows:(1)drdt=vsinγ,dvdt=Tcosα-Dm-gsinγ,dγdt=1v[Tsinα+Lmcosσ-gcosγ+v2rcosγ],dmdt=-TIsp,
where r, m, Isp, and g are the altitude, weight of the vehicle, impulse, and gravity acceleration, respectively.The lift forceL and the drag force D are given as follows:
(2)L=12ρv2SrefCL(Ma,α),D=12ρv2SrefCD(Ma,α),
where ρ and Sref are the air density of the current altitude and the reference area respectively, CL and CD are the lift coefficient and the drag coefficient, respectively, which are the nonlinear functions of the angle of attack α and Mach number Ma. For the air-breathing engine, the thrust is given by
(3)T=T(ρ,Ma,α,ϕ),
where T is the nonlinear function of ρ, α, Ma, and the throttle command ϕ.In [25], the equations were normalized to reduce the numerical calculation error. Variable substitutions are applied as follows:
(4)R=rr0,V=vg0r0,g=g0r02r2,
where r0 and g0 are Earth radius and the ground gravity acceleration, respectively. Considering the thrust T and aerodynamic forces as a composition of forces, we define the normalized accelerations aT and aL as follows:
(5)aT=Tcosα-Dmg0,aL=Tsinα+Lmg0.
We define
(6)τ=tr0/g0.
The normalized equations of motion are obtained from (1), (4), (5), and (6) as follows:
(7)dRdτ=Vsinγ,dVdτ=aT-sinγR2,dγdτ=1V[aL+(V2-1R)(cosγR)],dmdτ=-Tr0Ispg0.The mission of ascent is to send the vehicle to the required final state which is denoted as(Rf,Vf,γf). Define normalized-energy E as follows:
(8)E=V22-1R.
From (8) we obtain
(9)dEdτ=VaT,dRdE=sinγaT,dVdE=1aTV[aT-sinγR2],dγdE=1aTV2[aL+(V2-1R)(cosγR)],dmdE=-Tr0aTVIspg0.
The integration interval [E0,Ef] is fixed.


## 2.2. Optimal Feedback Control Problem for Ascent Guidance

In this section, the optimal feedback control problem for optimal ascent guidance is addressed. This problem differs from trajectory optimization off-board in that the initial state of the optimal control problem is obtained from the navigation system. We denote the state and the guidance command of the trajectory as(10)x(E)=[R(E)V(E)γ(E)m(E)]T,u(E)=[α(E)ϕ(E)]T.
With the feedback vehicle state xck obtained from navigation system, the optimal ascent guidance algorithm is used to generate new reference trajectory [x(E)u(E)]T and output u(E) from the following optimal control problem:
(11)minJ=ϕ(x(Ef))+∫Eck+1Ef‍g(x,u)dE,.s.t.00x˙=f(x,u),.s.t.00x(Eck)=xck,.s.t.00W(x(Ef))=0,.s.t.00C(x,u)≤0,
where Eck is the normalized energy for current vehicle state xck.


## 3. Optimal Ascent Guidance Algorithm

In [25, 26], the optimal control problem was solved on-board in every guidance cycle. However, it is difficult to solve the nonlinear and nondifferentiable TPBVP on-board for air-breathing launch vehicles. In this section, we propose an optimal guidance algorithm. This guidance algorithm updates the reference trajectory [xu]T in every guidance cycle to deal with the unknown disturbance. For the current state xck of the kth guidance cycle, the optimal feedback control problem is transformed into a linear TPBVP problem using optimality condition linear approximation. The searching direction of the correction values dk is obtained by solving the linear TPBVP. With the searching direction, the new reference trajectory is generated by linear search. For the kth guidance cycle, the trajectory update is performed as follows:
(12)[xk+1(E)uk+1(E)]T=[xk(E)uk(E)]T+lkdkT(E).The initial reference trajectory[x0(E)u0(E)]T is derived from the off-board from trajectory optimization by direct method. The flow chart of the algorithm is shown in Figure 2. Figure 3 shows the time sequence of algorithm.Figure 2
Flow chart of the optimal ascent guidance algorithm.Figure 3
Time sequence of the optimal ascent guidance algorithm.


### 3.1. Linear Approximation for the Optimality Condition

In this section, we linearize the optimality condition for optimality trajectory updating. The optimal control problem (11) without inequality constraints is given as follows:
(13)minJ=ϕ(x(Ef))+∫Eck+1Ef‍g(x,u)dE,.s.t.00x˙=f(x,u),.s.t.00x(Eck)=xck,.s.t.00W(x(Ef))=0.
The Hamiltonian of (13) is given as follows:
(14)H=g(x,u)+λTf(x,u),
where λ is the costate vector. The state equations, costate equations, and optimality condition for the optimal control problem are given by
(15)x˙=f(x,u),λ˙=-(∂H∂x)T,∂H∂u=0.
The initial and transversality conditions are given by
(16)x(Eck)=xck,W(x(Ef))=0,λ(Ef)=[∂ϕ∂x(Ef)]T+[ϖT∂W∂x(Ef)]T.
We denote the trajectory correction variables as follows:
(17)[Δx(E)Δu(E)Δλ(E)]=[x(E)-xk(E)u(E)-uk(E)λ(E)-λk(E)].
We denote that
(18)Huxk=∂2(g+λTf)∂u∂x(xk,uk,λk);Huk=[∂H∂u(xk,uk,λk)]T;Huuk=∂2(g+λTf)∂u2(xk,uk,λk);Hxk=[∂H∂x(xk,uk,λk)]T;Hxxk=∂2(g+λTf)∂x2(xk,uk,λk);fuk=∂f∂u(xk,uk);Hxuk=∂2(g+λTf)∂x∂u(xk,uk,λk);fxk=∂f∂x(xk,uk);guk=[∂g∂u(xk,uk)]T;gxk=[∂g∂x(xk,uk)]T.
The first-order Taylor expansion of (15) on the old reference trajectory [xkukλk]T is
(19)[∂H∂u]T≈Huk+HuxkΔx+HuukΔu+(fuk)TΔλ,λ˙≈-Hxk-HxxkΔx-Hxuk(xk,uk,λk)Δu-(fxk)TΔλ,x˙≈f(xk,uk)+fxkΔx+fukΔu.
The following equations are established:
(20)Huk+(fuk)TΔλ=guk+(fuk)Tλ,Hxk+(fxk)TΔλ=gxk+(fxk)Tλ.
Substituting (20) into (19) yields
(21)λ˙=-gxk-HxxkΔx-Hxuk(xk,uk,λk)Δu-(fxk)Tλ,Δx˙=fxkΔx+fukΔu-x˙k+f(xk,uk),guk+HuxkΔx+HuukΔu+(fuk)Tλ=0,Δx(Eck)+Δxk(Eck)=xck,∂W∂x(Ef)Δx(Ef)+W(xk(Ef))=0,λ(Ef)=[∂ϕ∂x(Ef)]T+[ϖT∂W∂x(Ef)]T.
This equation is a linear TPBVP about the variables [ΔxΔuλ]T. The linear TPBVP is solved using the Simpson rules in Section 6.


### 3.2. Backtracking Line Search for the Step Size

In this section, the step sizelk is determined by backtracking line search considering the solution of the linear TPBVP as the search direction. The penalty function for problem (13) is given by
(22)pk=J(x,u)+η{∥W(xf)∥+∫EckEf‍∥f(x,u)-x˙∥0000000000∫EckEf‍+∥x(Eck)-xck∥},
where η is the penalty parameter satisfied η>max|λ|. We denote
(23)y=[xkuk],dk=[ΔxΔu],
where [ΔxΔuλ]T is the solution of the linear TPBVP (21).The minimizer of penalty functionpk(x,u) corresponds to the solution of problem (11). Backtracking line search was used for large-scale nonlinear optimization problems in [31]. The step size lk is determined by backtracking line search that satisfies
(24)lk=max{l∈S:pk+1(yk+ldk)00000≤pk+1(yk)+lε1δpk+1(yk,dk)},S={l:l=ε2m,m=0,1,2,…},
where ε1∈(0,1), ε2∈(0,1), and δpk+1(yk,dk) is the variational of the functional pk+1(y) caused by variational dk of y on yk. Considering dk as the solution of (21), δpk+1(yk,dk)≤0 is established. The new state and guidance command of the reference trajectory are generated as follows:
(25)[xk+1uk+1]=[xkuk]+lk[ΔxΔu].


## 3.1. Linear Approximation for the Optimality Condition

In this section, we linearize the optimality condition for optimality trajectory updating. The optimal control problem (11) without inequality constraints is given as follows:
(13)minJ=ϕ(x(Ef))+∫Eck+1Ef‍g(x,u)dE,.s.t.00x˙=f(x,u),.s.t.00x(Eck)=xck,.s.t.00W(x(Ef))=0.
The Hamiltonian of (13) is given as follows:
(14)H=g(x,u)+λTf(x,u),
where λ is the costate vector. The state equations, costate equations, and optimality condition for the optimal control problem are given by
(15)x˙=f(x,u),λ˙=-(∂H∂x)T,∂H∂u=0.
The initial and transversality conditions are given by
(16)x(Eck)=xck,W(x(Ef))=0,λ(Ef)=[∂ϕ∂x(Ef)]T+[ϖT∂W∂x(Ef)]T.
We denote the trajectory correction variables as follows:
(17)[Δx(E)Δu(E)Δλ(E)]=[x(E)-xk(E)u(E)-uk(E)λ(E)-λk(E)].
We denote that
(18)Huxk=∂2(g+λTf)∂u∂x(xk,uk,λk);Huk=[∂H∂u(xk,uk,λk)]T;Huuk=∂2(g+λTf)∂u2(xk,uk,λk);Hxk=[∂H∂x(xk,uk,λk)]T;Hxxk=∂2(g+λTf)∂x2(xk,uk,λk);fuk=∂f∂u(xk,uk);Hxuk=∂2(g+λTf)∂x∂u(xk,uk,λk);fxk=∂f∂x(xk,uk);guk=[∂g∂u(xk,uk)]T;gxk=[∂g∂x(xk,uk)]T.
The first-order Taylor expansion of (15) on the old reference trajectory [xkukλk]T is
(19)[∂H∂u]T≈Huk+HuxkΔx+HuukΔu+(fuk)TΔλ,λ˙≈-Hxk-HxxkΔx-Hxuk(xk,uk,λk)Δu-(fxk)TΔλ,x˙≈f(xk,uk)+fxkΔx+fukΔu.
The following equations are established:
(20)Huk+(fuk)TΔλ=guk+(fuk)Tλ,Hxk+(fxk)TΔλ=gxk+(fxk)Tλ.
Substituting (20) into (19) yields
(21)λ˙=-gxk-HxxkΔx-Hxuk(xk,uk,λk)Δu-(fxk)Tλ,Δx˙=fxkΔx+fukΔu-x˙k+f(xk,uk),guk+HuxkΔx+HuukΔu+(fuk)Tλ=0,Δx(Eck)+Δxk(Eck)=xck,∂W∂x(Ef)Δx(Ef)+W(xk(Ef))=0,λ(Ef)=[∂ϕ∂x(Ef)]T+[ϖT∂W∂x(Ef)]T.
This equation is a linear TPBVP about the variables [ΔxΔuλ]T. The linear TPBVP is solved using the Simpson rules in Section 6.


## 3.2. Backtracking Line Search for the Step Size

In this section, the step sizelk is determined by backtracking line search considering the solution of the linear TPBVP as the search direction. The penalty function for problem (13) is given by
(22)pk=J(x,u)+η{∥W(xf)∥+∫EckEf‍∥f(x,u)-x˙∥0000000000∫EckEf‍+∥x(Eck)-xck∥},
where η is the penalty parameter satisfied η>max|λ|. We denote
(23)y=[xkuk],dk=[ΔxΔu],
where [ΔxΔuλ]T is the solution of the linear TPBVP (21).The minimizer of penalty functionpk(x,u) corresponds to the solution of problem (11). Backtracking line search was used for large-scale nonlinear optimization problems in [31]. The step size lk is determined by backtracking line search that satisfies
(24)lk=max{l∈S:pk+1(yk+ldk)00000≤pk+1(yk)+lε1δpk+1(yk,dk)},S={l:l=ε2m,m=0,1,2,…},
where ε1∈(0,1), ε2∈(0,1), and δpk+1(yk,dk) is the variational of the functional pk+1(y) caused by variational dk of y on yk. Considering dk as the solution of (21), δpk+1(yk,dk)≤0 is established. The new state and guidance command of the reference trajectory are generated as follows:
(25)[xk+1uk+1]=[xkuk]+lk[ΔxΔu].


## 4. The Global Convergence of the Algorithm

In [32, 33], the convergence of the off-board algorithm was analyzed. But the relationship between the global convergence and the guidance cycle τ has not been discussed for on-board optimization in the previous literatures. In this section, we state and prove the global convergence of the on-board algorithm proposed in the preceding section. The solution sequence [xkuk]T generated by the algorithm will converge to the solution of the feedback optimal control problem. Unlike off-board trajectory optimization, the initial condition of the feedback optimal control problem depends on the current state of the vehicle in real time. For the kth guidance cycle, the optimal feedback control problem is described as (13). The penalty function pk for this feedback optimal control problem is given by (22). For the subsequent guidance cycle, the optimal feedback guidance problem is given by
(26)minJ=W(x(Ef))+∫Eck+1Ef‍g(x,u)dE,.s.t.00x˙=f(x,u),.s.t.00x(Eck+1)=xck+1,.s.t.00W(x(Ef))=0,
where
(27)Eck+1=Eck+∫τkτk+Δτ‍dEdτdτ,Δτ is the guidance cycle  and xck+1 is the new state of vehicle from navigation system at the time τk+Δτ. The penalty function for optimal problem (24) is given by
(28)pk+1(x,u)=J(x,u)+η∥W(xk(Ef))∥+η∫Eck+1Ef‍∥f(x,u)-x˙∥dE+η∥x(Eck+1)-xck+1∥.
We define the entire penalty function as follows:
(29)p~k(x,u)=pk(x,u)+∑i=0k-1‍∫EciEci+1‍[η∥f(xi,ui)-x˙i∥+g(xi,ui)]dE.
The entire penalty function p~k(x,u) includes the penalty function on the whole integration interval [E0,Ef].Theorem 1.
If the minimizer ofp~k(xk,uk) is exited, g(x,u)≥0, and Δτ is selected and satisfies
(30)η∥xck+1-xck-∫EckEck+1‍x˙kdE∥-δpk+1(yk,dk)<lkε3,
where ε3∈(0,ε1). The solution sequence [xkuk]T generated from (24) will converge to the solution of the feedback optimal control problem.Proof.
From (22), (28), (29) and
(31)xk(Eck+1)=xk(Eck)+∫EckEck+1‍x˙kdE
we can obtain that
(32)p~k+1(xk,uk)-p~k(xk,uk)=η∥(xck+1-xck-∫EckEck+1x˙kdE‍)xk(Eck)-xck0000000-(xck+1-xck-∫EckEck+1x˙kdE)∥0000-η∥xk(Eck)-xck∥.
From (32), we obtain
(33)limΔτ→0[p~k+1(xk,uk)-p~k(xk,uk)]=0,limf(xk,uk)-x˙k→0[p~k+1(xk,uk)-p~k(xk,uk)]=0,(34)p~k+1(xk,uk)-p~k(xk,uk)=η{∥xk(Eck)-xck+1+∫EckEck+1‍x˙kdE∥0000000∥xk(Eck)-xck+1+∫EckEck+1‍x˙kdE∥-∥xk(Eck)-xck∥}≤η{∥xk(Eck)-xck∥∥xck+1-xck-∫EckEck+1‍x˙kdE∥0000000+∥xck+1-xck-∫EckEck+1‍x˙kdE∥0000000∥xck+1-xck-∫EckEck+1‍x˙kdE∥-∥xk(Eck)-xck∥}=η∥xck+1-xck-∫EckEck+1‍x˙kdE∥.
From conditions (24) and (30), it is obtained that
(35)p~k+1(xk+1,uk+1)≤p~k+1(xk,uk)+lkε1δpk+1(dk)=p~k(xk,uk)+lkε1δpk+1(dk)+η{∥xk(Eck)-xck(xck+1-xck-∫EckEck+1x˙kdE‍)000000000-(xck+1-xck-∫EckEck+1x˙kdE‍)∥00000000∥(xck+1-xck-∫EckEck+1x˙kdE‍)xk(Eck)-xck-∥xk(Eck)-xck∥}≤p~k(xk,uk)+η∥xck+1-xck-∫EckEck+1x˙kdE‍∥+lkε1δpk+1(yk,dk)≤p~k(xk,uk)+η∥xck+1-xck-∫EckEck+1x˙kdE‍∥+lkε3δpk+1(yk,dk)<p~k(xk,uk).
Thus, the sequence p~k(xk,uk), k=1,2,3,…. is monotonically decreasing. Considering the existence of a minimizer, p~k(xk,uk) converges. Thus, the following equation is established:
(36)limk→∞(lkε3δpk+1(yk,dk)-lkε1δpk+1(yk,dk))=0.
From (36), we derive
(37)limk→∞(lkδpk+1(yk,dk))=0.
If dk≠0, [xkuk]T is not the minimizer of penalty function, then δpk+1(yk,dk)≠0. Considering that lk is bounded away from zero, we derive
(38)limk→∞dk=0,limk→∞δpk+1(yk,dk)=0,limk→∞[x˙k-f(xk,uk)]=0,limk→∞(λ˙k+Hxk)=0,limk→∞Huk=0.
Thus, the solution sequence [xkuk]T generated by the algorithm will converge to the solution of the optimal control problem.


## 5. Smooth Approximation for Inequality Constraints

The nondifferentiable TPBVP from optimal control problem (11) with the inequality constraints brings difficulties for the proposed algorithm, which can only solve a smooth problem. Smooth approximations for nondifferentiable optimization problems have been studied in [5, 6, 34]. In this section we introduce the smooth approximation method for the inequality constraints. If the inequality constraint function is about the state and the guidance command
(39)C(x,u)≤0
the Hamiltonian for the optimal control problem (11) with inequality constraints (39) is given by
(40)H=g(x,u)+λTf(x,u)+μmax{C(x,u),0}.
However, the Hamiltonian (40) is a nondifferentiable function, which brings difficulty for the proposed algorithm. For the smooth Hamiltonian approximation, a differentiable function is substituted for max{C(x,u),0}. We denote the nondifferentiable function as follows:
(41)max{0,y}=∫-∞y‍w(z)dz,
where w(z) is the step function
(42)w(z)={1,ifz≥0,0,ifz<0.
Considering s(z,a) as the approximation function of w(z), the differentiable approximation to max{0,y} is obtained as
(43)fp(y,a)≈∫-∞y‍s(z,a)dz=y+1aIn(1+e-ay),
where s(z,a) is Sigmoid function
(44)s(z,a)=11+e-az,a>0.
The Properties of fp(y,a), a>0 are as follows [5]: (1)
fp(y,a) is n-times continuously differentiable for any positive integer n, with ∂fp(y,a)/∂y=1/(1+e-ay) and ∂2fp(y,a)/∂y2=ae-ay/(1+e-ay)2.(2)
fp(y,a) is strictly convex and strictly increasing on ℝ.(3)
fp(y,a)>max{0,y}.(4)
lima→∞fp(y,a)=max{0,y}.(5)
fp(y,a)>fp(y,a1) for a<a1, y∈ℝ. Substituting (43) into (40), the following differentiable Hamiltonian approximation is obtained:
(45)H^=g(x,u)+λTf(x,u)+μ[C+1aIn(1+e-aC)].
If the inequality constraint function is a qth-order state variable inequality constraint given by
(46)C(x)≤0.
If p times derivatives of C(x) are required before u appears explicitly in the result, the Hamiltonian with inequality constraints (46) is given by
(47)H=g(x,u)+λTf(x,u)+μC(q),
where
(48)μ=0ifC<0.
Considering the differentiable approximation to μ as follows:
(49)μ=μ^w(C)≈μ^11+e-aC
the following differentiable Hamiltonian approximation is obtained:
(50)H^=g(x,u)+λTf(x,u)+μ^11+e-aCC(q).
With the differentiable Hamiltonian approximation, the inequality constraints are considered in the proposed guidance algorithm.


## 6. Simulation Results

### 6.1. Numerical Calculation Based on Linear TPBVP Discretization

In this section, the linear TPVBP is transformed into linear equations using the Simpson rule, which corresponds to the three-point Newton-Cotes quadrature rule. The Simpson rule for differential equationdX/dE=y(E) is given by
(51)∫ab‍y(E)dE=b-a6[y(a)+4y(a+b2)+y(b)].
From (51), we obtain
(52)X(b)-X(a)=b-a6[y(a)+4y(a+b2)+y(b)].
We disperse the normalized energy E into N equal segments as [E0E1⋯EN], with E0=Eck and EN=Ef. The differential equation dX/dE=y(X,E) can be transformed into the following equation:
(53)X(Ei)-X(Ei-1)=Ei-Ei-16[y(Xi,Ei)(Xi+Xi-12,Ei+Ei-12)0000000000000+4y(Xi+Xi-12,Ei+Ei-12)0000000000000(Xi+Xi-12,Ei+Ei-12)+y(Xi-1,Ei-1)],i=1,2,3,…,N.
Based on (53), the linear TPBVP (21) can be transformed into linear equations described by
(54)AX=B.
The linear equations can be solved by Gaussian elimination.Finally, we can obtain the solutionΔx(E), Δu(E), and λ(E) of the linear TPBVP (21) from the solution of the linear equations through interpolating.


### 6.2. GHV Model

The proposed method is applied to the Generic Hypersonic Vehicle (GHV) to verify the effectiveness by flight simulation. In order to develop a complete set of aerodynamic coefficients, experimental longitudinal and lateral-directional aerodynamics were obtained for the GHV by using six Langley wind tunnels [35]. The aerodynamic database of the GHV is shown in [35]. The GHV model with scramjet propulsion system is used for the simulation, where Sref=334.73 m2 and the full-scale weight mini=110000 kg.According to the characterization of scramjet propulsion system described in [2, 27], the development of the propulsion system begins with the following equation:
(55)T=0.5ϕIsp(ϕ,Ma)ρvg0CT(α,Ma),
where ϕ is the throttle command which varies from 0 to 1. The coefficient CT depends on α and Ma as follows:
(56)CT=-0.0012Maα2+0.008α+1.4-0.1Ma.
The impulse Isp is given by the following equation:
(57)Isp=-81Ma2+1004Ma-52+100ϕ.


### 6.3. Simulation Results for Modeling and Initial State Error

Considering the high cost of the flight experimental test, researchers do flight numerical simulation, hardware-in-the-loop simulation and final flight experimental test to verify the guidance algorithm step by step. In this section, flight numerical simulation results for the GHV with modeling error and initial state error are given to show the effectiveness of the proposed guidance algorithm. For the modeling error caused by the wind tunnel experimental test, the aerodynamic coefficients bias is considered in the flight simulation. The simulations end once the final normalized energyEf is reached. The optimal guidance algorithm updates the trajectory every guidance cycle, which is 0.3 s. The initial and required final states of the air-breathing vehicle are shown in Table 1. The number of the discrete nodes N=10. The dynamic pressure inequality constraint is given by
(58)Q=12ρv2≤140000Pa.Table 1
Initial and required final states of the vehicle.
r
ini (m)
v
ini (m/s)
γ
ini (deg)
r
rf (m)
v
rf (m/s)
γ
rf (deg)
20000
1000
0
35000
2000
0To assess the capability of the proposed algorithm to deal with disturbance, the aerodynamic coefficients bias and the initial state error in the simulations are given by following cases:(1)
ΔCL=-0.1CL, ΔCD=+0.1CD.(2)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=+500 m.(3)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=-500 m.(4)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=+1000 m.(5)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=-1000 m.The simulation results of the optimal guidance algorithm are shown in Figures4, 5, 6, 7, and 8. The altitude, velocity, flight path angle, angle of attack, and dynamic pressure are, respectively, shown in Figures 4 to 8. As expected, the final state is achieved under aerodynamic bias and initial state error. Figure 8 shows that the state inequality constraint is satisfied. The dynamic pressure is less than 140000 Pa for different cases.Figure 4
The altitude of the simulation results.Figure 5
The velocity of the simulation results.Figure 6
The flight path angle of the simulation results.Figure 7
The angle of attack of the simulation results.Figure 8
The dynamic pressure of the simulation results.Table2 shows the terminal error and the fuel cost of the optimal guidance algorithm. The terminal accuracy is high for different cases. Moreover, the fuel cost of the optimal guidance algorithm is lower than that of the tracking guidance algorithm. A greater initial state error will result in more fuel savings. In summary, with lower fuel cost the results of the proposed optimal guidance algorithm satisfies all the equality and inequality constraints. The simulation results show the great potential for the final flight experimental test.Table 2
The simulation results for optimal guidance algorithm.
Case
Δ
r
f (m)
Δ
v
f (m/s)
Δ
γ
f (deg)
Fuel cost (kg)
Fuel cost (tracking method) (kg)
1
−28.20
0.136
0.20
4572.19
4574.62
2
8.59
−0.03
−0.15
4583.14
4628.39
3
−15.08
0.07
0.20
4584.74
4605.66
4
108.16
0.47
−0.23
4602.92
4776.00
5
−32.43
0.19
−0.16
4601.33
4688.94


## 6.1. Numerical Calculation Based on Linear TPBVP Discretization

In this section, the linear TPVBP is transformed into linear equations using the Simpson rule, which corresponds to the three-point Newton-Cotes quadrature rule. The Simpson rule for differential equationdX/dE=y(E) is given by
(51)∫ab‍y(E)dE=b-a6[y(a)+4y(a+b2)+y(b)].
From (51), we obtain
(52)X(b)-X(a)=b-a6[y(a)+4y(a+b2)+y(b)].
We disperse the normalized energy E into N equal segments as [E0E1⋯EN], with E0=Eck and EN=Ef. The differential equation dX/dE=y(X,E) can be transformed into the following equation:
(53)X(Ei)-X(Ei-1)=Ei-Ei-16[y(Xi,Ei)(Xi+Xi-12,Ei+Ei-12)0000000000000+4y(Xi+Xi-12,Ei+Ei-12)0000000000000(Xi+Xi-12,Ei+Ei-12)+y(Xi-1,Ei-1)],i=1,2,3,…,N.
Based on (53), the linear TPBVP (21) can be transformed into linear equations described by
(54)AX=B.
The linear equations can be solved by Gaussian elimination.Finally, we can obtain the solutionΔx(E), Δu(E), and λ(E) of the linear TPBVP (21) from the solution of the linear equations through interpolating.


## 6.2. GHV Model

The proposed method is applied to the Generic Hypersonic Vehicle (GHV) to verify the effectiveness by flight simulation. In order to develop a complete set of aerodynamic coefficients, experimental longitudinal and lateral-directional aerodynamics were obtained for the GHV by using six Langley wind tunnels [35]. The aerodynamic database of the GHV is shown in [35]. The GHV model with scramjet propulsion system is used for the simulation, where Sref=334.73 m2 and the full-scale weight mini=110000 kg.According to the characterization of scramjet propulsion system described in [2, 27], the development of the propulsion system begins with the following equation:
(55)T=0.5ϕIsp(ϕ,Ma)ρvg0CT(α,Ma),
where ϕ is the throttle command which varies from 0 to 1. The coefficient CT depends on α and Ma as follows:
(56)CT=-0.0012Maα2+0.008α+1.4-0.1Ma.
The impulse Isp is given by the following equation:
(57)Isp=-81Ma2+1004Ma-52+100ϕ.


## 6.3. Simulation Results for Modeling and Initial State Error

Considering the high cost of the flight experimental test, researchers do flight numerical simulation, hardware-in-the-loop simulation and final flight experimental test to verify the guidance algorithm step by step. In this section, flight numerical simulation results for the GHV with modeling error and initial state error are given to show the effectiveness of the proposed guidance algorithm. For the modeling error caused by the wind tunnel experimental test, the aerodynamic coefficients bias is considered in the flight simulation. The simulations end once the final normalized energyEf is reached. The optimal guidance algorithm updates the trajectory every guidance cycle, which is 0.3 s. The initial and required final states of the air-breathing vehicle are shown in Table 1. The number of the discrete nodes N=10. The dynamic pressure inequality constraint is given by
(58)Q=12ρv2≤140000Pa.Table 1
Initial and required final states of the vehicle.
r
ini (m)
v
ini (m/s)
γ
ini (deg)
r
rf (m)
v
rf (m/s)
γ
rf (deg)
20000
1000
0
35000
2000
0To assess the capability of the proposed algorithm to deal with disturbance, the aerodynamic coefficients bias and the initial state error in the simulations are given by following cases:(1)
ΔCL=-0.1CL, ΔCD=+0.1CD.(2)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=+500 m.(3)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=-500 m.(4)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=+1000 m.(5)
ΔCL=-0.1CL, ΔCD=+0.1CD, Δrini=-1000 m.The simulation results of the optimal guidance algorithm are shown in Figures4, 5, 6, 7, and 8. The altitude, velocity, flight path angle, angle of attack, and dynamic pressure are, respectively, shown in Figures 4 to 8. As expected, the final state is achieved under aerodynamic bias and initial state error. Figure 8 shows that the state inequality constraint is satisfied. The dynamic pressure is less than 140000 Pa for different cases.Figure 4
The altitude of the simulation results.Figure 5
The velocity of the simulation results.Figure 6
The flight path angle of the simulation results.Figure 7
The angle of attack of the simulation results.Figure 8
The dynamic pressure of the simulation results.Table2 shows the terminal error and the fuel cost of the optimal guidance algorithm. The terminal accuracy is high for different cases. Moreover, the fuel cost of the optimal guidance algorithm is lower than that of the tracking guidance algorithm. A greater initial state error will result in more fuel savings. In summary, with lower fuel cost the results of the proposed optimal guidance algorithm satisfies all the equality and inequality constraints. The simulation results show the great potential for the final flight experimental test.Table 2
The simulation results for optimal guidance algorithm.
Case
Δ
r
f (m)
Δ
v
f (m/s)
Δ
γ
f (deg)
Fuel cost (kg)
Fuel cost (tracking method) (kg)
1
−28.20
0.136
0.20
4572.19
4574.62
2
8.59
−0.03
−0.15
4583.14
4628.39
3
−15.08
0.07
0.20
4584.74
4605.66
4
108.16
0.47
−0.23
4602.92
4776.00
5
−32.43
0.19
−0.16
4601.33
4688.94


## 7. Conclusion

In this paper, we present an optimal guidance algorithm for air-breathing launch vehicle based on optimal trajectory correction. The proposed guidance algorithm corrects the reference trajectory in every guidance cycle to satisfy the optimality condition of the optimal feedback control problem. By linearizing the optimality condition, the linear TPBVP is obtained for optimal trajectory correction. The solution of the linear TPBVP is derived using the Simpson rule. The new trajectory is generated by a linear search for the step size of the solution. Smooth Hamiltonian approximation is implemented to the inequality constraints. The sufficiency condition for the global convergence of the guidance algorithm is given in this paper.Finally, simulations for 5 different cases of modeling and initial state error are presented. Compared with the tracking guidance method, the simulation results prove the low fuel cost and high precision of the proposed optimal guidance method.



---

*Source: 313197-2013-12-10.xml*