Receding Horizon Longitudinal Control Technology for Automatic Carrier Landing With Variable Reference Trajectory Based on Sliding Rate Information

In this paper, the longitudinal control of automatic carrier landing is studied. First, the carrier landing control problem is transformed into an optimal control problem of trajectory tracking. Considering the constraints of the control variables and the rate of change of control variables in the realistic landing process, the original linear small disturbance model is expanded. Based on the symplectic pseudospectral method and the adaptive regression prediction technology, a fast receding horizon carrier landing control technology with a variable reference trajectory is developed. Finally, the effectiveness of the control algorithm is verified by simulations at different sea states, initial deviations, and reference trajectory selection strategies. The simulation results demonstrate that the introduction of deck motion prediction can greatly reduce the phase delay of the control system and enhance the tracking ability of the carrier-based aircraft and improve the control effectiveness significantly. The proposed algorithm can precisely control the carrier landing trajectory under initial deviations, the external continuous wind disturbances, and random error of the state variables. Additionally, the calculation efficiency of the present control algorithm is sufficient for real-time online tracking.


I. INTRODUCTION
As the main combat weapon system of aircraft carrier formations, carrier-based aircraft greatly enhance the mastery of the seas, air superiority, and integrated attack-defense capabilities because of the high mobility. The problem of carrier landing is crucial to the effective cooperation between carrier-based aircrafts and the carrier. Carrier landing are significantly different from ground landing. Specifically, the landing area of the carrier's flight deck is limited, and the deck is constantly in motion, making the aircraft landing much more difficult. Changing meteorological conditions can cause the above adverse effects to be intensified. In order to assist the pilot in carrying out a secure carrier landing The associate editor coordinating the review of this manuscript and approving it for publication was Shihong Ding . in all weather conditions, the automatic carrier landing system (ACLS) of carrier-based manned aircraft has been greatly developed.
To ensure that the carrier-based aircraft can fly according to the ideal glide path in the ACLS system, a corresponding control algorithm must be designed to maintain the trajectory of the carrier-based aircraft. As a classical control method, Proportional Integration Differentiation (PID) control has been applied to the carrier landing control problem for many years [1], [2]. To improve the effectiveness of PID control, Zhiyuan Yang, Yimin Deng, and others have proposed a variety of parametric tuning methods utilizing intelligent optimization algorithms [3], [4]. By using these approaches, satisfactory control effectiveness has been obtained. Traditional PID control methods can achieve satisfactory landing performance in normal situations. How-ever, it is difficult to track the randomly changing reference glide path under the disturbances of deck motion [5]. Hengzhi Ding used the theory of nonlinear inverse dynamics combined with PID control to design an automatic carrier landing system, which can efficiently resist carrier air wake [6]. The effectiveness of the dynamic inverse method is heavily dependent on the accuracy of the model. Therefore, the effectiveness degrades significantly when there are random wake and sea wave disturbances. To achieve better anti-interference capabilities of the ACLS, Wang et al. [7] developed a stable adaptive control scheme based on LDU (lower-diagonal-upper) decomposition of the high-frequency gain matrix. This approach ensures closed-loop stability and asymptotic output tracking. Zhen et al. [8] studied a multivariable model reference adaptive control (MRAC) scheme for the automatic carrier-landing of unmanned aerial vehicles (UAVs) subject to nonlinear system dynamics, multivariable coupling, and parametric uncertainty. Lungu et al. [9] and Ju and Tsai [10] applied backstepping control to ground landing. The approach [10] addressed model parameter uncertainty in the path tracking process and it provided a reference for the glide path tracking problem of carrier landing. Wu et al. [11] studied robust carrier landing control by decoupling the height and velocity channels using an exact linearization method. Subsequently, robust controllers for each of the two channels were designed separately. In addition, many scholars have applied optimal control [12], [13], preview control [14], [15], sliding control [16], fuzzy control [17] and adaptive control [18], [19] technologies to the carrier landing trajectory control, and satisfactory control results are obtained.
The control algorithms discussed above play an important role in the research of automatic carrier landing control, and the control effect of PID algorithm has been verified in practice. However, the existing algorithms usually cannot deal with the control constraints explicitly, and when there are multiple control variables, there is no effective control allocation method. In addition, the above control algorithms generally divide the control system into multiple control loops to design separately, which is difficult to ensure the optimization of the whole control system. Carrier landing control is essentially a trajectory tracking problem, so some related control techniques proposed in recent years, such as quantized static output feedback control [20], [21], non-smooth control [22], homogeneous domination approach [23] and model predictive control (MPC) [24] can also be used for reference. Among the above methods, the MPC method has the unique advantage in its ability to treat constrained problem and utilize the forecast information. These abilities are very expected for automatic carrier landing. However, the basic idea of the MPC method is to solve a finite horizon open-loop optimal control problem at each sampling instant. Hence, the solving efficiency is the bottleneck that limits its wide application in automatic carrier landing. Various numerical methods for solving nonlinear optimal control have been developed, which can be roughly divided into direct methods and indirect methods [25]. In the last two decades, pseudospectral methods have become the most popular in the field of aeronautics and astronautics due to its high precision and sound robustness [26]. For years, pseudospectral methods were merely studied under the framework of direct methods. However, as a kind of discretization schemes, its application should not be limited only to direct methods. Recently, by utilizing the pseudospectral scheme, Peng and Wang et al. creatively develop a series of symplectic pseudospectral methods under the framework of indirect methods [27], [28]. These symplectic pseudospectral methods have excellent efficiency and accuracy due to the structure-preserving property [29]. And the successive convexification technique is integrated to further achieve excellent numerical robustness and fast convergence [30]. Various numerical tests demonstrate that the symplectic pseudospectral method is an appealing numerical method for solving trajectory planning problems [31], [32]. Owing to the good numerical characteristic of symplectic pseudospectral methods, together with the idea of receding horizon control, a fast receding horizon carrier landing control technology with the variable reference trajectory based on the glide rate information is designed in this paper. The contribution of this work is threefold.
First, the phase delay of the control system is greatly reduced by introducing prediction of the deck motion and the expected glide rate information into the selection progress of the reference trajectory. As a result, the following capability of the carrier-based aircraft with respect to deck motion can be significantly enhanced. Second, the derivative constraint problem of control variables in the actual control system is handled by augmenting the linear model to produce a control law that meets the requirements of practical engineering tasks. Third, by introducing Runge-Kutta integration into the original optimal tracking control method based on Receding Horizon Control (RHC), the influence of aerodynamic wake disturbances during the landing process are easily dealt with, and the simulation results are ensured to be realistic.
The remainder of this paper is structured as follows. In Section II, the problem definition and models are formulated. These models include the longitudinal linear model of small disturbances, the deck motion model, and the carrier air wake model. Section III describes the receding horizon optimal control algorithm for carrier landing based on symplectic pseudospectral algorithm. Section IV presents simulation results testing the control algorithm. Finally, Section V briefly summarizes the conclusions and discusses the further research orientation.

