Updating Finite Element Model Using Stochastic Subspace Identification Method and Bees Optimization Algorithm

Abstract

This study investigates the application of operational modal analysis along with bees optimization algorithm for updating the finite element model of structures. Bees algorithm applies instinctive behavior of honeybees as they look for nectar of flowers. The parameters that needed to be updated are uncertain parameters such as geometry and material properties of the structure. To determine these uncertain parameters, local and global sensitivity analyses have been performed. An objective function is defined based on the sum of the squared errors between the natural frequencies obtained from operational modal analysis and finite element method. The natural frequencies of physical structure are determined by stochastic subspace identification method which is considered as a strong and efficient method in operational modal analysis. To verify the accuracy of this method, the proposed algorithm is implemented on a three-story structure to update parameters of its finite element model. Moreover, to study the efficiency of bees algorithm, its results are compared with those of the particle swarm optimization, and Nelder and Mead methods. The comparison indicates that this algorithm leads more accurate results with faster convergence.

Keywords
Vibration analysis; Finite element model updating; Operational modal analysis; Stochastic subspace identification; Bees algorithm; Sensitivity analysis

1 INTRODUCTION

Due to increase in demand for improving efficiency and reducing the weight of structures in modern industries such as aviation and aerospace industries, developing an accurate understanding of dynamic and vibrational behavior of systems and providing a precise model to describe their behavior are necessary. Accurate dynamic model ensures accuracy of subsequent analyses done on the structure such as structural health monitoring, damage detection and many others ( He and Fu 2001 ). Over the last few years, numerous methods have been introduced to accurately model the dynamics of a real mechanical system, one of which is the vibrational experiment and measured data analysis method. This method also called as modal analysis is considered as one of the most potent methods for obtaining dynamic parameters of structures such as natural frequencies, damping ratios and modal shapes ( Brincker et al 2000 ). When modal analysis is implemented on large structures such as bridges, dams and buildings, measuring exciting forces is almost impossible and output responses are the sole practical information that can be used for system identification algorithms. Vibrational methods based only on output responses are called operational modal analysis ( Pioldi et al 2016 ). Operational modal analysis methods are categorized in frequency and time domains. Peak picking, transmissibility and frequency domain decomposition are the most demanding methods among the frequency domain category. On the other hand, stochastic subspace identification, linear autoregressive method and Ibrahim's time domain method are considered as the important methods in time domain. Increasing demands of various applications for precise calculation of modal parameters have resulted in expanding operational modal in the past two decades. Zaghbani and Songmene (2009) used operational modal analysis as a powerful tool for estimating modal parameters of a machine-tool during machining operations. They applied the autoregressive moving average method and the least square complex exponential method for extracting modal parameters of a machine-tool with high speed. Then, these modal parameters were applied to determine machine dynamic stability lobes. Ebrahimi et al. (2013) calculated the natural frequencies of a cutting blade platform of a harvest combine machine using real data by frequency domain decomposition (FDD) method. They modeled the blade using finite element method by Ansys software and calculate the dynamic behavior of the cutting blade platform in a harvest combine, and updated the model by natural frequencies obtained from FDD method. Then, by using mass change strategy and finite element model, a modification technique was introduced to reduce the vibration of cutting blade platform. Van Overschee and De Moor (1990) introduced the stochastic subspace identification (SSI) method which has been proposed as an alternative to classical methods. In this stochastic method, the state space model is calculated from system output measurements. The key step in the SSI method is to calculate projection of future outputs matrix on past outputs matrix. Noël and Kerchen (2013) developed a new method by combining SSI and frequency methods in order to calculate modal parameters of a nonlinear system. To verify their method, they calculated modal parameters of a 5 degree of freedom system with two nonlinear parts. Dӧhler and Mevel (2010) examined SSI method in large structures. They pointed out that, to obtain modal parameters of a large structure, a large number of sensors should be used to record vibration signals, which incurs higher costs. To solve this problem, they developed a sensor-integrating method which could extract modal parameters of a large structure using few sensors. Goursat et al. (2011) analyzed vibrational data of a spacecraft during a commercial mission. The analysis was performed using SSI algorithm, the direct data utilization and covariance functions calculation methods. It was found that natural frequencies of the structure show small changes during the mission; however, mode shapes were more stable before and after its operation. Kompalka et al. (2007) compared data obtained from the finite element model and empirical model in order to examine accuracy of SSI algorithm in detection of damage and its location in structures. They demonstrated that SSI algorithm can predict the presence and location of damage with an appropriate accuracy. Moreover, they confirmed their approach by performing experiments on damaged beam.

