Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

Lukas Einkemmer
email: lukas.einkemmer@uibk.ac.at

Erwan Faou
email: erwan.faou@inria.fr

HAMILTONIAN SPLITTING FOR THE VLASOV-MAXWELL EQUATIONS by

Keywords:

published or not. The documents may come

Hamiltonian splitting for the Vlasov-Maxwell equations

Introduction

The time evolution of an ensemble of charged particles in a plasma or a propagating beam is described by the Vlasov equation. This kinetic model describes the plasma response to electromagnetic fields. The unknown of the Vlasov equation is the distribution function f (t, x, v) of the considered plasma system, which depends on the time t, the space x, and the velocity v. The coupling of the Vlasov equation with Maxwell's equations takes into account the self-consistent fields, i.e. the electric E and magnetic field B generated by the particles in the plasma.

The so-obtained Vlasov-Maxwell system is nonlinear and thus analytical solutions are not available in general. Therefore, numerical simulations have to be conducted in order to study realistic physical phenomena (such as magnetic field generation).

Historically, Particle In Cell (PIC) methods have been widely used to numerically solve high dimensional Vlasov-Maxwell problems. The main advantage of this class of methods is their low computational cost [START_REF] Birdsall | Plasma physics via computer simulation[END_REF][START_REF] Bowers | Ultrahigh performance three-dimensional electromagnetic relativistic kinetic plasma simulation[END_REF]. In these methods, the trajectories of macro-particles are advanced by using the characteristics curves of the Vlasov equation whereas the electromagnetic fields are computed by gathering the charge and current densities of the simulated particles on a grid in physical space. Note that this grid is only up to three dimensional (as the electric and magnetic fields do not depend on velocity) whereas the phase space can consist of up to six dimensions.

More recently, deterministic methods have been developed [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF][START_REF] Filbet | Conservative numerical scheme for the Vlasov equation[END_REF][START_REF] Crouseilles | Numerical approximation of collisional plasmas by high order methods[END_REF][START_REF] Duclous | High order resolution of the Maxwell-Fokker-Planck-Landau model intend for ICF applications[END_REF]. In this approach, a phase space grid is used so that finite volumes, finite elements, or finite differences can be employed to approximate the differential operators. These methods usually require more memory when high dimensional problems are considered. However, they do not suffer from numerical heating and statistical noise. In most of the literature, a time splitting approach is employed; in such a scheme the transport in the spatial variable is split from the transport in the velocity variable ([23, 24, 20, 12, 4]). The electromagnetic fields are then advanced by approximating Maxwell's equations using a suitable time integrator. Finally, let us emphasize that it is well known that for both approaches (PIC and deterministic), the Poisson equation needs to be satisfied to machine precision (even if Poisson's equation is not solved directly in the numerical method). This is the charge conservation problem ( [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF][START_REF] Chen | An energy and charge conserving implicit electrostatic Particle In Cell algorithm[END_REF][START_REF] Langdon | On enforcing Gauss law in electromagnetic particle-in-cell codes[END_REF][START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF][START_REF] Barthelmé | Le problème de conservation de la charge dans le couplage des équations de Vlasov et de Maxwell[END_REF][START_REF] Crouseilles | Charge preserving scheme for the numerical solution of the Vlasov-Ampère equations[END_REF]).

In this work, we propose new time splitting schemes for the numerical solution of the Vlasov-Maxwell system using the deterministic approach outlined above. These splitting methods enjoy the following properties:

they are symplectic, in the sense that they preserve the mechanical Hamiltonian structure of the equation;

they can be generalized to arbitrary high order in time by composition (first order for the Lie splitting, second order for Strang splitting, fourth order for the triple Jump method, . . . ); -Poisson's equation is exactly satisfied at the semi-discrete level as well as at the fully discrete level for standard space discretization methods. they can be easily combined with arbitrary high order schemes in phase space (finite volume, semi-Lagrangian, spectral, . . . ).

These splitting methods are based on a decomposition of the Hamiltonian of the Vlasov-Maxwell equations. Splitting methods are widely used in the context of systems of ordinary differential equations governed by Hamilton's equations. In the case of velocity independent potentials for example, the application of the Strang splitting algorithm leads to the well known Störmer-Verlet method. Recently, these methods have been extensively studied in the literature (see, for example, [START_REF] Hairer | Geometrical Numerical Integration[END_REF] for an overview): as they are symplectic, they exhibit excellent long time behavior and energy preservation properties.

As far as partial differential equations are concerned, splitting methods are also well developed in the context of linear and semi-linear Hamiltonian equation such as the Schrödinger equation (see for instance [START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF][START_REF] Faou | Geometric numerical integration and Schrödinger equations[END_REF]). However, in the case of strongly nonlinear equation, where the highest order derivative depends on the solution itself, the formulation as a Hamiltonian system is significantly less common. Nevertheless, it is well known that the Vlasov-Maxwell equations can be considered as an infinite dimensional Hamiltonian system (see [START_REF] Morrison | The Maxwell-Vlasov equations as a continuous Hamiltonian system[END_REF][START_REF] Marsden | The Hamiltonian structure of the Maxwell-Vlasov equations[END_REF]) where the Hamiltonian structure is non canonical, and depends on the solution itself (Poisson structure). This is of interest, as contrary to the Vlasov-Poisson equation for which the natural splitting between spatial advection and velocity advection is also a Hamiltonian splitting (see [START_REF] Crouseilles | High order Runge-Kutta-Nyström splitting methods for the Vlasov-Poisson equation[END_REF]), it is not clear to decide from the structure of the Vlasov-Maxwell which choice of splitting scheme yields a priori good conservation properties. Furthermore, as we will discuss in more detail in section 5, the splitting methods for the Vlasov-Maxwell equations constructed in the literature have exclusively focused on two-term splittings (motivated by the two transport terms in space and in velocity) and therefore have used rather ad-hoc procedures in order to obtain methods of second order. This also impedes the construction of high order methods by composition (which is relatively straightforward for the Vlasov-Poisson equation, see [START_REF] Crouseilles | High order Runge-Kutta-Nyström splitting methods for the Vlasov-Poisson equation[END_REF]) as the resulting scheme are no longer symmetric.

In this paper, we outline the construction and numerical validation of these new splitting methods and apply them to the reduced 1+1/2 dimensional Vlasov-Maxwell system. Several numerical tests have been conducted and comparisons with two numerical schemes of the literature ( [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF][START_REF] Mangeney | A numerical scheme for the integration of the Vlasov-Maxwell system of equation[END_REF]) are performed for Landau damping, a Weibel type instability, and a magnetically induced two-stream instability.

The rest of the paper is organized as follows: first, the Vlasov-Maxwell model and its Hamiltonian structure are recalled. Then, the Hamiltonian splitting is presented and applied to the reduced 1+1/2 dimensional Vlasov-Maxwell system. Finally, numerical experiments are conducted with the new methods and comparisons with classical methods from the literature are performed.

Hamiltonian structure of the Vlasov-Maxwell system

We consider the Vlasov-Maxwell system that is satisfied by the electron distribution function f = f (t, x, v) and the electromagnetic fields (E, B) = (E(t, x), B(t, x)) ∈ R 3 × R 3 . Here, the spatial variable is denoted by x ∈ X 3 (X 3 being a three dimensional torus), the velocity variable is denoted by v ∈ R 3 , and the time is denoted by t ≥ 0. Using normalized units, the Vlasov-Maxwell system can be written as (2.1)

∂ t f + v • ∇ x f + (E + v × B) • ∇ v f = 0, ∂ t E = ∇ x × B - R 3 vf (t, x, v) dv, ∂ t B = -∇ x × E.
The two constraints on the electromagnetic fields (E, B) are given by (2.2)

∇ x • E = ρ(t, x) := R 3 f (t, x, v)dv -1, ∇ x • B = 0.
Note that if these constraints are satisfied at the initial time, they are satisfied for all times t > 0. Moreover, the total mass is preserved; that is,

X 3 R 3 f (t, x, v) dx dv = X 3 R 3 f (0, x, v) dx dv = 1
holds true for all time t > 0. Note, however, that this is not always the case when numerical approximation are considered. For a unique solution the Vlasov-Maxwel system has to be supplemented with initial conditions for the distribution function as well as the field variables; that is, we have to specify

f (t = 0, x, v) = f 0 (x, v), E(t = 0, x) = E 0 (x), B(t = 0, x) = B 0 (x).
The Hamiltonian associated with the Vlasov-Maxwell system is given by (see [START_REF] Morrison | The Maxwell-Vlasov equations as a continuous Hamiltonian system[END_REF] and [START_REF] Marsden | The Hamiltonian structure of the Maxwell-Vlasov equations[END_REF])

H = 1 2 X 3 |E| 2 dx + 1 2 X 3 |B| 2 dx + 1 2 X 3 ×R 3 |v| 2 f dv dx = H E + H B + H f . (2.3)
The three terms corresponding to electric energy, magnetic energy, and kinetic energy, respectively. For a given functional K(f, E, B), we denote by δK/δf , δK/δE and δK/δB the Fréchet derivatives of K with respect to f , E and B respectively (1) . The Poisson bracket of two functionals K(f, E, B) and G(f, E, B), is then defined as

[K, G] = X 3 ×R 3 f { δK δf , δG δf } dx dv (2.4) + X 3 δK δE • (∇ x × δG δB ) - δG δE • (∇ x × δK δB ) dx + X 3 ×R 3 δK δE • ∇ v f δG δf - δG δE • ∇ v f δK δf dx dv + X 3 ×R 3 δK δB • (∇ v f × v) δG δf - δG δB • (∇ v f × v) δK δf
∂ t K = [K, H] = [K, H E ] + [K, H B ] + [K, H f ],
for any functional K(t) = K(f (t), E(t), B(t)) evaluated along the solution of (2.5) (note that here, we assume that the solution (f, E, B) and the functional K are smooth enough to ensure the validity of the equations).

