Flight path reconstruction and parameter estimation using output-error method

This work describes the application of the output-error method using the Levenberg-Marquardt optimization algorithm to the Flight Path Reconstruction (FPR) problem, which constitutes an important preliminary step towards the aircraft parameter identification. This method is also applied to obtain the aerodynamic and control derivatives of a regional jet aircraft from flight test data with measurement noise and bias. Experimental results are reported, employing a real jet aircraft, with flight test data acquired by smart probes, inertial sensors (gyrometers and accelerometers) and Global Positioning Systems (GPS) receivers.


Introduction
Modeling and simulation has become an integral part of the aeronautical industry design and evaluation processes. One of its major parts is system identification and parameter estimation applied to complex aerodynamic systems such as airplane. This work focuses the system identification process which is a general procedure to match the observed input-output response of a dynamic system by a proper choice of an input-output model and its physical parameters. From this point of view, the aircraft system identification or inverse modeling comprises proper choice of aerodynamic models, the development of parameter estimation techniques by minimization of the mismatch error between the predicted and the real aircraft response. In recent years, several works have been written describing system identification applied to aircraft [8][9][10]12]. Another focus of this work is the problem of Flight Path Reconstruction (FPR), which arises naturally when the main goal is an accurate identification of the aircraft parameters, because, in this case, the proper characterization of the flight sensors constitutes a fundamental preliminary step. For example, if the bias of a certain sensor is not adequately estimated, the accuracy of the ensuing parameter identification may be degraded.
The FPR is especially useful in the validation of the instruments applied in a prototype. The interpretation of the results can furnish important information with respect to sources of problems. Additionally, it decreases the uncertainties about the quality of data, which is one of the main causes of poor flight tests results. One of the first approaches for FPR may be found in [5]. In this work, the FPR problem is investigated by parametric identification of a nonlinear model, based on output-error method and Levenberg-Marquardt algorithm. The results are reported for a real jet aircraft, considering 28 parameters and 6 outputs, and comparing the calibration results obtained with those determined by traditional methods [3,13]. The identification method used in this work, based on the optimization algorithm of Levenberg-Marquardt is of the output-error type, which is susceptible to process noise. However, this approach is justified in the present case, since the noise is not large and because the methods presented in [5,11] could lead to incorrect results: the EKF could mask instrumentation errors, mainly those arising from inadequate compatibilization between the INS and GPS coordinates.
This work is structured as follows: in Section 2, the kinematics model for FPR and the lateral-directional model for parameter estimation are presented. In Section 3, the parametric estimation method is described with special attention to the Gauss-Newton and Levenberg-Marquardt algorithms. The experimental results obtained in the FPR problem and parameter estimation are then analyzed.

