Hassan El Abiad
email: hassan.elabiad@supelec.fr

Laurent Le Brusquet
email: laurent.lebrusquet@supelec.fr

Marie-Ève Davoust
email: marie-eve.davoust@supelec.fr

MODEL-ROBUST DESIGN OF EXPERIMENTS FOR SEQUENTIAL IDENTIFICATION OF ODE PARAMETERS

This paper presents the idea of sequential model-robust Design of Experiments (DOE) for the identification of dynamic systems modeled with an Ordinary Differential Equation (ODE). The studied DOE problem consists in selecting sequentially the instants where the measures will be done in order to best estimate the system's parameter. The robustness is achieved by considering a statistical representation of the model error defined as the difference between the true ODE and the ODE used in the model. The idea of modeling the model error with a statistical representation has been widely explored in the DOE literature for the identification of static systems. However, there have been little previous works that apply this idea for the identification of dynamic systems. This paper initiates an exploration of this idea in the context of first-order ODE. The model error is modeled by using a kernel-based representation (Gaussian process). A new criterion for the instant selection is constructed and tested on an illustrative example. The design reached with the proposed sequential robust criterion is compared with the design reached with the non-robust version of criterion and with the classical uniform design.

INTRODUCTION

For parametric identification problems, the purpose of DOE is to adjust the selection of the experimental conditions to improve the quality of the parameter estimation. This has motivated many researches on DOE over almost a century [START_REF] Wald | On the efficient design of statistical investigations[END_REF][START_REF] Wolfowitz | The equivalence of two extremum problems[END_REF][START_REF] Wynn | Results in the theory and construction of d-optimum experimental designs[END_REF][START_REF] Federov | Theory of optimal experiments[END_REF][START_REF] Whittle | Some general points in the theory of optimal experimental design[END_REF]. The DOE technique was later adapted to the problem of dynamic systems identification [START_REF] Levadi | Design of input signals for parameter estimation[END_REF][START_REF] Arimoto | Optimal input test signals for system identification -an information theoretic approach[END_REF][START_REF] Mehra | Optimal inputs for system identification[END_REF]. Some of the recent works are summarized in [START_REF] Gevers | Identification for control: From the early acheivements to the revival of experiment design[END_REF].

There have been several works on robust DOE, especially for the identification of static models for engineering design problems. Design robustness with respect to the model error (misspecification) was first discussed by Box and Draper [START_REF] Box | A basis for the selection of a response surface design[END_REF] who studied the effect of taking a onedegree polynomial regression model when the target is twodegree. After that, many authors have further discussed and developed this idea with different assumptions about the misspecification [START_REF] Montepiedra | Minimum bias designs with constraints[END_REF][START_REF] Yue | Robust designs for fitting linear models with misspecification[END_REF][START_REF] Fang | Robust regression designs for approximate polynomial models[END_REF].

Several recent papers have considered these robust DOE techniques for the identification of ODE parameters. Aspery and Macchietto [START_REF] Asprey | Designing robust optimal dynamic experiments[END_REF] considered the problem of design robustness with respect to the parameter starting values. They considered that the design is affected by a poor starting values of the parameters and try to derive a criterion that is robust against these starting values. In [START_REF] Goodwin | Robust optimal experiment design for system identification[END_REF], Goodwin et al. propose a min-max design criterion for the parameter identification of a linear first-order ODE:

ẋ(t) = a * x(t) + u(t) + r(t) (1)
where x(t) is the state variable, u(t) is the control input, and r(t) is an unknown departure from the simple model ẋ(t) = a * x(t)+u(t) (misspecification). All these functions are scalar and defined on a finite interval [0, t max ]. Scalar a * is the parameter to be identified.

We have studied the same identification problem proposed by Goodwin et al. (equation ( 1)) but we have considered a statistical approach of the model error r(t) instead of the bound approach derived from the used min-max criterion.

The presence of a non-negligible misspecification has a great influence on the identification of the parameter a * and also on the selection of the instants t i where the state variable x(t) will be measured. Misspecified ODE may encounter frequently, especially in the biological field where the main reaction to be modeled is polluted by additional terms whose definition by a specific function would not be feasible.

We define a * as the value that minimizes the L 2 -norm of the misspecification r(t):

a * = arg min a tmax 0 ( ẋ(t) -ax(t) -u(t)) 2 dt (2)
For a practical reason justified in section 3.1, we choose to estimate τ * = -1 a * , the time constant of the system. The aim of designing experiments is to reach an accurate estimation of τ * with a small number of experiments. This is particularly important when obtaining experiments is difficult, for example for material or cost reasons.

The estimation of τ * is performed from a database of collected data. Suppose that the initial point (t 0 = 0, y 0 = x(t 0 )) is known and that a set of n collected data {(t i , y i ) ∈ ]0, t max ] × R, i = 1, ..., n} has already been collected. The t i 's form the initial design denoted by ξ n = [t 1 , ..., t n ] ⊤ , the y i 's are noisy observations of x i = x(t i ). Therefore, the observation equation is defined as:

y 0 = x 0 y i = x i + e i , i = 1 . . . n
where the observations error e i are normal and i.i.d. (e i ∼ N (0, σ 2 e ), σ 2 e is supposed to be known). Suppose that we desire to refine the estimation of τ * by adding a new design point t n+1 and its corresponding observation value y n+1 to the collected data. Then, the problem of sequential DOE is to choose the next design point t n+1 that will best refine the estimation.

An important issue in designing experiments for dynamical system identification is the fact that design criteria generally depend on the unknown parameter which has to be identified. To solve this problem, a sequential DOE strategy is generally implemented. This strategy consists in alternating between parameter estimation stages and designing stages, as shown in figure 1 [START_REF] Wolfowitz | The equivalence of two extremum problems[END_REF][START_REF] Wu | Asymptotic inference from sequential design in a nonlinear situation[END_REF][START_REF] Kulcsar | The interface between sequential experimental design and dual control[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF]. In the design stage, the estimated parameter τn is used in order to find the next design point t n+1 :

t n+1 = arg min t>tn J design (t, τn )
The paper is organized as follows. The next section illustrates the proposed approach and in particular the statistical representation used for r(t) and the derived design criterion J design (t, τn ). Mathematical developments for calculating the criteria J estim. (τ, ξ n , y n ) and J design (t, τn ) are not straightforward. They are detailed in section 3. In the fourth section, an example is presented in order to illustrate the improvements of the proposed sequential design over the classical uniform design and the non-robust version of the proposed sequential design. Section five concludes the paper by presenting the main contributions.

PROPOSED APPROACH SCHEME

This section presents the proposed sequential robust criterion J design (t, τn ). As mentioned in the introduction, the main contribution of this work is the presence of the model error r(t) which improves the robustness of our design criterion.

The model error r(t) is generally a relative smooth function. This is the reason why we model it by a Gaussian process as in [START_REF] Roger | A criterion for model-robust design of experiments[END_REF][START_REF] Abiad | Model-Robust Sequential Design of Experiments for Identification Problems[END_REF] (see [START_REF] Williams | Prediction with Gaussian processes: From linear regression to linear prediction and beyond[END_REF] for a Gaussian process tutorial). A Gaussian process is a random field defined by its mean and covariance function:

E r [r(t)] = 0, ∀ t ∈ [0, t max ] E r [r(t)r(t ′ )] = c(t -t ′ ), ∀ (t, t ′ ) ∈ [0, t max ] 2
where E[.] denotes the expected value and s 2 = c(0) denotes the variance of the model error. The relevance of modeling the misspecification by a Gaussian process rises because for some classes of covariance functions, Gaussian processes span a rather large space (infinite-dimensional). Therefore, this type of representation matches the robustness requirement: the design point t n+1 that we look for must lead to a good estimation performance whatever is the misspecification. The misspecification is unknown but our idea of robustness states that the chosen design point must guarantee a good level of performance (on average) over a wide range of potential misspecifications. The theoretical expression of the proposed criterion is directly derived from the use of the statistical representation of r(t) that was suggested:

J design (t, τn ) = E {en,en+1,r} (τ n+1 (t) -τ * ) 2 ( 3
)
where τn+1 is the updated estimation of τ * and e n = [e 1 , ..., e n ] ⊤ .

The use of the expected value over r (i.e. r(t), t ∈ [0, t max ])

assures the robustness with respect to the model misspecification.

MATHEMATICAL DEVELOPMENTS OF THE CRITERIA

This section presents the mathematical developments of the estimation and design criteria introduced in the previous section.

Estimation stage criterion

The first stage of the DOE procedure consists in estimating the unknown parameter τ * given n observations. This section presents a detailed formulation of the parameter estimation procedure. Two estimators of τ * are presented: the proposed robust estimator (considering model-error r(t) and measurement error e(t)) and its non-robust version (considering only measurement error e(t)). The non-robust version will be used for comparison in section 4.

First, let's use a numerical integration scheme in order to discretize the state equation (1) at the instants t i :

x(t i ) = x 0 + a * ti 0 x(θ)dθ + ti 0 u(θ)dθ + ti 0 r(θ)dθ (4)
where i = 1, ..., n. A numerical integration scheme is preferred over a derivative one for stability reasons.

The second term in ( 4) is approximated by the trapezoidal integration method. Thus, x n = [x(t 1 ), ..., x(t n )] ⊤ can be written in vector form:

x n = x 0 + a * W n x n + u n + r n ( 5
)
where u n = t1 0 u(θ)dθ, ..., tn 0 u(θ)dθ ⊤ is the control input vector defined by the input signal u(t) which is known. The n × n matrix W n contains the coefficients of the trapezoidal integration method.

The estimation of τ * is performed from the collected data y n = [y(t 1 ), ..., y(t n )]

⊤ which are noisy observations of (5):

y n = x 0 + a * W n (y n -e n ) + u n + e n + r n
The estimator τn must be coherent with the definition of τ * . In (2), τ * is obtained by minimizing the L 2 -norm, which is a natural Maximum Likelihood criterion when errors are i.i.d. The same assumptions are kept to construct the estimator τn . Thus, the obtained estimator will be consistent (τ n → τ * ), which is the main importance for an estimator. We will see in section 3.2 that a more realistic covariance function will be used to compute a precise value of E {en,r}

(τ n -τ * )
2 (main importance for the evaluation of the quality of an estimator). The τn estimator is thus obtained by considering the Maximum Likelihood of the vector e n + r n with r(t) supposed to be i.i.d. (r(t) ∼ N (0, s 2 )). When supposing r(t) i.i.d., the vector r n is a Gaussian vector. If k ≤ j, the (k, j) coefficient of its n × n covariance matrix D n is equal to s 2 t k . Thus, the Likelihood criterion is constructed with a covariance matrix equal to (D n + σ 2 e I n ) with I n a n × n identity matrix. Therefore, The estimator expression is given by:

τn = y ⊤ n W ⊤ n (D n + σ 2 e I n ) -1 W n y n y ⊤ n W ⊤ n (D n + σ 2 e I n ) -1 [y n -v n ] ⊤ (6)
where

v n = x 0 + u n .
The non-robust estimator is simply obtained by considering D n = 0 in the previous equation.

Design stage criterion

The expectation computation in (3) may be reached because, for a given misspecification r, a given design, and a given noise realization [e n , e n+1 ], it is possible to write analytically the expression of τn+1 . Indeed, the solution of (1) is:

x(t) = x 0 exp(a * t)+exp(a * t) tmax 0 (r(θ)+u(θ)) exp(-a * θ)dθ (7)
The integral terms in the previous equation are obtained by a fine discretization of the experimental domain [0, t max ]. Let {t ′ i } i=1...m be this time discretization. In the following, subscript ( ′ ) denotes vector calculated at these instants (e.g. r ′ m = [r(t ′ 1 )...r(t ′ m )] ⊤ ). Thus, equation ( 7) can be written in vector form:

x ′ m = P ′ m (r ′ m + u ′ ) + f ′ m
where P ′ m is the m × m matrix of integration coefficients (classical trapezoidal integration including exponential function) and

f ′ m = [x 0 exp(a * t ′ 1 ) . . . x 0 exp(a * t ′ m )] ⊤ is the function f (t) = x 0 exp(a * t) evaluated at the m discretiza- tion instants.
Following the same procedure, the n collected data can be written:

y n = P n (r ′ m + u ′ m ) + f n + e n ( 8
)
where P n is a n × m matrix containing the integration coefficients for the n collected data, and f n is the function f (t) evaluated at observation instants (t 1 , . . . , t n ).

The sequential DOE consists in choosing the next design point t n+1 which will refine the parameter estimation.

To estimate the unknown parameter τ * using n + 1 observations, we have to adopt the estimation expression in (6): 8), the expression of the parameter estimation ( 9) is rewritten in the following form:

τn+1 = y ⊤ n+1 W ⊤ n+1 (D n+1 + σ 2 e I n+1 ) -1 W n+1 y n+1 y ⊤ n+1 W ⊤ n+1 (D n+1 + σ 2 e I n+1 ) -1 [y n+1 -v n+1 ] ⊤ (9) Using (
τn+1 = r ′ ⊤ m M num r ′ m + r ′ ⊤ m v num + c num r ′ ⊤ m M den r ′ m + r ′ ⊤ m v den + c den = N um Den ( 10
)
where,

M den =P ⊤ n+1 W ⊤ n+1 (D n+1 + σ 2 e I n+1 ) -1 P n+1 v den =P ⊤ n+1 W ⊤ n+1 (D n+1 + σ 2 e I n+1 ) -1 (f n+1 + e n+1 -x 0 ) +P ⊤ n+1 (D n+1 + σ 2 e I n+1 ) -1 W n+1 (f n+1 + e n+1 ) c den = (f n+1 + e n+1 -x 0 ) ⊤ (D n+1 + σ 2 e I n+1 ) -1 W n+1 (f n+1 + e n+1 ) M num =P ⊤ n+1 W ⊤ n+1 (D n+1 + σ 2 e I n+1 ) -1 W n+1 P n+1 v num =P ⊤ n+1 W ⊤ n+1 ((D n+1 + σ 2 e I n+1 ) -1 + (D n+1 + σ 2 e I n+1 ) -1 ⊤ ) W ⊤ n+1 (f n+1 + e n+1 ) c num =W ⊤ n+1 (f n+1 + e n+1 ) ⊤ (D n+1 + σ 2 e I n+1 ) -1 W n+1 (f n+1 + e n+1 )
and x 0 is the n + 1 vector [x 0 , ..., x 0 ] ⊤ .

As the parameter estimation expression [START_REF] Box | A basis for the selection of a response surface design[END_REF] is derived, the criteria derivation is accomplished by taking the expectation of (τ n+1τ * )

2 with respect to the model and observation errors as shown in [START_REF] Wynn | Results in the theory and construction of d-optimum experimental designs[END_REF]. At this step, as τ * is unknown, it is replaced by its estimation τn . This is the purpose of the two-stage sequential DOE presented in the introduction. Expression [START_REF] Box | A basis for the selection of a response surface design[END_REF] can be written as follows:

τn+1 -τ * ≈ r ′ ⊤ m Mr ′ m + r ′ ⊤ m v + c r ′ ⊤ m M den r ′ m + r ′ ⊤ m v den + c den = A Den (11)
where, M = M numτn M den , v = v numτn v den and c = c numτn c den . Doing statistics on ( 11) is not evident because of the presence of stochastic variables in the numerator and denominator. Fortunately, for reasonable choice of the model error variance s 2 , the variations of

ǫ den = r ′ ⊤ m M den r ′ m + r ′ ⊤ m v den are negligible compared to c den : A Den = ǫ A + c ǫ den + c den ≈ c c den + 1 c den ǫ A - c c 2 den ǫ den
Furthermore, it was noticed that the variations of c c 2 den ǫ den are negligible compared to the variations of 1 c den ǫ A . Therefore, the denominator is supposed to be constant (equal to c den ). This point is illustrated in the next section. Now, doing statistics on [START_REF] Montepiedra | Minimum bias designs with constraints[END_REF] is possible because Den is constant and one has to do statistics only on A. This justifies our suggestion of working with τ * instead of a * . Knowing that Den is constant, the criterion in ( 3) is approximated by:

J design (t, τn ) = 1 Den E {en,en+1,r} A 2
The stochastic variable A = r ′ ⊤ m Mr ′ m + r ′ ⊤ m v + c consists of the products of correlated Gaussian variables. Therefore, it may be computed thanks to a singular value decomposition:

A = z ⊤ S z z + z ⊤ v z + c z
where c z is a constant value, z ∼ N (0, I n+1 ), and S z is a diagonal matrix.

The mean m A and the variance σ 2 A of A can be calculated as functions of S z , v z and c z :

m A = k λ k + c z σ 2 A =2 k λ 2 k + k λ k 2 + k v 2 k
where λ k is the k-th element in the diagonal of S z and v k is the k-th element of v z . The expression of J design is thus tractable:

J design (t, τn ) = 1 Den m 2 A + σ 2 A

ILLUSTRATIVE EXAMPLE

In this section, an illustrative example is given in order to show the advantage of performing DOE (comparison of designs reached by a sequential DOE approach and other designs reached with the classical uniform criterion) and the advantage of considering a model error (comparison of the proposed design with the non-robust design obtained by taking D n = 0 in the estimation stages and r n = 0 in the design stages).

In this example, the control input is taken to be u(t) = 0 and the system is only allowed for the initial condition x 0 = 1.

The Gaussian kernel is used because it is the most used kernel for the Gaussian process covariance [START_REF] Roger | A criterion for model-robust design of experiments[END_REF]:

c(t -t ′ ) = s 2 exp - t -t ′ λ 2 , ∀(t, t ′ ) ∈ [0, t max ] (12)
where, s 2 (Gaussian process variance) and λ (correlation distance) are the Gaussian process parameters. The kernel is used with the Gaussian process parameters values of s 2 = 0.01 and λ = 0.2. The approach used to choose these values is based on a maximin efficiency criterion (see [START_REF] Roger | A criterion for model-robust design of experiments[END_REF][START_REF] Abiad | Model-Robust Sequential Design of Experiments for Identification Problems[END_REF] for more details).

Having a set of initial points between the initial time value (t 0 = 0) and the intermediate time value (t int ), the goal is to select the time values between t int and t max in order to best estimate the parameter τ * . Let t max = 30 s the final time, t int = 1 s and ξ 3 = [0.33, 0.66, 1] ⊤ the initial design vector. The initial parameter value a * = -0.5 and the time constant corresponding value is τ * = 2 s.

The proposed approach has been tested with a Monte-Carlo method where 200 identification problems were created. Each problem is defined by a particular realization of the Gaussian process r(t) and by a particular realization of the observation error e n with σ 2 e = 0.05. In section 3, an assumption about Den was made. We have assumed that the variations of r ′ ⊤ m M den r ′ + r ′ ⊤ m v den are negligible in front of the constant term c den . Therefore, Den was considered to be constant Den ≈ c den . Figure 2 Figure 3 shows the mean of the square errors (τ n+1τ * ) 2 when 30 design points are added. The comparison is done among the proposed sequential approach, the nonrobust version of the sequential approach and the uniform approach where the design points are distributed uniformly over the experimental domain. The results show a faster Another way of comparison is to compare the number of added design points to reach a square error lower than a fixed threshold. Figure 4 compares the histograms of the added design points for the three criteria. The corresponding means are shown in Table 1. It is clear that the proposed approach is the most efficient because it needs the smallest number of added points. This result may be particularly interesting when the number of experiments has to be small.

CONCLUSION AND FUTURE WORK

This paper has discussed the idea of designing robust experiments for dynamic systems. The main contribution of this paper is the idea of modeling the misspecification by a statistical representation (Gaussian process). The paper has presented a detailed derivation of the proposed criterion showing the simplifications and main contributions. Finally, an illustrative example was presented showing the improvement in the speed of parameter convergence over other criteria.

The proposed idea is tested on a very simple context since it is restricted to first-order linear ODE. The goal of the paper is to show the relevance of using a statistical representation associated with a statistical criterion in order to reach robust estimators. Further works on more complex systems (non-scalar and non-linear ODE) are in progress. We also plan to exploit the proposed criterion to design for both t n and the value of the control input.

Fig. 1 .

 1 Fig. 1. Sequential DOE In the estimation stage, the unknown parameter τ * is estimated using the n collected observations: τn = arg min τ J estim. (τ, ξ n , y n )

  r n = t1 0 r(θ)dθ, ..., tn 0 r(θ)dθ ⊤ is unknown and its possible values depend on the used statistical representation.

2 denǫFig. 2 .

 22 Fig. 2. Histograms of Den and ǫ A × c den c

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison among the three criteria

Table 1 .

 1 Means of the number of added design points

		Robust Non-robust Uniform
	Mean Added Points	3.69	7.90	11.88


domain: sdv sde

text:
and tolerance to water deficit Control of root growth under water deficit ?

Roots Roots

-anchor plants -uptake water and nutrients

Under drought, growth of roots -allows the prospection of distant and wetter soil horizons -contributes to plant acclimation

Root growth under osmotic stress

osmotic stress (PEG 3350 g mol -1 ) :

-Moderate stress :

(90 mosmol kg -1 , 100g L -1 PEG)

Time (hour)

?

Distance from root tip (mm)

0 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 0.00 0.05 t 2
Mature organ :

? Dilution Temporal Convection Net accumulation ( ) ( ) ( ) ( ) x x x x x v ρ x ρ v t ρ D ∂ ∂ + ∂ ∂ + ∂ ∂ =

  photography

  The continuity equation (issued from a fluid mechanics formalism)gives access to the net accumulation rate of transcript, i.eExpression of aquaporins in the root apex -Root Symposium -August 2011 -16Distance from root tip (mm) in the root apex -Root Symposium -August 2011 -18Time before end of expansion (h)


Jacques Chauvin

Patrick Weber

Jean-Pierre Aubry

Frederic Lefebvre

Fabrice Sthal

Serge Galliou

Enrico Rubiola

Xavier Vacheret

A new generation of veryhigh stability BVA oscillators

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The main features targeted of that new "8607-C series" are to get significant improvements compared to the classical "state of the art" 8607-B design in a better short term stability @ 1 sec in Allan variance, a better-low phase noise and outstanding short term stability and a better-high isolation from "pressure and humidity" variations.

1) INTRODUCTION: MAIN OBJECTIVES

