C Imbert

A Mellet

Electrified thin films: Global existence of non-negative solutions

Keywords: Higher order equation, Non-local equation, Thin film equation, Non-negative solutions MSC: 35G25, 35K25, 35A01, 35B09

We consider an equation modeling the evolution of a viscous liquid thin film wetting a horizontal solid substrate destabilized by an electric field normal to the substrate. The effects of the electric field are modeled by a lower order non-local term. We introduce the good functional analysis framework to study this equation on a bounded domain and prove the existence of weak solutions defined globally in time for general initial data (with finite energy).

Introduction

In this paper, we construct solutions for a thin film type equation with a destabilizing singular integral term. This term models the effects of an electric field (see [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF]). From the analytical point of view, this paper belongs to the large body of literature devoted to the thin film equation with destabilizing terms such as long-wave unstable thin film problems [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF][START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF][START_REF] Witelski | Blowup and dissipation in a critical-case unstable thin film equation[END_REF] or the Kuramato-Sivashinsky equation in combustion and solidification [START_REF] Frankel | On the equation of a curved flame front[END_REF][START_REF] Frankel | On a free boundary problem associated with combustion and solidification[END_REF].

More precisely, we are considering the following equation, which is introduced by Tseluiko and Papageorgiou in [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF] (see also [START_REF] Papageorgiou | Gravity capillary waves in fluid layers under normal electric fields[END_REF]):

u t + u 3 (cu xx -αu -λI(u)) x x = 0 x ∈ [0, L], t > 0 (1)
(in [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF], (1) is supplemented with periodic boundary conditions). This equation models the evolution of a liquid thin film (of height u) wetting a horizontal solid

The existence of non-negative weak solutions for (2) was first established by F. Bernis and A. Friedman [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF] for f (u) = u n , n > 1. Further results (existence for n > 0 and further regularity results) were later obtained, by similar technics, in particular by E. Beretta, M. Bertsch and R. Dal Passo [START_REF] Beretta | Nonnegative solutions of a fourth-order nonlinear degenerate parabolic equation[END_REF] and A. Bertozzi and M. Pugh [START_REF] Bertozzi | The lubrication approximation for thin viscous films: the moving contact line with a "porous media" cut-off of van der Waals interactions[END_REF][START_REF]The lubrication approximation for thin viscous films: regularity and long-time behavior of weak solutions[END_REF]. Results in higher dimension were obtained in particular by Grün in [START_REF]On Bernis' interpolation inequalities in multiple space dimensions[END_REF][START_REF] Grün | Degenerate parabolic differential equations of fourth order and a plasticity model with non-local hardening[END_REF][START_REF] Passo | On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions[END_REF]. The thin film equation with lower order destabilizing terms has also received a lot of interest. In particular in [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF][START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF], the following equation

u t + (f (u)u xxx -g(u)u x ) x = 0 (3)
is considered. Such a destabilizing term (which, unlike that of [START_REF] Agranovich | On Fourier series in eigenfunctions of elliptic boundary value problems[END_REF], is a local term of order 2) models, for instance, the effects of gravity for a hanging thin film, or van der Waals type interactions with the solid substrate. In [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF], the nonlinearities f (u) = u n and g(u) = u m are considered and it is proved (among other things) that there is no blow-up for m < n + 2. In [START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF], for f (u) = u and g(u) = u m , it is proved that there is blow-up for m ≥ 3 and initial data in H 1 (R) with negative "energy". The reader is referred to [START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF] for a precise statement.

In our equation [START_REF] Agranovich | On Fourier series in eigenfunctions of elliptic boundary value problems[END_REF], the nonlinearities in front of the stabilizing and destabilizing terms are the same (f (u) = g(u) = u 3 ), but the destabilizing term is elliptic of order 3 and is nonlocal in space. It is known (see [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF]) that positive smooth solutions of (1) do not blow up. The goal of the present paper is to prove the existence of global in time weak solutions for [START_REF] Agranovich | On Fourier series in eigenfunctions of elliptic boundary value problems[END_REF].

Note that besides the existence of solutions, many important properties of the thin film equation [START_REF] Beretta | Nonnegative solutions of a fourth-order nonlinear degenerate parabolic equation[END_REF] have been investigated (finite speed of propagation of the support, waiting time phenomenon, existence of source-type solutions etc.).

The key tools in many of these studies are various delicate integral inequalities (in particular the so-called α-entropy inequalities and local entropy and energy inequalities. See [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF][START_REF]Nonnegative solutions of a fourth-order nonlinear degenerate parabolic equation[END_REF][START_REF]The lubrication approximation for thin viscous films: regularity and long-time behavior of weak solutions[END_REF][START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF]). It is not clear that similar functional inequalities holds for [START_REF] Agranovich | On Fourier series in eigenfunctions of elliptic boundary value problems[END_REF]. One reason is that the algebra involving the operator I(u) is considerably more difficult than that of the Laplace operator. Another reason, is the obvious difficulty in deriving local estimates (due to the nonlocal nature of the operator I(u)). For that reason, we only address the existence issue in this paper.

As in [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF][START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF], the main difficulty in proving the existence of solutions for [START_REF] Agranovich | On Fourier series in eigenfunctions of elliptic boundary value problems[END_REF] comes from the fact that the energy (see [START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF] below) can take negative values. In order to obtain H 1 a priori estimate, one thus has to use the conservation of mass which, for non-negative solutions, gives a global in time L 1 bound for the solution (see Lemma 1).

With such an estimate in hand, the existence of global in time solutions should follow from the construction of approximated solutions satisfying the right functional inequalities. Typically, one needs to regularize the mobility coefficient u 3 . One way to proceed is to replace the coefficient u 3 with u 3 + ε so that the equation becomes strictly parabolic. However, for such a regularized equation one cannot show the existence of non-negative solutions (the maximum principle does not hold for fourth order parabolic equations) and Lemma 1 is of no use. An alternative regularization is to replace the mobility coefficient u 3 with a function f ε (u) which satisfies in particular f ε (u) ∼ u 4 /ε. For such (more degenerate) mobility coefficient, solutions are expected to be strictly positive and therefore smooth. This second regularization procedure was first suggested by Bernis and Friedman [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF] and is used in particular by Bertozzi and Pugh [START_REF] Bertozzi | The lubrication approximation for thin viscous films: the moving contact line with a "porous media" cut-off of van der Waals interactions[END_REF][START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF]. However, the local in time existence for such a degenerate equation is not clear to us, since the corresponding proofs in [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF] rely on Schauder estimates which are not classical (and perhaps tedious) with our non-local singular term I(u). For this reason, we choose the first regularization approach; this requires us to pay attention to the lack of positivity of the approximated solutions. In particular, the L 1 norm is not controlled, and the H 1 norm will be controlled by combining the energy inequality with the entropy inequality. The idea of combining the energy together with the entropy in order to get a Lyapunov functional appeared previously in [START_REF] Shishkov | On the equation of the flow of thin films with nonlinear convection in multidimensional domains[END_REF] where a thin film equation with a nonlinear drift term is thouroughly studied in all dimensions.

The main contribution of this paper is thus to introduce the precise functional analysis framework to be used to treat the term I(u) and to provide a method for constructing weak solutions satisfying the proper a priori estimates.

Rather than working in the periodic setting, we will consider equation (1) on a bounded domain with Neumann boundary conditions (these Neumann conditions can be interpreted as the usual contact angle conditions and seem physically more relevant -the periodic framework could be treated as well with minor modifications). Further details about the derivation of (1) will be given in Section 2. Since the gravity term is of lower order than the electric field term, it is of limited interest in the mathematical theory developed in this paper. We will thus take α = 0 and c = λ = 1.

We thus consider the following problem:

     u t + f (u)(u xx -I(u)) x x = 0 for x ∈ Ω, t > 0 u x = 0, f (u)(u xx -I(u)) x = 0 for x ∈ ∂Ω, t > 0 u(x, 0) = u 0 (x) for x ∈ Ω. (4)
The domain Ω is a bounded interval in R; in the sequel, we will always take Ω = (0, 1). The mobility coefficient f

(u) is a C 1 function f : [0, +∞) → (0, +∞) satisfying f (u) ∼ u n as u → 0 (5)
for some n > 1. The operator I is a non-local elliptic operator of order 1 which will be defined precisely in Section 3 as the square root of the Laplace operator with Neumann boundary conditions (we have to be very careful with the definition of I in a bounded domain).

A priori estimates. As for the thin film equation ( 2), we prove the existence of solutions for (4) using a regularization/stability argument. The main tools are integral inequalities which provide the necessary compactness. Besides the conservation of mass, we will see that the solution u of (4) satisfies two important integral inequalities: We define the energy E(u) and the entropy e(u) by

E(u)(t) = 1 2 Ω (u 2 x (t) + u(t)Iu(t))dx and e(u)(t) = 1 2 Ω G(u(t))dx
where G is a non-negative convex function such that f G ′′ = 1. Classical solutions of (4) then satisfy:

E(u)(t) + t 0 Ω f (u) (u xx -I(u)) x 2 dx ds ≤ E(u 0 ), (6) e(u)(t)
+ t 0 (u xx ) 2 dx ds + t 0 u x I(u) x dx ds ≤ e(u 0 ). (7)
Similar inequalities hold for the thin film equation [START_REF] Beretta | Nonnegative solutions of a fourth-order nonlinear degenerate parabolic equation[END_REF]. However, we see here the destabilizing effect of the nonlocal term I(u): First, we note that as in [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF][START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF] the energy E(u) can be written as the difference of two non-negative quantities:

E(u)(t) = u(t) 2 Ḣ1 (Ω) -u(t) 2 Ḣ 1 2 (Ω) , (8)
and may thus take negative values. Similarly, the entropy dissipation can be written as

t 0 (u xx ) 2 dx ds + t 0 u x I(u) x dx ds = u(t) 2 Ḣ2 (Ω) -u(t) 2 Ḣ 3 2 N (Ω)
, so the entropy may not be decreasing. Nevertheless, it is reasonable to expect (4) to have solutions that exist for all times. Indeed, as shown in [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF], the conservation of mass, the inequality [START_REF]The lubrication approximation for thin viscous films: regularity and long-time behavior of weak solutions[END_REF] and the following functional inequality (see Lemma 1)

