A Trajectory Generation Algorithm for a Re-Entry Gliding Vehicle Based on Convex Optimization in the Flight Range Domain and Distributed Grid Points Adjustment

Featured Application

The method proposed can be used for accurate trajectory generation of a re-entry glide vehicle in the flight range domain with distributed grid points. It also provides ideas for trajectory generation based on the guidance mechanism.

Abstract

Optimal trajectory generation for the guidance of re-entry glide vehicles is of great significance. To realize a faster generation speed and consistency with the guidance mechanism, an improved convex optimization trajectory generation algorithm based on the flight range domain for the re-entry glide vehicles is proposed in this paper. Firstly, according to the definition of the range-to-go, the projected range-to-go of the re-entry glide vehicle is presented when the dynamic model is converted to the flight range domain. Then, the attack angle and bank angle are expanded to the state variables and the change rate, which is designed as a new control variable. When the dynamic models and constraints are convexificated and discretized, the vehicle trajectory discrete convex model in the flight range domain is proposed. In order to further improve the generation speed and accuracy, an initial trajectory generation method that is close to the guidance requirements is proposed by the landing points of different control laws. In addition, by analyzing the nonlinear illegal degree of grid points, the probability density of grid points and the adjustment strategy of grid points are proposed. Finally, the ablation experiment shows that the initial trajectory generation and distributed grid points method works. With different target points, different no-fly zones, different initial states, and the Monte Carlo experiment, our method can effectively and robustly complete the generation. The lateral and longitudinal generation error is less than 1 km. And compared with the Gaussian pseudo-spectral method, our method obtained comparable accuracy and faster speed.

1. Introduction

The re-entry glide vehicle [1] with a lift body structure has extremely fast re-entry and glide speeds and can make a jump maneuver of tens of kilometers in the longitudinal direction with hundreds of kilometers of maneuvering range in the lateral direction, which is difficult to predict and intercept [2], and thus has high research value.

Research on the trajectory generation mechanism of re-entry glide vehicles is beneficial to generating more effective re-entry trajectories and improving guidance accuracy. In general, the trajectory generation methods of re-entry glide vehicles are mainly divided into predictor–corrector guidance and nominal trajectory guidance. At present, the numerical predictor–corrector method is adopted within the mainstream predictor–corrector guidance methods, and the process is mainly divided into longitudinal and lateral guidance. Longitudinal guidance is primarily realized by setting the attack angle profile and using the secant method [3], the landing points method [4], or the flight range prediction network [5], to calculate the value of the bank angle. In the lateral direction, in order to consider the guidance and the avoidance of the no-fly zones, the avoidance logic [6] is usually combined with the heading deviation angle corridor [7] and cross-range corridor [8] when the adaptive cross-range corridor [5] integrating the two functions is also used to realize guidance. In general, although the trajectory generation speed of the predictor-corrector algorithm is fast, its guidance accuracy can be further improved.

At the same time, the nominal trajectory generation of a re-entry glide vehicle is a complex nonlinear problem. The primary method is transforming the problem of re-entry trajectory generation into the optimization problem of discrete points through the sequential convex optimization method [9]. In recent years, several types of research have been performed. Liu [10] et al. and Wang [11] et al. solved the multiple constraint re-entry optimization problem based on the sequential quadratic programming (SQP) method through linearization and variable relaxation. Wang [12] et al. treated the online trajectory generation problem of re-entry glider aircraft as a second-order cone programming problem and designed the corresponding optimal feedback tracking guidance design for effective tracking. Hong [13] et al. realized continuous guidance of re-entry vehicles based on the convex optimization algorithm and continuous closed-loop control. Liu [14,15] et al. introduced a nonlinear relaxation method, which can reduce the relaxation degree, and the effectiveness of the method is verified. Wang [16] et al. used the convexity method to provide a high-quality initial trajectory guess for the pseudo-spectral method. All the above methods realized the convexification of the re-entry dynamic model and constraints of the vehicle in the time domain and obtained the generated re-entry trajectory. Sandberg [17] et al. compared the optimal trajectories using Pontryagin’s method and the slew trajectories using sinusoidal functions, and approximately 1.5% lower control effort was obtained. Then, based on Sandberg’s approach, Raigoza [18] established an autonomous trajectory maneuver to de-orbit spacecraft back to Earth with a distributed waypoint for autonomous collision avoidance. However, according to the guidance mechanism, the planned spatial state of grid points is not directly reflected in the time domain, more states can be considered to describe the dynamic model.

In order to improve the convergence and accuracy of trajectory generation based on convex programming, many improved methods have been proposed. Wang [19] and Zhou [20,21] introduced new variables to adjust the trust region in the aircraft planning process adaptively, and this can realize faster convergence. Saglino [22] et al. converted the guidance problem of the re-entry glide vehicle to the energy domain, reduced the non-convexity of the equation, and solved it with the pseudo-spectral method. Liu [14] et al. realized the conversion of the guidance problem of the vehicle to the altitude domain and also realized the practical re-entry guidance. From the above methods, it can be seen that changing the described domain of dynamic equations may realize effective convexity and convergence. At the same time, the operation of the variable trust region is also worded to increase the convergence accuracy of the generated trajectory.

For the grid point adjustment method of re-entry trajectory generation, the Gaussian and Radau pseudo-spectral method are proven to be effective [23]. The hp pseudo-spectral method is used to select the root of the orthogonal polynomial as the grid points [24], which can be dynamically adjusted with iteration. Li [25] et al. proposed a multi-segment grid point modified Radau pseudo-spectral method to realize the generation and fast convergence of convex optimized trajectory. However, the pseudo-spectral method is better for linear problems and convex constraints and is not fully applicable to sequential planning problems. In addition to the pseudo-spectral method, some adaptive grid points adjustment methods are also proposed. Zhao [26] et al. proposed an adaptive grid point adjustment method, which solved the problem of Mars re-entry trajectory optimization by analyzing iteration errors to realize insertion and deletion points. Zhou [27] put forward the concept of nonlinear illegal degree and used the error threshold to be the insertion and deletion rules of grid points, thus realizing the dynamic adjustment of grid points. The concept of nonlinear illegal degree proposed by the above methods has an excellent guiding significance for grid point adjustment, but the grid point adjustment method is fixed, and the increase or decrease may be significantly affected by the convergence process.

To solve the problem that the description of re-entry trajectory generation in the time domain is not direct enough, the initial trajectory of iteration affects the convergence speed, and the grid point adjustment mode is fixed, a new convex optimized trajectory generation with the flight range domain and distributed grid points adjustment are proposed in this paper. Firstly, based on the definition of range-to-go, the concept of projected range-to-go is proposed, by which the dynamic model of the vehicle is transformed to the flight range domain for representation. At the same time, the states of the vehicle are expanded, and the dynamic model and constraints of the glider are convex and discrete. Then, the re-entry trajectory generation problem is transformed into a sequential programming problem. In order to improve the solving speed, the initial trajectory that basically meets the guidance requirements is obtained through transformation and interpolation according to the landing points of the vehicle under different constant control laws. In order to further improve the trajectory generation accuracy and speed up the trajectory convergence, a grid points probability density function is proposed according to the nonlinear illegal degree of the iterative trajectory, which is used to change the rules of generating new grid points and adjust the location of grid points. Through ablation experiments, the effectiveness of the improved part is proved. Additionally, by changing the initial and terminal positions and conducting Monte Carlo experiments, the robustness of our method is verified. The generation errors of our method are less than 1 km. Finally, through comparison with other mainstream methods, we prove the superiority of our method in search speed and accuracy.