II. LONGITUDINAL CONTROL MODEL OF CARRIER LANDING A. LONGITUDINAL SMALL DISTURBANCE EQUATION OF CARRIER LANDING
Compared with lateral carrier landing control, longitudinal control is more difficult and the control performance VOLUME 8, 2020 have a greater impact on the final landing appearance. So only the longitudinal control problem is discussed in this paper.
In this paper, the linear small disturbance equations of the F/A-18A in the longitudinal direction are used to describe the landing process of the carrier-based aircraft [3], [33]: where v, α, θ, q, h represent the velocity, angle of attack, pitch angle, pitch angular velocity, and altitude deviation of the aircraft relative to the nominal state, respectively. δ H , δ LEF , δ RT , δ PL represent the deviation of horizontal tail deflection, leading-edge flap deflection, rudder toe-in deflection, and engine throttle control angle, respectively. α g represents the deviation in attack angle caused by vertical wind disturbances. With the exception of δ PL , the above variables are all in international standard units. The unit of δ PL is deg. The trim state of carrier-based aircraft is: V 0 = 69.96m/s, α 0 = 8.3 • , γ 0 = −3 • . In actual flight control systems, there are some limits on the range and rate of change of the rudder deflection and throttle control angle. These limits result in the following limitations on the control variables in the carrier landing control model: In Eq. (2), the superscript * donates the nominal value of rudder deflection or throttle control angle, and * max , * min represent the maximum and minimum permissible values of the variable, respectively.

