An ED-PSO model updating algorithm for structure health monitoring of beam-like structures

2.1. Eigenvector difference approach

If the stiffness matrix $K$ and mass matrix $M$ of the structure are known, the generalized eigenvalues and generalized eigenvectors of the two matrices can be obtained from Eq. (1) [31]:

where ${\lambda }_i$ and ${\mathrm{\Psi }}_i$ represents the $i$th order generalized eigenvalues and generalized eigenvectors, respectively. The objective function formula of the eigenvector difference approach is shown in Eq. (2) [32]:

where $n_{mode}$represents the order of system vibration, ${\lambda }_i^{Exp}$ and ${\psi }_i^{Exp}$ represents the generalized eigenvalue and generalized eigenvector of the actual model, respectively. The generalized eigenvector represents the deflections of the $i$th vibration mode with different degrees of freedom when the system is free to vibrate. If the case under study is the deformation of the system under external loads, then ${\psi }_i$ in Eq. (2) can also represent the deflection vector of the system deformation.

However, for continuous systems, solving the mass stiffness matrix is very troublesome, so it is necessary to improve the objective function of the eigenvector difference approach. Yang Wang [33] et al. proposed the relationship between the generalized eigenvalue of a structure and its natural frequency, as shown in Eq. (3):

where $f$ represents the natural frequency of the structure, thus replacing the generalized eigenvalue in the eigenvector difference approach with the natural frequency. There is no need to manually solve the mass stiffness matrix of the structure. By introducing Eq. (3) into Eq. (2), an improved objective function formula for the eigenvector difference approach can be obtained, as shown in Eq. (4):

where, $f_{ai}$ and $f_{bi}$ represent the first n-order natural frequencies of the model to be updated and the actual model, ${\mathrm{\Psi }}_{bj}$, ${\mathrm{\Psi }}_{aj}$ represents the deflection calculated by the finite element model and the actual model for k points, respectively.

2.2. ED-PSO model updating algorithm

ED-PSO model updating algorithm optimization principle is the same as particle swarm algorithm, first assume that there is a swarm of particles in the space, to describe the characteristics of the particles in terms of position, velocity and adaptation value, the particles in the solution space through the individual extreme value and swarm extreme value to continuously update their own position velocity information, to derive a new adaptation value, and compare with the swarm extreme value, and finally find the assumption that in a D-dimensional search space, there are $N$ swarm particles where the $i$-th particle can be represented as an i-dimensional vector: $X_i=\left(x_{i1},x_{i2},\dots x_{iD}\right)$, the velocity of the $i$-th particle can also be expressed as an i-dimensional vector: $V_i=\left(v_{i1},v_{i2},\dots v_{iD}\right)$, the optimal position searched by the $i$-th particle so far is called the individual extreme value: $P_{best}=\left(p_{i1},p_{i2},\dots p_{iD}\right)$, the optimal position searched so far by the whole particle swarm is the global extreme value: $G_{best}=\left(g_{i1},g_{i2},\dots g_{iD}\right)$, and after finding these two optimal values, the particle swarm is updated by Eqs. (5) and (6) as follows:

where $\omega$ denotes the inertia weight (generally taken as 0.8-1.2), $i=$ 1,2,…,$N$, $j=$1, 2…,$D$, $k$ refers to the number of current iterations calculated, $v_{ij}$ refers to the velocity of the current particle, $x_{ij}$ refers to the position of the current particle, $c_1$, $c_2$ are acceleration factors, and $r_1$, $r_2$ are distributed as random numbers in [0, 1]. In order to prevent the particles from searching blindly in the solution space, the numerical intervals $[v_{min},v_{max}]$, $[x_{min},x_{max}]$ are generally set for the velocity and position of the particles, and the schematic diagram of the algorithm operation is shown in Fig. 1.

Fig. 1The schematic diagram of Particle Swarm Optimization

