# Adaptive Entry Guidance for Hypersonic Gliding Vehicles Using Analytic Feedback Control

Authors: Xunliang Yan; Peichen Wang; Shaokang Xu; Shumei Wang; Hao Jiang
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

This paper presents an adaptive, simple, and effective guidance approach for hypersonic entry vehicles with high lift-to-drag (L/D) ratios (e.g., hypersonic gliding vehicles). The core of the constrained guidance approach is a closed-form, easily obtained, and computationally efficient feedback control law that yields the analytic bank command based on the well-known quasi-equilibrium glide condition (QEGC). The magnitude of the bank angle command consists of two parts, i.e., the baseline part and the augmented part, which are calculated analytically and successively. The baseline command is derived from the analytic relation between the range-to-go and the velocity to guarantee the range requirement. Then, the bank angle is augmented with the predictive altitude-rate feedback compensations that are represented by an analytic set of flight path angle needed for the terminal constraints. The inequality path constraints in the velocity-altitude space are translated into the velocity-dependent bounds for the magnitude of the bank angle based on the QEGC. The sign of the bank command is also analytically determined using an automated bank-reversal logic based on the dynamic adjustment criteria. Finally, a feasible three-degree-of-freedom (3DOF) entry flight trajectory is simultaneously generated by integrating with the real-time updated command. Because no iterations and no or few off-line parameter adjustments are required using almost all analytic processing, the algorithm provides remarkable simplicity, rapidity, and adaptability. A considerable range of entry flights using the vehicle data of the CAV-H is tested. Simulation results demonstrate the effectiveness and performance of the presented approach.


---

## Body

## 1. Introduction

Atmosphere entry flight is a critical phase of operation for the unpowered lifting hypersonic flight vehicles such as reusable launch vehicles (RLVs) and hypersonic gliding vehicles (HGVs). Entry trajectory generation and guidance are challenging and responsible for the success of entry flight. Therefore, extensive studies can be found in recent years [1, 2]. Currently, entry guidance methods can be divided into two categories: the standard trajectory guidance and the predictor-corrector guidance [2]. The standard trajectory guidance that is more mature and widely used includes two parts: trajectory planning and tracking. The entry trajectory planning is usually based on numerical optimization or numerical iteration methods which are usually time-consuming and laborious [3]. More specifically, a well-known approach, i.e., planning aerodynamic drag acceleration profile as the reference trajectory, typically implemented in the shuttle entry guidance and lately extended to other instances (e.g., Evolved Acceleration Guidance Logic for Entry, EAGLE), has proven to be very effective and successful, which becomes the baseline approach for many entry vehicles [1].Even though the shuttle entry guidance is successful [4], there have been several promising extensions and applications of drag profile approach over the years. These studies strive to improve the accuracy and computation time of the reference drag profile by simplifying or automating the drag profile design [5–9], investigate linear and nonlinear full-state feedback tracking laws [10–12], improve on both the above two aspects [13, 14], or enhance the lateral maneuverability dealing with geographic constraints [15–17]. On the whole, these efforts are still considered as the variants of shuttle entry guidance and defined as the standard trajectory guidance.It is adequate to have the drag-based trajectory generator on the ground for the lifting vehicles with a limited flight envelope and focused mission; onboard trajectory generation is still necessary for the second generation RLVs or HGVs to achieve aircraft-like operation. A far-reaching contribution proposed by Shen and Lu [18] is a cornerstone of onboard 3DOF trajectory generation, which is not based on the drag acceleration profile. This benchmark effort uses the so-called quasi-equilibrium glide condition (QEGC) [19–21], a frequently observed phenomenon in the hypersonic lifting flight of vehicles with moderate to higher L/D ratios, as the foundation for the rapid online design of a feasible entry trajectory subject to all common conditions, and effective and efficient enforcement of the inequality constraints. On the basis of the trajectory generator presented by Shen, an adaptive lateral guidance logic for determining when to perform bank-angle reversals in the most stressful scenarios is investigated in [22]. Zang et al. [23] presented an on-line guidance algorithm for high L/D hypersonic reentry vehicles using a plane-symmetry bank-to-turn control method that can generate a feasible trajectory at each guidance cycle.Besides, the classical predictor-corrector algorithms have evolved and emerged to show significant potential to disengage from any dependence on the separate preplanned reference trajectory and tracking laws [24–31]. The predictor-corrector algorithms are aimed at iteratively determining a complete feasible entry trajectory onboard based on the current condition and the desired target condition. Despite many advantages, a long-standing weakness of the predictor-corrector algorithm is the lack of effective and broadly applicable means to enforce inequality trajectory constraints such as those on the heating rate and aerodynamic load [25–27]. In order to address this issue, Xue and Lu [28] presented a highly effective algorithm to enforce common inequality entry trajectory constraints in a predictor-corrector algorithm by employing the QEGC. Furthermore, Lu [29–31] presented a unified predictor-corrector method for both low and high lifting vehicles, in which the enforcement of common trajectory constraints is conducted by an augmentation of altitude-rate feedback to the baseline algorithm based on the natural time-scale separation of the trajectory dynamics and a nonlinear predictive control technique.Obviously, almost all aforementioned entry guidance algorithms require conducting several or more numerical iterations, in which repeated integrations of the equations of entry motion are involved, so as to generate a set of guidance commands and a feasible entry trajectory satisfying all common constraints. The main weakness of the numerical iterations, however, is the lack of the convergence guarantee of the numerical process. Moreover, such one or more repeated integrations involved in iterations add a so severe computation burden that the onboard capability and the terminal precision will both degenerate. Abandoning numerical iterations and repeated integrations, Xu et al. [32] presented a novel quasi-equilibrium glide adaptive entry trajectory generation algorithm based on the predictor-corrector principle for hypersonic lifting vehicles. The trajectory is converted into a special form to obtain the closed-form solution with the analytically calculated angle of attack and bank. Pan et al. [33] presented a three-dimensional guidance algorithm on the basis of analytical predictions for the trajectory using Lyapunov’s artificial small parameter method. However, this algorithm is essentially one of the standard trajectory guidance algorithms that numerical iterations cannot be avoided.In this paper, we present a rapid, relatively simple, and effective approach of trajectory planning for entry vehicles (such as HGVs) with a high L/D ratio. This approach is inspired by the contribution in [32] but owns an essentially distinct algorithmic principle. Novel utilization of the QEGC is the cornerstone for this rapid planning algorithm for fully constrained, three-dimensional feasible entry trajectories. The primary commands are the fixed velocity-dependent angle of attack and the adjustable bank angle which is calculated analytically. The magnitude of bank angle command consists of two parts: the baseline part derived from the analytical relation between the range-to-go and the velocity, and the augmented part that is generated by using the predictive objective-oriented altitude-rate feedback compensations required for the desired set of flight path angle. This set of flight path angle, treated as the pesudocontrol, is simply and readily deduced using the analytic expressions relating the range-to-go to the desired terminal altitude and relating the desired terminal altitude to the predicted terminal velocity, respectively. The inequality path constraints in the velocity-altitude space are dramatically translated into the velocity-dependent bounds for the magnitude of the bank angle by the QEGC. The sign of the bank command is determined by an automated bank-reversal logic based on the approximate linearity and proportional property between the crossrange and the range-to-go. The over-correct scheme is utilized with a constant parameter and conservative criterion to ensure that the crossrange and heading error requirements are all satisfied at an acceptable expense of one or more additional bank reversals. A feasible 3DOF entry trajectory is simultaneously generated by integrating the real-time updated command. No iterations are required, and few off-line parameter adjustments are necessary with only one time’s integration conducted along the trajectory. A considerable range of entry flights using the vehicle data of the CAV-H is tested. Simulation results demonstrate the effectiveness and performance of the presented approach.


## 2. Entry Guidance Problem

### 2.1. Entry Dynamics

The dimensionless 3DOF equations of motion of a HGV over a spherical, rotating Earth are given by(1)ṙ=Vsinθ,(2)λ̇=Vcosθsinσ/rcosϕ,(3)ϕ̇=Vcosθcosσ/r,(4)V̇=−D−sinθ/r2+rωe2cosϕcosϕsinθ−rωe2cosϕsinϕcosσcosθ,(5)Vθ̇=Lcosν+V2−1/rcosθ/r+2ωeVcosϕsinσ+rωe2cosϕcosϕcosθ+sinθsinϕcosσ,(6)Vσ̇=Lsinν/cosθ+V2/rcosθsinσtanϕ+2ωeVsinϕ−cosϕcosσtanθ+rωe2/cosθcosϕsinϕsinσ,where r is the radial distance from the Earth center to the HGV, λ the longitude, ϕ the latitude, V the Earth-relative velocity, θ the flight path angle, ν the bank angle defined such that a bank to the right is positive, and σ the velocity azimuth angle (i.e., heading angle) measured clockwise from the North. ωe is the self-rotation rate of Earth. In the nondimensional form, length and time are normalized by the radius of the Earth R0 and tscale=R0/g0 with g0=9.81m/s2, respectively, thus leading to dimensionless velocity Vscale=R0g0 and angular rate ωscale=g0/R0. The differentiation is with respect to the dimensionless time τ=t/tscale. The terms D and L are dimensionless aerodynamic accelerations (in g0), i.e.,
(7)D=ρVVscale2SrefCD/2mg0,(8)L=ρVVscale2SrefCL/2mg0,where Sref is the reference area of the vehicle and m is the mass of the vehicle. CD and CL are the aerodynamic drag and lift coefficients as functions of α and Mach number. The atmospheric density ρ is modelled using the exponential equation
(9)ρ=ρ0e−h/hs,where ρ0 is the atmospheric density at the sea level, h=R0r−1 is the altitude, and hs an altitude constant.The angle of attackα is assumed to be a fixed velocity-dependent profile determined synthetically by thermal protection, range capability, and control constraints, whereas it is slightly adjustable for entry tracking guidance not concerned in this paper. The only adjustable trajectory command ν is to be determined by the guidance approach in the following sections. Thus, the dimensionless entry dynamics can be rewritten as
(10)ẋ=dx/dτ=fx,u,xτ0=x0,where the state vector x=r,λ,ϕ,V,θ,σT and the control vector u=ν,αT. xτ0=x0 presents the initial conditions. Note that initial conditions will be denoted with a subscript “0,” and then target conditions will be denoted with subscript “f” in the following sections.


### 2.2. Trajectory Constraints

The entry trajectory should start with the initial conditions at the entry interface and terminate with the desired target conditions to ensure that the next phase can be successfully conducted. The typical terminal constraints for entry flight are specified so that the trajectory reaches to a location with a desired distancesf (sf can be zero) from the target point at a specified final altitude rf and velocity Vf. That is,
(11)rτf=rf,(12)Vτf=Vf,(13)stogoλτf,ϕτf=sf,where stogo denotes the value of range-to-go from the current point to the target location. Introducing the energy-like parameter e=1/r−V2/2, the first two conditions in Eqs. (11) and (12) maybe combined to define a specified final energy
(14)ef=1/rf−Vf2/2.Alternatively, under the conventions of entry flight, entry terminates at the specified final energyef instead of rf or Vf. Thus, the terminal conditions in Eqs. (11)–(13) are translated into
(15)eτf=ef,(16)stogoef=sf.Considering that the final velocity vector may be directed at the target point with a given toleranceΔσf, the heading error, which is the difference between the velocity azimuth angle and the line-of-sight angle from the vehicle to the target point, is limited by
(17)Δστf≤Δσf.The common entry trajectory inequality path constraints for hypersonic glide, including those on the heating rate at a stagnation pointQ̇, aerodynamic load n, and dynamic pressure q, are expressed as
(18)Q̇=kQρVVscale3.15≤Q̇max,(19)n=L2+D2≤nmax,(20)q=0.5ρVVscale2≤qmax,where kQ is a vehicle-dependent constant and Q̇max, nmax, and qmax are vehicle-dependent peak constants as well, respectively. These three constraints are considered “hard” constraints to be enforced strictly.For HGVs with a high L/D ratio, another path constraint is the equilibrium glide constraint withθ=0, ν=νEG, and the Earth self-rotation ignored. Then, it is expressed as
(21)LcosνEG+V2−1/r1/r≥0,where νEG is a specified constant. The steady flight could not be maintained when this condition is violated, because the vehicle would not have enough lift to maintain its flight path angle so that the phugoid oscillations in altitudes will be reduced. Nevertheless, the violation of this condition would not pose a risk to the vehicle not similar to the above three “hard” constraints. Thus, this condition is referred to as a “soft” constraint that does not need to be enforced strictly.Considering the attitude control system capability and the nominal angle of attack profile, limits are placed on flight control authority according to(22)ν≤νmax,ν̇≤ν̇max.The entry guidance is to determine the control historyu=ν,αT so that the corresponding entry flight should satisfy all of the aforementioned constraints in terms of the 3DOF entry dynamics, endpoint boundary conditions, typical path inequality constraints, and control authority constraints. Accordingly, a feasible trajectory is generated, and the rapidness and reliability are pursued subsequently.


