# Legendre Cooperative PSO Strategies for Trajectory Optimization

Authors: Lei Liu; Yongji Wang; Fuqiang Xie; Jiashi Gao
Journal: Complexity
 (2018)

---

## Abstract

Particle swarm optimization (PSO) is a population-based stochastic optimization technique in a smooth search space. However, in a category of trajectory optimization problem with arbitrary final time and multiple control variables, the smoothness of variables cannot be satisfied since the linear interpolation is widely used. In the paper, a novel Legendre cooperative PSO (LCPSO) is proposed by introducing Legendre orthogonal polynomials instead of the linear interpolation. An additional control variable is introduced to transcribe the original optimal problem with arbitrary final time to the fixed one. Then, a practical fast one-dimensional interval search algorithm is designed to optimize the additional control variable. Furthermore, to improve the convergence and prevent explosion of the LCPSO, a theorem on how to determine the boundaries of the coefficient of polynomials is given and proven. Finally, in the numeral simulations, compared with the ordinary PSO and other typical intelligent optimization algorithms GA and DE, the proposed LCPSO has traits of lower dimension, faster speed of convergence, and higher accuracy, while providing smoother control variables.


---

## Body

## 1. Introduction

Swarm intelligence is a collective dynamic behavior of distributed, self-organized systems, natural or artificial, employed in work on artificial intelligence. It introduces many simple agents with very general rules to achieve an “intelligent” global optimal behavior. Swarm intelligence-based techniques can be used in a number of applications on optimization. The US military is investigating the swarming techniques to control unmanned vehicles. The European Space Agency is thinking about an orbital swarm for self-assembly and interferometry. NASA is investigating the use of swarm technology for planetary mapping.In particular, trajectory optimization problem is one of the most important tasks in the preliminary design of the next generation of high speed vehicles, such as NASA’s X-43 unmanned hypersonic vehicle (HV), and has a great effect on the choice of conceptual design [1]. It is a daunting work to get the solution of nonlinear optimal control problem with the arbitrary final time and multiconstraints.Trajectory optimization with multiconstraints has been linked with some stochastic search algorithms. The typical approaches include genetic algorithms (GA) [2], differential evolution (DE) [3], and particle swarm optimization (PSO) [4–6]. Besides, several novel bionic optimization algorithms spring up in these years and show remarkable efficiency in the industrial domain, such as the honeybee mating optimization [7], harmony search algorithm [8], and ants swarm optimization [9]; however, they are rarely used in the aerospace field because of their excessive novelty. The PSO algorithm, as an optimization algorithm based on swarm intelligence, similar to GA and DE for their origination from population-based heuristic search mechanisms, has recently become more popular due to its simplicity and effectiveness. It is widely used because of its simple principle, small number of parameters to be adjusted, and easy realization. Considering that the dynamic parameter optimization needs to be transformed into the static format, it is validated to evaluate the trajectory parameters optimization using a spline function [10].Under the guidance of thought of the hybrid algorithm [11], a novel method of trajectory optimization is proposed to improve the global convergence ability in this paper, the Legendre cooperative PSO (LCPSO) method, which is a kind of cooperative PSO based on Legendre orthogonal polynomials.In LCPSO, the following improvements make it promising in solving complex problems. Firstly, the search space is divided into certain subspaces and different swarms are arranged to optimize the different parts of the search space. In this way, the optimized scale for each swarm can be reduced directly. The algorithm is suitable for the problems with larger scale and higher dimension.Secondly, for the optimal control problem by PSO, the discretization is necessary before solving the parameter optimization problem. Herein, the Legendre orthogonal polynomial approximation can achieve higher smoothness of control variables with lower dimensions. Additionally, a theorem on how to find the precise range of optimal parameters is given and proven, as well as how to find the boundaries of the coefficient of polynomials. Therefore, the proposed LCPSO is expected to realize the optimization with higher accuracy and efficiency.Thirdly, in order to solve the optimal control problems with arbitrary final time rather than the fixed one, the proposed Legendre orthogonal polynomial approximation method introduces an additional control variable to transcribe the original optimal problem to the one with fixed final time. Then, a traditional one-dimensional search method based on the interval analysis is proposed to optimize the additional control variable. This way, the specific optimal problem with single boundary can be solved.Finally, we use two typical trajectory optimization problems to illustrate efficiency of the proposed LCPSO algorithm. One is the ascent trajectory optimization of X-43 hypersonic vehicle, and the other one is the classic optimal orbit transfer problem. The simulation results demonstrate the advantages of the LCPSO in terms of solution accuracy and convergence rate by comparing with some traditional intelligent optimization algorithms.The organization of the remainder of this paper is as follows. Section2 formulates the optimal control problem by 3-DOF mass point dynamics. Section 3 describes the algorithm of novel Legendre cooperative PSO and some properties. Section 4 presents the numerical simulations on the performance of the proposed optimal algorithms. Section 5 draws conclusions to the paper.


## 2. Problem Description

