# Model Predictive Control for Automatic Carrier Landing with Time Delay

Authors: Kaikai Cui; Wei Han; Yujie Liu; Xinwei Wang; Xichao Su; Jie Liu
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

This paper focuses on the problem of automatic carrier landing control with time delay, and an antidelay model predictive control (AD-MPC) scheme for carrier landing based on the symplectic pseudospectral (SP) method and a prediction error method with particle swarm optimization (PE-PSO) is designed. Firstly, the mathematical model for carrier landing control with time delay is given, and based on the Padé approximation (PA) principle, the model with time delay is transformed into an equivalent nondelay one. Furthermore, a guidance trajectory based on the predicted trajectory shape and position deviation is designed in the MPC framework to eliminate the influence of carrier deck motion and real-time error. At the same time, a rolling optimal control block is designed based on the SP algorithm, in which the steady-state carrier air wake compensation is introduced to suppress the interference of the air wake. On this basis, the PE-PSO delay estimation algorithm is proposed to estimate the unknown delay parameter in the equivalent control model. The simulation results show that the delay estimation error of the PE-PSO algorithm is smaller than 2 ms, and the AD-MPC algorithm proposed in this paper can limit the landing height error within ±0.14 m under the condition of multiple disturbances and system input delay. The control accuracy of AD-MPC is much higher than that of the traditional pole assignment algorithm, and its computational efficiency meets the requirement of real-time online tracking.


---

## Body

## 1. Introduction

As an important symbol of national power, aircraft carriers play an indispensable role in maritime security [1]. As one of the key technologies in aircraft carrier systems, automatic carrier landing technology is of great importance to the navies of various countries. A well-designed automatic carrier landing system (ACLS) not only improves the landing accuracy of carrier-based aircraft but also reduces flight control difficulty and training costs for pilots. The small space on a carrier deck and significant marine environmental disturbances such as deck motion and carrier air wake impose severe limitations on landing performance [2]. According to the authors of [3], an ACLS can be divided into four layers for an inner loop, autopilot, guidance control, and guidance compensation. In [4], an ACLS was designed based on the H-infinity control technique to improve path control precision under the worst-case conditions of a vertical gust. Yu et al. proposed an active disturbance rejection control scheme for an ACLS in the final approach to achieve better tracking performance [5]. An ACLS design scheme with dynamic inversion techniques was also demonstrated to be promising and robust [6]. A stable adaptive control scheme was developed based on the LDU (lower-diagonal-upper) decomposition of the high-frequency gain matrix for carrier landing during the final approach [7]. Regarding the carrier landing problem for aircraft with short takeoff and landing capabilities, gain-scheduled linear optimal control and L-1 adaptive control were applied to an ACLS design in [8]. Additionally, several novel optimization algorithms have been proposed to optimize ACLS parameters and assist control algorithms to operate at full capacity [9, 10]. These optimization algorithms are mainly applied to controller parameter tuning. However, the influence of time delay in the control system is not considered in the above studies.Time delay is a difficult problem which cannot be ignored in an ACLS. During the process of automatic carrier landing, radar systems, shipboard computers, and flight control systems will introduce delay effects of different magnitudes for calculation or communication. Additionally, the response characteristics of the control actuator itself will also introduce delay effects. However, research on time delay in automatic carrier landing control is very sparse. The authors of [11] analyzed the sources and influences of delay in an ACLS, and a delay observer was designed to predict deck motion. However, the influence of control delay on ACLS was not considered in that article. Focusing on the manned aircraft carrier landing problem, Liu discussed the delayed control problem for carrier landing based on a pilot delay model [12]. In [13], the delay effects of the actuator and engine dynamics were eliminated via state space realization and state equation incorporation, but measurement signal delay and signal transmission delay were not considered. In general, these related studies were not comprehensive in terms of considering the delay.The main contributions for this paper are listed below: firstly, we proposed a model predictive control framework for the carrier landing control problem. Compared with other methods, the advantages of MPC are as follows: (1) In the longitudinal carrier landing trajectory control problem of an F/A-18, there are four control variables. The general control method often needs to consider the control allocation problem, while MPC can flexibly deal with the MIMO system with a coupling relationship. And the optimal control allocation scheme can be given according to the objective function. (2) The physical constraints usually exist in the controller actuator, and MPC can explicitly deal with the control problems with constraints. (3) The MPC control structure allows introducing model and environment prediction information into the control decision. So the predictable information of deck motion and the steady-state part of the ship air wake can be introduced into the current control decision-making, which is difficult to achieve by other methods. As a result, the influence of the ship air wake and deck motion on the control system can be well suppressed, and the landing control accuracy can be improved. Furthermore, there are few studies on the time delay control problem of carrier landing. Considering the time delay problem in the carrier landing progress, the MPC control framework is modified based on the Padé approximation principle, and to estimate the time delay, a PE-PSO delay estimate method is proposed in this paper. The modified MPC for carrier landing is suitable for time delay scenarios.This paper focuses on the problem of automatic carrier landing control with time delay, and the remainder of this paper is organized as follows. In Section2, the mathematical model for the carrier landing control with time delay is presented. In Section 3, an AD-MPC scheme for carrier landing based on the SP method and PE-PSO is designed. Experimental simulations and result analysis are presented in Section 4. Finally, Section 5 concludes this paper.


## 2. Mathematical Model for Carrier Landing with Time Delay

### 2.1. Small Disturbance Equation and Environment Model for Carrier Landing

In this paper, the following longitudinal linear small disturbance equation for an F/A-18A aircraft is presented to describe the carrier landing process [3, 13]:
(1)Ẋ=AMX+BMU+EMαg,X=ΔvV0,Δα,Δθ,Δq,ΔhV0T,U=ΔδH,ΔδLEF,ΔδRT,ΔδPLT, where Δv,Δα,Δθ,Δq, and Δh represent the velocity, angle of attack, pitch angle, pitch angle velocity, and altitude deviation relative to the nominal state, respectively. Additionally, ΔδH,ΔδLEF,ΔδRT, and ΔδPL represent deviations in terms of the horizontal tail deflection, leading-edge flap deflection, rudder toe-in deflection, and engine throttle control angle, respectively. Finally, αg represents the deviation of the attack angle caused by vertical wind disturbances. The nominal state for carrier landing is V0=69.96m/s, α0=8.3°, and γ0=−3°. In actual control systems, the following limits on the range and rate of change of the control surface deflection and throttle control angle must be considered:
(2)Uimin≤Ui≤UimaxU̇i≤U̇imax,i=1,2,⋯,4.A carrier is subjected to six-degree-of-freedom deck motions based on the wind and waves encountered at sea. In this study, the pitching and heaving motions of the carrier were taken into consideration based on their significant impact on landing accuracy. The deck motion models are defined as follows [14]:
(3)θSdeg=0.5sin0.6t+φ1+0.3sin0.63t+φ1+0.25,ZSft=4.0sin0.6t+φ2+1.0sin0.2t+φ2,where φ1 and φ2 represent random initial phases and θS and ZS represent the pitching and heaving motions of the carrier.The following engineering model for the carrier air wake from [14] is adopted in this paper:
(4)ug=u1+u2+u3+u4,vg=v1+v4,wg=w1+w2+w3+w4,where ug, vg, and wg represent the air wake in the x, y, and z directions, respectively, in the inertial frame. The subscript i (i=1, 2, 3, 4) represents atmospheric turbulence, the steady-state component, periodic component, and random component of the carrier air wake. Because the model in Equation (1) is only related to the air wake in the vertical direction, ug and vg were not considered in this study.


### 2.2. Sources of Time Delay in an ACLS

