Henri Gouin

Utilization of the second gradient theory in continuum mechanics to study motions and thermodynamics of liquid-vapor interfaces

Keywords: van der Waals fluids, second gradient theory, capillarity, surface tension, Marangoni effect PACS: 47.55.pf, 47.61.-k, 64.60.Bd, 64.60.Ht, 68.03.Cd, 68.03.Hj, 68.05.-n, 68.35.Gy 1991 MSC: 76A02, 76A20, 76B45

A thermomechanical model of continuous fluid media based on second gradient theory is used to study motions in liquid-vapor interfaces. At equilibrium, the model is shown to be equivalent to mean-field molecular theories of capillarity. In such fluids, conservative motions verify first integrals that lead to Kelvin circulation theorems and potential equations. The dynamical surface tension of liquid-vapor interfaces is deduced from viscous fluid equations. The result provides and explains the Marangoni effect.

Introduction

Usually, liquid-vapor interfaces are represented by material surfaces endowed with an energy related to the Laplace surface tension. Several studies conducted in the fields of fluid mechanics and thermodynamics represent the Email address: henri.gouin@univ-cezanne.fr (Henri Gouin).

Preprint submitted to ArXiv & Hal

interface as a surface separating two media and having its own energy density and characteristic behaviors [START_REF] Defay | Surface Thermodynamics[END_REF][START_REF] Sanfeld | Thermodynamics of Surfaces[END_REF][START_REF] Emschwiller | Chimie Physique[END_REF][START_REF] Davis | Rupture of thin liquid films in Waves on Fluid Interface[END_REF][START_REF] Landau | Fluids Mechanics[END_REF][START_REF] Sen | Steady thermocapillarity flows in two-dimensional slots[END_REF][START_REF] Ono | Molecular theory of surface tension in liquid[END_REF][START_REF] Levitch | Physicochemical Hydrodynamics[END_REF]. When working far from the critical temperature, detailed theoretical or experimental studies show that the capillary layer thickness ranges over a few molecular beams [START_REF] Emschwiller | Chimie Physique[END_REF][START_REF] Levitch | Physicochemical Hydrodynamics[END_REF][START_REF] Rocard | Thermodynamique[END_REF]. For fluid conditions defined by temperature and density, molecular models such as those used in gas kinetic theory lead to state laws associated with non-convex internal energies as the one associated with the van der Waals model for the pressure [START_REF] Rocard | Thermodynamique[END_REF][START_REF] Cahn | Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid[END_REF][START_REF] Dunn | On the thermodynamics of interstitial working[END_REF][START_REF] Brin | Contribution à l'étude de la couche capillaire et de la pression osmotique[END_REF][START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF][START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF]. Such models provide a direct verification of the Maxwell rule applied to isothermal liquid-vapor transition [START_REF] Aifantis | The mechanical theory of fluid interfaces and Maxwell's rule[END_REF]; nonetheless, they present two disadvantages: First, for density values ranging between the vapor and liquid bulk densities, the pressure may become negative; however, traction created by simple physical experiments can lead to such pressure values [START_REF] Rocard | Thermodynamique[END_REF][START_REF] Brin | Contribution à l'étude de la couche capillaire et de la pression osmotique[END_REF][START_REF] Bruhat | Cours de Physique Generale[END_REF]. Second, in the field between vapor and liquid, the internal energy distribution cannot be represented by a convex surface associated with the density and entropy and seems to be inconsistent with the existence of equilibrium state. To overcome these disadvantages, the thermodynamical investigation uses a plane domain in place of the portion corresponding to the non-convex internal energy. The fluid can no longer be considered as a continuous medium. The interface is represented by a material surface with a null thickness. In such a case, the only possible representation of the dynamic behavior of the interface is a discontinuous surface and its essential structure remains unknown. At equilibrium, the disadvantages can be eliminated by appropriately modifying the stress tensor of the capillary layer which is then expressed in an anisotropic form. As a consequence, the energy of the continuous medium must be modified [START_REF] Cahn | Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid[END_REF][START_REF]La theorie du second gradient et la capillarite[END_REF][START_REF] Bongiorno | Molecular theory of fluid interfaces[END_REF].

If the fluid interface must be represented as a continuous medium, how can the stress tensor of the capillary layer be written in a dynamic expression [START_REF] Davis | Rupture of thin liquid films in Waves on Fluid Interface[END_REF]? In previous articles, we proposed a dynamic theory named internal capillarity based on the second gradient of deformation of the medium [START_REF] Casal | Connection between the energy equation and the motion equation in korteweg theory of capillarity[END_REF][START_REF] Casal | Kelvin's theorems and potential equation in Korteweg's theory of capillarity[END_REF][START_REF] Gouin | Dynamical surface tension and Marangoni effect for liquid-vapor interfaces in internal capillarity theory[END_REF]. Elaborated in the sixties for the static case, the theory points out the advantage of using a three-dimensional approach of the capillarity in a continuous mechanical model [START_REF]La theorie du second gradient et la capillarite[END_REF][START_REF]Capillarité interne en mécanique des milieux continus[END_REF]. The deduced equations of equilibrium provide a satisfying representation of isothermal liquid-vapor equilibrium states. The approach is not new, and dates back to van der Waals and Korteweg [START_REF] Van Der Waals | Thermodynamique de la capillarité dans l'hypothèse d'une variation continue de densité[END_REF][START_REF] Korteweg | Sur la forme que prennent les équations du mouvement des fluides si l'on tient compte des forces capillaires[END_REF]. Cahn and Hilliard again used such an approach in 1959 in reference to the free energy [START_REF] Cahn | Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid[END_REF], and went back to a Landau-Ginsburg model [START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF]. The representation proposed in the present study is based on the notion of internal energy which is more useful when the temperature is not uniform. The apparent contradiction between the Korteweg classical stress theory and the Clausius-Duhem inequality makes the dynamic study of phase transformation more difficult [START_REF] Gurtin | Thermodynamics and the possibility of spatial interaction in elastic materials[END_REF], but the use of the capillary fluid model rectifies this apparent anomaly [START_REF] Dunn | On the thermodynamics of interstitial working[END_REF][START_REF] Casal | Connection between the energy equation and the motion equation in korteweg theory of capillarity[END_REF]. In the more general case of internal capillarity, by representing energy in terms of second gradient of deformation [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus[END_REF] and by applying a simple al-gebraic identity, it is possible to draw a relationship between energy equation, motion equation, mass conservation equation and entropy [START_REF] Casal | Connection between the energy equation and the motion equation in korteweg theory of capillarity[END_REF]. For a conservative fluid, an additional term with a heat flux dimension must be introduced into the energy equation. In the case of viscous fluid, previous results provide a set of equations which do not modify the Clausius-Duhem inequality and make them compatible with the second law of thermodynamics. In the non-viscous case, classical fluid motions and motions of fluids endowed with internal capillarity reveal a common structure that induces a thermodynamic form of equation of motion and potential equations [START_REF]La theorie du second gradient et la capillarite[END_REF][START_REF]Principes variationnels en fluide compressible et en magnétodynamique des fluides[END_REF][START_REF] Selinger | Variational principles in continuum mechanics[END_REF]. This leads to the same classification of motions, generalization of Kelvin theorems, Crocco-Vaszonyi equation and first integrals conserved across the interfaces [START_REF] Gouin | Contribution à une étude géometrique et variationnelle des milieux continus[END_REF][START_REF] Gouin | Exemples of non-conservative perfect fluid motions[END_REF]. By representing the internal energy as a function of entropy, density, density gradient and by using a single constant λ, the resulting equations are thereby identical to those obtained with molecular models of mean field theory in the isothermal case [START_REF] Rocard | Thermodynamique[END_REF].

For a surface area relatively large with respect to the capillary layer thickness, the surface tension is calculated by integration throughout the interface. Lengthwise, the surface tension is not constant and depends on the dynamical distribution of density and temperature. This dynamical distribution based on equations of motion, calls upon a Navier-Stokes viscosity. For interfaces in isothermal equilibrium, the results are classical [START_REF] Rocard | Thermodynamique[END_REF][START_REF] Cahn | Free energy of a nonuniform system. III. Nucleation in a two-component incompressible fluid[END_REF][START_REF] Brin | Contribution à l'étude de la couche capillaire et de la pression osmotique[END_REF][START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF][START_REF] Bongiorno | Molecular theory of fluid interfaces[END_REF]. The study of motion in the interface without mass transfer needs the surface tension gradient and the velocity gradients associated with the dynamic viscosity. The ensuing Marangoni effect has been interpreted by using a limit analysis wherein the approximate quantities correspond to the physical dimensions of the interface. When the mass flow across the interface is not null, a general dynamic form of the Laplace equation is given [START_REF] Gouin | Dynamical surface tension and Marangoni effect for liquid-vapor interfaces in internal capillarity theory[END_REF]. The method herein presented is completely different from classical calculations based on the balance equations established for both sides of a discontinuous surface, which takes account of density variations only as a difference across the interface [START_REF] Ishi | Thermo-fluid Dynamic. Theory of Two-phase Flow[END_REF]. In the particular case of isothermal liquid-vapor equilibria, an invariant integral of motions compatible with the interface coincides with a generalization of the rule advanced by Maxwell, associated with phase transitions [START_REF] Aifantis | The mechanical theory of fluid interfaces and Maxwell's rule[END_REF].

The model of a viscous fluid endowed with internal capillarity is therefore substantiated by the consequences verified in both equilibrium states and non-isothermal motions. The model provides a better understanding of the behavior of liquid-vapor interfaces in motion. It at least gives a partial answer to the question: Is the fluid at the interface rigid or moving ? [START_REF] Birkhoff | Numerical fluid dynamics[END_REF] and proposes a theory which takes account of the stress tensor and dynamics in the structure of the interface [START_REF] Davis | Rupture of thin liquid films in Waves on Fluid Interface[END_REF]. The resulting behaviors are not the laws governing classical fluids since they include an anisotropic stress tensor in the momentum equation and an additional heat flux term in the energy equation.

2 Equations of motion for a fluid endowed with internal capillarity

Case of conservative fluids

