Maximum likelihood identification and optimal input design for identifying aircraft stability and control derivatives

This section will discuss several of the previous identification techniques which have been used to process flight data. This discussion will bring out the similarities that exist among these methods, and mention the aircraft flight data to which the methods have been applied. The next section, then, will provide a detailed explanation of the Maximum Likelihood Technique to be applied later on to flight test data containing gusts.

3.1 Previous Identification Methods Although a large number of identification methods have been used in the past, only some of the more common methodswhich are currently in use will be described.

3.1.1 Time Vector Method The time vector methods for derivative identification is derived from the tlme-lnvariance of the amplitude and phase relations between the state variables (degrees of freedom) of an exponentially damped second order system and the derivatives and integrals of the state variables. This invarlance is used to determine the values of the amplltude-phase relations, thereby determining the aircraft stability and control derivatives (Ref. i).

When more than one state variable (degree of freedom) is involved in the transient response, and there is a common natural frequency, the instantaneous value of any one state may be readily determined if the characteristics of any one of the motions are known, along with the amplitude ratio and phase angle relative to the characteristic motion. The time invarlance of the amplitude ratios and their phase angles permits the representation of any one of the linearized equations of motion as vectors. The properties of these vectors, plus the requirement that the vector sum of the quantities in any one equation equal zero, makes possible the determination of two unknown derivatives in any one equation.

As Ref. 2 points out, the time vector method has the principal

disadvantage that it can only be applied to control-fixed, transient- oscillation aircraft responses with damping ratios less thanZ.3. Further- more, the successful application of the time-vector method is highly de- pendent on the operators' skill.

3.1.2 Analog-Matching Methods The analog matching technique is actually an output error method since it strives to iteratively minimize the errors of the various responses through operator manipulation of the values of the stability and control derivatives.

It is often used as a backup method for validating the more modern digital techniques. However, there are several disadvantages to the analog match- ing technique. First, the method works most successfully only when a single control surface is moved at a time and then only when the maneuvers are simple.(Ref. I). Second, when the maneuvers are made with a stability augmentation system or other form of dependent control input, the data is difficult to analyze. Finally, this method is extremely time consuming, even in face of the fact that recent procedures, through the use of hybrid computers, has reduced the time considerably. For example, the time in- volved in analyzing a lateral-directional flight maneuver, from receipt of flight data to final results, is approximately four hours (Ref. 2); the analog matching technique is also extremely susceptible to uniqueness problems since the success of a data analysis is very dependent upon the type of control maneuvers used. In such a case the skill and knowledge of the operat¢,r would play an important part in determining the ultimate success of the analys:is.

The analog matching technique has been used by the Air Force Flight Test Center (Ref. 3), the Naval Air Test Center (Ref. 4), and the NASA Fli_ht Test Center (Ref. 5) for the F-104. X-15, B-70, HL-IO, M2/F3, X-24 and PA-30 aircraft.

Most of the remaining identification techniques, almost all of which

require the use of a digital computer, can be classified as either

i. Equation error methods,

2. Output error methods, or

3. Advancedmethods.

These methodsdiffer by (i) the performance criterion that they are

developed from, (2) the kinds of estimates they produce, and (3) the problems

to which they can be applied.

3.1.3 Equation Error Methods

Equation error methods (Ref. 6) assume a performance criterion that minimizes the square of the equation error (process noise). All of these methods are basically least squares techniques and, in general, it is necessary to measure all the response variables and their derivatives. The procedure is to express the stability and control derivatives as functions of the measured responses using the equations of motion. This results in n or more linear equations in n unknowns. For those cases where the time derivatives are not measured, various "method functions" are used to operate on these equations (take time derivatives, Laplace or Fourier transforms, etc.) to obtain equations that are linear in the unknown stability and control derivatives. Since these methods do not allow for measurement errors (instru- mentation errors), they result in biased estimates when this type of error does exist. The principal use of these methods are as start-up techniques for the output error and advanced methods.

The equation error methods have been used or are being used by Cornell Aeronautical Laboratory (Ref. 6), Air Force Flight Test Center (Ref. 3), and Delft University of Technology (Ref. 7).

3.1.4 Output Error Nethods Output Error Methods (Refs. 8 through 17) minimize the square of the error between the actual system output and the output of a model used to represent the actual system. This method assumes measure- ment noise but no process noise. Typical output error methods include Newton-Raphson, Gradient methods, the Kalman Filter (without process noise), and modified Ne_on-Raphson, differential correction, and quasillnearization (all three of which are the same method).

The modified Newton-Raphson method has been used extensively in flight test applications for the past several years. It is the one method that has been used on an operational basis and for which the most experience exists.

This method has been or is being used by (among others): (1) the NASA Flight Test Center (Ref. 5) on the X-24, X-14, XB-70, 990, HL-10, M2/F3, Jet Star and PA-30 aircraft; (2) the NASA Ames Research Center (Reference 18 ) on the LearJet, XV-5, 990 and the C-8 aircraft; and (3) the NASA Langley Research Center (Refs. 19,20) on the XC-142, Navlon and F-4 aircraft.

(NASA Langley program has automatic update of the weighting matrix based on the maximum likelihood criterion.)

The principal disadvantages of the output error methods is that, because they do not include process noise in their performance criterion, the results degrade when process noise (gusts, modelling errors) exists. This may result in the computer program not converging or in estimates that have large vat- lances or poor estimates (Ref. 21). However, as long as these methods are applied to linear flight regions, or where the form of the equations is known, or where gusts do not exist, they work very well.

3.1.5 Advanced Methods The most general identification problem is one of extracting stability and control derivatives, for non-linear aircraft models, from flight data containing both measurement and process noise. The one advanced technique that has demonstrated the capability of extracting stability and control derivatives from flight data under these circumstances is an implementation of the maximum likelihood criterion. This numerical algorithm, developed by SCI, is a combination of three steps: (i) Kalman filtering to estimate the states and generate a residual sequence, (2) a modified Newton-Raphson algorithm for the parameter estimates, and (3) an algorithm to estimate the noise statistics (means and variances of the measurement and process noise).

The details of the numerical method are outlined in the next section.

The success of the SCI maximum likelihood technique can be attributed to several important attributes of this method: I. It does not require a priori knowledge of the process noise covariance, measurement noise covariance or the initial parameter estimate covarlance. These covarlances are determined as part of the identification procedure.

2. When process noise does not exist, this method simplifies to the modified Newton-Raphson output error method (although_rlth a spe- cific weighting matrix).

3. When no measurement noise exists (an unlikely event) this method simplifies to the least squares equation error method.

4. The Cramer-Rao lower bound on the covarlance of the error in the stability and control derivative estimates are obtained as part of the algorithm.

5. The minimum mean-square aircraft state variables (response variables) are obtained as an integral part of the algorithm. It is not re- quired, however, that initial state estimates be supplied.

The following section gives a detailed derivation of the Maximum Likelihood Identification Method.

IV MAXIMUM LIKELIHOOD (ML) IDENTIFICATION The notion of the maximum likelihood estimate which was introduced into statistics by R. A. Fisher in 1906 is based on a relatively simple idea. Assume that the outcome Z of an experiment depends on an unknown parameter e. We want to infer the best value of e from the observation Z. One answer is to choose that value of 8 which makes the observed value Z the most probable one to have occurred. This can be rigorously stated as: choose e to maximize the conditional probability of Z, given a value of e ; i.e.

^

e = max p(zlo)

A where 0 is the maximum likelihood estimate of 0 and p(ZlO) is the conditional probability of Z, given O. The same estimate is obtained by maximizing log p(ZlO) which is known as the likelihood function.

The above basic idea can be carried over to linear and nonlinear dynamic systems, with process and measurement noise, but the details of the application become quite involved. In practice, there are two major pro- blems in obtaining ML estimates for dynamic systems. These are: i. Deriving an expression for the likelihood function, and 2. Maximizing the likelihood function with respect to the unknown parameters.

These two problems are elaborated upon further. The likelihood function is the logarithm of the joint probability density of the observations given the parameters. If the observations are independent, the joint probability density function is easily written down since it is just the product of the 2O

probability densities of each observation given the parameters. The deriva-

tion of the likelihood function becomes muchmore difficult when the obser-

vations are correlated. This is necessarily the case for dynamic systems

with randominputs since the state at any time is correlated with the state

at all the previous times. In the next section, it is shownhow the likeli-

hood function for a dynamic system can be derived in a simple form using a

Kalmanfilter and the resulting white noise innovation sequence. This is

shownschematically in Figure 4.1.

The second problem in obtaining ML estimates is a computational one.

Generally, the likelihood function is highly nonlinear in terms of the para-

meters. For finite data lengths, it is also known to have several local

maxima. In the case of dynamic systems, certain differential equation

constraints have to be Satisfied. The choice of a suitable search algorithm

is very important for the successful application of ML identification.

The maximum likelihood identification method, as implemented by SCI, is an adaption and extension of the recent work of Astrom (Ref. 22),

Kashyap (Ref. 23) and Mehra (Ref. 24). It is capable of solving the most general

MeasurementNoise

Input J °a lOut u

'+ __J System I _+k._,l i I

Measure-f "_ Process[ ] I l +

Performance _2 Noise Output Measurement I / _ Inno-

ment _

Criterion

Noise

I ....... kJ va

I I Optimum Filter _ t- _) with z I - I Adjustable I Parameters

t

I Optimizer FIGURE 4.1 IMPLEMENTATION OF MAXIMUM LIKELIHOOD ESTD_TOR

identification problem, Including (I) systems governed by non-linear differential equations of motion, (2) the presence of additive random process noise in the equations of motion, and (3) random disturbances corrupting the measurements of the system inputs and outputs. In this section the development of the theory Justifying the use of the likelihood function as an optimizing criterion for identification is developed and the numerical algorithm for implementing the maximum likelihood identification method is outlined. The development is first given for linear systems and then extended to nonlinear systems.

4.1 Linear Systems: Consider the llnearized aircraft equations of motion

(4.1)

where x(t) ffi n x i state vector (p, q, r, u, v, w, etc.)

u(t) - p x I input vector (6e, 6 a, 6r) w(t) = q x i vector of random forcing functions Let the measurement equations be where y(t) = Hx(t) + Du(t) + v(t) y(t) = r x i output or measurement vector (4.27 v(t) = r x i vector of random measurement errors and E { w(t) } = O, E { w(t)wT(T) } = Q 6(t - z) where 8(t - z) is the Dirac Delta function.

Z { w(t)vT(z) } = 0 E { v(t) } = 0, E (v(t)vT(T) } = Z6t, T.

It is assumed that the structure of the model is known. The vector of unknown parameters in F_ G, F, H, Q, R and x(o) are denoted by 0. Thus 0 includes all the unknown stability and control derivatives, noise variances and the initial states.

If a sequence of observations y(1), ........... ,y(N) is made of the system state with noise, the maximum likelihood estimate of O, following the idea stated earlier, is given by ^

e = max p(YNIS) (4.3)

where YN = {y(1), ..... , y(N)}. With successive applications of Bayes Rule, an expression for p(YNI8) can be derived as

p(YN/e) ,,, p(y(1),.., y(N)/e)

'= p(yCN)[YN-I' 8)P(YN-1 IO) = p(y(N)IYN_I, ®)p(y(N-l) JYN-2' ' O)P(YN-2JS) N

" e)

J=l

Since the logarithm is a monotonic function, the maximum likelihood estimate can also be written as

(4.4)

0 " max [log P(XNl0)] = max log p(y(J)IYj_l,8 e 8 _=i where log p(YN[®) is the likelihood function.

If x(o), wCt) and v(t) are normally distributed, p(yCJ)iYj_l,e)'--

will also be normal and can be uniquely determined by computing the mean and covariance. Therefore define

(4.s)

E{y(J)IYj_1, e} A= y(jlj_l ) and

(4.6)

cov(y(J)IYj_l,O} = Z((y(J) - y(jlj-l)) (y (j ) -Y(jlJ -I) ) T}

B(J ]J-l)

With these assumptions, the term log p (y(J)[Yj_I,O) can be written as -I log p(y(j)IYj_l O) = Const.- -_-(y(j)-y(jlj-l))T B(j[j-l)(y(j)-y(jlj-l)) - log IB(J]j-'t)] (_.7) The problem of determining the maximum likelihood estimate has now become one of finding a way of calculating the conditional mean, y(J lJ-l), and the error covariance , B(JlJ-I). These quantities, however, are precisely the output of a Kalman Filter (Ref. 25) state estimator given 8. This filter is designed to recursively process measurements one at a time, and, at each point produce the minimum variance state estimate based on all the data received up to that point.

The Kalman filter prediction and update equations can be derived as follows : For a more rigorous derivation of the Kalman Filter, see either Kalman (Ref. 26) or Kailath (Ref. 27). Also note that the conditioning on G has been omitted from the equations to simply notation.

Initial Conditions: The Kalman filter is started with a priori state estimate 3(010) and covariance P(010 ).

The state prediction is done using the equations of motion and state update is done using the measurements.

Prediction Equations: At time (J-l), the Kalman filter has a state estimate _(j-llj-i ) and covariance P(J-IIJ-I). It is required to predict the state at time J. The resulting state estimate is denoted by _(J lJ-l) and has a covariance of e(j lj-l). The relationship be- tween the updated and the predicted estimates can be obtained by taking conditional expectations on both sides of Eq. (4.1) and interchanging the operations of expectation and differentiation. This gives d x(tlj_l) _ (tlj-l) + Gu(t) dt =

(4.8)

(J-l)_< t < j

where the predicted value of the white noise w(t) based on previous infor- mation is equal to zero.

A The covariance equation for the predicted estimate x(t/J-l) can be obtained by subtracting Eq. 44.8) from Eq. 44.1) and using the covariance propagation equation derived in Bryson and Ho (Ref. 25).

d

