# Steady Glide Dynamic Modeling and Trajectory Optimization for High Lift-to-Drag Ratio Reentry Vehicle

Authors: Liang Yang; Wanchun Chen; Xiaoming Liu; Hao Zhou
Journal: International Journal of Aerospace Engineering
 (2016)

---

## Abstract

Steady glide trajectory optimization for high lift-to-drag ratio reentry vehicle is a challenge because of weakly damped trajectory oscillation. This paper aims at providing a steady glide trajectory using numerical optimal method. A new steady glide dynamic modeling is formulated via extending a trajectory-oscillation suppressing scheme into the three-dimensional reentry dynamics with a spherical and rotating Earth. This scheme comprehensively considers all factors acting on the flight path angle and suppresses the trajectory oscillation by regulating the vertical acceleration in negative feedback form and keeping the lateral acceleration invariant. Then, a study on steady glide trajectory optimization is carried out based on this modeling and pseudospectral method. Two examples with and without bank reversal are taken to evaluate the performance and applicability of the new method. A comparison with the traditional method is also provided to demonstrate its superior performance. Finally, the feasibility of the pseudospectral solution is verified by comparing the optimal trajectory with integral trajectory. The results show that this method not only is capable of addressing the case which the traditional method cannot solve but also significantly improves the computational efficiency. More importantly, it provides more stable and safe optimal steady glide trajectory with high precision.


---

## Body

## 1. Introduction

Entry guidance plays an important role in generating the steering command to guide the vehicle from its initial condition to reach the destination safely and accurately. In general, traditional reentry guidance is divided into two parts. The first part is the generation of a feasible reference trajectory. The second part is the tracking of this reference trajectory [1]. This paper focuses on generating a feasible steady glide reference trajectory, especially for high lift-to-drag ratio reentry vehicle, using numerical optimal method. Previous researches in reentry trajectory optimization are summarized as follows. Scott applied the Legendre Pseudospectral method into the trajectory optimization of reentry vehicles. In Josselyn and Ross’s work [2], covector mapping theorem of Legendre Pseudospectral method was used to verify the first-order optimality condition arising in the path constraints trajectory optimization. Rao and Clarke [3] also studied the problem of reentry trajectory optimization using Legendre Pseudospectral method. The key features of the optimal trajectory and quality of trajectory obtained from the Legendre Pseudospectral method were discussed. Jorris and Cobb and Zhao and Zhou [4, 5] employed the Gauss Pseudospectral method to optimize the 2D and 3D reentry trajectory for Common Aero Vehicle (CAV), in which waypoint and no-fly zone constraints were considered as inner-point constraints in optimal process. Rahimi et al. [6] applied the particle swarm optimization into spacecraft reentry optimization. High-order polynomials were used to approximate the angle of attack and bank angle in problem formulation. The coefficients of both polynomials were considered as input variables in optimal process. It should be noted that, because of not considering the trajectory-oscillation suppressing scheme, the optimal trajectories generated from above methods are naturally oscillatory. In steady glide, the heating rate will not change sharply and the steady state will greatly release the burden of control system. Therefore, steady glide trajectory is the best reference trajectory for the reentry guidance. Actually, quasi-equilibrium-glide condition (QEGC) is a well-known “soft” path constraint that makes the trajectory change monotonously. However, complicated reentry dynamics, especially for high lift-to-drag ratio vehicle, are so sensitive to the “soft” path constraint that it is very difficult for numerical optimization method to converge when QEGC is considered in optimal process. Generally speaking, this constraint is suitable for the trajectory planning in which one or two parameters are searched by secant method so as to generate a feasible trajectory [7–9]. Therefore, steady glide trajectory optimization for the high lift-to-drag ratio vehicle is always a challenge for numerical optimization.The objective of this paper is to investigate the steady glide dynamic modeling and trajectory optimization for the high lift-to-drag ratio vehicle. A new steady glide dynamic modeling is formulated by extending a trajectory-oscillation suppressing scheme, which is presented by Yu and Chen in [10], into three-dimensional reentry dynamics. Firstly, a special fight path angle which is able to keep the vehicle flying in a steady glide is calculated from the command angle of attack, command bank angle, and current state. Then, the trajectory oscillation is suppressed via regulating the longitudinal acceleration in negative feedback form and keeping the lateral acceleration invariant. It should be noted that the negative feedback signal is the deviation between the special flight path angle and actual flight path angle. Simulation result shows that this scheme performs well in suppressing trajectory oscillation and guides the vehicle into steady glide as soon as possible. Additionally, a study on steady glide trajectory optimization is investigated based on this new modeling. The derivatives of command angle of attack and bank angle are chosen as the control variables. And the performance index is the weighted squares sum of those derivatives. The limits on actual angle of attack and bank angle are considered as the path constraints. In fact, steady glide trajectory optimization is a typical optimal control problem whose solutions change rapidly in certain regions. Therefore, Hp-adaptive Gaussian quadrature collocation method [11], which performs well in dealing with this kind of problem, is chosen to transfer the optimal control problem into a standard nonlinear programming problem and solve it. The notable difference from the Yu’s method is that the scheme is suitable for the three-dimensional reentry dynamics. Another notable difference is that the scheme comprehensively considers all factors (including the derivatives of reference angle of attack and reference bank angle) acting on the flight path angle. That makes it easy to integrate into the motion dynamics by the choice of those derivatives as the control variables. Two classical numeric optimal examples (with and without bank reversal) are taken to evaluate the performance of the steady glide trajectory optimization. In order to demonstrate the superior performance in applicability and computational efficiency, a comparison with the traditional method is also provided. Furthermore, a comparison between optimal trajectory and integral trajectory is carried out to verify the feasibility of the pseudospectral solution. The results show that the new method not only significantly improves the computational efficiency of trajectory optimization since using fewer nodes will achieve a higher accuracy for the steady glide reentry trajectory but also has an extensive applicability in considering more final constraints even with the bank reversal. Most importantly, it is capable of providing more stable and safe optimal steady glide trajectory with high precision, which would be a better choice for tracking guidance.This paper is organized as follows: entry dynamics including entry trajectory constraints and vehicle model are described in Section2; three-dimensional reentry trajectory-oscillation suppressing scheme is presented in Section 3; steady glide dynamic modeling and trajectory optimization are presented in detail in Section 4.


## 2. Dynamics and Vehicle Description

### 2.1. Entry Dynamics

The 3-DOF point mass dynamics of the reentry vehicle over a spherical, rotating Earth is described as follows [12]:(1)
r
˙
=
V
sin
⁡
γ
,
(2)
θ
˙
=
V
cos
⁡
γ
sin
⁡
ψ
r
cos
⁡
φ
,
(3)
φ
˙
=
V
cos
⁡
γ
cos
⁡
ψ
r
,
(4)
V
˙
=
-
D
m
-
g
sin
⁡
γ
+
ω
2
r
cos
⁡
φ
sin
⁡
γ
cos
⁡
φ
-
cos
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
(5)
γ
˙
=
1
V
L
cos
⁡
σ
m
+
V
2
r
-
g
cos
⁡
γ
+
2
ω
V
cos
⁡
φ
sin
⁡
ψ
+
ω
2
r
cos
⁡
φ
cos
⁡
γ
cos
⁡
φ
+
sin
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
(6)
ψ
˙
=
1
V
L
sin
⁡
σ
m
+
V
2
r
cos
⁡
γ
sin
⁡
ψ
tan
⁡
φ
-
2
ω
V
cos
⁡
φ
tan
⁡
γ
cos
⁡
ψ
-
sin
⁡
φ
+
ω
2
r
cos
⁡
γ
sin
⁡
φ
cos
⁡
φ
sin
⁡
ψ
,where r is the radial distance from the center of the Earth to the vehicle. In the latter, h denotes the altitude. The radius of the Earth is 6378145 m. θ and φ are the longitude and latitude, respectively. V is the Earth relative velocity. γ is the flight path angle of the Earth relative velocity. ψ is the azimuth angle of the Earth relative velocity. m is the mass of the vehicle. g
=
μ
/
r
2 is the gravity acceleration, where μ is Earth’s gravitational constant. ω denotes the Earth self-rotation rate. The aerodynamic lift L and drag D are given as follows:(7)
L
=
1
2
ρ
V
2
C
l
S
r
e
f
,
D
=
1
2
ρ
V
2
C
d
S
r
e
f
,where ρ
=
ρ
0
exp
⁡
(
-
h
/
H
) is the atmospheric density, where ρ
0 is the standard atmospheric pressure from the sea level. S
r
e
f is the reference area of the vehicle. C
l and C
d are lift and drag coefficients which are dependent on the vehicle configuration.


