Inclusion Identification by Inverse Application of Boundary Element Method, Genetic Algorithm and Conjugate Gradient Method

Characterization of the interior of an inhomogeneous body using displacement measurements obtained by conducting a simple tension test was investigated. A homogeneous elastic solid body is assumed to have another solid body with arbitrary shape hidden inside it. The shape and physical properties of this inclusion were unknown. The Boundary Element Method (BEM) coupled with a mete heuristic Genetic Algorithm (GA) and Conjugate Gradient Method (CGM) was used in this characterization problem. A fitness function, which is the summation of squared differences between the measured displacements and calculated displacements at the same locations on the boundary, is minimized using GA and CGM. GA is used to find a good initial estimate of the unknown parameters and the CGM was used as the hybrid function to get converged values of the unknown parameters. For the cases that the fitness function fluctuates severely, a regularization function was added to CGM in order to dampen the fluctuations.


INTRODUCTION
In this study, the feasibility of characterizing the internal structure of an inhomogeneous body solely on the basis of surface measurements will be explored. In particular we will examine the problem of determining the location, shape and physical properties i.e., the Young's modulus and Poisson's ratio of an inclusion in a body of arbitrary shape using measurements of displacement at a discrete number of points on the surface of the body only.
Application of the BEM to the inverse elasticity problems is not new. Estimation of the shape of an inclusion in two-dimensional region using impedance measured at the domain surface [1] , estimation of the diffusion constant of a homogeneous body as well as unknown internal heat sources utilizing surface temperature measurements [2] , estimation of the optimal sizes and locations of coolant flow passages for a userspecified steady distribution of surface temperatures and fluxes [3] , investigation of a spatially regularized solution of inverse elasticity problem using the BEM [4] , estimation of thermal and mechanical properties along with location and size of a circular inclusion [5] , nondestructive detection of cavities by an inverse elastostatics BEM [6] , damage detection and assessment of structures from static responses [7] , estimation of the geometric and material properties of an inclusion using BEM and implementation of regularization function [8] , Local optimization techniques in inverse elasticity problems using BEM [9] , estimation of geometric shape of an irregular boundary between a numbers of solid elastic bodies using surface temperatures [10] , are among some of the investigations done in this field.
In the present study, a body of some arbitrary but given shape is subjected to uniaxial tension. The body is assumed to contain an inclusion of unknown shape, location and elastic constants i.e. Young's modulus and Poisson's ratio. The simultaneous estimation of these parameters is to be accomplished by displacement measurement at surface locations. Since this inverse problem is highly ill posed, the GA and CGM are effective methods used to get the optimized solution.
Several questions arise in this inverse application of the BEM, which have not been adequately addressed in previous investigations. Among these are (1) what effect will the inevitable errors in experimental measurements have on the ability to estimate the sought parameters and (2) what is the influence of the inclusion size on the estimation process?
The present study assumes a two-dimensional body containing a single inclusion, but the method of analysis is not limited in principle. Extension to threedimensional problems and more complex internal structures is straightforward but the additional question of how many parameters one can hope to determine simultaneously requires further study. Nonetheless, the success of the technique demonstrated here shows great promise for application in the area of non-destructive evaluation.

MATERIALS AND METHODS
Direct problem: A linear elastic solid of uniform thickness h, loaded in some manner, is considered. The equations of equilibrium in terms of displacements are: For an isotropic material: ( ) ( ) for i = 1,2 , where u i and t i are the components of displacement and traction vectors in the i direction respectively. f i (s) and g i (s) are prescribed functions and s is a coordinate along the boundary as shown in Fig. 1. The displacement at any point x i in Ω can be determined from [11] as: where, T kj and U kj are known functions. smooth at x i and δ kj is the Kronecker delta. Subdividing the boundary Γ into N segments, Eq. 5 becomes: Linear elements are employed, so Eq. 6 is written as follows.
Note that the displacements are assumed to be piecewise linear and continuous whereas tractions are assumed to be piecewise linear and discontinuous. Here  [11] . Equation 7 is written in matrix form as: For an inhomogeneous body, the domain Ω is divided into two subdomains Ω 1 and Ω 2 , each having its own Poisson's ratio and Young's modulus and Eq. 8 reapplied to each subdomain. For two subdomains: where, continuity of traction's at the interface nodes is assumed. Finally after all boundary conditions are imposed, the system of equations is written in the form: Which is solved for {X}. The column vector {X} contains the values of unknown displacements and tractions at the outer boundary nodes as well as all displacements and tractions at the interface boundary nodes. A number of examples have been considered and solved using the BEM [5] .
Inverse problem: In the previous section the use of BEM to solve the direct linear elasticity problems for homogeneous and inhomogeneous bodies were demonstrated. In the direct problems, the governing equation, geometry, material properties and the boundary conditions are given and the unknown boundary data is computed. In the inverse problem the geometric features and material properties are unknown, but some of the unknown boundary data can be measured and used as additional information necessary to estimate the unknown input parameters. In this section, a body that is known to contain an inclusion with unknown shape and mechanical properties is considered. These where, E is the Young's modulus of the inclusion, ν is the Poisson's ratio of the inclusion, (x c , y c ) are the coordinates of the center of the inclusion and r i are the radial distance of n nodes on the boundary of inclusion to (x c , y c ) as shown in Fig. 2. The sum of the squared differences between measured and computed displacements is: where, « ∧ » denotes the estimated values of the unknown parameters.
Genetic algorithm: Genetic algorithm is a practical method of solving optimization problems on the basis of natural genetics. It is designed for large domain, nonlinear, discontinuous and not well defined solution domains where there is not enough information or modeling is difficult. It is also applicable to the problems where traditional estimation and optimization methods are not appropriate [12] .
In GA, each person is a possible solution to the problem, which include array of possible responses. Since the material constant (E) has broad range of values and the geometric parameters must be within a specified range, the reproduction operators must be defined specifically in order to get a correct solution to this problem.

