Abstract

Generation of optimal reentry trajectory for a hypersonic vehicle (HV) satisfying both boundary conditions and path constraints is a challenging task. As a relatively new swarm intelligent algorithm, an adaptive fireworks algorithm (AFWA) has exhibited promising performance on some optimization problems. However, with respect to the optimal reentry trajectory generation under constraints, the AFWA may fall into local optimum, since the individuals including fireworks and sparks are not well informed by the whole swarm. In this paper, we propose an improved AFWA to generate the optimal reentry trajectory under constraints. First, via the Chebyshev polynomial interpolation, the trajectory optimization problem with infinite dimensions is transformed to a nonlinear programming problem (NLP) with finite dimension, and the scope of angle of attack (AOA) is obtained by path constraints to reduce the difficulty of the optimization. To solve the problem, an improved AFWA with a new mutation strategy is developed, where the fireworks can learn from more individuals by the new mutation operator. This strategy significantly enhances the interactions between the fireworks and sparks and thus increases the diversity of population and improves the global search capability. Besides, a constraint-handling technique based on an adaptive penalty function and distance measure is developed to deal with multiple constraints. The numerical simulations of two reentry scenarios for HV demonstrate the validity and effectiveness of the proposed improved AFWA optimization method, when compared with other optimization methods.

1. Introduction

In recent decades, global strike and space transportation have spurred more and more interests in hypersonic vehicles (HVs) [1] for both military and civilian applications. The development of advanced guidance and control technologies for HV is promoted to meet the need for an effective and reliable access to the space. With aerodynamic control, the unpowered HV has the ability of reentering and gliding through the atmosphere. In order to steer an efficient and safety flight in the complex conditions, the reentry trajectory optimization problem of HV has been widely of concern. The aim of reentry trajectory optimization is to find an optimal solution under the reentry flight dynamics and physics, given the aerodynamics, structural strength, and the thermal protection system [2]. Besides, as the reentry dynamics is highly nonlinear, the reentry trajectory optimization is a nonconvex problem with multiple constraints, such as the control ability, heating rate, dynamic pressure, and aerodynamic load. Thus, it is difficult to solve these problems analytically, and numerical techniques are required to determine an approximation to the continuous solution.

There are two main kinds of numerical methods to solve trajectory optimization problems: indirect methods and direct methods [3]. In indirect methods [4], the original problem is transformed into a two-point boundary problem by applying the calculus of variations or the Pontryagin maximum principle. It has high accuracy, and it ensures that the solution satisfies the necessary optimality conditions. However, in indirect methods, it is quite complicated to derive the necessary optimality conditions. Besides, the radius of convergence is small, and it is difficult to guess the initial value of costate variable [5].

In order to overcome the disadvantages of indirect methods, direct methods have been developed, which are classified into two main types: direct shooting method and collocation methods. In the direct shooting method, only the control variables are parameterized and explicit numerical integration is used to satisfy the differential equation constraints. As an improvement, collocation methods [6–9] discretize both the states and control variables, and they use piecewise or global polynomials to approximate the differential equations at collocation points. Pseudospectral methods (PMs) are frequently used collocation methods, which include Gauss PM (GPM) [10, 11], Legendre PM (LPM), Radau PM (RPM) [12], and Chebyshev PM (CPM) [13]. By direct methods, the reentry trajectory optimization problem with infinite dimensions is transformed into the nonlinear programming problem (NLP) with finite dimensions, which is usually a high-dimensional multimodal nonsmooth problem. Some deterministic methods, such as sequential quadratic programming (SQP) [6, 14], second-order cone programming (SOCP) [15], and sequential convex programming method [16], are used to solve it. However, deterministic algorithms are sensitive, and they are easy to be stagnated at a local optimal point.