### 2.2. Entry Trajectory Constraints

Typical entry trajectory inequality path constraints for hypersonic vehicle include(8)
k
ρ
V
3.15
≤
Q
˙
m
a
x
,
(9)
L
cos
⁡
α
+
D
sin
⁡
α
≤
n
m
a
x
,
(10)
1
2
ρ
V
2
≤
q
m
a
x
,where (8) is the heating rate at a stagnation point on the surface of the vehicle. Equation (9) is the aerodynamic acceleration in the body-normal direction. Equation (10) is the dynamic pressure. The heating limit value is Q
˙
m
a
x, the load factor limit value is n
m
a
x, and the dynamic pressure limit value is q
m
a
x. These are dependent on the vehicle configuration and mission. Those three constraints are considered to be “hard” constraints that should be enforced strictly.


### 2.3. Vehicle Description and Model Assumption

Common Aero Vehicle (CAV) is one of the most representative hypersonic entry vehicles with high lift-to-drag ratio. Relying on aerodynamic control, this vehicle is able to glide without power through the atmosphere. There are two types of CAV in the report of Phillips [13], the low-lift CAV and the high-lift CAV. The high-lift CAV, namely, CAV-H, is modeled here to extend the trajectory optimization. The weight of CAV-H is 907 kg, the area reference is 0.4839 m2, and the maximum lift-to-drag ratio is about 3.5. In order to make derivation analysis more intuitive and easier to follow, it is assumed that the lift and drag coefficients are only dependent on the angle of attack. They can be expressed in the form as(11)
C
l
=
k
l
1
α
+
k
l
2
,
C
d
=
k
d
1
α
2
+
k
d
2
α
+
k
d
3
,where k
l
1
=
0.04675, k
l
2
=
-
0
.
10568, k
d
1
=
0
.
000508, k
d
2
=
0.004228, and k
d
3
=
0.0161. Moreover, because the angle of attack having the maximum lift-to-drag ratio is about 10 deg., the scope of angle of attack is extended to [5°, 20°]. The scope of bank angle is limited within [−60°, 60°]. The limiting values of heating rate, dynamic pressure, and normal load factor are 400 W/cm2, 60 kpa, and 2 g, respectively.


## 2.1. Entry Dynamics

The 3-DOF point mass dynamics of the reentry vehicle over a spherical, rotating Earth is described as follows [12]:(1)
r
˙
=
V
sin
⁡
γ
,
(2)
θ
˙
=
V
cos
⁡
γ
sin
⁡
ψ
r
cos
⁡
φ
,
(3)
φ
˙
=
V
cos
⁡
γ
cos
⁡
ψ
r
,
(4)
V
˙
=
-
D
m
-
g
sin
⁡
γ
+
ω
2
r
cos
⁡
φ
sin
⁡
γ
cos
⁡
φ
-
cos
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
(5)
γ
˙
=
1
V
L
cos
⁡
σ
m
+
V
2
r
-
g
cos
⁡
γ
+
2
ω
V
cos
⁡
φ
sin
⁡
ψ
+
ω
2
r
cos
⁡
φ
cos
⁡
γ
cos
⁡
φ
+
sin
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
(6)
ψ
˙
=
1
V
L
sin
⁡
σ
m
+
V
2
r
cos
⁡
γ
sin
⁡
ψ
tan
⁡
φ
-
2
ω
V
cos
⁡
φ
tan
⁡
γ
cos
⁡
ψ
-
sin
⁡
φ
+
ω
2
r
cos
⁡
γ
sin
⁡
φ
cos
⁡
φ
sin
⁡
ψ
,where r is the radial distance from the center of the Earth to the vehicle. In the latter, h denotes the altitude. The radius of the Earth is 6378145 m. θ and φ are the longitude and latitude, respectively. V is the Earth relative velocity. γ is the flight path angle of the Earth relative velocity. ψ is the azimuth angle of the Earth relative velocity. m is the mass of the vehicle. g
=
μ
/
r
2 is the gravity acceleration, where μ is Earth’s gravitational constant. ω denotes the Earth self-rotation rate. The aerodynamic lift L and drag D are given as follows:(7)
L
=
1
2
ρ
V
2
C
l
S
r
e
f
,
D
=
1
2
ρ
V
2
C
d
S
r
e
f
,where ρ
=
ρ
0
exp
⁡
(
-
h
/
H
) is the atmospheric density, where ρ
0 is the standard atmospheric pressure from the sea level. S
r
e
f is the reference area of the vehicle. C
l and C
d are lift and drag coefficients which are dependent on the vehicle configuration.


## 2.2. Entry Trajectory Constraints

Typical entry trajectory inequality path constraints for hypersonic vehicle include(8)
k
ρ
V
3.15
≤
Q
˙
m
a
x
,
(9)
L
cos
⁡
α
+
D
sin
⁡
α
≤
n
m
a
x
,
(10)
1
2
ρ
V
2
≤
q
m
a
x
,where (8) is the heating rate at a stagnation point on the surface of the vehicle. Equation (9) is the aerodynamic acceleration in the body-normal direction. Equation (10) is the dynamic pressure. The heating limit value is Q
˙
m
a
x, the load factor limit value is n
m
a
x, and the dynamic pressure limit value is q
m
a
x. These are dependent on the vehicle configuration and mission. Those three constraints are considered to be “hard” constraints that should be enforced strictly.


## 2.3. Vehicle Description and Model Assumption

Common Aero Vehicle (CAV) is one of the most representative hypersonic entry vehicles with high lift-to-drag ratio. Relying on aerodynamic control, this vehicle is able to glide without power through the atmosphere. There are two types of CAV in the report of Phillips [13], the low-lift CAV and the high-lift CAV. The high-lift CAV, namely, CAV-H, is modeled here to extend the trajectory optimization. The weight of CAV-H is 907 kg, the area reference is 0.4839 m2, and the maximum lift-to-drag ratio is about 3.5. In order to make derivation analysis more intuitive and easier to follow, it is assumed that the lift and drag coefficients are only dependent on the angle of attack. They can be expressed in the form as(11)
C
l
=
k
l
1
α
+
k
l
2
,
C
d
=
k
d
1
α
2
+
k
d
2
α
+
k
d
3
,where k
l
1
=
0.04675, k
l
2
=
-
0
.
10568, k
d
1
=
0
.
000508, k
d
2
=
0.004228, and k
d
3
=
0.0161. Moreover, because the angle of attack having the maximum lift-to-drag ratio is about 10 deg., the scope of angle of attack is extended to [5°, 20°]. The scope of bank angle is limited within [−60°, 60°]. The limiting values of heating rate, dynamic pressure, and normal load factor are 400 W/cm2, 60 kpa, and 2 g, respectively.


## 3. Trajectory-Oscillation Suppressing Scheme

The objective of trajectory-oscillation suppressing scheme is to make the vehicle flying in a steady glide so that the heating rate will not change sharply, which also will significantly release the burden of control system. In this section, a trajectory-oscillation suppressing scheme using the flight path angle feedback is extended. Firstly, the command angle of attack and bank angle are used to calculate the special flight path angle that keeps the two-order derivative of flight path angle zero. Then, the trajectory oscillation is suppressed via regulating the longitudinal acceleration in negative feedback form and keeping the lateral acceleration invariant. The negative feedback signal is the deviation between the special flight path angle and actual flight path angle. Simulation results show that the scheme is capable of guiding the vehicle into a steady condition in which the command angle of attack and bank angle can generate enough vertical lift to sustain the vehicle glide.


### 3.1. Generic Theory for the Oscillation Suppressing Scheme

