Abstract

Landing control of unmanned aerial vehicles (UAVs) is challenging because of the strong nonlinear dynamics, multivariable, model uncertainties, wind variations, and sensor noise. Motivated by this fact, this paper investigates an automatic landing system (ALS) that includes trajectory generation and guidance law for the first flight test of a turbine-based combined cycle technology demonstrator. Specifically, the control scheme increases the original model’s order to generate a reasonable monotone-decreasing throttle reference flare trajectory by the pseudospectral method. Subsequently, the guidance law based on innovative multivariable active disturbance rejection control is designed to robustly track the reference altitude and velocity simultaneously with high accuracy. The multivariable extended state observer (ESO) incorporated decoupling algorithm enhances the estimation capability and accuracy of potential problem in cross-coupling dynamics compared to the traditional ESO. It is proven that the closed-loop error dynamic has bounded-input bounded-output stability and an explicit upper bound is given. Numerical simulation verifies that the presented approach has better robustness and higher tracking accuracy for external disturbances and parametric uncertainties than the existing benchmark autolanding controller. Finally, flight tests show that the proposed ALS can land the vehicle effectively and safely under severe wind conditions.

1. Introduction

Despite the changing types of aircraft, fixed-wing aircraft still play the most important role both in civilian and military fields due to their higher flight speeds, larger payloads, and longer ranges [1]. With the widespread use of fixed-wing unmanned aerial vehicles (UAVs), their control problems have attracted the interest of researchers [2, 3]. High-speed UAVs have wider flight envelopes and higher dynamics than low-speed ones. Therefore, the former needs a more robust and better performance controller under finite actuator bandwidth. The automatic landing system (ALS) of UAVs is both a theoretical and practical problem because it is summarized as a multivariable nonlinear system with model uncertainty, external disturbances, and strong coupling [4]. This paper is aimed at developing an ALS for a high-speed fixed-wing turbine-based combined cycle (TBCC) technology demonstrator’s first flight test with robust and high tracking accuracy and fast soft landings.

Unlike manned aerial vehicles, most UAVs are only equipped with differential global position systems for current position information but do not have instrument landing systems. Therefore, it is necessary to design a nominal trajectory to guide the vehicle to land. Trajectory generation is essentially an optimal control problem with multiple constraints and has been widely studied in reusable launch vehicle landing control [5, 6]. However, UAVs are powered until the engine is turned off, which undoubtedly presents planning difficulties. Most literature prescribes the exponential reference height command using desired landing geometry [7, 8] or uses a few filters to soften the command [9]. Ignoring the speed variation may lead to too much speed and too little speed or even stall in the presence of wind shear and turbulence. There are many cases of landing accidents precisely due to unreasonable reference trajectories, resulting in UAVs with more significant thrusts touching the ground. Therefore, the trajectory of the flare phase is optimized based on the previous research work by combining landing geometry and the pseudospectral method [10].

In addition to traditional proportional integral derivative (PID) controllers [8], many advanced control algorithms have been investigated to guide and control UAVs to track the reference trajectory and touch down at the desired point. Dynamic inversion is adapted to cope with lightweight UAV landing and strong coupling, with assumptions known by an accurate model [11]. The discrete quantitative feedback theory (QFT) controller has been compared with a baseline PID controller for medium-sized UAV automatic landing [12]. However, flight results have been presented using QFT, where a significant number of controls need to be redesigned to achieve adequate experimental performance [13]. Lungu [14] designed an ALS controller that combines backstepping and dynamic inversion for a flying wing UAV subject to wind gusts and measurement sensor errors. Cui et al. [15] proposed an antidelay model predictive control scheme for carrier landing based on the symplectic pseudospectral method. In addition to the above control methods, ALSs also involve intelligent control methods, such as recurrent neural networks [16], deep learning-based fault estimation [17], fuzzy logic [18], swarm intelligence [19], and genetic algorithms [20]. Many of the aforementioned control methods are limited by the requirement of an accurate plant mathematical model, but such a requirement cannot be fulfilled in engineering practice. The dynamics of an automatic landing model are highly uncertain and suffer from substantially variable atmospheric disturbances. Therefore, the guidance law based on active disturbance rejection control (ADRC), requiring very little information about plant dynamics, is constructed to improve the antidisturbance capability [2123].

