Combining Stable Inversion and H ∞ Synthesis for Trajectory Tracking and Disturbance Rejection Control of Civil Aircraft Autolanding

: The landing phase during a ﬂight probably is the most dangerous part, as most of the accidents occur in this phase. A robust trajectory tracking controller is presented to autoland a civil aircraft subjected to severe wind disturbances to improve the aircraft’s safety. Firstly, the dynamic models of the aircraft and windshear are built. Secondly, a stable inversion (SI) based robust autolanding controller (SIRAC) is proposed. In this architecture, the SI algorithm is used to improve the output tracking precision, while the H ∞ synthesis is applied for enhancing robust stability against uncertainties caused by wind disturbances. Finally, two scenario simulations are carried out for the automatic landing control of a large civil aircraft. Signiﬁcant performances on the system have been achieved without any disturbance. In addition to that, the proposed SIRAC can also track the desired autolanding trajectory with high precision, even under large wind condition.


Introduction
For the civil aircraft, the automatic control system plays an important role in assisting the pilot in all flight phases. Usually, there are three phases in one flight process, that is, takeoff, cruise, and landing. The landing phase is the most challenging phase, during which many variables need to be controlled simultaneously and high airworthiness must be met. The altitude and velocity are rather low during the aircraft landing stage. Hence, accidents are more likely to happen in this phase. BAE firstly developed the Automatic Landing System (ALS) for commercial aircraft in 1965 to increase the safety of the landing maneuver [1]. It has been widely used since then, as it can provide safe and comfortable landing. However, the ALS is far from being robust to strong environmental disturbances (e.g., windshear, turbulence, etc.).
Typically, there are four segments in the landing maneuver [2] (see Figure 1): alignment, glide slope, flare, and taxiing. The glide slope and flare segments are mainly studied in this paper, as they are the most challenging segments. During landing, the aircraft needs to track the desired trajectory precisely; this desired trajectory is provided by the Instrumental Landing System (ILS). The ALS works in harmony with the ILS, which is commonly available in a rapidly growing list of airports. The ILS contains two beam transmitters to guide the aircraft during the landing procedure. One is for localizer and the other one is for glide slope. An aircraft aligns and localizes itself to a befitting altitude and approach angle, according to these transmitters. Subsequently, it starts landing and decreasing the altitude [3]. Hence, the autolanding control is inherently a trajectory tracking control problem. Over the past two decades, numerous research results of autolanding control design have been developed to cope with the civil aircraft landing problem. Classical control methods are popular methods due to their simplicity and reliability. Sadat-Hoseini et al. [4] carried out a linear quadratic regulator (LQR) feedforward closed-loop controller, which was simulated on both the longitudinal and lateral-directional channels through different loops simultaneously. Magni et al. [5] proposed a control design approach that was based on eigenstructure assignment while using dynamic feedback, which had two advantages, first was that it could be regarded as an efficient controller order reduction, second was the design methodology that could be initialized by H ∞ or μ synthesis. With the development of digital flight computers, modern control theories have been progressively developed in an autolanding system. Wang et al. [6] developed a multivariable model reference adaptive control method that was implemented with state feedback to enable a safety landing, despite the parametric variations. An ALS was designed using sliding mode control (SLC) technique in [7]. Lyapunov stability criteria was used to force the sliding function to reach the solution and converge the aircraft's landing path to the desired trajectory. Juang et al. [8] also used the SLC integrated with the cerebellar model articulation controller to improve the ability of disturbance rejection of landing system, but the parameters of the SLC were adjusted by different methods, such as the genetic algorithm, particle swarm optimization, and chaotic particle swarm optimization. Fuzzy-logic dynamic inversion was used to suppress hard landing and roll oscillation in turbulent air condition to increase the adaption of the landing system to different environments [9]. On the other hand, Juang et al. [10] proposed a controller that was based on a multi-layered fuzzy neural network structure, which provided the control signals during landing procedure. The backpropagation algorithm was used to train the network. In the meantime, H ∞ control has been widely used in the ALS. Tamkaya et al. [11] combined the model following the method with the H ∞ synthesis to construct a dynamic controller, which could improve the performance of the conventional ALS, even under severe weather conditions. Theis et al. [12] presented a comprehensive autopilot for crosswind landing. The individual control loops were designed while using robust control methods and classical loopshaping, which provided a complete qualitative design strategy. A flare control law was exploited via multi-channel H ∞ synthesis [13]. The controller controlled the vertical speed of an aircraft while minimizing the influence of windshear and ground effects. Although the autonomous technologies can improve flight operations and overall aircraft performance, pilots will remain at the heart of operations in practice. Autonomous technologies are paramount to supporting pilots, thus enabling them to focus less on aircraft operation and more on strategic decision-making and mission management during landing. The energy approach to the flight control was used to access the current and predict the future states of an aircraft in order to improve the informational and situational awareness of the aircrew [14]. The landing maneuver was simulated with the ahead obstacle and the engine failure, which showed the correction of the algorithm. Another good project was Clean Sky 2 of European Clean Sky [15]. One Over the past two decades, numerous research results of autolanding control design have been developed to cope with the civil aircraft landing problem. Classical control methods are popular methods due to their simplicity and reliability. Sadat-Hoseini et al. [4] carried out a linear quadratic regulator (LQR) feedforward closed-loop controller, which was simulated on both the longitudinal and lateral-directional channels through different loops simultaneously. Magni et al. [5] proposed a control design approach that was based on eigenstructure assignment while using dynamic feedback, which had two advantages, first was that it could be regarded as an efficient controller order reduction, second was the design methodology that could be initialized by H ∞ or µ synthesis. With the development of digital flight computers, modern control theories have been progressively developed in an autolanding system. Wang et al. [6] developed a multivariable model reference adaptive control method that was implemented with state feedback to enable a safety landing, despite the parametric variations. An ALS was designed using sliding mode control (SLC) technique in [7]. Lyapunov stability criteria was used to force the sliding function to reach the solution and converge the aircraft's landing path to the desired trajectory. Juang et al. [8] also used the SLC integrated with the cerebellar model articulation controller to improve the ability of disturbance rejection of landing system, but the parameters of the SLC were adjusted by different methods, such as the genetic algorithm, particle swarm optimization, and chaotic particle swarm optimization. Fuzzy-logic dynamic inversion was used to suppress hard landing and roll oscillation in turbulent air condition to increase the adaption of the landing system to different environments [9]. On the other hand, Juang et al. [10] proposed a controller that was based on a multi-layered fuzzy neural network structure, which provided the control signals during landing procedure. The backpropagation algorithm was used to train the network. In the meantime, H ∞ control has been widely used in the ALS. Tamkaya et al. [11] combined the model following the method with the H ∞ synthesis to construct a dynamic controller, which could improve the performance of the conventional ALS, even under severe weather conditions. Theis et al. [12] presented a comprehensive autopilot for crosswind landing. The individual control loops were designed while using robust control methods and classical loopshaping, which provided a complete qualitative design strategy. A flare control law was exploited via multi-channel H ∞ synthesis [13]. The controller controlled the vertical speed of an aircraft while minimizing the influence of windshear and ground effects. Although the autonomous technologies can improve flight operations and overall aircraft performance, pilots will remain at the heart of operations in practice. Autonomous technologies are paramount to supporting pilots, thus enabling them to focus less on aircraft operation and more on strategic decision-making and mission management during landing. The energy approach to the flight control was used to access the current and predict the future states of an aircraft in order to improve the informational and situational awareness of the aircrew [14]. The landing maneuver was simulated with the ahead obstacle and the engine failure, which showed the correction of the algorithm. Another good project was Clean Sky 2 of European Clean Sky [15]. One of the themes of the project was to make use of advanced autonomous technologies and cockpit-navigation for large passenger aircraft to make the aircraft more reliable.
Although the above works make great contribution to aircraft autolanding system, there are still some problems that have not been solved. During the landing phase, the velocity of the aircraft is lower than that under the least drag and the aircraft is in boundary flight condition. Hence, the aircraft easily becomes unstable and even becomes a non-minimum phase (NMP) system. Most of the aforementioned papers do not include analyses on this situation. Trajectory tracking problems of NMP systems are challenging. Chen et al. [16,17] proposed the SI algorithm to solve the tracking problem of the NMP systems and got bounded solutions based on the inverse system and differential geometry theory. Subsequently, various studies for the application expansion of the SI algorithm have been implemented. Olivier et al. [18] used the SI algorithm for feedforward control of flexible multibody systems. Maghzaoui et al. [19] applied a poles placement state feedback control of an induction machine system while using the stable dynamic inversion methodology. A new method to calculate the causal solution of SI was proposed to precisely track the airspeed and altitude for unmanned aerial vehicles [20].
This paper proposes an integrated control method that combines the SI algorithm and H ∞ synthesis for civil aircraft autolanding system in order to solve the trajectory tracking and disturbance rejection problem simultaneously, with the adaptive ways of thinking to propose some algorithms of control that are of wider application [21]. The main contributions of this paper are summarized, as follows: (1). A robust autolanding controller (RAC) is designed based on H ∞ synthesis, which can handle with the disturbances during landing procedure, such as windshear and noise. (2). A stable inversion (SI) based robust autolanding controller (SIRAC) is proposed to improve the RAC scheme. The SI algorithm is used to enhance the trajectory tracking ability of aircraft autolanding system, which calculates the desired input and state through the desired landing trajectory. While the disturbance rejection ability is also increased due to the integration of the SI algorithm and H ∞ synthesis.
The rest of the paper is organized, as follows: Section 2 deals with the models of the aircraft, actuator, windshear, and the desired landing trajectory used in this study. Section 3 describes the design of the SIRAC. Section 4 presents some simulations to illustrate the SIRAC performances. Section 5 summarizes the conclusions.