In the recent years, intelligent algorithms with the global searching ability, which are easy to implement, have been applied to solve the reentry trajectory optimization problems with the development of computational technology [14, 17–24]. In [17], a hybrid genetic algorithm (GA) collocation method was introduced for trajectory optimization, in which the initial guesses for the state and control variables are interpolated by the best solution of GA. Zhao and Zhou employed the particle swarm optimization (PSO) to obtain the end-to-end trajectory for hypersonic reentry vehicles in a quite simple formulation [19]. Zhang et al. applied the ant colony algorithm (ACO) to generate the optimal reentry trajectory for a reusable launch vehicle (RLV) [20]. Additionally, differential evolution (DE) [21, 22] that mimicked the process of nature evolution and simulated annealing (SA) inspired by annealing in metallurgy was also introduced to deal with reentry trajectory optimization problems.

Among these intelligent algorithms, the fireworks algorithm (FWA) is a relatively new swarm intelligence method firstly proposed by Tan and Zhu [25]. Inspired by the fireworks explosion, the algorithm selects some locations as the fireworks in space, each for exploding to generate a shower of sparks. Through the explosion procedure of the FWA, the diversity of population is enhanced. Besides, owing to the powerful local search capabilities and distributed parallel search mechanism, the FWA has a faster convergence speed compared with other intelligent algorithms.

The enhanced fireworks algorithm (EFWA) [26] is an improved version of the FWA in some operators. To improve the EFWA, the adaptive fireworks algorithm (AFWA) [27] was proposed with an adaptive explosion amplitude, which is calculated according to the fitness of individuals adaptively. Afterwards, some improved FWA were proposed and applied in solving various practical optimization problems. For example, Gao and Ming [28] combined the cultural algorithm (CA) with the FWA for the digital filter design. Zheng et al. [29] developed a hybrid FWA with DE operators to improve the diversity and avoid prematurity. Based on a self-adaption principle and bimodal Gaussian function, Xue et al. [30] proposed an advanced FWA to design the PID controllers with high-quality performances. From the previous research, it is recognized that the information interaction among all individuals of the FWA is relatively poor, whereas the interaction is vital in swarm intelligence algorithms. Thus, when solving complex multimodal problem of the optimal reentry trajectory generation, it may be easy to get trapped in a local optimum. As far as we are concerned, the application of the FWA for reentry trajectory optimization is scarce and has not been found in the published literature.

In this paper, we focus on the improvement of the AFWA and its application to the reentry trajectory optimization problems. Firstly, an improved version of the AFWA (I-AFWA) is developed by combining the standard AFWA (S-AFWA) with a new mutation strategy. In each iteration of the algorithm, the new mutation operator and auxiliary mutation individual selection strategy are applied to make Gaussian fireworks learn from more individuals (not only the best individual) and generate diverse sparks from all fireworks and explosion sparks. The interactions of the fireworks and sparks are enhanced to further improve the diversity of the population and avoid being trapped in local optima too early. Then, the I-AFWA is applied to the reentry trajectory optimization problems. The problem is transformed into NLP by using Chebyshev polynomial interpolation to discretize the angle of attack (AOA) and bank angle simultaneously, and the scope of the AOA is figured out by path constraints to simplify the problem. Next, a constraint-handling technique based on the adaptive penalty function and the distance measure is proposed to incorporate with the I-AFWA and used to deal with the multiconstrained parameter optimization problem. Finally, two reentry scenarios are given to illustrate the validity and effectiveness of the proposed I-AFWA method on the reentry trajectory optimization problem.

The remainder of this paper is organized as follows. The general reentry trajectory optimization problem is formulated in Section 2. Section 3 proposes the I-AFWA with a new mutation strategy, and a constraint-handling technique is incorporated with the I-AFWA for the reentry trajectory optimization. In Section 4, two different reentry tasks for HV are presented to evaluate the proposed algorithm. A few conclusions are made in Section 5.

2. Problem Formulation

2.1. Reentry Dynamics

