# Multiphase Return Trajectory Optimization Based on Hybrid Algorithm

Authors: Yi Yang; Ying Nan
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A hybrid trajectory optimization method consisting of Gauss pseudospectral method (GPM) and natural computation algorithm has been developed and utilized to solve multiphase return trajectory optimization problem, where a phase is defined as a subinterval in which the right-hand side of the differential equation is continuous. GPM converts the optimal control problem to a nonlinear programming problem (NLP), which helps to improve calculation accuracy and speed of natural computation algorithm. Through numerical simulations, it is found that the multiphase optimal control problem could be solved perfectly.


---

## Body

## 1. Introduction

Return trajectory, as an important part of satellite recovery mission, is flight trajectory of spacecraft from its working orbit to required reentry point. The purpose of return trajectory optimization is to obtain control variables of spacecraft that null the endpoint errors between the spacecraft final state vector and the target state vector while minimizing the total fuel mass consumption.Return trajectory optimization problem, which also could be considered as orbit transfer optimization problem, is an optimal control problem in the functional space in principle. Typically, orbit transfer optimization algorithm could be divided into two types according to their finite dimensional processing method for dynamic equation: direct methods and indirect methods.Indirect methods [1–4] transform the trajectory optimization problem into a two-point boundary value problem based on Pontryagin minimum principle and boundary condition. Indirect methods are of high solution accuracy and satisfy the optimality conditions strictly in theory, which have been widely used in orbit transfer optimization problem, especially in low-thrust, long-duration orbit transfer problem. However, indirect methods rely on correct estimation of unknown boundary condition and need good guess of initial value of consistent variable.Direct methods [5–8] usually discretize the optimal control problem into a nonlinear programming problem. Direct methods have broader convergence domain and lower requirements for the accuracy of initial value evaluation. However, compared to indirect methods, direct methods have lower accuracy and are easy to converge to local optimal solution, which indicates that global optimization algorithm is the main direction of direct methods.Pseudospectral method, as a well-known direct method, has been widely used to solve optimal control problem [9–12]. GPM transforms the original optimal control problem into a nonlinear programming problem (NLP) at the discrete points subject to continuous state equality and inequality constraints. Therefore, the dynamic equations can be transformed into the algebraic equations and the cost function can also be expressed as the function of the states and control. The NLP is usually solved by sequential quadratic programming method (SQP), but SQP easily falls into local optimal solution and is highly dependent on initial values, which make it hard to obtain global optimal solution.Natural computation algorithm, which was first proposed by Nan Ying in 2013, is a global optimization algorithm; it numerically simulates numerous substance units, combined substances, and reciprocities among the substances, including the organized running and decomposing and compounding of substance units and combined substances. The running purposes of substance units and combined substances are pursuing benefit performance indexes for themselves. Numerical simulation in [13, 14] indicates its global optimization ability, but it has shortcomings of numerous computation amounts and low computation accuracy according to its search principle.Therefore, in this paper, we proposed a hybrid algorithm to combine the advantages of these two methods, where GPM is adopted to transform the original optimal control problem into a nonlinear programming problem (NLP), and then solve this NLP based on natural computation algorithm (NCA).


## 2. Problem Formulation