ADRC is a technique that emerged in the 1990s for estimating and compensating uncertainties [24]. This technique also has the potential to address nonlinear cross-coupling in multivariable systems [25]. The vehicle longitudinal guidance system is typically multiple-input, multiple-output (MIMO), where altitude (or flight path angle) and airspeed are the primary control variables; throttle position and elevator deflection are used as two-channel inputs. Several papers have mentioned the application of ADRC to multivariable control systems, but they were considered underdriven models [26], facing coupling control problems by feedforward compensation with a known disturbance [27], or not accounting for special nonlinear coupling effects [21]. To address these challenges, this paper proposes a novel guidance law for ALS. The method generalizes classical ADRC for nonlinear cross-coupling dynamics and provides a general tool for the high-speed fixed-wing UAV autolanding problem. The main contributions are as follows. (1)Compared to the most existing autolanding control method, which only considers the SISO dynamic, a decoupling control approach based on a multivariable ADRC (MADRC) is proposed. The decoupling problem is reformulated as disturbance suppression, where the disturbance is defined as cross-channel interference under the assumption of less model information. MADRC solves the nonlinear cross-channel effects using predetermined input-output data in real time. In the MADRC framework, an extended state observer (MESO) with enhanced estimation and compensation capability cancels the total disturbance in the combined control law, which improves disturbance suppression and tracking accuracy(2)Rigorous theoretical analysis and flight test validation are carried out for the proposed guidance law. The upper bound of the closed-loop system error is also given. Some modifications are proposed to improve the robustness against model uncertainties and disturbances from the implementation point of view. The measured output variables are chosen as the controlled variables instead of the estimated output values in conventional ADRC to compensate for phase lag. Because of control saturation, the amplitude of the extended states should be limited to properly allocate the control source to each component in the feedback control law(3)The flare phase nominal trajectory is designed using the pseudospectral method, which is a numerical solution to optimal control problems. Monotonically decreasing throttle is achieved by increasing the order of the original model. The boundary conditions and path constraints are specified in detail and the cost function accounts for the flight trajectory smoothness and the controlled variable flatness. We use the ground distance to interpolate the altitude and airspeed as commands for the guidance system, which significantly eliminates the effect of ground velocity on touchdown point accuracy

The remainder of this article is structured as follows: in Section 2, we formulate the trajectory generation problem and plan a reasonable landing trajectory offline with a gradually decreasing throttle using the pseudospectral method. In Section 3, a decoupling trajectory tracking controller based on MADRC is introduced, and the stability analysis is provided. Section 4 presents the simulations, robustness analysis, and flight test results. Finally, concluding remarks and plans for future work are given in Section 5.

2. Trajectory Planning

The normal wheeled landing process consists of three stages: approach, glide slope, and flare. A schematic diagram of an autonomous landing is shown in Figure 1, where is the runway altitude, is the flare altitude, is the flare airspeed, is the flare ground distance, is the initial approach altitude, is the initial approach airspeed, is the initial approach ground distance, and is the flight path angle in glide slope phase. is the geodetic coordinate frame (north east down or NED), established with its origin at the theoretical touchdown so that for an approaching UAV.

2.1. Vehicle Model

The high-speed fixed-wing wheeled UAV studied in this paper is a wing body fusion design with a dual vertical tail, a midswept wing with slats, and a retractable landing gear, as illustrated in Figure 2. The parameters are shown in Table 1. is the body-fixed frame attached to the center of gravity (cg) of the UAV.

The UAV is considered a mass point, and its landing motion in the vertical plane is defined by the following ordinary differential equations with respect to a flat-Earth model [28]. where is the airspeed, is the flight path angle, is the down-track position along the runway centerline, is the altitude, and is the gravitational acceleration of Earth. represents the turbojet engine thrust, depending on throttle position , altitude , and airspeed . and denote lift and drag forces, respectively, which are defined in the usual manner: where is the dynamic pressure. The atmospheric density is computed using the U.S. 1976 Standard Atmosphere. The lift force coefficient and drag force coefficient are calculated using three-dimensional lookup tables with the angle of attack , Mach number , and elevator deflection as the independent variables.

The phugoid mode described in Equations (1a), (1b), (1c), and (1d) involves changes in the flight path and is a relatively slow mode. It also involves the translational degrees of freedom and is dependent on thrust, lift force, and drag force, as well as their variation with dynamic pressure.

2.2. Glide Slope Trajectory Design

This paper focuses on the glide slope and flare stages. Because of runway length limits and routine operation over civilian airspace, the glide phase is defined as beginning at  m and  m/s, and the flight path angle at the glide interface is set to -4 deg. Therefore, the reference altitude and reference velocity of the glide section are functions of the ground track distance:

In general, the velocity command is constant [29, 30]. Unexpectedly, the throttle will increase when the flight speed falls below the command, even when close to the ground. To avoid the common problem, a decreased glide phase airspeed is chosen to ensure the desired velocity at touchdown. In addition, the UAV attitude angle gradually increases to further reduce airspeed in the flare phase. Therefore, the pseudospectral method is chosen to generate a robust flight trajectory with a monotone-decreasing throttle.

2.3. Flare Trajectory Design Based on the Pseudospectral Method

The flare phase is activated at a predefined altitude . The traditional longitudinal trajectory for the flare phase is a trajectory with a smaller glide slope angle or an exponentially flattened trajectory, while the lateral trajectory is a straight line aligned with the runway centerline. By contrast, a landing trajectory planning problem with decreased throttle using the Radau pseudospectral method (RPM) is constructed as follows: the RPM is an orthogonal collocation method where the collocation points are the Legendre-Gauss-Radau points [31, 32]. Some of the interesting features of the RPM are as follows: (1) it is a Gauss quadrature implicit integration scheme; (2) the algorithm has exponentially fast convergence for problems where the solution is smooth; and (3) there is an elegant connection between the continuous-time optimal control problem and the discrete approximation. Increase the order of the original model by adding an equation to Equations (1a), (1b), (1c), and (1d).

Hence, and are treated as input variables. Minimize the cost function subject to the dynamic constraints in Equations (1a), (1b), (1c), (1d), and (4); the boundary conditions are as follows: and the inequality path constraints are as follows: where , , and are the weight coefficients. The first equation in Equation (7) constrains a monotone-decreasing throttle in the flare phase. The flare trajectory should be gentle and smooth in terms of control stability. The variation in energy and flight path angle should be as small as possible. In addition, it is necessary to minimize the amount of throttle position. Therefore, the cost function is formulated as Equation (5). It is critical to select reasonable weights based on mission requirements and flight conditions.

The values of boundary conditions and path constraints are explained as follows: the flare speed must be less than the maximum taxiing ground speed of 65 m/s allowed by the strength of the landing gear. The angle of attack at stall is 35 deg, and the tail-rub angle is 13.5 deg, so the tail-rub angle determines the lower limit of . The minimum speed for a mass of 430 kg at  deg is about 44 m/s. The normal touchdown speed is generally 1.1–1.2 times the minimum speed, or 50 m/s, to provide a safety margin. In a comprehensive analysis, the ideal touchdown speed is chosen to be 55 m/s.

From Equation (1d), the flight path angle is determined by the vertical sink rate and airspeed together. As mentioned earlier, the flight path angle in the glide phase is -4 deg. The flight path angle should be much larger than that in the glide phase to reduce the descent rate. According to test results, the maximum vertical sink rate that the landing gear can withstand is -1.5 m/s with a landing mass of 430 kg. A manned aircraft typically hits the ground with a vertical sink rate of -1.5 to -0.5 m/s, so the ideal flight path angle should be greater than -0.5 deg.

When the touchdown ground speed is 55 m/s and the average acceleration during taxiing is -2 m/s2, the taxi length is 756 m. Since the total length of the runway is 2500 m, the theoretical landing point is designed near the midpoint of the runway. It is a good practice to stay above the runway during the initial flare phase to prevent a touchdown in advance. Based on the direction of defined above, the maximum permissible ground trajectory position of the initial flare is  m and the touchdown position is  m. The preset flare altitude range is defined as 20–25 m above the runway, and the runway altitude is 1000 m.

The pitch angle is limited to 2–13.5 deg when touching the ground to avoid nose landing gear damage or the tail scraping the ground. The pitch angle is calculated approximately according to

According to trim results, when the flight path angle is -4 deg, the throttle position is 15% to 40%. Therefore, the lower and upper limits on the control inputs during the flare phase are defined in Equation (7).

2.4. Trajectory Planning Results

The initial flare height  m, the flare phase time  s, and the ground track distance m. The planning results of the state and control variables are shown in Figures 3 and 4. In the beginning, the altitude drops rapidly, and after 4 s, it decreases nearly linearly and slowly with a sink rate greater than -2.0 m/s. The velocity decreases approximately linearly and smoothly with the ground track distance.

