# Efficient Convex Optimization of Reentry Trajectory via the Chebyshev Pseudospectral Method

Authors: Chun-Mei Yu; Dang-Jun Zhao; Ye Yang
Journal: International Journal of Aerospace Engineering
 (2019)

---

## 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 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 time-consuming 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.


---

## Body

## 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 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 well-known 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–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 time-consuming 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–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.


## 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].
(1)ẋ=vcosθ,ẏ=vsinθ,ḣ=vγ,v̇=−Bv2e−βReh1+λ22E∗,γ̇=Bve−βRehλcosσ−1/v+v,θ̇=Bve−βRehλsinσ.Here, the independent variable is dimensionless timet which is normalized by Re/ge (Re is the earth’s radius, ge is the gravitational acceleration at sea level), states x and y are the horizontal positions normalized by the earth’s radius Re, 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=ρ0ReSrefCL∗/2m, in which ρ0 is the atmospheric density at sea level, Sref is the aerodynamic reference area, m is the vehicle mass, and CL∗ andCD∗, respectively, are the coefficients of lift and drag which produce the maximum lift-to-drag ratio E∗=CL∗/CD∗ for CAV. The control variables in (1) are the bank angle σ and the normalized coefficient of lift λ=CL/CL∗ where CL is the lift coefficient of the vehicle. Thus, the coefficient of lift and drag can be represented by CL=λCL∗ and CD=CL∗1+λ2/2E∗, 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
(2)ẋ=fx,u,t,where the state vector x=x,y,h,v,γ,θT and the control vector is u=λ,σT.


### 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, 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.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−Re/Re with the form
(3)ρ=ρ0e−βReh.Then, the normalized path constraints are defined by(4)P1xt,t=k¯Qe−0.5βRehv3−1k¯qe−βRehv2−1k¯ne−βRehv24E∗2λ2+1+λ22−1≤03×1,in which the normalized coefficients are defined by k¯Q=kQρ0geRe3/Q̇max, k¯q=0.5ρ0geRe/qmax, and k¯n=0.5ρ0ReSrefCL∗2/mnmax, where Q̇max (W/m2), qmax (N/m2), and nmax, respectively, are the maximum allowable heating rate, dynamic pressure, and loading factor.


#### 2.2.2. 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,nNFZ NFZs can be specified by their centre CjNFZ=xj,yjT and the corresponding radius RjNFZ for j=1,⋯,nNFZ. Then, NFZs constraints are formulated as
(5)Njxt,CjNFZ,RjNFZj=1,⋯,nNFZ=x−xjy−yjT≥RjNFZ.


#### 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(6)Γxt0,x0=0,Ψxtf,xf=0,where x0 and xf, 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
(7)Cu,t=σ−σmaxλ−2−λ≤0.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.
(8)u1=λcosσ,u2=λsinσ,u3=λ2=u12+u22,u=u1,u2,u3T.Then, the constraints in (7) are reformulated as
(9)C~u,τ=−u3u3−4−u2−u1tanσmaxu2−u1tanσmax≤0,with an auxiliary equality constraint
(10)Ce=u12+u22−u3=0.


### 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(11)P0:min∫0tf1dt,subjectto2,4,5,6,9,10.


## 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].
(1)ẋ=vcosθ,ẏ=vsinθ,ḣ=vγ,v̇=−Bv2e−βReh1+λ22E∗,γ̇=Bve−βRehλcosσ−1/v+v,θ̇=Bve−βRehλsinσ.Here, the independent variable is dimensionless timet which is normalized by Re/ge (Re is the earth’s radius, ge is the gravitational acceleration at sea level), states x and y are the horizontal positions normalized by the earth’s radius Re, 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=ρ0ReSrefCL∗/2m, in which ρ0 is the atmospheric density at sea level, Sref is the aerodynamic reference area, m is the vehicle mass, and CL∗ andCD∗, respectively, are the coefficients of lift and drag which produce the maximum lift-to-drag ratio E∗=CL∗/CD∗ for CAV. The control variables in (1) are the bank angle σ and the normalized coefficient of lift λ=CL/CL∗ where CL is the lift coefficient of the vehicle. Thus, the coefficient of lift and drag can be represented by CL=λCL∗ and CD=CL∗1+λ2/2E∗, 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
(2)ẋ=fx,u,t,where the state vector x=x,y,h,v,γ,θT and the control vector is u=λ,σT.


