Minimum-Fuel Ascent of Hypersonic Vehicle considering Control Constraint Using the Improved Pigeon-Inspired Optimization Algorithm

Trajectory optimization problem for hypersonic vehicles has long been recognized as a difficult problem. This paper brings control constraints into the trajectory optimization to make the optimal trajectory meet the requirements of control performance. The strong nonlinear characteristic of the ascent phase aerodynamics makes the trajectory optimization problem difficult to be solved by the optimal control theory. A trajectory optimization algorithm based on the improved pigeon-inspired optimization (PIO) algorithm is proposed to solve the complex trajectory optimization problem under multiple constraints. To overcome the obstacle of premature convergence and deceptiveness, the evolutionary strategy of qubit in quantum evolutionary algorithm (QEA) is introduced into the PIO to maintain population diversity and judge the optimal solution. To handle constraints, the penalty function is used to construct the fitness function. The optimal ascent trajectory is obtained by utilizing the improved PIO algorithm. Then, the trajectory inverse algorithm is used to verify the feasibility of the optimal trajectory to ensure that a feasible optimal trajectory is obtained. The comparison results show that the proposed algorithm outperforms particle swarm optimization (PSO) and standard PIO on trajectory optimization. Meanwhile, the simulation result shows that the performance of the optimal ascent trajectory with control constraints is improved and the trajectory is feasible. Therefore, the method is potentially feasible for solving the ascent trajectory optimization problem under control constraint for hypersonic vehicles.


Introduction
The hypersonic vehicle has received wide attention as it has high speed and large flight range. As an essential part, trajectory optimization for the hypersonic vehicle can help the design of the vehicle shape and a flight plan. The ascent trajectory optimization, due to its highly nonlinear dynamics and strong constraints, becomes a hot and difficult issue. Meanwhile, the reduction of fuel consumption in the ascent phase can help the hypersonic vehicle to increase the cruising distance.
A great deal of research has been done on the optimization of the hypersonic vehicle's trajectory, and a lot of new optimization methods have been proposed [1][2][3]. Kumar et al. presented a trajectory optimization of an aerodynamically controlled hypersonic boost-glide class of flight vehicles [4]. Xie et al. proposed a multiobjective gliding trajectory optimization scheme for a hypersonic vehicle with complicated constraints [5]. Chai et al. have done a lot of research on trajectory optimization and optimal control problem [6][7][8]; they considered the highly constrained trajectory optimization problems, applied a specific multiple-shooting discretization technique with the newest NSGA-III optimization algorithm, and constructed a new evolutionary optimal control solver to address the multiobjective trajectory planning problem [6]. Fu et al. solved the ascent trajectory optimization problem for hypersonic vehicles with the improved chicken swarm optimization (ICSO) algorithm [9]. Cheng et al. presented an iterative convex programming algorithm for the complex ascent trajectory planning problem, transcribed the dynamic equations into linearized algebraic equality constraints with a given initial guess based on the Newton-Kantorovich/Pseudospectral (N-K/PS) approach, and formulated the ascent trajectory planning problem as a convex programming problem [10]. Mahmoud et al. described the ascent and descent optimal trajectories with two different dynamics in each stage and the relation between the stages under the constraints of each phase [11].
However, the hypersonic vehicle is a strong coupling system; its plant performance, trajectory performance, and control performance are not independent but affect each other. Due to the limitation of control performance, the optimal trajectory that only considers trajectory performance may not be realized, so it is necessary to add control constraints to trajectory optimization.
At present, the numerical method is widely studied and applied in trajectory optimization [12][13][14], and it can be divided into indirect method and direct method. The indirect method uses optimal-control theory to obtain an analytical solution to the trajectory optimization problem. The merit of this method is that it can obtain an accurate analytical solution by the first-order necessary optimal condition; meanwhile, the initial costate is very difficult to guess. The direct method first transforms the trajectory optimization problem into nonlinear programming problem, then solves the nonlinear programming problem. Compared with the indirect method, the direct method has the advantages of low precision of the initial guess and simple derivation process. Recently, bioinspired heuristic algorithms have received widespread attention in the field of direct methods for solving the nonlinear programming problem. The advantage of the bioinspired heuristic algorithm is that the optimal solution is usually the global optimal solution, and these algorithms have strong robustness.
Due to the effectiveness and ease of implementation, more and more optimization algorithms based on swarm intelligence are proposed and researched nowadays [15][16][17]. The PIO algorithm performed better in a lot of popular benchmark optimization problems compared with DE [18] and PSO method and has proved effective in many areas. For example, three-dimensional path planning [19] controls parameter optimization [20] and predicts 3D protein structure [21]. Therefore, this paper chooses the PIO algorithm as the root-finding algorithm to solve the nonlinear programming problem.
Although the PIO algorithm is able to converge quickly, it is often inaccurate in dealing with complex optimization problems and easily falls into local optimization owing to premature convergence. Premature convergence mostly orig-inates from a loss of diversity and deceptiveness. The lack of diversity means that the difference between all solution candidates is small, which weakens the exploration. Meanwhile, a deceptive direction of convergence forestalls the exploration.
To overcome these obstacles and obtain better solutions, this paper introduces quantum representation and quantum rotation gate in the quantum evolutionary algorithm (QEA) to improve the PIO algorithm. The QEA is a probabilistic evolutionary algorithm that integrates concepts from quantum computing for robust search [22]. The QEA uses a qubit as the probabilistic representation, which represents a linear superposition of binary solutions [23]. The quantum rotation gate is the evolutionary strategy of qubit in the QEA [24], which is adopted as a mutation operator to update the pairs of probability amplitudes toward the one with the best fitness. In this paper, the position of pigeons is represented by qubit; the current best pigeon is considered to be a linear superposition of positive state and deceptive state. Other pigeon makes its own judgment by observation. If the current optimal solution still exists after iteration, the deceptive probability amplitude will decrease. In this way, the accuracy of the PIO algorithm is improved.
The rest of this paper is structured as follows. Section 2 briefly introduces the hypersonic vehicle model adopted in this paper; Section 3 establishes the optimization problem; constraints and performance index in trajectory optimization are determined; Section 4 proposes a new trajectory optimization based on the improved PIO algorithm and trajectory inverse algorithm in allusion to optimize ascent trajectory under complex constraints. Section 5 shows the simulation results. In this section, the comparison results of the proposed algorithm, PIO, and PSO on trajectory optimization are given. Meanwhile, the optimal trajectories with and without control constraints are compared; then, the control variables are obtained by using the trajectory inverse algorithm to determine whether the optimal trajectory is feasible. Finally, we conclude in the last section.

