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

Authors: Jiahao He; Yanbin Liu; Shuanglin Li; Yue Tang
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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–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–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–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–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 originates 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. Section2 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.


## 2. 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 (in the figure, α is the angle of attack; Lv 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.Figure 1
Geometry of the hypersonic vehicle.


### 2.1. Dynamic Equation

The dynamic equations of motion for the longitudinal dynamics of the hypersonic vehicle are as follows:(1)v̇=Tcosα−Dm−gsinθ−α,α̇=qp−Tsinα+Lmv+gv−vRe+hcosγ,ḣ=vsinθ−α,θ̇=qp,q̇p=MyIy,where v is the velocity; h is the altitude; α is the angle of attack; θ is the pitch angle; qp is the pitch rate; γ is the flight path angle; Iy is the moment of inertia; g is the gravitational acceleration; Re is the earth radius; m is the mass of the vehicle; and T, L, D, and My are thrust, lift, drag, and moment, respectively, expressed as
(2)L=0.5ρv2sCL,D=0.5ρv2sCD,T=0.5ρv2AeCT,My=0.5ρv2sc¯CM,where v is the velocity; s is the reference area; Ae is the air inlet area; ρ is the atmospheric density; c¯ is the mean aerodynamic chord; and CL, CD, CT, and CM 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:
(3)CL=CL0+CLα×α+CLδe×δe,CD=CD0+CDα×α+CDδe×δeCT=CT0+CTα×α+CTϕ×ϕCM=CM0+CMα×α+CMδe×δewhere CL0, CLα, CLδe, CD0, CDα,CDδe, CT0, CTα,CTϕ,CM0,CMα, andCMδeare 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].


### 2.2. 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(4)ṁ=−Tg0⋅Isp,where Isp is the specific impulse of the engine; T is thrust; and g0 is the normal acceleration of gravity. The fuel consumption between two adjacent points is
(5)dm=ṁ⋅△va,where a is the acceleration of vehicle and ∆vis the velocity difference between two adjacent points.


## 2.1. Dynamic Equation

The dynamic equations of motion for the longitudinal dynamics of the hypersonic vehicle are as follows:(1)v̇=Tcosα−Dm−gsinθ−α,α̇=qp−Tsinα+Lmv+gv−vRe+hcosγ,ḣ=vsinθ−α,θ̇=qp,q̇p=MyIy,where v is the velocity; h is the altitude; α is the angle of attack; θ is the pitch angle; qp is the pitch rate; γ is the flight path angle; Iy is the moment of inertia; g is the gravitational acceleration; Re is the earth radius; m is the mass of the vehicle; and T, L, D, and My are thrust, lift, drag, and moment, respectively, expressed as
(2)L=0.5ρv2sCL,D=0.5ρv2sCD,T=0.5ρv2AeCT,My=0.5ρv2sc¯CM,where v is the velocity; s is the reference area; Ae is the air inlet area; ρ is the atmospheric density; c¯ is the mean aerodynamic chord; and CL, CD, CT, and CM 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:
(3)CL=CL0+CLα×α+CLδe×δe,CD=CD0+CDα×α+CDδe×δeCT=CT0+CTα×α+CTϕ×ϕCM=CM0+CMα×α+CMδe×δewhere CL0, CLα, CLδe, CD0, CDα,CDδe, CT0, CTα,CTϕ,CM0,CMα, andCMδeare 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].


## 2.2. 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(4)ṁ=−Tg0⋅Isp,where Isp is the specific impulse of the engine; T is thrust; and g0 is the normal acceleration of gravity. The fuel consumption between two adjacent points is
(5)dm=ṁ⋅△va,where a is the acceleration of vehicle and ∆vis the velocity difference between two adjacent points.


## 3. 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(6)ẋt=fxt,ut,t,where xt is the vehicle states during the ascent phase and ut 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:
(7)minJ=Φxtf,tf,s.t.ẋt=Fxt,ut,t,Cxt,ut,t≤0,φxt0,xtf,t0,tf=0,where t0 and tf are the initial time and the terminal time, respectively; xt is the state variable; ut is input control variable; J=Φxtf,tf is performance index which is obtained by the objective function; φxt0,xtf,t0,tf=0 is given by boundary value constraint condition, which is converted to Cxt,ut,t≤0 by constraints during the flight.