As is outlined in the next section, this formulation provides the basis for the splitting methods proposed in this paper.

Hamiltonian splitting

In this section, we propose new splitting methods to compute the solution of the Vlasov-Maxwell system. These splitting are based on exact computations of the three parts of the Hamiltonian H E , H B and H f given by (2.3) respectively.

Let us start by detailing the equations associated with the different parts of the Hamiltonian. Note that these evolution equations can be derived by plugging in a representation for H E , H B , and H f into the right hand side of equation (2.5).

∂ t f + E(x) • ∇ v f = 0, ∂ t E = 0, ∂ t B = -∇ x × E.
For a given initial data (f 0 , E 0 , B 0 ) at time t = 0, the solution of this system at time t is given explicitly by

(3.2) f (t, x, v) = f 0 (x, v -tE 0 (x)), E(t, x) = E 0 (x), B(t, x) = B 0 (x) -t ∇ x × E 0 (x).
Formally we can write this solution as

F (t) := (f, E, B) T (t) = exp(H E t)(f 0 , E 0 , B 0 ) T . Moreover, if E 0 and f 0 satisfy the relation ∇ x •E 0 = R 3 f 0 (x, v) dv -1, then it holds that ∇ x •E(t, x) = R 3 f (t, x, v) dv-1;
this can be easily seen by the fact that E(t, x) is constant in time, and that the transformation v → v -tE 0 (x) preserves the volume.

Let us further remark that if ∇ x • B 0 = 0, then this property holds true for any later time as well. This is easily shown by considering the divergence of the last equation

∇ x • B(t, x) = ∇ x • B 0 (x) -t ∇ x • (∇ x × E 0 (x)) = ∇ x • B 0 (x) = 0.
The previous relations can be easily carried over to the fully discrete case by using for example a spectral discretization in x and an interpolation procedure in v to compute f (t, x, v) from f 0 (x, v). The volume preservation in this case can be easily ensured as the advection is just a translation with constant coefficients.

Equations for H

∂ t f + (v × B(x)) • ∇ v f = 0, ∂ t E = ∇ x × B, ∂ t B = 0.
Note that in the right hand side of the first equation the following cross product term appears

v × B =   v 2 B 3 -v 3 B 2 v 3 B 1 -v 1 B 3 v 1 B 2 -v 2 B 1   =   0 B 3 -B 2 -B 3 0 B 1 B 2 -B 1 0   v =: J B v.
Two important properties of this term have to be remarked. First, the k-th component of (v × B) does not depend on v k (for k = 1, 2, 3). Second, the 3 × 3 matrix J B is constant in time (which follows immediately from the third evolution equation).

For a given initial data (f 0 , E 0 , B 0 ) at time t = 0, the solution of this system at time t is therefore given explicitly by

(3.4) f (t, x, v) = f 0 (x, exp(-J B t)v), E(t, x) = E 0 (x) + t ∇ x × B 0 (x), B(t, x) = B 0 (x).
As before, we will denote this formally by

F (t) := (f, E, B) T (t) = exp(H B t)(f 0 , E 0 , B 0 ) T . Let us remark that if ∇ x • E 0 = R 3 f 0 (x, v) dv -1, then considering the divergence of the second equation leads to ∇ x • E(t, x) = ∇ x • E 0 (x).
Now as J B is skew-symmetric, the transformation v → exp(-J B t)v preserves the volume, and we thus follow that

(3.5) R 3 f 0 (x, v) dv = R 3 f 0 (x, exp(-J B t)v) dv
which ensures that the Poisson equation is satisfied after one integration step.

Note that in contrast with the previous case, the relation (3.5) is in general difficult to ensure by standard 3D (or 2D) interpolation at the fully discrete level. To remedy this difficulty, one possibility is to further approximate f (t, x, v) by directional splitting. Indeed, due the structure of the cross product, i.e. since

v × B =   v 2 B 3 -v 3 B 2 v 3 B 1 -v 1 B 3 v 1 B 2 -v 2 B 1   ,
the evolution equation for f (t, x, v) can be written as

∂ t f + (v 2 B 3 -v 3 B 2 )∂ v 1 f + (v 3 B 1 -v 1 B 3 )∂ v 2 f + (v 1 B 2 -v 2 B 1 )∂ v 3 f = 0,
which can be further split into the following three equations

∂ t f = -(v 2 B 3 -v 3 B 2 )∂ v 1 f ∂ t f = -(v 3 B 1 -v 1 B 3 )∂ v 2 f ∂ t f = -(v 1 B 2 -v 2 B 1 )∂ v 3 f.
Each of these equations can be solved explicitely by a linear advection with constant coefficients. For example the speed of the advection in the v 1 direction does only depend on the perpendicular velocities v 2 and v 3 . As proposed in [START_REF] Mangeney | A numerical scheme for the integration of the Vlasov-Maxwell system of equation[END_REF], a Strang splitting (or high order splitting) can be performed here.

3.3.

Equations for H f . -Finally, the equation associated with the Hamiltonian H f is given by (3.6)

∂ t f + v • ∇ x f = 0, ∂ t E = - R 3 vf (t, x, v) dv, ∂ t B = 0.
For a given initial data (f 0 , E 0 , B 0 ) at time t = 0, the solution at time t is given explicitly by

(3.7) f (t, x, v) = f 0 (x -tv, v), B(t, x) = B 0 (x), E(t, x) = E 0 (x) - t 0 R 3 vf 0 (x -sv, v) dv ds.
We will denote this solution formally by

F (t) := (f, E, B) T (t) = exp(H f t)(f 0 , E 0 , B 0 ) T .
Let us remark that the Poisson equation is propagated with time. Indeed, considering the divergence of the last equation leads to

∇ x • E(t, x) = ∇ x • E 0 (x) - t 0 R 3 v • ∇ x [f 0 (x -sv, v)] dv ds = R 3 f 0 (x, v) dv -1 + t 0 R 3 ∂ s [f 0 (x -sv, v)] dv ds = R 3 f 0 (x, v) dv -1 + R 3 [f 0 (x -tv, v) -f 0 (x, v)] dv = R 3 f 0 (x -tv, v) dv -1 = R 3 f (t, x, v) dv -1,
which implies that the splitting method proposed here satisfies the charge conservation property. The preservation of the charge conservation at the fully discrete level relies on a good approximation of the time integral in the previous formula, in combination with the interpolation procedure in the v variable. We will give an example in the next section where E(t, x) is easily computed using Fourier spectral approximations in the x variable.

3.4. Hamiltonian splitting. -The solution F (t) = (f, E, B) T (t) of the Vlasov-Maxwell system up to first order can be written as (this is the well known Lie splitting)

F (t) = exp(H E t) exp(H B t) exp(H f t)F (0).
The Strang splitting scheme (which is a second order accurate splitting method) is given by

F (t) = exp(H E t/2) exp(H B t/2) exp(H f t) exp(H B t/2) exp(H E t/2)F (0).
Higher order methods can be constructed using composition (see [START_REF] Hairer | Geometrical Numerical Integration[END_REF]III.3]). For example the well known triple jump scheme (see [16, II, Example 4.2]) based on the Strang splitting yields a method of order four, as the original Strang splitting is symmetric.

Application to the reduced 1 + 1/2 dimensional model

In this section, we apply the splitting method proposed in this paper to a fully discretized 1+1/2 model of the Vlasov-Maxwell system. This model has been studied in [START_REF] Palodhi | Nonlinear kinetic development of the Weibel instability and the generation of electrostatic coherent structures[END_REF][START_REF] Cheng | Discontinuous Galerkin Methods for Vlasov-Maxwell Equations[END_REF] and retains most of the properties of the full Vlasov-Maxwell system. All the numerical simulations in this paper will be conducted in the framework of this reduced model.

4.1. Reduced 1+1/2 model. -We consider the phase space (x 1 , v 1 , v 2 ) ∈ X × R 2 ,
where X is a one-dimensional torus, and the unknown functions

f (t, x 1 , v 1 , v 2 ), B(t, x 1 ) and E(t, x 1 ) = (E 1 , E 2 )(t, x 1
) which are determined by solving the following system of evolution equations (4.1)

∂ t f + v 1 ∂ x 1 f + E • ∇ v f + BJ v • ∇ v f = 0, ∂ t B = -∂ x 1 E 2 ; ∂ t E 2 = -∂ x 1 B - R 2 v 2 f (t, x 1 , v) dv, ∂ t E 1 = - R 2 v 1 f (t, x 1 , v) dv,
where v = (v 1 , v 2 ) and J denotes the symplectic matrix

J = 0 1 -1 0 .
This reduced system correspond to choosing an initial value of the form

E(x 1 , x 2 , x 3 ) =   E 1 (x 1 ) E 2 (x 1 ) 0   and B(x 1 , x 2 , x 3 ) =   0 0 B(x 1 )
  and a f depending on x 1 and (v 1 , v 2 ) only. Then it can be easily checked that this structure is preserved by the flow of (2.1) and that the equations reduce to (4.1). We refer the reader to [START_REF] Cheng | Discontinuous Galerkin Methods for Vlasov-Maxwell Equations[END_REF] for more details.

4.2. Phase space integration. -Consistent with most of the literature, we consider the periodic case in the x 1 direction. Thus, one natural method in this direction is the Fourier spectral approximation. It simplifies the computation of the evolution corresponding to H f , and the exact charge conservation can be easily recovered. In the evolutions corresponding to H E and H B , we compute an approximation to the transport of f by using a third order finite volume numerical scheme (see [START_REF] Filbet | Conservative numerical scheme for the Vlasov equation[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF]) whereas the equations for the electric field E and the magnetic field B are solved in Fourier space.