Vehicle Model
The hypersonic vehicle model used in this paper was developed by Bolender and Doman [25] and Parker et al. [26] as an attempt to extend earlier work done by Chavez and Schmidt [27]. The basic geometry of the vehicle model is shown in Figure 1, which contains only longitudinal dynamics Oblique shock Elevator Shear layer Figure 1: Geometry of the hypersonic vehicle.
2 International Journal of Aerospace Engineering (in the figure, α is the angle of attack; L v is vehicle length; δ e is elevator deflection; and M ∞ is freestream Mach number). A detailed explanation of the model derivation is given in [25]. Rather than using the relatively simple approach of Newtonian Impact Theory [27], the model was derived using a compressible flow theory. The combination of oblique shock wave and Prandtl-Meyer flow theory is used to determine the pressure within the range of possible angles of attack and structural bending conditions. The engine model is a scramjet taken directly from the paper by Chavez and Schmidt [27]. Since the angle of attack plays an important role in determining the characteristics of the flow characteristics into the scramjet, the thrust becomes very dependent on it.

Dynamic Equation.
The dynamic equations of motion for the longitudinal dynamics of the hypersonic vehicle are as follows: where v is the velocity; h is the altitude; α is the angle of attack; θ is the pitch angle; q p is the pitch rate; γ is the flight path angle; I y is the moment of inertia; g is the gravitational acceleration; R e is the earth radius; m is the mass of the vehicle; and T, L, D, and M y are thrust, lift, drag, and moment, respectively, expressed as L = 0:5ρv 2 sC L , D = 0:5ρv 2 sC D , where v is the velocity; s is the reference area; A e is the air inlet area; ρ is the atmospheric density; c is the mean aerodynamic chord; and C L , C D , C T , and C M are lift coefficient, drag coefficient, thrust coefficient, and moment coefficient, respectively. The coefficients determined by the angle of attack α and Mach number are as follows: where are determined by experimental and simulation data; α is the angle of attack; δ e is elevator deflection; and ϕ is fuel equivalent ratio. The exact analytical expressions of forces and moments are described in detail in [25].