### 3.2. 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.(8)minJopt=−mf,where Jopt is the performance index of the optimal trajectory in this paper; mf 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 functionP¯x constructed in this paper is as follows:
(9)P¯x=∑hi2x+∑min0,gix2,where hx and gx 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.(10)minJopt=−mf+σP¯x,where σ is the penalty matrix and x is the state variable.


### 3.3. Constraints for Ascent Phase

In this part, various constraints included inhx and gx 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.


#### 3.3.1. 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:(11)Q̇=Cρv3.15≤Q̇max,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:
(12)qmin≤q=12ρv2≤qmax,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:
(13)n=L2+D2mg≤nmax,where g is the gravitational acceleration; m is the mass of the vehicle; and L and D are lift and drag, respectively.


#### 3.3.2. 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:(14)vt=vf+△v,ht=hf+△h,γt=γf+△γ,where vt, ht, and γt are actual cruising state; vf, hf, and γf are the predetermined cruising state; and Δv, Δh, and Δγ are acceptable error.


#### 3.3.3. 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.(15)vmin≤v≤vmax,hmin≤h≤hmax,γmin≤γ≤γmax,where vmin and vmax are the lower and upper limits of the velocity; hmin and hmax 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.For the control constraint, the constraints of the angle of attack and the elevator deflection are granted.(16)αmin≤α≤αmax,δemin≤δe≤δemax,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.


## 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:
(7)minJ=Φxtf,tf,s.t.ẋt=Fxt,ut,t,Cxt,ut,t≤0,φxt0,xtf,t0,tf=0,where t0 and tf are the initial time and the terminal time, respectively; xt is the state variable; ut is input control variable; J=Φxtf,tf is performance index which is obtained by the objective function; φxt0,xtf,t0,tf=0 is given by boundary value constraint condition, which is converted to Cxt,ut,t≤0 by constraints during the flight.


## 3.2. 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.(8)minJopt=−mf,where Jopt is the performance index of the optimal trajectory in this paper; mf 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 functionP¯x constructed in this paper is as follows:
(9)P¯x=∑hi2x+∑min0,gix2,where hx and gx 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.(10)minJopt=−mf+σP¯x,where σ is the penalty matrix and x is the state variable.


## 3.3. Constraints for Ascent Phase

In this part, various constraints included inhx and gx 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.


### 3.3.1. 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:(11)Q̇=Cρv3.15≤Q̇max,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:
(12)qmin≤q=12ρv2≤qmax,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:
(13)n=L2+D2mg≤nmax,where g is the gravitational acceleration; m is the mass of the vehicle; and L and D are lift and drag, respectively.


### 3.3.2. 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:(14)vt=vf+△v,ht=hf+△h,γt=γf+△γ,where vt, ht, and γt are actual cruising state; vf, hf, and γf are the predetermined cruising state; and Δv, Δh, and Δγ are acceptable error.


### 3.3.3. 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.(15)vmin≤v≤vmax,hmin≤h≤hmax,γmin≤γ≤γmax,where vmin and vmax are the lower and upper limits of the velocity; hmin and hmax 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.For the control constraint, the constraints of the angle of attack and the elevator deflection are granted.(16)αmin≤α≤αmax,δemin≤δe≤δemax,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.


## 3.3.1. 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:(11)Q̇=Cρv3.15≤Q̇max,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:
(12)qmin≤q=12ρv2≤qmax,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:
(13)n=L2+D2mg≤nmax,where g is the gravitational acceleration; m is the mass of the vehicle; and L and D are lift and drag, respectively.


## 3.3.2. 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:(14)vt=vf+△v,ht=hf+△h,γt=γf+△γ,where vt, ht, and γt are actual cruising state; vf, hf, and γf are the predetermined cruising state; and Δv, Δh, and Δγ are acceptable error.


## 3.3.3. 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.(15)vmin≤v≤vmax,hmin≤h≤hmax,γmin≤γ≤γmax,where vmin and vmax are the lower and upper limits of the velocity; hmin and hmax 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.For the control constraint, the constraints of the angle of attack and the elevator deflection are granted.(16)αmin≤α≤αmax,δemin≤δe≤δemax,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.


## 4. 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.Figure 2
Ascent trajectory optimization algorithm.


### 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.


#### 4.1.1. Map/Compass Operator

