Reentry trajectory rapid optimization for hypersonic vehicle satisfying waypoint and no-fl zone constraints

School of Air and Missile Defense,Air Force Engineering University,Xi’an 710051,China

Reentry trajectory rapid optimization for hypersonic vehicle satisfying waypoint and no-fl zone constraints

Lu Wang*,Qinghua Xing,and Yifan Mao

School of Air and Missile Defense,Air Force Engineering University,Xi’an 710051,China

To rapidly generate a reentry trajectory for hypersonic vehicle satisfying waypoint and no-flzone constraints,a novel optimization method,which combines the improved particle swarm optimization(PSO)algorithm with the improved Gauss pseudospectral method(GPM),is proposed.The improved PSO algorithm is used to generate a good initial value in a short time, and the mission of the improved GPM is to fin the fina solution with a high precision.In the improved PSO algorithm,by controlling the entropy of the swarm in each dimension,the typical PSO algorithm’s weakness of being easy to fall into a local optimum can be overcome.In the improved GPM,two kinds of breaks are introduced to divide the trajectory into multiple segments,and the distribution of the Legendre-Gauss(LG)nodes can be altered,so that all the constraints can be satisfie strictly.Thereby the advantages of both the intelligent optimization algorithm and the direct method are combined.Simulation results demonstrate that the proposed method is insensitive to initial values,and it has more rapid convergence and higher precision than traditional ones.

hypersonic vehicle(HV),reentry trajectory optimization,waypoint,no-flzone,particle swarm optimization(PSO), Gauss pseudospectral method(GPM).

1.Introduction

Global strike(GS)and global persistent attack(GPA)are two of the seven united states air force concepts of operations[1].The hypersonic vehicle(HV)is very important due to its capability to reenterand glide without power through the atmosphere to attack the groundtarget by only relyingon aerodynamiccontrol.As one of the most important techniquesto supportthis ability,the reentrytrajectory optimizationmethodforthe HV has beenextensivelyin attention.

Generally,numerical methods for solving this problem fall into two categories:indirect method and direct method[2].Based on the first-orde optimal necessary conditions and the Pontryagin’s minimum principle,the indirect methods are converted into two point boundary value(TPBV)problems which are notoriously sensitive to initial conditions on the co-state variables[3].Moreover, the complex constrains during the reentry phase make the solving process more complicated.Compared with the indirect methods,the direct methods convert the original optimal control problem into the non-linear programming (NLP)problem,which can be solved more easier because it does not need the first-orde optimal necessary condition andcan directlyuse numericalmethodsto optimizethe objective function[4].

In recent years,the Gauss pseudospectral method (GPM),which belongs to the direct method,has been successfully applied to the ascent trajectory reconstruction[5,6]and the reentry trajectory optimization problem [7–12]with typical constraints due to its advantages of fewer nodes used and fast convergence speed.However, during the real reentry phase,more particular situations should be taken into consideration,such as waypoint and no-fl zoneconstraints.The waypointsare positions forreconnaissance or multiple payload deployments and no-fl zones are areas that the HV must overfl due to geopolitical factors.

In the GPM,the optimal control problem is discretized at a series of specifi discretizationnodes called Legendre-Gauss(LG)nodes,all the constraints are only satisfie at each discrete node.Because the distribution of LG nodes is uneven,sparse in middle and thick in both sides,the followingproblemsmayappearwhenthewaypointandno-fl zone constraints are taken into consideration.

(i)The trajectory between each two neighboring LG nodes may not meet the requirement of the no-fl zone or path constraints.

(ii)The waypoint constraints can be satisfie with high precision unless each waypoint locates at one of the LG nodes.

(iii)The iteration speed in the early stage is slow.Although the GPM is not sensitive to the initial value,with the number of constraints increasing,the iteration speed in the early stage will become slow,especially when the no-fl zone constraints are taken into consideration.

To solve Problem 1 and Problem 2,Jiang et al.[4]proposed a multi-phase technique based on the GPM,three kinds of specifi breaks were introduced to divide the full trajectoryinto multiple phases.However,byanalyzing,the positions of quasi-contact points and turn points in his paper are fi ed,so the terminal trajectory must go through those two kinds of breaks.As a result,the fina solution may not be the optimum.In addition,Problem 3 in his paper is not solved.To solve Problem 3,Zhang et al.[13] proposed a method combining the improved genetic algorithm with sequential quadratic programming for the design of reusable launch vehicle reentry trajectory,but the waypoint and no-fl zone constraints are not taken into account.Huang et al.[14]pointed out that the trend of improving both iteration speed and precision was to combine thedirect methodwiththe intelligentalgorithm,andthe intelligent algorithm could be used as an initial value builder for the direct method.Other authors,such as Jorris et al. [15–17]and Hu et al.[18]studied the trajectory optimizationof the commonaerovehicle(CAV)satisfyingthe waypoint and no-fl zone constraints.However,the three problems above are not solved.