Updating the finite element model is considered as an inverse method which involves reducing the difference between finite element model and empirical data. Methods based on gradient are extensively used in finite element model updating. Collins et al. (1974) have used inverse eigensensitivity method for updating finite element model. Failing to find global optimum point of the system is one of the fundamental problems of these methods. In addition, existence of optimum points in system boundaries results in reduction in efficiency of these methods. To cope with these problems, evolutionary intelligent optimization methods such as genetic, bees and particle swarm optimization algorithms were introduced, which inspired from natural evolution. These methods are not relying on the initial guess and do not need complicated mathematical concepts. Consequently, considering these characteristics, they can be used in finite element model updating. Dunn et al. (2005) applied genetic algorithm for updating finite element model of an F/A-18 aircraft based on experimental data. Moradi et al. (2010) updated a piping system with finite element method using bees algorithm and classical modal analysis data. Then, they compared their results with those obtained by genetic and PSO algorithms. Malekzehtab and Golafshani (2013) applied genetic algorithm for finite element model updating and damage detection of a jacket offshore platform. Their objective function was defined based on the natural frequencies and mode shapes of offshore platform. Chouksey et al. (2013) utilized experimental data for updating a rotating shaft with two journal bearing supports. They modeled the supports with linear and rotational springs and dampers. Moradi and Jamshidi moghadam (2015) used bees algorithm to find location and depth of crack in post-buckled beam-type structures. Bussetta et al. (2017) used discrete-time voltra series in order to update nonlinear finite element model of buckeled beam. Altunisik and Bayraktar (2017) updated model of Birecik Highway Bridge using operational modal analysis and finite element analysis. They changed some prominent parameters such as boundary conditions, material properties and section properties manually in order to reduce the difference between finite element and real models. In this study, the modal parameters of a system are calculated empirically using SSI method. Then, the effective parameters in model updating are obtained by performing a sensitivity analysis. Next, the sum of the squared errors between the natural frequencies obtained by the OMA and FEM is defined as the objective function, and will be minimized by the bees algorithm. In order to verify the integrity and the effectiveness of the proposed updating algorithm, the finite element model of a three-story frame is updated.

2 THEORY

2.1 Vibration Theory

Dynamic model of a linear system can be defined by a set of linear second-order differential equations with constant coefficients such as ( Goursat et al. (2011) ):

$$M\ddot{y}\left(t\right)+D\dot{y}\left(t\right)+Ky\left(t\right)=f\left(t\right)$$ (1)

where M, D and K represent the mass, damping and stiffness matrices, respectively. Also, y(t) is the vector of displacement and f(t) denotes vector of structure input forces. Equation (1) can be rewritten as a system of first order differential equations using various methods. Since measurement data are time-discrete, Equation (1) can be discretized such as the following steady-space model ( Dӧhler and Mevel (2010) ):

$$\begin{array}{l}{x}_{k+1}=A{x}_k+{\nu }_k\\ y_k=C{x}_k+w_k\end{array}$$ (2)

where x and y denote state and output vectors, respectively, and vectors v and w are input and output perturbations (system input is considered as a noise). A is state matrix and C is output matrix which is illustrated as a link between output and state variables. Moreover, v possesses noise due to modeling inaccuracy and w is the measurement noise related to data gathering system. Processing noise and data gathering noise are non-measurable signals and cannot be calculated by mathematical equations or signal processing, however, they are assumed as zero mean white noise signals between which the below relation exists ( Goursat et al. (2011) ).

$$E\left[\left[\begin{array}{c}{v}_r\\ w_r\end{array}\right]\left[\begin{array}{cc}{v}_s& w_s\end{array}\right]\right]=\left[\begin{array}{cc}Q& S\\ S^T& R\end{array}\right]{\delta }_{rs}$$ (3)

In Equation (3) , E denotes expected value operator and δ_rs is the Kronecker delta. Unscaled natural frequencies and mode shapes of system are obtained by eigenvalue and eigenvector analysis of matrix A and matrix C as:

$$\left(A-{\lambda }_{i}I\right){\varphi }_i=0\begin{array}{cc}& {\phi }_i\end{array}=C{\varphi }_i$$ (4)

where λ is unscaled natural frequency and φ is defined as mode shape of the system. Moreover, by using Equations (5) - (6) , scaled natural frequencies and damping ratios are estimated by the following relations:

$$a=\left|\arctan \frac{\mathrm{Im}\left(\lambda \right)}{\mathrm{Re}\left(\lambda \right)}\right|,\begin{array}{cc}& b=\ln \left|\lambda \right|\end{array}$$ (5)

$$\omega =\frac{a}{2\pi \tau }\begin{array}{cc}& \zeta =\frac{100\left|b\right|}{\sqrt{a^2+b^2}}\end{array}$$ (6)

where τ is considered as sampling rate. Therefore, vibration equation is converted to a system identification problem by solving which the desired modal parameters can be obtained.

2.2 Stochastic Subspace Identification Method

For obtaining modal parameters in time domain strategies, there is no need to transfer data to frequency domain, thus, eliminate any leakage error. Stochastic subspace identification method is one of the strong methods for estimating modal parameters in time domain. Making covariance functions of output data is regarded as the first step in SSI method ( Equation (7) ).

$$R_i=E\left(Y_k,Y_{k-i}^T\right)$$ (7)

where i is time delay and T is transpose symbol. Also, y_k is the output signals in time k which is defined in Equation (8) .

