B Kaviraj

O Coupiac

H Courtois

F Lefloch

Noise Correlations in Three-Terminal Diffusive Superconductor-Normal Metal-Superconductor Nanostructures

Keywords: numbers: 74.45.+c, 72.70.+m, 74.40.-n, 74.78.Na

We present measurements of current noise and cross-correlations in three-terminal Superconductor-Normal metal-Superconductor (S-N-S) nanostructures that are potential solid-state entanglers thanks to Andreev reflections at the N-S interfaces. The noise correlation measurements spanned from the regime where electron-electron interactions are relevant to the regime of Incoherent Multiple Andreev Reflection (IMAR). In the latter regime, negative cross-correlations are observed in samples with closely-spaced junctions.

Non-local entanglement -the emblematic ingredient to quantum physics -has been proposed and debated for nearly more than a century [1]. Having already been demonstrated in photons [2], it remains difficult to observe for massive particles like electrons. One interesting approach is to perform noise-correlation measurements using superconductors, as they are natural sources of entangled electrons. In mesoscopic systems, nonequilibrium current noise measurements provide information on the charge and statistics of current-carrying states [3]. For non-interacting electrons, the Pauli exclusion principle dictates the zero-frequency cross-correlations to be negative [4], whereas no such rule applies in the presence of interactions. In early experiments on electronic analogues of teh Hanbury-Brown and Twiss [5] experiment, noise correlations came out to be negative. For beams of partitioned electrons using high mobility GaAs two-dimensional electron gases (2DEG), anti-correlations were revealed, thus illustrating that fermions exclude each other [6]. Fermionic correlations have also been observed for free electrons [7] and neutrons [8]. Interestingly, positive correlations observed in a purely normal 2DEG were ascribed to different scattering mechanisms within the device [9]. Recently, cross-correlations with a bias-dependent sign were reported in a three-terminal superconducting hybrid nanostructure with tunnel contacts [10].

In a superconducting hybrid (N-S) beam splitter made of two normal metal (N) leads in contact with a superconductor (S), an incident electron (hole) from one lead can be reflected at the superconductor interface as a hole (electron) propagating into the other lead. This non-local crossed Andreev reflection (CAR) process corresponds to the transfer (creation) of a superconducting Cooper pair into (from) two entangled electrons in the two leads [11]. It creates positive correlations between currents flowing in each of the two leads [12]. On the contrary, elastic co-tunneling of an electron (hole) from one lead to the other contributes to negative correlations [13]. Whereas negative correlations are expected from a semi-classical approach neglecting proximity effect [14], positive correlations at sub-gap energies are predicted for intermediate values of transparencies [15]. Positive correlations could also arise from synchronized Andreev reflections [16]. In a normal metallic dot connected to all superconducting leads, positive or negative correlations are expected, depending upon properties of contacts between the dot and its leads [17].

In a S-N-S junction, charge transfer is mediated by Multiple Andreev Reflections (MAR), a process during which quasiparticles undergo successive Andreev reflections at both interfaces until their energy reaches the superconducting gap. At a bias above the Thouless energy E T h = hD/L 2 , where D is the diffusion constant and L is the junction length, MAR are not phase-correlated. In this incoherent MAR (IMAR) regime, the noise has been found to be very much enhanced compared to the normal case due to the confinement of the sub-gap electrons in the sandwiched normal metal [18]. In the same regime, cross-correlations are also predicted to be enhanced due to IMAR processes [17]. At low bias voltage, when the time taken by the quasiparticles to reach the superconducting gap exceeds the electron-electron interaction time, the IMAR cycle is interrupted. In this hot electron regime, a Fermi-Dirac like distribution with an elevated effective temperature is then restored [18].

In this Letter, we present measurements of noise and cross-correlations in three-terminal diffusive S-N-S nanostructures. In the IMAR regime, we demonstrate negative correlations in samples where the two junctions are closely spaced.

Samples comprising three-terminal double S-N-S junctions (see SEM pictures in Fig. 1) were fabricated by multiple angle evaporation through a PMMA-PMMA/MAA bilayer mask in an ultra-high vacuum chamber. Evaporation of Cu with a 50 nm thickness was followed immediately by the evaporation of 500 nm-thick Al electrodes, thereby forming diffusive S-N-S junctions across the Cu bridge with highly transparent S-N interfaces. The width of the Cu part was 0.9 µm. The samples have two different distances between junctions: 10 µm for the widely-separated junctions (sample W) and 0.5 µm for the closely-spaced ones (sample C). In the latter case, their distance is comparable with the superconducting Al coherence length, thus making CAR and elastic cotunneling probable. Here, we discuss results obtained on two samples W and C with a junction length L = 1.3 µm. Figure 1 depicts I-V characteristics of individual junctions of the two samples at 100 mK. The junction resistances of a given sample are very much symmetrical, although they differ slightly in the widely-separated geometry. All junctions show a superconducting branch with a small critical current (of the order of few µA), followed by a linear part corresponding to the normal-state resistance of the normal metal, of the order of 1.5 Ω. The related Thouless energy is estimated to be about 5 µeV. The abrupt transition in I-V curves at large currents (≥ 50 µA) is due to depairing effects in the superconducting Al electrodes. The temperature dependence of the I-V characteristics (not shown) does not show much of variation, apart from the depairing current, which decreases at higher temperatures and vanishing critical currents.

We have used a new experimental set-up especially designed to measure current fluctuations and noise crosscorrelations in three-terminal devices at low temperatures. The experiment operates down to 30 mK and is equipped with three commercial SQUIDs (Superconducting Quantum Interference Devices) as sensitive current amplifiers, see Fig. drop across each junction (V 1 , V 2 ) and across the reference resistor (

V ref ). We found that V 1 ≃ V 2 but V ref differs slightly from V 1(2)
. This difference, of about 5 %, is due to an additional resistor between the sample and the reference resistor. For the two samples discussed here, it corresponds to a resistance r ≃ 40 mΩ, which is not negligible compared to R ref and needs to be taken into account in the model discussed below. The intrinsic noise level of each SQUID expressed in equivalent current at its input coil is of few pA/ √ Hz. At frequencies above few hundreds of Hertz, 1/f noise contributions were negligible for all bias currents, so that a frequency domain between 800 Hz and 4 kHz was chosen for the noise and correlation measurements.

In accordance with Nyquist representation, each resistive part of the circuit is associated with a current source δI i . The two junctions of a sample are represented as two non-linear resistors R 1 and R 2 . Each SQUID measures partially each of the 3 current sources [19]. From the 3 fluctuating SQUID currents δI i sq , we can perform 3 auto-correlations AC i ≡ δI i sq δI i sq and 3 cross-correlations XC ij ≡ δI i sq δI j sq non-independent measurements. The spectral densities of noise AC i and correlations XC ij are related to the physical quantities S i ≡ δI i δI i and S ij ≡ δI i δI j through a 6 x 6 matrix. Here, S 1 and S 2 are the noise of each junction of a sample and S 3 the thermal noise of the resistor (R ref + r).

Figure 3 In a matrix form, the equations of the system reduce to:

  S 1 S 2 S 12   = M -1   AC 1 -AC 0 1 AC 2 -AC 0 2 XC 12 -XC 0 12   -N S 3 (1)
where AC 0 1 , AC 0 2 and XC 0 12 represent the background noise and correlation of the set-up. Finding S 1 , S 2 and S 12 thus relies on the determination of the noise and correlation backgrounds, as well as of the thermal noise S 3 . The elements of the 3 x 3 matrix [M -1 ] and of the vector [N ] depend only on the differential resistance of each junction and on the resistance R ref + r.

Let us first discuss the case of sample W, where zero correlations (S 12 = 0) are expected. For a S-N-S junction in the IMAR regime, the shot noise in the diffusive case 2eI/3 is amplified by the number 1 + 2∆/eV of times a quasi-particle is Andreev-reflected before it reaches the gap edge [20]:

S(V ) = 1 3 2eI[1 + 2∆ eV ] ( 2
)
where ∆ is the superconducting gap. Neglecting the contribution of XC achieve a reliable fit. Knowing AC 0 1 and AC 0 2 and choosing XC 0 12 = 0 first, we can plot the correlation S 12 . Our analysis shows that, in order to obtain a bias-independent S 12 , we need to consider the spurious resistance r as being split into two parts r , and r ,, sitting at the mixing chamber with a temperature T and 4.2 K respectively. The best result is obtained for r , = 15 mΩ and r ,, = 25 mΩ, see Fig. 4.

Finally, by choosing XC 0 12 , we get S 12 equal to zero over almost the entire voltage range. It differs from zero only at higher voltages approaching the depairing regime, where the common superconducting electrode becomes resistive and the model is no longer valid. From calibration measurements, we know that, due to the current bias, the electronic temperature of R ref can reach up to 300 mK from a bath temperature of 100 mK. However, this has only a minor effect here since most of the noise S 3 comes from the resistance r ,, sitting at 4.2 K [19].

For sample C with closely-spaced junctions, the measured contact resistance r has the same value as that of sample W. The results can therefore be analyzed exactly the same way as above. In order to adjust the correlation background XC 0 12 value, we consider the low bias regime. Here, the noise of each junction is thermal with an elevated effective electron temperature. The two noise sources δI 1(2) can then be assumed as uncorrelated and hence the cross-correlation noise is zero at low bias: S 12 (V ≃ 0) = 0. We used this criteria to determine XC 0 12 . The overall results of noise and correlations for the two different geometries are depicted in Fig. 5 for different temperatures from 100 to 700 mK. Clearly, the two samples show different cross-correlation behaviors above a voltage of about 40 µV, which corresponds to the cross-over between the hot electron regime and the IMAR regime. For sample C in the IMAR regime, correlations are negative, up to the depairing regime. When the temperature is increased, thermally activated quasiparticles generate additional thermal noise [18,20]. As these fluctuations are uncorrelated, correlations are expected to exhibit a negligible temperature dependence. This contrast in behavior is actually observed in Fig. 5. In addition to the two samples discussed here, other samples with the same geometry but a different junction length of 1 µm were studied and produced very similar results. For these shorter junctions, the noise level is larger since the resistance is lower, see Eq. 2. Again, correlations were found to be zero in the sample with widely separated junctions and negative in that with closely spaced junctions.

Negative correlations are expected for fermionic systems. Therefore, our results suggest that negative correlations arise from partition of quasiparticles injected above the superconducting gap as the result of IMAR processes. When IMAR processes are interrupted by inelastic collisions, the quasiparticle current is reduced and correlations vanish. This is in agreement with our findings. A quasiparticle current is known to vanish over the quasiparticle diffusion length, of the order of few µm in Al [21]. Zero correlations are thus expected in samples with widely separated junctions, as observed. To our knowledge, the role of the quasiparticle current on correlations has never been investigated theoretically. It may restrict the possibility to observe positive correlations in three-terminal devices with all superconducting contacts.

In conclusion, we have measured negative crosscorrelations in three-terminal diffusive S-N-S nanostructures in the incoherent multiple Andreev reflections regime. This experiment opens the way towards a better understanding of non-locality and entanglement in superconducting nano-devices with various interface transparencies and bias schemes. This work was supported by ANR ELEC-EPR project. Samples were fabricated using NanoFab facility at Néel Institute. We thank M. Houzet, J. Meyer and R. Mélin Supplemental Material.

SQUID currents equations.

In accordance with Nyquist representation, each resistive part of the circuit is associated with a current source δI i with i = 1, 2 or 3. The two junctions of a sample are represented as two non-linear resistors R 1 and R 2 , see Fig. 2 of the main text. The current fluctuations through each SQUID can thus be expressed as:

δI i sq = R i R i + R jk δI i - R jk R i + R jk (δI j + δI k ) (3)
where i, j and k ∈ {1, 2, 3} and R jk is the equivalent resistance of R j in parallel to R k . The third resistor R 3 in the model correspond to R ref + r. From the three fluctuating SQUID currents δI i sq , we can perform three autocorrelations AC i ≡ δI i sq δI i sq and three cross-correlations XC ij ≡ δI i sq δI j sq non-independent measurements given by:

AC i =< F F T * (δI i sq )F F T (δI i sq ) > XC ij =< F F T * (δI i sq )F F T (δI j sq ) > (4)
where F F T stands for the F astF ourierT ransf orm, F F T * its complex conjugate and < ... > the rms average at the spectrum analyzer.

The spectral densities of noise {AC i } and crosscorrelations {XC ij } are related to the physical quantities S i ≡ δI i δI i and S ij ≡ δI i δI j through a 6 x 6 matrix:

        AC 1 AC 2 AC 3 XC 12 XC 23 XC 13         = M 6×6         S 1 S 2 S 3 S 12 S 23 S 13         (5)
However, since R ref and r are macroscopic resistors sitting far apart from the sample, we consider that S 13 = S 23 = 0. The SQUID current equations reduce then to a 6x4 matrix. As an example, the 6x4 matrix for

R 1 = R 2 = 1.5Ω and R 3 = R ref + r = 0.1Ω is : M 6×4 =        0.
       (6)
As the six equations relating the six possible measurements, {AC i } and {XC ij }, to the four quantities, {S i } and S 12 , are not independent, the set of equations needs to be reduced to three independent equations. For this, we assume that the spectral density of noise S 3 generated by the resistance R 3 = R ref + r is only thermal.

The SQUID equations reduce then to :

  AC 1 AC 2 XC 12   = M   S 1 S 2 S 12   + N S 3 +   AC 0 1 AC 0 2 XC 0 12   (7)
where we have introduced the noise and correlation backgrounds AC 0 1 , AC 0 2 and XC 0 12 of SQUIDS 1 and 2. The matrix [M ] is now a 3x3 matrix and [N ] a vector, whose elements depend only on the differential resistance of each junction and on the resistance R ref + r.

To obtain the noise spectral densities S 1 and S 2 and the correlations S 12 we need to inverse the new 3x3 matrix. We then obtain :

  S 1 S 2 S 12   = M -1   AC 1 -AC 0 1 AC 2 -AC 0 2 XC 12 -XC 0 12   -N S 3 (8)
Again, as an example, the inverse matrix Effect of the temperature of the spurious resistor r

[M ] -1 for R 1 = R 2 = 1.5Ω and R 3 = R ref + r = 0.1Ω reads : M -1 =   1
As the additional resistor r is due to contact between distinct superconducting wiring elements, it can be split into two parts r , and r ,, sitting at the mixing chamber with a temperature T and 4.2 K respectively. The noise power S 3 is then

S 3 = 4k B [(R ref + r , )T + r ,, × 4.2]/(R ref + r) 2 (10)
To analyze the noise and correlations results of sample W, we proceed in two steps. As a first step, we only consider S 1 and S 2 and suppose that R ref and r are at the bath temperature T (r , = 40mΩ and r ,, = 0 in Eq 10). We adjust AC 0 1 , AC 0 2 and ∆ = 170µeV so that S 1 and S 2 fit the IMAR prediction (Eq. 2 of the main text), see Fig. 4 of the main text. Here, we neglect the contribution of XC 12 to S 1 (2) , which is supported by the facts that XC 12 is much smaller than AC 1 (2) and that the matrix elements coupling XC 12 to S 1(2) are small numbers.

Once we have estimated AC 0 1 , AC 0 2 , we can plot S 12 using equation 8 and choosing XC 0 12 = 0 first. At this stage, we see in Fig. 6 of this Supplemental Material, that the noise correlations S 12 is not constant as a function of the voltage (up triangle black symbols). We can also consider that the resistor r is entirely due to contact at the SQUID having therefore a temperature of 4.2K (r , = 0 and r ,, = 40mΩ in Eq 10). Again the S 12 is not constant (down triangle green symbols). We then adjust r , and r ,, so that the noise correlations S 12 becomes biasindependent. The new expression for S 3 is then fed back in Eq. 8 for iterative and fine adjustment of the parameters values. The best result is obtained for r , = 15 mΩ and r ,, = 25 mΩ. Finally, XC 0 12 is adjusted to get S 12 equal to zero almost over the entire voltage range (square red symbols). Since the contact resistance is found to be the same for sample C, we use the same split of the additional resistor r to analyze the data obtained on that sample with closely-spaced junctions.

FIG. 1 .

 1 FIG. 1. Current-voltage characteristics of each junction J1 and J2 of samples W and C at 100 mK. Insets: SEM pictures of both sample geometries. The distance between junctions is 10 µm for sample W and 0.5 µm for sample C.

2 .

 2 FIG. 2. Schematics of the circuit consisting of the sample, each junction being represented by differential resistors R1 and R2, and the three SQUIDS to measure the current fluctuations. The reference resistor R ref is used to voltage-bias the sample. According to Nyquist representation, each resistance in the model is associated with a current source in parallel. The two additional resistors r , and r ,, (r = r , + r ,, ) are due to contact between distinct superconducting elements. Their temperatures are the bath temperature and 4.2 K respectively. Their associated noise source is not shown for clarity.

FIG. 5 .

 5 FIG.5. Spectral density of noise S1 and cross-correlation S12 of the two samples W and C for various temperatures from 100 to 700 mK. The noise S2 would give very similar results to that of S1. The solid line reproduces the IMAR predictions of Eq. 2, whereas the dashed line enlightens that negative correlations appear in the IMAR regime.

FIG. 6 .

 6 FIG.6. Cross-correlation S12 between the two junctions of the sample W at 100 mK for various choices of the additional contact resistance split between the mixing chamber and the 4.2 K stages on the correlation amplitude S12.

  shows the 6 raw datas {AC i } and {XC ij } as a function of the voltage drop across sample W. The sum FIG. 3. Raw data of ACi's and XCij 's from sample W. The sum represented by the curve along about 20 pA 2 /Hz is constant, as expected, but non-zero due to the set-up background noise.AC 1 + AC 2 + AC 3 + 2(XC 12 + XC 13 + XC 23 ) is constant thanks to current conservation law, but non-zero due to the noise and correlation backgrounds of the experimental set-up. In order to extract the quantities of interest S 1 , S 2 and S 12 , one needs to choose 3 independent measurements. We chose to focus on AC 1 , AC 2 and XC 12 , which, up to the first order, are close to S 1 , S 2 and S 12 .

		60
	/ Hz)	40	Sample W
	2
	i and XC ij (pA	0 20	T=100mK
	AC	-20	AC1	AC2	AC3
			XC12	XC23	XC13	Sum
		-40
		0	20		40	60	80	100
				Voltage (µV)

  12 to S 1(2) [19], we adjusted AC 0 1 , AC 0 2 and ∆ = 170 µeV so that S 1 and S 2 fit this expression. As inelastic collisions are neglected here, Eq. 2 is valid in a voltage range that is limited but enough to

	/ Hz)	30		Sample W
	2
	Noise & Correlations (pA	0 10 20	0	20 T = 100 mK 40	60	80 S1 S2 S12 IMAR1 IMAR2	100
				Voltage (µV)
	FIG. 4. Spectral density of noise S1 and S2 and cross-
	correlation S12 of sample W at 100 mK, obtained using Eq.
	1. The lines correspond to the IMAR noise prediction Eq.
	2 for each individual junction. The dashed region indicates
	the interacting hot electron regime, where Eq. 2 is no longer
	valid.

  for discussions. . 99, 117002 (2007). [19] see supplementary document on line. [20] E. V. Bezuglyi, E. N. Bratus', V. S. Shumeiko, and G. Wendin, Phys. Rev. Lett. 83, 2050 (1999); K. E. Nagaev, Phys. Rev. Lett. 86, 3112 (2001). [21] F. Hübler, J. Camirand Lemyre, D. Beckmann and H. v. Löhneysen, Phys. Rev. B 81, 184524 (2010).

	Lett


Roland Lenain
email: roland.lenain@cemagref.fr

Benoit Thuilot
email: benoit.thuilot@lasmea.univ-bpclermont.fr

Oliver Hach

Philippe Martinet

High-speed mobile robot control in off-road conditions: a multi-model based adaptive approach

This paper is focused on the design of a control strategy for the path tracking of off-road mobile robots acting at high speed. In order to achieve high accuracy in such a context, uncertain and fast dynamics have to be explicitly taken into account. Since these phenomena (grip conditions, delays due to inertial and low-level control properties) are hardly measurable directly, the proposed approach relies on predictive and observer-based adaptive control techniques. In particular, the adaptive part is based on an observer loop, taking advantage of both kinematic and dynamic vehicle models. This multimodel based adaptive approach permits to adapt on-line the grip conditions (represented by cornering stiffnesses), enabling highly reactive sideslip angles observation and then accurate path tracking. The relevance of this approach is investigated through full scale experiments.

I. INTRODUCTION

Off-road mobile robots are a subject of increasing popularity, since their potential applications may meet social needs, especially concerning environmental issues. From surveillance to agricultural operations, the use of robots in natural environment may permit to increase the work efficiency, while reducing pollution as well as the risk for human personnel [START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF]. In order to make such applications feasible, autonomous path tracking operations must demonstrate high accuracy, even at high speed (see for instance [START_REF] Schäfer | Ravon -a robust autonomous vehicle for off-road navigation[END_REF]). This implies the design of advanced control strategies, since classical control approaches [START_REF] Micaelli | Trajectory tracking for unicycle-type and two-steering-wheels mobile robots[END_REF] relying on the assumption of rolling without sliding are no longer valid in the off-road context, especially at high speed. Robust approaches have been investigated by considering sliding as a perturbation (see for instance [START_REF] Eaton | Precision guidance of agricultural tractors for autonomous farming[END_REF] or [START_REF] Wang | An analysis of wheeled mobile robots in the presence of skidding and slipping: Control design perspective[END_REF]). Nevertheless, the efficiency of such approaches is limited, since they tend to be conservative, leading to an oscillating behavior of the robots even at limited speed.

Consequently, model-based control seems to be more attractive. However, simple kinematic approaches cannot reflect the entire robot behavior, while complete dynamic descriptions (as achieved in [START_REF] Ellouze | Control of unicycle-type robots in the presence of sliding effects with only absolute longitudinal and yaw velocities measurement[END_REF]) require the knowledge of numerous parameters. For example the ground-tire interaction model described in [START_REF] Bakker | Tyre modeling for use in vehicle dynamics studies[END_REF] employs 14 parameters, which are subject to changes in natural environment (as grip conditions are not constant). In [START_REF] Lenain | Multi-model based sideslip angle observer: Accurate control of high-speed mobile robots in off-road conditions[END_REF], an adaptive approach has been proposed, based on both extended kinematic and dynamic models. The control law is designed from the extended kinematic model, consisting in a Ackermann model that additionally accounts for the sideslip angles. Then, based on dynamical equations of motion, an observer is designed to estimate the grip conditions and the sideslip angles needed in the control law expression. However, the matrix inversion required in grip condition estimation becomes singular when the robot does not laterally slide (for instance, in straight line motion). In this case, the estimation process must be stopped. Most of the time, there is no serious consequence, since the mobile robot does not slide. Nevertheless, such switches can generate delay and transient lacks of accuracy, and exceptionally instabilities can occur.

In this paper, a new strategy for grip condition estimation is proposed, relying on a sensitivity-based gradient search algorithm, and its integration into the adaptive control method is shown. In contrast to the method presented in [START_REF] Lenain | Multi-model based sideslip angle observer: Accurate control of high-speed mobile robots in off-road conditions[END_REF], it allows a continuous estimation without any matrix inversion. Section II of this article pictures the modelling of an off-road mobile robot from kinematic and from dynamic points of view. Next, the control law, based on the extended kinematic model, is presented in Section III. The estimation of the sideslip angles, needed in the control law, is then described in Section IV. The estimation process is composed of several steps: first, a preliminary estimation of the sideslip angles relying on the extended kinematic model (slow-varying estimation), then the estimation of tire-soil interaction forces and the adaptation of the grip conditions (represented by cornering stiffnesses) and finally the observation of the sideslip angles based on the dynamic model (enabling a high reactive observation). In section V, results of experiments in off-road conditions show the efficiency of the approach in terms of path tracking accuracy at relatively high speed (up to 6 m s -1 ).

II. MOBILE ROBOT MODELLING

A. Scheme and notation

As the work described in this paper is dedicated to off-road mobile robot control, the classical "rolling without sliding" assumption can no longer be applied (it would lead to large tracking errors). Thus, a purely kinematic model cannot be used in such a context. Two alternative representations are therefore defined below, using the common framework depicted in Fig. 1.

Fig. 1. Path tracking parameters

The mobile robot is considered as a bicycle, whose motion is controlled with respect to a reference path Γ. The common modelling variables and parameters are listed below:

• L = a + b is the vehicle wheel base, where a and b are the distances of the front and rear wheel from the center of gravity G, • F F and F R are the lateral forces generated at the front and rear tires, respectively, • m and I z are dynamic parameters (mass and inertia along the vertical axis), • v is the vehicle linear velocity at the center of the rear axle (point to be controlled), • δ F is the front steering angle, it constitutes the control variable, • β F and β R are the front and rear side slip angles, while β denotes the global sideslip angle, • c(s) is the curvature of the path Γ, parametrized by its arc-length coordinate s, • θ is the vehicle orientation w.r.t. an absolute frame • θ is the angular deviation of the vehicle w.r.t. Γ, • y is the vehicle lateral deviation w.r.t. Γ.

B. Dynamic model

In order to take into account skidding effects and other dynamic phenomena, the natural approach consists in considering a dynamical model of motion (i.e. linking forces and accelerations). When considering only the lateral behavior of the mobile robot, the following dynamical model is commonly proposed in the literature (see [START_REF] Ellouze | Control of unicycle-type robots in the presence of sliding effects with only absolute longitudinal and yaw velocities measurement[END_REF]):

                   θ = 1 Iz (-aF F cos(δ F ) + bF R ) β = -1 v2m (F F cos(β -δ F ) + F R cos(β)) - θ β R = arctan tan β - b θ v2 cos(β) β F = arctan tan β + a θ v2 cos(β) -δ F v 2 = v cos(βR) cos(β) (1)
If lateral forces F F and F R are known, then system (1) accurately describes the robot dynamics. However, it does not rely directly on the controlled variable (lateral deviation y), and the required integrations make this model hardly tractable for control purposes.

C. Extended kinematic model

In order to design a relevant control law for high speed path tracking in off-road context, an alternative model is designed, comprising the lateral deviation as well as sliding effects. As introduced in [START_REF]Sideslip angles observers for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF] and generalized in [START_REF] Wang | An analysis of wheeled mobile robots in the presence of skidding and slipping: Control design perspective[END_REF], it is based on a kinematic description of the bicycle motion, to which the sideslip angles that appear in the dynamic model [START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF] are added. The motion of the robot can then be defined by the system of equations [START_REF] Schäfer | Ravon -a robust autonomous vehicle for off-road navigation[END_REF], which constitutes the extended kinematic model.

       ṡ = v cos( θ+βR) 1-c(s) y ẏ = v sin( θ + β R ) θ = v [cos(β R )λ 1 -λ 2 ] (2) with: λ 1 = tan(δF +βF )-tan(βR) L , λ 2 = c(s) cos( θ+βR)
1-c(s) y s, y and θ are assumed to be measured by appropriate sensors (an RTK GPS receiver and a gyrometer). If β F and β R are also available (sliding observation is described in section IV), then model (2) can describe accurately the mobile robot motion, including sliding effects.

III. CONTROL LAW

The structure of model (2) remains close to a purely kinematic description, as reported in [START_REF] Samson | Control of chained systems. Application to path following and time-varying point stabilization of mobile robots[END_REF]. In particular, it preserves the exact linearization properties of these descriptions. This attractive feature can then be used to derive a nonlinear control law, considering the sliding variables β F and β R as known: as described in [START_REF] Cariou | Automatic guidance of a four-wheel-steering mobile robot for accurate field operations[END_REF], once the exact linearization is achieved via nonlinear state and control transformations, a classical PD control can be proposed for the auxiliary inputs in order to ensure the convergence of the actual lateral deviation to zero. The reverse nonlinear transformations provide finally the nonlinear expression (3) for the steering control law.

δF = arctan tan(βR) + L cos(β R ) c(s) cos θ1 α + A cos 3 θ1 α 2 -βF
(3) with K p and K d strictly positive tunable gains and:

   θ1 = θ + β R α = 1 -c(s)y A = -K p y -K d α tan θ1 + c(s)α tan 2 θ1 (4)
In practice, steering actuators introduce some delay, which may damage the path tracking performances of control law (3). Consequently, Model Predictive Control techniques are here applied to address specifically curvature servoing. More precisely, the steering control law is split into two additive terms:

δ F = δ T raj + δ Deviation ( 5
)
δ Deviation is a term mainly concerned with control errors and sliding compensation, while δ T raj deals with the reference path's shape: it imposes that the robot's orientation follows the curvature of the path. As the future curvature of the path to be followed is known, as well as the steering actuator features, a model predictive algorithm can be derived: the value of δ T raj (called δ P red T raj in the sequel) is computed regarding an optimality criterion, in order to reach "at best" the future curvature on a fixed horizon of prediction. Substituting δ T raj by this term, leads to the adaptive and predictive control law:

δ F = δ P red T raj + δ Deviation (6)
IV. SIDESLIP ANGLE OBSERVATION Control law ( 6) requires indeed the knowledge of the sideslip angles β F and β R . As they are hardly measurable directly, their indirect estimation is addressed here. It is based on a algorithm, which relies on both kinematic and dynamic vehicle models, thus enabling accurate control even at high speed. The different steps of the hierarchical estimation process are detailed below.

A. Sideslip angle estimation from extended kinematic model

When considering kinematic model ( 2) from an observation point of view, it can be noticed that it is composed of two equations (the ṡ equation is indeed of no use, since only lateral motion is here considered), when only two unknowns have to be estimated (β F , β R ). Such an observation problem can be solved as described in detail in [START_REF] Cariou | Automatic guidance of a four-wheel-steering mobile robot for accurate field operations[END_REF]. The extended kinematic model (2) can be rewritten as:

Ẋobs = f (X obs , δ F , u) (7) =    v sin( θobs + u R ) v cos(uR)[tan(δF +uF )-tan(uR)] L -c(s) cos( θobs +uR) 1-c(s)y obs
where X obs = [y obs θobs ] T is the observed state vector and u

= [u F u R ] T = [β F β R ]
T are the sideslip angles to be estimated. The general principle is here to consider the sideslip angles u as control variables to be designed in order to achieve the convergence of X obs to the measured state X mes = [y mes θmes ] T . As sideslip angles do not exceed a few degrees in practice, this state equation can be linearized with respect to the control vector u in the vicinity of zero (i.e. no sliding), leading to

Ẋobs = f (X obs , δ F , 0) + B(X obs , δ F )u (8)
with B(., .) denoting the derivative of f with respect to u, evaluated at u = (0, 0).

Provided that θobs = π 2 [π] and v = 0, the matrix B is invertible. Let e = X obs -X mes be the observed error. Then, an exponential convergence ė = G e can be obtained by choosing:

u = B(X obs , δ F ) -1 G e + Ẋmes -f (X obs , δ F , 0) (9)
G has to be chosen as a Hurwitz matrix and constitutes the observer gain, assigning the settling times for the state observation. Since convergence of the observed state X obs to the measured one X mes is achieved, u can be regarded as a relevant estimate of the sideslip angles. This first observation step supplies a satisfactory sliding estimate in off-road context at limited speed: the results reported for instance in [START_REF] Cariou | Automatic guidance of a four-wheel-steering mobile robot for accurate field operations[END_REF] show very accurate results up to 4 m s -1 . Nevertheless, since this observer is designed from a kinematic model, the observation of dynamical variables such as sideslip angles necessarily presents some delays, that may generate oscillations when the robot's velocity is increased.

B. Dynamic sideslip angle estimation

In order to improve the reactivity of the sideslip angle estimation, the design of an observer based on dynamic model ( 1) is preferable. In order to use such a model, an analytical expression of the two forces F F and F R is required. A linear expression, quite standard in the literature (see for instance [START_REF] Ellouze | Control of unicycle-type robots in the presence of sliding effects with only absolute longitudinal and yaw velocities measurement[END_REF]), is here considered:

F F = C F β F F R = C R β R ( 10
)
where C F and C R are the cornering stiffnesses. While such a linear expression can be valid in on-road context, the motion in natural environment (moreover at high speed) implies non-linear tire-ground interaction forces. Nevertheless, as demonstrated in [START_REF] Lenain | Multi-model based sideslip angle observer: Accurate control of high-speed mobile robots in off-road conditions[END_REF], the linear model [START_REF] Samson | Control of chained systems. Application to path following and time-varying point stabilization of mobile robots[END_REF] can still be used, provided that cornering stiffnesses are on-line adapted, in order to be representative of the non-linearity and the variability of the grip conditions. On the condition that C F and C R can be properly online adapted, the integration of the tire forces [START_REF] Samson | Control of chained systems. Application to path following and time-varying point stabilization of mobile robots[END_REF] into the dynamic model (1) leads to the following state space representation (the third and fourth equations in [START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF] have been used to express β F and β R as functions of β):

Ẋ2 = A 2 X 2 + B 2 δ F (11)
where:

A 2 = -a 2 CF -b 2 CR v2Iz -aCF +bCR Iz -aCF -bCR v 2 2 m -1 -CF +CR v2m , B 2 = aCF Iz CF v2m (12) X 2 = [ θ β] T is the state vector, δ F is the control variable.
Since the steering anlge is limited to 20 • , it has been assumed that cos δ F ≈ 1. This allows to simplify the matrix A 2 and turns the dynamic model (1) into the linear form [START_REF] Astrom | Adaptive control (2nd edition)[END_REF].

It has also been assumed that the first state variable θ can be measured (for instance, with a gyrometer). An estimation of the second state variable β can be inferred from [START_REF]Sideslip angles observers for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF]. Let us denote this estimation as β and let us consider β as a virtual measurement of β. Then, a standard observer equation associated with model ( 11) is:

Ẋ2 = A 2 X2 + B 2 δ F + G 2 X2 (13)
where X2 = [ θ2 β2 ] T is the observed state, X = [ θ β] T is the measured state and X2 = X2 -X is the observer error.

From ( 11) and ( 13), it can be deduced that:

Ẋ2 = (A 2 + G 2 ) X2 (14)
The convergence of the observer error X2 to zero is then ensured, provided that G 2 is chosen such that A 2 + G 2 is Hurwitz. β has been chosen as the measurement associated with the observed state β2 , since its steady state value is always correct. However, since the values of β are not necessarily accurate during transient phases, the preference is given to the convergence of θ2 (since yaw rate measurement is reliable) with respect to the convergence of β2 . This can easily be achieved by tuning G 2 such that the settling time associated with θ2 is shorter than the one associated with β2 .

Finally, the front and rear sideslip angles to be injected into control law ( 6) can be obtained by reporting β2 into the third and fourth equations of (1):

   βR = β2 -b θ2 v2 βF = β2 + a θ2 v2 -δ F (15)
Equations ( 15) can be interpreted as a mixed kinematic and dynamic sideslip angle observer. As demonstrated in Section V, observer (15) improves the robustness and the reactivity of the sideslip angle estimation at high speed with respect to observer [START_REF]Sideslip angles observers for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF], and thus the performances of path tracking control.

C. Tire force observation and cornering stiffness adaptation

The previous observer indeed relies on relevantly adapted cornering stiffnesses C F and C R . When considering the dynamical model (1), together with the expression [START_REF] Samson | Control of chained systems. Application to path following and time-varying point stabilization of mobile robots[END_REF] for tire forces, the adaptation of the cornering stiffnesses can be achieved in the same way as for observer [START_REF]Sideslip angles observers for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF]. Such an approach (named "Direct Adaptation" in the sequel) has been presented in [START_REF] Lenain | Multi-model based sideslip angle observer: Accurate control of high-speed mobile robots in off-road conditions[END_REF]. However, it requires the inversion of a matrix, which becomes singular in the non-sliding case (i.e. when β = 0). This situation is likely to occur, even in off-road context (for instance, when following straight line path segments), and in such cases the adaptation must be stopped. Most of the time, these transient inaccuracies in cornering stiffness estimation have no serious consequence, since sideslip angles, and therefore lateral forces, are small. However, this switching behaviour can nevertheless generate delay and transient lacks of accuracy when the mobile robot leaves the straight line and enters into a curve, and, at worst, instabilities can occur.

This problem is addressed here by using the MIT Rule (as presented in [START_REF] Astrom | Adaptive control (2nd edition)[END_REF]), considering the cornering stiffnesses as slow-varying parameters of force model [START_REF] Samson | Control of chained systems. Application to path following and time-varying point stabilization of mobile robots[END_REF]. The first step is the estimation of the forces F F and F R . It is achieved by rewriting model (1) as follows (considering that the values of β are small):

Ẋc = A c Xc + B c (δ F ) [ FF FR ] T (16)
where:

A c = 0 0 -1 0 , B c (δ F ) = -a cos δF Iz b Iz -cos δF v2m -1 v2m
(17) Xc = [ θc βc ] T is the observed state, [ FF FR ] T are the forces to be estimated.

The vector X, introduced in section IV.B, can again be seen as a measurement of Xc . Then, the tire forces can be estimated by imposing the convergence of Xc to X. The matrix B c (δ F ) being invertible as soon as the velocity v 2 is non-zero, this convergence can be obtained by:

[ FF FR ] T = B c (δ F ) -1 G c Xc -A c Xc , ( 18
)
where Xc = Xc -X and G c is a gain matrix. The expression (18) indeed leads to the following error dynamics, provided that Ẋ is negligible:

Ẋc = G c Xc (19)
The estimates Fi i∈{F,R} can then be considered as the actual lateral tire forces applied to the robot. Consequently, in view of [START_REF] Samson | Control of chained systems. Application to path following and time-varying point stabilization of mobile robots[END_REF], and taking into account that cornering stiffnesses are slow-varying with respect to other dynamics, a relevant estimation Ĉi has to satisfy:

Fi = Ĉi u i ( 20
)
where u i are the sideslip angles estimated by [START_REF]Sideslip angles observers for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF]. Relying on a sensitivity-based gradient search algorithm, a cornering stiffness adaptation law is then:

Ċi = -γ( Fi -Ĉi u i ) ∂( Fi -Ĉi u i ) ∂ Ĉi = -γ( Fi -Ĉi u i )u i (21)
where γ is a strictly positive tunable gain.

Adaptation law (21) (named "MIT Adaptation" in the sequel) thus permits to adapt the cornering stiffnesses Ĉi i∈{F,R} , so that model ( 10) is consistent with both the estimation of the tire forces and the kinematic estimation of the sideslip angles. The dynamic observer (15) can then be run and supply the reactive estimates of the sideslip angles needed in control law [START_REF] Ellouze | Control of unicycle-type robots in the presence of sliding effects with only absolute longitudinal and yaw velocities measurement[END_REF].

D. Summary of the observer approach

The hierarchical observer proposed in this paper and composed of the previously described steps can be summarized by the scheme in Fig. 2.

The first step, shown as a red/dashed box in Fig. 2, consists in the preliminary observer [START_REF]Sideslip angles observers for vehicle guidance in sliding conditions: application to agricultural path tracking tasks[END_REF], derived from the extended kinematic model [START_REF] Schäfer | Ravon -a robust autonomous vehicle for off-road navigation[END_REF]. It supplies a first estimation of the sideslip angles, considered accurate enough to estimate slow-varying parameters. This estimation, considered as a virtual measurement, is then used together with the measured yaw rate in order to estimate lateral forces: it consists in observer (18), derived from dynamical model [START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF]. In turn, from these slow-varying estimates of the forces and of the sideslip angles, the cornering stiffnesses (also considered as slow-varying) are adapted via the expression (21), derived Fig. 2. Observer principle scheme also from dynamical model [START_REF] Siegwart | Introduction to Autonomous Mobile Robots[END_REF]. This second step is shown as green/plain boxes in Fig. 2. All the parameters of dynamical model ( 1) are then known (mass and inertial parameters are invariant, and cornering stiffnesses are supplied by expression ( 21)), so that a second sideslip angle observer, more reactive since it takes explicitly into account dynamics effects, can be derived: it consists in equations ( 13) and ( 15) and is shown as a blue/dashed box in Fig. 2. Finally, this last estimation of the sideslip angles is entered into the adaptive and predictive control law [START_REF] Ellouze | Control of unicycle-type robots in the presence of sliding effects with only absolute longitudinal and yaw velocities measurement[END_REF], in order to achieve accurate path tracking.

V. EXPERIMENTAL RESULTS

A. Experimental setup

The electric off-road vehicle depicted in Fig. 3 is used as an experimental platform. Designed for all-terrain mobility, it can climb slopes up to 45 • , its maximum speed is 8 m s -1 , and it has the following properties: The main exteroceptive sensor on-board is a "Magellan ProFlex 500" RTK-GPS receiver, which supplies absolute position measurements with an accuracy of 2 cm at a 10 Hz sampling frequency. The GPS antenna is located vertically above the center of the rear axle, so that the absolute position of the point to be controlled is straightforwardly obtained from the sensor. In addition, a MTi IMU manufactured by Xsens is on-boarded. In this framework, only the vertical axis of the gyrometer, supplying the yaw rate is used. It supplies this measure with an accuracy of 0.1 • /s and is fixed on the chassis.

B. Path tracking results

The results presented in this paper are illustrated by an extended video, available at ftp://ftp.clermont.cemagref.fr/pub/Tscf/Lenain/VideoIcra2011/ Fig. 3. Experimental platform 1) Relevance of cornering stiffness adaptation: The first experiment for the algorithm's validation consists in following the path depicted in Fig. 4. This path has been recorded beforehand, when the robot was steered manually at 1 m s -1 . It is composed of two circles: one performed on asphalt and the other one performed on dry grass. This trajectory has been followed twice, at 4 m s -1 , using control law (6) fed by observer (15):

• In the first run, the cornering stiffnesses are estimated via the "Direct Adaptation" method described in [START_REF] Lenain | Multi-model based sideslip angle observer: Accurate control of high-speed mobile robots in off-road conditions[END_REF].

The results are reported in red dotted line. Switches in the adaptation of the cornering stiffnesses, due to a matrix singularity (see discussion in Section IV.C), are encountered. • In the second run, the cornering stiffness adaptation is ensured by the "MIT Adaptation" approach (i.e. according to (21)). The results are reported in magenta solid line. Fig. 5 shows a comparison plot of the resulting tracking errors. It can first be noticed that both adaptation strategies allow accurate tracking: after initialization (i.e. after 5 s) the lateral deviation does not exceed 50 cm, while large oscillations with a constant deviation of more than 0.5 m during the circles are recorded when using a classical control law (i.e neglecting sliding, i.e. β F and β R are set to 0). Nevertheless, some significant differences can be observed. First, it can be noticed that transient deviations in the tracking error are satisfactorily less important when the "MIT Adaptation" is used. The main difference can be seen during the second circle (from 25 s to 35 s), as the error with the "MIT Adaptation" remains very close to zero, while the error with "Direct Adaptation" oscillates and converges to -25 cm. The explanation for this difference can be found by inspecting the recorded cornering stiffnesses. For the sake of clarity, only front cornering stiffnesses C F are compared in Fig. 6. Obviously, the estimation obtained with the "MIT Adaptation" is smoother than with the "Direct Adaptation", and moreover no delay is added. It thus permits, relying on (15), a smooth sideslip angle estimation, leading to a more stable behavior in path tracking control. Indeed, the major improvement of the proposed observer with respect to "Direct Adaptation" is its continuous operation: even during straight line path segments (from 5 s to 10 s and from 17.5 s to 25 s), the "MIT Adaptation" leads naturally to constant cornering stiffnesses (the estimated forces Fi and the sideslip angles u i are zero, leading to Ċi = 0 in equation ( 21)). On the contrary, such situations lead to singularities in the "Direct Adaptation", so that the adaptation process has to be stopped (before 12.5 s, from 19 s to 27 s and after 29 s). The cornering stiffnesses are then no longer adapted and stay at an arbitrary value, not necessarily representative of the actual grip conditions. As a result, an inaccurate estimation of the sideslip angles may be sent to the control law, explaining the large tracking errors recorded during the second circle.

This fact is highlighted in Fig. 7 that depicts a comparison of the front sideslip angles estimated via "Direct Adaptation" and via "MIT Adaptation", respectively. Differences in the cornering stiffness adaptation have indeed a direct impact on the sideslip angle estimation and consequently on the lateral deviation. They can be observed during both circular segments of the trajectory. During the first circle (10 s to 18 s), the cornering stiffness adapted with "Direct Adaptation" starts varying after the one adapted with "MIT Adaptation", leading to an overestimation of the front sideslip angles (especially at the beginning of the circle at 12.5 s where a difference of 5 • is recorded). During the second circle (from 25 s to 35 s), the "Direct Adaptation" method is interrupted because of a matrix singularity and the adapted cornering stiffness stays equal to 12500 N rad -1 . As a result, the corresponding estimation of the front sideslip angle is overestimated again (around 5 • instead of 2.5 • for the "MIT Adaptation"). These differences in sideslip angle estimation also explain the disparities in the tracking error depicted in Fig. 5 and already commented. At worst, in other situations, the "Direct Adaptation" algorithm may even lead to instability, as it generates discontinuities in the adapted cornering stiffnesses (e.g. in Fig. 6 at 27 s). The proposed "MIT Adaptation" process is thereby more robust, due to its continuous operation.

2) High-speed validation: The above results permit to investigate path tracking at higher speed, close to the controllability limits of the robot. A second reference trajectory, still manually recorded at 1 m s -1 , has been defined on the field shown in Fig. 3 and in the attached video. This trajectory is composed of a straight line, then a large circle and finally a quick inversion of the curvature (see Fig. 8).

Path tracking has then been achieved at 6 m s -1 , using control law [START_REF] Ellouze | Control of unicycle-type robots in the presence of sliding effects with only absolute longitudinal and yaw velocities measurement[END_REF], observer (15) and "MIT Adaptation". The resulting path is displayed as dashed line in Fig. 8, while Fig. 9 shows the recorded tracking error.

It can be noticed that the maximal tracking error of 1 m is recorded in the beginning of the curve (around 12.5 s), Therefore, in this experiment, in view of the actuator capabilities and the grip conditions, the mobile robot was close to its controllability limits, but nevertheless achieves successfully the path tracking mission. For slightly higher speed, this trajectory turns out to be non-admissible for this robot: at 7 m s -1 , the robot starts spinning and is thus no longer able to follow the reference path.

VI. CONCLUSION AND FUTURE WORK

This paper proposes a predictive and adaptive control law, dedicated to off-road path tracking of fast mobile robots. It is based on an observer scheme, mixing kinematic and dynamic models of the considered mobile robot. The interaction between several levels of modelling allows the on-line estimation of the unknown dynamic model parameters, and finally the achievement of a reactive and relevant estimation of the sideslip angles. The observed variables are subsequently used in an associated control law, permitting highly accurate path tracking with respect to the considered speeds and grip conditions variabilities.

However, at high speed, reference paths may become nonadmissible, depending on the grip conditions. Current work is then focused on controlling the robot speed, in order to ensure that the reference path is always achievable, with respect to robot capabilities and grip conditions. This will permit to ensure a minimal error and to avoid spinning situations.

  Total mass m = 350 kg Yaw inertia Iz = 270 kg m 2 Wheelbase L = 1.2 m Rear half-wheelbase b = 0.58 m

Fig. 4 .

 4 Fig. 4. Reference trajectory manually recorded

Fig. 5 .

 5 Fig. 5. Comparison of the tracking error at 4 m s -1

Fig. 6 .

 6 Fig. 6. Comparison of the adapted cornering stiffness

Fig. 7 .

 7 Fig. 7. Comparison of the estimated front sideslip angles

Fig. 8 .

 8 Fig. 8. Tracked pathwhile the lateral deviation afterwards converges to zero during the circle. A second transient overshoot occurs at 27 s, corresponding to the small right curve (implying a fast curvature inversion). The large transient error is the result of two cumulative delays: first, the settling time of the steering actuator is 0.4 s, meaning a settling distance of 2.4 m at the considered speed. Secondly, due to harsh grip conditions, variations in cornering stiffnesses are also reflected after a short delay: it can indeed be noticed, when examining the recorded cornering stiffnesses in Fig.10, that at the beginning of the circle (around 12 s) both ĈF and ĈR , need 2 s to converge to nearly constant values, namely 5000 N rad -1 and 4000 N rad -1 . This corresponds exactly to the duration of the transient deviation (from 12 s to 15 s).

Fig. 9 .

 9 Fig. 9. Tracking error

TABLE I EXPERIMENTAL

 I ROBOT DYNAMIC PARAMETERS

This work is supported by French National Research Agency (ANR), under the grant ANR-07-ROBO-0008 attributed to the FAST project, see https://projetfast.cemagref.fr/


Michel Broniatowski
email: michel.broniatowski@upmc.fr

Minimum divergence estimators, maximum likelihood and exponential families

Keywords: statistical divergence, minimum divergence estimator, maximum likelihood, exponential family

In this note we prove the dual representation formula of the divergence between two distributions in a parametric model. Resulting estimators for the divergence as for the parameter are derived. These estimators do not make use of any grouping nor smoothing. It is proved that all differentiable divergences induce the same estimator of the parameter on any regular exponential family, which is nothing else but the MLE.

1 Introduction

Context and scope of this note

This note presents a short proof of the duality formula for ϕ-divergences defined through differentiable convex functions ϕ in parametric models and discusses some unexpected phenomenon in the context of exponential families. First versions of this formula appear in [START_REF] Liese | Convex statistical distances[END_REF] p 33, in [START_REF] Broniatowski | Estimation of the Kullback-Leibler divergence[END_REF] in the context of the Kullback-Leibler divergence and in [START_REF] Keziou | Dual representation of ϕ-divergences and applications[END_REF] in a general form. The paper [START_REF] Broniatowski | An estimation method for the Neyman chi-square divergence with application to test of hypotheses[END_REF] introduces this form in the context of minimal χ 2 -estimation; a global approach to this formulation is presented in Broniatowski and Kéziou (2006) [START_REF] Broniatowski | Minimization of ϕ-divergences on sets of signed measures[END_REF]. Independently [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF] [START_REF] Liese | On divergences and informations in statistics and information theory[END_REF] have obtained a similar expression based on a much simpler argument as presented in all the above mentioned papers (formula (118) in their paper); however the proof of their result is merely sketched and we have found it useful to present a complete treatment of this interesting result in the parametric setting, in contrast with the aforementioned approaches.

The main interest of the resulting expression is that it leads to a wide variety of estimators, by a plug in method of the empirical measure evaluated on the current data set; so, for any type of sampling its estimators and inference procedures, for any ϕ-divergence criterion. In the case of the simple i.i.d. sampling resulting properties of those estimators and subsequent inferential procedures are studied in [START_REF] Broniatowski | Parametric estimation and tests through divergences and the duality technique[END_REF].

A striking fact is that all minimum divergence estimators defined through this dual formula coincide with the MLE in exponential families. They henceforth enjoy strong optimality under the standard exponential models, leading to estimators different from the MLE under different models than the exponential one. Also this result proves that MLE 's of parameters of exponential families are strongly motivated by being generated by the whole continuum of ϕ-divergences.

This note results from joint cooperation with late Igor Vajda.

Notation