HV is one of the main workhorses for most of the nations of the world for scientific studies and military and commercial applications [1]. Efforts are ongoing in the 21st century to enhance flexibility and reliability and reduce the overall cost of such systems. Here, without loss of generality, we provide the proposed algorithm to solve a trajectory optimal problem of HV. Its 3-DOF dynamics over a spherical earth are described by the following motion equations [12, 13]:(1)R˙=Vsin⁡γ,θ˙=Vcos⁡γsin⁡ψRcos⁡ϕ,ϕ˙=Vcos⁡γcos⁡ψR,V˙=Tcos⁡αm-D-sin⁡γR2,γ˙=1VTsin⁡αm+L+V2-1Rcos⁡γR,m˙=fMa,α,δ.The 3-DOF dynamics described in (1) are dimensionless equations of motion with six state variables and two control variables. The real variables are normalized as follows:(2)R=rr0,V=vg0r0,g=g0r02r2,where r is the radial distance from the earth’s center to the vehicle, θ is the longitude, ϕ is the latitude, v is the velocity of the vehicle, γ is the flight path angle, ψ is the velocity azimuth angle measured clockwise from the north, the control variables are the angle of attack α and the fuel throttle opening δ, m is the mass of the vehicle, and the terms D and L are the aerodynamic drag and lift acceleration, which are defined by(3)D=12ρv2CDScmg0,L=12ρv2CLScmg0.Here, Sc is the reference area of the vehicle. The terms of CD(α,Ma,δ) and CL(α,Ma,δ) are the coefficients of drag and lift, which are also functions of α and δ.The control variablesα and δ are approximated by Legendre orthogonal polynomial functions, respectively. Assuming LN(x) is the Legendre polynomial of degree N, the continuous function u(x) of variable x can be approximated as(4)ux=∑n=0NbnLnx.This way, the dynamic parameter optimization problem is transformed into a static one before the PSO algorithm performed.The evaluation of the performance index starts when the scramjet is launched. The vehicle releases from the boost phase and then enters into the later ascent and cruise phase. Here, we define the switching time astk. By analyzing the dynamic characteristics of HV, the fuel throttle opening δ shall be minimized during the later ascent stage.In this problem, the optimal time of ascent stage is free but satisfying some constraints. Thus, an additional control variable is introduced to transcribe the proposed problem to a fixed final time problem. The switching step is available and then the ascent trajectory can be divided into two phases. In the first phase, only the fuel throttle opening needs to be optimized to guarantee that it is a minimized constant value at the beginning of the second phase. In order to simplify the equations and problem, we give a reasonable assumption that the switching step satisfiestk=kt′f,  k∈[0,1], with a fixed time interval [t0,t′f], and the system’s continuous state does not jump at each switching point. Therefore, the arbitrary terminal time optimal control problem is turned into a two-point boundary value problem with the fixed terminal time. The optimal time tf would be obtained as long as all the terminal conditions are fulfilled.The design spaceS of HV ascent trajectory optimization using LCPSO can be defined by the following vector:(5)S=b1α,…,bNα,b1δ,…,bNδ,k.Here, the left two parts in the vectorS represent the two groups of coefficients of the orthogonal polynomials which approximate the control variables α and δ, while the last parameter k represents the ratio between the switching step and the supposed final time. Hence, the problem space can be divided into three subspaces bα, bδ, and k. Then, in LCPSO, we provide two groups of particles to optimize the subspaces bα and bδ, respectively, and use a one-dimensional search method based on interval analysis to optimize the variable k.To seek the solutions of both the continuous inputu(t) and the switching step tk according to the given initial state x(t0)=x0, we define a performance function as follows:(6)J=Φxtf+∫t0tkLxt,utdt+∫tktfL′xt,utdt.In this case, the above nonlinear optimal problem can be described as a basic optimal control problem of nonlinear Bolza type.Considering that the performance function here is selected to minimize the fuel expendable ratio (FER) through the control variablesα and δ, normally, we use the following relationship to describe it directly:(7)J=-∫t0tfm˙dt=m0-mf.Furthermore, the ascent trajectory terminate conditions can be represented by the following inequalities:(8)rf-Δrdown≤rtf,vf-Δvdown≤vtf,γf-Δγdown≤γtf≤γf+Δγup,where r(tf), v(tf), and γ(tf) are three state variables of the terminal point, while the tolerance is given by the upper bound Δγup and the lower bounds Δrdown,  Δvdown, and Δγdown.The terminal conditions are added to the performance function to ensure that all the constraints are satisfied. Then, the performance function of LCPSO is defined as follows:(9)Jα,δ=J+k1erf-rtf+k2evf-vtf+k3γ2tf.Here, in (9), the first part J represents the optimal index of minimum FER described in (7). The three terms k1e(rf-r(tf)),  k2e(vf-v(tf)), and k3γ2(tf) represent the terminal constraints index. r(tf) is the radial distance from the center of the earth to the vehicle and v(tf) represents the velocity at the final time, while the flight path angle index is γ(tf).Considering that the proposed LCPSO provides two groups of particles based on Legendre orthogonal polynomial approximation, the assumption of the fixed time interval[t0,t′f] can be transformed into a closed interval [-1,1] of Legendre orthogonal polynomials to facilitate subsequent processing. Meanwhile, the additional switching step tk can make the optimal interval of δ more accurate. After that, the PSO algorithm is employed to solve the optimal control problem about δ in the fixed interval [t0,t′f] and the next interval [tk, tf]. Furthermore, the 4th-order Runge-Kutta numeric integration is employed to evaluate the performance function, while the optimal time tf can be obtained as long as all the terminal conditions are satisfied.


## 3. Algorithm Formulation

### 3.1. LCPSO

