# Attainability Analysis for Entry Vehicles Based on Differential Evolution

Authors: Songtao Chang; Yongji Wang; Lei Liu; Xing Wei
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Attainable region provides crucial information on mission planning of entry vehicles. In order to obtain it, a series of
nonlinear optimal control problems which have similar formulations are needed to be solved. However, it is difficult
to compute due to severe nonlinearity of the dynamics and various constraints. In this paper, a novel method is
established to generate the attainable region at the end of the entry phase. It utilizes the parallel feature of differential
evolution (DE) and the high accuracy of Chebyshev polynomial interpolation. By using the Chebyshev polynomial
interpolation, the original problem is transformed to several nonlinear programming problems to facilitate employing
DE. Each individual in DE’s population represents a candidate point on the boundary of the attainable region. In order
to lead the population to the boundary simultaneously, a scheme is devised by exploiting the parallel feature of DE.
Different from conventional methods which generate one point of the boundary in each run, our proposed method
generates one side of the boundary of the attainable region. A scenario is presented to evaluate the designed method
and some analyses are conducted to evaluate the influence of the vehicle’s design parameters on the attainable region.


---

## Body

## 1. Introduction

Attainable region at the end of the entry phase (i.e., landing footprint) provides critical information on mission planning and evaluation of an entry vehicle's design. In order to obtain it, a series of nonlinear optimal control problems which have similar formulations are needed to be solved. The optimal goal is maximizing the terminal crossrange at each potential downrange under path constraints and terminal constraints. However, it is challenging to get accurate solutions of these problems due to severe nonlinearity of the vehicle's dynamics and various constraints. Some researchers derived approximate solution of the original problems by simplifying the vehicle's dynamics [1–3], in which path constraints were not considered, and the Earth's rotation was omitted. Therefore, the solution was inadequately accurate and poorly reliable. Saraf et al. [4] presented a rapid landing footprint computation algorithm, which was built around acceleration-based trajectory planner. Meanwhile, a near-optimal angle of attack control law was developed. But the footprint was inadequately accurate and the computed maximum downrange was much smaller than the one that the vehicle could achieve. Lu and Xue [5] developed a near-optimal closed-form control law by using quasiequilibrium glide condition (QEGC) [6]. It transformed the original problem to several closest-approach problems which were univariate root-finding ones. Path constraints were enforced through inflicting limitations on bank angle. However, only bank angle was optimized, and angle of attack was predetermined. Benito and Mease [7] used a direct method to solve the entry trajectory optimization problem. Direct method [8–11] transcribed the original optimal control problem to nonlinear programming problem (NLP) which could be solved easily. The method had the capacity to deal with path constraints and terminal constraints simultaneously. The accuracy of the solution can be improved by adjusting some parameters used in the transcription.It is noted that all of the aforementioned methods use determined optimization algorithm (e.g., sequential quadratic programming [12]) to solve the transcribed nonlinear programming problem. The performances of conventional optimization algorithms are highly affected by the initial guess of the solution since they are local optimization algorithms. To overcome this deficiency, some metaheuristic optimization algorithms [13–16] were proposed to approximate the global optimal solution. They were often inspired by some natural phenomena, such as genetic algorithm (GA) [13] which mimicked the process of natural selection, particle swarm optimization (PSO) [14] which was inspired from the movement of organisms like a flock of birds or a school of fishes, differential evolution (DE) [15] which mimicked the process of natural evolution, and simulated annealing (SA) [16] which was inspired by annealing in metallurgy. They made few or no assumptions on the problem to be optimized, could search very large spaces of candidate solutions, and did not need to have a good initial guess of the solution. Due to these advantages, they were widely applied in many fields, such as controller design [17], damage detection [18, 19], economic dispatch problem [20], and water distribution system optimization [21].In this paper, a parallel method is developed to calculate the landing footprint of entry vehicles based on DE and Chebyshev polynomial interpolation. The original optimal control problems are transcribed to NLPs by Chebyshev polynomial interpolation. After that, an improved DE is employed with each agent representing a candidate point on the boundary of the landing footprint. A new selection strategy is designed to drive the agents evolving to the boundary of the landing footprint simultaneously. It does not need to provide a good initial guess of the solution and can generate one side rather than one point of the boundary of the landing footprint by exploiting the parallel feature of DE. Furthermore, it does not rely on simplification of the system and considers all control signals which make it have high fidelity. A scenario is presented to evaluate the method. Through the analysis of the dynamic equations of the vehicle, design configuration can be translated to the lift and drag coefficients equivalently. In order to evaluate the effect of the vehicle's design parameters on the attainable region, various lift and drag configurations are used to generate the landing footprints.The remainder of this paper is organized as follows. Formulation of the landing footprint problem is given in Section2. Landing footprint computation algorithm is established in Section 3. In Section 4, a scenario is presented to validate the algorithm and some analyses are conducted. Finally, some remarks are concluded in Section 5.


## 2. Problem Formulation