Using a spherical rotating Earth model, the three-degree-of-freedom (3DOF) point mass dynamics equations [31] of a hypersonic reentry vehicle are described as follows: where is the radial distance from the center of the Earth to the vehicle. and are the longitude and latitude, respectively. is the Earth-relative velocity. is the velocity flight path angle and is the velocity heading angle measured from the North in a clockwise direction. is the bank angle, which is the angle between the vehicle longitudinal symmetry plane and the vertical plane [2]. is the Earth’s self-rotation angular velocity, and is the gravitational acceleration where is the gravitational constant. and are the aerodynamic lift force and drag force, respectively, which are given as follows: where is the density of atmosphere, with as the density of atmosphere at sea level, as the scalar height coefficient, and as the altitude from the sea level. , where denotes the Earth’s radius. and are the aerodynamic reference area and mass of the vehicle, respectively. and are the lift and drag coefficients determined by the AOA and Mach number . denotes the angle between the relative velocity and a reference line fixed with respect to the vehicle body [2], which does not occur explicitly in the equations and appear through the aerodynamic forces and . Given the altitude , velocity , and the AOA , the Mach number is calculated, and then the corresponding lift and drag coefficient and values are found by looking up the table functions with the value of and . The independent variable is the time , and control variables are and . The geometric sketch of reentry flight is shown in Figure 1.

Figure 1 

Geometric sketch of reentry flight.

2.2. Multiple Constraints

2.2.1. Heating Rate

indicates the heating rate at a specified point on the surface of the reentry vehicle, which normally chooses the stagnation point on the nose. It is determined by the atmospheric density and Earth-relative velocity [31], which must be limited as follows: where is the heating rate normalization constant. is the constant coefficient that typically takes the value of 3.0, 3.15, or 3.25. represents the maximum allowable value of heating rate.

2.2.2. Dynamic Pressure

The dynamic pressure is a typical path constraint, and it is limited to control the hinge moment of an actuator in a reasonable range. The peak value of the dynamic pressure must be less than the maximum value as follows:

2.2.3. Aerodynamic Load

In order to protect the mechanical system, the instantaneous value of aerodynamic load in the body-normal direction is given lower than the maximum value as follows:

2.2.4. Quasi-Equilibrium Glide Condition

As the flight path angle is small and varies relatively slowly in a major portion of the lifting reentry trajectory [32], we set and ignore the Earth’s self-rotation in (5), and it generates

This is the so called quasi-equilibrium glide condition (QEGC) [32], which guarantees that the reentry trajectory does not oscillate acutely for the convenience of the attitude controller design. It is a soft condition that forms the upper boundary of altitude-velocity reentry corridor. Additionally, hard constraints including heating rate, dynamic pressure, and aerodynamic load form the lower boundary of altitude reentry corridor. Therefore, the altitude-velocity reentry corridor is illustrated in Figure 2.

Figure 2 

The sketch map of altitude-velocity reentry corridor.

2.2.5. Control Boundary

To maintain the stability of the reentry vehicle, it is necessary to limit the two control variables of AOA and bank angle as follows:

2.2.6. Terminal Conditions

Given the different reentry flight missions and conditions, the terminal states, including terminal altitude, position, velocity, flight path angle, and heading angle, are restrained as follows: where is the terminal time of the reentry phase, and , and represent the ideal reentry terminal states, respectively.

2.3. Reentry Trajectory Optimization Problem Formulation

The objective of the reentry trajectory optimization problem is to design a reentry trajectory that can minimize (or maximize) the performance index, under multiple constraints. Hence, it is a multiconstrained nonlinear optimal control problem, and here, we consider it in the Bolza form [33].

The continuous Bolza problem is to determine the state and the control that minimize the Bolza objective function as follows: subject to the dynamics equations, boundary constraints, and path constraints, which are stated formally as follows: where is the endpoint performance index and is the integral performance index. In the reentry trajectory optimization problem, is the state vector of the vehicle and is the control vector. and are the initial and terminal time, respectively, of the reentry flight.

Generally, the objective function is defined corresponding to the specified mission. For the reentry vehicle, the objective function can be selected to minimize downrange, to maximize crossrange, or to minimize total heat and so on.

The reentry trajectory optimization is formulated as an optimal control problem by (14), (15), (16), and (17), in which a dynamics model is represented by (1), (2), (3), (4), (5), and (6), multiple constraints are given in (8), (9), (10), (11), (12), and (13), and the performance index function is described by (14). Thus, an improved AFWA is developed to effectively solve the nonconvex optimal control problem under constraints.

3. Trajectory Optimization Algorithm Description

