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

Authors: Shengnan Fu; Tianyu Lu; Jian Yin; Qunli Xia
Journal: International Journal of Aerospace Engineering
 (2021)

---

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


---

## Body

## 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 Section2. 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 energyE is an appropriate independent variable, which is determined by E=μ/r−V2/2, where r is the radial distance from the Earth’s center to the vehicle, V is the velocity of the vehicle, and μ is the gravitational constant designed as 3.98606e14. Based on a sphere and irrational earth mode, the 3-dimensional equations of motion [38] are shown as follows:
(1)θ′=−cosγcosψrcosϕD,(2)ϕ′=−cosγsinψrD,(3)r′=−sinγD,(4)ψ′=cosψtanϕcosγrD+LsinσV2cosγD,(5)γ′=g−V2/rcosγV2D−LcosσV2D,where θ is the longitude, ϕ is the latitude, ψ is the heading angle, γ is the flight path angle, and σ is the bank angle.The acceleration of gravityg is modeled as follows:
(6)g=μr2.The aerodynamic lift accelerationL and the aerodynamic drag acceleration D are given by
(7)L=qV2⋅S⋅CLα,Mam,(8)D=qV2⋅S⋅CDα,Mam,where q=ρV2/2 is the dynamic pressure, S is the reference area of the vehicle, CL is the aerodynamic lift coefficient, CD is the aerodynamic drag coefficient, α 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ρ based on the exponential atmosphere model is given by
(9)ρ=ρ0e−h/hs,where ρ0=1.225 kg/m3 is the atmospheric density at sea level, h=r−r0 is the altitude, r0=6378136.3 m is the equatorial radius of the Earth, and hs=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
(10)Q̇=kQρV3.15≤Q̇max,(11)q=ρV22≤qmax,(12)n=Lcosσ+Dsinα≤nmax,where Q̇ is the heating rate, kQ is the thermal model constant, and n is the normal acceleration. Moreover, Q̇max, qmax, and nmax are the maximum values.The terminal constraints for the altitude and the velocity are given by(13)hf−Δh≤hEf≤hf+Δh,(14)Vf−ΔV≤VEf≤Vf+ΔV,where Δh and ΔV are thresholds of the terminal constraint, rf and Vf are the expected terminal values, and Ef 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 Figure1. The center of the no-fly zone is designed on the surface of the earth at the position Nrnfz,θnfz,ϕnfz, where rnfz, θnfz, and ϕ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 1
Schematic of a no-fly zone model.Thus, in the geocentric coordinate system, coordinates of the positionN and the vehicle are given by
(15)xnfz=rnfzcosϕnfzcosθnfz,ynfz=rnfzcosϕnfzsinθnfz,znfz=rnfzsinϕnfz,(16)x=rcosϕcosθ,y=rcosϕsinθ,z=rsinϕ,where xnfz, ynfz, and znfz 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:
(17)x1y1z1=M1−90∘+αnfzM2φnfzM3θnfz−90∘xyz−0r00,where
(18)M1−90∘+αnfz=−sinαnfz0cosαnfz010−cosαnfz0−sinαnfz,(19)M2φnfz=1000cosφnfzsinφnfz0−sinφnfzcosφnfz,(20)M3θnfz−90∘=sinθnfz−cosθnfz0cosθnfzsinθnfz0001,x1, y1, and z1 are coordinates of the vehicle in the NUE coordinate system, and α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
(21)z12a2+x12b2≥1,y1≥0,where a and b are the long and short half-axes of the ellipse, respectively.


