Control of Drilling Vibrations:

A Time-Delay System-Based Approach Islam Boussaada * , * * Arben Cela * * * Hugues Mounier * Silviu-Iulian Niculescu * * Laboratoire des Signaux et Systèmes (L2S) Supélec-CNRS-Université Paris Sud, 3 rue Joliot-Curie 91192 Gif-sur-Yvette cedex (France) (e-mail: Silviu.Niculescu,Hugues.Mounier,islam.boussaada@lss.supelec.fr) * * Laboratoire de Modélisation et Calcul Scientifique (LMCS) Institut Polytechnique des Sciences Avancées, 7 rue Maurice Grandcoing 94200

Ivry-sur-Seine (France) * * * Department of Embedded Systems, UPE, 93162 Noisy-Le-Grand (France) (e-mail: celaa@esiee.fr)

Abstract: * The main purpose of this study is the control of both axial and torsional vibrations occurring along a rotary oilwell drilling system. This work completes a previous author's paper (Boussaada et al. [2012a]) which presents the description of the qualitative dynamical response of a rotary drilling system with a drag bit, using a model that takes into consideration the axial and the torsional vibration modes of the bit. The studied model, based on the interface bit-rock, contains a couple of wave equations with boundary conditions consisting of the angular speed and the axial speed at the top additionally to the angular and axial acceleration at the bit whose contain a realistic frictional torque. Our analysis is based on the center manifold theorem and normal forms theory whose allow us to simplify the model. By this way we design two control laws allowing to suppress the undesired vibrations guaranteeing a regular drilling process.

INTRODUCTION

Interconnected oscillatory systems often display what is called propagation phenomena, [START_REF] Hale | Introduction to functional differential equations[END_REF].

In general by Lossless propagation it is understood the phenomenon associated with long transmission lines for physical signals. In engineering, this problem is strongly related to electric and electronic applications, e.g. circuit structures consisting of multipoles connected through LC transmission lines; this can also be seen in steam or gaz flows or pressures and water pipes, [START_REF] Niculescu | Delay Effects on Stability, A Robust Control Approach[END_REF], [START_REF] Fu | On the Stability of Linear Delay-Differential Algebraic Systems: Exact Conditions via Matrix Pencil Solutions[END_REF], [START_REF] Rasvan | Oscillations in lossless propagation models: a liapunov-krasovskii approach[END_REF]. The mathematical model is described in all these cases by a mixed initial and boundary value problem for hyperbolic partial differential equations modeling the lossless propagation.