Let P := {P θ , θ ∈ Θ} an identifiable parametric model on R d where Θ is a subset of R s . All measures in P will be assumed to be measure equivalent sharing therefore the same support. The parameter space Θ need not be open in the present setting. It may even happen that the model includes measures which would not be probability distributions; cases of interest cover models including mixtures of probability distributions; see [START_REF] Broniatowski | Parametric estimation and tests through divergences and the duality technique[END_REF]. Let ϕ be a proper closed convex function from ] -∞, +∞[ to [0, +∞] with ϕ(1) = 0 and such that its domain domϕ := {x ∈ R such that ϕ(x) < ∞} is an interval with endpoints a ϕ < 1 < b ϕ (which may be finite or infinite). For two measures P α and P θ in P the ϕ-divergence between Q and P is defined by

φ(α, θ) := X ϕ dP α dP θ (x) dP θ (x).
In a broader context, the ϕ-divergences were introduced by [START_REF] Csiszár | Eine informationstheoretische Ungleichung und ihre Anwendung auf den Beweis der Ergodizität von Markoffschen Ketten[END_REF] as "fdivergences". The basic property of ϕ-divergences states that when ϕ is strictly convex on a neighborhood of x = 1, then φ(α, θ) = 0 if and only if α = θ.

We refer to [START_REF] Liese | Convex statistical distances[END_REF] chapter 1 for a complete study of those properties. Let us simply quote that in general φ(α, θ) and φ(, θ, α)are not equal. Hence, ϕ-divergences usually are not distances, but they merely measure some difference between two measures. A main feature of divergences between distributions of random variables X and Y is the invariance property with respect to common smooth change of variables.

Examples of ϕ-divergences

The Kullback-Leibler (KL), modified Kullback-Leibler (KL m ), χ 2 , modified χ 2 (χ 2 m ), Hellinger (H), and L 1 divergences are respectively associated to the convex functions ϕ

(x) = x log x-x+1, ϕ(x) = -log x+x-1, ϕ(x) = 1 2 (x -1) 2 , ϕ(x) = 1 2 (x -1) 2 /x, ϕ(x) = 2( √ x -1)
2 and ϕ(x) = |x -1|. All these divergences except the L 1 one, belong to the class of the so called "power divergences" introduced in [START_REF] Read | Goodness-of-fit statistics for discrete multivariate data[END_REF] (see also [START_REF] Liese | Convex statistical distances[END_REF] chapter 2), a class which takes its origin from Rényi [START_REF] Rényi | On measures of entropy and information[END_REF]. They are defined through the class of convex functions

x ∈]0, +∞[ → ϕ γ (x) := x γ -γx + γ -1 γ(γ -1) (1) if γ ∈ R \ {0, 1}, ϕ 0 (x) := -log x + x -1 and ϕ 1 (x) := x log x -x + 1.
So, the KL-divergence is associated to ϕ 1 , the KL m to ϕ 0 , the χ 2 to ϕ 2 , the χ 2 m to ϕ -1 and the Hellinger distance to ϕ 1/2 .

It may be convenient to extend the definition of the power divergences in such a way that φ(α, θ) may be defined (possibly infinite) even when P α or P θ is not a probability measure. This is achieved setting

x ∈] -∞, +∞[ → ϕ γ (x) if x ∈ [0, +∞[, +∞ if x ∈] -∞, 0[. (2)
when domϕ = R + / {0} . Note that for the χ 2 -divergence, the corresponding ϕ function φ 2 (x) := 1 2 (x -1) 2 is defined and convex on whole R.

We will only consider divergences defined through differentiable functions ϕ, which we assume to satisfy (RC)

There exists a positive δ such that for all c in [1 -δ, 1 + δ], we can find numbers c 1 , c 2 , c 3 such that ϕ(cx) ≤ c 1 ϕ(x) + c 2 |x| + c 3 , for all real x. Condition (RC) holds for all power divergences including KL and KL m divergences.

Dual form of the divergence and dual estimators in parametric models

Let θ and θ T be any parameters in Θ. We intend to provide a new expression for φ(θ, θ T ). By strict convexity, for all a and b the domain of ϕ it holds

ϕ(b) ≥ ϕ(b) + ϕ ′ (a)(b -a) (3)
with equality if and only if a = b.

Denote

ϕ # (x) := xϕ ′ (x) -ϕ(x).
For any α in Θ denote

a := dP θ dP α (x) . Define b := dP θ dP θ T (x) .
Inserting these values in (3) and integrating with respect to P θ T yields

φ(θ, θ T ) ≥ ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α dP θ T .
Assume at present that this entails

φ(θ, θ T ) ≥ ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α dP θ T (4)
for suitable α's in some set F θ included in Θ. When α = θ T the inequality in (4) turns to equality, which yields

φ(θ, θ T ) = sup α∈F θ ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α dP θ T (5)
Denote

h(θ, α, x) := ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α (6)
from which φ(θ, θ T ) = sup

α∈F θ h(θ, α, x)dP θ T . (7)
Furthermore by ( 4), for all suitable α

φ(θ, θ T ) -h(θ, α, x)dP θ T = h(θ, θ T , x)dP θ T . -h(θ, α, x)dP θ T ≥ 0 and the function x → h(θ, θ T , x) -h(θ, α, x) is non negative, due to (3). It follows that φ(θ, θ T ) -h(θ, α, x)dP θ T is zero only if h(θ, α, x) = h(θ, θ T , x)-P θ T a.
e. Therefore for any x in the support of

P θ T ϕ ′ dP θ dP θ T dP θ -ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α (x) +ϕ # dP θ dP θ T (x) = 0
which cannot hold for all x when the functions ϕ # dP θ dPα (x) , ϕ # dP θ dP θ T (x) and 1 are linearly independent, unless α = θ T . We have proved that θ T is the unique optimizer in [START_REF] Csiszár | Eine informationstheoretische Ungleichung und ihre Anwendung auf den Beweis der Ergodizität von Markoffschen Ketten[END_REF].

We have skipped some sufficient conditions which ensure that (4) holds.

Assume that

ϕ ′ dP θ dP α dP θ < ∞. (8)
Assume further that φ(θ, θ T ) is finite. Since

-ϕ # dP θ dP α (x) dP θ T ≤ φ(θ, θ T ) -ϕ ′ dP θ dP θ T dP θ ≤ φ(θ, θ T ) + ϕ ′ dP θ dP θ T dP θ < +∞ we obtain ϕ # dP θ dP α dP θ T > -∞
which entails (4). When ϕ # dP θ dPα dP θ T = +∞ then clearly , under ( 8)

φ(θ, θ T ) > ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α dP θ T = -∞.
We have proved that (5) holds when α satisfies [START_REF] Liese | Convex statistical distances[END_REF]. Sufficient and simple conditions encompassing (8) can be assessed under standard requirements for nearly all divergences. We state the following Lemma (see [START_REF] Liese | Convex statistical distances[END_REF] [START_REF] Liese | Convex statistical distances[END_REF]) and Broniatowski and Kéziou (2006) [START_REF] Broniatowski | Minimization of ϕ-divergences on sets of signed measures[END_REF], Lemma 3.2).

Lemma 1 Assume that RC holds and φ(θ, α) is finite. Then (8) holds.

Summing up, we state

Theorem 2 Let θ belong to Θ and let φ(θ, θ T ) be finite. Assume that RC holds.Let F θ be the subset of all α's in Θ such that φ(θ, α) is finite . Then

φ(θ, θ T ) = sup α∈F θ ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α dP θ T .
Furthermore the sup is reached at θ T and uniqueness holds.

For the Cressie-Read family of divergences with γ = 0, 1 this representation writes

φ γ (θ, θ T ) = sup α∈F θ 1 γ -1 dP θ dP α γ-1 dP θ - 1 γ dP θ dP α γ dP θ T - 1 γ(γ -1)
.

The set F θ may depend on the choice of the parameter θ. Such is the case for the χ 2 divergence i.e. ϕ(x) = (x -1) 2 /2, when p θ (x) = θ exp(-θx)1 [0,∞) (x). In most cases the difficulty of dealing with a specific set F θ depending on θ can be encompassed when exists a neighborhood U of θ T for which (A) φ(θ, θ ′ ) is finite whatever θ and θ ′ in U which for example holds in the above case for any θ T . This simplication deserves to be stated in the next result Theorem 3 When φ(θ, θ T ) is finite and RC holds, then under condition (A)

φ(θ, θ T ) = sup α∈U ϕ ′ dP θ dP α dP θ -ϕ # dP θ dP α dP θ T .
Furthermore the sup is reached at θ T and uniqueness holds.

Remark 4 Identifying F θ might be cumbersome. This difficulty also appears in the classical MLE case, a special case of the above statement with divergence function ϕ 0 ,for which it is assumed that log p θ (x)p θ T (x)dλ(x) is finite for θ in a neighborhood of θ T .

Under the above notation and hypotheses define

T θ (P θ T ) := arg sup α∈F θ h(θ, α, x)dP θ T . (9)
It then holds

T θ (P θ T ) = θ T
for all θ T in Θ. Also let

S (P θ T ) := arg inf θ∈Θ sup α∈F θ h(θ, α, x)dP θ T . (10)
which also satisfies

S (P θ T ) = θ T
for all θ T in Θ. We thus state Theorem 5 When φ(θ, θ T ) is finite for all θ in Θ and RC holds, both functionals T θ and S are Fisher consistent for all θ T in Θ.

Plug in estimators

From [START_REF] Keziou | Dual representation of ϕ-divergences and applications[END_REF] simple estimators for θ T can be defined, plugging any convergent empirical measure in place of P θ T and taking the infimum in θ in the resulting estimator of φ(θ, θ T ).

In the context of simple i.i.d. sampling, introducing the empirical measure

P n := 1 n n i=1 δ X i
where the X i 's are i.i.d. r.v's common unknown distribution P θ T in P, the natural estimator of φ(θ, θ T ) is

φ n (θ, θ T ) := sup α∈F θ h(θ, α, x) dP n (x) (11)
= sup

α∈F θ ϕ ′ dP θ dP α dP θ - 1 n n i=1 ϕ # dP θ dP α (X i ) . Since inf θ∈Θ φ(θ, θ T ) = φ(θ T , θ T ) = 0 the resulting estimator of φ(θ T , θ T ) is φ n (θ T , θ T ) := inf θ∈Θ φ n (θ, θ T ) = inf θ∈Θ sup α∈F θ h(θ, α, x) dP n (x) . ( 12
)
Also the estimator of θ T is obtained as

θ := arg inf θ∈Θ sup α∈F θ h(θ, α, x) dP n (x) . (13)
When A holds then F θ may be substituted by U in the above definitions. The resulting minimum dual divergence estimators (12) and (13) do not require any smoothing or grouping, in contrast with the classical approach which involves quantization. The paper [START_REF] Broniatowski | Parametric estimation and tests through divergences and the duality technique[END_REF] provides a complete study of those estimates and subsequent inference tools for the usual i.i.d. sample scheme. For all divergences considered here, these estimators are asymptotically efficient in the sense that they achieve the Cramer-Rao bound asymptotically. The case when ϕ = ϕ 0 leads to θ M L defined as the celebrated Maximum Likelihood Estimator (MLE), in the context of the simple sampling.

Minimum divergence estimators in exponential families

In this section we prove the following result Theorem 6 For all divergence φ defined through a differentiable function ϕ satisfying Condition (RC), the minimum dual divergence estimator defined by ( 13) coincides with the MLE on any full exponential families such that φ (θ, α) is finite for all θ and α in Θ.

Let P be an exponential family on R s with canonical parameter in

R d P := P θ such that p θ (x) = dP θ dλ (x) = exp [T (x) ′ θ -C(θ)] ; θ ∈ Θ where x is in R s and Θ is an open subset of R d ,
and λ is a dominating measure for P. We assume P to be full, namely that the Hessian matrix (∂ 2 /∂θ 2 ) C(θ) is definite positive for all θ in Θ.

Let X 1 , ..., X n be n i.i.d. random variables with common distribution

P θ T with θ T in Θ. Introduce M n (θ, α) := ϕ ′ dP θ dP α dP θ - 1 n n i=1 ϕ # dP θ dP α (X i ) We will prove that inf θ sup α M n (θ, α) = 0 (14)
whatever the function ϕ satisfying the claim. In (14) θ and α run in Θ. This result extends the maximum likelihood case for which inf

θ sup α M n (θ, α) = sup θ inf α 1 n n i=1 log p θ (X i ) -1 n n i=1 log p α (X i ) = 0. Direct substitution shows that for any θ, sup α M n (θ, α) ≥ M n (θ, θ) = 0 from which inf θ sup α M n (θ, α) ≥ 0 (15)
We prove that

α = θ M L is the unique maximizer of M n (θ M L , α) (16

Mϕ

  n (θ, α) ≤ sup α M n (θ M L , α) = M n (θ M L , θ M L ) = 0 (17)which together with (15) completes the proof.DefineM n,1 (θ, α) := ϕ ′ (exp A(θ, α, x)) exp B (θ, x) dλ(x) M n,2 (θ, α) := 1 n n i=1 exp (A (θ, α, X i )) ϕ ′ (exp A(θ, α, X i )) M n,3 (θ, α) (exp A(θ, α, X i ))withA(θ, α, x) := T (x) ′ (θ -α) + C(α) -C(θ) B(θ, x) := T (x) ′ θ -C(θ).It holdsM n (θ, α) = M n,1 (θ, α) -M n,2 (θ, α) + M n,3 (θ, α)with∂ ∂α M n,1 (θ, α) α=θ = -ϕ (2) (1) [∇C (θ) -∇C (α) α=θ ] = 0 for all θ, ∂ ∂α M n,2 (θ M L , α) α=θ M L = ϕ (2) X i ) + ∇C (α) α=θ M L X i ) + ∇C (α) α=θ M L = 0where the two last displays hold iff α = θ M L . Now∂ 2 ∂α 2 M n,1 (θ M L , α) α=θ M L = ϕ (3) (1) + 2ϕ (2) (1) ∂ 2 /∂θ 2 C(θ M L ) ∂ 2 ∂α 2 M n,2 (θ M L , α) α=θ M L = ϕ (3) (1) + 4ϕ (2) (1) ∂ 2 /∂θ 2 C(θ M L ) ∂ 2 ∂α 2 M n,3 (θ M L , α) α=θ M L = ϕ (2) (1) ∂ 2 /∂θ 2 C(θ M L ), whence ∂ ∂α M n (θ M L , α) α=θ M L = 0 ∂ 2 ∂α 2 M n (θ M L , α) α=θ M L = -ϕ (2) (1) ∂ 2 /∂θ 2 C(θ M L )which proves (16), and closes the proof.