In this section, return trajectory mission of spacecraft and dynamic equation applied will be precisely described.The spacecraft starts from a circular LEO with the altitude of 400 km and the deorbit window is shown in Table1. Spacecraft begins to decelerate under the backward thrust of engine and flies to the reference reentry point accurately under the control of attack angle and slide angle. The reference reentry point parameters are summarized in Table 2.Table 1
Initial orbit elements.
Universal time Semimajor axisa Eccentricitye Inclinationi RAANΩ Argument of perigeeω True anomalyf 2017/10/212:24:36 6778140 (m) 0.003 42° 122.722° 0° 227.298°Table 2
Reference reentry point parameters.
Parameters Reference value Reentry geographic height,h 120 km Reentry speed over ground,v2e 7586.2 m/s Reentry path angle,θ −1.235° Reentry heading angle,ψv 45.71° Reentry longitude,λ 34.00° Reentry latitude,ϕ 2.612°The initial mass of the spacecraft is assumed to be 3000 kg, and the engine which is fixed on the bottom of the spacecraft has stable thrust 1173.5 N and specific impulse 300 s.To minimize the performance index (i.e., fuel consumption):(1)J=∫t0tfm˙dtthe control variable u(t)=[α(t);β(t)], the state x(t)=[v(t);θ(t);ψv(t);h(t);ϕ(t);λ(t)], the initial time t0, and the terminal time tf have to be optimized.The states and control variables of spacecraft are determined by the following dynamic equations:(2)r˙=vsin⁡θ,ϕ˙=vcos⁡θcos⁡ψvr,λ˙=-vcos⁡θsin⁡ψvrcos⁡ϕ,v˙=Tcos⁡αcos⁡βm+grsin⁡θ+gϕcos⁡θcos⁡ψv-we2rcos⁡ϕsin⁡ϕcos⁡θcos⁡ψv-cos⁡ϕsin⁡θ,θ˙=Tsin⁡αmv+grcos⁡θ-gϕsin⁡λcos⁡ψvv+vcos⁡θr-2wecos⁡ϕsin⁡ψv+we2rvcos⁡ϕcos⁡θcos⁡ϕ+sin⁡ϕsin⁡θcos⁡ψv,-ψ˙v=Tcos⁡αsin⁡βmvcos⁡θ+gϕsin⁡ψvvcos⁡θ-vsin⁡ψvrtan⁡ϕcos⁡θ-2wecos⁡ϕtan⁡θcos⁡ψv-sin⁡ϕ-we2rvcos⁡θsin⁡ϕcos⁡ϕsin⁡ψv,where [r;ϕ;λ] denote the distance from earth’s core, altitude, and longitude of the spacecraft, respectively; [v;θ;ψv] denote the velocity, path angle, and heading angle of the spacecraft, respectively; [α;β] are attack angle and slide angle of the spacecraft; T is the thrust, which is fixed on 1173.5 N; m denotes the mass of spacecraft; we is the rotational-angular velocity of the earth; and [gϕ;gr] are the acceleration of gravity component on distance and altitude.Also, the initial and terminal constraints have to be considered, as described in (3) and (4). Moreover, the path constraints should be satisfied as (5):(3)v;θ;ψv;h;ϕ;λt0≡a;e;i;Ω;ω;ft0,(4)v;θ;ψv;h;ϕ;λtf≡vf;θf;ψvf;hf;ϕf;λf,(5)α˙2t+β˙2t≤0.3∘s.


## 3. Description of Hybrid Algorithm

This hybrid algorithm fully takes advantages of GPM and NCA; it converts the original multiphase trajectory optimization problem into a discrete NLP firstly by using Gauss pseudospectral method and then solves the NLP by using NCA. GPM improves accuracy and computation speed of algorithm and NCA can obtain the global optimal solution. GPM is one fitting method which uses global interpolation polynomial, which needs less input parameters and possesses higher accuracy compared to other methods. Compared with common SQP algorithm, NCA will not fall into local optimal solution, and, independent of initial values, numerical simulation in the next section will contrast optimization results of these two methods.


### 3.1. Gauss Pseudospectral Method