Let us pay attention to (5). While the flight path angle is small and varies relatively slow, it is assumed that cos
⁡
γ
=
1 and sin
⁡
γ
=
γ. Therefore, (5) without considering the earth rotation is formulated as(12)
γ
˙
=
1
V
L
cos
⁡
σ
m
+
V
2
r
-
g
,where L
=
1
/
2
ρ
V
2
C
l
S
r
e
f and g
=
μ
/
r
2. Substitute them into (12), then, (13) is the derivative of (12) with respect to time. Consider(13)
γ
¨
=
1
V
d
ρ
V
2
C
l
S
r
e
f
cos
⁡
σ
/
2
m
+
V
2
/
r
-
μ
/
r
2
d
t
+
D
m
+
g
γ
ρ
C
l
S
r
e
f
cos
⁡
σ
2
m
+
1
r
-
g
V
2
,where (14)
d
ρ
V
2
C
l
S
r
e
f
cos
⁡
σ
/
2
m
+
V
2
/
r
-
μ
/
r
2
d
t
=
V
2
C
l
S
r
e
f
cos
⁡
σ
2
m
ρ
˙
+
ρ
V
C
l
S
r
e
f
cos
⁡
σ
m
V
˙
+
ρ
V
2
S
r
e
f
2
m
C
˙
l
cos
⁡
σ
-
C
l
sin
⁡
σ
σ
˙
+
2
V
r
2
V
˙
-
V
2
r
2
r
˙
+
2
μ
r
3
r
˙
.Substitute (1), (4), and the derivative of the atmospheric density into (14). Then, substitute (14) into (13). The two-order derivative of flight path angle can be rewritten as(15)
γ
¨
=
-
ρ
V
2
C
l
S
r
e
f
cos
⁡
σ
2
m
H
γ
-
ρ
C
l
S
r
e
f
cos
⁡
σ
D
2
m
2
-
ρ
C
l
S
r
e
f
cos
⁡
σ
g
2
m
γ
+
ρ
V
S
r
e
f
2
m
C
˙
l
cos
⁡
σ
-
C
l
sin
⁡
σ
σ
˙
-
D
r
m
-
V
2
r
2
γ
-
D
g
m
V
2
+
g
r
γ
-
g
2
V
2
γ
,where the derivative of lift coefficient is formulated as follows:(16)
C
˙
l
=
∂
C
l
∂
α
·
∂
α
∂
t
+
∂
C
l
∂
M
a
·
∂
M
a
∂
r
·
∂
r
∂
t
.For reentry vehicle, the Mach number is much larger than five. Therefore, the aerodynamic coefficients are often assumed to be not dependent on the Mach number. Then, the last term in (16) can be neglected.Now, consider that the vehicle is in a certain condition of reentry process; if the command angle of attack and bank angle are given, the special flight path angle that keeps the two-order derivative of flight path angle zero can be formulated as follows:(17)
γ
m
=
D
∗
-
m
V
2
-
m
g
-
2
m
2
V
2
/
r
2
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
+
2
m
2
g
/
r
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
-
2
m
2
g
2
/
V
2
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
+
C
˙
l
∗
cos
⁡
σ
-
C
l
∗
sin
⁡
σ
σ
˙
V
C
l
∗
cos
⁡
σ
/
H
+
C
l
∗
cos
⁡
σ
g
/
V
+
2
m
V
/
r
2
ρ
S
r
e
f
-
2
m
g
/
r
ρ
V
S
r
e
f
+
2
m
g
2
/
V
3
ρ
S
r
e
f
+
D
∗
-
ρ
V
2
C
l
∗
S
r
e
f
cos
⁡
σ
r
/
2
H
-
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
g
r
/
2
-
V
2
/
r
m
+
m
g
-
g
2
r
m
/
V
2
+
D
∗
-
ρ
V
4
C
l
∗
S
r
e
f
cos
⁡
σ
/
2
H
g
-
ρ
V
2
C
l
∗
S
r
e
f
cos
⁡
σ
/
2
-
m
V
4
/
r
2
g
+
m
V
2
/
r
-
m
g
,where superscript ∗ denotes the lift and drag coefficients dominated by the command angle of attack and bank angle. The notable difference from Yu’s work is that the special fight path angle is derived from the three-dimensional reentry dynamics and comprehensively considers all factors including the derivatives of command angle of attack and bank angle. In the later section, the special flight path angle is easily calculated in the trajectory optimization when considering the derivatives of command angle of attack and bank angle as control variables. Then, trajectory oscillation is suppressed via regulating the longitudinal acceleration in negative feedback form and keeping the lateral acceleration invariant. The flight path angle also tends to the special flight path angle as soon as possible. Consider (18)
C
l
2
cos
⁡
σ
2
=
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
,
C
l
2
sin
⁡
σ
2
=
C
l
1
sin
⁡
σ
1
,where C
l
2 is the actual lift coefficient and C
l
1 is the command lift coefficient dominated by the command angle of attack. σ
2 is the actual bank angle. σ
1 is the command bank angle. K is a negative feedback gain; it should be noted that a better K will perform well in suppressing the trajectory oscillation. Then, a numerical simulation of entry process is carried out to evaluate the performance of the proposed scheme.


### 3.2. Performance of the Trajectory-Oscillation Suppressing Scheme

In this section, the proposed scheme is applied into reentry simulation with constant command angle of attack and bank angle. The simulation model, aerodynamic data, and some key parameters are stated in Section2.3. The command angle of attack is 10 deg. The command bank angle is also 10 deg. If the assumptions mentioned in Section 2.3 are held, the actual angle of attack and bank angle suppressing the trajectory oscillation can be calculated as follows: (19)
σ
2
=
arctan
⁡
C
l
1
sin
⁡
σ
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
,
α
2
=
1
k
l
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
cos
⁡
σ
2
-
k
l
2
.The initial states areh
0
=
60 km, θ
0
=
0 deg., φ
0
=
0 deg., V
0
=
7100 m/s, γ
0
=
0 deg., and ψ
0
=
90 deg. All programs run on a personal computer with a 3.3 Ghz processor and MATLAB 2008b. The solver of integral is ODE-45. The simulation stops when the altitude reduces to 30 km. Another worthy note is the negative feedback gain. After some attempts, it is easy to find that the trajectory oscillation will be suppressed perfectly when K is −16. Figure 1 shows the three-dimensional view of reentry trajectories. Note that the reference trajectory is the integral trajectory using the constant angle of attack and bank angle. It is obvious that the vehicle glides through the atmosphere in a great performance. Moreover, the glide trajectory is similar to the reference trajectory except that the trajectory oscillation is suppressed.Figure 1
Three-dimensional view of trajectory.Figures2 and 3 show the time histories of angle of attack and bank angle, respectively. The actual angle of attack and bank angle change substantially at the beginning of the flight and converge into the command angle of attack and bank angle quickly. Therefore, it is easy to conclude that applying the proposed scheme will guide the vehicle into a special condition in which the command angle of attack and bank angle will keep the vehicle steady glide. Another point should be noted is that the strict mathematical proof of convergence for this proposed method is absent. The future work will focus on this point.Figure 2
Histories of angle of attack and bank angle.
(a)
(b)Figure 3
Numeric results of Steady Glide Trajectory Optimization without Bank Reversal.
(a)
Altitude histories
(b)
Ground tracks
(c)
Velocity histories
(d)
Flight path angle histories
(e)
Azimuth angle histories
(f)
Azimuth angle histories
(g)
Bank angle histories
(h)
Heating rate histories
(i)
Dynamic pressure histories
(j)
Load factor histories


## 3.1. Generic Theory for the Oscillation Suppressing Scheme