$$Y_k=Y\left(k\Delta t\right)\begin{array}{cc}& k=\mathrm{1,2,...,}N\end{array}$$ (8)

here N is the total number of samples. Since, output data are discrete in time domain, the expected value operator of them is stated as Equation (9) .

$$R_i=\frac{1}{N}{\displaystyle \sum _{k=1}^N{Y}_k}{Y}_{k-i}^T$$ (9)

In second step, the obtained covariance functions are used to build the Henkel matrix.

$$H_{p+\mathrm{1,}q}=\left(\begin{array}{cccc}{R}_0& R_1& \mathrm{...}& R_{q-1}\\ R_1& R_2& \mathrm{...}& R_q\\ .& .& \mathrm{...}& .\\ R_p& R_{p+1}& \mathrm{...}& R_{p+\mathrm{1,}q}\end{array}\right)$$ (10)

Next, the Henkel matrix is decomposed into controllable and observable matrices using singular value decomposition factorization.

$$H=o\eta$$ (11)

In Equation (11) , o is the observable matrix and η denotes controllable matrix. Matrix C, that plays an important role in modal identification, equals to the first row block of observable matrix ( Equation (12) ).

$$o_p^{\uparrow }\left(C,A\right)\stackrel{def}{=}\left(\begin{array}{c}C\\ CA\\ .\\ C{A}^{P-1}\end{array}\right)$$ (12)

In addition, the state transfer matrix A is calculated by shift invariance property of the observable matrix o and least square solution which is presented in Equation (13) .

$$o_P^{\downarrow }\left(C,A\right)=o_p^{\uparrow }\left(C,A\right)A\begin{array}{ccc}& & o_p^{\downarrow }=\left(C,A\right)\stackrel{def}{=}\left(\begin{array}{c}CA\\ C{A}^2\\ ⋮\\ C{A}^P\end{array}\right)\end{array}\Rightarrow A=o_p^{\uparrow }{\left(C,A\right)}^{†}{o}_p^{\downarrow }\left(C,A\right)$$ (13)

where superscript † denotes pseudo inverse matrix. By identifying state matrix A and according to Equations (5) - (6) , the natural frequencies and damping ratios can be obtained. Furthermore, the mode shapes of the system can be computed using the C matrix.

In order to calculate correct natural frequencies, order selection is regarded as a prominent factor in SSI method. If user underestimates order of SSI method, the method will not be able to obtain all of the natural frequencies. In contrast, spurious natural frequencies can be calculated if order is overestimated. Therefore, to solve this problem stabilization diagram is applied to distinct real natural frequencies from spurious ones. In stabilization diagram, different selected orders versus the natural frequencies obtained at each order is plotted. Real natural frequencies show a consistent behavior at different system orders whilst superior natural frequencies do not have stable attitude. For identification and omitting spurious natural frequencies in stabilization diagram, three criteria must be implemented ( Cara et al. (2012) ).

$$\frac{\left|{\omega }_{pi}-{\omega }_{qi}\right|}{{\omega }_{pi}}\le {\epsilon }_f$$ (14)

$$\frac{\left|{\zeta }_{pi}-{\zeta }_{qi}\right|}{{\zeta }_{pi}}\le {\epsilon }_{\zeta }$$ (15)

$$1-MAC\left({\varphi }_{pi},{\varphi }_{qi}\right)\le {\epsilon }_{MAC}$$ (16)

where ɛ_f, ɛ_ζ and ɛ_Mac are three limiting criteria which indicate that i^th real mode obtained from order p differs slightly from the same mode from order q (in most cases, p and q are consecutive numbers i.e. q =p+1). As a result, by observing these conditions in SSI algorithm, spurious modes can be omitted. These limiting criteria are chosen by trial and error and in this research they are considered as 0.02, 0.05 and 0.05, respectively.

2.3 Selections of the Best Degrees of Freedom for Operational Modal Test

In order to transfer the excitation energy to all degrees of freedom a system, the excitation points should not be located close to the nodes of the mode shapes of the system. By using ODP parameter given in Equation (17) , the distance of the degrees of freedom to the modal nodes can be estimated.

$$ODP\left(j\right)={\displaystyle \prod _{r1=1}^m\Vert {\varphi }_{jr1}^2\Vert }$$ (17)

where j and m are the degree of freedom and number of mode shapes of the system, respectively. The points of which ODP values are zero or close to zero are not suitable to be used for stimulation of the system, since they are on the mode shape nodes or near them. By contrast, points with maximal ODP values are appropriate for excitation and mounting the accelerometers. While exciting a structure by a shaker, the possibility of interference between the shaker and structure establishes that should be minimized. Each shaker is composed of a system of mass, spring and damper and any interference between it and the structure causes changes in signal generated by the shaker. To reduce this effect, shaker should be installed in locations where average acceleration is of minimum value. Average acceleration is defined using ADDOFA as in Equation (18) .

$$ADDOFA\left(j\right)={\displaystyle \sum _{r1=1}^m\frac{{\varphi }_{jr1}^2}{{\omega }_{r1}^2}}$$ (18)