## 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 energyE is an appropriate independent variable, which is determined by E=μ/r−V2/2, where r is the radial distance from the Earth’s center to the vehicle, V is the velocity of the vehicle, and μ is the gravitational constant designed as 3.98606e14. Based on a sphere and irrational earth mode, the 3-dimensional equations of motion [38] are shown as follows:
(1)θ′=−cosγcosψrcosϕD,(2)ϕ′=−cosγsinψrD,(3)r′=−sinγD,(4)ψ′=cosψtanϕcosγrD+LsinσV2cosγD,(5)γ′=g−V2/rcosγV2D−LcosσV2D,where θ is the longitude, ϕ is the latitude, ψ is the heading angle, γ is the flight path angle, and σ is the bank angle.The acceleration of gravityg is modeled as follows:
(6)g=μr2.The aerodynamic lift accelerationL and the aerodynamic drag acceleration D are given by
(7)L=qV2⋅S⋅CLα,Mam,(8)D=qV2⋅S⋅CDα,Mam,where q=ρV2/2 is the dynamic pressure, S is the reference area of the vehicle, CL is the aerodynamic lift coefficient, CD is the aerodynamic drag coefficient, α 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ρ based on the exponential atmosphere model is given by
(9)ρ=ρ0e−h/hs,where ρ0=1.225 kg/m3 is the atmospheric density at sea level, h=r−r0 is the altitude, r0=6378136.3 m is the equatorial radius of the Earth, and hs=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
(10)Q̇=kQρV3.15≤Q̇max,(11)q=ρV22≤qmax,(12)n=Lcosσ+Dsinα≤nmax,where Q̇ is the heating rate, kQ is the thermal model constant, and n is the normal acceleration. Moreover, Q̇max, qmax, and nmax are the maximum values.The terminal constraints for the altitude and the velocity are given by(13)hf−Δh≤hEf≤hf+Δh,(14)Vf−ΔV≤VEf≤Vf+ΔV,where Δh and ΔV are thresholds of the terminal constraint, rf and Vf are the expected terminal values, and Ef 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 Figure1. The center of the no-fly zone is designed on the surface of the earth at the position Nrnfz,θnfz,ϕnfz, where rnfz, θnfz, and ϕ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 1
Schematic of a no-fly zone model.Thus, in the geocentric coordinate system, coordinates of the positionN and the vehicle are given by
(15)xnfz=rnfzcosϕnfzcosθnfz,ynfz=rnfzcosϕnfzsinθnfz,znfz=rnfzsinϕnfz,(16)x=rcosϕcosθ,y=rcosϕsinθ,z=rsinϕ,where xnfz, ynfz, and znfz 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:
(17)x1y1z1=M1−90∘+αnfzM2φnfzM3θnfz−90∘xyz−0r00,where
(18)M1−90∘+αnfz=−sinαnfz0cosαnfz010−cosαnfz0−sinαnfz,(19)M2φnfz=1000cosφnfzsinφnfz0−sinφnfzcosφnfz,(20)M3θnfz−90∘=sinθnfz−cosθnfz0cosθnfzsinθnfz0001,x1, y1, and z1 are coordinates of the vehicle in the NUE coordinate system, and α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
(21)z12a2+x12b2≥1,y1≥0,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σ. 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 Lcosσ by (7). The amplitude of ψ′ also can be obtained by (4), in which Lsinσ is given by
(22)Lsinσ=L2−Lcosσ2.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
(23)D′=−r0hs−2V2R2sinγ−2DV2+DCDCD′.Taking the derivative ofD′ with respect to E yields
(24)D″=4R2V2V′V+R′Rsinγ+−r0hs−2R2V2γ′cosγ+4DV3V′−2V2D′+CD′CDD′+DCD″CD−C′D2CD2.By substituting (5) into (24), we obtain Lcosσ=D″−b/a, which gives rise to
(25)σ=arccosD″−abL,where
(26)a=DCD″CD′−C′D2CD2+D′CD′CD+2V2−4DV4+1DV21hs+2gV2g−V2r,(27)b=−1V41hs+2gV2.A linear feedback control law is used to track drag acceleration profiles, as follows:(28)Dref″−D″+2ξωnDref′−D′+ωn2Dref−D=0,where ξ is the damp coefficient, ωn is the natural frequency of the tracking law, D is the real drag acceleration, and Dref is the reference drag acceleration.We can obtainD″ by (28) and substitute it into (25) to obtain the drag acceleration profile tracking law denoted as
(29)σ=arccosDref″+2ξωnDref′−D′+ωn2Dref−D−a/bL.In this paper, simulations are carried out on the Common Aero Vehicle (CAV) [43]. Let the input reference drag acceleration Dref=13 m/s2, the initial altitude h0=35 km, the attack angle α=20∘, 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 ξ and a constant natural frequency ω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 ξ=2 and ω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 2
Drag acceleration command tracking process withωn=5e−6.Figure 3
Drag acceleration command tracking process withξ=2.Figure 4
Drag acceleration command tracking process withξ=2 and ωn=5e−6.


### 3.2. Equilibrium Glide Condition

We assumed that the influence of the angle of attack on the range is ignored, and the angle of attack profile is determined by the Mach number, as follows:(30)α=αmax,Maα1≤MaαmaxL/D−αmaxMa−Maα1/Maα2−Maα1+αmax,Maα2≤Ma≤Maα1αmaxL/D,Ma≤Maα2,where αmax and αmaxL/D are the maximum angle of attack and the maximum lift drag ratio angle of attack. Maα1 and Maα2 are the turning points of the attack profile according to the vehicle parameters, as shown in Figure 5.Figure 5
Profile of the angle of attack.To achieve a smooth maximum downrange trajectory, we use the equilibrium glide condition by settingγ≤ε, where ε 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 γref, the feedback control law can be written as
(31)γref′−γ′+Kγγref−γ=0,where Kγ is a tracking coefficient.Substituting (5) into (31), we have
(32)LDcosσ=KγVDγref−γ+1Dg−V2rcosγ.To reduce bounces of the trajectory, we set the reference path angleγref=‐0.01∘. 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∘ at an appropriate switching time to satisfy the terminal velocity constraint.A simulation is carried out with parameters in section4 to study the effect of different values for the tracking coefficient Kγ 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γ=1 in this paper.Figure 6
Effect of differentKγ on the equilibrium glide condition.