In the basic PIO algorithm, the spatial position of the pigeon is the possible solution to the optimization problem. The position of the pigeon’s home is the global optimal solution. The fitness of the pigeon is the value of the objective function. The positionXi and the speed Vi of the ith pigeon are
(17)XiT=Xi,1,Xi,2,⋯,Xi,n,i=1,2,⋯,NP,ViT=Vi,1,Vi,2,⋯,Vi,n,i=1,2,⋯,NP,where NP 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:(18)Xlow≤Xii≤Xup,Vlow≤Vii≤Vup,where i=1,2,⋯NP; Xlow and Xup are the lower and upper limits of each element in the position, respectively; and Vlow and Vup 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.(19)Vit=Vit−1⋅e−Rt+rand⋅Xg−Xit−1,(20)Xit=Xit−1+Vit,where t is the current number of iterations, R is the map compass factor, Xg 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:(21)Xit=PiPgt±L2ln1u,where
(22)PiPgt=ft×Pit+1−ft×Pgt,where t is the number of iterations; u and f are uniformly distributed random numbers on [0,1]; Pit is the historical optimal position at tth iteration; Pgt is the global optimal position at the tth iteration; and L is defined as follows:
(23)L=2ωtmbestt−Xit,where ωt is inertia weight which has great influence on the convergence of the algorithm; mbestt is the average optimal position of all pigeons in the population at tth iteration and ωt and mbestt are defined as follows:
(24)ωt=ωmax−ωmax−ωmin×tT,mbestt=1NP∑i=1NPPit,where NP is the population size; ωmax and ωmin are lower and upper limits of inertia weight; and T is the maximum number of iterations.Hence, the improved PIO algorithm updates pigeons’ position according to the following principle:(25)Xit+1=PiPgt+ωt×mbestt−Xit×ln1ft,ft≥0.5,PiPgt−ωt×mbestt−Xit×ln1ft,ft<0.5.


#### 4.1.2. Landmark Operator

In the landmark operator, some individuals with better fitness are known as pigeons familiar with the landmarks, and the remaining pigeons may follow the elite pigeons or be abandoned by the population. In each iteration, pigeons are sorted according to the fitness, and half of the pigeons with poor fitness are eliminated. Then, the population centre pointCt is selected as the reference direction of the remaining pigeons and the position of the pigeon is updated.
(26)NPt=NPt−12,Xct=∑Xit⋅fitnessXitNP∑fitnessXit,Xit=Xit−1+rand⋅Xc−Xit−1,where NP is the population size; t is the number of iterations; Xc is the position of centre point; and fitness () is a fitness function, which is constructed by the objective function in Equation (10), i.e.,
(27)fitnessx=Joptmin=−mf+σP¯x.


#### 4.1.3. 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:(28)q=α1⋯αi⋯αmβ1⋯βi⋯βm,where q is qubit chromosomes and αi2+βi2=1, i=1,⋯,m,m, is the length of the chromosome. αi2 and βi2 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 αi2, if the random number is larger than αi2, 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(29)αit+1=cosΔθ−sinΔθαit1−αit2,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 optimalxbestt 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=2/2,i=1,2,⋯,n. If the current global optimal solution remains after iteration, run the mutation operator to increase αi; this means that xbestt is more likely to be a global optimal solution. Otherwise, the probability amplitude is reset to initial values to maintain vigilance against the deceptiveness.


#### 4.1.4. 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, setxbest(t) with the current best fitness.Step 3. Use the improved map compass operator to update; then use the landmark operator to convergence; obtaingbest(t).Step 4. Comparegbest(t) with xbest(t). If they are the same, go to (a); otherwise, go to (b).
(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 optimumxbest(t). Otherwise, return to Step 3.The above steps can be summarized in Figure3.Figure 3
The procedure of the improved PIO algorithm.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).Table 1
Test function optimization results.
Test functionDAlgorithmOptimal valueTime (s)Ackley5PSO7.1431E-040.11PIO0.01370.12The improved PIO5.0849E-150.09Rastrigin5PSO5.7730E-060.13PIO0.01210.12The improved PIO00.09Rosenbrock2PSO1.0878E-090.11PIO1.0900E-080.12The improved PIO7.9993E-330.09Table 2
Parameters of global optimization methods.
ParametersPSOPIOThe improved PIOPopulation size (NP)100100100Maximum number of iterations (T)505050Map factor (R)N/A0.20.2Learning factor (c1,c2)[2,2]N/AN/ADimension (D)101010Rotating angle (°)N/AN/A-11The 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.