Time delays in an ACLS can be divided into four categories according to their sources, namely, measurement signal delay, control signal delay, signal transmission delay, and actuator response delay. Regarding their effects on a control system, the first three types of delay effects can be regarded as pure time delays, whereas the actuator delay is a type of high-order dynamic response delay. Various methods can be adopted to handle time delay effects. The time delays for each component of the F/A-18 ACLS are listed in Table1 [11].Table 1
Time delays of the F/A-18 ACLS.
Source of delayDelay time (ms)α−β filter of the guidance radar10Asynchronous delay of signal transmission[0, 50]Receiver decoding time1Asynchronous delay of the task computer[0, 50]Task computer computation delay[12, 50]Asynchronous delay of the flight control system[0, 50]Flight control system foreground calculation delay[6, 56]Flight control system calculation delay6.25In Table1, there is a certain degree of delay in each component of the ACLS. Although the delay effect is not significant for a single component, delay will have a significant impact on the ACLS when the delay times of all components are combined.Figure1 presents a schematic of data transmission in the ACLS. The effects of pure time delay can be regarded as the control action exerted by the control system at the current time T, which is obtained according to the state of the control system at T-θ.Figure 1
Data transmission in the ACLS.After the position and speed of the carrier aircraft at timeT are measured and filtered by the shipboard radar, the guidance system on the carrier calculates the corresponding information at time T+dt1. The guidance and control system calculates the control law according to the position, speed, and flight status information and obtains a corresponding control law at T+dt1+dt2. After a series of coding, transmission, and decoding operations, the actuators receive the control law information and begin to work at T+dt1+dt2+dt3, resulting in a pure time delay for the entire closed loop of τ=dt1+dt2+dt3. It is typically necessary to limit τ to less than 200 ms in practice.During the process of carrier landing, there will be a time delay between the control actuator (rudder and engine) receiving control signals and achieving the required position or state. Generally, the delay for a control surface is relatively small; it has little influence on control effects and can be ignored [13]. However, the response time of the engine is sufficiently slow to have a significant impact on control system performance. Considering the dynamic response delay of the engine, the throttle input variable δpL in the longitudinal landing model given in Equation (1) must be replaced with the thrust response variable δT.Delay in the control law will result in an aircraft being unable to track the ideal glide path accurately. This is because the control inputs for the control system cannot eliminate the current output errors in time based on the delayed signal, which increases adjustment time. The resulting overshoot causes the system to oscillate [15]. If the time delay is too large, then the control system cannot make correct control actions according to the current system state, which could lead to control failure. And the follow-up simulation results will also verify the above analysis.


### 2.3. Carrier Landing Control Model with Time Delay

For the F/A-18 aircraft, the following transfer function model for the dynamic relationship between the thrust response and the throttle input was given as follows:(5)δTδpL=as+bs2+cs+b=2.6710s+1.1846s2+2.5336s+1.1846.For this type of time delay effect, the state space model of the transfer function is typically constructed using the state space realization method from the linear system theory. To eliminate the influence of time delay, the constructed state space model is combined with the original model to obtain a new system model as follows:(6)ẋ9ẋ10=Mx9x10+NδPL=−c1−b0x9x10+abδPL,where x9=δT and x10 are intermediate variables. Additionally, to introduce the rate of change information for ΔδH,ΔδLEF, and ΔδRT into the system model explicitly, it is necessary to augment Equation (1) in the following manner:
(7)x¯̇=A¯x¯+B¯u¯+E¯αg,x¯=X,ΔδH,ΔδLEF,ΔδRT,ΔδT,x10T,u¯=Δδ̇H,Δδ̇LEF,Δδ̇RT,ΔδPLT,where
(8)A-=AMBM05×1⋯⋯⋯03×10⋯⋯⋯02×8MB-=05×4⋯⋯⋯I3×303×1⋯⋯⋯02×3NE-=EM⋯05×1The main influence of the pure time delay is that it will generate control signal input lag. The measurement times of each state variable are different, which will lead to the desynchronization of various measurement signals. This delay difference is typically very small, and the influence on the control system can be ignored. However, it is assumed that the flight control and guidance signals are both calculated by the shipboard guidance system. Therefore, all state variable information required for the control law generation process is transmitted to the carrier through a data link after airborne measurement, except for the altitude and velocity information, which are measured by the shipboard radar. For the convenience of processing, it is assumed that the measurement delays of all states are approximately equal and denoted asdt1.Considering the time delay in the control system, the longitudinal small disturbance mathematical model for carrier landing can be defined as follows:(9)x¯̇=A¯x¯+B¯u¯t−τ+E¯αg,x¯=X,ΔδH,ΔδLEF,ΔδRT,ΔδT,x10T,u¯=Δδ̇H,Δδ̇LEF,Δδ̇RT,ΔδPLT,where τ represents the pure delay of the entire ACLS.In this study, to handle the time delay problem, the first-order Padé approximation [16] was incorporated by introducing virtual time delay variables, and the original control system with time delay was transformed into a delay-free system. First, an input term with a delay u¯t−τ is transformed using the Laplace transform. Then, by applying the first-order Taylor approximation, we have
(10)Lu¯t−τ=e−τsLu¯t≈1−τs/21+τs/2Lu¯t,resulting in
(11)Lu¯t−τLu¯t≈1−τs/21+τs/2=4/τs+2/τ−1.The intermediate variablexdt is selected to satisfy
(12)ẋdt=4τu¯t−2τxdt.Then, the original system can be rewritten as(13)ẋ=x¯̇ẋd=Ax+Bu+Eαg=A¯B¯0A~x¯xd+−B¯B~u¯t+E¯0αg,A~=−2τI4,B~=4τI4.The carrier landing problem can be regarded as a trajectory tracking problem for an ideal glide path under the constraints of flight dynamics, control variables, and state variables. Based on the analysis above, the equivalent mathematical model for the optimal control problem of carrier landing with time delay is defined as(14)minJ=12∫t0t0+Tx−xeTPx−xe+uTRudt,s.t.ẋ=Ax+Bu+Eαg,xmin≤x≤xmax,umin≤u≤umax,where xe is the ideal flight state vector of the carrier-based aircraft, x is the actual flight state vector, and u is the actual control input. P and R are diagonal coefficient matrices, where P is a semipositive definite matrix and R is a positive definite matrix.


## 2.1. Small Disturbance Equation and Environment Model for Carrier Landing

In this paper, the following longitudinal linear small disturbance equation for an F/A-18A aircraft is presented to describe the carrier landing process [3, 13]:
(1)Ẋ=AMX+BMU+EMαg,X=ΔvV0,Δα,Δθ,Δq,ΔhV0T,U=ΔδH,ΔδLEF,ΔδRT,ΔδPLT, where Δv,Δα,Δθ,Δq, and Δh represent the velocity, angle of attack, pitch angle, pitch angle velocity, and altitude deviation relative to the nominal state, respectively. Additionally, ΔδH,ΔδLEF,ΔδRT, and ΔδPL represent deviations in terms of the horizontal tail deflection, leading-edge flap deflection, rudder toe-in deflection, and engine throttle control angle, respectively. Finally, αg represents the deviation of the attack angle caused by vertical wind disturbances. The nominal state for carrier landing is V0=69.96m/s, α0=8.3°, and γ0=−3°. In actual control systems, the following limits on the range and rate of change of the control surface deflection and throttle control angle must be considered:
(2)Uimin≤Ui≤UimaxU̇i≤U̇imax,i=1,2,⋯,4.A carrier is subjected to six-degree-of-freedom deck motions based on the wind and waves encountered at sea. In this study, the pitching and heaving motions of the carrier were taken into consideration based on their significant impact on landing accuracy. The deck motion models are defined as follows [14]:
(3)θSdeg=0.5sin0.6t+φ1+0.3sin0.63t+φ1+0.25,ZSft=4.0sin0.6t+φ2+1.0sin0.2t+φ2,where φ1 and φ2 represent random initial phases and θS and ZS represent the pitching and heaving motions of the carrier.The following engineering model for the carrier air wake from [14] is adopted in this paper:
(4)ug=u1+u2+u3+u4,vg=v1+v4,wg=w1+w2+w3+w4,where ug, vg, and wg represent the air wake in the x, y, and z directions, respectively, in the inertial frame. The subscript i (i=1, 2, 3, 4) represents atmospheric turbulence, the steady-state component, periodic component, and random component of the carrier air wake. Because the model in Equation (1) is only related to the air wake in the vertical direction, ug and vg were not considered in this study.


## 2.2. Sources of Time Delay in an ACLS

Time delays in an ACLS can be divided into four categories according to their sources, namely, measurement signal delay, control signal delay, signal transmission delay, and actuator response delay. Regarding their effects on a control system, the first three types of delay effects can be regarded as pure time delays, whereas the actuator delay is a type of high-order dynamic response delay. Various methods can be adopted to handle time delay effects. The time delays for each component of the F/A-18 ACLS are listed in Table1 [11].Table 1
Time delays of the F/A-18 ACLS.
Source of delayDelay time (ms)α−β filter of the guidance radar10Asynchronous delay of signal transmission[0, 50]Receiver decoding time1Asynchronous delay of the task computer[0, 50]Task computer computation delay[12, 50]Asynchronous delay of the flight control system[0, 50]Flight control system foreground calculation delay[6, 56]Flight control system calculation delay6.25In Table1, there is a certain degree of delay in each component of the ACLS. Although the delay effect is not significant for a single component, delay will have a significant impact on the ACLS when the delay times of all components are combined.Figure1 presents a schematic of data transmission in the ACLS. The effects of pure time delay can be regarded as the control action exerted by the control system at the current time T, which is obtained according to the state of the control system at T-θ.Figure 1
Data transmission in the ACLS.After the position and speed of the carrier aircraft at timeT are measured and filtered by the shipboard radar, the guidance system on the carrier calculates the corresponding information at time T+dt1. The guidance and control system calculates the control law according to the position, speed, and flight status information and obtains a corresponding control law at T+dt1+dt2. After a series of coding, transmission, and decoding operations, the actuators receive the control law information and begin to work at T+dt1+dt2+dt3, resulting in a pure time delay for the entire closed loop of τ=dt1+dt2+dt3. It is typically necessary to limit τ to less than 200 ms in practice.During the process of carrier landing, there will be a time delay between the control actuator (rudder and engine) receiving control signals and achieving the required position or state. Generally, the delay for a control surface is relatively small; it has little influence on control effects and can be ignored [13]. However, the response time of the engine is sufficiently slow to have a significant impact on control system performance. Considering the dynamic response delay of the engine, the throttle input variable δpL in the longitudinal landing model given in Equation (1) must be replaced with the thrust response variable δT.Delay in the control law will result in an aircraft being unable to track the ideal glide path accurately. This is because the control inputs for the control system cannot eliminate the current output errors in time based on the delayed signal, which increases adjustment time. The resulting overshoot causes the system to oscillate [15]. If the time delay is too large, then the control system cannot make correct control actions according to the current system state, which could lead to control failure. And the follow-up simulation results will also verify the above analysis.


## 2.3. Carrier Landing Control Model with Time Delay

For the F/A-18 aircraft, the following transfer function model for the dynamic relationship between the thrust response and the throttle input was given as follows:(5)δTδpL=as+bs2+cs+b=2.6710s+1.1846s2+2.5336s+1.1846.For this type of time delay effect, the state space model of the transfer function is typically constructed using the state space realization method from the linear system theory. To eliminate the influence of time delay, the constructed state space model is combined with the original model to obtain a new system model as follows:(6)ẋ9ẋ10=Mx9x10+NδPL=−c1−b0x9x10+abδPL,where x9=δT and x10 are intermediate variables. Additionally, to introduce the rate of change information for ΔδH,ΔδLEF, and ΔδRT into the system model explicitly, it is necessary to augment Equation (1) in the following manner:
(7)x¯̇=A¯x¯+B¯u¯+E¯αg,x¯=X,ΔδH,ΔδLEF,ΔδRT,ΔδT,x10T,u¯=Δδ̇H,Δδ̇LEF,Δδ̇RT,ΔδPLT,where
(8)A-=AMBM05×1⋯⋯⋯03×10⋯⋯⋯02×8MB-=05×4⋯⋯⋯I3×303×1⋯⋯⋯02×3NE-=EM⋯05×1The main influence of the pure time delay is that it will generate control signal input lag. The measurement times of each state variable are different, which will lead to the desynchronization of various measurement signals. This delay difference is typically very small, and the influence on the control system can be ignored. However, it is assumed that the flight control and guidance signals are both calculated by the shipboard guidance system. Therefore, all state variable information required for the control law generation process is transmitted to the carrier through a data link after airborne measurement, except for the altitude and velocity information, which are measured by the shipboard radar. For the convenience of processing, it is assumed that the measurement delays of all states are approximately equal and denoted asdt1.Considering the time delay in the control system, the longitudinal small disturbance mathematical model for carrier landing can be defined as follows:(9)x¯̇=A¯x¯+B¯u¯t−τ+E¯αg,x¯=X,ΔδH,ΔδLEF,ΔδRT,ΔδT,x10T,u¯=Δδ̇H,Δδ̇LEF,Δδ̇RT,ΔδPLT,where τ represents the pure delay of the entire ACLS.In this study, to handle the time delay problem, the first-order Padé approximation [16] was incorporated by introducing virtual time delay variables, and the original control system with time delay was transformed into a delay-free system. First, an input term with a delay u¯t−τ is transformed using the Laplace transform. Then, by applying the first-order Taylor approximation, we have
(10)Lu¯t−τ=e−τsLu¯t≈1−τs/21+τs/2Lu¯t,resulting in
(11)Lu¯t−τLu¯t≈1−τs/21+τs/2=4/τs+2/τ−1.The intermediate variablexdt is selected to satisfy
(12)ẋdt=4τu¯t−2τxdt.Then, the original system can be rewritten as(13)ẋ=x¯̇ẋd=Ax+Bu+Eαg=A¯B¯0A~x¯xd+−B¯B~u¯t+E¯0αg,A~=−2τI4,B~=4τI4.The carrier landing problem can be regarded as a trajectory tracking problem for an ideal glide path under the constraints of flight dynamics, control variables, and state variables. Based on the analysis above, the equivalent mathematical model for the optimal control problem of carrier landing with time delay is defined as(14)minJ=12∫t0t0+Tx−xeTPx−xe+uTRudt,s.t.ẋ=Ax+Bu+Eαg,xmin≤x≤xmax,umin≤u≤umax,where xe is the ideal flight state vector of the carrier-based aircraft, x is the actual flight state vector, and u is the actual control input. P and R are diagonal coefficient matrices, where P is a semipositive definite matrix and R is a positive definite matrix.


## 3. AD-MPC Algorithm for Carrier Landing with Time Delay

MPC is a closed-loop rolling optimization control method based on the mathematical model of the controlled plant [17]. The main idea of MPC is to solve the optimal control problem in a small limited period under the condition of control constraints. The optimal control input sequence is obtained through the optimization of the objective function, and the first input value of the control sequence is selected as the actual input of the next step. Then, the actual state of the system is introduced into the next calculation step as the initial condition. Repeat the above process until the controlling task is finished. The control block diagram can be summarized as Figure 2.Figure 2
Schematic of the MPC algorithm.


### 3.1. Guidance Trajectory Design Based on Trajectory Shape and Tracking Error

A crucial problem for MPC tracking control is selecting a standard tracking trajectoryxe, which is similar to the guidance law design for traditional tracking control. A well-designed tracking trajectory can make an aircraft track the desired glide accurately and eliminate tracking error rapidly. In this study, by combining trajectory deviations with the predicted trajectories, a guidance trajectory design method based on predicted trajectory shapes and position deviations was developed.This guidance trajectory design method is illustrated in Figure3. In Figure 3, A represents the position of the aircraft at the current time, and he represents its height deviation from the ideal position A0. B0 is the expected position of the aircraft after ΔT, and A0B0⌢ is the ideal trajectory obtained based on deck motion prediction. AB⌢ is the reference trajectory obtained by moving A0B0⌢ in a parallel direction and is used to generate the guidance trajectory AB0⌢. The points on AB0⌢ satisfy the relationship ht=he×l1/l1+l2. Therefore, the guidance trajectory AB0⌢ contains the prediction information for the future ideal trajectory and the error correction information for the current position simultaneously. To improve the adaptability of the guidance trajectory in different environments, the adjustment coefficients Kα and Kβ are introduced into the guidance trajectory, and the corresponding height difference between AB0⌢ and AB¯ is adjusted as follows:
(15)Δh=Kαhs+Kβhel1l1+l2.When the values ofKα and Kβ are both 1.0, the guidance trajectory is represented by AB0⌢ in Figure 3.Figure 3
Schematic diagram of the guidance trajectory design.To obtain future reference trajectory heights, prediction information based on deck motion is incorporated. To improve the efficiency of the algorithm, an autoregressive (AR) prediction algorithm is adopted to predict the deck motion of an aircraft carrier. The details of this method can be found in [18].


### 3.2. Trajectory Tracking Optimal Control Based on the Symplectic Pseudospectral Algorithm

In this paper, the symplectic pseudospectral algorithm based on the second type of generation function [19] is used as the rolling optimization method in the MPC framework. It is an efficient and accurate computational optimal control technique and has the capability of handling constraints on state and control variables [1]. Under general conditions, Equation (14) is modified as follows:
(16)minJ=12∫t0t0+Tx−xeTPx−xe+uTRudt,s.t.ẋ=Ax+Bu+w,h=Cx+Du+V≤0,where w=Eαg, and by introducing a nonnegative relaxation vector α, the inequality constraint in Equation (16) can be rewritten as follows:
(17)Cx+Du+V+α=0.Then, by introducing the costate vectorλ and Lagrangian multiplier μ, Equation (16) can be transformed into an unconstrained optimal control problem. The objective function can be expressed as
(18)J=∫t0t0+TH−λTẋdt,where H is the Hamilton function:
(19)Hx,u,λ,β,α=12x−xeTPx−xe+12uTRu+λTAx+Bu+w+μTCx+Du+V+α.According to the parametric variational principle, the following equations should be satisfied when the objective functionJ reaches its minimum value:
(20)∂H∂u=Ru+BTλ+DTμ=0,ẋ=∂H∂λ=Ax−BR−1BTλ+DTμ+w,λ̇=−∂H∂x=−Px−xe−ATλ−CTμ.According to the inequality constraints and the Karush-Kuhn-Tucker condition, we have(21)Cx−DR−1BTλ+DTμ+V+α=0,αTμ=0,μ≥0,α≥0.Equation (16) defines an optimal control problem with a fixed terminal time. Therefore, when the terminal state is free, we have λt0+T=0 according to the transversality condition.The variablesx, λ, μ, and α are approximated using the Nj-degree Legendre-Gauss-Lobatto (LGL) pseudospectral method in the jth subinterval as follows:
(22)xjτ=∑i=0Njxijρijτ,λjτ=∑i=0Njλijρijτ,μjτ=∑i=0Njμijρijτ,αjτ=∑i=0Njαijρijτ,where ρijτ is the Lagrangian interpolation polynomial corresponding to the ith LGL node within the jth subinterval [19].The stagnation point condition for the second-generation function [20, 21] is applied to the jth subinterval, and we have
(23)Kjxj−1x¯jλ¯jλj+ξxjξλjμ∧j+ζxjζλj=λj−10dNj,10dNj,1xj.According to the constraint in Equation (21) and the complementarity condition, we obtain the following relationship for each subinterval Γj:
(24)Gjx∧j−Hjλ∧j−Mjμ∧j+v∧j+α∧j=0,α∧j≥0,μ∧j≥0,μ∧jα∧j=0.By assembling Equations (23) and (24) in every interval according to the boundary conditions, the two-point boundary value problem (TPBVP) in t0,t0+T can be rewritten as follows:
(25)a:Kσx^,λ^+ξμ^+ζ=r,b:Gx^−Hλ^−Mμ^+v^+α^=0,c:α^≥0,μ^≥0,μ∧Tα^=0,where x^, λ^, and μ^ contain the information regarding the state vectors, costate vectors, and Lagrange multipliers at the LGL collocation points, respectively. According to Equation (25)(a), the variables x^ and λ^ can be expressed as functions of μ^. Therefore, the TPBVP given in Equation (25) can be transformed into a standard linear complementarity problem based on Equations (25)(b) and (25)(c), as shown in Equation (26). Additionally, μ^ can be derived using the Lemke method [22], and x^ and λ^ can be derived according to Equation (25)(a). Finally, the control variable u can be obtained by substituting x^ and λ^ into Equation (20).
(26)Yμ^+q≥0,μ^≥0,μ∧TYμ^+q=0.The specific expressions of the variables involved in Equations (23) to (26) and the detailed derivation process of the symplectic pseudospectral algorithm based on the second type of generation function can be found in the appendix.It should be note that the atmospheric turbulence, periodic component, and random component of the carrier air wake in Equation (4) contain random factors. That means their strength information cannot be obtained in advance. However, the steady-state component w2 is usually only related to the strength of the deck wind, so this part information of the air wake field can be introduced into the rolling optimization process, that is, to let w=Eα2 in Equation (19), where α2=−w2/V0.


### 3.3. PE-PSO Estimation Algorithm for Pure Time Delay

The pure time delayτ in Equation (13) can be determined by experience. However, influenced by various factors, the time delay may vary in practical application. Therefore, it is necessary to estimate the time delay of the closed-loop system according to its input and output information. This problem lies in the category of parameter estimation, but traditional parameter estimation algorithms are largely based on the AR moving average model. In this study, system inputs are obtained based on the MPC control algorithm, and their probability distribution information cannot be obtained. Therefore, PE-PSO is utilized to identify the closed-loop time delay of the control system.In simulations and practical engineering, data sampling and calculations are conducted in discrete forms. Therefore, the following discrete prediction system is constructed to identify the time delay parameterτ:
(27)x¯^k+1=Gx¯^k+Hu¯k−τ^ΔT,G=eA¯ΔT,H=∫0ΔTeA¯ΔTdtB¯,where
(28)x¯k+1=Gx¯k+Nαgk,k∈N+&k<d,x¯k+1=Gx¯k+Hu¯k−d+Nαgk,k∈N+&k≥d,x¯k=xk,ΔδHk,ΔδLEFk,ΔδRTk,ΔδTk,x10kT,u¯k=Δδ̇Hk,Δδ̇LEFk,Δδ̇RTk,ΔδPLkT,N=∫0ΔTeA¯ΔTdtE¯.Here,x¯k and u¯k represent the values of x¯(kΔT) and u¯(kΔT), respectively. When d^=τ^/ΔT or d=τ/ΔT is a noninteger, u¯k−τ^/ΔT and u¯k−d are obtained via linear interpolation, where x¯^ and τ^ represent the estimators of x¯ and τ, respectively.The PE-PSO method is used to estimate the time delay variableτ, based on the deviations between the state outputs of the reconstructed and actual systems. Because the dummy variable x10 in x¯ cannot be measured, the observable state variable y in the system should be predicted according to the following prediction model to obtain its estimated value y^ at kΔT:
(29)y^k=fYk−1,U¯k,k,τ^=CMx¯^k=CMGk−k0xk0+∑i=k0k−k0−1Gk−i−1Hu¯i−τ^ΔT,where CM is the observation matrix, and the actual observable value yk can be expressed as
(30)yk=CMGk−k0xk0+∑i=k0k−k0−1Gk−i−1Hu¯i−τΔT+Nαgi+ek=CMGk−k0xk0+∑i=k0k−k0−1Gk−i−1Hu¯i−τΔT+∑i=k0k−k0−1Gk−i−1Nαgi+ek=fYk−1,U¯k,k,τ+εk,where εk=∑i=k0k−k0−1Gk−i−1Nαgi+ek. Yk−1 represents all historical datasets of observable variables y¯ before kΔT, and U¯k represents the datasets of u¯ at kΔT and before kΔT. The prediction error at kΔT is denoted as
(31)vkτ^=yk−y^k.When comparing the prediction model to the actual dynamic model, it can be seen thatεk is mainly composed of two parts: interference caused by atmospheric turbulence and measurement error in observable variables. Generally, measurement error can be regarded as white noise. When a carrier aircraft is located outside the influence area of the carrier air wake, atmospheric disturbance can also be regarded as white noise. Therefore, εk approximately obeys a zero-mean normal distribution. Additionally, based on the small time scale of the final glide process, the delay τ is not expected to change significantly during this period. Assuming that the observation sequence contains N observation values, the prediction error sequence can be expressed as vkτ^,k=k1,k2,⋯kN, and the estimated value of the delay variable τ is calculated as
(32)τ^=argminJτ^=argmin1N∑k=k1k=kNyk−y∧kTyk−y^k.Theorem 1.
When the sequenceεk follows a normal distribution with zero mean and the components are independent, the estimated value τ^ obtained by Equation (32) is strongly consistent when the objective function Jτ^ is minimized.This can be proven as follows:(33)Jτ^=1N∑k=k1k=kNyk−y∧kTyk−y^k=1N∑k=k1k=kNyk−fYk−1,U¯k,k,τ∧Txk−fYk−1,U¯k,k,τ^=1N∑k=k1k=kNεk+Δfτ∧,τTεk+Δfτ^,τ=1N∑k=k1k=kNεTkεk+1N∑k=k1k=kNεTkΔfτ^,τ+1N∑k=k1k=kNΔfTτ^,τεk+1N∑k=k1k=kNΔfTτ^,τΔfτ^,τ,where
(34)Δfτ^,τ=fYk−1,U¯k,k,τ−fYk−1,U¯k,k,τ^.Assuming that the statistical characteristics of each stochastic process do not change over time (i.e., stationary stochastic process), whenN⟶∞, each item in the formula above converges to its mean value with a probability of one.
(35)1N∑k=k1k=kNεTkεk⟶a.s.EεTkεk=E∑i=k1k−k1−1Gk−i−1Nαgi+ekT∑i=k1k−k1−1Gk−i−1Nαgi+ek=E∑i=k1k−k1−1∑j=k1k−k1−1CijαgTiαgj+2eTk∑i=k1k−k1−1Gk−i−1Nαgi+eTkek=∑i=k1k−k1−1∑j=k1k−k1−1CijRααi−j+∑i=1nσi2=C,where
(36)Cij=NTGk−i−1TGk−j−1N,Rααi,j=EαgTiαgj=Rααi−j.Here,Rααi,j and σi2 represent the cross-correlation function of atmospheric turbulence at different times and variance of state variable measurement errors, respectively. Under the assumption of a stationary random process, C is constant. Because εk only contains information regarding atmospheric disturbances and current measurement errors, its value is not affected by historical measurement data, the initial state, control inputs, or delay values. Therefore, εk and Δfτ^,τ are independent of each other, implying
(37)1N∑k=k1k=kNΔfTτ^,τεk⟶a.s.EΔfTτ^,τεk=EΔfTτ^,τEεk=0.Similarly,(38)1N∑k=k1k=kNΔfTτ^,τεk⟶a.s.EΔfTτ^,τεk=0,1N∑k=k1k=kNΔfTτ^,τΔfτ^,τ⟶a.s.EΔfTτ^,τΔfτ^,τ.WhenN⟶∞,
(39)Jτ^N⟶a.s.C+EΔfTτ^N,τΔfτ^N,τ.WhenJτ^ is taken as the minimum value, the following equation must be satisfied [23]:
(40)EΔfTτ^N,τΔfτ^N,τ=0.According to the definition ofΔfτ^N,τ, if τ^≠τ, then Δfτ^N,τ≠0 holds with a probability of one. Therefore, the above formula is equivalent to τ^N⟶a.s.τ and Theorem 1 is proven.To obtainτ^ that satisfies Equation (32), we transfer the problem of extremum determining to an optimation problem as follows:
(41)ϒ:minJτ^=1N∑k=k1k=kNyk−y∧kTyk−y^k,s.t.τ^∈R+,0≤τ^≤τmax,where τmax represents the maximum possible value of τ, which can be determined by practical experience. Then, the classical PSO algorithm is used to solve optimation problem ϒ, and the optimal estimation result τ^ of τ is obtained.The process of approaching and landing on an aircraft carrier is typically divided into two stages. In stage 1, the influence of deck motion is not considered, and the tracking trajectory is a straight-line glide path in space. Under the small disturbance linear model, the control process can be regarded as zero-trajectory control. In stage 2 (i.e., approximately 20 s before carrier deck contact), deck motion is introduced. The aircraft must adjust its height according to the motion of the ideal landing point. In stage 1, the carrier aircraft is still outside the influence range of the carrier air wake, and the atmospheric turbulence disturbance on the aircraft can be regarded as white noise. Therefore, the accuracy of time delay parameter identification in this stage has theoretical support.A flowchart for the carrier landing control algorithm with time delay developed in this study is presented in Figure4, where a dashed line or dotted box indicates that the corresponding step is only performed once.Figure 4
Flowchart for the carrier landing control algorithm with time delay.


## 3.1. Guidance Trajectory Design Based on Trajectory Shape and Tracking Error

A crucial problem for MPC tracking control is selecting a standard tracking trajectoryxe, which is similar to the guidance law design for traditional tracking control. A well-designed tracking trajectory can make an aircraft track the desired glide accurately and eliminate tracking error rapidly. In this study, by combining trajectory deviations with the predicted trajectories, a guidance trajectory design method based on predicted trajectory shapes and position deviations was developed.This guidance trajectory design method is illustrated in Figure3. In Figure 3, A represents the position of the aircraft at the current time, and he represents its height deviation from the ideal position A0. B0 is the expected position of the aircraft after ΔT, and A0B0⌢ is the ideal trajectory obtained based on deck motion prediction. AB⌢ is the reference trajectory obtained by moving A0B0⌢ in a parallel direction and is used to generate the guidance trajectory AB0⌢. The points on AB0⌢ satisfy the relationship ht=he×l1/l1+l2. Therefore, the guidance trajectory AB0⌢ contains the prediction information for the future ideal trajectory and the error correction information for the current position simultaneously. To improve the adaptability of the guidance trajectory in different environments, the adjustment coefficients Kα and Kβ are introduced into the guidance trajectory, and the corresponding height difference between AB0⌢ and AB¯ is adjusted as follows:
(15)Δh=Kαhs+Kβhel1l1+l2.When the values ofKα and Kβ are both 1.0, the guidance trajectory is represented by AB0⌢ in Figure 3.Figure 3
Schematic diagram of the guidance trajectory design.To obtain future reference trajectory heights, prediction information based on deck motion is incorporated. To improve the efficiency of the algorithm, an autoregressive (AR) prediction algorithm is adopted to predict the deck motion of an aircraft carrier. The details of this method can be found in [18].


## 3.2. Trajectory Tracking Optimal Control Based on the Symplectic Pseudospectral Algorithm

In this paper, the symplectic pseudospectral algorithm based on the second type of generation function [19] is used as the rolling optimization method in the MPC framework. It is an efficient and accurate computational optimal control technique and has the capability of handling constraints on state and control variables [1]. Under general conditions, Equation (14) is modified as follows:
(16)minJ=12∫t0t0+Tx−xeTPx−xe+uTRudt,s.t.ẋ=Ax+Bu+w,h=Cx+Du+V≤0,where w=Eαg, and by introducing a nonnegative relaxation vector α, the inequality constraint in Equation (16) can be rewritten as follows:
(17)Cx+Du+V+α=0.Then, by introducing the costate vectorλ and Lagrangian multiplier μ, Equation (16) can be transformed into an unconstrained optimal control problem. The objective function can be expressed as
(18)J=∫t0t0+TH−λTẋdt,where H is the Hamilton function:
(19)Hx,u,λ,β,α=12x−xeTPx−xe+12uTRu+λTAx+Bu+w+μTCx+Du+V+α.According to the parametric variational principle, the following equations should be satisfied when the objective functionJ reaches its minimum value:
(20)∂H∂u=Ru+BTλ+DTμ=0,ẋ=∂H∂λ=Ax−BR−1BTλ+DTμ+w,λ̇=−∂H∂x=−Px−xe−ATλ−CTμ.According to the inequality constraints and the Karush-Kuhn-Tucker condition, we have(21)Cx−DR−1BTλ+DTμ+V+α=0,αTμ=0,μ≥0,α≥0.Equation (16) defines an optimal control problem with a fixed terminal time. Therefore, when the terminal state is free, we have λt0+T=0 according to the transversality condition.The variablesx, λ, μ, and α are approximated using the Nj-degree Legendre-Gauss-Lobatto (LGL) pseudospectral method in the jth subinterval as follows:
(22)xjτ=∑i=0Njxijρijτ,λjτ=∑i=0Njλijρijτ,μjτ=∑i=0Njμijρijτ,αjτ=∑i=0Njαijρijτ,where ρijτ is the Lagrangian interpolation polynomial corresponding to the ith LGL node within the jth subinterval [19].The stagnation point condition for the second-generation function [20, 21] is applied to the jth subinterval, and we have
(23)Kjxj−1x¯jλ¯jλj+ξxjξλjμ∧j+ζxjζλj=λj−10dNj,10dNj,1xj.According to the constraint in Equation (21) and the complementarity condition, we obtain the following relationship for each subinterval Γj:
(24)Gjx∧j−Hjλ∧j−Mjμ∧j+v∧j+α∧j=0,α∧j≥0,μ∧j≥0,μ∧jα∧j=0.By assembling Equations (23) and (24) in every interval according to the boundary conditions, the two-point boundary value problem (TPBVP) in t0,t0+T can be rewritten as follows:
(25)a:Kσx^,λ^+ξμ^+ζ=r,b:Gx^−Hλ^−Mμ^+v^+α^=0,c:α^≥0,μ^≥0,μ∧Tα^=0,where x^, λ^, and μ^ contain the information regarding the state vectors, costate vectors, and Lagrange multipliers at the LGL collocation points, respectively. According to Equation (25)(a), the variables x^ and λ^ can be expressed as functions of μ^. Therefore, the TPBVP given in Equation (25) can be transformed into a standard linear complementarity problem based on Equations (25)(b) and (25)(c), as shown in Equation (26). Additionally, μ^ can be derived using the Lemke method [22], and x^ and λ^ can be derived according to Equation (25)(a). Finally, the control variable u can be obtained by substituting x^ and λ^ into Equation (20).
(26)Yμ^+q≥0,μ^≥0,μ∧TYμ^+q=0.The specific expressions of the variables involved in Equations (23) to (26) and the detailed derivation process of the symplectic pseudospectral algorithm based on the second type of generation function can be found in the appendix.It should be note that the atmospheric turbulence, periodic component, and random component of the carrier air wake in Equation (4) contain random factors. That means their strength information cannot be obtained in advance. However, the steady-state component w2 is usually only related to the strength of the deck wind, so this part information of the air wake field can be introduced into the rolling optimization process, that is, to let w=Eα2 in Equation (19), where α2=−w2/V0.