However, the natural frequency and deflection often need to be calculated, which means that if the natural frequency and deflection are directly used as the variables of the objective function, the range of the function variables cannot be directly determined, and it is not possible to randomly select a certain number of particles and iterative optimization within a certain range. However, the natural frequency and deflection are usually related to the material parameters of the structure, and the material parameters of the structure can be determined in a general range, so the ED-PSO model updating algorithm usually takes the material parameters of the structure as the object of the particle iterative optimization, and the ED-PSO model updating algorithm is implemented in the following steps:

(1) Select the studied structural material parameters and determine their ranges;

(2) Set the number of particle swarms, maximum number of iterations, inertia weights, learning factors, maximum and minimum particle positions, maximum and minimum values of motion velocity, and objective function tolerances for the ED-PSO model updating algorithm;

(3) Generate random particle swarms within the parameter setting range;

(4) Import the generated parameters into ANSYS for relevant calculations, obtain the natural frequency and deflection under each group of particles, and bring them into the objective function for comparison to obtain the local optimal particle;

(5) Enter into the iterative optimization process, as in the particle swarm algorithm, the particles are moved within the parameter setting range, and the real particles are continuously updated;

(6) Import the iterative optimization example into ANSYS for calculation, and get the corresponding natural frequency and deflection, and bring it into the objective function for calculation, and get the global optimal particle;

(7) Output the values of natural frequency and deflection at this point as the result of the model updating calculation.

3.1. Model updating of a reinforced concrete beam

In order to verify the superiority of the proposed ED-PSO model updating algorithm, this section takes a simply supported beam as the research object and performs the finite element model updating calculation by two methods respectively, and judges the advantages and disadvantages of the two methods by comparing the errors of the derived deflections and natural frequencies with the actual model.

The dimensions of the simply supported beam are 3 m long, 0.15 m wide and 0.3 m high with four reinforcement bars, and the cross-sectional dimensions of the reinforced concrete beam are shown in Fig. 2.

Fig. 2Sectional dimension picture of reinforced concrete beam

The modeling software is ANSYS APDL2021R1, which has strong finite element calculation capability and is very accurate in solving structural modal parameters and load problems. The concrete material is simulated by solid65 element, and the elastic modulus of concrete is set as 3.25×10⁴ MPa and the density is 2500 kg/m³; the reinforcement is simulated by Link180 element, and the elastic modulus of reinforcement is 2×10⁵ MPa and the density is 7850 kg/m³. The grid is divided as follows: 40 sections in the length direction of the reinforced concrete beam, 10 sections in the height and width direction, the grid shape is hexahedral cells, a total of 4000 cells, and the grid is divided by mapping, so that the grid is neatly arranged, and the constraint adopts simple support constraint, considering that the structure is a three-dimensional structure, so one end is displacement constrained in $x$, $y$, and $z$ directions, and the other end is displacement constrained in x and y directions with no moment constraint. The initial finite element model of the reinforced concrete beam is shown in Fig. 3.

Fig. 3Initial finite element model of reinforced concrete beam

For the purpose of finite element model updating work, the reinforced concrete beam is cut into ten equal parts, each part being numbered 1-10. Apply a force of 100 KN in the middle of the reinforced concrete beam, as shown in Fig. 4.

Fig. 4Schematic diagram of reinforced concrete beam

The elastic modulus of concrete of parts No. 2, 4 and 7 is reduced by 20 % as the actual model of the reinforced concrete beam. The middle deflection of beam bottom and the first third-order natural frequencies of beam vibration under free state at the beam length ($z$ direction) of 0.3 m, 0.6 m, 0.9 m, 1.2 m and 1.5 m under the load of 100 KN are taken as the reference data for the updating accuracy of the finite element model. According to the sectional structure of reinforced concrete beam and the calculation formula of the modal parameters of simply supported beam, it can be seen that the physical parameters that have a great influence on the modal parameters of beam structure are mainly the elastic modulus of concrete, concrete density, elastic modulus of reinforcement, section area, etc. The actual structure has equal cross-section and almost no change in cross-section, however, due to factors such as uneven mixing of concrete and impurities in reinforcement, the elastic modulus of concrete, concrete density and elastic modulus of reinforcement vary greatly, so the above three variables are selected here as the variables to be updated. Since concrete is the main material of the simply supported beam structure, the elastic modulus of concrete and density are set as local variables, according to the numbering in Fig. 3.7, the elastic modulus of concrete is set as $E_1,E_2,\dots ,E_{10}$,concrete density is set as ${\rho }_1,{\rho }_2,\dots ,{\rho }_{10}$, and the elastic modulus of reinforcement is a global variable set to $E_k$. In this way, there are 21 variables to be updated, and all parameters take values between 0.8 times and 1.2 times the design value of the parameters.