Fuel Consumption.
In the ascent phase of the hypersonic vehicle, fuel mass is time-varying. Fuel consumption is calculated to determine the fuel mass at each moment. Fuel consumption per unit time of the vehicle is where I sp is the specific impulse of the engine; T is thrust; and g 0 is the normal acceleration of gravity. The fuel consumption between two adjacent points is where a is the acceleration of vehicle and Δv is the velocity difference between two adjacent points.

Optimization Problem Establishment
The hypersonic vehicle model is a complex dynamic system owing to the coupling between subsystems. Therefore, the trajectory performance and control performance are also coupled with each other during the ascent phase. The fully integrated system model can be more simply expressed as where xðtÞ is the vehicle states during the ascent phase and uðtÞ is the control input. In this way, the hypersonic vehicle ascent trajectory becomes the multidisciplinary system model given in Equation (6) [28]. Hence, the control performance constraint has to be taken into account when optimizing the ascent trajectory.
3.1. Optimization Problem Formulation. The optimal control problem can be simply described as seeking the optimal control variable and the corresponding state variable, so that the performance index achieves the extreme value under the constraints of equality and inequality. The optimization problem of the multidisciplinary system in Equation (6) is as follows: where t 0 and t f are the initial time and the terminal time, respectively; xðtÞ is the state variable; uðtÞ is input control variable; J = Φðxðt f Þ, t f Þ is performance index which is obtained by the objective function; φðxðt 0 Þ, xðt f Þ, t 0 , t f Þ = 3 International Journal of Aerospace Engineering 0 is given by boundary value constraint condition, which is converted to CðxðtÞ, uðtÞ, tÞ ≤ 0 by constraints during the flight.

Construction of Performance
Index. The mission in the ascent phase is to make the vehicle move to the cruise window safely. The reduction of fuel consumption can help the vehicle to increase the cruise distance. In this paper, the maximum residual mass at the end of the ascent phase is set as the optimal objective function.
where J opt is the performance index of the optimal trajectory in this paper; m f is the mass of the vehicle at the trajectory termination point, which is obtained by integrating Equation (4).
To introduce equality constraints and inequality constraints into the optimization problem, this paper uses penalty function, which is usually used to solve constrained optimization problems, to construct objective function. The basic idea is that, according to the characteristics of constraint conditions, it can be converted into some punishment function and added to the objective function, to transform the constrained optimization problem into an unconstrained optimization problem to be solved, and the optimal solution can be obtained by solving the unconstrained optimization problem. The penalty function PðxÞ constructed in this paper is as follows: where hðxÞ and gðxÞ represent equality constraints and inequality constraints, respectively, and x is the state variable. Each term in the penalty function is multiplied by a certain penalty parameter and then combined with the original objective function to obtain the augmented objective function.
where σ is the penalty matrix and x is the state variable.

Constraints for Ascent
Phase. In this part, various constraints included in hðxÞ and gðxÞ are given. In addition to the common constraint, such as path constraints and terminal constraints, the impact of control constraints and state constraints on trajectory will also be taken into consideration.

Path Constraints.
During the flight, intense friction with the atmosphere will produce high temperatures. Therefore, the heat flux constraint must be considered to prevent the excessive surface temperature of the vehicle. The heat flux constraint is as follows: where C is constant; ρ is air density; and v is the velocity. Dynamic pressure generated by high-speed flight provides aerodynamic force and control torque to adjust attitude. However, if the dynamic pressure exceeds the limit, it will have a great impact on the vehicle. The dynamic pressure constraint is as follows: where ρ is air density and v is the velocity. In order to ensure the structural stability of the vehicle, the overload constraint should be considered in the ascent trajectory optimization, and the overload constraint is as follows: where g is the gravitational acceleration; m is the mass of the vehicle; and L and D are lift and drag, respectively.