The vehicle is pulled up smoothly from the initial vertical sink rate of -4.5 m/s, converges rapidly to within -1.0 m/s, and finally touches down at -0.5 m/s. The flight path angle changes in sync with the vertical sink rate. The pitch angle gradually increases to ensure that the attitude is effectively pulled up and the UAV lands in a two-wheel mode. The touchdown pitch angle is 10.1 deg, which has a large safety margin compared to the tail-rub angle. The throttle position decreases with velocity and is generally within the safe range as shown in Figure 4. In practice, the reference altitude and velocity trajectory in the flare phase are interpolation functions of the ground track distance : where and have been described in Figure 3. The initial height of the flare phase is used to determine whether Equation (3) or (9) should be implemented. The method generates a trajectory that satisfies all constraints, including a monotone-decreasing throttle position and an exponentially decreasing altitude profile.

3. Guidance Law

The control scheme for an automatic landing generally consists of two components. The first component is the trajectory, or the desired state during landing, which depends on the landing specification. The second component is reference trajectory tracking or how the vehicle follows the landing trajectory. While both components are intricately linked, each is discussed separately here. The guidance algorithm adopts a MADRC controller for the multivariable system to decouple the altitude and velocity channels while simultaneously maintaining the tracking performance. The throttle position and the angle of attack are obtained according to the reference velocity, sink rate, and their feedback variables, as shown in Figure 5.

3.1. Concept of Multivariable Extended State Observer

MESO is the core of the MADRC as shown in Figure 6. The MESO estimates the state variables and system disturbance, which is considered as an extended state variable including both internal dynamic uncertainties and external disturbances.

Consider the multivariable affine nonlinear system of m-inputs and m-outputs: where , , and are the state, input, and output of the system, respectively, while and are the input and output matrices, respectively. are functions of the state variables, time, and the external disturbance .

Suppose that is the best available estimate of . The extended state is defined by

It is treated as the total disturbance, including the external disturbance, the modeled dynamics, and the unmodeled dynamics, and then yields

where is the change rate vector of the extended state. It is assumed to be unknown but bounded.

Assumption 1. is an -order unity matrix, so the system given in Equation (12) is fully observable.
The observer for the system is described by where denotes the state estimation and denotes the MESO gain. , , and are the estimations of the original state , the extended state , and the output , respectively. where is the unity matrix. It is noted that the presented MESO is a linear system whether the original system is linear or nonlinear.

Let be the MESO estimation errors and be the system input. The error dynamic can be written as where

In particular, consider a special case where the observer gains are chosen as [33] where is defined as the MESO bandwidth. Hence, is a Hurwitz matrix, and the characteristic polynomial is written as follows:

Assumption 2. The external disturbances and their derivatives are bounded, and for , all partial derivatives of are bounded over .

Theorem 3. Suppose that Assumption 2 is satisfied, then there exists positive constant , such that .

Proof. Solving Equation (15) for , we have where is a Hurwitz matrix. According to Equation (18), the eigenvalues of matrix can be expressed in terms of multiple identical poles: . Sorting the eigenvalues by magnitude, , there exists an invertible real matrix, , that transforms into its Jordan form, yielding where is a Jordan block associated with the eigenvalue of . A Jordan block of order takes the form Therefore, Since is the maximum eigenvalue, the Frobenius norm of is given by Because of , the term is a monotonic decreasing function. For , .
The first derivative of is and only when . The second derivative of is negative; therefore, has a maximum extremum Note that for , the Frobenius norm of has supremum which implies that there exists such that According to the matrix norm inequality, and substituting Equation (28) yields

From Equation (30), if the derivative of the extended state is bounded; then, the estimation errors of the MESO are convergent. For Assumption 2, airspeed control is applied to explain its meaning. This means that the rate of change in acceleration is finite. The variation of the force acting on the vehicle is indeed bounded. More generally, it implies that there is a limited rate of change in the physical world or that no change occurs instantaneously.

Remark 4. Equation (17) parameterizes the observer bandwidth, reducing the number of parameters and simplifying the design process. Now, the tuning of becomes much simpler and more intuitive since the have explicit physical meanings. The observer errors can be as small as expected when is tuned to be large, according to Equation (30). In practice, the values should be tuned according to the variation of the extended state. Generally, the faster the extended state changes, the larger the value should be tuned. However, there are some cases where the observer bandwidth is limited by the noise frequency and the controller sampling frequency, leading to oscillations in the estimates. Therefore, the tuning of is a process of a tradeoff between performance and robustness.

3.2. Guidance Law Design

The use of the chain derivative rule for a compound function to take the time derivative of Equation (1d) yields

Substitution of the dynamic equations for flight path angle, Equation (1b); airspeed, Equation (1a); and pitch angle, Equation (8) into Equation (31) yields, after some algebra, the required reference vertical sink rate:

Equations (1a) and (32) constitute the dynamic model:

Let , , and be the system state, input, and output, respectively. Then, the dynamic model shown in Equation (33) is a two-input, two-output system, which has a similar form to Equation (10). The output matrix is a 2nd-order identity matrix satisfying Assumption 1, and the control matrix is a 2 by 2 square matrix.

Remark 5. is the best available estimation of matrix , which is obtained from known conditions. The elements of matrix have explicit physical meanings as follows: where subscript denotes the equilibrium point. and are the trim results of the lift coefficient and drag coefficient. is the partial derivative of thrust with respect to throttle, which can be obtained from the ground test engine test or user manuals. , namely, the lift-curve slope, determines how changes in due to turbulence translate into changes in lift. represents the slope of the drag coefficient versus . The difference between and is estimated and compensated as part of the extended state, including thrust slightly tilted with respect to the body -axis, aerodynamic data error, changes in aircraft mass properties, and ground effects. Since the matrix is defined artificially, it is guaranteed to be invertible.

The multivariable extended state dynamic equations are formulated from Equation (13):

The controller based on MADRC is illustrated in Figure 7. It is designed to drive the system output, defined in Equation (33), and track the reference signal, . The controller is then formulated as

In Equation (36), are the tunable control gains and are used for estimating and compensating the extended state.

Remark 6. The state variables are used in the feedback control law instead of the estimation used in conventional papers. One reason is that the state variables passing through MESO are equivalent to a low-pass filter, contributing to a certain phase delay. In addition, the sensor signals are redundant with low-measure noise and drift and have been used on other UAVs of the same class with the expected results. The airspeed is measured by the air data computer (ADC) with a sampling frequency of 50 Hz. Once the sensor fails, we must switch to ground speed, which is only 10 Hz from the real-time kinematic data in place of airspeed. Furthermore, the radar altimeter (RA), whose signal is filtered after differencing to accurately measure the sink rate , is single redundant. It is backed up by the integrated navigation system (INS), which measures the altitude relative to a theoretical reference while the actual runway might be curved or sloped.

3.3. Closed-Loop Stability Analysis

Assumption 7. The reference signal and its derivative are bounded.

Theorem 8. Suppose that Assumption 7 is satisfied, then there exist positive constants and such that the tracking error of the closed-loop system is bounded.

Proof. Define the tracking error as The dynamic equation of the error is described as , and substituting Equation (36) yields The solution is given by There exist and , such that is a Hurwitz matrix. Similar to the proof of Theorem 3, there exists such that the Frobenius norm of has supremum: where denotes the smallest eigenvalue of matrix and .

According to Theorem 3, the error estimate is bounded while Assumption 7 is satisfied. Thereby, the solution is guaranteed to be bounded.

3.4. Discussion

The is obtained from the guidance law in Equation (36). However, the pitch angle is treated as the controlled variable in practice. Equation (8) is an approximate relation between and , and it is easy to compute the pitch angle command from , because the pitch angle is derived from the INS and its measurement is less disturbed by wind disturbance than the aerodynamic angle, especially close to the ground. Noise from the sensor can make it challenging to achieve a precise landing. Hence, feedback is avoided, if possible, because of the difficulty of achieving an accurate, rapid responding, and noise-free measurement and because the sensor is vulnerable to mechanical damage. In addition, the controlled pitch angle needs to meet the landing specification to prevent a tail collision. The flight envelope of the vehicle covers a wide range of dynamic pressure . For example, the dynamic pressure is as low as 2.3 kPa during landing, but at Mach 0.6 at sea level, the dynamic pressure is 25.5 kPa. Large dynamic pressure changes can lead to large variations in the coefficients of the corresponding dynamic equations, as in Equations (1a), (1b), (1c), (1d), and (33). Other factors in Equation (34) also contribute to changes in vehicle dynamics. The dimensionless aerodynamic coefficients vary with Mach number and aerodynamic angles. This involves changing the measured amount of feedback as a function of a scheduling variable. Here, dynamic pressure is selected as the sole scheduling variable because it is the dominant term during autolanding.

4. Numerical Simulation and Flight Test Results

Trajectory generation and the guidance law for autolanding were discussed in the previous sections. The proposed method is implemented in MATLAB/Simulink using a full-state nonlinear numerical model. Moreover, a flight test based on the control scheme was conducted successfully in October 2019. Therefore, the superiority of the proposed method is fully demonstrated.

4.1. Implementation