## 2.1. Entry Dynamics

The dimensionless 3DOF equations of motion of a HGV over a spherical, rotating Earth are given by(1)ṙ=Vsinθ,(2)λ̇=Vcosθsinσ/rcosϕ,(3)ϕ̇=Vcosθcosσ/r,(4)V̇=−D−sinθ/r2+rωe2cosϕcosϕsinθ−rωe2cosϕsinϕcosσcosθ,(5)Vθ̇=Lcosν+V2−1/rcosθ/r+2ωeVcosϕsinσ+rωe2cosϕcosϕcosθ+sinθsinϕcosσ,(6)Vσ̇=Lsinν/cosθ+V2/rcosθsinσtanϕ+2ωeVsinϕ−cosϕcosσtanθ+rωe2/cosθcosϕsinϕsinσ,where r is the radial distance from the Earth center to the HGV, λ the longitude, ϕ the latitude, V the Earth-relative velocity, θ the flight path angle, ν the bank angle defined such that a bank to the right is positive, and σ the velocity azimuth angle (i.e., heading angle) measured clockwise from the North. ωe is the self-rotation rate of Earth. In the nondimensional form, length and time are normalized by the radius of the Earth R0 and tscale=R0/g0 with g0=9.81m/s2, respectively, thus leading to dimensionless velocity Vscale=R0g0 and angular rate ωscale=g0/R0. The differentiation is with respect to the dimensionless time τ=t/tscale. The terms D and L are dimensionless aerodynamic accelerations (in g0), i.e.,
(7)D=ρVVscale2SrefCD/2mg0,(8)L=ρVVscale2SrefCL/2mg0,where Sref is the reference area of the vehicle and m is the mass of the vehicle. CD and CL are the aerodynamic drag and lift coefficients as functions of α and Mach number. The atmospheric density ρ is modelled using the exponential equation
(9)ρ=ρ0e−h/hs,where ρ0 is the atmospheric density at the sea level, h=R0r−1 is the altitude, and hs an altitude constant.The angle of attackα is assumed to be a fixed velocity-dependent profile determined synthetically by thermal protection, range capability, and control constraints, whereas it is slightly adjustable for entry tracking guidance not concerned in this paper. The only adjustable trajectory command ν is to be determined by the guidance approach in the following sections. Thus, the dimensionless entry dynamics can be rewritten as
(10)ẋ=dx/dτ=fx,u,xτ0=x0,where the state vector x=r,λ,ϕ,V,θ,σT and the control vector u=ν,αT. xτ0=x0 presents the initial conditions. Note that initial conditions will be denoted with a subscript “0,” and then target conditions will be denoted with subscript “f” in the following sections.


## 2.2. Trajectory Constraints

The entry trajectory should start with the initial conditions at the entry interface and terminate with the desired target conditions to ensure that the next phase can be successfully conducted. The typical terminal constraints for entry flight are specified so that the trajectory reaches to a location with a desired distancesf (sf can be zero) from the target point at a specified final altitude rf and velocity Vf. That is,
(11)rτf=rf,(12)Vτf=Vf,(13)stogoλτf,ϕτf=sf,where stogo denotes the value of range-to-go from the current point to the target location. Introducing the energy-like parameter e=1/r−V2/2, the first two conditions in Eqs. (11) and (12) maybe combined to define a specified final energy
(14)ef=1/rf−Vf2/2.Alternatively, under the conventions of entry flight, entry terminates at the specified final energyef instead of rf or Vf. Thus, the terminal conditions in Eqs. (11)–(13) are translated into
(15)eτf=ef,(16)stogoef=sf.Considering that the final velocity vector may be directed at the target point with a given toleranceΔσf, the heading error, which is the difference between the velocity azimuth angle and the line-of-sight angle from the vehicle to the target point, is limited by
(17)Δστf≤Δσf.The common entry trajectory inequality path constraints for hypersonic glide, including those on the heating rate at a stagnation pointQ̇, aerodynamic load n, and dynamic pressure q, are expressed as
(18)Q̇=kQρVVscale3.15≤Q̇max,(19)n=L2+D2≤nmax,(20)q=0.5ρVVscale2≤qmax,where kQ is a vehicle-dependent constant and Q̇max, nmax, and qmax are vehicle-dependent peak constants as well, respectively. These three constraints are considered “hard” constraints to be enforced strictly.For HGVs with a high L/D ratio, another path constraint is the equilibrium glide constraint withθ=0, ν=νEG, and the Earth self-rotation ignored. Then, it is expressed as
(21)LcosνEG+V2−1/r1/r≥0,where νEG is a specified constant. The steady flight could not be maintained when this condition is violated, because the vehicle would not have enough lift to maintain its flight path angle so that the phugoid oscillations in altitudes will be reduced. Nevertheless, the violation of this condition would not pose a risk to the vehicle not similar to the above three “hard” constraints. Thus, this condition is referred to as a “soft” constraint that does not need to be enforced strictly.Considering the attitude control system capability and the nominal angle of attack profile, limits are placed on flight control authority according to(22)ν≤νmax,ν̇≤ν̇max.The entry guidance is to determine the control historyu=ν,αT so that the corresponding entry flight should satisfy all of the aforementioned constraints in terms of the 3DOF entry dynamics, endpoint boundary conditions, typical path inequality constraints, and control authority constraints. Accordingly, a feasible trajectory is generated, and the rapidness and reliability are pursued subsequently.


## 3. Entry Guidance Algorithms

This section presents a simple, adaptive, and autonomous guidance algorithm of constrained entry hypersonic flight for HGVs. The algorithm tackles the problem in two steps: the longitudinal guidance and the lateral guidance. The longitudinal guidance generates the feasible magnitude of bank angle in real-time, while the lateral guidance determines the sign of bank angle by a simple but efficient bank reversal logic. Conducting simultaneously these two channels with the successive states updated, the set of closed-loop commanded bank angle with analytic feedback control laws are easily deduced. Accordingly, a feasible and applicable entry trajectory is generated by integrating the whole trajectory only once with a pleasing computation cost.


### 3.1. QEGC and Translation of Inequality Path Constraints

As mentioned above, an ingenious utilization of QEGC is the cornerstone for this entry guidance algorithm and translation of inequality path constraints. Note that the flight path angle is small and varies relatively slowly in glide flight. Thus, the QEGC can be constructed by settingcosθ=1 and θ̇=0 in Eq. (5) and ignoring Earth self-rotation as follows
(23)Lcosν+V2−1/r/r=0.As obviously seen from Eq. ((23)), if two arbitrary terms of the states in terms ofr,V, andνare given, another time-varying parameter could be determined along the glide trajectory. Based on this principle, the altitude versus velocity profile can be determined by choosing a suitable bank angle ν. Combining the exponential density equation (9) with path constraint equations (18–20), a collective altitude versus velocity profile corresponding to the three path constraints, which constitutes the lower boundary of the so-called entry flight corridor, can be simply deduced and intuitively represented by lrmin,V. Obviously, rmin is the geocentric distance corresponding to the lower boundary of the entry corridor. Correspondingly, a velocity-dependent upper boundary of the bank-angle magnitude can be derived from the QEGC in Eq. (23) and denoted by νmax, that is,
(24)νmax=cos−11/rmin2−V2/rminLmaxrmin,V.On the other hand, the lower boundary of the bank-angle magnitude can be given intuitively asνmin=νEG, where νEG is a specified bank angle to enforce the equilibrium glide constraint as mentioned above. In this paper, νEG=0 is used to determine the lower boundary of the bank angle. It is advisable that an appropriate ν should be chosen within the admissible region specified by νmin and νmax to enforce all of the inequality path constraints. In other words, for any V in the glide phase where the QEGC is valid, the entry trajectory will stay inside the entry flight corridor if ν is chosen from the following simple box constraint
(25)νminV≤νV≤νmaxV.Note that the preceding equation and arguments are based on the QEGC which is only valid for the glide phase but not the initial descent phase. In [28], this issue has been addressed using a simple Newton-Secant method to solve the equation Fνdes‐max=Q̇−Q̇max=0, and as a result, a constant νdes‐max is derived as the upper bound of ν for the initial descent phase. In fact, this boundary is not activated in most cases, which will be described in the next section. In addition, a possible compensation term concerning the Earth self-rotation and the heating rate constraint is added to yield a modified QEGC so as to achieve higher accuracy for the upper boundary of the bank angle (cf. [28]). Unfortunately, only a compromise result will be achieved due to the supposition of r≈1. Thus, we still utilize the box constraint (25) and Eq. (24) to enforce the three inequality path constraints as well as other constraints expressed in the velocity-altitude space.


### 3.2. Longitudinal Subplanning and Guidance Algorithm

Taking into account the QEGC and the distinctive characteristics of entry flight mechanics, the algorithm tactically divides the longitudinal profiles into the two well-known phases: the initial descent phase and the quasi-equilibrium glide phase. In the initial phase, the dynamic pressure of the vehicle is inefficient for the aerodynamic lift to shape the trajectory; hence, the bank angle is forced to be a constant. The quasi-equilibrium glide (QEG) phase, which is distinctive and unique for HGVs with moderate to higher L/D ratios, starts from a transition point in which the rationality of QEGC is insured, covers the majority of the entry trajectory, and plays a crucial role in satisfying all path constraints and other terminal conditions.


#### 3.2.1. Initial Descent Phase

The objective of guidance for the initial descent phase is to determine the trajectory state and the corresponding control command which steers the vehicle flying from the entry interface to a transition point connecting to the quasi-equilibrium glide phase. An effective algorithm has been proposed in many literatures with slight differences [15, 16], [18]. For completeness, we briefly describe here how the algorithm can be adopted to the initial descent planning problem. The magnitude of the feasible constant bank angle, i.e., νdes, is determined by increasing the bank angle from zero at a fixed incremental (the sign is given by the lateral guidance in the later section) and numerically integrating the equations of motion until the following criteria are simultaneously satisfied
(26)dr/dV−dr/dVQEGC≤δ,(27)Q̇≤Q̇max,where δ is a small preselected positive value. The preceding criteria indicate that at the intersecting point inside the entry flight corridor, the slopes of the descent trajectory and the quasi-equilibrium glide trajectory closely match.Dividing Eq. (1) with Eq. (4), and ignoring Earth self-rotation, we can obtain the slope of the descent trajectory at the current point r,V(28)drdV=−VsinθD+sinθ/r2,The other slopedr/dVQEGC is obtained by differentiating the QEGC once with respect to V at r,V(29)drdVQEGC=2/V1−V2r+2VrβR01−V2r+V2−2/r,where β=1/hs is a constant and the other variables are all dimensionless. Note that this way of determining dr/dVQEGC is more efficient than that of [28] because of no need for solving the QEGC. Finally, the integrated initial descent trajectory can be obtained once the appropriate νdes is determined. Also determined is the transition point, in which the states xtrans are afforded to be the initial conditions for the next trajectory generation.


#### 3.2.2. Quasi-Equilibrium Glide (QEG) Phase