Aircraft Dynamics and Actuator Modelling
The dynamic model of a civil aircraft can be built via Newton law, and in normal operating conditions, such a model can be decoupled into longitudinal motion and lateral motion, as they have a slight impact on each other. For the purpose of autolanding, the longitudinal motion model is employed here [22,23], which is where m is the aircraft mass, V is the longitudinal speed, α is the angle of attack, q is the pitch rate, θ is the pitch rate, γ is the flight path angle, ε T is the thrust inclination angle, I yy is the principal moment of inertia in pitch axis, M is the pitching moment, and T, L, and D are the thrust, lift and drag, respectively. The aerodynamics forces and moment in Equation (1) can be described, as follows: where the aerodynamic coefficients can be described, as follows: Table 1 provides the corresponding parameters of the civil aircraft Boeing 747 [24], which are used for simulations and analysis. Equation (1) can be trimmed and linearized according to small perturbation linearization theory. Table 2 lists the trim conditions. The elevator actuator and engine models are listed, as follows: δ t δ tc = 0.25 s + 0.25 (5) where δ e is the elevator deflection angle, δ t is the throttle position changing, and δ ec and δ tc are commands. The maximum deflections and rates of elevator and engine are limited in order to ensure the aircraft's safety and comfortability, which are listed in Table 3 [25].

