Abstract
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-34 method, 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.
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 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 , 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 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 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 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 two-stage 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 fourth-order 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 variable-step 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 predictor-corrector 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.
2. 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 at this predicted value of to obtain at . 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 , , , , , , 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, 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 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 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:
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 three-step fourth-order ABM algorithms (12a), (12b), and (12c) to this equation using a fixed step-size , we have
Substituting into , 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 , , and 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 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.
4. 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 . 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 , , and . 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.
(a)
(b)
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 dynamics simulation requires the numerical solutions in real time.
Example 3. Due to the quality of being A-stable, the Rosenbrock method is a half-implicit 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 RTPEC-34 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 RTPEC-34 method (13c) and the Rosenbrock method in the interval . It is noted that the accuracy of RTPEC-34 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, body-axis rolling rate, pitch rate and yaw rate, x-position, y-position, and z-position 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 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)~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 F-16 aircraft at the altitude of 1000 m; however the curve of RTPEC-34 demonstrates smaller fluctuations than that of AB4 in Figures 5(a) and 5(b), 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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
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 predictor-corrector method achieves the highest accuracy.
5. 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 real-time 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 predictor-corrector 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.
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).