A Reentry Trajectory Planning Algorithm via Pseudo-Spectral Convexification and Method of Multipliers

: The reentry trajectory planning problem of hypersonic vehicles is generally a continuous and nonconvex optimization problem, and it constitutes a critical challenge within the field of aerospace engineering. In this paper, an improved sequential convexification algorithm is proposed to solve it and achieve online trajectory planning. In the proposed algorithm, the Chebyshev pseudo-spectral method with high-accuracy approximation performance is first employed to discretize the continuous dynamic equations. Subsequently, based on the multipliers and linearization methods, the original nonconvex trajectory planning problem is transformed into a series of relaxed convex subproblems in the form of an augmented Lagrange function. Then, the interior point method is utilized to iteratively solve the relaxed convex subproblem until the expected convergence precision is achieved. The convex-optimization-based and multipliers methods guarantee the promotion of fast convergence precision, making it suitable for online trajectory planning applications. Finally, numerical simulations are conducted to verify the performance of the proposed algorithm. The simulation results show that the algorithm possesses better convergence performance, and the solution time can reach the level of seconds, which is more than 97% less than nonlinear programming algorithms, such as the sequential quadratic programming algorithm.


Introduction
Hypersonic vehicles generally refer to near-space vehicles with flight speeds greater than Mach 5.They possess the advantages of strong maneuverability, flexible trajectory, they are difficult to intercept, and so on, and they have been increasingly valued by the major space powers due to their high flight speeds and vast airspace coverage.Among the related technologies, trajectory planning technology can provide important support for performance analysis regarding flight range, maneuverability, and ballistic characteristics.The hypersonic vehicle trajectory planning problem typically involves solving a nonlinear optimal control problem with various state and control constraints, including boundary conditions, no-fly-zone constraints, path constraints, etc. [1][2][3].
In this regard, many scholars and engineers have carried out a series of in-depth studies, and the proposed trajectory planning algorithms generally include indirect and direct methods.The indirect methods transform the trajectory planning problem into a Hamiltonian boundary value problem based on the Pontryagin maximum principle and solves it by employing the gradient method and other algorithms.On the other hand, many scholars have proposed and developed the collocation method, pseudo-spectral method, and other methods [4][5][6][7], demonstrating the advantages of direct methods in solving trajectory optimization problems.However, the direct methods still have some shortcomings under the requirement of rapid trajectory planning, such as high sensitivity to initial guess value and uncertain solving time and convergence, which restrict their efficient solving ability [8].
In recent years, convex-optimization-based methods have effectively met the demand for efficient solutions and attracted more and more attention in terms of spacecraft trajectory optimization.When a problem can be formulated within a convex optimization framework, its complexity is low and can be reliably solved to global optimality in the polynomial time by the primal-dual interior point method.The upper bound of the number of iterations required for convergence is also determined.Moreover, the primal-dual interior point method can be adopted to solve the convex problem without the initial guess value.Motivated by these preponderances, the convex-optimization-based methods have been applied to different aerospace problems, such as planetary reentry trajectory optimization [9], ascent trajectory optimization [10,11], Mars-landing trajectory planning [12,13], low-thrust orbit transfer [14], spacecraft rendezvous and proximity operations [15,16], and trajectory planning for satellite cluster reconfigurations [17].
Most aerospace problems are limited by nonlinear, nonconvex dynamics and path constraints, so they cannot be solved directly under the convex optimization framework.Therefore, convexification technologies that make the approximate error as small as possible are a significant research direction [18].Among them, the two mainstream convexification methods include lossless convexification and sequential convexification methods.In Refs.[12,13], Ackimese et al. employed the lossless convexification method to solve the Mars-landing trajectory planning problem by replacing nonconvex constraints with relaxed convex constraints without a loss of accuracy.But the lossless convexification method is only suitable for a few constraints with particular forms, which limits its wide application.On the other hand, the sequential convexification method can conduct highly nonlinear complex problems.The basic idea is to obtain a series of convex subproblems by approximating the nonlinear terms and then solve them iteratively until it converges to the expected precision [19,20].
However, the sequential convexification method can be further improved in terms of discretization and accelerating convergence.Firstly, the traditional trapezoidal discrete method [3,19,20] is often chosen to discretize the continuous optimization problem, which leads to a large deviation between the approximate model used in the solution procedure and the actual one.To obtain a precise-enough solution, the equidistant discrete nodes should be sufficiently numerous.Nevertheless, it also results in a dramatic increase in the number of optimization variables and takes lots of time to solve.In contrast, the pseudo-spectral discretization methods offer higher accuracy under the same number of discrete nodes and have been widely employed for solving the optimal control problem [5,8].Among the pseudo-spectral methods, the Chebyshev pseudo-spectral discretization method [21,22] is a special category, which can minimize the Runge phenomenon and supply the best polynomial approximation under the minimax norm.On the other hand, "artificial infeasibility" [23,24] caused by the convexification errors occurs when the original problem is feasible but the convex subproblem is not.Researchers address this issue by introducing slack variables to relax the feasible domain limited by various constraints, with large constant penalty parameters added to penalize these slack variables.As a result, the "artificial infeasible" gradually disappears and the slack variables tend towards zero as the iterative solution converges.Unfortunately, the fixed penalty parameters could cause the solution to converge to a stagnation point of the penalty problem rather than the original problem, according to Ref. [25].
To address the aforementioned issues, an improved sequential convexification algorithm is proposed to improve the performance of solving the trajectory planning problem of hypersonic vehicles in this paper.Firstly, the Chebyshev pseudo-spectral method with higher approximation precision is employed to discretize the continuous optimal control problem.And the flight terminal time is designed as an optimization variable, so that it can be applied to the optimal control problem with a fixed initial state and free terminal state.Then, by introducing automatically updated penalty parameters and Lagrange multipliers, the relaxed convex subproblem in the form of the augmented Lagrange function is constructed to improve the convergence and computational properties based on the method of multipliers (or the augmented Lagrange method).And it is solved iteratively by the interior point method in the framework of the improved sequential convexification algorithm.
The rest of this paper is organized as follows.In Section 2, the reentry trajectory optimization problem is formulated.In Section 3, the improved sequential convexification algorithm is detailed.In Section 4, numerical simulations are presented to verify the performance of the algorithm.Finally, the conclusion and discussion are provided in Section 5.

