3.1 Sensitivity analysis

Sensitivity analysis is normally used to determine the influence of the model’s parameters on its dynamic response. In the context of validation, the sensitivity analysis can be defined as the systematic investigation on the simulation responses to the model’s input^{(Reference Jategaonkar10,Reference Sargent26)} . In the environment of flight simulation, the behavior of flight response depends on the root location of the characteristic equations, which translates into the value of characteristic parameters such as frequency and damping ratio. This section presents the sensitivity analysis of the flight dynamic in longitudinal motion. The state equation and characteristic equation for the short-period approximation can be defined in the following formulas^{(Reference Nelson27)}:

(1) \begin{equation}\left[ \begin{array}{c}\Delta \dot{\alpha } \\[5pt] \Delta \dot{q} \end{array}\right]=\left[ \begin{array}{c@{\quad}c}\dfrac{Z_{\alpha }}{u_0} & 1 \\[12pt] M_{\alpha }+M_{\dot{\alpha }}\frac{Z_{\alpha }}{u_0} & M_q+M_{\dot{\alpha }} \end{array}\right]\left[ \begin{array}{c}\Delta \alpha \\[5pt] \Delta q \end{array}\right]\end{equation}

(2) \begin{equation}{\lambda }^2-\left(M_q+M_{\dot{\alpha }}+\frac{Z_{\alpha }}{u_0}\right)\lambda +M_q\frac{Z_{\alpha }}{u_0}-M_{\alpha }=0\end{equation}

The short-period roots can be obtained from the characteristic equation:

(3) \begin{equation}{\lambda }_{sp}=\frac{M_q+M_{\dot{\alpha }}+\frac{Z_{\alpha }}{u_0}}{2{{\omega }_n}_{sp}}\pm \frac{\sqrt{{\left(M_q+M_{\dot{\alpha }}+\frac{Z_{\alpha }}{u_0}\right)}^2-4\left(M_q\frac{Z_{\alpha }}{u_0}-M_{\alpha }\right)}}{2}\end{equation}

The frequency and damping ratio are calculated as:

(4) \begin{equation}{\omega }_{n_{sp}}=\sqrt{M_q\frac{Z_{\alpha }}{u_0}-M_{\alpha }}\end{equation}

(5) \begin{equation}{\zeta }_{sp}=-\frac{M_q+M_{\dot{\alpha }}+\frac{Z_{\alpha }}{u_0}}{2{{\omega }_n}_{sp}}=-\frac{M_q+M_{\dot{\alpha }}+\frac{Z_{\alpha }}{u_0}}{2\sqrt{M_q\frac{Z_{\alpha }}{u_0}-M_{\alpha }}}\end{equation}

Given the preceding modeling of longitudinal motion, the short period characteristic parameters are functions of a very large set of parameters. Some of these parameters are considered fixed quantities in the early design stage such as aircraft geometry dimensions, whereas some other parameters are known to change nominally with the aircraft flight envelope such as center of gravity (CG) position. Therefore, the short period sensitivity analysis focuses on the effects of parameters listed in Table 1. For example, Equations (4) and (5) denote that the primary aerodynamic coefficients involved in ${c_{{m_\alpha }}} = {C_{{L_\alpha }}}\left( {{x_{CG}} - {x_{AC}}} \right)$ have direct implications on ${M_\alpha }$, and then consequently on ${\omega _{{n_{sp}}}}$ and ${\zeta _{sp}}$. Therefore, the longitudinal locations of the center of gravity $\left( {{x_{CG}}} \right)$ and aerodynamic center $\left( {{x_{AC}}} \right)$ have major implications on the longitudinal dynamic stability. In a similar manner, this analysis is applied for all the influenced parameters.

Table 1

The influence of relevant parameters on affected derivatives

3.2 Parameter tuning technique