## 3.3. PE-PSO Estimation Algorithm for Pure Time Delay

The pure time delayτ in Equation (13) can be determined by experience. However, influenced by various factors, the time delay may vary in practical application. Therefore, it is necessary to estimate the time delay of the closed-loop system according to its input and output information. This problem lies in the category of parameter estimation, but traditional parameter estimation algorithms are largely based on the AR moving average model. In this study, system inputs are obtained based on the MPC control algorithm, and their probability distribution information cannot be obtained. Therefore, PE-PSO is utilized to identify the closed-loop time delay of the control system.In simulations and practical engineering, data sampling and calculations are conducted in discrete forms. Therefore, the following discrete prediction system is constructed to identify the time delay parameterτ:
(27)x¯^k+1=Gx¯^k+Hu¯k−τ^ΔT,G=eA¯ΔT,H=∫0ΔTeA¯ΔTdtB¯,where
(28)x¯k+1=Gx¯k+Nαgk,k∈N+&k<d,x¯k+1=Gx¯k+Hu¯k−d+Nαgk,k∈N+&k≥d,x¯k=xk,ΔδHk,ΔδLEFk,ΔδRTk,ΔδTk,x10kT,u¯k=Δδ̇Hk,Δδ̇LEFk,Δδ̇RTk,ΔδPLkT,N=∫0ΔTeA¯ΔTdtE¯.Here,x¯k and u¯k represent the values of x¯(kΔT) and u¯(kΔT), respectively. When d^=τ^/ΔT or d=τ/ΔT is a noninteger, u¯k−τ^/ΔT and u¯k−d are obtained via linear interpolation, where x¯^ and τ^ represent the estimators of x¯ and τ, respectively.The PE-PSO method is used to estimate the time delay variableτ, based on the deviations between the state outputs of the reconstructed and actual systems. Because the dummy variable x10 in x¯ cannot be measured, the observable state variable y in the system should be predicted according to the following prediction model to obtain its estimated value y^ at kΔT:
(29)y^k=fYk−1,U¯k,k,τ^=CMx¯^k=CMGk−k0xk0+∑i=k0k−k0−1Gk−i−1Hu¯i−τ^ΔT,where CM is the observation matrix, and the actual observable value yk can be expressed as
(30)yk=CMGk−k0xk0+∑i=k0k−k0−1Gk−i−1Hu¯i−τΔT+Nαgi+ek=CMGk−k0xk0+∑i=k0k−k0−1Gk−i−1Hu¯i−τΔT+∑i=k0k−k0−1Gk−i−1Nαgi+ek=fYk−1,U¯k,k,τ+εk,where εk=∑i=k0k−k0−1Gk−i−1Nαgi+ek. Yk−1 represents all historical datasets of observable variables y¯ before kΔT, and U¯k represents the datasets of u¯ at kΔT and before kΔT. The prediction error at kΔT is denoted as
(31)vkτ^=yk−y^k.When comparing the prediction model to the actual dynamic model, it can be seen thatεk is mainly composed of two parts: interference caused by atmospheric turbulence and measurement error in observable variables. Generally, measurement error can be regarded as white noise. When a carrier aircraft is located outside the influence area of the carrier air wake, atmospheric disturbance can also be regarded as white noise. Therefore, εk approximately obeys a zero-mean normal distribution. Additionally, based on the small time scale of the final glide process, the delay τ is not expected to change significantly during this period. Assuming that the observation sequence contains N observation values, the prediction error sequence can be expressed as vkτ^,k=k1,k2,⋯kN, and the estimated value of the delay variable τ is calculated as
(32)τ^=argminJτ^=argmin1N∑k=k1k=kNyk−y∧kTyk−y^k.Theorem 1.
When the sequenceεk follows a normal distribution with zero mean and the components are independent, the estimated value τ^ obtained by Equation (32) is strongly consistent when the objective function Jτ^ is minimized.This can be proven as follows:(33)Jτ^=1N∑k=k1k=kNyk−y∧kTyk−y^k=1N∑k=k1k=kNyk−fYk−1,U¯k,k,τ∧Txk−fYk−1,U¯k,k,τ^=1N∑k=k1k=kNεk+Δfτ∧,τTεk+Δfτ^,τ=1N∑k=k1k=kNεTkεk+1N∑k=k1k=kNεTkΔfτ^,τ+1N∑k=k1k=kNΔfTτ^,τεk+1N∑k=k1k=kNΔfTτ^,τΔfτ^,τ,where
(34)Δfτ^,τ=fYk−1,U¯k,k,τ−fYk−1,U¯k,k,τ^.Assuming that the statistical characteristics of each stochastic process do not change over time (i.e., stationary stochastic process), whenN⟶∞, each item in the formula above converges to its mean value with a probability of one.
(35)1N∑k=k1k=kNεTkεk⟶a.s.EεTkεk=E∑i=k1k−k1−1Gk−i−1Nαgi+ekT∑i=k1k−k1−1Gk−i−1Nαgi+ek=E∑i=k1k−k1−1∑j=k1k−k1−1CijαgTiαgj+2eTk∑i=k1k−k1−1Gk−i−1Nαgi+eTkek=∑i=k1k−k1−1∑j=k1k−k1−1CijRααi−j+∑i=1nσi2=C,where
(36)Cij=NTGk−i−1TGk−j−1N,Rααi,j=EαgTiαgj=Rααi−j.Here,Rααi,j and σi2 represent the cross-correlation function of atmospheric turbulence at different times and variance of state variable measurement errors, respectively. Under the assumption of a stationary random process, C is constant. Because εk only contains information regarding atmospheric disturbances and current measurement errors, its value is not affected by historical measurement data, the initial state, control inputs, or delay values. Therefore, εk and Δfτ^,τ are independent of each other, implying
(37)1N∑k=k1k=kNΔfTτ^,τεk⟶a.s.EΔfTτ^,τεk=EΔfTτ^,τEεk=0.Similarly,(38)1N∑k=k1k=kNΔfTτ^,τεk⟶a.s.EΔfTτ^,τεk=0,1N∑k=k1k=kNΔfTτ^,τΔfτ^,τ⟶a.s.EΔfTτ^,τΔfτ^,τ.WhenN⟶∞,
(39)Jτ^N⟶a.s.C+EΔfTτ^N,τΔfτ^N,τ.WhenJτ^ is taken as the minimum value, the following equation must be satisfied [23]:
(40)EΔfTτ^N,τΔfτ^N,τ=0.According to the definition ofΔfτ^N,τ, if τ^≠τ, then Δfτ^N,τ≠0 holds with a probability of one. Therefore, the above formula is equivalent to τ^N⟶a.s.τ and Theorem 1 is proven.To obtainτ^ that satisfies Equation (32), we transfer the problem of extremum determining to an optimation problem as follows:
(41)ϒ:minJτ^=1N∑k=k1k=kNyk−y∧kTyk−y^k,s.t.τ^∈R+,0≤τ^≤τmax,where τmax represents the maximum possible value of τ, which can be determined by practical experience. Then, the classical PSO algorithm is used to solve optimation problem ϒ, and the optimal estimation result τ^ of τ is obtained.The process of approaching and landing on an aircraft carrier is typically divided into two stages. In stage 1, the influence of deck motion is not considered, and the tracking trajectory is a straight-line glide path in space. Under the small disturbance linear model, the control process can be regarded as zero-trajectory control. In stage 2 (i.e., approximately 20 s before carrier deck contact), deck motion is introduced. The aircraft must adjust its height according to the motion of the ideal landing point. In stage 1, the carrier aircraft is still outside the influence range of the carrier air wake, and the atmospheric turbulence disturbance on the aircraft can be regarded as white noise. Therefore, the accuracy of time delay parameter identification in this stage has theoretical support.A flowchart for the carrier landing control algorithm with time delay developed in this study is presented in Figure4, where a dashed line or dotted box indicates that the corresponding step is only performed once.Figure 4
Flowchart for the carrier landing control algorithm with time delay.


