Predictability of Inverse Impact Force Location as Affected by Measurement Noise

Abstract

The impact force localization inverse problem is considered through a nonlinear optimization procedure. The objective function is derived in the particular case of elastic structures for which Maxwell-Betti theorem holds. Additional geometric constraints were introduced in order to stabilize optimum search. The solution of the constrained non linear mathematical problem was performed by means of two outstanding evolutionary algorithms that include Genetic Algorithm and Particle Swarm Optimization. Focus was done on the robustness aspect of force impact localization predictability when an additive white noise is assumed to perturbed strain measurement. It was found that the Genetic Algorithm fails to track the exact solution independently from the noise level as an error was systematically present in the solution. On the other hand, the Particle Swarm Optimisation based algorithm performed very well even for noise levels as high as 2% of the measured strain signal.

1. Introduction

Identification of impact force location for impact events occurring on elastic structures can be performed by various methods that were proposed in the literature [1–3]. To review briefly some of the important contributions in this filed, Martin and Doyle [4] have described how to find the location of an impact force using dynamic response measurements. They proposed a solution procedure using the spectral element method with a stochastic iterative search. Experimentally measured acceleration responses from two frame structures were used to achieve force localization by minimizing a fitness function. A Genetic Algorithm was used to guess iteratively the minimum through monitoring the actual error associated to a given sampling generation. The process enabled to discriminate between good and bad guesses and gave at convergence the correct impact location. An alternative technique which employs the arrival time of each frequency component of a pulse detected by means of wavelet transform was proposed by Inoue et al. [5]. But this approach suffers from the lack of accuracy in measurement of small arrival times of signals. Yen and Wu [6, 7] have used multiple strain responses along with a mutuality relationship based on Green’s functions and measured strains to achieve identification of force location on two-dimensional plate-like structures. Choi and Chang [8] minimized the error between measured strain responses in PZT sensors and numerically evaluated impact force locations. Shin [9] proposed a technique for identifying the force location using modal displacements and transient signal measured by accelerometers.

As mentioned by Doyle [10], the implicit character of force localization requires considering the inverse problem of force identification as being rather associated to two decoupled subproblems: localization of impact force and reconstruction of the time force signal. The reason for this is that the inverse problem of force identification cannot be handled by just adding the force point coordinates as extra unknown parameters to the discrete vector of time force signal values. These positions intervene in fact as implicit parameters, in contrast with the force history values which are explicit. Iterations are then required for solution of force positions unknowns in order to determine them at first. Treating all force history values to be implicit parameters is not a good choice since it would increase dramatically the number of unknowns and would consequently penalize the computational cost. When the force localization is obtained through the solution of a nonlinear mathematical program, the force time signal reconstruction can then be performed by means of a regularization deconvolution technique.

In this work, separation of the localization and reconstruction phases of the impact force inverse problem is adopted. Elastic structures subjected to nonpunctual impacts for which the resulting force field can be assumed to be uniformaly distributed over a finite domain of the structure are considered. The localization problem is solved by using minimization procedures that are of evolutionary type. Two methods are examined in the following: Genetic Algorithm (GA) based strategy [11] and the Particle Swarm Optimization (PSO) approach [12, 13]. Focus will be done on the particular role related to noise affecting strain measurement as to the resulting perturbations that produce and which may impede these algorithms to converge towards the exact problem solution.

The localization problem is derived straightforwardly from the reciprocity Maxwell-Betti theorem which is valid for any elastic structure. However, as the direct expression of the fitness function formed by this theorem is ill conditioned, because it admits a lot of trivial meaningless solutions which are associated to the fixed boundary conditions, the fitness function is modified [7]. The introduced modification transforms the nonlinear mathematical program to a unique solution problem and removes the trivial parasitic solutions which complicate optimization process convergence to the real solution. Some extra constraints are also introduced in order to guide the exploration of the optimal solution; these are associated to the geometric pieces of evidence that describe the domain containing the tracked unknowns.

The aim is to examine the predictability of impact location in this situation where noise is present in observation measurement. This means determining the amount of noise that can be tolerated and also selecting between the two proposed algorithms the most suited one that can be used to conduct solution of force localization problem with sufficient accuracy and minimal error.

2. Materials and Methods

2.1. Direct Problem Formulation

Although the problem can be stated for any elastic structural system, a simplified model having the form of a beam with a rectangular uniform section is considered, Figure 1. The beam is assumed to be simply supported on both ends. It has length 𝐿, width 𝑏, and height 𝑒. It is assumed to be made from a homogeneous and isotropic elastic material with Young’s modulus 𝐸 and density 𝜌. The applied force modelling impact is assumed to result from a uniformly distributed pressure, 𝑝, which is applied on a rectangular patch as shown in Figure 1. The pressure rectangle is assumed to be centred on 𝑥₀ and of length 2𝑢.