The innovation of this paper mainly includes the following three parts:

• According to the concept of range-to-go in the guidance of re-entry glide vehicles, the projected range-to-go is proposed. By the definition of projected range-to-go, the dynamic model of the vehicle is transformed from the time domain to the flight range domain. Then the dynamic model and constraints are convexification and discretization, and the final sequential convex optimization expression is proposed;
• According to the landing points under the different constant control laws, the initial trajectory generation problem of the vehicle in any spatial state is transformed into a similar initial state trajectory generation problem by rotation transformation. And the initial trajectory that can basically realize practical guidance and meet the dynamic model is obtained by interpolation;
• According to the concept of the nonlinear illegal degree of iteration trajectory and the distribution of the nonlinear illegal degree, the grid points probability density function is proposed, and the grid points adjustment law is proposed by the probability density function, which realizes efficient and fast grid points adjustment.

The section of this paper is arranged as follows: Section 1 describes the research status of convex trajectory generation of re-entry glide vehicles and the innovation of this paper; in the Section 2, the dynamic model and constraints of the trajectory generation in the flight range domain are described; in the Section 3, the convexity and discretization process of the trajectory programming in the flight range domain are described; in the Section 4, the initial trajectory generation method and grid point adjustment strategy are proposed; in the Section 5, operation of the simulation is used to verify the research content; the Section 6 summarizes the conclusion of this paper.

2. Dynamic Models and Constraints of Re-Entry Glide Vehicle

2.1. Dimensionless Dynamic Model in the Flight Range Domain

In order to facilitate the trajectory generation of the vehicle, the dynamic model of the vehicle is usually dimensionless in the time domain. The dimensionless dynamic equation is