B. DECK MOTION MODEL
When an aircraft carrier is at sea, it is subject to pitching, rolling, and heave motions due to both wind and waves. The pitching and heave motions of the carrier have the greatest impact on the longitudinal landing accuracy of the carrier-based aircraft. The pitching and heave motion model of the carrier deck for a carrier sailing at a typical speed of 15.4m/s can be found in Reference [34]: where ϕ 1 and ϕ 2 are the random initial phases.

C. CARRIER AIR WAKE MODEL
In the present work, the engineering model of carrier air wake in Reference [34] is adopted. Since the longitudinal landing dynamics model is only affected by the air wake of vertical direction in Eq. (1), only the simulation result of air wake in vertical direction is presented in Fig. 1.

III. RECEDING HORIZON OPTIMAL CONTROL ALGORITHM FOR CARRIER LANDING A. TRAJECTORY TRACKING OPTIMAL CONTROL PROBLEM
The optimal control problem with constraints can be written as follows: where S(x, u, t) is the objective function of the optimal control system. f (x, u, t) represents the dynamic model of the control system. The inequality h ≤ 0 refers to constraints of the actual control system, including state variable constraints and control variable constraints. In the present paper, a symplectic pseudospectral algorithm based on the second kind of generating function [28] is used to solve the optimal control problem in the receding time window. As an indirect numerical methods for optimal control problems [35], this approach is known to have high calculation efficiency and fast convergence speed and can deal with the optimal control problem with constraints on state and control variables. Without losing generality, Eq. (4) can be rewritten as: where, By introducing a non-negative relaxation vector α, the inequality constraint in Eq. (5) can be rewritten as an equality constraint as follows: Furthermore, by introducing the costate vector λ and Lagrangian multiplier µ, the optimal control problem with constraints can be transformed into an unconstrained optimal control problem. As a result, the objective function can be expressed as: where H is the Hamilton function given by: According to the parametric variational principle, when the objective function J achieves a minimum value, the following equations should be satisfied simultaneously: According to the inequality constraints and the KKT (Karush Kuhn Tucker) condition: Eq. (5) is an optimal control problem with a Mayer type objective function and a fixed terminal time. Therefore, when the terminal state is fixed, the boundary condition x(t 0 + T ) = x f should be added. When the terminal state is free, In the jth sub-interval, the variables x, λ, µ and α are discretized by the N (j) th Legendre-Gauss-Lobatto (LGL) nodes: where ρ (j) i (τ ) is the Lagrangian interpolation polynomial corresponding to the LGL nodes in the jth sub-interval. The specific expression for this polynomial can be found in Reference [28].
Applying the stagnation point condition of the second generating function [27], [36] in the jth sub-interval, then According to the constraint in Eq. (7) and the complementarity condition, the following relationship can be obtained in the subinterval j By assembling Eq. (13∼14) in every interval according to the boundary conditions, the two-point boundary value problem in [t 0 , t 0 + T ] can be rewritten as follows, where,x,λ andμ contain the information of state vectors, costate vectors, and Lagrange multipliers at each of the LGL collocation points, respectively. According to Eq. (15a), the variablesx andλ can be expressed as linear functions ofμ separately. Therefore, the two-point boundary value problem given in Eq. (15) can be transformed into a standard linear complementarity problem based on Eq. (15b) and Eq. (15c), as shown in Eq. (16). The variableμ can be solved using the Lemke method [37], andx andλ can be solved according to the linear functions. Finally, the required VOLUME 8, 2020 control variable u is obtained by substitutingx andλ into Eq. (10) A detailed derivation process of the symplectic pseudospectral algorithm based on the second kind of generating function, as well as the meaning and specific expressions of the variables involved in Eq. (13)-(16), can be found in the Reference [28], [36].