## 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, 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.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−Re/Re with the form
(3)ρ=ρ0e−βReh.Then, the normalized path constraints are defined by(4)P1xt,t=k¯Qe−0.5βRehv3−1k¯qe−βRehv2−1k¯ne−βRehv24E∗2λ2+1+λ22−1≤03×1,in which the normalized coefficients are defined by k¯Q=kQρ0geRe3/Q̇max, k¯q=0.5ρ0geRe/qmax, and k¯n=0.5ρ0ReSrefCL∗2/mnmax, where Q̇max (W/m2), qmax (N/m2), and nmax, respectively, are the maximum allowable heating rate, dynamic pressure, and loading factor.


### 2.2.2. 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,nNFZ NFZs can be specified by their centre CjNFZ=xj,yjT and the corresponding radius RjNFZ for j=1,⋯,nNFZ. Then, NFZs constraints are formulated as
(5)Njxt,CjNFZ,RjNFZj=1,⋯,nNFZ=x−xjy−yjT≥RjNFZ.


### 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(6)Γxt0,x0=0,Ψxtf,xf=0,where x0 and xf, 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
(7)Cu,t=σ−σmaxλ−2−λ≤0.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.
(8)u1=λcosσ,u2=λsinσ,u3=λ2=u12+u22,u=u1,u2,u3T.Then, the constraints in (7) are reformulated as
(9)C~u,τ=−u3u3−4−u2−u1tanσmaxu2−u1tanσmax≤0,with an auxiliary equality constraint
(10)Ce=u12+u22−u3=0.


## 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−Re/Re with the form
(3)ρ=ρ0e−βReh.Then, the normalized path constraints are defined by(4)P1xt,t=k¯Qe−0.5βRehv3−1k¯qe−βRehv2−1k¯ne−βRehv24E∗2λ2+1+λ22−1≤03×1,in which the normalized coefficients are defined by k¯Q=kQρ0geRe3/Q̇max, k¯q=0.5ρ0geRe/qmax, and k¯n=0.5ρ0ReSrefCL∗2/mnmax, where Q̇max (W/m2), qmax (N/m2), and nmax, respectively, are the maximum allowable heating rate, dynamic pressure, and loading factor.


## 2.2.2. 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,nNFZ NFZs can be specified by their centre CjNFZ=xj,yjT and the corresponding radius RjNFZ for j=1,⋯,nNFZ. Then, NFZs constraints are formulated as
(5)Njxt,CjNFZ,RjNFZj=1,⋯,nNFZ=x−xjy−yjT≥RjNFZ.


## 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(6)Γxt0,x0=0,Ψxtf,xf=0,where x0 and xf, 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
(7)Cu,t=σ−σmaxλ−2−λ≤0.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.
(8)u1=λcosσ,u2=λsinσ,u3=λ2=u12+u22,u=u1,u2,u3T.Then, the constraints in (7) are reformulated as
(9)C~u,τ=−u3u3−4−u2−u1tanσmaxu2−u1tanσmax≤0,with an auxiliary equality constraint
(10)Ce=u12+u22−u3=0.


## 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(11)P0:min∫0tf1dt,subjectto2,4,5,6,9,10.