The effect of the influenced parameters $\left( {{x_{CG}},{x_{AC}},{I_{yy}}} \right)$ on the characteristic parameters $\left( {{\zeta _{sp}},{\omega _{{n_{sp}}}}} \right)$ can be plotted based on the motion characteristic equations, as shown in Fig. 2. The most sensitive parameters are selected based on its effect on a specific motion. Assume that the motion characteristic can be explained as a function of the influenced parameters:

(6) \begin{equation}{\sigma }_i=f\left(p_j\right)\ \ \ \ \ with\ i=\left[\zeta ,{\omega }_n\right]\ \ \ and\ j=[x_{CG},\ x_{AC},I_{yy}]\end{equation}

Where ${\sigma _i}$ is the motion characteristic parameter such as damping ratio $\zeta $ or natural frequency ${\omega _n}$, and ${p_j}$ is the tuning parameters such as ${x_{CG}},\;{x_{AC}}$ and ${I_{yy}}$. The sensitive parameters can be tuned by applying the equation:

(7) \begin{equation}p^{tuned}_j=p^{init}_j+\dot{f}(p^{init}_j)\left[{\sigma }^{test}_i-{\sigma }^{sim}_i\right]\end{equation}

Where the superscript “tuned” presents the tuned parameters, “init” presents the initial parameter, “test” and “sim” is parameters obtained from the flight test and the simulation data, respectively.

Figure 2.

Longitudinal parameters sensitivity analysis result. a) Impact of $\textit{I}_{\textit{yy}}$, b) Impact of $\textit{X}_{\textit{CG}}$ and c) Impact of $\textit{X}_{\textit{AC}}$.

4.1 Maximum Likelihood method

As introduced earlier, the system identification process includes not only the model postulating and determining parameters, but performing suitable experiments, and gathering the system inputs and responses. This research utilises the classical time-domain Maximum Likelihood (ML) cost function and the Gauss-Newton (GN) computational algorithm to estimate the longitudinal dynamic derivatives from the light aircraft flight test data. The flow diagram presented in Fig. 3 shows the block schematic of the ML estimation method.

Figure 3.

The Maximum Likelihood (ML) parameter estimation concept^{(Reference Jategaonkar10)}.

Referring to Fig. 3, the dynamic system under investigation is simply assumed to be modeled by a set of dynamic mathematical equations containing unknown parameters. To determine the values of the unknown parameters, the system is excited by a suitable control input, and then the actual system response is measured. The difference between the measured response and estimated response, which is called the response error, comes to the ML cost function module. The GN computational algorithm is applied to minimise the cost function, which is a function of the difference between the measured and computed time histories^{(Reference Jategaonkar10,Reference Morelli and Klein13)} . The probability that the aircraft response time histories attain values near the observed values is defined for each estimate of the unknown parameters^{(Reference Iliff28)}. By maximising this probability, the unknown parameters are estimated. These unknown parameters are then used to update the mathematical model of the aircraft. The updated mathematical model is used to provide a new computed response and therefore a new response error. The updating continues iteratively until the mathematical model response matches the actual flight test response according to the convergence criteria. To clarify details of the ML method, the stepwise algorithm of the estimation program is presented in Fig. 4. First, the state and observation equations, Equations (9) and (10), of the dynamic system are defined. The control input $u\left( t \right)\;$and flight response $z\left( t \right)$ is loaded from flight test data, the initial value of unknown parameter $\xi $ is specified. Then, the state equations are integrated into the ML cost function $J$, the iteration process is applied to minimise the cost function by updating the unknown parameters. The ML estimation is achieved when the convergence criterion is satisfied.

Figure 4.

The program algorithm of the ML estimation method.

Since the aircraft is a continuous-time dynamic system, assume that the measurements are made at discrete time interval, the system equations are defined in the following^{(Reference Maine and Iliff29)}:

(8) \begin{equation}x\left(t_0\right)=x_0\left(\xi \right)\end{equation}

(9) \begin{equation}\dot{x}\left(t\right)=f\left[x\left(t\right),u\left(t\right),t,\xi \right]\end{equation}

(10) \begin{equation}z\left(t_i\right)=g\left[x\left(t_i\right),u\left(t_i\right),t_i,\xi \right]+{\eta }_i\end{equation}