The following dynamics are added to Equations (1a), (1b), (1c), and (1d); then, the full-state nonlinear longitudinal vehicle model is obtained. where is pitching moment coefficient and and are pitching moment dynamic derivatives. Also, and are lift force dynamic derivatives. denotes the actuator command,  rad/s denotes the natural frequency, and denotes the damping ratio. The other parameters have been shown in Table 1. In addition, the model consists of nonlinear components with a deflection limit of  deg and a rate limit of  deg/s. The throttle servo and engine response dynamic model is a simple-lag filter with a time constant of 0.2 s. The throttle position remains between idle (0%) and full power (100%), with a throttle rate of less than 20%/s.

The flight control system sample time is set to  s to represent sample times for real applications. A controller with integral action, like in Equation (36), that is combined with a saturated actuator may produce undesirable effects. The integrator may be an unstable system and the actuator remains saturated even if the process output changes, so the feedback path is broken. We avoid winding up the integrator by stopping the update of the integral when the actuator is saturated, so it is necessary to limit the amplitude of to allocate an amount of control to each component.

The pitch attitude control adopts a cascade ADRC controller based on pitch angle and pitch angle rate feedback, which is not discussed in this paper. As mentioned in the previous sections, the ADC, RA, and INS sampling times are 20 ms, 10 ms, and 10 ms, respectively. Except for the dynamic model introduced in Equations (1a), (1b), (1c), (1d), and (43), the simulation is delayed by 20 ms to include a simple representation of the accumulated time delay effects caused by both sensors and computations. In addition, there are limits on the controller in the simulation to more closely resemble the real physical system. When a state variable reaches a limit, a nonzero derivative is only allowed in the direction that takes the state variable out of the limit.

4.2. Comparison with Existing Design

We choose three existing benchmark methods to compare with the proposed control scheme. They are the (1) total energy control system (TECS), (2) ADRC, and (3) incremental nonlinear dynamic inverse (INDI). The design process of the baseline controller is briefly described here, and corresponding literature should be consulted for details.

4.2.1. Design 1. MADRC

The proposed method in Section 4 is designated as MADRC for the purpose of comparison. The MESO gains are selected based on the guideline that the observer bandwidth is larger than the controller bandwidth. Considering unmodeled dynamics, on the other hand, the observer bandwidth should be small enough to avoid stability problems in high frequency. The MESO gains are and . The baseline control gains are and .

4.2.2. Design 2. TECS

Lambregts and Faleiro and Lambregts [34, 35] proposed the TECS, where the total energy is the sum of kinetic energy and potential energy. Through throttle adjustment, energy is distributed by elevator deflection to achieve the control and distribution of energy. The basic regulation law is a PID controller with feedforward compensation.

4.2.3. Design 3. ADRC

The classical ADRC scheme does not have a special module to handle the coupling dynamics. It treats altitude and velocity as two separate SISO systems and uses two independent ESOs for estimation [21]. The controller gains are the same as in the MADRC scheme.

4.2.4. Design 4. INDI

INDI is a robust modification of NDI, which makes the controller less dependent on the model at the cost of sensitivity to sensor measurement delays and the use of measurement data that are not readily available. A polynomial prediction filter is presented in [36] to derive not-readily measured data. The baseline controller gain is consistent with the ADRC.

The chosen initial conditions were level flight at  m/s,  m, and  kg, with gear deployed for a nominal . The complete landing process curves for the four control methods are illustrated in Figure 8, including the approach, glide slope, and flare phase. Because of the accurate initialization of the controller, there are no transients in the level UAV flight. The aircraft maintains an airspeed of 80 m/s at 1325 m and is level with a trim throttle position of 55.8% and a trim of 5.1 deg. At the beginning of the descent, the airspeed increases in the ADRC and INDI approaches because of the conversion of potential energy into kinetic energy. The airspeed only tends to decrease when the excessive airspeed drives the throttle position reduction. The TECS, on the other hand, improves the transient response because of a throttle feedforward, which directly reduces the airspeed. The MADRC has a better transient response via its unique decoupling module, rather than feedforward compensation. It uses the MESO to feel the altitude drop in advance to generate a throttle reduction command. Owing to the tight control, the throttle position shows some rapid excursions that lead to possible rate limiting of the actuators. ADRC is the first to touch down, with a minimum sink rate of -0.25 m/s, but is still 332 m away from the designated point. INDI is the second to touch down, with a sink rate of -1.21 m/s, and is 215.7 m away from the ideal point. MADRC and ADRC have similar sink rates of -0.6 m/s and landing airspeeds of 57.7 m/s, but MADRC is closer to the designated point, within about 20 m.