B. LONGITUDINAL OPTIMAL CONTROL MODEL OF CARRIER LANDING
An aircraft must attempt to follow the ideal glide path during carrier landing process. The terminal point of the ideal glide path is located on the ideal landing point of the carrier deck. The inclination angle of the glide path is expressed as γ , as shown in Fig. 2. In Fig. 2, the coordinate systems O D and O P represent the deck coordinate system and the pitch coordinate system, respectively. The coordinate origin of O D is the ideal landing point on the deck, while the coordinate origin of O P is the pitch center of the aircraft carrier. The X -axes are both parallel to the straight section of the carrier deck and are directed forward. The Z -axes are both perpendicular to the deck and are directed vertically downward with respect to the deck. A virtual inertial coordinate system O V is defined to represent the average motion of the aircraft carrier. O V is primarily used to describe the base motion of the aircraft carrier. The average position of the carrier's pitch center during the navigation progress is taken as the coordinate origin, the X -axis is oriented to the average speed direction of the carrier's pitch center, and the Z -axis is perpendicular to the sea level and points downward. This coordinate system coincides with O P when the carrier is stationary.
During landing, it is necessary to track the ideal glide path under the constraints of flight dynamics, control variables, and state variables. The tracking error must be as small as possible, and the control process must be sufficiently stable during the trajectory tracking process. Based on the above analysis, the mathematical model of the optimal landing control problem of carrier-based aircraft can be described as Eq. (17), where X e represents the ideal flight path of the carrierbased aircraft and X represents the actual flight path of the carrier-based aircraft, which can be regarded as the actual state variable of the system here. U represents the actual control input.P andR are the diagonal coefficient matrices. P is required to be semi-positive definite whileR must be positive definite. The symbol≤ indicates that each component of the two vectors must satisfy an ≤ inequality. It can be seen from Eq. (5) According to Eq. (6), (17) and (2), it can be shown that: where the subscript (a: b) represents the a-th to the b-th components of the vector. It should be pointed out that the expression of W obtained from Eq. (6) is W =Ēα g . In Eq. (19), the physical quantity α g is caused by random wind disturbances. Additionally, the information of wind disturbance in the time interval [t 0 , t 0 + T ] cannot be obtained or predicted at t 0 . Therefore, the random disturbance cannot be handled when solving the optimal control model. To address this, this component of interference can be ignored in the model solving process, i.e. W = 0. The error caused by this simplification will be discussed in Section III-C. Inspired by Reference [38], a receding horizon carrier landing control algorithm with variable reference trajectory is designed in the present paper. The core idea of the receding horizon control method is to solve the optimal control problem over a specified future time interval and treat the resulting optimal control law at the current time as the actual input of the control system until the next calculation time point to repeat the process [39]. A time-series diagram of the VTGR-RHC algorithm is given in Fig. 3.
The symbols used in Fig. 3 are explained as follows: x e (t) : The desired state of the control system at time t.
x * (t) : The actual state of the control system at time t. Operation (n): Over the time interval n = [t 0 + nδ, t 0 + nδ + T ], according to the two-point boundary condition x (t 0 + nδ) = x * (t 0 + nδ), x (t 0 + nδ + T ) = x e (t 0 + nδ + T ) and the expected glide path x n e (t), Eq. (17) is solved using the symplectic pseudospectral algorithm introduced in Section III-A, where the matrices are selected according to Eq. (19)- (20). u opt (t 0 + nδ) : The optimal control law obtained by the calculation of Operation (n) over the time period [t 0 + nδ, t 0 + nδ + T ].
Taking u opt (t 0 + nδ), x * (t 0 + nδ) as the actual control input of the system and the initial condition, respectively, x * (t 0 + (n + 1) δ) is obtained by integrating the real control system of Eq. (17) in [t 0 + nδ, t 0 + (n + 1)δ]. The numerical integration method used in the present work is the Fourth-order-Runge-Kutta algorithm. By introducing the aerodynamic wake influence into the simulation process, the simulation error caused by the simplification of the matrix W in Eq. (19) can be eliminated.