u 2 Ḣ1 (Ω) ≤ αE(u) + β u 2 L 1 (Ω) , ∀u ∈ H 1 (Ω),
implies that non-negative solution remains bounded in L ∞ (0, T ; H 1 (Ω)) for all time T . Furthermore, the interpolation inequality

u(t) 2 Ḣ 3 2 N ≤ C u(t) Ḣ1 u(t) Ḣ2 ≤ 1 2 u(t) 2 Ḣ2 + C 2 u(t)
G ′′ (u) = 1 f (u)
for all u > 0.

Theorem 1. Let n > 1 and u 0 ∈ H 1 (Ω) be such that u 0 ≥ 0 and

Ω G(u 0 ) dx < ∞. (9)
For all T > 0 there exists a function u(t, x) ≥ 0 with

u ∈ C(0, T ; L 2 (Ω)) ∩ L ∞ (0, T ; H 1 (Ω)), u x ∈ L 2 (0, T ; H 1 0 (Ω))
such that, for all φ ∈ D([0, T ) × Ω) satisfying φ x = 0 on (0, T ) × ∂Ω,

Q uφ t -f (u)[u xx -I(u)]φ xx -f ′ (u)u x (u xx -I(u))φ x dt dx + Ω u 0 (x)φ(0, x) dx = 0. ( 10
)
Moreover, the function u satisfies for every t ∈ [0, T ],

Ω u(t, x) dx = Ω u 0 (x) dx, E(u(t)) + t 0 Ω f (u) (u xx -I(u)) x 2 ds dx ≤ E(u 0 ), ( 11
) Ω G(u(t)) dx + t 0 Ω (u xx ) 2 + u x I(u) x dsdx ≤ Ω G(u 0 ) dx. ( 12
)
We point out that the weak formulation [START_REF] Passo | On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions[END_REF] involve two integrations by parts. Our second main result is concerned with non-negative initial data whose entropies are possibly infinite (this is the case if u 0 vanishes on an open subset of Ω and n ≥ 2). In that case, only one integration by parts is possible, and the solutions that we construct are weaker than those constructed in Theorem 1. In particular, the equation is only satisfied on the positivity set of the solution and the boundary conditions are satisfied in a weaker sense.

Theorem 2. Assume n > 1 and let u 0 ∈ H 1 (Ω) be such that u 0 ≥ 0. For all T > 0 there exists a function u(t, x) ≥ 0 such that [START_REF] Grün | Degenerate parabolic differential equations of fourth order and a plasticity model with non-local hardening[END_REF] where P = {(x, t) ∈ Q : u(x, t) > 0, t > 0}. Moreover, the function u satisfies the conservation of mass and the energy inequality [START_REF] Frankel | On a free boundary problem associated with combustion and solidification[END_REF].

u ∈ C(0, T ; L 2 (Ω)) ∩ L ∞ (0, T ; H 1 (Ω)) ∩ C 1 2 , 1 8 (Ω × (0, T )) such that f (u)[u xx -I(u)] x ∈ L 2 (P ) and such that, for all φ ∈ D([0, T ) × Ω), Q uφ t dt dx + P f (u)[u xx -I(u)] x φ x dt dx + Ω u 0 (x)φ(0, x) dx = 0
Finally, u x vanishes at all points (x, t) of ∂Ω × (0, T ) such that u(x, t) = 0.

Comments. These results are comparable to those of [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF] when λ = 0. The reader might be surprised that they are presented in a different order than in [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF]. The reason has to do with the proofs; indeed, in contrast with [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF], weak solutions given by Theorem 2 are constructed as limits of the solutions given by Theorem 1. This is because the entropy is needed in order to construct the nonnegative solutions. See the discussion at the beginning of Section 5 for further details.

As pointed out earlier, the non-linearities in front of the stabilizing (u xxx ) and destabilizing ((I(u)) x ) are the same. This is in contrast with the work of Bertozzi and Pugh [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF][START_REF]Finite-time blow-up of solutions of some long-wave unstable thin film equations[END_REF]. By analogy with (3), one could consider the equation

u t + (u n u xxx -u m (I(u)) x ) x = 0
in which case a scaling analysis similar to that of [START_REF] Bertozzi | Long-wave instabilities and saturation in thin film equations[END_REF] suggests that blow up can only occur if m ≥ n + 1. However, to our knowledge, there is no physical motivation for such a generalization (in our case).

Organization of the article. In Section 2, we give more details about the physical model leading to (4). We gather, in Section 3, material that will be used throughout the paper. In particular, we detail the functional analysis framework and the definition of the non-local operator I (which is similar to that used in [START_REF] Imbert | Existence of solutions for a higher order non-local equation appearing in crack dynamics[END_REF]). Section 4, 5 and 6 are devoted to the proofs of the main results. Finally, we give in Appendix a technical result which is more or less classical.

Acknowledgements. The first author was partially supported by the French Ministry of Research (ANR project EVOL). The second author was partially supported by NSF Grant DMS-0901340.

Physical model

In this section, we briefly recall the derivation of (4) (see [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF] for further details). We consider a viscous liquid film which completely wets a solid horizontal substrate and is constrained between two solid walls (at x = 0 and x = 1), see Figure 1. The fluid is Newtonian and is assumed to be a perfect conductor. The substrate is a grounded electrode held at zero voltage. Thanks to the presence of another electrode (at infinity), an electric field E is created which is constant at infinity (in the direction perpendicular to the substrate):

E(x, y) -→ (0, E 0 )
as y → +∞.

The height of the fluid is denoted by u(t, x). Under the assumptions of the lubrication approximation, it is classical that the evolution of u is described by Poiseuille's law:

u t -∂ x u 3 3µ ∂ x p = 0 ( 14
)
where p is the pressure at the free surface of the fluid y = u(t, x). This pressure is the sum of three terms:

1. The capillary pressure due to surface tension, which can be approximated by p 1 ∼ -σu xx (replacing the mean curvature operator by the Laplacian).

2. The effect of gravity, given by p 2 = gu.

3. The additional pressure due to the action of the electric field E.

To compute the third term appearing in the pressure, we introduce the potential V such that E = -∇V , which satisfies ∆V = 0 for y ≥ u(x) and V (x, y) = 0 on y = u(x).

The condition at y → ∞ means that we can write

V ∼ E 0 (Y 0 -y)
with (using standard linear approximation)

   ∆Y 0 = 0 for y > 0, x ∈ Ω ∇Y 0 → 0 as y → ∞, x ∈ Ω Y 0 (x, 0) = u(x) for x ∈ Ω. ( 15
)
At the boundary of the cylinder, we assume that the electric field has no horizontal component:

∂ x V = 0, for x ∈ ∂Ω, y > 0.
The pressure exerted by the electric field is then proportional to

p 3 = γE y = -γ∂ y V (x, 0) = -γE 0 (∂ y Y 0 -1).
The application u → ∂ y Y 0 (x, 0) is a Dirichlet-to-Neumann map for the harmonic extension problem [START_REF] Imbert | Existence of solutions for a higher order non-local equation appearing in crack dynamics[END_REF]. We denote this operator by I(u). We will see in Section 3 that I(u) is in fact the square root of the Laplace operator on the interval Ω with homogeneous Neumann boundary conditions. We thus have

p = p 1 + p 2 + p 3 = -σu xx + gu -γE 0 I(u) + c 0
for some constant c 0 , and we obtain (1) with c = σ 3µ , α = g 3µ and λ = -γE0 3µ . Note that Poiseuille's law [START_REF]On Bernis' interpolation inequalities in multiple space dimensions[END_REF] is obtained under the no-slip condition for the fluid along the solid support. Other conditions, such as the Navier slip condition leads to f (u) = u 3 + Λu s with s = 1 or s = 2. This explains the interest of the community for general diffusion coefficient f (u).

Boundary conditions. Along the boundary ∂Ω, the fluid is in contact with a solid wall. It is thus natural to consider a contact angle condition at x = 0 and x = 1: Assuming that the contact angle is equal to π/2, we then get the boundary condition u x = 0 on ∂Ω.

In [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF], the authors derive their analytic results in a periodic setting which is obtained by considering the even extension of u to the interval (-1, 1) (recall that Ω = (0, 1)) and then taking the periodic extension (with period 2) to R.

Finally, since the equation is of order 4, we need an additional boundary condition. We thus assume that u satisfies the following null-flux condition

u 3 (u xx -I(u)) x = 0
on ∂Ω which will guarantee the conservation of mass.

Preliminaries

In this section, we recall how the operator I is defined (see [START_REF] Imbert | Existence of solutions for a higher order non-local equation appearing in crack dynamics[END_REF]) and give the functional analysis results that we will need to prove the main theorem. A very similar operator, with Dirichlet boundary conditions rather than Neumann boundary conditions, was studied by Cabré and Tan [START_REF] Cabré | Positive solutions of nonlinear problems involving the square root of the laplacian[END_REF].

Functional spaces

The space H s N (Ω). We denote by {λ k , ϕ k } k=0,1,2... the eigenvalues and corresponding eigenfunctions of the Laplace operator in Ω with Neumann boundary conditions on ∂Ω:

-

∆ϕ k = λ k ϕ k in Ω ∂ ν ϕ k = 0 on ∂Ω, (16)
normalized so that Ω ϕ 2 k dx = 1. When Ω = (0, 1), we have

λ 0 = 0 , ϕ 0 (x) = 1 and λ k = (kπ) 2 , ϕ k (x) = √ 2 cos(kπx) k = 1, 2, 3, . . .
The ϕ k 's clearly form an orthonormal basis of L 2 (Ω). Furthermore, the ϕ k 's also form an orthogonal basis of the space H s N (Ω) defined by

H s N (Ω) = u = ∞ k=0 c k ϕ k ; ∞ k=0 c 2 k (1 + λ s k ) < +∞ equipped with the norm ||u|| 2 H s N (Ω) = ∞ k=0 c 2 k (1 + λ s k )
or equivalently (noting that c 0 = Ω u(x) dx and λ k ≥ 1 for k ≥ 1):

||u|| 2 H s N (Ω) = u 2 L 1 (Ω) + u 2 Ḣs N (Ω)
where the homogeneous norm is given by:

u 2 Ḣs N (Ω) = ∞ k=1 c 2 k λ s k .
A characterisation of H s N (Ω). The precise description of the space H s N (Ω) is a classical problem.

