Research Article  Open Access
Efficient Convex Optimization of Reentry Trajectory via the Chebyshev Pseudospectral Method
Abstract
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 nofly zones. The ChebyshevGauss 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.
1. 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 infinitedimensional dynamic optimization problem to a finitedimensional 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 secondorder conic programming (SOCP) algorithm. Further, in [9], they improved their convex optimization algorithm based on the interiorpoint 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–11], the uniformly distributed discretized points are chosen to transform the original infinitedimensional optimization problem into a finitedimensional 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 highdimensional 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 RadauGauss pseudospectral method in [4, 5] and the Chebyshev method in [12–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 ChebyshevGauss 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.
2. 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 which is normalized by ( is the earth’s radius, is the gravitational acceleration at sea level), states and are the horizontal positions normalized by the earth’s radius , states and are the dimensionless altitude and velocity of the vehicle, and and are the flight path angle and heading angle, respectively. The vehiclespecific constant is defined by, in which is the atmospheric density at sea level, is the aerodynamic reference area, is the vehicle mass, and and, respectively, are the coefficients of lift and drag which produce the maximum lifttodrag ratio for CAV. The control variables in (1) are the bank angle and the normalized coefficient of lift where is the lift coefficient of the vehicle. Thus, the coefficient of lift and drag can be represented by and , respectively. 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 and the control vector is .
2.2. 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, nofly 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.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 with the form
Then, the normalized path constraints are defined by in which the normalized coefficients are defined by , , and , where (W/m^{2}), (N/m^{2}), and , respectively, are the maximum allowable heating rate, dynamic pressure, and loading factor.
2.2.2. NoFly Zone Constraints
Nofly zones (NFZs) are considered in this paper which are defined as circular exclusion zones in the horizontal plane with infinite altitude; thus, NFZs can be specified by their centre and the corresponding radius for . 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 and , respectively, are the given initial states and terminal states which should be satisfied in the optimized trajectory [16].
2.2.4. 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
2.3. Reentry Trajectory Optimization Problem
In the general trajectory optimization problem, there are several choices of performance indices to specify different optimization objectives such as maximum range, minimum heat load, and minimum time. In this research, minimum time is chosen as the performance index; thus, the reentry trajectory optimization problem is formulated as
3. 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 infinitedimensional trajectory optimization problem is discretized by the Chebyshev pseudospectral method and henceforth a finitedimensional 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 for , then the node points range from 1 to 1. Since such order is inconvenient for the trajectory optimization problem, then the modified CGL points proposed by [12] with the form, are employed in this work. The corresponding order Chebyshev polynomial is where can be mapped from the independent variable by the equation . If , then. Consequently, problem P0 is equivalent to the following problem P1. where the augment state and the control variable are used in (15) to reformulate the original timefree problem (11) as a timefixed 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), and . The constraints in (17) are the general form of the aforementioned nonlinear constraints (4) and (5). The augmented control variable is constrained by in (18). The state and control vectors of (15) can be approximated by where the basis order Lagrange interpolating polynomials for are given by with
Differentiating (20) at each CGL point yields the derivative approximation with the following form where are the entries of the differentiation matrix and can be obtained by
Thus, the dynamic constraint of (15) is transcribed into an algebraic constraint as where and . It is to be noted that if the function in (15) is a linear function of and , then (26) will be a series of linear equality constraints which are convex in nature. Unfortunately, is a strong nonlinear function; hence, the further convexification technique should be employed in order to formulate an SCvx problem.
3.2. 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.
3.2.1. Convexification of Dynamic Constraints
Due to the iterative nature of sequential optimization algorithms, we suppose a reference trajectory denoted by which is provided by the th iteration solution of the SCvx algorithm, then in the iteration, the dynamic equation (15) can be linearized along the reference trajectory as where and , the Jacobian matrices of with respect to and , 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 elementwise inequalities where and 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 .
Discretizing (27) on the CGL node points like (26) yields the following convex (linear) algebraic equality constraints: which will enforce the solution at CGL node points exactly satisfying the dynamic constraint (15).
3.2.2. Convexification of Path Constraints
Each entry of the concave path constraints (17) can be represented by a generalized inequality as where is the number of path constraints.
In the scenario of trajectory linearization and discretization, the abovementioned concave inequality constraints can be reformulated as where and , respectively, are the gradients of at and . 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.
3.2.3. Convexification of Control Constraints
The equality constraint on the control variables in (19) is obviously concave. We relax the strong equality constraints to a convex inequality constraint as follows [7, 18]: which is shown by the blue cone including its surface on the right side of Figure 1. Auxiliary inputs and are limited by and the maximum bank angle as illustrated by the left image in Figure 1. And the relaxed constraint illustrated by the right image in Figure 1 is obviously convex.
3.3. 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 :
3.4. SCvx Algorithm
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 , and choose an initial control profile on the CGL node points. An easy choice is set and bank angle . Then driven by , numerical integration is conducted according to the dynamical model in (1); we henceforth have the initial state profile .
Step 2. At the ^{th} iteration, problem P2 is set up by using, meanwhile pointwisely checking the violation states of the heat flux, dynamic pressure, load factor, and NFZ constraints defined in (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.
Step 3. Check whether the convergence condition as follows is satisfied: which consists of a series of componentwise inequalities; is a prescribed sufficiently small constant vector. Moreover, the constraints defined in (4) are pointwisely calculated. If all these constraints on the CGL node points and the convergence condition in (43) are satisfied, go to Step 4; otherwise, set and go back to Step 2.
Step 4. The solution for P0 is found to be . Stop.
4. Numerical Results
In this section, the hypersonic gliding vehicle (CAVH) 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 allsolution procedure is executed on a laptop computer with Intel Core i54200 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 , , and . The boundary conditions are as follows: , , , , , , and . The terminal flight heading angle is not constrained. In addition, the bank angle are limited by . The trust region in (42) is given as which is sufficiently large to satisfy Assumption 1 presented by [7]. The stopping criterion in (43) is set as

To verify the proposed SCvx algorithm via Chebyshev pseudospectral method (denoted by ChebSCvx), a comparison study between the proposed SCvx method in [7] (denoted by EqSCvx) as well as GPOPS is conducted. In the EqSCvx method, the equalspace discretized scheme is employed, while GPOPS uses the hpadaptive Radau pseudospectral method to solve optimal control problems. Similar to the proposed method in Section 3, we use the new controls , and to replace the original controls ( and ) in (1).
The comparison of numerical solutions obtained by ChebSCvx, EqSCvx, and GPOPS is presented in Table 2, which reveals that the proposed ChebSCvx method is more efficient and effective than EqSCvx and GPOPS while dealing with the highly constrained entry trajectory optimization problems. It is to be noted that, in the EqSCvx 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 ChebSCvx 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).
5. 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported in part by the National Key R&D Program of China (Grant No. 2018YFC0810202).
References
 A. V. Rao, “A survey of numerical methods for optimal control,” Advances in the Astronautical Sciences, vol. 135, pp. 1–32, 2010. View at: Google Scholar
 L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Interscience, New York, 1962.
 J. T. Betts, “Very lowthrust trajectory optimization using a direct SQP method,” Journal of Computational and Applied Mathematics, vol. 120, no. 12, pp. 27–40, 2000. View at: Publisher Site  Google Scholar
 D. Garg, M. Patterson, C. Darby et al., “Direct trajectory optimization and costate estimation of general optimal control problems using a Radau pseudospectral method,” in Proceedings of the AIAA Guidance, Navigation, and Control Conference and Exhibit, Chicago, Illinois, 2009. View at: Google Scholar
 G. T. Huntington, “Advancement and analysis of a gauss pseudospectral transcription for optimal control problems,” Tech. Rep., Ph.D dissertation, Massachusetts institute of Technology, 2007. View at: Google Scholar
 B. Açıkmeşe and S. R. Ploen, “Convex programming approach to powered descent guidance for Mars landing,” Journal of Guidance, Control, and Dynamics, vol. 30, no. 5, pp. 1353–1366, 2007. View at: Publisher Site  Google Scholar
 D.J. Zhao and Z.Y. Song, “Reentry trajectory optimization with waypoint and nofly zone constraints using multiphase convex programming,” Acta Astronautica, vol. 137, pp. 60–69, 2017. View at: Publisher Site  Google Scholar
 L. Blackmore, B. Acikmese, and D. P. Scharf, “Minimumlandingerror powereddescent guidance for Mars landing using convex optimization,” Journal of Guidance, Control, and Dynamics, vol. 33, no. 4, pp. 1161–1171, 2010. View at: Publisher Site  Google Scholar
 D. Dueri, B. Açıkmeşe, D. P. Scharf, and M. W. Harris, “Customized realtime interiorpoint methods for onboard powereddescent guidance,” Journal of Guidance, Control, and Dynamics, vol. 40, no. 2, pp. 197–212, 2017. View at: Publisher Site  Google Scholar
 X. Liu and P. Lu, “Solving nonconvex optimal control problems by convex optimization,” Journal of Guidance, Control, and Dynamics, vol. 37, no. 3, pp. 750–765, 2014. View at: Publisher Site  Google Scholar
 X. Liu, Z. Shen, and P. Lu, “Entry trajectory optimization by secondorder cone programming,” Journal of Guidance, Control, and Dynamics, vol. 39, no. 2, pp. 227–241, 2016, Articles in advance. View at: Publisher Site  Google Scholar
 F. Fahroo and I. M. Ross, “Direct trajectory optimization by a Chebyshev pseudospectral method,” Journal of Guidance, Control, and Dynamics, vol. 25, no. 1, pp. 160–166, 2002. View at: Publisher Site  Google Scholar
 S. Shahmirzaee Jeshvaghani, A. B. Novinzaddeh, and F. Pazooki, “Multiple stage satellite launch vehicle ascent optimization using Chebyshev wavelets,” Aerospace Science and Technology, vol. 46, pp. 321–330, 2015. View at: Publisher Site  Google Scholar
 X. Wei, L. Liu, Y. Wang, and Y. Yang, “Reentry trajectory optimization for a hypersonic vehicle based on an improved adaptive fireworks algorithm,” International Journal of Aerospace Engineering, vol. 2018, Article ID 8793908, 17 pages, 2018. View at: Publisher Site  Google Scholar
 G. N. Elnagar and M. A. Kazemi, “Pseudospectral Chebyshev optimal control of constrained nonlinear dynamical systems,” Computational Optimization and Applications, vol. 11, no. 2, pp. 195–217, 1998. View at: Publisher Site  Google Scholar
 T. R. Jorris and R. G. Cobb, “Threedimensional trajectory optimization satisfying waypoint and nofly zone constraints,” Journal of Guidance, Control, and Dynamics, vol. 32, no. 2, pp. 551–572, 2009. View at: Publisher Site  Google Scholar
 T. R. Jorris, “Common aero vehicle autonomous reentry trajectory optimization satisfying waypoint and nofly zone constraints,” Tech. Rep., Doctor, Engineering and Management, Air University, Ohio, US, 2007. View at: Google Scholar
 X. Liu et al., “Exact convex relaxation for optimal flight of aerodynamically controlled missiles,” IEEE Transactions on Aerospace and Electronic Systems, vol. 52, no. 4, pp. 1881–1892, 2016. View at: Publisher Site  Google Scholar
 J. Löfberg, “YALMIP : a toolbox for modeling and optimization in MATLAB,” in Proceedings of the CACSD Conference, Taipei, Taiwan, 2004. View at: Google Scholar
 A. Domahidi, E. Chu, and S. Boyd, “ECOS: an SOCP solver for embedded systems,” in European Control Conference, Zurich, Switzerland, 2013. View at: Google Scholar
Copyright
Copyright © 2019 ChunMei Yu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.