The proposed LCPSO is designed as a double iteration scheme: the inner iteration of dual cooperative PSO (CPSO) [14] and the outer iteration of interval analysis. The mathematical formulation of the CPSO algorithm is redesigned in (10) and (11).(10)xidrlt+1=xidrlt+vidrlt+1,with d∈{1,…,N},  l∈{1,…,S},  t<Tmax,  r∈{1,2,…,Rmax},(11)vidrlt+1=w1×vidrlt+c1×rand×pldrlt-xidrlt+c2×rand×pgdrlt-xidrlt,w1=wmin+wmax-wmin×Niter-tNiter,where xidrl(t) is the coordinate of particle i in dimension d at time t as subinterval l at interval iteration step r.The pseudocode of CPSO algorithm is reported in Pseudocode1.Pseudocode 1:Pseudocode of CPSO.
( 1 ) Input P (number of cooperation iteration), ( 2 )    K (number of cooperation iteration), ( 3 )    Nmax (subswarm size), ( 4 )    k, l and r (the parameters of interval analysis) are given from outer iteration. ( 5 ) Initialize (b01α,…,b0Nα) and (b01δ,…,b0Nδ) ( 6 ) For cooperation iteration j∈[1,…,P] ( 7 )  For each subswarm iteration t∈[1,…,Tmax] ( 8 )   For each particle i∈[1,…,Nmax] ( 9 )    If (j mod 2 > 0)  {xrlid(0)=(bαj1,…,bαjN)} ( 10 )    Else{xrlid(0)=(bδj1,…,bδjN)} ( 11 )   End If ( 12 )   To run Runge-Kutta numeric integration and evaluate performance function by using Equ. (9). ( 13 )   If (J(xirl(t)<pbestval)) ( 14 )     {pbest=xirl(t), pbestval=J(xirl(t))} ( 15 )   Else If(J(xirl(t)<gbestval)) ( 16 )   {gbest=xirl(t), gbestval=J(xirl(t))} ( 17 )   End If ( 18 )   To update the particle by using Equ. (10), Equ. (11). ( 19 )  End Fori ( 20 ) End Fort ( 21 ) If(j mod 2>0) ( 22 )   {To update (bαj1,…,bαjN)=gbest} ( 23 ) Else{To update (bδj1,…,bδjN)=gbest} ( 24 ) End If 25 End For jThe mathematical formulation of the one-dimensional search method based on interval analysis is designed as(12)krl=kr-1+l-0.5S+1wkr-1,withl∈1,…,S,r∈1,2,…,Rmax,wkr=wkr-1S.The one-dimensional search algorithm based on the interval analysis is reported in Pseudocode2.Pseudocode 2:Pseudocode of interval analysis.
( 1 ) Input T f, S, R, tk0, ( 2 ) Initialize t k 0 and w(k0) by using Equ. (12), ( 3 ) For interval analysis iteration r∈[1,…,R] ( 4 )  For each subinterval l∈[1,…,S] ( 5 )   To update the switch instant by using Equ. (11). ( 6 )   To perform CPSO and update J(krl) ( 7 )  End Forl ( 8 )   To update kr={krl∣l:min⁡(j(krl))} ( 9 )   To update trk=krTf. ( 10 ) End For r ( 11 ) To obtain the best position and the best fitness.The midpoint and the initial width of the switching step interval are defined as(13)k0=tk0Tf-1-tk0/TfS2,(14)wk0=S1-tk0Tf,where krl is the ratio between the switching step and the assumption final time in the interval iteration r and the subinterval l. The switch step interval is divided into S subintervals at each interval iteration step r.The method starts with the initial switch step interval, which is split into multiple subdivisions. Then, those subdivisions are either sent to the solution list which are considered later, or removed from further test list by certain cut-off condition. The above process is repeated by choosing a new switch step interval until no switch step could be considered or a global optimal point is found.


### 3.2. Theorem of Boundaries Selection of LCPSO

