Efficient Convex Optimization of Reentry Trajectory via the Chebyshev Pseudospectral Method

A novel sequential convex (SCvx) optimization scheme via the Chebyshev pseudospectral method is proposed for efficiently solving the hypersonic reentry trajectory optimization problem which is highly constrained by heat flux, dynamic pressure, normal load, and multiple no-fly zones. The Chebyshev-Gauss Legend (CGL) node points are used to transcribe the original dynamic constraint into algebraic equality constraint; therefore, a nonlinear programming (NLP) problem is concave and timeconsuming to be solved. The iterative linearization and convexification techniques are proposed to convert the concave constraints into convex constraints; therefore, a sequential convex programming problem can be efficiently solved by convex algorithms. Numerical results and a comparison study reveal that the proposed method is efficient and effective to solve the problem of reentry trajectory optimization with multiple constraints.


Introduction
The trajectory optimization problem for a hypersonic vehicle constrained by heat rating, dynamic pressure, normal load, and other constraints related to the specified mission is often a highly constrained nonlinear dynamic programming problem which, in general, can be solved by two types of methods: direct and indirect methods [1].Indirect methods rely on solving the necessary conditions which is analytically derived from the Pontryagin minimum principles [2].On the contrary, indirect methods, the analytical necessary conditions do not have to be derived, while the parameterization techniques are used to convert the original infinite-dimensional dynamic optimization problem to a finite-dimensional nonlinear programming problem (NLP) which can be solved by some nonlinear programming algorithms such as the wellknown sequential quadratic programming (SQP) [1,3].There are some software packages such as GPOPS [4] and GPOCS [5] based on direct methods for addressing the trajectory optimization problems.However, these mentioned nonlinear programming algorithms cannot provide an a priori guarantee on the convergence and acquisition of a global optimal solution [6].
In recent years, convex optimization methods have been introduced to solve complex trajectory optimization problems because of their unique theoretical advantages: (1) rapid convergent rate and (2) insensitivity to the initial guess solution [7].In [6,8], Açıkmeşe and his coauthors proposed a lossless convexification method for solving the soft landing problem in the Mars exploration, then the highly constrained nonlinear dynamic programming problem is converted into its convex version which is efficiently solved by the second-order conic programming (SOCP) algorithm.Further, in [9], they improved their convex optimization algorithm based on the interior-point methods and thereby an efficient online algorithm for the guidance of soft landing.In [10], an SCvx optimization framework is proposed for solving nonconvex optimal control problems, in which the concave inequality constraint is successively approximated by linearization on the iterated solution rendering a convex optimization problem suited to be solved by SOCP algorithms.This convex optimization method has been successfully used for addressing trajectory optimization problems of hypersonic vehicles [11].
Actually, as for the aforementioned SCvx optimization methods, the appropriate techniques of linearization and discretization are the key factors ensuring that the solution of the convexified problem is still the solution of the original problem.Hence, additional constraints of the trust region are applied to guarantee the validity of linearization in [10,11].Our previous study [7] reveals that the trust region should be carefully selected: if the trust region is given too large, the conditions of linearization will be violated and therefore the solution of the convexified problem is not that of the original problem; on the contrary, a smaller trust region will result in larger iterations and make the convergent rate slower.
In the procedure of discretization mentioned in [6][7][8][9][10][11], the uniformly distributed discretized points are chosen to transform the original infinite-dimensional optimization problem into a finite-dimensional parameter optimization problem.With such numerical discretization scheme, the discretized interval should be sufficiently small in order to obtain a sufficiently accurate solution, while this leads a high-dimensional transformed problem which is timeconsuming to be solved.Actually, the pseudospectral methods have been widely used to solve the trajectory optimization problems, such as the Radau-Gauss pseudospectral method in [4,5] and the Chebyshev method in [12][13][14].The Chebyshev pseudospectral method is a special case of the general spectral methods, in which the functions can be expanded in terms of interpolating polynomials and the expansion coefficients are the values of the function at the Gauss quadrature node points, thereby having the best accuracy in interpolation of a function [12].It has been shown that interpolation at the Chebyshev-Gauss Legend (CGL) nodes presents a closest result to the optimal polynomial approximation to a given function [15].
In this study, we develop a new SCvx optimization algorithm based on the Chebyshev pseudospectral method to improve the SCvx optimization method proposed in our previous work [7], in which the dynamic programming problem of reentry trajectory optimization is transcribed into a nonlinear programming problem by using the equispace discretizing technique, and then the convexification method and SCvx algorithm are employed to efficiently obtain the optimal trajectory.However, in order to obtain a sufficient well solution, the discretization points in the manner of equispace must be sufficiently large due to the requirement of discretizing accuracy.Based on the advantages of the Chebyshev psuedospectral method, CGL node points are used to discretize and approximate the state and control variables, therefore having more accurate approximations and less discretized points than using equispace discretized points.The less discretized points mean that less decision variables in the transformed SCvx problem and, consequently, the computational cost will be dramatically reduced, and this is validated by numerical study results in Section 4.