Terminal Constraints.
The terminal constraint is related to flight missions. In this paper, the flight height, the velocity, and the flight path angle are required to meet specific constraint as follows: where v t , h t , and γ t are actual cruising state; v f , h f , and γ f are the predetermined cruising state; and Δv, Δh, and Δγ are acceptable error.

Control Constraints and State Constraints.
Although the ascent trajectory without considering control constraints and state constraints is theoretically feasible, in practical application, due to the influence of ascending requirement and manoeuvring performance, it is necessary to control the state variables within a certain range and consider the control feasibility of the vehicle. As mentioned before, the control performance is affecting the ascent trajectory. Consequently, control constraints and state constraints are needed in trajectory optimization. In this paper, for the state constraint, the constraints of the flight altitude, the velocity, and the flight path angle are granted.
where v min and v max are the lower and upper limits of the velocity; h min and h max are the lower and upper limits of the flight altitude; and γ min and γ max are the lower and upper limits of the flight path angle. 4 International Journal of Aerospace Engineering For the control constraint, the constraints of the angle of attack and the elevator deflection are granted.
where α min and α max are the lower and upper limits of the angle of attack and δ emin and δ emax are the lower and upper limits of elevator deflection. Moreover, considering the attitude stability of the vehicle, the rate of change of the flight path angle cannot be too high or too low during the ascent phase. In this paper, the absolute value of the rate of change of the flight path angle never exceeds 0.002 rad/s.

Trajectory Optimization Algorithm
As mentioned in the previous section, trajectory performance and control performance are both considered in trajectory optimization. As the constraints increase, the complexity of trajectory optimization increases. An effective trajectory optimization algorithm is needed to solve the trajectory optimization problem under complex constraints.
The optimization can be divided into two steps: (1) identifying feasible designs and (2) identifying the optimal design from the set of feasible candidates [29]. This paper proposes an ascent trajectory optimization algorithm to solve a complex optimization problem under multiple constraints, as shown in Figure 2. The algorithm first establishes the constraints and fitness function of the ascent trajectory, optimizes the ascent trajectory using the optimization algorithm, and then performs a feasibility analysis on the optimal trajectory. When the trajectory is not feasible, the penalty value of fitness function is adjusted until the feasible optimal ascent trajectory is obtained. The specific methods included in the algorithm are introduced, respectively, in this section.
4.1. The Improved PIO Algorithm. In the global optimization problem, especially for multimodal problems, when the PIO algorithm cannot generate new children, it is easy for the algorithm to converge to local optimization prematurely. The main reason for premature convergence is the low diversity of the algorithm. Also, in the convergence stage, when the population aggregates at a certain point, even if this is not the correct global optimal solution, it will be regarded as the global optimal solution. Therefore, to overcome these weaknesses in the PIO algorithm, this paper represents the position of pigeons by qubit, the current best pigeon is considered to be a linear superposition of positive state and deceptive state. The other pigeon makes its own judgment by observation. Meanwhile, the quantum rotation gate is adopted in the PIO algorithm as a mutation operator to enhance positive probability.
where N P is the population size.
For the sake of the optimal solution to meet the multiple constraints presented in the previous section, strong constraints need to be applied to the value range of each possible solution element; i.e., the position and speed constraints of the pigeon are as follows: where i = 1, 2, ⋯N P ; X low and X up are the lower and upper limits of each element in the position, respectively; and V low and V up are the lower and upper limits of each element in the speed, respectively. In the map compass operator, each individual in the population updates its speed and position through geomagnetic, solar height information and the optimal information in the population.
where t is the current number of iterations, R is the map compass factor, X g is the global optimal solution in the current population, and rand is the uniform random value between the [0,1] regions. According to the map compass operator, each pigeon in the population adjusts its direction according to Equation (19). The map compass operator of the PIO algorithm is improved in this paper, so as to solve the problem that the basic PIO algorithm falls into local optimal due to premature convergence.
The improved PIO algorithm adopts qubit to describe the current state of the pigeon. To maintain population diversity, the position state of the pigeon is obtained according to Monte Carlo random simulation: where 5 International Journal of Aerospace Engineering where t is the number of iterations; u and f are uniformly distributed random numbers on [0,1]; P i ðtÞ is the historical optimal position at tth iteration; P g ðtÞ is the global optimal position at the tth iteration; and L is defined as follows: where ωðtÞ is inertia weight which has great influence on the convergence of the algorithm; m best ðtÞ is the average optimal position of all pigeons in the population at tth iteration and ωðtÞ and m best ðtÞ are defined as follows: where N P is the population size; ω max and ω min are lower and upper limits of inertia weight; and T is the maximum number of iterations.
where N P is the population size; t is the number of iterations; X c is the position of centre point; and fitness () is a fitness function, which is constructed by the objective function in Equation (10), i.e.,

