Finite element model updating using simulated annealing hybridized with unscented Kalman filter

Computers and Structures 177 (2016) 176–191

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Finite element model updating using simulated annealing hybridized with unscented Kalman filter Rodrigo Astroza a,⇑, Luan T. Nguyen b, Tamara Nestorovic´ b a b

Faculty of Engineering and Applied Sciences, Universidad de los Andes, Santiago, Chile Mechanics of Adaptive Systems, Department of Civil and Environmental Engineering, Ruhr-Universität Bochum, Germany

a r t i c l e

i n f o

Article history: Received 2 June 2016 2 August 2016 Accepted 1 September 2016

Keywords: Finite element model Model updating Simulated annealing Unscented Kalman filter Parameter estimation Global optimization

a b s t r a c t This paper proposes a method for finite element (FE) model updating of civil structures. The method is a hybrid global optimization algorithm combining simulated annealing (SA) with the unscented Kalman filter (UKF). The objective function in the optimization problem can be defined in the modal, time, or frequency domains. The algorithm improves the accuracy, convergence rate, and computational cost of the SA algorithm by local improvements of the accepted candidates though the UKF. The proposed methodology is validated using a mathematical function and numerically simulated response data from linear and nonlinear FE models of realistic three-dimensional structures. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Finite element model updating (FEMU) can be defined as the process of tuning a finite element (FE) model to minimize the discrepancy between the measured and FE predicted responses of the structure being modeled [18]. The tuning process is usually conducted in an off-line fashion by using batch processing methods and consists of seeking inaccurate or unknown parameters of the FE model assuming that the model structure is fixed. FEMU has attracted significant attention from the structural engineering community because of its applications in structural dynamics, mainly in damage identification (DID) and response prediction (e.g., [19,47,25,42,15,43]). One of the most popular approaches for DID makes use of a linear FE model which is calibrated using low amplitude vibration data recorded before and after the structure has suffered damage. Then, damage is identified as the reduction of effective stiffness in one or more regions of the structure. Also, approaches for DID based on nonlinear model updating have been proposed (e.g., [13,44,51]), however, all these studies have used simplified nonlinear models with lumped nonlinearities defined phenomenologically. Only in recent years

⇑ Corresponding author at: Faculty of Engineering and Applied Sciences, Universidad de los Andes, Mons. Álvaro del Portillo 12455, Las Condes, Santiago, Chile. E-mail address: [email protected] (R. Astroza). http://dx.doi.org/10.1016/j.compstruc.2016.09.001 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

high-fidelity mechanics-based structural nonlinear FE models, which are used for analysis and design, have been employed for DID (e.g., [3,14]). The use of advanced nonlinear structural FE models allows to describe the presence, location, type, and extent of damage in the structure, because these models provide valuable information about history of plastic deformations, residual deformations, loss of strength, and loss of ductility capacity. Methods for constrained nonlinear optimization are typically used to solve the inverse problem of FEMU considering an objective function describing the discrepancies between the FE predicted and measured responses or quantities derived therefrom (e.g., modal parameters or frequency response functions). Because of the complexity of the relationship between the model parameters to be identified and the objective function, the latter can include many local minima. Gradient-based methods are commonly used for FEMU (e.g., [48,6,37]), however they might be trapped in local minima and their solution highly depends on the starting point (i.e., initial guess of the model parameters). In order to avoid this issue, different global optimization algorithms (GOAs) have been used for FEMU. Teughels et al. [49] and Bakir et al. [7] investigated the use of coupled local minimizers (CLM) for linear FEMU using modal data. Shabbir and Omenzetter [41] combined Particle Swarm Optimization (PSO) with sequential niche technique to improve the finding of global minimum in FEMU and applied the approach to update a linear FE model of a full-scale pedestrian bridge using modal data. PSO was also used by Marwala [33] to update linear FE models of a suspended aluminum beam

177

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

and an unsymmetrical structure tested on a laboratory environment. Perera et al. [40] used PSO and modal data to update a linear FE model of a one-story one-bay reinforced concrete (RC) frame experimentally tested and compared the results with those obtained using Genetic Algorithms (GAs). Hasançebi and Dumlupinar [21] used artificial neural networks (ANNs) to updated FE models of a real RC bridge using the identified natural frequencies and measured static deflections. They employed linear and nonlinear FE models to generate datasets for network training and concluded that nonlinear models provide a better agreement between the updated FE model and the experimental data, especially for static deflection measurements. Betti et al. [9] combined ANNs and GAs to update a linear FE model of a reduced scaled three-story oneby-one bay steel frame at different damage levels using modal data. Other researchers have employed ANNs for FEMU and DID purposes (e.g., [4,16,52,33]). Hao and Xia [20] used modal data and GAs to identify damage in a cantilever beam and a twodimensional steel frame tested in laboratory conditions. Meruane and Heylen [35] used parallel GAs to identify damage in a RC beam and a simple and small aircraft prototype using model updating based on modal data. Chisari et al. [12] utilized GAs to update a linear FE model of an isolated bridge using experimentally collected static and dynamic data. Kang et al. [28] combined artificial bee colony algorithm with Nelder-Mead simplex method and applied the hybrid approach to identify the model parameters of concrete dam-foundation systems using numerically simulated data. A comprehensive review on system identification and model updating of civil structures using biologically-inspired (e.g., ANNs and GAs) methods can be found in Sirca and Adeli [43]. Methodologies based on modeling to generate alternative (MGA) have also been employed to provide multiple alternative solutions in linear FEMU (e.g., [53,10]). Simulated Annealing (SA) [29] has also been used to update FE models. Levin and Lieven [30] introduced the use of SA in FEMU and compared its performance with GAs using experimental data of a flat plate wing structure. They employed various perturbation schemes to minimize the objective function constructed based on experimental and FE simulated frequency response functions. The literature on the use of heuristic optimization methods for structural FEMU shows only a few works that have employed SA as an optimization method. The fact that SA has not received as much attention as other metaheuristic optimization methods is because the standard SA requires a large number of annealing cycles to converge to a satisfactory solution. This means that many simulations of the FE model need to be run during the sequential SA process. This issue limits the applicability of SA in large scale and physics-based complex FE models of geo-materials and structures such as tunnels, building, dams, and bridges. Zimmerman and Lynch [54] proposed the parallel SA implemented on a wireless sensor network for structural health monitoring (SHM) that reduces the runtime of the SA algorithm by breaking the annealing into a number of temperature steps, each run on a separate computer node supported with appropriate communication to other wireless nodes. Jeong and Lee [26] as well as He and Hwang [22] combined GAs with SA to create an adaptive SA-GA that can fix the poor hill-climbing capability of the GA for applications in generic system identification and damage detection, respectively. In this paper, an acceleration scheme recently proposed by some of the authors to reduce annealing steps of the SA by local improvements of the accepted candidates is employed to solve structural FEMU problems. The hybridized scheme has proved effective for a parameter identification problem of waveform inversion of the underground tunnel seismic waves in Nguyen and Nestorovic´ [39]. Similar to the principle in Martin and Otto [32], a method for improving an intermediate candidate locally is coupled to the main SA algorithm. However, in contrast to the idea

