I Clairand

J-M Bordy

J Daures

J Debroas

M Denozière

L Donadille

M Ginjaume

C Itié

C Koukorava

S Krim

A-L Lebacq

P Martin

L Struelens

M Sans-Mercé

M Tosic

F Vanhavere

ACTIVE PERSONAL DOSEMETERS IN INTERVENTIONAL RADIOLOGY: TESTS IN LABORATORY CONDITIONS AND IN HOSPITALS

Active personal dosemeters in interventional radiology:

Tests in laboratory conditions and in hospitals

INTRODUCTION

The optimization of the use of active personal dosemeters (APDs) in interventional radiology and cardiology (IR/IC) is performed by one of the work packages of the ORAMED project, a Collaborative Project (1) (2008-2011) supported by the European Commission within its 7 th Framework Program. APDs are used for monitoring of occupational exposure in many applications involving ionising radiation, especially the nuclear industry. In hospital environments, they are much less used [START_REF] Ginjaume | Overview of active personal dosemeters for individual monitoring in the European Union[END_REF] . In IR/IC, the possibility to assess the dose and/or dose rate in real time is particularly interesting since operators can receive relatively high doses while standing close to the primary radiation field. In addition, an attractive feature of the APD is the possibility to have an alarm when a particular dose rate or dose value is exceeded. Due to the specificity of the X-ray fields used in IR/IC (low energies and pulsed fields), the current technology of APDs can be inadequate. This problem was highlighted during two previous international intercomparisons [START_REF] Clairand | Intercomparison of active personal dosemeters in interventional radiology[END_REF][START_REF] Bordy | Evaluation of active personal dosemeters for interventional radiology[END_REF][START_REF]Intercomparison of personal dose equivalent measurements by active personal dosimeters[END_REF] . The work presented in this paper consisted in:  studying the real radiation field characteristics encountered in IR/IC in terms of energy, angular distribution, dose rate and pulse characteristics,  making a selection of commercial APDs deemed suitable for application in IR/IC according to several criteria, in particular the capacity to respond to photon energies down to 20 keV,  testing, under laboratory reference conditions, the dose, dose rate, energy and angular response of the selected APDs,  studying, under laboratory reference conditions, the effect of the dose rate, pulse frequency and pulse width on the APD response,  performing tests in several European hospitals in workplace conditions.

MATERIAL AND METHODS

Typical fields in IR/IC

The typical fields and parameters encountered in IR/IC were gathered through questionnaires sent to hospitals, a literature search and quality control outputs. Calculations of dose equivalent rate at specific points of interest and typical scattered spectra were performed using the Monte Carlo codes MCNPX and PENELOPE [START_REF]MCNPX User's manual, version 2.5.0[END_REF][START_REF] Salvat | A code system for Monte Carlo simulation of electron and photon transport[END_REF] .

Selection of APDs

The selection of commercial APD models was based on the results from international intercomparisons [START_REF] Clairand | Intercomparison of active personal dosemeters in interventional radiology[END_REF][START_REF] Bordy | Evaluation of active personal dosemeters for interventional radiology[END_REF][START_REF]Intercomparison of personal dose equivalent measurements by active personal dosimeters[END_REF] , and on their availability in different European countries. A pre-requisite for consideration was that each device should respond to photon energies down to 20 keV. Following those criteria, seven commercial APDs were selected (Figure 1

Tests of APDs under laboratory conditions

The tests with continuous X-ray fields were made in two calibration laboratories (IRSN in France and SCK•CEN in Belgium). These tests were performed to determine the dose, energy, dose rate and angle responses of the above mentioned APDs. Two devices of each type were always used. The following reference were used (N-15, N-20, N-25, N-30, N-40, N-60, N-80, N-100, N-120, S-Cs and S-Co) as defined in the ISO 4037-1 [START_REF]X and gamma reference radiation for calibrating dosemeters and doserate meters and for determining their response as a function of energy --Part 1: Radiation characteristics and production methods[END_REF] standard.

The tests [START_REF] Denozière | Pulsed X-rays for interventional radiology: test on active personal dosemeters (APD) European project FP7 ORAMED WP[END_REF] in pulsed mode were made at the French standard laboratory for ionizing radiation (Laboratoire National Henri Becquerel -LNHB, CEA LIST in France). The influence of several parameters on the response of the APD in pulsed mode was studied (70 kVp, HVL 5.17 mm Al):  the effect of the dose equivalent rate (i.e. the mean dose equivalent rate during one pulse) from 1 to 55 Sv.h -1 for a pulse duration = 20 ms and a pulse frequency = 10 s -1 -tests performed in multi-pulsed mode,  the effect of the pulse frequency from 1 to 20 s -1 for a dose equivalent rate = 1.8 Sv.h -1 and a pulse duration = 20 ms -tests performed in multi-pulsed mode,  the effect of the pulse width from 20 ms to 1000 ms for a dose equivalent rate = 1.8 Sv.h -1 -tests performed in single pulsed mode; for technical reasons, tests under 20 ms were not possible.

Tests of APDs under realistic conditions in hospitals

A series of tests were made in ten European hospitals during routine practice. The interventional radiologists and cardiologists were asked to wear, side by side and above their lead apron, an APD and an additional passive dosemeter during daily practice. The dosemeters were worn during several interventions to integrate doses of at least 300 Sv for several types of IR/IC procedures. The main objective of these tests was to compare the measurements performed by the active and passive dosemeters worn in routine practice in hospitals, where all kinds of procedures and parameter settings are used and without an accurate knowledge of the field parameters.

For practical reasons, only four dosemeters were tested in these realistic conditions: DMC 2000XB, EPD Mk2.3, EDM III and DIS-100.

RESULTS

Typical fields in IR/IC

The compilation of data presented in Table 1 gives an overview of typical fields encountered in IR/IC independently of the procedure considered.

The instantaneous dose equivalent rate obtained by quality control measurements using DAP-meters in the direct field at the level of the table ranges from 2 to 360 Sv.h -1 . The dose equivalent rate in the scattered beam, at the level of the operator for tube position of 0° and 90°, was found to range from 5.10 -3 to around 10 Sv.h -1 . In addition, Monte Carlo calculations showed that the influence of the filtration on the scattered spectra is very small. The energy of the scattered spectra ranges from 20 to 100 keV. Dose equivalent rate in the scattered beam (operator -above the lead apron) 5.10 -3 to 10 Sv.h -1 *Corresponding author: isabelle.clairand@irsn.fr

Response of APDs in continuous mode

In continuous mode, all APDs have a linear response with dose. The energy response (Figure 2) is within the interval [0.71 -1.67] as required in the IEC 61526 standard [START_REF]Radiation protection instrumentation. Measurement of personal dose equivalent Hp(10) and Hp(0.07) for X, gamma, neutron and beta radiation: Direct reading personal dose equivalent and/or dose equivalent rate dosemeters[END_REF] from 60 Co energy down to N-30 for all APDs except EDD30. For EDD30, these results are consistent with the fact that this APD is designed to work specifically at low energy. Figure 3 illustrates the response of the selected APDs as a function of dose equivalent rate. Most tested APDs provide a response for dose equivalent rates up to 10 Sv.h -1 , except PM1621A for which the response is diverging rapidly from 1 Sv.h -1 and EDD30 which saturates for dose equivalent rates above 2 Sv.h -1 . The dose equivalent rate range requirement specified in IEC 61526 [START_REF]Radiation protection instrumentation. Measurement of personal dose equivalent Hp(10) and Hp(0.07) for X, gamma, neutron and beta radiation: Direct reading personal dose equivalent and/or dose equivalent rate dosemeters[END_REF] is 1 Sv.h -1 thus, all the dosemeters fulfil the IEC standard. All devices showed under and over-responses for low energy photons and high angles, but these stayed within the limits of the IEC standard [START_REF]Radiation protection instrumentation. Measurement of personal dose equivalent Hp(10) and Hp(0.07) for X, gamma, neutron and beta radiation: Direct reading personal dose equivalent and/or dose equivalent rate dosemeters[END_REF] except AT3509C (Figure 4) for which the angle response is inside this interval at 60° only from N-80.

Response of APDs in pulsed mode

Effect of dose equivalent rate

For most APDs, the response decreases when the dose equivalent rate increases (Figure 5). For dose equivalent rates lower than 2 Sv.h -1 , the responses are, in general, close to 1 and fall down more or less rapidly for higher dose rates, except DIS-100 that gives correct response up to 55 Sv.h -1 . Figure 5. Dose equivalent rate response of APDs in pulsed mode for a pulse frequency equal to 10 s -1 .

Effect of pulse frequency

Table 2 sums up the effect of pulse frequency as a percentage of variation in the APD response between 1 and 20 s -1 . This variation is roughly equal to 30% for all devices, except PM1621A (no signal) and EDD30 for which a saturation was observed from 2 Sv.h -1 .

Effect of pulse width

When the pulse width is larger than 1 s, the responses in pulsed and continuous radiation field are quite similar. No significant effect of pulse width was observed. Mean dose equivalent rate (Sv.h -1 ) All results from the pulsed field tests show that the more continuous the field tends to, that is to say the longer the pulses and the higher the frequency, the more satisfactory the behaviour of the devices.

H p (10)m / H p(10)ref DMC2000XB EPDMk2.3 EDMIII
These results show that it is important to add tests in pulsed mode when type-testing APDs, and thus the IEC 61526 [START_REF]Radiation protection instrumentation. Measurement of personal dose equivalent Hp(10) and Hp(0.07) for X, gamma, neutron and beta radiation: Direct reading personal dose equivalent and/or dose equivalent rate dosemeters[END_REF] standard should be revised.

Response of APDs in hospitals

The results of tests performed in hospitals are presented in Figure 6 as the distribution of APD reading normalized to the passive dosemeter reading. With respect to passive dosemeters, DMC 2000XB (median 0.65), EPD Mk2.3 (median 0.69), DIS-100 (median 0.78) and EDM III (median 0.81) on average, present an under-response. The behaviour of APDs is globally more satisfactory in hospitals than in laboratories because devices are mainly exposed to scattered fields.

CONCLUSION

The tests performed with continuous X-ray beams showed that all APDs have a satisfactory response at lower energies typical of IR/IC. Most APDs a correct response for dose equivalent rates up to 10 Sv.h -1 , except PM1621A for which the response is diverging rapidly from 1 Sv.h -1 and EDD30 which saturates above 2 Sv.h -1 . However, the dose equivalent rates in the direct beam can be much higher than those tested here. So these tests cannot guarantee that the APDs will correctly measure the high dose equivalent rates in the direct beam.

The study in pulsed mode showed that, except PM1621A, whose display does not give any indication, all APDs provide a reading. The tests performed in laboratory conditions with pulsed X-ray fields determined the effect of the dose rate, pulse frequency and pulse width on the APD response. First, for most APDs, the response is generally equal to 1 for dose equivalent rates lower than 2 Sv.h -1 , and decreases for higher dose rates, except DIS-100 that gives correct reading for high dose rates. Second, for a pulse frequency ranging from 1 to 20 s -1 , a variation of 30% in average is observed for all APDs. Finally, no significant effect of pulse width was observed. The measurements in hospitals confirmed an underresponse of APDs (median ranging from 0.65 to 0.81) with respect to passive dosemeters. Since all selected APDs, except PM1621A, provide a reading in pulsed mode, this means they could be used in routine monitoring at hospitals, provided correction factors are introduced [START_REF] Ambrosi | Considerations concerning the use of counting active personal dosemeters in pulsed fields of ionising radiation[END_REF] . Some preliminary guidelines with recommendations about how to use APDs in practice were presented at IRPA2010 [START_REF] Clairand | Guidelines for the use of APDs in interventional radiology[END_REF] .

  ): DMC 2000XB (MGPi), EPD Mk2.3 (Thermo), EDM III (Dosilab), PM1621A (Polimaster), DIS-100 (Rados), EDD30 (Unfors) and AT3509C (Atomtex).

Figure 1 .

 1 Figure 1. Active personal dosemeters tested in this study

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Energy response of APDs in continuous mode (for dose equivalent rate Hp(10) around 10 mSv.h -1 , and integrated dose equivalent around 0.5 mSv)

Figure 6 .

 6 Figure 6. Distribution of APD response compared with passive dosemeter response in realistic conditions

Table 1

 1

	Typical field characteristics encountered in
	interventional radiology/cardiology
	Parameter	Range
	High peak voltage	50-120 kV
	Intensity	5-1000 mA
	Inherent Al equivalent filtration	4.5 mm
	Additional Cu filtration	0.1 -0.9 mm
	Pulse duration	1 -20 ms
	Pulse frequency	1 -30 s -1
	Dose equivalent rate in the direct beam (table)	2 to 360 Sv.h -1

Table 2 .

 2 Effect of the pulse frequency (1 to 20 s -1 ): percentage of variation on the APD response for a dose equivalent rate = 1.8 Sv.h -1 and a pulse duration = 20 ms

	APD	DMC 2000XB EPD MK2.3	EDM III	PM1621A	DIS-100	EDD30	AT3509C
	Variation on the APD response (%)	25-30	30-40	<10	NO SIGNAL	30	10 (1.8 Sv.h -1 ) saturation from 2 Sv.h -1	30 (10-20 s -1 ) No signal at 1 s -1

Acknowledgement

This study has received funding from the European Atomic Energy Community's 7 th Framework Program (FP7/2007-2011 -grant agreement n°211361).


Luc Mouton

Alban Leroyer

Bo Gan

Deng

Patrick Queutey

Thibaut Soler

TOWARDS UNSTEADY APPROACH FOR FUTURE FLUTTER CALCULATIONS

Keywords: fluid-structure interaction, flutter, frequency domain flutter analysis, CFD, modal approach 𝑅𝑒 Reynolds number [-] CFD Computational Fluid Dynamics FSI Fluid-Structure Interaction RANSE

In prior work to define an improved hydrodynamic approach to flutter calculations, Centrale Nantes, Bureau Veritas Marine & Offshore and Farr Yacht design investigated the possibility of defining a linearized unsteady hydrodynamic model using a fluid response database coming from a series of 2D unsteady RANSE computations. The approach is compared to the Theodorsen theory. The linearized hydrodynamic model was used in a strip theory model for frequency domain flutter analysis. In this latest work, the IMOCA 2006 keel which has been used previously in frequential domain flutter calculation is also analyzed using an alternative and more accurate solution, featuring a fully coupled FSI modal approach with CFD.

As described in the literature, the results present a large impact of the unsteadiness on the phase and module for both lift and moment with a fairly good match compared to Theodorsen theory. The implementation of the results on a frequency domain flutter analysis tool reduces the critical speed for the studied model. Thus the results are closer to the 3D modal CFD approach which gave an lower critical speed.

NOMENCLATURE

INTRODUCTION

At the 2015 HPYD conference, Burns Fallow from North Sails noted how quickly we become accustomed to extraordinary. Indeed five years ago, it was extremely rare for boats to actually 'fly'; only Hydroptere and a few foiling dinghies had this potential. After the 2013 America's Cup the number of fully foiling designs exploded. We are now witnessing the first fully foiling ocean going yachts. It is relevant to consider whether flutter will be an issue for larger and larger foils. Are foils not concerned by the hydro-structural instability, or is it simply that the design envelope is still below the dangerous area? And what is the margin between the actual designs and the flutter limit?

As a reminder, an unstable hydro-structural vibration or flutter caused the failure of the keel of the IMOCA 60 SILL in 2002 [START_REF] Balze | Flutter of racing yacht keels and appendages[END_REF]. Since then several IMOCA 60's were subject to the phenomenon which resulted in expensive and time consuming redesign efforts. [START_REF] Mouton | Exploratory study on the flutter behavior of modern yacht keels and appendages[END_REF] presented the promising results of a modal hydro-structural flutter analysis. This approach used a very basic hydrodynamic model based on the heave and pitch motion of 2D hydrofoil strips. For structures such as keels, the model gave fairly accurate results and improved understanding of the phenomenon. The previously developed model can be used successfully to predict flutter in the sailing domain of a canting keel. For example this model can be used to produce optimized designs for yachts such as mini Transat 6.50 and maxis where the class rules do not require one design keels. In contrast other canting keel classes, like the VOR65 or the present IMOCA 60's have avoided flutter by requiring one-design steel keels.

Unfortunately foils and keels are very different structures in terms of their dynamic response. Hydro-structural effect induced by lead carriers (keels) is very different from the response observed on a high modal frequency foil. Hence, the first step was to study the difference between the two appendages and the available models for unsteady flows.

The present paper relates the effort produced by the LHEEA lab and Bureau Veritas Marine & Offshore to enhance the flutter analysis possibilities from a quasi-static approach to an unsteady model. The first part of the work deals with the unsteady strip model of the hydrostructural modal calculation. The second part describes preliminary results using a fully coupled modal approach with CFD.

Flutter

The phenomenon of flutter is well known in the aviation industry. It was discovered at the beginning of the twentieth century on aircraft wings and was shown to be linked to a vibration problem. During flutter, a coupling of torsional and bending vibration modes of the structure with the aerodynamic forces leads to a transfer of energy from the fluid to the vibrating structure.

Below the flutter speed, and at low vibration frequencies, the hydrodynamic lift induced by a heave motion (bending) is in opposition of phase with the motions speed. It damps the vibration by transferring energy from the vibrating structure to the fluid. The lift induced by the torsional motion is in phase with the motion. It transfers energy from the structure to the fluid and reversely at each period. It is close to be energy conservative. In such condition, the damping is positive and any vibration of the structure reduces.

Approaching flutter speed, with low vibration frequencies, the coupling of the bending and torsional modes allows a mode to be combination of non-phased torsional and bending motions. During a vibration the lift of the torsional component may have a phase different from the bending motion. Thus it can bring energy to the bending component. The damping then becomes null or even negative. This transfer of energy, even for a few seconds, can be enough to induce vibration amplitudes that can lead to complete failure of the structure. The phenomenon is said to be explosive.

Lead bulbs and foils, quasi-static and unsteadiness

The main objective of a modern keel is to carry a heavy lead bulb that will produce a large part of the overturning moment necessary to keep the boat upright under the massive force of the sails. If we consider a typical 2006 IMOCA keel design that experienced the onset of flutter, the bulb inertia lowers the natural frequencies of the keel to:

-1.0 Hz for the bending mode, -2.15 Hz for the torsion mode.

The appendage modal vibration creates a harmonic oscillation of the fin's section. A vibration of the torsional mode will see a full period of fin section oscillation in 0.46 second. During this time, a particle of water travelling at 30 knots will travel a distance of 7.2 m. This means that after one oscillation, the particle has travelled 18 chord lengths downstream. In other words, the vortex wake induced by an oscillation of the profile will be very far aft of the profile at the end of the initiating oscillation. As the influence of the vortex is inversely proportional to the distance, it will have a small impact on the profile. The hydrodynamic response can be considered quasi-static. This is characterised by a low reduced frequency 𝑘 equal to 0.3 in this case.

On the other hand, a light and stiff foil will have natural frequencies classically around 20 Hz. For the same chord and speed, a particle of water will travel only 1.1 chord lengths downstream during an oscillation. The vortex wake induced by the oscillation will exist directly downstream of the section of the foil at the end of the period. This proximity creates a potential influence of the vortex wake on the hydrodynamic forces of the next oscillation.

The flow is said to be unsteady. The reduced frequency of this example is 5.3. The number of chord lengths travelled by a particle of water (or a vortex wake) during a period of the oscillation is . This explains why the quasi-static hydrodynamic model can give good results in the flutter calculation of a ballasted keel. However, to enable the faithful modelling of hydro-structural effect on a foil, it is necessary to take into account unsteady hydrodynamic effects.

Although turbulence models mainly designed from steady cases still need validation for such FSI configurations, the unsteady effects are intrinsically considered in a URANSE code. However, the use of complex CFD tools for flutter evaluation is cumbersome as the creation of such a model is complex and time consuming. Therefore, the frequency domain model remains interesting but requires enhancement of its hydrodynamic models to consider high reduced frequency problems.

In order to better take into consideration the unsteady effects, two options are considered and compared. The first is the analytical Theodorsen model. This model considers a 2D wing section harmonically oscillating in a flow. As a function of the reduced frequency, it evaluates the lift amplitude and phase lag. Momentum can also be expressed, but not drag.

The second option is to build a 2D unsteady model representing a foil profile pitching or heaving in a uniform flow produced by CFD. By simulating the fluid forces produced by different sets of harmonic motion in a database, the latter can be used as input to the frequency domain tool in order to incorporate unsteady effects into the fluid model.

The methodology is then based on two tools: the first one, developed by Bureau Veritas takes its background from the aeronautical industry. It is based on the frequency domain analysis of a hydro-structure model, the results of the 2D unsteady calculations are integrated via a strip theory model as hydrodynamic models of the code [START_REF] Mouton | Exploratory study on the flutter behavior of modern yacht keels and appendages[END_REF]. The second one, ISIS-CFD, is a Navier-Stokes solver developed at the LHEEA Lab., which is used here both to build a database as input to the previous tool and also to carry out a fully coupled time-marching resolution of this fluid-structure interaction problem.

THE FREQUENCY DOMAIN TOOL

Modal hydro-structural flutter analysis

The structural model is based on a truncated series of the n most energetic modes, without limitation on the number of modes. Accounting for two modes in a keel analysis is sufficient, but would be insufficient to analyze properly light structures such as a hydrofoil. The inputs to the flutter analysis software is a chosen number of n dry modes natural frequencies and shapes (normalized using the diagonalized mass matrix) provided as an output from Finite Element Analysis (FEA) of the structure.

Laws of dynamics give the following equation for the dry system with or without heel:

𝑀𝑋 ¨+ KX = 𝐹 Gravity (1)
As a reminder, and as detailed in [START_REF] Mouton | Exploratory study on the flutter behavior of modern yacht keels and appendages[END_REF], the strip model consists of several slices of the fin profile along the span. During any vibration of the keel or foil, the motion of each slice can be described with 6 degrees of freedom (see figure (2)) among which two are identified as having a sensible impact on the hydrodynamic flow.

In the local axis of a fin section (span wise, chord wise and perpendicular to the mean plane) the motion of the section profile can be named heave and pitch as described in figure (3). The hydrodynamic loads in the respective direction are the drag, lift and pitching moment.

For the sake of simplicity, only the lift will be shown here, but it should be noted that what follows can be reproduced for the other hydrodynamic loads.

2D flat foil section harmonically oscillating in pitch and heave

Let's consider a 2D profile section oscillating of small variations 𝑥 ̃ around an equilibrium position 𝑋 0 :

𝑋 = 𝑋 0 + 𝑥 ̃ (2)
If we consider harmonic motions around the equilibrium position, the motion can be written as equation ( 3).

𝑥 ̃(𝑡) = [ 𝜃 ̃(𝑡) = ℜ𝑒(𝜃 0 𝑒 𝑖𝜔(𝑡-𝑡 𝑝 ) ) ℎ ̃(𝑡) = ℜ𝑒(ℎ 0 𝑒 𝑖𝜔(𝑡-𝑡 ℎ ) ) ] (3)
Assuming the uncoupled hydrodynamic response we can write the unsteady lift per unit span through equation ( 4).

𝐿 ̃= 𝐿 ̃pitch (𝜃 ̃, 𝜃 , 𝜃 ) + 𝐿 ̃heave (ℎ ̃, ℎ ̇) (4)

In the hypothesis of small motions around the equilibrium position, we make the assumption that the hydrodynamic response is harmonic and its spectral content is purely in the first harmonic of the initial motion, leading to equation ( 5 .

Note that for very slow dynamic (𝑘 ≪ 1) the result should be in accordance with the potential theory. The phase of the pitch motion lift shall tend to 0 for slow motions (lift in phase with pitch). Then we should have for the derivative of the lift coefficient:

‖ 𝜕𝐶 𝐿 pitch 𝜕𝜃 ‖ , ‖ 𝜕𝐶 𝐿 heave 𝜕𝛼 ℎ ‖   0 k 2𝜋 and 𝜙 𝐿 𝑝 , 𝜙 𝐿 ℎ   0 k 0 (7)
It should be noted that the heave motion produces a lift which will tend to be in phase with the heave induce angle 𝛼 ℎ , therefore 𝜙 𝐿 ℎ   0 k 0. This leads to -𝜋 2 ⁄ of phase lag with the heave motion itself.

The unsteady effect will act on the amplitude and the phase lag of the lift response. It is important to note that any steady and unsteady effect (wake, added mass) can be modelled through this approach. It is therefore a practical way to compute a very complete hydrodynamic model for harmonic motions. Following the same approach we get for the pitching moment:

2D CFD Calculations

The objective of the 2D CFD analysis is to compute the amplitude 𝐿 0 pitch (𝑘,Re) and the phase 𝜙 𝐿 𝑝 (𝑘,Re) (or the equivalent for the heaving motion) of the hydrodynamic response.

Each motion (pitch or heave) is imposed to a profile in a constant infinite flow defining a couple (𝑘,Re). The unsteady calculation is run and the lift computed for each time step. Thus, the Fast Fourier Transformation of the resulting signal is achieved. Only the first harmonic is kept here assuming that it contains the major part of the response. This hypothesis is validated in the next section.

To build the database, the profile shape from the 2006 Farr Yacht Design IMOCA 60 fin is used (see figure (4)). Computations were carried out using the ISIS-CFD solver described in section Navier-Stokes solver. The 2D fluid domain of 51610 cells is 12 chord lengths and 10 chord heights. Far field conditions are used for all external boundaries, except the outlet where a constant pressure is imposed. A wall-function approach associated with the 𝑘 -𝜔 SST Menter turbulence model is used [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. Grid and time step independency were checked before running the complete database. The domain of calculation was 9.6e4 to 1.0e7 for the Reynolds number, and from 0.1 to 6 for the reduced frequency.

Theodorsen solution

The model of Theodorsen [START_REF] Theodorsen | General theory of aerodynamic instability and the mechanism of flutter[END_REF]Hodge and Pierce 2002) is derived from a theory of unsteady aerodynamics for a thin 2D airfoil undergoing small harmonic oscillations in incompressible and inviscid flow.

It contains added-mass and quasi-steady contributions as well as the effect of the wake, included in the Theodorsen's function 𝐶(𝑘), for both lift force and momentum. No information concerning the drag can be deduced from this model.

Since the publication of [START_REF] Theodorsen | General theory of aerodynamic instability and the mechanism of flutter[END_REF], this unsteady model of fluid loads and other variants have been widely used as pre-design for the flutter analysis and fatigue phenomenon of helicopter rotor blade [START_REF] Motta | Computational fluid dynamic analysis of a L-shaped Gurney flap for vibration control[END_REF] , wind turbine [START_REF] Leishman | Challenges in modelling the unsteady aerodynamics of wind turbines[END_REF][START_REF] Lobitz | Aeroelastic stability predictions for a MW-sized blade[END_REF]) and suspension bridges (Agar, T. J. A 1989 & Banerjee, J.R. 2001).

Hunsaker, Phillips ( 2015) also used it to analyze the efficiency of flapping foils and compare the results with CFD.

With the notations used in section Modal hydro-structural flutter analysis and in figure (3), the expressions of lift and momentum (expressed at the center of rotation) are given by equation ( 9).

𝐿(𝑡) = 2πρ𝑉 ∞ bC ( 𝑘 2 ) [-ℎ ˙+ 𝑉 ∞ 𝜃 + 𝑏 ( 1 2 -𝑎) 𝜃 ˙] + 𝜋.𝜌.𝑏 2 (-ℎ ¨+ 𝑉 ∞ .𝜃 ˙-𝑏.𝑎.𝜃 ¨) 𝑀 𝑂 (𝑡) = -πρ𝑏 3 [- 1 2 ℎ ¨+ 𝑉 ∞ 𝜃 ˙+ 𝑏 ( 1 8 - 𝑎 2 ) 𝜃 ¨] (9)
with

𝐶(𝑘) = 𝐻 1 (2) (𝑘) 𝐻 1 (2) (𝑘)+𝐻 0 (2) (𝑘) , 𝐻 𝑛 (2)
(𝑘) being the n th Hankel function of second species.

THE NAVIER-STOKES SOLVER

Brief description of the ISIS-CFD code

The ISIS-CFD code, developed by the METHRIC team of the LHEEA Lab. of Ecole Centrale Nantes, UMR-CNRS 6598, solves the incompressible Unsteady Reynolds-Averaged Navier-Stokes (URANS) equations in a strongly conservative way. It is based on a fully unstructured, cell-centered finite-volume method to build the spatial discretization of the conservation equations. Pressure-velocity coupling is obtained through a Rhie & Chow SIMPLE-type method: in each time step, the velocity update comes from the momentum equations and the pressure is given by the mass conservation, transformed into a pressure equation. The temporal discretization scheme is the Backward Difference Formula of order 2 (BDF2) when dealing with unsteady configurations. For each time step, an inner loop (denoted by non-linear loop) associated to a Picard linearization is used to solve the non-linearities of the system and converge all the sequential coupled equations.

In the case of turbulent flows, additional transport equations for modeled variables are solved in a form similar to the momentum equations and they can be discretized and solved using the same principles. An Arbitrary Lagrangian Eulerian (ALE) formulation is used to take into account modification of the fluid spatial domain due to body motion and deformation [START_REF] Leroyer | Experimental and numerical investigations of the flow around an oar blade[END_REF]. Only monofluid configurations are investigated here. However, free surface flowcan be addressed with an interface capturing method, by solving a conservation equation for the volume fraction of water, with specific compressive discretization schemes [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF]. The code is fully parallel using the MPI (Message Passing Interface) protocol.

Fluid-Structure Interaction coupling with modal approach

The modal approach in the ISIS-CFD code uses the n first natural modes of vibration of the dry structure as input. Since the eigen modal matrix makes the mass matrix and stiffness matrix orthogonal, the temporal resolution of the structure is then reduced to a set of uncoupled degree of freedom (DOF) governed by the n modal equations given by equation (10).

𝑞 ¨𝑖 + 2𝜀 𝑖 𝜔 𝑖 𝑞 ˙𝑖 + 𝜔 𝑖 2 𝑞 𝑖 = 𝜑 𝑖 𝑇 𝑓(𝑡),𝑖 ∈ [1, 𝑛] (10)
where 𝜑 𝑖 is the i th modal vector normalized by the mass, 𝑞 𝑖 the amplitude of the mode, 𝜔 𝑖 = 2𝜋𝑓 𝑖 the pulsation, 𝑓 𝑖 the frequency, and 𝜀 𝑖 is a possible damping coefficient assuming the Rayleigh damping hypothesis (damping matrix proportional to the mass and stiffness matrix). 𝑓(𝑡) refers to the fluid force acting on the structure. The total shape deformation is then given by equation ( 11).

∑ 𝑞 𝑖 (𝑡)𝜑 𝑖 (𝑥) 𝑛 𝑖=1 (11)
These equations are fully coupled with the resolution of the fluid flow through an internal implicit coupling. At each non-linear iteration of the fluid solver (within the same time step), all the DOF of the structure are solved and updated. A Radial Basis Function (RBF) approach is used here to compute the source term of equation ( 10) and as a mesh deformation technique to recover a body-fitted mesh after each resolution of the structure.

This approach was previously validated by extracting some configurations from Hollowel (1981) and [START_REF] Verma | Flutter prediction for a composite plate using unsteady CFD based FSI solver[END_REF], based on aerodynamic flutter.

2D RESULTS USED TO BUILD THE DATABASE

Verification of the harmonic content

Using the FFT, it was verified that all disregarded harmonics are less than 5% of amplitude of the fundamental mode.

Verification of the linearity

The first step was to validate the assumption of linearity. Additive and multiplicative linearities were checked. This was performed for various reduced frequencies. Results are very good as shown in figure (6). In this figure, we observe good agreement between the hydrodynamic force generated by a combined motion, and the reconstructed force 2L1+2L2 generated with separated pitch L1 and heave L2 motions at medium to high dynamic motion (k=2.51 et Re=3.96e5).

The slight difference can be attributed to the amplitude of the motion which becomes larger with the combination of both pitch and heave associated with the factor 2 for the amplitude. 𝜕𝐶 𝐿 𝜕𝜃 are identical for the quasi-static regime, large differences appear for high dynamic regime. The unsteady effect has a large impact on the lift curve coefficient as the minimum is more than 30% below quasi-static value, and for large value of 𝑘 it increases continuously, (figure ( 8)). Regarding the Theodorsen model, the body motion must produce small variations of velocity compared to the inflow velocity to fulfill the hypothesis of small disturbances. As a consequence, this approach becomes less and less valid when the reduced frequency increases. This explains why we can observe larger differences when the reduced frequency increases.

We observe that both approaches follow very similar trends. For 𝑘 = 0 both curves tend to the potential lift curve slope coefficient 2𝜋 of a thin plate assumed by the Theodorsen theory. It should be noted that the CFD calculation presents the actual foil profile and therefore tends to the actual value which is close but not equal to 2𝜋. In figure ( 9), the phases present a good agreement between the two methods. For low values of 𝑘 phases are in line with the quasi-static, 𝜙 𝐿 𝑝 , 𝜙 𝐿 ℎ tend to 0, which traduce no phase delay with 𝜃 ̃ the pitch motion, and α ̃ℎ the heave induced angle. Concerning the pitching moment, a good agreement can also be noticed between CFD and Theodorsen theory, figures (10,11). Again the quasi-static value is matched as the pitching moments tend to zero on this profile where the rotation center is defined at one quarter of the chord 𝐵 from the leading edge. This is verified for both pitch and heave. Similar to the lift, the quasi-static pitching moment is not a proper model to define a configuration involving high dynamic phenomena. Regarding the phase, figure (11), for large values of reduced frequency, the models agree fairly well. However, there are large differences at low frequencies. Even though the impact is probably limited as it happens as the amplitude tends to zero, the reason for these discrepancies is not known at this point. It may be related to the viscous friction that becomes more important than the pressure loads. The thickness of the profile may play a role since the Theodorsen model only considers a zero-thickness profile. Also, for high values of 𝑘 the presented results can be influenced by the lower Reynolds number used compared to the ones used for low values of 𝑘, see figure (7).