## 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 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 trajectory optimization problem, then the modified CGL points proposed by [12] with the form,
(12)τk=−cosπkN,k=0,⋯,N,are employed in this work. The corresponding ith order Chebyshev polynomial is
(13)Pi=cosicos−1τ,i=0,⋯,N,where τ∈−1,1 can be mapped from the independent variable t by the equation t=tf−t0τ/2+tf+t0/2. If t0=0, thent=tfτ+1/2. Consequently, problem P0 is equivalent to the following problem P1.
(14)P1:mint1,(15)subjecttodXdτ=tf2fx,u,τ1=FX,U,τ,τ∈−1,1,(16)ΓX−1,X0=0,ΨX1,Xf=0,(17)GX,U,τ=P1xτ,τ−Njxτ,CjNFZ,RjNFZj=1,⋯,nNFZ≤0,j=1,⋯,nNFZ,(18)CU,τ=C~u,τ−tf≤0,(19)Ce=0,where the augment state X=xT,tT and the control variable U=uT,tfT 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), X0=x0T,t0T and Xf=xfT,tfT. The constraints in (17) are the general form of the aforementioned nonlinear constraints (4) and (5). The augmented control variable tf is constrained by tf≥0 in (18). The state and control vectors of (15) can be approximated by
(20)Xτ=∑i=0NXτiΦiτ,(21)Uτ=∑i=0NUτiΦiτ,where the basis N-order Lagrange interpolating polynomials Φiτ for i=0,1,⋯,N are given by
(22)Φiτ=−1i+1N2ci1−t2ṖNττ−τi,with
(23)ci=2,i=0,N,1,1≤i≤N−1..Differentiating (20) at each CGL point yields the derivative approximation with the following form
(24)Ẋτk=∑i=0NXτiΦ̇iτk=∑i=0NDkiXτi,where Dki are the entries of the differentiation matrix D∈ℝN+1×N+1and can be obtained by
(25)D≔Dki≔−ckci−1i+kτk−τi,j≠k,τk21−τk2,1≤i=k≤N−1,−2N2+16,i=k=0,2N2+16,i=k=N.Thus, the dynamic constraint of (15) is transcribed into an algebraic constraint as
(26)∑i=0NDkiXτi−FXk,Uk,τk=0,k=0,⋯,N,where Xk≡Xτk and Uk≡Uτk. It is to be noted that if the function FX,U,τ in (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, FX,U,τ 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 byXnτ,Unτ 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 Xnτ,Unτ as
(27)dXdτ=AnXn,Un,τX+BnXn,Un,τU+CnXn,Un,τ,where An and Bn, the Jacobian matrices of F with respect to X and U, are given by
(28)An=∂F/∂XXnτ,Unτ=12000tfcosθ0−tfvsinθ0000tfsinθ0tfvcosθ0000tfγtfv0000tfβReBv2e−βReh1+u32E∗−tfBve−βReh1+u3E∗00000−tfβReBve−βRehu1tfBe−βr0hu1+1/v2+100000−tfβReBve−βRehu2tfBe−βRehu20000000000Xnτ,Unτ,Bn=∂F/∂UXnτ,Unτ=12000vcosθ000vsinθ000vγ00−tfBv2e−βReh2E∗−Bv2e−βReh1+u32E∗tfBve−βReh00Bve−βRehu1−1/v+v0tfBve−βReh0Bve−βRehu20001Xnτ,Unτ.The residue term can be calculated by(29)CnXn,Un,τ=FXn,Un,τ−AnXn,Un,τXn−BnXn,Un,τUn.Similar to references [7, 10], additional constraints of the trust region are denoted by the following element-wise inequalities
(30)X−Xn<εX,U−Un<εU,where εX∈R7 and εU∈R4 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 Xnτ,Unτ.Discretizing (27) on the CGL node points like (26) yields the following convex (linear) algebraic equality constraints:
(31)∑i=0NDkiX=AknXkn,Ukn,τX+BknXkn,Ukn,τU+CknXkn,Ukn,τ,k=0,⋯,N,which will enforce the solution at N+1 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
(32)gpX,U,τ≤0,p=1,⋯,m,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)gpXkn,Ukn+∇gpTXknX−Xkn+∇gpTUknU−Ukn≤0,k=0,1,2,⋯N,where ∇gpXkn and ∇gpTUkn, respectively, are the gradients of gp at Xkn and Ukn. 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]:
(34)u12+u22≤u3,which is shown by the blue cone including its surface on the right side of Figure 1. Auxiliary inputs u1 and u2 are limited by u3 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.Figure 1
[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).


