Rapid Algorithm for Generating Entry Landing Footprints Satisfying the No-Fly Zone Constraint

Abstract

An online estimation algorithm of landing footprints based on the drag acceleration-energy profile is proposed for an entry hypersonic vehicle. Firstly, based on the Evolved Acceleration Guidance Logic for Entry (EAGLE), drag acceleration-energy profiles are designed. To track the drag acceleration-energy profile obtained by the interpolation, a drag acceleration tracking law is designed. Secondly, based on the constraint model of the no-fly zone, flying around strategies are proposed for different conditions, and a reachable area algorithm is designed for no-fly zones. Additionally, by interpolating the minimum and maximum drag acceleration profiles, the terminal heading angle constraint is designed to realize the accurate calculation of the minimum and maximum downrange ranges by adjusting the sign of the bank angle. In this way, the distribution of landing footprints is more reasonable, and the boundary of a reachable area is more accurate. The simulation results under typical conditions indicate that the proposed method can calculate landing footprints for different situations rapidly and with the good adaptability.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Landing footprints provide key information for the mission planning, such as the terminal area energy management (TAEM) guidance for a shuttle entry [1]. Due to a vehicle’s high lift-drag ratio and the strong maneuverability, it forms a large range of landing footprints. The reachable area [2] formed by landing footprints is the maximum maneuvering range under given initial conditions, path constraints, and terminal constraints. During the entry process, the environment is complex, the flight state varies widely, and the shapes and sizes of reachable areas differ. Different from the case in a conventional mission planning, the initial condition has a great influence on a reachable area for an entry vehicle, and it has a great uncertainty. It is necessary to quickly obtain the reachable area based on landing footprints to select an appropriate landing site. Therefore, the reachable area needs to be calculated online, which requires the rapidity.

Hypersonic entry trajectory planning technology is the foundation of solving the reachable area. The main purpose of the entry trajectory planning system is to steer the vehicle flying from an entry point to a desired terminal area energy management interface. Since a well-planned trajectory is the key to the entry mission, a large amount of research on developing trajectory planning methods can be found. The sequential convex programming (SCP) method can be used to generate a reference trajectory under multiple constraints. Wang [3–5] used a convex quadratically constrained quadratic programming (QCQP) method to generate the optimal tracking guidance command. With the consideration of probabilistic constraints, a differentiable chance-constraint approximation method [6, 7] is proposed. Therefore, the transformed optimal control model becomes solvable for trajectory optimization methods. Based on the optimal trajectories generated off-line by an indirect method, a deep neural network (DNN) model [8–11] is designed to enable the capability of optimal action predictions. Lu et al. [12] and Cheng et al. [13] introduced the Full Numerical Predictor-corrector Entry Guidance (FNPEG) into entry missions without the need for a reference trajectory or the vehicle-dependent planning. Chai et al. [14] proposed an optimal trajectory generation method by constructing an integrated framework with an inner gradient solver.

As the missions become more and more diversified, the vehicle needs to avoid threat areas and politically forbidden areas [15]. Thus, the vehicle needs to consider the constraint of the no-fly zone. In the trajectory optimization, the constraint of the no-fly zone has been studied [16–19]. Yang et al. [20] used an additional bank reversal scheme to shape the entry trajectory at the right time to avoid no-fly zones. Gao et al. [21] improved the tentacle-based bank angle transient lateral strategy to avoid static and dynamic no-fly zones. Yu et al. [22, 23] proposed an analytical iteration scheme to plan the bank reversal sequence, for constraints of multiple no-fly zones and the flight-time. They also employed a new crossrange formula to the schedule proper bank reversal under the no-fly zone and terminal position constraints.

For a hypersonic vehicle, all possible trajectories can be obtained by trajectory optimization methods, and the reachable area can be formed. A well-designed solution to the reachable area should possess characteristics of the constraint satisfaction, the reliability, and the robustness. The existing reachable area algorithms can be divided into three categories: (1) Trajectory optimization techniques are employed directly in landing footprint generation algorithms, generally by a pseudospectral method (PSM) [24, 25]. Fahroo et al. [26, 27] employed the Legendre pseudospectral method for landing footprints. Other trajectory optimization methods such as the Gauss pseudospectral method [28–30] can also be applied to the landing footprint problem. The pseudospectral method has high accuracy, but the calculation takes a long time. (2) A simplified mathematical model based on the quasi-equilibrium glide condition (QEGC) can also generate landing footprints. Lu et al. [31, 32] proposed a closest-approach method to obtain the maximum crossrange trajectory at a particular downrange by setting up a nonreachable virtual target. The centerpiece of this method is a closed-form near-optimal bank angle control law that can solve landing footprints rapidly and reliably. However, this method is limited by a narrow model based on the assumption of the QEGC. Based on the same model, Liu et al. [33] obtained the maximum crossrange trajectory by the convex optimization. Li et al. [34] also proposed the selection scheme of a virtual target to make the method more applicable. Ngo and Doman [35] solved landing footprints under the failure of the control components based on the QEGC. The method only can obtain landing footprints corresponding to the maximum boundary due to the limitations of the simplified mathematical model of the quasi-equilibrium glide condition, and it takes a long time to discern the parameters. (3) Based on the EAGLE and the core idea that the trajectory range is inversely proportional to the drag acceleration, the entry corridor described by the drag acceleration profile can be designed. The maximum and minimum range trajectories and corresponding drag acceleration profiles can be obtained, respectively. Different drag acceleration profiles and different range trajectories can be obtained by linear interpolation. By keeping the sign of the bank angle positive or negative during a whole flight, Saraf et al. [36] can obtain terminal landing footprints. The method requires less calculation time, but its accuracy is lower than that of the pseudospectral method.

To solve landing footprints, a rapid algorithm for generating entry landing footprints satisfying the no-fly zone constraint is designed based on [36]. In this paper, we focus on the following three aspects of contribution: (1) we designed the drag acceleration profile tracking law and introduced the terminal heading angle constraint to realize the accurate calculation of initial landing footprints. (2) Considering the constraint of the no-fly zone, we proposed a fly-around strategy. Different relative trajectories were transformed into corresponding terminal heading angle constraints to realize the fast calculation of reachable areas under the constraint of the no-fly zone. (3) Based on the typical nonlinear model of an entry vehicle, the simulation of the proposed method was carried out. The results show that the method can obtain landing footprints based on the direct physical concept under the no-fly zone constraint. Compared with the traditional trajectory optimization, the calculation speed is faster and satisfies the requirement of the online calculation.

Although entry trajectory optimization problems for a hypersonic vehicle have been significantly addressed in the past years, the reachable area considered in this paper focuses on the generation algorithm of multiple trajectories. Different from [37] in which the no-fly zone constraint is transformed into the violation degree for the trajectory optimization, a fly-around strategy for the no-fly zone constraint is designed by adjusting the sign of the bank angle to avoid the threaten area in this paper. Compared with the reachable area algorithm [36], the no-fly zone constraint is added to make the application of the algorithm closer to the actual scene. In addition, the reachable area algorithm without the no-fly zone constraint is improved to distribute landing footprints more reasonable. Therefore, the results of the reachable area are more accurate compared with previous studies.