It can be noticed the present CFD results which uses a quite thick profile (see Fig. 4) match well with the trends obtained by [START_REF] Motta | Computational fluid dynamic analysis of a L-shaped Gurney flap for vibration control[END_REF] which study the influence of the thickness of the profile. Even if the study was restricted to a range of reduced frequency from 0 to 0.75, the results of [START_REF] Motta | Computational fluid dynamic analysis of a L-shaped Gurney flap for vibration control[END_REF] for quite thick NACA profile (NACA 0018 and NACA 0024) are very similar to the one presented here concerning the position of the curves relative to the Theodorsen ones for the amplitudes and for the phase. It can be noticed that the influence of the thickness of the profile is not negligible.

It seems that the hydrodynamic focal point at a quarter of the chord could be reconsidered at high reduced frequencies as the moment increases largely. However, this results should be further translated in terms of position of the focal point.

2D unsteady hydrodynamic calculations, by nature induce 2D unsteady wakes that tend to produce larger influence of the unsteady effects. In real 3D structure, the wake and vortices will fold in 3D and their influence tend to be reduced. Thus the 2D hydrodynamic model integrated along a fin to represent a 3D flow represents rather overestimated effect of the unsteadiness. 12), the slight differences change the flutter speed from 20.9 down to 18.9 knots.

As expected no tremendous differences are found on this keel. This is due to the low natural frequencies that induce an almost quasi-static hydrodynamic regime. Here the quasi-static model is sufficient to model correctly the physics of the flow. However, even for such a low reduced frequency the flutter speed is reduced by 10% which is not negligible.

It is different for a lighter structure, and the unsteady effects will have a larger influence. The method is now ready for this future step. Next steps are to include pitching moment into the code, and to allow the definition of curved foils. FSI coupled modal approach with CFD For this configuration, simulations have been done with a full hexahedron mesh of 3.8 million cells designed for a wall-function approach of the turbulence (y+=30). Concerning the structural part, the information of two first dry modes used for the frequential hydro-elastic analysis was simply converted into an input data file formatted for the CFD code. At the beginning of the simulation, the flow is at rest. The body is sped up to reach its nominal forward velocity within a short time (1s). FSI resolution is activated from the start. Perturbations generated during this violent transient state excite the structure modes. Thereafter, we can follow the evolution of their amplitudes in time and then the evolution of the displacement at some locations of the keel.

Results show that the approach succeeds in catching the flutter occurrence for a velocity around 8 m/s (figure ( 14)). For velocity of 7 m/s, it can be shown that the time derivative of the rear bulb tip displacement tends to slowly decrease, whereas at 9 m/s, the flutter velocity is exceeded, which results in an increase of the velocity in time. At 8 m/s, close to the flutter critical velocity found by this approach, the response tends to a limit periodic cycle. The range of the flutter velocity found here is a slightly lower than the value found by the hydro-elastic tool. However, additional simulations should be run to assess numerical convergence (grid, time step,...).

The fully coupled time marching approach between modal approach and RANSE solver is capable of detecting flutter phenomenon even if it requires large CPU time to accurately determine the range of the flutter velocity.

It is interesting to consider the differences between the methods presented in this work. As the dichotomic approach of the CFD 3D simply furnished upper and lower bound of the flutter critical speed we will consider the average value for this technique and compare to the frequency domain approach. Figure 15 summarizes the results. It can be seen that there is quite a large difference (around 26%) between the quasi-static hydrodynamic frequency approach and the 3D CFD modal calculation. This is reduced to 14% when considering the 2D unsteady linearized hydrodynamic in the frequency tool based on the CFD matrix response of the 2D profile.

The actual flutter speed on this boat was estimated around 18 to 20 knots. This suggests that the frequency domain approach is more precise than the modal CFD. However, in both analysis, the structural damping was not modelled. As explained in [START_REF] Mouton | Exploratory study on the flutter behavior of modern yacht keels and appendages[END_REF] the structural damping pushes the flutter speed higher as the negative damping from the hydrodynamics need to counter the positive structural damping. The impact on the critical speed was not assessed in the present work. Not considering the structural damping is conservative. Thus the quasi-static frequency domain approach is not conservative while CFD is probably the most accurate. The linearity of the model requires verification, and is only valid for small motions around the equilibrium position. However, it can include unsteady effects, such as phase delay, amplitude modifications, and added mass in a very simple linear approach. The comparison with Theodorsen's theory allows good confidence in the results. The presented database is sparse, and further work is required. However, the method allows the use of the frequency flutter analysis tool for foils with high dynamic regime even if only driven by the Theodorsen model maybe corrected for thickness as it is suggested in [START_REF] Motta | Computational fluid dynamic analysis of a L-shaped Gurney flap for vibration control[END_REF]. In the present test case the 2D unsteady model presents a lower flutter speed. However it is not sufficient to generalize.

The methodology is efficient to evaluate the onset of flutter which should be avoided in a design approach. In some cases limit cycles stabilize the flutter effect for large vibration amplitude. This constitutes a limit case and many further investigations will have to be performed prior to accept such effect in the sailing domain. On the contrary, a comfortable speed margin should be used.

The modal approach that is fully coupled with CFD, contains fewer assumptions and captures more of the physics, such as 3D effects. This approach requires running several unsteady simulations at different speeds via a dichotomic search to determine the critical flutter velocity range. More accurate estimations of the flutter phenomenon can be expected, at the expense of far larger computational power.

The tools are complementary; on the one hand the frequency analysis is able to investigate the onset of flutter for a large number of designs in a limited amount of time, on the other hand the FSI computation can provide a precise verification including a broader range of physical phenomenon. It can be noted that the use of linearized unsteady 2D hydrodynamic model in the frequency domain provides a critical speed closer to the URANSE calculated value.

  by pitch 𝜃 ̃ and heave 𝛼 ̃ℎ, [N/m] 𝑀 𝑂 (𝑡) Pitch moment at point O per unit of span [N] 𝑀, 𝐾 Mass and stiffness matrix O Rotational center of the profile Pdyn Dynamic pressure 1/2V∞² [N/m 2 ] 𝑞 𝑖 Amplitude of the i th dry mode of the structure [-] 𝑡 𝑝 , 𝑡 ℎ Delay of the pitch and heave motion compared to the time origin [s] 𝑉 ∞ Boat speed [m/s] 𝑋 Vector characterizing the position of the keel [m; rad] 𝑋 0 Equilibrium position [m; rad] 𝑥 ̃ Variation of 𝑋 around an equilibrium position 𝑋 0 , vector characterizing the deformation of the keel [m; rad] the lift induced by pitch and heave relative to the motion, [rad] 𝜑 𝑖 i th modal vector normalized by the mass 𝜃 Angle of incidence (pitch angle) of a foil section [rad] 𝜃 0 Amplitude of the harmonic pitch angle of a foil section [rad] 𝜌 Density of water [kg/m 3 ] 𝜔 Pulsation of the motion [rad/s]

Figure

Figure 1 .

 1 Figure 1. Illustration of the reduced frequency and indication of the potential impact of the vortex induced sheet

Figure 2 .

 2 Figure 2. Foil section with three possible translation (left) and three possible rotation (right). Motions with significant impact on the hydrodynamic are plotted in red.

Figure 3 .

 3 Figure 3. Section profile motions in the reference axes of the undeformed foil.

Figure 4 .

 4 Figure 4. Close view of the profile mesh

Figure 5 .

 5 Figure 5. Modes evolution in Reduced frequencies versus Reynolds number and associated 2D CFD calculations points,

Figure 6 .

 6 Figure 6. Example of verification of lift force multiplicative and additive linearity

Figure 7 .

 7 Figure 7. Reduced frequency versus Reynolds points computed in 2D unsteady CFD Comparisons for lift are given in figures (8, 9). Even though 𝜕𝐶 𝐿 𝜕𝛼 ℎ ;

Figure 8 .Figure 9 .

 89 Figure 8. Module of the lift curve slope as a function of the reduced frequency (left: pitch motion, right: heave motion )

Figure 10 .Figure 11 .

 1011 Figure 10. Module of the moment curve slope as a function of the reduced frequency (left: pitch motion, right: heave motion )