Where $x\left( t \right)$ is the state, $u\left( t \right)$ is the input and $z\left( t \right)$ is the response of the dynamic system. Assume that the measurement noise ${\eta _i}$ is a sequence of independent Gaussian random vector with zero mean. The forms of the $f$ and $g$ function are known, where the aerodynamic format of the f and g function is described detailed in the case study, Equations (16) and (17) in Section 5.4. The parameter estimation task is to estimate the value of the unknown parameter $\xi $. The ML estimates of $\xi $ is the value that minimise the cost function:

(11) \begin{equation}J\left(\xi \right)=\frac{1}{2}\sum^N_{i=1}{{\left[z\left(t_i\right)-{\tilde{z}}_{\xi }\left(t_i\right)\right]}^*}{\left(GG^*\right)}^{-1}\left[z\left(t_i\right)-{\tilde{z}}_{\xi }\left(t_i\right)\right]+\frac{1}{2}Nln\left|(GG^*)\right|\end{equation}

Where $ {\tilde{z}}_{\xi}(t_i)$ is the predicted response estimate of $z$ at ${t_i}$ for a given value of the unknown-parameter vector $\xi $. The ${\left( {G{G^*}} \right)^{ - 1}}$ is the measurement noise covariance matrix. The last term of the cost function is omitted by assuming there is no prior distribution to be considered. The Gauss-Newton optimisation algorithm is applied in order to minimise the cost function $J\left( \xi \right)$. Considering $k$ as the iteration number, the $\left( {k + 1} \right)$ estimate of $\xi $ is obtained by the $k$ estimate as:

(12) \begin{equation}{\widetilde{\xi }}_{k+1\ }={\widetilde{\xi }}_k+{\left[{\nabla }^2_{\xi }J\left({\widetilde{\xi }}_k\right)\right]}^{-1}\left[{{\nabla }^{{\rm *}}}^2_{\xi }J({\widetilde{\xi }}_k)\right]\end{equation}

If ${\left( {G{G^*}} \right)^{ - 1}}$ is assumed fixed, the first and second gradients are defined as:

(13) \begin{equation}{\nabla }_{\xi }J\left(\xi \right)=-\sum^N_{i=1}{{\left[z\left(t_i\right)-{\tilde{z}}_{\xi }\left(t_i\right)\right]}^*}{\left(GG^*\right)}^{-1}\left[{\nabla }_{\xi }{\tilde{z}}_{\xi }\left(t_i\right)\right]\end{equation}

(14) \begin{equation}{\nabla }^{{\rm 2}}_{\xi }J\left(\xi \right)=\sum^N_{i=1}{{\left[{\nabla }_{\xi }{\tilde{z}}_{\xi }\left(t_i\right)\right]}^*}{\left(GG^*\right)}^{-1}\left[{\nabla }_{\xi }{\tilde{z}}_{\xi }\left(t_i\right)\right]\end{equation}

The minimisation algorithm is used to find the derivatives that minimise the cost function. Each iteration of this algorithm provides a new estimate of the unknown derivatives on the basis of the response error. These new estimates are then used to update values of the derivatives of the mathematical model, providing a new estimated response and new response error. The updating of the mathematical model continues until a convergence criterion is satisfied.

4.2 Aerodynamic database correction

This part proposes a methodology to combine the existing aerodynamic DB with the estimated derivatives from flight test data. Initially, the baseline aerodynamic coefficients are in different format according to the analysis method. For example, the longitudinal aerodynamic coefficients obtained by W/T are mostly in 2D- and 3D-table format, whereas the dynamic coefficients calculated by AVL are in linear format (Table 3). However, the coefficients estimated from flight test data are all single values (or in linear format). Therefore, to update the existing aerodynamic DB, a correction technique is addressed to modify the linear region of the data. The linear correction formular is presented in Equation (15), where a 2D-table format coefficient is recalculated by two linear format coefficients. In terms of the dynamic coefficients such as ${C_{{L_q}}}$ and ${C_{{m_q}}}$, the initial coefficients computed by AVL are in linear format already. Therefore, these coefficients are simply updated by replacing the flight test estimations.

