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

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.


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.

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)): 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 nonmeasurable 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)).
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: 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: 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.

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)).
where i is time delay and T is transpose symbol. Also, yk is the output signals in time k which is defined in 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).
In second step, the obtained covariance functions are used to build the Henkel matrix.
Next, the Henkel matrix is decomposed into controllable and observable matrices using singular value decomposition factorization.
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)).
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).
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)).
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.

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.
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).
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)).

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)).
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)).
In this equation, Ymin and Ymax 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.

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 Nt random solutions, Nt1 solutions which have the highest fitness values are selected as the best solutions. Then, among the best solutions Nt2 solutions are selected as the elite ones. In order to find better solutions, the best solution neighborhoods are searched. nt1 and nt2 denote the number of neighborhoods searched around the best and elite solutions, respectively (nt1<nt2).
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.
where Nnew and Nold 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.
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.

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. 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. 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. 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. 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.

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.

Figure 6: a-Best locations of accelerometers installation in y direction b-Best locations of shaker installation in y direction
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 7 shows the best places for the installation of accelerometers and shaker in x direction and Table 4 presents their corresponding coordinates.

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)). Figures 8-9 present the local and global sensitivity results for the three-story structure, respectively. 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.

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 electrodynamic 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 modelling 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.  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.

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. 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. 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.

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 ns + 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. 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. 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.

Mode Shapes
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.

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.