GPM is one fitting method which uses global interpolation polynomial, which needs less input parameters and possesses higher accuracy compared to other methods.Since the states of return trajectory (x(t)=[v(t);θ(t);ψv(t);h(t);ϕ(t);λ(t)]) have different orders of magnitude, states with high order of magnitude need to be normalized as follows:(6)h-=hRe,v-=vg0Re,t¯=tRe/g0,m-=mm0,T-=Tm0g0,where m0 is the initial mass of the spacecraft and Re is the radius of the earth.Gauss pseudospectral method for solving multiphase optimal control problem is described as follows.Firstly, transform the independent variablet to the variable τ∈[-1,1], for every phase of the problem p∈[1,…,P], and ·(k) denotes information in the kth phase:(7)tp=tfp-t0p2τp+tfp+t0p2.Next, choose the Legendre-Gauss (LG) points as the collection points in each phase, which are the roots of theNth-degree Legendre polynomial PN(τ) as follows:(8)PN=12NN!dNdτNτ2-1N.The whole discrete points are LG points and 2 points (τ0=-1 and τN+1=1).The state of spacecraft trajectory is discretized atN+1 Legendre-Gauss points:(9)xpτ≈Xpτ=∑i=0NpLipτXpτip=∑i=0NpLipτXip,where Li(p)(τ) is defined as(10)Lipτ=∏j=0,j≠iNpτ-τjpτip-τjp.Differentials of state variables with respect to τ in (2) in each phase are(11)x˙pτ≈X˙pτ=∑i=0NpL˙ipτXip,p=1,2,…,P.The derivative of each Lagrange polynomial at the LG points can be represented compactly in the form of a differentiation matrix D∈RN×(N+1) as(12)Dkip=L˙ipτk=∑l=0N∏j=0,j≠i,lNτk-τj∏j=0,j≠iNτi-τj,where k=1,…,N and i=1,…,N.The dynamic constraint of (2) is converted into algebraic constraints in each phase p=[1,…,P] of the problem as(13)∑i=0NDkiXi-tf-t02fXk,Uk,τk;q,t0,tf=0,where Xk≡X(τk)∈Rn and Uk≡U(τk)∈Rm(k=1,…,N). Next, define additional variable X0 and Xf as(14)X0=Xτ0,XN+1=X0+tf-t02∑k=1NωkfXk,Uk,τk;q,t0,tf.Substituting (11) into (13), (13) will become a linear equation:(15)XN+1-X0-∑i=0N∑k=1NωkDkiXi=0.Similarly, the path constraints of (4) are discretized at LG points as(16)Cmin≤CXk,Uk,τk;t0,q,tf≤Cmax.Besides, the boundary conditions of (3) are given by(17)ϕmin≤ϕX0,t0,XN+1,tf≤ϕmax.Finally, the linking constraints are mapped using the values at the terminal and start points, respectively, of the phase pairs(pl,pr)∈[1,…,P](l,r=[1,…,P]) as (18)Lmins≤LsXN+1pls,tfpls;qpls,X0pus,t0pus;qpus≤Lmaxs.As the discretization process discussed above, the original multiphase continuous optimal control problem is transformed into a discrete NLP successfully, which could be solved by natural computation algorithm (NCA) later.


### 3.2. Natural Computation Algorithm (NCA)