### 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 forn=0,1,2⋯:
(35)P2:mint1,(36)subjectto∑i=0NDkiX=AknXkn,Ukn,τX+BknXkn,Ukn,τU+CknXkn,Ukn,τ,k=0,⋯,N,n=1,2,⋯,(37)ΓX−1,X0=0,ΨX1,Xf=0,(38)gpXkn,Ukn+∇gpTXknX−Xkn+∇gpTUknU−Ukn≤0,k=0,1,2,⋯N,(39)u12+u22≤u3,(40)0≤u3≤4,(41)−tanσmaxu1≤u2≤tanσmaxu1,(42)X−Xn<εX,U−Un<εU.


### 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 n=0, and choose an initial control profile U0τkk=0,1,2,N on the CGL node points. An easy choice is set λk=1 and bank angle σk=0. Then driven by U0τk, numerical integration is conducted according to the dynamical model in (1); we henceforth have the initial state profile X0τk.Step 2. At the n+1th iteration, problem P2 is set up by usingXnτk,Unτk, meanwhile point-wisely 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 byXn+1τk,Un+1τk.Step 3. Check whether the convergence condition as follows is satisfied:
(43)maxXn+1−Xn≤δ,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 Xn+1τk,Un+1τk. Stop.


## 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 trajectory optimization problem, then the modified CGL points proposed by [12] with the form,
(12)τk=−cosπkN,k=0,⋯,N,are employed in this work. The corresponding ith order Chebyshev polynomial is
(13)Pi=cosicos−1τ,i=0,⋯,N,where τ∈−1,1 can be mapped from the independent variable t by the equation t=tf−t0τ/2+tf+t0/2. If t0=0, thent=tfτ+1/2. Consequently, problem P0 is equivalent to the following problem P1.
(14)P1:mint1,(15)subjecttodXdτ=tf2fx,u,τ1=FX,U,τ,τ∈−1,1,(16)ΓX−1,X0=0,ΨX1,Xf=0,(17)GX,U,τ=P1xτ,τ−Njxτ,CjNFZ,RjNFZj=1,⋯,nNFZ≤0,j=1,⋯,nNFZ,(18)CU,τ=C~u,τ−tf≤0,(19)Ce=0,where the augment state X=xT,tT and the control variable U=uT,tfT 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), X0=x0T,t0T and Xf=xfT,tfT. The constraints in (17) are the general form of the aforementioned nonlinear constraints (4) and (5). The augmented control variable tf is constrained by tf≥0 in (18). The state and control vectors of (15) can be approximated by
(20)Xτ=∑i=0NXτiΦiτ,(21)Uτ=∑i=0NUτiΦiτ,where the basis N-order Lagrange interpolating polynomials Φiτ for i=0,1,⋯,N are given by
(22)Φiτ=−1i+1N2ci1−t2ṖNττ−τi,with
(23)ci=2,i=0,N,1,1≤i≤N−1..Differentiating (20) at each CGL point yields the derivative approximation with the following form
(24)Ẋτk=∑i=0NXτiΦ̇iτk=∑i=0NDkiXτi,where Dki are the entries of the differentiation matrix D∈ℝN+1×N+1and can be obtained by
(25)D≔Dki≔−ckci−1i+kτk−τi,j≠k,τk21−τk2,1≤i=k≤N−1,−2N2+16,i=k=0,2N2+16,i=k=N.Thus, the dynamic constraint of (15) is transcribed into an algebraic constraint as
(26)∑i=0NDkiXτi−FXk,Uk,τk=0,k=0,⋯,N,where Xk≡Xτk and Uk≡Uτk. It is to be noted that if the function FX,U,τ in (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, FX,U,τ 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 byXnτ,Unτ 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 Xnτ,Unτ as
(27)dXdτ=AnXn,Un,τX+BnXn,Un,τU+CnXn,Un,τ,where An and Bn, the Jacobian matrices of F with respect to X and U, are given by
(28)An=∂F/∂XXnτ,Unτ=12000tfcosθ0−tfvsinθ0000tfsinθ0tfvcosθ0000tfγtfv0000tfβReBv2e−βReh1+u32E∗−tfBve−βReh1+u3E∗00000−tfβReBve−βRehu1tfBe−βr0hu1+1/v2+100000−tfβReBve−βRehu2tfBe−βRehu20000000000Xnτ,Unτ,Bn=∂F/∂UXnτ,Unτ=12000vcosθ000vsinθ000vγ00−tfBv2e−βReh2E∗−Bv2e−βReh1+u32E∗tfBve−βReh00Bve−βRehu1−1/v+v0tfBve−βReh0Bve−βRehu20001Xnτ,Unτ.The residue term can be calculated by(29)CnXn,Un,τ=FXn,Un,τ−AnXn,Un,τXn−BnXn,Un,τUn.Similar to references [7, 10], additional constraints of the trust region are denoted by the following element-wise inequalities
(30)X−Xn<εX,U−Un<εU,where εX∈R7 and εU∈R4 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 Xnτ,Unτ.Discretizing (27) on the CGL node points like (26) yields the following convex (linear) algebraic equality constraints:
(31)∑i=0NDkiX=AknXkn,Ukn,τX+BknXkn,Ukn,τU+CknXkn,Ukn,τ,k=0,⋯,N,which will enforce the solution at N+1 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
(32)gpX,U,τ≤0,p=1,⋯,m,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)gpXkn,Ukn+∇gpTXknX−Xkn+∇gpTUknU−Ukn≤0,k=0,1,2,⋯N,where ∇gpXkn and ∇gpTUkn, respectively, are the gradients of gp at Xkn and Ukn. 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]:
(34)u12+u22≤u3,which is shown by the blue cone including its surface on the right side of Figure 1. Auxiliary inputs u1 and u2 are limited by u3 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.Figure 1
[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).