(15) \begin{equation}C^{mod}\left(x,y\right)=C^{init}\left(x,y\right)+\left[C^{est}_x\left(y\right)-\frac{\partial C^{init}\left(x,y\right)}{\partial x}\right]\times x+\left[C^{est}_y\left(x\right)-\frac{\partial C^{init}\left(x,y\right)}{\partial y}\right]\times y\end{equation}

Where the x, y represent the variables, the superscript “mod” presents the modified coefficient, the “init” presents the initial coefficient, and the “est” presents the flight test estimation.

5.1 Aircraft model specifications

In this paper, the KLA-100 light aircraft model is generated and simulated in the JSBSim Flight Dynamics Model (FDM). The KLA-100 is a two-seater light aircraft with low-wing configuration developed by the Light Aircraft Development Center of Konkuk University, Korea^{(6,Reference Tyan, Nguyen, Kim and Lee30,Reference Van Nguyen, Tyan, Lee and Kim31)} . The simulation model database consists of four modules: configuration DB, propulsion DB, mass and balance DB and aerodynamic DB.

The picture and three-view drawing with general specifications of the KLA-100 aircraft are shown in Fig. 5 and Table 2.

Figure 5.

The KLA-100 light aircraft – 3D view geometry^{(Reference Tyan, Nguyen, Kim and Lee30)}.

Table 2

Light aircraft KLA-100 specifications^{(Reference Tyan, Nguyen, Kim and Lee30)}

Table 3

The force components in aerodynamic database

5.2 Aerodynamic database construction

The aerodynamic DB is a comprehensive DB that consists of three force and three moment components. Table 3 shows an example of three force components together with their data format and analysis method.

The basic aerodynamic coefficients are classified into static and dynamic terms; the static components are obtained from CFD analysis or W/T static stability experiment data^{(Reference Van Nguyen, Tyan, Lee and Kim31,Reference Van Nguyen, Lee, Tyan, Lee and Kim32)} ; while the dynamic components are extracted from the dynamic stability analysis such as Athena Vortex Lattice (AVL) or empirical calculation.

5.3 Flight test data processing

The full-scale prototype of the light aircraft KLA-100 has been constructed for the flight test validation and certification compliance purposes. The flight test operation team determines the flight conditions and the control inputs based on the manoeuver description in the flight simulator evaluation handbook⁽⁷⁾. The flight test data is recorded by means of the on-board data acquisition system. The raw data includes all the parameters related to flight status and motion variables such as orientations, velocities, accelerations, and control surface deflections as well as throttle setting⁽³³⁾.

For the system identification and simulation validation purposes, the whole flight test data needs to be processed and classified into different manoeuver database packages. In this research, a flight test manoeuver classification process is introduced to identify flight motions based on their typical characteristics. Figure 6 presents the detailed steps of the flight test manoeuver classification process. First, a pre-processing module is applied to eliminate noise and extract the whole flight test sortie into segments which include several single manoeuvers^{(Reference Jategaonkar10)}. Here, a moving average filtering technique is utilised to smooth the raw data by replacing each data point with the average of the neighboring data points defined within the span. Then, using the extraction rules specified for each manoeuver such as load factor or control input, the full data is extracted into several single-manoeuver segments. The extracted manoeuvers are classified automatically into the specific type of motion using signal pattern similarity searching technique^{(Reference LV Thang, Maxim, Seung-Hyeog and Jae-Woo34)}.

Figure 6.

The flight test manoeuver classification process.