in Martin and Otto [32], the local improvement method employed in this work is run only if the suggested candidate is accepted following the Metropolis’ rule, thus reducing the number of FE model evaluations. The combined global-local optimization method presented in this work couples SA algorithm with the unscented Kalman filter (UKF). The UKF [27], which belongs to the sigma-point Kalman filter family, is a derivative-free estimator for nonlinear state-space models that can be efficiently adapted for solving parameter identification of static and dynamic models [50]. As shown later in the paper, the choice of the UKF as a local minimization method is based on the observation that the UKF can be very efficient to converge to a nearby local minimum. Application examples consisting of a test function (the Styblinski–Tang function) with six parameters and linear and nonlinear FE models are presented to verify the performance of the proposed method. For the FEMU examples, inverse problems involving different number of model parameters are analyzed and modal and time history data simulated numerically from realistic three-dimensional steel frame structures are used in the objective functions. Accurate results are obtained in the different application examples even when a limited number of response quantities (modal properties or time history responses) are used to define the objective function. The hybrid scheme improves the accuracy and reduces the computational cost of the standard SA algorithm, proving that its application to realistic civil structures is feasible. The effects of using modal parameters, acceleration time history responses, and heterogeneous time history response quantities on the estimation results are investigated. Based on the analyses presented in the paper, values for the parameters involved in the hybrid scheme are proposed to solve FEMU problems. It is noted that properly calibrated nonlinear FE models can be used to identify potential damage in the structure, providing an excellent tool for SHM. It is worth nothing that the hybrid scheme also has the capability to be employed with surrogate models. 2. Simulated annealing combined with the unscented Kalman filter 2.1. Forward model and the objective function The forward model that relates hidden model parameters to responses, for both linear and nonlinear finite element models, can be generally represented as

d ¼ hðmÞ þ v ;

ð1Þ

where vector m contains n hidden parameters of the model, d stores the model outputs obtained from the computational model hðÞ. Modeling errors v are caused by assumptions made in building the mathematical model and the numerical approximations involved. Throughout this work, v is assumed to be a stationary, zero-mean Gaussian white vector process with covariance Rm . Measuring the real life structural responses for the corresponding obs

observed quantities results in a set of measurement data d that are intrinsically contaminated by measurement noise. Measurement noise is often satisfactorily represented as having normal distribution with zero-mean and covariance Robs . Then, the total observation error covariance is formed by the addition rule R ¼ Robs þ Rm [45]. In this work the effects of modeling uncertainty are not considered since the same models are employed in the simulation and estimation phases. Thus, the total observation error covariance consists of the measurement noise covariance only, i.e., R ¼ Robs . In this work, different definitions to build the objective function that needs to be minimized are used. For the proposed algorithm, it is desired that the L2 norm of the residual is broken down into r

178

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

component-wise objectives s ¼ ðs1 ; s2 ; . . . ; sr Þ with obs

sj ¼ kdj

 dj k2

for j ¼ 1; . . . ; r

to

form

a

residual

vector

ð2Þ

For the FEMU problem tackled in this work, the residual components are based on the number of modes or the number of observation points depending on if the modal properties or time histories are used as estimated/observed quantities. The single scalar objective is then simply the summation of all residual components

S¼

r X sj

ð3Þ

j¼1

It is noted that the component-wise objective vector, Eq. (2), is used in place of the observation vector in the UKF iterations and the total objective, Eq. (3), is used to represent the energy in the simulated annealing process. In the definition of the total objective, different weighting coefficients can be employed for different component-wise objectives to take into account different resolution or accuracy levels of different response quantities. 2.2. SA algorithm SA is a probabilistic algorithm that aims to find near-optimal solutions of a given energy function to be minimized. The core of SA is the Metropolis sampling technique [36], which generates random samples from a proposal probability distribution function in a high dimensional space whose underlying mathematical description is unknown but its function values can be evaluated. The Metropolis procedure is based on the energy change (DE) of the system due to a random perturbation of the current configuration to decide if the proposed solution is accepted or rejected. Probability of acceptance follows the Boltzmann distribution   PðDEÞ ¼ exp  KDB ET c , in which K B = Boltzmann constant and T c = current temperature of the system. If DE 6 0, the proposal is accepted, otherwise it is accepted with probability PðDEÞ. The probabilistic acceptance of the worse candidate, which has higher energy than the current configuration, allows the algorithm to escape local minima and to explore other possible lower energy configurations. Mimicking the ‘physical’ annealing process of metals and glass, a cooling schedule is introduced to SA algorithm to regulate the acceptance rates along the ‘simulated’ annealing process. Typically, high temperature is set at the start of SA process and is decreased gradually to a very low temperature at the end. Accordingly, the acceptance rates are high in early annealing cycles to allow for more freely jumps of the states (configurations). As the annealing step increases, the corresponding acceptance rate decreases and becomes very low at the end, when the global minimum configuration is expected to be found, so that the random moves keep stabilized at the minimum energy configuration. In this paper, the energy function to minimize is the objective function defined in Eq. (3). The objective function can take values of any magnitude depending on the scale and the type of measurement data. To obtain an acceptance probability independent of the magnitude of the misfit, a normalization of the change of the objective value as proposed in Balling [8] is adopted here

  DS exp  DSa T c

ð4Þ

where DS is a change in objective function’s value of the proposed move compared with that from the previous accepted move and T c is the current temperature. The Boltzmann constant (K B ) is set to be the average DSa of the changes in objective function’s values of all the accepted moves up to the current annealing cycle. A fast

annealing schedule, in which the temperature is set to decrease in an exponential rate with a base a 2 ð0; 1Þ, is used in the implementation of SA. Temperature in the next annealing cycle c þ 1 is related to that in the current annealing cycle c by the following relation

T cþ1 ¼ aT c

ð5Þ

It is noted that the cooling rate a can be directly determined given the initial and end temperatures, T 0 and T N respectively, by the relation a ¼ ðT N =T 0 ÞNc 1 , with N c the number of annealing cycles. Random candidates are generated according to Ingber [24] under consideration of the upper bound U mi and the lower bound Lmi of each parameter mi , where i ¼ 1; . . . ; n. 1

mcþ1 ¼ mci þ yi ðU mi  Lmi Þ i

ð6Þ

where yi 2 ½1; 1 is generated by the uniform distribution variable ui 2 ½0; 1 as

# j2ui 1j   " 1 1 Tc 1 þ yi ¼ sgn ui  1 2 Tc

ð7Þ

It is noted that the random sampling Eq. (6) is repeated until 2 ½Lmi ; U mi , for i ¼ 1; . . . ; n. mcþ1 i 2.3. The UKF The UKF has two prominent advantages over the extended Kalman filter (EKF). First, statistical linearization used for the UKF suffers less from linearization error than analytical linearization (Taylor-series expansion) used in the EKF. Second, calculation of the gradients, which is easily prone to inaccuracy for a highly nonlinear function and can be impossible if the function is nonsmooth, is not required for the UKF. As such, the UKF has been used widely for state and parameter estimation applications in many disciplines. In the majority of works in structural parameter identification, the UKF is used for state or dual state-parameter estimation [31,5,11,1]. The UKF is used in this work for parameter-only identification to improve the SA random proposed candidates locally by minimizing an objective function. The use of the UKF for pure parameter estimation and minimization problems has been shown effective in Van Der Merwe and Wan [50], Astroza et al. [3], Nguyen and Nestorovic´ [38], etc. After the measurements are taken into account in the filtering (correction) step ðk  1Þ, the posterior mean and covariance are ^ þ and Pm denoted as m þ , respectively. In the next iterative step k, the posterior mean and covariance from the previous iterative step ^  and Pm are used to predict the prior mean and covariance m  by the use of the unscented transformation (UT) and the weighted average rule. This filtering procedure is repeated for a number of iterations. For details on the use of the UKF for FEMU, the reader can refer to Astroza et al. [3] and Nguyen and Nestorovic´ [38]. Here only a summary of the UT and a flow diagram of the UKF procedure are shown. With n the number of hidden model parameters, the UT requires ð2n þ 1Þ deterministic sample points, so-called sigmapoints (SPs), which are assigned as

^þ M0 ¼ m

ð8Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ^þþ ðn þ kÞPm ; Mi ¼ m þ

for i ¼ 1; . . . ; n

ð9Þ

i

^þ Mnþi ¼ m

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðn þ kÞPm ; þ

for i ¼ 1; . . . ; n

ð10Þ

i

The parameter k can be used to adjust the spread of the SPs about the predicted mean estimate. The notation ðÞi denotes the i-column of the matrix within parentheses. Each SP is associated

179

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

