An Accurate Method for Real-Time Aircraft Dynamics Simulation Based on Predictor-Corrector Scheme

Real-time aircraft dynamics simulation requires very high accuracy and stability in the numerical integration process. Nonetheless, traditional multistep numerical methods cannot effectively meet the new requirements. Therefore, a novel real-time multistep method based on Predict-Evaluate-Correct scheme of three-step fourth-order method (RTPEC-34) is proposed and developed in this research to address the gap. In addition to the development of a highly accurate algorithm based on predictor-corrector, the contribution of this work also includes the analysis of truncation error for real-time problems. Moreover, the parameters for the RTPEC-34 method are optimized using intelligent optimization algorithms. The application and comparison of the optimization algorithms also lead to general guidelines for their applications in the development of improved multistep methods. Last but not least, theoretical analysis is also conducted on the stability of the proposed RTPEC-34method, which is corroborated in simulation experiments and thus provides general guidelines for the evaluation of real-time numerical methods. The RTPEC-34 method is compared with other multistep algorithms using both numerical experiments and a real engineering example. As shown in the comparison, it achieves improved performance in terms of accuracy and stability and it is also a viable and efficient algorithm for real-time aircraft dynamics simulation.


Introduction
Aircraft dynamics simulation is a complex nonlinear process whereby engine, aerodynamic, and atmospheric models are solved simultaneously [1].Aircraft variables such as aircraft velocity, orientation angles, aerodynamic angles, and angular rates are assembled by ordinary differential equations (ODEs) with initial value problems (IVP).The real-time simulation of aircraft dynamics requires very accurate and stable numerical methods for solving the ODEs.
A numerical solution of an ODE is a table of approximated values of the variables for a set of discrete time points.Multistep methods, as a class of the most widely used numerical solution, use information at more than one previous point to estimate the solution at the next point.Linear multistep methods have the form shown in (1) where the parameters   and   are determined using polynomial interpolation [2] If  0 = 0, the method is explicit; and if  0 ̸ = 0, then the method is implicit.Since multistep methods require several values of variables and derivatives from previous calculations history, they cannot self-start before the values have been obtained.Using a single-step method, which does not require past history, is a common strategy to generate solution values at enough points for starting using a multistep method of a desired order.The Euler method and the Runge-Kutta (RK) method (or variants of them) are popular single-step methods to assist multistep methods in providing initial values.We know that implicit single-step methods are generally more 2 Mathematical Problems in Engineering accurate and stable than the explicit ones.Similarly, implicit multistep methods are usually more accurate and stable than the explicit ones but a nonlinear equation needs to be solved to obtain the derivative values  +1 .On the other hand, explicit linear multistep methods can conveniently make a guess of initial values and have simple and direct formulas.These advantages help meet the requirements of real-time input and calculation.As a result, the explicit and implicit methods are often used together as a predictor-corrector pair and the Predict-Evaluate-Correct (PEC) method has become a very useful scheme for solving ODEs with IVP.When computational complexity and the error constants relative to order are concerned [3], the most popular fourth-order Adams-Bashforth-Moulton (ABM) [2] of predictor-corrector method is widely used.This method, however, does not satisfy the needs of real-time simulation [4].In addition, although implicit methods have a much greater region of stability than explicit methods, they are still not necessarily unconditionally stable.Explicit methods cannot be A-stable due to the famous Dahlquist barrier, which means it is less useful for solving stiff ODEs [5] and a PEC scheme is actually explicit in a complicated way.On the contrary, a properly designed implicit multistep method is suitable for dealing with stiff problems.Since the accuracy restriction depends on the slowly varying components of the solution and the stability restriction depends on the rapidly varying ones, the Gear method [6] and the Rosenbrock method [7] aim to solve stiff ODEs with long intervals of stability regions, at the cost of accuracy.
According to the different requirements and constraints of a simulation, many researchers improved the multistep methods from different aspects.Huang [8] considered twostage hybrid methods of seventh order.Gottlieb et al. [9] considered a class of two-step and two-stage methods and Xu and Zhao [3] proposed an estimation of the longest stability interval for three-order four-step methods.In addition, Bulatov and Berghe [10] discussed two-step fourthorder methods of the second order and Bresten et al. [11] focused on strong stability preserving methods.Seong et al. [12] considered the fifth-order multistep method using constant step size and Wang et al. [13] discussed the variablestep interaction algorithm for multidisciplinary collaborative simulation.Although these methods attempted to improve the method by varying the order and number of integration steps, very few of them focused on optimizing parameters for a fixed order.When stability and accuracy of numerical solutions are considered at the same time, a predictorcorrector pair of order four is a very good choice based on which improved methods can be developed.
This research, with a particular focus on the aircraft dynamics ODE problems, aims to improve the accuracy of the real-time predictor-corrector methods and analyze the stability regions of the proposed methods.Both simulation experiments of numerical examples and an engineering example are used to evaluate the performance of the proposed method by comparing its results with those obtained from the classic RK method and the methods for stiff ODEs.It has been shown in the evaluation that these proposed methods are effective.The rest of this paper is organized as follows.
Section 2 details the development of a predictor-corrector method for highly accurate real-time aircraft dynamics simulation as well as the optimization of its key parameters.In Section 3, a theoretical analysis is conducted to evaluate the stability of the RTPEC-34 method.In Section 4, the RTPEC-34 method is evaluated by comparing it with other classic multistep algorithms in several simulation experiments.The main conclusions of this research are drawn in Section 5.

