Michel Pierre
email: michel.pierre@bretagne.ens-cachan.fr

Julien Vovelle
email: vovelle@math.univ-lyon1.fr

On the strong convergence of the gradient in nonlinear parabolic equations

We consider the Cauchy-Dirichlet Problem for a nonlinear parabolic equation with L 1 data. We show how the concept of kinetic formulation for conservation laws [LPT94] can be be used to give a new proof of the existence of renormalized solutions. To illustrate this approach, we also extend the method to the case where the equation involves an additional gradient term.

We consider the question of existence of solution to the nonlinear parabolic problem

u t -div(a(∇u)) = f in Ω × (0, T ), (1a) u = u 0 on Ω × {0}, (1b)
u = 0 on Σ, ( 1c
)
where Ω is a bounded subset of R N , N ≥ 1, T is positive and Σ = ∂Ω × (0, T ). Let p > 1 be given. In (1), the operator -div(a(∇u)) is assumed to be a Leray-Lions operator of exponent p (for example the p-Laplacian):

Assumption 1 The function a ∈ C(R × R n , R N ) satisfies: there exists α > 0, β ∈ C(R + , R + ) such that a(X) • X ≥ α|X| p , (2a)

|a(X)| ≤ β|X| p-1 , (2b)
(a(X) -a(Y )) • (X -Y ) > 0, (2c)
for all distinct X, Y ∈ R N , where X • Y is the canonical scalar product of two vectors of R N and |X| the associated euclidean norm of X.

The framework is L 1 :

Assumption 2 The data u 0 , f are L 1 functions on Ω and Ω × (0, T ) respectively.

Remark 1 The flux a may depend on x and u. More general problems also may be considered, with additional first-order terms div(Φ(u)), div(g) in Equation (1a), as in [START_REF] Blanchard | Existence and uniqueness of a renormalized solution for a fairly general class of nonlinear parabolic problems[END_REF] for example.

Introduction

The existence of solution (precisely, of renormalized solution, see Definition 1 below) to Problem (1) or quite more general problems has already been proved: we refer in particular to the paper by Blanchard, Murat, Redwane [START_REF] Blanchard | Existence and uniqueness of a renormalized solution for a fairly general class of nonlinear parabolic problems[END_REF]. Our purpose here is to give a new proof of this fact. The cornerstone in the proof of existence of solution (by means of a process of approximation) of such a nonlinear parabolic Problem as (1) is the proof of the strong convergence of the gradient. We give a new method (inspired from the kinetic formulation of conservation laws developed by Perthame and coauthors [START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF][START_REF] Perthame | Kinetic formulation of conservation laws[END_REF][START_REF] Chen | Well-posedness for non-isotropic degenerate parabolichyperbolic equations[END_REF]) to prove this result.

Let us briefly summarize how and in which context the question of strong convergence of the gradient occurs. First, as soon as the problem under consideration involves a nonlinear function of the gradient, for example (under the hypotheses above with p = 2), the nonlinear elliptic Problem -div(a(∇u)) = g in Ω, u = 0 on ∂Ω, for g ∈ L 2 (Ω). Indeed, in order to prove existence of a solution (in H 1 0 (Ω)), it is usual to prove existence by approximation (e.g. by Galerkin approximation), i.e. for a set of data g n converging to g. Then weak convergence in H 1 0 of (a subsequence of) u n , the solution with datum g n , although easily obtained by uniform estimate on u n H 1 (Ω) , is not enough to pass to the limit in the equation since a is nonlinear: one has to prove the strong convergence of the gradient ∇u n . This is done by use of monotonicity methods. We refer to [START_REF] Minty | on a "monotonicity" method for the solution of nonlinear equations in Banach spaces[END_REF][START_REF] Browder | Variational boundary value problems for quasi-linear elliptic equations of arbitrary order[END_REF][START_REF] Leray | Quelques résultats de Višik sur les problèmes elliptiques nonlinéaires par les méthodes de Minty-Browder[END_REF], and [START_REF] Evans | Partial differential equations[END_REF] for a brief explanation of the technique.

Nonlinear expressions of the gradient also occur after renormalization of an elliptic or parabolic equation. Actually, they occur even if the original equation is linear. Nevertheless, renormalization for elliptic or parabolic equation has been introduced to deal with nonlinear equations with data of low regularity, so that the renormalized equation involves (at least) two nonlinear expressions of the gradient (see, e.g. Eq. (4) below). In any case, it will be necessary to prove the strong convergence of the (truncates of) the gradient in order to get existence of a solution by approximation.

We give a new proof of the strong convergence of the gradient by use of an equation on the characteristic function on the level sets of the unknown, similar to the kinetic formulation for conservation laws introduced in [START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF] (see also [START_REF] Perthame | Kinetic formulation of conservation laws[END_REF] and [START_REF] Chen | Well-posedness for non-isotropic degenerate parabolichyperbolic equations[END_REF] concerning the kinetic formulation of secondorder conservation laws). We intend to use it to study certain systems of reaction-diffusion equations (a forthcoming paper).

Let us conclude this introduction by a few words about the concept of renormalized solutions. Introduced by DiPerna and Lions for the study of ordinary differential equations and Boltzmann Equation [START_REF]Ordinary differential equations, transport theory and Sobolev spaces[END_REF][START_REF] Diperna | On the Cauchy problem for Boltzmann equations: global existence and weak stability[END_REF], it has been extended to nonlinear elliptic equations in [START_REF] Boccardo | Existence and regularity of renormalized solutions for some elliptic problems involving derivatives of nonlinear terms[END_REF] (in parallel with the (equivalent) notion of entropy solution [BBG + 95]) and has been extended to nonlinear parabolic equations in [START_REF] Blanchard | Truncations and monotonicity methods for parabolic equations[END_REF][START_REF] Blanchard | Existence and uniqueness of a renormalized solution for a fairly general class of nonlinear parabolic problems[END_REF][START_REF] Lions | Mathematical topics in fluid mechanics[END_REF] (in parallel with the (equivalent) notion of entropy solution [START_REF] Prignet | Existence and uniqueness of "entropy" solutions of parabolic problems with L 1 data[END_REF]). It has also been extended to first-order conservation laws [START_REF] Bénilan | Renormalized entropy solutions of scalar conservation laws[END_REF][START_REF] Porretta | L 1 solutions to first order hyperbolic equations in bounded domains[END_REF].

The problem of strong convergence of the gradient (hence the question of existence of solution) has initially be solved by the method of Minty-Browder and Leray-Lions [START_REF] Minty | on a "monotonicity" method for the solution of nonlinear equations in Banach spaces[END_REF][START_REF] Browder | Variational boundary value problems for quasi-linear elliptic equations of arbitrary order[END_REF][START_REF] Leray | Quelques résultats de Višik sur les problèmes elliptiques nonlinéaires par les méthodes de Minty-Browder[END_REF], then extended to the case of nonlinear elliptic (then parabolic) equations with less and less regular data by several methods, see, e.g. [BG92a, BM92, BGM93, BGDM93, DMMOP97, BDGO99, DMMOP99, BMR01, BP05]. Notice that this list of references to some works in the field of renormalized solutions for elliptic and parabolic equations is far from being complete.

The paper is organized as follows : in Section 2.1, we introduce the notion of renormalized solution and state the equivalent formulation by the (so-called) level-set P.D.E. In Section 2.2, we analyze this formulation and explain how it can be relaxed (although still characterizing renormalized solutions), see Theorem 2 and Lemma 1. In Section 2.3, we apply our tools to prove the convergence of an approximation to Problem (1) and thus existence of a renormalized solution to (1) (of course, we focus on the strong convergence of the gradient). In Section 3, we give the proofs of various results, which are reported at the end of the paper to let the main arguments of Section 2 stand out. Eventually, in Section 4, we extend the method to prove the existence of a renormalized solution to the Cauchy-Dirichlet Problem for a nonlinear parabolic equation with a term with natural growth.

Notations :

We set Q T := Ω×(-1, T ) and U T := Q T ×R. Any measurable function v : Ω×(0, T ) → R m is implicitly extended to a measurable function Q T → R m still denoted by v, defined by v ≡ 0 on Ω×(-1, 0).

If ν is a Radon measure on U T , we denote by ν * be the push-forward of ν by projection on R ξ :

ν * (E) = ν(Q T × E), ∀E ∈ B(R),
where B(R) is the Borel σ-algebra of R. More generally, if E is a topological space, B(E) denotes the σ-algebra of the Borel subsets of E. If q ≥ 1 and V is an open subset of R q we denote by D(V ) the set of smooth (C ∞ ) functions on V compactly supported in V and we denote by D ′ (V ) the set of distributions on V .

2 Existence of a renormalized solution -strong convergence of the gradient 2.1 Renormalized solutions

Renormalized solutions

For k > 0, we let T k (u) be the truncate of a function u at level k:

T k (u) := min(u, k) if u ≥ 0, T k odd. Definition 1 A function u ∈ L ∞ (0, T ; L 1 (Ω)
) is said to be a renormalized solution of the problem (1) if

(Regularity of the truncates)

T k (u) ∈ L p (0, T ; W 1,p 0 (Ω)), ∀k > 0.

(3)

2. (Renormalized equation) For every function S ∈ W 2,∞ (R) with S(0) = 0 such that S ′ has compact support, the equation

S(u) t -div(S ′ (u)a(∇u)) = S(u 0 ) ⊗ δ t=0 + S ′ (u)f -S ′′ (u)a(∇u) • ∇u (4)
is satisfied in the sense of distributions in Q T .

(Recovering at infinity)

lim k→+∞ QT ∩{k<|u|<k+1} a(∇u) • ∇udxdt = 0. (5)
2.1.2 Level-set P.D.E.

For α ∈ R, ξ ∈ R, we set χ α (ξ) = 1 0<ξ<α -1 α<ξ<0 . This is the "equilibrium function" in the kinetic formulation of conservation laws [START_REF] Lions | A kinetic formulation of multidimensional scalar conservation laws and related equations[END_REF]. Let u ∈ L ∞ (0, T ; L 1 (Ω)) satisfy (3). Then we define the (vector-valued) distribution a(∇u)δ u=ξ on U T by its restriction to each space D K (U T ) N (the set of smooth vector-valued functions with support in the compact subset K of U T ) as

a(∇u)δ u=ξ , α = Q a(∇T k (u)) • α(x, t, T k (u))dxdt, (6)
where

α ∈ D K (U T ) N , K ⊂ Q T × [-k, k].
Similarly, we define the distribution a(∇u) • ∇u δ u=ξ on U T by