Julien Marzat

Hélène Piet-Lahanier

Frédéric Damongeot

Eric Walter

Tuning and comparing fault diagnosis methods for aeronautical systems via Kriging-based optimization

Many approaches address fault detection and isolation (FDI) based on analytical redundancy. To rank them, it is necessary to define performance indices and realistic sets of test cases on which these performance indices will be evaluated. For the ranking to be fair, each of the methods under consideration should have its internal parameters (often called hyperparameters) tuned optimally. However, no mathematical model linking hyperparameters and performance is available a priori. In this paper, we propose to use a combination of tools developed in the context of computer experiments to build such a model from a limited number of numerical evaluations of the performance indices at carefully chosen values of the hyperparameters. The optimal tuning of fault diagnosis methods may prove to be strongly sensitive to specifics of the test cases. This is why the methodology is extended so as to provide a tuning that is robust to variability in the conditions of use. The performance criteria are then replaced by their worst-case values when the sources of variability are assumed to belong to some predefined sets. This methodology is applied to tune fault diagnosis approaches on an aeronautical case study.

Introduction

Many approaches address fault detection and isolation (FDI) based on analytical redundancy. A review of the main methods can be, e.g., found in [START_REF] Marzat | Autonomous fault diagnosis: State of the art and aeronautical benchmark[END_REF][START_REF] Henry | Fault detection and diagnosis for aeronautic and aerospace missions[END_REF][START_REF] Ding | Model-based Fault Diagnosis Techniques: Design Schemes, Algorithms, and Tools[END_REF]. Ranking these methods makes it necessary to define performance indices and realistic sets of test cases on which these performance indices will be evaluated. For the ranking to be fair, each of the methods under consideration should have its internal parameters (often called hyperparameters) tuned optimally. This is a key issue, as the potential of each method derives directly from this tuning.

However, except in trivial cases, no mathematical model linking hyperparameters and performance is available a priori. In this paper, we propose to use a combination of tools developed in the context of computer experiments [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] to build such a model from a limited number of numerical evaluations of the performance indices at carefully chosen values of the hyperparameters. The methodology to be presented is particularly well-suited to computationally-intensive performance indices, which are the rule if the test cases are realistic enough. The approach is based on Kriging [START_REF] Matheron | Principles of geostatistics[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], which is used to build approximate models linking hyperparameters to performance indices, based on a small number of simulation results. These models are then used to optimize the hyperparameters using a procedure known as Efficient Global Optimization (EGO) [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF][START_REF] Marzat | Automatic tuning via Kriging-based optimization of methods for fault detection and isolation[END_REF].

Optimal tuning of fault diagnosis methods should not be too sensitive to specifics of the test cases. This is why the methodology is extended so as to provide a tuning that is robust to variability in the conditions of use. The performance criteria are then replaced by their worst-case values when the sources of variability are assumed to belong to some predefined set.

The procedure, in its basic and robust versions, is illustrated by the tuning and comparison of two classical fault diagnosis architectures applied on an aeronautical test case.

Problem formulation

This paper focuses on fault detection in aeronautics, for instance applied to an air-to-air missile (see Section 3.3). Several methods can be used for detecting faults occurring on this system such as an observer-based residual generator or other types of filters. The quality of the results obtained via each method depends strongly on the adequacy of its tuning to the problem to be treated. This tuning is parameterized by a vector of hyperparameters, defined as x c ∈ X c , where X c is assumed to be a known compact set. The performance measure of fault detection is assumed to be quantified by the real scalar y(x c ) to be minimized. The numerical value of y(x c ) for any feasible numerical value of x c is obtained by simulating fault scenarios for the aeronautical vehicle and evaluating the quality of detection using the method under consideration for this value of the hyperparameter vector. The complex simulation of such a test case is illustrated in Figure 1, including the dynamics of the aircraft, its control loop and a fault detection scheme composed of a residual generator and an evaluation test. Assuming that the output of the detection procedure is a Boolean signature (see Figure 2) on whether a specific fault has appeared, the performance measure y(•) can be defined as some trade-off between false-alarm and nondetection rates. The test is supposed to run from t on to t hor , and the fault starts at t from ∈]t on , t hor [. Then, the false-alarm rate is

r fd = i t i fd t from -t on , (1)
where t i fd is the i-th period of time where the decision is true in [t on ; t from ] and the non-detection rate

r nd = i t i nd t hor -t from , (2)
where t i nd is the i-th period of time where the decision is false in [t from ; t hor ].

The cost function to be used in this paper to measure performance of fault diagnosis methods is

y(x c ) = r fd + r nd . (3)
This is a simple trade-off between the conflicting objectives of minimizing false alarms and non detections. Tuning can then be formalized as the search for x c = arg min

x c ∈X c y(x c ). ( 4
)

Tuning via Kriging

As y(•) can only be evaluated via possibly very costly numerical simulation at sampled points, a black-box global optimization method is used, which has become very popular in the context of computer experiments. It has been applied in [START_REF] Marzat | Automatic tuning via Kriging-based optimization of methods for fault detection and isolation[END_REF] to tune several fault detection methods on the simple example of a change in the mean of a signal subject to Gaussian or uniformly distributed noise. The overall procedure is recursive and presupposes that a set of values of the performance index has already been computed, y c,n c = [y(x c,1 ), ..., y(x c,n c )] T , corresponding to an initial sampling of n c points in

X c , X c,n c = [x c,1 , ..., x c,n c ].
First, a Kriging model (Section 3.1) is fitted to approximate the mapping from the hyperparameter vector x c to the scalar performance measure y(•). Then, an iterative optimization procedure (Section 3.2), taking advantage of the Kriging model, is employed to look for a feasible global minimizer x c of the performance index.

Kriging

The unknown performance function y(•) is approximated by

Y(x c ) = f T (x c ) b + Z(x c ), (5)
where f T (x c ) b is a classical regression model with the entries of f (x c ) polynomial in the entries of x c and the regression coefficients b to be estimated. Z(•) is a zero-mean Gaussian process with known (or parametrized) covariance function k (•, •). Based on this modelling, Kriging is the search for the best linear unbiased predictor (BLUP) of Y(•) [START_REF] Matheron | Principles of geostatistics[END_REF][START_REF] Kleijnen | Kriging metamodeling in simulation: A review[END_REF]. The covariance function k (•, •) is usually considered as unknown, and expressed as

k Z x c,i , Z x c, j = σ 2 Z R x c,i , x c, j , (6)
where σ 2 Z is the process variance and R (•, •) some parametric correlation function. The parameters of R (•, •) and σ 2 Z must be fixed or estimated from available data. R x c,i , x c, j is often assumed to depend only on the displacement vector h = x c,ix c, j . In this paper, we adopt the classical power exponential correlation function [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF],

R (h) = exp         - d k=1 h k θ k 2         , (7)
where h k is the k-th component of h. For our application, empirical Kriging was used, which means that the covariance parameters are estimated in the maximum-likelihood sense. To simplify presentation, they will be considered as known in what follows and will not appear in the formulas. Define R as the n × n matrix such that

R(i, j) = R x c,i , x c, j , (8)
r(x c ) as the n vector r

(x c ) = R x c , x c,1 , ..., R x c , x c,n T (9)
and F as the n × dim b matrix

F = f x c,1 , ..., f x c,n T . ( 10
)
The maximum-likelihood estimate b of the vector of the regression coefficients b is

b = F T R -1 F -1 F T R -1 y c,n c . ( 11
)
The prediction of the mean of the Gaussian process at x c ∈ X c is then

Y(x c ) = f T (x c ) b + r (x c ) T R -1 y c,n c -F b . ( 12
)
Another crucial property of Kriging that will be exploited in the search for a global optimization of the performance index is the ability to compute the variance of the prediction error [START_REF] Schonlau | Computer Experiments and Global Optimization[END_REF] as

σ 2 (x c ) = σ 2 Z 1 -r (x c ) T R -1 r (x c ) . (13)
This variance can be used to assess the confidence in the prediction (12) at any point x c ∈ X c . In particular, it is equal to zero for all the sampled points in X c,n c if simulation errors are deemed negligible as here.

Efficient Global Optimization (EGO)

The idea of EGO [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] is to use the Kriging predictor ( 12) and the additional confidence information provided by [START_REF] Frank | Survey of robust residual generation and evaluation methods in observer-based fault detection systems[END_REF] to find an (n c + 1)-st vector in X c at which the complex simulation should be run. This point is chosen by maximizing a criterion that measures the interest of an additional evaluation at x c , given the past results y c,n c obtained at X c,n c . A possible choice for this criterion is the Expected Improvement (EI) [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF],

EI(x c , y n c min ) = σ (x c ) uΦ (u) + φ (u) , ( 14
)
where Φ is the cumulative distribution function and φ the probability density function of the normalized Gaussian distribution N (0, 1), with also

u = y n c min -Y (x c ) σ (x c ) , (15)
y n c min = min i=1...n c y x c,i . (16)
The maximization of EI achieves a trade-off between local search (numerator of ( 15)) and the exploration of unknown areas (where σ is high), and is thus well suited for global optimization. The EGO algorithm proceeds as follows.

1. Choose an initial sampling X c,n c and compute y c,n c .

2. Fit a Kriging predictor on those data points.

3. Compute ( 16) and find a new point x n c+1 where to evaluate y (•) by maximizing [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF] 4. Append x n c+1 to X c,n c , append y(x n c+1 ) to y c,n c and increment n c by one.

5. Go to step 2 until {max EI(•)} hits a lower threshold ε or until n c exceeds a predefined sampling budget n max .

The estimate of the global minimum is finally given by a last call to [START_REF] Chow | Analytical redundancy and the design of robust failure detection systems[END_REF], and the corresponding argument in X c,n c is an estimate of a global minimizer. To obtain the results presented in this paper, we used the toolbox SuperEGO [START_REF] Sasena | Flexibility and Efficiency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF]. Its parameters have been set to ε = 10 -4 and n max = 100.