This paper is organized as follows: The landing footprint problem is briefly described, and the simulation model is introduced in Section 2. The solution to landing footprints is designed in Section 3. The verification of the simulation model is presented in Section 4. Thereafter, this paper closes with a summary of the most important conclusions in Section 5.

2. Problem Formulation

2.1. Entry Dynamics

Since the time and the energy of an unpowered hypersonic entry vehicle are monotonically changing, they can be chosen as independent variables. Due to constraints of the terminal velocity and altitude, the initial energy and terminal energy are limited in a small range. Therefore, the drag acceleration-energy profile is designed for the trajectory planning in the following sections in this paper. If the time is chosen as the independent variable, the velocity can be calculated by integrating the velocity differential equation. However, considering the consistency with the drag acceleration-energy profile, the inverse mechanical energy $E$ is an appropriate independent variable, which is determined by $E=\mu /r-V^2/2$, where $r$ is the radial distance from the Earth’s center to the vehicle, $V$ is the velocity of the vehicle, and $\mu$ is the gravitational constant designed as 3.98606$e$14. Based on a sphere and irrational earth mode, the 3-dimensional equations of motion [38] are shown as follows: $$\begin{array}{}\text{(1)}& {\theta }^{\prime }=-\frac{\cos \gamma \cos \psi }{r\cos \varphi D},\text{(2)}& {\varphi }^{\prime }=-\frac{\cos \gamma \sin \psi }{rD},\text{(3)}& r^{\prime }=-\frac{\sin \gamma }{D},\text{(4)}& {\psi }^{\prime }=\frac{\cos \psi \tan \varphi \cos \gamma }{rD}+\frac{L\sin \sigma }{V^2\cos \gamma D},\text{(5)}& {\gamma }^{\prime }=\frac{g-V^2/r\cos \gamma }{V^{2}D}-\frac{L\cos \sigma }{V^{2}D},\end{array}$$where $\theta$ is the longitude, $\varphi$ is the latitude, $\psi$ is the heading angle, $\gamma$ is the flight path angle, and $\sigma$ is the bank angle.

The acceleration of gravity $g$ is modeled as follows: $$\begin{array}{}\text{(6)}& g=\frac{\mu }{r^2}.\end{array}$$

The aerodynamic lift acceleration $L$ and the aerodynamic drag acceleration $D$ are given by $$\begin{array}{}\text{(7)}& L=\frac{{qV}^2\cdot S\cdot C_L\alpha ,Ma}{m},\text{(8)}& D=\frac{{qV}^2\cdot S\cdot C_D\alpha ,Ma}{m},\end{array}$$where $q=\rho V^2/2$ is the dynamic pressure, $S$ is the reference area of the vehicle, $C_L$ is the aerodynamic lift coefficient, $C_D$ is the aerodynamic drag coefficient, $\alpha$ is the angle of attack, $Ma=V/c$ is the Mach number, $c$ is the velocity of sound, and $m$ is the mass of the vehicle. For the entry trajectory planning, the angle of attack is usually specified as a function of the Mach number.

The atmospheric density $\rho$ based on the exponential atmosphere model is given by $$\begin{array}{}\text{(9)}& \rho ={\rho }_0{e}^{-h/h_s},\end{array}$$where ${\rho }_0=1.225$ kg/m³ is the atmospheric density at sea level, $h=r-r_0$ is the altitude, $r_0=6378136.3$ m is the equatorial radius of the Earth, and $h_s=7600$ is the atmospheric scale height.

2.2. Constraints

The vehicle has maximum path constraints for the heating rate, dynamic pressure, and aerodynamic acceleration [39], which are given by $$\begin{array}{}\text{(10)}& \dot{Q}=k_Q\sqrt{\rho }{V}^{3.15}\le {\dot{Q}}_{\max },\text{(11)}& q=\frac{\rho V^2}{2\le q_{\max }},\text{(12)}& n=L\cos \sigma +D\sin \alpha \le n_{\max },\end{array}$$where $\dot{Q}$ is the heating rate, $k_Q$ is the thermal model constant, and $n$ is the normal acceleration. Moreover, ${\dot{Q}}_{\max }$, $q_{\max }$, and $n_{\max }$ are the maximum values.

The terminal constraints for the altitude and the velocity are given by $$\begin{array}{}\text{(13)}& h_f-\Delta h\le h{E}_f\le h_f+\Delta h,\text{(14)}& V_f-\Delta V\le V{E}_f\le V_f+\Delta V,\end{array}$$where $\Delta h$ and $\Delta V$ are thresholds of the terminal constraint, $r_f$ and $V_f$ are the expected terminal values, and $E_f$ is the terminal inverse mechanical energy.

When an area is not allowed to be passed due to regional factors, a no-fly zone should be established, as shown in Figure 1. The center of the no-fly zone is designed on the surface of the earth at the position $N{r}_{\text{nfz}},{\theta }_{\text{nfz}},{\varphi }_{\text{nfz}}$, where $r_{\text{nfz}}$, ${\theta }_{\text{nfz}}$, and ${\varphi }_{\text{nfz}}$ are the radial distance from the Earth’s center to the center of the no-fly zone, the longitude, and the latitude of the center of the no-fly zone, respectively.

[figure omitted; refer to PDF]

Thus, in the geocentric coordinate system, coordinates of the position $N$ and the vehicle are given by $$\begin{array}{}\text{(15)}& \begin{array}{c}{x}_{\text{nfz}}=r_{\text{nfz}}\cos {\varphi }_{\text{nfz}}\cos {\theta }_{\text{nfz}},\hfill \\ y_{\text{nfz}}=r_{\text{nf}z}\cos {\varphi }_{\text{nfz}}\sin {\theta }_{\text{nfz}},\hfill \\ z_{\text{nf}z}=r_{\text{nfz}}\sin {\varphi }_{\text{nf}z},\hfill \end{array}\text{(16)}& \begin{array}{c}x=r\cos \varphi \cos \theta ,\hfill \\ y=r\cos \varphi \sin \theta ,\hfill \\ z=r\sin \varphi ,\hfill \end{array}\end{array}$$where $x_{\text{nfz}}$, $y_{\text{nfz}}$, and $z_{\text{nfz}}$ are coordinates of the vehicle and $x$, $y$, and $z$ are coordinates of the vehicle in the geocentric coordinate system.