a(∇u) • ∇u δ u=ξ , α = Q a(∇T k (u)) • ∇T k (u)α(x, t, T k (u))dxdt, (7)
for all α ∈ D K (U T ). By (3) and assumption (2b), we have

| a(∇u)δ u=ξ , α | ≤ a(∇T k (u)) L p ′ (QT ) α L p (K) ≤ βC(K) T k (u) L p (0,T ;W 1,p 0 (Ω)) α L ∞ (K) ,
and

| a(∇u) • ∇u δ u=ξ , α | ≤ β T k (u) L p (0,T ;W 1,p 0 (Ω)) α L ∞ (K)
. This shows that the right-hand sides of (6) and (7) are distributions on U T of order 0. To prove that (6) and (7) makes sense, we must also show that their respective right-hand sides do not depend on the choice of k: suppose k < k ′ for example, with

K ⊂ Q T × [-k, k], then α(x, t, T k ′ (u)) = 0 for |u| ≤ k only, in which case T k (u) = T k ′ (u).
With this definitions at hand, we can give the "level-set" formulation of Definition 1.

Theorem 1 A function u ∈ L ∞ (0, T ; L 1 (Ω))
→ χ u(x,t) (ξ), denoted by χ u , is solution in D ′ (U T ) of the equation ∂ t χ u -div(a(∇u)δ u=ξ ) = χ u0 ⊗ δ t=0 + f δ u=ξ + ∂ ξ µ, ( 8
)
where µ is defined by µ := a(∇u) • ∇u δ u=ξ , (9)

(Recovering at infinity)

lim k→+∞ QT ∩{k<|u|<k+1} a(∇u) • ∇udxdt = 0. ( 10
)
Proof of Theorem 1 : see Section 3.1.

Relaxation of the definition of renormalized solution -analysis of µ

Analysis of µ

Since µ ≥ 0, µ is represented by a nonnegative Radon measure on U T . We study the properties of the push-forward µ * of µ: µ

* (E) = µ(Q T × E), E ∈ B(R). Fact 1. For every h ∈ C c (R), R h(ξ)dµ * (ξ) = UT h(ξ)dµ(x, t, ξ). (11)
Proof: by definition of µ * , (11) is satisfied if h = 1 E is the characteristic function of a Borel set E ⊂ R, and therefore if h is a simple function. There exists a pointwise converging sequence of bounded simple functions with limit h with the same compact support as h: the Lebesgue dominated convergence Theorem gives the result.

Fact 2. For every h ∈ C c (R) with, say, supp(h) ⊂ [-k, k], R h(ξ)dµ * (ξ) = QT a(∇T k (u)) • ∇T k (u)h(u)dxdt. ( 12
) Proof: let (ϕ n ) be a nonnegative sequence of C c (Q T ) such that ϕ n ↑ 1 everywhere on Q T . By definition of µ, we have UT ϕ n (x, t)h(ξ)dµ(x, t, ξ) = QT a(∇T k (u)) • ∇T k (u)ϕ n (x, t)h(u)dxdt.
The Lebesgue dominated convergence Theorem then gives, at the limit

[n → +∞], UT h(ξ)dµ(x, t, ξ) = QT a(∇T k (u)) • ∇T k (u)h(u)dxdt.
We conclude by (11).

Fact 3. The measure µ * has no atom.

Proof:

Given k > 0, set v = T k (u). For ξ * ∈ (-k, k), let (h n ) be a sequence of C c (-k, k) converging
monotonically to 1 {ξ * } (take the h n to be tent functions for example). For every n, we have, by (12),

R h n (ξ)dµ * (ξ) = QT a(∇v) • ∇vh n (v)dxdt.
At the limit [n → +∞], we obtain, by the Lebesgue dominated convergence Theorem,