Problem Formulation
The reentry trajectory optimization problem, including the reentry motion model, various constraints, and the performance index, is formulated in this section.Further, the corresponding convexification techniques are particularly demonstrated.
2.1.Reentry Dynamics of CAV.For simplicity, the nonrotating earth with constant gravitational acceleration is applied during modelling the motion of CAV; meanwhile, it assumes that CAV's motion is with small flight path angle and limited control input (bank angle and angle of attack).Then, the nondimensional equation of CAV is given as [16].
Here, the independent variable is dimensionless time t which is normalized by R e /g e (R e is the earth's radius, g e is the gravitational acceleration at sea level), states x and y are the horizontal positions normalized by the earth's radius R e , states h and v are the dimensionless altitude and velocity of the vehicle, and γ and θ are the flight path angle and heading angle, respectively.The vehicle-specific constant B is defined byB = ρ 0 R e S ref C * L / 2m , in which ρ 0 is the atmospheric density at sea level, S ref is the aerodynamic reference area, m is the vehicle mass, and C * L andC * D , respectively, are the coefficients of lift and drag which produce the maximum lift-to-drag ratio E * = C * L /C * D for CAV.The control variables in (1) are the bank angle σ and the normalized coefficient of lift λ = C L /C * L where C L is the lift coefficient of the vehicle.Thus, the coefficient of lift and drag can be represented by The further details of the aerodynamics of CAV can be found in the literature [17].For convenience, we rewrite (1) as the following general nonlinear system where the state vector x = x, y, h, v, γ, θ T and the control vector is u = λ, σ T .

Flight Constraints.
During the entry flight of CAV, the strong path constraints, such as heating flux, dynamic pressure, and load factor, should be satisfied for guaranteeing the safety of the vehicle.Moreover, in many cases, to avoid the enemy's detection and interception, no-fly zone constraints should be considered in the trajectory planning such that the vehicle keeps away from the area with the deployment of missile defense systems.
2 International Journal of Aerospace Engineering 2.2.1.Path Constraints.In this research, the path constraints including the maximal heating rate, dynamic pressure, and load factor are considered for the safe flight of the vehicle.
According to [7], the dimensional atmospheric density ρ is reasonably assumed be an exponential function of the nondimensional altitude h = R − R e /R e with the form ρ = ρ 0 e −βR e h 3 Then, the normalized path constraints are defined by in which the normalized coefficients are defined by 3 /Q max , k q = 0 5ρ 0 g e R e /q max , and / mn max , where Q max (W/m 2 ), q max (N/m 2 ), and n max , respectively, are the maximum allowable heating rate, dynamic pressure, and loading factor.

No-Fly Zone Constraints.
No-fly zones (NFZs) are considered in this paper which are defined as circular exclusion zones in the horizontal plane with infinite altitude; thus, n NFZ NFZs can be specified by their centre C NFZ j = x j, y j T and the corresponding radius R NFZ j for j = 1, ⋯, n NFZ .Then, NFZs constraints are formulated as 2.2.3.Boundary Constraints.According to the mission profile of CAV, the vehicle's entry start and the target point are prescribed; therefore, the following boundary constraints are where x 0 and x f , respectively, are the given initial states and terminal states which should be satisfied in the optimized trajectory [16].

Control Constraints.
During the entry of CAV, due to the vehicle's aerodynamic characteristics, the normalized lift coefficient λ is confined in a certain range (for example [0,2]).Further, the bank angle σ is bounded to guarantee the stability.Then, the inequality constraints imposed on the controls are formulated as It is noteworthy that if we directly use σ and λ as control variables in the followed SCvx algorithm, the sinusoidal function of σ in (1) will result in the chattering phenomenon in the solution.The detailed reason can be found in [11].Consequently, in accordance with the treatment proposed in [18], three new control variables are defined as follows to replace the original control variables in order to conveniently convexify the control constraints.
Then, the constraints in (7) are reformulated as with an auxiliary equality constraint

Sequential Convex Optimization Based on the Chebyshev Pseudospectral Method
In this section, the Chebyshev pseudospectral method is introduced to reformulate the original dynamic optimization problem (11) as an SCvx programming problem in order to be efficiently solved by the convex optimization algorithms.First, the infinite-dimensional trajectory optimization problem is discretized by the Chebyshev pseudospectral method and henceforth a finite-dimensional parameter optimization problem.Then, the sequential convex optimization problem is formulated to be solved by the SOCP algorithm.
3.1.Chebyshev Pseudospectral Method.In the classical Chebyshev pseudospectral method, the CGL points are given by τ k = cos πk/N for k = 0, ⋯, N, then the node points range from 1 to -1.Since such order is inconvenient for the 3 International Journal of Aerospace Engineering trajectory optimization problem, then the modified CGL points proposed by [12] with the form, are employed in this work.The corresponding i th order Chebyshev polynomial is where τ ∈ −1, 1 can be mapped from the independent variable t by the equation Consequently, problem P0 is equivalent to the following problem P1.
where the augment state X = x T , t T and the control variable U = u T , t f T are used in (15) to reformulate the original time-free problem (11) as a time-fixed problem for the convenience of the following discretized manipulations; meanwhile, Lagrange's problem P0 is converted into Mayer's problem, which is suited to be solved by convex programming since its performance index is linear in nature.In (16), X 0 = x T 0 , t 0 T and X f = x T f , t f T .The constraints in (17) are the general form of the aforementioned nonlinear constraints (4) and (5).The augmented control variable t f is constrained by t f ≥ 0 in (18).The state and control vectors of (15) can be approximated by where the basis N-order Lagrange interpolating polynomials Φ i τ for i = 0, 1, ⋯, N are given by with Differentiating (20) at each CGL point yields the derivative approximation with the following form where D ki are the entries of the differentiation matrix D ∈ ℝ N+1 × N+1 and can be obtained by Thus, the dynamic constraint of ( 15) is transcribed into an algebraic constraint as where 15) is a linear function of X and U, then (26) will be a series of linear equality constraints which are convex in nature.Unfortunately, F X, U, τ is a strong nonlinear function; hence, the further convexification technique should be employed in order to formulate an SCvx problem.