### 4.2. 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 variablex, control variable u, and an output vector y. The motion equation is shown as follows:
(30)ẋ=fx,u,x0=x0,y=gx,where f and g represent any functions related to system variables; x0 is the initial value of the state variable; and y is usually used to represent the known trajectory. The above 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 xn−1 at the current point and the current control input variable un−1,k. The state variable differential of discrete points can obtain the state variable of the next discrete point through the integration of time.
(31)ẋn−1,k=fxn−1,un−1,k,(32)xn,k=∫tn−1tnẋtdt+xn−1,k,where tn is the time at node n; tn−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.,
(33)Δt=tn−tn−1.WhenΔtis fairly small, Equation (32) can be approximated as
(34)xn,k≈ẋn−1,k△t+xn−1,k.The error vectoren,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).
(35)en,k=ẋn,k−x¯̇n,(36)un−1,k+1=un−1,k−J−1en,k,where J is the jacobian matrix, which is defined as follows:
(37)J=∂eẋv∂uα∂eẋv∂uδe∂eẋv∂uϕ∂eẋγ∂uα∂eẋγ∂uδe∂eẋγ∂uϕ∂eẋqp∂uα∂eẋqp∂uδe∂eẋqp∂uϕ,where uα,uδe, and uϕ are input of the angle of attack, the elevator deflection, and the fuel equivalent ratio, respectively; and eẋv,eẋγ, and eẋqp are the errors of the estimated values of the velocity, the path angle, and the pitch rate, respectively.Through the iterative solution of Equation (35), the control variable at the n−1 discrete point can be obtained, so that the trajectory can change from the n−1 discrete point to the point n.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:
(38)mn=mn−1+ṁΔt,where mn and mn−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 Figure4.Figure 4
The procedure of the trajectory inverse algorithm.


## 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.


### 4.1.1. Map/Compass Operator

In the basic PIO algorithm, the spatial position of the pigeon is the possible solution to the optimization problem. The position of the pigeon’s home is the global optimal solution. The fitness of the pigeon is the value of the objective function. The positionXi and the speed Vi of the ith pigeon are
(17)XiT=Xi,1,Xi,2,⋯,Xi,n,i=1,2,⋯,NP,ViT=Vi,1,Vi,2,⋯,Vi,n,i=1,2,⋯,NP,where NP 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:(18)Xlow≤Xii≤Xup,Vlow≤Vii≤Vup,where i=1,2,⋯NP; Xlow and Xup are the lower and upper limits of each element in the position, respectively; and Vlow and Vup 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.(19)Vit=Vit−1⋅e−Rt+rand⋅Xg−Xit−1,(20)Xit=Xit−1+Vit,where t is the current number of iterations, R is the map compass factor, Xg 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:(21)Xit=PiPgt±L2ln1u,where
(22)PiPgt=ft×Pit+1−ft×Pgt,where t is the number of iterations; u and f are uniformly distributed random numbers on [0,1]; Pit is the historical optimal position at tth iteration; Pgt is the global optimal position at the tth iteration; and L is defined as follows:
(23)L=2ωtmbestt−Xit,where ωt is inertia weight which has great influence on the convergence of the algorithm; mbestt is the average optimal position of all pigeons in the population at tth iteration and ωt and mbestt are defined as follows:
(24)ωt=ωmax−ωmax−ωmin×tT,mbestt=1NP∑i=1NPPit,where NP is the population size; ωmax and ωmin are lower and upper limits of inertia weight; and T is the maximum number of iterations.Hence, the improved PIO algorithm updates pigeons’ position according to the following principle:(25)Xit+1=PiPgt+ωt×mbestt−Xit×ln1ft,ft≥0.5,PiPgt−ωt×mbestt−Xit×ln1ft,ft<0.5.


### 4.1.2. Landmark Operator

In the landmark operator, some individuals with better fitness are known as pigeons familiar with the landmarks, and the remaining pigeons may follow the elite pigeons or be abandoned by the population. In each iteration, pigeons are sorted according to the fitness, and half of the pigeons with poor fitness are eliminated. Then, the population centre pointCt is selected as the reference direction of the remaining pigeons and the position of the pigeon is updated.
(26)NPt=NPt−12,Xct=∑Xit⋅fitnessXitNP∑fitnessXit,Xit=Xit−1+rand⋅Xc−Xit−1,where NP is the population size; t is the number of iterations; Xc is the position of centre point; and fitness () is a fitness function, which is constructed by the objective function in Equation (10), i.e.,
(27)fitnessx=Joptmin=−mf+σP¯x.