Figure 12 .

 12 Figure 12. On the left: Evolution of the natural frequencies of the keel as a function of: speed (top); corresponding reduced frequencies (bottom). On the right: Damping rate as a function of the boat speed, comparison between quasi-static and linearized unsteady hydrodynamic, (flutter speed is very close to 20 knots)

  Figure (13) show a global and a closer view of the mesh. The time step is set at 0.002 s. The Reynolds stress tensor is computed using the k-w (SST-Menter) model,[START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF].

Figure 13 .

 13 Figure 13. CFD mesh generated for the IMOCA's keel (left: global view, right: zoom on the keel)

Figure 14 .

 14 Figure 14. Evolution the transversal velocity at the rear bulb tip for different velocities

Figure 15 .

 15 Figure 15. Comparison of the flutter critical speed with the different approaches

Acknowledgements

This work was granted access to the HPC resources of CINES and IDRIS under the allocation 2B0129 made by GENCI (Grand Equipement National de Calcul Intensif).


Anne Auger
email: anne.auger@lri.fr

Marc Schoenauer
email: marc.schoenauer@inria.fr

Nicolas Vanhaecke
email: vanhaeck@fhi-berlin.mpg.de

LS-CMA-ES: a Second-order algorithm for Covariance Matrix Adaptation

Evolution Strategies, Evolutionary Algorithms based on Gaussian mutation and deterministic selection, are today considered the best choice as far as parameter optimization is concerned. However, there are multiple ways to tune the covariance matrix of the Gaussian mutation. After reviewing the state of the art in covariance matrix adaptation, a new approach is proposed, in which the covariance matrix adaptation method is based on a quadratic approximation of the target function obtained by some Least-Square minimization. A dynamic criterion is designed to detect situations where the approximation is not accurate enough, and original Covariance Matrix Adaptation (CMA) should rather be directly used. The resulting algorithm is experimentally validated on benchmark functions, performing much better than CMA-ES on a large class of problems.

Introduction

Among the class of Evolutionary Algorithms (EAs), Evolution Strategies (ESs), based on the so-called Gaussian mutations, are today considered the state-ofthe-art in parameter optimization, i.e. optimization of a function defined on R n for some n ≥ 1 [START_REF] Th | An Overview of Evolutionary Algorithms for Parameter Optimization[END_REF]. Basic ESs generate λ offspring by mutating the µ parents (without selection), and then deterministically replace the parents by the µ best individuals from either the offspring (this variant is then called the (µ, λ) -ES), or the merge of parents and offspring (and it is then called the (µ+λ)-ES). The Gaussian mutation operator for this large class of algorithm consists in adding some Gaussian random noise to the parent, and can be formally written as µ σ,C : X 0 → X 0 + σN (0, C).

(

) 1
where N (0, C) represents one realization of a multivariate Gaussian distribution with mean 0 and covariance matrix C, a symmetric positive definite matrix, and σ is a positive real value called the step-size. This distinction can seem arbitrary (the actual covariance matrix is σC), but it is often convenient to consider that C is somehow normalized, giving the "direction" in which sampling should occur, while σ indicates how far the mutation should go in that "direction". The performance of the EA is of course highly dependent on the choices of C and σ, that have to be tuned not only to the problem at hand, but also to the current state of the evolution process. Several techniques have been proposed, from the self-adaptation of the mutation parameters in Evolution Strategies (SA-ES) [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF] to the Covariance Matrix Adaptation (CMA-ES) [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] where the covariance matrix C is deterministically adapted from the last move of the algorithm. Some Estimation of Distribution Algorithms [START_REF] Larranaga | Estimation of Distribution Algorithms[END_REF] can also be considered as pertaining to the same category of algorithms that repeatedly sample some Gaussian distribution to generate next test points.

However, we argue in this paper (Section 2) that none of those algorithms does make an optimal use of all previously sampled points . We then discuss, in section 3, what the optimal covariance matrix should be, based on the case of elliptic functions. In Section 4 we propose the LS-CMA-ES algorithm, an ESlike algorithm that updates the covariance matrix by learning some curvature information of the fitness function. Finally, Section 5 presents and discusses experimental results on several benchmark functions from the literature, while section 6 derives from those results further work that must be done to extend the proposed method (e.g. to better handle multi-modal cases) and to complete its assessment.

State of the Art in Covariance Matrix Learning

Since the early work of Rechenberg and Schwefel, who first used Gaussian mutations to evolve real-valued parameters (see e.g. [START_REF] Bäck | A survey of Evolution Strategies[END_REF] for a detailed history of the field), it has been clear that the most critical issue of ESs was the tuning of the mutation parameters. Initial works were concerned mainly with the stepsize: Rechenberg's one-fifth rule [START_REF] Rechenberg | Evolutionstrategie: Optimierung Technisher Systeme nach Prinzipien des Biologischen Evolution[END_REF], derived from theoretical considerations on the progress rate, as well as Schwefel's first self-adaptive ES [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF] were only concerned with the optimal step-size, i.e. considered the covariance matrix to be the identity matrix. However, Schwefel soon proposed to self-adapt one step-size per variable, i.e. to use a covariance matrix that is diagonal with positive diagonal terms, and further extended the self-adapted parameters to the so-called correlated mutations, corresponding to the general case of Equation 1 where all parameters of the covariance matrix C are self-adapted.

The ground idea of self-adaptation for correlated mutations is to rely on mutations of the mutation parameter themselves to adjust the covariance matrix C. Whereas this clever idea removed the need to manually tune the mutation parameters, it does not take into account the history of evolution to direct the search -as did the one-fifth rule for the step-size. From thereon, Hansen and Ostermeier designed the Covariance Matrix Adaptation-ES (CMA-ES) algorithm [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF], later made more computationally efficient in [START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF]. The basic idea of CMA is to deterministically use the the most recent descent direction, i.e. the direction between the best parents at two consecutive generations: the covariance matrix C is gradually updated with rank-one matrices whose unique non-zero eigenvalue has as eigenvector the last descent direction (see the complete equations at lines 8 and 10 of Table 1). In some sense, whereas CMA has been thought as a derandomized correlated mutation from the ES point of view, it could also be viewed as a randomized steepest-descent algorithm from the point of view of standard numerical optimization. Note that the complete CMA algorithm also includes some derandomized step-size adaptation. However, it has been clearly experimentally demonstrated (not shown here for space reasons) that this stepsize adaptation is not crucial, as the performances of the algorithms are roughly the same if the step-size is self-adapted using the standard SA-ES mechanism (see [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF], or line 4 of Table 1, for the formal equation).

Estimation of Distribution Algorithms (EDAs) is a recent paradigm using similar ideas -at least in the continuous case (see [START_REF] Larranaga | Estimation of Distribution Algorithms[END_REF] for a survey of the whole field). EDAs evolve a probability distribution over the search space. The main loop of an EDA first samples that distribution, then selects some of the samples according to their fitness, and finally updates the distribution from those selected individuals. Hopefully, the distribution will gradually concentrate on the optima of the fitness function. Of course, the update step depends on the model that has been chosen for the distribution. A typical EDA for parameter optimization is the EM N A family [START_REF] Larranaga | Estimation of Distribution Algorithms[END_REF], where the distribution is sought as a multi-variate Gaussian distribution: the sampling phase of EM N A using such a distribution is then equivalent to the birth of offspring using Equation 1, and EMNA does adapt a covariance matrix, similarly to the ES algorithms described above.

On the one hand, EMNA uses many points to refine the covariance matrix -and hence one could think that it will get a better approximation of the "optimal" covariance matrix, if such a thing exists (see Section 3). However, the selection phase of EMNA (and of most EDAs) is merely binary: if for instance the best half of the population is selected, individuals get selected without any reference to their respective fitness. Moreover, in the case where a completely new covariance matrix is generated from the selected samples, using some Maximum Likelihood principle, such a covariance matrix only tries to "explain" the distribution of sample at hand, without paying any attention to other parts of the search space, even those that have been visited before. And this situation is only slightly improved if the covariance matrix is updated from the previous one at each generation.

Experiments have shown that the best-performing of those algorithms is clearly the CMA-ES algorithm (see [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] for comparison of SA-ES and CMA-ES, and [START_REF] Hansen | Learning Probability Distributions in Continuous Evolutionary Algorithms -A Comparative Review[END_REF] for a comparison of CMA-ES with some EDAs) -and the reason for that seems to be that the covariance matrix is better adapts to teh fitness landscape. The purpose of this paper is to try to make a better use of all the points that have been sampled by the algorithm, when updating the covariance matrix, by learning the local curvature of the target function. But before that, we need to clarify what could be the ideal covariance matrix within an algorithm using Gaussian sampling.

Rationale