2) SELECTION OF THE REFERENCE GLIDE PATH
According to the carrier landing optimal control model given in Eq. (17), there is a standard reference trajectory to be tracked in the process of carrier landing. However, when the aircraft carrier is at sea, it is accompanied by six degrees of freedom motion due to the wind and waves. As a result, the ideal landing point will move in response to the carrier deck motion. If the deck motion is ignored, the initial ideal landing trajectory will always be tracked as the straight trajectory. This can lead to significant landing errors. In order to solve this problem, two strategies can be adopted: a. The real-time motion information of the aircraft carrier can be monitored, and the straight-line trajectory between the current position of the carrier-based aircraft and the ideal landing point is regarded as the reference trajectory in the current time window.
b. The motion of the carrier can be predicted, and the future position of the ideal landing point is predicted. Subsequently, the straight-line trajectory between the current position of the carrier-based aircraft and the predicted position of the ideal landing point is used as the reference trajectory in the current time window.
The two strategies above correspond to two different ideal glide paths, as shown in Fig. 4. The hexagon stars in the figure represent the different positions of the ideal landing point according to the different strategies. Compared to Strategy (b), Strategy (a) is relatively simple and does not require a prediction of the deck motion. However, it can be seen from the later simulation results that the reference glide path generated by Strategy (a) results in a significant phase delay and poor control. Therefore, the present work is focused on Strategy (b).

3) PREDICTION OF AIRCRAFT CARRIER DECK MOTION
According to Strategy (b), an ideal glide slope based on predicted deck motion needs to be generated. In order to improve the efficiency of the algorithm, an autoregressive (AR) prediction algorithm is selected to predict the deck motion. The AR prediction model of p-th order can be simplified as [40]: ξ t = a 1 ξ t−1 + a 2 ξ t−2 + · · · + a p ξ t−p + ς t (21) where ξ is the variable to be predicted, ξ t−j , j = 0, 1, 2 · · · p represent the historical data of ξ , a j , j = 1, 2 · · · p are undecided variables, and ς t donates random white noise. Define ξ N = ξ p+1 , ξ p+2 , · · · , ξ N T and a N = a 1 , a 2 , · · · , a p T , and where N donates the number of historical data used for forecasting. In order to minimize the objective function

the least-squares
method is used to estimate the undetermined coefficients, and the optimal estimated value of a N is obtained aŝ To ensure that Eq. (23) is meaningful, when A N is not column full rank, the singular value decomposition A N = U V T can be applied and A T N A N −1 A T N in Eq. (23) can be replaced by V −1 U T . Then the data at the next l steps can be predicted by the following equation: When the AR algorithm is applied to predict the deck motion, the prediction step length t p , the maximum prediction step l m , the model order p and the number of historical data N should be determined first. Generally, t p and l m can be determined according to the problem. Numerical simulation experiments based on the actual problem can be carried out to determine p and N . The order of the model p can be also determined by the AIC criterion or SBC criterion [41].

