Parameter Uncertainty for Aircraft Aerodynamic Modeling using Recursive Least Squares Jared A. Grauer∗ and Eugene A. Morelli† NASA Langley Research Center, Hampton, Virginia 23681 A real-time method was demonstrated for determining accurate uncertainty levels of stability and control derivatives estimated using recursive least squares and time-domain data. The method uses a recursive formulation of the residual autocorrelation to account forcoloredresiduals, whichareroutinelyencounteredinaircraftparameterestimationand change the predicted uncertainties. Simulation data and flight test data for a subscale jet transport aircraft were used to demonstrate the approach. Results showed that the correcteduncertaintiesmatchedtheobservedscatterintheparameterestimates,anddidso moreaccuratelythanconventionaluncertaintyestimatesthatassumewhiteresiduals. Only smalldifferenceswereobservedbetweenbatchestimatesandrecursiveestimatesattheend of the maneuver. It was also demonstrated that the autocorrelation could be reduced to a small number of lags to minimize computation and memory storage requirements without significantly degrading the accuracy of predicted uncertainty levels. Nomenclature a , a , a body-axis accelerometer output, g y modeled output x y z b wing span, ft z measured output C , C , C body-axis moment coefficients α angle of attack, rad l m n C , C , C body-axis force coefficients β angle of sideslip, rad X Y Z c¯ mean aerodynamic chord, ft ∆ perturbation cov(.) covariance δ control surface deflection, rad D dispersion matrix (cid:15) measurement equation error E[.] expectation operator θ model parameter vector h altitude, ft σ standard deviation I identity matrix I moment of inertia, slug·ft2 Subscripts .. J cost function 0 bias or trim value K gain matrix a, e, r aileron, elevator, rudder m mass, slug cm center of mass N number of samples T thrust p, q, r body-axis angular rate, rad/s q¯ dynamic pressure, lbf/ft2 Superscripts S wing reference area, ft2 −1 inverse S power spectral density T transpose .. t time, s ˆ estimated value V true airspeed, ft/s ˙ time derivative X regressor matrix ∗ResearchEngineer,DynamicSystemsandControlBranch,MS308. SeniorMemberAIAA. †ResearchEngineer,DynamicSystemsandControlBranch,MS308. AssociateFellowAIAA. 1of20 AmericanInstituteofAeronauticsandAstronautics I. Introduction Least-squares parameter estimation is a commonly-used method for building models from measured data becauseofanumberofappealingqualitiesincludingsimplicity,efficiency,andextensibility. Therecursive formofthisapproachisanattractiveoptionforestimatingaircraftdynamicmodelsinrealtimeasflighttest data are being collected. Such a real-time analysis can enable rapid flight envelope expansion, control law verification, input type or amplitude selection, or aerodynamic database identification. Real-time estimates could instead update adaptive control laws during aircraft faults or loss-of-control type scenarios. These algorithms are simple enough to run on low-cost systems aboard unmanned air vehicles. However, recursive least squares is generally not favored and has instead been superseded by other methods, such as in Ref. [1], for real-time estimation. One reason is that recursive least squares produces inaccurateparameteruncertaintyestimates,whichisaprimaryindicatorofestimationquality. Onesourceof error is colored modeling residuals, which violate the underlying theory and result in parameter covariances that are too small. This error causes a misleadingly optimistic interpretation of the estimation results. For these reasons, parameter uncertainty calculations are delayed until the end of the maneuver2 or are omitted entirely.3 A residual sequence is said to be colored when it is correlated in time with itself. This occurs routinely with aircraft flight test data because of model structure error, which can arise from numerous sources. For example, turbulence, unsteady aerodynamics, and aeroelastic effects all create responses that are not considered in typical modeling applications. Other neglected dynamics can include nonlinearities, cross-axis coupling, and interactions with the propulsion system. Steps taken in the data reduction and analysis can also potentially color the residuals. Reference [4] identified colored residuals as the source of the discrepancy between predicted parameter covariance and observed scatter in flight test results. The phenomenon was well-known at the time, and the standard solution approach was to multiply the parameter stanadrd errors by an arbitrary “fudge factor” of 5 to10. In Ref. [5], abatchcorrection for theuncertaintylevelsbased on thecoloring of themodel residuals was presented for maximum likelihood estimators, with least squares as a special case. This correction was made recursive for real-time estimation in Ref. [6], and was applied to frequency response estimation from input and output measurements. Thispaperinvestigatesthereal-timemethodpresentedinRef.[6]butforadifferentapplication,namely, estimating accurate parameter uncertainties of aircraft stability and control derivatives. Estimation of sta- bility and control derivatives from flight test data is a common task, and accurate results are essential because they are used in important applications.7 The aerodynamic models containing these parameters are described in Section II, along with the batch and recursive least squares algorithms. That section then explains the recursive correction for colored residuals that enables accurate real-time parameter uncertainty prediction. Section III describes the T-2 subscale jet transport aircraft used to study the methods, as well as a linearized flight dynamics simulation model. Results using simulation data and flight test data for the T-2 aircraft are then presented in Section IV, followed by conclusions in Section V. Software for the input design, real-time smooth differentiation, and batch least-squares estimation used in this work are available in a MATLAB(cid:13)R toolbox called System IDentification Programs for AirCraft, or SIDPAC.2 II. Theoretical Development II.A. Aerodynamic Modeling Nondimensional aerodynamic force and moment coefficients cannot be directly measured in flight, and are instead computed from flight test data and known quantities. Under typical simplifying assumptions, 2of20 AmericanInstituteofAeronauticsandAstronautics the nondimensional aerodynamic coefficients can be computed as2 1 C = (ma −X ) (1a) X q¯S x T 1 C = (ma ) (1b) Y q¯S y 1 C = (ma ) (1c) Z q¯S z 1 C = [I p˙−I (r˙+pq)+(I −I )qr] (1d) l q¯Sb xx xz zz yy C = 1 (cid:2)I q˙+(I −I )pr+I (p2−r2)(cid:3) (1e) m q¯Sc¯ yy xx zz xz 1 C = [I r˙−I (p˙−qr)+(I −I )pq] (1f) n q¯Sb zz xz yy xx whichareexpressedinthebodyaxesandretainthefullnonlinearityoftheaircraftmotion. Theseequations require the aircraft geometry and mass properties, thrust, and measurements of the motion. Time histories of the rotational accelerations are often obtained by smoothly differentiating angular rate data.2,8 The aerodynamic coefficients can be modeled using linear expansions in the aircraft states and controls, assuming small excursions from a reference flight condition and short durations of time. These models can be simplified to include only the most dominant on-axis effects, for example9 C =C +C ∆α+C ∆δ (2a) X X0 Xα Xδe e b C =C +C ∆β+C r+C ∆δ (2b) Y Y0 Yβ Yr2V Yδr r C =C +C ∆α+C ∆δ (2c) Z Z0 Zα Zδe e b b C =C +C ∆β+C p+C r+C ∆δ +C ∆δ (2d) l l0 lβ lp2V lr2V lδa a lδr r c¯ C =C +C ∆α+C q+C ∆δ (2e) m m0 mα mq2V mδe e b C =C +C ∆β+C r+C ∆δ (2f) n n0 nβ nr2V nδr r The bias terms, for example C , subsume all the steady portions of the signals, including aerodynamic Z0 biases and non-zero trim contributions. These terms, sometimes referred to as “nuisance parameters,” are estimated but are not of primary interest.9 The quantities multiplying the normalized state and control variables are the nondimensional stability and control derivatives. II.B. Least Squares Parameter Estimation In the equation-error formulation, which is essentially the least-squares problem, the dependent or re- sponse variables are the aerodynamic coefficients computed using Eqs. (1), the regressor or explanatory variablesarethenormalizedstateandcontrolmeasurementsinEqs.(2), andtheunknownparameterstobe estimated are the corresponding stability and control derivatives. The estimation is decoupled in that each aerodynamic coefficient can be modeled separately, instead of all together. Each analysis has the form z=y+(cid:15) =Xθ+(cid:15) (3) wherezandy aretheN×1timehistoriesofthemeasuredandmodeleddependentvariable,respectively,X isaN×n matrixofindependentvariableorregressortimehistories,θ isann ×1arrayofunknownmodel θ θ parameters, and (cid:15) is an N ×1 array of the measurement equation errors. For example using the vertical force coefficient,       a xT 1 α δ z= mq¯Sg  azz...21 , X= x...1T2 = 1... α...12 δee...12 , θ = CCZZα0  (4)       C a xT 1 α δ Zδe zN N N eN 3of20 AmericanInstituteofAeronauticsandAstronautics Whentheregressorsareallmeasured,containnoerror,andformanadequatemodelstructure,leastsquares can be applied to match the aerodynamic coefficient time histories and accuratelyestimatethe stability and control derivatives. The least-squares cost function 1 J(θ)= (z−y)T (z−y) (5) 2 represents the total generalized distance between the measured and modeled dependent variable. This cost functionisminimizedbyequatingthecostgradienttozeroandsolvingfortheoptimalmodelparametersas θˆ=(cid:0)XTX(cid:1)−1XTz =DXTz (6) where D is the dispersion matrix and is related to the information content in the data. The difference between the dependent variable and the model output v=z−y (7) are called the model residuals. The uncertainty in the model parameter estimates is quantified by the parameter covariance matrix, which simplifies in this case as (cid:104) (cid:105) cov(θˆ)=E (θˆ−θ)(θˆ−θ)T =E(cid:104)(cid:0)XTX(cid:1)−1XT(z−y)(z−y)TX(cid:0)XTX(cid:1)−1(cid:105) =DXTE[vvT]XD (8) because the regressors are assumed to be deterministic. The square root of the diagonal entries of the parameter covariance are the Cram´er-Rao bounds, which provide a lower bound on the parameter standard error. The middle term in Eq. (8) is the autocorrelation matrix of the residual, defined as   R(0) R(1) ... R(N −1)  R(1) R(0) ... R(N −2)  E[vvT]= ... ... ... ...  (9)   R(N −1) R(N −2) ... R(0) which is symmetric about the diagonal. Entries in this matrix are approximated by the discrete sample autocorrelation function of the residual N−i Rˆ(i)= 1 (cid:88)v v , i=0,1,...,N −1 (10) N i+j j j=1 taken at different lag indices i. The autocorrelation function is even, so that R(i)=R(−i). The definition used in Eq. (10) is biased because it does not consider the statistical degrees of freedom. However, this form can have less mean square error than other definitions, and the bias is negligible for typical flight test samplingratesandrecordlengths. TheparametercovarianceinEq.(8)canalternativelybeexpressedusing summations as   N N cov(θˆ)=D(cid:88)(cid:88)xiR(i−j)xTjD (11) i=1j=1 If the model residuals can be assumed to have constant variance σ2 and are uncorrelated in time, which is the conventional assumption of white residuals, then  σ2, i=0 R(i)= (12)  0, otherwise 4of20 AmericanInstituteofAeronauticsandAstronautics and the parameter covariance reduces to   N N cov(θˆ)=D(cid:88)(cid:88)xi(σ2I)xTjD i=1j=1 =D(cid:2)σ2(cid:0)XTX(cid:1)(cid:3)D =σ2D (13) It is usually necessary to estimate the residual variance from the data because accurate values from prior or repeated experiments are difficult to obtain. One technique is to use the fit error variance 1 σˆ2 = vTv (14) N whichisbasedonthemodelresidualsandisthereforedependentupontheaccuracyofthemodelstructureand parameter estimates. When the residuals are significantly colored, Eq. (13) will underpredict the parameter uncertainties, and the more general expression in Eq. (11) is needed for accurate results.5 II.C. Recursive Least Squares The previous formulation of least squares is for batch estimation, after all the data from an experiment has been collected. In the recursive version, estimates are updated as each new set of measurements are obtained, viz.2 K =D x (cid:2)1+xTD xT(cid:3)−1 (15a) k k−1 k k k−1 k D =(cid:2)I−K xT(cid:3)D (15b) k k k k−1 (cid:104) (cid:105) θˆ =θˆ +K z −xTθˆ (15c) k k−1 k k k k−1 where the subscript k is used to denote the sample index. With each sample, a new set of regressors x k and a new measurement z are obtained. This information is used to compute a gain matrix K and the k k dispersion matrix D . The parameter estimate is then updated, based upon the previous estimate and the k weighted residual. Assumingwhiteresiduals, theparametercovarianceinEq.(13)isupdatedusingtherecursivedispersion matrix and a recursive estimate of the fit error variance σˆ2 = 1 (cid:0)v2+v2+...+v2 +v2(cid:1) k k 1 2 k−1 k (cid:18) (cid:19)(cid:18) (cid:19) = k−1 1 (cid:0)v2+v2+...+v2 (cid:1)+ 1v2 k−1 k 1 2 k−1 k k (cid:18) (cid:19) k−1 1 = σˆ2 + v2 (16) k k−1 k k Thereisasubtledifferenceinhowthebatchandrecursivemethodscomputetheparameteruncertainties. The batch algorithm uses all the data to estimate one estimate of the model parameters, from which the residuals and parameter uncertainties are computed. The recursive algorithm uses the most recent sample to update the parameter estimates, but does not recompute the residual history using these new estimates. Althoughtheparameterestimateswillmatchattheendofthemaneuverusingeitheralgorithm,theresiduals, and thus parameter uncertainties, are not expected to be the same. Additionally, when there is little information content in the data, such as early in the data record, parameter estimates are poor, which then make the residuals and parameter uncertainties also incorrect. Accuracy is sacrificed here for the sake of real-time estimation. Colored residuals will also increase the error in the recursive parameter uncertainties computed in this way. II.D. Recursive Parameter Covariance With Colored Residuals The batch algorithm for parameter covariances considering colored residuals given in Eq. (11) can be made recursive, as described in Ref. [6]. The first step is to write the residual autocorrelation in a recursive 5of20 AmericanInstituteofAeronauticsandAstronautics form, as (cid:18) (cid:19) (cid:18) (cid:19) k−1 1 Rˆ (i)= Rˆ (i)+ v v (17) k k k−1 k k−i k Equation (17) produces the same values as the analogous equation given in Ref. [6], but is expressed in a different, and perhaps simpler, form which was derived in a similar manner as Eq. (16). The residual autocorrelationfunctiondependsonpastautocorrelationvaluesandpastresiduals,whicharebothfunctions of past parameter estimates. This expression is not expected to give the same numerical values as a batch estimate, although they should be similar once parameter estimates begin to converge. The next step is to write the parameter covariance considering colored residuals in a recursive form. By expanding the summand and regrouping terms, Eq. (11) can be rewritten as (cid:34)N−1 (cid:35) cov(θˆ)=D (cid:88) R(i)Λ(i) D (18) i=0 where  N (cid:88)j=1xjxTj, i=0 Λ(i)= (19) N(cid:88)−ixi+jxTj +xjxTi+j, i>0 j=1 where the index i here represents the relative lags between samples. The first definition for Λ(0) represents thediagonalproductsinvolvingthezero-lagautocorrelation,andtheremainingdefinitionisfortheproducts involving non-zero lag indices. The recursive parameter covariance update is then (cid:34)k−1 (cid:35) cov(θˆ )=D (cid:88)R (i)Λ (i) D (20) k k k k k i=0 where  Λk−1(i)+xkxTk, i=0 Λ (i)= (21) k Λ (i)+x xT +x xT , i>0 k−1 k−i k k k−i and can be used for accurate real-time estimation of the parameter uncertainty levels. Note that as more data are collected, the summation in Eq. (20) operates over a growing number of indices. Morecomputationsandmorememorystoragearerequiredastheautocorrelationmatrixincreasesin size. Theautocorrelationfunctioncan,however,betruncatedafterafinitenumberoflagtermsbecauseonly proximateresidualsarecorrelatedintypicalflighttestdata. Computationandmemorystoragerequirements would then be reduced and bounded. This approximation can be visualized as cropping the autocorrelation function around the zero-lag index, or equivalently by thinning the autocorrelation matrix in Eq. (9) about its diagonal. Retaining only the first n lags, the truncated parameter covariance is τ cov(θˆ )=D (cid:34)(cid:88)nτ Rˆ (i)Λ (i)(cid:35)D (22) k k k k k i=0 Retaining only one term, i.e. n = 0, corresponds to the conventional uncertainty estimate which assumes τ white residuals. There is no procedure for selecting how many lags to maintain. For batch analysis, N/5 lags was suggested.2 For the real-time analysis, Ref. [6] used 11 lags. This number also varies according to aircraft scale, dynamics of interest, and sampling rate. III. T-2 Aircraft The estimation methods were applied to the airplane known as the T-2, which is a 5.5% dynamically scaled version of a generic transport aircraft. The T-2 has retractable tricyle landing gear, two jet engines 6of20 AmericanInstituteofAeronauticsandAstronautics Figure 1. T-2 airplane (credit: NASA Langley Research Center) mountedunderthewings,andaconventionaltailconfiguration. Figure1showsaphotographoftheairplane in flight, and Table 1 lists airplane geometry and nominal mass properties. Control surfaces for the T-2 are left and right ailerons, left and right inboard and outboard elevators, upper and lower rudders, left and right inboard and outboard trailing-edge flaps, and left and right inboard and outboard spoilers. This configuration amounts to 16 control surfaces, all of which can be moved inde- pendently. For the maneuvers examined here, only the elevators, ailerons, and rudders were deflected. The individual elevator surfaces were moved together as a single elevator surface, and similarly for the rudders, whereas the ailerons were moved asymmetrically, in the conventional way. Trailing-edge down is consid- ered positive deflection for wing and elevator surfaces, and trailing-edge left is positive for rudder surfaces. Definitions of control surface deflections are 1 δ = (δ +δ +δ +δ ) (23a) e 4 elo eli eri ero 1 δ = (δ −δ ) (23b) a 2 ar al 1 δ = (δ +δ ) (23c) r 2 ru rl The T-2 was outfitted with a variety of research-quality hardware. A micro-INS provided tri-axial translational accelerometer measurements, angular rate gyroscope measurements, estimated Euler angles, andGPS-derivedpositionandvelocity. Airdataprobesonbothwingtipsmeasuredangleofattack,sideslip angle, static pressure, and dynamic pressure. Measurements from static pressure sensors and outside air temperature sensors were used to compute air density and altitude. Engine speeds were measured and used in a thrust model which was identified from ground test data and augmented with adjustments for ram drag identified from flight data. Potentiometers on the rotation axes of the control surfaces measured controlsurfacedeflections. Masspropertieswerecomputedbasedonmeasuredfuelflow,pre-flightweightand balance,andbody-axisinertiameasurementsdoneonthegroundfortheaircraftwithoutfuel. Measurements were sampled at 200 Hz. The aircraft was flown using the NASA Langley AirSTAR (Airborne Subscale Transport Aircraft Re- search) flight test facility. This capability includes a mobile control station, which houses telemetry equip- ment, computational hardware, a pilot station, and stations for flight test personnel and researchers. A research pilot flies the aircraft from inside the control station, using synthetic vision drawn from telemetry dataandadatabaseofthelocalterrain. Controldeflectionsaregeneratedbythepilotandtheground-based flight control system. The AirSTAR flight control system has the capability to inject arbitrary inputs at the actuators, just before position and rate limiters. These inputs were engaged by the pilot pressing and holding a button on the thrust levers. Telemetered data, downsampled to 50 Hz, is available to researchers inside the control station through the internal ethernet network. For more information on AirSTAR, see for example Ref. [10]. Fortheresultsexaminedhere, theexcitationinputsappliedwereorthogonalphase-optimizedmultisines. This type of input was developed at NASA Langley11,12 and has been used for highly-efficient excitation of 7of20 AmericanInstituteofAeronauticsandAstronautics different types of vehicles in a variety of unusual flight conditions. Each multisine input has the form (cid:18) (cid:19) (cid:88) 2πk u(t)=a a sin t+φ (24) k T k k where a is the aggregate amplitude, a are the normalized relative sinusoid amplitudes, T is the excitation k record length, and φ are the phase angles. Choice of the record length determines the available excitation k frequencies ω =2πk/T. Selection of the sinusoid amplitudes and harmonic indices k can be used to tailor k the power spectrum of the inputs. The first few harmonic indices are usually discarded to ensure that the data record contains repeated cycles of each excitation frequency. Phase angles are optimized for minimum relativepeakfactor(RPF)ofthewaveformtocreatesmallperturbationresponses. Theaggregateamplitude is used to scale the resulting waveform to achieve good signal-to-noise ratios. Multisine input parameters used here are given in Table 2. IV. Results IV.A. Simulation Results AsimulationmodeloftheT-2longitudinalshort-perioddynamicswasfirstusedtotestthemethod. The state-space perturbation model was (cid:34) (cid:35) (cid:34) (cid:35)(cid:34) (cid:35) (cid:34) (cid:35) ∆α˙(t) = mq¯SVCZα 1 ∆α(t) + mq¯SVCZδe ∆δ (t) (25a) ∆q˙(t) qI¯Syyc¯Cmα qI¯Syyc¯2c¯VCmq ∆q(t) qI¯Syyc¯Cmδe e       ∆α(t) 1 0 (cid:34) (cid:35) 0 ∆α(t)  ∆q(t) = 0 1  + 0 ∆δ (t) (25b)     ∆q(t)   e ∆a (t) q¯SC 0 q¯SC z mg Zα mg Zδe The reference flight condition simulated was straight and level flight with 134 ft/s airspeed, 1370 ft altitude, and 4.8 deg angle of attack. Values for the stability and control derivatives were taken from a separate analysis using flight test data in calm air, and are given as the true values in Table 3. The simulation model was excited using the elevator multisine specified by Eq. (24) and given in Table 2. Angular accelerations were computed using a fixed-lag smooth differentiation method.8,2 The colored noise added to the measurements was simulated by summing together two noise sequences: a wide-band noise contribution from 0 Hz to the 25 Hz Nyquist frequency, and a band-limited contribution from 0–2 Hz. The wide-band component was realized using a normally-distributed pseudorandom number generator. Signal-to-noise ratios for the elevator, angle of attack, pitch rate, and vertical acceleration were 40, 12, 30, and 40, respectively. These values were also determined from flight test data in calm air. The band-limited component was realized by passing a different white Gaussian noise sequence through a fifth- order Chebyshev filter with a 2 Hz corner frequency. This method of generating colored noise is known to be representative of typical flight test data.4,5 This range is appropriate for this scale of aircraft because it extendsslightlybeyondtherigidbodydynamicsat1.1Hz. Thecolorednoisethereforecontainsthedynamics of interest and the bandwidth of the input, and is expected to color the modeling residuals. Simulations were run for band-limited noise set between 0% and 20%, in 5% increments, where the percentage pertains to the root mean square of the signal variation. Note that 0% band-limited noise means that only wide- bandGaussianwhitenoisewasaddedthemeasurements,andthatthemodelresidualsshouldthenbewhite. Band-limitednoiselevelsgreaterthan20%werenotattemptedbecausethesearegenerallynotrepresentative ofaircraftsystemidentificationscenarios. Simulationswererepeated250timesforeachlevelofband-limited noise using unique noise sequences each time. A summary of parameter estimation results, obtained at the end of the maneuver and averaged over all 250 simulation runs for each level of band-limited noise, is given in Table 3. For 0% band-limited noise, the parameterestimatesareaccurate, parameteruncertaintiesaresmall, andresultsareinstatisticalagreement withthetruevalues. Asthelevelofband-limitednoiseincreases,themeanestimatesbecomemorebiasedand themeanparameteruncertaintiesincrease. Forallcaseswithband-limitednoise,thecorrecteduncertainties match the observed scatter in the parameter estimates well, and does so more closely than the conventional uncertainty estimates that assume white residuals. 8of20 AmericanInstituteofAeronauticsandAstronautics In the remaining discussion, results will be illustrated for the 20% band-limited noise case, because it represents a somewhat exaggerated scenario of what is usually encountered. Of these results, the C Zα estimates will be highlighted, because it is an important parameter that is usually estimated well. For this case, theconventionalstandarderrorwastoosmallbyafactoroffour, whereasthecorrectedstandarderror was consistent with the observed scatter in the estimate. The measurements from the last simulation run are shown in Fig. 2. The recursive least-squares fits for C and C in this simulation example are shown in Fig. 3. The fits Z m hadcoefficientsofdeterminationof0.93and0.87,respectively,whichindicategoodfitsconsideringthislevel of noise. The model residuals are shown in Fig. 4, first as time histories and then as power spectra. The timehistoriesshowthattheresidualissmall,andcontainsbothcoherentandincoherentcontent. Thepower spectra show the the larger band-limited noise, its roll off at 2 Hz, and the wide-band noise floor at higher frequencies. These characteristics are particularly clear for the vertical force coefficient. Some amplification of the residual around 10 Hz is visible in the pitching moment residual power spectrum, which is due to the real-time differentiation of angular rate data. Overall, these plots agree with those shown later for the flight test results and are similar to those reported previously,4,5 which indicates a realistic simulation. The autocorrelation estimates for the C coefficient residual are shown in Fig. 5. For clarity, only Z positive lags are shown since autocorrelation is an even function. Figure 5(a) shows the time history, with a decreased temporal resolution for additional clarity. Figure 5(b) shows the final recursive estimate, along √ with the batch estimate and the ±2σ uncertainty bound of the autocorrelation, given2 as Rˆ(0)/ N. If the residual were white, the autocorrelation would have a large peak at the zero-lag index equal to σ2 and be zero otherwise. Because there are lagged autocorrelations outside the ±2σ bound, this is not the case and the residual is colored. Most of the residual coloring is due to the large autocorrelations near the zero-lag index. Note the similarity between the final recursive estimates and the batch estimate, particularly for the higher autocorrelations. This indicates that not much accuracy was scarified by using a recursive method. The evolving time history of estimation results for the model parameter C , averaged over all 250 Zα simulation runs, is shown in Fig. 6(a). On average, the parameter estimate converged by 2 s, which was 1.5 s after the elevator excitation began. By 6 s, there was enough data information content and good enough parameter estimates that the corrected parameter uncertainty envelope enclosed the observed scatter in the parameter estimate. This means that on average, the recursive estimate is slightly conservative, which is safer in modeling applications than having it too optimistic. The conventional estimate of the parameter uncertaintywastoosmallbyafactorof3.3,andwasnotinstatisticalagreementwiththetruevalue. Similar time histories were observed for the other model parameters. Batch and recursive estimates, obtained at the end of the maneuver, are compared for conventional and corrected parameter uncertainties in Fig. 6(b). In both cases, the corrected uncertainties are in statistical agreement with the true value, whereas the conventional estimates are not. There was a 1% difference in the size of the corrected error bounds, again meaning that not much accuracy is sacrificed by using a recursive technique. Figure 7(a) shows the impact of reducing the number of lags retained in the residual autocorrelation. Using only the zero lag, the parameter uncertainty is too small. Only two lags were required for statistical agreementwiththetruevalueinthisexample,althoughmorelagswereneededinotherrunshavingthesame level of band-limited noise. For roughly 10–200 reatined lags, the parameter covariance overshoots its final value by about 12%. A safe number of lags to retain would be about 50. This would reduce computation and memory, safeguard against selecting too few lags, and provide a conservative estimate of the parameter covariance. Figure 7(b) illustrates the autocorrelation time history with 50 lags. A timing test conducted on a 2.6 GHz Intel Core i7 laptop computer using MATLAB(cid:13)R R2015a showed that 20% of the available 50 Hz frame rate was needed to compute the full autocorrelation function in real time for this 12 s maneuver, but only 4% of the frame rate was needed when the autocorrelation was truncated to 50 lags. Note that thisnumberisdependentupontheaircraft, samplingrate, andthedynamicsofinterest. Forexample, afull scale aircraft and/or lower sampling rates would require many fewer lag terms, and thus less computation, to achieve the same level of accuracy. These results were obtained using degraded residuals from intermediate parameter estimates that were not fully converged, and an approximated recursive correction for colored residuals. In the past, these two factors have been primary obstacles to successfully applying recursive least squares in aircraft system identification. ThemethodpresentedinRef.[6]andusedhereisasimpleandeffectivemethodofmitigating those previous challenges. 9of20 AmericanInstituteofAeronauticsandAstronautics 2 δe, 0 deg -2 -4 8 α, 6 deg 4 2 0 20 q, 10 deg/s 0 -10 -20 -0.5 az, gunits -1 -1.5 0 2 4 6 8 10 12 Time,s Figure 2. Measurements from simulation data (20% band-limited noise, Run 250) 0 Data Model -0.2 CZ -0.4 -0.6 -0.8 0.2 0.1 Cm 0 -0.1 -0.2 0 2 4 6 8 10 12 Time,s Figure 3. Recursive fitting using simulation data (20% band-limited noise, Run 250) 10of20 AmericanInstituteofAeronauticsandAstronautics