Problem Formulation 2.1. Reentry Dynamics
In view of the dynamics characteristics of the reentry flight of hypersonic vehicles, the dimensionless three-degrees-of-freedom augmented dynamics model is established in the half-velocity coordinate system without considering Earth's flatness and rotation: where r is the radial distance from the Earth center to the vehicle, θ and ϕ are the longitude and latitude, respectively, v is the velocity, γ is the flight path angle, ψ is the heading angle, and σ is the bank angle.The above variables are defined as system state variables, i.e., x = [r, v, θ, ϕ, γ, ψ, σ] T , and the control variable is the bank angle rate u ∈ R. The variables r, v and time are scaled by R 0 , R 0 g 0 and R 0 /g 0 , respectively, where R 0 = 6371 km is the Earth radius and g 0 = 9.8 m/s 2 is the gravitational acceleration at sea level.Dimensionless lift and drag accelerations L and D are scaled by g 0 and calculated as: where m is the vehicle mass, S is the reference area of the vehicle, and ρ is the atmospheric density, ρ = ρ 0 e −H/H s , where H = rR 0 − R 0 is the height, H s is the atmospheric density scale height, and ρ 0 is the atmospheric density at sea level.Moreover, C L and C D are lift and drag coefficients, respectively.For reference to the aerodynamic parameters of the vehicle in [9], the aerodynamic coefficient can be expressed as: in which the angle of attack α is in degree and expressed as a function of velocity.The angle-of-attack velocity profile is preset as follows: In general, the design of the bank angle is the main means to change the trajectory of the vehicle when the angle of attack is preset.However, the bank angle rate is adopted here as the new control variable for the following reasons.Firstly, boundary constraints are applied to make the bank angle profile smoother.Secondly, the current system dynamics are transformed into affine control form, which expedite subsequent convexification operations and alleviate the high-frequency fluctuations in the bank angle [9].