4) DESIGN OF THE REFERENCE TRAJECTORY BASED ON GLIDE RATE INFORMATION
After the predicted aircraft carrier motion is obtained, the corresponding glide rate of the reference trajectory in time interval [t 0 + nδ, t 0 + nδ + T ] can be calculated. It can be seen from Fig.4 The position vector of the ideal landing point relative to its original position in O V can be expressed as where, P v D is the position vector of the ideal landing point in O V , Further, the expected glide rate of the carrier-based aircraft can be calculated as where t r is the expected remaining time for carrier-based aircraft to reach the touchdown point. h represents the height error between the aircraft and the virtual ideal landing point, is the regulatory factor. The reason for the introduction of κ is that when the carrier-based aircraft approaches the ideal landing point, the value of K gp becomes rather large due to the small value of t r . This can exceed the control ability of the actuator and result in a divergence of the calculation. The introduction of κ also increases the following capability of the carrier-based aircraft at the beginning of the landing phase. For t r can be calculated as: K gp can be rewrote as: The expected glide path is generated according to the expected glide rate information: where, t = t − t 0 − nδ. Letting t = T in Eq. (29), the terminal boundary condition can be obtained as h e (T ). The receding calculation ends when the carrier-based aircraft reaches the deck,

5) THE FRAMEWORK OF THE VTGR-RHC ALGORITHM
As described above, the VTGR-RHC algorithm is graphically summarized in Fig. 5. In order to obtain accurate longitudinal landing error, the dichotomy method is applied after RHC calculation. The constraints of the carrier-based aircraft control system are shown in Table 1 [12], [42].

B. THE CARRIER LANDING CONTROL EFFECTIVENESS UNDER TYPICAL SEA STATE
Under typical sea conditions, the deck motion of an aircraft carrier can be described by the model given in Section II-B. Based on the deck motion model and the simulation conditions in Section IV-A, the process of the carrier landing control is simulated. The receding time window length is taken as 1.5s, and the receding step length is taken as 0.05s. In addition, the aircraft carrier is assumed to travel along a straight line at a fixed speed V S = 15.4m/s. In order to simulate a realistic situation, a random error within ± 1% is introduced into the state variables to replicate uncertain interference in the simulation process. The parameter ζ is taken as 0.012 based on experimental test. In order to avoid the control signal given by the algorithm exceeding the control ability of the actual actuator, the absolute value of K gp in the simulation is limited to 1.5.
The control effectiveness is demonstrated in Fig. 6, where Fig. 6(a) shows the glide path of the carrier-based aircraft and Fig. 6(b) shows the variation in the height of the carrier-based aircraft during the landing process. It can be seen from Fig. 6 that by applying the VTGR-RHC algorithm the aircraft can overcome the initial error interference within 3 seconds, and then the aircraft can accurately track the deck motion and follow the ideal landing trajectory. The error between practical flight path and expected flight path is the same with the error between variation height of the plane and variation height of the ideal touch point numerically. And this error result can be found in Fig. 9 as the 'typical sea state', from which we can see that the altitude position error of the aircraft can be limited in 0.25 m during the steady-state control process. And the final longitudinal error of the landing point on the deck is 1.1147m.
The simulations were performed in MATLAB R2016b on an Intel Core i5-6200U personal computer with a 2.3 GHz processor and 8 GB of RAM. The single-step solution time of the algorithm was 18.7ms which was much shorter than the step length. It means the algorithm has high computational efficiency and can meet the requirements of real-time online tracking.
The rate of change of each control variable during landing is shown in Fig. 7. Fig. 7(a-d) show the rate of change of elevator deflection, leading-edge flap deflection, rudder toein deflection, and engine throttle control angle, respectively. It can be seen from Fig. 7 that the rate of change of elevator deflection, leading-edge flap deflection, and throttle control angle reach the performance limits of the actuators in the control process (as shown in the dotted line). This demonstrates that the algorithm can effectively limit the control variables.
The result of state variable θ is shown in Fig. 8. From  Fig. 8, we can see that during the carrier landing process, the pitch angle of the aircraft only change in a small scale (about 0.7 deg in the last 10 seconds). That means the aircraft can hold a satisfied attitude for hooking cable during the carrier landing process.