According to relations of the geocentric coordinate system and the north-up-east (NUE) coordinate system [40], we can determine the position of the vehicle in the NUE coordinate system as follows: $$\begin{array}{}\text{(17)}& \begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}={\mathbf{M}}_1-{90}^{\circ }+{\alpha }_{\text{nfz}}{\mathbf{M}}_2{\phi }_{\text{nfz}}{\mathbf{M}}_3{\theta }_{\text{nfz}}-{90}^{\circ }\begin{array}{c}x\\ y\\ z\end{array}-\begin{array}{c}0\\ r_0\\ 0\end{array},\end{array}$$where $$\begin{array}{}\text{(18)}& {\mathbf{M}}_1-{90}^{\circ }+{\alpha }_{\text{nfz}}=\begin{array}{ccc}-\sin {\alpha }_{\text{nfz}}& 0& \cos {\alpha }_{\text{nfz}}\\ 0& 1& 0\\ -\cos {\alpha }_{\text{nfz}}& 0& -\sin {\alpha }_{\text{nfz}}\end{array},\text{(19)}& {\mathbf{M}}_2{\phi }_{\text{nfz}}=\begin{array}{ccc}1& 0& 0\\ 0& \cos {\phi }_{\text{nfz}}& \sin {\phi }_{\text{nfz}}\\ 0& -\sin {\phi }_{\text{nfz}}& \cos {\phi }_{\text{nfz}}\end{array},\text{(20)}& {\mathbf{M}}_3{\theta }_{\text{nfz}}-{90}^{\circ }=\begin{array}{ccc}\sin {\theta }_{\text{nfz}}& -\cos {\theta }_{\text{nfz}}& 0\\ \cos {\theta }_{\text{nfz}}& \sin {\theta }_{\text{nfz}}& 0\\ 0& 0& 1\end{array},\end{array}$$

$x_1$, $y_1$, and $z_1$ are coordinates of the vehicle in the NUE coordinate system, and ${\alpha }_{\text{nfz}}$ is the angle between the short half axis and the north direction.

The model of the no-fly zone constraint [41] is established as an infinite elliptical cylinder in the NUE coordinate system, which is given as $$\begin{array}{}\text{(21)}& \frac{z_1^2}{a^2}+\frac{x_1^2}{b^2}\ge 1,y_1\ge 0,\end{array}$$where $a$ and $b$ are the long and short half-axes of the ellipse, respectively.

3. Algorithm Description

Lateral parameters of the trajectory are determined by the bank angle $\sigma$. Generally, the angle of attack is determined by a designed profile. Once the amplitude of the bank angle is determined, we can obtain $L$ and $L\cos \sigma$ by (7). The amplitude of ${\psi }^{\prime }$ also can be obtained by (4), in which $L\sin \sigma$ is given by $$\begin{array}{}\text{(22)}& L\sin \sigma =\sqrt{L^2-{L\cos \sigma }^2}.\end{array}$$

Therefore, lateral parameters of the trajectory are totally determined by the sign of the bank angle. Different drag acceleration profiles can be tracked by changing the magnitude of the bank angle. The sign of the bank angle should be kept consistent so that the landing footprint problem is translated into a maximum crossrange problem for different drag acceleration profiles. The terminal altitude and velocity constraints are satisfied by adjusting the sign of the terminal bank angle. Supplementary trajectories for maximum crossrange and minimum crossrange drag acceleration profiles make the reachable area more precise. For the no-fly zone constraint, lateral parameters need to be controlled by an appropriate strategy for the sign of the bank angle.

3.1. Drag Acceleration Profile Tracking Law

In this section, we design a law to track drag acceleration profiles [42]. For $D$, we take the derivative with respect to $E$, which gives rise to $$\begin{array}{}\text{(23)}& D^{\prime }=-\frac{r_0}{h_s}-\frac{2}{V^2{R}^2}\sin \gamma -\frac{2D}{V^2}+\frac{D}{C_D}{C}_D^{\prime }.\end{array}$$

Taking the derivative of $D^{\prime }$ with respect to $E$ yields $$\begin{array}{}\text{(24)}& D^{″}=\frac{4}{R^2{V}^2}\frac{V^{\prime }}{V}+\frac{R^{\prime }}{R}\sin \gamma +-\frac{r_0}{h_s}-\frac{2}{R^2{V}^2}{\gamma }^{\prime }\cos \gamma +\frac{4D}{V^3}{V}^{\prime }-\frac{2}{V^2}{D}^{\prime }+\frac{C_D^{\prime }}{C_D}{D}^{\prime }+D\frac{C_D^{″}}{C_D}-\frac{{C^{\prime }}_D^2}{C_D^2}.\end{array}$$

By substituting (5) into (24), we obtain $L\cos \sigma =D^{″}-b/a$, which gives rise to $$\begin{array}{}\text{(25)}& \sigma =\arccos \frac{D^{″}-a}{bL},\end{array}$$where $$\begin{array}{}\text{(26)}& a=D\frac{C_D^{″}}{C_D^{\prime }}-\frac{{C^{\prime }}_D^2}{C_D^2}+D^{\prime }\frac{C_D^{\prime }}{C_D}+\frac{2}{V^2}-\frac{4D}{V^4}+\frac{1}{{DV}^2}\frac{1}{h_s}+\frac{2g}{V^2}g-\frac{V^2}{r},\text{(27)}& b=-\frac{1}{V^4}\frac{1}{h_s}+\frac{2g}{V^2}.\end{array}$$

A linear feedback control law is used to track drag acceleration profiles, as follows: $$\begin{array}{}\text{(28)}& D_{ref}^{″}-D^{″}+2{\xi \omega }_n{D}_{ref}^{\prime }-D^{\prime }+{\omega }_n^2{D}_{ref}-D=0,\end{array}$$where $\xi$ is the damp coefficient, ${\omega }_n$ is the natural frequency of the tracking law, $D$ is the real drag acceleration, and $D_{\text{ref}}$ is the reference drag acceleration.

We can obtain $D^{″}$ by (28) and substitute it into (25) to obtain the drag acceleration profile tracking law denoted as $$\begin{array}{}\text{(29)}& \sigma =\arccos D_{\text{ref}}^{″}+2{\xi \omega }_n{D}_{\text{ref}}^{\prime }-D^{\prime }+{\omega }_n^2{D}_{\text{ref}}-D-a/bL.\end{array}$$

In this paper, simulations are carried out on the Common Aero Vehicle (CAV) [43]. Let the input reference drag acceleration $D_{\text{ref}}=13$ m/s², the initial altitude $h_0=35$ km, the attack angle $\alpha ={20}^{\circ }$, and the constant velocity $V=3000$ m/s. Drag acceleration command tracking processes are shown in Figures 2 and 3. The tracking system can depress the overshoot but slow down the tracking speed for a higher damp coefficient $\xi$ and a constant natural frequency ${\omega }_n$. If the damp coefficient is constant, the higher the natural frequency, the faster the tracking speed. However, if the natural frequency is too large, the high overshoot may be caused. Thus, tracking law coefficients are designed $\xi =2$ and ${\omega }_n=5e-6$. If there is a white noise with a mean of 0 and a standard deviation of 0.5 in the drag acceleration command, the tracking process is shown in Figure 4, and the tracking system can filter the noise. The overall tracking speed is slow because the initial state of the tracking system is far away from the reference drag acceleration command, so the vehicle will take some time to fly to the appropriate altitude to track the reference command.

[figure omitted; refer to PDF]

To achieve a smooth maximum downrange trajectory, we use the equilibrium glide condition by setting $\gamma \le \epsilon$, where $\epsilon$ is a small vector. The bank angle tracking law is designed to obtain the relationship between the drag acceleration $D$ and the inverse energy $E$ at the lowest boundary of the drag profiles through numerical integration of the equations of motion. We use the feedback linearization to keep the flight path angle close to zero during the flight. To track the reference flight path angle ${\gamma }_{\text{ref}}$, the feedback control law can be written as $$\begin{array}{}\text{(31)}& {\gamma }_{\text{ref}}^{\prime }-{\gamma }^{\prime }+K_{\gamma }{\gamma }_{\text{ref}}-\gamma =0,\end{array}$$where $K_{\gamma }$ is a tracking coefficient.