The Legendre orthogonal polynomials can be generated using Gram-Schmidt orthonormalization [15] in the interval [-1,1] with the function of the weight ρ(x)=1. LN(x) represents a Legendre polynomial with degree N. We provide the following transformation in (15) of the independent variable τ of 3-DOF dynamics when τ∈[τ0,τf] to guarantee that -1≤x≤1:(15)x=2τ-τf+τ0τf-τ0.This way, the continuous function u(x) of variable x can be approximated as(16)ux=∑n=0NbnLnx.Then, using Legendre polynomial approximation, the original optimal control problem can be transformed into a parameter optimization problem: to find the optimal coefficientsb=(b0,b1,…,bN) so as to minimize the fitness in (9) and satisfy the terminal constraints in (8).Obviously, the appropriate ranges of the coefficientsb can improve the convergence of the optimization procedure and prevent explosion. In our formal studies, a theorem to determine the coefficient of the Legendre orthogonal polynomial is proposed in [16]. However, the statement and proof of the theorem were incomplete and less rigorous. The complete expression of the theorem of how to determine the boundaries of the coefficients of the orthogonal polynomials is proposed and proven.Theorem 1.
Assume thatLN(x) represents an Nth degree Legendre polynomial. The continuous function u(x) of variable x can be approximated as in (16). And the range of u(x) belongs to the interval [Umin,Umax] with the variable x presented in (15). Then, the ranges of optimal coefficients b=(b0,b1,…,bN) can be calculated by(17)bn∈Umin,Umax,n=0,bn∈-n+12Umax-Umin∫-11Ln+xdx,n+12Umax-Umin∫-11Ln+xdx,n≥1,or(18)bn∈Umin,Umax,n=0,bn∈n+12Umax-Umin∫-11Ln-xdx,-n+12Umax-Umin∫-11Ln-xdx,n≥1,where (19)Lnx=1,n=012nn!dndxnx2-1n,n=1,2,…,(20)Ln+x=Lnx,Lnx≥0,0,otherwise,(21)Ln-x=Lnx,Lnx<0,0,otherwise.Proof.
According to the assumption, we have the following equations:(22)∫-11Lkxuxdx=∫-11Lkx∑n=0NbnLnxdx=bk∫-11LkxLkxdx+∑n=0,n≠kNbn∫-11LkxLnxdx.
According to the orthogonal relationship of the Legendre polynomials, we have(23)∫-11LkxLnxdx=22k+1δnkin which δnk is called Kronecker Delta and satisfies (24)δnk=1,n=k0,otherwise.
Substituting (23) and (24) into (22), we have(25)bk=k+12∫-11uxLkxdx,withk∈0,N.
Assume thatf(x) and g(x) are continuous in [a,b] and satisfy f(x)≤g(x); then, we have(26)∫abfxdx≤∫abgxdxexcept the finite numbers of points.
Since the range ofu(x) is [Umin,Umax] and from (20) and (21), we have(27)Lnx=Ln+x+Ln-x.Then, the following inequations should be achieved:(28)UminLn+x≤Ln+xux≤UmaxLn+x,UmaxLn-x≤Ln-xux≤UminLn-x.Hence,(29)UmaxLn-x+UminLn+x≤Ln+x+Ln-xux,Ln+x+Ln-xux≤UmaxLn+x+UminLn-x.
According to (26), (27), and (29), we can obtain(30)∫-11Lnxuxdx≥Umax∫-11Ln-xdx+Umin∫-11Ln+xdx,∫-11Lnxuxdx≤Umax∫-11Ln+xdx+Umin∫-11Ln-xdx.
Moreover, according to (25) and (30), we have the effective range of bn:(31)bn∈2n+12Umin∫-11Ln+xdx+Umax∫-11Ln-xdx,2n+12Umax∫-11Ln+xdx+Umin∫-11Ln-xdx.Therefore, the upper bound and lower bound of coefficients bn can be defined:(32)bnmax=2n+12Umax∫-11Ln+xdx+Umin∫-11Ln-xdx,bnmin=2n+12Umin∫-11Ln+xdx+Umax∫-11Ln-xdx,with n=0,1,2,….
The closed form for the orthogonal polynomials is given in (19). Assume that(33)gx=x2-1n,n=1,2,….Then,(34)∫-11Lnxdx=∫-11Ln+xdx+∫-11Ln-xdx=∫-1112nn!gnxdx=12nn!xxx-11=0.Therefore,(35)∫-11Ln+xdx=-∫-11Ln-xdx,n=1,2,….
Substituting (35) into (32), we have(36)bnmax=2n+12Umax-Umin∫-11Ln+xdx,bnmin=-bnmax,n≥1,or (37)bnmax=-2n+12Umax-Umin∫-11Ln-xdx,bnmin=-bnmax,n≥1.When n=0, importing L0(x)=1 and -1≤x≤1 into (32), we have(38)Umin≤b0≤Umax.The proof is achieved.Consider that the control valueu(x) always works in a symmetric feasible interval. Assuming that Umax>0 and Umin=-Umax, to simplify the results in practical application, we can obtain(39)n=0,b0∈Umin,Umax;n=1,b1∈1.5Umin,1.5Umax;n=2,b2∈1.9245Umin,1.9245Umax;n=3,b3∈2.275Umin,2.275Umax;n=4,b4∈2.5793Umin,2.5793Umax;n=5,b5∈3.8895Umin,3.8895Umax.


## 3.1. LCPSO

The proposed LCPSO is designed as a double iteration scheme: the inner iteration of dual cooperative PSO (CPSO) [14] and the outer iteration of interval analysis. The mathematical formulation of the CPSO algorithm is redesigned in (10) and (11).(10)xidrlt+1=xidrlt+vidrlt+1,with d∈{1,…,N},  l∈{1,…,S},  t<Tmax,  r∈{1,2,…,Rmax},(11)vidrlt+1=w1×vidrlt+c1×rand×pldrlt-xidrlt+c2×rand×pgdrlt-xidrlt,w1=wmin+wmax-wmin×Niter-tNiter,where xidrl(t) is the coordinate of particle i in dimension d at time t as subinterval l at interval iteration step r.The pseudocode of CPSO algorithm is reported in Pseudocode1.Pseudocode 1:Pseudocode of CPSO.
( 1 ) Input P (number of cooperation iteration), ( 2 )    K (number of cooperation iteration), ( 3 )    Nmax (subswarm size), ( 4 )    k, l and r (the parameters of interval analysis) are given from outer iteration. ( 5 ) Initialize (b01α,…,b0Nα) and (b01δ,…,b0Nδ) ( 6 ) For cooperation iteration j∈[1,…,P] ( 7 )  For each subswarm iteration t∈[1,…,Tmax] ( 8 )   For each particle i∈[1,…,Nmax] ( 9 )    If (j mod 2 > 0)  {xrlid(0)=(bαj1,…,bαjN)} ( 10 )    Else{xrlid(0)=(bδj1,…,bδjN)} ( 11 )   End If ( 12 )   To run Runge-Kutta numeric integration and evaluate performance function by using Equ. (9). ( 13 )   If (J(xirl(t)<pbestval)) ( 14 )     {pbest=xirl(t), pbestval=J(xirl(t))} ( 15 )   Else If(J(xirl(t)<gbestval)) ( 16 )   {gbest=xirl(t), gbestval=J(xirl(t))} ( 17 )   End If ( 18 )   To update the particle by using Equ. (10), Equ. (11). ( 19 )  End Fori ( 20 ) End Fort ( 21 ) If(j mod 2>0) ( 22 )   {To update (bαj1,…,bαjN)=gbest} ( 23 ) Else{To update (bδj1,…,bδjN)=gbest} ( 24 ) End If 25 End For jThe mathematical formulation of the one-dimensional search method based on interval analysis is designed as(12)krl=kr-1+l-0.5S+1wkr-1,withl∈1,…,S,r∈1,2,…,Rmax,wkr=wkr-1S.The one-dimensional search algorithm based on the interval analysis is reported in Pseudocode2.Pseudocode 2:Pseudocode of interval analysis.
( 1 ) Input T f, S, R, tk0, ( 2 ) Initialize t k 0 and w(k0) by using Equ. (12), ( 3 ) For interval analysis iteration r∈[1,…,R] ( 4 )  For each subinterval l∈[1,…,S] ( 5 )   To update the switch instant by using Equ. (11). ( 6 )   To perform CPSO and update J(krl) ( 7 )  End Forl ( 8 )   To update kr={krl∣l:min⁡(j(krl))} ( 9 )   To update trk=krTf. ( 10 ) End For r ( 11 ) To obtain the best position and the best fitness.The midpoint and the initial width of the switching step interval are defined as(13)k0=tk0Tf-1-tk0/TfS2,(14)wk0=S1-tk0Tf,where krl is the ratio between the switching step and the assumption final time in the interval iteration r and the subinterval l. The switch step interval is divided into S subintervals at each interval iteration step r.The method starts with the initial switch step interval, which is split into multiple subdivisions. Then, those subdivisions are either sent to the solution list which are considered later, or removed from further test list by certain cut-off condition. The above process is repeated by choosing a new switch step interval until no switch step could be considered or a global optimal point is found.