µ * ({ξ * }) = QT a(∇v) • ∇v1 {ξ * } (v)dxdt. ( 13
) For a.e. t, v(t) ∈ W 1,p (Ω). For such t's, we have ∇v(t) = 0 a.e. on {x ∈ Ω, v(x, t) = ξ * }. Indeed, we recall that, if w ∈ W 1,1 (Ω) and Z ⊂ R is a Borel negligible set, then the set {x ∈ Ω; w(x) ∈ Z, ∇w(x) = 0}
is negligible in Ω (the proof goes back to Stampacchia and can be found in [START_REF] Boccardo | Remarques sur l'homogénéisation de certains problèmes quasilinéaires[END_REF]). It follows therefore from (13) that µ * ({ξ * }) = 0.

Fact 4. For every l > k, R 1 (k,l) (ξ)dµ * (ξ) = QT ∩{k<u<l} a(∇u) • ∇udxdt. (14)
Proof: In the right hand-side of (14), u stands for

T m (u), m := max(|k|, |l|). Let (h n ) be a nonnegative sequence of C c (k, l) such that h n ↑ 1 { (k, l)}. For each n, we have by (12), R h n (ξ)dµ * (ξ) = QT a(∇u) • ∇uh n (u)dxdt.
At the limit [n → +∞], the dominated convergence Theorem gives the result.

Fact 5. For ϕ ∈ C c (Q T ), ϕ ≥ 0, define µ ϕ (A) := A ϕ(x, t)dµ(x, t, ξ), ∀A Borel subset of U T .
The measure µ ϕ has the same properties as µ and its analysis follows the same lines. In particular, µ ϕ, * has no atoms and, for every k > 0,

µ ϕ, * ([-k, k]) = µ ϕ, * ((-k, k)) = QT a(∇T k (u)) • ∇T k (u)ϕ(x, t)dxdt. ( 15
)
Remark 2 Note that the proof of the above Facts depends only on the property (3) of the truncates T k (u).

Actually, we may even replace a(∇u) by any measurable σ :

Q T → R N , such that σ1 |u|<k ∈ L p ′ (Q T ) N
for all k > 0. This will be used in paragraph 2.3.3.

Relaxation of the definition of renormalized solution

According to the above Facts (paragraph 2.2.1), the condition (10) may be rewritten in terms of the push-forward µ * uniquely as lim

k→±∞ µ * ((k, k + 1)) = 0, (16)
where we recall that µ is defined by (9). This simplifies the statement of Theorem 1 somewhat. However, what really makes plainer the characterization of renormalized solutions is the fact that, to some extent, it is not necessary to specify µ. This characterization is as follows.

Theorem 2 Let u be a function of L ∞ (0, T ; L 1 (Ω)) which has the regularity of the truncates (3) and satisfies the condition at infinity (5). Then u is a renormalized solution of the problem (1) if, and only if, there exists a nonnegative Radon measure µ on U T satisfying ( 16) and such that

∂ t χ u -div(a(∇u)δ u=ξ ) = χ u0 ⊗ δ t=0 + f δ u=ξ + ∂ ξ µ, ( 17
)
in the sense of distributions on U T .

The proof of Theorem 2 consists in showing that µ = a(∇u) • ∇u δ u=ξ . It is therefore a result of structure of µ: under the hypotheses of Theorem 2 and (17), µ has to be the measure a(∇u) • ∇u δ u=ξ . Theorem 2 has the virtue to give a plain characterization of renormalized solutions to (1). However, to prove the convergence of a sequence of approximate solutions to (1) and the existence of solution, we will need a slight generalization of Theorem 2 contained in the following lemma.

Lemma 1 Let u be a function of L ∞ (0, T ; L 1 (Ω)) which has the regularity of the truncates (3). Let σ be a measurable function

Ω × (0, T ) → R N such that ∀k > 0, σ1 |u|<k ∈ L p ′ (Q T ) N .
Suppose that there exists a nonnegative Radon measure µ on U T such that

lim k→±∞ µ * ((k, k + 1)) = 0, ( 18
)
and such that the following equation is satisfied in D ′ (U T )

∂ t χ u -div(σδ u=ξ ) = χ u0 ⊗ δ t=0 + f δ u=ξ + ∂ ξ µ. ( 19
)
Suppose also that: either u ≥ 0 a.e. and supp(µ

) ⊂ Q T × [0, +∞), ( 20
)
or the distribution σ • ∇u δ u=ξ satisfies the (sided) condition at infinity

lim sup k→+∞ σ • ∇u δ u=ξ , ϕ ⊗ 1 (k,k+1) ≤ 0, ∀ϕ ∈ C(Q T ), ϕ ≥ 0. ( 21
)
Then µ = σ • ∇u δ u=ξ .
Proof of Lemma 1 : see Section 3.2.

In Lemma 1 the definition of the distribution σ •∇u δ u=ξ is comparable to the definition of the distribution a(∇u) • ∇u δ u=ξ by (7):

σ • ∇u δ u=ξ , α = Q (σ1 |u|<k ) • ∇T k (u)α(x, t, T k (u))dxdt, ( 22
) for all α ∈ D K (U T ), K compact subset of Q T × [-k, k].
Equation ( 19) appears naturally when one considers limits of renormalized solutions, in particular of solutions of approximate equations u n tdiv(a(∇u n )) = f n : see Section 2.3.2 and Section 2.3.3. In the situation of Lemma 1, once the equality µ = σ • ∇u δ u=ξ has been proved, and thanks to Remark 2, we deduce the following corollary.

Corollary 1 Under the hypotheses of Lemma 1, and given ϕ ∈ C c (Q T ), ϕ ≥ 0, the measure µ ϕ, * has no atom and

µ ϕ, * ([-k, k]) = µ ϕ, * ((-k, k)) = QT σ • ∇T k (u)ϕ(x, t)dxdt, (23)
for all k > 0.

2.3 Existence of a renormalized solution -Strong convergence of the gradient

Approximation

Let (u n 0 ) and (f n ) be some approximating sequences of, respectively, u 0 and f in, respectively, L 1 (Ω) and

L 1 (Ω × (0, T )) such that u n 0 ∈ L p (Ω), f n ∈ L p ′ (Ω × (0, T )).
For each n, the problem

u n t -div(a(∇u n )) = f n in Ω × (0, T ), (24a) u n = u n 0 on Ω × {0}, (24b)
u n = 0 on Σ, (24c)
has a unique solution u n in the space

v ∈ L p (0, T ; W 1,p 0 (Ω)); v t ∈ L p ′ (0, T ; W -1,p (Ω)) .
(We refer to [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non-linéaires[END_REF] for example). The function u n is a weak solution to (24), hence a renormalized solution to (24), and therefore satisfies the equation

∂ t χ u n -div(a(∇u n )δ u n =ξ ) = χ u n 0 ⊗ δ t=0 + f n δ u n =ξ + ∂ ξ µ n , (25)
where µ n is defined by

µ n := a(∇u n ) • ∇u n δ u n =ξ . (26)

Estimates and limit equation

There are bounds independent on n on u n in L ∞ (0, T ;

L 1 (Ω)), on a(∇T k (u n )) • ∇T k (u n ) in L 1 (Q T ), on ∇T k (u n ) in L p (Q T ), on T k (u n ) t in L p ′ (0, T ; W -1,p (Ω))+L 1 (Ω×(0, T )).
They are obtained by multiplying the equation by T k (u n ) (see, e.g., [START_REF] Blanchard | Existence and uniqueness of a renormalized solution for a fairly general class of nonlinear parabolic problems[END_REF]). Aubin-Simon's compactness Theorem shows that there exists a function u ∈ L ∞ (0, T ; L 1 (Ω)) such that, up to a subsequence, u n → u a.e. and in L 1 (Q T ). By weak compactness of L p and L p ′ , we can suppose that

T k (u n ) → T k (u) in L p (0, T ; W 1,p 0 (Ω))-weak and a(∇T k (u n )) → σ k in L p ′ (Q T )-weak. Consider some common subsequences, still denoted (u n ), such that a(∇T k (u n )) → σ k and a(∇T k+1 (u n )) → σ k+1 in L p ′ (Q T )-weak. Since ∇T k (u n ) = ∇(T k • T k+1 )(u n ) = T ′ k • T k+1 (u n )∇T k+1 (u n ), and since a(0) = 0, we have a(∇T k (u n )) = 1 |u n |<k a(∇T k+1 (u n )). ( 27
)
As u n → u a.e., hence 1 |u n |<k → 1 |u|<k a.e. while being bounded, both sides of (27) converge in L p ′ (Q T )weak and, at the limit, we obtain σ k = 1 |u|<k σ k+1 . This shows (using a diagonal process) that there exists an additional subsequence still denoted (u n ) and a measurable function

σ : Ω × (0, T ) → R n such that ∀k > 0, σ1 |u|<k ∈ L p ′ (Q T ), a(∇T k (u n )) → σ1 |u|<k in L p ′ (Q T ) weak.
The bound on a(∇T

k (u n )) • ∇T k (u n ) in L 1 (Q T )
gives a uniform bound on µ n (K) for each compact subset K of U T . We can therefore suppose that (µ n ) converges weakly to a Radon measure µ in U T . Note that we have then

µ * (E) ≤ lim inf n→+∞ µ n * (E), (28)
for each E ⊂ R open. Indeed,

µ * (E) = µ(Q T × E) ≤ lim inf n→+∞ µ n (Q T × E) = lim inf n→+∞ µ n * (E), since Q T × E is open in U T .
With these results of convergence at hand, we let n → +∞ in (25) to obtain the limit equation

∂ t χ u -div(σδ u=ξ ) = χ u0 ⊗ δ t=0 + f δ u n =ξ + ∂ ξ µ. (29)
By multiplying (24a) by the function

(T k+1 -T k )(u n ) (k > 0) we obtain the estimate QT ∩{k<|u n |<k+1} a(∇u n ) • ∇u n dxdt ≤ Ω∩{|u n 0 |>k} |u n 0 |dx + QT ∩{|u n |>k} |f n |dxdt.
This is

µ n * ((k, k + 1)) + µ n * ((-k -1, -k)) ≤ Ω∩{|u n 0 |>k} |u n 0 |dx + QT ∩{|u n |>k} |f n |dxdt.
Up to a subsequence (and as a consequence of the strong convergence in L 1 ), there exists some functions

u 0 , u, f in L 1 (Ω) and L 1 (Q T ) respectively such that |u n 0 | ≤ u 0 , |u n | ≤ u, |f n | ≤ f a.
e. This implies the uniform estimates

µ n * ((k, k + 1)) + µ n * ((-k -1, -k)) ≤ Ω∩{u0>k} u 0 dx + QT ∩{u>k} f dxdt, (30)
from which we deduce by (28):

µ * ((k, k + 1)) + µ * ((-k -1, -k)) ≤ Ω∩{u0>k} u 0 dx + QT ∩{u>k} f dxdt.
In particular, lim

k→±∞ µ * ((k, k + 1)) = 0. ( 31
)
In the nonnegative case, i.e. u 0 , u n 0 ≥ 0 a.e., f, f n ≥ 0 a.e., the approximate solutions are nonnegative, and therefore u ≥ 0 a.e. and µ is supported in Q T × [0, +∞): Hypothesis (20) in Lemma 1 is satisfied. Let us show that, independently on any sign condition, Hypothesis (21) is satisfied: let ϕ ∈ C(Q T ), ϕ ≥ 0. For k > 0, n, m ∈ N, and by monotonicity of a, we have

0 ≤ QT a(∇v n k ) -a(∇v m k ), ∇v n k -∇v m k ϕdxdt, v n k := (T k+1 -T k ) + (u n ), i.e. QT (a(∇u n ) • ∇u m + a(∇u m ) • ∇u n )1 (k,k+1) (u n )1 (k,k+1) (u m )ϕdxdt ≤ QT (a(∇u n ) • ∇u n 1 (k,k+1) (u n ) + a(∇u m ) • ∇u m 1 (k,k+1) (u m ))ϕdxdt.
Denoting by ε k the right hand-side of (30), we deduce

QT (a(∇u n ) • ∇u m + a(∇u m ) • ∇u n )1 (k,k+1) (u n )1 (k,k+1) (u m )ϕdxdt ≤ 2ε k ϕ ∞ . Taking the limit [n → +∞], then [m → +∞] and [k → ∞], we obtain lim sup k→+∞ σ • ∇u δ u=ξ , ϕ ⊗ 1 (k,k+1) ≤ 0, (32)
which shows that Hypothesis (21) is satisfied.

Strong convergence of the gradient

Taking into account Eq. ( 29) and the estimates (31)-(32), we are in position to apply Lemma 1, which gives µ = σ • ∇u δ u=ξ .

We want to examine the weak convergence of the push-forward µ n * to µ * . We fix a test-function ϕ ∈ C c (Q T ), ϕ ≥ 0. We use the notations of Section 2.2.1, in particular

µ ϕ (A) := A ϕ(x, t)dµ(x, t, ξ), ∀A ∈ B(U T ). Then, if ψ ∈ C c (R), we have R ψdµ n ϕ, * = UT ϕ ⊗ ψdµ n , where ϕ ⊗ ψ(x, t, ξ) = ϕ(x, t)ψ(ξ) ∈ C c (U T ), hence R ψdµ n ϕ, * → UT ϕ ⊗ ψdµ = R ψdµ ϕ, * ,
and we conclude that (µ n ϕ, * ) converges weakly to µ ϕ on R. Let k > 0. By (23) the µ ϕ, * -measure of the boundary of [-k, k] is zero and, by weak convergence, we obtain

µ n ϕ, * ([-k, k]) → µ ϕ, * ([-k, k]). (33)
This identity (33) is the central result in the proof of the strong convergence of the gradient. Indeed, by ( 15) and ( 23), (33) reads

QT a(∇T k (u n )) • ∇T k (u n )ϕdxdt → QT σ • ∇T k (u)ϕdxdt (34)
and from (34) follows the strong convergence of the gradient

∇T k (u n ) → ∇T k (u) a.e. ( 35
)
Although the argument is classical, we give the proof of the implication (34) ⇒ (35) in Section 3.4. By (35), we have in particular σ = a(∇u) a.e. on Q T : therefore u is solution to the level-set p.d.e. associated to Problem (1):

∂ t χ u -div(a(∇u)δ u=ξ ) = χ u0 ⊗ δ t=0 + f δ u=ξ + ∂ ξ (a(∇u) • ∇u δ u=ξ ).
By Theorem 1, (u n ) converges to u, which is a renormalized solution to Problem (1).

3 Missing proofs 3.1 Proof of Theorem 1 8) is equivalent to: for all θ ∈ D(R),

Since {ϕ ⊗ θ; ϕ ∈ D(Q T ), θ ∈ D(R)} is dense in D(U T ), (
∂ t χ u -div(a(∇u)δ u=ξ ), θ D ′ (R ξ ),D(R ξ ) = χ u0 ⊗ δ t=0 + f δ u=ξ + ∂ ξ µ, θ D ′ (R ξ ),D(R ξ ) in D ′ (Q T )
. By definition of µ, this is equivalent to: for all θ ∈ D(R),

∂ t R χ u θdξ -div(θ(u)a(∇u)) = R χ u0 θdξ ⊗ δ t=0 + θ(u)f -θ ′ (u)a(∇u) • ∇u (36) in D ′ (Q T ).
The correspondence between (4) and ( 8) is obtained by taking θ = S ′ in (36), by the identity

R χ u (ξ)S ′ (ξ)dξ = S(u),
satisfied for all S ∈ W 2,∞ (R) such that S(0) = 0, and by a standard argument of density.

Proof of Lemma 1

Set ν := σ • ∇u δ u=ξ (see ( 22) for the definition of ν). We have to check that µ, ϕ ⊗ ψ = ν, ϕ ⊗ ψ for all ϕ ∈ D(Q T ), ψ ∈ D(R). We first suppose that ψ = ∂ ξ θ with θ ∈ D(R), so that µ, ϕ ⊗ ψ = -∂ ξ µ, ϕ ⊗ θ . By (19), µ, ϕ ⊗ ψ = ν, ϕ ⊗ ψ is then equivalent to the following identity

- T 0 Ω u u0 θ(ξ)dξ ϕ t + T 0 Ω (σ • ∇ϕ)θ(u) - T 0 Ω f ϕθ(u) = - T 0 Ω (σ • ∇u)ϕθ ′ (u). ( 37
)
By use of the rule of derivation of a product of functions in W 1,p ∩ L ∞ , we obtain the equivalent, more compact form of (37):

- T 0 Ω u u0 θ(ξ)dξ ϕ t + T 0 Ω σ • ∇(ϕθ(u)) - T 0 Ω f ϕθ(u) = 0. ( 38
)
Eq. ( 38) can be formally deduced from the chain-rule formula and from the equation

0 = ∂ t u -div(σ) -u 0 ⊗ δ t=0 -f. ( 39
)
Let us also remark that, formally, the equation ( 39) can be deduced from Eq. ( 19) by integrating with respect to ξ ∈ R. Indeed, that µ(ξ) → 0 when ξ → ±∞ is, still at the formal level, a consequence of the condition µ * ((k, k + 1)) → 0 when k → ±∞. Therefore, we begin with the derivation of an approximate form of Eq. ( 39): fix k > 0, let (ρ n ) n be an approximation of the unit on R (ρ

n having compact support in [-1/n, 1/n]), set α k := ρ k * 1 [k,k+1]
, and define

r k = r k (u) = ∞ |u| α k , v k := R χ u (ξ)r k (ξ)dξ, v k 0 := R χ u0 (ξ)r k (ξ)dξ. We have v k ∈ L p (-1, T ; W 1,p 0 (Ω)) ∩ L ∞ (Q T ), v k 0 ∈ L ∞ (Ω)
and, for l > 0,

r k → 1 a.e., T l (v k ) → T l (u) in L p (-1, T ; W 1,p 0 (Ω)), v k 0 → u 0 a.e.,
when k tends to +∞. Test Eq. ( 19) against ϕ(t, x)r k (ξ) to obtain

- T 0 Ω (v k -v k 0 )ϕ t + T 0 Ω σ • ∇ϕr k - T 0 Ω f ϕr k = T 0 Ω R ϕα k dµ. ( 40
)
This is the approximate form of (39). Now we want to use a kind of chain-rule formula to obtain an approximation of (38). To this purpose, we first infer from (40) the inequality

QT ϕ t (v k -v k 0 ) - T 0 G k , ϕ dt ≤ ϕ L ∞ ε k , (41)
where

G k := -(div(σr k (u)) + f r k (u)) ∈ L p ′ (0, T ; W -1,p ′ (Ω)) + L 1 (Q) and ε k := µ * ((k -1, k + 2)) + µ * ((-k -2, -k + 1)) → 0 when k → +∞.
We then consider the following lemma.

Lemma 2 Let ε > 0, v ∈ L p (0, T ; W 1,p 0 (Ω)) ∩ L 1 (Q T ), v 0 ∈ L ∞ (Ω) and G ∈ L p ′ (0, T ; W -1,p ′ (Ω)) + L 1 (Q T ) satisfy T 0 Ω ϕ t (v -v 0 ) - T 0 G, ϕ dt ≤ ϕ L ∞ ε, ( 42
)
for all ϕ ∈ D(Q T ). Then, for all ϕ ∈ D(R N × (-1, T )), for all h ∈ W 1,∞ (R) such that (h(v)ϕ)(t) = 0 on ∂Ω, for a.e. t ∈ (0, T ), ( 43
)
we have

T 0 Ω ϕ t v v0 h(ξ)dξ - T 0 G, h(v)ϕ dt ≤ ϕ L ∞ h L ∞ ε. ( 44
)
The Dirichlet condition (43) makes sense since h(v)ϕ ∈ L p (0, T ; W 1,p (Ω)). The proof of Lemma 2 is given in the following section. We apply Lemma 2 to (41

), with ϕ ∈ D(Q T ), h(v) = θ(v) to deduce T 0 Ω v k v k 0 θ(ξ)dξ ϕ t - T 0 Ω σr k (u) • ∇ ϕθ(v k ) + T 0 Ω f r k (u)ϕθ(v k ) ≤ ϕ ⊗ θ L ∞ ε k .
By use of the Lebesgue dominated convergence Theorem, we obtain (38) at the limit k → +∞. Recall that ψ = ∂ ξ θ, so that we actually proved that ∂ ξ (µν) = 0. By a classical Lemma in the theory of distributions, this shows that µν is constant with respect to ξ, or, more precisely, that for every κ ∈ D(Q T ) the distribution on R defined by ψ → µν, κ ⊗ ψ is represented by a constant c κ . There remains to show that c κ = 0.

In the case of nonnegative solution, i.e. under Hypothesis (20), this is straightforward since both µ and ν vanish on Q T × (-∞, 0) ξ . In the general case, i.e. under Hypothesis (21), we show that ν actually satisfies the condition at infinity

lim k→+∞ ν, ϕ ⊗ 1 (k,k+1) = 0, ∀ϕ ∈ C c (Q T ). ( 45
)
Since c κ = µν, κ ⊗ 1 (k,k+1) for all k, this will gives the result by (18).

To prove (45), we first observe that it is sufficient to obtain (45) for regular test-functions ϕ in the multiplicative form

ϕ(x, t) = ϕ 1 (t)ϕ 2 (x), ϕ 1 ∈ C 1 c (-1, T ), ϕ 2 ∈ C 1 c (Ω), ϕ i ≥ 0.
We then apply Lemma 2 to Eq. (40

) with ϕ(x, t) = ϕ 1 (t) ϕ 2 ∞ , h(v) = (T l+1 -T l )(v), l > 0 (observe that h ∈ W 1,∞ (R)
, and h(0) = 0 so that (43) is satisfied) and let k → +∞ to obtain Eq. (38) as above

with ϕ = ϕ 1 (t) ϕ 2 ∞ , i.e. T 0 Ω (σ • ∇u)ϕ 1 (t) ϕ 2 ∞ 1 (l,l+1) (u) = T 0 Ω u u0 (T l+1 -T l )(ξ)dξ ϕ ′ 1 (t) ϕ 2 ∞ + T 0 Ω f ϕ 1 (t) ϕ 2 ∞ (T l+1 -T l )(u)dξ.
Relabel l by k and take the limit

[k → +∞] to obtain lim k→+∞ ν, ϕ 1 ϕ 2 ∞ ⊗ 1 (k,k+1) = 0. ( 46
) Since -ϕ + ϕ 1 ϕ 2 ∞ ∈ C(Q T )
is nonnegative, we also have, by (21):

lim sup k→+∞ ν, (-ϕ + ϕ 1 ϕ 2 ∞ ) ⊗ 1 (k,k+1) ≤ 0.
This, combined with (46), shows that lim inf k→+∞ ν, ϕ ⊗ 1 (k,k+1) ≥ 0. Using (21) again, we obtain (45).

Proof of Lemma 2

It is a variation on the proof of Lemma 4.3 in [CW99] (Lemma 4.3 of [CW99] corresponds to the case ε = 0).

Step 1. Suppose that v 0 additionally satisfies v 0 ∈ W 1,p 0 (Ω). For t < 0, set v(t) = v 0 . Also first suppose h is non-increasing and ϕ nonnegative or h is non-decreasing and ϕ non-positive. We have

-ϕ L ∞ ε ≤ T 0 Ω ϕ t (v -v 0 ) - T 0 G, ϕ dt ≤ ϕ L ∞ ε (47)
for all ϕ ∈ D(Q T ) and thus, by regularity of v, G, for all ϕ ∈ L p (-1, T ;

W 1,p 0 (Ω)) ∩ L ∞ (Q T ) with ϕ t ∈ L p ′ (Q T ).
To use the function h(v) as a test-function in (47), we have first to regularize its dependence on t: for fixed ϕ ∈ D + (Q T ) and for η > 0 small enough (such that supp(ϕ)

⊂ Ω × (-1, T -2η]), we set ζ := ϕh(v), ζ η : (x, t) → 1 η t t-η ζ(x, s)ds.
In (47), this gives

T 0 G, ζ η dt ≤ ϕ L ∞ h L ∞ ε + T 0 Ω (ϕ η ) t (v -v 0 ) = ϕ L ∞ h L ∞ ε + T 0 Ω 1 η (ζ(x, t) -ζ(x, t -η))(v -v 0 )(x, t)dxdt = ϕ L ∞ h L ∞ ε + R Ω 1 η (v(x, t) -v(x, t + η))ζ(x, t)dxdt = ϕ L ∞ h L ∞ ε + R Ω 1 η (v(t) -v(t + η))h(v(t))ϕ(t)dxdt.
Since h is non-increasing and ϕ nonnegative or h is non-decreasing and ϕ non-positive, we have the inequality

(v(t) -v(t + η))h(v(t))ϕ(t) ≤ v(t+η) v(t) h(r)drϕ(t), t < T, hence T 0 G, ζ η dt ≤ ϕ L ∞ h L ∞ ε + R Ω ϕ(t) 1 η v(t+η) v(t) h(r)dr = ϕ L ∞ h L ∞ ε + R Ω 1 η (ϕ(t) -ϕ(t -η)) v(t) v0
h(r)dr.

At the limit η → 0, a first inequality is obtained

T 0 G, h(v)ϕ dt ≤ ϕ L ∞ h L ∞ ε + T 0 Ω ϕ t v(t) v0 h(r)dr.

By use of ζ

η : (x, t) → 1 η t+η t
ζ(x, s)ds as a test-function, we derive in a similar way the second inequality

T 0 G, h(v)ϕ dt ≥ -ϕ L ∞ h L ∞ ε + T 0 Ω ϕ t v(t) v0 h(r)dr,
which gives (44). In case h is non-decreasing and ϕ nonnegative or h is non-increasing and ϕ nonpositive, proceed similarly (just exchanging the order of the different time-regularizations) to prove (44), then decompose h as the sum of two monotone functions and ϕ as the sum of two signed functions to deduce the result in the general case.

Step 2. In the general case where

v 0 ∈ L ∞ (Q), regularize v 0 by v n 0 , v n 0 ∈ W 1,p 0 (Ω), v 0 -v n 0 L 1 (Ω) ≤ 1/n. Observe that, from (42), we deduce T 0 Ω ϕ t (v -v n 0 ) - T 0 G, ϕ dt ≤ ϕ L ∞ (ε + 1/n). (48)
Apply

Step 1. to get

T 0 Ω ϕ t v v n 0 h(ξ)dξ - T 0 G, h(v)ϕ dt ≤ ϕ L ∞ h L ∞ (ε + 1/n),
then pass to the limit [n → +∞] to achieve the proof of Lemma 2.

Proof of the strong convergence of the gradient

We start from (34) and prove the strong convergence of the gradient by the arguments of Minty, Browder and Leray, Lions [START_REF] Browder | Variational boundary value problems for quasi-linear elliptic equations of arbitrary order[END_REF][START_REF] Minty | on a "monotonicity" method for the solution of nonlinear equations in Banach spaces[END_REF][START_REF] Leray | Quelques résultats de Višik sur les problèmes elliptiques nonlinéaires par les méthodes de Minty-Browder[END_REF]. Let ϕ ∈ C c (Ω × (0, T )), ϕ ≥ 0 be given. Consider the sum

QT (a(∇T k (u n )) -a(∇T k (u))) • (∇T k (u n ) -∇T k (u))ϕdxdt. ( 49
)
We develop the product in this last term. The result (34) yields precisely the convergence of the term

QT a(∇T k (u n )) • ∇T k (u n )dxdt.
The other terms, which are linear with respect to ∇T k (u n ) or a(∇T k (u n )), converge by weak convergence. At the limit n → +∞ in (49), we obtain

lim n→+∞ QT (a(∇T k (u n )) -a(∇T k (u))) • (∇T k (u n ) -∇T k (u))ϕdxdt = 0. Since F n := (a(∇T k (u n )) -a(∇T k (u))) • (∇T k (u n ) -∇T k (u)
)ϕ is nonnegative (by monotonicity of a), this shows that F n → 0 in L 1 (Q T ). A subsequence of (F n ) (still denoted (F n )) therefore converges to 0 on a set A of full measure in Q T . Let (x, t) ∈ A and let q be an adherence value of (∇T

k (u n )) in R N .
Without loss of generality, we can suppose that ϕ(x, t) > 0. The vector q has finite-valued components as a consequence of the growth of a(∇T

k (u n )) • ∇T k (u n ), which gives (α|∇T k (u n )(x, t)| p -C|∇T k (u n )(x, t)|)ϕ(x, t) ≤ F n (x, t) → 0. At the limit [n → +∞] in F n (x, t) → 0, we thus obtain (a(q) -a(∇T k (u)(x, t))) • (q -∇T k (u)(x, t))ϕ(x, t) = 0.
By strict monotonicity of the flux a, q = ∇T k (u)(x, t). The sequence (∇T k (u n )(x, t)) has only one possible adherence value and is therefore convergent: ∇T k (u n ) → ∇T k (u) a.e. on Q T . Together with the uniform bound on ∇T k (u n ) in L p (Q T ), this shows the strong convergence of ∇T k (u n ) to ∇T k (u) in any L r (Q T ), r < p. Similarly, a(∇T k (u n )) converges to a(∇T k (u)) a.e. and in L r (Q T ), r < p ′ . In particular, σ = a(∇u) a.e.

To conclude, notice that we can recover the strong convergence ∇T k (u n ) → ∇T k (u) in L p loc (Q T ). Let K be compact subset of Q T . By the weak convergence of (a(∇T k (u n ))) and (∇T k (u n )) to a(∇T k (u)) and ∇T k (u) respectively, and by the convergence

(a(∇T k (u n )) -a(∇T k (u))) • (∇T k (u n ) -∇T k (u)) → 0 in L 1 (K), (a(∇T k (u n )) • ∇T k (u n )) converges to a(∇T k (u)) • ∇T k (u)) in L 1 (K)-weak. By Dunford-Pettis' Theorem, the family {a(∇T k (u n )) • ∇T k (u n )} is equi-integrable on K. By hypothesis (2a), |∇T k (u n )| p is dominated by α -1 a(∇T k (u n )) • ∇T k (u n ) and, therefore, {|∇T k (u n )| p } is also equi-integrable on K. Since ∇T k (u n ) → ∇T k (u) a.
(u n )| p → |∇T k (u)| p in L 1 (K). Besides, the weak convergence ∇T k (u n ) → ∇T k (u) in L p (K)-weak and the convergence ∇T k (u n ) → ∇T k (u) a.e. implies that lim n→+∞ ( ∇T k (u n ) L p (K) -∇T k (u n ) -∇T k (u) L p (K) ) = ∇T k (u) L p (K)
(this is a refinement of Fatou's Lemma by Brezis and Lieb [START_REF] Brézis | A relation between pointwise convergence of functions and convergence of functionals[END_REF]). By convergence of the L p (K)-norms, we conclude that lim

n→+∞ ∇T k (u n ) -∇T k (u) L p (K) = 0.
4 Parabolic equation with a term with natural growth

In this section, we briefly indicate how to adapt the arguments and proofs given above to solve the question of the strong convergence of the gradient (and, therefore, prove the existence of a renormalized solution) in the approximation by regularization and truncation of the following problem:

u t -div(a(∇u)) + γ(u)|∇u| p = f in Ω × (0, T ), (50a) u = u 0 on Ω × {0}, (50b)
u = 0 on Σ, (50c)
We keep the same assumptions on a and on the data: assumptions 1 and 2. The function γ ∈ C(R) is supposed to satisfies the sign condition

uγ(u) ≥ 0, ∀u ∈ R. (51)
This sign condition ensures a priori estimates for the additional term γ(u)|∇u| p , with a bound in L 1 (Q T ). More generally, we may consider a term γ(u)|∇u| r with a power r ∈ [1, p], instead of the term γ(u)|∇u| p .

Numerous works have been devoted to the study of Problem (50) (or to its elliptic version). Let us cite in particular [BMP83, BMP89, BG92b, BGM93, Por00, SdL03] and references therein.

In case p = 2, a = Id, there is a change of variables that transforms the equation in a classical Heat Equation:

v t -∆v = g, v = u 0 e -ξ 0 γ dξ, g = f e -u 0 γ .
It is this change of variables that we will adapt to the nonlinear case by use of the kinetic formulation (or level-set PDE).

A renormalized solution to (50) is defined as follows.

Definition 2 A function u ∈ L ∞ (0, T ; L 1 (Ω)) is a renormalized solution to (50) if T k (u) ∈ L p (0, T ; W 1,p 0 (Ω)), ∀k > 0,
where

R := γ(ξ)e -Γ+(ξ) {(α -1 1 ξ>0 + β -1 1 ξ<0 )µ n -|∇u n | p δ u n =ξ } (observe that the function ξ → γ(ξ)(α -1 1 ξ>0 + β -1 1 ξ<0 ) is continuous since γ(0) = 0). Since µ n = a(∇u n ) • ∇u n δ u n =ξ
, the hypotheses (2a) and (2b) on the flux a ensure that R ≥ 0 and, therefore, that

∂ t e -Γ+(ξ) χ u n -div(e -Γ+(ξ) a(∇u n )δ u n =ξ ) ≥ e -Γ+(ξ) (χ u n 0 ⊗ δ t=0 + f n δ u n =ξ ) + ∂ ξ (e -Γ+(ξ) µ n ). (56)
Similarly, we define

Γ -(ξ) =        1 β ξ 0 γ if ξ > 0, - 1 α 0 ξ γ if ξ < 0
and show the inequality

∂ t e -Γ-(ξ) χ u n -div(e -Γ-(ξ) a(∇u n )δ u n =ξ ) ≤ e -Γ-(ξ) (χ u n 0 ⊗ δ t=0 + f n δ u n =ξ ) + ∂ ξ (e -Γ-(ξ) µ n ). ( 57
)
It is then possible to pass to the limit [n → +∞] in ( 56) and (57) to obtain

∂ t e -Γ+(ξ) χ u -div(e -Γ+(ξ) σδ u=ξ ) ≥ e -Γ+(ξ) (χ u0 ⊗ δ t=0 + f δ u=ξ ) + ∂ ξ (e -Γ+(ξ) µ), (58)
∂ t e -Γ-(ξ) χ u -div(e -Γ-(ξ) σδ u=ξ ) ≤ e -Γ-(ξ) (χ u0 ⊗ δ t=0 + f δ u=ξ ) + ∂ ξ (e -Γ-(ξ) µ). (59)
What information do we extract from (58) and (59)? At a formal level, we can do the following computations: sum each inequality with respect to ξ ∈ R and use the condition at infinity (55) to obtain the (formal) weak equations

∂ t R e -Γ+(ξ) χ u dξ -div(e -Γ+(u) σ) ≥ R e -Γ+(ξ) χ u0 dξ ⊗ δ t=0 + e -Γ+(u) f, (60)
∂ t R e -Γ-(ξ) χ u dξ -div(e -Γ-(u) σ) ≤ R e -Γ l (ξ) χ u0 dξ ⊗ δ t=0 + e -Γ-(u) f . (61)
Multiply the first inequality by e Γ+(ξ)-Γ-(ξ) δ u=ξ and the second inequality by e -Γ+(ξ)+Γ-(ξ) δ u=ξ to obtain (still after formal computations)

∂ t e -Γ-(ξ) χ u -div(e -Γ-(ξ) a(∇u)δ u=ξ ) ≥ e -Γ-(ξ) (χ u0 ⊗ δ t=0 + f δ u=ξ ) -e -Γ-(ξ) σ • ∇δ u=ξ , ∂ t e -Γ+(ξ) χ u -div(e -Γ+(ξ) σδ u=ξ ) ≤ e -Γ+(ξ) (χ u0 ⊗ δ t=0 + f δ u=ξ ) -e -Γ+(ξ) σ • ∇δ u=ξ .
At last, use the identity e -Γ±(ξ) σ • ∇δ u=ξ = -∂ ξ (e -Γ±(ξ) ν), where

ν := σ • ∇u δ u=ξ ,
(this is also a very formal identity) to obtain ∂ t e -Γ-(ξ) χ udiv(e -Γ-(ξ) a(∇u)δ u=ξ ) ≥ e -Γ-(ξ) (χ u0 ⊗ δ t=0 + f δ u=ξ ) + ∂ ξ (e -Γ-(ξ) ν),

∂ t e -Γ+(ξ) χ udiv(e -Γ+(ξ) σδ u=ξ ) ≤ e -Γ+(ξ) (χ u0 ⊗ δ t=0 + f δ u=ξ ) + ∂ ξ (e -Γ+(ξ) ν). in D ′ (R). Consider the first of these inequalities. Since µ ϕ = ν ϕ = 0 on (-∞, 0), we have e -Γ+(ξ) µ ϕ ≤ e -Γ+(ξ) ν ϕ in D ′ (R). Similarly, using the second inequality, we obtain e -Γ-(ξ) µ ϕ ≥ e -Γ-(ξ) ν ϕ and conclude that µ ϕ = ν ϕ . This being true for every ϕ ∈ D + (Q T ), we have the desired result µ = ν.

Step 4. Strong convergence of the gradient. The identity µ = ν is the key point in the proof of the strong convergence of the gradient. Once this has been proved, we proceed as in Section 2.3.3. We prove in particular that ∇T k (u n ) → ∇T k (u) in L p loc (Q T ), and this allows to pass to the limit in Eq. ( 54) to obtain ∂ t χ udiv(a(∇u)δ u=ξ ) + γ(ξ)|∇u| p δ u=ξ = χ u0 ⊗ δ t=0 + f δ u=ξ + ∂ ξ (a(∇u) • ∇u δ u=ξ ), i.e. the fact that u is a renormalized solution.

Step 5. Rigorous proof of (64). This is a variation on the proof of Lemma 1 given in Section 3. We have v k ∈ L p (-1, T ; W 1,p 0 (Ω)) ∩ L ∞ (Q T ), v k 0 ∈ L ∞ (Ω) and T l (v k ) → T l (v) in L p (-1, T ; W 1,p 0 (Ω)) (l > 0), v k 0 → v 0 , r k → 1 a.e. when k tends to +∞. Test Eq. (58) against ϕ(t, x)r k (ξ) (with ϕ ∈ D + (Q T )), to obtain

- T 0 Ω (v k -v k 0 )ϕ t + T 0 Ω σ • ∇ϕr k - T 0 Ω f ϕr k ≥ T 0 Ω R ϕe -Γ+(ξ) α k dµ.
We deduce the inequality

QT ϕ t (v k -v k 0 ) - T 0 G k , ϕ dt ≤ ϕ L ∞ ε k ,
where G k := -(div(σr k (u)) + f rk (u)) ∈ L p ′ (0, T ; W -1,p ′ (Ω)) + L 1 (Q) and ε k := µ * ((k -1, k + 2)) → 0 when k → +∞. The analogue of Lemma 2 then shows that, for every h ∈ W 1,∞ (R), v k satisfies the following inequality:

QT ϕ t v k v k 0 h(ζ)dζ - T 0 G k , ϕh(v k ) dt ≤ ϕ L ∞ h L ∞ ε k .
Taking h with compact support, we obtain at the limit k → +∞ the inequality to obtain the weak form of (62). Similarly, we prove (63). As explained in Step 3., these two inequalities combined with (58) and ( 59) imply (64).

  e., and by Vitali's Theorem, this implies |∇T k

  the starting point (58)-(59) to deduce the inequalities∂ ξ (e -Γ+(ξ) µ) ≤ ∂ ξ (e -Γ+(ξ) ν), ∂ ξ (e -Γ-(ξ) ν) ≤ ∂ ξ (e -Γ-(ξ) µ).(64)Assume for the moment that (64) is satisfied in D ′ (U T ). A test-function ϕ ∈ D + (Q T ) being fixed, we consider the distributions on R defined byµ ϕ : ψ → µ, ϕ ⊗ ψ , ν ϕ : ψ → ν, ϕ ⊗ ψ .They satisfy the inequalities∂ ξ (e -Γ+(ξ) µ ϕ ) ≤ ∂ ξ (e -Γ+(ξ) ν ϕ ), ∂ ξ (e -Γ-(ξ) ν ϕ ) ≤ ∂ ξ (e -Γ-(ξ) µ ϕ )

  2. Let us explain the main arguments. Introduce α k := ρ k * 1 [k,k+1] , and definer k = r k (u) = ∞ |u| α k , v k := R e -Γ+(ξ) χ u (ξ)r k (ξ)dξ, v k 0 := R e -Γ+(ξ) χ u0 (ξ)r k (ξ)dξ.Set also rk = e -Γ+(u) r k and v := R e -Γ+(ξ) χ u (ξ)dξ, v 0 := R e -Γ+(ξ) χ u0 (ξ)dξ.

  )dζ -QT e -Γ+(u) σ • ∇(ϕh(v))dt + QT e -Γ+(u) f ϕh(v) ≤ 0. (65)We then fix θ ∈ D(R) and apply (65) withh(ζ) := e -(Γ--Γ+)(φ -1 (ζ)) θ(φ -1 (ζ)), φ(ξ) -(ξ) θ(ξ)dξ, h(v) = e -(Γ--Γ+)(u) θ(u),

and, for every function S ∈ W 2,∞ (R) such that S ′ has compact support and S(0) = 0, S(u) tdiv(S ′ (u)a(∇u)) + S ′ (u)γ(u)|∇u| p = S(u 0 ) ⊗ δ t=0 + S ′ (u)f -S ′′ (u)a(∇u) • ∇u and lim k→+∞ QT ∩{k<u<k+1} a(∇u) • ∇udxdt = 0.

We can also use directly the level-set PDE and define a renormalized solution to (50) as a function u ∈ L ∞ (0, T ; L 1 (Ω)) having the regularity of the truncates T k (u) ∈ L p (0, T ; W 1,p 0 (Ω)), ∀k > 0, which satisfies the equation:

where µ := a(∇u) • ∇u δ u=ξ , and satisfies the condition at infinity lim

We now explain how to prove the existence of a renormalized solution to Problem (50). For the sake of simplicity, we will suppose that the solution has a sign: we assume

Step 1. Approximation. Let (u n 0 ) and (f n ) be some nonnegative approximating sequences of, respectively, u 0 and f in, respectively, L 1 (Ω) and

has a unique solution u n in the space of functions v ∈ L p (0, T ; W 1,p 0 (Ω)) with v t ∈ L p ′ (0, T ; W -1,p (Ω)). The function u n is a weak solution to (24), hence a renormalized solution and therefore satisfies the equation

where µ n is defined by

Step 2. Estimates. As in Section 2.3.2, we show that, up to a subsequence,

|u|<k and µ n → µ weakly. We also prove, by the same technique as in Section 2.3.2, the conditions at infinity

Since u n ≥ 0 a.e., we also have u ≥ 0 a.e. and µ is supported in Q T × [0, +∞).

Step 3. Limit of the equation. To pass to the limit of Eq. ( 54), there is a difficulty in the fact that the term γ(ξ)|∇u n | p δ u n =ξ is uniformly bounded in L 1 and that no stronger a priori bound is available. We define

The function Γ + is continuous, not C 1 , on R, but a step of regularization shows that we have ∂ t e -Γ+(ξ) χ u ndiv(e -Γ+(ξ) a(∇u n )δ u n =ξ ) = e -Γ+(ξ) (χ u n 0 ⊗ δ t=0 + f n δ u n =ξ ) + ∂ ξ (e -Γ+(ξ) µ n ) + R


J Sadki

L Alfandari
email: alfandari@essec.fr

A Plateau
email: aplateau@cnam.fr

A Nagih

Space minimization in agricultural production planning by column generation

Keywords: Crop rotations, agricultural production planning, column generation, covering integer programming, dynamic programming

We deal in this paper with an agricultural production planning problem where crops must be scheduled on land plots so as to satisfy crop demands every period of time and to minimize the overall surface of land used for cultivation. This problem can be formulated as a covering integer program with a huge number of variables. A resolution scheme based on column generation is thus proposed, where the resulting pricing problem is efficiently solved by dynamic programming. The numerical experiments show that the method is all the more so efficient and robust as the planning horizon is long and plot sizes are small.

Introduction

This paper deals with a class of production problems in agriculture, that of planning crops over a given time horizon and given land plots so that the total production of each crop every season meets the needs of the farmer. The difficulty of these problems lies in the decreasing yields of cultivated lands over time, which enforce to plan fallow periods over the time horizon and to alternate crops in the best possible way, which is not a trivial issue. Whereas typical crop-rotation problems focus on building rotations that maximize a profit or yield function as in [START_REF] Detlefsen | Modelling optimal crop sequences using network flows[END_REF][START_REF] Klein Haneveld | Crop succession requirements in agricultural production planning[END_REF][START_REF] El-Nazer | The choice of crop-rotation: a modelling approach and case study[END_REF][START_REF] Santos | Crop rotation scheduling with adjacency constraints[END_REF], the objective of the specific problem studied in this paper is to minimize space needed to meet the demand. This problem was originally presented in [START_REF] Lemalade | A MIP flow model for crop rotation planning in a sustainable development context[END_REF][START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF] in a Madagascan context where the minimization of cultivated space helped to contribute in the long-term to the sustainable development of the primary forest around the farms. The model and software developed in [START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF] were part of a project of struggle against deforestation sponsored by the French-Madagascan Cooperation. Whereas the previous paper mainly focused on the formulation of the application and agricultural conclusions, this paper focuses on dedicated solving methods for a generic variant, where a single crop can be cultivated on a land plot every season when this plot is not left fallow. The proposed solving method is a column-generation method based on a Dantzig-Wolfe decomposition of the problem. A column generation approach was also applied in [START_REF] Santos | Crop rotation scheduling with adjacency constraints[END_REF], where the objective is to maximize the total time during which the land is cultivated (space is fixed contrary to our problem). Let us remark that the space-minimization problem studied in this paper has the particularity that the objective value is a multiple of the number of selected columns in the master program, which is also the number of plots used for production. Another particularity of this paper is that the planning is studied over a much longer horizon, up to 30 years, because the aim, inherited from the Madagascan project, was not only to satisfy the farmer constraints but also to measure the long term impact on the primary forest. This critical parameter in the computational tests motivated the use of a column generation approach which is known to be a powerful tool for solving large size covering problems.

The paper is organized as follows. Section 2 introduces the problem statement and data as well as the covering integer program formulation of the generic crop rotation problem. Section 3 describes the pricing problem and the dynamic programming method dedicated to its resolution. Section 4 presents computational experiments for various lengths of the time horizon and various sizes of the plots, comparing the column generation approach to direct solving using CPLEX solver. The numerical results show that the column generation algorithm is all the more so efficient as the time horizon is long. Section 5 concludes the paper.

The covering master problem for column generation

The Minimum-Space Crop Rotation Planning (MSCRP) problem is that of constructing crop rotations minimizing required space area for covering seasonal crop demands. We consider in this paper the case where all plots have the same unitary surface area and yield characteristics, which happens when the planting is done in the same area, so the objective amounts to minimizing the number of plots used. MSCRP can be formulated by an arc-flow model, as done by Alfandari et al. in [START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF] for a variant of the problem presented in this paper. However, a natural alternative attempt to model it consists in considering the set of feasible rotations over the scheduling time horizon, and then selecting a minimum number of rotations that meets the seasonal crops demand. Note that there is an equivalence between the arc-flow model and the model presented in this paper using the classical Dantzig-Wolfe decomposition. Let us define the notations that will be used throughout the paper. We consider:

• T = {1, .., T } the set of time periods of the planning horizon • R the set of feasible rotations of a plot over the planning horizon • C the set of crops that can be cultivated over the time horizon • b jt , j ∈ C, t ∈ T the demand (in tons) of crop j on season t • ρ r jt , j ∈ C, t ∈ T , r ∈ R the unitary yield of crop j (in tons) provided by rotation r on season t, which can be 0 if crop j is not cultivated on that period t in rotation r.

• S the surface area of a plot and let θ r be a decision variable representing the number of plots on which rotation r ∈ R is planned. The Minimum-Space Crop Rotation Planning (MSCRP) problem can be formulated as a Covering Integer Program where the objective function minimizes the total surface area used, while the covering constraints ensure the crop demand satisfaction. It is given as follows for a given surface area S:

(M SCRP ) S              z(S) = min r∈R Sθ r s.t r∈R Sρ r jt θ r ≥ b jt ∀j ∈ C, t ∈ T θ r ∈ N ∀r ∈ R. (1)
As problem (M SCRP ) S is untractable in this form because of the huge number of variables representing all feasible rotations, we apply the column generation approach suggested by Gilmore and Gomory [START_REF] Gilmore | A Linear Programming Approach to the Cutting-Stock Problem[END_REF]. It consists of generating new rotations (i.e. columns of the constraint matrix) only when needed instead of enumerating all feasible ones, in order to solve the continuous relaxation of the problem.

We denote the linear relaxation of (M SCRP ) S as the Master Problem (M P ) S and we note z(S) its optimal value:

(M P ) S              z(S) = min r∈R Sθ r s.t r∈R Sρ r jt θ r ≥ b jt ∀j ∈ C, t ∈ T θ r ≥ 0 ∀r ∈ R (2)
The well-known iterative principle of a column generation approach can be summarized as follows. We start with a limited number of crop rotations and solve a master problem reduced to this subset of columns. Then, given the dual variables obtained by this LP solving, we check whether there is any new feasible crop rotation of negative reduced cost that could be added to the restricted master problem in order to improve the LP bound. If no such negative reduced cost column is found, then the current solution is optimal for (M P ) S , otherwise we add a subset of negative reduced-cost columns to the Master Problem and reiterate the process. We keep the constant factor S in the objective because, as shown in the following proposition, this formulation will allow to solve the same pricing problem whatever the value of S.

Proposition 1. The continuous relaxation value of the master problem (M P ) S is a constant independent of S.

Proof. Let µ jt be the dual variables associated with the covering constraints, and let (D -M P ) S denote the dual problem associated with (M P ) S .

(D -M P ) S              z D (S) = max j∈C t∈T b jt µ jt s.t j∈C t∈T Sρ r jt µ jt ≤ S ∀r ∈ R µ jt ≥ 0 ∀j ∈ C, t ∈ T (3)
Dividing by S in the (D -M P ) constraints, we obtain a dual problem that does not depend on S. As z D (S) = z(S), we get the claimed result.

Pricing problem formulation

We can either generate new columns in the column generation process or verify optimality of a current solution by solving a pricing problem producing the minimum reduced cost rotation among all rotations of the original set R. Because of identical yields and surface areas of plots, a single pricing problem is solved at every iteration of the column generation process. We first present its formulation and then its efficient solving by dynamic programming.

Let µ jt be the dual variables associated with the covering constraints of the master problem (M P ) S , and let c r be the reduced cost of rotation r ∈ R. The pricing problem of minimizing the reduced cost of variables θ r is given as follows:

min r∈R cr = S(1 - t∈T j∈C ρ r jt µ jt ) (4)
Let us now model precisely a valid crop rotation of a plot. A rotation is defined as a sequence of cultivation and fallow periods over the time horizon. Let P (j) be the set of crops that can precede crop j ∈ C. We denote by (l, a) the state of a cultivated plot at period t, meaning that l fallow periods preceded a cultivation periods until current cultivated period t, where 1 ≤ l ≤ l and 1 ≤ a ≤ ā (l will be called the fallow length and a will be called the cultivation age of the plot). l is the number of periods beyond which no yield improvement is achievable when going back to cultivation, and ā is the number of cultivation periods beyond which the plot necessarily returns fallow because yield has decreased in too large amounts. If the plot has been left fallow at current period for l periods, l is also called the state of the plot. Taking into account the state of the plot enables to better model yields through plot history, as explained in [START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF]. Indeed, we denote by w la ij the yield of a plot in state (l, a) cultivated with crop j in current period preceded by crop i ∈ P (j) (decreasing function of a and increasing function of l), and by w l1 f j the yield of a plot that returns cultivated with crop j after l fallow periods (the yield is equal to zero when a plot is left fallow). In order to represent sequences that will form a rotation, we define decision variables associated with a given period t ∈ {2, ..T }, state (l, a), 1 ≤ l ≤ l, 1 ≤ a ≤ ā, crops j ∈ C and i ∈ P (j):

x la ijt =
1 if plot is in state (l, a), cultivated with crops j and i in periods t and t -1 respectively 0 otherwise we distinguish then the cases when a fallow of length l (1 ≤ l ≤ l) precedes or follows a crop j ∈ C:

x l f jt =
1 if plot returns to cultivation with crop j ∈ C in period t ∈ T after l fallow periods 0 otherwise

x la jf t =    1 if plot is left fallow in period t (t ≥ 2)
and was cultivated with crop j ∈ C in state (l, a) in period t -1 0 otherwise and finally the case of a fallow-fallow sequence without exceeding the maximal length l for 1 ≤ l ≤ l -1:

x l f f t =
1 if, in period t, plot has been left fallow for l + 1 periods 0 otherwise and the particular case where the maximal fallow length is exceeded on a plot:

F l t =
1 if plot has been left fallow at least l + 1 periods at period t 0 otherwise

These last two series of variables differ from [START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF] because going back to cultivation after a long fallow period is not mandatory in this paper, and lead to the following flow conservation constraints to define feasible crop rotations.

1 = j∈C x l f j1 + F l 1 (5) i∈P (j) x la ijt = k:j∈P (k)
x l,a+1 jk,t+1 + x la jf,t+1 ∀t ∈ {2, . . . , T -1}, j ∈ C, a ∈ {2, . . . , a -1}, l ∈ {1, . . . , l} [START_REF] Detlefsen | Modelling optimal crop sequences using network flows[END_REF] i∈P (j)

x la ijt = x la jf,t+1 ∀t ∈ {2, . . . , T -1}, j ∈ C, l ∈ {1, .., l} (7)
x l f jt = k∈C:j∈P (k)

x l2 jk,t+1 + x l1 jf,t+1 ∀t ∈ {1, . . . , T -1}, l ∈ {1, . . . , l}, j ∈ C (8)
l∈{1,...,l} a∈{1,...,a} j∈C

x la jf t = x 1 f f,t+1 + i∈C x 1 f i,t+1 ∀t ∈ {2, . . . , T -1} (9)
x l-1 f f t = j∈C x l f j,t+1 + x l f f,t+1 ∀t ∈ {1, . . . , T -1}, l ∈ {2, . . . , l -1} (10)
F l t + x l-1 f f t = j∈C x l f j,t+1 + F l t+1 ∀t ∈ {1, . . . , T -1} (11)
Assuming that the plot is initially in a fallow state of length l, constraint (5) establishes possible successions during the first period of the planning horizon. Constraints (6) are the flow conservation constraints when two crops follow one another without exceeding maximal cultivation age ā. When the maximal cultivation age is reached, then a fallow period is required to reset the plot; this is given by constraints [START_REF] Dogliotti | ROTAT, a tool for systematically generating crop rotations[END_REF]. Constraints [START_REF] Durr | Tiling with bars under tomographic constraints[END_REF][START_REF] El-Nazer | The choice of crop-rotation: a modelling approach and case study[END_REF] are the flow conservation constraints when a fallow state precedes or follows a crop ensuring the possibility to move from a cultivation to a fallow period and vice versa. Finally, constraints (10-11) define the flow conservation when there is a succession of fallows resulting either from a fallow of maximal age or a fallow of lower age.

Solving the pricing problem by dynamic programming

As established before, the objective of the pricing problem is to determine a rotation with minimum reduced cost according to the dual values provided by LP solving of the current Master Problem. The yield ρ r jt of a given crop j ∈ C at period t ∈ T in rotation r ∈ R can be written as

ρ r jt = l l=1 ( i∈P (j) a a=2
w la ij x la ijt +w l1 f j x l f jt ), where variables x la ijt , x l f jt satisfy constraints [START_REF] Gilmore | A Linear Programming Approach to the Cutting-Stock Problem[END_REF][START_REF] Detlefsen | Modelling optimal crop sequences using network flows[END_REF][START_REF] Dogliotti | ROTAT, a tool for systematically generating crop rotations[END_REF][START_REF] Durr | Tiling with bars under tomographic constraints[END_REF][START_REF] El-Nazer | The choice of crop-rotation: a modelling approach and case study[END_REF](10)(11). Thereby, the objective of the pricing problem (4) aiming at minimizing the reduced cost of rotations is given as follows:

min S -

t∈T j∈C l l=1 Sµ jt (w l1 f j x l f jt + i∈P (j) a a=2 w la ij x la ijt )
which is equivalent to

min 1 - t∈T j∈C l l=1 µ jt (w l1 f j x l f jt + i∈P (j) a a=2 w la ij x la ijt ) (12)
The following proposition provides a dynamic programming method for solving the pricing subproblem of finding the minimum reduced cost column at each iteration of the column generation scheme. Proposition 2. Let π jlat and π f lt denote the values computed by the following recurrence formula for t = 1, .., T :

π jlat = π f l,t-1 -µ jt w l1 f j if a = 1 min i∈P (j) {π i,l,a-1,t-1 -µ jt w la ij } if 2 ≤ a ≤ ā π f lt =      min j∈C,l ≤ l,a≤ā, {π jl a,t-1 } if l = 1 π f,l-1,t-1 if 2 ≤ l < l min(π f,l-1,t-1 , π f,l,t-1 ) if l = l
with initial conditions:

π f l0 = 1, π f l0 = π jla0 = ∞ ∀l ∈ {1, .., l -1}, a ∈ {1, .., ā}, j ∈ C Then, π * T = min(min j,l,a π jlaT , min l π f lT )
will be the value of a minimum reduced cost column. The corresponding column can be found by backtracking from t = T to t = 0.

Proof. We show that the recurrence formula of Proposition 2 amount to solve a shortest path in an acyclic graph G where paths represent all feasible columns and the value of a path is the reduced cost of a column. The acyclic graph is a connected graph G = ∪ T t=1 G t , where G t = (V t-1 ∪ V t , A t ) is a bipartite graph associated with each pair of periods (t -1, t), such that:

V t = {(j, l, a, t) j∈C,l∈{1,.., l},a∈{1,..,ā} } ∪ {(f, l, t) l∈{1,.., l} } where • (j, l, a, t), j ∈ C, l ∈ {1, .., l}, a ∈ {1, .., ā} represent cultivation state nodes, meaning that plot is in state (j, l, a) at period t

• (f, l, t), l ∈ {1, .., l} represent fallow state nodes, meaning that plot has been left fallow for l periods.

and

A t = {(u, v) ∈ V t-1 × V t : state u can precede state v}
represents feasible state sequences between each pair of nodes from t -1 to t according to the flow conservation constraints [START_REF] Gilmore | A Linear Programming Approach to the Cutting-Stock Problem[END_REF][START_REF] Detlefsen | Modelling optimal crop sequences using network flows[END_REF][START_REF] Dogliotti | ROTAT, a tool for systematically generating crop rotations[END_REF][START_REF] Durr | Tiling with bars under tomographic constraints[END_REF][START_REF] El-Nazer | The choice of crop-rotation: a modelling approach and case study[END_REF](10)(11). Values λ t uv on arcs (u, v) ∈ A t are equal to the opposite of the yield associated to state v after state u, multiplied by the dual variable associated with state v, i.e.

λ t uv =              -µ jt w la ij if u = (i, l, a -1, t -1), v = (j, l, a, t), -µ jt w l1 f j if u = (f, l, t -1), v = (j, l, 1, t), 0 if u = (j, l, a, t -1), v = (f, 1, t), 0 if u = (f, l, t -1), v = (f, l + 1, t), 0 if u = (f, l, t -1), v = (f, l, t) (13)
Figure 1 illustrates an example of a bipartite state subgraph G t defined on C = {Rice, P otatos} with (l, a) = (2, 2).

By construction, there is a 1-1 correspondence between paths in graph G and feasible columns of the Master Program, i.e. each path in G corresponds to a crop rotation on the plot, and vice versa. Also by construction, the value of a path in G is the reduced cost of the corresponding column. As graph G is acyclic, it can be solved by the Bellman-Ford algorithm computing for t = 1, . . . , T

π v = min u:(u,v)∈A t (π u + λ t uv ) ∀v ∈ V t (14) 0 0 season t -1 season t (P otato, 2, 2, t) (Rice, 2, 2, t) (Rice, 1, 2, t) (P otato, 1, 2, t) (Rice, 2, 1, t) (P otato, 2, 1, t) (Rice, 1, 1, t) (P otato, 1, 1, t) (f, 2, t) (f, 1, t) (P otato, 2, 2, t -1) (Rice, 2, 2, t -1)
(Rice, 1, 2, t -1)

(P otato, 1, 2, t -1) (Rice, 2, 1, t -1) (P otato, 2, 1, t -1) (Rice, 1, 1, t -1) (P otato, 1, 1, t -1) (f, 2, t -1) (f, 1, t -1) -µ R ic e ,t w 2 1 f R ic e
µ P o ta to ,t w 1 2

P o ta to ,P o ta to -µ R ic e ,t w 2 2 R ic e ,R ic e Figure 1: Bipartite subgraph G t = (V t-1 ∪ V t , A t ) with ( l, ā) = (2, 2)
with initial conditions π f, l,0 = 1, π v = ∞ for v = (f, l, 0), and output min v∈V T π v . As v ∈ V t describes all possible states (j, l, a, t) and (f, l, t), the above Bellman-Ford recurrence formula are equivalent to the dynamic programming formula of Proposition 2. The backtracking process to find the optimal column is the same as backtracking from t = T to t = 0 in graph G to find the path associated with the optimal value.

The next section is devoted to computational results comparing the column generation scheme to a direct solving of the arc-flow model by a commercial solver.

Computational experiments

To evaluate the contribution of the column generation approach, we have modeled the explicit arc-flow model whose constraints on plots are constraints ( 5)-( 11) with our problem specifications. In the following sections, direct solving means that the solution of the explicit arc-flow model is computed by the MIP-solver of cplex. We have generated (M SCRP ) S instances in the same way as in [START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF].

Regarding the column generation strategies, we opted to initialize the first restricted master problem with a heuristic set of columns provided by a greedy algorithm and to add, at each iteration, a set of columns of negative reduced cost obtained from solving the pricing problem by the dynamic programming method described above. Finally, to get integer solutions, we use an enumeration heuristic based on a depth-first search, within a specified time limit, based on the last master problem. We call this resolution scheme the Integer Column Generation Heuristic (ICGH).

The comparative tests are performed for time horizons T = 20, 30, 40, and 60 seasons (resp. 10, 15, 20, and 30 years) that correspond to a long term strategic planning in a context of sustainable development. On the other hand, for each class of instances of a fixed time horizon we evaluate the impact of varying the area S which gives another point of view that is more relevant for farmers than the long-term planning aspect. All the computations reported in this section have been carried out on a DELL Latitude D610 Personal Computer with Pentium M-1.73 Ghz processor and 512 Mo RAM. The code has been written in C++ and the solver is Cplex 12.0. About agricultural data, we considered a set C of two crops such as each one can precede the other, and we fixed the maximal cultivation and fallow ages to a = 4, l = 5 respectively as in [START_REF] Alfandari | A MIP flow model for crop-rotation planning in a context of forest sustainable developement[END_REF].

The first observation that can be made concerns the continuous relaxation value. As expected the value of the column generation relaxation, which is a lagrangian relaxation, is, in most cases, tighter (very slightly, however) than the classical continuous relaxation of direct solving and is obtained in very short time as shown in Table 1 where average CPU indicates the average running time over 10 instances of different S for each class of a fixed T . In addition, as the problem size increases (eg. T = 60 seasons) direct solving by Cplex is even unable to find a feasible solution for 3 instances within the allowed time of 3600s CPU. 2 presents the results of the Integer Column Generation Heuristic (ICGH) (i.e column generation followed by the enumeration heuristic) compared to direct solving. We compare the two methods in terms of: number of used plots (which, multiplied by corresponding surface S equals the objective value z(S)) and computational time CP U . We also report the time to the best integer solution found by both methods. We can see that, for most instances, (ICGH) provides integer solutions of high quality. More precisely, over the 40 instances we have that:

• for 95% of instances integer solutions found by (ICGH) are either optimal, or lower or equal to the best integer bound reached by direct solving on its time limit • for only 5% of instances, they are worse than direct solving with a gap of 1 unit corresponding to only one supplementary plot Furthermore, direct solving failed to find integer solutions in 25% of the cases.

Concerning the running time of the two methods, direct solving is clearly more timeconsuming than column generation. However, in the integer column generation heuristic, the enumeration heuristic dedicated to get integer solutions is still time consuming (see Total CPU), but as reported in column CPU-2-best, in most cases, this time is spent to prove optimality of the obtained integer solutions that are often found at the beginning of the enumeration process. On the other hand, direct solving is even unable to find feasible solutions when the surface value decreases and the time horizon increases, while column generation had the same trend for any value of S but increases only with T with a very low computational time. Plots of figures 2 -3 express statistical performance in terms of total running time to the best integer of each resolution method in function of the surface area in figure 2 and in function of the running time in figure 3. In order to have an homogeneous comparison, this average is computed over all T for a given S for figure 2, and it is computed over all S for a given T for figure 3.

We clearly see in figure 2 that the behavior of direct solving depends on the value of S, whereas (ICGH) is insensitive to this parameter. Also, we observe in figure 3 the dramatic increase of the running time while increasing the size of the problem, when column generation increases in a monotonic way with a small slope.

In conclusion, these tests show that our method based on column generation outperforms direct solving on all the following performance criteria: robustness, solution quality and running time.

Conclusion

In this paper, we presented a heuristic resolution scheme based on column generation for Minimum-Space Crop Rotation Planning. The master problem is modeled as a covering integer program and the pricing problem is efficiently solved by dynamic programming. Numerical tests demonstrate that this resolution scheme can be used on problems with up to 60 seasons, as well as problems of small values of surface area where direct solving failed to find solutions in the time limit. The method is quite robust then for long-term strategic planning in the framework of an environmental evaluation, and for small division of land as this is often the case in developing countries.

Figure 2 :Figure 3 :

 23 Figure 2: Evolution of CPU-time with surface S = 0.1, 0.2, .., 10

Table 1 :

 1 Continuous resolution of Master problem vs. Arc-flow modelTable

			Continuous resolution of:
	T	Master problem by CG	Arc-flow model by Cplex
		Optimal value average CPU Optimal value average CPU
	20	18.57	1.88	18.54	82.91
	30	19.07	7.23	19.05	247.21
	40	19.23	18.46	19.22	563.86
	60	19.32	99.37	19.32	***

Table 2 :

 2 Comparison between direct solving with Cplex and integer column generation heuristic ICGH ** a feasible integer solution was not obtained within one hour • for 37.5% of instances they are strictly better than the best integer bounds obtained by direct solving

			Integer Column Generation Heuristic		Direct Solving
	T	S # plots CG CPU Total CPU CPU-2-best # plots	CPU CPU-2-best
	20 0.1	187	1.34	>3600	1.34	188	>3600	3500
		0.2	94	1.31	>3600	1.31	94	>3600	2813
		0.5	38	1.51	1.52	1.52	38	477,46	477,46
		0.7	28	1.56	>3600	1.56	28	>3600	619,64
		1	20	3.66	>3600	3.66	20	>3600	270.23
		2	11	2.44	334.44	2.44	11	>3600	97.11
		4	6	1.84	2.11	2.11	6	2648	40.31
		6	4	1.65	1.69	1.69	4	11	11.07
		8	3	1.88	4.25	4.25	3	94	94.04
		10	3	1.64	1.89	1.89	3	3.95	3,95
	30 0.1	192	6.14	>3600	6.14	**	>3600	>3600
		0.2	97	8.01	>3600	8.01	**	>3600	>3600
		0.5	40	7.11	>3600	7.11	40	>3600	2282,64
		0.7	28	7.48	7.58	7.58	28	951,57	951,57
		1	20	9.30	9.41	9.41	21	>3600	755,29
		2	11	6.78	6.94	6.94	11	>3600	359,63
		4	6	1.96	2.02	2.02	6	62.4	62,4
		6	4	7.46	7.5	7.5	4	31.78	1.78
		8	4	8.09	34.84	8.09	3	550.85	550.85
		10	4	10.06	78	10.06	3	569.9	569.9
	40 0.1	194	19.51	>3600	19.51	**	>3600	>3600
		0.2	97	16.07	60.43	60.43	**	>3600	>3600
		0.5	40	16.3	601	16.3	**	>3600	>3600
		0.7	28	22.23	46.43	46.43	28	3469,1	3469,1
		1	20	15.37	15.9	15.9	20	1479.3	1479.3
		2	11	13.3	216.8	216.8	12	>3600	954.42
		4	6	13.98	61	61	7	>3600	279.3
		6	4	24.74	31.16	31.16	4 1476.48	1476.48
		8	4	22.15	103.83	22.15	4	>3600	140.13
		10	4	21	2048	21	4	>3600	71.32
	60 0.1	195	102.5	>3600	102.5	**	>3600	>3600
		0.2	98	128.8	>3600	128.8	**	>3600	>3600
		0.5	40	108.8	>3600	118.8	**	>3600	>3600
		0.7	28	123.14	123.87	123.87	**	>3600	>3600
		1	20	57.69	70.83	70.83	**	>3600	>3600
		2	11	67.21	>3600	615.21	13	>3600	2243.66
		4	6	60.21	640	640	8	>3600	1348.34
		6	4	119.4	363.92	363.92	4	1458.6	1458.6
		8	4	102	>3600	102	4	>3600	183.49
		10	4	124	3400	124	4	>3600	190.47