Figure 1 

An elastic beam with a uniform rectangular cross-section and loaded with a distributed uniform pressure; four sensors positions are indicated.

The differential equation for transverse vibrations of a beam writes 𝜕𝜌𝑆²𝑦(𝑥,𝑡)𝜕𝑡²()+𝑐𝜕𝑦𝑥,𝑡𝜕𝜕𝑡+𝐸𝐼⁴𝑦(𝑥,𝑡)𝜕𝑥⁴=𝑝(𝑥,𝑡), (1) where 𝑦 is the transverse displacement, 𝑥 is the axial position, 𝑡 is the time, 𝑆 is the cross-section area, 𝜌 is the density, 𝐸 is Young’s modulus, 𝐼 is the moment of inertia, 𝑐 is the viscous damping coefficient, and 𝑝(𝑥,𝑡) is the applied pressure over [𝑥₀−𝑢,𝑥₀+𝑢].

The considered boundary conditions are as follows: 𝑦()()𝜕0,𝑡=𝑦𝐿,𝑡=0,₂𝑦𝜕𝑡²()=𝜕0,𝑡²𝑦𝜕𝑡²(𝐿,𝑡)=0. (2) Using modal superposition and Duhamel’s integral, the transient dynamic solution of the impacted beam in terms of axial strain 𝜀 at location 𝑥_𝑖 is given by the following convolution problem: 𝜀𝑥_𝑖,𝑡,𝑥₀,𝑢=^𝑡₀^𝑥0_𝑥^+𝑢_0−𝑢𝑔𝑥_𝑖,𝜁,𝑡−𝜏𝑝(𝜁,𝜏)𝑑𝜁𝑑𝜏, (3) where 𝑔(𝑥_𝑖,𝜁,𝑡−𝜏) is the time response function between point 𝜁 and the sensor location 𝑥_𝑖.

Time discretisation of (3) under the assumption that the pressure 𝑝(𝜁,𝜏) is uniform yields the following linear system: 𝑌_𝑖=𝐺_𝑖𝑥₀,𝑢𝑃, (4) where 𝑌_𝑖 designates the vector containing the discrete values of axial strain at point location indicated by index 𝑖, 𝑃 is the vector containing discrete time values of the applied uniform pressure, and 𝐺_𝑖(𝑥₀,𝑢) is the Toeplitz-like matrix.

This last has the following form: 𝑔(1)𝑔(2)𝑔(1)𝑔(𝑁)𝑔()(1)𝐺=0⋯0⋱⋮⋮⋮⋱0𝑁−1⋯𝑔, (5) where 𝑁 is the number of time steps.

An explicit expression can be obtained for matrix 𝐺. Let us denote by Δ𝑡 the time step used in discretization and by 𝑀 the order of modal truncation; the impulse response function giving the deformation of the top fiber for a section having abscissa 𝑥_𝑖 is given by 𝐺_𝑘𝑗𝑥_𝑖,𝑥₀=,𝑢4𝜋𝜌𝑆𝐿^2𝑀_𝑚=1sin𝑚𝜋𝑥₀𝐿sin𝑚𝜋𝑢𝐿sin𝑚𝜋𝑥_𝑖𝐿×𝑔𝜔_𝑚,𝜉_𝑚,𝑗−𝑘Δ𝑇, (6) with 𝑔(𝜔_𝑚,𝜉_𝑚,(𝑗−𝑘)Δ𝑇)=(𝑚exp(−𝜉_𝑚𝜔_𝑚(𝑗−𝑘)Δ𝑇)/𝜔_𝑚1−𝜉²_𝑚)sin(𝜔_𝑚1−𝜉²_𝑚(𝑗−𝑘)Δ𝑇) if 𝑗≥𝑘, and 𝑔(𝜔_𝑚,𝜉_𝑚,(𝑗−𝑘)Δ𝑇)=0 if 𝑗<𝑘, where 𝜔_𝑚=(𝑚²𝜋²/𝐿²)𝐸𝐼/𝜌𝑆 and 𝜉_𝑚 are, respectively, the circular eigenfrequency and damping ratio for a given eigenmode 𝑚.

2.2. Localization Problem as Perturbed by the Presence of Measurement Noise