Three degrees of freedom dynamic equations of an entry vehicle are as follows:(1)
r
˙
=
v
sin
γ
,
(2)
θ
˙
=
v
cos
⁡
γ
sin
ψ
r
cos
⁡
ϕ
,
(3)
ϕ
˙
=
v
cos
⁡
γ
cos
⁡
ψ
r
,
(4)
v
˙
=
ω
2
r
cos
⁡
ϕ
(
sin
γ
cos
⁡
ϕ
-
cos
⁡
γ
sin
ϕ
cos
⁡
ψ
)
-
D
-
g
sin
γ
,
(5)
γ
˙
=
L
cos
⁡
σ
v
-
g
v
cos
⁡
γ
+
v
r
cos
⁡
γ
+
2
ω
cos
⁡
ϕ
sin
ψ
+
ω
2
r
cos
⁡
ϕ
(
cos
⁡
γ
cos
⁡
ϕ
+
sin
γ
cos
⁡
ψ
sin
ϕ
)
v
,
(6)
ψ
˙
=
L
sin
σ
v
cos
⁡
γ
+
v
r
cos
⁡
γ
sin
ψ
tan
ϕ
-
2
ω
(
tan
γ
cos
⁡
ψ
cos
⁡
ϕ
-
sin
ϕ
)
+
ω
2
r
v
cos
⁡
γ
sin
ψ
sin
ϕ
cos
⁡
ϕ
,
(7)
g
=
μ
r
2
,
ρ
=
ρ
0
e
-
(
r
-
r
0
)
/
h
s
,
L
=
ρ
v
2
S
C
L
(
2
m
)
,
D
=
ρ
v
2
S
C
D
(
2
m
)
,
where r is the radial distance from the earth center to the vehicle, θ and ϕ are the longitude and latitude, respectively, v is the earth-relative velocity, γ is the relative flight-path angle, ψ is the relative velocity heading angle measured clockwise from the north, ω is the earth self-rotation rate, m is the vehicle's mass, S is the reference area, ρ is the density of atmosphere, ρ
0 is the density of atmosphere at sea level, r
0 is the earth's radius, C
L and C
D are the lift and drag coefficients determined by angle of attack α and Mach M
a, g is the gravity acceleration, independent variable is the time t, and control variables are angle of attack α and bank angle σ. The geometric sketch of entry flight is illustrated by Figure 1.Figure 1
Geometric sketch of entry flight.Entry process must satisfy constraints of heating rate, dynamic pressure, and normal acceleration:(8)
Q
˙
=
c
ρ
v
k
q
≤
Q
˙
max
⁡
.
(9)
q
=
ρ
v
2
2
≤
q
max
⁡
.
(10)
n
=
L
cos
⁡
α
+
D
sin
α
≤
n
max
⁡
,
where Q
˙
max
⁡, q
max
⁡, and n
max
⁡ represent maximum allowable heating rate, dynamic pressure, and acceleration, respectively.In order to facilitate the attitude controller design, entry trajectory should not oscillate acutely and should satisfy QEGC as follows:(11)
L
cos
⁡
σ
-
g
+
v
2
r
≈
0
.Following the entry phase is Terminal Area Energy Management (TAEM) phase. Terminal state must be restrained as follows:(12)
|
r
(
t
f
)
-
r
f
|
≤
Δ
r
,
|
v
(
t
f
)
-
v
f
|
≤
Δ
v
,
|
γ
(
t
f
)
-
γ
f
|
≤
Δ
γ
,
where t
f is the terminal time of the entry phase and r
f, v
f, and γ
f represent the ideal terminal states of the entry phase, respectively.Constraints of heating rate, dynamic pressure, normal acceleration, and terminal states are hard constraints, which must be satisfied. QEGC is a soft constraint which may be violated a little.Landing footprint of an entry vehicle is the attainable region when it reaches TAEM. It can be formulated to a series of optimal control problems. Each has prescribed terminal downrange (the range along the norm-plan) and the optimal objective is maximizing or minimizing terminal crossrange (the range perpendicular to the norm-plan).LetR
(
t
f
), Z
(
t
f
) be the downrange and the crossrange at the end of the entry phase, respectively. To calculate the landing footprint of an entry vehicle, we have to solve (2
N
+
2) optimal control problems whose objectives are described as follows:(13)
R
min
⁡
=
min
⁡
(
J
1
)
=
min
⁡
(
R
(
t
f
)
)
.
R
max
⁡
=
max
⁡
(
J
1
)
=
max
⁡
(
R
(
t
f
)
)
.
min
⁡
(
J
2
)
=
min
⁡
(
Z
(
t
f
)
)
subject
to
R
(
t
f
)
=
R
min
⁡
+
k
R
max
⁡
-
R
min
⁡
N
-
1
,
k
=
0,1
,
…
,
N
-
1
.
max
⁡
(
J
2
)
=
max
⁡
(
Z
(
t
f
)
)
subject
to
R
(
t
f
)
=
R
min
⁡
+
k
R
max
⁡
-
R
min
⁡
N
-
1
,
k
=
0,1
,
…
,
N
-
1
,
where N is the number of the calculated sites to approximate the lateral boundary:
(14)
R
(
t
f
)
=
r
(
t
f
)
·
arccos
r
b
(
t
f
)
·
r
i
|
r
b
(
t
f
)
|
·
|
r
i
|
,
Z
(
t
f
)
=
r
(
t
f
)
·
(
π
2
-
arccos
r
(
t
f
)
·
n
|
r
(
t
f
)
|
·
|
n
|
)
,
r
b
(
t
f
)
=
r
(
t
f
)
-
(
r
(
t
f
)
·
n
)
·
n
,
r
i
=
r
i
·
[
cos
⁡
ϕ
i
cos
⁡
θ
i
,
cos
⁡
ϕ
i
sin
θ
i
,
sin
ϕ
i
]
,
r
(
t
f
)
=
r
(
t
f
)
·
[
cos
⁡
ϕ
(
t
f
)
cos
⁡
θ
(
t
f
)
,
c
c
c
c
c
c
c
i
cos
⁡
ϕ
(
t
f
)
sin
θ
(
t
f
)
,
sin
ϕ
(
t
f
)
]
,
n
=
v
i
×
r
i
|
v
i
×
r
i
|
,
v
i
=
v
i
·
[
cos
⁡
θ
i
cos
⁡
ϕ
i
sin
γ
i
-
sin
θ
i
cos
⁡
γ
i
sin
ψ
i
ccici
-
cos
⁡
θ
i
sin
ϕ
i
cos
⁡
γ
i
cos
⁡
ψ
i
,
sin
θ
i
cos
⁡
ϕ
i
sin
γ
i
ccici
+
cos
⁡
θ
i
cos
⁡
γ
i
sin
ψ
i
-
sin
θ
i
sin
ϕ
i
cos
⁡
γ
i
cos
⁡
ψ
i
,
ccici
sin
ϕ
i
sin
γ
i
+
cos
⁡
ϕ
i
cos
⁡
γ
i
cos
⁡
ψ
i
]
,
where r
i
,
θ
i
,
ϕ
i
,
v
i
,
γ
i, and ψ
i are the states at the initial time of the entry phase.