Reentry trajectory optimization problem is a continuous optimal control problem with infinite dimensions, so it is difficult to solve directly. To solve this problem, we first transform it to a finite dimensional NLP by the discretization of control curves with Chebyshev polynomial interpolation. However, it is also difficult to optimize the AOA and bank angle simultaneously, so we obtain the scope of AOA by path constraints ((8), (9), (10), and (11)). Subsequently, a constraint-handling technique based on an adaptive penalty function and a distance measure is proposed to deal with multiple constraints. Finally, a new improved AFWA with a new mutation strategy is developed to solve the NLP.

3.1. Discretization of the Control Curves

Substitute the expression of density into path constraints of (8) and (9), and we obtain the following equations:

Let . Based on the QEGC in (11), . Then, by substituting (7) into it, we obtain . Moreover, the derivation of , , and can be calculated as follows:

Since and , when and when . Lift coefficient increases as increases; that is, ; thus, . Consider ; then, when and when . Furthermore, is much larger than and , so they can be ignored. As and , we have

Denote , , and as the maximum allowable AOA, maximum allowable bank angle, and minimum allowable bank angle, respectively. Then, the scopes of control variables are determined. With the scopes, the discretization process is conducted in the following.

In the approximation of the original function, polynomial interpolation at the zeros of Chebyshev orthogonal polynomial can minimize the maximum error in interpolation intervals. Thus, Chebyshev polynomial interpolation is the optimal approximation of the original function in all polynomial interpolations in L_∞ norm sense. Hence, the control curves are approximated by Chebyshev polynomial interpolation as follows: where , is the interpolation basis function of -order Chebyshev interpolation polynomial, and is the interpolation basis function of -order Chebyshev interpolation polynomial. and are the orders of the interpolation polynomial and , respectively.

Set , in which is the decision vector of the transformed NLP. If is obtained, control curves are determined by (21), and then we can integrate (1), (2), (3), (4), (5), and (6) to figure out the state curves, objective function, and constraints. In this way, the NLP is solved.

3.2. Adaptive Fireworks Algorithm

In order to solve the NLP and obtain the optimal reentry trajectory, it is necessary to develop an effective intelligent optimization algorithm. The FWA [25] is a novel swarm intelligence algorithm inspired by fireworks explosion. Later, the AFWA is proposed with an adaptive explosion amplitude to improve the local search capability. It is efficient for unconstrained parameter optimization problem formulated to minimize D-dimensional function as follows: where and are the lower and upper boundaries, respectively, of the search space in dimension . Besides, , represents the position of the firework, where is the size of fireworks. The process of the AFWA is as follows.

3.2.1. Initialization

In the initialization stage, locations are randomly generated in the search space as the fireworks of the first generation as follows: where is the uniform distribution over interval .

3.2.2. Calculation of Spark Number and Explosion Amplitude

For each firework, a certain number of sparks are generated within a certain explosion amplitude. The number of sparks and explosion amplitude are calculated according to their fitness. The better the firework fitness is, the more sparks it will generate and the smaller its amplitude is, vice versa.

The number of explosion sparks for each firework is calculated as follows: where is the total number of explosion sparks generated by fireworks, stands for the firework with the worst (maximum) fitness among fireworks, refers to the fitness value of the firework, and denotes the smallest constant to avoid zero-division-error.

In order to avoid the overwhelming effects of splendid fireworks, the number of sparks is bounded as follows: where and are the lower bound and upper bound, respectively, for the spark numbers.

The explosion amplitudes of the fireworks (except for the optimal firework, i.e., the firework with the best fitness) are calculated as follows: where is the maximum explosion amplitude and stands for the firework with the best (minimum) fitness among the fireworks.

The initial explosion amplitude of the optimal firework is .

3.2.3. Generation of Explosion Sparks

The explosion sparks are generated randomly within the hypersphere with the firework as its center and the amplitude as its radius. In explosion, sparks may undergo the effects of explosion from random directions (dimensions). The affected dimensions are randomly selected by . For the selected dimensions (those dimensions, where equals to 1) of the firework, the different offset displacements are calculated by . The location of each spark generated by can be obtained by adding the displacement to , that is, where denotes the set of the preselected dimensions of the firework.