Let us consider the evolution corresponding to H f in more detail. By employing the Fourier transform in space (that is in the x 1 variable) and denoting the corresponding functions by fk , k ∈ Z, (that is, fk denotes the Fourier transform of f (x) and does still depend on v 1 and v 2 ) gives

∂ t f + v 1 ik f = 0, ∂ t Ê1 = - R 2 v 1 f dv.
The first equation can be solved exactly in time, given the initial condition f (t n ):

∀ k ∈ Z, f (t) = f (t n ) exp(-iv 1 k(t -t n )),
which inserted into the second equation (which we have integrated in time from t n to t n+1 ) yields the approximation

Ên+1 1 = Ên 1 - t n+1 t n R 2 v 1 f (t n ) exp(-iv 1 k(t -t n )) dv dt = Ên 1 - R 2 v 1 f (t n ) t n+1 t n exp(-iv 1 k(t -t n )) dt dv = Ên 1 - R 2 v 1 f (t n ) -1 ikv 1 (exp(-ikv 1 ∆t) -1) dv.
For the long time properties of the numerical scheme it is vital that the total charge is conserved exactly and that Poissons equation is satisfied. The charge is not part of the formulation of our numerical algorithm, but it can be shown that our scheme satisfied the charge conservation property. Indeed, assuming that Poisson's equation is satisfied initially, we can prove that it is satisfied for all times.

Assuming that Poisson equation is satisfied at time t n , that is ik

Ê1 n = ρn ,
where ρn is an approximation of the charge density ρ(t, x) in Fourier, see (2.2), we get from the last equation

Ên+1 1 = Ên 1 - R 2 v 1 f (t n ) -1 ikv 1 (exp(-ikv 1 ∆t) -1) dv = ρn ik + 1 ik R 2 f (t n )(exp(-ikv 1 ∆t) -1) dv = 1 ik ρn + R 2 f (t n ) exp(-ikv 1 ∆t) dv -ρn -1 = 1 ik ρn+1 ,
which shows that Poisson's equation is also true at time t n+1 , assuming a sufficient accuracy in the discretization of the integral in v. Note that k = 0 does not pose any additional difficulties as the average of the electric field is imposed to be zero by total mass conservation. This is a necessary condition to impose as periodic boundary conditions do not give a unique solution for the electric field.

Once the three steps are solved, the global time integration is given in subsection 3.4.

Numerical experiments

In this section, the Hamiltonian splitting is numerically studied by solving the 1+1/2 dimensional reduced Vlasov-Maxwell system introduced in the previous section. We also compared our numerical scheme to two methods of the literature: VALIS (from [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF], recalled in subsection 6.2) and Mangeney (from [START_REF] Mangeney | A numerical scheme for the integration of the Vlasov-Maxwell system of equation[END_REF], recalled in 6.1).

Landau damping.

-Landau damping is a classic test case for the Vlasov-Poisson equations. In this context a number of numerical schemes have been proposed (see e.g. [START_REF] Filbet | Conservative numerical scheme for the Vlasov equation[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF][START_REF] Crouseilles | High order Runge-Kutta-Nyström splitting methods for the Vlasov-Poisson equation[END_REF]). However, even though the equation we try to solve does not explicitly determine the electric field from the charge density (as is usually done for methods constructed in the context of the Vlasov-Poisson equations) it clearly reduces to the Vlasov-Poisson equation if no magnetic effects are present in a given configuration. Compared to the simulation discussed in sections 5.3 and 5.2 this problem does not rely on magnetic effects. Thus, it serves as a first, but rather limited, consistency check for the numerical schemes considered in this paper. However, the charge conservation property is crucial in this test; if it is not satisfied exactly, the solution of the Vlasov-Poisson system is not recovered.

For this test, the initial value for the distribution function f is given by

f 0 (x 1 , v 1 , v 2 ) = 1 2π e -1 2 |v| 2 (1 + α cos kx 1 ) ,
where x 1 ∈ X = [0, 2π/k], v ∈ R 2 and we have chosen α = 0.5 and k = 0.4, a configuration usually called strong (or nonlinear) Landau damping. The electric field, at the initial time, is determined by Poisson's equation, i.e.

∂ x 1 E 1 (t = 0, x 1 ) = R 2 f 0 (x 1 , v) dv -1,
whereas a vanishing magnetic and perpendicular electric field (i.e. E 2 ) is prescribed in order to completely specify the initial value.

The quantity of interest is the potential energy stored in the electric field

E pot (t) = 1 2 X |E(t, x 1 )| 2 dx 1 ,
and the total energy (sum of the potential and kinetic energy)

E tot (t) = E pot (t) + E kin (t),
where

E kin (t) = 1 2 X R 2 |v| 2 f (t, x 1 , v) dv dx 1 .
Note that the total energy is conserved with time.

In Figure 1, the time history of the electric energy is shown for a Strang Hamiltonian splitting, using 32 points in the spatial direction, 64 points in each velocity direction and a time step ∆t = 0.05. Note that compared to the linear case (i.e. α = 0.01, for example) the electric energy does only decay exponentially for a brief period of time before presenting an oscillating behavior for large times. The results is in a good agreement with results available in the literature (see [START_REF] Manfredi | Long time behavior of nonlinear Landau damping[END_REF]).

In Figure 2 the time history of the total energy error is plotted for a Lie and Strang Hamiltonian splitting (the scheme proposed in this paper) as well as the VALIS scheme described in [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF] (for more details see appendix 6). On the upper figure, we consider a phase space mesh of 32 × 64 × 64 and the time step is ∆t = 0.1. We can first observe that the three methods behave well as the error in energy is less than 1%. However, the Lie splitting leads to oscillations which is not the case when a second order charge preserving method (Strang splitting or VALIS) is used. When the numerical parameters are refined (lower figure), these two methods show virtually indistinguishable results (the error is about 0.1%) and we do not seem to experience a significant growth over the time period considered in this simulation. Let us remark that the charge conservation is very important with respect to energy conservation. If charge conservation is not exact, the accuracy of the algorithm is significantly reduced. This is more pronounced if we consider an increased number of grid points in space. The scheme proposed in [START_REF] Mangeney | A numerical scheme for the integration of the Vlasov-Maxwell system of equation[END_REF] for that reason only performs satisfactory if a significantly smaller time step is chosen (as compared to Hamiltonian splittings and VALIS). This phenomenon can be observed in all the simulations conducted in this section and we will therefore omit it from the plots. We will, however, show a comparison for a reduced step size in the fully magnetized case (section 5.3).

Weibel instability.

-In this section, we will discuss an incarnation of the Weibel instability. This instability includes a genuine magnetic effect. Following [START_REF] Palodhi | Nonlinear kinetic development of the Weibel instability and the generation of electrostatic coherent structures[END_REF], we impose the following initial particle density

f 0 (x 1 , v 1 , v 2 ) = 1 πv th √ T r e (v 2 1 +v 2 2 /Tr)/v th (1 + α cos(kx 1 )),
where

x 1 ∈ X = [0, 2π/k], v ∈ R 2
and we have chosen v th = 0.02, T r = 12, k = 1.25, and α = 10 -4 . The electric field E 1 is initialized by solving Poisson's equation and E 2 is set to zero at t = 0. The magnetic field is prescribe as

B 0 (x 1 ) = β cos(kx 1 ),
where we have chosen β = 10 -4 .

From Figure 3 we see that the electric field in the y-direction does show an exponential growth in magnitude. The (analytically derived) growth rate is proportional to k/ √ 1 + k 2 and, for the parameters, chosen here is approximately equal to 0.031. We observe that, similarly to [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF], the agreement to the theory is as well as one would expect given that the analytical growth rate is derived using the linearized Vlasov-Maxwell equations. Let us also note that even though the qualitatively features of the solution are well resolved on a grid with 32x64x64 points, quantitatively the result is wrong by at least an order of magnitude. As in the previous example, we observe that the scheme proposed by Mangeney does not produce accurate results with respect to energy conservation, except if the chosen step size is significantly smaller than what is necessary in case of the VALIS or the Hamiltonian splitting scheme (several orders of magnitude). However, compared to the simulation presented in the previous section, for sufficiently long integration times, we observe a better energy conservation for the Hamiltonian Strang splitting scheme as compared to the VALIS implementation (which is a scheme of second order as well).

Finally, let us investigate the behavior of the VALIS scheme more closely. As can be seen from Figure 5, if a time step close to the CFL condition is chosen we can observe (unphysical) oscillations in the solution. Let us duly note that such oscillations are not present in our Hamiltonian splitting approach. In fact, we can verify that these oscillations disappear for the VALIS scheme when the time step size is sufficiently small. This behavior implies that for the problem under considered in this section, the VALIS scheme has a significantly higher computational cost compared to the Hamiltonian Strang splitting scheme proposed in this paper. average velocity where x 1 ∈ X = [0, 2π], v ∈ R 2 and β = 2 • 10 -3 . The resulting instability is usually driven by a perturbation in the particle density. However, as we are more interested in magnetic effects we will drive the instability by a perturbation in the magnetic field only B 0 (x 1 ) = α sin(x 1 ), where we have chosen α = 10 -3 .

f 0 (x 1 , v 1 , v 2 ) = 1 2πβ e -v 2 2 /β (e -(v 1 -0.2) 2 /β + e -(v 1 +0.2) 2 /β ), 1e
The result of the simulation is shown in Figure 6. Even though we drive the instability by a magnetic perturbation only, we observe, after an initial oscillatory regime, an exponential increase and a saturation in the electric field.

