Multi-Objective Particle Swarm Optimization of Sensor Distribution Scheme with Consideration of the Accuracy and the Robustness for Deformation Reconstruction

For the inverse finite element method (iFEM), an inappropriate scheme of strain senor distribution would cause severe degradation of the deformation reconstruction accuracy. The robustness of the strain–displacement transfer relationship and the accuracy of reconstruction displacement are the two key factors of reconstruction accuracy. Previous research studies have been focused on single-objective optimization for the robustness of the strain–displacement transfer relationship. However, researchers found that it was difficult to reach a mutual balance between robustness and accuracy using single-objective optimization. In order to solve this problem, a bi-objective optimal model for the scheme of sensor distribution was proposed for this paper, where multi-objective particle swarm optimization (MOPSO) was employed to optimize the robustness and the accuracy. Initially, a hollow circular beam subjected to various loads was used as a case to perform the static analysis. Next, the optimization model was established and two different schemes of strain sensor were obtained correspondingly. Finally, the proposed schemes were successfully implemented in both the simulation calculation and the experiment test. It was found that the results from the proposed optimization model in this paper proved to be a promising tool for the selection of the scheme of strain sensor distribution.


Introduction
During the past few decades, deformation reconstruction has played an important role in real-time structural health monitoring (SHM) with the flourish of civil engineering, aerospace, and smart structure development. Shape sensing techniques can be employed to reconstruct the displacement field based on strain sensor values in a real-time manner [1,2]. The advantage of this method is that a priori knowledge, such as model shape, load form, and elastic-inertial material information, is not required in the process of reconstruction. This is particularly the case in scenarios where the applied loads are difficult to be accurately determined or measured. For example, the applied loads are the aerodynamic forces, temperatures, impact loads, or vibrating excitation [3].
The key technology in shape sensing was to accurately establish the mathematical model between the discrete surface strain and the displacement field. Previous researchers proposed many modeling methods. The adaptive fuzzy net was employed to establish the transformation relationship between
Based on the small-strain hypothesis, the strains can be obtained from the deviation of displacement based on the Timoshenko theory along the x axes, as given by Equation (1): ε x (x, y, z) = e 1 (x) + ze 2 (x) + ye 3 (x) γ xz (x, y) = e 4 (x) + ye 6 (x) γ xy (x, y) = e 5 (x) − ze 6 (x) (1) where e(u) = {e 1 , e 2 , e 3 , e 4 , e 5 , e 6 } describes the section strain of the theory, and it can be presented further in terms of nodal displacement as follows: where B i (i = 1, . . . ,6) is the deviation of the displacement shape function. u e contains the nodal displacement and rotation angle. The iFEM uses the weighted least square principle to minimize the difference between in situ section strain, e ε , and analytic section strains, e(u), which can be written as follows: The relationship between in situ strain measurement and deformed beam displacement can be obtained by substituting Equation (2) into Equation (3), which can be expressed as follows: k e u e =f e (4) The two terms k e and f e are defined as follows: with where n denotes the section number, x i is the section location along the x axis, and e ε i (x i ) is the in situ section strains. L denotes the element length.w i (k = 1,2, . . . ,6) are positive-valued constants for the axial stretching, bending, twisting, and transverse shearing, respectively, and their initial values have been set as 1 in this paper.
Once the scheme of strain sensor distribution has been determined, which means that the values of x = x i , θ = θ i , β = β i (i = 1,2, . . . ,6) are determined, the section strains can be calculated from the measured surface strains using strain-tensor transformations as follows: where µ is the Poisson ratio. R is the external radius of the section as shown in Figure 1.
where  is the Poisson ratio. R is the external radius of the section as shown in Figure 1. Then, based on the analysis and calculation results of Equation (12) in Reference [23], these unknown parameters (a1,..., a6), which are used to established the relationship between section strains and the position along the centroid axis, can be solved:  ()   Q T Q ε (9) where the Q matrix can be obtained by replacing variable xi by yi in the * Q matrix, and yi denotes an arbitrary point placement along the centroid axis direction. The transfer relationship between deformation displacement and surface strain can be constructed by substituting Equation (9) into Equation (4). This relationship can be expressed as follows: Then, based on the analysis and calculation results of Equation (12) in Reference [23], these unknown parameters (a 1 , . . . , a 6 ), which are used to established the relationship between section strains and the position along the centroid axis, can be solved: [a 1 , a 2 , a 3 , a 4 , a 5 , . . , T 6 ] T , i = 1,2, . . . ,6. Therefore, for an arbitrary point on the centroid axis, the section strains (e 1 ,e 2 , . . . e 6 ) can be determined by substituting Equation (8) into Equation (7): where the Q matrix can be obtained by replacing variable x i by y i in the Q * matrix, and y i denotes an arbitrary point placement along the centroid axis direction. The transfer relationship between deformation displacement and surface strain can be constructed by substituting Equation (9) into Equation (4). This relationship can be expressed as follows: where