The problem of finding the impact location for the beam considered in the present study consists in identifying the impact patch centre position 𝑥₀ and parameter 𝑢 defining the extent of the impacted zone. Using (4) to (6), the responses measured by strain sensors placed at points having the abscissa 𝑥_𝑖 and 𝑥_𝑗 can be expressed, respectively, under the following form: 𝑌_𝑖=𝐺_𝑖𝑥₀𝑌,𝑢𝑃,_𝑗=𝐺_𝑗𝑥₀,𝑢𝑃. (7) Equation (7) can be used to prove the following important commutativity property which does not contain the pressure vector 𝑃: 𝐺_𝑗𝑥₀,𝑢𝑌_𝑖−𝐺_𝑖𝑥₀,𝑢𝑌_𝑗=0. (8) In this way, (𝑥₀,𝑢) appears to be the solution of the following equation: 𝐺_𝑗(𝑠,𝑣)𝑌_𝑖−𝐺_𝑖(𝑠,𝑣)𝑌_𝑗²=0, (9) where ‖⋅‖ is the Euclidian norm.

There are a lot of other trivial solutions for (9) for which the factor sin(𝑚𝜋𝑠/𝐿)sin(𝑚𝜋𝑣/𝐿) appearing in the second hand side of (6) vanishes for all values of 𝑚∈[1,𝑀]. These are associated to the boundaries of the domain containing (𝑥₀,𝑢); that is, {𝑠=0;𝑠=𝐿;𝑣=0;𝑣=𝐿}. To get the exact solution of impact location, the parasitic solutions should be withdrawn from (9).

Denoting by 𝑁_𝑠 the number of sensors used, the fitness function to be minimized in order to find the impact zone parameters (𝑥₀,𝑢) is proposed under the following form: 𝑥₀,𝑢=Argmin_()𝑠,𝑣𝜙(𝑠,𝑣)=^𝑁𝑠𝑁_𝑖=1^𝑠_{𝑗=1𝑗≠𝑖}1𝛼_𝑖𝑗𝐺_𝑗(𝑠,𝑣)𝑌_𝑖−𝐺_𝑖(𝑠,𝑣)𝑌_𝑗², (10) with 𝛼_𝑖𝑗=‖𝐺_𝑖(𝑠,𝑣)‖²+‖𝐺_𝑗(𝑠,𝑣)‖² defining the weighting coefficients that are introduced in order to remove the parasitic solutions, as they vanish also for {𝑠=0;𝑠=𝐿;𝑣=0;𝑣=𝐿}.

In the presence of measurement noise, the perturbed fitness function takes the following form: 𝜙(𝑠,𝑣)=^𝑁𝑠𝑁_𝑖=1^𝑠_{𝑗=1𝑗≠𝑖}1𝛼_𝑖𝑗𝐺_𝑗(𝑠,𝑣)𝑌_𝑖−𝐺_𝑖(𝑠,𝑣)𝑌_𝑗² (11) with 𝑌_𝑖=1+𝜈_𝑖𝑧_𝑘𝑖𝑌_𝑘𝑖,𝑌_𝑗=1+𝜈_𝑗𝑧_𝑘𝑗𝑌_𝑘𝑗, (12) where 𝑧_𝑘 is a random number belonging to the interval [−1,1] and 𝜈_𝑖 and 𝜈_𝑗, respectively designate the noise level present in measurement delivered by sensor 𝑖 and sensor 𝑗, respectively.

To stabilize the minimization procedure and obtain the unique physical solution, the unconstrained mathematical program defined by (10) is constrained by adding the following geometrical bounds conditions: 𝐿40≤𝑥₀𝐿≤𝐿−,𝐿40𝐿40≤𝑢≤2. (13) To solve the mathematical program defined by (6) and (10)–(13), the PSO algorithm and the GA are considered. Their performance will then be assessed as function of the random noise level present in strain measurement.

2.3. PSO Algorithm

A PSO based method was proposed initially by Eberhart and Kennedy [12, 13]. This approach has gained since then considerable interest as being one of the most promising optimization methods that is able to provide high speed and high accuracy. PSO mimics the social behavior that a population of individuals adapts to its environment by returning to promising regions that were previously discovered [14]. This adaptation to the environment is a stochastic process that depends on both the memory of each individual, called particle, and the knowledge gained by the population, called swarm.

In the simplest numerical implementation of this method, each particle is characterised by four attributes: the position vector in the search space, the velocity vector, the best position achieved in its track, and the best position achieved by the swarm. The process steps can be outlined as follows.

Step 1. Generate the initial swarm involving 𝑁 given particles placed at random.

Step 2. Calculate the new velocity vector of each particle, based on its actual attributes.

Step 3. Calculate the new position of each particle from the current position and its new velocity vector.

Step 4. If the termination condition is satisfied, stop. Otherwise, go to Step 2.