The boundary conditions are of special type, being in feedback connection with some system described by ordinary differential equations. This leads to the so-called derivative boundary conditions considered in [START_REF] Cooke | Differential-difference equations and nonlinear initial boundary value problems for linear hyperbolic partial differential equations[END_REF], but also to the even more general boundary conditions of Abolina & Myshkis described by Volterra operators, see [START_REF] Rasvan | Oscillations in lossless propagation models: a liapunov-krasovskii approach[END_REF]. Integration along characteristics of the hyperbolic partial differential equations (here d'Alembert method) allows the association of certain system of functional equations to the mixed problem.

This paper is concerned by an application which can be modeled by such equations; therefore, the above idea is adopted, see [START_REF] Balanov | Bifurcation analysis of a neutral delay differential equation modelling the torsional motion of a driven drill-string[END_REF], [START_REF] Fridman | Bounds on the response of a drilling pipe model[END_REF], [START_REF] Rouchon | Flatness and stick-slip stabilization[END_REF], [START_REF] Saldivar | Stick-slip oscillations in oillwell drilstrings: Distributed parameter and neutral type retarded model approaches[END_REF]. The analysis and modeling of rotary drilling vibrations is a topic whose economical interest has been renewed by recent oilfields discoveries leading to a growing literature, see for instance [START_REF] Th | A simplified model to explore the root cause of stick-slip vibrations in drilling system with drag bits[END_REF], [START_REF] Germay | Nonlinear drilling dynamics analysis[END_REF], Navarro-López [2009], [START_REF] Navarro-López | Sliding-mode control of a multi-dof oilwell drillstring with stick-slip oscillations[END_REF], Navarro-López and Suárez [2004a], Boussaada et al. [2012a] and Navarro-López and Suárez [2004b].

Previous work by the authors Boussaada et al. [2012a] improved the modelling of vibrations of the drilling system, taking into account both axial and torsional vibrations and secondly by extending the qualitative analysis to the investigation of the nonlinear terms in the model. The use of the center manifold theorem and normal forms theory allow to the analysis of a finite dimensional approximation which conserves the main dynamics of the physical original system. Let us consider the following model for the axial vibrations U and torsional vibrations Φ:

   ∂ 2 t U (t, s) = c 2 ∂ 2 s U (t, s) E Γ ∂ s U (t, 0) = α∂ t U (t, 0) -H(t) M ∂ 2 t U (t, L) = -E Γ∂ s U (t, L) + F (∂ t U (t, L)) , (1)
and

     ∂ 2 t Φ(t, s) = c2 ∂ 2 s Φ(t, s) G Σ ∂ s Φ(t, 0) = β∂ t Φ(t, 0) -Ω(t) J ∂ 2 t Φ(t, L) = -G Σ ∂ s Φ(t, L) + F (∂ t U (t, L)) , (2)
where, in equation ( 1), H is the brake motor control and α∂ t U (t, 0) represents a friction force of viscous type. For equation ( 2), the right hand side of the second equation designates the difference between the motor speed and rotational speed of the first pipe. The physical parameters of the model ( 1)-( 2) are: G is the shear modulus of the drillstring steel and E the elasticity Young's modulus.

Then the wave speeds can be expressed by c = E/ρ and c = G/ρ and J the inertia J = M r 2 where r is taken as the averaged radius of drillpipe and Γ is the averaged section of the drillpipe and Σ is the quadratic momentum.

It is worth mentioning that the coupling term appears in the boundary condition of the torsional vibrations. Those parameters are taken following the numerical settings presented in the Appendix. The nonlinear aspect of the model is considered by taking functions F and F in the form: z → pk z/(k 2 z 2 + ζ) where the parameters p, k, ζ are some positive integer responsible of the sharpness of the top angle of the friction force graph and p is some parameter deciding the amplitude of the friction force such that 0 < ζ ≪ 1 and 0 < k < 1. Moreover, the behavior of the chosen friction model is close from the empirical model: the white friction force but is more handleable, which can be very useful in experimental identifications.

Note also that the proposed model can be expanded to Taylor sum, which is very important when the aim is to give accurate approximation at any fixed order. The chosen functions have a close behavior to the one used in [START_REF] Barton | Nonlinear dynamics of torsional waves in a drill string model with spacial extent[END_REF] for modelling the friction. In a previous work by the authors, Boussaada et al. [2012a] an analytical study of the uncontrolled drilling vibrations model as functional differential equations of neutral type is established, this is done by a qualitative theory method; center manifold theorem [START_REF] Carr | Application of Center Manifold Theory[END_REF] and normal forms theory [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcation of vector fields[END_REF]. Indeed, most of the references concerned by partial differential equations (PDE) or delay differential equations (DDE) models for the drilling problem have a numerical analysis character.

This work completes the contribution of Boussaada et al. [2012a] by including some design approaches. Thus we adopt the same qualitative methodology; we reduce the considered PDE model to a singularly perturbed system of ordinary differential equations (ODE), in the goal of designing appropriate control laws allowing to suppress these vibrations, which guarantee the desired drilling process (helicoidal evolution of the bit). First we propose a PID controller, the considered spectral projection allows to reduce the PDE system to a singularly perturbed system of ODE. Secondly, motivated by the technological constraints due to the use of wireless technology which induces delay, we propose a more realistic control law which consists of delayed feedback which might be more appropriate.

In this case a Bogdanov-Takens singularity (double zero eigenvalues) is considered. Furthermore, to the best of the authors knowledge, this type of singularity has never been studied for NDDE depending on parameters; thus we extend the methodology for computing the center manifold. Similar results can be found in [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF] with the analysis of a physiological control model of DDE with double-zero eigenvalue singularity. This study has the same spirit as the results of Bogdanov and Takens for ODEs. We refer the reader to [START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF], Gucken-heimer and Holmes [2002] for elements on the qualitative theory of differential equations and bifurcation theory.

The remaining part of the paper is organized as follows.

The second section is concerned by preliminaries, we describe the standard procedure for reducing the PDE drillstring model to a neutral delay differential equations (NDDE). In the third section, entitled Control of drilling vibrations, we present two controllers allowing to a stable drilling process. The methodological scheme described in Ait [START_REF] Babram | Computational scheme of a center manifold for neutral functional differential equations[END_REF], [START_REF] Campbell | Calculating centre manifolds for delay differential equations using maple[END_REF] is extended to the study of the parametrized model of neutral type. For the sake of self-containment, we report in the Appendix a table for the numerical settings for the parameters used in (1)-( 2). We refer the reader to Boussaada et al. [2012a] for the outlines of the methodology enabling to approximate a system of NDDE by a system of ordinary differential equations (center variety) and then the study of local bifurcations (normal forms) for for further insights in FDE see [START_REF] Campbell | Calculating centre manifolds for delay differential equations using maple[END_REF], [START_REF] Weedermann | Hopf Bifurcation calculations for scalar delay differential equations[END_REF], [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF].

PRELIMINARIES AND PREREQUISITES

To the best of the authors knowledge, the standard procedure allowing to transform the considered PDE model to a delay system of neutral type was presented for the first time in [START_REF] Cooke | Differential-difference equations and nonlinear initial boundary value problems for linear hyperbolic partial differential equations[END_REF], see also [START_REF] Balanov | Bifurcation analysis of a neutral delay differential equation modelling the torsional motion of a driven drill-string[END_REF] and [START_REF] Mounier | Propriétés structurelles des systèmes linéaires a retard : Aspects théorique et pratique[END_REF]. Indeed, by using d'Alembert theorem, in Boussaada et al. [2012a] the system of PDE (1)-( 2) is reduced to a system of NDDE. We adopt the normalization such that the units of length, time and torque the quantities L, T = L/c and EΓ/L, thus system (1)-( 2) is written

                                                 v(t) - α -1 α + 1 v(t -2) = - 1 M v(t) - α -1 M (α + 1) v(t -2) + 2 M (α + 1) H(t -1) + 1 M F ( v(t)) - α -1 M (α + 1) F ( v(t -2)) ẅ(t) - c E Γ β -c G Σ c E Γ β + c G Σ ẅ(t -2 τ ) = - c G Σ c E Γ J ẇ(t) - c G Σ c E Γ J c E Γ β -c G Σ c E Γ β + c G Σ ẇ(t -2τ ) + 1 J F ( v(t)) - c E Γ β -c G Σ J (c E Γ β + c G Σ) F ( v(t -2τ )) + 2c G Σ J (c E Γ β + c G Σ) Ω(t -τ )
(3) where τ is the ratio of the speeds τ = c c , and v and w are respectively the axial flat output v(t) = U (t, L) and the torsional flat output w(t) = Φ(t, L).

CONTROL OF DRILLING VIBRATIONS

Delayed feedback controller

Let denote by x 1 = v the flat output associated with the axial vibrations, x 2 = w the flat output associated with the torsional vibration, x 3 = v and by x 4 = ẇ and consider the matrix representation of the linear part of the above system where x = (x 1 , x 2 , x 3 , x 4 ) T and let set

H(t) = p H v(t -1) and Ω(t) = p Ω w(t -τ ).      ẋ(t) =D 1 ẋ(t -2) + D 2 ẋ(t - 2c c ) + A 0 x(t) + A 1 x(t -2) + A 2 x(t - 2c c ) + F(x(t), x(t -2), x(t - 2c c
))

(4) where F is the nonlinear part of the system (3)

D 1 =       0 0 0 0 0 0 0 0 0 0 d 1,1,1 0 0 0 0 0       , D 2 =       0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 d 2,2,2       , A 0 =       0 0 1 0 0 0 0 1 0 0 a 0,1,1 0 0 0 a 0,2,1 a 0,2,2       , A 1 =      0 0 0 0 0 0 0 0 p H 0 a 1,1,1 0 0 0 0 0      , A 2 =       0 0 0 0 0 0 0 0 0 0 0 0 0 p Ω a 2,2,1 a 2,2,2      
where the matrices coefficients are

d 1,1,1 = α-1 α+1 , d 2,2,2 = c1E(Γ) β-c2GΣ c1E(Γ) β+c2GΣ , a 0,1,1 = pk-ζ M ζ , a 0,2,1 = a2k Jζ , a 0,2,2 = -c2GΣ c1E(Γ) J , a 1,1,1 = -(α-1)(ζ+pk) (α+1)M ζ , a 2,2,1 = -pk(c1E(Γ) β-c2GΣ) ζ (c1E(Γ) β+c2GΣ)J , a 2,2,2 = -c2GΣ (c1E(Γ) β-c2GΣ) Jc1E(Γ) (c1E(Γ) β+c2GΣ
) . Recall that in the above quoted references (concerned by PDE models), the studies were concerned only by the torsional vibrations. Thus the associated NDDE (governing the speed of such vibrations) is scalar, which is easier to study compared with (4). And for the physiological model considered in [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF], A 2 = D i = 0 for i ∈ {1, 2} since the model is DDE with one delay. It is worth mentioning that delay of PD controller is chosen to be the normalized proper delay of the system. Setting the numerical values of the physical parameters given in the Appendix we have the following result Proposition 1. When p Ω = p H = 0:

• Zero is the only eigenvalue with zero real part and the remaining eigenvalues are with negative real parts. Moreover, zero is an eigenvalue of algebraic multiplicity 2 and of geometric multiplicity 1, that is zero eigenvalue is non-semisimple and the singularity is of Bogdanov-Takens, see [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcation of vector fields[END_REF]. • The system (4) is formally stable but not asymptotically stable (although there are no characteristic roots with positive real parts).

When p H = -24δ r and p Ω = 45µr 2 10 for a small parameter r:

• The dynamics of (4) reduces on a cubic center manifold to ż1 = z 2 ż2 = δz 1 + µz 2 -3z 2 z 2 1

(5)

for which the function

I(z) = 1 2 z 1 2 - 1 2 z 2 2 δ - 1 2 δ z 1 2 z 2 2 + 1 4 δ 2 z 1 4 + 1 4 z 2 4
is a Lyapunov function when δ < 0 and µ < 0 and then the system is globally asymptotically stable.

Sketch of the proof:

The first assertion is obtained by establishing a linear analysis of the system. Indeed, it can be easily checked by computing the associated characteristic equation and substituting the physical values. Numerical tools as QPMR Vyhlídal & Zítek [2003] are also very useful for locating those spectral values. Let us stress on the second part of the proposition which concerns the nonlinear analysis. We follow the approach described in [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF] that considers a singular delay system linearly dependent on a parameter, and in the same spirit to the decomposition established in [START_REF] Faria | Normal forms for retarded functional-differential equations with parameters and applications to Hopf bifurcation[END_REF] in the goal of computing the normal form for delay systems depending on a parameter, we extend the scheme of computing the center manifold to the case of NDDE depending on parameters and thus look for the system (4) as a perturbation of

d dt D x t = L 0 x t , where L 0 = L| {p H =0, pΩ=0} (6)
Indeed, system (4) can be written as

d dt D x t : = L 0 x t + F(x t ) = L 0 x t + (L -L 0 ) x t + F(x t ) (7)
such that

F µ,p = -0.0405 x 3 1 (t) + 0.0377 x 3 3 (t -2) + p H x 1 (t -2) -1.875 px 3 1 (t) + 1.874998 px 3 1 (t -1.264911064)
Here we follow the theoretical schemes briefly presented in [START_REF] Campbell | Calculating centre manifolds for delay differential equations using maple[END_REF], [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF] and give computations steps for the equation of the evolution of the problem's solutions on the center variety for system (6). First, we compute the basis of the generalized eigenspace corresponding to the double eigenvalue λ 0 = 0.

Φ(θ) T = 1 + θ 1 + 2 θ 1 2 1 2 0 0 , where θ ∈ [-2, 0]. Recall that the adjoint linear equation associated to (4) is      u(t) =D 1 u(t + 2) + D 2 u(t + 2c c ) -A 0 u(t) -A 1 u(t + 2) -A 2 u(t + 2c c ) (8)
with a basis for the generalized eigenspace associated to the double eigenvalue zero is given by

Ψ(θ) = -1 -3 -7 -13
ξ + 1 3 ξ + 2 7 ξ + 3 13 ξ + 4 .

Let us consider the bilinear form, see [START_REF] Hale | Introduction to functional differential equations[END_REF] (ψ, ϕ) =ψ(0

)(ϕ(0) -D 1 ϕ(-2) -D 2 ϕ(-1.264911)) + 0 -2 ψ(ξ + 2)A 1 ϕ(ξ)dξ + 0 -1.264911 ψ(ξ + 1.264911)A 2 ϕ(ξ)dξ - 0 -2 ψ ′ (ξ + 2)D 1 ϕ(ξ)dξ - 0 -1.264911 ψ ′ (ξ + 1.264911)D 2 ϕ(ξ)dξ. (9)
By using (9) we can easily normalize Ψ such that (Ψ, Φ) = I d , thus the space C can be decomposed as

C = P Q, where P = {ϕ = Φz; z ∈ R 2 } and Q = {ϕ ∈ C; (Ψ, ϕ) = 0}.
Recall that each of those subspaces is invariant under the semigroup T (t) and that the matrix B (introduced in the previous section concerned by the theoretical settings) satisfying AΦ = ΦB is given by

B = 0 0 1 0 . ( 10
)
Let us first set the following decomposition x t = Φy(t) + z(t) where z(t) ∈ Q and y(t) ∈ R 2 , z(t) = h(y(t)) and h is some analytic function h : P → Q. Thus the explicit solution on the center manifold can be obtained by the use of the proven formula in [START_REF] Campbell | Calculating centre manifolds for delay differential equations using maple[END_REF], [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF] that is

ẏ(t) = By(t) + Ψ(0)F[Φ(θ)y(t) + h(θ, y(t))] (11)
∂h dy {By + Ψ(0)F[Φ(θ)y + h]} + Φ(θ)Ψ(0)F[Φ(θ)y + h] =    ∂h dθ , -2 ≤ θ ≤ 0 L(h(θ, y)) + F[Φ(θ)y + h(θ, y)], θ = 0 ( 12
) where h = h(θ, y) and F is defined in (7).

Since our aim is to establish the parameter values of p H and p Ω guaranteeing an asymptotic suppression of the vibrations after some fixed time t 0 we study the parameter bifurcations, the computation of the evolution equation of of the problem's solutions on the center variety for system (7) is required. In the next step, we introduce a small parameter r as a scaling parameter for making a zoom into the neighborhood of the singularity. We introduce the following changes of coordinates p H = -24δ r, p Ω = 45µr 2 10 , y 1 = rz 1 , y 2 = r 2 z 2 and we scale the time by t old = r t new which allows us to the following cubic normal form reduction of (4), ż1 = z 2 ż2 = δz 1 + µz 2 -z 3 1 for which a normal form is given by (5). Under the negativity condition of δ and µ positive function beyond (0, 0) and if δ < 0 and µ > 0 then İ(z 1 , z

2 ) = -z2 2 µ δ + 3 z2 2 z1 2 δ - δ z 1 2 z 2 2 µ + 3 δ z 1 4 z 2 2 + z 2 4 µ -3 z 2 4 z 1 2
which is always negative. Thus the system is globally asymptotically stable and the undesired vibrations are suppressed. The proposed scheme offset the computation of a Lyapunov function for a system of PDE with nonlinear boundary conditions.

PID controller

Let denote by x 3 = v the flat output associated with the axial vibrations, x 4 = w the flat output associated with the torsional vibration, x 5 = v, by x 6 = ẇ and x 1 = t 0 x 3 (s) ds and x 2 = t 0 x 4 (s) ds and let us consider the matrix representation of the linear part of the above system where x = (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ) T and let set H(t) = H p v(t -1) + H d v(t -1) and Ω(t) = Ω p w(tτ ) + Ω i t-τ 0 w(s)ds, then the system

     ẋ(t) =D 1 ẋ(t -2) + D 2 ẋ(t - 2c c ) + A 0 x(t) + A 1 x(t -2) + A 2 x(t - 2c c ) + F(x(t), x(t -2), x(t - 2c c
))

(13) where F is the nonlinear part of the system (3) with matrices

D 1 =        
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 d 1,1,1 0 0 0 0 0 0 0         , D 2 =      
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 d 2,2,2       , A 1 =        
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a 1,1,1 0 0 0 0 0 0 0

        , A 2 =            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a 2,2,1 a 2,2,2            A 0 =            0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 H p 0 H d + a 0,1,1 0 0 Ω i 0 Ω p a 0,2,1 a 0,2,2           
where the matrices coefficients are

d 1,1,1 = α-1 α+1 , d 2,2,2 = c1E(Γ) β-c2GΣ c1E(Γ) β+c2GΣ , a 0,1,1 = pk-ζ M ζ , a 0,2,1 = a2k Jζ , a 0,2,2 = -c2GΣ c1E(Γ) J , a 1,1,1 = -(α-1)(ζ+pk) (α+1)M ζ , a 2,2,1 = -pk(c1E(Γ) β-c2GΣ) ζ (c1E(Γ) β+c2GΣ)J , a 2,2,2 = -c2GΣ (c1E(Γ) β-c2GΣ) Jc1E(Γ) (c1E(Γ) β+c2GΣ
) . Setting the numerical values of the physical parameters given in the Appendix we have the following result Proposition 2. When Ω p = Ω i = H p = H d = 0:

• Zero is the only eigenvalue with zero real part and the remaining eigenvalues are with negative real parts. Moreover, zero is an eigenvalue of algebraic multiplicity 4 and of geometric multiplicity 1, that is a generalized Bogdanov-Takens singularity, see [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcation of vector fields[END_REF]. • The system (13) is formally stable but not asymptotically stable (although there are no characteristic roots with positive real parts).

When H d = 13.27 r 11 δ 3 , H p = 26.55 r 12 δ 2 , Ω i = -11.47 r 13 δ 1 , Ω p = 11.47 r 10 δ 4 for a small parameter r:

• The dynamics of (13) reduces on a degree 15 center manifold to ż1

= z 2 , ż2 = z 3 , ż3 = z 4 ż4 = δ 1 z 1 + δ 2 z 2 + δ 3 z 3 + δ 4 z 4 -z 3 4 (14)
for which the the linear part is written as a Companion matrix,

A =       0 1 0 0 0 0 1 0 0 0 0 1 δ 1 δ 2 δ 3 δ 4      
thus there exist values for δ 1 , δ 2 , δ 3 , and δ 4 such that A is Hurwitz, which guaranteeing local asymptotic stability.

Sketch of the proof: By the same way as for the previous control law, we use the spectral projection methodology. The singularity here is zero with algebraic multiplicity 4 and geometric multiplicity 1. A basis for the generalized eigenspace M 0 here is

Φ(θ) =               1 6 θ 3 + θ 2 + θ -1 1 2 θ 2 + 2 θ + 1 θ + 2 1 - 1 6 θ 3 + 1 2 θ 2 -θ + 2 - 1 2 θ 2 + θ -1 -θ + 1 -1 θ 2 + θ 2 θ + 1 2 0 θ + 1 1 0 0 4 θ 4 0 0 1 0 0 0              
The matrix B associated to this type of singularity and staisfying AΦ = BΦ is given by Using the following changes of coordinates (r is a sufficiently small parameter) H d = 13.27 r 11 δ 3 , H p = 26.55 r 12 δ 2 , Ω i = -11.47 r 13 δ 1 , Ω p = 11.47 r 10 δ 4 , y = r 5 z 1 , y 2 = r 4 z 2 , y 3 = r 3 z 3 , y 4 = r 2 z 4 and an appropriate scaling of time allows us to the reduced system (14).

B =       0
For the complete proofs see for instance Boussaada et al. [2012b] the full version of the present paper.

CONCLUDING REMARKS

The main purpose of this paper is the design of control laws allowing to suppress the undesired vibrations using a FDE qualitative analysis based approach. We establish the appropriate gains of a PID controller as well as the gains of a delayed feedback controller. Moreover, in the case of delayed feedback controller we arrive to establish a Lyapunov function guaranteeing global stability for the central dynamics, which guarantees local asymptotic stability for the initial PDE and permits to achieve the desired dynamic for the drilling process. From a realistic point of view we suggest the delayed feedback controllers since the data collection comes delayed due to the use of wireless technology which induces a delay τ = 2.2 s.

Ψ

  basis of the adjoint space satisfying (Ψ, Φ) = I d , where evaluated at zero gives

ACKNOWLEDGEMENTS

The authors would like to thank the reviewers for valuable comments.

APPENDIX

Graphical Illustration

The projection of the dynamics on the center manifold for the critical value δ = µ = 0 is given in Figure 1 and


Quentin Maillet
email: quentin.maillet@mines-paristech.fr

Huan Xu

Necmiye Ozay

Richard M Murray

Dynamic State Estimation in Distributed Aircraft Electric Control Systems via Adaptive Submodularity

We consider the problem of estimating the discrete state of an aircraft electric system under a distributed control architecture through active sensing. The main idea is to use a set of controllable switches to reconfigure the system in order to gather more information about the unknown state. By adaptively making a sequence of reconfiguration decisions with uncertain outcome, then correlating measurements and prior information to make the next decision, we aim to reduce the uncertainty. A greedy strategy is developed that maximizes the one-step expected uncertainty reduction. By exploiting recent results on adaptive submodularity, we give theoretical guarantees on the worst-case performance of the greedy strategy. We apply the proposed method in a fault detection scenario where the discrete state captures possible faults in various circuit components. In addition, simple abstraction rules are proposed to alleviate state space explosion and to scale up the strategy. Finally, the efficiency of the proposed method is demonstrated empirically on different circuits.

I. INTRODUCTION AND MOTIVATION

The increasing focus on more-electric aircraft within the aerospace industry signifies progress in the direction of more energy efficient vehicles. Electric systems are replacing pneumatic, mechanical, and hydraulic subsystems, thus reducing weight, easing maintenance, and improving aircraft controllability and configurability [START_REF] Moir | Aircraft Systems: Mechanical, Electrical and Avionics Subsystems Integration[END_REF]. Yet as more subsystems rely on electric power, the flight-criticality of an aircraft becomes more dependent on the electric power system as well. Because the state of the system is determined solely from sensor measurements, the problem of state estimation from sensor readings is crucial to the safety of the entire aircraft.

Estimation of electric power systems using optimizationbased techniques is a well-established area [START_REF] Abur | Power system state estimation: theory and implementation[END_REF], [START_REF] Bose | Real-time modeling of power networks[END_REF], [START_REF] Monticelli | Electric power system state estimation[END_REF]. A large body of work exists on diagnostics of electric power systems focusing on AC systems [START_REF] Davidson | A toolset for applying model-based reasoning techniques to diagnostics for power systems protection[END_REF], as well as large vehicle systems. [START_REF] Liu | Fault detection, diagnostics, and prognostics: software agent solutions[END_REF] examines the diagnostics for the international space station, [START_REF] Keller | Aircraft electrical power systems prognostics and health management[END_REF] for an aircraft electric system, and [START_REF] Gonzalez | Model-based, real-time control of electrical power systems[END_REF] for a marine vehicle power system. For a DC system, [START_REF] Gorinevsky | Estimation of faults in dc electrical power system[END_REF] uses an optimization-based approach to estimate fault states.

Previous work in electric power system state estimation has focused on static, centralized estimation problems with continuous states. We perform discrete state estimation using active control of switches within a distributed control architecture. The system reconfigures itself through a set of controllable contactors (i.e., electrically controlled switches).

Once reconfigured, new sensor measurements are taken to gain more information about the unknown state. We adaptively sequence switching actions by use of a greedy strategy that maximizes the one-step expected uncertainty reduction. By exploiting recent results in adaptive submodularity [START_REF] Krause | Submodular function maximization[END_REF], [START_REF] Golovin | Adaptive submodularity: Theory and applications in active learning and stochastic optimization[END_REF], we provide theoretical bounds for the worst-case performance of the greedy strategy. Such dynamic state estimation techniques have been proposed in the context of Markov jump linear systems [START_REF] Blackmore | Active estimation for jump markov linear systems[END_REF], information gathering in robotics [START_REF] Mihaylova | A comparison of decision making criteria and optimization methods for active robotic sensing[END_REF], [START_REF] Singh | Efficient planning of informative paths for multiple robots[END_REF], active hypothesis testing [START_REF] Naghshvar | Active sequential hypothesis testing[END_REF], and active learning [START_REF] Guillory | Average-case active learning with costs[END_REF]. To the best of our knowledge, these ideas have not been applied before in electric power system state estimation and fault diagnosis problems.

Recently, correct-by-construction control synthesis has been applied to power allocation and distribution in aircraft electric power systems [START_REF] Wongpiromsarn | Formal synthesis of embedded control software for vehicle management systems[END_REF], [START_REF] Ozay | Distributed power allocation for vehicle management systems[END_REF], [START_REF] Xu | A case study on reactive protocols for aircraft electric power distribution[END_REF]. A critical assumption in these papers is that the high-level reactive control protocol has an accurate knowledge of the system states, including fault states, so that it can reroute the power accordingly. An expensive, hence undesirable, solution to achieve accurate state estimates is to equip the system with a large number of sensors. Software, however, is cheaper and more amenable to change than hardware.

The goal of the current paper is to obtain high-accuracy state estimates with a limited number of sensors by utilizing software-based dynamic estimation strategies. We are particularly interested in detecting and localizing faults in the system. As it is common to use discrete models for fault diagnosis [START_REF] Sampath | Failure diagnosis using discrete-event models[END_REF], continuous values of voltage and current, as well as health statuses of components in the system, are discretized before performing state estimation. Additionally, a discrete framework is well-suited for combining the proposed estimation strategy with control synthesis results described in [START_REF] Xu | A case study on reactive protocols for aircraft electric power distribution[END_REF] in future work.

II. PROBLEM SETUP

A. General problem description

Consider an aircraft electric power system topology, which can be represented by a graph data structure G = (N , E). Figure 1 shows a representative single-line diagram. The set N of nodes in the graph contains the following components: generators (G), rectifier units (R), and voltage sensors (S).

The set E of edges contains all contactors (and solid wire links) between components. The status of contactors C ⊆ E can either be open or closed. A node corresponding to a rectifier unit has no outgoing edges on the AC side and no incoming edges on the DC side to reflect the fact that they contain a diode. The rest of the edges in the graph are bidirectional.

Elements in the set of generators G ⊆ N and rectifier units R ⊆ N are uncontrollable, and can take values of unhealthy (i.e., the component is online but outputting a voltage not in admissible range), healthy (i.e., the component is online and outputting the correct voltage), or offline (i.e., no power output, open circuit). Measurements read from the sensors S ⊆ N will depend on the status of generators, rectifier units, and contactors. We say that there is a live path between two components if there exists a simple path in the graph G that connects the two nodes corresponding to these components, there is no offline component along the path including end nodes, and the contactors along this path are all closed. The readings of a sensor s ∈ S can then take the following values (i) improper voltage: if there is a live path between s and some g ∈ G (not offline by definition of live path), and either g or some rectifier r ∈ R along such a path is unhealthy; (ii) admissible voltage: for all g ∈ G that have a live path to s, both g and r along such paths are healthy; or (iii) no voltage: there is no live path between s and any generator g ∈ G.

On top of the circuit and sensing topology is a distributed control architecture with a dynamic state estimation mechanism. We assume that one of the embedded controllers is responsible for dynamic state estimation, hereafter referred to as the fault detection controller. The fault detection controller is able to control a subset C \ C ′ of contactors (e.g., those labelled with blue in Fig. 1). The state x of the system is defined as a valuation on all components n ∈ G ∪ R and uncontrollable contactors e ∈ C ′ ⊆ C. We define Ω as the set of all states, i.e all the different valuations of the components. The state x is unknown and hence modeled as a random variable X that can only be determined by sensor measurements mapped back to a set of possible states in which the circuit may be. The overall goal is to design a strategy the fault detection controller runs to adaptively estimate the discrete state of the circuit by taking "actions" (i.e., closing and opening controllable contactors), and then reading voltage sensor measurements.

B. Mathematical formulation

In this section we introduce the relevant notations used throughout the rest of the paper and mathematically formulate the estimation problem.

The state X of the circuit is modeled as a random variable. Data on component types and reliability levels can be used to build a probability measure P[x] on Ω. At the beginning of the state estimation process the system is in the (unknown) state x 0 ∈ Ω. We assume that faults in the system are independent, and that x 0 remains fixed during the estimation process. This is a reasonable assumption because the timescale of the estimation process is meant to be much smaller than the failure rates of the components and the timescales of the other controllers in the system.

For the controllable subset of contactors, there exists a set V of actions v that can be performed and a set Y of measurements y that can be observed. For an action v ∈ V, y = µ(v, x) ∈ (Y ) is the unique outcome of performing action v if the system is in the state x. The actions {v 0 , ..., v t } performed and outcomes {y 0 , ..., y t } observed up until step t are represented by the partial realization

ψ t = {(v i , y i )} i∈{0,...,t} . Given two partial realizations ψ t and ψ t ′ , we say that ψ t is a subrealization of ψ t ′ if ψ t ⊆ ψ t ′ .
At each step t, the probability measure P[x] can be updated by conditioning it on ψ t to obtain

P[x | ψ t ].
We are interested in an estimation process adaptively eliminating "invalid" states to get to the actual state x 0 . We define D(y, v), with y = µ(v, x 0 ), to be the set of states x ∈ Ω that are indistinguishable from x 0 under the action v. Formally,

D(µ(v, x 0 ), v) = {x ∈ Ω | µ(v, x) = µ(v, x 0 )}.
We further extend this concept by defining h(v 0:t , x 0 ), the set of states that produce the same set of outcomes {µ(v 0 , x 0 ), . . . , µ(v t , x 0 )} as x 0 under the same set of actions {v 0 , . . . , v t }. In the remainder of the paper, we use S t as a shorthand for h(v 0:t , x 0 ). If, at step t, we perform a new action v ′ / ∈ ψ t , there exists a recursive relation between the two sets of states:

h(v 0:t ∪ {v ′ }, x 0 ) = h(v 0:t , x 0 ) ∩ D(µ(v ′ , x 0 ), v ′ ), (1)
which leads immediately to

S t = ∩ i∈{0,...,t} D(µ(v i , x 0 ), v i ).
(

) 2
As only intersections are taken, the order of actions v i does not matter.

To represent the uncertainty in the state estimate, we define an objective function f : 2 V×Y × Ω → R + that maps the the set of actions A ⊆ V under state x 0 to reward f (A, x 0 ). A strategy π is a function from partial realizations to actions such that π(ψ t ) is the action v t+1 taken by π when observing ψ t . We denote Ṽ(π, x 0 ) ⊆ V the set of all the actions performed under the strategy π, the state of the system being x 0 . In the general case, Ṽ(π, x 0 ) = V.

The fault detection controller is assigned a budget k ≪ |V|, indicating the number of steps within which the estimation process should terminate. The system is initially in the state x 0 , which is fixed and unknown, and the controlled contactors are in some initial configuration v 0 . Initial configuration v 0 and the corresponding measurement y 0 constitute ψ 0 . Then, for i = 1, . . . , k, we consider the following process:

v i = π(ψ i-1 ) (3a) y i = µ(v i , x 0 ) (3b) ψ i = ψ i-1 ∪ (v i , y i ). (3c)
Equations ( 3a) -(3c) represent the decision making, measurement, and update in the estimation process, respectively. The goal is to reduce the uncertainty of X represented by the probability distribution P[x] through performing k actions. To that end, the following reward function is considered:

f (v 0:k , x 0 ) = -P[S k ] = - x∈S k P[x]. (4)
The behavior driven by the maximization of f is to remove as much probability mass from Ω as possible in k steps. It is also worth noting that when the underlying probability distribution on Ω is uniform, f is just proportional to the size of S k and so maximizing f is equivalent to minimizing the number of indistinguishable states.

The goal of estimation is to find the strategy that allows the "best expected estimate" for the state, i.e, the strategy π * such that

π * ∈ arg max π E[f ( Ṽ(π, X), X)], (5)
subject to | Ṽ(π, x)| k for all x, and with expectation taken with respect to P[x].

III. STRATEGY

In this section, we describe the algorithm used to solve the state estimation problem and give performance guarantees on the worst-case execution.

A. Greedy strategy

The optimal strategy for the fault detection controller would plan ahead for k steps. Complexity, however, scales up exponentially with k. To address the problem efficiently we develop a greedy strategy that selects, at each step, the action maximizing the expected one-step gain in uncertainty reduction. At step t, the greedy strategy uses the available information ψ t to compute the probability measure P[x | ψ t ] on the set S t , using a Bayesian update:

P[x | ψ t ] = P[ψ t | x] P[x] P[ψ t ] , ∀ x ∈ Ω. (6)
As the measurement process is deterministic, for a given x ∈ Ω we have P[x | ψ t ] = 1 {x∈St} , meaning that P[x | ψ t ] = 1 if x belongs to S t , and P[x | ψ t ] = 0 otherwise. From (6) we then get:

P[x | ψ t ] = P[x] P[ψt] ∀ x ∈ S t 0 elsewhere (7)
The term P[ψ t ] is the same for all x. It is a normalization coefficient that can be computed using

x∈St P[x | ψ t ] = 1 to obtain P[ψ t ] = x∈St P[x]. (8)
At each step t, the strategy consists of choosing the next action v t+1 that maximizes the gain in uncertainty reduction. Our measure of uncertainty comes from the value of the function f , established in Eq. ( 4), and therefore the benefit is expressed in terms of the change in f as we choose the action v. Consistent with our goal, we choose to maximize in mean the benefit at each step, the expectation taken with respect to the updated probability measure P[x | ψ t ]. We obtain the greedy strategy:

v t+1 ∈ arg max v∈V E[f (v 0:t ∪ {v}, X) -f (v 0:t , X) | ψ t ]. (9)

B. Worst-case performance guarantees

Greedy strategies in general can perform arbitrarily bad [START_REF] Bang-Jensen | When the greedy algorithm fails[END_REF]. However, by exploiting recent results on adaptive submodularity, we give a lower bound on the performance of the proposed strategy. For a brief overview of adaptive submodularity and related definitions, see Appendix A. We next show that the function f defined in Eq. ( 4), is adaptive monotone and adaptive submodular (Def. 2 and 3).

Proposition 1: The function f defined in Eq. ( 4) is adaptive monotone.

Proof: Given an action v ∈ V and partial realization ψ t at step t, we need to show the expected marginal benefit ∆(v|ψ t ) (see Def. 1) is nonnegative. For the cost function f , ∆(v|ψ t ) can be written as:

∆(v | ψ t ) = E[f (v 0:t , X) | ψ t ] -E[f (v 0:t ∪ {v}, X) | ψ t ].
(10) By Eq. ( 7), we get

∆(v | ψ t ) = x∈h(v0:t,x0) P[x|ψ t ] φ(x), (11)
with

φ(x) = x∈h(v0:t,x) P[x] - x∈h(v0:t∪{v},x) P[x]. (12)
By Eq. ( 1), h(v 0:t ∪ {v}, x) is a subset of h(v 0:t , x) for every x ∈ Ω. Thus, φ(x) 0, all the terms in the sum in Eq. ( 11) are non-negative, and ∆(v|ψ t ) 0.

Proposition 2: The function f defined in Eq. ( 4) is adaptive submodular.

Proof: Given in Appendix B. Theorem 1: For any true state x 0 ∈ Ω, the uncertainty reduction achieved in k steps by the greedy strategy given in Algorithm 1 is no worse than (1 -1/e) of what can be achieved in k steps by any other strategy, including the best possible strategy.

Proof: Follows directly from Propositions 1 and 2 and Theorem 2 given in Appendix A.

IV. IMPLEMENTATION

In this section, we give implementation details on the dynamic estimator employing the greedy strategy on some typical aircraft electric power system topologies. In order to reduce online computation, the inverse mapping from sensor measurements to compatible states of the circuit is conducted offline. Additionally, we propose some abstraction rules to reduce the size of the circuit as well as computation time.

A. Implementation details

The overall estimation process is summarized in Algorithm 1. In this algorithm, some items can be precomputed to improve run time. In particular, the inverse mapping from sensor measurements to compatible states does not have a closed form expression and the computation of the inverse map involves searching for paths on the graph G = (N , E) representing the circuit topology. Therefore, for all measurements y ∈ Y and all actions v ∈ V, the sets D(y, v) of states consistent with the action-measurement pairs (v, y) are computed offline to achieve a faster implementation. This collection is then accessed on the fly to significantly reduce the computation time as it is the most costly part of the algorithm.

Algorithm 1 Adaptive greedy strategy

Assumptions about the components and circuit can be easily incorporated in our framework. In particular, because these circuits are designed to achieve certain reliability levels, one common assumption is that at least one generator and one rectifier unit are online (delivering correct or improper voltage). These assumptions render certain states impossible, which are removed from the initial state set Ω.

B. Model reduction via abstraction

Although the greedy strategy provides an efficient way (with performance guarantees) to solve the dynamic state estimation problem, the offline computation for complex topologies can be very demanding as the number of possible states is exponential in the number of components whose states are being estimated. In this section, we give a set of rules that can be recursively applied to reduce the size of the circuit by clustering certain components together into metacomponents.

Components (generators, rectifier units, contactors) are connected through their ports to form the circuit, and sensors are placed on some of these ports. The main reduction idea is that when two uncontrolled components are connected together and there is no sensor on their internal connecting port, some of the individual states of the components may become indistinguishable from what can be measured with the available sensors. Therefore, they can be treated as a single basic component, called a metacomponent, having the same global overall behavior. It is then possible to hierarchically estimate the system state, first by estimating the state of the metacomponent, and then mapping this state to possible states of individual components forming the metacomponent. When running the greedy algorithm on the reduced circuit, the probabilities of metacomponent states should be adjusted accordingly to ensure a lossless abstraction.

The rules we use to simplify the circuits are summarized in Fig. 2. Figure 2(a), for example, shows how the combination of generator and contactor can be abstracted into a single "generator" metacomponent. For the original combination of components, the contactor can either be open (o) or closed (c), and the generator can either be healthy (h), unhealthy (u), or offline (o). Thus, the set Ω has six possible states, represented as a tuple of contactor status and generator health:

x 1 = (c, h), x 2 = (c, u), x 3 = (c, o), x 4 = (o,

!"

(a) Generator Metacomponent.

!"# !"# !"# !" !" (b) Rectifier Unit Metacompo- nent.

!" !"

(c) Contactor Metacomponent.

V. EXAMPLES

To assess the performance of the greedy strategy, we have systematically tested the dynamic estimator described in Algorithm 1 on multiple circuits typical of those found in electric power systems. For these experiments we have taken a uniform probability distribution over Ω, the reward function f defined in Eq. ( 4) becoming the size of the feasible set.

In many cases it is not possible to completely eliminate the uncertainty on the state of the system when there is a limited number of sensors. In order to evaluate the performance of the greedy strategy, we compare it with a brute force strategy, which exhaustively tries every action v ∈ V. Hence, no strategy can gather more information than the brute force strategy. Although the brute force strategy is not practically applicable, as |V| can be very large, it gives an upper bound on achievable performance, and can be used as a benchmark. Overall test methodology is summarized in Algorithm 2 and has been performed for every or some of the initial configurations v 0 ∈ V of the uncontrolled contactors.

Algorithm 2 Test methodology

Input: Initial configuration of the controlled contactors v0 ∈ V 1: for x0 ∈ Ω do 2:
Set the circuit in the state x0

3:

Put the controlled contactors in the configuration v 4:

Run the strategy tested (Greedy or brute force strategy) 5:

Record the computation time and the value of f at the end for statistics. 6: end for

A. Tests on a small-size circuit

The test small-size circuit, shown in Fig. 1 is comprised of 12 components. Six components are unknown (G1, G2, R1, R2, C2, C5), four contactors are controlled (C1, C3, C4, C6), depicted in blue in Fig. 1, and two voltage sensors are available (S1 and S2). The size of the state-space generated, taking into account the assumptions on faults, is 1600. A more precise description of the actual hardware circuit can be found in [START_REF] Rogersten | An aircraft electric power testbed for validating automatically synthesized reactive control protocols[END_REF]. On this particular example there are four controlled contactors, so the brute force strategy performs |V| = 2 4 = 16 actions. Both strategies have been run on the same Intel Core T M i3-2310M 64 bits CPU 2.10 GHz, 6.00 Gb RAM.

Simulation result: The greedy strategy with a horizon length of k = 6 performs as well as the brute force strategy, i.e., the value of the objective function f at the end of the 6 steps using the greedy strategy is the same as after the brute force strategy with 16 steps.

1) Average execution time: The average execution time for the greedy strategy is shown in Fig. 3. It takes on the order of milliseconds to compute the next best action to perform. On the other hand, the offline computation takes 30 seconds for this circuit.

2) Final value of the reward function f : Using our metric of performance and comparing greedy and brute force strategies, Fig. 4 shows that the greedy strategy performs as well as the brute force. A point at coordinates (n, m) signifies that in m% of the cases, there are n or fewer indistinguishable states after the strategy (greedy or brute force) stops. Starting from 1600 possible states, the greedy strategy reduces the number of candidates to less than 20 in 100% of the cases. In approximately half of the cases, there are four states or fewer that are still indistinguishable after the k = 6 steps. While the greedy strategy shows potential, real-world problems are bigger and more complex. To that end, the lossless abstraction process designed in Subsection IV-B proves useful for scaling up the strategy.

B. Tests on larger circuit

We also test the greedy strategy on a larger circuit representative of more-electric aircraft power distribution systems with multiple generators and demonstrate how abstraction can reduce the offline computation time. The circuit topology is shown in Fig. 5. Contactors controlled by the fault detection controller are depicted in blue.

Applying the lossless abstraction established in IV-B leads to a reduced circuit that eliminates four uncontrolled contactors. Comparing the offline computation for the full and reduced circuit, abstraction reduces the offline computation time from 3795 to 378 seconds, i.e, by a factor of 10.

We have tested the performance of the greedy strategy on this circuit with the methodology described in Algorithm 2. The results obtained on this test set by the greedy strategy (with k = 6 actions) were again equivalent to the results for the brute-force strategy (with |V| = 32 actions).

VI. CONCLUSIONS AND FUTURE DIRECTIONS

We proposed a greedy strategy for dynamic state estimation and fault detection in aircraft electric power systems. The proposed strategy was shown to have theoretical worstcase performance guarantees. Moreover, the practical performance of the strategy is well above the theoretical lower bound as demonstrated by examples.

The output of the overall estimation process is the set of all states consistent with observations, or the unique feasible state if possible. Future work will integrate these results with partial information games and synthesis of reactive control protocols. We are also interested in incorporating safety requirements while performing estimation. Although one can incorporate such information to disallow potentially unsafe actions, the main difficulty here is to obtain performance bounds when the set of available actions are changing in time. Currently, placement of sensors on the circuit topology is a given. By changing the number and locations of sensors, however, it may be possible to improve state estimation performance. Exploring this design space and trade-offs therein are likewise topics for future work.

APPENDIX

A. Background results in submodularity

We give some definitions and results on adaptive submodularity that follow the exposition provided in [START_REF] Golovin | Adaptive submodularity: Theory and applications in active learning and stochastic optimization[END_REF] and [START_REF] Krause | Submodular function maximization[END_REF]. Notations used here are defined in Subsections II-A and II-B.

Definition 1: Given an objective function f , an action v ∈ V, and a partial realization ψ t , ∆(v|ψ t ) is the conditional expected marginal benefit of v conditioned on having observed ψ t , defined as

∆(v|ψ t ) . = E[f (v 0:t ∪ {v}, X) -f (v 0:t , X)|ψ t ],
and the expectation taken with respect to P[x|ψ t ]. Definition 2: The function f : 2 V ×Y × Ω → R + is adaptive monotone with respect to distribution P[x] if the conditional expected marginal benefit of any action is nonnegative. Thus, for all v ∈ V and ψ t with P[ψ t ] > 0, ∆(v|ψ t ) ≥ 0. Definition 3: The function f : 2 V ×Y × Ω → R + is adaptive submodular with respect to distribution P[x] if the conditional expected marginal benefit of any fixed action v does not increase as more actions are performed and measurements are taken. Thus, f is adaptive submodular with respect to distribution P[x] if for all ψ t , ψ t ′ such that ψ t is a subrealization of ψ t ′ , and for all v ∈ V \ {v 0 , . . . , v t ′ }, ∆(v|ψ t ) ≥ ∆(v|ψ t ′ ). The adaptive greedy algorithm, a generalization of the greedy algorithm [START_REF] Krause | Submodular function maximization[END_REF], is a strategy that selects the action maximizing the conditional expected marginal benefit, conditioned on outcomes from all previous actions.

Theorem 2 (Theorem 1.14 in [START_REF] Krause | Submodular function maximization[END_REF]): Let π greedy l be a greedy strategy run for l iterations (so that it selects l actions). Let π * k be any policy selecting at most k actions for any realization x. Then,

f avg (π greedy l ) ≥ 1 -e -l/k f avg (π * k ),
where f avg (π) .

= E[f ( Ṽ (π, X), X)] is the expected reward of π.

In particular, by setting k = l we see that the greedy strategy selecting k items step by step obtains at least (1 -1/e) of the value of the optimal strategy that selects k items step by step.

B. Proof of proposition 2

We first state a lemma that will be useful in the proof.

(τ 1 , τ 2 , . . . , τ Y ) = Y i=1 τ i - Y i=1 τ 2 i Y i=1 τ i , (13)
is increasing on the positive orthant, i.e., b(τ

1 , τ 2 , . . . , τ Y ) ≥ b(s 1 , s 2 , . . . , s Y ) if τ i ≥ s i ≥ 0 for all 1 ≤ i ≤ Y .
Proof: Note that because b is symmetric, i.e., permutation invariant with respect to its arguments, it is enough to show that it is increasing in one of its arguments.

Let k 1 . = Y i=2 τ i and k 2 . = Y i=2 τ 2 i . Define b(x) . = b(x, τ 2 , . . . , τ Y ) = k 1 + x -k2+x 2
k1+x . The partial derivative of b with respect to state x is ∂ b/∂x = k 2 1 +k2 (k1+x) 2 , which is non-negative by definitions of k 1 and k 2 . Now, we are ready to prove Proposition 2. Consider two partial realizations ψ t and ψ t ′ s.t ψ t ⊆ ψ t ′ and the corresponding sets S t and S t ′ . Fix an action v ∈ V \ v 0:t ′ . To prove adaptive submodularity, ∆(v, ψ t ) can be expressed as a function dependent on the size of S t . We examine the variation of ∆ between S t and S t ′ .

Since the probability measure is non-uniform and can take values in some set {p 1 , . . . , p N }, we define the subsets of Ω where P[x] is constant:

We then separately compute the two terms in Eq. [START_REF] Keller | Aircraft electrical power systems prognostics and health management[END_REF]. First term becomes:

E[f (v 0:t , X) | ψ t ] = x0∈St P[x 0 | ψ t ]
x∈h(v0:t,x0)

P[x]. (17)
For x 0 ∈ S t , h(v 0:t , x 0 ) = S t , we obtain

E[f (v 0:t , X) | ψ t ] = y∈Y τ y . ( 18
)
For the second term in Eq. ( 10), we first get f (v 0:t ∪ {v}, x) = x∈h(v0:t,x)∩D(µ(v,x),x)

P[x]

= τ µ(v,x) .

From Eq. ( 14) and Eq. ( 16), we obtain:

E[f (v 0:t ∪ {v}, X) | ψ t ] =
where Y . = |Y|. This expression of ∆(v|ψ t ) in terms of the variables τ i is similar for the partial realization ψ t ′ ; the only change is the set S t , which is represented in the function b by a different value of the τ i denoted τ ′ i . Since ψ t ⊆ ψ t ′ and S t ′ ⊆ S t , τ i and τ ′ i satisfy τ ′ i τ i for all i. Therefore, adaptive submodularity is equivalent to showing that b is increasing on the positive orthant, and Lemma 1 concludes the proof.

Fig. 1 .

 1 Fig. 1. A single-line diagram of a simple circuit with AC components (in black) and DC components (in red).

Input: 6 : 8 :

 68 Probability measure P[x] on Ω, number of actions to perform k. The system is in the state x0 ∈ Ω, fixed and unknown, and the controlled contactors are in some configuration v0. Output: Partial realization ψ k and the set S k of compatible states after k actions are taken based on the strategy π greedy 1: Take the measurement y0 = µ(v0, x0). 2: ψ0 = {(v0, y0)} 3: for t ∈ {1, . . . , k} do 4: vt = π greedy (ψt-1) Take the measurement yt = µ(vt, x0) 7: ψt = ψt-1 ∪ {vt, yt} St = St-1 ∩ D(yt, vt) 9: Compute P[x | ψt] (Bayesian update) 10: end for 11: return (ψ k , S k , P[x | ψ k ])

  h), x 5 = (o, u), and x 6 = (o, o). The "generator" metacomponent, however, has three possible states, corresponding to healthy, unhealthy, and offline: x1 = h, x2 = u, and x3 = o. These metacomponent states can be mapped back to the corresponding original components, such that x1 = {x 1 }, x2 = {x 2 }, and x3 = {x 3 , x 4 , x 5 , x 6 }.

Fig. 2 .

 2 Fig. 2.Metacomponents used for abstraction. In terms of possible external behaviors (i.e., what can be measured from the external ports), two-component circuit units (shown in black) are equivalent to the single component units (shown in red).

Fig. 3 .

 3 Fig. 3. Histogram of execution time for the greedy strategy.

Fig. 4 .

 4 Fig. 4. Performance comparison between greedy and brute-force strategies.

Fig. 5 .

 5 Fig. 5. A single-line diagram of a larger circuit with AC and DC components.

Lemma 1 :

 1 The function b : R Y → R, defined as b

  F n = {x ∈ Ω | P[x] = p n } for n ∈ {1, . . . , N }. The collection F 1:N is trivially a partition of Ω. It is possible to show that the sets {D(y, v * ) ∩ F n |y ∈ Y, n ∈ 1 : N } form a partition of Ω and thus a partition of S t . Let α n,y . = S t ∩ D(y, v) ∩ F n . Then for all x ∈ α n,y , we have µ(v, x) = y and P[x] = p n .(14)By Eq. (8) , we get a new expression for P[ψ t ]:P[ψ t ] = x∈St P[x] = y∈Y n∈1:N p n |S t ∩ D(y, v) ∩ F n |. (15) Let τ y . = n∈1:N p n |α n,y |. Then, conditional probabilities on F n can be rewritten as ∀x ∈ F n , P[x | ψ t ] = p n y∈Y τ y .

x∈Stf.

  (v 0:t ∪ {v}, x)P[x | ψ t ]Finally, putting the two terms of Eq. (10) leads to∆(v|ψ t ) = b(τ 1 , τ 2 , . . . , τ Y ) =

ACKNOWLEDGMENTS

The authors would like to thank Jean-Michel Maillet and Eric Wolff for very useful and enlightening discussions. This work was supported in part by IBM and UTC through the iCyPhy consortium.


P Salzenstein
email: patrice.salzenstein@femto-st.fr

E Pavlyuchenko

A Hmima

N Cholley

M Zarubin

S Galliou

Y K Chembo

L Larger

Estimation of the uncertainty for a phase noise optoelectronic metrology system

A configuration of phase noise measurement system operating in X-band (8.2-12.4 GHz) using a photonic delay line as a frequency discriminator is presented in this paper. This system doesn't need any excellent frequency reference and works for any frequency in this band. Oscillator frequency fluctuation is converted to phase frequency fluctuation through the delay line. The measured phase noise includes the Device Under Test (DUT) noise and the instrument background. Then the use of a cross correlation decreases the cross spectrum terms of uncommon phase noise as √(1/m), where m is the average number. Using cross correlation on 500 averages, noise floor of the instrument £(f) becomes respectively -150 and -170 dBc/Hz at 10 1 and 10 4 Hz from the 10 GHz carrier (-90 and -170 dBc/Hz including 2 km delay lines). We then focus on determining uncertainty. There are two categories of uncertainties terms : "type A", statistic contribution such as repeatability and experimental standard deviation; "type B" due to various components and temperature control. Elementary term of uncertainty for repeatability is found to be equal to 0.68 dB. Other elementary terms still have lower contributions. Its leads to a global uncertainty of 1.58 dB at 2σ.

Introduction.

Recent progress in optoelectronic oscillators based on delay lines [START_REF] Yao | High frequency optical subcarrier generator[END_REF] or optics resonators with modulation frequency in microwaves [START_REF] Ilchenko | High-Q microsphere cavity for laser stabilization and optoelectronic microwave oscillator[END_REF][START_REF] Volyanskiy | Compact Optoelectronic Microwave Oscillators using Ultra-High Q Whispering Gallery Mode Disk-Resonators and Phase Modulation Optics[END_REF][START_REF] Tavernier | Magnesium Fluoride Whispering Gallery Mode Disk-Resonators for Microwave[END_REF][START_REF] Schliesser | Cavity Optomechanics with Whispering-Gallery Mode Optical Micro-Resonators[END_REF] require low phase noise system. A convenient method is to lock the oscillator to be measured on a reference standard. However it is only possible to determine the phase noise for oscillators compared to another one with the same frequency if noise is expected to be better than a synthesizer. That's why there is a need to develop methods based on the use of two delay lines and cross correlation to be able to characterize an oscillator without referring to another one with the same delivered frequency [START_REF] Salzenstein | Realization of a Phase Noise Measurement Bench Using Cross Correlation and Double Optical Delay Line[END_REF]. We present in this short paper the development made on an instrument dedicated to such measurement for metrology applications and the evaluation of the uncertainty.

Principle of the phase noise measurement system.

A quasi-perfect RF-microwave sinusoidal signal can be written as :

υ(t)=V 0 [1+α(t)]cos(2πν 0 t+φ(t))
Eq.( 1) where V 0 is the amplitude, ν 0 is the frequency, α(t) is the fractional amplitude fluctuation, and φ(t) is the phase fluctuation. Equation (1) defines α(t) and φ(t) in low noise conditions: |α(t))| <<1 and |φ(t)| <<1. Short-term instabilities of signal are usually characterized in terms of the single sideband noise spectral density PSD S(f). S is typically expressed in units of decibels below the carrier per hertz (dBc/Hz) and is defined as the ratio between the one-side-band noise power in 1 Hz bandwidth and the carrier power:

£(f) = ½ Sφ(f) Eq.
(2) This definition given in equation (2) includes the effect of both amplitude and phase fluctuations. However we must know the amplitude and phase noise separately because they act differently in the circuit. For example, the effect of amplitude noise can be reduced by amplitude limiting mechanism and mainly suppressed by using a saturated amplifier. Phase noise of microwave oscillators can usually be characterized by heterodyne measurement. Whereas, for such a system, we need a reference oscillator operating exactly at the frequency of the DUT with lower phase noise. Reference is no more required for homodyne measurement with a delay line discriminator. At microwave frequencies, electrical delay is not suitable because of its high losses. However photonic delay line offers high delay and low attenuation equal to 0.2 dB/km at the wavelength λ=1.55 μm. Optoelectronic phase noise measurement system is schematically represented on figure 1. It consists on two equal and fully independent channels. The phase noise of the oscillator is determined by comparing phase of the transmitted signal to a delayed replica through optical delay using a mixer. It converts the phase fluctuations into voltage fluctuations. An electro-optic modulator allows modulation of the optical carrier at microwave frequency. The length of the short branch where microwave signal is propagating is negligible compared to the optical delay line. Mixers are used as phase detectors with both saturate inputs in order to reduce the amplitude noise contribution. The low pass filters are used to eliminate high frequency contribution of the mixer output signal. DC amplifiers are low flicker noise. The oscillator frequency fluctuation is converted to phase frequency fluctuation through the delay line. If the mixer voltage gain coefficient is K φ (volts/radian), then mixer output rms voltage can be expressed as :