In this phase, the magnitude of the bank angle command consists of two parts: the baseline part and the augmented part.(1) Range Control and Determination of the Baseline Bank Angle. The baseline command is derived from the analytical relation between the range-to-go and the velocity. As shown in the preceding section, we let stogo denote the range-to-go along the great circle connecting the current location of the vehicle and the final site on the surface of a spherical Earth. The time derivative for stogo is
(30)ṡtogo=−VcosθcosΔσ/r,where Δσ again is the offset between the heading angle and the azimuth of this great circle. Under the great circle assumption, the offset is so small that the usual approximation cosΔσ≈1 holds. Thus, Eq. (30) is simplified as
(31)ṡtogo=−Vcosθ/rDividingV̇ in equation (4) by ṡtogo and ignoring Earth self-rotation, we get the differential equation as follows
(32)dVdstogo=rD+sinθ/rVcosθ.Note thatθ≈0 and cosθ≈1 are acceptable when the QEGC is valid. Thus, in the QEG phase, Eq. (32) can be simplified to
(33)dVdstogo=rDV.ReplacingD with LCD/CL and substituting L from the QEGC in Eq.(23) lead to
(34)dVdstogo=1/r−V2CD/CLVcosν,which can be further rewritten as
(35)dstogo=CLCDcosν⋅V1/r−V2dV.Note that the dimensionless radial distancer varies so slowly in the QEG phase that it can be approximated as a constant value r~=rtrans+rf/2 (i.e., the average radial distance of the QEG phase). Since the angle of attack α is preselected to maintain the gliding flight with the highest L/D ratio, the relational term CL/CD could be assumed to be a constant too. We also assume that the bank angle ν is independent of the flight velocity V. Based on all the above assumptions, both sides of Eq. (35) can be analytically integrated into the corresponding interval for the range-to-go and the velocity, respectively. For the current state of arbitrary point, the integrated form of the preceding equation can be expressed as
(36)sf−stogo=−12CLCDcosν⋅lnV2−1/r~VVf,where V is the current velocity. The current range-to-go stogo is computed by spherical trigonometric functions and is given below to get the analytic solution
(37)stogo=cos−1sinϕfsinϕ+cosϕfcosϕcosλf−λ.Hence, it is easy to have(38)cosν=2CDCLstogo−sflnV2−1/r~VVf.Finally, the magnitude of the bank angle can be analytically obtained from Eq. (38). In fact, the assumptions of the constant r~ and CL/CD may not be sufficiently accurate. Also, the precision of stogo along the great circle is insufficient. Thus, all approximations are only used for the above analytic derivations in each guidance cycle, but continuously updated along the trajectory propagation. The accumulated errors of the final states can be reduced by the above successively updates. And it can be ensured by some augmented terms given in the following sections.(2) Command Augmentation and Altitude Control. Similar to the concepts and principles in [29], let νbase denote the analytically calculated bank angle at the current time by range control from Eq. (38). Note that some additional needs for trajectory shaping can be accomplished by augmenting the bank command νcmd by an altitude-rate compensation, expressed as [29].
(39)Lcosνcmd=Lcosνbase−kḣ−ḣref,where Lcosνbase is the vertical component of the baseline aerodynamic lift acceleration, ḣ is the current altitude rate, ḣref is the corresponding reference value with different forms for different purposes, and k>0 is a gain.In the rest of this section, a suitableḣref is designed and deduced analytically to eliminate the terminal altitude error. Consider the relation between the range-to-go and the radial distance. Dividing ṙ in equation (1) with ṡtogo, we get the variational equation as follows
(40)dstogo=−1tanθdrr.Because the fight path angle is very small in the QEG phase when the QEGC is valid, the approximationtanθ≈θ is acceptable.Similarly, for the current state of arbitrary point, the integrated form of the preceding equation can be expressed as(41)sf−stogo=−1θ⋅lnrrrf.Equation (41) gives rise to the required flight path angle to ensure the relation between the range-to-go and the radial distance
(42)θalt=1stogo−sf⋅lnrrrf,where subscript “alt” depicts the value accounting only for the terminal altitude constraint.Define the altitude rate required to altitude control by(43)ḣalt=Vsinθalt,where V is the current velocity. Substituting Eq. (43) into Eq. (39) gives
(44)Lcosνcmd‐alt=Lcosνbase−kaltḣ−ḣalt,where ḣ is the current altitude rate and νcmd‐alt is the commanded bank-angle magnitude required by the altitude control. The constant gain kalt>0 can be determined by simulations. Now, the magnitude of the commanded bank angle is calculated from Eq. (44), in which the altitude control is considered to reduce the terminal altitude error.(3) Command Augmentation Based on Velocity Control. Except for the need of altitude control, velocity control is also a crucial issue to be addressed so that the error of terminal velocity is tolerable. To meet the desired velocity at the final altitude (i.e., the final radial distance), the terminal velocity can be separated into two parts: one part due to atmospheric drag and the other from the gravity.Now, we will first determine the loss of velocity due to the aerodynamic drag. By ignoring the gravity term and Earth self-rotation, Eq. (4) can be rewritten as
(45)V̇=−D.Dividing the above equation byṙ in equation (1) and substituting D from Eq. (7) and ρ from Eq. (9) yields
(46)dVV=−kBρ0⋅1sinθ⋅R0exp−R0r−1/hs⋅dr,where kB is the ballistic coefficient with the form of
(47)kB=SrefCD2m.By treatingθ and CD all as constants, the analytical integration of Eq. (46) from the current state to the terminal state yields
(48)lnVVVDf=hskBsinθρ0exp−R0r−1/hsrrf.This result gives rise to the predicted terminal velocity accounting only for aerodynamic drag(49)VDf=exphskBsinθρ0exp−R0r−1/hsrrf=VexphskBsinθρf−ρ,where the subscript “D” denotes the effect of the aerodynamic drag and V is the current velocity. Hence, the loss of velocity due to the aerodynamic drag only is
(50)ΔVaero=VDf−V.Moreover, taking into account the orbital dynamics of the vehicle under the Earth gravitational field, the terminal velocity can easily be obtained using Keplerian laws(51)Vgf=V2+21/rf−1/r,where the subscript “g” denotes the effect of Earth gravity. Based on the law of conservation of energy, ΔVaero can also be expressed as
(52)ΔVaero=Vf−Vgf.Note that the true terminal velocity should be equal to the desired value, that is, Eq. (12) should be expected. To this done, the desired flight path angle is determined by substituting Eq. (49) into Eq. (50) and is expressed as follows
(53)sinθvel=hskBlnΔVaero/V+1ρ0exp−R0r−1/hsrrf=hskBlnΔVaero/V+1ρf−ρ,where ΔVaero is computed using Eq. (52) and (51) and the subscript “vel” denotes the derived value considering only the terminal velocity constraints. θvel is used to attain the desired velocity loss given by Eq. (52); hence, the specified terminal velocity can be reached when the flight terminates at the required final altitude.Similarly, define the altitude rate required only to the velocity control by(54)ḣvel=Vsinθvel.Substituting Eq. (54) into Eq. (39) gives
(55)Lcosνcmd‐vel=Lcosνbase−kvelḣ−ḣvel,where the constant gain kvel>0 can be determined by simulations and νcmd‐vel is the commanded bank-angle magnitude required by the terminal velocity constraints. Now, the magnitude of the commanded bank angle is calculated from Eq. (55), in which the velocity control is considered to reduce the terminal velocity error.Finally, accounting for the terminal constraints in terms of the range, altitude, and velocity depicted by Eqs. (11)–(13), the commanded vertical component of aerodynamic lift acceleration Lcosνcmd is taken as a weighted combination of that obtained from Eq. (44) and (55) as follows
(56)Lcosνcmd=ϖLcosνcmd‐alt+1−ϖLcosνcmd‐vel,where ϖ is a weighted value and ϖ∈0,1. Substituting Eq. (44) and (55), the above equation can be rewritten as
(57)Lcosνcmd=Lcosνbase−ϖkaltḣ−ḣalt−1−ϖkvelḣ−ḣvel.Setting new feedback gainsKalt=ϖkalt≥0 and Kvel=1−ϖkvel≥0, which should be scheduled synthetically by simulations, we get
(58)Lcosνcmd=Lcosνbase−Kaltḣ−ḣalt−Kvelḣ−ḣvel.It is worth noting thatνcmd should be limited by Eq. (25) to observe the inequality path constraints. That is,
(59)νcmd=νminV,ifνcmdV<νminVνcmdV,ifνminV≤νcmdV≤νminVνmaxV,ifνcmdV>νmaxV.Up to now, the magnitude of the commanded bank angle is ultimately calculated from Eq. (58), in which the terminal constraints in terms of the range, altitude, and velocity are all accounted for. Therefore, the remained task is to determine the sign of the bank angle that is presented in the next section.


### 3.3. Lateral Guidance Algorithm

With the longitudinal subplanning and guidance accomplished in the preceding section, we proceed to the lateral guidance problem to specify the sign of the bank angleν, so that the terminal heading error and crossrange are nullified or kept within specified tolerances, respectively.In the initial descent phase, the sign ofνdes is chosen to be opposite from that of the heading error Δψ. As mentioned previously, Δψ denotes the difference between the velocity azimuth angle and the line-of-sight angle from the vehicle to the target point and is expressed as
(60)Δψ=σ−ψLOS,where the line-of-sight to the final destination can be computed using spherical trigonometric functions as follows
(61)ψLOS=sin−1sinλf−λcosλf/sinstogo.Hence, the sign ofνdes is given by
(62)signνdes=−signΔψ0=−signσ0−ψLOS0,where the subscript “0” denotes the initial value of the trajectory parameters similar to the preceding section.In the QEG phase, the sign of bank angleνcmd is determined using an automatic, simple but efficient bank reversal logic to be discussed later in this section. Different from the lateral logic used by the Apollo and the Shuttle, we define a crossrange parameterχin radian instead of the heading error by
(63)χ=sin−1sinstogosinΔψ,which denotes the angle between the line-of-sight vector and its projection on the current flight plane. As demonstrated in [22], the appealing features of the crossrange parameter are the approximate piecewise linearity and slow variation with respect to the range-to-go for different vehicles/missions, which the heading error lacks in contrast. In this paper, the lateral logic presented automatically regulates the crossrange and corrects the heading error using an overcorrect scheme based on the dramatic feature as mentioned above. In the rest of this section, some parts of the principle in [22] are extended and revised to determine the reversal moment of the bank angle for lateral guidance.Considering the approximate linearity ofχ, we differentiate Eq. (63) with respect to stogo, and express the slope of dχ/dstogo at the current point as
(64)χ′=dχdstogo=cosstogosinΔψ+Δψ′sinstogocosΔψcosχ,where the prime denotes the derivative with respect to stogo, and similarly, we get
(65)Δψ′=−rV2cosθcosΔψLsinνcosθ+V2rcosθsinσtanϕ−ψLOS′.Assumingstogo≪1, we can also obtain ψLOS′ by the simplified expression
(66)ψLOS′=tanΔψ/stogo.Obviously, the value ofχ′ is a function of ν. It changes whenever the sign of ν is reversed as well and is depicted as χ′−ν. Furthermore, the bank angle criterion is depicted as follows.Suppose that the bank reversal to be immediately performed is the last one in the entry flight, and the reversal takes place at pointR. Let stogoR denote the range-to-go at point R. As mentioned above, we assume that the crossrange were truly linear with respect to stogo as long as stogo<stogoR. To ensure that the terminal heading error and crossrange are ideally nullified at the terminal range sf, the following relationship should be satisfied
(67)χR=χR′−νstogoR−sf,where the subscript “R” denotes the value at point R. The geometric meaning of χR and Eq. (67) is shown in Figure 1. It should be noted that R is not a fixed point. It can be clearly seen that the bank reversal occurs once the current χ satisfies Eq. (67), which leads to zero crossrange at sf. That is to say, as long as Eq. (67) is true at any point along the trajectory, the sign of ν should be reversed. If not, the crossrange error will exceed the specified tolerance. For example, if the reversal took place at point R2 when χR2>χR, this undercorrected reversal will be too late to meet the terminal constraints. That is, χsf>0 as seen in Figure 1. Conversely, the bank reversal at point R1 is overcorrected.Figure 1
Geometric illustration of the principle for the bank reversal.Unfortunately, the crossrange is not exactly linear, and the control constraints could cause some crossrange error. Hence, the bank reversal should take place no later than the instant when Eq. (67) is true. To do so, a margin is added multiplying the right side of Eq. (67) by a coefficient ε∈0,1. Therefore, ignoring the subscript “R,” the bank reversal will be performed when the following criterion is violated
(68)χ≤χthreshold=εχ′−νstogo−sf.Obviously, the above more conservative criterion could command the bank reversal so early that additional reversals may be needed later, which is preferred to improve the lateral precision. In essence, the dynamicχthreshold plays the role of a range-dependent threshold similar to that of heading error in the Apollo lateral guidance. The smaller the ε is, the tighter the χthreshold is. A tighterχthresholdleads to the bank reversal so early that excessive bank reversals should be performed later. An appropriate ε could strike a balance and hold a favorable precision without using excessively many bank reversals.For each guidance cycle in the QEG phase, we determine the magnitude of the commanded bank angle using the longitudinal guidance algorithm. Meanwhile, the sign of the bank angle is given by the lateral guidance algorithm. Then, Eqs. (1)–(6) are integrated using the commanded bank angle and the angle of attack in each guidance cycle. Once the above steps are accomplished, a feasible entry trajectory is generated as well as the closed-loop control commands. As seen from the guidance steps, any planning trajectory obtained is perfectly flyable.