### 3.3. Path Constraint Translation

The heat rate, dynamic pressure, and load constraints are translated into the corresponding drag profiles. By combining these drag acceleration profiles, the minimum downrange drag acceleration profile can be determined. We take the dynamic pressure constraint as an example to explain the transformation process.(1)
According to the initial and terminal states, the initial and terminal energy are obtained byE=μ/r−V2/2. We discretize the independent variable E into N intervals by ΔE=Ef−E0/N and obtain N+1 node points denoted as E0,E1,⋯,EN(2)
The dynamic pressure constraint is transformed into an equality constraint:f=q−qmax(3)
The values ofr and V satisfying E=Ei and f=0 are determined by the secant method [44]. The detailed iterative process of the secant method is as follows:
(a)
In thekth iteration, rk−1 is the radial distance from the center of the Earth to the vehicle for the last iteration, Vk−1 is obtained by E=μ/rk−1−Vk−12/2, and ρk−1 is obtained by (6). Especially in the first iteration, we can calculate these parameters according to the initial state of the vehicle. Based on the above definition, the (k−1)th equality constraint becomesfk−1=qk−1−qmax(b)
We setrk=rk−1−100, and the kth equality constraint fk is given in the same way as fk−1(c)
We use a simple secant method to achieverk+1=rk−fkrk−rk−1/fk−fk−1(d)
We check whether the convergence condition is met as follows:(33)rk+1−rk−1<ς,where ζ is a sufficiently small vector. If the condition in step d is met, we go to step e; otherwise, we go back to step a and start the next iteration using k=k+1(e)
Oncerk+1 is determined, Vk+1 is also determined by the energy formula. A converged solution rk+1,Vk+1 satisfying E=Ei and f=0 is found(4)
By performing iterative processes for each nodeE0,E1,⋯,EN, we can obtain the minimum downrange drag acceleration profileThe magnitude of the bank angle is constrained to0∘ at an appropriate switching time to satisfy the terminal velocity constraint. The path constraints are translated into the minimum downrange drag acceleration profile, as shown in Figure 7. Since the minimum downrange drag acceleration profile cannot be differentiable at switching points where Dref′ and Dref″ cannot be obtained, a cubic spline is used to realize the second derivative of the drag acceleration profile continuously at the switching points.Figure 7
Maximum downrange and minimum downrange drag profiles.For profiles between the maximum downrange and minimum downrange drag acceleration profiles, the interpolation equation is given by(34)D=Dmin+i2n2Dmax−Dmin,where Dmax and Dmin are the maximum and minimum drag acceleration profiles corresponding to the minimum and maximum range, respectively. i presents the ith 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ψmax1=ψ0 and ψmax2=ψegcf, where ψ0 is the initial entry heading angle and ψegcf is the terminal heading angle of the maximum downrange trajectory; thus, the supplementary terminal constraint ψmax is given by
(35)ψmax=ψmax1+iψ2nψ2ψmax1−ψmax2,where iψ represents the ith supplementary trajectory and nψ is the number of supplementary trajectories. Moreover, another supplementary terminal constraint is given by ψmin=ψmax−δψ, and δψ is the threshold of the heading angle.Similarly, for minimum downrange trajectories, we apply an interpolation between the heading anglesψmin1=ψ0 and ψmin2=ψconsf, where ψconsf is the terminal heading angle of the minimum downrange trajectory. Thus, the supplementary terminal constraint ψmin is given by
(36)ψmin=ψmin1+iψ2nψ2ψmin1−ψmin2.Moreover, another supplementary terminal constraint is given byψmax=ψmin+δψ.For the supplementary terminal constraints, supplementary trajectories are designed by adjusting the sign of the bank angle, as given by(37)signσi=sgnψmax−ψ,ψ>ψmaxsgnψmin−ψ,ψ<ψminsgnσi−1,ψmin≤ψ≤ψmax,where σi−1 is the bank angle at the previous moment and σi−1 is the bank angle at this moment. δψ is the threshold of the heading angle. If the value of δψ 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 Section2.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=Rz and b=Rz in (21), in which Rz is the radius of a circular no-fly zone. Assuming that the center of the no-fly zone Nθnfz,ϕnfz and the position of the vehicle Mθ,ϕ are in the same horizontal plane, the relationship between N and M is shown in Figure 8.Figure 8
Relative relationship between the center of the no-fly zone and the vehicle.MA is the tangent from M to the no-fly zone, and A is the tangent point. ψ¯A is the angle between the tangent MA and the north, and ψ¯Z is the angle between MN and the north. ψ¯ is the redefined heading angle, where ψ¯=‐ψ. ψ¯Z is solved by the following equation:
(38)ψ¯Z=arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ,ϕ<ϕnfz,ψ¯Z=π−arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ,ϕ≥ϕnfz.If the vehicle passes the no-fly zone from the right side, the boundary of the heading angle is given by(39)ψ¯min=ψ¯A,ψ¯max=ψ¯A+δψ¯,where
(40)ψ¯A=arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ+arcsinRZMN,δψ¯ 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 δψ¯ is the same as δψ 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(41)ψ¯min=ψ¯A−δψ¯,ψ¯max=ψ¯A,where
(42)ψ¯A=arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ‐arcsinRZMN.The sign of the bank angle is determined similarly to the equation in Section3.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 asJ, 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 J1 and J4. The trajectory is denoted as j2 or j4 when it corresponds to the positive or negative bank angle tracking from the maximum downrange drag acceleration profile. Similarly, the trajectory is denoted as j1 or j5 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 ψ0 corresponding to the maximum downrange drag acceleration profile is denoted as j3 and that corresponding to the minimum downrange drag acceleration profile is denoted as j6(b)
Supplementary landing footprints: we can obtain the set of interpolated landing footprintsJ2 or J3 constrained to the supplementary terminal heading angle by (35), in which the value of ψegcf is ψ2 (the terminal heading angle of j2) or ψ4 (the terminal heading angle of j4). The sign of the bank angle is given by (37). Similarly, we can obtain the set of interpolated landing footprints J5 or J6 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 footprintsJ1−J2−J3−J4−J5−J6 successively. The sets of landing footprints correspond to Figure 9Case 2.
(a)
Initial landing footprints: sets of landing footprintsJ2, J3, J4, and J5 are obtained by steps a and b in Case 1. For trajectories not affected by the no-fly zone in J1, we denote the trajectory closest to the no-fly zone from the left as j7 and closest from the right as j11. Landing footprints calculated from j7 to j2 are denoted as J7, and landing footprints from j11 and j1 are denoted as J16(b)
Modified landing footprints for the left side: use the left flying-around strategy in Section3.5 for trajectories between j7 and j1 and denote these landing footprints as J8. The terminal heading angles of the flying-around trajectories j8 and j9 are denoted as ψ8 and ψ9, corresponding to the original trajectories j7 and j1. Where the trajectory j9 and the no-fly zone intersect at the point of tangency, we denote the heading angle here as ψ10. The interpolation between ψ9 and ψ10 is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints J9 by the same method used in Section 3.4. Similarly, we can obtain set J10 by the interpolation of ψ0 and ψ10(c)
Modified landing footprints for the right side: landing footprints calculated fromj11 to j1 are denoted as J16. Use the right flying-around strategy in Section 3.5 for the trajectories between j11 and j1 and denote these landing footprints as J13. The terminal heading angle of flying-around trajectory j13 or j14 is denoted as ψ13 or ψ14, corresponding to the original trajectory j11 or j1. The terminal heading angle of the trajectory j11 or j1 is denoted as ψ11 or ψ1. Where the trajectory j13 or j14 and the no-fly zone intersect at a point of tangency, denote the heading angle here as ψ12 or ψ15. The interpolation between ψ11 and ψ12 is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints J11 by the same method used in Section 3.4. Similarly, we can obtain set J12, J14, or J15 by the interpolation of ψ12 and ψ13, ψ14, and ψ15, or ψ15 and ψ1(d)
Final landing footprint: the area of landing footprints consists of two closed figures constructed by connecting setsJ2−J3−J4−J5−J10−J9−J8−J7 and J16−J11−J12−J13−J14−J15 successively. The sets of landing footprints correspond to Figure 10Case 3.
(a)
Initial landing footprints: sets of landing footprintsJ2, J3, J4, J5, J6, J7, and J16 are obtained through steps a and b in Case 2(b)
Modified landing footprints for the left side: we can obtain setJ8 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 J8. In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in J8 is denoted as j9(c)
Modified landing footprints for the right side: the trajectoriesj11, j12, j13, set J11, and set J12 are given by the same step in Case 2. We can obtain set J13 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 J13. In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in J13 is denoted as j14(d)
Part of the no-fly zone: select pointsp1 and p2 closest to the landing footprints of j9 and j14 on the edge of the no-fly zone and denote points between p1 and p2 as P1(e)
Final landing footprint: we can obtain the final landing footprints by connecting sets of landing footprintsJ2−J3−J4−J5−J6−J16−J11−J12−J13−P1−J8−J7 successively. The sets of landing footprints correspond to Figure 11Case 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 Figure12.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 asj17 and that closest from the right side is denoted as j18. The landing footprints of the trajectories between j6 and j17 are denoted as J17, and the landing footprints of the trajectories between j18 and j1 are denoted as J18(e)
Part of the no-fly zone: select pointsp1 and p3 closest to the landing footprints of j9 and j17 on the edge of the no-fly zone and denote the points between p1 and p3 as P1. Select points p2 and p4 closest to the landing footprints of j14 and j18 on the edge of the no-fly zone and denote the points between p2 and p4 as P2(f)
Final landing footprint: the area of the landing footprints consists of two closed figures constructed by connecting setsJ2−J3−J4−J5−J17−P1−J8−J7 and J16−J11−J12−J13−P2−J18 successively
The above solution is based on the assumption that the no-fly zone affects sets of landing footprintsJ1, J4, J5, and J6. If the no-fly zone has an impact on J2 or J3, 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 9
Landing footprints for Case 1.Figure 10
Landing footprints for Case 2.Figure 11
Landing footprints for Case 3.Figure 12
Landing footprints for Case 4.Figure 13
Landing footprints for Case 5.