4.3. Simulations with Wind Variation

The wind is the most significant factor affecting safety during the landing process. We divided the wind model into two parts that are simulated separately [37]. The first part is mean wind, discrete gusts, and wind shear. Near the ground, the horizontal component of the wind is about 5 m/s. And there are two cases in the first part: tailwinds and headwinds. The other simulated part is turbulence. Since this paper only covers the design in the longitudinal mode, only the wind components in the -axis and -axis directions are considered.

Table 2 indicates the performance requirements at touchdown used to verify the success probability of the landing guidance and control law. We focus on four metrics: landing sink rate, ground distance, pitch angle, and ground speed. The ground distance indicates the deviation from the designated point, and a positive value indicates a touchdown past the ideal touchdown point. Once the results exceed one of the given metrics ranges, the landing is considered a damaging landing.

Table 3 summarizes the simulation results of the four design methods for different wind conditions. We pay the most attention to the state variations at the time of touchdown, so the altitude and airspeed profiles near the flare are plotted in Figure 9. The pentagram indicates the planned start point of the flare. MADRC is the best performing method with soft landing for all three deviation states. It is clear that the MADRC is more concentrated at the beginning of the flare in terms of altitude and airspeed. During the flare phase, there is almost no effect of tailwind or headwind on the MADRC altitude profile. The TECS performance is the worst among the four methods due to its inaccurate feedforward compensation. The landing is so early under tailwinds that the ground speed is more than 72.6 m/s. The other nonnominal cases are also hard landings caused by excessive sink rates. Since the sink rate exceeds the limit, the ADRC is a hard landing in the turbulent case. This indicates that the independent ESO is unable to cope with the fast time-varying turbulence, which introduces stochastic angular rate noise. When landing in headwinds, INDI fails to slow down the speed because of poor airspeed control, despite with the same control parameters as ADRC, leading to damage. Also similar to ADRC, the turbulent condition is a hard landing in the INDI method. Therefore, these results demonstrate that MADRC has better antidisturbance performance in the presence of wind variation.

4.4. Monte Carlo Simulation

Monte Carlo simulations [38] containing 1000 autolanding trajectories were performed for each control scheme with random wind dispersions and uncertainties in drag coefficient, lift coefficient, atmospheric density, engine thrust, thrust-line skew and eccentricity, and vehicle mass. Turbulence is already a stochastic process, and the other wind components are modeled as Gaussian distributions with the local average wind speed at the airport as the mean value and a 3 of 20%. The local wind field is a function of the flight height, so it contains mean wind, wind shear, and gusts. Random biases of lift coefficient, drag coefficient, air density, engine thrust, and mass are also applied with Gaussian distributions having a zero mean and 3 of 10%, 20%, 5%, 10%, and 30 kg, respectively. It is noted that the worst combination of these parameters changes the aerodynamic acceleration by 15%. The thrust tilt and eccentricity generate force and moment disturbances that separately impact the guidance and attitude control systems. Their biases are defined as normally distributed mean values of -2 deg and 20 mm with a 3 of 10%.

The Monte Carlo simulation results are shown in Figure 10 and are summarized in Table 4. The model uncertainty simulations show that the MADRC controller has a 99.8% soft landing rate compared to 71.9% for TECS, 98.9% for ADRC, and 85.1% for INDI. The mean touchdown vertical sink rate for the MADRC controller is -0.61 m/s, with a standard deviation of 0.14 m/s. ADRC has similar statistical results. The mean value of for TECS and INDI is about -1.2 m/s, and their standard deviations are 0.59 m/s and 0.22 m/s, respectively. Among the four methods, the hard landings due to the landing sink rate exceeding -1.5 m/s are 2 for MADRC, 247 for TECS, 8 for ADRC, and 69 for INDI. Except for the negative control performance of TECS, the dispersion in touchdowns is small for the other three methods because the control commands are interpolated with respect to the ground trajectory distance, so the effect of ground speed is nearly eliminated. The attraction of the interpolating approach is evident in situations where the runway distance is limited, and the landing speed is relatively large. The statistical results of pitch angle all meet the specification between 2 and 13.5 deg. 80 trajectories of INDI are damaging landings, the majority of which are caused by exceeding ground speed, supported by the evidence of the wind variation simulations. However, many of the damaging landings of TECS are due to early touchdown, as illustrated in Figure 10(a).

