Design of Optimal Attack-Angle for RLV Reentry Based on Quantum Particle Swarm Optimization

Abstract

The attack-angle optimization is a key problem for reentry trajectory design of a gliding type reusable launch vehicle (RLV). In order to solve such a problem, the equations of motion are derived first. A physical programming (PP) method is briefly presented and the preference function is reflected in mathematical representation. The attack-angle optimization problem with four criteria (i.e., downrange, total heat, heat rate, and trajectory oscillation) is converted into a single-objective optimization problem based on the PP method. A winged gliding reentry RLV is chosen as a simulation example and the transformed single-objective problem is solved by the quantum-behaved particle swarm optimization (QPSO) algorithm based on two types of preference structures, longer range preference and smaller total heat preference. The constraints of maximizing heating rate, normal load factor, and dynamic pressure and minimizing terminal velocity are handled by a penalty function method. The simulation results demonstrate the efficiency of these methods. The physical causation of the optimal solution and the typical profiles are presented, which reflect the designer's preference. At last, the feasibility and advantages of QPSO are revealed by comparison with the results of genetic algorithm (GA) and standard particle swarm optimization (PSO) algorithm on this optimization problem.

1. Introduction

A reusable launch vehicle (RLV) is a launch system which is capable of launching payloads, astronauts, or other equipment into space more than once. It plays an important role in the international space station (ISS) construction, materials transportation, satellite maintenance, and space experiments [1]. However, no true orbital RLV is currently in use. The closest example is the partially reusable Space Shuttle. Because of its expensive cost, the Shuttle retired in 2011. The development of RLV is driven by the desire to reduce expenditure. In order to provide better space transportation products and services, it is necessary to accelerate the development of new RLVs and related technologies and improve the reliability of RLVs during launch, on-orbit flying, and reentry. Moreover, this efficient transport system would lay the foundations for space technology applications in commercial and civilian areas.

RLV is a kind of multipurpose spacecraft that can be reused and flown freely between ground and space. Its total flying course can be divided into several phases: ascend, orbit injection, on-orbit, de-orbit, and reentry. A RLV is by far the most complex aircraft, involving advanced propulsion technology, lightweight materials and structure technology, aerodynamic heating technology, thermal protection, and guidance and control technology. The reentry trajectory and guidance design technologies play a vital role in the reentry process [2].

In conventional reentry trajectory optimization and tracking guidance strategy design of RLV, the law of attack-angle is generally given according to engineering experience. The vehicle usually flies at large attack-angle in early reentry phase to ease the burden of thermal protection structure system and then turn to attack-angle of maximum lift-to-drag ratio (L/D) to increase its lateral and longitudinal mobility. When the sideslip angle is ignored and the magnitude of lift force is fixed, the direction of total lift can be regulated by changing the bank angle to control the trajectory. This method has been applied in the Space Shuttle, X-33, and some other spacecrafts [3, 4].

Reentry control is presented by Roenneke and Markl [5] by tracking the drag-energy curve, a reference trajectory generated by sequential quadratic programming (SQP) algorithm. The vehicle's attack-angle is controlled to a Mach-dependent schedule. A methodology, which considers common equality and inequality constraints, is developed for fast design of reentry trajectory [6]. This approach makes a novel use of the quasiequilibrium glide phenomenon in lifting entry for effective enforcement of the inequality constraints. The application of orthogonal arrays and related statistical analysis of the results to reentry trajectory optimization are addressed by Mooij [7]. The main focus is put on the sensitivity of selected performance indices to the variation in the control parameters attack-angle and bank-angle. The integration of guidance and control system of an unpowered, winged reentry vehicle is discussed by Mooij et al. [8] and the entry guidance system commands attack-angle and bank-angle to steer the vehicle to the prescribed final state. To generate a reference trajectory, a reference attack-angle profile is provided as a function of normalized energy. The results of a feasibility study are presented by Chu et al. [9] about how internal analysis can be applied to the nonlinear reentry trajectory optimization problem of minimizing the heat load while returning to a desired landing point. The two control inputs are attack-angle and bank-angle. In order to minimize the maximum stagnation heat flux, the first part of reentry should be flown at high attack-angle (maximum lift configuration), and to maximize the cross-range and downrange, the flight should be continued at the condition of maximum L/D. The optimal trajectory planning for RLVs is studied by Zhang et al. [10] based on SQP. The cost function is the minimum accumulated heat load or the maximum downrange under multiple constraints. A migrant particle swarm optimization (PSO) algorithm, which employs stochastic search method and adaptive linear search method, respectively, for PSO search spaces, is used by Xie et al. [11] to solve the trajectory optimization problem in the presence of path constraints. A modified ant colony optimization (ACO) algorithm is applied by Zhang et al. [12] to optimize the drag-energy profile and its cost function is to minimize the total heat load for the RLV reentry trajectory. In [13], genetic algorithm (GA) also serves as a useful optimization tool in solving the trajectory optimization problem with minimum terminal states error. All of these studies choose bank-angle as the control variable and fix the profile of attack-angle.

However, presetting the pattern of attack-angle can result in intense air friction at the beginning of reentry. This will decrease the velocity at the expense of increasing the total amount of heat. Moreover, it is difficult for the control system to change the direction of total lift to steer the vehicle to the destination under fixed attack-angle profile. The conclusion can be drawn that the traditional design of attack-angle has severely limited the development of RLV—it is a prerequisite of RLV reentry trajectory optimization. The main task of this paper is to design a more rational attack-angle sequence under the constraints of heating rate, normal load factor, and dynamic pressure.

Although the optimization design of attack-angle has been proposed by some researchers, most of them focus just on one optimization criterion such as maximum downrange or minimum total heat load. Though some researchers have taken several criteria into account simultaneously, they mostly converted the multiobjective optimization (MOO) problem into single-objective optimization (SOO) problem by using conventional weight-based methods. The motivation and major concern of this work is that these simple approaches bring little flexibility and cannot meet the engineering requirements. In the initial design of reentry trajectory for RLVs, the engineers or designers are usually required to choose the sequence of attack-angle according to their preference based on multiple criteria.

This paper focuses on the optimization of attack-angle for gliding-type RLV reentry trajectory in the concept design phase and studies the optimal reentry attack-angle scheme, which can be applied directly to the study on the three-dimensional trajectory optimization and guidance design. When comparing to the ballistic reentry case, the gliding reentry can increase the cross-range and downrange. However, the aerodynamic heating is quite serious and the total heat and peak heat flux also bring enormous challenges to the structure design. Furthermore, the oscillation of the reentry trajectory is hoped not to be too large in the view of trajectory stability. So the studies on multiobjective optimization should be carried out. Along with the development of nonlinear programming (NLP) methods, a variety of MOO techniques have been developed over the past few years, such as the methods based on weighted approach, multiobjective evolutionary algorithms, multiobjective simulated annealing algorithms, and multiobjective PSO algorithms [14 –19]. In addition, MOO methods have been applied in the optimization design of reentry trajectory. A nondominated sorting genetic algorithm (NSGA) is used by Chen et al. [20] to the RLV multiobjective reentry optimization design with minimum heat load and maximum maneuverable range. In [21], the traditional skipping trajectory optimal control problem is discretized into a NLP problem using a direct parameter optimization method. Then the NLP problem is solved by a direct optimization algorithm. The constraints of the terminal state, heating rate, overload, and dynamic pressure are all taken into account. The performance index is selected to minimize the total heat, overload, and dynamic pressure. The control law of attack-angle is also given in advance.

In this paper, four optimization criteria (maximum down range, minimum heat load, minimum peak heat flux, and minimum trajectory oscillation) are considered, and the attack-angle optimization problem is brought into physical programming (PP) framework. Section 1 establishes the equations of motion for RLV reentry. The PP approach is introduced in Section 2. By using the PP method, the multiobjective attack-angle optimization problem is transformed into a single-objective optimization problem in Section 3. Then the transformed single-objective problem is solved by the QPSO algorithm in Section 4, in which the simulation results are also presented and compared with the standard PSO algorithm and GA. Section 5 concludes the work.