d7 P(t[j-1) = F P(t[j-1) + P(tlJ-I)F T + rQrT

(4.9) Update Equations: The update equations for the Kalman filter can be derived using a well-known property of the conditional normal distribu- tions (Ref. 28), viz.,

E{ alb } ffi 7+ PabPb (b - b) (4.10)

pT

c°v(alb) =Paa - PabPb ab 44.11) where a and b are normal random variables with E{a} = a , E{b} = coy(a) = Paa ' coy(b) = Pbb E((a-_) (b-b) T } = Pab Replacing a by x(J) with mean a = x(JlJ-l) and covariance p = p(J l j-l) and replacing b by y(J) with mean aa = y(J[J-1) = Hx(J [J-l) and covariances

(4.12)

Pbb = HP(J[J-1) HT + R Pab = P(j[j-1) H T (4.13) we obtain, A

(4.14)

x(J[J) = x (J[j-1) + K(J) (y(J) - R_(j[J-1))

(4.15)

K(J) - P(JlJ-I)H T (HP(JlJ-I)H T + R) -I and (4.16) P(JlJ) = (I - K(j)H) P(JlJ-I) The quantity (y_) - y_[J-l) ) represents the new information brought forth by the measurement y_). It is known as the "innovation" sequence and has been shown to be zero mean, Gaussian and white (Ref. 29) Denoting the innovations by u_), the likelihood function can be written as N (4.17) _ 12 I { vT(j)B-IQ[j-1) u(j) + ipglBQ[j-1)[} log p(YNle) ffi J=l where B(J[j-1) = HP(J[J-I)H T + R (4.18) A The maximum likelihood estimate 0 is obtained by maximizing (4.17) with respect to O, subject to the constraints in equation (4.8)-(4.9), (4.14)-(4.16). This is a very difficult optimization problem. An approximation suggested in Ref. 24 simplifies the problem tremendously. It is assumed that the filter gain K(J) and covariance B(j lj-l) have reached constant values K and B. The vector 0 of un- known parameters is now defined to include (in addition to F,G) K and B instead of Q and R. Reference 57 gives a detailed derivation of the relation between K, B and Q, R. Then N log p(YN]8) = 1 l(vT_)B-iv(j) + log[B[]

- _ (4.19)

j=l

Maximizing (4.19) over B, produces N (4.20)

J'=!

A where e is the ML estimate of unknowns of F, G and K.

It is given by the root of the equation N I vT(j)B I _ = O.

(4.21) J=l where _(J)_a is calculated from Eq. (4.8) - (4.18). The root of (4.21) ^ is found by a Newton-Raphson iteration. Once a is obtained, R and Q can be obtained from equations (4.9) and (4.18). In this way the non-linear constraints imposed by the equations (4.9), (4.15), (4.16) and (4.18) are avoided during optimization.

4.2 Nonlinear Systems The approach to obtaining the maximum likelihood parameter estimates for nonlinear models is conceptually similar to that for linear models.

Consider a nonlinear dynamic systemmodel of the form

x(t) = f(x(t), e, u(t)) + Fw(t) (4.22)

y(t) = h(x(t)) + v(t) (4.23) where f(') and h(') are n x i and r x i vectors of nonlinear functions.

Also, w(t) and v(t) are Gaussian white noise sequences with zero mean and covariances Q and R.

The evaluation of the exact maximum likelihood estimate involves the calculation of the conditional probability p(y(J)IYj_I,8) as in the linear model case. This would require an optimal nonlinear filter, which, to date, is computationally unfeasible since a complete description of p(y(J)[Yj_1,8) requires computing all its moments. As a result, it is proposed to use an Extended Kalman Filter (Ref. 30) of the following form: (4.24)

= f(x(tlj-l), e, u(t))

(4.25',)

_(JlJ) = x(Jlj-l) + K(j)v(j)

(4.26)

v(j) = y(j) - h(i(J[j-1))

The Kalman gain K(t) is calculated from equations (4.13)-(4.16) by using the time varying matrices H and F, defined by H(t) = _x[ (4.27)

x = lj-l)

F(t) = _xl (4.28)

= : ;:(J IJ)

Notice that the Extended Kalman filter linearizes the equations around the latest best estimate of the state. More advanced filters such as Second Order Filters, Single Stage Smoothing Filters, etc. (Ref. 31) can be used for state estimation, but for the aircraft parameter identification problem, when all the states and accelerations are being measured accurately, an Extended Kalman filter comes quite close to the optimal nonlinear filter in accuracy.

Kailath (Ref. 27) has shown that the density of the innovation u(t) tenas to a Gausslan density as the sampling rate is increased. Thus, for high sampling rates the likelihood function can again be written as N

_0T(J) B-I(J) v(j) + loglB(j) I

J -= log(YNlO) I _!2 I (4.29)

J=l

The validity of the above two assumptions viz: high sampllng rates and accurate measurements should be checked in practice for each application of this method.

Remark: The use of an Extended Kalman filter here is for state estimation only.

It is also possible to use an Extended Kalman Filter for simultaneous state and parameter estimation (Refs. 21, 32). In the authors' opinion, this is not desirable since the uncertainties in the states are much smaller than the uncertainties in the parameters. Therefore, the assumptions of lineariza- tion which are valid for state estimation are generally not valid for para- meter estimation in the aircraft parameter identification problem. Moreover, _he Extended Kalman Filter for simultaneous estimation of the state and the parameters assumes knowledge of the a priori covariances which are unknown for the parameters. This is one of the reasons why an Extended Kalman filter typically gives unreliable confidence limits on the parameter estimates (Ref. 21). The maximum likelihood method described here will be shown to provide realistic estimates of confidence limits on a test case.

It has also been found to converge in several cases where the Extended Kalman Filter for simultaneous state and parameter estimation failed to converge properly due to poor a priori values for the parameters.

4.3 Numerical Optimization Algorithm

The optimization algorithm described here for obtaining the maximum of the likelihood function is the Modified Newton-Raphson or Quasilinearizatlon Method. Determining an update to a set of parameter estimes 8 , which will decrease the value of the likelihood function (cost), J, using this method in- volves computing two matrices: the gradient of the cost with respect to the aJ and the information matrix, M.

unknown parameters, a-_' For the case of a nonlinear system with process noise, the likelihood function, J, is computed using an extended Kalman filter. The gradient and information matrix computation must therefore include at least the first order partials of the Kalman gain with respect to the parameters. With the nonlinear system model given by Equations (4.22)-(4.23), and the extended Kalman filter by Equations (4.24)-(4.28), the gradient of J with respect to 8 is given by n aJ _)T(j) B-I(j) ao - _ a8 (j) v(j) av(j) 1 T(j) B-I(j) aB(J)B-I

ae

J,=l

1 (j)_--_) (4.30)

+ _ tr (B-I av(j)) -H(t) _(jlj-l) ah (4.31) where

ae = as ae

aB(j) = aN(J) p(j [J-l) HT(j) + N(J) aP(j lJ-I)HT(j) and ao k ae k ao k aR

+ H(J)P(jlj-l) --

ae k aek

Based on prior experience (Ref. 33), the convergence of other optimization algorithms, including conjugate gradient and Davldon method, has been found to be slower than that of the quasilinearization method.

3O

The recurslve equations for @P(J l J-l)

are obtained directly from the B0 k Riccati equations, with the separation into update at the measurement times and prediction between measurement times again being made.

Prediction: Prom Equation (4.9), the prediction equation for (j-l) <t <J becomes @P(tlJ-l) = B__F p(tlj_l) + F @P(tlJ-l) + BP(tlJ-I) FT @0k B0 k @0k @O k

_F T @r _q rT @r T

(4.32)

+ _(_lJ _ _k+_ Q_+_ _°_ + _Q_

Update: The update equation, at the jth measurement point, is obtained from Equation (4.161 @P(JlJ) = (I-K(J)R) SP(JlJ-I) _ @k(j) HP(jlJ-1) Be k so k so k @H (4. 331

K (J)_k P(jlj i)

where @P(j lJ-l) HT [HP(J]J-1) H T + R] -I @O k @H T H T + P(j]j-l) _k [HP(j[J-1) + R] -I FB_, HT BP(J 'J-l)H T - K(J) L._P(J lj-l) + H 30 k

(4.34)

+ HP(j l j-l) _k + _k (HP(j IJ-I)HT+R) -1 _H T _R _

The recursive equations for the term a_(J !J-l) defined as the sensi-

a8 tlvity equations and appearing in Equation (4.31), are obtained from the update and prediction filter equations (4.24)-(4.25).

Prediction: a_ af a_ (4.35)

_e-- k = _e-- k + F(t)_k ; j-1<_ t < J

Update: = Be k ae k ae k These same sensitivity equations are used to compute the information matrix, which is given by

avT(J)

av(J)_ T(j) B-I(j) BB(J) B-I(j av ) _,& A a2j ffi E a@k

k (J

.1=1

av(J)

_B(J) B-I(j) __ vT (J)B-I (J) ao_ a@ k (4.37) i I B-I -_tr aB(j) aB(j) I Note that second order partial terms and several first order partials (of matrix inverses) have been neglected. All the partial derivative terms appear- ing in Equation 4.37 can be obtained via Equation 4.31.

The update to the parameter estimates e, or step size Ae, is then com-

puted using the following equation

(4.38)

_/

Since M is the Fisher Information Matrix, M -I provides the Cram_r-Rao Lower Bound on the covariance of the e estimates. The ML method approaches this lower bound asymptotically.

4.4 Relationship to Output Error and Equation Error Methods As stated earlier, one of the principal advantages of the maximum likelihood method is that, under special circumstances, it reduces to the output error or equation error method, both of which have been widely used for extracting stability and control derivatives from flight test data.

For the case where there is no process noise present, i.e., w(t) = O, the process noise covariance, Q(t), is identically zero. With P(0) either equal to zero, (if the initial state estimates are known) or small, this implies P(tlt-l)EO for all t after some initial transient (see Equations (4.9)-(4.16)). The Kalman gain will then also be identically zero (Equa- tion (4.15) and the innovations sequence reduces to

v(t) = y(t) - Hx(t)

for the linear case, and _(t) = y(t) - h(x(t)) in the nonlinear case. In both cases _(t) is exactly the output error.

The only difference then between the maximum likelihood method and the more classical output error method is the choice of the weighting matrix. In the maximum likelihood method, it is given as R, the measurement noise covariance, or its estimate. In the maximum likelihood method, is is chosen as N v(j,_) T(j,8).

J=l For the case where no measurement noise exists, the measurement noise covariance R(t) is identlcally zero. For the case in which all the states and their time derivatives are measured, the likelihood function is the sum of squares of the equation error at sampling times. Thus, the ML estimates are identical to the equation error estimates.

4.5 Identifiability and Uniqueness Problems in Extraction of Stability and Control Derivatives Although the maximum likelihood method discussed in the previous section represents one of the most advanced identification techniques developed to date, there still remain some basic problems associated with extracting stability and control derivatives from flight test data. Most of these pro- blems can be classified under the heading of "identifiability," which is re- lated to the degree of excitation for the particular modes of the system under investigation and the ability to identify the associated parameters. Identi- liability also relates to whether the parameters themselves can be identified or whether they can only be identified as part of a linear combination. This section will discuss some of the symptoms and causes of identifiability pro- blems, and a few of the methods which have been used to solve them.

4.5.1 Symptoms and Causes of Identifiability Problems The most obvious symptoms of identlfiability problems are physically nonrealizable parameter estimates and large associated error covariances.

Either of these symptoms may arise for a number of different reasons. If the input sequence does not adequately excite some of the modes, or if the Stability Augmentation System is operating, thereby suppressing some of the aircraft modes, the associated parameters may not be identifiable. If the model chosen to get the input-output data is inadequate, the parameters of that model may be forced to account for some major unmodeled effects. The estimated parameter values may, therefore, be quite different from what aero- dynamic theory and previous results may indicate. If there are large, un- accounted for instrumentation errors or errors in the location of the e.g.

and the sensors, again non-physlcal parameter values may result. Finally, such additional factors as too short a data length, local minima in the cost functional and poor initial parameter estimates may also result in non- physical parameter values.

Large error covariances principally result from poor input sequences and attempts at identifying too many parameters. The first factor reduces the sensitivity of the ouput to variations in some parameter values, and the second factor causes linear dependencies between parameter estimates.

Since an extraneous parameter in the model does not, by definition, improve or degrade the fit to the observed data, its estimated value will be of no significance and the error covariance of the estimated value will be large.

Probably the most common identlfiability problem encountered in processing flight data results from parameter dependencies. This may occur through a pair of parameters which always appear in the equations of motion together, as with CM and CM. , or through a poor choice of q inputs such that some of the aircraft response variables are linearly correlated, or it may occur through an overspeciflcation of the number of parameters to be identified. In each case, the result of the dependencies is a nearly singular information matrix, which when inverted to obtain the step size in the parameter estimates, causes numerical problems.

Additional numerlcal problems associated with a nearly singular information matrix arise when the control input is expressible as a linear combination of the aircraft response variables. This matrix singularity results from the linear dependence between the partial of a response variable with respect to a parameter in the input matrix amd the partials with respect to a parameter of the dynamics matrix. Since the same singularity exists at each data point, the resulting information matrix will also be singular. A second input related problem arises when the input is of such a nature that the time histories of several of the aircraft response variables appear highly correlated. All the elements of the partial derivative of the output vector with respect to any one of the parameters will be the same, introducing a singularity.

4.5.2 Approaches to Identlfiabillty Problems Four different approaches have been used to alleviate identifiability problems. These are: i. Fixing Parameters - The usual remedy for parameter dependencies has often been to fix some of the dependent parameters during iden- tification. While this generally improves the numerical convergence the choice of a particular parameter to fix and the value at which it is fixed are generally not clear. Although it is possible to fix the parameters at the wind tunnel or theoretical (DATCOM) values, the estimated parameter values may depend upon these fixed values. In those instances where the wind tunnel values are inaccurate or DATCOM doesn't apply and no other a priori information is available, a better way of dealing with the parameter dependency is needed.

2. A priori Weighting - Whenever a priori values exist for certain parameters in a given model structure, they can be included in the maximum likeli- hood method by using a Bayesian formulation. If the a priori values have a Gaussian distribution, a quadratic term involving the weighted dif- ference between the estimated parameter values and the a priori para- meter values is added to the likelihood function. Depending on the weights given to these differences, it is possible to force any of the parameter values to the a priori ones. The a priori values for the aircraft stability and control derivatives are usually derived from the wind tunnel estimates and theoretical calculations.

The weights, which are the most subjective part of this technlquep signify the confidence in the a priori values. An alternate pro- cedure is to successively reduce the weights at each iterationp thereby discounting the dependence on a priori information. The main advantage of this procedure is numerical since the information matrix with a priori weighting is generally better conditioned than the one without it. This procedure is a special case of Tychonov Regularlzation used for solving an ill-conditloned set of equations (Ref. 34).

B Constrained Optimization - If, from practical or theoretical con- siderations, a range of allowable values or relationships between the stability and control derivatives can be specified, they can be used as constraints on the parameter estimates to avoid non-physical estimates. Such a procedure would require a constrained optimization teehnlque in lieu of the Newton-Raphson optimization method normally used (for the output error criterion). Including such parameter value constraints will most probably also reduce the convergence rate.

o Rank-Deflcient Solutions - Without any of the above remedies, the parameter identlflability problems will usually appear as a difficulty with inverting the information matrix and obtaining accurate parameter estimates and error covariances. This numerical problem can be related to the spread in the eigenvalues of the information matrix. A perfect dependency among the parameters should, strictly speaking, result in a zero eigenvalue. However, since round-off and other numerical errors prevent the matrix from being exactly singular, all the eigenvalues will be non-zero with a spread between the smallest and largest elgen- value being many orders of magnitude. In such a case, it might be better to use a rank deficient solution for the inverse rather than a full rank solution (Ref 35). That is, the inverse to the information matrix should be computed leaving out one or more of the smallest eigenvalues. Each eigenvalue which is left out relates to a singular direction in parameter space and, therefore, indicates a combination of parameters which cannot be identified uniquely. (see section 5.3.6) The maximum likelihood identification program described below has options to use the above methods for solving identifiability problems. Further research in this important area is badly needed if identification programs are to be used on a routine basis for extracting stability and control derivatives. It should be mentioned that two other topics related to identiflability (one of which is discussed elsewhere in this report) are those of input design and model structure determination.

4.6 Maximum Likelihood Identification Program This section describes the computer program that was designed to implement the maximum likelihood identification method for extracting stability and control derivatives from flight test data. Three options are provided for dealing with the identifiability problem: (I) a priori weighting, (2) fixing parameters at a priori values, and (3) rank-deflcient solution for the information matrix inverse. At the outset of an identi- fication run one of these three options is indicated (including the weight- ing matrix if a priori weighting is specified) and the program thereafter runs automaticallyz Step size cutting (in the event of a cost increase) and parameter bounding routines are always included in the algorithm, al- though they can be easily rendered inactive, if it is so desired.

The flowchart for the maximum likelihood program is shown in Fig. 4.2.

The principal steps of the algorithm are all blocked out, omitting the numerical procedures used to compute such quantities as the solution to differential equations, matrix inverses, etc. The following paragraphs briefly outline the functions carried in each of the numbered blocks and how the logic progresses from block to block.

READ IN: (I) initial parameter values, Po (NP parameters)

(i)

(2) parameter bounds PL' PU (3) observation and input control sequence Yl, ul (N data pts) (4) initial process and measurement noise covarlance estimates , Q and R, and state error covarlanee.

SPECIFY: (i) which parameters are to be fixed (if any)

(2)

(2) whether a priori weighting is used (if yes, supply weighting matrix).

(3) whether rank deficient solution is to be used.

Compute cost associated with a priori parameter estimates: J = I vT(j)BI(j)u(J) -log [ B(J)

,)

j=l

I

'Solve the following equations simultaneously (using updated parameter estimates) (i) Time history of aircraft states via Kalman or extended Kalman filter (equations 4.24 - 4.26) (2) Kalman gain time history (equations 4.9, 4.15, 4.16, 4.27, and 4.28) ^ (3) _x (3) Time history of sensitivity equations (equations 4.35, 4.36) • (4) Time history of additional partials _P(J[J-I) _Pk ' _ (equations 4.32 - 4.34) _Pk _Pk FIGURE 4.2 FLAXIHUM LIKELIHOOD PROCRAM FLOWCHART Compute

(4)

(1) gradient vector IXI (equation 4.30), (2) information matrix M (equation 4.37) YES rows and columns of the information matrix

'(5)

I Eliminate those rows of the gradient, and those associated with the fixed parameter.

of M (I i is largest, A N is smallest) I ompute eigenvalues h i and eigenvectors V i

(6)

I

lqO I

IS RANK DEFICIENT INVERSE FOR M DESIRED? i i I YES

i

K ffi NP I I SPECIFY K ffi minimum no. of eigenvalues to be retained

t

i

lillgil_l IIIIl'illll II_oIOlill_ Illll_lll IIIOla _ I I

i

i

Compute : t (1) Rank deficient information matrix inverse --M -1 ffi E )'i I- VivtT

!

if1

I

(2) Parameter step size, APt = -M-I DJ

(7) i

Cost associated with new parameter values (3)

I

*" P + APt Pt

i

N

i

I

___ uT(j)B-I(j)u(J) -log [ B(J)

J=_

i

J=1

i

+ (P£ Po )T W(P£ - Po )

i

i

(W -= 0 if a priori weighting NO__T_T used)

i

I

m II, FIGbT_ 4.2 (CONTINUED) 4O

l

(8)

i - NO ! IS Cost, J, greater than cost from previous iteration? 1 YES

! - - 1

Do minimization of J £ times: NC - I, ...

* 1 (9) i - (2) recompute cost =_ IF cost, J, still greater than cost of previous iteration, | increase NC by i.

_|I lilmil_lilililimlil__ mlili (10) From amongst (PN ' APN )' (PN-I' APN-I) .... (P_, APK)

I *

.retain pair with lowest cost; define as P , AP For new parameters p*, are any parameters constalnts violated I YES

(I)

if any component of P (previous parameter estimate set) was on a constraint and AP took it into infeasible region, set that component of _P E 0 i

(2)

compute new P = P + AP

(3)

For each component of P* beyond constraint boundaries,

(ii)

compute ratio DSTI =I API - distance beyond boundary API I if none are beyond boundary, go to step (12)

(4)

find smallest DST i _ DST r replace parameter estimates P with new set of estimates (12) I

I P + AP'DST

i

return to step 3 to begin new iteration

]

FIGURE 4,2 (CONTINUED) Tb_ initialization procedure for the maximum likelihood identification method, indicated in blocks (1) and (2), consists of specifying a set of a priori parameter values, including the measurement and process noise covarlances and, if desired, a set of upper and lower bounds for each parameter. The observations and input control time histories are then read in and stored. Since the maximum likelihood method is a batch pro- cessor, it will use the entire data record for each iteration. The ini- tialization concludes by specifying which of the several options are to be used: (i) fixing parameters, (2) a priori weighting, (3) or rank deficient solution for the information matrlx inverse.

With block (3) the first iteration begins. Using the equations given in Section 4.3 for the nonlinear system equations, the extended Kalman filter, the sensitivity functions and all the required partial derivatives, the time history for each of these quantities at each data point is com- puted. These differential equations can be solved using any one of a number of numerical techniques, e.g., Runge-Kutla. However, the majority of the computer time for each iteration will be consumed in solving these equations. When there is no process noise, Denery (Ref. 36) has shown that by using a transformation some of the sensitivity equations need not be evalu- ated, but rather can be expressed as a function of the others. The number of differential equations _¢hich need to be solved is thereby reduced.

In block (4) the time histories of the quantities computed in block (2) are combined to form the gradient, D J, and the information matrix, M.

Up to this point, all the computations can be performed considering a very complete and general set of parameters to be identified. This set may, however, be more general than needed for a particular application.

For example, it may not be desirable to identify the rudder derivatives if there is no rudder input. If this is the case, the components of the information matrix and gradient due to the parameters which are not to be identified (thereby being considered fixed) must he removed. Follow- ing through the computation of DJ and M, this can easily be done by simply deleting the rows of DJ and' the rows and columns of M associated with the fixed parameters. This is performed in block (5).

Many of the computational problems associated with the Newton-Raphson optimization technique are involved with the large spread in the elgen- values of M. Perhaps the most exact way of computing M -I is therefore to use the eignvalue - eigenvector decomposition. This decomposition is performed in block (6). The eigenvalues and eigenvector are also required if a rank deficient inverse for M is desired.

The second program option consists of specifying if a rank dificient inverse is to be used, and if so, specifying, in addition, the minimum number of eigenvalues which are to be retained in computing M -I . Note that if a rank deficient inverse is not desired, this minimum number is Just set equal to the total number of parameters.

The logic for determining the rank deficient M is given in blocks (7) - (9). For each rank from the minimum to the full rank, the appropriate number of smallest eigenvalues are discarded and the information matrix inverse is computed. The associated parameter step is calculated and the likelihood function value is determined using the new set of parameter estimates. This involves computing the aircraft state and observation time histories and deriving the innovation sequence. Note that the third option enters in block (7) in the specification of whether a nonzero a priori weighting matrix is to be used.

Blocks (8) and (9) are concerned with the solution where the cost determined from the new parameter set is greater than the cost of the previous iteration. (If this is the first iteration, the previous cost is that associated with the a priori parameter estimates). When the new cost is higher, the parameter step size is cut in half and the cost reevaluated. If the new cost is still larger than the cost from the previous iteration, the step size is cut in half again. This same procedure is repeated a given number of times. The reason for this step size cutting is the nonquadratic nature of the likelihood surface.

This same step size cutting routine can be carried out for the rank deficient solution procedure. The end result is (k) sets of parameter estimates, each one resulting in some value of the likelihood function. In block (I0), the set of parameter values and step sizes associated with the lowest cost is picked out and retained. The other sets of parameter values need not be saved.

The final option of the maximum likelihood program is to alter the parameter step size if any of the parameter constraints are violated.

If this option is not desired, the parameter bounds are simply set to very large values. The routine for computing the optimal step size without exceeding the parameter constralntslnvolves four calculations.* The first calculation checks the individual parameter values to see which ones are on a constraint. If the parameter step associated with any of these parameters results in violation of the constraint boundary, that step size is set equal to zero. In the second calculation, the new set of parameter estimates are computed, using the modified step size (some elements are zero).

In the third calculation, each component of the new set of parameter estimates is compared with the constraint boundaries. For any individual parameter value which is beyond the boundary, the absolute value of the ratio of the allowable parameter step to the actual parameter step is computed. This ratio is exactly the factor needed to have that particular parameter value fall on the constraint boundary. The smaller that factor, the farther beyond the constraint boundary the new parameter estimate would have been. In the last calculation, the smallest factor from among those computed for the individual parameter estimates is determined and retained.

* This procedure is based on the Generalized Reduced Gradient method of Abadle (Ref.37).

The final block (12) of the algorithm involves multiplying all the parameter step sizes by the smallest factor determined in block (11) or by I, if the constrained otlmlzation option was not chosen. If the option was used_only one additional parameter estimate will be on the constraint boundary. All other parameter estimates besides the ones with zero step sizes will be within the constraint boundaries.

The computation of a new step size for parameters marks the end of an iteration. To begin another iteration, these parameter values are used in the computations of block (3), and the cycle is restarted.

The original cost now becomes the cost associated with these new parameter estimates.

V RESULTS OF IDENTIFYING AIRCRAFT STABILITY AND CONTROL DERIVATIVES This section discusses in detail, the experience and results of applying the maximum likelihood identification techniaue to simulated and real flizht test data from three different aircraft. Included will be a discussion of the problems that were encountered and all the techniques that were used to alleviate them. Wherever possible, the cause of the problems is also spelled out along with possible implications for flight test procedures.

The first data that was used was from a computer simulation of X-22 VTOL aircraft. The aircraft model was highly non-linear and the data included process noise as well as measurement noise. Experiments were run with different input sequences and different measurement noise levels to investigate their effects on the parameter estimates. In all, 23 parameters were identified, excluding the measurement and process noise covariances, which were assumed known.

The second case involved actual flight data from an HL-10 lifting body. The digitized data, comprising approximately six and one-half seconds of flight, was supplied to SCI by NASA-Edwards FRC. A linear aircraft model was assumed in fitting the data and, in all, 20 parameters were identified, including the measurement noise covariance and the initial flight conditions. The data was assumed not to contain any wind gust (process noise) effects.

The third set of data, also supplied by Edwards FRC, was from an M2/F3 lifting body. This data, covering approximately eight seconds of the flight test, did contain wind gust effects and represented a test of the maximum likelihood technique in reducing flight data which had not been successfully reduced by the output error technique. In all 22 parameters were identified, including the measurement noise covariance, the parameters (time constant and driving noise covarlance) of a wind gust model and the initial flight conditions.

5.1 X-22 Simulated Data At the time this contract began, the maximum likelihood method had been applied to simulated X-22 data, containing the effects of gusts, with very promising results. However there remained several important problem areas which needed further investigation and improvements to be made to the existing program. This section will outline these problem areas, including the method of approach, results, and conclusions.

5.1.1 Generation of X-22 Simulated Data The model of the longitudinal motion of the X-22 is given in Appendix A. These equations can be put in the nonlinear form

x_'= f (x,c,p)

+ g (x, p) v

where x = [q, @, u, w] T is the 4-dimensional state vector (q_ pitch rate, 6_ A pitch angle, u-A longitudinal velocity, w _ vertical velocity) p is the 23 x i vector of unknown parameters (consisting of the coefficients of the polynominal expansion in u of the derivatives Mo, Mw, Mq, M_, Xo, Xw, X_, Zo, Zw, Z_) c is the vector of deterministic control surface deflections and biases = [i, _es IT (_es _ elevator deflection) v is a 3-dimensional white, Gausslan process noise with mean O and covariance Q.

The elements of the g matrix are obtained from the matrix of 3f first partials _x and, therefore, the parameters and states appearing in f also appear in g.

The measurement equations are X x

."I

n 5 q + Z 6 + qw + g sin O n 6

nx I

- qu - g cos 0 n 7 m . i m where n i , i = i,. • •,7 are independent, white, gaussian measurement noise samples with the properties E{n_t} = 0 and E{_I = R6ts Substituting for &, _, and 5, however, introduces, process noise into the measurement equations. The measurement equation can then be rewritten as _x- +n+ Z _ I f' where f' and g' are made up of specific rows of f and g, respectively• This gives rise to a correlation between the process noise and the measurement noise now consisting of the sum of the vector _ and _ X.

With the specification of an elevator deflection sequence, _es' and

the process and measurement noise covarlances, Q and R, the data, z, could be generated (using 4th-order Runga-Kutta integration of the nonlinear equations of motion). For each trial approximately i0 secs. of data was used, with a sampling rate of 20 per sec.

5.1.2 Program Description The program that was initially used to extract the stability and control derivatives from the simulated data consisted of basically two parts. The first part was a least squares start-up routine (Ref. 13) which generated an initial estimate of the parameter values. This least squares technique is an equation error method which, in one pass through the data, obtains parameters estimates that minimize the following criterion rain (x--i - fi(x' P)) 2 i=i,..,4 -- , for each derivative xi which is measured or derived. Since, in the X-22 simulation, _, n and n were measured, it was first necessary to x y express these quantities as linear functions of the parameters to be identi- fied. From Appendix A it is possible to write _ as a linear function of the parameters (polynomial coefficients) in the derivatives Mo(U), Mw(U) , M (u) and M_ (u); n as a linear function of the parameters in X (u), q es x o Xw(U) and X 6 (u); and n 6 as a linear function of the parameters in Z (u), es o Zw(U) and Z6 (u). Since no parameter appears in more than one expression, e8 a unique least squares estimate can be obtained for all the parameters.

The second part of the program was the maximum likelihood identification technique in a form designed to identify the parameters of a non-linear model when both the process and measurement noise co- variances, Q and R, are known. The least squares parameter estimates were used as the initial conditions for the first iteration through the data of the maximum likelihood routine. These estimates are updated with each itera- tion, untll the algorlthm converges. However, since one iteration through the data of the maximum llkellhood technique required I minute of UNIVAC 1108 CPU time, only a few iterations were used.

5.1.3 Limitations of Previous Results The first limited trial of the maximum likelihood identification technique applied to the problem of extracting aircraft stability and control derivatives was on simulated X-22 VTOL data supplied to SCI by Cornell Aeronautical Laboratory. The complete data set consisted of four cases; two without process noise and two with process noise. In each case a single step input was used to generate the data. For the no process noise cases 2A and 2C, satisfactory estimates were obtained for all parameters except for X 6 and Z 6 derivatives. For both the low process noise (2B) and moderate process noise (2D) cases, the errors in all the identified parameters were much larger.

However, when a multistep input sequence, which supplied much more excitation, was used, the results for the process noise case improved greatly. It soon became apparent that the quality of the parameter estimates were very "input" dependent.

As outlined in Section4.3 , the calculation of the update in the parameter estimates involves the computation of the gradient ___JJ, and the information matrix _p _2j, where J is the likelihood function (Eq.4.29). These quantities in turn in- volve solving a differential equation for the sensitivity matrix _ (i/i-l), _P where _ (i/i-l) is the output of a Kalman filter. Since the equation for the state estimate error covariance, P (i/i-l), does not reach steady state and involves the unknown parameters, the partial derivatives of P (i/i-l), and the Kalman gain W i, with respect to p should be included in the computation of _J and _2j . (See Appendix B) These were neglected in the earlier X-22 _p _p2 5O identification work. It was orglnally thought that the lack of monotonic convergence of the identification algorithm could be attributed to these extra partials being neglected in the computation of the partial derivatives.

As noted earlier the accelerometer measurements introduce process noise into the measurement equation, thereby correlating the total effective measurement noise with the process noise. This correlation which effects the equations of the Kalman filter and, therefore, the computation of the sensitivity matrix, was not accounted for in the earlier application. There is an additional correlation between the g(.) function in the dynamics equation and the g'(.) function in the measurement equation since both are a function ^ Of the state estimate, x (i/i-l). It was originally thought that this might also have a significant effect on the parameter estimates and on the standard deviations supplied by the Cramer-Rao lower bound.

In both the data supplied by Cornell and generated at SCI the process noise and control were kept constant over an integration step. It was important to distinguish the cases where the process noise changed value before or after a measurement. In one case there would be a correlation between the measurement noise at a sampling point and the process noise during the pre- ceeding integration step while in the other case the ocrrelation would be with the succeeding integration step. This difference, though subtle, is important.

All these areas were investigated with the objective of determining the effects, on the parameter estimates and standard deviations, of different modifications. The following description of the results is broken up into separate sections, each one involving a different area of investigation.

Included in an accompanying table are the parameter estimates and standard deviations resulting from each change in the algorithm. These standard de- viations are actually lower bounds on the actual values and are obtained from the diagonal elements of the matrix/_2j_ -I. Also noted, in each case, is the number of iterations of the algorithm used in obtaining the results.

5.1.4 Comparison of Results with Sinsle Step and Multi-Step Input Sequences The initial processing of the data supplied to SCI by Cornel Aero- nautical Labs resulted in unsatisfactory parameter estimates both for the low and moderate process noise cases, as shown in Columns i and 2 of Table 5.1. Since it was already known that the input sequence shown in Fig. 5.1, used to generate the Cornell data did not sufficiently +1.0 TIME (SEC)

t I I '1 I I I I I>

2 3 4 S 6 7 8 9 10 0.0 _ ![RIM -1.0 FIGURE 5.1 INPUT SEQUENCE USED IN GENERATING COP, NELL DATA excite all the modes of the system to allow adequate identification and since there was considerable uncertainty already as to how this data was generated, SCI programmed its own data generator using the equations of motion in Appendix A. The elevator deflection sequence used by SCI in generating the X-22 simulated flight data is shown in Fig. 5.2 and the process and measurement noise standard deviations are given in Table 5.2.

TABLE 5.1 X-22 ZDEHTTFICATION RESULTS IParaameterl I | Case 20 2: Case 2B 3 : Case 2B 4 : Case 23 5 : Case 20 6 _ Case 211 After After ull in8 SCI WZth With With In|flaX |nlt fal J'r,ul ¢ |eta p Forvord Coaltlnt 8 Constant • COet J Ch_n_ee Ch.ln es Correlation and Cot relat|oa 3,3L27 3.)7_$ 3.5388 3. $/,99 3.5.355 3 _17 31.47 )2.02 -. 176 -1.632 -°835 -. 159 -.206 -.209 -.224 -.000355 .00429 .O02U -.O004.6_ -.000284 -.G00292 +. 000174 -.148 -.166 My(: ) -.1 -1.268 -.671 -. 09t, 7 -.137 -.00317 .00558 .000009 -.00321 -.0027& -.00268 -.00235 -.SO0 "-.436 Hq(1) -.497 -2.0721 -.99$ -.,89 -.S20 u -.00103 .00434 .00278 -.00114 -.000814 -.000975 -.00182 18.15 18.11 (:) 18.66 -11.85 ,,,5 18.38 18.17 K6 .0669 • 311 • 179 .0708 • 0734 .0733 .0758 17.87 17.97 X u -.0917 -.352 -.160 -.103 -.110 o 0.0869 0.0874 / 12/ 18.3 3.5.80 23.01 _ 18.88 19.41 v -.G003 .000659 -.0000539 -.00024 -.000228 -.000`3: -.000.316 .2233 -.223 (_) .2211 .0162 .176 .220 ,2085 rv -.00159 -.0000_1 -.0OJ22 -.00159 -.00143 - • 0015 - • 00160 -,896 -,875 (:) -.778 8.$86 -.976 -.691 -1.011 x4 .0164 -.0521 .0202 .0171 .0212 • 0199 .0196 /1_-32.17 '5.95 -31.22 -35.15 -2,.08 -30.19 -28.0.3 ,_u ,) .91 -.1367 .511 .969 .825 • 879 .813 _tu2/ -.007 -.0017 .00568 -.00728 -. 00669 -.O0692 -.0066 -..3.397 -..356 (_) -.2939 -1..320 -,$99 -.272 -.361 Iv -.00287 .00516 --000399 -.00314 -.00213 -.00233 -,00213 -1.136 -1.11 (:) -.351 20.6 -11.26 -.67`3 -1.143 s4 -.0167 -.122 .106 .0200 .O269 • 02697 +. 0271 Mo. of' lteratLons .1 2 2 2 2 2 Cost ji .81S8 .9061 .$296 Standard L Dev/at Ion8 2.M) 3.$01 1.104 .254 .0352 ,0371 • 0418 .0573 • O0/+O$ .000996 • 0OO145 .0001_4 • 000171 .00023$ • 0299 .0399

_(1) ,9, .8 +2,1 +2,9

00,,, 00_. +., _0., • 000279 .000371 • 0164 .0157 ,,(_) -, .., 0., .0.

00.371 001.3, 000,,2 0.1.

• 000157 .00017`3 • `3359 ..3368 ,,(:) 1_.6 18.. .279 .2.

.,,2 .1+ .,_6 .00,1, • 00361 .00356 • 7890 1,010 XO .276 .0688 .O1.32 .0135

(i,) ++ "+ ++ ++ • 0127 .0163

.0009_;0 .000241 .O000S .0000_,31 .0000508 .0000660 .00919 .0116 X(_) .17`3 .0422 .0105 .OlO] • 00142 .000.151 .0000899.0000893 • 0000857 ,000109 • 1073 .224 s., (:) 41.0, '.686 .11)9 .11212 .3173 .0596 .001134 .0011, • 00U6 .00134 &.667 6.533 so+ 1.276 .283 .0427 .04486 • 0767 .I057 • 00,28 .000972 .000172 .000181 • 000299 ,000428 .0_45 .0749 'v(_) .779 .16, .032`3 .0334 • 00657 .001,2 .000321 .00q330 • 000507 .000696 • 636 ,692 • 4(1 t 271.0 '1,10 .`388 .398 u 1.21 • )189 .004`31 .0046 • 00688 .00737 TABLE 5.1 CONTII_UED 8 z Cane 2D ! 9 s CaN 2D _1. Cne 2I !

Vich With With Constant I Constant $ Some Added No Added lud Add _JctlIXe rllCtllll IPlct tll 3.$368A 3.6949 $.$76Y 47,291 76,69 _3.58 -.427 -1.049 -.261 +.000295 .00345 -.000035g -.421 -.583 -.157 .0000073 .00193 °.00247 -.5068 -.321 -.502 -.000826 -.00345 -.00102 15.92 17.79 18.04 .1003 .0807 .0752 16. _g 13.$0 19,17 -.0_,93 -.0319 " -.)09 - .000449 -.000449 -.000221 .228 .2764 .215 -.00159 -.00191 -.00150 -1.402 -.£61 -°$$2 .0261 .0181 .0207 -14.18 2S.5i_, -|5.88 .536 -.179 .196 -. O061g -.0024$ -o006$6 -.634 -og65g -.:M_9 .000741 +.00389 -.00196 -4o618 -1.146 -1,213 .0703 .0278 .0289 2 2 :t 22.893& 23.167 .5_29 14.06 .2289 .000534 .1(;28 ,001$1 ,0760 .000824 1.706 .0181 4.270 .0688 .000217 .0498 .000462 .S4cl ..

.00581 2$.27 .406 .00163L ,2e#51 ,00273 3.231 .03443 | _r. (sin) 0,0 I, ', I I I 3 4, _"

I

FIGURE 5.2 MULTISTEP INPUT TABLE 5.2 STANDARD DEVIATION OF PROCESS AND MEASUREMENT NOISE Standard Deviation Low Moderate gusts (process noise) 1.0 fps 5.0 fps 1.0 fps 5.0 fps .2 deg/sec 1.0 deg/sec measurement noise n 0.5 fps 2.5 fps u n 0.075 fps .375 fps w .03 deg .15 deg n O n .01 deg/see .05 deg/sec q n ,001 g .005 g n x n 005 g .025 g n z no D025 deg/sec 2 .0125 deg/see 2 q The initial state estimates used in the Kalman filter were set equal to the state measurements observed at the first data point, i.e., _i (I/i) = z i. Since these measurements consist of the true state _i(1) plus noise, the error covariance P (i/i) is given directly as P (i/i) = E {(x (I) = __ (i/i)) (x (I) - __ (i/i)) T }= dlag{ RII, R22, } where Rii is the i th diagonal element of the measurement noise R33, R44 covarlance matrix, R.

The effect of using the SCI data generation program with a multistep input sequence instead of a single step input sequence can be seen by comparing Cols. 2 and 3 of Table 5.1. The parameter estimates are greatly improved and the standard deviations are reduced. This enhanced ability to identify the parameters is attributed to the fact that the more varied the input sequence, the more the system modes are excited and the higher is the signal-to-noise ratio at the output.

It is important to realize that the results as shown were obtained for only one noise sequence, and therefore, although the parameter esti- mates improved considerably, the parameter estimates by themselves are not sufficient for comparison. Neither are the costs, themselves, since changes in the noise sequence will influence the costs.

As will be seen in the later areas of investigation, some of the para- meter estimates improved as the result of some change in the algorithm and some did not. This is almost always the case, and unless the individual relative effect of the parameters on the cost is known, it is very difficult to say on the basis of only a few of the parameter estimates improving, that the algorithm itself is improved by any change. The criterion that is more suitable for comparison is the standard deviation of the parameter estimates.

A lower bound on the standard deviations is obtained from the inverse of the in- formation matrix and this is adequate in many cases. However, if the differences in standard deviations are small, the Cramer-Rao lower bound may not reflect these differences.

5.1.5 Comparison of Forward and Backward Correlation The next area of investigation involved the effects on the parameter estimates and standard deviations of the type of correlation between the input and output noise sequences. The original Cornell data specified that the process noise was kept constant over an integration step. This meant that whatever correlation existed at the measurement times (note that the measurements are taken at discrete instants) also existed throughout the entire integration interval. The SCI multlstep data was first generated with the accelerations 4, n x and n being calculated using the process noise from Y the previous integration step and the control, which was also held constant over an integration step, from the next integration step. Figure 5.3 below graphically shows when the values of vl, the process noise, and 61, the control, were changed in relation to the measurement instances.

v i vi+2 _ Process Noise 6i+i 6i+2 -_-- Control

/

Measurements FIGURE 5.3 v AND _ ORIGINALLY USED IN THE CALCULATION OF q, n , n (DOTS INDICATE THE INTEGRATION TIME POINTS) x y

This seemedinconsistent with the way correlation betweenprocess and

measurement noise is usually represented in state space models. For example,

the discrete analog of the continuous time representation

_(t) = Fx(t) + G v(t)

y(t) = }Ix(t) + n(t)

is where _,F and G,P are related by T, the sampling interval. In this model, v k and n k are correlated, i.e. the process noise during (tk, tk+l) is correlated with the measurement noise at t k. The correlation between n k and v k effects Xk+l, not x k. Similarly, Xk+ I is calculated using 6k (the control at time tk). In continuous time this means that the correlation between v(t) and n k must be during the interval between t k and tk+ I. In addition, the derivative _(t) at time t k must be calculated using the control that existed between tk_ I and t k. Therefore, _, n and n should be calculated using the x Y process noise that will exist in the next integration step and the control that existed in the previous integration step, as in Fig. 5.4.

vi+2 • :_=_- Process Noise i Vi+l 61+i _--_ Control AFTER CHANGE n FIGURE 5.4 v and _ USED IN CALCULATING OF q, n x, Y Note that the correlation between measurement and process noise has changed from a backward correlation to a forward correlation, i.e. the measurement noise at the sampling instant is now correlated with the process noise during the succeeding integrationinterval_ This forward correlation, if unmodeled with effect the (forward) Kalman filter operation while the backward correlation will effect a second (backwards), smoothing run through the data. In an actual application with continuous time dynamics and measurements, this problem will not arise. However, if discrete measurements are recorded, the type of correlation that exists will be an issue, and for data generated by a physical system the forward correlation is the more natural.

The effect of using the forward correlation in the data generation and then identifying the parameters, although not accounting for this forward correlation in the Kalman filter, is shown in Cols. 3 and 4 of Table 5.1 (for the low process noise case, 2B). Almost all the parameter estimates have degraded, offset by an accompanying slight increase in standard deviation.

These results indicate that if forward correlation (the type normally used in computer models of discrete systems) is not modeled, it can have a detrimental effect on the quality of the parameter estimates.

5.1.6 Additional Performance Index In Col. 4 of Table 5.1 a new performance measure is introduced, labelled J'. This represents the unweighted mean square error in estimating the output, given by N j, _ 1 ^ 2-N i_l[Zi - h(xi/i-l' _)]T [z i - h(_i/i_l ' _)] J' was introduced as another means of comparing the results of the different runs and differs from J by the fact that changes in the Kalman filter only effect it through the state estimates xi/i-l" It was not substituted for J in the identification algorithm since it does not weight

the various state estimates and is not the likelihood function. As was

pointed out earlier, the weighting matrix used in J was a function of

the Kalman filter covariance and therefore varied as different changes

or additions were madein the filter.

Experience with the algorithm and the two costs, J and J', indicated

that J' was more sensitive to the parameter estimates than J, which

always appeared to be in the range 3.5 to 4. However, this is to be

expected since, if the parameter estimates are bad, the state estimates

will likewise be bad and the associated state estimate error covariances

will be large. The weighted residuals, which are inversely proportional

to the state estimate error covariance, may therefore changevery little.

J', on the other hand, has no weighting, and therefore reflects the

absolute accuracy of the state estimates.

5.1.7 Accounting for Correlation Between Process and Measurement Noise

The initial attempts to account for the correlation between the measurement and process noise were not successful either in reducing the cost J or in improving the parameter estimates. It was decided that part of this problem was due to the fact that the noise term g(xi).v i appearing both in the system equations and the measurement equations, depends on the state x i and gives a long-term correlation. For this reason, it was decided to alter the system equations to include a constant g matrix, calculated from the initial control values, the nominal state values and the actual parameter values.

The identification program was first run without accounting for the process noise-measurement noise correlation, in order to get a new standard for comparison. As shown in Column 5 of Table 5.1, many of the parameter estimates were worsened and the J' cost increased slightly. The fact that the J cost decreased slightly can be attributed to the fact that the weighting matrix is a function of the g matrix.

Also the worsened parameter estimates were not,in all cases, offset by increased standard deviations.

The improved accuracy of the parameter estimates with the non-

constant g matrix can be attributed to the fact that the parameters

of the g matrix account for the system gain factors and are therefore

easily identifiable. Since, in the original system model, the same

parameters appeared in the f and g matrices, the parameter estimates

were overall improved.

It was agreed that, although the constant g assumption was a

large change from the original problem, it did not represent a

departure from reality. As can be seen from the system equations, the

g matrix was originally constructed from the linearized f (',')

matrix, the motivation being that the process noise would then

enter the dynamical equation linearly.

With the constant g assumption, the correlation between process and

measurement noise was accounted for by adding the following terms (indicated

by |_| !I) to the identification algorithm. (See Section 4.3)

T

Defining S i = E{v. n.} to be the measurement noise and process noise

_th m I

correlation at the m measurement time: Kalman filter state prediction:

- I R-I I

x = f(_) +I g Si_ 1 [zi-i - h(xi-i/i_l)] I

L _!

Kalman filter covariance prediction: = T gT Pi/i-i _i-i Pi-I/k-i _i-i + g Q (AT)2 R-I T I - g Si_ 1 Si_ 1 gT(AT) 2 I

L l

"- 7 where 0i_ I = I +_x _AT- I g Si-i _P + 1 I L x"xi- i / i-i x=xi-i / i-I Sensitivy equations:

I

ti_ I < t < t i pj' Vp+Tp - ,gsi-i Tx +Tp , I __I

x--&i_i/i_ 1

The resulting parameter estimates and cost are given in Col. ,6 of Table 5.1. Comparing these results with those of Col. 5 it is seen that, while some parameter estimates improved, others did not. The overall cost J remained the same, while J' decreased slightly. An important point is that the standard deviations for almost all the parameters increased.

This implies that, with the initial set of least squares parameter estimates as good as they are, the inclusion of the terms accounting for the input/output noise correlation does not gain much by way of the parameter estimates.

However, with these correlation terms included, the standard deviations come out to be more realistic, in view of the differences between the actual and estimated parameter values. All this is not to say, however, that for a less exact set of initial parameter estimates, the correlation terms won't improve the algorithm performance.

5.1.8 Inclusion of Additional Partial Derivatives There were two principle motives for adding the additional first order partial derivative telnns of the state covariance matrix to the identification algorithm. The first was that they would be required if Q and R were to be identified, since they both appear explicitly in the equations for the state estimate error covariance. The second motive was that the cost J, instead of monotonically decreasing with each iteration, was oscillating.

A possible cause of this was that the gradient direction was being

calculated incorrectly due to the fact that the neglected partial

derivatives have an appreciable effect on the gradient of the likelihood

function in the vicinity of the minimizing set of parameter estimates.*

It was necessary, therefore, to investigate the importance of these

extra partial terms in identifying the system paramters.

The inclusion of the additional first order partial derivatives into

the identification algorithm presented special problems due to the nature

of the X-22 model. Since the g matrix is a function of both the states

and of the parameters, the derivative of the g Q gT term (appearing in

the covariance equations) with respect to p must be included. This is a

particularly lengthy computation. A first attempt to include all the

additional partials except those of the g Q gT term is shownin Col. 7 of

Table 5.1. Comparingthis with Col. 4 (since constant g was not assumed), it

is seen that both costs J and J' increased slightly and someof the parameter

estimates, themselves, are slightly degraded.

Twochangeswere then decided upon. The first was that case 2D

(moderate process noise) would be used instead of 2B (low process noise).

The second was that the g would again be modeledas constant.

The first changewas motivated by the desire to see larger variations

in the costs. More process noise would makethe initial least squares

parameter estimates worse and therefore the effect of the additional

partials would, potentially, be the greatest. The second changewas

motivated by the fact that with the constant g assumption, partials

of the g Q g terms are identically zero.

A second possible cause was too large a step size, which was corrected

by halving the step size along the calculated gradient direction whenever

the cost J, increased.

Column 8 of Table 5.1 gives the parameter estimates, cost and standard deviations for the 2D case without the added partial terms, and with the constant g assumptions. Column 9 gives the same results with the additional partials included in the algorithm. Once again the cost, J, increased slightly ( _ 5%). This can be attributed to the fact that the convergence of the algorithm with the added partial terms may be slower.

The standard deviations of these parameter estimates (for both Cols. 7 and 9 ) are not given since their calculation requires another full iteration of the algorithm. Since the added partials quadruple the run time per iteration, it was decided not to calculate these values. The important point is that, considering the slight variation in J, even for this worst case, and the possible benefits of the added partials in terms of the vastly increased run time, it is not necessary to include these added partials in the identification algorithm.

5.1.9 Aerodynamic Derivative Estimates As was noted at the beginning of this section, the aerodynamic deri- vatives themselves were not identified. Rather, the coeffients of first or second order polynomial expansions in the longitudinal velocity, u, of these derivatives were identified. Using these identified coefficients, it was then possible to reconstruct the time histories of the total derivatives, and compare these estimates with the actual values. These comparisons are shown in Fig. 5.5 for the 2B data and the original model structure. The fits to most of the derivatives was good. This indicated that, although some of the estimates of the polynomial coefficients had relatively large uncer- tainties,their influence on determining the total derivative behavior was small.

X22 STrlOILITY FIND CONTROL DERIVrtTIVES X2_ STABILITY _ c'ONTROL DERTYRTIYE$ II ! • i • i - i • i • i I " i • i - . i - i • -.Ill • I - i • i • i • I " ' i • i • i - i " i - ",Ill o,III _.lYI p.

MII III{_ iii III ._L,Itl k- l.- I I l-,IO0 o I_L,lol r 31: 31131 II KIII Ill I Ill Ill II IK l! + _xx x -.ill -.Ill 't'llt,l-t,tll,l,111, • _.III 1.0 I.I l.l I.I +.0 S.0 1.0 Y.I 1.0 1.0 I0,0 11.0 l.+ TIME - SECOndS I.O I,I I.I t,I I.I • l.I ?.I I.I lol II.t II.I I I ! f ! I t l t l ; TIME - SECONDS C_ U7 114 $I _ N 16 _ M i) )+ I I f 1 I I I 1 _f I J ul u1 iii II . 12 Ii i_ ii II 4l IP letl_tlllllX_ I'll_ ?r _2 _maILttY FIND CONTRCL DERtVRTtVE$ -,24 • X'21 STRBILITy RND CONTROL DERIYRTIVES • i - i • i - i • i i • i • i • i • i • • I • i • ! i • i • I • i • i • i • i • ; I_ M i{ -.IT R x -0|i tIP II I( + Ill M |I ++.

p. |I x • J¢ i( x x iI ii M ii l( x I _C ql II 1" M " .

II -.I!

IJ l.l • I,I + 1.0 $.0 ?.0 l.l l,l 10.O ]l.I IIlc_ • W , ! , t , ! , t , ! ,. t _ ! _ _ I . I , _ T|HE - SECONDS ! °0 t.B "I,l I .I °°0 ¢ .I T.0 |.I |.I II.l | l.l I , I l I I " I I I I I TIP_ - _COmOS II_ III H Irl II M 10 M I_ St " . I . I I I l I I I l I .i Ul Illi II m II II . II 14 II II I{i_I_IIIM. IIIIJ_IIT FIGURE 5.5 X-22 ESTIMATED AND ACTUAL (SIMULATED) STABILITY AND CONTROL DERIVATIVE TIME HISTORIES X22 b"TRDZLZTY P_O CONTROl. DERIVRTIV1E$ X22 STF_L_ITY _ CON'TROt DERTVRTTVEg Jtp • • i - | • i i . i - i - i • nu • ! - | " u " tll • I " ! " I _ • I 1 • ' i - ! " i - e " .ll .H .It • X illll IM HIll • _ .It .el .80 _¢ °e4 • .Or • ! . f • I • I m * • I • ! ' f • | " f * O.i)e • ! . ! * l * I * 1 • f * e • I |.qP N.O I,II ,,O 0.0 I°0 l.o IB.O q, .Ill $.ql $,g 1'o0 II°O I.O I*1 l.O l.O I.I 5.il I.I ToO I*(J I.l |O*l IJ*l TtNE - SECONOS Tit'S. - So.COt',lOS f I I I ,, I l I " I f I" I I , ! I ! I ! I ! ; ! '" u I11 I1_ N _ N II 8o 64 -- IL,It 114, H iS N It* ii Ikb i:l II lamaCllmDmlJ. _tLec¢lrlr O_ X22 STNBILITY _ CONTROL. DERTVflTIVE$ Cl Ill Iii ]11 ..

l .:It |.0 S[ W[ oil .8 *4 . ! . I , I1 , | • I • 1 • I , I , t ."He .

.t I.I JJ _,1 • |.1 |.O 5.11J I.O T.O l.§ J.O 10.0 J|.al rlr,_ - SECONDS I I l I t l I 'l f t J i11 ul fl 1o 141 16 to 1,4 II lip MI_II_WlL VlUCI_ FIGURE 5.5 (CONT.) X-22 ESTIMATED AND ACTUAL (SIMULATED) STABILITY AND CONTROL DERIVATIVE TIME HISTORIES X22 $TADIL|TY _D CONTR_ OERZVNT|VES _2 STRBItTI"Y _ CONTROl. DERIWtTTV_9 -.au o • | - ; • l • i - ! - 1 '. ! • i • i - | - X_XXXXXXXX • i - i • | - u • i _ ! - n - u - | - i -

!

!-i

-,| !

_-°$q *Z0 I -.St I r "N NN • *.¢1 -,¢$ MX _ll;t * I . ! P • I t * | . t • , . a 1 f * I a _$ tfl *-t.t-e.l-l.f-w-9- 0.0 I.O I,g I.e 4.0 1.0 t.O T,0 IJ 1.0 10.0 3l',O 8,0 1,0 t*l I.I t°O I.| 1.0 ?.1 |,O $°0 ||,0 I1°1 TJ_ - SECO_S • • TIP_ - $[C0_S I,I f I I I I 1 f I "!

I l l "'f I I I I f I I u; st4 go _ _ e! ,_ _ 63 $1 lit' • 114 tl t_ _ 16 IO 14 I$ JJ| ILIII_I 1 lll41b lllLdbCl'Irff ,,,..,I X22 _OILI_ f_DC_ DERIVNT[YES II°t _xkKl - i - z - I " l • l • l - J • i - l - t.| |,e . t.| a.t i '$.t 9.1} N .I .$ lllllll ,I.O.l,l',llflt_f,?.,,_.

.4.

|.0 f,0 t,t l.l rl.O" l.O l,O f,I I.I l.O ll.l llJ Tl_ -_CGNOS I f f |' t f I I f I I Iz; 114 tl _ u lo _ 13 J_ FIGURE 5.5 (CONT.) X-22 ESTIMATED AND ACTUAL (SIMULATED) STABILITY AND CONTROL DERIVATIVE TIME HISTORIES 5.2 HL-10 Flight Test Data This data was used mainly for checking and validating the maximum likelihood program. Flight data for the HL-10 lifting body was supplied to SCI by NASA FRC (Flight Research Center) Edwards along with as much information about the flight condition as was available at the time. FRC also supplied SCI the results of their stability and control derivative extraction program along with the specific measurement and a priori parameter weights that were used. The HL-10 data did not contain gust effects and a linear model for the lateral dynamics was used.

It should be noted for the case of unknown measurement and process noise covariances an additional term, N_n(det(HP Ht+ R)) is R_ added to the cost criterion, where Pss is the steady-state Kalman filter error covariance matrix. This is because the weighted mean square error N will always have a constant value of _ x (no. of states) since the weighting matrix is the sample covariance.

5.2.1 Dynamical Equations of Motion and Observation Equations The linearized lateral equations of motion, including the effect of the wind gusts, are: -- u

0 " rp"

Lp L r L¢

P 0 rr

i

Np Nr NB

cose .B

C ) =

Yp+sin a Yr -c°s a YB V

1 tan e 0

F

m

L6 Lo _a Le

L6 a

r

N6 NO 6r

N B

N6a

+ r

(S,:l)

B n

Ye

Y6 Yo 1

Y6a

r

0 0 0

G u

r where p is roll rate (°Isec) r is yaw rate ('Isec) is sideslip angie (°) is roll angle (°) 6 is aileron deflection (') a _r is rudder deflection (o) _n is wind gust in equivalent sideslip angle (o) C is a transformation matrix

Ixz

1 -IT o o

. Ixz

I 0 0

I z

C

0 0 1 0

0 0 0 1

All the quantities are in the body axes system.

For both the HL-10 and the M2/F3 flight data, (which is discussed in Section 5.3), a , the angle of attack, y, the flight path angle, and V, the velocity, were assumed constant over the data record and their values were supplied The same was true of I I and I . In • XZ ' X Z most cases, Y and Y were assumed to be zero, which implied that there p r were nominally 20 parameters to identify (7 in F, 9 in G and 4 initial conditions), excluding the wind gusts and biases in measurements.

The observation equations, again assuming the existence of gusts, are given below: m m

0 0

l 0 0 0 0

P

0 6a

Yl

0 0

r

0 l 0 0 0

0 _r

Y2

0 0 0

0 0 1 0 0

+

l

Y3

0 0

0 0 0 l 0

L-C-

Y4

Y6a Y6 r Yo

0 0 YB 0 0

LYs_

_ - __

H D u

Y

0 nl

0 n2

-- -]

+ l Bn + n3 (5.2)

,!0 n4J

_YB _ n5

F n where Y5 is the lateral acceleration ni, i = i,...,5 are independent white noise Gaussian measurement errors.

Notice that Bn always appears with B when gusts are assumed to exist.

In the processing of the HL-10 data, 8n=0, thereby leaving the measurement noise, _, as the only random noise source. The Kalman filter equations reduce to the original state equations and the Maximum Likelihood Method is essentially similar to the Output Error Method With the exception of the weighting matrix.

5..2.2 Characteristics of HL-10, Fli_ht 19-2 The HL-10 data was assumed to contain no gusts and was, therefore, processed by using the generalized output error criterion. The results, however, were somewhat unexpected and, in all, eight different runs, each with a variation on the original output error run, were made to solve the problems which were encountered. It was later learned that these same problems had been encountered by FRC. Each of the different runs are described in this section, along with the objectives, observations and conclusions particular to each.

Along with the observations and control time histories (see Fig. 5.6) wind tunnel derived parameter estimates were supplied to SCI by FRC. There were 327 data points at a sampling rate of 50 per second, for a total of 6.54 seconds of data. The angle of attack of this flight was 16.8 ° and the mach number was 1.22.

It was evident from the data that a substantial amount of clipping and quantization had occurred during data collection. An accurate model for the observations would include the dynamics of the instrumentation system, but no information was available. One effect of not including a model of the instrumentation dynamics and quantization effects would be to have correlated nongaussian residuals, in each of the observations.

_L ROLL RnTE(DEG/$EC) ROLL(OE6) I .-t I-4 I t_ t11 I I YR_ RRI'E ([_G/SEC) SIDESt._P(D£G) o_ 0' ., L._ ! i D | i ! !

LATER_. RCCELERRTION-TrHE .11 • | , i • | - i • i • f • i - ! - i • ! • i • ! - i - ! - i _ .St .I|

• l

.1¢ • -°Or

V

_J *.0t _ -.01 • _ -,11| -* " J f i ! * | • | i _ i | , _ • | • ! • | , | , 1 I | |,O .q ,l I,t |.| ' It,l 1*I Z.| 3.1 11.| 1.| 1.4 1.8 11,t $.$ I.# $,1 TIHE(SEC) R*J'J2_J_ POSlTZ O_-TIH_ R _LF.RC/¢ POS _ T |01¢-7|ME || " i - ! * i • ; * ! - j . ! . ! . ! - i - i - | - ! - t - ! -

.)

.

o -| ../ -9 * ! ,tf , ! • ! • 1 , ! • t j • |t • I • | * | , t , | m t • _ , t • ! * ! * ! * t . _ • t • !

| .I .I .1 | ._) I.lr I.t 1.$ |,1 3.1 |.$ 1.11 q,$ S,l 11,1 _,11 I,l T_ME(SEC) 1).| .l ,$ 1,11 |,? |.t t.| |,S |,11 J.| I.t l.t |,It |.Jl !.| I1._) T_MEISEC) FIGURE 5.6 (CONT.) HL-10 OBSERVED DATA AND CONTROL SEQUE_;CETIME HISTORIES 5.2.3 Results of Flight 19-2 The first processing of the HL-10 data was with the maximum likelihood identification program in an output error mode. As such, it differed from the FRC Newton-Raphson identification program in only one way, the weight- ings of the fit errors. FRC uses a constant weighting matrix, based on an idea of the instrumentation accuracies. The SCl approach estimates the measurement noise covariance matrix and weights the fit errors by the inverse of this matrix. It is shown in Section 4.4 that the ML estimate of the measurement noise covariance matrix is given by the sample covar- lance of the observations. Edwards also uses a priori weighting on the difference between the parameter estimates and the wind tunnel estimates, which was not Inltlally included in the SCI maximum likelihood program.

The time histories of the five response variables along with the estimated values obtained after ii iterations of the data are given in Figure 5.7. Included also in Figure 5.8 is the fit error in the p and r observations. A comparison of the parameter values themselves and the wind tunnel values is given in Table 5.3 along with the associates standard devia- tions in the parameter estimates. The parameter values obtained from the FRC Output error method with fixed weights after 7 iterations and a priori weighting are also shown in Table 5.3, along with the associated confidence bounds.

As Figure 5.7 shows, the fits in all the observations were very good.

However, as is often the case, the fit error alone does not indicate an acceptable set of parameter values. The two major problems that appeared were that (i) the signs of the Lp, Lr, Np and N r derivatives had all changed from those of the wind tunnel values and (2) the fit error in * In this investigation, it was assumed that the wind-tunnel and theoretical values had correct signs. This may not necessarily be the case for lifting bodies flown at transonic speeds due to limitations of wind-tunnel testing and theoretical calculations. The question of how much confidence can really be placed in these values has not been resolved.

TABLE 5.3 HL-IO PARAMETER ESTIMATES AND STANDARD DEVIATIONS

FRC Values wlth Wind Max. ilk. Estimates IT&nnel Max. Lik. Estimates Max. Lik. gstlmates a priori weighting Max. Llk. _st_tee (with St'd. day.)

Parameter and Theo- (with confidence With a priori Weighting with Yp and Yr With a priori weighting and Biam (with St d dev.)

bounds) (with Stld dev.)

(with St'd dev.)

L -0.3435 -0.3436 _0.196) 0.915 (0.025) 0.395 (0.022) P -0.295 (0.0114) -0.271 (0.0114) L 0.2723 1.188 (0.0196) -1.363 (0.138) r 0.0671 (0.iii) 1.574 (0.0728) 1.349 (0.0747) -52.073 (1.18) L B -30;75 -56.489 (0.305) -46.94 (0.417) -47.179 (0.285) -50.32 (0.330) N 0.0245 0.0326 (o.o157) -0.160 (0.00429) -0.187 (0.00370) P 0.0380 (0.00352 0.0550 (o.o.o323) N r -0.1290 -0.1114 (0.0157: 0.432 (0.0187) 0.548 (0.0152) -0.111 (0.0114) -o. o89 6 (o. o168) N 8 6.8411 7.0496 (0.118) 8.523 (o.719) 7.474 (0.0978) 6.292 (0.0574) 6.845 (0.0570) Y -0.0617 -0.0584 (0.0039_) 0.329* (0.00307) 0_335" (0.00192) P 0.310 (0.00192) Y -0.0120 -0.0122 (0) "4 r -1.091' (0.0122) -1.064" (0.00754) -1.019 (0.00762) tn -0.0855 (0.00392) Y8 -0.0916 -1.471 (0.0202) -0.i458 (0.0161) -0.0949 (0.00383) -0.0918 (0.00382) L6a ;11.2464 11.996 (1.18) 12.415 (0.0731) 12.494 (0.0775) 12.124 (0.0586) 12.282 (0.631) 5.877 (1.18) L6a 5.665 6.288 (0.173) 6.544 (0.148) 6.252 (0.111) 6.429 (0.114) io 21.03 (0.302) 17.341 (0.222) 0.247 (0.117) -0.0137 (0.0951) N6a 0.8135 1.456 (0.118) 1.262 (o.o136) 1.245 (0.0201) 1.404 (0.0131) 1.357 (0.123) -3.178 (0_118) N6r -3.617 -3.186 (0.033) -3.313 (0.0326) -3.257 (0.0229) -3.194 (0.0212) N O -1.633 (0.0479) -1.206 (0.0401) 1.561 (0.0205) 1.562 (0.0,63) -0.0018 (0) Y6o -0.00180 0.0623 (0.00390) •0515 (0.00373 -0.0231 (0.00435) -0.0482 (0.00443) Y_r 0.0111 -0.00427 (0) 0.919 (0.00992) 0.0629 (0.00895) 0.0717 (0.00782) 0.519 (0.00768) Y o.513 (0.00898) o 0.412 (0.00778) -0.751 (0.0112) -0.639 (0.00829) Ptnitial 1.76 (0.158) 3.325 (0.146) 2.225 (0.116) 2.348 (0.117) Rlnltlal 0.117 (0.322) 0.0251 (0.0263) -o.71o (o:o21) -0.599 (0.0201) 0.398 (0.0188) _Inltlal 0.1305 (0.0152) -0.226 (0.0125) -0.219 (0.0111) 1.939 (0.0697) 1.767 (0.0511) -1.748 (0.0663) -1.429 (0.0571) Likelihood Function Value -2243 -2359 -1552 -2264 The identified qumntitiea are Yp + sin •, Yr- 508 • OL ROLL RArE (C_G/SEC) $1DESLIP(DEG) i I I i o

+i "

' ! I I I TLI I I I I yR_I RRTECOEG/SEC) -,d !

C, M O • I " I " I " I ' ! ' I " I " I ' i ' I " i • i " I ' I " I " $ o $ S $ S i °| -It . I , ! _ ! , I , I , I I I , I J I , ! i | . l , ! , I , ! .

• 8.0 .qi .8 !._ l .f, |.! 2.5 9.S 3.3 ].'/ q.] 4.S 4.$ 5.3 5.8 I;.Z T II'tE(SEC) I " i ' I ' I ' I ' l,i'l'l'l'l'l i'l_ i | I o -11 v z c) -4 e.4 i-.- G: I l . I , l , I , Ir . t , ! , ) , I , t , | i I . l , t , t .

-T ,11 1.2 1.¢ 2.) Z.5 2.$ 1.3 3.T 4.1 q.5 S.| .S.] 5.$ £.2 0.0 .S T II'IE($EC) FLg. 5.8 HL-10 FIT ERRORS IN p AND r MEASUREMENTS - OUTPUT ERROR the observations (although only p and r are shown) exhibited a sinusoidal characteristic. Although the precise reasons for these problems are not known, there are several contributory factors.

The first factor is that the parameters with the opposite signs are weak parameters -- that is, they cannot be accurately identifed from flight test data. This is the approach FRC has taken with the L and r Np parameters and in the a priori weighting, they weighed their wind tunnel values very strongly. There is also an identifiability problem because the stability augmentation system (SAS) was used on this flight_ Such a system would tend to suppress certain modes of the system, while emphasizing others. The parameters of the suppressed mode are, therefore, hard to identify.

A second factor could be that the linearlzed dynamics are not accurate enough for the flight conditions of this data. Also, there may be coupling between the longitudinal and lateral modes, which is not included in the model.

A third factor, which may account for the slnusoidal characteristic in the fit error, is that due to the instrumentation dynamics, the measure- ment noise is actually correlated. This hypothesis could be verified by reprocessing the data with the measurement noise modeled by a second order linear system. Final verification of these possibilities would have to be based on processing additional data under similar flight conditions.

In many instances, the solution to such problems is to adjoin to the likelihood function other measures of performance, usually indicating some a priori knowledge of the parameter values. Several examples of these are the a priori weighting and constrained parameter values which are discussed later. A less direct approach, of reducing the dlmensionality of the parameter space, is also presented.

5.2.4 Output-Error with Yp and Yr Identified The remaining series of runs were all aimed at solving the problems encountered with the first processing of the HL-10 data. First, the Y P and Y derivatives were considered as two additional parameters to be r identified.

In the previous run, both Yp and Yr were considered zero. However, by examining the equations of motion, p can be expressed as a function of (Yp + sin _) and _ can be expressed as a function of (Yr - cos _)r. This would introduce previously neglected second order effects into the estimate of p and r, and possibly account for the sinusoidal characteristic in the fit error.

The results indicated that the fit in each of the observations was about the same as in the straight output error case although, as shown in Table 5.2, L r does have the same sign as the wind tunnel value and Lp is less positive. However, on the other hand, both Np and N r are worse.

In addition, the sinusoidal characteristic of the fit error remained, diminished only slightly.

5.2.5 Output Error With Constrained Parameter Values Since it seemed clear that the opposite signs on the four parameters were a result of trying to minimize the fit error, and until additional 8O terms were added to the model, these signs were likely to remain opposite, the next run constrained the values of Lp and N r to remain negative. This would answer the questions of whether there was a set of parameter values which would minimize the cost criterion (although not globally) with the indicated parameters having the same sign as the wind tunnel values. If these values remained on the constraints, no such minimizing set exists.

The results of this run were that the L and N values remained on p r the constraints and the L r and Np values again had the opposite signs. In addition, the fit in the observations was drastically degraded. Only the fit on r was of equal quality as in the two previous trials.

5.2.6 Output Error With Different Initial Conditions One remaining possible cause of the changed values could be the presence of local minima having the opposite signs on L p, Lr, Np and N r.

The next run used initial parameter estimates of L and N which were p r more negative than the wind tunnel values, while keeping the control derivatives the same. The results of this run were that the signs of all four parameters (Lp, Lr, Bp, Nr) were again reversed and, if more itera- tions had been performed, the final parameter values would have, more than likely, been equal to the initial set of output error values.

Although the results from many sets of initial parameter values would be necessary to conclusively determine if the values from the initial output were truly the global minimum, it appears that this might be the case. If the signs on the four parameters are to be the same as the wind tunnel values, an additional cost must be put on the difference between the parameter estimates and the wind tunnel values. This is precisely the reason for the "a priori weighting" mentioned earlier.

5.2.7 Output Error With A Priori Weighting The values for the parameter weights used in this run were obtained

dlrectly from the FRC's runs supplied to SCI. Figure 5.9 shows the time

histories of the observations and the resultlng estimates (except for _) for this weighting. It is clear that the fits to the observatlonsD except for p, have been degraded. However, as shown in Table 5.3, the values of the four indicated parameters have the same sign as the wind tunnel values.

It was found, however, that except for p and r there was appreciable bias in the fits to the observed data. Not accounting for this bias in the computation of themeasurement noise covariance will cause incorrect weights to be assinged to the different observation residuals.

This will effect both the computation of the gradient and the information matrix, resulting in incorrect parameter step sizes. Another run was made with the sample bias of the observation residuals computed, at each iteration, and accounted for in the sample covariance calculation. The fits to the observed data for this second processing of the data with a priori weighting are shown in Figure 5.10. The fits to r, B and $ are much improved over the previous case. Only the fit to the lateral acceleration data has not improved. As shown in Table 5.3, many of the parameter estimates for the run are closer to the original wind tunnel values than for the previous run, without considering the biases. It, therefore, appears that when using a priori weighting, consideration must be given to the possibility of having biased residuals, which must be used in computing the sample covariance.

A final processing of the data with a priori weighting and including the identification of the output biases was made with the additional feature of retaining only the diagonal elements of the sample covariance for the estimation of the instruction noise variances. All the off-diagonal terms were set to zero. The rationale for this was that each of the measuring instruments on board the aircraft operate independently and there- fore the errors would be uncorrelated. The fit to the observed data did not improve over the previous run and many of the parameters were now

C9

$IDESLIP(DCG) | i e E _ "LNTERflL NCCELERRI ION(G) -( ,,] I '1 i _ I I I I I i i I I ROLl. RKI'E-TII_ " i • i • I • 6 • i - 6 • I " I • I " ! " I " | " i • L • i -'" N le II SB I I

\

,q o|l 4[ °1| , 4 I: ,.41

V

I,I . .I oi I I I.I I.I I.I I I I l I.I loO I I I.I I I I.I I.I i.l T il'lr (SEC) YRGI RnTE TIME II • ii I - • I "/__" i • i • I • i • I • i • i • i • i • $ " I ' b # i " V I I I c iiI IiI iiI I IiI Ii I III III iii iii iii I I I $ I I I I I II S TIMEISEC) SIDESL IP-TIHE l.l • i • i • ! " i • i • i • i - i • i " I • I • i - i • i • i - I.I I,I I.I I.I J .I l!l Iii ..... _ ....

I II II II II II II II I I II II II II I,I .I ,I TIHFCSEC) OUPUT ERROR WlTli A PRIORI HEIGtlTING AIID BIASES llL-I 0 : Fig. 5.10 ROLt.-TI_ ,.J ,J -Jr/ - _ ! i I , I • | • I • | • I _ I • | J ".' , ! • I • | _ ! • I , I 0.0 .q ol 1.2 1.6 too t._ 8.8 J.| 3.6 q.O t.4 !.0 $.2 |.$ i,0 ¢.q TIHE($EC) LAIERIb, L I_,CCELERRT l ON-T I_ • i , l , i , ! , l .

i ' i ' l ' I " _ ' l ' i ' I '. I " l ' .10 .o8 .Ol ^ .02 v • .oo F- w_" -.ot u - .OI t.--

,.. +

-.oI , i , I . J ! , I , ! , ! • t . I , ! . I , t . I , I . I . i O.I .q .I I.t 1.i 8.0 |.q l.i $.8 I.I q.I 4.4 t.i $.1 l.I $.0 I.q T IPIE <,SEC )

Fig. 5.1o HL-IO : OUTPUT ERROR WITH A PRIORI WEIGHTING AND BIASES

farther removed from the wind tunnel values.

No additional processings with Varied a priori weights were made, since sufficient data by which these variations could be Justified were not available.

5.2.8 Parameter Estimates Used for Prediction It has been often stated that using a set of estimates for the stability and control derivatives to predict the measurements from a flight test, under similar conditions and with similar instruments as the one used to identify the derivatives, would be the most valid test of the accuracy of the parameter estimates. Since another set of flight data for the HL-10, under similar conditions, was not available, an experiment was run in which only the first 227 points of data were used to identify the parameter estimates and these results were used to predict the final i00 points (2 seconds) of data. The identification algorithm which was used included the a priori weighting and the identification of the output biases. As the fits to the observed data, given in Figure 5.11, indicate, there is some divergence at the end, especially for r. However, the divergence in _ and _ was anticipated since the observed data suffers heavily from clipping during the final 1 second. The fit to the lateral acceleration was as good as might be expected considering the fit to the first 227 points, 5.3 M2/F3 Flight Test Data The data supplied to Systems Control, Inc. on flight No. 21, case 6 of the M2/F3 lifting body is sho_m in Figure 5.12. The influence of wind gusts is evident in the time histories of the sideslip angle and the lateral acceleration. Referring to Section 5.2.1, the wind gusts were assumed to enter the dynamical equations of motion in exactly the same manner as the sideslip angle _ • Nothing was known, a priori, about the ItOt_t. I_T(-T II_ -II

i"

o.e .q .e t_ t.| 8.6 z,e |,8 3.1 J.8 o.e i.t _.e $.| T|I'E(SEC) S.I a.o ¢.4 YPa_f RRT£ Tr_ | ' • i • i ....... .......

e i -4 -tl

a

-J

4 " ' I • I i I , I . I . I , I i t , I . I • I . I • l . I • I , I,I ,I .IT |I'$E ( ,__C )I'I I.I I.I.I.I I.I I,I I I I,I I,I I .I I.I I ,. _,_._I I I I I ROLL-T IlIE I . ! !

• i . i i t . t , t . i , i . i , i . i . I . l ! . i . | i I "ill ,I .I .I i+! 11 II Ii !1 I.i I I i I 41 i1 I II'E(SEC) I I Is e e e q FIGURE 5.11 HL-10: PREDICTION OF FINAL 2 SECONDS OF DATA

V_V_ _0 S(IN03_S Z XVNIH _0 NOI_DIa_BH :OI-IH ('IRO3) II'_ _flOlH

S'S 0"$ $'$ I'g I'te b'b O'b $'t" Z'¢ I1"1 F| 0"1 5"1 l'l I' b' 0"0 |0"- 7 ' I ' l ' I ' I ' I ' | ' | ' I " ' I ' | ' I ' I ' I ' | ' I ,'

• N

$0 "° _ r bO'- ('1 O.

F" prl " I _]D

,VA 00..°

- " I , z

I,O" SO" I|° ¢ I • I . I * I . I i I I i * I , I . I _ I • I , I • I , I • _I" 3N 11-NO111d8373DOU -I_I_13161"I (c)3g)Z_tl I |" I," 0"0 O'g ,1'$ i'g $'1, b't O'b $'1 1'¢ ll'l I,'Z 0'l S'I I*!

b'S , _ ' U '" IP'- 18) I ' I [_I ' I ' i ' i ' i ' I o I. ' ! ' I ' i ' ! ' ' i °- _ "0 0"0 O ÷ I- v ¢- |, $.

I'1 $oi |'!

I'| . I • I • I • I • I • I, • I • I • I n I • I • I • I * I • I * 31411-d17S3015 , I • ll|_||t|t|l,lllll,tl|ll|lll,l,l,l.i I .I 1,1 |.1 l.g I,i I.I I.I I,I t*l i.i t.1 I,1 Iol 1.1 I,_ I,1 ?*1 Y.I I,I TIMECSEC) FIGURE 5.12 O6 RUOOC._R POS | 1' I ON(O_G) LRTERRL RCCELERRT ION(G) .- - - - _ _ :.. :.- _. :.. _ .,...

: .... _, _____, , ,

i iI,,l i

I 1 I I I I == i i I i i I I . t ! I statistics (correlation time, mean square value) of the wind gusts.

The only information supplied SCI was that the output error program used by FRC had failed to match the time histories adequately. A total of 401 data points were supplied, representing 8.02 seconds of data. Once again, the effects of the instrumentation (quantization, clipping) were ignored as were the dynamics of the boom which measures sideslip angle (B vane). It was also interesting to note that the flight conditions were appreciably different than for the HL-10. The angle of attack was only 1.57 ° and the Mach number was .468.

Seven separate runs were made with the M2/F3 data, indicating a succession of possible model representations for the equations of motion.

Since neither the measurement noise nor wind gust statistics were known a priori, these were included, where called for, in the list of parameters to be identified, along with the stability and control derivatives and initial conditions.

5.3.1 Output Error - No Wind Gusts Included The maximum likelihood algorithm, in the output error mode, with the wind gusts assumed zero, was first used in trying to process the M2/F3 flight data. It was intended that from such a run, it would become apparent where the wind gusts were having the most impact and also the results would serve as a standard against which the identification algorithm performance with the wind gusts included, could be measured.

The time histories of the fit in each of the five measurements are given in Figures 5.13. As these figures indicate, the worst fits were obtained on thesideslip angle and lateral acceleration measurements, although none of the fits were as good as with the HL-10 data. These results also indicated that the model for including the wind gusts, suggested in Section 5.3.3 is appropriate, since the measurements involving the sideslip angle show the most random fluctuation when compared to the data from the HL-IO flight.

I-I I',o o _0 I-4 !

I ! i !

V

LRTERRL RCCELERAT I ON- .OC4 .OSG .@q8 .OqO .032

G

.O_C Z o .DOS LLJ O *OOO I -.OO@ -.0.15

g

.,J - .0|4 O.O ._ *@ 1.2 l.G Z.O 2.4 s'.8 |.t |,¢ t.0 q.q q,t $,t S,i_ S.@ s.q S @ I',t I'.E 0.0 TI_E(SEC)

Flg. 5.13(cont'd) M2/F3: OBSERVATIONS AND ESTIMATES -OUTPUT ERROR

The parameter values obtained for this processing of the data are given in Table 5.4.

5.3.2 Perfect Measurement of Sideslip Angle For this processing of the M2/F3 data, it was assumed that the measure- ment noise on the sideslip angle measurement is much smaller than the gust noise. With this assumption there is a perfect correlation between the process noise and the sideslip angle measurement noise, both being _n" The state and sideslip angle measurement equations now appear as C_ = Fx + Cu + rB n Y3 ffi_ + Sn + nB _ B +B n The Kalman filter for the complete four state, five output model must account for the perfect Bn correlation. The most direct method for doing this is to first construct an equivalent four state model which is uncorrelated with the sideslip angle measurement. This is done by adding the quantity y3-B-Bn, which has value zero, to the dynamics, i.e., C_ = Fx + Gu + FB n + 0 (y3-B-Bn) and solving for 0 such that E (FB n n { _ 0Bn) B T} ffi 0. O = r is seen to be the solution and the equivalent model has the resulting from C_ = Fx + Gu + r(Y3-¢) This equation is in the form of the Kalman filter, and it can further be shown that is the exact Kalman gain. Since is the third column of F, the dependence of _ on B is eliminated and the equations of motion and the measure- ment equations can be rewritten as TABLE 5.4 M2/F3 PARAMETER ESTIMATES AND STANDARD DEVIATIONS Wlnd Max. ilk. Estimate Mak. Lik Estimate Max. Ilk. Estimate -Max. Ilk. Estimates Wax. ilk. est_atas Parameter Tunnel -output error mode Assuming perf. B Meas.

Directly ident, of B With a priori veightln E with dependent parsms.

& Theo- (with St'd dew.)

(with 8t'd dew.) (with St'd dew.) n (with Sttd dev.) fixed, fixed.

retical (with St'd dev.)

0.679 (0.035) -i. 779 Lp -0.4673 -1.548 (0.0935) (0.214) -0.461 (0.0182) L r 0.8878 2.008 (Z.187) 10.49 (0.547) 25.46 (1.908) 4.154 (0.140) -97.79 (1.615) -135.38 L_ 75.140 -54.49 (2.45) (2.238) -67.95 (1.o2) -.0203 (0.00393 -.142 Np .0802 .i02 (0.006) (0.0147) .00475 (0.00349) N r -.6876 -.0307 (0.078) -1.675 (0.0590) 1. 628 (0.199) -.764 e (0.0134) 7.324 (0.152) -9. 890 N_ 7.5342 2.876 (0.136) (0.349) 6.763 (.0876) 4.435 (.113) y * P ¥ * r YB -.2001 -.0476 (0.125) -1.249 (.0597) -1.466

(0.0386) -.202 (o.oo392)

-1.36 (.0594) L6a 14.04 14.82 (0.301) 9.804 (0.109) 16.022 (0.3017) 10.96 (0.161) 9.66 (.169) L6r 10.03 73.97 (8.59) -109.28 (5.519) -157.130 (17.88) -42.18 e (3.13) L o 0 11.14 (1.828) -10.46 (0.328) .145 (2.11) -.572 (o.115) -9.004 (_141) N6a .83 .596 (.0223) .719 (.0104) .762 2.128 (0.0456) (o.ln) _D .756 (.0134) _n N6r -4.06 -12.874 (0.578) 6.844 (0.643) -Ii.754 (1.467) -4.37 e (0.106) N o 0 -.345 (0.121) .177 (0.0357) .427 (0.198) -.233 (0.0433) -.00239 (.0320) Y6a 0 -.00033 (0.0151) -.0363 (0.00669) -.0125 (0.00689 -.0847 (0.00867) -.0273 (.00634) ¥6r 0 .0301 (0.363) -.874 (0.222) -.926 (0.227) -1.932 (0.286) To 0 .0179 (0.354) .378 (0.0299) .295 (0.0313 -.0974 (0.37_ .456 (.0189) _b/ae -.281 (0.0531) 2.936 (0.0629) -6.01 (0.0933) -.667 (.108) Plnltlal 3.807 (0.521) •359 (0.188) 1.657 (0.852) 4.846 (0.296) -3.125 (.239) rlnltlel! -2.262 (0.0785) -1.66 (0.0280) -1.604 (0.0556) -2.22 (0.0453) -1.061 (.0597) _inltlal -.558 (.0251) - .565 (.0279) ¢i._Itial "34.44 (.175) -32.69 (.158) !-33.630 (.0796) -31.91 (.223) -31.52 (.224) -44.147 (.856) q s.231 (.107) 8_nltial 2.937 (.256) m Likelihood -1502 -2237 -2038 -1122 -1051 Function Value TABLE 5.4 (CONT'D) Max. Lik. with Rank Deficient Solution Parameter -.531 (.0189) L P 4.268 (.144) L r _103.35 (.105) L 8 N .0397 (.00682) P -.989 (.0672) N r 7.568 (.306) N B Y P Y r

-1.19 (.0590)

Y8 10.25 (.0845) L6a -5.539 (.0257) ,.o L_r

-10.89 (.280)

L O .561 (.0254) N6 a -.512 (.651) N6r .587 (.0833) N o -.0360 (.00660) Y6a Ygr -.737 (.219) Y .408 (.0296) O -.164 (.0428) ¢bias P -1.029 (.139) initial

-2.054 (.0552)

rinitial -31.576 (.0705) ¢initial Likelihood -1689 Function Value m w _ P L L 0 P 6 L 6 L 6 L o L8 p r a a r 4% C t _ ' r r + N_ N_ N O N Np N r 0 r r 8 a

1; 1 tan8 0 0

0 _o 0 m _Y3.

w n I I 'Yl "_ n 2 i y2 I Ir + n 4 + Y_ .6 +Y r o

rib_

I_Y5 J 8 " Y3 + Y_a_a r where Y3 is treated as a deterministic control. The order of the dynamical system has been reduced to 3 and the number of measurements to 4.

Once a complete set of parameters has been obtained for this reduced order system, the time history of 8n can be recovered. This is important since the identification of the statistics of the wind gusts is also possible using identification. Sideslip angle estimate 8 can be found by substituting the parameter values of the three state model into the original four state model and solving for its time history. Then subtracting from the sideslip angle measurement gives the time history of 8n + n 8 (n 8 was originally assumed small). A flrst-order linear model of the form _n ffi aBn ÷ vB (Y3 - _) = Bn ÷ n8 where v 8 is the process noise with covariance q, can be fit to this data and the time constant a, the process noise covarlance q, and the covariance of the measurement noise n 8 identified.

The time histories of the fit to the four observations (not in- cluding sideslip angle) are given in Figs. 5.14. The parameter estimates along with the estimates of the process noise covariance and the (recip- rocal of the) time constant for the wind gust are given in Table 5.3.

The time history of the wind gust 8n (including the neglible measurement noise) is shown in Fig. 5.15.

The fit in each of the four measurements is very good, although time histories of the fit error indicate that there is still the same slnu- soidal variation, especially in p, that was observed in the HL-IO fit errors. Only the fit error in the lateral acceleration, ay, approached being white noise, which is the indication of the best possible fit. The value of the covarlance of the noise on the sideslip angle measurement was almost two orders of magnitude smaller than the process no_se covarlance which supports the original assumption of this run. One surprising result, ROI.L(OE6) ROLL R_TE(DEG/_EC)

, _. _, .-

P

P

r vb u =

¢

"I ! I -I S , i , I , | ' i'' | I ' I ' I ' I ' '| _ | ' I ' | ' | ' j ' 1 ' | ' | ° | ' (degs) -3 -2 -3 , | t I, I i I, I , I I I, I t ] I I I I' I ' I , I , I , I , 1 , 1 t '* ,- -So.o .4 .8 1.2 ].G 2.0 2,_. 2,8 3.2 3.6 _..0 _._, _,.8 5.2 5.6 G.O _._ 6.8 7.2 7.G 8.0 TIIIE(SEC) Fig. 5.15 TIME HISTORY of B + n13 n I00 however, is that the L N_ and L_ parameters changed sign from the p, wind tunnel value, r r There appeared to be two principle reasons for the N5 and L_ r r parameters having the wrong signs. The first was that the magnitude of the aileron variation was much larger than the rudder variation, unlike in the HL-10 case. Since the effect of the controls is additive in determining p and r, there is an identifiability problem with respect to N5 and L_ . This is substantiated by the small values of the terms r r of the sensitivity matrix corresponding to the N_ and L5 parameters.

r r The second factor contributing to the incorrect signs is the opera- tion of the yaw damper. This causes a feedback loop which activates the rudder as a result of yaw rate. The time histories of r and _ appear r in phase, therefore, in the M2/FS time histories. In such a situation, unless the control d r is modeled as a linear combination of the states, there is a uniquesness problem as to whether the actual aircraft dynamics on the feedback loop is being identified. With the yaw rate and _ in r phase, the N_ parameter, at least, will appear with a positive sign.

r The problem of incorrect signs was of major concern and was the motivating factor for many of the remaining processing of the M2/F3 data. The experience with the HL-10 data indicated that constr_ning those parameters with wrong signs to have the same signs as the wind tunnel values would not correct the problem. The solution had to lie either in a more complete aircraft model or in dealing directly with the numerical problems causing the incorrect signs.

5.3.3 Wind Gusts Included: Direct Identification of Process Nois_C_v_rlance and Time Constant of Correlated Gusts For the third processing of the M2/F3, the gusts were included dir- rectly in the dynamical and measurement equations requiring that a full Kalman filter be used in the maximum likelihood identification algorithm in order to obtain the sensitivities. For this processing of the data, the model of the wind gusts obtained from the previous run was used

(5.3)

_n = a_+V6 is an unknown where a is the reciprocal of the time constant and v 8 in the original system disturbance with covariance q. Replacing 8n equations by equation (5.3) results in

-L u L 0 O- p

6a ! 0

L_ a L_ r L o

P r

&r 0

_[ N N N_ 0 0 r

N No

N6 a 6r

p r

l ÷ 0

v8

C _l = Y + Yr + YB _ 0 B

+ "Y_a Y_r Yo

s_n_ cosa v

_i 1 tane 0 0 0 ¢ 0

0 0 ,x, o

0 0 l

_n L 0 0 0 0 a B n

t_ has been added to the _ equation.

Note that an additional bias term o With t he inclusion of the covariance of v_ and the _n state, there are now 24 parameters to identify: 17 stability and control derivatives (Yp and Y are assumed zero), 5 initial conditions, q and a.

r The difficulty involved in setting up the identification algorithm in such a case is that both the measurement and process noise covariances, R and q respectively, are unknown. However, both R and q are needed in establishing the Kalman filter gain, which is assumed to be in steady state.

To begin the identification, therefore, some initial estimates of both q and R are necessary. R is assumed to be the diagonal elements of the sample covarlance matrix obtained from the output error method. An initial value for q is obtained from the results of the previous run.

Once the initial iteration is completed, the value of q is updated like any other parameter and the measurement noise covariance, R, is ob- tained from the sample covariance. This last fact is derived from the A property that with xili_ I defined as the Kalman estimate of the state at time t I given data up to time tl_l, N lim l )T.

_ (Yi-H_i l i-I " Dui ) (yi-Hxi I i- ]-Dui : HPss HT + R

i=l where Pss is the steady state error covariance matrix, obtained from solving a discrete Ricatti equation. The above expression is only approximate for finite data lengths.

The time histories of the observations and the estimates are given in Fig. 5.16, and the final parameter values are given in Table 5.4. Although the fits to the p, 8, _ and a measurements Y obtained from this run improved over those obtained from the output error method, they are not totally acceptable. It is interesting to note that the signs of the parameters L and L have retained the p r same sign as the wind tunnel values.

The time histories in Fig. 5.16 also indicate that most of the fits to the observed data are biased. Inclusion of measurement biases in the list of parameters to be identified did not, however, improve the performance.

STI3ZS',. T P(OIZO ) t i t | -'___'-,-__ _, a .=t e.. ° ,J_2_-'-- l.J

-e':,e

,_,_'_

o ROLL-TIlE -.IY • I • I • I • I • I • I " I • i • I • I • I • I ' I • I • I • I ' l • I • I • -Z!

-49 -SO -|!

-J!

-|q -95 • _ W D .J _ I, I, t,lil, l, I, I , I , | m I , I , I,I, r, I , I , I , f , -18 O.o .q .l 1.2 I.I l.O l,q 2,1 ).2 3.6 q.o 4.q q.l I.I S.l 1,0 l,q I+.I 'r,t Y,I 1.0 TIME(SEC) I._TERFII_ RCCELERRT ION-TIME .OE4 • I " I " l " l ' I ' I ' I ' I " l ' I " I " I " I " I " I " l ' I " I " I " ,0$| .Oq$ .010 .OSt .024 .011 b-n p_ .005 0.000 ILl IUl U O: --.008 -.011 O_ . .olq , I , I = I , l , | o l , l , l , I , | j l , f • I , l , I , l , I , | _ f , 0.0 .t .I t.t t.¢ 2.0 |.q t,ll I,Z 3.1 q.O 'I,II q.l i.t 5.I 1.0 6.q I.I 'r.l 'r.l 1.0 TIME(SEC) FIG. 5.16 (CONT'D): DIRECT IDENTIFICATION OF WIND GUST MODEL The reason the maximum likelihood algorithm with the Kalman filter cannot reduce the fit error any further is basically numerical. What has occurred is that the diagonal element of the measurement noise covariance matrix, R, associated with the measurement of the sideslip angle has become very small when compared to the covariance of the wind gust disturbance. Indeed, the wind gust itself is practically white.

As a result, the measurement noise cannot be distinguished from the wind gust. The alternative was to restructure the model so that the measurements of _ + 8n' the total sideslip angle are perfect. This is precisely what was done in the previous processing of the data.

5.3.4 Three State Model With A Priori Weightin_ The results of the two previous processings of the M2/F3 data indicated that the assumption of perfect measurements of the sideslip angle was reasonable and produced the best fits to the data. However, as stated earlier, the three state model resulted in wrong signs for many of the parameters. The first processing of the data in an attempt to correct these incorrect signs used the priori weighting technique.

The same weights as for the HL-IO data were used for the M2/F3. Measure- ment biases were included in the list _ of parameters to be identified since the use of a priori weighting on the HL-10 data indicated the need for bias estimation.

As shown in Fig. 5.17 the fits to the observed data resulting from a priori weighting were quite poor, especially in roll angle and lateral acceleration, althoughthe parameter estimates themselves, as given in Table 5.4, were quite close to the wind tunnel values. It is interesting to note, however, that L 6 still has a wrong sign.

r 5.3.5 Three State Model With Fixed Parameters The basic causes of the incorrect signs for some of the parameters are that either the sensitivity of the output to changes in that parameter are small, as indicated by a relatively small diagonal element in the information matrix, or that there is a correlation, with respect to the sensitivity, between two or more of the parameters being identified.

Such a situatlon would be indicated by an off-diagonal element of the normalized information matrix being close to _i. If this were the case the correlated parameters could not be individually identified.

Both these problems existed with the M2/F3 data.

One technique which has been used for treating both these problems is to fix one or more of a set of parameters that are correlated.

The results of the identification run with fixed parameters indicated (i) the convergence of the algorithm to the final set of parameters estimates is more rapid and monotonic than when all the parameters are being identified, and (li), the final fit to the observed data is degraded to a certain degree. This latter characteristic is due to the fact that the number of degrees of freedom (equal to the number of parameters to be identified) for fitting the observed data has been reduced. In comparing the value of the likelihood function or cost for two cases with different numbers of parameters being identified (measure- ment noise covariance R being identified in both cases) the comparison should be made between a corrected cost, given by

f(N,k)AnlR l

y_ _ T1_ ..i.l.k-l'l'l'l" I,Inl.l,l,i,l'l'l'l'l'l'll ..,i.l,l-i'l'l'l" .$.1.i.i-i,1"|'1"1"111 .I .q 1.1 °|.t v ,,l.q

'".'.'.'.i.'.'.'.'-'.'.'.'-'-'-

_'1"111.| .q o! I,:{ I.I I.I I qq 1 |1 |1! |11 11| qlq ¢11! i1| |11 IIi iii iiIi t I I| I I • • ! I ! ! I ! I ! ! ! I | TIME(SEC) "151.1 .q .I 1.| I.I I,I l,q I.I 1.1 |.1 t,I q.q _l.I |*l I,I Iol I*q Iol lf.t 1Pol I.I TIt'_(SEC) U_TER_. RCCELERAT I OhI-T | I'E POLL-TIME ,Ill • i " i - 0 - 1 , I • i " i " i • i • i - i - i - i , i i i i "" OD .1_ . i - i - i - i - i - I - i - i - i • i i i i i i i o111 -Is .Iql -II .a_e "11 .111 -I 1 *oI1¢1 ,Ill III " ! "Oil .lit , t, lit. I,I, I, I'1"1"1" II I'!" I" I'!'1" ..,,, tol ,i .I l.l I.I I.I l.I l,l I.I |,1 q,l q,i ¢ol I.t I,I I.I I,! I,i f,I f.I I,I il.|,l,li|.l. I, |11 ||11"111"111"1|11|11" n, -II *1 LI i*1 l,I I.q I*1 I,I I.! 4ol q.q q,I I*1 I,I I.I I.q I.I t*l Iol I,t TII_aE(SF.C) l.I *q T U'_(SEC) .--------- ol)sez'ved (!.al:a FIG. 5.17 _/F3 TI_ HISTORIES WITH A PRIORI WEIGHTING + + + estimate where N is the number of data points, and k is the number of parameters being identified and f(N,k) is a monotonically increasing function of k which depends upon the objective of identification. For certain types of systems, Akaike (Ref. 38) has shown that the one-step ahead final prediction error using the identlfied model is given by J = N-k _n IR_ N+k A plot of J vs k for a typical model is shown in Fig. 5.18. The important thing to note is that the predictive qualities of a model do not improve monotonically with the number of parameters, even though the fit error decreases monotonically with the number of parameters It was indicated earlier that there was very little variation in the rudder during the M2/F3 flight, which would make identification of r the _ derivatives very difficult. In fact, the identified values of the L 6 and N_ derivatives with the three state model were physically unreasonable, _eing opposite in sign from the wind tunnel values. The first processing of the M2/F3 data in this set of runs was therefore made with the rudder derivatives fixed at the wind tunnel values, with measurement noise biases being included in theunknown parameter set.

The results indicated a strong correlation between L and N and P r almost all the other parameters, and a fairly poor fit to the data.

Fixing the same parameters but including a priori weighting did not improve the performance.

The second processing of the data in this series included fixed P r .

L and N derivatives as well as fixed L 6 and N 6 The results showed only a slight improvement over the _reviou_ processing, and a strong correlation still existed between L r and Lp; N and and P N_; Value of Criterion Akaike's / , Fit Error Criterion !

Optimal No. Number of of Parameter Parameters FIG. 5.18 PERFORMANCE CRITERION AS A FUNCTION OF THE NUMBERS OF MODEL PARAMETERS ii0 Np and YS" It appeared, at this point, that a problem of identifiability was developing between the measurement biases and the initial conditions.

deri- The next processing of the data therefore included fixed L Nr, d r p' vatives and initial conditions. The fits to the observed data were still very bad and correlation still existed between Lr and LS, Lr and Lo, and Np and N .

After several more experiments it was decided that L L r, LB, Np p' N and the _ derivatives should all be fixed and the initial conditions r r be identified instead of the measurement biases. The results, shown in Fig. 5.19, were the best fits to the observed data obtained with the technique of fixing parameters at the wind tunnel values. The values of the parameters which were identified were all of the same sign as the wind tunnel values, as shown in Table 5.4.

5.3.6 Three State Model With Rank Deficient Solution The results of the data processing, using the technique of fixing selected parameters of a correlated set at fixed values, showed that the convergence rate was improved due to better conditioning of the informa- tion matrix. The basic reasons why the parameter fixing technique does not always work are: (i) the correlation is usually not simply between pairs of parameters, but may involve the entire set of unknown parameters, and (ii) it is not usually possible to correctly choose a set of para- meters that should be fixed and the values at which they should be fixed.

It was decided at this point to investigate, more fully, the problem of possible correlation between more than just pairs of unknown parameters.

The solution of this type of dependency problem is to find the directions in parameter space corresponding to combinations of parameters which cannot be identified. A perfect dependency among the parameters would, strictly speaking, result in a zero eigenvalue of the information matrix, causing it to be singular. However, since round-off and other numerical / 4" "4" -I- + errors prevent the information matrix from being exactly singular, all the eigenvalues will be non-zero with a spread between the smallest and largest eigenvalue being many orders of magnitude. In such a case it is better to use a rank deficient solution for the inverse rather than a full rank solution. That is, the inverse to the information matrix should be computed leaving out one or more of the smallest eigenvalues.

The reason for neglecting the smallest eigenvalues can be explained in terms of the parameter step. The eigenvalues of the information matrix are the dimensions of the uncertainty ellipsoid, associated with the parameter estimates, in parameter space. The smaller eigenvalues indi- cating a larger dimension and therefore more uncertainty. Since the LxL information matrix M can be expressed in terms of its eigenvalues and eigenvectors, %i and Vi, i = I, . . , L L M = %i Vt viT L=I the parameter step is given by (see Section 4.3) L m Ap = M -I (D J) = _L-1%1 i Vi viT where DJ is the gradient of the likelihood function. Therefore, the smaller elgenvalues also contribute the largest proportions to the parameter step.

This implies that the largest components of the parameter step are in the direct- ion of the most parameter uncertainty. Therefore, the information matrix inverse is computed neglecting a certain number of the smaller eigenvalues, i.e. L-K M -I = i_l= %1 -I V i viT Each eigenvalue which is left out relates to a singular dlrection in parameter space, and, therefore indicates a combination of parameters which cannot be identified uniquely. Rather than fix the value of one or more of the parameters, as was necessary with the a priori weighting technique, the rank deficient solutio_ fixes combinations of parameters corresponding to nearly zero eigenvalues. Thus the dimension of the space in which the set of parameter values that minimizes the cost are sought, is reduced. It is important to realize, however, that at each iteration the values for all the specified un- known parameters are assigned updated estimates.

The number of eignevalues to be neglected depends on which order rank deficient solution produces a parameter step resulting in a set of parameters with the lowest associated cost. The procedures is as follows.

Starting from some inimimum number of eigenvalues (i0 in the HL-10 case), a reduced rank inverse is computed, the parameter step isdetermined and the associated cost evaluated. One more eignvalue is then added and the procedure repeated. This same thing is done until all the eigenvalues are added in, with the last inverse being a full rank inverse. The same procedure, starting from a minimum number of eigenvalues and progressing to the full rank, is repeated every iteration.

The fits to the observed data are shown in Fig. 5.20. Comparing these with Fig. 5.14, it can be seen that the fits are only slightly degraded. The fits to roll angle and lateral acceleration are much better than those obtained for either a priori weighting or fixing of correlated parameters. The parameter values obtained for this third order rank deficient solution (3 eigenvalues neglected) are given in Table 5.4. Although several of the parameter still have opposite signs from the wind tunnel values, many of them are much closer to the wind tunnel values and are more reasonable than the full-rank, 3 state parameter es- timates. Some, such as N 6r, now have the correct sign from physical considerations, where before they did not. It is clear that further development work on this rank-deficient solution approach will improve the estimates even more.

It should be mentioned that the basic identifiability problem in the M2/F3 data is due to the stability augmentation system (SAS) providing all the sudden movement. The methods described above (viz. a priori weighting, fixing parameters and rank-deficient solutions) are indirect means of handling this problem. It would bemore exact to model the relevent characteristics of SAS in order to alleviate the problems.

ROU. l_qTE-TIl_ m tm_nr T_ " i- i. i-l_l-l, i-i_i, l.l.l.l.l.l_ i . i I. i - i .

_ II-1- i i i -1 - i - i - l- i - i- i - l - i - li. i*- i- i-i- l _ its Io eJ i -o_ -I.II W C] v _I,Q -18 w Q: ..I Q: "qI "I'l* t,I, I, t. I.I, t, t,I. t. t. 1, I, I, OI|** t.

)- -e.I " t' t'_'_I" I" I' t, I. I. t, t, I, I* t. l. t- I- _l Ij • • I.J .q .| |.l !.$ Z.I |._ m.I |._ $.g _,I 9,9 q.I S.E g.I I.I O._ • $ _._ ?.I I I TXI_(S_) • • TIn_(S(C) I--" ,110 .lie ul| e.Ill "i o41$ .=1.6 -,led v .J _ *.I¢4 • f_t, 1.t,t,f. f. T . t,T, I.t.l, l,l.l. 9at_t, • I,t, I, 1,1. f_l.l.t*l.l,l_l.I.l.l. !,_ i.

I.I ,q *J I.I 2.8 _,e I,l |,1 |*1 I.I 0,0 q.q q,$ 1.8 IJ I,I I.Q I.I I.I _*i t,I I*l .e *l I*1 |*l £.1 I.t I*l I.I I*1 I*t _*t _*1 _.1 I*1 I*_ I*e I*_ _,1 _.1 D,I TI/'_(S_¢) T|r_(¢(¢) observed da_a + + ÷ estimate FIGURE 5.20 M2/F3 TIME HISTORIES WITH RANK DEFICIENT SOLUTION VI BACKGROUND FOR LINEAR SYSTEM INPUT DESIGN The importance of choosing appropriate inputs (i.e., control surface deflections) for exciting specific modes of an aircraft or executing spe- cific maneuvers has been recognized for a long time. Several conslderations which enter into the selection of inputs for an aircraft are: i. Pilot Acceptability - The inputs should be capable of being im- plemented easily by a pilot and the resulting response of the aircraft should not endanger pilot safety.

2. Parameter Sensitivity - The measured response of the aircraft should be sensitive to the parameters that are being identified.

This is necessary for obtaining good estimates of the parameters from the flight test data during the inverse computation or the identification process.

3. Instrumentation Limitations - The dynamic range of the instruments and their slgnal-to-noise characteristics impose limitations on the types and magnitudes of aircraft maneuvers. The relationship between input design and instrumentation specification has been emphasized in (Ref. 40).

4. Derivative Extraction Method - In the past, the choice of control inputs has often been dictated by the desire to use a particular method for derivative extraction. For example, sinusoidal inputs were used initially to obtain the transfer function of an aircraft at specified frequencies (Ref. 40). However, it was soon realized that this was very expensive in terms of the total flight test time required to obtain the aircraft stability and control derlvatives.(Ref. 40) Next, the step and the doublet type of inputs were used and specialized

methodssuch as Prony's Method (Ref. 41) and the Time Vector Method(Ref. l)

were devised to extract derivatives. With the more powerful digital techniques available today such as the Newton-Raphson (Ref. 8) and the Maximum Likelihood Methods, (Ref. 33) arbitrary inputs can be handled and it is no longer necessary to limit the inputs for the success of the derivative extraction method.

, Modeling Assumptions - The slx-degree-of-freedom equations of motion and the nonlinear aerodynam/c model for an aircraft contain a large number of parameters (over 20C). The simultaneous esti- matlon of all these parameters from a single maneuver is not attempted since this would lead to nonuniqueness and identifl- billty problems. Generally, llnearized decoupled equations of motion are used for the extraction of longitudinal and lateral stability and control derivatives. The inputs selected for exciting these modes should be such that the assumptions of llnearlty and decoupllng are not violated. The inputs currently in use are mostly of the doublet type. The resulting aircraft response is an impulse-type of response about a given trim condition. Generally no attempt is made to optimize the frequency, the shape, or the timing of the impulses in order to make the air- craft response sensitive to the parameters that are being identified.

The motivation for the present study comes from a simulation of the X-22 VTOL Aircraft perfomed in 1970, and described in Section 5.1.4. The multlstep input gave parameter estimates which a;e an order of magnitude more accurate than the estimates obtained using the Cornell input. At about the same time, one of

the authors developed a general theory of optimal input design

based on the techniques of modern control theory. (Ref. 42) In the followlng section, the salient features of the theory and the computation aspects of input design are presented. The results of applying it to the C-8 aircraft parameter identification problem are considered in section VIII.

6.1 Related Work on Input Design in System:Identlfication The importance of input selection for system identification has also been recognized for a long time, though a unified mathematical treatment has emerged only recently. Some of the earlier attempts at input design were based on fre- quency domain methods and engineering judgment. A large amount of literature exists on Pseudo Random Binary Sequence (PRBS) inputs which have been found to provide improved identification for a large number of systems. (Ref. 43-45) However, PRBS inputs use very little information about the known properties of the system.

Since in a number of physical systems some a priori information is available about the modes of the system (e.g., short period mode, phugoid mode, etc. of an aircraft's longitudinal motion), one can use this information to design inputs for identifying these modes more precisely.

The work described in this report is most closely related to that of Aoki and Staley, (Ref. 46) Levadi, (Ref. 47) Nahi and Wallis, (Ref. 48) and Levin (Ref. 49) on input signal design for system identification and to that of McAulay (Ref. 50) and Esposito (Ref. 51) for signal synthesis. Aoki and Staley (Ref. 46) consider single-input, slngle-output discrete-time systems. Levadi's results; (Ref. 47) are only applicable to the case in which the unknown parameters enter linearly in the system impulse response. Levine's results (Ref. 49) are applicable when linear regression is used to estimate the unknown parameters.

The results presented here are applicable to multi-input, multi-output, continuous time systems. Tile computational algorithms proposed are new and and have not been used earlier for input design purposes.

VII THEORY OF INPUT DESIGN FOR LINEAR SYSTEM IDENTIFICATION* The problem of inputdesign for linear system identification is formu- lated here as an optimal control problem• The performance criterion used is the sensitivity of the system response to the unknown parameters. The other criteria for input design such as pilot acceptability, instrumentation characteristics and state deviations described in Section VI are considered indirectly through an energy or power type of constraint on the input and through modifications of the final results. It is assumed that an output error or maximum likelihood method which can handle arbitrary inputs is used for derivative extraction. In these methods, the measured response of the system z(t) is expanded in terms of the system parameters as follows: z(t) = y(8,t) + n(t)

(7.l)

where y(t,e) is the true system response, 8 is an Nxl vector of unknown parameters and n(t) represents measurement noise• Let 80 be the best a priori guess of 8. By a Taylor series expansion of Equation (7•i), z(t) = y(O o, t) + V 0 Y(0o, t) (@ - Oo ) + n(t) + higher order (7.2) terms where In the output error method, the step (0-00) is determined by a least squares solution of Eq. (7.2) over the time interval [0,T]. It is intui- tively clear that the sensitivity function V e y(80,t) must be sufficiently large to achieve high accuracy in determining the parameter step (8-80) .

Furthermore, the sensitivity functions must not be linearly dependent over 0 _ t < T for a unique solution of (0-00) The theoretical aspects of this work were supported in part under AFOSR contract No. F44620-71-C-0077. See "Dual Control and Identification Methods for Avionic Systems - Part II, Optimal Input Design for Linear System Identi- fication".

For input design, we now define a scalar performance index in terms of the sensitivity functions. For a single parameter, a suitable performance index is the weighted mean square value of the sensitivity function over the time interval [0,T].

P

dt (7.3)

Jl _ .__12/ T 2 i,,i where o_ is the variance of noise in the measurement Yi" The perform- ance index J1 favors the measurements which are more accurate and makes them more sensitive to 81 compared with the mesurements that are less accurate. It can also be shown that Jl represents the inform- ation in the measured response about the parameter 01 • For the multiparameter case, the choice of a scalar performance index is much more complicated due to the conflicting nature of the individual sensitivity functions. One possible choice is a weighted sum of the individual sensitivity measures.

N J -- _ wj J] (7.4)

j=l

where the weighting coefficients wj reflect the relative importance of parameters. It is clear that the criterion (7.4) can also be written as the weighted trace of the Fisher Information Matrix M.

J = tr (_) (7.5)

where W=

(7.6)

T T

and M=

dt

_0 _0 (7.7)

I a2 2 P_.

(7.8)

° o° I

Other criteria such as the weighted trace of M -I the smallest elgenvalue of M or the determinant of M lead to nonlinear optimi- zation problems which are much more difficult to solve.

An alternative method for choosing a criterion function is to con- sider the objective of identification. For example, if the identified parameters are to be used for control system design, the inputs should be selected so as to minimize the control error. Similarly, if the parameters are used for response prediction, the input design should be based on minimizing the mean square response error. However, these cri- teria generally lead to optimization problems which are mathematically intractable or extremely difficult to solve. It is also felt that the advantage gained by solving these difficult problems may be more than offset by the basic uncertainty about the initial parameter values 80" The input design process described here should be viewed as a sequential design process in which tile inputs are selected based on the current best knowledge of the parameters.

7. I Problem Formulation

Consider the linearized aircraft equations of motion (7.9) _- FX+ Gu where x is nxl state vector and u is a mxl control vector. F and G are n x n and n x m matrices of unknown parameters (stability and control derivatives). The measurements are denoted by a p x i vector z(t) which is contaminated with white noise :n(t) (7.10) z(t) = Hx(t) + n(t).

H is a p x n matrlx and n(t) is a zero mean Gaussian white noise process (7.11) Z{n} = 0, E{n(t)nT(T)} = R_(t - z).

Let O denote the N x i vector of unknown parameters in the above equations. It is required to select the input {u(t), 0 _ t _ T} to maximize the weighted trace of the information matrix M subject to an energy type constraint. The information matrix M for the unknown parameter set e can be written as (7.12) (V e x) T H T R -I H(Vex ) dt The energy constraint on u(t) is (7.13) uudt = E o T T The performance index J can be written as, J = hr{WM} = tr{WI/2MW I/2}

(7.14)

# = tr{ (Vsx " WI/2) T RTR -I H(V8x • W I/2) dt}

o In order to use this criterion, the quantity (Vox) must be calculated as a function of time. A differential equation for (Vsx) can Be easily found from the system equation. Since the multiplication by W I/2 represents a column operation, this is equivalent to calculating successive equations for 1/2 ax wl ao-_ ' i = I,...,N This can be accomplished by the following: @ , 1/2 ax , 1/2 aF , 1/2 ax , 1/2 aG

(7.15)

_w i tw i _w i _i -_i ) = x +

+

with _ = Fx + Gu.

An equivalent way of formulating the problem which makes use of the simultaneous computation of w_/2 ax/a8, involves the specification of aug- l mented FA, GA, and H A matrices and x A vector. These are defined as follows (with m the number of inputs and p the number of outputs): N

-F I0 I lO-

aG

aF I F I I 0

wl_ I 1 I

Wl

q (7.16) : I : I I: aF aG

_ o I IF

WN_N I I I

(N+I) nx (N+I) n (N+l)nxm eo a

-- T T = 0 R__

• t o (7.17) H A = -r- "7 "..,%.

H 0

o o

Np x Np Np x (N+l)n and m X _X _x (7.18) x A =

w2-_2

_x (N+l)n x I With these definitions and using Eq. (92), it is possible to write: (7.19) XA = FAXA + GAU The performance criterion is now redefined as T (7.20) j ffi XAHARA HAX A d t Io T T -i In the next section, we obtain a set of necessary conditions for the optimal input using the Pontryagin's Maximum Principle (Ref. 52).

7.2 Optimal Engergy-Constralnedlnput UslngMax_mum Prlncip!o Let A(t) denote the costate vector for XA(t) and _ be the constant multiplier associated with the constraint (7.13). The Hamiltonlan of the augmented system is, 1 (xA)THAT IHACXA) _ E CAU ] (7.21) _= 2 [- RA + _(uTu _)] + AT[FAXA + the necessary conditions of optimality are or (7.22) and C/e = 0 u or u* = -1 (G A )T X (7.23) The boundary conditions are homogeneous.

XA(0) = 0 , ),¢r) = 0

Substituting for u* in Eq. 7.19, we obtain the two point boundary value problem, Since the boundary conditions are homogeneous, the solution is trival vlz. xA E O, k E O, u E 0 except for certain values of _ which are the elgenvalues of the two point boundary value problem. In other words, the problem is of the Sturm-Liouville type.(Ref- 52) The eigenvalues and the optimal input can be determined in a number of ways. Two possible methods are (i) the transition matrix method and (ii) the Riccati equation method.

7.2.1 Transition Matrix Method denote the transition matrix of 7.24 for a particular _.

Let _(t,0;_) (7.25) _(t,O;_) = exp Partition _(t,0;_) into X A and X parts as follows: (7. 265 Then (7.275 %(T) L_xx(T,O; _) The second equation in (7.275along with the boundary conditions gives (7.28) X(T) = _%k(T,0:W) X(O) = 0 For a nontrivial solution (7.29)

I+xx(T,0; )I = o

Eq.(7.29)is the elgenvalue equation for the Hamiltonian system (7.24). It is a nonlinear algebraic equation in _ and can be solved by a Newton-Raphson iteration. In general, there is an infinite set of eigenvalues, but we will be only interested in the largest value of p which Will be shown to maximize S (Section 7.3).

7.2.2 Rlccati E_uation Method

The eigenvalues _ are functions of the interval length T. Therefore, one can fix p and determine T for which ¢%%(T,0;p) becomes singular.

Another way is to use the Riccati matrix P(t) defined by the relatiQnship xA(t) = P(t)%(t) (7.30) An equation for P(t) is obtained by differentiating both sides of Eq. (7.30) and substituting from Eq. (7.24).

x A = ;>+ P; or 1 T [FAP-_ O A G A or I GA (7.31) P(0) = 0.

(7.32) The Riccati Eq. (7.31) differs from the usual Riccati-equation of the Llnear-Quadratic problem in that the forcing term (last term) in Eq. (7.31) enters negatively. Eq. (7.30) can also be written as X(t) = p-1 (t) xA(t) X(T) = o, t = T, since whenever p-1 exists. At final time p-1 (T) ffi O. (7.33) which means that a conjugate point exists at t = T.

Eq. (7.33) provides us with a method to obtain the critical interval length T corresponding to an eigenvalue _. The Riccati Eq. (7.31) is integrated forward in time for a particular _ using initial conditions(7.32). When the elements of P(t) become very large, the critical length T correspond- ing to an eigenvalue is being reached. Now p-l(t) is integrated using the equation d p-I p-i

(p-l)= _

or (7.34) d I p-i GA p-i = - p-IFA- FIp -I + HTRAIHA + _ GA T (p-l) At the critical interval length T, all the elements of p-i go to zero.

It follows from the Sturmian property(Ref. 52) that the smallest T corresponds to the largest eigenvalue _.

After the critical length T corresponding to the largest value of has been determined, Eq. (7"24) is solved forward in time using _(0) obtained from Eq.(7"28)and(7.29) as an eigenvector of #xI(T,0;_) corresponding to

the zero eigenvalue. Thereby, the boundary condition %(T) = 0 is auto-

matieally satisfied. A unique value of %(0) is found by using the normal- ization condition of Eq. (7.13).

7.3 Application of Functional Anal_s%S I In the last section, the optimal input u was characterized in terms of the solution to a two point boundary valueproblem. In this section, we show that the optimal u is an eigenfunction of a positive self-adjoint operator corresponding to the largest eigenvalue p.

Let A denote the operator corresponding to Eq. (7.19) viz.

t A[u] =i eFA(t-T) GA U(T)dT (7.35) Let A*[-] denote the adjoint operator to A[.]

A*[w] = GAT ; TeFAT(s-t) w(s)ds (7.36) "t Let <u,w> denote the inner product <u,w> = uT(t)w(t)dt (7.37) The performance index J can be written as I = <.xA HTRAIHA XA> T-I

- <Au, HAAS>

(7.38) The energy constraint of Eq. (7.13) is written as <UtU> " E It is well known that J is maximized subject to the above constraint by u which is an elgenfunction corresponding to the largest eigenvalue of the operator A*HIRAIH¢" Furthermore, since A*HTR-1HA A is a positive self- adJoint operator, all its eigenvalues are real and positive (Ref. 52). For finite T, the operator is also compact and has a finite maximum eigenvalue. The * H_RA RAA corresponding to this eigen- optimal u is the eigenfunction of A* T -i value and normalized according to <u,u> = E.

(7.39) A A Also, (7.40) Max Jffi vE U To show the relationship of the above elgenvalue problem with the two point boundary value of Eq. (7.24), define (7.41) Z = Au and (7.42) m ft T eFJ (s-t) w(s) ds A , Then, using the definition of A and (7.43)

= F Z+GAU , z(0) = 0

A T -1 • (7.44)

. _ FTn + .AR, .,z.n(T) = o

From Eq. (7.39) and (7.42)

(7.45)

or

1 _T n (7.46)

P Therefore

z (0) = 0

- 1 GA%T ,(7.47)

n(T) = 0

_ - - FTn + HTR-AIHA z , A comparison of Eq. (7.46) and (7.47) with Eq. (7.24) shows that Z = XA 7.4 Examples We now apply the above results to two first o.rder and a second order example.

7.4.1 First Order System with Unknown Gain Consider the system = -x + Ou (7.48) where x and u are scalars and 8 is the unknown gain.

y = x+v s (7.49) where E{v} = 0, E{v(t)v(r)} = r _(t - 3).

From Eq. (7.23) and (7.24) (7.50) u = ---1 tl

d__ e_ "-1 - _ (7.5l)

dt Ei__genvalues Equation 7.51 gives

=! v0x+ x

r _ IX]+ 1

_: i [_ve x _-_ [¥ re. + x]

r (7.52)

= C1-!-) x •

_r It 2 a > O, a = 0 and _ < O.

Let a = I - i/_r. Three cases arise: 2 2 2 can be easily shown that only a < 0 leads to a nontrivial solution.

Let (7.53) w = - 1 Then k(t) = C I sin wt + C 2 cos wt Vox(t) = r[Clw cos wt - C2w sin wt - C1 sin wt - C 2 cos wt]

vex(O) = o=> clw- c 2 = o

k(T) = 0 => Cl[sin wT + w cos wT] = 0 For nontrivial solution,

tan wT = -w

(7.54)

Eq.(7.54) is the eigenvalue equation. The smallest root w of Eq.(7.54) o corresponds to the largest value of _. The optimal input u is obtained from Eq. (7.50) as , C 1 u = [sin w t + w cos w t] _o o o o (7.55) = 2 where _o i/r(l + Wo). Notice that u satisfies the same second order dif- ferential equation as _ viz. Eq.(7.52). C 1 is determined from the condition I T u2dt = E.

o This gives E Cl = /2 - 1/4 sin 2T (7.56) 7.4.2 Levadi's Example Levadi considers the following example in his paper (Ref. 47).

_-- -x+ bu

(7.57)

y = X+ V

(7.5S)

where v is a correlated noise process with autocorrelation function.

E{v(t)v(T)} = c e -'alt-Tl It is required to estimate b only.

Levadi's (Ref. 47) results can be easily derived as follows: v(t) can be represente d as (7.59) where e is a white noise process and E{E} = 0 , E{c(t)c(_)} = 2ac 6(t - T) A new measurement can be generated by differentiating Eq.(7.58). (This procedure is similar to that of Bryson and Johansen, Ref. 53).

= (a - l)x + bu - ay + s Now 9 can be regarded as a new measurement which has white noise _ in The new information matrix is: it.

T r --!-i [(a - l)VbX + u] 2dt (7.60) M = Jo 2ac where _b x = -VbX + u (7.61) Now the problem is in the same form as example 7 .4.1 except that the performance index is slightly different.

The following equations of optimality are easily derived.

u =- 1 [_ _ a- 1

- z "VbX] (7.62)

2 (p

= .(a - 1) 2

ac VbX + a - 1 u + A (7.63) ac An equation only in terms of u can also be obtained from Eq. (7.61) - (7.63) u- [1 a -1 2pac- 1 -] u = 0 . (7.64) The eigenvalue equation is tan(wT + _) =w a (7.65) a - 1 ]1/2 where w = [-1 + 2 ac - (7.66) = tan -I w The optimal input u = A sin(wt + @) . (7.67) Notice that the results for example 7 .4.1 can be obtained by letting a ÷ _ and 2c/a _ r, where 2c/a represents the area under the autoeorrelation function of v.

The optimum value of w is chosen to maximize p.

From Eq. [7.66), 1 a - i P = 2a---_ [1 + 2 ] " (7.68) l+w It is seen from Eq.(7.68) that when a > i, the maximum of p is attained for the smallest value of w. This corresponds to the case when the noise is wide-band. For the narrow band noise case viz. a 2 < i, the second term in Eq.

(7.68)is negative and the maximum of _ viz. i/2ac is reached at w = =. The

practical implication of this result is that the input should be of as high

frequency as possible. Since the noise is narrow band, this increases the

high frequency signal to no_se ratio at the output.

7.4.3 Second Order Example

The following system represents the short-period longitudinal dynamics of an aircraft.

(7.69) e

°jLx j

Y=Xl+V x I represents pitch rate, x 2 is a linear combination of pitch rate and angle- of-attack and 6 is the elevator deflection. The transfer function of the e system from elevator to pitch rate is Xl(S) blS + b 2 (7.70)

6 (s) 2

e s - als - a 2 Optimal Input for Identifying bl: Sensitivity equations are (7.7].)

e Optimal input, 6 e = - Ii/2_ (7.72) 2 Xl_ llal_ _2a2

_I ° ¥ Vb1

(7.73) _'2 = -tl " Both the optimal input and X 2 satisfy fourth order linear differential e equations of the form z= [D4 + (p_r- 2a2- al)D2 + a2]l 2 0 (7.74) where D denotes the differential operator d/dt.

A solution for 12 is easily written as X2(t) = C 1 sin Wlt+ C 2 cos Wlt+ C3 sin w2t + C4 cos w2t (7.75) where Wl = 7_ - 4a 2

"2 + - 'a2-

Jo

2 2 2 a = (_ - 2a 2 _ al ) = Wl + w2 (7.76) Also w I w 2 = a2 .

It is assumed here that o > 2a 2 or i/_r > (4a 2 + al ) since other cases lead to hyperbolic functions which become unbounded for large T. They are rejected as possible solutions using the same reasoning as used in examples 7.4.1 and 7.4.2.

The expressions for Xl(t), VblXl, u(t) and Vb2X 2 are easily obtained from Eq.(7.75)using Eq.(7.71)-(7.73). The elgenvalue equation, assuming Wl _ w2 ' is obtained as sin wlT cos wiT sin w2T cos w2T - 0 (7.77) -w I cos wIT w I sin wiT -w 2 cos w2T w 2 sin w2T a 2 + w 2 -alw 2 a 2 + w 2 wlal w I 3 w2 2 _ -- wla I w 2 - w2a 2 - _-_ 3 2 w2a I w I - wla2 _r Equation(7.76)is used along with (7.77) to select w I and w 2 where w 2 = - a21w I.

From Eq. (7.76) which maximize _.

I__= a2 2 2 - w2 )2 + a_ • pr (Wl + _i ) + al = (wl The minimum of I/gr is attained at w_ = -a 2 or Wl = w 2, i.e., at the undamped natural frequency of the system. However since w I = w 2 is ruled out by the , solution considered here, the root of Eq.(7.77) closest to the undamped natural frequency should be chosen. Since WlW 2 = -a 2, the two frequencies will bracket the natural frequency of the system.

The optimal input 6 e is ,C7.78) * 1 6e = 2-_ [ClWl cos wit - C2w I sin Wlt + C3w 2 cos w2t - C4w 2 sin w2t] CI , C2, C3, and C 4 are determined as the eigenvector of Eq.(7;77)corresponding to the root w I. They are normalized using the condition T 62 dt = E e o 7.5 Extension to S_stems with Process Noise Consider the linearized aircraft equations of motion = Fx + Gu + rn (7.79) !

Z =Hx+v

where n(t) is a gaussian white noise forcing function representing random gusts

and modeling errors, z{n(t)} = o , z(n(t)n_(_D = q_(t - _) .

The information matrix M in this case is given in terms of the Kalman filter for the above system.

R = FR + Gu + K(y - HR) (7.80) K = EHTR -I (7.81)

E = FT + IF T + FQF - EHTR-IHE

(7.82) where R denotes the best filtered estimate of x and E denotes the covariance of 3. The Kalman filter provides a linear causally-invertible whitening transformation for the process y since the innovation sequence (y - H_) is a gaussian white noise sequence. The likelihood function is easily written in terms of the innovation sequence (Refs. 24, 27). The information matrix M is given as fT HTR-IH(vo_)}dt M = E{(Vo_)T (7.83) where Vs_ denotes the sensitivity function of the filtered estimate _ with respect to the unknown parameter vector 8. Note that both K and E are func- tions of 8 so that the sensitivity equations are much more complicated than for the no process noise case. Moreover M, in general, depends on the random quantities q and v so that its expected value needs to be maximized.

A special case arises when 8 contains parameters from G only and the initial state is known exactly. Since K and P do not depend upon G, the sensitivity equation has a simple form V8 _ = (F - KH)VS_ + V8C u (7.84) K is, in general, time-varying, but if the system is completely controllable and observable, K reaches a constant steady-state value (Ref. 54). Then Eq.(7.84) is essentially similar to Eq. (7.19) and most of the theory developed in Sections 7.2 and 7.3 carries over to this case.

7.5.1 Example: Let us consider example 7.4.1 with additive process noise.

= -x + eu + n (7.85) y = X+V where E{B} = 0 , E{_(t)_(_)} = q_(t - T) , x(O) = 0 • The filter sensitivity equation for 8 under steady-state filter gain, k > 0 is (7.86) VOi = -(i + k)Vsi + u , Vsi(O) = 0 Proceeding as in example 7.4.1, and defining I (7.87) 1 _ ]112 = [3 (I + k) 2 or P = r[,,_ 2 + (i + k) 2] it is seen that the optimal input u obeys Eq. (7.54)-(7.56). Notice that by increasing process noise q, the gain increases and _ decreases. Thus the information M = _E for the same input energy E decreases. The frequencies m, however, remains unchanged.

7.6 State-Variable Constraints

Linear state variable constraints can be handled either directly by adding a quadratic penalty function to the performance index or indirectly by adjusting the total input energy E. The examples 7.4.1-7.4.3 show that the amplitude of the input u is determined by E. Thus by adjusting E, the amplitude of the input u and the state x can be bounded. Of course, the inputs obtained in this fashion are not strictly optimal.

7.7 Steps in Optimal Input Program As outlined in Sections 7.2 and 7.3, there are several possible com- putational techniques which can be used to solve for a consistent pair of interval length T and largest elgenvalue _max" The method that has been implemented in this program is to numerically find the first time instant at which the solution of the ricatti equation, relating XA(t ) and %(t), becomes infinite. The instant at which this singularity occurs can be found with arbitrary accuracy by continually using a finer and finer step size and noting the instant at which the diagonal elements of the Riccati solution change sign.* This change in sign of the diagonal elements is one key indication that the ricatti solution has blown up through one direction (e.g., -_) and returned from the other (e.g. +_).

The complete flowchart of the computer program designed to calculate optimal inputs is given in Fig. 7.1. Many of the detailed steps have been combined into a single descriptive step since their description is beyond the scope of this report. For example, the actual technique used to inte- grate the ricatti equation is explained elsewhere (Ref. 55).

The only instant in the computational algorithm where the theoretical development was not followed exactly was in the calculation of the eigen- The alternative technique of calculating the largest eigenvalue of the Hamiltonian system is much more difficult since the determinant of the transition matrix for the llamiltonian system exhibits a very sharp zero.

1NPUT: W, F, G, HI n, N, p. m, Q, R,

aF aF

"_1' " ' " ' aes

_G _G CONSTRUCT:

I

F A, GA, HA, RA P m.O INTEGRATE RICCATI EQUATION UNTIL SIGN CHANGE: o iterate _ accuracy TI4AX CALCULATE : !

F A I-c Qc T

I A A Z m T-1 I T

! '_R A HA I-FA

.i

I

CALCULATE ¢_ " EXP(Z'A)

*T " ExP(Z.T_)

PAY, TITION

| 7

I

CT " [@22 I

I I(N+l)n

(N+l)n 1]4 2 CALCULATE: Ao = elgenvector assoc, with smallest.elgenvalue of IAo_ specified by input energy requirement I. ,, DEFINE: CALCULATE : Sz(t+A ) ° _ASz(t) control: u(t) = -QGATA (t) energy: E = E+(uTu) * A CONSTRUCT : i

,""

(VoX)' =6 I " a • ; I _-_N] e_ _62 I i

$

CALCULATE

, , 'THTR-IH(VoX ) , M = M +(VoX ) ._

M+W-I/2M'W-I/2 iterate until t _ THA X Fig. 7.1 FLOW-CHART OF OPTIMAL INPUT COMPUTER PROGRAM vector of _xx(T,0;_) associated with the zero elgenvalue. The problem en- countered was basically numerical. Very seldom, if at all, would the matrix _x(T,0;_) exhibit an exactly zero elgenvalue. This was caused, to a large extent, by the fact that T is never the exact instant of singularity of the ricatti solution. This difficulty was resolved by choosing _(0) such that the product _xx(T,0;_) %(0) had the smallest norm. If _%%(T,0;_) did, in- deed have a zero eigenvalue, this would be the associated eigenvalue and _XX _(0) E 0. The solution of this minimum norm problem is to choose %(0) to be the (normalized) eigenvector associated with the smallest eigenvalue of _k T_kx. This was the technique incorporated into the optimal input program.

Two additional items should be mentioned concerning the steps in the computational algorithm. The first involves specifying some additional fact about the eigenvector _(0) so that it can be uniquely specified. It was shown in Section 7.2 that the contro% at any time point, t, is a linear fune- tlonal of l(t). Therefore, if a particular input energy is required, the norm of _(0) is scaled to the proper amount.

The second item concerns the reconstruction of the M matrix as the last step in the computer algorithm. To just calculate the performance index, tr{WM}, it would be sufficient to use T T T -I tr{ XAHAR A HAXAdt} • o However, the information matrix M, the Cramer-Rao lower bound M -1, and the det(M) give important information about the identification of the parameters which is not reflected in tr{WM}. For this reason, the augmented state vector, x A, is rearranged to construct the matrix (Vox)WI/2, which is used to compute the information matrix M. In addition, the eigenvalues of M and a flgure-of-merit tr[M}• tr{M -I} is also computed.

7".8 Specialized Algorithms This section describes the two specialized algorithms which have been incorporated into the main program. These are: (i) the algorithm for com- puting the optimal set of weights when using the product of diagonal elements Of M criterion, and (il) the algorithm for computing _max when the data length, T, is specified.

An Algorithm for Choosing Weights: As outlined in Appendix B, one purpose of the weighted trace is to obtain a closer approximation to the determinant of M as the performance criterion. This is done by bringing the ratio of the largest to the smallest diagonal element of M as close to 1 as possible. The computer algorithm takes the form of an iteratlve sequence of choosing weights, calculating the resulting M and then updating the weights and repeating.

The formula for the updating of the weights, once an optimal input has been computed along with the accompanying information matrix, is as follows: W new i Woldi + e[_ -!-I - ] mii Woldi (7.88) where m''ll iS the i th diagonal element of M. Additional logic was subsequently added to the program to reduce the factor e by successive factors of 2 if the new set of weights failed to reduce the ratio of the largest to smallest diagonal element of M. This was made necessary since the equation (7.88) for updating the weights is by no means optimal.

An Algorithm for Determining _ for a Given T: max The second special algorithm built into the program enables a user to determine the Umax for a specific length of data, T. The most direct, but costly, way of determining a value of Umax would be to pick several values of _max and run through the program, finding the associated values of T.

These pairs (Umax,T) could then be used to construct a _max from this curve, the _max for a desired T could be determined. However, this procedure would require a great many (_max,T) pairs. A more direct method is to employ the optimal input program with associated zero crossing logic which would converge onto the correct _max' For the simple case of (_I,TI) and (_2,T2) as known associated pairs and T 1 > T d > T 2 where Td is the desired data length, the iterative equation for successive choices of _d is (_d2 E _2 ) (Ref. 56).

1 = (Ti - Td) 1 (T 2 - T d) _ . (7.89) Ti - TI _i Ti T1 _di+l - _di With a new, updated choice of _di, the program is run and an associated T i is found. This new pair is used in Eq. (7.89) to find an updated Bd ' and so on. Once the change in values for _d becomes smaller than some 7, the procedure is stopped. For c of 1%, this procedure usually requires only four or five repetitions. Of course, if the value of T d is outside the range of the given initial pair (_1,T1) and (_2,T2), the logic of Eq. (7.89) is altered appropriately.

Examples of both these specialized algorithms are given in Section VIII.

VIII

APPLICATION OF OPTIMAL INPUTDESIGNTO C-8 AIRCRAFT

This section discusses the results of applying the optimal input

design technique described in Section _II to a two and four state model

of the longitudinal dynamics of a C-8 aircraft. In both cases the input was designed to optimize the identification of the five parameters associated with the short period mode. Additional constraints were put on the input in the form of limiting the maximum excursion of the angle of attack and pitch rate hy adjusting the total energy of the optimal input. In the first part of this section the results of the two-state example are described, including a frequency domain analysis of the optimal input, a comparison of the input performance versus a suboptimal doublet elevator input and use of the weighted trace criterion. The second part of this section describes the results of a Monte Carlo simulation of the identification process for the five short period parameters from simulated flight test data for both optimal input and a suboptimal doublet input.

8.1 Short Period (Two-State) Longitudinal Dynarics of C-8 Aircraft 8.1.1 Optimal Input for two-state model This investigation involved finding the optimal elevator deflections 6 e to identify the parameters in the two dimensional short period longitudinal equations of motion for a C-8 aircraft. The state variables are the pitch rate q and angle-of-attack._. The equations for the short period dynamics of the C-8 aircraft were = + 6 - 0.737_] [_0.005_] e and the measurement equations were and the values given In determining the power spectral densities of nq n e, in Ref. 13 (i deg./sec, error in q, and 2 ° in _) were multiplied by two times the correlation time of the noise sources, which is assumed to be 0.01 secs. The measurement noise spectral density matrix is therefore given as R Ii 02 _.041 For this investigation the data length T was fixed at 4 secs. The ° appropriate _max was found from a _max-T curve shown in Fig.8.1.

The value of _max associated with a T of 4 secs is 0.015. The shape of the _max-T curve in Fig. 8.1 is characteristic of the general relationship between these two variables.

For the _max and R values indicated above the optimal input with the respect to the three parameters in the F matrix and the two parameters in G is given in Fig. 8.2. The energy of the input was 311 and tr{M} = 20,460. The check value of _max E was approxi _ mately 20,200, indicating a numerical error of 0.1%. The determinant of M was computed to be i x 1015 , with the ratio of the largest to smallest .eigenvalue of M being almost three orders of magnitude. The eigenvalues themselves indicated that much of the relative uncertainty in the parameter estimates was concentrated in two of the five dimensions.

The standard deviations of the parameter esti-,ates were .167 0.063 Standard for G = Standard for F =

E 0 l r°°l

Deviation 0.035 J Deviation L0. 025J The time histories of the states _ and q, resulting from the optimal T I0 %0 | I 8 I I I I l 2 3 4 6 7 8 9 10' 5O 5O 5O 5O 5O 50 50 50 50 50 8,1 -- Fig.

VS. T CURVE FOR A 2 STATE/5 PARAMETER MODEL TIME, SEC.

1 2 3 4

+-0 -4 -6 Elevator Deflection, _ 8 • 8 e - 10 (Deg.)

- 12 - 14 't OPTIMAL INPUT E = 311 Tr{M} ffi20460 OPTIMAL INPUT FOR SHORT PERIOD LONGITUDINAL DYNAMICS Fig. 8.2 - PITCH RATE - q /DEGS_ 1 2 3 4

TIMZ ~(SEC)

ANGLE OF ATTACK a

(DEC)

6 .

i ! !

I 2 3 4

TIME ~(SEC)

Fig. 8.3 - PITO! RATE AND ANGLE-OF-ATTACK TLXIE HISTORIES WITH OPTIMAL INPUT

input are shownin Fig. 8.3. The energy of the input was constrained

@ so that _ does not exceed I0 . This method of energy limitation is the most direct way of applying state constraints; although, as mentioned in Section VII, the penalty function approach can also be used.

8.1.2 Fourier Transform of the Optimal Input Since it was specified that the input be designed with respect to the parameters both in F and G, it would be reasonable to infer that the input.

would have a low frequency component for identifying the parameters in G and a higher frequency component for identifying the parameters in F. The actual Fourier transform of the input signal is given in Fig. 8.4. The vertical scale has been reduced to I/i0 its actual height in order to illustrate the smaller variations. The DC component is 0.98. The small peak in the transform occurs at a frequency of 0.375 cyc/sec, which is close to the short period frequency of the aircraft.

The most important point demonstrated by the frequency domain anaiysls is that the sinusoidal llke component of the input signal occurred at a very specific frequency. It is well known that the maximum signal power can be obtained from a second order system if it is excited at its natural frequency.

Therefore, in order to maximize the sensitivity of the output signal to the parameters in the F matrix (which is given by the M matrix), or in other words, to maximize the component of the output signal due to the parameters in F, the input signal had a specific component set at the natural system frequency.

8.1.3 Comparison with a Doublet Inut_al Duration and Energy The third part of the investigation was to compare the performance of the optimal input to that of a doublet input of equal duration and energy. The doublet input used for this comparison is shown in Fig. 8.5.

•05

Ln (.m [-_ Z O

1 2 3 4

5 6 7 8 9 I0 II

FREQUENCY (I/SEC)

Fig. 8.4 - FOURIER TRANSFORM OF OPTIMAL INPUT

Elevator Deflection (DEGS) 22.7

• 3 .6 4.0 T_(SEC)

-22.

FIG. 8.5 SUBOPTIMAL DOUBLET INPUT The determinant of M was 7 x 108 , seven orders of magnitude less than for the optimal, while the trace of M was only 1600. The standard deviations of the parameter estimates were Standard for F = for G =

[.17 271 stadard [ °°3

Deviation .164 ; Deviation 0346 Two of these values are smaller than those obtained for the optimal input, however the standard deviations in the parameter estimates are not explicitly in the optimization criterion.

_1.4 Effect of Small Parameter Value Chan_es on Optimal Input Since, in an operational application, the actual parameter values are not known, it was important to investigate the effect that changes in the parameter values had on the optimal input shape. These changes might reflect, for example, the difference between wind tunnel or theoretical estimates of an aircraft's stability and control derivative and the actual stability and control derivative values. It would be these estimates, however, that would be used to compute the optimal input for the identification of those same derivatives. The modified F matrix which was used in this investigation was -2 -4 -6 -8 -10 -12 -14 ELEVATOR DEFLECTION (DEG) Fig. 8.6 - OPTI_L INPUT FOR SYSTEM NITH 10% PARAmeTER VARIATION F -1.429 -.613 1 -.663 The optimal input resulting from this approximately 10% change in the parameter values is shown in Fig. 8.6. The most Obvious difference is that the input length is now 3.21 seconds instead of 4.0 secs. Otherwise, the input has the same qualitative shape as the optimal input. The expected stan- dard deviations in the parameter estimates resulting from this input when it is applied to the original system are as follows: ".255 .0886 Standard for G = Standard for F = ; Deviation .026

l ['14 I

Deviation .036 [ J Three of these values are very close to the standard deviations obtained for the optimal input, while the other two represent increases of approximately 50%. The improvement over the suboptimal input, however, still exists.

8.1.5 Weighted Trace Criterion: This part of the investigation involved using the weighted trace criterion to derive the optimal input and choosing the weights to make the diagonal elements of the information matrix equal. As detailed in Section VII the performance criterion is tr{WM 1 where W is a diagonal matrix of weights chosen to set _imll = ... = _ m .

P PP When all the diagonal elements of WM are equal, maximizing the trace of WM is equivalent to maximizing the product of the diagonal elements which is a better approximation to det(M).

Since the input which maximizes the performance criterion depends on the values of the weights, which in turn affect the input, an iterative scheme was used to update the weights until convergence was achieved.

ELEVATOR DEFLECTION .I .2 .3 .4 .5 .6 .7 (OEG) 02 • I , I -4 -6 - 8 - 1o Fig. 8.7 - OPTIHAL ELEVATOR DEFLECTION WITH UNITY NEIGIITS

J

5- STATES qs a (DEG/SEC, DEC) 1.

•1 .2 .3 .4 .5 .6 .7 TX_ (s_c) Fig. 8.8 OPTIHAL STATE TIHE HISTOILIES FOR U_IITY IdEIGHTS The optimal input and the state time histories for a T of 0.77 sec and an energy of 62.61, and unity weights are given in Figs. 8.7 and 8.8. The trace of M for this input was 611 and the determinant of M was 8.11 x 103 . The ratio of the largest to smallest eigenvalue of M was 1270. The eigenvalues themselves indicated that the parameter uncertainty was quite disproportionate along two of the eigenvector directions. The standard deviations in the parameter estimates are given below.

Deviation for F = Standard for G = D Standard 12.03 5.57_ V0.381 _ i. 43_ Deviation L0.129_ Notice the increase in standard deviations due to a shorter data length (0.77 sec. vs 4 sec.) and smaller input energy (62.6 vs 311).

Using the weighted trace criterio_ and 12 iterations to bring the ratio of the largest to the smallest element of WM down to 1.14, the optimal input and state time histories given in Fig. 8.9 were obtained.

The determinant of M was calculated to be 2.27 x 104 which is five times greater than the unity weights determined. The volume of the uncertainty ellipsoid decreased by the same factor. Another indication of this was the fact that the ratio of the largest to smallest eigenvalue of M was reduced by a factor of 2 to 640, with the largest eigenvalue itself being reduced by a factor of 2. The standard deviations for the parameter estimates are as follows: Deviation for F = Standard I Standard 11.56 4.04_ F0. 304 _ 1.23_ Deviation for G = L 0"122_ Notice on Fig. 8.8 that the parameters which are poorly estimated with unity weights were assigned higher weights. As a result, the lengths of the error ellipsoid along each of the axes has become more uniform and the total volume of the uncertainty ellipsoid has been decreased.

O .2 ELEVATOR 4 INPUT (DEG) 6 - • X.O " / eelghts = [29,4 ) I0 • s .273 # q STATES q,Q (DEC/SEC,

VEG)

TL'_E(SEC)

OPTI,_b%L INPUT AND STATE TI.X[E I[ISTORIES Fig. 8.9 - WIT}[ NEI(;ilTED TRACE It is clear from this example that using a weighted trace criterion does result in an input which can, in an overall sense, identify the unknown parameters with improved accuracy. This improvement is measured by the increase in the value of the determinant of M.

8.2 C-8 Monte Carlo Simulation A more realistic test for verifying the advantages of the optimal input over a sub-optimal doublet input of equal energy was to perform a monte carlo simulation of the identification process using simulated flight test data. To make the simulation representative of an actual flight test the following were included. First, both phugoid and the short period modes in the lonzitudinal eouations of motion were used to generate the data, although only the short period derivatives were to be identified. The optimal input had the characteristic of suppressln_, as much as Dosslble, the Dhu=oid mode in order to maximize the sensitivity of the output to the values of the short period parameters.

Second, short period parameters of the four state model that were used to compute the optimal input were changed by approximately 50% from the equi- valent parameter in the model that generated the data. In this way the situation of designing the input based on the wind tunnel values of the stability and control derivative was simulated.

There were several criteria for comparing the performance of the optimal input with that of the suboptimal doublet input. Since the optimal input was computed on the basis of maximizing the trace of the expected information matrix, this is an obvious candidate. Maximizing the determinant of the information matrix, or equivalently minimizing the volume of the uncertainty ellipsoid, in parameter space, is another. How- ever, by the nature of a monte carlo simulation, the parameter estimates themselves can be used to calculate the sample covarlance, and its trace, determinant or the eigenvalues can be used as performance criteria.

Finally, histograms of the parameter estimates themselves and the associated probability distributions can also be used in determining input performance.

8.2.1 Optimal and Subopt_m-I Inputs for Monte Carlo Simulation The four state longitudinal equations of motion for the C-8 aircraft which were used for computing the optimal input are m N U 33.74 -.02 -32.2 0 0 0 1 6 e + -1.658 -.562 q .003 0 -1.588 .005 -.737- Gt -.004 0 1 n u Zlk Uk

qe

z2 k 8 k + n q z3 k qk n .125 .125

61j

where R = E{nln_ j } .125 .25C is the discrete measurement noise covariance. The optimal input itself was designed with respect to enhancing the ability to identify the five short _parameters, given in the locations marked by an x in F : |x x

x

These parameters are Mq, Ma, Za, M5 and Z_ . Considering the e amount of computer time needed to generate the optimal input for five parameters and four state variables, a 2 second data length was decided upon. The optimal input, for an energy of 58.0 is shown in Fig. 8.10, along with the suboptimal doublet input, of equal energy, which was used for comparison.

8.2.2 Generation of Simulated Flight Data The different stability and control derivatives which were used in the generation of the simulated flight data are illustrated below.

!

Same !

-f F _ !

Same -2238-28 G=1829!

area I

i 368 L007

!

The random noise added to each of the four measurements were derived from a Gaussian random number generator. These noise sequences presented a slight problem since, with a 2 second data length and a .02 sampling period, the i00 samples were sometimes insufficient for the noise statistics to have the desired mean and covariance. However, since the performance of the optimal and suboptimal inputs were always compared for the same measurement noise sequence, the problem of incorrect noise statistics should be of m_nor importance for comparison purposes.

• 5 1.0 1.5 2.0

TIME (sFc)

ELEVATOR -5 DEFLECTION (DE_) -10 Optlm_l Input ELEVATOR DEFLECTION 2.0 (DEGS)

TIME (SEC)

-10 Suboptimal Inpu.t.

Fig. 8.10 - OPTIMAL AND SUBOPTIMAL INPUT FOR MONTE CARLO SIMULATION _2.3 Desc_ion of Monte Carlo Identification Simulation A flow diagram of the Monte Carlo identification simulation is given In Fig. 8.11 . Fifty different sets of random numbers are added to the computed time histories in order to simulate a statistically significant set of runs at the same flight condition. Each of the 50 sets of simulated flight test data was processed in exactly the same manner using the SCI maximum likelihood identification program. The parameters to be identified consisted of the five unknown stability and control derivatives and the measurement noise covariance matrix. The initial stability and control derivative estimates were set equal to the values used in computing the optimal input, while the initial state estimates were set equal to the value of the first data point, since all four variables were being measured.

The estimate of the measurement noise covariance matrix is obtained, at each iteration, from the sample error covarlance (fit error covarlance).

The four diagonal elements of the sample error covarianee are taken to be the estimates of the four independent measurement nolsesequences.

For each set of data the identification procedure is carried out until either the change in cost, between two iterations, or the norm of parameter step size, or the norm of the gradient becomes smaller than some present tolerances. For both the optimal and suboptimal inputs, convergence, indicated by one of these three conditions being satisfied, was usually obtained after 4 iterations through the data. At the end of processing a set of data, the final parameter estimates and the information matrix would be stored on tape for use in the latter compilation of the results.

Since there wasn't any process noise being considered in the simu- lation, the maximum likelihood identification technique reduces to an output error method, with the weighting matrix being the inverse of the sample error covariance. The required calculations consist of, for a given set of parameter estimates, integrating the system equations and the sensitivity equations. Since there are 4 states and 5 parameters, this makes a total of 24 simultaneous differential equations. The information matrix is computed from the time history of the sensitivity functions.

A priori

Generate Parameter Es timate_ Suboptimal Input Optimal Input 4-State C-8 Dynamical Model Generate State Actual Parameter Time History Values Actual Noise Add Noise Covarianees to Data

sEqUeNCE

OF STEPS IS Initialize Identifi- A priori Paramet REPEATED 50 cat/on Routine TIMES - FOR Estimates _____.._._....___, EACH INPUT I Identify Stability I & control deriva- I t_ves with maximum | likelihood identi- Matrix Parameter Estlmles_rmati°n Compile Results I for each input ] Fig. 8.11 BLOCK DIAGRAM OF MONTE CARLO SIMULATION

It is interesting to note that since the weighting matrix is the

sample error covarlance, the value of the likelihood function will be determined completely by the log of the determinant of the sample error covarlance. There are also two different methods for choosing the weighting matrix. One is that it should be considered a diagonal matrix since the measurement noise sources are all independent. The other is that, non zero off diagonal terms be allowed in the weighting matrix. For this Monte Carlo simulation, the former method gave better results. In the latter case, the error covarlance matrix had a tendency to become singular at times.

For this application, none of the options discussed in Section 4.6 , including parameter fixing, rank-deflclent solution for M -I or bounding the parameter estimates were used. Neither was it necessary during any of the iterations to cut the step size in order to improve the convergence properties.

8.2.4 Results of Monte Carlo Simulation The monte carlo simulation consisted of 50 runs of the identification procedure, with the parameter estimates, information matrix and its elgenvalues, and the parameter covariance matrix and its elgenvalues being saved at the end of each run. The ensemble results are tabulated in Table 8.1.

The theoretical values of the trace of the information matrix, using the actual values for the stability and control derivatives, were computed to be 2.12 x I0 7 and 4.74 x 105 for the optimal and suboptimal inputs respectively. The average values, from Table 8.1, for the 50 runs were 2.15 x 107 and 4.79 x 105 , indicating that 50 runs were sufficient for obtaining accurate parameter estimate and information matrix statistics. In addition, the trace of the sample covarlance, computed from Sample Covarlance =____I (Apj-A_) (Apj-A_) T

J=i

MONTE CARLO RESULTS BASED ON IDENTIFICATION TABLE 8.1 FOR 50 SETS OF SIMULATED DATA Suboptimal Input Optimal Input .315 .242 Trace of Sample Covariance 1.501 x 10 -15 1.62 x 10 -19 Determinant of Sample Covariance .262 Eigenvalues of Sample Covariance .514 x 10-1 .725 x 10 -2 .115 x 10 -2 .252 x 10 -3 .766 x 10 -4 .188 x 1Q -4 I .234 .126 x 10 -5 .202 x 10 -7 .307 .40?

Parameter Standard Deviations .491 .295 .0235 .00925 .0537 .0771 .00168 .000570 Average Trace of Information 4.79 x 105 2.15 x 107 Hatrix 4.79 x 105 2.14 x 107 Eigenvalues of Average 8.46 x 103 2.95 x 104 Information Matrix 4.77 x 102 6.56 x 103 2.18 x 101 1.39 x 102 4.14 1.12 1.95 x 1014 4.70 x 1018 Average Determinant of Information Matrix .312 .182 Average Trace of the Covariance Matrix (Cramer-Rao Lower Bound) .303 Lower bound on parameter .441 .234 .0311 standard deviations (from .00876 .0568 I .351 .0665 Cramer-Rao Lower Bound) .00262 .000247

Apj

where is the error in the parameter estimates for the jth run and Ap is its

mean. As stated, it was found that the sample covariance matrix was fairly close to the Cramer-Rao lower bound.

By almost all measures of performance the optimal input performed better than the suboptimal input. The determinant of the sample covariance, giving a measure of the overall parameter uncertainty based on the actual derived parameter estimates, was four orders of magnitude smaller for the optimal input. The eigenvalues of the sample covariance were smaller, on a one-to-one basis, for the optimal input, indicating a smaller dimension for each axis of the uncertainty ellipsoid.

The histograms of parameter estimates are shown in Figs. 8.12 - 8.1.6 For the M , Z and Z 6 parameters, the optimal input definitely pro- duced a better ensemble of parameter estimates. The mean estimate value was much closer to the actual parameter value and the standard deviations and mean square errors were smaller. The performance for the two inputs was about the same for the Mq, while the suboptimal input did outperform the optimal input on the fourth parameter. (M s ) I{owever, e it should be kept in mind that the accuracy of the parameter estimates themselves was not a direct performance objective. Rather, the overall input performance, as measured by the determinant or trace of the covariance matrix was the criterion of interest.

It is also interesting to note that the standard deviations for the

parameter estimates obtained from the inverse of the information matrix (Cramer-Rao Lower Bound) did in fact bound the parameters standard deviations obtained from compiling the individual parameter estimates. For the optimal input, a comparison of these standard deviations is given below: Actual Standard Deviation Lower Bound From of Parameter Estimates Information Matrix Inverse (Cramer-Rao Lower Bound) .407 .351 .295 .234 .00925 .00876 .0771 .0665 .000570 .000247 The histograms of the error In estimating two of the four components of the observation noise covarlance are shown in Flg. 8.1,7 and 8.18,. The errors In estimating the covariance are plotted rather than the covariance estimates themselves because the actual (sample) covarlance of the noise varied from run to run due to the finite data length. As these histograms indi- cate, the error in estimating the 2nd component was consistently less than 5% for both the optimal and sub-optlmal inputs and most often within 1% for the optimal input. For the third component, the error with the optimal input is consistently less than 7% and with the sub-optlmal input, 9%. Overall, it is accurate to state that the performance of the optimal input in identifying the measurement noise covariances was only slightly improved over that of the sub- optimal input.

Optimal Input Suboptimal Input

z {p }- -2.68

Z (p} - 2.63 St'd. Vev. - .4069 SEtd Dey. 1 .3068 m E{p 2 }- 7.35 p2 } . 7.018

7-

m Real Value - -2.238 m

F-

I.-, I I

I i t j, i , !

i ]

I

_i I

II

; i

L , t

i !

.... i!

I!

I

- slr

II

-3.4 -3.0 .6 :2.2 -1.8 -3.4 -3 Fig. 8.1_ PARAMETER ESTIMATE HISTOGRAMS FOR M (Z Suboptimal Input Optimal Input.

E ( p } = -.0738 = -.122 St'd Dev. - .49!4 - .295 E( p2} . .2469 - .1019 Real Value = -.28 1.6 -.5 -.3 -.i .i .3 -.9 -.7 Fig. 8.1_ PARAMETER ESTIMATE HISTOGRAMS FOR M O_tlmal Input Suboptimal In ut

E{ p }= -.364

E { p } - -.397 St'd Dev. - .0092 St'd Dev. " .0235 E{_2} . .1328 _. {p2) . .1583

!

p-J ",4 Real Value - -.368

I

f

L

I

i i

il]

I

II ,_[

I Jl ! I __U_

-.39 -.37 -.35 -.45 -.43 -.41 -.39 -.37 -.35 -.33 Fig. 8.14 PARAMETER ESTIMATE HISTOGRAMS FOR Z Suboptimal Input -- . Optimal Input

E{p)- -.859

E {p) - -.933 St'd Dev. - .0537 St'd Bey. - .0771

p_ . .74l

E ( p2 } . .876

F-

F---

m Real Value

- -.829 l-

m F-

!

F-

--I

m m p u m -I.0 -.8 -1.2 -I.0 -.8 Fig. 8.]_ PARA)_TER ESTIMATE HISTOGRAMS FOR M 6 e .Optimal Input Suboptimal Input Z{ p} - .008Z9 E{ p }- .00732 St'd Dev. - ,00057 St'd Dev = .00168 E{P 2} = .000067 E{p 2} " .G00056 Real Value - .0076 m m

-i

L

t

,007 .009 ,003 .005 .007 .009 .010 Fig. 8.16 PARAMETER ESTImaTE HISTOGRAMS FOR Z 6 e .OptimaY Input Suboptimal Input I i !

_===_ I

t----

I !

I I

I--1-- _

i!

-.003 .001 .005 .009 .013 .007 .011 .015 -.005 -.001 .003 Fig. 8.17 HISTOGRAM OF ESTIMATION ERRORS FOR n 8 Optimal Input Suboptimal Input _4

I

_u "-4

_ IF-

.i f, k

i

N

M JI, l • II

-.003 •001 .005 .009 •013 -.017 -.013 -.009 -.005 -.001 .003 .007 .011 FiE. 8.18 HISTOGRAM OF ESTIMATION ERRORS FOR n q 8.3 Opt_m-I Input Through First-Order Filter One possible objection to the optimal inputs which are derived from the optimization technique is that they start from nonzero value which cannot be realized in practice. For example, in the case of the C-8 Monte Carlo simulation, the initial elevator deflection for the optimal input was I0 °. Although the suboptimal doublet input had an initial discontinuity of the same magnitude, a more realistic model of the situation would include the lag action of the control servos. This would prevent an instantaneous change in elevator deflection for both the optimal and suboptimal inputs. In addition, the optimal input, for identifying the short period parameters, has the disadvantage that the longitudinal velocity does not return to the nominal level after the input has been applied. For a nominal steady state value, the input must be two sided, i.e. both positive and negative values.

This section compares the performance of a two sided optimal input, with a total energy of 101.25 and length of 4 seconds to the performance of a suboptimal doublet input, of equal energy_ The 2 second optimal input which was used to generate the 4 - second input was derived for the C-8 four state longitudinal equations of motion given in section 8.2.1. A first order dynamical system with a time constant of .2 seconds and unity gain was used to simulate the control surface servo mechanism.

The 4 second optimal input after passed through the model of the servo dynamics, appears as in Fig. 8.19 . Note that the actual control surface deflection has no discontinuities, starts at 0 ° and has no violent maneuvers.

The suboptimal input, on the other hand, which is shown in Fig. 8w20 , while also beginning at 0 °, entails some much more drastic deflections.

Although a Monte Carlo indentlfication simulation was not performed for these two inputs, their performance, based on the charateristics of the information matrix, can still be compared. The determinant of the ELEVATOR DEFLECTION (DEG) .4 4.0 I Tr_(S_C) Fig. 8.19 TWO-SIDED oPTIMAL I1TPUT FROM OUTPUT OF FIRST ORDER SERVO v _D l--I CO ll4 c_ I-I I...4 £-1 !

A I-I !

CO ,,.4 !

I ! ! : " ,0 ,...4 I ! I ! ,-.4 ,,.4 ,-4 !

I I I Z I,.l I::1 information matrix, representing the reciprocal of the volume of the uncertainty ellipsoid, was 1.6 x 1026 for the optimal input and 7.3 x 1024 for the suboptimal doublet input. Three out of five eigenvalues of the information for the optimal input were larger, on a one-to-one basis, than for the suboptimal input, and the trace of M was 1.4 x 108 versus 9.7 x 106 , with the one for the optimal input again being larger.

In summary, the performance of the 4 second optimal input versus the suboptimal input including the effects of the making the input two-sided and the control surface servo delay, changed only slightly from that observed during the extensive Monte Carlo simulation. Further work along these lines is needed to make optimal inputs realizable in practice.

IX

CONCLUSIONS In this report, the problem of identifying aircraft stability and control derivatives from flight test data was discussed. It was shown that the three important elements of the identlflcatlon process are (i) the identification method, (li) the design of control inputs, and (ill) the instrumentation. The first part of the report described various methods that have been used in the past for identification. The special character- istics and limitations of these methods were discussed. A new technique, based on the Maximum Likelihood Criterion, was then derived and discussed in detail. The technique is applicable to nonlinear models containing both process and measurement noise. As such, the Maximum Likelihood Method presented here is the most general technique that has been developed for identifying stability and control derivatives from flight test data.

It consists of a combination of a Kalman filter (linear) or an Extended Kalman filter (nonlinear) for estimating the state and a Modified Newton- Raphson iterative procedure for estimating the parameters.

The Maximum Likelihood Identification Method was applied successfully to three different problems : (i) Identification of parameters from simulated data for a nonlinear model of X-22 VTOL aircraft containing both process and measure- ment noise.

(ii) HL-10 Lifting Body flight test data without gusts.

(ill) M2/F3 Lifting Body flight test data containing gusts.

The Output Error and Equation Error methods which were also tried gave either poor results or failed to converge on problems (i) and (ill).

Additional problems of identifiability were encountered in the processing of the HL-10 and M2/F3 flight test data. Those problems were manifested as a) opposite signs on the parameter estimates compared to wind-tunnel and theoretical values, b) large parameter covariances, and c) high correlation between the parameter estimates. Three different methods were tried for alleviating these problems: (i) A priori weighting (ii) Fixing certain parameters and (iil) Rank-deflclent solutions.

The last method seems to hold the greatest promise for further developments.

The second part of the report discusses the theory and the computation of optimal inputs for linear system identification. The criterion used for input design is the maximization of the output sensitivity to para- meter variations. This criterion is related to the Fisher Information matrix and the Cramer-Rao lower bound for the covariance of the parameter estimates. The objective of this effort was to determine the control surface deflections that enhance the effectiveness of identifying stability and control derivatives from flight test data. The specific application considered was the longitudinal dynamics of a C-8 aircraft. The optimal elevator surface deflection time histories were computed and compared with a doublet input of the same energy and time duration. A Monte Carlo simu- lation was done using the optimal and the doublet inputs and identifying short-period parameters from noisy output data. The optimal inputs were shown to provide more accurate estimates of aircraft parameters compared to a doublet input.

X AREAS FOR FURTHER INVESTIGATION The followlng areas are proposed for further research on and practical application of the maximum likelihood identification method and the design of optimal inputs: Complete preflight simulation of the Integrated Aircraft

(i)

Identification process for a future flight test program.

Identification of stability and control derivatives from (ii) flight test data for specific important (e.g., high angle-of- attack , transoni_flight regimes. This would require the development of methods for model structure determination since the aero- dynamic models in such flight regimes are not well known.

Correlation of the flight test results with the wind-tunnel

(iii)

results for high performance aircraft.

Determination of the effect of SAS on the identifiability

(iv)

of the parameters.

Further developments of the Rank-Deficient Method for solving

(v)

identifiability problems.

Extensions of the Input Design procedure to nonlinear models

(vi)

and process noise.

Modifications of the optimal inputs for pilot acceptability (vii) and ease of implementation.

Flight Test validation of optimal inputs.

(viii)

APPENDIX A EQUATION OF MOTION FOR X-22 VTOL The model for the equations of motion for the X-22 was designed to represent the longitudinal-vertical three degree of freedom motions of the vehicle in the body axis coordinate system. Defining u = change in velocity along the x axis from trim condition w = change in velocity along the z axis from turn condition 8 = change in pitch angle (fuselage attitude relative to the horizontal) from the trim value q = change in pitch rate from zero trim value; since vehicle is restricted to longitudinal-vertical, wings level motion only, q = = change in elevator stick deflection (positive _ gives es es positive pitching moment) from trim condition n x accelerometer signal along the x axis n = accelerometer signal along the z axis z the nonlinear model, derived from an examination of wind tunnel data, is given by m mi_-- X(u)w-qw-g sin 8 U Xo(u) X6es(U) w Zw(u)w + qu + gcos 8 Zo(u) Z,6es(U) + o o

(u)w +Mq (u)q

Mo(U)M_es(U) i m Xu(U) Xw(U) 0 " Vl z u (u) z (u) o V 2 + 0 0 0 V 3 m Mu(U Hw(U) Mq(U) The stability and control derivatives are all expressed as first order functions of u (e.g. M =A M + _ • u) except for Mo(U), Xo(U) and w w w o u Zo(U ) which are second order functions. The derivatives of Xu(U), Zud(u) d (Xo(U) + X(u) w + X6es(U) ; (u) and Mu(U) are given as du " 6es) du (Zo + Z (u)w + Z6es(u) and d Hw(u ) w q M_e s w " 6es) _u (Mo(U) + . + M (u) " q + (u) • _es ) are independent, zero mean white noise respectively and v i' v2' and v 3 sequences with (diagonal matrix) covariance Q. There are therefore a total of 23 parameters to identify, enumerated as H ; M u M _u! ; u o w ; q M6e s X o ; Xw ; X6 es and Z ; Z ; Z_e s u o 2 w There are seven measurements being made of the aircraft state. These include the four state variables u,w,8 and q , the pitch acceleration, q, and the two accelerometer outputs nv'"and n z . The equations for the accelerometer outputs are = 1 (u + qw) + sin e

g

.I (w- qu) - cos e n z = g Substituting in u and w, these expressions can be written as 0 u ..I- -Xo(U) X_e _ (u)] 0 w n z Zw(u)

_Z'o_(U) z6e_(u)_J

e q u m

X(u) "i

'u(U) ÷

u(u) z.w(u) v2

i v 3

APPENDIX B GRADIENT AND INFORMATION MATRIX CALCULATION WITH ADDED PARTIAL DERIVATIVE TERMS Since the presence of the process noise requires the use of a Kalman filter in computing the sensitivity functions, the exact equations for the gradient and information matrix will include the additional partials of the Kalman gain and state estimate error covariance with respect to the parameters. In those cases where the Kalman gain reaches a steady state, this steady state value can often be included as unknown parameters (wlth the only error being during the transient period). However, whenever the system equations are non-linear, as with the X-22 model, the Kalman filter will not reach a steady state and added partials will appear.

Using the notation of Section 5.1, the gradient and information computations as well as the sensitivity function computations are given below. The additional terms included in these calculations, due to the time varying Kalman gain, are indicated by a I i The term pj I _J" L-- -- is used to denote the jth entry in the vector of unknown parameters.

Gradient: ^ T _J _ I l i/i-i N E i/i-]

I 1 1

__ff _ xi/i- i _f'

i

_PJ i=l Yi _i/i-i _ 5p. + 3 3 ^ T I - xi/i- 1 -I _i/i-I

- _ _2 _ _/_-_

_P.

<i-xi,i_,)

where = H. H_ + R i/i-I i Pi/i-I l 5H.

H I +H i H_ + H.

_P. }Sp.

- Pi/i-I " BP L z Pi/i-i _p.

3 3

J

+ _P.

In£ormatfon Matrix: T -I _3^

_._:/°_^ _,_ +

i=l _ B _'Xili-i + Of' - B aXili-1 af' I °_il_-1 1 i=l Yl -i o O _ ili-I _I/i-i OPj I zi - _li-llT -i i=l @PK ^ B i _ + OP K Sermttivity Function: Between measurements: - + dt _p At a measurement:

(t:___ aP. ti - w i

ON i zi - kill_ I] i A

/ °xi/i-_. )

_f' OXl/i-i Of' --'_x - 0 5 - + O_j where

T

i/i-I Pi/i-i _ i/i-i H T -i O_i/i-1 -1

-Pi/i-i i _/i-i 05 _i/i-1

= (I + F i At) _Pi/i

_5

_ (I + r i At)T + _pj(hQgZ) At2 _t + _pj At Pi/i (I + F i At)T+ (I + Fi t_t) Pi/i \_--_j/ bPi/--i = (I - w i Hi) bPi/i-1 f_wi _H. \ _Pj apj apj and Fi "= _x_fl t=t; LIST OF REFERENCES i. C.H. Wolcowicz, "Considerations in the Determination of Stability and Control Derivatives and Dynamic Characteristics from Flight Data", AGARD Report 549, Part I 2. C.H. Wolcowlcz, K. W. lliff, G. B. Gilyard, "Flight Test Experience in Aircraft Parameter Identification", presented at AGARD Symposium on Stability and Control, Braunschwelg, Germany, April 1972 3. R. Turley, Private Communication, U.S. Air Force Flight Test Center, April 1971.

4. T.M. Kastner, J. A. Eney and J. J. McCue, "Flight Evaluation of Various Short Period Dynamics in the Variable Stability F-80", NATC TR FT-13R-TO.

5. H. Redless, Private Communication, NASA Flight Test Center, October 1971.

6. Dante DiFranco, "In-FlightParameter Identification by the Equation- of-Motlon Technique -- Application to the Variability Stability T-33 Airplane", Cornell Aeronautical Laboratory Report No. TC-1921-F-3, 15 December 1965.

7. O.H. Gerlach, "Determination of Performance and Stability Parameters from Non-Steady Flight Test Maneuvers", SAE National Business Air- craft Meeting, Wichita, Kansas, March 1970.

8. G.C. Goodwin, "Application of Curvature Methods to Parameter and State Estimation", Proc. IEEE, Vol. 16, No. 6, June 1969.

9. L. Taylor,et. al., "A Comparison of Newton-Raphson and Other Methods for Determining Stability Derivatives from Flight Data", Third Techni- cal Workshop on Dynamic Stability Problems, Ames Research Center, 1968. Also presented at AIAA Thlrd Flight Test, Simulation and Sup- port Conference, Houston, Texas, March 1969.

i0. R. Bellman, et. al., "Quasillnearization, System Identification, and Prediction, RAND Corporation RM-3812, August 1963.

ii. K.S.P. Kumar and R. Shridhar, "On the Identification of Control Systems by the Quasilinearization Method", IEEE Trans., Vol. AC-10, pp 151- 154, April 1964.

12. D. Larson, "Identification of Parameters by Method of Quasilinear- ization", CAL Report 164, May 1968.

13.

D. G. Denery, "An Identification AlgorlthmWhich is Insensitive to Initial Parameter Estimates", AIAA Eighth Aerospace Science Conference, January 1970.

14.

P.C. Young, "Process ParameterEstimatlon and Adaptive Control", In Theory of Self-Adaptive Control Systems, P. Hammond, ed., Plen_m Press, New York, 1966.

15.

R. D. Schalow, "Quasillnearizatlon and Parameter Estimation Accuracy", Ph.D. Thesis, Syracuse University, 1967.

16.

B. Dolbln, "A Differential Correction Method for the Identification of Airplane Parameters from Flight Test Data", University of Buffalo, Masters Thesis, December 1968.

17.

L. S. Lason, et. al., "The Conjugate Gradient Method for Optimal Control Problems", IEEE Trans., G-AL, Vol. 12, No. 2, April 1967.

18.

R. Wingrove, Private Communication, NASA Ames Research Center, October 1971.

19.

W.T. Suit, "Aerodynamic Parameters of the Navion Airplane Extracted from Flight Data", NASA TN D-6643, March 1972.

20.

G.S. Steinmetz, R.V. Parrish, and R.L. Bowles, "Longitudinal Stabil- ity and Control Derivatives of a Jet Fighter Airplane Extracted from Flight Test Data by Utilizing Maximum Likelihood Estimation", NASA TN D-6532, March 1972.

21.

J.S. Tyler, J.D. Powell, R. K. Mehra, "The use of Smoothing and Other Advanced Techniques for VTOL Aircraft Parameter Identification", Final Report to Cornell Aeronautical Laboratory under Naval Air Systems Command Contract No. N00019-69-C-0534, June 1970.

22.

K.J. Astrom and S. Wenmark, "Numerical Identification of Stationary Time Series", Sixth International Instruments and Measurements Congress, Stockholm, Sept. 1964.

23.

R.L. Kashyap, "Maximum Likelihood Identification of Stochastic Linear Dynamic Systems", IEEE Trans. AC, Feb. 1970.

24.

R.K. Mehra, "Identification of Stochastic Linear Dynamic Systems Using Kalman Filter Representation", AIAA Journal, Vol. 9, No. I, Jan. 1971 25.

A. E. Bryson and Y.C. Ho, Applied Optimal Control, Blalsdell Publishing Co., Waltham, Mass., 1969.

26.

R. E. Kalman, "A New Approach to Linear Filtering and Prediction Problems," Trans. ASME, J. Basic Eng. Vol. 82, March 1960.

27. T. Kailath, "A General Likelihood - Ratio Formula for Random Signals in Gausslan Nolse", IEEE Tzans. Info. Theory, Vol. IT-15, May 1969.

28. A. Papoulls, 'Probabillty_ RandomVarlableS andStochastlc Processes , McGraw-Hill, New York, 1965.

29. T. Kailath, "An Innovations Approach to Least-Squares Estimation, Part I", IEEE Trans. on AC, Vol AC-13, No. 6, Dec. 1968.

30. R. E. Larson, R. M. Dressler, and R.S. Rather, "Application of the Extended Kalman Filter to Balllstic Trajectory Estimation", Final Report, SRI, ProJ. 5188-103, Jan. 1967.

31. R. K. Mehra, "A Comparison of Several Non-llnear Filters for Radar Target Tracking", IEEE T-AC, Vol. AC-16, No. 4, Aug. 1971.

32. R.T. Chen, B.J. Eulrich and J.V. Lebacgz, "Development of Advanced Techniques for the Identification of V/STOL Aircraft Stability and Control Parameters", CAL Report No. BM-2820-F-I, Aug. 1971.

33. R.K. Mehra, '_aximum Likelihood Identification of Aircraft Parameters", 1970 JACC, Atlanta, Georgia R. Bellman, Introduction to Matrix Analysis, McGraw-Hill, New York, 1970 R. J. Hanson and C. L. Lawson, "Extensions and Applications of the Householder Algorithm for Solving Linear Least Squares Problems," Math. of Comp., October 1969.

36. D.G. Denery, "Simplification in the Computation of the Sensitivity Functions for Constant Coefficient Linear System", IEEE Trans. on AC, Aug. 71.

37. R.K. Mehra and R.E. Davis, "A Generalized Gradient Method for Optimal Control Problems with Inequality Constraints and Singular Arcs", SCI Research Report No. 2, Sept. 1970.

38. H. Akaike, "Statistical Predictor Identification", Ann. Inst. Statist.

Math., Vol. 22, 1970.

39. J. Sorensen, "Analysis of Instrumentation Error Effects on the Identi- fication Accuracy of Aircraft Parameters", Contract NAS 1-10791, May 1972.

NASA CR-I12121 40.

C. L. Muzzey and E.A. Kidd, "Measurement and Interpretation of Flight Test Data for Dynamic Stability and Control", Chapter ii, Vol II, AGARD Flight Test Manual, Pergammon Press, 1959.

41.

H. Greenberg, "A Survey of Methods for Determining Stability Parameters of an Airplane from Dynamic Flight Measurements", NASA TN 2340, April 1951.

42.

R.K. Mehra, "Optlmallnputs for Linear System Identification, Part I - Theory", JACC, Stanford, Californla 1972.

43.

P.A.N. Briggs, K.R. Godfrey and P.H. Hammond, "Estimation of Process Dynamic Characteristics by Correlation Methods Using Pseudo-Random Signals", I.F.A.C. Symposium on Identification and Process Parameter Estimation, Prague, June 1967.

44.

K.R. Godfrey, "The Appplicatlon of Pseudo-Random Sequences to Industrial Processes and Nuclear Power Plants", I.F.A.C. Symposium on Identification and Process Parameter Estimation, Prague, 1970.

45.

I.G. Cu_mlng, "Frequency of Input Signal in Identification", I.F.A.C.

Symposium on Identification and Process Parameter Estimation, Prague 1970.

46.

M. Aokl and R.M. Staley, "On Input Signal Synthesis in Parameter Identi- fication", Automatlca, Vol. 6, 1970.

47.

V.C. Levadi, "Design of Input Signal for Parameter Estimation," IEEE G-AC, Vol. AC-II, No. 2, April 1966.

48.

N.E. Nahi and D.E. Wallls, Jr., "Optimal Inputs for Parameter Estimation in Dynamic Systems with White Noise Observation Noise," Preprlnts, JACC Boulder, Colo. Aug. 1969.

49.

M.J. Levlne, "Estimation of a System Pulse Transfer Function in the Presence of Noise," IEEE AC, 1964, pp 229-235.

50.

R.J. McAulay, "Optimal Control Techniques Applied to PPM Signal Design", Information and Control 12, 1968, pp 221-235.

51.

R. Esposlto and M.A. Schumer, "Probing Linear Filters -- Signal Design for the Detection Problem", IEEE Trans Information Theory, Vol T-16, No 2, March 1970.

52. L.S. Pontryagin, V. Boltyanskil and E. Mishehenko, The Mathematical Theory of Optimal Processes, New York, Intersclence, 1962.

53.

A.E. Bryson and D.E. Johansen, "Linear Filtering for Time-Varying Systems Using Measurements Containing Colored Noise," IEEE Trans. AC Vol. AC-10 pp 4-10, January 1965.

54.

R.E. Kalman, "New Methods and Results in Linear Prediction and Filtering Theory," Proc. Symp. on Engineering Applications of Random Function Theory and Probability," New York, Wiley, 1961.

55.

R.E. Kalman and T.S. Enslar, "A User's Manual for the Automatic Synthesis Program," NASA CR-475, June 1966.

56. A. Ralston, A First Course inNumerical Analysis, McGraw-Hill, New York, 1965.

57. R. K. Mehra, "On-Line Identification of Linear Dynamic Systems with Applications to Kalman Filtering", IEEE Trans. Auto. Control, April 1970 U.S. GOVERNMENT PAINTING OFFLC£: ]_i _ - ?_ '_I', R_]]<N _. _-]_