Reentry Trajectory Optimization for a Hypersonic Vehicle Based on an Improved Adaptive Fireworks Algorithm

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.


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][7][8][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], secondorder 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][18][19][20][21][22][23][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. International Journal of Aerospace Engineering dynamics equations [31] of a hypersonic reentry vehicle are described as follows:

Problem Formulation
where r is the radial distance from the center of the Earth to the vehicle.θ and ϕ are the longitude and latitude, respectively.v 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].ω e is the Earth's self-rotation angular velocity, and g = μ/r 2 is the gravitational acceleration where μ is the gravitational constant.L and D are the aerodynamic lift force and drag force, respectively, which are given as follows: where ρ = ρ 0 exp −h/h s is the density of atmosphere, with ρ 0 as the density of atmosphere at sea level, h s as the scalar height coefficient, and h as the altitude from the sea level.h = r − r e , where r e denotes the Earth's radius.S ref and m are the aerodynamic reference area and mass of the vehicle, respectively.C L and C D are the lift and drag coefficients determined by the AOA α and Mach number M a = v/V s r .α 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 L and D.  [31], which must be limited as follows: where k Q is the heating rate normalization constant.c is the constant coefficient that typically takes the value of 3.0, 3.15, or 3.25.Q s max represents the maximum allowable value of heating rate.
2.2.2.Dynamic Pressure.The dynamic pressure q 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 3 International Journal of Aerospace Engineering dynamic pressure must be less than the maximum value q max 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 n T is given lower than the maximum value n T max 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 γ ≈ 0, γ ≈ 0 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.

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 t f is the terminal time of the reentry phase, and r f , θ f , ϕ f , v f , γ f , and ψ f represent the ideal reentry terminal states, respectively.

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 x t ∈ R n and the control u t ∈ R m 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 g is the integral performance index.In the reentry trajectory optimization problem, x = r, θ, ϕ, v, γ, ψ T is the state vector of the vehicle and u = α, σ T is the control vector.t 0 and t f 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.

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: Based on the QEGC in (11), f r, v, α, σ = 0.Then, by substituting (7) into it, we obtain α = α v, r, σ .Moreover, the derivation of r, α, and σ can be calculated as follows: Furthermore, L/h s is much larger than ∂C L /∂r and 2g/r, so they can be ignored.As ∂f /∂r ≈ −v 2 /r 2 − L cos σ/h s < 0 and ∂α/∂r = − ∂f /∂r / ∂f /∂α > 0, we have Denote α max , σ max , and σ min 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: k τ is the kth interpolation basis function of N α -order Chebyshev interpolation polynomial, and l N σ k τ is the kth interpolation basis function of N σ -order Chebyshev interpolation polynomial.N α and N σ are the orders of the interpolation polynomial α and σ, respectively. Set in which X is the decision vector of the transformed NLP.If X 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.

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 x k,min and x k,max are the lower and upper boundaries, respectively, of the search space in dimension k.Besides, , n, represents the position of the ith firework, where n is the size of fireworks.The process of the AFWA is as follows.
3.2.1.Initialization.In the initialization stage, n locations are randomly generated in the search space as the fireworks of the first generation as follows:
5 International Journal of Aerospace Engineering

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 m is the total number of explosion sparks generated by n fireworks, X # stands for the firework with the worst (maximum) fitness among n fireworks, F X i refers to the fitness value of the ith 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 S min and S max 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 A is the maximum explosion amplitude and X * stands for the firework with the best (minimum) fitness among the n fireworks.
The initial explosion amplitude of the optimal firework is x k,max − x k,min , k = 1, 2, … , D.

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 z k = round rand 0, 1 , k = 1, 2, … , D. For the selected dimensions (those dimensions, where z k equals to 1) of the ith firework, the different offset displacements are calculated by Δx i k = A i ⋅ rand −1, 1 .The location of each spark X j generated by X i can be obtained by adding the displacement Δx i k to x i k , that is, where Z denotes the set of the preselected dimensions of the ith firework.

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 X i , the location of the Gaussian mutation spark X j is generated as follows: where the offset displacement e = Gaussian 0, 1 is the Gaussian distributed with mean 0 and standard deviation 1, and it is utilized to define the coefficient of the mutation.m is the number of Gaussian mutation sparks and x * k is the kth dimensional position of the best firework in the current generation.
3.2.5.Mapping Rule.When the location of a new generated spark X j (explosion spark or mutation spark) exceeds the search range in dimension k, 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, where X * G+1 denotes the best individual among all sparks and fireworks in generation G, namely, the optimal firework of generation G + 1, and X * G is the optimal firework of generation G. X i here represents all individuals in this generation.d is a certain measurement of distance.
Then, the adaptive amplitude of the optimal firework in generation G + 1 is calculated as follows: where A * G and A * G+1 represent the adaptive amplitude in generation G and G + 1, 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 6 International Journal of Aerospace Engineering 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.

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 v X 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 c i,max = max X c i X is the maximum value of the ith constraint violation.n is the total number of constraints, l is the number of inequality constraints, and n − l is the number of equality constraints.If the individual X satisfies the ith constraint, then the constraint violation c i X is set to zero; otherwise, c i X will be the actual constraint value, which is greater than zero.Particularly, if X is a feasible individual, namely, all the constraints of X are satisfied, then c i X = 0, i = 1, … , n, and c i,max = 0, which make the constraint violation in (34) run into the zero-divisionerror.Hence, we define v X = 0 in this case.Thus, the distance measure for each individual is formulated with the objective value and constraint violation as follows: where r f is the proportion of feasible solutions in the current population, and it is described as r f = the number of f easible individuals in the current population the size of the current population , 37 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: and The weight r f 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 X 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 International Journal of Aerospace Engineering 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 X j is generated by the selected firework X i in all dimensions as follows: where a i = a i 1 , a i 2 , … , a i D defines which individual the selected firework X i should follow; that is, X i stretches along the direction between the selected individual and itself.
is the corresponding kth dimension of any individual (all fireworks and explosion sparks, including itself), and the decision depends on the mutation probability P i , which can take different values for different fireworks.
For each dimension of the firework X i , a random number is generated.If the random number is greater than P i , the corresponding dimension will learn from another individual; otherwise, it will learn from itself.We employ a selection procedure of a i k when the dimension of the firework learns from another individual as follows.
Firstly, two individuals (excluding the selected firework X i ) 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., a i k = i in all dimensions; due to that, all generated random numbers above are less than P i ), we will randomly choose one dimension to learn from the corresponding dimension of another individual.The detailed procedure of choosing a i in (42) is given in Figure 3, where the ceil operator obtains the smallest integer greater than or equal to the specified expression and ps = n + m − 1 is the population size.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  14) and multiple constraints as ( 8)~( 13) to obtain the fitness value of each firework as (41) Initialize the adaptive amplitude of the optimal firework A ⁎ = X max − X min Calculate the number of explosion sparks S i of the firework X i as ( 24)~( 25) Calculate the explosion amplitude A i (except for A ⁎ ) of the firework X i as (25) Do all fireworks generate explosion sparks?
Generate S i explosion sparks of the firework X i as ( 27) and ( 29) Randomly select a firework and choose the corresponding mutation individuals as Figure 3 Generate a mutation spark X j for the selected firework as ( 42) and ( 29) Are all mutation sparks generated?
Evaluate the fitness of all individuals (fireworks, explosion sparks, and mutation sparks) as described above and calculate A ⁎ as ( 30)~( 31) Choose the best individual and other n − 1 randomly selected individuals as the next-generation fireworks Is the termination criteria met?
Choose the location of the optimal individual to be the control parameters of the reentry vehicle dynamics system  In order to address optimization problems in a general manner, each firework has a different mutation probability P i .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.