## 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
(23)D′=−r0hs−2V2R2sinγ−2DV2+DCDCD′.Taking the derivative ofD′ with respect to E yields
(24)D″=4R2V2V′V+R′Rsinγ+−r0hs−2R2V2γ′cosγ+4DV3V′−2V2D′+CD′CDD′+DCD″CD−C′D2CD2.By substituting (5) into (24), we obtain Lcosσ=D″−b/a, which gives rise to
(25)σ=arccosD″−abL,where
(26)a=DCD″CD′−C′D2CD2+D′CD′CD+2V2−4DV4+1DV21hs+2gV2g−V2r,(27)b=−1V41hs+2gV2.A linear feedback control law is used to track drag acceleration profiles, as follows:(28)Dref″−D″+2ξωnDref′−D′+ωn2Dref−D=0,where ξ is the damp coefficient, ωn is the natural frequency of the tracking law, D is the real drag acceleration, and Dref is the reference drag acceleration.We can obtainD″ by (28) and substitute it into (25) to obtain the drag acceleration profile tracking law denoted as
(29)σ=arccosDref″+2ξωnDref′−D′+ωn2Dref−D−a/bL.In this paper, simulations are carried out on the Common Aero Vehicle (CAV) [43]. Let the input reference drag acceleration Dref=13 m/s2, the initial altitude h0=35 km, the attack angle α=20∘, 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 ξ and a constant natural frequency ω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 ξ=2 and ω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 2
Drag acceleration command tracking process withωn=5e−6.Figure 3
Drag acceleration command tracking process withξ=2.Figure 4
Drag acceleration command tracking process withξ=2 and ωn=5e−6.