Initial population:
The initial population should be chosen such that all the unknown parameters i.e. Eq. 12 are represented. For this reason, chromosomes which include Young's modulus, Poisson's ratio and x-y coordinates of the center of inclusion and a number of radiuses needed to estimate the shape of the inclusion are chosen as the initial population. The first four genes don't have any limit in their values and after being normalized, they pick a number between zero and one. But the radius should be such that the boundary of the inclusion remains within the domain of the problem. Therefore if the inclusion is assumed to be inside a rectangular domain with dimensions of a and b as shown in Fig. 2, then the radius should be bounded by the following equations: If the radius is within the above limits, it will be chosen as the candidate answer; otherwise another random value for radius will be tested.

Crossover operator:
The crossover operator is the most important tool in GA, which produces new solutions (individual). It chooses randomly two solutions from the current population as parents. They could reproduce and make new offspring, which are better than parents. If the offspring have advantage over their parents, they remain in the population and can be chosen to reproduce. If offspring don't fit within the problem's geometry or the physical properties of the inclusion, then they are kicked out of the population and would not be reconsidered again. The crossover operator is chosen according to the following criteria: where, R is a random vector with dimension equal to the number of parameters (genes), i.e. dimension of [x] in Eq. 12. It doesn't have any limit for the first four genes, but after the forth gene, i.e. the radiuses, this random number should be such that it specifies the geometric limitations. An important factor in crossover operator is the method of defining crossover fraction (p c ). It specifies the fraction of the next generation, other than elite individuals, which is produced by crossover. The remaining individuals, other than elite, in the next generation are produced by mutation. Crossover fraction gets a value between 0 and 1. With regard to this problem and the nature of parameters, p c is defined as a variable with the number of generations. For the first 500 generations, it decreases linearly from 0.8-0.5, then it takes a constant value of 0.5.

Mutation operator:
The mutation operator is very effective in finding the absolute optimum point without getting stock in local optimum points.
According to this problem, uniform mutation operator is considered up to the forth gen, but after the forth gen, it should satisfy the following criteria: rand is a random number between [0,1], p m is the mutation rate. Since the Young's modulus has a broad range of values, p m is chosen dynamically and variable with the number of generations. In order to cover all possible answers, p m gets a bigger value for the initial generations. To have convergence of the middle generations, p m gets a lower value. To skip local optimum points and to find an answer, p m increases again at the last generations.
Penalty function: With all the efforts made in order to find suitable chromosome, sometimes at different stages of mutation or crossover operations a wrong chromosome would appear as answer. To avoid this, a penalty function is defined to check all geometric or physical conditions of the chromosome. The chromosomes that pass all the checks are chosen as suitable answers.
Hybrid function: Convergence of GA is not always guarantied, but it could be achieved by having suitable conditions [13] . In complicated problems with many variables, it takes a long time to converge to the right answer. To converge faster, hybrid functions are used.
After a suitable answer is found for the initial estimation of unknown parameters, the GA stops and a faster converging method i.e. CGM takes over until convergence is achieved.

Conjugate gradient method:
CGM is based on minimizing the fitness function i.e. Eq. 13. To minimize this function we need the direction of descent kp (x) and the search step size k β in the following equation [14] : The direction of descent is given by: where, kĝ (x) are gradient directions at iteration k. According to Polack-Ribiere [14] the conjugate coefficients γ k are: where, γ 0 = 0. In cases where algorithm doesn't converge after n step, γ gets the value of zero.
The gradient function kĝ (x) is computed using finite difference method.  [8] , a special regularization function is introduced: In this equation E 0 , ν 0 and L 0 are the first guesses of E, ν and L (circumference of the inclusion). λ is a weighting factor which is computed by trail and error. At the first iteration, it is computed with regard to the minimum of normalized fitness function. The regularization function is added to the function f (x) as: The convergence criterion is: where, ε 1 and ε 2 are chosen such that we have a stable solution to the problem. Details of implementing GA and CGM could be found in [15] .

