UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
AN ANALYSIS OF METHODS FOR EXTRACTING AERODYNAMIC COEFFICIENTS FROM TEST DATA By DONALD CLIFTON DANIEL A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMNT OF THE REQUIRPF.IEENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY. UNIVERSITY OF FLORIDA 1973 To Pat. Digitized by the Internet Archive in 2010 with funding from University of Florida, George A. Smathers Libraries with support from Lyrasis and the Sloan Foundation http://www.archive.org/details/analysisofmethod00dani ACKNOWLEDGEMENTS Many people contribute in a variety of ways during the course of most graduate programs. I would like to express my sincere appreciation at this time to some of those who have helped me. First, I want to thank the members of my supervisory committee for the assistance and encouragement they have offered. In particular, this work would not have been possible without financial aid procurred from the Air Force Armament Laboratory by Dr. M. H. Clarkson. I also want to thank Dr. T. E. Bullock for the time he has spent teaching me to appreciate some aspects of modern control theory as well as the use and convenience of state variable notation. I would also like to acknowledge the teaching efforts of Dr. U. H. Kurzweg for his excellent courses in applied mathematics. Finally, Dr. W. W. Menke is due special thanks, not only for giving me a great deal of encouragement, but also for helping me to appreciate the importance of technological management. iii TABLE OF CONTENTS Page ACKNOWLEDGMENTS iii LIST OF TABLES V LIST OF FIGURES vii NOMENCLATURE LIST xii ABSTRACT xiv CHAPTER I. Introduction 1 II. The Coefficient Extraction Technique of Chapman and Kirk 7 III. Analysis of the ChapmanKirk Coefficient Extraction Technique 13 IV. Development of the Extended Kalman Filter for Estimating Parameters and Their Uncertainties 45 V. Use of the Extended Kalman Filter for Estimating Parameters and Their Uncertainties 60 VI. Conclusion 111 APPENDIX I. ChapmanKirk Coefficient Extraction Nomenclature List and Program Listing 119 II. Extended Kalman Filter Nomenclature List and Program Listing 142 REFERENCES 153 BIOGRAPHICAL SKETCH 158 LIST OF TABLES Page I. Initial Conditions and Nonhomogeneous Terms for Parametric Differential Equations 15 II. Program Options 19 III. True Values of Constants Used in Generating Data 22 IV. Summary of Results for Extracting Aero dynamic Coefficients from Dynamic Data Containing Random Measurement Errors 24 V. Extracted Coefficient Uncertainty Ratios 27 VI. Summary of Results for Extracting Aero dynamic Coefficients from Dynamic Data Containing Random Measurement Errors and System Noise 30 VII. Effects of Variations in the Initial Parameter Variances on Near Steady State Parameters and Their Variances (Linear System, Measurement Errors Only) 67 VIII. Effects of Variations in the Initial Parameter Variances on Near Steady State Parameters and Their Variances (Linear System, Measurement Errors and System Noise) 70 IX. Effects of Variations in the Initial Parameter Variances on Near Steady State Parameters and Their Variances (Nonlinear System, Measurement Errors Only) 77 X. Effects of Errors in the Estimate of the Measurement Error Variance on Near SteadyState Parameters and Their Variances (Nonlinear System, Measure ment Errors Only) 79 List of Tables, continued Page XI. Effects of Variations in the Initial Parameter Variances on Near Steady State Parameters and Their Variances (Nonlinear System, Measurement Errors and System Noise) 83 XII. Effects of Errors in the Estimates of Noise Variances on Near SteadyState Parameters and Their Variances (Nonlinear System, Measurement Errors and System Noise) 84 LIST OF FIGURES Page 1. Example of Program Fit to Data Containing Measurement Errors 34 2. Variation of Percent Error in Extracted C3 with Measurement Error 35 3. Variation of Percent Error in Extracted C4 with Measurement Error 35 4. Variation of Percent Error in Extracted C5 with Measurement Error 36 5. Variation of Percent Error in Extracted C6 with Measurement Error 36 6. Variation of RMS Residual with Measurement Error 37 7. Variation of Normalized Estimated Standard Deviation of C3 with Measurement Error 38 8. Variation of Normalized Estimated Standard Deviation of C4 with Measurement Error 38 9. Variation of Normalized Estimated Standard Deviation of C5 with Measurement Error 39 10. Variation of Normalized Estimated Standard Deviation of C6 with Measurement Error 39 11. Example of Program Fit to Data Containing System 'oise and Measurement Errors 40 12. Variation of Percent Error in Extracted C3 withI System Noise 41 13. Variation of Percent Error in Extracted C4 with System Noise 41 14. Variation of Percent Error in Extracted C5 with System Noise 42 vii List of Figures, continued Page 15. Variation of Percent Error in Extracted C6 with System Noise 42 16. Variation of Normalized Estimated Standard Deviation of C3 with System Noise 43 17. Variation of Normalized Estimated Standard Deviation of C4 with System Noise 43 18. Variation of Normalized Estimated Standard Deviation of C5 with System Noise 44 19. Variation of Normalized Estimated Standard Deviation of C6 with System Noise 44 20. Variation of Percent Error in X3 with Time (Linear System, Measurement Errors Only) 87 21. Variation of Percent Error in Xi with Time (Linear System, Measurement Errors Only) 87 22. Variation of Near SteadyState Parameter Error with Initial Parameter Variance (Linear System, Measurement Errors Only) 88 23. Variation of Normalized Near SteadyState Parameter Uncertainty with Initial Parameter Variance (Linear System Measurement Errors Only) 89 24. Variation of X3 Variance with Time for Three Initial Estimates (Linear System, Measurement Errors Only) 90 25. Variation of X4 Variance with Time for Three Initial Estimates (Linear System, Measurement Errors Only) 90 26. Variation of Percent Error in X3 with Time (Linear System, Measurement Errors and System Noise) 91 27. Variation of Percent Error in X4 with Time (Linear System, Measurement Errors and System Noise) 91 viii List of Figures, continued Page 28. Variation of Near SteadyState Parameter Error with Initial Parameter Variance (Linear System, Measurement Errors and System Noise) 92 29. Variation of Normalized Near SteadyState Parameter Uncertainty with Initial Parameter Variance (Linear System, Measurement Errors and System Noise) 93 30. Variation of Percent Error in X4 with Time for Two Estimates of the Initial Parameter Variance (Linear System, Measurement Errors and System Noise) 94 31. Variation of X3 Variance with time for Three Initial Estimates (Linear System, Measurement Errors and System Noise) 95 32. Variation of X4 Variance with Time for Three Unitial Estimates (Linear System, Measurement Errors and System Noise) 95 33. Variation of Percent Error in X4 with Time for Large Initial Parameter Estimate Error and Two Initial Parameter Variance Estimates (Linear System, Measurement Errors and System Noise) 96 34. Variation of Percent Error in X3 with Time (Nonlinear System, Measurement Errors Only) 97 35. Variation of Percent Error in Xi with Time (Nonlimear System, Measurement Errors Only) 97 36. Variation of Percent Error in Xs with Time (Nonlinear System, Measurement Errors Only) 98 37. Variation of Percent Error in X6 with Time (Nonlinear System, Measurement Errors Only) 98 38. Variation of Near SteadyState Parameter Error with Initial Parameter Variance (Nonlinear System, Measurement Errors Only) 99 List of Figures, continued Page 39. Variation of Normalized Near SteadyState Parameter Uncertainty with Initial Para meter Variance (Nonlinear System, Measurement Errors Only) 100 40. Variation of Near SteadyState Parameter Error with Error in the Estimate of the Measurement Error Variance (Nonlinear System, Measurement Errors Only) 101 41. Variation of Normalized Near SteadyState Parameter Uncertainty with Error in the Estimate of Measurement Error Variance (Nonlinear System, Measurement Errors Only) 102 42. Variation of Percent Error in X3 with Time for Large Error in Initial Estimate of X5 (Nonlinear System, Measurement Errors Only) 103 43. Variation of Percent Error in X4 with Time for Large Error in Initial Estimate of Xc (Nonlinear System, Measurement Errors Oniy) 103 44. Variation of Percent Error in Xs with Time for Large Error in Initial Estimate of Xc (Nonlinear System, Measurement Errors Only) 104 45. Variation of Percent Error in Xg with Time for Large Error in Initial Estimate of Xc (Nonlinear System, Measurement Errors Only) 104 46. Variation of Percent Error in X with Time (Nonlinear System, Measurement Errors and System Noise) 105 47. Variation of Percent Error in X with Time (Nonlinear System, Measurement Errors and System Noise) 105 48. Variation of Percent Error in X with Time (Nonlinear System, Measurement Errors and System Noise) 106 List of Figures, continued Page 49. Variation of Percent Error in X6 with Time (Nonlinear System, Measurement Errors and System Noise) 106 50. Variation of Near SteadyState Parameter Error with Initial Parameter Variance (Nonlinear System, Measurement Errors and System Noise) 107 51. Variation of Normalized Near SteadyState Parameter Uncertainty with Initial Parameter Variance (Nonlinear System, Measurement Errors and System Noise) 108 52. Variation of Near SteadyState Parameter Error with Error in the Estimate of the Measurement Error Variance (Nonlinear System, Measurement Errors and System Noise) 109 53. Variation of Normalized Near SteadyState Parameter Uncertainty with Error in the Estimate of Measurement Error Variance (Nonlinear System, Measurement Errors and System Noise) 110 NOMENCLATURE Symbol Definition A Reference area d Reference length I Vehicle moment of inertia about an axis through the center of gravity and normal to its pitch plane q Dynamic pressure, 1/2 pv2 V Freestream velocity C Static pitching moment coefficient moo derivative, rad1 C Static pitching moment coefficient ma2 derivative, rad3 C MStatic pitching moment coefficient Cmt4 derivative, rad5 C mqo Pitch damping coefficient mqo  C Pitch damping coefficient, rad2 C Static pitching moment coefficient at a=0 a Pitch angle p Freestream density B1 The inverse of the matrix B BT The transpose of the matrix B E[g(X)] Expected value of g(X) E[g(X)] = f. g(c)Fx(O)d where Fx () is the probability density function of the random variable X xii Symbol Definition 6.. Kronecker delta tt') Dirac delta function 6(tt') Dirac delta function xiii Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy AN ANALYSIS OF METHODS FOR EXTRACTING AERODYNAMIC COEFFICIENTS FROM TEST DATA By Donald Clifton Daniel March, 1973 Chairman: M. H. Clarkson Major Department: Engineering Science, Mechanics, and Aerospace Engineering An analysis of numerical methods for extracting aerodynamic coefficients from dynamic test data has been conducted. The emphasis of the analysis is on the effects that random measurement errors in the data and random disturbances in the system have on the accuracy with which the coefficients for linear and nonlinear systems can be determined. Both deterministic and stochastic methods for extracting the coefficients and determining their uncertain ties are considered. The deterministic technique considered, due to Chapman and Kirk, provides excellent estimates of both linear and nonlinear static pitching moment coefficients for the range xiv of measurement errors and system noise considered. Somewhat less accurate estimates of linear damping coefficients are obtained. Nonlinear pitch damping coefficients extracted using this deterministic technique are affected considerably by both measurement errors and system noise. The estimated standard deviations of the extracted coefficients obtained using standard techniques are generally adequate when the data being analyzed contain only measurement errors. The stochastic approach considered demonstrates the feasibility of using an extended Kalman filter, with a parameter augmented state vector, for determining the values of the aerodynamic coefficients and their uncertainties from noisy dynamic test data. The specific filter used generally reaches near steady state conditions in its estimates of the parameters in less than one second. Variations in the initial parameter variances or in the estimates of the noise statistics essentially affect only the determination of the nonlinear damping parameter. Parameter estimates obtained with the extended filter compare favorably with previously obtained results using deterministic techniques. Estimates of the parameter uncertainties provided by the filter are generally superior to those obtained with deterministic techniques particularly when system noise has corrupted the data. CHAPTER I INTRODUCTION This dissertation presents an analysis of methods for determining aerodynamic coefficients of flight vehicles from dynamic test data. This determination is accomplished by finding numerical values of the aerodynamic coefficients appearing in the equations of motion such that the solutions to these equations are adequate representations of the given test data. When this determination has been made, the coefficients are said to have been extracted from the dynamic data. The emphasis of the analysis is placed on the effect that random measurement errors in the data and random disturbances in the system have on the accuracy with which the coefficients for linear and nonlinear systems can be determined. Both deterministic and stochastic system models for extracting the coefficients and determining their uncertainties are considered. Historical Background The determination of aerodynamic coefficients from dynamic data is generally agreed to have begun with the work of Fowler, Gallop, Lock and Richmond1 in the first quarter of this century. Their basic technique was con cerned with determining moment and damping characteristics of artillery shells by firing these projectiles through spaced cards and reconstructing the pitch and yaw angle time histories by observing the obliqueness of the holes that resulted when the, shells passed through the cards. Their technique of data measurement is still in use at some ballistic ranges to this date.2 Nielsen and Synge3 later clarified the linear theory of Fowler et al. in their work during, and immediately following, World War II. Coefficient extraction techniques that are currently in use at many ballistic and wind tunnel facilities evolved from the work of Murphy4'5 and Nicolaides,6 both of whom have considered various combinations of degrees of freedom for a variety of flight vehicles. Their aerodynamic co efficient extraction techniques feature a least squares fit to dynamic data of the exact solutions of linear equations of motion, or approximate closed form solutions of slightly nonlinear equations of motion using the quasi linearization technique of Kryloff and Bogoliuboff.7 Ex tensions of the work of Murphy and Nicolaides, particularly for more complex nonlinear systems, have been made by Eikenberry and Ingram. The requirement for closed form solutions of the equations of motion has recently been eliminated by the formulation of least squares techniques which fit numerical solutions of the equations of motion to dynamic data. The minimization of the least squares criterion function involves differential corrections, as do the techniques of Murphy and Nicolaides; the required partial derivatives are determined by numerically integrating parametric differ ential equations which are derived from the equations of motion. Contributions to this numerical coefficient extraction technique have been made independently by Chapman and Kirk,10 Knadler,11 Goodman12'13 and 14 Meissinger, although the method is usually referred to as the ChapmanKirk technique. The computational require ments of this technique are sometimes extensive,15 but this is usually outweighed by the fact that it can be used to analyze highly nonlinear aerodynamic forces and moments.16 All of the coefficient extraction techniques that have been discussed to this point are deterministic in nature in that the modeling of the equations of motion does not account for random disturbances in the system. Most angular or translational motion data obtained from ballistic rac:ges or wind tunnel dynamic tests have, never theless, been affected by random system noise and random measurement errors. The effects of noise of these types on the accuracy of coefficients extracted from dynamic data using deterministic techniques are of current interest.17 With the exception of Ingramn's partial analysis of noise 8 effects on Eikenberry's "Wobble" program, however, little has been reported concerning these effects. One of the primary goals of the research reported here is to show some effects of random measurement errors and system noise on the accuracy of aerodynamic coefficients extracted from dynamic data using the most general of the various determi nistic techniques, that of Chapman and Kirk. In recent years increasing attention has been given to parameter and state variable estimation through the use of stochastic modeling of the physical system of interest. Most of this work has been done by optimal control specialists and is a direct result of pioneering contributions in linear filtering theory of Kalman 8'19 19 and Bucy. The Kalman filter provides estimates of the states of noisy, linear dynamic systems as well as estimates of the state variable uncertainties. Bryson and Ho,20 among others, have extended the Kalman filter theory to include estimates of the states for nonlinear systems. Mehra21 has recently proposed a maximum likelihood technique for the determination of aerodynamic coefficients from dynamic data using the Kalman filter for linear systems 5 and the extended Kalman filter when the system is non linear. Results verifying his proposals are, as yet, unavailable. An additional maximum likelihood technique has been developed by Grove et al.22 and has been used with modest success by Suit.23 A stochastic approach to the coefficient extraction problem using an extended Kalman filter with parameter augmented state vector is developed in this dissertation. Results obtained with this filter are also presented. Scope The following chapters are devoted to an analysis of the effects of random system noise and measurement errors on coefficients extracted from dynamic data using the method of Chapman and Kirk, as well as the development and evaluation of an extended Kalman filter for solving essentially the same problem. The equation of motion used in the analysis is one describing the onedegree offreedom pitching motion of a rigid body. Pitching moments that are both linear and nonlinear functions of angleofattack are considered. The theory of the ChapmanKirk technique is recounted briefly in Chapter II, followed by the analysis of noise effects on this method in Chapter III. 6 Developments leading to the formulation of the extended Kalman filter to be used are given in Chapter IV. Results obtained with the extended filter are presented and discussed in Chapter V. Some concluding remarks are offered in Chapter VI. Computer program listings and definitions of program variables are given in the two appendices following the body of the paper. CHAPTER II THE COEFFICIENT EXTRACTION TECHNIQUE OF CHAPMAN AND KIRK A brief reconstruction of the ChapmanKirk algorithm for extracting aerodynamic coefficients from dynamic data is given in this chapter. The standard technique for estimating uncertainties in coefficients determined from a least squares fit of a given function to experimental data is also presented. Development of the Extraction Algorithm The extraction technique of Chapman and Kirk has two very basic requirements: (1) the different al equation of motion for the body of interest must be given and (2) a set of experimental data based on the observed motion of the body must be available. As an example, consider the nonlinear equation of pitching motion for a rigid body a + (C4 + CGa2)& + (C3 + C5a2)a = 0 (1) subject to the initial conditions a(0) = C1 & (0) = C2 where a indicates the derivative of a with respect to time. Suppose, also, that the time history of the pitch angle as recorded during some experiment, ae(ti), i=l,2,...m, is also available. The technique of Chapman and Kirk is used to determine the values of the C. (j=1,2,...,6) in Equation (1) which result in the solution to this equation of motion being a best fit to the test data in a least squares sense. Thus it is necessary to minimize the least squares cri terion function m 2 S = ae(ti) ac(ti) (2) i=l where ac(ti) is obtained from the solution to Equation (1). Now it is well known that in order to determine para meters directly by a least squares fit of a given function to test data, the parameters must appear linearly in the function. This requirement is met in the problem at hand by expanding aC in a truncated Taylor series about the numerical solution resulting from some initial estimates of the parameters of interest, Cjo, i.e., "c(ti) = aco(ti) + AC. (3) j=l 3ji " where Da are evaluated for Cj=C. and AC.=C.C acjj 30 o J 3jo Substituting (3) into (2) yields M 6 2 S = e (t) a (ti) ( a AC. (4) i= co1 i j= aCJ 3 Now, assuming that the, [a are known, Equation (4) is a function of the AC.'s only. Therefore, to determine the values of these coefficients that will minimize this equation, it is necessary to take the partial derivative of Equation (4) with respect to each AC., set each of the resulting equations to zero and solve for the AC.'s. Carrying out these operations yields the following matrix equation [A] [AC] = [B] (5) where A is, in this case, a 6 X 6 matrix with elements given by A C J Ci (6) jk Xi=l ac i Ck AC is a 6 X 1 column matrix, or vector, and B is a 6 X 1 column matrix with elements given by m ( B. = a (ti) ac t.i) 1, (7) 3 il i co I i The solution to Equation (5) is [AC] = [A]1 [B] .(8) The solutions for the AC.'s obtained from Equation (8) are exactly correct only if a is a linear function of the C.'s as assumed by Equation (3). This condition of linearity is seldom the case and the process must be repeated with new initial guesses, C = C. + ACo (9) until the change in the criterion function [Equation (2)] from one iteration to the next is sufficiently small. The algorithm just presented requires that time histories of the influence coefficients, a, be available. These time histories are determined by numerically inte grating parametric differential equations which are derived by differentiating the equation of motion with respect to each of the parameters of interest. As an example, the parametric differential equation for C1, the initial condition of a is L[a + (C4+C6a2)& + (C3+Csa2)a] = 0 aLn as aCs 3Cs Ba a 2Cs 3a + aC4 + 2C6 a2& + 2Ca& a + C6a2 3& + 3C a C1 aC, aC 6 ac, aC1 aC1 3a 3C5 aa + C3 a +  a3 + 3C5d2 0 3R+1 a CI Assuming that the parameters are independent of each other and that the order of differentiation can be reversed, the final form of the desired equation is d2 fa d fa at dt2 C a2 d 1 + (C +3C a2+2C a&) J 0 (10) dt2 to th dtinitial cons subject to the initial conditions aa(0) aC, aC aC1 ac(O) aC2 aC1 aC1 The complete set of parametric differential equations for the given equation of motion [Equation (1)] is given in Chapter III. In summary, the process for extracting numerical coefficients from test data given the system model [Equa tion (1)] and criterion function [Equation (2)] is 1. Estimate the numerical values of the C.'s. 2. Integrate Equation (1) to obtain aco(ti). (11) 3. Determine L 4. Solve Equation (8) for the AC 's. 5. Repeat the process with C. =C. +AC. until the 31 Jo Jo change in Equation (2) is sufficiently small. Estimation of Extracted Parameter Uncertainties The estimation of the uncertainties, or standard deviations, of parameters that have been determined by the least squares fitting of a given function to test data is a wellknown result available in a variety of references (see, for example, References 24 or 25). The least square parameter uncertainties are estimated in this paper by A.. S (12) mK j jj K (12) where A. is the jth diagonal element of the inverse Grammian matrix [Equations (6) and (8)], S is the sum of the squares of the residuals as given by Equation (2), m is the total number of data points and K is the total number of parameters being determined by the fit. CHAPTER III ANALYSIS OF THE CHAPMANKIRK COEFFICIENT EXTRACTION TECHNIQUE Coefficient Extraction Computer Program This section provides a detailed description of a onedegreeoffreedom coefficient extraction computer program based on the previously described iterative process of Chapman and Kirk. This program considers the pitching motion of a symmetric missile about a fixed point. It is used to determine values of static pitching moment coefficient derivatives, pitch damping coefficients, and a trim term so that the solution to the appropriate differ ential equation of motion is a best fit to test data in a least squares sense. The program was developed to be used as an economical tool in the evaluation of the sensitivity to noise and convergence sensitivity of the ChapmanKirk technique when operating in its least complex mode. Computational Equations Equation of motion The complete equation of motion that is used in the program is 14 S+ (C +C a2)& + (C3+C a2+C7 a)a + C = 0 with initial conditions caO) = C1 where (Cmao)q 3 I &(0) = C2 Ad C (Cmqo)qAd2 S2VI 2VI C (Cma2)qAd sI S(Cma4)qAd C7  (Cmq2)qAd2 6 2VI (CM6 )qAd 8 I Parametric differential equations The eight parametric differential equations for the given equation of motion [Equation (13)] and initial con ditions [Equation (14)] are of the form P. + APj + BP. = F. 3 j = 1,2,...,8 (15) 3 a3C 3j S d 3a ' J dt 3Cj 3 d2 'ca 3 dt32 9 A = (C4 + C6a2) (13) (14) where and B = (C3 + 3C5a2 + 2C6a& + 5C7a4). The values of the initial conditions P.(0) and P (0) as J J well as the functional forms of the nonhomogeneous term F. are give in Table I for j=1,2,...,8. Table I. INITIAL CONDITIONS AND NONHOMOGENEOUS TERMS FOR PARAMETRIC DIFFERENTIAL EQUATIONS Pj (0) P (0) 0 Ci 3 a a5 1 Program Description The program is written in Fortran IV for use primarily on an IBM 360/65. The paragraphs that follow describe the functions of the main program and its four subroutines, the required input data, and the program options. A listing of the complete program is given in Appendix I. Main program functions The functions of the main program in their approximate order of occurrence are as follows: (a) Reads and writes input data. (b) Computes initial parameter values (these are the C.'s that appear in Equation (13). (c) Calls the numerical integration subroutine ADDUM. (d) Writes current parameter values. (e) Computes the matrix elements required for incre menting the parameter values. (f) Calls the matrix inversion subroutine MINV. (g) Computes the sum of the squares of the residuals between the calculated and experimental data points. (h) Computes the rootmeansquare residual and root meansquare error (Reference 24) of the current fit to date. (i) Computes the estimated standard deviations (Reference 24) of the current parameter values. (j) Tests the difference between the current and previous values of the rootmeansquare error to determine if the iteration process has converged. (k) Computes the incremental changes for the para meters (if convergence has not occurred). (1) Computes the updated parameter values (if convergence has not occurred). (m) Returns to (c) (if convergence has not occurred). (n) Computes the values of the extracted coefficients from the current parameter values (after con vergence is achieved). (o) Computes the estimated standard deviations of the extracted aerodynamic coefficients (after convergence is achieved). (p) Writes extracted aerodynamic coefficients, estimated standard deviations, and the pitch angle output that represents the final fit to the experimental data (after convergence is achieved). Subroutine functions The names and functions of each of the four subroutines that are used in the coefficient extraction program are given below. ADDUM.Subroutine ADDUM integrates the equation of motion and parametric differential equations using a fourth order RungeKutta method for starting and a fourth order AdamsBashforth predictorcorrector method for running. It calls subroutines XDOT and OUT. This sub routine is described in detail in Reference 26. XDOT.Subroutine XDOT computes current values of first derivatives that are required by ADDUM. OUT.Subroutine OUT stores the results of the numerically integrated equation of motion [Equation (13)] and parametric differential equations [Equation (15)]. MINV.Subroutine MINV inverts an K x K matrix using a standard GaussJordan technique and is described in detail in Reference 27. Required input data The program reads six categories of input data. These categories and the specific data in each are de lineated in Appendix I. The format and units of the entries on a specific data card can be determined from the program listing and nomenclature list provided in this appendix. Program options This program has options for extracting various combinations of aerodynamic coefficients from the given test data, in addition to the option of extracting no coefficients and merely integrating the equation of motion. These options are controlled by the numerical value of the number of first order ,equations to be integrated, N, which is read by the program on the first data card. The various options are given in Table II. Table II. PROGRAM OPTIONS N Coefficients to be Extracted 2 No coefficients are extracted. The program integrates the equation of motion for the given input data and prints the results. 8 ao, &,o Cmao 10 ao, o', Cmao, Cmqo 12 ao, &o, Cmao, Cmqo, Cma2 14 ao, ao, Cmao, Cmqo, Cma2, Cmq2 16 ao &0o Cmao' Cmqo' Cma2' Cmq2' Cma4 18 ao, &o, Cmao' Cmqo' Cma2' Cmq2 Cma4' Cm6a .1.alysis of the Sensitivity to Noisej of the ChlapmanKirl Technique The evaluation of the sensitivity to noise of the ChapmanKirk coefficient extraction technique has been approached from several directions. Coefficients have been extracted from numerous sets of computer program generated data containing only measurement errors as well as data containing system noise and measurement errors. In addition to the use of artificially generated data, some very preliminary work has been done with noisy wind tunnel data obtained from actual experimentation. The design of the experiment used to produce the dynamic data is described by Turner in Reference 28 along with some preliminary results obtained with the onedegreeoffreedom program described earlier in this chapter. Since the wind tunnel experimentation is still in a developmental stage, however, the results presented in the remainder of this chapter are those obtained from arti ficially generated data. Generation of Noisy Dynamic Data The computer program UFNOISE described in detail in Reference 29 was used to generate the dynamic data from which the aerodynamic coefficients were extracted. This program simulates the pitching motion of symmetric missile oscillating in a wind stream about a fixed point, for any given initial pitch angle displacement and initial pitch angle rate, by numerically integrating the equation of motion. The program has two options for simulating system noise: it considers the magnitude and direction of the freestream velocity vector as normally distributed random variables, with programmer set mean values and standard deviations, to determine the system noise perturbation accelerations or it simply selects random perturbation accelerations which have zero mean, are normally distributed and have a programmer set standard deviation. Both methods are essentially equivalent. New values of the system noise perturbations are randomly selected at each numerical integration step. The program also has the option of simulating random measurement errors of the pitch angle by superimposing normally distributed random noise on the output of the numerical integration. The basic equation of motion used to generate the noisy dynamic data was a slightly simplified form of Equation (13): a + (C4 + C6a2)& + (C3 + C5a2)a = w(t) (16) subject to a(0) = Ci (O) = C2 (17) The true values of the physical and aerodynamic con stants used in generating the data are given in Table III. The various noise level standard deviations and mean values are discussed in the following paragraphs to which they are pertinent. Table III. TRUE VALUES OF CONSTANTS USED IN GENERATING DATA Constant Value d(Ft) A(Ft2) I(slu.g ft2) V (ft/sec) m q(lb/ft2) Cmoo(rad 1) C (rad 3) ma2 C mqo Cmq(rad2) a0 (rad) a (rad/sec) 0 0.333 0.0873 0.1080 500.0 297.0 2.00 24.5 60.0 163.0 0.5235 0.0 Effects of Random Measurement Errors Aerodynamic coefficients were extracted from a total of nine sets of data containing only measurement errors. Each data set was made up of 201 discrete points which re sult from integrating the equation of motion [Equation (16)] in increments of 0.005 second for a total of 1.000 second. The desired mean value of the measurement errors for each of the mine data sets was 0.0 radian. The desired standard deviations of the random errors were 0.00146 radian, 0.00582 radian, and 0.01745 radian,with three different sets of data being generated at each noise level. These standard deviations are typical of measurement un certainty for data of this type (Reference 9). The above standard deviations correspond to percent noise levels of 1.1, 4.4, and 13.2, respectively where percent noise level is defined in Reference 30 as Percent Noise = (18) Approximate average peak amplitude with a being the standard deviation of interest. The approximate average peak amplitude can be determined in a variety of ways. The value used here is the mean positive peak amplitude for the three cycles that result when the nominal equation of motion is integrated for 1.000 second. The results of this portion of the analysis are given in Table IV and on Figures 1 through 10. Table IV is primarily a summary of the percent error in the extracted aerodynamic coefficients together with their normalized estimated standard deviations for the various actual measurement errors, om. The percent error in the extracted coefficients is defined by Percent Error (C) = (CC) 100 (19) J r rf Ci : J Id t, I Cd " E Z 0 0 U J M a0 U E Z U Sw ;2Fo P4d u < u O M U O m oltfr I O O 9*# V )i *iJ 0) 4) S.rt 4 ,.4 0 i u d a3 u n tJ OV) md X c 0 rO S$ a z :: O .0 ro co 0) %O tn '0 a 0 0 N O 0 f 4 r c o I V, r 'd  in N N * 0 0 * * r a 0 0 0 0 oo oo 0 0 0 0 * ra 0 r4 0 0 * *0 0 0 CD' C) iCl C) * L/ *= C14 OC cM cn 0 0 0 0 0 0 C14 1 l il N t ) 00 * OH . 4 4  C4 4q 0 r4 c M v v , 0 a cc 0 N . v * * 0 CM 0 0 0 0 0 0 C> C C) C oi 0 O. fr 4) cd U) m where C. is the coefficient estimate determined with the 3 program in the fitting of the data and C. is the true value of the coefficient of interest. The normalized estimated standard deviations are the estimated standard deviations calculated from Equation (12) divided by the true value of the parameter of interest, Percent = 100 (20) Jn Cj Table IV also contains the RMS residual for each fit to the noisy data and the percent noise levels based on the values of am. An ex3 ple of one of the program fits to a noisy data set is shown on Figure 1. Figures 2 through 5 show the variation of the percent error in the individual extracted coefficients with measurement error. The variation of RMS residual with measurement error is shown on Figure 6. The variations of the normalized estimated standard deviations of the extracted coefficients with measurement error are shown on Figures 7 through 10. An analysis of Table IV and Figures 2 through 10 reveals the following. (1) The extracted static pitching moment coefficient derivatives show little or no error for the entire range of measurement errors considered. ( 'e Figures 2 and 3.) (2) The extracted pitch damping coefficients show some error for the lower noise values (om 0.005 radian) and deviate significantly for am=0.018 radian. (See Figures 4 and 5.) (3) The variation of RMS residual with measurement noise is essentially linear with a slope of unity (Figure 6). ,This indicates that the RMS residual, as calculated by the coefficient extraction pro gram, is a good estimate of the amount of measure met noise in the data when this is the only type of noise present. (4) All variations of the normalized estimated standard deviations of the extracted coefficients, ajn. in increase linearly with measurement error (Figures 7 through 10). Approximate values of the slopes of these variations are given in Table V. These normalized uncertainty ratios give the relative (to each other) uncertainty which can be expected when extracting coefficients from data containing measurement noise using the given coefficient extraction program. (5) The true value of a given coefficient is contained within the interval difined by its extracted value 3ac for every coefficient extracted from the nine sets of data. 27 Table V. EXTRACTED COEFFICIENT UNCERTAINTY RATIOS Coefficient (Cj) Approximate Normalized Uncertainty Ratio (Aj /A m) Effects of Random Measurement Errors and System Noise For this portion of the analysis, attempts were made to extract aerodynamic coefficients from twelve basic sets of dynaunic data. Each data set was made up of 201 discrete points resulting from integrating the equation of motion ina increments of 0.005 second for 4.00 seconds; the pitch augle was output every four integration steps so that the time between data points was 0.020 second. The time between data points used in this portion of the analysis is different from that previously used. This change was made so that the time between data points would correspond to that used by Turner28 in actual experiment ation. Of the twelve basic data sets, nine contain system noise and measurement errors, whereas three contain only measurement errors. Coefficients were extracted from these latter three to determine if the abovementioned change in the time increment between data points had any appreciable effect on the extracted coefficients or their uncertainties. The estimates of the uncertainties were improved due pri marily to the fact that a longer data record was being analyzed. The desired mean values of the freestream velocity vector direction fluctuations and measurement errors for all data sets were zero. The desired standard deviation of the measurement errors, am, was 0.00582 radian (30m = 1.0 degree) for all sets; the desired standard deviation of the velocity magnitude fluctuations, av, was 5.0 ft/sec for all nine data sets containing system noise. The desired standard deviations of the velocity direction, oa, were 0.02910 radian (3 a = 5 degrees), 0.05820 radian (3aa = 10 degrees), and 0.11640 radian (3oa = 20 degrees), with three different sets of data being generated at each noise level. The random velocity fluctuations act as a forcing function for the equation of motion and cause oscillations even if the vehicle has no initial displacement. The resultant maximum amplitudes of these forced oscillations are usually within a certain magnitude or noise band width. The widths of the noise bands for the data used in this analysis are generally equivalent to the 3aa values, and these have been used to calculate the percent system noise values. The approximate average peak amplitude used in the percent noise calculations is the mean positive peak ampli tude for the 10 cycles that result when the nominal equation of motion is integrated for approximately four seconds. The results of this portion of the analysis are given in Table VI and on Figures 11 through 19. Table VI is a summary of the percent error in the extracted aerodynamic coefficients and their normalized estimated standard deviations along with the various noise Table VI. SUMMARY OF RESULTS FOR EXTRACTING AERODYNAMIC COEFFICIENTS FROM DYNAMIC DATA CONTAINING RANDOM MEASUREMENT ERRORS AND SYSTEM NOISE a (Rad) 0.00 0.00 0.00 0.02L . 0.03038 0.02828 0.005516 0.06120 0.05662 0.11262 0.11706 0.12058 Percent System Noise 0.00 0.00 0.00 41.4 43.1 40.2 78.4 86.8 80.3 152.9 158.9 163.7 Data Set a m (Rad) 0.30527 0.00589 0.00574 0.00611 0.00597 0.00587 0.00603 0.00562 0.00600 0.00572 0.00547 0.00571 ov (Ft/sec) 0.00 0.00 0.00 5.05 5.18 4.88 5.01 4.90 4.96 5.08 4.85 4.96 ~ (Extended) Percent Error in Extracted Coefficients and Normalized Estimated Standard Deviations C3(03) C4(04) C5(05) C (as) RMS Residual (Rad) 0.27(0.25) 0.34(0.28) 0.42(0.27) 0.55(0.64) 4.19(0.63) 2.87(0.55) 2.71(1.24) 14.44(1.32) 1.41(0.70) 0.84(0.93) 0.29(1.04) 0.36(1.01) 8.25(2.75) 5.84(2.41) 20.73(2.10) 42.80(3.95) 0.90(5.35) 35.67(2.33) 0.68(0.65) 1.51(0.73) 0.94(0.69) 9.55(1.88) 1.35(1.63) 2.80(1.24) 6.54(3.64) 5.44(3.87) 14.59(2.30) DIVERGED DIVERGED DIVERGED 5.72(12.92) 3.60(14.43) 4.40(13.88) 329.59(43.58) 44.80(34.40) 247.20(26.39) 605.77(66.79) 197.10(90.61) 588.34(44.19) 0.00520 0.00574 0.00563 0.01258 0.01244 0.01026 0.02806 0.02627 0.01488 Table VI. level standard deviations, percent system noise levels, and the RMS residual of each fit to a noisy data set. The percent system noise levels are similar to those encountered experimentally by Turner in Reference 28. Figure: 11 is an example of one of the program fits to a noisy data set. Figures 12 through 15 show the varia tion of percent error for the extracted coefficients with system noise. The variations of the normalized estimated standard deviations of the extracted coefficients with system noise level are shown on Figures 16 through 19. An analysis of Table VI and Figures 12 through 19 reveals the following: (1) The errors in the linear static pitching moment coefficient derivative parameter, C3, are gen erally less than 5 percent for the entire range of system noise considered. The errors in the corresponding nonlinear term, C5, are generally less than 10 percent. (See Figures 12 and 13.) (2) ThMe pitch damping parameters show significant error for all nonzero system noise levels. The e:xtracted values of the nonlinear term, C6, are sometimes in error by several orders of magnitude at the higher noise levels. (See Figures 14 ard 15.) (3) The estimated standard deviations of the extracted coefficients generally increase as the system noise increases. The estimated standard deviations of coefficients extracted from dynamic data containing both system noise and measurement errors are generally too small and do not reflect the true uncertainty of the extracted coeffi cients as was true in the previous case when the data contained only measurement errors. (See Table VI and Figures 16 through 19.) (4) All attempts at extracting the complete set of coefficients from data with a system noise band of approximately 20 degrees failed (see Table VI). The failures resulted when a value of C6 was eventually calculated by the iterative process which caused the solution to the equation of motion to diverge. (5) The divergence problem can sometimes be circum vented by not attempting to extract Cm from extremely noisy data. The resulting coefficients that are extracted, Cmo, Cm2, and Cmqo have accuracies comparable to those extracted from noisy data with a 10 degrees system noise band. N4 N I + to It Ia II E E 0 C. 00 4) *4 c+ 0 0 0 0 I I 30 200 10 0 Figure 2. 0.005 0.010 0.010 om(Rad) Variation of Percent Error with Measurement Error I 0.015 000 0.020 in Extracted C3 30 20  10 Figure 3. 0 0.605 0.i10 om(Rad) Variation of Percent Error with Measurement Error o.d1s 0 0.020 in Extracted C4 30 20 10 0 Figure 4. 0 100 200 300 Figure 5. _~h rO0 0.010 om(Rad) 0.015 Variation of Percent Error in Extracted C5 with Measurement Error 0.005 I 0.010 0.015 I 0.020 o (Rad) .m Variation of Percent Error in Extracted C6 with Measurement Error 0 r Q Oi o.oosO 0.0050 0.020 0.015  0.010  0.005  0.000 GP 0.005 0.005 0.010 0.015 0.015 0.020 a (Rad) Variation of RMS Residual with Measurement Error I __ I __ _ Figure 6. C, 15  10 Figure 7. 15  10  5 0.005 0.005 I 0.010 m (Rad) 0.015 0.015 0.020 0.020 i I E i da Variation of Normalized Estimated Standard Deviation of C3 with Measurement Error o@ 0.005 0.005 0.010 0.010 0.015 0.015 0.020 Figure 8. Variation of Deviation of om(Rad) Normalized Estimated Standard C4 with Measurement Error 15 10 Figure 9. 120  80 40 0.005 0.005 0.010 om(Rad) 0.015 0.020 0 0 Varitio ofNraie stmtdSadr Variation of Normalized Estimated Standard Deviation of C5 with Measurement Error O 0.005 0.005 0.010 0.015 0.020 Figure 10. m (Percent) Variation of Normalized Estimated Standard Deviation of C6 with Measurement Error cq 0 S4) sU r4 '4 0r '4 (p~U) 70 c0 t N9 00 S N O cq C> to I 0 II U i 0 U 0 41 0 o. & cM O 0 ' C9i I 0.04 a (Rad) a Figure 12. 10 10 20 30 40 Figure 13. Variation of Percent Error in Extracted C3 with System Noise 0.02 0.04 0.06 0.08 o (Rad) a 0 Variation of Percent Error in Extracted C4 with System Noise 20 10 0 0 0.02 0.02 I 0.06 1 0.08 e p 20 10 Figure 14. 600  S0 0 in 0 4J U Pc) 0.04 0.06 0.08 o (Rad) Variation of Percent Error in Extracted C5 with System Noise 400 I 200  0 200  0.02 Q 0.04 0.06 o.b8 a (Rad) a Figure 15. Variation of Percent Error in Extracted CG with System Noise 0.02 0.02 ci) ci) 1 _ Y 0.02 Figure 16. 10 4J U r &S 0.04 0.06 0.08 a (Rad) Variation of Normalized Estimated Standard Deviation of C3 with System Noise O 0.02 0.04 0.06 0.08 Figure 17. a (Rad) Variation of Normalized Estimated Deviation of C4 with System Noise Standard L.    . 1 7 10  0 Figure 18. 80 40  0.02 Variation Deviation 0.04 o (Rad) a 0.06 0.06 of Normalized Estimated of C5 with System Noise 0.08 Standard O I 0.02 0.04 0.06 I 0.08 Figure 19. a (Rad) Variation of Normalized Estimated Deviation of C6 with System Noise Standard _ __ CHAPTER IV DEVELOPMENT OF THE EXTENDED KALMAN FILTER FOR ESTIMATING PARAMETERS AND THEIR UNCERTAINTIES This chapter presents an abridged derivation of the Kalman filter for discrete and continuous linear systems, followed by a general statement of the extended Kalman filter for continuous nonlinear systems. The extended filter is then used as a base for the development of a specific algorithm for estimating states and parameters of the second order equation of pitching motion for a missile being forced by random disturbances. Development of the Kalman Filter The original derivation of the Kalman filter was presented by Kalmanl8 in 1960 for multistage systems making discrete linear transitions from one stage to another. Kalman and Bucy9 gave the analogous development for continuous linear systems approximately one year after the first work was published. The purpose of the linear filter is to provide estimates of the state of a system by making use of measurements of all, or some, of the state vector components of the system. The system is assumed to be operating in the presence of random distur bances, the statistical properties (i.e., mean and variance) of which are known. The measurements of the state vector components of the system are also assumed to have random errors of known statistics. In addition to the original derivations of Kalman and Bucy, a number of other methods offering various degrees of insight but leading to the same results are available. A brief development taken primarily from Bryson and Ho20 is presented here. Other derivations or developments of 31 32 the filter equations are given by Jazwinsky,3 Sorenson3 and Barham.33 The treatment given here starts with a static system, is extended for a singlestage linear transition which leads directly to linear multistage processes. Finally, by making use of a limiting process, the desired form of the equations for a continuous linear dynamic system is given. Static System The problem at hand is to estimate the ncomponent state vector X of a static system using the pcomponent measurement vector, z, containing random errors, v, which are independent of the state. The measurement vector can be represented as z = HX + v (21) where H is a known p x n matrix. Conditions on the measure ment errors are E(v) = 0 (22) E(vvT) = R (23) where R is a known matrix of dimension p. It is assumed that a prior estimate of the state, designated as X, is available and also that the covariance of the prior esti mate, M, is known. Thus T M = E[(XX)(XX) T] (24) where M is of dimension n. The desired estimate of X, taking into account the measurement, z, is the weightedleastsquares estimate, X, which minimizes J = [(X))T M1 (X() + (zHX)T R1 (zHX)]. (25) Differentiating Equation (25) with respect to X, setting the resulting equation to zero and solving for X yields S= X + PHTR1 (zHX) (26) where p1 = M1 + HTR1H (27) The quantity P is the covariance matrix of the error in the state estimate after measurement, X, i.e., P = E[(XX)(XX)T] (28) SingleStage Transitions It is now desired to estimate the state of a system that makes a discrete transition from state 0 to state 1 according to the linear relation X1 = 0oX + ro (29) where o is a known transition matrix of dimension n and r is a known n x r matrix. The mean and variance of the random forcing vector are assumed to be known and are given by E(w ) = wo E[(w wo)(wow ) Qo (30) The statistical properties of the random state vector are assumed to be known initially as E(Xo) = X0 E[(XoX )(X X )] = Po (31) It is also assumed that Xo and w0 are independent. From this information X1 is a random vector whose mean value and covariance are x = OoXo + ro.w (32) M1 = T + rQ T (33) Suppose now that measurements of the state are made after the transition to state 1; then from Equations (26) and (27) the best estimate of X is X, = X + PiHIR1 (zHTX) (34) where P, = (M 1 + HT1R1H )1 = M M H (HM HT + R1) HM, (35) Linear Multistage Processes The developments of the two preceding sections can be extended for linear, stochasitc, multistage processes. Given the following difference equation system model Xi+ = iXi + r.w. (36) where E(X ) = X (37) E(wi) = wi (38) E[(X X )(X o o)T = M (39) 00 0 E[(wiwi)(wjw.j)] = Qi6ij E[(wiwi)(XoXo) ] = 0 and measurements of the state zi = HiXi + vi 1 1 1 1 where E(vi) = 0 E(v v) = Ri6ij E[(w.wi)vj] = 0 11 , E[(XoRo)V] = 0 00 1 the estimate of the state is Xi = X + Ki(ziHiX) where .i+ = ii + riwi K. = P.HTR. 1 i l 1 P. = (Mi + HTR. lHi)1 1 i 11 1 = Mi M.HT(H.M.HT + R.1)H.M. 1 1 1 1 1 1 1 i 1 Mi+ = iPii + iQir i+1 1 1r1.Q 1. (40) (41) (42) (43) (44) (45) (46) (47) (48) (49) (50) Equations ('46), (47) and (48) are the discrete Kalman filter with the state variable covariances given by Equations (49) and (50). As can be seen from the above equations, the Kalman filter is essentially the same as the system model [Equation (36)]. The differences are (1) the actual system noise, which is random from one stage to the next in Equation (36), is replaced by its mean or expected value, and (2) there is a correction term based on the difference between the actual measurement of the state and its predicted value. The difference term is multiplied by a gain, Ki, which is essentially the ratio of the uncertainty in the state to the uncertainty in the measurement. If the covariances of the measurement errors are large, the gain will be relatively small and the corresponding difference term will have little effect on the estimate of the state; if, however, the system noise is relatively large or the measurement errors are very small, the gain will be large and differences between the actual and predicted measurements of the state at a given stage have increased significance in the state variable estimates. The Continuous Kalman Filter By applying a limiting process to the discrete filter with the time between stages tending to zero, the linear system model becomes X = F(t)X + G(t)w(t) (51) and the continuous Kalman filter is given by SHT ^ X = FX + Gw + PH Rl(zHX) (52) P = FP + PFT + GQGT PHTRHP (53) The Extended Kalman Filter The extension of the Kalman filter for estimating the states of nonlinear systems in the presence of noise has been given by Bryson and Ho20 and Jazwinski,31 among others, as X = f(X,t) + G(t)w(t) + P[ R1[zt)h(,t)] (54) ) aa af (hah S P + P fX + GQG' PTJ R1JP (55) 1x 1X axj P, ax] for the nonlinear system X = f(X,t) + G(t)w(t) (56) with measurements z(t) = h(X,t) + v(t) (57) where E[w(t)] = w(t) (58) E[w(t)w(t)][w(t')] = Q(t)6(tt') (59) E[X(t )] = o (60) E[X(t) j[X(to)Xo I = P0 (61) E[X(t )X ][w(t)w(t)]T = 0 (62) E[v(t)] = 0 (63) E[v(t)vT(t')] = R(t)6(tt') (64) E[X(t )Xo lv(t)]T 0 (65) E[w(t)w(t)][v(t')] = 0 (66) The Parameter Estimation Algorithm In this section a specific parameter and state esti mation algorithm using the extended Kalman filter is developed for the basic equation of pitching motion pre viously analyzed with the ChapmanKirk technique. The nonlinear equation of motion, or system model, is thus a + (C4+C6a2)& + (C3+C5a2)a = w(t) (67) where E[w(t)] = 0 E[w(t)w(t')] = q6(tt') (68) The measurements of the state of the system are assumed to be given by z(t) = a(t) + v(t) (69) where 1 ;, Ct) =o0 E[v(t)vT(t')] = r6(tt') S (70) Now to reduce Equation (67) to the required system of first order differential equations, let a = X1 & = X2 At this point, to estimate the aerodynamic parameters appearing in the equation of motion in addition to the state variables, the state vector is augmented by setting C3 = X3 C4 = X4 C5 = X5 Ce = X6 (71) (72) with the constraints = 0 = 0 = 0 = 0 (73) Making use of Equations (71), (72) and (73), Equation (67) now reduces to the following nonlinear system of first order equations: Xi = X2 12 = (X4+X6X12)X2(X3+X5X12)X1 + w(t) X3 = 0 X, = 0 is = 0 Xi = 0 with linear measurements z = Xi + v Comparing Equations (74) and (75) with (56) and (57) the following matrices may be identified: C(t) = H(t) = X1 0 0 0 0 0 f (X, t) = X2 (X4+X6Xi2)X2 (X3+X5Xi2)XI 0 0 0 0 4 (76) Now, since the mean value of the system noise is assumed to be zero JEquation (68)] and because the measurements (74) (75) of the state are related linearly to the state [Equation (75)], the extended Kalman filter, previously given by Equations (54) and (55), can be simplified somewhat to X= (X,t) + PH [zHX] (77) P f P + PI fHT q T X= p + +1 GqGT PHT 1HP (78) P x + where H = [ 1 0 Applying Equations (76) Xl X2 X2 (X4+XX12) X2  X3 0 X4 0 Xc 0 Xs 0 So 0 0 ] and (77) yields (X3+XsX1 2)X P13 P14 P15 P16 P23 P24 P25 P26 P33 P34 P35 P36 P34 P44 P45 P46 P35 P45 P55 P56 P36 P46 P56 P66 f(z[1 0 0 0 0 0] r PI1 P12 P12 P22 P13 P23 P14 P24 PIs P25 P16 P26 Carrying out the prescribed multiplications and equating like components of the resulting two column matrices yields the desired filter equations, = P + r (z) 2 r1 = (4+X6, X2 (X3+Xi25) + r (zX) P13  (zX1) r (zX) Sp15 r (zX1) r 16 (79) The necessary covariance equations are obtained from the matrix Ricatti equation [Equation (78)] making use of Equations (76) and the fact that xx= 0 2XiX2X6X33XiX5 0 0 0 0 1 (X4+X6eXi) 0 0 0 0 0 X2 0 0 0 0 0 _3 xl 0 0 0 0 0 XiX 0 0 0 0 J (80) The results of substituting Equations (76) and (80) into (78), carrying out the rather tedious matrix multipli cations and equating like terms are P11 = 2PI2 P11 r 1 12 = P22AP1iBP12XIP3iXz2PlCPIsDPiS P zlP12 P13 = P23P1P13 1 P14 = P24P11P14 P15 = P25 llP15 1 P16 = P26_P11P16 22 = q2(AP12+BP22++X1P23+X2P24+CP25+DP26)P12 P23 = AP13BP23X1P33X2P34CP35DP36ppl2p13 P24 = AP14BP24X1P34X2P44CP4SDP46P1l2Pl4 P25 = AP15BP25XIP35X2P45CP55DP56 P12P15 P26 = API6BP26X1P36X2P46CP56DP661pl2P16 P33 = P13 P34 = P13P14 1 P35 = 1I1315 P36 = P13P16 1 2 P44 = Pl4 P45 = IP14P15 I P46 = PI4P16 Ps 1 = Pl5 P56 = PP15P16 P66 1P162 (81) where ^2^ A = 2XIX2X6 + X3 + 3XIX5 B = X4 + XiX6 ^3X c = Xi D = X1X2 (82) The desired estimates of the states and parameters are obtained by numerically integrating the filter equations [Equations (79)] and their covariances [Equations (81)]. Initial estimates of the states and their covariances as well as the variances of the system noise and measurement errors are assumed to be available. In the event that the data consist of discrete measurements, which is usually the case, the constant time between data points should be equal to the numerical integration step size. This is necessary because the state estimates are updated at each integration increment and in doing this the filter requires knowledge of the difference between measured and estimated state values. If the filter successfully adjusts the true states and the imbedded parameters so that XI is an adequate match to the data, z, the time derivatives of the filter equations for the parameters become small [see Equations (79)] and the parameters reach steadystate or near steadystate conditions. CHAPTER V USE OF THE EXTENDED KALMAN FILTER FOR ESTIMATING PARAMETERS AND THEIR UNCERTAINTIES This chapter presents an analysis using the extended Kalman filter with parameter augmented state vector to determine aerodynamic coefficients and their uncertainties from noisy dynamic data. A description of the digital computer program used in the analysis is given first, followed by results obtained with this program from a variety of noisy data sets. The Extended Kalman Filter Program The extended Kalman filter program is written in Fortran IV for use primarily on an IBM model 360/65 digital computer. The basic function of the program is to integrate numerically the 27 first order differential equations which are given in the previous chapter, and which comprise the extended Kalman filter with parameter augmented state vector. The specific functions of the main program and its subroutines, as well as the required input data, are described in the following paragraphs. A listing of the complete program is given in Appendix II. Main Program Functions The main program reads and writes the input data and calls the numerical integration subroutine, ADDUM. Subroutine Functions The names and functions of each of the three sub routines that are used in the extended Kalman filter program are given below. ADDUM Subroutine ADDUM integrates the 27 filter and covariance equations using a fourth order RungeKutta method for starting, followed by a fourth order Adams Bashforth predictorcorrector method for running. It is essentially identical to the subroutine of the same name used in the ChapmanKirk coefficient extraction program and is described in detail in Reference 26. Subroutines XDOT and OUT are called from ADDUM. XDOT Subroutine XDOT computes values of the time deri vatives for the differential equations being integrated by ADDUM. OUT This subroutine writes the output of the numerical integration. Required Input Data The program reads four categories of input data. These categories and the specific data in each are delineated in Appendix II. The format and units of the entries on a specific data card can be determined from the program listing and nomenclature list, which are also given in Appendix II. Analysis of the Extended Kalman Filter The use of the extended Kalman filter for estimating the parameters of interest is analyzed for linear and non linear systems with data containing only measurement errors as well as data containing both system noise and measurement errors. The analysis begins with a linear system and data containing only measurement errors and progresses through increasing stages of difficulty up to nonlinear systems and data which have been corrupted by both system noise and measurement errors. The previously discussed computer program and algorithm require no modifications to consider the linear case; by initially setting the parameters that are the numerical coefficients of the nonlinear terms in the equation of motion (and their variances) equal to zero, the extended filter for a nonlinear system reduces to the one required for a linear system. All data used in the analysis were generated by the previously mentioned computer program, UFNOISE (Reference 29). The true values of the constants used in the data generation are those previously given in Table III with the exception of the initial value of the pitch angle for the linear system. This initial displacement is a more realistic 0.1745 radian for the linear cases considered. Linear Systems with Measurement Errors Only After the extended Kalman filter computer program had been constructed and checked, it became apparent that the technique would probably be best understood by considering relatively simple cases at first and then progressing to more difficult ones as confidence in the technique was gained. To this end, the first success with the program was realized when analyzing data containing only measurement errors for a linear system. The pitch angle data used in this part of the analysis were generated by integrating the equation of motion a + Cq4 + C3a = 0 (83) and superimposing random Gaussian errors on the output. The standard deviation of the errors is 0.00588 radian, which corresponds to 3c measurement errors of approximately 1.0 degree. The numerical integration step size used in generating the data was 0.005 second and the pitch angle was output every integration increment. The equation was integrated for a total of 2.00 seconds. Basic filter performance The results of using the extended Kalman filter to identify the correct values of the parameters of interest, X3 and X4, are shown on Figures 20 and 21. These figures show time histories of the percent errors in the estimated values of the parameters that are computed by the extended filter program when analyzing the noisy data described previously. The percent error is defined by (X.X.).O00 PERCENT ERROR (X.) = jJ (84) X. The input for the measurement error variance was the one previously computed by the simulation program based on the errors that were actually put into the data. The initial values of the parameter variances were computed from knowledge of the true value of the parameters and the arbitrarily selected initial values of the parameters. These sample variances are defined by P.jj(0) = [X.(O)Xj(O)]2 (85) This method of initializing the parameter variances was chosen merely to generate the required initial numerical values. In actual experimentation, the initial variance values would depend on the method of selecting the initial parameter estimates and prior knowledge as to how accurately these initial parameter estimates were in relation to the true values of the parameters. The initial values of all covariances were chosen as zero. This implies that errors in the estimates of the individual state vector components are initially uncor related. The initial values of the variances for the pitch angle and pitch angle rate were also chosen as zero. For wind tunnel dynamic experimentation where the model is initially held rigid at some given displacement to the wind stream, this seems to be a valid assumption. As can be seen from Figures 20 and 21, the filter does an excellent job of identifying the parameters of interest. The X3 parameter, which corresponds to the Cm term, and which is initially in error by 25 percent, is identified with approximately zero percent error in less than 0.3 second. The correct identification of the damping parameter takes slightly more time but the results are of equal quality. Effects of variations in the initial parameter variances This section presents some effects on the near steady state estimates of the parameters of interest and their near steadystate variances for various values of the initial parameter variances. The parameter variances were initialized to values that were 25 percent higher than those used in the previous section and also to values that were 25 percent lower than the referenced values. The effects of these initializations are shown on Figures 22 through 25 and in Table VII for the linear system of interest. Table VII is a summary of the percent errors in the parameter variances and the normalized parameter uncertain ties for the above mentioned variations in the initial parameter variances. The same results are depicted graphically on Figures 22 and 23. The numerical results in Table VII are based on the near steadystate parameter and variance values that result after the filter has integrated for 1.5 seconds. These correspond to the last point shown on the supporting figures. The normalized parameter uncertainties are defined by P?. PERCENT (P/P) =  100 (86) J3 n x.j Time histories of the parameter variances, P33 and P44, are shown on Figures 24 and 25. All of the information presented on these figures indicates that initial values Table VII. EFFECTS OF VARIATIONS IN THE INITIAL PARAMETER VARIANCES ON NEAR STEADY STATE PARAMETERS AND THEIR VARIANCES (LINEAR SYSTEM, MEASUREMENT ERRORS ONLY) Variation in Initial Parameter Variance* (Percent) 25 0 25 Errors in Near SteadyState Parameters (Percent) X3 X4 0.17 0.19 0.22 0.55 0.46 0.38 Normalized Near SteadyState Variances (Percent) 0.12 0.11 0.11 0.91 0.91 0.91 * relative to sample variance I of parameter variances in the range investigated have no effect on final parameter estimates or their variances when using the extended Kalman filter to analyze data for a linear system containing only measurement errors. Linear Systems with Measurement Errors and System Noise Data for this portion of the analysis were generated with the computer program UFNOISE by numerically integrating the following equation of motion, a + C4& + C3a = w(t) (87) and superimposing random Gaussian errors on the output. The forcing function, w(t), is also random in nature and Gaussian with zero mean. The standard deviation of the system noise was 4.80 rad/sec2, which is of sufficient magnitude to drive the oscillations in a noise band of approximately 0.0872 radian (5 degrees) for zero initial displacement. The standard deviation of the measurement errors was 0.00556 radian. Also, as before, the equation of motion was integrated in increments of 0.005 second for a total of 2.00 seconds. Basic filter performance The basic performance of the filter is demonstrated on Figures 26 and 27 which show time histories of the percent error in X3 and X4. As can be seen from these figures, the initial estimates of X3 and X4 are 25*percent in error and these estimates are corrected to within approximately 1 per cent and 3 percent of their respective true values. The noise variances used are those given by the simulation which generated the pitch angle data. The initial para meter variances are calculated from Equation 85. Effects of variations in the initial parameter variances The sensitivity of the filter to variations in the initial parameter variances is given in Table VIII and on Figures 28 through 32. Table VIII is a summary of the percent error in the parameters of interest and their normalized uncertainties for variations in the initial parameter variances of 25 percent, 0.0 percent and 25 percent relative to the sample variances. Figures 28 and 29 show the variation of the near steadystate errors in the parameters and their normalized uncertainties as functions of the percent variations in the initial parameter variances. Figure 30 shows the two time histories of the error in the estimate of X4 that result for two different initial values of the parameter variances. One of the initial variances is the sample variance based on the initial and true values of the parameters [equation (85)] and the other is 25 percent less than this value. The X3 error time histories are essentially identical for both cases and are not presented. The variations of Table VIII. EFFECTS OF VARIATIONS IN THE INITIAL PARAMETER VARIANCES ON NEAR STEADY STATE PARAMETERS AND THEIR VARIANCES (LINEAR SYSTEM, MEASUREMENT ERRORS AND SYSTEM NOISE) Variation in Initial Parameter Variance* (Percent) 25 0 25 Errors in Near SteadyState 'Parameters (Percent) 0.72 0.88 0.98 0.47 2.66 4.95 Normalized Near SteadyState Variances (Percent) 2.42 2.42 2.42 14.39 15.18 15.81 * relative to sample variance the parameter variances with time, for the three initial values considered, are shown on Figures 31 and 32. From the information given on the figures mentioned above, it is obvious that neither X3 nor P33 are affected by variations in the initial variances within the 25 percent range considered. The estimates of the damping parameter as well as its variance are affected, however, by these variations. The near steadystate damping para meter estimates vary almost linearly with initial variance values from a low of 0.47 percent to a high of 4.95 percent (see Figure 28). The variation in the normalized uncertainty is less pronounced, ranging from a low of 14.39 percent to a high of 15.81 percent (see Figure 29). Nevertheless, the near steadystate uncertainties are of sufficient magnitude so that the true values of the coefficients are within less than one standard deviation of their estimated values regardless of the initial parameter variance, when the standard deviation is taken as the square root of the near steadystate parameter variance. Attention is also called to the fact that since system noise is now present in the data, the near steadystate parameter variances are of greater magnitude than was the case when only measurement errors were present (see Figures 24, 25, 31 and 32). Filter response to large errors in the initial parameter estimates The results of using the extended Kalman filter to identify the damping parameter, X4, when the initial esti mate of the parameter is in error by 100 percent are shown on Figure 33 for two initial parameter variance estimates. When the initial parameter variance is based on the initial estimate of X4, the near steadystate parameter estimate is approximately 14 percent high; for an initial parameter variance that is greater by approximately 25 percent, the resulting near steadystate value of X4 is nearly 19 percent high. The corresponding parameter variances 1.5 seconds after the filter begins integrating are such that the estimated values of X4 are within approximately one standard deviation of the true value. Comparison of filter performance to ChapmanKirk technique The parameters of interest and their uncertainties were extracted from the given noisy data set using the previously described technique of Chapman and Kirk in order to be compared to the filter results. The percent error in the parameters and their normalized estimated standard deviations obtained with the ChapmanKirk program are X3:12.5 percent (1.8 percent) X4: 3.1 percent (0.3 percent) The above results are similar to those presented in Chapter III (Table VI) in that the estimated parameter uncertain ties are not of sufficient magnitude to reflect the true error in the extracted parameters. Whereas the extended filter estimates the parameters and variances so that the true value of the parameter is generally within one standard deviation of the estimate, it is shown that such is not the case when using the ChapmanKirk technique. The errors in the above parameters are also slightly greater than those obtained with the extended filter, although this is probably of less consequence than the differences in the uncertainties computed by the two techniques. Nonlinear System with Measurement Errors Only The pattern used in the previous two sections of generating a data set, investigating the basic response of the filter in identifying the parameters of interest and checking the sensitivity of the filter to variations in the initial parameter variances is essentially repeated here for the more complicated system model a + (C4+Ca2)E + (C3+C5a2)a = 0 (88) In addition, some effects of errors in the measurement error variance are also included and discussed. Some comparisons of results achieved with the extended filter and the Chl: aLnnKirk technique are also made. The standard deviation of the measurement errors in the generated data is 0.00569 radian. The mean value of the errors is zero. Basic filter performance The time histories of the four parameters of interest as computed by the extended filter while analyzing the data just described are shown on Figures 34, 35, 36 and 37. The initial errors in the parameters are 25 percent for X3, X4 and X, and 5 percent for X5. The measurement error variance is that obtained from the simulation. The initial values of the parameter variances are based on the difference between the true values of the parameters and the above mentioned initial estimates. The reason the initial estimate of Xs is only 5 percent low instead of 25 percent, as is the case with the other parameters, is due to its magnitude and a weakness in the numerical scheme. The true value of X5 is 1960 based on the data in Table III of Chapter III. A 25 percent initial error in this parameter results in an initial parameter sample variance of Pss(0) = (1960(0.25)1960)2 = 240,000 which is of sufficient magnitude to cause the filter to diverge for the numerical integration step size of 0.005 second. The divergence occurs because the large initial value of P55 results in large values of P25 and P15 which in turn cause the values of P55 that are fed into the numerical integration subroutine to be very large [Equa tions (81) and (82)]. The result of this cascading effect is a rapid decrease in P55 in large increments until it becomes negative (a physical impossibility since this term is the variance of a parameter). The structure of the covariance equations is such that a negative variance causes numerical divergence in the solutions to several of the equations. To circumvent the problem, an initial error of 10 percent in X5, with a corresponding adjustment in P55(0), was tried. The result was again divergence. Finally an initial estimate of 5 percent was found to be successful. The results shown on Figures 34, 35 and 36 indicate that the filter does an excellent job in correctly identi fying the linear and nonlinear static restoring moment parameters, X3 and X5, as well as the linear pitch damping term, X4. The results for the nonlinear damping parameter, Xg, are not as gratifying as can be seen by referring to Figure 37. The error in this term is still approximately 10 percent after 1.5 seconds. The reason for this problem is related to the magnitude of this parameter as compared to others in the equation of motion. The true value of X3 and XS are 160 and 1960 respectively, although the magnitude of X5 is effectively reduced by at least an order of magnitude when it is multiplied by a3. X4 and X6, on the other hand, have true values of 1.60 and 4.35 respec tively. After X6 is multiplied by a2& it is effectively the smallest parameter in the equation of motion by an order of magnitude and it is naturally more difficult to identify, since it has the least effect on the trajectory. Effects of variations in the initial parameter variances The effects of variations in the initial parameter variances are summarized in Table IX and on Figures 38 and 39. As can be seen from these results, the near steadystate values of the parameters and their variances are relatively insensitive over the range of initial para meter variances considered. Only the X6 parameter estimate and its variance show a variation of more than 1 percent. As is obvious from the referenced table and figures, no data are available for initial parameter variances of 50 percent above the sample variances; the reason for this is, again, that the magnitude of the P55 variance caused the filter to diverge. Attention is called to the fact that for the range of initial variances considered, the near steadystate estimates of the X3, X4 and Xs parameters are all within one standard deviation of their true values when the standard deviation is taken as the square root of the near S9 Zf to to o ; (t) H Cd a) O 4 Uo man (d *n * <  0 0 0 0 UL4 p1 r 4 Lnr 00 nI r* H 0 "H $ , Cr M r1 4 " Z SH rHt 0 . H 0 0< t H 4 ,4 o A ;2 P0 o oD 0 m n H Hc H * U 4 00 0 0 0 0 C a ( 00 E 4 04U u< wU el CI) C) C PZ z 1 *.dC C%: 1dr O r0 )O C Q ( ; 0 cd i4 ( C r* r4 Ft;_ d J () c ) r) rl r tl *K steadystate variance for the parameter of interest. The Xg parameter is always within two standard deviations of its true value. Effects of errors in the estimate of the measurement error variance Up to this point, exact values of the measurement error variance, r, have been used in all runs made with the extended filter. Exact knowledge of this quantity is avail able because the pitchangle data being used with the extended filter are computer generated with specified noise statistics. In actual test situations, true values of the noise para meters may not be known exactly and, in fact, methods for determining both measurement error variances and system noise variances from noisy dynamic data are of current research interest.3 Some effects of errors in the estimates of the measure ment error variance on the near steadystate parameter estimates and their uncertainties are given in Table X and on Figures 40 and 41. These results show that errors in the measurement noise variance of 25 percent have essentially no effect on estimates of the four parameters of interest. The uncertainties of the X3, X4 and X5 parameters are changed by less than 1 percent over the range considered. The normalized Xg uncertainty varies from approximately 5 percent to 12 percent. SV) O CO0  C/) 2 0 03 CO) O O CO MH o E 0a DiM K; U ca te r,.  4 c1 rI 4. IJ C .t I  0a or.c) F: DD 0I 04 t4 t C! I " in ro ^ z c! >< s) m >. (a) *a 4 ;V4 m c m *r4 P $4 i O  o~g ce W: r Md wf 0P CLn M cn Filter response to larger error in the initial estimate of Xs The difficulty associated with.numerically large initial parameter variance estimates has been discussed earlier. The sample variance of P55 which corresponds to a 5 percent error in the initial X5 parameter estimate appears to be an approximate upper limit for the integration step size and physical constants used in this analysis. A run was made, however, with the previously described nonlinear data set containing only measurement errors where all initial parameter estimates were in error by 25 percent; the P55 variance, however, was set equal to the previously used sample variance which corresponds to a 5 percent initial estimate of X5. The parameter estimate time histories resulting from this initialization are shown on Figures 42 through 45. These results indicate that the near steadystate values of the parameters are very nearly equal to those presented when the sample variances based on the initial parameter estimates were used initially. Comparison of filter performance to ChapmanKirk technique The linear and nonlinear static restoring moment and damping parameters were extracted from the noisy data set being considered in this section using the ChapmanKirk program. Th:e percent errors in the parameters and their normalized estimated uncertainties are given below: X3:0.36 percent (0.35 percent) X4:0.69 percent (0.43 percent) X5:0.41 percent (0.49 percent) X6:4.60 percent (11.27 percent) Comparing the above results with those previously given in Table IX reveals that the estimates obtained using the extended filter are slightly more accurate for all but the X4 parameter. Since the errors are, for the most part, less than 1 percent the differences are essentially insignificant. The normalized uncertainties obtained with both techniques are equally good. Nonlinear System with Measurement Errors and System Noise In this section, as before, the basic filter performance, sensitivity of the filter to initial parameter variance vari ations and sensitivity of the filter to errors in the esti mates of the noise variances are investigated. The data used in the analysis were generated by numerically inte grating the following nonlinear equation with a random forcing function, E + (C4+C6a2)& + (C3+C5a2)a = w(t) (89) Approximately the same noise variances used previously are repeated once again; the variance of the measurement errors is 0.00576 radian and the variance of the system noise is 4.75 rad/sec2. Basic filter performance The time histories of the errors in the parameters of interest obtained with the extended filter are given on Figures 46 through 49. The response time required for the parameters to reach their near steadystate values is slightly greater than in the simpler cases that have been considered previously. The accuracy of the parameters is similar to that which was achieved for the nonlinear system with only measurement errors in the data with the exception of X4, which is approximately 5 percent higher. Effects of variations in the initial parameter variances The effects of variations in the initial parameter variances are given in Table XI and are shown on Figures 50 and 51. Only the nonlinear damping parameter and its uncertainty vary by more than 1 percent for initial variance variations of up to 50 percent relative to the sample variances. Effects of errors in the estimates of the noise variances The effects of errors in the estimates of r, the measure ment error variance, and q, the system noise variance, are given in Table XII and are shown on Figures 52 and 53. These figures show the variation of the near steadystate parameter errors and normalized parameter uncertainties as functions of the error in the estimate of the measurement error variance for various values of error in the system U 1w c* *r > 1d (1 41 a *rH C Q) r4 4J U Z ci) 0 ( cI m (^1 ' t I 14 CO 2 P 0 U Z 2Z < HZ < z I  U) CO MH PCu)> dHZ Sr4 ~C ci)D cn H ? 4: ,1 O(UN *8 cd O 4 04 c u r: + cd rd U Cd c l * *H q N 4 ) N Cd Cd a4 cd *H > *H '4< 9 43 cdCd 0 Cd c) 24 4J 12 cd e 0) *rl X +> (1, U) > ( (M W N N r4 11 fN 00 4* * . *Co 4 co (* * 0ct 00 Sl 0 0 cd 00 0 rt C)l CO 4 In cd 0 0 N *r4 ci 4C ov, 0 11 zWZ P4 m P4 :2: S< P4 = :1 En ZO4 E< 0 E OSH 0 0< W 0 0Z2 O (. m 0 c oo p 0 i P< d 0 U) 0 cd F4 r4 4 U rt 0r 3 ! U) cU a 9 0 Z c $ 0 CS 0 0 0 to COd W Cd o ri C U)  4 4l4 LO (M 0 t) o LC ) L)n o in i noise variance estimate. From these results it is obvious that only X6 is affected by errors of 25 percent in the estimate of the system noise variance or measurement error variance. Comparison of filter performance to ChapmanKirk technique The percent errors in the parameters and their normalized estimated standard deviations that are obtained using the ChapmanKirk program with the data currently being analyzed are given below: X3:0.62 percent (0.40 percent) X4:2.63 percent (1.23 percent) Xs:0.41 percent (0.61 percent) X6:46.3 percent (14.8 percent) Comparing these results with the results in Table XI indicates that both techniques yield essentially the same accuracy in their determination of X3 and X5. The filter does considerably better in estimating X6 but is slightly worse in its estimate of X4. The uncertainties computed using the filter are of sufficient I.::gnitude so that the true value of the X3, X5 and X6 parameters are within one standard deviation of the filter estimates; the true value of X4 is within two standard deviations. 