3.2.4. Generation of Gaussian Mutation Sparks

To keep the diversity of sparks and avoid falling into the local optimum, the Gaussian mutation process is introduced. Each selected firework stretches along the direction between the current location of the firework and the location of the best firework. For the selected dimensions of the randomly selected firework , the location of the Gaussian mutation spark is generated as follows: where the offset displacement is the Gaussian distributed with mean 0 and standard deviation 1, and it is utilized to define the coefficient of the mutation. is the number of Gaussian mutation sparks and is the dimensional position of the best firework in the current generation.

3.2.5. Mapping Rule

When the location of a new generated spark (explosion spark or mutation spark) exceeds the search range in dimension , this spark will be mapped to another location within the search space with uniform distribution according to

3.2.6. Calculation of the Adaptive Amplitude of the Optimal Firework

The optimal firework in the next generation is the best individual found (could be a spark generated or a firework) in this generation. To calculate the adaptive amplitude, we choose an individual and use its distance to the best individual as the amplitude of the optimal firework in the next explosion. The chosen individual should subject to the following conditions: (1)Its fitness is worse than that of the optimal firework in this generation.(2)Its distance to the best individual is minimal among all individuals satisfying 1.

Namely, with constraint , where denotes the best individual among all sparks and fireworks in generation , namely, the optimal firework of generation , and is the optimal firework of generation . here represents all individuals in this generation. is a certain measurement of distance.

Then, the adaptive amplitude of the optimal firework in generation is calculated as follows: where and represent the adaptive amplitude in generation and , respectively. denotes the infinity norm, which is used as the distance measure, namely, the maximum difference among all dimensions. is a certain coefficient (usually bigger than 1) that is used to further slow down the convergence rate and improve the global search.

3.2.7. Selection

Elitism-random selection method [34] is used for the selection. The best individual will be selected first as the optimal firework in the next generation. Then, the other individuals are selected randomly.

3.3. Constraint Handling

In this subsection, we consider how to deal with the various constraints of the optimal reentry trajectory generation. Typically, a penalty function is added to the objective function to solve the constrained parameter optimization problem. However, it is difficult to choose appropriate weight coefficient for each constraint. In this case, a distance measure method and an adaptive penalty function are incorporated to form a fitness function to deal with multiple constraints and check dominance in the population. The values of the two functions vary with the sum of constraint violations and the objective function value of an individual.

3.3.1. Distance Measure

The distance measure is obtained by introducing the effect of an individual’s constraint violation into its objective function. Firstly, the equality constraints defined in (16) are transformed into the following inequality constraints: where is the small constant that represents the tolerance value for equality constraints.

Since the magnitude of an individual’s constraint value and its objective function are different, it is necessary to normalize the two indexes. With the maximum and minimum objective function values of the current generation, the normalized objective function is described as follows:

The constraint violation of an individual is then calculated as the summation of the normalized violations of each constraint divided by the total number of constraints. with the constraint violation that is defined as follows: where is the maximum value of the constraint violation. is the total number of constraints, is the number of inequality constraints, and is the number of equality constraints. If the individual satisfies the constraint, then the constraint violation is set to zero; otherwise, will be the actual constraint value, which is greater than zero. Particularly, if is a feasible individual, namely, all the constraints of are satisfied, then , and , which make the constraint violation in (34) run into the zero-division-error. Hence, we define in this case.

Thus, the distance measure for each individual is formulated with the objective value and constraint violation as follows: where is the proportion of feasible solutions in the current population, and it is described as

3.3.2. Adaptive Penalty Function

Additionally, an adaptive penalty function is added to the fitness value of infeasible individuals, which contains two parts: one is based on the constraint violation and the other is based on the objective function. The weight between the two components is controlled by the number of feasible individuals in the current population. The adaptive penalty function is formulated as follows: with and

The weight in the adaptive penalty function allows the tradeoff between finding more feasible solutions and finding better solutions during the evolutionary process.

Thus, the final fitness function of each individual is formulated as the sum of the distance measure and penalty function.

