Athmane Bakhta
email: athmane.bakhta@cermics.enpc.fr

Virginie Ehrlacher
email: virginie.ehrlacher@enpc.fr

David Gontier
email: gontier@ceremade.dauphine.fr

Numerical reconstruction of the first band(s) in an inverse Hill's problem

This paper concerns an inverse band structure problem for one dimensional periodic Schrödinger operators (Hill's operators). Our goal is to find a potential for the Hill's operator in order to reproduce as best as possible some given target bands, which may not be realisable. We recast the problem as an optimisation problem, and prove that this problem is well-posed when considering singular potentials (Borel measures). We then propose different algorithms to tackle the problem numerically.

Introduction

The aim of this article is to present new considerations on an inverse band structure problem for periodic one-dimensional Schrödinger operators, also called Hill's operators. A Hill operator is a self-adjoint, bounded from below operator of the form A V := -d 2 dx 2 + V , acting on L 2 (R), and where V is a periodic real-valued potential. Its spectrum is composed of a reunion of intervals, which can be characterised using Bloch-Floquet theory as the reunion of the spectra of a family of self-adjoint compact resolvent operators A V q , indexed by an element q ∈ R called the quasimomentum or k-point (see [22, Chapter XIII] and Section 2.1). The m th band function associated to a periodic potential is the function which maps q ∈ R to the m th lowest eigenvalue of A V q . The properties of these band functions are well-known, especially in the one-dimensional case (see e.g. [22, Chapter XIII]).

The inverse band structure problem is an interesting mathematical question of practical interest, which can be roughly formulated as follows: is it possible to find a potential V so that its first bands are close to some target functions?

A wide mathematical literature answers the question when the target functions are indeed the bands of some Hill's operator, corresponding to some V ref . In this case, we need to recover a potential V that reproduces the bands of V ref . We refer to [START_REF] Eskin | Inverse spectral problem for the Schrödinger equation with periodic vector potential[END_REF][START_REF] Eskin | On isospectral periodic potential in R n . I[END_REF][START_REF] Eskin | On isospectral periodic potential in R n . II[END_REF][START_REF] Pöschel | Inverse spectral theory[END_REF][START_REF] Freiling | Inverse Sturm-Liouville problems and their applications[END_REF][START_REF] Veliev | Multidimensional periodic Schrödinger operator[END_REF] for the case when V ref is a regular potential, and to [START_REF] Hryniv | Inverse spectral problems for Sturm-Liouville operators with singular potentials[END_REF][START_REF] Hryniv | Half-inverse spectral problems for Sturm-Liouville operators with singular potentials[END_REF][START_REF] Hryniv | Inverse spectral problems for Sturm-Liouville operators with singular potentials. II. Reconstruction by two spectra[END_REF][START_REF] Hryniv | Inverse spectral problems for Sturm-Liouville operators with singular potentials. III. Reconstruction by three spectra[END_REF][START_REF] Hryniv | Inverse spectral problems for Sturm-Liouville operators with singular potentials. IV. Potentials in the Sobolev space scale[END_REF] when V ref is singular (see also the review [START_REF] Kuchment | An overview of periodic elliptic operators[END_REF]). The main ideas of the previous references are as follows. First, the band structure of a Hill's operator can be seen as the transformation of an analytic function. In particular, the knowledge of any band on an open set is enough to recover theoretically the whole band structure. A potential is then reconstructed from the high energy asymptotics of the bands.

The previous methods use the knowledge of the behaviour of the high energy bands, and therefore are unsuitable for practical purpose (material design) since we usually have no accurate and numerically stable information about these high energy bands. Moreover, in practice, only the low energy bands are usually of interest. The fact that there exists no explicit characterisation of the set of the first band functions associated to a given admissible set of periodic potentials is an additional numerical difficulty. For applications, it is therefore interesting to know how to construct a potential such that only its first bands are close to some given target functions, which may not be realisable (for instance not analytic). In this present work, we therefore adopt a different point of view, which, up to our best knowledge, has not been studied: we recast the inverse problem as an optimisation problem.

The outline of the paper is as follows. In Section 2, we recall basic properties about Hill's operators with singular potentials. and we state our main result (Theorem 2.3). Its proof is given in Section 3. Finally, we present in Section 4 some numerical tests and propose an adaptive optimisation algorithm, which is observed to converge faster than the standard one. This adaptive algorithm relies on the use of an a posteriori error estimator for discretised eigenvalue problems, whose computation is detailed in the Appendix.

2 Spectral decomposition of periodic Schrödinger operators, and main results

In this section, we recall some properties of Hill's operators with singular potentials. Elementary notions on the Bloch-Floquet transform [START_REF] Reed | Methods of modern mathematical physics. IV: Analysis of operators[END_REF] are gathered in Section 2.1. The spectral decomposition of one-dimensional periodic Schrödinger operators with singular potentials is detailed in Section 2.2, building on the results of [START_REF] Kato | Schrödinger operators with singular potentials[END_REF][START_REF] Hryniv | 1-D Schrödinger operators with periodic singular potentials[END_REF][START_REF] Gesztesy | On spectral theory for Schrödinger operators with strongly singular potentials[END_REF][START_REF] Mikhaelets | One-dimensional Schrödinger operators with singular periodic potentials[END_REF][START_REF] Dias | One-dimensional Schrödinger operators with singular potentials: A Schwartz distributional formulation[END_REF]. We state our main results in Section 2.3.

Bloch-Floquet transform

We need some notation. Let D ′ denotes the Schwartz space of complex-valued distributions, and let D ′ per ⊂ D ′ be the space of distributions that are 2π-periodic. In the sequel, the unit cell is Γ := [-π, π), and the reciprocal unit cell (or Brillouin zone) is

Γ * := [-1/2, 1/2]. For u ∈ D ′ per and k ∈ Z, the k th normalised Fourier coefficient of u is denoted by u(k). For s ∈ R, we denote by H s per := { u ∈ D ′ per , ∥u∥ 2 H s per := ∑ k∈Z (1 + |k| 2 ) s | u(k)| 2 < +∞ }
the complex-valued periodic Sobolev space, which is a Hilbert space when endowed with its natural inner product. We write H s per,r for the real-valued periodic Sobolev space, i.e.

H s per,r := { u ∈ H s per , ∀k ∈ Z, u(-k) = u(k) } .
We also let L Lastly, we denote by C 0 per the space of 2π-periodic continuous functions, and by C ∞ c the space of C ∞ functions over R, with compact support.

To introduce the Bloch-Floquet transform, we let H := L 2 (Γ * , L 2 per ). For any element f ∈ H, we denote by f q (x) its value at the point (q, x) ∈ Γ * × Γ. The space H is an Hilbert space when endowed with its inner product ∀f, g ∈ H, ⟨f, g⟩ H := ˆΓ * ˆΓ f q (x)g q (x)dx dq.

The Bloch-Floquet transform is the map R+x) .

B : L 2 (R) → H defined, for smooth functions φ ∈ C ∞ c (R), by ϕ q (x) := (Bφ) q (x) := ∑ R∈Z φ(x + R)e -iq(
It is an isometry from L 2 (R) to H, whose inverse is given by

( B -1 ϕ ) (x) := ˆΓ * ϕ q (x)e iqx dq = φ(x).
The Bloch theorem states that if A is a self-adjoint operator on L 2 (R) with domain D(A) that commutes with Z-translations, then BAB -1 is diagonal in the q-variable. More precisely, there exists a unique family of self-adjoint operators (A q ) q∈Γ * on L 2 per such that for all φ ∈ L 2 (R)∩D(A),

(Aφ)(x) = ˆΓ * (A q ϕ q )(x) dq.
In this case, we write

A = ˆ⊕ Γ * A q dq.

Hill's operators with singular potentials

Giving a rigorous mathematical sense to a Hill's operator of the formd 2 dx 2 + V on L 2 (R), when the potential V is singular is not an obvious task. In the present paper, we consider V ∈ H -1 per,r , which is a case that was first tackled in [START_REF] Kato | Schrödinger operators with singular potentials[END_REF] (see also [START_REF] Hryniv | 1-D Schrödinger operators with periodic singular potentials[END_REF][START_REF] Dias | One-dimensional Schrödinger operators with singular potentials: A Schwartz distributional formulation[END_REF][START_REF] Gesztesy | On spectral theory for Schrödinger operators with strongly singular potentials[END_REF][START_REF] Mikhaelets | One-dimensional Schrödinger operators with singular periodic potentials[END_REF] for recent results).

The results which are gathered in this section are direct corollaries of results which were proved in these earlier works, particularly in [START_REF] Hryniv | 1-D Schrödinger operators with periodic singular potentials[END_REF].

Proposition 2.1. [Theorem 2.1 and Lemma 3.2 of [START_REF] Hryniv | 1-D Schrödinger operators with periodic singular potentials[END_REF]] For all V ∈ H -1 per,r , there exists

σ V ∈ L 2 per and κ V ∈ R such that V = σ ′ V + κ V in D ′ per . (2.1) Moreover, if a V : H 1 (R) × H 1 (R) → C is the sesquilinear form defined by ∀v, w ∈ H 1 (R), a V (v, w) = ˆR dv dx dw dx + ˆR κ V vw -ˆR σ V ( dv dx w + v dw dx ) , (2.2)
then a V is a symmetric, continuous sesquilinear form on H 1 (R) × H 1 (R), which is closed and bounded from below. Besides, a V is independent of the choice of σ V ∈ L 2 per and κ V ∈ R satisfy- ing (2.1).
Remark 2.2. The expression (2.2) makes sense whenever v, w ∈ H 1 (R). This can be easily seen with the Cauchy-Schwarz inequality, and the embedding H 1 (R) → L ∞ (R). It is not obvious how to extend this result to higher dimension.

A direct consequence of Proposition 2.1 is that one can consider the Friedrichs operator on L 2 (R) associated to a V , which is denoted by A V in the sequel. The operator A V is thus a densely defined, self-adjoint, bounded from below operator on L 2 (R), with form domain H 1 (R) and whose domain is dense in L 2 (R). Formally, it holds that

A V = - ∂ 2 ∂x 2 + V.
The spectral properties of the operator A V can be studied (like in the case of regular potentials) using Bloch-Floquet theory.

The previous result, together with Bloch-Floquet theory, allows to study the operator A V via its Bloch fibers

( A V q )
q∈Γ * . For q ∈ Γ * , it holds that A V q is the self-adjoint extension of the operator

-i d dx + q 2 + V.
It holds that A V q is a bounded from below self-adjoint operator acting on L 2 per , whose form domain is H 1 per , and with associated quadratic form a V q , defined by (recall that H 1 per is an algebra)

∀v, w ∈ H 1 per , a V q (v, w) := ˆΓ [ ( -i d dx + q ) v ( -i d dx + q ) w ] + ⟨V, vw⟩ H -1 per ,H 1 per . (2.3)
In other words, we have

A V = ˆ⊕ Γ * A V q dq.
The fact that L 2 per is compactly embedded in H 1 per implies that A V q is compact-resolvent. As a consequence, there exists a non-decreasing sequence of real eigenvalues

( ε V q,m
) m∈N * going to +∞ and a corresponding orthonormal basis

(u V q,m ) m∈N * of L 2 per such that ∀m ∈ N * , A V q u V q,m = ε V q,m u V q,m . (2.4)
The map Γ * ∋ q → ε V q,m is called the m th band. Since the potential V is real-valued, it holds that A V -q = A V q , so that ε V -q,m = ε V q,m for all q ∈ Γ * and m ∈ N * . This implies that it is enough to study the bands on [0, 1/2]. Actually, we have

σ(A V ) = ∪ q∈[0,1/2] ∪ m∈N * {ε V q,m }.
In the sequel, we mainly focus on the first band. We write ε V q := ε V q,1 for the sake of clarity. Thanks to the knowledge of the form domain of A V q , we know that

ε V q := min v∈H 1 per ∥v∥ L 2 per =1 a V q (v, v). (2.5)
This characterisation will be the key to our proof. When the potential V is smooth (say V ∈ L 2 per ), then the map Γ * ∋ q → ε V q,m is analytic on (-1/2, 1/2). Besides, it is increasing on [0, 1/2] if m is odd, and decreasing if m is even (see e.g. [22, Chapter XIII]).

Main results

The goal of this article is to find a potential V so that the bands of the corresponding Hill's operator are close to some given target functions. In order to do so, we recast the problem as a minimisation one, of the form

V * ∈ arg min V ∈V J (V ).
Unfortunately, we were not able to consider the full setting where the minimisation set V is the whole set H -1 per,r . The problem was that we were unable to control the negative part of V . To bypass this difficulty, we chose to work with potentials that are bounded from below. Such a distribution is necessary a measure (see e.g. [19]). Hence measure-valued potentials provide a natural setting for band reconstruction. We recall here some basic properties about measures.

We denote by M + per the space of non-negative 2π-periodic regular Borel measures on R, in the sense that for all ν ∈ M + per , and all Borel set S ∈ B(R), it holds that ν(S) = ν(S + 2π) ≥ 0, and ν(Γ) < ∞. For all ϵ > 0, from the Sobolev embedding H

1/2+ε per → C 0 per , we deduce that M + per → H -1/2-ε per → H -1
per , where the last embedding is compact. For ν ∈ M + per , we denote by

V ν ∈ H -1
per,r the unique corresponding potential, which is defined by duality through the relation:

∀ϕ ∈ H 1 per , ˆΓ ϕdν = ⟨V ν , ϕ⟩ H -1
per ,H 1 per .

For B ∈ R, we define the set of B-bounded from below potentials

V B := { V ∈ H -1 per,r | ∃ν ∈ M + per , V = V ν -B } ⊂ H -1 per,r .
This will be our minimisation space for our optimisation problem. Note that

V B1 ⊂ V B2 for B 1 ≥ B 2 .
We now introduce the functional J to minimise. First, we introduce the set T of allowed target functions:

T := { b ∈ C 0 (Γ * ), b is even and b is increasing on [0, 1/2] } . (2.6) Of course, for all V ∈ H -1 per,r , it holds that Γ * ∋ q → ε V q ∈ T .
Finally, in order to quantify the quality of reconstruction of a band b ∈ T , we introduce the error functional

J b : H -1 per,r → R defined by ∀V ∈ H -1 per,r , J b (V ) := 1 2 ˆΓ * |b(q) -ε V q | 2 dq = ˆ1/2 0 |b(q) -ε V q | 2 dq. (2.7)
The main result of the present paper is the following.

* := ffl Γ * b(q) dq ∈ R. Then, for all B > 1/4 -b * , there exists a solution V b,B ∈ V B to the minimisation problem V b,B ∈ arg min V ∈V B J b (V ).
(2.8)

The proof of Theorem 2.3 relies on the following proposition, which is central to our analysis. Both the proofs of Theorem 2.3 and Proposition 2.4 are provided in the next section.

Proposition 2.4. Let B ∈ R and let (V n ) n∈N * ⊂ V B . For all n ∈ N * , let ν n ∈ M + per such that V n := V νn -B. Let us assume that the sequence ( ε Vn 0 ) n∈N *
is bounded and such that

ν n (Γ) -→ n→+∞ +∞. Then, up to a subsequence (still denoted n), the functions q → ε Vn q converge uniformly to a constant function ε ∈ R, with ε ≥ 1 4 -B. In other words, there is ε ≥ 1 4 -B such that max q∈[0,1/2] ε Vn q -ε ----→ n→∞ 0.
(2.9)

Conversely, for all ε ≥ 1 4 -B, there is a sequence (V n ) n∈N * ⊂ V B such that (2.9) holds.
This result implies that the first band of the sequence of operators

)
n∈N * that converges uniformly to a constant function. However, as the first band of any Hill's operator must be increasing and analytic, the limit is not the first band of a Hill's operator.

