Dynamic Modeling from Flight Data with Unknown Time Skews Eugene A. Morelli * NASA Langley Research Center, Hampton, Virginia, 23681 A method for estimating dynamic model parameters from flight data with unknown time skews is described and demonstrated. The method combines data reconstruction, nonlinear optimization, and equation-error parameter estimation in the frequency domain to accurately estimate both dynamic model parameters and the relative time skews in the data. Data from a nonlinear F-16 aircraft simulation with realistic noise, instrumentation errors, and arbitrary time skews were used to demonstrate the approach. The approach was further evaluated using flight data from a subscale jet transport aircraft, where the measured data were known to have relative time skews. Comparison of modeling results obtained from time-skewed and time-synchronized data showed that the method accurately estimates both dynamic model parameters and relative time skew parameters from flight data with unknown time skews. Nomenclature a ,a ,a = body-axis translational accelerometer measurements, g x y z b = wing span, ft c = mean aerodynamic chord, ft C ,C ,C = body-axis nondimensional aerodynamic force coefficients X Y Z C ,C ,C = body-axis nondimensional aerodynamic moment coefficients l m n Ε[] = expectation operator [] = Fourier transform operator g = gravitational acceleration, ft/s2 Im[] = imaginary part I ,I ,I ,I = mass moments of inertia, slug-ft2 x y z xz j = imaginary number = −1 m = aircraft mass, slug M = number of frequencies M = body-axis pitching moment from engine thrust, ft-lbf T p,q,r = body-axis roll, pitch, and yaw rates, rad/s or deg/s P = ambient static pressure, lbf/ft2 a q = dynamic pressure, lbf/ft2 Re[] = real part S = wing reference area, ft2 T = maneuver length, s T ,T = body-axis engine thrust components, lbf x z u,v,w = body-axis air-relative velocity components, ft/s V = airspeed, ft/s * Senior Research Engineer, Dynamic Systems and Control Branch, MS 308, Associate Fellow 1 American Institute of Aeronautics and Astronautics α = angle of attack, rad or deg β = sideslip angle, rad or deg δ,δ,δ ,δ = stabilator, elevator, aileron, and rudder deflections, rad or deg s e a r φ,θ,ψ = Euler roll, pitch, and yaw angles, rad or deg θ = parameter vector Σ = covariance matrix superscripts T transpose ˆ estimate  time derivative  Fourier transform –1 matrix inverse † complex conjugate transpose subscripts cg = center of gravity I = inboard L = left r = reconstructed R = right o = reference value or base term O = outboard I. Introduction EV ERY commonly-used method for dynamic modeling based on measured data makes an implicit assumption that all measured signals are sampled at the same time1,2. In practice, this is never exactly true, due to issues such as multiplexing in the analog-to-digital conversion electronics, and different sampling rates for various signals, for example. In high-quality aircraft instrumentation systems designed to collect data for dynamic modeling, such issues are well known, and steps are taken to minimize the differences in the sampling times for the measured data. The usual approaches involve accurately recording the time when each sample was taken (time tagging), and/or linear interpolation to approximate simultaneous data sampling. Unfortunately, not every dynamic modeling task is undertaken with data that have been carefully time tagged or interpolated. Some examples are data from an aircraft accident, and data from a low-budget instrumentation system, such as might be found on munitions, unmanned air vehicles, or remotely-piloted vehicles. Past approaches to solving this problem have included estimating the time skews separately, as well as estimating model parameters and time skews simultaneously in a combined nonlinear parameter estimation formulation3-5. In some practical situations, the effort necessary to time synchronize the sampled data is prohibitively expensive, and sometimes not even possible. This paper reports the results of an investigation focused on finding a dynamic modeling method that can be used in cases where the sampled signals may have unknown relative time skews. This can simply mean that the assumption of simultaneous data sampling has been violated, and the exact sampling times of the measured signals are not known. But it can also mean that there are practical or instrumentation issues that produce a time skew in the signal before it is sampled. Common examples are lags in air flow angle data due to sensor position on the aircraft or pressure tubing lengths, and control surface data time skews from low sampling rates or the use of control surface command data rather than control surface position data. For conventional dynamic modeling methods, the presence of relative time skews in the measured data has known detrimental effects on dynamic modeling accuracy6. 2 American Institute of Aeronautics and Astronautics A real-world example of the importance of this issue is the first flight of the X-43A hypersonic research vehicle, or Hyper-X, which was aborted during the launch phase7. Investigation of the event focused (among other things) on the dynamic model of the booster and Hyper-X stack (see Fig. 1) in the transonic region. Uncertainty in the sample times of the measured signals used in the investigation led to an extensive effort to accurately determine what those sample times were. It was very important to the mishap investigation to determine the sample times accurately, so that the dynamic modeling based on measured flight data would produce accurate results, and the root cause of the mishap could be identified. A dynamic modeling method that could be applied regardless of the actual relative time skews in the data would have saved a large amount of Figure 1. Hyper-X stack being carried time and expense in the mishap investigation. to altitude on the NASA B-52B Credit: NASA Armstrong Flight Research Center The purpose of this paper is to explain and demonstrate a method for accurately identifying dynamic models based on flight data with unknown time skews. The next section describes the general idea and gives some background information. Following this, the method is applied to data from an F-16 nonlinear simulation, using realistic noise, instrumentation errors, and known time skews in the data. Finally, flight data from a subscale transport aircraft were used to demonstrate that the method produces accurate modeling results for data with unknown time skews, and can also serve as a diagnostic tool to identify data that have significant relative time skews. II. Method In the time domain, it is difficult to accurately determine time skews in sampled data, because of the quantization related to the sampling interval and the fact that implementing any time skew value that is not an integer multiple of the sampling interval involves interpolation of the data. This causes convergence problems in time-domain parameter estimation algorithms, as well as resolution limitations for the estimated values of the time skews. However, if the analysis is done in the frequency domain, the time skews are continuous-valued parameters, and these problems do not exist. Consequently, the approach described here uses frequency-domain methods. The problems that arise from time-skewed data are related to relative time skews, or time skews of sampled time series relative to one another. If all of the time series are time skewed in exactly the same way, then for dynamic modeling purposes, there is no problem, and the modeling can proceed normally. When the data have relative time skews, these time skews can be falsely interpreted as phase lag due to the system dynamics, for example, and this causes errors in the model parameter estimates. A good approach to the problem is to first define a reference for the time skews in the data. Body-axis translational accelerations and angular rates typically come from a single Inertial Measurement Unit (IMU), so that the time skew between these measured quantities is negligible, and the data from the IMU provide a good reference for the time skews. Time skews for other sampled time series can be defined relative to this reference. With this timing reference established, there are typically two categories of data whose time skew must be determined relative to the reference: 1) Air flow data: airspeed V , sideslip angle β, and angle of attack α 2) Control surface deflection data: elevator δ , aileron δ , and rudder δ (and possibly other controls) e a r For example, aircraft longitudinal short-period dynamics can be described by the following linear equations1: α(t−τ )=Z α(t−τ )+(1+Z )q(t)+Z δ (t−τ ) (1a) α α α q δe e δe q(t) =M α(t−τ )+M q(t)+M δ (t−τ ) (1b) α α q δe e δe 3 American Institute of Aeronautics and Astronautics a (t)=Vo Z α(t−τ )+Z q(t)+Z δ (t−τ ) (1c) z g  α α q δe e δe  where the dimensional stability and control derivatives Z , Z , Z ,M ,M ,M , and the relative time skews α q δe α q δe τ ,τ , are unknown parameters to be estimated from the data. The aircraft response and control variables in α δ e Eqs. (1) are perturbation quantities, as usual for linearized dynamic equations1. Equations (1) can be readily transformed into the frequency domain1, jωαe−jωτα =Z αe−jωτα+(1+Z )q+Z δ e−jωτδe (2a) α q δe e jωq =M αe−jωτα+M q+M δ e−jωτδe (2b) α q δe e a =Vo Z αe−jωτα+Z q+Z δ e−jωτδe (2c) z g  α q δe e  In the transformed Eqs. (2), the relative time skews are continuous-valued parameters that can be estimated in the same manner as the stability and control derivatives. At first glance, it would appear that all of the unknown parameters in Eqs. (2) could be estimated from the data using standard equation-error or output-error parameter estimation in the frequency domain1,8. However, this is not true, because the relative time skew parameter for angle of attack is highly correlated with the stability derivatives, which leads to inaccurate parameter estimates. The reason for this can be understood physically - any time skew in the angle of attack data can be interpreted as phase shift in the dynamic system, which can be modeled nearly equally well by adjusting the values of the stability derivatives or by adjusting the angle of attack time skew parameter. This is called model parameter correlation, which gives rise to an ill-posed parameter estimation problem that cannot be solved accurately. Note that the control time skew parameter does not have this problem, because the control has nothing to do with the phase shift associated with the dynamic system. Consequently, control time skew parameters are well-conditioned, identifiable parameters that can be estimated accurately along with the stability and control derivatives, using standard parameter estimation algorithms in the frequency domain. The correlation between the stability derivatives and the angle of attack time skew can be broken if an accurate independent estimate of the angle of attack time skew parameter can be found. This can be done using data reconstruction, which in some contexts is called flight path reconstruction. In this approach, air flow data are reconstructed using IMU data (defined as having zero time skew) for translational accelerations, angular rates, and Euler attitude angles. The reconstructed air flow data are transformed into the frequency domain, as explained in the next subsection. Then a simple nonlinear optimization in the frequency domain can be used to accurately determine the relative time skews for the air flow data (airspeed, sideslip angle, and angle of attack). Once these relative time skew quantities are estimated accurately and separately, the resulting values are substituted into the dynamic equations (e.g., Eqs. (2)), and the remaining unknown parameters, consisting of stability and control derivatives and control time skew parameters, can be estimated accurately using a nonlinear parameter estimation algorithm in the frequency domain. The next section describes transformation of measured flight data into the frequency domain. This is followed by an explanation of the data reconstruction and nonlinear optimization in the frequency domain used to accurately estimate time skews in the air-flow data. Then, the application of parameter estimation in the frequency domain to estimate stability and control parameters and the control time skews is explained, which completes the modeling. Because all of the unknown model parameters are estimated using standard parameter estimation algorithms in the frequency domain, the uncertainties in the estimated parameters can also be computed using standard methods1. 4 American Institute of Aeronautics and Astronautics A. Transforming Data from the Time Domain into the Frequency Domain Measured flight data are transformed from the time domain into the frequency domain using the finite Fourier transform. For an arbitrary scalar signal x(t) on the time interval [0,T], the finite Fourier transform is defined by x(t)≡x(ω)≡∫Tx(t)e−jωtdt (3) 0 The finite Fourier transform can be computed very accurately for arbitrary frequencies ω using a numerical method based on Filon quadrature and cubic interpolation1,9. When a time series has non-zero endpoints, the finite Fourier transform of the time derivative of that time series requires endpoint correction terms. Applying integration by parts to Eq. (3), x(t)≡∫Tx(t)e−jωtdt 0 =x(T)e−jωT −x(0) + jω∫Tx(t)e−jωtdt (4) 0 =x(T)e−jωT −x(0) + jωx(ω) Applying high-pass filtering to a time series or assuming perturbations from an initial reference condition imposes a zero initial condition. This leaves only one term for the endpoint corrections, so that x(t)=x(T)e−jωT + jωx(ω) (5) Equation (5) was used to compute the finite Fourier transform for the time derivative of a time series that has been high-pass filtered or conditioned as a perturbation from an initial reference condition. B. Data Reconstruction and Time Skew Estimation in the Frequency Domain Data from the IMU can be used in the nonlinear translational kinematic equations of motion to reconstruct body-axis velocity components u,v,w, from which V,β,α can be computed. The required equations are1: u=rv−qw−gsinθ+ga (6a) x v= pw−ru+gcosθsinφ+ga (6b) y w=qu−pv+gcosθcosφ+ga (6c) z Airspeed, sideslip angle, and angle of attack can be computed from u,v, and w using1 V = u2+v2+w2 (7a) r β =sin−1(v V ) (7b) r r α =tan−1(w u) (7c) r Equations (6)-(7) can be used to reconstruct airspeed, sideslip angle, and angle of attack data using IMU measurements. Consequently, the reconstructed V ,β,α time series do not have relative time skew. r r r 5 American Institute of Aeronautics and Astronautics A practical problem with this data reconstruction is that IMU data for angular rates p,q,r, and translational accelerations a ,a ,a , all typically have bias errors. This causes a time-dependent drift in the V ,β,α x y z r r r reconstructions from Eqs. (6)-(7), even when the bias errors are small, because of the additive effect of the time integration. Furthermore, the initial conditions for the integrations in Eqs. (6) must be specified, and any errors in these values will bias the reconstructed data. Scale factor errors for angular rate and translational accelerometer measurements are negligible in practice, and any errors in the Euler angles φ and θ will be mitigated by the trigonometric functions. Most IMU provide Euler angle time series (or equivalent quaternion data) that have been corrected for drift, using magnetometer measurements, for example. For the data reconstruction in Eqs. (6), initial conditions must be specified. The initial conditions can be computed from initial measured values for the air flow data, or a mean of measured values in a steady initial flight condition. Initial conditions for the body-axis velocity component reconstructions are then computed from u(0)=V(0)cosα(0)cosβ(0) (8a) v(0)=V(0)sinα(0) (8b) w(0)=V(0)sinα(0)cosβ(0) (8c) In practical situations, particularly when the initial flight condition is not steady, this approach to finding the initial conditions for Eqs. (6) will be approximate. However, any errors in the initial conditions will be shown to have little effect on the modeling results. From the foregoing discussion, it is clear that reconstructed V ,β,α time series can have significant bias and r r r drift errors, from the approximate initial conditions used in Eqs. (6), as well as the biases in the IMU data for angular rates p,q,r, and translational accelerations a ,a ,a . However, if the modeling is done in the frequency x y z domain, where the bias and drift are removed for other reasons (discussed next), then this method of reconstructing air flow data can be used very effectively. When transforming time-domain data into the frequency domain, the constant bias and drift (also called the trend, or the best-fit linear function of time) are always removed prior to applying the finite Fourier transform. This processing is done to avoid leakage from relatively large low-frequency components that can pollute the frequency- domain data at low frequencies of interest1,10,11. It follows that bias and drift errors in the air flow data reconstructions are of no concern when the analysis is done in the frequency domain, because the bias and drift of every time series are removed anyway. In essence, the reconstruction is done in the time domain, and the errors incurred in doing that are discarded prior to transformation into the frequency domain. At that point, standard frequency-domain modeling methods can be applied. Note that this approach conforms with the practical reality that angular rate measurements and translational accelerometer measurements typically have small but non-zero bias errors, and negligible slope or scale factor errors1. To estimate the time skew in the angle of attack measurement, a simple nonlinear optimization can be done in the frequency domain. The model equation relating the measured angle of attack and the reconstructed angle of attack in the frequency domain is simply α =α e−jωτα (9) r A nonlinear optimizer such as Gauss-Newton implemented in the frequency domain1 can be used to find the best estimate of the single unknown parameter τ . Equivalently, the parameter estimation problem can be set up as an α output-error parameter estimation problem in the frequency domain, which is also solved using a nonlinear optimizer. This parameter estimation problem is very well conditioned and the solution is obtained very rapidly using Gauss-Newton optimization in the frequency domain. A similar approach can be used to find the time skews 6 American Institute of Aeronautics and Astronautics for airspeed and sideslip angle, individually. The angle of attack corrected for time skew in the frequency domain is computed by inverting Eq. (9) using the angle of attack time skew estimate τˆ found by the nonlinear optimizer, α α =αejωτˆα (10) c and similarly for airspeed and sideslip angle. Note that this approach retains data information from the air flow data measurements and only estimates relative time skews. Once air flow angle data are corrected using the estimated time skews, the remaining aerodynamic modeling and control time skew estimation can be done in the frequency domain, as described in the next subsection. Time skews in the Euler angle data from the IMU can be estimated in an analogous way using the nonlinear rotational kinematic equations of motion to reconstruct Euler angles φ,θ,ψ, from body-axis angular rates p, q, r. The rotational kinematic equations are1: φ = p+tanθ (qsinφ +rcosφ) (11a) r r r r θ =qcosφ −r sinφ (11b) r r r qsinφ +rcosφ ψ = r r (11c) r cosθ r Reconstructed φ,θ,ψ time series have no relative time skew, because the data reconstruction is based on r r r IMU body-axis angular rate data. Equation (11c) can be omitted because Euler yaw angle does not affect dynamic modeling, and does not appear in Eqs. (11a)-(11b). Initial conditions for Eqs. (11a)-(11b) can be set to initial measured values or a mean of measured values in a steady initial flight condition. As mentioned previously, most IMU provide Euler angle time series that have been corrected for drift. However, filtering and data processing involved in producing the Euler angle time series can produce time skews relative to the angular rate data. The time skew in the Euler pitch angle data can be estimated in the same manner as described earlier, using a simple nonlinear optimization in the frequency domain. For example, the model equation relating the Euler pitch angle data from the IMU and the reconstructed Euler pitch angle data in the frequency domain is simply θ=θ e−jωτθ (12) r The Euler pitch angle data corrected for time skew in the frequency domain is computed by inverting Eq. (12) using the time skew estimate τˆ found by the nonlinear optimizer, θ θ =θejωτˆθ (13) c The same approach can be applied for the Euler roll angle. The rotational kinematic analysis can be done separately and independently from the translational kinematic analysis. Therefore, it is convenient to do the rotational analysis first, then correct the Euler angle data for the estimated time skews before doing the translational kinematic analysis to estimate the air flow data time skews. However, because the Euler angles appear only as arguments in trigonometric functions in the translational reconstruction Eqs. (6), and are not explanatory variables for the aerodynamic modeling, the rotational time skew estimation can be skipped. C. Aerodynamic Modeling and Control Time Skew Estimation in the Frequency Domain Continuing with the aircraft longitudinal short-period dynamics example described earlier, and using the angle of attack data corrected for time skew, the frequency-domain equations are: 7 American Institute of Aeronautics and Astronautics jωα =Z α +(1+Z )q+Z δ e−jωτδe (14a) c α c q δe e jωq =M α +M q+M δ e−jωτδe (14b) α c q δe e a =Vo Z α +Z q+Z δ e−jωτδe (14c) z g  α c q δe e  As mentioned earlier, the remaining unknown parameters, consisting of stability and control derivatives and control time skew parameter, can now be estimated using standard parameter estimation methods in the frequency domain. Using an equation-error formulation, Eqs. (14) can be analyzed either individually or together. The equation- error approach produces parameter estimates that minimize the least squares fit of the model to derivative information in the frequency domain for each equation1. For the pitching moment equation (14b), the equation error in the frequency domain is ε(θ)=z−y(θ) (15) where z = jωq (16a) y(θ)=M α +M q+M δ e−jωτδe (16b) α c q δe e T θ=M M M τ  (16c)  α q δe δe The Fourier transform quantities based on the data, i.e., α ,q,δ , are computed at M selected frequencies c e ω, k =1,2,,M. The number of frequencies M is selected to be much greater than the number of unknown k parameters n (n =4 in Eq. (16c)), and the frequencies are chosen to cover the frequency band associated with p p the dynamic model, e.g., [0.1, 1.5] Hz for rigid-body dynamics of a typical fighter or subscale aircraft. Note that the transform frequencies can be chosen arbitrarily (subject only to data information content limitations) using the high- accuracy numerical method described in Refs. [1] and [9]. This results in an overdetermined set of equations with complex data, z =y(θ)+ε (17) where z =jω1q(ω1) jω2q(ω2)  jωMq(ωM)T (18a) αc =αc(ω1) αc(ω2)  αc(ωM)T (18b) q=q(ω1) q(ω2)  q(ωM)T (18c) δe =δe(ω1) δe(ω2)  δe(ωM)T (18d) 8 American Institute of Aeronautics and Astronautics ω=[ω ω  ω ]T (18e) 1 2 M y(θ)=M α +M q+M δ e−jωτδe (18f) α c q δe e and ε is a vector of complex equation errors. Equation (18f) shows that the model output is nonlinearly dependent on the parameters, because of the time skew parameter τ in the exponential. This simply means that the solution δ e cannot be obtained with the algebraic pseudo-inverse calculation normally used for equation-error parameter estimation, but rather must be solved using a nonlinear optimizer, such as Gauss-Newton in the frequency domain. The solution is obtained in a straightforward manner using nonlinear optimization (typically applied for other nonlinear parameter estimation, such as output-error problems) in the frequency domain, with the least squares cost function: J(θ)= 1 z−y(θ)†z−y(θ) (19) 2 The process produces accurate estimates for the stability and control parameters, as well as the control time skew. For these estimated parameters, as well as the air flow time skew parameters (which also come from a nonlinear parameter estimation solution), the estimated parameter covariance matrix can be computed from1 Σ(θˆ)≡E(θˆ−θ)(θˆ−θ)T=σˆ2{ReS†S}−1 (20a)       ∂y(θ) S≡ (20b) ∂θ θ=θˆ where the equation-error variance estimate σˆ2 can be found from the model residuals, σˆ2 = 1 z−y(θˆ)†z−y(θˆ) (20c) ( )     M −n   p and n is the number of unknown parameters, i.e., the number of elements in parameter vector θ. The standard p errors of the estimated parameters are given by the square roots of the diagonal elements of the covariance matrix Σ(θˆ), s(θˆ )= Var(θˆ ) = Σ j=1,2,,n (20d) j j jj p Explanations of why the estimated parameter standard errors are computed in this way, and why this calculation in the frequency domain produces parameter error measures that are consistent with the scatter in parameter estimates from repeated maneuvers, can be found in Ref. [1]. It is also possible to solve this parameter estimation problem using an output-error formulation in the frequency domain, which is equivalent to finding model parameter estimates that match the aircraft response information in a least squares sense, rather than the derivative information. However, that approach is more complicated, and it was found that the equation-error formulation provided better parameter sensitivity for the control time skew parameters, and produced modeling results with high accuracy. Because the nonlinear optimization is being done in the frequency domain, all of the calculations are algebraic, and the nonlinear optimization can be carried out very rapidly on a desktop computer. 9 American Institute of Aeronautics and Astronautics Using the equation-error approach described here, each equation (for example, Eq. (14b) or (14c)) can be analyzed individually. However, model parameters can also be estimated in a combined equation-error formulation. For example, Eqs. (14b) and (14c) could be combined, and the nonlinear optimization applied to the combined frequency-domain data. This allows the estimator to find a single best value for the control time skew, which appears in both Eqs. (14b) and (14c). Experience has shown that using the output equation (14c) for a gives superior modeling results compared to z using the state equation (14a) for α, when estimating the model parameters that appear in both equations. This is related to improved parameter sensitivity and signal-to-noise ratio using the output equation (14c). Similar statements apply for lateral/directional modeling, relative to using the output equation for a versus the state y equation for β. Although the example shown here is for longitudinal short-period dynamics using dimensional stability and control derivatives, the same approach can be used for other models, including lateral/directional dynamic models. Furthermore, the approach also applies if the modeling is done using nondimensional aerodynamic coefficients and nondimensional stability and control derivatives. In that case, the measured data to be matched by the model using the equation-error formulation are nondimensional aerodynamic force and moment coefficient data, computed from flight measurements as follows1: C =(ma −T ) qS (21a) X x x C =(ma −T ) qS (21b) Z z z C = I q +(I −I )pr+I (p2−r2)−M  qSc (21c) m  y x z xz T C =ma qS (22a) Y y C =I p −I (pq+r)+(I −I )qr qSb (22b) l  x xz z y  C =I r−I (p−qr)+(I −I )pq qSb (22c) n  z xz y x  These expressions retain the full nonlinear dynamics in the equations of motion for a symmetric rigid aircraft. For local modeling, the aerodynamic force and moment coefficients computed from Eqs. (21) and (22) can be modeled using linear expansions in aircraft states and controls, e.g., qc C =C α+C +C δ (23a) X Xα Xq 2V Xδ qc C =C α+C +C δ (23b) Z Zα Zq 2V Zδ qc C =C α+C +C δ (23c) m mα mq 2V mδ pb rb C =C β+C +C +C δ (24a) Y Yβ Yp 2V Yr 2V Yδ pb rb C =C β+C +C +C δ (24b) l lβ lp 2V lr 2V lδ 10 American Institute of Aeronautics and Astronautics