Attainability Analysis for Entry Vehicles Based on Differential Evolution

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.


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][2][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][9][10][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][14][15][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 Section 2. 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.

Problem Formulation
Three degrees of freedom dynamic equations of an entry vehicle are as follows: where  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,  is the vehicle's mass,  is the reference area,  is the density of atmosphere,  0 is the density of atmosphere at sea level,  0 is the earth's radius,   and   are the lift and drag coefficients determined by angle of attack  and Mach ,  is the gravity acceleration, independent variable is the time , and control variables are angle of attack  and bank angle .The geometric sketch of entry flight is illustrated by Figure 1.Entry process must satisfy constraints of heating rate, dynamic pressure, and normal acceleration: where Q max ,  max , and  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: Following the entry phase is Terminal Area Energy Management (TAEM) phase.Terminal state must be restrained as follows: where   is the terminal time of the entry phase and   , V  , and   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).
Let (  ), (  ) 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+2) optimal control problems whose objectives are described as follows: where  is the number of the calculated sites to approximate the lateral boundary: where   ,   ,   , V  ,   , and   are the states at the initial time of the entry phase.

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.
DE involves two stages: initialization and evolution.Initialization stage uses (19) to generate the initial population  0 randomly.Then  0 evolves to  1 ,  1 evolves to  2 , . .., and so on, until the terminal conditions are reached.Each evolving phase includes three operations, namely, differential mutation, crossover, and selection: 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.
Let V , = [V ,,1 , V ,,2 , . . ., V ,, ] denote the th mutant individual.Classical mutation operator is described as follows: where  1 ,  2 , and  3 are mutually exclusive indices randomly chosen in the range [1, NP], which are different from the base index . > 0 is the mutant factor.The mutant individual V , exchanges its components with the base vector X , through crossover operation to get the trial vector U , = [ ,,1 ,  ,,2 , . . .,  ,, ].Common crossover operation is described as follows: where  rand is a random index in the range [1, 𝐷].It ensures that the trial vector has at least one component from the mutant vector.  > 0 is the crossover probability.By using one to one greedy selection operation to decide whether the trial vector U , substitutes X , to the next iteration.This strategy guarantees that the population evolves to superior region.
Step 1. Set the population scale NP, the mutant factor , the crossover probability   , the maximum generations  max , the orders of the interpolation polynomials of  and  ( 1 and  2 ), the number of subintervals for downrange splitting , and let  = 0. Initialize the population  0 using (19), and calculate the downrange, the crossrange and the constraints of each individual in  0 by integrating the state equations ( 1)-( 6) using the control curves described as (18) with decision vector X = X , .
Step 2. For  = 1 to NP Do mutation and crossover with the th individual X , to generate the offspring U , .Calculate the downrange, the crossrange and the constraints of U , by integrating the state equations using the control curves described as (18) with decision vector X = U , .End For Step 3. Calculate the individuals' fitness of   and its offspring.

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.
Let   represent the population constituted by   and its offspring.At the first stage, sort all individuals in   according to the violation magnitude of the constraints (magnitude of the normalized constraints vector) in ascending order.Suppose it has  sets that are denoted as ϝ  ,  = 1, 2, . . ., .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  classifications.We set the fitness of each individual in the th level classification in ϝ 1 as  +  −  ( = 1, 2, . . ., ), the fitness of each individual in ϝ  as  + 1 −  ( = 2, 3, . . ., ).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  subintervals evenly according to the span of the downrange of ϝ 1 , and then fill all individuals in each subinterval.Denote all individuals in the th subinterval as   (suppose it has   individuals); then  = max(  ).The th level classification includes each subinterval's individual corresponding to the th 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 Figure 2 ( = 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 = {A, B, C, D},  2 = {E, F, G, H}, and  3 = {0, 1}.Therefore,  1 =  2 = 4,  3 = 2, = 4.The first level classification includes B, E, 0, A, 1; the second level classification includes C, H; the third level classification includes D, F; and the fourth level classification includes G.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.
The pseudocode of the landing footprint computation method is described as Pseudocode 1. 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.

Simulation
An entry scenario of X-33 [23,24] is used to validate the algorithm described in this paper.

Landing Footprint of X-33 in Nominal Configuration.
The landing footprint of X-33 in nominal configuration is described in Figure 3, 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.
The time consumption to calculate the landing footprint is 323.2 s; that is, the mean cost to calculate each landing point is 3.A conventional method is used to verify our method.This method solves (2 + 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.3s 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.

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 ( and ).From (7),  and  are determined by four design parameters: the mass of the vehicle , the reference area of the vehicle , the lift coefficient   , and the drag coefficient   .The deviations of the mass  and the reference area  can be transformed to the deviations of the lift and drag coefficients   and   equivalently.For     example, increasing the reference  is equivalent to increase   and   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.

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 Mathematical Problems in Engineering 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.

Step 4 .Figure 2 :
Figure 2: Schematic diagram for classification of all feasible individuals.

Figure 3 :
Figure 3: Landing footprint with ground tracks of X-33 in nominal configuration.

Figure 4 :
Figure 4: Control curves according to the boundary of the landing footprint in nominal configuration.

Figure 5 :
Figure 5: Landing footprints with ground tracks in various lift and drag configurations of X-33.
cos   cos   sin   − sin   cos   sin   − cos   sin   cos   cos   , sin   cos   sin   + cos   cos   sin   − sin   sin   cos   cos   , sin   sin   + cos   cos   cos   ] , is the th function of basis functions of  1 th order Chebyshev polynomial interpolation, and Υ  () is the th function of basis functions of  2 th order Chebyshev polynomial interpolation. 1 and  2 are the orders of the interpolation polynomials of  and , respectively.Let X = [ 1 ,  2 , . . .,   1 ,  1 ,  2 , . . .,   2 ]

Table 2 :
Deviations of the lift and drag coefficients.