## 3.1. QEGC and Translation of Inequality Path Constraints

As mentioned above, an ingenious utilization of QEGC is the cornerstone for this entry guidance algorithm and translation of inequality path constraints. Note that the flight path angle is small and varies relatively slowly in glide flight. Thus, the QEGC can be constructed by settingcosθ=1 and θ̇=0 in Eq. (5) and ignoring Earth self-rotation as follows
(23)Lcosν+V2−1/r/r=0.As obviously seen from Eq. ((23)), if two arbitrary terms of the states in terms ofr,V, andνare given, another time-varying parameter could be determined along the glide trajectory. Based on this principle, the altitude versus velocity profile can be determined by choosing a suitable bank angle ν. Combining the exponential density equation (9) with path constraint equations (18–20), a collective altitude versus velocity profile corresponding to the three path constraints, which constitutes the lower boundary of the so-called entry flight corridor, can be simply deduced and intuitively represented by lrmin,V. Obviously, rmin is the geocentric distance corresponding to the lower boundary of the entry corridor. Correspondingly, a velocity-dependent upper boundary of the bank-angle magnitude can be derived from the QEGC in Eq. (23) and denoted by νmax, that is,
(24)νmax=cos−11/rmin2−V2/rminLmaxrmin,V.On the other hand, the lower boundary of the bank-angle magnitude can be given intuitively asνmin=νEG, where νEG is a specified bank angle to enforce the equilibrium glide constraint as mentioned above. In this paper, νEG=0 is used to determine the lower boundary of the bank angle. It is advisable that an appropriate ν should be chosen within the admissible region specified by νmin and νmax to enforce all of the inequality path constraints. In other words, for any V in the glide phase where the QEGC is valid, the entry trajectory will stay inside the entry flight corridor if ν is chosen from the following simple box constraint
(25)νminV≤νV≤νmaxV.Note that the preceding equation and arguments are based on the QEGC which is only valid for the glide phase but not the initial descent phase. In [28], this issue has been addressed using a simple Newton-Secant method to solve the equation Fνdes‐max=Q̇−Q̇max=0, and as a result, a constant νdes‐max is derived as the upper bound of ν for the initial descent phase. In fact, this boundary is not activated in most cases, which will be described in the next section. In addition, a possible compensation term concerning the Earth self-rotation and the heating rate constraint is added to yield a modified QEGC so as to achieve higher accuracy for the upper boundary of the bank angle (cf. [28]). Unfortunately, only a compromise result will be achieved due to the supposition of r≈1. Thus, we still utilize the box constraint (25) and Eq. (24) to enforce the three inequality path constraints as well as other constraints expressed in the velocity-altitude space.


## 3.2. Longitudinal Subplanning and Guidance Algorithm

Taking into account the QEGC and the distinctive characteristics of entry flight mechanics, the algorithm tactically divides the longitudinal profiles into the two well-known phases: the initial descent phase and the quasi-equilibrium glide phase. In the initial phase, the dynamic pressure of the vehicle is inefficient for the aerodynamic lift to shape the trajectory; hence, the bank angle is forced to be a constant. The quasi-equilibrium glide (QEG) phase, which is distinctive and unique for HGVs with moderate to higher L/D ratios, starts from a transition point in which the rationality of QEGC is insured, covers the majority of the entry trajectory, and plays a crucial role in satisfying all path constraints and other terminal conditions.


### 3.2.1. Initial Descent Phase

The objective of guidance for the initial descent phase is to determine the trajectory state and the corresponding control command which steers the vehicle flying from the entry interface to a transition point connecting to the quasi-equilibrium glide phase. An effective algorithm has been proposed in many literatures with slight differences [15, 16], [18]. For completeness, we briefly describe here how the algorithm can be adopted to the initial descent planning problem. The magnitude of the feasible constant bank angle, i.e., νdes, is determined by increasing the bank angle from zero at a fixed incremental (the sign is given by the lateral guidance in the later section) and numerically integrating the equations of motion until the following criteria are simultaneously satisfied
(26)dr/dV−dr/dVQEGC≤δ,(27)Q̇≤Q̇max,where δ is a small preselected positive value. The preceding criteria indicate that at the intersecting point inside the entry flight corridor, the slopes of the descent trajectory and the quasi-equilibrium glide trajectory closely match.Dividing Eq. (1) with Eq. (4), and ignoring Earth self-rotation, we can obtain the slope of the descent trajectory at the current point r,V(28)drdV=−VsinθD+sinθ/r2,The other slopedr/dVQEGC is obtained by differentiating the QEGC once with respect to V at r,V(29)drdVQEGC=2/V1−V2r+2VrβR01−V2r+V2−2/r,where β=1/hs is a constant and the other variables are all dimensionless. Note that this way of determining dr/dVQEGC is more efficient than that of [28] because of no need for solving the QEGC. Finally, the integrated initial descent trajectory can be obtained once the appropriate νdes is determined. Also determined is the transition point, in which the states xtrans are afforded to be the initial conditions for the next trajectory generation.


### 3.2.2. Quasi-Equilibrium Glide (QEG) Phase

In this phase, the magnitude of the bank angle command consists of two parts: the baseline part and the augmented part.(1) Range Control and Determination of the Baseline Bank Angle. The baseline command is derived from the analytical relation between the range-to-go and the velocity. As shown in the preceding section, we let stogo denote the range-to-go along the great circle connecting the current location of the vehicle and the final site on the surface of a spherical Earth. The time derivative for stogo is
(30)ṡtogo=−VcosθcosΔσ/r,where Δσ again is the offset between the heading angle and the azimuth of this great circle. Under the great circle assumption, the offset is so small that the usual approximation cosΔσ≈1 holds. Thus, Eq. (30) is simplified as
(31)ṡtogo=−Vcosθ/rDividingV̇ in equation (4) by ṡtogo and ignoring Earth self-rotation, we get the differential equation as follows
(32)dVdstogo=rD+sinθ/rVcosθ.Note thatθ≈0 and cosθ≈1 are acceptable when the QEGC is valid. Thus, in the QEG phase, Eq. (32) can be simplified to
(33)dVdstogo=rDV.ReplacingD with LCD/CL and substituting L from the QEGC in Eq.(23) lead to
(34)dVdstogo=1/r−V2CD/CLVcosν,which can be further rewritten as
(35)dstogo=CLCDcosν⋅V1/r−V2dV.Note that the dimensionless radial distancer varies so slowly in the QEG phase that it can be approximated as a constant value r~=rtrans+rf/2 (i.e., the average radial distance of the QEG phase). Since the angle of attack α is preselected to maintain the gliding flight with the highest L/D ratio, the relational term CL/CD could be assumed to be a constant too. We also assume that the bank angle ν is independent of the flight velocity V. Based on all the above assumptions, both sides of Eq. (35) can be analytically integrated into the corresponding interval for the range-to-go and the velocity, respectively. For the current state of arbitrary point, the integrated form of the preceding equation can be expressed as
(36)sf−stogo=−12CLCDcosν⋅lnV2−1/r~VVf,where V is the current velocity. The current range-to-go stogo is computed by spherical trigonometric functions and is given below to get the analytic solution
(37)stogo=cos−1sinϕfsinϕ+cosϕfcosϕcosλf−λ.Hence, it is easy to have(38)cosν=2CDCLstogo−sflnV2−1/r~VVf.Finally, the magnitude of the bank angle can be analytically obtained from Eq. (38). In fact, the assumptions of the constant r~ and CL/CD may not be sufficiently accurate. Also, the precision of stogo along the great circle is insufficient. Thus, all approximations are only used for the above analytic derivations in each guidance cycle, but continuously updated along the trajectory propagation. The accumulated errors of the final states can be reduced by the above successively updates. And it can be ensured by some augmented terms given in the following sections.(2) Command Augmentation and Altitude Control. Similar to the concepts and principles in [29], let νbase denote the analytically calculated bank angle at the current time by range control from Eq. (38). Note that some additional needs for trajectory shaping can be accomplished by augmenting the bank command νcmd by an altitude-rate compensation, expressed as [29].
(39)Lcosνcmd=Lcosνbase−kḣ−ḣref,where Lcosνbase is the vertical component of the baseline aerodynamic lift acceleration, ḣ is the current altitude rate, ḣref is the corresponding reference value with different forms for different purposes, and k>0 is a gain.In the rest of this section, a suitableḣref is designed and deduced analytically to eliminate the terminal altitude error. Consider the relation between the range-to-go and the radial distance. Dividing ṙ in equation (1) with ṡtogo, we get the variational equation as follows
(40)dstogo=−1tanθdrr.Because the fight path angle is very small in the QEG phase when the QEGC is valid, the approximationtanθ≈θ is acceptable.Similarly, for the current state of arbitrary point, the integrated form of the preceding equation can be expressed as(41)sf−stogo=−1θ⋅lnrrrf.Equation (41) gives rise to the required flight path angle to ensure the relation between the range-to-go and the radial distance
(42)θalt=1stogo−sf⋅lnrrrf,where subscript “alt” depicts the value accounting only for the terminal altitude constraint.Define the altitude rate required to altitude control by(43)ḣalt=Vsinθalt,where V is the current velocity. Substituting Eq. (43) into Eq. (39) gives
(44)Lcosνcmd‐alt=Lcosνbase−kaltḣ−ḣalt,where ḣ is the current altitude rate and νcmd‐alt is the commanded bank-angle magnitude required by the altitude control. The constant gain kalt>0 can be determined by simulations. Now, the magnitude of the commanded bank angle is calculated from Eq. (44), in which the altitude control is considered to reduce the terminal altitude error.(3) Command Augmentation Based on Velocity Control. Except for the need of altitude control, velocity control is also a crucial issue to be addressed so that the error of terminal velocity is tolerable. To meet the desired velocity at the final altitude (i.e., the final radial distance), the terminal velocity can be separated into two parts: one part due to atmospheric drag and the other from the gravity.Now, we will first determine the loss of velocity due to the aerodynamic drag. By ignoring the gravity term and Earth self-rotation, Eq. (4) can be rewritten as
(45)V̇=−D.Dividing the above equation byṙ in equation (1) and substituting D from Eq. (7) and ρ from Eq. (9) yields
(46)dVV=−kBρ0⋅1sinθ⋅R0exp−R0r−1/hs⋅dr,where kB is the ballistic coefficient with the form of
(47)kB=SrefCD2m.By treatingθ and CD all as constants, the analytical integration of Eq. (46) from the current state to the terminal state yields
(48)lnVVVDf=hskBsinθρ0exp−R0r−1/hsrrf.This result gives rise to the predicted terminal velocity accounting only for aerodynamic drag(49)VDf=exphskBsinθρ0exp−R0r−1/hsrrf=VexphskBsinθρf−ρ,where the subscript “D” denotes the effect of the aerodynamic drag and V is the current velocity. Hence, the loss of velocity due to the aerodynamic drag only is
(50)ΔVaero=VDf−V.Moreover, taking into account the orbital dynamics of the vehicle under the Earth gravitational field, the terminal velocity can easily be obtained using Keplerian laws(51)Vgf=V2+21/rf−1/r,where the subscript “g” denotes the effect of Earth gravity. Based on the law of conservation of energy, ΔVaero can also be expressed as
(52)ΔVaero=Vf−Vgf.Note that the true terminal velocity should be equal to the desired value, that is, Eq. (12) should be expected. To this done, the desired flight path angle is determined by substituting Eq. (49) into Eq. (50) and is expressed as follows
(53)sinθvel=hskBlnΔVaero/V+1ρ0exp−R0r−1/hsrrf=hskBlnΔVaero/V+1ρf−ρ,where ΔVaero is computed using Eq. (52) and (51) and the subscript “vel” denotes the derived value considering only the terminal velocity constraints. θvel is used to attain the desired velocity loss given by Eq. (52); hence, the specified terminal velocity can be reached when the flight terminates at the required final altitude.Similarly, define the altitude rate required only to the velocity control by(54)ḣvel=Vsinθvel.Substituting Eq. (54) into Eq. (39) gives
(55)Lcosνcmd‐vel=Lcosνbase−kvelḣ−ḣvel,where the constant gain kvel>0 can be determined by simulations and νcmd‐vel is the commanded bank-angle magnitude required by the terminal velocity constraints. Now, the magnitude of the commanded bank angle is calculated from Eq. (55), in which the velocity control is considered to reduce the terminal velocity error.Finally, accounting for the terminal constraints in terms of the range, altitude, and velocity depicted by Eqs. (11)–(13), the commanded vertical component of aerodynamic lift acceleration Lcosνcmd is taken as a weighted combination of that obtained from Eq. (44) and (55) as follows
(56)Lcosνcmd=ϖLcosνcmd‐alt+1−ϖLcosνcmd‐vel,where ϖ is a weighted value and ϖ∈0,1. Substituting Eq. (44) and (55), the above equation can be rewritten as
(57)Lcosνcmd=Lcosνbase−ϖkaltḣ−ḣalt−1−ϖkvelḣ−ḣvel.Setting new feedback gainsKalt=ϖkalt≥0 and Kvel=1−ϖkvel≥0, which should be scheduled synthetically by simulations, we get
(58)Lcosνcmd=Lcosνbase−Kaltḣ−ḣalt−Kvelḣ−ḣvel.It is worth noting thatνcmd should be limited by Eq. (25) to observe the inequality path constraints. That is,
(59)νcmd=νminV,ifνcmdV<νminVνcmdV,ifνminV≤νcmdV≤νminVνmaxV,ifνcmdV>νmaxV.Up to now, the magnitude of the commanded bank angle is ultimately calculated from Eq. (58), in which the terminal constraints in terms of the range, altitude, and velocity are all accounted for. Therefore, the remained task is to determine the sign of the bank angle that is presented in the next section.