BVA oscillators are mostly used as local oscillators in ground atomic clocks (Cs fountain, primary references,..). Then the main goal in this application is a phase noise lower than -130 dBc/Hz @ 1 Hz offset, and sigma tau lower than 8 10 -14 @ 1 s, 5 10 -14 on floor up to 30s. For metrology purpose, our target is to get a phase noise floor better than -160 dBc / Hz and thermal stability better than 1.10

-10 in -15°C to + 60 °C temperature range. For specific application, medium term stability is a strong requirement. In DORIS application, the goal is to get, from a set of 90 samples (τ =10s) an average slope (linear regression) better than 2 10 -13 per min, with a residual noise (distribution around the average slope) lower than 1 10 -13 . In such case, ageing thermal stability and Allan variance must be at their best level.

2) THERMO-MECHANICAL STRUCTURE

2-1) GENERAL PRINCIPLE

In order to satisfy new exigencies expressed in §1, we have chosen an original thermo-mechanical structure with the main following characteristics, described on next graph: a) A double oven structure placed at the neck and inside a Dewar-glass a) The "internal oven" housing both BVA resonator and oscillator circuit hermetically sealed in order to reduce "pressure & humidity" sensitivity c) The "internal oven" is mechanically realized in a heavy copper block to reduce the thermal transients d) The Dewar-glass neck is closed by a temperature controlled copper plate (external Oven) e) Both "Ovens" are temperature controlled by 2 separated P.I loops circuits (proportional-integral type) f) A rigid central "composite Beam" links the two ovens together and fixes the double oven structure on the external case g) Due to the low thermal conductivity coefficient ( §2.4) of the beam material, there is no major thermal link between the copper plate (at the neck of the Dewar-glass) and the internal oven (at the extremity of the beam)

2.2) INTERNAL OVEN

The "internal oven" is the most important piece regarding thermal performances. Its main features are: -Perfect thermal homogeneity & high thermal inertia -Hermetic enclosure & electrical links insured by glass-beads -Heating ensured by power transistors -Must support a connexion PCB to drive the heating transistors Machined in a heavy copper block, the "internal oven" mass is # 370g, providing a thermal time constant of the internal oven around 2000s.

2.3) EXTERNAL OVEN

The "external oven" function is mainly filtering the external temperature variations by closing the "Dewarglass neck" by a quasi-isothermal cape. In these conditions we have chosen a simple circular copper plate heated by a single power MOSFET transistor. Added to "Dewar-glass" natural insulation properties and "Internal oven" large thermal inertia, the significant reduction of the temperature variation magnitude seen by the internal oven must allow to fulfil "static" and "dynamic" high performances under temperature gradient.

2.4) CENTRAL COMPOSITE BEAM

The first function of the central composite beam is to fix strongly both "external" and "internal" ovens to external structure. The second function is to allow a high thermal insulation between the 2 ovens to avoid thermal leakage. Two physical properties of the selected material are of main interest in our application:

-The tensile modulus of elasticity: 18000 N/mm2 -The thermal conductivity: 0.30 W/K.m

2.5) COMPLETE THERMO-MECHANICAL STRUCTURE

The complete thermo-mechanical structure is presented on Fig1.

3) DOUBLE OVEN THERMAL REGULATION

Both ovens are "temperature controlled" by two separated P.I loops (proportional-integral). "Heating modules" are realized with MOSFET power transistors and "Temperature sensors" are "glass sealed" precision small size thermistors.

3.1) GENERAL P.I.D TRANSFER FUNCTION

The fig 2 and fig 3 show a general P.I.D (Proportional-Integral-Derivative) electrical scheme equivalent to our thermal loop and its equivalent transmittance. The internal oven is considered as equivalent to a "first order low pass filter" associated to a pure delay time (∆t # 1s) with a very high time constant (τ # 2000s). The thermistor assembly is also equivalent to a "low pass filter" but with a small time constant (τ capt # 6s). In a first approximation we can consider than the two "first order" filters are connected in tandem. In our specific application, the "derivative function" will not be connected because not useful regarding the 8607's very long stabilisation times. Fig 2 In order to avoid frequency transients when the external power supply voltage is perturbed, the internal oven power chain is powered by a regulated "power supply" placed outside the Dewar-glass. Both Gain and Phase of equivalent "open loop thermal transmittance" are shown on fig 5. The internal oven temperature reference will be adjusted at "ToT" (Resonator Turn over Temperature). For example, the average "ToT" for our BVA / 5MHz / SC cut (overtone 3) resonators is # 85°C. The "internal oven" real thermal gain (ratio between external temperature deviation & internal temperature deviation), must be ≥ 2500. 1262

P.I OPEN LOOP TRANSMITANCE

3.3) EXTERNAL OVEN P.I THERMAL LOOP

The same electrical scheme is used for "external oven" but with only one heating MOSFET transistor (instead of two in the internal oven).

The copper plate objective is to filter the external thermal perturbations at the neck of the "Dewarglass". The copper plate reference temperature "Tcp" will be given by the following expression: If Top max = max operating temperature, the best solution is to have Tcp > Top max in order to keep the 2 thermostats both operating even if the external temperature is at the maximum of the operating temperature range.

Tcp ≤ ToT -(Rth x
In our particular case, the internal temperature increasing due to internal dissipative power (Rth x Pint) is about 10°C" For example, when the max operating temperature is 50°C, then "Tcp" will be fixed ≥ 60°C. Concerning the "external oven" real thermal gain and contrarily to § 3.2, this point is not critical and a low value is widely sufficient (e.g: ≥ 30)

Remark:

The combination of "external oven" and "internal oven" placed in "tandem" must insure a total real thermal gain ≥ 75000. E.g.: the maximum temperature variation measured on BVA resonator and oscillator circuit will be ≤ 0.001°C when external temperature is varying from -20°C to + 50°C.

4) OSCILLATOR & BUFFER TOPOLOGY

4.1) GENERAL PRINCIPLE

We have chosen a classical feedback transmission "Colpitts type" oscillator allowing a good optimisation of all parameters able to affect the oscillator's noise performances.

In particular, the selected configuration is well adapted to adjust: -The BVA resonator excitation level -The "Input" and "Output" load impedances seen by the resonator -The oscillator bandwidth "High Input impedance" Output Buffer allows keeping a large S/N ratio. The general principle scheme is given in fig 6.

In our particular application Zx in # Zx out # 20Ω; Zb in ≥ 2KΩ; Zb out = 50Ω

4.2) "ADS" SIMULATION RESULTS

The computations results obtained with ADS software (Advanced Design Software) are shown in the table hereunder.

In that simulation configuration, the intrinsic BVA crystal noise isn't taken into account but the transistor noise factor is represented by two different values in order to view the consequences on oscillator's phase noise spectrum. As seen in that table, the main differences are only obvious below 10Hz.

4.3) LEESON'S MODEL

