Ascent phase trajectory optimization for vehicle with multi-combined cycle engine based on improved particle swarm optimization

Accepted Manuscript Ascent phase trajectory optimization for vehicle with multi-combined cycle engine based on improved particle swarm optimization Hongyu Zhou, Xiaogang Wang, Yuliang Bai, Naigang Cui PII:

S0094-5765(17)30815-9

DOI:

10.1016/j.actaastro.2017.08.024

Reference:

AA 6443

To appear in:

Acta Astronautica

Received Date: 14 June 2017 Revised Date:

1 August 2017

Accepted Date: 21 August 2017

Please cite this article as: H. Zhou, X. Wang, Y. Bai, N. Cui, Ascent phase trajectory optimization for vehicle with multi-combined cycle engine based on improved particle swarm optimization, Acta Astronautica (2017), doi: 10.1016/j.actaastro.2017.08.024. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Ascent phase trajectory optimization for vehicle with multi-combined cycle engine based on improved particle swarm optimization Hongyu Zhou , Xiaogang Wang, Yuliang Bai, Naigang Cui School of Astronautics, Harbin Institute of Technology, 150001 Harbin, China. Abstract: An improved particle swarm optimization (IPSO) algorithm is proposed to optimize the ascent phase trajectory for vehicle with multi-combined cycle engine. Aerodynamic and thrust models are formulated in couple with flying states and environment. Conventional PSO has advantages in solving complicated optimization problems but has troubles in constraints handling and premature convergence preventing. To handle constraints, a modification in the fitness function of infeasible particles is executed based on the constraints violation and a comparation is executed to choose the better particle according to the fitness. To prevent premature, a diminishing number of particles are chosen to be mutated on the velocity by random times and directions. The ascent trajectory is

RI PT

divided into sub-phases according to engine modes. Different constraints, control parameters and engine models are considered in each sub-phase. Though the proposed algorithm is straightforward in comprehension and implementation, the numerical examples demonstrate that the algorithm has better performance than other PSO variants. In comparation with the commercial software GPOPS, the performance index of IPSO is almost the same as GPOPS but the results are less oscillating and dependent on initial values.

Keywords: Trajectory optimization; Combined cycle engine; Particle swarm optimization; Constraint handling; Turbulence operator

SC

1 Introduction

In recent years, the concept of reusable vehicle has obtained steadily increasing scientific and engineering attention for its convenience, flexibility and cheapness. To realize this, utilizing reusable engines instead of traditional rocket engine is a good choice. The scheme of using multi-combined cycle engine consisting of turbine based combined cycle (TBCC) and rocket based combined cycle (RBCC) is one of the best options. The vehicle takes off horizontally with the help of TBCC and the engine mode converts to

M AN U

RBCC when the predefined switch condition is met. By this way, the merits of both TBCC and RBCC are fully utilized and the defects be conquered. To this point, the vehicle with multi-combined cycle engine is of great importance for reusable vehicle and the trajectory optimization problem is of great significance.

In consideration of the complex mathematical models combined with dynamic, aerodynamic and engine models, various constraints such as dynamic pressure, heating rate, normal load, heat load, attack angle, etc. during the flight and the coupling in aerodynamics, power and flying states, the ascent phase trajectory optimization problem is of extreme complication and the optimal solution is difficult to be obtained unless an appropriate optimization method is utilized. Apart from these, the opportunity of engine mode changing is another important decision to be made. Because of the complexity of the problem, analysis optimization methods are not applicative so only numerical optimization methods are considered here.

Generally, numerical optimization methods can be classified as deterministic methods and heuristic methods. Deterministic

TE D

methods such as single shooting method, multiple shooting method, pseudospectral method, etc. are gradient-based and the gradient information can lead to fast convergence and high accuracy. Kirches [1] utilized direct multiple shooting method to discretize the optimal control problem and then solve it by quadratic programs. Chai [2] solved fuzzy physical programming problem for space maneuver trajectory optimization with hp-adaptive pseudospectral. Tian [3] proposed an indirect Legendre pseudospectral method to yield the optimal guidance law for reentry vehicles. These deterministic methods are based on the assumption of continuity and differentiability of objective function and the appropriate choice of initial values which lie in feasible region. The assumption may not be satisfied and the suitable values are usually unknown a priori. Moreover, the gradient calculation relays on numerical methods and

EP

an inappropriate method will lead to incorrect gradient information or even ill conditioning [4]. What’s more, for optimization problem in this paper, the convert opportunity of engine mode should be optimized, resulting in discontinuity in control and deficiency of

AC C

gradient information. In fact, there exist some deterministic derivative-free methods [5] and deterministic derivative-free particle

ACCEPTED MANUSCRIPT

methods [6,7]. The dependence on gradient information is well solved in these methods. Because they are rarely used in trajectory optimization of vehicles, they are not included into deterministic methods. In contrary to deterministic methods, heuristics methods such as genetic algorithm (GA), colony optimization (ACO), simulated annealing (SA) and PSO don’t need any gradient information or feasible solution knowledge [8]. Between these algorithms, Rahimi [9] demonstrates that PSO is the most outstanding one because of its convergence, robust, convenience and effectiveness. It’s introduced by Kennedy in 1995, mimicking the social behavior of flocks of birds while searching for food and takes advantage of the phenomenon of information sharing between the swarm. Each individual in the swarm, referred to as particle in PSO, is abstracted as a mess-less and volume-less point and is associated with a position vector and a velocity vector. Each particle represents a possible solution and corresponds to a performance index. Despite its simplicity and efficiency, conventional PSO has some shortcomings such as lack of consideration in constraints and the nature of premature which means getting stuck in local optimality. To conquer these and obtain better solutions, many PSO variants are proposed by researchers [10]. Most of them concentrate on the diversity keeping of the swarm, the promotion of exploitation and the treatment of constraints [11].