In Figure 7 the error in energy for different schemes is shown. Note that for the scheme proposed by Mangeney, even for a relatively coarse space discretization, the step size has to be reduced by at least an order of magnitude in order to get comparable results to the VALIS and Hamiltonian Strang splitting schemes. In addition, we observe that, especially for a sufficiently fine space discretization, the error in energy for the VALIS and our Hamiltonian Strang splitting scheme is almost indistinguishable. scheme. The order of the triple jump scheme is of particular interest as the results provided in this section demonstrate that we can in fact use composition to obtain a method that is of order four. The numerical results for the Weibel instability are shown in Table 1 for a space discretization using 64 grid points and in Table 2 for a space discretization using 256 grid points. In both cases we compute the l 1 error of the field variables (i.e. the two components of the electric field and the single component of the magnetic field) by using a reference solution for which a sufficiently small time step has been chosen.

Appendix: Other algorithms

We present in this section some algorithms found in the literature [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF]. The algorithms only differ from the computation of the Maxwell equations. The computation of the distribution function through the Vlasov equation is always performed using a Strang splitting 1. Order determination for the Hamilton Lie, Strang, and triple jump splitting scheme in case of the Weibel instability introduced in section 5.2. In these simulations 64 grid points are used in each dimension and the system is integrated up to final time t = 1.

-f = f n exp(-ikv 1 ∆t/2) solution of ∂ t f + v 1 ∂ x 1 f = 0 on ∆t/2, -compute f by solving ∂ t f +(E +BJ v)•∇ v f = 0 on ∆t, using a splitting procedure, -f n+1 = f exp(-ikv 1 ∆t/2) solution of ∂ t f + v 1 ∂ x 1 f = 0 on ∆t/2.
compute Note that this algorithm conserves the charge.

B n+1/2 = B n -∆t/2 ∂ x 1 E n 2 -compute E n+1/2 2 = E n 2 -∆t/2 ∂ x 1 (B n-1/2 + B n+1/2 )/2 -∆t/2J n 2 -compute E n+1/2 1 = E n 1 -∆t/2 J n 1 -compute f n+1 with Strang splitting, the fields (E n+1/2 , B n+1/2 ) being known at time t n+1/2 , -compute the current J n+1 1,2 = R 2 v 1,2 f n+1 (x 1 , v) dv -compute E n+1 2 = E n 2 -∆t∂ x 1 B n+1/2 -∆t/2(J n 2 + J n+1 2 ) -compute E n+1 1 = E n 1 -∆t/2(J n 1 + J n+1 1 ) 6
compute the current In [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF], the current is extracted from the two x 1 -advection steps. Using Fourier transforms, it is written f = f n exp(-ikv 1 ∆t/2) = f n -ik ∆t 2 f n [exp(-ikv 1 ∆t/2) -1] /(ik∆t/2), which becomes, after integration in v ∈ R 2 , ρ = ρn -ik ∆t 2 Ĵ 1 with

J n 1,2 (x 1 ) = R 2 v 1,2 f n (x 1 , v) dv -compute B n+1/2 = B n -∆t/2 ∂ x 1 E n 2 -compute E n+1/2 2 = E n 2 -∆t/2 ∂ x 1 (B n-1/2 + B n+1/2 )/2 -∆t/2J n 2 -compute E n+1/2 1 = E n 1 -∆t/2 J n 1 -compute f n+1 and the current • from f , compute and store Ĵ 1,2 = R 2 v 1,2 f dv = R 2 f n (exp(-ikv 1,2 ∆t/2) -1)/(ik∆t/
Ĵ 1 = R 2
f n [exp(-ikv 1 ∆t/2) -1] /(ik∆t/2) dv.

Similarly the second part of the current is denoted by

Ĵ 1 = R 2
f [exp(-ikv 1 ∆t/2) -1] /(ik∆t/2) dv.

Gathering these two steps, we can compute the density at time t n+1 by integrating in v each step of the splitting for f which is exactly the Poisson equation at time t n+1 . Note that this can be generalized for finite volume method following the ideas of [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF].

Conclusion

In this work, a new time splitting is proposed for the numerical solution of the Vlasov-Maxwell system. This splitting is based on a decomposition of the Hamiltonian of the Vlasov-Maxwell system. Each step can be solved exactly in time so that the error is only due to the splitting procedure (first order for Lie, second order for Strang, fourth order for the triple jump scheme). In addition to the fact that methods of arbitrary high order can be constructed by composition, we prove that when one uses a spectral method for the spatial discretization, the charge is exactly preserved. Several numerical tests show the very good behavior of the splitting compared to methods from the literature, especially regarding the conservation of total energy.

Several further perspectives of this work can be envisaged. First, a generalization to finite volumes methods for the spatial discretization have to be performed, taking care of the charge conservation property. We are confident this can be achieved by combining ideas of [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D Vlasov-Maxwell system[END_REF] and [START_REF] Crouseilles | Charge preserving scheme for the numerical solution of the Vlasov-Ampère equations[END_REF]. Second, a more detailed study of high order splittings (such as the triple Jump scheme) is necessary, since the use of small steps due to the CFL condition imposed by the Maxwell equations is a limiting factor. A substepping methodology would alleviate this problem. Finally, considering more realistic problems in 4 or 5 dimensions will require the parallelization of the code in the context of distributed memory architectures (i.e. using MPI). Due to the fact that the VALIS method requires more temporary variables than the splitting method proposed in this paper, a more efficient MPI parallelization of our Hamiltonian splitting seems to be possible and this will constitute further research.

3. 1 .

 1 Equations for H E . -The equations associated with the Hamiltonian H E are given by(3.1)

  B . -The equations associated with the Hamiltonian H B are given by(3.3)

Figure 1 .

 1 Figure 1. Time evolution of the electric energy E pot (semi-log scale) for the Landau damping test using the Hamiltonian Strang splitting with 32x64x64 grid points and ∆t = 0.05. Top figure: nonlinear Landau damping (α = 0.5); bottom figure: linear Landau damping (α = 0.01). The analytic decay rate γ ≈ 0.066 is shown for comparison in case of the linear Landau damping.

Figure 2 .

 2 Figure 2. Time evolution of the error in energy (semi-log scale) for the Landau damping test |E tot (t) -E tot (0)|. Top figure: 32x64x64 grid points and ∆t = 0.1; bottom figure: 256x256x256 grid points and ∆t = 0.0125.

E 1 Figure 3 .

 13 Figure 3. Time evolution of the k = 1.25 Fourier mode of the electric and magnetic fields (semi-log scale) for Weibel instability using the Hamiltonian Strang splitting method. Top figure: 32x64x64 grid points and ∆t = 0.2; bottom figure: 256x256x256 grid points and ∆t = 0.0125.

Figure 4 .

 4 Figure 4. Weibel test: time evolution of the total energy error |E(t) -E(0)| (semi-log scale). Top figure: 32x64x64 grid points and ∆t = 0.2; bottom figure: 256x256x256 grid points and ∆t = 0.0125.

Figure 5 .

 5 Figure 5. Weibel test: time evolution of the first mode of E 2 (semilog scale). Comparison between VALIS, and Strang splitting scheme using 32x64x64 grid points.

5. 4 .Figure 6 .

 46 Figure 6. Time evolution of the electric, potential and total energies (semi-log scale) for the two-stream instability test using the Strang splitting scheme with 32x64x64 grid points and a time step of ∆t = 0.1.

Figure 7 .

 7 Figure 7. Time evolution of the error in energy |E tot (t) -E tot (0)| (semi-log scale) of the two-stream instability test. Top figure: 32x64x64 grid points and ∆t = 0.1 for VALIS, Strang and Lie, and ∆t = 0.01 for Mangeney; bottom figure: 256x256x256 grid points and ∆t = 0.0125.

ρn+1 = ρ -ik ∆t 2 Ĵ 1 = ρ -ik ∆t 2 Ĵ 1 ( 2 Ĵ 1 -ik ∆t 2 Ĵ 1 = ρn -ik ∆t 2 Ĵ 1 + Ĵ 1 = ρn -ik∆t Ĵn+1/ 2 1 2 1to 1 - 1 ) = ik Ên+1 1 ,

 2121212121122111 since the advection in v is conservative) = ρn -ik ∆t Assuming the Poisson equation satisfied at time t n , we get ik Ên 1 = ρn and since the current is used Ĵn+1/Ên

Table

  6.1. Predictor-Corrector. -At time t n , the electric field E n and the magnetic field B n-1/2 are known, together with f n .compute the currentJ n 1,2 (x 1 ) = R 2 v 1,2 f n (x 1 , v) dv

		Lie splitting scheme
	stepsize l 1 error	order
	0.2	1.881556551139132e-06
	0.1	8.769892305391039e-07 1.1013
	0.05	4.237597343971976e-07 1.04931
	0.025	2.083331679902264e-07 1.02435
	0.0125 1.033058984634212e-07 1.01197
	0.00625 5.144678908873264e-08 1.00577
		Strang splitting scheme
	stepsize l 1 error	order
	0.2	2.124324456920513e-07
	0.1	5.27651798019974e-08 2.00935
	0.05	1.315400004930548e-08 2.00408
	0.025	3.276015657686228e-09 2.00549
	0.0125 8.068242472080718e-10 2.02162
	0.00625 1.86093775964047e-10 2.11622
	Triple jump splitting scheme
	stepsize l 1 error	order
	0.2	1.321445314225445e-08
	0.1	8.095000243664432e-10 4.02894
	0.05	5.162111972126581e-11 3.971
	0.025	2.031275496535852e-12 4.6675
	0.0125 4.297134131610342e-13 2.24094

Table 2 .

 2 .2. Valis. -At time t n , the electric field E n and the magnetic field B n-1/2 are known, together with f n . Order determination for the Hamilton Lie, Strang, and triple jump splitting scheme in case of the Weibel instability introduced in section 5.2. In these simulations 256 grid points are used in each dimension and the system is integrated up to final time t = 1.

		Lie splitting scheme
	stepsize l 1 error	order
	0.2	1.823383817066824e-06
	0.1	8.465392041182712e-07 1.10697
	0.05	4.082184424389359e-07 1.05224
	0.025	2.004799316967084e-07 1.02588
		Strang splitting scheme
	stepsize l 1 error	order
	0.2	2.271250317272867e-07
	0.1	5.640193594261031e-08 2.00967
	0.05	1.407771489982866e-08 2.00233
	0.025	3.51820025143404e-09 2.0005
	Triple jump splitting scheme
	stepsize l 1 error	order
	0.05	5.922461122797243e-11
	0.025	3.371093950592525e-12 4.13491

  2) dv• compute f by solving the (conservative) advection in v 1 , v 2 , the fields being known at time t n+1/2 , • from f n+1 , compute and store Ĵ

	1,2 =	R 2 v 1,2	f n+1 dv =	R 2 f (exp(-ikv 1,2 ∆t/2)-1)/(ik∆t/2) dv
	-compute the current	Ĵn+1/2 1,2	= ( Ĵ 1,2 + Ĵ 1,2 )/2,
	-compute E n+1 2 -compute E n+1 1	= E n 2 -∆t∂ x 1 B n+1/2 -∆tJ 2 n+1/2 = E n n+1/2 1 -∆tJ 1