The second gradient theory [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus[END_REF], conceptually more straightforward than the Laplace theory, can be used to construct a theory of capillarity. In the present text the only addition with respect to compressible fluids is an internal mass energy which is a function of density ρ, entropy s and also grad ρ. The specific internal energy ε characterizes both the compressibility and capillarity properties of the fluid [START_REF] Pratz | Contribution à la théorie du second gradient pour les milieux isotropes[END_REF], ε = f (ρ, s, β), where β = (grad ρ) 2 . The equation of motion for a conservative fluid writes

ρ a = div σ -ρ grad Ω, (1)
where a denotes the acceleration vector, Ω the extraneous force potential and σ the general stress tensor

σ = -p I -λ (grad ρ)(grad ρ) T (2)
or

σ ij = -p δ ij -λ ρ ,i ρ ,j , i, j ∈ {1, 2, 3} where λ = 2 ρ ε ′ β and p = ρ 2 ε ′ ρ -div(λ grad ρ) (3) It should be noted that: θ = ε ′ s ( 4
) is the Kelvin temperature expressed as a function of ρ, s and β.

Proof : The equation of motion is obtained by using the virtual work principle. The virtual displacement, denoted by δx, has been defined by J. Serrin ([34], p. 145). For a fluid endowed with internal capillarity, the virtual work principle ( [START_REF] Serrin | Mathematical principles of classical fluid dynamics[END_REF], iv, section 14) is stated as follows: the motion of a fluid is such that:

δ V ρ (ε + Ω) dv - V ρ a δx dv = 0,
where V is an arbitrary material volume. The variation of entropy ( [START_REF] Serrin | Mathematical principles of classical fluid dynamics[END_REF], p. 148) must satisfied the condition δs = 0; while taking account of

δ ∂ρ ∂x = ∂δρ ∂x - ∂ρ ∂x ∂δx ∂x ,
we obtain

δβ = 2 ∂δρ ∂x - ∂ρ ∂x ∂δx ∂x grad ρ and δε = ε ′ ρ δρ + ε ′ s δs + ε ′
β δβ, as well as eqs. (14.5) and (14.6) in [START_REF] Serrin | Mathematical principles of classical fluid dynamics[END_REF]. For virtual displacement δx null at the edge of V and integration by parts, the integrals calculated on the edge of V have a null contribution; we obtain:

V (ρ a -div σ + ρ grad Ω) δx dv = 0,
where σ is expressed by Eq. ( 2). Classical methods of the calculus of variations lead to Eq. ( 1). It appears that a single term λ accounts for the effects of the second gradient deformation in the equation of motion. Scalar λ depends on ρ, s and β as ε does. In a study of surface tension based on gas kinetic theory, Rocard obtained the same expression as Eq. ( 2) for the stress tensor but with λ constant. If λ is constant, specific internal energy ε writes:

ε(ρ, s, β) = α(ρ, s) + λ 2ρ β.
The second gradient term λ β/(2ρ) is simply added to the specific internal energy α(ρ, s) of the classical compressible fluid. The pressure of the compressible fluid is denoted P ≡ ρ 2 α ′ ρ and the temperature is denoted Θ ≡ α ′ s .

Consequently,

p = P -λ β 2 + ρ ∆ρ .
For the pressure P , Rocard uses either the van der Waal pressure

P = ρ R Θ 1 -bρ -aρ 2
or other laws of which he provides a comparison. It should be noted that if λ is constant, this implies that θ = Θ and there exits a relationship, independent of the second gradient terms, between θ, ρ and s.

Case of viscous fluids

If the fluid is endowed with viscosity, the equation of motion includes not only the stress tensor σ, but also the viscous stress tensor σ v :

σ v = ν tr D + 2 µ D,
where D is the deformation tensor, symmetric gradient of the velocity field u

D = 1 2    ∂u ∂x + ∂u ∂x T    .
Of course in second gradient theory, it would be coherent to add terms accounting for the influence of higher order derivatives of the velocity field to the viscous stress tensor σ v ; this has not been done in the present case. Equation ( 2) is modified by adding the virtual work of the forces of viscosity and we obtain:

ρ a = div(σ + σ v ) -ρ grad Ω.

Equation of energy for viscous fluids endowed with internal capillarity

Let e be total volumic energy, e ≡ ρ 1 2 u 2 + ε + Ω , q the heat flux vector, r the heat supply and h = ε + p ρ the specific enthalpy. We define the following quantities:

M = ρ a -div(σ + σ v ) + ρ grad Ω, B = • ρ + ρ div u, S = ρ θ • s + div q -r -tr(σ v D), E = ∂e ∂t + div [(e -σ -σ v )u] -div(λ • ρ grad ρ) + div q -r -ρ ∂Ω ∂t .
Theorem: For an internal energy written as in Eq. ( 1) and for any motion in the fluid, the relation

E -M • u - 1 2 u 2 + h + Ω B -S ≡ 0 (5)
is an identity.

Proof : In the first member, the dissipative terms q, r and σ v cancel out; this also is the case for the extraneous force potential and the inertia terms. After having replaced σ, p and θ by their respective values in Eqs (3 -6), it remains to prove that the terms from internal energy ε also cancel out. These terms include the following expressions:

a) in E : ∂(ρε) ∂t + div(ρ ε u) -div(σu) -div(λ • ρ grad ρ) ≡ ρ • ε + ε( • ρ + ρ div u) -div(σ)u -tr σ ∂u ∂x -div(λ • ρ grad ρ),
or finally,

≡ ρ(ε ′ s • s + ε ′ ρ • ρ + ε ′ β • β) + ε( • ρ + ρ div u) -div(σ)u + p div u +λ (grad ρ) T ∂u ∂x grad ρ -div(λ • ρ grad ρ); b) in M • u : -div(σ)u; c) in 1 2 u 2 + h + Ω B : ε + p ρ ( • ρ + ρ div u); d) in S : ρ ε ′ s • s.
So, the first member of Eq. ( 5) writes:

ρ ε ′ ρ • ρ + λ 2 • β - p ρ • ρ + λ (grad ρ) T ∂V ∂x grad ρ -div(λ • ρ grad ρ).
Finally, by taking Eq. ( 5) into account, we obtain that

λ (grad ρ) T d(grad ρ) dt + • ρ div λ (grad ρ) + λ (grad ρ) T ∂u ∂x grad ρ -div(λ • ρ grad ρ)
is identically null.

Corollary 1 : For a conservative motion of a fluid endowed with internal capillarity, the conservation of entropy along the trajectories is equivalent to

∂e ∂t + div[(e -σ)u] -div(λ • ρ grad ρ) -ρ ∂Ω ∂t = 0, (6)
which is the equation of conservation of energy. This is derived from Eq. ( 5) when σ v = 0, div q = r = 0 and M = 0, B = 0, S = 0. It leads us to add the additional term div(λ • ρ grad ρ) to the equation of energy. The vector λ

• ρ grad ρ has the dimension of a heat flux vector and even occurs in the conservative case.

Corollary 2 : For a motion of a viscous fluid endowed with internal capillarity, the energy equation

∂e ∂t + div[(e -σ -σ v )u] -div(λ • ρ grad ρ) + div q -r -ρ ∂Ω ∂t = 0
is equivalent to the entropy equation

ρθ • s + div q -r -tr(σ v D) = 0. (7)
Equation ( 7) corresponds to the classical version of the entropy variation expressed by the function of dissipation of viscous stresses ψ = tr(σ v D).

The Planck inequality and the Clausius -Duhem inequality

For any motion in a viscous fluid endowed with capillarity, it has been assumed that ρ θ

• s + div q -r 0 . (8)
Inequality ( 8) represents the Planck inequality [START_REF] Truesdell | Rational Thermodynamics[END_REF]. Let generally write the Fourier principle as q • grad θ 0 .

Then, the Clausius-Duhem inequality can be directly deduced as

ρ • s + div q θ - r θ 0 .
When in motion, the fluid is endowed with internal capillarity, the law of thermodynamics leads to the existence of a heat flux vector. This heat flux vector introduces an additional term into the classical equation of energy, even if the fluid is non-viscous. It is possible to extend these results to continuous media with an internal energy that contains gradient terms of deformation.

Transformation of the motion equations for fluids endowed with internal capillarity

In this section, we rewrite Eq. ( 2) in other forms.

Case λ constant

Let us note ω = Ω -λ ∆ρ. Then, Eq. ( 2) can be written

ρ a + grad P + ρ grad ω = 0 . (9)
This is a representation of a perfect fluid where P is the van der Waals pressure; term ω contains all capillarity terms. Equations ( 3) and ( 5) yield

σ ij = -P δ ij + λ 1 2 ρ ,k ρ ,k + ρρ, kk δ ij -λ ρ ,i ρ ,j .
Then, σ ij,j = -P ,i + λ ρρ ,ijj i.e. div σ = -grad P + λ ρ grad(∆ρ) . We can note that for viscous fluids, Eq. ( 9) writes ρ a + grad P + ρ grad(Ω -λ ∆ρ) -div σ v = 0 .

(10)

Thermodynamic form of the equation of motion

Commonly -and not only when λ is constant -the equation of motion (1) can be written in the following form

a = θ grad s -grad(h + Ω), (11)
which is the thermodynamic form of the equation of motion. In the case without capillarity (ε = α(ρ, s), ε ′ β = 0), Eq. ( 11) is well-known ( [START_REF] Serrin | Mathematical principles of classical fluid dynamics[END_REF], p. 171). The equation remains valid for fluids endowed with internal capillarity.

The verification simply consists on writing out the second member of Eq. ( 11). Let us denote

A = θ grad s -grad h i.e. A i = θs ,i -h ,i .
According to Eq. ( 4), we get:

A i = ε ′ s s ,i -ε + p ρ ,i
. Then,

A i = ε ′ s s ,i -ε ′ s s ,i + ε ′ ρ ρ ,i + ε ′ β ρ ,k ρ ,k ,i - p ρ 2 ρ ,i - 1 ρ p ,i . (12)
Using Eq. (3), Eq. ( 12) writes σ ij,j /ρ and consequently A = div σ/ρ, which proves Eq. ( 11). Corollary 1 : For a homentropic flow [START_REF]Principes variationnels en fluide compressible et en magnétodynamique des fluides[END_REF], the velocity circulation along a closed fluid curve is constant.

Theorem 2 : The velocity circulation along a closed isothermal fluid curve is constant.

Corollary 2 : For an isothermal flow, the velocity circulation along a closed fluid curve is constant.

Potential equations for conservative fluids and classification of motions

The results of [START_REF]La theorie du second gradient et la capillarite[END_REF][START_REF]Principes variationnels en fluide compressible et en magnétodynamique des fluides[END_REF][START_REF] Selinger | Variational principles in continuum mechanics[END_REF][START_REF] Gouin | Contribution à une étude géometrique et variationnelle des milieux continus[END_REF] can be applied to conservative motions of fluids endowed with internal capillarity. To any motion of fluids endowed with internal capillarity, there correspond scalar potentials ϕ, ψ, τ and χ verifying

                           • ϕ = 1 2 u 2 -h -Ω, • τ = 0, • ψ = θ, • χ = 0, • s = 0. (13)
The fluid velocity is given by

u = grad ϕ + ψ grad s + τ grad χ, (14)
together with ∂ρ ∂t + div(ρ u) = 0.

The equations can be used to classify motions in the same way as in the case of non-capillary perfect fluids [START_REF]Principes variationnels en fluide compressible et en magnétodynamique des fluides[END_REF][START_REF] Gouin | Exemples of non-conservative perfect fluid motions[END_REF].

Homentropic motions

This case corresponds to s constant in all the fluid. Eq. ( 14) writes u = grad ϕ + τ grad χ, which leads to the Cauchy theorem [START_REF] Serrin | Mathematical principles of classical fluid dynamics[END_REF] d dt

rot u ρ = ∂u ∂x rot u ρ .

Oligotropic motions

Equation ( 14) writes u = grad ϕ + ψ grad s. The flows verify the relation rot u • grad s = 0. Surfaces with equal entropy are eddy surfaces. The velocity circulation along a closed isentropic fluid curve is null.

Remark: the form of Eqs. ( 9) and ( 11) allows to generalize the results of [START_REF] Gouin | Exemples of non-conservative perfect fluid motions[END_REF].

Generalized Crocco -Vazsonyi equation

The equation of energy [START_REF] Sen | Steady thermocapillarity flows in two-dimensional slots[END_REF]

The first member is equal to ρ

• H+H(

• ρ+ρ div u). Because of mass conservation, the definitive form of the equation of energy writes

ρ • H = ρ ∂Ω ∂t + ∂p ∂t + div λ ∂ρ ∂t grad ρ .
By taking account of the identity

a = ∂u ∂t + rot ∧ u + grad 1 2 u 2 ,
Equation ( 11) can be written

