- Email: [email protected]
Fundamentals of system identification in structural dynamics
REVIEW PAPER Fundamentals of system identification in structural dynamics H. Imai
Setsunan University, Osaka, Japan C.-B. ¥un
Korea Advanced Institute of Science and Technology, Seoul, South Korea O. Maruyama
Musashi Institute of Technology, Tokyo, Japan M. Shinozuka
Princeton University, Princeton, NJ 08544, USA Methods of identification for structural dynamic systems are reviewed in this report. Emphasis is placed on methods relevant to the identification of linear and nonlinear behavior of structures subjected to such environmental loads as ground motion due to earthquakes, wind-generated pressure and wind-induced ocean wave forces. The methods considered are the least square, instrumental variable, maximum likelihood and a method utilizing the extended Kalman filter. In order to verify the validity of these methods, numerical simulation studies are carried out utilizing mathematical models of a suspension bridge, offshore tower and building structure. On the basis of such simulation studies, the efficiency of these methods are investigated under several conditions of observational noise. 1. I N T R O D U C T I O N This study reviews existing methods of system identification relevant to the identification of dynamic behavior of structures under various environmental loads. The problem of system identification has become increasingly important in the area of structural engineering, particularly in connection with the prediction of structural response to adverse environmental loadings such as earthquakes, wind and wave forces 2°'2s'9'18 and also with respect to estimation of the existing conditions of structures for the assessment of damage and deteriorationlV'9'2. The general subjects of system identification originally began in the area of electrical engineering and later extended to the field of mechanical/control engineering. Various techniques have been developed. One can find general surveys on the subject in Hart and Yao 7, IFAC Symposium 1°, Kozin and Natke ~2 and Ljung 14. However, those methods available may not be readily or directly applicable to problems of structural engineering systems for the following reasons: (i) these systems are generally much larger in size and much more complex in behavior so that accurate mathematical idealization is not easy, (ii) the options for input-output observational data are usually very limited, (iii) observational data are generally heavily contaminated by measurement noise, and (iv) in the case of a damaged or deteriorated system, the behavior may be highly nonlinear. Therefore, for the
purpose of effective structural engineering applications, specialized techniques of system identification need to be developed. In general, system identification methods are classified as parametric and nonparametric. Parametric identification involves estimation of system parameters, while nonparametric identification determines the transfer function of the system in terms of analytical representation. Identification methods can also be categorized as those in the time domain and the frequency domain. In the time domain method, system parameters are determined from observational data sampled in time. On the ,Jther hand, in the frequency domain method, modal quantities such as natural frequencies, damping ratio and modal shapes are identified using measurements in the frequency domain. In this study, mainly parametric identification techniques in the time domain are investigated, since it is more relevant to the identification of structural parameters related to environmental loading and/or nonlinear behavior. In Section 2, methods for modelling structural systems are discussed with an aim to system identification. In Section 3, several identification techniques based on least squares, maximum likelihood and extended Kalman filter are described. In Section 4, numerical simulation analyses are given for different structural systems, i.e., a suspension bridge, offshore structure and structure with nonlinear hysteresis. The efficiency of various identification
1989 ComputationalMechanics Publications 162
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
Fundamentals o f system identification in structural dynamics: H. Imai et al.
methods is investigated under several noise conditions. One of these simulation analyses is unique in that the method for system identification is applied to obtain an equivalent linear system representing a nonlinear offshore structure• The result obtained by this new method of identifying the equivalent linear system has been compared with those of conventional methods. Computer programs have been developed for the various system identification methods discussed in this study and used for example studies• Documentation of these programs is currently underway and will be made available in the form of a technical report in the near future• Presently, an experimental study is also being carried out using laboratory models in order to verify the validity of the identification techniques demonstrated in this study. Upon such verification, field experiments will follow• The results of these experimental studies will also be reported in the form of technical reports in the near future•
3 4
state space model ARMAX model.
In the following, it is shown that most of the structural systems can be modelled by means of ordinary differential equations (ODE) and that the remaining three models are derived from the ordinary differential equation. Example 1: Multi-story building ~9
Assuming a shear building model, an n-story building can be modeled as shown in Fig. 1. The equation of motion of this structure under ground excitation is written as: mi~ ~+ c , ( 4 , - ~ _ , ) - C,+ l ( ~ i + l
- - ~ i ) -~- k i ( ~ i - -
-k~+1(~+1-~)= -mi2 o
2. STRUCTURAL SYSTEM AND MATHEMATICAL MODEL Civil engineers deal with many types of structural dynamic systems: for example, multi-story buildings, suspension bridges, nuclear power plants, offshore structures, etc. The dynamic characteristics of these structures can be described by mathematical models. A variety of models have been developed for different purposes. Models commonly used in structural and systems engineering are the following:
~i-1)
(i= 1,2 . . . . . n)
(1)
where ~, 4i and ~'i are the relative horizontal displacement, velocity and acceleration of the ith floor to the ground, 20 is the horizontal ground acceleration, and m~, % ki are the mass, damping, stiffness coefficients, respectively• In vector-matrix notation, equation (1) can be written as follows: M 2 + C i + Kz = Lu
(2)
where z={~l,~: . . . . . (,}T and u = - 2 o. M is an n x n diagonal matrix with Mii = ml. L is an n x 1 matrix with L~= ml. C and K are both n x n symmetric matrices n
1 2
ordinary differential equation (ODE) transfer function C=
~n
m
Cl -}-£2
--C2
0
- - C2
C 2 "Jr-C 3
-- C3
0
- c3
c 3 + c4
......
0
. . . . . .
0
......
0
Ill . . . . . . . . . . . .
0
0
0
C N _ I -~-C N
....
--CN_
cN- 1
1
cN
(3a)
~n-i ~n-2
K= i
"k1 +k 2
-k 2
0
......
0
-k 2
k2+k 3
-k 3
......
0
0
- k3
k 3 q- k 4
......
0
. . . . . . . . . . . . 0
0
0
kN_l + k N ....
"
--kN_l
kN- 1
kN
(3b) 2
•
Example 2: Suspension bridge 2°
To investigate the dynamic characteristics of a suspension bridge such as aerodynamic instability, an idealized two-dimensional model of the bridge section as shown in Fig. 2 is frequently employed. The equation of motion of the model can be expressed as M ~ + C2 + Kz = L f
x
Fig. 1.
g
Shear beam structure model o f multi-story building
(4)
where z={h,a} z, f = { u , v } r with a=pitching motion, h = heaving motion, and u and v fluctuating components of the wind velocity in the horizontal and vertical directions, respectively. M is the mass matrix. C and K are
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
163
Fundamentals o1" system identification in structural dynamics: H. lmai et al.
/
where
a{t
wind
7 +
/
v(t)~/~. U
+ u(t)
O : IO- 11 (r)
The measurement y of the output z may be subjected to a measurement noise v, i.e.
"///////////////////////////////////////////////~
y=z+v=Hx+v
Fig. 2. Bridge deck model the damping and stiffness matrices aerodynamic effects.
including
the
Example 3: Offshore structure 25,1s An idealized model for the dynamic analysis of an offshore structure is shown in Fig. 2. The equation of motion of the structure is written by
where H = [I, 0]. The set of equations (9) and (10) is called the state-space model in continuous form. Usually, ~ and v are assumed to be white noise vectors with zero mean and variances Q and R, respectively. In recent years, most signal processing is accomplished by digital computers which cannot handle continuous time signals. Therefore, the output signal has to be in discrete form ( t = 0 , At, 2At . . . . ) where At denotes the sampling interval. Thus, for a linear system, the state space model is discretized as follows 2°
(Mo + C~t)2 + C2 + K o z : C~t~ + C o ( + - ~ ) l + - ~ I (5) where z is the vector of horizontal displacement, ?, ~ are the vector of horizontal wave particle velocity and acceleration, Mo, C, K o are the matrices of structural mass, damping, and stiffness, and C M and C o are the diagonal matrices containing, respectively, the inertia and drag coefficients associated with the wave force acting on the structure. In the case of linear structural systems as in equation (2), the equation of motion can also be represented by using the transfer function as Z(s) = G(s)U(s)
(6a)
G(s) = ( M s 2 + C s + K)- 1L
(6b)
where Z(s), U(s) are the Laplace transforms of z(t) and u(t), and G(s) is the transfer function of the system. The state space model in continuous form can be derived for linear and nonlinear systems. For example, by defining the state vector as
X: ~XI'~ {.X2) : {:}
(10)
x(i + 1 )= Fx(i) + Gu(i) + w(i)
(1 l a)
y(i)=Hx(i)+v(i)
(1 lb)
where F = e Am
G:
e At dt B
and
x(i)-x(iAt),
y(i)=-y(iAt)
u(i)-~u(iAt),
v(i)~v(iAt)
where it is assumed that the input u(t) is piecewise constant within each sampling interval, i.e.
u(t)=u(iAt),
iAt<=t<(i+l)At
(12)
The system noise w(i) which is defined by
(7)
w(i)=
~
(i+ 1)At CA[(/+ I)AT q D ~ ( z )
¢
dr
(13)
d i At
Equation 2 can be transformed into a state equation as
Xl}
{ f ( X l , X2 ' U)} +
is a white noise sequence with zero mean and the covariance matrix
w
x2
Qw=E[w(i)wr(i)] = x2
)
~0A!eA~DQDT[eA~]T dr
(14)
F 0 q The state space model may also be represented in innovation form 15, i.e.
where w and ~ are system noises which account for the unmeasurable input disturbances and errors in modelling. If the system is linear as in equation (2), equation (8) can be rewritten as
~=Ax+Bu+D{ 164
(9)
x(i+ l)=Fx(i)+Gu(i)+Fe(i)
(15a)
y(i) = Hx(i) + e(i)
(15b)
where x(i) is the estimation at t = i At, F is the so-called steady state Kalman gain and e(i) the innovation process.
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
Fundamentals of system ident~cation in structural dynamics: H. Imai et al. The innovation model is preferable over the general state-space model, since it has fewer parameters than the latter. When the system is controllable and observable, we can eliminate x(i) from equation (15) to obtain the following ARMAX model 2°'15
~F= [0(3), 0(4) . . . . . O(N)] r
(21c)
Then, equation (18) can be rewritten as
y(i)=oro(i)+~(i)
for i = 3 , 4 . . . . . N
(22a)
and further into matrix form as y(i) = F l y ( i - 1) + F2y(i--2) + G~u(i- 1 ) + G 2 u ( i - 2 ) +e(i)+J~e(i- l)+J2e(i-2 )
Y = ~F0+ It(3), 5(4) . . . . . ~(N)] r
(16)
If the measurable input u(-) is missing, the model is called an ARMA model. The mathematical models introduced above are used for different purposes. For conventional problems of the analysis and design of a structural system, ordinary differential equations and transfer functions have commonly been used for a long time. For the identification of modal quantities, the transfer function has been also widely employed. However, as for the identification of structural parameters, state space equations and ARMAX models are preferable because the identification techniques recently developed are based on sampled data in time.
3. ALGORITHMS FOR PARAMETER
(22b)
Then, from equation (22), the least squares estimate 0 which minimizes equation (19) can be obtained as OLs = [ v T v ] - ' [ v r v ]
(23)
It is noted that in general the least squares estimate is biased; i.e., OLs does not converge to the true value 0 as N approaches infinity. This is a major drawback of the least squares method. To overcome this difficulty several improved methods have been developed. Among them, an instrumental variable method 24'4 and a generalised least squares method should be mentioned. The least squares method is also available in sequential form ~s. Let 0~ denote an estimate of 0 based on the data {y(1), y(2) . . . . . y(i), u(1) . . . . . u(i)}, then O, can be estimated sequentially by
ESTIMATION There exists many different algorithms of parameter estimation. Most commonly used algorithms are least squares method, maximum likelihood method, extended Kalman filter, and variations of them, which will be briefly discussed in the following.
with
3.1. Least squares method 4'2°'~4
Pi = Pi-1
Consider a system described by an ARMAX model, equation (16), and define the equation error vector as
priori.
e(i)=e(i)+Jle(i-1)+J2e(i-2 )
-
-
O,=Oi_~ +K,(yr(i)-or(i)O~-~)
(24)
Ki=P,_ 10(i)(1 +Or(i)Pi_,O(i)) -~
(25)
P,- 10(i)( 1 + 0r(i)P~ - ~0(i))-'0r(i)P,-~
(26)
where the initial values 0 o, Po, 0(1) should be given a
(17)
3.2. Maximum likelihood method 1'4'2°'14 Then equation (16) can be rewritten as y(i) = F l y ( i - 1)+Fzy(i-2)+G~u(i- 1)
+ G2u(i- 2) + e(i)
(18)
Let P(y[u, 0) denote the conditional probability density, where u = [u(1 ), u(2) . . . . . u(N)] r. For given data [Y, u], one may regard P(ylu, 0) as a function of 0, which is called the likelihood function. An estimate which maximizes the likelihood function or equivalently
The least squares method is based on minimization of N
L = log P(y[u, 0)
(27)
N
J= E ~(i)%(i)= ~ [l~(i)ll 2 i=3
is called the maximum likelihood estimate and denoted by
OML,
i=3
N
As for the linear discrete-time system, the maximization of equation (27) corresponds to the minimization al
= ~ Ily(i)-Ely(i- 1)-E2y(i-2) i=3
- G l U ( i - 1)-G2u(i-2)[] 2
(19)
e(i)er(i)
J=det
with respect to
Li
0 = [F~, F2, G1, G2] r
(20)
where N is the number of data points. To solve the minimization problem, one may define a vector 0(i) and matrices Y and W as
0(i)= [y(i-- 1) r, y(i--2) r, u ( i - 1) r, u(i--z)r] r (21a) Y = [y(3), y(4) . . . . . y(N)] r
(21b)
(28)
= 1
where e(i) is a one step ahead prediction error which is defined by e(i) = y(i) - y(i[i- 1)
(29)
and y(i[i- 1)-E{y(i)ly(1), y(2) . . . . . y ( i - 1),
u(1), u(2) . . . . . u ( i - 1), 0}
(3o)
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4 165
Fundamentals of system ident~cation in structural dynamics: H. lmai et al. The prediction error e(i) can be computed sequentially according to e(1)=y(l)
P(k +
l lk+ l) = [I- K(k + 1)H]P(k + 1Ik)[l- K(k + 1)H] r + K(k + 1)Vt,K(k + 1)T
(38)
where K(k + 1) is the Kalman gain matrix which is defined as
e(2)=y(2)-Fly(l)-Gxu(1)-J~e(1)
e(i)=y(i)-Fly(i)-F2y(i-2)-Gxu(i- 1 ) - G 2 u ( i - 2 ) -Jle(i- 1)-J2e(i-2) i=3,4 ..... N (31) The problem of minimization of equation (28) may not be solved analytically in general. Thus nonlinear programming algorithms have to be used, e.g., Newton-Raphson method, Davison method 3, etc. While the maximum likelihood estimate has a nice property that it is consistent and asymptotically efficient, it usually requires a substantial amount of computational time. To circumvent this disadvantage, several approximate methods are proposed, e.g., recursive maximum likelihood method s , approximate maximum likelihood method 6, etc. Most of them are given in sequential form.
K(k+ 1 ) = P ( k +
1]k)HT[HP(k+ llk)HT+V,,] - ' (39)
The extended Kalman filtering technique described above can be applied to the identification of system parameters as follows. Consider a dynamic system described by a state equation as in equation (8), /~= f(x, u, 0; t) + w(t)
(8a)
where x is the state (response) of the system, u is the input excitation and 0 denotes system parameters. Then, by defining a new augmented state vector as
x={;}
3.3. Extended Kalman filter 6'25"14
,40,
The basic algorithm of the extended Kalman filter is a recursive process for estimating the optimal state of a nonlinear system based on observed data for the input (excitation) and output (response). It can be summarized as follows. Consider a general continuous state equation described by
a new nonlinear state equation can be constructed from equation (8a) as
:~(t) = g(X, u; t) + w(t)
(32)
The observation equation corresponding to equation (10) can be obtained as
(33)
y(k)=[H,O]{o(k)}+v(k)
X--{f(X'OoU;t)}+{W(ot) }
(41)
with observations at t = k At x(k)
y(k) = HX(k) + v(k)
in which X(k)=state vector at t=k At, y(k)=observational vector at t = k At, v(k) = observational noise vector with covariance of V~, w(t)= systems noise vector with covariance of Vw and H--matrix associated with observations. The predicted state X(k + 1[k) and its error covariance matrix P(k +l[k) can be evaluated as X(k +llk) = E{X(k + 1)[y(1), y(2) . . . . . y(k)} =X(klk)+
~
(42)
By applying the extended Kalman filter to equations (41) and (42), system parameter 0 can be estimated recursively as part of the state vector. The estimates obtained by the extended Kalman filter method may be in general biased or divergent 22'23. Some modified algorithms are proposed to improve the convergence property 13'2t'8. In this study, the weighted global iteration procedure proposed by Hoshiya and Saito 8 is also used in example studies.
( k + 1)At
g(X(tlk),u;t)dt
(34)
d k At
P(k+ l l k ) = O ( k +
1,k)P(klk)OT(k+ 1,k)+V w (35)
where E{AIB } is the expected value of A conditional to B and O(k + 1,k) is the state transition matrix which can be approximately obtained as
O(k+ l,k)~-I + At [[-0gi(x(t)' ~ t) lxl,):x/klk)
(36)
for small At. Then, the filtered state X ( k + l [ k + l ) and its error covariance matrix P(k + l lk +1) can be estimated as X(k+ ilk+ 1)= E{X(k + 1)ly(1 ), y(2) . . . . . y(k+ 1)} = X ( k + llk)+ K(k + l) x [y(k+ 1 ) - H X ( k + llk)]
166
(37)
4. N U M E R I C A L EXAMPLES AND DISCUSSION In order to investigate the accuracy and efficiency of the various system identification techniques discussed in the preceding sections, example analyses have been carried out for four different cases: (i) an idealized suspension bridge model for wind loading, (ii) a simplified model of an offshore structure for wave forces, (iii) a singledegree-of-freedom structure with bilinear hysteric characteristics under seismic excitations, and (iv) an equivalent linear system for an offshore structure model. In each case, simulated data for input and output time histories are utilized for parameter identification purposes. The estimated system parameters and responses are compared with exact ones which have been assumed a prioiri. In the last example, the system identification method is applied for determination of the equivalent linear system of an offshore structure model.
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
Fundamentals of system identification in structural dynamics: H. lmai et al. 4.1. Suspension bridge model 2° The idealized structural model shown in Fig. 2 is uscd. The equation of motion for wind loading is given in equation (4). Assumcd values of thc systcm parameters are shown in Table 1. Simulated time histories of the wind velocity fluctuations, u(t) and v(t), are shown in Fig. 4. Thc (simulatcd) observation timc histories of the bridgc motion, Yh(t) and Y~(t) are shown in Figs 5 and 6. The level of observational noise included in Yh(t) and Y~(t) is assumed to be 10% of the structural response in the root mean square (RMS) value in both cases. For the purpose of system identification, the ARMAX model has been utilized in this example. Parameter estimation is carried out by using the least square (LS) method, instrumental variable (IV) method and maximum likelihood (ML) method. The estimated parameters arc compared with (assumed) exact ones in Table 1. The response time histories rccomputed using estimated parameters, ~(t) and ~(t), are also compared with those of the observation in Figs 5 and 6. From the results in Table 1 and Figs 5 and 6, one can scc that the IV and ML methods give good estimations of the parameters as well as structural responses. On thc other hand, the results from thc LS method are found to be unsatisfactory. The poor rcsults of the LS estimates are due to the fact that the observational
5.00 +00
(~/~)
u(t)
0.00 +00
~
-5.00 +00
-
-
.
-
-
.
i
.
-
-
.
-
-
-
-
.
-
~
L
~
5.00 *00
(~/8~¢) 0.00
*OO
-S.OO +OO
.--i.-.---J-___ 0.00 25.50
50.00
75.00
TIME
(SEC)
IC0.00
125.00
Fig. 4. Simulated time histories of wind velocity fluctuations (mean horizontal wind speed = 10 m/sec)
II MASS 1 x
M.SL.
Z.O0 -Of
~7
P MASS 2
ty
F3---->
t MASS 3
r~--i->
|
% '
I MASS
MASS
(m)
o.o
-2.00
.oo
zhC~)
v,,wvv.vvvvv.vw,
VV/
-01
2.00 -01
(m)
I l MASS 6
v
~(t)LS
0.00 *00
v "
~
V -V
-v-
" "V%,,-'
--V xs,,'-,
/\ "/// X "///////
.',//////
"//////////
a) OFFSHORE TOWER
Fig. 3.
I
-2.00 -01
I
T
p
2.00 -0t
(m)
~(t)iv
b) M O D E L
Fixed offshore tower and model
0.~8 -0C
' vvvvvv vvVVV'vV vwvv -VVV/
Table 1. Exact and estimated parameters of suspension bridge P ....
ters
Exact Values LS
ML
M-tK
~1
M-'C
M-'L
{
oz
3.70 -0.12
0.95 18.25
0,128 -0.020
0.050 0.060
0.0020 0.0003
-0.057 0.012
0,0309 0.0022
5.50 0.14
-0.75 19.40
3.194 -0.013
-0.563 0.606
0.0037 0.0003
-0.122 -0.013
0.0160 0.0009
3,65 0.13
1.36 18,31
0.128 ,.0,020
0.148 0.058
0,0015 0,0004
-0.053 0.012
0.0300 0.0023
3.72 -0.11
-1,62 18.19
0.122 -0.020
1.325 0.080
0.0019 0.0011
-0.062 0.012
0.0327 0.0021
Note : i. Observational noise is assumed a.s I0 % in R M S value 2. o. = standard deviation of Zl and z2 3. U n i t : [ M - t K ] = 1/see; [ M - ' C ] = 1/see; I M - ' L ] = llsecCfirst row), X/m/see(second row); [ a , , ] = m and [t%~ I = rad/sec.
-2.0C -0: 2.00 -01
(m)
E(t)~u "
0 . 0 0 *00
-2.00 -01 ~ O.Ol
25.00
50.00
' 75.20
IOO.O0
125.00
TIME (sec)
Fig. 5. Time histories of observation Yh(t) and estimated responses h(t) for heaving motion of suspension bridge Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
167
Fundamentals of" system identification in structural dynamics: H. lmai et al. I.O0 -02
It is assumed that time histories of the displacement at the two discrete masses are available. Identification of the system parameters is carried out by using the extended Kalman filtering with the weighted global iteration procedure s,9. Table 2 shows the (assumed) exact and estimated values of the parameters for two sea states corresponding to two wind speeds of 25 and 75 ft/sec. Two different conditions of the observational noise (5 and 10% in the RMS values) are considered. The estimated values have been obtained through five global iterations. Figure 7 shows input time histories of wave particle velocities and
(t)
(r~a) 0.00 -00
-[.O0 -02 1.00 -02
(tad)
~(t)us
O.O0 -00
,,
.,J..A,h,~lL,,,h^..,, ..i..d..,h.
,~,~,"q. ] r ',, ~ " " H " ~ " " ~ ) W .
,ll,,Aill,. ,,,~)~O ,'qlr~'" I....
Table 2.
Exact and estimated parameters o/" +~[hhore structure model
~
l't/sec 2
-1.00 -02
JLt Js:
J
(~,,)
~.00 -02
(tad)
~(t)£7
D (~)
0.187 -0.022
0.4 0.0
J~2
-0.1u o.o
JT~'
0.150
0.4
Dtt
@060
0.1
@060
0.1
Ktl
3.7,50
5.0
K21
-0.750
0.0
D,,
K
(.-fi~) x,~
0.CC -,3G
L
-3.7so
0.226 .0.008 -0.147 0.127
oo
K.;
2.500
4.0
Ltl L~2
0,400 0.500
0.2
75
25
0.3
0.166 .0.005 -0.073 0.132
0.049 @064
0.261 0.028 -0.168 0.103 0.037 0.071
0.062 O.OSg
0.150 0.005 -0.043 0,125 O.Oro~ 0.058
3.743 -0.721 -3.780 2.465
3.717 -0.705 -3.779 2.449
3.779 -0,753 -3,772 2.503
3.794 -0.765 -3.765 2.520
0.408 0.494
0.419 0.486
0.422 0.405
0,446 0,492
5
10
5
10
Observatlon~ Noise Level
(% of response in R.M.S,)
OC -'3Z
- i.
1.00
-02
(~ad)
Note : #o = initial guesses of parameters 0, = Estimates by the E x ~ n d e d Kalman Filter after Fifth Global Iteration
~(t)ML
(a) Node t 0.00 -00
~0,00
(ftl~¢) O.CO
-~.00 -02 ~ 0.00
25.00
, ZO.O0
75.00
a___ lO0.CO
125.00 -~O.CO
TIME (SEC)
I0.00
Fig. 6. Time histories of observation Y,(t) and estimated responses ~(t) .for pitching motion of suspension bridge
(#/~) O.CO
vector qt(i) in equation (21a) is not statistically independent of the moving average noise term r(i) in equation (17). Further numerical investigation indicates that the results from the ML method are less sensitive to observational noise than those from the IV method.
~
I
r
I
20.00
40.00
50.00
'
f
I
90.00 time (see)
I00.00
(b) Node 2 '.O.CO i
'
(#/se~2) |
,
'
~
'
J
h (t)
'
~
'
.
--
4.2. Offshore structure 25.1s In this example, an offshore structure idealized as a two-degrees-of-freedom system is utilized. The equation of motion for wave loading, equation (5), is simplified as -[0.00
i + J~- D{(¢-i)l¢-
~1} + K z = L~
(43)
where J = [ M o + C M ] ~C, D=[Mo+CM]-ICD, K= [Mo+CM] 1Ko and L = [ M 0 + C M ] - I C M . For system identification, equation (43) is transformed into a state-space model by using the augmented state vector defined as {X}
=
{Z l Z2 Z1 -72 J l l J 2 l J12 J22 D l l D2z
Kll K21 K12 K22 LI~ 168
L22} T
(44)
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
IO.O0
I
~
'
'
I
[
,
,
,
L
,
!
'
)
,
,
~
,
4
/ O.CO ~
-I0.00
1
0.00
.
,
)
)
)
20.00
~.0.O0
60.00
~
)
,
80.00 ICO.OC time (see)
Fig. 7. Simulated time histories o,1"wave particle velocities and accelerations (wind speed = 75 ft/sec)
Fundamentals of system identification in structural dynamics: H. Imai et al. accelerations. Figure 8 shows the exact response z l(t ) and z2(t ) obtained using the assumed (exact) parameters and simulated response observation, }I1(t) and Y2(t), as well as the estimated responses, ~(t) and ~2(t), computed based on the identified parameters. Figure 9 demonstrates the convergence of recursive estimation of system parameters by the extended Kalman filtering algorithm. The results in Table 2 indicate that the extended Kalman filtering technique yields fairly good estimates even under relatively severe nonlinear hydrodynamic
(a) First Iteration
(~/s ~t) O. 30
' I
'
2g(k/k)
,
'
~
~
,
0.!5
O. O0
I
r
I
(a) Node 1 7.50
,
~
(it)
,
,
,
:
,
,
,
,
~ O. 00
20. O0
40. O0
60.00
0.00
80. OO
I00.00
time (~ec)
(b) Fifth Iteration -7.50
I
7. m.O
'
'
I
0.30
'
I
'
I
0//~)
,
i
f
fit)
D~1(t)
=
2g(~/k)
i , i
0.15 0.00 I
0.00
O/.ed)
-7.50 7. ~ 0
"
,
I
(/t)
,
I0.00
~
i
,
I
'
l
,
I
t)=2~3i~tk)
5.00 0.00
0.00
-7.50 0.00
•
20.00
40.00
60.00
80.00 time
(Set~)
,
,
,
,
,
l
,
~
,
~
,
,
,
,
,
,
(/t) O. 00
-5.
O0
r
5.00
,
~
(/t ) 0.00
-5.00
I
(It)
t
'
I
'
I
i
'
i
"
'
I
,
I
I
'
i
'
i
i
i 20.00
60.00
80.00
I00,00
time (see)
Fig. 9. Estimated parameters of offshore structure as augmented state variable (wind speed= 75 ft/sec; noise level= 10% of response in RMS)
+ 2~co~+ (oZg(ze. cO= - 2o(t)
~
i 40.00
|
80.00
A single-degree-of-freedom structure which exhibits bilinear hysteretic behavior under seismic excitation is identified. The equation of motion is written as
0.00
-5.00 0.00
!
60.00
4.3. Structure with bilinear hysteresis9
i
5.00
i 40. O0
loading conditions for both of the observational noise conditions. The initial estimates of the parameters were chosen to be quite far from the exact values. The error covariance matrix of the initial estimate has been taken fairly arbitrarily as a diagonal matrix with large values of the diagonal elements. However, it has been found that the estimated values converge to reasonable ones as the time step increases and as the iteration proceeds (Fig. 9). The estimated response obtained by using the identified parameters is found to be virtually identical to the exact response (Fig. 8).
(b) Nod~ 2 O0
20. O0
ICO.OC
Fig. 8-a. Exact, observed and estimated responses of offshore structure (wind speed= 75 ft/sec; noise level= 10% of response in RMS)
5.
' O. O0
I00.00
time (see)
Fig. 8-b. Exact, observed and estimated responses of offshore structure (wind speed= 75 fi/sec; noise level= 10% of response in RMS)
(45)
where co and { are the natural frequency and damping ratio, respectively and g(Ze,C0 defines the bilinear hysteresis with z e being the yielding displacement and c( being the ratio of post-yielding stiffness to pre-yielding stiffness as shown in Fig. 10. For the purpose of system identification, equation (45) is rewritten into a nonlinear state equation by using the state vector defined as {X} = {z ~ ~ co z e ~}r
(46)
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4 169
Fundamentals of system identification in structural dynamics: H. lmai et al.
g(Ze,a)
The results in Table 3 as well as Figs 12 14 indicate that extended Kalman filtering, particularly with the weighted global iteration procedure, estimates the system parameters and hysteretic behavior remarkably well even for the case of very severe material nonlinearity as shown in Fig. 10.
D
Fig. 10.
z
4.4. Identification of equivalent linear system 18 In this example, equivalent linearization by means of a system identification technique is discussed. In general, the nonlinear characteristics of actual structures under severe environmental loading, such as earthquakes, wind and waves are very difficult to model and analyze. Hence approximations in modelling are inevitable. For instance,
Bilinear hysteresis
400.00
~,(:)
(,~/,.,=)
f
I
T
C.OO
wlJ
"
O.nO 4 -!5.C,.0
I .=.30
0.00
I 10.SO
; 1.=.00
2O.C0
t~.e (s~) t5.CO
-40C.O0
o. oo
'=.oo
to. oc
15.00
2O.C0
Fig. II. Input ground acceleration (1940 El Centro earthquake, N S component)
O.CO
Table 3. Exact and estimated parameters of structure with bilinear hysteresis Parameters
-~5.00 0.00
Estimates
Exact
Initial
Va|ues
G uesse~
First Iteration
Fifth Iteration
o.1
o.s
0.096
0.098
w(rad/sec)
3.14
1.0
3.157
3.147
z,(cm)
3.0
t.o
2.979
3.025
a
0.1
0.5
0.090
0.088
Fig. 12. structure
5.00
Observed
50.00
'l
tO.O0
15.00
20.C0
and estimated displacement
I"
of
T
(c=Is~) Identification of parameters is carried out using the extended Kalman filtering technique with the weighted global iteration procedure. As an input ground acceleration, the 1940 E1 Centro earthquake record (N-S component, Fig. 11 ) is utilized. It is assumed that time histories of the structural displacement and velocity are available (Figs 12 and 13). The observational noise levels are taken as 10% of the structural response in the RMS values. Table 3 shows the exact and estimated values of the system parameters. Figures 12 and 13 also show time histories of the estimated responses obtained after the fifth global iteration. Finally, Fig. 14 compares four different hysteretic characteristics. They are (a) true, (b) that obtained from observational data without the extended Kalman filtering procedure by g(Ze, 0~)-
( - j~(t) - ~ - 2 ~o~) (1)2
-50.00 0.00
5.00
10.00
i5.00
20.00
tLme @~)
50.00
(c=/'s~)
5(t)
0.00
(47)
and the remaining two, (c) and (d), obtained by using extended Kalman filtering with weighted global iteration.
170
0,00
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
-50.00 0.00
Fig. 13.
I
1
I
5. O0
I0. O0
IS. GO
20. GO
Observed and estimated velocity of structure
Fundamentals of system identification in structural dynamics: H. lmai et al.
(a) 12.00
Q
(6) 12.00
L
0.00
0.00
-12.00 -12.00
0.00
12.00
-12.00 -12.00
D i s p l a c e m e n t (cm)
O. O0
12.00
Displacement (cm)
(c)
(d)
12.00
IZ.O0
0.00
0.00
A
-12.00 -12.50
0.00
12.00
-12.00 -t2.00
Displacement (cm)
(a) True Hysteresis (e) E s t i m a t e d by First Iteration
Fi9.14.
O. O0
12.00
Displacement (cm)
(b) Estimated by Conventional Method (d) Estimated by Fifth Iteration
True and estimated hystereses
in the case of offshore structures, the nonlinear hydrodynamic drag force is represented using drag coefficients which are determined semi-empirically based on representative values of the wave particle velocity and size, shape and surface roughness of the structural members. Because of the uncertainty and computational complexity associated with the nonlinear terms, the equivalent linearization technique has frequently been used in many analyses. The conventional method performs linearization by minimizing the error resulting from linearization and requires considerable numerical effort. In the present example, the idealized offshore structure model used in Section 4.2 is utilized. From the nonlinear equation of motion, equation (43), a linearized equation is
derived as
+ J*~ + K*z = L*~ + D*(~ - 2)
(48)
where J*, K*, L* and D* are the parameters of the linearized system. As in the previous example, time histories of the wave particle and acceleration and structural displacement are assumed to be available at each node. The measurement noises are assumed to be 5% of the structural responses in the RMS values. The linearized parameters for the hydrodynamic drag force as well as for the other system parameters are identified by using the extended Kalman filtering technique with weighted global iteration. Two sea states corresponding to two wind speeds of 25 and
Probabilistic Engineerin9 Mechanics, 1989, Vol. 4, No. 4
171
Fundamentals o/" system ident!fication in structural dynamics: H. lmai et al.
75 ft/sec are considered. Table 4 shows the assumed exact and estimated parameters. It should be noted here that the estimated values under each wind speed are the ensemble average of five sets of estimations arising from five sets of statistically identical but individually different wave particle motions, b(t) and/)(t), and observation time histories Y(t) generated by simulation. The estimated coefficient matrix D* has also been compared with the values computed by the conventional equivalent linearization procedure ~6. Figure 15 compares the observed response and estimated responses by two linearization methods. From the results in Fig. 15, it can be seen that the equivalent linear system obtained by using system identification can simulate the responses of the actual nonlinear system exceptionally well. The results in Table 4 indicate that for the case of small wave conditions where the nonlinear effect is not so significant, the linearized parameters obtained by the two linearization methods are quite different. However, as the wave condition becomes more severe and the nonlinear term becomes more important, the estimated values by two methods are found to be in good agreement. The above results indicate that the present method of equivalent linearization by means of system identification is considered to provide an efficient alternative to the conventional linearization method, in view of the fact that all the coefficient matrices are assumed to be unknown in the present method, while only the linearized drag coefficient matrix is considered unknown in the conventional method.
Table 4. Estimated parameters of equivalent linear system Parameters
J
(,,~)
(.-~-.) L D
(~) D"
( ,-,u, )
L
}
0.00
r
0.00
,
2 0 . O0
[
I
[
6 0 . O0
40. O0
80.30
O.187 -0.022
-0.111 0.150
0.295 -0.063 -0.172 0.249
0.258 0,080
J,2
0.400 0.000 0,000 0.400
3.73O -0,737 -3.766 2.488
3.944 -0.744 -4.234 2.578
K,j, K2,
3.750
0.097
K22
-3.750 2.54)0
5.0(}0 0.000 0,000 4.000
Ll L:
0.400 0.500
0.200 0.300
0.405 0.506
0,353 0.597
Dz|
0.060 0.060 0.300 0.300
o.o47(o.tm) o.o~(o.os0
o.3s'qo.,L,u) o.',so(o.33o)
z,,
D.
-0.750
D;
D;
In this paper, methods of parameter identification for structural dynamic systems are reviewed with emphasis on the identification of system parameters which vary with environmental loadings. For linear systems, the least squares, instrumental variable and maximum likelihood methods are reviewed. For nonlinear systems, a method utilizing the extended Kalman filter is discussed. Example numerical analyses are carried out for identification of the aerodynamic coefficients of a suspension bridge under wind loading, drag coefficients of an offshore structure under wind-induced wave forces and yield displacement and stiffness ratio of a structure with bilinear hysteresis subjected to seismic excitation. From the numerical results, it has been found that the instrumental variable and maximum likelihood methods provide good estimates for a linear system, while the extended Kalman filtering technique with weighted global iteration yields excellent estimates for nonlinear cases. A method for developing an equivalent linear system by means of system identification is also presented and a numerical simulation study indicates that this method offers an efficient alternative to the conventional linearization method.
(P)
-5.C0
Jtz J2z
E a t i m a t ~ Valu4 gstlmated Value W = l S ft/s~: W=TSft/sec
5. C O N C L U S I O N S
I
{
Initial Gum
* T h e values in parenthesis Were computed by the con~ntional linearisation method assuming D i and D~ are the only unknowns.
(a) Node 1 5.C0
Linear Syatem
System
Jll K
Nonlinear
lCO. O~
ACKNOWLEDGEMENT (b) Node 2 5.C0
This work was supported by the National Center for Earthquake Engineering Research under grant numbers 87-3012 and 88-3011 under the auspices of the National Science Foundation master contract number ECE86-07591.
i
fit) ,
i
O.CO
REFERENCES
1 -5.00
~ 0 . O0
t 2 0 . O0
,
r qO. O0
~
l
r
riO. O0
I
,
80. O0
t 00. OC
2
time (see) .
..... ....
g'xac
~.
~.ar.J.mar.~.d by P r e ~ e n ~ L L n e a r i z a r . l . o a
- F. , ~t t ~a ~e ~ b y C o n v e n e . £ o . a l
L£,qe~r-Lza~Lon
Fig. 15. Exact and estimated displacement of offshore structure (wind speed=75 ft/sec; noise l e v e l = 5 % of response in R M S )
172
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
3 4 5
AstrSm,K. I. and Eykhoff,P. System identification - a survey, Automatica, 1971, 7, 123 162 Chen,J.-C. and Garba, J. A. Structural damage assessmentusing a system identification technique, Proceedings of the Workshop on Structural Safety Evaluation Based on System Identification Approaches, Lambrecht, West Germany, June 29 July l, 1987, F. Vieweg& Sohn, Braunschweig, West Germany, pp. 474-492 Dahlquist,G. and Bj6rck, A. Numerical Methods, Prentice-Hall, Englewood Cliffs, New Jersey, 1974 Eykhoff, P. System Identification - Parameter and State Estimation, John Wiley & Sons, New York, 1974 Gertler, J. and Banyasz, C. A recursive (on-line) maximum likelihood identification method, IEEE Trans. on Automatic Control, 1974, AC-19, 816-820
Fundamentals o f system ident!~cation in structural dynamics: H. lmai et al. 6
7 8 9
Goodwin, G. C. and Payne, R. L. Dynamic System ldemifcation: Experimental Design and Data Analysis, Academic Press, New York, 1977 Hart, G. C. and Yao, J. T.-P. System identification in structural dynamics, Journal of Engineering Mechanics, ASCE, 1977, 106(EM6), 1089 1104 Hoshiya, M. and Saito, E. Structural identification by extended Kalman filter, Journal of Engineering Mechanics, ASCE, 1984, 110(12), 1757 1770 Hoshiya, M. and Maruyama, O. Identification of nonlinear structural systems, Proceedings of ICASP-5, Vancouver, Canada, 1987, Vol. I, pp. 182-189
10
Proceedingsof 6th IFAC Symposiumon Identification and System Parameter Estimation. Arlington, Virginia, June 1982
11
Kashyap, R. L. Maximum likelihood identification of stochastic linear systems, IEEE Trans. on Automatic Control, 1970, AC-15, 25-34 Kozin, F. and Natke, H. G. System identification techniques, Structural Safety, 1986, 3, 269 316 Ljung, L. Asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems, IEEE Trans. Automatic Control, 1979, AC-24, 36-50 Ljung, L. System Identification - Theory for the User, Prentice-Hall, Englewood Cliffs, New Jersey, 1987 Ljung, L. and S6derstr6m, T. Theory and Practice of Recursive Identification, MIT Press, Cambridge, Massachusetts, 1983 Malhotra, A. and Penzien, J. Nondeterministic analysis of offshore structures, Journal of Engineering Mechanics, ASCE, 1970, 96(EM6), 985-1003
12 13 14 15 16
17
18 19 20 21 22 23
24 25
Natke, H. G. and Yao, J. T.-P. System identification approach in structural damage evaluation, ASCE Structures Congress '86, 1986, Preprint 17-1 Paliou, C. and Shinozuka, M. Identification of equivalent linear systems, accepted for publication in the Journal of Engineering Mechanics, ASCE, 1988 Shinozuka, M., Itagaki, H. and Hakuno, M. Dynamic safety analysis of multistory buildings, Journal of Structural Engineering, ASCE, 1968, 94(ST1), 309-330 Shinozuka, M., Yun, C.-B. and Imai, H. Identification of linear structural dynamic systems, Journal of Engineering Mechanics, ASCE, 1982, 108(EM6), 1371 1390 Song, T. L. and Speyer, J. L. The modified gain extended Kalman filter and parameter identification in linear systems, Automatica, 1986, 22, 59 75 Urain, B. Asymptotic convergence properties of the extended Kalman filter using filtered state estimates, IEEE Trans. on Automatic Control, 1980, AC-25, 1207 1211 Westerlund, T. and Tyss, T. Remarks on asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems, IEEE Trans. Automatic Control, 1980, AC-25, 1011-1012 Wong, K. Y. and Polack, E. Identification of linear discrete time systems using the instrumental variable method, IEEE Trans. on Automatic Control, 1967, AC-12, 707 718 Yun, C.-B. and Shinozuka, M. Identification of nonlinear structural dynamic systems, Journal of Structural Mechanics, 1980, 8, 187-203
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
173
Setsunan University, Osaka, Japan C.-B. ¥un
Korea Advanced Institute of Science and Technology, Seoul, South Korea O. Maruyama
Musashi Institute of Technology, Tokyo, Japan M. Shinozuka
Princeton University, Princeton, NJ 08544, USA Methods of identification for structural dynamic systems are reviewed in this report. Emphasis is placed on methods relevant to the identification of linear and nonlinear behavior of structures subjected to such environmental loads as ground motion due to earthquakes, wind-generated pressure and wind-induced ocean wave forces. The methods considered are the least square, instrumental variable, maximum likelihood and a method utilizing the extended Kalman filter. In order to verify the validity of these methods, numerical simulation studies are carried out utilizing mathematical models of a suspension bridge, offshore tower and building structure. On the basis of such simulation studies, the efficiency of these methods are investigated under several conditions of observational noise. 1. I N T R O D U C T I O N This study reviews existing methods of system identification relevant to the identification of dynamic behavior of structures under various environmental loads. The problem of system identification has become increasingly important in the area of structural engineering, particularly in connection with the prediction of structural response to adverse environmental loadings such as earthquakes, wind and wave forces 2°'2s'9'18 and also with respect to estimation of the existing conditions of structures for the assessment of damage and deteriorationlV'9'2. The general subjects of system identification originally began in the area of electrical engineering and later extended to the field of mechanical/control engineering. Various techniques have been developed. One can find general surveys on the subject in Hart and Yao 7, IFAC Symposium 1°, Kozin and Natke ~2 and Ljung 14. However, those methods available may not be readily or directly applicable to problems of structural engineering systems for the following reasons: (i) these systems are generally much larger in size and much more complex in behavior so that accurate mathematical idealization is not easy, (ii) the options for input-output observational data are usually very limited, (iii) observational data are generally heavily contaminated by measurement noise, and (iv) in the case of a damaged or deteriorated system, the behavior may be highly nonlinear. Therefore, for the
purpose of effective structural engineering applications, specialized techniques of system identification need to be developed. In general, system identification methods are classified as parametric and nonparametric. Parametric identification involves estimation of system parameters, while nonparametric identification determines the transfer function of the system in terms of analytical representation. Identification methods can also be categorized as those in the time domain and the frequency domain. In the time domain method, system parameters are determined from observational data sampled in time. On the ,Jther hand, in the frequency domain method, modal quantities such as natural frequencies, damping ratio and modal shapes are identified using measurements in the frequency domain. In this study, mainly parametric identification techniques in the time domain are investigated, since it is more relevant to the identification of structural parameters related to environmental loading and/or nonlinear behavior. In Section 2, methods for modelling structural systems are discussed with an aim to system identification. In Section 3, several identification techniques based on least squares, maximum likelihood and extended Kalman filter are described. In Section 4, numerical simulation analyses are given for different structural systems, i.e., a suspension bridge, offshore structure and structure with nonlinear hysteresis. The efficiency of various identification
1989 ComputationalMechanics Publications 162
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
Fundamentals o f system identification in structural dynamics: H. Imai et al.
methods is investigated under several noise conditions. One of these simulation analyses is unique in that the method for system identification is applied to obtain an equivalent linear system representing a nonlinear offshore structure• The result obtained by this new method of identifying the equivalent linear system has been compared with those of conventional methods. Computer programs have been developed for the various system identification methods discussed in this study and used for example studies• Documentation of these programs is currently underway and will be made available in the form of a technical report in the near future• Presently, an experimental study is also being carried out using laboratory models in order to verify the validity of the identification techniques demonstrated in this study. Upon such verification, field experiments will follow• The results of these experimental studies will also be reported in the form of technical reports in the near future•
3 4
state space model ARMAX model.
In the following, it is shown that most of the structural systems can be modelled by means of ordinary differential equations (ODE) and that the remaining three models are derived from the ordinary differential equation. Example 1: Multi-story building ~9
Assuming a shear building model, an n-story building can be modeled as shown in Fig. 1. The equation of motion of this structure under ground excitation is written as: mi~ ~+ c , ( 4 , - ~ _ , ) - C,+ l ( ~ i + l
- - ~ i ) -~- k i ( ~ i - -
-k~+1(~+1-~)= -mi2 o
2. STRUCTURAL SYSTEM AND MATHEMATICAL MODEL Civil engineers deal with many types of structural dynamic systems: for example, multi-story buildings, suspension bridges, nuclear power plants, offshore structures, etc. The dynamic characteristics of these structures can be described by mathematical models. A variety of models have been developed for different purposes. Models commonly used in structural and systems engineering are the following:
~i-1)
(i= 1,2 . . . . . n)
(1)
where ~, 4i and ~'i are the relative horizontal displacement, velocity and acceleration of the ith floor to the ground, 20 is the horizontal ground acceleration, and m~, % ki are the mass, damping, stiffness coefficients, respectively• In vector-matrix notation, equation (1) can be written as follows: M 2 + C i + Kz = Lu
(2)
where z={~l,~: . . . . . (,}T and u = - 2 o. M is an n x n diagonal matrix with Mii = ml. L is an n x 1 matrix with L~= ml. C and K are both n x n symmetric matrices n
1 2
ordinary differential equation (ODE) transfer function C=
~n
m
Cl -}-£2
--C2
0
- - C2
C 2 "Jr-C 3
-- C3
0
- c3
c 3 + c4
......
0
. . . . . .
0
......
0
Ill . . . . . . . . . . . .
0
0
0
C N _ I -~-C N
....
--CN_
cN- 1
1
cN
(3a)
~n-i ~n-2
K= i
"k1 +k 2
-k 2
0
......
0
-k 2
k2+k 3
-k 3
......
0
0
- k3
k 3 q- k 4
......
0
. . . . . . . . . . . . 0
0
0
kN_l + k N ....
"
--kN_l
kN- 1
kN
(3b) 2
•
Example 2: Suspension bridge 2°
To investigate the dynamic characteristics of a suspension bridge such as aerodynamic instability, an idealized two-dimensional model of the bridge section as shown in Fig. 2 is frequently employed. The equation of motion of the model can be expressed as M ~ + C2 + Kz = L f
x
Fig. 1.
g
Shear beam structure model o f multi-story building
(4)
where z={h,a} z, f = { u , v } r with a=pitching motion, h = heaving motion, and u and v fluctuating components of the wind velocity in the horizontal and vertical directions, respectively. M is the mass matrix. C and K are
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
163
Fundamentals o1" system identification in structural dynamics: H. lmai et al.
/
where
a{t
wind
7 +
/
v(t)~/~. U
+ u(t)
O : IO- 11 (r)
The measurement y of the output z may be subjected to a measurement noise v, i.e.
"///////////////////////////////////////////////~
y=z+v=Hx+v
Fig. 2. Bridge deck model the damping and stiffness matrices aerodynamic effects.
including
the
Example 3: Offshore structure 25,1s An idealized model for the dynamic analysis of an offshore structure is shown in Fig. 2. The equation of motion of the structure is written by
where H = [I, 0]. The set of equations (9) and (10) is called the state-space model in continuous form. Usually, ~ and v are assumed to be white noise vectors with zero mean and variances Q and R, respectively. In recent years, most signal processing is accomplished by digital computers which cannot handle continuous time signals. Therefore, the output signal has to be in discrete form ( t = 0 , At, 2At . . . . ) where At denotes the sampling interval. Thus, for a linear system, the state space model is discretized as follows 2°
(Mo + C~t)2 + C2 + K o z : C~t~ + C o ( + - ~ ) l + - ~ I (5) where z is the vector of horizontal displacement, ?, ~ are the vector of horizontal wave particle velocity and acceleration, Mo, C, K o are the matrices of structural mass, damping, and stiffness, and C M and C o are the diagonal matrices containing, respectively, the inertia and drag coefficients associated with the wave force acting on the structure. In the case of linear structural systems as in equation (2), the equation of motion can also be represented by using the transfer function as Z(s) = G(s)U(s)
(6a)
G(s) = ( M s 2 + C s + K)- 1L
(6b)
where Z(s), U(s) are the Laplace transforms of z(t) and u(t), and G(s) is the transfer function of the system. The state space model in continuous form can be derived for linear and nonlinear systems. For example, by defining the state vector as
X: ~XI'~ {.X2) : {:}
(10)
x(i + 1 )= Fx(i) + Gu(i) + w(i)
(1 l a)
y(i)=Hx(i)+v(i)
(1 lb)
where F = e Am
G:
e At dt B
and
x(i)-x(iAt),
y(i)=-y(iAt)
u(i)-~u(iAt),
v(i)~v(iAt)
where it is assumed that the input u(t) is piecewise constant within each sampling interval, i.e.
u(t)=u(iAt),
iAt<=t<(i+l)At
(12)
The system noise w(i) which is defined by
(7)
w(i)=
~
(i+ 1)At CA[(/+ I)AT q D ~ ( z )
¢
dr
(13)
d i At
Equation 2 can be transformed into a state equation as
Xl}
{ f ( X l , X2 ' U)} +
is a white noise sequence with zero mean and the covariance matrix
w
x2
Qw=E[w(i)wr(i)] = x2
)
~0A!eA~DQDT[eA~]T dr
(14)
F 0 q The state space model may also be represented in innovation form 15, i.e.
where w and ~ are system noises which account for the unmeasurable input disturbances and errors in modelling. If the system is linear as in equation (2), equation (8) can be rewritten as
~=Ax+Bu+D{ 164
(9)
x(i+ l)=Fx(i)+Gu(i)+Fe(i)
(15a)
y(i) = Hx(i) + e(i)
(15b)
where x(i) is the estimation at t = i At, F is the so-called steady state Kalman gain and e(i) the innovation process.
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
Fundamentals of system ident~cation in structural dynamics: H. Imai et al. The innovation model is preferable over the general state-space model, since it has fewer parameters than the latter. When the system is controllable and observable, we can eliminate x(i) from equation (15) to obtain the following ARMAX model 2°'15
~F= [0(3), 0(4) . . . . . O(N)] r
(21c)
Then, equation (18) can be rewritten as
y(i)=oro(i)+~(i)
for i = 3 , 4 . . . . . N
(22a)
and further into matrix form as y(i) = F l y ( i - 1) + F2y(i--2) + G~u(i- 1 ) + G 2 u ( i - 2 ) +e(i)+J~e(i- l)+J2e(i-2 )
Y = ~F0+ It(3), 5(4) . . . . . ~(N)] r
(16)
If the measurable input u(-) is missing, the model is called an ARMA model. The mathematical models introduced above are used for different purposes. For conventional problems of the analysis and design of a structural system, ordinary differential equations and transfer functions have commonly been used for a long time. For the identification of modal quantities, the transfer function has been also widely employed. However, as for the identification of structural parameters, state space equations and ARMAX models are preferable because the identification techniques recently developed are based on sampled data in time.
3. ALGORITHMS FOR PARAMETER
(22b)
Then, from equation (22), the least squares estimate 0 which minimizes equation (19) can be obtained as OLs = [ v T v ] - ' [ v r v ]
(23)
It is noted that in general the least squares estimate is biased; i.e., OLs does not converge to the true value 0 as N approaches infinity. This is a major drawback of the least squares method. To overcome this difficulty several improved methods have been developed. Among them, an instrumental variable method 24'4 and a generalised least squares method should be mentioned. The least squares method is also available in sequential form ~s. Let 0~ denote an estimate of 0 based on the data {y(1), y(2) . . . . . y(i), u(1) . . . . . u(i)}, then O, can be estimated sequentially by
ESTIMATION There exists many different algorithms of parameter estimation. Most commonly used algorithms are least squares method, maximum likelihood method, extended Kalman filter, and variations of them, which will be briefly discussed in the following.
with
3.1. Least squares method 4'2°'~4
Pi = Pi-1
Consider a system described by an ARMAX model, equation (16), and define the equation error vector as
priori.
e(i)=e(i)+Jle(i-1)+J2e(i-2 )
-
-
O,=Oi_~ +K,(yr(i)-or(i)O~-~)
(24)
Ki=P,_ 10(i)(1 +Or(i)Pi_,O(i)) -~
(25)
P,- 10(i)( 1 + 0r(i)P~ - ~0(i))-'0r(i)P,-~
(26)
where the initial values 0 o, Po, 0(1) should be given a
(17)
3.2. Maximum likelihood method 1'4'2°'14 Then equation (16) can be rewritten as y(i) = F l y ( i - 1)+Fzy(i-2)+G~u(i- 1)
+ G2u(i- 2) + e(i)
(18)
Let P(y[u, 0) denote the conditional probability density, where u = [u(1 ), u(2) . . . . . u(N)] r. For given data [Y, u], one may regard P(ylu, 0) as a function of 0, which is called the likelihood function. An estimate which maximizes the likelihood function or equivalently
The least squares method is based on minimization of N
L = log P(y[u, 0)
(27)
N
J= E ~(i)%(i)= ~ [l~(i)ll 2 i=3
is called the maximum likelihood estimate and denoted by
OML,
i=3
N
As for the linear discrete-time system, the maximization of equation (27) corresponds to the minimization al
= ~ Ily(i)-Ely(i- 1)-E2y(i-2) i=3
- G l U ( i - 1)-G2u(i-2)[] 2
(19)
e(i)er(i)
J=det
with respect to
Li
0 = [F~, F2, G1, G2] r
(20)
where N is the number of data points. To solve the minimization problem, one may define a vector 0(i) and matrices Y and W as
0(i)= [y(i-- 1) r, y(i--2) r, u ( i - 1) r, u(i--z)r] r (21a) Y = [y(3), y(4) . . . . . y(N)] r
(21b)
(28)
= 1
where e(i) is a one step ahead prediction error which is defined by e(i) = y(i) - y(i[i- 1)
(29)
and y(i[i- 1)-E{y(i)ly(1), y(2) . . . . . y ( i - 1),
u(1), u(2) . . . . . u ( i - 1), 0}
(3o)
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4 165
Fundamentals of system ident~cation in structural dynamics: H. lmai et al. The prediction error e(i) can be computed sequentially according to e(1)=y(l)
P(k +
l lk+ l) = [I- K(k + 1)H]P(k + 1Ik)[l- K(k + 1)H] r + K(k + 1)Vt,K(k + 1)T
(38)
where K(k + 1) is the Kalman gain matrix which is defined as
e(2)=y(2)-Fly(l)-Gxu(1)-J~e(1)
e(i)=y(i)-Fly(i)-F2y(i-2)-Gxu(i- 1 ) - G 2 u ( i - 2 ) -Jle(i- 1)-J2e(i-2) i=3,4 ..... N (31) The problem of minimization of equation (28) may not be solved analytically in general. Thus nonlinear programming algorithms have to be used, e.g., Newton-Raphson method, Davison method 3, etc. While the maximum likelihood estimate has a nice property that it is consistent and asymptotically efficient, it usually requires a substantial amount of computational time. To circumvent this disadvantage, several approximate methods are proposed, e.g., recursive maximum likelihood method s , approximate maximum likelihood method 6, etc. Most of them are given in sequential form.
K(k+ 1 ) = P ( k +
1]k)HT[HP(k+ llk)HT+V,,] - ' (39)
The extended Kalman filtering technique described above can be applied to the identification of system parameters as follows. Consider a dynamic system described by a state equation as in equation (8), /~= f(x, u, 0; t) + w(t)
(8a)
where x is the state (response) of the system, u is the input excitation and 0 denotes system parameters. Then, by defining a new augmented state vector as
x={;}
3.3. Extended Kalman filter 6'25"14
,40,
The basic algorithm of the extended Kalman filter is a recursive process for estimating the optimal state of a nonlinear system based on observed data for the input (excitation) and output (response). It can be summarized as follows. Consider a general continuous state equation described by
a new nonlinear state equation can be constructed from equation (8a) as
:~(t) = g(X, u; t) + w(t)
(32)
The observation equation corresponding to equation (10) can be obtained as
(33)
y(k)=[H,O]{o(k)}+v(k)
X--{f(X'OoU;t)}+{W(ot) }
(41)
with observations at t = k At x(k)
y(k) = HX(k) + v(k)
in which X(k)=state vector at t=k At, y(k)=observational vector at t = k At, v(k) = observational noise vector with covariance of V~, w(t)= systems noise vector with covariance of Vw and H--matrix associated with observations. The predicted state X(k + 1[k) and its error covariance matrix P(k +l[k) can be evaluated as X(k +llk) = E{X(k + 1)[y(1), y(2) . . . . . y(k)} =X(klk)+
~
(42)
By applying the extended Kalman filter to equations (41) and (42), system parameter 0 can be estimated recursively as part of the state vector. The estimates obtained by the extended Kalman filter method may be in general biased or divergent 22'23. Some modified algorithms are proposed to improve the convergence property 13'2t'8. In this study, the weighted global iteration procedure proposed by Hoshiya and Saito 8 is also used in example studies.
( k + 1)At
g(X(tlk),u;t)dt
(34)
d k At
P(k+ l l k ) = O ( k +
1,k)P(klk)OT(k+ 1,k)+V w (35)
where E{AIB } is the expected value of A conditional to B and O(k + 1,k) is the state transition matrix which can be approximately obtained as
O(k+ l,k)~-I + At [[-0gi(x(t)' ~ t) lxl,):x/klk)
(36)
for small At. Then, the filtered state X ( k + l [ k + l ) and its error covariance matrix P(k + l lk +1) can be estimated as X(k+ ilk+ 1)= E{X(k + 1)ly(1 ), y(2) . . . . . y(k+ 1)} = X ( k + llk)+ K(k + l) x [y(k+ 1 ) - H X ( k + llk)]
166
(37)
4. N U M E R I C A L EXAMPLES AND DISCUSSION In order to investigate the accuracy and efficiency of the various system identification techniques discussed in the preceding sections, example analyses have been carried out for four different cases: (i) an idealized suspension bridge model for wind loading, (ii) a simplified model of an offshore structure for wave forces, (iii) a singledegree-of-freedom structure with bilinear hysteric characteristics under seismic excitations, and (iv) an equivalent linear system for an offshore structure model. In each case, simulated data for input and output time histories are utilized for parameter identification purposes. The estimated system parameters and responses are compared with exact ones which have been assumed a prioiri. In the last example, the system identification method is applied for determination of the equivalent linear system of an offshore structure model.
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
Fundamentals of system identification in structural dynamics: H. lmai et al. 4.1. Suspension bridge model 2° The idealized structural model shown in Fig. 2 is uscd. The equation of motion for wind loading is given in equation (4). Assumcd values of thc systcm parameters are shown in Table 1. Simulated time histories of the wind velocity fluctuations, u(t) and v(t), are shown in Fig. 4. Thc (simulatcd) observation timc histories of the bridgc motion, Yh(t) and Y~(t) are shown in Figs 5 and 6. The level of observational noise included in Yh(t) and Y~(t) is assumed to be 10% of the structural response in the root mean square (RMS) value in both cases. For the purpose of system identification, the ARMAX model has been utilized in this example. Parameter estimation is carried out by using the least square (LS) method, instrumental variable (IV) method and maximum likelihood (ML) method. The estimated parameters arc compared with (assumed) exact ones in Table 1. The response time histories rccomputed using estimated parameters, ~(t) and ~(t), are also compared with those of the observation in Figs 5 and 6. From the results in Table 1 and Figs 5 and 6, one can scc that the IV and ML methods give good estimations of the parameters as well as structural responses. On thc other hand, the results from thc LS method are found to be unsatisfactory. The poor rcsults of the LS estimates are due to the fact that the observational
5.00 +00
(~/~)
u(t)
0.00 +00
~
-5.00 +00
-
-
.
-
-
.
i
.
-
-
.
-
-
-
-
.
-
~
L
~
5.00 *00
(~/8~¢) 0.00
*OO
-S.OO +OO
.--i.-.---J-___ 0.00 25.50
50.00
75.00
TIME
(SEC)
IC0.00
125.00
Fig. 4. Simulated time histories of wind velocity fluctuations (mean horizontal wind speed = 10 m/sec)
II MASS 1 x
M.SL.
Z.O0 -Of
~7
P MASS 2
ty
F3---->
t MASS 3
r~--i->
|
% '
I MASS
MASS
(m)
o.o
-2.00
.oo
zhC~)
v,,wvv.vvvvv.vw,
VV/
-01
2.00 -01
(m)
I l MASS 6
v
~(t)LS
0.00 *00
v "
~
V -V
-v-
" "V%,,-'
--V xs,,'-,
/\ "/// X "///////
.',//////
"//////////
a) OFFSHORE TOWER
Fig. 3.
I
-2.00 -01
I
T
p
2.00 -0t
(m)
~(t)iv
b) M O D E L
Fixed offshore tower and model
0.~8 -0C
' vvvvvv vvVVV'vV vwvv -VVV/
Table 1. Exact and estimated parameters of suspension bridge P ....
ters
Exact Values LS
ML
M-tK
~1
M-'C
M-'L
{
oz
3.70 -0.12
0.95 18.25
0,128 -0.020
0.050 0.060
0.0020 0.0003
-0.057 0.012
0,0309 0.0022
5.50 0.14
-0.75 19.40
3.194 -0.013
-0.563 0.606
0.0037 0.0003
-0.122 -0.013
0.0160 0.0009
3,65 0.13
1.36 18,31
0.128 ,.0,020
0.148 0.058
0,0015 0,0004
-0.053 0.012
0.0300 0.0023
3.72 -0.11
-1,62 18.19
0.122 -0.020
1.325 0.080
0.0019 0.0011
-0.062 0.012
0.0327 0.0021
Note : i. Observational noise is assumed a.s I0 % in R M S value 2. o. = standard deviation of Zl and z2 3. U n i t : [ M - t K ] = 1/see; [ M - ' C ] = 1/see; I M - ' L ] = llsecCfirst row), X/m/see(second row); [ a , , ] = m and [t%~ I = rad/sec.
-2.0C -0: 2.00 -01
(m)
E(t)~u "
0 . 0 0 *00
-2.00 -01 ~ O.Ol
25.00
50.00
' 75.20
IOO.O0
125.00
TIME (sec)
Fig. 5. Time histories of observation Yh(t) and estimated responses h(t) for heaving motion of suspension bridge Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
167
Fundamentals of" system identification in structural dynamics: H. lmai et al. I.O0 -02
It is assumed that time histories of the displacement at the two discrete masses are available. Identification of the system parameters is carried out by using the extended Kalman filtering with the weighted global iteration procedure s,9. Table 2 shows the (assumed) exact and estimated values of the parameters for two sea states corresponding to two wind speeds of 25 and 75 ft/sec. Two different conditions of the observational noise (5 and 10% in the RMS values) are considered. The estimated values have been obtained through five global iterations. Figure 7 shows input time histories of wave particle velocities and
(t)
(r~a) 0.00 -00
-[.O0 -02 1.00 -02
(tad)
~(t)us
O.O0 -00
,,
.,J..A,h,~lL,,,h^..,, ..i..d..,h.
,~,~,"q. ] r ',, ~ " " H " ~ " " ~ ) W .
,ll,,Aill,. ,,,~)~O ,'qlr~'" I....
Table 2.
Exact and estimated parameters o/" +~[hhore structure model
~
l't/sec 2
-1.00 -02
JLt Js:
J
(~,,)
~.00 -02
(tad)
~(t)£7
D (~)
0.187 -0.022
0.4 0.0
J~2
-0.1u o.o
JT~'
0.150
0.4
Dtt
@060
0.1
@060
0.1
Ktl
3.7,50
5.0
K21
-0.750
0.0
D,,
K
(.-fi~) x,~
0.CC -,3G
L
-3.7so
0.226 .0.008 -0.147 0.127
oo
K.;
2.500
4.0
Ltl L~2
0,400 0.500
0.2
75
25
0.3
0.166 .0.005 -0.073 0.132
0.049 @064
0.261 0.028 -0.168 0.103 0.037 0.071
0.062 O.OSg
0.150 0.005 -0.043 0,125 O.Oro~ 0.058
3.743 -0.721 -3.780 2.465
3.717 -0.705 -3.779 2.449
3.779 -0,753 -3,772 2.503
3.794 -0.765 -3.765 2.520
0.408 0.494
0.419 0.486
0.422 0.405
0,446 0,492
5
10
5
10
Observatlon~ Noise Level
(% of response in R.M.S,)
OC -'3Z
- i.
1.00
-02
(~ad)
Note : #o = initial guesses of parameters 0, = Estimates by the E x ~ n d e d Kalman Filter after Fifth Global Iteration
~(t)ML
(a) Node t 0.00 -00
~0,00
(ftl~¢) O.CO
-~.00 -02 ~ 0.00
25.00
, ZO.O0
75.00
a___ lO0.CO
125.00 -~O.CO
TIME (SEC)
I0.00
Fig. 6. Time histories of observation Y,(t) and estimated responses ~(t) .for pitching motion of suspension bridge
(#/~) O.CO
vector qt(i) in equation (21a) is not statistically independent of the moving average noise term r(i) in equation (17). Further numerical investigation indicates that the results from the ML method are less sensitive to observational noise than those from the IV method.
~
I
r
I
20.00
40.00
50.00
'
f
I
90.00 time (see)
I00.00
(b) Node 2 '.O.CO i
'
(#/se~2) |
,
'
~
'
J
h (t)
'
~
'
.
--
4.2. Offshore structure 25.1s In this example, an offshore structure idealized as a two-degrees-of-freedom system is utilized. The equation of motion for wave loading, equation (5), is simplified as -[0.00
i + J~- D{(¢-i)l¢-
~1} + K z = L~
(43)
where J = [ M o + C M ] ~C, D=[Mo+CM]-ICD, K= [Mo+CM] 1Ko and L = [ M 0 + C M ] - I C M . For system identification, equation (43) is transformed into a state-space model by using the augmented state vector defined as {X}
=
{Z l Z2 Z1 -72 J l l J 2 l J12 J22 D l l D2z
Kll K21 K12 K22 LI~ 168
L22} T
(44)
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
IO.O0
I
~
'
'
I
[
,
,
,
L
,
!
'
)
,
,
~
,
4
/ O.CO ~
-I0.00
1
0.00
.
,
)
)
)
20.00
~.0.O0
60.00
~
)
,
80.00 ICO.OC time (see)
Fig. 7. Simulated time histories o,1"wave particle velocities and accelerations (wind speed = 75 ft/sec)
Fundamentals of system identification in structural dynamics: H. Imai et al. accelerations. Figure 8 shows the exact response z l(t ) and z2(t ) obtained using the assumed (exact) parameters and simulated response observation, }I1(t) and Y2(t), as well as the estimated responses, ~(t) and ~2(t), computed based on the identified parameters. Figure 9 demonstrates the convergence of recursive estimation of system parameters by the extended Kalman filtering algorithm. The results in Table 2 indicate that the extended Kalman filtering technique yields fairly good estimates even under relatively severe nonlinear hydrodynamic
(a) First Iteration
(~/s ~t) O. 30
' I
'
2g(k/k)
,
'
~
~
,
0.!5
O. O0
I
r
I
(a) Node 1 7.50
,
~
(it)
,
,
,
:
,
,
,
,
~ O. 00
20. O0
40. O0
60.00
0.00
80. OO
I00.00
time (~ec)
(b) Fifth Iteration -7.50
I
7. m.O
'
'
I
0.30
'
I
'
I
0//~)
,
i
f
fit)
D~1(t)
=
2g(~/k)
i , i
0.15 0.00 I
0.00
O/.ed)
-7.50 7. ~ 0
"
,
I
(/t)
,
I0.00
~
i
,
I
'
l
,
I
t)=2~3i~tk)
5.00 0.00
0.00
-7.50 0.00
•
20.00
40.00
60.00
80.00 time
(Set~)
,
,
,
,
,
l
,
~
,
~
,
,
,
,
,
,
(/t) O. 00
-5.
O0
r
5.00
,
~
(/t ) 0.00
-5.00
I
(It)
t
'
I
'
I
i
'
i
"
'
I
,
I
I
'
i
'
i
i
i 20.00
60.00
80.00
I00,00
time (see)
Fig. 9. Estimated parameters of offshore structure as augmented state variable (wind speed= 75 ft/sec; noise level= 10% of response in RMS)
+ 2~co~+ (oZg(ze. cO= - 2o(t)
~
i 40.00
|
80.00
A single-degree-of-freedom structure which exhibits bilinear hysteretic behavior under seismic excitation is identified. The equation of motion is written as
0.00
-5.00 0.00
!
60.00
4.3. Structure with bilinear hysteresis9
i
5.00
i 40. O0
loading conditions for both of the observational noise conditions. The initial estimates of the parameters were chosen to be quite far from the exact values. The error covariance matrix of the initial estimate has been taken fairly arbitrarily as a diagonal matrix with large values of the diagonal elements. However, it has been found that the estimated values converge to reasonable ones as the time step increases and as the iteration proceeds (Fig. 9). The estimated response obtained by using the identified parameters is found to be virtually identical to the exact response (Fig. 8).
(b) Nod~ 2 O0
20. O0
ICO.OC
Fig. 8-a. Exact, observed and estimated responses of offshore structure (wind speed= 75 ft/sec; noise level= 10% of response in RMS)
5.
' O. O0
I00.00
time (see)
Fig. 8-b. Exact, observed and estimated responses of offshore structure (wind speed= 75 fi/sec; noise level= 10% of response in RMS)
(45)
where co and { are the natural frequency and damping ratio, respectively and g(Ze,C0 defines the bilinear hysteresis with z e being the yielding displacement and c( being the ratio of post-yielding stiffness to pre-yielding stiffness as shown in Fig. 10. For the purpose of system identification, equation (45) is rewritten into a nonlinear state equation by using the state vector defined as {X} = {z ~ ~ co z e ~}r
(46)
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4 169
Fundamentals of system identification in structural dynamics: H. lmai et al.
g(Ze,a)
The results in Table 3 as well as Figs 12 14 indicate that extended Kalman filtering, particularly with the weighted global iteration procedure, estimates the system parameters and hysteretic behavior remarkably well even for the case of very severe material nonlinearity as shown in Fig. 10.
D
Fig. 10.
z
4.4. Identification of equivalent linear system 18 In this example, equivalent linearization by means of a system identification technique is discussed. In general, the nonlinear characteristics of actual structures under severe environmental loading, such as earthquakes, wind and waves are very difficult to model and analyze. Hence approximations in modelling are inevitable. For instance,
Bilinear hysteresis
400.00
~,(:)
(,~/,.,=)
f
I
T
C.OO
wlJ
"
O.nO 4 -!5.C,.0
I .=.30
0.00
I 10.SO
; 1.=.00
2O.C0
t~.e (s~) t5.CO
-40C.O0
o. oo
'=.oo
to. oc
15.00
2O.C0
Fig. II. Input ground acceleration (1940 El Centro earthquake, N S component)
O.CO
Table 3. Exact and estimated parameters of structure with bilinear hysteresis Parameters
-~5.00 0.00
Estimates
Exact
Initial
Va|ues
G uesse~
First Iteration
Fifth Iteration
o.1
o.s
0.096
0.098
w(rad/sec)
3.14
1.0
3.157
3.147
z,(cm)
3.0
t.o
2.979
3.025
a
0.1
0.5
0.090
0.088
Fig. 12. structure
5.00
Observed
50.00
'l
tO.O0
15.00
20.C0
and estimated displacement
I"
of
T
(c=Is~) Identification of parameters is carried out using the extended Kalman filtering technique with the weighted global iteration procedure. As an input ground acceleration, the 1940 E1 Centro earthquake record (N-S component, Fig. 11 ) is utilized. It is assumed that time histories of the structural displacement and velocity are available (Figs 12 and 13). The observational noise levels are taken as 10% of the structural response in the RMS values. Table 3 shows the exact and estimated values of the system parameters. Figures 12 and 13 also show time histories of the estimated responses obtained after the fifth global iteration. Finally, Fig. 14 compares four different hysteretic characteristics. They are (a) true, (b) that obtained from observational data without the extended Kalman filtering procedure by g(Ze, 0~)-
( - j~(t) - ~ - 2 ~o~) (1)2
-50.00 0.00
5.00
10.00
i5.00
20.00
tLme @~)
50.00
(c=/'s~)
5(t)
0.00
(47)
and the remaining two, (c) and (d), obtained by using extended Kalman filtering with weighted global iteration.
170
0,00
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
-50.00 0.00
Fig. 13.
I
1
I
5. O0
I0. O0
IS. GO
20. GO
Observed and estimated velocity of structure
Fundamentals of system identification in structural dynamics: H. lmai et al.
(a) 12.00
Q
(6) 12.00
L
0.00
0.00
-12.00 -12.00
0.00
12.00
-12.00 -12.00
D i s p l a c e m e n t (cm)
O. O0
12.00
Displacement (cm)
(c)
(d)
12.00
IZ.O0
0.00
0.00
A
-12.00 -12.50
0.00
12.00
-12.00 -t2.00
Displacement (cm)
(a) True Hysteresis (e) E s t i m a t e d by First Iteration
Fi9.14.
O. O0
12.00
Displacement (cm)
(b) Estimated by Conventional Method (d) Estimated by Fifth Iteration
True and estimated hystereses
in the case of offshore structures, the nonlinear hydrodynamic drag force is represented using drag coefficients which are determined semi-empirically based on representative values of the wave particle velocity and size, shape and surface roughness of the structural members. Because of the uncertainty and computational complexity associated with the nonlinear terms, the equivalent linearization technique has frequently been used in many analyses. The conventional method performs linearization by minimizing the error resulting from linearization and requires considerable numerical effort. In the present example, the idealized offshore structure model used in Section 4.2 is utilized. From the nonlinear equation of motion, equation (43), a linearized equation is
derived as
+ J*~ + K*z = L*~ + D*(~ - 2)
(48)
where J*, K*, L* and D* are the parameters of the linearized system. As in the previous example, time histories of the wave particle and acceleration and structural displacement are assumed to be available at each node. The measurement noises are assumed to be 5% of the structural responses in the RMS values. The linearized parameters for the hydrodynamic drag force as well as for the other system parameters are identified by using the extended Kalman filtering technique with weighted global iteration. Two sea states corresponding to two wind speeds of 25 and
Probabilistic Engineerin9 Mechanics, 1989, Vol. 4, No. 4
171
Fundamentals o/" system ident!fication in structural dynamics: H. lmai et al.
75 ft/sec are considered. Table 4 shows the assumed exact and estimated parameters. It should be noted here that the estimated values under each wind speed are the ensemble average of five sets of estimations arising from five sets of statistically identical but individually different wave particle motions, b(t) and/)(t), and observation time histories Y(t) generated by simulation. The estimated coefficient matrix D* has also been compared with the values computed by the conventional equivalent linearization procedure ~6. Figure 15 compares the observed response and estimated responses by two linearization methods. From the results in Fig. 15, it can be seen that the equivalent linear system obtained by using system identification can simulate the responses of the actual nonlinear system exceptionally well. The results in Table 4 indicate that for the case of small wave conditions where the nonlinear effect is not so significant, the linearized parameters obtained by the two linearization methods are quite different. However, as the wave condition becomes more severe and the nonlinear term becomes more important, the estimated values by two methods are found to be in good agreement. The above results indicate that the present method of equivalent linearization by means of system identification is considered to provide an efficient alternative to the conventional linearization method, in view of the fact that all the coefficient matrices are assumed to be unknown in the present method, while only the linearized drag coefficient matrix is considered unknown in the conventional method.
Table 4. Estimated parameters of equivalent linear system Parameters
J
(,,~)
(.-~-.) L D
(~) D"
( ,-,u, )
L
}
0.00
r
0.00
,
2 0 . O0
[
I
[
6 0 . O0
40. O0
80.30
O.187 -0.022
-0.111 0.150
0.295 -0.063 -0.172 0.249
0.258 0,080
J,2
0.400 0.000 0,000 0.400
3.73O -0,737 -3.766 2.488
3.944 -0.744 -4.234 2.578
K,j, K2,
3.750
0.097
K22
-3.750 2.54)0
5.0(}0 0.000 0,000 4.000
Ll L:
0.400 0.500
0.200 0.300
0.405 0.506
0,353 0.597
Dz|
0.060 0.060 0.300 0.300
o.o47(o.tm) o.o~(o.os0
o.3s'qo.,L,u) o.',so(o.33o)
z,,
D.
-0.750
D;
D;
In this paper, methods of parameter identification for structural dynamic systems are reviewed with emphasis on the identification of system parameters which vary with environmental loadings. For linear systems, the least squares, instrumental variable and maximum likelihood methods are reviewed. For nonlinear systems, a method utilizing the extended Kalman filter is discussed. Example numerical analyses are carried out for identification of the aerodynamic coefficients of a suspension bridge under wind loading, drag coefficients of an offshore structure under wind-induced wave forces and yield displacement and stiffness ratio of a structure with bilinear hysteresis subjected to seismic excitation. From the numerical results, it has been found that the instrumental variable and maximum likelihood methods provide good estimates for a linear system, while the extended Kalman filtering technique with weighted global iteration yields excellent estimates for nonlinear cases. A method for developing an equivalent linear system by means of system identification is also presented and a numerical simulation study indicates that this method offers an efficient alternative to the conventional linearization method.
(P)
-5.C0
Jtz J2z
E a t i m a t ~ Valu4 gstlmated Value W = l S ft/s~: W=TSft/sec
5. C O N C L U S I O N S
I
{
Initial Gum
* T h e values in parenthesis Were computed by the con~ntional linearisation method assuming D i and D~ are the only unknowns.
(a) Node 1 5.C0
Linear Syatem
System
Jll K
Nonlinear
lCO. O~
ACKNOWLEDGEMENT (b) Node 2 5.C0
This work was supported by the National Center for Earthquake Engineering Research under grant numbers 87-3012 and 88-3011 under the auspices of the National Science Foundation master contract number ECE86-07591.
i
fit) ,
i
O.CO
REFERENCES
1 -5.00
~ 0 . O0
t 2 0 . O0
,
r qO. O0
~
l
r
riO. O0
I
,
80. O0
t 00. OC
2
time (see) .
..... ....
g'xac
~.
~.ar.J.mar.~.d by P r e ~ e n ~ L L n e a r i z a r . l . o a
- F. , ~t t ~a ~e ~ b y C o n v e n e . £ o . a l
L£,qe~r-Lza~Lon
Fig. 15. Exact and estimated displacement of offshore structure (wind speed=75 ft/sec; noise l e v e l = 5 % of response in R M S )
172
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
3 4 5
AstrSm,K. I. and Eykhoff,P. System identification - a survey, Automatica, 1971, 7, 123 162 Chen,J.-C. and Garba, J. A. Structural damage assessmentusing a system identification technique, Proceedings of the Workshop on Structural Safety Evaluation Based on System Identification Approaches, Lambrecht, West Germany, June 29 July l, 1987, F. Vieweg& Sohn, Braunschweig, West Germany, pp. 474-492 Dahlquist,G. and Bj6rck, A. Numerical Methods, Prentice-Hall, Englewood Cliffs, New Jersey, 1974 Eykhoff, P. System Identification - Parameter and State Estimation, John Wiley & Sons, New York, 1974 Gertler, J. and Banyasz, C. A recursive (on-line) maximum likelihood identification method, IEEE Trans. on Automatic Control, 1974, AC-19, 816-820
Fundamentals o f system ident!~cation in structural dynamics: H. lmai et al. 6
7 8 9
Goodwin, G. C. and Payne, R. L. Dynamic System ldemifcation: Experimental Design and Data Analysis, Academic Press, New York, 1977 Hart, G. C. and Yao, J. T.-P. System identification in structural dynamics, Journal of Engineering Mechanics, ASCE, 1977, 106(EM6), 1089 1104 Hoshiya, M. and Saito, E. Structural identification by extended Kalman filter, Journal of Engineering Mechanics, ASCE, 1984, 110(12), 1757 1770 Hoshiya, M. and Maruyama, O. Identification of nonlinear structural systems, Proceedings of ICASP-5, Vancouver, Canada, 1987, Vol. I, pp. 182-189
10
Proceedingsof 6th IFAC Symposiumon Identification and System Parameter Estimation. Arlington, Virginia, June 1982
11
Kashyap, R. L. Maximum likelihood identification of stochastic linear systems, IEEE Trans. on Automatic Control, 1970, AC-15, 25-34 Kozin, F. and Natke, H. G. System identification techniques, Structural Safety, 1986, 3, 269 316 Ljung, L. Asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems, IEEE Trans. Automatic Control, 1979, AC-24, 36-50 Ljung, L. System Identification - Theory for the User, Prentice-Hall, Englewood Cliffs, New Jersey, 1987 Ljung, L. and S6derstr6m, T. Theory and Practice of Recursive Identification, MIT Press, Cambridge, Massachusetts, 1983 Malhotra, A. and Penzien, J. Nondeterministic analysis of offshore structures, Journal of Engineering Mechanics, ASCE, 1970, 96(EM6), 985-1003
12 13 14 15 16
17
18 19 20 21 22 23
24 25
Natke, H. G. and Yao, J. T.-P. System identification approach in structural damage evaluation, ASCE Structures Congress '86, 1986, Preprint 17-1 Paliou, C. and Shinozuka, M. Identification of equivalent linear systems, accepted for publication in the Journal of Engineering Mechanics, ASCE, 1988 Shinozuka, M., Itagaki, H. and Hakuno, M. Dynamic safety analysis of multistory buildings, Journal of Structural Engineering, ASCE, 1968, 94(ST1), 309-330 Shinozuka, M., Yun, C.-B. and Imai, H. Identification of linear structural dynamic systems, Journal of Engineering Mechanics, ASCE, 1982, 108(EM6), 1371 1390 Song, T. L. and Speyer, J. L. The modified gain extended Kalman filter and parameter identification in linear systems, Automatica, 1986, 22, 59 75 Urain, B. Asymptotic convergence properties of the extended Kalman filter using filtered state estimates, IEEE Trans. on Automatic Control, 1980, AC-25, 1207 1211 Westerlund, T. and Tyss, T. Remarks on asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems, IEEE Trans. Automatic Control, 1980, AC-25, 1011-1012 Wong, K. Y. and Polack, E. Identification of linear discrete time systems using the instrumental variable method, IEEE Trans. on Automatic Control, 1967, AC-12, 707 718 Yun, C.-B. and Shinozuka, M. Identification of nonlinear structural dynamic systems, Journal of Structural Mechanics, 1980, 8, 187-203
Probabilistic Engineering Mechanics, 1989, Vol. 4, No. 4
173