The fitness function formulation of the optimal reentry trajectory generation is flexible, and it helps to solve the multiconstrained optimization problem more efficiently and effectively.

3.4. Improved Adaptive Fireworks Algorithm

Although there are some variants for the AFWA, premature convergence when solving complex multimodal problems of the optimal reentry trajectory generation is still the main deficiency of the AFWA. The search process of the AFWA indicates that the diversification mechanism is not very flexible, and in particular it does not utilize more information about other better quality solutions in the swarm. Especially, the interaction of fireworks in the explosion operator is not sufficient when obtaining the number and amplitude of the fireworks, which can be easily polarized. In other words, the best firework may generate an overwhelming number of sparks, while the other fireworks produce very few sparks. In the mutation operator, the interaction between the explosion sparks is ignored, and thus the diversity of sparks is reduced. Besides, the mutation sparks can hardly be passed down to the next generation [35]. For a D-dimensional optimization problem, the fitness of an individual may be determined by values in all dimensions, and an individual that has discovered the region corresponding to the global optimum in some dimensions may have a low fitness value because of the poor solutions in the other dimensions [36]. Thus, it is vital to generate more diverse sparks in all dimensions.

In order to make use of the beneficial information in the population more effectively to generate better quality solutions, a new mutation strategy is proposed to improve the S-AFWA, which can enhance the diversity of the population and improve the global search capability. In the new mutation strategy, the location of the mutation spark is generated by the selected firework in all dimensions as follows: where defines which individual the selected firework should follow; that is, stretches along the direction between the selected individual and itself. is the corresponding dimension of any individual (all fireworks and explosion sparks, including itself), and the decision depends on the mutation probability , which can take different values for different fireworks.

For each dimension of the firework , a random number is generated. If the random number is greater than , the corresponding dimension will learn from another individual; otherwise, it will learn from itself. We employ a selection procedure of when the dimension of the firework learns from another individual as follows.

Firstly, two individuals (excluding the selected firework ) are randomly selected among the population. Subsequently, compare the fitness values of the two individuals and select the better one. At last, the winner is used as the exemplar to learn from for that dimension.

After the above selection procedure, if all exemplars in different dimensions of the firework are itself (i.e., in all dimensions; due to that, all generated random numbers above are less than ), we will randomly choose one dimension to learn from the corresponding dimension of another individual. The detailed procedure of choosing in (42) is given in Figure 3, where the ceil operator obtains the smallest integer greater than or equal to the specified expression and is the population size.

Figure 3 

Selection procedure of auxiliary mutation individuals for the firework in all dimensions.

The new mutation strategy can generate new spark in the search space using the information from all the individuals, and thus the interaction between the fireworks and sparks is enhanced, not limited to the selected firework and the best firework. The proposed mutation strategy exhibits a larger potential search space for mutation sparks than that of the S-AFWA, and thus the diversity is also enhanced. Especially when the best firework falls into the local optimum region, the mutation sparks have the increasing ability to jump out of the local optimum via the cooperative behavior of the whole population.

In order to address optimization problems in a general manner, each firework has a different mutation probability . Therefore, fireworks have different levels of exploration and exploitation ability in the population, and they are able to solve diverse problems. We empirically set the mutation probability for each firework as follows:

In detail, the procedure of reentry trajectory optimization based on the I-AFWA with a new mutation strategy is depicted in Figure 4.

Figure 4 

Reentry trajectory optimization based on the I-AFWA with a new mutation strategy.

4. Simulation and Results

In this section, one HV model [32] is used to validate the proposed algorithm. The reentry scenario is set as follows: where and is the Earth’s radius. Hence, and .

In order to verify the effectiveness of the proposed algorithm, the I-AFWA is tested in comparison with the S-AFWA, DE, PSO, GA, and FWA. These algorithms all run 20 times repeatedly based on a similar simulation environment. The population size (spark size in the I-AFWA, S-AFWA, and FWA and particle number in DE, PSO, and GA) is set as 50, and the maximum number of function evaluations for all algorithms. The AOA discrete point number , and bank angle discrete point number . Besides, other main parameters are shown in Table 1.

The hard path constraint limits and control boundaries remain fixed throughout all simulations: