Abstract: Simulated annealing algorithms are evolutionary processes that may well be applied to solve optimization problems under some hypotheses such as non-convexity and non-differentiability, etc. However, premature convergence is an intrinsic trait of such simulated annealing that makes them incompetent of searching best solution of the problem domain. In order to evade this problem, a new simulated annealing method for solving inverse electromagnetic problem is proposed. Adaptive simulated annealing is used as the new method in electrical engineering for solving the parameters identification problem. With the target as the metal under test, the calculation of the eddy current densities produced between the coil and the material is also detailed. Here the aim of the research explains the calculation procedure for the reactance and resistance of the eddy current sensor exited by voltage supply using a fast semi analytical method called coupled circuits method. Both the Coupled Electromagnetic Circuits Method (CECM) and the new approach (ASA) are used to estimate the electric conductivity and magnetic permeability of a circular material under test. To have the quality of this method, the performances of ASA are compared with other algorithm such as Genetic Algorithm (GA) in term of accuracy of the solution and computation time.

INTRODUCTION

The Simulated Annealing Algorithm (SA) (Yang et al., 2000) and Genetic Algorithm (GA) (Chen et al., 2001) have been extensively used to solve the optimization of electromagnetic problems. Among the principals advantages of these methods are the ability to avoid many difficulties such as the calculated gradient and local minima solution. However, it should be expected that GA may be better suited for some problems than SA (Glover and Kochenberger, 2003).

Also, in optimal electromagnetic field where the functions are frequently nonlinear and multimodal, the simulated annealing cannot afford the adequate fidelity of the electromagnetic inverse problem because the convergence to an optimal solution cannot theoretically be guaranteed after a number of iterations and function evaluations.

Generally, simple simulated annealing always takes much computation time in finding a globally optimal solution which is not acceptable for many engineering applications (Glover and Kochenberger, 2003). To improve the convergence of SA to some area, many supplementary methods have been proposed (Renyuan et al., 1996). SA using Gaussian distribution had been shown that the convergence is fairly slow (Drago et al., 1999). Tasllis and Stariolo (1996) planned a rapid SA using Cauchy algorithm. If the field computation using finite element modelling of an optimization problem takes several minutes, it is difficult to apply these methods to it for their long running time (Glover and Kochenberger, 2003). To overcome these difficulties for the best minimum, adaptive simulated annealing has been developed for solving the inverse electromagnetic problem. Recently, the adaptive simulated annealing algorithm is widely used because they are efficient and have both the ability to find large-scale solution and evade premature convergence. The algorithm is planned so that facilitates a global search and escapes the local minima. This improved algorithm can be worked adequately when the cost function is multimodal and under some hypotheses such as non-convexity and non-differentiability.

This study deals with the identification of material parameters using an inverse strategy. The method based on the use of coupled electromagnetic circuits method (Mohdeb et al., 2009a, b; Cheriguen and Mekideche, 2003) and adaptive simulated annealing algorithm. The identified material parameters are the relative permeability and electric conductivity. In order to quantify the quality of the agreement between the measured and calculated responses, an objective or cost function has to be defined.

MATERIALS AND METHODS

Sensor description and equations
Description of the system: The test configuration chosen for the evaluation of the new optimization method is shown in Fig. 1. The simplified 2D configuration is a circular conductive plate placed underneath a flat spiral coil where the optimization target is to identify the electric conductivity and magnetic permeability from the coil impedance measurement. The flat spiral coil is constituted of a bobbin with Nw wires (in our problem two wires). The bobbin is supplied by a sinusoidal voltage source with constant amplitude U and pulsation w. The forward problem predicts the coil impedance calculation wite more excitation frequency using the coupled electromagnetic circuit (Mohdeb et al., 2009a, b; Cheriguen and Mekideche, 2003).

Electromagnetic field computation: The axsiymmetrical model is based on the Coupled Electromagnetic Circuits Method (CECM) which permits to calculate the coil impedance. The CECM analysis is used to calculate the coil impedance and the eddy currents from the magnetic potential vector. Generally, this method is suitable for some hypothesis:

In this research, researchers exploitd the CECM for solving materials properties determination inverse problem because this model is fast-running. The CECM consists in associating the integral form of the solution to a subdivision in elementary coils. In CECM only the conductive regions are meshed. The current densities are the unknown vector in the inductive coil and the load. The materials are discretized in elementary loops for determine these unknowns. The bobbin and load are represented by Nb and Np elementary coaxial loops, respectively (Nb; Bobbin and Np; Plate). Each elementary loop is in magnetic interaction with itself and with the other ones. Each loop is discretized in several elements (Ne discretizations). So, the total number of unknowns is Nb + Np. The CECM mesh generated for the regions conductive is shown in Fig. 2b for 2 wires. The relations between 2 loops can be explicated with the electric transformer model (Fig. 2a). The discretization scheme is an extension of this electric transformer model (Cheriguen and Mekideche, 2003; Mohdeb et al., 2009a, b).

Fig. 1:	Example of typical joining application

Fig. 2:	Model of coupled electromagnetic circuits method (a) Equivalent electric circuits, (b) CECM discritization

Maxwell’s equations are a set of equations stating the relationship between the fundamental electromagnetic quantities. We are going to search from Maxwell’s equations and the Ohm’s law an equation that describes the electromagnetic phenomena in an elementary circular loop (Fig. 3). This equation is the basis formulation of the coupled electromagnetic circuits method (CECM).

The magnetic and induction effects between the different domains in axisymmetric devices are represented by the Maxwell’s equations. Then The mathematical model with 2D axisymmetrical case for quasi-static problem can be written from Maxwell’s equations as follow:

Fig. 3:	Representation of an elementary circular loop

Fig. 4:	Distribution of eddy current in load at frequency 3 kHz

(1)

Where:

The gradient of the potential in the circular reference is shown by Gradd = -U/2πr. In the problem the excitation varies sinusoidally with time then ∂/∂t can be written as jω. From both the Maxwell’s Equations and the Ohm’s law and considering the simplified notation of the gradient of the potential in the circular reference, the combination of electromagnetic system and the voltage supply in elementary loop p is shown by:

(2)

To calculate the potential vector magnetic A at point p generated by the current densities J (q) at q point, researchers used Biot Savart’s law. This law is written as follow:

(3)

By simplification of this equation, the sum magnetic vector potential A in a known point p with the current densities, is expressed as:

(4)

Note that:

(5)

Where S is the gross section of the elementary loop p, r and z are cylindrical coordinates. The functions E₁ and E₂ are the Legendre function of the first and second kinds. When applying circuit’s laws the equation above is simplified (after elimination magnetic potential in our system) in a linear system as:

(6)

Where:

Eddy currents calculation: Once the magnetic vector potential has been determined, all electromagnetic field quantities can be calculated. The current densities in the conducting plate are expressed as:

(7)

The total impedance of the exciting coil can be calculated from the voltage supply and the current densities in the loops. In that case, the expression is:

(8)

where, Z_coil is the coil impedance of the eddy current sensor.

Adaptive simulated annealing method: Simulation optimization by simulated annealing was first described

by Glover and Kochenberger (2003) and is based on research by Yang (2000) in the area of statistical mechanics. A simulated annealing optimization starts with a metropolis Monte Carlo simulation for state-space variables at a high temperature.

This means that a relatively large percentage of the random steps that result an increase in energy will be accepted. After a sufficient number of Monte Carlo steps or attempts, the temperature is decreased. The acceptance of the novel result is according to the Metropolis’s condition based on the Boltzman’s probability (Glover and Kochenberger, 2003). SA algorithm contains two steps: the 1st, perform search while the temperature is decreasing. The 2nd determine the acceptance. The acceptance probability of solution point is defined by:

(9)

Where:

If, Ej’-Ej≤0,j’ is accepted as a starting point for the next iteration; otherwise solution j’ is accepted with Boltzman’s probability. The above procedure is repeated time until temperature T is reduced. The aim of the Metropolis’s succession is to authorize the system to attain thermal equilibrium. It should be noted that classical optimization algorithm only accept improved design and never accept a worse design. In simulated annealing, the condition Ej’-Ej≥0 gives the algorithm a chance of get out of a local minimum.