To be more specific, the new velocity vector of the 𝑖th particle at time 𝑡+1, denoted 𝑣_𝑖^𝑡+1, is calculated according to the following Shi and Eberhart [15] formula: 𝑣_𝑖^𝑡+1=𝜔^𝑡𝑣^𝑡_𝑖+𝑐₁𝑅^𝑡₁𝑝^𝑡_𝑖−𝑥^𝑡_𝑖+𝑐₂𝑅^𝑡₂𝑝^𝑡_𝑔−𝑥^𝑡_𝑖. (14) In (14), 𝑅^𝑡₁ and 𝑅^𝑡₂ are random numbers between 0 and 1, 𝑝^𝑡_𝑖 is the best position of the 𝑖th particle in its track, and 𝑝^𝑡_𝑔 is the best position of the swarm. There are three problem dependent parameters that fix performance of this algorithm, namely, the inertia of the particle 𝜔^𝑡 and the two trust parameters 𝑐₁ and 𝑐₂.

The new position of the 𝑖th particle at time 𝑡, denoted 𝑥_𝑖^𝑡+1, is then calculated as follows: 𝑥_𝑖^𝑡+1=𝑥^𝑡_𝑖+𝑣_𝑖^𝑡+1, (15) where 𝑥^𝑡_𝑖 is the current position of the 𝑖th particle at time 𝑡. The 𝑖th particle actual position enables to determine the best position in its track 𝑝^𝑡_𝑖. When considering all the particles, the global best position of the swarm 𝑝^𝑡_𝑔 is then obtained.

PSO algorithm works such that particles concentrate on the best search position of the swarm. They cannot easily escape from the local optimal solution since the search direction vector 𝑣_𝑖^𝑡+1 calculated by (14) always includes the direction vector to the best search position of the swarm. This shows the major feature of PSO algorithm as being a robust process of continuous enhancement for optimum search.

In the presence of constraints, a particle move should be restricted in order to remain in the feasible solution space by examining the given constraints. A modified PSO version was introduced for constrained problems in order to manage this situation [16].

2.4. Genetic Algorithm

Genetic Algorithm (GA) was firstly introduced by Holland [11]. It is a probabilistic optimization method that is able to achieve global search by mimicking natural biological evolution. GA operates on a population of individuals called the set of potential solutions. Each individual is represented by an encoded string (chromosome) that contains the decision variables (genes). Traditionally, GA uses binary strings as chromosome representation.

The GA has an iterative procedure structure that comprises generally the following five main steps.

Step 1. Creating an initial population (𝐺₀).

Step 2. Evaluation of the performance of each individual or chromosome (𝑐_𝑘) of the population, by means of a fitness function to be maximized.

Step 3. Selection of individuals for the reproduction of a new population.

Step 4. Application of genetic operators: Crossover and Mutation.

Step 5. Iteration of Steps 2 to 4 until a termination criterion is fulfilled.

In the localization problem considered in this work, the candidate solution is the centre position 𝑠₀ and the extent of the impact zone 𝑢. These variables are then coded in a chromosome using a binary coding scheme.

To start the algorithm, an initial population of individuals (chromosomes) is defined. The GA is configured, so that it creates a fixed number of initial individuals at random from the whole feasible solution space. An important parameter in initialization is the population size. In general, the population size affects both the ultimate performance and the efficiency of GA and should be determined in a case by case study.

3. Results and Discussion

A pinned-pinned beam having the following material and geometric properties is considered: 𝐸=7.06×10¹⁰ Pa; 𝐿=0.5 m; 𝑏=𝑒=5×10⁻³ m; 𝜌=2660 kg·m⁻³; 𝜉_𝑚=2%; 𝑥₀=0.417 m; 𝑢=0.0417 m. The time interval considered has the duration 𝑇_𝑐=1 s. The beam is assumed to be subjected to a half sine pulse pressure having the shape and spectral content depicted in Figure 2. The maximum pressure is taken to be 10⁵ Pa.

(a)

(b)

Figure 2 

Characteristics of the impulse impact pressure; (a) time history and (b) frequency content.

The first five modes were retained. Their frequencies are given by 𝑓₁=46.72Hz, 𝑓₂=186.9Hz, 𝑓₃=420.5Hz, 𝑓₄=747.5Hz, and 𝑓₅=1168Hz. The time step used was Δ𝑡=3.4246×10⁻⁴ s which satisfies largely Shannon condition with regards to 𝑓₅.

Four gauge strain sensors were used in this study. Their labeling and positions are indicated in Table 1. Previous studies conducted by means of PSO algorithm have shown that, while considering the noise free problem, localization needs at least three sensors to be achieved with adequate accuracy. Here, a fourth sensor has been added in order to enhance the performance of the localization inverse problem in the presence of measurement noise.