C. SIMULATION RESULTS UNDER DIFFERENT SEA STATES
In this section, the effectiveness of the VTGR-RHC algorithm for carrier landing is analyzed under different sea states. The control effectiveness under typical (normal) sea state was shown in Fig. 6. The landing error results under low, medium-high and high sea states (corresponding to 0.7, 1.3, VOLUME 8, 2020 FIGURE 6. The control effect of VTGR-RCH algorithm. and 1.6 times of deck motion and air wake intensity) are shown in Fig. 9. From Fig. 9, it can be seen that the tracking error between the carrier-based aircraft and the ideal glide path is strongly dependent on the severity of sea states.
In order to make a statistical comparison of the landing control under different sea states, 50 different initial phases of waves were selected for simulation. The longitudinal deviation of the landing position and calculation time were then calculated for each initial state. The 50 individual cases with different initial phases were simulated under each sea state, and the statistical average results are shown in Table 2. Here, to improve the performance of the control method under the complex sea states and initial phases, the parameter ζ is taken as 0.013. From the data given in Table 2, it can be seen that the landing error increases with increased severity of the sea states. Table 2 provides the average statistics of the landing error results. In order to further analyze the distribution of landing errors, Fig. 10 shows the distribution of the 50 landing points on the deck under each sea state. It can be seen from Fig. 10 that with the increase of the sea state severity, the errors in landing accuracy as well as the error dispersion range increase gradually. Under low and normal sea states, all landing points are located between the second and the third arresting cables, which ensures the ideal situation of the third cable hooking. In most cases of the medium-high sea states, the carrier-based aircraft can achieve the third cable hooking, and the second cable hooking may occur on rare occasions. While under the   high sea state, the fourth cable hooking may occur in some cases, which is barely acceptable in typical scenarios. These results indicate that a further increase in the severity of the sea state may lead to a landing failure. For this reason, the navies of all countries have strict rules on the sea state of carrier landing.

D. THE CARRIER LANDING CONTROL EFFECTIVENESS UNDER DIFFERENT INITIAL DEVIATION CONDITIONS
Similarly, considering the control performance of the algorithm under different initial deviation conditions, three groups of different initial errors (1.5, 0.5 and -1 times of x 0 , respectively) were selected for simulation under a typical sea state.  The resulting control effectiveness are shown in Fig. 11(a). Fig. 11(b) shows a comparison of landing error under four different initial deviation conditions. From Fig. 11, it can be seen that the altitude control errors in the last 10 seconds under different initial error conditions are all less than 0.2m and the different initial error conditions has little influence on the final landing control effect.
A total of 50 carrier landing simulations were carried out under each initial error condition in order to produce statistics of longitudinal landing deviation and calculation time. The statistical average results are shown in Table 3. From these results, it can be seen that the average longitudinal landing  deviations under different initial conditions are all less than 1.3m and there is no obvious correlation between the relative size of the initial error and the final landing error. What's more, the average computation times of the whole landing process are around 5.5s, which is only about 1/3 of the real time range.
In order to further explore the influence of different sea states and different initial deviations on the final landing control effect, the two cases described above were simulated again with the error caused by random interference being ignored. The results are shown in Fig. 12. Fig.12 further confirms that the tracking error between the aircraft and the ideal glide path is positively correlated to the severity of sea states. Conversely, the relative size of the initial error has no significant impact on the final landing control error.

