Automatic Approach and Landing Trajectory Planner for Unpowered Reusable Launch Vehicle

A new guidance scheme for the approach and landing (A & L) phase of an unpowered reusable launch vehicle (RLV) has been developed. The main advantage of the new guidance is the use of glide-efficiency factor as the independent variable to compute the geometrical flare parameters by a set of analytical functions. The trajectory-planning algorithm generates its reference geometry based on the steep and shallow subphases, respectively. During the steep segment, the quasi-equilibrium glide (QEG) solution, which assumes a constant dynamic pressure and flight-path angle during the flight, is used to create the flight reference while the shallow segment is defined by polynomial functions for altitude and dynamic pressure profiles. Standard linearization methods are used to design a closed-loop command in order to track the QEG profile. Furthermore, proportion-derivative (PD) control is used to modulate the lift coefficient during the flare flight. Once the reference trajectory is created, a closed-loop simulation is obtained to track the reference. Off-nominal conditions, in terms of change in initial glide-efficiency factor, dynamic pressure, flight-path angle, and altitude are tested using a Monte-Carlo simulation. The simulated results demonstrate the effectiveness of the proposed algorithm to land the vehicle successfully under large dispersions of glide-efficiency factor.


Introduction
The National Aeronautics and Space Administration (NASA) has sought to de- velop advanced guidance and controls methods in an effort to improve safety and reliability of future reusable launch vehicle (RLV) trajectories [1] [2] [3].
These innovative technologies may apply to an unpowered winged lifting body or an unmanned booster vehicle.In either case, the RLV's aerodynamic parameters are limited, such as low lift-to-drag (L/D) ratio.Hence, gliding from the current state to the desired state is a critical task for a successful landing, especially under large initial conditions dispersions.The RLV needs to replan or reshape its trajectory to achieve the desired runway touchdown state.
Traditionally, the RLV descent can be divided into three segments: Entry, terminal area energy management (TAEM), and approach and landing (A & L) [4].The approach and landing phase represents a critical flight phase because the RLV needs to manage its energy in order to glide and meet runway touchdown conditions.In general, the A & L phase begins at the end of the TAEM phase at a low subsonic Mach number (M = 0.5) and an altitude of 10,000 ft above the runway.In the past, the A & L flight used reference altitude profiles which consist of steep and shallow glideslopes [5].This two-phase flight trajectory has been proven to be successful for low L/D unpowered vehicles such as the US Space Shuttle.The Shuttle, however, relied on a small set of fixed reference profiles, and therefore the shuttle guidance may not be well-suited in the presence of large trajectory, atmospheric, aerodynamics, and vehicle-mass dispersions.
Several researchers have developed guidance systems for the A & L flight phase.Some of these algorithms have been proposed to replan the reference trajectory under initial condition dispersions.Kluever [6] developed an A & L predicator-corrector guidance system that reshapes the reference trajectory by adjusting the glide-efficiency factor during the mission.Kluever [7] also presented an A & L guidance method based on limited normal acceleration capabilities.
The reference path with minimum load factor is generated by iterating on the initial flight-path angle until achieving the desired touchdown sink rate.Kluever [8] developed an onboard trajectory-planning algorithm that recomputes a new reference trajectory based on the wind conditions, the energy state, and the aerodynamic performance.Harl et al. [9] and Zhao et al. [10] presented an A & L guidance method that computes an online reference path based on sliding mode control.
Other studies have considered optimization methods to design the A & L guidance system.Schierman et al. [11] [12] [13] proposed an A & L guidance system that creates a set of optimal trajectories offline based on indirect optimization techniques.Then, the best trajectory is selected from the family of optimum paths to guide the RLV from the current state to the desired state.Trent et al. [14] developed a trajectory-planning algorithm that applies Pontryagin's minimum principle to solve a two-point boundary value problem (TPBVP).

Equations of Motion
The unpowered RLV is assumed to be a point mass, and its governing equations of motion are where V is the Earth-relative velocity, γ is the flight-path angle, h is the altitude above the runway, x represents the downtrack position of the vehicle, g is the gravitational acceleration, and m is the mass.The aerodynamic lift force L and drag force D can be written as where C L and C D are lift and drag coefficients, and S is the vehicle's reference area.The dynamic pressure is Time t is the independent variable in equations of motion (1-4); however, our Advances in Aerospace Science and Technology algorithm uses altitude as the independent variable.Thus, the chain-rule method is applied to replace time with altitude as the independent variable.Dividing Equations ( 1), (2), and (4) by Equation ( 3) we obtain