RI PT

Lots of effort has been taken by researchers to handle the constraints. In [12], the method based on preserving feasible solutions is adopted. The particles are initialized repeatedly until all of them lie in feasible region. Only the particles lying in feasible region are used to update the swarm, which needs a large number of particles and affects the effectiveness. In [13], penalty coefficients are incorporated into the objective function to penalize constraint violations, transforming the constrained optimization problem into an unconstrained one [14]. Tong [15] separates the objective function from the constraints and sets two fitness values to choose the leading particles. In [16], particles in infeasible region are picked out and their positions are solved again by exterior point method.

Premature convergence always exists in PSO, leading the particles to converge to a position or some specific ones. This

SC

phenomenon is mainly due to loss of diversity. The most common solution is to use a changing inertia weight instead of a fixed one. Linearly decreasing inertia weight [17], exponential decreasing inertia weight [18] and adaptive inertia weight [19] are typical ones. These methods aim to balance of the global exploration and local exploitation of the search space but the effect is unsatisfactory. An alternative method is to change the cognitive factors which decide the motion of the particles. To balance the overall exploration and local exploitation, Clerc [20] introduced a constriction factor while Zhang [21] varies the cognitive factors just like the inertia weight is

M AN U

changed. Another method is to utilize the mutation operator. Dang [22] introduced Cauchy mutation into PSO and demonstrated its effectiveness in dealing with multimodal functions. In [23], Gauss mutation is operated on every particle’s best position and the stability condition is introduced.

In this paper, the scheme of using vehicle with multi-combined cycle engine is considered and the ascent trajectory is optimized. Different constraints and control parameters are considered in accordance with the engine mode and the optimal trajectory is generated by a PSO-based method with the switch opportunity and some other variables as the unknown parameters. To resolve the problem effectively, a modification in the fitness function of particles in infeasible region and a comparation mechanism to choose the better particle are introduced and a diminishing number of particles are chosen to be mutated, solving the constraints handling problem and premature convergence, respectively. When comparing two particles to choose whose position is better, the concept of the closeness to feasible region is utilized. Under this concept, the infeasible particles are still useful so the population has not to be so large, leading to and GPOPS. 2 Description of IPSO 2.1 Basic PSO

TE D

fast convergence. Numerical examples are carried out to demonstrate the superior of IPSO when compared with other PSO variants

In PSO, each particle represents a possible solution to the optimization problem of interest and is associated with a position vector

EP

and a velocity vector. The motion of a particle is decided by the experience obtained by itself and its companions. Consider M is the size of the swarm and D is the dimension of the search space. The position and the velocity of the i th particle are donated by X i = ( xi ,1 , xi ,2 ,⋯ , xi , D ) and Vi = ( vi ,1 , vi ,2 ,⋯ , vi , D ) , respectively. Meanwhile, the best experience of the i th

AC C

particle representing the cognitive components is notated as p-best and is donated as Pi = ( pi ,1 , pi ,2 ,⋯ , pi , D ) while the best experience

ACCEPTED MANUSCRIPT

of the total swarm representing the social components is notated as g-best and is donated as Pg = ( pg ,1 , pg ,2 ,⋯ , pg , D ) . Note that in the defined swarm, there are M p-best and only a single g-best. Let the superscript t be iteration number, then the swarm is manipulated according to the flowing two equations [24]:

(

)

(

vit,+j1 = wvit, j + c1r1 pi , j − xit, j + c2 r2 pg , j − xit, j   t +1 t t +1  xi , j = xi , j + vi , j

) (1)

where w is the inertia weight and is set before the operation; c1 and c2 are positive cognitive factors and are also set before the operation; r1 and r2 are independent random numbers between 0 and 1, they are generated in every iteration. The flowchart of the basic PSO algorithm is shown in Fig. 1. For complete presentation of the PSO algorithm, the reader is

M AN U

SC

RI PT

referred to [24].

Fig. 1

Flowchart of basic PSO

TE D

The number of particles is decided according to the complexity of the problem, for simple problems dozens of particles is enough while for some complicated problems hundreds of particles may be need. The inertia weigh reflects the succession of the current velocity and influences the exploitation of the search space. The cognitive factors give a particle the ability of summarizing and learning thus leading the particle to the best position in the swarm. The settings of w , c1 and c2 are wildly studied for their enormous implications on the performance of PSO and will be mentioned below. To handle the constraints, methods based on preserving feasibility seem not so appropriate. Firstly, preserving feasible solutions needs a large number of particles, thus increasing the CPU time. Secondly, though the particle in infeasible regions can’t be used to

EP

guide the motion in the next step, the closest distance from its current position to the boundary of the feasible region can help guide the particle to the feasible region. So it’s very important to keep the percentage of infeasible particles at a certain level, which is very

AC C

difficult to realize. However, in methods based on penalty functions, the balance between fitness function and constraints function and