Constraint Conditions
The trajectory planning problem of hypersonic vehicles belongs to a class of optimal control problems, and its purpose is to determine the optimal control variables and optimize the performance index function under various constraints.
Generally, the vehicle must meet different constraints, such as heat rate, dynamic pressure, load, and control margin, during the ultra-high-speed flight to ensure the safety of the vehicle structure, thermal protection, guidance, and control systems.To meet the mission requirements, the initial conditions and terminal states of the vehicle also need to be limited.In addition, the vehicle should avoid certain no-fly zones due to radar detection, geopolitics, or other considerations.In summary, all of these constraint conditions can be divided into two categories: equality constraints and inequality constraints.
First of all, the equality constraints mainly include dynamic equations and initial and terminal constraints.Here, the system dynamics Equation ( 1) is abstractly expressed as where t is the system time variable, and f (•, •, •) is the right function of the dynamic equations.And The initial and terminal constraints are determined by the flight mission, which contains the requirements of the reentry start point and the target point, and they can be represented as follows: where t 0 and t f are the initial and terminal time, respectively.Φ(•, •) and Ψ(•, •) are the initial and terminal state constraints, respectively.Then, inequality constraints can be divided into the following three categories.First, the path constraints, including the maximal heat rate, dynamic pressure, and load, are expressed as: Q max , q max and n max are the corresponding maximum values, respectively.Second, the bounded constraints about the state and control variables are given by: where • min and • max are the lower and upper bound values, respectively.
The last class of inequality constraints is the no-fly-zone (NFZ) constraint.In general, an NFZ is modeled as a circular exclusion zone in the horizontal place with infinite height and limits the longitude and latitude of the flight trajectory by where θ c and ϕ c are the longitude and latitude at the center of the NFZ, and d is the radius.

Optimal Control Problem
In the trajectory optimization problem, the performance indexes can be reasonably selected according to the different flight tasks and design requirements.Common performance indexes include minimum flight time, maximum range, minimum heat load, etc.In this paper, the minimum flight time is chosen as the performance index, considering that the vehicle needs to reach the anticipative target point quickly.Henceforth, the original reentry trajectory optimization problem can be formulated as a highly constrained optimal control problem: 5), ( 6), ( 7), ( 8), (9) .

Improved Sequential Convexification Algorithm
In this section, an augmented Lagrange-based Chebyshev pseudo-spectral form improved sequential convexification (AL-CP-ISC) algorithm is proposed to solve the above continuous and nonlinear reentry trajectory optimization problem.The advantages of the pseudo-spectral method with high discrete accuracy and the method of multipliers with good convergence performance and stable numerical computation are synthesized to empower the AL-CP-ISC algorithm.Firstly, the discretization and convexification process of the vehicle's nonlinear system are given on the strength of the Chebyshev pseudospectral method and first-order Taylor expansion.Subsequently, the original problem is transformed into a series of relaxed convex subproblems by introducing the slack variables, penalty parameters, and Lagrange multipliers.Finally, an algorithm solution procedure is presented.