Kinematic model for FPR
The equations that constitute the kinematics model of an aircraft can be grouped in 3 sets of first-order differential equations, providing translational velocities, angular velocities and attitude angles. The axis systems are shown in Fig. 1. Using the standard body-fixed reference frame F B , the equations for the components u, v, w of true air speed V T along the body axes X B , Y B and Z B are: where p, q er denote the rates of rotation about the axes of F B ; θ, φ denote pitch and roll angle, respectively; m denotes aircraft mass and g denotes the local acceleration due to gravity. X, Y and Z represent the components of the total aerodynamic force, including the aerodynamic effects of propulsion systems.
For an aircraft with a geometrical plane of symmetry, the rotational dynamics are given by, where L, M, N denote the total aerodynamic moments, including any aerodynamic effects of the propulsion system; I x , I y and I z denote the moments of inertia and I xz the only non-zero product of inertia in F B (due to symmetry).
The orientation of F B with respect to the earth-fixed vertical reference frame F E is governed by the following equations for the Euler angles φ, θ, ψ.
To integrate the Eqs (1) and (3), it is necessary to determine X, Y, Z in Eq. (1). This is done by assuming that these accelerations are measured, giving. (4) in which a x , a y ea z denote the specific aerodynamic forces along the body axes X B , Y B and Z B , respectively.
By replacing Eq. (4) into Eq. (1) and dividing by m leads tȯ Once mass m (or any other physical properties) has been eliminated from Eqs (5) and (3), these equations can be integrated. More precisely, the solution of these equations can be obtained using the acceleration components (a x , a y , a z ) and the angular rates (p, q, r) as input variables, since they are measured by sensors installed on the aircraft, which are part of the inertial system. It is precisely the measure of these components that allow the realization of the FPR before the parameter identification of the aircraft is carried out.
Aiming to use GPS readings (geographical coordinates), we must characterize the position of the aircraft relative to the earth-fixed reference frame. To improve the quality of this signal, it was used the differential technique, namely the Differential Global Positioning System (DGPS), which is an extension of the GPS system that use terrestrial radio signals to transmit corrections of the position to the receptors. This position is obtained from Eqs (5) and (3), where L EB denotes an orthogonal matrix of reference frame transformation, determined by the roll, pitch and yaw angles. The three vectors that form this matrix are and denotes the components of a constant atmospheric wind vector W E along the axes of F E .
To summarize, the aircraft motion can be described by the nonlinear model Eqs (5), (3) and (6), which can be rewritten in the form with state and input vectors given by Basically, the observation models take the form of nonlinear algebraic relations between the observed variables and the state and input vector components. Aerodynamic forces and moments arise in the wind or aerodynamic axes systems. The wind axes system is obtained by successive rotations the body axes system over the sideslip angle and the angle of attack, see Fig. 2. In this work the models are derived for observations of true air speed V T , angle of attack α, side slip angle β and geographical position measurements.
True air speed V T can be derived from differential and absolute barometric and temperature transducers, resulting where K v is a scale factor and ΔV T the bias term. By definition, V T is the absolute value of the resultant of the air velocity components u, vew along the axes of F B , i.e., Also by definition, the angle of attack and the side slip angle are given by, respectively, These values differ from the measured angles, due to many effects, like velocities induced by aircraft rotational motion and modification of the air flow due to disturbs of the air near the aircraft, resulting the following measurement equations where K α , K α β and K β are scale factors, Δα and Δβ denote bias terms and the parameters x α , x β , y α and z β denote the position of the sensors. More details can be found in [11]. Finally, the geographical coordinates are obtained by DGPS. Therefore, the observation model takes the form where Δx E , Δy E and Δz E denote bias terms. Based in Eqs (12), (15) and (16), the observation vector is defined as From Eqs (5), (3), (6), (12), (15)-(17), and adding bias terms in the measurements of accelerations and angular velocities, the following dynamic model is obtained: state equations: control signals: output signals: where the last terms in Eq. (20) stand for sensor noise. In the dynamic model given by Eqs (18)-(20), the parameter vector to be estimated is formed by 28 components, namely where the last 9 terms in Eq. (21) denote initial conditions of the state vector in Eq. (11). Further details on the modeling procedure can be found in [7]. Therefore, from Eqs (18)-(21) we conclude that the FPR constitutes a parametric identification problem applied to a dynamic system of the form,

Dynamic model of lateral-directional movement of aircraft
The aircraft dynamic system is described by a stochastic nonlinear hybrid model in the form of Eq. (22). In this section the inverse problem formulation is applied to the lateral-directional movement of the aircraft, for which the linear state and output equations can be written as, Equation (23) has 14 unknown parameters that need to be estimated, i.e., As usually formulated in the aeronautical literature, the components of the parameter vector are the dimensional aerodynamic derivatives, which in turn can be written in term of non-dimensional coefficients by proper choice of flight parameters, all assumed known a priori. More details about this conversion can be found in [1].

Parametric estimation method
In this section, the parametric identification, in particular the parameter estimation applied to a linear causal model of an aircraft, in space state formulation according to Eq. (22). The output-error method is one of the most used estimation methods in aircraft identification and aerodynamic parameter estimation. It has several desirable statistical properties, including its application to nonlinear dynamical systems and the proper accounting of measurements noise [5].
The structure of the model is considered to be known, and the identification process consists in determining the parameter vector, which gives the best prediction of the output signal y(t), using some sort of optimization criteria. The attainment of an estimate through optimization of a cost function based on the prediction error of the plant requires, usually, the minimization of a nonlinear function. Thus, the Levenberg-Marquardt method is used here to estimate the parameters in model Eq. (23). Therefore, the cost function to be minimized involves the prediction error, whereŷ(k + 1,Θ) is the output prediction based on the actual estimate of the parameter vector.

Maximum likelihood estimation criteria
Consider a dynamic system, identifiable, with model structure M (Θ) defined and output y. Suppose that p(y|Θ) is the conditional probability Gaussian distribution of the random variable y with dimension m, mean f (Θ) and covariance F F T , with dimension m × m. p(y|Θ) is known as the likelihood functional, and in [2] the authors attribute its name due to the fact that it is a measure of the probability of occurrence of the observation y for a given parameter Θ. The Maximum Likelihood Estimate is defined as the value of Θ which maximizes this functional, in such a way that the best estimate of Θ, according to the MLE criteria is.
Thus, the likelihood functional is whose maximization is equivalent to the minimization of since, in the optimization process, J(Θ) is equivalent to − ln(p(y|Θ)), except for a constant term.