## 3.2. Theorem of Boundaries Selection of LCPSO

The Legendre orthogonal polynomials can be generated using Gram-Schmidt orthonormalization [15] in the interval [-1,1] with the function of the weight ρ(x)=1. LN(x) represents a Legendre polynomial with degree N. We provide the following transformation in (15) of the independent variable τ of 3-DOF dynamics when τ∈[τ0,τf] to guarantee that -1≤x≤1:(15)x=2τ-τf+τ0τf-τ0.This way, the continuous function u(x) of variable x can be approximated as(16)ux=∑n=0NbnLnx.Then, using Legendre polynomial approximation, the original optimal control problem can be transformed into a parameter optimization problem: to find the optimal coefficientsb=(b0,b1,…,bN) so as to minimize the fitness in (9) and satisfy the terminal constraints in (8).Obviously, the appropriate ranges of the coefficientsb can improve the convergence of the optimization procedure and prevent explosion. In our formal studies, a theorem to determine the coefficient of the Legendre orthogonal polynomial is proposed in [16]. However, the statement and proof of the theorem were incomplete and less rigorous. The complete expression of the theorem of how to determine the boundaries of the coefficients of the orthogonal polynomials is proposed and proven.Theorem 1.
Assume thatLN(x) represents an Nth degree Legendre polynomial. The continuous function u(x) of variable x can be approximated as in (16). And the range of u(x) belongs to the interval [Umin,Umax] with the variable x presented in (15). Then, the ranges of optimal coefficients b=(b0,b1,…,bN) can be calculated by(17)bn∈Umin,Umax,n=0,bn∈-n+12Umax-Umin∫-11Ln+xdx,n+12Umax-Umin∫-11Ln+xdx,n≥1,or(18)bn∈Umin,Umax,n=0,bn∈n+12Umax-Umin∫-11Ln-xdx,-n+12Umax-Umin∫-11Ln-xdx,n≥1,where (19)Lnx=1,n=012nn!dndxnx2-1n,n=1,2,…,(20)Ln+x=Lnx,Lnx≥0,0,otherwise,(21)Ln-x=Lnx,Lnx<0,0,otherwise.Proof.
According to the assumption, we have the following equations:(22)∫-11Lkxuxdx=∫-11Lkx∑n=0NbnLnxdx=bk∫-11LkxLkxdx+∑n=0,n≠kNbn∫-11LkxLnxdx.
According to the orthogonal relationship of the Legendre polynomials, we have(23)∫-11LkxLnxdx=22k+1δnkin which δnk is called Kronecker Delta and satisfies (24)δnk=1,n=k0,otherwise.
Substituting (23) and (24) into (22), we have(25)bk=k+12∫-11uxLkxdx,withk∈0,N.
Assume thatf(x) and g(x) are continuous in [a,b] and satisfy f(x)≤g(x); then, we have(26)∫abfxdx≤∫abgxdxexcept the finite numbers of points.
Since the range ofu(x) is [Umin,Umax] and from (20) and (21), we have(27)Lnx=Ln+x+Ln-x.Then, the following inequations should be achieved:(28)UminLn+x≤Ln+xux≤UmaxLn+x,UmaxLn-x≤Ln-xux≤UminLn-x.Hence,(29)UmaxLn-x+UminLn+x≤Ln+x+Ln-xux,Ln+x+Ln-xux≤UmaxLn+x+UminLn-x.
According to (26), (27), and (29), we can obtain(30)∫-11Lnxuxdx≥Umax∫-11Ln-xdx+Umin∫-11Ln+xdx,∫-11Lnxuxdx≤Umax∫-11Ln+xdx+Umin∫-11Ln-xdx.
Moreover, according to (25) and (30), we have the effective range of bn:(31)bn∈2n+12Umin∫-11Ln+xdx+Umax∫-11Ln-xdx,2n+12Umax∫-11Ln+xdx+Umin∫-11Ln-xdx.Therefore, the upper bound and lower bound of coefficients bn can be defined:(32)bnmax=2n+12Umax∫-11Ln+xdx+Umin∫-11Ln-xdx,bnmin=2n+12Umin∫-11Ln+xdx+Umax∫-11Ln-xdx,with n=0,1,2,….
The closed form for the orthogonal polynomials is given in (19). Assume that(33)gx=x2-1n,n=1,2,….Then,(34)∫-11Lnxdx=∫-11Ln+xdx+∫-11Ln-xdx=∫-1112nn!gnxdx=12nn!xxx-11=0.Therefore,(35)∫-11Ln+xdx=-∫-11Ln-xdx,n=1,2,….
Substituting (35) into (32), we have(36)bnmax=2n+12Umax-Umin∫-11Ln+xdx,bnmin=-bnmax,n≥1,or (37)bnmax=-2n+12Umax-Umin∫-11Ln-xdx,bnmin=-bnmax,n≥1.When n=0, importing L0(x)=1 and -1≤x≤1 into (32), we have(38)Umin≤b0≤Umax.The proof is achieved.Consider that the control valueu(x) always works in a symmetric feasible interval. Assuming that Umax>0 and Umin=-Umax, to simplify the results in practical application, we can obtain(39)n=0,b0∈Umin,Umax;n=1,b1∈1.5Umin,1.5Umax;n=2,b2∈1.9245Umin,1.9245Umax;n=3,b3∈2.275Umin,2.275Umax;n=4,b4∈2.5793Umin,2.5793Umax;n=5,b5∈3.8895Umin,3.8895Umax.