## 3.2.1. Initial Descent Phase

The objective of guidance for the initial descent phase is to determine the trajectory state and the corresponding control command which steers the vehicle flying from the entry interface to a transition point connecting to the quasi-equilibrium glide phase. An effective algorithm has been proposed in many literatures with slight differences [15, 16], [18]. For completeness, we briefly describe here how the algorithm can be adopted to the initial descent planning problem. The magnitude of the feasible constant bank angle, i.e., νdes, is determined by increasing the bank angle from zero at a fixed incremental (the sign is given by the lateral guidance in the later section) and numerically integrating the equations of motion until the following criteria are simultaneously satisfied
(26)dr/dV−dr/dVQEGC≤δ,(27)Q̇≤Q̇max,where δ is a small preselected positive value. The preceding criteria indicate that at the intersecting point inside the entry flight corridor, the slopes of the descent trajectory and the quasi-equilibrium glide trajectory closely match.Dividing Eq. (1) with Eq. (4), and ignoring Earth self-rotation, we can obtain the slope of the descent trajectory at the current point r,V(28)drdV=−VsinθD+sinθ/r2,The other slopedr/dVQEGC is obtained by differentiating the QEGC once with respect to V at r,V(29)drdVQEGC=2/V1−V2r+2VrβR01−V2r+V2−2/r,where β=1/hs is a constant and the other variables are all dimensionless. Note that this way of determining dr/dVQEGC is more efficient than that of [28] because of no need for solving the QEGC. Finally, the integrated initial descent trajectory can be obtained once the appropriate νdes is determined. Also determined is the transition point, in which the states xtrans are afforded to be the initial conditions for the next trajectory generation.


## 3.2.2. Quasi-Equilibrium Glide (QEG) Phase

In this phase, the magnitude of the bank angle command consists of two parts: the baseline part and the augmented part.(1) Range Control and Determination of the Baseline Bank Angle. The baseline command is derived from the analytical relation between the range-to-go and the velocity. As shown in the preceding section, we let stogo denote the range-to-go along the great circle connecting the current location of the vehicle and the final site on the surface of a spherical Earth. The time derivative for stogo is
(30)ṡtogo=−VcosθcosΔσ/r,where Δσ again is the offset between the heading angle and the azimuth of this great circle. Under the great circle assumption, the offset is so small that the usual approximation cosΔσ≈1 holds. Thus, Eq. (30) is simplified as
(31)ṡtogo=−Vcosθ/rDividingV̇ in equation (4) by ṡtogo and ignoring Earth self-rotation, we get the differential equation as follows
(32)dVdstogo=rD+sinθ/rVcosθ.Note thatθ≈0 and cosθ≈1 are acceptable when the QEGC is valid. Thus, in the QEG phase, Eq. (32) can be simplified to
(33)dVdstogo=rDV.ReplacingD with LCD/CL and substituting L from the QEGC in Eq.(23) lead to
(34)dVdstogo=1/r−V2CD/CLVcosν,which can be further rewritten as
(35)dstogo=CLCDcosν⋅V1/r−V2dV.Note that the dimensionless radial distancer varies so slowly in the QEG phase that it can be approximated as a constant value r~=rtrans+rf/2 (i.e., the average radial distance of the QEG phase). Since the angle of attack α is preselected to maintain the gliding flight with the highest L/D ratio, the relational term CL/CD could be assumed to be a constant too. We also assume that the bank angle ν is independent of the flight velocity V. Based on all the above assumptions, both sides of Eq. (35) can be analytically integrated into the corresponding interval for the range-to-go and the velocity, respectively. For the current state of arbitrary point, the integrated form of the preceding equation can be expressed as
(36)sf−stogo=−12CLCDcosν⋅lnV2−1/r~VVf,where V is the current velocity. The current range-to-go stogo is computed by spherical trigonometric functions and is given below to get the analytic solution
(37)stogo=cos−1sinϕfsinϕ+cosϕfcosϕcosλf−λ.Hence, it is easy to have(38)cosν=2CDCLstogo−sflnV2−1/r~VVf.Finally, the magnitude of the bank angle can be analytically obtained from Eq. (38). In fact, the assumptions of the constant r~ and CL/CD may not be sufficiently accurate. Also, the precision of stogo along the great circle is insufficient. Thus, all approximations are only used for the above analytic derivations in each guidance cycle, but continuously updated along the trajectory propagation. The accumulated errors of the final states can be reduced by the above successively updates. And it can be ensured by some augmented terms given in the following sections.(2) Command Augmentation and Altitude Control. Similar to the concepts and principles in [29], let νbase denote the analytically calculated bank angle at the current time by range control from Eq. (38). Note that some additional needs for trajectory shaping can be accomplished by augmenting the bank command νcmd by an altitude-rate compensation, expressed as [29].
(39)Lcosνcmd=Lcosνbase−kḣ−ḣref,where Lcosνbase is the vertical component of the baseline aerodynamic lift acceleration, ḣ is the current altitude rate, ḣref is the corresponding reference value with different forms for different purposes, and k>0 is a gain.In the rest of this section, a suitableḣref is designed and deduced analytically to eliminate the terminal altitude error. Consider the relation between the range-to-go and the radial distance. Dividing ṙ in equation (1) with ṡtogo, we get the variational equation as follows
(40)dstogo=−1tanθdrr.Because the fight path angle is very small in the QEG phase when the QEGC is valid, the approximationtanθ≈θ is acceptable.Similarly, for the current state of arbitrary point, the integrated form of the preceding equation can be expressed as(41)sf−stogo=−1θ⋅lnrrrf.Equation (41) gives rise to the required flight path angle to ensure the relation between the range-to-go and the radial distance
(42)θalt=1stogo−sf⋅lnrrrf,where subscript “alt” depicts the value accounting only for the terminal altitude constraint.Define the altitude rate required to altitude control by(43)ḣalt=Vsinθalt,where V is the current velocity. Substituting Eq. (43) into Eq. (39) gives
(44)Lcosνcmd‐alt=Lcosνbase−kaltḣ−ḣalt,where ḣ is the current altitude rate and νcmd‐alt is the commanded bank-angle magnitude required by the altitude control. The constant gain kalt>0 can be determined by simulations. Now, the magnitude of the commanded bank angle is calculated from Eq. (44), in which the altitude control is considered to reduce the terminal altitude error.(3) Command Augmentation Based on Velocity Control. Except for the need of altitude control, velocity control is also a crucial issue to be addressed so that the error of terminal velocity is tolerable. To meet the desired velocity at the final altitude (i.e., the final radial distance), the terminal velocity can be separated into two parts: one part due to atmospheric drag and the other from the gravity.Now, we will first determine the loss of velocity due to the aerodynamic drag. By ignoring the gravity term and Earth self-rotation, Eq. (4) can be rewritten as
(45)V̇=−D.Dividing the above equation byṙ in equation (1) and substituting D from Eq. (7) and ρ from Eq. (9) yields
(46)dVV=−kBρ0⋅1sinθ⋅R0exp−R0r−1/hs⋅dr,where kB is the ballistic coefficient with the form of
(47)kB=SrefCD2m.By treatingθ and CD all as constants, the analytical integration of Eq. (46) from the current state to the terminal state yields
(48)lnVVVDf=hskBsinθρ0exp−R0r−1/hsrrf.This result gives rise to the predicted terminal velocity accounting only for aerodynamic drag(49)VDf=exphskBsinθρ0exp−R0r−1/hsrrf=VexphskBsinθρf−ρ,where the subscript “D” denotes the effect of the aerodynamic drag and V is the current velocity. Hence, the loss of velocity due to the aerodynamic drag only is
(50)ΔVaero=VDf−V.Moreover, taking into account the orbital dynamics of the vehicle under the Earth gravitational field, the terminal velocity can easily be obtained using Keplerian laws(51)Vgf=V2+21/rf−1/r,where the subscript “g” denotes the effect of Earth gravity. Based on the law of conservation of energy, ΔVaero can also be expressed as
(52)ΔVaero=Vf−Vgf.Note that the true terminal velocity should be equal to the desired value, that is, Eq. (12) should be expected. To this done, the desired flight path angle is determined by substituting Eq. (49) into Eq. (50) and is expressed as follows
(53)sinθvel=hskBlnΔVaero/V+1ρ0exp−R0r−1/hsrrf=hskBlnΔVaero/V+1ρf−ρ,where ΔVaero is computed using Eq. (52) and (51) and the subscript “vel” denotes the derived value considering only the terminal velocity constraints. θvel is used to attain the desired velocity loss given by Eq. (52); hence, the specified terminal velocity can be reached when the flight terminates at the required final altitude.Similarly, define the altitude rate required only to the velocity control by(54)ḣvel=Vsinθvel.Substituting Eq. (54) into Eq. (39) gives
(55)Lcosνcmd‐vel=Lcosνbase−kvelḣ−ḣvel,where the constant gain kvel>0 can be determined by simulations and νcmd‐vel is the commanded bank-angle magnitude required by the terminal velocity constraints. Now, the magnitude of the commanded bank angle is calculated from Eq. (55), in which the velocity control is considered to reduce the terminal velocity error.Finally, accounting for the terminal constraints in terms of the range, altitude, and velocity depicted by Eqs. (11)–(13), the commanded vertical component of aerodynamic lift acceleration Lcosνcmd is taken as a weighted combination of that obtained from Eq. (44) and (55) as follows
(56)Lcosνcmd=ϖLcosνcmd‐alt+1−ϖLcosνcmd‐vel,where ϖ is a weighted value and ϖ∈0,1. Substituting Eq. (44) and (55), the above equation can be rewritten as
(57)Lcosνcmd=Lcosνbase−ϖkaltḣ−ḣalt−1−ϖkvelḣ−ḣvel.Setting new feedback gainsKalt=ϖkalt≥0 and Kvel=1−ϖkvel≥0, which should be scheduled synthetically by simulations, we get
(58)Lcosνcmd=Lcosνbase−Kaltḣ−ḣalt−Kvelḣ−ḣvel.It is worth noting thatνcmd should be limited by Eq. (25) to observe the inequality path constraints. That is,
(59)νcmd=νminV,ifνcmdV<νminVνcmdV,ifνminV≤νcmdV≤νminVνmaxV,ifνcmdV>νmaxV.Up to now, the magnitude of the commanded bank angle is ultimately calculated from Eq. (58), in which the terminal constraints in terms of the range, altitude, and velocity are all accounted for. Therefore, the remained task is to determine the sign of the bank angle that is presented in the next section.


