Home

User's Manual for MMLE3, A General FORTRAN Program

1. EM 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 G ds ds gs sei 2 a g mg mg mg gi JR 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 The FM matrix is filled with 1 s AL 0 sina COS Q cos q cos 9 g 0 0 0 0 0 0 0 0 0 1 cos q tan 9 0 The BL SL RL CL DL HL EL and FL matrices are filled with 0 s The vectors v and w in the lateral directional equations depend on TIMVAR sec 3 3 8 20 and MX sec 3 3 8 11 All of the terms in v and w are computed using only measured data The average measured V is used in the v 1 equation if TIMVAR is FALSE otherwise V is redefined at each time point si g v 1 sin cos 9 v i ALi 4 9 YOFFA where YOFF is the vector of biases in the measurement equation If USERIC sec 3 3 8 26 is TRUE YOFF is a 0 vector otherwise YOFF is the vector of initial measurements If MX is less than 4 the y YOFF term is omitted This implies that the program uses the measured alone rather than linearizing about it 99 4 2 DEL Laz VASTET g PA R Xx Ty y Lez 6Nq ee VOD P TEMA R T z y 64 v 4 q sin tan 9 oe a p r2 R g 4 2 Files This section describes the file formats used by the standard aircraft routines Only details specific to the standard aircraft routines are mentioned here General descriptions of the files are found in section 3 2 4 2 1 Input Time History File The standard aircraft routine READTH reads time history data
2. eR A c X tis ER 1s p u t 1 40 ER 1s H 4 ER v t w t 1 where o R AAt a y RAs ds 42 0 and where throughout this section an index of i indicates the average of the value at the beginning of a time interval index i and the value at the end of a time interval index i 1 One special case must be noted The predicted state X at the beginning of a time interval is the same as the corrected state at the end of the previous interval Thus T tipi is the average of t and ti41 not x t and t The same principle applies to Lo Correction step Be t 1 Ea tins K 2 tis1 Ze ti 1 43 Gradients of matrices are involved in several of the equations below The gra dient of a matrix is a third order tensor Formal tensor notation is avoided because the operations involved should be unambiguous without the multiplicity of indices formally needed For instance the term Y5 represents the tensor inner product Y B pu Those uncomfortable with tensors can substitute partial deri vatives for the gradients Each gradient equation then becomes a set of partial derivative equations involving no tensors The first gradient of the filter equations is as follows 24 Prediction step Vee tins OVER ti VR OA He ting VEB Cina 85 72 ti e Era c y A ER V A where we have defined Correction step where eR AX C 1 Vets
3. second gradient transpose predicted estimate corrected estimate at flight or reference center of gravity corrected general indices measured rotary derivatives per rad static derivatives per deg control derivatives with respect to indicated quantity function of bias Prefix to matrix names APR a priori weighting 1 0 Suffixes to matrix names inverse L dimensionalization addition M dimensionalization ratio N nondimensional V variation 1 0 PARAMETER ESTIMATION THEORY The parameter estimation problem can be defined quite simply in general terms The system investigated is assumed to be modeled by a set of dynamic equations containing unknown parameters The system is excited by a suitable input The input and the actual system response are measured The values of the unknown parameters are then inferred from the requirement that the model response to the same input should match the actual system response Formulated in this manner the identification of the unknown parameters could be easily accomplished by many methods however complicating factors arise when considering application to real systems The first complication results from the impossibility of obtaining perfect measure ments of the response of any real system The inevitable sensor errors are usually included as additive measurement noise in the dynamic model Once this noise is introduced the theoretical nature of the problem ch
4. 89 1 Report No 2 Government Accession No 3 Recipient s Catalog No NASA TP 1563 4 Title and Subtitle 5 Report Date November 1980 6 Performing Organization Code 8 Performing Organization Report No H 1084 10 Work Unit No 9095 36 24 11 Contract or Grant No 13 Type of Report and Period Covered Technical Paper 14 Sponsoring Agency Code USER S MANUAL FOR MMLE3 A GENERAL FORTRAN PROGRAM FOR MAXIMUM LIKELIHOOD PARAMETER ESTIMATION 7 Author s Richard E Maine and Kenneth W Iliff 9 Performing Organization Name and Address NASA Dryden Flight Research Center P O Box 273 Edwards California 93523 12 Sponsoring Agency Name and Address National Aeronautics and Space Administration Washington D C 20546 15 Supplementary Notes 16 Abstract This report is a user s manual for the FORTRAN IV computer program MMLE3 a maximum likelihood parameter estimation program capable of handling general bilinear dynamic equations of arbitrary order with measurement noise and or state noise process noise This volume describes the theory and use of the pro gram The program listings and discussion of the inner structure and coding of the program are found in the Programmer s Manual for MMLE3 A General FORTRAN Program for Maximum Likelihood Parameter Estimation by Richard E Maine NASA TP 1690 The basic MMLE3 program is quite general and therefore applicable to a wide variety of problems
5. P cos a tr sin a AL 3 8 YOFF my Sin a where YOFF is the vector of biases in the measurement equation If USERIC sec 3 3 8 26 is TRUE YOFF is a 0 vector otherwise YOFF is the vector of initial measurements If MX is less than 3 the 0 YOFF term is not included This implies that the program uses measured 6 alone rather than linearizing about it 91 4 1 3 2 li I _ Z_ 2 2 xz A ae LR rQ p a v 3 r sin q Pp w 1 YALF V VAN ZAN 2 2 w 4 gO Pt gt a p g R g _XAX DCGFT 2 2 _ YAX p T w 5 5 Hog R I p Xe I y 6N 63 Rg mg 4 1 3 2 Lateral directional The nondimensional lateral directional matrices are 52 Y Cy Cy 0 B p r o Co Co 0 B p r C C 0 n n n B p r 0 0 0 0 SN Gu Ep Ya n C Lo Ch 0 Po 4 1 3 2 RN CN 1 0 0 0 KB KBXZB KBX XB DCGFT 0 Do 0 1 0 0 ie ae Le 0 0 0 B x 0 0 0 1 XZ Ca TC 1 0 om Cy z 0 B Z B p 0 0 0 1 0 0 0 0 0 0 DN HN 0 0 0 0 0 0 0 0 0 0 0 0 C C C Es a Yo a r EN 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ZAY XAY DCGFT 0 0 1 0 0 0 0 1 0 93 4 1 3 2 The dimensionalization matrices defined by the standard aircraft routines for a lateral directional case are BM AM S gQ Q Q ote on ah 7 S S S Q GE as qa 7 AAA S al ala dg 2 8 al el als als RM SM DM CM ls ca ot H H 94 4 1 3 2 HM
6. 11 ZALF ZB ZAN ZAX ZAY distances of the angle of attack angle of side slip normal acceleration longitudinal acceleration and lateral acceleration sensors respectively below the center of gravity ftor m The default values are 0 4 3 2 3 Angle of attack break points The angle of attack break points are listed on one to three cards in 8F10 format 4 3 2 4 Mach number break points The Mach number break points are listed on one to three cards in 8F10 format 4 3 2 5 Param break points The break points for the extra distinguishing parameter param are listed on one to three cards in 8F10 format 4 3 2 6 Predicted derivative data This section describes the predicted data input format for one derivative The format is repeated for as many derivatives as are desired no particular order of the derivatives is required For each derivative there is a derivative header card followed by a derivative data table The header card contains the derivative name type location and functional dependence The derivative name is four characters starting in column 1 the name is not used by the program but is solely for ease of user identification The type should be either LONG or LATR Starting in column 11 to indicate a longitudinal or lateral directional derivative if anything else is specified for type the derivative will effectively be ignored The derivative location is specified by matrix name row and column
7. FULL1 19 INCH 38 ITAPR 22 ITDFAC 21 ITG 24 MAXREC 5 MB 14 MU 13 MX 11 MZ 12 NCASE 1 NEAT 15 NEXPLT 37 NOITER 16 NREC 6 NUPLT 36 PAPER 39 PLOTEM 33 PLTFAC 38 PLTMAX 34 PRINTI 29 PRINTO 30 PRINTX Gl PRINTY 31 PUNCH 45 RELAB 10 RLXG 24 SPS 4 TAPE 2 TEST 28 THIN 3 TIMESC 40 84 Page 39 NAMELIST INPUT Continued TIMVAR 20 UBIAS 9 UCHAN 7 UMAX 42 UMIN 42 USCALE 8 USERIC 26 VARIC 27 WAPR 22 XMAX 43 XMIN 43 XPLOT 35 ZBIAS 9 ZCHAN 7 ZMAX 41 ZMIN 41 ZSCALE 8 NAMELIST USER ection 4 3 3 ALPHA 17 AREA 4 CG 10 CHORD 6 IX 8 IXE 9 IXZ 8 IY 8 IZ 8 KALF 11 KB 11 LATR 1 LONG 1 MACH 20 METRIC 3 PARAM 21 PHI 19 QBAR 15 SHIFT 2 SPAN 5 THETA 18 UVAR 22 V 16 W XALF 12 Page NAMELIST USER Continued XAN 12 XAX 12 XAY 12 XB 12 x YALF 13 YAN 13 YAX 13 YAY 13 YB 13 ZALF 14 ZAN 14 ZAX 14 ZAY 14 ZB 14 NAMELIST WIND Section 4 3 2 2 AREA 7 CG CHORD 7 KALF 8 KB 8 Page NAMELIST WIND Continued NABP 2 NMBP 3 NPBP 4 PRINT 6 SPAN 7 STAB 1 XALF 9 XAN 9 XAX 9 XAY 9 XB 9 YALF 10 YAN 10 YAX 10 YAY 10 YB 10 ZALF 11 ZAN 11 ZAX 11 ZAY 11 ZB 11
8. MAXZ MAXU and LEX word vectors respec tively biases for the observations controls and extra signals respectively Elements of ZBIAS UBIAS or EXBIAS are added to corresponding measured data to compensate for known bias errors These additions are performed after any scale or sign corrections item 8 but before any other operations on the data The default values are all 0 10 RELAB logical option to read labels If RELAB is set to TRUE labels for the observations states controls and extra signals will be read as described in section 3 3 9 If RELAB is FALSE the default labels used are Zi Xi Ui or EXi for the ith observation state control and extra signal respectively These default labels are changed by the standard aircraft routines The default value of RELAB is FALSE 11 MX length of the state vector All matrix dimensions will be forced to conform with MX states The maximum value for MX is given by the variable MAXX sec 2 currently 7 The minimum value for MX is 1 If MX is specified as negative it will be obtained from the number of rows of the AN matrix sec 3 3 11 The standard aircraft routines define a default number of rows of AN sec 4 3 5 The default value for MX is 1 12 MZ length of the observation vector All matrix dimensions will be forced to conform with MZ observations The maximum value for MZ is given by the variable MAXZ sec 2 currently 8 The minimum value for MZ is 1
9. RLXG G determination variables If ITG is not 0 G determination will start after ITG normal iterations or convergence item 17 whichever occurs first Normal iterations are defined as iterations without any iteration dependent options Iteration dependent options include DFAC item 21 and the linear terms only option item 19 IfITAPR item 22 is nonzero then a priori is considered an iteration dependent option If BOUND item 17 indicates convergence after G determination has started iteration terminates If ITG is 0 G determination is not used RLXG is a relaxation factor used in the G determination The default values are 0 for ITG and 1 0 for RLXG 25 DIAGG logical diagonal G matrix option If G determination item 24 is not used DIAGG is ignored If G determination is used DIAGG specifies whether a diagonal or full matrix will be determined A value of TRUE for DIAGG specifies determination of a diagonal G matrix If a nondiagonal G matrix is used for a starting estimate DIAGG is forced to FALSE the printout of the options may be incorrect in this case because the printout precedes the matrix input The default for DIAGG is TRUE 26 USERIC logical user initial condition option If USERIC is TRUE user routine UINIT will be called to determine the initial conditions and fixed biases for each maneuver item 1 Some implications of the use of UINIT are discussed in section 3 1 UINIT should define the i
10. The matrix name can be AN BN SN RN CN DN HN or EN starting in card column 21 A left parenthesis separates the matrix name and row number a comma separates the row and column numbers and a right parenthesis follows the column number The row and column numbers can be either one or two digits No blanks are allowed before the closing parenthesis except as the first character of a two digit row or column number Examples of permissible forms of the derivative location are AN 1 1 AN 01 1 and ANC 1 1 The functional dependence for each derivative is specified by one two or three of the characters A and M and P in any order starting in column 31 The character A stands for angle of attack dependence M for Mach number depend ence and P for param dependence Blanks in columns 31 to 33 are equivalent to AMP any character other than A or M or P or a blank in column 31 indicates the derivative is a constant independent of all three variables The form of the data table depends on the functional dependence indicated on the header card If a constant derivative was indicated the data table is a single card with the constant value in F10 format If the derivative depends on one parameter the derivative values at the break points for that parameter are listed on 61 4 3 2 6 one to three cards in 8F10 format If the derivative depends on two parameters the one parameter form is repeated for each break point o
11. a q 0 0 0 BN SN T n Cr a 0 e C mg 0 6 i 0 KALF X XALF DCGFT 0 1 0 0 1 ES 0 q Es 0 q 0 0 HN EN 0 0 0 0 0 0 0 0 0 C 0 XAN DCGFT No 0 ZAX C A 0 1 OS O oOo OO o DN S 49 4 1 3 1 The longitudinal dimensionalization matrices defined by the user routines are AM _ WS y qs e m2 myy qse qsc Cc TR Fw y y 1 1 1 SM dsg Us 38 my mv A my qsc _qSC qsc I R I R I R y y y 1 1 1 CM 1 1 y 1 1 1 1 1 1 1 ds gc mg mg 2VR gs _qs 1 mg mg 2VR 90 BM _ gs 98 y 98 g mv A mv XK my I R I R I y y 1 1 1 RM _ gs 1 mvt qsc qsc I LER y y 1 1 DM 1 1 1 1 1 1 1 1 1 qs qs qS mg mg mg 4S qS _ 48 mg mg mg ale af wae 4 1 3 1 EM HM 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 5 1 qs qs qs q5 g mg mg mg mg 1 qS qS qs qs IR mg mg mg mg The FM matrix is filled with 1 s AL q _ 0 1 7 cos sin cos a cos 0 sina C 0 0 0 0 coso 0 The BL SL CL DL HL and EL matrices are filled with 0 s The remaining quantities to be defined for the longitudinal equations are the known forcing functions v and w All of the terms in v and w are computed using only measured data The definition of V depends on TIMVAR sec 3 3 8 20 and MX sec 3 3 8 11 The average V is used if TIMVAR is FALSE otherwise V is redefined at each time point v 1 y A cos O cos q cos a sin 6 sin A B
12. c sin a 2 cos a 66 _ sin a _ cos a _ sin a qc _ cosa qc _ sina _ cos a RN Ay CN IVR R CA VR R NOR ee where the 5 terms are summed over all controls The average values of a q and are used If MZ sec 3 3 8 12 is 4 the above computation is replaced by the low a approximations Cr CN and Cr CN q q a a 68 4 3 5 The lateral directional AN default remains Cy Cy Cy 0 B p r Co Cy Cy 0 B p r Ca Ca om 0 B p r 0 0 0 0 The default size of AN is set to 3 by 3 for a longitudinal case or 4 by 4 for a lateral directional case The longitudinal BN default is Ls e On e 0 where MATDEF computes C C ecosa C sina If MZ sec 3 3 8 12 is Le Ng As e e e 4 the approximation C C is used re Ng e e The lateral directional BN default remains C C Ys Ys a r C C Y Le a r Ci Ci 8 8 a r 0 0 69 4 3 9 The longitudinal SN default is where MATDEF computes C C cosa C sina If MZ sec 3 3 8 12 is 4 h No S0 the approximation Cr Cy is used The default size of the longitudinal SN is set to 3 by 1 0 0 The lateral directional SN default remains The default size of the lateral directional SN is set to 4 by 1 The RN matrix defaults are Longitudinal Lateral directional 1 0 0 1 0 0 0 C 1 0 _IXZ Ms da 1 ix 0 0 0 1 _IXZ Cm a P 0 0 0 1 where IX IZ and IXZ are described in section 4 3 3 8 70 4 3 5 The CN matrix defau
13. continuation cards The standard deviations of the measured observations controls and extra signals follow The format is the same as for the averages except that SIG replaces AVG in the labels The signal minima follow in the same format labeled MIN then the maxima are punched labeled MAX Following the above header information the relevant nondimensional matrices and Cram r Rao bound matrices are punched in standard matrix format sec 3 3 11 Any nondimensional matrices that contain no independent unknowns are omitted along with the corresponding Cram r Rao bound matrices The last card punched contains ENDCASE in the first eight columns 38 4 2 4 4 2 4 Output Time History File The standard aircraft routine THOUT is intended for use in creating a simulated data file An unformatted FORTRAN file is written which can be used as an input time history file in a later run The file has one record per time point with time in integer hours minutes seconds and milliseconds 1n the first four words of the record The rest of the record contains the final iteration predicted observation time history the controls and the extra signals The complete dimensioned lengths of the observations controls and extra signals are written These lengths are MAXZ MAXU and LEX respectively sec 2 4 3 Input Description This section describes the card input required by the standard aircraft routines The input description o
14. ref 48 is used for the solution of equation 29 The Hamiltonian matrix is defined as L cr 66 y lo A The QR algorithm ref 49 is used to find the eigenvalues and normalized eigen vectors of the Hamiltonian matrix The matrix of normalized eigenvectors is par titioned into four equal size matrices X111 12 3 X21 X22 where the eigenvectors corresponding to eigenvalues with positive real parts are in the left partition Controllability by the state noise and observability is sufficient to guarantee that exactly one half of the eigenvalues will have positive real parts and that the solution to equation 29 is given by PSK oy 39 The gradient of P is the solution to a set of Lyapunov equations one for each partial derivative of the general form A V P t Ye M AF 31 The eigenvectors of 4 are computed and used to transform into 4 as follows A T AT 32 where 4 is block diagonal with 1 by 1 blocks for real eigenvalues and 2 by 2 blocks for complex eigenvalues T is a real transformation matrix the columns of which are the real eigenvectors and the real and imaginary parts of the complex eigen vector pairs Then defining 1 1 V P T V P T 33 ee ee a Ll A1 T C T 34 we have the block diagonal Lyapunov equation SV P V P L 35 l ee a 19 2 0 Equation 35 is separable into 1 by 1 2 by 2 and 1 by 2 par
15. 0 The o and are included to allow for instrument biases The q equation al ready has the freedom to correct such biases using Co the Cr in the equation 0 0 does not allow similar freedom because it is related to the Cn and Ca terms in the observation equations below Measured data are used for B p r and q in these equations Measured a is also used wherever a appears explicitly above The gravity term in the equation is normally linearized about the measured 9 however it can be evaluated at the measured 9 if it is not desired to integrate the 9 equation 44 4 1 2 1 The longitudinal observation equations are T Ya Tm Ka t q4 p Gm q O 9 m 96 x z X a a a _ gs ns h q la ns a Cy zd q p RP Nn mg N g2 22g IR Za lt q Ya _ _ Qs Xa KY 2 eS are Y a Ca tq gt r mr Xn mg A JRI mi IR mg Im 7 17 d The do is the instrument bias on q Biases on a and a are included in Ca and 0 Ca those in a q and 9 were handled in the state equations Measured data are 0 used for p r p and in these equations The terms involving q are also eval uated with measured data to eliminate the nonlinearity The expansions of the longitudinal force and moment coefficients are C qe Cn Cy 2 Cy YE Cy 8 Cy a q 6 0 7 qc ae Cin ae a VR Cm ha C IWR a q 0 57 qe Ca ie j lt A ZUR A b En Cy Cy cosa C sin a where the 6 term is sum
16. 2x 6 1 1 The last term is constant and can be neglected in the minimization If Gis known the next to last term is also constant Multiplying by 1 to get a minimization rather than a maximization problem results in the following simple and intuitively attractive cost function J 22 RO 2 t gr 1 2 7 1 1 2 Measurement Noise and State Noise When both measurement and state noise are present there is no known practical method of computing the likelihood ratio for nonlinear systems Attention is therefore restricted to linear and bilinear systems The assumed equations of the actual system are x t Ax t But Fn t z t Cx t Du t n ti where n is the state noise vector The matrices A B C D and F are functions of In general A B C and D can include known functions of time the time depend ence is omitted from the notation for simplicity The state noise vector n is assumed zero mean Gaussian and white with unity spectral density The measure ment noise vectors n t gt are assumed to be zero mean Gaussian independent 8 random vectors with identity covarience The measurement and state noise are assumed to be independent The log likelihood ratio is then z log p Z 3 Ze ti z t GG 1 eC z t 5 log IGG m log 2x 9 1 1 where Z is the predicted estimate of z and GG is the variance of the filter resid uals Note that this e
17. 9 60 Z X a Ya a Wo Ys Pah once pay Ym mg Y gR GR Rg Pm P Po A The p and are instrument biases on p and Biases on a are included in C 0 0 y Y those on B p r and q were included in the state equations Measured data are used for p and po in the a equation to eliminate the nonlinearity The f observation equation does not include the division by cos a which is needed if B is measured by a fixed vane The expansions of the force and moment coefficients are L b rb Cy Cy B Cy apa Cy gwa Cy B Cy B p r 5 0 Cy Cy B C goby Cy yuyg Cy t Co Cy in Nhat p p r 6 0 B pb ah _Bb_ En ag id n IVR j Cn WVA id li gt ng ng 2VR where the 6 term is summed over all controls 47 4 1 3 4 1 3 Matrix Equations This section describes the linearized equations of motion in the matrix form exactly as used by the standard aircraft routines The variables TIMVAR sec 3 3 8 20 and MZ sec 3 3 8 12 affect the exact forms of the linearized equations used The general form of the equations is R A P x Bmut SIDAD v Fn t 36 z t C x D u t HHIH E w Gn 1 where A B S R C D H and E are defined by relationships such as A AM t X AN os t 37 Further discussion of the general equations and gradients is found in section 3 1 The N matrices for the standard aircraft routines are further discussed in sections 4 3 2 6 and 4 3 5 The
18. BN SN RN CN DN HN EN nondimensional starting matrices There are two types of default values for these matrices First the data on the predicted derivative file sec 4 2 2 second the specific values set in routine MATDEF The values set in MATDEF will override the values from the predicted derivative file if they refer to the same matrix element Data from the predicted derivative file are used only if the predicted derivative control card sec 3 3 4 specifies OLD or NEW These data can contain tables for any location in AN BN SN RN CN DN HN or EN The table values are inter plotted using ALPHA MACH and PARAM sec 4 3 3 17 20 and 21 Only tables with a type of LONG will be used for longitudinal cases or LATR for lateral directional cases Details of the predicted derivative data are given in sections 4 2 2 and 4 3 2 The remainder of this section discusses the computations made by MATDEF and shows the resulting default matrices It is assumed that the predicted derivative data specify the derivatives and matrix locations shown in section 4 3 2 6 If SHIFT sec 4 3 3 2 is TRUE the moment derivatives from the predicted derivative file are corrected from the reference to the flight center of gravity location The default longitudinal AN is Cr Cy 0 o q C C 0 m m 0 q 0 0 0 where MATDEF computes AN 1 2 C C cosa C sin a N A q q q B 7 P a a AN 1 1 ga Cy cos a 2 sin a
19. CN DN and HN matrices instead of C N and C A derivatives Subroutine WTTRAN will convert the Cy and Ch derivatives to the required Cy and Ca derivatives Note that the Cy derivatives should not be placed directly in the AN BN and SN matrices the Cr derivatives in these matrices are computed by MATDEF regardless of the axis system of the predicted data For a lateral directional case the following matrix locations are those normally used RN BN Cy Cy C 0 0 0 0 0 Oy Cy p p r Cy 0 0 0 a r oF Co 0 c c B p r Ca 0 0 0 5 5 C C 0 B p r 0 0 0 0 ns ns 0 0 0 0 7 0 0 63 4 3 2 7 4 3 2 7 End card The end of the predicted derivative file is indicated by a card with the characters END starting in column 1 4 3 3 User Input User input sec 3 3 7 for the standard aircraft routines consists entirely of a NAMELIST called USER A brief discussion of NAMELIST format is contained in section 3 3 8 for complete details consult the FORTRAN reference manuals for specific computer systems All values are reset to the defaults at the beginning of each case regardless of any values used in previous cases If user routines are bypassed sec 3 3 2 user input is not included in the deck The following variables are included in NAMELIST USER 1 LONG LATR logical type of aerodynamic mode to be analyzed Either LONG or LATR but not both should be set to TRUE to indicate a longitudinal or lateral directio
20. D 7831 1975 Suit William T Aerodynamic Parameters of the Navion Airplane Extracted From Flight Data NASA TN D 6643 1972 27 28 29 30 31 32 33 34 35 36 37 38 39 40 82 Klein V On the Adequate Model for Aircraft Parameter Estimation Cranfield Report Aero No 28 Cranfield Inst Technol Mar 1975 Steers Sandra Thornberry and Iliff Kenneth W Effects of Time Shifted Data on Flight Determined Stability and Control Derivatives NASA TN D 7830 1975 Iliff Kenneth W and Maine Richard E Further Observations on Maximum Likelihood Estimates of Stability and Control Characteristics Obtained From Flight Data AIAA Paper 77 1133 Aug 1977 Park Gary D Parameter Identification Technology Used in Determining In Flight Airload Parameters J Aircraft vol 14 no 3 Mar 1977 pp 251 256 Iliff Kenneth W Maine Richard E and Montgomery T D Important Factors in the Maximum Likelihood Analysis of Flight Test Maneuvers NASA TP 1459 1979 Maine Richard E Aerodynamic Derivatives for an Oblique Wing Aircraft Estimated From Flight Data by Using a Maximum Likelihood Technique NASA TP 1336 1978 Iliff Kenneth W Maximum Likelihood Estimates of Lift and Drag Characteristics Obtained From Dynamic Aircraft Maneuvers Proceedings of AIAA 3rd Atmos pheric Flight Mechanics Conference c 1976 pp 137 150 Kirsten Paul W and Ash Lawrence G A Comparison a
21. FALSE plots are sized for centimeter grid paper Other values of PLTFAC multiply the plot size proportionally The default value is 1 for PLTFAC and FALSE for INCH 39 PAPER paper grid width centimeters or half inches Plots will be scaled to fit on paper with a grid width specified by PAPER The units for PAPER depend on INCH item 38 The default value for PAPER is 25 centimeters if INCH is FALSE or 20 half inches if INCH is TRUE 40 TIMESC minimum time scale for plots seconds per centimeter or seconds per half inch The program automatically increases the time scale above the value TIMESC as necessary to fit the plots on the paper size specified by PAPER item 39 The length available for the time axis is PAPER 2 centimeters or half inches The units of TIMESC depend on INCH item 38 The default value of TIMESC is 0 5 41 ZMIN ZMAX MAXZ word vectors minimum and maximum ordinate values for observation time history plots The ordinate axes are 4 centimeters or 4 half inches long depending on INCH item 38 If corresponding elements of 38 3 3 9 ZMIN and ZMAX are equal for any signal automatic scaling will be used for that signal The automatic scaling algorithm used does not force 0 to be included in the seale If automatic scaling is used for some signal and both the measured and com puted time histories of that signal have the same constant value for an entire maneuver item 1 the plot for
22. M and L matrices are defined by user routines MAKEM and MAKEL respectively The v and w vectors are defined by user routine MAKEVW In the following matrix definitions the standard aircraft routines use average measured quantities unless TIMVAR sec 3 3 8 20 is TRUE in which case the matrices are redefined at each time point with the current measured values The a and Bo used are obtained from a and Bn by the equations a Ta q XALF DCGFT _ pYALF c KALF V V 62 B g mM pZB _ r XB DCGFT c KB V V The variables KALF XALF YALF KB XB and ZB are defined in section 4 3 3 If these quantities are allowed to vary only the starting values are used to compute a and Bae The Bo computation does not include the multiplication by cos a gt which is needed if B is measured by a fixed vane If SHIFT sec 4 3 3 2 is TRUE DCGFT is defined as the flight center of gravity position sec 4 3 3 10 minus the wind tunnel reference center of gravity position sec 4 3 2 5 times CHORD sec 4 3 3 6 If SHIFT is FALSE DCGFT is 0 User routine UINIT defines the initial conditions and biases The initial condition is set equal to the measured observation for all states except a and B The a and B initial conditions are defined using equation 62 The biases UOFF and YOFF sec 3 1 are defined as 0 48 4 1 3 1 4 1 3 1 Longitudinal The nondimensional longitudinal matrices are AN C C 0 La he v vn 0
23. The basic program can interact with a set of user written problem specific routines to simplify the use of the program on specific systems A set of user routines for the aircraft stability and control derivative estimation problem is provided with the program 17 Key Words Suggested by Author s Maximum likelihood Parameter estimation Aircraft stability and control Computer programs System identification Aircraft flight testing 18 Distribution Statement Unclassified Unlimited STAR Category 59 22 Price 8 00 19 Security Classif of this report 20 Security Classif of this page Unclassified Unclassified 21 No of Pages 89 For sale by the National Technical Information Service Springfield Virginia 22161 NASA Langley 1980
24. be obtained by maximizing the likelihood functional shown in equation 9 This would be easy except that Ze is a complicated funciton of G This section describes a two step estimator that estimates G and the other unknowns separately When 16 1 5 state noise is not present this procedure obtains the maximum likelihood estimates With state noise the estimates obtained are very close to the maximum likelihood estimates The basic G estimator is derived by maximizing equation 9 with respect to G ignoring the dependence of l on G The solution to this maximization problem can be written explicitly se 2 E OE ORKO es To estimate G and the other unknowns the following two step algorithm is used each iteration 1 Use one iteration of the Newton Balakrishnan algorithm to estimate all of the unknowns except G 2 Use equation 24 with Z evaluated at the revised to obtain a new estimate of G Repeat steps 1 and 2 until convergence The best convergence is usually obtained by doing a few iterations with G fixed before starting the above algorithm This is because the sample residual power in the first iteration or two is often quite large and likely to give a worse estimate of G than the starting value An overrelaxation algorithm can sometimes speed the convergence of G deter mination The MMLE3 program includes an option to use logarithmic overrelaxation on the diagonal elements of Gas If the previous e
25. change from run to run For large or complicated models the use of this option can considerably simplify the input requirements 3 3 3 User Initialization Input The user initialization input contains any input required by the user routine ONCE If user routines are bypassed sec 3 3 2 the user initialization input should be omitted Subroutine ONCE should perform any initialization or input required by the user which is not to be repeated for each case Predicted derivative input is not included here see secs 3 3 4 and 3 3 5 The input required by the standard aircraft subroutine ONCE is described in section 4 3 1 3 3 4 Predicted Derivative Control Card This card controls the program call to the user routine WTIN If user routines are bypassed sec 3 3 2 this card should be omitted WTIN reads predicted derivative data from cards and writes the predicted derivative file The first three characters of the predicted derivative control card must be NEW or OLD or NO or ONL for only any other first three characters will result in an error message The entire card is printed on the output listing so the remaining 77 columns can be used for comments on the listing In addition to controlling the call to WTIN the first three characters of this card determine the value of the logical variable WI FILE WTFILE is intended to inform the user routines whether or not data are available on the predicted derivative file 30 3 3
26. direct series evaluation of geet may become computationally unstable for low sample A At NEAT rates In such cases e and its integral are first evaluated using 10 terms of the power series The desired transition matrices are then obtained by recursive applications of the formulae t 2 A pes L 2 51 Des i Aa 52 This process provides improved computational stability without significantly in creasing complexity or computer time ref 50 NEAT 0 implies direct series computation In general NEAT should be large enough so that magnitudes of the eigenvalues of A are less than 0 1 The default value of NEAT is 0 At NEAT 2 16 NOITER maximum total number of iterations The convergence criterion specified by BOUND item 17 or the error criterion specified by ERRMAX item 18 can cause iteration to terminate early otherwise NOITER iterations will be done If NOITER is 0 the parameter identification step is skipped entirely and the program computes the final time histories using the initial coefficient estimates Program dimensions limit NOITER to 49 or less The default value for NOITER is 6 34 3 3 8 17 BOUND convergence bound If for any iteration the error sum changes by less than BOUND times the error sum of the previous iteration the algorithm is assumed to have converged One of three actions is then taken depending on ITAPR item 22 and ITG item 24 status a priori is turned off G determ
27. ee ee 20 30 BASIC PROGRAM a a ll a AR a a AL 3 1 Equations of Motion Filter and Gradients a lt 22 3 2 Files A KE c l vat Ge ZO 3 2 1 Card Reader and Line Printer Files Be mks Bie ths oo t i dc od At ees sehr oa XS c 3 2 2 Input Time History File s s s a sos s ea w ee ee e ane m 246 3 2 3 Predicted Derivative File 1 eee ee ee ee 26 3 2 4 Punch File LA A A a A a ss s EN 3 2 5 Output Time History File E ra a TOD et A e o sl 3 2 6 Plot File Be abe Se 6 tes ES A cae i A on Tai 3 2 7 Internal Scratch Files bi t de oh S are e fo on ae n a v 0 ee SO V a 3 3 Input Description s sos oe 27 Syntax Check Card se Sete s c o se Ge 7 2 Ge Fy eet Gee d O User Routines Control Card Be gg Ok cis oh A to as ie ee es Se a L C 41 O User Initialization Input ir toatl id ee 099g Ge HE A Saget s25 3 Predicted Derivative Control Card iG s u ed OA TO AG See Sen tao e a ae Se Predicted Derivative Input 656 8 ee ee eee Sl Tie Card ee y dot hr Aw er Bev ch ae ae Ee a E RL User Input LHe i eda c y SE i NAMELIST INPUT St Me ies J he Me a c l ee Wh es te es Gd t c t Boe ee OE Signal Labels 4 v se eo ht AO oe ee EE 6 9096 S 3 OR ee ar ee 10 Time Cards e s s e 4 6 amp E GEO ee c t oe e 2U Ai Mattix puta 4 4 a amp l 0 166 5 s2 29287248 Se s a A a VU 12 Endease Card s
28. from file number UDATA only if CARD sec 3 3 8 2 is FALSE otherwise time history data are read from cards and file number UDATA is ignored Several details about the time history data are described in section 3 3 8 items 3 to 9 and in section 4 3 6 The input time history file is a standard FORTRAN unformatted file with one record for each time point The first four words of each record contain the time as integer hours minutes seconds and milliseconds The rest of the record contains the data words as real FORTRAN variables The default order of the channels is a q V 9a a 8 amp n x e eee longitudinal 5 phMqBpra pr65 8 8 5 N T This longitudinal y lateral lateral order can be changed with ZCHAN UCHAN and EXCHAN sec 3 3 8 7 The last record on the file should have a time greater than or equal to the ma neuver stop time sec 3 3 10 otherwise an end of file error will occur Cases need not be run in the order that they appear on the file 4 2 2 Predicted Derivative File FORTRAN file number UWT is used by the standard aircraft routines for predicted derivatives and related data If the predicted derivative control card sec 3 3 4 specifies NO this file is not used if it specifies OLD an existing file is read If the predicted derivative control card specifies NEW or ONL data are read from cards sec 4 3 2 to create a predicted derivative file If the ONL option was used the program then stops if the N
29. gradient like algorithm is generally discouraged because of its slow convergence An alternate approach to initial convergence problems is to use an a priori weighting sec 1 3 for the first few iterations 1 2 3 Inequality Constraints When state noise is present the minimization is subject to the inequality con straint of equation 17 As mentioned previously this constraint is approximated by constraining some diagonal elements of KC to be less than or equal to 1 Although simpler than the eigenvalue constraint this is still a nonlinear inequality constraint 14 1 3 Minimizing a nonlinear functional subject to a nonlinear inequality constraint is a problem in nonlinear programming The most common solution technique for such problems is to solve a series of quadratic programming problems that locally approxi mate the nonlinear problem and converge to the solution of the nonlinear problem ref 46 This approach is used by MMLE3 1 3 Inclusion of A Priori Information Predicted derivative estimates are often available from previous cases or other sources It is sometimes desirable to have the estimation algorithm consider this a priori information in addition to the new information obtained from the maneuver This can be heuristically accomplished by adding to the cost functional a quadratic penalty function for departure from the predicted values The resulting cost fuctional is J gt 2 00 z t GG 1 z
30. least one data point with time greater than or equal to the end time must be read from the input file in order for the program to detect the end of the maneuver end of file errors can result from neglecting this fact the program has no end of file checks because such checks are system dependent 3 3 11 Matrix Input This section describes the input of the data matrices Any number including 0 of the matrices listed below can be read in any order Each of the matrices has a default value which is used if that matrix is not read in Any of these defaults can be changed by user routine MATDEF unless user routines are bypassed sec 3 3 2 The default values used by the standard aircraft routines are specified in section 4 3 5 All values are reset to the defaults at the beginning of each case regardless of any values used in previous cases The same format is used for all of the matrices except for the HARD and SOFT matrices the format for HARD and SOFT is described in item 6 Each matrix is preceded by a header card with the matrix name number of rows and number of columns in A4 16 110 format The matrix name always starts in column 1 of the card The body of the matrix follows one row to a card in 8F10 format A diagonal matrix can be input in a more compact form by specifying 0 as the number of columns the diagonal elements are then read from one card in 8F10 format 1 AN BN SN RN CN DN HN EN nondimensional starting
31. mg o 1 2 1 VK V P ER A C GG P S C WeE R 1 1 1 y A ER YA ER Y asa E E ti VP a E Jg ta 9P tenn VES VE e ti 1 R Bu R S R v Yee Hier h t i PEC V g Ml Y E Ze 8 41 c ti 1 R R A GG 1 the gradient Ver is computed as the solution to the set of Lyapunov equations where 2 9 P UPA 6 wf R LA L K ER A c At 4 xs r R v A P vi 1 aq 9 c 1 P K V R 1 AP R V g FDA RT we R F RE gt KE 1 R s R R AP 1 3 1 44 45 46 47 48 49 50 29 3 2 3 2 Files The MMLE3 program uses several files for input and output The general usage of these files is discussed in this section Specific file formats are generally defined by the user and thus are not discussed here The file formats defined by the standard aircraft routines are discussed in section 4 2 To facilitate compatibility with different computer systems the following program variables are used for all file numbers The Programmer s Manual ref 38 describes how the values of these variables can be changed Variable Name Description Current Value UCARD Card reader 1 UPRINT Line printer 3 UDATA Input time history data 4 UWT Predicted derivative data 10 UPUNCH Card punch 2 UTHOUT Output time history data 9 UPLOT Plot file on some systems 13 UTI Tem
32. or w 2 2 4 46 Bd ea a ew ew A 13 Time History Cards lt ece 2 3s 44 6 a4 As a as AA coaoanoob WN a G2 G2 www G2 G2 G ww Www A w G2 G2 G2 G2 G2 G2 G2 amp G9 G2 G2 G2 ili 4 0 STANDARD AIRCRAFT ROUTINES 4 1 Equations of Motion 4 1 1 Nonlinear Six Degree of Freedom Equations 4 1 2 Uncoupled Linearized Equations 4 1 2 1 Longitudinal 4 1 2 2 Lateral directional 4 1 3 Matrix Equations 4 1 3 1 Longitudinal 4 1 3 2 Lateral directional 4 2 Files a 16 16 s l o s 4 2 1 Input Time History File 4 2 2 Predicted Derivative File 4 2 2 1 Header information 4 2 2 2 Predicted derivative data 4 2 2 3 End card Er 4 2 3 Punch File 4 2 4 Output Time History File 4 3 Input Description i 4 3 1 User Initialization Input 4 3 2 Predicted Derivative mins 1 Title card 2 NAMELIST WIND yA E 3 Angle of attack break points 4 Mach number break points 5 Param break points 6 Predicted derivative data 1 End card User Input Time History Card Input Matrix Defaults Altered Basic Program Defaults Be Eo Dee ee IS 5 0 CONCLUDING REMARKS REFERENCES INDEX OF NAMELIST VARIABLES 1V USER S MANUAL FOR MMLE3 A GENERAL FORTRAN PROGRAM FOR MAXIMUM LIKELIHOOD PARAMETER ESTIMATION Richard E Maine and Kenneth W Iliff Dryden Flight Research Center INTRODUCTION This report is a user s manual for the FORTRAN IV computer program MMLE3 a maximum likelihood parameter estim
33. program sec 3 3 8 unchanged variables are not included The matrix default changes are described in section 4 3 5 and are not repeated here Many of these defaults depend on the variables LONG and LATR sec 4 3 3 1 If LONG is TRUE the case is longitudinal if LATR is TRUE the ease is lateral directional The default value of NREC sec 3 3 8 6 is changed to 25 The defaults for ZCHAN UCHAN and EXCHAN sec 3 3 8 7 are Longitudinal ZCHAN 1 5 3 0 6 UCHAN 8 9 10 11 EXCHAN 15 16 17 18 12 14 13 3 20 21 26 27 28 29 30 31 32 33 34 35 Lateral directional ZCHAN 16 17 18 12 19 20 21 UCHAN 22 23 24 29 EXCHAN 15 1 2 5 4 14 13 3 6 7 26 27 28 29 30 31 32 33 34 35 The default value of RELAB sec 3 3 8 10 remains FALSE but the default signal labels are changed as follows Longitudinal OBSERVATIONS ALPHA Q THETA AN AX Q DOT STATES ALPHA Q THETA CONTROLS DELTA E DELTA C DELTA C1 DELTA C2 EXTRA SIGNALS Q BAR BETA P R PHI MACH ALT V P DOT R DOT RPM THRUST 8 9 0 Lateral directional OBSERVATIONS BETA P R PHI AY P DOT R DOT STATES BETA P R PHI CONTROLS DELTA A DELTA R DELTA CI DELTA C2 EXTRA SIGNALS Q BAR ALPHA Q AN THETA MACH ALT V Q DOT AX RPM THRUST The default value of USERIC see 3 3 8 26 is changed to TRUE Furthermore USERIC should not be set to FALSE or the equations of motion of the standard aircraft routines may be incorrect 9 0 CONCLUDING REMARKS A digital computer
34. program written in FORTRAN IV is available for maximum likelihood parameter estimation This program is capable of handling general linear or bilinear equations of arbitrary order with or without state noise The basic program is quite general and therefore applicable to a wide variety of problems The basic program can interact with a set of user written problem specific routines to simplify the use of the program on specific systems A set of user routines for the standard aircraft stability and control problem is provided with the program The program incorporates features that have been found useful for analysis of actual flight data at the NASA Dryden Flight Research Center based on experience with over 9000 maneuvers from 35 different aircraft Dryden Flight Research Center National Aeronautics and Space Administration Edwards Calif August 24 1979 79 10 11 REFERENCES Parameter Estimation Techniques and Applications in Aircraft Flight Testing NASA TN D 7647 1974 Gerlach O H High Accuracy Instrumentation Techniques for Nonsteady Flight Measurements Flight Test Instrumentation Vol 3 M A Perry ed Pergamon Press 1964 pp 77 100 Taylor Lawrence W Jr and Iliff Kenneth W A Modified Newton Raphson Method for Determining Stability Derivatives From Flight Data Paper presented at 2nd International Conference on Computing Methods in Opti mization Problems Sanremo Italy Sept 1968 Goodwin Graham
35. several independent variables In this case the changes implied by each constraint are added If the constraint ratio is entered as 0 it will be automatically computed as the ratio of the starting values if the independent variable has a starting value of 0 a warning message is printed and a ratio of 1 is used The maximum number of hard constraints is MAXHRD 1 sec 2 The default uses no hard constraints 7 SOFT soft constraints The input format for soft constraints is the same as that for hard constraints except that the matrix name on the header card is SOFT For a soft constraint both variables must be specified as independently varying by the variation matrices item 3 or the constraint will be ignored as immaterial in specific neither may be a dependent variable in a hard constraint The soft con straint means simply that a priori weighting will penalize the constrained variable for departure from the constraint instead of from its initial condition Thus soft con straints are immaterial if a priori weighting is not used or if the constrained variable has an a priori weighting of 0 Multiple or chained soft constraints are allowed in any combination The maximum number of soft constraints is MAXSFT 1 sec 2 The default uses no soft constraints 3 3 12 Endcase Card The end of the matrix input section is signaled by a card with END or ENDCASE starting in column 1 ENDCASE indicates that more cases follow
36. that signal will be omitted for that maneuver The default values for ZMIN and ZMAX are all 0 implying automatic scaling 42 UMIN UMAX MAXU word vectors minimum and maximum ordinate values for control time history plots The ordinate axes are 4 centimeters or 4 half inches long depending on INCH item 38 If corresponding elements of UMIN and UMAX are equal for any signal automatic scaling will be used for that signal The automatic scaling algorithm used forces 0 to be included in the scale If automatic scaling is used for some signal and the time history of that signal is 0 for an entire maneuver item 1 the plot for that signal will be omitted for that maneuver The default values for UMIN and UMAX are all 0 implying automatic scaling 43 XMIN XMAX MAXX word vectors minimum and maximum ordinate values for state time history plots The interpretations of XMIN and XMAX are identical to those of UMIN and UMAX item 42 The default values of XMIN and XMAX are all 0 implying automatic scaling 44 EXMIN EXMAX LEX word vectors minimum and maximum ordinate values for extra signal time history plots The interpretations of EXMIN and EXMAX are iden tical to those of UMIN and UMAX item 42 The default values of EXMIN and EXMAX are all 0 implying automatic scaling 45 PUNCH logical punched output option If PUNCH is TRUE user routine OUTPUN will be called to write data on the punch file sec 3 2 4 If NOIT
37. used the same basic algorithm as that of reference 6 but incor porated many features found useful for routinely processing large amounts of flight Stability and control data The program was widely circulated among industry and government agencies Extensive worldwide experience has been obtained using maximum likelihood estimation on actual flight data refs 13 to 26 The NASA Dryden Flight Research Center alone has analyzed over 5000 maneuvers from 35 different aircraft with these techniques Several studies have been made of the factors to be considered in care fully applying maximum likelihood estimation to flight data refs 27 to 29 This background of experience and careful study has resulted in maximum likelihood estimation becoming recognized as the leading technique for estimation of aircraft stability and control derivatives ref 30 Reference 31 discusses the Dryden Flight Research Center s approach to the gathering and analysis of flight data for maximum likelihood estimation The success of maximum likelihood estimation of standard aircraft stability and control derivatives has prompted interest in expanding the scope of its application Most of the work with actual flight data has been done in stabilized flight at low to moderate angles of attack a standard linear set of uncoupled equations with two degrees of freedom has been used for longitudinal data or three degrees of freedom for lateral directional data Investigation o
38. 1 2 ti 2 s ws 50 21 where So is the vector of a priori values and W is the a priori weighting matrix The maximum likelihood method with the quadratic penalty function can be interpreted as maximizing the unconditional probability density instead of the conditional prob ability density function The a priori weighting option can be used to improve the convergence of ill conditioned cases The resulting estimates however are biased When a priori weighting is used to aid convergence it is recommended that the a priori weighting be removed after initial convergence and the last few iterations run without a priori weighting The converged values from the algorithm with a priori weighting will usually be close enough to the maximum likelihood estimates that the program will rapidly converge to the asymptotically unbiased maximum likelihood estimates if the a priori weighting is then removed When a priori weighting is used without this final step care should always be taken that the resulting bias in the estimates does not lead the user to false conclusions 1 4 Cram r Rao Bound With any parameter estimation method it is important to have a measure of the accuracy of the estimates Reference 47 discusses the evaluation of the accuracy of the estimates including detailed treatment of the subjects briefly mentioned here In the absence of modeling error or bias the scatter of the estimates is a reasonable indication of
39. 5 If the first three characters OLD or NO are used subroutine WTIN will not be called With the NO option WTEILE will be set to FALSE indicating no data on the predicted derivative file With the OLD option WTFILE will be set to TRUE indicating data available on the predicted derivative file created by a previous job If the first three characters NEW are used WTIN will be called and WTFILE will be set to TRUE The user subroutine WTIN should create the predicted derivative file If the first three characters ONL are used WTIN will be the last subroutine called After this call the program will stop This option is used if WTIN is to create a predicted derivative file to be saved for later use but no cases are to be run at present 3 3 9 Predicted Derivative Input The predicted derivative input contains any input required by the user routine WTIN The input required by the standard aircraft routine WTIN is described in section 4 3 2 If the user routines control card sec 3 3 2 specifies NO WTIN is not called so the predicted derivative input should be omitted If the user routines control card is anything else the call to WTIN is controlled by the predicted derivative control card sec 3 3 4 If the predicted derivative control card specifies NO or OLD WTIN is not called and the predicted derivative input should be omitted if the predicted derivative control card is NEW or ONL the input is included 3 3 6 Title Ca
40. 9 Gupta Narendra K Hall W Earl Jr and Trankle Thomas L Advanced Methods of Model Structure Determination From Test Data AIAA Paper 77 1170 Aug 1977 Jazwinski Andrew H Stochastic Processes and Filtering Theory Academic Press 1970 Luenberger David G Optimization by Vector Space Methods Wiley 1969 Dixon L C W Nonlinear Optimisation Crane Russak and Co Inc 1972 Maine Richard E and Iliff Kenneth W Estimation of the Accuracy of Dynamic Flight Determined Coefficients AIAA Paper 80 0171 Jan 1980 Potter James E Matrix Quadratic Solutions J SIAM Appl Math vol 14 May 1966 pp 496 501 Wilkinson J H The Algebraic Eigenvalue Problem Clarendon Press London 1978 Moler Cleve and Van Loan Charles Nineteen Dubious Ways to Compute the Exponential of a Matrix SIAM Review vol 20 no 4 Oct 1978 pp 801 836 Gainer Thomas G and Hoffman Sherwood Summary of Transformation Equa tions and Equations of Motion Used in Free Flight and Wind Tunnel Data Reduction and Analysis NASA SP 3070 1972 83 INDEX OF NAMELIST VARIABLES The variables in the three NAMELISTS of MMLE3 are listed alphabetically below The numbers in parentheses refer to the sections and item numbers where definitions are found NAMELIST INPUT Section 3 3 8 BOUND 17 CARD 2 DFAC 21 DIAGG 25 ERRMAX 18 ERRTH 32 EXBIAS 9 EXCHAN 7 EXMAX 44 EXMIN 44 EXSCAL 8 FREQCR 23
41. C and Payne Robert L Dynamic System Identification Experiment Design and Data Analysis Academic Press 1977 Balakrishnan A V Communication Theory McGraw Hill Book Co Inc c 1968 Taylor Lawrence W Jr and Iliff Kenneth W Systems Identification Using a Modified Newton Raphson Method A FORTRAN Program NASA TN D 6734 1972 lliff Kenneth W and Taylor Lawrence W Jr Determination of Stability Derivatives From Flight Data Using a Newton Raphson Minimization Technique NASA TN D 6579 1972 Larson Duane B Identification of Parameters by the Method of Quasilineariza tion CAL Rep 164 Cornell Aeronautical Lab Inc May 1968 Grove Randall D Bowles Roland L and Mayhew Stanley C A Procedure for Estimating Stability and Control Parameters From Flight Test Data by Using Maximum Likelihood Methods Employing a Real Time Digital System NASA TN D 6735 1972 Ross A Jean and Foster G W Fortran Programs for the Determination of Aerodynamic Derivatives From Transient Longitudinal or Lateral Responses of Aircraft C P No 1344 British A R C 1976 lliff Kenneth W and Maine Richard E Practical Aspects of Using a Maximum Likelihood Estimation Method to Extract Stability and Control Derivatives From Flight Data NASA TN D 8209 1976 12 Maine Richard E and Iliff Kenneth W A FORTRAN Program for Determining 13 80 Aircraft Stability and Control Derivatives From Flight Data NASA TN
42. END is used for the last case to be processed 3 3 13 Time History Cards If the user routine READTH reads the time history data from cards these cards should be placed in this location Various aspects of the time history data are de scribed in section 3 3 8 items 2 to 9 The card format required by the standard aircraft routine READTH is described in section 4 3 4 All of the time history data cards should be read by READTH if more cases follow otherwise the first unread card will be interpreted as the title card for the next case 4 0 STANDARD AIRCRAFT ROUTINES This section describes a particular set of user routines referred to as the standard aircraft routines These routines are intended for the aircraft longitudinal or lateral directional stability and control problem with no turbulence 42 4 1 4 1 Equations of Motion 4 1 1 Nonlinear Six Degree of Freedom Equations This section presents the nonlinear aircraft equations of motion from which the linearized equations used in MMLE3 are derived The coupled nonlinear equations are presented first These equations assume a rigid vehicle and a flat nonrotating earth The time rate of change of mass is assumed negligible and fuel sloshing effects are ignored No small angle approximations are used but the absolute values of B and must be less than 90 The aircraft velocity must not be 0 No symmetry assumptions are made Engine inertia and thrust terms are included assum
43. ER item 16 is 0 PUNCH is forced to FALSE If the ERRMAX limit Gtem 18 is exceeded OUTPUN will be skipped regardless of the value of PUNCH The format of the data written is defined by the user The default value of PUNCH is FALSE 3 3 9 Signal Labels The cards described in this section contain the labels to be used for the obser vations states controls and extra signals Signal labels are included only if RELAB sec 3 3 8 10 is set to TRUE If RELAB is FALSE the default labels used are Zi Xi Ui or EXi for the ith observation state control or extra signal respectively These default labels are changed by the standard aircraft routines The default labels defined by the standard aircraft routines are specified in section 4 3 6 If RELAB is TRUE all of the labels must be read there is no provision for re defining only a portion of the labels The labels are read left justified in fields of 10 up to 8 labels per card Each label is limited to eight characters The labels for MAXZ observations are read from the first card s MAXX states from the next card s MAXU controls from the following card s and LEX extra signals from the last card s 39 3 3 10 3 3 10 Time Cards One time card is required for each of the NCASE sec 3 3 8 1 maneuvers to be analyzed The time cards contain the start and end times for each maneuver expressed as hours minutes seconds and milliseconds in the format 2 312 13 1X At
44. EW option was used execution continues using the 26 4 2 2 1 predicted derivative file just created The predicted derivative file is composed of 80 column card images The contents are described below The units must be consistent with METRIC sec 4 3 3 3 4 2 2 1 Header information Card 1 Card 2 Card 3 Card 4 Card 5 Card 6 Card 7 Card 8 Card 9 Card 10 80 column title card identifying the data set Reference area reference span reference chord reference center of gravity and axis system The axis system is either STAB or BODY The data format is 3F10 3 F10 5 A4 Angle of attack and angle of sideslip flow amplification factors KALF and KB respectively sec 4 3 3 11 in 2F10 3 format Instrument locations XALF XB XAN XAX and XAY sec 4 3 3 12 1n 5F10 3 format Instrument location YALF YB YAN YAX and YAY sec 4 3 3 13 in 5F10 3 format Instrument locations ZALF ZB ZAN ZAX and ZAY sec 4 3 3 14 in 5F10 3 format Number of angle of attack Mach number and extra parameter break points in the predicted derivative data tables The format of this eard is 3110 Angle of attack break points in 8G10 3 format The break points are continued to further cards in the same format if required Mach number break points in 8G10 3 format The break points are continued to further cards in the same format if required Extra parameter break p
45. If MZ is specified as negative it will be obtained from the number of rows of the GGI matrix sec 3 3 11 The standard aircraft routines define a default number of rows of GGI sec 4 3 5 The default value for MZ is 1 33 3 3 8 13 MU length of the control vector All matrix dimensions will be forced to conform with MU controls The maximum value for MU is given by the variable MAXU sec 2 currently 4 The minimum value of MU is 1 If MU is specified as negative it will be obtained from the number of columns of the BN matrix sec 3 3 11 If this number is 0 MU will be set to the largest column number that contains a nonzero element in BN DN BV or DV The standard aircraft routines define a default num ber of columns of BN The default value for MU is 1 14 MB length of the bias vector All matrix dimensions will be forced to con form with MB elements in the bias vector sec 3 1 The maximum value for MB is given by the variable MAXB sec 2 currently 4 The minimum value of MB 1s 1 If MB is specified as negative it will be obtained from the number of columns of the SN matrix sec 3 3 11 If this number is 0 MB will be set to the minimum of NCASE item 1 or MAXB sec 2 MB is forced to be less than or equal to NCASE The default value for MB is 1 15 NEAT number of time reductions in the computation of the transition matrix and its integral The transition matrix is eet where in general Ae Rta The
46. M and L dimensionalization matrices will be recomputed at each time point The time varying option is only partially implemented for the state noise case in that K and its gradients are considered time invariant time varying state noise cases should therefore be critically examined If user routines are bypassed sec 3 3 2 the time varying option is meaningless and will be overridden if attempted USERIC item 26 is sometimes required in time varying cases depending on the system model Use of the time varying option should be avoided if possible because it increases the com puter time used by approximately a factor of three The default for TIMVAR is FALSE 21 DFAC ITDFAC diagonal convergence factor option The diagonal elements of the second gradient matrix will be multiplied by DFAC before computing parameter changes for ITDFAC iterations This option starts when a priori is turned off if ITAPR item 22 is not 0 If ITAPR is 0 the DFAC option starts on iteration 2 iteration 1 if FULL1 item 19 is TRUE NOITER item 16 should usually be larger than ITDFAC so that the final iterations will be done without the diagonal convergence factor see sec 1 2 2 This option sometimes improves the convergence properties of marginally stable maneuvers or cases where the starting estimates are far from the correct values Convergence using this option tends to be monotone but very slow Because of the slow convergence the use of this op
47. NASA Technical Paper 1563 User s Manual for MMLE3 a General FORTRAN Program for Maximum Likelihood Parameter Estimation Richard E Maine and Kerineth W Iliff NOVEMBER 1980 NASA NASA Technical Paper 1563 User s Manual for MMLE3 a General FORTRAN Program for Maximum Likelihood Parameter Estimation Richard E Maine and Kenneth W Iliff Dryden Flight Research Center Edwards California NASA National Aeronautics and Space Administration Scientific and Technical Information Branch 1980 CONTENTS Page INTRODUCTION c 404 2 44 Oe Bee he AAA Se eee eee SYMBOLS e da Sb kh ROE Se ER EER EOS eS ee dc OS 1 0 PARAMETER ESTIMATION THEORY ee eee eee 8 1 1 Maximum Likelihood Estimation e eee ees 9 1 1 1 Measurement Noise Only A i dc AO DEA s4 Ve 1 1 2 Measurement Noise and State Noise A A E 1 2 Minimization of the Cost Functional 24 13 1 2 1 Newton Balakrishnan Algorithm 13 1 2 2 Improvement of Convergence Properties 14 1 2 3 Inequality ConstralMts s y mas cea 9 3 wee s t eee A 1 3 Inclusion of A Priori Information 6 ew ee ee 1 4 Cram r Rao Bound ur an dE Aa e2 a a 0 1 5 Estimation of Residual POWER A A E AS Se Se sec E 1 6 Solution of the Riccati Equation s s 18 2 0 MMLE3 PROGRAM STRUCTURE e ee
48. RMAX item 18 or PLTMAX item 34 bounds are exceeded The default value of PLOTEM is TRUE 34 PLTMAX maximum error for plotting If the final error sum is greater than PLTMAX plots will not be produced regardless of PLOTEM item 33 In this case the fits will probably be uninformative and may exceed reasonable plotter limits If ERRTH item 32 is TRUE and PLTMAX is exceeded measured time histories will be printed as a diagnostic aid The default value for PLTMAX is 10 35 XPLOT MAXX word logical vector states to be plotted For each element of XPLOT that has the value TRUE the corresponding corrected state estimate time history will be plotted provided PLOTEM item 33 is TRUE The default values for XPLOT are all FALSE 36 NUPLT number of controls to be plotted The first NUPLT control time histories will be plotted if PLOTEM item 33 is TRUE The program limits the maximum value of NUPLT to MAXU sec 2 currently 4 The minimum value is 0 The default value for NUPLT is MAXU 37 NEXPLT number of extra signals to be plotted The first NEXPLT extra signal time histories will be plotted if PLOTEM item 33 is TRUE The program limits the maximum value of NEXPLT to LEX sec 2 currently 20 The minimum value is 0 The default value of NEXPLT is 0 38 INCH logical PLTFAC plot sizing eontrols If PLTFAC is 1 and INCH is TRUE plots will be sized for inch grid paper If PLTFAC is 1 and INCH is
49. alue of ALPHA is 999 18 THETA average pitch angle deg If THETA is set to 999 a value is obtained from the average of the measured pitch angle time history The default value of THETA is 999 19 PHI average bank angle deg If PHI is set to 999 a value is obtained from the average of the measured bank angle time history The default value of PHI is 999 20 MACH real average Mach number If MACH is set to 0 a value is obtained from the average of the measured Mach number time history The default value of MACH is 0 21 PARAM extra identifying parameter This parameter is used in the predicted derivative table lookup sec 4 3 2 6 it is also included on the punched output sec 4 2 4 The default value of PARAM is 0 22 UVAR MAXU word integer vector vector of control flags UVAR is a convenient means of defining the BV and DV matrices for the standard aircraft routines without reading in the entire matrices UVAR affects the default values of the BV and DV matrices therefore if BV or DV matrices are read from cards the matrices read in will override those defined by UVAR For each element of UVAR that is set to 1 the standard derivatives for the cor responding control are allowed to vary for the elements of UVAR that are 0 the control derivatives are fixed The C y C Q and Ch control derivatives are standard for lateral cases the C C A and Cn derivatives are standard for longitudinal N ca
50. anges drastically It is no longer possible to exactly identify the values of the unknown parameters rather the values must be estimated by some statistical criterion At this point of development the problem should precisely be referred to as parameter estimation rather than para meter identification although both terms are often used in the literature Since the resulting parameter values can no longer be considered perfect an entire branch of theory ref 39 is introduced to answer questions about the accuracy of the estimates For discrete time observations the theory of estimation in the presence of measure ment noise is relatively straightforward requiring only basic probability The second complication of real systems is the presence of state noise also re ferred to as process noise or input noise State noise is the random excitation of the system from unmeasured sources the standard example for the aircraft stability and control problem being atmospheric turbulence If state noise is present but measure ment noise is neglected the standard line of analysis results in the regression algorithm ref 40 The MMLE3 program does not currently include regression When both state and measurement noise are considered the algorithm is more complicated than in the state noise only or measurement noise only case However with reasonable care in the implementation the computer costs need not be much higher than in the measurement noise
51. ased on the assumptions of Gaussian white noise Reference 47 discusses the problems of applying this equation to actual data with band limited noise An approximate correction for the effects of band limited noise is implemented in MMLE3 The residuals 7 2 are filtered with a first order filter the break frequency of which can be defined by the user sec 3 3 8 23 The cost functional eq 10 is then evaluated using the filtered residuals in place of the unfiltered residuals and is multiplied by the Nyquist frequency ref 47 divided by the filter break frequency For white residuals the result should approximately equal the original cost functional value which should approximately equal the length of the observation vector MMLE3 multiplies the variances of equation 23 by the adjusted filtered cost functional divided by the length of the observation vector The MMLE3 program prints both the Cram r Rao bounds and the estimated cor relations ref 47 The Cram r Rao bounds are the best available indicators of overall accuracy The smaller the Cram r Rao bound the more confidence that can be placed in the estimates The estimated correlations can sometimes be used to help find the sources of problems indicated by the Cram r Rao bounds 1 5 Estimation of Residual Power The previous discussions have largely neglected the question of estimating G Theoretically maximum likelihood estimates of G as well as the other unknowns can
52. ation program capable of handling general bilinear dynamic equations of arbitrary order with or without state noise Contin uous time system equations with discrete sampled observations are used Many of the advances in practical application of estimation methodology have been in the area of aircraft stability and control derivatives ref 1 The well defined dynamic model the high quality automated data acquisition systems ref 2 and the presence of analysts oriented towards computer techniques have contributed to this trend In 1966 NASA Dryden Flight Research Center began investigating ref 3 the use of maximum likelihood estimation ref 4 based on the theoretical work of Balakrishnan ref 5 Taylor and Iliff ref 6 refer to the computer prog ram developed as the Newton Raphson program because it uses a modified Newton Raphson algorithm Newton Balakrishnan to implement the maximum likelihood estimation Reference 7 compares the maximum likelihood method to previously used techniques of estimating stability and control derivatives There have been numerous parallel and subsequent developments refs 8 to 10 of techniques identical or similar to the Iliff Taylor Balakrishnan approach The Newton Raphson program underwent a gradual evolution based on expe rience in application to flight data ref 11 Reference 12 documents the resulting program released in 1973 This program named MMLE modified maximum likeli hood estimator
53. ave nondimen sional forms In some cases the same unknown parameter can appear in two or more matrix locations This results in hard constraints that must be satisfied as the matrices are updated The handling of such constraints is discussed in section 3 3 11 6 Soft constraints those implemented by a penalty function approach are discussed in section 3 3 11 7 Two general restrictions are necessary R must be nonsingular and ER r must be 0 for all values of the unknown parameters in some neighborhood of the true value The restriction on ER T insures that the state noise and measurement noise are independent it is not needed if the x in the Ex term is considered to be the x computed with F 0 The residual covariance matrix see below must be non singular The steady state Riccati covariance matrix P is defined by 1 r 14 p p R la p ER lA o GE ers C P RFR 1p 0 38 with the requirement that P be symmetric The steady state Kalman filter gain matrix is defined by K P ER_ la c ery 39 When the system matrices are time varying P K and their gradients are computed using average values in the system matrices All of the other expressions are re evaluated at each time point Naturally if F 0 the no state noise case then P 0and K 0 23 3 1 The filter equations for the state estimation are as follows Prediction step tia OE 1 R Bu t p YR 44 UR Cv ty 2 t 1
54. c procedures MMLE3 does not incorporate automatic model structure determination 1 1 Maximum Likelihood Estimation This section describes the maximum likelihood estimation method used by MMLE3 Maximum likelihood estimation has many desirable statistical characteristics for example it yields asymptotically unbiased consistent and efficient estimates refs 4 5 and 36 The maximum likelihood estimates are obtained by choosing the vector of un known parameters to maximize the likelihood ratio p 2 8 where Z is the measured response of the system The crucial element of the theory is the definition of the likelihood ratio This will be discussed separately for the cases with measure ment noise only and with both state and measurement noise 1 1 1 Measurement Noise Only When only measurement noise is present the likelihood ratio is easily computed even for nonlinear systems The MMLE3 program is written specifically for linear or bilinear sec 3 1 systems but the theory is valid in general Most nonlinear problems can be reformulated as bilinear systems so the restriction on MMLE3 is not severe The actual system is assumed to be described by x t fan u t E He Z g e t u t Etna Ot 1 1x24 where x state vector u 1nput vector Z measured observation vector n measurement noise vector e true value of vector of unknown parameters The measurement z is obtained at the discrete time poin
55. dicted moment derivatives if SHIFT item 2 is TRUE The default value of CG is 0 25 11 KALF KB real flow amplification factors for angle of attack and angle of sideslip respectively The computed observations for these angles are multi plied by KALF or KB These multiplication factors are used to linearly model local flow and upwash effects The default values of KALF and KB are obtained from the predicted derivative file sees 4 2 2 and 4 3 2 2 8 if present Other wise the defaults are 1 12 XALF XB XAN XAX XAY longitudinal locations of the angle of attack angle of sideslip normal acceleration longitudinal acceleration and lateral accel eration sensors respectively ft or m These locations are all given as offsets from the center of gravity positive for sensors forward of the center of gravity The center of gravity location referenced is the wind tunnel reference center of gravity if SHIFT item 2 is TRUE or the actual flight center of gravity if SHIFT is FALSE The units used for the sensor locations depend on METRIC item 3 The default values are obtained from the predicted derivative file secs 4 2 2 and 4 3 2 2 9 if present otherwise the defaults are 0 13 YALF YB YAN YAX YAY lateral locations of the angle of attack angle of sideslip normal acceleration longitudinal acceleration and lateral acceleration sensors respectively ft or m These locations are all given as offsets f
56. ds in each input record This variable is passed to user routine READTH to define the number of data words required at each time point The maximum value of NREC is limited by program dimensions to 100 The normal default value for NREC is given by the variable LORD sec 2 currently 32 This default is changed by the standard aircraft routines 7 ZCHAN UCHAN EXCHAN MAXZ MAXU and LEX word integer vectors respectively channel numbers of the observations controls and extra Signals respectively ZCHAN UCHAN and EXCHAN specify the channels of the signals in the record returned from user routine READTH The values normally lie between 1 and NREC item 6 Values greater than NREC but less than or equal to the dimension limit of 100 can be used to obtain data values of 0 The normal default values are increasing integers starting from 1 in the vector formed by concatenating ZCHAN UCHAN and EXCHAN thus the normal default values depend on the lengths of these vectors The default is changed by the standard aircraft routines 8 ZSCALE USCALE EXSCAL MAXZ MAXU and LEX word vectors respec tively scale factors for observations controls and extra signals respectively The measured data are multiplied by corresponding elements of ZSCALE USCALE or EXSCAL to compensate for known scaling or sign errors This multiplication is done before any other operations on the data The default values are all 1 9 ZBIAS UBIAS EXBIAS
57. e variable initial conditions are often equivalent to combinations of variables in the SN and HN matrices resulting in identifiability problems The default values for VARIC are all FALSE 28 TEST logical extra output for debugging If TEST is TRUE extra out put useful for debugging purposes is printed Internal location maps are printed below the list of unknowns and several matrices are printed each iteration The matrices printed include the following A B S R C D H E P K RIA Ra RIB R YB RIS R Ss ERIAC ER lA C ERIBD ER B D ERISH ER H R JAAt PHI e PSI Y PSIB WR B PSIS wR S and DK Y K where At pl S y f e sds The SUM matrix requires special mention The lower triangular 0 part of SUM contains the second gradient of the cost function The upper triangular part of this symmetric matrix is not stored The first gradient is stored and printed as an additional column of the SUM matrix The gradients printed do not include terms attributable to a priori and are not multiplied by DFAC item 21 The default value of TEST is FALSE 29 PRINTI logical option to print input time history If PRINTI is TRUE the measured time histories of the observations controls and extra signals are printed All scale factors item 8 biases item 9 and modifications by user routine THMOD are applied before the data are printed The default value of PRINTI is FALSE 30 PRINTO logica
58. equiv alent to allowing bilinear equations The unknown coefficients however cannot be treated as time functions The time varying system matrices are related to the time invariant unknowns by equations such as A AM t X eat ot t 37 22 Note that the operation used is element by element multiplication not matrix multi plication The time invariant unknown coefficients are contained in the AN matrix The form of this equation is motivated by the aircraft equations of motion sec 4 where the unknown nondimensional coefficients are multiplied by known time functions to get the dimensional coefficients We will refer to AN as a nondimensional matrix AM and AL as dimensionalization matrices and A as a dimensional matrix The use of this terminology based on the aircraft equations does not imply any program requirement that the elements of AN actually be nondimensional AM and AL are arbitrary known functions of time defined by the user routines Variations of AM and AL with time should be continuous and slow compared with the system response time these matrices are often time invariant which results in a considerable saving of computer time If the user routines are not called each element of AM t is de fined as 1 and AL t as 0 so that the A and AN matrices are equivalent In this same manner the B S R C D H and E matrices are defined in terms of time invariant nondimensional matrices The F and G matrices do not h
59. f SHIFT is TRUE 3 METRIC logical metric units option If METRIC is TRUE SI MKS units are used for all data otherwise U S Customary EGS units are used The default value of METRIC is FALSE 64 4 3 3 4 AREA reference wing area cet or m The units of AREA depend on METRIC item 3 The default value is obtained from the predicted derivative file secs 4 2 2 and 4 3 2 2 7 if present otherwise the default is 0 5 SPAN reference span ft or m The units of SPAN depend on METRIC item 3 The default value is obtained from the predicted derivative file secs 4 2 2 and 4 3 2 2 7 if present otherwise the default is 0 6 CHORD reference chord ft or m The units of CHORD depend on METRIC item 3 The default value is obtained from the predicted derivative file secs 4 2 2 and 4 3 2 2 7 if present otherwise the default is 0 7 W aircraft gross weight lb or N The units of W depend on METRIC item 3 The default value is 0 8 IX IY IZ IXZ real aircraft moments of inertia slug ft or kg m The units depend on METRIC item 3 The default values are all 0 9 IXE real X axis moment of inertia of rotating mass of the engine The units depend on METRIC item 3 The default value is 0 10 CG aircraft center of gravity in fraction of the reference chord The center of gravity is used to adjust the longitudinal instrument offsets item 12 and the pre
60. f the basic program given in section 3 3 still applies when the standard aircraft routines are used The subsections below supplement section 3 3 with details applicable to the standard aircraft routines In several cases the standard routines affect the card input of the basic routines for instance they change default values such effects are also described here 4 3 1 User Initialization Input The standard aircraft user initialization input consists of default GGI and F matrices The matrices are read in standard matrix format sec 3 3 11 except that the matrix names are five and seven characters long The matrix names used are LONGGGI LATRGGI LONGF and LATRE for the longitudinal and lateral directional matrices If not read in the defaults for LONGF and LATRF are 0 matrices and the defaults for GGI are diagonal matrices of size 5 by 3 with the following values LONGGGI 10 60 30 200 200 LATRGGI 150 0 5 300 10 5000 The last card of the user initialization input should have END starting in column 1 If no F or GGI matrices are to be read in this can be the only card 4 3 2 Predicted Derivative Input The standard aircraft version of the predicted derivative input sec 3 3 5 is described herein The predicted derivative control card sec 3 3 4 should specify NO if predicted derivatives are not used or OLD if the predicted derivative file is already available in either of these cases the predicted derivati
61. f the second parameter Finally if the derivative depends on all three parameters the two parameter form is repeated for each break point of the third parameter Two examples should aid understanding of the data table organization For both examples we assume two angle of attack break points denoted by Al and A2 two Mach number break points denoted by M1 and M2 and two param break points denoted by P1 and P2 Example 1 Column 1 11 21 31 Y CLDA LATR BN 2 1 Header card Al A2 ee a A aa M1 Data card 1 P1 M2 2 Data table Data card M1 Data card 3 P2 M2 Data card 4 Example 2 Column 1 11 21 31 CMQ LONG AN 2 2 PA Header card P1 P2 A A N ngat E amena A A Al Data card 1 Data table A2 Data card 2 Although the predicted derivative file can contain values for any of the non dimensional matrix locations the following coefficients are those normally used For a longitudinal case AN RN BN SN 0 0 0 0 0 0 0 0 C C 0 C 0 0 C C m m Me m m a q C 0 0 0 0 0 0 0 62 4 3 2 6 CN DN HN 0 0 0 0 0 0 0 0 0 0 0 0 90 0 0 Ce Ge 0 ae on a q 9 C C 0 Ca Ay Ag Ca 9 Note that total trimmed CN and Ca are required in the HN matrix See section 4 3 5 1 for further discussion of this point If the predicted longitudinal derivatives are in stability axes secs 4 3 2 2 1 and 4 2 2 then Cr and Cp derivatives should be placed in the fourth and fifth rows of the
62. f unusual aircraft ref 32 and extreme flight conditions requires more flexibility in the choice of equations of motion Application to problems such as performance estimation ref 33 or structural mode identification ref 29 also requires this flexibility Nonaircraft applications have also been considered The MMLE program developed primarily for the standard stability and control equations lacks the versatility to conveniently handle all of the new investigations Kirsten ref 34 notes that MMLE requires the dynamic pressure to remain constant during a maneuver and Park ref 35 says that it is too special ized to analyze structural responses Neither MMLE nor any of the other widely available programs combine sufficient generality with those features DFRC has found useful in routinely analyzing actual data In response to the requirements for a more versatile maximum likelihood estima tion program DFRC began the development of the MMLE3 computer program MMLE3 was designed to handle a general set of linear or bilinear sec 3 1 dynamic equations of arbitrary order Although greatly refined in application this estimation problem is solved using the same theoretical tools developed by Balakrishnan in 1966 ref 5 However the final capability added to MMLE3 estimation in the presence of state noise turbulence required further theoretical development ref 36 Iliff ref 37 demonstrated its application to flight data in turbule
63. he basic program structure to be very general while retaining the simplicity of input possible in programs designed for specific systems The use of the user routines to simplify the required input deck is most desirable when large numbers of cases are being processed The user routines also can be coded to perform computations or operations unique to a specific application Examples of this function include reading data from special formats normalizing or correcting the estimates to standard conditions and punching the results in a form suitable to auxiliary programs for purposes such as plotting esti mates or updating simulators The following is a list of the user subroutines and a very brief description of their functions Naturally lower level subroutines called from these routines can be defined the standard aircraft routines include the lower level subroutines INTERP WTDEF and WTTRAN User Subroutines Purpose ONCE Initialization WTIN Predicted derivative input USERIN User input READTH Read time histories THMOD Modify time history data AVERAG Access averages MATDEF Matrix defaults UINIT Define initial conditions 20 3 0 MAKEL Define dimensionalization addition L matrices MAKEM Define dimensionalization ratio M matrices MAKEVW Define v and w vectors THOUT Time history output OUTPUN Punched output TITPLT Plot titling The basic program forms the structure within which the user routines fit there fore the descri
64. ignores the TAPE variable The default values are FALSE for CARD and TRUE for TAPE 3 THIN integer thinning factor for input time histories If THIN is 1 every time point is used If THIN is 2 every second time point is used and so forth The default value is 1 4 SPS sampling rate of the time histories before thinning samples per second The sampling rate of the thinned data is obtained by dividing SPS by THIN item 3 If SPS is 0 the sampling rate of the thinned data will be determined from the times of the thinned data using the following algorithm The times of the first two thinned data points are subtracted and the difference is rounded to the nearest 5 milliseconds The reciprocal of this rounded difference is used for the thinned sampling rate The rounding process corrects for slippages of up to 2 milliseconds in the recorded times as long as the actual thinned sample interval is a multiple of 5 milliseconds If the thinned sample interval is not a multiple of 5 milliseconds this algorithm will compute an incorrect sampling rate and thus SPS should be specified The default value of SPS is 0 5 MAXREC maximum number of input records If user routine READTH is called more than MAXREC times for a single case the program will terminate with an error message This feature prevents infinite loops caused by bugs in READTH The default value is 100 000 32 3 3 8 6 NREC number of data words exclusive of time wor
65. ilter gain matrix flow amplification factors for angle of attack and angle of sideslip Mach number mass kg slug number of observations engine speed positive clockwise viewed from rear rpm number of time points state noise vector Riccati covariance matrix roll rate deg sec likelihood ratio pitch rate deg sec dynamic pressure N m2 Ibt ft state equation matrix degrees per radian 57 2958 relaxation factor yaw rate deg sec state equation matrix reference area m cet thrust N bf total time sec transformation matrix time sec control vector velocity m sec ft sec forcing function in state equation a priori weighting matrix forcing function in observation equation state vector longitudinal instrument offsets from center of gravity m ft lateral instrument offsets from center of gravity m ft time history of measured response measured observation vector vertical instrument offsets from center of gravity m ft bias vector angle of attack deg angle of sideslip deg time interval sec control deflection deg aileron deflection deg extra control deflection elevator deflection deg rudder deflection deg extra controls measurement noise vector Superscripts Subscripts CG CG fit ref pitch angle deg vector of unknowns vector of a priori values transition matrix bank angle deg integral of the transition matrix gradient row vector
66. ination is started or iteration is stopped Convergence may occur several times in one case for different a priori and G determination conditions If the FULL1 option item 19 is FALSE the bound will not be checked on the first iteration BOUND is also not checked on iterations that use DFAC item 21 The default value of BOUND is 0 001 18 ERRMAX maximum allowable error sum If at any time the error sum becomes greater than ERRMAX it is assumed that the algorithm is not converging properly Iteration will then terminate and depending on ERRTH item 32 the measured time histories might be printed to provide debugging clues ERRMAX Should be less than the largest single precision floating point number defined on the computer divided by the total number of time points in the case If ERRMAX is less than PLTMAX item 34 it will be set equal to PLTMAX The default value of ERRMAX is 102 19 FULL1 logical full first iteration option If FULL1 is FALSE only initial conditions and B S D and H terms are determined in the first iteration This eliminates convergence problems caused by large initial bias errors Since these terms affect the response in a linear manner the Newton Balakrishnan algorithm converges in one iteration if they are the only unknowns If FULL1 is true the first iteration is full normal iteration The default value of FULL1 is FALSE 20 TIMVAR logical time varying option If TIMVAR is true the
67. ing the engine alignment and thrust vector are along the X axis The equa tions are written in body axes referenced to the center of gravity More complete forms of the equations of motion and transformations are found in reference 51 All angles are in degrees The V equation is not included mm mV cos B Qe 11 q tan B p cos a r sin a RC T ox cos 8 cos q cos a sin 9 sin a sina V cos B mg 3 p sin a 3 Go 2 i B p sina r cos a R cos B C J cos 8 sin R sin B EA 2 cos 8 cos y sin a sin 9 cos a m cos a 53 Ep c see a 1 _ Mi aa pl Gey e qsbKC o Re ar 1 1 a r Jiz pq FPI y p cata 7 a 1 D re al oz Bly ISCRC m 2 rp t 1 r p Y E Pal 6N o F ZE 1 E 2 2 M _ rl pl alos qsb C y 23 pa 1 p q Vary pel ari GNL q cos rsin p r cos tan 9 q sin q tan 9 4 1 2 Uncoupled Linearized Equations In order to divide the equations into longitudinal and lateral directional sets small angle approximations are needed for B The sin f term in the f equation is ignored and cos f is replaced by 1 in the and B equations The tan f in the equation is replaced by B 4 it is not ignored because it is multiplied by p which can be quite large Symmetry about the XZ plane is assumed so Luz and ey are 0 The remaining cross coupling terms are computed using measured ata the assump tions inherent in th
68. is usage are discussed in reference 26 Whenever measured a and B are used they are corrected for upwash and center of gravity offset 43 4 1 2 1 Measured data are also used to linearize the other nonlinear terms in the equations There is some leeway for engineering judgment in the use of the measured data Strictly speaking the linearized longitudinal equations should be obtained by a first order expansion about the measured longitudinal data similarly for lateral directional equations However many of the nonlinear terms are satisfactory when simply evaluated at the measured data that is using a zero order expansion A value judgment must be made for each nonlinear term to decide whether the additional fidelity of the first order expansion justifies the additional complexity The implementation in the standard aircraft routines represents the authors judg ment of these tradeoffs modification for other choices is relatively simple 4 1 2 1 Longitudinal The longitudinal state control observation and extra signal vectors respectively are x aq 9 u 8 Pe oj 59 3 m Im on an a m 94 extra F B p r MHhVprN T The nonlinear longitudinal state equations are Y 48 gt g i i B i T j VRC dy q Z cos cos cos a sin a sin 9 p cos a r sin a z sin a y p2 Bz ONY 55 1 b n hi 1 4 SCRC n 2 Er R lie I zZz x q cos r sin
69. its true value The second term can therefore be neglected if the algorithm starts near enough to the converged solution The Newton Raphson algo rithm with this term neglected is referred to as the Newton Balakrishnan algorithm 1 2 2 Improvement of Convergence Properties Since Newton Raphson type methods converge only if the starting estimates are close enough to the minimizing values it is common to use special startup algorithms to improve the probability of convergence For this purpose MMLE3 contains an option sec 3 3 8 19 to estimate only control derivatives biases and initial con ditions for the first iteration The cost functional J amp is quadratic in these unknowns when the other unknowns are constrained to the starting values thus the Newton Balakrishnan algorithm attains the constrained minimum in one iteration In subse quent iterations all unknowns are estimated A program option is available to aid convergence by multiplying the diagonal elements of the second gradient matrix by a constant With this option invoked the algorithm becomes similar to the gradient algorithm ref 45 convergence becomes more nearly monotone but requires many more iterations than with Newton Raphson The gradient like algorithm tends to stall before reaching the minimum Since this stall is easily confused with convergence to the minimum the final iterations should always be done with the full Newton Balakrishnan algorithm Use of the
70. ively uncorrelated with K and thus with G F however enters the cost function only through K therefore unknowns in F are strongly affected by changes in G A heuristic adjustment for F is derived assuming that the optimum K is independent of G therefore when G is changed F is adjusted to minimize the resulting change in K This is effected to a first approximation by multiplying the ith row of F by og J gd gt On g 7 0 E ADA 5 27 C g d j If F was at the optimum for the previous G this adjustment puts F near enough to the optimum for the new G that convergence to the new optimum is rapid Only elements of F that are varying are adjusted using this algorithm 1 6 Solution of the Riccati Equation This section discusses the methods used for the computation of the Riccati covariance matrix P and its gradients MMLE3 uses only the steady state solution of the Riccati equation The discrete time Riccati equation is usually written as a function of the measurement noise covariance matrix G9 MMLE3 however works in terms of the residual covariance matrix GG The discrete time Riccati equation in terms of GG is DPO P PC GG TCPO At OFF O 0 28 Equation 28 is closely approximated by the continuous time Riccati equation 1 E PC GG CP FF 0 29 AP PAS MMLE3 uses equation 29 because its solution is more convenient than that of equation 28 18 1 6 Potter s metnod
71. l option to print final output time history If PRINTO is TRUE the time histories of the corrected observation estimate are printed in the final iteration PRINTO is irrelevant if PRINTY item 31 is TRUE The default value of PRINTO is FALSE 31 PRINTX PRINTY logical options to print states and observations for all iterations If PRINTX or PRINTY is TRUE time histories of the predicted state esti mates or corrected observation estimates respectively will be printed each iteration If both PRINTX and PRINTY are TRUE the predicted state at each time point will be printed on one line and the corrected observation on the next line The printout columns of the states and observations are staggered so they can be easily distin guished The default values of PRINTX and PRINTY are FALSE 32 ERRTH logical option for time history printout after errors If ERRTH is TRUE and a case terminates as a result of exceeding the ERRMAX item 18 or 37 3 3 8 PLTMAX item 34 bounds the measured observations controls and extra signals will be printed This printout provides a check to see if any obvious data problems are present If PRINTI item 29 is TRUE ERRTH is forced to FALSE since the printout would be redundant in this case The default value of ERRTH is FALSE 33 PLOTEM logical time history plot option If PLOTEM is TRUE time history plots comparing measured and predicted time histories are produced unless the ER
72. les are also listed The default values can be altered by the user subroutine USERIN unless user routines are bypassed sec 3 3 2 The changes in default values made by the standard aircraft routines are listed in section 4 3 All values are reset to the defaults at the beginning of each case regardless of any values used in previous cases 1 NCASE number of maneuvers to be analyzed as a single case If two or more maneuvers were performed at disjoint times they may be analyzed together to obtain a single set of estimates Program dimensions limit NCASE to 15 or less One time card sec 3 3 10 is required for each maneuver The use of multiple maneuvers in this manner presumes that the values of the unknown coefficients are the same for all of the maneuvers in the case The default value of NCASE is 1 2 CARD TAPE logical input source for the time history data These var iables are passed to the user subroutine READTH to allow selection of alternate sources for the input time history data Either CARD or TAPE may be set to TRUE A third selection is possible by specifying both CARD and TAPE to be FALSE Both cannot be TRUE if CARD is TRUE TAPE will be forced to FALSE READTH can assign any data sources to the three available combinations the sources need not correspond to the variable names The READTH subroutine supplied reads from cards if CARD is TRUE or from tape otherwise see secs 4 2 1 and 4 3 4 for the data formats It
73. lts are Longitudinal Lateral directional KALF KALF X XALF DCGFT 0 KB KBXZB KBX XB DCGFT 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 Cy CN 0 0 0 0 1 a S Cy Cy 0 A A p l q 0 0 0 0 0 0 0 0 0 0 where Es ES and a are obtained from AN Y AN 2 and AN 3 KALF and KB are described in section 4 3 3 11 XALF and XB in 4 3 3 12 and ZB in 4 3 3 14 DCGFT is discussed in section 4 1 3 The defaults for DN are Longitudinal Lateral directional 0 0 0 0 0 0 0 C Ne 0 0 e Cy Cy 5 5 C A a r 5 e 0 The longitudinal values in DN are taken directly from the predicted derivative file The lateral directional values are obtained from C BN and C BN Y Fl Ys 1 2 n r 71 4 3 5 The defaults for HN are Longitudinal Lateral directional 0 0 0 C N 0 0 Cy C A 0 0 Although the total trimmed C N and C A should be placed in HN from the predicted derivative file the equations of motion require C N and C A in these locations Subroutine WTDEF computes 0 0 C SOC as No N Ny Ns 67 Ca C TC aC 6 0 a 5 and replaces C andC by C andC in the HN matrix N A No Ay The defaults for EN are Longitudinal _Lateral directional 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 XAN DCGFT 0 0 0 0 0 0 ZAX 0 O ZAY XAY DCGFT 0 0 1 0 0 1 0 0 0 0 1 0 XAN and XAY are described in section 4 3 3 12 and ZAX and ZAY in 4 3 3 14 DCGFT is discussed in section 4 1 3 2 F s
74. lumns like that shown for each element of UVAR that is 1 For elements of UVAR that are 0 the corresponding column is 0 74 4 3 5 The defaults for HV are Longitudinal Lateral directional Rh e CC oo e e O C e e CO oo o h e e oOo o o e e e O o o o Rh e e C o C e COC o o o e O o o o If the weighting diagonal GGI element for any signal is 0 the elements in the corresponding row of HV are forced to 0 this is true whether HV was read in or defaulted The defaults for EV and FV are 0 4 APRA APRB APRS APRR APRC APRD APRH APRE APRF a priori weighting matrices The defaults for APRA are Longitudinal Lateral directional 100 0 0 300 0 0 0 m 1 0 1000 10 10 0 0 0 0 1000 10 10 0 0 0 0 0 The defaults for APRB are Longitudinal Lateral directional 300 300 300 300 1000 1000 0 1000 1000 0 0 15 4 3 5 The default for APRC is 0 in a lateral directional case in a longitudinal case the default is 0 0 0 0 0 0 0 0 0 100 0 0 100 0 0 0 0 0 The default for APRD in a lateral directional case is 0 in a longitudinal case the default is 300 300 The defaults for APRS APRR APRH APRE and APRF are 0 5 GGI inverse of GG The GGI default is specified in the user initialization input sec 4 3 1 If not otherwise read in the default is a diagonal matrix of size 5 by 5 The diagonal values are Longitudinal 10 60 30 200 200 Lateral direc
75. ly tied to the treatment of the initial condition The basic program uses pertur bation equations The state initial condition is defined as 0 UOFF and YOFF are set to the measured control and observation at the first time point Other treatments of the initial conditions and biases are possible through user routine UINIT If infor mation is available to define the initial state in absolute as opposed to perturbation terms it would be natural to define UOFF and YOFF as 0 The vector 7 t is used to allow unknown bias terms to be included in the model The first element of 7 t is defined as a constant value of 1 for the entire case in order to allow constant bias terms The remaining elements of 7 t are used to allow the biases for multiple maneuvers sec 3 3 8 1 to be independent The second element of 7 t is defined to be 0 during the first maneuver and 1 for all subsequent maneuvers the third element is defined to be 0 during the first and second maneuvers and 1 for all subsequent maneuvers and soon Consequently the biases for the first maneuver are the derivatives for the first element of 7 t the biases for the second maneuver are the sums of the derivatives for the first and second elements of a t and so on The first MAXB maneuvers can therefore have independent biases the biases for any subsequent maneuvers are equal to the biases of the MAXBth maneuver The MMLE3 program allows the system matrices to be time varying this is
76. matrices The default values of these matrices are all 0 except for RN which defaults to an identity matrix The default matrix sizes are all set to 0 The sizes of the AN BN and SN input matrices are sometimes used to determine MX MU and MB sec 3 3 8 11 13 and 14 All matrix sizes are adjusted by the program to be consistent with MX MU MB and MZ sec 3 3 8 12 Maximum matrix sizes are discussed in section 2 2 F state noise power spectral density matrix The default value for F is 0 3 AV BV SV RV CV DV HV EV FV variation matrices These matrices indicate which elements in the nondimensional matrices item 1 and F item 2 are allowed to vary A value of 1 indicates that the corresponding nondimensional or F matrix element is allowed to vary independently a value of 0 indicates that it is held fixed or is constrained to another element see HARD item 6 The defaults for the variation matrices are all 0 40 3 3 11 4 APRA APRB APRS APRR APRC APRD APRH APRE APRF a priori weighting matrices These matrices contain relative weighting data for the a priori feature The value of each element in the a priori weighting matrices should be pro portional to the inverse of the a priori standard deviation of the corresponding non dimensional coefficient A priori weightings are only implemented for independent unknowns that is matrix locations corresponding to values of 1 in the variation
77. matrices item 3 A priori weightings on dependent unknowns in hard constraints item 6 are ignored The overall a priori weighting is adjusted by WAPR sec 3 3 8 22 The actual a priori weighting for each coefficient is WAPR times the square of the corresponding APR element The default values for the a priori weighting matrices are all 0 5 GGI inverse of the residual covariance matrix The GGI matrix is often diagonal If the GGI used is diagonal the program will take advantage of the diagonal condition to reduce computation time If GGI is not symmetric an error message will be printed and the upper triangular part will be replaced by the transpose of the lower triangular part to force symmetry The size of the GGI input matrix is sometimes used to determine MZ sec 3 3 8 12 Maximum sizes are discussed in section 2 The default value for GGI is the 0 matrix 6 HARD hard constraints Constraint matrices are read in the special form described here The header card contains the matrix name HARD as usual and the number of constraints replaces the number of rows A special provision is made to control whether constraints from cards replace or supplement any default con straints If any nonzero value is punched for the number of columns on the header card of the HARD matrix the default constraints will be used in addition to the ones read in If the number of columns is left blank on the header card of the HARD matrix inp
78. med over all controls The equation for C is linearized E about the average measured a If MZ sec 3 3 8 12 is 4 the approximation Cr Cy is used This approximation is good for low angles of attack 45 4 1 2 2 4 1 2 2 Lateral directional The lateral directional state control observation and extra signal vectors respectively are x 6B pr q w 8 Sy 53 58 A A Bn Pm m Pm Yn Pm L n da 1 e extra Taqa MhV a NT The nonlinear lateral directional state equations are E qs g 1 i B Ay KR Cy Bp y Rcos 8 sin p sina rcosa I qr 7 XZ pl e SDRC o R 1 1 pq oe c _ _ lez 6N PI BL TSDC R RL I wr at Lae p p r cos tan 9 q sin tan 9 gt The Bo and Po are included to allow for instrument biases The p and equations already have the freedom to correct for such biases using Co and an the Cy in 0 0 0 the equation does not allow similar freedom because it also appears in the ay equation below Measured data are used for a q and 0 in these equations Measured is used for the trigonometric terms in the equation to eliminate the nonlinearity The gravity term in the B equation is normally linearized about the measured q however it can be evaluated at the measured if it is not desired to integrate the equation 46 4 1 2 2 The observation equations are Z X Bm Kap p Pr Pm 7P ro r m
79. mpossible constraints are not imposed MMLE3 constrains only those diagonal elements of KC which correspond to unknown diagonal elements of F 1 2 Minimization of the Cost Functional The maximum likelihood method consists of choosing the vector of unknowns to maximize the likelihood functional The likelihood functional has been defined in the previous sections It remains to describe an algorithm for the maximization The problem is usually restated as minimizing the negative log likelihood functional of equation 7 or 10 referred to as the cost functional 1 2 1 Newton Balakrishnan Algorithm The Newton Raphson algorithm is a well known technique for iterative mini mization of nonlinear functions An initial estimate 0 is required The estimate is then revised iteratively using the equation e a iZ a Sai Si vz7 5 V J ate The properties of this algorithm are well documented in the basic textbooks ref 45 13 a v The Newton Raphson algorithm requires values of the first and second gradients of J Computation of the first gradient is straightforward N VIE 2 40 2 t GG vez 19 The second gradient is N 2 P r A _ E gest 22 r US ALEA GG Vel Ze 1 Z t 2 t GG AKG 20 i i Computation of the first term of equation 20 is straightforward The second term requires inordinate amounts of computer time and core however it has expected value 0 if is at
80. nal case respectively The default case type is lateral directional 2 SHIFT logical option to shift center of gravity reference from wind tunnel to flight If SHIFT is TRUE all predicted moment derivatives are corrected by the program for the longitudinal distance between the wind tunnel reference center of gravity sec 4 3 2 2 5 and the flight center of gravity location CG item 10 The equations used for the moment corrections are C n Cr CG q CGref CN cc Sca o fit ref 65 7 CHORD Ch Cn CGn Cer SPAN CY ScG ng fit ref where 53 represents any a B control or bias derivative In addition all longi tudinal instrument locations item 12 are assumed to be given as offsets from the wind tunnel reference center of gravity the program adds the center of gravity difference to these offsets to obtain offsets from the flight center of gravity location If SHIFT is FALSE all predicted data are assumed to be referenced to the actual flight center of gravity location No corrections are made to the moment derivatives or instrument locations regardless of any values specified for the wind tunnel ref erence or flight centers of gravity No corrections for vertical or lateral center of gravity shifts are included in the program whether SHIFT is TRUE or not If the predicted derivative control card sec 3 3 3 specifies NO SHIFT is forced to FALSE since no wind tunnel reference is defined The default value o
81. nce The state noise algorithm is similar in form to the algorithm used without state noise indeed it reduces to the same algorithm if the state noise power is 0 thus implementation in the same computer program was relatively easy MMLE3 has the option to use either the state noise or the no state noise algorithm Miff s original implementation of the state noise algorithm used the pure contin uous time formulation discussed in references 36 and 37 Recent unpublished studies by the authors conclude that a continuous system model with discrete sampled observations is superior in the state noise case the formulations result in identical algorithms when there is no state noise MMLE3 therefore uses a mixed continuous discrete time formulation Section 1 of this report outlines the theory behind the MMLE3 program Sec tion 2 defines the overall structure and coding philosophy of the program the most important item discussed is the division of MMLE3 into the basic program and the user routines Section 3 documents the use of the basic MMLE3 program which is a general maximum likelihood estimator for use with any linear system A set of user routines designed for the standard aircraft stability and control problem is described in section 4 Program listings and reference maps along with detailed descriptions of the purpose and operation of each subroutine are found in the Programmer s Manual ref 38 Particular attention is paid to the desc
82. nd Evaluation of Two Methods of Extracting Stability Derivatives From Flight Test Data AFFTC TD 73 5 Air Force Flight Test Center Edwards Air Force Base May 1974 Park Gary D Determination of Tail Off Aircraft Parameters Using Systems Identification Proceedings of AIAA 3rd Atmospheric Flight Mechanics Con ference c 1976 pp 128 136 Balakrishnan A V Stochastic Differential Systems I Filtering and Control A Function Space Approach Lecture Notes in Economics and Mathematical Systems 84 M Beckmann G Goos and H P Kunzi eds Springer Verlag Berlin 1973 lliff K W Identification and Stochastic Control With Application to Flight Control in Turbulence UCLA ENG 7340 School of Engineering and Applied Science Univ Calif Los Angeles Calif May 1973 Maine Richard E Programmer s Manual for MMLE3 A General FORTRAN Program for Maximum Likelihood Parameter Estimation NASA TP 1690 1981 Cramer Harald Mathematical Methods of Statistics Princeton Univ Press 1946 Klein V Parameter Identification Applied to Aircraft Cranfield Report Aero No 26 Cranfield Inst Technol Oct 1973 41 42 43 44 49 46 47 48 49 90 91 Schweppe Fred C Uncertain Dynamic Systems Prentice Hall Ine 1977 Taylor Lawrence W Jr Application of a New Criterion for Modeling Systems Methods for Aircraft State and Parameter Identification AGARD CP 172 May 1975 pp 4 1 to 4
83. ne will be read as the title card for the next case If too few cards are present a read error or end of file error will occur 4 3 5 Matrix Defaults As previously mentioned all of the input matrix defaults can be changed by the user routines This section describes the input matrix defaults for the standard aircraft routines Although in some cases these defaults are unchanged from the basic program defaults described in section 3 3 11 they are repeated here for completeness descriptions of input format and use of the matrices are not repeated The dimension alization matrices used to define the equations of motion are not described in this section as they are defined internally in routines MAKEM and MAKEL rather than being read in The descriptions of these matrices for the standard aircraft routines are found in section 4 1 3 In general the matrix sizes are determined by the variables MX MZ MU and MB sec 3 3 8 11 to 14 This section describes the default values of the matrix elements but does not discuss how large a partition of each matrix will actually be used Any locations not shown herein are assumed to default to 0 Many of the matrix defaults are affected by LONG and LATR sec 4 3 3 1 A longitudinal default will refer to a default used if LONG is TRUE and conversely the lateral directional defaults are used if LATR is TRUE where neither type is specified the default applies to all cases 67 4 3 9 1 AN
84. new case Any number of cases can be run in one program execution The last case is signified by END on the endcase card The input flow chart on page 29 shows all of the input card groups and the order of their appearance Each card group is denoted as basic program or user routine Card groups which may be omitted under some circumstances are labeled as op tional Of course all of the user routine card groups are omitted if user routines are not used furthermore a specific set of user routines may require no cards in a particular group A user routine card group is not marked as optional unless the basic program has an option to omit the card group while user routines are being used Details are given in the appropriate subsections of this section 3 3 1 Syntax Check Card MMLE3 has an option to check the syntax of the input cards without executing the estimation This option can be useful on runs with many cases which will require large amounts of computer time In order to activate this option a card with SYNTAX in the first seven columns is added to the front of the deck After all the syntax errors are found the syntax check card should be removed and the job re submitted for actual execution When the syntax check option is used the input time history file is not used Therefore the input time history file sec 3 2 2 need not be available and the time history card input sec 3 3 13 should be omitted If syntax checki
85. ng is not desired the input deck begins with the user routines control card 3 3 2 User Routines Control Card This card controls the calls to most of the user routines Only the first three characters of the card are tested but the entire card is printed on the output listing so it may be used for comments on the listing The card must not have SYNTAX in the first seven columns or it will be mistaken for a syntax check card sec 3 3 1 If the first three characters are NO the user routines defining the equations of motion are bypassed and the card groups described in sections 4 3 3 4 3 4 4 3 5 and 4 3 7 should be omitted In this case the equations are defined solely by the matrices in the card input The default values of nondimensional matrices N matrices 28 3 3 2 Input flow chart Syntax check card Basic Optional User routines control card Basic Once per User initialization input run User Predicted derivative control card Basic Optional Predicted derivative input User Optional Title card Basic User input User NAMELIST INPUT Basic Signal labels Repeated Basic Optional for each case Time cards Basic Matrix input Basic Optional Endcase card Basic Time history cards User Optional 29 3 3 3 are all 0 except for RN which defaults to the identity the dimensionalization ratio matrices M matrices are all 1 s and the dimensionalization addition ma
86. nitial state the initial control and an ob servation bias to be added to the computed observations The program will compute a measured control bias by subtracting the measured initial control from that specified by UINIT This control bias will be added to the measured controls for each time point USERIC is sometimes necessary when TIMVAR item 20 is used If USERIC is FALSE the initial state and control are set to 0 and the observation bias is set to the initial measured observation vector The program prints out the initial state and control vectors used for each maneuver The control time histories printed and plotted are the measured controls without any biases added The default value for USERIC is FALSE This default is changed by the standard aircraft routines 36 3 3 8 27 VARIC MAXX word logical vector variable initial condition increments for states For each element of VARIC that is TRUE an increment to the initial condition of the corresponding state is estimated as one of the unknowns This increment is added to any initial condition defined by user routine UINIT to obtain the initial condition used to compute the time histories If USERIC item 26 is FALSE the increments are added to the value 0 since UINIT is not called If VARIC is used in conjunction with NCASE gt 1 item 1 the same increment from the values defined by UINIT or 0 is used for each maneuver Caution must be taken with the use of VARIC becaus
87. nnot be used on the predicted derivative file subroutine WTIN expands and reorders the data before writing them to the file 4 2 2 3 End card The end of the predicted derivative file is indicated by a card with the characters END starting in column 1 This card is necessary be cause FORTRAN end of file checks are machine specific 4 2 3 Punch File Subroutine OUTPUN is called to punch the estimates if PUNCH sec 3 3 8 45 is TRUE The first card punched is the title card sec 3 3 6 The second card contains the case type Mach number sec 4 3 3 20 angle of attack sec 4 3 3 17 extra parameter sec 4 3 3 21 and center of gravity sec 4 3 3 10 in format A4 6X 4G10 3 The case type is either LONG or LATR sec 4 3 3 1 If the date and time are available to the program see the Programmer s Manual ref 38 they will be punched in the last two fields of 10 columns on the second card The third card contains the value MAXZ sec 2 followed by MAXZ channel averages for the measured observations This card is in format Z AVG 13 7G10 3 If MAXZ is greater than 7 continuation cards in format 8G10 3 are used The fourth card contains U AVG the value MAXU and the MAXU channel averages for the controls The fifth card contains EX AVG the value LEX and the LEX channel averages for the extra signals The formats of the fourth and fifth cards are the same as for the third card including the possibility of
88. oints in 8G10 3 format The break points are continued to further cards in the same format if required 4 2 2 2 Predicted derivative data This section describes the format of the predicted data for one derivative The format is repeated as many times as there are derivatives No particular order of the derivatives is required For each derivative there is a derivative header card followed by a derivative data table The header card contains the derivative name type and location in format A4 6X A4 6X A2 1X I2 1X I2 The name is not used by the program but is solely for ease of user identification The type should be either LONG or LATR to indicate a longitudinal or lateral directional derivative if any thing else is specified for type the derivative will effectively be ignored The 37 4 2 2 3 derivative location is specified by matrix name AN BN SN RN CN DN HN or EN row and column The derivative data table contains the derivative values at the angle of attack break points on one card in 8G10 3 format further cards in the same format are added if there are more than eight angle of attack break points This card or cards is repeated for each Mach number break point Finally all of the above cards are repeated for each parameter break point This data table corresponds to that described in section 4 3 2 for AMP functional dependence The other compressed and reordered forms described in section 4 3 2 ca
89. only case The MMLE3 program accounts for both state and measurement noise The measurement noise only algorithm can also be used 8 1 1 The final problem to be faced for real systems is modeling It has been assumed throughout the above discussion that for some value called the true value of the unknown parameter vector the system is correctly described by the dynamic model Physical systems are seldom described exactly by simple dynamic models so the question of modeling error arises There is no comprehensive theory of modeling error available The most common approach taken amounts to ignoring it Any modeling error is simply treated as state noise and or measurement noise in spite of the fact that the modeling error may be deterministic rather than random The assumed noise statistics can then be adjusted to include the contribution of the modeling error This procedure has not been rigorously justified but combined with careful choice of the model is probably the best approach available Schweppe ref 41 discusses this question in heuristic terms Several methods have been advocated which purport to determine the structure of the dynamic model from the data refs 42 and 43 This basic concept however is fraught with problems for real systems It is likely that model structure determin ation will continue to require the engineer s careful consideration of the phenom enology of the physical system to supplement the automati
90. ossible distinctions in clude aircraft configuration altitude and power setting NPBP cannot be greater than 20 or less than 1 The default value is 1 5 CG reference center of gravity in fraction of chord The default value is 0 25 6 PRINT logical option to print predicted data If PRINT is TRUE all of the predicted data will be printed otherwise only the title card reference center of gravity axis system break points and header cards will be printed The de fault value of PRINT is FALSE 7 AREA SPAN CHORD reference area span and chord respectively at and ft or m and m The default values are all 0 8 KALF KB real flow amplification factors for angle of attack and angle of sideslip The default values are both 1 9 XALF XB XAN XAX XAY distances of the angle of attack angle of sideslip normal acceleration longitudinal acceleration and lateral acceleration sensors respectively forward of the center of gravity ft orm The center of gravity location referenced is the wind tunnel reference center of gravity if SHIFT sec 4 3 3 2 is TRUE or the actual flight center of gravity if SHIFT is FALSE The default values are 0 60 4 3 2 3 10 YALF YB YAN YAX YAY distances of the angle of attack angle of sideslip normal acceleration longitudinal acceleration and lateral acceleration sensors respectively right of the center of gravity ftorm The default values are 0
91. porary scratch file 1 7 UT2 Temporary scratch file 2 8 3 2 1 Card Reader and Line Printer Files The card reader and line printer files are assigned FORTRAN device numbers UCARD and UPRINT respectively A detailed description of the data for the card reader file is given in section 3 3 The line printer output data are strongly depend ent on the input data 3 2 2 Input Time History File The FORTRAN file number for the input time history file is UDATA The format of this file can be defined by the user READTH is the user routine to read the input time history file no reference to the file is made in any other routine READTH can be written to read the time history data from the card reader sec 3 3 13 or to compute simulated time history data in these cases file UDATA is not used Various aspects of the data from the input time history file are described in section 3 3 8 items 2 to 9 3 2 3 Predicted Derivative File The FORTRAN device number for the predicted derivative file is UWT This is a miscellaneous file which can be used for any special purpose required by the user This file is not restricted to predicted derivative tables it can be used to store any data required by the user routines The file can be used for input output or 26 3 2 4 scratch storage Since the purpose of this file can be defined by the user it can be referenced in any of the user routines No reference to this file is made outside of the use
92. ption of the standard aircraft user routines in section 4 should be regarded as a supplement to the description of the basic program in section 3 The MMLE3 program is coded to facilitate changing the maximum matrix dimen sions Systems with matrix sizes less than the maximum can be analyzed without changing the code however it may prove useful to reduce the maximum dimensions in order to lower core requirements If the required matrix sizes are larger than the current program dimensions the code must be changed Details of an automated system for easily changing the maximum matrix dimensions of the program are des cribed in the Programmer s Manual ref 38 The names used in the program and in this report to refer to the maximum matrix dimensions are listed below along with the current values of the maximum dimensions Variable Name Description Current Value LEX Maximum length of extra signal vector 20 LORD MAXZ MAXU LEX 32 MAXB Maximum length of bias vector 4 MAXHRD Maximum number of hard constraints 1 36 MAXSFT Maximum number of soft constraints 1 11 MAXFV Maximum number of independent unknowns 4 in F MAXKV Maximum number of independent unknowns 15 affecting K MAXTV Maximum total number of unknowns 90 including constraints MAXU Maximum length of control vector 4 MAXX Maximum length of state vector 1 MAXZ Maximum length of observation vector 8 NI Maximum number of independent unknowns 2 39 The maximum number of time point
93. quation is identical to equation 6 except that 7 is now com puted using the Kalman filter ref 4 and G replaces Y Naturally equation 6 can be viewed as a special case of equation 9 with 0 for the Kalman filter gain 11 1 1 2 There are two subtle aspects of the differences between equations 6 and 9 First in equation 9 Z is a function of the measured data z because of the Kalman filter For the most part this is of no concern but it creates some problems in the exact computation of the Cram r Rao bound sec 1 4 Second G is a complicated function of the coefficients even if 8 is known MMLE3 avoids this problem by determining G directly instead of GS For known G the cost functional to be mini mized is N JE T Z 1 2897 Gc 1 KOKO 10 i 1 which corresponds to equation 7 The equation for Ze is as follows Prediction step Fe tina 0 g 1 VB y Z a 11 tina Cty tuy tPU tu where ge AA 12 At Y f a ds J0 Correction step e tina y ai E 2 tier T Ze tir 13 The Kalman filter gain matrix is defined as K PC GG 14 where P is the solution to the discrete time Riccati equation MMLE3 uses the steady state continuous time Riccati equation for an approximate solution for P 1 AP PA Pc GG CP FF 0 15 12 1 2 Although this algorithm incorporates a Kalman filter it should not be confused with the extended Kalman fil
94. r characteristics to the operating system File UT1 is used during iteration for the measured time history data for a case Scale factors and biases sec 3 3 8 8 and 9 and any corrections made by subroutine THMOD have been applied to the data on this file Each record of file UT1 contains LORD 4 FORTRAN variables File UT2 is used for measured and estimated time history data for plotting Each record of file UT2 contains 2 X MAXZ MAXX MAXU LEX FORTRAN variables 3 3 Input Description This section describes the card input required to run the basic program Each subsection describes a group of input cards The subsections are presented in the order in which the card groups appear in the input deck 27 3 3 1 Subsections are included for input to the user routines The input cards re quired for these subsections depend on the user routines therefore the details of such input are not documented in this section For the basic program such cards are omitted If user routines are called sec 3 3 2 any input required for the user routines should be placed in the locations indicated by these subsections The special input requirements for the standard aircraft routines are described in sec tion 4 3 The syntax check card user routines control card user initialization input predicted derivative control card and predicted derivative input occur only once per program execution The remaining input card groups are repeated for each
95. r routines 3 2 4 Punch File The FORTRAN device number for the punch file is UPUNCH The format of this file can be defined by the user Subroutine OUTPUN is the user routine to write the punch file no reference to this file is made in any other routine The information written on the punch file would normally include the final derivative estimates and Cram r Rao bounds 3 2 5 Output Time History File The FORTRAN device number for the output time history file is UTHOUT The format of this file can be defined by the user Subroutine THOUT is the user routine to write the output time history file no reference to this file is made in any other routine Both the measured and last iteration predicted response time histories are available to subroutine THOUT as well as the last iteration corrected state estimates and the control and extra signal time histories Any of these signals or data com puted from them can be written to the output time history file 3 2 6 Plot File The FORTRAN device number for the plot file is UPLOT The software on some systems may ignore the device number and create a file with a special name defined by the system System dependent details of plotting are discussed in the Programmer s Manual ref 38 3 2 7 Internal Scratch Files FORTRAN file numbers UT1 and UT2 are used for unformatted internal scratch files The user generally does not need to be concerned with these files except possibly to describe thei
96. rd The first card of each case is an 80 column title card This card is used on page headings plotted output and often in the punched output created by user routine OUTPUN 3 3 7 User Input Any input required by the user routine USERIN should be inserted at this point in the input deck The input required by the standard aircraft routine USERIN is described in section 4 3 3 If user routines are bypassed sec 3 3 2 USERIN will not be called and the user input should be omitted 3 3 8 NAMELIST INPUT The input for this section is in the form of a FORTRAN NAMELIST called INPUT The FORTRAN manuals for specific computer systems should be consulted for the format required by these systems The variables in the NAMELIST are listed below Most of these variables are scalar with type real or integer depending on the standard FORTRAN type convention A first letter of I J K L M or N indicates an integer variable any other first letter indicates a real variable Vector variables and excep tions to the standard type convention are indicated parenthetically after the variable 31 3 3 8 names Vector lengths are defined by the variables described in section 2 Some closely related variables are listed together under a single item number The variables are grouped into those relating to the input items 1 to 10 the esti mation process items 11 to 27 and the output items 28 to 45 The default values of all of the variab
97. ription of the user routines so that the engineer can understand how to program or modify a set of user routines for his specific problem SYMBOLS All data are referenced to fuselage body axes according to right handed sign conventions A state equation matrix A dummy matrix a dummy variable a normal acceleration g a longitudinal acceleration g a lateral acceleration g B state equation matrix b reference span m ft b dummy variable C observation matrix C A axial force coefficient Cp drag coefficient CG center of gravity fraction of chord C L lift coefficient C Q rolling moment coefficient C En pitching moment coefficient normal force coefficient yawing moment coefficient lateral force coefficient reference chord m ft observation matrix observation matrix expected value state noise power spectral density matrix arbitrary function residual covariance matrix measurement noise covariance matrix acceleration of gravity m sec ft sec diagonal element of 99 gt arbitrary function observation matrix altitude m ft identity matrix moment of inertia about roll axis N m slug ft 2 2 engine moment of inertia N m slug ft l i 2 2 cross products of inertia N m slug ft moment of inertia about pitch axis N m slug ft moment of inertia about yaw axis N m slug ft zZz 3 B KE 3 p Z cost functional Kalman f
98. rom the actual flight center of gravity positive for sensors right of the center of gravity The units used for these offsets depend on METRIC item 3 The defaults are obtained from the predicted derivative file secs 4 2 2 and 4 3 2 2 10 if present otherwise the defaults are 0 65 4 3 3 14 ZALF ZB ZAN ZAX ZAY vertical locations of the angle of attack angle of sideslip normal acceleration longitudinal acceleration and lateral acceleration sensors respectively ft orm These locations are all given as offsets from the actual flight center of gravity positive for sensors below the center of gravity The units used for these offsets depend on METRIC item 3 The default values are obtained from the predicted derivative file secs 4 2 2 and 4 3 2 2 11 if present otherwise the defaults are 0 15 QBAR average dynamic pressure Ib ft or Nim The units of QBAR depend upon METRIC item 3 If QBAR is set to 0 a value is obtained from the average of the measured dynamic pressure time history The default value of QBAR is 0 16 V average velocity ft sec or m sec The units of V depend on METRIC item 3 If V is set to 0 a value is obtained from the average of the measured velocity time history The default value of V is 0 17 ALPHA average angle of attack deg If ALPHA is set to 999 a value is obtained from the average of the corrected measured angle of attack time history The default v
99. s that can be plotted for each maneuver is 1200 plots will automatically be thinned to avoid exceeding this limit There is no soft ware limit on the number of time points that can be analyzed in a maneuver 3 0 BASIC PROGRAM This section describes the basic program and the input required to execute it In some cases the input variables of the basic program interact with the user 21 3 1 routines the most common interaction is for the user routines to change the default values of the basic program These interactions are described where appropriate In addition this section points out where the input cards for the user routines would be placed in the input deck 3 1 Equations of Motion Filter and Gradients The general bilinear equations of motion used in the MMLE3 program are R AXLA B u SID v Fn t 36 z t C t x t D u H EHL w Gn t Although equation 36 contains more terms than equation 8 in the theory discussion the differences do not affect the theory in any important manner Equation 36 can be recast in the form of equation 8 by appropriate substitutions The theoretical discussion of section 1 uses the simpler form to keep the details from obscuring the theory In this section the exact equations used by MMLE3 are given The vector u is the measured control minus a known bias UOFF and z is the measured response minus a Known bias YOFF The biases UOFF and YOFF are close
100. ses The BV and DV matrix defaults defined by UVAR are explicitly shown in section 4 3 5 2 For a longitudinal case UVAR defaults to1 0 0 0 6 66 4 3 4 derivatives will be determined For a lateral directional case UVAR defaults to 1 1 0 0 aileron and rudder derivatives will be determined 4 3 4 Time History Card Input This section describes the time history card input data sec 3 3 13 required by the standard aircraft routines The time history data are read from cards only if CARD sec 3 3 8 2 is TRUE otherwise time history data are read from tape or disk files and this card input should be omitted Other details of the time history data are described in sections 3 3 8 items 3 to 9 4 1 and 4 3 6 The format of the first data card for each time point is 312 13 1X 7F10 The time in hours minutes seconds and milliseconds is read in the integer format the last seven fields of the first card contain data for the first seven data channels Sub sequent cards contain the data for the remaining channels in 8F10 format The default order of the channels is a q V 9 a q a 6 8 5 n e ne one longitudinal 8 phMqppra pro 8 8 5 N T This longitudinal y lateral lateral order can be changed with ZCHAN UCHAN and EXCHAN sec 3 3 8 7 The last record of data on the cards should have the first time greater than or equal to the maneuver end time sec 3 3 10 If further time history data cards are present o
101. stimate of a diagonal element is ggo and the estimate from equation 24 is J a further revised estimate Jo is obtained from the equation Jo ag 25 where r 1 J1 a 5 75 26 gp 91 where b is the geometric average of the Ia values for all of the weighted signals and 0 r is a relaxation factor Values of r greater than 1 correspond to overrelaxation values smaller than 1 correspond to underrelaxation If r equals 1 there is no relaxation This relaxation is used only for diagonal elements of Geyt Each off diagonal element is multiplied by the square roots of the a values for the diagonal 17 1 6 elements in the corresponding row and column If the state noise algorithm is used the two step procedure described for esti mating G and the other coefficients results in slow sometimes erratic convergence The reasons for this problem are discussed and a heuristically based fix is described below The two step procedure does not compute correlations between G and the other unknowns It therefore converges best when such correlations are small When state noise is not present G affects the other coefficients only through the relative signal weighting in the cost functional In general the coefficients are quite insen sitive to changes in the weighting therefore the two step procedure converges very well When state noise is present G also affects the Kalman filter gain matrix K Most of the coefficients are relat
102. tate noise spectral density The default for F is specified by the user initialization input sec 4 3 1 If not otherwise read in the default for F is 0 12 4 3 9 3 AV BV SV RV CV DV HV EV FV variation matrices The defaults for AV are Longitudinal Lateral directional 0 0 0 1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 The defaults for BV depend on UVAR sec 4 3 3 22 For the default UVAR the BV defaults are Longitudinal Lateral directional 0 1 1 1 i 1 0 1 1 0 0 If UVAR is changed the BV defaults contain columns like those shown for each element of UVAR that is 1 For elements of UVAR that are 0 the corresponding BV column is 0 The defaults for SV are Longitudinal Lateral directional Ll lo th I 1 1 1 1 1 1 1 1 r dd xx 1 1 1 1 1 1 1 1 1 1 1 1 1 If the weighting on 6 GGI 3 3 is 0 the default of the third row of the longi tudinal SV is 0 Similarly if the weighting on GGI 4 4 is 0 the default of the fourth row of the lateral directional SV is 0 3 4 3 9 The defaults for RV are 0 The default for CV is 0 in a lateral directional case In a longitudinal case the CV default is e C o O OO o CC o O CC O The defaults for DV depend on UVAR sec 4 3 3 22 For the default UVAR the DV defaults are Longitudinal Lateral directional 0 0 0 1 1 0 0 0 0 0 0 0 0 If UVAR is changed the longitudinal DV default contains co
103. ter algorithm ref 44 The filter is not introduced arbi trarily into the formulation it appears in the implementation of the spectral factor ization of the noise The Kalman filter is used here only for the linear problem of state estimation The extended Kalman filter algorithm applies the Kalman filter in an ad hoc linearized manner to the nonlinear problem of combined state and param eter estimation The likelihood functional in MMLE3 is expressed in terms of the residual co variance matrix GG instead of the measurement noise covariance matrix GG This results in several simplifications in the algorithm It does however raise one complication which must be addressed F and are unrelated but there are some combinations of F and G that are not physically meaningful note that GG CPC 16 The left hand side must be positive semidefinite to be meaningful but the right hand side can be made negative for any GG by making F and thus P large enough The constraint GG CPC 2 0 17 must be satisfied to insure meaningful results This constraint is equivalent to re quiring that the eigenvalues of KC be less than or equal to 1 A necessary condition for this is that the diagonal elements of the matrix KC are less than 1 Because KC generally is nearly diagonal the diagonal constraint may reasonably approximate the eigenvalue constraint MMLE3 uses the diagonal constraint because of its sim plicity To insure that i
104. the overall accuracy If only one or two cases are available at a given condition however the scatter cannot be evaluated an alternate measure of the accuracy is necessary Even when many cases are used it is useful to have an indication of which individual estimates are the most reliable The Cram r Rao bound provides the best known analytical measure of the accuracy of the estimates 15 1 5 The Cram r Rao bound for unbiased estimators refs 4 5 and 39 is Variance gt GIK log p 2 8 Vg log p 2 0 JN 22 This equation gives only a lower bound for the variance of the estimates The maxi mum likelihood estimates however are asymptotically efficient this means that for large time intervals the variance is approximately equal to the expression in the above inequality provided that the system and noise are correctly modeled For the no state noise case equation 22 can be evaluated as N Variance y KON 66 v 2 1 23 i 1 When state noise is present the exact Cram r Rao bound is awkward to compute but Balakrishnan has shown that equation 23 approaches the bound for large time intervals ref 36 Equation 23 is seen to be in inverse of the portion of the second gradient matrix used in the Newton Balakrishnan algorithm eq 20 The program already computes this matrix so the Cram r Rao bounds are available with negligible extra computational effort The derivation of equation 23 is b
105. tion is not recommended unless convergence is not achieved with the standard algorithm The convergence check item 17 is not used on iterations with DFAC The minimum value of DFAC is 1 Very small increases in DFAC are usually sufficient to affect convergence significantly Typically used values range between 1 000001 and 1 1 The default values are 1 01 for DFAC and 0 for ITDFAC 35 3 3 8 22 WAPR ITAPR a priori control variables WAPR is the overall weighting factor for a priori information The a priori term in the cost functional sec 1 3 is proportional to WAPR A value of 0 for WAPR implies that the a priori feature is not used If WAPR is 0 ITAPR is also forced to 0 ITAPR is the maximum number of iterations for the a priori option If ITAPR is 0 the a priori weighting WAPR is used for all iterations If ITAPR is not 0 a priori is used only for ITAPR iteration or until convergence item 17 Iteration then continues without a priori The default value of WAPR is 0 and the default value of ITAPR is 0 23 FREQCR break frequency for first order residual filter Hz The power of the filtered residual for each observation is printed out each iteration The sum of the filtered residual powers weighted by GGI is used to adjust the Cram r Rao bounds sec 1 4 If FREQCR is 0 the residuals are not filtered and unfiltered residuals are used for the Cram r Rao bound adjustment The default value of FREQCR is 0 24 ITG
106. tional 150 5 300 10 5000 6 HARD hard constraints For a lateral directional case the default hard constraints are 76 CN C5 1 DN 5 1 oe AN 1 1 BN 1 1 4 3 5 1 i 1 2 3 1 i l For a longitudinal case the default hard constraints are BN 1 1 DN 4 i cosa 1 1 MU BN 1 i DN 5 i sina i 1 MU AN 1 1 CN 4 1 cos a 5 sin a AN 1 1 CN 5 1 sin a 2 cos a AN 1 1 HN 4 1 lt sin a AN 1 1 HN 5 1 cos a AN 1 1 DN 4 i 3 sin a 8 Tak y MU AN 1 1 DN 5 i cos a 8 i 1 MU _ P 1 4c AN 1 1 CN 4 2 a sin e s A qc AN 1 1 CN 5 2 E cos a 25 AN 1 2 CN 4 2 cosa AN 1 2 CN 5 2 sina MU is described in section 3 3 8 13 If MZ sec 3 3 8 12 is 4 a The only operative constraints are then the above constraints tudinal constraints BN 1 i AN 1 1 AN 1 2 H II DN 4 1 CN 4 1 CN 4 2 Average values of a q and are used in 0 is used in all of the longi 1 i 1l MU 1 1 Section 3 3 11 6 describes how the user can contro whether constraints read from cards supplement or replace the default constraints t 4 3 6 7 SOFT soft constraints The default includes no soft constraints 4 3 6 Altered Basic Program Defaults This section summarizes the changes made by the standard aircraft routines to the default values of the basic
107. titions the solutions to which can be written explicitly The inverse of the transformation of equation 33 is then used to obtain Y P S i 2 0 MMLE3 PROGRAM STRUCTURE Because of its generality MMLE3 sometimes requires a large amount of input to specify a given problem In addition some program features that are useful for specific problems have little or no meaning for other problems To satisfy the con flicting requirements of generality and ease of use the MMLE3 program is divided into two levels The basic level consists of a general maximum likelihood estimation program applicable to any linear system The basic program can be run by itself using input data to completely describe the system to be analyzed and the program options to be used The use of the basic program is documented in section 3 At the second level the basic program is used as the core of a program adapted to a particular application This adaptation is accomplished by a set of user routines so called because the user can write or modify them for a particular application Section 4 describes a particular set of user routines designed for the standard aircraft stability and control problem When the user routines are written for a particular application the program input does not have to contain the detailed system specifi cations only those items that change from case to case would need to be input This concept of a modular set of user routines allows t
108. trices L matrices are all 0 s thus the dimensional matrices are identical to the nondimen sional matrices The time varying option sec 3 3 8 20 is meaningless when the user routines are bypassed it will be overridden if attempted The subroutines bypassed by this option are ONCE WTIN USERIN AVERAG MATDEF MAKEL MAKEM and MAKEVW The remaining six user routines are not affected Subroutine READTH is always called when a case is executed since the program must obtain a set of time histories Subroutines THMOD THOUT and TITPLT are also always called though they perform no functions inherently necessary to the execution of the program The call to UINIT is controlled by the variable USERIC sec 3 3 8 26 and the call to OUTPUN is controlled by the variable PUNCH sec 3 3 8 45 The NO option on the user routines control card is useful for simple systems or occasional runs of complicated systems where it is more convenient to read as input all of the relevant matrices than to code routines to define defaults If any first three characters other than NO are punched on the user routines control card the user subroutines ONCE WTIN USERIN AVERAG MATDEF MAKEL MAKEM and MAKEVW are called The use of this option is suggested when many runs are made with the same or similar models The user routines would be coded to define most of the features of the models used input would only be required for those items in the models that
109. ts t The vectors n t are assumed to be zero mean Gaussian independent random vectors with identity covariance For arbitrary the model or calculated system is O EOL u t 2 2 ti g Xe t u t 5 The tilde notation for the calculated quantities is used for consistency with the state noise and measurement noise case The observation Z is measured at a finite number of time points Therefore it seems natural to maximize p Z the probability of given Z Unluckily p Z is difficult to define in general However we can take advantage of Bayes rule to write p S Z p Z p 2 p Z 8 p 8 3 and thus p E p Z PIOS 4 A specific Z has been measured so p Z is a constant independent of Further more we assume that there is no a priori preference for any particular value of so p is a constant section 1 3 discusses the case where a priori information is avail able Thus p Z differs from p Z only by a constant factor The maximum likelihood method works with p Z 8 instead of p Z because the former is easier to define Using the Gaussian assumption for the measurement noise we can write directly 1 N 5N oZy L ax l exp 2 gt 2 ali ggr 1 k ti i 5 i 1 10 1 1 2 The logarithm of the likelihood ratio is usually used because of its simpler form N log p Z 2g ti eG 1 z 1 a log 19 m log
110. ut the constraints read in will replace the default constraints The constraint instructions follow the header card one to a card in the format of the sample below CNG 1 AN 1 1 1 S ee Constrained Independent variable variable Valid matrix names here are AN BN SN RN CN DN HN and EN The names of the constrained variable and the independent variable start in columns 1 and 11 re spectively Either one or two digits can be used for the row and column numbers The constraint ratio is in columns 21 to 30 The constraint ratio in this example is 1 For a hard constraint the constrained variable should not be specified as varying in the variation matrices item 3 else the constraint will not work properly the variation matrices define independently varying elements only Conversely the independent variable must be specified as varying in the variation matrices else the constraint will be ignored as immaterial Thus hard constraints may not be chained A constrained to B constrained to C as B is not independently varying However an equivalent form can be obtained by constraining A and B to C 41 3 3 12 It is extremely important to note that only the change from the starting value is constrained the starting values may or may not satisfy the constraint It is solely the user s responsibility to enter the desired starting value As a consequence of this formulation it is valid to constrain one variable to
111. ve input is omitted This input is included if either NEW or ONL is specified on the predicted 99 4 3 2 1 derivative control card 4 3 2 1 Title card The first card of the predicted derivative input is an 80 column title card It is suggested that this card indicate the aircraft name data source initials of person preparing the data date prepared and any other iden tifying information necessary 4 3 2 2 NAMELIST WIND The input described in this section is in the form of a FORTRAN NAMELIST called WIND NAMELIST format is briefly discussed in section 3 3 8 For further details see the FORTRAN reference manuals for specific computer systems The parameters in NAMELIST WIND are described below 1 STAB logical stability axes If STAB is TRUE the longitudinal deriv atives are in stability axes otherwise they are in body axes The lateral directional data must be in body axes regardless of STAB The default value of STAB is FALSE 2 NABP number of angle of attack break points in the predicted data NABP cannot be greater than 20 or less than 1 The default value is 1 3 NMBP number of Mach break points in the predicted data NMBP cannot be greater than 20 or less than 1 The default value is 1 4 NPBP number of extra parameter param break points in the predicted data Param is used to represent any parameter other than angle of attack or Mach number that is used to distinguish derivative estimates P

User's Manual for MMLE3, A General FORTRAN Program

Contents

Download Pdf Manuals

Related Search

Related Contents