## 3.2. Equilibrium Glide Condition

We assumed that the influence of the angle of attack on the range is ignored, and the angle of attack profile is determined by the Mach number, as follows:(30)α=αmax,Maα1≤MaαmaxL/D−αmaxMa−Maα1/Maα2−Maα1+αmax,Maα2≤Ma≤Maα1αmaxL/D,Ma≤Maα2,where αmax and αmaxL/D are the maximum angle of attack and the maximum lift drag ratio angle of attack. Maα1 and Maα2 are the turning points of the attack profile according to the vehicle parameters, as shown in Figure 5.Figure 5
Profile of the angle of attack.To achieve a smooth maximum downrange trajectory, we use the equilibrium glide condition by settingγ≤ε, where ε 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 γref, the feedback control law can be written as
(31)γref′−γ′+Kγγref−γ=0,where Kγ is a tracking coefficient.Substituting (5) into (31), we have
(32)LDcosσ=KγVDγref−γ+1Dg−V2rcosγ.To reduce bounces of the trajectory, we set the reference path angleγref=‐0.01∘. 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∘ at an appropriate switching time to satisfy the terminal velocity constraint.A simulation is carried out with parameters in section4 to study the effect of different values for the tracking coefficient Kγ 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γ=1 in this paper.Figure 6
Effect of differentKγ on the equilibrium glide condition.


## 3.3. Path Constraint Translation