## 3.3. Lateral Guidance Algorithm

With the longitudinal subplanning and guidance accomplished in the preceding section, we proceed to the lateral guidance problem to specify the sign of the bank angleν, so that the terminal heading error and crossrange are nullified or kept within specified tolerances, respectively.In the initial descent phase, the sign ofνdes is chosen to be opposite from that of the heading error Δψ. As mentioned previously, Δψ denotes the difference between the velocity azimuth angle and the line-of-sight angle from the vehicle to the target point and is expressed as
(60)Δψ=σ−ψLOS,where the line-of-sight to the final destination can be computed using spherical trigonometric functions as follows
(61)ψLOS=sin−1sinλf−λcosλf/sinstogo.Hence, the sign ofνdes is given by
(62)signνdes=−signΔψ0=−signσ0−ψLOS0,where the subscript “0” denotes the initial value of the trajectory parameters similar to the preceding section.In the QEG phase, the sign of bank angleνcmd is determined using an automatic, simple but efficient bank reversal logic to be discussed later in this section. Different from the lateral logic used by the Apollo and the Shuttle, we define a crossrange parameterχin radian instead of the heading error by
(63)χ=sin−1sinstogosinΔψ,which denotes the angle between the line-of-sight vector and its projection on the current flight plane. As demonstrated in [22], the appealing features of the crossrange parameter are the approximate piecewise linearity and slow variation with respect to the range-to-go for different vehicles/missions, which the heading error lacks in contrast. In this paper, the lateral logic presented automatically regulates the crossrange and corrects the heading error using an overcorrect scheme based on the dramatic feature as mentioned above. In the rest of this section, some parts of the principle in [22] are extended and revised to determine the reversal moment of the bank angle for lateral guidance.Considering the approximate linearity ofχ, we differentiate Eq. (63) with respect to stogo, and express the slope of dχ/dstogo at the current point as
(64)χ′=dχdstogo=cosstogosinΔψ+Δψ′sinstogocosΔψcosχ,where the prime denotes the derivative with respect to stogo, and similarly, we get
(65)Δψ′=−rV2cosθcosΔψLsinνcosθ+V2rcosθsinσtanϕ−ψLOS′.Assumingstogo≪1, we can also obtain ψLOS′ by the simplified expression
(66)ψLOS′=tanΔψ/stogo.Obviously, the value ofχ′ is a function of ν. It changes whenever the sign of ν is reversed as well and is depicted as χ′−ν. Furthermore, the bank angle criterion is depicted as follows.Suppose that the bank reversal to be immediately performed is the last one in the entry flight, and the reversal takes place at pointR. Let stogoR denote the range-to-go at point R. As mentioned above, we assume that the crossrange were truly linear with respect to stogo as long as stogo<stogoR. To ensure that the terminal heading error and crossrange are ideally nullified at the terminal range sf, the following relationship should be satisfied
(67)χR=χR′−νstogoR−sf,where the subscript “R” denotes the value at point R. The geometric meaning of χR and Eq. (67) is shown in Figure 1. It should be noted that R is not a fixed point. It can be clearly seen that the bank reversal occurs once the current χ satisfies Eq. (67), which leads to zero crossrange at sf. That is to say, as long as Eq. (67) is true at any point along the trajectory, the sign of ν should be reversed. If not, the crossrange error will exceed the specified tolerance. For example, if the reversal took place at point R2 when χR2>χR, this undercorrected reversal will be too late to meet the terminal constraints. That is, χsf>0 as seen in Figure 1. Conversely, the bank reversal at point R1 is overcorrected.Figure 1
Geometric illustration of the principle for the bank reversal.Unfortunately, the crossrange is not exactly linear, and the control constraints could cause some crossrange error. Hence, the bank reversal should take place no later than the instant when Eq. (67) is true. To do so, a margin is added multiplying the right side of Eq. (67) by a coefficient ε∈0,1. Therefore, ignoring the subscript “R,” the bank reversal will be performed when the following criterion is violated
(68)χ≤χthreshold=εχ′−νstogo−sf.Obviously, the above more conservative criterion could command the bank reversal so early that additional reversals may be needed later, which is preferred to improve the lateral precision. In essence, the dynamicχthreshold plays the role of a range-dependent threshold similar to that of heading error in the Apollo lateral guidance. The smaller the ε is, the tighter the χthreshold is. A tighterχthresholdleads to the bank reversal so early that excessive bank reversals should be performed later. An appropriate ε could strike a balance and hold a favorable precision without using excessively many bank reversals.For each guidance cycle in the QEG phase, we determine the magnitude of the commanded bank angle using the longitudinal guidance algorithm. Meanwhile, the sign of the bank angle is given by the lateral guidance algorithm. Then, Eqs. (1)–(6) are integrated using the commanded bank angle and the angle of attack in each guidance cycle. Once the above steps are accomplished, a feasible entry trajectory is generated as well as the closed-loop control commands. As seen from the guidance steps, any planning trajectory obtained is perfectly flyable.


## 4. Numerical Examples

### 4.1. Vehicle and Missions

The simulations presented in this section use the model of Lockheed-Martin’s CAV-H, which is a typical hypersonic gliding vehicle with a high lifting and lift-to-drag ratio. The CAV-H has a mass of 907 kg, with the reference area of 0.4839m2, and the maximum lift-to-drag ratio of 3.5 corresponding to an angle of attack of about 10∘. The nominal α profile is fixed and given by a piece of linear function of velocity
(69)α=20deg,V≥V110−20V2−V1V−V1+20deg,V2≤V<V110deg,V<V2,where V1=4800m/s and V2=2500m/s. The aerodynamic lift and drag coefficients are fitted by the functions of the angle of attack and Mach number using the tabulated data. The flight control authority is restricted by the following conditions: ν∈−80,80deg, ν̇≤20deg/s.To evaluate the guidance algorithm, several mission scenarios are set up and tested with the same initial conditions of entry interface and different terminal conditions for different flight missions. The uniform initial conditions areh0=80km, λ0=0deg, ϕ0=0deg, V0=7000m/s, θ0=0deg, and σ0=60deg. The terminal conditions for different missions are listed in Table 1, including the specified final velocity, longitude, latitude, and the peak heating rate limit. The range-to-go and crossrange at the entry interface were also computed and listed. Negative values indicate the left crossranges for lateral motions. The first mission is the nominal case. The other terminal conditions are hf=25km, sf=50km, and Δψf=5deg. In addition, the other same peak path constraints on entry trajectories for all cases are nmax=3 and qmax=150kPa.Table 1
Entry mission scenarios.
MissionVf (m/s)λf (deg)ϕf (deg)stogo0(km)χ0 (km)Q̇max(kW/m2)1220065257503-28515002220065257503-28512503240065257503-28515004200065257503-28515005220060106727-178212506220060407503144615007220070107819-20241250822007040831912621500922007040831912621250The guidance algorithm is coded and implemented in MATLAB on an ordinary laptop computer. The update rate of the commanded bank angle in all tests is set to be 1 Hz for the QEG phase. The integration step is set to be 0.1 s for the initial descent phase and 0.01 s for the QEG phase.


### 4.2. Simulation Results

#### 4.2.1. Preliminary Testing and Discussions

Firstly, the test for mission 1 was implemented to verify the principles of the guidance algorithm and assess how well the algorithm works. A conventional bank reversal logic used by the Apollo and the Shuttle is demonstrated and validated. As can be seen in [29], the sign of bank angle should be maintained until the following criterion is violated
(70)Δψ≤ΔazmthV.An appropriate design ofΔazmthV, which avoids too many bank reversals, depends on the vehicle performance and the mission scenarios and is time-consuming. For comparison, a velocity-dependent ΔazmthV with a piecewise linear form, which is similar to that in [29], is carefully chosen to the above lateral logic by trial simulations.The tests of mission 1 were completed for three cases in terms of the baseline algorithm described in Eq. (38) (represented by BA), the augmented algorithm described in Eq. (58) with new lateral logic in Eq. (68) (represented by AANLL), and the same augmented algorithm but with conventional lateral logic (represented by AACLL). The comparison of the testing results for the three cases is shown in Figure 2.Figure 2
Trajectory comparison of the three cases for mission 1.
(a)
Entry trajectories(b)
Ground tracks(c)
Flight path angles(d)
Bank angles(e)
Crossranges(f)
Heading errorsFigures2(a) and 2(b) show the altitude versus range-to-go profiles and the ground tracks. Although the baseline algorithm has poor precision and large phugoid oscillations in the altitude, the augmented algorithm performs well for two lateral logics. The large phugoid oscillations are mostly eliminated with a high terminal precision by the feedback augmentation, which can also be seen from the flight path angle plotted in Figure 2(c). The AANLL gives a terminal altitude error of 73.94 m, a velocity error of 1.54 m/s, and a range error of 121.93 m. Based on the carefully chosen ΔazmthV, the AACLL gives the same level of precision but different ground track as shown in Figure 2(b) because of the different bank reversals. Also seen in Figure 2(c) is that the flight path angles of the AANLL are small negative values and vary rather smoothly except for the initial descent phase. Hence, the corresponding hypothesis for the QEGC is reasonable and effective as well as an approximately equilibrium glide.Figure2(d) shows significantly the different bank angle histories. The bank angle is maintained to be a constant zero till the initial descent phase terminates. In the QEG phase, the AANLL gives four bank reversals, which is less than the AACLL and leads to different ground tracks as can be seen in Figure 2(b). Figure 2(e) and 3(f) give the comparison of heading errors and crossranges for the above three cases. With the feedback augmentation in Eq. (58), the corresponding algorithms render the approximately piecewise linearity crossrange parameter acceptable and reasonable, in contrast to the volatilization and the nonlinearity of heading errors. The final heading error is only about 0.002 deg for AANLL, dramatically better than that of 3.08 deg for AACLL. This validates the high lateral precision that the AANLL can offer. Because the first bank reversal of the AANLG occurs considerably later than that of the AACLL, the maximum of the crossrange parameter is larger than that of the AACLL.Figure 3
Crossranges and bank angles of the AANLL for differentε.
(a)
Crossranges(b)
Bank anglesFigure3 compares the crossranges and bank angles with different scaling factors ε for the AANLL. Obviously, as ε is increased, the maximal lateral excursion increases while the bank reversals decrease, aside from the rearward shift of the reversal point in time and loss of accuracy. Hence, an appropriate ε should be chosen to balance a preferred terminal precision and bank reversals.


#### 4.2.2. Adaptability Testing and Simulations

