Cop~rig· 111
IF:\( : IdClllili(;lIillll Estilllatioll 1~IX :-., York . t ' K , I qK :-.
;tl l ci S',I(,[II 1':lr:IlI1t'!( ' !
FREQUENCY DOMAIN PARAMETER ESTIMATION OF AERONAUTICAL SYSTEMS WITHOUT AND WITH TIME DELAY M. Marchand and Fu K-h ii / 1,11//- 111It! Utl/IIII/flin/ (IJFI' IR J, III.I/i/II/ /Ii/ FlIlp,'/11I'(llIlIIik, ) )()() Ii UII 111,11 1111'I'ig, Fl'rll'ml UI'/iII;'/i{ of' (;I'I'IlItIlIY
[)('///vhl' FII/,1111/1II,!!.,I -
11111/ \ ' /'I ,III(h,IIIII ,I/II// /
Abstract. The paper considers the application of frequency domain methods to flight test data from rotorcraft and fixed wing aircraft. In particular, a modified frequency domain output error method is presented, which can be applied to linear systems with time delays. In contrast to existing methods, the new approach enables the combination of data from several manoeuvres for one evaluation, the use of test records with arbitrary initial state conditions x(t=O), and the direct identification of time delays (instead of using Pade-approximation). Results are presented as well from the identification of a six-degrees-of-freedom helicopter model (Bo-I05 Research Helicopter of the DFVLR) as from the identification of an equivalent low order system with time delay (HFB 320-In-Flight-Simulator of the DFVLR). Keywords. Identification; multivariable linear system; hood; frequency domain; aircraft.
I. INTRODUCTION System identification technique s are widely used for the evaluation of flight test data from fixed wing aircraft and rotorcraft. During the last ten years, the time domain methods have been shown to be particularly effective (Ref. I, 2). In most cases, a reduced mathematical model representing up to 6 degrees of freedom of the rigid body motion was used. For various applications, however, it is necessary to extend this model and to explicitly describe additional effects (for instance elastic or rotor degrees of freedom, dynamics of control system and actuators, time delays). In this case the evaluation in the time domain may reach its limits of applicability due to - high system order - large spread of smallest and largest eigenvalue - large number of data due to long data record, high sampling rate and large number of recorded input/output variables - large number of parame ters to be identified (aerodynamic parameters and additive constants within the equations). An approach to alleviate the numerical difficulties of time domain evaluat ion is the frequ e ncy domain evaluation. This approa c h enables a reduction of both the number of data to be evalua ted and the number of parame ters to be identified. Frequency domain identification methods have been used since about 195 0 (Ref. 3, 4). The approa ch i n the first methods was to minimize either th e t rans fer func t ion e rrors or th e equa t ion error s (Ref. 5 through 8). During the last few years, the theoretical backg round for the implementation of more advanced method s ha s bee n developed and output-error-methods as well a s maximum-likelihood methods have been proposed (Ref. 9, 10).
time-delay;
maximum likeli-
In this paper, a Modified Frequency Domain Method is presented and its extension to linear systems with time delay is developed. The paper is completed by examples dealing with helicopter and fixed wing aircraft identifications from flight test data.
2. THE IDENTIFICATION METHOD Basic approach. The basi c approach of the frequenc y domain identification method corresponds to the maximum likelihood method (V. Klein, Ref. 9). In this method, the error cost function J( 0 ) =
L £*("'k)
S-I £("'k) + log
I si
(I)
k
whi c h is derived from the log-likelihood-function, is minimized by a modified Newton-Raphson-iteration. In our institute, this approach was appli ed to a frequency domain model with the structure j<1.X(W) y( w)
=
Ax(w)
+
Bu(w)
= Cx( w) + Du( w)
(2a) (2b)
to identify the unknown elements of the system matric es A, B, C, D and of the error covariance matr i x S from the Fourier transformed input and output data u(w) and y( w). The system is assumed to be free of proc ess noise. In contrast to time domain methods, additive constants (bias) have not to be estimated, since initial values x(t=O) and measurement zero shifts do not influence the tr a nsformed data for frequencies "'k > O. First identification runs performed using this ba si c approach le a d to disappointing results and showed the need of a number of modifications.
1;711 Modified Method. The main modifications applied to the basic approac h consider the following items: - Extension of model equations for nonperiodic signals - Multi-Run-Evaluation-Technique - Combination of equation-error and maximum-likelihood-methods - Verification of results in time domain.
lOO
lO
HAX. (RROR Of O£RtVA livES.
MEAN ERROR OF O(RIVA TlVES %
%
100
10
100
10
The modified method is shown schematically in Fig. 1. Some characteristic details of this method are discussed in the following.
DATA EXTRACTION FOURIER TRANSFORMATION
1l"
SELECTION OF FREQUENCIES
1l"
MODE~
WITHOUT
!
t.x
WITH
1l"
MQDE~
t.x
WITHOUT
MODE~
t.x
WITH
t.x
FIRST IDENTIFICATION EQUATlON-ERROR-METHOD
l
Fig. 2
Maximum and mean error of 21 identif ied derivatives for two mathematical models (Equation (2) without c'x, Equation (5) with ill<.).
jtJ. x(w)
=
FINAL IDENTIFICATION M- l - OR OUTPUT-ERROR - METHOD
y( w) = Cx( w)
!
ii
where
PLOT OF OUTPUT SPECTRA
Concept of the DFVLR-Frequency Identification Method.
Domain
Extension of model equations. For th e basic proach, it was assumed that the signals can regarded as pe riodic, so that x(T) = x(O) ~(",) = jwx(w). In general, real test data do meet this condition. As shown in the Appendix,
Fourier transform of
~(w) = jwx + ill< =
[x(T -
Co n sequently modified to
x is
ill< e jwC,t/2
with
the state Equations
jtJ.x(w) = Ax(w)
+ Bu(w) -
apbe and not the
(3a)
C,t/2) - x(-C,t/2) ]/T
(3b)
(2a) have
ill< e jWC,t/2
=
+ S;(w)
(5)
+ Bil( w)
lB, -ill<. ],
D=
lD,
0], u
=
lu, e jwC,t/2]T.
During testing of the method with simulated helicopter flight data, both the unmodified Equation (2) and the modified Equation (3) were used. The accuracy of the identified models was compared using the errors of the 21 identified derivatives. As shown in Fig. 2, a satisfactory accuracy could be achieved only by using the modified Equation
BIAS IDENTIFICATION AND SIMULATION OF TIME HISTORIES
Fig. 1
Ax(w)
to be
(3) •
Multi-Run-Evaluation. For the identification of complex systems it is necessary to use data which contain information about each of the eigenmotions and about the effectiveness of each of the controls. In aircraft identification, it is not possible to obtain all this information from one single manoeuvre, because the test pilot is not able to excite several controls Simultaneously in a prescribed manner. Therefore it is necessary to use information from different runs that were flown independently. For this, the time records suitable for identification are selected (with equal length) and their Fourier transformed data are combined in one identification by minimizing the multi-run cost function
(4) J =
The additional term in Equation (4) can be treated in a simple manner by adding a fictitious control variable and setting the val ue of this variable to e j ",C,t/2 in the frequency domain. The corresponding additional column of the control matrix B has to be filled with the elements of the vector c'x or with unknown parameters to be identified. The extended Equations (2) then are
l. L £*("'k' i
i) S-1 £("'k' 1)
+ log I S I
(6)
k
(i = run number) instead of (1). (Note, c,x-term of Equation (3) depends on the ber.)
that the run num-
(i71
Fn:qul'ncv iloll1ain Two-Step-Identification. The advantages of two identification methods are combined. First, a robust Least-Squares-Equation-Error-Method is applied, which quickly provides preliminary results without requiring a-priori-values of the parameters. Then, a M-L-method is used to further improve these results and to obtain unbiased estimates.
jux
with state vector
x
=
tu, v, w, p, q, r,
4>,
control vector
u
=
[6 co ll '
6
]I
6y '
x'
6 TR
MAGNITUDE
I y(w) I
10-4
Fig. 3 presents the magnitudes of output spectra from the measured data and the identified model, Fig. 4 the verification of the model by time history plots. From both, it can be seen tha t the frequency domain method applied works satisfactory. Some discrepancies, especially in the plots of yaw rate could be observed also in time domain identifications and could be traced back to deficiencies of the flight test instrumentation system (Ref. 11, 12).
I PITCH RATE 0
10- 2
..
10- 6 10-3 10-4
RUN 3 10- 5
10- 6 10- 2
YAW RATE
~
RUN 2
10-7 10- 2
~ ..
10- 6 10-3
~
10- 4
10- 6
~. . .. . ~ . ..
10- 3
10- 4
r--\
10-5
10-7 10- 2
10- 4
1010-3
10- 5
.~
10- 4
~
~..
10- 5 10- 6 10- 2
10- 5
10- 7
~
10- 3 10- 3 10- 4
RUN 1 10-4 10- 5 0
Fig. 3
15 VOL l-W
.~ 10- 6 w
2
rid/lee
4
el T
performed using data from four manoeuvres flown at a trim speed of V = 150 km/h. These manoeuvres consist of different inputs into collective, longitudinal, lateral and tail rotor control (6 011' 6 , 6 , 6TR ). For the evaluation, only data c for t~e lirst eleven frequencies were used, since the power spectra of these data decrease rapidly for higher frequencies.
10-' RUN 4
el T
A multi-run-evaluation was
The frequency domain equations are
0
(7)
lax , ay, a z , u, v, w, p, q, r, 4>,
Data from a flight program with the Bo-105 research helicopter of the DFVLR were used for the identification of a six-degrees-of-freedom model containing 40 unknown state and control derivatives.
RATE
6x e jwllt / 2
measurement vector y
3. EXAMPLE 1
I ROLL
+ Bu{w) +
yew) = Cx(w) + Du(w)
Comparison of time histories. To verify the validity of the identified model in time domain, it is necessary to perform an integration of the differential equations, taking into account unknown initial states and zero shifts. For this, after identification of the system parameters (elements of A, B, C, D) an additional identification of bias-parameters only (with system parameters fixed) is performed using a time domain M-Lmethod. In this way, the time response plots of Figures 4 and 6 are obtained.
10- 2
Ax{",)
+ nonlinear gravity and inertia terms
0
~ w
2
rid ,"ee 4
w
2
rld'".c 4
Example I. Magnitude of output spectra from Bo-105 flight test data (--) and identified model (xxx) for the rotational degrees of freedom (combined identification from the 4 data run s shown in Fig. 4).
1\ 1. I\larchalld alld I'll K-h
RUN 1
RUN 2
RUN 3
RUN 4
20 YAW RATE
deg/sec
o
-20-L--------~~----~----L---~--~~~
______~
70 PITCH RATE
deg/sec
o
-70-L--~~--~~----~----L---~----~~------~
140 ROLL RATE
deg/sec
o
-140 +
YAW CONTROL
-]
_:~N .
LATERAL CONTROL
LONGITUDINAL CONTROL
COLLECTIVE PITCH CONTROL
Fi g. 4
FLIGHT TEST DATA + MODEL OUTPUT
I
~
I
I
-::~ j :]-, , L30sec
+
~: ~
I I
t ,
TIME
: ~
=
I
Example I. Time his tori es of i npu t signal s , meas ur ed output signals a nd comput ed outpu ts of th e identified mode l (Bo 105 Research Helicopter)_
fi73
Frequ e nc y Domain
INPUT VECTOR
LINEAR MODEL
TIME-DELAY MATRIX
Tou
~
To(w.e,.e 2• ... e m )
u
Fi g. 5
M(w.T,.T2 •.. ·T,)
OUTPUT
=
VECTOR y-MTou
Fr equ e ncy domain i denti fica t ion model with time-delays.
4. IDENTIFICATION OF SYSTEMS WITH TIKEDELAY
(8)
y - M To u aT
~-= M --11. u a0 i
a0 i
In contrast to time domain representation. the frequency domain representation allows a simple way of treating models with time-delays. As shown in Fig. 5 an output time-delay matrix M can be defined. describing the influence of one or more delay pa~ameters Ti • M is a diagonal matrix with
WTi T
dX
dTi- -
0
u
( i
m)
(9a)
( i -
r)
(9b)
As example. if only one time-delay is to be ident if ied and with M - e - jWT I. the modified equations are
-JOlT
Mkk - e i. i f Yk is to be delayed by Ti and Mkk - l. if Yk is not to be delayed. In a similar way, an input time-delay matrix can be defined, if the input variables have to be delayed instead of the output variables. The identification method described earlier can be used to identify the extended par,meter vector l01' 02' ••• 0. Tl' T2' ••• ] , when the model equations and t~e sensitivity equations are modified as follows
y - e-
ro-- i a
°-
j WT
To u
= e -j WT
~ _ -j aT
W
aT 'did\!. u
(i -
i
1 ••• m)
Y
PITCH STICK INPUT
r PITCH RATE
deg/sec
IDENTIFICATION WITHOUT TIME-DELAY
0
..... .... .. 4~. . .. + . . . . . ....
o
-4~--~----~--~--~~--~--~--
.
__~__~__~__- - J
[ IDENTIFICATION WITH TIME-DELAY Tq
0
PITCH RATE
deg/sec
o
Fi g. 6
6
12
18 TIME
24
sec
Example 2. Identification of a seco nd order model with time-delay from flight test data (DFVLR-HFB 320 In-Fl.ight Simu lat or. --- measur ed data. +++ output of identified mod el ; identified time=~elay Tq = .8 sec) . Model transfer function = K
I +
q I + a
l
Ts s +
2 e 8
2
s
q
30
(10)
(i71 5. EXAMPLE 2 Flight tests with the In-Flight-Simulator HFB 320 of the DFVLR were conducted for handling qualities investigations with defined time-delays between the commanded control inputs (pitch stick, roll stick) and the corresponding control surfaces (elevator, rudder). The method described above was used to identify a so-called equivalent model describing the pitch rate response q due to pitch control o. Fig. 6 shows the outputs of two identified models in comparison with the measured data. In both cases, four parameters of the secondorder-model have been identified. A good fit of pitch rate could be obtained when a time-delay Tq was identified as fifth parameter.
6. CONCLUSIONS
A frequency domain identification method for linear systems was developed using the maximum-likelihood-approach. The practical applicability of the method was improved by implementing a number of modifications to the basic approach and could be demonstrated during the evaluation of helicopter flight test data.
8. Gupta, N.K.: New Frequency Domain Methods for System Identification. Joint Automatic Control Conference, 1977. 9. Klein, V., and D.A. Keskar: Frequency Domain Identification of a Linear System Using Maximum Likelihood Estimation. Proceedings of the 5th IFAC Symposium, 1979. 10. Klein, V.: Maximum Likelihood Method for Estimating Airplane Stability and Control Parameters from Flight Data in Frequency Domain. NASA TP 1637, 1980. 11. Kaletka, J.: Practical Aspects of Helicopter Parameter Identification. AIAA 11th Atmospheric Flight Mechanics Conference, Seattle, 1984. 12. Fu, K.-H., and M. Marchand: Helicopter System Identification in the Frequency Domain. Ninth European Rotorcraft Forum, Stresa, 1983. -----
8. APPENDIX In addition, the method was extended to the identification of systems with one or more time delay parameters. This approach avoids the disadvantages associated with the Pade-approximation of time delays (increased model order, limited frequency range). First applications of the extended method to flight test data have been successful.
The discrete Fourier Transform of x(t) is defined as x(w) = l
x(w) P., et a1.: Parameter Identification. AGARD Lecture Series No. 104, 1979.
2. Maine, R.E., and K.W. Iliff: Identification of Dynamic Systems. AGARDograph, 1984. 3. Greenberg, H.: A Survey of Methods for Determining Stability Parameters of an Airplane from Dynamic Measurement. NACA TN 2340, 1951. 4. Shinbrot, M.: A Least-Squares Curve Fitting Method with Application to the Calculation of Stability Coefficients from Transient Response Data, NACA TN 2341, 1951. 5. Levy, E.C.: Complex Curve Fitting. IRE-Trans. Autom. Control 4, pp. 37-43, May 1959. 6. Marchand, M.: The Identification of Linear Multivariable Systems from Frequency Response Data. Proceedings of the 3rd IFAC Symposium, 1973. 7. Marchand, M., and R. Koehler: Determination of Aircraft Derivatives by Automatic Parameter Adjustment and Frequency Response Methods. AGARD Conference Proceedings No. 172, 1974.
1.
xk ;;jukllt
(AI)
The sum can be approximated by an integral and evaluated by partial integration
7. REFERENCES 1. Hamel,
N-l
N k=o
=
t T-lI~/2 -lI /2
x(t) e- jwt dt
(A2 )
. T-lI~/2 fl!!
x(t) e- jwt dt -lI /2 + lX(T-lIt/2)e- jw(T-lIt/2) _ x(-lIt/2)e- j w(-lIt/2) l/T.
=
jwT Since e - 1 for all frequencies used in Eq.(Al), the Fourier Transform of x(t) is x(w) = jux(w) + !Ix ejwllt/2
(A3)
where !Ix = [x(T-lIt/2) - x(-lIt/2) lIT !Ix can be calculated from the sampled data by the linear interpolation
(A4)
which requires two extra data points, x_I and xN' not used in Eq.(Al).