## 4. Simulation Results

### 4.1. Ascent Trajectory Optimization of the X-43 Vehicle

Firstly, we use the model of the X-43 vehicle [13] as the test case. Particularly, the vector of the initial states is [r0,θ0,φ0,v0,γ0,m0]=[6398000,120,0,1200,0,1000], while the terminal states are [rf,vf,γf]=[6408000,1820,0]. The tolerance of state variables in the terminal point is [Δr,Δv,Δγ]=[±1000,±50,±1.5]. Some other experiment parameters are as follows: the number of LCPSO iterations that occurred T=35, the population size of a single particle swarm N=40, the dimension of the particle D=6, the number of swarm iterations of the cooperated particle P=4, the acceleration coefficients c1=c2=2, and the range of the inertia weight wmax=0.9 and wmin=0.4. In addition, the weights of the performance function are [k1,k2,k3]=[109,104,103], and the parameters of interval analysis iteration are [t′f,tk,S,R]=[1,0.9,5,3].In Table1, during the optimal processing, five subintervals are introduced and the optimal costs vary with the switch steps in each iteration step of interval analysis, and the results converged at the global optimal point only through three steps of iteration. It takes about 580 seconds under the hardware environment of Intel Core i3 3220 CPU and DDR4G 1600 MHz memory and the software environment of Windows 7 operating system with VC++ 6.0 compiling environment.Table 1
Results after 1 time numerical simulation.
Iteration timet Switching step ratiok Corresponding steptk Optimal costJ 1 0.450 90.0 320.849 1 0.550 110.0 314.069 1 0.650 130.0 291.451 1 0.750 150.0 433.650 1 0.850 170.0 491.392 2 0.610 122.0 293.675 2 0.630 126.0 254.625 2 0.650 130.0 275.729 2 0.670 134.0 294.822 2 0.690 138.0 348.983 3 0.622 124.4 285.098 3 0.626 125.2 259.907 3 0.630 126.0 254.625 3 0.634 126.8 262.411 3 0.638 127.6 278.003The optimal flight path angle profile is shown in Figure1, with the corresponding optimal control variables profile in Figures 2 and 3.Figure 1
The profile of opening of fuel valve.Figure 2
The profile of the angle of attack.Figure 3
The profile of the angle of flight path.In Table2, according to the ratios of the feasible results, the best and the mean results of 200 times simulations using interval analysis of LCPSO with different coefficients’ ranges are presented. Apparently, the proposed method with the coefficients’ ranges selection has the best performance, while all state variables meet the state terminal constraints with higher precision. Furthermore, the feasible results distribution with standard range and half-range selection is shown in Figure 4, while the feasible results distributions with standard, two, and five times range selection are displayed in Figure 5. From the comparison, we can find that if the optimal solutions existed in this reductive area, the reduction of the search range to its half will probably decrease the search complexity. However, if the performance index function corresponding to the reductive area is of bad quality, the optimized solution is likely to be lost. On the other hand, the multifold coefficient will provide the unwanted search space, which will increase the complexity and reduce efficiency of the optimization.Table 2
Coefficient range selection results.
Standard range 0.5 times range 2 times range 5 times range Ratio of feasible results 72.0 59.1 63.5 24.5 Least consumption of fuel 51.8% 55.9% 53.6% 56.4% Mean consumption of fuel 63.0% 61.1% 64.5% 70.5%Figure 4
Result distribution with standard and half times range.Figure 5
Result distribution with standard, double times, and five times range.


### 4.2. Earth-Mars Transfer Problem