The possibility of interface excitation in points with high ADDOFA is higher than the other points. Therefore, the points with higher ODP/ADDOFA are considered more appropriate points for excitation ( Imamovic (1998) ).

2.4 Sensitivity Analysis

Various mathematical models have been implemented for estimation of complex phenomena in different areas such as engineering, economy and physics. Identification of effective parameters in these mathematical models is the major challenge which users often deal with. Sensitivity analysis is one of the most important methods capable of identifying the parameters that have the most impact on results. In other words, any small changes in sensitive parameters would result in significant changes in output.

Once at a time index (OAT) is considered as one of the most applicable criteria for determination of sensitive parameters. Equation (19) states OAT index for identification of sensitive parameters. According to Equation (19) , a dimensionless index is considered in order to remove effects of various units of parameters ( Hamby (1994) ).

$$OAT=\frac{\Delta Y}{\Delta X}\frac{X}{Y}$$ (19)

where X and Y are input and output parameters of the model, respectively. Moreover, factor X/Y is defined as a normalized coefficient to omit effects of the units. The OAT Index defined in Equation (19) measures the local sensitivity. Global sensitivity index (GSI) is defined by Equation (20) in which overall sensitivity can be estimated ( Hamby (1994) ).

$$GSI=\frac{Y_{\max }-Y_{\min }}{Y_{\max }}$$ (20)

In this equation, Y_min and Y_max are the minimum and maximum output of the model using upper and lower limit bounds of the input parameters. According to Equations (19) - (20) , parameters that play significant role in output can be identified. In this research, input parameters are physical properties of the structure, whereas the natural frequencies are considered as outputs. By using sensitivity analysis, physical parameters which are suitable for optimization algorithm can be determined.

2.5 Bees Algorithm

Bees optimization algorithm is classified in evolutionary algorithms. There is an organized social behavior among bees which can be used for solving complex optimization problems. There are scout bees in each swarm whose main task is to find food sources for their hives. As scout bees find new food sources, they return to their hives and evaluate various discovered gardens based on specific parameters. Then, by performing a toggle dance, they provide the information of the direction, distance and amount of nectar in these gardens for worker bees. Next, the worker bees fly to the detected locations. The number of worker bees sent to these locations is proportional to the available nectar amount in the detected gardens. In other words, more worker bees are sent to gardens which have more nectars and shorter distance to the hive. Therefore, this strategy enables bees swarm to obtain food sources in an efficient manner. Figure 1 shows the flowchart of the bees algorithm. From N_t random solutions, N_t1 solutions which have the highest fitness values are selected as the best solutions. Then, among the best solutions N _t2 solutions are selected as the elite ones. In order to find better solutions, the best solution neighborhoods are searched. n_t1 and n_t2 denote the number of neighborhoods searched around the best and elite solutions, respectively (n_t1<n _t2).

Figure 1
Flowchart of bees optimization algorithm

Next, the remaining solutions are chosen randomly in the search space to find other solutions. Equation (21) indicates formation of a new generation in bees algorithm.

$$N_{new}=\underset{N_{t1}\begin{array}{cc}selected& \end{array}}{\underbrace{\left((N_{t1}-N_{t2})n_{t1}+(N_{t2})n_{t2}\right)}}+\left(N_{old}-N_{t1}\right)$$ (21)

where N_new and N_old are the new and the previous generation, respectively. Additionally, the number of population in each generation is fixed. These steps continue until the convergence criteria have been reached.

In this research, an objective function is defined to update the finite element model and bring it closer to the experimental model. It is expressed as the summation of the errors between the natural frequencies obtained by the operational modal analysis and the finite element method as represented by the following equation.

$$Errorfunction\left(\left\{Z\right\}\right)={\displaystyle \sum _{j=1}^s{\left({\Omega }_{c,j}^{}-{\Omega }_{uc,j}^{}\right)}^2}$$ (22)

where s is the number of natural frequencies used in optimization problem, Ω_c denotes natural frequencies calculated from the finite element method and Ω_uc are the natural frequencies obtained using the operational modal analysis. Moreover, vector Z includes design parameters which are identified by the sensitivity analysis (section (2-4)). By minimization of the objective function ( Equation (22) ) using bees optimization algorithm, design parameters are obtained and a precise finite element model based on real structure is designed.

3 RESULTS AND DISCUSSION

The main objective of this research is to optimize the finite element model of structures by using a combination of SSI, sensitivity analysis and finite element methods. The algorithm of the proposed method is described in Figure 2 .

Figure 2
Finite element model updating algorithm

To verify the proposed algorithm, a three-story structure is built numerically and experimentally, the updating algorithm is applied on it, and the results are presented in the following sections.

3.1 Finite Element Model

Figure 3 shows sketch of the three-story structure. It is built by steel bars with all connections welded. Bars are cut carefully and welded together using templates.

Figure 3
Plan of three views of the three-story structure

Moreover, anchor bolts are used to connect the structure to the ground. The dimensions of the three-story structure are presented in Table 1 .

Table 1
Dimensions of the three-story frame