with a weight. These weights are defined in Eq. (11) and their sum is equal to one.

W0 ¼

k 1 ; W i ¼ W iþn ¼ nþk 2ðn þ kÞ

for i ¼ 1; . . . ; n

ð11Þ

The mapping of the SPs through the objective function (Eq. (2)) calculated for each of the r multiple residuals s ¼ ðs1 ; s2 ; . . . ; sr Þ; s : Rn ! Rr , is done for each SP

S i ¼ sðMi Þ for i ¼ 0; . . . ; 2n

SA accepted move

UKF

1 1*

ð12Þ

This nonlinear mapping of the SPs helps to preserve second order accuracy of the Gaussian distribution as opposed to the first order approximation involved in the EKF. The flowchart of the implemented UKF procedure is shown in Fig. 1. 2.4. The hybrid method The proposed hybrid algorithm is a global-local optimization scheme that combines the merits of both optimization methods – SA algorithm and the UKF – to achieve the global minimum solution at reduced computation time. The mechanism for achieving fast and accurate convergence of the combined algorithm is illustrated in Fig. 2. As can be seen in Fig. 2, the UKF is executed whenever random candidate solutions (represented as 1 and 2 in Fig. 2) are accepted (according to the Metropolis sampling). The execution of the UKF in the whole SA cycles is designed in this manner to drive the accepted candidate solutions (1 and 2) proposed randomly to positions having lower objective values (1⁄ and 2⁄). It is noted that the improved solution in local minimum area (1⁄) is

UKF

2

2* Fig. 2. Global-local optimization concept of the new hybridized algorithm (Adapted from Nguyen and Nestorovic´ [39]).

more biased than that in the global minimum area (2⁄) because the residual at the local minima is much larger than that at the global minimum. The acceptance rates in later annealing cycles become very small if near-global or global solution is found owing to the Metropolis sampling that is used in the SA algorithm. The socalled unscented hybrid simulated annealing (UHSA) algorithm [39] is based on early convergence to the minima to achieve fast convergence to the global minimum. Thus, applying the UHSA algorithm to solve inverse problems does not require very long annealing cycles compared with the standard SA. Fig. 3 shows the complete flowchart of the proposed UHSA algorithm. The proposed hybrid algorithm is derivative-free, therefore it does not require computation of Jacobian matrices. Although it is simple to implement, the statistical convergence of the SA is assured while the convergence speed and accuracy of the solution

Fig. 1. The UKF procedure used for local minimization.

180

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

Start

Configure settings

c =1

Calculate crosscovariance

Calculate Kalman gain

Predict model responses

Correct estimates

Propose a new candidate

Calculate acceptance propability

c=c+1

Predict estimates

Y k =1

Generate and transfer the sigma-points

Y k = k +1

k < Nk ?

N

UKF

Store intermediate results

Lower temperature

Y

N

c < Nc ?

N

SA

End Fig. 3. Proposed hybridized algorithm composed of SA (left block) and the UKF (right block).

are improved [39]. As for the standard SA, the UHSA requires a starting point to set its current energy state. It is suggested to set the cooling rate (a) very close to one and configure the algorithm parameters related to the UKF such that the acceptance rates are kept as high as between 20% and 25%. 3. Convergence behavior with test function 3.1. The multi-dimensional Styblinski-Tang function The multi-dimensional Styblinski-Tang function is chosen to test the convergence of the proposed UHSA algorithm and to compare its convergence behavior to that of the standard SA. To have the non-negative function, which is the case of general objective functions, a constant equal to the global minimum value 39:166 n (with n the number of dimensions) is subtracted from the original function to form

f ðm1 ; m2 ; . . . ; mn Þ ¼

n 1X ðm4  16m2i þ 5mi Þ  ð39:166 nÞ 2 i¼1 i

of dimensions increases. As such, to solve this problem properly, a global optimization method is preferable. 3.2. Convergence compared with standard SA The Styblinski-Tang function in the 6-dimensional parameter space limited from 5 to 5 in each dimension is subjected to be minimized. The UHSA is set to run N c ¼ 2000 annealing cycles with initial temperature T 0 ¼ 6 and final temperature T N ¼ 0:01. The UKF is configured with N k ¼ 8 iterations. The measurement error is set to be approximately 2% deviated from the function value

ð13Þ

The function in Eq. (13) has global minimum f min  0 at point mmin ¼ ð2:903; 2:903; . . . ; 2:903Þ. Fig. 4 plots an elevated surface of the bivariate function (n ¼ 2) described by Eq. (13). The projected contour curves show that the test function has a global minimum and three local minima located at far distances. The problem might look not very challenging because a gradient-based linear optimization method may still have some chance to converge to the global minimum if the starting point is accidentally chosen close to the global minimum. However, the chance of choosing an appropriate starting point decreases significantly as the number

Fig. 4. Surface and contour plots of the bivariate multi-minima Styblinski-Tang test function.

181

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

chosen is that the objective becomes less than 2.5% of the initial objective value. It is observed that in many cases the UHSA needs fewer numbers of function evaluations to attain the same convergence criterion. As shown in Fig. 6, in the best case the UHSA can achieve one order of magnitude less number of function evaluations than that required by the best standard SA. Even in the worst case of UHSA run, the number of function evaluations needed is still slightly less than that of the best standard SA run. As for the stochastic global optimization, not only is important that the global minimum is found, but the distribution of each parameter over the predefined feasible range also provides an insight into the multi-minima objective surface as well as an estimate of how influential each parameter is on that surface. Fig. 7 plots the histograms of such distribution for the test problem considered. The frequencies of the accepted candidates show that aside from the global minimal region, where high concentration of the accepted candidates locates, there exists one local minimal region with much scarcer appearance of the accepted candidates at a far distance for each parameter.

evaluated at the randomly chosen initial point ^ 0 ¼ ð1; 1; 1; 1; 1; 1Þ. The initial estimation error covariance Pm m 0 is set to have diagonal structure whose standard deviations are of one fortieth of the feasible ranges. The covariance matrix Q is set to be a small fraction of Pm 0 . In particular, the UHSA is set numerically as below:  Settings of SA: a ¼ 0:9968, Lm ¼ ð5; 5; 5; 5; 5; 5Þ, ^ 0 ¼ ð1; 1; 1; 1; 1; 1Þ Um ¼ ð5; 5; 5; 5; 5; 5Þ, m  m m 2  U L  Settings of the UKF: N k ¼ 8, Pm , R ¼ 4:02 , 0 ¼ diag 40 Q ¼ 0:12 Pm 0 , k ¼ 0:001 In order to compare the performance of the UHSA to that of the standard SA (without UKF local improvement), the latter is set to run N c ¼ 30; 000 annealing cycles. In case of standard SA, the same implementation as the UHSA is used, but with number of UKF iterations N k ¼ 0. Five pairs of annealing processes (five standard SA runs and five UHSA runs), each with the same random seed, are run for a comparative study. Fig. 5 shows that for all the random runs performed, the UHSA needs fewer annealing cycles than the standard SA to attain satisfactory small objective values. If it is aimed to estimate the best single solution, the annealing process can be set to stop when a certain residual is met in order to save the runtime. In this example, the convergence criterion

10

3.3. Influence of the algorithm parameters Assuming that the settings for SA algorithm are standard, it is aimed to study how the UKF influences the convergence behavior of the hybrid optimization method.

3

Best objective

UHSA Standard SA 10

2

10

1

10

0

10

−1

10

0

10

1

10

2

10

3

10

4

Annealing cycle

Number of function evaluations

Fig. 5. Best objectives along annealing cycles for five pairs of random annealing processes.

10

5

10

4

10

3

10

2

10

1

10

0

UHSA Standard SA 10

0

10

1

10

2

10

3

10

4

Annealing cycle Fig. 6. Number of function evaluations to achieve convergence. Each red circle indicates the cycle at which convergence criterion is met for each annealing process. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

182

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

Accepted candidates