Substituting (5) into (31), we have $$\begin{array}{}\text{(32)}& \frac{L}{D}\cos \sigma =\frac{K_{\gamma }V}{D}{\gamma }_{\text{ref}}-\gamma +\frac{1}{D}g-\frac{V^2}{r}\cos \gamma .\end{array}$$

To reduce bounces of the trajectory, we set the reference path angle ${\gamma }_{\text{ref}}=‐0{.01}^{\circ }$. The reference path angle can be tracked in the initial phase, but the phugoid motion cannot be avoided during a whole flight with a fixed bank angle. The magnitude of the bank angle has been constrained to be ${90}^{\circ }$ at an appropriate switching time to satisfy the terminal velocity constraint.

A simulation is carried out with parameters in section 4 to study the effect of different values for the tracking coefficient $K_{\gamma }$ on the equilibrium glide condition, as shown in Figure 6. The excessive coefficient may cause the divergence of the tracking system. According to Figure 6, we set $K_{\gamma }=1$ in this paper.

[figure omitted; refer to PDF]

For profiles between the maximum downrange and minimum downrange drag acceleration profiles, the interpolation equation is given by $$\begin{array}{}\text{(34)}& D=D_{\min }+\frac{i^2}{n^2}{D}_{\max }-D_{\min },\end{array}$$where $D_{\max }$ and $D_{\min }$ are the maximum and minimum drag acceleration profiles corresponding to the minimum and maximum range, respectively. $i$ presents the $i$th drag acceleration profile, and $n$ is the number of drag acceleration profiles.

3.4. Supplementary Trajectories

Since the bank angle is kept positive or negative during flight, the interval of landing footprints between the maximum crossrange trajectories is considerable. For maximum downrange trajectories, we apply an interpolation between the heading angles ${\psi }_{\max 1}={\psi }_0$ and ${\psi }_{\max 2}={\psi }_{\text{egcf}}$, where ${\psi }_0$ is the initial entry heading angle and ${\psi }_{\text{egcf}}$ is the terminal heading angle of the maximum downrange trajectory; thus, the supplementary terminal constraint ${\psi }_{\max }$ is given by $$\begin{array}{}\text{(35)}& {\psi }_{\max }={\psi }_{\max 1}+\frac{{i_{\psi }}^2}{{n_{\psi }}^2}{\psi }_{\max 1}-{\psi }_{\max 2},\end{array}$$where $i_{\psi }$ represents the $i$th supplementary trajectory and $n_{\psi }$ is the number of supplementary trajectories. Moreover, another supplementary terminal constraint is given by ${\psi }_{\min }={\psi }_{\max }-\delta \psi$, and $\delta \psi$ is the threshold of the heading angle.

Similarly, for minimum downrange trajectories, we apply an interpolation between the heading angles ${\psi }_{\min 1}={\psi }_0$ and ${\psi }_{\min 2}={\psi }_{\text{consf}}$, where ${\psi }_{\text{consf}}$ is the terminal heading angle of the minimum downrange trajectory. Thus, the supplementary terminal constraint ${\psi }_{\min }$ is given by $$\begin{array}{}\text{(36)}& {\psi }_{\min }={\psi }_{\min 1}+\frac{{i_{\psi }}^2}{{n_{\psi }}^2}{\psi }_{\min 1}-{\psi }_{\min 2}.\end{array}$$

Moreover, another supplementary terminal constraint is given by ${\psi }_{\max }={\psi }_{\min }+\delta \psi$.

For the supplementary terminal constraints, supplementary trajectories are designed by adjusting the sign of the bank angle, as given by $$\begin{array}{}\text{(37)}& \mathrm{sign}{\sigma }^i=\begin{array}{cc}\mathrm{sgn}{\psi }_{\max }-\psi ,\hfill & \psi >{\psi }_{\max }\hfill \\ \mathrm{sgn}{\psi }_{\min }-\psi ,\hfill & \psi <{\psi }_{\min }\hfill \\ \mathrm{sgn}{\sigma }^{i-1},\hfill & {\psi }_{\min }\le \psi \le {\psi }_{\max }\hfill \end{array},\end{array}$$where ${\sigma }^{i-1}$ is the bank angle at the previous moment and ${\sigma }^{i-1}$ is the bank angle at this moment. $\delta \psi$ is the threshold of the heading angle. If the value of $\delta \psi$ is too large, the result may be inaccurate. On the contrary, if the value is too small, it will cause the rapid reversal of the bank angle, which may cause the out of control.

3.5. Flying-Around Strategy for the No-Fly Zone Constraint

For the no-fly zone constraint, according to the relationship among the no-fly zone, the position of the vehicle, and the heading angle of the vehicle, we designed a strategy to prevent the vehicle from entering the no-fly zone. In Section 2.2, we established a model of a no-fly zone. In this paper, the elliptical no-fly zone is simplified as a circular zone by the assignment of $a=R_z$ and $b=R_z$ in (21), in which $R_z$ is the radius of a circular no-fly zone. Assuming that the center of the no-fly zone $N{\theta }_{\text{nfz}},{\varphi }_{\text{nfz}}$ and the position of the vehicle $M\theta ,\varphi$ are in the same horizontal plane, the relationship between $N$ and $M$ is shown in Figure 8.

[figure omitted; refer to PDF]

$MA$ is the tangent from $M$ to the no-fly zone, and $A$ is the tangent point. ${\overline{\psi }}_A$ is the angle between the tangent $MA$ and the north, and ${\overline{\psi }}_Z$ is the angle between $MN$ and the north. $\overline{\psi }$ is the redefined heading angle, where $\overline{\psi }=‐\psi$. ${\overline{\psi }}_Z$ is solved by the following equation: $$\begin{array}{}\text{(38)}& {\overline{\psi }}_Z=\begin{array}{c}arc\tan \frac{\sin {\theta }_{\text{nfz}}-\theta }{\cos \varphi \tan {\varphi }_{\text{nfz}}-\sin \varphi \cos {\theta }_{\text{nfz}}-\theta },\varphi <{\varphi }_{n\text{fz}},\hfill \\ {\overline{\psi }}_Z=\pi -arc\tan \frac{\sin {\theta }_{\text{nfz}}-\theta }{\cos \varphi \tan {\varphi }_{\text{nfz}}-\sin \varphi \cos {\theta }_{\text{nfz}}-\theta },\varphi \ge {\varphi }_{\text{nfz}}.\hfill \end{array}\end{array}$$

If the vehicle passes the no-fly zone from the right side, the boundary of the heading angle is given by $$\begin{array}{}\text{(39)}& \begin{array}{c}{\overline{\psi }}_{\min }={\overline{\psi }}_A,\hfill \\ {\overline{\psi }}_{\max }={\overline{\psi }}_A+\delta \overline{\psi },\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(40)}& {\overline{\psi }}_A=arc\tan \frac{\sin {\theta }_{\text{nfz}}-\theta }{\cos \varphi \tan {\varphi }_{\text{nfz}}-\sin \varphi \cos {\theta }_{\text{nfz}}-\theta }+\arcsin \frac{R_Z}{MN},\end{array}$$