To overcome the three problems above,in this paper, a method combining the intelligent algorithm with the direct method is proposed.The mission of the intelligent algorithm is to generate an initial solution rapidly with low precision for the direct method,based on this initial value, the direct method can fin the fina solution with high precision quickly.

Because the precision of the fina solution is guaranteed by the direct method,so the improvement on the traditional direct method is the key point in this paper.To deal with Problem 1 and Problem 2 above,an improved GPM named multi-segment GPM(MGPM)is proposed. Two kinds of breaks are introduced to divide the whole trajectory into multiple segments.In this way,the distribution of LG nodes is altered without increasing their number.Thus,for the whole trajectory,all the constraints can be satisfied

Because the speed is the most concerned problem for a good initial value builder of the MGPM,so the intelligent algorithm should not be very complex.Among all intelligent algorithms,the particle swarm optimization(PSO) algorithm[19]has been successfully applied to wide range of areas[20–23]and many scholars from all over the world have done a lot of work to enhance the ability of the traditional PSO algorithm[24–26].In this paper,slight improvement on the traditional PSO algorithm has been done to keep its rapid convergence advantage and to overcome its weakness of being easy to fall into a local optimum.

Thereby,the advantages of both the direct method and theintelligentoptimizationalgorithmarecombined.Simulationresults show thatthe proposedmethodcan overcome all three problems above.

2.Problem formulation

In this section,the trajectory optimization model for the HV is described,which includes dynamic equations,multiple constraints and objective function.

2.1 Dynamic equations

During the reentry phase,the HV is modeled as a point mass and the earth is modeled as a spherical non-rotating sphere.The dynamic equations of the HV[27]are as follows:

where r is the radial distance from the earth center to the vehicle,V is the velocity,θ is the flight-pat angle,ψ is the heading angle,λ is the longitude,? is the latitude,γ is the bank angle,m is the mass of the vehicle,D and L are the aerodynamic drag force and lift force respectively.The equations of D and L are given by

where S is the reference area.The atmospheric density ρ=ρ0e(?h/β),where h is the altitude from the sea level, ρ0=1.225 kg/m3and β=7 100.The drag coefficien CDand lift coefficien CLare both assumed to be continuous function of the angle of attack α and Mach number Ma[28].The equations of CDand CLare as follows:

whereCLiandCDi(i=0,...,3)are the constant coeffi cient[28,29].

2.2Constraints

2.2.1 Path constraints

To meet the requirements of the thermal protection system and structure,the typical path constraints,which include dynamic pressure,heating rate and overload,should be limited during the reentry phase.The path constrains’models are as follows:

whereqL,˙QLandnLrepresent the maximum allowable values of dynamic pressure,heating rate and aerodynamic load respectively.Cis a constant andg0is the gravity coefficien at Earth surface.

2.2.2 Control constraints

To have a goodstability duringthe reentryphase,the angle of attack and the bank angle should be limited as follows:

where the subscripts“max”and“min”are the maximum and minimum allowable values,respectively.

2.2.3 Terminal constraints

In order to obtain the ability of hitting the target,the terminal state should be limited as

where the subscripts“f”means the terminal state.

2.2.4 Waypoints constraints

To accomplish the given tasks such as reconnaissance or multiple payload deployments,it is necessary for the HV to fl over some special positions which can be described by(λi,?i).The model is given by

whereNpointis the number of waypoints.

2.2.5 No-fl zone constraints

Due to the geopolitical restrictions or threat avoidance,the HV should avoid some special areas during the reentry phase.In general,the no-fl zones are specifie as cylinder zoneswith infinit altitude[16],the cross section ofthejth no-fl zone is pressed by radiusRjand center(λj,?j)[4]. The model can be described as

whereNzoneis the number of no-fl zones.

2.3Objective function

In general, the objective function is different under different missions.To have a good thermal protection,the objective function can be set as the minimum total heat or minimum peak heat.If the mission is urgent and the target’s position is fi ed,the objective functioncan beset as theminimumarrivingtime.Inaddition,it alsocan be set as the maximum crossrange,maximum downrange, etc.