3 Proof of Theorem 2.3 and Proposition 2.4

Preliminary lemmas

We first prove some intermediate useful lemmas before giving the proof of Proposition 2.4 and Theorem 2.3. We start by recording a spectral convergence result.

Proposition 3.1. [Theorem 4.1 [11]] Let (V n ) n∈N * ⊂ H -1 per,r be a sequence such that (V n ) n∈N * converges strongly in H -1 per to some V ∈ H -1 per,r . Then, ∀m ∈ N * , max q∈[0,1/2] ε Vn q,m -ε V q,m ----→ n→∞ 0.
In our case, since we are working with potentials that are measures, we deduce the following result.

Proposition 3.2. Let B ∈ R and (V n ) n∈N * ⊂ V B be a bounded sequence, in the sense sup n∈N ⟨V n , 1 Γ ⟩ H -1 per ,H 1 per < ∞. For all n ∈ N * , let ν n ∈ M + per such that V n = V νn -B.
Then, there exists ν ∈ M + per such that, up to a subsequence (still denoted n), (ν n ) n∈N converges weakly-* to ν in M per , and

(V n ) n∈N * converges strongly in H -1 per to V := V ν -B ∈ V B . Moreover, it holds that ∀m ∈ N * , max q∈[0,1/2] ε Vn q,m -ε V q,m ----→ n→∞ 0.
Proof. The fact that we can extract from the bounded sequence (ν n ) n∈N * a weakly-* convergent sequence in M + per is the Prokhorov's theorem applied in the torus Γ * . The second part comes from the compact embedding M per → H -1 per . The final part is the direct application of Proposition 3.1.

Remark 3.3. This proposition explains our choice to consider measure-valued potentials. Note that a similar result does not hold in the L 1 per setting for instance.

We now give a lemma which is standard in the case of regular potentials V (see [START_REF] Evans | Partial Differential Equations[END_REF]).

Lemma 3.4. Let V ∈ V B for some B ∈ R. The first eigenvector u V q=0 ∈ H 1 per of A V q=0 is

unique up to a global phase. It can be chosen real-valued and positive.

Proof. We use the min-max principle (2.5), and the fact that, for u ∈ H 1 per , the following holds

d dx |u| ≤ d dx u a.e.
We see that if u is an eigenvector corresponding to the first eigenvalue, then so is |u|. We now consider a non-negative eigenvector u ≥ 0, and prove that it is positive. The usual argument is Harnack's inequality. However, it is a priori unclear that it works in our singular setting. To prove it, we write V = V ν -B for ν ∈ M + per , and consider the repartition function F ν of ν: F ν (x) := ν((0, x]). This function is not periodic, but the function

f ν (x) := F ν (x) -ν(Γ) x |Γ| is. Since F ν is an non decreasing, right-continuous function, we deduce that f ν ∈ L ∞ per . Moreover, it holds, in the H -1 per sense, that f ′ ν = V ν -|Γ| -1 ν(Γ) = V + B -|Γ| -1 ν(Γ).
As a result, we see that u is solution to the minimisation problem

u ∈ arg min v∈H 1 per,r ∥v∥ L 2 per =1 { ˆΓ dv dx 2 + ( ν(Γ) |Γ| -B ) -2 ˆΓ f ν ( v dv dx ) } .
There exists λ ∈ R so that the corresponding Euler-Lagrange equations can be written in the weak-form:

div F (x, u, u ′ ) + G(x, u, u ′ ) = 0, with F (x, u, p) = p -f ν u and B(x, u, p) = f ν p + λu.
We are now in the settings of [23, Theorem 1.1], and we deduce that u > 0. The rest of the proof is standard.

Proof of Proposition 2.4

We now prove Proposition 2.4. Let B ∈ R and let

V n = V νn -B ∈ V B with ν n ∈ M +
per , be a sequence such that the sequence

( ε Vn q=0 ) n∈N * is bounded and ν n (Γ) goes to +∞. Since ( ε Vn 0 )
n∈N * is bounded, then up to a subsequence (still denoted by n), there exists ε ∈ R such that ε Vn 0 converges to ε. Our goal is to prove that the convergence also holds uniformly in q ∈ Γ * .

Let u Vn

0 ∈ H 1 per be the L 2 per -normalised positive eigenvector of A Vn 0 associated to the eigenvalue ε Vn 0 (see Lemma 3.4). We denote by α n := min x∈Γ u Vn 0 (x) > 0. Let us first prove that the following convergences hold:

α n ˆΓ u Vn 0 dν n -----→ n→+∞ 0 and α 2 n ν n (Γ) -----→ n→+∞ 0. (3.1)
From the equality ˆΓ d dx

( u Vn 0 ) 2 + ˆΓ |u Vn 0 | 2 dν n = ε Vn 0 + B, we get α 2 n ν n (Γ) ≤ α n ˆΓ u Vn 0 dν n ≤ ˆΓ |u Vn 0 | 2 dν n ≤ ε Vn 0 + B. (3.2)
As the right-hand side is bounded, and ν n (Γ) → +∞ by hypothesis, this implies α n → 0. Moreover, we have

0 ≤ ˆΓ u Vn 0 dν n = a Vn 0 (u Vn 0 , 1 Γ ) + B ˆΓ u Vn 0 = (ε Vn 0 + B) ˆΓ u Vn 0 ≤ (ε Vn 0 + B)|Γ| 1/2 ,
where we used the Cauchy-Schwarz inequality for the last part. As a result, we deduce that the sequence

( ´Γ u Vn 0 dν n ) n∈N *
is bounded. The first convergence of (3.1) follows. The second convergence is a consequence of the first inequality in (3.2).

Let

x n ∈ Γ = [0, 2π) be such that α n = u Vn 0 (x n ). The fact that α n → 0 implies that l n := ∥u Vn 0 (x n + •) -α n ∥ 2 L 2
per → 1 and we can thus define for n large enough

v n := u Vn 0 (x n + •) -α n ∥u Vn 0 (x n + •) -α n ∥ L 2 per . It holds that v n ∈ H 1 per , ∥v n ∥ L 2 per = 1.
Besides, it holds that v n (0) = 0. For q ∈ Γ * , we introduce the function v q,n defined by: ∀x ∈ R, v q,n (x) := v n (x)e -iq [x] , where we set [x] := x mod 2π.

Thanks to the equality v n (0) = 0, it holds that v q,n ∈ H 1per , and that ∥v q,n ∥ L 2 per = 1. This function is therefore a valid test function for our min-max principle 1 .

From the min-max principle (2.5) and the expression (2.3), we obtain

B + ε Vn q ≤ B + a Vn q (v q,n , v q,n ) = ˆΓ ( -i d dx + q ) v q,n 2 + ˆΓ |v q,n | 2 dν n = ˆΓ dv n dx 2 + ˆΓ |v n | 2 dν n = 1 l n ( ˆΓ d dx ( u Vn 0 (x n + •) ) 2 + ˆΓ |u Vn 0 (x n + •) -α n | 2 dν n ) = 1 l n ( ˆΓ d dx ( u Vn 0 ) 2 + ˆΓ |u Vn 0 | 2 dν n -2α n ˆΓ u Vn 0 dν n + α 2 n ν n (Γ) ) = 1 l n ( B + ε Vn 0 -2α n ˆΓ u Vn 0 dν n + α 2 n ν n (Γ)
) .

We infer from these inequalities, and from (3.1) that

0 ≤ max q∈Γ * ε Vn q -ε Vn 0 ≤ ( B + ε Vn 0 ) ( 1 l n -1 ) + 1 l n ( -2α n ˆΓ u Vn 0 dν n + α 2 n ν n (Γ) ) -----→ n→+∞ 0.
This already proves the convergence (2.9).

To see that ε ≥ 1 4 -B, we write, for

V = V ν -B with ν ∈ M + per that ∀q ∈ [-1/2, 1/2], A V q = -i d dx + q 2 + V ν -B ≥ -i d dx + q 2 -B ≥ q 2 -B,
where we used the fact that the lowest eigenvalue of -i d dx + q 2 is q 2 for q ∈ [-1/2, 1/2] (this can be seen with the Fourier representation of the operator). As a consequence, for q = 1 2 , we obtain that for all

V ∈ V B , ε V q=1/2 ≥ 1 4 -B.
The result follows. To prove the converse, we exhibit an explicit sequence of measures

(ν n ) n∈N * ⊂ M + per such that ε Vν n q → 1 4 .
The general result will follow by taking sequences of the form