Windshear Model
Windshear is a rapid variation in the velocity and direction of air flows. If the diameter of the windshear is less than 4 km, then it is called a microburst, otherwise a macroburst. Although windshear might last only a few minutes, it is one of the most dangerous factors in the takeoff and landing of an aircraft at low altitude because of its extreme speed and variation.
Researchers have developed many different models of windshear. Woodfield and Wood developed the vortex-ring model [26]. This paper uses the simplified vortex-ring downburst model [27] in simulations. As the lateral motion is ignored, the horizontal and vertical wind speeds are given, as follows: where W x is the horizontal wind speed, W h is the vertical wind speed, D = V 0 T/2, x = V 0 t, V 0 is the approaching speed, T is the total time during which the aircraft flies in the windshear, and f x and f h are the strengths of the windshear.
As the windshear speeds will influence the aircraft landing maneuver, Equation (7) needs to be embedded into the aircraft dynamics in the flight simulation. Additionally, the longitudinal motion equations can be written, as follows: Appl. Sci. 2020, 10, 1224 6 of 22

The Desired Landing Trajectory
In the alignment phase, the aircraft altitude will be around 500 m. Additionally, it keeps at a constant speed. In this case, the approaching speed is selected as 67.4 m/s. After that, the aircraft will enter the glide slope mode, when the approach path reaches the desired glide path. In this mode, the flight path angle should be kept at −3 deg and the gliding velocity should be maintained at a constant value. The aircraft enters the flare phase when the altitude is around 15 m. The aircraft will gradually approach the ground and the vertical speed will be reduced to 0.3 m/s at the touch down point.
From above analysis, two variables need to be controlled, such that the aircraft follows the desired trajectory: one variable is the longitudinal speed V and another is the altitude h. According to the civil aircraft landing requirements, the desired longitudinal speed V d is almost constant during the landing stage. As the vertical speed is very small, the horizontal speed V x is assumed to be equal to The desired altitude trajectory h d includes two parts: glide slope and flare segment. In the glide slope, the flight path angle is −3 deg, so the desired trajectory is In the flare segment, the desired trajectory is chosen as [28] where the k 1 , · · ·, k 4 meet the following equations: .
where t a is the total time of the flare segment, which is chosen as t a = 10s. We choose the initial k 0 = k 1 0 , k 2 0 , k 3 0 , k 4 0 = [0.0001, 0.01, 0, 5] and the allowable error is 10 −4 in order to find k 1 , · · ·, k 4 .
Using the function f solve(·) in Matlab to get the following results:

The SI Algorithm
While considering the NMP linear system: then the H ∞ controller can stabilize the output trajectory.

The SI Algorithm
While considering the NMP linear system: .
where x ∈ R n , u ∈ R p , y ∈ R m , z(t) is the measurement output, and A, B, C, D, C m , and D m are suitable dimension matrices. If y d (t) is the desired trajectory, the desired input u d (t) and state x d (t) meet the following equations: then the H ∞ controller can stabilize the output trajectory.
where ξ(t) includes the output y(t) and its time derivatives, that is, where r = [r 1 , r 2 , · · · , r m ] is the relative degree vector and C p is the pth row of C. For the aircraft, the state According to [29], the relative degree r = [3,2], so the ξ(t) can be selected as: Additionally, the system internal state η(t) = [q(t), θ(t)] T . Now, Equation (15) can be rewritten as: In order to find T, assume Substitute Equation (21) into Equation (17) yields According to Equations (20) and (22), we get Currently, the system Equation (13) can be rewritten into the new coordinates, as: . . where Without the loss of general, assume that the desired output trajectory y d (t) and its time derivatives are specified, that is, y (r) (t) = y (r) d (t). Subsequently, the following steps are used to find the desired input u d (t) and desired state x d (t).
Step 1: Find the inverse system.
This is the inverse system.
A transformation M is used to decouple the system Equation (31) into a stable subsystem σ s (t) and an unstable subsystem σ u (t), that is, If A s and A u are eigenvalues ofÂ η , then diag(A s , A u ) = MÂ η M −1 . Additionally, Equation (31) can be rewritten as: The stable subsystem is integrated from t = 0 to t → ∞ , while the unstable subsystem is integrated from t = 0 to t → −∞ . Subsequently, the bounded solution η d (t) can be found.
Step 3: Calculate u d (t) and x d (t).
If a bounded solution of η d (t) for the Equation (31) can be found, then the desired input u d (t) can be obtained through Equation (29), and the associated desired state x d (t) is obtained, as: Figure 3 shows the whole process of the SI algorithm.
If a bounded solution of ( ) d for the Equation (31) can be found, then the desired input ( ) d u t can be obtained through Equation (29), and the associated desired state ( ) d x t is obtained, as: Figure 3 shows the whole process of the SI algorithm.  Figure 3. Flowchart of the SI algorithm.

Design of RAC
H ∞ control is applied to process the wind disturbance. The system equations belonging to ( ) P s in Figure 4 are given, as follows: where ( ) u t is the control signal, ( ) w t is the exogenous input, ( ) z t is the exogenous output, and ( ) y t is the sensed output. Additionally, the RAC ( ) K s can be represented, as follows:

Design of RAC
H ∞ control is applied to process the wind disturbance. The system equations belonging to P(s) in Figure 4 are given, as follows: .
where u(t) is the control signal, w(t) is the exogenous input, z(t) is the exogenous output, and y(t) is the sensed output. Additionally, the RAC K(s) can be represented, as follows: and The desired performance criterion to be minimized can be taken as a H ∞ norm of the closed loop transfer function. Hence, the cost function takes the following form [30] ( ) ( , ) where ( ) K s should make the closed-loop system stable and satisfy the following H ∞ norm constraints ( ) , 0 Linear Matrix Inequalities (LMIs) are used to solve the H ∞ optimal control problem for a minimum γ without losing convexity in order to find the controller ( ) K s . Next, we give three steps to find the controller ( ) K s [31]. Step 1: Find symmetric matrices R and S .  Substitute u(t) in Equation (36) into Equation (35) and y(t) in Equation (35) into Equation (36), the closed-loop system can be written as: and The desired performance criterion to be minimized can be taken as a H ∞ norm of the closed loop transfer function. Hence, the cost function takes the following form [30] where K(s) should make the closed-loop system stable and satisfy the following H ∞ norm constraints Linear Matrix Inequalities (LMIs) are used to solve the H ∞ optimal control problem for a minimum γ without losing convexity in order to find the controller K(s). Next, we give three steps to find the controller K(s) [31].
Step 1: Find symmetric matrices R and S.
where N R is the basis of the null space of B 2 D 12 T , and N S is the basis of the null space of Step 2: Find symmetric positive definite matrix X.
When the R and S are found, we solve MN T = I − RS, and let where Step 3: Find the controller K(s).
When the X is found, according to the bounded real lemma, if S(P, K) is stable and J ∞ (K) < γ, then X should satisfy the LMI Let Inequality (48) left multiplies diag F T 1 , I, I and right multiplies diag(F 1 , I, I) and substitute A cl , B cl , C cl , and D cl into Inequality (48) yields where diag(·) represents diagonal matrix and Now, the controller K(s) can be constructed, as follows: The above method has been expanded into RAC design. Figure 5 shows the interconnected structure of the general plant for the RAC.    Where n is the sensors noises, W in , W p , W e , W act , and W n are the weighting functions. In this case, h and V are the variables to be controlled. The weighting functions are selected, as follows, according to the specifications to achieve the desired performances. There are two main tuning criteria that are used for selecting the weighting functions of the H ∞ controller: one is the weighting functions of some components of a vector signal should be bigger if they are more important than others, the other is each component of the signal should be scaled to the same units according to the weighting functions to make these components comparable.
The altitude h d and speed V d commands are scaled by W in in order to normalize the reference inputs. Select the average altitude h ave = 250 m and average speed V ave = 67.4 m/s, and yields: The gain of W e should be big enough at low frequency to enable the controlled system to track ramp commands with a very small steady state error, as the altitude error e h and the speed error e V should be small, and the altitude and speed do not change quickly. Hence, the W e is selected as: where 20,000 and 12,000 are proportional gains and 1/0.007 and 1/0.002 are time constants. It is obvious that the magnitude of W e is big at low frequency, which is similar with a proportional integral (PI) element, as the integral element can eliminate the steady state error.
As the maximum pitch rate is ±0.52 rad/s [25], the weighting function of W p is selected as: The weighting function of W act is applied to confine the control deflections and their rates. According to the maximum deflections and rates of elevator and engine in Table 3, the W act is selected as: As measurements are often corrupted with some noises, according to the measurement noise error [25], the measurement noise weighting function W n is selected, as: Based on the above weighting functions, RAC is constructed following the preceding three steps and the matrices are given, as follows:    Figure 6 shows the architecture of SIRAC, which includes SI algorithm and RAC. The SIRAC can simultaneously meet trajectory tracking and disturbance rejection requirements. The desired input u d and state x d can be calculated off-line based on the desired landing trajectory while using the SI algorithm. The RAC is used to robust against the wind disturbances. The SI drives the output trajectory of the aircraft to the desired trajectory. Now, the algorithm of the SIRAC is designed, as follows.
Step  Step 7 Knowing u , solve the system equation of aircraft for z .
Step 8 Iterate over Step 3 to 7 until the landing process is done.