$\delta \overline{\psi }$ is the threshold of the redefined heading angle, and $MN$ is the distance between the point $M$ and point $N$. The selection of the value for $\delta \overline{\psi }$ is the same as $\delta \psi$ in Section 3.4.

If the vehicle passes the no-fly zone from the left side, the boundary of the heading angle is given by $$\begin{array}{}\text{(41)}& \begin{array}{c}{\overline{\psi }}_{\min }={\overline{\psi }}_A-\delta \overline{\psi },\hfill \\ {\overline{\psi }}_{\max }={\overline{\psi }}_A,\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(42)}& {\overline{\psi }}_A=arc\tan \frac{\sin {\theta }_{\text{nfz}}-\theta }{\cos \varphi \tan {\varphi }_{\text{nfz}}-\sin \varphi \cos {\theta }_{\text{nfz}}-\theta }‐\arcsin \frac{R_Z}{MN}.\end{array}$$

The sign of the bank angle is determined similarly to the equation in Section 3.4.

3.6. Landing Footprint Solution

In this section, the reachable method is designed based on the condition that the no-fly zone is located on the right side of the initial heading direction of the vehicle. Under the condition of a left no-fly zone, a similar strategy can be designed but is not described in detail in this paper.

According to the relationship between the no-fly zone and landing footprints, the solution is divided into five cases, as follows: Case 1, without no-fly zones; Case 2, the no-fly zone is located entirely outside the area of landing footprints; Case 3, the no-fly zone is located entirely inside the area of landing footprints; Case 4, the no-fly zone intersects the edges of the supplementary minimum downrange landing footprints, while the center of the no-fly zone is located outside the area of landing footprints; and Case 5, the no-fly zone intersects the edges of the supplementary minimum downrange landing footprints, while the center of the no-fly zone is located inside the area of landing footprints. For the cases in which the no-fly zone intersects the edges of the supplementary maximum downrange or interpolated maximum crossrange landing footprints, the no-fly zone has a small effect on landing footprints, which is not described in this paper. The set of landing footprints is denoted as $J$, and the trajectory is denoted as $j$. The steps based on EAGLE for solving the landing footprints for each case are described as follows:

Case 1.

(a) Initial landing footprints: interpolated drag acceleration profiles between the maximum downrange and minimum downrange drag acceleration profiles determined by (34) are tracked by the drag acceleration profile tracking law presented in Section 3.1. By keeping the bank angle positive or negative during the whole flight, we can obtain sets of interpolated landing footprints $J_1$ and $J_4$. The trajectory is denoted as $j_2$ or $j_4$ when it corresponds to the positive or negative bank angle tracking from the maximum downrange drag acceleration profile. Similarly, the trajectory is denoted as $j_1$ or $j_5$ when it corresponds to the positive or negative bank angle tracking from the minimum downrange drag acceleration profile. The trajectory with the terminal heading angle ${\psi }_0$ corresponding to the maximum downrange drag acceleration profile is denoted as $j_3$ and that corresponding to the minimum downrange drag acceleration profile is denoted as $j_6$

(b) Supplementary landing footprints: we can obtain the set of interpolated landing footprints $J_2$ or $J_3$ constrained to the supplementary terminal heading angle by (35), in which the value of ${\psi }_{\text{egcf}}$ is ${\psi }_2$ (the terminal heading angle of $j_2$) or ${\psi }_4$ (the terminal heading angle of $j_4$). The sign of the bank angle is given by (37). Similarly, we can obtain the set of interpolated landing footprints $J_5$ or $J_6$ in the same way, corresponding to the minimum downrange drag acceleration profile

(c) Final landing footprints: we can obtain final landing footprints by connecting sets of landing footprints $J_1-J_2-J_3-J_4-J_5-J_6$ successively. The sets of landing footprints correspond to Figure 9

Case 2.

(a) Initial landing footprints: sets of landing footprints $J_2$, $J_3$, $J_4$, and $J_5$ are obtained by steps a and b in Case 1. For trajectories not affected by the no-fly zone in $J_1$, we denote the trajectory closest to the no-fly zone from the left as $j_7$ and closest from the right as $j_{11}$. Landing footprints calculated from $j_7$ to $j_2$ are denoted as $J_7$, and landing footprints from $j_{11}$ and $j_1$ are denoted as $J_{16}$

(b) Modified landing footprints for the left side: use the left flying-around strategy in Section 3.5 for trajectories between $j_7$ and $j_1$ and denote these landing footprints as $J_8$. The terminal heading angles of the flying-around trajectories $j_8$ and $j_9$ are denoted as ${\psi }_8$ and ${\psi }_9$, corresponding to the original trajectories $j_7$ and $j_1$. Where the trajectory $j_9$ and the no-fly zone intersect at the point of tangency, we denote the heading angle here as ${\psi }_{10}$. The interpolation between ${\psi }_9$ and ${\psi }_{10}$ is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints $J_9$ by the same method used in Section 3.4. Similarly, we can obtain set $J_{10}$ by the interpolation of ${\psi }_0$ and ${\psi }_{10}$

(c) Modified landing footprints for the right side: landing footprints calculated from $j_{11}$ to $j_1$ are denoted as $J_{16}$. Use the right flying-around strategy in Section 3.5 for the trajectories between $j_{11}$ and $j_1$ and denote these landing footprints as $J_{13}$. The terminal heading angle of flying-around trajectory $j_{13}$ or $j_{14}$ is denoted as ${\psi }_{13}$ or ${\psi }_{14}$, corresponding to the original trajectory $j_{11}$ or $j_1$. The terminal heading angle of the trajectory $j_{11}$ or $j_1$ is denoted as ${\psi }_{11}$ or ${\psi }_1$. Where the trajectory $j_{13}$ or $j_{14}$ and the no-fly zone intersect at a point of tangency, denote the heading angle here as ${\psi }_{12}$ or ${\psi }_{15}$. The interpolation between ${\psi }_{11}$ and ${\psi }_{12}$ is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints $J_{11}$ by the same method used in Section 3.4. Similarly, we can obtain set $J_{12}$, $J_{14}$, or $J_{15}$ by the interpolation of ${\psi }_{12}$ and ${\psi }_{13}$, ${\psi }_{14}$, and ${\psi }_{15}$, or ${\psi }_{15}$ and ${\psi }_1$

(d) Final landing footprint: the area of landing footprints consists of two closed figures constructed by connecting sets $J_2-J_3-J_4-J_5-J_{10}-J_9-J_8-J_7$ and $J_{16}-J_{11}-J_{12}-J_{13}-J_{14}-J_{15}$ successively. The sets of landing footprints correspond to Figure 10

Case 3.

(a) Initial landing footprints: sets of landing footprints $J_2$, $J_3$, $J_4$, $J_5$, $J_6$, $J_7$, and $J_{16}$ are obtained through steps a and b in Case 2