Judgment on Global Optimal Solution.
To avoid the problem that the wrong global optimal solution was obtained by reason of the aggregation of pigeons, herein, this paper adopts the quantum rotating gate in the QEA [24] to solve this problem.
Qubit chromosomes can be represented as follows: where q is qubit chromosomes and jα i j 2 + jβ i j 2 = 1, i = 1, ⋯, m, m, is the length of the chromosome. jα i j 2 and jβ i j 2 represent the probability that the state of the quantum bit is 0 and 1, respectively. By comparing the random number generated within the interval with jα i j 2 , if the random number is larger than jα i j 2 , the corresponding qubit of the chromosome is set as 1, otherwise set as 0 [30]. In the QEA, the probability amplitude update is calculated by where Δ + θ is the rotation angle and t is the number of iterations. Unlike the quantum bit evolution strategy in QEA, it is used as a mutation operator to enhance positive probability. The improved PIO algorithm will get a global optimal x best ðtÞ after each iteration, but for the reasons mentioned above, this global optimal may not be the accurate global optimal. Probability amplitudes are initialized as α i = β i = ffiffi ffi 2 p /2, i = 1, 2, ⋯, n. If the current global optimal solution remains after iteration, run the mutation operator to increase α i ; this means that x best ðtÞ is more likely to be a global optimal International Journal of Aerospace Engineering solution. Otherwise, the probability amplitude is reset to initial values to maintain vigilance against the deceptiveness.

The
Procedure of the Improved PIO Algorithm. The procedure of the improved PIO algorithm is as follows.
Step 1. Initialize population information and parameters of the improved PIO algorithm. Set a random speed and path for each pigeon.
Step 2. Compare the fitness of each pigeon, set x best (t) with the current best fitness.
Step 3. Use the improved map compass operator to update; then use the landmark operator to convergence; obtain g best (t).
Step 4. Compare g best (t) with x best (t). If they are the same, go to (a); otherwise, go to (b).

11
International Journal of Aerospace Engineering (a) Conduct the mutation operator by Equation (29) (b) Conduct the mutation operator by Equation (29) Step 5. If the maximum iteration number is reached, get global optimum x best (t). Otherwise, return to Step 3.
The above steps can be summarized in Figure 3. To validate the improved PIO algorithm, three test functions are used to compare three algorithms, i.e., PSO [31], PIO [18], and the improved PIO. In order to rule out the effect of a few special values on the results, all experimental data in this section are the average of the values obtained from 40 independent runs of each algorithm. The optimized results for the test functions of these algorithms are listed in Table 1, where D is the dimension of the design variables. To avoid the difference in optimized results originating from the selection of control parameters, all equivalent control parameters are set to be the same (the specific parameters are listed in Table 2).
The improved PIO algorithm shows great potential for solving the multimodal optimization problem in terms of accuracy. The optimal value of the improved PIO algorithm is much smaller than that of other algorithms; the improved PIO algorithm also has the shortest running time. Therefore, the improved PIO algorithm outperforms PSO and PIO on benchmark functions.