## 4. Simulation Experiments and Analysis

### 4.1. Simulation Results for AD-MPC Carrier Landing Control Algorithm

The matrices of the longitudinal linear small disturbance equation in Equation (1) are listed as follows:
(42)AM=‐0.07050.0475‐0.14030‐5.8×10−5‐0.3110‐0.343000.99131.02×10−3000100.0218‐1.16600‐0.254400‐1100,BM=0.01210.002480.16900.2316‐0.07210.01400.0128‐0.03380000‐1.8150‐0.07900.16810.00230000,EM=0.0475‐0.3430‐1.1660T.The initial state of the aircraftx0=10−2×1−0.170.35−0.17−0.2901×9T, which means the initial error of Δh is -0.2 m. The control parameter matrices are P=103×221510001×9T and R=20201020T. The control parameters are Kα=0.1 and Kβ=1. In the PE-PSO algorithm, the number of particles was set to 20, and the maximum number of iterations was set to 30. In the time delay estimation stage, the time delay in the control system in Equation (13) is set to 200 ms by default. All simulation experiments are performed using MATLAB R2016b on a Personal Computer with a 2.3 GHz CPU and 8 G of RAM.Firstly, assume that the pure time delay of the landing systemτ is 200 ms and the proposed AD-MPC carrier landing control method is applied. Figure 5 presents the simulation results. The glide path tracking result of aircraft is shown in Figure 5(a), where the dotted line indicates the time when the deck motion is introduced. It can be analyzed from Figure 5(a) that the AD-MPC carrier landing control method can accomplish accurate tracking of the ideal glide path under the conditions of carrier air wake, deck motion, random interference, and time delay. Then, the simulation is calculated for 10 times, and the average calculation time of a single step is 35 ms, which is significantly less than the simulation step of 50 ms. It indicates that the AD-MPC algorithm meet the requirements of online tracking calculation efficiency.Figure 5
Carrier landing control effectiveness.
(a)
Glide path tracking result(b)
Deck motion tracking result(c)
Height error of the aircraft(d)
Delay estimation resultFurther, the deck motion tracking result is given in Figure5(b). It is shown in Figure 5(b) that the control scheme can eliminate the initial error in about 3 seconds when the initial height error is -0.2 m and the aircraft can accurately track the deck motion within the last 15 seconds before landing. The height error of aircraft is shown in Figure 5(c), which demonstrates that within the last 15 seconds before landing, the maximum height tracking error of the aircraft is only 0.13 m, which is far less than the allowable height error of 1.5 m during the carrier landing process [24]. Meanwhile, the delay estimation results of PE-PSO are given in Figure 5(d). It can be seen that the PE-PSO algorithm can estimate the time delay in ACLS system effectively, and the final estimation error is less than 1 ms.The state variables and control variables are presented in Figure6. One can see that the control variables obtained using the AD-MPC carrier landing method are relatively stable, and there is no significant jitter. Additionally, the control variables are all constrained by the performance range of the corresponding actuator, which ensures the feasibility of the control law.Figure 6
Results of the state and control variables.


### 4.2. Control Performance Comparison

In order to verify the control performance of the AD-MPC carrier landing control algorithm, the traditional pole assignment method is selected for comparison. To get better comparison result, the initial value is increased to -0.35 m, and the initial conditionx0 is set as 10−2×0.43−0.17−0.30−0.501×9T.Firstly, setτ as 100 ms and the deck motion tracking results of the two algorithms are given in Figure 7. It is shown in Figure 7(a) that the aircraft can accomplish the glide path tracking under the control of pole assignment method when the deck motion is not introduced into the control system. However, there will be a significant phase lag on the aircraft height tracking result after the introduction of deck motion. Further, Figure 7(b) demonstrates that in the initial control stage, the response speed of pole assignment method is slightly faster than that of AD-MPC algorithm, but during the last 15 seconds before landing, the altitude tracking error range of pole assignment method is about ±0.55 m, while this error range of AD-MPC is ±0.13 m, which is only 1/4 of that of the pole assignment method.Figure 7
Carrier landing control effectiveness.
(a)
Deck motion tracking results(b)
Height errors of the aircraftFurthermore, when the time delay is set as 150 ms, the comparison results of the two algorithms are given in Figure8. As can be seen from Figure 8(a), due to the increase of time delay, the pole assignment method cannot effectively control the height of aircraft even if the deck motion is not considered. The height of aircraft will fluctuate continuously. After the introduction of deck motion, the phase lag still exists, and the tracking performance is worse than that of Figure 7(a) (τ=100ms). Further, it can be seen from Figure 8(b) that within the last 15 seconds before landing, the height error range of pole placement method is about ±0.6 m, while this error range of AD-MPC is ±0.14 m. When the time delay is increased to 200 ms, the pole assignment method will diverge, and it can be seen from Figure 5 that AD-MPC algorithm can still effectively control the aircraft.Figure 8
Carrier landing control effectiveness.
(a)
Deck motion tracking results(b)
Height errors of the aircraftThe above analysis shows that the AD-MPC carrier landing control method designed in this paper has stronger antidelay capability and higher control accuracy compared with the traditional pole assignment method.


### 4.3. Simulation Results for Time Delay Estimation Algorithm

The time delay estimation results of the two examples in Section4.2 are shown in Figures 9(a) and 9(b), respectively. It can be seen from Figure 9 that the estimation errors of the two cases are both less than 2 ms, which shows that the PE-PSO algorithm can estimate the time delay accurately.Figure 9
Time delay estimation results for the ACLS.
(a)
τ=100ms(b)
τ=150msTo further verify the effectiveness of the PE-PSO time delay estimation algorithm, control delays of 50, 100, 150, and 200 ms are introduced into the landing control system successively. And for comparison, the artificial fish school algorithm (AFS) and pigeon-inspired optimization (PIO) algorithm are selected to replace PSO in PE-PSO. Considering the four time delay scenarios, Monte Carlo simulations are carried out 500 times for every situation. The numbers of population sizes are all set as 20. And the statistical results are show in Figure10 and Table 2.Figure 10
Statistics of Monte Carlo simulation results.
(a)(b)(c)Table 2
Statistics of Monte Carlo simulation results.
Real time delay (ms)RMSE (ms)Mean fitness function (10-7)Mean CPU time (ms)AFSPIOPSOAFSPIOPSOAFSPIOPSO501.97291.97311.97718.40338.40338.4038551.83340.15147.431002.08682.08692.08758.55298.55288.5530565.67352.86148.081501.23721.24111.24378.96808.96768.9683541.63325.35137.822000.86660.86540.86229.23049.22989.2319614.15328.56141.37It can be seen from the statistical results that the PIO algorithm yields the smallest mean fitness function value (MFFV) among the three methods. The RMSE of AFS is the minimal except when time delay is 200 ms. The RMSE and the MFFV of the three methods are quite close to each other (the largest difference is less than 0.52 %), while the CPU time of AFS is almost 4 times of that of PSO for the unique trial mechanism of AFS. Meanwhile, the CPU time of PIO is about 2.3 times of that of PSO.In fact, there is only one variable to be optimized in time delay parameter estimation, and the optimization problem is not complicated. Comparing with the absolute accuracy of parameter estimation, we are more concerned about how to obtain the enough accurate delay estimate result in the shortest time. Shorter time consumption means faster updating frequency of the delay estimate result. And from the Monte Carlo simulations, it can be seen that the PE-PSO proposed in this paper meets this design requirement.


## 4.1. Simulation Results for AD-MPC Carrier Landing Control Algorithm