[START_REF] Barthelmé | Le problème de conservation de la charge dans le couplage des équations de Vlasov et de Maxwell[END_REF] Note that in the absence of domain boundaries, this definition of functional derivative is not ambiguous.

This work is supported by the Fonds zur Förderung der Wissenschaften (FWF) -project id: P25346. This work was partly supported by the ERC Starting Grant Project GEOPARDI.


Marie Duhamel

Prof Philippe Vandenkoornhuyse
email: philippe.vandenkoornhuyse@univ-rennes1.fr

Sustainable agriculture: possible trajectories from mutualistic symbiosis and plant neodomestication

Food demand will increase concomitantly with human population.

A sustainable food production ?

Feeding the world and securing access to food are both major social and scientific issues. Food crises have often occurred in the past. In recent years, the rapidly increasing demand for food (i.e., for human populations and livestock) along with biofuels has led to food price volatility [START_REF] Battisti | Historical warnings of future food insecurity with unprecedented seasonal heat[END_REF]. Furthermore, recent work suggests that food crises are exacerbated by global warming. It has been clearly indicated that agricultural productivity has declined worldwide as a consequence of the hottest summers experienced in the recent past, and according to different global warming scenarios "… the hottest seasons on record will represent the future norm in many locations ..." [START_REF] Battisti | Historical warnings of future food insecurity with unprecedented seasonal heat[END_REF]. Human population is continuously increasing, but expected to peak before the end of the century, with 10 billion people before 2100 [START_REF] Lutz | The end of world population growth[END_REF]. Contrary to common assumption, non-linearities between population expansion and environmental degradation are likely to increase disproportionately rapidly [START_REF] Harte | Human population as a dynamic factor in environmental degradation[END_REF]. Human population expansion will be coupled with an increased demand for space, water and food. These demands will therefore be accompanied by urban and cropland expansion, and more than 10 9 hectares of natural ecosystems are likely to be lost by 2050 [START_REF] Tilman | Forecasting agriculturally driven global environmental change[END_REF]. This represents collateral damage for the environment because cropland expansion can only be achieved by replacing non agricultural, mainly forested areas. According to recent studies, agricultural production will have to expand by about 100% during the 21 st century to satisfy forecasted world demands [START_REF] Cirera | Income distribution trends and future food demand[END_REF]. At the same time, agriculture is a major threat to the environment eventually leading to a decline in biodiversity and related ecosystem services, including degradation of soil and water quality [START_REF] Foley | Global consequences of land use[END_REF].

A fundamental issue for agriculture during this century is thus to confront two contradictory goals, (i) the need to produce enough food to minimize human malnutrition and support world population expansion and (ii) the need to limit collateral damage to the environment, which can in turn negatively impact agriculture. Based on recent findings about strategies in plant mutualisms and plant selection, we develop new ideas in this paper and suggest trajectories for a more sustainable agricultural development.

Intensive vs. extensive agriculture ?

The aim in intensive agriculture is to maximize productivity per unit of surface, whereas in extensive agriculture, lower productivity yields are accepted as a counterpart to less potential ecosystem damage. The main advantage of extensive agriculture is that no or few inputs are required. However this is often countered by a need for a larger soil area to obtain comparable production. It has been shown that agricultural intensification with high yield production eventually increases greenhouse gas emissions per unit surface. However, much higher carbon emissions can be expected if the same production is obtained by expanding low-yield farming (e.g. [7;8]). Similarly, the need to increase agricultural productivity to limit adverse effects on the environment has also been underlined by modeling land use/land cover changes [START_REF] Nelson | Projecting global land-use change and its effect on ecosystem service provision and biodiversity with simple models[END_REF] and by projecting possible improvements of productivity in existing agricultural areas [START_REF] Foley | Solution for a cultivated planet[END_REF]. One key element which has emerged is the necessity for agricultural intensification to preserve biodiversity and the related ecosystem services. As developed below, new ideas for maintaining high crop productivity with lower inputs have recently been put forward.

Crop selection from traits?