No	Parameter	Symbol	Length (m)
1	Height of vertical members in the y-z plane	La1	0.600
2		La2	1.100 1.600
3		La3
4	Height of vertical members in the x-z plane	Lb1	0.600
5		Lb2	1.100 1.600
6		Lb3
7	Length of horizontal members in the x-y plane	Lc1	0.500
8		Lc2	0.400
9		a	0.016
10	Cross section	b	0.016

In Figure 4 , the finite element model of three-story structure in ANSYS software is displayed. Solid 186 element is used for building this model. Additionally, in this model, welds are taken into account as a change in Young's modulus in connections.

Figure 4
Finite element model of the three-story structure

Next, an eigenvalue analysis was performed on the model to obtain the natural frequencies and the results for the first six natural frequencies are tabulated in Table 2 . Also, the corresponding mode shapes are displayed in Figure 5 .

Table 2
Natural frequencies obtained by FEM

Frequency no.		Finite element method (rad/sec)
1	73.052
2	84.700
3	92.862
4	234.564
5	276.033
6	326.483

Figure 5
Modal shapes of the three-story structure

3.2 Appropriate Points for Shaker and Accelerometer Installation

Acquiring modal parameters requires an accurate planning for conducting experiments. Appropriate locations of accelerometers and shaker can lead to more precise modal parameters.

According to section (2-3), by using ODP and ODP/ADDOFA, the best points for accelerometers placement and structure stimulation by shaker can be detected. Figure 6 displays these points for installation of accelerometers and shaker in y direction. The points are also tabulated in Table 3 .

Figure 6
a- Best locations of accelerometers installation in y direction b- Best locations of shaker installation in y direction

Table 3
Coordinates of the best locations of accelerometers and shaker in y direction

Point	Length (cm)	point	Length (cm)
A1	95	A6	12
A2	110	A7	46
A3	122	A8	44
A4	144	A9	86
A5	160	A10	105

Figure 7 shows the best places for the installation of accelerometers and shaker in x direction and Table 4 presents their corresponding coordinates.

Figure 7
a-The Best locations for accelerometers installation in x direction b-The Best location for shaker installation in x direction

Table 4
Coordinates of the best locations of accelerometers and shaker in x direction

Point	Length (cm)	Point	Length (cm)
B1	91	B5	13
B2	120	B6	33
B3	148	B7	50
B4	160	B8	72

3.3 Sensitivity Analysis Results

In order to determine effective parameters for finite element model updating, a sensitivity analysis is performed on the three-story frame according to section (2-4). Table 5 presents the lower and upper limits of design parameters for sensitivity analysis (see section (2-4)).

Table 5
Upper and lower bounds of design parameters

no	Parameters	Symbol (unit)	Lower limit	Upper limit
1	Density	(kg/m3) ρ	7600	8000
2	Young's modulus	(Gpa) E	180	220
3	Poisson's ratio	([]) θ	0.2	0.4
4	Length of vertical members	Da1 (m)	1.55	1.65
5		Da2 (m)	1.55	1.65
6		Da3 (m)	1.55	1.65
7		Da4 (m)	1.55	1.65
8	Height of horizontal members in x-z plane	Dd1 (m)	0.58	0.62
9		Dd2 (m)	1.07	1.14
10		Dd3 (m)	1.55	1.65
11	Height of horizontal members in y-z plan	De1 (m)	0.58	0.62
12		De2 (m)	1.07	1.14
13		De3 (m)	1.55	1.65
14-29	Width of cross section	a1….a16 (m)	0.0155	0.0165
30-45	Length of cross section	b1….b16 (m)	0.0155	0.0165
46	Length of horizontal bars in x-z plane	Dc (m)	0.48	0.52
47	Length of horizontal bars in y-z plane	Dd (m)	0.38	0.42
48-71	Weld Young's modulus	E1…E24 (Gpa)	180	220

Figures 8 - 9 present the local and global sensitivity results for the three-story structure, respectively.

Figure 8
Results of local sensitivity analysis

Figure 9
Results of global sensitivity analysis

Figure 8 shows that the maximum local sensitivity belongs to design parameters related to length of the members, and physical characteristics of the structure such as Young's modulus and density. However, the minimum values are related to Young's modulus of welded connections and the cross-section of the members.

Moreover, in global sensitivity, approximately the same parameters that considered in local sensitivity can have direct impact on the obtained natural frequencies ( Figure 9 ).

3.4 Operational Modal Analysis

Operational modal analysis is regarded as a subset of modal analysis that only depends on output responses. In this research, SSI method is applied to identify dynamic parameters of the three-story frame. Random inputs are the main assumption of SSI method. Therefore, for random excitation of the three-story structure, an electro-dynamic shaker is used. This electro-dynamic shaker can produce random, burst random, pseudo random, sweep random and periodic random signals. Then, by using accelerometers attached on the structure, the output signals are captured and send to the time recorder software. Figure 10 shows equipment used in data acquisition of the structure.

Figure 10
Equipment used in empirical modeling