Discretization and Convexification
The Chebyshev pseudo-spectral discretization method with unique time-domain mapping is adopted to discretize the continuous reentry trajectory optimization problem P 0 with the free terminal time.In the Chebyshev pseudo-spectral discretization method, the domain of the Chebyshev-Gauss-Lobatto (CGL) points is τ ∈ [−1, 1], but the flight time interval is t ∈ [t 0 , t f ] in the practical problem.Hence, the time variable is transformed into The CGL points are unevenly distributed on the interval [−1, 1]: Taking the real state and control variables at N + 1 nodes above, the Lagrange interpolation polynomials are constructed, respectively, as approximations of continuous state and control variables.The approximate expressions of the real state variable x and the control variable u are where ϕ j (τ) is the N th Lagrange interpolation basis function where T N (τ) is the Nth − order Chebyshev polynomial and Deriving the approximate expression for the state variables x(τ) yields: .
where D kj is the k th row and j th column element of the differential matrix D (N+1)×(N+1) of the Chebyshev pseudo-spectral discretization method, and the calculation of each element's value of the matrix D is shown in Ref. [21].
The derivatives of the state variables in the dynamic equations are replaced using the right-hand term in Equation ( 15) and discretized at the nodes, so that the original differential dynamic equation constraints are transformed into discrete algebraic constraints: Next, the above nonlinear algebraic constraints are linearized based on first-order Taylor expansion as follows: where x * , u * and t * f are the reference values of the optimization variables [x, u, t f ], respectively, and One can see Appendix A for more details on the matrix A. Similarly, the nonlinear path constraints are given by: where The NFZ constraint is a nonconvex function and also needs to be convexified by the first-order Taylor expansion, as shown below: where Finally, to place limits on the deviation of the state variables between the linearized and original system, a trust region constraint is introduced, so as to reduce the linearization error and improve the convergence of the sequential linear approximation, denoted as where δ ∈ R 7 is a constant vector, and the inequality is expressed in components.This is a second-order conic constraint that is, itself, convex.Adding the trust region constraint is necessary to guarantee the linearized constraints to legitimately approximate the original constraints.Meanwhile, proposition 2 in Ref. [9] theoretically explains that a feasible solution of the linearized problem satisfies the linearized path constraints, and it also satisfies the original path constraints.

Problem Transformation
Combining the discretization and convexification processes, the original problem is transformed into an augmented Lagrange formal convex problem.The AL-CP-ISC algorithm can profit from this conversion, resulting in convergence-rate promotion and numerical difficulty avoidance according to the following analysis [26][27][28].
The linearized error causes the feasible region of the original problem to shrink to a certain extent in the process of convexification, resulting in the "artificial infeasible" situation in which the original problem is feasible but the linearized problem is not feasible.Hereon, slack variables are introduced to relax the "hard constraints", such as dynamic equations, path constraints, and terminal conditions, to compensate for the linearization errors [8,18].In the meantime, penalty parameters are introduced to punish the slack variables.When the iterative solution converges, the slack variables also gradually approach zero due to the penalty imposed.However, the penalty parameters are usually selected as a larger constant value in many studies of spacecraft trajectory optimization.Unfortunately, the constant penalty parameters can cause the iterative solution to converge to the stationary point of the penalty problem instead of the original problem [25].
Therefore, to resolve this matter, the slack variables are penalized by penalty parameters, which are automatically updated incrementally.Meanwhile, to prevent the condition number of the Hessian matrix from getting worse and worse when the penalty parameters are updated to infinity, which leads to numerical difficulties in the algorithm, Lagrange multipliers are introduced to transform the penalty function problem into the augmented Lagrange function one.According to the augmented Lagrange function method, by embedding the multiplier update mechanism, the satisfaction degree of constraints and the optimization of the objective function can be considered in each iteration, thus reducing the number of iterations required to achieve the same precision.Additionally, the augmented Lagrange function method can also effectively control the condition number of the Hessian matrix through its unique construction, which avoids the instability of numerical computation.In short, the iterative solution of the algorithm can converge to the solution of the original problem, and the penalty parameters do not need to go to infinity.
Here, the original problem is abstracted into a general nonlinear programming problem form to illustrate the relaxed convex problem of an augmented Lagrange form, as shown below: , where x = [x, u, t f ] T are augmented optimization variables, all equality constraints are denoted as h( x), all inequality constraints are expressed as g( x), and s( x) represents the second-order cone constraints.And E, I, and C are the corresponding feasible domain sets, respectively.Then, the slack variables ζ h , ζ g , penalty parameters p, and Lagrange multipliers µ, λ are introduced to construct the relaxed problem of augmented Lagrange function form as follows: , where, for further transformation, the penalty term of the inequality function is not treated as a standard cutoff function because the slack variables ζ g are nonnegative but as an absolute value term.The similar treatment method can also be found in Ref. [25].
Finally, the problem P 2 is convexified by the first-order Taylor expansion: The objective function of the problem P 3 contains an absolute value term, so it is not a convex second-order cone problem.For this purpose, auxiliary control variables η can be introduced to convert the objective function with absolute values into a combination of the linear objective function and the second-order cone constraints [29].Specifically, an equivalent unconstrained optimization problem is given below min ∑ p|x|, p > 0.
The minimization problem is equivalent to min ∑ pη s.t.|x| ≤ η , in which the constraint condition is a second-order cone function, and it is naturally convex.
At this point, a standard convex second-order cone problem P SOCP is represented as: in which it contains linear objective function, affine equality constraints, and second-order cone constraints.And T is the set of the augmented optimization variables, slack variables, and auxiliary control variables.Suppose N x is the dimension of x after the discretization.The variables N C and N I are the number of the second-order cone and inequality constraints, respectively, and let N E be the number of the equality constraints.
The parameters , and d i ∈ R are calculated based on the previous iteration solution T .It should be noted that the linear inequality constraints can be included in equation (*) because it is a particular case of the second-order cone constraints.And the final convex subproblem P SOCP can be solved iteratively by the advanced interior point method.
For the penalty parameters and Lagrange multipliers in a discrete convex problem P 3 , the specific updating methods are as follows: where ρ p > 1 is the update multiple of the penalty parameters.Additionally, to accelerate the convergence of the algorithm, an update approach is designed for the trust region radius in the second-order cone constraints.Let k be the number of iterations and k sc > 1 be the number of iterations in which the trust region starts to update.When where 0 < b 0 < 1 < b 1 are the trust region contraction factors.Hereon, the trust region is updated in the middle step of the iterations, rather than in the first step.It is to ensure that there is a large manually set scope of trust regions in the early iterations when the system has a certain approximation error.In this way, the algorithm can easily find a feasible optimization direction and perform continuous iterations.