2.4Trajectory optimization

Supposeu(t)∈Rmandx(t)∈Rnrepresent the control vector and the state vector respectively.The reentry trajectory optimization problem can be described as follows: finu(t)andx(t)to minimize the objective function

At the same time,the dynamics of the system,which can be described by state equations,path constraints and boundary constraints[30],should be satisfie as follows:

3.Methodology

In this section,firstl,the traditio nal PSO algorithmis analyzed,and then the EPSO algorithm for rapidly generating the initial feasible trajectory is proposed.Secondly,by introducing and analyzing the traditional GPM,the MGPM is to fin the fina solution with high accuracy.Finally, the EPSO/MGPM is presented and its advantages are discussed.

3.1Traditional PSO and analysis

3.1.1 Traditional PSO

PSO is an evolutionary optimization technique based on the competition cooperation and among individuals[31] called particles to fin the optimal solution in a particular search space.In the PSO algorithm,theith particle can be denoted by a velocity vectorVi=(vi1,vi2,...,vin)anda position vectorPi=(pi1,pi2,...,pin),wherenis the dimension size.During the evolution,each particle is updated by the personal historical best position called pbest andthe globalhistorical best positioncalled gbest.For particlei,the velocity and position can be updated by

where the superscripts“t”represents thetth iteration,ωis the inertia weight of particlei,c1is the accelerationcoeffi cient which represents the self-confidenc of particleiandc2is the other accelerationcoefficien which representsthe social confidenc oftheswarm,r1andr2are randomnumbers between 0 and 1.

3.1.2 Analysis

Although the PSO algorithm has good iteration speed in the early stage,it is easy for swarm to fall into local optimum.This problem usually happens in two phases:population initialization phase and evolution phase.

To describe the problem occurringin the population initialization phase,we present an example of ten particles’distribution.All the two kinds of distribution of swarm in Fig.1 may happen after initialization.Compared with the swarm in the left condition,the swarm in the right one is excessivelyconcentratedandcannotoverlaysolutionspace very well,so the swarm in left may have more chances to fin the global optimum.Therefore,it is important to keep the discrete degree of the initial swarm.

Fig.1 Two kinds of initial swarm’s distribution

The problemmay also happen in the evolutionphase,so we also give an example.Suppose that the dimension size of the solution is four and particleican be described byPi=(pi1,pi2,pi3,pi4).As shown in Fig.2,the distribution in each dimension of all the particles is sparse except that in dimension three.This phenomenonmay give rise to a result that all the particles in the third dimension fall into a local optimum too fast.

Fig.2 Distribution of particle swarm in each dimension

Thus if this phenomenon happens,it is necessary to make all the particles jump out of the local optimum in dimension three.

3.2EPSO

In this section,the EPSO algorithm is proposed to insure the discrete degree of the initial swarm and avoid prematureconvergencein the early stage of evolution.Firstly,the definition of code entropy and system entropy are given; secondly,the main idea of EPSO is introduced;finall,the detail algorithm fl w of EPSO is given.

3.2.1 Related definition

Entropy is the disorder degree measurement of the system. The bigger the entropy is,the more disordered the system is.In this paper,two kinds of entropy,which include code entropy and system entropy,are introduced.By adjusting the relationship between code entropy and system entropy in real-time,the requirements of keeping the discreteness of the initial swarm and avoiding premature convergence in the early stage of evolution can be satisfied The defi nitions of code entropy and system entropy are given as follows.

Suppose that the vectorDj=(p1j,...,pMj)(j= 1,2,...,n)represents the distribution of swarm in dimensionj,whereMmeans the number of particles.

Definitio1The code entropy in dimensionjis indicated byHj,which means the entropy ofDj(j= 1,2,...,n).

Suppose thatpjminandpjmaxare the minimal and maximal elements ofDj,respectively,rj=pjmax?pjminpresents the length ofDj.Then divide the interval[jmin,jmax]intosequal parts and choosenij(i= 1,2,...,s)to present the number of elements which locate in the intervali.ForDj,the density of the intervalican be calculated bybij=nij/m.Then the calculation equation ofHjis given by

Definitio2The system entropy is indicated byH. The calculation equation ofHis given as

Based on the two definition above,the proportion ofHjtoH,denoted byQj,can be calculated by

3.2.2 Main idea of EPSO

BecauseQj(j=1,2,...,n)represents the proportion ofHjtoH,ifQjis toosmall,it meansthat distributionoftheparticle swarm in dimensionjis too dense but it is sparse in other dimensions.