Figure 11 demonstrates the time response obtained from the three-story structure under random excitation. These data were recorded using a sampling rate of 16328 samples per second.

Figure 11
Response of the three-story structure under random excitation

Table 6 presents the natural frequencies obtained by SSI method under different excitations, plus those frequencies obtained by classical modal analysis.

Table 6
Experimental natural frequencies calculated by SSI and classical methods

Method	Input signal	Frequency no.
		1st	2nd	3rd	4th	5th	6th
SSI (method)	Random	62.230	71.830	90.120	205.230	239.230	293.970
	Burst random	63.140	73.430	91.250	208.380	240.460	294.080
	Pseudo random	62.450	71.210	89.640	205.090	237.900	292.010
	Sweep random	64.170	74.120	95.420	209.450	241.340	294.140
	Periodic random	59.230	78.230	96.150	210.130	244.250	296.340
Classical modal	Hammer	62.328	71.341	89.380	206.240	238.010	292.210

Figure 12 shows the relative error of natural frequencies of the structure obtained from operational modal analysis and classical modal analysis (hammer test). According to this figure, natural frequencies calculated by pusedo random signals are more consistent with those obtained from classic modal analysis, therefore, only this signal will be used in the following sections of this research.

Figure 12
The natural frequencies relative error between SSI and classical modal analysis

Stability diagram of SSI method for pusedo random excitation signal in y direction is depicted in Figure 13 .

Figure 13
Stability diagram obtained from SSI method for the three-story structure in y direction

As stated in section (2-2), points deployed next to each other in rows are an approximation of natural frequencies and can be seen in the figure clearly. As can be seen from Figure 13 , system order varies between 40 and 200. Table 7 compares the natural frequencies obtained by SSI method and those computed by the finite elements method. As can be seen from the table the relative error between the two methods is comparable.

Table 7
Natural frequencies (rad/sec) obtained from SSI and FEM methods

Frequency no.	SSI method	FEM method	Relative error (%)
1	62.450	73.052	14.513
2	71.210	84.700	15.927
3	89.640	92.862	3.470
4	205.090	234.564	12.566
5	237.900	276.033	13.815
6	292.010	326.483	10.559

3.5 Bees Algorithm Results

According to Table 7 , the relative error between the natural frequencies obtained by SSI and FE methods is high and therefore the finite element model needs to be updated. The mismatch may be due to the values of physical properties such as density, Young’ module or can be related to geometric properties of structure such as length, thickness and size of the members.

In order to update the finite element method and reduced the discrepancy between the computed and measured natural frequencies, an objective function is defined and minimized by the bees algorithm. Three different sets of design parameters are defined for bees algorithm.

The first set includes parameters obtained from sensitivity analysis, second set includes parameters not important for sensitivity analysis and material properties, and the third set consists of all design parameters. Table 8 summarizes all the design parameter sets in this study. As shown in the table, these sets consist of 14, 59 and 71 design variables, respectively.

Table 8
Design parameter sets used in FEM updating

Design Parameters	Set 1	Set 2	Set 3
ρ	√	√	√
E	√	√	√
Da1	√	×	√
Da2	√	×	√
Da3	√	×	√
Da4	√	×	√
Dd1	√	×	√
Dd2	√	×	√
Dd3	√	×	√
De1	√	×	√
De2	√	×	√
De3	√	×	√
a1….a16	×	√	√
b1….b16	×	√	√
Dc	√	×	√
Dd	√	×	√
E1…E24	×	√	√
ν	×	√	√

Figure 14 shows the convergence of bees algorithm for different sets of design parameters. As it is clear from this figure the convergence of the first and third sets are faster and the results are more precise than the second set. However, the convergence occurred with 14 design parameters by set 1 and 71 design parameters by set 3. Additionally, the run time of the optimization process in set 1 is lower than that in set 3. This shows that optimization using set 1 could achieve optimal response in minimum time. Moreover, since all of design parameters in set 1 were obtained by sensitivity analysis, it is clear that sensitivity analysis is able to identify effective parameters in optimization process.

Figure 14
Convergence of bees algorithm

Table 9 presents control parameters utilized in bees optimization algorithm. These values have been obtained by empirical studies.

Table 9
Control parameters of bees algorithm

Control Parameter	Values
Number of objective functions	40000
Neighborhood radius	0.01
Nt1	20
Nt2	12
nt1	10
nt2	15
Nt	100
Number of iterations	100

The natural frequencies obtained from updated finite element model are tabulated in Table 10 . It can be deduced from Tables 7 and 10 that relative error between natural frequencies calculated by finite element and experimental methods decreased remarkably after model updating. Total relative error before model updating was %70.850 that reduced to %0.021 after model updating

Table 10
Design parameter sets used in FEM updating

Mode number	Natural frequencies (rad/sec)		Relative error (%)
	Experimental model	Updated finite element model
1	62.450	62.451	0.002
2	71.210	71.213	0.004
3	89.640	89.647	0.008
4	205.090	205.098	0.004
5	237.900	237.904	0.002
6	292.010	292.008	0.001

3.6 Validation of Bees Algorithm