The orbit transfer problem is another kind of hotspot issues on trajectory optimization research [17]. Here, we use the proposed LCPSO method to solve the Earth-Mars trajectory transfer problem and compare it with the traditional DE and GA method.The Earth-Mars transfer orbit can be divided into three stages: geocentric escape section, heliocentric transition section, and areosynchronous capture section. Each segment has different constraints. In the geocentric escape section, the initial state constraints are[r0,vr0,vθ0]=[6.6107RE,0,0.38893RE/TUE] and the terminal state constraint is rf=145RE. On the other hand, the constraints required in the third segment are the terminal radius ρf and the terminal speeds vρf and vδf, which satisfy [ρf,vρf,vδf]=[6.0236RM,0,0.40745RM/TUM]. And the initial radius is the radius of the sphere of influence for Mars, ρ0=170RM.The LCPSO algorithm parameters are set as follows: the population sizem=20 and the maximum number of iterations n=200. Meanwhile, the parameters of DE and GA algorithm are set as follows: the population size m=20, the maximum number of iterations n=200, the cross probability pc=0.4, and the mutation probability pm=0.6. The optimization variables contain the thrust angle of the aircraft, the initial and terminal orbital angles, and the time of the whole three-segment orbits flight. Currently, the range of the thrust angle and the initial/terminal orbital angles of the vehicle is set in [0,π]. The problem is to transfer the orbit within the shortest time to minimize the fuel consumption, so the fitness function can be described as J=t1+t2+t3, where t1, t2, and t3 are the transfer periods of three sections, respectively. Figure 6 shows the best fitness process using LCPSO, DE, and GA algorithms. It is shown that the LCPSO has a faster speed of convergence and higher accuracy than DE and GA algorithms.Figure 6
Convergence processes of the fitness functions.The final resulting candidate solutions by LCPSO, DE, and GA are shown in Figures7–9. Because the intelligent algorithm obtains the feasible solution to meet the accuracy requirements, the state curves of any two algorithms are not necessarily completely consistent. But both of them satisfy the constraints and accuracy requirements. For instance, as shown in Figure 7, although the final radius of the geocentric orbit is about 145RE by LCPSO, 145.1RE by GA, and 143RE by DE, the orbital radius and the velocity are all within the permitted constrained ranges. And in this case, the optimization result of the final geocentric orbital radius by LCPSO is closer to the constraint requirements. Hence, the LCPSO exhibits better efficiencies on the constraints requirements than the other two algorithms. Moreover, the results of the whole transfer time by different optimization algorithms are shown in Table 3. It can be found that, in the case of satisfying constraints, the optimal trajectory by LCPSO will spend the shortest flight time among the three algorithms, which has the best performance.Table 3
Final times for the geocentric, heliocentric, and areocentric segments.
Segment Final time (days) LCPSO GA DE 1 (geocentric) 37.1 34.6 34.5 2 (heliocentric) 171.4 182.4 174.7 3 (areocentric) 19.3 20.3 19.6 Whole flight time 227.8 237.3 228.8Figure 7
State variables in geocentric segment.Figure 8
State variables in heliocentric segment.Figure 9
State variables in areocentric segment.


## 4.1. Ascent Trajectory Optimization of the X-43 Vehicle

Firstly, we use the model of the X-43 vehicle [13] as the test case. Particularly, the vector of the initial states is [r0,θ0,φ0,v0,γ0,m0]=[6398000,120,0,1200,0,1000], while the terminal states are [rf,vf,γf]=[6408000,1820,0]. The tolerance of state variables in the terminal point is [Δr,Δv,Δγ]=[±1000,±50,±1.5]. Some other experiment parameters are as follows: the number of LCPSO iterations that occurred T=35, the population size of a single particle swarm N=40, the dimension of the particle D=6, the number of swarm iterations of the cooperated particle P=4, the acceleration coefficients c1=c2=2, and the range of the inertia weight wmax=0.9 and wmin=0.4. In addition, the weights of the performance function are [k1,k2,k3]=[109,104,103], and the parameters of interval analysis iteration are [t′f,tk,S,R]=[1,0.9,5,3].In Table1, during the optimal processing, five subintervals are introduced and the optimal costs vary with the switch steps in each iteration step of interval analysis, and the results converged at the global optimal point only through three steps of iteration. It takes about 580 seconds under the hardware environment of Intel Core i3 3220 CPU and DDR4G 1600 MHz memory and the software environment of Windows 7 operating system with VC++ 6.0 compiling environment.Table 1
Results after 1 time numerical simulation.
Iteration timet Switching step ratiok Corresponding steptk Optimal costJ 1 0.450 90.0 320.849 1 0.550 110.0 314.069 1 0.650 130.0 291.451 1 0.750 150.0 433.650 1 0.850 170.0 491.392 2 0.610 122.0 293.675 2 0.630 126.0 254.625 2 0.650 130.0 275.729 2 0.670 134.0 294.822 2 0.690 138.0 348.983 3 0.622 124.4 285.098 3 0.626 125.2 259.907 3 0.630 126.0 254.625 3 0.634 126.8 262.411 3 0.638 127.6 278.003The optimal flight path angle profile is shown in Figure1, with the corresponding optimal control variables profile in Figures 2 and 3.Figure 1
The profile of opening of fuel valve.Figure 2
The profile of the angle of attack.Figure 3
The profile of the angle of flight path.In Table2, according to the ratios of the feasible results, the best and the mean results of 200 times simulations using interval analysis of LCPSO with different coefficients’ ranges are presented. Apparently, the proposed method with the coefficients’ ranges selection has the best performance, while all state variables meet the state terminal constraints with higher precision. Furthermore, the feasible results distribution with standard range and half-range selection is shown in Figure 4, while the feasible results distributions with standard, two, and five times range selection are displayed in Figure 5. From the comparison, we can find that if the optimal solutions existed in this reductive area, the reduction of the search range to its half will probably decrease the search complexity. However, if the performance index function corresponding to the reductive area is of bad quality, the optimized solution is likely to be lost. On the other hand, the multifold coefficient will provide the unwanted search space, which will increase the complexity and reduce efficiency of the optimization.Table 2
Coefficient range selection results.
Standard range 0.5 times range 2 times range 5 times range Ratio of feasible results 72.0 59.1 63.5 24.5 Least consumption of fuel 51.8% 55.9% 53.6% 56.4% Mean consumption of fuel 63.0% 61.1% 64.5% 70.5%Figure 4
Result distribution with standard and half times range.Figure 5
Result distribution with standard, double times, and five times range.


