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.

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 Section 2, 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]: where , and represent the velocity, angle of attack, pitch angle, pitch angle velocity, and altitude deviation relative to the nominal state, respectively. Additionally, , and represent deviations in terms of the horizontal tail deflection, leading-edge flap deflection, rudder toe-in deflection, and engine throttle control angle, respectively. Finally, represents the deviation of the attack angle caused by vertical wind disturbances. The nominal state for carrier landing is , , and . 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:

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]: where and represent random initial phases and and 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: where , , and represent the air wake in the , , and directions, respectively, in the inertial frame. The subscript (, 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, and 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 Table 1 [11].

In Table 1, 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.

Figure 1 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 , which is obtained according to the state of the control system at -.

After the position and speed of the carrier aircraft at time are measured and filtered by the shipboard radar, the guidance system on the carrier calculates the corresponding information at time . 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 . After a series of coding, transmission, and decoding operations, the actuators receive the control law information and begin to work at , resulting in a pure time delay for the entire closed loop of . 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 in the longitudinal landing model given in Equation (1) must be replaced with the thrust response variable .

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:

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: where and are intermediate variables. Additionally, to introduce the rate of change information for , and into the system model explicitly, it is necessary to augment Equation (1) in the following manner: where

The 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 as .

Considering the time delay in the control system, the longitudinal small disturbance mathematical model for carrier landing can be defined as follows: 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 is transformed using the Laplace transform. Then, by applying the first-order Taylor approximation, we have resulting in

The intermediate variable is selected to satisfy

Then, the original system can be rewritten as

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 where is the ideal flight state vector of the carrier-based aircraft, is the actual flight state vector, and is the actual control input. and are diagonal coefficient matrices, where is a semipositive definite matrix and 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.

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

A crucial problem for MPC tracking control is selecting a standard tracking trajectory , 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 Figure 3. In Figure 3, represents the position of the aircraft at the current time, and represents its height deviation from the ideal position . is the expected position of the aircraft after , and is the ideal trajectory obtained based on deck motion prediction. is the reference trajectory obtained by moving in a parallel direction and is used to generate the guidance trajectory . The points on satisfy the relationship . Therefore, the guidance trajectory 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 and are introduced into the guidance trajectory, and the corresponding height difference between and is adjusted as follows:

When the values of and are both 1.0, the guidance trajectory is represented by in Figure 3.

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: where , and by introducing a nonnegative relaxation vector , the inequality constraint in Equation (16) can be rewritten as follows:

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 where is the Hamilton function:

According to the parametric variational principle, the following equations should be satisfied when the objective function reaches its minimum value:

According to the inequality constraints and the Karush-Kuhn-Tucker condition, we have

Equation (16) defines an optimal control problem with a fixed terminal time. Therefore, when the terminal state is free, we have according to the transversality condition.

The variables , , , and are approximated using the -degree Legendre-Gauss-Lobatto (LGL) pseudospectral method in the th subinterval as follows: where is the Lagrangian interpolation polynomial corresponding to the th LGL node within the th subinterval [19].

The stagnation point condition for the second-generation function [20, 21] is applied to the th subinterval, and we have

According to the constraint in Equation (21) and the complementarity condition, we obtain the following relationship for each subinterval :

By assembling Equations (23) and (24) in every interval according to the boundary conditions, the two-point boundary value problem (TPBVP) in can be rewritten as follows: where , , 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 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 and can be derived according to Equation (25)(a). Finally, the control variable can be obtained by substituting and into Equation (20).

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 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 in Equation (19), where .

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 : where

Here, and represent the values of () and (), respectively. When or is a noninteger, and are obtained via linear interpolation, where and represent the estimators of 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 in cannot be measured, the observable state variable in the system should be predicted according to the following prediction model to obtain its estimated value at : where is the observation matrix, and the actual observable value can be expressed as where . represents all historical datasets of observable variables before , and represents the datasets of at and before . The prediction error at is denoted as

When comparing the prediction model to the actual dynamic model, it can be seen that 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, 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 observation values, the prediction error sequence can be expressed as , and the estimated value of the delay variable is calculated as

Theorem 1. When the sequence 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 is minimized.

This can be proven as follows: where

Assuming that the statistical characteristics of each stochastic process do not change over time (i.e., stationary stochastic process), when , each item in the formula above converges to its mean value with a probability of one. where

Here, and 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, is constant. Because 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, and are independent of each other, implying

Similarly,

When ,

When is taken as the minimum value, the following equation must be satisfied [23]:

According to the definition of , if , then holds with a probability of one. Therefore, the above formula is equivalent to and Theorem 1 is proven.

To obtain that satisfies Equation (32), we transfer the problem of extremum determining to an optimation problem as follows: where 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 Figure 4, where a dashed line or dotted box indicates that the corresponding step is only performed once.

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:

The initial state of the aircraft , which means the initial error of is -0.2 m. The control parameter matrices are and . The control parameters are and . 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.

Further, the deck motion tracking result is given in Figure 5(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 Figure 6. 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.

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 condition is set as .

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.

Furthermore, when the time delay is set as 150 ms, the comparison results of the two algorithms are given in Figure 8. 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) (). 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.

The 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 Section 4.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.

To 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 Figure 10 and Table 2.

It 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.

Appendix

A. Reference [21]

The second type of generating function on the time interval is expressed as

By taking the variation of Equation (A.1), and considering Equation (20), we have

Then, the time interval is divided into subintervals (). And each subinterval is transformed into the time interval by the linear transformation

The generating function in the th subinterval is expressed as

The variables , , , and are approximated using the -degree Legendre-Gauss-Lobatto (LGL) pseudospectral method in the th subinterval as follows: where

In (A.6), is the -degree Legendre polynomial of the interval [-1, 1]. are the LGL points, which are the roots of equation . And ; .

Because and are regard as “independent variables” in the th subinterval, the remaining components of the state and costate in the th subinterval are donated as

Define the first-order differentiation matrix as

Then, Equation (A.2) can be written as where

Based on Equation (A.2), by applying the variational principle to Equation (A.9) for each subinterval, one has where

Detailed expressions of , , , , , , , and are listed as follows:

Furthermore, Equation (A.11) can be rewritten into a compact form as where

Define

Equation (A.11) within all subintervals can be rewritten into a more compact form as where

In Equation (A.18), is symmetric since Equation (A.17) is obtained by variational principle.

Consider the th subinterval and impose Equation (21) at all LGL nodes, we have where

All relationships of Equation (A.21) within the th subinterval can be rewritten as where

Equation (A.23) within all subintervals can be rewritten into a more compact form as where

Similarly, the inequality in Equation (21) within all subintervals can be rewritten as follows:

The boundary conditions of free terminal state are applied to Equation (A.17). And Equation (A.17) can be modified as where the detailed expressions of , , , and can be found in reference [19].

By deducing Equation (A.28), can be written as where and .

Then, and can be expressed as the linear functions of as where where refers to or and the symbol refers to the submatrix which contains rows of the matrix .

Substituting Equation (A.30) into Equation (A.25), one has where

Then, we get a standard LCP as follows:

The LCP can be solved by the classical Lemke’s algorithm. Then, and can be obtained by Equation (A.30). Finally, the control variables can be obtained based on the first equation of Equation (20).

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.