Benoît Puel

Dominique Lesselier

Sylvain Chatillon

Pierre Calmon
email: pierre.calmon@cea.fr

Optimization of ultrasonic arrays design and setting using a Differential Evolution algorithm

Keywords: ultrasonics, array, optimization, simulation, transducer

The optimization of phased-array design and setting can be a delicate task when done manually through parametric studies. An optimization method based on an Evolutionary Algorithm and numerical simulation is proposed and evaluated. The Randomized Adaptive Differential Evolution has been adapted to meet the specificities of the addressed applications and in particular to solve multi-objective problems with the implementation of the concept of Pareto's optimal set of solutions. The algorithm has been implemented and connected to the ultrasonic simulation modules of the CIVA software used as forward model. The efficiency of the method is illustrated on two realistic cases of application: The optimization, considered as a mono-objective problem, of the position and delay laws of a flexible array inspecting a nozzle and the optimization of the design of a surrounded array and its delay laws, considered as a constrained bi-objective problem.

Introduction

Over the last decade, the use of ultrasonic arrays for Non Destructive Evaluation has continuously increased and phasedarray techniques play now a major role in the field. Phasedarrays offer numerous advantages in terms of coverage, sensitivity and imaging capabilities. They allow to overpass some limitations of conventional ultrasonic inspections and to tackle increasingly complex inspections. The price to pay for these enhanced capabilities is the complexity increase of the design of the array (number and arrangement of elements) and the definition of the inspection setup (T-R functionalities, electronic delays, transducer displacements, etc.). Nowadays, efficient ultrasonic models are available which can help the expert in these tasks. Up to now this is done through computational parametric studies, which can be time-consuming, fastidious, and somehow delicate. Moreover, due to the necessarily small number of tested set of parameters such manual optimization may lead to one local optimum and miss a better global solution. The motivation of the present work is to propose, prototype and evaluate an automated optimization method based on the use of ultrasonic simulation tools connected to an Evolutionary Algorithm.

Evolutionary Algorithms (EAs) have been used by several authors to solve array optimization problems [START_REF] Rocca | Evolutionary optimization as applied to inverse scattering problems[END_REF]. Chen et al. [START_REF] Chen | An evolution strategy for improving the design of phased array transducers[END_REF] used an Evolution Strategy to design a transducer by looking for the emission surface shape that minimizes the number of rings of an annular phased-array. This algorithm has been used to design an array in order to inspect titanium bullets [START_REF] Lupien | Improved titanium billet inspection sensitivity through optimized phased array design, part I: Design technique, modeling and simulation[END_REF].