The Leeson's equation takes into account significant parameters [see remark hereunder] that determine the oscillator's single-sided phase noise density (including the "flicker corner" (fc) of the active component) like expressed hereunder:

                +         + + = m c l m m c m f f Q f f f f Pavs FKT Log f L 1 2 1 2 10 ) (

Remark: Even if "fc" term can be considered as large enough and general, it is not obvious to integer in this model some added noise (proper to XO's) like the intrinsic crystal's "flicker noise".

Some other "deterministic" additional noise types caused by external phenomenon's like:

-The OCXO's thermal and pressure sensitivity -The dynamic thermal behaviour proper to the "Crystal Cut" (AT or SC) are not taken into account. In these conditions, the real measurements must be theoretically worst than those given by the Leeson's model especially for long integration times (e.g. ≤ 0.1Hz offset frequencies). For information, the first "8607-C" prototype's real measurements have been obtained between N2-A and N3-A units.

It is easy to see that the N2-A real measurements are in good correlation with the ADS simulations results above 2Hz offset, but the model becomes not correct below 2Hz. The little difference on the floor is only due to the excitation level difference.

The comparison with Leeson's model shows a good correlation with the OSA reference144 (8607-B previous design) up to 100Hz even if the low excitation level applied in that oscillator can explain a significant floor level difference. An important gap can be observed with the N2-A prototype in the 2Hz---1000Hz offset range where the real results are significantly better showing a clear improvement. Another remark is the perfect correlation between Leeson's model (intending a probably non realistic fc=280Hz ?...) and our N2-A measurements in the 10Hz---200KHz bandwidth.

Below 10Hz, the apparent limitation could be mainly attributed to crystal's "intrinsic flicker noise" (up to 5dB below 2Hz).

The most significant point can be expressed by the fact that below 2Hz, there is no difference between our new prototypes (N2-A) and the "previous version" reference (Nr144) showing probably the intrinsic crystal's noise limit.

The Leeson's model, for which this crystal intrinsic noise isn't taken into account, shows better results (up to -5dB between 0.1Hz and 2Hz).

The particular point -134dBc/Hz @ 1Hz could be considered as a good target if we want to valid a new low noise BVA resonator (progres in design or manufacturing process).

4.5) FIRST PROTOTYPES NOISE; A PHYSICAL TENTATIVE INTERPRETATION

According to noise theory it is easy to write the noise frequency spectrum as the sum of some physical types:

α α f h f S y ∑ - = 0 3 ) (
A well known relationship allows expressing the phase noise spectrum as:

α α ν ν f h f f S f f S y 2 2 0 3 2 2 ) ( ) ( ∑ - Φ = =
And it is possible simplifying by using α b coefficients as: According to E.Rubiola theory [ref 2] a fine graphical analyse shows some interesting improvements in term of flicker phase noise (see table hereunder).

α α f f b f S ). ( ) (
According to graphical observations we can deduct the main physical characteristics: -The apparent loaded Q is relatively high # 1.5E6 (Ql # Qx60% according to design) -The sustaining amplifier "Noise Factor" is # 2.2 (3.5dB) -The parametric "flicker noise" contribution specially generated by Output Buffers is # 6dB with probably a 3dB potential improvement margin.

-The intrinsic "flicker noise" contribution generated by resonator only can be estimated at 4dB. -The potential electronic performance for such oscillator without any other design improvements is # -134dBc/Hz @ 1Hz offset (in compliance with §4.4).

-The ultimate electronic performance (assuming -3dB flicker noise Output Buffers improvements) can be estimated # -137dBc/Hz @ 1Hz offset (in compliance with ADS simulations).

5) CHOICE OF BVA RESONATOR'S DESIGN

In many oscillators phase noise measurements, the results obtained by exchanging only "electrically equivalent resonators" can show significant differences (up to 10dB @ 1Hz). The same observation has been made at Femto-ST Institute with the passive method ( §6) on resonators only. This problem widely discussed in some theoretical works [ref 2] is probably due to resonator's "intrinsic flicker noise" generally not known and not measured before utilization in oscillators.

In these conditions, an important work has been started to optimise "BVA design or manufacturing process parameters" able to reduce the crystal's intrinsic flicker noise. In parallel, measurements methods allowing knowing quickly resonator's intrinsic noise performances (before utilization in oscillators) must be developed to evaluate "design or process" progress.

In order to get directions to discussion's supports, we have compared in the scope of this work two different BVA SC 3rd Overtone resonators. The idea behind is, one hand to optimise design maximising the Q value (assuming that the 1/Q 4 law derived by F Walls, JJ Gagnepain & all.. [Ref 3] is still valid), and on the other hand to get resonators with similar Q while slightly different C1 parameter, in order trying to quantify the electronic contribution.

6.1) SHORT DESCRIPTION OF THE PASSIVE METHOD

The intrinsic flicker frequency noise of quartz crystal resonators can be measured by means of passive methods. In this case, the resonator noise is observed in term of phase noise without the noise usually associated with an active oscillator. In 1975, the first technique using passive measuring system was developed by Walls and Wainwright [Ref 4]. We will call it "simple bridge method". Two crystals as identical as possible are driven through a π transmission network by a unique low-noise source and amplified. Then, both signals are mixed with 90° phase difference by a double balanced mixer (DBM). The DBM output, used as a phase detector produces a DC voltage which is proportional to the instantaneous phase difference between the two signals. A good balancing of each arm in term of resonator loaded Q's and resonant frequencies gives a DBM output signal insensitive to the source noise. The reduction of the source noise authorizes the noise detection of both resonators. This kind of measurement system is well adapted for numbers of resonators but the noise floor of the system is limited about -140 dBc/Hz @ 1 Hz carrier offset. This is due to the added flicker noise of the amplifier needed to increase the resonator signals before the DBM. Thus, the detection of inherent resonator noise of the best crystal is not possible (Pair measurement). Another technique is now available to solve this problem. At the end of nineties, crystal resonators testers were designed to assist in the PM noise characterization of quartz crystal resonators in the 1 to 200 MHz domain. We will call this method "carrier suppression technique". Recently, FEMTO-ST Institute has developed an improved carrier suppression system in order to find the origin of intrinsic flicker noise of acoustic wave resonators [Ref 9].

6.2) MEASUREMENTS PRINCIPLE

Fig. 10 shows the principle of the carrier suppression technique.

Figure 10: Carrier suppression principle.

The DUT's bridge is quasi identical to the simple bridge method. The carrier signal of the source is splitted into equal parts to drive both devices under test (DUT). The DUT's can be resistors to measure the noise floor of the system or crystal resonator pairs to measure their inherent phase noise. The resonant frequency of each arm of the bridge is tuned to the source frequency with a serial capacitor. The difference between both methods is that the crystal output signals are not 90° mixed together but combined 180° out of phase. In this case, the carrier signal is subtracted. Since phase noise is defined relative to the carrier power, reducing the carrier level has the effect of amplifying the phase noise of the DUT. The combiner output signal is increase by about 60 dB by the HF amplifier and then mixed by the phase noise detector with a 90° phase shift parts of the carrier signal. Then, the LF amplifier pushes the output signal to a level compatible with the Fast Fourier Transform (FFT) analyzer. Thus, the carrier suppression technique gives a better gain than the simple bridge method. Moreover the flicker noise of the HF amplifier is very low because of the low carrier level. Calibration of the measurement system is obtained by injecting a known side band signal or a known white noise source (KNS) on one of the arms of the bridge. The noise of the DUT, as seen on the FFT analyzer, is corrected using the calibration factor determined with the side band or the KNS. The resonator noise is given through the single sideband power spectral density of phase fluctuations L(f). The measured noise is attributed to both resonators if they can be considered identical S φ (1 Hz) = L(f). If one resonator noise is known and significantly better, the measured noise can be attributed to only one resonator S φ (1 Hz) = L(f) + 3dB. S φ (1 Hz) is converted in Allan standard deviation by means of loaded quality factor QL of the resonator.

) 1 ( 2 1 ) 2 ln( 2 ) ( 2 Hz S Q L y φ τ σ         = QL = ν0 /2fL is computed from the Leeson cut-off frequency fL. f L represent the f
-1 to f -3 slope change of the L(f) curve. But, more precision is usually given on fL by a measure from the transfer function of the resonator obtained with the KNS source. In the case of Fig. 11, fL is equal to 1.8 Hz Thus, if the BVA 4304 029 36 is considered as the reference, we have:

13 10 / 130 2 6 10 34 . 1 10 10 5 8 . 1 ) 2 ln( 2 ) ( - - ⋅ =       ⋅ = τ σ y
The Allan standard deviation calculated this way represents the Allan standard deviation of an oscillator using this resonator in which the only source of flicker frequency noise is the resonator under test. Thus, a correspondence between this value and the classical measurement of the oscillator flicker floor gives the proof that the resonator noise limits the oscillator.

6.3) RESULTS with 5MHz "STANDARDS" BVA

The measurements results (table hereunder) are obtained with "passive method" @ 60µW by using BVA 4304 029 36 as noise reference (no correction). These homogeneous results can be considered as representative of the standard 5MHz BVA design.

6.4) RESULTS with 5MHz "SPECIFIC" BVA

The measurements results (Table hereunder) are obtained with "passive method" @ 60µW by using BVA 4304 035 27 as noise reference (no correction). Other tests are in progress @ 100µW to evaluate the excitation level impact. The measurements results show the phase noise improvement obtained with the new design compared to the actual state of the art. Despite the values measured at one Hz, around -128 dBc/ Hz, which are a little bit disappointing regarding the Allan variance target at 1 Hz, we can notice the significant improvement, # -4 dB, around 10 and 100 Hz.

The noise floor improvement was expected by design, but the flicker contribution improvement is a prime result.

7.2) ALLAN VARIANCE DEVIATION (5 MHz

Standards BVA) "Allan variance deviation" measurements are issued from a "5110A Time Interval Analyser" (Timing Solutions).

Remark:

The results presented hereunder are given as the sum of two pieces. A correction factor equivalent @ 0.707 can be applied in the specific case where the two pieces can be considered as strictly equivalent.

In the first line, the grey case (

9) CONCLUSION

The results obtained with the "New Design 8607-C" first prototypes shows that our oscillator's electrical scheme allows significant progress in term of phase noise above 1Hz offset frequencies (up to -5dBc/Hz @ 100Hz compared to previous 8607-B design). Below 1Hz, the noise performances are mainly given by the resonator. The important works started in collaboration with Femto-ST Institute show a good correlation between results obtained on the resonator only (passive method) then on oscillator (active method). The phase noise measurements differences, always included in the "error margin", show a "not significant" contribution of the oscillator's electronic (at least very close of the "Passive Method" residual contribution).

In order to still progress and reach noise performances widely better than -130dBc/Hz @ 1Hz (objective §4.4: -134dBc/Hz), an important effort must be done in term of resonator's design or manufacturing process to reduce their intrinsic flicker noise. Even if all the works aren't yet finished, the first "noise measurements" obtained on pieces 4304 035 27 & 4304 035 28 ( §6.4) with the last "Specific BVA" resonators are a stimulating factor which merits to be continued.

Another interesting result to be signalled concerns our last series of standard BVA resonators for which a particular manufacturing process has been applied. A first outcome of this work is the successful comparison between passive and active method of noise characterisation even if the results are still obtained with a little number of resonators.

A second outcome of this work is the good sensitivity of the passive method (still to be confirmed on oscillator) showing a remarkable "6 10 -14 @ 1s" Allan variance.

A third main outcome is the successful attempt to drastically reduce the flicker noise on the active oscillator, between 10Hz and 1KHz with an electronic flicker noise limitation estimated at -134dBc/Hz @ 1Hz offset. A new "tests campaign" with the most performing resonators recently characterised will be started to approach that potential limit.