## 3.2.1. Convexification of Dynamic Constraints

Due to the iterative nature of sequential optimization algorithms, we suppose a reference trajectory denoted byXnτ,Unτ 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 Xnτ,Unτ as
(27)dXdτ=AnXn,Un,τX+BnXn,Un,τU+CnXn,Un,τ,where An and Bn, the Jacobian matrices of F with respect to X and U, are given by
(28)An=∂F/∂XXnτ,Unτ=12000tfcosθ0−tfvsinθ0000tfsinθ0tfvcosθ0000tfγtfv0000tfβReBv2e−βReh1+u32E∗−tfBve−βReh1+u3E∗00000−tfβReBve−βRehu1tfBe−βr0hu1+1/v2+100000−tfβReBve−βRehu2tfBe−βRehu20000000000Xnτ,Unτ,Bn=∂F/∂UXnτ,Unτ=12000vcosθ000vsinθ000vγ00−tfBv2e−βReh2E∗−Bv2e−βReh1+u32E∗tfBve−βReh00Bve−βRehu1−1/v+v0tfBve−βReh0Bve−βRehu20001Xnτ,Unτ.The residue term can be calculated by(29)CnXn,Un,τ=FXn,Un,τ−AnXn,Un,τXn−BnXn,Un,τUn.Similar to references [7, 10], additional constraints of the trust region are denoted by the following element-wise inequalities
(30)X−Xn<εX,U−Un<εU,where εX∈R7 and εU∈R4 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 Xnτ,Unτ.Discretizing (27) on the CGL node points like (26) yields the following convex (linear) algebraic equality constraints:
(31)∑i=0NDkiX=AknXkn,Ukn,τX+BknXkn,Ukn,τU+CknXkn,Ukn,τ,k=0,⋯,N,which will enforce the solution at N+1 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
(32)gpX,U,τ≤0,p=1,⋯,m,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)gpXkn,Ukn+∇gpTXknX−Xkn+∇gpTUknU−Ukn≤0,k=0,1,2,⋯N,where ∇gpXkn and ∇gpTUkn, respectively, are the gradients of gp at Xkn and Ukn. 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]:
(34)u12+u22≤u3,which is shown by the blue cone including its surface on the right side of Figure 1. Auxiliary inputs u1 and u2 are limited by u3 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.Figure 1
[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).