### 4.1.3. 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:(28)q=α1⋯αi⋯αmβ1⋯βi⋯βm,where q is qubit chromosomes and αi2+βi2=1, i=1,⋯,m,m, is the length of the chromosome. αi2 and βi2 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 αi2, if the random number is larger than αi2, 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(29)αit+1=cosΔθ−sinΔθαit1−αit2,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 optimalxbestt 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=2/2,i=1,2,⋯,n. If the current global optimal solution remains after iteration, run the mutation operator to increase αi; this means that xbestt is more likely to be a global optimal solution. Otherwise, the probability amplitude is reset to initial values to maintain vigilance against the deceptiveness.


### 4.1.4. 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, setxbest(t) with the current best fitness.Step 3. Use the improved map compass operator to update; then use the landmark operator to convergence; obtaingbest(t).Step 4. Comparegbest(t) with xbest(t). If they are the same, go to (a); otherwise, go to (b).
(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 optimumxbest(t). Otherwise, return to Step 3.The above steps can be summarized in Figure3.Figure 3
The procedure of the improved PIO algorithm.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).Table 1
Test function optimization results.
Test functionDAlgorithmOptimal valueTime (s)Ackley5PSO7.1431E-040.11PIO0.01370.12The improved PIO5.0849E-150.09Rastrigin5PSO5.7730E-060.13PIO0.01210.12The improved PIO00.09Rosenbrock2PSO1.0878E-090.11PIO1.0900E-080.12The improved PIO7.9993E-330.09Table 2
Parameters of global optimization methods.
ParametersPSOPIOThe improved PIOPopulation size (NP)100100100Maximum number of iterations (T)505050Map factor (R)N/A0.20.2Learning factor (c1,c2)[2,2]N/AN/ADimension (D)101010Rotating angle (°)N/AN/A-11The 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.


## 4.1.1. Map/Compass Operator

In the basic PIO algorithm, the spatial position of the pigeon is the possible solution to the optimization problem. The position of the pigeon’s home is the global optimal solution. The fitness of the pigeon is the value of the objective function. The positionXi and the speed Vi of the ith pigeon are
(17)XiT=Xi,1,Xi,2,⋯,Xi,n,i=1,2,⋯,NP,ViT=Vi,1,Vi,2,⋯,Vi,n,i=1,2,⋯,NP,where NP 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:(18)Xlow≤Xii≤Xup,Vlow≤Vii≤Vup,where i=1,2,⋯NP; Xlow and Xup are the lower and upper limits of each element in the position, respectively; and Vlow and Vup 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.(19)Vit=Vit−1⋅e−Rt+rand⋅Xg−Xit−1,(20)Xit=Xit−1+Vit,where t is the current number of iterations, R is the map compass factor, Xg 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:(21)Xit=PiPgt±L2ln1u,where
(22)PiPgt=ft×Pit+1−ft×Pgt,where t is the number of iterations; u and f are uniformly distributed random numbers on [0,1]; Pit is the historical optimal position at tth iteration; Pgt is the global optimal position at the tth iteration; and L is defined as follows:
(23)L=2ωtmbestt−Xit,where ωt is inertia weight which has great influence on the convergence of the algorithm; mbestt is the average optimal position of all pigeons in the population at tth iteration and ωt and mbestt are defined as follows:
(24)ωt=ωmax−ωmax−ωmin×tT,mbestt=1NP∑i=1NPPit,where NP is the population size; ωmax and ωmin are lower and upper limits of inertia weight; and T is the maximum number of iterations.Hence, the improved PIO algorithm updates pigeons’ position according to the following principle:(25)Xit+1=PiPgt+ωt×mbestt−Xit×ln1ft,ft≥0.5,PiPgt−ωt×mbestt−Xit×ln1ft,ft<0.5.


## 4.1.2. Landmark Operator