The signal pattern similarity searching technique recognises a time-varying segment, which is closely matching with the sample signal by using the distance matric. The basic concept of this method is to sweep the sample signal across the data and then looking for the lowest difference between the sample signal and the data locally at each location. The Euclidian distance has been used to measure the difference. In addition, the sequences of the signals need to be warped in a nonlinear fashion to match each other. The dynamic time warping algorithm^{(Reference Li35)} has been used to find an optimal alignment between two given time-dependence sequences under certain restrictions. The sample signal is generated regarding their typical characteristic description. For instance, short period motion is a fast-responding longitudinal mode that is initiated from horizontal level flight by applying an elevator or other longitudinal control input, while phugoid mode can be excited through either an elevator pulse or thrust variation^{(Reference Ward, Strganac and Niewoehner36)}. In this case, the elevator single doublet has been applied for short-period longitudinal motion. Figure 7 shows an example for searching a single doublet signal in one control input segment of the flight test data, where the dot line presents the segment signal that closely matches with the sample signal.

Figure 7.

Single doublet control input searching result.

5.4 Longitudinal derivatives estimation

In this subsection, an application of the ML estimation method is demonstrated in a basic longitudinal motion, particularly the short-period mode. First, the state equations for short-period mode are defined as follows:

(16) \begin{equation}\dot{\alpha }=Z_0+Z_{\alpha }\alpha +\left(Z_q+1\right)q+Z_{\delta e}{\delta }_e\end{equation}

(17) \begin{equation}\dot{q}=M_0+M_{\alpha }\alpha +M_qq+M_{\delta e}{\delta }_e\end{equation}

The observation equations:

(18) \begin{equation}{\alpha }_m=\alpha \end{equation}

(19) \begin{equation}q_m=q,\end{equation}

where subscript m denotes the measured variable from the flight test. The unknown parameter vector ${\rm{\Theta }}$ consisting of the dimensional derivatives and the bias parameters is given by:

(20) \begin{equation}{\Theta }^T=[Z_0\ Z_{\alpha }\ Z_q\ Z_{\delta e}\ M_0\ M_{\alpha }\ M_q\ M_{\delta e}]\end{equation}

The postulated model is defined in Table 4.

Table 4

Definition of variable parameters

The unknown parameter vector ${\rm{\Theta }}$ consisting of the dimensional derivatives and the bias parameters is:${\rm{\;}}{{\rm{\Theta }}^T} = \left[ {{Z_0}\;{Z_\alpha }\;{Z_q}\;{Z_{\delta e}}\;{M_0}\;{M_\alpha }\;{M_q}\;{M_{\delta e}}} \right]$. The MATLAB program for the ML method⁽¹⁰⁾ has been used to estimate the longitudinal derivatives from flight test data of the light aircraft KLA-100. Figure 8 shows the convergence of the estimated longitudinal derivatives for the KLA-100 aircraft. As discussed in Section 4.1, the ML estimation is achieved when the convergence criteria are satisfied. In this case, the iteration process is stopped when the convergence criterion in terms of the relative change in the cost function meets the tolerance at ${10^{ - 5}}$ deg/s after 20 iterations.

Figure 8.

Convergence of longitudinal derivatives for KLA-100 aircraft.

From the estimated longitudinal derivatives, the related coefficients can be calculated using the Equations (21)–(26). The calculated aerodynamic coefficients are shown in Table 5.

(21) \begin{equation}C_{L_{\alpha }}=-\left(\frac{mZ_{\alpha }}{QS}+C_{D_0}\right)\end{equation}

(22) \begin{equation}\hspace*{-19pt}C_{m_{\alpha }}=\frac{u_0I_{yy}M_{\alpha }}{QS\overline{c}}\hspace*{19pt} \end{equation}

(23) \begin{equation}C_{L_q}=-\frac{2\mu_0Z_q}{QS\overline{c}}\end{equation}

(24) \begin{equation}C_{m_q}=\frac{2u_0I_{yy}M_q}{QS{\overline{c}}^2}\end{equation}

(25) \begin{equation}C_{L_{\delta e}}=-\frac{2\mu_0Z_{\delta e}}{QS\overline{c}}\end{equation}

(26) \begin{equation}C_{m_{\delta e}}=\frac{I_{yy}M_{\delta e}}{S\overline{c}}\end{equation}

Table 5