The heat rate, dynamic pressure, and load constraints are translated into the corresponding drag profiles. By combining these drag acceleration profiles, the minimum downrange drag acceleration profile can be determined. We take the dynamic pressure constraint as an example to explain the transformation process.(1)
According to the initial and terminal states, the initial and terminal energy are obtained byE=μ/r−V2/2. We discretize the independent variable E into N intervals by ΔE=Ef−E0/N and obtain N+1 node points denoted as E0,E1,⋯,EN(2)
The dynamic pressure constraint is transformed into an equality constraint:f=q−qmax(3)
The values ofr and V satisfying E=Ei and f=0 are determined by the secant method [44]. The detailed iterative process of the secant method is as follows:
(a)
In thekth iteration, rk−1 is the radial distance from the center of the Earth to the vehicle for the last iteration, Vk−1 is obtained by E=μ/rk−1−Vk−12/2, and ρk−1 is obtained by (6). Especially in the first iteration, we can calculate these parameters according to the initial state of the vehicle. Based on the above definition, the (k−1)th equality constraint becomesfk−1=qk−1−qmax(b)
We setrk=rk−1−100, and the kth equality constraint fk is given in the same way as fk−1(c)
We use a simple secant method to achieverk+1=rk−fkrk−rk−1/fk−fk−1(d)
We check whether the convergence condition is met as follows:(33)rk+1−rk−1<ς,where ζ is a sufficiently small vector. If the condition in step d is met, we go to step e; otherwise, we go back to step a and start the next iteration using k=k+1(e)
Oncerk+1 is determined, Vk+1 is also determined by the energy formula. A converged solution rk+1,Vk+1 satisfying E=Ei and f=0 is found(4)
By performing iterative processes for each nodeE0,E1,⋯,EN, we can obtain the minimum downrange drag acceleration profileThe magnitude of the bank angle is constrained to0∘ at an appropriate switching time to satisfy the terminal velocity constraint. The path constraints are translated into the minimum downrange drag acceleration profile, as shown in Figure 7. Since the minimum downrange drag acceleration profile cannot be differentiable at switching points where Dref′ and Dref″ cannot be obtained, a cubic spline is used to realize the second derivative of the drag acceleration profile continuously at the switching points.Figure 7
Maximum downrange and minimum downrange drag profiles.For profiles between the maximum downrange and minimum downrange drag acceleration profiles, the interpolation equation is given by(34)D=Dmin+i2n2Dmax−Dmin,where Dmax and Dmin are the maximum and minimum drag acceleration profiles corresponding to the minimum and maximum range, respectively. i presents the ith 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ψmax1=ψ0 and ψmax2=ψegcf, where ψ0 is the initial entry heading angle and ψegcf is the terminal heading angle of the maximum downrange trajectory; thus, the supplementary terminal constraint ψmax is given by
(35)ψmax=ψmax1+iψ2nψ2ψmax1−ψmax2,where iψ represents the ith supplementary trajectory and nψ is the number of supplementary trajectories. Moreover, another supplementary terminal constraint is given by ψmin=ψmax−δψ, and δψ is the threshold of the heading angle.Similarly, for minimum downrange trajectories, we apply an interpolation between the heading anglesψmin1=ψ0 and ψmin2=ψconsf, where ψconsf is the terminal heading angle of the minimum downrange trajectory. Thus, the supplementary terminal constraint ψmin is given by
(36)ψmin=ψmin1+iψ2nψ2ψmin1−ψmin2.Moreover, another supplementary terminal constraint is given byψmax=ψmin+δψ.For the supplementary terminal constraints, supplementary trajectories are designed by adjusting the sign of the bank angle, as given by(37)signσi=sgnψmax−ψ,ψ>ψmaxsgnψmin−ψ,ψ<ψminsgnσi−1,ψmin≤ψ≤ψmax,where σi−1 is the bank angle at the previous moment and σi−1 is the bank angle at this moment. δψ is the threshold of the heading angle. If the value of δψ 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 Section2.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=Rz and b=Rz in (21), in which Rz is the radius of a circular no-fly zone. Assuming that the center of the no-fly zone Nθnfz,ϕnfz and the position of the vehicle Mθ,ϕ are in the same horizontal plane, the relationship between N and M is shown in Figure 8.Figure 8
Relative relationship between the center of the no-fly zone and the vehicle.MA is the tangent from M to the no-fly zone, and A is the tangent point. ψ¯A is the angle between the tangent MA and the north, and ψ¯Z is the angle between MN and the north. ψ¯ is the redefined heading angle, where ψ¯=‐ψ. ψ¯Z is solved by the following equation:
(38)ψ¯Z=arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ,ϕ<ϕnfz,ψ¯Z=π−arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ,ϕ≥ϕnfz.If the vehicle passes the no-fly zone from the right side, the boundary of the heading angle is given by(39)ψ¯min=ψ¯A,ψ¯max=ψ¯A+δψ¯,where
(40)ψ¯A=arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ+arcsinRZMN,δψ¯ 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 δψ¯ is the same as δψ 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(41)ψ¯min=ψ¯A−δψ¯,ψ¯max=ψ¯A,where
(42)ψ¯A=arctansinθnfz−θcosϕtanϕnfz−sinϕcosθnfz−θ‐arcsinRZMN.The sign of the bank angle is determined similarly to the equation in Section3.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 asJ, 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 J1 and J4. The trajectory is denoted as j2 or j4 when it corresponds to the positive or negative bank angle tracking from the maximum downrange drag acceleration profile. Similarly, the trajectory is denoted as j1 or j5 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 ψ0 corresponding to the maximum downrange drag acceleration profile is denoted as j3 and that corresponding to the minimum downrange drag acceleration profile is denoted as j6(b)
Supplementary landing footprints: we can obtain the set of interpolated landing footprintsJ2 or J3 constrained to the supplementary terminal heading angle by (35), in which the value of ψegcf is ψ2 (the terminal heading angle of j2) or ψ4 (the terminal heading angle of j4). The sign of the bank angle is given by (37). Similarly, we can obtain the set of interpolated landing footprints J5 or J6 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 footprintsJ1−J2−J3−J4−J5−J6 successively. The sets of landing footprints correspond to Figure 9Case 2.
(a)
Initial landing footprints: sets of landing footprintsJ2, J3, J4, and J5 are obtained by steps a and b in Case 1. For trajectories not affected by the no-fly zone in J1, we denote the trajectory closest to the no-fly zone from the left as j7 and closest from the right as j11. Landing footprints calculated from j7 to j2 are denoted as J7, and landing footprints from j11 and j1 are denoted as J16(b)
Modified landing footprints for the left side: use the left flying-around strategy in Section3.5 for trajectories between j7 and j1 and denote these landing footprints as J8. The terminal heading angles of the flying-around trajectories j8 and j9 are denoted as ψ8 and ψ9, corresponding to the original trajectories j7 and j1. Where the trajectory j9 and the no-fly zone intersect at the point of tangency, we denote the heading angle here as ψ10. The interpolation between ψ9 and ψ10 is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints J9 by the same method used in Section 3.4. Similarly, we can obtain set J10 by the interpolation of ψ0 and ψ10(c)
Modified landing footprints for the right side: landing footprints calculated fromj11 to j1 are denoted as J16. Use the right flying-around strategy in Section 3.5 for the trajectories between j11 and j1 and denote these landing footprints as J13. The terminal heading angle of flying-around trajectory j13 or j14 is denoted as ψ13 or ψ14, corresponding to the original trajectory j11 or j1. The terminal heading angle of the trajectory j11 or j1 is denoted as ψ11 or ψ1. Where the trajectory j13 or j14 and the no-fly zone intersect at a point of tangency, denote the heading angle here as ψ12 or ψ15. The interpolation between ψ11 and ψ12 is regarded as the terminal heading angle constraint. Thus, we can obtain the set of landing footprints J11 by the same method used in Section 3.4. Similarly, we can obtain set J12, J14, or J15 by the interpolation of ψ12 and ψ13, ψ14, and ψ15, or ψ15 and ψ1(d)
Final landing footprint: the area of landing footprints consists of two closed figures constructed by connecting setsJ2−J3−J4−J5−J10−J9−J8−J7 and J16−J11−J12−J13−J14−J15 successively. The sets of landing footprints correspond to Figure 10Case 3.
(a)
Initial landing footprints: sets of landing footprintsJ2, J3, J4, J5, J6, J7, and J16 are obtained through steps a and b in Case 2(b)
Modified landing footprints for the left side: we can obtain setJ8 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 J8. In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in J8 is denoted as j9(c)
Modified landing footprints for the right side: the trajectoriesj11, j12, j13, set J11, and set J12 are given by the same step in Case 2. We can obtain set J13 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 J13. In particular, the trajectory corresponding to the landing footprint closest to the no-fly zone in J13 is denoted as j14(d)
Part of the no-fly zone: select pointsp1 and p2 closest to the landing footprints of j9 and j14 on the edge of the no-fly zone and denote points between p1 and p2 as P1(e)
Final landing footprint: we can obtain the final landing footprints by connecting sets of landing footprintsJ2−J3−J4−J5−J6−J16−J11−J12−J13−P1−J8−J7 successively. The sets of landing footprints correspond to Figure 11Case 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 Figure12.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 asj17 and that closest from the right side is denoted as j18. The landing footprints of the trajectories between j6 and j17 are denoted as J17, and the landing footprints of the trajectories between j18 and j1 are denoted as J18(e)
Part of the no-fly zone: select pointsp1 and p3 closest to the landing footprints of j9 and j17 on the edge of the no-fly zone and denote the points between p1 and p3 as P1. Select points p2 and p4 closest to the landing footprints of j14 and j18 on the edge of the no-fly zone and denote the points between p2 and p4 as P2(f)
Final landing footprint: the area of the landing footprints consists of two closed figures constructed by connecting setsJ2−J3−J4−J5−J17−P1−J8−J7 and J16−J11−J12−J13−P2−J18 successively
The above solution is based on the assumption that the no-fly zone affects sets of landing footprintsJ1, J4, J5, and J6. If the no-fly zone has an impact on J2 or J3, 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 9
Landing footprints for Case 1.Figure 10
Landing footprints for Case 2.Figure 11
Landing footprints for Case 3.Figure 12
Landing footprints for Case 4.Figure 13
Landing footprints for Case 5.