The simplest problem, on which the behavior of optimization algorithm can be studied in great detail, is the sphere model, where the goal is to minimize the function f S defined by

f s (x) = x T x. (2)
It has been shown numerically [START_REF] Th | An Overview of Evolutionary Algorithms for Parameter Optimization[END_REF] that the optimal covariance matrix in that case is the identity matrix I d . Moreover, some theoretical studies have proved the convergence of the on (1, λ) -ES, either with a dynamic step-size [START_REF] Auger | Dimension-independent Convergence Rate for Non-isotropic (1, λ)-ES[END_REF] or in the case of self-adaptive step-size [START_REF] Auger | Convergence Results for (1,λ)-SA-ES using the Theory of ϕ-irreducible Markov Chains[END_REF]. From those results, plus the naive consideration of the isotropy of the fitness function itself, we shall from now on suppose that indeed I d is a good, if not optimal, covariance matrix when optimizing the sphere function.

But suppose that we now want to minimize the elliptic function

f E f E (x) = 1 2 x T Hx, (3)
where H is a symmetric positive definite matrix. The choice of the optimal covariance matrix to solve that problem with Evolution Strategies becomes obvious after a change of variables that will turn this problem into a sphere problem. Consider the eigen decomposition of H in an orthonormal basis:

H = P T ∆∆P ( 4
)
where P is an orthonormal matrix (P -1 = P T ) and ∆ a diagonal matrix whose diagonal terms are the square roots of the eigenvalues of H. Now let W = ∆ -1 P X. Then simple calculations show that f E (X) = f S (W ) and that the mutation operator given by equation ( 1) with C = ( 1 2 H) -1 transforms W 0 into W 0 +σN (0, I d ). Hence, if we consider that I d is the best choice for the covariance matrix for the sphere problem, then ( 1 2 H) -1 is the best choice for the covariance matrix for the elliptic problem.

For arbitrary functions, But on the other hand, thanks to Taylor expansion, all regular functions can be locally approximated by an elliptic function:

f (X) = f (X 0 )+(X -X 0 ) T ∇f (X 0 )+ 1 2 (X -X 0 ) T H(X 0 )(X -X 0 )+o(||X -X 0 || 2 )
(5) where ∇f (X 0 ) denotes the gradient of f at point X 0 and H(X 0 ) the Hessian matrix of f at point X 0 . It hence makes sense to try to use within ESs an approximation of the Hessian matrix of the function at hand within the mutation operator. But before proposing a deterministic algorithm using points previously sampled during evolution to approximate that Hessian matrix, let us first look at existing methods to derive a covariance matrix in the evolutionary community.

The LS-CMA-ES algorithm

The goal of this section is to introduce an algorithm that will try to use an approximation of the Hessian matrix of the target function as advocated in section 3: This approximation will be deterministically built using the available sample points -and later used within the Gaussian mutation 1.

Whereas such strategy is optimal for elliptic functions, it can fail on functions that are "more flat" (e.g. the famous Rosenbrook function), because the quadratic approximation then becomes very inaccurate. A criterion needs to be designed to detect such cases and to alternatively use the CMA "steepest descent update rule" for the covariance matrix.

Approximating the Hessian matrix

This step of the proposed algorithm aims at learning the local Hessian matrix of the target function. Starting from the Taylor expansion (see equation 5), and supposing that we have some sample points of values taken by f , i.e. a set of points X k together with their actual fitness f (X k ), located "not too far" from a given point X 0 , approximations g and H for both ∇f (X 0 ) and H(X 0 ) can be easily obtained by solving the linear least-square minimization problem:

min g∈R d ,H∈S(R d ) N k=1 f (X k )-f (X 0 )-(X k -X 0 ) T g- 1 2 (X k -X 0 ) T H(X k -X 0 ) 2 (6)
The unknowns of this problem are g (n unknown parameters), and H (n(n + 1)/2 unknown parameters). As soon as more than n(n + 3)/2 sample points are available, the linear system corresponding to equation 6 is overdetermined: the solution is found by evaluating the pseudo-inverse of this linear system. The cost of the direct computation of this pseudo-inverse by standard numerical methods (e.g. directly available in Matlab) is scaling as n 6 . Note that if the function f is elliptic, the solution of problem 6 is exactly given by ∇f (X 0 ) and H(X 0 ), and the least-square value reaches 0. However, in the general case, denoting by ĝ and Ĥ the solutions of problem 6, the minimum is not 0. Moreover, a measure of the quality of the approximation of the gradient and the Hessian matrix can be given by

N k=1 f (X k ) -f (X 0 ) -(X k -X 0 ) T ĝ - 1 2 (X k -X 0 ) T Ĥ(X k -X 0 ) 2 .
More precisely, a criterion to detect the cases of poor approximations where matrix Ĥ should not be used in the mutation will be based on the following normalized error:

Q = 1 N N k=1 f (X k ) -f (X 0 ) -(X k -X 0 ) T ĝ -1 2 (X k -X 0 ) T Ĥ(X k -X 0 ) f (X k ) -f (X 0 ) -(X k -X 0 ) T ĝ 2 . (7)
This Q factor is invariant under any dilatation or offset of the target function f .

Gathering examples

The approximation method described above requires at least n(n 2 +3)/2 samples (i.e. points of the search space together with their actual fitnesses). Those samples will of course be gathered from the points used during evolution. A trade-off has to be found, however, between the accuracy of the approximation and the computational effort, as drawing exactly n(n 2 + 3)/2 points at every generation would obviously be far too costly. Hence only points that need to be evaluated during the normal course of the algorithm will be used, and the approximation will be based on the most recent visited points. The approximation could nevertheless be computed at every generation, using some sliding window of recent points. But again, this has a cost, and another trade-off has to be made, computing a new approximation only at given times, using the same covariance matrix in between (this is similar to what is done in most pseudo-Newton methods).

Adapting the step-size

So far, we have only discussed the covariance matrix of the Gaussian mutation described by Equation 1. However, even with the optimal covariance matrix, an ES algorithm will never work properly if the step-size is not correctly adapted too. As mentioned in section 2, CMA-ES has a deterministic update rule for its step-size, but performs as well with self-adaptive step-size. Moreover, the transition between Least-Square phases and CMA-phases of the algorithm makes it difficult to mix the different strategies for the step-size. It was hence decided to use the standard self-adaptive rule of ES algorithms for the step-size throughout the LS-CMA-ES algorithm.

The (1, λ)-LS-CMA-ES

The first complete algorithm based on the ideas developed above is the (1, λ)-LS-CMA-ES algorithm, whose pseudo-code is given in Table 1, where g is the generation counter, x (g) the current parent, σ (g) the current step-size, C (g) the current covariance matrix, and c cov a real parameter that will tune the covariance matrix update rule. We shall now discuss the different steps in turn.

A (1, λ)-SA-ES. The basis of the LS-CMA algorithm is a standard self-adaptive Evolution Strategy: Line 4 of Table 1 is the generation of λ offspring using the Gaussian mutation of equation 1 where-the step-size is a log-normal mutation of the step-size of the parent. Lines 5 and 6 are the usual evaluation and selection steps of (1, λ)-ES. Covariance matrix update. Lines 8 and 10 are the standard covariance matrix adaptation of the CMA method [START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF]. However, the computation of the steepest direction p c must be performed even in LS mode to be accurate if mode must be switched. In LS mode, the covariance matrix of the mutation is unchanged (line 9). Quadratic approximation. Every n upd generations (line 11), the approximation of the Hessian is computed from equation 6 (line 12). The sample points have been gathered along evolution (line 7).

If the error on the approximation is below the threshold, then the new approximation replaces the old one, whether it came from CMA or LS modes, and mode is turned to LS (line 13). Otherwise, mode is turned to CMA (line 14). And the main loop resumes (line 16).