(b) Modified landing footprints for the left side: we can obtain set $J_8$ through a process which is similar to step b in Case 2 but slightly different. If a trajectory will not reach the no-fly zone at the point of tangency, the landing footprint of this trajectory will no longer belong to $J_8$. In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in $J_8$ is denoted as $j_9$

(c) Modified landing footprints for the right side: the trajectories $j_{11}$, $j_{12}$, $j_{13}$, set $J_{11}$, and set $J_{12}$ are given by the same step in Case 2. We can obtain set $J_{13}$ by a process that is similar to step c in Case 2 but slightly different. If a trajectory will not reach the no-fly zone at the point of tangency, the landing footprint of this trajectory will no longer belong to $J_{13}$. In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in $J_{13}$ is denoted as $j_{14}$

(d) Part of the no-fly zone: select points $p_1$ and $p_2$ closest to the landing footprints of $j_9$ and $j_{14}$ on the edge of the no-fly zone and denote points between $p_1$ and $p_2$ as $P_1$

(e) Final landing footprint: we can obtain the final landing footprints by connecting sets of landing footprints $J_2-J_3-J_4-J_5-J_6-J_{16}-J_{11}-J_{12}-J_{13}-P_1-J_8-J_7$ successively. The sets of landing footprints correspond to Figure 11

Case 4.

The steps used to solve the landing footprints for Case 4 are the same as the steps used for Case 2. The sets of landing footprints correspond to Figure 12.

Case 5.

(a)–(c) The steps are the same as steps a to c for Case 3.

(d) Modified landing footprints: the trajectory corresponding to the landing footprint closest to the no-fly zone from the left side is denoted as $j_{17}$ and that closest from the right side is denoted as $j_{18}$. The landing footprints of the trajectories between $j_6$ and $j_{17}$ are denoted as $J_{17}$, and the landing footprints of the trajectories between $j_{18}$ and $j_1$ are denoted as $J_{18}$

(e) Part of the no-fly zone: select points $p_1$ and $p_3$ closest to the landing footprints of $j_9$ and $j_{17}$ on the edge of the no-fly zone and denote the points between $p_1$ and $p_3$ as $P_1$. Select points $p_2$ and $p_4$ closest to the landing footprints of $j_{14}$ and $j_{18}$ on the edge of the no-fly zone and denote the points between $p_2$ and $p_4$ as $P_2$

(f) Final landing footprint: the area of the landing footprints consists of two closed figures constructed by connecting sets $J_2-J_3-J_4-J_5-J_{17}-P_1-J_8-J_7$ and $J_{16}-J_{11}-J_{12}-J_{13}-P_2-J_{18}$ successively

The above solution is based on the assumption that the no-fly zone affects sets of landing footprints $J_1$, $J_4$, $J_5$, and $J_6$. If the no-fly zone has an impact on $J_2$ or $J_3$, the key algorithm is the same as that in the above study and is not described in this paper. The set of landing footprints corresponds to Figure 13.

[figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF]

4. Numerical Examples

In this section, a 3-degrees model of the CAV is used for the landing footprint simulation. To verify the algorithm, MATLAB is used to model the problem. The mass of the vehicle is 907.2 kg, and the aerodynamic reference area is 0.484 m². The path constraints are set as ${\dot{Q}}_{\max }=2000$ kW/m², $q_{\max }=100$ kPa, $n_{\max }=3$ g, and $k_Q=1\text{e}‐4$. The initial states of the vehicle are shown in Table 1. The terminal constraints for the altitude and the velocity are 3000 m/s and 30 km, respectively. And thresholds of the terminal constraints are designed as $\Delta h=100$ m and $\Delta V=50$ m/s. The angles of attack profile parameters are designed as ${\alpha }_{\max }={20}^{\circ }$, ${\alpha }_{\max L/D}={20}^{\circ }$, ${Ma}_{\alpha 1}=12$, and ${Ma}_{\alpha 2}=12$. The tracking law coefficients for the drag acceleration command tracking process are designed $\xi =1$ and ${\omega }_n=5e-6$. Similarly, the tracking law coefficient for the equilibrium glide condition is designed $K_{\gamma }=1$. The threshold of the heading angle in Section 3.4 and Section 3.5 is designed as $\delta \psi =2^{\circ }$. The number of drag acceleration profiles is designed as $n=10$, and the number of supplementary trajectories is designed as $n_{\psi }=10$. Thus, for Case 1, the number of total trajectories is 80. For comparison, the number of total trajectories for the Gauss pseudospectral method is also set as 80.

Table 1

Initial states of the entry vehicle.

$\theta {}^{\circ }$	$\varphi {}^{\circ }$	$h_{0}km$	$\psi {}^{\circ }$	$\gamma {}^{\circ }$	$Vm/s$
0	0	60	0	-1	5000

4.1. Without No-Fly Zones

In this section, the no-fly zone constraint is not added, corresponding to Case 1, and the effectiveness of the method proposed in this paper based on EAGLE is verified. Figures 14 and 15 show the tracking ability of the minimum and maximum acceleration profiles. At the initial point and the turning point of the attack angle profile, there are little lags in the tracking process, but the overall tracking is still good.

[figure omitted; refer to PDF]

Table 2 shows the computation times of different methods. The solving time for the method proposed in this paper is 13.42 s, which is suitable for the online rapid engineering, in contrast to the 429.93 s by the Gauss pseudospectral method, which is suitable for accurate offline engineering.

Table 2

Computation times of different methods.

Computation time of the EAGLE method	Computation time of the Gauss pseudospectral method
Equilibrium glide condition: 0.14 sPath constraint translation: 0.23Each tracking trajectory: 0.16 s	–
Total time: 13.42 s	Total time: 429.93 s

Figure 16 compares landing footprints obtained by the EAGLE and Gauss pseudospectral methods, and we can see that landing footprints obtained by EAGLE are significantly smaller. This is because the trajectory based on the drag acceleration profile has not been optimized, the calculation process of landing footprints obtained by EAGLE has some approximate simplification steps, and the equilibrium glide condition is added. Figure 17 shows altitude profiles from the equilibrium glide condition to transformed path constraints. Due to the flight capacity and terminal constraints, the maximum range trajectory maintains the equilibrium glide condition only for a while, and the bouncing of the trajectory cannot be avoided during the whole process. Figures 18 and 19 show drag acceleration profiles and corresponding lift acceleration profiles. Drag acceleration profiles are obtained by tracking the designed interpolation profiles. Figures 20 and 21 show bank angle profiles and the angle of attack profiles. And Figures 22 and 23 show flight path angle profiles and heading angle profiles.

[figure omitted; refer to PDF]

Figure 24 shows the dispersion of the terminal altitude within 12 m, and Figure 25 shows the dispersion of the terminal velocity within 10 m/s, which satisfy terminal constraints. Figures 26–28 show all path constraint profiles, which also satisfy proposed path constraints.

[figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF]

4.2. With No-Fly Zone

