Bogdan Robu

Lucie Baudouin
email: baudouin@laas.fr

Christophe Prieur
email: christophe.prieur@gipsa-lab.grenoble-inp.fr.

Denis Arzelier
email: arzelier@laas.fr

Simultaneous H ∞ vibration control of fluid/plate system via reduced-order controller

Keywords: Partial differential equation, flexible system, fluid/plate system, H ∞ control, HIFOO

We consider the problem of the active reduction of structural vibrations of a plane wing induced by the sloshing of large masses of fuel inside partly full tank. This study focuses on an experimental device composed of an aluminum rectangular plate equipped with piezoelectric patches at the clamped end and with a cylindrical tip-tank, more or less filled with liquid at the opposite free end. The control is performed through piezoelectric actuators and the main difficulty comes from the complex coupling between the flexible modes of the wing and the sloshing modes of the fuel. First, a partial derivative equation model and then a finite -dimensional approximation, calculated using the first 5 structural modes of the plate and the first 2 liquid sloshing modes, is established. Second, after a model matching procedure, a simultaneous H∞ control problem associated to the vibration attenuation problem for two different fillings of the tank is stated. Due to the large scale of the synthesis model and to the simultaneous performance requirements, a reduced-order H∞ controller is computed with HIFOO 2.0 package and is compared with individual designs for different filling levels. Experimental results are finally provided illustrating the relevance of the chosen strategy.

Introduction

As commercial transport aircraft designs become larger and more flexible, the impact of aeroelastic vibration of the flight dynamics grows in prominence. See e.g. [START_REF] Winther | Reduced-order dynamic aeroelastic model development and integration with nonlinear simulation[END_REF], [START_REF] Becker | Active buffeting vibration alleviation demonstration of intelligent aircraft structure for vibration & dynamic load alleviation[END_REF], [START_REF] Buttrill | The impact of structural vibration on flying qualities of a supersonic transport[END_REF] or even [START_REF] Kubica | European patent office: "procedé et dispositif pour réduire les mouvements vibratoires d'un fuselage d'un aéronef[END_REF] for airplanes and see [START_REF] Shaw | Active control of the helicopter rotor for vibration reduction[END_REF] for helicopter dynamics. Moreover for space applications, the interaction of flexible modes and sloshing modes disturbs the dynamics (see [START_REF] Deyst | Effects of structural flexibility on entry vehicle control systems[END_REF]).

For applications in aeronautics, it is crucial to attenuate the plane wing's vibrations when the wing is in interaction with the movement of the fuel inside of it. Smart materials are used for many applications for instance in civil engineering. Thus flexible structures, which are equipped with piezoelectric patches, occupy a major place in the control research area. Their capability of attenuating the vibrations and measuring the deformation is described in [START_REF] Becker | Comparison of piezoelectric systems and aerodynamic systems for aircraft vibration alleviation[END_REF][START_REF] Bhikkaji | Multivariable integral control of resonant structures[END_REF][START_REF] Fleming | Optimal impedance design for piezoelectric vibration control[END_REF] and [START_REF] Sun | Vibration suppression of a hard disk driver actuator arm using piezoelectric shunt damping with a topology-optimized PZT transducer[END_REF] among other references. In the present paper, it is shown how piezoelectric devices can be useful for vibration control in aeronautics.

We are studying the experiment depicted by Figure 1 which is an example of a coupled fluidflexible structure system and has the same first flexible modes as a plane wing [START_REF] Pommier-Budinger | Active control of a structure with sloshing phenomena[END_REF]. For vibration's Figure 1: Plant description: the rectangular plate and the horizontal cylinder attenuation of flexible structures, some studies investigate the use of piezoelectric patches to effectively suppress the vibrations (see [START_REF] Baillargeon | Active vibration suppression of sandwich beams using piezoelectric shear actuators: Experiments and numerical simulations[END_REF], [START_REF] Denoyer | Dynamic modelling and vibration suppression of a swelling stucture utilizing piezoelectric sensors and actuators[END_REF], [START_REF] Hwang | Robust LQG control of an all-clamped thin plate with piezoelectric actuators/sensors[END_REF], [START_REF] Yamada | Hybrid vibration suppression of multiple vibration modes of flexible structures using piezoelectric elements and analog circuit[END_REF], [START_REF] Zhang | Active vibration control of flexible structures using piezoelectric materials[END_REF] among others). However, only a few results are already available in the literature for fluid-structure systems. Reference [START_REF] Lasiecka | Boundary feedback control in fluid-structure interactions[END_REF] gives a recent theoretical result and [START_REF] Terasawa | Adaptive identification of MR damper for vibration control[END_REF] validates the method by means of experimental results. To the best of our knowledge, there are less studies of fluid-structure system dedicated to aerospace applications, except [START_REF] Pommier-Budinger | FEM design of a piezoelectric active control structure. application to an aircraft wing model[END_REF] or [START_REF] Pommier-Budinger | Fractional robust control with iso-damping property[END_REF] where controllers are designed using a numerical model.

To take the large numbers of degrees of freedom in the dynamics into account, we focus on an infinite-dimensional model described by partial differential equations. Many results are available about the suppression of vibrations using a distributed parameter model. Let us cite [START_REF] Demetriou | An adaptive control scheme for a class of second order distributed parameter systems with structured perturbations[END_REF] or [START_REF] Potami | Actuator switching for vibration control of spatially distributed systems[END_REF] among other references.

The first contribution of this paper is the construction of an infinite-dimensional model for the coupled fluid-flexible structure system of Figure 1. We consider a plate equation to describe the dynamics of the flexible wing, as in [START_REF] Lee | Vibration control experiments on a piezoelectric laminate plate using spatial feedforward control approach[END_REF], in an infinite-dimensional context, and Bernoulli equation for the dynamics of the fluid. Other dynamical models and associated control problems have been considered in the fluid control literature. In [START_REF] Coron | Local controllability of a 1-D tank containing a fluid modeled by the shallow water equations[END_REF], the author studies the Shallow Water equations for modeling of the dynamics of the fluid while using the return method and directly controlling the tank's displacements. In [START_REF] Prieur | Stabilization of a 1-D tank containing a fluid modeled by the shallow water equations[END_REF], a Lyapunov method is applied to the same nonlinear equation for the fluid while [START_REF] Petit | Dynamics and solutions to some control problems for water-tank systems[END_REF] proposes to compute a flat output for the linearized model. An identical control problem is also considered in [START_REF] Grundelius | Constrained iterative learning control of liquid slosh in an industrial packaging machine[END_REF][START_REF] Grundelius | Iterative optimal control of liquid slosh in an industrial packaging machine[END_REF] where the iterative learning control is applied to the linear Bernoulli model. In the present paper, the setting is quite different since water-tank oscillations are indirectly controlled through a flexible structure (more precisely through a flexible plate) by means of piezoelectric actuators. Moreover, in addition to the stabilization problem, we are interested in minimizing a performance criterion based on the H ∞ norm of some transfer.

The geometry of the tank has been chosen so that the mock-up has the same first flexion/torsion undamped natural frequencies as a real airplane wing filled with liquid. Considering the difficulty to properly derive a model for a cylindrical tank, an approximation by a rectangular tip-tank and by a system of masses-pendulum is used. Moreover, due to the physical interaction of the flexible plate with the fluid, a coupling between the plate equation and the dynamical fluid model needs to be considered. It is detailed here in both infinite and finite dimensions. The finite dimension approximation is then modified in order to match with the measures observed on the experimental setup.

The second contribution is to define and solve a robust control problem for this fluid-flexible structure system. The structure is subject to a disturbance created by the voltage applied to one piezoelectric actuator and the control problem aims at rejecting the perturbation that may create some vibrations of the coupled fluid/plate system. Robust H ∞ controllers are then computed and it is shown that we need to consider reduced order controllers. To do so, a non smooth optimization approach proposed in the HIFOO package is used [START_REF] Gumussoy | Multiobjective robust control with HIFOO 2.0[END_REF]. It has also the interest of allowing the computation of a unique controller for several fillings of the tank, by solving a simultaneous H ∞ optimal control problem. Finally, some experiments on the real setup are presented and the relevance of the stabilizing controller to suppress the vibrations of the fluid-structure system is shown. To avoid a spillover effect when closing the loop, a suitable filter is also defined in the standard control problem (see [START_REF] Hwang | Robust LQG control of an all-clamped thin plate with piezoelectric actuators/sensors[END_REF]).

The paper is organized as follows. The plant under consideration, equipped with piezoelectric patches (sensors and actuators) along with the control objectives are detailed in Section 2. In Section 3 the fluid/plate model of the system is established. In Section 4, first the model matching problem is solved, second a full-order robust controller and then a reduced-order one using HIFOO package are computed. Finally, the effectiveness of the robust controller of reduced order is checked on various experiments. Section V comprises some concluding remarks and Section VI some technical developments.

2 Problem statement and control objectives

Plant description

The plant to be controlled is located at ISAE-ENSICA, Toulouse, France and has been constructed to have the vibration frequencies of a real plane wing filled with fuel in its tip-tank [START_REF] Pommier-Budinger | Active control of a structure with sloshing phenomena[END_REF].

The device is composed of an aluminium rectangular plate and a plexiglas horizontal cylindrical tip-tank filled with liquid (see Figures 1 and2). The two piezoelectric actuators made from PZT (Lead zirconate titanate) are bonded next to the plate clamped side. Two sensors (made from PVDF -Polyvinylidene fluoride), are located on the opposite side of the plate with respect to the actuators. The characteristics of the collocated sensors and actuators are given in Table 2.

The tank is centered at 1.28 m from the plate clamped side and is symmetrically spread along the horizontal axis. Due to the configuration of the whole system (see Figure 2), the tank undergoes The experimental setup is depicted in Figure 3. The system (plate + tank) is connected to the computer via a high voltage amplifier delivering ±100 V and a charge amplifier 2635 from Brüer & Kjaer. Moreover, the controller is implemented in a DSpace c board.

First, the voltage delivered by the DSpace c card is amplified by the high voltage amplifier and then applied to the piezoelectric actuator. The deflection of the beam is measured by the piezoelectric sensor and then transmitted to the charge amplifier which will deliver a voltage to the DSpace c card.

In order to obtain the Bode plots, tests are made using a sampling time of 0.004 seconds on the DSpace c card and a Dynamic Signal Analyzer.

Control objectives

The main control objective is to attenuate the vibrations of the plate and the sloshing modes of the liquid in the tank. The control input is the voltage applied to one of the two piezoelectric actuators. The flexible structure is subject to a disturbance created by a voltage of frequency varying between 0 Hz and 50 Hz applied to the other piezoelectric actuator. This perturbation is a source of vibrations of the fluid-structure system and is modeled by a low-pass filter of order 1 with a bandwidth of 50 Hz. The filter before the piezoelectric actuator used as a disturbance actuator has therefore a transfer function given by H 1 (s) = 100π s+100π . The output of the system is the voltage measured with a piezoelectric sensor collocated with the piezoelectric control actuator. The to-be-controlled output is composed of the output of the system and the control input. The objective of the controller is to minimize the maximum of the frequency response of the transfer function between the control input and this output. This specification is formulated as minimizing the H ∞ norm of the transfer between the input disturbances and the controlled output.

Modeling of the system

In this section, a mathematical model of the system presented above is derived. We first describe an infinite-dimensional model using partial derivative equations and then explain how a finitedimensional model is obtained for numerical implementation by considering only the first modes of vibration and sloshing.

Concerning the horizontal tank, it is known (see e.g. in Chapters 1 and 2 of [START_REF] Ibrahim | Liquid sloshing dynamics[END_REF]) that the solution of the sloshing problem depends on the geometry of the tank. In the case of an horizontal cylindrical tank, the difficulty resides in the fact that the shape of a partially filled horizontal cylinder does not fit into any standard coordinate system and thus the velocity potential, which quantifies the sloshing, cannot be derived using separation of variables method. In our situation, where the tank undergoes longitudinal movement, the mode frequencies are simply curves faired through experimental data since there is no analytical result available (see [START_REF] Dodge | The new "dynamic behavior of liquids in moving containers[END_REF]Chapter 1.6]) for this particular case.

In order to overcome this difficulty, we use a geometrical approximation. The cylindrical tank is replaced by an "equivalent" rectangular one with the same sloshing frequencies (where length and width are respectively denoted a and b and are respectively along the x-axis and y-axis). Analytical results are available in the literature for the calculation of modes and forces/moments (see [START_REF] Ibrahim | Liquid sloshing dynamics[END_REF]Chapter 1.6]). The parameters of the rectangular "virtual" tank are computed as follows:

The length is equal to the one of the cylindrical tank while the width is equal to the width of the free surface of the liquid for the corresponding tank filling. Finally, the height of the liquid in the tank is defined for the liquid volume in both containers to be the same. The height of the rectangular tank is chosen in order to keep the same tank filling ratio (e) in both containers. Furthermore, one can see from [START_REF] Ibrahim | Liquid sloshing dynamics[END_REF]Chapter 1.5] that the tank height has no influence on the computation of the sloshing frequencies. Technical details about the approximation method can be found in Section 6.1 of the Appendix.

After comparing the sloshing frequencies of the rectangular "virtual" tank (using the method detailed earlier) to the ones calculated from the experimental data, a difference of only about 10 -2 Hz is observed.

Infinite-dimensional model

On the one hand, considering the plate, we assume that the kinematic hypothesis given by [START_REF] Blevins | Formulas for natural frequency and mode shape[END_REF]Chapter 11.1] is respected. This assumption allows us to derive the partial derivative equation of the plate (see [START_REF] Géradin | Mechanical Vibrations: theory and application to structural dynamics[END_REF]Chapter 4.6]):

m s ∂ 2 w ∂t 2 + ζ(w) ∂w ∂t + Y I s ∆ 2 w = ∂ 2 m y ∂y 2 + ∂ 2 m z ∂z 2 (1)
where ζ(w) is an operator quantifying the damping, w = w(y, z, t) is the displacement, m s the mass per unit plate area, Y the Young modulus and I s =

h 3 12(1-ν 2 )
the moment of inertia of the plate. ∆ is the Laplace operator with ∆ 2 being equal to

∂ 2 ∂y 2 + ∂ 2 ∂z 2 2
. The notations m y , m z correspond to external moments, along the y and z-axis, delivered to the plate by the actuators (see [START_REF] Dimitriadis | Piezoelectric actuators for distributed vibration excitation of thin plates[END_REF] or [START_REF] Crépeau | Control of a clamped-free beam by a piezoelectric actuator[END_REF]) and by the sloshing modes of the liquid in the tank. The moments delivered by the actuators are denoted by m a y and m a z and those delivered by the fluid by m f y and m f z . In other words, we have:

m y = m a y + m f y , m z = m a z + m f z
Equation ( 1) is to be solved using the boundary conditions given in Section 6.2 for given initial conditions, as in [5, Chapter 8.1.1] or [START_REF] Crépeau | Control of a clamped-free beam by a piezoelectric actuator[END_REF].

On the other hand, we now consider the longitudinal movement of the liquid along the x-axis. Because the liquid motion is starting from rest, there is a velocity potential φ(x, y, z, t) such that (see [START_REF] Sir | Hydrodynamics[END_REF]Chapter 1.12]) the equation of liquid continuity may be written as:

∂ 2 φ ∂x 2 + ∂ 2 φ ∂y 2 + ∂ 2 φ ∂z 2 = 0 (2)
The linearized (Bernoulli) equation of liquid motion is given by [START_REF] Sir | Hydrodynamics[END_REF]Chapter 2.20] or [START_REF] Khandelwal | A mechanical model for liquid sloshing in a rectangular container[END_REF]:

∂φ ∂t + p ρ + g(z -h) -C 0 x = 0 (3)
where C 0 stands for the acceleration along the x-axis, g for the gravitational acceleration and h for the liquid height in the container at rest position; ρ is the density of the liquid and p(x, y, z, t) is the pressure in the liquid. The equation is linearized by omitting the square products of the velocity potential (see [START_REF] Sir | Hydrodynamics[END_REF]Art. 20] for the complete equation before linearisation). Further on, (2) and (3) are to be solved using the boundary conditions detailed in Section 6.2 of the appendix (more details are also given in [25, Chapter 1.9] for the rectangular tank).

As in [13, Chapter 1.2], we denote the height of the oscillating free surface z s by:

z s = h + δ(x, y, t) (4)
where δ(x, y, t) is the small displacement of the free surface above the undisturbed mean level z = h.

As in [START_REF] Khandelwal | A mechanical model for liquid sloshing in a rectangular container[END_REF], the external acceleration C 0 can be calculated at the gravity center G = (y G , z G ) of the tank for steady motion:

C 0 = ẅ(y G , z G , t). (5)
In order to complete the description of the coupling between the two pde's, let us now compute the moment m f = (m f x , m f y , m f z ) applied by the liquid on the plate. Since it is assumed that there is no displacement along the y-axis, we get m f x = m f z = 0, and

m f y = b 0 h 0 (z - h 2 ) (p| x=a -p| x=0 ) dydz - a 0 b 0 (x - a 2 )p| z=0 dxdy. (6)
The second contribution to the whole moment is due to the voltage applied by the piezoelectric patches (see references [START_REF] Dimitriadis | Piezoelectric actuators for distributed vibration excitation of thin plates[END_REF][START_REF] Kim | Influence of piezo-actuator thickness on the active vibration control of a cantilever beam[END_REF]). As our experimental setup is not symmetric with respect to the x-axis (the actuators are only on one side of the plate), the position of the neutral fiber needs to be recalculated. We finally get that:

m a y (y, z, t) = m a z (y, z, t) = K a V a (t) if (y, z) ∈ A = 0 else ( 7
)
where A = (y a1 , y a2 ) × (z a1 , z a2 ) is the square delimited by the surface of the actuator and K a is a parameter depending on the actuator/plate geometry and on the position of the neutral fiber (see [START_REF] Kim | Influence of piezo-actuator thickness on the active vibration control of a cantilever beam[END_REF]).

The piezoelectric sensor delivers a voltage V c given by:

V c = K c yc2 yc1 zc2 zc1 ( ∂ 2 w ∂y 2 + ∂ 2 w ∂z 2 )dydz (8)
where K c is a parameter depending on the piezoelectric characteristics and on the capacity of the charge amplifier. y c1 , y c2 , z c1 , z c2 define the extremities of the piezoelectric sensor.

The perturbation entry, defined as a noise on the perturbation actuator glued next to the controlled actuator, creates a momentum similar to [START_REF] Crépeau | Control of a clamped-free beam by a piezoelectric actuator[END_REF]:

m p y (y, z, t) = m p z (y, z, t) = K p W p (t) (9)
where W p is a noise with a frequency of 50Hz.

Finite-dimensional approximation

In order to compute an approximation of the previous infinite-dimensional model, it is needed to describe an appropriate basis in which the solution of the partial derivative equation may be written. This issue is well covered in the literature (see for example [START_REF] Géradin | Mechanical Vibrations: theory and application to structural dynamics[END_REF]Chapter 5.2] for the plate and [START_REF] Khandelwal | A mechanical model for liquid sloshing in a rectangular container[END_REF] for the liquid sloshing) thus it will be detailed only in the appendix (see Section 6.3 below). This allows us to compute a finite-dimensional approximation of the model. Knowing that the influence of system's modes is inversely proportional to the mode's frequency, the approximation of the displacement by the first N modes of the system may be justified. The derivation of the finite-dimensional approximation of our model is now detailed. We aim at describing a state space model where the state vector should gather the first N vibration modes of the plate and the first M = N sloshing modes of the tank.

Plate model

In this section, we calculate the finite-dimensional approximation of size N of (1) for the rectangular plate, doing a truncation of the equation [START_REF] Prieur | Stabilization of a 1-D tank containing a fluid modeled by the shallow water equations[END_REF] at n = N .

The approximation of the influence of the tank on the plate will be considered in Section 3.2.3. We write:

Ẋp (t) = A p X p (t) + B p u(t) y(t) = C p X p (t) (10)
where

X p = q1 ω 1 q 1 • • • q N ω N q N -
the q n being the same as in (33) -is the state space vector. u(t) stands for the voltage applied to the piezoelectric patch and y(t) is the voltage delivered by the piezoelectric sensor.

In [START_REF] Deyst | Effects of structural flexibility on entry vehicle control systems[END_REF], the dynamic matrix A p is written as:

A p =     A p1 0 • • • 0 0 A p2 • • • 0 • • • 0 0 • • • A p N     (11)
with A pn = -2ζ n ω n -ω n ω n 0 for every n from 1 to N .

The angular frequency of the n th mode is ω n in [rad s -1 ] given in [5, Chapter 11] and ζ pn is the damping factor. The experiments have shown that the damping is not constant for every vibration mode but it depends on the quality factor which is different for every structure mode. For a given input voltage, we measure the quality vector and then identify the damping by analogy to a second order system with damping [START_REF] Siebert | Circuits, Signals, and Systems[END_REF].

For the control matrix B p , we use ( 1) and ( 7), and we integrate on the surface A of the actuator. Thus we get:

B p = (b p1 , 0, ..., b pn , 0, ..., b p N , 0) T ( 12
)
where

b pn = K b (Y ′ in (y a2 ) -Y ′ in (y a1 )) za2 za1 Z jn (z)dz +K b (Z ′ jn (z a2 ) -Z ′ jn (z a1 )) ya2 ya1
Y in (y)dy

Here i n and j n are the specific deformation modes of the plate along y and z-axis. Moreover, the Heaviside step is used here to suggest that the moment generated by the actuator is different than zero only under the actuator. More details of the derivation are given in Section 6.4 of the appendix. The output matrix C p is calculated with [START_REF] Demetriou | An adaptive control scheme for a class of second order distributed parameter systems with structured perturbations[END_REF] and the new position of the neutral fiber:

C p = (0, c p1 , ..., 0, c p k , ..., 0, c p N ) (13)
where

c pn = Kc ωp n (Y ′ in (y c2 ) -Y ′ in (y c1 )) zc2 zc1 Z jn (z)dz +(Z ′ jn (z c2 ) -Z ′ jn (z c1 )) yc2 yc1
Y in (y)dy and K c is the same parameter depending on the sensor/plate geometry and on the position of the neutral fiber as in [START_REF] Demetriou | An adaptive control scheme for a class of second order distributed parameter systems with structured perturbations[END_REF]. Here also, the Heaviside step is used to infer that the sensor senses the plate deformation which is only under it. The computation method is similar to the one of the control matrix from Section 6.4.

Liquid model

In this section, we calculate the finite-dimensional approximation of the tank's sloshing using a mechanical approach.

With the geometrical parameters of the rectangular tank, one can calculate the movement equations of the liquid and the boundary conditions (see [START_REF] Sir | Hydrodynamics[END_REF]Chapter 17]). The approximation given by [START_REF] Khandelwal | A mechanical model for liquid sloshing in a rectangular container[END_REF] may be used to model the sloshing. Indeed we consider a finite number of mass pendulum systems, each corresponding to a sloshing mode, which will exert on the tank the same force and moment as the sloshing liquid does (one can see for instance that for each odd integer i of ( 36) or [START_REF] Rubin | Improved component-mode representation for structural dynamic analysis[END_REF] in the appendix corresponds to a mass-pendulum systems of index k = 1

• • • M ).
The state-space representation is then easily obtained using the pendulum equation under external acceleration [START_REF] Ibrahim | Liquid sloshing dynamics[END_REF]:

θk + 2ξ θ g l k θk + g l k θ k = - 1 l k C 0 ( 14
)
where θ m is the angle of the m th pendulum compared to its equilibrium position, l k its length and ξ θ is the value of the damping fixed at 0.01. Choosing the state space vector for the liquid sloshing equal to

X θ = θ1 g l1 θ 1 • • • θ M g l M θ M T ,
the dynamic equation, for the general case of an input u acc , is:

Ẋθ = A θ X θ + B θ u acc ( 15
)
where the matrix A θ computed from ( 14) for each k is:

A θ =     A θ1 0 • • • 0 0 A θ2 • • • 0 • • • 0 0 • • • A θ M     (16)
with

A θ k =   -2ξ θ g l k -g l k g l k 0   .
In the case of the control matrix B θ the construction is also straightforward by considering u acc = C 0 as the control variable:

B θ = (b θ1 , 0, ..., b θ k , 0, ..., b θ M , 0) T ( 17
)
where b θ k = -1 l k 0 .

Model coupling

Plate influence on the rectangular tank The influence of the plate on the movement of the fluid is given by ( 5) which is expressed using the Ritz functions as (see [START_REF] Prieur | Stabilization of a 1-D tank containing a fluid modeled by the shallow water equations[END_REF]):

C 0 = N n=1 Y in (y G )Z jn (z G ). qn (t) (18)
Taking

K G = (Y i1 (y G )Z j1 (z G ), 0, • • • , Y i N (y G )Z j N (z G ), 0)
this equation is equivalent to:

C 0 = K G q1 (t) ω 1 q1 (t) • • • q N (t) ω N q N (t) T . ( 19
) Since q1 (t) ω 1 q1 (t) • • • q N (t) ω N q N (t) T =
Ẋp is the derivative of the state space vector of the rectangular plate, using [START_REF] Deyst | Effects of structural flexibility on entry vehicle control systems[END_REF] we obtain

C 0 = K G A p X p + K G B p u. (20)
Using now ( 16), ( 17) and [START_REF] Hwang | Robust LQG control of an all-clamped thin plate with piezoelectric actuators/sensors[END_REF], equation ( 15) becomes:

Ẋθ = A θ X θ + B θ (K G A p X p + K G B p u). (21)
This equation describes the coupling between the movement of the plate and the movement of the pendulum mass systems, that is to say the movement of the liquid in the tank.

Tank's influence on the rectangular plate The liquid sloshing is viewed by the plate as a perturbation that comes into the state space representation by means of a matrix A θp :

Ẋp = A p X p + B p u + A θp X θ y = C p X p (22)
(see paragraph 3.2.1). For the calculation of the perturbation matrix A θp , an identical approach as in the case of the control matrix B p has been followed. The external moment applied to the tank due to the movement of the fluid is computed in [START_REF] Schotte | Various modelling levels to represent internal liquid behavior in the vibration analysis of complex structures[END_REF]. Due to our choice of the physical parameters of the pendulum, the moment from the liquid sloshing equals the moment due to the movement of all mass-pendulum systems. The external moment written for M mass-pendulum system is

M θp = M k=1 m k L k l k θk
where l k θk is the total acceleration due to pendulum oscillations and plate movement. Using ( 16), the expression of M θp can be rewritten as:

M θp = (m 1 l 1 L 1 , 0, • • • , m M l M L M , 0) A θ X θ . (23)
In similar way as the computation of the control matrix B p given by ( 12), using [START_REF] Coron | Local controllability of a 1-D tank containing a fluid modeled by the shallow water equations[END_REF] and integrating on a small square (y 1G , z 1G ) × (y 2G , z 2G ) around the point G, we get that the matrix A θp ∈ R 2N ×2M is given by:

A θp = a T θp1 0 . . . a T θpn 0 . . . a T θp N 0 T (24)
where a θpn ∈ R 2M is given by

a θpn = K θp (Y ′ in (y 2G ) -Y ′ in (y 1G )) z 2G z 1G Z jn (z)dz +(Z ′ jn (z 2G ) -Z ′ jn (z 1G )) y 2G y 1G
Y in (y)dy and

K θp = (m 1 l 1 L 1 , 0, • • • , m k l k L k , 0, • • • , m M l M L M , 0) A θ .

Complete Model

From ( 21) and ( 22), one can write the complete state-space model of the experimental setup.

As previously stated, in the presentation of the complete model, we consider the case of N Ritz functions and M mass-pendulum systems. For the complete system, the state space vector is:

X = X p X θ .
Combining all the matrices leads to the following finite-dimensional approximation using the statespace representation:

           Ẋ = A p A θp A pθ A θ X + B p B pθ u + B w 0 w z = C p 0 0 0 X + 0 1 u y = C p 0 X. (25)
where

A pθ = B θ K G A p , B pθ = B θ K G B p
and 0 denotes null matrices of appropriate dimensions. B w is the perturbation matrix computed from (9) using the same methodology as B p .

As stated earlier in [START_REF] Sir | Hydrodynamics[END_REF], y stands for the measured output (the voltage measured through the piezoelectric sensor), u is the control input (the voltage applied to the piezoelectric actuator), z stands for the output to be controlled and is composed in this case of y and u. Finally, w is the disturbance input (not to be confused with the plate's displacement) and is generated by the voltage applied to the piezoelectric actuator used as source of disturbances.

Controller synthesis

The system model ( 25) is established for the first N modes of the plate and for the first M modes of the liquid sloshing. As it is observed in Section 4.1 below, this choice allows to consider the flexional and torsional deformation of the plate. Besides, it is shown in [START_REF] Halim | An optimization approach to optimal placement of collocated piezoelectric actuators and sensors on a thin plate[END_REF] that the first modes contain the main part of the energy of the deformation of the flexible structure. Moreover, using the energy approach from [START_REF] Tliba | H∞ control design for active vibration damping of flexible structure using piezoelectric transducers[END_REF], it is possible to check that the first five modes contain almost all the energy of the structure. Since the objective is to control the flexional deformation but also the torsional modes of the plate, the number of modes considered is also chosen so that torsion modes are included. Therefore the controller will be computed with N = 5 and M = 2. So X p ∈ R 2N and X l ∈ R 2M and the total state space vector X belongs to R 14 .

Model matching

The complete system model from [START_REF] Sir | Hydrodynamics[END_REF] was validated in [START_REF] Robu | A controlled distributed parameter model for a fluid-flexible structure system: numerical simulations and experiment validations[END_REF] by a comparison of a time-response for a given initial deformation of the plate. However, in order to obtain a model which provides a good match of the measured frequency response, some adjustments are required by considering the Bode plots. This adaptation is done following a trial-and-error method (first the frequencies are matched and then is the damping). Other methods are possible for flexible structures (see e.g., [START_REF] Roos | Generation of flexible aircraft lft models for robustness analysis[END_REF] and references therein).

This model matching is necessary since some mechanical elements are not well known and have not been taken into account in the modeling. These elements include the circular ring used to attach the tank on the plate, the dynamical behavior of the piezoelectric patches, the non-homogeneity of the plate and the weight of the tank.

As a first step of the model matching, the frequencies of the plate are adjusted. A bigger amount of liquid will be sensed by the plate as an increase of the total mass of the plate and this will lead to (see [START_REF] Blevins | Formulas for natural frequency and mode shape[END_REF]Chapter 11]) a shift towards zero of all the mode frequencies. A second step of the model matching is the adding of a static gain that corresponds to the high frequency modes neglected during the model reduction static correction [START_REF] Rubin | Improved component-mode representation for structural dynamic analysis[END_REF]. This allows to get a more realistic model at low frequencies. The comparison of the bode plots for e = 0.7 and e = 0.9 on Figures 4 and5 shows that the model, for e = 0.7 is quite accurate with respect to the real data while there is some discrepancy in the amplitude of the first sloshing mode for e = 0.9. In the considered figures, the first peak corresponds to the first flexional mode of the plate (0.625 Hz) and the second peak to the first sloshing mode (1.25 Hz for e = 0.7 and 1.415 Hz for e = 0.9) in the tank. The next four peaks are respectively representing: The first torsional mode (the third peak) (6.38 Hz) and the second (8.75 Hz), third (14.45 Hz) and fourth (21.50 Hz) flexional modes of the plate. The second mode of the liquid sloshing (1.99 Hz for e = 0.7 and 2. [START_REF] Géradin | Mechanical Vibrations: theory and application to structural dynamics[END_REF] Hz for e = 0.9) cannot be identified on the Bode diagrams due to its very small amplitude.

Robust H ∞ control problem

In this section a robust controller is computed and some experiments are performed. The controller is calculated using the standard H ∞ problem given in Figure 6. In order to take disturbances into account, the low-pass filter H 1 (s) is included in the design scheme. The residual modes divergence describing the spillover phenomenon is a common problem when working with a truncation of an infinite-dimensional systems. In order to avoid this non desirable effect, a high-pass filter H 2 (s) is added on the controlled output (see [START_REF] Tliba | Active vibration damping of a smart flexible structure using piezoelectric transducers: H∞ design and experimental results[END_REF]). This filter has a transfer function

H 2 (s) = (1 + s 2π25 ) 3 (1 + s 2π160 )
3 that allows to get a 60 dB attenuation above the cut-off frequency of 25 Hz. The cut-off frequency is slightly greater than the frequency of the last considered mode in the controller synthesis.

The H ∞ controller is designed and is first tested through simulations and on the experimental setup afterwards. The simulation model is a system of larger dimension, with the first six modes of the plate and the first three modes of the liquid sloshing, in order to test the existence of the spillover effect. Two different levels of tank filling are considered (e = 0.7 and e = 0.9).

When looking at the synthesis plant model augmented by the filters H 1 and H 2 and keeping in mind the number of modes retained for design, it is clear that a full-order H ∞ controller given by the Matlab c Robust Control Toolbox will have severe restrictions for its numerical implementation on the experimental set-up. More clearly, the Robust Control Toolbox c fails in finding a controller when the number of considered modes is increasing and when the two filters are added. Moreover, it cannot directly tackle the problem of different filling levels in the tip-tank. It is therefore reasonable to resort to algorithms for the synthesis of reduced-order controllers for the simultaneous H ∞ control problem. This is the case of the HIFOO package based on non smooth optimization from [START_REF] Gumussoy | Multiobjective robust control with HIFOO 2.0[END_REF].

Order of the reduced-order HIFOO controllers

At first, one specific controller is considered for each tank filling level. In order to choose the suitable order of these HIFOO controllers, H ∞ controllers of different orders for a fixed tank filling of 0.7 are computed using the Simulink c diagram of Figure 6. The computations show that controllers of order 4 and 1, have almost the same H ∞ norm: 4.28 for a 4 th order and 4.24 for a 1 st order. Consequently, a 4 th order controller and a 1 st order controller for the same tank filling e = 0.7 are tested on the plant. The idea behind this is to see if greater order controllers are really more efficient than a very simple first-order controller. The experimental results are plotted in Figure 7. The Bode plots show the closed-loop attenuation in the case of a 4 th order and a 1 st order controller computed using HIFOO. One can notice a slightly better attenuation for the first sloshing mode in case of the 4 th order controller and a better attenuation (approx. 4 dB) for the first flexion mode in the case of the first-order controller.

It can be seen that the complexity of a 4 th order controller is not justified. Therefore, from now on, only first-order controllers will be computed with HIFOO.

Experimental results are given in Figures 8 and9. The transfer functions of both controllers are given below. Therefore, the transfer function of the first order HIFOO controller for the case e = 0.9 is: K 09 = 0.1112s + 0.5257 s -0.0504 and the one for the case e = 0.7 is:

K 07 =
-1.826s + 7.155 s + 6.52 .

It may be observed that the first peak is well attenuated for the different considered tank fillings. An attenuation of 14 dB is measured when e = 0.9 and of 11.7 dB when e = 0.7.

Concerning the first twisting mode (3 rd peak on the Bode plots) the attenuation is very small for e = 0.7 and quite good for e = 0.9 (1.5 dB). For higher order modes, it is seen that the controller for e = 0.9 is also quite efficient.

Further on, in order to analyze time behavior of the mock-up, the plate is allowed to freely oscillate after release from a deflected position of 10 cm at the plate free end, corresponding to the excitation of the first flexion mode. The tank filling level is chosen to be equal equal to 0.9 and HIFOO controller K 09 is tested for this new configuration. The time response of the closed-loop system, compared to the open-loop one is given in Figure 10, while the input voltage delivered by the first order HIFOO controller K 09 is depicted in Figure 11. We notice, that the voltage remains in the saturation limits ±100V of the amplifier delivering voltage to the piezoelectric actuator.