ACCEPTED MANUSCRIPT

the appropriate setting of penalty coefficients are difficult to yield. Neither under-penalization nor over- penalization is expected. To avoid premature, methods based on changing inertia weight or cognitive factors have to exactly catch the evolution course of the swarm which needs a lot of experiments and is difficult to be realized. If the value is too large the feasible solution will be hard to get while a too small may lead to premature convergence to a local minimum. According to the analyses above, a mechanism based on the relationship between feasible and infeasible solutions to handle constraints and a turbulence operator (some form of mutation) to the velocity to handle premature convergence are proposed and then the improved particle swarm optimization is constituted. 2.2 Mechanism to Handle Constraints In general, optimization problems can be summarized into a nonlinear programming problem: f ( x)

(2)

Subject to:  gi ( x ) ≤ 0,   h j ( x ) = 0,

i = 1, 2,⋯ , n j = 1, 2,⋯ , p

(3)

RI PT

min

where x is the vector of unknown parameters, i.e. the control variables for the optimization problem; n is the number of inequality constraints and p is the number of equality constraints. It’s common to transform equality constraints into inequality constraints: hj ( x ) − ε ≤ 0

(4)

SC

where ε is the threshold. Now only inequality constraints are considered in our study. What we do to handle constraints is to perform a little change in the fitness function when comparing two particles. The detail of the mechanism is as follows: if we compare two feasible particles, the fitness function won’t be changed and the particle with the better fitness value wins; if one is infeasible and the other is feasible, then

M AN U

the feasible particle wins; if both particles are infeasible, then the fitness function is changed (described below) and the particle with the better fitness value wins. According to this method, the p-best and g-best are generated and then the position and velocity of a particle is update.

In contrary to the method proposed by [8], in our method, the fitness of an infeasible particle isn’t set to an infinite value and the velocity of it is kept. This method utilizes the information of the closeness to the feasible region of the infeasible particles, thus avoid wasting the information of the infeasible particles. The new fitness function is as follows: m

∑u

ij

f new =

j =1 n m

(5)

∑∑ u

ij

TE D

i =1 j =1

where i donates the i th particle, n is the number of particles, j donates the j th constraint that the i th particle violates, m is the number of constraints violated by the i th particle.  hj ( x ) − ε  uij =  ε g ( x )  j

if h j ( x ) − ε > 0

(6)

if g j ( x ) > 0

EP

By this way, a smaller fitness value indicates that the particle is closer to the feasible region, despite it violates more kinds of constraints than others.

The principle is graphically simply depicted in Fig. 2. In Fig. 2, three particles and two constraints are considered: particle 1

AC C

violates two constraints, particle 2 only violates the second constraint and particle 3 only violates the first constraint. By our method,

ACCEPTED MANUSCRIPT

Particle 1 is the best because it’s the closest to the feasible region though it violates two constraints while the other two only violate one,

Fig. 2

Representation of the constraints handling mechanism

2.3 Turbulence operator

RI PT

respectively.

The turbulence operator can expend the search space of a particle, let the particle detect minor and isolated feasible regions and lead the particle to a more favorable direction. But at the end of searching, excessive mutation may cause unwanted turbulence on velocity and influence the convergence. So the number of the mutated particles should be decreased over iteration. A dynamic tm =

tc

SC

probability of turbulence operator is adopted. (7)

tmax

Pm = tm1.65 − 2tm + 1 pm = 100 × Pm

M AN U

(8) (9)

where tm and Pm are temporary variables, tc is the current iteration number (generation), tmax is the maximum iteration number, pm is the probability of mutation. The variable pm means that the top pm percent particles with better fitness value are chosen to be mutated. Thus the leader particles are lead to more favorable directions then the rest ones are guided towards better positions. Assume that tmax equals to 200, then the pm with respect to generation is shown in Fig. 3. At the beginning, most of the particles are mutated to enhance the exploration ability and prevent premature convergence, while at the end of evolution, only a littlie number of particles are mutated to enhance the exploitation ability and to accelerate the convergence.

TE D

100

Probability

80

60

40

20

EP

0 0

50

100 Generation

150

200

AC C

Fig. 3 Probability with respect to generation

ACCEPTED MANUSCRIPT

The mutation equation of velocity is as follows: vit,+j1 = 2rvθ ⋅ vit, j

(10)

where rv is a random number between 0 and 1, θ is a integer which decides the direction of turbulence and is assigned either -1 or 1. During the updating, only a part of particles’ velocity is chosen to be mutated according to eq.(10). In eq.(10), θ is an integer which decides the direction of turbulence and is assigned either -1 or 1. The random number rv is generated by uniform distribution between 0 and 1. So, the velocity is multiplied by a random number between -2.0 and 2.0. The maximum absolute value of multiplier is 2, while the average absolute value of the multiplier is 1. Therefore, the turbulence operator can be carried out without velocity climbing problem. Furthermore, in some specified problems the velocities will be bounded in order to restrict the position which will be discussed hereinafter. 2.4 The Algorithm Flow of IPSO

algorithm of the improved PSO is as follows:

RI PT

According to the process of basic PSO, the constraints handling mechanism and the turbulence operator mentioned before, the 1) Set the population of the swarm, the total generation and the threshold, then initialize the swarm;

