A Convex Programming Method for Rocket Powered Landing With Angle of Attack Constraint

Real time trajectory planning is vital in rocket precision powered landing guidance. Due to the nonconvex angle of attack (AOA) constraint and other constraints, including nonlinear dynamics and thrust constraint, the powered landing trajectory planning problem is highly nonconvex, which makes it difficult to be solved in real time via existing nonconvex optimization algorithms. As the main contribution of this paper, the AOA constraint is taken into account in the problem. A convex feasible set method is presented to handle it, in which a quadratic concave function is introduced to find a convex feasible subset in the original AOA constraint and the second order term of this function is estimated by Gersgorin disc theorem. For the remaining nonconvexities, the lossless convexification is employed to address the thrust constraint, and the successive linearization is performed to handle the nonlinear dynamics. Thus, a convex optimization problem with AOA constraint for landing trajectory planning is built. The optimal solution to the original problem is obtained by iteratively solving convex problems. Numerical simulations show that the proposed method can find a feasible landing trajectory where the AOA constraint is satisfied and has a better convergence performance compared with the traditional linearization method.


I. INTRODUCTION
Reusable rocket is becoming an important space vehicle due to its low launch cost [1]. One of the key technologies to make rockets practically reusable is called as precision landing [2]. To achieve precision landing under large initial state deviation, the landing trajectory should be replanned in real time when new information is obtained in flight. The fast trajectory planning method, which can meet the real time requirement, is demanded in the future.
Due to nonlinear dynamics and various nonconvex state and control constraints, the landing trajectory planning problem is highly nonconvex, which makes it difficult to be solved in real time via existing nonconvex optimization algorithms. In the study of fast trajectory optimization methods [3]- [7], the methods based on convex optimization have a great advantage in solution efficiency. Theoretically, as long as the problem can be formulated in the form of convex optimization, classical convex algorithms (such as primal-dual interior point algorithms) can be utilized to obtain the optimal solution in polynomial time [8]. Various methods have been The associate editor coordinating the review of this manuscript and approving it for publication was Zhiwei Gao . developed to address the nonconvexity. One popular way is lossless convexification. This concept was first introduced by Acikmese in solving Mars landing problem [9], in which the key idea is to relax nonconvex equality constraints by convex inequality constraints. The lossless convexification method has been applied into powered landing [10]- [12], low-thrust transfers [13], UAV-UGV rendezvous [14], and entry problems [15]. This method highly depends on the inherent characteristic of planning problem and may only handle specific control constraints like the thrust constraint in powered landing [9] and the constraint related to bank angle in Ref. [15]. For most other constraints, there does not exist a general lossless convexification method. Another way is the successive linearization. This method is mainly utilized to deal with nonlinear objective function [16], nonlinear dynamics [17], and concave inequality constraints [18]. Mao et al. [19] and Bonalli et al. [20] proposed a guaranteed successive linearization method. However, for nonconvexnonconcave inequalities, the convex set obtained by the successive linearization method may contain infeasible region that may adversely affect the convergence performance (See the numerical simulations in Sec. V). In addition to the above two methods, Sun et al. [21] utilized reasonable VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ approximation and equivalent conversion method to convexify the nonlinear dynamics and nonconvex energy-optimal objective function respectively. However, it is difficult to eliminate the nonlinearity caused by aerodynamic force via approximation and the lossless convexification of the thrust constraint may not be compatible with the energy -optimal objective function [7]. In Ref. [2], [22], the powered landing trajectory planning problem in the presence of angle of attack (AOA) constraint has been studied, but the problem is described in the velocity coordinate frame, and the AOA is treated as a control input. The motion dynamics equations are highly nonlinear in this frame which may adversely affect the convergence performance of successive linearization algorithm. In the landing site coordinate frame [9], [23], the motion dynamics equations have better linearity. However, to our best knowledge, the literature on addressing the AOA constraint in this frame remains limited. One of the main reasons is that the AOA constraint is a nonconvex-nonconcave inequality constraint which is difficult to be handle by lossless convexification and successive linearization method.
As the main contribution of this paper, the AOA constraint is taken into account. A convex feasible set (CFS) method is proposed to handle it. The CFS method was first proposed to address the obstacle avoidance problem in robot motion planning [24]- [26], in which the key idea is to find a CFS of the original problem. Using this method, we find out a convex feasible subset in the AOA constraint by introducing the quadratic concave function. The second order term of this function is estimated by Gersgorin disc theorem. For the remaining nonconvexities, the lossless convexification is employed to address the thrust constraint, and the successive linearization is performed to handle the nonlinear dynamics. Finally, a convex optimization problem with AOA constraint for landing trajectory planning is built. The optimal solution to the original problem is obtained by iteratively solving convex problems.
The rest of this paper is organized as follows: Sec. II describes the original continuous-time nonconvex rocket powered landing trajectory planning problem. In Sec. III, the method to address the nonconvexity of the original problem is discussed. In Sec. IV, we describe the successive solution procedure. Numerical simulations are conducted to validate the proposed algorithm in Sec. V, and Sec. VI concludes the whole work.