Let us pay attention to (5). While the flight path angle is small and varies relatively slow, it is assumed that cos
⁡
γ
=
1 and sin
⁡
γ
=
γ. Therefore, (5) without considering the earth rotation is formulated as(12)
γ
˙
=
1
V
L
cos
⁡
σ
m
+
V
2
r
-
g
,where L
=
1
/
2
ρ
V
2
C
l
S
r
e
f and g
=
μ
/
r
2. Substitute them into (12), then, (13) is the derivative of (12) with respect to time. Consider(13)
γ
¨
=
1
V
d
ρ
V
2
C
l
S
r
e
f
cos
⁡
σ
/
2
m
+
V
2
/
r
-
μ
/
r
2
d
t
+
D
m
+
g
γ
ρ
C
l
S
r
e
f
cos
⁡
σ
2
m
+
1
r
-
g
V
2
,where (14)
d
ρ
V
2
C
l
S
r
e
f
cos
⁡
σ
/
2
m
+
V
2
/
r
-
μ
/
r
2
d
t
=
V
2
C
l
S
r
e
f
cos
⁡
σ
2
m
ρ
˙
+
ρ
V
C
l
S
r
e
f
cos
⁡
σ
m
V
˙
+
ρ
V
2
S
r
e
f
2
m
C
˙
l
cos
⁡
σ
-
C
l
sin
⁡
σ
σ
˙
+
2
V
r
2
V
˙
-
V
2
r
2
r
˙
+
2
μ
r
3
r
˙
.Substitute (1), (4), and the derivative of the atmospheric density into (14). Then, substitute (14) into (13). The two-order derivative of flight path angle can be rewritten as(15)
γ
¨
=
-
ρ
V
2
C
l
S
r
e
f
cos
⁡
σ
2
m
H
γ
-
ρ
C
l
S
r
e
f
cos
⁡
σ
D
2
m
2
-
ρ
C
l
S
r
e
f
cos
⁡
σ
g
2
m
γ
+
ρ
V
S
r
e
f
2
m
C
˙
l
cos
⁡
σ
-
C
l
sin
⁡
σ
σ
˙
-
D
r
m
-
V
2
r
2
γ
-
D
g
m
V
2
+
g
r
γ
-
g
2
V
2
γ
,where the derivative of lift coefficient is formulated as follows:(16)
C
˙
l
=
∂
C
l
∂
α
·
∂
α
∂
t
+
∂
C
l
∂
M
a
·
∂
M
a
∂
r
·
∂
r
∂
t
.For reentry vehicle, the Mach number is much larger than five. Therefore, the aerodynamic coefficients are often assumed to be not dependent on the Mach number. Then, the last term in (16) can be neglected.Now, consider that the vehicle is in a certain condition of reentry process; if the command angle of attack and bank angle are given, the special flight path angle that keeps the two-order derivative of flight path angle zero can be formulated as follows:(17)
γ
m
=
D
∗
-
m
V
2
-
m
g
-
2
m
2
V
2
/
r
2
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
+
2
m
2
g
/
r
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
-
2
m
2
g
2
/
V
2
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
+
C
˙
l
∗
cos
⁡
σ
-
C
l
∗
sin
⁡
σ
σ
˙
V
C
l
∗
cos
⁡
σ
/
H
+
C
l
∗
cos
⁡
σ
g
/
V
+
2
m
V
/
r
2
ρ
S
r
e
f
-
2
m
g
/
r
ρ
V
S
r
e
f
+
2
m
g
2
/
V
3
ρ
S
r
e
f
+
D
∗
-
ρ
V
2
C
l
∗
S
r
e
f
cos
⁡
σ
r
/
2
H
-
ρ
C
l
∗
S
r
e
f
cos
⁡
σ
g
r
/
2
-
V
2
/
r
m
+
m
g
-
g
2
r
m
/
V
2
+
D
∗
-
ρ
V
4
C
l
∗
S
r
e
f
cos
⁡
σ
/
2
H
g
-
ρ
V
2
C
l
∗
S
r
e
f
cos
⁡
σ
/
2
-
m
V
4
/
r
2
g
+
m
V
2
/
r
-
m
g
,where superscript ∗ denotes the lift and drag coefficients dominated by the command angle of attack and bank angle. The notable difference from Yu’s work is that the special fight path angle is derived from the three-dimensional reentry dynamics and comprehensively considers all factors including the derivatives of command angle of attack and bank angle. In the later section, the special flight path angle is easily calculated in the trajectory optimization when considering the derivatives of command angle of attack and bank angle as control variables. Then, trajectory oscillation is suppressed via regulating the longitudinal acceleration in negative feedback form and keeping the lateral acceleration invariant. The flight path angle also tends to the special flight path angle as soon as possible. Consider (18)
C
l
2
cos
⁡
σ
2
=
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
,
C
l
2
sin
⁡
σ
2
=
C
l
1
sin
⁡
σ
1
,where C
l
2 is the actual lift coefficient and C
l
1 is the command lift coefficient dominated by the command angle of attack. σ
2 is the actual bank angle. σ
1 is the command bank angle. K is a negative feedback gain; it should be noted that a better K will perform well in suppressing the trajectory oscillation. Then, a numerical simulation of entry process is carried out to evaluate the performance of the proposed scheme.


## 3.2. Performance of the Trajectory-Oscillation Suppressing Scheme

In this section, the proposed scheme is applied into reentry simulation with constant command angle of attack and bank angle. The simulation model, aerodynamic data, and some key parameters are stated in Section2.3. The command angle of attack is 10 deg. The command bank angle is also 10 deg. If the assumptions mentioned in Section 2.3 are held, the actual angle of attack and bank angle suppressing the trajectory oscillation can be calculated as follows: (19)
σ
2
=
arctan
⁡
C
l
1
sin
⁡
σ
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
,
α
2
=
1
k
l
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
cos
⁡
σ
2
-
k
l
2
.The initial states areh
0
=
60 km, θ
0
=
0 deg., φ
0
=
0 deg., V
0
=
7100 m/s, γ
0
=
0 deg., and ψ
0
=
90 deg. All programs run on a personal computer with a 3.3 Ghz processor and MATLAB 2008b. The solver of integral is ODE-45. The simulation stops when the altitude reduces to 30 km. Another worthy note is the negative feedback gain. After some attempts, it is easy to find that the trajectory oscillation will be suppressed perfectly when K is −16. Figure 1 shows the three-dimensional view of reentry trajectories. Note that the reference trajectory is the integral trajectory using the constant angle of attack and bank angle. It is obvious that the vehicle glides through the atmosphere in a great performance. Moreover, the glide trajectory is similar to the reference trajectory except that the trajectory oscillation is suppressed.Figure 1
Three-dimensional view of trajectory.Figures2 and 3 show the time histories of angle of attack and bank angle, respectively. The actual angle of attack and bank angle change substantially at the beginning of the flight and converge into the command angle of attack and bank angle quickly. Therefore, it is easy to conclude that applying the proposed scheme will guide the vehicle into a special condition in which the command angle of attack and bank angle will keep the vehicle steady glide. Another point should be noted is that the strict mathematical proof of convergence for this proposed method is absent. The future work will focus on this point.Figure 2
Histories of angle of attack and bank angle.
(a)
(b)Figure 3
Numeric results of Steady Glide Trajectory Optimization without Bank Reversal.
(a)
Altitude histories
(b)
Ground tracks
(c)
Velocity histories
(d)
Flight path angle histories
(e)
Azimuth angle histories
(f)
Azimuth angle histories
(g)
Bank angle histories
(h)
Heating rate histories
(i)
Dynamic pressure histories
(j)
Load factor histories


## 4. Steady Glide Dynamic Modeling and Trajectory Optimization

In this section, the steady glide trajectory optimization is carried out to find the optimal controls satisfying path and final constraints. Firstly, a new steady glide dynamic modeling is formulated via integrating the trajectory-oscillation suppressing scheme into the motion dynamics while the intrinsic properties remain. Secondly, the steady glide trajectory optimization is formulated based on this new modeling. The derivatives of command angle of attack and bank angle used to determine the special flight path angle are considered as control variables, and a performance index used to provide the stable controls is selected. Because steady glide trajectory optimization is a typical control problem whose solutions change rapidly in certain regions or are discontinuous, the Hp-adaptive Gaussian quadrature collocation method is naturally selected to transfer it into a nonlinear programming problem. Finally, two numeric examples with and without bank reversal are used to evaluate the performance and applicability of the new method. In order to demonstrate its superior performance in computational efficiency, a comparison with the traditional method is also provided. Moreover, a comparison between the optimal trajectory and integral trajectory is carried out to further verify the feasibility of the solution provided by the pseudospectral solver.


### 4.1. Steady Glide Dynamic Modeling