2) Calculate the fitness values of the particles, if a particle violates at least one constraint, calculate its fitness value according to eq.(5);

3) Update the swarm with eq.(1), and then operate turbulence on the velocities of particles with better fitness value according to eq.(7) to eq.(10); 4) Go to step 2) until the maximum iteration number is reached.

TE D

M AN U

SC

The flowchart of the IPSO algorithm is shown in Fig. 4.

Fig. 4

3 Formulation of the Problem

Flowchart of IPSO

A multi-combined cycle engine constituted by TBCC and RBCC is used in this paper. The vehicle takes off horizontally, climb

EP

and accelerate under the action of TBCC until the height and velocity predefine according to the characters of TBCC and RBCC are

AC C

reached (i.e. the switch opportunity is get), then it climbs and accelerates continually under the action of RBCC till the terminal point.

ACCEPTED MANUSCRIPT

3.1 Mathematical Models 3.1.1 Mathematical Models of Vehicle

Regardless of the rotation and eccentricity of the earth, the three-degree (3DOF) point-mass dynamics of the vehicle in the longitudinal plane are given by the following equations:  dV P cos α − X − g sin γ  dt = m   dγ = P sin α + Y −  g − V  cos γ    dt mV V r  . .  dr = V sin γ  dt   dm = − mɺ  dt

RI PT

(11)

where V is the velocity scalar of the vehicle, r is the distance from the vehicle to the earth's core, γ is the flight path angle, P is the thrust, g is the gravitational acceleration, α is the attack angle, m is the mass of the vehicle, mɺ is the mass flow rate, X is the drag force and Y is the lift force. The drag force and lift force are calculated by the following equations:  X = cx qS  Y = c y qS

(12)

SC

where q is the dynamic pressure, S is the aerodynamic reference area, cx and c y are drag coefficient and lift coefficient, respectively and can be represented by functions of attack angle and Mach in the form of polynomial. ρV 2 q= (13) 2

(donated by Ma ): cx = 0.059+0.0183α +0.1332exp ( −0.2Ma ) c y = 0.088 + 0.0015α 2 + 0.54 exp ( −0.407 Ma )

M AN U

where ρ is the density of air and is the function of height. In this paper, the drag coefficient and the lift coefficient are calculated by experiential functions of attack angle and Mach

(14)

(15)

Note that in eq.(14) and (15), the unit of α is deg. In the rest of this paper, unless particularly illustrated, the units of the variables are standard units.

The trust and the mass flow rate represent the performance of the engine, which should be calculated according to the character of the engine. For vehicle engines, a variable representing the character of the fuel called specific impulse and is donated by I s instead of the mass flow rate is used to describe the engine model.

TE D

To TBCC, the trust and the specific impulse can be regarded as the function of height and Mach:  P = 109 − 12Ma − 4.7 H + 70 Ma 2 − 0.88MaH − 6 Ma3 −  1.7 Ma 2 H   I = 3790 − 211Ma + 20 H  s

(16)

where H donates the height. Note that in eq.(16), the unit of height is km and the unit of trust is kN. The RBCC has mainly four modalities: rocket ejection modality, ramjet modality, scramjet modality and rocket modality. Neglecting the mechanism that is out of the region of the research, the performance can be regarded as the outcome of the role the

EP

rocket component plays in RBCC. In fact, the more crucial role the rocket component plays, the bigger is the trust and the less is the specific impulse. So to RBCC, if the role that the rocket component plays is represented by a real number between 0 and 1 and is

AC C

donated by K , then the trust and the specific impulse can be described as follows:

ACCEPTED MANUSCRIPT

 P = Tmax K   I s = I max (1 − K )

(17)

If K equals to 0 it means the rocket component doesn’t work, while K equals to 1 means RBCC completely depends on the rocket component; Tmax is the maximum trust of RBCC and I max is the maximum specific impulse, they are set constants according to practical issues. K is called relative flow in this paper. Since TBCC and RBCC are of air-breathing style, the trust calculated by eq.(16) or eq.(17) should be multiplied by capture area donated by Sc . 3.2.2 Control Parameters and Objective Function Base on the models, the control parameters are set according to the engine mode. In TBCC the control parameter is attack angle while in RBCC the control parameters are attack angle and relative flow. The trajectory and the corresponding fitness are yield by integrating the dynamics equations with initial states and control

RI PT

parameters. To realize this, the time histories of these parameters are essential. Here, we consider the attack angle in TBCC phase as the cubic polynomial of time, the attack angle in RBCC phase as the linear polynomial of time and the relative flow in RBCC as the quadratic polynomial of time. The attack angle in TBCC phase is donated by α1 and the attack angle in RBCC phase by α 2 , then the expressions of these parameters are as follows: α1 = a1t 3 + a2 t 2 + a3t + a4  α 2 = b1t + b2  2  K = c1t + c2 t + c3

(18)

SC

where ai , bi , ci are coefficients of the polynomials and t is the fly time since each phase. Furthermore, there is another important decision to be made i.e. the switch opportunity. It’s represented by the switch height and velocity, donated by H q and Vq , respectively. vector of a particle is of 11 dimensions:

X = ( a1 , a2 , a3 , a4 , b1 , b2 , c1 , c2 , c3 , H q ,Vq )

M AN U

By far, we can get the known parameters that construct the position of a particle and will be decided by PSO. So the position

(19)

The objective is the least fuel consumption and is as follows: f =