II. ROCKET POWERED LANDING TRAJECTORY PLANNING PROBLEM DESCRIPTION
The landing trajectory planning problem has been described in Ref. [23] and the AOA constraint is added in this paper. The specific form can be described as the following optimal control problem. P1: subject to where Eq. (1) is the objective function, Eqs. Eq. (7) is AOA constraint, η max represents the max AOA limit, The AOA is the angle between the velocity and negative thrust direction. The two-dimensional diagram of rocket landing is shown in Figure 1. Although the AOA constraint (7) (represented by dotted lines) looks like a convex cone, it is nonconvex because its axis v will change. To the best of our knowledge, there are very little researches on the AOA constraint. This constraint will be the biggest challenge in this paper.
In a convex optimization problem, the objective function and the inequality constraints are required to be convex or linear, and the equality constraints are required to be linear. In problem P1, the dynamics Eqs. (3)-(5) are nonlinear caused by R, 1/m, T . The minimum-thrust constraint (the left inequality of Eq.(6) is concave, and the AOA constraint is nonconvex-nonconcave. Therefore, the problem P1 cannot be solved by convex optimization directly. In the next section, we will discuss how to formulate it into a convex optimization problem.

III. CONVEXFICATION FORMULATION A. CONVEX FEASIBLE SET METHOD FOR AOA CONSTRAINT
If the inequality constraint is concave, a feasible half space can be obtained directly by the simple linearization method mentioned in the Ref. [23]. However, for nonconvexnonconcave constraints (See Figure 2), the intersection exists between the half space (the red region in Figure 2) and the infeasible region (the gray region in Figure 2) which will adversely affect the numerical solution. In this paper, we will  propose a novel second-order method to find the CFS (blue region in Figure 2). The The The gradient and Hessian do not exist at v = 0 or T = 0 (the feasible T can be zero after the relaxation in Sec. III B). To do this the function (8) is changed into the following form: The gradient off η (v, T ) is (12), as shown at the bottom of the next page. Define where v r and T r are the reference values. Now, we introduce a quadratic function where E = diag(e 1 , . . . , e n ), e i ≥ 0 for i = 1, . . . , n and n = 6. The Proposition 1 shows that there exists a diagonal matrix E such thatf η (v, T ) is a concave function.
Proposition 1: There exists such a diagonal matrix E such thatf η (v, T ) is a concave function. Proof: To make sure thatf η (v, T ) is concave, the eigenvalues ofH should be less than or equal to zero.
According to Gersgorin disc theorem [27], the eigenvalues are contained in the union of D i where z represents the eigenvalue ofH,H ij represents the elements inH, the subscript represents their position in the matrix, and i, j = 1, . . . , n. Eq Eq. (20) can make sure the eigenvalues ofH are less than or equal to zero.
Substituting Eq. (17) into Eq. (20), we get BecauseH ij is bounded, e i exists. This completes the proof. Define where the diagonal elements of E are set as Eq. (21). The set is equivalent to the AOA constraint. The proposition 2 shows that the set is a convex feasible subset of .
Proof: From the properties of concave function, we havē Substitute (16) into (24) to get This completes the proof.

2) ESTIMATION OF E
Due to the large range of v and T , the elements of E will be very large which may adversely affect the numerical solution. This paper presents an estimation method. Firstly, define the following new variables to shrink the range of the elements ofH. v Then, we can use the Gersgorin disc theorem (18)- (21) and The estimation (29) requires only a small amount of addition and absolute value operations and need no calculation of eigenvalues. The numerical evidence from our tests suggests the estimation (29) is effective, even though no rigorous proof is available yet.
The second order term in the set will shrink the AOA constraint, and increases the sensitivity to the initial guess. To reduce the sensitivity to initial guess, we introduce the following coefficient κ in the second term: where the κ is a number increasing from 0 to 1 with iteration, and