The excellent results obtained with 2 specific BVA resonators (pieces 4304 035 27 & 4304 035 28) and with the last pieces of series 34 show that it's possible to reduce significantly "intrinsic noise" by design and manufacturing process even if the physical mechanisms behind aren't yet known and if the reproducibility is still far to be insured. In the next steps, understanding the crystal's intrinsic noise physical mechanism could be an interesting subject to be studied with the funding support of Femto-ST institute [ Ref 5].

Concerning the oscillator and buffer noise contributions, E.Rubiola theoretical and graphical approach seems a good way understanding and identifying noise sources in order improving electronic performances. Some other works concerning thermal and humidity characterisation remain to be done before starting the pre-serial phase.

A

  NEW GENERATION OF VERY HIGH STABILITY BVA OSCILLATORS Dr Jacques Chauvin*, Patrick Weber*, Dr Jean-Pierre Aubry*, Frederic Lefebvre*, Dr Fabrice Sthal**, Dr Serge Galliou **; Dr Enrico Rubiola**, Xavier Vacheret** * Oscilloquartz S.A, rue des Brévards 16, CH-2002 Neuchâtel -Switzerland ** Femto-ST Institute, Chemin de l'Epitaphe 26, 25030 Besançon Cedex -France chauvin@oscilloquartz.com SUMMARY A third generation of "Oscilloquartz" OCXO's using the technique of housing a BVA SC-cut crystal resonator and its associated oscillator components in double oven technology has been developed with the funding support of European Space Operations Centre (E.S.O.C). The main purpose is to provide a local oscillator for high performances ground clock [ref 1].

Fig 3 3

 3 Fig 3 3.2) INTERNAL OVEN P.I THERMAL LOOP Fig 4 shows the electrical scheme of the "internal oven P.I thermal loop".

Fig 5

 5 Fig 5

  ) = Phase noise (in dBc/Hz) Ql = Loaded Q fm = Carrier offset frequency (Hz) fo = Carrier center frequency (Hz) fc = Flicker corner frequency of the active device (Hz) T = Temperature (°K) Pavs = Average power trough the resonator (W) F = Noise factor of the active device (F=Nout/(G*Nin)) K = Boltzman constant (J/°K)

Fig 9

 9 Fig 9 According to Fig 9, our "electronic configuration" seems highly performing compared to previous "8607-B" version by reducing drastically the flicker noise in the 2Hz---1000Hz bandwidth.

=

  In these conditions a simple measurement of α f intercept point with the 1Hz axis gives directly the α b value.

  Fig. 11 gives an example of L(f) measurement @ 60 µW for BVA 4304 029 36 and BVA 4304 027 47. L(1Hz) = -133 dBc/Hz

Figure 11 :

 11 Figure 11: L(f) @ 60 µW of BVA 4304 029 36 and BVA 4304 027 47.

5 MHz SSB Phase Noise (fc=1000Hz; NF=3.5dB)

					In our particular case, it is necessary to reduce
					fc@280Hz (instead of 1KHz ?...) to obtain a good
					correlation with our prototype's real measurements,
					showing a significant "flicker noise" reduction in the
					"2Hz---1000Hz" bandwidth.
					Even if "fc" can't be considered as a relevant physical
					parameter [ref 2], it's obvious that such theoretical
					value seems not physically realistic.
					To view this questioning, the Fig 9 shows the results
					issued from:
					-LEESON's model computations (but with fc reduced
					to 280Hz)
					-N2-A prototype (8607-C) real measurements
					-Nr 144 5MHz OSA reference (previous design).
					In that particular case, according to previous remark,
					all the parameter's values applied in Leeson's model
					remain the same except "fc" (280Hz):
					5 MHz SSB Phase Noise (fc=280Hz; NF=3.5dB)
					Leeson model	ADS simulation	N2-A prototype	OSA Ref 144
					-90
					-100
					-110
					Lf (dBc/Hz)	-140 -130 -120
					-150
					-160
					-170
					1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06
					Fourier offset frequencies (Hz)
		Leeson model	ADS simulation	N2-A prototype	OSA Ref 144
		-90
		-100
		-110
	Lf (dBc/Hz)	-140 -130 -120
		-150
		-160
		-170
		1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+05 1.E+06
			Fourier offset frequencies (Hz)
	Fig 8

Remark:

Regarding the Leeson's model, the most significant parameter able to change the SSB Phase Noise shape is the flicker noise level given by "fc" (oscillating transistor "flicker corner").

  These units use carrier suppression based on the bridge technique [Ref 6-8]. The noise floor of these systems obtained by means of resistors is # -155 dBc/Hz for a 70 µW carrier power.

	Noise Type	Freq slope	bα	Ref 144 SФ(1Hz)	N2-A proto. SФ(1Hz)
				dB	rad2 /Hz	dB	rad2 /Hz
	P.white	f0	b0	-155	3.2 e-16	-161	8 e-17
	P.flicker	f-1	b-1	-131.5	7.1 e-14	-136.5	2.2 e-14
	F.flicker	f-3	b-2	-127.5	1.8 e-13	-127.5	1.8 e-13

dBc / Hz) for 1 piece (-3dB correction)

					BVA /4304		BVA /4304	BVA /4304
					035 29		035 30	035 31
	Lf(dBc/Hz@1Hz)	-136		-123		-129	-133
	SФ(dBrad^2/Hz)	-133		-120		-126	-130
	Q(10^6)		2.7		2.7		2.7	2.7
	Ql (10^6)		1.5		1.56		1.56	1.45
	Ql / Q		55.44%	57.87%		57.87%	53.83%
	SФ(rad^2/Hz)		5E-14		1E-12		2.5E-13	1E-13
	σy(-ε) (E-13)		0.65		2.88		1.41	0.958
	σy(1s) (E-13)		0.88		3.77		1.89	1.28
	σy(+ε) (E-13)		1.17		5.04		2.53	1.71
	As opposed to previous standard design ( §6.3), these
	first "specific" BVA show an important dispersion in
	term of phase noise (e.g.: up to 13dB @ 1Hz between
	Nr 28 and 29 pieces).
	That significant difference isn't explainable by the
	quasi-identical motional parameters and seems
	demonstrate a wide intrinsic noise level dispersion.
	Another interesting point is the very low noise level
	measured on resonators Nr 4304 035 27 & 4304 035
	28 (σy = 8.8E-14).
	By applying a "0.707" correction factor (assuming that
	the 2 pieces are strictly identical), the average Allan
	variance @1s for one piece is remarkable # 6.23E-14
	(+/-2E-14)!
	7) FIRST OSCILLATORS PROTOTYPES RESULTS
	(ACTIVE METHOD)
	7.1) PHASE NOISE RESULTS (5 MHz Standards
	BVA)
	In order to avoid any "error risks" particularly for very
	low offset frequencies (inside the PLL bandwidth) all
	phase noise measurements are issued from 2
	different systems:
	In Frequency domain (>1Hz): PN9000-Phase noise
	measurement system /Aeroflex
	In Time domain (>1s): 5110A-Time interval analyser
	system / Timing Solutions
	Remark: The "References pieces" Nr 144 & 303 are
	5MHz "noise references" in use @ "Oscilloquartz S.A"
	(previous 8607-B type).
	Pieces "N1-A2", "N2-A" and "N3-A" are the 3 first
	prototypes @ 5MHz ("New Generation" 8607-C Type).
	PROTOTYPES & REFERENCES PHASE NOISE PERFORMANCES
	Values ( DOMAIN TIME (5110A)		SPECTRAL (PN9000)
	Fourier	0.1	1	10	100	1	10	100	1
	Frequencies	Hz	Hz	Hz	Hz	KHz	KHz	KHz	MHz
	OSA Ref. ( 8607-B) 144 or 303	-102	-129	-145	-153	-155	-157	-158	-161
	N1-A2 (8607-C)	-97	-127	-148	-158	-162	-164	-164	-166
	N2-A (8607-C)	-103	-129	-149	-158	-162	-164	-164	-167
	N3-A (8607-C)	-100	-128	-149	-158	-162	-164	-164	-167

ALLAN variance" for 1 piece (E-13)

							3.59 N.C) shows a "not
	correct" measurement caused by a beat frequency
	pollution problem.
			PROTOTYPES & REFERENCES SHORT TERM STABILITY Allan Variance (10-13) for 2 pieces (no correction)
	Time τ (s)	1	2	4	10	20	40	100	200	400	1k
	Ref 144 / 303	3.59 N.C	1.81	1.18	1.04	0.99	1.12	1.40	1.67	2.22	3.7
	144 / N1A2	2.80	2.18	2.02	2.15	2.49	3.23	4.01	5.3	7.3	13.2
	144 / N2-A	1.84	1.74	1.91	2.59	1.41	1.48	1.86	2.64	4.0	7.8
	144 / N3-A	2.29	2.02	2.06	3.01	3.28	2.19	2.96	3.6	4.8	10.1
	8) COMPARISON BETWEEN PASSIVE AND
	ACTIVE METHOD
	We have reprinted on the same table the results
	obtained with the same BVA resonators measured
	with the two methods (Passive and Active).
	Except for "N1-A prototype" where instabilities seen on
	oscillator get the measurements not significant, a good
	correlation can be observed between the two methods
	on the 3 other pieces.
		SPECTRAL DOMAIN				TIME DOMAIN
		PASSIVE METHOD				ACTIVE METHOD
	(Resonator Nr ; Sigma@1s)	(Osc. Nr ; Sigma @1s ; Best Sigma)
	BVA 4304 029 36	≤ 1.34	N1-A			N.A (Instabilities)
	BVA 4304 025 14	1.66 (+/-0.4)	N1-A2		1.97	1.41(τ = 4s)
	BVA 4304 027 47	1.34 (+/-0.4)	N2-A		1.30	0.996 (τ = 20s)
	BVA 4304 027 57	1.75 (+/-0.4)	N3-A		1.61	1.50 (τ = 4s)

PASSIVE & ACTIVE METHOD COMPARISON

Average "