∂u ∂t + rot u ∧ u = θ grad s -grad H. ( 17
)
If the motion is a steady flow, Eq. ( 17) makes it possible to conclude that H is constant along the stream lines (

• H = 0) and Eq. ( 16) shows that the divergence of the partial energy flux ρHu is null. Equation ( 17) leads to the Crocco-Vazsonyi equation generalized to fluids endowed with internal capillarity; in stationary motion,

rot u ∧ u = θ grad s -grad H,
where s and H are constant along every stream line.

Remark : Equations [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF] show that along each trajectory, and consequently across interfaces, τ , χ and the entropy s are constant scalars. They represent first integrals of the motion. Other integrals can obviously be found as the Kelvin integrals

J = Γ τ (grad χ) T dx ≡ Γ τ dχ
are constant along any fluid curve Γ. Thanks to the Noether theorem, any law of conservation can be represented with an invariance group [START_REF] Gouin | Noether theorem in fluid mechanics[END_REF]. The law of conservation expressed with the Kelvin theorems corresponding to isentropic fluid curves is related to the group of permutations associated with particles of equal entropy; the group keeps the equations of motion invariant for both classical perfect fluid and fluid endowed with internal capillarity, but also for any perfect fluid endowed with an internal energy depending on ρ, s and their gradients of any order. It is even tempting to define a general perfect fluid by identifying it with an invariant group or with a continuous medium whose motions verify the Kelvin theorems [START_REF] Gouin | Thermodynamic form of the equation of motion for perfect fluid of grade n[END_REF].

Dynamic surface tension of liquid-vapor interfaces

Far from the critical conditions, the thickness of a liquid-vapor interface is very small [START_REF] Rocard | Thermodynamique[END_REF][START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF]. Outside the capillary layer, density and its spatial derivatives have smooth variations. The density in each phase is reached at points located within the immediate vicinity of the layer. In our study, bubbles and drops of size of a few molecular beams have not been considered. Surfaces of equal density materializing the interface are parallel surfaces and are used to define a system of orthogonal coordinates [START_REF] Valiron | Equations fonctionnelles: Applications[END_REF]. Notations are those used in [START_REF] Germain | Mécanique des Milieux Continus[END_REF]. The subscript 3 refers to the direction normal to surfaces of equal density. Let us denote by e 3 the unit vector of this index oriented in the increasing density direction. The mass flow through the interface is assumed to be null such that surfaces of equal density are material surfaces. In the capillary layer, div u = 0.

(

Extraneous forces being neglected, Eq. ( 18) represents velocities compatible with the interface. Equation ( 10) yields:

ρ a 3 + 1 h 3 ∂P ∂x 3 = λ ρ 1 h 3 ∂∆ρ ∂x 3 + 2µ (div D) 3 , (19)
ρ a tg + grad tg P = λ ρ grad tg ∆ρ + 2µ (div D) tg , where the subscript tg denotes the tangent component to the surfaces of equal density. The different linear sizes of the interface should be taken in consideration. The capillary layer is measured in Angström and the surface curvature radii are of non-molecular dimensions. The deduced relations result from a limit analysis where the parameter related to the thickness of the capillary layer tends to zero.

The subscripts v and l designating the vapor and the liquid bulks, respectively, integration of Eq. ( 19) along the third coordinate line yields, when we assume that a 3 is bounded through the interface such that its integral is negligible,

P -P v = λ ρ ∆ρ -λ x 3 x v 3 ∆ρ ∂ρ ∂x 3 dx 3 + 2 x 3 x v 3 µ (div D) 3 h 3 dx 3 . (20)
The partial derivatives of velocity with respect to coordinates x 1 and x 2 are assumed to be bounded. Taking account of Eq. ( 18), the last term in Eq. ( 20) is negligible. Furthemore, we have the relation:

∆ρ = - 2 R m 1 h 3 ∂ρ ∂x 3 + 1 h 3 1 h 3 ∂ρ ∂x 3 ,3 , (21)
where R m denotes the mean curvature of the surfaces of equal density in the capillary layer, oriented by e 3 [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus[END_REF]. It follows:

P -P v = λ ρ∆ρ - 1 2 (grad ρ) 2 + 2 λ R m x 3 x v 3 (grad ρ) 2 h 3 dx 3 . (22)
Denoting dn ≡ h 3 dx 3 , we finally obtain:

P -P v = 2 H R m where H = λ n l nv (grad ρ) 2 dn,
which expresses the Laplace equation for motions. We interpret H as the dynamical surface tension.

Remark: Interpretation of the surface tension at equilibrium. The eigenvalues of the internal capillary stress tensor are deduced from Eq.

(2): λ 1 = -p + λ (grad ρ) 2 is the eigenvalue associated with the plane perpendicular to grad ρ. λ 2 = -p is the eigenvalue associated with the direction grad ρ.

In the system of curvilinear coordinates related to the interface, the stress tensor writes

σ =        λ 1 0 0 0 λ 1 0 0 0 λ 2       
.

The equilibrium equation for the plane interface is drawn from Eq. ( 1); neglecting extraneous forces we get:

λ 2 = -P 0 ,
where P 0 denotes the pressure in the liquid and vapor bulks. The force per unit length along the edge of the interface (see Fig. 1) is:

F = x l 3 x v 3 λ 1 h 3 dx 3 = -P 0 L + x l 3 x v 3 λ (grad ρ) 2 h 3 dx 3 ,
where L denotes the capillary layer thickness. Obviously, P 0 L is negligible. Let us denote

H = x l 3 x v 3 λ (grad ρ) 2 h 3 dx 3 ,
the force per unit length; then H represents the surface tension of the plane interface at equilibrium.

Practical calculation of the surface tension

By integrating Eq. ( 19) on the third coordinate line, the viscosity term being assumed negligible, we obtain

λ ∆ρ = ∂ ∂ρ ρ ρ ρ v P -P v ρ 2 dρ . (23)
Let us notice that ∂P ∂θ ∂θ ∂x 3 is negligible with respect to ∂P ∂ρ ∂ρ ∂x 3 .

Taking Eq. ( 21) into account, a further integration on the third coordinate leads to

λ 2 (grad ρ) 2 = 2 λ R m x 3 x v 3 (grad ρ) 2 h 3 dx 3 + ρ ρ ρ v P -P v ρ 2 dρ . (24)
By denoting x i 3 the third coordinate of the surface of equal density ρ = ρ i , then for

x 3 ∈ [x v 3 , x i 3 ], the term 2 λ R m x 3 x v 3 (grad ρ) 2 h 3 dx 3 is negligible with respect to λ 2 (grad ρ) 2 ; that is the case for ρ i = 1 2 (ρ v + ρ l ). For ρ ∈ [ρ v , ρ i ],
Eq. ( 24) writes λ 2 (grad ρ) 2 = ρ ρ ρ v P -P v ρ 2 dρ. We obtain the same result for ρ ∈ [ρ i , ρ l ] and finally,

H = √ 2λ ρ i ρ v u u ρ v P -P v ρ 2 dρ du + ρ l ρ i u u ρ i P -P l ρ 2 dρ du . (25)
Relation [START_REF] Gurtin | Thermodynamics and the possibility of spatial interaction in elastic materials[END_REF] expresses the surface tension of a liquid-vapor interface for motions compatible with the interface. Here P is not only a function of density but also temperature and varies according to the location of the point within the interface. The viscosity has explicitly disappeared in the expression of H (the expression of H does not explicitly take the viscosity into account). In the case of plane interface in isothermal equilibrium, the expression given by Rocard [START_REF] Rocard | Thermodynamique[END_REF] is found again. The value of the internal capillarity constant can be numerically calculated on the basis of experimental values of H and expressions of P . By injecting the value of λ in Eq. ( 25), we can calculate surface tension for any dynamic temperature distribution.

Remark: Eq. ( 23) yields

λ ∆ρ = P ρ - P v ρ v + ρ ρ v P ρ 2 dρ
and we deduce

ρ l ρ v P ρ 2 dρ = P v ρ v - P l
ρ l which represents an invariant corresponding to motions compatible with the interface. This invariant leads to the equations (4-11) of [START_REF] Aifantis | The mechanical theory of fluid interfaces and Maxwell's rule[END_REF] demonstrated in the specific case of an interface in isothermal equilibrium. The case of plane interface leads to the Maxwell equal area rule [START_REF] Rocard | Thermodynamique[END_REF][START_REF] Aifantis | The mechanical theory of fluid interfaces and Maxwell's rule[END_REF].

Marangoni effect in liquid-vapor interfaces

Let us denote

H j (x 3 ) = x 3 x j 3 (grad ρ) 2 h 3 dx 3 ,
where j ∈ {v, l}.

By using calculations as in Section 10, for example for

x 3 ∈ [x v 3 , x l 3 ], we obtain H v (x 3 ) = √ 2λ ρ ρ v u u ρ v P -P v ρ 2 dρ du .
By injecting Eq. ( 22) into Eq. ( 23), we obtain:

ρ a tg = λ 2 grad tg (grad ρ) 2 -grad tg 2 R m H j (x 3 ) -grad tg P j +2µ(div D) tg . (26)
Equation ( 18) and conditions on partial derivatives of velocity imply:

x l 3 x v 3 (div D) tg h 3 dx 3 = [De 3 ] x l 3 x v 3 .
Integration of Eq. ( 26) yields:

λ 2 x l 3 x v 3 grad tg (grad ρ) 2 h 3 dx 3 + 2[µDe 3 ] x l 3 x v 3 = 0 .
Taking account of λ 2

x l 3 x v 3 grad tg (grad ρ) 2 h 3 dx 3 = grad tg λ x l 3 x v 3 (grad ρ) 2 h 3 dx 3 ,
we can conclude grad tg H + 2[µDe 3 ] l v = 0 . If the viscosity stress tensor of the vapor bulk is supposed to be negligible, we get:

grad tg H + 2[µD l e 3 ] = 0 . This is the usual Marangoni condition for free boundary problems [START_REF] Sen | Steady thermocapillarity flows in two-dimensional slots[END_REF]. In the limit case when the viscosity coefficient is null (as for superfluid helium), the problem must be proposed in another way: the momentum associated with the interface can no longer be neglected and other physical effects must be taken into account [START_REF] Landau | Fluids Mechanics[END_REF]. We must note there exist other phenomenal presentations of the Marangoni effect (for example [START_REF] Scriven | Dynamics of a fluid interface. Equation of motion for Newtonian surface fluids[END_REF]), but to our knowledge, they all consider the interface as a discontinuous surface of the fluid medium. The calculations performed in the present text do not call upon the use of any linear approximations, and only take account of the various physical quantities describing liquid-vapor interfaces while working far from the critical temperature. The use of second gradient theory for representing interfaces has, of course, been considered by several authors. In [START_REF] Dunn | On the thermodynamics of interstitial working[END_REF][START_REF] Slemrod | Admissibility criteria for propagating phase boundaries in a van der Waals fluid[END_REF][START_REF] Slemrod | Dynamic phase transitions in a van der Waals fluid[END_REF][START_REF] Dunn | Interstitial working and a nonclassical continuum thermodynamics in New perspectives in thermodynamics[END_REF], their source is found to be in free energy which is more directly useable in case of isothermal flows.

Conclusion and further developments

As we have seen in this revisited version of [START_REF] Gouin | [END_REF], the second gradient theory, conceptually more straightforward than the Laplace theory can be used to build a theory of capillarity. Such a theory is able to take account of systems in which fluid interfaces are present. The internal capillarity associated with fluids is the simplest case. A mathematical limit analysis associated with the thickness of the interface when the size of the layer goes toward zero and the behavior of the layer between fluid phases yield the model of material surfaces. Such a theory is able to calculate the superficial tension as well in the case of thin interfaces as thick interfaces. The static model in continuum mechanics of second gradient theory is extended to dynamics. The equation of motion is able to induce a stress tensor. In the fluid case, the theory does not lead to an isotropic stress tensor.

Many developments have been done since paper [START_REF] Gouin | [END_REF] by the author:

It is possible to obtain the radius of nucleation of microscopic drops and bubbles and to develop a macroscopic theory as the Laplace theory of capillarity for curved interfaces [START_REF] Dell'isola | Radius and surface tension of microscopic bubbles by second gradient theory[END_REF]. The stability of interfaces is investigated with differential or partial derivative equations [START_REF] Gouin | Stability of spherical isothermal liquid-vapor interfaces[END_REF]. Classification of fluids endowed with internal capillarity is the extension of the classical fluid one's [START_REF] Casal | Invariance properties of perfect fluids of grade n[END_REF].

In fact, contact forces are of a different nature than the ones associated with the Cauchy stress tensor. Classical conditions with the tetrahedron construction due to Cauchy are not efficient to study the non-linear behavior of fluids endowed with internal capillarity. We deduced contact forces concentrated on edges representing the boundaries of surfaces of separation. For example, such conditions are necessary to study the stability of thin films in contact with solid walls and the connection with the mean field molecular theories [START_REF] Gouin | Energy of interaction between solid surface and liquids[END_REF][START_REF] Gouin | Boundary conditions for a capillary fluid in contact with a wall[END_REF][START_REF] Gouin | Wetting problem for multicomponent fluid mixtures[END_REF][START_REF] Gouin | Elastic effects of liquids on surface physics[END_REF][START_REF] Gouin | A mechanical model for the disjoining pressure[END_REF].

The model dynamically interprets phenomena in the vicinity of the critical point and authorizes an investigation -at least qualitative -of the dynamic change of phases between bulks in fluids or fluids mixtures [START_REF] Gouin | Dynamic effects in gradient theory for fluid mixtures[END_REF][START_REF] Gouin | Mixtures of fluids involving entropy gradients and accelerations waves in interfacial layers[END_REF]. New mathematical equations associated with the hyperbolicity in continuous media can be considered in second gradient theory and consequently in internal capillarity [START_REF] Gavrilyuk | A new form of governing equations of fluids arising from Hamilton's principle[END_REF][START_REF] Gavrilyuk | Symmetric form of governing equations for capillary fluids in Trends in applications of mathematics to mechanics[END_REF].

More recently, we must notice that fluids endowed with internal capillarity are able to model fluid layers at a scale of some nanometers and to recognize the disjoining pressure concept for very thin liquid films [START_REF] Gouin | A mechanical model for the disjoining pressure[END_REF]. The model can be applied in vegetal biology to interpret the ascent of the crude sap in very high trees as sequoia and giant eucalyptus [START_REF] Gouin | A new approach for the limit to tree height using a liquid nanolayer model[END_REF][START_REF] Gouin | Liquid-solid interaction at nanoscale and its application in vegetal biology[END_REF]

6

  Generalized Kelvin theoremsLet J be the circulation of the velocity vector along a closed fluid curve C convected by the fluid flow J = and consequently yields the following theorems which are valid for fluids endowed with internal capillarity: Theorem 1 : The velocity circulation along a closed isentropic fluid curve is constant.

Fig. 1 .

 1 Fig. 1. Interpretation of the surface tension in liquid-vapor interfaces as a layer of fluid endowed with internal capillarity.


Marco Kletting

Michel Kieffer

Eric Walter

Two Approaches for Guaranteed State Estimation of Nonlinear Continuous-Time Models

HAL is

Introduction

When building knowledge-based models, for instance models based on the laws of physics, one frequently ends up with a continuous-time state-space model, which may depend on a possibly time-varying vector of parameters p ∈ R n p : ẋ(t) = f x (x (t) , p (t)) ,

Marco Kletting

Multi-Function Airborne Radars (OPES22), Cassidian Electronics,Woerthstr. 85, 89077 Ulm, Germany e-mail: marco.kletting@cassidian.com

Michel Kieffer

Laboratoire des Signaux et Systèmes -CNRS -SUPELEC -Univ Paris-Sud, 3 rue Joliot-Curie, F-91192 Gif-sur-Yvette cedex, on leave at LTCI -CNRS -Telecom ParisTech, 46 rue Barault, F-75013 Paris e-mail: michel.kieffer@lss.supelec.fr

Eric Walter

Laboratoire des Signaux et Systèmes -CNRS -SUPELEC -Univ Paris-Sud, 3 rue Joliot-Curie, F-91192 Gif-sur-Yvette cedex e-mail: eric.walter@lss.supelec.fr 1 where x ∈ R n x is the state vector. Provided that an equation is specified for the dynamics of p, such as

ṗ(t) = ∆ p with ∆ p ∈ [∆ p, ∆ p],
one can concatenate the state and parameter vectors into an extended state vector z(t) = [x(t) T , p(t) T ] T . Then ż(t) = f(z(t)),

where

f = f x (x(t), p(t)) ∆ p with f : D → R n , D ⊂ R n = R n x × R n p .
Any time-invariant parameter p i should satisfy ∆ p i = 0. This paper is devoted to the estimation of the state of such a model from discretetime measurements. Measurement times t 1 < t 2 < • • • < t k max may not be regularly spaced. The step from t k to t k+1 , with size h k = t k+1 -t k , is referred to as the (k + 1)st step. The vector y(t k ) ∈ R m of the values measured at time t k is assumed to satisfy

y(t k ) = y(t k ) + δ (t k ), (2)
where δ (t k ) is an additive measurement error vector and

y(t k ) = h (z (t k ) ,t k ) (3)
is the result that would have been obtained from ideal measurements. Usually, n < m. The absolute value of the measurement error is assumed to be bounded, with a known upper bound δ , which implies that δ (t k ) ∈ [-δ ; δ ] for all k. As a result, we have y(t k ) ∈ [y(t k )] = [y(t k )δ , y(t k ) + δ ] .

If there are uncertain parameters in the measurement equation, then they can be incorporated in the extended state vector, just as we have done with the uncertain parameters in the state equation. In what follows, we shall therefore only be concerned with the computation of guaranteed estimates of the (extended) state vector, which contain guaranteed estimates of the parameters of the state and observation equations, if any. This is why the extended state will simply be called state in what follows, unless we need to distinguish the parameters from the state proper.

Guaranteed nonlinear state estimation in this context of bounded-measurement errors have been addressed by a number of authors, see e.g. [START_REF] Chernousko | Optimal guaranteed estimates of indeterminacies with the aid of ellipsoids[END_REF][START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF][START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][START_REF] Kieffer | Guaranteed parameter estimation for cooperative systems[END_REF][START_REF] Kühn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF][START_REF] Meslem | Interval observers for uncertain nonlinear systems. application to bioreactors[END_REF][START_REF] Moisan | Near optimal interval observers bundle for uncertain bioreactors[END_REF][START_REF] Raissi | Set membership state and parameter estimation for systems described by nonlinear differential equations[END_REF][START_REF] Raissi | Interval observer design for consistency checks of nonlinear continuous-time systems[END_REF][START_REF] Veres | Geometric bounding toolbox[END_REF][START_REF] Walter | Guaranteed nonlinear parameter estimation in knowledge-based models[END_REF]. All these approaches enclose the set of all state values consistent with the model, and the measurements and noise bounds. They differ first by the wrappers used to perform enclosure. Ellipsoids are used in [START_REF] Chernousko | Optimal guaranteed estimates of indeterminacies with the aid of ellipsoids[END_REF][START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF][START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF][START_REF] Veres | Geometric bounding toolbox[END_REF], zonotopes in [START_REF] Kühn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF], boxes in [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF][START_REF] Kieffer | Guaranteed parameter estimation for cooperative systems[END_REF][START_REF] Meslem | Interval observers for uncertain nonlinear systems. application to bioreactors[END_REF][START_REF] Moisan | Near optimal interval observers bundle for uncertain bioreactors[END_REF][START_REF] Raissi | Set membership state and parameter estimation for systems described by nonlinear differential equations[END_REF][START_REF] Raissi | Interval observer design for consistency checks of nonlinear continuous-time systems[END_REF], and union of boxes in [START_REF] Walter | Guaranteed nonlinear parameter estimation in knowledge-based models[END_REF]. A second important difference is in the hypotheses about the measurements. In [START_REF] Meslem | Interval observers for uncertain nonlinear systems. application to bioreactors[END_REF][START_REF] Moisan | Near optimal interval observers bundle for uncertain bioreactors[END_REF][START_REF] Raissi | Interval observer design for consistency checks of nonlinear continuous-time systems[END_REF], continuous-time measurements are assumed to be available. This rather unrealistic assumption allows nice convergence properties of the estimators to be obtained.