In previous section, the trajectory-oscillation suppressing scheme is presented in Section3.1. The performance of this scheme is evaluated by the simulation in Section 3.2. Now, the steady glide dynamic modeling is formulated via integrating the trajectory-oscillation suppressing scheme into the three-dimensional reentry dynamics in Section 2.1. Consider (20)
r
˙
=
V
sin
⁡
γ
;
θ
˙
=
V
cos
⁡
γ
sin
⁡
ψ
r
cos
⁡
φ
;
φ
˙
=
V
cos
⁡
γ
cos
⁡
ψ
r
,
V
˙
=
-
D
2
m
-
g
sin
⁡
γ
+
ω
2
r
cos
⁡
φ
sin
⁡
γ
cos
⁡
φ
-
cos
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
γ
˙
=
1
V
L
2
cos
⁡
σ
2
m
+
V
2
r
-
g
cos
⁡
γ
+
2
ω
V
cos
⁡
φ
sin
⁡
ψ
+
ω
2
r
cos
⁡
φ
cos
⁡
γ
cos
⁡
φ
+
sin
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
ψ
˙
=
1
V
L
2
sin
⁡
σ
2
m
+
V
2
r
cos
⁡
γ
sin
⁡
ψ
tan
⁡
φ
-
2
ω
V
cos
⁡
φ
tan
⁡
γ
cos
⁡
ψ
-
sin
⁡
φ
+
ω
2
r
cos
⁡
γ
sin
⁡
φ
cos
⁡
φ
sin
⁡
ψ
,
σ
2
=
atan
⁡
C
l
1
sin
⁡
σ
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
;
α
2
=
1
k
l
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
cos
⁡
σ
2
-
k
l
2
,
γ
m
=
f
r
,
V
,
α
1
,
σ
1
,
α
˙
1
,
σ
˙
1
;
α
˙
1
=
u
1
;
σ
˙
1
=
u
2
,where L
2 is the actual aerodynamic lift dominated by the command angle of attack, α
˙
1, σ
˙
1 are the derivatives of command angle of attack and bank angle, and α
1, σ
1, and σ
2 are defined in previous sections. If α
˙
1 and σ
˙
1 are given, there exists a specific steady glide trajectory determined by the steady glide dynamic modeling. Thus, steady glide trajectory optimization is continued based on this modeling in the next section so as to find the optimal control satisfying path and final constraints. Another point that should be noted is that there is existing analytical solution for the actual angle of attack and bank angle because of the simplification of aerodynamic coefficients. However, if complex aerodynamic coefficients are considered in the steady glide trajectory optimization, it only needs to regard the complex function of aerodynamic coefficients as the equality constraint in optimal process.


### 4.2. Hp-Adaptive Gaussian Quadrature Collocation Method

Hp-adaptive Gaussian quadrature collocation method is an efficient tool for solving multiple-phase optimal control problems using variable-order Gaussian quadrature collocation methods [14, 15]. Because an adaptive mesh refinement scheme [16, 17] is implemented to achieve a specified accuracy, this method performs well in solving the OCP whose solutions change rapidly in certain regions or are discontinuous. In each mesh interval, the dynamics and constraints are transferred into a set of nonlinear algebraic constraints by employing a Legendre-Gauss-Radau quadrature collocation method. Therefore, the continuous-time optimal control problem is transferred to a large sparse nonlinear programming problem that is solved using well-established techniques. This method takes the advantage of the exponential convergence in regions where the solution is smooth and places mesh points only near potential discontinuities or in regions where the solution changes rapidly.Steady glide trajectory optimization problem for high lift-drag ratio reentry vehicle is a complex constrained continuous-time optimal control problem. What is more, there exist several regions where the state and control variables change rapidly. At the beginning of the reentry, because the vehicle has no enough control capacity at high altitude, the altitude decreases rapidly in this phase. When the vehicle decreases to an appropriate height in which the vehicle has enough vertical lift to sustain the vehicle glide, the state and control variables also change rapidly. At the end of flying, it only takes a short time to steer the vehicle to meet the final requirements. The state and control variables also change rapidly. However, the vehicle glides in a steady and smooth condition for most of the flight duration. Naturally, Hp-adaptive Gaussian quadrature collocation method is suitable for solving the steady glide trajectory optimization problem. The later numerical results also show that Hp-adaptive Gaussian quadrature collocation method performs well in various glide trajectory optimizations.


### 4.3. Numerical Example of Trajectory Optimization without Bank Reversal

In this subsection, Hp-adaptive Gaussian quadrature collocation method is applied to solve the steady glide trajectory optimization problem. Meanwhile, a comparison with the traditional entry trajectory optimization using the dynamics of motion in (1)~(6) is also provided so as to further demonstrate the superior performance of the new method. Some of the simulation parameters are the same as those mentioned in above sections. For steady glide trajectory optimization, the motion dynamics are formulated in (18). In optimization process, the derivatives of command angle of attack and bank angle are chosen as the control variables, and the angle of attack and bank angle are considered as procedure variables. In order to make the trajectory state and optimal control as smooth and stable as possible, the performance index is selected as follows: (21)
J
=
∫
t
0
t
f
K
1
α
˙
1
2
+
K
2
σ
˙
1
2
d
t
,where K
1 and K
2 are control weighting coefficients. The limits on actual angle of attack and actual bank angle are considered as path constraints and presented as follows:(22)
σ
m
i
n
≤
arctan
⁡
C
l
1
sin
⁡
σ
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
≤
σ
m
a
x
,
α
m
i
n
≤
1
k
l
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
cos
⁡
σ
2
-
k
l
2
≤
α
m
a
x
.For the traditional trajectory optimization, the control variables are the derivatives of angle of attack and bank angle, and the performance index is the weighted square sum of those derivatives. The initial and final conditions for both methods are listed in Table 1. It should be noted that the traditional trajectory optimization will fail to converge due to considering too many final constraints. So the final longitude is set to be free. However, the new method is capable of addressing the problem with all final constraints. Therefore, the final longitude for such optimization is fixed at the optimal longitude provided by the traditional method. Another point that should be noted is that both optimal results come from the same initial guess.Table 1
Initial and final conditions for different methods.
Initial condition
θ
0 (deg.)
φ
0 (deg.)
h
0 (m)
V
0 (m/s)
γ
0 (deg.)
ψ
0 (deg.)
Traditional method
0
0
70000
6900
0
65
New method
0
0
70000
6900
0
65
Final condition
θ
f (deg.)
φ
f (deg.)
h
f (m)
V
f (m/s)
γ
f (deg.)
ψ
f (deg.)
Traditional method
—
0
30000
2400
0
—
New method
156.6
0
30000
2400
0
—In the procedure of optimization, finite difference method is used to provide the derivatives of nonlinear programming problem for both optimizations. Because the flight path angle is so sensitive, it is necessary to set the step of finite difference method to be a small value. The step is then set at 10−9. In addition, the feasible and optimal tolerances are set at 10−8, respectively.Numeric results for the case without bank reversal are presented in Figure3. As seen from those figures, it is apparent that a perfect steady glide optimal trajectory satisfying all constraints is provided via the new method, around which the entry trajectory provided by the traditional method shows some damped oscillation with a decreasing period. Moreover, the angle of attack and bank angle provided by the new method are very smooth and stable except at the beginning and end of the flight. The angle of attack and bank angle provided by the traditional method also appear to oscillate around that provided by the new method. It should be noted that because the angle of attack, bank angle, and the flight path angle are plotted over a large time scale (it is from 0 to 4000 sec), they seem to vary very sharply at the end of flight. In fact, those angles vary very slowly and completely meet the control requirements. The derivative of flight path angle is less than 1 deg./s. Another significant phenomenon depicted in Figures 3(h), 3(i), and 3(j) is that the path constraints generated by the traditional method are much larger than that provided by the new method due to the trajectory oscillation. It is obvious that the heating rate for the traditional method is much larger than the limit (the maximum is up to 430 W/cm2), but the one for the new method is within the safe range. The distribution of mesh nodes for both methods is shown in Figure 4. It is clear that the initial mesh nodes are the same for both methods. After several iterations, the mesh nodes for the traditional method are distributed with a high density, but the ones for the new method are distributed sparsely.Figure 4
Distribution of mesh nodes (blue: new method, red: traditional method).Table2 shows the statistics of optimal results by the NLP solver (SNOPT [18]). It is apparent that it only takes 4 times of mesh refinement to achieve the required feasible and optimal tolerances for both methods. However, the number of total nodes for the traditional method is much more than that for the new method, which is 498 for the traditional method but 285 for the new method. Moreover, because the time consumption will increase exponentially with the increase in the number of total nodes, the computational efficiency has been significantly improved for the new method. It only costs 82.06 sec to finish the optimization for the new method, which is only one-third of that for traditional method. Therefore, it is concluded that the new method not only is capable of solving the problem that the traditional method cannot solve but also has high computational efficiency since it will use less nodes to achieve a higher accuracy. Importantly, it will also provide more stable and safe optimal entry trajectory. It should be noted that it only takes half the time if “forward-difference” or “back-difference” is used in optimization. Applying more efficient resources is another way to reduce the computing time.Table 2
Optimal result by SNOPT.
Mesh iteration
Node
Feasible
Optimal
Maximum relative error for each mesh
Time consumption in each iteration
Time (sec)
Traditional method
1
41
9.2
E
-
15
6.8
E
-
9
0.1483900
37.17
222.04
2
287
3.4
E
-
15
3.2
E
-
9
1.8234
E
-
5
88.02
3
473
2.3
E
-
15
6.7
E
-
9
5.6194
E
-
7
82.58
4
498
2.3
E
-
15
7.5
E
-
9
9.3151
E
-
9
14.27
New method
1
41
5.1
E
-
15
3.7
E
-
9
0.0193110
1.61
82.06
2
209
2.5
E
-
15
4.5
E
-
9
0.0015244
24.91
3
278
2.2
E
-
15
8.7
E
-
9
3.9624
E
-
7
41.56
4
285
1.9
E
-
15
1.0
E
-
9
9.0488
E
-
9
13.98