Intuitively, for s < 3/2, the boundary condition u ν = 0 does not make sense, and one can show that (see Agranovich and Amosov [START_REF] Agranovich | On Fourier series in eigenfunctions of elliptic boundary value problems[END_REF] and references therein):

H s N (Ω) = H s (Ω) for all 0 ≤ s < 3 2 .
In particular, we have H

1 2 N (Ω) = H 1 2
(Ω) and we will see later that

u 2 Ḣ 1 2 (Ω) = Ω Ω (u(y) -u(x)) 2 ν(x, y)dxdy
where ν(x, y) is a given positive function; see (20) below.

For s > 3/2, the Neumann condition has to be taken into account, and we have in particular

H 2 N (Ω) = {u ∈ H 2 (Ω) ; u ν = 0 on ∂Ω}
which will play a particular role in the sequel. More generally, a similar characterization holds for 3/2 < s < 7/2. For s > 7/2, additional boundary conditions would have to be taken into account, but we will not use such spaces in this paper. In Section 4, we will also work with the space H 3 N (Ω) which is exactly the set of functions in H 3 (Ω) satisfying u ν = 0 on ∂Ω.

The case s = 3/2 is critical (note that u ν | ∂Ω is not well defined in that space) and one can show that

H 3 2 N (Ω) = u ∈ H 3 2 (Ω) ; Ω u 2 x d(x) dx < ∞
where d(x) denotes the distance to ∂Ω. A similar result appears in [START_REF] Cabré | Positive solutions of nonlinear problems involving the square root of the laplacian[END_REF]; more precisely, such a characterization of H 3 2

N (Ω) can be obtained by considering functions u such that u x ∈ V 0 (Ω) where V 0 (Ω) is defined in [START_REF] Cabré | Positive solutions of nonlinear problems involving the square root of the laplacian[END_REF] as the equivalent of our space H 1/2 N (Ω) with Dirichlet rather than Neumann boundary conditions. We do not dwell on this issue since we will not need this result in this paper.

The operator I

As it is explained in the Introduction, the operator I is related to the computation of the pressure as a function of the height of the fluid. Spectral definition. With λ k and ϕ k defined by ( 16), we define the operator

I : ∞ k=0 c k ϕ k -→ - ∞ k=0 c k λ 1 2 k ϕ k ( 17
)
which clearly maps H 1 (Ω) onto L 2 (Ω) and

H 2 N (Ω) onto H 1 (Ω).
Dirichlet-to-Neuman map. We now check that this definition of the operator I is the same as the one given in Section 2, namely I is the Dirichletto-Neumann operator associated with the Laplace operator supplemented with Neumann boundary conditions:

We consider the following extension problem:

   -∆v = 0 in Ω × (0, +∞), v(x, 0) = u(x) on Ω, v ν = 0 on ∂Ω × (0, ∞). (18)
Then, we can show (see [START_REF] Imbert | Existence of solutions for a higher order non-local equation appearing in crack dynamics[END_REF]):

Proposition 1 ([15]). For all u ∈ H 1 2
N (Ω), there exists a unique extension v ∈ H 1 (Ω × (0, +∞)) solution of [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF].

Furthermore, if u(x) = ∞ k=1 c k ϕ k (x), then v(x, y) = ∞ k=1 c k ϕ k (x) exp(-λ 1 2 k y). ( 19
)
and we have:

Proposition 2 ([15]
). For all u ∈ H 2 N (Ω), we have

I(u)(x) = - ∂v ∂ν (x, 0) = ∂ y v(x, 0) for all x ∈ Ω,
where v is the unique harmonic extension solution of [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF].

Furthermore I • I(u) = -∆u.

Integral representation. Finally, the operator I can also be represented as a singular integral operator:

Proposition 3 ([15]
). Consider a smooth function u : Ω → R. Then for all x ∈ Ω,

I(u)(x) = Ω (u(y) -u(x))ν(x, y)dy
where ν(x, y) is defined as follows: for all x, y ∈ Ω,

ν(x, y) = π 2 1 1 -cos(π(x -y)) + 1 1 -cos(π(x + y)) . ( 20
)

Functional equalities and inequalities

Equalities. The semi-norms

|| • || Ḣ 1 2 (Ω) , || • || Ḣ1 (Ω) , || • || Ḣ 3 2 N (Ω)
and

|| • || Ḣ2 N (Ω)
are related to the operator I by equalities which will be used repeatedly.

Proposition 4 (The operator I and several semi-norms - [START_REF] Imbert | Existence of solutions for a higher order non-local equation appearing in crack dynamics[END_REF]).

For all u ∈ H 1 2 (Ω), we have

-uI(u) dx = 1 2 Ω Ω (u(x) -u(y)) 2 ν(x, y)dxdy = ||u|| 2 Ḣ 1 2 (Ω)
where ν is defined in [START_REF] Witelski | Blowup and dissipation in a critical-case unstable thin film equation[END_REF].

For all u ∈ H 2 N (Ω), we have

- Ω u x I(u) x dx = ||u|| 2 Ḣ 3 2 N (Ω)
.

For all k ∈ N and u ∈ H k+1 N (Ω), we have

Ω (∂ k x I(u)) 2 dx = u 2 Ḣk+1 N (Ω) . (21)
Inequalities. First, we recall the following Nash inequality:

u L 2 (Ω) ≤ C u 1 3 H 1 (Ω) u 2 3
L 1 (Ω) . It implies in particular that,

u 2 L 2 (Ω) ≤ 1 2 u 2 Ḣ1 (Ω) + C u 2 L 1 (Ω) . (22)
This inequality will allow us to control the H 1 norm by the energy E(u) and the L 1 norm. Indeed, we recall that the energy is defined by

E(u) = 1 2 Ω |u x | 2 + u I(u) dx = 1 2 u 2 Ḣ1 (Ω) - 1 2 u 2 Ḣ 1 2 (Ω) . (23)
We then have:

Lemma 1. There exist positive constants α, β such that for all u ∈ H 1 (Ω),

u 2 Ḣ1 (Ω) ≤ αE(u) + β u 2 L 1 (Ω) .
Remark 1. See also Lemma 4.1 in [START_REF] Tseluiko | Nonlinear dynamics of electrified thin liquid films[END_REF].

Proof. We have

u 2 Ḣ1 (Ω) = 2E(u) + u 2 Ḣ 1 2 (Ω)
, and using (21) with k = 0 and (22), we get:

u 2 Ḣ 1 2 (Ω) = -u I(u) dx ≤ u L 2 (Ω) I(u) L 2 (Ω) ≤ u L 2 (Ω) u Ḣ1 (Ω) ≤ 1 2 u 2 L 2 (Ω) + 1 2 u 2 Ḣ1 (Ω) ≤ 3 4 u 2 Ḣ1 (Ω) + C 2 u 2 L 1 (Ω) ,
hence the result.

A regularized problem

We can now turn to the proof of Theorem 1 and 2. As usual, we introduce the following regularized equation:

     u t + (f ε (u)(u xx -I(u)) x ) x = 0 for x ∈ Ω, t > 0 u x = 0, f ε (u)(u xx -I(u)) x = 0 for x ∈ ∂Ω, t > 0 u(x, 0) = u 0 (x) for x ∈ Ω (24)
where the mobility coefficient f (u) is approximated by

f ε (u) = min(max(ε, f (|u|)), M ) which satisfies ε ≤ f ε (u) ≤ M for all u ∈ R.
Ultimately, we will show that the solution u satisfies 0 ≤ u(t, x) ≤ M 0 for some constant M 0 independent of M , so that we do not have to worry about M (provided we take it large enough). The ε is of course the most important parameter in the regularization since it makes (24) non-degenerate. The proof of Theorem 1 consists of two parts: First, we have to show that the regularized equation (24) has a solution (which may take negative values). Then we must pass to the limit ε → 0 and show that we obtain a non-negative solution of (4).

In this section, we prove the first part. Namely, we prove (to be compared with Theorem 1.1 in [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF]): Theorem 3. Let u 0 ∈ H 1 (Ω). For all T > 0 there exists a function u ε (t, x) such that

u ε ∈ C(0, T ; L 2 (Ω)) ∩ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 3 N (Ω)) such that, for all φ ∈ D([0, T ) × Ω), Q u ε φ t + f ε (u ε )[u ε xx -I(u ε )] x φ x dt dx + Ω u 0 (x)φ(0, x) dx = 0. ( 25
)
Moreover, the function u ε satisfies for every t ∈ [0, T ],

Ω u ε (t, x) dx = Ω u 0 (x) dx, E(u ε (t)) + t 0 Ω f ε (u ε ) (u ε xx -I(u ε )) x 2 ds dx ≤ E(u 0 ), (26)
and

Ω (u ε x ) 2 dx + t 0 Ω f ε (u ε )(u ε xxx ) 2 ds dx ≤ Ω ((u 0 ) x ) 2 dx + t 0 Ω f ε (u ε )u ε xxx (I(u ε )) x ds dx (27)
and

Ω G ε (u ε (t)) dx + t 0 Ω (u ε xx ) 2 + u ε x I(u ε ) x dsdx ≤ Ω G ε (u 0 ) dx. ( 28
)
where G ε is a non-negative function such that f ε G ′′ ε = 1. Finally, u ε is 1 2 -Hölder continuous with respect to x and 1 8 -Hölder continuous with respect to t; more precisely, there exists a constant C 0 only depending on

Ω and u ε L ∞ (0,T ;H 1 (Ω)) and f ε (u ε )[u ε xx -I(u ε )] x L 2 (Q) such that u ε C 1 2 , 1 8
t,x (Q) ≤ C 0 . ( 29
)
Remark 2. This theorem is very similar to Theorem 1.1 in [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF], and some steps in our proof follow along the lines of [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF]. For instance, getting the Hölder estimates from the L ∞ H 1 estimate is done in the same way. However, the main difficulty in proving this theorem is precisely to get the L ∞ H 1 estimate; this step is not straightforward at all and this is a significant difference with [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF].

Proof of Theorem 3. Theorem 3 follows from a fixed point argument: For some T * , we denote V = L 2 (0, T * ; H 2 N (Ω)) and we define the application

F : V → V such that for v ∈ V , F (v) is the solution u of      u t + (f ε (v)(u xxx -I(v) x )) x = 0 for x ∈ Ω, t > 0 u x = 0, f ε (v)(u xx -I(v)) x = 0 for x ∈ ∂Ω, t > 0 u(x, 0) = u 0 (x) for x ∈ Ω. ( 30
)
The fact that F is well defined follows from the observation that for v ∈ V , we have