In the population initialization phase,Qjis too small, it means that the distribution in dimensionjis not dense enough.Themethodthat wechoosetoovercomethis problem is to reinitialize the swarm in dimensionj.

When the phenomenon above happens during the early evolution phase,it means that the iteration speed of all the particles in dimensionjis much faster than that in other dimensions.This phenomenon may bring a result that the swarm in dimensionjmay fall into local optimum.The method to overcome this problem is to change thegbestin dimensionj.According to(11),as shown in Fig.3,if thegbestis changed,all the particles in this dimensionmay searchtowardsa newdirection.Inthisway,allthe particles have bigger probability to jump out of the local optimum and fin the global optimum.The equation to changegbestis given as wherer3is a randomnumberbetween–1 and1.gjmaxandgjminare the maximum and minimum allowable values in dimensionj.

Fig.3 Sketch map of gbest updating

3.2.3 EPSO application

(i)Encoding method for particle

Bydividingtheinterval[t0,tf]inton?1equalparts,the number of nodes isn.Then,the particle can be described as

where for particlei,tifrepresents the total flyin time,αijmeans the angleof attack at nodej,γijmeans the bankangleatnodej.Thevelocity,pbestandgbestcanbedescribed as(17).

(ii)Objective function

Inthis paper,we use thepenaltyfunctionmethodto deal with the inequality constraints in the model.

(iii)Detail fl w of EPSO

The detail fl w of EPSO is as follows:

Step 1Initialize necessary parameters for the algorithm;

Step 2Initialize the swarm randomly;

Step 3For each dimensionj,calculateQj,ifQj＜Qmin,reinitialize the swarm randomly in dimensionjand then return to Step 3,otherwise,go to Step 4;

Step 4For each particlei,calculate the fitnes value ofpiand updateli;

Step 5Updateg;

Step 6If the terminal condition is satisfiedgis the solution of the problem,otherwise,go to Step 7;

Step 7For each particlei,updateviandpiaccording to(9);

Step 8For each dimensionj,calculateQj,ifQj＜Qmin,updategjaccording to(13)and updateviandpiaccording to(9),then calculateQjand return to Step 8; otherwise,return to Step 4.

3.3Traditional GPM and analysis

3.3.1 Traditional GPM

In the GPM,the optimal control problem is discretized at LG nodes and then converted into a nonlinear programming problemby approximatingthe states and controls using Lagrange interpolating polynomials[32].

Since the total flyin time for the HV lasts fromt0totfbut the LG nodes belong to the range[?1 1],so in order to describe the problem inτ∈[?1 1],it is necessary to transformt∈[t0tf]toτ∈[?1 1]bythetransformation

Furthermore,theNdiscrete nodes contain the initial nodeτ0=?1,the terminal nodeτf=1 andK=N?2 interior LG nodes which are the roots of theKth-degree Legendre polynomial.Then the state and control variables can be approximated as

Convert(19)to its derivative

In the GPM,the derivatives of each state at the LG node shouldsatisfy(21).ThedifferentialofeachLagrangepolynomial at the LG nodes can be expressed in a differential approximation matrixD∈RN×N+1and the elements of the matrix can be calculated by the following equation [33].

wherePN(τ)is theNth degree Lagendre polynomial,i=1,2,...,N,k=1,2,...,N.

Then(10)are transformed as

wherek=1,2,...,N,because only the state at each LG node should satisfy(23).Therefore,the state at each boundary node[30]should be added as

whereωkis the Gauss weights.

Then(9)can be approximated as follows:

with boundary and path constraints

Finally,the optimal control problem is converted to an NLP which is described as follows:fin the NLP variables

to minimizes(25)and subject to(23)and(26).

3.3.2 Analysis

As shown in Fig.4,the distribution of LG nodes is uneven, sparse in middle and dense in both sides.

Fig.4 Distribution of LG nodes

The distribution of Gauss points may give rise to some problems as follows:

(i)Because the passage time for each waypoint is uncertain but the distribution of LG nodes is certain in the interval[–1 1],so the waypoints constraints can be satisfie unless each waypoint locates at one of the LG nodes, as shown in Fig.5.

Fig.5 Example of trajectory near the waypoints

(ii)As shown in Fig.6,although the no-fl zone constraints are satisfie at each LG node,the segment between each two neighboringLG nodes may not meet the requirements of the no-fl zone constraints.