In practice, a geometric cooling schedule is generally utilized to have SA settle down at some solution in a finite amount of time. It has been proved by some reseachers that by carefully controlling the rate of cooling of the temperature, SA can find the global optimum. However, this requires infinite time. Fast annealing and Very Fast Simulated Reannealing (VFSR) or Adaptive Simulated Annealing (ASA) are each in turn exponentially faster and overcome this problem.

The first simulated annealing employed Gaussian distribution as a generator and was proposed by Kirkpatrick. Glover and Kochenberger (2003) proposed a fast simulated annealing by using Cauchy/Lorentzian distribution. Another modification of the SA, the so-called adaptive simulated annealing was proposed by Ingber (1993, 2003) and was designed for optimization problem in a constrained search space. For x^k a parameter in dimension i at annealing time k with rang x^kε(x_i^min, x_i^max) the new value is generated by:

(10)

Where:

(11)

The i values T_i and are identifies the parameter index and temperature. To find λ_i one most find the normalized cumulative probability distribution of g (λ_i). The cumulative probability distribution can be defined as:

(12)

To simplify this generating function λ_i for a uniform distribution is preferred. A normal uniform distribution is defined as follow:

(13)

where, u_i (ε [0, 1]) is the uniform distribution function and is the cumulative of a uniform distribution. Each parameter is generated using a cumulative function. In this case by the idea of Ingber it can be seen to choose g (λ_i) = u_i. Then, to calculate λ_i according to the preceding distribution, researchers applied this formulation:

(14)

The new generation distribution function in ASA has much fatter trails than Gaussian and Cauchy generation function. Temperature T is a key element in the cooling system in the ASA algorithm. After every generated points, annealing takes place with a new annealing schedule. A global optimum can be obtained statistically if the annealing schedule is:

(15)

Where:

Reannealing in ASA algorithm periodically rescales the generating temperature in terms of the sensitivities s_i calculated at the most current minimum values of the cost function. After every acceptance points, reannealing takes place by the first calculating the sensitivities:

(16)

The annealing time is adjusted according to s_i based on the heuristic concept that the generating distribution used in the relatively insensitive dimension should be wider than that of the distribution produced in a dimension more sensitive to change.

Genetic algorithm: The Genetic Algorithm (GA) is an optimization and search technique based on the principles of genetics and natural selection (Haupt, 1995). The method was developed by John Holland over the course of the 1960s and 1970s and finally popularized by one of his students, David Goldberg (Rahmat-Samii and Michielssen, 1999), who was able to solve a difficult problem involving the control of gas-pipeline transmission for his dissertation. Genetic operators are used to explore a solution space through the generation of new individuals from existing ones. There are two major genetic operators: mating and mutation (Glover and Kochenberger, 2003). Mating operators generate new solutions by exchanging the genetic information among two or more individuals from the population. Mutation operators generate new solutions by randomly changing the genetic information of singles individuals. The genetic operators deal with the individuals in a population over several generations to improve their fitness gradually. Individuals standing for possible solutions are often compared to chromosomes and represented by strings of binary numbers. Like the other method, GA is also expected to find the global minimum solution even in the case where the objective function has several extrema including local maxima, saddle points as well as local minima. The details of genetic algorithm are discussed by Chen et al. (2001), Glover and Kochenberger (2003) and Rahmat-Samii and Michielssen (1999).

RESULTS AND DISCUSSION

The geometry of the problems considered is shown schematically in Fig. 3b. Considering the symmetry, the model is only designed for the half.

Table 1:	The results from the coupled circuits method

Fig. 5:	Distribution of eddy curent in load at frequency 10 kHz

The aim of the simulation study is to calculate the coil impedance from the eddy current densities in the all system (bobbin). The parameters used in the computation for harmonic excitation are inner radius of the coil (r₀ = 2.5 mm), outer radius of coil (r₁ = 3.5 mm), coil width (l₀ = 1 mm), the distance between coils is 1 mm (r2), relative permeability (μ_r = 1 for coil and air, 100 for load), conductivity of half space (σ_s = 5.7e7 S m^-1 for coil (copper) and = 7.6e6 S m^-1 for load (steel)), lift off (e = 1 mm). Frequencies used for excitation of the coil are 0.1-30 kHz. The piece to test is a cylinder of 0.80 mm thickness (l₁) and the dimension of it is 14 mm length (l). The coil with two wires has been energized by voltage supply with U = 30V at several frequencies. The model of the probe coil is shown in Fig. 2b.