a(t, x) = f ε (v(t, x)) ∈ [ε, M ] and g(t, x) = I(v) x ∈ L 2 (Q)
and the following proposition:

Proposition 5. Consider u 0 ∈ H 1 (Ω) and a(t, x) ∈ L ∞ (Q) such that ε ≤ a(t, x) ≤ M a.e. in Q. If g ∈ L 2 (Q), then there exists a function u ∈ C(0, T ; L 2 (Ω)) ∩ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 3 N (Ω))
and

Q uφ t + a(u xxx -g)φ x dt dx + Ω u 0 (x)φ(0, x) dx = 0 for all φ ∈ D([0, T ) × Ω). Moreover, for every t ∈ [0, T ], u satisfies Ω u(t, x) dx = Ω u 0 (x) dx
and

Ω (u x ) 2 (t) dx + 1 2 t 0 Ω a(u xxx ) 2 ds dx ≤ M 2 t 0 Ω g 2 ds dx + Ω (u 0 ) 2 x dx. (31)
Furthermore, u is 1 2 -Hölder continuous with respect to x and 1 8 -Hölder continuous with respect to t; more precisely, there exists a constant C 0 only depending on Ω and u L ∞ (0,T ;H 1 (Ω)) and

u xxx -g L 2 (Q) such that u C 1 2 , 1 8
t,x (Q) ≤ C 0 . ( 32
)
This proposition is a very natural existence result for the fourth order linear parabolic equation

u t + (au xxx ) x = (ag) x .
Its proof is fairly classical, we give some details in Appendix A for the interested reader.

Next, we show the following result:

Lemma 2. There exists a (small) time T * > 0, depending only on ε, M and Ω, such that F has a fixed point u in V = L 2 (0, T * ; H 2 N (Ω)) for any initial data u 0 ∈ H 1 (Ω). Furthermore, u satisfies

u V ≤ R||u 0 || Ḣ1 (Ω) and u L ∞ (0,T * ; Ḣ1 (Ω)) ≤ √ 2||u 0 || Ḣ1 (Ω) . (33)
Before proving this lemma, let us complete the proof of Theorem 3.

Construction of a solution for large times. Lemma 2 gives the existence of a solution u ε 1 of (24) defined for t ∈ [0, T * ]. Since T * does not depend on the initial condition, we can apply Lemma 2 to construct a solution u ε 2 in [T * , 2T * ] with initial condition u ε 1 (T * , x) which is H 1 (Ω) by ( 33). This way, we obtain a solution u ε of (24) on the time interval [0, 2T * ]. Note that we also have

u ε L ∞ (0,2T * ; Ḣ1 (Ω)) ≤ √ 2 2 ||u 0 || Ḣ1 (Ω) .
Iterating this argument, we construct a solution u ε on any interval [0, T ] satisfying in particular, for all k ∈ N such that kT * ≤ T ,

||u ε || L ∞ (0,kT * , Ḣ1 (Ω)) ≤ √ 2 k R||u 0 || Ḣ1 (Ω) .
Energy and entropy estimates. The conservation of mass follows from Proposition 5, but we need to explain how to derive (26), ( 27) and ( 28) from (25). Formally, one has to choose successively φ = -u ε xx + I(u ε ), φ = -u ε xx and φ = G ′ ε (u ε ). Making such a formal computation rigourous is quite standard; details are given in Appendix B for the reader's convenience. Finally, (29) follows from (32).

Proof of Lemma 2. We need to check that the conditions of Leray-Schauder's fixed point theorem are satisfied:

F is compact. Let (v n ) n be a bounded sequence in V and let u n denote F (v n ). The sequence (I(v n )) x is bounded in L 2 (Q),
and so

g n = f ε (v n )∂ x I(v n ) is bounded in L 2 (Q). Estimate (31) implies that u n is bounded in L 2 (0, T * ; H 3 N (Ω)). In particular ∂ xxx u n is bounded in L 2 (0, T * ; L 2 (Ω)) and Equation (30) implies that ∂ t (u n ) is bounded in L 2 (0, T * , H -1 (Ω)). Using Aubin's lemma, we deduce that (u n ) n is pre-compact in V = L 2 (0, T * ; H 2 N (Ω)). F is continuous. Consider now a sequence (v n ) n in V such that v n → v in V and let u n = F (v n ). We have in particular v n → v in L 2 (Q) and, up to a subsequence, we can assume that v n → v almost everywhere in Q. Hence, f ε (v n ) → f ε (v) almost everywhere in Q. We also have that (I(v n )) x converges to (I(v)) x in L 2 (Q), and since |f ε (v n )| ≤ M a.
e., we can show that

g n = f ε (v n )∂ x I(v n ) → f ε (v)∂ x I(v) = g in L 2 (Q).
Next, the compacity of F implies that (u n ) n is pre-compact in the space L 2 (0, T ; H 2 N (Ω)), and so u n converges (up to a subsequence) to U in V . In particular, u n → U in L 2 (Q) and (up to a another subsequence), u n → U almost everywhere in Q. We thus have f ε (u n ) → f ε (U ) in L 2 (Q), and passing to the limit in the equation, we conclude that U = u = F (v) (by the uniqueness result in Proposition 5). Since this holds for any subsequence of u n , we deduce that the whole sequence u n converges to u hence

F (v n ) → F (v) in V as n → ∞
and F is continuous.

A priori estimates. It only remains to show that there exists a constant R > 0 such that for all functions u ∈ V and σ ∈

[0, 1] such that u = σF (u), we have u V ≤ R.
This is where the smallness of T * will be needed.

Using (31), we see that