Simulation Analysis
While considering the linearized aircraft model of Equation (13), matrices A and B have been already given in Section 2.1, and the other matrices are listed, as follows:  Now, the algorithm of the SIRAC is designed, as follows.
Step 1 Calculating Step 2 Knowing y d , solve Equation (29) and Equation (34) using SI algorithm to find the u d and x d .
Step 3 Knowing u d and x d , solve z d = C m x d + D m u d .
Step 6 Knowing [δ ec , δ tc ] T , solve u d + [δ ec , δ tc ] T for the input of the aircraft: u.
Step 7 Knowing u, solve the system equation of aircraft for z.
Step 8 Iterate over Step 3 to 7 until the landing process is done.

Simulation Analysis
While considering the linearized aircraft model of Equation (13), matrices A and B have been already given in Section 2.1, and the other matrices are listed, as follows: To find the u d and x d in step 2 of SIRAC algorithm, we need to know y d and its derivatives. According to Equations (10) and (11), in the glide slope, In the flare segment, The proposed method is also compared with the method in [4] for the windshear effect, which developed a linear quadratic regulator (LQR) feedforward closed-loop controller for aircraft landing. The brief introduction is as follows. More details can be found in [4].
Defining a new variable ε as: where t f is the final time. While combining Equation (59) with the linearized aircraft model, it can be shown that where d is the measurable disturbances. Using the assumption of constant steady-state values of y d , the Equation (60) can be written, as: where While combining the optimal LQR with the integrator and feedforward controller, the form of the final control law is where It is evident that there are three terms in Equation (62). The first is the LQR, the second term is the integral control, and the third term is the feedforward of the disturbances and the desired trajectory. Figure 7 shows the block diagram of the LQR feedforward closed-loop controller.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 23 It is evident that there are three terms in Equation (62). The first is the LQR, the second term is the integral control, and the third term is the feedforward of the disturbances and the desired trajectory. Figure 7 shows the block diagram of the LQR feedforward closed-loop controller.
Now, all of the models and methods that are employed in this paper have been established. Two different scenarios are simulated while using these models and methods to demonstrate the trajectory tracking and disturbance rejection capacity of SIRAC. In the first scenario, the landing Now, all of the models and methods that are employed in this paper have been established. Two different scenarios are simulated while using these models and methods to demonstrate the trajectory tracking and disturbance rejection capacity of SIRAC. In the first scenario, the landing procedure of the aircraft is simulated under no wind disturbance. The results are compared with the LQR and RAC to demonstrate the trajectory tracking capacity of SIRAC. In the second scenario, the windshear effect is taken into consideration in the landing process to demonstrate the disturbance rejection capacity of SIRAC.

Scenario 1
In this scenario, the landing procedure is simulated under no wind disturbance. Figure 8a compares the actual trajectory and desired trajectory. Additionally, Figure 8b shows the errors between actual trajectory and desired trajectory. It is clearly seen that, when compared with the LQR and RAC, the actual trajectory more closely follows the desired trajectory in glide slope and flare segment while using SIRAC. Additionally, the trajectory error of SIRAC is less about 3.5 m and 3 m than the LQR and RAC, respectively. The velocity of the SIRAC is also more close to the desired velocity, which is 67.4 m/s, see Figure 8c. Based on this result, it is also expected that the descent rate is slow. As seen from Figure 8d, the descent rates of three controllers are both approximately −3.5 m/s in the glide slope and reduce to zero in the flare segment. However, at the transition, the descent rate of SIARC has a larger variation, which will enable the aircraft to track the trajectory more quickly and closely. Figure 8e shows the response of pitch angle and Figure 8f shows the response of pitch rate. It can be seen that the pitch angle and pitch rate of three controllers are smoother in the glide slope. Similarly, the variation of SIRAC is larger in the transition, as it more closely tracks the trajectory.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 23 procedure of the aircraft is simulated under no wind disturbance. The results are compared with the LQR and RAC to demonstrate the trajectory tracking capacity of SIRAC. In the second scenario, the windshear effect is taken into consideration in the landing process to demonstrate the disturbance rejection capacity of SIRAC.