B. CONVEXFICATION OF THRUST CONSTRAINT AND DYNAMICS 1) THRUST CONSTRAINT
Lossless convexification and successive linearization methods are utilized to handle the nonconvex thrust constraint and nonlinear dynamics. Both methods are very classic, so only a brief process is given below.
According to the lossless convexification [9], a slack variety is introduced to augment the original space, the thrust constraint (6) can be relaxed to and replace T in Eq.(5) with . Proposition 3 will show that the relaxation (32)-(33) is lossless under AOA constraint.
Proposition 3: For problem P1, the relaxation (32)-(33) is lossless, i.e., where superscript * represents the optimal value. Proof: See the Appendix. where The implicit free terminal time t f in Eq. (16) is also a nonlinear term. Ref. [28] gives a way to express it explicitly. Define τ t/t f , τ ∈ [0, 1]. Apply the chain rule to the left side of (35): Eq. (35) can be rewritten in term of τ as follows: Eq. (18) is approximated by a first order Taylor expansion at a given reference point (x r , u r , t f ,r ) as follows where ∇ x f and ∇ u f are the gradients of f with respect to x and u respectively.
To ensure the validity of the linearization, trust region constraints are introduced as follows where p could be r, v, m, T , and t f .

C. DISCRETIZATION
The above convex optimization problem is still an infinite dimensional optimal control problem. To apply numerical methods to solve it, the continuous dynamics Eq. (19) must be discretized. In this paper, the trapezoidal discretization method is used. Eq. (19) can be rewritten as where i = 1, . . . , N − 1 and N is the number of discretized points.

IV. A SUCCESSIVE SOLUTION PROCEDURE
After convexification and discretization, the problem P1 can be transformed into a convex problem, donated as problem P2, which is expressed as follows P2: subject to (41) and where i = 1, . . . , N − 1.
The objective function, dynamics equation constraint, and state and control constraints are all convex in problem P2, which can be readily solved by the convex optimization algorithm with high efficiency.
The optimal solution of the original problem P1 can be obtained by iteratively solving convex problem P2. The flow diagram of the successive solution method is shown in Figure 3. Compared with the successive convex method in Ref. [2], the main difference between the two methods is that a step is added to calculate the matrix E before solving the problem P2 every time in the proposed method.
The stop criterion shown in Figure 3 is as follows: where ζ could be r, v, m , and ζ * represents the optimal solution obtained in the current iteration.   Proposition 4 explains the relationship between the solutions to problem P1 and P2.
Proposition 4: Suppose that the problem P2 is strictly feasible at the converged solution. If the optimal solution y * = [x * T u * T ] T of the problem P2 is y r itself, then y * = y r is a KKT solution to the problem P3. The problem P3 is the discrete form of problem P1.

V. NUMERICAL SIMULATIONS
The numerical simulation is divided into two parts. The first is to verify the effectiveness of the CFS method. The second is to analyze the convergence performance. The dynamics parameters are listed in Table 1 and the constraint parameters are listed in Table 2.
The number of discretized points N is 40. The trust region constraint boundary is The convergence criterion is where s represents the number of iteration and s min = 3 and s max = 6. The initial r r , v r , m r , r is chosen as a linear function of time from its initial value to the final value, and the thrust direction is simply set as the opposite direction of v r , i.e., To enhance the numerical stability, the parameters in problem P2 and simulation parameters are all normalized by the Earth's average radius R 0 = 6371km, g 0 , and initial mass m 0 of the rocket.
At each iteration, the MOSEK solver based on MATLAB is used to solve the problem P2 on a laptop computer with Intel Core i5-3320 CPU (2.60GHz) [29].