V n = V νn + ( ε -1 4
) -B. We denote by δ x the Dirac mass at x ∈ R, and consider, for λ > 0, the measure

ν λ := λ ∑ k∈Z δ 2πk ∈ M + per . (3.3)
From the first part of the Proposition, it is enough to check the convergence for q = 0. We are looking for a solution to (we denote by

ω 2 λ := ε Vν λ 0 ≥ 0 for simplicity) -u ′′ + λδ 0 u(0) = ω 2 λ u, u ≥ 0, u(2π) = u(0). (3.4)
On (0, 2π), u satisfies the elliptic equation -u ′′ = ω 2 λ u, hence is of the form

u(x) = Ce iω λ x + De -iω λ x ,
for some C, D ∈ R. The continuity of u at 2π implies Ce 2iπω λ + De -2iπω λ = C + D. Moreover, integrating (3.4) between 0 -and 0 + leads to the jump of the derivative -u ′ (0)+u ′ (2π)+λu(0) = 0,

or iω λ (D -C) + iω λ ( Ce 2iπω λ -De -2iπω λ ) + λ(C + D) = 0.
We deduce that (C, D) is solution to the 2 × 2 matrix equation

( 1 -e 2iπω λ 1 -e -2iπω λ -iω λ ( 1 -e 2iπω λ ) + λ iω λ ( 1 -e -2iπω λ ) + λ ) ( C D ) = ( 0 0
) .

The determinant of the matrix must therefore vanish, which leads to

1 = cos(2πω λ ) + λ 2 sin(2πω λ ) ω λ .
As λ → ∞, one must have ω λ → 1/2, or equivalently ε Vν λ 0

→ 1/4. The result follows.

Proof of Theorem 2.3

We are now in position to give the proof of Theorem 2.

* := ffl Γ * b(q) dq. Let V n = V νn -B ⊂ V B
be a minimising sequence associated to problem (2.8). Let us first assume by contradiction that ν n (Γ) → ∞. Then, according to Proposition 2.4, up to a subsequence (still denoted by n), there exists ε ≥ 1 4 -B such that ε Vn q converges uniformly in q ∈ Γ * to the constant function ε. Also, from the second part of Proposition 2.4, the fact that B > 1 4 -b * and the fact that b * is the unique minimiser to

inf c∈R K b (c), (3.5)
where

K b (c) := ´[0,1/2] |b(q) -c| 2 dq for all c ∈ R, it must hold that ε = b * .
We now prove that inf

V ∈V B J b (V ) ̸ = inf c∈R K b (c) = K b (b * ).
To this aim, we exhibit a potential

W ∈ V B such that J b (W ) < K b (b * ).
Since b is continuous and increasing on [0, 1/2], there exists a unique q * ∈ (0, 1/2) such that b(q * ) = b * . We choose δ > 0 small enough such that 0 < q * -δ < q * + δ < 1/2, and set

η ext := ˆq * -δ 0 |b(q) -b * | 2 dq + ˆ1/2 q * +δ |b(q) -b * | 2 dq and η int := ˆq * +δ q * -δ |b(q) -b * | 2 dq, so that K b (b * ) = η ext + η int .
Since b is increasing and continuous, it holds that η int > 0 and η ext > 0, and that b(q * -δ) < b * < b(q * + δ).

We now choose a constant σ > 0 such that

0 < σ < min { η int 8δ , B + b * - 1 4 , b * -b(q * -δ), b(q * + δ) -b * } .
Let ν n be the measure defined in (3.3) for λ = n ∈ N, and let

W n := V νn + b * - 1 4 .
Since ε Wn q converges to b * uniformly in Γ * , there exists n 0 ∈ N * large enough such that

∀q ∈ Γ * , ε Wn 0 q -b * < σ/2.
We then define

W := W n0 + b * -ε Wn 0 q * = V νn + [( B + b * - 1 4
) -

( ε Wn 0 q * -b * )] -B. Since σ < B + b * -1/4, it holds that W ∈ V B .
Moreover, it holds that b * -σ < ε W q < b * + σ for all q ∈ Γ * . Finally, for q = q * , we have ε W q * = b * .

Let us evaluate J b (W ). We get

J b (W ) = ˆq * -δ 0 |b(q) -ε W q | 2 dq + ˆq * +δ q * -δ |b(q) -ε W q | 2 dq + ˆ1/2 q * +δ |b(q) -ε W q | 2 dq.
For the first part, we notice that for 0 ≤ q < q * -δ, we have

b(q) < b(q * -δ) < b * -σ < ε W q < ε W q * = b * .
This yields that

∀ 0 ≤ q < q * -δ, |b(q) -ε W q | = ε W q -b(q) < b * -b(q) = |b(q) -b * |.
Integrating this inequality leads to ˆq * -δ

0 |b(q) -ε W q | 2 dq < ˆq * -δ 0 |b(q) -b * | 2 dq.
Similarly, we obtain that ˆ1/2

q * +δ |b(q) -ε W q | 2 dq < ˆ1/2 q * +δ |b(q) -b * | 2 dq.
Lastly, for the middle part, we have ˆq * +δ

q * -δ |b(q) -ε W q | 2 dq < 2δ [ ε W q * +δ -ε W q * -δ ] ≤ 4δσ ≤ η int 2 < ˆq * +δ q * -δ |b(q) -b * | 2 dq.
Combining all these inequalities yields that J b (W ) < K b (b * ). This contradicts the minimising character of the sequence

(V n ) n∈N * .
Hence the sequence (ν n (Γ)) n∈N * is bounded. The proof of Theorem 2.3 then follows from Proposition 3.2.

Numerical tests

In this section, we present some numerical results obtained on different toy inverse band structure problems. We propose an adaptive optimisation algorithm in which the different discretisation parameters are progressively increased. Such an approach, although heuristic, shows a significant gain in computational time on the presented test cases in comparison to a naive optimisation approach.

In Section 4.1, we present the discretised version of the inverse band problem for multiple target bands. We present the different optimisation procedures used for this problem (direct and adaptive) in Section 4.2. Numerical results on different test cases are given in Section 4.3. The reader should keep in mind that although the proof given in the previous section only works for the reconstruction of the first band, it is possible to numerically look for methods that reproduce several bands.

Discretised inverse band structure problem

For k ∈ Z, we let e k (x) := 1 √ 2π e ikx be the k-th Fourier mode. For s ∈ N * , we define by

X s := Span {e k , k ∈ Z, |k| ≤ s} (4.1)
the finite dimensional space of L 2 per consisting of the N s := 2s + 1 lowest Fourier modes. We denote by Π Xs : L 2 per → X s the L 2 per orthogonal projector onto X s . In practice, the solutions of the eigenvalue problem (2.4) are approximated using a Galerkin method in X s . We denote by ε V,s q,1 ≤ • • • ≤ ε V,s q,Ns the eigenvalues (ranked in increasing order, counting multiplicity) of the operator A V,s q := Π Xs A V q Π * Xs . We also denote by (u V,s q,1 , • • • , u V,s q,Ns ) an orthonormal basis of X s composed of eigenvectors associated to these eigenvalues so that

∀1 ≤ j ≤ N s , A V,s q u V,s q,j = ε V,s q,j u V,s q,j . (4.2)
An equivalent variational formulation of (4.2) is the following:

∀1 ≤ j ≤ N s , ∀v ∈ X s , a V q ( u V,s q,j , v ) = ε V,s q,j ⟨ u V,s q,j , v ⟩ L 2 per .
As s goes to +∞, it holds that ε V,s q,m -→

s→+∞ ε V q,m .
In order to perform the integration in (2.7), we discretise the Brillouin zone. We use a regular grid of size Q ∈ N * , and set

Γ * Q := { - 1 2 + j Q , j ∈ {0, • • • , Q -1} } .
We emphasise that since the maps q → ε q,m are analytic and periodic, the discretisation error coming from the integration will be exponentially small with respect to Q. In practice, we fix

Q ∈ N * .
Let M ∈ N * be a desired number of targeted bands and b 1 , • • • , b M ∈ C 0 per be real-valued even functions, and such that b m is increasing when m is odd and decreasing when m is even. Our cost functional is therefore J : H -1 per,r → R, defined by

∀V ∈ H -1 per,r , J (V ) := 1 Q ∑ q∈Γ * Q M ∑ m=1 |b m (q) -ε V q,m | 2 .
Its discretised version, when the eigenvalues problems are solved with a Galerkin approximation, is

∀s ∈ N * , ∀V ∈ H -1 per,r , J s (V ) := 1 Q ∑ q∈Γ * Q M ∑ m=1 |b m (q) -ε V,s q,m | 2 .
Recall that our goal is to find a potential V ∈ H -1 per,r which minimise the functional J s . In practice, an element V ∈ H -1 per,r is approximated with a finite set of Fourier modes. For p ∈ N * , we denote by

Y p := Span    ∑ k∈Z, |k| 2 ≤p V k e k , ∀k ∈ Z, |k| ≤ p, V -k = V k    . (4.3)
Altogether, we want to solve V s,p := arg min

V ∈Yp J s (V ).

Algorithms for optimisation procedures

Naive algorithm

We first present a naive optimisation procedure, using a gradient descent method, where the parameters s and p are fixed beforehand. We tested three different versions of the gradient descent algorithm: steepest descent (SD), conjugate gradient with Polak Ribiere formula (PR) and quasi Newton with the Broyden-Fletcher-Goldfarb-Shanno formula (BFGS). We do not detail here these classical descents and corresponding line search routines for the sake of conciseness and refer the reader to [START_REF] Bakhta | Mathematical models and numerical simulation of photovoltaic materials[END_REF][START_REF] Bonnans | Numerical Optimization[END_REF].

For all V ∈ H -1 per,r , there exists real-valued coefficients

( c V k ) k∈N and ( d V k ) k∈N * such that V (x) = c V 0 + ∑ k∈N * c V k cos(kx) + d V k sin(kx), and
∑ k∈N * (1 + |k| 2 ) -1 ( |c V k | 2 + |d V k | 2 ) < +∞.
For all k ∈ N (respectively k ∈ N * ), we can express the derivative

∂ c V k J s (V ) (respectively ∂ d V k J s (V )
) exactly in terms of the Bloch eigenvectors u V,s q,m . Indeed, it holds that

∂ c V k J s (V ) = 1 Q ∑ q∈Γ * Q M ∑ m=1 2 ( ε V,s q,m -b m (q) ) ∂ c V k ( ε V,s q,m
) .

On the other hand, from the Hellman-Feynman theorem, it holds that

∂ c V k ( ε V,s q,m ) = ⟨ u V,s q,m , ∂ c V k A V q , u V,s q,m ⟩ = ⟨u V,s q,m , cos(k•)u V,s q,m ⟩ L 2 per .
Similarly, for all k ∈ N * ,

∂ d V k ( ε V,s q,m ) = ⟨ u V,s q,m , ∂ d V k A V q , u V,s q,m ⟩ = ⟨u V,s q,m , sin(k•)u V,s q,m ⟩ L 2 per .
In the rest of the article, for all p ∈ N * , we will denote by ∇J s (V )| Y p the 2p + 1-dimensional real-valued vector so that

∇J s (V ) Y p = ( ∂ d V p J s (V ), ∂ d V p-1 J s (V ), • • • , ∂ d V 1 J s (V ), ∂ c V 0 J s (V ), ∂ c V 1 J s (V ), • • • , ∂ c V p J s (V )
) .

In order for the reader to better compare our adaptive algorithm with this naive one, we provide its pseudo-code below (Algorithm 1).

Algorithm 1: Naive optimisation algorithm

Input: p, s ∈ N * ; W 0 ∈ Y p : initial guess; ε > 0: prescribed global precision; ν > 0: tolerance for the norm of the gradient;

Output: W * ∈ Y p such that ∥∇J s (W * ) Yp ∥ ≤ ν; Instructions: n = 0, W = W 0 ; while ∥∇J s (W ) Yp |∥ > ν do compute a descent direction D ∈ Y p at J s (W ) (using SD / PR / BFGS); choose t ∈ R so that t ∈ argmin t∈R J s (W + tD); set W ← W + tD; return W * = W .
Although this method gives satisfactory numerical optimisers as shown in Section 4.3, its computational time grows very quickly with the discretisation parameters p and s. Besides, it is not clear how these parameters should be chosen a priori, given some target bands. This motivates the design of an adaptive algorithm.