Figures 11(a)11(d) each plot 1000 history curves of the Monte Carlo simulations for different controllers. It is noted that MADRC has the best control performance of airspeed, both in the early glide slope and in the flare phase. Moreover, the airspeed curve gathers around the command, even though wind deviation is added to the simulation. TECS also has positive control performance of airspeed in the flare section because of the feedforward compensation, but the standard deviation of increases due to early touchdown. ADRC has similar sink rate curves to MADRC, but the airspeed control is worse because the independent ESO cannot perceive the altitude drop in the throttle channel in advance. From the curves, the airspeed control of INDI is the worst, meaning that it is more sensitive to control gain and wind disturbance.

In summary, the evaluation results including wind variation and Monte Carlo simulations show that the MADRC controller provides significantly better tracking performance and robustness with plant uncertainty and disturbances when considering landing safety, sink rate, roll-out distance, and touchdown pitch attitude.

4.5. Flight Data

In addition to the simulation and verification completed in the preceding sections, we performed extensive hardware-in-the-loop semiphysical simulations with an actual flight control computer and servosystem actuators with loads. The vehicle, as shown in Figure 2, was designed to demonstrate a turbine-based combined cycle engine. The first flight test verified the vehicle’s takeoff and landing performance, with only the turbojet engine working. The test was conducted in October 2019 with a flight duration of 8 min and a maximum flight speed of Mach 0.24. The flight trajectory is shown in Figure 12. The planned five-sided route is a green curve, each side is 6 km long, and the yellow line represents the real flight path of the vehicle in the - plane. WP6 is both a theoretical touchdown site and a takeoff point. So the taxiing heading during takeoff and landing is opposite. And the pentagram represents the initial glide descent point. A range of 5.1 km is reserved for the glide slope and flare phase. The controller decoupled the altitude and velocity channels well during the landing process and tracked the ramp signal with simultaneous altitude and velocity changes, guaranteeing high tracking precision. The attitude during the flight was smooth and the controlled attitude angle was tracked with high accuracy. The touchdown point was precise, less than 10 m from the desired point. The landing ground speed was 58.8 m/s and the airspeed was 55.7 m/s. In addition, the main landing gear touched the ground in a comfortable position with a pitch angle of 8.8 deg. The glide and flare process achieved a successful and gentle soft landing, and the landing performance was excellent and attributed to the control scheme proposed in this paper.

Figure 13 depicts the complete ALS flight data of the first test. The flight data are sampled at 10 Hz and corrected by the time scale. In Figure 13(b), the dashed line refers to the ground speed and the solid line to the airspeed. It is noted that the ground speed meets the requirement, even for a tailwind landing. During the flare phase, the sink rate in the flight data exhibits small amplitude oscillations. It is inferred that this is due to a lag in pitch angle command tracking in Figure 13(d). The overshoot of the sink rate control resulted in a long glide time in the flare phase, with the touchdown point actually being only 5 m from the ideal point. This is significantly better than the numerical simulation. After some time, the sink rate control error converges to zero, as shown in the stability analysis. The terminal throttle does not monotone decrease as expected but acceptable oscillates slightly from 20% to 30%. The reason is that the tailwind near the ground causes a decrease in airspeed and an increase in throttle command to maintain the airspeed. The flight test, which agreed with the numerical simulation, verified the stability and feasibility of the proposed method against model inaccuracy and wind disturbance.

5. Conclusions

This paper proposes a longitudinal ALS consisting of a reference trajectory generator and a guidance law to meet the high accuracy and strong robustness requirements for the high-speed fixed-wing UAV first flight test. The Radau pseudospectral method is used to generate a flare trajectory with a monotone-decreasing throttle, and the guidance law uses MADRC as the baseline algorithm. In the MADRC control scheme, the control parameters have physical meanings, thereby making the design simple, intuitive, and efficient. Moreover, some of the implementation issues encountered during the design, such as sensor redundancy, feedback signal selection, and control gain tuning, were discussed and illustrated in detail. The Monte Carlo simulation gives a soft landing rate of 99.8%, compared to 71.9% for TECS, 98.9% for ADRC, and 85.1% for INDI, which establishes that the proposed ALS has high tracking accuracy, strong robustness, and effective decoupling capability on the altitude loop and velocity loop, even in the presence of severe wind variation and significant parameter uncertainties. The successful flight test demonstrates that the method has promising potential in advancing the state of the art in the engineering application of fixed-wing UAV autolanding.

Future work involves (1) eliminating underdamping phenomena occurring in the flight test to further improve landing performance, (2) extending to the lateral-directional to replace existing PID controller, and (3) investigating high-order MESO to enhance rejection capability to fast time-varying turbulence.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. U21B6003).