## 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 m2. The path constraints are set as Q̇max=2000 kW/m2, qmax=100 kPa, nmax=3 g, and kQ=1e‐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 Δh=100 m and ΔV=50 m/s. The angles of attack profile parameters are designed as αmax=20∘, αmaxL/D=20∘, Maα1=12, and Maα2=12. The tracking law coefficients for the drag acceleration command tracking process are designed ξ=1 and ωn=5e−6. Similarly, the tracking law coefficient for the equilibrium glide condition is designed Kγ=1. The threshold of the heading angle in Section 3.4 and Section 3.5 is designed as δψ=2∘. The number of drag acceleration profiles is designed as n=10, and the number of supplementary trajectories is designed as nψ=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.
θ∘ϕ∘h0kmψ∘γ∘Vm/s00600-15000


### 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. Figures14 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 14
Minimum drag acceleration profile tracking ability.Figure 15
Maximum drag acceleration profile tracking ability.Table2 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 methodComputation time of the Gauss pseudospectral methodEquilibrium glide condition: 0.14 sPath constraint translation: 0.23Each tracking trajectory: 0.16 s–Total time: 13.42 sTotal time: 429.93 sFigure16 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 16
Landing footprints for Case 1 (method 1: the Gauss pseudospectral method; 2: the EAGLE method).Figure 17
Altitude profiles for Case 1.Figure 18
Drag acceleration profiles for Case 1.Figure 19
Lift acceleration profiles for Case 1.Figure 20
Bank angle profiles for Case 1.Figure 21
Angle of attack profiles for Case 1.Figure 22
Flight path angle profiles for Case 1.Figure 23
Heading angle profiles for Case 1.Figure24 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 24
Dispersion of the terminal altitude.Figure 25
Dispersion of the terminal velocity.Figure 26
Normal acceleration profiles for Case 1.Figure 27
Dynamic pressure profiles for Case 1.Figure 28
Heating rate profiles for Case 1.