Fig.6 Example of trajectory near the no-fl zone

(iii)The peak value of all the constraints may be less accurate even though all the constraints are satisfie at each LG node.As shown in Fig.7,the peak values exceed the given range.

Fig.7 Example of constraints near peak value

3.4MGPM

In this section,the MGPM is proposed to overcome the problem described above.Firstly,the definition of the fi ed node and virtual node are given;secondly,the main ideaofEPSOispresented;finall,thedetailalgorithmfl w of EPSO is given.

3.4.1 Related definition

Definitio3The fi ed node is define as the node which onlyhas positioninformationbut does not have passage time information,such as the waypoint.

Definitio4The virtual node is define as the node which onlyhas passage time informationbut does not have position information.

3.4.2 Main idea of MGPM

The main idea of the MGPM can be described as the following steps:

Step 1Divide trajectory by fi ed nodes

As shown in Fig.8,the whole trajectory is divided into multiple segments by fi ed nodes(waypoints).In this way,the waypoints constraints can be satisfied Then go to Step 2.

Fig.8 Trajectory divided by fixe nodes

Step 2Solve model

Suppose that the whole trajectory is divided intoZsegments and the number of LG nodes in each segment isN. The objective function in segmentzcan be described as

Similarly,all the constraints in segmentzcan be transformed as

Then the objective function for the whole trajectory is transformed as follows

Toconnecteachtwoneighborsegments,theboundaryconditions for each segment should be added as

Thus,the GPM is transformed into the MGPM which is described by(28)–(31).The processing mode of the MGPM is similar as that of GPM.In this paper,the sequential quadratic programming(SQL)algorithm,which is one of the best methods to solve the nonlinear optimization problem,is chosen to solve the finall converted NLP. Then solve the problem by SQL and go to Step 3.

Step 3Check solution

As shown in Fig.5–Fig.7,although the constraints are satisfie at each LG node,it does not mean the segment of trajectory between each two neighboring LG nodes could meet the requirements of all the constraints.Here,it is necessary to check the feasibility of the solution by calculating all the constraints at the checkpoints which locate in the middle of each two neighboring LG nodes.Suppose thatrepresents the location of the checkpointi.The distribution of LG nodes and checkpoints is shown in Fig.9.The calculation equation ofis given as

If all the constraints are satisfie at each checkpoint,the algorithm stops and the solution is the fina solution of the problem,or else go to Step 4.

Fig.9 Distribution of LG nodes and checkpoints

Step 4Divide trajectory by additional virtual nodes

After Step 3,the checkpoints that do not meet the requirements of any constraints can be regarded as virtual nodes.Then,as shown in Fig.10,the whole trajectory is divided into multiple segments by The fl w of the EPSO/ MGPM is as follows:firstl,run the EPSO to generate a reentry those virtual nodes.Because of the distribution character of LG nodes,there will be more LG nodes located near the virtual nodes.In this way,all the constraints may be satisfied Then return to Step 2.

Fig.10 Sketch map of the influenc by inserting virtual nodes

3.5EPSO/MGPM and analysis

trajectory with low precision;secondly,obtain the initial values of the NLP at each discrete node based on the result of EPSO;finall,solve the problem by the MGPM based on the initial values above.

Because the proposed method combines the direct method with the intelligent algorithm,so the advantages of the EPSO-MGPM can be summarized as follows:

(i)The speed of the algorithmhas a significan improvementcomparedwiththetraditionaldirectmethod.Because the improved PSO algorithm is insensitive to initial values and more rapidly convergence,the EPSO algorithm can give a feasible solution near the optimum rapidly.Based on the initial values,the fina solution with high accuracy can be found quickly.

(ii)Theaccuracyofthealgorithmis improvedcompared with the traditional direct method.

Via dividing the trajectory into multiple segments by fi ed nodes,each waypoint will become an LG node.Because all the constraints can be satisfie strictly at each LG node,so the waypoint constraints can also be satisfie with high precision.

Via dividing the trajectory into multiple segments by virtual nodes,the distribution of LG nodes is altered and more LG nodes will locate near the no-fl zone or other appointedplaces.As a result,the no-fl zoneand path constraints can be satisfie with high precision.

4.Experiments and results