### 4.4. Numerical Example of Trajectory Optimization with Bank Reversal

Generally speaking, the objective of steady glide trajectory optimization is to provide reference trajectory for traditional entry guidance such as linear-quadratic regulator tracking laws. Those laws focus on tracking only the reference longitudinal profile. The bank reversal is used to null the lateral errors. Meanwhile, reentry flight with bank reversal will significantly increase the trajectory shaping capability for the vehicle and provide more smooth control commands except when the bank reversal happened. Therefore, in order to provide the available reference trajectory, steady glide trajectory optimization with bank reversal is considered in this subsection. It is assumed that there is only one bank reversal during the reentry flight. Therefore, the optimal control problem is divided into two phases. And the signs of the bank angles before and after the bank reversal are opposite. The interior point constraints used to maintain continuity in the state at the phase boundaries are presented as(23)
t
+
=
t
-
,
r
t
+
=
r
t
-
,
θ
t
+
=
θ
t
-
,
φ
t
+
=
φ
t
-
,
v
t
+
=
v
t
-
,
γ
t
+
=
γ
t
-
,
ψ
t
+
=
ψ
t
-
,
α
t
+
=
α
t
-
,
σ
t
+
=
-
σ
t
-
,where t denotes the switching time for two phases. The performance index is stated in (21), and the motion dynamics are also stated in (20). All simulation environments are the same as that in Section 4.1. Initial states for glide vehicle are h
0
=
70000 m, θ
0
=
0 deg., φ
0
=
0 deg., V
0
=
7100 m/s, γ
0
=
0 deg., and ψ
0
=
90 deg. In addition, various cases for different cross-ranges are done to evaluate the applicability of the new method. The final longitude, latitude, velocity, and flight path angle are listed in Table 3.Table 3
Different terminal states for various cases.
Cases
Longitude (deg.)
Latitude (deg.)
Velocity (m/s)
Flight path angle (deg.)
Case1
155
0
2400
0
Case2
165
0
2400
0
Case3
175
0
2400
0
Case4
185
0
2400
0
Case5
195
0
2400
0Numeric results for the cases with bank reversal are presented in Figure5. Obviously, the new method is applicable for the cases with the bank reversal in different cross-ranges. All trajectories are optimal steady glide trajectories and meet all constraints. And it is very easy to find the bank reversal in the time history of bank angle. And another point that should be noted is that the magnitudes of the bank angle vary more stably and linearly than that without bank reversal. The meshes for various cases are placed perfectly in regions where the state and control change substantially. Additionally, the maximum heating rate in case 1 reaches the limit. The statistics for optimal results are presented in Table 4. It has been seen that the numbers of mesh iterations are 7 for case 1, case 2, and case 3, 4 for case 4, and 5 for case 5. The total node for various cases is about 190. All tolerances satisfy the stop criteria that are set at the beginning of optimization. From the statistics of CPU time, it only takes 2-3 min to successfully finish the steady glide optimization with high numerical precision.Table 4
Optimal results by SNOPT.
Mesh iterations
Total node
Feasible
Optimal
Maximum relative error for each mesh
Time (sec)
7
179
2.6
E
-
15
7.8
E
-
9
9.5388
E
-
7
170.77
7
188
2.5
E
-
15
9.2
E
-
8
6.8891
E
-
7
168.96
7
196
3.4
E
-
15
6.1
E
-
8
1.9260
E
-
7
204.70
4
192
3.2
E
-
15
6.4
E
-
8
8.0681
E
-
7
96.56
5
187
2.6
E
-
15
8.0
E
-
8
6.5696
E
-
7
135.68Figure 5
Numeric results of Steady Glide Trajectory Optimization with Bank Reversal.
(a)
Altitude histories
(b)
Ground tracks
(c)
Velocity histories
(d)
Flight path angle histories
(e)
Azimuth angle histories
(f)
Angle of attack histories
(g)
Bank angle histories
(h)
Heating rate histories
(i)
Dynamic pressure histories
(j)
Load factor histories


### 4.5. Verification of Feasibility for the Pseudospectral Solution

According to the results obtained above, it is not difficult to conclude that the proposed scheme not only performs well in suppressing the trajectory oscillation but also has high efficiency in three-dimensional steady glide optimization with and without bank reversal. However, there are three reasons that an elementary check should be conducted to validate the solution provided by the pseudospectral method. First, the feasible tolerance listed in Table3 is for NLP [19]. Second, according to the research in [20], Hp-adaptive Gaussian quadrature collocation method has some deficiency that it is unable to solve some special simple problems [20]. Third, the solution provided by the pseudospectral method meets dynamics only at a limited number of nodes. Therefore, a comparison between the optimal trajectories and integral trajectories is carried out to verify the feasibility of the resulting optimal trajectories.Figure6 shows the results comparing between the optimal trajectories provided by the pseudospectral solver and integral trajectories generated by integrating the equations of motion using the optimal controls presented in Section 4.3 with variable step (the min step size is 0.001 sec and the max step size is 0.1 sec). The comparison includes altitude, flight path angle, velocity, and ground footprint. It is apparent that those trajectories are practically the same even though the flight path angle varies within the range from −0.5 to 0.1 deg. The statistic on final errors for various integral trajectories is presented in Table 5. The maximal error for altitude is 2.5193 m, the maximal error for velocity is 0.1714 m/s, the maximal error for flight path angle is −0.003006 deg., the maximal error for longitude is 0.002742 deg., and the maximal error for latitude is 0.002329 deg. Obviously, the steady glide trajectory optimal control problem formulated in this paper can be solved by Hp-adaptive pseudospectral solver with a very high efficiency. And the solution is of high precision that it is suitable to be considered as the reference trajectory for the tracking reentry guidance law.Table 5
Final errors for various integral trajectories.
Cases
Altitude (m)
Velocity (m/s)
Flight path angle (deg.)
Longitude (deg.)
Latitude (deg.)
Case1
0.2858
0.1398
−0.001026
0.002742
0.002329
Case2
2.5193
0.1714
0.001847
0.002007
0.001763
Case3
1.3820
0.1538
−0.001858
0.002772
0.001763
Case4
2.1876
0.0058
0.001411
0.000452
0.001763
Case5
−1.1173
0.1686
−0.003006
0.002534
0.001690Figure 6
Comparison between numeric optimal trajectories and integral trajectories.
(a)
Altitude histories
(b)
Ground tracks
(c)
Velocity histories
(d)
Flight path angle histories


## 4.1. Steady Glide Dynamic Modeling