Here, discrete-time measurements will be considered, which appears more realistic, but makes convergence analysis much more difficult, especially for nonlinear systems. Techniques using boxes or union of boxes as wrappers usually rely on interval analysis and guaranteed integration of ODEs.

We shall assume in this chapter that bounds are available on the possible value of the initial state, and that estimates of the present state are to be produced based on the past measurements only. As for most state estimators including the celebrated Kalman filter, the (k + 1)-st step of the procedure will then consist of two steps: a prediction step that predicts the evolution of the state between two instants of time at which measurements are obtained, and a correction step during which the newly acquired data are taken into account to reduce the uncertainty in the result of the prediction step. In the context of guaranteed estimation, the prediction step involves guaranteed integration under consideration of all uncertainties, and the correction step must eliminate any part of the predicted set that can be proved to be inconsistent with the new measurements given the bounds on the measurement errors. For the sake of simplicity, no state perturbation has been considered here. On how state perturbations may be considered, see, e.g., [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF].

Section 2 overviews bounded-error state estimation in an idealized context. Section 3 considers two approaches that rely on interval analysis and try to counteract the pessimism of guaranteed integration when the initial conditions are uncertain, as is the case in an estimation context. The first of these approaches is based on Müller's theorem and presented in Section 3.1, while the second, based on the use of high-order Taylor models is described in Section 3.2. Various ways to implement correction steps are then described in Section 4. The resulting algorithms are compared in Section 5 on a simple test-case.

Idealized State Estimation

Let Z(t) be the set of all state values that are consistent with the information available up to time t. An idealized bounded-error counterpart of the Kalman filter for nonlinear discrete-time systems, alternating prediction and correction steps [START_REF] Jaulin | Applied Interval Analysis[END_REF], may be considered to build Z(t) at each t.

If z(t k ) at t k is only known to belong to Z(t k ), the set of solutions at time

t > t k of (1) is z(t;t k , Z(t k )) = {z(t;t k , z(t k )) | z(t k ) ∈ Z(t k )} .
For the (k + 1)-st prediction step, one has thus to compute the predicted set

Z + (t k+1 ) = z(t k+1 ;t k , Z(t k )).
For the (k + 1)-st correction step, one has to take into account the measurement available at time t k+1 to update Z + (t k+1 ) and obtain

Z(t k+1 ) = z ∈ Z + (t k+1 ) | h (z,t k+1 ) ∈ [y(t k+1 ) -δ , y(t k+1 ) + δ ] , (4)
thus

Z(t k+1 ) ⊆ Z + (t k+1 ).
Provided that all hypotheses on the state equation and the measurement noise are satisfied, Z(t k+1 ) does contain z (t k+1 ). The main difficulty in this idealized algorithm comes from the fact that Z(t k ) may have a quite complex shape. Outerapproximations of the sets Z and Z + have thus to be evaluated. These outerapproximations may consist of interval vectors (boxes), union of non-overlapping boxes (subpavings), or may be described using Taylor models. Implementable counterparts of the idealized prediction and correction steps are now described.