Adaptive algorithm

In order to improve on the efficiency of the numerical optimisation procedure, we propose an adaptive algorithm, where the discretisation parameters s or p are increased during the optimisation process. To describe this procedure, we introduce two criteria to determine whether s or p need to be increased during the algorithm. As the parameter s is increased, the approximated eigenvalues ε V,s q,m becomes more accurate, and the discretised cost functional J s gets closer to the true one J . Our criterion for s relies on the use of an a posteriori error estimator for the eigenvalue problem (4.2). More precisely, assume we can calculate at low numerical cost an estimator ∆ V,s q,m ∈ R + such that

|ε V m,q -ε V,s m,q | ≤ ∆ V,s q,m ,
(see Appendix A), then we would have that

|J (V ) -J s (V )| = 1 Q ∑ q∈Γ * Q M ∑ m=1 ( |b m (q) -ε V q,m | 2 -|b m (q) -ε V,s q,m | 2 ) = 1 Q ∑ q∈Γ * Q M ∑ m=1 ( 2b m (q) -ε V q,m -ε V,s q,m ) ( ε V,s q,m -ε V q,m ) ≤ 1 Q ∑ q∈Γ * Q M ∑ m=1 ( 2 b m (q) -ε V,s q,m + ∆ V,s q,m ) ∆ V,s q,m =: S s V .
The quantity S s V estimates the error between J (V ) and J s (V ) and therefore gives information on the necessity to adapt the value of the discretisation parameter s.

We now derive a criterion for the parameter p. When this parameter is increased, the minimisation space Y p gets larger. A natural way to decide whether or not to increase p is therefore to consider the gradient of J s , at the current minimisation point W ∈ Y p , but calculated on a larger subspace Y p ′ ⊃ Y p with p ′ > p.

In practice, the natural choice p ′ = p + 1 is inefficient. This is not a surprise, as there is no reason a priori to expect a sudden change at exactly the next Fourier mode. We therefore took the heuristic choice p ′ = 2p. More specifically, we define

P p V := ∇ V J s (V ) Y2p .
Note that this estimator needs to be computed only when V is a local minimum of J s on Y p . When this estimator is larger than some threshold, we increase p so that the new space Y p contains the Fourier mode which provides the highest contribution in

(∇ V J s (V )) Y2p .
The adaptive procedure we propose is described in details in Algorithm 2:

Algorithm 2: Adaptive optimisation algorithm Input: p 0 , s 0 ∈ N * : initial discretisation parameters; W 0 ∈ Y p0 : initial guess; η > 0: global discretisation precision; ν > 0: gradient norm precision;

Output: p ≥ p 0 , s ≥ s 0 : final discretisation parameters; W * ∈ Y p such that ∥∇J s (W * ) Yp ∥ ≤ ν, S s W * ≤ η and P p W * ≤ η; Instructions: n = 0, W = W 0 ; while ∥∇J s (W ) Yp ∥ > ν or S s W > η or P p W > η do while ∥∇J s p (W ) Yp ∥ > ν do compute a descent direction D ∈ Y p at J s (W ) (using SD / PR / BFGS); choose t ∈ R so that t ∈ arg min t∈R J s (W + tD); set W ← W + tD; if S s W > η then set s ← s + 1; else if P p W > η then set p ← arg max p<p≤2p max ( ∂ d V p J s (W ) , ∂ c V p J s (W )
) ;

return W * = W .

Numerical results

In this section, we illustrate the different algorithms presented above.

We consider the case where the target functions come from a target potential V t ∈ Y pt , whose Fourier coefficients are randomly chosen for some p t ∈ N * . We therefore take b m (q) := ε Vt,st q,m , and try to recover the first M functions b m . The numerical parameters are M = 3, Q = 25, ν = 10 -5 , η = 10 -6 and s t = 20. The initial guess is W 0 = 0. The naive algorithms are run with s = s t and p = p t , while the adaptive algorithms start with s 0 = p 0 = 1. In addition, the a posteriori estimator is obtained with s ref = 250 and θ = 0.01 (see Appendix A). All tests are done with the naive and adaptive algorihms, with steepest descent (SD), conjugate gradient with Polak Ribiere formula (PR) and quasi Newton with the Broyden-Fletcher-Goldfarb-Shanno formula (BFGS).

In our first test, we try to recover a simple shifted cosine function (i.e. p t = 1). Results are shown in Figure 1. We observe that the bands and the potential are well reconstructed. We also notice that the adaptive algorithm takes more iterations to converge. However, as we will see later, most iterations are performed for low values of the parameters s and p, and therefore are usually faster in terms of CPU time (see Table 1

below).

In the second test case, we try to recover a more complex potential with p t = 8 (see Figure 2). In this case, all the algorithms reproduce well the first bands, but fail to recover the potential. Actually, we see how different methods can lead to different local minima for the functional J . This reflects the complex landscape of this function.

We end this section by reporting results obtained with the different algorithms, and for different target potential V t ∈ Y pt with p t = 1, 4, 8, 12 (see Table 1 We notice that although the adaptive approach requires more iterations to converge, it is usually faster than the naive one. As we already mentioned, this is due to the fact that most of the iterations are performed with small values of p and s, and are therefore faster. Moreover, we notice that the adaptive algorithms tend to find an optimised potential which p N ≤ p t , i.e. a less oscillatory potential than the target one.

Ackowledgements

The authors heartily thank Éric Cancès, Julien Vidal, Damiano Lombardi and Antoine Levitt for their great help in this work and for inspiring discussions. The IRDEP institute is acknowledged for funding.

A A posteriori error estimator for the eigenvalue problem

We present in this appendix the a posteriori error estimator for eigenvalue problems that we use in Section 4. [START_REF] Bonnans | Numerical Optimization[END_REF]. More details about this estimator are given in [START_REF] Bakhta | An a posteriori error estimator based on shifts for positive hermitian eigenvalue problems[END_REF].

Let H be a finite dimensional space of size N ref and let A be a self-adjoint operator on H. In our case, H is some X s ref for some large s ref ≫ 1, and A = A V,s ref q . The eigenvalues of A, counting multiplicities are denoted by

ε 1 ≤ ε 2 ≤ • • • ≤ ε N ref .
For N ≪ N ref , we consider X N a finite dimensional subspace of H. We denote by Π X N the orthogonal projection on X N , and by A

N := Π X N AΠ * X N . The eigenvalues of A N are denoted by ε N 1 ≤ ε N 2 ≤ • • • ≤ ε N N .
Let us also denote by

( u N m ) 1≤j≤N a corresponding orthogonal basis of X N , so that ∀1 ≤ m ≤ N, A N u N m = ε N m u N m .
We recall that, from the min-max principle, it holds that ε m ≤ ε N m . A certified a posteriori error estimator for the m-th eigenvalue is a non-negative real number

∆ N m ∈ R + such that ε N m -ε m ≤ ∆ N m .
We also require that the expression of ∆ N m only involves the approximated eigenpair ε N m and u N m (and not ε m ).

Proposition A.1. Assume that ε m (resp. ε N m ) is a non-degenerate eigenvalue of A (resp. A N ), and that 0 < ε N m -ε m < dist ( ε N m , σ(A) \ {ε m } ) . (A.1)
Let λ m < ε m . Then there exists δ m > 0 such that, for all 0 ≤ δ < δ m , we have

ε N m -ε m ≤ ⟨ r N m , (A -c δ ) -1 (A -d δ ) (A -c δ ) -1 r N m ⟩ , (A.2)
where we set c δ := ε N m + δ, d δ := λ m + δ, and where r N m :=

( A -ε N m ) u N m is the residual. Proof. Assumption (A.1) implies that ε N m / ∈ σ (A), so that ( A -ε N m ) is invertible. From the fact that ⟨u N m , Au N m ⟩ = ε N m
, and the definition of the residual, it holds that

ε N m -ε m = ⟨ r N m , ( A -ε N m ) -1 (A -ε m ) ( A -ε N m ) -1 r N m ⟩ . (A.3)
Thus, a sufficient condition for (A.2) to hold is that

(A -c δ ) -1 (A -d δ ) (A -c δ ) -1 ≥ ( A -ε N m ) -1 (A -ε m ) ( A -ε N m ) -1 .
Thanks to the spectral decomposition of A, this is the case if and only if,

∀1 ≤ m ≤ N ref , ε m -d δ (ε m -c δ ) 2 ≥ ε m -ε m (ε m -ε N m ) 2 .
Denoting by η :

= dist ( ε N m , σ(A) \ {ε m } ) - ( ε N m -ε m )
, this holds true as soon as δ ≤ δ m := min (ε m -λ m , η). The result follows.

In order to use the left-side of (A.2) as an a posteriori estimator, we need to choose λ m < ε m and δ m > 0. For the choice of λ m , we follow [START_REF] Wolkowicz | Bounds on eigenvalues using traces[END_REF], and notice that

ε m ≥ λ m := µ - ( N ref -m -1 m + 1 ) 1/2 σ,
where we set

µ := 1 N ref Tr A and σ 2 := 1 N ref Tr A 2 -µ 2 .
For the choice of δ m , we chose the simple rule

δ m = θ ( ε N m -κ ) with 0 < θ ≪ 1 and κ ∈ R independent of m.
The real number κ is chosen to be an a priori lower bound of the lowest eigenvalue ε 1 of A. This choice is heuristic in the sense that we cannot guarantee that the assumptions of Proposition A.1 are satisfied. However, the encouraging numerical results we obtain below motivated our choice to use such an estimator (see Section A).

Numerical test To illustrate the efficiency of our heuristic, we tested it to compute the first bands of the Hill's operator A V with

V (x) = 3 ∑ k=-3
Vk e k , where V0 = 2 and V-1 = V-2 = V1 = V2 = 1 + 0.5 i.

The reference operator is A := A V,s ref q with s ref = 250, and the first three bands are computed on the space X s defined in 4.1 with s = 6. We plot in Figure 3 the true error ε V,s q,m -ε V,s ref q,m for m = 1, 2, 3, and the corresponding a posteriori error with κ = 0 and different values of θ (namely θ = 0.1, 0.5, 1). We observe that our estimator is sharp for a large range of θ.

Theorem 2 . 3 .

 23 Let b ∈ T , and denote by b

3 .

 3 Let b ∈ T and B > 1/4 -b * where b

  Convergence of the algorithms

Figure 1 :

 1 Figure 1: Recovery of the cosine potential

  Convergence of the algorithms

Figure 2 :

 2 Figure 2: Recovery of an oscillating potential.

Figure 3 :

 3 Figure 3: Numerical validation of the a posteriori error estimator proposed in Appendix A.

  ( A Vn ) n∈N * , where (V n ) n∈N * satisfies the assumptions of Proposition 2.4, becomes flat.

	Remark 2.5. Here we have a sequence of first bands	(	ε Vn q

Table 1 :

 1 ). In this table, N denotes the number of iterations, s N and p N are the values of the parameters s and p at the last iteration (in particular, for the naive algorithms, we have s N = s t = 20 and p N = p t ). Lastly, for each algorithm algo, we define a relative CPU timeτ algo = t algo t SD ,where t algo is the CPU time consumed by the algorithm algo and t SD is the CPU time consumed by the classical steepest descent. In particular, τ SD = 1. Results for recovery test with different algorithms. Red values are reference values.

	p t	--	BFGS naive adaptive naive adaptive naive adaptive PR SD
	1	τ N	0.259 8	1.176 31	0.929 21	1.320 154	1 24	1.255 90
		s N	20	3	20	4	20	3
		p N	1	3	1	2	1	3
	4	τ N	0.070 54	0.009 1424	0.464 1927	0.281 7091	1 8453	0.259 19095
		s N	20	8	20	7	20	5
		p N	4	5	4	3	4	3
	8	τ N	0.470 553	0.151 1041	1.090 1023	0.144 1515	1 7326	0.519 26783
		s N	20	6	20	7	20	6
		p N	8	4	8	4	8	4
	12	τ N	0.007 765	0.001 2474	0.054 2413	0.004 2727	1 50312	0.044 34865
		s N	20	9	20	9	20	9
		p n	12	8	12	8	12	8

This construction only works in one dimension. We do not know how to construct similar test functions in higher dimension.


A N T R Monteiro

F Garcia-Launay

L Brossard

A Wilfart

J.-Y Dourmad

Effect of precision feeding on environmental impact of fattening pig production

Keywords: fattening pig, precision feeding, environmental evaluation, life cycle assessment Precision Livestock Farming '17  565

Various studies have shown that improved feed efficiency is an efficient lever to reduce the environmental impact of pig production. The development of feeders which allow the distribution of different diets, together with new low-cost animal identification technologies and sensors for high-throughput collection of information on animals (body weight, feed intake), should allow the development in practice of individual feeding of fattening pigs in the near future. With this in mind, a modelling approach was used to evaluate the potential of this strategy to reduce the environmental impact of pig fattening. Eight populations of 2000 pigs, fed according to either a conventional two-phase (2P) or a precision feeding (PR) strategy, were simulated using the population version of InraPorc. This was performed for two geographical production contexts (Brazil and France) with two soyabean meal origins (Centre West Brazil, a region with deforestation, and South Brazil, without deforestation). Environmental impacts were evaluated using a cradle-to-farm gate life cycle assessment and were expressed per kg of weight gain. Compared to 2P, PR improved average daily gain and feed efficiency by about 4-5% and reduced feed costs by 7-10%. It reduced nitrogen excretion by about 20%. On average, compared to 2P, the PR feeding strategy reduced environmental impacts for climate change, energy use, eutrophication, acidification and land occupation by 6.1%, 6.9%, 10.3%, 12.7% and 3.5%, respectively. It was concluded that PR appears to be a very interesting strategy for improving economic and environmental dimensions of sustainability in pig production.

Introduction

The environmental impacts of pig production have been debated increasingly in recent years, resulting in a greater focus on identifying and mitigating the environmental degradation they may cause. Better adjustment of nutrient supply to animal requirements [START_REF] Dourmad | Impact of nutrition on nitrogen, phosphorus, Cu and Zn in pig manure, and on emissions of ammonia and odours[END_REF]) may be a key factor in improving the efficiency of nutrient retention, reducing excretion and consequently increasing the sustainability of pig production. The development of innovative feeders which allow the distribution of different diets, together with new low-cost animal identification technologies and sensors for high-throughput collection of animal information (body weight, feed intake), should allow the development of individual pig feeding in practice in the near future.