Convexification.
In the above subsection, the dynamic constraint is transformed into an algebraic constraint via the Chebyshev pseudospectral method, and the trajectory optimization problem can be reformulated as an NLP problem which cannot be directly solved by convex optimization algorithms.Hence, appropriate convexification techniques are required to reformulate problem P1 as an SCvx problem.It is obvious that the constraints in ( 15) and ( 17)-( 19) are concave and should be convexified. 4 International Journal of Aerospace Engineering 3.2.1.Convexification of Dynamic Constraints.Due to the iterative nature of sequential optimization algorithms, we suppose a reference trajectory denoted by X n τ , U n τ which is provided by the n-th iteration solution of the SCvx algorithm, then in the n + 1 −th iteration, the dynamic equation ( 15) can be linearized along the reference trajectory X n τ , U n τ as where A n and B n , the Jacobian matrices of F with respect to X and U, are given by The residue term can be calculated by Similar to references [7,10], additional constraints of the trust region are denoted by the following element-wise inequalities where ε X ∈ R 7 and ε U ∈ R 4 are imposed to make the linearized system sufficiently approximate to the original nonlinear system.Actually, the constraints of the trust region (30) confine the deviated trajectory in a prescribed neighbourhood about the reference trajectory X n τ , U n τ .Discretizing (27) on the CGL node points like (26) yields the following convex (linear) algebraic equality constraints: which will enforce the solution at N + 1 CGL node points exactly satisfying the dynamic constraint (15).