Establishing the Multi-Objective Particle Swarm Optimization Model
The section strains can be computed from the strain sensor values which can be measured by Fiber Bragg Grating (FBG) sensors. However, during the FBG sensor installation, installation errors occur between the installation position and the theoretical position of the sensor along the axial length x i , angle parameter θ i, and β i . Therefore, these position parameters with errors are substituted into Equation (7), which will affect the accuracy of the computed section strain. Moreover, in Equation (10), the transfer relationship (TR) is the matrix about the scheme of strain sensor distribution, and acts as a bridge between the surface strain and the deformation displacement. However, it was found that an inappropriate strain sensor distribution scheme could lead to a singular or ill-condition TR matrix [23]. For example, a tiny surface measured strain change might cause a significant reconstructed displacement change using the iFEM, which further deteriorates the accuracy or even leads to the process failure. Thus, the robustness of transfer relationship is one optimization objective ( f 1 ), which can be constructed with well-separated eigenvalues [27]: where λ i represents the eigenvalues of transformation matrix TR(ξ i , θ i , β i ). Parameter ξ i represents the sensor location along the centroid axis direction. It is difficult to determine the strain sensor locations on both end nodes in engineering applications. Therefore, the installed placements are selected in the range of ξ i ∈ [ L 10 , 9L 10 ]. In order to minimize the positioning errors for a strain sensor distribution scheme, one sensor is placed at β i = 45 • , whereas the other five sensors are located at β i = 0 [23].
The accuracy is the key reference to estimate the performance of deformation reconstruction using the iFEM. Generally, structural beam deformation is varied due to different load forms. In order to eliminate the effect of working conditions on deformation, the relative root mean square (RRMS) is proposed to act as another optimization objective function ( f 2 ) for estimating reconstruction accuracy, which can be expressed as follows: where disp(x i ) is the deformation displacement along the centroid axis in one direction. The superscript measurement denotes the actual deformation value from the simulation software or a measurement device. The iFEM is the predicted value computed from the iFEM. RMS indicates deformation reconstruction root-mean-square errors. RRMS is the ratio of RMS and maximum deformation value from the actual deformation value.
The purpose of strain sensor distribution scheme optimization is to find an optimal scheme, which can make the optimal balance between robustness and accuracy throughout a limited beam surface measurement. Thus, the optimization model can be described as follows: For multiple objective optimization problems, Coello et al. proposed an extending particle swarm optimization (PSO) approach to solve this problem, referred to as an MOPSO [28]. Compared with other algorithms, the advantages of this method are a high convergence speed and a simple algorithm structure [29].
According to the PSO algorithm theory [30][31][32], a population size is 50 and randomly initialized in the limited beam surface. The maximum iteration number is 100. Each particle position can be expressed by a vector 1×36 , which determines the direction that particle can move for the search of the optimal solutions. Then, each particle is introduced into an optimal model to calculate a potential solution (13). The optimal solutions are continuously searched by updating the velocity and position of particle according to Equation (14): where t indicates the tth iteration in the PSO algorithms. Parameter w denotes inertia weight, which is decreased from 0.9 to 0.4. Acceleration constant parameters are c 1 = c 2 = 2. Parameters r 1 and r 2 are random values uniformly distributed in [0,1]. Parameters p i and p g indicate the element of personal best solution (pbest) and global best solution(gbest). The maximum velocity (v max ) is equal to the dynamic range in each dimension of the particle, and v t+1 . PSO will be stopped if the tth iteration satisfies the predefined maximum number of iterations. The concrete specifications of this optimization algorithm are explained in Figure 2.
swarm optimization (PSO) approach to solve this problem, referred to as an MOPSO [28]. Compared with other algorithms, the advantages of this method are a high convergence speed and a simple algorithm structure [29].
According to the PSO algorithm theory [30][31][32], a population size is 50 and randomly initialized in the limited beam surface. The maximum iteration number is 100. which determines the direction that particle can move for the search of the optimal solutions. Then, each particle is introduced into an optimal model to calculate a potential solution (13). The optimal solutions are continuously searched by updating the velocity and position of particle according to Equation (14): where t indicates the tth iteration in the PSO algorithms. Parameter w denotes inertia weight, which is decreased from 0.9 to 0.4. Acceleration constant parameters are 1   The pbest is continuously updated in the following process: in the search space, if a particle searches a new solution that dominates its pbest, the new solution can be selected as the new pbest. Otherwise, the old solution has a probability of 0.5 to be replaced by a new solution. The criterion of crowding distance [29], which is used to estimate the density of solutions surrounding a particular solution in the Pareto front, is proposed to select gbest. For example, when the value of crowding distance is large, the solution in the Pareto front will be deviated from its surrounding solutions. Additionally, the reflecting wall proposed in Reference [30] is used to improve the exploration of the particle.