Precision feeding appears to be the most promising approach to improving the efficiency of nutrient use. This technique has been successfully tested on growing-finishing pigs [START_REF] Pomar | Applying precision feeding techniques in growing-finishing pig operations[END_REF][START_REF] Andretta | The impact of feeding growing-finishing pigs with daily tailored diets using precision feeding techniques on animal performance, nutrient utilization, and body and carcass composition[END_REF]. The strategy is to predict the expected performance of each individual pig according to previous real-time measurements of feed intake and body weight, and to feed it with a ration which provides the amount of nutrient required to achieve the expected performance [START_REF] Hauschild | Development of sustainable precision farming systems for swine: estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF].

In recent years, life cycle assessment (LCA) has been widely used in agriculture [START_REF] Guinée | Handbook on life cycle assessment: Operational guide to the ISO standards[END_REF] and several studies of swine production chains have been conducted [START_REF] Nguyen | Fossil energy and GHG saving potentials of pig farming in the EU[END_REF][START_REF] Garcia-Launay | Evaluation of the environmental implications of the incorporation of feeduse amino acids in pig production using Life Cycle Assessment[END_REF]. LCA makes it possible to perform an integrated environmental evaluation of the whole production chain, considering the potential impacts associated with the raising of pigs as well as those related to the production of inputs (e.g. the feed) and the disposal of waste (e.g. the manure produced).

The aim of this study was thus to use a simulation approach to evaluate the potential for precision feeding of fattening pigs, compared to phase feeding, in two different production contexts, in southern Brazil and western France.

Materials and methods

Feeding strategies and feed specifications

The study considered the growing-finishing pig production system, with three different feeding programmes: two phases (2P), four phases (4P) and precision feeding (PR). Two experimental feeds (named A and B) were independently formulated at least cost for each sex (females and castrated males) in each country. Feeds A and B differed in their amino acid and mineral concentrations, with feed A being formulated with a high nutrient density to meet the estimated nutrient requirements at the beginning of the growing period, and feed B formulated with a low nutrient density to meet the estimated nutrient requirements at the end of the finishing period. Feeds A and B were combined according to the feeding programme. In the 2P feeding programme, feeds A and B were combined in the following proportions: 100% and 0% in phase 1, and 42% and 58% in phase 2 for feeds A and B, respectively. In the 4P feeding programme, feeds A and B were combined in the following proportions: 100% and 0% in phase 1, 65% and 35% in phase 2, 42% and 58% in phase 3 and 15% and 85% in phase 4 for feeds A and B, respectively. In the precision feeding strategy, the proportion of feed A and B was calculated each day for each pig according to its expected amino acid requirement [START_REF] Hauschild | Development of sustainable precision farming systems for swine: estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF]. Feeds contained the usual feed ingredients, including industrial amino acids, used in each country and were formulated at least cost using OpenSolver for Excel®, using the mean prices of feed ingredients in Brazil and France for the year 2014.

Pig production

Performance data from experimental studies in Brazil [START_REF] Monteiro | The impact of feeding growing-finishing pigs with reduced dietary protein levels on performance, carcass traits, meat quality and environmental impacts[END_REF] and France [START_REF] Brossard | Comparison of in vivo and in silico growth performance and variability in pigs when applying a feeding strategy designed by simulation to control the variability of slaughter weight[END_REF] were used to adjust average animal profile parameters for growth and feed intake using InraPorc® software. These profiles were used to calculate, according to InraPorc®, the average nutritional requirement curves for each sex (females and castrated males), these requirements being used for diet formulation. To take account of the variability, the nutrient requirement of the population was calculated as 110% of the mean requirement as generally recommended [START_REF] Pomar | Applying precision feeding techniques in growing-finishing pig operations[END_REF]Brossard et al., 2009). Parameters for growth and feed intake profiles were thus defined for a population of 1000 castrated males and 1000 females for each country, according to the method described by Brossard et al. (2009) using a variance-covariance matrix. Simulations for 2000 pigs (50% females, 50% castrated males) were performed for each feeding scenario in each country in order to determine animal performance, and nutrient balance and excretion.

Life Cycle Assessment LCA was performed according to [START_REF] Nguyen | Fossil energy and GHG saving potentials of pig farming in the EU[END_REF] and [START_REF] Garcia-Launay | Evaluation of the environmental implications of the incorporation of feeduse amino acids in pig production using Life Cycle Assessment[END_REF]. The LCA considered the entire pig farming activity, including crop production, grain drying and processing, production and transport of feed ingredients, feed production at the feed factory, transport of the feed to the farm, growing to finishing pig production, and manure storage, transport and spreading. The pig production system considered was typical of conventional growing-finishing pig farms located in Brittany and southern Brazil. The environmental consequences of manure utilisation were evaluated using system Precision Livestock Farming '17  567 expansion as described by [START_REF] Nguyen | Fossil energy and GHG saving potentials of pig farming in the EU[END_REF]. It was assumed, for both countries, that soyabean was produced in Brazil either in the Centre West (CW), a region with recent deforestation, or in the South (SO) where there is no deforestation. For Brazilian crops, the life cycle inventory (LCI) came from Prudêncio da Silva et al. ( 2010), taking into account the effect of land-use change on carbon release due to conversion of Brazilian forest to cropland. For French crops, the LCI came from a national database developed by French research institutes with data for the environmental impacts of all main ingredients used in animal feeds [START_REF] Wilfart | LCIA results of feedstuffs for French livestock[END_REF]. Emissions to air during swine production and management of manure were estimated according to [START_REF] Rigolot | Modelling of manure production by pigs and NH 3 , N 2 O and CH 4 emissions. Part II: effect of animal housing, manure storage and treatment practices[END_REF] and IPCC (2006). The following potential impacts of pig production were considered: climate change (CC, kg CO 2 -eq.), eutrophication potential (EP, g PO 4 -eq.), acidification potential (AP, g SO 2 -eq.), terrestrial ecotoxicity (TE, g 1,4-DCB-eq.), cumulative energy demand (CED, MJ), and land occupation (LO, m 2 .year). The CC was calculated according to the 100-year global warming potential factors in kg CO 2 -eq. Impacts were calculated at the farm gate and the functional unit considered was one kg of body weight gain (BWG) over the fattening period.

Results and discussion

Animal performance Simulated pig performance and the results for nitrogen excretion are presented in Table 1. Average dietary crude protein (CP) content was significantly affected by country and feeding programme. Compared with 2P, the PR feeding programme reduced dietary CP by 7% and 13% for Brazil and France, respectively. The reduction of CP content was accompanied by a decrease in feeding cost (-7 and -8% between 2P and PR programmes for Brazil and France, respectively). In both countries, ADG was affected by the feeding programme: the highest growth performance was obtained for PR (946 g/d on average) and the lowest (914 g/d on average) for 4P. Similar effects were observed for feed conversion ratio, which was lower for PR than for 4P (2.66 vs. 2.75 kg/kg). Compared with the 2P programme, the PR programme reduced nitrogen excretion by 19 and 21% in Brazil and France, respectively, 4P feeding programmes being intermediate.

Acidification and eutrophication potential

With soyabean from SO, the AP values for the different feeding programmes ranged from 53.9 to 60.2 g SO 2 -eq. per kg BWG in Brazil and from 40.8 to 47.1 g SO 2 -eq. per kg BWG in France (Table 1). With soyabean from CW the values increased slightly to 54.7 to 61.3 g SO 2 -eq. per kg BWG in Brazil and 41.3 to 48.0 g SO 2 -eq. per kg BWG in France. The lowest AP impact was obtained for PR, both for soyabean from SO and CW, with on average about a 12% reduction in AP impact for PR compared to 2P and 4P, which did not differ from each other. The EP values for the feeding programmes ranged from 15.6 to 17.5 g PO 4 -eq. per kg BWG, with no difference between countries and with similar results for soyabean from SO and CW (Table 1). The lowest EP impact was obtained for PR in both countries (mean of 15.6 g PO 4 -eq. per kg BWG) and the highest for 2P and 4P which did not differ (mean of 17.2 g PO 4 -eq. per kg BWG). Since nitrogen contributes to eutrophication and to acidification through ammonia emissions [START_REF] Guinée | Handbook on life cycle assessment: Operational guide to the ISO standards[END_REF], the AP and EP impacts were reduced in both countries by increasing the number of feeding phases and to a higher extent by precision feeding. This was not surprising because these strategies reduce nitrogen excretion and, consequently, also reduce NH 3 emissions from animal housing, manure management and field application. The distribution of EP impact among pigs is presented in Figure 1 which shows a non-normal distribution for phase feeding, indicating that the highest performing pigs (i.e. those with the lowest impact) are not able to achieve their expected performance, probably because of insufficient amino acid supplies during the Precision Livestock Farming '17  569 transition periods. With the change in diet composition, this explains the lower average EO impact obtained with precision feeding.

Climate change potential

With soyabean from SO, the average CC values for the different feeding programmes ranged from 2.30 to 2.40 kg CO 2 -eq. per kg BWG in Brazil and from 2.32 to 2.33 kg CO 2 -eq. per kg BWG in France (Table 1). When soyabean meal from CW was used, CC values increased to 2.61 to 2.81 kg CO 2 -eq. per kg of BWG in Brazil and to 2.48 to 2.65 kg CO 2 -eq. per kg BWG in France. The impact on CC did not differ between the 2P and 4P feeding programmes.

Conversely, the PR feeding programme significantly reduced the CC impact, with values on average 3% and 6% lower than for the two other strategies, for soyabean from SO or from CW, respectively. The effect of the PR feeding programme on CC was thus more pronounced with soyabean from CW. This indicated that the effect of phase feeding on CC may depend on the origin and the amount of soyabean, in agreement with the studies of [START_REF] Eriksson | Environmental systems analysis of pig production -the impact of feed choice[END_REF] and Garcia-Launay et al. ( 2014).