Prediction Step

A naive approach to obtain an outer-approximation for z(t k ;t 0 , [z(0)]), k = 1 . . . k max , would be to use one of the guaranteed ODE solvers based on interval analysis AWA [START_REF] Lohner | Computation of guaranteed enclosures for the solutions of ordinary initial and boundary value-problem[END_REF], VNODE [START_REF] Nedialkov | Methods for initial value problems for ordinary differential equations[END_REF], COSY IV [START_REF] Hoefkens | Efficient high-order methods for ODEs and DAEs[END_REF], VSPODE [START_REF] Lin | Validated solution of ODEs with parametric uncertainties[END_REF], or ValEncIA-IVP [START_REF] Rauh | Extensions of ValEncIA-IVP for reduction of overestimation, for simulation of differential algebraic systems, and for dynamical optimization[END_REF][START_REF] Rauh | VALENCIA-IVP: A comparison with other initial value problem solvers[END_REF]. The main difficulty is to obtain accurate enclosures for the solutions, when there are uncertain parameters, bounded state perturbation, or uncertain initial conditions.

One may enclose the solutions of (1) with uncertain initial conditions between a pair of coupled system of ODEs with deterministic initial conditions using Müller's theorem [START_REF] Müller | Über das Fundamentaltheorem in der Theorie der gewöhnlichen Differentialgleichungen[END_REF]. Other types of uncertainty may be taken into account as well, such as unknown but bounded inputs. Any guaranteed tool for solving ODEs may then be used to solve this system, see Section 3.1.

When only the initial conditions are undetermined, but known to belong to some box, guaranteed ODE solvers such as COSY IV, VSPODE, or ValEncIA-IVP perform quite well, since they are evaluating a Taylor development of the solution with interval remainder, this developments being made also with respect to the initial condition, see Section 3.2.

Using M üller's Theorem

Here, the solution is obtained by bounding the solutions of dynamical systems with uncertain parameters or initial conditions using deterministic dynamical systems. This approach has been previously presented in [START_REF] Kieffer | Interval analysis for guaranteed nonlinear parameter and state estimation[END_REF][START_REF] Walter | Interval analysis for guaranteed nonlinear estimation[END_REF] in the context of cooperative dynamical models, i.e., models such as [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF] for which the off-diagonal terms of the Jacobian matrix of f are positive. These results were inspired by the interval observer proposed in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF]. Müller's theorems [START_REF] Müller | Über das Fundamentaltheorem in der Theorie der gewöhnlichen Differentialgleichungen[END_REF], which have recently been used in the context of guaranteed simulation [START_REF] Gennat | Simulation of uncertain systems with guaranteed bounds[END_REF], make it possible to bound the solutions of more general dynamical models, see also [START_REF] Kieffer | Guaranteed estimation of the parameters of nonlinear continuoustime models: contributions of interval analysis[END_REF].

M üller's Theorem

The following theorem [START_REF] Müller | Über das Fundamentaltheorem in der Theorie der gewöhnlichen Differentialgleichungen[END_REF] may be directly applied to bound the solutions of dynamical models such as [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF]

D ± ω i (t) min T i (t) f i (z) and D ± Ω i (t) max T i (t) f i (z) ,
where T i (t) is the subsets of D defined by

T i (τ) :    z i = ω i (τ) , ω j (τ) z j Ω j (τ) , j = i, t = τ,
and T i (t) is the subset of D defined by

T i (τ) :    z i = Ω i (τ) ω j (τ) z j Ω j (τ) , j = i, t = τ.
Then, for any given z (t k ) ∈ [z(t k )], a solution to (1) exists, such that

ω (t) z(t) Ω (t) ∀t ∈ [t k ,t k+1 ],
and this solution is equal to

z (t k ) at t = t k . Moreover, if for any t ∈ [t k ,t k+1 ] , f (z,t) is Lipschitz with respect to z, then for any given z (t k ) ∈ [z(t k )] this solution is unique. ♦
The main idea of this theorem is to bracket the solutions of (1) between the solution of two deterministic ODEs. The initial conditions for these ODEs are given by 1. and the conditions that have to be satisfied by each solution are given by 2., see Section 5 for an example.

Prediction Step Using M üller's Theorem

Müller's theorem allows the evaluation of lower and upper bounds for the solution of (1) provided that two functions ω (t) and Ω (t) are available. The interval function [Φ] (t) = [ω (t) , Ω (t)] can then be seen as an inclusion function for all solutions z (t) of the state equation in [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF].

[Φ] (t) provides a box containing the state at each time instant. More accurate descriptions of the predicted set may be obtained by splitting the box corresponding to the initial conditions into non-overlapping subboxes, and to apply Müller's theorem on each of the resulting subboxes. A list of overlapping boxes is obtained, which may be merged into a subpaving using the ImageSp algorithm [START_REF] Jaulin | Applied Interval Analysis[END_REF]. This subpaving may then be used by any implementable correction step described in Section 4, or used to perform a new prediction until a measurement is available. The accuracy of the description of Z + (t k ) at time t k depends on some precision parameter ε, which determines the size of the boxes obtained after splitting the initial box or subpaving used by the ImageSp algorithm.

The construction of ω (t) and Ω (t) is usually easy on a case-by-case basis, as illustrated in Section 5.1. For more complex systems, hybrid automata may be put at work to build these functions, as detailed in [START_REF] Meslem | Guaranteed state bounding estimation for uncertain non linear continuous systems using hybrid automata[END_REF].

An inclusion function for h may similarly be obtained using the sensitivity functions of the output with respect to the state, see Section 4.2.3.

Verified Integration Based on Taylor Models

Verified integration techniques such as VNODE [START_REF] Nedialkov | Methods for initial value problems for ordinary differential equations[END_REF] are based on a Taylor series expansion in time. COSY-VI [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF] performs, in addition to this expansion in time, an expansion with respect to the initial state vector, denoted by z in what follows. The box to which z is assumed to belong is given by [z]. The expansion point with respect to the initial state vector z is some ẑ ∈ [z], and the expansion point with respect to time is t k . The flow of the differential equation in a given time interval

[t k ,t k+1 ] is enclosed by a n-dimensional Taylor model T ρ (z -ẑ,t -t k ) := P ρ (z -ẑ,t -t k ) + I ρ,k+1 , with z ∈ [z] and t ∈ [t k ,t k+1 ] ,
where P ρ (zẑ,tt k ) is the multivariate polynomial part of order ρ and I ρ,k+1 is the remainder box. The i-th entry of the n-dimensional vector

T ρ (z -ẑ,t -t k ) is denoted by T ρ,i (z -ẑ,t -t k ).
The Taylor model at t = t k+1 is written as

T ρ,k+1 (z -ẑ) := P ρ,k+1 (z -ẑ) + I ρ,k+1 . The i-th entry of T ρ,k+1 (z -ẑ) is given by T ρ,i,k+1 (z -ẑ).
Verified integration methods that use a single box or a single parallelepiped for the state enclosure may suffer from large overestimation, especially for nonlinear systems. The flow representation by Taylor models makes it possible to obtain tight enclosures of non-convex sets and leads to less overestimation.

The integral form of the differential equation ( 1) is given by

O (z(t)) := z(t k ) + t t k f(z(t ),t )dt
Applying O to a Taylor model for the integration in the time-interval [t k ,t k+1 ] yields

O(P ρ (z -ẑ,t -t k ) + I ρ,k+1 ) = z(t k ) + t t k f(P ρ (z -ẑ,t -t k ) + I ρ,k+1 )dt ,
where z(t k ) is represented by its Taylor model enclosure at t = t k .

T ρ,k = P ρ,k (z -ẑ) + I ρ,k .
This leads to

O(P ρ (z -ẑ,t -t k ) + I ρ,k+1 ) = P ρ,k (z -ẑ) + I ρ,k + t t k f(P ρ (z -ẑ,t -t k ) + I ρ,k+1 )dt .
The goal for the integration from t k to t k+1 consists in determining a Taylor model

T ρ (z -ẑ,t -t k ) such that O(P ρ (z -ẑ,t -t k ) + I ρ,k+1 ) ⊂ P ρ (z -ẑ,t -t k ) + I ρ,k+1 ∀z ∈ [z] and ∀t ∈ [t k ,t k+1 ].
The polynomial part and the interval remainder are determined in separate steps. A detailed description of these steps is given in [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF][START_REF] Kletting | Verified Methods for State and Parameter Estimators for Nonlinear Uncertain Systems with Applications in Engineering[END_REF].

For numerical and implementation reasons it is advantageous to have the unit box [-1; 1] n as the domain box in each integration step [START_REF] Makino | Suppression of the wrapping effect by Taylor model-based verified integrators: Long-term stabilization by shrink wrapping[END_REF]. Thus the initial enclosure [z(0)] of the extended state vector z(t) is expressed as a Taylor model according to

[z(0)] = T(z) = c + D z with z i ∈ [-1; 1], i = 1 . . . n, where c is the midpoint of [z(0)] and D is a diagonal matrix with d i,i = rad([z(0)]).
The expansion in initial state reduces the overestimation that may occur during integration. To limit the long-term growth of the remainder error and to reduce overestimation the following strategies can be applied:

• Shrink Wrapping: the interval remainder is absorbed in the polynomial part [START_REF] Makino | Suppression of the wrapping effect by Taylor model-based verified integrators: Long-term stabilization by shrink wrapping[END_REF].

• Preconditioning: the ODE solution is studied in a more suitable coordinate system [START_REF] Makino | Suppression of the wrapping effect by Taylor model-based verified integrators: Long-term stabilization by preconditioning[END_REF]. • Splitting: the domain box [z] is split into subboxes and the enclosure of z(t) is given by a list of Taylor models [START_REF] Kletting | Verified Methods for State and Parameter Estimators for Nonlinear Uncertain Systems with Applications in Engineering[END_REF]. This is described in the following section.

Splitting the Domain Box

As when using Müller's theorem, splitting of the domain box into subboxes [START_REF] Kletting | Verified Methods for State and Parameter Estimators for Nonlinear Uncertain Systems with Applications in Engineering[END_REF] may be useful to reduce overestimation. The state vector z k+1 at t = t k+1 is then enclosed by a list T k+1 of Taylor models

z k+1 ∈ T k+1 = T (1) ρ,k+1 (z), T (2)
ρ,k+1 (z) . . .

T (L k+1 ) ρ,k+1 (z) with z i = [-1; 1], i = 1 . . . n and L k+1 ≤ L max .
where L max is the maximum allowed number of Taylor models. Consider a Taylor model T ρ,k+1 (z) with the domain box

[z], z ∈ [z] . The domain box of this Taylor model is split into subboxes [z (l) ], l = 1 . . . L, L l=1 [z (l) ] = [z] .
To obtain again the unit box as a domain box, [z (l) ] is expressed as a Taylor model according to

[z (l) ] = T (l) (z) = c (l) + D (l) z with z i ∈ [-1; 1] , i = 1 . . . n ,
where c (l) is the midpoint of [z (l) ] and D (l) is a diagonal matrix with d

(l) i,i = rad([z (l) i ]
). The components of the original initial state vector z of T ρ,k+1 (z) are replaced by the components of T (l) (z) by substituting T (l) i (z) for z i , which results in a modified Taylor model

T (l) ρ,k+1 (z) = T ρ,k+1 ( T (l) (z))
for each subbox [z (l) ].