2. Problem Statement

Since the attack-angle optimization problem is generally carried out in the concept design phase of RLVs, the Earth can be considered as a symmetrical sphere and its own rotation effect can be ignored to simplify the problem. As discussed in [22], the amount of loss of accuracy is about 4.1% when neglecting the oblateness and rotation of the Earth. The motion parameters of the RLV are displayed in Figure 1. The magnitude of vector r indicates the distance from the Earth's center to the RLV, while λ and ϕ denote the longitude and latitude, respectively. All of the above three parameters determine the position of the RLV in the coordinates. The velocity can be characterized by the velocity vector v , heading angle ψ, and flight path angle θ. The last two parameters can be derived by projection of vector v to the local level.

Figure 1:

Depiction of RLV reentry flight.

The sideslip angle β which relates to the rotation of the RLV centerline from the relative wind is generally neglected in the reentry trajectory design and the thrust is assumed zero as the vehicle considered in this work is a gliding vehicle; then a nondimensional representation of the longitudinal equations of motion is as follows [23]:

\begin{matrix} \dot{v} = - D - \frac{\sin θ}{r^{2}}, \\ \dot{θ} = \frac{L}{v} - \frac{\cos θ}{v r^{2}} + \frac{v \cos θ}{r}, \\ \dot{r} = v \sin θ, \\ \dot{Φ} = \frac{v \cos θ}{r}, \end{matrix}

(1)

where v and r are the nondimensional velocity and distance from the Earth center. θ and Φ denote the path angle and range angle, respectively. The equations are normalized so that a change in one unit of each state element is approximately of equal significance. The dimensional units are divided by scale length R₀, scale time $\sqrt{R_{0} / g_{0}}$ , scale velocity $V_{c} = \sqrt{g_{0} R_{0}}$ , and scale acceleration g₀. R₀ is the average radius of the Earth and g₀ is the gravity acceleration of the sea level. The nondimensional time $τ = t / \sqrt{R_{0} / g_{0}}$ . Here L and D denote the nondimensional lift and drag acceleration which can be expressed as follows:

\begin{matrix} L = \frac{ρ S C_{L} {(v \cdot V_{c})}^{2}}{2 m g_{0}} = \frac{ρ S C_{L} v^{2} R_{0}}{2 m}, \\ D = \frac{ρ S C_{D} {(v \cdot V_{c})}^{2}}{2 m g_{0}} = \frac{ρ S C_{D} v^{2} R_{0}}{2 m}, \end{matrix}

(2)

where m and S are the mass and reference area of the vehicle and ρ denotes the atmospheric density. The aerodynamic coefficients, C_L and C_D, can be seen as functions of attack-angle approximately [24]:

\begin{matrix} C_{L} = - 0.041065 + 0.016292 α + 0.00026024 α^{2}, \\ C_{D} = 0.080505 - 0.03026 C_{L} + 0.86495 C_{L}^{2}, \end{matrix}

(3)

where α is the attack-angle which is ready to be optimized and determined.

During the reentry, a RLV is facing up with severe constraints of aerodynamic heating, load, dynamic pressure, and downrange. How to search for an optimal reentry trajectory which not only meets the design requirements but also satisfies multiple constraints is a difficult but meaningful work to be carried out.

In conventional trajectory optimization, as mentioned before, the attack-angle command is usually chosen as function of velocity or Mach number according to engineering experience, which lacks theoretical support. Considering the defect of traditional design, a new approach based on the PP method and QPSO algorithm is presented in this paper. In the next several sections, the PP method is introduced first and then applied to the transformation of the multiobjective optimization problem into a single-objective optimization model, which can reflect the preference of the designer. Later the QPSO algorithm is used to solve the optimization model and search for the optimal attack-angle sequence, which can satisfy the performance indexes and the constraints, simultaneously. Simulation studies and result comparisons are shown at last. The PP method for multiobjective problem will be firstly illustrated in the next section.

3. Physical Programming Method

The physical programming method that is suitable for generating a preferred compromise solution during multiobjective system design and optimization is described here. The initial development of the PP methodology is presented by Messac [25]. Since then, the PP method has been used in structure design, interactive design, robust design, and some other areas [26 –29].

3.1. Outline of Physical Programming Method

Physical programming is not a novel optimization algorithm, but a new framework dealing with multiobjective optimization problems. This new approach is intended to substantially reduce the computational burden for large problems and to place the design process into a more flexible and natural framework.

The PP method has two key advantages.

In stark contrast to conventional weight-based methods, once the designer's preferences are articulated, obtaining the corresponding optimal design is a noniterative process.

The PP method can find the optimal solution without much computational burden in contrast to Pareto-based multiobjective optimization methods. So the PP method is appropriate for the attack-angle optimization problem mentioned above in the concept design of RLV.

3.2. Mathematical Representation of Reference Function

In PP methods, the designer classifies the objective functions into four different classes, each comprising two cases: the Soft (S) case and the Hard (H) case:

Class 1-S: Smaller-is-Better, that is, Minimization,

Class 2-S: Larger-is-Better, that is, Maximization,

Class 3-S: Center-is-Better, that is, Seek Center Value,

Class 4-S: Range-is-Better, that is, Seek Range.

These classes refer to the behavior of the preference function with respect to each generic criterion. A lower value of the class function is considered better than a higher value. The ideal value of the class function during optimization is zero.

For each criterion, regions are defined to delineate the degrees of desirability: unacceptable, highly undesirable, undesirable, tolerable, desirable, and highly desirable. Figure 2 illustrates the case of Class 1-S [25, 30].

Figure 2:

Class 1-S regions for ith generic criterion.

The parameters g_i1 through g_i5 are physically meaningful values that are specified to quantify the preference functions associated with the ith objective. The mathematical implications of the above definitions and the development of the class function for each objective are given in [25]. Once the range parameters are defined for each objective function, the preference functions will be constructed. The following equations define the class function (Class 1-S) for the ith criterion in the domain of the kth region:

{\bar{g}}_{i} = {\begin{cases} {\bar{g}}_{i 1} \cdot \exp [(\frac{S_{i 1}}{{\bar{g}}_{i 1}}) (g_{i} - g_{i 1})], (k = 1) \\ T_{0} (ξ_{i k}) {\bar{g}}_{i (k - 1)} + T_{1} (ξ_{i k}) {\bar{g}}_{i k} + {\bar{T}}_{0} (ξ_{i k}, λ_{i k}) S_{i (k - 1)} \\ + {\bar{T}}_{1} (ξ_{i k}, λ_{i k}) S_{i k} (k = 2,3, 4,5), \end{cases}

(4)

where

\begin{matrix} ξ_{i k} = \frac{g_{i} - g_{i (k - 1)}}{g_{i k} - g_{i (k - 1)}}, λ_{i k} = g_{i k} - g_{i (k - 1)} (k = 2,3, 4,5), \\ T_{0} (ξ) \equiv \frac{1}{2} ξ^{4} - \frac{1}{2} {(ξ - 1)}^{4} - 2 ξ + \frac{3}{2}, \\ T_{1} (ξ) \equiv - \frac{1}{2} ξ^{4} + \frac{1}{2} {(ξ - 1)}^{4} + 2 ξ - \frac{1}{2}, \\ {\bar{T}}_{0} (ξ, λ) \equiv λ [\frac{1}{8} ξ^{4} - \frac{3}{8} {(ξ - 1)}^{4} - \frac{1}{2} ξ + \frac{3}{8}], \\ {\bar{T}}_{1} (ξ, λ) \equiv λ [\frac{3}{8} ξ^{4} - \frac{1}{8} {(ξ - 1)}^{4} - \frac{1}{2} ξ + \frac{1}{8}], \end{matrix}

(5)

where ${\bar{g}}_{i k}$ and S_ik (k = 1, 2, 3, 4, 5) denote the displacements and slopes at the region boundaries, respectively. The calculation of the displacements and slopes can be found in [25, 30].

From (4), the class function regions are represented by two generic functions. The first one is represented by a decaying exponential, which is fully defined by its value and slope at a point (region 1). For regions 2, 3, 4, and 5, the generic function takes the form of a spline segment that can be defined by its value and slope at its left and right boundaries.

It is obvious that the effectiveness of the PP method rests substantially on its ability to create class functions that truly reflect the priorities of the designer and the preceding classification offers significantly more flexibility than the typical weighted-criterion approach. For each criterion, a class function is formed that constitutes a component of the aggregate preference function to be minimized.

3.3. Physical Programming Problem Model

Building on previous development, the PP problem model with aggregate preference function is as follows [25]:

\begin{matrix} \min f (x) = lg {\frac{1}{n_{s c}} \sum_{i = 1}^{n_{s c}} {\bar{g}}_{i} [g_{i} (x)]} \\ s . t . g_{i} (x) \leq g_{i 5}, x_{\min} \leq x \leq x_{\max}, \end{matrix}

(6)

where x is the variable vector to be determined and x _max and x _min denote the maximum and minimum values of x , respectively.

The preceding problem model conforms to the framework of most NLP codes, with possible minor rearrangements. The PP method involves converting a multiobjective problem into a single-objective problem by using preference function that captures the designer's preferences. Before searching for a tradeoff solution, a single-objective optimization model for attack-angle should be established based on the PP method, which is the main work of next section.

4. Establishment of Optimization Model

As mentioned before, a proper optimization model for attack-angle is a precondition and foundation for optimal attack-angle solving. The single-objective model with multiple constraints and multiple criteria is derived as follows using PP method.

4.1. Optimization Criteria

In traditional trajectory optimization and guidance design, the attack-angle sequence is generally given based on the consideration of aerodynamic heating and horizontal mobility. However, more constraints and criteria should be treated theoretically in designing the attack-angle for reentry trajectory optimization. The unpowered RLVs usually adopt lifting body configurations to increase the range. In addition, severe conditions such as heating rate, total heat, and trajectory oscillation should also be considered. Here, the optimization criteria are chosen as

maximum down range: min f₁ = − Φ(τ_f),

minimum total heat: min $f_{2} = Q = \int_{τ_{0}}^{τ_{f}} \dot{Q} d τ$ ,

minimum heating rate: min $f_{3} = {\dot{Q}}_{\max}$ ,

minimum trajectory oscillation: min $f_{4} = {\dot{θ}}_{\max}$ .

4.2. Constraint Conditions

For the thermal, structural, and operational consideration, such severe constraints as heating rate $\dot{Q}$ , normal load factor n, and dynamic pressure q should be also taken into account [24]:

\begin{matrix} \dot{Q} (τ) = \frac{k_{s}}{\sqrt{R_{n}}} {(\frac{ρ}{ρ_{0}})}^{0.5} {(v)}^{3.15} \leq {\dot{Q}}_{\max}, \\ n (τ) = \sqrt{L^{2} + D^{2}} \leq n_{\max}, \\ q (τ) = \frac{1}{2} ρ {(v \cdot V_{c})}^{2} \leq q_{\max}, \end{matrix}

(7)

where k_s is a constant of the vehicle, R_n is the curvature radius of the stagnation point area, and ρ₀ is the atmospheric density at sea level. ${\dot{Q}}_{\max}$ , n_max, and q_max denote the corresponding maximum of the constraints of heating rate, normal load factor, and dynamic pressure, respectively. Their values are given in the section on simulation study.

The terminal conditions of the reentry are chosen to meet the specific terminal height and minimum terminal velocity:

\begin{matrix} h (τ_{f}) = h_{f}, \\ v (τ_{f}) \geq v_{f \min}, \end{matrix}

(8)

where h_f is the expected terminal height and v_fmin the minimum terminal velocity.

By using a penalty function method, this constrained optimization problem can be transformed into an unconstrained optimization problem and then the constraints with different physical meanings are converted to the values sharing the same order of magnitude [31]. The penalty term should be added into the performance index as long as the constraints are unsatisfied. The penalty terms for each constraint are as follows:

\begin{matrix} p_{1} (x) = {\begin{cases} 0 & v (τ_{f}) < v_{f \min} \\ {(\frac{v (τ_{f}) - v_{f}}{1000})}^{2} & v (τ_{f}) < v_{f \min}, \end{cases} \\ p_{2} (x) = {\begin{cases} 0 & {\bar{\dot{Q}}}_{\max} \leq {\dot{Q}}_{\max} \\ {(\frac{{\bar{\dot{Q}}}_{\max} - {\dot{Q}}_{\max}}{1000})}^{2} & {\bar{\dot{Q}}}_{\max} > {\dot{Q}}_{\max}, \end{cases} \\ p_{3} (x) = {\begin{cases} 0 & {\bar{n}}_{\max} \leq n_{\max} \\ {({\bar{n}}_{\max} - n_{\max})}^{2} & {\bar{n}}_{\max} > n_{\max}, \end{cases} \\ p_{4} (x) = {\begin{cases} 0 & {\bar{q}}_{\max} \leq q_{\max} \\ {(\frac{{\bar{q}}_{\max} - q_{\max}}{1 0^{5}})}^{2} & {\bar{q}}_{\max} > q_{\max}, \end{cases} \end{matrix}

(9)

where ${\bar{\dot{Q}}}_{\max}$ , ${\bar{n}}_{\max}$ , and ${\bar{q}}_{\max}$ are the maximum values of heating rate, normal load factor, and dynamic pressure during reentry, respectively.

4.3. Design Variables

It can be seen from (3) that the control variable is attack-angle α(τ). In order to find the optimal solution, the optimal control problem should be converted into a parameter optimization problem. A direct shooting method is chosen to perform the conversion [32].

The conversion begins with the definition of N + 1 fixed times τ₀ < τ₁ < ··· < τ_{N − 1} < τ_N = τ_f. The unknown control parameters α₀, α₁,…, α_{N − 1}, α_N at each point are to be determined. In each interval ∈ [τ_i, τ_{i + 1}] (i = 0, 1,…, N − 1), the control history α(τ) is formed by linear interpolation; that is,

α (τ) = α_{i} + \frac{τ - τ_{i}}{τ_{i + 1} - τ_{i}} (α_{i + 1} - α_{i}),

(10)

where α_i and α_{i + 1} denote the values of attack-angle at τ_i and τ_{i + 1}. By using the piecewise linear interpolation function, the unknown variable α(τ) is substituted by the discrete sequence α₀, α₁,…, α_{N − 1}, α_N. In addition, the final time is free, so the variable t_f should also be seen as a design parameter. Based on the conversion, the design parameters of the attack-angle optimization make up the vector u = (α₀, α₁,…, α_{N − 1}, α_N, t_f), which is going to be solved by optimization method.

4.4. Optimization Model

In the former development, four optimization criteria (i.e., maximum down range, minimum total heat, minimum heat rate, and minimum trajectory oscillation) are placed into the PP framework by the preference function. Based on the aggregate preference function (6) and equality/inequality constraints of heating rate, normal load factor, dynamic pressure, terminal height, and velocity, the integrated optimization index can be formed by adding the penalty terms of each constraint into the aggregate preference function. Then the optimization model of attack-angle is established as follows:

\begin{matrix} \min F (u) = lg {\frac{1}{n_{s c}} \sum_{i = 1}^{n_{s c}} {\bar{g}}_{i} [g_{i} (u)]} + \sum p_{j} (u) \\ s . t . g_{i} (u) \leq g_{i 5}, u_{\min} \leq u \leq u_{\max}, \end{matrix}

(11)

where u = (α₀, α₁,…, α_{N − 1}, α_N, t_f) is the design parameter vector to be solved. u _min and u _max are the boundary values of attack-angle.

As it is seen from preceding process, the multiobjective optimization problem of attack-angle design has been transformed into a single-objective optimization problem based on the PP method. The nonlinear optimization model in (11) will be solved by QPSO algorithm in the next section to obtain the optimal attack-angle sequence and flying time. The results will satisfy the designer's preference.

5. Attack-Angle Optimization Based on QPSO

The optimization model of attack-angle derived in previous section belongs to the class of large-scale nonlinear optimization problems, which can be solved by NLP methods. However, the traditional gradient-based approaches are very sensitive to initial guess and converge rapidly but may get stuck in local minima of multimodal functions [33, 34]. To improve the optimization performance, stochastic optimal approaches such as GA, ACO, PSO, and simulated annealing (SA) algorithms are proposed and some of which have been applied in trajectory optimization problems [35 –38]. Although stochastic optimization algorithms waste computational effort by doing a random search, they have better global search ability during the offline trajectory optimization and attack-angle design, which do not have much requirement on real-time environment.

5.1. PSO Algorithm

PSO algorithm, as a newly developed swarm intelligence paradigm proposed by Kennedy and Eberhart [39], has seen increased popularity recently in the field of optimization for many problems. PSO does not mimic Darwinian natural selection but rather it imitates the social behavior of birds.

PSO optimizes a problem by having a population (swarm) of candidate solutions (particles). These particles are moved around in the search space according to a few simple formulae. The movements of the particles are guided by their own best known position in the search space as well as the entire swarm's best known position. When improved positions are being discovered, the movements of the swarm will be guided. The process is repeated and by doing so it is hoped that a satisfactory solution will eventually be discovered.

PSO algorithm does not use the gradient of the problem being optimized, which means PSO does not require that the optimization problem be differentiable as it is required by classic optimization methods such as gradient descent and quasi-Newton methods. When compared to GA, ACO, SA, and some other stochastic optimization algorithms, the PSO algorithm presents relatively short optimization process and fewer control variables. Therefore the PSO algorithm is used for many kinds of optimization problems [40 –43].

5.2. QPSO Algorithm

PSO algorithm makes few or no assumptions about the problem being optimized and can search very large spaces of candidate solutions. However, PSO cannot guarantee a global optimal solution, which has been proved by van den Bergh [44, 45]. In order to enhance the chance to find out better solutions in search space, the individual particle of a PSO system moving in a quantum multidimensional space is studied by Sun et al. [46] and a quantum delta potential well model for PSO is established, and then the quantum delta-potential-well-based particle swarm optimization (QDPSO) algorithm is proposed. Based on the QDPSO algorithm, the mean best (m best) position is introduced by Sun et al. [47] to evaluate the search scope of a particle and put forward a novel parameter control method and thus formulate the philosophy of QPSO.

In PSO, the particle must be in bound state to guarantee collectiveness of the particle swarm and converge to optima or suboptima, whereas, in QPSO, a particle can appear at any position in the whole feasible search space with a certain probability and such position may be superior to the current best position in the population. In contrast to the traditional PSO, the QPSO has even simpler model and less control variables, which promote its application in many research fields [48, 49].

In a classical PSO system with the population size of M, the position vector of particle i in N-dimensional space in the tth iteration is denoted as

X_{i} (t) = {X_{i, 1} (t), X_{i, 2} (t), \dots, X_{i, N} (t)}, i = 1,2, \dots, M .

(12)

The particle in quantum space has no velocity vector. For minimization problems, a smaller value of the performance index corresponds to a better fitness value [21], so the best position (the position giving the best fitness value) of particle i can be updated by

P_{i} (t) = {\begin{cases} X_{i} (t) & if f [X_{i} (t)] < f [P_{i} (t - 1)] \\ P_{i} (t - 1) & if f [X_{i} (t)] > f [P_{i} (t - 1)] \end{cases}

(13)

and the best particle among all particles in the population can be determined by

\begin{matrix} g = \arg \min_{1 \leq i \leq M} {f [P_{i} (t)]}, \\ G (t) = P_{g} (t), \end{matrix}

(14)

where g is the suffix of the best particle.

Actually, before every update of the ith particle P_i(t), the best particle across the whole population G(t) should be known. Therefore, the current fitness value of the ith particle f[P_i(t)] should be compared with the fitness value of the ith particle's best position f[G(t)] and then update G(t) if f[P_i(t)] is better.

Each particle in the PSO system converges to its only local attractor, of which the coordinates are

p_{i, j} (t) = φ_{j} (t) \cdot P_{i, j} (t) + [1 - φ_{j} (t)] \cdot G_{j} (t), 1 \leq j \leq N,

(15)

where φ_j(t)∼U(0, 1) is a random number between 0 and 1.

In QPSO, the state of the particle in quantized search space is depicted by a wave function or by a probability function of the position, which cannot reveal any certain information about the position of a particle that is vital to evaluate its fitness. So it is necessary to transform the quantum state to classical state. A Monte Carlo method can realize the process of transformation. The iterative equation of position in QPSO is

X_{i, j} (t + 1) = p_{i, j} (t) \pm η | C_{j} (t) - X_{i, j} (t) | \cdot \ln [\frac{1}{u_{i, j} (t)}],

(16)

where u_i,j(t)∼U(0, 1) is a random number between 0 and 1. Here, η is called contraction-expansion (CE) coefficient, which is the only parameter in the algorithm except for the population size and iteration times, and works on individual particle convergence speed and the performance of the algorithm. In the simulation, the value of η decreases from 1.2 to 0.5 linearly, and C_j(t) is the mean best position (m best) defined as follows:

C_{j} (t) = \frac{1}{M} \sum_{i = 1}^{M} P_{i, j} (t) .

(17)

The QPSO algorithm is described in Algorithm 1.

Algorithm 1:

5.3. Simulation Procedure

Based on the PP method and QPSO algorithm, the attack-angle optimization design problem can be solved by the following procedure.

Design the preference structure of each optimization criterion according to the designer's preference and divide the total domain of each criterion into six subintervals as illustrated in Figure 2.

Determine the preference function of each criterion based on the preference structure. Get the boundary conditions (displacements and slopes) of the regions for each criterion, and then obtain the preference function value ${\bar{g}}_{i} (i = 1,2, 3,4)$ according to (4).

Substitute the preference function values into the aggregate preference function (6) to convert the multiobjective problem into a single-objective problem.

Establish the integrated optimization model of attack-angle by adding the penalty terms of constraints into the aggregate preference function, and then solve the model as shown in (11) using the QPSO algorithm.

5.4. Simulation Studies

In order to investigate the proposed method for the attack-angle optimization problem, a gliding type reentry RLV is chosen as a simulation example. The aerodynamic coefficients are generated based on the common aero vehicle (CAV) data [30, 50]. These coefficients are functions of attack-angle shown in (3) and the effect of the Mach number can be neglected according to the data. The initial conditions and model parameters are given in Table 1. This simulation task is carried out on a PC with 2 GB memory and Intel Core 2 Duo CPU. The version of MATLAB is R2007b.

Table 1:

Numerical values for simulation.

Variable	Sign	Value (unit)

Mass	m	907 (kg)
Reference area	S	0.35 (m²)
Maximum attack-angle	α_max	20 (deg)
Attack-angle of maximum L/D	α_{L/D max}	8.5 (deg)
Initial velocity	V ₀	7000 (m/s)
Initial height	h ₀	100 (km)
Initial path angle	θ₀	−0.2 (deg)
Terminal height	h _f	20 (km)
Minimum terminal velocity	V _fmin	1000 (m/s)
Maximum heating rate	${\dot{Q}}_{\max}$	2000 (kW/m²)
Maximum normal load factor	n _max	4 (g)
Maximum dynamic pressure	q _max	800 (kPa)

According to engineering experience and the simulation results of [30], the vehicle could achieve maximum range when flying at attack-angle of maximum L/D, but the trajectory oscillation is severe; when flying at the maximum attack-angle, the total time-cost is reduced and the heat load is minimum, but the peak heat flux increases and the range is too short. Thus, it is noticed that there exist multiple objectives in the attack-angle optimization problem which conflict with each other. Generally, it is difficult or impossible to find the best solution that all of the criteria achieve the best values simultaneously. Under this consideration, it is preferred to obtain the compromised solution based on the designer's expectation or preference.

Traditionally, the sequence of attack-angle is usually chosen as a simple profile of velocity or Mach number according to the flight demonstration and engineering experience in the trajectory optimization problems. In order to compare with the conventional design method, two groups of preference structure (Preferences I and II) that represent longer downrange expectation and smaller total heat expectation are chosen. The two preference structures are shown in Tables 2 and 3.

Table 2:

Preference I: longer downrange expectation.

Criterion (unit)	Preference type	Region boundaries
Criterion (unit)	Preference type	g _i1	g _i2	g _i3	g _i4	g _i5

f₁ (deg)	Class 1-S	−90	−85	−80	−75	−70
f₂ (10⁵ kJ/m²)	Class 1-S	6	10	14	18	22
f₃ (kW/m²)	Class 1-S	1200	1400	1600	1800	2000
f ₄	Class 1-S	1.2	2.0	2.8	3.6	4.4

Table 3:

Preference II: smaller total heat expectation.

Criterion (unit)	Preference type	Region boundaries
Criterion (unit)	Preference type	g _i1	g _i2	g _i3	g _i4	g _i5

f₁ (deg)	Class 1-S	−70	−65	−60	−55	−50
f₂ (10⁵ kJ/m²)	Class 1-S	2	5	8	11	14
f₃ (kW/m²)	Class 1-S	1200	1400	1600	1800	2000
f ₄	Class 1-S	1.2	2.0	2.8	3.6	4.4

The total time space is divided into 20 and 10 equal intervals, respectively, to illustrate the feasibility of the optimization method for the two preferences. Before simulation, the parameters of QPSO are chosen as follows: the population size is 40, the maximum iteration is limited to 500, the terminating value of the integrated index is 0.5, and the CE coefficient η decreases from 1.2 to 0.5 linearly. Since QPSO is a stochastic optimization algorithm, 50 independent optimization trials are carried out to guarantee the accuracy of simulation.

In order to provide complete information about the distribution of the observations, the five-number summary is introduced to give a robust description of the location (from the median), range (from the minimum and maximum), and spread (from the percentile) of the observations. The five-number summary of a data set consists of the minimum observation, 25th percentile, median, 75th percentile, and maximum observation of the set. The observations of the 50 independent trials are put in ascending order. Then the values of the five-number summary can be obtained. The results and typical curves are shown as follows.

N = 20. The simulation results and the spread of the fitness are shown by Table 4 and Figure 3; the typical iterative processes for two preferences are depicted by Figure 4; the typical curves are presented by Figures 5, 6, 7, 8, 9, and 10.

Table 4:

Simulation results of 50 trials (N = 20).

Results	Preference I	Preference II

Best value (min)	0.3081	0.3517
25th percentile	0.3749	0.4429
Median	0.4342	0.4679
75th percentile	0.4776	0.4919
Worst value (max)	0.4974	0.4985
Arithmetic mean	0.4318	0.4691
Typical results
f₁ (deg)	85.41	76.05
f₂ (10³ kJ/m²)	968.6	539.1
f₃ (kW/m²)	1.3429	1.0883
f₄	1.42	1.80

Figure 3:

The spread of the fitness values of 50 trials for the two preferences (N = 20).

Figure 4:

The typical iterative processes for the two preferences (N = 20).

Figure 5:

Typical profiles of height-range for reentry (N = 20).

Figure 6:

Typical profiles of velocity-range for reentry (N = 20).

Figure 7:

Typical profiles of optimal attack-angle for reentry (N = 20).

Figure 8:

Typical profiles of heat rate for reentry (N = 20).

Figure 9:

Typical profiles of normal load factor for reentry (N = 20).

Figure 10:

Typical profiles of dynamic pressure for reentry (N = 20).

N = 10. The simulation results and the spread of the fitness are shown by Table 5 and Figure 11; the typical iterative processes for two preferences are depicted by Figure 12; the typical curves are presented by Figures 13, 14, 15, 16, 17, and 18.

Table 5:

Simulation results of 50 trials (N = 10).

Results	Preference I	Preference II

Best value (min)	0.2048	0.2001
25th percentile	0.3188	0.3079
Median	0.3732	0.3664
75th percentile	0.4382	0.4347
Worst value (max)	0.4948	0.4949
Arithmetic mean	0.3818	0.3693
Typical results
f₁ (deg)	84.89	78.50
f₂ (10³ kJ/m²)	991.2	572.0
f₃ (kW/m²)	1.1843	1.1004
f₄	1.49	1.61

Figure 11:

The spread of the fitness values of 50 trials for the two preferences (N = 10).

Figure 12:

The typical iterative processes for the two preferences (N = 10).

Figure 13:

Typical profiles of height-range for reentry (N = 10).

Figure 14:

Typical profiles of velocity-range for reentry (N = 10).

Figure 15:

Typical profiles of optimal attack-angle for reentry (N = 10).

Figure 16:

Typical profiles of heat rate for reentry (N = 10).

Figure 17:

Typical profiles of normal load factor for reentry (N = 10).

Figure 18:

Typical profiles of dynamic pressure for reentry (N = 10).

Tables 4 and 5 show the simulation results of the 50 independent optimization trials for Preferences I and II, respectively. The smallest fitness value (the best results), 25th percentile, median, 75th percentile, and the largest fitness value (the worst results) are measured according to the observations of 50 independent trials. The arithmetic mean and the median are all given as measures of the average result. In the first case of Preference I, we have a median of 0.4342 and an arithmetic mean of 0.4318, so they are quite comparable. Even the worst fitness values are less than 0.5, which is the terminating value of the integrated index. The same holds for other cases and preferences. So the convergence of every iterative process for the attack-angle optimization problem can be guaranteed and all the results lie in an ideal scope.

Figures 3 and 11 represent the spread of the results. In every trial, the population size, the maximum iteration, and the terminating value of the integrated index are specified. So the simulation process stops and the fitness value (the integrated index) is obtained at the end. Nearly all of the trials got a different result under different initial values. As can be seen from the results, the QPSO algorithm shows good performance during the optimization design of attack-angle.

Figures 4 and 12 depict the typical iterative processes of QPSO for the two preferences. It is obvious that the terminating value of the integrated index can be reached within several steps by doing a random search, which indicates that QPSO has a fast convergence speed for attack-angle optimization problem.

Figures 5 to 10 and Figures 13 to 18 present the typical curves of the 50 independent optimization trials. It is noticed in these figures that all of the optimization criterion values are in the tolerable or better regions. By comparing the simulation results of the two groups of preference structures, the optimal attack-angle profiles are trade-off solutions under different considerations of preference structures. Figures 5 and 13 present the damped oscillation trajectories of gliding reentry, which are compromise results of every single optimization object as discussed by Yong et al. [30]. The trajectories with amplitude of damped oscillation are response to the consideration of increasing the range, while, in the view of trajectory stability, the oscillation is expected to be not so severe. Meanwhile, the total flying time should be reduced to meet the constraint of thermal protection system (TPS) at the expense of range decrease. Figures 5 and 13 can be understood easily based on these physical analyses. Figures 7 and 15 show the history of optimal attack-angle sequence corresponding to the optimal reentry trajectory. Because the peak heating flux usually occurs at the first trough, the RLV usually flies at high value of attack-angle in the initial part of the reentry trajectory. In the later segment of the trajectory, when a designer prefers to obtain an ideal trajectory with smaller amount of total heat, the attack-angle stays at the close region of maximum value to reduce the heat; while a designer has a bias towards longer downrange, the attack-angle varies near the vicinity of maximum L/D to increase the downrange. These characteristics coincide with the results of single-objective optimization problem. The histories of heating rate, normal load factor, and dynamic pressure under different preferences are depicted in Figures 8 to 10 and Figures 16 to 18, which satisfy the constraint conditions.

5.5. Results Comparison

The defects of the traditional gradient-based approaches to NLP problems are discussed in the beginning of this section. Meanwhile, many stochastic optimization algorithms are developed, such as GA, ACO, and PSO. Compared with gradient-based methods, these stochastic approaches are not sensitive to initial guess and present strong global search ability when solving complex NLP problems.

As a novel attempt, the PP method is firstly introduced into the process of solving multiobjective trajectory optimization by Yong et al. [30]. The multiobjective optimal reentry trajectory is obtained by GA. The similar optimization problem is solved by Arora et al. [51] based on a trajectory optimization software which uses the GA as an optimization tool. The total time consumed is about 18 minutes. The PP method and the traditional PSO algorithm are introduced into an aircraft concept design problem by Wang and Li [52]. The efficiency and feasibility of PSO are revealed by comparing with the results of GA. The optimal trajectory planning for RLV based on ACO is studied by Zhang et al. [53]. The cost function is the minimum accumulated heat load under multiple constraints such as heating rate, normal load factor, dynamic pressure, and equilibrium glide. A distributed architecture to deal with global optimization of engineering problems in a parallel computing environment is presented by Mooij and H $\ddot{a}$ nninen [54], in which a particular implementation with PSO algorithm is discussed and used for reentry trajectory optimization that minimizes the heat load and maximizes the total flight range simultaneously. A new maximum cross-range guidance method is presented by Chao et al. [55] and divided into two subproblems, the trajectory design and trajectory following subproblem. The optimal reference trajectory is found by GA base on the analysis of heat flux, overload factor, and dynamic pressure constraints.

Through comparison, GA is sensitive to the initial population and has poor ability for local search, and “prematurity” usually occurs. The distribution of the initial population and the selection of genetic operators directly concern the convergence and searching efficiency of GA. Unreasonable initial population will lead to local optimal solution, so the initial population of GA should be chosen scientifically and will reflect the information in the whole solution space. ACO also has an initial guess sensitivity and unavoidable stagnation phenomenon. When the population is initialized, the pheromone on each path is equal and it is difficult for the ants to find a better path among the rambling paths within short time. So the convergence speed is slow. Besides, Differential Evolution (DE) and Artificial Bee Colony (ABC) algorithms also have good performance in addressing numerical optimization problems and both of them are good candidates for attack-angle optimization problems. However, the research of DE is relatively scattered and lacks systematic and theoretical foundation when compared with GA and PSO. Moreover, appropriate control parameters like crossover rate and mutation rate should also be selected, which is quite similar to that of GA. The ABC algorithm is a newly emerging evolutionary optimization algorithm and the study of ABC algorithm for constrained optimization problems is still in primary stage. So much work needs to be done to improve the capabilities of classical ABC algorithm before applications widely, such as algorithm complexity, theoretical proofs, convergence capability, and easily sinking in local optima. A traditional PSO algorithm has many advantages like having a simple search policy and fewer control parameters, which are suitable for the attack-angle optimization design problem in the initial concept design phase. However, the search processes of these traditional intelligent optimization algorithms usually require more time. In order to test the performance of QPSO on attack-angle optimization, simulation studies are carried out and the results are compared with that of GA and PSO.

For the attack-angle optimization problem and optimization model that have been established before, total time domain is divided into 10 equal intervals, and then the traditional GA, PSO, and QPSO algorithms are used to solve this optimization problem for 50 times, based on Preference I (longer downrange). The population size is 40; the maximum iteration is limited to 500; the terminating value of the integrated index is 0.5.

It is well known that the choice of parameter values for GA, PSO, and QPSO is critical to the success of these algorithms. For example, the crossover rate controls the capability of GA in exploiting the local optima. The higher the crossover rate is, the quicker the exploitation proceeds. However, if the crossover rate is too large, it would disrupt individuals faster than they could be exploited. In contrast, the mutation rate controls the speed of GA in exploring a new area. Small mutation rate values are commonly adopted in GA. Typical values of crossover rate are in the range0.5∼1.0, while typical values of mutation rate are in the range 0.001∼0.05 [56]. In addition, PSO is also very sensitive to efficient parameter setting and many approaches have been proposed to select the values of inertia weight, cognitive, and social parameters for PSO to enhance and improve its performance. However, choosing the optimal parameter values for these algorithms is not the main work of our paper, so typical values are used in the simulation studies by setting p_c = 0.8 and p_m = 0.02 for GA, respectively. The inertia weight and learning effects of PSO are ω = 0.7298, c₁ = 1.49618, and c₂ = 1.49618. The CE coefficient η of QPSO decreases from 1.2 to 0.5 linearly.

As can be seen from previous discussion, the five-number summary provides a concise summary of the distribution of the observations. In addition, it is possible to quickly compare several sets of observations by comparing their five-number summaries, which can be represented graphically using a boxplot. A boxplot is a visual representation of the five-number summary. It can be used to provide important information about the location and dispersion of the data. So it is particularly useful for comparing distributions between different groups. The simulation results are reported in Table 6. The boxplots of the three algorithms are drawn vertically in Figure 19.

Table 6:

Results comparison of GA, PSO, and QPSO.

Results	GA	PSO	QPSO

Best value (min)	0.4383	0.3142	0.2048
25th percentile	0.6681	0.4187	0.3188
Median	0.8766	0.5468	0.3732
75th percentile	1.1839	0.6697	0.4382
Worst value (max)	1.5134	0.8346	0.4948
Arithmetic mean	0.9344	0.5537	0.3818
Typical results
f₁ (deg)	82.07	83.13	84.89
f₂ (10³ kJ/m²)	878.3	953.2	991.2
f₃ (kW/m²)	1.4736	1.2942	1.1843
f₄	2.43	1.87	1.49
Convergence time (s)	744.36	134.82	32.67
Average iteration	386.58	65.74	12.46

Figure 19:

Boxplots comparison of GA, PSO, and QPSO.

The fitness value is obtained in every trial when the population size, maximum iteration, and the terminating value of the integrated index are all specified. Almost all of the trials will get a different result under different initial values and different optimization algorithms. Table 6 shows the five-number summaries of the observations of GA, PSO, and QPSO algorithms. The best value, 25th percentile, median, 75th percentile, worst value, and arithmetic mean result under each algorithm are all measured according to the distribution of the results of 50 independent trials to illustrate the different performances of the three algorithms.

Figure 19 depicts the boxplots of the data in different groups. On each box, the central mark is the median and the edges of the box are the 25th and 75th percentiles. The upper and lower extremes are the greatest and least numbers that occur in the data groups. The spacings between the different parts of the box indicate the degree of dispersion and skewness in the data. Given the maximum iteration, most of results of GA are larger than 0.5, while the fitness values of PSO and QPSO are all around or smaller than 0.5. More importantly, all the observations of QPSO lie in a narrow scope, which indicates a better performance of the QPSO than GA and PSO on the attack-angle optimization problem.

In addition, a statistical study based on the Mann-Whitney U test is also carried out to highlight the differences in performance of the three different algorithms. The significant level is chosen as 0.05 and the P values of the tests are all approaching zero, which indicate a rejection of the null hypothesis at the 5% significance level and significant differences in performance of the algorithms.

As can be seen from the comparison of the results, all of the three optimization algorithms can achieve good solutions, but the QPSO algorithm is superior to GA and PSO algorithm in search ability, convergence speed, and the accuracy of the solution on the particular problem discussed in this paper. These characteristics can also be found in other quantum-based algorithms [57, 58]. First, when adopting the same data, the simulation results of QPSO are apparently better than the results of GA and PSO; this is because the particles in QPSO can appear at any position in the whole feasible search space. Second, the convergence speed of QPSO is faster and its optimization process is more efficient than GA and PSO, because the quantum computation was introduced into the stochastic procedure. Due to the characteristics of wave function of quantum delta potential well model, the global search ability of the algorithm is enhanced. Considering the advantages of QPSO and the requirement of fast design of attack-angle in engineering practice, QPSO algorithm could serve as a novel trial on attack-angle optimization problem in concept design of RLVs.

It is easy to see from the simulation results that the QPSO performs better on the attack-angle optimization problem, but sometimes the QPSO is sensitive to the control parameters, which will result in unstable convergence speed of the particle. In addition, fast convergence speed is sometimes acquired in sacrifice of global search ability of the algorithm, so effective methods should be found to improve the performance of QPSO algorithm.

6. Conclusions

The reentry trajectory optimization problem for a horizontal-landing reusable launch vehicle (RLV) is based on the successful choosing of attack-angle sequence. This paper studies the attack-angle optimization problem using physical programming (PP) method and quantum-behaved particle swarm optimization (QPSO) algorithm and obtains the optimal solution corresponding to the designer's preference.

First, the longitudinal equations of motion were introduced to depict the reentry course of the RLV, followed by the problem and its motivation. Then the PP method was briefly introduced and the mathematical representation of preference function was presented to deal with the multiobjective optimization problems. The attack-angle optimization problem with four criteria (maximum down range, minimum total heat, minimum heat rate, and minimum trajectory oscillation) was transformed into a single-objective optimization problem based on the PP method and the integrated optimization model was established at last. The attack-angle discrete sequence based on direct shooting method acts as the design parameters of the model; the constraints of maximum heating rate, maximum normal load factor, maximum dynamic pressure, and minimum terminal velocity were handled by penalty function method, and then the integrated optimization index was formed by adding the penalty terms of each constraint into the aggregate preference function. Finally, this single-objective optimization problem of attack-angle was solved by QPSO algorithm based on two types of preference structures (longer range preference and small total heat preference). In addition, the simulation results of QPSO are compared with the results of other algorithms, especially GA and standard PSO. The QPSO shows strong search ability and high efficiency on the attack-angle optimization problem because of the application of quantum computation.

The combination of the PP method and QPSO algorithm serves as a useful and efficient mean on the attack-angle optimization design problem for RLV reentry. According to this method, the engineers can establish the optimization model and obtain the trade-off solution that satisfies their preference within several minutes, which could further enhance the accuracy and performance of the reference trajectory optimization and guidance system design. In the future studies, the three-dimensional reentry trajectory optimization problem can be directly solved using this approach and obtain the optimal solution of attack-angle and bank-angle, which can be treated as the reference trajectory for a reentry guidance system. In addition to the algorithms discussed in this paper, some other algorithms like DE and ABC are all good candidates for the reentry attack-angle optimization and trajectory optimization problems. In the future, more advanced algorithms will be tried to handle this kind of problems.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Footnotes

Acknowledgement

This work is supported in part by the National Key Technology Research and Development Program (no. 2011BAK16B03).

References

Zhao

Cui

N. G.

Guo

J. F.

, and Fu

, “Reentry trajectory planning based on the secondary reversal dynamics for the second generation reusable launch vehicles,” Advanced Materials Research, vol. 625, pp. 109–112, 2012.

Fasano

and Pinter

J. D.

, Modeling and Optimization in Space Engineering, Springer, New York, NY, USA, 2013.

Lafleur

J. M.

and Saleh

J. H.

, “Survey of mission evolution and flexibility in the Space Shuttle program,” Space Policy, vol. 26, no. 4, pp. 236–245, 2010.

D. W.

, “Reentry guidance based on parametric optimization,” Applied Mechanics and Materials, vol. 278–280, pp. 1496–1499, 2013.

Roenneke

A. J.

and Markl

, “Re-entry control to a drag-vs-energy profile,” Journal of Guidance, Control, and Dynamics, vol. 117, no. 5, pp. 916–920, 1994.

Shen

Z. J.

and Lu

, “Onboard generation of three-dimensional constrained entry trajectories,” Journal of Guidance, Control, and Dynamics, vol. 26, no. 1, pp. 111–121, 2003.

Mooij

, “Parametric control-variable analysis in support of reentry trajectory optimisation,” in Proceedings of the Atmospheric Fright Mechanics Conference and Exhibit (AIAA '04), pp. 1–12, Providence, RI, USA, August 2004.

Mooij

Mease

K. D.

, and Benito

, “Robust re-entry guidance and control system design and analysis,” in Proceedings of the Guidance, Navigation, and Control Conference and Exhibit (AIAA '07), pp. 1–14, Hilton Head, SC, USA, August 2007.

Chu

Mooij

van Kampen

, and Chu

Q. P.

, “A feasibility study to the application of interval analysis to re-entry trajectory optimization,” in Guidance, Navigation, and Control Conference and Exhibit (AIAA '08), pp. 1–18, Honolulu, Hawaii, USA, August 2008.

10.

Zhang

Q. Z.

Gao

Guo

, and Ren

, “Reentry trajectory planning optimization based on sequential quadratic programming,” in Proceedings of the 2nd International Symposium on Systems and Control in Aerospace and Astronautics (ISSCAA '08), pp. 1–5, Shenzhen, China, December 2008.

11.

Xie

F. Q.

Wang

Y. J.

Zheng

Z. C.

, and Zhang

, “Constrained trajectory optimization using migrant particle swarm optimization algorithm,” in Proceedings of the 4th International Conference on Bio-Inspired Computing, pp. 1–7, Beijing, China, October 2009.

12.

Zhang

Guan

S. Y.

, and Chen

S. L.

, “Modified ant colony algorithm-based reentry trajectory optimization for suborbital launch vehicle,” Journal of Astronautics, vol. 32, no. 6, pp. 1244–1248, 2011.

13.

Tang

Chen

, and Xu

, “Reentry trajectory planning based on genetic optimization of terminal matching,” in Proceedings of the International Conference on Electronic and Mechanical Engineering and Information Technology (EMEIT '11), pp. 1151–1155, chn, August 2011.

14.

Boonlong

Kittipong

, “Multiobjective optimization of a vehicle vibration model using the improved compressed-objective genetic algorithm with convergence detection,” Advances in Mechanical Engineering, vol. 2013, Article ID 121495, 14 pages, 2013.

15.

Zhang

A. A.

and Yang

H. G.

, “Pareto-set based multi-objective dynamic reactive power and voltage control,” Asian Journal of Control, vol. 14, no. 4, pp. 1028–1035, 2012.

16.

Xie

, “Multi-objective H2/L2 performance controller synthesis for LPV systems,” Asian Journal of Control, vol. 14, no. 5, pp. 1273–1281, 2012.

17.

Ganguly

Sahoo

N. C.

, and Das

, “Multi-objective planning of electrical distribution systems using dynamic programming,” International Journal of Electrical Power & Energy Systems, vol. 46, pp. 65–78, 2013.

18.

Sadollah

Bahreininejad

Eskandar

, and Hamdi

, “Optimum material gradient for functionally graded dental implant using particle swarm optimization,” Advanced Materials Research, vol. 647, pp. 30–36, 2013.

19.

Shivasankaran

Kumar

P. S.

Nallakumarasamy

, and Raja

K. V.

, “Repair shop job scheduling with parallel operators and multiple constraints using simulated annealing,” International Journal of Computational Intelligence Systems, vol. 6, no. 2, pp. 223–233, 2013.

20.

Chen

Wan

Z. M.

, and Chen

S. L.

, “RLV reentry trajectory multi-objective optimization design based on NSGA-II algorithm,” in Proceedings of the Atmospheric Plight Mechanics Conference and Exhibit (AIAA '05), pp. 1–6, San Francisco, Calif, USA, August 2005.

21.

L. H.

Yuan

J. P.

Fang

, and Luo

J. J.

, “Reentry skipping trajectory optimization using direct parameter optimization method,” in Proceedings of the 14th AIAA/AHI International Space Planes and Hypersonics Systems Technologies Conference, pp. 1–14, Canberra, Australia, November 2006.

22.

Zheng

Y. Y.

Liu

, and Li

J. F.

, “Spacecraft gliding trajectory optimization,” Journal of Astronautics, vol. 30, no. 3, pp. 858–862, 2009.

23.

Vinh

N. X.

, Optimal Trajectory in Atmospheric Flight, Elsevier, New York, NY, USA, 1981.

24.

, “Entry guidance and trajectory control for reusable launch vehicle,” Journal of Guidance, Control, and Dynamics, vol. 20, no. 1, pp. 143–149, 1997.

25.

Messac

, “Physical programming: Effective optimization for computational design,” AIAA Journal, vol. 34, no. 1, pp. 149–158, 1996.

26.

Kovach

Cho

B. R.

, and Antony

, “Development of an experiment-based robust design paradigm for multiple quality characteristics using physical programming,” International Journal of Advanced Manufacturing Technology, vol. 35, no. 11–12, pp. 1100–1112, 2008.

27.

Sanchis

Martínez

M. A.

, and Blasco

, “Integrated multi-objective optimization and a priori preferences using genetic algorithms,” Information Sciences, vol. 178, no. 4, pp. 931–951, 2008.

28.

Ilgin

M. A.

and Gupta

S. M.

, “Physical programming: a review of the state of the art,” Studies in Informatics and Control, vol. 21, no. 4, pp. 349–366, 2012.

29.

Tian

Z. G.

Lin

D. M.

, and Wu

B. R.

, “Condition based maintenance optimization considering multiple objectives,” Journal of Intelligent Manufacturing, vol. 23, no. 2, pp. 333–340, 2013.

30.

Yong

Chen

, and Tang

, “Trajectory optimization of hypersonic gliding reentry vehicle based on the physical programming,” Acta Aeronautica et Astronautica Sinica, vol. 29, no. 5, pp. 1091–1097, 2008.

31.

Chen

Wan

Z. M.

, and Chen

S. L.

, “Parameterization method and constraints transformation method for RLV reentry trajectory optimization using genetic algorithm,” Journal of System Simulation, vol. 17, no. 11, pp. 2737–2740, 2005.

32.

Zhang

and Chen

, “Trajectory optimization for hypersonic vehicle satisfying maneuvering penetration,” Applied Mechanics and Materials, vol. 110–116, pp. 5223–5231, 2011.

33.

Zhang

Q. Z.

Wang

Z. B.

Tao

, and Sarker

B. R.

, “On-line optimization design of sliding mode guidance law with multiple constraints,” Applied Mathematical Modeling, vol. 37, no. 14–15, pp. 7568–7587, 2013.

34.

Zhang

Q. Z.

Wang

Z. B.

, and Tao

, “Optimal guidance law design for impact with terminal angle of attack constraint,” Optik, vol. 125, no. 1, pp. 243–251, 2014.

35.

Tao

Laili

Liu

, “Concept, principle and application of dynamic configuration for intelligent algorithms,” IEEE Systems Journal, vol. 8, no. 1, pp. 28–42, 2014.

36.

Dai

Tang

D. B.

Zheng

, and Cai

Q. X.

, “An improved genetic-simulated annealing algorithm based on a hormone modulation mechanism for a flexible flow-shop scheduling problem,” Advances in Mechanical Engineering, vol. 2013, Article ID 124903, 13 pages, 2013.

37.

Tao

Laili

, and Zhang

, “FC-PACO-RM: a parallel method for service composition optimal-selection in cloud manufacturing system,” IEEE Transactions on Industrial Informatics, vol. 9, no. 4, pp. 2023–2033, 2013.

38.

Kim

and Lee

, “Damage identification of trusses with elastic supports using FEM an genetic algorithm,” Advances in Mechanical Engineering, vol. 2013, Article ID 584901, 12 pages, 2013.

39.

Kennedy

and Eberhart

R. C.

, “Particle swarm optimization,” in Proceedings of the IEEE International Conference on Neural Networks, pp. 1942–1948, Perth, Australia, December 1995.

40.

Tao

Zhao

, and Zhou

, “Resource service composition and its optimal-selection based on particle swarm optimization in manufacturing grid system,” IEEE Transactions on Industrial Informatics, vol. 4, no. 4, pp. 315–327, 2008.

41.

Sahoo

N. C.

Ganguly

, and Das

, “Multi-objective planning of electrical distribution systems incorporating sectionalizing switches and tie-lines using particle swarm optimization,” Swarm and Evolutionary Computation, vol. 3, pp. 15–32, 2012.

42.

Talbi

and Belarbi

, “Fuzzy takagi sugeno system optimization using hybrid particle swarm optimization and tabu search learning algorithm,” International Journal of Tomography and Simulation, vol. 22, no. 1, pp. 4–16, 2013.

43.

Safari

Ardehali

M. M.

, and Sirizi

M. J.

, “Particle swarm optimization based fuzzy logic controller for autonomous green power energy system with hydrogen storage,” Energy Conversion and Management, vol. 65, pp. 41–49, 2013.

44.

van den Bergh

, An analysis of particle swarm optimizers [Ph.D. thesis], University of Pretoria, Pretoria, South Africa, 2001.

45.

Van den Bergh

and Engelbrecht

A. P.

, “A new locally convergent particle swarm optimiser,” in Proceedings of the IEEE International Conference on Systems, Man and Cybernetics, pp. 94–99, Hammamet, Tunisia, October 2002.

46.

Sun

Feng

, and Xu

, “Particle swarm optimization with particles having quantum behavior,” in Proceedings of the 2004 Congress on Evolutionary Computation (CEC '04) vol. 1, no. 4, pp. 325–331, San Diego, Calif, USA, June 2004.

47.

Sun

, and Feng

, “A global search strategy of quantum-behaved particle swarm optimization,” in Proceedings of the IEEE Conference on Cybernetics and Intelligent Systems, pp. 111–116, Singapore, December 2004.

48.

Pan

J. H.

Wang

, and Yang

X. G.

, “A random particle swarm optimization algorithm with application,” Advanced Materials Research, vol. 634, pp. 3940–6383944.

49.

Jamalipour

Sayareh

Gharib

Khoshahval

, and Karimi

M. R.

, “Quantum behaved particle swarm optimization with differential mutation operator applied to WWER-1000 in-core fuel management optimization,” Annals of Nuclear Energy, vol. 54, pp. 134–140, 2013.

50.

Jorris

T. R.

Schulz

C. S.

Friedl

F. R.

, and Rao

A. V.

, “Constrained trajectory optimization using pseudospectral methods,” in Proceedings of the Atmospheric Fright Mechanics Conference and Exhibit (AIAA '08), pp. 1–33, Honolulu, Hawaii, August 2008.

51.

Arora

R. K.

George

J. K.

, and Xavier

A. G.

, “Trajectory design for a flyback mission of a reusable launch vehicle,” in Proceedings of the Atmospheric Fright Mechanics Conference and Exhibit (AIAA '03), pp. 1–5, Austin, Texas, 2003.

52.

Wang

Y. L.

and Li

W. J.

, “Physical programming and its application in aircraft concept design,” Acta Aeronautica et Astronautica Sinica, vol. 26, no. 5, pp. 562–566, 2005.

53.

Zhang

Q. Z.

Liu

C. J.

Yang

, and Ren

, “Reentry trajectory planning optimization based on ant colony algorithm,” in Proceedings of the IEEE International Conference on Robotics and Biomimetics, pp. 1064–1068, Sanya, China, December 2007.

54.

Mooij

and Hänninen

P. G.

, “Distributed global trajectory optimization of a moderate lift-to-drag re-entry vehicle,” in Proceedings of the Guidance, Navigation, and Control Conference and Exhibit (AIAA '09), pp. 1–15, Chicago, Illinois, USA, August 2009.

55.

Chao

Wang

, and Yang

, “Maximum crossrange guidance under multiple constraints for lifting body reentry vehicle,” in Proceedings of the 3rd International Symposium on Systems and Control in Aeronautics and Astronautics (ISSCAA '10), pp. 1265–1270, Harbin, China, June 2010.

56.

Lin

W. Y.

Lee

W. Y.

, and Hong

T. P.

, “Adapting crossover and mutation rates in genetic algorithms,” Journal of Information Science and Engineering, vol. 19, no. 5, pp. 889–903, 2003.

57.

Tao

Zhang

Z. H.

, and Nee

A. Y. C.

, “A quantum multi-agent evolutionary algorithm for selection of partners in a virtual enterprise,” CIRP Annals-Manufacturing Technology, vol. 59, no. 1, pp. 485–488, 2010.

58.

Tao

Laili

Y. J.

Zhang

Z. H.

, and Nee

A. Y. C.

, “QMAEA: a quantum multi-agent evolutionary algorithm for multi-objective combinatorial optimization,” SIMULATION, vol. 90, no. 2, pp. 182–204.