- Email: [email protected]
THE EXTENDED KALMAN FILTER IN THE FREQUENCY DOMAIN FOR THE IDENTIFICATION OF MECHANICAL STRUCTURES EXCITED BY SINUSOIDAL MULTIPLE INPUTS
Mechanical Systems and Signal Processing (2000) 14(3), 327}341 doi:10.1006/mssp.1999.1256, available online at http://www.idealibrary.com on
THE EXTENDED KALMAN FILTER IN THE FREQUENCY DOMAIN FOR THE IDENTIFICATION OF MECHANICAL STRUCTURES EXCITED BY SINUSOIDAL MULTIPLE INPUTS RICCARDO PROVASI AND GIAN ANTONIO ZANETTA ENEL, Struttura Ricerca, Area Generazione, Via Reggio Emilia 39-20090 Segrate (Mi), Italy. E-mail: [email protected]; [email protected] AND
ANDREA VANIA Politecnico di Milano-Dipartimento di Meccanica-P. za Leonardo da Vinci, 32-20133 Milano, Italy. E-mail: [email protected] (Received 3 February 1999, accepted 28 July 1999) A modal parameter identi"cation method applied to mechanical structures excited by correlated sinusoidal multiple inputs was developed. The algorithm is based on the same formulation of the extended Kalman "lter, applied as a system parameter identi"er in the frequency domain to mechanical structures subject to excitations characterised by an inherently high degree of correlation. The algorithm was validated by using simulated data on a multi-degree-of-freedom system. The tests demonstrate that the proposed technique is of practical application value. The method is devoted to the identi"cation of the modal parameters of supporting structures of rotating machinery, using data obtained during the normal operation of the machines. Nonetheless, it can have an interest for more general applications in the "eld of the experimental modal analysis. This study was undertaken in the framework of the BRITE EURAM III project MODIAROT (MOdel based DIAgnosis of ROTors in power plants). 2000 Academic Press
1. INTRODUCTION
The dynamic behaviour of large and complex mechanical structures can be studied using "nite element (FE) methods. But for such structures, reliable FE models are costly and di$cult to set up, because of the complexity of the structure and, often, due to lack of necessary information. Moreover, they could be sometimes unnecessarily detailed and, in any case, they should be adequately tuned. Therefore, a reasonable goal is to model the dynamics of these structures by means of simpli"ed and yet accurate equivalent models, identi"ed on the basis of experimental data obtained from tests with arti"cial excitations or better still during the normal operation of the machine of which the structure could be a component. The identi"cation methods in common use require the knowledge of the transfer functions of the structure at its measurement points with respect to any single excitation point. Large structures may require arti"cial multiple-input excitations to increase the energy input. In order to evaluate the transfer functions, the exciting forces must be generated by uncorrelated random signals. In this case, special algorithms of signal analysis have to 0888}3270/00/050327#15 $35.00/0
2000 Academic Press
328
R. PROVASI E¹ A¸.
be applied for separating the contribution of each single source. Alternatively, multiple sinusoidal excitations can introduce more energy in any single-frequency band and can better point out non-linearities in the structure, but require appropriate control systems to drive two or more exciters in order to obtain predetermined force shapes. However, this requires rather sophisticated equipment for the control of the relative phases and amplitudes of the generated forces and the synchronisation with the data acquisition system. On the other hand, during the normal operation of the machines, often the excitations generated by some component of the machine or exchanged between di!erent components at the connection points are intrinsically correlated, thus preventing a straight determination of the transfer functions of the system. Most times these forces are sinusoidal or can be reduced to the harmonic components of the working frequency of the machine, as in the case of rotating machinery. Whenever the operational forces exciting a structure can be directly measured or indirectly evaluated, or when only simple excitation systems can be used, the development of identi"cation procedures based on the analysis of the frequency responses of structures excited by correlated multiple inputs can be very useful. In view of applying some identi"cation procedure, the dynamics of a structure or a substructure could be described by the sti!ness, damping and mass matrices, condensed at a reduced number of physical degrees of freedoms (dofs), possibly associated, in the case of substructures, with the connection points with adjacent components. However, serious theoretical problems are encountered and, in any case, the number of parameters to be identi"ed may become exceedingly large. This problem arose in particular in the "eld of rotating machinery, where the identi"cation of the supporting structures has been pursued by several authors for the last decade and di!erent identi"cation schemes were proposed [1}4]. All these techniques could be classi"ed as direct system parameter estimation methods [5]. The advantage of identifying the direct parameters is due to the fact that the inverse problem of computing the forces in di!erent conditions from the measured vibrations could be easily solved, once the direct parameter model is identi"ed. But, apart from the comparative robustness of these identi"cation procedures, it is authors' experience that the identi"ed matrices are in general non-unique and they are likely to be faced with a set of incomplete data [6]. The natural alternative is to resort to modal models, that have the additional advantage of making it possible to use also incomplete models for the substructures in order to build up the total system model. On the other hand, the possibility of solving the inverse problem is lost. Having stated the problem in terms of modal parameter estimation in the frequency domain, it was necessary to devise a method capable of carrying out the identi"cation with vibration data generated even by correlated multiple inputs, in particular sinusoidal inputs, without any a priori constraints between them. In fact, the commercially available modal analysis software does not deal with such a problem. In the literature, applications of sinusoidal multiple-input excitations are described in connection with the well-known appropriation techniques [7] or, more recently, with multiphase-step-sine methods [8], where the sources are separated by appropriate excitation sequences with force patterns having pre-established amplitudes and phases. In a di!erent approach, the modal analysis of the response to sinusoidal multiple-input excitations is performed by a gradient method, but once again, the phase and amplitude relationships of the forces are "xed [9]. In order to overcome these restrictions, an attempt was made to generalise the last method allowing for a generic set of sinusoidal forces [10], but the results obtained were
EXTENDED KALMAN FILTER
329
erratic and often the algorithm gave divergent solutions, being exceedingly sensitive to initial parameter values and noise. Additonally, some preliminary trials with other approaches provided by a MATLAB] toolbox did not give satisfactory results. Much more promising was an approach based on the same formulation of the extended Kalman "lter (EKF) theory, applied in the frequency domain [11]. In fact, the problem stated above requires the identi"cation of the modal parameters of a linear, time-invariant mechanical system. The non-linearity arises in the equation of the observations as a function of the state vector to be identi"ed, that is to say the unknown modal parameters. In general, one of the possible application of the Kalman "lter is to use it simply as a system identi"er. In the paper, it is clearly shown that the problem of the identi"cation in the frequency domain can be solved by a formulation coincident with the one of the EKF in the time domain, by a simple substitution of variables (frequency instead of time), where the usual assumptions for the application of the "lter about the observation and state noise hold. This method has been developed in the framework of the BRITE-EURAM research project MODIAROT (MOdel based DIAgnosis of ROTors in power plants) and it was particularly devoted to the above-mentioned problem of identifying rotor foundations. In the present paper, the formulation of the algorithm and some practical aspects of the implementation are discussed. The results obtained by the validation tests on simulated data for a multi-degree-of-freedom system are then discussed.
2. THEORETICAL DEVELOPMENTS
In what follows, the fundamentals of the extended Kalman "lter (EKF) theory are shortly recalled "rst. Then the application of the "lter in the frequency domain to the modal model identi"cation of mechanical structures excited by correlated sinusoidal multiple inputs is illustrated. 2.1. THE EXTENDED KALMAN FILTER (EKF) With a notation familiar to the control theory, the dynamics of a discrete-time system is described by the following equations [12]: x
"f (x , u )#g (x )w I> I I I I I I z "h (x )#v . (1) I I I I In general, the state vector x is unknown and the observation variables z are known I> I at the discrete times t . The state and observation noise, in the vectors w and v , are random I I I processes, i.e. it is assumed that they are independent, white, Gaussian, with zero mean and non-negative covariance matrices. Moreover, let us assume that the state is independent of the input u . I The extended Kalman "lter is a recursive algorithm developed to estimate the state vector on the basis of noisy measurements of the input and output signals. After linearization of equations (1) with the "rst order of the Taylor series in the proximity of the estimate xL , it is possible to demonstrate that the state estimate at time t from the observations I I I> up to time t , is I xL "f (xL )#K (z !h (xL )) (2) I> I I I I\ I I I I I\ and the observation estimate is zL "h (xL ) I I I I\
(3)
330
R. PROVASI E¹ A¸.
Figure 1. Flow chart of extended Kalman "lter.
where the "lter gain matrix K is de"ned as follows: I K "F P H2 (H P H2#R )\. I I I I\ I I I I\ I I
(4)
In equation (4), P is the covariance matrix of the state estimate error at time k, based I I\ on the estimate up to time k!1 and R is the covariance matrix of the observation noise at I time k. In Fig. 1, a block diagram shows the "lter working scheme, with the measurable response z from the real world and its estimate zL computed by the "lter. I I The predicted and "ltered state together with its error covariance matrix, can be computed at any time via a set of recursive formulas that give the "lter correction: xL "xL #P H2 (H P H2#R )\(z !zL ) I I I I\ I I\ I I I I\ I I I I P "P !P H2 (H P H2#R )\H P I I I I\ I I\ I I I I\ I I I I I\ xL "f (xL ) I> I I I I P "F P F2#G Q G2 . I> I I I I I I I I
(5)
The state estimate at time t from the observation up to time t is computed starting I> I from an initial guess of the state vector xL , considered as a random variable with given mean value and covariance matrix and uncorrelated to the noise on the state and on the observations. At the end of the process, the "lter provides a suboptimal estimate of the actual state. 2.2.
STATEMENT OF THE PROBLEM OF THE MODAL MODEL IDENTIFICATION OF MECHANICAL
STRUCTURES
For a mechanical structure excited by sinusoidal correlated multi-point inputs, the steady-state equation of motion can be written as !Muz#jCuz#Kz"F.
(6)
EXTENDED KALMAN FILTER
331
When in equation (6), the viscous damping matrix C is assumed to be proportional to the mass and sti!ness matrices, the resulting modes are real. This hypothesis is not very restrictive in practice and, anyway, can be removed, without a!ecting the kernel of the developed method. The response at the ith dof of the structure for any particular frequency in a given frequency range can be expressed as a modal parameter expansion:
t2 F(u) t2F(u) , t2 F(u) P #tG #tG S . zG(u)" tG !u S u P (u!u#2jm u u) P P P S P
(7)
In equation (7), the modes are real and normalised to the unit real mass. The last two right-hand side residual terms are used to approximate the contributions of the lower and upper modes at frequencies below or above the range of analysis. For the upper modes, the values of the pseudo-eigenfrequency u , is chosen somewhat arbitrarily and it is set high S enough so that u and any possible term 2jm u u can be neglected. S S If the response zG(u) is measured at least at all the points where a force is applied and the forces are known, directly or indirectly, from measurements, the unknown modal parameters u , m , t , t , t can be identi"ed from equation (4) without the need to start P P P J S from a transfer function, which cannot be computed in case of correlated multiple inputs. It should be emphasised that it is not necessary to have forces on all the response points in order to identify completely t , t , t , so the identi"cation procedure can be applied to P J S a more general experimental modal analysis problem. Another point to be outlined is that the excitation forces need not necessarily be sinusoidal. Both the forces and the responses in equation (7) could be the result of a simple FFT of random signals performed on each single channel of acquisition, without any need to separate the contribution of multiple-input random forces. 2.3.
APPLICATION OF THE EKF TO MODAL PARAMETERS IDENTIFICATION IN THE
FREQUENCY DOMAIN
Equation (7) is non-linear in the unknown parameters. The extended Kalman "lter algorithm can be applied in the frequency domain [11, 13] to the modal parameter identi"cation problem stated in equation (7), with the state vector de"ned as x"+t2,2, t2,2, t2 , m ,2, m ,2, m , u ,2, u ,2, u , t2 , t2,2. S P , P , P , J
(8)
The system equations are x "x I> I z "h (x )#v I I I I
(9)
where z , x and v are referred to the kth forcing frequency u . The function h has the I I I I I same expression as the term on the right-hand side of equation (7) and, therefore, is a non-linear function of the state vector x , the input forces and the frequency u . The I I usual hypotheses of random, white, Gaussian, zero mean and independent noise can be applied as well to the noise on the observations in the frequency domain, while no noise is applied to the state vector, because it is made up by deterministic, constant modal parameters.
332
R. PROVASI E¹ A¸.
The identi"cation technique based on the EKF is a predictor-corrector method of the state estimate and of its corresponding error covariance matrix, working through the following steps: E prediction at the frequency u , using the frequency responses measured at frequencies I> up to u : I xL "xL I> I I I P "P (10) I> I I I E correction using the frequency responses up to u : I> xL "xL #K (z !h (xL )) I> I> I> I I> I> I> I> I P "(I!K H )P (I!K H )2#K R K2 (11) I> I> I> I> I> I I> I> I> I> I> where the "lter gain is given by K "P H2 (H P H2#R )\. (12) I I I\ I I I I\ I I At the end of the iteration cycle, one obtains the estimate xL of state vector and the associated con"dence intervals for each component of xL through its covariance matrix P. 2.4. NUMERICAL PROBLEMS In applying the "ltering procedure to a speci"c system, a dynamic model must be de"ned and noise statistics and initial guesses for the state must be speci"ed. In general, the model is only an approximation of a physical system; the model parameters and the noise statistics may be uncorrect. It has already been outlined that the initial guesses for the unknown parameters can be a rough approximation of the true values and are often speci"ed somewhat arbitrarily when starting the "ltering process. The factors can a!ect the "lter performance and make the process to diverge or, even worse, to converge to unreliable results. The last case can be detected by inspecting the values in the covariance matrix P I I or applying some of the traditional modal model validation techniques. Considering the identi"cation problem stated above, three factors primarily in#uence the convergence of the iterative procedure: E The values of the initial system parameters. E The condition number of the initial error covariance matrix, determined by the relative order of magnitude of the di!erent parameters to be identi"ed. As a matter of fact, the "lter could become insensitive to the parameters with the smallest values. E The condition number of the matrix (HPH2#R), which is updated at each step during the "ltering process. This matrix must be non-negative in order to obtain reliable results from the numerical inversion. As for the initial value of the parameters, it is in general not di$cult to give a rough approximation for the eigenvalues, the error being usually much higher for the damping than for the frequency. It is more di$cult to give a reasonable guess for the order of magnitude of the modes, since they are normalised to the unit mass. Therefore, a devoted algorithm capable of working out automatically from the measured frequency responses a "rst approximation for the mode shapes was developed. To solve the convergence problem, several mathematical techniques have been proposed [14]. In particular, the weighted global iteration technique [15] appears to be more e!ective in reducing the in#uence of the initial values. According to this method, a "rst iteration is carried out with the assigned initial guesses to obtain the converged values after the
EXTENDED KALMAN FILTER
333
processing of all the observations. Any subsequent iteration is performed starting from the "nal values of the previous step, but overweighting the error covariance matrix, until the convergence is obtained or the number of iterations is exhausted. While it was veri"ed that the algorithm can stand a high noise rate on the measured displacements, the noise rejection on the measured forces can be a problem. In its current implementation, the algorithm does not have an explicit formulation for counteracting the noise in the input forces, although it was veri"ed that it can nonetheless stand a signi"cant random noise rate on the forces, as will be discussed later. The practical implementation of the algorithm allows for the choice of a subset of measurement points to be considered in a "rst identi"cation step, i.e. those for which the modes are more clearly excited. The process is then repeated for the remaining points, possibly holding some parameters to the values estimated in the previous step. In particular, but not necessarily, this also applies to the identi"cation of the components of the modes for the points where no forces are acting.
3. EXAMPLE OF APPLICATION
In order to validate the proposed method, a "rst series of tests was performed on a lumped mass}spring system with eight dofs, with viscous damping proportional to the mass and sti!ness matrices so as to obtain real normal modes. The spring connections of each mass with up to "ve other masses were included. A drawing of the system is shown in Fig. 2. This system was set up with the aim of simulating the behaviour of the supporting structure of a complex rotating machine, such as a turbogenerator unit. The resulting mass and sti!ness matrices are reported in Table 1. The damping rate was assigned by proportionality coe$cients of 0.03 and 0.0033 with respect to the mass and the sti!ness values. The reference eigenvalues and eigenvectors were computed by standard routines (Table 2). The synchronous forced response was then obtained over a frequency range 5%65 Hz, thus encompassing all the eigenfrequencies; the frequency resolution was 0.25 Hz. In Fig. 3 it is possible to have an overview of the exciting forces and the relevant frequency responses. The resulting curves are similar to those provided by the monitoring systems during run-ups and shut-downs of turbogenerator units. The exciting forces were essentially a combination of components constant and proportional to the square of the frequency, in order to simulate the e!ect of rotating forces on a real rotor. No particular care was taken in order to obtain the best modal excitation on all modes. The correlation coe$cients between the forces are very close to one.
Figure 2. Drawing of the eight dofs lumped mass}spring system.
334
R. PROVASI E¹ A¸.
TABLE 1 Mass and sti+ness matrices of the eight dofs system Diagonal mass matrix (kg) 2.0
89 400
1.8
!36 000 80 000
Symmetric
1.6
!12 000 !30 000 99 850
1.2
1.7
0.8
Sti!ness matrix (N/m) !3600 !1200 !600 !10 000 !3000 !1000 !39 000 !13 000 !3900 101 000 !33 000 !11 000 107 200 !40 000 92 500
1.2
!1300 !3300 !13 000 !27 000
1.0
!650 !1100 !4000 !9000 30 000 64 750
With these data it was possible to validate the method against problems such as high correlation between the input forces, a comparatively high number of modes to be identi"ed, signi"cant modal coupling due to the damping and presence of local modes, i.e., modes with signi"cant components only on two or three dofs. Moreover, the e!ect of the noise was investigated by addition of a uniformly distributed random error on the measurements and the forces, proportional to the respective maximum amplitudes among all the dofs. Up to eight modes were identi"ed at a time. The MATLAB] version of the code took 2900 s of CPU time with noisy data on a 200 MHz Windows NT PC computer. The number of modes did not seem to be a problem for the algorithm. In some cases, working on a reduced frequency range, modes external to the frequency range of analysis could still be identi"ed with su$cient, sometimes even good, accuracy. The initial values of the eigenfrequencies were taken with a 10% error and the initial modal damping with approximately 100% error (see Table 3). For these tests, the initial values of the mode shapes were generated as a random combination of the true ones, with errors on the single component even greater than 100% and resulting shapes of the modes far away from the correct ones (see Fig. 4). With no noise on the data, almost the exact computational values were obtained for all the parameters. More than acceptable approximations were obtained also with 10% of random error, as de"ned above, on the vibration response. As expected, the method was much more sensitive to the noise on the input forces. At most a 4% error could be added to the forces on this set of data, in order to obtain acceptable results. The results for the identi"ed eigenvalues with a 4% error on both input and output, as in the curves of Fig. 2, are reported in Table 2, where they are compared with the computed values and the values assigned at the "rst iteration step. The only remarkable deviation is on the damping of the "rst eigenfrequency. The amplitudes of the "rst mode components were correspondingly underestimated. As observed from Fig. 3, the signal-to-noise ratio for the forces is very poor at low frequencies, where it is di$cult even to guess a phase value and the energy of the noise is comparable with that of the forces themselves. The identi"ed mode shapes were well correlated with the reference modes obtained by computation, as illustrated in Table 3 by the Modal Assurance Criterion (MAC) coe$cients.
Computed u (Hz) A Identi"ed u (Hz) A Initial u (Hz) A Computed f (%) A Identi"ed f (%) A Initial f (%) A 23.10 23.10 20.79 2.43 2.49 5.40
1.02 1.13 2.27
Mode C 2
9.51 9.51 8.56
Mode C 1
3.53 3.48 7.85
33.64 33.65 30.28
Mode C 3
4.14 4.18 9.20
39.42 39.39 35.48
Mode C 4
4.50 4.39 10.00
42.82 42.82 38.53
Mode C 5
TABLE 2 Computed, identi,ed and initial eigenvalues
5.03 5.00 11.18
47.87 47.87 43.09
Mode C 6
5.51 5.49 12.24
52.40 52.42 47.16
Mode C 7
6.22 6.28 13.82
59.12 59.11 53.21
Mode C 8
EXTENDED KALMAN FILTER
335
336
R. PROVASI E¹ A¸.
Figure 3. Exciting forces and frequency response for the eight dofs system with 4% random error.
In Fig. 4 the whole set of the reference eigenvectors, given by numerical computations, are compared with the initial values assigned at the "rst iteration step and with the identi"ed mode shapes. The overall quality of the results can be evaluated from Fig. 5, where the simulated measurements are compared with the synthesized curves obtained from the identi"ed modal parameters and from their initial values. The curves are shown for three di!erent dofs, as amplitude and phase of the frequency response in a range not encompassing the "rst eigenfrequency. Apart from the "rst mode which was signi"cantly underestimated for the reasons discussed above, for the other modes the curve "tting was very good. As a crosscheck, the curve "tting was then repeated with a set of data generated by forces di!erent from those used in the identi"cation process, with similar results. In order to ascertain the quality of the results, an error index was computed for every single dof in the system as the per cent ratio of the rms values of the error between the curves synthesized by the identi"ed parameters and the rms values of the simulated measurements. A global error index was then obtained as the average error among all the dofs. The same error indices were de"ned with respect to the curves obtained from the initial values assigned to the parameters, in order to appreciate the initial distance from the simulated measurements. The indices are reported in Table 4. Separate values are given in case the frequency range of the "rst mode is taken into account or neglected. Finally, it can be mentioned that the e!ects of having pseudo-resonances in the measurements were also tested. This is the case where there are forces acting at connection points with other components, having their own resonances. Successful results were obtained also with these data sets.
4. CONCLUSIONS
The application of a method based on the same formulation of the extended Kalman "lter and applied in the frequency domain proved to be a suitable tool for solving the problem of the time-invariant modal parameter identi"cation of a mechanical structure excited by sinusoidal multiple inputs having an inherently high degree of correlation.
Mode Mode Mode Mode Mode Mode Mode Mode
C C C C C C C C
Identi"ed
1 2 3 4 5 6 7 8
1.000 0.041 0.001 0.001 0.001 0.000 0.006 0.029
Mode C 1 0.036 1.000 0.035 0.000 0.000 0.001 0.000 0.019
Mode C 2 0.000 0.034 1.000 0.017 0.001 0.006 0.007 0.002
Mode C 3 0.001 0.000 0.014 1.000 0.010 0.007 0.000 0.011
Mode C 4 0.001 0.000 0.001 0.011 1.000 0.013 0.003 0.014
Mode C 5 0.000 0.000 0.006 0.008 0.014 1.000 0.000 0.000
Mode C 6
0.005 0.000 0.007 0.000 0.003 0.000 1.000 0.004
Mode C 7
TABLE 3 Model assurance criterion coe.cients between identi,ed and computed modes
0.028 0.016 0.002 0.011 0.013 0.000 0.003 1.000
Mode C 8
EXTENDED KALMAN FILTER
337
Figure 4. Mode shapes of the eight dofs system: 䊊*䊊*䊊, reference values; # *# *#, identi"ed values; *****, initial values.
338 R. PROVASI E¹ A¸.
339
EXTENDED KALMAN FILTER
Figure 5. Curve "tting with identi"ed and initial modal parameters at dofs C5}7.
TABLE 4 Error values in the curve ,tting with the initial and the identi,ed modal parameters Point 1 Point 2 Point 3 Point 4 Point 5 Point 6 Point 7 Point 8 Mean Initial % rms error} "rst mode included
89.1
84.8
99.8
99.6
92.4
105.6
101.9
88.1
95.2
Final % rms error} "rst mode included
81.1
26.1
29.1
21.9
30.1
18.0
20.7
13.4
22.2
Initial % rms error} "rst mode excluded
88.1
82.6
99.4
97.2
91.2
105.1
101.7
84.9
93.8
Final % rms error} "rst mode excluded
4.7
6.1
6.6
4.5
7.3
3.8
5.5
3.7
5.3
The present paper was aimed at illustrating the formulation of the algorithm and to stress its potential applications, rather than proving its theoretical bases, although it was clearly shown that the proposed formulation in the frequency domain is identical to the more familiar one in the time domain. It was also discussed that the usual hypotheses for the application of the EKF in the time domain are equally satis"ed at least for the considered identi"cation problem dealt with in the frequency domain.
340
R. PROVASI E¹ A¸.
The method was validated so far on di!erent sets of data simulated by a multi-degree-offreedom lumped parameter model, well representative of possible problems to be dealt with in a real situation. These tests, though not exhaustive of any possible occurrence, demonstrated that this approach gave considerably better results than other alternative approaches, in terms of convergency, sensitivity to initial guesses and consistency of the results with respect to di!erent data sets used for the identi"cation. In its current implementation, the method was successfully tested against such problems like strong coupling between di!erent modes, presence of local modes, strong di!erences between mode components measured on di!erent directions, pseudo-resonances in the exciting forces, and poor modal excitation. Up to eight modes were identi"ed at a time. The program was not particularly fast, but the software programming was not yet optimized. The initial values needed to such a recursive algorithm can be worked out from the experimental data ignoring any a priori knowledge, as it would be in a real case. No divergence problems were experienced with such initial guesses, although sometimes they were pretty distant from the true values, but more extensive tests should be performed. In some cases, the ability of the algorithm to identify modes outside the frequency range of analysis was also veri"ed. The random noise rejection was good for noise on the measured vibrations and still acceptable for noise on the input forces, although in the present revision there are no speci"c tools to reduce the e!ect of the last kind of noise. Future developments of the method should include tests against deterministic errors on the input forces, the study of algorithms to "lter the noise on the forces, a generalization of the algorithm for the identi"cation of complex modes, a revision of the iteration scheme that should further increase stability and speed of convergence and the optimisation of the software in order to reduce the computation time. Even in its current revision, the method seems to be appropriate for the application to the modal parameter identi"cation of supporting structure of rotating machines, for which it was primarily developed. Moreover, it could be applied as an additional tool in experimental modal analysis packages, to be used whenever multiple excitations are required in order to introduce more energy into the structure and when sophisticated input controls are not available.
ACKNOWLEDGEMENTS
This work was partially funded by the EU in the project MODIAROT BE95-2015 Contract No. BRPR-CT95-0022 of the research programme BRITE EURAM III.
REFERENCES 1. C. P. FRITZEN 1986 Journal of
EXTENDED KALMAN FILTER
341
6. G. A. ZANETTA 1992 IMechE International Conference on
APPENDIX: NOMENCLATURE M C K F z x I u I w I v I f I g I h I xL I I zL I I u P m P t P u S t J t S F I H I G I K I P I Q I R I
structure mass matrix structure damping matrix structure sti!ness matrix exciting force vector measured vibration vector state vector at time k input vector at time k noise vector on the state at time k noise vector on the observations at time k non-linear function of the state at time k non-linear function on the state noise at time k non-linear function of the observ. at time k state estimate of x based on time k I observation estimate of z based on time k I rth-mode eigenfrequency rth-mode modal damping rth-mode shape vector pseudo-eigenfrequency of the equivalent upper mode shape vector of the lower-mode residuals shape vector of the upper mode residuals matrix of derivatives of f at x I I matrix of derivatives of h at x I I matrix of derivatives of g at x I I "lter gain matrix at time k covariance matrix of the state estim. error at time k covariance matrix of the state noise at time k covariance matrix of the observ. noise at time k
THE EXTENDED KALMAN FILTER IN THE FREQUENCY DOMAIN FOR THE IDENTIFICATION OF MECHANICAL STRUCTURES EXCITED BY SINUSOIDAL MULTIPLE INPUTS RICCARDO PROVASI AND GIAN ANTONIO ZANETTA ENEL, Struttura Ricerca, Area Generazione, Via Reggio Emilia 39-20090 Segrate (Mi), Italy. E-mail: [email protected]; [email protected] AND
ANDREA VANIA Politecnico di Milano-Dipartimento di Meccanica-P. za Leonardo da Vinci, 32-20133 Milano, Italy. E-mail: [email protected] (Received 3 February 1999, accepted 28 July 1999) A modal parameter identi"cation method applied to mechanical structures excited by correlated sinusoidal multiple inputs was developed. The algorithm is based on the same formulation of the extended Kalman "lter, applied as a system parameter identi"er in the frequency domain to mechanical structures subject to excitations characterised by an inherently high degree of correlation. The algorithm was validated by using simulated data on a multi-degree-of-freedom system. The tests demonstrate that the proposed technique is of practical application value. The method is devoted to the identi"cation of the modal parameters of supporting structures of rotating machinery, using data obtained during the normal operation of the machines. Nonetheless, it can have an interest for more general applications in the "eld of the experimental modal analysis. This study was undertaken in the framework of the BRITE EURAM III project MODIAROT (MOdel based DIAgnosis of ROTors in power plants). 2000 Academic Press
1. INTRODUCTION
The dynamic behaviour of large and complex mechanical structures can be studied using "nite element (FE) methods. But for such structures, reliable FE models are costly and di$cult to set up, because of the complexity of the structure and, often, due to lack of necessary information. Moreover, they could be sometimes unnecessarily detailed and, in any case, they should be adequately tuned. Therefore, a reasonable goal is to model the dynamics of these structures by means of simpli"ed and yet accurate equivalent models, identi"ed on the basis of experimental data obtained from tests with arti"cial excitations or better still during the normal operation of the machine of which the structure could be a component. The identi"cation methods in common use require the knowledge of the transfer functions of the structure at its measurement points with respect to any single excitation point. Large structures may require arti"cial multiple-input excitations to increase the energy input. In order to evaluate the transfer functions, the exciting forces must be generated by uncorrelated random signals. In this case, special algorithms of signal analysis have to 0888}3270/00/050327#15 $35.00/0
2000 Academic Press
328
R. PROVASI E¹ A¸.
be applied for separating the contribution of each single source. Alternatively, multiple sinusoidal excitations can introduce more energy in any single-frequency band and can better point out non-linearities in the structure, but require appropriate control systems to drive two or more exciters in order to obtain predetermined force shapes. However, this requires rather sophisticated equipment for the control of the relative phases and amplitudes of the generated forces and the synchronisation with the data acquisition system. On the other hand, during the normal operation of the machines, often the excitations generated by some component of the machine or exchanged between di!erent components at the connection points are intrinsically correlated, thus preventing a straight determination of the transfer functions of the system. Most times these forces are sinusoidal or can be reduced to the harmonic components of the working frequency of the machine, as in the case of rotating machinery. Whenever the operational forces exciting a structure can be directly measured or indirectly evaluated, or when only simple excitation systems can be used, the development of identi"cation procedures based on the analysis of the frequency responses of structures excited by correlated multiple inputs can be very useful. In view of applying some identi"cation procedure, the dynamics of a structure or a substructure could be described by the sti!ness, damping and mass matrices, condensed at a reduced number of physical degrees of freedoms (dofs), possibly associated, in the case of substructures, with the connection points with adjacent components. However, serious theoretical problems are encountered and, in any case, the number of parameters to be identi"ed may become exceedingly large. This problem arose in particular in the "eld of rotating machinery, where the identi"cation of the supporting structures has been pursued by several authors for the last decade and di!erent identi"cation schemes were proposed [1}4]. All these techniques could be classi"ed as direct system parameter estimation methods [5]. The advantage of identifying the direct parameters is due to the fact that the inverse problem of computing the forces in di!erent conditions from the measured vibrations could be easily solved, once the direct parameter model is identi"ed. But, apart from the comparative robustness of these identi"cation procedures, it is authors' experience that the identi"ed matrices are in general non-unique and they are likely to be faced with a set of incomplete data [6]. The natural alternative is to resort to modal models, that have the additional advantage of making it possible to use also incomplete models for the substructures in order to build up the total system model. On the other hand, the possibility of solving the inverse problem is lost. Having stated the problem in terms of modal parameter estimation in the frequency domain, it was necessary to devise a method capable of carrying out the identi"cation with vibration data generated even by correlated multiple inputs, in particular sinusoidal inputs, without any a priori constraints between them. In fact, the commercially available modal analysis software does not deal with such a problem. In the literature, applications of sinusoidal multiple-input excitations are described in connection with the well-known appropriation techniques [7] or, more recently, with multiphase-step-sine methods [8], where the sources are separated by appropriate excitation sequences with force patterns having pre-established amplitudes and phases. In a di!erent approach, the modal analysis of the response to sinusoidal multiple-input excitations is performed by a gradient method, but once again, the phase and amplitude relationships of the forces are "xed [9]. In order to overcome these restrictions, an attempt was made to generalise the last method allowing for a generic set of sinusoidal forces [10], but the results obtained were
EXTENDED KALMAN FILTER
329
erratic and often the algorithm gave divergent solutions, being exceedingly sensitive to initial parameter values and noise. Additonally, some preliminary trials with other approaches provided by a MATLAB] toolbox did not give satisfactory results. Much more promising was an approach based on the same formulation of the extended Kalman "lter (EKF) theory, applied in the frequency domain [11]. In fact, the problem stated above requires the identi"cation of the modal parameters of a linear, time-invariant mechanical system. The non-linearity arises in the equation of the observations as a function of the state vector to be identi"ed, that is to say the unknown modal parameters. In general, one of the possible application of the Kalman "lter is to use it simply as a system identi"er. In the paper, it is clearly shown that the problem of the identi"cation in the frequency domain can be solved by a formulation coincident with the one of the EKF in the time domain, by a simple substitution of variables (frequency instead of time), where the usual assumptions for the application of the "lter about the observation and state noise hold. This method has been developed in the framework of the BRITE-EURAM research project MODIAROT (MOdel based DIAgnosis of ROTors in power plants) and it was particularly devoted to the above-mentioned problem of identifying rotor foundations. In the present paper, the formulation of the algorithm and some practical aspects of the implementation are discussed. The results obtained by the validation tests on simulated data for a multi-degree-of-freedom system are then discussed.
2. THEORETICAL DEVELOPMENTS
In what follows, the fundamentals of the extended Kalman "lter (EKF) theory are shortly recalled "rst. Then the application of the "lter in the frequency domain to the modal model identi"cation of mechanical structures excited by correlated sinusoidal multiple inputs is illustrated. 2.1. THE EXTENDED KALMAN FILTER (EKF) With a notation familiar to the control theory, the dynamics of a discrete-time system is described by the following equations [12]: x
"f (x , u )#g (x )w I> I I I I I I z "h (x )#v . (1) I I I I In general, the state vector x is unknown and the observation variables z are known I> I at the discrete times t . The state and observation noise, in the vectors w and v , are random I I I processes, i.e. it is assumed that they are independent, white, Gaussian, with zero mean and non-negative covariance matrices. Moreover, let us assume that the state is independent of the input u . I The extended Kalman "lter is a recursive algorithm developed to estimate the state vector on the basis of noisy measurements of the input and output signals. After linearization of equations (1) with the "rst order of the Taylor series in the proximity of the estimate xL , it is possible to demonstrate that the state estimate at time t from the observations I I I> up to time t , is I xL "f (xL )#K (z !h (xL )) (2) I> I I I I\ I I I I I\ and the observation estimate is zL "h (xL ) I I I I\
(3)
330
R. PROVASI E¹ A¸.
Figure 1. Flow chart of extended Kalman "lter.
where the "lter gain matrix K is de"ned as follows: I K "F P H2 (H P H2#R )\. I I I I\ I I I I\ I I
(4)
In equation (4), P is the covariance matrix of the state estimate error at time k, based I I\ on the estimate up to time k!1 and R is the covariance matrix of the observation noise at I time k. In Fig. 1, a block diagram shows the "lter working scheme, with the measurable response z from the real world and its estimate zL computed by the "lter. I I The predicted and "ltered state together with its error covariance matrix, can be computed at any time via a set of recursive formulas that give the "lter correction: xL "xL #P H2 (H P H2#R )\(z !zL ) I I I I\ I I\ I I I I\ I I I I P "P !P H2 (H P H2#R )\H P I I I I\ I I\ I I I I\ I I I I I\ xL "f (xL ) I> I I I I P "F P F2#G Q G2 . I> I I I I I I I I
(5)
The state estimate at time t from the observation up to time t is computed starting I> I from an initial guess of the state vector xL , considered as a random variable with given mean value and covariance matrix and uncorrelated to the noise on the state and on the observations. At the end of the process, the "lter provides a suboptimal estimate of the actual state. 2.2.
STATEMENT OF THE PROBLEM OF THE MODAL MODEL IDENTIFICATION OF MECHANICAL
STRUCTURES
For a mechanical structure excited by sinusoidal correlated multi-point inputs, the steady-state equation of motion can be written as !Muz#jCuz#Kz"F.
(6)
EXTENDED KALMAN FILTER
331
When in equation (6), the viscous damping matrix C is assumed to be proportional to the mass and sti!ness matrices, the resulting modes are real. This hypothesis is not very restrictive in practice and, anyway, can be removed, without a!ecting the kernel of the developed method. The response at the ith dof of the structure for any particular frequency in a given frequency range can be expressed as a modal parameter expansion:
t2 F(u) t2F(u) , t2 F(u) P #tG #tG S . zG(u)" tG !u S u P (u!u#2jm u u) P P P S P
(7)
In equation (7), the modes are real and normalised to the unit real mass. The last two right-hand side residual terms are used to approximate the contributions of the lower and upper modes at frequencies below or above the range of analysis. For the upper modes, the values of the pseudo-eigenfrequency u , is chosen somewhat arbitrarily and it is set high S enough so that u and any possible term 2jm u u can be neglected. S S If the response zG(u) is measured at least at all the points where a force is applied and the forces are known, directly or indirectly, from measurements, the unknown modal parameters u , m , t , t , t can be identi"ed from equation (4) without the need to start P P P J S from a transfer function, which cannot be computed in case of correlated multiple inputs. It should be emphasised that it is not necessary to have forces on all the response points in order to identify completely t , t , t , so the identi"cation procedure can be applied to P J S a more general experimental modal analysis problem. Another point to be outlined is that the excitation forces need not necessarily be sinusoidal. Both the forces and the responses in equation (7) could be the result of a simple FFT of random signals performed on each single channel of acquisition, without any need to separate the contribution of multiple-input random forces. 2.3.
APPLICATION OF THE EKF TO MODAL PARAMETERS IDENTIFICATION IN THE
FREQUENCY DOMAIN
Equation (7) is non-linear in the unknown parameters. The extended Kalman "lter algorithm can be applied in the frequency domain [11, 13] to the modal parameter identi"cation problem stated in equation (7), with the state vector de"ned as x"+t2,2, t2,2, t2 , m ,2, m ,2, m , u ,2, u ,2, u , t2 , t2,2. S P , P , P , J
(8)
The system equations are x "x I> I z "h (x )#v I I I I
(9)
where z , x and v are referred to the kth forcing frequency u . The function h has the I I I I I same expression as the term on the right-hand side of equation (7) and, therefore, is a non-linear function of the state vector x , the input forces and the frequency u . The I I usual hypotheses of random, white, Gaussian, zero mean and independent noise can be applied as well to the noise on the observations in the frequency domain, while no noise is applied to the state vector, because it is made up by deterministic, constant modal parameters.
332
R. PROVASI E¹ A¸.
The identi"cation technique based on the EKF is a predictor-corrector method of the state estimate and of its corresponding error covariance matrix, working through the following steps: E prediction at the frequency u , using the frequency responses measured at frequencies I> up to u : I xL "xL I> I I I P "P (10) I> I I I E correction using the frequency responses up to u : I> xL "xL #K (z !h (xL )) I> I> I> I I> I> I> I> I P "(I!K H )P (I!K H )2#K R K2 (11) I> I> I> I> I> I I> I> I> I> I> where the "lter gain is given by K "P H2 (H P H2#R )\. (12) I I I\ I I I I\ I I At the end of the iteration cycle, one obtains the estimate xL of state vector and the associated con"dence intervals for each component of xL through its covariance matrix P. 2.4. NUMERICAL PROBLEMS In applying the "ltering procedure to a speci"c system, a dynamic model must be de"ned and noise statistics and initial guesses for the state must be speci"ed. In general, the model is only an approximation of a physical system; the model parameters and the noise statistics may be uncorrect. It has already been outlined that the initial guesses for the unknown parameters can be a rough approximation of the true values and are often speci"ed somewhat arbitrarily when starting the "ltering process. The factors can a!ect the "lter performance and make the process to diverge or, even worse, to converge to unreliable results. The last case can be detected by inspecting the values in the covariance matrix P I I or applying some of the traditional modal model validation techniques. Considering the identi"cation problem stated above, three factors primarily in#uence the convergence of the iterative procedure: E The values of the initial system parameters. E The condition number of the initial error covariance matrix, determined by the relative order of magnitude of the di!erent parameters to be identi"ed. As a matter of fact, the "lter could become insensitive to the parameters with the smallest values. E The condition number of the matrix (HPH2#R), which is updated at each step during the "ltering process. This matrix must be non-negative in order to obtain reliable results from the numerical inversion. As for the initial value of the parameters, it is in general not di$cult to give a rough approximation for the eigenvalues, the error being usually much higher for the damping than for the frequency. It is more di$cult to give a reasonable guess for the order of magnitude of the modes, since they are normalised to the unit mass. Therefore, a devoted algorithm capable of working out automatically from the measured frequency responses a "rst approximation for the mode shapes was developed. To solve the convergence problem, several mathematical techniques have been proposed [14]. In particular, the weighted global iteration technique [15] appears to be more e!ective in reducing the in#uence of the initial values. According to this method, a "rst iteration is carried out with the assigned initial guesses to obtain the converged values after the
EXTENDED KALMAN FILTER
333
processing of all the observations. Any subsequent iteration is performed starting from the "nal values of the previous step, but overweighting the error covariance matrix, until the convergence is obtained or the number of iterations is exhausted. While it was veri"ed that the algorithm can stand a high noise rate on the measured displacements, the noise rejection on the measured forces can be a problem. In its current implementation, the algorithm does not have an explicit formulation for counteracting the noise in the input forces, although it was veri"ed that it can nonetheless stand a signi"cant random noise rate on the forces, as will be discussed later. The practical implementation of the algorithm allows for the choice of a subset of measurement points to be considered in a "rst identi"cation step, i.e. those for which the modes are more clearly excited. The process is then repeated for the remaining points, possibly holding some parameters to the values estimated in the previous step. In particular, but not necessarily, this also applies to the identi"cation of the components of the modes for the points where no forces are acting.
3. EXAMPLE OF APPLICATION
In order to validate the proposed method, a "rst series of tests was performed on a lumped mass}spring system with eight dofs, with viscous damping proportional to the mass and sti!ness matrices so as to obtain real normal modes. The spring connections of each mass with up to "ve other masses were included. A drawing of the system is shown in Fig. 2. This system was set up with the aim of simulating the behaviour of the supporting structure of a complex rotating machine, such as a turbogenerator unit. The resulting mass and sti!ness matrices are reported in Table 1. The damping rate was assigned by proportionality coe$cients of 0.03 and 0.0033 with respect to the mass and the sti!ness values. The reference eigenvalues and eigenvectors were computed by standard routines (Table 2). The synchronous forced response was then obtained over a frequency range 5%65 Hz, thus encompassing all the eigenfrequencies; the frequency resolution was 0.25 Hz. In Fig. 3 it is possible to have an overview of the exciting forces and the relevant frequency responses. The resulting curves are similar to those provided by the monitoring systems during run-ups and shut-downs of turbogenerator units. The exciting forces were essentially a combination of components constant and proportional to the square of the frequency, in order to simulate the e!ect of rotating forces on a real rotor. No particular care was taken in order to obtain the best modal excitation on all modes. The correlation coe$cients between the forces are very close to one.
Figure 2. Drawing of the eight dofs lumped mass}spring system.
334
R. PROVASI E¹ A¸.
TABLE 1 Mass and sti+ness matrices of the eight dofs system Diagonal mass matrix (kg) 2.0
89 400
1.8
!36 000 80 000
Symmetric
1.6
!12 000 !30 000 99 850
1.2
1.7
0.8
Sti!ness matrix (N/m) !3600 !1200 !600 !10 000 !3000 !1000 !39 000 !13 000 !3900 101 000 !33 000 !11 000 107 200 !40 000 92 500
1.2
!1300 !3300 !13 000 !27 000
1.0
!650 !1100 !4000 !9000 30 000 64 750
With these data it was possible to validate the method against problems such as high correlation between the input forces, a comparatively high number of modes to be identi"ed, signi"cant modal coupling due to the damping and presence of local modes, i.e., modes with signi"cant components only on two or three dofs. Moreover, the e!ect of the noise was investigated by addition of a uniformly distributed random error on the measurements and the forces, proportional to the respective maximum amplitudes among all the dofs. Up to eight modes were identi"ed at a time. The MATLAB] version of the code took 2900 s of CPU time with noisy data on a 200 MHz Windows NT PC computer. The number of modes did not seem to be a problem for the algorithm. In some cases, working on a reduced frequency range, modes external to the frequency range of analysis could still be identi"ed with su$cient, sometimes even good, accuracy. The initial values of the eigenfrequencies were taken with a 10% error and the initial modal damping with approximately 100% error (see Table 3). For these tests, the initial values of the mode shapes were generated as a random combination of the true ones, with errors on the single component even greater than 100% and resulting shapes of the modes far away from the correct ones (see Fig. 4). With no noise on the data, almost the exact computational values were obtained for all the parameters. More than acceptable approximations were obtained also with 10% of random error, as de"ned above, on the vibration response. As expected, the method was much more sensitive to the noise on the input forces. At most a 4% error could be added to the forces on this set of data, in order to obtain acceptable results. The results for the identi"ed eigenvalues with a 4% error on both input and output, as in the curves of Fig. 2, are reported in Table 2, where they are compared with the computed values and the values assigned at the "rst iteration step. The only remarkable deviation is on the damping of the "rst eigenfrequency. The amplitudes of the "rst mode components were correspondingly underestimated. As observed from Fig. 3, the signal-to-noise ratio for the forces is very poor at low frequencies, where it is di$cult even to guess a phase value and the energy of the noise is comparable with that of the forces themselves. The identi"ed mode shapes were well correlated with the reference modes obtained by computation, as illustrated in Table 3 by the Modal Assurance Criterion (MAC) coe$cients.
Computed u (Hz) A Identi"ed u (Hz) A Initial u (Hz) A Computed f (%) A Identi"ed f (%) A Initial f (%) A 23.10 23.10 20.79 2.43 2.49 5.40
1.02 1.13 2.27
Mode C 2
9.51 9.51 8.56
Mode C 1
3.53 3.48 7.85
33.64 33.65 30.28
Mode C 3
4.14 4.18 9.20
39.42 39.39 35.48
Mode C 4
4.50 4.39 10.00
42.82 42.82 38.53
Mode C 5
TABLE 2 Computed, identi,ed and initial eigenvalues
5.03 5.00 11.18
47.87 47.87 43.09
Mode C 6
5.51 5.49 12.24
52.40 52.42 47.16
Mode C 7
6.22 6.28 13.82
59.12 59.11 53.21
Mode C 8
EXTENDED KALMAN FILTER
335
336
R. PROVASI E¹ A¸.
Figure 3. Exciting forces and frequency response for the eight dofs system with 4% random error.
In Fig. 4 the whole set of the reference eigenvectors, given by numerical computations, are compared with the initial values assigned at the "rst iteration step and with the identi"ed mode shapes. The overall quality of the results can be evaluated from Fig. 5, where the simulated measurements are compared with the synthesized curves obtained from the identi"ed modal parameters and from their initial values. The curves are shown for three di!erent dofs, as amplitude and phase of the frequency response in a range not encompassing the "rst eigenfrequency. Apart from the "rst mode which was signi"cantly underestimated for the reasons discussed above, for the other modes the curve "tting was very good. As a crosscheck, the curve "tting was then repeated with a set of data generated by forces di!erent from those used in the identi"cation process, with similar results. In order to ascertain the quality of the results, an error index was computed for every single dof in the system as the per cent ratio of the rms values of the error between the curves synthesized by the identi"ed parameters and the rms values of the simulated measurements. A global error index was then obtained as the average error among all the dofs. The same error indices were de"ned with respect to the curves obtained from the initial values assigned to the parameters, in order to appreciate the initial distance from the simulated measurements. The indices are reported in Table 4. Separate values are given in case the frequency range of the "rst mode is taken into account or neglected. Finally, it can be mentioned that the e!ects of having pseudo-resonances in the measurements were also tested. This is the case where there are forces acting at connection points with other components, having their own resonances. Successful results were obtained also with these data sets.
4. CONCLUSIONS
The application of a method based on the same formulation of the extended Kalman "lter and applied in the frequency domain proved to be a suitable tool for solving the problem of the time-invariant modal parameter identi"cation of a mechanical structure excited by sinusoidal multiple inputs having an inherently high degree of correlation.
Mode Mode Mode Mode Mode Mode Mode Mode
C C C C C C C C
Identi"ed
1 2 3 4 5 6 7 8
1.000 0.041 0.001 0.001 0.001 0.000 0.006 0.029
Mode C 1 0.036 1.000 0.035 0.000 0.000 0.001 0.000 0.019
Mode C 2 0.000 0.034 1.000 0.017 0.001 0.006 0.007 0.002
Mode C 3 0.001 0.000 0.014 1.000 0.010 0.007 0.000 0.011
Mode C 4 0.001 0.000 0.001 0.011 1.000 0.013 0.003 0.014
Mode C 5 0.000 0.000 0.006 0.008 0.014 1.000 0.000 0.000
Mode C 6
0.005 0.000 0.007 0.000 0.003 0.000 1.000 0.004
Mode C 7
TABLE 3 Model assurance criterion coe.cients between identi,ed and computed modes
0.028 0.016 0.002 0.011 0.013 0.000 0.003 1.000
Mode C 8
EXTENDED KALMAN FILTER
337
Figure 4. Mode shapes of the eight dofs system: 䊊*䊊*䊊, reference values; # *# *#, identi"ed values; *****, initial values.
338 R. PROVASI E¹ A¸.
339
EXTENDED KALMAN FILTER
Figure 5. Curve "tting with identi"ed and initial modal parameters at dofs C5}7.
TABLE 4 Error values in the curve ,tting with the initial and the identi,ed modal parameters Point 1 Point 2 Point 3 Point 4 Point 5 Point 6 Point 7 Point 8 Mean Initial % rms error} "rst mode included
89.1
84.8
99.8
99.6
92.4
105.6
101.9
88.1
95.2
Final % rms error} "rst mode included
81.1
26.1
29.1
21.9
30.1
18.0
20.7
13.4
22.2
Initial % rms error} "rst mode excluded
88.1
82.6
99.4
97.2
91.2
105.1
101.7
84.9
93.8
Final % rms error} "rst mode excluded
4.7
6.1
6.6
4.5
7.3
3.8
5.5
3.7
5.3
The present paper was aimed at illustrating the formulation of the algorithm and to stress its potential applications, rather than proving its theoretical bases, although it was clearly shown that the proposed formulation in the frequency domain is identical to the more familiar one in the time domain. It was also discussed that the usual hypotheses for the application of the EKF in the time domain are equally satis"ed at least for the considered identi"cation problem dealt with in the frequency domain.
340
R. PROVASI E¹ A¸.
The method was validated so far on di!erent sets of data simulated by a multi-degree-offreedom lumped parameter model, well representative of possible problems to be dealt with in a real situation. These tests, though not exhaustive of any possible occurrence, demonstrated that this approach gave considerably better results than other alternative approaches, in terms of convergency, sensitivity to initial guesses and consistency of the results with respect to di!erent data sets used for the identi"cation. In its current implementation, the method was successfully tested against such problems like strong coupling between di!erent modes, presence of local modes, strong di!erences between mode components measured on di!erent directions, pseudo-resonances in the exciting forces, and poor modal excitation. Up to eight modes were identi"ed at a time. The program was not particularly fast, but the software programming was not yet optimized. The initial values needed to such a recursive algorithm can be worked out from the experimental data ignoring any a priori knowledge, as it would be in a real case. No divergence problems were experienced with such initial guesses, although sometimes they were pretty distant from the true values, but more extensive tests should be performed. In some cases, the ability of the algorithm to identify modes outside the frequency range of analysis was also veri"ed. The random noise rejection was good for noise on the measured vibrations and still acceptable for noise on the input forces, although in the present revision there are no speci"c tools to reduce the e!ect of the last kind of noise. Future developments of the method should include tests against deterministic errors on the input forces, the study of algorithms to "lter the noise on the forces, a generalization of the algorithm for the identi"cation of complex modes, a revision of the iteration scheme that should further increase stability and speed of convergence and the optimisation of the software in order to reduce the computation time. Even in its current revision, the method seems to be appropriate for the application to the modal parameter identi"cation of supporting structure of rotating machines, for which it was primarily developed. Moreover, it could be applied as an additional tool in experimental modal analysis packages, to be used whenever multiple excitations are required in order to introduce more energy into the structure and when sophisticated input controls are not available.
ACKNOWLEDGEMENTS
This work was partially funded by the EU in the project MODIAROT BE95-2015 Contract No. BRPR-CT95-0022 of the research programme BRITE EURAM III.
REFERENCES 1. C. P. FRITZEN 1986 Journal of
EXTENDED KALMAN FILTER
341
6. G. A. ZANETTA 1992 IMechE International Conference on
APPENDIX: NOMENCLATURE M C K F z x I u I w I v I f I g I h I xL I I zL I I u P m P t P u S t J t S F I H I G I K I P I Q I R I
structure mass matrix structure damping matrix structure sti!ness matrix exciting force vector measured vibration vector state vector at time k input vector at time k noise vector on the state at time k noise vector on the observations at time k non-linear function of the state at time k non-linear function on the state noise at time k non-linear function of the observ. at time k state estimate of x based on time k I observation estimate of z based on time k I rth-mode eigenfrequency rth-mode modal damping rth-mode shape vector pseudo-eigenfrequency of the equivalent upper mode shape vector of the lower-mode residuals shape vector of the upper mode residuals matrix of derivatives of f at x I I matrix of derivatives of h at x I I matrix of derivatives of g at x I I "lter gain matrix at time k covariance matrix of the state estim. error at time k covariance matrix of the state noise at time k covariance matrix of the observ. noise at time k