To determine the component in which the domain box has to be split, splitting criteria have to be evaluated for the considered Taylor model T ρ,k+1 (z). Splitting is carried out perpendicularly to the direction which has been calculated by the spliting criteria. Approaches to determine the splitting direction are described in [START_REF] Kletting | Verified Methods for State and Parameter Estimators for Nonlinear Uncertain Systems with Applications in Engineering[END_REF].

If several Taylor models are already present, the most appropriate Taylor model for the splitting has to be selected. This is done by calculating the interval enclosure of each Taylor model and the corresponding pseudo volume of the resulting box. The pseudo volume of a n-dimensional interval vector is calculated by the multiplication of the interval diameters of all its components. The Taylor model with the largest pseudo volume is selected. Alternatively, the Taylor model with the largest interval remainder could be selected.

How splitting of the domain box is combined with preconditioning is described in detail in [START_REF] Kletting | Verified Methods for State and Parameter Estimators for Nonlinear Uncertain Systems with Applications in Engineering[END_REF].

Prediction Step Using Taylor Models

At time t k , the extended state vector is enclosed by a list T k of Taylor models

z k ∈T k = T (1) ρ,k (z), T (2) ρ,k (z) . . . T (L k ) ρ,k (z) with z i = [-1; 1], i = 1 . . . n and L k ≤ L max .
First a Taylor model is selected for splitting, then a splitting criterion is evaluated, and the Taylor model is split by splitting the domain box in subboxes. Taylor models are split until a pre-specified number of Taylor models or number of splittings is reached. Next, for each Taylor model a verified integration is performed.

The resulting enclosure of the extended state vector at time t k+1 (after the prediction step) is then given by

T pr k+1 = T (pr,1) ρ,k+1 (z), T (pr,2) ρ,k+1 (z) . . . T (pr,L pr k+1 ) ρ,k+1 (z)
with

z i = [-1; 1], i = 1 . . . n and L pr k+1 ≤ L max .
The prediction step is repeated until measurements become available. The result of the last prediction step before measurements become available is used as an initial enclosure of the next correction step.

Correction Step

Assume that the prediction step at time t k+1 has produced a set Z + (t k+1 ) that contains z(t k+1 ). This set may consist of a single box, a list of potentially overlapping boxes, or may be a subpaving.

The measurement vector y(t k+1 ) obtained at time t k+1 has now to be taken into account. Several practical implementations of the idealized correction step (4) are now presented.

Using Set Inversion Via Interval Analysis

The aim of the Set Inverter Via Interval Analysis (SIVIA) algorithm [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF] is to eliminate parts of Z + (t k+1 ) that can be proved to be inconsistent with the measurements, the measurement equation and the bounds on the measurement noise.

Consider a box

[z] ⊂ Z + (t k+1 ). If h ([z],t k+1 ) ⊂ [y(t k+1 )-δ , y(t k+1
)+δ ], then all z ∈ [z] are consistent with the measurements, model, and noise bounds. Therefore, [z] has been proved to belong to Z(t k+1 ). If h ([z],t k+1 ) ∩ [y(t k+1 )δ , y(t k+1 ) + δ ] = / 0, then no z ∈ [z] is consistent with the measurement, model and noise bounds.

Thus, [z] has an empty intersection with Z(t k+1 ) and can be rejected. If none of the previous conditions is satisfied, [z] is said undetermined, as parts of [z] may belong to Z(t k+1 ).

The SIVIA algorithms iteratively performs selection, elimination, or bisection of boxes, starting from a large initial search box, list of boxes or subpaving. Undetermined boxes are bisected until their width is smaller than some precision parameter ε, which helps to trade-off complexity and accuracy of representation of the solution set, Z(t k+1 ) here. The solution provided by SIVIA is a subpaving, consisting of inner boxes and undetermined boxes deemed too small to be further bisected. This subpaving may be fed to the next prediction step. See [START_REF] Jaulin | Applied Interval Analysis[END_REF] for more details and implementation issues.

Using Contractors

Bounded-error measurements translate into vector inequality constraint k (z) 0, to be understood componentwise. A contractor C k for z is an algorithm to compute a box

C k ([z]) such that C k ([z]) ⊂ [z] , {z ∈ [z] |k (z) 0} ⊂ C k ([z]) . (5)
The first relation in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF] ensures that [z] is contracted, while the second guarantees that no value of z satisfying the constraints is lost. Contractors can be similarly defined in the case of equality constraints.

Improving the State Estimate

Given the measurement equations ( 2), (3) and the fact that δ k ∈ -δ , δ , two constraints may be obtained that have to be satisfied by the state vector at time t k , namely

k 1 (z) = y (t k ) -h (z,t k ) + δ 0 and k 2 (z) = -y (t k ) + h (z,t k ) + δ 0.
Various types of contractors may be considered [START_REF] Jaulin | Applied Interval Analysis[END_REF], depending on the structure of h (z,t k ). If m = n, the interval Newton or the Krawczyk contractors may be employed.

Improving the Initial Conditions Estimate

If one is interested in obtaining a better estimate of the initial conditions, one may write h (z,t k ) as a function of z. For that purpose, consider the i-th entry of h (z,t k ). We have

h i ([z k+1 ] ,t k ) ⊂ h i ([z k+1 ] ( z) ,t k ) + ([z] -z) T ∂ h i (z,t k ) ∂ z [z k+1 ] , (6)
where z ∈ [z] and [z k+1 ] ( z) is the box obtained when integrating the system with known initial conditions taken as z. In ( 6),

∂ h i (z,t k ) ∂ z =     ∂ z 1 ∂ z 1 . . . ∂ z n ∂ z 1 . . . . . . . . . ∂ z 1 ∂ z n . . . ∂ z n ∂ z n         ∂ h i ∂ z 1 . . . ∂ h i ∂ z n    
where ∂ z i ∂ z j is the sensitivity function of the i-th state component with respect to the initial condition of the j-th state component.

Taking into account the measurement at time t k ,

h i ([z k+1 ] ( z) ,t k ) + ([z] -z) T ∂ h i (z,t k ) ∂ z [z k+1 ] ⊂ y i (t k ) -δ , y i (t k ) + δ , i = 1 . . . m.
Thus each [z i ], i = 1 . . . n, has to satisfy

[z i ] ⊂ y (t k ) --δ , δ -h i ([z k+1 ] ( z) ,t k ) -∑ j =i ([z j ] -z j ) ∂ h i ∂ z j ([z k+1 ]) / ∂ h i ∂ z i ([z k+1 ]) + z i , leading to a contracted box [z] new = C i ([z]
), the components of which are defined as

[z i ] new = [z i ] ∩ y (t k ) --δ , δ -h i ([z k+1 ] ( z) ,t k ) -∑ j =i ([z j ] -z j ) ∂ h i ∂ z j ([z k+1 ]) / ∂ h i ∂ z i ([z k+1 ]) + z i for i = 1 . . . n.
This contractor requires the computation of an inclusion function for the sensitivity functions of each state component with respect to the initial condition; see the next section.

Sensitivity Functions

Denote the first-order sensitivity of z j with respect to z k by s jk

s jk (z,t) = ∂ z j ∂ z k (z,t) .
Assume for simplicity that the model output is linear in the state and it is given by

h (z (t) ,t) = Cz (t) ,
where

C = (c 1 . . . c n ) T is a known matrix. Then ∂ h (z (t) ,t) ∂ z =     c T 1 ∂ z ∂ z 1 . . . c T n ∂ z ∂ z n     .
Differentiate once the jth row of (1) with respect to z k to obtain the differential equation

ṡ jk = ∂ f j (z) ∂ z j s jk . ( 7
)
Since z is a scaled version of z (t 0 ), the initial conditions for the sensitivity functions are

s jk (t 0 ) = ∂ z j (t 0 ) ∂ z k = d j j if j = k, 0 else.
The sensitivity functions may then be obtained by considering a new extended statespace model consisting of ( 1) and of all differential equations (7) satisfied by the sensitivity functions. Müller's theorem turns out to be especially useful, as the extended state-space model is seldom cooperative, even if this is the case of the initial state-space model.

Using Taylor Models

In the case of Taylor models, the correction step is quite different from what has been described in Section 4.1. The measurement equation ( 2), ( 3) is rewritten as

h (z k+1 ) + δ -y(t k+1 ) = 0. ( 8
)
Each Taylor model T

(pr,l) ρ,k+1 (z) of T (pr)
k+1 is now considered separately and substituted for z k+1 in (8) to obtain

h T (pr,l) ρ,k+1 (z) + δ -y(t k+1 ) = 0. ( 9
)
For nonlinear systems, the left hand side of ( 9) is nonlinear even for a linear measurement equation, since the Taylor model T (pr,l) ρ,k+1 (z) is nonlinear. One knows that δ ∈ [-δ , δ ], thus, (9) may be solved for z with an interval Newton method, which leads to a tightened domain interval [z], hence z ∈ [z]. Here, the Krawczyk method is used. For n > m, (9) is under-determined and cannot be inverted. A simple approach is to solve [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF] only for m variables (components of z), while considering the remaining nm variables as constant intervals.

This procedure is illustrated for a very simple linear example with n = 2 and m = 1. Consider the case when ( 9) is given by

z 1 + 0.5z 2 + δ 1 -y 1 (t k+1 ) = 0, ( 10
)
with 10) is first solved for z 1 , the interval [-1; 1] being used for z 2 . We have