As a first step in testing and assessing the efficiency and adaptability of the guidance algorithm represented as AANLL, we present and discuss the results of all mission scenarios set up and listed in Table1. For demonstration and comparison, all initial conditions and guidance parameters were kept the same as mission 1 for all missions.The nominal terminal conditions for all the above mission scenarios are listed in Table2. Note that all missions have a rather high accuracy. The terminal altitude errors are all less than 1 km, the velocity errors are less than 5 m/s, and range errors are all less than 3 km. The heading errors are significantly less than the 5 deg requirement, the maximum of which is a rather small value of -0.52 deg. The peak heating rates for all missions are no more than the corresponding limits, respectively. In addition, the other two peak path constraints are also well satisfied. In fact, the peak aerodynamic loads and dynamic pressures for all missions did not exceed 115 kPa and 2.5, which are no more than the corresponding limits, respectively. Therefore, Table 2 only gives a focus on the comparison of peak heating rates.Table 2
Terminal condition precision for all missions.
MissionΔhf(m)ΔVf (m/s)Δsf (km)Δψf (deg)MaximumQ̇ (kW/m2)173.941.540.120.001294.582-268.470.010.180.261249.913991.72-7.522.880.031292.364-486.171.84-0.110.181295.705-666.201.982.42-0.521250.006-152.78-0.39-0.440.121293.217328.33-2.77-1.00-0.241249.788853.92-4.530.38-0.031278.489405.52-4.48-0.320.071249.64The computation time used for generating the entry trajectory for each mission with flight time of about 1500 s is only about 3-4 s. It should be noted that most time is consumed by integrating Eqs. (1)–(6) throughout the entire entry trajectory. In every guidance cycle, the computation time required to generate the commanded bank angle is dramatically less than 1 ms. Note that all computations were implemented on a laptop computer and the algorithm is coded in MATLAB without any optimization. Predictably, improvements in software and hardware could provide large room for improvement in computation speed. Thus, the guidance algorithm has an indubitable potential for onboard application.Taking missions 5, 6, 7, 8, and 9 as examples, Figure4 shows the comparison of altitude, ground track, bank angle, and heating rate histories for the QEG phase only. It is evident from Figure 4(a) that the large phugoid oscillations in the initial part are gradually eliminated along the trajectory propagations, which ensures the validity of the QEGC. The ground tracks shown in Figure 4(b) demonstrate the high position accuracy for different missions. The characteristic of phugoid oscillations is also verified by the flight path angles shown in Figure 2(c). It can be seen from Figure 2(d) that the entry trajectories reach the final site with different heading angle (i.e., in different directions) for different missions.Figure 4
Trajectory comparison for all missions.
(a)
Longitudinal profiles(b)
Ground tracks(c)
Flight path angles(d)
Heading errors(e)
Bank angles(f)
Heating ratesCombining with the heating rate profiles is shown in Figure4(f), and Figure 4(e) illustrates the validity of the translation of inequality path constraints. Look closely at the comparison of the bank angle histories for missions 8 and 9. One of the main differences lies in the initial bank angle when the QEG phase initiates. As expected, the magnitude of the bank angle changes considerably as the heating rate constraint is imposed on. Once the trajectory has entered the QEG phase, the magnitude of the bank angle of mission 8 will increase immediately to focus on the range requirement due to the loose heating rate limit. The path constraint is not active for this case. However, the magnitude of bank angle of mission 9 is still maintained to be zero for the initial small time range, so that the trajectory would be forced to drive shallower into the dense atmosphere initially. It is also confirmed by the altitude profiles as seen in Figure 4(a). Then, the peak heating rate, often found on the first trough point, is induced effectively as shown in Figure 4(f).To further test and assess the efficiency and adaptability of the guidance algorithm, 100 dispersion cases were studied for mission 1 using the Monte Carlo simulations. The dispersed initial conditions are considered and modeled by the zero-mean Gaussian dispersions with 3-sigma values. The 3-sigma value of initial entry condition dispersions are as follows: the dispersed altitude of 10 km, the longitude and latitude of 2 deg, the velocity of 300 m/s, the flight path angle of 0.2 deg, and the heading angle of 3 deg, respectively. Table3 summarizes the statistics on the final conditions for 100 dispersed trajectories. Obviously, the small final errors are all derived for each final state concerned in terms of small means and standard deviations. The efficiency, robustness, and adaptability are dramatically confirmed and approved again. However, in the presence of significant aerodynamic modeling uncertainties and atmosphere modeling dispersions (the atmosphere density, aerodynamic coefficients all have the 3-sigma values of 15% in respective nominal values), this algorithm performs not well with large scatters especially in the terminal range error and heading error. Hence, no results are given here. When the guidance algorithm serves as an entry guidance approach, some techniques (e.g., aerodynamics filter, compare [27]) can be used to improve the guidance performance. This comment will be carried out in future works about entry guidance.Table 3
Statistics on the final conditions for 100 dispersed trajectories.
hf (km)Vf (m/s)sf (km)Δψf (deg)Average25.062198.1949.670.12Max26.002202.6050.771.03Min24.272194.3146.50-0.60Standard deviation0.0401.770.830.23


## 4.1. Vehicle and Missions

The simulations presented in this section use the model of Lockheed-Martin’s CAV-H, which is a typical hypersonic gliding vehicle with a high lifting and lift-to-drag ratio. The CAV-H has a mass of 907 kg, with the reference area of 0.4839m2, and the maximum lift-to-drag ratio of 3.5 corresponding to an angle of attack of about 10∘. The nominal α profile is fixed and given by a piece of linear function of velocity
(69)α=20deg,V≥V110−20V2−V1V−V1+20deg,V2≤V<V110deg,V<V2,where V1=4800m/s and V2=2500m/s. The aerodynamic lift and drag coefficients are fitted by the functions of the angle of attack and Mach number using the tabulated data. The flight control authority is restricted by the following conditions: ν∈−80,80deg, ν̇≤20deg/s.To evaluate the guidance algorithm, several mission scenarios are set up and tested with the same initial conditions of entry interface and different terminal conditions for different flight missions. The uniform initial conditions areh0=80km, λ0=0deg, ϕ0=0deg, V0=7000m/s, θ0=0deg, and σ0=60deg. The terminal conditions for different missions are listed in Table 1, including the specified final velocity, longitude, latitude, and the peak heating rate limit. The range-to-go and crossrange at the entry interface were also computed and listed. Negative values indicate the left crossranges for lateral motions. The first mission is the nominal case. The other terminal conditions are hf=25km, sf=50km, and Δψf=5deg. In addition, the other same peak path constraints on entry trajectories for all cases are nmax=3 and qmax=150kPa.Table 1
Entry mission scenarios.
MissionVf (m/s)λf (deg)ϕf (deg)stogo0(km)χ0 (km)Q̇max(kW/m2)1220065257503-28515002220065257503-28512503240065257503-28515004200065257503-28515005220060106727-178212506220060407503144615007220070107819-20241250822007040831912621500922007040831912621250The guidance algorithm is coded and implemented in MATLAB on an ordinary laptop computer. The update rate of the commanded bank angle in all tests is set to be 1 Hz for the QEG phase. The integration step is set to be 0.1 s for the initial descent phase and 0.01 s for the QEG phase.


## 4.2. Simulation Results

### 4.2.1. Preliminary Testing and Discussions

Firstly, the test for mission 1 was implemented to verify the principles of the guidance algorithm and assess how well the algorithm works. A conventional bank reversal logic used by the Apollo and the Shuttle is demonstrated and validated. As can be seen in [29], the sign of bank angle should be maintained until the following criterion is violated
(70)Δψ≤ΔazmthV.An appropriate design ofΔazmthV, which avoids too many bank reversals, depends on the vehicle performance and the mission scenarios and is time-consuming. For comparison, a velocity-dependent ΔazmthV with a piecewise linear form, which is similar to that in [29], is carefully chosen to the above lateral logic by trial simulations.The tests of mission 1 were completed for three cases in terms of the baseline algorithm described in Eq. (38) (represented by BA), the augmented algorithm described in Eq. (58) with new lateral logic in Eq. (68) (represented by AANLL), and the same augmented algorithm but with conventional lateral logic (represented by AACLL). The comparison of the testing results for the three cases is shown in Figure 2.Figure 2
Trajectory comparison of the three cases for mission 1.
(a)
Entry trajectories(b)
Ground tracks(c)
Flight path angles(d)
Bank angles(e)
Crossranges(f)
Heading errorsFigures2(a) and 2(b) show the altitude versus range-to-go profiles and the ground tracks. Although the baseline algorithm has poor precision and large phugoid oscillations in the altitude, the augmented algorithm performs well for two lateral logics. The large phugoid oscillations are mostly eliminated with a high terminal precision by the feedback augmentation, which can also be seen from the flight path angle plotted in Figure 2(c). The AANLL gives a terminal altitude error of 73.94 m, a velocity error of 1.54 m/s, and a range error of 121.93 m. Based on the carefully chosen ΔazmthV, the AACLL gives the same level of precision but different ground track as shown in Figure 2(b) because of the different bank reversals. Also seen in Figure 2(c) is that the flight path angles of the AANLL are small negative values and vary rather smoothly except for the initial descent phase. Hence, the corresponding hypothesis for the QEGC is reasonable and effective as well as an approximately equilibrium glide.Figure2(d) shows significantly the different bank angle histories. The bank angle is maintained to be a constant zero till the initial descent phase terminates. In the QEG phase, the AANLL gives four bank reversals, which is less than the AACLL and leads to different ground tracks as can be seen in Figure 2(b). Figure 2(e) and 3(f) give the comparison of heading errors and crossranges for the above three cases. With the feedback augmentation in Eq. (58), the corresponding algorithms render the approximately piecewise linearity crossrange parameter acceptable and reasonable, in contrast to the volatilization and the nonlinearity of heading errors. The final heading error is only about 0.002 deg for AANLL, dramatically better than that of 3.08 deg for AACLL. This validates the high lateral precision that the AANLL can offer. Because the first bank reversal of the AANLG occurs considerably later than that of the AACLL, the maximum of the crossrange parameter is larger than that of the AACLL.Figure 3
Crossranges and bank angles of the AANLL for differentε.
(a)
Crossranges(b)
Bank anglesFigure3 compares the crossranges and bank angles with different scaling factors ε for the AANLL. Obviously, as ε is increased, the maximal lateral excursion increases while the bank reversals decrease, aside from the rearward shift of the reversal point in time and loss of accuracy. Hence, an appropriate ε should be chosen to balance a preferred terminal precision and bank reversals.


### 4.2.2. Adaptability Testing and Simulations

As a first step in testing and assessing the efficiency and adaptability of the guidance algorithm represented as AANLL, we present and discuss the results of all mission scenarios set up and listed in Table1. For demonstration and comparison, all initial conditions and guidance parameters were kept the same as mission 1 for all missions.The nominal terminal conditions for all the above mission scenarios are listed in Table2. Note that all missions have a rather high accuracy. The terminal altitude errors are all less than 1 km, the velocity errors are less than 5 m/s, and range errors are all less than 3 km. The heading errors are significantly less than the 5 deg requirement, the maximum of which is a rather small value of -0.52 deg. The peak heating rates for all missions are no more than the corresponding limits, respectively. In addition, the other two peak path constraints are also well satisfied. In fact, the peak aerodynamic loads and dynamic pressures for all missions did not exceed 115 kPa and 2.5, which are no more than the corresponding limits, respectively. Therefore, Table 2 only gives a focus on the comparison of peak heating rates.Table 2
Terminal condition precision for all missions.
MissionΔhf(m)ΔVf (m/s)Δsf (km)Δψf (deg)MaximumQ̇ (kW/m2)173.941.540.120.001294.582-268.470.010.180.261249.913991.72-7.522.880.031292.364-486.171.84-0.110.181295.705-666.201.982.42-0.521250.006-152.78-0.39-0.440.121293.217328.33-2.77-1.00-0.241249.788853.92-4.530.38-0.031278.489405.52-4.48-0.320.071249.64The computation time used for generating the entry trajectory for each mission with flight time of about 1500 s is only about 3-4 s. It should be noted that most time is consumed by integrating Eqs. (1)–(6) throughout the entire entry trajectory. In every guidance cycle, the computation time required to generate the commanded bank angle is dramatically less than 1 ms. Note that all computations were implemented on a laptop computer and the algorithm is coded in MATLAB without any optimization. Predictably, improvements in software and hardware could provide large room for improvement in computation speed. Thus, the guidance algorithm has an indubitable potential for onboard application.Taking missions 5, 6, 7, 8, and 9 as examples, Figure4 shows the comparison of altitude, ground track, bank angle, and heating rate histories for the QEG phase only. It is evident from Figure 4(a) that the large phugoid oscillations in the initial part are gradually eliminated along the trajectory propagations, which ensures the validity of the QEGC. The ground tracks shown in Figure 4(b) demonstrate the high position accuracy for different missions. The characteristic of phugoid oscillations is also verified by the flight path angles shown in Figure 2(c). It can be seen from Figure 2(d) that the entry trajectories reach the final site with different heading angle (i.e., in different directions) for different missions.Figure 4
Trajectory comparison for all missions.
(a)
Longitudinal profiles(b)
Ground tracks(c)
Flight path angles(d)
Heading errors(e)
Bank angles(f)
Heating ratesCombining with the heating rate profiles is shown in Figure4(f), and Figure 4(e) illustrates the validity of the translation of inequality path constraints. Look closely at the comparison of the bank angle histories for missions 8 and 9. One of the main differences lies in the initial bank angle when the QEG phase initiates. As expected, the magnitude of the bank angle changes considerably as the heating rate constraint is imposed on. Once the trajectory has entered the QEG phase, the magnitude of the bank angle of mission 8 will increase immediately to focus on the range requirement due to the loose heating rate limit. The path constraint is not active for this case. However, the magnitude of bank angle of mission 9 is still maintained to be zero for the initial small time range, so that the trajectory would be forced to drive shallower into the dense atmosphere initially. It is also confirmed by the altitude profiles as seen in Figure 4(a). Then, the peak heating rate, often found on the first trough point, is induced effectively as shown in Figure 4(f).To further test and assess the efficiency and adaptability of the guidance algorithm, 100 dispersion cases were studied for mission 1 using the Monte Carlo simulations. The dispersed initial conditions are considered and modeled by the zero-mean Gaussian dispersions with 3-sigma values. The 3-sigma value of initial entry condition dispersions are as follows: the dispersed altitude of 10 km, the longitude and latitude of 2 deg, the velocity of 300 m/s, the flight path angle of 0.2 deg, and the heading angle of 3 deg, respectively. Table3 summarizes the statistics on the final conditions for 100 dispersed trajectories. Obviously, the small final errors are all derived for each final state concerned in terms of small means and standard deviations. The efficiency, robustness, and adaptability are dramatically confirmed and approved again. However, in the presence of significant aerodynamic modeling uncertainties and atmosphere modeling dispersions (the atmosphere density, aerodynamic coefficients all have the 3-sigma values of 15% in respective nominal values), this algorithm performs not well with large scatters especially in the terminal range error and heading error. Hence, no results are given here. When the guidance algorithm serves as an entry guidance approach, some techniques (e.g., aerodynamics filter, compare [27]) can be used to improve the guidance performance. This comment will be carried out in future works about entry guidance.Table 3
Statistics on the final conditions for 100 dispersed trajectories.
hf (km)Vf (m/s)sf (km)Δψf (deg)Average25.062198.1949.670.12Max26.002202.6050.771.03Min24.272194.3146.50-0.60Standard deviation0.0401.770.830.23