Since the beginning of agriculture, crops have been selected for different traits, including plant productivity. The main current approach to modern plant breeding is to maximize the fitness of individual plants. However other contrasting breeding strategies have been suggested. One of the most exciting of these new solutions would be to base plant breeding on group selection rather than on individual plant fitness [START_REF] Weiner | Evolutionary agroecology, the potential for cooperative, high density, weed-suppressing cereals[END_REF] (where group selection refers to the selection '… for attributes that increase total crop yield but reduce plants' individual fitness...' [START_REF] Weiner | Evolutionary agroecology, the potential for cooperative, high density, weed-suppressing cereals[END_REF]). This would imply a completely new approach to selection criteria. For example, selecting for cooperative shading, which would allow a passive control of weeds, seems promising to improve yield and sustainability [START_REF] Weiner | Evolutionary agroecology, the potential for cooperative, high density, weed-suppressing cereals[END_REF].

In all these breeding approaches (i.e. individual selection and group selection), however, plants are always considered as standalone entities, which is arguably a mistake. Plants are deeply dependent on mutualist microorganisms for their growth, and these can be damaged by conventional agricultural practices and current plant breeding strategies.

Arbuscular mycorrhiza and consequences of agricultural practices

The arbuscular mycorrhizal symbiosis is responsible for massive global nutrient transfer (Box 1). It is a mutualism 'that helps feed the world' [START_REF] Marx | The roots of plant-microbe collaborations[END_REF]. Arbuscular mycorrhizal fungi, because of their functions, can be considered as key microorganisms for soil productivity.

Intensive agricultural management (i.e., conventional agriculture in Europe and North America) has exerted a high selection pressure on microorganisms through profound modification of their habitats and niches, notably brought about by tillage, the high increase of mineral nutrients, and low plant diversity (i.e., crops). Tillage, ploughing and ripping, for example, represent an intense form of soil disruption. In natural habitats, AM mutualism is not subjected to perturbations of this intensity. Such disruption leads to degradation of the hyphal network, ecological functions, and AM fungal diversity [START_REF] Verbruggen | Evolutionary ecology of mycorrhizal functional diversity in agricultural systems[END_REF]. Soil nutrient availability is a strong driving influence for producing an evolved geographic structure in AM mutualism (i.e., a coevolutionary selection mosaic) [START_REF] Johnson | Resource limitation is a driver of local adaptation in mycorrhizal symbioses[END_REF]. As a result, soil fertilization in agricultural ecosystems has had a negative impact on AM fungal functions [START_REF] Johnson | Can fertilization of soil select less mutualistic mycorrhizae[END_REF] and diversity [START_REF] Egerton-Warburton | Mycorrhizal community dynamics following nitrogen fertilization, a cross-site test in five grasslands[END_REF]. Thus confounding factors, related to conventional agricultural trajectories, likely act synergistically against mycorrhizal symbiosis.

Mutualistic strategy and agriculture

From a theoretical point of view, mutualisms (i.e. cooperative interactions among different species) can exhibit instability: individuals potentially benefit from defecting from cooperation if cooperation is costly. Organisms will increase their own fitness, even if this comes to a cost of others. Kiers et al, [START_REF] Kiers | Reciprocal rewards stabilize cooperation in the mycorrhizal symbiosis[END_REF] have demonstrated the capacity of plants to sanction less-cooperative strains (i.e. 'cheaters') through a carbon embargo. The gain in fitness for the cheater is therefore reduced by this plant trait. This in itself can explain the stability of this symbiosis. A similar sanction of carbon allocation has been observed in the case of nitrogen-fixing nodules in leguminous plants to control Rhizobium cheaters [START_REF] Kiers | Human selection and the relaxation of legume defences against ineffective rhizobia[END_REF]. The most cooperative AM fungal symbionts transfer more phosphorus to the roots when they receive more carbon [START_REF] Kiers | Reciprocal rewards stabilize cooperation in the mycorrhizal symbiosis[END_REF]. Such mutualism is therefore bilaterally controlled because both partners can enforce the cooperation and any possible enslavement strategy is also limited. This fairly explains the stability of arbuscular mycorrhizal symbiosis. In addition, the main advantage for the plant to not enslave its symbionts is the access offered to numerous potential functions harbored by the reservoir of soil AM fungi into which the plant can tap depending on its nutritional requirements. For the fungi, the main advantage of not being enslaved is to be able to maintain a high level of diversity. This symbiosis is one reason for the success of plants in terrestrial ecosystems.

Less cooperative AM fungi do exist in nature. We can expect them to become more abundant as the diversity of AM fungi decreases because the symbiotic options offered to the plants are more limited. It has been shown that AM fungi cheaters can develop 'dealer' behavior by keeping phosphorus in polyphosphate chains and delivering it at an expensive cost for the host plant [START_REF] Kiers | Reciprocal rewards stabilize cooperation in the mycorrhizal symbiosis[END_REF]. The plant's capacity to sanction cheaters is a tremendously important trait to maintain, given the fact that most mineral nutrients (~70% of the phosphorus for example) are delivered to plants by AM fungi [START_REF] Smith | Mycorrhizal fungi can dominate phosphate supply to plants irrespective of growth responses[END_REF] (Box 1) in 'natural' environments.

New ideas for more sustainable agricultural practices by promoting mutualisms

Ecosystem productivity has been shown to be driven by AM symbiosis diversity [e.g. 20]. Thus, AM fungi constitute a key compartment of soil fertility.

The plant can be colonized by a variety of AM fungi (i.e., no host-specificity).

However, the recent findings suggest that plants can choose to reward and enroll some fungal colonizers in order to ensure access to particular functions related to their needs [START_REF] Kiers | Reciprocal rewards stabilize cooperation in the mycorrhizal symbiosis[END_REF]. This selective rewarding is likely to lead to the exclusion of certain colonizers and culminate in an observed 'host-plant preference' e.g. [21;22]. This leads to the idea that a plant can filter soil AM fungi depending on its requirements, the season and location. Conventional field-based agriculture makes use of very limited crop plant diversity, fungicides, soil tillage and fertilizer. The pressure exerted by agricultural practices leads to a reduction in AM fungal diversity compared to more natural ecosystems e.g. [23;24]. Breeders generally select crop cultivars from rich soils which have been under conventional agriculture for many years. In fact, the ultimate result of this selection strategy is to produce a plant that is best adapted to current agricultural practices and the related agrosystems anthropization. Agricultural soils have been enriched with fertilizers for decades and the ecological function of AM fungi as plant phosphorus providers is less important in these enriched soils. This, together with the breeding trajectory, will have relaxed the plant sanction trait in modern crops, as is the case in soybean (Glycine max (L.) Merr.) where ancient varieties are better able to control Rhizobia cheaters than modern ones [START_REF] Kiers | Human selection and the relaxation of legume defences against ineffective rhizobia[END_REF]. From an interesting meta-analysis performed from 39 publications it appears that there is '...no evidence that new crops plant genotypes lost their ability to respond to mycorrhiza due to agricultural and breeding practices...' [START_REF] Lehman | Mycorrhizal responsiveness trends in annual crop plants and their wild relativesa meta-analysis on studies from 1981 to 2010[END_REF].

Two alternative hypotheses for AM symbiosis can be put forward. First, we can hypothesize that the same trend as in the Rhizobium/legume mutualism will have already occurred for AM mutualism with a resulting loss of the sanction trait against AM fungal cheaters. As a consequence, an increase in AM fungal cheaters can be expected in agricultural soils. Because AM fungi constitute a fundamental component of soil fertility, solutions for a more ecologically intensive agriculture should focus on this compartment (Box 2.).

Plant breeders could also imagine new selection trajectories where the sanction trait is considered as a major selection target (i.e. the capacity of plants to punish bad cooperators by a carbon embargo [START_REF] Kiers | Reciprocal rewards stabilize cooperation in the mycorrhizal symbiosis[END_REF]). In this way the possibilities offered by AM functional efficiency could be restored and agricultural practices modified by reducing soil inputs and tillage (Box 2). The alternative hypothesis is that plant breeders have selected cultivars that are very efficient for mineral foraging through soil AM fungal mutualists. This apparently optimistic hypothesis is worse than that of a loss of the sanction trait in crops, because of the lack of long term sustainability. Furthermore, one important component of soil fertilizer, phosphorus, is known to rely on high quality rock phosphate, which is a finite resource. More than 85% of the global phosphate resources are dominated by only 3 countries which is far fewer than the number of countries controlling the world's oil reserves [START_REF] Elser | A broken biogeochemical cycle[END_REF]. Phosphorus (P) supply is thus of strategic importance for many countries, and "...many food producers are in danger of becoming completely dependent on this trade..." [START_REF] Elser | A broken biogeochemical cycle[END_REF]. Major agricultural regions such as India, America, and Europe are already dependent on P imports. Phosphate market prices can soar, as shown by the 700% increase in 2008 [START_REF] Elser | A broken biogeochemical cycle[END_REF], especially as phosphate mining production is predicted to attain a peak in 2030 [START_REF] Cordell | The story of phosphorus, global food security and food for thought[END_REF].

Other plant mutualisms, in addition to arbuscular mycorrhiza, should potentially have a synergistic impact on plant productivity and plant resistance against stresses (Box 2). For example, infection of barley (Hordeum vulgare) with an endophytic fungus, Piriformosa indica, increases resistance to stresses including salinity and systemic resistance of the crop to root and leaf pathogens, and a concomitant increase in yield production [START_REF] Waller | The endophytic fungus Piriformospora indica reprograms barley to salt-stress tolerance,disease resistance,and higher yield[END_REF]. Native plants in coastal environments and geothermal habitats require fungal endophytes in order to grow [START_REF] Rodriguez | Stress tolerance in plants via habitat-adapted symbiosis[END_REF]. Thus a passive adaptation of the plant is observed, with the endophytic fungus providing a selective advantage to the colonized plant.

Infection of the tomato (Solanum lycopersicum) plant with these endophytes, for example, confers salt or heat resistance [START_REF] Rodriguez | Stress tolerance in plants via habitat-adapted symbiosis[END_REF]. It can thus be argued that solutions, which support a more productive and sustainable agriculture and involve the use of endophytic microorganisms, do exist but have as yet been little explored.

Concluding remarks

The Green Revolution that started about 50 years ago, allowed food shortages to be limited. Given the stocks of resources and human population growth, this Green Revolution can continue for only a few more decades. The counterpart of this Green Revolution is a high cost to the environment and global environmental changes [START_REF] Tilman | Forecasting agriculturally driven global environmental change[END_REF]. If nothing is done to counteract these changes, thresholds will be exceeded, with dramatic consequences [START_REF] Harte | Human population as a dynamic factor in environmental degradation[END_REF] and indeed the impossibility for natural ecosystems to regenerate. A more sustainable agriculture and a plant neodomestication has to emerge to guarantee food supply over the next 50 years. One way of achieving a more ecologically intensive agriculture would be to consider and protect the ecological functions displayed by AM fungi, which have been effective for more than 400 million years, whatever the ecosystem. This will not only improve natural plant mineral nutrition but also water supply and other ecological functions that have already been clearly documented [START_REF] Smith | Mycorrhizal Symbiosis[END_REF]. Research efforts must also stimulate/accompany this possible plant neodomestication.

Box 1. Arbuscular mycorrhizal symbiosis

Among plant mutualistic symbioses, the arbuscular mycorrhiza relationship has been evolving for more than 400 million years [START_REF] Redecker | Glomalean fungi from the Ordovician[END_REF]. This symbiosis is really mundane and widespread with approximately 80% of land plants colonized by AM fungi [START_REF] Smith | Mycorrhizal Symbiosis[END_REF] across a huge diversity of ecosystems. In this symbiosis, plants provide carbohydrates to the arbuscular mycorrhizal fungi in exchange for minerals, drought resistance and protection against pathogens e.g. [30;32].

The fungus in this mutualistic relationship is an obligate biotroph, its transmission is horizontal and there is no genetic uniformity between fungal symbionts. Several different fungal symbionts colonize the same plant roots.

Box 2. Future of agricultural trajectories guidelines

Forests represent important carbon stocks which, when converted into agrosystems, have a huge impact on CO2 emission to the atmosphere [33] as well as a collateral effect on biodiversity [6;9]. In the context of global changes, it seems fundamental to limit agricultural expansion [START_REF] Foley | Solution for a cultivated planet[END_REF]. The key point seems to be to improve crop yields within existing agrosystems. However, conventional agricultural practices and plant breeding strategies have arguably entered a 'cul-de-sac' because they are "...unlikely to improve attributes already favored by millions of years of natural selection..." [START_REF] Weiner | Evolutionary agroecology, the potential for cooperative, high density, weed-suppressing cereals[END_REF] while underexplored natural keys to crop yield improvement, such as AM fungi, exist but are ignored and maltreated.

To maintain or restore this essential component of soil fertility, conventional agricultural practices need to be modified. The following are suggested guidelines to improve the sustainability of human land use and crop productivity:

(i) Because AM diversity is positively correlated with plant diversity e.g. [START_REF] Van Der Heijden | Mycorrhizal fungal diversity determines plant biodiversity, ecosystem variability and productivity[END_REF],

agriculture will need to make use of greater plant diversity.

(ii) Tillage, if employed, will need to be restricted to maintain hyphal networks and functional efficiency and also to preserve soil aggregates and limit water losses [34].

(iii) Plant breeders should select plants in poor soils, taking into account the two previous aspects, the aim being to maximize the efficiency of AM fungi symbiosis (i.e., plants able to take full advantage of the AM fungi available in soils). These new selected plants might also be able to restore effective AM fungi in the field (iv) Additional mutualist microorganisms such as endophytic fungi should also be considered as important targets to improve plant resistance and productivity.

This should facilitate a passive promotion of AM fungal mutualism and, at the same time, reduce the use of fertilizers, biocides and water. These guidelines have the potential to enhance crop yields and reduce the problems associated with conventional agriculture in both developed and developing countries.

Food production therefore needs

  to be high enough and, at the same time, minimize damage to the environment. This equation cannot be solved with current strategies. Based on recent findings, new trajectories for agriculture and plant breeding which take into account the below-ground compartment and evolution of mutualistic strategy, are proposed in this opinion article. In this context, we argue that plant breeders have the opportunity to make use of native Arbuscular Mycorrhizal symbiosis in an innovative ecologically intensive agriculture.

Acknowledgements

The work is supported by a grant from 'l'Agence Nationale de la Recherche' (ANR-10-STRA-0002). We acknowledge E.T. Kiers and D. Warwick for comments and suggestions of modifications on previous versions of the manuscript.

Hyphal network of arbuscular mycorrhiza. This dense network propagated from the plant roots explores a high volume of soil and capture mineral nutrients and water which are transfered to the roots to the benefit of the host-plant. In return host-plant provides photosynthesized sugars and polylosides to sustain the mutualistic fungal compartment.


In this paper, we propose as a first step a software solution to measure the electrical power consumed in an industrial furnace intended essentially for heat treatments. The soft sensor is constructed from the power physical measurement taken as the output of the set (dimmer + resistances), and the control signal measurement provided by a controller with an unknown structure. The second step consists in a detection of faults like a resistance disconnection, for instance. This phase requires the knowledge of the controller model and the furnace system. An overparametrization method was chosen for the controller estimation. An indirect closed-loop Input-Output (IO) identification approach was used for the furnace model estimation through a Tailor-Made and a decomposition of the closed-loop algorithms. A validation with two other experimental tests concludes the paper.

I. INTRODUCTION

In France, the energy consumption of the industrial sector represents roughly 28% of the french energy consumption and 23% of its total CO 2 emissions [START_REF] Mouel | Fostering Energy Efficiency in manufacturing plants through economical breakthroughs in power and flow rate measurement[END_REF]. This situation is becoming unacceptable in a context of increasingly demands with respect to energy efficiency of industrial and household uses, especially when a great part of this energy can be saved.

That is in this context that the project ANR CHIC (CHaînes de mesures Innovantes à bas Coût) emerged. The objective is to develop and test low-cost sensors to monitor and analyze the energy consumption of major fluids used in industrial sites (electricity, gas, compressed air). These sensors should allow the consumption monitoring and the detection of consumption deviations. The objectives are to achieve a measurement accuracy of about 5%. The project involves the design of new low-cost sensors (both physical and software sensors) in the following areas: current sensors, voltage, power, and gas flow. The work presented in this paper only concernes the study of power measurement.

The objective of the study presented in this paper is, in a first step, to design of a power "soft" sensor. The concept of soft sensors is to combine measurement available or easily achievable, representing the evolution of the studied process, and mathematical models relating the measure and the quantities to determine. This concept is used in different fields and especially in chemical or biological processes [START_REF] Chéruy | Software sensors in bioprocess engineering[END_REF]. Modeling is a key step to determine the quality of the measurement and can be based on a physical or empirical approach, or a combination of both (grey box models). The design of a physical model is excluded since it is poorly adapted to an industrial environment. For this purpose, we propose to build behavioral models.

The second step consists in a diagnosis of faults usually noticed in this kind of equipment. When it is the case, an overshoot in power is observed and may be caused by a door that was poorly closed or by a failure in one of the resistances of the furnace. The diagnosis procedure depends on the kind of occurred failure. We will see that a failure presence can be detected only with the control signal. Therefore, it is important to have a good estimation of the controller. However, this criterion is not sufficient in the sense that we can not distinguish what sort of failure has occurred. In this case, the study of the furnace model parameters such as its gain and time constants is essential.

The I/O furnace modeling in a closed-loop working environment was rarely discussed [START_REF] Bitschnau | Modeling and control of an industrial continuous furnace[END_REF], [START_REF] Stankovski | Industrial furnace control: an experiment in iterative identification and design[END_REF]. There are two types of identification in such case: identification in open-loop and identification in closed-loop. The first one is in practice difficult or simply not feasible. For example, some plants have an integrator or are instable in open-loop, and others can undergo long term drift in open-loop operation and non-stationary disturbances which favor a closed-loop experiment for data collection [START_REF] Landau | Recursive algorithms for identification in closed-loop: A unified approach and evaluation[END_REF]. Also, it has been shown that a closed-loop identification experiment is often the optimal experimental setup [START_REF] Agüro | A virtual closed loop method for closed-loop identification[END_REF]. In the remainder of this paper, two different outputerror algorithms are discussed.

II. PLANT DESCRIPTION

The furnace, shown in figure1, contains a two zones refractory brick chamber. Each zone includes 6 resistances (12kW, 380V) which leads to 144kW of total power. The temperature in each zones is regulated by a digital controller actuated by a dimmer following a setpoint profile imposed by the operator (figure 2), and measured by a thermocouple. Measurements (electric power [kW], furnace temperature [ • C], control signal [%]) are collected each 10s. The closed-loop scheme is summarized in the figure 3. The maximum temperature is of 1600 • C. The figure 4 presents some of the data collected on the furnace during a free-load common test. We can see that power is strongly related to the control signal.

III. POWER MEASUREMENT

The power soft sensor development comes down to model the set dimmer and resistances defined by the relationship P = f (u) between the power and the control signal as shown in the figure 3. Knowing that the furnace is of power equal to 144 kW and that the control signal varies between 0 and 100%, the theoretical model of the dimmer and resistances is then deduced: P (t) = 1, 44 u(t). The plot (Figure 5) of the experimental data P = f (u) reveals a linearity. Thus, a simple leastsquares algorithm (LS) is sufficient to model this relationship and gives : P (t) = 1.506 u(t) -0.9117. I), without taking into account the outlier points, shows that its maximum is less than the fixed objective of 5% with a clearly advantadge for the LS-model. The interest of the estimated model is its simplicity to handle. It can be easily implemented with a microcontroller, for instance. It also gives the possibility to calibrate or recalibrate the model online.