In order to the verify the performance of bees optimization algorithm in model updating, the results obtained by this algorithm are compared with those of particle swarm optimization (PSO) algorithm, one of the most popular evolutionary algorithms, and Nelder and Mead simplex method. PSO was introduced by Kennedy and Eberhart (1995) . This method is inspired by migration behavior of birds or fish schooling. PSO algorithm is based on the particles movement in which each particle can be considered as a possible optimal solution. Particles in PSO method follow a very simple behavior: trying for neighboring particles success and their success. Finding optimal solution area in the search space with a large number of dimensions is the result of mass movement behavior. The Nelder-Mead algorithm is a direct method for finding minimal value of an unconstrained ns-dimensional objective function presented by Nelder and Mead (1965) .

This method is based on the comparison of the function value in the n_s + 1 vertices of a simplex and replacing the worse vertex in terms of the objective function with a new point. Moreover, the method works only by using the function values and does not need any function derivatives.

Figure 15 portrays the convergence of BA, PSO and Nelder-Mead methods for finite element model updating of the three-story structure using first set of design parameters.

Figure 15
Evolution of the best function values

As can be seen from the Figure 15 , BA algorithm converges faster than other algorithms and its objective function has the lowest value. Relative errors between natural frequencies obtained from the finite element updated model and the experimental modal analysis using different optimization methods are tabulated in Table 11 . The table shows that the natural frequencies calculated by BA are more accurate than other optimization algorithms. The optimized values of the first set of design parameters are tabulated in Table 12 .

Table 11
Relative error (%) between natural frequencies of updated finite element model with experimental frequencies

Optimization Method	Frequency no
	1st	2nd	3rd	4th	5th	6th
Bees	0.002	0.004	0.008	0.004	0.002	0.001
PSO	0.034	0.038	0.056	0.095	0.053	0.081
Nelder-Mead	1.725	1.106	1.034	0.503	0.448	1.001

Table 12
Design parameter sets used in FEM updating

Design parameters	BA	Nelder-Mead	PSO
(kg/m3) ρ	7710	7650	7765
(Gpa) E	198	203	201
Da1	1.61	1.63	1.65
Da2	1.61	1.63	1.65
Da3	1.61	1.62	1.65
Da4	1.61	1.62	1.64
Dd1	0.59	0.61	0.60
Dd2	1.08	1.13	1.08
Dd3	1.64	1.63	1.63
De1	0.61	0.58	0.62
De2	1.09	1.11	1.09
De3	1.64	1.63	1.64
Dc	0.51	0.48	0.50
Dd	0.39	0.38	0.39

To evaluate the efficiency of optimization algorithms, two criteria are defined, practical reliability index and normalized price value. Practical reliability index is defined as the probability of the solution to reach a practical optimum. A practical optimum is stated as an optimal solution within 0.1% of the global optimum. In this study, the best solution found after 80000 objective functions evaluations were used to define global optimum. Moreover, the normalized price value is defined as the ratio of the number of evaluated objective functions to the practical reliability.

The values of practical reliability and normalized price for the three optimization algorithms are shown in Table 13 . According to this table, BA algorithm has the highest practical reliability amongst all the methods. It also has the lowest cost to find the global optimum.

Table 13
Practical reliabilities and normalized prices obtained from BA, PSO and Nelder-Mead methods

Method	practical reliability	normalized price
BA	0.9	44444.4
PSO	0.7	57142.9
Nelder-Mead	0.6	66666.7

3.7 Mode Shapes

Figure 16 portrays the mode shapes obtained from the updated finite element model and SSI method.

Figure 16
Mode shapes using FEM and SSI methods

To compare the mode shapes obtained from SSI and FEM methods, the mode assurance criterion (MAC) is used which is defined by the following equation:

$$MAC(i,j)=\frac{{\left|{\left\{\varphi \right\}}_{SSI}^T{\left\{\varphi \right\}}_{FEM}^T\right|}^2}{\left|{\left\{\varphi \right\}}_{SSI}^T{\left\{\varphi \right\}}_{SSI}^T\right|\times \left|{\left\{\varphi \right\}}_{FEM}^T{\left\{\varphi \right\}}_{FEM}^T\right|}$$ (23)

where ϕ is the mode shape. MAC values greater than 0.7 shows good correlation between the two methods. The MAC values for the mode shapes of SSI and FEM is shown in Figure 17 . As can be seen from the Figure 17 , the mode shapes calculated by these methods are compatible with each other and the value of MAC criteria is more than 0.9 for the same mode shapes.

Figure 17
MAC criteria between SSI and FEM mode shapes

4 CONCLUSIONS