$$\left\{\begin{array}{l}\dot{r}=vsin\theta \\ \dot{\varphi }=vcos\theta sin\psi /rcos\phi \\ \dot{\phi }=vcos\theta cos\psi /r\\ \dot{v}=-D-sin\theta /r^2\\ \dot{\theta }=Lcos\beta /v+vcos\theta /r+cos\theta /\left(v{r}^2\right)+C_{\theta }+{\tilde{C}}_{\theta }\\ \dot{\psi }=Lsin\beta /\left(vcos\theta \right)+vcos\theta sin\psi tan\phi /r+C_{\psi }+{\tilde{C}}_{\psi }\end{array}\right.$$

(1)

where, $r$ represents the dimensionless geocentric distance of the vehicle, $\varphi$, $\phi$ represent the longitude and latitude coordinates, $v$ represents the dimensionless speed, $\theta$ represents the flight path angle, $\psi$ represents the flight heading angle of the target, $L$ and $D$ represent the dimensionless lift and drag of the vehicle, $C_{\theta }$, ${\tilde{C}}_{\theta }$ and $C_{\psi }$, ${\tilde{C}}_{\psi }$ represent the coriolis acceleration item and the involved acceleration item caused by the earth’s rotation.

$$\left\{\begin{array}{l}L=\left(\rho {\left(v\cdot v_c\right)}^{2}S{C}_L\right)/2m{g}_0\\ D=\left(\rho {\left(v\cdot v_c\right)}^{2}S{C}_D\right)/2m{g}_0\end{array}\right.$$

(2)

where, $C_L$ represents the lift coefficient, $C_D$ represents the drag coefficient, $S$ represents the reference area of the vehicle, $m$ represents the mass of the vehicle, and $g_0$ represents the gravity acceleration at zero altitude, $v_c=\sqrt{g_0{R}_e}$ represents the dimensional velocity where $R_e$ is the radius of earth, $\rho$ is the atmosphere density. The dimensionless operation of other quantities can be referred to in the literature [28].

Define the initial range-to-go of the re-entry glide vehicle as $S_{f0}$.

$$S_{f0}=\arccos \left(cos{\phi }_{f}cos{\phi }_{0}cos({\varphi }_f-{\varphi }_0)+sin{\phi }_{f}sin{\phi }_0\right)$$

(3)

where, ${\varphi }_0,{\phi }_0$ represent the longitude and latitude coordinates of the initial point, and ${\varphi }_f,{\phi }_f$ represent the longitude and latitude coordinates of the target point.

Then the speed of the vehicle between the initial point and the target point is $v_s=vcos\theta cos(\psi -{\psi }_p)$, as shown in Figure 1. When the initial point and target point are given, the problem of effective trajectory generation for the vehicle can be regarded as a problem of decreasing the projected range-to-go of the vehicle between the initial point with the target point. If the projected range-to-go is $s$, then

$$ds/dt=-vcos\theta cos(\psi -{\psi }_p)/r$$

(4)

where ${\psi }_p$ represents the flight heading angle of the target point relative to the initial point, ${\psi }_p=\arcsin \left(sin({\varphi }_f-{\varphi }_0)cos{\phi }_f/sin{S}_f\right)$.

Operating Equation (1) with Equation (4), the dynamic model of the vehicle in the flight range domain can be expressed as

$$\left\{\begin{array}{l}dr/ds=-rtan\theta /cos(\psi -{\psi }_p)\\ d\varphi /ds=-sin\psi /\left(cos\phi cos(\psi -{\psi }_p)\right)\\ d\phi /ds=-cos\psi /cos(\psi -{\psi }_p)\\ dv/ds=Dr/vcos\theta cos(\psi -{\psi }_p)+tan\theta /vrcos(\psi -{\psi }_p)\\ d\theta /ds=-Lcos\beta r/\left(v^{2}cos\theta cos(\psi -{\psi }_p)\right)-1/cos(\psi -{\psi }_p)+1/\left(v^{2}rcos(\psi -{\psi }_p)\right)+\left(C_{\theta }+{\tilde{C}}_{\theta }\right)\\ d\psi /ds=-Lsin\beta r/\left(v^2{cos}^2\theta cos(\psi -{\psi }_p)\right)-sin\psi tan\phi /\left(cos(\psi -{\psi }_p)\right)+\left(C_{\psi }+{\tilde{C}}_{\psi }\right)\end{array}\right.$$

(5)

The corresponding acceleration term can be converted into

$$\left\{\begin{array}{l}{C}_{\theta }=-2r{\omega }_{e}sin\psi cos\phi /\left(vcos\theta cos(\psi -{\psi }_p)\right)\\ C_{\theta }^{\prime }=-{\omega }_e^2{r}^2(cos\phi sin\phi cos\psi sin\theta +{cos}^2\varphi cos\theta )/\left(v^{2}cos\theta cos(\psi -{\psi }_p)\right)\\ C_{\psi }=-2r{\omega }_e(sin\phi -cos\phi tan\theta cos\psi )/\left(vcos\theta cos(\psi -{\psi }_p)\right)\\ C_{\psi }^{\prime }=-{\omega }_e^2{r}^{2}sin\phi cos\phi sin\psi /v^2{cos}^2\theta cos(\psi -{\psi }_p)\end{array}\right.$$

(6)

where ${\omega }_e$ represents the earth rotation constant.

2.2. Constraints Settings

There are process constraints, no-fly zone constraints, and terminal constraints in the flight of the vehicle. Where the terminal constraint is expressed as

$$\left\{\begin{array}{l}r(s=0)=r_f\\ v(s=0)=v_f\\ \varphi (s=0)={\varphi }_f\\ \phi (s=0)={\phi }_f\end{array}\right.$$

(7)

where, $r_f$, $v_f$, ${\varphi }_f$, ${\phi }_f$ represent the terminal altitude, terminal speed and longitude and latitude, respectively.

The process constraint includes the control variables constraint, heat flux constraint ${\dot{Q}}_{\max }$, dynamic pressure constraint $q_{\max }$, and overload constraint $n_{\max }$, of the vehicle. The control variables generally include the attack angle $\alpha$ and bank angle $\beta$ constraints

$$\left\{\begin{array}{l}{\alpha }_{\min }\le \alpha \le {\alpha }_{\max }\\ {\beta }_{\min }\le \beta \le {\beta }_{\max }\end{array}\right.$$

(8)

where, ${\alpha }_{\min },{\alpha }_{\max },{\beta }_{\min },{\beta }_{\max },$ represent the minimum, the maximum value of attack angle, the minimum and the maximum value of the bank angle, respectively.

The heat flow density, dynamic pressure constraint and overload calculation formulas are shown as

$$\left\{\begin{array}{l}\dot{Q}=C{\rho }^{0.5}{(v\cdot v_c)}^{3.15}\le {\dot{Q}}_{\max }\\ q=0.5\rho {(v\cdot v_c)}^2\le q_{\max }\\ n=\sqrt{L^2+D^2}{r}^2\le n_{\max }\end{array}\right.$$

(9)

where, ${\dot{Q}}_{\max }$ represents the maximum value of heat flow density, $q_{\max }$ represents the maximum value of dynamic pressure, $n_{\max }$ represents the maximum value of overload, $C$ represents the aerodynamic heat coefficient.

In this paper, the no-fly zone is set as a circular area, and the re-entry glide vehicle shall not pass through the no-fly zone. It can be expressed as

$${\left(\varphi -{\varphi }_{bi}\right)}^2+{\left(\phi -{\varphi }_{bi}\right)}^2-r_{bi}^2\ge 0$$

(10)

where, ${\varphi }_b$ and ${\phi }_b$ represent the central longitude and latitude of the center of no-fly zone, $r_b$ represent the longitude and latitude radius of the no-fly zone, and $i=1,\dots ,n$ represent the number of no-fly zones.

2.3. Description of Re-Entry Generation Problem in the Flight Range Domain

The purpose of trajectory generation is to realize terminal constraints under several constraints. Thus, the optimization objective is

$$J_0={\kappa }_1\left({({\varphi }_f^s-{\varphi }_f)}^2+{({\phi }_f^s-{\phi }_f)}^2\right)+{\kappa }_2{(r_f^s-r_f)}^2$$

(11)

where, ${\varphi }_f^s,{\phi }_f^s,r_f^s$ represent longitude and latitude and geocentric distance when the projected range-to-go of the vehicle is 0. ${\kappa }_1$, ${\kappa }_2$ are optimization coefficients.

The optimization problem of trajectory generation in the flight range domain can be expressed as

$$\begin{array}{l}{P}_0:\min J_0={\kappa }_1\left({({\varphi }_f^s-{\varphi }_f)}^2+{({\phi }_f^s-{\phi }_f)}^2\right)+{\kappa }_2{(r_f^s-r_f)}^2\\ \mathrm{subject} \mathrm{to}: \left(5\right), \left(8\right), \left(9\right), \left(10\right),v(s=0)=v_f\end{array}$$

(12)

3. Convexification and Discretization of the Re-Entry Trajectory Generation Problem

The optimization problem in the Section 2 is a nonlinear and non-convex problem, which is difficult to solve directly. Therefore, the problem needs to be convex.

3.1. State, Control Variable Settings and Model Convexity

According to the equation dynamic model (5), the control variables are the attack angle and bank angle, while the lift and drag coefficients of the CAV [29], like the vehicle simulated in this paper generally meet the following equations:

$$\left\{\begin{array}{l}{C}_L=c_{l2}\alpha +c_{l1}\\ C_D=c_{d3}{\alpha }^2+c_{d2}\alpha +c_{d1}\end{array}\right.$$

(13)

where, $c_{l1}$, $c_{l2}$ represent the lift parameters, and $c_{d1}$, $c_{d2}$, $c_{d3}$ represent the drag parameters.

Set the state variable $X=[r,\varphi ,\phi ,v,\theta ,\psi ,\alpha ,\beta ]$ and control variable $U=[\dot{\alpha },\dot{\beta }]$ of the vehicle. With the new control quantity, the dynamic model (5) is converted into

$$\left\{\begin{array}{l}dr/ds=-rtan\theta /cos(\psi -{\psi }_p)\\ d\varphi /ds=-sin\psi /\left(cos\phi cos(\psi -{\psi }_p)\right)\\ d\phi /ds=-cos\psi /cos(\psi -{\psi }_p)\\ dv/ds=C_{D}rqs/m{g}_{0}vcos\theta cos(\psi -{\psi }_p)+tan\theta /\left(vrcos(\psi -{\psi }_p)\right)\\ d\theta /ds=-C_{L}cos\beta rqs/\left(m{g}_0{v}^{2}cos\theta cos\left(\psi -{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em} -1/cos\left(\psi -{\psi }_p\right)+1/\left(v^{2}rcos\left(\psi -{\psi }_p\right)\right)+\left(C_{\theta }+{\tilde{C}}_{\theta }\right)\\ d\psi /ds=-C_{L}sin\beta rqs/\left(m{g}_0{v}^{2}cos\theta cos\left(\psi -{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em}\hspace{1em} -sin\psi tan\phi /\left(cos\left(\psi -{\psi }_p\right)\right)+\left(C_{\psi }+{\tilde{C}}_{\psi }\right)\\ d\alpha /ds=\dot{\alpha }\\ d\beta /ds=\dot{\beta }\end{array}\right.$$

(14)

Equation (14) can be expressed as

$$\dot{X}=g(X,S)+BU\left(S\right)+h(X,S)=F\left(X,U,S\right)$$

(15)

$$g\left(X,S\right)=\left[\begin{array}{l}-rtan\theta /cos(\psi -{\psi }_p)\\ -sin\psi /\left(cos\phi cos(\psi -{\psi }_p)\right)\\ -cos\psi /cos(\psi -{\psi }_p)\\ C_{D}rqs/m{g}_{0}vcos\theta cos(\psi -{\psi }_p)+tan\theta /\left(vrcos(\psi -{\psi }_p)\right)\\ -C_{L}cos\beta rqs/m{g}_0{v}^{2}cos\theta cos(\psi -{\psi }_p)-\frac{r}{v^{2}cos(\psi -{\psi }_p)}(\frac{v^2}{r}-\frac{1}{r^2})\\ -C_{L}sin\beta rqs/m{g}_0{v}^2{cos}^2\theta cos(\psi -{\psi }_p)-sin\psi tan\phi /cos(\psi -{\psi }_p)\\ 0\\ 0\end{array}\right]$$

(16)

$$B={\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]}^T$$

(17)

$$h\left(X,S\right)={\left[0\hspace{1em}0\hspace{1em}0\hspace{1em}0\hspace{1em}{C}_{\theta }+C_{\theta }^{\prime }\hspace{1em}{C}_{\psi }+C_{\psi }^{\prime }\hspace{1em}0\hspace{1em}0\hspace{1em}\right]}^T$$

(18)

Since $g\left(X,S\right)$ is a nonlinear function of $X$, further linearization of the dynamic equation is required. The first order Taylor expansion is operated on the last iterative trajectory, which can be expressed as

$$\dot{X}={\left.\frac{\partial g\left(X\right)}{\partial X}\right|}_{x=x^{(k)}}X+BU+g\left(X^{(k)}\right)-{\left.\frac{\partial g\left(X\right)}{\partial X}\right|}_{X=X^{(k)}}{X}^{(k)}+h(X^{\left(k\right)})$$

(19)

Let $A={\left.\frac{\partial g\left(x\right)}{\partial x}\right|}_{x=x^{(k)}}$, $C=g\left(X^{(k)},S\right)-A{X}^{(k)}+h(X^{(k)},S)$, then

$$\dot{X}=AX+BU+C$$

(20)

$$A=\left[\begin{array}{cccccccc}{a}_{11}& 0& 0& 0& a_{15}& a_{16}& 0& 0\\ 0& 0& a_{23}& 0& 0& a_{26}& 0& 0\\ 0& 0& 0& 0& 0& a_{36}& 0& 0\\ a_{41}& 0& 0& a_{44}& a_{45}& a_{46}& a_{47}& 0\\ a_{51}& 0& 0& a_{54}& a_{55}& a_{56}& a_{57}& a_{58}\\ a_{61}& 0& a_{63}& 0& a_{65}& a_{66}& a_{67}& a_{68}\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$

(21)

$$\left\{\begin{array}{l}{a}_{11}=-tan{\theta }^{\left(k\right)}/cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\\ a_{41}=D^{\left(k\right)}\left(1-{\mathrm{R}}_e{r}^{\left(k\right)}/H_0\right)/\left(v^{\left(k\right)}cos{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em} -tan{\theta }^{\left(k\right)}/v^{\left(k\right)}{r}^{\left(k\right)2}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\\ a_{51}=-L^{\left(k\right)}cos{\beta }^{\left(k\right)}\left(1-{\mathrm{R}}_e{r}^{\left(k\right)}/H_0\right)/\left(v^{2\left(k\right)}cos{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em} -1/v^{\left(k\right)2}{r}^{\left(k\right)2}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\\ a_{61}=-L^{\left(k\right)}sin{\beta }^{\left(k\right)}\left(1-{\mathrm{R}}_e{r}^{\left(k\right)}/H_0\right)/\left(v^{2\left(k\right)}cos{\theta }^{2\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\end{array}\right.$$

(22)

$$\left\{\begin{array}{l}{a}_{23}=-sin{\psi }^{(k)}sin{\phi }^{(k)}/\left({cos}^2{\phi }^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\\ a_{63}=-sin{\psi }^{(k)}/\left({cos}^2{\phi }^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\end{array}\right.$$

(23)

$$\left\{\begin{array}{l}{a}_{44}=D^{(k)}r/\left(v^{2(k)}cos{\theta }^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\\ \hspace{1em}\hspace{1em} -tan{\theta }^{(k)}/\left(v^{(k)}{}^2{r}^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\\ a_{54}=-2/\left(r^{(k)}{v}^{(k)3}cos({\psi }^{(k)}-{\psi }_p)\right)\end{array}\right.$$

(24)

$$\left\{\begin{array}{l}{a}_{15}=-r^{\left(k\right)}/\left({cos}^2{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ a_{45}=D^{\left(k\right)}{r}^{\left(k\right)}tan{\theta }^{\left(k\right)}/\left(v^{\left(k\right)}cos{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em} +1/\left({cos}^2{\theta }^{\left(k\right)}{v}^{\left(k\right)}{r}^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ a_{55}=-L^{\left(k\right)}cos{\beta }^{\left(k\right)}{r}^{\left(k\right)}tan{\theta }^{\left(k\right)}/\left(v^{2\left(k\right)}cos{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ a_{65}=-2L^{\left(k\right)}sin{\beta }^{\left(k\right)}{r}^{\left(k\right)}tan{\theta }^{\left(k\right)}/\left(v^{2\left(k\right)}{cos}^2{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\end{array}\right.$$

(25)

$$\left\{\begin{array}{l}{a}_{16}=-tan{\theta }^{\left(k\right)}{r}^{\left(k\right)}sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)/{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)\\ a_{26}=\left(-cos{\psi }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)-sin{\psi }^{\left(k\right)}sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)/\left(cos{\phi }^{\left(k\right)}{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ a_{36}=\left(sin{\psi }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)-cos{\psi }^{\left(k\right)}sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)/{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)\\ a_{46}=D^{\left(k\right)}{r}^{\left(k\right)}tan\left({\psi }^{\left(k\right)}-{\psi }_p\right)/\left(v^{\left(k\right)}cos{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em} +tan{\theta }^{\left(k\right)}sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)/\left(v^{\left(k\right)}{r}^{\left(k\right)}{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ a_{56}=-L^{\left(k\right)}cos{\beta }^{\left(k\right)}{r}^{\left(k\right)}tan\left({\psi }^{\left(k\right)}-{\psi }_p\right)/\left(v^{2\left(k\right)}cos{\theta }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em} -sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)/{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)+sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)/\left(r^{\left(k\right)}{v}^{\left(k\right)2}{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ a_{66}=-L^{\left(k\right)}sin{\beta }^{\left(k\right)}{r}^{\left(k\right)}tan\left({\psi }^{\left(k\right)}-{\psi }_p\right)/\left(v^{2\left(k\right)}cos{\theta }^{2\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)\right)\\ \hspace{1em}\hspace{1em} +\left(-cos{\psi }^{\left(k\right)}cos\left({\psi }^{\left(k\right)}-{\psi }_p\right)tan{\phi }^{\left(k\right)}-sin{\psi }^{\left(k\right)}sin\left({\psi }^{\left(k\right)}-{\psi }_p\right)tan{\phi }^{\left(k\right)}\right)/{cos}^2\left({\psi }^{\left(k\right)}-{\psi }_p\right)\end{array}\right.$$

(26)

$$\left\{\begin{array}{l}{a}_{47}=\left(2c_{d3}{\alpha }^{(k)}+c_{d2}\right)r^{(k)}{q}^{(k)}s/\left(m{g}_0{v}^{(k)}cos{\theta }^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\\ a_{57}=-c_{l2}cos{\beta }^{(k)}{r}^{(k)}{q}^{(k)}s/\left(m{g}_0{v}^{2(k)}cos{\theta }^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\\ a_{67}=-c_{l2}sin{\beta }^{(k)}{r}^{(k)}{q}^{(k)}s/\left(m{g}_0{v}^{2(k)}cos{\theta }^{2(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\end{array}\right.$$

(27)

$$\left\{\begin{array}{l}{a}_{58}=L^{(k)}sin{\beta }^{(k)}{r}^{(k)}/\left(v^{2(k)}cos{\theta }^{(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\\ a_{68}=-L^{(k)}cos{\beta }^{(k)}{r}^{(k)}/\left(v^{2(k)}cos{\theta }^{2(k)}cos({\psi }^{(k)}-{\psi }_p)\right)\end{array}\right.$$

(28)

3.2. Constraint Convexity and Relaxation

(1)

Constraints of state variables

The model needs to establish certain trust region constraints on the latest iterative trajectory to ensure that the deviation is not too large, which can reduce the error between the convex model and the real model.

$$\left\{\left|X-X^{\left(k\right)}\right|\le {\xi }_X\right.$$

(29)

where ${\xi }_X$ is the trust region of state.

(2)

Convexification of process constraint

If the process constraint (8) of the re-entry glide vehicle is expanded on the R-V (H-V) profile

$$\left\{\begin{array}{l}r\ge \left(2\ast H_0\ln \left({\rho }_0{\left(v\ast vc\right)}^{3.15}C/{\dot{Q}}_{\max }\right)+R_e\right)/R_e=r_{\dot{Q}}\\ r\ge \left(H_0\ln \left({\rho }_0{\left(v\ast vc\right)}^2/2q_{\max }\right)+R_e\right)/R_e=r_q\\ r\ge \left(H_0\ln \left(S{\rho }_0{\left(v\ast vc\right)}^2\sqrt{C_D^2+C_L^2}/\left(2n_{\max }m{g}_0\right)\right)+R_e\right)/R_e=r_n\end{array}\right.$$

(30)

where, $r_{\dot{Q}},r_q,r_n$ represents thermal density, dynamic pressure and overload height.

Expand Equation (30) in the last iterative trajectory

$$\left\{\begin{array}{l}{r}_{\dot{Q}}=r_{\dot{Q}}^{(k)}+{\left.\frac{\partial r_{\dot{Q}}}{\partial v}\right|}_{X^{(k)}}(v-v^{(k)})\\ r_q=r_q^{(k)}+{\left.\frac{\partial r_q}{\partial v}\right|}_{X^{(k)}}(v-v^{(k)})\\ r_n=r_n^{(k)}+{\left.\frac{\partial r_n}{\partial v}\right|}_{X^{(k)}}(v-v^{(k)})+{\left.\frac{\partial r_{nj}}{\partial \alpha }\right|}_{X^{(k)}}(u_3-u_3^{(k)})\end{array}\right.$$

(31)

$$\left\{\begin{array}{l}{\left.\frac{\partial r_{\dot{Q}}}{\partial v}\right|}_{X^{(k)}}=6.3H_0/\left(R_e{v}^{(k)}\right)\\ {\left.\frac{\partial r_q}{\partial v}\right|}_{X^{(k)}}=2H_0/\left(R_e{v}^{(k)}\right)\\ \left[{\left.\frac{\partial r_n}{\partial v}\right|}_{x^{(k)}},{\left.\frac{\partial r_n}{\partial \alpha }\right|}_{x^{(k)}}\right]=\left[2H_0/\left(R_e{v}^{(k)}\right),\frac{H_0}{R_{\mathrm{e}}}\frac{C_D^{(k)}(2c{d}_3{\alpha }^{(k)}+c{d}_2)+C_L^{(k)}c{l}_2}{C_D^{2(k)}+C_L^{2(k)}}\right]\end{array}\right.$$

(32)

where $H_0$ represents the air density calculation constant.

(3)

Convexification of the no-fly zones constraint

The no-fly zones constraint (9) is non-convex and cannot be solved directly and effectively, so it needs to be convex.

$$\begin{array}{l}{\left({\varphi }^{(k)}-{\varphi }_{bi}\right)}^2+{\left({\phi }^{(k)}-{\varphi }_{bi}\right)}^2+2\left({\varphi }^{(k)}-{\varphi }_{bi}\right)\left(\varphi -{\varphi }^{(k)}\right)\\ +2\left({\phi }^{(k)}-{\phi }_{bi}\right)\left(\phi -{\phi }^{(k)}\right)-r_{bi}^2\ge 0\end{array}$$

(33)

3.3. Improved Description of Convex Optimization Problem

According to the settings in Section 3.1 and Section 3.2, the new optimization problem can be transformed into

$$J_1=J_0-C_3{\displaystyle {\int }_{s_f}^0\left|\dot{\alpha }\right|+\left|\dot{\beta }\right|dS}$$

(34)

Then convert problem $P_0$ into problem $P_1$ through convexity.

$$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{P}_1:\min \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} J_1\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}subject to:\hspace{1em}\hspace{1em} (20),(29),(31),(33)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} X\left(S_{f0}\right)=[r_0,{\varphi }_0,{\phi }_0,v_0,{\theta }_0,{\psi }_0], v(s=0)=v_f\end{array}$$

(35)

3.4. Discretization of Dynamic Model

The problem $P_1$ is a representation in the continuous flight range domain, which is a model with infinite dimension. In order to further solve the problem, it needs to be discretized. First, set grid point locations of the generated trajectory in the flight range domain as $S_i$, $i=1,\dots ,m$ represents the serial number of the grid point, $m$ represents the number of the grid points, $S_1=S_{f0}\ge S_{i+1}\ge S_i\ge S_m=0$.

Then the dynamic Equation (20) can be discretized using the Euler method.

$$X_{i+1}=X_i+\Delta S_i{\dot{X}}_i=X_i+\Delta S_i\left(A_i^{(k)}{X}_i+B_i^{(k)}{U}_i+C_i^{(k)}\right)$$

(36)

where, $\Delta S_i$ represents the projected range-to-go difference between grid point $i$ and $i+1$.

At the same time, due to different iterations, the trust region size required of optimization may be different, so the variable trust region method is adopted in Equation (30) to improve the search range of solution

$$\left\{\left|X_j-X_j^{(k)}\right|\le \eta {\xi }_x\right.$$

(37)

where $\eta$ represents variable of trust region

On this basis, the optimization objective $J_1$ is further discretized to obtain a new optimization objective $J_2$

$$J_2=J_0-{\kappa }_3\left({\displaystyle \sum _{j=0}^N\left|{\dot{\alpha }}_j\right|\Delta S_j}+{\displaystyle \sum _{j=0}^N\left|{\dot{\beta }}_j\right|\Delta S_j}\right)+{\kappa }_4\eta$$

(38)

where ${\kappa }_3$, ${\kappa }_4$ are optimization coefficients.

Then discretize the constraints (31), (33) to further obtain the final optimization problem $P_2$

$$\begin{array}{l}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{P}_2:\min \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{J}_2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}subject to:\hspace{1em}\hspace{1em}\hspace{1em}(31),(32),(36),(37)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}X\left(S_{f0}\right)=[r_0,{\varphi }_0,{\phi }_0,v_0,{\theta }_0,{\psi }_0], v(s=0)=v_f\end{array}$$

(39)

3.5. Termination Condition of Solving

On the basis of problem $P_2$, when the convexity value of each grid point for two consecutive times iteration meets the convergence threshold $\epsilon$, the convexity iteration is terminated. The iteration termination conditions are as follows

$$\left\{\begin{array}{l}\underset{1\le i\le n}{\max }\left|r^{\left(k+1\right)}\left(S_i\right)-r^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_r,\underset{1\le i\le n}{\max }\left|{\varphi }^{\left(k+1\right)}\left(S_i\right)-{\varphi }^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_{\varphi }\\ \underset{1\le i\le n}{\max }\left|{\phi }^{\left(k+1\right)}\left(S_i\right)-{\phi }^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_{\phi },\underset{1\le i\le n}{\max }\left|v^{\left(k+1\right)}\left(S_i\right)-v^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_v\\ \underset{1\le i\le n}{\max }\left|{\theta }^{\left(k+1\right)}\left(S_i\right)-{\theta }^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_{\theta },\underset{1\le i\le n}{\max }\left|{\psi }^{\left(k+1\right)}\left(S_i\right)-{\psi }^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_{\psi }\\ \underset{1\le i\le n}{\max }\left|{\alpha }^{\left(k+1\right)}\left(S_i\right)-{\alpha }^{\left(k\right)}\left(S_j\right)\right|\le {\epsilon }_{\alpha },\underset{1\le i\le n}{\max }\left|{\beta }^{\left(k+1\right)}\left(S_i\right)-{\beta }^{\left(k\right)}\left(S_i\right)\right|\le {\epsilon }_{\beta }\end{array}\right.$$

(40)

where $\left[{\epsilon }_r,{\epsilon }_{\varphi },{\epsilon }_{\phi },{\epsilon }_{\nu },{\epsilon }_{\theta },{\epsilon }_{\psi },{\epsilon }_{\alpha },{\epsilon }_{\beta }\right]$ represent the convergence threshold of the corresponding state and control variables.

4. Fast Initial Trajectory Setting and Distributed Grid Points Adjustment

4.1. Fast Initial Trajectory Setting

In order to improve the solution seeking speed of the method and ensure that the initial trajectory meets a large number of constraints as much as possible, according to the landing points of the vehicle, a fast initial trajectory setting method that meets the dynamic models was designed. Since the terms of Coriolis acceleration and implicated acceleration caused by the earth rotation are small, the flight path of the re-entry glide vehicle in any direction, at any longitude and latitude coordinate can be approximated to the equatorial origin, and with the same initial speed and same control law in the longitudinal direction by rotating and moving, which is shown in Figure 2a,b. In this paper, the landing points map of the constant angle of attack and bank angle control strategy of the re-entry glide vehicle is established offline, and suitable control variables and initial mapping trajectory are obtained by interpolation, as shown in Figure 2b. Finally, the initial generated trajectory satisfying the constraint can be obtained by integrating the corresponding control law in the original state.

Reachable landing point locations of the re-entry glide vehicle are shown in Figure 3a,b. Among them, the change of constant attack angle and constant bank angle has regularity. Generally, with the fixed attack angle and the different bank angle, the lateral landing points locations are mainly affected. With the same bank angle and different attack angle, the flight range is affected. When the setting of the attack angle is less than the maximum lift drag ratio attack angle, the flight range is proportional to the attack angle; when it is greater than the maximum lift drag ratio attack angle, the flight range is inversely proportional to the attack angle.

4.2. Distributed Grid Points Adjustment

According to the literature [28], with the last iterative grid points, the grid points can be effectively adjusted by calculating the non-linear illegal degree of the grid points. The non-linear illegal degree of the grid point $i$ is

$$\Delta n_i={\displaystyle \sum _{i=2}^m{‖X^{(k)}\left(S_i\right)-{\displaystyle {\int }_{S_{i-1}}^{S_i}F\left(X,U,S\right)dS}‖}_2}$$

(41)

Let the number of grid point intervals be $K$, $\Delta S_j=(S_i-S_{j-1})/K$

$${\displaystyle {\int }_{S_{i-1}}^{S_i}F\left(X,U,S\right)dS}={\displaystyle \sum _{j=1}^{K}F\left(X,U,S_i+j\Delta s_j\right)}$$

(42)

In order to adjust grid points more quickly and dynamically, this paper defines the estimated nonlinear illegal degree distribution probability density function $P$, and defines increasing the value threshold of the grid points’ nonlinear illegal degree $\Delta n_{max}$. When the nonlinear illegal degree of a point exceeds this value, it means that the nonlinearity illegal degree of this point is high, and grid points need to be increased; At the same time, the deleting value threshold $\Delta n_{\min }$ is also defined. The grid point supplement function $P^{\prime }$ of a point is expressed as

$$P^{\prime }(S_i)=\left\{\begin{array}{l}0,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}if \Delta n_i\le \Delta n_{\min }\\ 1,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}if \Delta n_{\min }\le \Delta n_i\le \Delta n_{\max } or i=1 or i=n\\ \Delta n_i/\Delta n_{\max }, if \Delta n_i\ge \Delta n_{\max }\end{array}\right.$$

(43)

At the same time, in order to ensure the iteration speed of the method, set the upper number of grid points $n_{\max }$, and the adjusted number of grid points $n^{(k+1)}$ meets

$$n^{(k+1)}=\left\{\begin{array}{l}{\displaystyle \sum _{i=1}^n{P}^{\prime }(S_i)},\hspace{1em}if {\displaystyle \sum _{i=1}^n{P}^{\prime }(S_i)}\le n_{\max }\\ n_{\max },\hspace{1em}if {\displaystyle \sum _{i=1}^n{P}^{\prime }(S_i)}>n_{\max }\end{array}\right.$$

(44)

The number of grid points to be adjusted is $\Delta n^{(k+1)}=Count\left(P^{\prime }=0\right)+n^{(k+1)}-n^{(k)}$.

Let the probability density function of the adjusted grid points be $P$, and it can be expressed as

$$P(S_i)=P^{\prime }(S_i)/\left({\displaystyle \sum _{i=1}^n{P}^{\prime }(S_i)}\left|\Delta S_i\right|\right)$$

(45)

For any distance domain location $S$, the probability density of this point can be obtained by interpolating by the calculated value of Equation (44). Finally, the corresponding cumulative distribution function is inversely solved by (45) to obtain the location, when the corresponding cumulative probability value is equal to $i/\left(\Delta n^{(k+1)}+1\right),i=1,\dots ,n^{(k+1)}$ as the additional grid point location.

As shown in Figure 4, it is assumed that after the generation of initial trajectory, the grid points will be a uniform distribution according to the flight range domain in the first iteration. After the first iteration, the probability density function distribution of grid points and the degree of nonlinear illegal degree will be calculated. In the first iteration, the red point location is of lower nonlinear illegal degree, so it can be deleted for adjustment. According to the probability density function distribution, the adjustment position is the blue dotted point of the second iteration; After the second optimization, find the location where the degree of nonlinear illegal is lower again, and adjust them to the new grid points location according to the second probability distribution, and continue to the third iteration.

5. Simulation

5.1. Experimental Subjects and Parameter Settings

In this paper, CAV-H is selected as the simulation vehicle, and its specific parameters can be referred to in the literature [29], where its process constraints are set ${\dot{Q}}_{\max }=4\times {\mathrm{10}}^6 \mathrm{W}/{\mathrm{m}}^2$, $q_{\max }=100 \mathrm{Kpa}$, $n_{\max }=4$. The trust region $[{\xi }_r,{\xi }_{\varphi },{\xi }_{\phi },{\xi }_v,{\xi }_{\theta },$${\xi }_{\psi },{\xi }_{\alpha },{\xi }_{\beta }]$ are $[2\times {\mathrm{10}}^4/R_e,500\times {\mathrm{10}}^4\pi /R_e,500\times {\mathrm{10}}^4\pi /R_e,1000/v_c,20\pi /180,30\pi /180,10\pi /180,$$90\pi /180]$, convergence error threshold $\left[{\epsilon }_r,{\epsilon }_{\varphi },{\epsilon }_{\phi },{\epsilon }_v,{\epsilon }_{\theta },{\epsilon }_{\psi },{\epsilon }_{\alpha },{\epsilon }_{\beta }\right]$ are $[200/R_e,5\times {\mathrm{10}}^4\pi /R_e,5\times {\mathrm{10}}^4\pi /R_e,$$100/v_c,1\pi /180,5\pi /180,1\pi /180,5\pi /180]$, maximum number of grid points $n_{\max }=200$. The simulation software is python with Win10, deploying ECOS-BB solver for sequential programming and Matlab for plotting. The simulation CPU is Inter-I7.

5.2. Verification of the Initial Trajectory Generation

In order to verify the effectiveness of the initial trajectory generation method proposed in this paper, the three types of relationship between initial point and target point are selected to be tested, respectively, to analyze the guidance effects and compare with the iterative trajectories. Where the initial states of all trajectories are set as $[r_0,{\varphi }_0,{\phi }_0,v_0,{\theta }_0,{\psi }_0]=[\left({\mathrm{R}}_e+70 \mathrm{km}\right)/R_e,0^{°},0^{°},6000/v_c,-2^{°},{45}^{°}]$. Then, three guidance location are selected to generate trajectories on different sides in the velocity direction of the initial state $[r_{f1},{\varphi }_{f1},{\phi }_{f1}]=[\left(R_e+30 \mathrm{km}\right)/R_e,{30}^{\mathrm{o}},{15}^{\mathrm{o}}]$, $[r_{f2},{\varphi }_{f2},{\phi }_{f2}]=$$[\left(R_e+30 \mathrm{km}\right)/R_e,{34}^{\mathrm{o}},{24}^{\mathrm{o}}]$, $[r_{f3},{\varphi }_{f3},{\phi }_{f3}]=[\left(R_e+30 \mathrm{km}\right)/R_e,{24}^{\mathrm{o}},{34}^{\mathrm{o}}]$. It can be seen from Figure 5a,b that the longitudinal and lateral landing point of the initial trajectory can basically reach the target point by the initial trajectory generation algorithm in this paper. Then, target point 2 was taken as the example to compare the iterative trajectory. As can be seen from Figure 5c,d, by our method the vehicle could effectively converge to the optimal trajectory after solving for five iterations with a total time of 18.46 s. When the maximum attack angle strategy method was adopted as the initial trajectory in this paper, the number of iterations of the algorithm increased to nine. The total time is 40.21 s, which indicates our initial trajectory generation method has a better effectiveness for the optimization algorithm.

5.3. Validity Verification

In order to verify the effective generation of the re-entry trajectory under constraints by the proposed method, three cases are selected for verification. The initial states of the vehicle are the same in all three cases. $[r_0,{\varphi }_0,{\phi }_0,v_0,{\theta }_0,{\psi }_0]=[\left({\mathrm{R}}_e+70 \mathrm{km}\right)/R_e,0^{°},0^{°},6000/v_c,0^{°},{60}^{°}]$ when the terminal states are also the same $[r_f,{\varphi }_f,{\phi }_f,v_f]=[\left(R_e+30 \mathrm{km}\right)/R_e,{30}^{°},{25}^{°},2000/v_c]$. In the three cases, the vehicle will have a flight with the no no-fly zone, flying- around the no-fly zones and passing through the no-fly zones in three cases. In case 2,3, the no-fly zones are set as follows:

$$\mathrm{Case} 2: \left[{\varphi }_{b1},{\phi }_{b1},r_{b1}\right]=\left[{15}^{\mathrm{o}},{8.5}^{\mathrm{o}},2^{\mathrm{o}}\right], \left[{\varphi }_{b1},{\phi }_{b1},r_{b1}\right]=\left[{24}^{\mathrm{o}},{16}^{\mathrm{o}},2^{\mathrm{o}}\right]$$

$$\mathrm{Case} 3:\left[{\varphi }_{b1},{\phi }_{b1},r_{b1}\right]=\left[{15}^{\mathrm{o}},{11}^{\mathrm{o}},2^{\mathrm{o}}\right], \left[{\varphi }_{b1},{\phi }_{b1},r_{b1}\right]=\left[{24}^{\mathrm{o}},{14}^{\mathrm{o}},2^{\mathrm{o}}\right]$$

The simulation effect is shown in Figure 6, Figure 7 and Figure 8. As can be seen from Figure 6, the convexation algorithm proposed in this paper can realize effective trajectory generation after a variety of situations for all three cases. As can be seen from case 1, although the longitudinal guidance effect of the initial trajectory is not very good, the method can realize accurate lateral and longitudinal trajectory generation. On the other hand, it can be seen from cases 2 and 3 that although there is crossover between the initial trajectory and the no-fly zones, the method can avoid the no-fly zone through effective adjustment and satisfy the longitudinal guidance at the same time.

On the other hand, according to the effect of Figure 6, Figure 7 and Figure 8, all the trajectories can meet the final constraint requirements. In Case 2 and Case 3, the peak values of several constraints are high due to the avoidance of the no-fly zones, especially at low flight altitudes. In order to effectively meet the constraints, at the initial stage of flight, the vehicle usually employs a high attack angle. In case 3, in order to pass through and between the no-fly zones and increase the range, a smaller attack angle value is adopted; The change of the corresponding bank angle represents the avoidance maneuver process of the vehicle. In Case 1, there are no no-fly zones, so the change of the bank angle is small; However, in Case 2 and Case 3, considering the influence of no-fly zones, there are many variation of the bank angle; this can be seen from the grid point changes when the uniform grid points are used in the initial trajectory. For the last iteration trajectory, due to the small aerodynamic force in the initial stage of vehicle, where the nonlinear illegal degree is low, and the grid points are relatively sparse; at the end of flight, the flight height is low, the aerodynamic force changes greatly, and the grid points are obviously dense.

The optimized trajectories’ accuracy, iteration times and the waste time of the three cases are shown in

Table 1. Generation effect with different cases.

Case	${\mathit{E}}_{\mathit{l}} \left(\mathbf{km}\right)$	${\mathit{E}}_{\mathit{h}} \left(\mathbf{km}\right)$	Waste Time (s)	Iteration
1	0.001	−6.43 × 10−4	22.93	5
2	0.041	0.005	27.20	5
3	0.016	−0.001	110.67	19

. Where $E_l$ represents the lateral error and $E_h$ represents the longitudinal error. There is a no-fly zone in Case 1, so its lateral and longitudinal trajectories generation accuracy is very high; due to the influence of the no-fly zone, Case 2 and Case 3 have poor guidance accuracy, but the lateral error is not more than 1.6 km, and the longitudinal error is not more than 100 m. For the number of iteration convergence, Case 3 needs to pass through and between the no-fly zones, so it has many iterations and takes a longer time, but the average iteration time is not more than 6 s.

5.4. Simulation Experiment with Different Target Points

The simulations are operated with the same initial and final states and different no-fly zones. The simulated results of different target points are analyzed in this section. The target points and simulation results are shown in

Table 2. Generation effect with different target points.

Case	${\mathit{\varphi }}_{\mathit{f}}{(}^{\mathbf{o}})$	${\mathit{\phi }}_{\mathit{f}}{(}^{\mathbf{o}})$	${\mathit{E}}_{\mathit{l}} \left(\mathbf{km}\right)$	${\mathit{E}}_{\mathit{h}} \left(\mathbf{km}\right)$	Waste Time (s)	Iteration
1	28	23	0.167	0.009	135.84	23
2	30	20	0.082	−0.002	38.63	7
3	30	10	0.079	−0.009	34.87	6
4	27	27	0.139	−0.019	122.70	21

. The setting of the no-fly zones is the same as that in Case 3 in Section 5.3. The simulation effect is shown in Figure 9. It can be seen from Figure 9 that for different target points, this method can realize lateral and longitudinal guidance with a high accuracy. When effective avoidance is achieved, reasonable control angles are adopted. Among them, a trajectory with a long range and requiring the avoidance of no-fly zones has more iterations, relatively poor accuracy and more frequent changes in control angles.

5.5. Simulation Experiment with Initial State Disturbance

To further demonstrate the effectiveness and robustness of our method, the simulation is operated when the initial state of the vehicle is disturbed with the same target point and the no-fly zones. The target point and no-fly zones are set as the same as Case 2 in Section 5.3, the initial altitude $\Delta h_0=+2 \mathrm{km}$, initial speed $\Delta v_0=-100 \mathrm{m}/\mathrm{s}$, initial flight path angle $\Delta {\theta }_0=+{0.5}^{\mathrm{o}}$ and initial flight heading angle $\Delta {\psi }_0=-{10}^{\mathrm{o}}$ are, respectively, simulated. The optimization effect is shown in Figure 10 and

Table 3. Generation effect with disturbance of the initial state.

Case	${\mathit{E}}_{\mathit{l}} \left(\mathbf{km}\right)$	${\mathit{E}}_{\mathit{h}} \left(\mathbf{km}\right)$	Waste Time (s)	Iteration
$\Delta h_0=+2 \mathrm{km}$	0.002	−0.001	24.92	5
$\Delta v_0=-100 \mathrm{m}/s$	0.005	−9 × 10−4	19.63	4
$\Delta {\theta }_0=+{0.5}^{\mathrm{o}}$	0.002	−0.002	25.08	5
$\Delta {\psi }_0=-{10}^{\mathrm{o}}$	0.004	−6 × 10−4	45.54	9

. It can be seen from Figure 10 that with certain disturbance of the initial state, the vehicle can still effectively avoid the no-fly zone and reach the target points with high accuracy in both the lateral and longitudinal directions. The control angle obtained by changing the height, speed and the flight angle is similar, which has a certain impact on the longitudinal trajectory. When the initial flight heading angle is changed, the optimization trajectory changes obviously, and the number of iterations in this case increases significantly. This is because the adjusted initial flight heading angle is larger than the LOS angle. In general, the proposed method is robust to initial state disturbances.

5.6. Monte Carlo Robustness Simulation

In order to further reflect the robust effect of the method proposed in this paper on process disturbance, Monte Carlo simulations are operated 100 times with the same conditions of Case 2 in Section 5.2, where the lift and drag coefficient, vehicle mass and reference area have 5% random error each time, and the atmospheric density has 20% random error, which are used to verify the robustness of our method. The results are shown in Figure 11. As can be seen in Figure 11, the convexity trajectory generation algorithm proposed in this paper is very accurate in lateral and longitudinal accuracy in the Monte Carlo test, which proves the robustness of our method with process disturbance conditions.

5.7. Comparison of Mainstream Methods

In order to reflect the fast-planning speed of the method proposed in this paper, the method in this paper is compared with the Gauss pseudo spectral method. The optimization results of different methods are shown in Figure 12 and

Table 4. Generation effect of different methods.

Method	${\mathit{E}}_{\mathit{l}} \left(\mathbf{km}\right)$	${\mathit{E}}_{\mathit{h}} \left(\mathbf{km}\right)$	Waste Time (s)	Iteration	Number of Grid Points
Our method	0.016	0.001	110.67	21	200
Gaussian Spectral method	0.010	1 × 10−4	154.20	25	341

. The trajectory obtained by using the Gauss pseudo spectral method, and the method in this paper, can effectively realize trajectory generation, and meet various constraints when the accuracy of the method is comparable. The trajectory generated by the Gauss pseudo spectral method is shorter, but the control amount changed greatly. In terms of the number of grid points, the method in this paper uses 200 grid points due to the upper limit of grid points, while the Gaussian pseudo spectral method finally uses 341 grid points; in terms of planning time, the method in this paper takes 110.67 s, the Gauss pseudo spectral method takes 154.2 s, and the programing time is longer.

Through several simulations in this section, the effectiveness and robustness of our method have been verified. But the method of changing the maximum number of gird points adaptively could be researched in the future. Inspired by the multi-agent system guidance and control technology [30,31], convex trajectory generation algorithm of multi-re-entry glide targets will also be researched in our future work.

6. Conclusions

In this paper, according to the guidance mechanism of the vehicle, the concept of the projected range-to-go of the vehicle is first proposed. Then, the dynamic model is converted to the flight range domain when the model is convexated and discretized. In order to improve the generation speed and accuracy, a fast and accurate initial trajectory-generating method is proposed according to the landing points of the vehicle under different control laws. According to the nonlinear illegal degree of iterative trajectory, a grid point probability density function is proposed to dynamically change grid points. Through the simulation experiments against various disturbances, all of our final guidance errors are less than 1 km, and waste times are less than 135.62 s, which prove the effectiveness and robustness of our method. Compared with the conventional Gaussian pseudo-spectral method, we can obtain comparable accuracy and a faster generation speed.