Variance (10-13) for 2 pieces (no correction)

  The following table shows exceptional "Allan variance deviations" measured on series 34 (active method): SERIES 34 SHORT TERM STABILITY (Active method) Allan

	Time τ (s)	1	2	4	10	20	40	100	200	400	1k
	144 / 460	2.27	1.41	1.07	0.96	0.94	1.02	1.23	1.38	1.61	2.40
	144 / 172	2.5	1.39	1.09	0.95	0.95	1.00	1.29	1.55	1.93	2.90
	144 / 472	2.27	1.45	1.14	1.02	1.05	1.11	1.19	1.37	1.66	1.62
	460 / 172	2.03	1.39	1.10	0.96	0.98	1.01	1.08	1.13	1.29	2.00

Remark:

In this batch the best results obtained between 10s and 20s averaging time are near of 6.8E-14 (average value for one piece after application of the 0.707 correction factor) showing a clear progress in the manufacturing process.

& 21th European Frequency and Time Forum, Jun 2007, Geneva, Switzerland. pp.1261-1268. hal-


est destinée au dépôt

Data entering, data management and QC Matieu Henry, Laurent Saint-André, Nicolas Picard To cite this version: Matieu Henry, Laurent Saint-André, Nicolas Picard. Data entering, data management and QC. National technical staff training on allometric equations (ae) under the UN-REDD Programme, Jun 2012, Hanoi, Vietnam. 7 p. hal-02805406 DATA FILTERING Field Data Sheet http://g-redd.cartong.org/aboutInventory.html Software http://nrs.fs.fed.us/inventory_monitoring/inventory/data_tools_services/ DATA STRUCTURE One line for one individual of the sampled population One column -> one descriptive variable

How data resulting from inventory and monitoring efforts will be acquired, managed, processed, checked, analyzed, distributed, and archived.

In order to avoid errors:

To check all the data several times: double entry


A new Physio-Demo-Genetics (PDG) model accounting for complex interactions among

genes, functional traits and climate when environment or demography are unstable -Account simultaneously for plasticity, adaptation and gene flow/migration -Provide short-time scale predictions (few generations), out of equilibrium -Use biophysical model : take climate as input, no assumption on trait optimum, trait-fitness relationship… -Both components of fitness decreased with increasing Tsum Results : strength of selection on tsum Francois Lefèvre  Aurore Bontemps, Maxime Cailleret, Julie Gaüzere  The field team of URFM (Nicolas, William) and UEFM (Norbert, Olivier, Frank)