Convexification of Path Constraints.
Each entry of the concave path constraints ( 17) can be represented by a generalized inequality as where m is the number of path constraints.
In the scenario of trajectory linearization and discretization, the abovementioned concave inequality constraints can be reformulated as

33
where ∇g p X n k and ∇g T p U n k , respectively, are the gradients of g p at X n k and U n k .Note that the trust region defined by (30) is necessary in order to guarantee the linearized constraints in (33) reasonably approximate to the original constraints in (17).Lemma 1 of [10] presented a theoretical illustration in which the feasible solution of P1 comes from the linearized inequality constraints.Equation ( 33) is also the feasible solution of the original problem.19) is obviously concave.We relax the strong equality constraints to a convex inequality constraint as follows [7,18]:

Convexification of Control Constraints. The equality constraint on the control variables in (
which is shown by the blue cone including its surface on the right side of Figure 1.Auxiliary inputs u 1 and u 2 are limited by u 3 and the maximum bank angle σ max as illustrated by the left image in Figure 1.And the relaxed constraint illustrated by the right image in Figure 1 is obviously convex.

Sequential Convex Optimization
Problem.After the convexifications of path constraints and the control constraints, the optimal solution of P1 is obtained by solving the following relaxed sequential convex optimal problems defined on the CGL node points for n = 0,1,2 ⋯ : According to the solution procedure of the sequential convex programming method, the SCvx algorithm to find the original problem P0 is given as follows.
Step 1. Set n = 0, and choose an initial control profile U 0 τ k k = 0, 1, 2, N on the CGL node points.An easy choice is set λ k = 1 and bank angle σ k = 0. Then driven by U 0 τ k , numerical integration is conducted according to the dynamical model in (1); we henceforth have the initial state profile X 0 τ k .
Step 2. At the n + 1 th iteration, problem P2 is set up by using X n τ k , U n τ k , meanwhile point-wisely checking [18] Relaxation of control constraints.(Left) original nonconvex set: blue surface (10).(Right) relaxed convex set: blue cone includes its surface (34) Reproduced from Liu et al. [18] (under the Creative Commons Attribution License/public domain).4).If all or some of these constraints are violated, the corresponding convexified constraints in (38) are set up.Then solving problem P2 by the SOCP algorithm yields the solution denoted by X n+1 τ k , U n+1 τ k .
Step 3. Check whether the convergence condition as follows is satisfied: which consists of a series of component-wise inequalities; δ is a prescribed sufficiently small constant vector.Moreover, the constraints defined in (4) are point-wisely calculated.If all these constraints on the CGL node points and the convergence condition in (43) are satisfied, go to Step 4; otherwise, set n = n + 1 and go back to Step 2.
Step 4. The solution for P0 is found to be X n+1 τ k , U n+1 τ k .Stop.