Yang et al. [START_REF] Yang | A novel method to design sparse linear arrays for ultrasonic phased array[END_REF] used a Genetic Algorithm to design a sparse array transducer. Gueudre et al. [START_REF] Gueudre | Direct model optimisation for data inversion. Application to ultrasonic characterisation of heterogeneous welds[END_REF] used a Genetic Algorithm to characterize heterogeneous welds based on ultrasonic data inversion. EAs have also been applied in many other domains, e.g., optimization of focus pattern delay laws for ultrasound surgery using a Genetic Algorithm [START_REF] Lu | Design and experiment of 256-element ultrasound phased array for noninvasive focused ultrasound surgery[END_REF] or the retrieval of spheres buried in subsoil using a Differential Evolution [START_REF] Bréard | Hybrid differential evolution and retrieval of buried spheres in subsoil, Geoscience and Remote Sensing Letters[END_REF]. In this work, among many available EAs, which makes non-trivial a choice, the so-called Randomized Adaptive Differential Evolution (RADE) [START_REF] Nobakhti | A simple self-adaptive differential evolution algorithm with application on the alstom gasifier[END_REF] has been chosen in view of its good performance in solving engineering optimization problems and its small number of tuning parameters. RADE has been extended in order to solve multi-objective problems, as well as to handle constraints and criteria of accuracy on the variables. These algorithms, described in section 2, have been implemented and connected to the ultrasonic module of CIVA [9,[START_REF] Mahaut | Recent advances and current trends of ultrasonic modelling in CIVA[END_REF]. This module allows to compute beams radiated by arrays in different materials and to simulate signals arising from defects in complex geometries [START_REF] Benoist | Simulation and applications of NDT ultrasonic array techniques[END_REF]. Through this connection, it has been possible to evaluate the method on realistic NDE cases. In sections 3 and 4 we give two examples of such cases, illustrating respectively mono-objective and multi-objective problems. The first one concerns the optimization of the positioning of a flexible array on a complex part (a nozzle) and of the applied electronic delays. The second one (constrained bi-objective problem) concerns the design and set up of an array for the inner inspection of a pipe. A brief conclusion follows in section 5.

Optimization algorithm

The algorithms that have been used to optimize ultrasonic inspections are described in this section. Since RADE is able to solve only unconstrained mono-objective problems, extensions are proposed to handle constraints, accuracies, and multiobjective strategies.

Optimization problem

Whenever an engineering problem is to be solved by an EA, it should be transformed into an optimization problem that is given, in most general form, as

Minimize f l (X) l = 1, . . . , N o Subject to g m (X) ≤ 0 m = 1, . . . , N c x L j ≤ x j ≤ x U j j = 1, . . . , N (1)
A candidate solution X = [x 1 ; . . . ;

x N ] is a vector of N variables tested to evaluate how well it solves the problem. The problem is defined with N o objective functions f l that have to be representative of the problem at hand; and N c constraint functions g m . The lower and upper bounds are respectively given by the vectors X L = [x L 1 ; . . . ; x L N ] and X U = [x U 1 ; . . . ; x U N ] and define the search range of all variables.

RADE

This algorithm is derived from differential evolution [START_REF] Price | Differential Evolution: A Practical Approach To Global Optimization[END_REF] (DE) with auto-regulation of one tuning parameter (the mutation factor). It makes a population of candidate solutions to converge to the optimum using three evolutionary operators: mutation, crossover, and selection. One loop of those operators is called a generation. The ith candidate solution of the kth generation is noted X (k) i and its jth variable x (k) i j . The strategy used is denoted as DE/rand/1/bin [START_REF] Price | Differential Evolution: A Practical Approach To Global Optimization[END_REF].

The first step is the initialization of the population by a uniform distribution in the search domain defined by

x (0) i j = x L j + rand(x U j -x L j ) (2)
where rand is a random real between 0 and 1. Then, evolutionary operators are used to create new candidate solutions. Mutation creates new vectors (V (k+1

) i = [v (k+1) i1 ; . . . ; v (k+1) iN ]
) by adding perturbations to a reference solution:

V (k+1) i = X (k) a + F i (X (k) b -X (k) c ) (3)
where a, b and c are randomly chosen in the parent population such that a b c i, and where F i is the mutation factor, a self-adapted tuning parameter linked to the ith candidate solution of the population (X i ), as follows. F i is initialized with a random value between 0.1 and 0.9; then, each five generations, F i is reinitialized as well for half of the solutions X (k) i with the lowest evolution of their objective function Γ (see [START_REF] Nobakhti | A simple self-adaptive differential evolution algorithm with application on the alstom gasifier[END_REF] for more details):

Γ(X (k) i ) = f (X (k-5) i ) -f (X (k) i ) (4)
Crossover yields vectors

U (k+1) i = [u (k+1) i1 ; . . . ; u (k+1)
iN ] by combining variables between the new and parent solutions:

x L j x U j u (k+1) i j u (k+1) i j 1 Figure 1: If a new variable u (k+1) i j
is generated out of the search range, it is put back inside by axis symmetry on the closest boundary.

u (k+1) i j =        x (k) i j if rand > C r or j rand i = j v (k+1) i j otherwise ( 5
)
where j rand i is randomly chosen in interval [1, N] for each solution i and C r ∈ [0.1; 0.7] is a user-defined parameter called crossover rate. Selection keeps solutions that improve in the objective function sense:

X (k+1) i = U (k+1) i if f (U (k+1) i ) < f (X (k) i ) X (k) i otherwise ( 6
)
This evolution is performed until stopping criteria are satisfied. Many [START_REF] Zielinski | Examination of stopping criteria for differential evolution based on a power allocation problem[END_REF] can be introduced but the most common one is a maximum number of generations (k max ), which will be the one chosen here.

Constraint handling

Two kinds of constraints have to be handled. First, since evolutionary operators provide diversity, some U (k+1) i can be out of the search range defined by X L and X U . So, variables u (k+1) i j are reintroduced by symmetry on the closest boundary (figure 1):

u (k+1) i j =              2x U j -u (k+1) i j if u (k+1) i j > x U j 2x L j -u (k+1) i j if u (k+1) i j < x L j u (k+1) i j otherwise (7)
The other constraints are inequalities on some variables. The method used here is similar to the one in [START_REF] Huang | Self-adaptive differential evolution algorithm for constrained real-parameter optimization[END_REF]. It is based on the notion of feasible solutions that respect all inequalities and non-feasible solutions that do not respect at least one. This information is given for a candidate solution X at the generation k by the constraint violation function (g (k) (X)), defined by

g (k) (X) = n c m=1 max(0, g m (X)) g (k) m,Max (8)
where g (k) m,Max is the maximum value observed on the constraint g m until generation k. Afterwards, to generate the least modifications in the algorithm, an equivalent objective function (η (k) (X)) is defined by

η (k) (X) =                1 + g(X) n c if g (k) (X) > 0 f (X) -f (k) Min f (k) Max -f (k) Min otherwise (9)
where f (k) Min and f (k) Max are, respectively, the minimum and the maximum values observed on the objective function f until generation k.

This equivalent objective has interesting properties

: if 0 ≤ η (k) (X) ≤ 1 then X is feasible; if 1 < η (k) (X) ≤ 2 then X is non-feasible.
Thus, it can replace the objective function for the selection (equation 6) while it respects Deb's proposition [START_REF] Deb | An efficient constraint handling method for genetic algorithms[END_REF], as well as for the calculation (equation 4) of F i since it does not depend on the order of magnitude due to normalization.

Accuracy handling

RADE had been designed to solve real-variable problems. Since it is interesting to take into account a precision step that determines the accuracy beyond which the solution is of no use to the end-user, a modification (accuracy handling) has been made as described now. The first task is to sample the search space with using the precision step. Then, whenever a new candidate solution (U (k+1) i ) is created, the algorithm finds the closest solution on the chosen grid (called formatted solution) given by

u (k+1) i j = x L j + x P j floor         u (k+1) i j -x L j x P j         (10)
where U (k+1) i is the formatted solution of U (k+1) i , x P j the precision of jth variable, and floor(a) computes the integer part of a.

Yet, one still has to care for diversity preservation since loss of diversity leads to being unable to generate useful new solutions. That is, the mutation operator (equation 3) cannot add a perturbation if X (k) b and X (k) c are too close or equal. So, if the formatted solution has already been computed, one is simply looking for the closest, non-computed solution on the grid within the neighborhood. In brief, computation time is spared since solutions too close to one another are not computed; and the optimal solution found makes sense, from the design point of view, since it has no unusable digits.

Multi-objective

The main difficulty faced with multi-objective problems is that two solutions cannot be easily compared, e.g., since one objective might be better reached and another one less or viceversa. Then, so-called non-dominated sorting is often employed, as is implemented herein. The start point is that a solution

X 1 dominates X 2 if F (X 1 ) F (X 2 ) [16, chap. 8] where F (X) = [ f 1 (X), . . . , f n o (X)], this being defined as F (X 1 ) F (X 2 ) ⇔ ∀i ∈ [1, n o ], f i (X 1 ) ≤ f i (X 2 ) and F (X 1 ) F (X 2 ) (11)
Obviously, solutions can be sorted out in order to produce sets of solutions with an attached rank. Rank one is defined by S 1 , which is the set of all solutions that are not dominated: the Pareto optimal set. Rank two is defined by S 2 , which is the set of all solutions only dominated by some solutions from S 1 . And so on for the next ranks. Obviously, S 1 yields the

f 2 (X) S 1 f 1 (X) S 5 S 4 S 3 S 2 1
Figure 2: The Pareto optimal solutions (red circles) are all non-dominated solutions of all computed feasible candidate solutions (S 1 ∪ S 2 ∪ S 3 ∪ S 4 ∪ S 5 ) in objective functions frame ( f 1 and f 2 ).

best compromises between all objectives and is our main target. Figure 2 shows an example with five ranks of candidate solutions. The rank of a candidate solution X is noted as χ(X). This non-dominated sorting is performed on a database containing all feasible candidate solutions that have been evaluated since the initialization.

Since candidate solutions can be ordered with nondominated sorting, solutions from the same rank should be ordered so as to favor the most isolated ones, a good spread of solutions being aimed at in the Pareto optimal set. To achieve it, a diversity preservation method is applied, involving the normalized crowding-distance, as inspired from the crowding-distance method [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-II[END_REF], properly modified to avoid infinite values. That is, the density of solutions surrounding X in its rank set is estimated via the average distance of two solutions on either side of X for each objective, and is denoted δ (k) (X). The modification considers the distance between the first and last solutions of the rank to estimate the diversity of those solutions.

Then the selection is designed for multi-objective problems to first favor the solution with the best rank, then with the highest normalized crowding-distance [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-II[END_REF]. An equivalent objective function (η(X)) for constrained multi-objective optimization that makes sense to perform the self-adaptation of F i , is then defined as

η (k) (X) =              χ (k) Max + g(X) n c if g (k) (X) > 0 χ (k) (X) - δ (k) (X) n o otherwise ( 12
)
where χ (k) Max is the total number of ranks at generation k.

One example of mono-objective problem

Problem description

This first example given here concerns the inspection of a nozzle with a 8×8 flexible array operating at 2 MHz. The sought defects are breaking cracks on the inner surface at the junction of the two pipes and oriented radially relatively to the secondary pipe as depicted in figure 3 detected thanks to the corner echo involving the reflexion on the inner surface of the nozzle. The objective of the optimization is to find both the position of the array and the focusing point (inputted in the delays computation) that maximize the amplitude of this echo for the different possible angular positions Θ of the defect. The problem is not trivial since this amplitude both depends on the deviation of the beam and on the angle of incidence on the defect. A manual optimization study had been previously performed by an expert for different possible angular positions of the defect. Solutions given by the presently proposed method can be compared to the solution obtained by the expert. Results presented here are obtained for a defect oriented along the primary pipe axis x (Θ = 0 • ).

Variables and objective function

The variables of the problem are the position and orientation of the transducer and of the focal point witch is the input of the delay laws computation. They are listed in Table 1 with their search intervals and precision below which the information is considered useless.

Each candidate solution X generated by the algorithm is evaluated through one objective function f defined as where A max (X) is the amplitude of the corner echo computed by the defect located at Θ = 0 • being assumed to be planar with a height H = 5 mm.

f (X) = -A max (X) (13

Results and discussion

For all the defect orientations (Θ), which have been investigated the solutions given by the automated method are, at least, as good as the solutions manually obtained by the expert in terms of signal magnitude. In the particular case Θ = 0 • , the amplitude of the corner echo has been improved of +5 dB. This improvement is illustrated in figure 4 by the two corresponding simulated Bscan. The parameters of the two solutions are listed in Table 2. The following values of tuning parameters have been used to perform the optimization: Cr = 0.7, NP = 30 and k max = 100. In this example then 3030 candidate solutions have been generated and evaluated by the algorithm. The mean computation time required for one simulation being 42 seconds, the optimization process asked for about 36h total running time, on a standard PC. The evolution of the amplitude corresponding to the solution selected at each generation is reported in figure 5 for two versions of the algorithm: with and without accuracy handling. The amplitudes displayed are given in decibels, the 0 dB level being the manual solution. It can be seen that the accuracy handling significantly reduces the number of generations and therefore the corresponding computation time required to converge to the optimized solution: 33 generations and 12 hours vs. 50 generations et 18 hours.

One example of multi-objective problem

Problem description

The second example given in this contribution concerns the inner inspection of a stainless steel pipe by a surrounded array, depicted on figure 6. It is solved as a constrained multiobjective problem. The inner diameter of the pipe is 50 mm and its thickness 5 mm. The aim of the control is to detect cracks oriented along the pipe axis with 45 • compressional waves. The array is moved along the axis z while electronic commutation is used to inspect the whole section of the pipe. The pipe is filled by water. The goal is to optimize the design of the array, and the electronic setting: the electronic commutation and applied delay laws. The number of available channels, and therefore the maximal number of elements of the array, is fixed at 256. The objective is to guarantee the good characteristics of the beam, direction (45 • on the defect) and focalization on the outer surface, and to minimize grating-lobes.

Variables

The optimization concerns the design of the array (dimensions, arrangement and central frequency of the elements) the definition of the commuting aperture and the angular deviation defining the electronic delays. The variables, search intervals and accuracy are listed in Table 3.

Objective functions

The definition and evaluation of the objective functions are based on the computation by CIVA of the ultrasonic beams radiated by the array in the pipe and corresponding to the different candidate solutions. More precisely, only the longitudinal component is calculated and the reflection on the outer surface is not taken into account. Computation parameters are sketched in figure 7. The computation zone is a rectangular grid of 161 × 61 points, distant of 0.2 mm along both axis. The steel is assumed to be isotropic with compressional wave velocity of 5650 m/s and negligible attenuation. The sound wave velocity in water is set at 1483 m/s in water. For each deviation angle θ and active aperture, the delay law applied to the array is computed by using the geometrical approach available in the CIVA package. Two objective functions are defined: ( f 1 ) ensures a beam orientation as close as possible to 45 • on the defect, and ( f 2 ) maximizes the ratio between useful amplitude distribution inside the beam and unwished amplitude distribution outside the beam. The later mainly corresponds to the possible apparition of the grating lobe, which is to be avoided.

Beam orientation

The objective function f 1 is defined as the difference between 45 • and the average orientation in the beam. This average orientation is evaluated on a set of pixel A i of highest amplitudes.

f 1 (X) = 1 N N i=1 |45 • -α i (X)| (14)
The orientation (angle α i ) assigned to one pixel A i is issued from the local direction of propagation d i output of the CIVA computation. α i measures the angle between d i and the normal of the outer surface at B i (see figure 8).

Amplitude distribution

The spatial extension of the useful beam is defined for each candidate solution from the -6 dB width ∆ and the direction of propagation d at the focal point F (see figure 9), which is an output of the field computation. The second objective function f 2 , related to the amplitude distribution, is expressed as

f 2 (X) = - Amp in (X) Amp out (X) (15)
where Amp out (X) is the average amplitude of all the pixels laying outside the useful beam:

Amp out (X) = Amp out i out [START_REF] Angell | Optimization Methods in Electromagnetic Radiation[END_REF] and Amp in (X) is the weighted average amplitude of all the pixels laying inside the useful beam:

Amp in (X) = w i Amp in i in (17)
The weights w i are introduced in order to favor solutions with highest amplitude on the outer surface. w i decreases with increasing distance to the outer surface l i .

Constraints

One constraint function is introduced to guarantee that the number of elements of the array is lower or equal to 256. The expression of the constraint is:

g 1 (X) = 40π g(X) + w(X) -256 ≤ 0 ( 18

Results and discussion

The following values of tuning parameters have been used to perform the optimization: Cr = 0.7, NP = 25, and k max = 50. 1530 candidate solutions have been generated among which 1422 were feasible. The evaluation of feasible solutions (implying a field computation) took 50h total running time on a standard PC. The Pareto optimal sets obtained at generations 0, 10, 20 and 50 are depicted in figure 10. At the 50th generation, the Pareto set is composed of 94 solutions corresponding to different compromises between f 1 and f 2 . Along the front, f 1 and f 2 range respectively from 1 • to 25 • and from 11 dB to 19 dB. All the solutions have approximately the same values of array length (l = 4 ± 0.5 mm), element width (w = 0.2 mm) and gap between elements (g = 0.3 mm). Unsurprisingly, these values correspond to the allowed maximum number of elements (256) and minimum width of the elements. The three parameters that mostly influence the compromise between the two objectives are the number of elements in the active aperture n, the frequency f , and the deviation angle θ. We give three representative solutions of the Pareto optimal front corresponding to reasonable compromises between the two objectives in the Table 4. The computed fields associated to these three solutions can be seen in figure 11. Solution 1 favors the beam orientation while solution 3 maximizes the focalization on the outer surface, the solution 2 being an intermediate solution. This example illustrates the interest of the algorithm which proposes in such multi-objective problem the best set of compromises among with the expert can select the final solution.

Conclusion

A method based on the use of a Differential Evolution algorithm, the Randomized Adaptive Differential Evolution (RADE) has been proposed to optimize array design and setting. The RADE has been adapted to meet the specificities of the addressed NDT applications and in particular to: i) manage the constraints on the variables in the problem: bounding domain and inequality constraints, ii) consider the finite precision of the variables, iii) solve multi-objective problems with the implementation of the concept of Pareto's optimal set of solutions. The algorithm has been implemented and connected to the ultrasonic simulation modules of the CIVA software platform. For one given solution (expressed as a set of parameters input into the simulation) the value of the objective function is inferred from the result of the corresponding simulation. The method has been evaluated on two realistic cases of application considered respectively as one mono-objective problem and one constrained multi-objective problem: the optimization of the 4.

positioning and electronic delays applied to a flexible array inspecting a nozzle and the optimization of the array design and electronic delays for a pipe inspection. The results confirmed the ability of the algorithm to converge towards the desired solution within a reasonable time. The comparison with manual reference solutions when available highlighted the equivalent or better quality of the solution provided by the algorithm.

1 Figure 3 :

 13 Figure 3: Schematic view of the nozzle inspection. The array probe is located on the outer surface of the secondary pipe. The sought defect is schematized at the junction of the two pipes. The focal point is materialized by the green dot.

1 Figure 4 :

 14 Figure 4: B-Scans comparison between solutions obtained by the expert and the RADE. Y Θ α r θ ϕ X exp 91 49 20 153 44 0 X RADE 94.5 54.6 22.5 188 39.4 0.4

Generations 1 Figure 5 : 1 Figure 6 :

 1516 Figure 5: Evolution of the objective function (the maximum of the echo amplitude) during the optimization process with (in blue) and without (in red) accuracy handling. The manual solution corresponds to the 0 dB level.

3 : 1 Figure 7 :

 317 Figure 7: Sketch of the computation zone and array setting. The electronic delays applied to the rotating active aperture are plotted in red.

1 Figure 8 : 1 Figure 9 :

 1819 Figure 8: Definition of the parameters related to the objective function f 1 , sketched on the map of the computed ultrasonic field corresponding to one candidate solution.

Figure 10 : 1 Figure 11 :

 10111 Figure 10: Pareto optimal sets obtained each ten generations during the optimization of the surrounded array, where f dB 2 (X) = 20log(-f 2 (X)).

Table 1 :

 1 Variables defined to optimize position and setting of a flexible array for a nozzle inspection.

	name description	interval	accuracy	unit
	Position of the array on the nozzle
	Y	vertical position [30; 100]	0.5	mm
	Θ	angular position	[0; 90]	0.3	degree
	α	probe rotation	[0; 90]	0.5	degree
	Focal point coordinates relative to the probe center
	r	distance	[100; 250]	1	mm
	θ	polar angle	[-80; 80]	0.3	degree
	ϕ	azimuthal angle [-80; 80]	0.3	degree

Table 2 :

 2 Parameters of the solutions obtained by the expert (X exp ) and the RADE (X RADE ).

Table 4 :

 4 Three interesting solutions from the Pareto front.

	)

Acknowledgement

This work has been supported by the Région Ile-de-France in the framework of the CAPVERS project.