E. CONTROL EFFECTS OF DIFFERENT REFERENCE TRAJECTORY SELECTION STRATEGIES AND CONTROL METHODS
As mentioned before, there are two strategies for the selection of the reference trajectory used in the control algorithm (described in Section III-C). Fig. 13(a) shows the control effects of Strategy (a) and Fig. 13(b) shows the control effect comparison of these two selection strategies. Compare to Fig. 6, it can be seen that the phase delay of the control system is greatly reduced by introducing prediction of the carrier deck motion in the control process. As a result, the following capability of the carrier-based aircraft to the deck motion can be enhanced. As well, the height error of the aircraft can be reduced from ± 1.5m to ± 0.25m.
In the online optimal tracking control method based on RHC proposed in the literature [38], interpolation is used to obtain the state variables at the current time. This operation makes it impossible to deal with the influence of the aerodynamic wake disturbances on the carrier-based aircraft, and therefore, the simulation results do not accurately reflect the real physical scenario. The Runge-Kutta integration technique applied in this paper can easily deal with the influence of the aerodynamic wake disturbances.
Accounting for aerodynamic wake disturbances, the two methods above were used to simulate the carrier landing control process, respectively. The results are shown in Fig. 14(a), and the height deviation of the two methods is shown in Fig. 14(b). From Fig. 14, it can be seen that the interpolation operation will lead to extra height error by a maximum of 0.25m during the simulation due to the absence of aerodynamic wake disturbance. It should be noted that the random disturbance factors were not considered in Fig.13 and Fig. 14 in order to obtain a clearer conclusion.
Finally, to confirm the algorithm's effectiveness, comparison experiments are done using typical LQ (Linear Quadric) feedback optimal control algorithm and VTGR-RCH algorithm. The simulation result of LQ optimal control is shown in Fig. 15(a), and the height deviations of the two methods are shown in Fig. 15(b). From Fig. 15, it can be seen that the typical LQ algorithm cannot effectively track the movement of the carrier, and the steady-state error is far greater than VTGR-RCH algorithm. Within the first 2.5 seconds of the control process, the height error of typical LQ algorithm is less than that of VTGR-RCH. This is due to that the control variables cannot be constrained effectively by typical LQ algorithm and the initial height error can be restrained rapidly under the over-limit control variables. While the over-limit control variables cannot be obtained in the real physical scenario. This also illustrates that the VTGR-RCH algorithm is more practical in the actual situation of engineering. It should be noted that the random disturbance factors were not considered in this simulation, either.

V. CONCLUSION
In this paper, the problem of longitudinal automatic carrier landing control is studied. First, the carrier landing control problem of aircraft was transformed into an optimal control problem of trajectory tracking. Further, the original linear small disturbance model was extended to account for constraints on control variables and derivatives of control variables in the landing process. With the help of a symplectic pseudospectral algorithm based on the second kind of generating function and adaptive regression prediction technology, a receding horizon carrier landing control scheme with a variable reference trajectory based on the information of the glide rate is designed. Finally, simulations were performed to demonstrate the control effectiveness of the algorithm under different sea states, initial deviations, and reference trajectory selection strategies.
The simulation results show that the proposed algorithm has the following advantages: (1) The phase delay of the control system is greatly reduced by introducing prediction of deck motion into the reference trajectory selection. This also enhances the following ability of the aircraft and improves the landing accuracy significantly.
(2) Through the augmentation of the linear small disturbance model, the constraint problem of the control variable derivations in the control system is solved. This allows the control law generated by the algorithm to satisfy actual engineering requirements in real scenarios.
(3) The algorithm proposed in this paper can effectively control the landing trajectory of carrier-based aircraft with high accuracy in the case of initial deviation, external continuous wind disturbances, and random errors of state variables.
(4) The present algorithm has better performance comparing with the typical LQ algorithm, and the computational efficiency meets the requirements of real-time online tracking.
In this paper, the selection of the coefficient matrix, time window length, and receding step length were all obtained by experiments. Further, an intelligent optimization algorithm may be applied to further optimize these parameters so as to obtain better control effectiveness. And aiming at the carrier landing control problem with mismatched model term, some robust control methods [43], [44] can also serve as reference and guidance.