Natural computation algorithm is a global optimization method, which simulates running rules of all radical element units (REUs) and their combinations. In NCA, all units and combinations are pursuing for optimal performance index. Original NCA adopts average discrete method, which may lead to some shortcomings such as low accuracy and long computation time because equidistance interpolating causes fluctuations near initial and terminal points. To overcome this issue, GPM is used here to convert continuous trajectory optimization problem to NLP, which could be solved by NCA with higher accuracy and computation speed.Computational steps of NCA are described as follows.Initial SettingSetting  1. Structure a Ns=4-dimensional physical space, which consists of time dimension and 3-space dimension. Define the sizes of 3-space dimension: iSX_max, iSY_max, and iSZ_max. Dimension of control and strategy’s storage space is iCl_max,l=1,2,…,lmax, where lmax is the number of control variables.Setting  2. Transform the trajectory optimization problem into NLP by using GPM, which was described above in Section 3.1. The number of calculation steps is N, which is the sum of each phase’s discrete points (N=∑p=1PNp, where NP equals degree of Legendre polynomial in pth phase).The detailed calculation methods for each step of natural computation are described as follows.Step  1. Calculate the running states from time τϒ,min to τϒ+1,min(τϒ+1,min=τϒ,max), where ϒ=1. Eliminate the combinations with poor development according to the evolution and elimination law (discussed in below). The first level’s optimal control U→∗(1)=[u→∗(1),μ→∗(1)] would be obtained as follows.Retain the suitable combinations on the basis of (23) which is posed below and also the optimal control:(19)u→∗X→∗1,J1∗X→∗1,1.Save the states of REUs and combinations in storagePsi[1][iSX][iSY][iSZ][0], save the performance index in storage Psi[1][iSX][iSY][iSZ][1], save the interrelation function in storage Psi[1][iSX][iSY][iSZ][2], save the optimal control in storage Usi[1][lC,max], and save the history information iCτ_Opt in storage †i[1][lC,max].Stepκ. Calculate the running states from time τϒ,min to τϒ+1,min(τϒ+1,min=τϒ,max), where ϒ=κ. Eliminate the combinations with poor development according to evolution and elimination law. The κth level’s optimal control U→∗(κ)=[u→∗(κ),μ→∗(κ)] would be obtained as follows:(20)Jκ∗X→∗κ,κ=minu→k,iκ-1→u→k,iκLk,i,jX→κ,u→k,iκ,κ+Jκ-1∗x→∗κ-1,κ-1,x→k,iκ=f→k,iX→κ-1,u→k,iκ-1,κ-1.Retain the suitable combinations on the basis of (23) and also the optimal control:(21)u→∗X→∗1,J1∗X→∗1,1.Save the states of REUs and combinations in storagePsi[κ][iSX][iSY][iSZ][0], save the performance index in storage Psi[κ][iSX][iSY][iSZ][1], save the interrelation function in storage Psi[κ][iSX][iSY][iSZ][2], save the optimal control in storage Usi[κ][lC,max], and save the history information iCτ_Opt in storage †i[κ][lC,max].Step  N. Calculate the running states from time τϒ,min to τϒ+1,min(τϒ+1,min=τϒ,max), where ϒ=N. Eliminate the combinations with poor development according to (23). The Nth level’s optimal control U→∗(N-1)=[u→∗(N-1),μ→∗(N-1)] would be obtained as follows:(22)JN∗X→∗N,N=minu→k,iN-1→u→k,iNLk,i,jX→N,u→k,iN,N+JN-1∗x→∗N-1,N-1,x→k,iN=f→k,iX→N-1,u→k,iN-1,N-1.The final states have to satisfy the terminal constraints (17). Similar to Step κ, save the calculation data in storage corresponding to N.After the calculations, all the pieces of information for each step, including the state variables, control variables, and performance indexes, have been saved successfully and could be read freely. The final optimal combination and its history are the optimal control plan.


#### 3.2.1. The Evolution and Elimination Law

(  1) Survival of the Fittest. By the end of each time interval [τϒ,min,τϒ+1,min], numerous combinations’ running states are generated in the physical space (each REU generates Qk,i2 states at the next time point, and each combination generates Qα,δ2 states). In the natural computation algorithm, all REUs and combinations need to follow the law of nature, which means that the survival space is finite and the total number of REUs and combinations cannot exceed the size of physical space. In each space unit, there exists only one fittest REU, and in the physical space, the number of combinations is also restricted.Suppose that℘Agv is the capacity of one space unit; for all ℘ combinations in the physical space, only ƛ combinations can survive:(23)ƛ=0,if℘=01,if1≤℘<℘Agvξ,if℘Agv·ξ≤℘<℘Agv·ξ+℘Agv-1.


## 3.1. Gauss Pseudospectral Method