Trajectory Feasibility Analysis.
To verify whether the control variables required by the trajectory conforms to constraints, this paper uses a trajectory inverse algorithm, which can solve the control variables needed according to the state variables at the current time node and the state variables at the next node.
For the vehicle, the known output is an ascent trajectory varying with time, and the algorithm can work out the value of the actual control variable required by the vehicle following a given trajectory. By analyzing the change curve of the control variable, it can be estimated whether the required control variable exceeds the physical limit when the vehicle follows the reference trajectory [32].
In the trajectory inverse algorithm, the dynamic system to be analyzed can be defined by state variable x, control variable u, and an output vector y. The motion equation is shown as follows: where f and g represent any functions related to system variables; x 0 is the initial value of the state variable; and y is usually used to represent the known trajectory. The above

12
International Journal of Aerospace Engineering  13 International Journal of Aerospace Engineering formula can associate the control variable u with the known trajectory to solve the required control variable.
As shown in Equation (31), the state differential value of the trajectory at the discrete point n − 1 is determined by the state variable x n−1 at the current point and the current control input variable u n−1,k . The state variable differential of dis-crete points can obtain the state variable of the next discrete point through the integration of time.
where t n is the time at node n; t n−1 is the time at node n − 1; and k is the number of iterations. Let Δt be the time interval between two nodes, i.e., When Δt is fairly small, Equation (32) can be approximated as    The error vector e n,k is defined as the error between the state variable differential estimation value at the nth discrete point and the target trajectory point _ x n after the trajectory is solved in the kth iteration, as shown in Equation (35). According to the defined error vector, the control variable of the next iteration can be obtained by Newton's method, as shown in Equation (36).
where J is the jacobian matrix, which is defined as follows: where u α , u δ e , and u ϕ are input of the angle of attack, the elevator deflection, and the fuel equivalent ratio, respectively; and e _ x v , e _ x γ , and e _ The trajectory inverse algorithm needs to consider the mass change caused by vehicle fuel consumption. The mass of the vehicle needs to be calculated according to the vehicle mass during the last iteration and the instantaneous fuel flow rate. Assuming that the fuel consumption rate remains unchanged within the time interval Δt, the following equation can be obtained: where m n and m n−1 are the mass of the vehicle at the n discrete point and the n − 1 discrete point, respectively, and _ m is fuel consumption per unit time of the vehicle given in Equation (4).
The entire procedure of trajectory inverse algorithm is shown in Figure 4.

Results
The results in this section are mainly divided into two parts. The first part is the comparison of the optimal trajectories obtained by different methods. The second part is to verify the feasibility of the optimal trajectory obtained by the improved PIO algorithm using the trajectory inverse algorithm.
This section uses the Doman model described in Section 2 and studies its optimal ascent trajectory from 6 Mach, altitude 22379 m, to 7.8 Mach, altitude 26495 m. The optimization variables were obtained by using the Chebyshev    16 International Journal of Aerospace Engineering interpolation method to discretize the angle of attack and elevator deflection. Then, the optimization variables are expressed as follows: where α is the angle of attack; δ e is the elevator deflection; X is the position vector of pigeons; D is the dimension of the position vector; and t f is the terminal time; 5.1. Optimal Ascent Trajectory. To invest the feasibility and optimality of the proposed method, comparative studies with other global optimization methods such as PSO and PIO are presented. Besides, to analyze the impact of control constraints, the optimal trajectories with and without control constraints are obtained by the proposed method. The control parameters of global optimization methods are given in Table 2; to avoid the influence of parameters on the effectiveness of the optimization algorithm, all equivalent control parameters are set to be the same. According to the actual situation, select the initial state as ½v, h, r, γ, m = ½1780 m/s, 22379 m, 0, 5°, 16817 kg. For constraints, set the terminal constraints as v f = 2337 m/s, h f = 26495 m, and γ f = 0; the path constraints are _ Q max = 600 kW/m 2 , 10 kPa ≤ q ≤ 100 kPa, and n max = 3:5; the state constraints are 1500 m/s ≤ v ≤ 2500 m/s and 20000 m ≤ h ≤ 30000 m. The above constraints are the basic constraints, and this paper will optimize a standard ascent trajectory based on the basic constraints. This paper adds control constraints to the basic constraints so that the control performance limitation can be considered into trajectory optimization; the control constraints are as follows: where α is the angle of attack; δ e is the elevator deflection; ϕ is the fuel equivalent ratio; and γ is the flight path angle. Besides, if the flight path angle too large or too small, the difficulty of trajectory control will increase. Therefore, this paper also constrains the flight path angle, i.e., −2°≤ γ ≤ 5°. The optimal trajectories generated using different algorithms are plotted in Figures 5 and 6. Table 3 lists the terminal states of optimal trajectories obtained by different algorithms.
When control constraint is considered, it can be seen from Table 3 and Figure 5 that the optimal trajectory obtained by the proposed algorithm is most in line with terminal constraints and has the lowest fuel consumption. Figure 6 shows that compared to the trajectories obtained by other algorithms, the dynamic pressure of the trajectory obtained by the proposed algorithm is largest in the constrained range and the higher dynamic pressure resulted in better performance. The heat flux of the trajectory obtained by the proposed algorithm is smallest, which can reduce the requirement of aerodynamic heating for the vehicle materials