True value

Best estimate

2000

m1

1000 0 2000

m2

1000 0

Frequency

2000

m3

1000 0 2000

m4

1000 0 1000

m5

500 0 1000

m6

500 0 −5

−4

−3

−2

−1

0

1

2

3

4

5

m Fig. 7. Histograms resulting from a UHSA run.

Number of function evaluations

3.3.1. Influence of the covariance matrices Covariance matrices R, Q , and Pm 0 together decide the convergence path of an UKF iteration. The initial covariance Pm 0 needs not be close to the true estimation error covariance but it is preferably set small enough to reduce too big jumps in the UKF iterations because the UKF is aimed only at attaining the local minima in the vicinity of the randomly accepted candidates. Practically, Pm 0 is a diagonal matrix whose square roots of the elements are set inversely proportional to the degree of multimodality of the objective landscape, but typically not exceeding one tenth of the range for each dimension. Elements in the covariance matrix Q can be set so that their standard deviations should be around ten percent the corresponding standard deviations in Pm 0 . Data uncertainty matrix R is assigned according to the noise level of the data recorded at the receivers. It is noted that in general, R includes the effects of measurement and modeling uncertainty, however, modeling uncertainty is out of the scope of this paper.

10

6

10

4

10

2

3.3.2. Influence of the number of UKF iterations Because the UT required by the UKF needs ð2n þ 1ÞN k function evaluations, it is important to set a proper number of UKF iterations (N k ). Five random annealing processes with N k being set to 4, 8, and 24 are run. The results plotted in Fig. 8 show that neither a very small nor a very large number of UKF iterations is good. Very small number of UKF iterations (N k ¼ 4) generally requires less number of function evaluations but causes the UHSA to meet the convergence criterion at very late annealing cycles. Very large UKF iterations (N k ¼ 24) does not help the hybrid algorithm to achieve better results either in the number of function evaluations or the annealing cycles needed. 4. FEMU: verification examples The proposed algorithm is applied to parameter estimation and model updating of mechanics-based linear and nonlinear FE

UHSA − 4UKF UHSA − 8UKF UHSA − 24UKF 10

0

10

0

10

1

10

2

10

3

Annealing cycle Fig. 8. Effect of the number of UKF iterations N k on the number of function evaluations and annealing cycles required to achieve convergence (denoted by circles).

183

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

models of steel building structures using numerically simulated data. A given FE model is initially used to simulate the response of the structure and then some parameters of the FE model are assumed unknown (representing parameters with values unknown or highly uncertain in practice) and estimated using the presented hybrid method. All the FE models used in this study were developed in the software OpenSees [34]. It is noted that in all the application examples the same FE model is used for response simulation and parameter estimation, i.e., the effects of model uncertainty are not considered in this study. The FEMU were run in a desktop workstation with an Intel Xeon 2.66 GHz processor and 48 GB RAM. 4.1. Linear finite element model updating A three-dimensional (3D) 5-story 3-by-2 bay steel frame building (Fig. 9) is considered as verification example for linear FEMU. Three case studies are considered in this section. First, in order to better graphically illustrate the results, only two model parameters are considered unknown and the model is updated based on the modal properties of the building. Then, six model parameters are identified using both the modal properties and the time history response data. In the latter, the effects of output measurement noise in the estimation results are investigated. The building is designed as an intermediate moment-resisting frame according to the 2012 International Building Code [23] for a location in downtown Seattle, WA, with Site Class D soil conditions and a short-period spectral acceleration SMS = 1.37g and a one-second spectral acceleration SM1 = 0.53g. The frame has a story height of 3.5 m and bay width of 5.0 m and 6.0 m in the longitudinal (X) and transverse (Z) directions, respectively (Fig. 9). The columns located on axis X = 0 m and X = 15 m have W14  61 cross

section, and those on axis X = 5 m and X = 10 m have cross section W14  90. Longitudinal beams on second and third levels have W21  62 cross section, while on fourth, fifth, and roof levels have W21  55 cross section. Transverse beams on second and third levels have W18  40 cross section, while on fourth, fifth, and roof levels have W18  35 cross section. Columns and beams are made of A992 and A36 steel, respectively. Fig. 10 depicts the FE model of the linear steel frame. Modal masses, gravity loads acting on beams, and geometry are shown. Each beam and column member of the frame is modeled with a single Bernoulli–Euler beam element. The properties of the cross-sections (area, bending moments of inertia, and torsional moment of inertia) are taken as the nominal values and the Young’s and shear modulus of the material are assumed as Es ¼ 200 GPa and Gs ¼ Es =½2ð1 þ mÞ, where m ¼ 0:3 is the Poisson’s ratio for steel. 4.1.1. Verification examples 1 and 2: linear FEMU using modal data The natural frequencies and mode shapes of the first r modes of the frame are used to define the objective function

SðmÞ ¼

2 2 r X ðxID;i Þ  ðxFE;i ðmÞÞ þ ð1  MACð/ID;i ; /FE;i ðmÞÞÞ ðxID;i Þ2 i¼1 ð14Þ

where k  k denotes absolute value, xID;i = identified natural frequency of mode ith, xFE;i ðmÞ = natural frequency of mode ith of the FE model, MACð/p ; /q Þ is the modal assurance criterion [2] between mode shapes /p and /q , /ID;i = ith identified mode shape, /FE;i ðmÞ = ith mode shape of the FE model. It is noted that other objective functions involving natural frequencies and mode shapes can be considered to define the objective function SðmÞ when using

Y

X

W14×90

W14×61

3.5 m (typ)

Z Fig. 9. Steel frame building used for verification examples of linear FE models.

184

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

Fig. 10. FE model of the steel frame considered in verification examples of linear models.

modal data (e.g., [18,33]). Investigation on the effects of the selection of the objective function is out of the scope of this paper. A first verification example (VE) considers r ¼ 3 and only two material parameters to be identified, which consist of the Young’s modulus of the steel of two sets of columns of the first story. The first set comprises the first story columns located on axis X = 0 and X = 5 m (marked in red1 in Fig. 9) and the second set includes the first story columns located on axis X = 15 m (marked in blue in Fig. 9). These sets of columns were chosen to have an objective function with two minima. To obtain the ‘‘true” dynamic characteristics of the structure, the following material parameter values, referred to as true values, are assumed mtrue ¼ mtrue ¼ 200 GPa. The initial esti1 2 mate of the FE model parameters is chosen as true ^ 0 ¼ ½m ^ 1;0 ; m ^ 2;0  ¼ ½0:5mtrue m 1 ; 1:5m2 , with a value of the objective function Sjm^ 0 ¼ 1:024. The feasible search domain in the estimation process is assumed to be Lm ¼ 0:3 mtrue 6 m 6 3:0 mtrue ¼ Um . In the implementation of the optimization algorithm presented in Section 2, the following parameters are chosen: T 0 ¼ 4:0 (see Eq. (5)),  m m 2  U L , k ¼ 0:0001, and T N ¼ 0:01 (see Eq. (5)), Pm 0 ¼ diag 100 Q ¼ 0:01Pm 0 . Two cases for the data uncertainty matrix (R), which is assumed to be diagonal with entries Rjj ðj ¼ 1; . . . ; rÞ, are analyzed. Table 1 reports the estimation results for VE1 considering different values of Nc , Nk , and both values of Rjj . As the number of iterations of the UKF (Nk ) increases, the estimation results (in terms of parameter estimates and value of the objective function) improve, confirming the good performance of the proposed global optimization algorithm combining SA and the UKF. For fixed values of Nc and Rjj , the estimation results improve as Nk increases because a better local search is provided by the iterations of the UKF. In all the cases the objective function is considerably reduced from Sjm^ 0 ¼ 1:024 to Sjm^  6 0:1506 and Sjm^  6 0:0104 for Rjj ¼ 1e4 and Rjj ¼ 1e10 , respectively. Better estimation results are obtained in the case of Rjj ¼ 1e10 because there is no uncertainty in the simulation and therefore a low value of Rjj in the estimation process provides a better guess of the error involved in the simulated modal data. When a wrong estimate of the diagonal entries of R is considered (Rjj ¼ 1e4 ), the proposed algorithm is robust, and good estimates of the model parameters are obtained for cases VE1_06, VE1_08, VE1_09, VE1_10, VE1_11, and VE1_12, i.e., when large enough values of Nc