The matrices of the longitudinal linear small disturbance equation in Equation (1) are listed as follows:
(42)AM=‐0.07050.0475‐0.14030‐5.8×10−5‐0.3110‐0.343000.99131.02×10−3000100.0218‐1.16600‐0.254400‐1100,BM=0.01210.002480.16900.2316‐0.07210.01400.0128‐0.03380000‐1.8150‐0.07900.16810.00230000,EM=0.0475‐0.3430‐1.1660T.The initial state of the aircraftx0=10−2×1−0.170.35−0.17−0.2901×9T, which means the initial error of Δh is -0.2 m. The control parameter matrices are P=103×221510001×9T and R=20201020T. The control parameters are Kα=0.1 and Kβ=1. In the PE-PSO algorithm, the number of particles was set to 20, and the maximum number of iterations was set to 30. In the time delay estimation stage, the time delay in the control system in Equation (13) is set to 200 ms by default. All simulation experiments are performed using MATLAB R2016b on a Personal Computer with a 2.3 GHz CPU and 8 G of RAM.Firstly, assume that the pure time delay of the landing systemτ is 200 ms and the proposed AD-MPC carrier landing control method is applied. Figure 5 presents the simulation results. The glide path tracking result of aircraft is shown in Figure 5(a), where the dotted line indicates the time when the deck motion is introduced. It can be analyzed from Figure 5(a) that the AD-MPC carrier landing control method can accomplish accurate tracking of the ideal glide path under the conditions of carrier air wake, deck motion, random interference, and time delay. Then, the simulation is calculated for 10 times, and the average calculation time of a single step is 35 ms, which is significantly less than the simulation step of 50 ms. It indicates that the AD-MPC algorithm meet the requirements of online tracking calculation efficiency.Figure 5
Carrier landing control effectiveness.
(a)
Glide path tracking result(b)
Deck motion tracking result(c)
Height error of the aircraft(d)
Delay estimation resultFurther, the deck motion tracking result is given in Figure5(b). It is shown in Figure 5(b) that the control scheme can eliminate the initial error in about 3 seconds when the initial height error is -0.2 m and the aircraft can accurately track the deck motion within the last 15 seconds before landing. The height error of aircraft is shown in Figure 5(c), which demonstrates that within the last 15 seconds before landing, the maximum height tracking error of the aircraft is only 0.13 m, which is far less than the allowable height error of 1.5 m during the carrier landing process [24]. Meanwhile, the delay estimation results of PE-PSO are given in Figure 5(d). It can be seen that the PE-PSO algorithm can estimate the time delay in ACLS system effectively, and the final estimation error is less than 1 ms.The state variables and control variables are presented in Figure6. One can see that the control variables obtained using the AD-MPC carrier landing method are relatively stable, and there is no significant jitter. Additionally, the control variables are all constrained by the performance range of the corresponding actuator, which ensures the feasibility of the control law.Figure 6
Results of the state and control variables.


## 4.2. Control Performance Comparison

In order to verify the control performance of the AD-MPC carrier landing control algorithm, the traditional pole assignment method is selected for comparison. To get better comparison result, the initial value is increased to -0.35 m, and the initial conditionx0 is set as 10−2×0.43−0.17−0.30−0.501×9T.Firstly, setτ as 100 ms and the deck motion tracking results of the two algorithms are given in Figure 7. It is shown in Figure 7(a) that the aircraft can accomplish the glide path tracking under the control of pole assignment method when the deck motion is not introduced into the control system. However, there will be a significant phase lag on the aircraft height tracking result after the introduction of deck motion. Further, Figure 7(b) demonstrates that in the initial control stage, the response speed of pole assignment method is slightly faster than that of AD-MPC algorithm, but during the last 15 seconds before landing, the altitude tracking error range of pole assignment method is about ±0.55 m, while this error range of AD-MPC is ±0.13 m, which is only 1/4 of that of the pole assignment method.Figure 7
Carrier landing control effectiveness.
(a)
Deck motion tracking results(b)
Height errors of the aircraftFurthermore, when the time delay is set as 150 ms, the comparison results of the two algorithms are given in Figure8. As can be seen from Figure 8(a), due to the increase of time delay, the pole assignment method cannot effectively control the height of aircraft even if the deck motion is not considered. The height of aircraft will fluctuate continuously. After the introduction of deck motion, the phase lag still exists, and the tracking performance is worse than that of Figure 7(a) (τ=100ms). Further, it can be seen from Figure 8(b) that within the last 15 seconds before landing, the height error range of pole placement method is about ±0.6 m, while this error range of AD-MPC is ±0.14 m. When the time delay is increased to 200 ms, the pole assignment method will diverge, and it can be seen from Figure 5 that AD-MPC algorithm can still effectively control the aircraft.Figure 8
Carrier landing control effectiveness.
(a)
Deck motion tracking results(b)
Height errors of the aircraftThe above analysis shows that the AD-MPC carrier landing control method designed in this paper has stronger antidelay capability and higher control accuracy compared with the traditional pole assignment method.


## 4.3. Simulation Results for Time Delay Estimation Algorithm

The time delay estimation results of the two examples in Section4.2 are shown in Figures 9(a) and 9(b), respectively. It can be seen from Figure 9 that the estimation errors of the two cases are both less than 2 ms, which shows that the PE-PSO algorithm can estimate the time delay accurately.Figure 9
Time delay estimation results for the ACLS.
(a)
τ=100ms(b)
τ=150msTo further verify the effectiveness of the PE-PSO time delay estimation algorithm, control delays of 50, 100, 150, and 200 ms are introduced into the landing control system successively. And for comparison, the artificial fish school algorithm (AFS) and pigeon-inspired optimization (PIO) algorithm are selected to replace PSO in PE-PSO. Considering the four time delay scenarios, Monte Carlo simulations are carried out 500 times for every situation. The numbers of population sizes are all set as 20. And the statistical results are show in Figure10 and Table 2.Figure 10
Statistics of Monte Carlo simulation results.
(a)(b)(c)Table 2
Statistics of Monte Carlo simulation results.
Real time delay (ms)RMSE (ms)Mean fitness function (10-7)Mean CPU time (ms)AFSPIOPSOAFSPIOPSOAFSPIOPSO501.97291.97311.97718.40338.40338.4038551.83340.15147.431002.08682.08692.08758.55298.55288.5530565.67352.86148.081501.23721.24111.24378.96808.96768.9683541.63325.35137.822000.86660.86540.86229.23049.22989.2319614.15328.56141.37It can be seen from the statistical results that the PIO algorithm yields the smallest mean fitness function value (MFFV) among the three methods. The RMSE of AFS is the minimal except when time delay is 200 ms. The RMSE and the MFFV of the three methods are quite close to each other (the largest difference is less than 0.52 %), while the CPU time of AFS is almost 4 times of that of PSO for the unique trial mechanism of AFS. Meanwhile, the CPU time of PIO is about 2.3 times of that of PSO.In fact, there is only one variable to be optimized in time delay parameter estimation, and the optimization problem is not complicated. Comparing with the absolute accuracy of parameter estimation, we are more concerned about how to obtain the enough accurate delay estimate result in the shortest time. Shorter time consumption means faster updating frequency of the delay estimate result. And from the Monte Carlo simulations, it can be seen that the PE-PSO proposed in this paper meets this design requirement.


## 5. Conclusion

Automatic carrier landing control is a challenging problem based on a variety of disturbances and high control accuracy requirements. Additionally, the impact of time delays in an ACLS further increases the difficulty of the carrier landing control. In this study, the problem of automatic carrier landing control with time delay was examined. The main contributions of this paper can be summarized as follows.Based on the SP algorithm and PE-PSO time delay estimation algorithm, an AD-MPC carrier landing control algorithm is proposed in this paper. In the framework of MPC, the algorithm makes full use of the prediction information of deck motion and the steady-state component of carrier air wake to improve the control performance. And Padé approximation principle is also applied to transform the time delay control system into a nondelay one to avoid the control performance degradation caused by time delay.Simulation results demonstrated that the PE-PSO algorithm designed in this study can accurately estimate the pure time delay in an ACLS. Additionally, the AD-MPC algorithm can effectively overcome the adverse effects of time delays in the process of landing control. Under typical sea conditions, the steady-state height error is less than 0.14 m, which is far less than the allowable height error during the carrier landing process.



---

*Source: 8613498-2021-08-18.xml*