- Email: [email protected]
Identification of Bus Dynamics from Test Data
Copyright
©
IFAC Idclltii'i(.'alioll "lId SySI('1lI Paralllt'tcr
Estimation l'lHC" York , UK, I"H r,
IDENTIFICATION OF BUS DYNAMICS FROM TEST DATA P. Michelberger, A. Keresztes,
J.
Bokor and P. Varlaki
/)I'j}{/rlllll'lll o/IH('('hlllli('s, TI'l'illlim/ Ulli7ll'rsit\' , HIlr/apl'sl, Hllng{l/ )
Abstract. The paper discusses the identification of multivariable models for vehicle structures excited by stochastic road profiles. A time domain approach is used to the identification of the transfer matrix and the modal characteristics /natural frequencies and eigenvalues/ from accelerations and stress measurements. It will be shown that the application of ESS-representation arise naturally in modelling vibrating structures. The methods are uc.ed to identify models of an air suspension bus for various speed levels and road profile variances. Keywords. Commercial vehicle, stiffness, damping, identification, models. INTRODUCTION
modal parameters. In the widely spread frequency domain approach the input force is sinusoidal and the frequency response is measured, or in case of T-transformable inputs,a transfer matrix is computed from the l'- transfonns of input and O:.ltput functions. The modal parameters are usually estimated by a multi-degree freedom curve fit over the frequency range of interest.
The effect of road profile excitation on the design and operation of automotive vehicles has been subject of research activity for a number of years. In the design of commercial road vehicles including buses the determination of dynamic stresses arising in certain structural elements of the body is one of the most important problem. The dynamic analysis requires a mathemathical model of the vehicle and the road profile. The latter one is usually modeled as a stationary stochastic process. The role of identification in vehicle system modeling can be characterised as follows. The body of the vehicle is basically a distributed parameter system excited by the road profile through the wheels, rear and front axles. These last elements can be considered as concentrated parameter systems. The dynamic behaviour of the distributed parameter system can be described by partial differential equations, where the parameters are in general functions of the spatial variable. In vibrating systems these parameters represent mass, stiffness and damping distributions. In the design of vehicles the common practice is to discretize the model e.g. by finite element method, resulting in a second order vector difference equation, thus the requirement to modify the finite element model so as to better correlate with measured data occurs very often. In this case the structure of the model is given apriori and the parameters for subset of parameters/ appearing in the mass and stiffness matrices has to be identified from the measured data. The other situation that turns out very often is the requirement to determine the transfer matrix relating the input excitation to the stress response in given structural elements. Modal parameters play an important role in the dynamic analysis thus the determination of the natural frequencies and eigenvectors of the discrete system can also be the goal of the identification. Various methods are known for identification of the system transfer matrix and
This paper applies the time domain approach to the identification of the transfer ~atrix and the modal parameters from measured stress response and accelerations arising from stochastic road profile excitation. It will be shown, that the application of Elementary Sub-System /ESS/ representation of multivariable systems arise naturally in the identification of vibrating structures. This representation ensures a direct relationship between the transfer matrix and the modal parameters. Using the structure and parameter estimation procedure suggested by Bokor and Keviczky, 1983, not only the transfer matrix but also the number of dominant eigenfrequencies and the elements of eigenvectors can be identified in the time domain. The proposed methods were applied to identify models for an air suspension bus at various speed levels and road profile variances, thus the random vibration analysis of the vehicle can be performed in the whole range of operation. THE EQUATION OF MOTIONS FOR THE VEHICLE The body of the vehicle shown on Fig.l. can be considered as a distributed parameter system excited by the road profile through the wheels, rear and front axles. These last elements can be modelled as concentrated parameter systems. The dynamic behaviour of the body can be described by partial differential equations, but following the common practice, these equations are converted to a set of ordinary differential equations by applymg spatial discretisation and e.g. finite element methods. The bending vibration of the body associated with an orthotrope road nodel can be described by the follOWing equation of the F.Otion:
IH ,I
P. Michelbe rgcr , 1'1 al.
184
M V (t)+ K V(t)+ C V(t)= f g (t); fg(t)= G u(t)
(1)
where V(~~Rn is the vector of vertical displacement of the pOints determined by the method the body is divided into finite elements, and of the vertical axle displacements. The matrices M, K, C are the mass, damping and stiffness matrices respectively, and their structure is determined by the finite-element method applied, but generally these are symmetric matrices and M is positive definite. The elements of u(t) are the road profile displacements under the front and rear tires, usually it is assumed that u 2 (t)= ul(t-T), and T is defined by the speed and axle distance of the vechile. The elements of G are the tire-stiffnesses. In many cases the dynamic stresses arising in certain elements of the body are measured. With spatial discretization the form of equations for motions is similiar to that in Eq. 1, but applying finite element method with special appoximating functions, the first n-l elements of vet) can be made proportional to the second derivatives of the displacements according to the spatial variable, and the last two are the axle displacements. For the construction of the appropriate function space of global approximating polynomials and the structure of M, K, C matrices see Michelberger et-al, 1984. The direct relation between the derivatives in vet) and the measured stress vector m y (t) E R is given by yet) = H V(t) (2) Where the matrix H is of full rank m, and its rows are usually unit vectors in Rn. The direct relationship between v(t) and fg(t) is given by the transfer matrix T(s) obtained from the Laplace-transform of Eq.l with zero initial conditions: T(S)
(M
=
D-1
lM
s2+
s +
K
c]
D(S) = det[M s2 +
K
s +
c] ;
1t
i=l
nc (s+A.) l.
'Jt'
j =1
(s+).
+ .)(s+).,x + .) nr J nr J
where
(5)
n , nc are the number of real and r complex pair of poles in T(s) respectively. In Eq.4 it was assumed that all the roots of D(s) have mUltiplicity one, which is the case very often in the practice, but all the foregoing discussion can be extended for multiple poles too. It is known, that the second order differential equation in Eq.l can be rewritten as Me ~ (t) + C e z (t) = f e (t ) where ZT(t)=[VT(t), vTlt)}, z, feER2n and
f~(t)=[f~(t),
C e = diag { C, - MS'
(6 )
u
+ ' ... ; nr l
u~
r
+n
u~r+l"'"
J.
c
(7)
Let the matrix X be constructed from U as (8) XT = [ UT: /\ UT Utilysing the simmetric structure of matrices in Eq.l and Eq.6., one can conlude, that the following normalising conditions hold:
1.
(9)
XT M X = I e 2n and the extended transfer matrix associated wi~h Eq.6 is given by
Te(S)
2n
Te (s)=
L
(10)
i=l
As the eigenvectors u consist of the i first n elements of x.l. the transfer matrix T{s) nxn minor of 2n T
L i=l
(s)
~
can be obtained as the first Tees) , as T nr T Uj u i ui ui + s +).i i=l s+ ;A. i
n ( _u_..::r:...+_J_·
deg D(S)= 2n (4)
Let the polynomial D(s) be factored over the field of complex numbers as nr
U = [ u l ' ... , u nr
j=nr+l
The link between the transfer matrix and the modal characteristics /natural frequencies, eigenvectors/ can be shown as follows.
D(S) =
)., ~n ] , the diagonal matrix of the eigenvalues, and let the matrix of the associated eigenvectors be denoted by U
+
where
n
equations of Eq. 6. Denote /\.= diag[ \ , '"
=
s2 + K s + c)-l
(s ) ad j
Eq.l can be reobtained from the first
2:'
_u_~.!:r_+_j
. (ll)
+
s +An +j r
Eq.ll shows, how to obtaine the transfer matrix T(s) if the modal parameters /eigenvectors, eigenvalues/ are given. In case of T{s) is given e.g. as a result of identification, the following procedure can be applied. The following result is known from linear algebra. Let in the nonsingular matrix pencil s Me + Ce be Me nonsingular, and assume that the polynomial + C
/
De(S)= det/sM e +
has roots with multiplicity one.
e Then utilysing the symmetric structure of Me' C , it can be proved that e adj[s Me+ Cel/D'(s)!s=
A.
= xixr ;
l.
i =1, ... ,2n
(12)
where D'(~= d D(~ /d s , x. is the eigenvector associated with tfie eigenvalue ~ , and xi satisfies the normalising conditions in Eq.9. The eigenvectors are obtained from the first of xi'
n
ui elements
185
Identification of Bus Dynamics
As adj[s Me + C~1/D'(s)IS-.Ai are the residue matrices at s = Ai' the following procedure can be used to obtain the modal parameters. Assume, that T(S) is given.
of T(s). In the identification of vibrating structures, however, this model has the disadvantage that the direct relationship with the modal characteristics is lost, and the estimation of the degree of polynomials can be quite involved.
(i) Construct the least common denominator D (s) of the elements of T (s), and factor it as shown in Eq.5. The roots of DeS) are the eigenvalues Ai ' i=l, ... , 2n of T(s).
As it was shown in the previous paragraph, the structure of the transfer matrix T(s) in Eq.ll is given by the sum of partial fraction matrices, where the numerator matrices Pi are given as the dyadic pro-
(ii)Decompose each elements of T(S) into patial fractions, and collect each terms associated with Ai into the partial
ducts of the eigenvectors. Thus the application of the Elementary SubSystem IESSI representation of multivarialbe system is straightforward for the identification of vibrating systems. The ESS representation as defined in Bokor and Keviczky, 1983, is based on the decomposition of the discrete transfer matrix into partial fraction matrices. The partial fraction matrices are further decomposed, using a special dyadic decomposition procedure, given in Bokor and Peto, 1981. The number of dyads is equal to the rank of the numerator matrix associated with a given pole. The multivariable ESS-s associated with a given pole Ai can be represented by a transfer
fraction matrix
Ti(S), where
Ti(S)= Pi/(s + Ai)'
rank Pi = 1, i = 1, ... , 2n
(iii) Construct the Iminimall position of Pi as
(13)
dyadic decom-
i=l, ... , 2n , to obtain the eigenvectors
u . i
This procedure can be generalised for multiple poles, but in this case rank p, ~ 1 and in step (iii) a minimal dyadi~ decomposition procedure has to be used. In case of stress ~easurements, the transfer matrix relating the road profile excitation u(t) to the measured stress vector can be written from Eq.1-2 as 2n
L i=l
It can be seen, that either the second order derivatives in vet) or the stresses in yet) = H v(s) are measured, the structure of the transfer matrix arising in the dynamic modelling of vehicles is given by the sum of partial fraction matrices associated with the poles li.e. with natural frequencies I of the system. This stresses that this structure should also be used for identification as it will be done in the next paragraph. IDENTIFICATION OF THE TRANSFER MATRIX AND MODAL CHARACTERISTICS Various methods are known for identification of the system transfer matrix. In the frequecy domain approach the input force is sinusoidal and the frequency response is measured, or in case of the input is Fourier transformable function, the elements of the transfer matrix are computed from the Fourier-transforms of input and output functions. Using the frequency response or transfer functions, the modal parameters are usually estimated by multidegree of freedom curve fit over the frequency range of interest. The transfer matrix can also be identified in the domain, where the application of well developed and powerful discrete time structure and parameter estimation procedures become possible. The discrete-time equivalent of T(s) is also a rational matrix thus the identification of T(s) means the estimation of the degrees and coefficients in the numerator and denominator polynomi.a ls, and possible the dead times in the elements
matrix with a single pole
Ai
and numera-
tor matrix consisting of the dyadic products of the vectors obtained in the above dyadic decomposition. One can conclude, that the form of representing T(s) by the eigenvectors in Eq. 11 is just a special case of ESS representation, thus the ESS model structures naturally arise in modelling vibrating sys~ tern. Applying ESS-representations in the identification procedure, not only the estimation of the transfer matrix, but also the determination of the number of natural frequencies Ipoles I and the estimation of the elements of the associated eigenvectors can be solved, as it will be shown in the sequel. In the discrete-time identification, the discrete time equivalent of T(s) is needed, which can be obtained certain assumtions on the approximations of input process during the sampling time interval. An example is the zero order approximation, u(kh+t) =u(kh} if o=t=h where h is the sampling interval. Denote tij(s) the i,j - th element in T (s) : ., n
p,lJ
~ t-'k (s)= L -- + 1.J k=l s +~
t"
+
i ( f>~~+k k=l
s +
A. n
+ +k r
(15)
Cl.kij = u t-' u ' or t a k'1.ng th e oVserik jk vation in Eq.2 into consideration.
where
p,ij = hT, T 1. u k u k gj
I'"' k
(16)
It can be shown, that using zero order approxima tion for u (t), the i, j - th element in T (z) , the discrete time equivalent of vT(S) is given by
P. Michelberger, 1'1 al.
186
where
Yt
denote the vector of measured,
sampled noisy stress responses, TN(z),
+
(
bij
TNri(z), TNcj(z) denote the noise transfer matrix and its partial fraction matrices, and et denotes anm-dimensional
') x
n r +k ) --=....:..:x~
+
(17)
Z + an +k r
zut=u t + l ' denotes the sampled input process, and
where Ut
Z
is the shift operator,
zero-mean white noise vector sequence with diagonal covariance matrix. The properties of multivariable stochastic models matrix have been discussed in Bokor and Keviczky, 1983, where also an effective structure arid parameter estimation procedure has been developed. It is shown that the parameter vector T ' e T = [T Qv, eN where @v' @N consist of the unknown coefficients in the partial fraction matrices of Tv(Z) and TNlz) res-
1
i
In order to avoid the estimation of complex coefficients, and to have Tv(Z) in a form directly related to Eq.ll, we can write: n Tv(Z)
n
r
L
A.
--~- +
i=l
Z + a.
~
c
BljZ + B2j
L
2
j=l Z + P ij z+P2j
n
n
t
pectively, can be estimated by the method of maximum likelihood /ML/. The ML estimation of e has the following properties.
(22a) c
L
Tvr~. (z) +
i=l
(l9)
TVC]. (z)
and let the estimated value of
j=l
noted by where partial fraction matrices
I3
9
k
V
be de-
v '
Tvrilz),
T . (z) are the transfer matrices of ESS-s aK~6ciated with i-th real and j-th complex poles respectively, and z2+PljZ+P2j =(z+
where k denotes the number of samples used in the estimation and ~ n , ~
)"n+j)(Z+A~+j)' r
nr
r
j=l, .•. ,n
c
'
~
BljZ + B 2j
(20)
constructed analogously to Eq.22a and 22b. It has been prooved in Keviczky et-al., 1984, that under certain conditions the following limit relations hold w.p.l. as
mod R[z] /(z +Pljz+P2j) if
f
v(t) and T
H bibiG ;
g
k-.oo
[g l z+g2] [g l z+g2]T
2
nc -
indicating that the number of poles i~ overestimated. The vectors eN' 9~ are
The numerator matrices A., B ·, B . are l 2 obtained from Eq.17 using~step~(i) 1 liii) in the procedure described in the previous paragraph, and have the properties, that Ai = bibi
r
n
i=l, ... ,n r
(21a)
'"Ap( i)·
~=nr +
1
°
~
, ... , nr
(t) are measured, and
Bljz+B2j~H
mod R[Z]
[
T
j=l, ... , nc
g l z+g2] [g l z+g2] G
2 / (z +Plj+P2j)
(21b)
Hv(t) and ult) are measured ( R(Z) / Y (z) denotes the residue class ring according to the polynomial Y ~z) ) . The measurements are usually carried out in noisy environment because of neglected nonsignificant dynamic modes and the measurement noise. This situation can be modelled by the inclusion of an output noise model into the system model if
mr
Tv (z) Ut + TN(Z) et' TN (Z) =
I
i=l
T Nri (z) +
"Bl'O(j)O;
"'-
B 2 'OUro
j=nc+l,·
"'~c'
where p(i), O(j) denote permutations of the first n a n d n integers, respectively. Similaf limit rglation hold for the k
parameters in the noise model. Eq.23 N intuitively states that in the model the numerator coefficients in the partial fraction matrices associated with a pole / eigenfrequency/ that do not occur in the system converge to zero. These properties utilised in the structure and parameter e s timation procedures can be used to determine the number and estimates of dominant poles /natural frequencies/ in the frequency range of interest. As the numerator matrices are parametrized e.g. as A. b.b T for real Q
~
~
~
187
Identification of Bus Dynamics
poles and the vectors
b
are directly esi timated, the previous results show, that using the ESS-representations, the identification of the eigenvectors or some of its eleme nts can also be solved in the time domain. EXPERIMENTAL CONDITIONS AND IDENTIFICATION RESULTS Experiments were carried out on airsuspension bus, using a computer-controled Shenk Hydropulse /SH/ apparatus. The input excitation of the wheels were generated by the vertical displacement of the cylinders of the SH apparatus. The cylinders were forced to moove by computer control in accordance with the previously registered stochastic road profile realisation. The left and right wheels had the same excitation /orthotrope road/, and the front wheel excitation ul(t) was simple delayed to get the rear wheel exitation u 2 (t)= Ul(t-T s )' where Ts was computed from axle distance and the speed of the vehicle. The variance of the road profile process and the speed of the vehicle could be made to vary continuously between given lower and upper bounds, The possible location of stress gauges are shown on Fig.l. Two accelerometers were also used to measure the front and rear axle accelerations. The input and the measured output processes were continuously registred on magnetic tape and sampled for use as input data to the digital computer programs de v eloped to perform the identification procedure. The sampling inter v al h=0,02 sec. proved to be satisfactory. The numb e r of samples used in the identification was N=500. Measurements were obtained at speed levels 20,30, ... ,80 km/h, and five road profile variances. The usual statistical functions like auto- and crosscorrelation functions, auto- and cross spectral densities e.t.c. were first computed and used in the preliminary data analysis. In the for e going discussion, the model identified from two stress processes Yl (t), Y2(t), and from the front and rear axle displace~ ments Y3(t), Y4(t) as outputs and from th e road profile displacement under the front and real tires as inputs will be investigated. Because of limited space, only the mod e l idcp.tified at speed 60 km/h and road profilevariance 1 cm will be discussed. Omitting the details of structure and parame t e r estimation procedure, the final estimate of the vehicle and noise transfer matric e s are the following:
D. l ...
t
O.ot
.2~O .h.U.)'
~L .I
-O.1 0 1." O. ~.
+o.l 2;1 • 0. 1 j
.2-o . 4 7 .t-O, h
.,
o. lla - O .• } I. -o.09J.0.41
-O . ilK .. 0. 10
.. 2_ 0 ,o'l+o , n
As the elements of T (~ are directly identified in the form de~cribed in Eq. 19, the parameters of the continuous time model can be computed from Eqs. 15, 17, 18. When the modal characteristics associated with a finite element model are of interest, Eq.14 has to be taken into consideration also. The partial fractions in the elements of the transfer matrix Tv(Z) can be related to the natural frequencies of the vehicle. Denote f . a natural frequency as~ociated with the denominator polynomial Z +Pijz+ + Pzj· Then a direct relation can be given to compute fj from the discrete time model parameters,
Computing the frequencies f. associated with the identified transferJmatrix, the model can be explained as follows. The estimated frequencies can be clustered into four groups. The values around 11.5-12 Hz can be related to the second bending frequency of the body, and this. frequency appears in each outputs. The frequencies around 10.3-10.4 Hz appear in stress outputs Y2t and in axle accelerations Y3t' Y4t. These can be related to the natural frequencies of the front and rear axles. The frequency 9.1 Hz occurs only in the stress output Ylt' and it approximates the naturalfrequency of the engine. It c an be observed, that the input excites the system in the frequency range 9 Hz: f~ -12.5 Hz. The frequencies associated with the noise model poles lie in the range 3.8 Hz~ f ~ 10.5 Hz. Thus it can be expecN ted, that the model of vehicle given by the transfer matrix T (z), gives a good approximation in the fr~quency range 9Hz~ f~ £ 12.5 Hz, and for lower frequencies the effect of noise is dominant. This noise effect can be rather explained by the neglec ted nonsignificant low-frequency dynamic modes then by measurement errors. To analyse the above phenomena the logarithmic amplitude, 20 log I s .. (f)1 and the 1.J
amplitude-phase diagrams
Sij (jf)
were
computed by the use of continuous transfer matrix T (s) and were also estimated by cross-spe~tral method using 2000 samples. For the illustration of the results these diagrams for a stress output Y2t and for an acceleration output Y4t are shown on Fig.3. The comparison of ehe computed and estimated diagrams supports the above explanations. In order to illustrate the fit of the rrodel in the time domain, 100 samples of the measured Y2t' Y4t processes and of the corresponding outputs computed by the identified model are plotted on Fig.2.
P.
188
Mi c helbc r~e r , 1'/
ni.
CONCLUSION The paper discussed the identification of transfer matrix model for air suspension bus structures. The identified models can be used to determine the dynamic stress e s and stress st~tistics for arbitrary road profile excitation and speed levels within the examined range. REFERENCES 1. Michelberger,P., Bokor, J.,Keresztes A. and Varlaki P. (198~. Multivariable identification for commercial vehicle stress analysis. American Control Conference, San Francisco, California,
Fiq.Z-Plots of measurements and of identified model outputs. model ---- measurements
outputs
2. Ross, G.R. (1971). Synthesis of stiffness and mass matrices from experimental vibration modes. SAE paper. No 710787. 3. Natke, H.G. (1983). Einfilhrung in Theorie und Praxis der zeitreichen - und Modalanalyse. Braunschweig ; Wiesbaden: Vieweg. 4. Bokor,J. and Keviczky,L. \1984). Structure and parameter estimation of MIMO system using elementary sub-system representation. Int. Journal of Control., Vol 39., No.S. 5. Keviczky,L. Bokor,J. and Veres S. \1984). Strong consistency of ML estimators using partial fraction and elementary subsystem representation of multivariable system. 9 th • IFAC World Congress, Budapest 6. Bokor,J. and Pete, J. (1981).Minimal realisation procedure from partial fraction representation of a rational transfer matrix. Problems of Control and Int. Theory, Vol. 10, 7. Meirovitch, A. and Baruh,H. (1982). Identification of distributed-parameter vibrating systems. Proc . 6. IFAC Symp. on Ident, and System Par.Est., Washington 8. Michelberger,P. Keresztes, A. Bokor,J. and Varlaki, P. (1984) Dynamic Modelling of Commercial Road Vehicle Structures from test Data. XX FISITA Congress, Vol. 6-12. pp.4.96-4 . 103 . Wien.
Y, Y2 stresses
Y3 accelerations Y4 Fiq.l .- Locations of stress gauges acceleraneters.
am
Pig.) Oro . . epeotra - - - oOlDJlUte4 117 tlw data - - . .t1mate4 n- tile dat.
©
IFAC Idclltii'i(.'alioll "lId SySI('1lI Paralllt'tcr
Estimation l'lHC" York , UK, I"H r,
IDENTIFICATION OF BUS DYNAMICS FROM TEST DATA P. Michelberger, A. Keresztes,
J.
Bokor and P. Varlaki
/)I'j}{/rlllll'lll o/IH('('hlllli('s, TI'l'illlim/ Ulli7ll'rsit\' , HIlr/apl'sl, Hllng{l/ )
Abstract. The paper discusses the identification of multivariable models for vehicle structures excited by stochastic road profiles. A time domain approach is used to the identification of the transfer matrix and the modal characteristics /natural frequencies and eigenvalues/ from accelerations and stress measurements. It will be shown that the application of ESS-representation arise naturally in modelling vibrating structures. The methods are uc.ed to identify models of an air suspension bus for various speed levels and road profile variances. Keywords. Commercial vehicle, stiffness, damping, identification, models. INTRODUCTION
modal parameters. In the widely spread frequency domain approach the input force is sinusoidal and the frequency response is measured, or in case of T-transformable inputs,a transfer matrix is computed from the l'- transfonns of input and O:.ltput functions. The modal parameters are usually estimated by a multi-degree freedom curve fit over the frequency range of interest.
The effect of road profile excitation on the design and operation of automotive vehicles has been subject of research activity for a number of years. In the design of commercial road vehicles including buses the determination of dynamic stresses arising in certain structural elements of the body is one of the most important problem. The dynamic analysis requires a mathemathical model of the vehicle and the road profile. The latter one is usually modeled as a stationary stochastic process. The role of identification in vehicle system modeling can be characterised as follows. The body of the vehicle is basically a distributed parameter system excited by the road profile through the wheels, rear and front axles. These last elements can be considered as concentrated parameter systems. The dynamic behaviour of the distributed parameter system can be described by partial differential equations, where the parameters are in general functions of the spatial variable. In vibrating systems these parameters represent mass, stiffness and damping distributions. In the design of vehicles the common practice is to discretize the model e.g. by finite element method, resulting in a second order vector difference equation, thus the requirement to modify the finite element model so as to better correlate with measured data occurs very often. In this case the structure of the model is given apriori and the parameters for subset of parameters/ appearing in the mass and stiffness matrices has to be identified from the measured data. The other situation that turns out very often is the requirement to determine the transfer matrix relating the input excitation to the stress response in given structural elements. Modal parameters play an important role in the dynamic analysis thus the determination of the natural frequencies and eigenvectors of the discrete system can also be the goal of the identification. Various methods are known for identification of the system transfer matrix and
This paper applies the time domain approach to the identification of the transfer ~atrix and the modal parameters from measured stress response and accelerations arising from stochastic road profile excitation. It will be shown, that the application of Elementary Sub-System /ESS/ representation of multivariable systems arise naturally in the identification of vibrating structures. This representation ensures a direct relationship between the transfer matrix and the modal parameters. Using the structure and parameter estimation procedure suggested by Bokor and Keviczky, 1983, not only the transfer matrix but also the number of dominant eigenfrequencies and the elements of eigenvectors can be identified in the time domain. The proposed methods were applied to identify models for an air suspension bus at various speed levels and road profile variances, thus the random vibration analysis of the vehicle can be performed in the whole range of operation. THE EQUATION OF MOTIONS FOR THE VEHICLE The body of the vehicle shown on Fig.l. can be considered as a distributed parameter system excited by the road profile through the wheels, rear and front axles. These last elements can be modelled as concentrated parameter systems. The dynamic behaviour of the body can be described by partial differential equations, but following the common practice, these equations are converted to a set of ordinary differential equations by applymg spatial discretisation and e.g. finite element methods. The bending vibration of the body associated with an orthotrope road nodel can be described by the follOWing equation of the F.Otion:
IH ,I
P. Michelbe rgcr , 1'1 al.
184
M V (t)+ K V(t)+ C V(t)= f g (t); fg(t)= G u(t)
(1)
where V(~~Rn is the vector of vertical displacement of the pOints determined by the method the body is divided into finite elements, and of the vertical axle displacements. The matrices M, K, C are the mass, damping and stiffness matrices respectively, and their structure is determined by the finite-element method applied, but generally these are symmetric matrices and M is positive definite. The elements of u(t) are the road profile displacements under the front and rear tires, usually it is assumed that u 2 (t)= ul(t-T), and T is defined by the speed and axle distance of the vechile. The elements of G are the tire-stiffnesses. In many cases the dynamic stresses arising in certain elements of the body are measured. With spatial discretization the form of equations for motions is similiar to that in Eq. 1, but applying finite element method with special appoximating functions, the first n-l elements of vet) can be made proportional to the second derivatives of the displacements according to the spatial variable, and the last two are the axle displacements. For the construction of the appropriate function space of global approximating polynomials and the structure of M, K, C matrices see Michelberger et-al, 1984. The direct relation between the derivatives in vet) and the measured stress vector m y (t) E R is given by yet) = H V(t) (2) Where the matrix H is of full rank m, and its rows are usually unit vectors in Rn. The direct relationship between v(t) and fg(t) is given by the transfer matrix T(s) obtained from the Laplace-transform of Eq.l with zero initial conditions: T(S)
(M
=
D-1
lM
s2+
s +
K
c]
D(S) = det[M s2 +
K
s +
c] ;
1t
i=l
nc (s+A.) l.
'Jt'
j =1
(s+).
+ .)(s+).,x + .) nr J nr J
where
(5)
n , nc are the number of real and r complex pair of poles in T(s) respectively. In Eq.4 it was assumed that all the roots of D(s) have mUltiplicity one, which is the case very often in the practice, but all the foregoing discussion can be extended for multiple poles too. It is known, that the second order differential equation in Eq.l can be rewritten as Me ~ (t) + C e z (t) = f e (t ) where ZT(t)=[VT(t), vTlt)}, z, feER2n and
f~(t)=[f~(t),
C e = diag { C, - MS'
(6 )
u
+ ' ... ; nr l
u~
r
+n
u~r+l"'"
J.
c
(7)
Let the matrix X be constructed from U as (8) XT = [ UT: /\ UT Utilysing the simmetric structure of matrices in Eq.l and Eq.6., one can conlude, that the following normalising conditions hold:
1.
(9)
XT M X = I e 2n and the extended transfer matrix associated wi~h Eq.6 is given by
Te(S)
2n
Te (s)=
L
(10)
i=l
As the eigenvectors u consist of the i first n elements of x.l. the transfer matrix T{s) nxn minor of 2n T
L i=l
(s)
~
can be obtained as the first Tees) , as T nr T Uj u i ui ui + s +).i i=l s+ ;A. i
n ( _u_..::r:...+_J_·
deg D(S)= 2n (4)
Let the polynomial D(s) be factored over the field of complex numbers as nr
U = [ u l ' ... , u nr
j=nr+l
The link between the transfer matrix and the modal characteristics /natural frequencies, eigenvectors/ can be shown as follows.
D(S) =
)., ~n ] , the diagonal matrix of the eigenvalues, and let the matrix of the associated eigenvectors be denoted by U
+
where
n
equations of Eq. 6. Denote /\.= diag[ \ , '"
=
s2 + K s + c)-l
(s ) ad j
Eq.l can be reobtained from the first
2:'
_u_~.!:r_+_j
. (ll)
+
s +An +j r
Eq.ll shows, how to obtaine the transfer matrix T(s) if the modal parameters /eigenvectors, eigenvalues/ are given. In case of T{s) is given e.g. as a result of identification, the following procedure can be applied. The following result is known from linear algebra. Let in the nonsingular matrix pencil s Me + Ce be Me nonsingular, and assume that the polynomial + C
/
De(S)= det/sM e +
has roots with multiplicity one.
e Then utilysing the symmetric structure of Me' C , it can be proved that e adj[s Me+ Cel/D'(s)!s=
A.
= xixr ;
l.
i =1, ... ,2n
(12)
where D'(~= d D(~ /d s , x. is the eigenvector associated with tfie eigenvalue ~ , and xi satisfies the normalising conditions in Eq.9. The eigenvectors are obtained from the first of xi'
n
ui elements
185
Identification of Bus Dynamics
As adj[s Me + C~1/D'(s)IS-.Ai are the residue matrices at s = Ai' the following procedure can be used to obtain the modal parameters. Assume, that T(S) is given.
of T(s). In the identification of vibrating structures, however, this model has the disadvantage that the direct relationship with the modal characteristics is lost, and the estimation of the degree of polynomials can be quite involved.
(i) Construct the least common denominator D (s) of the elements of T (s), and factor it as shown in Eq.5. The roots of DeS) are the eigenvalues Ai ' i=l, ... , 2n of T(s).
As it was shown in the previous paragraph, the structure of the transfer matrix T(s) in Eq.ll is given by the sum of partial fraction matrices, where the numerator matrices Pi are given as the dyadic pro-
(ii)Decompose each elements of T(S) into patial fractions, and collect each terms associated with Ai into the partial
ducts of the eigenvectors. Thus the application of the Elementary SubSystem IESSI representation of multivarialbe system is straightforward for the identification of vibrating systems. The ESS representation as defined in Bokor and Keviczky, 1983, is based on the decomposition of the discrete transfer matrix into partial fraction matrices. The partial fraction matrices are further decomposed, using a special dyadic decomposition procedure, given in Bokor and Peto, 1981. The number of dyads is equal to the rank of the numerator matrix associated with a given pole. The multivariable ESS-s associated with a given pole Ai can be represented by a transfer
fraction matrix
Ti(S), where
Ti(S)= Pi/(s + Ai)'
rank Pi = 1, i = 1, ... , 2n
(iii) Construct the Iminimall position of Pi as
(13)
dyadic decom-
i=l, ... , 2n , to obtain the eigenvectors
u . i
This procedure can be generalised for multiple poles, but in this case rank p, ~ 1 and in step (iii) a minimal dyadi~ decomposition procedure has to be used. In case of stress ~easurements, the transfer matrix relating the road profile excitation u(t) to the measured stress vector can be written from Eq.1-2 as 2n
L i=l
It can be seen, that either the second order derivatives in vet) or the stresses in yet) = H v(s) are measured, the structure of the transfer matrix arising in the dynamic modelling of vehicles is given by the sum of partial fraction matrices associated with the poles li.e. with natural frequencies I of the system. This stresses that this structure should also be used for identification as it will be done in the next paragraph. IDENTIFICATION OF THE TRANSFER MATRIX AND MODAL CHARACTERISTICS Various methods are known for identification of the system transfer matrix. In the frequecy domain approach the input force is sinusoidal and the frequency response is measured, or in case of the input is Fourier transformable function, the elements of the transfer matrix are computed from the Fourier-transforms of input and output functions. Using the frequency response or transfer functions, the modal parameters are usually estimated by multidegree of freedom curve fit over the frequency range of interest. The transfer matrix can also be identified in the domain, where the application of well developed and powerful discrete time structure and parameter estimation procedures become possible. The discrete-time equivalent of T(s) is also a rational matrix thus the identification of T(s) means the estimation of the degrees and coefficients in the numerator and denominator polynomi.a ls, and possible the dead times in the elements
matrix with a single pole
Ai
and numera-
tor matrix consisting of the dyadic products of the vectors obtained in the above dyadic decomposition. One can conclude, that the form of representing T(s) by the eigenvectors in Eq. 11 is just a special case of ESS representation, thus the ESS model structures naturally arise in modelling vibrating sys~ tern. Applying ESS-representations in the identification procedure, not only the estimation of the transfer matrix, but also the determination of the number of natural frequencies Ipoles I and the estimation of the elements of the associated eigenvectors can be solved, as it will be shown in the sequel. In the discrete-time identification, the discrete time equivalent of T(s) is needed, which can be obtained certain assumtions on the approximations of input process during the sampling time interval. An example is the zero order approximation, u(kh+t) =u(kh} if o=t=h where h is the sampling interval. Denote tij(s) the i,j - th element in T (s) : ., n
p,lJ
~ t-'k (s)= L -- + 1.J k=l s +~
t"
+
i ( f>~~+k k=l
s +
A. n
+ +k r
(15)
Cl.kij = u t-' u ' or t a k'1.ng th e oVserik jk vation in Eq.2 into consideration.
where
p,ij = hT, T 1. u k u k gj
I'"' k
(16)
It can be shown, that using zero order approxima tion for u (t), the i, j - th element in T (z) , the discrete time equivalent of vT(S) is given by
P. Michelberger, 1'1 al.
186
where
Yt
denote the vector of measured,
sampled noisy stress responses, TN(z),
+
(
bij
TNri(z), TNcj(z) denote the noise transfer matrix and its partial fraction matrices, and et denotes anm-dimensional
') x
n r +k ) --=....:..:x~
+
(17)
Z + an +k r
zut=u t + l ' denotes the sampled input process, and
where Ut
Z
is the shift operator,
zero-mean white noise vector sequence with diagonal covariance matrix. The properties of multivariable stochastic models matrix have been discussed in Bokor and Keviczky, 1983, where also an effective structure arid parameter estimation procedure has been developed. It is shown that the parameter vector T ' e T = [T Qv, eN where @v' @N consist of the unknown coefficients in the partial fraction matrices of Tv(Z) and TNlz) res-
1
i
In order to avoid the estimation of complex coefficients, and to have Tv(Z) in a form directly related to Eq.ll, we can write: n Tv(Z)
n
r
L
A.
--~- +
i=l
Z + a.
~
c
BljZ + B2j
L
2
j=l Z + P ij z+P2j
n
n
t
pectively, can be estimated by the method of maximum likelihood /ML/. The ML estimation of e has the following properties.
(22a) c
L
Tvr~. (z) +
i=l
(l9)
TVC]. (z)
and let the estimated value of
j=l
noted by where partial fraction matrices
I3
9
k
V
be de-
v '
Tvrilz),
T . (z) are the transfer matrices of ESS-s aK~6ciated with i-th real and j-th complex poles respectively, and z2+PljZ+P2j =(z+
where k denotes the number of samples used in the estimation and ~ n , ~
)"n+j)(Z+A~+j)' r
nr
r
j=l, .•. ,n
c
'
~
BljZ + B 2j
(20)
constructed analogously to Eq.22a and 22b. It has been prooved in Keviczky et-al., 1984, that under certain conditions the following limit relations hold w.p.l. as
mod R[z] /(z +Pljz+P2j) if
f
v(t) and T
H bibiG ;
g
k-.oo
[g l z+g2] [g l z+g2]T
2
nc -
indicating that the number of poles i~ overestimated. The vectors eN' 9~ are
The numerator matrices A., B ·, B . are l 2 obtained from Eq.17 using~step~(i) 1 liii) in the procedure described in the previous paragraph, and have the properties, that Ai = bibi
r
n
i=l, ... ,n r
(21a)
'"Ap( i)·
~=nr +
1
°
~
, ... , nr
(t) are measured, and
Bljz+B2j~H
mod R[Z]
[
T
j=l, ... , nc
g l z+g2] [g l z+g2] G
2 / (z +Plj+P2j)
(21b)
Hv(t) and ult) are measured ( R(Z) / Y (z) denotes the residue class ring according to the polynomial Y ~z) ) . The measurements are usually carried out in noisy environment because of neglected nonsignificant dynamic modes and the measurement noise. This situation can be modelled by the inclusion of an output noise model into the system model if
mr
Tv (z) Ut + TN(Z) et' TN (Z) =
I
i=l
T Nri (z) +
"Bl'O(j)O;
"'-
B 2 'OUro
j=nc+l,·
"'~c'
where p(i), O(j) denote permutations of the first n a n d n integers, respectively. Similaf limit rglation hold for the k
parameters in the noise model. Eq.23 N intuitively states that in the model the numerator coefficients in the partial fraction matrices associated with a pole / eigenfrequency/ that do not occur in the system converge to zero. These properties utilised in the structure and parameter e s timation procedures can be used to determine the number and estimates of dominant poles /natural frequencies/ in the frequency range of interest. As the numerator matrices are parametrized e.g. as A. b.b T for real Q
~
~
~
187
Identification of Bus Dynamics
poles and the vectors
b
are directly esi timated, the previous results show, that using the ESS-representations, the identification of the eigenvectors or some of its eleme nts can also be solved in the time domain. EXPERIMENTAL CONDITIONS AND IDENTIFICATION RESULTS Experiments were carried out on airsuspension bus, using a computer-controled Shenk Hydropulse /SH/ apparatus. The input excitation of the wheels were generated by the vertical displacement of the cylinders of the SH apparatus. The cylinders were forced to moove by computer control in accordance with the previously registered stochastic road profile realisation. The left and right wheels had the same excitation /orthotrope road/, and the front wheel excitation ul(t) was simple delayed to get the rear wheel exitation u 2 (t)= Ul(t-T s )' where Ts was computed from axle distance and the speed of the vehicle. The variance of the road profile process and the speed of the vehicle could be made to vary continuously between given lower and upper bounds, The possible location of stress gauges are shown on Fig.l. Two accelerometers were also used to measure the front and rear axle accelerations. The input and the measured output processes were continuously registred on magnetic tape and sampled for use as input data to the digital computer programs de v eloped to perform the identification procedure. The sampling inter v al h=0,02 sec. proved to be satisfactory. The numb e r of samples used in the identification was N=500. Measurements were obtained at speed levels 20,30, ... ,80 km/h, and five road profile variances. The usual statistical functions like auto- and crosscorrelation functions, auto- and cross spectral densities e.t.c. were first computed and used in the preliminary data analysis. In the for e going discussion, the model identified from two stress processes Yl (t), Y2(t), and from the front and rear axle displace~ ments Y3(t), Y4(t) as outputs and from th e road profile displacement under the front and real tires as inputs will be investigated. Because of limited space, only the mod e l idcp.tified at speed 60 km/h and road profilevariance 1 cm will be discussed. Omitting the details of structure and parame t e r estimation procedure, the final estimate of the vehicle and noise transfer matric e s are the following:
D. l ...
t
O.ot
.2~O .h.U.)'
~L .I
-O.1 0 1." O. ~.
+o.l 2;1 • 0. 1 j
.2-o . 4 7 .t-O, h
.,
o. lla - O .• } I. -o.09J.0.41
-O . ilK .. 0. 10
.. 2_ 0 ,o'l+o , n
As the elements of T (~ are directly identified in the form de~cribed in Eq. 19, the parameters of the continuous time model can be computed from Eqs. 15, 17, 18. When the modal characteristics associated with a finite element model are of interest, Eq.14 has to be taken into consideration also. The partial fractions in the elements of the transfer matrix Tv(Z) can be related to the natural frequencies of the vehicle. Denote f . a natural frequency as~ociated with the denominator polynomial Z +Pijz+ + Pzj· Then a direct relation can be given to compute fj from the discrete time model parameters,
Computing the frequencies f. associated with the identified transferJmatrix, the model can be explained as follows. The estimated frequencies can be clustered into four groups. The values around 11.5-12 Hz can be related to the second bending frequency of the body, and this. frequency appears in each outputs. The frequencies around 10.3-10.4 Hz appear in stress outputs Y2t and in axle accelerations Y3t' Y4t. These can be related to the natural frequencies of the front and rear axles. The frequency 9.1 Hz occurs only in the stress output Ylt' and it approximates the naturalfrequency of the engine. It c an be observed, that the input excites the system in the frequency range 9 Hz: f~ -12.5 Hz. The frequencies associated with the noise model poles lie in the range 3.8 Hz~ f ~ 10.5 Hz. Thus it can be expecN ted, that the model of vehicle given by the transfer matrix T (z), gives a good approximation in the fr~quency range 9Hz~ f~ £ 12.5 Hz, and for lower frequencies the effect of noise is dominant. This noise effect can be rather explained by the neglec ted nonsignificant low-frequency dynamic modes then by measurement errors. To analyse the above phenomena the logarithmic amplitude, 20 log I s .. (f)1 and the 1.J
amplitude-phase diagrams
Sij (jf)
were
computed by the use of continuous transfer matrix T (s) and were also estimated by cross-spe~tral method using 2000 samples. For the illustration of the results these diagrams for a stress output Y2t and for an acceleration output Y4t are shown on Fig.3. The comparison of ehe computed and estimated diagrams supports the above explanations. In order to illustrate the fit of the rrodel in the time domain, 100 samples of the measured Y2t' Y4t processes and of the corresponding outputs computed by the identified model are plotted on Fig.2.
P.
188
Mi c helbc r~e r , 1'/
ni.
CONCLUSION The paper discussed the identification of transfer matrix model for air suspension bus structures. The identified models can be used to determine the dynamic stress e s and stress st~tistics for arbitrary road profile excitation and speed levels within the examined range. REFERENCES 1. Michelberger,P., Bokor, J.,Keresztes A. and Varlaki P. (198~. Multivariable identification for commercial vehicle stress analysis. American Control Conference, San Francisco, California,
Fiq.Z-Plots of measurements and of identified model outputs. model ---- measurements
outputs
2. Ross, G.R. (1971). Synthesis of stiffness and mass matrices from experimental vibration modes. SAE paper. No 710787. 3. Natke, H.G. (1983). Einfilhrung in Theorie und Praxis der zeitreichen - und Modalanalyse. Braunschweig ; Wiesbaden: Vieweg. 4. Bokor,J. and Keviczky,L. \1984). Structure and parameter estimation of MIMO system using elementary sub-system representation. Int. Journal of Control., Vol 39., No.S. 5. Keviczky,L. Bokor,J. and Veres S. \1984). Strong consistency of ML estimators using partial fraction and elementary subsystem representation of multivariable system. 9 th • IFAC World Congress, Budapest 6. Bokor,J. and Pete, J. (1981).Minimal realisation procedure from partial fraction representation of a rational transfer matrix. Problems of Control and Int. Theory, Vol. 10, 7. Meirovitch, A. and Baruh,H. (1982). Identification of distributed-parameter vibrating systems. Proc . 6. IFAC Symp. on Ident, and System Par.Est., Washington 8. Michelberger,P. Keresztes, A. Bokor,J. and Varlaki, P. (1984) Dynamic Modelling of Commercial Road Vehicle Structures from test Data. XX FISITA Congress, Vol. 6-12. pp.4.96-4 . 103 . Wien.
Y, Y2 stresses
Y3 accelerations Y4 Fiq.l .- Locations of stress gauges acceleraneters.
am
Pig.) Oro . . epeotra - - - oOlDJlUte4 117 tlw data - - . .t1mate4 n- tile dat.