Ω (u x (t)) 2 dx + ε 2 T * 0 Ω (u xxx ) 2 dx dt ≤ M 2 T * 0 Ω ((I(u) x ) 2 dx dt + Ω (u 0 ) 2 x dx, (34)
and using (21) and the interpolation inequality

u x L 2 (Ω) ≤ C u 1 2 L 2 (Ω) u xx 1 2 L 2 (Ω) , (35)
we get:

T * 0 Ω ((I(u) x ) 2 dx dt ≤ 2C ε T * 0 I(u)(t) 2 2 dt + ε 2M T * 0 (I(u)) xx (t) 2 2 dt ≤ 2C ε T * 0 u x (t) 2 2 dt + ε 2M T * 0 u xxx (t) 2 2 dt ≤ 2C ε T * sup t∈[0,T * ] u x (t) 2 2 + ε 2M T * 0 u xxx (t) 2 2 dt (36)
where C ε only depends on the constant C in (35) and the parameters ε and M . Combining (34) and (36), we conclude that

(1 -M C ε T * ) sup t∈[0,T ] u x (t) 2 2 + ε 4 T * 0 u xxx (t) 2 2 dt ≤ Ω (u 0 ) 2 x dx.
Therefore, choosing T * := 1 2MCε , we get the following estimates

u L ∞ (0,T * ; Ḣ1 (Ω)) ≤ √ 2||u 0 || Ḣ1 (Ω) and u L 2 (0,T * ; Ḣ3 N (Ω)) ≤ 2 √ ε ||u 0 || Ḣ1 (Ω)
Since we also have

Ω u(t) dx = Ω u 0 dx,
we deduce that u V ≤ R for some constant R depending on ε, which completes the proof.

Proof of Theorem 1

As pointed out in the introduction, one of the main difficulties in the proof of Theorem 1 is that the natural energy estimate (31) does not give any information by itself, since E(u) may be negative. Even if Lemma 1 implies that the quantity

αE(u ε ) + β||u ε || L 1 (Ω)
is bounded below by the H 1 norm of u, the mass conservation only allows us to control the L 1 norm of u ε if we know that u ε is non-negative. Unfortunately, it is well known that equation (24) does not satisfy the maximum principle, and that the existence of non-negative solutions of ( 4) is precisely a consequence of the degeneracy of the diffusion coefficient, so that while we can hope (and we will prove) to have lim ε→0 u ε ≥ 0 we do not have, in general, that u ε ≥ 0. Lemma 3 below will show that it is nevertheless possible to derive some a priori estimates that are enough to pass to the limit, provided the initial entropy is finite [START_REF] Cabré | Positive solutions of nonlinear problems involving the square root of the laplacian[END_REF].

Proof of Theorem 1. Consider the solution u ε of (24) given by Theorem 3. In order to prove Theorem 1, we need to show that lim ε→0 u ε exists and solves [START_REF] Passo | On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions[END_REF].

Since we cannot use the energy inequality to get the necessary estimates on u ε , we will need the following lemma: Lemma 3. Let H ε denote the following functional:

H ε (v) = Ω [v 2 x + 2M G ε (v)] dx.
Then the solution u ε given by Theorem 3 satisfies

H ε (u ε (t)) + M 2 t 0 Ω (u ε xx ) 2 dx ds + 1 2 t 0 Ω f ε (u ε )(u ε xxx ) 2 dx ds ≤ H ε (u 0 )e t/2
for every t ∈ [0, T ].

Proof of Lemma 3. Using ( 27) and (28) we see that

Ω G ε (u ε (t)) dx + 1 2 t 0 Ω (u ε xx ) 2 dx ds ≤ Ω G ε (u 0 ) dx + 1 2 t 0 Ω (u ε x ) 2 dx ds
and

Ω (u ε x ) 2 dx + 1 2 t 0 Ω f ε (u ε )[u ε xxx ] 2 dx ds ≤ Ω ((u 0 ) x ) 2 dx + M 2 t 0 Ω (u ε xx ) 2 dx ds.
This implies

H ε (u ε (t)) ≤ H ε (u 0 ) + t 0 H ε (u ε (s))ds,
and Gronwall's lemma yields the desired result.

Sobolev and Hölder bounds. We now gather all the a priori estimates: Using the conservation of mass, Lemma 3 and inequality (9), we see that there exists a constant C independent of ε such that:

sup t∈[0,T ] Ω (u ε x (t)) 2 dx ≤ C, (37)
sup

t∈[0,T ] Ω G ε (u ε (t)) dx ≤ C, (38)
T 0 Ω f ε (u ε )[u ε xxx ] 2 dx dt ≤ C, (39)
T 0 Ω (u ε xx ) 2 dx dt ≤ C. (40)
Next, we note that (37) yields

E(u ε ) ≥ -||u ε || L ∞ (0,T ; Ḣ1/2 (Ω)) ≥ -C||u ε || L ∞ (0,T ;H 1 (Ω)) ≥ -C
and so (26) gives

T 0 Ω f ε (u ε ) u ε xxx -I(u ε ) xx 2 ds dx ≤ C. (41)
Finally, estimates (29), ( 37) and (41) yield that u ε is bounded in

C 1/2,1/8 x,t (Q).
Limit ε → 0. The previous Hölder estimate implies that there exists a function u(x, t) such that u ε converges uniformly to u as ε goes to zero (up to a subsequence). Inequality (40) also implies that u ε ⇀ u in L 2 (0, T ; H 2 (Ω))-weak and Aubin's lemma gives

u ε -→ u in L 2 (0, T ; H 1 (Ω))-strong.
After integration by parts, (25) can be written as

Q u ε φ t -f ε (u ε )[u ε xx -I(u ε )]φ xx -f ′ ε (u ε )u ε x (u ε xx -I(u ε ))φ x dt dx = Ω u 0 (x)φ(0, x) dx
and passing to the limit ε → 0 gives [START_REF] Passo | On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions[END_REF].

Non-negative solution. It only remains to show that u is non-negative. This can be done as in [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF], using (38) (and the fact that f satisfies (5) with n > 1).

L ∞ a priori estimate. Finally, (37) and Sobolev's embedding implies that there exits a constant M 0 depending only on

||u 0 || H 1 (Ω) such that 0 ≤ u(t, x) ≤ M 0 . (42)
Choosing M > M 0 , we deduce that u solves (4).

Proof of Theorem 2

In order to get Theorem 2, we need to derive the following corollary from Theorem 1.

Corollary 1. The solution u constructed in Theorem 1 satisfies for all φ ∈ D((0, T ) × Ω),

Q uφ t dt dx + P f (u)[u xx -I(u)] x φ x dt dx = 0 ( 43
)
where P = {(x, t) ∈ Q : u(x, t) > 0, t > 0}.

Proof. In view of the proof of Theorem 1, u is the uniform limit of a subsequence of (u ε ) ε>0 where u ε is given by Theorem 3. Since u ε satisfies (25), it is thus enough to pass to the limit in this weak formulation as ε → 0 in order to get the desired result. Let

h ε denote f ε (u ε )[u ε xx -I(u)] x . Estimates (41) and (42) imply Q h 2 ε dx dt ≤ C. (44)
In other words, (h ε ) ε is bounded in L 2 (Q). Hence, up to a subsequence,

h ε ⇀ h in L 2 (Q)-weak.
Furthermore, we recall that there exists a continuous function u(x, t) such that u ε converges uniformly to u as ε goes to zero (up to a subsequence).

Passing to the limit in (25), we deduce that the function u satisfies

Q uφ t dt dx + Ω hφ x dt dx = 0.
We now have to show that

h = 0 in {u = 0}, f (u) u xx -I(u) x in P = {u > 0}.
First we note that for any test function φ and η > 0, we have

T 0 {u≤η} f ε (u ε ) u ε xxx -I(u ε ) x φ dx dt ≤ C(φ) f ε (3η/2) 1/2 T 0 {u≤η} f ε (u ε ) u xxx -I(u) x 2 dx dt 1/2
for ε small enough (so that |u ε -u| ≤ η/2). Inequality (41) thus implies lim sup

ε→0 T 0 {u≤2η} f ε (u ε ) u ε xxx -I(u ε ) x φ dx dt ≤ C(φ)f (η/2) 1/2 .
We deduce (since f (0) = 0)

h = 0 on {u = 0}. (45)
Next, (44) yields (for ε > 0 and η small enough)

{u>2η} |u ε xxx -I(u ε ) x | 2 dx ds ≤ C(η). This implies that, if Q η denotes {u > 2η}, (u ε xxx -I(u ε ) x ) is bounded in the space L 2 (Q η ). Hence, we can extract from (u ε xxx -I(u ε ) x ) ε>0 a subsequence converging weakly in L 2 (Q η ). Moreover, remark that Q η is an open subset of Q (recall that u is Hölder continuous) and u ε xxx -I(u ε )
x converges in the sense of distributions to u xxx -I(u) x (use the integral representation for I(•)). We thus conclude that,

u ε xxx -I(u ε ) x ⇀ u xxx -I(u) x in L 2 (Q η ). This yields h = f (u) u xxx -I(u) xx in {u > 0}
which concludes the proof of Corollary 1.

We now turn to the proof of Theorem 2.

Proof of Theorem 2. When u 0 does not satisfy (9), we lose the L 2 (0, T ; H 2 (Ω)) bound on u ε , and the previous analysis fails. However, we can introduce u δ 0 = u 0 + δ which satisfies [START_REF] Cabré | Positive solutions of nonlinear problems involving the square root of the laplacian[END_REF]. Theorem 1 then provides the existence of a non-negative solution u δ of [START_REF] Passo | On a fourth-order degenerate parabolic equation: global entropy estimates, existence, and qualitative behavior of solutions[END_REF]. In view of Corollary 1, u δ satisfies:

Q u δ φ t dt dx + P f (u δ )[u δ xx -I(u δ )] x φ x dt dx = 0. (46)
Since u δ is non-negative, the conservation of mass gives a bound in the space L ∞ (0, T ; L 1 (Ω)) and allows us to make use of the energy inequality: Indeed, using [START_REF] Frankel | On a free boundary problem associated with combustion and solidification[END_REF] and Lemma 1 we see that there exists a constant C independent of δ such that

||u δ || L ∞ (0,T ;H 1 (Ω)) ≤ C and T 0 Ω f (u δ ) u δ xxx -I(u δ ) x 2 ds dx ≤ C. (47)
We now define the flux

h δ = f (u δ ) u δ xxx -I(u δ )
x . Inequality (47) implies that h δ is bounded in L 2 (Q), and that there exists a function h ∈ L 2 (Q) such that

h δ ⇀ h in L 2 (Q)-weak.
Proceeding as in the proof of Theorem 1, we deduce that u δ is bounded in C 1/2,1/8 (Ω×(0, T )) and that there exists a function u(x, t) such that u δ converges uniformly to u as δ goes to zero (up to a subsequence). We can now argue (with minor changes) as in the proof of Corollary 1 and conclude.

A Proof of Proposition 5

Our goal here is to prove the existence of a weak solution of u t + (au xxx ) x = (ag) x .

We first prove the following proposition.

Proposition 6. For all h ∈ H 1 (Ω), there exists v ∈ V 0 := H 1 ∩ H 3 N such that for all φ ∈ D( Ω),

- Ω v -h τ φ dx + av xxx φ x dx = Ω agφ x dx. (48)
In particular,

Ω v dx = Ω h dx, 1 2 Ω v 2 x + τ Ω av 2 xxx ≤ 1 2 Ω h 2 x + τ Ω agv xxx . (49)
Proof. In order to prove this proposition, we have to reformulate the equation. More precisely, instead of choosing test functions φ ∈ D( Ω), we choose φ = -ψ xx + Ω ψ dx where ψ is given by the following lemma We check that it is continuous and coercive:

|A(v, w)| ≤ v x 2 w x 2 + M τ v xxx 2 w xxx 2 + v 1 w 1 ≤ C v V0 w V0 A(v, v) ≥ Ω [(v x ) 2 + ετ (v xxx ) 2 ] dx + Ω v dx 2 ≥ v 2 V0 .
We now consider the following linear form L in V 0 : for all w ∈ V 0 ,

L(w) = Ω h x w x dx + Ω h dx Ω w dx + τ Ω agw xxx dx
Since 0 ≤ a ≤ M and g ∈ L 2 , L is continuous as soon as h ∈ H 1 (Ω). Lax-Milgram theorem thus implies that there exists v ∈ V 0 such that (50) holds true for all w ∈ V 0 . Eventually, remark that conservation of mass and (49) are direct consequences of (48). The proof of Proposition 6 is now complete.

We can now prove Proposition 5.

Proof of Proposition 5. For any τ > 0, we consider N τ = ⌈ T τ ⌉. We then define inductively a sequence (u n ) n=0,...,Nτ of V 0 as follows: u 0 = u 0 and u n+1 is obtained by applying Proposition 6 to h = u n . We then define u τ : [0, N τ τ ) × Ω as follows:

u τ (t, x) = u n (x) for t ∈ [nτ, (n + 1)τ ).

We have Ω u τ (t, x) dx = Ω u 0 (x) dx for all t. We also derive from (49) that we have In particular, (u τ ) τ is bounded in L ∞ (0, T ; H 1 (Ω)) and (S τ u τ -u τ ) τ is bounded in L 2 (0, T -τ ; H -1 (Ω)) where S τ v(t, x) = v(t + τ, x). We derive from [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]Theorem 5] that (u τ ) τ is relatively compact in C(0, T ; L 2 (Ω)).

We now have to pass to the limit in (48). Since (u τ xxx ) τ is bounded in L 2 (Q) and we can find a sequence τ n → 0 such that u τn → u in C(0, T, L 2 (Ω)) and u τn xxx → u xxx in L 2 (Q). This is enough to conclude. We next explain how to get (32). Sobolev's embedding imply that there exists a constant K (depending on u L ∞ (0,T ;H 1 (Ω)) ) such that

|u(x 1 , t) -u(x 2 , t)| ≤ K|x 1 -x 2 | 1/2

Figure 1 :

 1 Figure 1: A viscous thin film submitted to an electric field E and gravity g

Lemma 4 .= v xxx 2 L 2 +Ω

 422 For all φ ∈ D( Ω), there exists ψ ∈ D( Ω) such that-ψ xx + Ω ψ dx = φ.Hence, we considerV 0 := H 1 ∩ H 3 N equipped with the norm v 2 V0 ( v dx) 2 and we look for v ∈ V 0 such that for all ψ ∈ D( Ω), Ω v x ψ x dx + τ av xxx ψ xxx dx + agψ xxx dx. (50)We thus consider the bilinear form A in V 0 defined as follows: for all v, w ∈ V 0 , A(v, w) = Ω v x w x dx + τ av xxx w xxx dx +

Ω (u τ x ) 2

 2 (T, x) dx + T 0 Ω a(u τ xxx ) 2 (t, x) dtdx ≤ Ω ((u 0 ) x ) 2 dx + T 0 Ω (agu xxx )(t, x) dtdx

for all x 1 , x 2 ∈ Ω and for all t ∈ (0, T ). Since u satisfies u t = h x with h ∈ L 2 (Q), it is a fairly classical result that Hölder regularity in space implies Hölder regularity in time. More precisely, we have (see [START_REF] Bernis | Higher order nonlinear degenerate parabolic equations[END_REF], Lemma 2.1 for details): Lemma 5. There exists a constant C such that for all x 1 , x 2 in Ω and all t 1 ,

The proof of Proposition 5 is thus complete.

B Proof of (26), ( 27) and (28)

We have to derive (26), ( 27) and ( 28) from (25). Using the fact that u ε lies in C([0, T ], H 1 (Ω)), we first state the following lemma

The proof of such a lemma is fairly classical. It relies on mollifiers that are decentered in the time variable. More precisely, one considers a smooth even function ρ : R → [0, 1] compactly supported in [-1, 1] and such that ρ = 1. Then for α > 0 and δ ∈ R, one can define

If now a function f is defined in [0, T ], it can be extended by 0 to R; in other words, it can be replaced with f 1 Ω where 1 Ω (x) = 1 if x ∈ Ω and 1 Ω (x) = 0 if not; then the convolution product in R:

Consider a smooth function ρ such as in the proof of Lemma 6. Consider α > 0 and define ρ α (x) = αρ( x α ). Then consider θ(x, t) = ρ α (x)ρ α (t).

We next apply Lemma 6 with φ = v ⋆ θ where v is chosen to be successively u ε xx ⋆ θ, I(u ε ⋆ θ) and G ′ ε (u ε ⋆ θ). After direct computations, we can let α → 0 and get the desired estimates.


Emmanuel Witrant
email: emmanuel.witrant@gipsa-lab.inpg.fr.

Pangun Park
email: pgpark@ee.kth.se.

Mikael Johansson ⋆⋆
email: mikaelj@ee.kth.se.

Time-delay estimation and finite-spectrum assignment for control over multi-hop WSN

Keywords: Network controlled systems, wireless networks, multi-hop networks, time-delay estimation, predictive control, finite spectrum assignment

Recent technological developments and the need for global control strategies to meet with stringent performance requirements rendered the use of wireless sensor networks (WSN) of prime importance for today's automation. Such communication devices heavily affect the information transport, especially when a multi-hop configuration is introduced to minimize the energy consumption of the network. The associated time-delays are strongly varying and necessitate the design of a dedicated control approach for remote regulation, as proposed in this chapter. First, the communication process is analyzed and illustrated with the Breath protocol, which results in a description of the delay in the frequency domain. Such analysis motivates the use of a CUMSUM Kalman filter to estimate the delay from round-trip-time (RTT) measurements while selecting the desired frequencies. This estimation then allows for the setup of a state-predictor that ensures a finite spectrum assignment (FSA) on the closed loop spectrum. The efficiency of this approach and the impact of packet losses is finally investigated on an inverted pendulum using experimental network measurements.

Introduction

There are major advantages in terms of increased productivity and reduced installation costs in the use of wireless communication technology in industrial control systems [START_REF] Willig | Wireless technology in industrial networks[END_REF][START_REF] Park | Protocol design for control applications using wireless sensor networks[END_REF]. Several wireless technologies are currently reaching the plant and workshop floor addressing a range of applications. A multi-hop wireless network is one such technology that offers great flexibility, but also challenges when it comes to predictability.

Wireless sensor networks (WSN) are used in emerging applications in industrial automation and dedicated algorithms, such as randomized multi-hop routing protocols [START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF], were recently shown to be useful in these applications. From the perspective of designing control systems utilizing wireless multi-hop networks, one of the most significant characteristics is the need of hiding the system complexity through suitable abstraction of the networks and their nodes. One simple but important way of doing this is by interpreting the network as a communication channel with time-varying delay.

Experimental results on telerobotics over TCP/IP networks illustrate the importance of varying transmission delays in the control loop [START_REF] Lelevé | Telerobotics over IP networks: towards a lowlevel real-time architecture[END_REF], where the authors combine a generic estimation/prediction control scheme with a buffering of the received packets, which compensates for the delays jitter. While the buffer may decrease the transient performances of the closed-loop system [START_REF] Witrant | On the use of state predictors in networked control systems[END_REF], this work motivates the use of estimation/prediction schemes to achieve robust and efficient control over networks. Another interesting experimental result is provided in [START_REF] Ploplys | Closed-loop control over wireles networks[END_REF], where a peer-to-peer wireless network (802.11b) and a UDP protocol are used to transmit the sensors and actuators signals. An inverted pendulum is controlled through this network using LQR techniques without compensating the delay (classical state feedback), which results in relatively large oscillations in the system response (complex poles induced by the delay) even for small latencies (average value less then 3.6 ms). The packet losses (significant effect above 15%) are efficiently compensated using a dynamic adjustment of the sampling interval. These experimental results illustrate the fact that, when the transmissions are fairly deterministic and the control/sensors signals have a significant effect on the network load, simple solutions can be used to compensate for the packet losses and timing issues while the delay effects need a more careful consideration.

The main focus of this chapter is a new scheme for control over wireless multi-hop networks based on time-delay and finite spectrum assignment (FSA). Time-varying delays for the transmission of sensor and control data over the wireless network are caused by a randomized multi-hop routing protocol called Breath [START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF]. Breath is designed for control applications using WSN and it ensures a desired packet delivery and delay probabilities while minimizing the energy consumption of the network. We also consider the modification of Breath by using an Automatic Repeat reQuest (ARQ) mechanism [START_REF] Leon-Garcia | Counmmunication Netorks[END_REF] to improve the packet reception rate (reliability). The characteristics of the routing protocol together with lower-layer network mechanisms give rise to a delay process with high variance and stepwise changing mean. Estimation of the transport delay over multi-hop networks is particularly difficult, since detailed modeling from first principles is prevented by the high complexity and multitude of traffics over the network. Here we use the low-order model proposed in [START_REF] Jacobsson | Some modeling and estimation issues in traffic control of heterogeneous networks[END_REF] to set a predictive control scheme with a delay estimator. The estimator is based on a Kalman filter with a change detection algorithm [START_REF] Jacobsson | Some modeling and estimation issues in traffic control of heterogeneous networks[END_REF]. It is able to track the delay mean changes but efficiently attenuate the high frequency jitter. Compared to similar works such as [START_REF] Seuret | Implementation of an internet-controlled system under variable delays[END_REF], where an output feedback scheme uses the maximal value of the delay and GPS synchronization for remote control over Internet, we explicitly compensate for the time-varying delay effect. The control scheme is analyzed and its implementation detailed. Network data from an experimental setup, where the randomized multi-hop routing protocol [START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF] is implemented on sensor nodes, are used to illustrate the efficiency of the approach. This work is an extended version of [START_REF] Witrant | Control over wireless multi-hop networks[END_REF], with additional insights on the delay estimation and WSN aspects.

The outline of the paper is as follows. Section 6.2 defines the considered problem of control over a wireless multi-hop network. The wireless network and its randomized routing protocol is described in Section 6.3. The timedelay estimator is presented in Section 6.4. Section 6.5 presents the controller and its delay compensator. Experimental results are given in Section 6.6, where the tradeoff between lossless communications and time-delays is discussed.

Problem description

The problem considered is depicted in Figure 6.1, where a linear system is remotely controlled over a wireless multi-hop network. We suppose that roundtrip time or end-to-end delay measurements (i.e., between the sensors and controller, which also implies a proper clock synchronization) are available. A predictive control law uses the filtered measurements to compensate for variations induced by the multi-hop network, which induces communication delays (time-varying average latency and jitter), packet losses and communication bandwidth constraints. More precisely, we consider the remote control of linear systems that write as:

ẋ(t) = Ax(t) + Bv(t), x(0) = x 0 (6.1)
y(t) = Cx(t) (6.2)
where x ∈ R n is the internal state, v ∈ R is the system input, y ∈ R m is the system output, and A, B, C are matrices of appropriate dimensions. The pairs (A, B) and (A, C) are assumed to be controllable and observable, respectively, but no assumption is made on the stability of A. The delays introduced by the network dynamics are defined as:

• τ 1 (t): delay between the sensors and the controller,

• τ 2 (t): delay between the controller and the actuator.

The actuator is event-driven and uses the last received signal. Packet losses are included in this description by considering that the last received and most recent controller output u(•) is taken as the input for the system, which is then:

w(t k ) = u(t k -τ 2 (t k )) if packet k is received w(t k ) = w(t k-1 ) else (6.3)
where w(•) is the signal received by the actuator from the network, u(•) is the controller output and v(t) = w(t k ) ∀ t ∈ [t k , t k+1 ). The index k is used to differentiate the discrete data packets from the continuous physical variables. The chronology of the signals is indicated in the transmitted data in order to reconstruct their history on the receiver site. The receivers consider only the most recent information, meaning that a packet is also considered as lost if a more recent one is received first.

Wireless multi-hop networks

In many scenarios of relevant interest, wireless multi-hop networks do not have fixed communication paths [START_REF] Karaki | Routing technques in wireless sensor networks: A survey[END_REF], but the end-to-end path followed by packets happens according to a dynamic selection of hops [START_REF] Greuen | Adaptive sleep discipline for energy conservation and robustness in dense sensor networks[END_REF][START_REF] Xu | Geography-informed energy conservation for ad hoc routing[END_REF]. Indeed, it is often impossible to build fixed routing tables for these networks due to the time-varying communication channel and network topology. In such networks, a packet is routed to a particular next-hop node as long as the node selected for a hop ensures progress toward destination. Furthermore, a major technique to ensure power savings and longer network lifetime is to turn off a node whenever its presence is not strictly required for the correct operation of the network [START_REF] Park | Protocol design for control applications using wireless sensor networks[END_REF][START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF][START_REF] Xu | Geography-informed energy conservation for ad hoc routing[END_REF][START_REF] Chen | Span: An energy-efficient coordination algorithm for topology maintenance in ad hoc wireless networks[END_REF]. In this case, each node goes to sleep for a random amount of time depending on traffic and network conditions, which means that the network topology is changing randomly. Consider the multi-hop network reported in Figure 6.1, where the transmitter to receiver path is composed of a time-varying number of hops h(t). Therefore we can consider the following sources of delay:

• the time α i to wait before sending a data packet, which is typically a random variable; • the time F to forward a data packet once the connection with the next node is set up, which is constant and contains the propagation and transmission delay; • possibly the time β i induced by an ARQ mechanism, which retransmits a packet for a predefined number of time if no acknowledgement is received.

The end-to-end delay can then be expressed as follows:

τ (t) = h(t)F + h(t) i=1 (α i + β i ) (6.4)
For the Breath protocol [START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF] used to set the experimental setup, the routing, medium access control and sleeping discipline are randomized in such a way that the energy consumption of the communication network is minimized while meeting the constraints on packet delay and reliability. Since it is representative of a class of protocol for wireless automation, we take it as a reference for the experimental results. With this protocol, α i depends on the exponentially distributed wake up rate, which average value µ c,i is constrained by a function of the traffic rate, channel condition and application requirements. Note that, since h(t), α i and β i are random variables, the end-to-end delay clearly results in a random variable.

If one looks at the outcomes of (6.4), there are three possible behaviors:

1. Low frequency delays, due to changes in the number of hops h(t). The end-to-end number of hops is expected to change when an obstacle changes the channel conditions for long time (i.e. the shadow fading of the wireless propagation). Since shadowing changes with low frequency, this causes low frequency changes in the delay. 2. Middle frequency delays, which depend on the random selection of the next awaken hop for a fixed number of clusters (related to µ c,i in the example considered). This is induced by the fact that a transmitting node has to wait for some node to wake up in the forwarding region. For example, in the randomized protocol, this phenomenon is directly caused by the wake up rate µ c,i , which is higher than the rate with which the shadow fading changes. 3. High frequency delays, related to the transmissions between nodes and packet loss probability. Since the wireless transmission may be lossy due to fast fading (which causes the execution of the ARQ) there might be highly time-varying delays. Hence, this delay is directly related to the packet loss probability. Note that the frequency of the fast fading is very high in comparison with the shadow fading. Example 1. Considering the WSN multi-hop network presented in Section 6.6, which describes the specific experimental setup, time-delay measurements are obtained for a three-hops configuration. The power spectral density of these delays is obtained thanks to the periodograms computation and depicted in Figure 6.2. We can clearly notice the dominent low-frequency component of the delays (0 dB/rad/sample), followed by the medium (two peaks above -15 dB/rad/sample) and high (below -18 dB/rad/sample, close to white noise) frequency behaviors.

The measurements can also be used to investigate the class of models that can be associated with such delays. As an example, a spectral model (empirical transfer function estimate) and a time-domain model (auto regressive moving average -ARMA) are presented in Figure 6.3. The ARMA model is set with a second order polynomial on the delay and a first order polynomial on the white noise. While providing for a rough approximation of the signal, we can see that the ARMA model gives a reasonable approximation of the amplitudes at the three frequency levels. This motivates the structure of the delay estimation presented in the next section.

Time-delay estimator

A critical component of the proposed control strategy is the delay estimator, which output is used to set the time-varying horizon of the predictive controller. Examining the real delay traces from our network (introduced in the previous example), shown in dashed lines in Figure 6.4, we see that the delay exhibits a large variability and persistent changes in the mean. Although these changes can be traced back to a particular feature of our protocol (discussed in the previous section) similar effects have been reported for a wide range of network quantities and network technologies [15]. Since the controller-actuator latency is not immediately available to the controller, but has to be deduced from the (delayed information about the) controlleractuator round-trip time, we will look for a delay estimator that allows to smooth out the fast variations while being able to react quickly to the changes in mean.

To this end, we propose to use a combination of Kalman filter and CUM-SUM change-detection (cf., [START_REF] Jacobsson | Some modeling and estimation issues in traffic control of heterogeneous networks[END_REF][START_REF] Gustafsson | Filtering and Change Detection[END_REF]). The underlying signal model is:

τ k+1 = τ k + v k r k = cτ k + w k
where τ k is the one-way latency, r k is the measured time, c is the number of delays included in r k (i.e., 1 if the end-to-end delay of the channel is measured and 2 if it is the round trip time), w k is a white-noise sequence with covariance R k and: v k = 0 with probability 1 -q, ν with probability q, where Cov(ν k ) = q -1 Q k , models the abrupt changes in the mean. The associated Kalman filter updates are given by:

τk+1 = τk + K k (r k -cτ k ), with K k = cP k-1 R k + c 2 P k-1 , P k = P k-1 - c 2 P 2 k-1 R k + c 2 P k-1 + Q k
The tracking ability of the filter is proportional to the variance of the process noise Q k , but increasing Q k also decreases the smoothing effect that we are looking for.

To improve performance, we run a CUMSUM detector in parallel with the Kalman filter. The detector looks at the prediction errors ε k = r k -cτ k , sums up the negative and positive prediction errors:

g + k = max(g k-1 + ε k -κ, 0) g - k = max(g k-1 -ε k -κ, 0)
and monitors if g + k or g - k exceeds a threshold g thr . Here, the parameter κ is a design parameter that specifies the negative drift (see [START_REF] Gustafsson | Filtering and Change Detection[END_REF] for design guidelines). At these threshold crossings, Q k is set to a large value, g + and g -are reset to zero, τk is set to c -1 r k and the Kalman filter iterations resumed. The combined Kalman-filter and CUMSUM test allows to achieve both fast tracking of the changes in the mean and good suppression of the high-frequency variations, as illustrated in Figure 6.4 (solid line).

Predictive control approach

The control approach proposed here is based on the state predictor with a time-varying horizon first proposed in [START_REF] Olbrot | On controllability of linear systems with time delays in the control[END_REF]. The main advantage of such method compared to other delay compensation strategies [START_REF] Richard | Time delay systems: An overview of some recent advances and open problems[END_REF] is to allow for a finite spectrum assignment on systems with time-varying delays [START_REF] Manitius | Finite spectrum assignment problem for systems with delays[END_REF][START_REF] Artstein | Linear systems with delayed control: a reduction[END_REF]. An overview on the use of state predictors in networked control systems is provided in [START_REF] Witrant | On the use of state predictors in networked control systems[END_REF], where it is compared with more classical control approaches.

Ideal delay compensation

We use the results established in [START_REF] Witrant | Remote stabilization via time-varying communication network delays: Application to TCP networks[END_REF] and subsequent works, where the delay induced by the network dynamics is explicitly taken into account in the control setup. For simplicity and to focus on the network effect on the closed-loop system, we suppose that the full delayed state3 x(t-τ 1 ) is available to set the control law. This hypothesis can be easily removed using an observer-based control setup as the one proposed in [START_REF] Witrant | Remote output stabilization under two channels time-varying delays[END_REF]. When there is no packet loss and the delay is fully known (and can be predicted), the main theorem in [START_REF] Witrant | Remote stabilization via time-varying communication network delays: Application to TCP networks[END_REF] can be extended to the case where both communication channels experience a delay as follows.

Theorem 1. Consider the system:

ẋ(t) = Ax(t) + Bu(t -τ 2 ), x(0) = x 0 (6.5)
with (A, B) a controllable pair. Assume that the network dynamics is such that the following holds for τ 1 and τ 2 :

A1) τ 1,2 ∈ C(R + , [0, τ max ]) ∀t, A2) τ1,2 < 1 ∀t ≥ t 0 .
Then the feedback control law:

u(t) = -Ke Aδ e Aτ1 x(t -τ 1 ) + e At t+δ t-τ1 e -Aθ Bu(θ -τ 2 (θ))dθ (6.6) δ(t) = - λ 1 -dτ 2 (ζ)/dζ δ + dτ 2 (ζ)/dζ + λτ 2 (ζ) 1 -dτ 2 (ζ)/dζ (6.7)
with ζ = t + δ, λ a positive constant, δ(0) = δ 0 , ensures that the system trajectories x(t) converge to zero.

Proof. (Outline) First note that (6.6) is obtained by setting the feedback law as:

u(t) = -Kx(t + δ(t))
where x(t + δ(t) is obtained from the fundamental of (6.5) as:

x(t + δ) = e Aδ e Aτ1 x(tτ 1 ) + e At t+δ t-τ1

e -Aθ Bu(θτ 2 (θ))dθ

Noting that with (A1) -(A2) the dynamics (6.7) ensures the exponential convergence towards zero of δ(t) -τ (t + δ(t)) the closed-loop dynamics is set as:

d dt x(t + δ) = (A -BK)x(t + δ)
Such dynamics is stable if (A -BK) is Hurwitz and the convergence of x(t) is a direct consequence of the one of x(t + δ).

⋄ ⋄ ⋄

If the communication delay is perfectly known and satisfies A1-A2 (which correspond to a lossless network description), then the closed-loop system behavior can be inferred from the time-shifted dynamics:

dx dζ = (A -BK)x(ζ)
and its spectrum coincides with the spectrum of the matrix:

A - 1 + δ τ (t + δ) e -Aδ BK
Considering the worst-case delay properties to set the closed loop spectrum, which correspond to a maximum delay τ max and a maximum delay variation τmax , the feedback gain has to be chosen such that the matrix:

A - 1 1 -τmax e -Aτmax BK (6.8)
is Hurwitz.

Use of the delay estimator

The control law is set using the delay estimator described in section 6.4. Both delays τ 1 and τ 2 are then replaced by τ in the control setup description introduced in Theorem 1. This is motivated by the fact that the predictor architecture in used to compensate for the time-variations of the latency (i.e., the slow and medium dynamics of the delay), which are estimated with the CUMSUM Kalman filter. It is reasonable to suppose that the average values of the delays in both communication channels are close and provided by the estimator. The control law is then given as:

u(t) = -Ke Aδ e Aτ x(t -τ 1 ) + e At t+δ t-τ
e -Aθ Bu(θτ (θ))dθ (6.9)

The predictor horizon is computed using the actual values of the estimated delay and delay variation directly instead of the predicted ones (at time ζ).

We then approximate δ with δ, which dynamics is given by:

δ(t) = τ (t) + λ(τ (t) -δ(t))
1 -τ (t) (6.10)

Control algorithm

Special care has to be taken in the computation of the control law (6.9)-(6.10), especially concerning the integral term used in the predictor's part. Indeed, the discretization of the integral may introduce some numerical instabilities, as illustrated in [START_REF] Van Assche | Some problems arising in the implementation of distributed-delay control laws[END_REF] and explained in [START_REF] Mondié | Finite spectrum assignment of unstable time-delay systems with a safe implementation[END_REF]. The effect of these instabilities is further increased if the communication channel experiences packet losses and bandwidth limitation. A solution to this problem is to include the discretized integral in a dynamic control law, as described in this section. First, considering the integral term:

I(t) . = e At t+ δ t-τ
e -Aθ Bu(θτ (θ))dθ computed at the time instant denoted by k, we have that:

I k = e Aτ k τk + δk 0 e -Aι Bu(µ + t k -τk -τ (µ + t k -τk ))dµ Introducing the function f (µ) . = e -Aµ Bu(µ + t k -τk -τ (µ + t k -τk )
) and considering that the delay evolution can not be predicted, the integral term is computed with:

f (µ) = e -Aµ Bu(min(µ + t k -τk -τ (min(µ + t k -τk , t k )), t k ))
to ensure the causality of the predictor computation.

The second step is to discretize the integral term, which is done according to the guidelines proposed in [START_REF] Mondié | Finite spectrum assignment of unstable time-delay systems with a safe implementation[END_REF] for the single delay case. More precisely, the integral term is discretized with a trapezoidal rule and included into the dynamic controller:

     ż(t) = -az(t) + bKe A( δk +τ k ) x(t -τ 1 ) + t s 2 n k i=1 f ((i -1)t s ) + f (it s ) u k = z(t k ) (6.
11) where n k . = ( δk + τk )/t s , a and b are some positive design parameters, t s is the sampling time set by the controller clock and n k is the number of steps used to discretize the integral, which is time-varying. Note that we chose the controller sampling time as the integrand step, since it sets the buffering of the past estimated delays and controller outputs.

Finally, an upper bound τthr is set on the maximum estimated delay variation that can be used for the computation of δ(t) since the delay estimator may have some fast variations, especially when it is triggered by the alarm signal. The dynamics of the predictor horizon is then set with:

δ(t) = τk + λ(τ k -δ(t)) 1 -τk (6.12)
where:

   τk = τk -τk-1 t s if τ < τthr τk =

τthr else

The issues considered in this section are all critical for a safe implementation of the predictor with a time-varying horizon, especially to control a system over the network considered. They are summarized with the following result.

Result 6.5.1 Using the delay estimation provided by the CUMSUM Kalman filter described in section 6.4, the proposed predictive control approach is set with the dynamic controller (6.11), where the predictor horizon is set by the dynamics (6.12) and the controller gain is such that the closed-loop matrix (6.8) is Hurwitz.

Experimental and simulation results

In this section, the network setup is presented along with some experimental results, which are used to estimate the efficiency of the proposed control setup using appropriate simulation tools. More details on the experimental setup and protocol specifications can be found in [START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF].

Communication protocol

We use the Breath protocol proposed in [START_REF] Park | Breath: a self-adapting protocol for wireless sensor networks in control and automation[END_REF] and briefly discussed in Section 6.3. Breath is designed for WSNs where source nodes attached to a plant must transmit information via multi-hop routing to a sink. Breath ensures a desired packet delivery and delay probabilities while minimizing the energy consumption of the network. The design approach relies on a constrained optimization problem, whereby the objective function is the energy consumption and the constraints are the packet delay and reliability. The optimal working point of the protocol is achieved by a simple algorithm, which adapts to traffic variations and channel conditions with negligible overhead. The protocol has been implemented and experimentally evaluated on a testbed with off-the-shelf wireless sensor nodes [START_REF]Tmote Sky Data Sheet[END_REF], and it has been compared with a IEEE 802.15.4 standard solution [26]. Analytical and experimental results show that Breath is tunable and meets delay and reliability requirements. Breath exhibits a good distribution of the working load, thus ensuring a long lifetime of the network. In addition, since this protocol does not include any specific error control mechanism that can ensure a reliable communication, we also introduce a stop-and-wait ARQ algorithm to provide for a more reliable data transfer service. ARQ combines error detection and retransmission to ensure that the data is delivered accurately to the user, despite the errors that may occur during transmission [START_REF] Leon-Garcia | Counmmunication Netorks[END_REF]. However there is packet loss due to the limited number of retransmissions with ARQ mechanism (2 in our experiment). In the following, we consider the Breath protocol without and with ARQ mechanism.

Network setup and experimental measurements

Our experiment is based on the Breath protocol using Tmote nodes [START_REF]Tmote Sky Data Sheet[END_REF], which are equipped with the radio controller Chipcon CC2420 operating at 2.4 GHz and supporting the IEEE 802.15.4 standard. The source is located at 30 m from the destination node, which is connected to the computer and generates 20 pckts/s. Eighteen nodes are placed at regular intervals (1.58 m) along a straight line without obstacles. Surrounding objects are static, with minimal time-varying changes in the wireless channel due to multi-path fading effects. At each location the sensor node is placed in the same orientation to avoid non-linearities in the transmission pattern. Figures 6.5 and 6.6 show the average end-to-end delay and reliability of the Breath protocol without and with ARQ, respectively. In the figure, note that "NH" refers to the number of hops of the network. We recall that the number of hops is set depending on the channel condition, traffic load and application requirements of the network. It is clear to observe that both the average end-to-end delay and reliability increase as the number of hops increases. In addition, ARQ mechanism increases the end-to-end delay because of the waiting time to receive an acknowledgement. We remark here that there is tradeoff between the packet delay and reliability [START_REF] Park | Protocol design for control applications using wireless sensor networks[END_REF]. Furthermore, increasing the number of hops gives an outstandingly bad reliability if more than four hops are used.

Physical system and control setup

The closed-loop system is presented as a block diagram in Figure 6.7, which depicts the plant, network, controller, and their interconnection. As an example, the linear system considered is similar to the T-shape ECP inverted pendulum presented in [START_REF] Witrant | Remote stabilization via time-varying communication network delays: Application to TCP networks[END_REF], an unstable open-loop system with no-minimum phase. Its dynamics is set with the state-space matrices:

A =     0 1 0 0 0.80 0 12.56 0 0 0 0 1 -2.42 0 -8.33 0     , B =     0 4.57 0 0.38    
The eigenvalues of A are λ 1,2 = 0.74 ± 2.08i, λ 3,4 = -0.74 ± 2.08i. The choice of an unstable open-loop plant emphasizes the efficiency of the controller to compensate for the transmission perturbations and the performance limitations.

Experimental delay measurements provide for τ 1 and τ 2 by splitting the initial measured data into two sets, one every other sample being attributed to each set. This means that the channels are not symmetric but belong to the same network (same global configuration). Only the last received packets are considered and the packet losses are handled as described in (6.3). The bandwidth limitation is set by zero-order hold with a 50 ms sampling time. The controller includes:

• the time-delay estimator (based on end-to-end delay measurements); • a reference trajectory block that sets the control objective;

• the predictive control setup. The CUMSUM Kalman filter is set with R = 10 000, Q = 1, ν = 0.01 and g thr = 35. The reference trajectory, for the pendulum angle x(3), is a filtered square signal. The predictive control is set with λ = 10, a = b = 100 (the introduced controller dynamics is equivalent to a low-pass filter) and K = [0.48 0.36 2.20 0.57] sets the poles of (6.8) to (-3 + 0.1i; -3 -0.1i; -4; -4).

Simulation results

Figures 6.8 and 6.9 (without and with ARQ) present the actual and estimated time-delays (top), the reference and the system output (middle), and the control output (bottom). In both cases the transmissions are very noisy, in terms of jitter and packet losses (up to 60% of the signal). These figures illustrate the controller efficiency and the proper handling of the high frequency noise induced by the numerical integration problem. The effect of the number of hops on the system response clearly appears in both cases, particularly when a change in the reference sign occurs when there are four hops. Indeed Figure 6.8 depicts a situation where the bad quality of the communication channel introduces some oscillations, which are reduced as soon as the number of hops is set back to two. The benefit of ARQ appears at the beginning of the simulations but is discussable when the decrease in the average reliability is not compensated by the communication algorithm. Comparing the different simulation results, it appears that the system oscillations are induced by the low reliability rather than high endto-end delay (i.e., comparing the different 4-hops cases). This is due to the fact that the control setup explicitly compensate the delay effects while the robustness with respect to packet losses is not specifically addressed in the feedback design.

These simulations also show that an increase of the number of hops has few effects on the system response when it is established in a steady state. A control-oriented multi-hop wireless network protocol could efficiently use this information to minimize the emission energy (increase the number of hops) according to some measurements of the physical state while ensuring satisfying performances. Fig. 6.9 Delays, system output and control signal with ARQ.

Conclusions

This chapter addressed the problem of remote control over a multi-hop wireless network. The constraints of such a network on the communication channel, such as time-varying delays (jitter), packet losses and bandwidth limitations are considered. The time-delay induced by the network is estimated thanks to a CUMSUM Kalman filter that provides for a fast tracking of the mean and good suppression of the high frequency variations of the delay (jitter). This estimator is used to design a dynamic predictive control law that compensates for the estimated part of the delay explicitly. The causality of this controller is set according to the estimator and specific considerations on safe implementation are detailed. A necessary condition on the assigned closed-loop spectrum is provided but the particular poles location is not specified, which provides for a degree of freedom for optimal or robust control approaches. Based on experimental delay and packet loss, some simulations illustrate the proposed control scheme and the choice of a communication protocol.

Fig. 6 . 1

 61 Fig. 6.1 Overview of the system and control setup.

  Fig. 6.2 Power spectral density of the delays induced by a three-hops network.

Fig. 6 . 3

 63 Fig. 6.3 Empirical transfer function estimate (ETFE) and auto regressive moving average (ARMA) model of the network delays.

Fig. 6 . 4

 64 Fig. 6.4 Delay traces and delay estimates computed with the combined Kalman-filter and CUMSUM detector.

Fig. 6 . 5 Fig. 6 . 6

 6566 Fig. 6.5 Average end-to-end network-induced delay and reliability of Breath protocol.

Fig. 6 . 7

 67 Fig. 6.7 Block diagram of the closed loop system with plant, network, measurements and controller.

Fig. 6 . 8

 68 Fig. 6.8 Delays, system output and control signal without ARQ.

The time dependency of the variables δ(t), τ 1 (t) and τ 2 (t) are omitted in the notations from this point.

This work was supported by the EU project FeedNetBack, the Swedish Research Council, the Swedish Strategic Research Foundation, and the Swedish Governmental Agency for Innovation System. 1