In this section, the no-fly zone constraint is added. Centers of designed no-fly zones are set as $1{.4}^{\circ },10^{\circ }$, $3^{\circ },12{.5}^{\circ }$, $2^{\circ },10{.6}^{\circ }$, and $2^{\circ },11{.1}^{\circ }$, corresponding to Cases 2-5 in Section 3.6, respectively. The radius of no-fly zones is designed as 40 km, and the angle between the short half axis and the north direction is designed as ${\alpha }_{\text{nfz}}=0^{\circ }$. The number of the interpolated trajectories is designed as 5, which are supplementary trajectories for the flying-around strategy in Section 3.6. By applying the flying-around strategy described in Section 3.5 and the landing footprint solution presented in Section 3.6, we can obtain final modified landing footprints shown in Figures 29–32. To see the effect more clearly, landing footprints near the no-fly zone are partially magnified. Except for the Case 3, final landing footprints consist of two closed figures. Table 3 shows the computation times of different cases. Trajectories affected by the no-fly zone mean flying-around trajectories and supplementary trajectories for the flying-around strategy. The solving time for these cases is short, which is suitable for the online rapid engineering.

[figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF]

Table 3

Computation times of different cases.

Different cases	Number of trajectories not affected by the no-fly zone	Computation time of trajectories not affected by the no-fly zone	Number of trajectories affected by the no-fly zone	Computation time of trajectories affected by the no-fly zone	Total computation time
Case 2	69	11.21 s	46	7.96 s	19.18 s
Case 3	77	12.51 s	21	3.64 s	16.15 s
Case 4	73	11.86 s	33	5.71 s	17.57 s
Case 5	73	11.86 s	30	5.19 s	17.06 s

5. Conclusions

In this paper, a fast method to generate the reachable area for entry vehicles is developed, which is subjected to inequality entry path constraints, terminal constraints, and the no-fly zone constraint. The reachable area algorithm follows the basic reference profile tracking scheme and consists of three modules: drag acceleration-energy profiles and the corresponding tracking law are designed to achieve initial landing footprints; a flying-around strategy is designed to avoid the no-fly zone; the solution of the reachable area is proposed by the reasonable distribution of landing footprints. Compared with the Gauss pseudospectral method, the rapidity and the effectiveness of the proposed solution are comprehensively analyzed. The results show that the solution has the good applicability and dynamic characteristics. Therefore, the online reachable area solution developed in this paper provides the key information for the trajectory planning of hypersonic vehicles. In addition, the method has some shortcomings in terms of accuracy. Future studies will try to find a highly accurate online landing footprint solution for no-fly zones.

References

[1] L. Mu, X. Yu, Y. M. Zhang, P. Li, X. Wang, "Onboard guidance system design for reusable launch vehicles in the terminal area energy management phase," Acta Astronautica, vol. 143 no. 1, pp. 62-75, DOI: 10.1016/j.actaastro.2017.10.027, 2018.

[2] S. R. Ploen, H. Seraji, C. E. Kinney, "Determination of spacecraft landing footprint for safe planetary landing," IEEE Transactions on Aerospace and Electronic Systems, vol. 45 no. 1,DOI: 10.1109/TAES.2009.4805259, 2009.

[3] Z. Wang, "Optimal trajectories and normal load analysis of hypersonic glide vehicles via convex optimization," Aerospace Science and Technology, vol. 87 no. 4, pp. 357-368, DOI: 10.1016/j.ast.2019.03.002, 2019.

[4] Z. Wang, M. J. Grant, "Constrained trajectory optimization for planetary entry via sequential convex programming," Journal of Guidance Control and Dynamics, vol. 40 no. 10, pp. 2603-2615, DOI: 10.2514/1.g002150, 2017.

[5] Z. Wang, M. J. Grant, "Autonomous entry guidance for hypersonic vehicles by convex optimization," Journal of Spacecraft and Rockets, vol. 55 no. 4, pp. 993-1006, DOI: 10.2514/1.a34102, 2018.

[6] R. Chai, A. Savvaris, A. Tsourdos, S. Chai, Y. Xia, "Stochastic spacecraft trajectory optimization with the consideration of chance constraints," IEEE Transactions on Control Systems Technology, vol. 28 no. 4, pp. 1550-1559, DOI: 10.1109/TCST.2019.2908938, 2019.

[7] R. Q. Chai, A. Tsourdos, A. Savvaris, S. C. Chai, "Trajectory planning for hypersonic reentry vehicle satisfying deterministic and probabilistic constraints," Acta Astronautica, vol. 77 no. 12, 2020.

[8] R. Q. Chai, A. Tsourdos, A. Savvaris, S. Chai, Y. Xia, C. L. P. Chen, "Six-dof spacecraft optimal trajectory planning and real-time attitude control: a deep neural network-based approach," IEEE Transactions on Neural Networks and Learning Systems, vol. 31 no. 11, pp. 5005-5013, DOI: 10.1109/TNNLS.2019.2955400, 2020.

[9] L. Cheng, F. Jiang, Z. Wang, J. Li, "Multi-constrained real-time entry guidance using deep neural networks," IEEE Transactions on Aerospace and Electronic Systems, vol. 9,DOI: 10.1109/TAES.2020.3015321, 2020.

[10] R. Chai, A. Tsourdos, A. Savvaris, Y. Xia, S. Chai, "Real-time reentry trajectory planning of hypersonic vehicles: a two-step strategy incorporating fuzzy multi-objective transcription and deep neural network," IEEE Transactions on Industrial Electronics, vol. 67 no. 8, pp. 6904-6915, DOI: 10.1109/tie.2019.2939934, 2020.

[11] Y. Shi, Z. B. Wang, "A deep learning-based approach to real-time trajectory optimization for hypersonic vehicles," AIAA Scitech 2020 Forum, .

[12] P. Lu, C. W. Brunner, S. J. Stachowiak, G. F. Mendeck, M. A. Tigges, C. J. Cerimele, "Verification of a fully numerical entry guidance algorithm," Journal of Guidance and Control Dynamics, vol. 40 no. 2, pp. 230-247, DOI: 10.2514/1.g000327, 2017.

[13] L. Cheng, Z. Wang, Y. Cheng, Q. Zhang, K. Ni, "Multi-constrained predictor-corrector reentry guidance for hypersonic vehicles," Proceedings of the Institution of Mechanical Engineers Part G Journal of Aerospace Engineering, vol. 232 no. 16, pp. 3049-3067, DOI: 10.1177/0954410017724185, 2018.

[14] R. Q. Chai, A. Tsourdos, A. Savvaris, S. C. Chai, "Trajectory optimization of space maneuver vehicle using a hybrid optimal control solver," IEEE Transactions on Cybernetics, vol. 49 no. 2, 2017.

[15] E. Besada-Portas, L. de la Torre, J. M. de la Cruz, B. de Andrés-Toro, "Evolutionary trajectory planner for multiple UAVs in realistic scenarios," IEEE Transactions on Robotics, vol. 26 no. 4, pp. 619-634, DOI: 10.1109/tro.2010.2048610, 2010.

[16] N. Zhang, C. Wang, F. Li, "A rapid landing footprint calculation method for entry vehicles," Proceedings of IEEE Chinese Guidance, Navigation and Control Conference, pp. 553-557, .