The finite element model updating numerical calculation of reinforced concrete beam is realized by ED-PSO model updating algorithm, setting the swarm size as 50, the maximum number of iterations as 200, each particle corresponds to 21 updating variables, the inertia weight w is set as 0.729, the learning factor $c_1=c_2=$2, the particle position interval is the interval where the parameters take values,and the maximum and minimum value of velocity is ±0.4 times of the parameter design value.The errors of the updated frequencies and deflections with the actual model are obtained as a reference.

First, determine the value interval of the elastic modulus of concrete ($E_{concrete}$), density of concrete (${\rho }_{concrete}$) and elastic modulus of reinforcement ($E_k$), as shown in Table 1.

Since the first three orders of natural frequencies and five measured points of deflection are used as a reference for the updating results, Eq. (4) is rewritten as Eq. (7):

Because each swarm represents a 21-dimensional vector, in order to visualize the initial value of each variable in the vector, here the concrete elastic modulus, concrete density, and steel elastic modulus are combined, i.e., combined into three-dimensional coordinate points $\left(E_1,{\rho }_1,E_k\right)$, $\left(E_2,{\rho }_2,E_k\right)$,…, $\left(E_{10},{\rho }_{10},E_k\right)$, and the resulting random particle coordinate points are shown in Fig. 5.

By solving the particle swarm local optimum, the individual local optimum particle is obtained, at which time $E_k=$2.27×10¹¹ Pa, and similarly, the individual local optimum particle is expressed in two-dimensional coordinates $\left(E_1,{\rho }_1\right)$, $\left(E_2,{\rho }_2\right)$,…, $\left(E_{10},{\rho }_{10}\right)$, the individual local optimum particles are shown in Fig. 6, individual optimal particle values are shown in Table 2.

Fig. 5Initial particle distribution

Fig. 6Individual optimal particle distribution

Import the above data into ANSYS as a text file, solve the structure’s natural frequency and deflection, and compare with the actual structure. The graphs of the first three orders of modal analysis of the reinforced concrete beam under load are shown in Fig. 7.

Fig. 7Diagram of first third order vibration and load deformation of beam

a) First order vibration

b) Second order vibration

c) Third order vibration

d) Load deformation

Then set the parameters of the actual structure, calculate the corresponding deflection and natural frequency in the case of the actual structure, and finally calculate the deflection and natural frequency of the structure to be updated in the individual optimal case as the result before the model updating, and calculate the error of the deflection and natural frequency of the model to be updated with the actual structure by Eq. (8) and Eq. (9):

The final relevant calculation results are shown in Table 3.

The global optimal particle is then solved, at which point $E_k=$1.85×10¹¹ Pa, and the two-dimensional coordinates of the variables taken are shown in Fig. 8. The global optimal particle results are shown in Table 4.

Similarly, save the calculated global optimal particles to a text file, and then import the data into ANSYS for the solution of deflection and natural frequency, and then carry out the error calculation of the deflection under the actual structure, the results of which are shown in Table 5. The updated error is obtained and compared with the pre-updating error, as shown in Fig. 9.

It is easy to see from the obtained data that the ED-PSO model updating algorithm has a high updating accuracy for the finite element model updating method.

Fig. 8Global optimal particle distribution

Fig. 9Error comparison before and after updating

3.2. Comparison with CSD-PSO

In order to verify the advantages of the proposed method, simply supported beams in the previous section are used as the research object in this section, and performs the model updating calculation by the previous study - the combination of static and dynamic methods and particle swarm algorithm (CSD-PSO), and compares the results with the proposed method, and the results show that the proposed method has higher accuracy.