Numerical Results
In this section, the hypersonic gliding vehicle (CAV-H) described by ( 1) is used to validate the proposed algorithm in the above section.The MATLAB modelling tool YALMIP [19] is used to formulate the SOCP problem P2, and a light embedded SOCP solver ECOS [20] is used to obtain the solution.The all-solution procedure is executed on a laptop computer with Intel Core i5-4200 at 2.5 GHz.The CAV mission is described in Table 1, which defines the horizontal positions of the initial and target points, as well as two waypoints during the flight.The radius of the first NFZ is chosen to be much smaller than the turn capability of the CAV, while the second NFZ has a large enough radius.The path constraints imposed on the trajectory are given by Q max = 4000 kW/m 2 , q max = 5 × 10 5 N/m 2 , and n max = 2 5g e .The boundary conditions are as follows: which is sufficiently large to satisfy Assumption 1 presented by [7].The stopping criterion δ in (43) is set as δ = 0 03 0 03 0 003 0 04 0 06 0 03 T To verify the proposed SCvx algorithm via Chebyshev pseudospectral method (denoted by Cheb-SCvx), a comparison study between the proposed SCvx method in [7] (denoted by Eq-SCvx) as well as GPOPS is conducted.In the Eq-SCvx method, the equal-space discretized scheme is employed, while GPOPS uses the hp-adaptive Radau pseudospectral method to solve optimal control problems.Similar to the proposed method in Section 3, we use the new controls u 1 , u 2 , and u 3 to replace the original controls (λ and σ) in (1).
The comparison of numerical solutions obtained by Cheb-SCvx, Eq-SCvx, and GPOPS is presented in Table 2, which reveals that the proposed Cheb-SCvx method is more efficient and effective than Eq-SCvx and GPOPS while dealing with the highly constrained entry trajectory optimization problems.It is to be noted that, in the Eq-SCvx algorithm, when the discretized point number is less than 250, the algorithm cannot converge, and no solution can be obtained.The reason lies in the fact that less discretized points result in a larger distance between two adjacent points; thus, the linearization condition in (27) will be violated.And after multiple numerical experiments, we can ascertain that the appropriate number of CGL node points will extremely reduce the computational time and improve the performance index.The ground track of CAV constrained by NFZs is presented in Figure 2, which reveals that both target and NFZ constraints are satisfied by three methods.However, the proposed method renders lower cost than the other two methods do.Altitude and velocity profiles are given in Figure 3, while the flight path angle and heading angle histories are illustrated in Figure 4.Although the terminal altitude, velocity, and heading angle are not assigned a fixed value but limited in a prescribed range, the overall trend of the solutions provided by the three methods is similar.The optimized control variables' histories are presented in Figure 5, which shows that the bank angles obtained by the Cheb-SCvx method are with a more chattering phenomenon than those obtained by the other methods during the final phase, but the constraints on the bank angle are well satisfied.Further, the constraints on the heat rate, dynamic pressure, and normal load are all satisfied (as shown in Figure 6).

Conclusions
Inspired by the application of convex optimization in aerospace trajectory generation and optimization, a novel sequential convex programming algorithm based on the classical Chebyshev pseudospectral method is proposed to solve a highly constrained entry trajectory optimization problem with free final time.By using sequential linearization, convexification, and discretization on the CGL node points, the original concave optimization problem is approximated by a series of SOCP problems, which are solved by opensource solver ECOS.In this work, we concentrate on the conversion of a nonconvex problem to a convex space, so that the converted problem can be efficiently solved by the SOCP method.An efficient and accurate discretization method based on the Chebyshev interpolating polynomials are proposed to facilitate transcribing the dynamic constraint into an algebraic equality constraint; then, the convexification technique based on the linearization is used to set up the convex version of the original trajectory optimization problem.
The numeric results reveal that the proposed method can dramatically reduce computational cost by appropriately choosing the number of discretized points and will be very competitive to fulfill the onboard optimization in the future.
2000 m/s, and γ f = −4 ∘ .The terminal flight heading angle is not constrained.In addition, the bank angle are limited by σ max = π/ 3. The trust region in (42) is given as|X − X n | ≤ 0

Figure 4 :
Figure 4: Flight path angle and heading angle histories.

Table 2 :
Comparison of the minimum-time solution for CAV.