Brazil France

Figure 1. Frequency of distribution of pigs (%) and normal curve for eutrophication potential, according to the feeding programme and country (soyabean from SO).

Cumulative energy demand

With soyabean from SO, the CED values for the different feeding programmes ranged from 12.9 to 13.4 MJ-eq. per kg of BWG in Brazil and from 12.1 to 12.7 MJ-eq. per kg of BWG in France. With soyabean from CW, the values increased to 14.5 to 15.5 MJ-eq. per kg BWG in Brazil and to 12.7 to 13.8 MJ-eq. per kg of BWG in France. The lowest CED impact was calculated for PR, with a 4% and 7% reduction in CED impact compared to 2P and 4P, for soyabean from SO and CW, respectively. The possibility of reducing the CED impact by increasing the number of feeding phases was confirmed for diets based on soyabean meal from CW, but not for soyabean from SO. Precision feeding only resulted in reduced CED impact in that situation.

Land occupation

The LO values for the different feeding programmes ranged from 2.22 to 2.39 m 2 .year per kg BWG in Brazil and from 3.77 to 3.82 m 2 .year per kg BWG in France, with similar results for soyabean from both origins (Table 1). The feeding programme affects

LO in Brazil, with the lowest impact for PR which was 6% lower than 2P and 4P, but not in France.

Conclusions

The results of this study indicate that precision feeding would be the most efficient approach for reducing the life cycle impact of pig fattening, whereas the potential of phase feeding programmes depends on the impact considered, soyabean origin and the geographical context of pig production. The benefit of phase feeding for reducing the climate change impact is limited with soyabean from South Brazil, whereas it appears to be an efficient strategy with soyabean from the Centre West. Conversely, potential eutrophication and acidification impacts are largely reduced by phase feeding in a rather similar way in all situations.

Table 1 .

 1 Performance and environmental impacts of growing-finishing pigs fed with different feeding strategies, in 2-phase (2P), 4-phase (4P) or individual precision feeding (PR)

			Brazil 1			France 1		Stat Sign. 2
		2P	4P	PR	2P	4P	PR	C	S	CxS
	CP, g/kg	145 a	139 b	135 b	158 a	149 b	137 c	**	*** ***
	ADG, g/d	842 ab	830 b	878 a	909 ab	899 b	915 a	***	*
	FCR, kg/kg	2.75 b	2.79 a	2.64 c	2.69 b	2.72 a 2.69 b	**	*** ***
	Cost, €/kg ADG	0.54 a	0.53 a	0.50 b	0.60	0.58	0.55	*** ***
	N exc., kg/pig	3.39 a	3.25 a	2.75 b	3.69 a	3.43 a 2.91 b *** ***
	LCA per kg BWG (SO) 3
	CC, kg CO 2 eq. CED, MJ	2.40 a 13.4 b	2.43 a 13.5 a	2.30 b 12.9 b	2.33 ab 2.35 c 12.7 a 12.6 a 12.1 b *** *** 2.32 ** ***	*
	AP, g SO 2 eq. EP g PO 4 eq. LO, m 2 .year	60.2 a 17.5 a 2.39 a	59.6 a 17.3 a 2.40 a	53.9 b 15.7 b 2.26 b	47.1 a 17.2 a 3.78	45.3 a 40.8 b *** *** 16.8 a 15.7 b *** 3.82 3.79 *** **	**
	LCA per kg BWG (CW) 3
	CC, kg CO 2 eq. CED, MJ	2.81 a 15.5 a	2.80 a 15.4 a	2.61 b 14.5 b	2.65 a 13.8 a	2.61 a 2.48 b *** *** 13.5 a 12.7 b *** ***
	AP, g SO 2 eq. EP g PO 4 eq. LO, m 2 .year	61.3 a 17.4 a 2.35 a	60.5 a 17.2 a 2.36 a	54.7 b 15.6 b 2.22 b	48.0 a 17.1 a 3.75	46.0 a 41.3 b *** *** 16.8 a 15.6 b *** 3.80 3.77 *** **	**
	1 Within-country means followed by same or no letter do not differ (P > 0.05)
	2 Statistical significance of effects of Country (C), Feeding Strategy (S) and their
	interaction (CxS), * P<0.05., ** P< 0.01, *** P<0.001
	3 Life Cycle Assessment with soyabean from Centre West (CW) or South (S) Brazil

 Precision Livestock Farming '17

 Precision Livestock Farming '17

Precision Livestock Farming '17  573

Acknowledgements

This project was conducted within a sandwich PhD thesis programme supported by the CAPES Foundation, Brazil.

Precision Livestock Farming '17  571


J Y Dourmad

L Brossard

C Pomar

J Pomar

P Gagnon

L Cloutier

Development of a decision support tool for precision feeding of pregnant sows

Keywords: sow, gestation, precision feeding, decision support tool

Nutritional studies indicate that nutrient requirements for pregnancy largely differ among sows and according to the stage of pregnancy, whereas in practice the same diet is generally fed to all sows from a given herd. In this context, the availability of new technologies for high throughput phenotyping of sows and their environment, and of innovative feeders that allow the distribution of different diets, offers opportunities for a renewed and practical implementation of prediction models of nutrient requirements, in the perspective of improving feed efficiency and reducing feeding costs and environmental impacts. The objective of this study was thus to design a decision support tool that could be incorporated in automated feeding equipment. The decision support tool was developed on the basis of InraPorc® model. The optimal supply for a given sow is determined each day according to a factorial approach considering all the information available on the sow: genotype, parity, expected prolificacy, gestation stage, body condition (i.e. weight and backfat thickness), activity and housing (i.e. type of floor and ambient temperature). The approach was tested using data from 2500 pregnancies on 540 sows. Energy supply was calculated for each sow to achieve, at farrowing, a target body weight established based on parity, age at mating and backfat thickness (18 mm). Precision feeding (PF) with the mixing of two diets was then simulated in comparison with conventional (CF) feeding with a single diet. Compared to CF, PF reduced protein and aminoacid intake, N excretion and feeding costs. At the same time, with PF, amino acid requirement was met for a higher proportion of sows, especially in younger sows, and a lower proportion of sows, especially older sows, received excessive supplies.

Introduction

Nutrient requirements for pregnant sows largely differ among animals according to their body condition at mating, their parity, their expected reproductive performance, their stage of pregnancy, their physical activity and the housing conditions [START_REF] Dourmad | InraPorc: a model and decision support tool for the nutrition of sows[END_REF]. In practice, the feeding level of pregnant sows is to some extent adapted to take account of these variations, but generally the same diet is fed to all sows from a given herd. Moreover, the group housing of pregnant sows, for welfare issues, makes it sometime difficult to feed each animal individually, especially when sows are raised in small groups with a common feeding trough. Conversely, the group housing of pregnant sows also favoured the development of innovative technologies allowing the distribution of feed individually, for sows raised in large groups, using automated electronic feeders and animal identification.

The objective of this study was thus to design a decision support tool to be incorporated in automated feeding equipment, as already developed for fattening pigs [START_REF] Cloutier | Evaluation of a method estimating real-time individual lysine requirements in two lines of growing-finishing pigs[END_REF], and to test it using a set of data already available from an experimental farm.

Development of the decision support tool

Description of the general approach

The description of the general approach is illustrated in Figure 1. The final objective is to send a command to the automated feeder to proceed with feed distribution. This command informs the feeder with the amount of each of the different diets, generally two diets differing in their nutrient content, to be fed to a given sow over a given day or period.

Figure 1. General description of the approach

To take that decision, the decision support system (DSS) uses information relative to the individual sow to be fed, her housing conditions and the general feeding strategy in the farm. This information is stored in a herd database that will not be described in detail in this paper. Different types of information are stored in this database: (i) the description of herd profile and performance, and the general strategy of management of sow body reserves, as described in InraPorc tool [START_REF] Dourmad | InraPorc: a model and decision support tool for the nutrition of sows[END_REF], (ii) information about each individual sow at mating, especially their age, parity, body weight (BW) and backfat thickness (BF), and their history of performance, and (iii) real time data collected either automatically by different sensors about the sows (e.g. BW, physical activity, feeding or drinking activity…) or their environment (ambient temperature, humidity…).

From the available information, which may vary according to the equipment available on the farm and the data management system, the DSS build the "best guess" decision to be transmitted to the automated feeder. This involves mainly two steps: (i) the determination of energy, amino acid and mineral requirements and (ii) the determination of the amount and composition of the ration to be fed. This ration is prepared from the mixing of different diets, generally two diets, available in the automated feeder.

Determination of energy and amino acid requirements

Energy and nutrient requirements are determined according to a factorial approach. Metabolizable energy (ME) requirement is calculated as the sum of the requirements for maintenance, physical activity and thermoregulation, growth and constitution of body reserves, and development of foetuses and uterine contents (Table 1,eq. 1). depending on housing conditions (eq. 4 and 5). The cumulated amount of energy retained in sow body reserves over pregnancy (ER m ) is determined according to the amount of energy in maternal body at mating and the targeted amount after farrowing. These amounts are calculated from sow BW and BF according to the equations proposed by Dourmad et al (1997) (eq. 7). The corresponding metabolizable energy (ME) requirement is calculated from ER m and the efficiency of energy retention in maternal tissues (k m ). ME requirement for conceptus is calculated according to energy retention in conceptus (ER c , eq. 6) and the efficiency of energy retention in conceptus (k c ).

Total nitrogen retention (eq. 9) is calculated as the sum of N retained in conceptus (eq. 8) and the nitrogen retained in maternal tissues. Standardized ileal digestible (SID) lysine requirement is the calculated assuming 6.5% lysine in retained protein (NR x 6.25) with an efficiency of retention of 65% (eq.10).

Utilisation of the decision support system (DSS)

Description of the database A database issued from an experimental farm (Table 2) was used to simulate the utilization of the DSS. This database contains the data from 2511 gestating sows with information about their body condition at mating (i.e. body weight, BW, and backfat thickness, P2) and at farrowing. These data were used to calibrate InraPorc parameters for this phenotype. Litter size at farrowing averaged 14.1 for a mean piglet birth weight of 1.48 kg. Sows BW at mating increased from 163 to 251 kg from parity 1 to parity 8, whereas P2 tended to be higher in first and second parity and then remained rather constant. Determination of energy and lysine requirement A target of maternal BW (i.e. total BW minus uterus contents) was calculated for each sow according to her age and BW at mating (Table 2). Target of P2 at farrowing was fixed to 18 mm for all sows according to the usual recommendation for this farm.

The DSS was then used to calculate the average ME and feed requirement over pregnancy (Table 3). Average ME requirement varied according to parity, from 31 to 36.8 MJ /d, and it was highly variable between sows with a coefficient of variation of about 7%. Figure 2. Evolution of SID lysine requirement (g/kg feed) of primiparous (mean), multiparous (mean) of all sows (mean ± 2sd) sows, and minimum and maximum requirement, according to gestation stage.

The dynamic of SID lysine requirement, expressed in g/kg feed, according to gestation stage is presented in Figure 2. Average SID lysine requirement increases with gestation stage with a large variability among sows, the highest value being 3-fold higher that the lowest. The requirement is also affected by sows parity with much higher values in primiparous than in multiparous sows.

Evaluation of a precision feeding strategy These simulated data were used to evaluate the interest of precision feeding. A conventional 1-phase feeding strategy (CF) was compared to a precision feeding (PF) strategy consisting in the mixing of two diets with either a low (L) or a high (H) nutrient content. SID lysine content was assumed to 4.8, 3.0 and 6,0 g/kg feed and protein content to 14%; 9% and 16% in diets CF, L and H, respectively.