Solution Procedure
The solution procedure of the AL-CP-ISC to find the solution to the original problem is given.As shown in Figure 1, the convex subproblem can be solved iteratively until predetermined convergence precision is reached.x by propagating the dynamical Equation ( 1) with the fixed control variables.

Set initial value
4. Define the value of constraint violation v : ( ) When k v κ ≤ , go to 6, otherwise, go to 5, where κ is a sufficiently small posi- tive number.
5. Update penalty parameters k p and Lagrange multipliers , Then, 1 k k = + , and go to 3.
6. Obtain the optimal solution of the original problem: A more detailed solution procedure of the AL-CP-ISC algorithm is shown as follows (Algorithm 1): Let k = 0, set the initial reference trajectory x0 by propagating the dynamical Equation ( 1) with the fixed control variables.2. Assign initial values to the following parameters: penalty parameters p 0 , penalty parameter update multiple ρ p > 1, initial Lagrange multipliers λ 0 , µ 0 , initial trust region radius δ 0 , the iteration number of the trust region starts to update k sc , the trust region contraction factor 0 < b 0 < 1 < b 1 , and the number of discrete points N. 3. k ≥ 1, solve the convex subproblem P SOCP by the interior point method, and find solution pairs 4. Define the value of constraint violation v: When v k ≤ κ, go to 6, otherwise, go to 5, where κ is a sufficiently small positive number.5. Update penalty parameters p k and Lagrange multipliers λ k , µ k .Then, k = k + 1, and go to 3. 6. Obtain the optimal solution of the original problem: Remark 1.In Step 1, the initial reference trajectory generated by the numerical integration method meets the dynamic constraints, and the change in the bank angle and its rate are smooth.It is beneficial to enable linearization and facilitate the iteration.Nevertheless, notice that the quality of this initial trajectory may affect the convergence effect, such as the number of iterations and convergence accuracy.Therefore, taking into account factors such as the flight time, mission characters, and other constraint conditions, it is necessary to judiciously select the fixed control variable.