Optimization of the Accuracy of Multistep Numerical Methods Using Predictor-Corrector Scheme
Aircraft simulation involves a computation process that is subject to real-time constraints as it is a real-time system.This means the responses or results of the simulation process have a deadline that must be met, regardless of system load, in order for the system considered to be correct [14,15].The predictor-corrector multistep formula is rearranged as the real-time form discussed in [16].In particular, the predictor equation is obtained as follows: Evaluate (( +1 ),   +1 ) at this predicted value of   +1 to obtain  +1 at  +1 .The corrector equation is then obtained as In order to improve the accuracy of the real-time predictor-corrector multistep methods, the step of the predictor equation needs to be reduced.Consider the fourth-order three-step algorithm as follows: where parameters  0 ,  1 ,  2 ,  ] ,  0 ,  1 , and  2 are yet to be determined.Numerical convergence can be achieved for formula (3) when the following constraint equations are established [17]: Additionally, based on the method of undetermined coefficients [18], formula (3) can be further constrained by five equations as follows: Thus, (4a) and (4b) include seven undetermined equations which contain eight unknown variables.The multistep numerical methods for solving an ODE suffer from two distinct sources of error.Compared to truncation error, round-off error generally plays only a minor role [19,20].Therefore, truncation error is the main focus of this research and the estimation of the p-order truncation error can be done using the following formula: The local truncation error coefficient is then represented as For (3), the estimation of truncation error coefficient can be represented as Based on the above discussion, the problem has now become how to figure out the minimum value of formula (7) with the constraints described in (4a) and (4b).
There are many optimization methods to find out the minimum value with nonlinear constraints, such as feasible direction method (FDM) [21], penalty function method, and quadratic programming (QP) [22].However, these methods are not accurate enough to solve multivalue equality.The genetic algorithm (GA) [23] has the potential to fix the problem of falling into local optimum compared with the classic methods mentioned above.While the equality nonlinear constraints in the above model are strict, the infeasible solution will become a large proportion in the populations as the genetic algorithm is used separately [24].
Hence, some other optimization methods should be utilized to assist the GA [25].The penalty function has the advantages of considering the points out of the feasible regions when solving nonlinear constraints problems [26].In this research, the GA with a punishment strategy is chosen to find solutions outside the feasible regions.
A penalty function is generally constructed using the addition form below: where  is the chromosome, () is the objective function,   () is th penalty element, and   is the changeable penalty coefficient of the   ().The element ∑  =1     () in formula ( 8) is for multiplying constraint violations, which will be used to alter feasible region unless the original function is unsolvable.The complete expression of penalty function then can be established as follows: To simplify the constraint, use V to represent  0 ,  1 , and  2 : Thus, the objective function (7) can be rearranged as Then, the following three penalty elements are used for the GA optimization process, which are derived from (4a), (4b), and (10): The setting for the GA running is as follows: binary code is used; the number of individuals in the population is 1000; the crossover operator is 0.9; and mutation operator is 0.08.A simulation is done using MATLAB software, and then the parameters can be obtained for the minimum truncation error of the RTPEC-34 method.
Compared with the classic optimization methods, it can be found that different parameters are obtained as shown in Table 1.
Apply the parameters to (3); a family of the RTPEC-34 method with the highest accuracy can be obtained.The RTPEC-34 method using FDM optimization is shown as follows: The RTPEC-34 method using QP optimization is shown as follows: The RTPEC-34 method using GA optimization with a punishment strategy is shown as follows:

Stability Analysis
A widely used approach to determining stability region is to apply the method to the linear ODE / =  with initial condition (0) =  0 , and its analytical solution is given as () =  0   .Applying the three-step fourth-order ABM algorithms (12a), (12b), and (12c) to this equation using a fixed step-size ℎ, we have Substituting  +V into  +1 , we can obtain the following stability polynomial: where  = ℎ.
Obtaining RTPEC-34 stability polynomial with the parameters shown in Table 1, we can then use root-locus method to obtain the diagram for the stability regions of this method as shown Figure 1.
The curves  1 ,  2 , and  3 are the root locus of RTPEC-34 (13a), (13b), and (13c) stability polynomial, respectively.Since no multistep methods of greater than second order are unconditionally stable, an explicit linear multistep method Im u Re u cannot be A-stable if it is higher than order two due to the famous Dahlquist barrier theorem.Even if it is implicit, it is still the case.Thus the fourth-order method can be stable only in the finite stability region.From Figure 1, it can be concluded that the highest accuracy RTPEC-34 method (13c) using GA with a punishment strategy has the largest stability region among the three methods.

Evaluation and Discussion
4.1.Numerical Experiments.In this section, the numerical results obtained by using the highest accuracy RTPEC-34 method are presented in Example 1.In addition, the method proposed is compared with the Runge-Kutta method of fourth-order (RK4) method to demonstrate its good performance in supporting real-time simulation in Example 2. Furthermore, the Rosenbrock method for stiff ODEs is used to evaluate the accuracy of the proposed method.Example 1.The ODE with IVP ( 16) is used to compare the proposed methods (13a), (13b), and (13c) using different parameters.Consider The analytical solution of this ODE is () = √ 1 + 2.Applying (13a), (13b), and (13c) to the ODE ( 16), the numerical calculation errors in different simulation settings with different step-size ℎ are shown in Table 2.According to Figure 2, the accuracy of (13c) is better than that of (13a) and (13b) under the condition of ℎ = 0.005, ℎ = 0.01, and ℎ = 0.02.The RTPEC-34 methods using GA optimization with a punishment strategy (13c) are much more stable than the other two cases.Thus, we choose (13c) as the method with the highest accuracy after the optimization and comparison.It is used to demonstrate the efficiency of the proposed methods in other simulation experiments.
Example 2. To compare the proposed method (13c) with the RK4 method, the real-time simulation procedure of an aircraft propulsion system [27] is considered in which the motion ODE with IVP has the following form: Applying (13c) and the RK4 to system (17), numerical solutions can be obtained, as shown in Figure 3.It is shown in the curves in Figure 3 that the solutions obtained using the RTPEC-34 method come earlier in time (for about one) than the RK4 method, which is a significant phenomenon in the real-time simulation.This means the RTPEC-34 method (13c) is much more adaptive than RK4 when the aircraft