GPM is one fitting method which uses global interpolation polynomial, which needs less input parameters and possesses higher accuracy compared to other methods.Since the states of return trajectory (x(t)=[v(t);θ(t);ψv(t);h(t);ϕ(t);λ(t)]) have different orders of magnitude, states with high order of magnitude need to be normalized as follows:(6)h-=hRe,v-=vg0Re,t¯=tRe/g0,m-=mm0,T-=Tm0g0,where m0 is the initial mass of the spacecraft and Re is the radius of the earth.Gauss pseudospectral method for solving multiphase optimal control problem is described as follows.Firstly, transform the independent variablet to the variable τ∈[-1,1], for every phase of the problem p∈[1,…,P], and ·(k) denotes information in the kth phase:(7)tp=tfp-t0p2τp+tfp+t0p2.Next, choose the Legendre-Gauss (LG) points as the collection points in each phase, which are the roots of theNth-degree Legendre polynomial PN(τ) as follows:(8)PN=12NN!dNdτNτ2-1N.The whole discrete points are LG points and 2 points (τ0=-1 and τN+1=1).The state of spacecraft trajectory is discretized atN+1 Legendre-Gauss points:(9)xpτ≈Xpτ=∑i=0NpLipτXpτip=∑i=0NpLipτXip,where Li(p)(τ) is defined as(10)Lipτ=∏j=0,j≠iNpτ-τjpτip-τjp.Differentials of state variables with respect to τ in (2) in each phase are(11)x˙pτ≈X˙pτ=∑i=0NpL˙ipτXip,p=1,2,…,P.The derivative of each Lagrange polynomial at the LG points can be represented compactly in the form of a differentiation matrix D∈RN×(N+1) as(12)Dkip=L˙ipτk=∑l=0N∏j=0,j≠i,lNτk-τj∏j=0,j≠iNτi-τj,where k=1,…,N and i=1,…,N.The dynamic constraint of (2) is converted into algebraic constraints in each phase p=[1,…,P] of the problem as(13)∑i=0NDkiXi-tf-t02fXk,Uk,τk;q,t0,tf=0,where Xk≡X(τk)∈Rn and Uk≡U(τk)∈Rm(k=1,…,N). Next, define additional variable X0 and Xf as(14)X0=Xτ0,XN+1=X0+tf-t02∑k=1NωkfXk,Uk,τk;q,t0,tf.Substituting (11) into (13), (13) will become a linear equation:(15)XN+1-X0-∑i=0N∑k=1NωkDkiXi=0.Similarly, the path constraints of (4) are discretized at LG points as(16)Cmin≤CXk,Uk,τk;t0,q,tf≤Cmax.Besides, the boundary conditions of (3) are given by(17)ϕmin≤ϕX0,t0,XN+1,tf≤ϕmax.Finally, the linking constraints are mapped using the values at the terminal and start points, respectively, of the phase pairs(pl,pr)∈[1,…,P](l,r=[1,…,P]) as (18)Lmins≤LsXN+1pls,tfpls;qpls,X0pus,t0pus;qpus≤Lmaxs.As the discretization process discussed above, the original multiphase continuous optimal control problem is transformed into a discrete NLP successfully, which could be solved by natural computation algorithm (NCA) later.


## 3.2. Natural Computation Algorithm (NCA)