In previous section, the trajectory-oscillation suppressing scheme is presented in Section3.1. The performance of this scheme is evaluated by the simulation in Section 3.2. Now, the steady glide dynamic modeling is formulated via integrating the trajectory-oscillation suppressing scheme into the three-dimensional reentry dynamics in Section 2.1. Consider (20)
r
˙
=
V
sin
⁡
γ
;
θ
˙
=
V
cos
⁡
γ
sin
⁡
ψ
r
cos
⁡
φ
;
φ
˙
=
V
cos
⁡
γ
cos
⁡
ψ
r
,
V
˙
=
-
D
2
m
-
g
sin
⁡
γ
+
ω
2
r
cos
⁡
φ
sin
⁡
γ
cos
⁡
φ
-
cos
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
γ
˙
=
1
V
L
2
cos
⁡
σ
2
m
+
V
2
r
-
g
cos
⁡
γ
+
2
ω
V
cos
⁡
φ
sin
⁡
ψ
+
ω
2
r
cos
⁡
φ
cos
⁡
γ
cos
⁡
φ
+
sin
⁡
γ
sin
⁡
φ
cos
⁡
ψ
,
ψ
˙
=
1
V
L
2
sin
⁡
σ
2
m
+
V
2
r
cos
⁡
γ
sin
⁡
ψ
tan
⁡
φ
-
2
ω
V
cos
⁡
φ
tan
⁡
γ
cos
⁡
ψ
-
sin
⁡
φ
+
ω
2
r
cos
⁡
γ
sin
⁡
φ
cos
⁡
φ
sin
⁡
ψ
,
σ
2
=
atan
⁡
C
l
1
sin
⁡
σ
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
;
α
2
=
1
k
l
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
cos
⁡
σ
2
-
k
l
2
,
γ
m
=
f
r
,
V
,
α
1
,
σ
1
,
α
˙
1
,
σ
˙
1
;
α
˙
1
=
u
1
;
σ
˙
1
=
u
2
,where L
2 is the actual aerodynamic lift dominated by the command angle of attack, α
˙
1, σ
˙
1 are the derivatives of command angle of attack and bank angle, and α
1, σ
1, and σ
2 are defined in previous sections. If α
˙
1 and σ
˙
1 are given, there exists a specific steady glide trajectory determined by the steady glide dynamic modeling. Thus, steady glide trajectory optimization is continued based on this modeling in the next section so as to find the optimal control satisfying path and final constraints. Another point that should be noted is that there is existing analytical solution for the actual angle of attack and bank angle because of the simplification of aerodynamic coefficients. However, if complex aerodynamic coefficients are considered in the steady glide trajectory optimization, it only needs to regard the complex function of aerodynamic coefficients as the equality constraint in optimal process.


## 4.2. Hp-Adaptive Gaussian Quadrature Collocation Method

Hp-adaptive Gaussian quadrature collocation method is an efficient tool for solving multiple-phase optimal control problems using variable-order Gaussian quadrature collocation methods [14, 15]. Because an adaptive mesh refinement scheme [16, 17] is implemented to achieve a specified accuracy, this method performs well in solving the OCP whose solutions change rapidly in certain regions or are discontinuous. In each mesh interval, the dynamics and constraints are transferred into a set of nonlinear algebraic constraints by employing a Legendre-Gauss-Radau quadrature collocation method. Therefore, the continuous-time optimal control problem is transferred to a large sparse nonlinear programming problem that is solved using well-established techniques. This method takes the advantage of the exponential convergence in regions where the solution is smooth and places mesh points only near potential discontinuities or in regions where the solution changes rapidly.Steady glide trajectory optimization problem for high lift-drag ratio reentry vehicle is a complex constrained continuous-time optimal control problem. What is more, there exist several regions where the state and control variables change rapidly. At the beginning of the reentry, because the vehicle has no enough control capacity at high altitude, the altitude decreases rapidly in this phase. When the vehicle decreases to an appropriate height in which the vehicle has enough vertical lift to sustain the vehicle glide, the state and control variables also change rapidly. At the end of flying, it only takes a short time to steer the vehicle to meet the final requirements. The state and control variables also change rapidly. However, the vehicle glides in a steady and smooth condition for most of the flight duration. Naturally, Hp-adaptive Gaussian quadrature collocation method is suitable for solving the steady glide trajectory optimization problem. The later numerical results also show that Hp-adaptive Gaussian quadrature collocation method performs well in various glide trajectory optimizations.


## 4.3. Numerical Example of Trajectory Optimization without Bank Reversal

In this subsection, Hp-adaptive Gaussian quadrature collocation method is applied to solve the steady glide trajectory optimization problem. Meanwhile, a comparison with the traditional entry trajectory optimization using the dynamics of motion in (1)~(6) is also provided so as to further demonstrate the superior performance of the new method. Some of the simulation parameters are the same as those mentioned in above sections. For steady glide trajectory optimization, the motion dynamics are formulated in (18). In optimization process, the derivatives of command angle of attack and bank angle are chosen as the control variables, and the angle of attack and bank angle are considered as procedure variables. In order to make the trajectory state and optimal control as smooth and stable as possible, the performance index is selected as follows: (21)
J
=
∫
t
0
t
f
K
1
α
˙
1
2
+
K
2
σ
˙
1
2
d
t
,where K
1 and K
2 are control weighting coefficients. The limits on actual angle of attack and actual bank angle are considered as path constraints and presented as follows:(22)
σ
m
i
n
≤
arctan
⁡
C
l
1
sin
⁡
σ
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
≤
σ
m
a
x
,
α
m
i
n
≤
1
k
l
1
C
l
1
cos
⁡
σ
1
+
K
γ
-
γ
m
cos
⁡
σ
2
-
k
l
2
≤
α
m
a
x
.For the traditional trajectory optimization, the control variables are the derivatives of angle of attack and bank angle, and the performance index is the weighted square sum of those derivatives. The initial and final conditions for both methods are listed in Table 1. It should be noted that the traditional trajectory optimization will fail to converge due to considering too many final constraints. So the final longitude is set to be free. However, the new method is capable of addressing the problem with all final constraints. Therefore, the final longitude for such optimization is fixed at the optimal longitude provided by the traditional method. Another point that should be noted is that both optimal results come from the same initial guess.Table 1
Initial and final conditions for different methods.
Initial condition
θ
0 (deg.)
φ
0 (deg.)
h
0 (m)
V
0 (m/s)
γ
0 (deg.)
ψ
0 (deg.)
Traditional method
0
0
70000
6900
0
65
New method
0
0
70000
6900
0
65
Final condition
θ
f (deg.)
φ
f (deg.)
h
f (m)
V
f (m/s)
γ
f (deg.)
ψ
f (deg.)
Traditional method
—
0
30000
2400
0
—
New method
156.6
0
30000
2400
0
—In the procedure of optimization, finite difference method is used to provide the derivatives of nonlinear programming problem for both optimizations. Because the flight path angle is so sensitive, it is necessary to set the step of finite difference method to be a small value. The step is then set at 10−9. In addition, the feasible and optimal tolerances are set at 10−8, respectively.Numeric results for the case without bank reversal are presented in Figure3. As seen from those figures, it is apparent that a perfect steady glide optimal trajectory satisfying all constraints is provided via the new method, around which the entry trajectory provided by the traditional method shows some damped oscillation with a decreasing period. Moreover, the angle of attack and bank angle provided by the new method are very smooth and stable except at the beginning and end of the flight. The angle of attack and bank angle provided by the traditional method also appear to oscillate around that provided by the new method. It should be noted that because the angle of attack, bank angle, and the flight path angle are plotted over a large time scale (it is from 0 to 4000 sec), they seem to vary very sharply at the end of flight. In fact, those angles vary very slowly and completely meet the control requirements. The derivative of flight path angle is less than 1 deg./s. Another significant phenomenon depicted in Figures 3(h), 3(i), and 3(j) is that the path constraints generated by the traditional method are much larger than that provided by the new method due to the trajectory oscillation. It is obvious that the heating rate for the traditional method is much larger than the limit (the maximum is up to 430 W/cm2), but the one for the new method is within the safe range. The distribution of mesh nodes for both methods is shown in Figure 4. It is clear that the initial mesh nodes are the same for both methods. After several iterations, the mesh nodes for the traditional method are distributed with a high density, but the ones for the new method are distributed sparsely.Figure 4
Distribution of mesh nodes (blue: new method, red: traditional method).Table2 shows the statistics of optimal results by the NLP solver (SNOPT [18]). It is apparent that it only takes 4 times of mesh refinement to achieve the required feasible and optimal tolerances for both methods. However, the number of total nodes for the traditional method is much more than that for the new method, which is 498 for the traditional method but 285 for the new method. Moreover, because the time consumption will increase exponentially with the increase in the number of total nodes, the computational efficiency has been significantly improved for the new method. It only costs 82.06 sec to finish the optimization for the new method, which is only one-third of that for traditional method. Therefore, it is concluded that the new method not only is capable of solving the problem that the traditional method cannot solve but also has high computational efficiency since it will use less nodes to achieve a higher accuracy. Importantly, it will also provide more stable and safe optimal entry trajectory. It should be noted that it only takes half the time if “forward-difference” or “back-difference” is used in optimization. Applying more efficient resources is another way to reduce the computing time.Table 2
Optimal result by SNOPT.
Mesh iteration
Node
Feasible
Optimal
Maximum relative error for each mesh
Time consumption in each iteration
Time (sec)
Traditional method
1
41
9.2
E
-
15
6.8
E
-
9
0.1483900
37.17
222.04
2
287
3.4
E
-
15
3.2
E
-
9
1.8234
E
-
5
88.02
3
473
2.3
E
-
15
6.7
E
-
9
5.6194
E
-
7
82.58
4
498
2.3
E
-
15
7.5
E
-
9
9.3151
E
-
9
14.27
New method
1
41
5.1
E
-
15
3.7
E
-
9
0.0193110
1.61
82.06
2
209
2.5
E
-
15
4.5
E
-
9
0.0015244
24.91
3
278
2.2
E
-
15
8.7
E
-
9
3.9624
E
-
7
41.56
4
285
1.9
E
-
15
1.0
E
-
9
9.0488
E
-
9
13.98