Initialization: x (0) ← x0; σ (0) ← σ0; C (0) ← I d ; p (0) c ← 0; archive ← ∅ Mode ← LS 1 while (not stopping condition ) do 2 g ← g + 1 3 Create λ offspring: x (g+1) j ← x (g) + σ (g) exp(τ Ñj(0, 1))Nj (0, C (g) ), j ∈ [1, λ] 4 Evaluate offspring: Compute f (x (g+1) j ) for all 1 ≤ j ≤ λ 5 Select best offspring x (b) : x (g+1) ← x (g+1) b ; σ (g+1) ← σ (g) exp(τ Ñb (0, 1)) 6 Store offspring in archive 7 p (g+1) c ← (1 -cc)p (g) c + √ cc(2-cc) σ (g) (x (g+1) -x (g) ) 8 if in LS mode C (g+1) ← C (g) 9 else (CMA mode) C (g+1) ← (1 -ccov)C (g) + ccovp (g+1) c (p (g+1) c
) T 10 if (g mod n upd =0) 11 Compute ĝx (g) , Ĥx (g) from n2 recent archived samples solving Eq. 6 The parameters. The values of some parameters of the LS-CMA-ES algorithm have to be manually fixed. The following values, either from [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF] or [START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF] when relevant, or adjusted by trial-and-error during preliminary runs, have been used in all experiments of Section 5. The number of offspring λ was set to 10 for the LS-CMA-ES and to 4 + 3 * log(N ) for the CMA algorithm, c c and c cov , the relaxation parameters, to respectively 4 n+4 and that will decide to switch back and forth between LS mode and pure CMA mode, was set to 10 -3 , after intensive monitoring on the elliptic and Rosenbrook functions. Finally, n upd was set to 100 with of course a possibly delay in the first update in order to have n 2 samples available for the first update.

12 if (Q(ĝ x (g) , Ĥx (g) ) < Q th mode ← LS; C (g+1) ← ( 1 2 Ĥx (g) ) -

Experimental Results

This section presents the first results obtained using the (1, λ)-LS-CMA-ES algorithm. It will be compared to the best-performing to-date algorithm in the class of ES-like algorithms (i.e. relying on repeated Gaussian mutations), the original CMA algorithm, as described in [START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF], using the original implementation that the first author kindly makes available [START_REF] Hansen | CMA-ES for Noisy Optimization: Implementations in Matlab[END_REF]. Two variants of CMA are used, the best performing CMA, that uses λ 2 parents as default value, and a (1, λ)-CMA, for a fair comparison, as no experiment has yet been performed on multi-membered LS-CMA-ES. The analytical form of all test function is given in Table 5.

Results: First, the elliptic function f elli and the Rosenbrock function f Ros have been intensively tested. A priori, if our hypothesis of section 3 about the optimality of ( 12 H) -1 as covariance matrix is correct, the LS-CMA-ES method should perform very well on the elliptic function. Rosenbrock function, on the other hand, is quite different, and behaves more like |x| 6 close to the optimum. The approximation error Q defined in Equation 7should hence be larger.

All tests have been performed 100 times independently, from the same initial position (5, 5, . . . , 5). At each generation, the minimum, median and maximum values out of the 100 runs have been gathered, and are shown on Figure 1.

The results for f elli are very clear, and match our expectations: the three groups of three different lines, from left to right are respectively the min, median and max for LS-CMA-ES, then the min, median and max for pure CMA and finally the min, median and max for (1, λ)-CMA. The variance is very small for all three algorithms. For both CMA-ES algorithms, the flat plateau before actually diving toward 0 is the time necessary to learn the correct covariance matrix -while LS-CMA-ES has it all right from the beginning, whatever the dimension (only the minimum number of sample points before the first matrix can be learned makes a small difference). Note, however, that such plateaus can result in a huge difference of actual results in practice, when the stopping criterion is likely to be driven by the computational cost when the actual value of the minimum is generally unknown.

The picture is somewhat different with the Rosenbrock function -which is not unimodal, and has a very flat region between the local optimum and the global optimum. First, only the min and median values have been plotted. For both 20 and 50 dimensions, the first three lines from left to right are the min values of, respectively, LS-CMA-ES, pure CMA and (1, λ)-CMA. Next three lines are the median value for, respectively, the pure CMA, the (1, λ)-CMA and the LS-CMA-ES. As can be seen in dimension 50, the median run of LS-CMA-ES did not converge to the right minimum and the fitness value stayed above 1. Indeed, for all three algorithms, some of the 100 runs did not converge to the global optimum. More precisely, this concerns, for LS-CMA-ES, pure CMA and (1, λ)-CMA, 25 (resp. 12 and 22) runs in dimension 20, and 65 (resp. 15 and 20) runs in dimension 50. Furthermore, a detailed investigation of the approximation error shows that LS is active at the beginning of evolution, far from the flat optimum, and boosts the search there. As the individuals get closer to the minimum, the approximation error increases, and at some point the algorithm switches to CMA, sometimes keeping its advance (e.g. for the runs who gave the min values), sometimes not. Some better tuning of the error threshold might prevent this to happen, but more detailed investigations are needed to understand why CMA catches up.

Several other functions, also used in [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF][START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF] (except for f exp ), have been investigated to further validate those first conclusions. Some off-line results are summarized in table 5. The tendencies observed on the elliptic and Rosenbrock functions are indeed confirmed: f cigar , f tablet and f cigtab are elliptic functions, on which LS-CMA-ES definitely outperforms even the complete CMA; On the other hand, f dif f -pow is and f exp are not elliptic, f exp is even "infinitely flat" (all derivatives are 0 at the minimum), and the results are similar to those on the Rosenbrock function in terms of switching between LS mode and CMA mode. However, as these functions are unimodal, no premature convergence takes place. (or -when it never did within 10 5 evaluations.

f dif f -pow = n 1 |xi| 2+10 i-1 n-1
Discussion: Those experimental results show an outstanding improvement of performance for LS-CMA-ES over pure CMA-ES on elliptic functions. This should not be too surprising if referring to classical methods of numerical optimization: Whereas the CMA-ES can be seen as a steepest descent technique, where the covariance matrix is updated using the previous descent direction (see Section 2), the LS-CMA-ES can be considered similar to second-order quasi-Newton techniques, i.e. relying on local curvature information to find out next points to sample.

Yet, the improvement has only be demonstrated on unimodal or quasi-unimodal functions -and the niche for Evolutionary Algorithms is global optimization in multi-modal context. However, the exploitation phase of EAs is important too, and many refinements of ES algorithms have concentrated on that. As for exploration, it is generally emphasized in the early generations, when the step-sizes are still large enough, by using more than one parent. At the moment, only (1, λ)-LS-CMA-ES have been experimented with, and for few values of λ. It is expected that using multiple parents within LS-CMA-ES will increase its performance, for multi-modal problems, of course, including the Rosenbrock function, but also on non-elliptic unimodal problems, where at the moment if performs better than the corresponding (1, λ)-CMA.

One particular aspect of the LS-CMA-ES that cannot be seen on the plots of Section 5 is the computational cost: as mentioned in Section 4.1, direct solution of the least-square linear system theoretically scales like n 6 . Thanks to Matlab efficient implementation, it was nevertheless possible to go up to dimension 70 (overhead costs for one approximation on a 3.0GHz Pentium roughly take 1, 280 and 1320 seconds for dimensions 20, 50 and 70). A first remark is that in high dimensions, this additional cost could be greatly reduced by using an iterative method to solve Equation 6. Furthermore, when solving very costly real-world problems, e.g. in CFD domain, where the computation of a single fitness value can take up to a few hours, the approximation overhead will be negligible compared to the gain of a few fitness computation. This paper has introduced LS-CMA-ES, the first second-order evolutionary algorithm for parameter optimization. Initial results have demonstrated the efficiency of this algorithm on a large class of (quasi-)unimodal functions, where the best CMA-ES algorithm is outperformed by several orders of magnitude on elliptic problems. Moreover, a dynamic criterion allows the algorithm to switch back to CMA when the second-order approximation is poor, resulting in better results than the one-parent CMA-ES on "more flat" functions, as far as the best runs are concerned. There remains a lot of room for improvement of the proposed algorithm. First, its parameters (e.g. the number of offspring) needs to be finetuned. Next, there might exist some CMA-like update rule for the step-size to tackle premature convergence problems like observed on Rosenbrock function. Finally, the multi-membered version of the algorithm will be investigated, to address multi-modal problems. We do believe that this algorithm is a significant step forward in Evolutionary Parameter Optimization, and will be able to combine the accuracy of the best deterministic methods with the global search capabilities of Evolutionary Algorithms.

Fig. 1 .

 1 Fig. 1. Comparative on-line results for the elliptic function in dimensions 50 and 70 (above) and for the Rosenbrook function in dimensions 20 and 50 (below). See text for line labels.

Table 1 .

 1 Pseudo code for the (1, λ)-LS-CMA-ES algorithm

	1	13

Table 2 .

 2 Off-line results for different test functions in dimension 20. The figures are the number of function evaluations (unit = 10 4 ) before the algorithm reached 10 -10

	1.53 2.14 3.66 1.18 1.36 1.58 1.53 2.14 2.66

(n+ √2) 2[START_REF] Hansen | Reducing the Time Complexity of the Derandomized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES)[END_REF]. The update parameter of the self-adaptive log-normal mutation of the step-size, τ , was set to 1 √ n[START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF]. The threshold Q th on the approximation error for the Hessian matrix