FrostFrost

  1700 m on North face -840 and 1615 m on the South face Ecological data (growth, fructification, dispersal, phenotypic and genetic diversity) the studied case -No overlapping generations (sequential logging of reproductive trees) -No selection on date of budburst during the regeneration step -No mutation (only standing genetic variation + recombination) -No competition for light among adults -Climate (2002-2006 repeated) with elevation effect on T, Rainfall, and relative humidity Results: population dynamics Year simulated mortality at Alt2 → Alt4 -Expansion slower at Alt1 -Colonization with mortality at Alt5; Treeline at 1620 m (versus observed =1700 m) Results : validation Frost effect on LAI = 0 IUFRO meeting, 9-12/10/2012, Bordeaux -The model reproduces well patterns of variation of physiological and ecological variables with elevation.Results : spatial patterns of budburst



  offspring at next generation -Higher individual fitness in population Alt5 due to mortality combined with density-dependent selection -Strong selection gradient on Tsum in population Alt5 Back to the real world… Frost effect on LAI = 0 Tsum Population -In the case of beech, genetic cline shows counter-gradient variation (Tsum is lower for populations of high elevation) mortality at Alt2 → Alt4 -Expansion slower at Alt1 -Colonization with mortality at Alt5; Treeline at 1620 m (versus observed =1700 Population i + Generation j +ε i (+ interaction(Population i , Generation i ) )Results: variation of frost risk with elevation Frost effect on LAI = (Alt5) early genotypes were not more affected by late frost than late genotypes → early genotypes were still selected At middle elevation ( Alt 3-4) risks of late frost were maximum → early genotypes were counter-selected At low elevation, risks of late frost are lower→ early genotypes are still selected high gene flow, diversifying selection is capturing first beneficial allelic associations (θB) distributed among the loci contributing to the trait, prior to modifying allele frequencies → decoupling of Q ST and F STqInvestigation of this prediction with PDG modelWhat is the power of F ST -based methods (outlier approaches) to detect selection in such cases ? Mechanistic biophysical models → valid for other environment conditions  Modularity → adapted to other species 2. Evolution can be quick and strong  At high elevation Tsum decreased by 4° in 5 generations (mortality) François de Coligny  Christian Pichot, Philippe Dreyfus, Marianne Alleaume-Benharira &


Ali Chahine

Sylvain Dupont

Yves Brunet

Carole Sinfort

Numerical simulation of pesticide dispersal at the scale of rows and vineyards

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Grégory Thibault
email: gthibaul@bat710.univ-lyon1.fr

Alexandre Aussem
email: aaussem@univ-lyon1.frbonnevay@univ-lyon1.fr

Stéphane Bonnevay

Analyse critique des algorithmes EDA dans le cadre de l'apprentissage de structure de réseaux Bayésiens

Keywords: réseaux Bayésiens, apprentissage de structure, EDA, GES Bayesian networks, structure learning, EDA, GES

Récemment, une nouvelle classe d'algorithmes d'optimisation baptisée EDA (pour estimation of distribution algorithms) a montré son efficacité sur le problème de l'apprentissage de structure de réseaux Bayésiens. Pour cela, les chercheurs ont comparé les EDA avec des algorithmes génétiques. Dans cet article nous proposons de comparer les EDA avec un algorithme à base de score reconnu comme performant : GES. L'objectif de ce travail est de réaliser une analyse critique des EDA dans ce cadre.

ABSTRACT. Recently, a new class of optimization algorithms called EDA (stand for estimation of distribution algorithms) has shown its effectiveness on the problem of learning Bayesian networks structures. For this, the researchers compared the EDA with genetic algorithms. We propose here to compare the EDA with an algorithm-based scoring recognized as powerful: GES. The aim of this work is to achieve a critical analysis of EDA in this context.

Introduction

La recherche sur les réseaux Bayésiens a connu un grand essor ces dernières années, en particulier sur l'apprentissage automatique de la structure à partir de données [START_REF] Buntine | A guide to the literature on learning probabilistic networks from data[END_REF]. Ce problème très documenté est un problème NP-difficile [START_REF] Chickering | Learning Bayesian Networks is NP-Complete[END_REF].

Les EDA (ou Estimation of Distribution Algorithms) [START_REF] Blanco | Learning Bayesian networks in the space of structures by estimation of distribution algorithms[END_REF] sont des algorithmes qui ont fait leurs preuves dans de nombreux domaines [START_REF] Larrañaga | Estimation of Distribution Algorithms. A new tool for evolutionary computation[END_REF], y compris sur l'apprentissage de réseaux Bayésiens. Ainsi, [START_REF] Blanco | Learning Bayesian networks in the space of structures by estimation of distribution algorithms[END_REF], [START_REF] Romero | Learning Bayesian networks in the space of orderings with estimation of distribution algorithms[END_REF] et [START_REF] Thibault | Learning Bayesian network structures by estimation of distribution algorithms : An experimental analysis[END_REF] présentent plusieurs applications de ces algorithmes et montrent empiriquement leur efficacité en les comparant avec des approches génétiques [START_REF] Larrañaga | Structure Learning of Bayesian Networks by Genetic Algorithms : A Performance Analysis of Control Parameters[END_REF]. Dans une précédente étude [START_REF] Thibault | Learning Bayesian network structures by estimation of distribution algorithms : An experimental analysis[END_REF], nous avons réalisé une optimisation du paramétrage de l'algorithme afin d'en augmenter la performance.

Nous nous proposons ici d'évaluer la qualité des résultats de cette approche en comparant avec l'algorithme GES (ou Greedy Equivalent Search) [START_REF] Chickering | Optimal structure identification with greedy search[END_REF], un algorithme très performant utilisant l'espace des équivalents de Markov.

2. Apprentissage à base de score [START_REF] Heckerman | A Bayesian approach to causal discovery[END_REF] discutent de l'intérêt des méthodes à base de score. Elles permettent notamment l'incorporation aisée de connaissances a priori [START_REF] Cheng | Learning Bayesian networks from data : an information-theory based approach[END_REF], un travail plus efficace avec les données manquantes, et la possibilité de combiner plusieurs modèles pondérés (model averaging) [START_REF] Friedman | Being Bayesian about Network Structure[END_REF] rend l'étape d'inférence bien moins dépendante des erreurs éventuellement commises pendant l'apprentissage de la structure et permet de mixer des modèles provenant de plusieurs sources.

Scores

Les scores AIC (Akaike Information Criterion) [START_REF] Akaike | Statistical predictor identification[END_REF] et BIC (Bayesian Information Criterion) [START_REF] Schwarz | Estimating the dimension of a model[END_REF], deux scores dérivés du PML (Penalized Maximum Likelihood) sont composés de deux critères inconciliables : le maximum de vraisemblance (par rapport aux données) et la complexité du modèle. Les structures maximisant ces deux critères sont respectivement le graphe complet (au sens orienté), et le graphe vide.

P M L(B | D) = log L(D | θ M V , B) -f (N ) • Dim(B) [1] Où f (N ) = 1 pour AIC, log(N ) 2
pour BIC. Dans le cas de données discrètes finies 1 :

log L(D | θ M V , B) = n i=1 qi j=1 ri k=1 N ijk • log( N ijk N ij ) [2] Dim(B) = n i=1 (r i -1) • q i [3]
Les scores dits « Bayésiens », par exemple BD (Bayesian Dirichlet) [START_REF] Cooper | A bayesian method for the induction of probabilistic networks from data[END_REF], estiment la probabilité d'une structure a posteriori. Dans le cas des données statiques, on peut simplifier le calcul en utilisant la formule de Bayes :

p(B|D) = p(B, D) p(D) ∝ p(B, D) = θ L(D | θ, B) • p(θ | B) • dθ [4]
Une dérivation est possible lorsque les paramètres suivent une distribution de Dirichlet et que les données sont IID (indépendants et uniformément répartis) :

BD(B | D) = p(B) • n i=1 qi j=1 Γ(α ij ) Γ(N ij + α ij ) ri k=1 Γ(N ijk + α ijk ) Γ(α ijk ) [5]
Plus tard ce critère est adapté aux équivalents de Markov pour obtenir BDe [START_REF] Heckerman | Learning Bayesian networks : The combination of knowledge and statistical data[END_REF], puis avec des a priori uniformes, ce qui donne BDeu [START_REF] Buntine | Theory refinement on Bayesian networks[END_REF].

Il existe d'autres scores, comme MDL (Minimum Description Length) [START_REF] Bouckaert | Probabilistic network construction using the Minimum Description Length principle[END_REF] dont le principe est de minimiser la taille de codage nécessaire à la description du modèle et des données relativement au modèle. Notons enfin l'utilisation de scores adaptés aux recherches par arcs, i.e., qui utilisent des poids sur les arcs (au lieu des noeuds), comme l'information mutuelle [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]. D'une façon générale, le score définit l'objectif de la recherche : les critères PML sont utilisés lorsque le modèle recherché doit être efficace en inférence tout en minimisant sa complexité, les scores Bayésiens (BD, BDe etc.) favorisent les structures offrant le meilleur a posteriori, etc.

En pratique, les scores sont utilisés lorsqu'ils réunissent certaines propriétés.

-La DÉCOMPOSABILITÉi.e., la possibilité de calculer un score local en chaque noeud -garantit des temps de calcul moindres lors de modifications locales.

-La propriété d'ÉQUIVALENCE garantit que deux structures d'une même classe d'équivalence de Markov auront les mêmes scores. Cette propriété est particulièrement utile lorsque la recherche utilise l'espace des équivalents de Markov.

-La CONSISTANCE fait qu'un modèle contenant la distribution -contenue dans les données -obtiendra un meilleur score qu'un modèle ne la contenant pas. De plus cette propriété discrimine les modèles contenants tous deux la distribution mais dont le nombre de paramètres est différent. Il existe une variante locale de cette propriété : la consistance locale qui permet à GES d'améliorer localement les structures.

Méthodes

Certains des algorithmes de recherche utilisent des heuristiques de limitation de l'espace de recherche, comme le font les algorithmes K2 [START_REF] Cooper | A bayesian method for the induction of probabilistic networks from data[END_REF] ou MWST (Maximum Weight Spanning Tree) [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]. D'autres effectuent leur recherche directement dans l'espace complet à l'aide de métaheuristiques. L'algorithme glouton, ou Greedy Search [START_REF] Chickering | learning bayesian networks : search methods and experimental results[END_REF], progresse à l'aide d'opérations locales, mais converge dans des maxima locaux. Des techniques plus évoluées ont alors été utilisées :

iterated hill-climbing [START_REF] Campos | An iterated local search algorithm for learning Bayesian networks restarts based on conditional independence tests[END_REF], simulated annealing [START_REF] Janzura | A simulated annealing-based method for learning Bayesian networks from statistical data : Research Articles[END_REF], tabu search [START_REF] Munteanu | The EQ Framework for Learning Equivalence Classes of Bayesian Networks[END_REF], genetic algorithm [START_REF] Larrañaga | Structure Learning of Bayesian Networks by Genetic Algorithms : A Performance Analysis of Control Parameters[END_REF], etc.

En 2003, [START_REF] Blanco | Learning Bayesian networks in the space of structures by estimation of distribution algorithms[END_REF] puis [START_REF] Romero | Learning Bayesian networks in the space of orderings with estimation of distribution algorithms[END_REF] proposent d'utiliser les EDA (Estimation of Distribution Algorithms) [START_REF] Mülhenbein | From recombination of genes to the estimation of distributions I. Binary parameters[END_REF] pour palier les problèmes rencontrés avec les algorithmes génétiques et décrits dans [START_REF] Larrañaga | Combinatorial optimisation by learning and simulation of Bayesian networks[END_REF]. De nombreux travaux montrent l'efficacité des EDA dans de nombreux domaines [START_REF] Larrañaga | Estimation of Distribution Algorithms. A new tool for evolutionary computation[END_REF], y compris pour l'apprentissage de réseaux Bayésiens où les EDA surpassent les algorithmes génétiques [START_REF] Blanco | Learning Bayesian networks in the space of structures by estimation of distribution algorithms[END_REF][START_REF] Romero | Learning Bayesian networks in the space of orderings with estimation of distribution algorithms[END_REF][START_REF] Thibault | Learning Bayesian network structures by estimation of distribution algorithms : An experimental analysis[END_REF].

Algorithmes EDA

Les EDA sont des algorithmes « évolutionnaires » à base de populations. Ils reposent sur une connaissance explicite de la distribution des individus dans l'espace des solutions n'utilisant ainsi aucun opérateur de croisement ou mutation. Les EDA repartissent les individus d'une génération selon la distribution manipulée, puis réestiment cette distribution selon les meilleurs individus. Dans notre précédent travail expérimental [START_REF] Thibault | Learning Bayesian network structures by estimation of distribution algorithms : An experimental analysis[END_REF] dist ← distribution uniforme initiale 3:

répéter 4: population ← GÉNÈRE_POPULATION_SELON_DIST(dist, M ) 5: meilleurs ← SÉLECTIONNE_MEILLEURS(population, N ) 6: dist ← RÉESTIME_DISTRIBUTION(meilleurs) 7: jusqu'à CONVERGENCE(dist) == 1 8:
retourne meilleur_individu 9: fin fonction Illustrons l'algorithme sur One Max, un problème jouet dont le principe est de maximiser la somme des éléments d'un vecteur binaire. La fitness 2 choisie ici est le nombre de 1 par ligne.

ind. x1 x2 x3 fit x(1) 0 1 0 1 x(2) 0 1 0 1 x(3) 1 0 1 2 x(4) 1 0 1 2 x(5) 0 1 1 2 x(6) 1 0 0 1 dist 0.5 0.5 0.5 dist ′ 0.66 0.33 1 ⇒ ind. x1 x2 x3 fit x(1) 1 1 1 3 x(2) 1 1 1 3 x(3) 1 0 1 2 x(4) 1 0 1 2 x(5) 1 0 1 2 x(6) 0 0 1 1 dist 0.7 0.3 1 dist ′ 1 0.66 1 ⇒ . . .
Partant d'une distribution uniforme 6 individus sont générés (1.4). A partir d'une sélection (1.5) des 3 meilleurs (en gras), la distribution est ré-estimée (1.6) en calculant la probabilité fréquentielle de chaque partie de la solution (formule 6). Ce processus recommence jusqu'à un critère d'arrêt (1.7) -généralement la convergence.

dist ′ i = 1 N k∈meilleurs x i (k) [6]
2. La fonction à maximiser qui représente la capacité d'un individu à survivre et se reproduire.

Notre application de l'algorithme aux réseaux Bayésiens utilise une matrice qui associe à chaque arc X → Y une probabilité, ce qui suppose que la distribution des structures peut se calculer à partir d'une combinaison des probabilités de chaque arc. La règle de mise à jour est similaire à celle de l'équation 6. Notons qu'une telle structure nécessite l'utilisation d'un opérateur de coupure des cycles orientés.

Algorithme GES

En 2002, [START_REF] Chickering | Optimal structure identification with greedy search[END_REF] publient l'algorithme GES qui travaille directement dans l'espace des équivalents. Pour comprendre son fonctionnement, il faut revenir sur quelques notions des réseaux Bayésiens.

Une CARTE D'INDÉPENDANCE est un graphe G tel que (G, P ) respecte la condition de Markov [START_REF] Verma | Causal networks : Semantics and expressiveness[END_REF] (où P une distribution). On parle de carte d'indépendance minimale lorsqu'aucun graphe partiel de G n'est une carte d'indépendance de P. A l'inverse, une CARTE DE DÉPENDANCE est un graphe G ′ tel que toute indépendance dans P implique une séparation dans G ′ . Parallèlement, on peut définir la notion de carte maximale de dépendance. Une CARTE PARFAITE est l'unique carte qui soit à la fois une carte d'indépendance (minimale) et une carte de dépendance (maximale). Lorsque le nombre d'exemples est suffisamment grand, la distribution P de la base d'apprentissage est une carte parfaite du réseau original G (⇔ (G,P) est "faithful"). Notons que la séparation est une notion générale des graphoïdes [START_REF] Pearl | Probabilistic reasoning in intelligent systems : Networks of plausible inference[END_REF], sa version « orientée » est la d-séparation -qui s'applique donc aux réseaux Bayésiens.

Un arc couvert est un arc X → Y tel que les parents de X et Y sont identiques (à part X qui n'est pas un parent de lui-même). S'il existe une réprésentation DAGfaithful H de la distribution P, alors le théorème de Meek [START_REF] Meek | Graphical Models : Selecting causal and statistical models[END_REF] (théorème 1) prouve qu'il existe une suite finie d'opérations pour trouver la carte H. Cette conjecture a été démontrée par [START_REF] Chickering | Optimal structure identification with greedy search[END_REF]. Plusieurs algorithmes utilisent le résultat de Meek [START_REF] Meek | Graphical Models : Selecting causal and statistical models[END_REF] mais GES travaille dans l'espace des équivalents grâce à des opérateurs définit pour cet espace (Chickering, 2002a). Chaque classe d'équivalence est réprésentée par un CPDAG (Completed Partially Directed Acyclic Graph), ou pattern depuis que [START_REF] Verma | Equivalence and synthesis of causal models[END_REF] retourne G ⊲ carte minimale d'indépendance ⇔ carte parfaite 10: fin fonction 3. Comparaisons L'objet de cette étude est d'évaluer la qualité des résultats des approches par EDA en comparant avec l'algorithme GES (ou Greedy Equivalent Search) [START_REF] Chickering | Optimal structure identification with greedy search[END_REF], un algorithme très performant utilisant l'espace des équivalents de Markov. Pour cela, nous comparons les résultats obtenus sur plusieurs benchmarks à l'aide de plusieurs critères de qualité.

Benchmarks

Nous avons utilisés 5 réseaux tests classiques de la littérature : ASIA [START_REF] Lauritzen | Local computations with probabilities on graphical structures and their application on expert systems[END_REF], INSURANCE [START_REF] Binder | Adaptive Probabilistic Networks with Hidden Variables[END_REF], INSULIN [START_REF] Le | Using prior knowledge to improve genetic network reconstruction from microarray data[END_REF], ALARM [START_REF] Beinlich | The ALARM monitoring system : a case study with two probabilistic inference techniques for belief networks[END_REF], HAILFINDER [START_REF] Abramson | Hailfinder : A Bayesian system for forecasting severe weather[END_REF]. Pour chacun, une base d'apprentissage de taille arbitraire (selon le nombre de variables et de modalités) a été générée par la méthode "logic sampling" [START_REF] Henrion | Propagating uncertainty in Bayesian networks by probabilistic logic sampling[END_REF]. Chaque réseau a servi à la génération de 10 bases de d'apprentissage. Les résultats sont des moyennes sur 10 lancements différents avec les écarts-types correspondants.

Criteres de comparaison

Tandis que les EDA effectuent leur recherche dans B -l'espace des réseaux Bayésiens -GES utilise E -l'espace des équivalents. -l'a priori sur les structures p(B) = n i=1 0.001 (ri-1)qi , -l'a priori sur les paramètres

α ijk = 10 ri•qi , -donc, α ij = ri k=1 α ijk = 10
qi . Ce score est donné dans [START_REF] Chickering | Optimal structure identification with greedy search[END_REF] pour l'implémentation de GES. Dans nos expérimentations, il est implémenté avec un cache LRU (Least Recently Used) afin d'accélerer les calculs.

Afin de pouvoir comparer les structures, les DAG obtenus par l'EDA sont transformés en CPDAG par la méthode décrite par [START_REF] Chickering | A Transformational Characterization of Equivalent Bayesian Network Structures[END_REF]. Ainsi, il est possible de compter le nombre de différences structurelles par rapport au représentant de la classe d'équivalence du benchmark. Ces différences structurelles peuvent être de trois types : arc en trop, mauvaise orientation (inversion, absence ou présence abusive d'orientation), ou arc manquant. La distance d'édition (ED) est la somme de ces trois mesures, c'est à dire le nombre de changements élémentaires nécessaires pour retrouver la structure originale.

Enfin, le nombre d'itérations de chaque algorithme est présenté. A chaque itération de l'EDA tous les individus sont scorés, ce qui n'est pas le cas dans GES puisque seules les modifications locales applicables à l'individu sont scorées. A titre indicatif, nous donnons le temps total d'exécution, directement proportionnel au nombre d'appels effectifs à la fonctions de scorei.e., lorsque la valeur ne se trouve pas en cache.

Résultats

Le tableau 2 présente les écarts de score par rapport à la structure originale (référence). Afin de rendre compte de ces écarts, ils sont également donnés de façon relative. En terme de score, GES obtient de bien meilleurs résultats qu'EDA. Mieux, il obtient de meilleurs scores que celui du réseau original pour les 2 benchmarks IN- SULIN et HAILFINDER. Le En général, l'objectif des méthodes à base de score n'est pas de retrouver la structure initiale, mais de trouver un modèle à la fois simple et précis en inférence. Toutefois, la CONSISTANCE de notre score permet de distinguer les structures « proches » de la structure originale. Un score consistant augmente lorsque le nombre d'indépendances conditionnelles retrouvées augmente, mais diminue si trop d'indépendances supplémentaires sont encodées dans la structure. Nous pouvons donc observer les résultats en terme de défauts structurels (tableau 4) où f p représente le nombre de faux-positifs, f n les faux-négatifs et f o les fautes-d'orientation par rapport au réseau original 3 , la somme formant la distance d'édition (ed). Le

Conclusion

Dans ce travail, nous avons mis en évidence les limites des EDA en montrant que l'algorithme GES est plus performant en qualité et en temps d'exécution. Pour cela, nous avons comparé nos deux implémentations sur un ensemble de benchmarks connus, en utilisant les critères usuels. La démonstration de la conjecture de Meek offre un résultat théorique important, qui permet d'élaborer des algorithmes plus efficaces que les méthodes généralistes telles qu'EDA.

  Théorème 1 (rappel) Soit G et H deux DAG tels que H est une carte d'indépendance de G. Soit r le nombre d'arcs dans H qui ont une orientation opposée dans G, et m le nombre d'arcs de H qui n'existent dans aucune orientation dans G. Il existe une séquence d'au plus r + 2m opérations de réversion ou d'ajout d'arc dans G telle que : -chaque reversion concerne un arc couvert, -après chaque opération, G est un DAG et H une carte d'indépendance de G, -après toutes les réversions ou ajouts, G = H.

  ont montrés que tous les graphes d'une même classe possèdent le même squelette et les mêmes V-structures, i.e., les mêmes liaisons convergentes. L'opération qui oriente d'une façon quelconque les arcs non-orientés -appelée INSTANCIATION -est close dans l'espace des équivalents. GES procède en deux étapes : une étape gloutonne d'insertion d'arcs (dans E) qui permet d'augmenter le graphe jusqu'à l'obtention d'une carte d'indépendance, puis une étape de suppression d'arcs permettant de minimiser la carte d'indépendance, donc de trouver la carte parfaite. L'utilisation d'un score consistant est donc nécessaire pour découvrir les opérations d'améliorations. tant que arc ← TROUVER_AMELIORATION("deletion",G) faire 7:G ← RETIRE_ARC(G, arc) ⊲ minimisation de la carte d

  nous avons sélectionné la variante UMDA de l'algorithme, qui donnait de meilleurs résultats que la variante PBIL. Nous reprenons ici nos conclusions, en paramétrant l'algorithme avec 1000 individus, une sélection élitiste de 100 individus et sans conservation des individus d'une population sur la suivante.

	Algorithm 1 Algorithme EDA (variante UMDA)

1: fonction EDA(M, N ) ⊲ M taille de la population, N taille de la sélection 2:

  Réseaux tests : ASIA, INSURANCE, INSULIN, ALARM et HAILFINDER la propriété d'équivalence assure de pouvoir comparer les scores des structures obtenues par les algorithmes. Le score utilisé, DÉCOMPOSABLE, ÉQUIVALENT et CONSIS-TANT est BDeu avec :

	benchmarks	variables	type	arcs données
	ASIA	8	binaires	8	2 000
	INSURANCE	27	discrètes	52	10 000
	INSULIN	35	binaires	52	5 000
	ALARM	37	discrètes	46	10 000
	HAILFINDER	56	discrètes	66	30 000
	Tableau 1.
			L'utilisation d'un score possédant

  score utilisé étant CONSISTANT, cela signifie que le réseau Ecart absolu des scores (ecart -score_ref ) et écart relatif (| ecart score_ref |)obtenu par l'algorithme représente mieux la distribution que le benchmark lui-même. Remarquons que l'écart-type est généralement plus faible pour GES, ce qui montre l'aspect très stochastique des EDA.Le nombre d'itérations rend compte de la faculté de l'algorithme à converger, i.e., trouver une solution en un temps « acceptable ». Chaque génération d'EDA correspond ici à l'évaluation complète de 1000 individus (hors utilisation du cache). La méthode est donc plus gourmande que GES qui n'évalue en général qu'un nombre assez réduit d'opérateurs. On remarquera encore l'écart-type plus important dans les EDA (tableau 3). La grande différence de temps observée est une conséquence du mode Tableau 3. Nombre moyen d'itérations et temps moyen d'exécution avant convergence de fonctionnement de chaque algorithme. Quand l'un procède par « petites » modifications locales, améliorant de fait considérablement la solution, l'autre dissémine de nombreux individus dans l'espace de recherche afin de le « fouiller ».

					écart absolu	relatif (%)
	scores	benchmark		EDA	GES	EDA GES
	ASIA	-4741.8	-13.66 ± 8.45	-1.86 ± 9.66 0.29	0.04
	INSURANCE	-149113	-5128 ± 1701	-3829 ± 963	3.44	2.57
	INSULIN	-87253.6	-970.4 ± 376	+152.4 ± 127	1.11 -0.17
	ALARM	-110637	-7311 ± 1049	-2925 ± 120	6.61	2.64
	HAILFINDER	-1500330 -67328 ± 10434	+4364 ± 92	4.49 -0.29
	Tableau 2. nombre d'itérations	temps d'exécution (en s.)
	temps		EDA		GES	EDA	GES
	ASIA		11.5 ± 2.55	7.0 ± 0.47	0.4 ± 0.1	0.0 ± 0.1
	INSURANCE	26.9 ± 1.91 42.1 ± 1.91	103.4 ± 1.4	1.2 ± 0.1
	INSULIN		28.2 ± 2.86 47.1 ± 1.20	77.3 ± 1.2	1.8 ± 0.1
	ALARM		30.8 ± 1.75 42.9 ± 1.52	168.0 ± 2.3	2.3 ± 0.3
	HAILFINDER	34.8 ± 4.87 65.5 ± 0.71 1794.7 ± 16.7 27.8 ± 0.5

  nombre d'erreurs d'orien-± 3.2 20.0 ± 1.9 23.6 ± 3.0 54.9 ± 6.1 EDA INSULIN 8.4 ± 4.0 13.5 ± 2.8 10.0 ± 2.1 31.9 ± 7.0 ALARM 21.4 ± 2.8 11.5 ± 2.0 28.1 ± 2.3 61.0 ± 4.2 HAILFFINDER 26.9 ± 4.3 33.4 ± 2.8 28.3 ± 3.1 88.6 ± 6.0 ± 2.5 20.3 ± 2.1 10.5 ± 2.0 46.6 ± 4.9 Tableau 4. Qualité de la structure obtenue avec les algorithmes EDA et GES tation est légèrement biaisé : calculant ce critère sur les des CPDAG, une mauvaise orientation d'un arc entrainera souvent une mauvaise orientation de ses voisins par propagation. Malgrès cet obstacle, nous avons conservé cette méthode de comparaison, car chaque différence représente un changement de classe d'équivalence.

	erreurs	fp	fn	fo	ed
	ASIA	1.5 ± 0.7	1.1 ± 0.6	4.6 ± 1.7	7.2 ± 2.1
	INSURANCE 11.3 ASIA 0.2 ± 0.6	1.2 ± 0.8	0.3 ± 0.7	1.7 ± 1.8
	INSURANCE	6.1 ± 1.7 16.0 ± 1.3 14.8 ± 0.8 36.9 ± 2.8
	GES INSULIN	1.0 ± 0.5	5.9 ± 1.1	7.5 ± 1.8 14.4 ± 2.8
	ALARM	5.1 ± 0.3	8.4 ± 0.7	4.4 ± 0.5 17.9 ± 0.3

HAILFFINDER

15.8