RTPEC-34 (13c) RK4
Figure 3: The numerical solutions obtained using the RTPEC-34 method and the RK4 method.dynamics simulation requires the numerical solutions in real time.

Aircraft Dynamics Experiment.
The aircraft dynamics experiment is focused on the simulation of the motion model, which is a system of 12 scalar order differential equations [1].To avoid the singularities of derivatives of roll angular rate ṗ and yaw angular rate ṙ when pitch angle  passes through ±(/2), quaternion methods are used for the aircraft orientation representation [28].We have got an aircraft system representation consisting of 13 scalar first order differential equations as follows: where To demonstrate the effectiveness and efficiency of the proposed method (13c), it has been applied to a real aircraft dynamics simulation problem.A number of comparisons with other classic linear multistep methods such as the fourth-order Adams-Bashforth method (AB4) have been made in the simulation which has high real-time requirements The nonlinear dynamics model of the F-16 aircraft [29] is used in the experiment as the simulation in such an engineering problem requires a detailed model and practical data.For the atmospheric data relative to the coefficient of formulas (20a), (20b), (20c), and (20d) an approximation is used based on the guidelines by the International Standard Atmosphere (ISA) [30].In the experiment, the flight status of the F-16 aircraft with the altitude at 1000 m and the velocity at 200 m/s and the initial values of variables are calculated by using the ISA atmospheric model, aerodynamics model, and engine model.The simulation is performed for the F-16 aircraft dynamical model using different numerical methods, the RTPEC-34 (13c) and the AB4 method.The configuration for the simulation is shown in Table 4.
As demonstrated in these simulation experiments, the aircraft variables curves in Figure 5 can illustrate the real state of a whole fight, where the roll rate, pitch rate, and yaw rate [, , ] tend to be stable in Figures 5(a  , which means that the convergence speed of RTPEC-34 (13c) is faster than that of AB4.With the number of iteration increased, the accuracy of RTPEC-34 will be higher than AB4.The example of the F-16 aircraft simulation demonstrates that the proposed method is adaptive in the field of dynamic ODEs.The comparison with the classic multistep methods of the same order shows the good accuracy and speed of convergence of the RTPEC-34 (13c) and confirms that the proposed real-time fourth-order three-step predictorcorrector method achieves the highest accuracy.

Conclusions
In this paper, a novel real-time multistep method based on Predict-Evaluate-Correct scheme of three-step fourth-order (RTPEC-34) method is proposed and developed for the realtime simulation of aircraft dynamics, which is improved by obtaining the predictor-corrector parameters using optimization algorithms.The method of using GA optimization with a punishment strategy is more stable and accurate than others both theoretically and practically.This work involves theoretical analysis of the truncation error of the predictorcorrector for real-time problems as well as of the stability of the proposed RTPEC-34 method.The analysis work can be used by other applications that involve improvement of accuracy and stability.Both numerical examples and engineering examples are used in simulation experiments to evaluate the performance of the RTPEC-34 method by comparing it with the RK4 and the Rosenbrock methods.It is shown in the evaluation that the proposed method achieves improved performance in terms of accuracy, stability, and support of real-time simulation.Moreover, the successful application of the proposed method in the F-16 aircraft experiment has shown that the proposed method is adaptive to solve the multivariable ODEs popular in aircraft dynamics simulation.

Figure 1 :
Figure 1: The region of absolute stability of the RTPEC-34 method.

Figure 4 :
Figure 4: Comparison of the numerical solutions obtained by the RTPEC-34 and Rosenbrock methods.
)∼5(c).This is also the case for [, ] in Figure5(d) through to Figure 5(f) and the coordinates [  ,   ,   ] of the fight position in the

Figure 5 :
Figure 5: Simulation results of different numerical methods for the F-16 aircraft model.

Table 1 :
Parameter values obtained using different optimization methods.

Table 3 :
Comparison of numerical solutions obtained with different step sizes.