Natural computation algorithm is a global optimization method, which simulates running rules of all radical element units (REUs) and their combinations. In NCA, all units and combinations are pursuing for optimal performance index. Original NCA adopts average discrete method, which may lead to some shortcomings such as low accuracy and long computation time because equidistance interpolating causes fluctuations near initial and terminal points. To overcome this issue, GPM is used here to convert continuous trajectory optimization problem to NLP, which could be solved by NCA with higher accuracy and computation speed.Computational steps of NCA are described as follows.Initial SettingSetting  1. Structure a Ns=4-dimensional physical space, which consists of time dimension and 3-space dimension. Define the sizes of 3-space dimension: iSX_max, iSY_max, and iSZ_max. Dimension of control and strategy’s storage space is iCl_max,l=1,2,…,lmax, where lmax is the number of control variables.Setting  2. Transform the trajectory optimization problem into NLP by using GPM, which was described above in Section 3.1. The number of calculation steps is N, which is the sum of each phase’s discrete points (N=∑p=1PNp, where NP equals degree of Legendre polynomial in pth phase).The detailed calculation methods for each step of natural computation are described as follows.Step  1. Calculate the running states from time τϒ,min to τϒ+1,min(τϒ+1,min=τϒ,max), where ϒ=1. Eliminate the combinations with poor development according to the evolution and elimination law (discussed in below). The first level’s optimal control U→∗(1)=[u→∗(1),μ→∗(1)] would be obtained as follows.Retain the suitable combinations on the basis of (23) which is posed below and also the optimal control:(19)u→∗X→∗1,J1∗X→∗1,1.Save the states of REUs and combinations in storagePsi[1][iSX][iSY][iSZ][0], save the performance index in storage Psi[1][iSX][iSY][iSZ][1], save the interrelation function in storage Psi[1][iSX][iSY][iSZ][2], save the optimal control in storage Usi[1][lC,max], and save the history information iCτ_Opt in storage †i[1][lC,max].Stepκ. Calculate the running states from time τϒ,min to τϒ+1,min(τϒ+1,min=τϒ,max), where ϒ=κ. Eliminate the combinations with poor development according to evolution and elimination law. The κth level’s optimal control U→∗(κ)=[u→∗(κ),μ→∗(κ)] would be obtained as follows:(20)Jκ∗X→∗κ,κ=minu→k,iκ-1→u→k,iκLk,i,jX→κ,u→k,iκ,κ+Jκ-1∗x→∗κ-1,κ-1,x→k,iκ=f→k,iX→κ-1,u→k,iκ-1,κ-1.Retain the suitable combinations on the basis of (23) and also the optimal control:(21)u→∗X→∗1,J1∗X→∗1,1.Save the states of REUs and combinations in storagePsi[κ][iSX][iSY][iSZ][0], save the performance index in storage Psi[κ][iSX][iSY][iSZ][1], save the interrelation function in storage Psi[κ][iSX][iSY][iSZ][2], save the optimal control in storage Usi[κ][lC,max], and save the history information iCτ_Opt in storage †i[κ][lC,max].Step  N. Calculate the running states from time τϒ,min to τϒ+1,min(τϒ+1,min=τϒ,max), where ϒ=N. Eliminate the combinations with poor development according to (23). The Nth level’s optimal control U→∗(N-1)=[u→∗(N-1),μ→∗(N-1)] would be obtained as follows:(22)JN∗X→∗N,N=minu→k,iN-1→u→k,iNLk,i,jX→N,u→k,iN,N+JN-1∗x→∗N-1,N-1,x→k,iN=f→k,iX→N-1,u→k,iN-1,N-1.The final states have to satisfy the terminal constraints (17). Similar to Step κ, save the calculation data in storage corresponding to N.After the calculations, all the pieces of information for each step, including the state variables, control variables, and performance indexes, have been saved successfully and could be read freely. The final optimal combination and its history are the optimal control plan.


### 3.2.1. The Evolution and Elimination Law

(  1) Survival of the Fittest. By the end of each time interval [τϒ,min,τϒ+1,min], numerous combinations’ running states are generated in the physical space (each REU generates Qk,i2 states at the next time point, and each combination generates Qα,δ2 states). In the natural computation algorithm, all REUs and combinations need to follow the law of nature, which means that the survival space is finite and the total number of REUs and combinations cannot exceed the size of physical space. In each space unit, there exists only one fittest REU, and in the physical space, the number of combinations is also restricted.Suppose that℘Agv is the capacity of one space unit; for all ℘ combinations in the physical space, only ƛ combinations can survive:(23)ƛ=0,if℘=01,if1≤℘<℘Agvξ,if℘Agv·ξ≤℘<℘Agv·ξ+℘Agv-1.


## 3.2.1. The Evolution and Elimination Law

(  1) Survival of the Fittest. By the end of each time interval [τϒ,min,τϒ+1,min], numerous combinations’ running states are generated in the physical space (each REU generates Qk,i2 states at the next time point, and each combination generates Qα,δ2 states). In the natural computation algorithm, all REUs and combinations need to follow the law of nature, which means that the survival space is finite and the total number of REUs and combinations cannot exceed the size of physical space. In each space unit, there exists only one fittest REU, and in the physical space, the number of combinations is also restricted.Suppose that℘Agv is the capacity of one space unit; for all ℘ combinations in the physical space, only ƛ combinations can survive:(23)ƛ=0,if℘=01,if1≤℘<℘Agvξ,if℘Agv·ξ≤℘<℘Agv·ξ+℘Agv-1.