## 3. Landing Footprint Computation Algorithm

In order to calculate the landing footprint of an entry vehicle, we must solve a series of multiconstraints optimal control problems whose objectives are described as (13), state equations are described as (1)–(6), and constraints are described as (8)–(12). These optimal control problems have similar formulations; the only difference among them is the way to deal with the terminal crossrange and downrange. By applying this characteristic and the parallel feature of DE, an algorithm is designed to solve half of these problems (maximizing or minimizing) simultaneously. First, the angle of attack's scope is determined by path constraints to simplify the problem. And the QEGC constraint (11) is guaranteed by inflicting the low bound of the angle of attack. Then, control profiles are parameterized by Chebyshev polynomial interpolation to facilitate employing DE algorithm. After that, DE is given. At last, an improved DE is applied to calculate the landing footprint. Each run of the algorithm generates one side of the landing footprint. To get the full landing footprint, it just needs to run the algorithm twice.


### 3.1. Discretization of the Control Profiles

From (8)-(9) and the expression of density ρ in (7), we have
(15)
r
≥
r
0
-
2
h
s
ln
⁡
(
Q
˙
max
⁡
ρ
0
c
v
k
q
)
≜
r
Q
(
v
)
,
r
≥
r
0
-
h
s
ln
⁡
2
q
max
⁡
ρ
0
v
2
≜
r
q
(
v
)
.
Let y
(
α
,
v
,
r
,
σ
)
=
L
cos
⁡
σ
-
g
+
v
2
/
r
=
(
S
/
2
m
)
ρ
v
2
cos
⁡
σ
C
L
(
α
,
v
,
r
)
-
μ
/
r
2
+
v
2
/
r, according to (11), y
≈
0. Solve the equation y
=
0 with dependent variable being α; we have α
=
α
(
v
,
r
,
σ
),
(16)
∂
y
∂
σ
=
-
L
sin
σ
,
∂
y
∂
α
=
ρ
v
2
S
2
m
∂
C
L
∂
α
cos
⁡
σ
,
∂
y
∂
r
=
2
g
r
-
v
2
r
2
-
L
h
s
cos
⁡
σ
+
L
C
L
∂
C
L
∂
r
cos
⁡
σ
,
σ
∈
[
-
9
0
∘
,
9
0
∘
]
,
L
>
0; thus, ∂
y
/
∂
σ
>
0 when σ
∈
[
-
9
0
∘
,
0
∘
) and ∂
y
/
∂
σ
<
0 when σ
∈
(
0
∘
,
9
0
∘
]. Lift coefficient C
L increases as α increases; that is, ∂
y
/
∂
α
>
0. Consider ∂
α
(
v
,
r
,
σ
)
/
∂
σ
=
-
(
∂
y
/
∂
σ
)
/
(
∂
y
/
∂
α
); then, ∂
α
(
v
,
r
,
σ
)
/
∂
σ
<
0 when σ
∈
[
-
9
0
∘
,
0
∘
) and ∂
α
(
v
,
r
,
σ
)
/
∂
σ
>
0 when σ
∈
(
0
∘
,
9
0
∘
]; therefore, α
(
v
,
r
,
σ
)
≥
α
(
v
,
r
,
0
).∂
C
L
/
∂
r and 2
g
/
r are significantly smaller than L
/
h
s, so they can be ignored; then ∂
y
/
∂
r
<
0, ∂
α
(
v
,
r
,
σ
)
/
∂
r
=
-
(
∂
y
/
∂
r
)
/
(
∂
y
/
∂
α
)
>
0; therefore
(17)
α
=
α
(
v
,
r
,
σ
)
≥
α
(
v
,
r
,
0
)
≥
α
(
v
,
max
⁡
(
r
Q
(
v
)
,
r
q
(
v
)
)
,
0
)
≜
α
min
⁡
(
v
)
.
Suppose the maximum angle of attack and bank angle are α
max
⁡ and σ
max
⁡, respectively, and the minimum bank angle is σ
min
⁡. Chebyshev interpolation polynomial is an optimal approximation of the original function in all interpolation polynomials in L
∞ norm-sense [22]. Herein, control profiles are approximated by Chebyshev interpolation polynomials as follows:
(18)
η
(
τ
)
=
∑
i
=
1
N
1
‍
η
i
Φ
i
(
τ
)
,
ξ
(
τ
)
=
∑
i
=
1
N
2
‍
ξ
i
Υ
i
(
τ
)
.
α
(
v
)
=
(
1
-
η
(
τ
)
)
α
min
⁡
(
v
)
+
η
(
τ
)
α
max
⁡
,
σ
(
v
)
=
(
1
-
ξ
(
τ
)
)
σ
min
⁡
+
ξ
(
τ
)
σ
max
⁡
,
where τ
=
2
(
v
-
v
0
)
/
(
v
f
-
v
0
)
-
1
∈
[
-
1,1
], Φ
i
(
τ
) is the ith function of basis functions of N
1th order Chebyshev polynomial interpolation, and Υ
i
(
τ
) is the ith function of basis functions of N
2th order Chebyshev polynomial interpolation. N
1 and N
2 are the orders of the interpolation polynomials of α and σ, respectively.LetX
=
[
η
1
,
η
2
,
…
,
η
N
1
,
ξ
1
,
ξ
2
,
…
,
ξ
N
2
]
,
η
k
,
ξ
k
∈
[
0,1
], and X is the decision vector. When X is known, control profiles are determined by (18); then we can integrate (1)–(6) to figure out the objective function (13), and the constraints (8)–(10) and (12). The nonlinear programming problems are well defined.


### 3.2. Differential Evolution

