Multiphase Return Trajectory Optimization Based on Hybrid Algorithm

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.


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][2][3][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][6][7][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][10][11][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).

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 Table 1.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.
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): the control variable () = [(); ()], the state () = [V(); ();  V (); ℎ(); (); ()], the initial time  0 , and the terminal time   have to be optimized.The states and control variables of spacecraft are determined by the following dynamic equations: where [; ; ] 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;  is the thrust, which is fixed on 1173.5 N;  denotes the mass of spacecraft;   is the rotational-angular velocity of the earth; and [  ;   ] 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):

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.

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 (() = [V(); ();  V (); ℎ(); (); ()]) have different orders of magnitude, states with high order of magnitude need to be normalized as follows: where  0 is the initial mass of the spacecraft and   is the radius of the earth.Gauss pseudospectral method for solving multiphase optimal control problem is described as follows.
Firstly, transform the independent variable  to the variable  ∈ [−1, 1], for every phase of the problem  ∈ [1, . . ., ], and (⋅) () denotes information in the th phase: Next, choose the Legendre-Gauss (LG) points as the collection points in each phase, which are the roots of the th-degree Legendre polynomial   () as follows: The whole discrete points are LG points and 2 points ( 0 = −1 and  +1 = 1).
The state of spacecraft trajectory is discretized at  + 1 Legendre-Gauss points: where  ()  () is defined as Differentials of state variables with respect to  in (2) in each phase are The derivative of each Lagrange polynomial at the LG points can be represented compactly in the form of a differentiation matrix  ∈ R ×(+1) as where  = 1, . . .,  and  = 1, . . ., .

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 Setting
Setting 1. Structure a   = 4-dimensional physical space, which consists of time dimension and 3-space dimension.Define the sizes of 3-space dimension:   max ,   max , and   max .Dimension of control and strategy's storage space is   max ,  = 1, 2, . . .,  max , where  max 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 , which is the sum of each phase's discrete points ( = ∑  =1   , where   equals degree of Legendre polynomial in th phase).
The detailed calculation methods for each step of natural computation are described as follows.
Retain the suitable combinations on the basis of (23) which is posed below and also the optimal control: Save the states of REUs and combinations in storage , save the optimal control in storage   [1][ ,max ], and save the history information   Opt in storage †  [1][ ,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 ⃗  * () = [ ⃗  * (), ⃗  * ()] would be obtained as follows: Retain the suitable combinations on the basis of ( 23) and also the optimal control: Save the states of REUs and combinations in storage The final states have to satisfy the terminal constraints (17).Similar to Step , save the calculation data in storage corresponding to .
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.

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  2 , states at the next time point, and each combination generates  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:

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 Figures 1-12, the red lines are optimal trajectory obtained by the hybrid algorithm, and the black lines are optimal trajectory obtained by GPM and SQP.
From the simulation results shown in Figures 1-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 is 2885.23 kg, which means that the performance index (fuel consumption) is 114.77kg.
( (3) The state errors between terminal state of optimal trajectory and reference reentry point state are listed in Table 3, which satisfy required precision.

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.

Figure 1 :
Figure 1: Time history of velocity of spacecraft in return trajectory.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Time history of path angle of spacecraft in return trajectory.

Figure 5 :
Figure 5: Time history of latitude of spacecraft in return trajectory.

Figure 6 :
Figure 6: Time history of longitude of spacecraft in return trajectory.

Figure 7 :
Figure 7: Satellite ground track of the return trajectory.

Figure 8 :
Figure 8: Time history of attack angle  of spacecraft in return trajectory.

Figure 9 :
Figure 9: Time history of slide angle  trajectory of spacecraft in return trajectory.

Figure 10 :Figure 11 :
Figure 10: Time history of mass of spacecraft in return trajectory.

Figure 12 :
Figure 12: Time history of slide angle rate β of spacecraft in return trajectory.