m0 − m f

(20)

mf

where m0 and m f are the initial mass and final mass of the vehicle, respectively. Note that only minimum value of objective (fitness) is considered and optimized in this paper. 3.2 Constraints in Optimization

TE D

In this problem, the equality constraints come from the boundary conditions in both phases while the inequality constraints come from the flying course which can be divided into two aspects: state constrains and control constrains. The equality constraints contain the boundary conditions in both phase and will be assigned in next section. The inequality constraints are show in Table 1. The subscript means the phase: 1 represents the TBCC phase and 2 represents the RBCC phase, respectively. N = P sin α + Y / ( mg ) is the normal load. The symbol ‘-’ means the value is not need or doesn’t exist. The attack angle constraint in RBCC phase is more rigorous than that in TBCC phase because the RBCC may shut down if the attack angle is too big.

AC C

EP

The constraints in dynamic pressure and nominal load are derived from the endurance of thermal protection and structural strength.

ACCEPTED MANUSCRIPT Table 1 The inequality constraints Variables

Maximum value

Minimum value

Unit

α1

30

-30

deg

α2

10

-10

deg

N1

4.0

—

—

N2

3.0

—

—

q1

100

—

kPa

q2

70

—

kPa

4 Simulation and Results

RI PT

4.1 Parameters Setting

The problem is solved by employing canonical units. The distance unit is the mean radius of the earth, which equals to 6371004m and is donated by Re ; the velocity unit is the first cosmic velocity of the earth and is calculated by gravitational parameter of the earth.

µe / Re where µe is the

The equality constraints are shown in Table 2. It means that the values of the state variables should be satisfied at appointed points with errors less than the specified threshold. The symbol ‘-’ means that the value is not constrained and is derived by

Table 2

The equality constraints

V (m/s)

Switch point Terminal

150

(deg) 0

H (km) 0

m (kg) 60000

M AN U

Take off

γ

SC

optimization. The final mass can not be less than 10000kg in consideration of payload and vehicle structure.

Vq

—

Hq

—

5000

5.0

110

>10000

In PSO, usually the unknown parameters should be constrained by bounds according to the nature of the problem. They are constrained to their respective ranges: xLj ≤ xi , j ≤ xUj

( j = 1, 2,⋯ , D )

(21)

where xLk is the least values while xUk the biggest. As the positions are bounded, the corresponding velocity components are also constrained because if the velocity value is too large then, starting from any position, the updated position will beyond the border.

( j = 1, 2,⋯ , D )

(22)

TE D

xLj − xUj ≤ vi , j ≤ xUj − xLj

To this end, H q and Vq are bounded in this paper. According to eq.(16), the trust vs. height and Mach is shown in Fig. 5. It can be detected that the trust decreases as the height rises; when the height is over 20km, the trust will be quite low. At the same time, the trust rises as the Mach rises when the height is below 15km, while if the height approach 18km the trust won’t change apparently as the Mach changes and the trust will decrease as when Mach exceeds 4.0. The trust can be quite large if TBCC work at a low altitude, but low altitude means dense air, leading to enormous dynamic pressure which doesn’t benefit for climbing and work for RBCC. So the

EP

initial values of H q and Vq are restrict to [15km, 20km] and [3.0, 4.0], respectively, thus the useless search regions are ignored and

AC C

the convergence is accelerated.

ACCEPTED MANUSCRIPT 800

P / kN

600 400 200 0 0

4 10

Fig. 5

2 Mach

20 30

0

Trust of TBCC vs. height and Mach

RI PT

H /km

The position and the velocity of a particle are initialized randomly, so the initial values of the unknown parameters are generated randomly without any priori knowledge. For the unknown parameters that are bounded, their corresponding dimensions of the position are initialized randomly in the boundaries while the rest dimensions and all the velocities are initialized completely randomly. These can be done by MATLAB code easily.

The character area is set to 20m2, the capture area is set to 9m2, the maximum trust in eq.(17) is set to 200kN, the maximum specific impulse in eq.(17) is set to 2500s, and the threshold value in eq.(4) is set to 0.001.

SC

PSO with linearly decreasing inertia weight [18] (w-PSO) and PSO with constriction factor (c-PSO) [20] are used to compare with the proposed IPSO, and the constraints handling mechanism in w-PSO and c-PSO is the penalty function method described in [13]. In w-PSO, the initial inertia weight varies from 0.9 to 0.1 and the cognitive factors are both 2.0. In c-PSO, c1 and c2 are both 20 independent runs are applied for each algorithm.

M AN U

2.05. In IPSO, w is 0.5 and the cognitive factors are both 2.0. The population size is 50 and the maximum iteration is 150. A total of The best result of IPSO is also compared with the result generated by GPOPS [25]. In GPOPS, the switch height and the switch velocity are not optimized and are set according to the best result of IPSO. 4.2 Results and Discussion

Table 3 presents the results of three algorithms. It can be seen that IPSO performs better than others in all items, proving the effectiveness of IPSO. Note that in most runs the three algorithms converge before the maximum iteration number is reached. For the beat results of these algorithms, the objective value versus generation is shown in Fig. 6. Table 3

The objective values of each algorithm (20 runs) IPSO

w-PSO

c-PSO

Best value

0.4629

0.4629

0.4663

Worst value

0.5750

0.5856

0.5890

Means value