Simultaneous reduced-order HIFOO controller

In practice, the liquid in the plane tanks is varying during flight. Therefore, one controller must be valid for different fillings.

The first-order controller previously computed for the 90% filled tank is tested on the 70% filled tank. One can notice from Figure 12 that the controller increases the amplitude of the first vibration mode of the plate and does not attenuate the other modes. Therefore a simultaneous first order robust controller is computed for all different levels using again the HIFOO package under Matlab c .

In Figures 13 and14, the Bode plot of the experimental setup in open-loop and the Bode plot of the experimental setup in closed-loop using the simultaneous HIFOO controller are compared It is observed that the first and most important mode in terms of plate displacement is very well attenuated for e = 0.9 (10 dB) and e = 0.7 (5.7 dB). Regarding the twisting mode and higher order modes they are also well attenuated especially for e = 0.9.

Conclusion

After describing the experimental setup used to test the obtained theoretical results, we computed an infinite dimensional model for the plate bending and liquid sloshing using partial derivative equations. A model approximation using Ritz method was then performed and a state -space representation of the system was realized. Before the computation of robust controllers, a model matching is performed, in order to have a better similarity in terms of amplitude and frequency between the model and the plant. Finally, two controllers are calculated for the model considering 5 modes for the plate and 2 modes for liquid sloshing using HIFOO package. The Robust Control Toolbox from Matlab c is not used since it is not well adapted to the system under study. Tests on the experimental setup are achieved. With a simple first-order controller computed using HIFOO, vibrations alleviation is obtained for the first modes on one tested tank filling but performance is not so good when applying this controller for another tank filling demonstrating some lack of robustness. This last issue is greatly improved by designing a simultaneous first-order controller that shows much better performance on both tip-tank fillings. Work is in progress to compute a controller that is more robust to the filling of the tank by using different weightings for the simultaneous HIFOO controller.

Appendix

Tank transformation method

To compute the dimensions of the virtual rectangular tank we propose the following method. The method was chosen after comparing the computed natural frequencies to the experimental ones on one side (the precision of the method) and after analyzing the computation complexity plus the time employed in implementing the method.

This method keeps the length b and the width a of the rectangular tank equal to the length of the cylinder and the width l s of the free surface respectively.

l s = R 2 1 -(2e -1) 2 .
where R is the radius of the internal diameter of the tank and e = liquid height 2R is the tank filling level in the cylindrical tank.

The liquid height is chosen so that the volume of liquid in both tanks is the same. By equating the two volumes, the liquid height in the rectangular container is:

h = R(θ -sin(θ)) 1 -(2e -1) 2 . ( 26
)
with θ = 2 arccos(1 -2e). For more details about this method and for alternative approximation methods see [34, Chapter 2.3.2].

Boundary and initial conditions

Plate The boundary conditions for the three free sides of the plate are given by:

∂ 3 w ∂y 3 = ∂ 3 w ∂z 3 = ∂ 2 w ∂y 2 = ∂ 2 w ∂z 2 = 0 ∀(y, z) ∈ {L} × [0, l], ∂ 3 w ∂y 3 = ∂ 3 w ∂z 3 = ∂ 2 w ∂y 2 = ∂ 2 w ∂z 2 = 0 ∀(y, z) ∈ (0, L] × {0, l}. (27)
and for the clamped side:

w = ∂w ∂y = ∂w ∂z = 0 ∀(y, z) ∈ 0 × [0, l]. (28)
The plate initial conditions are:

w(y, z, 0) = w 0 (y, z) ∀(y, z) ∈ [0, L] × [0, l] ∂w ∂t (y, z, 0) = w 1 (y, z) ∀(y, z) ∈ [0, L] × [0, l] (29)
where w 0 and w 1 stand for the initial deformation and velocity respectively.

Tank with liquid The boundary conditions for the tank with liquid are: ∂φ ∂x = 0, ∂φ ∂y = 0, ∂φ ∂z = 0 on the tank walls ∂δ ∂t = ∂φ ∂z for z = z s [START_REF] Pommier-Budinger | Fractional robust control with iso-damping property[END_REF] which state that the speed of the liquid on the tank walls are zero and the speed of the liquid on the free surface is only the time derivative of the liquid movement along the mean level h. Moreover the initial condition of the liquid is:

φ(x, y, z, 0) = φ 0 (x, y, z) ∀(x, y, z) in the liquid.

(31)

Basis functions for the infinite-dimensional model

For the plate equation ( 1), we consider the homogeneous pde in the absence of external forces:

m s Y I s ∂ 2 w ∂t 2 + ∆ 2 w = 0 ( 32
)
and use the separation of variables method. The displacement w can then be expressed in a Hilbertian basis:

w(y, z, t) = ∞ n=1 η n (y, z)q n (t) = ∞ n=1
Y in (y)Z jn (z)q n (t) [START_REF] Prieur | Stabilization of a 1-D tank containing a fluid modeled by the shallow water equations[END_REF] where Y in (y) and Z jn (z) are the Ritz functions of two beams along the y-axis and z-axis. For a more detailed description of the choice of Ritz functions see [START_REF] Géradin | Mechanical Vibrations: theory and application to structural dynamics[END_REF]Chapter 5.2]. φ(x, y, z, t) and δ(x, y, t) are then calculated by solving (2), (3) with the proper boundary conditions [START_REF] Pommier-Budinger | Fractional robust control with iso-damping property[END_REF] and the initial conditions [START_REF] Pommier-Budinger | Active control of a structure with sloshing phenomena[END_REF]. We first write them in the Hilbertian basis of L 2 ([0, a] × [0, b]) composed of the eigenfunctions of ∆ x,y with Newman homogeneous boundary conditions. Thus, they are written as:

φ(x, y, z, t) = ∞ i=0 ∞ j=0 g ij (t)f ij (z)S ij (x, y) (34)
δ(x, y, t)

= ∞ i=0 ∞ j=0 r ij (t)S ij (x, y) (35)
Because the tank is rectangular, the separation of variables method can be used to write S ij (x, y) as a product of cosine functions depending separately on x and y: S ij (x, y) = cos(π i a x) cos(π i b y)

Using ( 30) and (2) we find:

g ij (t) = dr ij (t) dt f ij (z) = cosh(Υ ij z) Υ ij sinh(Υ ij h)
where Υ ij = π i a + π j b . As the movement of the liquid is along the x-axis, it can be proven that the free surface displacement and the velocity potential are not depending on the variable y and thus [START_REF] Robu | Active vibration control of a fluid/plate system[END_REF] and [START_REF] Robu | A controlled distributed parameter model for a fluid-flexible structure system: numerical simulations and experiment validations[END_REF] can be written as (see [START_REF] Khandelwal | A mechanical model for liquid sloshing in a rectangular container[END_REF]): δ(x, y, t) = δ(x, t) = ∞ i=0 r i (t) cos(Υ i x) [START_REF] Roos | Generation of flexible aircraft lft models for robustness analysis[END_REF] φ(x, y, z, t) = φ(x, z, t)

= ∞ i=0 ṙi (t) cosh(Υ i z) Υ i sinh(k i h) cos(Υ i x) (37)
where Υ i = π i a . Replacing ( 4), ( 36) and ( 37) in (3), we get a differential equation for r i which is to be solved at the free surface. This allows us to compute the r i functions. After tedious calculations, the whole moment exerted by the liquid sloshing is given by: where w(y, z, t) = N k=1 Y i k (y)Z j k (z)q k (t) is the deformation of the plate given by equation ( 33) and that the momentum along both axes is given by (6.4), we obtain from the earlier equation the components b p k of the input matrix B p . where the Heaviside step function is defined as:

m y = -2ρb
H(r) = 0 if r > 0, 1 if r < 0 (40)
and K b is a coefficient depending on the physical parameters of the plate and piezoelectric actuator.

For more details about this computation [34, Chapter 2.2.3.2] can be consulted.

Figure 2 :

 2 Figure 2: Deformation of the rectangular plate (1 st mode)The length of the plate is along the horizontal axis and the width along the vertical one. The plate is clamped on one end and free on the three other sides. The characteristics of the aluminium plate are given in Table1.

Figure 3 :

 3 Figure 3: Experimental setup

Figure 4 :

 4 Figure 4: Frequency matching for the tank filling level e = 0.9 (numerical model -plain line and experimental setup -dotted line)

Figure 5 :

 5 Figure 5: Frequency matching for the tank filling level e = 0.7 (numerical model -plain line and experimental setup -dotted line)

Figure 6 :

 6 Figure 6: Standard H ∞ problem

Figure 7 :

 7 Figure 7: Comparison between a 4 th and a 1 st order HIFOO controllere = 0.7

Figure 8 :

 8 Figure 8: Open-loop and closed-loop transfers using HIFOO controllere = 0.9

Figure 9 :

 9 Figure 9: Open-loop and closed-loop transfers using HIFOO controllere = 0.7

Figure 10 : 9 Figure 11 :

 10911 Figure 10: Open-loop and closed-loop time response using K 09e = 0.9

Figure 12 :

 12 Figure 12: HIFOO controller calculated for the tank 90% filled and tested on the tank 70% filled

Figure 13 :

 13 Figure 13: Open loop and closed loop using simultaneous HIFOO controllere = 0.9

Figure 14 :

 14 Figure 14: Open loop and closed loop using simultaneous HIFOO controllere = 0.7

6 . 4 ∂ 2 m a z ∂z 2

 642 Computation of the control matrixIn order to obtain the input matrix B p , we compute the total bending momentum Γ generated by the piezoelectric patch along both axes y and z: )w(y, z, t)dydz[START_REF] Shaw | Active control of the helicopter rotor for vibration reduction[END_REF]

m a y = m a z =

 z K b V a [H(yy a1 ) -H(yy a2 )][H(zz a1 ) -H(zz a2 )]

Table

Table 1 :

 1 . Plate characteristics

	Plate length	L	1.36 m
	Plate width	l	0.16 m
	Plate thickness	h	0.005 m
	Plate density	ρ 2970 kg m -3
	Plate Young modulus	Y	75 GPa
	Plate Poisson coefficient ν	0.33

Table 3 :

 3 Characteristics of the cylindrical tank

	Actuator length/width/thickness	0.14/0.075/5e -4 m
	Sensor length/width/thickness	0.015/0.025/5e -4 m
	Actuator/Sensor density	7800 kg m -3
	Actuator/Sensor Young modulus	67 GPa
	Actuator piezoelectric coefficient	-210e -12 m V -1
	Sensor piezoelectric coefficient	-9.6 N (Vm) -1
	Actuator/Sensor Poisson coefficient	0.3
	Table 2: Characteristics of the piezoelectric patches
	a longitudinal movement when the plate has a flexion movement and a pitch movement if the plate
	has a torsion movement. It has the dimensions given in Table 3 and it can be filled with water or
	ice up to an arbitrary level. If the tank is filled with ice, it can be easily modeled by a steady mass
	[38]. When the tank is filled with water up to a level close to 0 or close to the cylinder diameter
	(near empty tank empty or near full tank), there is no sloshing behavior, and the modeling process
	is similar to the case of frozen water.
	Tank exterior diameter	0.11 m
	Tank interior diameter	0.105 m
	Tank length	0.5 m
	Tank density	1180 kg m -3
	Tank young modulus	4.5 GPa

Acknowledgements

The authors gratefully acknowledge Michael Overton and the Courant Institute of Mathematical Sciences of NYU, New York City, New York, USA, for hospitality, where this work was done. Support for this work was provided in part by the grant DMS-0714321 from the U.S. National Science Foundation. The authors would like to acknowledge Valérie Budinger and ISAE for the helpful support concerning the experimental setup.


Frédéric Cérou
email: frederic.cerou@irisa.frteddy.furon@inria.fr

Pierre Del Moral
email: pierre.delmoral@inria.fr

Teddy Furon

Arnaud Guyader
email: arnaud.guyader@uhb.fr

SEQUENTIAL MONTE CARLO FOR RARE EVENT ESTIMATION *

Keywords: Primary 65C05; secondary 65C35, 60F05, 62G30 Sequential Monte Carlo Methods, Rare Event Analysis