Vehicle Model
A low L/D gliding vehicle is used for guidance algorithm development, and the vehicle data is taken from [6].The specified RLV has wing area S = 2370 ft 2 , mass 4211.5 slugs m = ( 1355 bf 00 I W = ), and the wing loading is 57.2 psf W S = .
The aerodynamic drag coefficient is computed using the standard drag polar equation: where 0 D C and K are the zero-lift and lift-induced drag coefficients, respec- tively.Figure 1 in Ref [6] shows that the aerodynamic coefficients The best fit of the zero-lift coefficient can be described by a rational (5 th -order/4 th -order) function for 0.
The best fit of the induced-drag coefficient can be described by a rational (2 th -order/5 th -order) function for 0.5 M ≤ : Although the zero-lift 0 D C and the induced-drag K coefficients are obtained by complicated equations, there is very little change in the drag polar parameters with Mach number at low subsonic flight.Table 1 presents the drag polar parameters at three discrete Mach numbers during the A & L phase.
The guidance algorithm estimates the drag coefficient D C using Equation ( 11) by using the commanded value of L C , where 0 D C and K are computed using Equations ( 12) and ( 13) for Mach number M .

Reference Trajectory
Unlike most trajectory-planning algorithms, our new guidance does not require online numerical integration of the equations of motion.Instead, onboard computation is performed based on the available vehicle glide-efficiency factor η to generate a new reference trajectory.The glide-efficiency factor is defined as a lift-to-drag ratio (L/D) divided by the maximum lift-to-drag ratio ( ) max

L D .
Once a reference trajectory is created, a closed-loop simulation is propagated to track the reference trajectory relying on the steep and shallow flight controls.At this point, it is important to describe the steps for generating the new reference during the steep and shallow segments:

Steep Subphase
The QEG condition serves as the reference profile for dynamic pressure and flight-path angle.It is unrealistic to assume that velocity and flight-path angle are constant; however, it is realistic to expect that the dynamic pressure and flight-path angle remain constant during the majority of A & L. For the QEG condition, the gliding flight is steady with no change in dynamic pressure and flight-path angle.The time derivative of dynamic pressure is Substituting Equation (1) and Equation (3) into Equation ( 14), we obtain Note that d dh ρ ρ ′ = is the change in density with altitude; it can be computed using the U.S 1976 Standard Atmospheric model.
Substituting Equation ( 5) into Equation ( 2), we obtain Hence, the QEG condition consists of two nonlinear equations, Equations ( 15) and ( 16) with four free variables: Maximum L/D can be computed from the drag polar parameters: The glide-efficiency factor is used to adjust the vehicle range, in which flight at higher glide-efficiency produces lower dynamic pressure and shallower flight-path angle.Hence 1 η = results the greatest range.In contrast, flight at lower glide-efficiency produces higher dynamic pressure and steeper flight-path angle and reduces range.
The only way to solve the QEG condition of two equations [Equations ( 17)-( 18)] with four unknown variables is to fix two variables and compute the remaining two unknown variables. In

Shallow Subphase
Before discussing how the shallow glideslope reference trajectory is generated, it is instructive to compute the flare geometrical parameters such as the flare altitude, h F , and the downtrack location of the start of flare, s F .Both are functions of glide-efficiency factor [6].As with ALI x , the initial flare altitude h F and the downtrack location of the start of flare phase s F are fitted also using the cftool function in MATLAB by a rational (4 th -order/5 th -order) and (2 th -order/5 th -order) functions, respectively.The h F and s F profiles are shown in Figure 3 and Figure 4, respectively.1 ft 0 Figure 2. Downtrack approach and landing interface location vs. glide efficiency.A fourth-order polynomial altitude profile defines the flare maneuver phase where x is the downtrack position along the runway centerline.The derivative of altitude with respect downtrack is the tangent of the flare maneuver flight-path angle: Fourth-order coefficients 0 a and 1 a are computed from Equation (27) and Equation (28) and the RLV's altitude and flight-path angle at the start of flare phase where 1400 ft A quadratic polynomial profile defines dynamic pressure during the flare maneuver phase: The derivative of dynamic pressure with respect downtrack is: The quadratic coefficients o c and 1 c are computed from Equation (29) and Equation (30) and the RLV's dynamic pressure and its derivative at the start of flare phase where 1400 ft that QEG q is constant during the steep phase and hence QEG 0 q =  .The remaining coefficient 2 c is computed from Equation (29) at the touchdown point where TD 1400 ft x x = = and 2 TD SSL TD 0.5 Where the standard sea level air density ( )