In this section,we evaluate the performance of the proposed EPSO/MGPM approach by two experiments.The firs one is about trajectory optimization with typical constraints by EPSO and EPSO/MGPM to demonstrate the contribution of EPSO in the EPSO/MGPM approach.The second one is on trajectory optimization with complex constraints to show the good performance of the MGPM compared with the GPM and the physical programming approach[34].The related parameters throughout all the simulations are asS=0.35 m2,m=907 kg,αmax= 90°,αmin=?90°,γmax=90°,γmin=?90°,˙Qmax= 2 000 kW/m2,qmax=500 kPa,nmax=4 g.The optimization results are found through MATLAB 7.14.

4.1Trajectory optimization with typical constraints by EPSO and EPSO/MGPM

To be an initial value builder,EPSO should have fast itera-tion speed and the ability to approach global optimum. In order to prove the contribution of EPSO in the EPSO/MGPM,we make a simulation of maximum downrange trajectory with typical constraints,which does not include waypoints or no-fl zones,by using the EPSO and EPSO/MGPM.The initial and terminal conditions of the simulation are shown in Table 1.The objective function is described as J=maxλf.Moreover,the number of particles is 20 and the number of LG nodes is 60.

Table 1 Initial and terminal condition of maximum downrange trajectory

Fig.11 Results of optimal trajectory by EPSO and EPSO/MGPM

As shown in Fig.11,the results of optimal trajectory using EPSO is close to the fina results of the optimal trajectory using the EPSO/MGPM.Although the accuracy of EPSO is poorer than that of the MGPM because the maximum longitude is 101.92°by EPSO and 103.04°by EPSO/MGPM(Fig.11(d)),the result by EPSO is close to the result by EPSO/MGPM.It should be noted that the running time is about 20 s by EPSO and 30 s by EPSO/MGPM,but when used by the GPM,it is about 4 min.With the constraints increased,the runningtime difference between the GPM and the EPSO/MGPM will become bigger.That is the reasonwhy the EPSO algorithmis chosen to be an initial value builder for the MGPM in the EPSO/MGPM approach.

Fig.11(a)–Fig.11(e)show that the trajectory is smooth,and the boundary constraints such as the altitude and velocity are satisfied The path constraints,such as the dynamic pressure,heating rate,and aerodynamic load,are also in the allowable range Fig.11(f)–Fig.11(h).

4.2 Trajectory optimization with complex constraints by EPSO/MGPM and other traditional methods

In this section,the comparison with other two traditional methods is implemented.The traditional methods we choose are the physical programmingapproachand the GPM,which are representatives for the numerical method and the non-numericalmethodrespectively.The initial and terminal conditions of the simulation are shown in Table 2. Theparametersofthewaypointandno-fl zoneconstraints are shown in Table 3.The objective function is described asJ=maxtf.

Table 2 Initial and terminal condition of minimum time trajectory

Table 3 Parameters of waypoint and no-fl zone constraints

Fig.12 and Fig.13 show the optimal trajectories by the physical programming approach and the traditional GPM respectively.And the result of the optimal trajectoryby the EPSO/GPM is shown in Fig.14.

Fig.12 Results of optimal trajectory by physical programming approach

Fig.13 Results of optimal trajectory by GPM

Fig.14 Results of optimal trajectory which divided by fixenodes and virtual nodes

Fig.12 shows that the optimal trajectory generated by the physical programming approach is smooth(Fig.12 (a)and Fig 12(b))and no-fl zone constraints are satisfie(Fig.13(c)and Fig.13(d)).However,the precision of the physical programming approach is lower than that of the EPSO/MGPM by comparing with Fig.14(c)and Fig.14(d).It also demonstrates that the numerical method is better than the non-numerical method on precision.

Fig.13 shows that the optimal trajectory is also smooth (Fig.13(a)and Fig.13(b)),but all the no-fl zone constraints can only be satisfie at each LG node(Fig.13(c) and Fig.13(d)).Thus the trajectorygeneratedby the GPM is not available.However,the trajectory generated by the EPSO/MGPM satisfie all the no-fl zone constraints between each two neighboring LG nodes(Fig.14(c)and Fig.14(d)).The path constraints,such as the dynamic pressure,heating rate,and aerodynamic load,are also in the allowable range(Fig.14(h)–Fig.14(j)).

Furthermore,the objective function values are 1 820 s by the physical programming approach,1 730 s by the GPM and 1 762 s by the EPSO/MGPM.And the running time is 2 min,5 min and 40 s,respectively.

4.3 Discussion

What we can conclude from the simulation results are as follows:

(i)To be an initial value builder for the MGPM in the EPSO/MGPM approach,EPSO can generate a feasible solution near the optimum quickly and the speed of the algorithm has a significan improvement compared with othertraditional methods.