In the landmark operator, some individuals with better fitness are known as pigeons familiar with the landmarks, and the remaining pigeons may follow the elite pigeons or be abandoned by the population. In each iteration, pigeons are sorted according to the fitness, and half of the pigeons with poor fitness are eliminated. Then, the population centre pointCt is selected as the reference direction of the remaining pigeons and the position of the pigeon is updated.
(26)NPt=NPt−12,Xct=∑Xit⋅fitnessXitNP∑fitnessXit,Xit=Xit−1+rand⋅Xc−Xit−1,where NP is the population size; t is the number of iterations; Xc is the position of centre point; and fitness () is a fitness function, which is constructed by the objective function in Equation (10), i.e.,
(27)fitnessx=Joptmin=−mf+σP¯x.


## 4.1.3. 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:(28)q=α1⋯αi⋯αmβ1⋯βi⋯βm,where q is qubit chromosomes and αi2+βi2=1, i=1,⋯,m,m, is the length of the chromosome. αi2 and βi2 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 αi2, if the random number is larger than αi2, 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(29)αit+1=cosΔθ−sinΔθαit1−αit2,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 optimalxbestt 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=2/2,i=1,2,⋯,n. If the current global optimal solution remains after iteration, run the mutation operator to increase αi; this means that xbestt is more likely to be a global optimal solution. Otherwise, the probability amplitude is reset to initial values to maintain vigilance against the deceptiveness.


## 4.1.4. 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, setxbest(t) with the current best fitness.Step 3. Use the improved map compass operator to update; then use the landmark operator to convergence; obtaingbest(t).Step 4. Comparegbest(t) with xbest(t). If they are the same, go to (a); otherwise, go to (b).
(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 optimumxbest(t). Otherwise, return to Step 3.The above steps can be summarized in Figure3.Figure 3
The procedure of the improved PIO algorithm.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).Table 1
Test function optimization results.
Test functionDAlgorithmOptimal valueTime (s)Ackley5PSO7.1431E-040.11PIO0.01370.12The improved PIO5.0849E-150.09Rastrigin5PSO5.7730E-060.13PIO0.01210.12The improved PIO00.09Rosenbrock2PSO1.0878E-090.11PIO1.0900E-080.12The improved PIO7.9993E-330.09Table 2
Parameters of global optimization methods.
ParametersPSOPIOThe improved PIOPopulation size (NP)100100100Maximum number of iterations (T)505050Map factor (R)N/A0.20.2Learning factor (c1,c2)[2,2]N/AN/ADimension (D)101010Rotating angle (°)N/AN/A-11The 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.


## 4.2. 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 variablex, control variable u, and an output vector y. The motion equation is shown as follows:
(30)ẋ=fx,u,x0=x0,y=gx,where f and g represent any functions related to system variables; x0 is the initial value of the state variable; and y is usually used to represent the known trajectory. The above 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 xn−1 at the current point and the current control input variable un−1,k. The state variable differential of discrete points can obtain the state variable of the next discrete point through the integration of time.
(31)ẋn−1,k=fxn−1,un−1,k,(32)xn,k=∫tn−1tnẋtdt+xn−1,k,where tn is the time at node n; tn−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.,
(33)Δt=tn−tn−1.WhenΔtis fairly small, Equation (32) can be approximated as
(34)xn,k≈ẋn−1,k△t+xn−1,k.The error vectoren,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).
(35)en,k=ẋn,k−x¯̇n,(36)un−1,k+1=un−1,k−J−1en,k,where J is the jacobian matrix, which is defined as follows:
(37)J=∂eẋv∂uα∂eẋv∂uδe∂eẋv∂uϕ∂eẋγ∂uα∂eẋγ∂uδe∂eẋγ∂uϕ∂eẋqp∂uα∂eẋqp∂uδe∂eẋqp∂uϕ,where uα,uδe, and uϕ are input of the angle of attack, the elevator deflection, and the fuel equivalent ratio, respectively; and eẋv,eẋγ, and eẋqp are the errors of the estimated values of the velocity, the path angle, and the pitch rate, respectively.Through the iterative solution of Equation (35), the control variable at the n−1 discrete point can be obtained, so that the trajectory can change from the n−1 discrete point to the point n.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:
(38)mn=mn−1+ṁΔt,where mn and mn−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 Figure4.Figure 4
The procedure of the trajectory inverse algorithm.