## 4.2. Earth-Mars Transfer Problem

The orbit transfer problem is another kind of hotspot issues on trajectory optimization research [17]. Here, we use the proposed LCPSO method to solve the Earth-Mars trajectory transfer problem and compare it with the traditional DE and GA method.The Earth-Mars transfer orbit can be divided into three stages: geocentric escape section, heliocentric transition section, and areosynchronous capture section. Each segment has different constraints. In the geocentric escape section, the initial state constraints are[r0,vr0,vθ0]=[6.6107RE,0,0.38893RE/TUE] and the terminal state constraint is rf=145RE. On the other hand, the constraints required in the third segment are the terminal radius ρf and the terminal speeds vρf and vδf, which satisfy [ρf,vρf,vδf]=[6.0236RM,0,0.40745RM/TUM]. And the initial radius is the radius of the sphere of influence for Mars, ρ0=170RM.The LCPSO algorithm parameters are set as follows: the population sizem=20 and the maximum number of iterations n=200. Meanwhile, the parameters of DE and GA algorithm are set as follows: the population size m=20, the maximum number of iterations n=200, the cross probability pc=0.4, and the mutation probability pm=0.6. The optimization variables contain the thrust angle of the aircraft, the initial and terminal orbital angles, and the time of the whole three-segment orbits flight. Currently, the range of the thrust angle and the initial/terminal orbital angles of the vehicle is set in [0,π]. The problem is to transfer the orbit within the shortest time to minimize the fuel consumption, so the fitness function can be described as J=t1+t2+t3, where t1, t2, and t3 are the transfer periods of three sections, respectively. Figure 6 shows the best fitness process using LCPSO, DE, and GA algorithms. It is shown that the LCPSO has a faster speed of convergence and higher accuracy than DE and GA algorithms.Figure 6
Convergence processes of the fitness functions.The final resulting candidate solutions by LCPSO, DE, and GA are shown in Figures7–9. Because the intelligent algorithm obtains the feasible solution to meet the accuracy requirements, the state curves of any two algorithms are not necessarily completely consistent. But both of them satisfy the constraints and accuracy requirements. For instance, as shown in Figure 7, although the final radius of the geocentric orbit is about 145RE by LCPSO, 145.1RE by GA, and 143RE by DE, the orbital radius and the velocity are all within the permitted constrained ranges. And in this case, the optimization result of the final geocentric orbital radius by LCPSO is closer to the constraint requirements. Hence, the LCPSO exhibits better efficiencies on the constraints requirements than the other two algorithms. Moreover, the results of the whole transfer time by different optimization algorithms are shown in Table 3. It can be found that, in the case of satisfying constraints, the optimal trajectory by LCPSO will spend the shortest flight time among the three algorithms, which has the best performance.Table 3
Final times for the geocentric, heliocentric, and areocentric segments.
Segment Final time (days) LCPSO GA DE 1 (geocentric) 37.1 34.6 34.5 2 (heliocentric) 171.4 182.4 174.7 3 (areocentric) 19.3 20.3 19.6 Whole flight time 227.8 237.3 228.8Figure 7
State variables in geocentric segment.Figure 8
State variables in heliocentric segment.Figure 9
State variables in areocentric segment.


## 5. Conclusion

In this paper, a novel interval analysis based Legendre cooperative PSO algorithm is proposed and applied to solve the trajectory optimization problems. The Legendre orthogonal polynomial approximation is synthesized with the dual cooperative particle swarms to transfer the arbitrary final time optimal problem into a two-point boundary value problem with fixed terminal one. Then, a fast one-dimensional interval search method is provided in each selected interval to reduce the search space of the particles in the iterations. Furthermore, a theorem that determines the range of the parameters of the Legendre polynomial is investigated, and the problem can be solved to get closer to the global optimal solution.Lastly, in numeral simulation, the results demonstrate that the LCPSO algorithm can solve the unsmooth trajectory optimization problem of X-43 effectively and obtain a smooth control variable. And the appropriate range of parameter values will significantly reduce the complexity of the optimization search. Moreover, in the solution to the orbit transfer problem, the comparisons with the existing GA and DE algorithms represent the notion that the proposed LCPSO method has better performance in the speed of convergence, final accuracy, and constraints satisfaction.The flaw of the proposed LCPSO algorithms is that the parameters of the Legendre polynomial are sensitive in adjustment. Nevertheless, in the trajectory optimization problem of HV, the variety of parameters is a relatively gentle process because of the dynamic characters of HV, and the parameters are adjusted less frequently in iterations. Normally, the parameters will remain the same after they are adjusted once before. The optimization between parameter sensitivity and optimal solution will be studied to make the method more practical in engineering in our future work.



---

*Source: 5036791-2018-04-12.xml*