Numerical Verification
In this section, numerical simulations are employed to testify the effectiveness and convergence of the AL-CP-ISC algorithm.To compare with the simulation results of the AL-CP-ISC algorithm, CPM and P-CP-ISC algorithms are used to solve the same reentry trajectory optimization problem of the hypersonic vehicle.The CPM is the Chebyshev pseudo-spectral method that transforms the original problem into a general nonlinear programming problem and adopts the sequential quadratic programming algorithm to solve the problem.The difference between P-CP-ISC and AL-CP-ISC lies in different transformation approaches, and the objective function of the P-CP-ISC algorithm adopts the form of a penalty function without Lagrange multipliers.
According to Ref.
[30], it is considered that the mass of the reentry vehicle is m = 104, 035 kg, and the aerodynamic reference area is S = 391.22m 2 .In this paper, the fight mission of the vehicle is to plan an optimal trajectory with the minimum flight time under various constraints.Correspondingly, the initial and terminal constraints are: Q max = 1500 kW/m 2 , q max = 18, 000 N/m 2 , and n max = 2.5g 0 , respectively.In addition, to keep enough control margin, the magnitude of the bank angle and bank angle rate are limited by 80 • and 10 • /s.The longitude and latitude of the center of the NFZ constraint are θ c = 2 • and ϕ c = 50 • , respectively, and the radius is about 222 km.The initial trust region size in Equation ( 26) is given as: T .
In the solution procedure of the AL-CP-ISC algorithm, the termination condition is set to be v k ≤ 10 −9 , and the number of CGL points is 80.In Step 1 of the algorithm procedure, the initial reference trajectory is generated by the numerical integration method, in which the fixed control variable is set to 0.015 • /s, and the estimated terminal time is 1610 s.All optimization algorithms in this section are implemented by MATLAB 2020b on a desktop computer equipped with Intel Core i7-10700K/3.80GHz CPU.For the convex subproblems P SOCP , the professional software MOSEK (Version 9. 3. 7.) [31] is used to solve them.
First, Figures 2-7 present the trajectories of the vehicle obtained by the three algorithms, and Table 1 shows all the optimal solutions and solve time of each algorithm.The solid red curves are the solution of the AL-CP-ISC, the blue double lines represent the results of the CPM, and the other black dotted lines are the ones of the P-CP-ISC.In Figure 2, the profiles of each iteration of the AL-CP-ISC algorithm are depicted in different colors, from dark blue to warm red, to make the progression of the convergence clearer.The sequence solution converges in the 25th iteration, and the corresponding minimum flight time is 1636.68 s, which is quite close to the optimal solution from the other two algorithms.The trajectory curves almost overlap and are hard to distinguish in late iterations according to the zoom-in view.Moreover, Figures 3-6 show a comparison of results between state variables and control variables, including altitude, velocity, longitude, latitude, flight-path angle, heading angle, bank angle, and its rate.It is obvious that all state and control profiles tend to be relatively consistent, and the initial, terminal, and bounded conditions required by the original problem are satisfied.In particular, it can be seen from Figure 4 that the trajectory of the vehicle successfully evades the no-fly zone.Meantime, the vehicle trajectories obtained by the three algorithms all meet the conditions of path constraints in Figure 7.In Table 1, in addition to the similar optimal solutions, it is worth noting that the solve time of AL-CP-ISC and P-CP-ISC is only 4.20 s and 5.30 s, far less than the solve time required by CPM.Given these analyses, it can be fully explained that the improvement measures of the sequential convexification algorithm do not affect the optimal solution, and the AL-CP-ISC algorithm can vastly reduce the computational cost; that is, the effectiveness and efficiency of the proposed algorithm are verified.
ures are the value of constraint violations and the terminal flight time in each it respectively.The only difference between the two algorithms is the form of the o function.In a numerical sense, it indicates that the AL-CP-ISC algorithm with t mented Lagrange function formal cost function has better convergence properties paring the trajectory results of these two algorithms.Then, the state trajectories of the AL-CP-ISC and P-CP-ISC algorithms are almost the same, and some slight differences are seen in the optimal flight time and the solve time.This shows that both algorithms can solve the approximate convex problem well.However, the number of iterations required by both algorithms is quite different according to Figures 8  and 9, and they converge at the 25th and 32th iterations, respectively.Both figures are the value of constraint violations and the terminal flight time in each iteration, respectively.The only difference between the two algorithms is the form of the objective function.In a numerical sense, it indicates that the AL-CP-ISC algorithm with the augmented Lagrange function formal cost function has better convergence properties by comparing the trajectory results of these two algorithms.