18
International Journal of Aerospace Engineering to a certain extent. The overload amplitude of the trajectory obtained by the proposed algorithm decreases, which ensures better control performance. Therefore, the proposed trajectory optimization algorithm based on the improved PIO algorithm is feasible and optimal.
The optimal trajectories obtained by the proposed algorithm with and without control constraints in Figures 5  and 6 and Table 3 show that with control constraints having been taken into account, the fuel consumption increases by 0.11% while the range increases by 4.50%, which is of great

19
International Journal of Aerospace Engineering practical significance for the long-distance flight. Besides, the ascent trajectory with control constraints avoids the subduction segment and the flight path angle changes more gently; that is, the trajectory control performance is better. From the perspective of dynamic pressure, heat flux, and overload, the trajectory performance with control constraints is also improved.
In order to further analyze the computational complexity of the proposed method, Table 4 illustrates the results of the three methods in terms of the computational time and values of the fitness function, while Figure 7 shows the fitness value versus iteration time. Table 4 shows that the proposed trajectory optimization algorithm has much less computational time and the minimum fitness value. Besides, it can be seen from Figure 7 that the proposed trajectory optimization algorithm can converge to the global optimal solution fastest. These results further prove that the proposed algorithm is feasible and optimal.

Trajectory Feasibility.
The optimal trajectory is input into the trajectory inverse algorithm to verify the feasibility of trajectory with control constraints. The simulation results are shown in Figure 8.
The control variables obtained by the trajectory inverse algorithm are all within the constraint range. The control variables are brought into the dynamic equations of the vehicle, and the actual ascent trajectory is obtained by numerical integration of state variables by the Runge-Kutta method and compared with the optimal trajectory solved by the proposed algorithm; the results are shown in Figure 9.
The actual trajectory obtained by the trajectory inverse algorithm is consistent with the optimal trajectory solved by the proposed algorithm, thus verifying the feasibility of the optimal trajectory, and the effectiveness of the algorithm is also proved.

Conclusion
In this paper, a new trajectory optimization algorithm based on the improved PIO algorithm is proposed to solve the trajectory optimization problem under control constraints.
To get rid of premature convergence, a quantum-based approach is incorporated into the PIO algorithm to perceive deceptiveness and preserve the diversity of the population. To be specific, the position of pigeons is represented by qubit, the quantum representation of the current best solution can effectively maintain the diversity in exploration. Then, the quantum rotation gate is used to decrease the deceptive probability.
By comparing the optimization result with trajectory optimization algorithm based on PSO and the standard PIO, it can be seen that although the trajectories are quite different, the trajectory obtained by the proposed algorithm has a better performance index, and the terminal states are more accurate. Meanwhile, compared with the optimal trajectory without control constraints, the optimal trajectory with control constraints has better control performance. Also, the proposed algorithm has shorter computational time and smaller fitness value and can converge to the optimal solution faster. Therefore, the proposed algorithm is feasible in trajectory optimization.
In the following research, two problems can be further studied. The first point is how to choose the appropriate interpolation method for the given optimal control problem. Different interpolation methods and interpolation nodes affect the effectiveness and efficiency of the algorithm. The second point is that the efficiency of the PIO algorithm is affected by parameters. How to select appropriate parameters to optimize the efficiency of the algorithm is worth studying.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.