1 For interpretation of color in Figs. 9 and 16, the reader is referred to the web version of this article.

and Nk are considered. In terms of computational cost, the wallclock time increases as N c and Nk increase, as expected. However, only when Nc ¼ 500 and Nk ¼ 8 the wall-clock time is higher than 10 min, confirming a feasible computational cost of the hybrid algorithm when applied to realistic civil structures. It is observed that comparing cases with similar values of Sjm^  , the computational cost is lower if Nk increases. Figs. 11 and 12 plot the estimation results for cases VE1_08 (N c ¼ 300 and N k ¼ 4) and VE1_10 (N c ¼ 500 and N k ¼ 1) with Rjj ¼ 1e4 , respectively. Both cases have similar computational cost, however a better performance is obtained when the UKF is used to improve locally the accepted candidates of the SA algorithm. A considerable smaller value of the objective function eval^  ) is obtained in case VE1_08 uated at the final estimate (m compared to VE1_10 (see Table 1). As observed from the contour plots, which depict the value of the objective function in terms of the parameter values, the accepted candidates are much more concentrated in the zone of low values of the objective function in case VE1_08 (see Fig. 11a), showing the improvement due to the iterations of the UKF. In addition, in both figures some accepted candidates are located around the local minimum ½m1 ; m2  ¼ ½0:7; 0:4 (see Figs. 11a and 12a), showing that the algorithm is able to escape from local minima due to the global search capabilities of the SA. The second VE considers six material parameters to be estimated, which are the Young’s modulus of: columns located on axis X = 0 and X = 15 m (section W14  61); columns located on axis X = 5 m and X = 10 m (section W14  90); longitudinal beams at 2nd and 3rd levels (section W21  62); longitudinal beams at 4th, 5th, and roof levels (section W21  55); transverse beams at 2nd and 3rd levels (section W18  40); transverse beams at 4th, 5th, and roof levels (section W18  35). The true model parameters are

chosen

as

T

true true true true true mtrue ¼ ½mtrue 1 ; m2 ; m3 ; m4 ; m5 ; m6  ¼ T

½200; 200; 200; 200; 200; 200 GPa and the initial estimate is true ^ 0 ¼ ½m ^ 1;0 ; m ^ 2;0 ; m ^ 3;0 ; m ^ 4;0 ; m ^ 5;0 ; m ^ 6;0 T ¼ ½0:7mtrue taken as m 1 ; 1:5m2 ; true true true T 0:5mtrue 3 ; 0:8m4 ; 1:1m5 ; 1:3m6  . The same values for the feasible search domain and parameters T 0 , T N , Pm 0 , Q , and k as in VE1 are used. Based on the results obtained for VE1, here only N c ¼ 800 and Rjj ¼ 1e10 are analyzed. Two different number of modes are considered in this verification example, r ¼ 3 and r ¼ 6. Table 2 summarizes the estimation results for VE2. The objective function decreases from Sjm^ 0 ¼ 0:34 (objective function evaluated at the initial estimate) to Sjm^  6 0:0031 implying the

185

R. Astroza et al. / Computers and Structures 177 (2016) 176–191 Table 1 Estimation results for linear FE model with two unknown parameters using modal properties. Case

Nc

Nk

VE1_01

50

1

VE1_02

50

4

VE1_03

50

8

VE1_04

100

1

VE1_05

100

4

VE1_06

100

8

VE1_07

300

1

VE1_08

300

4

VE1_09

300

8

VE1_10

500

1

VE1_11

500

4

VE1_12

500

8

^ S at m

Wall-clock time (min)

1.53 1.00 1.13 1.00 1.10 1.00

0.1506 0.0104 0.0415 0.0001 0.0314 0.0008

0.3 0.5 0.4 0.6 0.7 0.7

1.02 1.00 1.03 1.00 1.03 1.00

1.08 0.99 1.05 1.00 1.04 1.00

0.0397 0.0048 0.0169 0.0006 0.0134 0.0014

0.5 0.8 1.0 1.3 1.6 1.0

1e4 1e10 1e4 1e10 1e4 1e10

1.03 1.00 1.01 1.00 1.01 1.00

1.07 0.99 1.02 1.00 1.01 1.00

0.0256 0.0053 0.0072 4.0e6 0.0030 4.0e6

2.2 3.7 3.0 5.0 5.7 13.4

1e4 1e10 1e4 1e10 1e4 1e10

0.98 1.00 0.98 1.00 1.00 1.00

0.99 1.02 0.97 1.00 1.00 1.00

0.0275 0.0095 0.0084 0.0001 0.0006 3.0e6

3.5 6.8 8.9 7.4 9.8 12.3

Rjj

Parameter estimate ^ 1 =mtrue m 1

^ 2 =mtrue m 2

1e4 1e10 1e4 1e10 1e4 1e10

1.31 0.99 1.08 1.00 1.06 1.00

1e4 1e10 1e4 1e10 1e4 1e10

(a)

(b)

(c)

(d)

Fig. 11. Estimation results for VE1_08 with Rjj ¼ 1e4 : (a) contour plot, (b) history of the objective function, (c) convergence history of estimates of m1 , and (d) convergence history of estimates of m2 .

success update of the FE model. The wall-clock time for the four cases is around 20–30 min and the computational cost is similar if 3 or 6 modes are considered to update the FE model. Considerably better parameter estimates are obtained when r increases. In particular, it is observed that model parameters m3 and m4 (Young’s modulus of longitudinal beams) are identified with relatively high errors when r ¼ 3 but they are better estimated when r ¼ 6. This is due to the fact that the natural frequencies and mode shapes of the first three modes are not very sensitive to these parameters, however, higher modes are more sensitive to them. This can be observed by comparing the modal parameters (natural frequencies

^ 0 , and m ^  for the cases and mode shapes) obtained with mtrue , m with r ¼ 3 and r ¼ 6 (Table 3). In agreement with the results obtained for VE1, estimation results slightly improve as N k increases. Fig. 13 shows the estimation results for case VE2_01 with r ¼ 6. The objective function decreases as the annealing cycles progress and the estimates of the model parameters converge toward their true values. It is noted that standard SA (N k ¼ 0) with N c ¼ 5000 provided a worse estimation of the model parameters reaching a value of the objective function at the final estimate Sjm^  ¼ 0:15 and a wall-clock time 68.8 min, confirming the superior performance of the hybrid approach.

186

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

(a)

(b)

(c)

(d)

Fig. 12. Estimation results for VE1_10 with Rjj ¼ 1e4 : (a) contour plot, (b) history of the objective function, (c) convergence history of estimates of m1 , and (d) convergence history of estimates of m2 .

Table 2 Estimation results for linear FE model with six unknown parameters using modal properties. Case

Nk

# of modes (r)

VE2_01

8

VE2_02

12

3 6 3 6

^ i =mtrue Parameter estimate (m ) i i=1

i=2

i=3

i=4

i=5

i=6

0.95 0.99 0.93 1.01

0.96 0.98 0.93 1.01

1.24 1.05 1.40 0.97

1.12 1.07 1.37 0.96

1.04 1.00 1.04 1.00

1.01 1.02 1.05 0.99

^ S at m

Wall-clock time (min)

0.0006 0.0031 0.0006 0.0028

21.0 24.4 25.0 31.5

Table 3 ^ 0 ) and updated (with m ^  ) models with respect to the true model (with mtrue ). Comparison between modal parameters of initial (with m Mode

r

Natural frequency (xFE =xtrue ) 1

2

3

4

5

6

1

2

3

4

5

6

^ 0) Initial (m ^ ) VE2_01 (m

– 3 6 3 6

0.971 1.000 1.000 1.000 1.000

1.059 1.000 0.999 1.000 0.999

0.987 1.000 1.000 1.000 1.000

1.023 – 1.000 – 1.000

1.011 – 1.000 – 1.000

1.079 – 1.000 – 1.000

0.999 1.000 1.000 1.000 1.000

0.928 1.000 1.000 1.000 1.000

0.940 1.000 1.000 1.000 1.000

0.999 – 1.000 – 1.000

0.971 – 1.000 – 1.000

0.965 – 1.000 – 1.000

^ ) VE2_02 (m

Mode shapes (MACð/FE ; /true Þ)

4.1.2. Verification example 3: linear FEMU using time history response data In order to investigate the influence of the response quantities of the objective function on the identification results, acceleration time history response data are used, instead of modal properties, in the verification example in this section. The following definition is used to construct the objection function

SðmÞ ¼

r X ^est RRMSEðyi  y i ðmÞÞ

ð15Þ

i¼1

^est where RRMSEðyi  y i Þ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi PN 2 1 1 ^ est ðyi;k Þ2 is the k¼1 ðy i;k  y i;k Þ N N

relative root-mean-square error between the ith recorded time history response yi (obtained using mtrue and polluted by noise) and ^ est the ith FE predicted time history response y i , r denotes the number of measured responses, and N is the number of data samples in the time histories.

The same structure, six model parameters, and variables ^ 0 , T 0 , T N , Pm (mtrue , m 0 , Q , N c , Rjj , and k) employed in VE2 are used here, therefore, the only difference is the objective function used in the optimization problem (time history responses instead of modal parameters). Before running the seismic base excitation, the gravity loads are applied quasi-statically. Rayleigh damping is employed to model the energy dissipation mechanisms by assigning a damping ratio of 2% to the first and fourth modes of the structure (f 1 ¼ 1:05 Hz and f 2 ¼ 3:13 Hz). Ground acceleration recorded at the Sylmar County Hospital station during the 1994 Northridge earthquake (Fig. 14) is used as input base excitation. The 360° and 90° components are applied in the longitudinal and transverse directions of the building, respectively. The earthquake ground motions were recorded at a sampling rate of 50 Hz and have a N ¼ 600 data samples. The recorded peak ground acceleration were 0.84g and 0.60g for components 360° and 90°, respectively. Absolute acceleration response time histories in both horizontal

187

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

(a)

(b)

Acceleration (g) Acceleration (g)

Fig. 13. Estimation results for VE2_01 with n ¼ 6: (a) objective function and (b) convergence history of parameter estimates.

1 Sylmar 360° (Long.)

0 −1 1 Sylmar 90° (Trans.)

0 −1 0

2

4

6

8

10

12

Time (sec) Fig. 14. Ground acceleration records used as seismic input motions.

directions at the 3rd, 5th, and roof levels (see blue arrows in Fig. 9), i.e. r ¼ 6, are polluted by additive white Gaussian noises with a root-mean-square of 2.0%g and used as measured responses (y). In the estimation, a diagonal R matrix assuming a standard deviation of 1.0%g for each response is assumed. The diagonal entries of R are purposely taken different to the RMS values used to pollute the true response to mimic a real world case, where only an estimate of the covariance matrix of the observation error is possible. Table 4 summarizes the estimation results for VE3 and Fig. 15 plots the convergence history of the objective function and parameter estimates for VE3_02. Precise estimation of the model parameters is obtained and the FE model is accurately updated, reducing the objective function from Sjm^ 0 ¼ 486 to Sjm^  6 11:0. All the

parameters converge to their true values except m3 that converges when N k ¼ 8. When N k increases from 8 to 16, the estito 1:03mtrue 3 mation of the model parameters improves slightly further and all the model parameters converge to their true values while the objective function decreases to 7.3. Since the only difference between VE2 and VE3 is the objective function used in the optimization problem, it is concluded that better estimation results are obtained when using time histories instead of modal parameters, but at a higher computational cost. By comparing the wallclock time of VE3_01 (see Table 4) and VE2_01 (see Table 2), both with N k ¼ 8, it is noted that the computational cost increases by two to three times if acceleration time history responses are used in the objective function instead of modal parameters.

Table 4 Estimation results for linear FE model with six unknown parameters using acceleration time histories. Case

VE3_01 VE3_02

Nk

8 16

^ i =mtrue Parameter estimate (m ) i i=1

i=2

i=3

i=4

i=5

i=6

1.00 1.00

1.00 1.00

1.03 1.00

1.00 1.00

1.00 1.00

1.00 1.00

^ S at m

Wall-clock time (min)

11.0 7.3

66.0 127.2

188

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

(a)

(b)

Fig. 15. Estimation results for VE3_02: (a) objective function and (b) convergence history of parameter estimates.

4.2. Verification example 4: nonlinear FEMU using time history response data A 3D four-story two-by-one bay steel frame building (Fig. 16a) subjected to bi-directional horizontal seismic excitation is considered as a verification example of updating nonlinear FE models. Nonlinear seismic response of the structure is simulated using a mechanics-based nonlinear FE model with force-based fibersection beam-column elements (distributed plasticity). In these nonlinear models, FEs with fiber-discretized cross sections are used to simulate the spread of plastic areas along the element length. The nonlinear stress-strain behavior of a material is described by a uniaxial material constitutive law of the corresponding fiber, which are defined by time-invariant parameters

and can capture stiffness and strength degradation (see [46] for more details). The structure is designed for the same location as the 4-story building used for linear FEMU (see Fig. 9). The building has a story height of 3.5 m and bay width of 7.0 m and 8.0 m in the longitudinal (X) and transverse (Z) directions, respectively (Fig. 16a). The columns located on axis X = 0 m and X = 14 m have W14  61 cross section, and those on axis X = 7 m have cross section W14  90. Longitudinal beams on second and third levels have W21  62 cross section, while on fourth and roof levels have W21  55 cross section. Transverse beams on second and third levels have W18  40 cross section, while on fourth and roof levels have W18  35 cross section. Columns and beams are made of A992 and A36 steel, respectively.

(a) Y

(b)

nbf = 12 ntf = 2

Monitored fiber

(c) X 3.5 m (typ)

Stress (σ)

W14×90

ndw = 12

bE

f E

Strain (ε )

Z Fig. 16. FE model of the steel frame considered in the verification example of nonlinear models.

189

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

Fig. 17. Steel frame building used for verification examples of nonlinear FE models.

Fig. 17 shows the FE model of the nonlinear steel frame. Modal masses, gravity loads acting on beams, and geometry are shown. Fiber section force-based structural finite elements [46] are used to model the beams and columns and Gauss–Lobatto quadrature is used for numerical integration along the elements. Material nonlinearity can spread over several cross sections monitored along the element (integration points or IPs). The IPs are discretized into fibers, whose stress-strain behavior is characterized by nonlinear uniaxial material constitutive laws. A set of material model parameters (e.g., elastic modulus, yield stress, etc.) define material constitutive models. Cross-sections of beams and columns are discretized into longitudinal fibers as shown in Fig. 16b. The nonlinear uniaxial stress-strain behavior of the steel fibers is modeled using the modified Giuffré-Menegotto-Pinto model [17]. This material model is governed by eight parameters, three of which are primary parameters, while the other five are secondary parameters controlling the curvature of the hysteresis branches between consecutive strain reversal points. The primary material parameters consist of the elastic modulus (Es ), initial yield strength (f y ), and strain hardening ratio (b). Fig. 16c shows the uniaxial material model used for the reinforcing steel fibers together with its corresponding parameters assumed to be unknown in the estimation phase. The true material parameter values used to simulate the true true f ycol ¼ 345 MPa, bcol ¼ 0:04, true true ¼ 200 GPa, f ycol ¼ 250 MPa, and bcol ¼ 0:03, then true true true true true true ¼ ½Etrue m scol ; f ycol ; bcol ; Esbeam ; f ybeam ; bbeam . Mass and tangent

response

are

Etrue scol ¼ 200 GPa,

Etrue sbeam

stiffness-proportional Rayleigh damping, with a critical damping ratio of 2% for the first and second longitudinal modes of the building (f1 = 0.68 Hz and f2 = 1.32 Hz, after application of the gravity loads), is used to model the sources of energy dissipation distinct to material hysteretic energy dissipation. Components 360° and 90° recorded at the Sylmar Station during the 1994 Northridge earthquake (Fig. 14) are applied in the longitudinal (X) and transverse (Z) directions of the building, respectively. Different response measurements are considered in the estimation phase. First, true absolute acceleration response time histories in both horizontal directions at 2nd, 4th, and roof levels (see blue arrows in Fig. 16a) are polluted by additive white Gaussian noises with a root-mean-square of 0.5%g and used as measured responses,

y (i.e., ykþ1 2 R61 ). Two cases for matrix R are considered here, consisting of diagonal entries equal to (0.5%g)2 and (0.3%g)2. These cases are referred to as VE4_01 and VE4_03, respectively. In order to investigate the effect of heterogeneous response measurements, the estimation problem is also solved considering the relative displacement responses at the fourth and roof levels (see green circles in Fig. 16a) in addition to the contaminated absolute acceleration time histories used in previous case, i.e., ykþ1 2 R101 . The displacement responses are polluted by additive zero-mean white Gaussian noises with a RMS of 1.0 mm. In the estimation problem two cases for R are considered. In the first case, the diagonal entries of R corresponding to acceleration and displacement responses are assumed equal to (0.5%g)2 and (1.0 mm)2, respectively, while in the second they are taken as (0.3%g)2 and (0.7 mm)2, respectively. These cases are referred to as VE4_02 and VE4_04, respectively. In all cases, the following values for the other variables involved ^ 0 ¼ 0:7mtrue , T 0 ¼ 4, in the estimation problem are considered: m  m m 2  m m true true U L , T N ¼ 0:01, L ¼ 0:5 m , U ¼ 2:0 m , Pm 0 ¼ diag 100 k ¼ 0:0001, Q ¼ 0:01Pm 0 , N c ¼ 800, and N k ¼ 8. It is noted that N k ¼ 8 is chosen because, based on the results obtained in Sections 3 and 4.1, accurate parameter estimates are achieved with this number of iterations of the UKF. Using higher values of N k increases the computational cost of the estimation process with no significant improvement in the estimation results. Table 5 shows the estimation results for VE4. The objective function evaluated at the initial mode parameter estimates is Sjm^ 0 ¼ 568 and Sjm^ 0 ¼ 1001 for cases of acceleration only (y 2 R61 ) and acceleration + displacement (y 2 R101 ) response measurement, respectively. The objective function evaluated at the final model ^  ) is lower than or equal to 32.46 for parameter estimates (S at m all cases, thus confirming the successful updating of the FE model by minimizing the objective function. Good estimation results are true obtained for the initial stiffness (Etrue scol and Esbeam ) and yield-related true

true

parameters (f ycol and f ybeam ), with relative errors less than or true

true

equal to 6%. Post-yield related parameters (bcol and bbeam ) are only accurately estimated in case VE4_02, in which acceleration and displacement responses are considered and the true level of mea-

190

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

Table 5 Estimation results for the case of nonlinear model with six model parameters using acceleration and heterogeneous response time histories. pffiffiffiffiffiffi Rjj

Case

VE4_01 VE4_02 VE4_03 VE4_04

^ i =mtrue Parameter estimate (m ) i

0.5%g 0.5%g 1.0 mm 0.3%g 0.3%g 0.7 mm

i=1

i=2

i=3

i=4

i=5

i=6

1.01 0.99

0.59 1.02

0.99 0.99

0.99 1.00

1.10 1.00

28.26 22.82

585.1 601.2

1.00 0.99

1.05 1.06

1.09 0.50

1.00 1.00

0.95 0.95

1.78 1.88

24.94 32.46

593.9 598.6

ini

y4L (m/s2)

10 0

y2T (m/s2)

−10 20 0 −20 10

y5T (m/s2)

y2L (m/s2) y5L (m/s2)

Wall-clock time (min)

1.00 1.00

true

y4T (m/s2)

^ S at m

0 −10 2

3

4

5

6

Time (s)

fin

20 0 −20 10 0 −10 10 0 −10

2

3

4

5

6

Time (s)

Fig. 18. Comparison of true responses and estimated responses based on the final estimate of model parameters for case VE4_03 (time window 2.0–6.0 s).

surement noise is considered in R. The measured responses are not sensitive enough to post-yield parameters, and therefore limited true

true

information is contained in y to estimate bcol and bbeam . The wall-clock time is always higher than 580 min (9 h), implying a significant computational cost when updating detailed nonlinear mechanics-based FE models. However, it is important to emphasize that the use of structural nonlinear FE models requires estimating a limited number of model parameters because a large structural system can be characterized by a relatively small number of structural materials. Fig. 18 compares the true and estimated acceleration responses for case VE4_03. The excellent agreement between the true measured responses and the estimated responses based on the final estimate of the model parameters is clearly observed, with RRMSEs lower than 5.1%. It is noted that for all applications examples in Section 4, standard SA (N k ¼ 0) with a large number of annealing cycles (N c ) were also run, and worse estimation results at a considerably higher computational cost than those obtained with the hybrid approach were obtained.

5. Conclusions This paper presents a new method for finite element model updating (FEMU) of civil structures. The inverse problem of solving the unknown model parameters of structural FE models is solved using a global optimization approach, which combines the global search capabilities of simulated annealing (SA) with local improve-

ments of the accepted candidates through the unscented Kalman filter (UKF). The performance of the proposed method is investigated and the effects of different variables involved in the method are assessed. The proposed method improves the accuracy, convergence rate, and computational cost of the SA algorithm. Application examples consisting of a test function (the Styblinski–Tang function) with six parameters and linear and nonlinear FE models are presented to verify the performance of the proposed approach. For FEMU examples, numerically simulated data of realistic three-dimensional (3D) steel frame structures are employed. Different verification examples with objective functions defined in the modal (linear models) and time (linear and nonlinear models) domains involving different number of model parameters to be identified (from two to six) are presented. From the updating of the linear FE model with six unknown parameters, it is concluded that considering accelerations time histories to build the objective function provides better estimation results than using modal responses of the building. In addition, the advantages of using heterogeneous response measurements is shown through the updating of a mechanics-based nonlinear FE model of a 3D steel frame subjected to bidirectional horizontal seismic excitation. The updated FE models could be potentially used for damage identification purposes. Acknowledgements The first author acknowledges the support provided by the Universidad de los Andes through the research grant Fondo de Ayuda a la Investigacion (FAI) Iniciación 2016.

R. Astroza et al. / Computers and Structures 177 (2016) 176–191

References [1] Al-Hussein A, Haldar A. Novel unscented Kalman filter for health assessment of structural systems with unknown input. ASCE J Eng Mech 2015;141(7). [2] Allemang RJ, Brown DL. A correlation coefficient for modal vector analysis. In: Proc, 1st international modal analysis conference (IMAC I), Orlando, FL. p. 110–6. [3] Astroza R, Ebrahimian H, Conte JP. Material parameter identification in distributed plasticity FE models of frame-type structures using nonlinear stochastic filtering. ASCE J Eng Mech 2015;141(5):04014149. [4] Atalla MJ, Inman DJ. On model updating using neural networks. Mech Syst Signal Process 1998;12(1):135–61. [5] Azam SE, Ghisi A, Mariani S. Parallelized sigma-point Kalman filtering for structural dynamics. Comput Struct 2012;92:193–205. [6] Bakir PG, Reynders E, De Roeck G. Sensitivity-based finite element model updating using constrained optimization with a trust region algorithm. J Sound Vib 2007;305(1–2):211–25. [7] Bakir PG, Reynders E, De Roeck G. An improved finite element model updating method by the global optimization technique ‘Coupled Local Minimizers’. Comput Struct 2008;86(11–12):1339–52. [8] Balling RJ. Optimal steel frame design by simulated annealing. ASCE J Struct Eng 1991;117(6):1780–95. [9] Betti M, Facchinim L, Biagini P. Damage detection on a three-storey steel frame using artificial neural networks and genetic algorithms. Meccanica 2015;50:875–86. [10] Caicedo JM, Yun G. A novel evolutionary algorithm for identifying multiple alternative solutions in model updating. Struct Health Monit 2011;10 (5):491–501. [11] Cha YJ, Buyukozturk O. Structural damage detection using modal strain energy and hybrid multiobjective optimization. Comput-Aided Civ Infrastruct Eng 2015;30(5):347–58. [12] Chisari C, Bedon C, Amadio C. Dynamic and static identification of baseisolated bridges using Genetic Algorithms. Eng Struct 2015;102:80–92. [13] Distefano N, Rath A. System identification in nonlinear structural seismic dynamics. Comput Methods Appl Mech Eng 1975;5:353–72. [14] Ebrahimian H, Astroza R, Conte JP. Extended Kalman filter for material parameter estimation in nonlinear structural finite element models using direct differentiation method. Earthquake Eng Struct Dynam 2015;44 (10):1495–522. [15] Fan W, Qiao P. Vibration-based damage identification methods: a review and comparative study. Struct Health Monit 2011;10(1):83–111. [16] Fang X, Luo H, Tang J. Structural damage detection using neural network with learning rate improvement. Comput Struct 2005;83(25–26):2150–61. [17] Filippou FC, Popov EP, Bertero VV. Effects of bond deterioration on hysteretic behavior of reinforced concrete joints, UCB/EERC-83/19. Berkeley, CA: EERC report 83-19, Earthquake Engineering Research Center; 1983. [18] Friswell MI, Mottershead JE. Finite element model updating in structural dynamics. Dordrecht: Kluwer Academic Publishers; 1995. [19] Fritzen CP, Jennewein D, Kiefer T. Damage detection based on model updating methods. Mech Syst Signal Process 1998;12(1):163–86. [20] Hao H, Xia Y. Vibration-based damage detection of structures by genetic algorithm. J Comput Civ Eng ASCE 2002;16(3):222–9. [21] Hasançebi O, Dumlupinar T. Linear and nonlinear model updating of reinforced concrete T-beam bridges using artificial neural networks. Comput Struct 2013;119:1–11. [22] He RS, Hwang SF. Damage detection by an adaptive real-parameter simulated annealing genetic algorithm. Comput Struct 2006;84(31):2231–43. [23] International Code Council (ICC). International building code, Falls Church, VA, 2012. [24] Ingber L. Simulated annealing: practice versus theory. Math Comput Model 1993;18(11):29–57. [25] Jafarkhani R, Masri SF. Finite element model updating using evolutionary strategy for damage detection. Comput-Aided Civ Infrastruct Eng 2011;26:207–24. [26] Jeong I, Lee J. Adaptive simulated annealing genetic algorithm for system identification. Eng Appl Artif Intell 1996;9(5):523–32. [27] Julier SJ, Uhlmann JK. New extension of the Kalman filter to nonlinear systems. In: Proc SPIE 3068, signal processing, sensor fusion, and target recognition VI. p. 182–93.

191

[28] Kang F, Li J, Xu Q. Structural inverse analysis by hybrid simplex artificial bee colony algorithms. Comput Struct 2009;87(13):861–70. [29] Kirkpatrick S. Optimization by simulated annealing: quantitative studies. J Stat Phys 1984;34(5):975–86. [30] Levin RI, Lieven NAJ. Dynamic finite element model updating using simulated annealing and genetic algorithms. Mech Syst Signal Process 1998;12 (1):91–120. [31] Mariani S, Ghisi A. Unscented Kalman filtering for nonlinear structural dynamics. Nonlinear Dyn 2007;49(1–2):131–50. [32] Martin OC, Otto SW. Combining simulated annealing with local search heuristics. Ann Oper Res 1996;63(1):57–75. [33] Marwala T. Finite element model updating using computational intelligence techniques. Norwell: Kluwer Academic Publishers Group; 2010. [34] McKenna F, Fenves GL, Scott M. Open system for earthquake engineering simulation. Berkeley, CA: University of California; 2000. Available from: . [35] Meruane V, Heylen W. Damage detection with parallel genetic algorithms and operational modes. Struct Health Monit 2010;9(6):481–96. [36] Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. J Chem Phys 1953;21 (6):1087–92. [37] Moaveni B, He X, Conte JP, Restrepo JI. Damage identification study of a sevenstory full-scale building slice tested on the UCSD-NEES shake table. Struct Saf 2010;32(5):347–56. [38] Nguyen LT, Nestorovic´ T. Nonlinear Kalman filters for model calibration of soil parameters for geomechanical modeling in mechanized tunneling. ASCE J Comput Civ Eng 2015:04015025. [39] Nguyen LT, Nestorovic´ T. Unscented hybrid simulated annealing for fast inversion of the tunnel seismic waves. Comput Methods Appl Mech Eng 2016;301(1):281–99. [40] Perera R, Fang SE, Ruiz A. Application of particle swarm optimization and genetic algorithms to multi-objective damage identification inverse problems with modelling errors. Meccanica 2010;45(5):723–34. [41] Shabbir F, Omenzetter P. Particle swarm optimization with sequential niche technique for dynamic finite element model updating. Comput-Aided Civ Infrastruct Eng 2015;30:359–75. [42] Shiradhonkar SR, Shrikhande M. Seismic damage detection in a building frame via finite element model updating. Comput Struct 2011;89:2425–38. [43] Sirca GF, Adeli H. System identification in structural engineering. Sci Iran 2012;19(6):1355–64. [44] Song W, Dyke SJ. Real-time dynamic model updating of a hysteretic structural system. J Struct Eng ASCE 2014;140(3):04013082. [45] Tarantola A. Inverse problem theory and methods for model parameter estimation. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM); 2005. [46] Taucer FF, Spacone E, Filippou FC. A fiber beam-column element for seismic response analysis of reinforced concrete structures Report 91/17, EERC. Berkeley: Earthquake Engineering Research Center (EERC), University of California; 1991. [47] Teughels A, De Roeck G. Damage detection and parameter identification by finite element model updating. Arch Comput Methods Eng 2005;12 (2):123–64. [48] Teughels A, De Roeck G. Structural damage identification of the highway bridge Z24 by finite element model updating. J Sound Vib 2004;278 (3):589–610. [49] Teughels A, De Roeck G, Suykens JAK. Global optimization by coupled local minimizers and its application to FE model updating. Comput Struct 2003;81 (24–25):2337–51. [50] Van Der Merwe R, Wan EA. Efficient derivative-free Kalman filters for online learning. In: 9th European symposium on artificial neural networks (ESANN’2001), Bruges, Belgium. p. 205–10. [51] Yang J, Xia Y, Loh CH. Damage detection of hysteretic structures with pinching effect. J Eng Mech ASCE 2014;140(3):462–72. [52] Zapico JL, Gonzalez-Buelga A, Gonzalez MP, Alonso RA. Finite element model updating of a small steel frame using neural networks. Smart Mater Struct 2008;17(4):1–11. [53] Zarate B, Caicedo JM. Finite element model updating: multiple alternatives. Eng Struct 2008;30(12):3724–30. [54] Zimmerman AT, Lynch JP. A parallel simulated annealing architecture for model updating in wireless sensor networks. IEEE Sens J 2009;9(11):1503–10.