RESULTS AND DISCUSSION
The obtained results are presented in two parts. In the first part, the validation of the numerical code for GA is checked and in the second part, a number of example problems are solved and some questions related to this inverse problem are investigated.
Analysis model: An inclusion of unknown physical properties i.e. the Young's modulus and Poisson's ratio and also unknown shape is located some where inside a rectangular solid body. To simulate the experimental measurements of displacement, the body Ω 1 with E s = 210 Gpa and ν = 0.3 containing an Aluminum (E a = 70 Gpa, ν = 0.34) or Tungsten (E t = 380 Gpa, ν = 0.2) inclusion with known shape, under uniaxial tension as shown in Fig. 3 is solved by BEM. The dimensions of rectangular outer boundary a and b are one meter and the boundary condition shown in Fig. 3 are hold throughout the investigation. The outer boundary is divided into 64 linear elements and the interface boundary is divided into 36 elements. The boundary displacement on the left and right sides at this model are computed and used as experimental results.
To investigate convergence process, an error coefficient is defined as below: In this error coefficient, x act is the value of real unknown parameter and x est is the final estimated value of the parameter in inverse algorithm. NUP is the number of unknown parameters.
Physical properties of the inclusion are chosen within the following range: 1 GPa E 400 GPa

Fig. 3: Plane inhomogeneous body with boundary conditions
Investigating the validation of numerical code: To validate the proposed numerical scheme, the first problem is to estimate the shape and location of a circular inclusion i.e. (x-0.5) 2 +(y-0.5) 2 = 0.2 2 . The x-y coordinates of the center of inclusion and its radius are estimated. The outer and inner boundaries are divided into 16 and 12 linear elements, respectively. The initial population is 50, crossover fraction is p c = 0.8, the mutation rate is p m = 0.01 and the selection operator is Roulette Wheel. According to Eq. 27, up to 97% of cases converged to the correct solution.
The second problem investigated is to estimate the location and shape of an elliptic inclusion. The same BEM model as previous example is used. The equation of ellipse is 9x 2 +36y 2 = 1 with its center at x c = y c = 0.5. Fourteen parameters, i.e. 12 radiuses and the coordinate of the center of ellipse are estimated. The results are shown in Fig. 4. According to Eq. 27, after 600 generations, up to 90% of the cases converged to the correct solution.
When the two material properties i.e. Young's modulus and Poisson's ratio are added to the geometric parameters, it is observed that the estimation process gets rough and the ability of GA diminishes considerably.
Inclusion identification: After the validity of the proposed numerical scheme is checked, the next problem investigated is to estimate all parameters i.e. geometric and material properties of a pear shaped inclusion, represented by: The GA gives the best estimate of the initial guess of unknown parameters. Then the CGM is employed to get better results. The first step is to find regularization factor λ by trail and error. λ is found and is shown in Fig. 6. The best value of λ for Aluminum is found to be 0.01 and for Tungsten is found to be 0.05. The convergence criteria for this problem are according to Eq. 26. ε 1 and ε 2 are set equal to 10 −11 and 10 −7 respectively. These numbers are chosen with regard to the rang of values of computed outputs (displacements). The results are shown in Fig. 7 and 8. It is observed that error coefficient for this example using GA and CGM together, is about 10%.

The effect of inclusion size on converges:
The first question investigated is the influence of the inclusion size on the estimation process. The investigation is done on an elliptic shape inclusion with equation 2 2 A (x 0.5) B(y 0.5) 1 − + − = for Tungsten. As the percent area of inclusion decrease, the error coefficient in converged solution increases. The given results as shown in Table 1 are the mean value from 10 different runs.
The effect of experimental measurement errors on the estimation process: The next question investigated is the effect of inevitable experimental errors on the estimation problem. Random errors are added to the computed displacements and these are taken to be the measured data. The statistical assumptions regarding the introduced errors are that they are additive, noncorrelated, normally distributed and have zero mean and constant variance. These errors are generated according to: where, σ is the error. At first the GA uses all measured displacements with experimental errors to estimate the inclusion's shape, then the CGM is used in two stages. In the first stage 10 displacements measured on left and right sides  of the boundary as shown in Fig. 5. With symbol are used and in the second stage all the boundary displacements are used. This is done because the experimental measurements have errors and if all the measurements are used, the estimation process diverges after one or two iterations. Therefore after some iteration the changes in the value of unknown parameters become more stable and then the second stage starts where all the displacements are used and convergence is achieved. Since CGM depends on the gradient of the function but GA is only a direct search method, the experimental errors have more effect on CGM than GA. The unknown parameters are estimated by experimental measurements with 1 and 3% errors. The results are shown in Fig. 9 and 10.
For errors up to 3% we get acceptable shape of inclusion but after 3%, the estimated shape would be completely wrong. It is interesting to observe that as the number of erroneous experimental measurement increases, the error in estimated shape also increases but the error in estimated material property decreases.

CONCLUSION
The inverse application of BEM to estimate material properties and the shape of an inclusion using GA and CGM proves to be an effective optimizing • Estimation doesn't need initial guess of unknown parameters.
• Simultaneous estimation of the geometric and physical properties of an inclusion is possible.
• Convergence is achieved even with inevitable measurement errors.
• If the number of erroneous experimental measurements increases, the error in estimated shape also increases, but the error in estimated material property decreases.