Automatic Approach and Landing Trajectory Planner for Unpowered Reusable Launch Vehicle

Show more

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

$\stackrel{\dot{}}{V}=-\frac{D}{m}-g\mathrm{sin}\gamma $ (1)

$\stackrel{\dot{}}{\gamma}=\frac{L}{mV}-\frac{g}{V}\mathrm{cos}\gamma $ (2)

$\stackrel{\dot{}}{h}=V\mathrm{sin}\gamma $ (3)

$\stackrel{\dot{}}{x}=V\mathrm{cos}\gamma $ (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

$L=\stackrel{\xaf}{q}S{C}_{L}$ (5)

$D=\stackrel{\xaf}{q}S{C}_{D}$ (6)

where C_{L} and C_{D} are lift and drag coefficients, and S is the vehicle’s reference area. The dynamic pressure is
$\stackrel{\xaf}{q}=0.5\rho {V}^{2}$ . The atmospheric density
$\rho $ 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

$\frac{\text{d}V}{\text{d}h}=\frac{-D}{mV\mathrm{sin}\gamma}-\frac{g}{V}$ (7)

$\frac{\text{d}\gamma}{\text{d}h}=\frac{L}{m{V}^{2}\mathrm{sin}\gamma}-\frac{g}{{V}^{2}}\mathrm{cot}\gamma $ (8)

$\frac{\text{d}x}{\text{d}h}=\mathrm{cot}\gamma $ (9)

$\frac{\text{d}t}{\text{d}h}=\frac{1}{V\mathrm{sin}\gamma}$ (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 ft^{2}, mass
$m=4211.5\text{\hspace{0.17em}}\text{slugs}$ (
$W=135500\text{\hspace{0.17em}}\text{Ibf}$ ), and the wing loading is
$W/S=57.2\text{\hspace{0.17em}}\text{psf}$ . The aerodynamic drag coefficient is computed using the standard drag polar equation:

${C}_{D}\left(M,{C}_{L}\right)={C}_{D0}+K{C}_{L}^{2}$ (11)

where ${C}_{D0}$ and $K$ are the zero-lift and lift-induced drag coefficients, respectively. Figure 1 in Ref [6] shows that the aerodynamic coefficients ${C}_{D0}$ and $K$ are functions of Mach number. We fit ${C}_{D0}\left(M\right)$ and $K\left(M\right)$ with analytical functions using MATLAB’s curve fitting toolbox (cftool).

The best fit of the zero-lift coefficient can be described by a rational (5^{th}-order/4^{th}-order) function for
$M\le 0.5$ :

$\begin{array}{l}{C}_{D0}=\frac{0.0152{M}^{5}+0.0245{M}^{4}-0.1456{M}^{3}+0.225{M}^{2}-0.1616M+0.047}{{M}^{4}-3.329{M}^{3}+4.76{M}^{2}-3.376M+0.9794},\text{\hspace{0.17em}}\\ M\le 0.5\end{array}$ (12)

The best fit of the induced-drag coefficient can be described by a rational (2^{th}-order/5^{th}-order) function for
$M\le 0.5$ :

$K=\frac{4.427{M}^{2}-6.368M+2.633}{{M}^{5}-3.996{M}^{4}+2.722{M}^{3}+20.33{M}^{2}-30.2M+12.51},\text{\hspace{0.17em}}M\le 0.5$ (13)

Although the zero-lift ${C}_{D0}$ 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 ${C}_{D}$ using Equation (11) by using the commanded value of ${C}_{L}$ , where ${C}_{D0}$ and $K$ are computed using Equations (12) and (13) for Mach number $M$ .

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 $\eta $ 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 ${\left(L/D\right)}_{\mathrm{max}}$ . 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

$\stackrel{\dot{}}{\stackrel{\xaf}{q}}=\frac{1}{2}\frac{\text{d}\rho}{\text{d}h}\frac{\text{d}h}{\text{d}t}{V}^{2}+\rho V\frac{\text{d}V}{\text{d}t}$ (14)

Substituting Equation (1) and Equation (3) into Equation (14), we obtain

Table 1. RLV Drag polar parameters.

$\stackrel{\dot{}}{\stackrel{\xaf}{q}}=V\mathrm{sin}\gamma \left[\left(\frac{{\rho}^{\prime}}{\rho}-\frac{\rho S{C}_{D}}{m\mathrm{sin}\gamma}\right)\stackrel{\xaf}{q}-\rho g\right]=0$ (15)

Note that ${\rho}^{\prime}=\text{d}\rho /\text{d}h$ 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

$\stackrel{\dot{}}{\gamma}=\frac{\stackrel{\xaf}{q}S{C}_{L}}{mV}-\frac{W\mathrm{cos}\gamma}{mV}=0$ (16)

Hence, the QEG condition consists of two nonlinear equations, Equations (15) and (16) with four free variables:

$\stackrel{\dot{}}{\gamma}={f}_{1}\left(\stackrel{\xaf}{q},\gamma ,h,{\stackrel{\xaf}{C}}_{L}\right)=0$ (17)

$\stackrel{\dot{}}{\stackrel{\xaf}{q}}={f}_{2}\left(\stackrel{\xaf}{q},\gamma ,h,{\stackrel{\xaf}{C}}_{L}\right)=0$ (18)

where ${\stackrel{\xaf}{C}}_{L}$ is the open-loop lift or reference coefficient. The required lift coefficient ${C}_{L}^{*}$ and drag coefficient ${C}_{D}^{*}$ during the QEG flight are determined using the lift and drag coefficients corresponding to maximum $L/D$ , which can be computed by:

${C}_{L}^{*}=\sqrt{{C}_{D0}/K}$ (19)

${C}_{D}^{*}=2{C}_{D0}$ (20)

Maximum L/D can be computed from the drag polar parameters:

${\left(L/D\right)}_{\mathrm{max}}=\frac{1}{2\sqrt{{C}_{D}{}_{0}K}}$ (21)

Table 2 presents maximum L/D at three discrete values of Mach number. The lift and drag coefficients of the flight with ${\left(L/D\right)}_{\mathrm{max}}$ (or $\eta =1$ ) 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

$\eta =\frac{L/D}{{\left(L/D\right)}_{\mathrm{max}}}$ (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 $\eta =1$ 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 ${\stackrel{\xaf}{C}}_{L}$ so that the RLV flies on the “front side” of the $L/D$ curve, we obtain

${\stackrel{\xaf}{C}}_{L}=\frac{1}{2}\left\{\frac{1}{\eta K/\left(2\sqrt{{C}_{D0}K}\right)}-\sqrt{{\left(\frac{1}{\eta K/\left(2\sqrt{{C}_{D0}K}\right)}\right)}^{2}-\frac{4{C}_{D0}}{K}}\right\}$ (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
$h$ is fixed and the lift coefficient
${\stackrel{\xaf}{C}}_{L}$ is computed from Equation (23) for a fixed value of
$\eta $ . The QEG solutions are obtained using a Newton-Raphson method that typically needs 6 - 7 iterations to determine the
${\stackrel{\xaf}{q}}_{\text{QEG}}$ and
${\gamma}_{\text{QEG}}$ values. 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,
${x}_{\text{ALI}}$ . 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
${x}_{\text{ALI}}$ 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.

${x}_{\text{ALI}}={10}^{4}\left(\frac{9.69{\eta}^{3}+11.19{\eta}^{2}-6.526\eta -14.93}{{\eta}^{5}-1.086{\eta}^{4}-4.926{\eta}^{3}+6.52{\eta}^{2}-7.216\eta +5.819}\right)\text{\hspace{0.17em}}\text{ft}$ (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, h_{F}, and the downtrack location of the start of flare, s_{F}. Both are functions of glide-efficiency factor [6] . As with
${x}_{\text{ALI}}$ , 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.

${h}_{F}={10}^{3}\left(\frac{{P}_{1}{\eta}^{4}+{P}_{2}{\eta}^{3}+{P}_{3}{\eta}^{2}+{P}_{4}\eta +{P}_{5}}{{\eta}^{5}+{q}_{1}{\eta}^{4}+{q}_{2}{\eta}^{3}+{q}_{3}{\eta}^{2}+{q}_{4}\eta +{q}_{5}}\right)\text{\hspace{0.17em}}\text{ft}$ (25)

${s}_{F}={10}^{4}\left(\frac{{P}_{1}{\eta}^{2}+{P}_{2}\eta +{P}_{3}}{{\eta}^{5}+{q}_{1}{\eta}^{4}+{q}_{2}{\eta}^{3}+{q}_{3}{\eta}^{2}+{q}_{4}\eta +{q}_{5}}\right)\text{\hspace{0.17em}}\text{ft}$ (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 ${V}_{\text{TD}}=280\text{\hspace{0.17em}}\text{ft}/\text{s}$ , ${\stackrel{\dot{}}{h}}_{\text{TD}}=-5\text{\hspace{0.17em}}\text{ft}/\text{s}$ , and ${\gamma}_{\text{TD}}={\mathrm{sin}}^{-1}\left({\stackrel{\dot{}}{h}}_{\text{TD}}/{V}_{\text{TD}}\right)$ $=-1.023\text{\hspace{0.17em}}\text{deg}$ .

Figure 4. Start of flare location vs. glide efficiency factor.

A fourth-order polynomial altitude profile defines the flare maneuver phase

${h}_{f}={a}_{0}+{a}_{1}x+{a}_{2}{x}^{2}+{a}_{3}{x}^{4}$ (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:

$\frac{\text{d}{h}_{f}}{\text{d}x}=\mathrm{tan}{\gamma}_{f}={a}_{1}+2{a}_{2}x+4{a}_{3}{x}^{3}$ (28)

Fourth-order coefficients ${a}_{\text{0}}$ and ${a}_{1}$ are computed from Equation (27) and Equation (28) and the RLV’s altitude and flight-path angle at the start of flare phase where $x={s}_{F}+1400\text{\hspace{0.17em}}\text{ft}$ , ${h}_{f}={h}_{F}$ and ${\gamma}_{f}={\gamma}_{\text{QEG}}$ (note that ${\gamma}_{\text{QEG}}$ is the flight-path angle at the end of the steep phase). The second-and third-order polynomial coefficients ${a}_{2}$ and ${a}_{3}$ are also computed from Equation (27) and Equation (28) and the RLV’s altitude and flight-path angle at the touchdown point where $x=1400\text{\hspace{0.17em}}\text{ft}$ , ${h}_{f}=0$ , and ${\gamma}_{f}={\gamma}_{\text{TD}}$ .

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:

${\stackrel{\xaf}{q}}_{f}={c}_{0}+{c}_{1}x+{c}_{2}{x}^{2}$ (29)

The derivative of dynamic pressure with respect downtrack is:

$\frac{\text{d}{\stackrel{\xaf}{q}}_{f}}{\text{d}x}={c}_{1}+2{c}_{2}x$ (30)

The quadratic coefficients ${c}_{o}$ and ${c}_{1}$ are computed from Equation (29) and Equation (30) and the RLV’s dynamic pressure and its derivative at the start of flare phase where $x={s}_{F}+1400\text{\hspace{0.17em}}\text{ft}$ , ${\stackrel{\xaf}{q}}_{f}={\stackrel{\xaf}{q}}_{\text{QEG}}$ and $\left(\text{d}{\stackrel{\xaf}{q}}_{f}/\text{d}x\right)={\stackrel{\dot{}}{\stackrel{\xaf}{q}}}_{f}={\stackrel{\dot{}}{\stackrel{\xaf}{q}}}_{\text{QEG}}$ (note that ${\stackrel{\xaf}{q}}_{\text{QEG}}$ is constant during the steep phase and hence ${\stackrel{\dot{}}{\stackrel{\xaf}{q}}}_{\text{QEG}}=0$ . The remaining coefficient ${c}_{2}$ is computed from Equation (29) at the touchdown point where $x={x}_{\text{TD}}=1400\text{\hspace{0.17em}}\text{ft}$ and ${\stackrel{\xaf}{q}}_{f}={\stackrel{\xaf}{q}}_{\text{TD}}=0.5{\rho}_{\text{SSL}}{V}_{\text{TD}}^{2}$ .

Where the standard sea level air density ${\rho}_{\text{SSL}}=0.002377\text{\hspace{0.17em}}\text{slug}/\text{f}{\text{t}}^{3}$

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:

$\stackrel{\dot{}}{\gamma}=\frac{\text{d}{\gamma}_{f}}{\text{d}x}\cdot \frac{\text{d}x}{\text{d}t}=\frac{\stackrel{\xaf}{q}S{C}_{L}}{mV}-\frac{g}{V}\mathrm{cos}\gamma $ (31)

Substituting Equation (4) into Equation (31) and solving for lift coefficient, we obtain the open-loop lift coefficient:

${\stackrel{\xaf}{C}}_{L}=\frac{\mathrm{cos}\gamma}{\stackrel{\xaf}{q}\left(S/m\right)}\left({V}^{2}\frac{\text{d}{\gamma}_{f}}{\text{d}x}+g\right)$ (32)

where $V=\sqrt{2\stackrel{\xaf}{q}/\rho}$ . The downtrack derivative of the flight-path angle $\text{d}{\gamma}_{f}/\text{d}x$ can be determined from Equation (28)

$\frac{\text{d}{\gamma}_{f}}{\text{d}x}=\frac{12{a}_{3}{x}^{2}+2{a}_{2}}{{\left(4{a}_{3}{x}^{3}+2{a}_{2}x+{a}_{1}\right)}^{2}+1}$ (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 $\stackrel{\xaf}{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

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 ${h}_{f}$ during the shallow phase. During the QEG and flare phases, the lift coefficient is determined using the closed-loop command:

${C}_{L}={\stackrel{\xaf}{C}}_{L}+\delta {C}_{L}$ (34)

During the steep phase, the closed-loop lift coefficient ${C}_{L}$ 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 ${\stackrel{\xaf}{C}}_{L}$ is determined using Equation (23) and the feedback term is

$\delta {C}_{L}=-{K}_{\gamma}\left(\gamma -{\gamma}_{\text{QEG}}\right)-{K}_{\stackrel{\xaf}{q}}\left(\stackrel{\xaf}{q}-{\stackrel{\xaf}{q}}_{\text{QEG}}\right)$ (35)

The closed-loop control ${C}_{L}$ guides the RLV so that it tracks the reference dynamic pressure and flight-path angle. The feedback gains ${K}_{\gamma}$ and ${K}_{\stackrel{\xaf}{q}}$ 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 ${\stackrel{\xaf}{C}}_{L}$ is evaluated using Equation (32); however the closed-loop feedback term can be formulated by:

$\delta {C}_{L}=-{K}_{H}\left(h-{h}_{\text{ref}}\right)-{K}_{HD}\left(\stackrel{\dot{}}{h}-{\stackrel{\dot{}}{h}}_{\text{ref}}\right)$ (36)

where the feedback gains ${K}_{H}$ and ${K}_{HD}$ 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 $\eta $ 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 ${x}_{\text{ALI}}$ , flare altitude ${h}_{F}$ , and downtrack location of the start of flare ${s}_{F}$ ) 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 $\eta =0.835$ , $\stackrel{\xaf}{q}=216.98\text{\hspace{0.17em}}\text{psf}$ , $\gamma =-12.5\text{\hspace{0.17em}}\text{deg}$ , and $h=10000\text{\hspace{0.17em}}\text{ft}$ .

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 $\eta =0.835$ 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 ( $h=0$ and $x=1400\text{\hspace{0.17em}}\text{ft}$ ).

Figure 7(c) presents altitude versus downrange for several initial ranges, starting from ( ${x}_{\text{ALI}}=-4.987\times {10}^{4}\text{\hspace{0.17em}}\text{ft}$ ) corresponded to flight with $\eta =0.74$ to ( ${x}_{\text{ALI}}=-5.281\times {10}^{4}\text{\hspace{0.17em}}\text{ft}$ ) corresponded to flight with $\eta =0.93$ . 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 ${C}_{L}$ at runway touchdown matches well with the reference lift coefficient ${\stackrel{\xaf}{C}}_{L}$ 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.

References

[1] Hanson, J. (2000) Advanced Guidance and Control Project for Reusable Launch Vehicles. AIAA Guidance, Navigation, and Control Conference and Exhibit, 1-10.

https://doi.org/10.2514/6.2000-3957

[2] Hanson, J.M. (2002) A Plan for Advanced Guidance and Control Technology for 2nd Generation Reusable Launch Vehicles. AIAA Paper, 4557, 5-8.
https://doi.org/10.2514/6.2002-4557

[3] Hanson, J.M. and Jones, R.E. (2004) Test Results for Entry Guidance Methods for Space Vehicles. Journal of Guidance Control and Dynamics, 27, 960-966.

https://doi.org/10.2514/1.10886

[4] Moore, T.E. (1991) Space Shuttle Entry Terminal Area Energy Management. NASA Technical Memorandum 104744, November 1991.

[5] Tsikalas, G.M. (1982) Space Shuttle Autoland Design. AIAA Paper 82-1604, AIAA Guidance and Control Conference, San Diego, 9-11 August 1982.

[6] Kluever, C.A. and Neal, D.A. (2015) Approach and Landing Range Guidance for an Unpowered Reusable Launch Vehicle. Journal of Guidance, Control, and Dynamics, 38, 2057-2066.

https://doi.org/10.2514/1.G000909

[7] Kluever, C.A. (2007) Unpowered Approach and Landing Guidance with Normal Acceleration Limitations. Journal of Guidance Control and Dynamics, 30, 882.

https://doi.org/10.2514/1.28081

[8] Kluever, C.A. (2004) Unpowered Approach and Landing Guidance Using Trajectory Planning. Journal of Guidance Control and Dynamics, 27, 967-974.
https://doi.org/10.2514/1.7877

[9] Harl, N. and Balakrishnan, S.N. (2010) Reentry Terminal Guidance through Sliding Mode Control. Journal of guidance, control, and dynamics, 33, 186.
https://doi.org/10.2514/1.42654

[10] Zhao, Y., Sheng, Y. and Liu, X. (2014) Unpowered Landing Guidance with Large Initial Condition Errors. Guidance, Navigation and Control Conference (CGNCC), 2014 IEEE Chinese, IEEE, 1862-1867.
https://doi.org/10.1109/CGNCC.2014.7007465

[11] Schierman, J.D., Ward, D.G., Monaco, J.F. and Hull, J. (2001) A Reconfigurable Guidance Approach for Reusable Launch Vehicles. Proceedings of the AIAA, Guidance, Navigation and Control Conference and Exhibit, Montreal, 6-9 August 2001.

https://doi.org/10.21236/ADA436263

[12] Schierman, J., Hull, J.R. and Ward, D. (2003) On-Line Trajectory Command Reshaping for Reusable Launch Vehicles. Proceedings of the 2003 AIAA Guidance, Navigation, and Control Conference, Austin, 11-14 August 2003.

[13] Schierman, J.D., Ward, D.G., Hull, J.R., Gandhi, N., Oppenheimer, M.W. and Doman, D.B. (2004) Integrated Adaptive Guidance and Control for Re-Entry Vehicles with Flight-Test Results. Journal of Guidance Control and Dynamics, 27, 975-988.
https://doi.org/10.2514/1.10344

[14] Trent, A.D., Doman, D.B. and Iyer, R.V. (2007) Trajectory Planning for a Reentry Vehicle under Failure Conditions. American Control Conference, ACC’07, New York, 9-13 July 2007, 3868-3873. https://doi.org/10.1109/ACC.2007.4282306

[15] Heydari, A. and Balakrishnan, S.N. (2013) Path Planning Using a Novel Finite Horizon Suboptimal Controller. Journal of Guidance, Control, and Dynamics, 36, 1210-1214. https://doi.org/10.2514/1.59127

[16] Kluever, A. and AL Bakri, F. (2017) Semi-Analytical Terminal Trajectory-Planer for Unpowered Reusable Launch Vehicle. Journal of Guidance Control and Dynamics. (Under Review)