Application of proper orthogonal decomposition and radial basis functions for crack size estimation using particle swarm optimization

Complex engineering problems require simulations, which are computationally expensive in cases of inverse identification tasks since they commonly requires hundreds of thousands of simulations. This paper propose a method based on model reduction for crack size estimation, combining the proper orthogonal decomposition method with radial basis functions. The reduced model is validated by comparing the obtained boundary displacements with the corresponding results from a finite element model. This inverse procedure is formulated as the minimization of the difference between the measured and computed values of displacement at selected boundary nodes, called sensor points, using particle swarm optimization algorithm. Convex and a non-convex specimens have been considered for investigations of crack presence, and identification of its size, different crack sizes have been tested to demonstrate the efficiency of the proposed approach.


Introduction
Crack initiation and propagation is an omnipresent fact in all structures undergoing cyclic loads due to the fatigue phenomenon. In most cases, cracks are engaged in a predictable location. Thus maintenance measures give big importance to the crack size, trying to follow its state to prevent reaching the dangerous level. Damage detection using vibration data has been the focus of a large number of studies in the literature [1][2][3][4][5][6][7][8]. There are several numerical methods for crack detection [9][10][11], that employ different theoretical bases. Many of these methods are dedicated to the use of parameters which are not accessible experimentally.
Inverse problems are defined as the problems, in which the output is known and the input or source of output remains to be determined. They are opposite to the direct problems, in which output or response are determined using information from input [12]. In the case of the Inverse Elastostatics Problem of internal flaw detection, the geometric parameters of the flaw are unknown, but the displacements along the boundaries are known. In order to analyze this kind of problems, the boundary displacements, usually called "experimental data", are obtained under known boundary conditions and compared with the calculated ones.
Inverse crack identification problems can be stated as an optimization task [4]. There are several optimization techniques summoned in [13]. The particle swarm optimization (PSO) is a very popular algorithms used in diverse range of applications. The most utilized methods in the calculation of the mechanical behaviour of structures are boundary element method (BEM), the finite element method (FEM) [14][15][16][17][18][19][20][21][22] and recently isogeometric analysis (IGA) [23][24][25][26][27][28][29], mainly used to obtain the displacement field. The FEM and BEM was employed to solve inverse methods in structural analysis [30,31]. Generally speaking, the weak point of FEM based inverse methods is their very high computational cost. Model reduction is an alternative to solve this FEM difficulty.
The proper orthogonal decomposition (POD) is a model reduction techniques proceed by the approximation of the problem solution using the appropriate set of approximation functions [32], which contributes to the huge acceleration of the procedure since, once a trained model is built, the system response are computed by means of it, in a time shorter by about five orders of magnitude compared to FEM [33]. This leads to a very quick alternative in inverse problems, which provides simplicity and a considerably lower computational time.
Boundary measurement are employed to determine cracks since several decades [34]. The proposed crack size estimation procedure uses the same principle of classical approaches. Its main contribution lies in the way the structural response is obtained, which opposed to existing methods that employs simulation methods, they are calculated through a reduced model. Section 2 and 3 present the basic blocks of this approach, which are respectively the model reduction method POD-RBF and optimization method PSO, section 4 is dedicated to the formulation of the inverse problem by combining the POD-RBF and PSO, section 5 studies the ability of model reduction in the estimation of boundary displacement in both convex and non-convex specimen. Finally in section 6, the efficiency of this approach have been tested in identifying small, medium and large cracks in both specimens.

POD-RBF as a model reduction method
POD is a powerful statistical method for data analysis, used as model order reduction technique in different fields [35,36]. In our study, the POD was applied to determine the boundary displacement field of a two dimensional elastic structure containing an unknown crack. Used input data were finite element boundary displacements corresponding to various known cracks, called snapshots. They were first stored in matrix U, which is expressed in equation 1, as shown in figure 1.  Where N is the number of sensor points and S represents the number of snapshot vectors that represent the boundary displacement field of each crack configuration. The matrix P stores the crack parameter sets P i of all simulations, considered in our study as the crack size.
The main purpose of the POD method is to propose a set of orthogonal vectors Φ called POD basis vectors, to reassemble the snapshot matrix U in an optimal way. Φ is related to U by the following linear relationship: In Equation (2), A is the amplitude matrix collecting the coefficients of the new basis combination and, according to the orthogonality of Φ , it can be computed from: Optimal basis vectors are defined by the performance of the POD method, also known as the singular value decomposition operation [37,38]: Where V is the matrix storing the normalized eigenvectors of the covariance matrix C , and Λ the diagonal matrix storing its eigenvalues. The matrix C is given by the following equation: Due to the optimality of the new system Φ constructed as a POD basis, a low dimensional approximation Φ of high accuracy is extracted from it by preserving only K (K ≪ S) columns that correspond to the largest eigenvalues. Since the eigenvalues of the covariance matrix C , called the energy of the system, are stored in a descending order, POD directions that hold little information are then discarded without influencing the accuracy of the representation. This is known as the truncation of the POD basis, and is accomplished by choosing the fraction of system energy that will be neglected in later calculations. Consequently, the amplitude matrix Â is given by: To determine the boundary displacement field of a two dimensional elastic structure containing an unknown crack, RBF interpolation was used. This method can generate different sets of parameters, which were not included in the initial selection in the matrix P. The amplitude matrix A is defined as a multiplicative form of the function G, defined as the matrix of interpolation parameters, and the matrix B containing the unknown coefficients: The interpolation functions are expressed by [32,37,39]: P i is the parameter corresponding to U i (i=1,2,…,S). The argument of the i-th RBF is the distance |P-P i |, P and Pi being respectively current and reference parameters. c is the RBF smoothing factor defined in the range from 0 to 1. If all or some of the the knot points P i are relatively close to each other, the matrix G could be singular, which is circumvented by reducing the c value. After the evaluation of the coefficient matrix B , a low-dimensional model issued from (8) is written in the following vector form: a(P) = B ⋅ g(P), (10) Equation (7) can be expressed as an approximation of the snapshot u corresponding to a new parameter vector P : This model will now be referred as the trained POD-RBF network. It is capable of reproduce the unknown boundary displacement field of the structure that corresponds to any set of crack parameters P . It must be noted that extrapolation outside the range of P leads probably to poor precision of the model. Increasing the value of the smoothing factor c leads to a better interpolation. But it can make the matrix G singular, depending on the closeness of knot points. In the present work, the parameter c was chosen to be constant for all functions, and equal to the mean value of normalized parameters.