Table 1 shows the variation of the coil impedance with several frequencies: f1 = 0.1 kHz, f1 = 0.4 kHz, f1 = 1.5 kHz, f1 = 9.0 kHz, f1 = 30.0 kHz. The Fig. 5 and 6 show the distribution of the current densities in the load. It is important to select an adequate mesh to represent correctly the electromagnetic phenomena and then to reduce the numerical errors that can influence the convergence of the identification process.

In this case, every coil of this sensor contains 64 elementary loops distributed in 8 following the axial direction and 8 following the radial direction and researchers have considered 320 loops in the load (40 along the radial direction and 8 in the axial direction). The experiment was run on a Dell which contain an Intel Pentium 4, 3.6 GHz CPU and 256 Mb RAM. The program was implemented in MATLAB 7.1.

Fig. 6:	Methodology of inverse problem procedure

When a voltage supply is used to excite a coil, a magnetic field is produced and magnetic lines of flux are concentrated at the center of the coil. Then as the coil is brought near an electrically conductive material, the magnetic field penetrates the material and generates continuous, circular eddy currents. Larger eddy currents are produced near the test surface. As the penetration of the induced field increases, the eddy currents become weaker. The induced eddy currents produce an opposing (secondary) magnetic field. This opposing magnetic field, coming from the material has a weakening effect on the primary magnetic field and the test coil can sense this change. In effect, the impedance of the test coil is reduced proportionally as eddy currents are increased in the test piece.

The model in voltage source driven gives the possibility to compute the inductive current and therefore the impedance of the system that is interesting for solving material properties determination inverse problem while using as observables the coil impedance of the eddy current sensor.

For identification parameters of the real parameters, a lot reiterate are required for predicting the permeability magnetic and the electric conductivity. In this step, researchers described the implementation of the ASA algorithm in parameters identification approach. Figure 6 shows the scheme of electromagnetic inverse problem. The proposed methodology can be summarized as follows: Choose a true experimental test and used for the identification procedure and saving the measured Z^mes. Generate randomly the input parameters.

The forward problem or the model program is applied to simulate the output vectors Z^cal. Optimization analysis by ASA algorithm simulation. This step is performed by the calculation of the error function for new materials parameters X^sim. Verification of the ASA results with original measured parameters X^mes. The new results Z^cal is calculated by using the model code CECM and compared its values with original measured results Z^mes.

The Fig. 6 shows the identification procedure in the case of non destructive evaluation of a parameter σ_c and μ_r from the measured (or true) impedance Z_mes and the calculated impedance Z_cal. If the calculated values are generally different from than that of measured values the algorithm are used for minimized this difference margin. Researchers performed the estimation approach of the electric conductivity and the relative permeability of the material under test while solving an inverse problem. Here, the inverse problem to analysis, is expressed as follow: To find

(17)

where, Z_cal is the impedance calculate and Z_mes is the impedance measured. In this study, the measured values are replaced by those gotten, while using the direct model (Coupled Electromagnetic Circuits Method CECM).

Researchers defined the error function as the difference between the properties measured and calculated values by the coupled circuits method. Then, the identification is considered as nonlinear problem to find a solution x that minimizes the function, as follows:

(18)

where, R_cal and X_cal denote the resistance and reactance calculated by the CECM approach, R_mes and X_mes are the resistance and reactance measured at the mth point due to a frequency. This function (Eq. 18) is minimized by using the new hybrid CECM-ASA. The values of μ_r and σ_c are predicted through minimization of this cost function. The first experiment involved the use of additive noisy data to examine the proposed method. Two cases are investigated; the same parameters geometric and physic are used. The material properties of cases are respectively set as:

It should be mentioned that these material properties are used only to provide the simulated measurements of impedance responses using CECM code and to check the accuracy of the material properties determined by the inverse analysis (CECM-ASA) using these simulated measurements.

Minimal and maximal initial values of the materials were respectively for σ_c (1e6-15e6) and μ_r (50-200). The bounds on the materials parameters are required to define a finite search space for the ASA. In engineering practice, a narrower range is always preferred for accuracy in inverse solution and for computational efficiency.

Table 2:	The parameters for the new approach

Table 3:	ASA results with different noise levels

The settings of parameters used in this method for the two tests are shown in Table 2. Noise is inevitably involved in the measured data in practice. A few test cases were also run considering noisy data.

The results are shown in Table 3. In order to simulate the measured impedance, a Gaussian noise defined by traditional equation (τ, δ where τ is a random number in the range [-1, 1] and represents the standard deviation of the measurement errors) is directly added to the coil impedance responses and then the noise-contaminated responses are used as inputs for the identification. To investigate the sensitivity and stability of the present inverse procedure to the noise level, two noise levels of 4 and 9% are considered.

The presence of noise can make the identification much more complicated compared to the noise-free case. If the noise is too large, a local optimum could be found as the true results rather than the global optimum. It has been found that when the noise is >9%, the true results were not identified. If the noise is>4%, the true results can still be found as shown in Table 3 and the characterized result remains stable regardless of the presence of the noise.

With ASA optimization, the convergence to an optimal solution can theoretically be guaranteed after a number of iterations. Interestingly, when a combination of adaptive simulated annealing and the semi analytical method was applied and even better result was achieved. This can be explained with the fact that the ASA method has different strength. The adaptive simulated annealing is very good at finding the correct area of the solution, tolerant of local maxima and minima and the new generation function is excellent at refining a solution systematically to the nearest maximum or minimum (best solution).

Table 4:	Simulation results for various search algorithms

The new algorithm is better equipped for global optimization because it is more aggressive in the exploration of the search space. This algorithm can be worked adequately when the cost function is multimodal and not derived for the design parameters. Now the researchers compared this new method with other optimization method such as Genetic Algorithm (GA) that determines the identification parameters in any problem. Table 4 shows the values of parameters of the material under test and the results of optimal computations, using the ASA and GA for the noise-free case.

For this optimization method, the code has been programmed in MATLAB 7.0 with the Genetic algorithm Toolbox (MATLAB 04). The parameters used for the genetic algorithm are selected as follow: The test start with number of population equals to 60 with 20 generations. Each generation stores the best fitness string and at the end gives us the best candidate. A binary encoding is used. The crossover probability is equal 0.61 for the test. The mutation probability was 0.001. Also, the method of tournament selection is used.

As observed for this problem, the ASA method resulted in the optimal solution close to global optimal solution but the GA solution was stuck in local optima. The performance of the ASA was the best among genetic algorithm. As such, not only can the global optima be ensured but results can also be obtained at a reasonably fast speed (Table 4). The other advantages of ASA are the capability to escape from the local optima (Table 3). When compared the CECM-ASA with the CECM-GA, the numerical results show that the adaptive simulated annealing gives us an excellent convergence in a minimal CPU time (Table 4). It is evident from the above results that adaptive simulated annealing is superior to on this problem both in terms of optima found and speed convergence.

CONCLUSION

A new stochastic optimization based on adaptive simulated annealing algorithms was carried out for solving inverse electromagnetic problem. A new Simulated Annealing Algorithm (ASA) is an extension of the traditional simulated annealing algorithm. It uses a new function generation to reduce the likelihood of the premature convergence.

When used to solve the optimization problem in the parameters identification of a eddy current sensor, which its objective function is under some hypotheses such as non-convexity and non-differentiability, adaptive simulated annealing can not only obtain the global optimal solution but also the convergence history showed that the ASA converged to the optima faster than the Genetic Algorithm (GA). The results proved satisfactory in terms of the best identification parameters for eddy current sensor showing a good behaviour of the Adaptive Simulated Annealing Algorithm (ASA) in optimizing the analyzed problem.

Finally, the new ASA algorithm can be extensively used in any other situation to solve different optimization problems of electromagnetic devices.