0.4994

0.5183

0.5263

Standard deviation

0.0586

0.0636

0.0647

Avg. converge generation

49

69

74

AC C

EP

TE D

Items

ACCEPTED MANUSCRIPT 0.49

IPSO w-PSO c-PSO

0.485

gbest

0.48 0.475 0.47 0.465 0.46 0

50

100

150

Fig. 6

Best objective values vs. generation

gbests worst best mean

IPSO 0.8 0.6 0.4

0

2

4

6

8

10

12

14

16

12

14

16

12

14

16

18

20

w-PSO

0.6 0.4

0

2

4

6

8

10 c-PSO

0.8 0.6

Fig. 7

0

2

4

6

8

10 runs

18

20

18

20

M AN U

0.4

SC

gbest

0.8

RI PT

generation

Results of every independent run by different algorithms

The results of every independent run, the mean value, the best value and the worst value of the three algorithms can be seen in Fig. 7. It’s obvious that compared with w-PSO and c-PSO, the proposed IPSO is more stable because it always derives a fitness near the mean value.

The swarm diversity represents the exploration and can be calculated by eq.(23); the kinetic energy represents the ability to escape local optimality and can be calculated by eq.(24).

(23)

TE D

 2 1 M D t xi , j − x jt )  Dv = ( ∑ ∑ M  i =1 j =1  M 1  t t  x j = M ∑ xi , j i =1 

where Dv donates the swarm diversity and x jt donates the average position of j th dimension in t th generation. M

D

K e = ∑∑ i =1 j =1

(x )

2 t i, j

2

(24)

EP

where K e donates the kinetic energy. For the beat result of each algorithm, the diversity and kinetic energy versus generation until convergence can be seen in Fig. 8 and Fig. 9. According to eq.(23), the diversity is indeed the standard deviation of the positions of particles, representing the disparity of

AC C

the swarm. A large diversity means big difference between particles. According to eq.(24), kinetic energy means the dispersion degree

ACCEPTED MANUSCRIPT

of particles. If the kinetic energy is large, the positions of particles will not aggregate together. When convergence occurs, diversity is close to zero and kinetic energy tends to be consistent, which can be seen in Fig. 8 and Fig. 9. Premature convergence means that the particles get together before the minimum fitness is found, i.e. diversity and kinetic energy decrease to zero and a low level too early, respectively. As a result, the swarm lost the capacity of exploitation and will converge to a local minimum. It can be seen in Fig. 8 and Fig. 9 that the diversity and kinetic energy of each algorithm are at a high level at the beginning. However, as the generation increases, the diversity of w-PSO and c-PSO decrease to a pretty small number while the diversity of IPSO keeps at a high level at the same time untie the swarm finds the minimum fitness. Analogously, the kinetic energy of w-PSO and c-PSO is much lower than IPSO, meaning that IPSO is superior in minimum exploitation. Furthermore, compared with the other two algorithms, the diversity and kinetic energy of IPSO keep fluctuating during the evolution, enhancing the exploitation capacity and preventing premature convergence further. According to the relevance between premature convergence and diversity and kinetic energy

RI PT

and the corresponding results shown in Fig. 8 and Fig. 9, the proposed IPSO is superior in premature convergence reducing.

4000 IPSO w-PSO c-PSO

3500

2500 2000 1500 1000 500 1

Fig. 8

7

x 10

5

5

Swarm diversity vs. generation

6

IPSO w-PSO c-PSO

6 kinetic energy

2 3 4 generation (log)

M AN U

0 0

SC

diversity

3000

5

4

TE D

3

2 0

50

100

150

generation

Fig. 9

Kinetic energy vs. generation

The inspection of the results generated from the best result of IPSO and GPOPS can be seen in Table 4. The time histories of

AC C

EP

some variables (e.g. height, velocity, attack angle, etc.) are shown in Fig. 10 to Fig. 16. Fig. 17 shows the profile of the trajectory.

ACCEPTED MANUSCRIPT 5000

IPSO GPOPS

Velocity /m/s

4000

3000

2000

1000

0 0

50

100 t/s

150

200

12

x 10

RI PT

Fig. 10 Velocity vs. time

4

IPSO GPOPS

10

6 4 2

50

100 t/s

Fig. 11

40 30

α /deg

20 10

150

200

M AN U

0 0

SC

Height /m

8

Height vs. time

IPSO GPOPS

TE D

0 -10

-20 0

50

AC C

EP

Fig. 12

100 t/s

150

Attack angle vs. time

200

ACCEPTED MANUSCRIPT 50

IPSO GPOPS 40

γ /deg

30

20

10

0 0

50

100 t/s

150

200

RI PT

Fig. 13 Flight path angle vs. time

1

0.8

K

0.6

SC

0.4 IPSO GPOPS

0.2

0 0

100

150

M AN U

50

t/s

Fig. 14

100

Relative flow of RBCC vs. time

IPSO GPOPS

60

40

TE D

q /kPa

80

20

0 0

50

AC C

EP

Fig. 15

100 t/s

150

Dynamic pressure vs. time

200

ACCEPTED MANUSCRIPT 3

IPSO GPOPS

Nominal Load

2

1

0

-1

50

Fig. 16

100 t/s

150

200

Nominal load vs. time

RI PT

-2 0

120 IPSO GPOPS

100

60

SC

Height /km

80

40 20 0 0

10 Mach

15

20

M AN U

5

Fig. 17 Flight profile: Mach vs. height Table 4

The best result of IPSO and the result of GPOPS Value (IPSO)

Value (GPOPS)

Fuel consumption in TBCC (kg)

2470.41

2534.21

Flight time of TBCC (s)

56.10

56.17

Fuel consumption in RBCC (kg)

37176.18

36814.46

Flight time of RBCC (s)

104.80

124.77

Switch height (km)

19.03

—

Switch velocity (m/s) / Mach

885.21 / 3.01

—

Total fuel consumption (kg)

39646.59

39348.67

TE D

Item

It can be seen that although the result of IPSO is far different from that of GPOPS, the performance index of IPSO is nearly the same as GPOPS (only 0.76% in discrepancy). In fact, the result generated by GPOPS may be violently oscillating if the setting of the be superior.

EP

package or the initial value is not appropriate, which won’t appear in IPSO. To this point, the algorithm proposed in this paper seems to

AC C

Note that in Fig. 12 and Fig. 16, the attack angle is discontinuous at the switch point, so as the nominal load which is calculated

ACCEPTED MANUSCRIPT

by attack angle, because the continuity which seems unnecessary is not considered here. In real application, the vehicle can turn on RBCC after the attitude has been changed to meet the attack angle. The constraints of dynamic pressure, nominal load and attack angle presented in Table 1 are not active and the terminal state constrains predefined for orbit injection are completely satisfied. To find the effects of number of particles, weight and factors, various population, weight and factors are utilized in IPSO. In general, the two cognitive factors are of the same value between 0 and 4, so they are assigned the same value. The best and mean fitness values (mean value after best value) under different population, weight and factor in independent 20 runs are shown in Table 5. In Table 5, the cognitive factors are represented by c for short. Analogously, the population is represented by p and weight by w. It can be seen that the best fitness value decreases as the population increases. This is because a large number of particles means greater capability of exploitation. The mean value derived by 20 particles is worse than that of 50 particles and 80 particles, while there is no obvious difference between 50 and 80. This means too little particles can’t solve the complicated trajectory problem well. Meanwhile, a big swarm can slower the convergence, so a swarm of 40 to 60 particles is appropriate for this problem. For the effect of

RI PT

weight, the best value doesn’t change very obviously under different weights. In the contrary, with respect to mean value, the result derived from the weight of 0.8 is better than that of 0.2 but worse than that of 0.5. If the weight is too small, the exploitation capability will be very weak and it’s easy to get into local minimum while a too large weight makes it difficult to converge at the end of the evolution process because the particles will jump around the minimum and finally also converge to a local minimum. For the factor, it has little impact on the performance of IPSO, so its effect can be negligible. It needs to be emphasized that other variants with respect to population, weight and factor exist in many literatures. They are not discussed here and will be researched in the future work.

SC

Results generated by IPSO under different population, weight and factor p=20

p=50

w=0.5/c=1.5

0.4870/0.5271

0.4636/0.4997

0.4574/0.4982

w=0.5/c=2.0

0.4875/0.5257

0.4629/0.4994

0.4574/0.4968

w=0.5/c=2.5

0.4863/0.5263

0.4622/0.5008

0.4972/0.4582

w=0.2/c=2.0

0.4884/0.5512

0.4627/0.5265

0.4582/0.5190

w=0.8/c=2.0

0.4871/0.5399

0.4619/0.5122

0.4575/0.5084

5 Conclusion

p=80

M AN U

Table 5

In this paper, an idea of using a reusable vehicle which takes off and lands horizontally to send payloads into space with the help of multi-combined cycle engine is proposed, utilizing the merits of TBCC and RBCC and offsetting their defects. The 3DOF point-mass dynamics of the vehicle are given so are the mathematical models of aerodynamic and dynamic system in relevant with state variables. Different constraints and control parameters are considered in accordance with the engine mode. The switch

TE D

opportunity of engine modes is another control parameter and is optimized based on an analysis of mathematical model of TBCC. In consideration of the shortcomings of PSO in constraints handling and premature which leads to local minimum, a novel mechanism based on the concept of the closeness to the feasible region and a turbulence operator that is operated according to the generation are proposed, constructing the IPSO.

A scenario in consideration of boundary constraints and course constraints with the objective of least fuel consumption is presented and is solve by IPSO, GPOPS and other two PSO variations. By analyzing and comparing the convergence, diversity and kinetic energy, it illustrates that the effectiveness and competitiveness of IPSO is better than other PSO variants. By comparing the

EP

results derived from the best result of IPSO and the optimization package GPOPS, it can be seen that though the results are far different, the performance index of IPSO is nearly the same as GPOPS, while in comparation of the oscillation of results and the dependence on

AC C

initial values, IPSO is superior to GPOPS.

ACCEPTED MANUSCRIPT

References

AC C

EP

TE D

M AN U

SC

RI PT

[1] Kirches C, Bock H G, Schloder J P, et al. Block-structured quadratic programming for the direct multiple shooting method for optimal control [J]. Optimization Methods & Software, 2011, 26(2):239-257. [2] Chai R, Savvaris A, Tsourdos A. Fuzzy physical programming for Space Manoeuvre Vehicles trajectory optimization based on hp-adaptive pseudospectral method [J]. Acta Astronautica, 2016, 123:62-70. [3] Tian B, Zong Q. Optimal guidance for reentry vehicles based on indirect Legendre pseudospectral method [J]. Acta Astronautica, 2011, 68(7):1176-1184. [4] Poustini M J, Esmaelzadeh R, Adami A. A new approach to trajectory optimization based on direct transcription and differential flatness [J]. Acta Astronautica, 2015, 107:1-13. [5] Jones D R, Perttunen C D, Stuckman B E. Lipschitzian optimization without the Lipschitz constant [J]. Journal of Optimization Theory and Applications, 1993, 79(1): 157-181. [6] Serani A, Leotardi C, Iemma U, et al. Parameter selection in synchronous and asynchronous deterministic particle swarm optimization for ship hydrodynamics problems [J]. Applied Soft Computing, 2016, 49: 313-334. [7] Pellegrini, R., Serani, A., Leotardi, C, et al. Formulation and parameter selection of multi-objective deterministic particle swarm for simulation-based optimization [J]. Applied Soft Computing, 2017, 58: 714-731. [8] Pontani M, Ghosh P. Particle Swarm Optimization of Multiple-Burn Rendezvous Trajectories [J]. Journal of Guidance Control & Dynamics, 2012, 35(4): 1192-1207. [9] Rahimi A, Kumar K D, Alighanbari H. Particle Swarm Optimization Applied to Spacecraft Reentry Trajectory[J]. Journal of Guidance Control and Dynamics, 2013, 36(1):307-310. [10] Koziel S, Michalewicz Z. Evolutionary algorithms, homomorphous mappings, and constrained parameter optimization [J]. Evolutionary Computation, 2002, 7(1): 19- 44. [11] Micheli D, Pastore R, Gradoni G, et al. Reduction of satellite electromagnetic scattering by carbon nanostructured multilayers [J]. Acta Astronautica, 2013, 88(3):61–73. [12] Pulido G T, Coello C. A constraint-handling mechanism for particle swarm optimization [C]//Evolutionary Computation, 2004. CEC2004. Congress on. IEEE, 2004: 1396-1403 Vol.2. [13] Melton R G. Hybrid methods for determining time-optimal, constrained spacecraft reorientation maneuvers [J]. Acta Astronautica, 2014, 94(1):294-301. [14] Banks A, Vincent J, Anyakoha C. A review of particle swarm optimization. Part II: hybridization, combinatorial, multi-criteria and constrained optimization, and indicative applications [J]. Natural Computing, 2008, 7(1): 109-124. [15] Tong S, Powell D J. Genetic Algorithms: A Fundamental Component of an Optimization Toolkit for Improved Engineering Designs [J]. Polyhedron, 2016, 107(2): 124- 135. [16] Yi W, Cao C, Zhang C. The geometric constraint solving based on hybrid differential evolution and particle swarm optimization algorithm [C] // International Conference on Intelligent Control and Information Processing. IEEE, 2010:692-695. [17] Sakamoto S, Oda T, Ikeda M. Implementation and evaluation of a simulation system based on particle swarm optimization for node placement problem in wireless mesh networks [J]. International Journal of Communication Networks & Distributed Systems, 2016, 17(1):1. [18] Xue W, Ma J J, Sheng W. Distributed Particle Swarm Optimization and Simulated Annealing for Energy-efficient Coverage in Wireless Sensor Networks [J]. Sensors, 2007, 7(5): 628-648. [19] Zhang P, Li J, Baoyin H, et al. A low-thrust transfer between the Earth–Moon and Sun–Earth systems based on invariant manifolds [J]. Acta Astronautica, 2013, 91(10): 77-88. [20] Clerc M, Kennedy J. The particle swarm explosion, stability, and convergence in a multi dimensional complex space [J]. Evolutionary Computation IEEE Transactions on, 2002, 6(1): 58-73. [21] Zhang Y, Cao J. 3D Human Motion Key-Frames Extraction Based on Asynchronous Learning Factor PSO [C] // Fifth International Conference on Instrumentation and Measurement, Computer, Communication and Control. IEEE, 2015: 1617- 1620. [22] Dang C T, Wu Z, Deng C. An improved approach of particle swarm optimization and application in data clustering [J]. Intelligent Data Analysis, 2015, 19(5): 1049- 1070. [23] Assarzadeh Z, Nilchi A R N. Chaotic PSO with Mutation for Classification [J]. International Journal of Tourism Research, 2015, 4(1): 29–38. [24] Wang M M, Luo J J and Walter U. Trajectory planning of free-floating space robot using Particle Swarm Optimization (PSO) [J]. Acta Astronautica, 2015, 112:77-88.

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

[25] Patterson M A, Rao A V. GPOPS-II: A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using hp-Adaptive Gaussian Quadrature Collocation Methods and Sparse Nonlinear Programming [J]. Acm Transactions on Mathematical Software, 2014, 41(1):1-37.

ACCEPTED MANUSCRIPT The highlight of the paper is as follows:

1. Trajectory optimization for vehicle with combined cycle engine is researched. 2. An improved particle swarm optimization (IPSO) is proposed solve the problem

RI PT

3. A modification in the fitness function is done to handle the constraints. 4. A turbulence operator is proposed to conquer premature convergence.

AC C

EP

TE D

M AN U

SC

5. Superiority of IPSO is certified by results in comparison with other algorithms.