A. THE EFFECTIVENESS
To verify the effectiveness of the proposed method, the simulation results are compared with the solutions of problem P1 obtained via GPOPS that is a MATLAB software for solving optimal control problems [30]. GPOPS employs the Legendre-Gauss orthogonal collocation method to transform the continuous optimal control problem into a nonlinear programming (NLP) problem. The NLP solver used by GPOPS is set IPOPT [31]. The initial sates are listed in Table 3. Figures 4-8 and Table 4, from which it is found that the two methods produce very similar results. The slight difference may be caused by the differences between the discrete strategy of the two methods. Figure 7 verifies the relaxation on the thrust constraint is lossless.
As can be seen from Table 4, the difference of final time, fuel consumption obtained by the two methods is very small, which verifies the accuracy and optimality of the proposed method. The CPU time taken by the CFS method is 1.18s, which is only 5.3% of that by the GPOPS method, indicating the high efficiency of the CFS method.

B. THE CONVERGENCE 1) NOMINAL CASE
The convergence performance of the CFS method is analyzed under the nominal initial states listed in Table 4. The converged solution is obtained at the 7 th iteration. Figure 9 shows the trajectories in all iterative steps. The green trajectories are obtained in intermediate iterations. Figure 9 indicates that the trajectory can quickly get convergent. The convergence   process of the AOA profile is depicted in Figure 10. At the 4 th iteration, the AOA profile almost converges and completely meets the constraint requirement. Figure 11 shows that the fuel consumption gets convergent at the 4 th iteration. In general, Figures 9-11 verify the rapid convergence of the proposed method.

2) MONTE CARLO CAMPAIGN
In this subsection, we show the results associated with a Monte Carlo campaign (50 cases) to test the sensitivity of the VOLUME 8, 2020   method to the initial guess and the convergence performance. Perturbations in terms of the initial position, velocity, and mass are included. The ranges and the type of uncertainty are listed in the Table 5. States and controls are depicted in Figure 13. The simulation results in Figures 14-15 are compared with the solutions obtained via the successive   linearization (SL) method [18], i.e., All the cases in the CFS method use the initial guess (58) and can get convergent well. This shows that the proposed  method is insensitive to the initial guess. Figure 14 shows the comparison of the difference of the states and control between consecutive iterations between the CFS method and SL method. In about the 6 th -9 th iteration, all cases in the CFS method can meet the stop criterion while most cases in the SL method cannot. The oscillations or jumps occur after iteration 10 are quite small (about 10 −5 with dimension or 10 −12 without dimension). Those may be caused by the error of optimization solver. Figure 15 shows the comparison of AOA between CFS and SL. It can be seen that the AOA profiles obtained by CFS completely meet the constraint requirement while most of those obtained by SL cannot. This verifies that the proposed CFS method has a better convergence performance.

VI. CONCLUSION
In this paper, we propose a novel convex programming method to solve the rocket powered landing problem with the AOA constraint. The main contribution of the study is to propose a convex feasible method for the AOA constraint. We construct a quadratic concave function and use Gersgorin disc theorem to find a convex feasible set of the AOA constraint. Compared with the traditional linearization method for the AOA constraint, the main advantage is that the proposed method has a better convergence performance. The numerical simulations are performed to verify the effectiveness and convergence of the proposed method. Consequently, the proposed method has a potential for engineering applications.

A. PROOF OF PROPOSITION 3
Define Hamiltonian and Lagrangian as follows: Obtaining by the maximum principle in optimal control theory [32]: (1) The nontriviality condition: (2) The pointwise maximum condition: i.e., T * = * p 1 / p 1 , p 1 = 0 Singular arc, where R * v represents the partial derivative of aerodynamic drag R with respect to v.
Similarly, P3 can be converted into the following form. P5: Thus Y * satisfy all constraints (Eqs. (B10)-(B13)). Define the Lagrangian of P4 as [4]: The dual problem to P4 is given by [8] max ι≥0,υ≥0, The problem (B16) can be rewritten in the following form. D4: Because the problem P4 is a strictly feasible convex problem, the strong duality holds which means that the optimal duality gap is zero. Denote the optimal values of the objective functions of the primal and dual problems by p * and d * , respectively. For Y * = Y r , the zero duality gap implies: Because Eq.(B19) and P4 is strictly feasible, we have where i = 1, . . . , N − 1. Because all terms in Eq. (B21) are nonnegative, we have