Minimization of the cost function by Levenberg-Marquardt
The identification algorithms based on the Gauss-Newton method is of second order. This method, although complex, is suitable for a quadratic cost function, and is expected to converge quickly. First, we approximate J(Θ) by a parabolic function (retaining only the 3 first Taylor series terms), The optimization condition is obtained when, which can be used to find the minima of the original cost function through the recursion, The complexity in the calculation of the Hessian matrix, in Eq. (31), is avoided through the Gauss-Newton method, which uses the approximation, where the terms involving the second derivative are discarded. The gradient of the estimated output, is called Sensibility Function. The Levenberg-Marquardt algorithm is an extension of the Gauss-Newton. The idea is to modify ∇ 2 Θ J(Θ) to ∇ 2 Θ J(Θ) + λI , and the inversion of the matrix is not done in a explicit way, but solving by Singular Value Decomposition (SVD) the following expression The inclusion of λI in Eq. (33) solves the problem of an ill conditioned approximated Hessian. The Levenberg-Marquardt algorithm can be interpreted in the following way: for small values of λ it behaves like the Gauss-Newton algorithm, while for high values of λ it behaves like the gradient method.

Flight path reconstruction
A flight test was performed and data was gathered with sampling time T = 0.09 s. Based on the input signals and in the measured variables according to the output vector, the identification algorithm was used to determine the parameter vector containing the 28 parameters indicated in Eq. (21). The identification algorithm was executed many times, aiming to investigate the influence of the design parameters. The influence of the integration method used to solve the state equation was also investigated, concluding that the Euler method is not suitable, but the 4th order Runge-Kutta method produces adequate results.
After executing the identification, we must evaluate its performance. The first quality measure is the mean square prediction error. The plots of these errors in the present case indicate small values, and these plots are omitted,  except for the air speed, which is shown in Fig. 3. The relative vertical scale indicates a total range of 6 m/s. Hence, the difference between the measured and the predicted values of air speed is small.
Next, it is considered the main performance measure here, the comparison of the estimated values with those obtained by the aerospace industry via traditional procedures. Two of these variables are considered here: the angle of attack and the side slip angle. In Fig. 4 is presented the relative values of angle of attack. Three plots compose Fig. 4: the values measured by the aircraft sensors, the values calibrated by the aerospace industry through traditional  procedures, and the values calibrated by the method proposed in this paper. The vertical scale omits the absolute values of the angles, but presents the total variation. Based on Fig. 4, we conclude that the procedure proposed here for the flight path reconstruction presents results compatible with that obtained by traditional techniques employed by the aeronautical industry.

Matching of flight test data for lateral-directional movement
The aerodynamic derivatives associated with the lateral-directional model, as shown in Eq. (23), were estimated by matching the real flight test data with the model predicted simulation. A Dutch-roll maneuver of a regional transport aircraft was used to investigate the effectiveness of the discussed output-error method (the Levenberg-Marquardt), applied to estimate the aerodynamic parameter vector.
The aircraft input signals are the aileron and rudder deflections, and the output signals are five attitude parameters: sideslip angle, roll rate, yaw rate, bank angle, and lateral acceleration. The experimental input signals are shown in  Table 1 shows the final values of the non-dimensional aerodynamic derivatives obtained by the Levenberg-Marquardt algorithm. We use the values achieved by the Nelder-Mead method to initialize the Levenberg-Marquardt algorithm, since this method is more computationally demanding and a good initial estimate     can speed up its convergence. A maximum likelihood cost was used, in which case the weighting factor was the estimated covariance matrix associated to the prediction errors.
Since the flight data employed to generate Table 1 was obtained experimentally and no wind tunnel tests are available, an indirect measure of performance is used, based on the prediction error. So, the main focus of the present inverse aerodynamic modeling is to check that this local minimization procedure can provide good matching to the experimental flight data and stable input-output modeling for the aircraft. This prediction capability, as obtained by the output-error method, can be accessed from the model validation results shown in Figs 8 to 12, where the estimation error are small for most of the output variables.

Conclusions
Based on the small prediction error and good agreement between the calibrated values, we can conclude that the proposed procedure for FPR, based on parametric identification via optimization using the Levenberg-Marquardt method, exhibited satisfactory performance. Thus, the proposed procedure can be added to the repertory of FPR techniques and constitute a relevant alternative for practical applications. This work also presented the estimation of aircraft linear aerodynamic derivatives, which presented good convergence properties and good matching to the experimental flight data. The results obtained with the Levenberg-Marquardt algorithm demonstrate the feasibility of the method.