## 5. 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 Section2 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 interpolation method to discretize the angle of attack and elevator deflection. Then, the optimization variables are expressed as follows:
(39)X=α1,α2,⋯,αD,δe1,δe2,⋯δeD,tf,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 tf 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 Table2; 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 asv,h,r,γ,m=1780m/s,22379m,0,5°,16817kg. For constraints, set the terminal constraints as vf=2337m/s, hf=26495 m, and γf=0; the path constraints are Q̇max=600kW/m2, 10kPa≤q≤100kPa, and nmax=3.5; the state constraints are 1500m/s≤v≤2500m/s and 20000m≤h≤30000m. 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:(40)−6∘≤α≤8∘,−10∘≤δe≤10∘,0≤ϕ≤1,γ̇≤0.002rad/s,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.Figure 5
State variables of optimal ascent trajectories.
(a)
Velocity(b)
Height(c)
Flight range(d)
Flight path angle(e)
MassFigure 6
Dynamic pressure, heat flux, and overload of optimal ascent trajectories.
(a)
Dynamic pressure(b)
Heat flux(c)
OverloadTable 3
Terminal states of trajectories.
AlgorithmTerminal statesPSO (with control constraints)v,h,r,γ,m=2338m/s,26485m,287863m,0.11°,15776kgPIO (with control constraints)v,h,r,γ,m=2337m/s,26495m,285903m,0.03°,15833kgThe improved PIO (with control constraints)v,h,r,γ,m=2337m/s,26495m,287800m,0°,15893kgThe improved PIO (without control constraints)v,h,r,γ,m=2337m/s,26495m,274864m,0°,15894kgWhen control constraint is considered, it can be seen from Table3 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 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 Figures5 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 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, Table4 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
Performance for different methods.
AlgorithmComputational time (s)Values of the fitness functionPSO504-15237PIO521-15612The improved PIO327-15894Figure 7
Fitness value versus iteration time.Table4 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.


### 5.2. 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 Figure8.Figure 8
Control variables obtained by trajectory inverse algorithm.
(a)
Angle of attack(b)
Rudder deflection(c)
Fuel equivalent ratioThe 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 Figure9.Figure 9
Comparison between actual trajectory and optimal trajectory.
(a)
Velocity(b)
Height(c)
Flight path angle(d)
MassThe 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.


## 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 Table2; 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 asv,h,r,γ,m=1780m/s,22379m,0,5°,16817kg. For constraints, set the terminal constraints as vf=2337m/s, hf=26495 m, and γf=0; the path constraints are Q̇max=600kW/m2, 10kPa≤q≤100kPa, and nmax=3.5; the state constraints are 1500m/s≤v≤2500m/s and 20000m≤h≤30000m. 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:(40)−6∘≤α≤8∘,−10∘≤δe≤10∘,0≤ϕ≤1,γ̇≤0.002rad/s,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.Figure 5
State variables of optimal ascent trajectories.
(a)
Velocity(b)
Height(c)
Flight range(d)
Flight path angle(e)
MassFigure 6
Dynamic pressure, heat flux, and overload of optimal ascent trajectories.
(a)
Dynamic pressure(b)
Heat flux(c)
OverloadTable 3
Terminal states of trajectories.
AlgorithmTerminal statesPSO (with control constraints)v,h,r,γ,m=2338m/s,26485m,287863m,0.11°,15776kgPIO (with control constraints)v,h,r,γ,m=2337m/s,26495m,285903m,0.03°,15833kgThe improved PIO (with control constraints)v,h,r,γ,m=2337m/s,26495m,287800m,0°,15893kgThe improved PIO (without control constraints)v,h,r,γ,m=2337m/s,26495m,274864m,0°,15894kgWhen control constraint is considered, it can be seen from Table3 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 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 Figures5 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 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, Table4 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
Performance for different methods.
AlgorithmComputational time (s)Values of the fitness functionPSO504-15237PIO521-15612The improved PIO327-15894Figure 7
Fitness value versus iteration time.Table4 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.


## 5.2. 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 Figure8.Figure 8
Control variables obtained by trajectory inverse algorithm.
(a)
Angle of attack(b)
Rudder deflection(c)
Fuel equivalent ratioThe 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 Figure9.Figure 9
Comparison between actual trajectory and optimal trajectory.
(a)
Velocity(b)
Height(c)
Flight path angle(d)
MassThe 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.


## 6. 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.



---

*Source: 3024607-2020-06-08.xml*