(ii)The way points are used as the segment breaks,in this way,the way point constraints are transformed into the terminal constraints in each segment.Thus the waypoint constraints can be satisfied

(iii)Via dividing the trajectory into multiple segments by virtual nodes,there will be more LG nodes locating near the no-flzones or other appointed places.By this technique,the no-fl zone and path constraints can be satisfied

In addition,because the MGPM only adjusts the distribution of LG nodes but not increases the number of LG nodes,so the EPSO/MGPM can obtain the optimal solution with fewer LG nodes than other multi-segment methods,and can decrease the computational complexity remarkably.

5.Conclusions

In this paper,the EPSO/MGPM is proposed to solve the reentry trajectory optimization problem which includes waypoint and no-fl zone constraints.The MGPM is devised to ensure the terminal solution to fulfil the requirements of the waypoint,no-flzone and path constraints with high precision.And the EPSO is used to generate a goodinitialfeasiblesolutioninashorttimefortheMGPM. So thattheadvantagesofboththedirectmethodandthe intelligent optimization algorithm are combined.Simulation results show that by using EPSO,the speed of the algorithm has a significan improvement;by using the MGPM, all the complex constraints between each two neighboring LG nodes can be satisfied

[1]USAF.Air force handbook-109th congress.Washington,D.C.: Department of the Air Force,2005.

[2]G.T.Huntington,D.Benson,A.V.Rao.A comparison of accuracy and computational efficien y of three pseudospectral methods.Proc.of the AIAA Guidance,Navigation,and Control Conference and Exhibit,2007,AIAA 20–23.

[3]G.Q.Huang,Y.P.Lu,Y.Nan.A survey of numerical algorithms for trajectory optimization of fligh vehicles.Science China Technological Sciences,2012,55(9):2538–2560.

[4]Z.Jiang,Z.Rui.Reentry trajectory optimization for hypersonic vehicle satisfying complex constraints.Chinese Journal of Aeronautics,2013,26(6):1544–1553.

[5]Q.Zong,B.L.Tian,L.Q.Dou.Ascent phase trajectory optimization for near space vehicle based on Gauss pseudospectral method.Journal of Astronautics,2010,31(7):1775–1781.(in Chinese)

[6]D.Zhang,X.F.Lu,L.Liu,et al.An online ascent phase trajectory reconstruction algorithm using Gauss pseudospectral method.Proc.of the International Conference on Modeling, Identificatio and Control,2012:1184–1189.

[7]X.X.Yang,Y.Yang,J.H.Shi.A survey of boost-glide aircraft trajectory design.Journal of Naval Aeronautical and Astro-nautical University,2012,27(3):245–252.(in Chinese)

[8]P.C.Zhang,X.G Li.Based on trajectory optimization of glide phase with GPM for hypersonic glide vehicle.Electronic Design Engineering,2013,21(15):105–108.(in Chinese)

[9]L.Tang,J.M.Yang,F.Y.Chen.Trajectory optimization with multi-constraints for a lifting reentry vehicle.Missiles and Space Vehicles,2013,1(1):1–5.(in Chinese)

[10]H.J.Cao,Z.L.Ge.A study on optimizing design of trajectory for glide extend-range projectile.Journal of Sichuan Ordnance,2013,34(11):8–10.(in Chinese)

[11]K.Li,W.S.Nie,B.M.Feng.Research on multi-phase trajectory optimizationfor boost-glide vehicle.Command Control& Simulation,2012,34(5):21–25.(in Chinese)

[12]E.M.Yong,G.J.Tang,L.Chen.Rapid trajectory optimization for hypersonic reentry vehicle via Gauss pseudospectral method.Journal of Astronautics,2008,29(6):1766–1771.(in Chinese)

[13]D.N.Zhang,Y.Liu.Reentry trajectory optimization based on improved genetic algorithm and sequential quadratic programming.Journal of Zhejiang University(Engineering Science), 2014,48(1):161–167.(in Chinese)

[14]C.Q.Huang,H.F.Guo,D.L.Ding.A Survey of trajectory optimization and guidance for hypersonic gliding vehicle.Journal of Astronautics,2014,35(4):369–379.(in Chinese)

[15]T.R.Jorris,R.G.Cobb.Multiple method 2-D trajectory optimization satisfying waypoints and no-flzone constraints.Journal of Guidance,Control,and Dynamics,2008,31(3): 543–552.

