- Email: [email protected]
Parameter identification of vehicles moving on continuous bridges
ARTICLE IN PRESS
Journal of Sound and Vibration 269 (2004) 91–111
JOURNAL OF SOUND AND VIBRATION www.elsevier.com/locate/jsvi
Parameter identification of vehicles moving on continuous bridges F.T.K. Au*, R.J. Jiang, Y.K. Cheung Department of Civil Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, People’s Republic of China Received 5 June 2002; accepted 19 December 2002
Abstract This paper describes a method for the identification of parameters of vehicles moving on multi-span continuous bridges. Each moving vehicle is modelled as a 2-degree-of-freedom system that comprises four components: an unsprung mass and a sprung mass, which are connected together by a damper and a spring. The corresponding parameters of these four components, namely, the equivalent unsprung mass and sprung mass, the damping coefficient and the spring stiffness are identified based on dynamic simulation of the vehicle–bridge system. The identification process makes use of acceleration measurements at selected stations on the bridge. In the study, the acceleration measurements are simulated from the solution to the forward problem of a continuous beam under moving vehicles, together with the addition of artificially generated measurement noise. The identification is carried out through a robust multi-stage optimization scheme based on genetic algorithms, which searches for the best estimates of parameters by minimizing the errors between the measured accelerations and the reconstructed accelerations from the moving vehicles. This multi-stage optimization scheme reduces the variable search domains stage by stage using the identified results of the previous stage. Besides the basic operators in simple genetic algorithms, some advanced genetic operators and techniques are adopted here. Therefore, it makes the proposed identification procedure much more efficient than other traditional optimization methods. The identification procedure is then verified with a few test cases. The estimated vehicle parameters can also be used to get the time varying contact forces between the vehicles and bridge surface. r 2003 Elsevier Ltd. All rights reserved.
1. Introduction Parameter identification is an important kind of problem relevant to many fields, and various researchers have experimented with different optimization methods and techniques. For the *Corresponding author. Tel.: +1-2859-2650; fax: +1-2559-5337. E-mail address: [email protected] (F.T.K. Au). 0022-460X/03/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-460X(03)00005-1
ARTICLE IN PRESS 92
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
identification of vehicle parameters, one of the traditional tools is the least squares method [1–3] and its amended versions [2]. The other methods for identification of vehicle parameters include the extended Kalman filter method [4,5], the maximum likelihood method [6], and the recursive prediction error method [7]. A comparative study [8] has been carried out on the last three methods of identification of vehicle parameters, and the study has discovered that the most common problem with these optimization methods is convergence. Different models have been created to represent the actual behaviour of vehicles with different emphasis. For models of vehicle–tyre system with different numbers of degrees-of-freedom (d.o.f.s) [1,3,6,8], different numbers of system parameters are identified. In some cases, attention is paid to the driver–vehicle interaction analysis [9] and the vehicle parameters for riding quality optimization [4,10]. At other times, the aim of parameter identification is to get the tyre forces and the road friction coefficients for vehicles on different road surfaces [11]. Most of the identifications mentioned above are realized based on the simulated or experimental test data of the vehicles. From the point of view of dynamic coupling between the vehicle and the road, the motion of a bridge is largely related to the motion of a vehicle on it. This is well verified by the results of forward analysis of bridge dynamics under moving loads [12] and moving vehicles [13]. All the forward dynamic analysis of bridges under moving vehicles is based on the assumption that the vehicle parameter values are known correctly. However, for the normal operation and structural health monitoring of existing bridges open to public, it is unrealistic to expect that all the correct vehicle parameters are known beforehand. Therefore the identification of vehicle parameters is an important step in the study of dynamic behaviour of existing bridges. This paper describes a method for the identification of parameters of vehicles moving on multispan continuous bridges. Instead of just modelling the vehicles as moving masses [14], each moving vehicle is modelled as a 2-d.o.f. system [13] that comprises four components: an unsprung mass and a sprung mass, which are connected together by a damper and a spring. The corresponding parameters of these four components, namely, the equivalent unsprung mass and sprung mass, the damping coefficient and the spring stiffness are identified based on the acceleration measurements at various stations along the continuous beam. In the study, the acceleration measurements are simulated from the solution to the forward problem of a continuous beam under moving vehicles based on the modified vibration functions [12] and modal superposition, together with the addition of artificially generated measurement noise. The identification is carried out through a robust multi-stage optimization scheme based on genetic algorithms [15], which searches for the best estimates of parameters by minimizing the errors between the measured accelerations and the reconstructed accelerations from the moving vehicles. This multi-stage optimization scheme reduces the variable search domains stage by stage using the identified results of the previous stage. Based on the mechanics of natural selection and natural genetics, genetic algorithms combine survival of the fittest among string structures with a structured yet randomized information exchange to form a search algorithm with some of the innovative flair of human search. Their special search mechanics leads to their special robust function of finding the global optimal values of the variables definitely without worrying about the convergence problem nor the restrictive requirements of continuity and derivative existence, as compared with the traditional optimization techniques. Apart from the basic strategies of simple genetic algorithms, other advanced genetic techniques, such as niching, elitism, fitness scaling, are also adopted to improve the efficiency. In the present parameter identification problem, a natural
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
93
search domain that contains all the possible options for each variable does not exist, at least at the initial stage of solution. In order to overcome this problem and to improve the efficiency, a multistage optimization procedure based on genetic searching is put forward by reducing the variable search domains stage by stage using the identified results of the previous stage. A comprehensive case study is then carried out to verify the robustness of the proposed identification method. The effects of different noise levels, various parameters used in the genetic algorithms, the number of measurement stations, and the number of vehicles on the results are studied. The time-dependent contact forces between the moving vehicles and the bridge surface can also be calculated from the identified vehicle parameters. The work described in this paper will be useful to the operation and structural health monitoring of existing bridges.
2. Simulation of acceleration measurements 2.1. Modelling of vehicles and bridge In this paper, a multi-span continuous bridge is modelled as a linear elastic Bernoulli–Euler beam with totally (P+1) point supports as shown in Fig. 1. The bridge may have varying sections characterized by the density rðxÞ; cross-sectional area AðxÞ; Young’s modulus EðxÞ and moment of inertia IðxÞ at location x: The N moving vehicles on the bridge are modelled as 2-d.o.f. systems each with four parameters {Ms1, Ms2, cs, ks, s=1,2, y, N}. The typical sth moving vehicle is modelled by the unsprung mass Ms1 and the sprung mass Ms2 ; which are interconnected by a spring of stiffness ks and a dashpot of damping coefficient cs : The convoy of vehicles moves from left to right along the bridge at a speed vðtÞ: Setting the position of the leftmost point support of the bridge as the origin O; the horizontal position of the sth vehicle is defined as xs ðtÞ; a function of t. The deflection of the bridge at position x and at time t is denoted by yðx; tÞ where upward deflection is taken as positive, as shown in Fig. 1. The vertical displacements of the unsprung mass and sprung mass of sth vehicle are, respectively, represented by fys1 ðtÞ; ys2 ðtÞ; s ¼ 1; 2y; Ng; which are also measured vertically upwards with reference to their respective vertical positions before coming onto the bridge.
v(t) ys2(t)
Ms2 y
O
cs
ks Ms1
0
1
2
ys1(t) x P
xs(t) L
Fig. 1. A multi-span continuous bridge with N moving vehicles.
ARTICLE IN PRESS 94
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
If the surface roughness of the bridge and the separation of the vehicle from the bridge are neglected, the vertical displacement ys1 ðtÞ of the unsprung mass of the sth vehicle model is the same as the deflection of the beam yðxs ðtÞ; tÞ at the corresponding location xs ðtÞ: Then the vertical displacement ys1 ðtÞ and its first two derivatives are as follows: ys1 ðtÞ ¼ yðx; tÞjx¼xs ðtÞ ; dys1 ðtÞ @yðx; tÞ @yðx; tÞ ¼ þ v ; dt @t @x x¼xs ðtÞ 2 2 d2 ys1 ðtÞ @ yðx; tÞ @2 yðx; tÞ @yðx; tÞ 2 @ yðx; tÞ ¼ þ 2v þ a ; þ v dt2 @t2 @x @t @x2 @x x¼xs ðtÞ
ð1Þ ð2Þ
ð3Þ
where a is the horizontal acceleration of the group of vehicles. The contact force fcs ðtÞ between the sth unsprung mass and the beam can be expressed as fcs ðtÞ ¼ ðMs1 þ Ms2 Þg þ Dfcs ðtÞ;
ð4Þ
in which g is the acceleration due to gravity and Dfcs ðtÞ is the time-dependent variation of contact force. The equation of vertical motion for the sth unsprung mass Ms1 and sprung mass Ms2 can be written as d2 ys1 ðtÞ dys2 ðtÞ dys1 ðtÞ ½ Ms1 ¼ Df ðtÞ þ k y ðtÞ y ðtÞ þ c ; ð5Þ cs s s2 s1 s dt2 dt dt d2 ys2 ðtÞ dys2 ðtÞ dys1 ðtÞ Ms2 ¼ ks ½ys2 ðtÞ ys1 ðtÞ cs : ð6Þ dt2 dt dt From Eqs. (4)–(6), the contact force between the sth vehicle and the bridge can be written as d2 ys1 ðtÞ d2 ys2 ðtÞ þ M : ð7Þ fcs ðtÞ ¼ ðMs1 þ Ms2 Þg þ Ms1 s2 dt2 dt2 For the modelling of the continuous bridge, the modified vibration functions [12] are chosen, which are built up from the vibration modes of a hypothetical single-span simply supported beam having the same end supports and total length as the real beam, and cubic spline expressions that ensure the boundary conditions at the supports. Therefore the assumed vibration modes can be written as * i ðxÞ; % i ðxÞ þ F ð8Þ Fi ðxÞ ¼ F % i ðxÞ; i ¼ 1; 2; y; ng are the first n vibration modes of a hypothetical beam of total length where fF * i ðxÞ; i ¼ 1; 2; y; ng L with the same end supports but without the intermediate supports, and fF are the augmenting cubic spline expressions which are so chosen that each Fi ðxÞ satisfies the boundary conditions at the two ends and the zero deflection conditions at the intermediate point % i ðxÞ are the supports. For a hypothetical beam that is simply supported, the vibration modes F * Fourier sine series. The augmenting cubic spline expressions Fi ðxÞ can be determined by using the zero deflection conditions at the intermediate supports and the boundary conditions at end supports.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
95
By modal superposition and separation of variables, the vibration of the beam can be expressed as yðx; tÞ ¼
n X
qi ðtÞFi ðxÞ
ði ¼ 1; 2; y; nÞ;
ð9Þ
i¼1
where qi ðtÞ are the generalized coordinates. 2.2. Simulation of noise-free acceleration measurements Based on the vehicle and bridge models introduced above, dynamic analysis of the system can be carried out and the simulated noise-free acceleration measurements at specified stations of the bridge can be obtained for the subsequent parameter identification of the vehicles. This part is in fact a conventional forward problem. The velocity of vibration and the curvature of the beam at position x are, respectively, n @yðx; tÞ X ¼ q’ i ðtÞ Fi ðxÞ; @t i¼1
ð10Þ
n @2 yðx; tÞ X ¼ qi ðtÞF00i ðxÞ; @x2 i¼1
ð11Þ
in which the dot represents differentiation with respect to time t, and the prime represents differentiation with respect to the abscissa x. Therefore the kinetic energy T, the bending energy V, and the work done by external forces W are, respectively, Z n X n 1 L @yðx; tÞ 2 1X rðxÞAðxÞ dx ¼ ð12Þ T¼ q’ i ðtÞmij q’ j ðtÞ; 2 0 @t 2 i¼1 j¼1 1 V¼ 2
Z
L 0
W ¼
2 2 n X n @ yðx; tÞ 1X EðxÞIðxÞ dx ¼ qi ðtÞkij qj ðtÞ; @x2 2 i¼1 j¼1
Z 0
N LX
fcs yðx; tÞd½x xs ðtÞ dx ¼
s¼1
n X
qi ðtÞfi ðtÞ;
ð13Þ
ð14Þ
i¼1
where the work done by external forces W is defined in terms of the Dirac delta function ( N ðx ¼ x0 Þ ð15aÞ dðx x0 Þ ¼ 0 ðxax0 Þ Z
N
dðx x0 Þ dx ¼ 1; N
ð15bÞ
ARTICLE IN PRESS 96
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
while mij ; kij and fi ðtÞ are, respectively, the generalized mass, stiffness and force, which can be written as Z L mij ¼ rðxÞAðxÞFi ðxÞFj ðxÞ dx; ð16Þ 0
kij ¼
Z
L
EðxÞIðxÞF00i ðxÞF00j ðxÞ dx;
0
fi ðtÞ ¼
N X
ð fcs Fi ½xs ðtÞ Þ ¼
s¼1
N X
fis ðtÞ:
ð17Þ
ð18Þ
s¼1
The component fis ðtÞ in the generalized force fi ðtÞ can be further written as fis ðtÞ ¼ ðMs1 þ Ms2 ÞgFi ðxs ðtÞÞ Ms2 y. s2 ðtÞFi ðxs ðtÞÞ n X Ms1 Fi ðxs ðtÞÞ fq. j ðtÞFj ðxs ðtÞÞ þ 2vq’ j ðtÞF0j ðxs ðtÞÞ þ qj ðtÞ½v2 F00j ðxs ðtÞÞ þ aF0j ðxs ðtÞÞ g:
ð19Þ
j¼1
The Lagrangian function for the vehicle-bridge system is G ¼ T ðV W Þ; where ðV W Þ is the total potential energy, and the Euler-Lagrange equation appears as d @G @G ¼ 0: ð20Þ dt @q’ i @qi Substituting Eqs. (16–17) and (18) into Eq. (20) gives n n n X X X # ij q. j ðtÞ þ m c#ij q’ j ðtÞ þ k#ij qj ðtÞ ¼ p# i ðtÞ; j¼1
j¼1
i ¼ 1; 2; y; n;
ð21Þ
j¼1
where # ij ðtÞ ¼ mij þ m
N X
Ms1 Fi ðxs ðtÞÞFj ðxs ðtÞÞ;
ð22Þ
s¼1
c#ij ðtÞ ¼
N X
2vMs1 Fi ðxs ðtÞÞF0j ðxs ðtÞÞ;
ð23Þ
s¼1
k#ij ðtÞ ¼ kij þ
N X
h i Ms1 Fi ðxs ðtÞÞ v2 F00j ðxs ðtÞÞ þ aF0j ðxs ðtÞÞ ;
ð24Þ
s¼1 N X ðMs1 þ Ms2 ÞgFi ðxs ðtÞÞ: p# i ðtÞ ¼
ð25Þ
s¼1
Should a particular vehicle be outside the bridge, the corresponding terms under the summation signs should be omitted. Eq. (21) can be solved by the Newmark-b method. Therefore fqi ðtÞg and fq. i ðtÞg can be obtained accordingly. The simulated ‘‘noise-free’’ acceleration anf ðx; tÞ of the beam at measurement
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
97
location x can then be expressed as anf ðx; tÞ ¼ yðx; . tÞ ¼
n X
q. i ðtÞFi ðxÞ:
ð26Þ
i¼1
The vector containing the simulated ‘‘noise-free’’ acceleration time histories at the measurement stations then forms the noise-free acceleration vector Anf : 2.3. Simulation of measurement noise In practical experiments or fieldwork, the measurements are usually polluted by measurement noise. To evaluate the sensitivity of the results to such measurement noise, noise-polluted measurements are simulated by adding to the noise-free acceleration vector a corresponding noise vector whose root-mean-square (RMS) value is equal to a certain percentage of the r.m.s. value of the noise-free data vector. The components of all the noise vectors are of Gaussian distribution, uncorrelated and with a zero mean and unit standard deviation. Then the noise-polluted acceleration Anp of the bridge at location x can be simulated by Anp ¼ Anf þ RMS ðAnf Þ Nlevel Nunit ;
ð27Þ
where RMSðAnf Þ is the root-mean-square value of the noise-free acceleration vector Anf ; Nlevel is the noise level, and Nunit is the randomly generated noise vector with zero mean and unit standard deviation. In the present study, besides the case of zero noise, five levels of noise pollution will be investigated, namely 5, 10, 15, 20% and 25%.
3. Multi-stage parameter identification based on genetic algorithms 3.1. Genetic algorithms Genetic algorithms [15] are stochastic global search techniques based on the mechanics of natural selection and natural genetics. They combine survival of the fittest among string structures with a structured yet randomized information exchange to form a search algorithm with some of the innovative flair of human search. The central theme of genetic algorithms is robustness, which results from their special search ways, namely, working with a coding of the parameter set, not the parameters themselves; searching from a population of points, not a single point; using objective function information, not derivatives or other auxiliary knowledge; and using probabilistic transition rules, not deterministic rules. Because of these intrinsic characteristics, the search and optimization direction of genetic algorithms is definitely guided to global optimal values of the variables instead of the local optimal solution in traditional search and optimization methods. Convergence is not a problem for genetic algorithms at all. Genetic algorithms have found applications in a wide spectrum of areas, covering problems with discrete variables and continuous variables. In the present study, the variables to be identified are the equivalent unsprung mass and sprung mass, the damping coefficient and the spring stiffness for each vehicle. Assuming that there are all together N vehicles on a bridge, the number of parameters to be identified is therefore 4N. All the
ARTICLE IN PRESS 98
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
{M11, M12, …, M1N} {M21, M22, …, M2N} {c1, c2, …, cN} {k1, k2, …, kN} N
N
N
N
4N
Fig. 2. The structure of an individual with 4N variables in the present study.
parameters are coded as binary strings with different string lengths, where the string length is defined as the number of the binary bits for one variable. For the purpose of optimization using genetic algorithms, these strings are assembled together in prescribed order to form individuals, and the structure of an individual adopted in the present study is shown in Fig. 2. To kick off the search and optimization process, a specified number of individuals must be randomly created to form an initial population. 3.2. Genetic techniques adopted The three basic genetic operations in simple genetic algorithms are adopted here, namely reproduction of individuals, crossover between individuals in a population, and the mutation of individuals. Besides these, a few advanced genetic techniques are used, which improves the identification efficiency in each stage. After an initial population is randomly created with the specified population size, the tournament selection scheme is used with a shuffling technique for choosing random pairs for mating. In order to keep the fittest individuals in the later populations, elitism is also included, which means that the best individuals in a generation will be reproduced to the next generation. Uniform crossover is adopted for the proposed identification process. Jump mutation or creep mutation scheme is used to realize mutation. Niche theory is also introduced. A sharing function used to induce niche in this genetic algorithm limits the uncontrolled growth of particular species within a population. Another important improvement scheme—linear fitness scaling [15] is implemented to keep appropriate levels of competition throughout a simulation. Without scaling of the objective function values, there is a tendency for a few super-individuals to dominate the selection process early on during a run, which is a leading cause of premature convergence. In this case, the objective function values must be scaled down to prevent a takeover of the population by these super-strings. Later on during the run, when the population is largely converged, competition among population members is less strong and the simulation tends to wander. In this case, the objective function values must be scaled up to accentuate differences between population members to continue to reward the best performers. In both cases, at the beginning of the run and as the run matures, the linear fitness scaling can help. The fitness value is related to the objective function to reflect the desirability of the individual. Linear scaling requires a linear relationship between the scaled fitness f 0 and the raw fitness f as follows: f 0 ¼ af þ b:
ð28Þ
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
99
The coefficients a and b are chosen based on two criteria. The first is to enforce equality of the 0 of the scaled fitness values. average value favg of the raw fitness values and the average value favg The second criterion is to control the number of offspring given to the population member with maximum raw fitness. In particular, the best population member is expected to give Cmult copies in the next generation, i.e., 0 ¼ Cmult favg ; fmax
ð29Þ
0 where fmax is the maximum scaled fitness. Here Cmult is taken as 2.0, which means that the number of expected copies desired for the best population member is 2. In the present identification of vehicle parameters, the errors at the measurement stations between the measured or simulated accelerations and the reconstructed accelerations from the identified vehicle parameters should be minimized. For a bridge with Q measurement stations and a total of R time points in each measured time history, the objective function F to be maximized in the present study can be written in terms of the error vector ei as !1=2 Q Q R X X X 2 F ¼ jjei jj2 ¼ ðeij Þ ; ð30Þ i¼1
i¼1
j¼1
where jjei jj2 is the Euclidean vector norm of the error vector ei defined as ei ¼ faij;simu g faij;iden g;
ð31Þ
where aij;simu is the element of the simulated acceleration vector and aij;iden is the element of the acceleration vector calculated from the trial identified vehicle parameters in each generation, both at the ith measurement station and the jth time point. Linear scaling works well except when negative fitness calculation prevents its use [15]. To overcome this problem associated with Eq. (30), a mapping from the objective function F to the raw fitness f is carried out as f ¼ F þ Cmin ;
ð32Þ
where Cmin is chosen to enforce that the raw fitness f is always positive. In the present study, the value of Cmin is taken to be the absolute value of the worst objective function value in each current generation, which prevents the occurrence of the negative fitness in every generation. 3.3. Multi-stage optimization In the present parameter identification problem, a natural search domain that contains all the possible options for each variable does not exist obviously, at least at the initial stage of solution. Although the identification of the contact forces provides useful clues to the identification of moving masses [14], this may not be as useful here. In order not to miss the global optimal set of parameters and converge on a wrong solution, the initial search domain for each variable must be set sufficiently large. However within these very large search domains, it is often difficult to find the ideal global optimal values by a single optimization run, and it also takes an extremely long running time. Without any idea of the likely magnitudes of the optimal parameters, it is difficult to determine the appropriate string lengths, which are essential for accurate results. In order to overcome this problem, the parameter identification is carried out through a multi-stage
ARTICLE IN PRESS 100
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
optimization scheme involving the repeated use of genetic algorithms and gradual refinement of search domains as well as accuracy. Using reasonable string lengths and a relatively small number of generations, it is possible for the initial stage of optimization to produce rough estimates of the parameters. Such estimates are then used to map out the refined search domains for the second stage. It is known that the genetic optimization gradually guides the trial solutions to the optimal one. Therefore through this staged optimization process, the search domains will be reducing in each successive stage in the light of the identification results of the previous stage until the desired identification precision is reached. In the present study, all parameters to be identified are non-negative values, and so the lower limits of their initial search domains are all set as zero. The upper limits of the initial search domains for the spring stiffness and the damping coefficient are set to 1020 or even larger values in standard units. For both the unsprung and sprung masses for normal bridges, a convenient upper limit of the initial search domain can be taken as 106 kg. These search domains can then be narrowed down in subsequent stages of optimization. For example, assuming the identified unsprung mass of the sth vehicle in a stage to be Ms1;iden ; the search domain D of this unsprung mass for the next stage can be set as D ¼ ðDl ; Du Þ ¼ Ms1;iden ðCl ; Cu Þ;
0oCl o1
1oCu ;
ð33Þ
where Cl and Cu are, respectively, coefficients to define the lower limit Dl and the upper limit Du of the search domain D. The same strategy can be applied to the search for parameters for the sprung mass Ms2 ; the spring stiffness ks and the damping coefficient cs : The values of Cl and Cu may be varied at different stages for improved efficiency. The termination of the search at each stage is controlled both by the changes of parameters between successive generations as well as the maximum number of generations allowed. However the termination of the whole optimization procedure is controlled by the differences between the identified parameters of consecutive stages. The efficiency of this multi-stage optimization process based on genetic algorithms at each stage will be verified by a comprehensive case study.
4. Case study 4.1. General description To evaluate the efficiency and effectiveness of the proposed method, a comprehensive case study covering various aspects have been carried out. The effects of population sizes, number of measurement stations, levels of measurement noise, and number of vehicles on the accuracy are examined. Table 1 summarizes the different aspects examined in various cases. The contact forces between the moving vehicles and the bridge surface are also estimated using the identified vehicle parameters. The identification error Erroriden of, for example, the unsprung mass of sth vehicle is defined as Erroriden ¼
jMs1;iden Ms1;true j
100%; Ms1;true
ð34Þ
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
101
Table 1 Cases for investigation Case
No. of vehicles
No. of stations
Noise level (%)
Population size
Effects to investigate
1 2 3 4 5
1 1 1 1 3
1 1 1 1, 2, 3, 4, 6 2, 6
0 5 5, 10, 15, 20, 25 10 10
30 5, 10, 20, 30, 40, 60, 80, 100 30 30 60
Basic strategy Population size Measurement noise Number of stations Number of vehicles
36 m
48 m
36 m
L = 120 m (a) 4m 0.2 m
0.4 m
0.4 m
3.0 m 0.2 m
0.75 m
2.5 m
0.75 m
(b)
Fig. 3. A three-span continuous bridge. (a) Idealized beam model; (b) idealized cross-section.
where Ms1;iden and Ms1;true are the identified and true values of unsprung mass Ms1 ; respectively. A similar definition is also applied to other parameters such as the sprung mass, spring stiffness and damping coefficient. A three-span continuous box girder bridge is used in the case study, and the idealized beam model and cross-section are shown in Fig. 3. The density r and Young’s modulus E of concrete are 2400 kg/m3 and 30 000 Mpa, respectively. The span lengths from left to right are 36, 48, 36 m, respectively. The cross-sectional area A and the second moment of inertia I are 3.46 and 3.854 m4, respectively. For both the acceleration simulation and parameter identification, the first three
ARTICLE IN PRESS 102
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
vibration modes of the bridge are adopted. All the vehicles are assumed to travel at a horizontal speed v of 17 m/s across the bridge. Both for the simulation of acceleration measurements and the subsequent analysis, the total traversing time T, which is counted from the instant when the first vehicle arrives at the starting end of the bridge to that when the last vehicle leaves the bridge, is divided into 200 time intervals. The histories of acceleration and contact forces are therefore presented accordingly with respect to time ratio t/T. The contact forces between the vehicles and the bridge surface are calculated from fcs ðtÞ ¼ ðMs1 þ Ms2 ÞgFi ðxs ðtÞÞ Ms2 y. s2 ðtÞFi ðxs ðtÞÞ n X Ms1 Fi ðxs ðtÞÞ fq. j ðtÞFj ðxs ðtÞÞ þ 2vq’ j ðtÞF0j ðxs ðtÞÞ þ qj ðtÞ½v2 F00j ðxs ðtÞÞ þ aF0j ðxs ðtÞÞ g:
ð35Þ
j¼1
The parameters for the present genetic algorithm with uniform crossover are as follows: crossover probability Pc ¼ 0:5; creep mutation probability Pcm ¼ 0:04; and jump mutation probability Pjm ¼ 1=Npop ; where Npop denotes the population size. As seen in the following cases, where the signals are not contaminated with measurement noise, essentially the true parameters are obtained. 4.2. Case study 1: basic strategy In Case 1, there is only one moving vehicle with simulated noise-free accelerations at the measurement station located at the middle of the central span, i.e., x ¼ 60 m: This is mainly to establish the validity and effectiveness of this multi-stage optimization procedure. The true parameters for the vehicle used are the unsprung mass M1=4800 kg, the sprung mass M2=27000 kg, the damping coefficient c=86000 Ns/m, and the spring stiffness k=9.12 106 N/m. These parameters are applicable to the other cases unless otherwise stated. In the first optimization stage, the initial search domains of the four vehicle parameters are set as M1 Að0; 10 000Þ kg; M2 Að0; 60 000Þ kg; cAð0; 200 000Þ Ns/m and kAð0; 2 107 Þ N=m; and their string lengths are 7, 8, 8 and 8, respectively. The population size is set as 30. By generation 300 of the first stage, the identified parameters M1, M2, c and k are 4961, 26 824 kg, 89 412 Ns/m and 9.098 106 N/m, respectively. Based on these identified values, and taking the coefficients Cl and Cu for Eq. (33) as 0.9 and 1.1, respectively, the search domains are narrowed down to M1A(4465, 5457) kg, M2A(24 142, 29 506) kg, cA(80 471, 98 353) Ns/m and kA(8.188 106, 10.0078 106) N/ m. Similarly, the genetic search is carried out again within these refined search domains. The identified parameters in this stage are then used to map out the search domains for the next stage. In this manner, the correct results are gradually identified as shown in Table 2. Only three stages are necessary to get the essentially true values of the parameters. The redefinition of search domains for successive stages enables the identification process to be both accurate and efficient. Table 3 shows the evolution of results for the first stage. The maximum number of generations for this stage is set as 800. It can be seen that among the first 90 generations, the results change a lot, while the results in generations 244–800 are essentially the same. Therefore the chosen number of generations is enough for the first stage, and it can be safely terminated.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
103
Table 2 Case 1: summary of results at different stages noise-free measurements, station at x=60 m, population size=30 Stage
1 2 3
Limit coefficients for search domain
Identified parameters
Cl
Cu
M1 (kg)
M2 (kg)
c (Ns/m)
k(N/m)
— 0.9 0.95
— 1.1 1.05
4961 4729 4802
26 824 27 292 27 002
89 412 88 167 86 038
9 098 039 9 224 698 9 120 620
4800 0.04
27 000 0.01
86 000 0.04
9 120 000 0.01
True value Error (%)
Table 3 Case 1: evolution of results in the first stage (noise-free measurements, station at x=60 m, population size=30) Generation number
Identified parameters M1 (kg)
M2 (kg)
c (Ns/m)
k (N/m)
1–5 6–7 8 9–11 12 13–17 18 19–24 25–26 27–28 29 30–31 32–33 34–39 40 41–54 55 56–69 70–73 74–88 89–243 244–800
7874 2992 6457 2835 4094 4094 3780 3780 4409 3937 4724 3937 3780 3937 4252 4094 4094 4330 4252 4252 4252 4961
25647 29412 25647 29882 29412 29412 29412 29412 29176 28941 28470 28470 28470 28470 28470 28471 28471 28235 28235 28235 28235 26824
71373 144314 69804 96470 94902 96471 97255 99608 99608 99608 99608 99608 99608 99608 99608 99608 98039 98039 98039 94902 96471 89412
9 254 902 9 725 490 9 333 333 9 725 490 9 882 353 9 725 490 9 725 490 9 725 490 9 725 490 9 647 059 9 568 627 9 568 627 9 411 765 9 490 196 9 490 196 9 490 196 9 490 196 9 490 196 9 411 765 9 411 765 9 411 765 9 098 039
Error (%)
3.35
0.65
3.97
0.24
4.3. Case study 2: effect of population size This is essentially the same as Case 1 except that 5% measurement noise is introduced and that the population size is varied to study its effect on the accuracy and efficiency of the method.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
104
Table 4 Case 2: effects of population size on accuracy. (5% measurement noise, station at x=60 m) Population size
5 10 20 30 40 60 80 100
Number of stages
4 3 3 2 2 2 2 2
Identified values
Errors (%)
M1 (kg)
M2 (kg)
c (Ns/m)
k (N/m)
M1
M2
c
k
4778 4758 4790 4775 4774 4775 4792 4780
27052 27035 26942 27042 27035 26957 26963 26963
86 728 86 147 85 267 85 668 85 806 86 229 85 845 86 192
9 118 394 9 095 260 9 106 979 9 107 513 9 118 959 9 109 327 9 141 072 9 113 324
0.46 0.88 0.21 0.52 0.55 0.52 0.17 0.42
0.19 0.13 0.22 0.15 0.13 0.15 0.14 0.14
0.85 0.17 0.85 0.39 0.23 0.27 0.18 0.22
0.02 0.27 0.14 0.14 0.01 0.12 0.23 0.07
325
Contact force (kN)
320
315
310
305
300
0
0.2
0.4 0.6 Time ratio t / T
0.8
1
Fig. 4. Case 2: contact forces calculated from the true and identified vehicle parameters using different population sizes: —— true; – – – size = 5; — — size=30; — – — size=60.
Starting from the initial search domains M1A(0, 10 000) kg, M2A(0, 60 000) kg, cA(0, 200 000) Ns/m and kA(0.2 107) N/m, the proposed identification method is applied with eight different kinds of population sizes, namely 5, 10, 20, 30, 40, 60, 80 and 100. The results are summarized in Table 4. First and foremost, no matter what population size is used, extremely good results are obtained with errors less than 1%. For population sizes larger than 20, results of comparable accuracy are obtained after the second stage. The computation time required for the larger population sizes of 40, 60, 80 and 100 is of course much longer than that for a smaller population size. However for a population size of 5, at least four stages are needed to get results comparable to those obtained from two stages with population sizes larger than 20. In other words, a population size that is too small reduces the speed by requiring more optimization stages, while a population size that is too large increases the computation time in each optimization stage. A population size of 30 is considered appropriate to the parameter identification of a single moving vehicle.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
105
Table 5 Case 3: effects of measurement noise on accuracy (station at x=60 m, population size=30) Noise level (%)
5 10 15 20 25
Number of stages
2 3 3 3 3
Identified values
Errors (%)
M1 (kg)
M2 (kg)
c (Ns/m)
k (N/m)
M1
M2
c
k
4775 4812 4774 4783 4762
27 042 27 136 26 962 26 783 27 185
85 668 86 474 85 855 85 491 86 006
9 107 513 9 113 098 9 162 722 9 127 093 9 147 911
0.52 0.25 0.54 0.35 0.79
0.15 0.50 0.14 0.80 0.69
0.39 0.55 0.17 0.59 0.01
0.14 0.08 0.47 0.08 0.31
325
Contact force (kN)
320
315
310
305
300
0
0.1
0.2
0.3
0.4
0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 5. Case 3: contact forces calculated from the true and identified vehicle parameters using measurements of different noise levels: —— true; – – – 5% noise; — — 10% noise.
The histories of contact forces calculated from the identified parameters using different population sizes are plotted in Fig. 4 and compared with the true values. Very good agreement is observed. 4.4. Case study 3: effect of measurement noise The effect of measurement noise on the accuracy and efficiency of the proposed method is next investigated. This case is essentially the same as Case 1 except that five measurement noise levels of 5%, 10%, 15%, 20% and 25% are introduced. Following the conclusions of the studies for Case 2, a population size of 30 is adopted for this identification problem with a single moving vehicle. Table 5 shows the identification results, from which it can be seen that the identification values vary little with the increase of the measurement noise level. This demonstrates that the proposed identification procedure is not particularly sensitive to measurement noise levels. For the five noise levels, the errors are still well below 1%. However the more measurement noise there is, the more optimization stages are required. Fig. 5 shows the identified contact forces based on noise-free measurements, and measurements with noise levels of 5% and 10%. Basically, very
ARTICLE IN PRESS 106
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
Table 6 Case 4: arrangement of measurement stations Numbers of measuring stations
Locations of measuring stations, x (m)
1 2 3 4 6
60 18, 18, 18, 12,
60 60, 102 52, 68, 102 24, 52, 68, 96, 108
Table 7 Case 4: effects of number of measurement stations Number of measurement stations Number of stages Identified values
Errors (%)
M1 (kg) M2 (kg) c (Ns/m) k (N/m)
M1
M2
c
k
5% measurement noise, population size=30 1 2 2 2 3 2 4 2 6 2
4775 4813 4808 4814 4781
27 042 27 019 26 982 26 875 26 967
85 668 86 424 86 149 85 887 85 777
9 107 513 9 110 873 9 106 786 9 087 109 9 116 734
0.52 0.27 0.17 0.30 0.40
0.15 0.07 0.07 0.09 0.12
0.39 0.49 0.17 0.07 0.26
0.14 0.10 0.15 0.36 0.04
10% measurement noise, population size=30 1 3 2 3 3 3 4 2 6 2
4812 4793 4827 4773 4799
27 136 27 038 26 882 27 012 26 962
86 474 86 464 85 798 85 767 85 777
9 113 098 9 110 588 9 097 493 9 099 936 9 115 070
0.25 0.15 0.56 0.56 0.02
0.50 0.14 0.44 0.04 0.14
0.55 0.54 0.23 0.27 0.26
0.08 0.10 0.25 0.22 0.05
20% measurement noise, population size=30 1 3 2 3 3 3 4 2 6 2
4783 4844 4821 4819 4814
26 783 26 948 27 024 27 106 27 011
85 491 86 114 85 743 85 862 85 777
9 127 093 9 088 694 9 106 786 9 123 056 9 092 955
0.35 0.92 0.44 0.39 0.30
0.80 0.19 0.09 0.39 0.04
0.59 0.13 0.30 0.16 0.26
0.08 0.34 0.15 0.03 0.30
good agreement is observed although the presence of noise in the measurements does give rise to certain discrepancies. 4.5. Case study 4: effect of number of measurement stations Very often more than one accelerometer is used in field measurements. It not only improves the accuracy in most cases, but also caters for the possibility of occasional breakdown of certain instruments. The case of an increased number of measurement stations is studied. For comparison, the identification problem is repeated using five patterns of measurement stations, namely, 1, 2, 3, 4 and 6 measurement stations, respectively. The distribution and locations of the
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
107
325
Contact force (kN)
320
315
310
305
300
0
0.1
0.2
0.3
0.4 0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 6. Case 4: contact forces calculated from the true and identified vehicle parameters using different number of measurement stations: —— true; – – – One station; — — three stations. Table 8 Case 5: assume true vehicle parameters Vehicle
Unsprung mass (kg)
Sprung mass (kg)
Damping coefficient (Ns/m)
Spring stiffness (N/m)
1 2 3
4800 3900 5400
27 000 24 000 32 000
86 000 72 000 100 000
9.12 106 8.68 106 10.0 106
measurement stations for each arrangement are summarized in Table 6. A population size of 30 is used and three levels of measurement noise, namely 5%, 10% and 20%, are considered. The identification results as shown in Table 7 are very good. The errors are well below 1% even for 20% measurement noise. Similar to what is seen before, the higher the measurement noise, the more stages are needed. With the same level of measurement noise, the efficiency of optimization in general increases with the number of measurement stations, as seen in the number of stages required. However this is achieved at the expense of computation time. The time histories of the contact forces calculated from the true and identified values of the vehicle parameters using 1 and 3 measurement stations with 10% measurement noise are plotted in Fig. 6. It is obvious that the time histories of contact forces calculated from identified parameters both agree well with the true time histories, but that identified based on three measurement stations is closer to the true one. 4.6. Case study 5: effect of number of vehicles In real life, there are always many vehicles on a bridge at the same time. In order to test the accuracy and efficiency of the proposed method for more practical situations, a convoy comprising three vehicles at a spacing of 10 m is considered and their true parameters are shown in Table 8. Both the two and six-station measurement arrangements are used, and their locations are shown in Table 6. The measurement noise for acceleration simulation is taken to be 10%. The
ARTICLE IN PRESS 108
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
Table 9 Case 5: effect of number of vehicles (10% measurement noise, population size=60) Parameters
Identified values
Errors (%)
Number of stages
Two stations
six stations
Two stations
Six stations
Vehicle 1
M11 (kg) M12 (kg) c1 (Ns/m) k1 (N/m)
4757 27266 85790 9237977
4807 27065 86196 9131049
0.90 0.99 0.24 1.29
0.15 0.24 0.23 0.12
Vehicle 2
M21 (kg) M22 (kg) c2 (Ns/m) k2 (N/m)
3861 23988 73212 8711342
3881 23979 71482 8709654
1.00 0.00 1.68 0.36
0.49 0.09 0.72 0.34
M31 (kg) M32 (kg) c3 (Ns/m) k3 (N/m)
5465 32276 101615 10068262
5429 31901 100165 10062176
1.20 0.86 1.62 0.68
0.54 0.31 0.17 0.62
Vehicles 3
4
340
Contact force (kN)
330
320
310
300
290
0
0.1
0.2
0.3
0.4 0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 7. Case 5: contact forces of the first vehicle calculated from the true and identified vehicle parameters using two measuring stations: —— true; – – – identified.
same initial search domains as in Case 1 are used. In view of the larger number of parameters to be identified in this case, the population size is set as 60. The identification results as shown in Table 9 are very good. For example, in the two-station measurement arrangement, identification errors below 2% can be obtained after only four optimization stages even though there are 12 parameters to be identified. In the six-station measurement arrangement, identification error below 1% can be obtained after only four optimization stages. It demonstrates that the proposed method can also cope with a convoy of vehicles with reasonable accuracy although longer computation time is required for more vehicles. It is also seen that the accuracy improves when more measurement stations are used. For the present three-span continuous bridge, six measurement stations are sufficient for parameter identification of the vehicle convoy.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
109
300
Contact force (kN)
290
280
270
260
250
0
0.1
0.2
0.3
0.4 0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 8. Case 5: contact forces of the second vehicle calculated from the true and identified vehicle parameters using two measuring stations: —— true; – – identified.
410
Contact force (kN)
395 380 365 350 335 320
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time ratio t / T
Fig. 9. Case 5: contact forces of the third vehicle calculated from the true and identified vehicle parameters using two measuring stations: —— true; – – – identified.
The time histories of contact forces of the three vehicles calculated from the true and identified parameters using two measurement stations are plotted in Figs. 7–9, respectively. Note that the total traversing time T is counted from the instant when the first vehicle arrives at the starting end of the bridge to that when the last vehicle leaves the bridge. Therefore the time history for each vehicle only occupies part of the total traversing time. Good agreement is again observed.
5. Discussions In Case 1 when noise-free simulated accelerations are used, effectively the exact results are obtained. It shows that the proposed procedure is theoretically correct and accurate. Then in the other cases, the procedure is tested on various aspects including the population size, noise level of
ARTICLE IN PRESS 110
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
acceleration measurements, number of measurement stations, and number of vehicles. Examination of the results confirms that the identification method is very robust. Even for a larger number of variables and a reasonably high noise level, very acceptable results can be obtained after only a few optimization stages with an appropriate population size. In all the five cases studied above, even only a few measurement stations are chosen along the whole bridge, very acceptable results can be obtained no matter how many vehicles are on the bridge at the same time.
6. Conclusions A procedure for the identification of vehicles moving on multi-span continuous bridges is described. The moving vehicles are modelled as 2-d.o.f. systems each comprising an unsprung mass and a sprung mass interconnected by a damper and a spring. The aim of the study is to identify the associated vehicle parameters based on the acceleration measurements at selected stations. Here the acceleration measurements are simulated by the numerical solution to the associated forward problem with random noise introduced to account for the likely errors in practice. The search for the vehicle parameters is formulated as an optimization problem in which the error between the measured accelerations and the accelerations reconstructed from trial vehicle parameters is minimized. A robust multi-stage optimization scheme based on genetic algorithms has been proposed in which the search domains are systematically reduced to converge on the optimal solution. A comprehensive study has been carried out to test the proposed procedure. It shows that in spite of the presence of measurement noise, accurate results are still obtained. The method can also cope with the presence of multiple vehicles and the use of instruments at multiple sections.
Acknowledgements The work described in this paper has been supported by the Committee on Research and Conference Grants, The University of Hong Kong, China. The present study has made use of the code of a genetic algorithm written by Dr. D.L. Carroll at the University of Illinois, and this is gratefully acknowledged.
References [1] C.A. Julien, H.M. Hung, Technique in system identification for dynamic mechanical systems, SAE-Paper 690497, 1969, p. 9. [2] R. Schlegemilch, M. Railey, System identification comparison study, Proceedings of the 21st Annual Pittsburgh Conference Part 5 (of 5), Pittsburgh, PA, May 3–4, 1990, pp. 2269–2273. [3] Y. Kyongsu, K. Hedrick, Observer-based identification of nonlinear system parameters, Journal of Dynamic Systems, Measurement and Control, Transactions of the American Society of Mechanical Engineers 117 (2) (1995) 175–182.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
111
[4] K. Wise, R.E. Reid, Modeling and identification of a light truck engine mounting system for ride quality optimization, International Congress & Exposition—Society of Automotive Engineers, Warrendale, PA, 1984, 841142P. [5] D.A. Derradji, N. Mort, Multivariable adaptive control using artificial neural networks, Proceedings of the 1996 UKACC International Conference on Control. Part 2 (of 2), Exeter, UK, September 2–5 1996, 427/2, pp. 889–893. [6] P. Michelberger, J. Bokor, A. Keresztes, P. Varlaki, Identification of a multivariable linear model for road vehicle dynamics from test data, International Journal of Vehicle Design 8 (1) (1987) 96–114. [7] K.R. Goheen, E.R. Jefferys, The application of alternative modelling techniques to ROV dynamics. Proceedings of the 1990 IEEE International Conference on Robotics and Automation, Cincinnati, OH, May 13–18 1990, pp. 1302–1309. [8] D.H. Lee, J.G. Lee, Comparison of parameter identification algorithms for flight vehicles, Transactions of the Japan Society for Aeronautical and Space Sciences 36 (114) (1994) 249–256. [9] E. Bisimis, Use of system identification methods for the determination of vehicle parameters from measured time histories, Vehicle System Dynamics 6 (1977) 182–184. [10] O. Kropac, J. Sprinc, Identification of the system vehicle-road parameters, Vehicle System Dynamics 11 (4) (1982) 241–249. [11] L.R. Ray, Nonlinear tire force estimation and road friction identification: simulation and experiments, Automatica 33 (10) (1997) 1819–1833. [12] D.Y. Zheng, Y.K. Cheung, F.T.K. Au, Y.S. Cheng, Vibration of multi-span non-uniform beams under moving loads by using modified beam vibration functions, Journal of Sound and Vibration 212 (3) (1998) 455–467. [13] Y.K. Cheung, F.T.K. Au, D.Y. Zheng, Y.S. Cheng, Vibration of multi-span non-uniform bridges under moving vehicles and trains by using modified beam vibration functions, Journal of Sound and Vibration 228 (3) (1999) 611–628. [14] R.J. Jiang, F.T.K. Au, Y.K. Cheung, Identification of masses moving on multi-span beams based on genetic algorithms, Computers and Structures 81 (22–23) (2003) 2137–2148. [15] D.E. Goldberg, Genetic Algorithms in Search Optimization, and Machine Learning, Addison-Wesley Reading, MA, 1989.
Journal of Sound and Vibration 269 (2004) 91–111
JOURNAL OF SOUND AND VIBRATION www.elsevier.com/locate/jsvi
Parameter identification of vehicles moving on continuous bridges F.T.K. Au*, R.J. Jiang, Y.K. Cheung Department of Civil Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong, People’s Republic of China Received 5 June 2002; accepted 19 December 2002
Abstract This paper describes a method for the identification of parameters of vehicles moving on multi-span continuous bridges. Each moving vehicle is modelled as a 2-degree-of-freedom system that comprises four components: an unsprung mass and a sprung mass, which are connected together by a damper and a spring. The corresponding parameters of these four components, namely, the equivalent unsprung mass and sprung mass, the damping coefficient and the spring stiffness are identified based on dynamic simulation of the vehicle–bridge system. The identification process makes use of acceleration measurements at selected stations on the bridge. In the study, the acceleration measurements are simulated from the solution to the forward problem of a continuous beam under moving vehicles, together with the addition of artificially generated measurement noise. The identification is carried out through a robust multi-stage optimization scheme based on genetic algorithms, which searches for the best estimates of parameters by minimizing the errors between the measured accelerations and the reconstructed accelerations from the moving vehicles. This multi-stage optimization scheme reduces the variable search domains stage by stage using the identified results of the previous stage. Besides the basic operators in simple genetic algorithms, some advanced genetic operators and techniques are adopted here. Therefore, it makes the proposed identification procedure much more efficient than other traditional optimization methods. The identification procedure is then verified with a few test cases. The estimated vehicle parameters can also be used to get the time varying contact forces between the vehicles and bridge surface. r 2003 Elsevier Ltd. All rights reserved.
1. Introduction Parameter identification is an important kind of problem relevant to many fields, and various researchers have experimented with different optimization methods and techniques. For the *Corresponding author. Tel.: +1-2859-2650; fax: +1-2559-5337. E-mail address: [email protected] (F.T.K. Au). 0022-460X/03/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-460X(03)00005-1
ARTICLE IN PRESS 92
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
identification of vehicle parameters, one of the traditional tools is the least squares method [1–3] and its amended versions [2]. The other methods for identification of vehicle parameters include the extended Kalman filter method [4,5], the maximum likelihood method [6], and the recursive prediction error method [7]. A comparative study [8] has been carried out on the last three methods of identification of vehicle parameters, and the study has discovered that the most common problem with these optimization methods is convergence. Different models have been created to represent the actual behaviour of vehicles with different emphasis. For models of vehicle–tyre system with different numbers of degrees-of-freedom (d.o.f.s) [1,3,6,8], different numbers of system parameters are identified. In some cases, attention is paid to the driver–vehicle interaction analysis [9] and the vehicle parameters for riding quality optimization [4,10]. At other times, the aim of parameter identification is to get the tyre forces and the road friction coefficients for vehicles on different road surfaces [11]. Most of the identifications mentioned above are realized based on the simulated or experimental test data of the vehicles. From the point of view of dynamic coupling between the vehicle and the road, the motion of a bridge is largely related to the motion of a vehicle on it. This is well verified by the results of forward analysis of bridge dynamics under moving loads [12] and moving vehicles [13]. All the forward dynamic analysis of bridges under moving vehicles is based on the assumption that the vehicle parameter values are known correctly. However, for the normal operation and structural health monitoring of existing bridges open to public, it is unrealistic to expect that all the correct vehicle parameters are known beforehand. Therefore the identification of vehicle parameters is an important step in the study of dynamic behaviour of existing bridges. This paper describes a method for the identification of parameters of vehicles moving on multispan continuous bridges. Instead of just modelling the vehicles as moving masses [14], each moving vehicle is modelled as a 2-d.o.f. system [13] that comprises four components: an unsprung mass and a sprung mass, which are connected together by a damper and a spring. The corresponding parameters of these four components, namely, the equivalent unsprung mass and sprung mass, the damping coefficient and the spring stiffness are identified based on the acceleration measurements at various stations along the continuous beam. In the study, the acceleration measurements are simulated from the solution to the forward problem of a continuous beam under moving vehicles based on the modified vibration functions [12] and modal superposition, together with the addition of artificially generated measurement noise. The identification is carried out through a robust multi-stage optimization scheme based on genetic algorithms [15], which searches for the best estimates of parameters by minimizing the errors between the measured accelerations and the reconstructed accelerations from the moving vehicles. This multi-stage optimization scheme reduces the variable search domains stage by stage using the identified results of the previous stage. Based on the mechanics of natural selection and natural genetics, genetic algorithms combine survival of the fittest among string structures with a structured yet randomized information exchange to form a search algorithm with some of the innovative flair of human search. Their special search mechanics leads to their special robust function of finding the global optimal values of the variables definitely without worrying about the convergence problem nor the restrictive requirements of continuity and derivative existence, as compared with the traditional optimization techniques. Apart from the basic strategies of simple genetic algorithms, other advanced genetic techniques, such as niching, elitism, fitness scaling, are also adopted to improve the efficiency. In the present parameter identification problem, a natural
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
93
search domain that contains all the possible options for each variable does not exist, at least at the initial stage of solution. In order to overcome this problem and to improve the efficiency, a multistage optimization procedure based on genetic searching is put forward by reducing the variable search domains stage by stage using the identified results of the previous stage. A comprehensive case study is then carried out to verify the robustness of the proposed identification method. The effects of different noise levels, various parameters used in the genetic algorithms, the number of measurement stations, and the number of vehicles on the results are studied. The time-dependent contact forces between the moving vehicles and the bridge surface can also be calculated from the identified vehicle parameters. The work described in this paper will be useful to the operation and structural health monitoring of existing bridges.
2. Simulation of acceleration measurements 2.1. Modelling of vehicles and bridge In this paper, a multi-span continuous bridge is modelled as a linear elastic Bernoulli–Euler beam with totally (P+1) point supports as shown in Fig. 1. The bridge may have varying sections characterized by the density rðxÞ; cross-sectional area AðxÞ; Young’s modulus EðxÞ and moment of inertia IðxÞ at location x: The N moving vehicles on the bridge are modelled as 2-d.o.f. systems each with four parameters {Ms1, Ms2, cs, ks, s=1,2, y, N}. The typical sth moving vehicle is modelled by the unsprung mass Ms1 and the sprung mass Ms2 ; which are interconnected by a spring of stiffness ks and a dashpot of damping coefficient cs : The convoy of vehicles moves from left to right along the bridge at a speed vðtÞ: Setting the position of the leftmost point support of the bridge as the origin O; the horizontal position of the sth vehicle is defined as xs ðtÞ; a function of t. The deflection of the bridge at position x and at time t is denoted by yðx; tÞ where upward deflection is taken as positive, as shown in Fig. 1. The vertical displacements of the unsprung mass and sprung mass of sth vehicle are, respectively, represented by fys1 ðtÞ; ys2 ðtÞ; s ¼ 1; 2y; Ng; which are also measured vertically upwards with reference to their respective vertical positions before coming onto the bridge.
v(t) ys2(t)
Ms2 y
O
cs
ks Ms1
0
1
2
ys1(t) x P
xs(t) L
Fig. 1. A multi-span continuous bridge with N moving vehicles.
ARTICLE IN PRESS 94
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
If the surface roughness of the bridge and the separation of the vehicle from the bridge are neglected, the vertical displacement ys1 ðtÞ of the unsprung mass of the sth vehicle model is the same as the deflection of the beam yðxs ðtÞ; tÞ at the corresponding location xs ðtÞ: Then the vertical displacement ys1 ðtÞ and its first two derivatives are as follows: ys1 ðtÞ ¼ yðx; tÞjx¼xs ðtÞ ; dys1 ðtÞ @yðx; tÞ @yðx; tÞ ¼ þ v ; dt @t @x x¼xs ðtÞ 2 2 d2 ys1 ðtÞ @ yðx; tÞ @2 yðx; tÞ @yðx; tÞ 2 @ yðx; tÞ ¼ þ 2v þ a ; þ v dt2 @t2 @x @t @x2 @x x¼xs ðtÞ
ð1Þ ð2Þ
ð3Þ
where a is the horizontal acceleration of the group of vehicles. The contact force fcs ðtÞ between the sth unsprung mass and the beam can be expressed as fcs ðtÞ ¼ ðMs1 þ Ms2 Þg þ Dfcs ðtÞ;
ð4Þ
in which g is the acceleration due to gravity and Dfcs ðtÞ is the time-dependent variation of contact force. The equation of vertical motion for the sth unsprung mass Ms1 and sprung mass Ms2 can be written as d2 ys1 ðtÞ dys2 ðtÞ dys1 ðtÞ ½ Ms1 ¼ Df ðtÞ þ k y ðtÞ y ðtÞ þ c ; ð5Þ cs s s2 s1 s dt2 dt dt d2 ys2 ðtÞ dys2 ðtÞ dys1 ðtÞ Ms2 ¼ ks ½ys2 ðtÞ ys1 ðtÞ cs : ð6Þ dt2 dt dt From Eqs. (4)–(6), the contact force between the sth vehicle and the bridge can be written as d2 ys1 ðtÞ d2 ys2 ðtÞ þ M : ð7Þ fcs ðtÞ ¼ ðMs1 þ Ms2 Þg þ Ms1 s2 dt2 dt2 For the modelling of the continuous bridge, the modified vibration functions [12] are chosen, which are built up from the vibration modes of a hypothetical single-span simply supported beam having the same end supports and total length as the real beam, and cubic spline expressions that ensure the boundary conditions at the supports. Therefore the assumed vibration modes can be written as * i ðxÞ; % i ðxÞ þ F ð8Þ Fi ðxÞ ¼ F % i ðxÞ; i ¼ 1; 2; y; ng are the first n vibration modes of a hypothetical beam of total length where fF * i ðxÞ; i ¼ 1; 2; y; ng L with the same end supports but without the intermediate supports, and fF are the augmenting cubic spline expressions which are so chosen that each Fi ðxÞ satisfies the boundary conditions at the two ends and the zero deflection conditions at the intermediate point % i ðxÞ are the supports. For a hypothetical beam that is simply supported, the vibration modes F * Fourier sine series. The augmenting cubic spline expressions Fi ðxÞ can be determined by using the zero deflection conditions at the intermediate supports and the boundary conditions at end supports.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
95
By modal superposition and separation of variables, the vibration of the beam can be expressed as yðx; tÞ ¼
n X
qi ðtÞFi ðxÞ
ði ¼ 1; 2; y; nÞ;
ð9Þ
i¼1
where qi ðtÞ are the generalized coordinates. 2.2. Simulation of noise-free acceleration measurements Based on the vehicle and bridge models introduced above, dynamic analysis of the system can be carried out and the simulated noise-free acceleration measurements at specified stations of the bridge can be obtained for the subsequent parameter identification of the vehicles. This part is in fact a conventional forward problem. The velocity of vibration and the curvature of the beam at position x are, respectively, n @yðx; tÞ X ¼ q’ i ðtÞ Fi ðxÞ; @t i¼1
ð10Þ
n @2 yðx; tÞ X ¼ qi ðtÞF00i ðxÞ; @x2 i¼1
ð11Þ
in which the dot represents differentiation with respect to time t, and the prime represents differentiation with respect to the abscissa x. Therefore the kinetic energy T, the bending energy V, and the work done by external forces W are, respectively, Z n X n 1 L @yðx; tÞ 2 1X rðxÞAðxÞ dx ¼ ð12Þ T¼ q’ i ðtÞmij q’ j ðtÞ; 2 0 @t 2 i¼1 j¼1 1 V¼ 2
Z
L 0
W ¼
2 2 n X n @ yðx; tÞ 1X EðxÞIðxÞ dx ¼ qi ðtÞkij qj ðtÞ; @x2 2 i¼1 j¼1
Z 0
N LX
fcs yðx; tÞd½x xs ðtÞ dx ¼
s¼1
n X
qi ðtÞfi ðtÞ;
ð13Þ
ð14Þ
i¼1
where the work done by external forces W is defined in terms of the Dirac delta function ( N ðx ¼ x0 Þ ð15aÞ dðx x0 Þ ¼ 0 ðxax0 Þ Z
N
dðx x0 Þ dx ¼ 1; N
ð15bÞ
ARTICLE IN PRESS 96
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
while mij ; kij and fi ðtÞ are, respectively, the generalized mass, stiffness and force, which can be written as Z L mij ¼ rðxÞAðxÞFi ðxÞFj ðxÞ dx; ð16Þ 0
kij ¼
Z
L
EðxÞIðxÞF00i ðxÞF00j ðxÞ dx;
0
fi ðtÞ ¼
N X
ð fcs Fi ½xs ðtÞ Þ ¼
s¼1
N X
fis ðtÞ:
ð17Þ
ð18Þ
s¼1
The component fis ðtÞ in the generalized force fi ðtÞ can be further written as fis ðtÞ ¼ ðMs1 þ Ms2 ÞgFi ðxs ðtÞÞ Ms2 y. s2 ðtÞFi ðxs ðtÞÞ n X Ms1 Fi ðxs ðtÞÞ fq. j ðtÞFj ðxs ðtÞÞ þ 2vq’ j ðtÞF0j ðxs ðtÞÞ þ qj ðtÞ½v2 F00j ðxs ðtÞÞ þ aF0j ðxs ðtÞÞ g:
ð19Þ
j¼1
The Lagrangian function for the vehicle-bridge system is G ¼ T ðV W Þ; where ðV W Þ is the total potential energy, and the Euler-Lagrange equation appears as d @G @G ¼ 0: ð20Þ dt @q’ i @qi Substituting Eqs. (16–17) and (18) into Eq. (20) gives n n n X X X # ij q. j ðtÞ þ m c#ij q’ j ðtÞ þ k#ij qj ðtÞ ¼ p# i ðtÞ; j¼1
j¼1
i ¼ 1; 2; y; n;
ð21Þ
j¼1
where # ij ðtÞ ¼ mij þ m
N X
Ms1 Fi ðxs ðtÞÞFj ðxs ðtÞÞ;
ð22Þ
s¼1
c#ij ðtÞ ¼
N X
2vMs1 Fi ðxs ðtÞÞF0j ðxs ðtÞÞ;
ð23Þ
s¼1
k#ij ðtÞ ¼ kij þ
N X
h i Ms1 Fi ðxs ðtÞÞ v2 F00j ðxs ðtÞÞ þ aF0j ðxs ðtÞÞ ;
ð24Þ
s¼1 N X ðMs1 þ Ms2 ÞgFi ðxs ðtÞÞ: p# i ðtÞ ¼
ð25Þ
s¼1
Should a particular vehicle be outside the bridge, the corresponding terms under the summation signs should be omitted. Eq. (21) can be solved by the Newmark-b method. Therefore fqi ðtÞg and fq. i ðtÞg can be obtained accordingly. The simulated ‘‘noise-free’’ acceleration anf ðx; tÞ of the beam at measurement
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
97
location x can then be expressed as anf ðx; tÞ ¼ yðx; . tÞ ¼
n X
q. i ðtÞFi ðxÞ:
ð26Þ
i¼1
The vector containing the simulated ‘‘noise-free’’ acceleration time histories at the measurement stations then forms the noise-free acceleration vector Anf : 2.3. Simulation of measurement noise In practical experiments or fieldwork, the measurements are usually polluted by measurement noise. To evaluate the sensitivity of the results to such measurement noise, noise-polluted measurements are simulated by adding to the noise-free acceleration vector a corresponding noise vector whose root-mean-square (RMS) value is equal to a certain percentage of the r.m.s. value of the noise-free data vector. The components of all the noise vectors are of Gaussian distribution, uncorrelated and with a zero mean and unit standard deviation. Then the noise-polluted acceleration Anp of the bridge at location x can be simulated by Anp ¼ Anf þ RMS ðAnf Þ Nlevel Nunit ;
ð27Þ
where RMSðAnf Þ is the root-mean-square value of the noise-free acceleration vector Anf ; Nlevel is the noise level, and Nunit is the randomly generated noise vector with zero mean and unit standard deviation. In the present study, besides the case of zero noise, five levels of noise pollution will be investigated, namely 5, 10, 15, 20% and 25%.
3. Multi-stage parameter identification based on genetic algorithms 3.1. Genetic algorithms Genetic algorithms [15] are stochastic global search techniques based on the mechanics of natural selection and natural genetics. They combine survival of the fittest among string structures with a structured yet randomized information exchange to form a search algorithm with some of the innovative flair of human search. The central theme of genetic algorithms is robustness, which results from their special search ways, namely, working with a coding of the parameter set, not the parameters themselves; searching from a population of points, not a single point; using objective function information, not derivatives or other auxiliary knowledge; and using probabilistic transition rules, not deterministic rules. Because of these intrinsic characteristics, the search and optimization direction of genetic algorithms is definitely guided to global optimal values of the variables instead of the local optimal solution in traditional search and optimization methods. Convergence is not a problem for genetic algorithms at all. Genetic algorithms have found applications in a wide spectrum of areas, covering problems with discrete variables and continuous variables. In the present study, the variables to be identified are the equivalent unsprung mass and sprung mass, the damping coefficient and the spring stiffness for each vehicle. Assuming that there are all together N vehicles on a bridge, the number of parameters to be identified is therefore 4N. All the
ARTICLE IN PRESS 98
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
{M11, M12, …, M1N} {M21, M22, …, M2N} {c1, c2, …, cN} {k1, k2, …, kN} N
N
N
N
4N
Fig. 2. The structure of an individual with 4N variables in the present study.
parameters are coded as binary strings with different string lengths, where the string length is defined as the number of the binary bits for one variable. For the purpose of optimization using genetic algorithms, these strings are assembled together in prescribed order to form individuals, and the structure of an individual adopted in the present study is shown in Fig. 2. To kick off the search and optimization process, a specified number of individuals must be randomly created to form an initial population. 3.2. Genetic techniques adopted The three basic genetic operations in simple genetic algorithms are adopted here, namely reproduction of individuals, crossover between individuals in a population, and the mutation of individuals. Besides these, a few advanced genetic techniques are used, which improves the identification efficiency in each stage. After an initial population is randomly created with the specified population size, the tournament selection scheme is used with a shuffling technique for choosing random pairs for mating. In order to keep the fittest individuals in the later populations, elitism is also included, which means that the best individuals in a generation will be reproduced to the next generation. Uniform crossover is adopted for the proposed identification process. Jump mutation or creep mutation scheme is used to realize mutation. Niche theory is also introduced. A sharing function used to induce niche in this genetic algorithm limits the uncontrolled growth of particular species within a population. Another important improvement scheme—linear fitness scaling [15] is implemented to keep appropriate levels of competition throughout a simulation. Without scaling of the objective function values, there is a tendency for a few super-individuals to dominate the selection process early on during a run, which is a leading cause of premature convergence. In this case, the objective function values must be scaled down to prevent a takeover of the population by these super-strings. Later on during the run, when the population is largely converged, competition among population members is less strong and the simulation tends to wander. In this case, the objective function values must be scaled up to accentuate differences between population members to continue to reward the best performers. In both cases, at the beginning of the run and as the run matures, the linear fitness scaling can help. The fitness value is related to the objective function to reflect the desirability of the individual. Linear scaling requires a linear relationship between the scaled fitness f 0 and the raw fitness f as follows: f 0 ¼ af þ b:
ð28Þ
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
99
The coefficients a and b are chosen based on two criteria. The first is to enforce equality of the 0 of the scaled fitness values. average value favg of the raw fitness values and the average value favg The second criterion is to control the number of offspring given to the population member with maximum raw fitness. In particular, the best population member is expected to give Cmult copies in the next generation, i.e., 0 ¼ Cmult favg ; fmax
ð29Þ
0 where fmax is the maximum scaled fitness. Here Cmult is taken as 2.0, which means that the number of expected copies desired for the best population member is 2. In the present identification of vehicle parameters, the errors at the measurement stations between the measured or simulated accelerations and the reconstructed accelerations from the identified vehicle parameters should be minimized. For a bridge with Q measurement stations and a total of R time points in each measured time history, the objective function F to be maximized in the present study can be written in terms of the error vector ei as !1=2 Q Q R X X X 2 F ¼ jjei jj2 ¼ ðeij Þ ; ð30Þ i¼1
i¼1
j¼1
where jjei jj2 is the Euclidean vector norm of the error vector ei defined as ei ¼ faij;simu g faij;iden g;
ð31Þ
where aij;simu is the element of the simulated acceleration vector and aij;iden is the element of the acceleration vector calculated from the trial identified vehicle parameters in each generation, both at the ith measurement station and the jth time point. Linear scaling works well except when negative fitness calculation prevents its use [15]. To overcome this problem associated with Eq. (30), a mapping from the objective function F to the raw fitness f is carried out as f ¼ F þ Cmin ;
ð32Þ
where Cmin is chosen to enforce that the raw fitness f is always positive. In the present study, the value of Cmin is taken to be the absolute value of the worst objective function value in each current generation, which prevents the occurrence of the negative fitness in every generation. 3.3. Multi-stage optimization In the present parameter identification problem, a natural search domain that contains all the possible options for each variable does not exist obviously, at least at the initial stage of solution. Although the identification of the contact forces provides useful clues to the identification of moving masses [14], this may not be as useful here. In order not to miss the global optimal set of parameters and converge on a wrong solution, the initial search domain for each variable must be set sufficiently large. However within these very large search domains, it is often difficult to find the ideal global optimal values by a single optimization run, and it also takes an extremely long running time. Without any idea of the likely magnitudes of the optimal parameters, it is difficult to determine the appropriate string lengths, which are essential for accurate results. In order to overcome this problem, the parameter identification is carried out through a multi-stage
ARTICLE IN PRESS 100
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
optimization scheme involving the repeated use of genetic algorithms and gradual refinement of search domains as well as accuracy. Using reasonable string lengths and a relatively small number of generations, it is possible for the initial stage of optimization to produce rough estimates of the parameters. Such estimates are then used to map out the refined search domains for the second stage. It is known that the genetic optimization gradually guides the trial solutions to the optimal one. Therefore through this staged optimization process, the search domains will be reducing in each successive stage in the light of the identification results of the previous stage until the desired identification precision is reached. In the present study, all parameters to be identified are non-negative values, and so the lower limits of their initial search domains are all set as zero. The upper limits of the initial search domains for the spring stiffness and the damping coefficient are set to 1020 or even larger values in standard units. For both the unsprung and sprung masses for normal bridges, a convenient upper limit of the initial search domain can be taken as 106 kg. These search domains can then be narrowed down in subsequent stages of optimization. For example, assuming the identified unsprung mass of the sth vehicle in a stage to be Ms1;iden ; the search domain D of this unsprung mass for the next stage can be set as D ¼ ðDl ; Du Þ ¼ Ms1;iden ðCl ; Cu Þ;
0oCl o1
1oCu ;
ð33Þ
where Cl and Cu are, respectively, coefficients to define the lower limit Dl and the upper limit Du of the search domain D. The same strategy can be applied to the search for parameters for the sprung mass Ms2 ; the spring stiffness ks and the damping coefficient cs : The values of Cl and Cu may be varied at different stages for improved efficiency. The termination of the search at each stage is controlled both by the changes of parameters between successive generations as well as the maximum number of generations allowed. However the termination of the whole optimization procedure is controlled by the differences between the identified parameters of consecutive stages. The efficiency of this multi-stage optimization process based on genetic algorithms at each stage will be verified by a comprehensive case study.
4. Case study 4.1. General description To evaluate the efficiency and effectiveness of the proposed method, a comprehensive case study covering various aspects have been carried out. The effects of population sizes, number of measurement stations, levels of measurement noise, and number of vehicles on the accuracy are examined. Table 1 summarizes the different aspects examined in various cases. The contact forces between the moving vehicles and the bridge surface are also estimated using the identified vehicle parameters. The identification error Erroriden of, for example, the unsprung mass of sth vehicle is defined as Erroriden ¼
jMs1;iden Ms1;true j
100%; Ms1;true
ð34Þ
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
101
Table 1 Cases for investigation Case
No. of vehicles
No. of stations
Noise level (%)
Population size
Effects to investigate
1 2 3 4 5
1 1 1 1 3
1 1 1 1, 2, 3, 4, 6 2, 6
0 5 5, 10, 15, 20, 25 10 10
30 5, 10, 20, 30, 40, 60, 80, 100 30 30 60
Basic strategy Population size Measurement noise Number of stations Number of vehicles
36 m
48 m
36 m
L = 120 m (a) 4m 0.2 m
0.4 m
0.4 m
3.0 m 0.2 m
0.75 m
2.5 m
0.75 m
(b)
Fig. 3. A three-span continuous bridge. (a) Idealized beam model; (b) idealized cross-section.
where Ms1;iden and Ms1;true are the identified and true values of unsprung mass Ms1 ; respectively. A similar definition is also applied to other parameters such as the sprung mass, spring stiffness and damping coefficient. A three-span continuous box girder bridge is used in the case study, and the idealized beam model and cross-section are shown in Fig. 3. The density r and Young’s modulus E of concrete are 2400 kg/m3 and 30 000 Mpa, respectively. The span lengths from left to right are 36, 48, 36 m, respectively. The cross-sectional area A and the second moment of inertia I are 3.46 and 3.854 m4, respectively. For both the acceleration simulation and parameter identification, the first three
ARTICLE IN PRESS 102
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
vibration modes of the bridge are adopted. All the vehicles are assumed to travel at a horizontal speed v of 17 m/s across the bridge. Both for the simulation of acceleration measurements and the subsequent analysis, the total traversing time T, which is counted from the instant when the first vehicle arrives at the starting end of the bridge to that when the last vehicle leaves the bridge, is divided into 200 time intervals. The histories of acceleration and contact forces are therefore presented accordingly with respect to time ratio t/T. The contact forces between the vehicles and the bridge surface are calculated from fcs ðtÞ ¼ ðMs1 þ Ms2 ÞgFi ðxs ðtÞÞ Ms2 y. s2 ðtÞFi ðxs ðtÞÞ n X Ms1 Fi ðxs ðtÞÞ fq. j ðtÞFj ðxs ðtÞÞ þ 2vq’ j ðtÞF0j ðxs ðtÞÞ þ qj ðtÞ½v2 F00j ðxs ðtÞÞ þ aF0j ðxs ðtÞÞ g:
ð35Þ
j¼1
The parameters for the present genetic algorithm with uniform crossover are as follows: crossover probability Pc ¼ 0:5; creep mutation probability Pcm ¼ 0:04; and jump mutation probability Pjm ¼ 1=Npop ; where Npop denotes the population size. As seen in the following cases, where the signals are not contaminated with measurement noise, essentially the true parameters are obtained. 4.2. Case study 1: basic strategy In Case 1, there is only one moving vehicle with simulated noise-free accelerations at the measurement station located at the middle of the central span, i.e., x ¼ 60 m: This is mainly to establish the validity and effectiveness of this multi-stage optimization procedure. The true parameters for the vehicle used are the unsprung mass M1=4800 kg, the sprung mass M2=27000 kg, the damping coefficient c=86000 Ns/m, and the spring stiffness k=9.12 106 N/m. These parameters are applicable to the other cases unless otherwise stated. In the first optimization stage, the initial search domains of the four vehicle parameters are set as M1 Að0; 10 000Þ kg; M2 Að0; 60 000Þ kg; cAð0; 200 000Þ Ns/m and kAð0; 2 107 Þ N=m; and their string lengths are 7, 8, 8 and 8, respectively. The population size is set as 30. By generation 300 of the first stage, the identified parameters M1, M2, c and k are 4961, 26 824 kg, 89 412 Ns/m and 9.098 106 N/m, respectively. Based on these identified values, and taking the coefficients Cl and Cu for Eq. (33) as 0.9 and 1.1, respectively, the search domains are narrowed down to M1A(4465, 5457) kg, M2A(24 142, 29 506) kg, cA(80 471, 98 353) Ns/m and kA(8.188 106, 10.0078 106) N/ m. Similarly, the genetic search is carried out again within these refined search domains. The identified parameters in this stage are then used to map out the search domains for the next stage. In this manner, the correct results are gradually identified as shown in Table 2. Only three stages are necessary to get the essentially true values of the parameters. The redefinition of search domains for successive stages enables the identification process to be both accurate and efficient. Table 3 shows the evolution of results for the first stage. The maximum number of generations for this stage is set as 800. It can be seen that among the first 90 generations, the results change a lot, while the results in generations 244–800 are essentially the same. Therefore the chosen number of generations is enough for the first stage, and it can be safely terminated.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
103
Table 2 Case 1: summary of results at different stages noise-free measurements, station at x=60 m, population size=30 Stage
1 2 3
Limit coefficients for search domain
Identified parameters
Cl
Cu
M1 (kg)
M2 (kg)
c (Ns/m)
k(N/m)
— 0.9 0.95
— 1.1 1.05
4961 4729 4802
26 824 27 292 27 002
89 412 88 167 86 038
9 098 039 9 224 698 9 120 620
4800 0.04
27 000 0.01
86 000 0.04
9 120 000 0.01
True value Error (%)
Table 3 Case 1: evolution of results in the first stage (noise-free measurements, station at x=60 m, population size=30) Generation number
Identified parameters M1 (kg)
M2 (kg)
c (Ns/m)
k (N/m)
1–5 6–7 8 9–11 12 13–17 18 19–24 25–26 27–28 29 30–31 32–33 34–39 40 41–54 55 56–69 70–73 74–88 89–243 244–800
7874 2992 6457 2835 4094 4094 3780 3780 4409 3937 4724 3937 3780 3937 4252 4094 4094 4330 4252 4252 4252 4961
25647 29412 25647 29882 29412 29412 29412 29412 29176 28941 28470 28470 28470 28470 28470 28471 28471 28235 28235 28235 28235 26824
71373 144314 69804 96470 94902 96471 97255 99608 99608 99608 99608 99608 99608 99608 99608 99608 98039 98039 98039 94902 96471 89412
9 254 902 9 725 490 9 333 333 9 725 490 9 882 353 9 725 490 9 725 490 9 725 490 9 725 490 9 647 059 9 568 627 9 568 627 9 411 765 9 490 196 9 490 196 9 490 196 9 490 196 9 490 196 9 411 765 9 411 765 9 411 765 9 098 039
Error (%)
3.35
0.65
3.97
0.24
4.3. Case study 2: effect of population size This is essentially the same as Case 1 except that 5% measurement noise is introduced and that the population size is varied to study its effect on the accuracy and efficiency of the method.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
104
Table 4 Case 2: effects of population size on accuracy. (5% measurement noise, station at x=60 m) Population size
5 10 20 30 40 60 80 100
Number of stages
4 3 3 2 2 2 2 2
Identified values
Errors (%)
M1 (kg)
M2 (kg)
c (Ns/m)
k (N/m)
M1
M2
c
k
4778 4758 4790 4775 4774 4775 4792 4780
27052 27035 26942 27042 27035 26957 26963 26963
86 728 86 147 85 267 85 668 85 806 86 229 85 845 86 192
9 118 394 9 095 260 9 106 979 9 107 513 9 118 959 9 109 327 9 141 072 9 113 324
0.46 0.88 0.21 0.52 0.55 0.52 0.17 0.42
0.19 0.13 0.22 0.15 0.13 0.15 0.14 0.14
0.85 0.17 0.85 0.39 0.23 0.27 0.18 0.22
0.02 0.27 0.14 0.14 0.01 0.12 0.23 0.07
325
Contact force (kN)
320
315
310
305
300
0
0.2
0.4 0.6 Time ratio t / T
0.8
1
Fig. 4. Case 2: contact forces calculated from the true and identified vehicle parameters using different population sizes: —— true; – – – size = 5; — — size=30; — – — size=60.
Starting from the initial search domains M1A(0, 10 000) kg, M2A(0, 60 000) kg, cA(0, 200 000) Ns/m and kA(0.2 107) N/m, the proposed identification method is applied with eight different kinds of population sizes, namely 5, 10, 20, 30, 40, 60, 80 and 100. The results are summarized in Table 4. First and foremost, no matter what population size is used, extremely good results are obtained with errors less than 1%. For population sizes larger than 20, results of comparable accuracy are obtained after the second stage. The computation time required for the larger population sizes of 40, 60, 80 and 100 is of course much longer than that for a smaller population size. However for a population size of 5, at least four stages are needed to get results comparable to those obtained from two stages with population sizes larger than 20. In other words, a population size that is too small reduces the speed by requiring more optimization stages, while a population size that is too large increases the computation time in each optimization stage. A population size of 30 is considered appropriate to the parameter identification of a single moving vehicle.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
105
Table 5 Case 3: effects of measurement noise on accuracy (station at x=60 m, population size=30) Noise level (%)
5 10 15 20 25
Number of stages
2 3 3 3 3
Identified values
Errors (%)
M1 (kg)
M2 (kg)
c (Ns/m)
k (N/m)
M1
M2
c
k
4775 4812 4774 4783 4762
27 042 27 136 26 962 26 783 27 185
85 668 86 474 85 855 85 491 86 006
9 107 513 9 113 098 9 162 722 9 127 093 9 147 911
0.52 0.25 0.54 0.35 0.79
0.15 0.50 0.14 0.80 0.69
0.39 0.55 0.17 0.59 0.01
0.14 0.08 0.47 0.08 0.31
325
Contact force (kN)
320
315
310
305
300
0
0.1
0.2
0.3
0.4
0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 5. Case 3: contact forces calculated from the true and identified vehicle parameters using measurements of different noise levels: —— true; – – – 5% noise; — — 10% noise.
The histories of contact forces calculated from the identified parameters using different population sizes are plotted in Fig. 4 and compared with the true values. Very good agreement is observed. 4.4. Case study 3: effect of measurement noise The effect of measurement noise on the accuracy and efficiency of the proposed method is next investigated. This case is essentially the same as Case 1 except that five measurement noise levels of 5%, 10%, 15%, 20% and 25% are introduced. Following the conclusions of the studies for Case 2, a population size of 30 is adopted for this identification problem with a single moving vehicle. Table 5 shows the identification results, from which it can be seen that the identification values vary little with the increase of the measurement noise level. This demonstrates that the proposed identification procedure is not particularly sensitive to measurement noise levels. For the five noise levels, the errors are still well below 1%. However the more measurement noise there is, the more optimization stages are required. Fig. 5 shows the identified contact forces based on noise-free measurements, and measurements with noise levels of 5% and 10%. Basically, very
ARTICLE IN PRESS 106
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
Table 6 Case 4: arrangement of measurement stations Numbers of measuring stations
Locations of measuring stations, x (m)
1 2 3 4 6
60 18, 18, 18, 12,
60 60, 102 52, 68, 102 24, 52, 68, 96, 108
Table 7 Case 4: effects of number of measurement stations Number of measurement stations Number of stages Identified values
Errors (%)
M1 (kg) M2 (kg) c (Ns/m) k (N/m)
M1
M2
c
k
5% measurement noise, population size=30 1 2 2 2 3 2 4 2 6 2
4775 4813 4808 4814 4781
27 042 27 019 26 982 26 875 26 967
85 668 86 424 86 149 85 887 85 777
9 107 513 9 110 873 9 106 786 9 087 109 9 116 734
0.52 0.27 0.17 0.30 0.40
0.15 0.07 0.07 0.09 0.12
0.39 0.49 0.17 0.07 0.26
0.14 0.10 0.15 0.36 0.04
10% measurement noise, population size=30 1 3 2 3 3 3 4 2 6 2
4812 4793 4827 4773 4799
27 136 27 038 26 882 27 012 26 962
86 474 86 464 85 798 85 767 85 777
9 113 098 9 110 588 9 097 493 9 099 936 9 115 070
0.25 0.15 0.56 0.56 0.02
0.50 0.14 0.44 0.04 0.14
0.55 0.54 0.23 0.27 0.26
0.08 0.10 0.25 0.22 0.05
20% measurement noise, population size=30 1 3 2 3 3 3 4 2 6 2
4783 4844 4821 4819 4814
26 783 26 948 27 024 27 106 27 011
85 491 86 114 85 743 85 862 85 777
9 127 093 9 088 694 9 106 786 9 123 056 9 092 955
0.35 0.92 0.44 0.39 0.30
0.80 0.19 0.09 0.39 0.04
0.59 0.13 0.30 0.16 0.26
0.08 0.34 0.15 0.03 0.30
good agreement is observed although the presence of noise in the measurements does give rise to certain discrepancies. 4.5. Case study 4: effect of number of measurement stations Very often more than one accelerometer is used in field measurements. It not only improves the accuracy in most cases, but also caters for the possibility of occasional breakdown of certain instruments. The case of an increased number of measurement stations is studied. For comparison, the identification problem is repeated using five patterns of measurement stations, namely, 1, 2, 3, 4 and 6 measurement stations, respectively. The distribution and locations of the
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
107
325
Contact force (kN)
320
315
310
305
300
0
0.1
0.2
0.3
0.4 0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 6. Case 4: contact forces calculated from the true and identified vehicle parameters using different number of measurement stations: —— true; – – – One station; — — three stations. Table 8 Case 5: assume true vehicle parameters Vehicle
Unsprung mass (kg)
Sprung mass (kg)
Damping coefficient (Ns/m)
Spring stiffness (N/m)
1 2 3
4800 3900 5400
27 000 24 000 32 000
86 000 72 000 100 000
9.12 106 8.68 106 10.0 106
measurement stations for each arrangement are summarized in Table 6. A population size of 30 is used and three levels of measurement noise, namely 5%, 10% and 20%, are considered. The identification results as shown in Table 7 are very good. The errors are well below 1% even for 20% measurement noise. Similar to what is seen before, the higher the measurement noise, the more stages are needed. With the same level of measurement noise, the efficiency of optimization in general increases with the number of measurement stations, as seen in the number of stages required. However this is achieved at the expense of computation time. The time histories of the contact forces calculated from the true and identified values of the vehicle parameters using 1 and 3 measurement stations with 10% measurement noise are plotted in Fig. 6. It is obvious that the time histories of contact forces calculated from identified parameters both agree well with the true time histories, but that identified based on three measurement stations is closer to the true one. 4.6. Case study 5: effect of number of vehicles In real life, there are always many vehicles on a bridge at the same time. In order to test the accuracy and efficiency of the proposed method for more practical situations, a convoy comprising three vehicles at a spacing of 10 m is considered and their true parameters are shown in Table 8. Both the two and six-station measurement arrangements are used, and their locations are shown in Table 6. The measurement noise for acceleration simulation is taken to be 10%. The
ARTICLE IN PRESS 108
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
Table 9 Case 5: effect of number of vehicles (10% measurement noise, population size=60) Parameters
Identified values
Errors (%)
Number of stages
Two stations
six stations
Two stations
Six stations
Vehicle 1
M11 (kg) M12 (kg) c1 (Ns/m) k1 (N/m)
4757 27266 85790 9237977
4807 27065 86196 9131049
0.90 0.99 0.24 1.29
0.15 0.24 0.23 0.12
Vehicle 2
M21 (kg) M22 (kg) c2 (Ns/m) k2 (N/m)
3861 23988 73212 8711342
3881 23979 71482 8709654
1.00 0.00 1.68 0.36
0.49 0.09 0.72 0.34
M31 (kg) M32 (kg) c3 (Ns/m) k3 (N/m)
5465 32276 101615 10068262
5429 31901 100165 10062176
1.20 0.86 1.62 0.68
0.54 0.31 0.17 0.62
Vehicles 3
4
340
Contact force (kN)
330
320
310
300
290
0
0.1
0.2
0.3
0.4 0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 7. Case 5: contact forces of the first vehicle calculated from the true and identified vehicle parameters using two measuring stations: —— true; – – – identified.
same initial search domains as in Case 1 are used. In view of the larger number of parameters to be identified in this case, the population size is set as 60. The identification results as shown in Table 9 are very good. For example, in the two-station measurement arrangement, identification errors below 2% can be obtained after only four optimization stages even though there are 12 parameters to be identified. In the six-station measurement arrangement, identification error below 1% can be obtained after only four optimization stages. It demonstrates that the proposed method can also cope with a convoy of vehicles with reasonable accuracy although longer computation time is required for more vehicles. It is also seen that the accuracy improves when more measurement stations are used. For the present three-span continuous bridge, six measurement stations are sufficient for parameter identification of the vehicle convoy.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
109
300
Contact force (kN)
290
280
270
260
250
0
0.1
0.2
0.3
0.4 0.5 0.6 Time ratio t / T
0.7
0.8
0.9
1
Fig. 8. Case 5: contact forces of the second vehicle calculated from the true and identified vehicle parameters using two measuring stations: —— true; – – identified.
410
Contact force (kN)
395 380 365 350 335 320
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time ratio t / T
Fig. 9. Case 5: contact forces of the third vehicle calculated from the true and identified vehicle parameters using two measuring stations: —— true; – – – identified.
The time histories of contact forces of the three vehicles calculated from the true and identified parameters using two measurement stations are plotted in Figs. 7–9, respectively. Note that the total traversing time T is counted from the instant when the first vehicle arrives at the starting end of the bridge to that when the last vehicle leaves the bridge. Therefore the time history for each vehicle only occupies part of the total traversing time. Good agreement is again observed.
5. Discussions In Case 1 when noise-free simulated accelerations are used, effectively the exact results are obtained. It shows that the proposed procedure is theoretically correct and accurate. Then in the other cases, the procedure is tested on various aspects including the population size, noise level of
ARTICLE IN PRESS 110
F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
acceleration measurements, number of measurement stations, and number of vehicles. Examination of the results confirms that the identification method is very robust. Even for a larger number of variables and a reasonably high noise level, very acceptable results can be obtained after only a few optimization stages with an appropriate population size. In all the five cases studied above, even only a few measurement stations are chosen along the whole bridge, very acceptable results can be obtained no matter how many vehicles are on the bridge at the same time.
6. Conclusions A procedure for the identification of vehicles moving on multi-span continuous bridges is described. The moving vehicles are modelled as 2-d.o.f. systems each comprising an unsprung mass and a sprung mass interconnected by a damper and a spring. The aim of the study is to identify the associated vehicle parameters based on the acceleration measurements at selected stations. Here the acceleration measurements are simulated by the numerical solution to the associated forward problem with random noise introduced to account for the likely errors in practice. The search for the vehicle parameters is formulated as an optimization problem in which the error between the measured accelerations and the accelerations reconstructed from trial vehicle parameters is minimized. A robust multi-stage optimization scheme based on genetic algorithms has been proposed in which the search domains are systematically reduced to converge on the optimal solution. A comprehensive study has been carried out to test the proposed procedure. It shows that in spite of the presence of measurement noise, accurate results are still obtained. The method can also cope with the presence of multiple vehicles and the use of instruments at multiple sections.
Acknowledgements The work described in this paper has been supported by the Committee on Research and Conference Grants, The University of Hong Kong, China. The present study has made use of the code of a genetic algorithm written by Dr. D.L. Carroll at the University of Illinois, and this is gratefully acknowledged.
References [1] C.A. Julien, H.M. Hung, Technique in system identification for dynamic mechanical systems, SAE-Paper 690497, 1969, p. 9. [2] R. Schlegemilch, M. Railey, System identification comparison study, Proceedings of the 21st Annual Pittsburgh Conference Part 5 (of 5), Pittsburgh, PA, May 3–4, 1990, pp. 2269–2273. [3] Y. Kyongsu, K. Hedrick, Observer-based identification of nonlinear system parameters, Journal of Dynamic Systems, Measurement and Control, Transactions of the American Society of Mechanical Engineers 117 (2) (1995) 175–182.
ARTICLE IN PRESS F.T.K. Au et al. / Journal of Sound and Vibration 269 (2004) 91–111
111
[4] K. Wise, R.E. Reid, Modeling and identification of a light truck engine mounting system for ride quality optimization, International Congress & Exposition—Society of Automotive Engineers, Warrendale, PA, 1984, 841142P. [5] D.A. Derradji, N. Mort, Multivariable adaptive control using artificial neural networks, Proceedings of the 1996 UKACC International Conference on Control. Part 2 (of 2), Exeter, UK, September 2–5 1996, 427/2, pp. 889–893. [6] P. Michelberger, J. Bokor, A. Keresztes, P. Varlaki, Identification of a multivariable linear model for road vehicle dynamics from test data, International Journal of Vehicle Design 8 (1) (1987) 96–114. [7] K.R. Goheen, E.R. Jefferys, The application of alternative modelling techniques to ROV dynamics. Proceedings of the 1990 IEEE International Conference on Robotics and Automation, Cincinnati, OH, May 13–18 1990, pp. 1302–1309. [8] D.H. Lee, J.G. Lee, Comparison of parameter identification algorithms for flight vehicles, Transactions of the Japan Society for Aeronautical and Space Sciences 36 (114) (1994) 249–256. [9] E. Bisimis, Use of system identification methods for the determination of vehicle parameters from measured time histories, Vehicle System Dynamics 6 (1977) 182–184. [10] O. Kropac, J. Sprinc, Identification of the system vehicle-road parameters, Vehicle System Dynamics 11 (4) (1982) 241–249. [11] L.R. Ray, Nonlinear tire force estimation and road friction identification: simulation and experiments, Automatica 33 (10) (1997) 1819–1833. [12] D.Y. Zheng, Y.K. Cheung, F.T.K. Au, Y.S. Cheng, Vibration of multi-span non-uniform beams under moving loads by using modified beam vibration functions, Journal of Sound and Vibration 212 (3) (1998) 455–467. [13] Y.K. Cheung, F.T.K. Au, D.Y. Zheng, Y.S. Cheng, Vibration of multi-span non-uniform bridges under moving vehicles and trains by using modified beam vibration functions, Journal of Sound and Vibration 228 (3) (1999) 611–628. [14] R.J. Jiang, F.T.K. Au, Y.K. Cheung, Identification of masses moving on multi-span beams based on genetic algorithms, Computers and Structures 81 (22–23) (2003) 2137–2148. [15] D.E. Goldberg, Genetic Algorithms in Search Optimization, and Machine Learning, Addison-Wesley Reading, MA, 1989.