IV. FAULT DIAGNOSIS

An overshoot in power may be caused by a door that was poorly closed or by a fault in one of the resistances of the As it has been said in introduction, this step involves the identification of the controller and the furnace system. Indeed, the feedback is taken into acount by identifying a closed-loop transfer function and determining the normal-operating using the knowledge of the linear controller.

A. Controller estimation

It might be obvious to know exactly the control law, particularly when the systems are numerically controlled. But in practice, the implemented controller can be very different from the real one, because of programming mode, nonlinearities, antisaturation devices, etc. A solution consists to identify the controller. In order to avoid the choice of the controller model, an overparametrization technique is used. Then a criterion based on discrete moment is used to choose the model order [START_REF] Poinot | Contribution à l'identification des systèmes par la méthode de surparamétrisation en traitement des eaux[END_REF].

Consider the system given by the closed-loop setting of the figure 9. y k is the measurement of x k added to measurement

R(z) S(z) H(z) + - + + y k b k x k u k e k Fig. 9. Closed loop noise b k (y k = x k + b k with x k unknown).
Once the measurement is obtained, y k is well known. Thus, y k is a certain value that we can use in a processing control algorithm. The controller identification is therefore a free-noise identification problem. With the prior knowledge of the control signal u k and the setpoint e k , we can apply an ordinary least-squares algorithm [START_REF] Ljung | System identification. Theory for the user[END_REF]. Let us define for example the controller to be estimated by

R(z) S(z) = r 0 + r 1 z -1 1 + s 1 z -1 (1)
Then

u k + s 1 u k-1 = r 0 (e k -y k ) + r 1 (e k-1 -y k-1 ) (2)
Knowing the control signal, the setpoint and the noisy output, we can write ûk = ϕ T k θ

with ϕ T k = e k -y k e k-1 -y k-1 -u k-1 and θT = r0 r1 ŝ1 . If the exact structure is unknown, an overparametrization principle can be used. Thus, C (z) = Rs(z) Ss(z) is chosen such as

deg[R s (z)] > deg[R(z)], deg[S s (z)] > deg[S(z)] (4)
Consider the overparametrized estimated controller

C s (z) = r 0 + r 1 z -1 + . . . + r 2 z -S 1 + s 1 z -1 + . . . s S z -S
where S is the degree of overparametrization, and ûk = ϕ T Sk θS ,with

ϕ T Sk = [e k -y k . . . e k-S -y k-S -u k-1 . . . -u k-S ] θT S = [r 0 . . . r S s 1 . . . s S ]
(5) It is necessary to be sure that the S degree is sufficiently high in order to avoid modeling error and at the same time not uselessly high. For that purpose, a characterization test based on system invariants called "discrete moments" is used [START_REF] Trigeassou | Contribution à l'extension de la méthode des moments en automatique[END_REF].