## 4. Numerical Simulation and Analysis

Five phases are assumed for the whole return trajectory, and [OFF, ON, OFF, ON, OFF] are engine work states corresponding to five phases separately. The simulation results of optimal return trajectory are shown in Figures1–12, the red lines are optimal trajectory obtained by the hybrid algorithm, and the black lines are optimal trajectory obtained by GPM and SQP.Figure 1
Time history of velocity of spacecraft in return trajectory.Figure 2
Time history of path angle of spacecraft in return trajectory.Figure 3
Time history of heading angle of spacecraft in return trajectory.Figure 4
Time history of height of spacecraft in return trajectory.Figure 5
Time history of latitude of spacecraft in return trajectory.Figure 6
Time history of longitude of spacecraft in return trajectory.Figure 7
Satellite ground track of the return trajectory.Figure 8
Time history of attack angleα of spacecraft in return trajectory.Figure 9
Time history of slide angleβ trajectory of spacecraft in return trajectory.Figure 10
Time history of mass of spacecraft in return trajectory.Figure 11
Time history of attack angle rateα˙ of spacecraft in return trajectory.Figure 12
Time history of slide angle rateβ˙ of spacecraft in return trajectory.From the simulation results shown in Figures1–12, the hybrid algorithm gets a better return trajectory with lower fuel consumption compared to SQP algorithm, which indicates that NCA has a better ability for global optimization.The optimal return trajectory can be described as follows: at the beginning of the second phase, the attack angle is almost 180°, which means that the thrust has been turned to the same direction as that of the velocity; then the engine works to decline the velocity of spacecraft until the end of second phase; in the third phase, the spacecraft flies without power in a transition track; the attack angle and slide angle are adjusted to optimal values; and, in fourth phase, the engine works to accelerate reaching the reentry point.The final optimal results are as follows:(1)
The maximum final mass is2885.23kg, which means that the performance index (fuel consumption) is 114.77 kg.(2)
The terminal state of flight trajectory isvt;θt;ψvt;ht;ϕt;λt=[7586.20,-1.2351,314.3132,119992.011,2.5951,43.0032].(3)
The state errors between terminal state of optimal trajectory and reference reentry point state are listed in Table3, which satisfy required precision.Table 3
Terminal state errors.
Trajectory state Errors V   ( m / s ) 0.0107501755680642 θ   ( r a d ) - 7.66425969788223 e - 06 ψ v   ( r a d ) - 5.64008839631924 e - 07 h   ( m ) - 8.8399759697844 ϕ   ( r a d ) 3.40221324294371 e - 06 λ   ( r a d ) 4.11743958561583 e - 06


## 5. Conclusion

We have proposed a novel hybrid algorithm, consisting of Gauss pseudospectral method and natural computation algorithm, to solve the multiphase optimal control problem. First, Gauss pseudospectral method (GPM) is used to discretize the original optimal control problem to NLP; and then, natural computation algorithm is used to search the optimal solution of NLP. GPM is characterized as high-grade precision fitting method while NCA is capable of searching the optimal solution globally. We take advantages of both GPM and NCA in this hybrid algorithm. To verify the effectiveness, a multiphase return trajectory design has been illustrated, where the fuel consumption is minimized; meanwhile all the terminal constraints and path constraints are satisfied.As a new global optimization algorithm, NCA has some pros and cons. The main drawbacks are low fitting accuracy and large computation amount. The newly proposed hybrid algorithm uses GPM to replace the average discrete method to discretize the original continuous problem, which improves both calculation accuracy and speed.However, there is still a lot of work that can be done to improve the results. In this paper, we set the whole trajectory as a REU; as a result, the length of the REU is the sum of five phases’ discrete points, which may cause dimension disaster and large calculation amount. The possible solution is to set every phase as a REU and use combination of REUs to represent the whole trajectory. We believe that it will reduce the dimension of calculation greatly.



---

*Source: 6257050-2016-10-24.xml*