This paper discusses a novel strategy for simulating rare events and an associated Monte Carlo estimation of tail probabilities. Our method uses a system of interacting particles and exploits a Feynman-Kac representation of that system to analyze their fluctuations. Our precise analysis of the variance of a standard multilevel splitting algorithm reveals an opportunity for improvement. This leads to a novel method that relies on adaptive levels and produces estimates with optimal variance.

The motivation for this theoretical work comes from problems occurring in watermarking and fingerprinting of digital contents, which represents a new field of applications of rare event simulation techniques. Some numerical results show the performance of our technique for these practical applications.

1. Introduction. Monte Carlo approach is a common tool to estimate the expectation of any function of a random object when analytical or numerical methods are not available. However, if the function is non-zero on a set which has a very small probability, then the naive Monte Carlo method will return the estimate zero, unless the sample size is so large that it becomes intractable. Typically we want to estimate precisely and in a reasonable time the small probability, say 10 -9 or below, of an extreme event. A naive Monte Carlo method is impractical as it would require an excessively large sample. To circumvent this difficulty, importance sampling (see e.g. Bucklew [4]) changes the law of the simulated random objects, and reweighs the sample consequently. The difficulty is then to choose the appropriate change of probability in order to achieve a good estimate. This is not always obvious, especially when there is no large deviation result to consider. Importance splitting, or multilevel splitting, is another approach that is well adapted when the random object is a Markovian process. The basic idea is to reinforce trajectories that approach the targeted set by splitting (or branching) them and discarding the others. This very powerful approach in fact dates back to Kahn and Harris [START_REF] Kahn | Estimation of particle transmission by random sampling[END_REF] and Rosenbluth and Rosenbluth [START_REF] Rosenbluth | Monte Carlo Calculation of the Average Extension of Molecular Chains[END_REF]. We refer the reader to the paper by Glasserman et al. [START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF] which contains a precise review on these methods as well as a detailed list of references. Recently, the connection between importance splitting for Markovian processes and particle methods for Feynman-Kac models has lead to some improvements and to a rigorous framework for mathematical analysis (see Del Moral [START_REF] Del Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF]). Unlike previous works concerning rare event estimation and simulation, the present work deals with rare events for a fixed probability distribution. We are simply concerned here with events of the type {X ∈ A} for some random vector X, with p = È(X ∈ A) = È(Φ(X) > L) ≪ 1, where Φ is a mapping from Ê d to Ê, and where there is no dynamical model for X, i.e. X is not a process indexed by the time. In order to use the framework developed for Markov processes (see Cérou et al. [START_REF] Cérou | Genetic Genealogical Models in Rare Event Analysis[END_REF], Cérou and Guyader [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF], Del Moral and Lezaud [START_REF] Del Moral | Branching and interacting particle interpretation of rare event probabilities[END_REF]), we construct a family of Markov transition kernels M k whose invariant measures are the successive laws of X restricted on smaller and smaller sets, the smallest being A. As usual when using a splitting technique in rare event simulation, we decompose the rare event in not so rare nested events, whose product of probabilities equals the probability of the rare event. Del Moral et al. [START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF] and Johansen et al. [START_REF] Johansen | Sequential Monte Carlo samplers for rare events[END_REF] were first to use fixed-levels algorithms for static rare events. These articles were written in a different framework, and thus do not deal with the practical details of our precise setting. In the present article, we detail both a fixed and an adaptive multilevels algorithm, the adaptive one consisting in optimally placing the levels on the fly. Recently and independently from the authors of the present work, Botev and Kroese [START_REF] Botev | An Efficient Algorithm for Rare-event Probability Estimation, Combinatorial Optimization, and Counting[END_REF] proposed the same approach. These authors work on a similar algorithm, including the use of quantiles of the random variable Φ(X) on the swarm of particles in order to estimate the next level. The main difference is their two stage procedure (like in Garvels [START_REF] Garvels | The splitting method in rare event simulation Thesis[END_REF]): they first run the algorithm just to compute the levels, and then they restart from the beginning with these proposed levels. Actually we prove that by computing the levels on the fly (within the same run as the one to compute the rare event probability), we only pay a small bias on the estimate. Note also that [START_REF] Botev | An Efficient Algorithm for Rare-event Probability Estimation, Combinatorial Optimization, and Counting[END_REF] does not address the general construction of the transition kernels M k , since the authors only tackle examples where they can derive a Gibbs sampler at each step. This is mainly possible because their function Φ is linear, which is a severe restriction. Another related approach is the recent work on combinatorial counting of Rubinstein [START_REF] Rubinstein | The Gibbs Cloner for Combinatorial Optimization, Counting and Sampling[END_REF]. This article deals with optimization and counting problems in which X has a uniform distribution over a discrete but very large state space. The author uses what he calls a cloning procedure, where the number of offsprings is fixed (i.e. the same for all the particles in the sample) but adaptive to keep the number roughly constant, while removing redundant particles after the MCMC step. This is a main difference since we use a resampling with replacement procedure. Finally, we would like to mention that these last two papers [3][24] have demonstrated the performance of their algorithms via an extensive simulation study, to which we now lay out the mathematical foundations. The paper is organized as follows. Section 2 describes and analyses the fixedlevels algorithm. Section 3 provides and studies the adaptive levels version, which proves to be optimal in terms of asymptotic variance of the estimator. Section 4 deals with the tuning of the algorithm and especially the choice and the iteration of the transition kernel which is at the core of the method. Section 5 shows the relevance of our algorithm for watermarking and fingerprinting, which constitute a new application area of rare event simulation techniques. Finally, all the proofs are gathered in the appendix.

2. The fixed-levels method.

Assumptions and ingredients.

We assume that X is a random element on Ê d for some d > 0, and denote by µ its probability distribution on the underlying probability space (Ω, F, È). We denote by A the rare set of interest, and we assume that A = {x ∈ Ê d s.t. Φ(x) ≥ L} for some function Φ : Ê d → Ê and some real number L. We also assume that we know how to draw i.i.d. samples from µ. Our algorithm makes use of the following ingredients. An increasing sequence {L 0 , . . . , L n } in Ê, with L 0 = -∞ and L n = L defines a sequence of corresponding sets A k = {x ∈ Ê d , Φ(x) ≥ L k }. These sets are thus nested:

Ê d = A 0 ⊃ A 1 ⊃ • • • ⊃ A n = A. We now need to choose the sequence {L 0 , . . . , L n } in such a way that p k = È(X ∈ A k |X ∈ A k-1 ) is not too small. For indices k > n, we assume that L k = L n and A k = A n . We also need to choose a Markov transition kernel K on Ê d which is µ-symmetric, that is ∀(x, y) ∈ Ê d × Ê d , µ(dx)K(x, dy) = µ(dy)K(y, dx).
As a consequence, K has µ as an invariant measure. As we will see in the sequel, the choice of the L k 's can be made adaptive and is thus not an issue. However, the choice of the kernel K is crucial. Even if any µ-symmetric kernel would eventually do the job, we need to carefully choose it to make the algorithm efficient as discussed in section 4. Consider now a Markov chain (X k ) k≥0 defined by: L(X 0 ) = µ and the inhomogeneous transitions kernels È(

X k ∈ dy|X k-1 = x) = M k (x, dy), with M k (x, dy) = ½ A c k (x) δ x (dy) + ½ A k (x)(K(x, dy)½ A k (y) + K(x, A c k ) δ x (dy))
. Moving a particle according to M k is then twofold: firstly a new transition according to K is proposed, and secondly we accept this transition only if it stays in A k , keeping the old position otherwise. For k ∈ {0, . . . , n}, denote

µ k (dx) = 1 µ(A k ) ½ A k (x) µ(dx)
the normalized restriction of µ on A k . The following stationarity property holds for µ and µ k .

Proposition 1. The measures µ and µ k are both invariant by the transition kernel M k .

From a general point of view, a Feynman-Kac representation for µ k is a formula of the form

µ k (ϕ) = [ϕ(X k ) k-1 m=0 G m (X m )] [ k-1 m=0 G m (X m )]
,

where the potentials G m are positive functions, and (X k ) k≥0 is a non homogeneous Markov chain with transitions M k . If we know how to draw realizations of the Markov chain, then we can compute µ k (ϕ) with a Monte Carlo approach. But naive Monte Carlo is not efficient, because most of the realizations of the chain have small values for the product of the potentials. Anyway, in this form a much nicer Monte Carlo algorithm can be used. It mainly consists in keeping a cloud of particles (ξ j k ), with time 0 ≤ k ≤ n and particle index 1 ≤ j ≤ N . Then for each time step k, discard those with small potential G k , and branch the others, with a rate proportional to G k (ξ j k ). Then apply the Markov transition M k to all the surviving particles, and iterate on the time step. This approach has given birth to a huge amount of literature, and is often referred to as Sequential Importance Sampling (SIS) or Sequential Monte Carlo (SMC). See the monograph by Del Moral [START_REF] Del Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF] for a theoretical overview and Doucet et al. [START_REF]Sequential Monte Carlo Methods in Practice[END_REF] for examples of applications. In our context, the Feynman-Kac representation for µ k has the following form. Proposition 2. For every test function ϕ, for k ∈ {0, . . . , n}, the Feynman-Kac representation is as follows where (X k ) k≥0 is a Markov chain given by the following conditions: X 0 ∼ µ and the inhomogeneous transition kernels (M k ) k≥1 .

µ k (ϕ) = [ϕ(X k ) k-1 m=0 ½ A m+1 (X m )] [ k-1 m=0 ½ A m+1 (X m )]
2.2. The fixed-levels algorithm. Proposition 2 shows that the framework of Feynman-Kac formulae does apply, and thus this grants access to the approximation of the associated measures using an interacting particle method as studied by Del Moral in [START_REF] Del Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF]. Basically, at each iteration k, it consists in selecting the particles according to the potentials, here ½ A k+1 , and then in propagating the particles according to the transitions given by M k+1 .

The approximation of the rare event probability stems from the following obvious property

p = È(X ∈ A n ) = n-1 k=0 È(X ∈ A k+1 |X ∈ A k ) = n-1 k=0 µ k (A k+1 )
and finally

p = n-1 k=0 [½ A k+1 (X k ) k-1 m=0 ½ A m+1 (X m )] [ k-1 m=0 ½ A m+1 (X m )]
,

where the last equality comes from Proposition 2. We approximate at each

stage p k = È(X ∈ A k+1 |X ∈ A k
) by the proportion of the particles already in the next set, and the total probability is estimated as the product of those. This gives the algorithm 1.

Algorithm 1.

Parameters

N the number of particles, the sequence {L 0 , . . . , L n } of levels.

Initialization

Draw an i.i.d. N -sample (ξ j 0 ) 1≤j≤N , of the law µ.

Iterations The last set of particles is a (non independent) sample that provides a good approximation of the law µ n of the rare event.