V² out (f) = K² φ │H φ (jf)│²S φ (f) Eq. ( 3
)
Where │H φ (jf)│² = 4.sin²(πfτ) is the transfer function of optical delay line, and f is the offset frequency from the microwave carrier. Equation [START_REF] Volyanskiy | Compact Optoelectronic Microwave Oscillators using Ultra-High Q Whispering Gallery Mode Disk-Resonators and Phase Modulation Optics[END_REF] shows that the sensitivity of the bench depends directly on K² φ and │H φ (jf)│. The first is related to the mixer and the second essentially depends on the delay τ . In practice, we need an FFT analyzer to measure the spectral density of noise amplitude V² out (f)/B, where B is the bandwidth used to calculate V out (f)/B. The phase noise of the DUT is finally defined by Eq. ( 4) and taking into account the gain of DC amplifier G DC as :

£(f)=[V² out (f)] / [2K² φ .│H φ (jf)│² G² DC B]
Eq. ( 4)

Validation of the performances.

The measured phase noise includes the DUT noise and the instrument background. The cross correlation method allows to decrease the cross spectrum terms of uncommon phase noise as √(1/m), where m is the average number. Thereby uncorrelated noise is removed and sensitivity of measure is improved. To validate the measure of our phase noise bench, we need to compare datasheet of the commercial frequency synthesizer Anritsu/Wiltron 69000B [7] with the phase noise we measure using our system. Figure 2 shows the result of this measure. We can see that our bandwidth is limited to 100 kHz (τ = 10 μs) and the measured phase noise corresponds to the datasheet. Figure 3 represents the background phase noise of the bench after performing 500 averaged with cross-correlation method, when removing the 2 km optical delay line. In this case, phase noise of the 10 GHz synthesizer is rejected. The solide curve shows noise floor (without optical transfer function) respectively better than -150 and -170 dBc/Hz at 10 1 and 10 4 Hz from the 10 GHz carrier. Dotted curve is the noise floor when optical fiber is introduced.