## 4.2.1. Preliminary Testing and Discussions

Firstly, the test for mission 1 was implemented to verify the principles of the guidance algorithm and assess how well the algorithm works. A conventional bank reversal logic used by the Apollo and the Shuttle is demonstrated and validated. As can be seen in [29], the sign of bank angle should be maintained until the following criterion is violated
(70)Δψ≤ΔazmthV.An appropriate design ofΔazmthV, which avoids too many bank reversals, depends on the vehicle performance and the mission scenarios and is time-consuming. For comparison, a velocity-dependent ΔazmthV with a piecewise linear form, which is similar to that in [29], is carefully chosen to the above lateral logic by trial simulations.The tests of mission 1 were completed for three cases in terms of the baseline algorithm described in Eq. (38) (represented by BA), the augmented algorithm described in Eq. (58) with new lateral logic in Eq. (68) (represented by AANLL), and the same augmented algorithm but with conventional lateral logic (represented by AACLL). The comparison of the testing results for the three cases is shown in Figure 2.Figure 2
Trajectory comparison of the three cases for mission 1.
(a)
Entry trajectories(b)
Ground tracks(c)
Flight path angles(d)
Bank angles(e)
Crossranges(f)
Heading errorsFigures2(a) and 2(b) show the altitude versus range-to-go profiles and the ground tracks. Although the baseline algorithm has poor precision and large phugoid oscillations in the altitude, the augmented algorithm performs well for two lateral logics. The large phugoid oscillations are mostly eliminated with a high terminal precision by the feedback augmentation, which can also be seen from the flight path angle plotted in Figure 2(c). The AANLL gives a terminal altitude error of 73.94 m, a velocity error of 1.54 m/s, and a range error of 121.93 m. Based on the carefully chosen ΔazmthV, the AACLL gives the same level of precision but different ground track as shown in Figure 2(b) because of the different bank reversals. Also seen in Figure 2(c) is that the flight path angles of the AANLL are small negative values and vary rather smoothly except for the initial descent phase. Hence, the corresponding hypothesis for the QEGC is reasonable and effective as well as an approximately equilibrium glide.Figure2(d) shows significantly the different bank angle histories. The bank angle is maintained to be a constant zero till the initial descent phase terminates. In the QEG phase, the AANLL gives four bank reversals, which is less than the AACLL and leads to different ground tracks as can be seen in Figure 2(b). Figure 2(e) and 3(f) give the comparison of heading errors and crossranges for the above three cases. With the feedback augmentation in Eq. (58), the corresponding algorithms render the approximately piecewise linearity crossrange parameter acceptable and reasonable, in contrast to the volatilization and the nonlinearity of heading errors. The final heading error is only about 0.002 deg for AANLL, dramatically better than that of 3.08 deg for AACLL. This validates the high lateral precision that the AANLL can offer. Because the first bank reversal of the AANLG occurs considerably later than that of the AACLL, the maximum of the crossrange parameter is larger than that of the AACLL.Figure 3
Crossranges and bank angles of the AANLL for differentε.
(a)
Crossranges(b)
Bank anglesFigure3 compares the crossranges and bank angles with different scaling factors ε for the AANLL. Obviously, as ε is increased, the maximal lateral excursion increases while the bank reversals decrease, aside from the rearward shift of the reversal point in time and loss of accuracy. Hence, an appropriate ε should be chosen to balance a preferred terminal precision and bank reversals.


## 4.2.2. Adaptability Testing and Simulations

As a first step in testing and assessing the efficiency and adaptability of the guidance algorithm represented as AANLL, we present and discuss the results of all mission scenarios set up and listed in Table1. For demonstration and comparison, all initial conditions and guidance parameters were kept the same as mission 1 for all missions.The nominal terminal conditions for all the above mission scenarios are listed in Table2. Note that all missions have a rather high accuracy. The terminal altitude errors are all less than 1 km, the velocity errors are less than 5 m/s, and range errors are all less than 3 km. The heading errors are significantly less than the 5 deg requirement, the maximum of which is a rather small value of -0.52 deg. The peak heating rates for all missions are no more than the corresponding limits, respectively. In addition, the other two peak path constraints are also well satisfied. In fact, the peak aerodynamic loads and dynamic pressures for all missions did not exceed 115 kPa and 2.5, which are no more than the corresponding limits, respectively. Therefore, Table 2 only gives a focus on the comparison of peak heating rates.Table 2
Terminal condition precision for all missions.
MissionΔhf(m)ΔVf (m/s)Δsf (km)Δψf (deg)MaximumQ̇ (kW/m2)173.941.540.120.001294.582-268.470.010.180.261249.913991.72-7.522.880.031292.364-486.171.84-0.110.181295.705-666.201.982.42-0.521250.006-152.78-0.39-0.440.121293.217328.33-2.77-1.00-0.241249.788853.92-4.530.38-0.031278.489405.52-4.48-0.320.071249.64The computation time used for generating the entry trajectory for each mission with flight time of about 1500 s is only about 3-4 s. It should be noted that most time is consumed by integrating Eqs. (1)–(6) throughout the entire entry trajectory. In every guidance cycle, the computation time required to generate the commanded bank angle is dramatically less than 1 ms. Note that all computations were implemented on a laptop computer and the algorithm is coded in MATLAB without any optimization. Predictably, improvements in software and hardware could provide large room for improvement in computation speed. Thus, the guidance algorithm has an indubitable potential for onboard application.Taking missions 5, 6, 7, 8, and 9 as examples, Figure4 shows the comparison of altitude, ground track, bank angle, and heating rate histories for the QEG phase only. It is evident from Figure 4(a) that the large phugoid oscillations in the initial part are gradually eliminated along the trajectory propagations, which ensures the validity of the QEGC. The ground tracks shown in Figure 4(b) demonstrate the high position accuracy for different missions. The characteristic of phugoid oscillations is also verified by the flight path angles shown in Figure 2(c). It can be seen from Figure 2(d) that the entry trajectories reach the final site with different heading angle (i.e., in different directions) for different missions.Figure 4
Trajectory comparison for all missions.
(a)
Longitudinal profiles(b)
Ground tracks(c)
Flight path angles(d)
Heading errors(e)
Bank angles(f)
Heating ratesCombining with the heating rate profiles is shown in Figure4(f), and Figure 4(e) illustrates the validity of the translation of inequality path constraints. Look closely at the comparison of the bank angle histories for missions 8 and 9. One of the main differences lies in the initial bank angle when the QEG phase initiates. As expected, the magnitude of the bank angle changes considerably as the heating rate constraint is imposed on. Once the trajectory has entered the QEG phase, the magnitude of the bank angle of mission 8 will increase immediately to focus on the range requirement due to the loose heating rate limit. The path constraint is not active for this case. However, the magnitude of bank angle of mission 9 is still maintained to be zero for the initial small time range, so that the trajectory would be forced to drive shallower into the dense atmosphere initially. It is also confirmed by the altitude profiles as seen in Figure 4(a). Then, the peak heating rate, often found on the first trough point, is induced effectively as shown in Figure 4(f).To further test and assess the efficiency and adaptability of the guidance algorithm, 100 dispersion cases were studied for mission 1 using the Monte Carlo simulations. The dispersed initial conditions are considered and modeled by the zero-mean Gaussian dispersions with 3-sigma values. The 3-sigma value of initial entry condition dispersions are as follows: the dispersed altitude of 10 km, the longitude and latitude of 2 deg, the velocity of 300 m/s, the flight path angle of 0.2 deg, and the heading angle of 3 deg, respectively. Table3 summarizes the statistics on the final conditions for 100 dispersed trajectories. Obviously, the small final errors are all derived for each final state concerned in terms of small means and standard deviations. The efficiency, robustness, and adaptability are dramatically confirmed and approved again. However, in the presence of significant aerodynamic modeling uncertainties and atmosphere modeling dispersions (the atmosphere density, aerodynamic coefficients all have the 3-sigma values of 15% in respective nominal values), this algorithm performs not well with large scatters especially in the terminal range error and heading error. Hence, no results are given here. When the guidance algorithm serves as an entry guidance approach, some techniques (e.g., aerodynamics filter, compare [27]) can be used to improve the guidance performance. This comment will be carried out in future works about entry guidance.Table 3
Statistics on the final conditions for 100 dispersed trajectories.
hf (km)Vf (m/s)sf (km)Δψf (deg)Average25.062198.1949.670.12Max26.002202.6050.771.03Min24.272194.3146.50-0.60Standard deviation0.0401.770.830.23


## 5. Conclusions

A simple, adaptive, and autonomous guidance algorithm is developed for entry vehicles with a high L/D ratio. The novel utilization of the QEGC is the cornerstone of this guidance algorithm for fully constrained, three-dimensional feasible entry flight. The algorithm for tackling the problem contains two parts: the longitudinal profile guidance and the lateral guidance. The longitudinal guidance generates the feasible magnitude of the bank angle analytically and successively in real-time, while the lateral guidance determines the sign of the bank angle by a simple but efficient bank reversal logic. Conducting simultaneously these two channels with the successive states updated, the set of closed-loop commanded bank angle with analytic feedback laws are easily deduced. The inequality path constraints in the velocity-altitude space are also analytically translated into the velocity-dependent bounds for the magnitude of the bank angle by the QEGC. Because no iterations and few off-line parameter adjustments are necessary, the algorithm provides remarkable simplicity, rapidity, and adaptability. A considerable range of entry flights using the vehicle data of the CAV-H is tested. Simulation results demonstrate the effectiveness and performance of the presented approach. Accordingly, a feasible and applicable entry trajectory is generated by integrating the whole trajectory only once with a pleasing computation cost so that the guidance algorithm can also serve as an analytical trajectory planning method. In the future, entry trajectory generation with waypoints and no-fly zones will be carried out using this presented algorithm. The performance will also be tested and assessed when the algorithm serves as an entry trajectory planning approach.



---

*Source: 8874251-2020-11-18.xml*