[16]T.R.Jorris,R.G.Cobb.Three-dimensional trajectory optimization satisfying waypoint and no-flzone constraints.Journal of Guidance,Control,and Dynamics,2009,32(2): 551–572.

[17]T.R.Jorris.Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and no-fl zone constraints.Ohio,USA:Air Force Institute of Technology,2007.

[18]Z.D.Hu.Research on trajectory planning and guidance for space-based strike weapon.Changsha:National University of Defense Technology,2009.

[19]J.Kennedy,R.C.Eberhart.Particle swarm optimization.Proc. of the IEEE International Conference on Neural Networks, 1995:1942–1948.

[20]Y.F.Liao,D.H.Yau,C.L.Chen.Evolutionary algorithm to traveling salesman problems.Computer and Mathematics with Applications,2012,64:788–797.

[21]W.C.Yeh,Y.C.Lin,Y.Y.Chung,et al.A particle swarm optimization approach based on Monte Carlo simulation for solving the complex network reliability problem.IEEE Trans. on Reliability,2010,59(1):212–221.

[22]H.Cho,D.Kim,F.Olivera,et al.Enhanced speciation in particle swarm optimization for multi-modal problems.European Journal of Operational Research,2011,213:15–23.

[23]Y.Wang,Y.Yang.Particle swarm with equilibrium strategy of selection for multi-objective optimization.European Journal of Operational Research,2010,200:187–197.

[24]H.J.Li,Y.Yi,X.X.Li,etal.Humanactivityrecognition based on HMM by improved PSO and event probability sequence.Journal of Systems Engineering and Electronics,2013,24(3): 545–554.

[25]M.W.Li,H.G.Kang,P.F.Zhou,et al.Hybrid optimization algorithm based on chaos,cloud and particle swarm optimization algorithm.Journal of Systems Engineering and Electronics,2013,24(2):324–334.

[26]J.J.Qi,B.Kuo.H.T.Lei,et al.Solving resource availabilitycost problem in project scheduling by pseudo particle swarm optimization.Journal of Systems Engineering and Electronics, 2014,25(1):69–76.

[27]N.X.Vinh.Optimal trajectories in atmospheric fligh.New York:Elsevier Scientifi Publishing Company,1981.

[28]Y.Sun,G.R.Duan,M.R.Zhang.Modifie aerodynamic coefficien models of hypersonic vehicle in reentry phase.System Engineering and Electronics,2011,33(1):134–137.(in Chinese)

[29]T.H.Phillips.A common aero vehicle(CAV)model,description,and employment guide.Schafer Corporation for AFRL and AF-SPC,2003.

[30]B.L.Tian,Q.Zong.3DOF ascent phase trajectory optimization for aircraft based on adaptive gauss pseudospectral method.Proc.of the 3rd International Conference on Intelligent Control and Information Processing,2012:431–435.

[31]X.C.Zhao,W.Q.Lin,Q.F.Zhang.Enhanced particle swarm optimization based on principal component analysis and line search.Applied Mathematics and Computation,2014,229: 440–456.

[32]G.T.Huntington.Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems.Cambridge,United States:Massachusetts Institute of Technology, 2007.

[33]Y.Liu,Y.J Qian,J.Q.Li,et al.Mars exploring trajectory optimization using gauss pseudo-spectral method.Proc.of the International Conference on Mechatronics and Automation, 2012:2371–2377.

[34]E.M.Yong,L.Chen,G.J.Tang.Trajectory optimization of hypersonic gliding reentry vehicle based on the physical programming.Acta Aeronautica et Astronautica Sinica,2008, 29(5):1092–1097.

Biographies

Lu Wangwas born in 1987.He received his M.S.degree from Air Force Engineering University (AFEU)in 2013.He is currently pursuing his Ph.D. in military operational research at AFEU.His research interests include trajectory optimization,target tracking,information fusion and sensor resource management.

E-mail:18049689239@163.com

Qinghua Xingwas born in 1966.She received her Ph.D.degree from Air Force Engineering University in 2003.She is a professor.Her research interests include systems modeling and simulation,operation decision analysis in air and missile defense.

E-mail:liuxqh@126.com

Yifan Maowas born in 1988.She received her M.S.degree from Air Force Engineering University (AFEU)in 2013.She is currently pursuing her Ph.D. in military operational research at AFEU.Her research interests include trajectory optimization,target tracking and information fusion.

E-mail:myf1210@126.com

10.1109/JSEE.2015.00140

Manuscript received January 27,2015.

*Corresponding author.

This work was supported by the National Natural Science Foundation of China(61272011).