Illustrative example

The test case is the reduced longitudinal model of a missile flying at a constant altitude of 6000 m. The state vector, consisting of the angle of attack, angular rate and Mach number, is x = [α, q, M] T . The control input is the rudder angle, u = δ, and the available measurement is the normal acceleration, γ = a z . The linearized model around the operating point x 0 = [ ᾱ, q, M] T = [20 deg, 18.4 deg/s, 3] T is given by the following state-space model, after discretization with a time step of 0.02s,

       x k+1 = Ax k + Bu k γ k = Cx k + Du k + w k + f k (17)
where

A =          
0.9163 0.0194 0.0026 -5.8014 0.9412 0.5991 -0.0485 -0.005 0.996

          , B =           -0.0279 -2.5585 -0.0019           , C = -2.54 0 -0.26 , D = -0.204 . ( 18
)
This model is simulated on a time horizon of 50 seconds. Measurement γ = a z suffers at time = 25 s from a progressive fault f steadily increasing with a slope equal to 0.1. This measurement is also corrupted by a Gaussian noise w, with zero mean and standard deviation 10 -3 .

Fault diagnosis methods

Two methods are considered to detect this sensor fault, in line with the architecture described in Figure 1. The first one uses a Luenberger observer and the other is based on a Kalman filter [START_REF] Frank | Survey of robust residual generation and evaluation methods in observer-based fault detection systems[END_REF]. A residual is produced by each filter as the difference between the available measurement and the predicted output value. In both schemes, the residuals are analyzed using a CUSUM test [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF].

The performance of a Luenberger observer depends on the settings of its three poles [p 1 , p 2 , p 3 ], while the tuning of a Kalman filter relies on the choice of an initial covariance matrix W for the state vector perturbation (dimension 3) and of a covariance matrix V for measurement noise (dimension 1), both matrices are assumed to be diagonal. This corresponds to 4 hyperparameters [w 1 , w 2 , w 3 , v 1 ], such that

W =           w 1 0 0 0 w 2 0 0 0 w 3           , V = v 1 (19)
The initial value of the covariance matrix P of the state estimation error is arbitrarily set to identity (dimension 3). Initial estimated state vector is identical for both filters and taken as [17, 21, 2.5] T . The CUSUM test has two hyperparameters, namely the minimum size of change to be detected µ and a threshold λ. The tuning hyperparameters of each scheme are those of the residual generator and of the statistical test. Search spaces for these values for each of the competing methods are as indicated in Table 1.

Table 1: Space search of hyperparameters for the two diagnostic methods

Results

One hundred initialisations with the tuning methodology have been effected. The number of simulations required for obtaining a suitable tuning remains low (Table 2), and so is the variation on the results. The use of a Kalman filter proves to be more efficient in this setup than that of a Luenberger observer. The resulting decisions are illustrated on Figure 3 .

Extension to robustness to environmental conditions

The previous approach does not take into account the variability of environmental conditions (model and measurement uncertainty, variation in atmospheric conditions or magnitude of faults, etc). However, the simulation response can vary greatly with these conditions at the risk of making the optimal tuning inadequate. Environmental conditions are assumed to be described by a vector of environmental variables x e belonging to a known compact set X e . The tuning of a fault detection method should remain valid for a reasonable set of possible values of the environmental variables. Most of the studies on computer experiments in this context use a probabilistic modeling of the environmental variables [START_REF] Lehman | Designing computer experiments to determine robust control variables[END_REF]. We assume instead that bounds of the values taken by these variables are available and look for an optimal tuning in the worst-case sense. Worst-case optimality is a concern that has been raised in many studies on fault detection, since [START_REF] Frank | Survey of robust residual generation and evaluation methods in observer-based fault detection systems[END_REF][START_REF] Chow | Analytical redundancy and the design of robust failure detection systems[END_REF].

Since no algorithm seemed to have been reported to deal with environmental variables for the robust optimization of black-box functions evaluated by costly computer experiments, we have proposed such a procedure in [START_REF] Marzat | Min-max hyperparameter tuning with application to fault detection[END_REF]. The problem to be solved can be formalized as that of finding x c and x e such that x c , x e = arg min

x c ∈X c max x e ∈X e y(x c , x e ), (20)
which is a difficult problem since both spaces X c and X e are now continuous. This corresponds to the search for the best hyperparameters for the worst values of the environmental variables. The optimization procedure proposed to solve (20) combines the relaxation algorithm described in [START_REF] Shimizu | Necessary conditions for min-max problems and algorithms by a relaxation procedure[END_REF] with EGO (Section 3.2). A detailed description of the resulting robust tuning algorithm can be found in [START_REF] Marzat | Min-max hyperparameter tuning with application to fault detection[END_REF]. The main principles are to transform (20) into

         min x c ∈X c ,τ τ subject to y(x c , x e ) ≤ τ, ∀x e ∈ X e , (21)
and then to solve this minimization problem (that has an infinite number of constraints) by an iterative relaxation procedure that computes a finite number of constraints from values in X e . The resulting procedure involves two intertwined global-optimization steps, which are addressed with EGO.

Illustrative example

The test case and the methods to be compared remain the same as in Section 3. ]. The hyperparameter spaces for the two schemes are as before (see Table 1).

One hundred runs of the entire procedure have again been performed to assess convergence, repeatability and dispersion of its results. Performance level, number of evaluations and values of the worst-case environmental variables and of the best hyperparameter tuning for both schemes are reported in Table 3. It is important to note that the number of evaluations is quite low with an average sampling of approximately 33 points per dimension, leading to a quick robust tuning. Moreover, the repetition of the procedure suggests that on this example an acceptable value for the tuning is always obtained, and that the worst-case is correctly identified. Indeed, for both strategies, the worst environmental conditions are located near the smallest value of the fault and highest value of the noise, which is what common sense suggests. Figure 4 shows the decision functions obtained for both diagnosis methods tuned at the best hyperparameters under the conditions of the worst-case environmental conditions. It clearly appears on this example that the Kalman filter scheme is much better at detecting small incipient faults than the observer scheme, which is not conceived to take into account such noise level. Figure 5 shows the value of the objective function y over X e for the estimated worst-case optimal tuning of the hyperparameters of the two methods. This shows that the optimal worst-case tuning leads to good performance for other environmental conditions.

Conclusions and perspectives

This article has presented a framework for automatically tuning fault diagnosis methods. The setup only requires a simulation of a collection of test cases, where the internal parameters of candidate methods may be adjusted and a measure of performance computed. The tuning methodology relies on fitting the mapping from hyperparameters to performance level via Kriging, which is then used to find the optimal hyperparameters. Since any simulation should take into account uncertainty, a robust version of this tuning procedure has been developed in the worst-case sense. It allows one to find at the same time the optimal hyperparameters and the worst environmental conditions. These strategies have been illustrated on an academic version of an aeronautical case study, where two competing fault diagnosis schemes comprising a residual generator and a statistical test were applied. The number of evaluations of the simulation to achieve an acceptable tuning value is always very low, which makes it practicable when the associated computational cost is expensive. The proposed procedures may be readily extended to other more realistic aerospace applications.

Figure 1 :

 1 Figure 1: Simulation of the test case

Figure 2 :

 2 Figure 2: Boolean decision function

Observer and CUSUM p 1 ∈

 1 [0; 1] p 2 ∈ [0; 1] p 3 ∈ [0; 1] µ ∈ [0.01; 5] λ ∈ [0.1; 20] Kalman and CUSUM w 1 ∈ [0; 2] w 2 ∈ [0; 2] w 3 ∈ [0; 2] v 1 ∈ [0; 10] µ ∈ [0.01; 5] λ ∈ [0.1; 20] (a) Observer-based scheme (b) Kalman-based scheme

Figure 3 :

 3 Figure 3: Boolean decision functions for both schemes with best hyperparameter tuning

3 .

 3 The sources of variability are the magnitude of the incipient sensor fault f and the noise level on w. The slope of the fault s and the standard deviation ζ of the Gaussian white noise are then the environmental variables to which the tuning should be robust, x e = [ζ, s] T . The search space X e is such that ζ ∈ [10 -7 , 10 -3 ] and s ∈ [10 -3 , 10 -1

Fault slope s 1 • 10 -3 1 • 10 - 3

 1101103 Hyperparameter vector x c Hyperparameter values p 1 = 0.73 p 2 = 0.726 w 1 = 1.57 w 2 = 1.11 for the median performance p 3 = 0.72 w 3 = 1.06 v 1 = 2.04 µ = 0.065 λ = 4.553 µ = 0.12 λ = 3.33 (a) Observer-based scheme (b) Kalman-based scheme

Figure 4 :

 4 Figure 4: Boolean decision function for both schemes with robust tuning and the worst-case environmental variables

  (a) Observer-based scheme (b) Kalman-based scheme

Figure 5 :

 5 Figure 5: Value of the objective function over X e for the minimax-optimal hyperparameters

Table 2 :

 2 Tuning of complete decision approaches

		Observer and CUSUM	Kalman and CUSUM
	Ranking	2	1
	Median of performance	0.0455	0.0184
	standard deviation of performance	0.009	0.0024
	Average number of simulations	102	136
	standard deviation of number of simulations	20	21
	False alarm rate	0	0
	Non detection rate	0.0455	0.0184
	Hyperparameter values	p

1 = 0.7675 p 2 = 0.8498 w 1 = 0.5428 w 2 = 0.2593 for median of performance p 3 = 0.8951 w 3 = 1.1015 v 1 = 2.2291 µ = 0.4537 λ = 2.91 µ = 0.134 λ = 4

Table 3 :

 3 Results for 100 replications of the minimax tuning procedure

		Observer and CUSUM	Kalman and CUSUM
	Ranking	2	1
	Median minimax performance	0.114	0.0312
	standard deviation of performance	4.7 • 10 -2	1.97 • 10 -2
	Average number of simulations	168	199
	standard deviation nb simulations	26	8
		Environmental vector x e
	Noise level ζ	9.8 • 10 -4	9.81 • 10 -4