### 4.2. With No-Fly Zone

In this section, the no-fly zone constraint is added. Centers of designed no-fly zones are set as1.4∘,10∘, 3∘,12.5∘, 2∘,10.6∘, and 2∘,11.1∘, 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 αnfz=0∘. 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 29
Landing footprints for Case 2.Figure 30
Landing footprints for Case 3.Figure 31
Landing footprints for Case 4.Figure 32
Landing footprints for Case 5.Table 3
Computation times of different cases.
Different casesNumber of trajectories not affected by the no-fly zoneComputation time of trajectories not affected by the no-fly zoneNumber of trajectories affected by the no-fly zoneComputation time of trajectories affected by the no-fly zoneTotal computation timeCase 26911.21 s467.96 s19.18 sCase 37712.51 s213.64 s16.15 sCase 47311.86 s335.71 s17.57 sCase 57311.86 s305.19 s17.06 s


## 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. Figures14 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 14
Minimum drag acceleration profile tracking ability.Figure 15
Maximum drag acceleration profile tracking ability.Table2 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 methodComputation time of the Gauss pseudospectral methodEquilibrium glide condition: 0.14 sPath constraint translation: 0.23Each tracking trajectory: 0.16 s–Total time: 13.42 sTotal time: 429.93 sFigure16 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 16
Landing footprints for Case 1 (method 1: the Gauss pseudospectral method; 2: the EAGLE method).Figure 17
Altitude profiles for Case 1.Figure 18
Drag acceleration profiles for Case 1.Figure 19
Lift acceleration profiles for Case 1.Figure 20
Bank angle profiles for Case 1.Figure 21
Angle of attack profiles for Case 1.Figure 22
Flight path angle profiles for Case 1.Figure 23
Heading angle profiles for Case 1.Figure24 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 24
Dispersion of the terminal altitude.Figure 25
Dispersion of the terminal velocity.Figure 26
Normal acceleration profiles for Case 1.Figure 27
Dynamic pressure profiles for Case 1.Figure 28
Heating rate profiles for Case 1.


## 4.2. With No-Fly Zone

In this section, the no-fly zone constraint is added. Centers of designed no-fly zones are set as1.4∘,10∘, 3∘,12.5∘, 2∘,10.6∘, and 2∘,11.1∘, 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 αnfz=0∘. 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 29
Landing footprints for Case 2.Figure 30
Landing footprints for Case 3.Figure 31
Landing footprints for Case 4.Figure 32
Landing footprints for Case 5.Table 3
Computation times of different cases.
Different casesNumber of trajectories not affected by the no-fly zoneComputation time of trajectories not affected by the no-fly zoneNumber of trajectories affected by the no-fly zoneComputation time of trajectories affected by the no-fly zoneTotal computation timeCase 26911.21 s467.96 s19.18 sCase 37712.51 s213.64 s16.15 sCase 47311.86 s335.71 s17.57 sCase 57311.86 s305.19 s17.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.



---

*Source: 8827377-2021-01-21.xml*