Automatic Approach and Landing Trajectory Planner for Unpowered Reusable Launch Vehicle ()
1. Introduction
The National Aeronautics and Space Administration (NASA) has sought to develop 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). Heydari et al [15] developed an A & L guidance algorithm that introduces a state-dependent Riccati equation (SDRE) to solve the finite-horizon optimal problem.
This paper presents an A & L guidance methodology that utilizes a trajectory planning algorithm. This work differs from the prior work in two ways: a) planning a new reference trajectory does not rely on numerical integration of the governing equations of motion; b) our algorithm is able to successfully guide the RLV despite simultaneous large variations in initial downrange, altitude, dynamic pressure, and flight-path angle. Our proposed algorithm generates its reference based on the quasi-equilibrium glide (QEG) solutions and a set of analytical expressions. In the steep phase, the QEG solutions (dynamic pressure, flight-path angle, and open-loop lift coefficient) are stored to generate a two-dimensional look-up table with respect to altitude. The shallow phase is defined by a fourth-order polynomial and a quadratic polynomial for altitude and dynamic pressure profiles, respectively.
All reference calculations are performed online and the open- and closed-loop commands are readily available after the reference trajectory is updated. Finally, the effectiveness of the proposed algorithm is validated using a Monte-Carlo simulation method.
2. System Models
2.1. Equations of Motion
The unpowered RLV is assumed to be a point mass, and its governing equations of motion are
(1)
(2)
(3)
(4)
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
(5)
(6)
where CL and CD are lift and drag coefficients, and S is the vehicle’s reference area. The dynamic pressure is
. 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. Figure 1 presents the force vectors and flight-path angle during the A & L phase.
Time t is the independent variable in equations of motion (1-4); however, our
Figure 1. Approach and landing coordinate frame.
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
(7)
(8)
(9)
(10)
2.2. 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 ft2, mass
(
), and the wing loading is
. The aerodynamic drag coefficient is computed using the standard drag polar equation:
(11)
where
and
are the zero-lift and lift-induced drag coefficients, respectively. Figure 1 in Ref [6] shows that the aerodynamic coefficients
and
are functions of Mach number. We fit
and
with analytical functions using MATLAB’s curve fitting toolbox (cftool).
The best fit of the zero-lift coefficient can be described by a rational (5th-order/4th-order) function for
:
(12)
The best fit of the induced-drag coefficient can be described by a rational (2th-order/5th-order) function for
:
(13)
Although the zero-lift
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
using Equation (11) by using the commanded value of
, where
and
are computed using Equations (12) and (13) for Mach number
.
3. 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
. 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:
3.1. 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
(14)
Substituting Equation (1) and Equation (3) into Equation (14), we obtain
Table 1. RLV Drag polar parameters.
(15)
Note that
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
(16)
Hence, the QEG condition consists of two nonlinear equations, Equations (15) and (16) with four free variables:
(17)
(18)
where
is the open-loop lift or reference coefficient. The required lift coefficient
and drag coefficient
during the QEG flight are determined using the lift and drag coefficients corresponding to maximum
, which can be computed by:
(19)
(20)
Maximum L/D can be computed from the drag polar parameters:
(21)
Table 2 presents maximum L/D at three discrete values of Mach number. The lift and drag coefficients of the flight with
(or
) can be computed 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
(22)
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
results the greatest range. In contrast, flight at lower glide-efficiency produces higher dynamic pressure and steeper flight-path angle and reduces range.
Table 2. Maximum L/D vs. Mach number.
Substituting Equation (21) into Equation (22) and solving for the reference lift coefficient
so that the RLV flies on the “front side” of the
curve, we obtain
(23)
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 our work, the altitude
is fixed and the lift coefficient
is computed from Equation (23) for a fixed value of
. The QEG solutions are obtained using a Newton-Raphson method that typically needs 6 - 7 iterations to determine the
and
values. These profiles are obtained for altitude ranging from 10,000 ft (beginning of A & L) to the initial flare altitude hF (beginning of shallow glide subphase). The initial flare altitude hF 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,
. 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
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 (3th-order /5th-order) function as shown in Figure 2.
(24)
3.2. 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, hF, and the downtrack location of the start of flare, sF. Both are functions of glide-efficiency factor [6] . As with
, the initial flare altitude hF and the downtrack location of the start of flare phase sF are fitted also using the cftool function in MATLAB by a rational (4th-order/5th-order) and (2th-order/5th-order) functions, respectively. The hF and sF profiles are shown in Figure 3 and Figure 4, respectively.
(25)
(26)
Figure 2. Downtrack approach and landing interface location vs. glide efficiency.
Figure 3. Initial flare altitude vs. glide efficiency factor.
while it is important to know the initial conditions at the start of A & L, it is also important to know the touchdown conditions. The desired touchdown state can be defined as
,
, and
.
Figure 4. Start of flare location vs. glide efficiency factor.
A fourth-order polynomial altitude profile defines the flare maneuver phase
(27)
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:
(28)
Fourth-order coefficients
and
are computed from Equation (27) and Equation (28) and the RLV’s altitude and flight-path angle at the start of flare phase where
,
and
(note that
is the flight-path angle at the end of the steep phase). The second-and third-order polynomial coefficients
and
are also computed from Equation (27) and Equation (28) and the RLV’s altitude and flight-path angle at the touchdown point where
,
, and
.
Although A & 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.
A quadratic polynomial profile defines dynamic pressure during the flare maneuver phase:
(29)
The derivative of dynamic pressure with respect downtrack is:
(30)
The quadratic coefficients
and
are computed from Equation (29) and Equation (30) and the RLV’s dynamic pressure and its derivative at the start of flare phase where
,
and
(note that
is constant during the steep phase and hence
. The remaining coefficient
is computed from Equation (29) at the touchdown point where
and
.
Where the standard sea level air density
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:
(31)
Substituting Equation (4) into Equation (31) and solving for lift coefficient, we obtain the open-loop lift coefficient:
(32)
where
. The downtrack derivative of the flight-path angle
can be determined from Equation (28)
(33)
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
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
Figure 5.Reference trajectory profiles for different glide-efficiency factor; (a) flight-path angle vs. altitude; (b) dynamic pressure vs. altitude; (c) altitude vs. downrange; (d) open-loop lift coefficient vs. altitude.
is initiated at the beginning of A & L phase, all the reference polynomial coefficients are calculated to establish the A & L trajectory to touchdown.
4. 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
during the shallow phase. During the QEG and flare phases, the lift coefficient is determined using the closed-loop command:
(34)
During the steep phase, the closed-loop lift coefficient
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
is determined using Equation (23) and the feedback term is
(35)
The closed-loop control
guides the RLV so that it tracks the reference dynamic pressure and flight-path angle. The feedback gains
and
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 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.
During the shallow phase, the reference lift coefficient
is evaluated using Equation (32); however the closed-loop feedback term can be formulated by:
(36)
where the feedback gains
and
are linear profiles with altitude which are determined by trial and error.
5. 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 altitude). 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)]. 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
, flare altitude
, and downtrack location of the start of flare
) 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)].
5.1. 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
,
,
, and
.
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 to track the rapid pulled-up reference trajectory. Finally, the sink velocity has a small change along the steep phase due to a small change in the RLV velocity, and eventually it increases significantly to meet the touchdown state as shown in Figure 6(f).
5.2. 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.
Figure 6. Nominal approach and landing profiles for glide-efficiency for
with quasi-equilibrium glide initial conditions; (a) flight-path angle vs. altitude; (b) dynamic pressure vs. altitude; (c) altitude vs. downrange; (d) closed-loop lift coefficient vs. altitude; (e) load factor vs. altitude; (f) sink velocity vs. altitude.
Table 3. Statistics for initial state errors and glide efficiency.
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 (
and
).
Figure 7(c) presents altitude versus downrange for several initial ranges, starting from (
) corresponded to flight with
to (
) corresponded to flight with
. Figure 7(d) and Figure 7(e) show also the 1000 closed-loop lift coefficient and load factor histories, respectively. These figures show that the actual trajectories track the references without excess the closed-loop lift coefficient effort. In addition, the actual lift-coefficient
at runway touchdown matches well with the reference lift coefficient
for all glide efficiency factors. Finally, Figure 7(f) presents the actual sink velocity histories which also satisfy the touchdown condition.
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.
Figure 7. Approach and landing profiles for 1000 Monte-Carlo simulations; (a) flight-path angle vs. altitude; (b) dynamic pressure vs. altitude; (c) altitude vs. downrange; (d) closed-loop lift coefficient vs. altitude; (e) load factor vs. altitude; (f) sink velocity vs. altitude.
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
Figure 8. Histograms of touchdown errors from 1000 Monte-Carlo simulations; (a) touchdown flight-path angle errors; (b) touchdown dynamic pressure errors; (c) touchdown downrange errors; (d) touchdown sink velocity errors.
Table 4. Statistics for runway touchdown state errors.
deviations are small. Hence, these results show the effectiveness and robustness of the designed algorithm to guide the RVL from the A & L initiation point to the runway touchdown point.
6. Conclusions
A new guidance scheme for the approach and landing phase of an unpowered reusable launch (RLV) has been developed. The guidance system employs a trajectory-planning algorithm that creates the reference trajectory by relying on the initial glide-efficiency factor. The approach and landing trajectory consists of two-phase (steep and shallow) reference flight. During the steep segment, the quasi-equilibrium glide (QEG) solution is used to create the flight reference with altitude as the independent variable while the shallow segment is defined by a fourth-order polynomial altitude and a quadratic dynamic pressure profiles.
Several closed-loop simulations using a Monte-Carlo method were obtained to validate the proposed guidance algorithm. Trajectories with random off-nominal conditions, (with changes in initial glide-efficiency factor, dynamic pressure, flight-path angle, and altitude) were simulated, and the results demonstrate the robust performance of the proposed algorithm.
Acknowledgements
The first author thanks Mhawesh Rafal for her suggestions and help.