Particle swarm algorithm-based finite element model updating techniques combined with a variety of objective functions, model updating accuracy reference data are mainly natural frequency, deflection, etc.[34], static and dynamic methods are commonly used to measure the above data. The data measured by the static method are more accurate and easier to measure, but the measured data cannot reflect the dynamic changes of the structure. The data measured by dynamic method is less accurate and the testing process is tedious, but it can reflect the dynamic changes of the structure. In order to combine the advantages of the two, the combination of static and dynamic methods of finite element model updating was produced and widely used in the study of finite element model updating. However, while the combination of static and dynamic methods combines the advantages of both, it also amplifies the disadvantages of both, and the combination of static and dynamic methods requires more objective functions and is more complicated to calculate. The commonly used objective function of the combination of static and dynamic methods is shown in Eq. (10) [35]:

where $f_m$ and $f_s$ are the dynamic and static objective functions of the structure, respectively. The objective function of the combination of static and dynamic methods is shown in 3.4, where $f_m$ and $f_s$ are the dynamic and static objective functions of the reinforced concrete beam, respectively, as shown in Eqs. (11 and 12):

where ${\omega }_{bi}$ and ${\omega }_{ai}$ denote the first three orders of natural frequencies calculated by the actual structure and the first three orders of natural frequencies calculated by the model to be updated, respectively. $f_{bi}$, $f_{ai}$ denote the middle deflection of the bottom of the beam at 0.3 m, 0.6 m, 0.9 m, 1.2 m and 1.5 m in the direction of the length of the extended beam under the load action of the actual structure and the structure to be updated, respectively.

As in the previous subsection, the range of values of concrete modulus of elasticity, concrete density, and steel modulus of elasticity is still the same as the range of values in Table 1.Similarly,50 sets of particle vectors are randomly generated, and each variable will be expressed in three-dimensional coordinate points $\left(E_1,{\rho }_1,E_k\right)$, $\left(E_2,{\rho }_2,E_k\right)$,…,$\left(E_{10},{\rho }_{10},E_k\right)$, the resulting random particle coordinate points are shown in Fig. 10.

Fig. 10Initial particle distribution

Fig. 11Individual optimal particle distribution

Similarly, the optimal value among the 50 particle vectors is found by performing a local optimization search calculation, at which time $E_k=$ 2.11×10¹¹ Pa, and the individual optimal particles are expressed in two-dimensional coordinates $\left(E_1,{\rho }_1\right)$, $\left(E_2,{\rho }_2\right)$,…, $\left(E_{10},{\rho }_{10}\right)$, the results are shown in Fig. 11. The specific values of each parameter variable are shown in Table 6. Save the calculated individual optimal particles to a text file, and then import the data into ANSYS for the solution of the deflection and natural frequency, and the final calculation results are shown in Table 7.

The global optimal particle distribution is obtained by the global iterative optimization search of the particle swarm algorithm, at which time $E_k=$ 2.37×10¹¹ Pa, similarly, the individual optimal particles are represented by two-dimensional coordinates $\left(E_1,{\rho }_1\right)$, $\left(E_2,{\rho }_2\right)$,…, $\left(E_{10},{\rho }_{10}\right)$, the global optimal particle distribution is shown in Fig. 12. The specific values of the global optimal particles are shown in Table 8.

Fig. 12Global optimal particle distribution

Save the calculated global optimal particles to a text file, and then import the data into ANSYS for the solution of deflection and natural frequency, and the error is calculated with the deflection under the actual structure, and the calculated results are shown in Table 9. The error after updating is compared with the error before updating, as shown in Fig. 13.

Compare the error of the updated data of the model updating method ED-PSO model updating algorithm and CSD-PSO model updating algorithm, and the result is shown in Fig. 14.

As can be seen from Fig. 14, ED-PSO model updating algorithm is obviously superior to CSD-PSO model updating algorithm.

Fig. 13Error comparison before and after updating

Fig. 14Error comparison of two methods