Evaluation of the uncertainty.

There are two categories of uncertainties terms [START_REF] Gum | Guide to the Expression of Uncertainty in Measurement, fundamental reference document[END_REF] : the first category is called "type A", statistic contribution such as reproducibility, repeatability, special consideration about Fast Fourier Transform analysis, and the experimental standard deviation. Second family of uncertainties contributions is "type B" and due to various components and temperature control. Experience with or general knowledge of the behavior and properties of relevant materials and instruments, manufacturer's specifications, data provided in calibration and other certificates (noted BR), uncertainties assigned to reference data taken from handbooks. By the way, uncertainty on £(f) strongly depends on propagation of uncertainties through the transfer function as deduced from equation (4). Global uncertainty is strongly related to repeatability of the measurements. Repeatability is the variation in measurements obtained by one person on the same item and under the same conditions. Repeatability conditions include: the same measurement procedure, the same observer, the same measuring instrument, used under the same conditions, repetition over a short period of time, the same location.

We switch on all the components of the instrument and perform a measurement keeping the data of the curve. Then we need to switch them off and switch them on again to obtain another curve. We must repeat this action several times until we have ten curves. Elementary term of uncertainty for repeatability is found to be equal to 0.68 dB. Other elementary terms still have lower contributions. For instance, temperature effects are less than 0.1 dB as optical fiber regulation system of temperature in on. Resolution of instruments is determined by the value read on each voltmeter when we need to search minimum and maximum for the modulator but also for wattmeter. Resolution is then no worse than 0.1 dB. Uncertainty on the determination of the coefficient K φ dependant for the slope expressed in Volt/rad is found to be lower than 0.08dB. For the contribution of the use of automatic/manual range, we can deduce from the curves that this influence is no more than 0.02 dB. In our case, all these terms were found lower than repeatability. Other elementary terms still have lower contributions. Its leads to a global uncertainty of 1.58 dB. Table 1 summarizes each contribution of elementary term of uncertainty.