Conclusions
In this paper, an improved sequential convexification algorithm is proposed to efficiently solve the trajectory planning problem of hypersonic vehicles based on the Chebyshev pseudo-spectral method and the method of multipliers.By employing the Chebyshev pseudo-spectral discretization method, the differential dynamical equations are transformed into the algebraic equation constraint.The original problem is then discretized into a finite-dimensional nonlinear programming problem.Next, through linearization and relaxation techniques, the slack variables, penalty parameters, and multipliers are introduced to transform the discrete problem into a series of relaxed convex subproblems in the form of the augmented Lagrange function, which are iteratively and efficiently solved by the MOSEK solver until the expected convergence precision is satisfied.The

Conclusions
In this paper, an improved sequential convexification algorithm is proposed to efficiently solve the trajectory planning problem of hypersonic vehicles based on the Chebyshev pseudo-spectral method and the method of multipliers.By employing the Chebyshev pseudo-spectral discretization method, the differential dynamical equations are transformed into the algebraic equation constraint.The original problem is then discretized into a finite-dimensional nonlinear programming problem.Next, through linearization and relaxation techniques, the slack variables, penalty parameters, and multipliers are introduced to transform the discrete problem into a series of relaxed convex subproblems in the form of the augmented Lagrange function, which are iteratively and efficiently solved by the MOSEK solver until the expected convergence precision is satisfied.The numerical simulation results verify the effectiveness, efficiency, and convergence performance of the proposed AL-CP-ISC algorithm.In the future, the proposed algorithm will be competitive for use in realizing onboard optimization due to its excellent optimization efficiency after code optimization.

Appendix A
All elements of the matrix A in Equation ( 17

Figure 1 .Algorithm 1 .
Figure 1.Algorithm solution flow.A more detailed solution procedure of the AL-CP-ISC algorithm is shown as follows (Algorithm 1):

2 .δ
Assign initial values to the following parameters: penalty parameters 0 , the iteration number of the trust region starts to update sc k , and the number of discrete points N .3. 1 k ≥ , solve the convex subproblem SOCP P by the interior point method, and find solution pairs: ˆ[ , , , , ]
The maximum values of the path constraints are .

Figure 2 .
Figure 2. Convergence of the trajectories using AL-CP-ISC.

Figure 5 .
Figure 5.Comparison of fight-path angle and heading angle profile.Figure 5. Comparison of fight-path angle and heading angle profile.

Figure 5 .
Figure 5.Comparison of fight-path angle and heading angle profile.Figure 5. Comparison of fight-path angle and heading angle profile.

Figure 5 .
Figure 5.Comparison of fight-path angle and heading angle profile.

Figure 6 .
Figure 6.Comparison of bank angle and bank angle rate profile.

Figure 7 .
Figure 7.Comparison of the path constraints.

Figure 6 .
Figure 6.Comparison of bank angle and bank angle rate profile.

Figure 5 .
Figure 5.Comparison of fight-path angle and heading angle profile.

Figure 6 .
Figure 6.Comparison of bank angle and bank angle rate profile.

Figure 7 .
Figure 7.Comparison of the path constraints.Figure 7. Comparison of the path constraints.

Figure 7 .
Figure 7.Comparison of the path constraints.Figure 7. Comparison of the path constraints.

Figure 8 .
Figure 8. Value of constraints violation for each iteration.Figure 8. Value of constraints violation for each iteration.

Figure 8 .
Figure 8. Value of constraints violation for each iteration.Figure 8. Value of constraints violation for each iteration.

Figure 8 .
Figure 8. Value of constraints violation for each iteration.

Figure 9 .
Figure 9. Flight terminal time for each iteration.

Figure 9 .
Figure 9. Flight terminal time for each iteration.

Table 1 .
Comparison of the simulation results.

Table 1 .
Comparison of the simulation results.