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.
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 accelerationbased trajectory planner. Meanwhile, a nearoptimal 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 nearoptimal closedform control law by using quasiequilibrium glide condition (QEGC) [6]. It transformed the original problem to several closestapproach problems which were univariate rootfinding 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 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.
2. 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, is the earthrelative velocity, is the relative flightpath angle, is the relative velocity heading angle measured clockwise from the north, is the earth selfrotation rate, is the vehicle's mass, is the reference area, is the density of atmosphere, is the density of atmosphere at sea level, 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 , , and 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 , , 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 normplan) and the optimal objective is maximizing or minimizing terminal crossrange (the range perpendicular to the normplan).
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 () optimal control problems whose objectives are described as follows: where is the number of the calculated sites to approximate the lateral boundary: where , and 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 Let , according to (11), . Solve the equation with dependent variable being ; we have , ; thus, when and when . Lift coefficient increases as increases; that is, . Consider ; then, when and when ; therefore, .
and are significantly smaller than , so they can be ignored; then , ; therefore Suppose the maximum angle of attack and bank angle are and , respectively, and the minimum bank angle is . Chebyshev interpolation polynomial is an optimal approximation of the original function in all interpolation polynomials in normsense [22]. Herein, control profiles are approximated by Chebyshev interpolation polynomials as follows: where , is the th function of basis functions of th order Chebyshev polynomial interpolation, and is the th function of basis functions of th order Chebyshev polynomial interpolation. and are the orders of the interpolation polynomials of and , respectively.
Let , and is the decision vector. When 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 dimensional function , . Let represent the th generation of the population, where the th individual's position is , ; is the population's scale.
DE involves two stages: initialization and evolution. Initialization stage uses (19) to generate the initial population randomly. Then evolves to , evolves to , and so on, until the terminal conditions are reached. Each evolving phase includes three operations, namely, differential mutation, crossover, and selection: where is a uniform random number in .
Mutation and crossover create offsprings of the population, and the selection operation controls the direction of evolution. The following will describe these operations.
Let denote the th mutant individual. Classical mutation operator is described as follows: where , and are mutually exclusive indices randomly chosen in the range , which are different from the base index . is the mutant factor.
The mutant individual exchanges its components with the base vector through crossover operation to get the trial vector . Common crossover operation is described as follows: where is a random index in the range . It ensures that the trial vector has at least one component from the mutant vector. is the crossover probability.
By using one to one greedy selection operation to decide whether the trial vector substitutes 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.
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 , . All individuals who do not violate any constraint (these individuals also called feasible individuals) have the highest priority, and they are denoted as . At the second stage, downrange and crossrange are used to classify . Suppose it has classifications. We set the fitness of each individual in the th level classification in as , the fitness of each individual in as . In this manner, the fitness is higher when the violation magnitude of the constraints is smaller.
The fitness of feasible individuals (individuals in ) is determined by the classification according to the downrange and crossrange. We will detail the classification approach of the individuals in . First, split the whole interval into subintervals evenly according to the span of the downrange of , and then fill all individuals in each subinterval. Denote all individuals in the th subinterval as (suppose it has individuals); then . 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 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 ( 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, , , and . Therefore, , . 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.
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.

4. Simulation
An entry scenario of X33 [23, 24] is used to validate the algorithm described in this paper.
The parameters of the entry problem are as follows: The parameters of the algorithm: , , , , , .
4.1. Landing Footprint of X33 in Nominal Configuration
The landing footprint of X33 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.
(a) Angle of attack
(b) Bank angle
The time consumption to calculate the landing footprint is ; that is, the mean cost to calculate each landing point is , 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 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 under the same environment; that is, it takes to calculate one point of the boundary. It demonstrates that our method is very efficient (nearly 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 ( 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 are in ascending order according to the ratio between the lift and drag coefficients (lifttodrag 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., lifttodrag ratio does not change). It is indicated that the lifttodrag ratio is the most important factor to the attainability of the entry vehicle and the attainability enhances as the lifttodrag becomes larger.
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 X33 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 lifttodrag ratio is the most important factor to the attainability of the entry vehicle and the attainability enhances as the lifttodrag gets larger.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is partially supported by National Nature Science Foundation of China no. 61203081 and no. 61174079, Doctoral Fund of Ministry of Education of China no. 20120142120091, and Precision Manufacturing Technology and Equipment for Metal Parts no. 2012DFG70640.