Table 1 

Positions of the gauge sensors considered for the measurement of axial strain at the upper beam fiber.

Three noise levels were considered in order to test robustness of both PSO and GA based algorithms. They correspond to the following.

Case 1. Noise free case where noise amplitude is set to the value 𝜈_𝑖=0% for the four sensors.

Case 2. Noise level is fixed at 𝜈_𝑖=2% for the four sensors.

Case 3. Noise level is fixed at 𝜈_𝑖=5% for the four sensors.

The parameters used for GA were as follows:(i)stopping test: 10⁻⁶;(ii)population size: 100;(iii)probability of intersection: 1;(iv)probability of mutation: 0.05;(v)maximum number of generation: 200.

For GA algorithm the first values of the unknown parameters impact centre and extent were initialized with 𝑥₀=𝑢=1.25×10⁻².

For PSO algorithm, the following stability parameters were used: 𝑤=0.4; 𝑐₁=1.25; 𝑐₂=0.5; size population: 100.

Figure 3 gives evolution of impact location characteristics as a function of iterations in case of noiseless problem. One can see the difference existing between these two algorithms; PSO algorithm gives the exact solution of the impact location. On the opposite, GA has not converged to the right solution.

(a)

(b)

Figure 3 

Noise free case; evolution of impact zone characteristics. (a) Impact centre and (b) impact zone extent.

Figure 4 gives evolution of impact location characteristics as function of iterations in case of noisy measurement conditions with a noise level 𝜈_𝑖=2%. One can observe that GA has not converged at all, as the calculated solution is too far from the exact solution. The relative error reached 92% for the impact zone extent. Meanwhile, PSO based algorithm has continued to give the exact solution with only a moderate error, representing 13% of the impact zone extent.

(a)

(b)

Figure 4 

Noise level 𝜈𝑖=2%; evolution of impact zone characteristics. (a) Impact centre and (b) impact zone extent.

Figure 5 gives the evolution of impact location characteristics as function of iterations in case of noise level value given by 𝜈_𝑖=5%. Unexpectedly, GA is better in this condition than for noise level 𝜈_𝑖=2% as the relative error for impact zone extent has decreased from 92% to 68%. Moreover, one can see that both GA and PSO algorithms fail to give the exact solution as the error on the impact zone extent is too large, 61% and 68%, respectively.

(a)

(b)

Figure 5 

Noise level 𝜈𝑖=5%; evolution of impact zone characteristics. (a) Impact centre and (b) impact zone extent.

Tables 2 and 3 recall the obtained relative error as the function of the algorithm used and noise level. One can see that GA fails to predict the exact solution, with a behavior that is not monotonous, while PSO based algorithm performs well for small measurement noise levels not exceeding 2%. The performance decreases after that when identifying the impact zone extent, while it remains sufficiently good for the impact zone centre as the relative error is smaller than 5%.

Table 2 

Relative error affecting impact zone extent 𝑢 as compared with the exact solution.

Table 3 

Relative error affecting impact zone centre 𝑥₀ as compared with the exact solution.

To emphasize the irregular behavior of GA, a comparison has been made between the results of two runs of GA. Figure 6 shows a comparison between the obtained evolutions. They are not the same because of the stochastic nature of GA. GA is found thus to be inappropriate for the inverse force location problem as it is too sensitive to noise level and its repeatability is too poor.

(a)

(b)

Figure 6 

Noise level 𝜈𝑖=5% for the four sensors; evolution of impact centre and extent as function of iterations (a) first test GA 1 and (b) second test GA 2.

On the opposite, PSO algorithm runs always in the same manner with perfect repeatability.

4. Conclusion

Based on the separation approach that decouples force location from force signal reconstruction in an inverse impact problem occurring on an elastic beam, robustness of a particular localization procedure was analyzed. This uses a modified fitness function derived from Maxwell-Betti theorem by applying some filtering coefficients that enable to remove parasitic solutions. Solution of the obtained constrained nonlinear mathematical program that provides the impact zone location was performed by GA and PSO based algorithms. Predictability of force location was studied as a function of measurement noise level. It was found that PSO algorithm continues to achieve exact prediction of the solution even in the presence of 2% of noise intensity, while GA algorithm fails at that noise as the associated error was too large. Both of these algorithms fail to predict correct impact extent when noise level reaches 5%, even if they continue to predict the position of the impact centre with rather acceptable accuracy.

To assess predictability in a large sense, system model noise should also be integrated in the future in conjunction with measurement noise.