The aerodynamic derivatives and relative coefficients estimation

5.5 Aerodynamic database correction

Table 6 shows the value of the estimated parameters and the baseline aerodynamic DB format. The modification of aerodynamic coefficients can be obtained by applying the developed linear correction technique in Section 4.2.

Table 6

The estimated parameters and the baseline aerodynamic DB format

As observed from Fig. 9(a) and (b), the surface presents the baseline aerodynamic DB, while the scatter points present the updated results using database correction method. The result figures depict that the gradient of the lift coefficient is slightly changed due to modification according to the flight test result.

Figure 9.

Aerodynamic coefficient modification results. a) Modification of $\textit{C}_{\textit{L}_{\alpha\beta}}$ and b) Modification of $\textit{C}_{\textit{L}_{\alpha\delta e}}$.

6.1 Simulation results

To evaluate the efficiency of the proposed enhancement approach, the flight simulation response and the flight test data are compared side-by-side regarding the evaluation criteria. The simulation results are conducted by putting the same control input and flight condition into the 6 DOF JSBSim FDM. A control input segment that including the series of four elevator doublet deflection has been tested. Figure 10 shows the comparison between the actual flight test (continuous lines), the simulation responses of the initial model (dotted lines) and the flight responses of the enhanced model (dashed line) in time histories. Three evaluation criteria have been plotted for comparison of the short-period mode are the pitch angle $\theta $, the pitch rate $p$ and the normal acceleration ${a_z}.$

Figure 10.

Simulation longitudinal response comparison. a) Test 1, b) Test 2, c) Test 3 and d) Test 4.

It can be noticed that the initial simulation model depicts a major margin of error in the flight response, which violates the regulation limits, compared to the flight test measured. After applying the model tuning methodology in stage 1 and the corrected aerodynamic DB in stage 2, the accuracy of the simulation model significantly improves and closely matches the actual result from the flight test. To provide a more detailed evaluation of the enhanced simulation model, a quantitative and statistical analysis is addressed in the following validation section.

6.2 Validation of flight simulation

The simulation model is validated by the tolerances defined in the simulator regulations. The simulator regulations are the documents concerning for aircraft simulator qualification, published by the FAA⁽²⁴⁾ and the RAeS⁽⁷⁾. The validation procedure is applied for short-period dynamic of longitudinal motion regarding the flying qualities section. Referring to Fig. 10, the shaded area represents the tolerance band according to the flight measurement, which is $ \pm {1.5^\circ }$ for pitch angle, $ \pm {2^\circ }/{\rm{s}}$ for pitch rate and $ \pm 0.10{\rm{g}}$ for normal acceleration. The flight simulator is expected to be validated if the simulation response does not violate the upper and lower limitations. It can be inferred that the response of the initial model does violate the tolerance band whereas the response of the enhanced model lies within the limitations.

Table 7 presents a comparison between the initial and the enhanced simulation model with the flight test data according to the regulation criteria. The comparison indicates that the accuracy of the simulation response is significantly improved in the enhanced model. As an illustration, the average maximum pitch angle error is reduced from ${1.59^\circ }$ for the initial model, to ${1.08^\circ }\;$for the enhanced model. In the similar manner, the average maximum errors also show a considerable reduction in terms of pitch rate and normal acceleration, from 3.60$\;$ to ${1.36^\circ }/{\rm{s}}$ and from $0.16$ to $0.08{\rm{g}}$, respectively. The satisfaction of these regulation criteria shows that the initial model does not satisfies the regulation, whereas the enhanced model satisfies all the requirements. In addition to the regulation criteria, the major characteristic parameters for the short-period motion are also taken into comparison. The results indicate that the errors are greatly reduced from 10.3% to 6.2% for damping ratio and from 14.3% to 8.0% for undamped natural frequency. To summarise, the enhanced simulation model shows significant improvement from the initial model to the enhanced model in stage 1 and stage 2, in both the regulation criteria and the motion characteristics.

Table 7

The evaluation criteria comparison between the initial and the enhanced simulation model