Simulation Calculation
In the simulation, the aluminum hollow beam using 188 elements was modeled based on the Timoshenko beam theory. The beam's Young's modulus is E = 73000 MPa with the Poisson ratio of v = 0.3 and a density of p = 2712.63 kg/m 3 . The beam length is L = 660 mm, with an external diameter R = 13 mm as well as an inner diameter r = 11.5 mm. The model was discreted into 200 elements and the cross-section of the beam was discreted into 120 segments as shown in Figure 3.
The free-end of the cantilever beam was subjected to static loads with a combination of two forces in three conditions (Table 1) for obtaining these parameters in nodes, containing deformation displacements, rotation angles, strains, and shear strains. Then, based on the obtained data, the MOPSO was established to optimize the strain sensor distribution scheme. The Pareto front and bi-objective function values ( f 1 , f 2 ) are shown in Figure 4.
crowding distance [29], which is used to estimate the density of solutions surrounding a particular solution in the Pareto front, is proposed to select gbest. For example, when the value of crowding distance is large, the solution in the Pareto front will be deviated from its surrounding solutions. Additionally, the reflecting wall proposed in Reference [30] is used to improve the exploration of the particle.

Simulation Calculation
In the simulation, the aluminum hollow beam using 188 elements was modeled based on the Timoshenko beam theory. The beam's Young's modulus is E = 73000 MPa with the Poisson ratio of v = 0.3 and a density of p = 2712.63 kg/m 3 . The beam length is L = 660 mm, with an external diameter R = 13 mm as well as an inner diameter r = 11.5 mm. The model was discreted into 200 elements and the cross-section of the beam was discreted into 120 segments as shown in Figure 3.
The free-end of the cantilever beam was subjected to static loads with a combination of two forces in three conditions (Table 1) for obtaining these parameters in nodes, containing deformation displacements, rotation angles, strains, and shear strains. Then, based on the obtained data, the MOPSO was established to optimize the strain sensor distribution scheme. The Pareto front and bi-objective function values ( 12 , ff ) are shown in Figure 4.  According to the Pareto front as shown in Figure 4a, the better front can be found in a small range. Figure 4b describes the objective function values of all positions visited by the swarm during optimization, which can be divided into the accurate region and the stable region, namely, C1 region and C2 region. However, in engineering, taking into consideration the strain senor distribution location errors, a good scheme should be searched in the C2 region, where a high accuracy and good robustness can be identified. On the other hand, the higher accuracy C1 region was also selected to compare the robustness of the C2 region. The specific strain sensor distribution schemes from ranges C1 and C2 are described in Table 2.  According to the Pareto front as shown in Figure 4a, the better front can be found in a small range. Figure 4b describes the objective function values of all positions visited by the swarm during optimization, which can be divided into the accurate region and the stable region, namely, C1 region and C2 region. However, in engineering, taking into consideration the strain senor distribution location errors, a good scheme should be searched in the C2 region, where a high accuracy and good robustness can be identified. On the other hand, the higher accuracy C1 region was also selected to compare the robustness of the C2 region. The specific strain sensor distribution schemes from ranges C1 and C2 are described in Table 2.    The C1 and C2 schemes were tested on the same hollow beam which was subjected to static loads with a two-force combination. The comparisons of iFEM deformation reconstruction of the hollow beam with two different schemes are presented in Table 3. The C1 and C2 schemes were tested on the same hollow beam which was subjected to static loads with a two-force combination. The comparisons of iFEM deformation reconstruction of the hollow beam with two different schemes are presented in Table 3. Table 3. Comparison of reconstruction accuracy under two schemes (mm). In this paper, the estimation criteria of reconstruction accuracy can be defined as follows:

Heavy Load (N) Middle Load (N) Slight Load (N)
where MER denotes the maximum errors. Max-disp is the maximum deformation displacement in theory. Table 3 shows that the reconstruction accuracy of the C1 scheme is higher than C2 on the condition of not considering the position errors of strain sensor distribution.
The robustness and the accuracy of the C1 and C2 schemes are further compared with the consideration of these errors from the strain sensor positions and measuring equipment, which may affect the accuracy of deformation reconstruction computed from Equation (10). It is assumed that these errors are as follows: It is assumed that these noises are randomly distributed in the set, and these tests are performed 500 times in total. Then, the highest error can be selected for the reconstruction accuracy, which is shown in Tables 4 and 5. Table 4. Comparison of the maximum deformation errors with adding disturbance (mm).