44
where h = r − r e and r e = 6378 km is the Earth's radius.Hence, r 0 = 6448 km and r f = 6408 km.
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 G max = 100000 for all algorithms.The AOA discrete point number N α = 5, and bank angle discrete point number N σ = 10.Besides, other main parameters are shown in Table 1.
The hard path constraint limits and control boundaries remain fixed throughout all simulations: Q s max = 5000 kW/m 2 , q max = 30 kPa, In order to deal with multiple constraints, these six algorithms all adopt the constraint-handling technique proposed in this paper.The simulation of reentry trajectory optimization is conducted on two cases using different objective functions: minimum downrange and minimum total heat.4.1.Minimum Downrange Trajectory.For convenience of description, the definition of the standard reentry longitudinal plane and related concepts is given.The schematic diagram of the downrange is presented in Figure 5, where 10 International Journal of Aerospace Engineering r 0 and v 0 are the position vector and velocity vector, respectively, of the reentry point.O E is the center of the Earth, and r f is the position vector of the reentry terminal point.The standard reentry longitudinal plane is a plane composed of r 0 and v 0 , and r b is the projection of r f on this plane.The vector n is the unit normal vector of the plane.The motion in the standard longitudinal plane is called the longitudinal motion, and the motion deviating from the plane is called the lateral motion.Therefore, the reentry downrange Z is defined as the arc length of a circle with the radius of r f corresponding to the angle δ Z between vectors r b and r 0 , which is formulated as follows: where the relevant parameters are defined as r 0 = r 0 ⋅ cos ϕ 0 cos θ 0 , cos ϕ 0 sin θ 0 , sin ϕ 0 , − cos θ 0 sin ϕ 0 cos γ 0 cos ψ 0 , sin θ 0 cos ϕ 0 sin γ 0 + cos θ 0 cos γ 0 sin ψ 0 − sin θ 0 sin ϕ 0 cos γ 0 cos ψ 0 , sin ϕ 0 sin γ 0 + cos ϕ 0 cos γ 0 cos ψ 0 , 47 where r 0 , θ 0 , ϕ 0 , v 0 , γ 0 , and ψ 0 represent the initial states of the reentry phase.
In the minimum downrange reentry trajectory optimization, the terminal constraints of the longitude, latitude, and heading angle need not be considered, and the objective is given The average convergence curves of the six comparative algorithms are illustrated in Figure 6.As we can see, in the minimum downrange reentry trajectory optimization, both the convergence speed and the final optimization result of our I-AFWA are better, when compared with those of the other five algorithms.Figure 7 presents the simulation results under running the I-AFWA 20 times, including the 3D trajectory, ground trajectory, control curves, and path constraint curves.It can be found that all constraints are satisfied and the reentry trajectory and control profiles are smooth throughout the flight, which indicate that our algorithm is quite efficient in handling multiple constraints.
The mean and standard deviation of the optimization performance index, path constraints, and terminal constraints of minimum downrange reentry trajectory obtained by the six intelligent algorithms are presented in Table 2.The performance index of the minimum downrange generated by the I-AFWA is much better than that generated by the S-AFWA, DE, PSO, GA, and FWA, which demonstrates that the optimization precision of the I-AFWA is higher and our improvement strategy is effective.Besides, the distribution range of the minimum downrange is very small, which shows that the I-AFWA achieves better performance and robustness.
The minimum downrange of the reentry vehicle is 1151.8km, and the corresponding trajectory has a larger AOA, which can increase the aerodynamic drag force and decrease the lift-drag ratio, thus consuming more energy and reducing the maneuver capacity of reentry vehicles, resulting in smaller flight downrange.