Particle swarm optimization
The Particle Swarm Optimization (PSO) is a population based optimization method inspired from the behaviour of bird flocks that is characterized by distinct social and psychological principles. Large attention has been paid to this method in few last decades. The algorithm was initially proposed by Kennedy and Eberhart [40]. Its implementation requires a small number of parameters, which facilitates its application and reduces the computational cost.
The main idea of PSO is to consider the potential solution as a particle moving through the space, looking for the global optimum position. Initiated as a group of random particles, each particle is characterized by its position in the multidimensional space and by its movement speed. These particles, each other, cooperate to achieve the solution, based on their personal previous experience and the experience of other particles. The speed and the position of the particles are calculated, respectively, as follows: The weight inertia w is multiplied by the particle speed value at every iteration to maintain the particle acceleration in its original direction. c 1 is a positive constant, called cognitive parameter and controlling the step size toward the particle's personal best position. c 2 is social parameter that controls the step size toward the global best position. {r 1 } and {r 2 } are vectors containing random numbers within the interval [0,1]. {x j (t)} is the vector of the current positions of particles. {x Pb,j } is the vector of the personal best position found by the particle j and {x Gb } is the vector of the global best position found by the entire swarm.

Inverse problem formulation
The existence of a crack changes the behavior of the plate when put under traction, therefore the deformation of the structure, which is also affected by the changes of the crack length. Benefiting from this fact, the deformation of the structure's border is measured using deformation sensors for inverse crack size estimation.
The crack estimation algorithm consists of two main stages. In the first stage, the identification problem is defined, and the response data corresponding to the unknown crack is chosen. In the second stage, the optimization algorithm is executed. All fitness function values are obtained from calculation on the reduced model, unlike classical methods where this operation needs a full analysis of the whole structure. Figure 2 summarizes the proposed approach. The fitness function was evaluated from the following equation: By introducing the crack parameters, corresponding to each possible solution, in the trained POD-RBF network, the resulting boundary displacement vector u(P) is generated. Then, the fitness function value is the norm error between this vector and the reference displacement vector u(P 0 ) caused by the real crack parameters. The PSO algorithm is operating in the following manner: Step 1. Initialization of the algorithm by randomly generating the position vectors of the particles within the design space (x 1 ⃗⃗⃗ , x 2 ⃗⃗⃗⃗ , … , x NP ⃗⃗⃗⃗⃗⃗⃗ ) and calculating their corresponding speed vectors (v 1 ⃗⃗⃗ , v 2 ⃗⃗⃗⃗ , … , v NP ⃗⃗⃗⃗⃗⃗⃗ ), where NP the number of particles.
Step 2. Analysis and evaluation of the fitness value for current positions of the particles ( (x 1 ⃗⃗⃗ ), (x 2 ⃗⃗⃗⃗ ), … , (x NP ⃗⃗⃗⃗⃗⃗⃗ )) , where NP the number of particles. Equation (14) in our study.
Step 3. Personal best calculation, if the current fitness value is better than the best fitness value (Pbest) in the particle's history then this value is set as the new Pbest and the current position as the new x Pb,j for each j particle.
Step 4. Global best calculation, the best fitness value of all the particles Pbest is set as Gbest and the corresponding position as the new x Gb .
Step 6. Check for any dimension i x i ≤ x L or x i ≥ x U set x i = x L or x i = x U respectively and ν i = 0.
Step 7. If the maximum number of generations or a defined fitness level is reached, the algorithm is terminated; else, the steps 2-6 are repeated.