Approach and Landing Simulation
In this section, we will derive the guidance command required to track the dynamic pressure and flight-path angle during the steep phase and to track the flare altitude profile f h during the shallow phase.During the QEG and flare phases, the lift coefficient is determined using the closed-loop command: The closed-loop control L C guides the RLV so that it tracks the reference dynamic pressure and flight-path angle.The feedback gains K γ and q K are designed based on the pole-placement technique so that the system has very good damping with fast response.Equation (15), Equation ( 16), and Equation (3), have been linearized with respect to QEG path by considering the dynamic During the shallow phase, the reference lift coefficient L C is evaluated using Equation (32); however the closed-loop feedback term can be formulated by: where the feedback gains H K and HD K are linear profiles with altitude which are determined by trial and error.

Numerical Results
The purpose of dispersion simulations is to evaluate the performance of the guidance algorithm in the presence of significant deviations in glide-efficiency factor η and trajectory states (dynamic pressure, flight-path angle, and alti- tude).Nominal and off-nominal guided trajectories are obtained by numerically integrating Equations (7-10) using a variable-step, variable-order solver (MATLAB's ode45).The lift coefficient is calculated using the open-and closed-loop commands in order track the reference trajectory.
In summary, the QEG solution is used to create the steep flight reference for a known initial glide-efficiency factor and its corresponded dynamic pressure and flight-path angle [Equations ( 15)-( 16

Nominal A & L Trajectory
The nominal initial conditions for the A & L phase are taken from [16].These initial A & L states corresponded to the nominal end-state from our TAEM trajectory studies in [16].These nominal values are 0.835

Off-Nominal A & L Trajectory
Many dispersed A & L trajectories are tested in order to demonstrate the robustness of our trajectory-planning algorithm.In this study, we impose random off-nominal initial conditions to changes in initial glide-efficiency factor, flight-path angle, dynamic pressure, and the altitude.The statistics of the random A & L initial conditions used in this paper are based on the end-conditions of the guided TAEM phase [16], which are summarized in Table 3.The statistics of all runway touchdown state errors are summarized in Table 4.
Figure 8 and Table 4 present a summary of the Monte-Carlo simulation results.These statistical results reveal that the minimum and maximum values are close to the target conditions at the runway touchdown, and all the standard

.
The atmospheric density ρ is computed using the US 1976 Standard Atmosphere.The governing equations of motion (1-4) are formulated with respect to a "flat-Earth" model, where the + x axis points along the runway centerline and the runway threshold is the origin point.Figure1presents the force vectors and flight-path angle during the A & L phase.

Figure 1 .
Figure 1.Approach and landing coordinate frame.
18) where L C is the open-loop lift or reference coefficient.The required lift coefficient * L C and drag coefficient * D C during the QEG flight are determined using the lift and drag coefficients corresponding to maximum L D , which can be computed by:

Figure 4 .
Figure 4. Start of flare location vs. glide efficiency factor.
angle at the end of the steep phase).The second-and third-order polynomial coefficients 2 a and 3 a are also computed from Equation (27) and Equation (28) and the RLV's altitude and flight-path angle at the touchdown L has two separate altitude profiles, the RLV has a smooth transition from the steep phase to the shallow phase due to a continuous and differentiable reference path.
factor, η Initiial flare altitude,h f , ft Advances in Aerospace Science and Technology At this point, we can substitute Equation (5) into Equation (2) and rewrite the result using the chain-rule for the time rate of flight-path angle: into Equation (31) and solving for lift coefficient, we obtain the open-loop lift coefficient:

.Figures 5
Figures 5(a)-5(d) show the flight-path angle, dynamic pressure, altitude, and lift coefficient reference profiles for A & L altitude ranging from 10,000 ft to zero, respectively, (note that A & L progress right-to-left in Figure 5.We expect the RLV to complete the steep phase and begin the shallow phase at altitudes below the initial flare altitude.The reference profiles are obtained for four values of glide-efficiency factor, where corresponds to flight with is that 74% of the maximum lift-to-drag ratio, etc. Figure 5(a) and Figure 5(b) indicate that the flight-path angle and dynamic pressure profiles are nearly constant during the steep phase due to the QEG solutions.However and significantly decrease during the shallow flare phase in order to achieve the touchdown states.As expected, the flight at higher glide-efficiency factor produces a shallower flight-path angle and lower dynamic pressure.Figure 5(c) shows that the altitude vs. downrange profile has a linear behavior during the steep segment due to a constant flight-path angle; however altitude has a fourth-order profile during the flare maneuver according to Equation (27).As we can see from Figure 5(c), higher glide-efficiency produces higher range.Because q is nearly constant, the open-loop lift coefficient is approximately constant during the steep phase.Lift coefficient increases significantly to pull up the RLV toward the touchdown state as shown in Figure 5(d).The open-loop lift coefficient for the pull-up flare maneuver should be calculated in order to find the drag coefficient [Equation (11)].All reference calculations are performed online for an initial glide-efficiency factor and desired runway touchdown conditions.Once a glide efficiency factor
During the steep phase, the closed-loop lift coefficient L C is designed to track the reference dynamic pressure and flight-path angle profiles which are obtained by a two-dimensional interpolation of the stored QEG solutions.The open-loop lift coefficient L C is determined using Equation (23) and the feedback term is

η
= 0.93 Advances in Aerospace Science and Technology pressure and flight-path angle as the state variables with altitude as the independent variable.The gains are scheduled with altitude by solving a sequence of pole-placements problems along the steep phase.
)].Then a shallow reference profile is defined by a fourth-order polynomial profile for altitude and quadratic polynomial profile for dynamic pressure [Equations (27)-(28)].All the A & L geometrical parameters (the downtrack location of the start of A & L phase ALI x , flare alti- tude F h , and downtrack location of the start of flare F s ) are computed based on the initial glide-efficiency factor [Equations (24)-(26)].All of the reference polynomial coefficients are calculated based on the A & L geometrical parameters and the desired runway touchdown conditions [Equations (27)-(30)].Finally, open-and closed-loop commands are obtained in order to track the reference profiles under nominal and off-nominal conditions [Equation (23), Equation (32), Equation (34), Equation (35), and Equation (36)].