Conclusion and further work.

We detailed performances and consideration about estimation of the uncertainty to show the main advantage of the instrument developed for metrology applications. With high performance better than -150 dBc/Hz at 10 Hz from the 10 GHz carrier, integration in the accredited laboratory is expected and the new field of capabilities in terms of calibration of single source in X-band will be soon achieved to extend the Calibration Metrology Capabilities (CMC) of the LNE (french national metrology institute). It is interesting to underline that it is possible to determine the phase noise of a single oscillator in X-band with a global uncertainty set to be better than ±2 dB.

Figure 1 .

 1 Figure 1. Phase noise bench.

Figure 2 .Figure 3 .

 23 Figure 2. Phase noise (dBc/Hz) of the synthesizer measured at 10 GHz with K φ =425 mV/rad and G DC = 40dB versus Fourier Frequency

Table 1 .

 1 Budget of uncertainties

	Uncerta inty	Designation	Value (in dB)
	A1	Repeatability	0.68
	A2	Reproductibility	0
	A3	Uncertainty term due to the number of sample	0.1
	A	(ΣAi²) 1/2	0.69
	BR	Not applicable	0
	BL1	Gain of the DC amplifier	0.04
	BL2	Influence of the temperature	0.1
	BL3	Influence of the resolution of the instrument	0.1
	BL4	Influence of the power of the DUT	0
	BL5	Uncertainty on the determination of K Φ	0.08
	BL6	Contribution of automatic/manual range	0.02
	BL7	Influence of the variation of the input power	0.02
	BL	ΣBLi	0.38
	uC	(A²+BR²+BL²) 1/2	0.79