Heavy Load (N)
Middle Load (N) Slight Load (N)  Table 5. Comparison of the maximum RRMS with adding disturbance (%).  Table 3 shows that the C1 scheme has a higher deformation reconstruction accuracy than C2 using iFEM without considering these disturbances. In particular, for the Y direction in heavy load, the error of the C2 scheme is 0.19 mm. On the other hand, the error of the C1 scheme drops to 0.08 mm.

Heavy Load (N) Middle Load (N) Slight Load (N)
When the noise is introduced, the maximum errors in the X, Y, Z directions are 1.01 mm, 20.6 mm, and 23.7 mm, respectively, for the C1 scheme. On the other hand, the maximum errors in the X, Y, Z directions drop to 2.64 mm, 2.61 mm, and 2.71 mm, respectively, for the C2 scheme, as shown in Table 4. The RRMS in the Y and Z directions are within 34% for the C1 scheme, and the RRMS in the Y and Z directions are stable within 4.8% for the C2 scheme, as shown in Table 5. It is indicated that the accuracy and the robustness of the C2 scheme is better than the C1 scheme. Figure 5 shows the experiment rig and its measurement equipment, which includes a 3-D optical measurement device (NDI Optrotrak Certus, Northern Digital INC, Waterloo, ON, Canada) and Fiber Bragg Grating (FBG) strain sensors. The two-strain sensor distribution schemes and six-position sensors are used for the deformation reconstruction of the hollow beam under the end-node loadings, which are presented in Table 6.

Experiment Verification
In the experiments, the surface strains are measured by two schemes of the FBG strain sensor, respectively. In addition, the beam shapes are recorded by the NDI (in Figure 5b), which can accurately capture the position of the position sensors in Figure 5c. The comparison between the measured deformation captured from the NDI and the reconstructed deformation computed using the iFEM is shown in Table 7. Moreover, the deformation displacement values in measurement points are also shown in Tables 8-10.
From Table 7, MD NDI denotes the maximum deformation measured by NDI. MDiFEMC1 and MDiFEMC2 denote the maximum deformation computed using the iFEM with the sensor distribution C1 and C2 schemes. MER C1 (MER C2 ) and RRMS C1 (RRMS C2 ) can be computed from Equations (17) and (20), respectively. Table 7 compares the results of reconstructed displacement using the iFEM with the deformation displacement measured with NDI. It is found that the C2 scheme result gives a better robustness and a higher accuracy than the C1 scheme in most loading scenarios. In particular, the maximum errors (MER) in the Y and Z directions are 2.57 mm and 3.95 mm, and 0.2 mm and 0.16 mm, respectively, in most loading scenarios. The corresponding relative root mean square (RRMS) values are 22.8% and 20.9%, and 3.1% and 2.6%, respectively. However, in the X direction, the deformation displacement should be zero in theory, but there are slight deviations due to the imperfect distribution of the actual end-load in the YOZ plane. Under a heavy load, the deformation displacement of the hollow beam is only 0.49 mm, which gives the same magnitude order with the 3-D measurement device accuracy (NDI). Therefore, the maximum error and relative root mean square (RRMS) of reconstructed deformation displacement is comparatively high. Based on the aforementioned discussion, the optimal model proposed in this paper can allow the selection of a proper strain sensor distribution scheme which could result in a strain-displacement transfer relationship with better robustness and higher reconstructed deformation displacement accuracy.  Figure 5 shows the experiment rig and its measurement equipment, which includes a 3-D optical measurement device (NDI Optrotrak Certus, Northern Digital INC, Waterloo, ON, Canada) and Fiber Bragg Grating (FBG) strain sensors. The two-strain sensor distribution schemes and six-position sensors are used for the deformation reconstruction of the hollow beam under the end-node loadings, which are presented in Table 6.

Heavy Load (N) Middle Load (N) Slight Load (N)
Load (kg) 9 6 3 In the experiments, the surface strains are measured by two schemes of the FBG strain sensor, respectively. In addition, the beam shapes are recorded by the NDI (in Figure 5b), which can accurately capture the position of the position sensors in Figure 5c. The comparison between the measured deformation captured from the NDI and the reconstructed deformation computed using the iFEM is shown in Table 7. Moreover, the deformation displacement values in measurement points are also shown in Tables 8, 9, and 10.

Conclusions
The robustness of the strain-displacement relationship and the accuracy of the reconstructed deformation displacement are two key factors to evaluate the strain sensor distribution scheme. In this paper, the MOPSO was employed to reach a balance between robustness and accuracy. In the algorithm, the well-separated eigenvalues of transfer matrix was employed as an optimal objection, and the relative root mean square of reconstructed deformation displacement was used as another optimal objection. The optimization model was then tested on a hollow cantilever beam subjected to various loads. In the test, the Pareto front and the objective function values of all visited positions were obtained. Next, two schemes were chosen from the solution space based on practical application conditions and