Minimum Total Heat
Trajectory.The objective function of minimum total heat reentry trajectory optimization problem is given as follows: Figure 8 presents the average convergence curves of the comparative algorithms.It can be seen that the I-AFWA also has the fastest convergence speed in the six algorithms.The simulation results using the I-AFWA, including the state and control profiles, are depicted in Figure 9.The minimum total heat and constraints of minimum total heat reentry trajectory with the six algorithms are compared in Table 3.The terminal constraints, such as the final altitude, longitude, latitude, velocity, flight path angle, and heading angle, are satisfied accurately.The heating rate, dynamic pressure, and aerodynamic load are less than the maximum allowable values strictly.The minimum total heat of the I-AFWA is much smaller than that of the other five algorithms, which proves once again that the optimization results of the I-AFWA are better.Additionally, the variance of the performance index is also relatively small.
From the simulation results, it can be seen that the minimum total heat of this reentry scenario is about 1.47 × 10 6 kJ, maybe smaller in fact.The reentry trajectory with minimum total heat corresponds to larger AOA and smaller bank angle, where the aerodynamic lift force of the vehicle is relatively large, the flight altitude is higher, and the atmospheric density is lower, thus reducing the heating rate (in terms of ( 8)) and decreasing the total heat; when the velocity is lower (close to the end of the reentry phase), pull down the AOA, so that the flight height drops rapidly, and meanwhile lateral maneuver is carried out to consume excess energy, making the vehicle reach the reentry terminal states and complete the flight mission.
Tables 2 and 3 also give the run time of each algorithm in the last row.On the two objective functions, we can see that the run time of the I-AFWA is typically  12 International Journal of Aerospace Engineering a little more than that of the S-AFWA, which is mainly because the I-AFWA takes more computational time on the selection of auxiliary mutation individuals.Note that the run time of PSO is the smallest and that of DE is the largest.However, for reentry trajectory optimization problem, compared with the long flight time, the little gap between the run times of the I-AFWA and S-AFWA is almost negligible.Although PSO runs much faster, the optimization result is relatively poor; DE is just the opposite.The running speed of the FWA is second only to that of PSO, but its optimization result is the worst.GA presents relatively poor performances in both two aspects.In summary, the results of the numerical simulations demonstrate that the proposed new mutation strategy, which enhances the information interactions between the fireworks and sparks to a great extent, can greatly improve the population diversity and thus avoid premature convergence during the search.Consequently, the I-AFWA can achieve a much better balance of exploration and exploitation than S-AFWA, DE, PSO, GA, and FWA and thus improve convergence speed and solution accuracy effectively.