In this study, a method to update the finite element model of structures has been developed using operational modal analysis and bees optimization algorithm. An objective function based on the summation of the squared errors between the natural frequencies obtained from the numerical method and experimental analysis was defined. Numerical and experimental frequencies obtained by the finite element method and operational modal analysis, respectively. By minimization of this function, the values of uncertain design parameters were calculated. To validate the proposed method, it was implemented on numerical and physical model of a three-story frame. Numerical frequencies were calculated by eigenvalue analysis of Ansys software and experimental frequencies were obtained by stochastic subspace identification. Results show that the combined errors between the first six natural frequencies were reduced from %70.850 to %0.021. In addition, a very good compatibility was found between the mode shapes obtained by the finite element and operational modal analysis methods where the MAC values were more than 0.9 for all the same modes. Moreover, to investigate the efficiency of the bees algorithm, its results compared with those obtained by particle swarm optimization and Nelder and Mead methods. The results indicated that the bees algorithm is faster and more efficient compare with those methods and can successfully bring the finite element model closer to the physical model.

• http://dx.doi.org/10.1590/1679-78254189

References

• Altunisik, A.C., Bayraktar, A., (2017). Manual model updating of highway bridges under operational condition. Jou-rnal of Structures and Systems, 19: 39-46.
• Brincker, R., Zhang, L., Andersen, P., (2000). Modal identification from ambient responses using frequency domain decomposition. Proceedings of the International Modal Analysis Conference (IMAC), San Antonio, TX, USA.
• Bussetta, P.,Shiki, S.B., da Silva, S., (2017). Updating of a nonlinear finite element model using discrete-time voltera series. Latin American Journal of Solid and Structures, 14:1183-1199.
• Cara, F.J, Juan, J., Alarco’n, E., Reynders, E., De Roeck, G., (2012). Modal contribution and state space order selection in operational modal. Journal of Mechanical Systems and Signal Processing, 38: 276–298.
• Chouksey, M., Dutt, J.K., Modak S.V., (2013). Model updating of rotors supported on ball bearings and its application in response prediction and balancing. Journal of Measurement, 46: 4261-4273.
• Collins, J.D, Hart, G.C., Hasselman, T.K., Kennedy, B., (1974). Statistical Identification of Structures. AIAA Journal, 12: 185-190.
• Dӧhler, M., Mevel, L., (2010). Modular subspace-based system identification from multi-setup measurements. IEEE Transactions on Automatic Control, 5: 2951-2956.
• Dunn, S., Peucker, S., Perry, J., (2005). Genetic algorithm optimisation of mathematical models using distributed computing. Journal of Applied Intelligence, 23: 21-32.
• Ebrahimi, R., Esfahanian, M. and Ziaei-rad, S., (2013). Vibration modeling and modification of cutting platform in a harvest combine by means of operational modal analysis (OMA). Journal of Measurement, 46: 3959-3967.
• Goursat, M., Dӧhler, M., Mevel, L., Anderson, P., (2011). Crystal clear SSI for operational modal analysis of aerospace vehicles. 28th International Modal Analysis Conference, (Florida).
• Hamby, D.M., (1994). A review of techniques for parameter sensitivity analysis of environmental models. Journal of Environmental Monitoring and Assessment, 32: 135-154.
• He, J., Fu, Z. F., (2001). Modal analysis, Oxford, (London).
• Imamovic, N., (1998). Validation of large structural dynamics models using modal test data Ph.D. Thesis, Imperial College, London, England.
• Kennedy, J., Eberhart, RC., (1995). Particle swarm optimization. In: Proceedings of the IEEE international Conference on Neural Networks, 4: 1942–1948.
• Kompalka, A.S, Reese S., Bruhns, O.T., (2007). Experimental investigation of damage evolution by data-driven stochastic subspace identification and iterative finite element model updating. Journal of Applied Mechanics, 77: 559-573.
• Malekzehtab, H., Golafshani, A.A., (2013). Damage detection in an offshore jacket platform using genetic algorithm based finite element model updating with noisy modal data. Procedia Engineering, 54: 480-490.
• Moradi, S., Fatahi, L., Razi, P., (2010). Journal of Finite element model updating using bees algorithm. Structural and Multidisciplinary Optimization. 42:283-291.
• Moradi, S, Jamshidi Moghadam, P., (2015). Crack identification in post-buckled beam-type structures. Journal of Smart Structures and Systems, 15:1233-1252.
• Nelder, J.A, Mead, R., (1965). A simplex method for function minimization. Journal of Compute,7: 308-313.
• Noël, J.P., Kerchen, G., (2013). Frequency-domain subspace identification for nonlinear mechanical systems. Journal of Mechanical System and Signal Processing, 40: 701–717.
• Pioldi, F., Ferrari, R, Rizzi, E., (2016). Output-only modal dynamic identification of frames by a refined FDD algori-thm at seismic input and high damping. Journal of Mechanical Systems and Signal Processing, 68-69: 265-291.
• Van Overschee, P., De Moor, B., (1990). Subspace identification for linear system. Journal of Mechanical System and Signal Processing, 5: 101-135.
• Zaghbani, I., Songmene, V., (2009). Estimation of machine-tool dynamic parameters during machining operation trough operational modal analysis. International Journal of Machine Tools and Manufacture, 49: 947-957.

Publication Dates

• Publication in this collection
  2018

History

• Received
  01 July 2017
• Reviewed
  12 Oct 2017
• Accepted
  18 Nov 2017