## 4.4. Numerical Example of Trajectory Optimization with Bank Reversal

Generally speaking, the objective of steady glide trajectory optimization is to provide reference trajectory for traditional entry guidance such as linear-quadratic regulator tracking laws. Those laws focus on tracking only the reference longitudinal profile. The bank reversal is used to null the lateral errors. Meanwhile, reentry flight with bank reversal will significantly increase the trajectory shaping capability for the vehicle and provide more smooth control commands except when the bank reversal happened. Therefore, in order to provide the available reference trajectory, steady glide trajectory optimization with bank reversal is considered in this subsection. It is assumed that there is only one bank reversal during the reentry flight. Therefore, the optimal control problem is divided into two phases. And the signs of the bank angles before and after the bank reversal are opposite. The interior point constraints used to maintain continuity in the state at the phase boundaries are presented as(23)
t
+
=
t
-
,
r
t
+
=
r
t
-
,
θ
t
+
=
θ
t
-
,
φ
t
+
=
φ
t
-
,
v
t
+
=
v
t
-
,
γ
t
+
=
γ
t
-
,
ψ
t
+
=
ψ
t
-
,
α
t
+
=
α
t
-
,
σ
t
+
=
-
σ
t
-
,where t denotes the switching time for two phases. The performance index is stated in (21), and the motion dynamics are also stated in (20). All simulation environments are the same as that in Section 4.1. Initial states for glide vehicle are h
0
=
70000 m, θ
0
=
0 deg., φ
0
=
0 deg., V
0
=
7100 m/s, γ
0
=
0 deg., and ψ
0
=
90 deg. In addition, various cases for different cross-ranges are done to evaluate the applicability of the new method. The final longitude, latitude, velocity, and flight path angle are listed in Table 3.Table 3
Different terminal states for various cases.
Cases
Longitude (deg.)
Latitude (deg.)
Velocity (m/s)
Flight path angle (deg.)
Case1
155
0
2400
0
Case2
165
0
2400
0
Case3
175
0
2400
0
Case4
185
0
2400
0
Case5
195
0
2400
0Numeric results for the cases with bank reversal are presented in Figure5. Obviously, the new method is applicable for the cases with the bank reversal in different cross-ranges. All trajectories are optimal steady glide trajectories and meet all constraints. And it is very easy to find the bank reversal in the time history of bank angle. And another point that should be noted is that the magnitudes of the bank angle vary more stably and linearly than that without bank reversal. The meshes for various cases are placed perfectly in regions where the state and control change substantially. Additionally, the maximum heating rate in case 1 reaches the limit. The statistics for optimal results are presented in Table 4. It has been seen that the numbers of mesh iterations are 7 for case 1, case 2, and case 3, 4 for case 4, and 5 for case 5. The total node for various cases is about 190. All tolerances satisfy the stop criteria that are set at the beginning of optimization. From the statistics of CPU time, it only takes 2-3 min to successfully finish the steady glide optimization with high numerical precision.Table 4
Optimal results by SNOPT.
Mesh iterations
Total node
Feasible
Optimal
Maximum relative error for each mesh
Time (sec)
7
179
2.6
E
-
15
7.8
E
-
9
9.5388
E
-
7
170.77
7
188
2.5
E
-
15
9.2
E
-
8
6.8891
E
-
7
168.96
7
196
3.4
E
-
15
6.1
E
-
8
1.9260
E
-
7
204.70
4
192
3.2
E
-
15
6.4
E
-
8
8.0681
E
-
7
96.56
5
187
2.6
E
-
15
8.0
E
-
8
6.5696
E
-
7
135.68Figure 5
Numeric results of Steady Glide Trajectory Optimization with Bank Reversal.
(a)
Altitude histories
(b)
Ground tracks
(c)
Velocity histories
(d)
Flight path angle histories
(e)
Azimuth angle histories
(f)
Angle of attack histories
(g)
Bank angle histories
(h)
Heating rate histories
(i)
Dynamic pressure histories
(j)
Load factor histories


## 4.5. Verification of Feasibility for the Pseudospectral Solution

According to the results obtained above, it is not difficult to conclude that the proposed scheme not only performs well in suppressing the trajectory oscillation but also has high efficiency in three-dimensional steady glide optimization with and without bank reversal. However, there are three reasons that an elementary check should be conducted to validate the solution provided by the pseudospectral method. First, the feasible tolerance listed in Table3 is for NLP [19]. Second, according to the research in [20], Hp-adaptive Gaussian quadrature collocation method has some deficiency that it is unable to solve some special simple problems [20]. Third, the solution provided by the pseudospectral method meets dynamics only at a limited number of nodes. Therefore, a comparison between the optimal trajectories and integral trajectories is carried out to verify the feasibility of the resulting optimal trajectories.Figure6 shows the results comparing between the optimal trajectories provided by the pseudospectral solver and integral trajectories generated by integrating the equations of motion using the optimal controls presented in Section 4.3 with variable step (the min step size is 0.001 sec and the max step size is 0.1 sec). The comparison includes altitude, flight path angle, velocity, and ground footprint. It is apparent that those trajectories are practically the same even though the flight path angle varies within the range from −0.5 to 0.1 deg. The statistic on final errors for various integral trajectories is presented in Table 5. The maximal error for altitude is 2.5193 m, the maximal error for velocity is 0.1714 m/s, the maximal error for flight path angle is −0.003006 deg., the maximal error for longitude is 0.002742 deg., and the maximal error for latitude is 0.002329 deg. Obviously, the steady glide trajectory optimal control problem formulated in this paper can be solved by Hp-adaptive pseudospectral solver with a very high efficiency. And the solution is of high precision that it is suitable to be considered as the reference trajectory for the tracking reentry guidance law.Table 5
Final errors for various integral trajectories.
Cases
Altitude (m)
Velocity (m/s)
Flight path angle (deg.)
Longitude (deg.)
Latitude (deg.)
Case1
0.2858
0.1398
−0.001026
0.002742
0.002329
Case2
2.5193
0.1714
0.001847
0.002007
0.001763
Case3
1.3820
0.1538
−0.001858
0.002772
0.001763
Case4
2.1876
0.0058
0.001411
0.000452
0.001763
Case5
−1.1173
0.1686
−0.003006
0.002534
0.001690Figure 6
Comparison between numeric optimal trajectories and integral trajectories.
(a)
Altitude histories
(b)
Ground tracks
(c)
Velocity histories
(d)
Flight path angle histories


## 5. Conclusion

In this paper, a new steady glide dynamic modeling for high lift-to-drag ratio reentry vehicle is presented via extending the trajectory-oscillation suppressing scheme into three-dimensional reentry dynamics. Simulation shows that this scheme performs well in trajectory-oscillation suppressing. Based on this new modeling, a study on steady glide trajectory optimization with multiple final constraints is investigated. The derivatives of command angle of attack and bank angle are chosen as the control variables. And the weighted square sum of those derivatives is selected as the performance index so as to achieve the smoothest controls. Then, because the steady glide trajectory optimization is a typical optimal control problem whose solutions change rapidly in certain regions or are discontinuous, Hp-adaptive Gaussian quadrature collocation method is selected to solve it. Two classical optimization examples (with and without bank reversal) are conducted to show that the new method performs well in providing optimal steady glide trajectory even considering multiple final constraints. A comparison with the traditional method is done to demonstrate that the new method significantly improves the computational efficiency since the steady glide trajectory can be discretized with fewer nodes. And a comparison between optimal trajectory and integral trajectory is also carried out to verify that the solution provided by the pseudospectral method is of very high accuracy.



---

*Source: 3527460-2016-08-02.xml*