[17] Z. Liang, S. Liu, Q. Li, Z. Ren, "Lateral entry guidance with no-fly zone constraint," Aerospace Science and Technology, vol. 60 no. 12, pp. 39-47, DOI: 10.1016/j.ast.2016.10.025, 2017.

[18] J. A. Englander, B. A. Conway, "An automated solution of the low-thrust interplanetary trajectory problem," Journal of Guidance Control and Dynamics, vol. 40 no. 1, pp. 15-27, DOI: 10.2514/1.g002124, 2017.

[19] R. Z. He, L. H. Liu, G. J. Tang, W. M. Bao, "Rapid generation of entry trajectory with multiple no-fly zone constraints," Advances in Space Research, vol. 60 no. 7, pp. 1430-1442, DOI: 10.1016/j.asr.2017.06.046, 2017.

[20] Y. Liang, Y. Jin, C. Wanchun, "Entry guidance with no-fly zone avoidance using linear pseudospectral model predictive control," IEEE Access, vol. 7, pp. 119246-119258, 2019.

[21] Y. Gao, G. Cai, X. Yang, "Improved tentacle-based guidance for reentry gliding hypersonic vehicle with no-fly zone constraint," IEEE Access, vol. 7, pp. 98589-98602, 2019.

[22] W. B. Yu, W. C. Chen, Z. G. Jiang, "Analytical entry guidance for coordinated flight with multiple no-fly-zone constraints," Aerospace Science and Technology, vol. 84 no. 10, 2018.

[23] W. Yu, W. Chen, Z. Jiang, M. Liu, D. Yang, M. Yang, Y. Ge, "Analytical entry guidance for no-fly-zone avoidance," Aerospace Science and Technology, vol. 72 no. 1, pp. 426-442, DOI: 10.1016/j.ast.2017.11.029, 2018.

[24] R. Chai, A. Savvaris, A. Tsourdos, "Violation learning differential evolution-based hp-adaptive pseudospectral method for trajectory optimization of space maneuver vehicle," IEEE Transactions on Aerospace and Electronic Systems, vol. 53 no. 4, pp. 2031-2044, DOI: 10.1109/TAES.2017.2680698, 2017.

[25] B. Wu, D. Wang, E. K. Poh, G. Xu, "Nonlinear optimization of low-thrust trajectory for satellite formation: Legendre pseudospectral approach," Journal of Guidance Control and Dynamics, vol. 32 no. 4, pp. 1371-1381, DOI: 10.2514/1.37675, 2009.

[26] F. Fahroo, D. Doman, A. Ngo, "Footprint generation for reusable launch vehicles using a direct pseudospectral method," Proceedings of American Control Conference, pp. 2163-2168, .

[27] F. Fahroo, D. Doman, "A direct method for approach and landing trajectory reshaping with failure effect estimation," Proceedings of AIAA Guidance, Navigation, and Control Conference and Exhibit, .

[28] Y. Liu, J. Heiligers, M. Ceriotti, "Loosely-displaced geostationary orbits with hybrid sail propulsion," Aerospace Science and Technology, vol. 79 no. 9, pp. 105-117, DOI: 10.1016/j.ast.2018.05.034, 2018.

[29] F. Han, Z. Wang, L. He, "Trajectory plan for an ultra-short distance on-orbit service based on the Gaussian pseudo-spectral method," IEEE/CAA Journal of Automatica Sinica, vol. 79 no. 2, 2018.

[30] B. Tian, W. Fan, R. Su, Q. Zong, "Real-time trajectory and attitude coordination control for reusable launch vehicle in reentry phase," IEEE Transactions on Industrial Electronics, vol. 62 no. 3, pp. 1639-1650, DOI: 10.1109/TIE.2014.2341553, 2015.

[31] P. Lu, S. Xue, "Rapid generation of accurate entry landing footprints," Journal of Guidance Control and Dynamics, vol. 33 no. 3, pp. 756-767, DOI: 10.2514/1.46833, 2010.

[32] P. Lu, "Asymptotic analysis of quasi-equilibrium glide in lifting entry flight," Journal of Guidance Control and Dynamics, vol. 29 no. 3, pp. 662-670, DOI: 10.2514/1.15789, 2006.

[33] X. Liu, Z. Shen, P. Lu, "Solving the maximum-crossrange problem via successive second-order cone programming with a line search," Aerospace Science and Technology, vol. 47 no. 12, pp. 10-20, DOI: 10.1016/j.ast.2015.09.008, 2015.

[34] H. F. Li, R. Zhang, Z. Y. Li, "Footprint problem with angle of attack optimization for high lifting reentry vehicle," Chinese Journal of Aeronautics, vol. 25 no. 2, pp. 105-113, 2012.

[35] A. D. Ngo, D. B. Doman, "Footprint determination for reusable launch vehicles experiencing control effector failures," Proceedings of AIAA Guidance, Navigation and Control Conference and Exhibit, .

[36] A. Saraf, J. A. Leavitt, M. Ferch, "Landing footprint computation for entry vehicles," Proceedings of AIAA Guidance, Navigation and Control Conference and Exhibit, .

[37] R. Chai, A. Savvaris, A. Tsourdos, Y. Xia, S. Chai, "Solving multi-objective constrained trajectory optimization problem by an extended evolutionary algorithm," IEEE Transactions on Cybernetics, vol. 50 no. 4, pp. 1630-1643, DOI: 10.1109/tcyb.2018.2881190, 2020.

[38] M. SagliaNo, E. Mooij, S. Theil, "Adaptive disturbance-based high-order sliding-mode control for hypersonic-entry vehicles," Journal of Guidance Control and Dynamics, vol. 40 no. 3, pp. 521-536, DOI: 10.2514/1.G000675, 2017.

[39] J. Ji, A. Khajepour, W. W. Melek, Y. Huang, "Path planning and tracking for vehicle collision avoidance based on model predictive control with multiconstraints," IEEE Transactions on Vehicular Technology, vol. 66 no. 2, pp. 952-964, DOI: 10.1109/TVT.2016.2555853, 2017.

[40] Y. Xie, L. Liu, G. Tang, W. Zheng, "Highly constrained entry trajectory generation," Acta Astronautica, vol. 88 no. 1, pp. 44-60, DOI: 10.1016/j.actaastro.2013.01.024, 2013.

[41] W. Yu, W. Chen, "Guidance law with circular no-fly zone constraint," Nonlinear Dynamics, vol. 78 no. 3, pp. 1953-1971, DOI: 10.1007/s11071-014-1571-2, 2014.

[42] A. Saraf, J. A. Leavitt, D. T. Chen, K. D. Mease, "Design and evaluation of an acceleration guidance algorithm for entry," Journal of Spacecraft and Rockets, vol. 41 no. 6, pp. 986-996, DOI: 10.2514/1.11015, 2004.

[43] Y. Wang, W. Zhou, "An improved reentry trajectory planning method for Common Aero Vehicle," Proceedings of IEEE Chinese Guidance, Navigation and Control Conference, pp. 369-374, .

[44] R. G. Gottlieb, W. T. Fowler, "Improved secant method applied to boost trajectory optimization," Journal of Spacecraft and Rockets, vol. 14 no. 2, pp. 96-101, DOI: 10.2514/3.57166, 1977.