## 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 forn=0,1,2⋯:
(35)P2:mint1,(36)subjectto∑i=0NDkiX=AknXkn,Ukn,τX+BknXkn,Ukn,τU+CknXkn,Ukn,τ,k=0,⋯,N,n=1,2,⋯,(37)ΓX−1,X0=0,ΨX1,Xf=0,(38)gpXkn,Ukn+∇gpTXknX−Xkn+∇gpTUknU−Ukn≤0,k=0,1,2,⋯N,(39)u12+u22≤u3,(40)0≤u3≤4,(41)−tanσmaxu1≤u2≤tanσmaxu1,(42)X−Xn<εX,U−Un<εU.


## 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 n=0, and choose an initial control profile U0τkk=0,1,2,N on the CGL node points. An easy choice is set λk=1 and bank angle σk=0. Then driven by U0τk, numerical integration is conducted according to the dynamical model in (1); we henceforth have the initial state profile X0τk.Step 2. At the n+1th iteration, problem P2 is set up by usingXnτk,Unτk, meanwhile point-wisely 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 byXn+1τk,Un+1τk.Step 3. Check whether the convergence condition as follows is satisfied:
(43)maxXn+1−Xn≤δ,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 Xn+1τk,Un+1τk. Stop.


## 4. 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 Table1, 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=4000kW/m2, qmax=5×105N/m2, and nmax=2.5ge. The boundary conditions are as follows: h0=122km, v0=7315.2m/s, γ0=−1.5∘, θ0=4∘, hf≥20km, vf≥2000m/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
(44)∣X−Xn∣≤0.20.210000Re5000geReπ6π6T,∣U−Un∣≤0.50.510.5T,which is sufficiently large to satisfy Assumption 1 presented by [7]. The stopping criterion δ in (43) is set as δ=0.030.030.0030.040.060.03T.Table 1
CAV mission description.
Latitude (deg)
Longitude (deg)
Radius (km)
Initial point
N 28.5
W -80.6
—
NFZ1
N 20.1
W -3.4
177.7920
NFZ2
N 55.6
E 58.5
277.8000
Target
N 31.2
E 65.7
—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 u1,u2, and u3 to replace the original controls (λ and σ) in (1).The comparison of numerical solutions obtained by Cheb-SCvx, Eq-SCvx, and GPOPS is presented in Table2, 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.Table 2
Comparison of the minimum-time solution for CAV.
Method
Number of discretized points
Iteration number
Optimal index (s)
CPU time (s)
Cheb-SCvx
25
7
3056.4
5.59
30
18
2860.5
14.20
40
8
2862.6
6.75
60
6
2909.5
5.68
80
6
2968.4
12.21
Eq-SCvx
300
7
2869.2
11.15
250
10
2893.8
15.42
GPOPS
-
-
3100.0
20.430The ground track of CAV constrained by NFZs is presented in Figure2, 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).Figure 2
Ground tracks of CAV.Figure 3
Altitude and velocity profile.Figure 4
Flight path angle and heading angle histories.Figure 5
Control variables (normalized lift coefficient and bank angle) histories.Figure 6
Heat flux, dynamic pressure, and normal load histories.


## 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 open-source 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.



---

*Source: 1414279-2019-05-02.xml*