DE is a simple, efficient, and robust evolutionary algorithm that deals with unconstrained parameter optimization problem [15]. Unconstrained parameter optimization can be formulated to minimize a single D-dimensional function f
(
X
)
,
X
=
[
x
1
,
x
2
,
…
,
x
D
], x
i
∈
[
x
l
,
i
,
x
u
,
i
]. Let P
G
=
{
X
G
,
1
,
X
G
,
2
,
…
,
X
G
,
NP
} represent the Gth generation of the population, where the ith individual's position is X
G
,
i
=
[
x
G
,
i
,
1
,
x
G
,
i
,
2
,
…
,
x
G
,
i
,
D
], i
=
1,2
,
…
,
NP; NP is the population's scale.DE involves two stages: initialization and evolution. Initialization stage uses (19) to generate the initial population P
0 randomly. Then P
0 evolves to P
1, P
1 evolves to P
2
,
…, and so on, until the terminal conditions are reached. Each evolving phase includes three operations, namely, differential mutation, crossover, and selection:
(19)
x
0
,
i
,
j
=
x
l
,
j
+
rand
[
0,1
]
×
(
x
u
,
j
-
x
l
,
j
)
,
i
=
1,2
,
…
,
NP
,
j
=
1,2
,
…
,
D
,
where rand
[
0,1
] is a uniform random number in [
0,1
].Mutation and crossover create offsprings of the population, and the selection operation controls the direction of evolution. The following will describe these operations.LetV
G
,
i
=
[
v
G
,
i
,
1
,
v
G
,
i
,
2
,
…
,
v
G
,
i
,
D
] denote the ith mutant individual. Classical mutation operator is described as follows:
(20)
V
G
,
i
=
X
G
,
r
1
+
F
×
(
X
G
,
r
2
-
X
G
,
r
3
)
,
where r
1
,
r
2, and r
3 are mutually exclusive indices randomly chosen in the range [
1
,
NP
], which are different from the base index i. F
>
0 is the mutant factor.The mutant individualV
G
,
i exchanges its components with the base vector X
G
,
i through crossover operation to get the trial vector U
G
,
i
=
[
u
G
,
i
,
1
,
u
G
,
i
,
2
,
…
,
u
G
,
i
,
D
]. Common crossover operation is described as follows:
(21)
u
G
,
i
,
j
=
{
v
G
,
i
,
j
,
j
=
j
rand
or
rand

[
0,1
]
≤
C
R
x
G
,
i
,
j
,
otherwise
,
where j
rand is a random index in the range [
1
,
D
]. It ensures that the trial vector has at least one component from the mutant vector. C
R
>
0 is the crossover probability.By using one to one greedy selection operation to decide whether the trial vectorU
G
,
i substitutes X
G
,
i to the next iteration. This strategy guarantees that the population evolves to superior region.


### 3.3. Landing Footprint Generation Algorithm

When coming to calculate the landing footprint, the algorithm must guarantee that none of the constraints is violated, and it should have extreme crossrange and maintains diversity in downrange distribution. DE uses fitness to drive the evolving of the population. Here, we design a scheme to determine the fitness of each agent to lead the agent to the boundary of the landing footprint.LetT
G represent the population constituted by P
G and its offspring. At the first stage, sort all individuals in T
G according to the violation magnitude of the constraints (magnitude of the normalized constraints vector) in ascending order. Suppose it has M sets that are denoted as Ϝ
k, k
=
1,2
,
…
,
M. All individuals who do not violate any constraint (these individuals also called feasible individuals) have the highest priority, and they are denoted as Ϝ
1. At the second stage, downrange and crossrange are used to classify Ϝ
1. Suppose it has L classifications. We set the fitness of each individual in the jth level classification in Ϝ
1 as M
+
L
-
j
(
j
=
1,2
,
…
,
L
), the fitness of each individual in Ϝ
k as M
+
1
-
k
(
k
=
2,3
,
…
,
M
). In this manner, the fitness is higher when the violation magnitude of the constraints is smaller.The fitness of feasible individuals (individuals inϜ
1) is determined by the classification according to the downrange and crossrange. We will detail the classification approach of the individuals in Ϝ
1. First, split the whole interval into N subintervals evenly according to the span of the downrange of Ϝ
1, and then fill all individuals in each subinterval. Denote all individuals in the ith subinterval as χ
i (suppose it has κ
i individuals); then L
=
max
⁡
(
κ
i
). The jth level classification includes each subinterval's individual corresponding to the jth most favorite (largest/smallest) crossrange. The individuals corresponding to the minimum downrange and the maximum downrange of Ϝ
1 are in the first classification too. By picking up favorite individuals from each subinterval, this method balances the optimality in crossrange and the diversity in the distribution of downrange.In order to make it easier to understand, a schematic diagram is shown in Figure2 (N
=
3 in this case). The whole interval is divided into three subintervals evenly by the two vertical dashed lines. Fill all individuals in each subinterval; then, χ
1
=
{
①
,
②
,
③
,
④
}, χ
2
=
{
⑤
,
⑥
,
⑦
,
⑧
}, and χ
3
=
{
⑨
,
⑩
}. Therefore, κ
1
=
κ
2
=
4, κ
3
=
2
,
L
=
4. The first level classification includes ②
,
⑤
,
⑨
,
①
,
⑩; the second level classification includes ③
,
⑧; the third level classification includes ④
,
⑥; and the fourth level classification includes ⑦. This scheme for assigning the fitness considers the constraints firstly and then leads the feasible individuals to the boundary of the landing footprint. Through splitting the potential downrange interval into several subintervals evenly and making each subinterval's agent have the same probability to survive, the algorithm balances the optimality in crossrange and the diversity in the distribution of downrange.Figure 2
Schematic diagram for classification of all feasible individuals.The pseudocode of the landing footprint computation method is described as Pseudocode1. First, initialize the population and calculate the corresponding downrange, crossrange, and constraints. Then, evolve the population of DE until the maximum generations are reached. And, the evolving direction is determined by superior agents (i.e., have large fitness) of the population. Run the algorithm to generate one side of the landing footprint according to the maximum crossrange due to different downrange with positive bank angle. Then, use the same algorithm to generate the other side of the landing footprint with negative bank angle.Pseudocode 1:Pseudocode of the landing footprint computation algorithm.
Step  1. Set the population scale NP, the mutant factor F, the crossover probability C
R, the maximum generations G
max
⁡,
the orders of the interpolation polynomials of α and σ (N
1 and N
2), the number of subintervals for downrange
splitting N, and let G
=
0. Initialize the population P
0 using  (19), and calculate the downrange, the crossrange
and the constraints of each individual in P
0 by integrating the state equations  (1)–(6) using the control curves
described as  (18) with decision vector X
=
X
G
,
i.
Step  2. For i
=
1 to NP
Do mutation and crossover with the ith individual X
G
,
i to generate the offspring U
G
,
i. Calculate the
downrange, the crossrange and the constraints of U
G
,
i by integrating the state equations using the control
curves described as  (18) with decision vector X
=
U
G
,
i.
End For
Step  3. Calculate the individuals' fitness of P
G and its offspring.
Step  4. Choose NP most favorite individuals (have larger fitness) from P
G and its offspring to form the next population
P
G
+
1 and set G
=
G
+
1.
Step  5. If G
≤
G
max
⁡
Go toStep  2.
Else
Print the boundary of the landing footprint and the according trajectories, and then terminate the procedure.
End If


## 3.1. Discretization of the Control Profiles

From (8)-(9) and the expression of density ρ in (7), we have
(15)
r
≥
r
0
-
2
h
s
ln
⁡
(
Q
˙
max
⁡
ρ
0
c
v
k
q
)
≜
r
Q
(
v
)
,
r
≥
r
0
-
h
s
ln
⁡
2
q
max
⁡
ρ
0
v
2
≜
r
q
(
v
)
.
Let y
(
α
,
v
,
r
,
σ
)
=
L
cos
⁡
σ
-
g
+
v
2
/
r
=
(
S
/
2
m
)
ρ
v
2
cos
⁡
σ
C
L
(
α
,
v
,
r
)
-
μ
/
r
2
+
v
2
/
r, according to (11), y
≈
0. Solve the equation y
=
0 with dependent variable being α; we have α
=
α
(
v
,
r
,
σ
),
(16)
∂
y
∂
σ
=
-
L
sin
σ
,
∂
y
∂
α
=
ρ
v
2
S
2
m
∂
C
L
∂
α
cos
⁡
σ
,
∂
y
∂
r
=
2
g
r
-
v
2
r
2
-
L
h
s
cos
⁡
σ
+
L
C
L
∂
C
L
∂
r
cos
⁡
σ
,
σ
∈
[
-
9
0
∘
,
9
0
∘
]
,
L
>
0; thus, ∂
y
/
∂
σ
>
0 when σ
∈
[
-
9
0
∘
,
0
∘
) and ∂
y
/
∂
σ
<
0 when σ
∈
(
0
∘
,
9
0
∘
]. Lift coefficient C
L increases as α increases; that is, ∂
y
/
∂
α
>
0. Consider ∂
α
(
v
,
r
,
σ
)
/
∂
σ
=
-
(
∂
y
/
∂
σ
)
/
(
∂
y
/
∂
α
); then, ∂
α
(
v
,
r
,
σ
)
/
∂
σ
<
0 when σ
∈
[
-
9
0
∘
,
0
∘
) and ∂
α
(
v
,
r
,
σ
)
/
∂
σ
>
0 when σ
∈
(
0
∘
,
9
0
∘
]; therefore, α
(
v
,
r
,
σ
)
≥
α
(
v
,
r
,
0
).∂
C
L
/
∂
r and 2
g
/
r are significantly smaller than L
/
h
s, so they can be ignored; then ∂
y
/
∂
r
<
0, ∂
α
(
v
,
r
,
σ
)
/
∂
r
=
-
(
∂
y
/
∂
r
)
/
(
∂
y
/
∂
α
)
>
0; therefore
(17)
α
=
α
(
v
,
r
,
σ
)
≥
α
(
v
,
r
,
0
)
≥
α
(
v
,
max
⁡
(
r
Q
(
v
)
,
r
q
(
v
)
)
,
0
)
≜
α
min
⁡
(
v
)
.
Suppose the maximum angle of attack and bank angle are α
max
⁡ and σ
max
⁡, respectively, and the minimum bank angle is σ
min
⁡. Chebyshev interpolation polynomial is an optimal approximation of the original function in all interpolation polynomials in L
∞ norm-sense [22]. Herein, control profiles are approximated by Chebyshev interpolation polynomials as follows:
(18)
η
(
τ
)
=
∑
i
=
1
N
1
‍
η
i
Φ
i
(
τ
)
,
ξ
(
τ
)
=
∑
i
=
1
N
2
‍
ξ
i
Υ
i
(
τ
)
.
α
(
v
)
=
(
1
-
η
(
τ
)
)
α
min
⁡
(
v
)
+
η
(
τ
)
α
max
⁡
,
σ
(
v
)
=
(
1
-
ξ
(
τ
)
)
σ
min
⁡
+
ξ
(
τ
)
σ
max
⁡
,
where τ
=
2
(
v
-
v
0
)
/
(
v
f
-
v
0
)
-
1
∈
[
-
1,1
], Φ
i
(
τ
) is the ith function of basis functions of N
1th order Chebyshev polynomial interpolation, and Υ
i
(
τ
) is the ith function of basis functions of N
2th order Chebyshev polynomial interpolation. N
1 and N
2 are the orders of the interpolation polynomials of α and σ, respectively.LetX
=
[
η
1
,
η
2
,
…
,
η
N
1
,
ξ
1
,
ξ
2
,
…
,
ξ
N
2
]
,
η
k
,
ξ
k
∈
[
0,1
], and X is the decision vector. When X is known, control profiles are determined by (18); then we can integrate (1)–(6) to figure out the objective function (13), and the constraints (8)–(10) and (12). The nonlinear programming problems are well defined.


## 3.2. Differential Evolution

DE is a simple, efficient, and robust evolutionary algorithm that deals with unconstrained parameter optimization problem [15]. Unconstrained parameter optimization can be formulated to minimize a single D-dimensional function f
(
X
)
,
X
=
[
x
1
,
x
2
,
…
,
x
D
], x
i
∈
[
x
l
,
i
,
x
u
,
i
]. Let P
G
=
{
X
G
,
1
,
X
G
,
2
,
…
,
X
G
,
NP
} represent the Gth generation of the population, where the ith individual's position is X
G
,
i
=
[
x
G
,
i
,
1
,
x
G
,
i
,
2
,
…
,
x
G
,
i
,
D
], i
=
1,2
,
…
,
NP; NP is the population's scale.DE involves two stages: initialization and evolution. Initialization stage uses (19) to generate the initial population P
0 randomly. Then P
0 evolves to P
1, P
1 evolves to P
2
,
…, and so on, until the terminal conditions are reached. Each evolving phase includes three operations, namely, differential mutation, crossover, and selection:
(19)
x
0
,
i
,
j
=
x
l
,
j
+
rand
[
0,1
]
×
(
x
u
,
j
-
x
l
,
j
)
,
i
=
1,2
,
…
,
NP
,
j
=
1,2
,
…
,
D
,
where rand
[
0,1
] is a uniform random number in [
0,1
].Mutation and crossover create offsprings of the population, and the selection operation controls the direction of evolution. The following will describe these operations.LetV
G
,
i
=
[
v
G
,
i
,
1
,
v
G
,
i
,
2
,
…
,
v
G
,
i
,
D
] denote the ith mutant individual. Classical mutation operator is described as follows:
(20)
V
G
,
i
=
X
G
,
r
1
+
F
×
(
X
G
,
r
2
-
X
G
,
r
3
)
,
where r
1
,
r
2, and r
3 are mutually exclusive indices randomly chosen in the range [
1
,
NP
], which are different from the base index i. F
>
0 is the mutant factor.The mutant individualV
G
,
i exchanges its components with the base vector X
G
,
i through crossover operation to get the trial vector U
G
,
i
=
[
u
G
,
i
,
1
,
u
G
,
i
,
2
,
…
,
u
G
,
i
,
D
]. Common crossover operation is described as follows:
(21)
u
G
,
i
,
j
=
{
v
G
,
i
,
j
,
j
=
j
rand
or
rand

[
0,1
]
≤
C
R
x
G
,
i
,
j
,
otherwise
,
where j
rand is a random index in the range [
1
,
D
]. It ensures that the trial vector has at least one component from the mutant vector. C
R
>
0 is the crossover probability.By using one to one greedy selection operation to decide whether the trial vectorU
G
,
i substitutes X
G
,
i to the next iteration. This strategy guarantees that the population evolves to superior region.


## 3.3. Landing Footprint Generation Algorithm

When coming to calculate the landing footprint, the algorithm must guarantee that none of the constraints is violated, and it should have extreme crossrange and maintains diversity in downrange distribution. DE uses fitness to drive the evolving of the population. Here, we design a scheme to determine the fitness of each agent to lead the agent to the boundary of the landing footprint.LetT
G represent the population constituted by P
G and its offspring. At the first stage, sort all individuals in T
G according to the violation magnitude of the constraints (magnitude of the normalized constraints vector) in ascending order. Suppose it has M sets that are denoted as Ϝ
k, k
=
1,2
,
…
,
M. All individuals who do not violate any constraint (these individuals also called feasible individuals) have the highest priority, and they are denoted as Ϝ
1. At the second stage, downrange and crossrange are used to classify Ϝ
1. Suppose it has L classifications. We set the fitness of each individual in the jth level classification in Ϝ
1 as M
+
L
-
j
(
j
=
1,2
,
…
,
L
), the fitness of each individual in Ϝ
k as M
+
1
-
k
(
k
=
2,3
,
…
,
M
). In this manner, the fitness is higher when the violation magnitude of the constraints is smaller.The fitness of feasible individuals (individuals inϜ
1) is determined by the classification according to the downrange and crossrange. We will detail the classification approach of the individuals in Ϝ
1. First, split the whole interval into N subintervals evenly according to the span of the downrange of Ϝ
1, and then fill all individuals in each subinterval. Denote all individuals in the ith subinterval as χ
i (suppose it has κ
i individuals); then L
=
max
⁡
(
κ
i
). The jth level classification includes each subinterval's individual corresponding to the jth most favorite (largest/smallest) crossrange. The individuals corresponding to the minimum downrange and the maximum downrange of Ϝ
1 are in the first classification too. By picking up favorite individuals from each subinterval, this method balances the optimality in crossrange and the diversity in the distribution of downrange.In order to make it easier to understand, a schematic diagram is shown in Figure2 (N
=
3 in this case). The whole interval is divided into three subintervals evenly by the two vertical dashed lines. Fill all individuals in each subinterval; then, χ
1
=
{
①
,
②
,
③
,
④
}, χ
2
=
{
⑤
,
⑥
,
⑦
,
⑧
}, and χ
3
=
{
⑨
,
⑩
}. Therefore, κ
1
=
κ
2
=
4, κ
3
=
2
,
L
=
4. The first level classification includes ②
,
⑤
,
⑨
,
①
,
⑩; the second level classification includes ③
,
⑧; the third level classification includes ④
,
⑥; and the fourth level classification includes ⑦. This scheme for assigning the fitness considers the constraints firstly and then leads the feasible individuals to the boundary of the landing footprint. Through splitting the potential downrange interval into several subintervals evenly and making each subinterval's agent have the same probability to survive, the algorithm balances the optimality in crossrange and the diversity in the distribution of downrange.Figure 2
Schematic diagram for classification of all feasible individuals.The pseudocode of the landing footprint computation method is described as Pseudocode1. First, initialize the population and calculate the corresponding downrange, crossrange, and constraints. Then, evolve the population of DE until the maximum generations are reached. And, the evolving direction is determined by superior agents (i.e., have large fitness) of the population. Run the algorithm to generate one side of the landing footprint according to the maximum crossrange due to different downrange with positive bank angle. Then, use the same algorithm to generate the other side of the landing footprint with negative bank angle.Pseudocode 1:Pseudocode of the landing footprint computation algorithm.
Step  1. Set the population scale NP, the mutant factor F, the crossover probability C
R, the maximum generations G
max
⁡,
the orders of the interpolation polynomials of α and σ (N
1 and N
2), the number of subintervals for downrange
splitting N, and let G
=
0. Initialize the population P
0 using  (19), and calculate the downrange, the crossrange
and the constraints of each individual in P
0 by integrating the state equations  (1)–(6) using the control curves
described as  (18) with decision vector X
=
X
G
,
i.
Step  2. For i
=
1 to NP
Do mutation and crossover with the ith individual X
G
,
i to generate the offspring U
G
,
i. Calculate the
downrange, the crossrange and the constraints of U
G
,
i by integrating the state equations using the control
curves described as  (18) with decision vector X
=
U
G
,
i.
End For
Step  3. Calculate the individuals' fitness of P
G and its offspring.
Step  4. Choose NP most favorite individuals (have larger fitness) from P
G and its offspring to form the next population
P
G
+
1 and set G
=
G
+
1.
Step  5. If G
≤
G
max
⁡
Go toStep  2.
Else
Print the boundary of the landing footprint and the according trajectories, and then terminate the procedure.
End If


## 4. Simulation

An entry scenario of X-33 [23, 24] is used to validate the algorithm described in this paper.The parameters of the entry problem are as follows:(22)
[
r
i
,
θ
i
,
ϕ
i
,
v
i
,
γ
i
,
ψ
i
]
=
[
6499518
m
,
-
117.0
1
∘
,
-
18.25
5
∘
,
c
7622
m
/
s
,
-
1.437
9
∘
,
30
.
0
∘
]
[
r
f
,
v
f
,
γ
f
]
=
[
6408427
m
,
908
m
/
s
,
-
7
.
5
∘
]
,
[
Δ
r
,
Δ
v
,
Δ
γ
,
]
=
[
2000
m
,
20
m
/
s
,
1
.
0
∘
]
,
Q
max
⁡
=
431259
W
/
m
2
,
q
max
⁡
=
11970
N
/
m
2
,
n
max
⁡
=
2.5
g
,
5
∘
≤
α
≤
4
5
∘
,
|
α
˙
|
≤
5
∘
s
,
|
σ
|
≤
9
0
∘
,
|
σ
˙
|
≤
5
∘
s
.
The parameters of the algorithm: NP
=
200, F
=
0.7, C
R
=
0.5, G
max
⁡
=
500, N
=
50, N
1
=
5
,
N
2
=
12.


### 4.1. Landing Footprint of X-33 in Nominal Configuration

The landing footprint of X-33 in nominal configuration is described in Figure3, and the center of each circle in the figure is the terminal site of each trajectory. Control variables are described in Figures 4(a) and 4(b), and the constraints are described in Table 1. It can be found that all constraints are satisfied and the control profiles are smooth. By observing the control profiles, the angle of attack varies a lot to achieve optimality which indicates that it is necessary to take the angle of attack under consideration.Table 1
Constraints of the trajectories according to the landing footprint.
max
⁡
(
Q
˙
)
max
⁡
(
q
)
max
⁡
(
n
)
|
Δ
r
(
t
f
)
|
|
Δ
v
(
t
f
)
|
|
Δ
γ
(
t
f
)
|
(W/m2)
(N/m2)
(g)
(m)
(m/s)
(degree)
Maximum
431249
9916
2.41
1952
19.7
0.96
Mean
406557
8614
1.83
955
9.6
0.44
Minimum
346819
7242
1.56
5
0.8
0.01
Standard
26432
584
0.13
442
5.9
0.28Figure 3
Landing footprint with ground tracks of X-33 in nominal configuration.Control curves according to the boundary of the landing footprint in nominal configuration.
(a)
Angle of attack
(b)
Bank angleThe time consumption to calculate the landing footprint is323.2
s; that is, the mean cost to calculate each landing point is 3.2
s, under the environment:Compiler: Microsoft Visual Studio 2010
CPU: AMD Athlon(tm)Π X3 435
RAM: 4 G.A conventional method is used to verify our method. This method solves (
2
N
+
2
) optimal control problems (optimal objectives are described as (13)). Each run generates one point in the boundary of the landing footprint. This method discretizes the control profiles as (18) and applies the conventional DE algorithm to solve each optimal control problem. The landing footprint calculated by this method is shown in Figure 3 (the red line with +) too. It shows that the two boundaries are almost the same, except the near boundary of the landing footprint. Near boundary corresponds to large path constraints and these trajectories are dangerous to operate, so the conservatism of our method has little effect in engineering. The time consumption of the conventional method is 15852.3
s under the same environment; that is, it takes 155.4
s to calculate one point of the boundary. It demonstrates that our method is very efficient (nearly 50 times fast) compared with the conventional method.


### 4.2.  Influences of Vehicle's Design Parameters on Landing Footprint

By analyzing the dynamic equations of entry vehicles (1)–(6), the vehicle is maneuvered by the lift and drag acceleration (L and D). From (7), L and D are determined by four design parameters: the mass of the vehicle m, the reference area of the vehicle S, the lift coefficient C
L, and the drag coefficient C
D. The deviations of the mass m and the reference area S can be transformed to the deviations of the lift and drag coefficients C
L and C
D equivalently. For example, increasing the reference S is equivalent to increase C
L and C
D simultaneously at the same rate. To analyze the effect of the design parameters on the landing footprint of the entry vehicle, nine different deviations of the lift and drag coefficients described in Table 2 are used to generate the landing footprints. Cases (
1
)
–
(
9
) are in ascending order according to the ratio between the lift and drag coefficients (lift-to-drag ratio). The landing footprint is described by Figure 5. From the figure, we can find out that the landing footprint is larger when the lift coefficient becomes bigger or the drag coefficient becomes smaller, but the landing footprint is almost the same when the deviations of the lift and drag coefficients are in the same manner (i.e., lift-to-drag ratio does not change). It is indicated that the lift-to-drag ratio is the most important factor to the attainability of the entry vehicle and the attainability enhances as the lift-to-drag becomes larger.Table 2
Deviations of the lift and drag coefficients.
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
Δ
C
L
-
5.0
%
-
5.0
%
0.0
%
-
5.0
%
0.0
%
+
5.0
%
+
5.0
%
0.0
%
+
5.0
%
Δ
C
D
+
5.0
%
0.0
%
+
5.0
%
-
5.0
%
0.0
%
+
5.0
%
0.0
%
-
5.0
%
-
5.0
%
(
1
+
Δ
C
L
)
(
1
+
Δ
C
D
)
-
1
-
9.5
%
-
5.0
%
-
4.8
%
0.0
%
0.0
%
0.0
%
5.0
%
5.3
%
10.5
%Figure 5
Landing footprints with ground tracks in various lift and drag configurations of X-33.


## 4.1. Landing Footprint of X-33 in Nominal Configuration

The landing footprint of X-33 in nominal configuration is described in Figure3, and the center of each circle in the figure is the terminal site of each trajectory. Control variables are described in Figures 4(a) and 4(b), and the constraints are described in Table 1. It can be found that all constraints are satisfied and the control profiles are smooth. By observing the control profiles, the angle of attack varies a lot to achieve optimality which indicates that it is necessary to take the angle of attack under consideration.Table 1
Constraints of the trajectories according to the landing footprint.
max
⁡
(
Q
˙
)
max
⁡
(
q
)
max
⁡
(
n
)
|
Δ
r
(
t
f
)
|
|
Δ
v
(
t
f
)
|
|
Δ
γ
(
t
f
)
|
(W/m2)
(N/m2)
(g)
(m)
(m/s)
(degree)
Maximum
431249
9916
2.41
1952
19.7
0.96
Mean
406557
8614
1.83
955
9.6
0.44
Minimum
346819
7242
1.56
5
0.8
0.01
Standard
26432
584
0.13
442
5.9
0.28Figure 3
Landing footprint with ground tracks of X-33 in nominal configuration.Control curves according to the boundary of the landing footprint in nominal configuration.
(a)
Angle of attack
(b)
Bank angleThe time consumption to calculate the landing footprint is323.2
s; that is, the mean cost to calculate each landing point is 3.2
s, under the environment:Compiler: Microsoft Visual Studio 2010
CPU: AMD Athlon(tm)Π X3 435
RAM: 4 G.A conventional method is used to verify our method. This method solves (
2
N
+
2
) optimal control problems (optimal objectives are described as (13)). Each run generates one point in the boundary of the landing footprint. This method discretizes the control profiles as (18) and applies the conventional DE algorithm to solve each optimal control problem. The landing footprint calculated by this method is shown in Figure 3 (the red line with +) too. It shows that the two boundaries are almost the same, except the near boundary of the landing footprint. Near boundary corresponds to large path constraints and these trajectories are dangerous to operate, so the conservatism of our method has little effect in engineering. The time consumption of the conventional method is 15852.3
s under the same environment; that is, it takes 155.4
s to calculate one point of the boundary. It demonstrates that our method is very efficient (nearly 50 times fast) compared with the conventional method.


## 4.2.  Influences of Vehicle's Design Parameters on Landing Footprint

By analyzing the dynamic equations of entry vehicles (1)–(6), the vehicle is maneuvered by the lift and drag acceleration (L and D). From (7), L and D are determined by four design parameters: the mass of the vehicle m, the reference area of the vehicle S, the lift coefficient C
L, and the drag coefficient C
D. The deviations of the mass m and the reference area S can be transformed to the deviations of the lift and drag coefficients C
L and C
D equivalently. For example, increasing the reference S is equivalent to increase C
L and C
D simultaneously at the same rate. To analyze the effect of the design parameters on the landing footprint of the entry vehicle, nine different deviations of the lift and drag coefficients described in Table 2 are used to generate the landing footprints. Cases (
1
)
–
(
9
) are in ascending order according to the ratio between the lift and drag coefficients (lift-to-drag ratio). The landing footprint is described by Figure 5. From the figure, we can find out that the landing footprint is larger when the lift coefficient becomes bigger or the drag coefficient becomes smaller, but the landing footprint is almost the same when the deviations of the lift and drag coefficients are in the same manner (i.e., lift-to-drag ratio does not change). It is indicated that the lift-to-drag ratio is the most important factor to the attainability of the entry vehicle and the attainability enhances as the lift-to-drag becomes larger.Table 2
Deviations of the lift and drag coefficients.
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
Δ
C
L
-
5.0
%
-
5.0
%
0.0
%
-
5.0
%
0.0
%
+
5.0
%
+
5.0
%
0.0
%
+
5.0
%
Δ
C
D
+
5.0
%
0.0
%
+
5.0
%
-
5.0
%
0.0
%
+
5.0
%
0.0
%
-
5.0
%
-
5.0
%
(
1
+
Δ
C
L
)
(
1
+
Δ
C
D
)
-
1
-
9.5
%
-
5.0
%
-
4.8
%
0.0
%
0.0
%
0.0
%
5.0
%
5.3
%
10.5
%Figure 5
Landing footprints with ground tracks in various lift and drag configurations of X-33.


## 5. Conclusion

This paper developed a parallel algorithm to generate the attainable region at the end of the entry phase based on DE and Chebyshev polynomial interpolation. By exploiting the parallel feature of DE, a scheme is designed to lead the individuals of DE's population evolving to the boundary of the landing footprint. It generates one side of the boundary of the attainable region while conventional methods generate one point of the boundary in each run. Furthermore, the algorithm takes all control signals into consideration and enforces path and terminal constraints, which make it have high fidelity. An entry scenario of X-33 is used to validate the algorithm and some analyses are conducted. It is indicated that the algorithm has the capacity to deal with multiconstraints and to maintain the diversity in the distribution of downrange in the landing footprint. Through the analyses, we can find out that the lift-to-drag ratio is the most important factor to the attainability of the entry vehicle and the attainability enhances as the lift-to-drag gets larger.



---

*Source: 306947-2014-03-17.xml*