[z 1 ] = [-1; 1], [z 2 ] = [-1; 1], y 1 (t k+1 ) = 0.9, and [δ 1 ] = [-0.1; 0.1]. Now, (
[z 1 ] = y 1 (t k+1 ) -[δ 1 ] -0.5[z 2 ] = 0.9 -[-0.1; 0.1] -0.5[-1; 1] = [0.3; 1.5].
Now this result is intersected with the initial interval enclosure [-1; 1] resulting in

[z 1 ] = [0.3; 1.5] ∩ [-1; 1] = [0.3; 1].
Next, ( 10) is solved for z 2 with the new

[z 1 ] [z 2 ] = 2(y 1 (t k+1 ) -[δ 1 ] -[z 1 ]) = 2(0.9 -[-0.1; 0.1] -[0.3; 1]) = [-0.4; 1.4].
An intersection with the initial interval enclosure [-1; 1] results in

[z 2 ] = [-0.4; 1.4] ∩ [-1; 1] = [-0.4; 1].
Another possibility is to consider a sufficient number of previous measurements y(t ≤ t k ) and the corresponding Taylor models of the right hand side of (9) to obtain the missing nm equations.

If no solution in [z] can be found, then the corresponding Taylor model is inconsistent with the data and can be deleted.

As previously, the Taylor model is computed for a normalized domain box and written as

[z] ∈ T(z) = c + D z with z i ∈ [-1; 1], i = 1 . . . n,
where c is the midpoint of [z] and D is a diagonal matrix with di,i = rad([z i ]). The Taylor model T (c,l) k+1 after the correction step is then given by

T (c,l) ρ,k+1 (z) = T (pr,l) ρ,k+1 ( T(z)),
which defines the l-th Taylor model for the next prediction step:

T (l) ρ,k+1 (z) = T (c,l) ρ,k+1 (z).
The list of Taylor models at k + 1 is then given by

z k+1 ∈T ρ,k+1 = T (1) ρ,k+1 (z), T (2)
ρ,k+1 (z) . . .

T (L k+1 ) ρ,k+1 (z) with z i = [-1; 1], i = 1 . . . n with L k+1 ≤ L pr k+1 ≤ L max .
Interval Newton methods, like the Krawczyk method, involve the computation of the derivatives of h with respect to the initial states z similar to [START_REF] Hoefkens | Efficient high-order methods for ODEs and DAEs[END_REF] in the correction step for the method based on Müller's theorem. Since the Taylor model performs an expansion in the initial states z, the derivatives are easily obtained by calculating them for the Taylor model obtained from [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF].

The correction step in combination with preconditioning is explained in [START_REF] Kletting | Verified Methods for State and Parameter Estimators for Nonlinear Uncertain Systems with Applications in Engineering[END_REF], together with a correction step that involves consistency tests.

Simulation Results

For a comparison of the two approaches, consider the following nonlinear system

   ẋ1 = -p 3 x 1 - p 1 x 1 1 + p 2 x 1 + p 4 x 2 ẋ2 =p 3 x 1 -p 4 x 2
with the initial value of the state x 1 (0) = 1 and x 2 (0) = 0. Only x 2 is measured and the measurement equation at t = t k is

y k = x 2 (t k ) + δ k .
The parameters p 1 , p 2 and p 3 are given by p 1 = 1, p 2 = 1.2 and p 3 = 0.5. The parameter p 4 is uncertain with p 4 ∈ [0.1, 0.5]. Thus, the extended state satisfies

       ż1 = -p 3 z 1 - p 1 z 1 1 + p 2 z 1 + z 3 z 2 ż2 =p 3 z 1 -z 3 z 2 ż3 =0 (11)
with z 1 = x 1 , z 2 = x 2 , and z 3 = p 4 . Measurements are assumed to take place every 2 s. The measurement uncertainty is given by δ (t k ) ∈ [-0.005, 0.005] for all t k .

Results with Prediction Based on M üller's Theorem

Correction using Sivia

We used the Sivia algorithm combined with ImageSp to estimate the extended state vector. For this purpose, a coupled system of ODEs is built from ( 11)

                           ż1 = -p 3 z 1 - p 1 z 1 1 + p 2 z 1 + z 3 z 2 ż2 =p 3 z 1 -z3 z 2 ż3 =0 ż1 = -p 3 z1 - p 1 z1 1 + p 2 z1 + z3 z2 ż2 =p 3 z1 -z 3 z2 ż3 =0 (12)
with z 1 (0) = z1 (0) = 1, z 2 (0) = z2 (0) = 0, and z 3 (0) , z3 (0) dependent on the iteration in Sivia. Using [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF], an inclusion function for h (z

(t k )) = z 2 (t k ) is ob- tained.
The estimates for p 4 obtained at t = 10 s using the Sivia algorithm for various values of the precision parameter ε are provided in Table 1. As expected, reducing ε increases the computing time and the accuracy of the estimate. These results have been obtained using a single core of a two-processor, quadri-core Intel Xeon CPU E5462 at 2.80 GHz with 6 144 KB cache and 64 GB RAM.

Correction by Constraint Propagation

Since only the third initial condition is unknown, only the first-order sensitivity functions with respect to z 3 have to be evaluated. For that purpose, each ODE in [START_REF] Kieffer | Interval analysis for guaranteed nonlinear parameter and state estimation[END_REF] is derived with respect to z 3 to get

           ṡ13 = -p 3 - p 1 1 + p 2 z 1 + p 1 p 2 (1 + p 2 z 1 ) 2 s 13 + s 33 z 2 + z 3 s 23 ṡ23 =p 3 s 13 -s 33 z 2 -z 3 s 23 ṡ33 =0 (13)
with s 13 (0) = s 23 (0) = 0 and s 33 (0) = 1. Müller's theorem may be used with [START_REF] Kieffer | Guaranteed estimation of the parameters of nonlinear continuoustime models: contributions of interval analysis[END_REF] to compute an enclosure of the sensitivity function of the state with respect to the third initial condition, which is also

p 4                                  ṡ13 = -p 3 - p 1 1 + p 2 z 1 + p 1 p 2 z 1 (1 + p 2 z 1 ) 2 s 13 + z 2 s 33 + z 3 s 23 ṡ23 =p 3 s 13 -s 33 z 2 -z 3 s 23 ṡ33 =0 s13 = -p 3 - p 1 1 + p 2 z 1 + p 1 p 2 z 1 1 + p 2 z 1 2 s13 + s33 z 2 + z 3 s23 s23 =p 3 s13 -s 33 z 2 -z 3 s 23 s33 =0
This enclosure uses the fact that all quantities are positive, except s 23 (t) ≤ 0 for t ≥ 0. The contractor described in Section 4.2 may then be employed in conjunction with Sivia to reduce the uncertainty on the initial value of the third component of the state. Table 2 describes the evolution with time of the estimate for p 4 as a function of the precision parameter ε. Reducing ε again increases the accuracy at which the estimate is obtained, but the price to be paid is an increased complexity, since more ODEs are solved. Contrary to Sivia, a decent estimate is obtained even with the largest value of ε. The first and second measurements provide the most information about p 4 , since for these measurements, the best decrease in the size of [z 3 ] is observed. See also Figure 1. The same processor as in Section 5.1.1 has been used here.

ε = 0.1, =21 s t ε = 0.05, =42 s t ε = 0.025, =83 s t ε = 0.1, =21 s t ε = 0.05, =42 s t ε = 0.025, =83 s t ε = 0.1, =21 s t ε = 0.05, =42 s t ε = 0.025, =83 s t (c) Estimation of p 4
Fig. 1 Comparison of the estimation results with Sivia for different values of the precision parameter ε used in the Müller method

Prediction and Correction Involving Taylor Models

Estimation results for the approach based on Taylor models are shown for different orders and numbers L max of Taylor models. In Figure 2 results for ρ = 4, ρ = 6, and ρ = 12 for a single Taylor model (L max = 1, hence without splitting of the domain box) are depicted together with one result for ρ = 5 with 4 Taylor models (L max = 4), hence with splitting of the domain box. In this section, all computations have been done on a Intel Centrino Core2 Duo T7300 at 2 GHz. In Table 3 the interval enclosures for the evolution of the estimate for p 4 as a function of time and total computing time are given. In Table 4 the interval enclosures for the estimated parameter p 4 at t = 10 s for various orders and numbers of Taylor models are presented. Increasing the order from 4 to 6 leads to an improvement with a slightly higher computation time. However if the order is increased to 12, the computation time increases drastically (from 1.12s to 4.83s) without tightening the enclosures in any significant manner. Further improvement of the estimation results is achieved only if splitting is applied as the results for ρ = 5 and L max = 4 indicate. Even the computation time is lower than for ρ = 12 and L max = 1.

Conclusions and Perspectives

Physics, chemistry (and most other experimental sciences) tend to produce continuoustime models whose outputs are nonlinear in their unknown parameters. When these models are in state-space form, they are almost always nonlinear in the extended state vector obtained by concatenating the state and parameter vectors. Most of the methods available for estimating this extended state vector are based on linearization or random exploration and cannot provide any guarantee as to their results. Since each of the parameters and state variables of such models has a concrete physical meaning, this is unfortunate. One would rather like to be able to characterize in some guaranteed way the set of all acceptable estimates given what is known of the uncertainty in the experimental data.

For quite some time, this has been completely out of reach, but a combination of advances in interval analysis, constraint propagation, and guaranteed integration of ODEs, together with a massive increase in computing power have now made this achievable, at least for small-scale models. This paper has presented, in a coordinated manner, a variety of tools that can be used when bounds are available on the acceptable difference between the data and corresponding model output.

Guaranteed set estimation makes it possible to bypass the usual requirements of identifiability, observability, and persistency of excitation. These properties should nevertheless definitely contribute to the quality of the parameter and state estimates. More generally, the problem of experiment design for guaranteed parameter or state estimation is an interesting but still largely open question.

The main challenge is to increase the complexity of the models that can be studied with this type of guaranteed approach. To this end, it is necessary to use and possibly combine tools that struggle as efficiently as possible with the pessimism inherent to interval analysis and guaranteed integration and the curse of dimensionality. Contractors, which make it possible to eliminate parts of the search region without the need for bisections, and high-order Taylor models are among the most promising avenues for research.

1 ρ

 1 = 4, 1 T. mod., time = 0.87 s ρ = 6, 1 T. mod., time = 1.12 s ρ = 12, 1 T. mod., time = 4.83 s ρ = 5, 4 T. mod., time = 3.51 s (a) Estimation of x 1

Fig. 2

 2 Fig. 2 Comparison of the estimation results for different orders and number of Taylor models.

  in the presence of uncertain initial conditions z (t k ) ∈ [z(t k )], where [z(t k )] is a box in state space. (t) and Ω i (t) , i = 1 . . . n x , are continuous on [t k ,t k+1 ] and satisfy 1. ω (t k ) = z (t k ) and Ω (t k ) = z (t k ), 2. the left derivatives D -ω i (t) and D -Ω i (t) and right derivatives D + ω i (t) and D + Ω i (t) of ω i (t) and Ω i (t) satisfy, for i = 1 . . . n and all t ∈ [t k ,t k+1 ] ,

	Theorem 1. Assume that f (z(t)) in (1) is continuous on
	T :	ω (t) x Ω (t) t k t t k+1
	where ω i

Table 1

 1 Estimates for p 4 obtained with Sivia at t = 10 s for various values of ε

	ε	time (s) estimate for p 4
	0.05	1.28	[0.125,0.35]
	0.025	1.85 [0.1625,0.3125]
	0.01	4.3 0.2[1718,8750]
	0.005	8.28 0.2[3303,6875]
	0.0025 15.77 0.2[4063,6251]
	0.001 31.95 0.2[4546,5804]

Table 2

 2 Evolution of the estimate for p 4 with Sivia for various values of the precision parameter ε as a function of the time and total computing time

	Time	ε = 0.2	ε = 0.1	ε = 0.05	ε = 0.01
	2 s	0.2[20488,57674] 0.2[24304,57668] 0.2[23906,57668] 0.2[24323,56828]
	4 s	0.2[30813,57674] 0.2[34575,57668] 0.2[35411,57668] 0.2[40556,56828]
	6 s	0.2[30813,57674] 0.2[36483,57668] 0.2[37597,57668] 0.2[44613,56828]
	8 s	0.2[30813,57674] 0.2[36483,57668] 0.2[37597,57668] 0.2[44613,56828]
	10 s 0.2[30813,57674] 0.2[36483,57668] 0.2[37597,57668] 0.2[44613,56250]
	Comp.	1.0 s	2.13 s	3.47 s	17 s

Table 3

 3 Evolution of the estimate for p 4 for various values of ρ and L max as a function of time, and total computing time

	max = 4

Time t ρ = 4, L max = 1 ρ = 6, L max = 1 ρ = 12, L max = 1 ρ = 5, L

Table 4

 4 Estimates obtained with Talyor models at t = 10 s for various values of ρ and numbers of Taylor models.


Nesrine Changuel
email: nesrine.changuel@alcatel-lucent.fr

Bessem Sayadi
email: bessem.sayadi@alcatel-lucent.frmichel.kieffer@lss-supelec.fr

Michel Kieffer

Statistical multiplexing of distributed video streams

Nous considérons un élément du réseau sensible au contenu multimédia fourni par plusieurs serveurs. Son rôle est de stocker le contenu vidéo encodé et le multiplexer pour être diffusé sur une liaison sans fil. Nous proposons une technique de contrôle décentralisée capable de satisfaire une contrainte d'équité en qualité vidéo entre les programmes. Contrairement à la plupart des systèmes de multiplexage statistique, notre système est en partie décentralisé. L'allocation de bande passante entre les programmes est centralisée et faite dans le réseau en prenant en compte la contrainte de l'équité en qualité. Chaque serveur vidéo est contrôlé indépendamment des autres, ne nécessitant aucun échange avec les autres serveurs. Le système de contrôle de l'allocation de bande renvoie aux serveurs vidéo le niveau de remplissage de la mémoire tampon à laquelle il est associé. Les encodeurs utilisent ces retours d'information pour adapter leurs compromis débit-distorsion pour que la mémoire tampon atteigne un certain niveau de remplissage référence. Les résultats expérimentaux montrent que, dans le cas de sources Gaussiennes, un équilibre est atteint, et que la contrainte de l'équité de la qualité est satisfaite.

Introduction

In multimedia broadcasting systems, multiple video services are transmitted over a fixed bandwidth channel. Therefore, the data rates required by these multimedia services are expected to vary over a wide range, depending on their multimedia content.

Statistical multiplexing of compressed video signals consists in adapting jointly the rate and the distortion of the encoded video to use the available bandwidth in an optimal way. Usually, Variable Bit Rate (VBR) [START_REF] Chengsheng | An efficient two-pass VBR encoding algorithm for H.264[END_REF] video coders are used in order to satisfy several constraints linked to the quality of the encoded streams [START_REF] Tagliasacchi | Minimum variance optimal rate allocation for multiplexed H.264/AVC bitstreams[END_REF][START_REF] Changuel | Predictive control for efficient statistical multiplexing of digital video programs[END_REF] and to the transmission delay [START_REF] Rezaei | Implementing statistical multiplexing in DVB-H[END_REF][START_REF] Changuel | Joint encoder and buffer control for statistical multiplexing of multimedia contents[END_REF]. The video quality fairness is an essential constraint that should be considered in the statistical multiplexing system in order to provide programs with quality levels of the same order of magnitude.

When considering a statistical multiplexer of H.264/AVCencoded video, the rate-distortion (R-D) control may be performed by adjusting dynamically the encoding param-eters of the source coders as in [START_REF] Changuel | Joint encoder and buffer control for statistical multiplexing of multimedia contents[END_REF]. This may also be performed by re-encoding or transcoding the encoded video to match the R-D target fixed by the control scheme. With H.264/SVC, the encoded streams are organized into quality layers and the number of transmitted layers for each frame may be adjusted [START_REF] Jacobs | Statistical multiplexing using svc[END_REF][START_REF] Changuel | End-to-end stochastic scheduling of scalable video overtime-varying channels[END_REF]. In both cases, to satisfy video quality fairness among programs, it is necessary to share quality information between encoders, transcoders, or layer filtering units. Since video servers are usually located in separate places, this may introduce some communication delays. Therefore, a decentralized rate control scheme is desirable. The problem of rate control, buffer management and video quality control for video streaming system is addressed in [START_REF] Huang | A control-theoretic approach to rate control for streaming videos[END_REF][START_REF] Huang | Analysis and design of a proportional-integral rate controller for streaming videos[END_REF]. The video streaming problem is considered by jointly controlling the encoder rate and the network congestion. The control system is formulated via a feedback control system using first a Proportional and second a Proportional-Integral (PI) controller. The robustness and the stability of the proposed controltheoretic approach are derived.

In this paper, we consider a media aware network element (MANE) fed by several remote digital media servers. The MANE is in charge of buffering the encoded streams and to build a centralized bandwidth allocation among programs. It takes into account quality information provided by each video coder to try to satisfy the fairness constraint. Each video server is controlled independently from the others, requiring no exchange between servers. The MANE feeds back to each video server the level of its associated buffer to help the remote video servers adapt their R-D trade-off so that the buffer within the MANE reaches some reference level. We adopt a control-theoretic approach that performs rate allocation via a feedback control system. In the case of Gaussian sources, we show experimentally that the proposed control system allows an equilibrium state to be reached.

The paper is organized as follows. Section 2 introduces the statistical multiplexing problem and defines some notations. Section 3 presents the way we propose to solve the statistical multiplexing in a distributed manner. In Section 4 we study the equilibrium and the stability of the proposed solution.

Problem statement

We consider a typical broadcast system in which N video programs have to be encoded and transmitted in parallel. The compressed bitstreams are multiplexed and transmitted over a communication channel allowing a constant transmission rate R c as represented in Figure 1. We assume that all video contents are provided with the same frame period ∆t. At time index j, the j-th encoded frames are fed by all video coders. For each video stream i ∈ {1, . . . , N } an individual Rate-Distortion (R-D) control process is performed per frame generating an encoding rate R ij and a distortion D ij for each flow i ∈ {1, . . . , N } at time j. Each controlled frame is then stored in its corresponding buffer of the MANE, see Figure 1. The available transmission rate R c is distributed among the N programs, (α 1 R c , . . . , α N R c ), in such a way that a fair video quality is maintained.

Proposed solution

In this section, we propose a distributed solution for multi video programs transmission over a shared and fixed channel rate. The idea is to regulate the video flow in two steps. The first regulation is the R-D control done for each program individually without sharing the encoding information of the other programs. The second regulation is the centralized buffer regulation during which the video quality of each stored flow in the buffers, considered located in a MANE, are analyzed and compared in order to adjust the amount of data to extract from the buffers such that comparable video quality is maintained in the N buffers. The two regulations processes are performed in a closed loop since the rate regulation feds back to each video server the level of its associated buffer to help the R-D controllers adapt their R-D trade-off so that the buffer within the MANE reaches some reference level. These two regulations processes are detailed in the following.

Rate allocation: Quality fairness

The available transmission rate R c is allocated among video programs in a centralized way within the MANE to meet the fairness constraint. For that purpose, we assume that each video server i provides to the MANE the distortion of the j-th frame D ij and the corresponding Peak Signal-to-Noise Ratio (PSNR) P ij is

P ij = 10log 10 ( 255 2 D ij ). (1)
The average PSNR Pj among the N program is then computed. We adopt a control-theoretic approach to regulate the video quality and model the system as a feedback control system. Using a a proportional controller (P), the rate α ijj R c allocated to the j-th frame in the i-th program between time j and time j + 1 is evaluated as

α ij = 1 N + K P α .( Pj -P ij ), 0 α ij 1, i = 1 . . . N (2)
where K P α is a proportional correction gain. More (less) rate is allocated to programs with a PSNR smaller (larger) than Pj . With this allocation, the evolution of the level of the buffer of the i-th program between time j and time j + 1 is

B ij+1 = B ij + (R ij -α ij R c ) ∆t. (3)
As a consequence, the level of the buffer of programs which produce encoded video with smaller PSNR than the average PSNR Pj decreases faster than programs producing video with a larger PSNR. For those programs, the buffer level may even decrease.

R-D control

For each video coder an individual R-D control process is assumed to be performed at a frame level. This control may be done, for example, by adjusting the encoding parameters with H264-AVC encoders or the number of SNR layer to transmit with H264-SVC encoders. A reference coding rate R 0 = R c /N is provided to each server. The rate R ij at which the j-th frame is adjusted depending on the discrepancy between the buffer level of the i-th program and some reference buffer level B 0 as follows

R ij = R 0 -K P R B ij-1 -B 0 ∆T , (4)
where K P R is a proportional correction factor. PI controllers are widely used in feedback systems, which the integral term is proportional to both the magnitude of the error and the duration of the error. PI controller is that its integral part may drive the steadystate error to zero, which enhances the steady-state performance of the system. This is a perfect match for the objective of regulating the encoding rate to reach stabilized buffer fullness for the N buffers in the MANE. Using a PI controller, the allocated rates fir the j-th frame and considered between time j -1 and time j are

R ij = R 0 -K P R B ij-1 -B 0 ∆T -K I R j k=1 B ik-1 -B 0 ∆T , (5)
where K I R is a Integral correction factor. One sees from (4) and ( 5) that the encoding rate increases when the buffer is below the reference level. This control technique only requires the MANE to feed back the buffer level B ij-1 to each remote video encoder.

Fig. 1 :

 1 Fig. 1: Statistical multiplexing system of video streams provided by remote servers

  The advantage of

			Proxy: Media Aware Network Element
	Video server	Buffer state
	Video Encoder 1	R-D control process 1			R	t 1
			Buffer Program 1
	Video server	Buffer state
	Video Encoder 2		Buffer Program 2	Rate allocation: scheduling	R	t 2	R c	Broadcast Channel
				decision
	Video server	Buffer state
	Video				R	t N
	Encoder N
			Buffer Program N

Analysis and results

In this section, we consider Gaussian sources clipped to [0, 255], with mean 128, and variance σ 2 ij with i = {1, . . . , N } that may vary with time. From time to time, there is a jump in the variance, small variations are modeled by uniformly distributed additional variances z ij distributed in the interval [-1, 1], so that

) Temporal variations of the source variance correspond to scene changes in the video sequences. The R-D function of each source is