Conclusion
In this paper, based on the AFWA and a new mutation strategy, a new optimization algorithm called the I-AFWA is proposed to generate an optimal trajectory for a hypersonic reentry vehicle under multiple constraints.By incorporating the new mutation operator and auxiliary mutation individual selection strategy to the S-AFWA, the information sharing among the fireworks and sparks is enhanced to generate more diverse individuals, and thus the global search capability is improved.Additionally, based on the distance measure and adaptive penalty function, a constraint-handling technique is developed to deal with all path and terminal constraints.The I-AFWA-based optimal reentry trajectory generation method takes all control variables into consideration, which makes it have high fidelity.
The simulation results of two different reentry missions are used to validate the reentry trajectory optimization algorithm, and some analyses are conducted.It is indicated that the I-AFWA exhibits better search efficiency and achieves better optimization performance than S-AFWA, DE, PSO, GA, and FWA.Besides, the convergence speed of the I-AFWA is also the fastest.The I-AFWA can also be applied to other optimization functions and other optimization problems.However, due to the increased computational cost caused by the selection process of auxiliary mutation individuals, the running speed of the I-AFWA is not the fastest one among all the algorithms.In the future research, the proposed algorithm will be further studied and tested on more complicated reentry scenarios involving waypoint and no-fly zone constraints.

Figure 2 :
Figure 2: The sketch map of altitude-velocity reentry corridor.

Figure 3 :
Figure 3: Selection procedure of auxiliary mutation individuals for the firework X i in all dimensions.
Figure out the objective function as (14) and multiple constraints as (8)~(13) to obtain the fitness value of each firework as (41)

Figure 4 :
Figure 4: Reentry trajectory optimization based on the I-AFWA with a new mutation strategy.
de θ (d eg ) L a t i t u d e 휙 ( d e g ) Altitude h (km)

Figure 7 :
Figure 7: Simulation results of 20 times minimum downrange reentry trajectory optimization using the I-AFWA.

Figure 8 :
Figure 8: The average convergence curves of the fitness function for each algorithm in minimum total heat reentry trajectory optimization.

Figure 9 :
Figure 9: Simulation results of 20 times minimum total heat reentry trajectory optimization using the I-AFWA.

Table 2 :
Performance index and constraints in minimum downrange reentry trajectory optimization with the six comparative algorithms.

Table 3 :
Performance index and constraints of minimum total heat reentry trajectory optimization with the six comparative algorithms.