On average the level of incorporation of L diet in the PF strategy was 84%, the value being lower in first parity sow (67%). The level or incorporation of L diet decreased with gestation stage from almost 100% in the first week to less than 30% in the last week. Compared to the 1-phase strategy, PF strategy resulted in 27% decrease in total SID lysine supply and 28% decrease in total crude protein supply. Moreover, the proportion of sows that were underfed in the last two weeks of lactation decreased from more than 60% with 1CF to less than 5% with PF. For first parity sow the difference was even more marked, with almost all primiparous sows receiving deficient diets over the last 10 days of pregnancy with CF, compared to about 10% with PF (Figure 3). Conversely the proportion of sows that were overfed was drastically reduced (Figure 3).

Conclusions

The results from this study indicate that, in the same way as in fattening pigs [START_REF] Pomar | Applying precision feeding techniques in growing-finishing pig operations[END_REF], precision feeding of gestating sows appears a win-win strategy which allows improving nutritional supplies of sows whilst reducing total protein supply and consequently reducing N excretion. The effect on Precision Livestock Farming '17  591 feeding cost was not evaluated but it may be expected that it will also be reduced [START_REF] Dourmad | Nutrition-Environment interaction in the sow : interest of modelling approach[END_REF]. The DSS developed in this study allows adapting the amount and composition of feed to each sow according to her body condition at mating and expected prolificacy, and to stage pregnancy. This DSS will also allow taking account of information collected by sensors during gestation, such as BW, backfat thickness or physical activity, on the environment, such as ambient temperature.

Figure 3 .

 3 Figure 3. Effect of the feeding strategy (1-phase vs precision feeding) during gestation on the proportion of sows (among all sows or among first-parity sows) that received adequate, deficient or in excess lysine supplies.

Table 1 .

 1 Main equations describing nutrient utilization (from[START_REF] Dourmad | InraPorc: a model and decision support tool for the nutrition of sows[END_REF]

	Energy utilisation ME = ME m + ER c / k c + ER m / k m	[1]
		ME m :	ME for maintenance
		ER c :	energy retention in conceptus
		ER m :	energy retained in maternal tissues
		k c = 0.50	efficiency of ME retention in conceptus
		k m = 0.77 efficiency of ME for maternal
	ME for	in thermoneutral conditions
	maintenance	ME m = 440 kJ.BW -0.75 .d -1 (for 240 min.d -1 standing activity)	[2]
	and effect of	physical activity = 0.30 KJ. kg BW -0,75 .d -1 .min -1 standing	[3]
	activity and	below lower critical temperature (LCT)
	ambient	In individually housed sows
	temperature	LCT = 20°C and HP increases by 18 kJ.kg BW -0.75 .d -1 .°C -1 [4]
		In group-housed sows
		LCT = 16°C and HP increases by 10 kJ.kg BW -0.75 .d -1 .°C -1 [5]
	Energy retention	ERc : Energy in conceptus (kJ)
		Ln(ERc) = 11.72 -8.62 e -0,0138 t + 0.0932 Litter size	[6]
		ERm : Energy in maternal tissues (MJ)
		ER m = 13.65 BW gain + 45.94 BF gain	[7]
	Nitrogen retention NR : total N retention (g.d -1 ),
		NR c : N in conceptus (g)
		Ln(6.25 NR c ) = 8.090 -8.71 e -0,0149 t + 0.0872 Litter size	[8]
		NR = 0.85 (d(NR c)/dt -0.4 + 45.9 (t/100) -105.3 (t/100) 2 +
		64.4 (t/100) 3 + a (ME -ME mm ) )	[9]
		where a = f(BW at mating) and ME mm = ME m at mating
	SID lysine
	requirement g/d	SID Lys = 0.036 x BW -0.75 + 6.25 NR x 0.065 /0.65	[10]
	ME requirement for maintenance is calculated according to BW (eq. 2) and
	modulated according to physical activity (eq. 3), and ambient temperature

Precision Livestock Farming '17  587

Table 2 .

 2 Description of the data base BW sows net body weight; P2 sows backfat thickness 2 A target of BW is calculated for each sow according to BW and age at mating. The same target of P2 is used for all sows.

			Average	Average	Average	Average target
		n°	Litter	Piglet	at mating	after farrowing 2
	Parity sows	size	weight, g BW 1 , kg P2 1 , mm	BW 1 , kg	P2 1 , mm
	1	392	13.3	1405	163	16.9	203	18
	2	389	13.5	1557	192	15.9	227	18
	3	413	14.1	1523	211	15.0	243	18
	4	384	14.9	1480	227	14.4	255	18
	5	335	15.0	1472	234	14.1	260	18
	6	253	14.8	1438	241	14.1	263	18
	7	187	13.9	1445	246	14.6	265	18
	8	158	13.6	1455	251	14.9	267	18
	all	2511	14.1	1478	214	15.2	244	18

588  Precision Livestock Farming '17 1

Table 3 .

 3 Calculated ME (MJ/d) and SID lysine requirement (g/kg feed) and supplies per parity, and % of low nutrient density feed (L) and % of reduction of lysine in precision feeding (PF) compared to conventional feeding (CF) strategy.

				Av. lysine req	Lysine supply in PF strategy 1
	Parity	ME MJ/d	Feed. kg/d	30 d	114 d	Average g/kg feed	L feed	Reduction PF vs CF,%
	1	31.0	2.4	3.63	6.23	4.01	67	17%
	2	34.0	2.6	3.20	5.80	3.62	78	24%

Precision Livestock Farming '17  585

 Precision Livestock Farming '17

 Precision Livestock Farming '17

 Precision Livestock Farming '17

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme, grant agreement No 633531. The data used for the simulations were issued from a project conducted by the Centre de Développement du Porc du Québec inc (CDPQ) in collaboration with INRA, within the AgriInnovation Program from Agriculture and Agri-food Canada.

This project has received funding from the European Union's Horizon 2020 research and innovation programme, grant agreement No 633531. The data 584  Precision Livestock Farming '17 used for the simulations were issued from a project conducted within the AgriInnovation Program from Agriculture and Agri-food Canada.


Ludovic Brossard
email: ludovic.brossard@inra.fr

I Taoussi

J Van Milgen

J Y Dourmad

Selection of methods to analyse body weight and feed intake data used as inputs for nutritional models and precision feeding in pigs

Keywords: pig, precision feeding, body weight, feed intake, real-time analysis

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

As other livestock production systems, pig production is continuously facing the challenge of sustainability. To contribute to the growing demand for animal protein, feed efficiency in pig production has to be optimised. This optimisation will contribute to reduce the environmental impact and to improve competitiveness of a production system where feed represents a major part of the production costs (typically 60 to 70%). The progress of technologies (sensors, automates) in precision livestock farming enables the development of potentially novel feeding techniques such as precision feeding. Precision feeding is based on the dynamic adjustment (if possible day by day) of the nutrient supply to the requirement, at a group or at an individual level [START_REF] Pomar | Applying precision feeding techniques in growing-finsihing pig operations[END_REF]. Recent studies have shown that precision feeding is a promising way to improve feed efficiency (e.g., [START_REF] Andretta | The impact of feeding growing-fisnihing pigs with daily tailored diets using precision feeding techniques on animal performance, nutrient utilization, and body and carcass composition[END_REF]. Compared to common phase-feeding programs applied to groups of pig, precision feeding allows a better consideration of the change in requirements during growth and of the variability among pigs of the same age, sex, and genotype [START_REF] Brossard | Modelling the variation in performance of a population of growing pig as affected by lysine supply and feeding strategy[END_REF].

In addition to the design of adapted feeders (e.g., [START_REF] Pomar | Applying precision feeding techniques in growing-finsihing pig operations[END_REF] for ad libitum feeding, [START_REF] Marcon | Precision feeding based on individual daily body weight of pigs with an automatic feeder developed to allow for restricting feed allowance[END_REF] for restricted feeding), and a precise evaluation of feed ingredients, precision feeding requires decision-support tools to manage data and apply nutritional models that calculate the optimal feed composition and allowance. Current models developed to simulate pig growth and to determine nutrient requirements (e.g. van Milgen et al., 2008) are difficult to be used in real-time decision-support tools as they require historical data on body weight (BW) and daily feed intake (DFI) data to characterize the animal. However, in precision feeding, nutrient requirements for individuals need to be determined in real time on the basis of their own growth and FI patterns. [START_REF] Hauschild | Development of sustainable precision farming systems for swine: Estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF] proposed a model for pigs fed ad libitum with empirical and mechanistic components, where the empirical model allowed the estimation of DFI and BW at time t+i from historical information measured for each individual animal up to time t. [START_REF] Hauschild | Development of sustainable precision farming systems for swine: Estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF] used the double exponential smoothing (DES) forecasting time series method on DFI and on weekly BW. However, with the technological development of sensors, BW can now also be obtained daily and individually (e.g., [START_REF] Marcon | Precision feeding based on individual daily body weight of pigs with an automatic feeder developed to allow for restricting feed allowance[END_REF].

The aim of this study was to test the most accurate forecasting method for BW and DFI for ad libitum feed allowance among three statistical methods: the DES method, the multivariate adaptive regression splines (MARS) method, and the k-nearest neighbours (kNN) method. These methods were tested in relation to data sampling frequency (i.e., daily or weekly measurements) and data availability.

Material and Methods

Dataset

A dataset on 119 pigs from Topigs was used for the analyses and calculations. During collection of data, animals had ad libitum access successively to two diets formulated to meet or exceed nutritional requirements (National Research Council, 1998), with a diet change at 65 kg BW: a growing diet (18.7% crude protein (CP), 10.51 MJ net energy (NE)/kg, 1.06 g standard ileal digestible lysine (SID Lys) / MJ NE on an as-fed basis), and a finishing diet (15.6% CP, 10.24 MJ NE/kg and 0.89 g SID Lys / MJ NE on an as-fed basis). During the test, DFI was recorded for each pig using automatic feeder systems. Animals were weighed daily by automatic weighing devices. Mean BW at the beginning and end of the data collection period were 34.0 ± 4.7 kg (75.9 ± 6.6 d of age) and 139.9 ± 9.8 kg (176.0 ± 9.2 d of age), respectively. Average DFI observed during the data collection was 2.32 ± 0.73 kg/d. On average, 99 and 101 data were available per animal for BW and DFI, respectively.

Calculation methods

To estimate the BW or DFI at a time t+i, using data up to time t, three forecasting methods were tested: the DES method, the MARS method, and the kNN method.

Double exponential smoothing (DES) method

In the panel forecasting time series methods, exponential smooting techniques are appropriate to reduce fluctuations from irregular observations in the studied time series [START_REF] Claycombe | Fundamentals of forecasting[END_REF], with a goal of short-time forecasting. As indicated by [START_REF] Hauschild | Development of sustainable precision farming systems for swine: Estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF], the DES is well adapted to the study of DFI and BW in pigs as they show long-run trends without seasonal components and because the method works with a limited number of observations (at least 3). At a given time t, the DES forecasting of the series to i days ahead is given by: 𝑋 ̂𝑡(𝑖) = 𝑎 𝑡 𝑖 + 𝑏 𝑡 where the coefficients 𝐚 𝐭 and 𝐛 𝐭 vary with time. In the DES method, the smoothed series 𝐒 𝟏 is smoothed a second time to obtain 𝐒 𝟐 with

𝑆 1 (𝑡) = 𝛼𝑦 𝑡 + (1 -𝛼) 𝑆 1 (𝑡 -1) 𝑎𝑛𝑑 𝑆 2 (𝑡) = 𝛼𝑆 1 (𝑡) + (1 -𝛼)𝑆 2 (𝑡 -1)
where α is the smoothing constant comprised between 0 and 1 and used to weigh differently past and recent observations. Values of α close to 1 give a greater weight to recent observations, while values of α close to 0 have a greater smoothing effect and are less responsive to recent changes.

Multivariate adaptive regression splines (MARS) method

The MARS method is an adaptation of techniques developed by [START_REF] Friedman | Multivariate Adaptive Regression Splines[END_REF] to resolve regression problems, with the aim to predict values of one or several continuous variables using a set of explicative variables. This nonparametric method has been largely used in data mining because it does not require a hypothesis on residues or relationships. The MARS method can be seen as an extension of linear regression to model automatically interactions and non-linearity. The method uses a database to establish functional relationships between explicative and predicted variables, even if the relationship between variables is not continuous and difficult to approximate with parametric models. In the MARS method, these relationships are approximate using linear-based functions such as:

𝐵(𝑥) = {
(𝑥 -𝑡) 𝑖𝑓 𝑥 ≥ 𝑡 𝑒𝑙𝑠𝑒 0 where t points are nodes connecting two regression segments. The general equation of the model for the variable t depending on the variable x is obtained by combining the linear-based functions estimated through the least squares method:

𝑦 = 𝑓(𝑡) = 𝛽 0 + ∑ 𝛽 𝑚 𝐵 𝑚 (𝑡) 𝑀 𝑚=0
where 𝐵 𝑚 are linear-based functions with associated coefficients β m .

K-nearest neighbours (kNN) method

The kNN method is an intuitive and non-parametric method used for classification and regression [START_REF] Altman | An introduction to kernel and nearest-neighbor nonparametric regression[END_REF]. The method is based on the determination of the k-nearest neighbours in a population of training values. In the regression application, the predicted value is the average of the k-nearest neighbours. To forecast a variable value at time t, the kNN method requires the use of historical data from the individual to be forecasted, and a learning database with values for previous days and time t for a population of individuals. Different parameters have to be determined for the calculation: the number k of neighbours, the type of distance to be calculated between individuals, the possible weighing of distance in calculation of the average (to assign weights to the contributions of the neighbours, so that the nearer neighbours contribute more to the average than the more distant ones), the number of previous days to be used to determine the nearest neighbours, and the weight to give to recent observations compared to past observations. For continuous variables such as BW or DFI, the Euclidian distance is commonly Precision Livestock Farming '17  577 used. Preliminary tests we made showed that a triweight kernel and k = 3 could be considered for calculations on BW and DFI.

Calculations

For daily measurements of FI and BW, the forecasting performance at day+1 of the DES and MARS methods were compared. For the DES method, values for the smoothing constant α ranging from 0.1 to 0.9 by a 0.1 increment were tested. As DFI can vary considerably from one day to another, we considered that the forecasting of DFI may be too sensitive to forecast these variations. Consequently, and similarly to BW, the cumulative FI was forecasted. The DFI can be determined from the cumulative FI.

For weekly available data, we considered that forecasting of DFI is required for the application of precision feeding. Consequently only forecasting of BW at day+7 was tested. The MARS method requires at least 8 historical data to perform forecasting whereas the DES method requires 3 data. Considering that this limits the usefulness of the MARS method, only the DES method was tested on BW for weekly data, with α ranging from 0.1 to 0.9. BW data at 7 day intervals were extracted from the original dataset to perform calculations.

The kNN method was used to illustrate the possibility to forecast BW or DFI at day+1 when only 1 or 2 historical data are available. This occurs for example at the beginning of data collection where the MARS or DES methods cannot be used. In the original dataset, data of 83 animals were used to create a learning database to forecast BW and DFI of the 36 other pigs.

All calculations were performed every day (or week) for each pig using the R software (version 3.3.2). The following functions and R packages were used: the earth function from the earth package [START_REF] Milborrow | Earth: multivariate adaptive regression spline models[END_REF] for the MARS method; the HoltWinters function from the stats package for the DES method; the kkNN function from the FNN [START_REF] Beygelzimer | FNN: Fast Nearest Neighbor Search Algorithms and Applications[END_REF] and kknn [START_REF] Schliep | InraPorc: A model and decision support tool for the nutrition of growing pigs[END_REF] packages.

Number of previous data used in calculation

The DES and MARS methods were tested on daily BW and cumulative DFI data using the 8, 13 or 20 historical (i.e., latest) data (the 8 data corresponds to the minimal number of data for the MARS method). Moreover, the 8 historical data refer approximatively to one week of data recording. The 13 and 20 historical data were chosen as they refer to a BW or DFI forecasting at 14 th or 21 st days, i.e. at the end of the 2 or 3 last weeks of data collection. Thus, the calculation started at day 9 and was performed on 8 historical data or integrated 578  Precision Livestock Farming '17 progressively up to 13 or 20 historical data depending on the targeted number of historical data. For the weekly BW, the calculation started at day 4 and a maximum of 8 historical data was used.

Missing data

For some pigs, data for BW or DFI were missing in the dataset. As the methods used cannot deal with missing values, missing data were created to obtain a specific completed datasets for each tested method. For data before day 4, BW and DFI were created by adding 0.75 kg and 27 g to previous BW and DFI, respectively. From day 4, data were created using the tested method with corresponding number of historical data.

Statistics

The residual mean square error of prediction (RMSEP) was calculated between forecasted values and measured value for each pig, excluding the first 8 days for tests on DFI and daily BW forecasting, and the first 3 days for weekly BW forescasting. The RMSEP were submitted to an analysis of variance (proc MIXED, SAS v9.4, Inst. Inc. Cary, NC). Least-square means were compared. For the daily BW or cumulative FI, the main effects were the forecasting method (with MARS and DES with α ranging from 0.1 to 0.9), the number of historical data (8, 13 or 20) and their interaction. For the weekly BW, the main effect was the α value (ranging from 0.1 to 0.9).

Results and Discussion

Daily forecasting of BW The RMSEP of daily BW decreased with increasing number of historical data used (Table 1). However, this decrease was significant only from 8 to 13 data for MARS method and α = 0.1 and 0.3 for DES method and from 13 to 20 data for α = 0.3 (method x number of historical data interaction, P < 0.001). The lowest RMSEP (1.1 kg) were obtained with the MARS method (13 or 20 historical data). This corresponds to 3% and 0.7% of BW at the beginning and the end of test period, respectively. For the DES method, the lowest RMSEP were obtained with α = 0.3 to 0.6. These α values give an intermediate weight between recent and less recent data. This means that most recent or oldest historical data should not be given too much weight to forecast BW. The results indicate that the MARS method used with 13 to 20 historical data could be preferred to the DES method for daily forecasting of BW, avoiding the choice of an α value. [START_REF] Quiniou | Assessment of the dynamic growth of the fattening pigs from body weight measured daily and automatically to elaborate precision feeding strategies In[END_REF] observed that the DES method with α = 0.6 to obtain the best RMSEP for data from pigs restrictively fed. [START_REF] Hauschild | Development of sustainable precision farming systems for swine: Estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF] used α = 0.1 to forecast BW (BW range = 25 to 105 kg, without reporting the RMSEP), giving a high weight to less recent data. This indicates that such comparisons could be influenced in part by the BW range or by the feeding level that can affect BW.

Weekly forecasting of BW The RMSEP of weekly BW were lowest with α = 0.5 to 0.6 (2.14 kg; Table 2) and increased for lower or higher α values. The highest RMSEP were obtained with α = 0.1 and 0.2. As for daily BW forecasting, giving a quite balanced weight to recent and less recent data allowed to better forecast BW, especially compared to privileging oldest data. The RMSEP values observed for weekly BW were 1 to more than 2 kg higher than for daily BW forecasting. The BW varies from day to day because of growth but also due to eating, defecating and urinating patterns. Consequently the precision of forecasting is affected when using weekly data where the difference between two successive points can be sensitive to the conditions of BW measurement.

Table 1: RMSEP (kg) of daily BW forecasting using the double exponential smoothing (α value ranging from 0.1 to 0.9) or MARS methods and 8, 13, or 20 historical data 1 . 1. Least-square means. The main effects of method, number of historical data and their interaction were significant at P < 0.001 (residual standard deviation of the model = 0.31 kg). Within a row, values followed by common letters are not significantly different for the method effect (P < 0.05). Daily forecasting of cumulated FI With the MARS method, the RMSEP of cumulated FI significantly increased with an increasing number of historical data used (Table 3). With the DES method, the number of historical data did not influence the RMSEP (P > 0.05), except for α = 0.2 where the RMSEP was significantly lower for 8 historical data compared to 13 or 20 historical data. The lowest RMSEP values (0.49 kg) were obtained with the DES method and α = 0.6 to 0.8, even though the difference between the MARS and DES methods was not significant for 8 historical data. These RMSEP values are lower than those obtained for daily BW, despite the fact that the daily increase in FI is higher than in BW (ADFI = 2.3 kg/d vs average daily gain = 1.06 kg/d) and the variation in DFI can be important from a day to another. This potential less smooth evolution of cumulative FI could also explain the increase in RMSEP with increasing number of historical data observed with the MARS method that uses combination of linear regressions. For the cumulative FI, the best forecasting was obtained by giving a higher weight to more recent data (α > 0.6). In contrast, [START_REF] Hauschild | Development of sustainable precision farming systems for swine: Estimating real-time individual amino acid requirements in growing-finishing pigs[END_REF] used α = 0.1 to forecast DFI (rather than cumulated FI). The present results indicate that the DES method used with 8 to 20 historical data could be preferred to the MARS method for daily forecasting of the cumulative FI, with α = 0.6 to 0.8, taking advantage of using a larger number of available historical data (more than 8) with a better RMSEP. 1. Least-square means. The main effects of method, number of historical data and their interaction were significant at P < 0.001 (residual standard deviation of the model = 0.17 kg). Within a row, values followed by common letters are not significantly different for the method effect (P < 0.05).

The kNN method for obtaining initial data We forecasted the BW and DFI with the kNN method at day 2 and 3 of data collection, on the basis of the first day or first 2 days of available data. The RMSEP of BW and DFI were 0.86 kg (± 0.73 kg) (n = 35 due to an outlier) and 0.33 kg (± 0.16 kg) (n = 36), respectively. This RMSEP was lower than for the two other methods although obtained with fewer data. The RMSEP of BW was higher and more variable than that for DFI, probably due to differences in the absolute value between BW and DFI at this stage (38 kg vs 1.27 kg/d, respectively). This method can be useful for forecasting when other methods cannot be used. However it requires a database on BW and DFI obtained in similar pigs reared in similar conditions. It is likely to be more sensitive to dayto-day variation than the DES or MARS methods that can smooth variation by using a larger number of historical data for the same pig.

Conclusions

The results of this study indicated that the MARS method used with 13 to 20 historical data is to be preferred to the DES method for daily forecasting of BW. Conversely the DES method is preferred to the MARS method for daily forecasting of cumulative FI, and can be used to forecast weekly BW. The kNN method can be useful to forecast BW or DFI at the start of data collection, when the two other methods cannot be used. The results about DES and MARS methods have to be confirmed on larger datasets and on different rearing conditions (e.g., with restricted feeding). However, the methods can be integrated for an efficient forecasting of BW and FI in a decision-support tool for precision feeding.

Table 2 :

 2 RMSEP (kg) of weekly BW forecasting using the double exponential smoothing method on 8 historical data with α values ranging from 0.1 to 0.9 1 .

					α value
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
	RMSEP (kg) 4.18 a 2.79 b 2.43 c 2.23 cd 2.14 d	2.14 d 2.19 d 2.28 cde 2.45 cde
	1. Least-square means. Effect of α value significant at P < 0.001 (residual
	standard deviation = 0.90 kg). Values followed by common letters are not
	significantly different (P < 0.05).

Table 3 :

 3 RMSEP (kg) of cumulated FI forecasting using the double exponential smoothing (α value ranging from 0.1 to 0.9) or MARS methods and 8, 13, or 20 historical data 1 .

	Nb.					Method
	of	MARS			Double exponential smoothing (α value)
	data		0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
	8	0.52 a	1.06 b 0.67 c 0.60 d 0.54 ae 0.50 ae 0.48 a 0.48 a	0.49 a 0.51 a
	13	0.60 a	1.04 b 0.73 c 0.61 a 0.54 d 0.50 de 0.49 e 0.49 e	0.49 e 0.52 de
	20	0.70 a	1.02 b 0.75 c 0.61 e 0.54 f	0.51 fg 0.49 g 0.49 g 0.50 fg 0.52 fg

 Precision Livestock Farming '17

Precision Livestock Farming '17  575

 Precision Livestock Farming '17

Precision Livestock Farming '17  579

Precision Livestock Farming '17  581

 Precision Livestock Farming '17

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant N° 633531. The authors gratefully acknowledge N. Quiniou and M. Marcon (IFIP) for methodological support.

This study was performed in the Feed-a-Gene Project funded by the European Union's H2020 Programme (grant agreement no 633531).