Figures 6
Figures 6(a)-6(f), show the nominal flight-path angle, dynamic pressure, altitude, closed-loop lift command, load factor, and sink velocity profiles, respectively.Figures 6(a)-6(d) show that the flight-path angle, dynamic pressure, ground track range, and lift coefficient references are tracked accurately without excessive closed-loop lift.Figure 6(e) shows that at the beginning of the circular flare (at altitude h = 1,028 ft), the RLV requires a load-factor maneuver in order
in Aerospace Science and Technology A Monte-Carlo simulation is used to test the proposed trajectory-planning algorithm by considering 1000 guided trajectories based on the random variations in the initial states as shown in Table 3. Figures 7(a)-7(c) present the 1000 flight path-angle, dynamic pressure, and ground-track range histories.Although the dynamic pressure, the flight-path angle and the altitude profiles show different initial values, they reach the same touchdown conditions ( 0 h = and 1400 ft x =).

Figures
Figures (8a)-(8d) shows the statistical results of the flight-path angle, dynamic pressure, downrange, and touchdown sink velocity errors, respectively, at the runway touchdown (TD).It can be seen from Figure 8 that all runway touchdown state errors are reasonably small.

Table 2
presents maximum L/D at three discrete values of Mach number.The lift and drag coefficients of the flight with ( ) max L D (or 1 η = ) can be com- puted directly using Equations (19) and (20), however, these equations are not suitable for all glide efficiency factors.Therefore, we need to derive a general equation to find the open-loop lift coefficient for all values of L/D.As previously mentioned, the glide-efficiency factor can be describe as

Table 2 .
Maximum L/D vs. Mach number.
our work, the altitude h is fixed and the lift coefficient L [6]alues.These profiles are obtained for altitude ranging from 10,000 ft (beginning of A & L) to the initial flare altitude h F (beginning of shallow glide subphase).The initial flare altitude h F is a function of glide-efficiency factor; it will be described in the next section.At this point, it is appropriate to compute the downtrack location of the start of the A & L phase, ALIx .The location represents the matching point between the TAEM and A & L phases, which is also a function of glide-efficiency factor.Following the methods of Ref[6], the ALI x profile is computed by numerically integrating the velocity differential equation backward from the touchdown state to the start of the A & L phase for different glide-efficiency factors.The downtrack location of the start of the A & L interface (ALI) is fitted using MATLAB's curve fitting toolbox (cftool) by a rational (3 th -order /5 th -order) function as shown in Figure 2.

Table 3 .
Statistics for initial state errors and glide efficiency.

Table 4 .
Statistics for runway touchdown state errors.