Recent Advances on Modeling, Control, and Optimization for Complex Engineering Systems
View this Special IssueResearch Article  Open Access
Jiaxin Zhao, Hongwei Wang, Heming Zhang, "An Accurate Method for RealTime Aircraft Dynamics Simulation Based on PredictorCorrector Scheme", Mathematical Problems in Engineering, vol. 2015, Article ID 193179, 10 pages, 2015. https://doi.org/10.1155/2015/193179
An Accurate Method for RealTime Aircraft Dynamics Simulation Based on PredictorCorrector Scheme
Abstract
Realtime 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 realtime multistep method based on PredictEvaluateCorrect scheme of threestep fourthorder method (RTPEC34) is proposed and developed in this research to address the gap. In addition to the development of a highly accurate algorithm based on predictorcorrector, the contribution of this work also includes the analysis of truncation error for realtime problems. Moreover, the parameters for the RTPEC34 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 RTPEC34 method, which is corroborated in simulation experiments and thus provides general guidelines for the evaluation of realtime numerical methods. The RTPEC34 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 realtime aircraft dynamics simulation.
1. 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 realtime 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 , the method is explicit; and if , then the method is implicit. Since multistep methods require several values of variables and derivatives from previous calculations history, they cannot selfstart before the values have been obtained. Using a singlestep 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 RungeKutta (RK) method (or variants of them) are popular singlestep methods to assist multistep methods in providing initial values. We know that implicit singlestep methods are generally more 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 . 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 realtime input and calculation. As a result, the explicit and implicit methods are often used together as a predictorcorrector pair and the PredictEvaluateCorrect (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 fourthorder AdamsBashforthMoulton (ABM) [2] of predictorcorrector method is widely used. This method, however, does not satisfy the needs of realtime 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 Astable 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 twostep and twostage methods and Xu and Zhao [3] proposed an estimation of the longest stability interval for threeorder fourstep methods. In addition, Bulatov and Berghe [10] discussed twostep fourthorder methods of the second order and Bresten et al. [11] focused on strong stability preserving methods. Seong et al. [12] considered the fifthorder 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 realtime predictorcorrector 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 predictorcorrector method for highly accurate realtime 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 RTPEC34 method. In Section 4, the RTPEC34 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.
2. Optimization of the Accuracy of Multistep Numerical Methods Using PredictorCorrector Scheme
Aircraft simulation involves a computation process that is subject to realtime constraints as it is a realtime 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 predictorcorrector multistep formula is rearranged as the realtime form discussed in [16]. In particular, the predictor equation is obtained as follows:
Evaluate at this predicted value of to obtain at . The corrector equation is then obtained as
In order to improve the accuracy of the realtime predictorcorrector multistep methods, the step of the predictor equation needs to be reduced. Consider the fourthorder threestep algorithm as follows: where parameters , , , , , , and 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, roundoff error generally plays only a minor role [19, 20]. Therefore, truncation error is the main focus of this research and the estimation of the porder truncation error can be done using the following formula:The local truncation error coefficient is then represented asFor (3), the estimation of truncation error coefficient can be represented asBased 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 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 to represent , , and :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 RTPEC34 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 RTPEC34 method with the highest accuracy can be obtained. The RTPEC34 method using FDM optimization is shown as follows:
The RTPEC34 method using QP optimization is shown as follows:
The RTPEC34 method using GA optimization with a punishment strategy is shown as follows:
3. Stability Analysis
A widely used approach to determining stability region is to apply the method to the linear ODE with initial condition , and its analytical solution is given as . Applying the threestep fourthorder ABM algorithms (12a), (12b), and (12c) to this equation using a fixed stepsize , we have
Substituting into , we can obtain the following stability polynomial:where .
Obtaining RTPEC34 stability polynomial with the parameters shown in Table 1, we can then use rootlocus method to obtain the diagram for the stability regions of this method as shown Figure 1.
The curves , , and are the root locus of RTPEC34 (13a), (13b), and (13c) stability polynomial, respectively. Since no multistep methods of greater than second order are unconditionally stable, an explicit linear multistep method cannot be Astable 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 fourthorder method can be stable only in the finite stability region. From Figure 1, it can be concluded that the highest accuracy RTPEC34 method (13c) using GA with a punishment strategy has the largest stability region among the three methods.
4. Evaluation and Discussion
4.1. Numerical Experiments
In this section, the numerical results obtained by using the highest accuracy RTPEC34 method are presented in Example 1. In addition, the method proposed is compared with the RungeKutta method of fourthorder (RK4) method to demonstrate its good performance in supporting realtime 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 . Applying (13a), (13b), and (13c) to the ODE (16), the numerical calculation errors in different simulation settings with different stepsize 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 , , and . The RTPEC34 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.
(a)
(b)
Example 2. To compare the proposed method (13c) with the RK4 method, the realtime 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 RTPEC34 method come earlier in time (for about one) than the RK4 method, which is a significant phenomenon in the realtime simulation. This means the RTPEC34 method (13c) is much more adaptive than RK4 when the aircraft dynamics simulation requires the numerical solutions in real time.
Example 3. Due to the quality of being Astable, the Rosenbrock method is a halfimplicit RK method used for solving stiff ODEs. The Rosenbrock method of order three has the following form:where , , , , , and . However, it is difficult to estimate the truncation error directly. Hence, the following ODE (19) is used to compare the RTPEC34 method (13c) with the Rosenbrock method in terms of accuracy:
The analytical solution of ODE (18) is .
Table 3 shows the exact solutions and numerical solutions using the RTPEC34 method (13c) and the Rosenbrock method in the interval . It is noted that the accuracy of RTPEC34 method (13c) is much better than that of the Rosenbrock method at in Figure 4. It is indicated in this phenomenon that ODEs can be solved more accurately by the proposed method unless they are extremely stiff.

(a)
(b)
4.2. 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 , 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:
whereThe 12 aircraft variables, [], [], [], and [], represent velocity, sideslip angle, attack angle, bodyaxis rolling rate, pitch rate and yaw rate, xposition, yposition, and zposition with respect to earth, rolling angle, pitch angle, and yaw angle, respectively. For the sake of calculation simplicity, 7 parameters are introduced, where [] are velocity in axis, axis, and axis direction and [] are quaternion components which do not have any physical meaning.
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 fourthorder AdamsBashforth method (AB4) have been made in the simulation which has high realtime requirements
The nonlinear dynamics model of the F16 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 F16 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 F16 aircraft dynamical model using different numerical methods, the RTPEC34 (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)~5(c). This is also the case for [] in Figure 5(d) through to Figure 5(f) and the coordinates [] of the fight position in the whole air route in Figure 5(g) through to Figure 5(i). Both of the numerical methods are effective in the simulation of the F16 aircraft at the altitude of 1000 m; however the curve of RTPEC34 demonstrates smaller fluctuations than that of AB4 in Figures 5(a) and 5(b), which means that the convergence speed of RTPEC34 (13c) is faster than that of AB4. With the number of iteration increased, the accuracy of RTPEC34 will be higher than AB4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
The example of the F16 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 RTPEC34 (13c) and confirms that the proposed realtime fourthorder threestep predictorcorrector method achieves the highest accuracy.
5. Conclusions
In this paper, a novel realtime multistep method based on PredictEvaluateCorrect scheme of threestep fourthorder (RTPEC34) method is proposed and developed for the realtime simulation of aircraft dynamics, which is improved by obtaining the predictorcorrector 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 realtime problems as well as of the stability of the proposed RTPEC34 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 RTPEC34 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 realtime simulation. Moreover, the successful application of the proposed method in the F16 aircraft experiment has shown that the proposed method is adaptive to solve the multivariable ODEs popular in aircraft dynamics simulation.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research is supported by the National Natural Science Foundation of China (no. 61374163), the National High Technology Research and Development Program (863 Program) of China (no. 2013AA041302), and the Aviation Science Foundation of China (no. 20124060179).
References
 B. L. Stevens and F. L. Lewis, Aircraft Control and Simulation, vol. 2, John Wiley & Sons, New York, NY, USA, 2003.
 M. T. Heath, Scientific Computing, McGrawHill, New York, NY, USA, 1997.
 Y. Xu and J. J. Zhao, “Estimation of longest stability interval for a kind of explicit linear multistep methods,” Discrete Dynamics in Nature and Society, vol. 2010, Article ID 912691, 18 pages, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 Z. A. Majid, N. Z. Mokhtar, and M. Suleiman, “Direct twopoint block onestep method for solving general secondorder ordinary differential equations,” Mathematical Problems in Engineering, vol. 2012, Article ID 184253, 16 pages, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 G. G. Dahlquist, “A special stability problem for linear multistep methods,” BIT Numerical Mathematics, vol. 3, pp. 27–43, 1963. View at: Publisher Site  Google Scholar  MathSciNet
 L. F. Shampine and C. W. Gear, “A user's view of solving stiff ordinary differential equations,” SIAM Review, vol. 21, no. 1, pp. 1–17, 1979. View at: Publisher Site  Google Scholar  MathSciNet
 H. H. Rosenbrock, “Some general implicit processes for the numerical solution of differential equations,” The Computer Journal, vol. 5, no. 4, pp. 329–330, 1963. View at: Google Scholar  MathSciNet
 C. Huang, “Strong stability preserving hybrid methods,” Applied Numerical Mathematics, vol. 59, no. 5, pp. 891–904, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 S. Gottlieb, C.W. Shu, and E. Tadmor, “Strong stabilitypreserving highorder time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001. View at: Publisher Site  Google Scholar  MathSciNet
 M. V. Bulatov and G. V. Berghe, “Twostep fourth order methods for linear ODEs of the second order,” Numerical Algorithms, vol. 51, no. 4, pp. 449–460, 2009. View at: Publisher Site  Google Scholar  MathSciNet
 C. Bresten, S. Gottlieb, Z. Grant, D. Higgs, D. I. Ketcheson, and Németh, “Strong stability preserving multistep RungeKutta methods,” http://arxiv.org/abs/1307.8058. View at: Google Scholar
 H. Y. Seong, Z. A. Majid, and F. Ismail, “Solving secondorder delay differential equations by direct AdamsMoulton method,” Mathematical Problems in Engineering, vol. 2013, Article ID 261240, 7 pages, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 H. Wang, H. Mao, and H. Zhang, “A variablestep interaction algorithm for multidisciplinary collaborative simulation,” Integrated ComputerAided Engineering, vol. 21, no. 3, pp. 263–279, 2014. View at: Google Scholar
 F. Bergero and E. Kofman, “PowerDEVS: a tool for hybrid system modeling and realtime simulation,” Simulation, vol. 87, no. 12, pp. 113–132, 2011. View at: Publisher Site  Google Scholar
 X. Wang, J. Zhang, and M. Scalia, “Parallel motion simulation of largescale realtime crowd in a hierarchical environmental model,” Mathematical Problems in Engineering, vol. 2012, Article ID 918497, 15 pages, 2012. View at: Publisher Site  Google Scholar
 O. B. Widlund, “A note on unconditionally stable linear multistep methods,” BIT Numerical Mathematics, vol. 7, pp. 65–70, 1967. View at: Publisher Site  Google Scholar  MathSciNet
 K. E. Atkinson, An Introduction to Numerical Analysis, John Wiley & Sons, New York, NY, USA, 2008. View at: MathSciNet
 Q. Y. Li, Z. Guan, and F. S. Bai, Numerical Calculation Principle, Tsinghua University Press, Beijing, China, 2000.
 T. E. Hull, W. H. Enright, B. M. Fellen, and A. E. Sedgwick, “Comparing numerical methods for ordinary differential equations,” SIAM Journal on Numerical Analysis, vol. 9, no. 4, pp. 603–637, 1972. View at: Publisher Site  Google Scholar
 H. M. Zhang, S. Liang, S. J. Song, and H. W. Wang, “Truncation error calculation based on Richardson extrapolation for variablestep collaborative simulation,” Science China Information Sciences, vol. 54, no. 6, pp. 1238–1250, 2011. View at: Publisher Site  Google Scholar
 G. Zoutendijk, Methods of Feasible Directions; A Study in Linear and Nonlinear Programming, Elsevier, New York, NY, USA, 1960.
 M. Frank and P. Wolfe, “An algorithm for quadratic programming,” Naval Research Logistics Quarterly, vol. 3, no. 1, pp. 95–110, 1956. View at: Publisher Site  Google Scholar  MathSciNet
 C. R. Houck, J. Joines, and M. Kay, “GA genetic algorithm for function optimization: a Matlab implementation,” NCSUIE TR, 1995. View at: Google Scholar
 O. Kramer and H.P. Schwefel, “On three new approaches to handle constraints within evolution strategies,” Natural Computing, vol. 5, no. 4, pp. 363–385, 2006. View at: Publisher Site  Google Scholar  MathSciNet
 F. Zong, H. Lin, B. Yu, and X. Pan, “Daily commute time prediction based on genetic algorithm,” Mathematical Problems in Engineering, vol. 2012, Article ID 321574, 20 pages, 2012. View at: Publisher Site  Google Scholar
 C. Coello and C. Artemio, “Constrainthandling techniques used with evolutionary algorithms,” in Proceedings of the 14th International Conference on Genetic and Evolutionary Computation Conference Companion, ACM, 2012. View at: Google Scholar
 J. S. Litt, D. L. Simon, S. Garg et al., “A survey of intelligent control and health management technologies for aircraft propulsion systems,” Journal of Aerospace Computing, Information and Communication, vol. 1, no. 12, pp. 543–563, 2004. View at: Google Scholar
 L. Stevens, Aircraft Control and Simulation, John Wiley & Sons, 1992.
 R. S. Russel, “Nonlinear f16 simulation using Simulink and matlab,” Tech. Rep., University of Minnesota, 2003. View at: Google Scholar
 J. A. Mulder, W. H. J. J. van Staveren, and J. C. van der Vaart, “Flight dynamics,” Tech. Rep., Delft University of Technology, 2000. View at: Google Scholar
Copyright
Copyright © 2015 Jiaxin Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.