with i = 1, . . . , N . Combining [START_REF] Changuel | End-to-end stochastic scheduling of scalable video overtime-varying channels[END_REF] and [START_REF] Chengsheng | An efficient two-pass VBR encoding algorithm for H.264[END_REF], the PSNR can be expressed as a the following linear function of the rate

with

) and γ = 6.

Equilibrium

To study the equilibrium of the proposed system, we assume that the level of the buffers B ij is constant with time. Thus, using a P controller and combining (2), (3), and ( 9), one gets the following system of equations

corresponding to the rates at equilibrium.

Experiments

The performance of the proposed statistical multiplexing system is evaluated with three Gaussian sources each one is characterized by a time-varying variance. The transmission rate allocation is performed according to (2) using a P controller with K P α = 0.07. The encoding rates are obtained from (4) when using P controller with K P R = 0.05 and from (5) when using PI controller with K P R = 0.05 and K I R = 0.006. The instantaneous variance σ 2 of the source, the buffer fullness, the PSNR, the rate allocation α, and the encoding rate of each program are represented in Figure 2. These curves have been obtained with B 0 = 30 bits, R c = 10 bits/s and ∆t = 1 s.

From Figure 2, one can see that the control system reaches equilibrium after short transient modes to satisfy the quality fairness constraint. At equilibrium, thanks to the rate allocation control in (2) based on the quality fairness constraint, almost equal PSNR are obtained for all programs. Using a P controller in the feedback control system of the R-D control, a significant discrepancy between the level of the buffers at equilibrium and the refer- Study of the equilibrium of the proposed system using a P controller ence buffer level is not reached. This discrepancy are fully compensated as shown in Figure 3 by using a PI control as in ( 5) instead of a P control in (4).

Conclusion

In this paper, we propose a control system that allow the statistical multiplexing of multiple video sources in a distributed way. A MANE is in charge of buffering the encoded streams and building a centralized bandwidth allocation among programs. This first regulation system should determine the appropriate amount of data to extract from each buffer based on quality information provided by each video coder in order to try to satisfy the fairness constraint. The buffer fullness is then fed back to the R-D control system of each video server to adapt the encoding rates and the video quality so that a reference buffer level is achieved. This regulation process involves a P controller for the rate allocation process and a P or PI controller for the R-D control process. Experimental tests using Gaussian signals show that the targeted quality fairness constraint is achieved. In addition, thanks to the PI controller, results show that the proposed control system allows an equilibrium state of the buffer fullness to be reached. This control system has now to be adapted to perform the control of actual video coders. The main difficulty comes from the fact that not all R-D points may be reached in a continuous way (this would require very fine grain scalability). The R-D point closest to the requested one will then have to be applied.