Let us define the controller impulse response g k with a finite sum in [0, ∞[. A Taylor series expansion of C S (z) in neighborhood of z -1 = 1 gives

Z {g k } = C S (z) = ∞ n=0 (z -1 -1) n n! M n (g k ) (6)
where

M n (g k ) = ∞ k=n A n k g k with A n k = k! (k-n)!
, is the n th order discrete moment of the impulse response g k . In the case of a controller with an integrator, we deduce the moments of order 0, 1 and 2 as follows

   M 0 = β 0 /α 1 M 1 = (β 1 -C 0 α 2 )/α 1 M 2 = 2(β 2 -C 0 α 3 -C 1 α 2 )/α 1 (7)
with

β n = M k=n A n k r k , α n = N k=n
A n k s k and M, N the numerator and denominator orders.

If an overparametrized structure includes the exact structure, then all its moments are equivalent to those of the real system. In practice, the exact moments are unknown, but we can increase the structure complexity and calculate its moments. When a stability of moments is reached, the considered structure includes certainly the real system.

This algorithm is now applied to the experimental data. Because of the lack of knowledge about this structure, an overparametrization to a degree S = 3 was done. In order to find the good S degree, the discrete moments are tested to find the degree from which the system invariants does not change as shown in table II. We can clearly observe a negligible variation of the discrete moment values for models higher or equal to S = 2. Consequently, the estimated controller has the Estimated models 10.

0 th moment 1 st moment 2 nd moment F IT [%] S = 1 2, 2 × 10 -2 5,
Ĉ(z) = r 0 + r 1 z -1 + r 2 z -2 1 + s 1 z -1 + s 2 z -2

B. Furnace model

The parameters are estimated using a nonlinear optimization algorithm which minimizes a quadratic criterion according to the iterative Levenberg-Marquardt algorithm [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]. The criterion is based on the output error, i.e. the error between the measured system output y and the model output ŷ( θ), defined by

ε OE = y -ŷ( θ) = y -H(s, θ)u (8)
where The Levenberg-Marquardt algorithm is defined by the following iterative equation

ŷ( θ) = B(s, θ) Â(s, θ) u(s) (9)
θi+1 = θi -[J ′′ θθ + λI n ] -1 J ′ θ θ= θi ( 10
)
where J ′ θ and J ′′ θθ are the gradient and the approximated Hessian, respectively defined by

J ′ θ = - 2 K K k=1 ε OE ( θ)σ( θ), J ′′ θθ = 2 K K k=1 σ( θ)σ( θ) T (11) σ( θ) is the vector of sensitivity function: σ( θ) = ∂ ŷ( θ) ∂ θ , λ
is a scalar used to control the convergence, and K the number of samples.

The sensitivity functions are computed as follows σ y,âj = -s j ŷ( θi )

A(s, θi ) , σ y, bj = s j u( θi )

A(s, θi ) (12)
The property of asymptotic convergence in this algorithm is achieved at the cost of the minimization of a nonlinear quadratic criteria, which can lead to a local optimum [START_REF] Landau | Unbiased recursive identification using model reference adaptive techniques[END_REF]. To avoid the local optimum, a solution is to use an equationerror algorithm to initialize near the global optimum. The routine IVCTRPM (Instrumental Variable Continuous-Time Reinitialized Partial Moments) can be used [START_REF] Ouvrard | On embedded FIR filter models for identifying CT and discrete-time transfer functions: the RPM approach[END_REF], [START_REF] Tohme | A pseudo-output error algorithm to improve global convergence[END_REF]. Two different methods are compared: the OE identification method based on a closed-loop decomposition (OE-CLD) [START_REF] Grospeaud | Contribution à l'identification en boucle fermée par erreur de sortie[END_REF] and the Tailor-Made first OE method (TM) [START_REF] Grospeaud | Contribution à l'identification en boucle fermée par erreur de sortie[END_REF], [START_REF] Van Donkelaar | Analysis of closedloop identification with a tailor-made parametrization[END_REF].

1) OE-CLD:

We can use a DT controller with a continuoustime (CT) system as shown in Figure 11 [START_REF] Grospeaud | Contribution à l'identification en boucle fermée par erreur de sortie[END_REF]. The solution provided by the OE-CLD consists in simulating y from H(s, θ). The output sensitivity functions σ y, θ should be calculated with taking account of σ u, θ (the sensitivity of the predictive control signal û to θ).

Let us define σ y, θi . We know that

ŷ(s) = H(s, θ)û(s) + - r C(z) u H(s) y b + + + + - C(z) û H(s, θ) - ŷ Optimization Algorithm θ
σ y, θ = H(s, θ) ∂ ûk ∂ θ + ∂H(s, θ) ∂ θ û(s) (13)
with ∂ ûk ∂θi obtained from a difference equation like (2), for example, which gives

∂ ûk ∂ θi + s 1 ∂ ûk-1 ∂ θi = -r 0 ∂ ŷk ∂ θi -r 1 ∂ ŷk-1 ∂ θi in other terms σ ûk , θi + s 1 σ ûk-1 , θi = -r 0 σ y k , θi -r 1 σ y k-1 , θi (14)
With this sensitivity function formulation, it is possible to estimate the CT system model parameters controlled by a digital controller.

2) TM: An OE method is used to globally identify the DT closed loop of the figure 12 and it was widely treated in [START_REF] Grospeaud | Contribution à l'identification en boucle fermée par erreur de sortie[END_REF], [START_REF] Van Donkelaar | Analysis of closedloop identification with a tailor-made parametrization[END_REF] . Let us consider C(z) = R(z) S(z) and H(z, θ) = B(z) Â(z) .

H(z, θ) + -

y M k x k u k r k C(z)
S(z) Â(z) + R(z) B(z) r(z) (15)
or y M (z) P (z) = R(z) B(z)r(z) with

P (z) = [S(z) Â(z) + R(z) B(z)]
The calculation of the sensitivity function

σ k = ∂y M k ∂ θ gives σ k = S(q -1 ) P (q -1 ) ϕ M k ( 16
)
with q is the forward-shift time operator, and

ϕ M k =           -y M (k -1) . . . -y M (k -n a ) û(k -1) . . . û(k -n b )          
The equation ( 16) is used in the OE algorithm. Then, the CT transfer function Ĥ(s) is obtained from the DT one Ĥ(z).

On experimental data, we first initialize the OE algorithm with the IVCTRPM algorithm estimated in open-loop by taking as input the control signal and the measured temperature as output. The results shown in table IV are the estimated parameters of a transfer function of orders n a = 2 and n b = 1 defined by H(s, θ) = G(1+τzs) (1+τp 1 )(1+τp 2 ) . F IT denotes the fitting criterion [START_REF] Ljung | System identification. Theory for the user[END_REF] which is defined by

F IT = 100 × 1 - y(t) -ŷ(t) y(t) -mean(y(t)) (17)
The choice of orders has been made after testing several orders. It can be seen that the IVCTRPM initialized parameters are very close to the TM and OE-CLD estimated ones. However, the difference could be more important for a more disturbed application. The fitting values show that the OE methods improve the quality of estimation. The temperature outputs are drawn in the figure 13.

F This identification is followed by a validation procedure on two other data sets by simulating the estimated closed-loop system (see Figures 14 and15). We can easily observe that the simulated outputs are close to measurement and this result is confirmed by table V.

C. Fault diagnosis results

The furnace model is now estimated and validated. The final step consists of a diagnosis of a resistance disconnection. We have two possibilities:

• First, we can simulate the closed-loop identified below and then deduce the control signal with taking as input the setpoint temperature of the experimental data with a fault. A comparison to the measured control signal

V. CONCLUSION

In this paper, a complete study of a furnace plant was introduced with the objectives of designing a soft power sensor and of a fault diagnosis. It has also been shown that we can easily identify the model using the existing identification approaches. An indirect closed-loop estimation using two OE algorithms were investigated. For that, a previous step to estimate the controller consisting of an overparametrization technique was exploited. A successful validation is also given and the results show that the proposed methods are effective. However, only a resistance disconnection fault diagnosis was investigated but future directions as a door opening or poorly closed fault diagnosis will be explored. This diagnosis deals with building a grey box model in order to obtain a more flexible model which takes into account the mass and the load nature.

Fig. 1 Fig. 3 .

 13 Fig. 1. Experimental plant

Fig. 4 .

 4 Fig. 4. Typical experimental data

Fig. 5 .

 5 Fig. 5. Least-squares dimmer/resistances model estimation comparison between the theoretical model and the LS-model. To verify the consistency of the estimated model, a validation was done on another experimental test as shown in the figures 7 and 8. The calculation of the relative error (see tableI), without taking into account the outlier points, shows that its maximum is less than the fixed objective of 5% with a clearly advantadge for the LS-model.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Comparison between the estimated LS model output and the theoretical dimmer/resistances model output

Fig. 10 . 1 j=0

 101 Fig. 10. Measured and estimated control signals for S = 2 with θ = [â 0 , . . . , âna-1 , b0 , . . . , bn b ] and B(s, θi ) = n b j=0

Fig. 11 .

 11 Fig. 11. OE algorithm based on a closed-loop decomposition

Fig. 12 .

 12 Fig. 12. DT closed-loop system

Fig. 13 .

 13 Fig. 13. Identification results

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Validation on a 2nd free-load test

Fig. 16 .

 16 Fig. 16. Control signal divergence • We simulate all the closed loop which includes the estimated furnace model. With this procedure, a separation can be done between a fault in the set (dim-mer+resistances) (e.g. a resistance disconnection) from a fault in the furnace itself (e.g. a door opening).

TABLE I

 I

	MEAN AND MAXIMUM RELATIVE ERROR WITH ESTIMATION AND
	VALIDATION DATA

  The estimated parameters are given in table III. The measured

	r0	r1	r2	s1	s2
	0, 02	-1, 7 × 10 -3	0.011	-0, 677	-0, 321

TABLE III ESTIMATED

 III PARAMETERS OF THE DT CONTROLLER and estimated control signals are plotted in figure

TABLE V VALIDATION FITTINGS

 VFITTINGS

978-1-4799-0223-1/13/$31.00 ©2013 IEEE

ACKNOWLEDGMENT

This work has been supported by the French Research National Agency (ANR) through Efficacité énergétique et réduction des émissions de CO 2 dans les systèmes industriels program (project CHIC n • ANR-10-EESI-02).