Scenario 1
In this scenario, the landing procedure is simulated under no wind disturbance. Figure 8a compares the actual trajectory and desired trajectory. Additionally, Figure 8b shows the errors between actual trajectory and desired trajectory. It is clearly seen that, when compared with the LQR and RAC, the actual trajectory more closely follows the desired trajectory in glide slope and flare segment while using SIRAC. Additionally, the trajectory error of SIRAC is less about 3.5 m and 3 m than the LQR and RAC, respectively. The velocity of the SIRAC is also more close to the desired velocity, which is 67.4 m/s, see Figure 8c. Based on this result, it is also expected that the descent rate is slow. As seen from Figure 8d, the descent rates of three controllers are both approximately −3.5 m/s in the glide slope and reduce to zero in the flare segment. However, at the transition, the descent rate of SIARC has a larger variation, which will enable the aircraft to track the trajectory more quickly and closely. Figure 8e shows the response of pitch angle and Figure 8f shows the response of pitch rate. It can be seen that the pitch angle and pitch rate of three controllers are smoother in the glide slope. Similarly, the variation of SIRAC is larger in the transition, as it more closely tracks the trajectory.

Scenario 2
In this scenario, the landing procedure is simulated in windshear to examine the robustness of the SIRAC. Windshear is modeled, as in Section 2.2, and the approaching speed 0 V is 67.4 m/s, the time = 60 s T during which the aircraft will fly in the windshear and the strengths of the windshear u f and w f are both set as 1.5. Figure 9 shows the horizontal speed and vertical speed of the windshear, respectively. The aircraft will meet an increasing horizontal headwind and the maximum speed is about −14 m/s. Subsequently, the horizontal headwind will become horizontal downwind and the maximum speed is about 14 m/s. In the meantime, there is a downward flow, and the maximum speed is about −15 m/s.  Figure 10a shows a comparison of the altitude response between the actual trajectory and desired trajectory in a windshear environment. Additionally, Figure 10b shows the altitude error of LQR, RAC, and SIRAC. It can be seen that the largest altitude error of SIRAC is less about 10 m and 15 m than RAC and LQR in the glide slope and both 3 m in the flare segment. The velocity response of SIRAC is also smoother and more closed to desired velocity than LQR and RAC, as shown in Figure 10c. The largest velocity amplitude of SIRAC is less about 5 m/s and 6m/s than RAC and LQR, respectively. In the windshear, it is also expected that the descent rate is slow and smooth. The descent rate of SIRAC is more closed to −3.5 m/s, as the largest error of SIRAC is about 1 m/s, RAC

Scenario 2
In this scenario, the landing procedure is simulated in windshear to examine the robustness of the SIRAC. Windshear is modeled, as in Section 2.2, and the approaching speed V 0 is 67.4 m/s, the time T= 60 s during which the aircraft will fly in the windshear and the strengths of the windshear f u and f w are both set as 1.5. Figure 9 shows the horizontal speed and vertical speed of the windshear, respectively. The aircraft will meet an increasing horizontal headwind and the maximum speed is about −14 m/s. Subsequently, the horizontal headwind will become horizontal downwind and the maximum speed is about 14 m/s. In the meantime, there is a downward flow, and the maximum speed is about −15 m/s. (e) (f)