for k = 0 to n -1 /* level number */ Let I k = {j : ξ j k ∈ A k+1 }. Let pk = |I k | N . for j = 1 to N /*
2.3. Fluctuations analysis. Del Moral has extensively studied in a very general context the asymptotic behavior of the interacting particle model as the number N of particles goes to infinity [START_REF] Del Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF]. For example, it is well known that the estimate p is unbiased. The next proposition presents a precise fluctuation result in our context of rare event analysis. Proposition 3. Let p denote the estimate given by the fixed-levels algorithm, then

√ N p -p p D -----→ N →+∞ N (0, σ 2 ),
with

σ 2 = n-1 k=0 1 -p k p k + n-1 k=0 1 p k È(X n-1 ∈ A n |X k ) È(X n-1 ∈ A n |X k-1 ∈ A k ) -1 2 X k-1 ∈ A k .
The samples are not independent due to the splitting of successful particles. In fact the variance is lower bounded

σ 2 ≥ n-1 k=0 1 -p k p k ,
with equality if and only if for all k = 0, . . . , n -1 and knowing that

X k-1 ∈ A k , one has È(X n-1 ∈ A n |X k ) ⊥ X k .
This means that equality holds if, between step k and step n -1, the algorithm forgets the initial position X k . In order to reach this goal, a possible imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010 route is to begin step k by applying an infinite number of times (and not only one time as is the case in Algorithm 1) the transition kernel M k with stationary distribution µ k = L(X|X ∈ A k ). We will discuss this point in details in section 4.

Anyway, from now on, we assume that at each step k it is possible to draw an i.i.d. sample of the law of X conditionally on the event

{X ∈ A k } = {Φ(X) > L k }.
Then the relative variance of the estimator reduces to:

σ 2 = n-1 k=0 1 -p k p k .
Thus, for a fixed value of p and a fixed number n of levels, this asymptotic variance would be minimal if p k ≡ p 0 for all k. This is indeed a simple constrained optimization problem:

arg min p 0 ,...,p n-1 n-1 k=0 1 -p k p k s.t. n-1 k=0 p k = p.
In this case, the minimal asymptotic variance is simply n 1-p 0 p 0 , with p 0 = p 1 n . This optimal situation corresponds to the case where the levels are evenly spaced in terms of probability of success: as far as multilevel splitting for Markov processes is concerned, this point was also mentioned in Glasserman et al. [START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF], Lagnoux [START_REF] Lagnoux | Rare event simulation[END_REF] and Cérou et al. [START_REF] Cérou | Genetic Genealogical Models in Rare Event Analysis[END_REF]. The following section addresses this crucial issue for the adaptive version of the algorithm. Before this, this section ends with two remarks.

Remarks:

1. If one's particular interest is the variance of p rather than a convergence in distribution like the CLT-type result of Proposition 3, then we can turn to the recent non asymptotic results obtained in Cérou et al. [5, corollary 5.2]. Under some regularity conditions (mainly about the mixing property of the kernel K), there exist positive constants α k , for 0

≤ k ≤ n -1, such that for all N ≥ N 0 = n-1 k=0 α k p k , p -p p 2 ≤ 4 N 0 N .
If we assume an i.i.d. sample, then all the α k 's are all equal to 1, and

N 0 = n-1 k=0 1 p k .
imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010 2. Finally, there is a maybe small, but non-zero, probability that the particle system dies at some stage. This may typically happen when two consecutive levels are too far apart, or when the number of particles is too small. A first solution to this problem is given in Le Gland and Oudjane [START_REF] Le Gland | A sequential algorithm that keeps the particle system alive[END_REF]. The idea is to go on sampling new particles until a given number of them have reached the next level. The price to pay is a possibly very long computation time. A second solution is proposed in the next section.

3. The adaptive method.

3.1. The algorithm. As we may not have a great insight about the law µ and/or the mapping Φ, typically when Φ is a 'black box', the choice of the levels L 1 , . . . , L n-1 might prove to be quite problematic. We propose from now on to adaptively choose the levels, ensuring not only that the particle system never dies but also that the asymptotic variance of the estimate p is minimized.

The method is indeed simple to implement. We choose a prescribed success rate p 0 between two consecutive levels (e.g. p 0 = 0.75). After application of the kernel M k , the algorithm sorts the particles ξ j k according to their scores Φ(ξ j k ). Then it sets the next level to the (1p 0 ) empirical quantile Lk+1 , which means that a proportion p 0 of the particles scores are above it. Starting from this sample of p 0 N particles which are independently and identically distributed according to the law L(X|Φ(X) > Lk+1 ), an i.i.d. sample of size N is drawn with the same distribution, and the rest of the algorithm is unchanged. The algorithm then stops when some Ln 0 +1 ≥ L, and the probability is estimated by p = r0 p n0 0 , where r0 denotes the number of particles in the last iteration being above level L. The number n0 of steps is random, but if N is large enough, then Lemma 1 in Appendix proves that, outside an event of exponentially small probability, n0 is actually fixed by the ratio of the logarithms

(3.1) n 0 = log È(X ∈ A) log p 0 = log p log p 0 .
As mentioned above, this variant enforces evenly spaced levels in terms of probability of success, and therefore a minimal asymptotic variance for the estimate p of p. The pseudo-code for the adaptive algorithm is given in Algorithm 2 below.

-p 0 ) quantile of Φ(ξ j k ), j = 1, . . . , N . k = k + 1; endwhile Let N L the number of particles ξ j k-1 , j = 1, . . . , N , such that Φ(ξ j k-1 ) ≥ L.

Output

Estimate the probability of the rare event by p = N L N p k-1 0 . The last set of particles is a (non independent) sample that provides a good approximation of the law µ n of the rare event.

Remarks:

1. In this algorithm, the step drawing an N -sample starting from a p 0 Nsample is of course the trickiest one. The analytical study of the next subsection assumes it can be done perfectly. In section 4, we propose a way to implement it in practice, at least approximately, and section 5.1 shows its practical efficiency. 2. The costs of adaptive levels is a higher complexity by a factor log N (due to the quick sort), and a slight loss of accuracy due to a bias. Yet, Proposition 4 proves that this bias becomes negligible compared to the standard deviation as N increases. Experimental results of section 5.1 illustrate this. 3. Estimation of quantiles: Some applications require the estimates of quantiles of the random variable Φ(X). This can be done at no additional cost within the previous algorithm. For α ∈ (0, 1), define the α-quantile by q α = sup{x : P (Φ(X) ≤ x) ≤ α}. When the algorithm is in step k, with a set of particles {ξ j k , j = 1, . . . , N }, such that

p k+1 0 ≤ 1 -α < p k 0 , then let r = 1 -(1 -α)p -k .
An estimate of the quantile q α is then given by the r quantile of the sample {Φ(ξ j k ), j = 1, . . . , N }.

Bias and variance.

The assumption of a continuous cumulative distribution function (cdf) F of Φ(X) is now required to derive the properties of the adaptive algorithm. Let us write the rare event probability as p = r 0 p n 0 0 , with n 0 = log p log p 0 and r 0 = p p -n 0 0 , so that r 0 ∈ (p 0 , 1]. In the same way we write p = r0 p n0 0 , with n0 the number of steps before the algorithm stops. A first theorem shows a CLTtype convergence.

Theorem 3.1. If F is continuous, then we have √ N (p -p) D -----→ N →+∞ N (0, σ 2 ),
where

σ 2 = p 2 n 0 1 -p 0 p 0 + 1 -r 0 r 0 .
Unlike the fixed-levels version of the algorithm, the adaptive version is biased. Nevertheless, the next result shows that the bias is of order 1/N , and is thus negligible compared to the standard deviation given in Theorem 3.

N [p] -p p -----→ N →+∞ n 0 1 -p 0 p 0 .
Thus the bias is positive and of order 1 N when N goes to infinity:

[p] -p p ∼ 1 N n 0 (1 -p 0 ) p 0 .
Putting all things together, we can write the following expansion:

p = p 1 + 1 √ N n 0 1 -p 0 p 0 + 1 -r 0 r 0 Z + 1 N n 0 1 -p 0 p 0 + o È 1 N ,
where Z is a standard Gaussian variable.

Remarks:

1. If one does not want any bias, then the solution is to make a first run to compute the levels, and a second run to actually compute the probability. But it is clear from the above formula that, for the same computational cost, it is better in term of relative error to use directly our algorithm 2 and compute the levels on the fly. 2. In this regard the remark 3.2 page 489 of [START_REF] Botev | An Efficient Algorithm for Rare-event Probability Estimation, Combinatorial Optimization, and Counting[END_REF] might be misleading: if the levels are chosen from start, or using a preliminary run, then the resulting probability estimate is unbiased, even if the level crossing probabilities are indeed dependent (see Cérou et al. [START_REF] Cérou | Genetic Genealogical Models in Rare Event Analysis[END_REF]). 3. It is worth mentioning that the bias is always positive, giving a slightly overvalued estimate. As rare event analysis usually deals with catastrophic events, it is not a bad thing that the real value be a bit lower than the provided estimate.

4. Tuning the algorithm.

Choice of the kernel K.

There is no completely general method for finding the best transition kernel K because it depends on the application. But in the very classical case of a Gibbs measure given by a bounded potential, we can use the Metropolis algorithm, as first proposed by Metropolis et al. [START_REF] Metropolis | Equation of State Calculations by Fast Computing Machines[END_REF], or the more general version later proposed by Hastings [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF].

Less dependent sample.

As mentioned in section 2.3, for the fixedlevels version of the algorithm we always have

σ 2 ≥ n-1 k=0 1 -p k p k .
The equality holds if and only if for all k, knowing that

X k-1 ∈ A k , one has È(X n-1 ∈ A n |X k ) ⊥ X k .
To reach this goal, a simple idea is to iterate the transition kernel M k at each step as it provides more independence among particles.

The application of the kernel over the particles altogether can be seen as a kernel (M k ) ⊗N , and the joint law of an i.i.d. sample (µ k ) ⊗N is an invariant measure for (M k ) ⊗N according to Proposition 1. Then, from Meyn and Tweedie [22, Prop. 13.3.2], the total variation distance between the law of the sample and (µ k ) ⊗N is non increasing as we iterate the kernel. So these iterations can only improve the independence of the sample, and thus decrease the variance.

We now focus on Metropolis-Hastings type kernels. We define these kernels as Tierney [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF]: Denote λ the Lebesgue measure Ê d and assume that µ has a density with respect to λ: µ(dy) = f µ (y)λ(dy). We also need an instrumental kernel Q with the same density assumption: Q(x, dy) = q(x, y)λ(dy). Now define the Metropolis-Hastings ratio α by

α(x, y) =          min f µ (y)q(y, x) f µ (x)q(x, y) , 1 if f µ (x)q(x, y) > 0 1 if f µ (x)q(x, y) = 0
Let us now denote by ℓ the off-diagonal density

ℓ(x, y) =      q(x, y)α(x, y) if x = y 0 if x = y.
and r(x) = 1 -Ê d ℓ(x, y)λ(dy). We say that K is a Metropolis kernel if it is of the form K(x, dy) = ℓ(x, y)λ(dy) + r(x)δ x (dy).

Note that in the framework of section 2.1, M k is a Metropolis-Hastings kernel as soon as K is.

We further assume that K is µ-irreducible, and that µ({x ∈ Ê d , r(x) > 0}) > 0. This is usually the case, especially whenever ∀(x, y) ∈ Ê d ×Ê d , q(x, y) > 0.

Remark that if K is µ-irreducible, then M k is obviously µ k -irreducible1 . Then, using Tierney [26, Corollary 2], they are also Harris recurrent, and by Theorem 13.3.3 in Meyn and Tweedie [START_REF] Meyn | Markov chains and stochastic stability[END_REF], any initial distribution ν such that2 ν(A k ) = 1 yields

ν(dx)(M k ) m (x, .) -µ k -----→ m→+∞ 0,
where the norm is in total variation. Then, for any initial cloud of particles Ξ = (ξ 1 , . . . , ξ N ) and any test function (ϕ 1 , . . . , ϕ N ),

δ Ξ ((M k ) ⊗N ) m   N j=1 ϕ j   - N j=1 µ k (ϕ j ) = N j=1 (M k ) m (ϕ j )(ξ j ) - N j=1 µ k (ϕ j ) -----→ m→+∞ 0.
For all test functions ϕ on (Ê d ) N , a standard density argument gives

δ Ξ ((M k ) ⊗N ) m (ϕ) -µ ⊗N k (ϕ) -----→ m→+∞ 0.
This means that the more the kernel is iterated, the closer (in distribution) we get to an independent sample. Thus, at each step of Algorithm 2 (adaptive levels), we can think of iterating the kernel a fixed number of times (for example 20 times in the simulations of section 5.1).

4.3.

Mixing property of the kernel K. We have written the algorithms using a unique kernel K for all the iterations. Usually it is quite easy to construct not only one, but a family of kernels that are all µ-symmetric, but with different mixing properties. This is useful when applying K to the current particles, most of the transitions are refused (their scores are below the current threshold). In this case, we propose a change to another kernel which is less mixing, i.e. in some way with "smaller steps", and thus with a lower probability of going below the current level L k . On the other hand, when almost all the transitions are accepted, it means that the kernel is poorly mixing, and that we could decrease the variance by choosing a kernel K that is more mixing, i.e. with "larger steps". For example, this is tuned by the parameter α in section 5.1.

Applications.

Our motivation comes from problems occurring in the protection of digital contents. Here the term watermarking refers to a set of techniques for embedding/hiding information in a digital file (typically audio or video), such that the change is not perceptible, and very hard to remove. See the web site of the Copy Protection Working Group [8] for details. In order to be used in an application, a watermarking technique must be reliable. Here are two application scenarii where a wrong estimation of the probability of error could lead to a disaster.

Copy protection.. Assume commercial contents are encrypted and watermarked and that future consumer electronics storage devices have a watermark detector. These devices refuse to record a watermarked content since it is copyrighted material. The probability of false alarm is the probability that the detector considers an original piece of content (which has not been watermarked) as protected. The movie that a user shot during his holidays could be rejected by his storage device. This absolutely non user-friendly behavior really scares consumer electronics manufacturers. In the past, the Copy Protection Working Group of the DVD forum evaluated that at most one false alarm should happen in 400 hours of video [8]. As the detection rate was one decision per ten seconds, this implies a probability of false alarm in the order of 10 -5 . An accurate experimental assessment of such a low probability of false alarm would demand to feed a real-time watermarking detector with non-watermarked content during 40,000 hours, i.e. more than 4 years! Proposals in response of the CPTWG's call were, at that time, never able to guarantee this level of reliability.

Fingerprinting.. In this application, users' identifiers are embedded in a purchased content. When this content is found in an illegal place (e.g. a P2P network), the copyright holders decode the hidden message, find an identifier, and thus they can trace the traitor, i.e. the customer who has illegally broadcast his copy. However, the task is not that simple because dishonest users might collude. For security reason, anti-collusion codes have to be employed. Yet, these solutions (also called weak traceability codes, see Barg et al. [START_REF] Barg | Digital fingerprinting codes: problem statements, constructions, identification of traitors[END_REF]) have a non-zero probability of error (defined as the probability of accusing an innocent). This probability should be, of course, extremely low, but it is also a very sensitive parameter: anti-collusion codes get longer (in terms of the number of bits to be hidden in content) as the probability of error decreases. Fingerprint designers have to strike a tradeoff, which is hard to conceive when only rough estimation of the probability of error is known. The major issue for fingerprinting algorithms is the fact that embedding large sequences implies also assessing reliability on a huge imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010 amount of data, which may be practically unachievable without using rare event analysis. 5.1. Zero-bit watermarking. In this example of zero-bit watermarking, X is a Gaussian vector in Ê d , with zero mean and identity covariance matrix,

Φ(X) = | X,u | X
where u is a fixed normalized vector (see Merhav and Sabbag [START_REF] Merhav | Optimal watermarking embedding and detection strategies under limited detection resources[END_REF]). Then the region A = {x ∈ Ê d s.t. Φ(x) ≥ L} is a cone as shown on figure 1. For a Gaussian distribution, the obvious choice for the kernel is the following: if we start from any point x, then the new position is given by

x ′ = x + αW √ 1 + α 2 ,
where W is a N (0, I d ) Ê d valued random vector and α a positive number.

This simple setup allows us to compare our estimates of the rare event probability with the result of a numerical integration. For example, in our simulations d = 20 and L = 0.95, so that p = È(Φ(X) ≥ L) ≈ 4.70 • 10 -11 .

For such a low probability, it is of course out of question to run a classical Monte Carlo algorithm. Our algorithm with adaptive levels was run with 20 iterations of the kernel M k at each step. The choice of the mixing parameter α = 0.3 has experimentally proved to be a good trade-off, see discussion in section 4.3 above for details. The proportion of particles surviving from one step to another has been fixed to p 0 = 0.75, so that n 0 = log p log p 0 = 82 and r 0 = p p -n 0 0 ≈ 0.83.

For several numbers of particles, ranging from N = 100 to N = 5, 000, we have run the algorithm 100 times in order to estimate the variance. shows in log-log plots the convergence of the normalized standard deviation to minimum achievable, which is that of i.i.d. samples at each stage, that is (see Theorem 3.1):

Var p -p p ∼ 1 √ N n 0 1 -p 0 p 0 + 1 -r 0 r 0 .
Clearly this indicates that even if we use the mixing kernel a finite number of times (here only 20 times), the empirical variance is on the same level as it would be for an independent sample (that is in the limit of infinite applications of the kernel), and therefore our theoretical results in this setting give a good picture of the accuracy of the actual algorithm. Anyway, from a practical point of view, one would like to obtain an estimation of p by running only one time the algorithm. In this aim, Theorem 3.1 and Proposition 4 allow us to construct confidence intervals. Indeed, we have

p -p p ≈ N n 0 (1 -p 0 )/p 0 N , n 0 (1 -p 0 )/p 0 + (1 -r 0 )/r 0 N ,
so that an approximate 95% confidence interval for p is given by I = [p -, p+ ], where

p+ -= p   1 - n0 (1 -p 0 )/p 0 N ± 2 n0 (1 -p 0 )/p 0 + (1 -r0 )/r 0 N   .
imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

This is illustrated on figure 3, where 100 such confidence intervals have been drawn for N = 500 particles. In this example, 2 of them do not contain the true value p = 4.704 • 10 -11 .

Tardos probabilistic codes.

We are interested here in embedding an identifier in each copy of a purchased content. Then a copy, which is the result of a collusion, is found on the web, and we want to decide whether or not it can be originated from a certain user. The rare event will be to consider an innocent as guilty.

The embedded message, called a fingerprint, consists of a sequence of bits X = (X 1 , . . . , X m ), where each X i is independent from the others, and drawn from a Bernoulli's B(p i ). The p i 's are themselves i.i.d. random variables, drawn from a given distribution with density f on [0, 1]. Then we find a copy with fingerprint y = (y 1 , . . . , y m ) ∈ {0, 1} m . We conclude that a user is guilty if the score

Φ(X) = m i=1 y i g i (X i )
is larger than some value L, for some given functions g i 's. This approach was proposed by Tardos in [START_REF] Tardos | Optimal probabilistic fingerprint codes[END_REF], where he derives optimal choices for f and the g i 's.

To apply our algorithm, we need to choose the kernel K. As the X i 's are independent, we randomly choose r indices {j 1 , . . . , j r } ∈ {1, . . . , m}, with r being a fixed parameter. Then for each j ℓ , we draw a new X ′ j ℓ independently from the Bernoulli distribution B(p j ℓ ). Using our adaptive algorithm, we made some numerical experiments on such codes. More precisely, we can easily estimate the probability of false detection (false positive) for some code length m, and collusion size c. The collusion strategy is to randomly pick up the symbols of pirated copy among the c colluders' sequences. We can also estimate the probability of not accusing someone guilty (false negative). The results for m = 200, and c = 2, 3, 4 are shown on figure 4. From these curves, one can then decide how to set the threshold L to minimize the total error.

APPENDIX A: PROOFS

A.1. Proof of Proposition 1. We first prove that µ k is invariant by the transition kernel M k . We mainly use here the fact that the transition kernel imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010 K is µ-symmetric:

x µ k (dx)M k (x, dy) = x µ k (dx)(K(x, dy)½ A k (y) + K(x, A c k ) δ x (dy)) = x z µ k (dx)K(x, dz)(½ A k (z)δ z (dy) + ½ A c k (z)δ x (dy)) = x 1 µ(A k ) ½ A k (x)µ(dx)K(x, dy)½ A k (y) + z 1 µ(A k ) ½ A k (y)µ(dy)K(y, dz)½ A c k (z) = x 1 µ(A k ) ½ A k (x)µ(dx)K(x, dy)½ A k (y) + z 1 µ(A k ) ½ A k (y)µ(dz)K(z, dy)½ A c k (z) = µ k (dy).
Now, in order to prove that µ is also invariant by the transition kernel M k , we first remark that

µ(dx) = ½ A c k (x)µ(dx) + ½ A k (x)µ(dx) = ½ A c k (x)µ(dx) + µ(A k )µ k (dx). Since we have x ½ A c k (x)µ(dx)M k (x, dy) = x ½ A c k (x)µ(dx) δ x (dy)) = ½ A k (y)µ(dy),
and by the proof above

x µ(A k )µ k (dx)M k (x, dy) = µ(A k )µ k (dy),
we can conclude that

x µ(dx)M k (x, dy) = ½ A k (y)µ(dy) + µ(A k )µ k (dy) = µ(dy),
which ends the proof of Proposition 1.

A.2. Proof of Proposition 2. We use induction to show that

ϕ(X k ) k-1 m=0 ½ A m+1 (X m ) = µ(A k ) µ k (ϕ).
The case k = 0 is obvious. Then assume the property be true for k. We write, using the Markov property and Proposition 1,

ϕ(X k+1 ) k m=0 ½ A m+1 (X m ) = [ϕ(X k+1 )|X 0 , . . . , X k ]½ A k+1 (X k ) k-1 m=0 ½ A m+1 (X m ) = µ k (M k+1 (ϕ)½ A k+1 ) × µ(A k ) = (½ A k µ)(M k+1 (ϕ)½ A k+1 ) = (½ A k ½ A k+1 µ)(M k+1 (ϕ)) = µ k+1 (M k+1 (ϕ)) × µ(A k+1 ) = µ k+1 (ϕ) × µ(A k+1 ).
Then taking the case ϕ = ½ we have

k m=0 ½ A m+1 (X m ) = µ(A k+1 ),
which concludes the proof.

A.3. Proof of Proposition 3. Adopting the notations of Proposition 9.4.1 page 301 of [START_REF] Del Moral | Feynman-Kac formulae, Genealogical and interacting particle systems with applications[END_REF], the application of the first formula of page 304 leads to

Var(p) = [W γ n-1 (1) 2 ] = γ 2 n-1 n-1 k=0 η k   Q k,n-1 (1) η k Q k,n-1 (1) -1 2   ,
In our context this can be rewritten as follows

(A.1) Var(p) p 2 = n-1 k=0 È(X n-1 ∈ A n |X k ) È(X n-1 ∈ A n |X k-1 ∈ A k ) -1 2 X k-1 ∈ A k . imsart-
p 2 = n-1 k=0 È(X n-1 ∈ A n |X k ) 2 X k-1 ∈ A k È(X n-1 ∈ A n ) 2 -1 . Now one can readily check that È(X n-1 ∈ A n |X k ) 2 X k-1 ∈ A k = È(X n-1 ∈ A n |X k ) 2 ½ A k+1 (X k ) X k-1 ∈ A k ,
and this is equivalent to

È(X n-1 ∈ A n |X k ) 2 X k-1 ∈ A k = È(X n-1 ∈ A n |X k ) 2 ½ A k+1 (X k )½ A k (X k-1 ) È(X k-1 ∈ A k ) . Since {X k ∈ A k+1 } implies {X k-1 ∈ A k }, this last expression can be sim- plified È(X n-1 ∈ A n |X k ) 2 X k-1 ∈ A k = È(X n-1 ∈ A n |X k ) 2 ½ A k+1 (X k ) È(X k-1 ∈ A k )
, and rewritten as

È(X n-1 ∈ A n |X k ) 2 X k-1 ∈ A k = È(X n-1 ∈ A n |X k ) 2 X k ∈ A k+1 × È(X k ∈ A k+1 ) È(X k-1 ∈ A k ) . Since p k = È(X k ∈ A k+1 |X k-1 ∈ A k ) and {X k ∈ A k+1 } ⇒ {X k-1 ∈ A k }, we deduce È(X n-1 ∈ A n |X k-1 ∈ A k ) = p k È(X n-1 ∈ A n |X k ∈ A k+1 ), so that È(X n-1 ∈ A n |X k ) 2 X k-1 ∈ A k È(X n-1 ∈ A n |X k-1 ∈ A k ) 2 = È(X n-1 ∈ A n |X k ) 2 X k ∈ A k+1 p k È(X n-1 ∈ A n |X k ∈ A k+1 ) 2 .
Now it remains to notice that

È(X n-1 ∈ A n |X k ) 2 X k ∈ A k+1 È(X n-1 ∈ A n |X k ∈ A k+1 ) 2 = 1 + È(X n-1 ∈ A n |X k ) È(X n-1 ∈ A n |X k ∈ A k+1 )
so that coming back to (A.1) gives the desired result Var(p)

p 2 = n-1 k=0 1 -p k p k + n-1 k=0 1 p k È(X n-1 ∈ A n |X k ) È(X n-1 ∈ A n |X k-1 ∈ A k ) -1 2 X k-1 ∈ A k .
A.4. Proof of Theorem 3.1. In order to simplify the writings, we will suppose that p 0 N is an integer. Then, for all real numbers L and L ′ such that L < L ′ , let us denote

F (L, L ′ ) = È(Φ(X) ≤ L ′ |Φ(X) > L) = F (L ′ ) -F (L) 1 -F (L)
,

with the convention that F (L, L ′ ) = 0 if F (L) = 1.
We first notice the following crucial point: given Lk , the random vectors ξ j k , for j ∈ {1, . . . , N }, are i.i.d., and thus so are the random variables Φ(ξ j k ). Since F is continuous, given Lk , the random variable F ( Lk , Lk+1 ) has the same distribution as the random variable U ((1-p 0 )N ) , where (U i ) 1≤i≤N is a triangular array of identically distributed random variables with uniform law on [0, 1] and row-wise independent, and for all N ≥ 1

U (1) ≤ U (2) ≤ . . . ≤ U (N ) .
Let us denote G N the empirical cdf of (U i ) 1≤i≤N and G(x) = x the cdf of the uniform law on [0, 1]. Then from the basic identities G

N -G ∞ = G -1 N -G ∞ and U ((1-p 0 )N ) = G -1 N (1 -p 0 ), we can deduce that |U ((1-p 0 )N ) -(1 -p 0 )| ≤ G N -G ∞ .
Using Dvoretzky-Kiefer-Wolfowitz (DKW) inequality (see for example [START_REF] Van Der Vaart | Asymptotic statistics. Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]), we have then for all ε > 0

È( G N -G ∞ > ε) ≤ 2 exp(-2N ε 2 ), hence È(|F( Lk , Lk+1 ) -(1 -p 0 )| > ε) = È(|U ((1-p 0 )N ) -(1 -p 0 )| > ε) ≤ 2 exp(-2N ε 2 )
Using Borel-Cantelli lemma, we conclude that

F ( Lk , Lk+1 ) a.s.
----→ From the theory of order statistics (see for example [START_REF] Arnold | A first course in order statistics[END_REF], Theorem 8.5.1, p. 223), we also deduce the convergence in distribution

√ N (1 -F ( Lk , Lk+1 ) -p 0 ) D ----→ N →∞ N (0, p 0 (1 -p 0 )). (A.3)
To prove the result of Theorem 3.1, we proceed by induction, assuming that:

√ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 D ----→ N →∞ N (0, σ 2 k ).
For the next step, we use the decomposition

√ N k+1 m=1 (1 -F ( Lm-1 , Lm )) -p k+1 0 = √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 (1 -F ( Lk , Lk+1 ) -p 0 ) (A.4) +p 0 √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 + p k 0 √ N (1 -F ( Lk , Lk+1 ) -p 0 ).
The almost sure convergence of equation (A.2) and the induction hypothesis ensure that the first term converges in probability to 0 when N goes to infinity. To prove the convergence in distribution of the other terms of (A.4), let us introduce the characteristic function

φ N (t) = exp it p 0 √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 +p k 0 √ N (1 -F ( Lk , Lk+1 ) -p 0 .
Conditioning with respect to L1 , . . . , Lk gives

φ N (t) = exp it p 0 √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 × exp it p k 0 √ N (1 -F ( Lk , Lk+1 ) -p 0 | L1 , . . . , Lk .
Thanks to the strong Markov property of the Lk 's, this can be reduced to

φ N (t) = exp it +p 0 √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 × exp it p k 0 √ N (1 -F ( Lk , Lk+1 ) -p 0 | Lk ,
and we can remark that

exp it p k 0 √ N (1 -F ( Lk , Lk+1 ) -p 0 | Lk = exp it p k 0 √ N (1 -U ((1-p 0 )N ) -p 0 ,
where U ((1-p 0 )N ) is independent of L1 , . . . , Lk . This leads to

φ N (t) = exp it +p 0 √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 × exp it p k 0 √ N (1 -U ((1-p 0 )N ) -p 0 .
Thanks to the induction hypothesis and to equation (A.3), it comes

p 0 √ N k m=1 (1 -F ( Lm-1 , Lm )) -p k 0 +p k 0 √ N 1 -F ( Lk , Lk+1 ) -p 0 D ----→ N →∞ N 0, p 2 0 σ 2 k + p 2k+1 0 (1 -p 0 ) .
Putting all pieces together gives finally

√ N k+1 m=1 (1 -F ( Lm-1 , Lm )) -p k+1 0 D ----→ N →∞ N (0, σ 2 k+1 ) (A.5) with σ 2 k+1 = p 2 0 σ 2 k + p 2k+1 0 (1 -p 0 ). From this recursion we deduce that for all k ≥ 1 σ 2 k = kp 2k-1 0 (1 -p 0 ).
It remains to deal with the last step. For the sake of simplicity, we suppose that log p/ log p 0 is not an integer. Let us first consider an alternative algorithm defined as follows: we run algorithm 2 with the deterministic number of steps n 0 and denote pd = rd p n 0 0 the corresponding estimator, where

rd = 1 N N j=1 ½ {Φ(ξ j n 0 )≥L} . (A.6)
In this expression, knowing Ln 0 , the random variables (½ {Φ(ξ j n 0 )≥L} ) 1≤j≤N are i. As a consequence n0 a.s.

----→ N →∞ n 0 .

Proof. Let us denote B = {n 0 = n 0 } the event for which the algorithm stops after the correct number of steps. The following equalities are straightforward:

B = Ln 0 ≤ L < Ln 0 +1 = 1 -F ( Ln 0 +1 ) < 1 -F (L) ≤ 1 -F ( Ln 0 ) = n 0 +1 m=1 (1 -F ( Lm-1 , Lm )) < p ≤ n 0 m=1
(1 -F ( Lm-1 , Lm )) .

For all m = 1, . . . , n 0 + 1, if we denote

B m = p

,

  imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

N

  the number of particles, the number N 0 < N of succeeding particles, and let p 0 = N 0 /N . Initialization Draw an i.i.d. N -sample (ξ j 0 ) 1≤j≤N of the law µ. Compute L1 , the (1p 0 ) quantile of Φ(ξ j 0 ), j = 1, . . . , N . k = 1; Iterations while Lk < L do Starting from an i.i.d. p 0 N -sample with law L(X|Φ(X) > Lk ), draw an i.i.d. N -sample (ξ j k ) 1≤j≤N with the same law. Compute Lk+1 , the (1

Fig 1 .

 1 Fig 1. Detection region for zero-bit watermarking.

Figure

2

 2

Fig 2 .

 2 Fig 2. Theoretical and empirical relative standard deviations with 100 simulations.

11 Fig 3 .

 113 Fig 3. 95% confidence intervals for p = 4.704 • 10 -11 with 100 simulations and N = 500 particles.

Fig 4 .

 4 Fig 4. Mappings of the false positive probability (blue) and false negative probability (red) against the threshold. m = 200, c ∈ {2, 3, 4}. The score of a particle is the mean of the c colluders scores.

1 -r 0 p n 0 0 1 -√ N 1 -ÈLemma 1 . 1 n 0 , p 1 n

 1011111 i.d. Bernoulli trials with parameter r = È ½ {Φ(ξ j n 0 )≥L} = 1 Ln 0 = 1-F ( Ln 0 , L) = p F ( Ln 0 ) = F ( Ln 0) . imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010Then we can write√ N (r 0r) = √ N 1 -F ( Ln 0 )p n 0 0 r 0 1 -F ( Ln 0 ). Now we use the almost sure convergence1 -F ( Ln 0 ) in distribution from equation (A.5) F ( Ln 0 )p n 0Coming back to the "true" estimator p, we have√ N (pp) = √ N (pp)½ n0 =n 0 + √ N (pp)½ n0 =n 0 ,but one can readily see that p½ n0 =n 0 = pd ½ n0 =n 0 almost surely, so that√ N (pp) = √ N (p dp) + √ N (ppd )½ n0 =n 0 ,and the proof of Theorem 3.1 will be complete if we show that√ N (ppd )½ n0 =n 0 , let us first remark that for all ε > 0 È(| √ N (pp)½ n0 =n 0 | > ε) ≤ È(n 0 = n 0 ),then the following lemma allows us to conclude.imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010 Denoting c = min p 0p 0 +1p 0 , we have È(n 0 = n 0 ) ≤ 2(n 0 + 1)e -2N c 2 .

1 n 0 1 n 1 -( 1 - 1 n 0 , p 1 nA. 5 .p n 0 0 1 - 1 =F- 1 k=0 1 -

 1011111501111 p 0 < 1p 0 -F ( Lm-1 , Lm ) < p 0 +1p 0 , we have È(B) ≥ È(B 1 ∩ • • • ∩ B n 0 +1 ) ≥ È(B m )). Denoting c = min p 0p 0 +1p 0 , the DKW inequality implies 1 -È(B m ) ≤ È(|1p 0 -F ( Lm-1 , Lm )| > c) ≤ 2e -2N c 2 ,so that the result of Lemma 1 is provedÈ(B) = È(n 0 = n 0 ) ≥ 1 -2(n 0 + 1)e -2N c 2 .imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010 Proof of Proposition 4. As for the analysis of the standard deviation, we first notice thatN (pp) = N (p dp) + N (ppd )½ n0 =n 0 , (A.7)and applying Lemma 1 yieldsN [|ppd |½ n0 =n 0 ] ≤ N È(n 0 = n 0 ) ----→ N →∞ 0,so that only the first term of the right-hand-side of equation (A.7) is worth considering for the convergence of N [pp]. Recall that the estimate is then pd = rd p n 0 0 , where rd is defined as in equation (A.6), so that[r d ] = rd Ln 0 = [r] = p 1 -F ( Ln 0 ) .Then the normalized bias is[p]p p = F ( Ln 0 ) -( Ln 0 ) -F (L n 0 ) 1 -F ( Ln 0 ) = W a -W with W = F ( Ln 0 )-F (L n 0 ) = F ( Ln 0 )-(1-p n 0 0 ),and a = 1-F (L n 0 ) = p n 0 0 . If we remark that 1 -F ( Ln 0 ) = n 0 F ( Lk , Lk+1 ) with the convention L0 = -∞, then the result of equation (A.2) implies W a = F ( Ln 0 ) -(1p n 0 0 ) expansion (1x) -1 = 1 + x + o(x) leads to [p]use the following lemma. imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

  aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

	This leads to
	Var(p)

imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

As discussed in Tierney [26, page 1705], K and M k are also aperiodic.

This is only to avoid uninteresting pathological cases. imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

X k ∈ A k+1 , imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

* This work was partially supported by the French Agence Nationale de la Recherche (ANR), project Nebbiano, number ANR-06-SETI-009.

Lemma 2.

[W ] = [F ( Ln 0 ) -(1p n 0 0 )] = 0 and

The proof of this lemma is left to the next subsection. Coming back to (A.8), we have obtained

which ends the proof of Proposition 4.

A.6. Proof of Lemma 2. First of all, some notation. Let (U i ) 1≤i≤N be an i.i.d. family of random variables uniformly distributed on (0, 1). We denote by

For simplicity, we will assume that p 0 = N 0 N for some 1 ≤ N 0 ≤ N . Then it is well known from the theory of order statistics (see for example e.g. [1, formula (2.2.20) page 14]) that

Expectation of W .. We will prove that it is equal to zero by induction on n 0 . For n 0 = 1, F ( L1 ) has the same law as U (N -N 0 ) , thus the result is obvious by equation (A.9). Now assume that [F ( Ln 0 -1 )] = 1p n 0 -1 0 . From the decomposition

which proves that W has zero mean.

imsart-aos ver. 2010/04/27 file: cdfg.tex date: July 2, 2010

Variance of W .. From the proof of Theorem 3.1, we know that