POD-RBF for the computation of boundary displacement field
The POD in this stage has been used to build a reduced model describing the effect of crack size on the vertical boundary displacements, this model is based on a snapshot matrix containing boundary displacement vectors collected from 11 scenarios, corresponding to crack length s varying in the range: 0 (no crack) to 2.5 mm, in two structures where the largest value of the crack (2.5 mm) is half the width of the specimen.
The first specimen (Specimen 1) consists of a rectangular plate, subjected to traction force from the upper and lower sides simulated as a to displacement of 0.1 mm. Young modulus and Poisson coefficient of the material are E = 210 GPa and ν = 0.3, respectively. The crack is located at the center of the plate. The displacement data considered for validating the POD-RBF model is obtained from the higher half of the vertical boundary as describes in Figure 3, which contains 40 nodes. The second (Specimen 2) specimen is a non-convex structure, with same material properties and subjected to the boundary conditions as described in Figure 4. The displacement results obtained from the 49 boundary nodes, highlighted by red color in descriptive figure, are considered for validating the POD-RBF. The choice of the best K value is important, as a too low value will lead to poor precision, and too large value will cost more computational time. To test the effect of the chosen truncation point K on the accuracy of the model, we compared the POD-RBF results, corresponding to crack with size equal to 1.4 mm, with equivalent results from FEM. Figure 5 and Figure 6 study respectively the convex and the non-convex examples by means of the ratio error of the boundary displacement field, produced by POD-RBF model, issue from models constructed based on different K values, and equivalent boundary displacement field from FEM.  From Figure 5, we notice that maximum error for all K values is located at node number 40. This can be explained by the fact that this node is located on the centre of the structure, where the displacement is equal to zero. Therefore, a small difference between estimated results and FEM results affect largely the error value (POD-RBF/FEM). We can also see that the largest error value for all K values is about 2% found for K=1 and K=2. On the other hand, more than 4 modes ( ≥ 4) have led to better results with error around 1%. Moreover, when the first three modes are used, the results were even better with error less than 0.5%. Figure 6 shows that node 49 has the largest error value for all modes, similarly to the convex problem, this is due to the fact that displacement value in this node is close to null. It is also noticed that error in the interval between nodes 21 and 31 has a considerable error values for K=1 and K=2, which is not present for K values larger than 3 with error less than 0.1%. Therefore, the first three modes are chosen to represent the full model in the inverse calculations.

Crack size estimation
The inverse problem solved by PSO, minimizing the cost function, which is the norm error between the vertical displacement vector caused by the crack that we want to estimate and the one proposed by PSO.
The displacement data considered in the fitness evaluation is obtained from only 4 sensor point in both cases, convex and non-convex, where their position, as described in Figures 3 and 4, are chosen based on early study [40], in which, it has been found that better results are acquired when the sensors are dispersed uniformly respecting a near distance between every sensor point. Tables 1 summarize the PSO parameters, and Table 2 and 3 present crack size estimation results, for convex and non-convex specimen, respectively, the results presented in both table represent the best over 10 runs. In each table, three levels of crack size are studied (small, medium and large), as well as, the case of the absence of the crack, which is represented by crack size value equal to zero. The fitness value of each case is shown along with the error of the estimated results.  Tables 2 and 3, we can see that the presented approach could estimate the absence of the crack in both specimens with high accuracy, leading, for specimen 1, to a very small crack of 0.04 mm. It is noted that the fitness value does not follow any order. This is due to using only four sensors, as the displacement results in those sensors varies from one crack size to another. For the three levels of crack sizes, the maximum error value is found in the smallest crack (0.35mm), with a difference between the estimated and real crack equal to 0.02 mm, not the error is large because is calculated compared to a From both tables, we can see the estimation results are very close, which indicates that geometry of the specimen doesn't affect the inverse estimation method. The average computation time in PC with Intel Dual-Core Processor 3.0 GHz and 2 GB RAM, is about 4 seconds.
For the sake of illustration, figures 7 and 8 depict the crack size evolution and the fitness convergence, for the crack size equal to 1.4 in the case of specimen 2.   8 shows that the algorithm approached the final result very early, i.e. before the 10th iteration, and reach the optimal result around the 40 th iteration. This means that good precision can be conserved even though a low number of iteration is adopted as stopping criteria.

Conclusion
In this numerical study, we presented POD-RBF method as model reduction of cracked specimens under traction. The finite element model of the structure was created for different crack lengths, and the results were obtained from the new model was compared to the original ones to insure the efficiency. Crack size was investigated based on vertical boundary displacement data using the PSO for the inverse calculation. The results have clearly shown that the developed algorithm was capable of estimating crack size accurately, and proved its effectiveness even with a very low number of sensors, in both convex and non-convex structures. The crack size estimation using PSO showed high accuracy, even with a very low number of sensors, and showed that the inverse calculation on the reduced model by POD-RBF was practical and helpful in avoiding computational time problems typical for the simulation based inverse problems.