Scenario 2
In this scenario, the landing procedure is simulated in windshear to examine the robustness of the SIRAC. Windshear is modeled, as in Section 2.2, and the approaching speed 0 V is 67.4 m/s, the time = 60 s T during which the aircraft will fly in the windshear and the strengths of the windshear u f and w f are both set as 1.5. Figure 9 shows the horizontal speed and vertical speed of the windshear, respectively. The aircraft will meet an increasing horizontal headwind and the maximum speed is about −14 m/s. Subsequently, the horizontal headwind will become horizontal downwind and the maximum speed is about 14 m/s. In the meantime, there is a downward flow, and the maximum speed is about −15 m/s.  Figure 10a shows a comparison of the altitude response between the actual trajectory and desired trajectory in a windshear environment. Additionally, Figure 10b shows the altitude error of LQR, RAC, and SIRAC. It can be seen that the largest altitude error of SIRAC is less about 10 m and 15 m than RAC and LQR in the glide slope and both 3 m in the flare segment. The velocity response of SIRAC is also smoother and more closed to desired velocity than LQR and RAC, as shown in Figure 10c. The largest velocity amplitude of SIRAC is less about 5 m/s and 6m/s than RAC and LQR, respectively. In the windshear, it is also expected that the descent rate is slow and smooth. The descent rate of SIRAC is more closed to −3.5 m/s, as the largest error of SIRAC is about 1 m/s, RAC  Figure 10a shows a comparison of the altitude response between the actual trajectory and desired trajectory in a windshear environment. Additionally, Figure 10b shows the altitude error of LQR, RAC, and SIRAC. It can be seen that the largest altitude error of SIRAC is less about 10 m and 15 m than RAC and LQR in the glide slope and both 3 m in the flare segment. The velocity response of SIRAC is also smoother and more closed to desired velocity than LQR and RAC, as shown in Figure 10c. The largest velocity amplitude of SIRAC is less about 5 m/s and 6 m/s than RAC and LQR, respectively. In the windshear, it is also expected that the descent rate is slow and smooth. The descent rate of SIRAC is more closed to −3.5 m/s, as the largest error of SIRAC is about 1 m/s, RAC is 2.5 m/s, and LQR is 3 m/s, as seen from Figure 10d. The pitch angle and pitch rate responses are also shown in Figure 10e,f, respectively. It is clearly seen that when the aircraft encounters a windshear, the pitch angle response has a biggest variation 8deg while using SIRAC, 15 deg using RAC, and 18 deg using LQR. Additionally, the pitch rate also has a larger variation in the transition from initial approach to glide slope and glide slope to flare segment while using SIRAC, which means that the transition time of SIRAC is shorter. As a result, it is evident that the SIRAC provides a more precise and steady landing.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 20 of 23 is 2.5 m/s, and LQR is 3m/s, as seen from Figure 10d. The pitch angle and pitch rate responses are also shown in Figure 10e and Figure 10f, respectively. It is clearly seen that when the aircraft encounters a windshear, the pitch angle response has a biggest variation 8deg while using SIRAC, 15 deg using RAC, and 18 deg using LQR. Additionally, the pitch rate also has a larger variation in the transition from initial approach to glide slope and glide slope to flare segment while using SIRAC, which means that the transition time of SIRAC is shorter. As a result, it is evident that the SIRAC provides a more precise and steady landing.

Conclusions
This paper presented a new approach to the ALS of the civil aircraft. Differently from the existing control methods for the autolanding system, this approach exploited the integration of the SI algorithm and H ∞ synthesis, which could precisely track the autolanding trajectory, even in large wind disturbance. The windshear disturbance was predefined. Subsequently, the desired input and state were calculated while using SI algorithm according to the desired altitude trajectory and the desired velocity. The general plant for RAC was created. Weighting functions were selected according to the desired performances. RAC was constructed by solving the LMIs. Finally, two scenarios were simulated and studied. The simulation results demonstrated the trajectory tracking and disturbance rejection capacity of the SIRAC.
In this paper, it is assumed that the system parameters are accurate and do not change over time. However, a real system generally has system parameter uncertainties. Thus, future work will be concentrated on the robustness of the system by considering system uncertainties and turbulence. Although the paper is based on mathematical models and just shows the simulation results computed by computer, the production of prototypes has been basically completed, which indicates that our laboratories and cooperative industrial companies have reached the Technology Readiness Levels 5 (TRL5) [32].  Acknowledgments: Authors would like to thank the editors and the reviewers for their constructive suggestions and comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature m
aircraft mass δ e , δ t elevator deflection angle and throttle position changing, respectively δ ec , δ tc elevator command and throttle command, respectively h, h d altitude and desired altitude, respectively V, V d longitudinal speed and desired longitudinal speed, respectively α, θ, γ angle of attack, pitch angle and flight path angle, respectively q pitch rate ε T thrust inclination angle I yy principal moment of inertia in pitch axis M pitching moment T, L, D thrust, lift and drag, respectively C L , C D , C M aerodynamic coefficients of lift, drag and pitching moment, respectively