Fatima-Zahra Tani

Alain Rapaport

Térence Bayen

Periodic controls for discriminating density dependent growth in the chemostat

We show the benefit of considering periodic dilution rates in the chemostat model for discriminating between a growth function which does not depend on the density of the micro-organisms population (such as the Monod law) and another one which does depend (such as the Contois law). This goal is achieved thanks to the measurement of the abiotic resource only. We then present a simple procedure for a robust discrimination between the two types of kinetics using a single experiment in three phases. Finally, the shape of the best periodic excitation is discussed and the method is illustrated on numerical simulations.

I. INTRODUCTION

Continuous culture has been invented simultaneously by Monod [START_REF] Monod | La technique de culture continue: Théorie et applications[END_REF] and Novick & Szilard [START_REF] Novick | Experiments with the chemostat on spontaneous mutations of bacteria[END_REF] in the fifties as a means to measure accurately growth rates of microorganisms. Typically, the so-called chemostat device consists in feeding a culture vessel of volume V with a nutritive solution of concentration s in at a constant flow rate Q, and then in extracting micro-organisms and substrate at the same rate Q from the culture vessel (so that the volume remains constant). The concentrations of micro-organisms and substrate, denoted respectively by b and s are then solutions of the following differential equations, that reflect the writing of the mass balance

ṡ = - 1 Y µ(s)b + D(s in -s), ḃ = µ(s)b -Db, (1)
where D := Q/V denotes the dilution rate, Y the conversion rate and µ(•) the specific growth rate. Each bacterial species is characterized by the positive constant Y and the nonnegative function s → µ(s) satisfies µ(0) = 0. The usual way to identify Y and µ(•) are as follows.

1. The parameter Y is obtained from the batch mode which consists in taking D = 0 (see for instance [START_REF] Corman | Bacterial growth measurement using an automated system: Mathematical modelling and analysis of growth kinetics[END_REF]). One can check that solutions of (1) with positive initial conditions (s 0 , b 0 ) converge asymptotically to the steady state (s , b ) := (0, b 0 +Y b 0 ). The parameter Y is thus determined by

Y = b -b 0 s 0 , s
which is usually considered as a robust method.

F.-Z. Tani and A. Rapaport are with MISTEA, Univ. Montpellier, INRA, Montpellier SupAgro, France alain.rapaport@inra.fr, fatima-zahra.tani@supagro.fr T. Bayen is with Avignon Université, Laboratoire de Mathématiques d'Avignon (EA 2151) F-84018 Avignon, France terence.bayen@univ-avignon.fr 2. There are two general approaches for estimating the growth rate function µ(•): either by adjusting the dynamical model to growth curves (in batch or in continuous culture) [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF], [START_REF]Dochain State and parameter estimation in chemical and biochemical processes : a survey[END_REF], either by a series of steady-states. The first method is often considered less accurate because it deals with transient dynamics and biomass sampling can perturb the culture. The second method measures steady-states only, for which cell density no longer change with time (historically this method owes its name to "chemostat"), and is considered more reliable [START_REF] Hoskisson | Continuous culture -making a comeback[END_REF]. The growth function µ(•) is then identified step by step, considering a series of experiments with different constant values of D (see for instance [START_REF] Ziv | The use of chemostats in microbial systems biology[END_REF]). One can easily check that positive equilibria (s , b ) of ( 1) are given by µ(s ) = D and b = Y (s ins ), when s is below s in . Therefore each value of D provides a point (s , D) of the graph of µ(•). This method implicitly assumes that there exists a unique solution s to the equation µ(s) = D and that the corresponding steady state is asymptotically stable. One can easily show (see for instance [START_REF] Harmand | The Chemostat: mathematical theory of microorganisms cultures[END_REF]) that this is fulfilled when µ(•) is monotone and D is chosen less than µ(s in ). For non-monotonic growth functions, other methods can be used (see for instance [START_REF] Sieber | A new method for the reconstruction of unknown non-monotonic growth function in the chemostat[END_REF]).

Most often, the growth functions that suit data are concave increasing, such as the Monod law [START_REF] Monod | La technique de culture continue: Théorie et applications[END_REF]:

µ M (s) := µ max s K s + s . (2)
However, the function µ, differently to the parameter Y , might depend on environmental growth conditions (temperature, pH,...). Parameters µ max , K s could then depend on experiments and measurements. It also happens that the growth is inhibited by large concentrations b of the biomass, reflecting a crowding effect of micro-organisms. This is particularly met in bioreactors of industrial wastewater plants. Then, one has to look for density dependent functions, i.e., of the form (s, b) → µ(s, b), increasing with respect to s and non-increasing with respect to b. The most popular one is the Contois law [START_REF] Contois | Kinetics of Bacterial Growth: Relationship between Population Density and Specific Growth Rate of Continuous Cultures[END_REF], given by the following expression

µ C (s, b) := µ max s s + kb . (3)
Notice that Contois law can be interpreted as a Monod expression but applied to the ratio s/b instead of s. However it does not generalize the Monod function for small density effect, but one can consider the generalized Contois law

µ CG (s, b) := µ max s K s + s + kb , (4)
to recover Monod law (taking k = 0) and Contois law (taking K s = 0) expressions. Here the parameter k > 0 measures the density effect. Recently, another expression of density dependent function has been proposed and experimentally validated [START_REF] Krichen | High biomass density promotes density-dependent microbial growth rate[END_REF]:

µ K (s, b) := µ 0 s b α , (5)
where α is a non-negative parameter.

Before setting up a model identification and control strategies, it is worth to mention that the choice between Monod law or density dependent extensions has important impacts on the predictions of the behavior of bacterial species (the greater the parameter k or α is, the greater the impacts are). First of all, steady state concentrations of substrate do not depend on the input concentration s in in the case of Monod law, while it does in the case of Contois law [START_REF] Harmand | The Chemostat: mathematical theory of microorganisms cultures[END_REF]. This is clearly an issue for the piloting of waste-water plants for which the input concentration s in fluctuates with time. Secondly, the Competitive Exclusion Principle [START_REF] Smith | The Theory of Chemostat -Dynamics of microbial competition[END_REF], [START_REF] Harmand | The Chemostat: mathematical theory of microorganisms cultures[END_REF] applies when the growth functions depend on s only: it asserts that only one species coexists in the long run in the chemostat model. On the contrary, a dependency on biomass concentration in the kinetics can allow coexistence of several species [START_REF] Lobry | Persistence in ecological models of competition for a single resource[END_REF], [START_REF] Leenheer | Crowding effects promote coexistence in the chemostat[END_REF], [START_REF] Lobry | Persistence in Ratio-Dependent Models of Consumer-Resource Dynamics[END_REF], which can impact the performances of bioreactors. Therefore, the discrimination between substrate and ratio dependency is a question of prime importance, in particular for the choice of sensors and control laws.

The classical identification method we recalled previously imposes many experiments to reconstruct accurately the graph of µ (as one experiment provides only one point of its graph) and can be then lengthy and costly, whereas control laws do not necessarily require a perfect knowledge of the function µ (when it depends on s only, see for instance [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF], [START_REF] Rapaport | Robust regulation of a class of partially observed nonlinear continuous bioreactors[END_REF], [START_REF] Antonelli | Set point regulation of an anaerobic digestion process with bounded output feedback[END_REF]). Notice that the identification of density dependent functions requires the additional measurement of the biomass, which is often more difficult and expensive to obtain than the abiotic component. Ratio-dependency has received a great attention in the ecological literature (see for instance [START_REF] Arditi | How species interact? Altering the Standard View on Trophic Ecology Oxord[END_REF]), and a conceptual experiment of several tanks in series has been designed to qualitatively discriminate ratiodependency in the context of prey-predators [START_REF] Arditi | Saiah Empirical Evidence of the Role of Heterogeneity in Ratio-Dependent Consumption[END_REF]. Inspired by this approach, a similar experiment has been investigated in the context of bioreactors [START_REF] Harmand | Density-dependent kinetics models for a simple description of complex phenomena in macroscopic massbalance modeling of bioreactors[END_REF]. However, it imposes to be able to retain perfectly the biomass from one reactor to the following one, which can be a practical issue.

The purpose of the present work is to investigate a robust methodology for discriminating between the two kinds of growth functions, on the basis of a simple experiment with a single tank and the single measurement of the substrate. Instead of steady state operations, as in the methods described above, we propose to control the experiment with periodic dilution rate D(•) instead of constant ones. Periodic operations of chemostat have already been investigated in the literature, but for optimizing performances when the growth function is known (see for instance [START_REF] Abulesz | Periodic Optimization of Microbial Growth Processes[END_REF], [START_REF] Abulesz | Periodic Operation of a Continuous Culture of Baker's Yeast[END_REF]).

The paper is organized as follows. In Sec. II, we give results about periodic solutions of the chemostat model with periodic controls. In particular, we show in Proposition 2 the role played by the convexity of µ(•) in the average value of s(•) over a period, compared to a steady state s . This result is the core of the procedure that is presented in Sec. III. In this section, we first define "weak" and "strong" density dependent effects, and we then give an experiment procedure based on bang-bang controls for discriminating between these two aforementioned cases with the single measurement of the variable s only. In Sec. IV, we justify the choice of bangbang controls as the optimal ones. Finally, Sec. V illustrates this methodology on numerical simulations showing also the robustness of the method when environmental conditions could suddenly change.

II. THE CHEMOSTAT MODEL UNDER PERIODIC CONTROL

In the sequel, we shall consider non-negative measurable control functions D : R + → R + such that: -the control D(•) is periodic of period T (to be chosen), -the control D(•) is persistently exciting ( [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF]) and satisfies

1 T T 0 D(t) dt = D > 0, (6)
where D > 0 is given. Considering the total mass concentra-

tion m = Y s + b, (1) implies that ṁ = D(t)(Y s in -m).
Then, any periodic solution of (1) has to fulfill m(t) = Y s in for any time t. Therefore, we shall consider the reduced dynamics on the invariant line Y s + b = Y s in , that is, the scalar dynamics

ṡ = F(s, D) := (s in -s)(D -ν(s)), (7)
on the invariant interval I := (0, s in ), where the function ν is defined as follows

ν(s) := µ(s) if µ = µ M , µ(s,Y (s in -s)) if µ = µ CG .
Notice that in any case, the function ν is increasing. Therefore, for D < ν(s in ), there exists a unique solution s of and(s , x ) is thus the unique asymptotically stable equilibrium of (1) for the constant control D = D. Let us define for convenience the function

ν(s ) = D, with x := Y (s in -s ) > 0,
λ (D) := ν -1 (D) if D < ν(s in ), s in otherwise.
We start by giving a result about the convergence of solutions of [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] to a unique T -periodic solution, using the usual argumentation for scalar dynamics based on the Poincaré mapping (see for instance [START_REF] Smith | The Theory of Chemostat -Dynamics of microbial competition[END_REF]).

Lemma 1: Suppose that D < ν(s in ). Then, for any Tperiodic control D(•) that satisfies [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF] and any initial condition s 0 ∈ I, the associated solution of ( 7) converges asymptotically to the unique T -periodic solution of (7) associated with D(•).

Proof: For a given T -periodic control D(•), consider the C 1 function φ : I → I defined by φ (s 0 ) := s(T, s 0 ) where s(•, s 0 ) denotes the solution of [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] for the initial condition s 0 at time 0. By differentiating φ , one has

φ (s 0 ) = ∂ s ∂ s 0 (T, s 0 ) = X(T ),
where X(•) is the fundamental solution of the variational equation ż = ∂ f ∂ s (s(t, s 0 ), D(t))z with z(0) = 1. One gets:

φ (s 0 ) = e -T 0 [D(t)-ν(s(t,s 0 ))+ν (s(t,s 0 ))(s in -s(t,s 0 ))] dt .
Clearly, φ > 0 and thus φ is increasing. For s 0 = 0, the solution cannot be identically equal to 0 because of condition [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF] and F(0, D) ≥ 0. Therefore one has φ (0) = s(0, T ) > 0.

For s 0 = s in , the solution of ( 7) is constant equal to s in . Then φ (s in ) = s in and

φ (s in ) = e -T 0 [D(t)-ν(s in )] dt .
Since D < ν(s in ), we deduce that φ (s in ) > 1. Consequently, the C 1 function s 0 → ψ(s 0 ) := φ (s 0 )s 0 satisfies ψ(0) > 0, ψ(s in ) = 0, and ψ (s in ) > 0, thus it has necessarily a zero in I which proves the existence of a fixed point s0 of φ as well as the existence of a T -periodic solution s(•) := s(•, s0 ) of [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF]. Moreover, at s 0 = s0 , one finds

ψ ( s0 ) = e -T 0 [D(t)-ν( s(t))+ν ( s(t))(s in -s(t)] dt -1.
Observe now that for a periodic solution s(•), one has

0 = ln s in -s(T ) s in -s0 = - T 0 [D(t) -ν( s(t))] dt, (8)
from which one obtains

ψ ( s0 ) = e -T 0 ν ( s(t))(s in -s(t)) dt -1 < 0.
Since at every point s0 ∈ I such that ψ(s 0 ) = 0, one has ψ ( s0 ) < 0, we then conclude that ψ has exactly one zero s0 ∈ I which shows the uniqueness of the T -periodic solution s(•). Notice that s(•) is asymptotically stable since φ ( s0 ) < 1 (see [START_REF] Smith | The Theory of Chemostat -Dynamics of microbial competition[END_REF]). Finally, for any s 0 ∈ I \ { s0 }, the sequence (φ n (s 0 )) n is monotone since φ is increasing, and thus it converges to a fixed point of φ . The only fixed point at boundary of I is s in with φ (s in ) < 1 which implies that φ (s 0 ) < s 0 for s 0 close to s in . One can verify that φ n (s 0 ) necessarily converges to s0 when n goes to infinity which implies that |s(t, s 0 )s(t)| → 0 when t goes to infinity.

The next result will play an important role in the following.

) 9
When ν is convex on the interval J, applying Jensen's inequality to ν yields

ν 1 T T 0 s(t) dt ≤ 1 T T 0 ν( s(t)) dt. (10)
Similarly, when ν is concave on J, one obtains the reversed inequality. Since equality in Jensen's inequality holds true if and only if ν is affine over J, the inequality in ( 10) is strict. Combining ( 9)-( 10) then gives

ν 1 T T 0 s(t) dt < ν( s),
and, as ν is increasing, we get

1 T T 0 s(t) dt < s .
The reversed inequality is obtained in the same way when v is concave over J. Lemma 3: The solution s verifies the inequalities:

min t∈[0,T ] s(t) ≤ s ≤ max t∈[0,T ] s(t).
Proof: Equality ( 9) can be equivalently written

1 T T 0 [ν( s(t)) -ν(s )] dt = 0.
Then, the function t → ν( s(t))ν(s ) has to change its sign over [0, T ]. Since ν is monotone over I, we deduce that there exists t -and t + in [0, T ] such that s(t -) < s < s(t + ).

Without any loss of generality, we can then consider periodic solutions of [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] with periodic controls satisfying [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF] and such that s(0) = s . The constant solution s = s for constant D = D is one of them. We show how to use Prop. 2 to discriminate density-dependency.

III. A DISCRIMINATING PROCEDURE

Consider the growth functions ν associated to Monod and (generalized) Contois laws (following the notation of Sec. II):

ν M (s) := µ max s K s + s ; ν CG (s) := µ max s K s + s + kY (s in -s)
.

One can easily check that these functions are increasing. A straightforward computation also gives

ν M (s) = -2K s µ max (K s + s) 3 ; ν CG (s) = µ max 2(kY -1)(K s + kY s in ) (K s + s + kY (s in -s)) 3 .
One can observe that ν M is always concave, while ν CG is convex when kY > 1. We shall say that the growth rate ν has a weak (or null) density effect if it follows the laws of Monod or (generalized) Contois with k ≤ 1/Y . Conversely, we shall say that it has a strong density effect if it follows the law of (generalized) Contois with k > 1/Y .

We investigate now how to use the results of Proposition 2 to discriminate between convex and concave functions using periodic controls. Let us underline that for a given value of s , determining periodic non constant functions D(•) such that the solution s(•) with s(0) = s is periodic requires the perfect knowledge of the function µ and the constants Y and s in . On another hand, waiting for the convergence to a periodic solution could be long and difficult to test in practice. We propose below an experiment procedure to discriminate between weak and strong density effects.

Procedure : It consists in conducting a single chemostat experiment with three phases: constant, bang-bang and constant. It requires the single measurement of the substrate s.

• Phase 1. Conduct an experiment for a constant dilution rate D with 0 < D < ν(s in ) (to avoid the washout equilibrium) and wait the chemostat to be at (quasi) steady state. Let s eq be the measurement of s, and t 0 the current time.

• Phase 2. Choose two values D min , D max such that 0 ≤ D min < D < D max , and two time durations δ 1 > 0, δ 2 > 0. Apply first D = D max on the time interval [t 0 ,t 0 + δ 1 ). Next, apply D = D min until the time t such that s(t) = s eq , and carry on applying D = D min during the duration δ 2 . Finally, change the flow rate value to D = D max for t > t + δ 2 . Store the measurements history {s(t)} t∈[t 0 ,t 0 +T ] until the first time t 0 + T > t + δ 2 such that s(t 0 + T ) = s eq (which exists as 0 ≤ D min < D < D max ). To summarize, one applies

D(t) :=    D max , t ∈ [t 0 ,t 0 + δ 1 ), D min , t ∈ [t 0 + δ 1 , t + δ 2 ), D max , t ∈ [t + δ 2 ,t 0 + T ]. (11)
Let D be the mean value of the dilution rate over [t 0 ,t 0 + T ]:

D := D max (T + t 0 -t + δ 1 -δ 2 ) + D min (t -t 0 + δ 2 -δ 1 )
T .

• Phase 3. Apply for t > t 0 + T the constant dilution rate D := D and wait for the (quasi) steady state, denoted by s . With the data stored on [t 0 ,t 0 + T ], determine the average value of s on [t 0 ,t 0 + T ]:

ŝ := 1 T t 0 +T t 0 s(t) dt.
Conclusion. If ŝ < s , we validate the strong density effect.

Phase 1 allows the chemostat to be at steady state, so that dynamics ( 7) is considered to be valid at time t 0 . The phase 2 generates a T -periodic solution with a periodic control, over a single period, without having to wait the convergence to the periodic solution. Phase 3 compares the average value ŝ of the substrate concentration on the periodic solution with the steady state one s obtained for a constant control D equal to the average value of the periodic control. Proposition (2) allows to conclude.

Remark 1: Note that s has no a priori reason to be equal to s eq , as the initial value of D is arbitrary.

Remark 2: This procedure is a way to discriminate between Monod and (generalized) Contois laws. It allows also to discriminate between convex or concave functions on the invariant interval (λ (D min ), λ (D max )) and can then be applied to other situations whose growth is expected to follow expressions such as [START_REF] Arditi | Saiah Empirical Evidence of the Role of Heterogeneity in Ratio-Dependent Consumption[END_REF]. For this function, one has

ν K (s) = α µ 0 2s in + (α -1)s (s in -s) α+2 .
Therefore, ν K is convex when α is larger than one (which has been observed in [START_REF] Krichen | High biomass density promotes density-dependent microbial growth rate[END_REF]) or when s in is large enough.

The proposed procedure considers bang-bang periodic controls for discriminating convex functions ν. One may wonder if the periodic control that satisfies [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF] allows to obtain the largest difference between ŝ and s . We justify this choice in the next section.

IV. THE BEST SHAPE OF PERIODIC CONTROLS

The goal of this section is to justify the construction of D in [START_REF] Contois | Kinetics of Bacterial Growth: Relationship between Population Density and Specific Growth Rate of Continuous Cultures[END_REF] as a bang-bang control, in terms of optimal control. Given bounds D min , D max on the control D(•) with 0 ≤ D min < D < D max , and given a convex function ν, we look for Tperiodic controls D(•) that satisfy [START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF] with

D(t) ∈ [D min , D max ] a.e. t ≥ 0, (12)
such that the associated solution s with s(0) = s is Tperiodic and gives the smallest average value ŝ: min

D(•) ŝ := 1 T T 0 s(t) dt. (13)
Similarly, when ν is concave, we would look for controls D(•) such that ŝ is the largest. Problem ( 13) is an optimal periodic control problem. Several contributions about the theory of optimal control with periodic controls are available in the literature (see for instance [START_REF] Bittanti | Status of periodic optimization of dynamical systems[END_REF], [START_REF] Gilbert | Optimal periodic control: A general theory of necessary conditions[END_REF], [START_REF] Gilbert | Necessary and Sufficient Conditions for Local Optimality of a Periodic Process[END_REF] and references herein). However, in the present case, it has to fulfill the additional integral constraint (6), which has not been yet studied, up to our knowledge, expect in the recent works [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF], [START_REF] Bayen | Optimal periodic control of the chemostat with Contois growth function[END_REF].

Let K := (λ (D min ), λ (D max )) and introduce the function

η : K → R * + as η(s) := 1 (s in -s)(D max -ν(s))(ν(s) -D min )
, s ∈ K.

One has the following optimality result. Proposition 4: Assume ν to be convex increasing on K, and take arbitrarily D ∈ (D min , D max ) and T > 0. Then: (i) There are exactly two T -periodic optimal solutions of (13) under the constraints ( 6) and ( 12), and such that s(0) = s . (ii) The controls D ± generate the two optimal trajectories (with same cost):

D + (t) := D max , t ∈ [0,t + 1 ) ∪ [t + 2 , T ) mod T, D min , t ∈ [t + 1 ,t + 2 ) mod T, (14)
D -(t) := D min , t ∈ [0,t - 1 ) ∪ [t - 2 , T ) mod T, D max , t ∈ [t - 1 ,t 2
-) mod T, where times t ± i , i = 1, 2 are given by:

t + 1 := inf{t ≥ 0 ; s(t) = s M }, t + 2 := inf{t ≥ t + 1 ; s(t) = s m }, t - 1 := inf{t ≥ 0 ; s(t) = s m }, t - 2 := inf{t ≥ t - 1 ; s(t) = s M }, (15)
and where (s m , s M ) is the unique pair in K 2 that satisfies

s M s m η(s) ds = T D max -D min , s M s m η(s)ν(s) ds = DT D max -D min . (16)
Proof: We apply a result of [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] that we recall in the Appendix. One can easily check that ( 13) is a particular instance of [START_REF] Harmand | The Chemostat: mathematical theory of microorganisms cultures[END_REF] taking and (s) = s, and verify that the required assumptions to apply the result of the Appendix are satisfied. The optimal solutions u ± of (18) provided by [START_REF] Krichen | High biomass density promotes density-dependent microbial growth rate[END_REF] once translated into [START_REF] Desmond-Le Qumner1 | A thermodynamic theory of microbial growth[END_REF], lead to the optimal solutions D ± of [START_REF] Desmond-Le Qumner1 | A thermodynamic theory of microbial growth[END_REF] given by ( 14)- [START_REF] Gilbert | Optimal periodic control: A general theory of necessary conditions[END_REF]. As well, [START_REF] Gilbert | Necessary and Sufficient Conditions for Local Optimality of a Periodic Process[END_REF] straightforwardly follows from [START_REF] Hoskisson | Continuous culture -making a comeback[END_REF].

I = K, η = η, f (s) = (s in - s) D max +D min 2 -ν(s) , g(s) = (s in -s) D max -D min
This result gives the existence of δ 1 , δ 2 that give the largest value | ŝs | in the procedure of Sec. III. Let us underline that s is not known in advance (as the function ν is unknown) but is learned in phase 2 of the procedure.

V. NUMERICAL ILLUSTRATIONS

We have considered Monod and Contois laws (see the graph of the associated functions ν on Fig. 1) that provide (for the same values of D and s in ) the same steady state s eq = 1. It is then not possible to distinguish between these two functions without changing the value of D. Fig. 1: In blue: Monod's kinetics for K s = 2; Y = 1; s in = 2; µ max = 2 (plain) and µ max = 1.4 (dashed). In red: Contois's kinetics for K s = 0; Y = 1; k = 2; s in = 2; µ max = 3 (plain); µ max = 2 (dashed) and s in = 3 (dotted).

illustrates the three phases generated by the procedure (phase 1: blue, phase 2: green, phase 3: red) with the average value ŝ (in black). Then Tables I and II (left columns) gives The interest of the procedure, compared to classical growth identification, relies on its robustness with respect to disturbances. It is well known that continuous cultures are often subject to uncontrolled environmental changes, such as variations of s in , temperature or pH which directly impact the maximum growth rate µ max . We have simulated such a sudden change during the phase 2 of the procedure. These changes correspond to the functions ν in dashed line on Fig. 1 and modify the s-trajectory as depicted in dashed line on Fig. 2. Tables I and II present the results under two kinds of disturbances: change of s in and change of µ max . In any case, one can observe that the sign of ∆s is preserved. Notice that under these changes, the duration T is modified and the ssolution between t 0 and t 0 + T is no longer periodic, neither for the old or the new functions ν, but ∆s is a continuous function of the parameters s in , µ max , justifying the robustness of the test.

VI. CONCLUSION

We have proposed an experiment procedure to test a density dependency in the chemostat, using non-constant dilution rates. This procedure consists in a sequence of a steady state, a periodic operation over only one period, and a last steady-state, with the single measurement of the substrate concentration. We have shown its robustness with respect to disturbances on the growth curves during the procedure, demonstrating its interest compared to classical identification methods. For more complex growths for which the function µ could be neither convex or concave on their whole domain (such as the Hill function [START_REF] Moser | The Dynamics of Bacterial Populations Maintained in the Chemostat[END_REF] or the one from the Microbial Transition State theory [START_REF] Desmond-Le Qumner1 | A thermodynamic theory of microbial growth[END_REF]), one would need to choose adequately the bounds D min , D max in the procedure. This is an open problem that could be the matter of a future work.

VII. APPENDIX

We recall results from [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] about optimal periodic control. Consider the following dynamical system: ẋ = f (x) + u(t)g(x), u(t) ∈ [-1, 1] and x(0) = x, [START_REF] Harmand | Density-dependent kinetics models for a simple description of complex phenomena in macroscopic massbalance modeling of bioreactors[END_REF] where f : I → R, g : I → R * + and I is an open interval of R. f and g are assumed to be C 1 with f + g > 0 and fg < 0 over I . Consider now the following optimal control problem:

min u(•) 1 T T 0 (x(t)) dt, (18)
where : I → R and x(•) is a T -periodic solution of ( 17) associated with a measurable control u(•) satisfying

1 T T 0 u(t) dt = ū,
where ū :=f ( x)/g( x). Define ψ : I → R as ψ := -f g and γ := ψ • -1 whenever -1 is well-defined, as well as η(x) := 1 f (x) + g(x) -1 f (x)g(x)

, x ∈ I .

It is supposed in [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] that I is invariant by the dynamics [START_REF] Harmand | Density-dependent kinetics models for a simple description of complex phenomena in macroscopic massbalance modeling of bioreactors[END_REF] with ψ(I ) ⊂ [-1, 1], and that there is a unique x ∈ I s.t.

(ψ(x)ψ( x))(xx) > 0, x ∈ I , x = x.

It is also shown in that when is increasing and γ strictly convex increasing, there is a unique (x m , x M ) ∈ I 2 s.t.

Theorem 3.6 in [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input[END_REF] shows that the bang-bang control whose restriction over [0, T ] is defined as

u + (t) := +1 if t ∈ [0, t1 ], -1 if t ∈ [t 1 , t2 ], +1 if t ∈ [t 2 , T ], (20)
(where t1 is the first time x reaches x M and t2 the first time t > t1 where x reaches x m ) is optimal (up to a unique time translation defining an optimal control u -with same cost).

  The key point is to consider functions D(•) that fulfill[START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF] for the same value D. Hereafter, for a given non constant T -periodic control D(•) that fulfills[START_REF] Bastin | On-Line Estimation and Adaptive Control of Bioreactors[END_REF], we denote by s(•) the unique periodic solution associated with D(•). Proposition 2: If the restriction of ν(•) to the interval J := {ν( s(t)), t ∈ [0, T ]} is convex, respectively concave, and does not coincide with a linear function on this interval, then ŝ := 1 T T 0 s(t) dt < s , respectively ŝ := 1 T T 0 s(t) dt > s . Proof: Let s(•) be a non-constant T -periodic solution with values in I. From (8), we deduce the equality 1 T T 0 ν( s(t)) dt = D. (

2

 2

Fig. 2 : 4 δ 1

 241 Fig. 2: Solution s(•) generated by the procedure ( for the Monod law)with s(0) = 0.2; K s = 1; µ max = 2; Y = 1; s in = 2; D max = 1.3; D min = 0.2 (in dashed under the perturbation) the difference ∆s := ŝs provided by the procedure, for different values of the durations δ 1 , δ 2 . One can see that a

2 δ 1 = 10 δ 2 = 8 ∆s 22

 212822 s in = 2 and 3, s in = 2, µ max = 3 µ max = 3 µ max = 3 and = -0.20 ∆s = -0.55 ∆s = -0.32 T = 19.5 T = 19.73 T =

  )ψ(x) dx = T ū.

TABLE I :

 I Results for Monod law: K S = 2; Y = 1; D max = 1.3; D min = 0.2.

	Without	Perturbation	Perturbation
	perturbation:	on s in :	on µ max :
	s in = 2,

TABLE II :

 II Results

for Contois law: k = 2; Y=1; D max = 2.7; D min = 0.1.

VIII. ACKNOWLEDGMENTS

The authors thank the Algeria-France program AVER-ROES for the PhD grant of F.-Z. Tani, the support of the LabEx NUMEV incorporated into the I-Site MUSE and the ANR project THERMOMIC, that gave us the idea to develop this procedure.


Fatima Zahra
email: fatima-zahra.tani@supagro.fr

Alain Rapaport
email: alain.rapaport@inra.fr

Térence Bayen
email: terence.bayen@umontpellier.fr

F.-Z Tani

Mistea

A hybrid control against species invasion in the chemostat

A hybrid control against species invasion in the chemostat

Fatima-Zahra Tani, Alain Rapaport and Térence Bayen Abstract-In this work, we consider the classical chemostat model with the objective to limit the invasion of a new species having negative effect on the resident one, playing with the removal rate. We study the resilience of the system to the apparition of an invader, and propose a new concept of weak resilience when the system cannot return and stay at its original state (whatever is the removal rate). The weak resilience guarantees that the measure of the time for which the density of the resident species is above a given threshold is infinite. We show that this can be achieved by a hybrid controller with very few knowledge on the dynamics of the system. As it is not possible to eradicate totally the invasive species, the controller makes the resident species return indefinitely many times above the desired threshold, and the solutions converge asymptotically to periodic solutions. We illustrate the effectiveness of the proposed controller on numerical simulations.

I. INTRODUCTION

The well known chemostat model [START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF], [START_REF] Harmand | The chemostat, mathematical theory of continuous culture of micro-organisms[END_REF] describes the microbial growth in continuous culture [START_REF] Panikov | Microbial Growth Kinetics[END_REF] and is also often used to represent continuously stirred tank reactors [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]:

ṡ = -1 Y µ(s)x + D(s in -s), ẋ = (µ(s) -D)x. (1)
The variables s and x denote the concentrations of substrate and biomass, respectively. The operating parameters are the input concentration of substrate s in and the removal rate D (we assume that the effective removal rate is the same on both concentrations, which amounts to neglect mortality or attachment of bacteria). The parameter Y and the function µ(•) are the yield conversion and the specific growth rate, that we assume here to be increasing with µ(0) = 0 (as it is often the case in microbiology [START_REF] Panikov | Microbial Growth Kinetics[END_REF]). Without loss of generality, one can take Y = 1 (at the price to change the biomass unit). In many industrial applications, bacteria are expected to provide ecosystem services (such as water purification) or produce substances of economic interest (in pharmacology for instance). It is often of primer importance to maintain the concentration of the biomass x at or above a given threshold. One can easily check that a positive equilibrium (s , x ) of (1) verifies µ(s ) = D and x = Y (s ins ). Typically, bioreactors are operated at steady state choosing the removal rate D to be equal to a nominal value D := µ(s inx /Y ) where x is the desired biomass concentration (see for instance [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF], [START_REF] Harmand | The chemostat, mathematical theory of continuous culture of micro-organisms[END_REF]). However, it is also well known that open systems are subject to bacterial contamination or genetic evolution. Therefore, another bacterial species can appear in the chemostat system at any time. The chemostat model with more than one species (here two, assuming without loss of generality that the yield coefficients Y i are equal to 1):

   ṡ = -µ 1 (s)x 1 -µ 2 (s)x 2 + D(s in -s), ẋ1 = (µ 1 (s) -D)x 1 , ẋ2 = (µ 2 (s) -D)x 2 , (2)
has been widely studied in the literature [START_REF] Humphrey | Dynamics of a chemostat in which two organisms compete for a common substrate[END_REF], [START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF], [START_REF] Grover | Resource Competition[END_REF], [START_REF] Grivet | Nonlinear population dynamics in the chemostat[END_REF], [START_REF] Harmand | The chemostat, mathematical theory of continuous culture of micro-organisms[END_REF]. We also assume that the removal rates of bacteria are identical). The main result is the so-called "Competitive Exclusion Principle" [START_REF] Hsu | A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms[END_REF], [START_REF] Armstrong | Competitive exclusion[END_REF], [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], which states that (generically) at most one species survives at steady state. Then, depending on the characteristics of the invading species, this one can settle and eradicate the resident species, which can be catastrophic for operators if the new species does not provide similar ecosystem services or by-products. However, competing species often present compromises. Typically, a species could be specialized for low substrate concentrations while another one could be for larger concentrations. In such cases, it has been shown that applying a periodic timevarying removal rate D(•) could allow coexistence of both species [START_REF] Smith | Competitive coexistence in an oscillating chemostat[END_REF], [START_REF] Butler | A mathematical model of the chemostat with periodic washout rate[END_REF], [START_REF] Wolkowicz | N-species competition in a periodic chemostat[END_REF]. It appears that the design of such time-varying controls requires the perfect knowledge of the growth functions µ i (•). In addition, to our best knowledge, the performances in terms of services that continue to be provided by the resident species, have not been studied in the literature. This is precisely the objective of the present work. In particular, we shall see that hybrid controllers can be defined without a precise knowledge of the functions µ i (•) allowing coexistence of both species. In addition, for the corresponding solution, the time spent by the density of the resident species above a given threshold is infinite. The paper is organized as follows. In Section II, we setup assumptions and study the resilience of the chemostat model [START_REF] Aubin | Viability Theory[END_REF] in presence of an invading species, considering system (2). We introduce here a new concept of weak resilience.

In Section III, we present recent results about the design of time-varying removal rates that make the system weak resilient. We then show that such open-loop controls can indeed be synthesized by a hybrid feedback controller. Section IV illustrates with numerical simulations the behavior of the proposed controller and discusses the role of its parameters.

II. RESILIENCE ANALYSIS

For system (2), let us consider the following (classical) assumptions:

Assumption A1: The functions µ i (•) (i = 1, 2) are C 1 in- creasing functions with µ i (0) = 0.
Let us also introduce the break-even concentrations:

λ i (D) := sup {s ∈ [0, s in ] ; µ i (s) < D} , i = 1, 2.
Then, the Competitive Exclusion Principle (CEP) [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF], [START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF], [START_REF] Harmand | The chemostat, mathematical theory of continuous culture of micro-organisms[END_REF] states as follows.

Theorem 1 (CEP): Under Assumption A1, if λ i (D) < min(λ j (D), s in ) for i = j, i, j ∈ {1, 2}, then any positive solution of (2) converges asymptotically to the equilibrium (s , x 1 , x 2 ) with s = λ i (D), x i = s inλ 2 (D) and x j = 0. Therefore, several cases appear: -if µ 2 (s) > µ 1 (s) for any s > 0, then one has λ 2 (D) < λ 1 (D) for any D ∈ (0, µ 1 (s in )), and the resident species never survives under invasion by species 2.

-On the opposite, if µ 2 (s) < µ 1 (s) for any s > 0, the resident species is safe in presence of such an invader (with D < µ 1 (s in ) to avoid the washout equilibrium).

-However, as mentioned in the introduction, it often happens that species are complementary: if the invader can settle and eradicate the resident species for small values of s, the issue of the competition is reversed under large concentrations of s, what we shall consider in the sequel with the following assumption:

Assumption A2: There exists s ∈ (0, s in ) such that (µ 1 (s) -µ 2 (s))(s -s) > 0, s ∈ [0, s in ] and s = s. (3) Denote then D := µ 1 ( s) = µ 2 ( s).
For any D ∈ ( D, µ 1 (s in )), the resident species 1 wins the competition and its density at steady state is equal to x * 1 := s inλ 1 (D). Therefore the number

x1 := s in -λ 1 ( D) = s in -s > 0.
is the largest density level that the system can maintain in presence of the invader 2 under constant removal rate. However, one may wonder if it is possible to maintain higher levels with time varying removal rate

D(•) taking values in an interval [D m , D M ] with 0 < D m < D < µ 1 (s in ) < D M . (4)
One can easily check that the domain

{(s, x 1 , x 2 ) ∈ R 3 + ; s + x 1 + x 2 = s in },
is invariant by ( 2) and attractive for any persistently exciting control D(•). By persistently exciting control, we mean any non-negative function D(•) such that +∞ 0 D(t)dt = +∞ (see [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]). Considering that ( 1) is already at steady state before invasion, and that the quantity x 2 at invasion is very small, we shall consider in the sequel the reduced dynamics on this domain, that is, the planar dynamics

ẋ = f (x, D) := (µ 1 (s in -x 1 -x 2 ) -D)x 1 (µ 2 (s in -x 1 -x 2 ) -D)x 2 , (5)
defined on the invariant set

S := {x ∈ R 2 + ; x 1 + x 2 ≤ s in }. Consider then a threshold x r 1 ∈ ( x1 , s in ) satisfying 0 < D m < µ 1 (s in -x r 1 ), (6)
and define a subset K (x r 1 ) of S as:

K (x r 1 ) := {x ∈ S ; x 1 ≥ x r 1 and x 2 > 0}.
Let us now formulate the invasion problem in terms of resilience. Assume that, at the initial time, species 1 is alone with a density level equal to x r 1 ensured by the choice of the constant dilution rate D r := µ 1 (s in -x r 1 ) < D, and that species 2 (satisfying Assumption A2) appears. We then say that system [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] with the constant removal rate D r is not resilient with respect to the set K (x r 1 ) because the trajectory escapes from K (x r 1 ) and never comes back (indeed, it converges to the equilibrium point (λ 2 (D r ), 0, s inλ 2 (D r )).

Let us now study resilience when one allows time-varying removal rate, in terms of viability analysis. Let us recall that the viability kernel of a closed subset K of R n for a controlled dynamics ẋ = f (x, u(t)) (with f : R n × R m → R n and u(t) ∈ U ⊂ R m ), is the largest closed subset V of K such that for any x 0 ∈ V there exists u(•) for which the solution with x(0) = x 0 verifies x(t) ∈ K for any t ≥ 0 (see [START_REF] Aubin | Viability Theory[END_REF]). 4) and ( 6)) is given by

Viab(cl K (x r 1 )) = [x r 1 , s in ] × {0}. Proof: From (6), x r 1 is such that x r 1 ∈ ( x1 , s in -λ 1 (D m ))
. Moreover, Assumptions A1, A2 and condition (4) also imply the inequality 0 < λ 2 (D m ) < λ 1 (D m ) < s. Take an initial condition (x 1,0 , x 2,0 ) in cl K (x r 1 ). If x 2,0 = 0, the solution of (5) verifies x 2 (t) = 0 for any t > 0 and any control D(•). From (6), we can choose a constant control D such that D m ≤ D < µ 1 (s inx r 1 ), and we observe that (5) with D(t) = D satisfies

x 1 = x r 1 ⇒ ẋ1 = (µ 1 (s in -x r 1 ) -D)x r 1 > 0, x 1 = s in ⇒ ẋ1 = -Ds in < 0.
Thus, the constant control D allow the trajectory to stay in the desirable set K (x r 1 ). Assume now that x 2,0 > 0. Then, a solution x(•) of ( 5) associated with an admissible control D(•) verifies x 2 (t) > 0 for any time t ≥ 0. Suppose by contradiction that x(•) stays in K (x r 1 ) for any time t ≥ 0. Then one has s(t) < s inx r 1 < s for any time t ≥ 0. By Assumption A2, one has for any t ≥ 0, µ 1 (s(t))µ 2 (s(t)) < 0 and thus, we deduce the inequality

ṡ(t) > -µ 2 (s(t))x 1 (t) -µ 2 (s(t))x 2 (t) + D m (s in -s(t)) = (D m -µ 2 (s(t)))(s in -s(t)),
for any time t ≥ 0. Notice that any positive solution

ζ (•) of ζ = (D m -µ 2 (ζ (t)))(s in -ζ (t)) converges to λ 2 (D m ) when t → +∞. From (6) one has λ 2 (D m ) < s in -x r
1 , hence there exist t 1 ≥ 0 and s ∈ (0, λ 2 (D m )) such that one has ζ (t) > s for any t ≥ t 1 . From the preceding inequality, we deduce that s(•) satisfies s(t) ≥ ζ (t) for any time t ≥ 0. We thus deduce that for any time t ≥ t 1 , one has µ 1 (s(t))µ 2 (s(t)) ≤ c with

c := min{µ 1 (σ ) -µ 2 (σ ) ; σ ∈ [s, s in -x r 1 ]} < 0. If we differentiate the function q 1 := x 1 /x 2 w.r.t. t, we find that q1 = µ 1 (s(t)) -µ 2 (s(t)) q 1 ,
with s(t) = s inx 1 (t)x 2 (t). One then obtains q1 < c q 1 . Therefore q 1 decreases to zero and x 1 as well, leading to a contradiction. We conclude that the only solutions of ( 5) that stay in cl K (x r 1 ) for any time are the ones starting with x 2,0 = 0 as was to be proved.

This lemma shows that no control strategy allows the resident species to stay in the desirable set K (x r 1 ) in presence of the invader. The dynamics is thus not resilient for the domain K (x r 1 ). We are now in a position to introduce a weaker concept of resilience.

Definition 2: Let x r 1 ∈ ( x1 , s in -λ 1 (D m )).
The system ( 5) is said to be weakly resilient for the set K (x r 1 ) if for any initial condition in K (x r 1 ), there exists a control function D(•) with values in [D m , D M ] such that the corresponding solution of (5) satisfies meas {t ≥ 0 ; x(t) ∈ K (x r 1 )} = +∞. Such a function D(•) will be called a weakly resilient control.

The purpose of the next section is to construct explicitly such control functions.

III. THE HYBRID CONTROL

In this section, we shall consider bang-bang controls with values D m or D M (parameters D m and D M satisfying conditions (4) and ( 6)).

A. Construction of a weakly resilient control

Let us first give properties of the dynamical system (5) when D is constant (D = D m or D = D M ).

1. For D = D M , the "washout" equilibrium E 0 := (0, 0) is the single equilibrium in S , which is globally asymptotically stable (see Fig. 1).

2. For D = D m : there are three equilibria in S :

E 0 , E 1 := (x m 1 , 0) and E 2 := (0, x m 2 ) with x m i := s in -λ i (D m ).
One can easily check that E 2 is a stable node while E 1 is a saddle point. By the Theorem of the stable and unstable manifolds, we deduce that the unstable manifold W u of E 1 connects E 1 to E 2 , and that the stable manifold of E 1 belongs to the x 2axis (see Fig. 2).

One can also easily check that solutions of (5) with D = D M reach the origin tangentially to the x 2 -axis. Therefore, if D(t) = D M for a sufficiently long time, the state will get so close to the x 2 -axis that a switching to D = D m will produce a further trajectory that will remain close from the stable manifold of E 1 for a long time and that can then reach the set K (x r 1 ), before moving towards the equilibrium E 2 . Then, another switch to D = D M when the state leaves K (x r 1 ) allows to come back again close to the origin, and so on.

The construction of a weakly resilient control precisely relies on asymptotic properties of (5) with a constant control D = D m or D = D M . The next theorem is our main result and formalizes this methodology.

Theorem 2: Assume that D m and D M satisfy (4) and [START_REF] Bastin | On extremum seeking in bioprocesses with multivalued cost functions[END_REF]. For any ε > 0, define the set

E ε := (0, x1 ] × (0, ε].
For ε small enough and any initial condition x(0) ∈ K (x r 1 ), there exists a piecewise constant function D(•) which alternates the values D M , D m on time intervals [T i , T i+1 ), i ∈ N satisfying: [START_REF] Bayen | About Moreau-Yosida regularization of the minimal time crisis problem[END_REF] where x(•) is the unique solution of (5) associated with

T 0 = 0 < T 1 < • • • < T i < T i+1 < • • • and lim i→∞ T i = +∞,
D(•) such that (i) if x(T i ) / ∈ E ε , one has D(t) = D M for t ∈ [T i , T i+1 ) with T i+1 defined as the first next entry time into E ε . (ii) if x(T i ) ∈ E ε , one has D(t) = D m for t ∈ [T i , T i+1
) and the trajectory x(•) enters into the set K (x r 1 ) in finite time ; T i+1 is defined as the first next exit time from K (x r 1 ). Moreover, D(•) is a weakly resilient control.

Du to lack of space, we do not provide the proof of this result (which is long and technical). It mainly consists in studying deep properties (such as image, continuity) of the operator:

O : (0, s in -x r 1 ] → (0, s in -x r 1 ] x 2,0 → x 2 (T 2 )
where x(•) = (x 1 (•), x 2 (•)) denotes the unique solution of ( 5) for the initial condition (x r 1 , x 2,0 ) and the time-varying control D(•) given by Th. 2, parameters D m , D M , ε being fixed. More details can be found in [START_REF] Bayen | Weak resilience of the chemostat model to a species invasion with non-autonomous removal rates[END_REF].

Let us underline the fact that the control law given by Th. 2 is an open-loop control, and that the computation of the switching times T i requires the perfect knowledge of the functions µ j (•). In practice, an open-loop control is not robust with respect to the knowledge of the initial condition and the switching times. Moreover the characteristics of the invasive species is usually not known in advance. However, this control strategy cannot be written as a pure state feedback because when x 1 reaches the value x r 1 with x 2 > ε, one cannot decide if D = D m or D = D M without knowing the past, i.e., if x 1 is increasing or decreasing.

B. Synthesis with a hybrid control

Instead, we propose a hybrid controller associated with a logic-based switching algorithm represented by a piecewise signal σ (•) that can take three possible states denoted by χ 1 , χ 2 and χ 3 . The transitions between these three states is illustrated on Fig. 3.

x 1 < x r 1 D = D M D = D m χ 1 χ 3 χ 2 D = D m x 1 > x r 1
x ∈ E Fig. 3: The logic-based algorithm associated with the hybrid control.

One has the following result.

Proposition 3: Let σ (•) : R + → {χ 1 , χ 2 , χ 3 } be such that σ (t) := g(x 1 (t), σ (t -)) =        χ 2 if σ (t -) = χ 1 and x 1 (t) > x r 1 , χ 3 if σ (t-) = χ 2 and x 1 (t) < x r 1 , χ 1 if σ (t -) = χ 3 and x(t) ∈ E ε , σ (t -) otherwise,
with σ (0) = χ 3 . Then, the σ -feedback:

D(σ (t)) = D m if σ (t) ∈ {χ 1 , χ 2 }, D M if σ (t) = χ 3 . ( 8
)
guarantees the weak resilience of the dynamics with the control [START_REF] Bayen | About the minimal time crisis problem[END_REF].

Proof: This automata is a simple way to memorize if the state x 1 is entering or leaving the set K (x r 1 ). One can easily see that solutions of the coupled dynamics

ẋ(t) = f (x(t), D(σ (t))), x(0) ∈ K (x r 1 ), σ (t) = g(x 1 (t), σ (t -)), σ (0) = χ 3 ,
(see for instance [START_REF] Liberzon | Switching in Systems and Control[END_REF] for the solution concept) are such that the times T i given by Th. 2 correspond exactly to the ones of σ (when switching from χ 2 to χ 3 or from χ 3 to χ 1 ). Provided that ε is small enough and that D m and D M satisfy (4) and ( 6), Th. 2 then guarantees weak resilience of the hybrid controller [START_REF] Bayen | About the minimal time crisis problem[END_REF].

In [START_REF] Bayen | Weak resilience of the chemostat model to a species invasion with non-autonomous removal rates[END_REF], it has been shown that there exist periodic solutions of (5) associated with the control given in Th. 2. It has been conjectured that for ε sufficiently small, the operator O : x 2,0 → x 2 (T 2 ), where T 2 is given by Th. 2 for the initial condition (x r 1 , x 2,0 ), is contractive. Then, under this condition, it has been proved that any trajectory with initial condition in K (x r 1 ) and the weakly resilient control given by Th. 2 converges asymptotically to a unique periodic solution x † (•) with a period T † , up to a time shift.

Remark 1:

The synthesis of the hybrid controller requires the single choice of the parameters D m and D M that satisfy (4) and ( 6) (which are quite loose conditions) and ε > 0 small enough. Therefore, it does not require the precise knowledge of the functions µ i to ensure weak resilience against a species that fulfills Assumption A2, in contrast with more sophisticated approaches, such as model predictive control, which rely on more knowledge on the growth functions. Let us also underline that the controller switches between an environment that is unfavorable to both species (D = D M ) and an environment which is favorable to the invasive species (D = D m ). It is then not intuitive that the resident species could be dominant most of the time under such a switching. This property is strongly linked to the choice of ε that has simply to be sufficiently small, although we are not able to provide an explicit bound.

IV. NUMERICAL SIMULATIONS Numerical simulations have been performed to illustrate the behavior of the trajectories generated by the hybrid controller and the role of its parameters on the proportion of time spent in the desired set K (x r 1 ), as well as their impact on the productivity of the resident species.

We consider two growth functions of Monod's type: Let us first notice that in absence of species 2, the hybrid controller gives D(t) = D m . We have simulated an invasion of species 2 at time t I = 50 with a sudden input of a small amount of concentration x 2 (t + I ) = 0.2. Before the invasion, one can observe that the system is at quasi-steady state and after the invasion, the control D(t) alternates between D m and D M . The solution finally converges asymptotically to a periodic solution x † (•). A. Proportion of time spent in K (x r 1 ) On Figs. 5 and6, two different values of ε have been chosen. In Table I, we also present for different values of ε, the proportion of time spent by the periodic solution x † (•) the desirable set K (x r 1 ), during one period (other parameters being fixed). When ε decreases, one observes that the period increases, as it requires a longer time for the trajectory to reach the set E ε . Interestingly, one can also observe that the proportion of time spent in K (x r 1 ) increases when ε decreases. Indeed, when ε is small, the trajectory gets close to the stable manifold of the saddle equilibrium E 1 when D = D m (which is the x 1 -axis) and therefore the state stays a long time in the vicinity of E 1 which belongs to K (x

B. Productivity of the resident species

Let us consider the productivity P of the species 1 alone, as the quantity produced per unit of volume and time when the system is at steady steady state. It is given by P := Dx eq 1 (D), where x eq 1 (D) is the steady state associated with a constant control D (or equivalently P = µ 1 (s eq )(s in -s eq )). The largest value of P is then obtained when the substrate concentration at steady state that maximizes the function s → µ 1 (s)(s ins) at a certain s ∈ (0, s in ). Then, the maximal productivity is obtained for the dilution rate D = µ 1 (s ) and one has x 1 = s in -s . For the chosen value of s in , the maximal productivity of species 1 is obtained for P max 0.4289 (see Fig. 7). We have then chosen a threshold x r 1 equal to x 1 and a value of D m close to D (such that inequality (6) fulfills), so that it does not impact too much the productivity of species 1 as long as species 2 is absent (see Fig. 4). In absence of species 2, as the controller takes the value D m , the productivity of species 1 alone at steady state is P m := D m x m 1 0.4216 (which is quite close to the maximal one).

In presence of species 2, we have considered the average productivity P := 1

T † T † 0 D(t)x † 1 (
t) dt of species 1 over a period of the asymptotic periodic solution x † (•), to be compared with the productivity of species 1 alone (see Table I). The percentage of loss of productivity P m - P P m is also indicated. It is clear that the presence of the species 2 impacts the productivity of species 1, as there is a new consumer of the resource. However, one can see that it is possible to maintain a relatively low decrease of productivity, under the condition that the parameter ε is not too small.

In conclusion, as already known in the literature in the context of single species (see for instance [START_REF] Bastin | On extremum seeking in bioprocesses with multivalued cost functions[END_REF]), one can appreciate the trade-off between having a high productivity and maintaining a high conversion. The parameter ε is thus a lever of choice for the practitioners.

V. CONCLUSION

In this work we have proposed a hybrid controller that allows weak resilience in the chemostat model. This controller switches between bounds D m and D M , so that the corresponding solution enters infinitely many times into a desirable subset (while each bound is unfavorable to resilience when considered alone). One of the main advantages of this controller is that it does not require the precise knowledge of the growth characteristics of the invader and can cope with a large variety of unknown invading species. In practice, one can simply apply this controller even in absence of invading species, guaranteeing a robustness of the performances against possible future invasions. We have also shown its features in terms of time spent in the desirable subset and productivity thanks to an adequate choice of the bounds D m and D M and of the parameter ε, although we do not know if the bang-bang strategy is optimal for the so-called "time of crisis" (see [START_REF] Bayen | About Moreau-Yosida regularization of the minimal time crisis problem[END_REF], [START_REF] Bayen | About the minimal time crisis problem[END_REF], [START_REF] Bayen | Minimal time crisis versus minimum time to reach a viability kernel: a case study in the prey-predator model[END_REF]): This amounts to minimize w.r.t. admissible controls the time spent by trajectories outside the desirable set K (x r 1 ) either over a finite given horizon [0, T ] (which is more relevant in the present context), or over [0, +∞) using the notion of finitely optimal control as in [START_REF] Aseev | Maximum principle for infinite-horizon optimal control problems under weak regularity assumptions[END_REF]. This could be the matter of a future work. When practitioners possess more information on the growth functions (such as average and variability, especially on the possible invasive species), it would be interesting to compare the performances of the proposed controller with more classical approaches such as model predictive control. This could be also studied in a future work. Finally, an extension of this work could deal with different removal rates of species.

Lemma 1 :

 1 The viability kernel Viab(cl K (x r 1 )) of the closure of the set K (x r 1 ) for (5) with controls D(•) ∈ [D m , D M ] (where D m , D M fulfill (

Fig. 1 :Fig. 2 :

 12 Fig. 1: Phase portrait of (5) with the constant control D = D M .

µ 1 (

 1 + s , with s in = 5, see Fig.4.

λ 1 (Fig. 4 :

 14 Fig.4: Graphs of the functions µ 1 and µ 2 with s in = 5.

Fig. 5 : 1 εFig. 6 :

 516 Fig. 5: Time courses for ε 0.1, D m = 0.13 and D M = 2.

Fig. 7 :

 7 Fig. 7: Graph of the function s → µ 1 (s)(s ins) with s in = 5.

  inf D(•) meas {t ∈ [0, T ] ; x(t) / ∈ K (x r 1 )} .

TABLE I :

 I Asymptotic performances with D m = 0.13 and D M = 2.

					r 1 ).
	ε	period	% of time in K (x r 1 ) productivity average	% of loss
	0.001	293	78.50	0.3562	15.51
	0.01	177	75.20	0.3596	14.71
	0.07	80.3	63.65	0.3674	12.86
	0.15	40.5	36.67	0.3755	10.93


Hibridni Sol-Gel-Nanosi

S Dopirani

Z A Cerijem

Korozijsko Za

^ito Magnezijevih

Zlitin Noé

Verner Murillo-Gutiérrez

Florence Ansart

Jean-Pierre Bonino

Marie-Jöelle Menu

Marie Gressier

HYBRID SOL-GEL COATINGS DOPED WITH CERIUM TO PROTECT MAGNESIUM ALLOYS FROM CORROSION

Keywords: magnesium, coating, sol-gel, corrosion inhibitor, EIS

Hybrid coatings produced via the sol-gel route were deposited onto an Elektron 21 magnesium alloy. The sol consisted of tetraethyl-orthosilicate (TEOS) and 3-(trimethoxysilyl)propylmethacrylate (MAP) to which corrosion inhibitors were added. The influence of the cerium concentration on the anti-corrosion properties of the hybrid coating is presented. Furthermore, the morphology of the organic/inorganic coatings deposited on the magnesium alloy was determined with scanning electron microscopy (SEM). In parallel, the electrochemical behavior during the immersion in a 0.05 M NaCl corrosive solution was studied with electrochemical impedance spectroscopy (EIS). It was proven that the hybrid films exhibit a high impedance modulus during the first hours of the immersion and that an addition of cerium to the sol with a concentration of 0.01 M considerably increases the durability of the film, delaying its degradation during the immersion. In addition, this project especially focusses on determining the critical concentration of the cerium salt at which the impedance modulus of the hybrid coating strongly decreases during the immersion.

Hibridni nanosi, pripravljeni s sol-gel-postopkom, so bili naneseni na magnezijevo zlitino Elektron 21. Osnova je bila tetraetil-ortosilikat (TEOS) in 3-(trimetoksisilil)propilmetakrilat (MAP), ki so ji bili dodani inhibitorji korozije. Predstavljen je vpliv koncentracije cerija na protikorozijske lastnosti hibridnega nanosa. Poleg tega je bila dolo~ena morfologija organskih/ neorganskih nanosov na magnezijevo zlitino z vrsti~nim elektronskim mikroskopom (SEM). Vzporedno je bilo preu~evano elektrokemijsko vedenje med potopitvijo v korozijsko raztopino 0,05 M NaCl, z uporabo elektrokemijske impedan~ne spektroskopije (EIS). Dokazano je bilo, da izkazujejo hibridni nanosi visok impedan~ni modul med prvimi urami namakanja in da dodatek cerija osnovi v koncentraciji 0,01 M mo~no pove~a zdr`ljivost nanosa z zadr`anjem njegove degradacije med namakanjem. Ta projekt je bil usmerjen v dolo~anje kriti~ne koncentracije cerijeve soli, pri kateri se impedan~ni modul hibridnega nanosa mo~no zmanj{a med namakanjem. Klju~ne besede: magnezij, nanos, sol-gel, inhibitor korozije, EIS

METHODOLOGY

Preparation of the materials and coatings

Samples with the dimensions of 40 mm × 20 mm × 6 mm were obtained by making cuttings from a cast Elektron 21 (El21) alloy. The chemical composition of this alloy is shown in Table 1. The samples were first mechanically polished with abrasive papers with a grit of up to grade 4000, then with alumina paste (3 μm and 1 μm) and finally they were rinsed with ethanol and dried under a flux of cold air. The sols were produced by mixing the starting precursors consisting of tetraethyl-orthosilicate (TEOS) and

INTRODUCTION

With a density equivalent to 2/3 of that of aluminium, magnesium and its alloys are interesting weight-saving materials for the automotive and aeronautics industries. However, compared to steel and aluminium alloys, magnesium alloys have a very low corrosion resistance. In order to prevent this problem, various surface treatments and coatings have been developed with different techniques over the last few years 1,2 . However, most of these processes make use of chromium (Cr VI) compounds, nowadays forbidden by international regulations since these are classified as carcinogen, mutagenic and reprotoxic compounds. The sol-gel route is an efficient method to produce "green" coatings and their anti-corrosion performances have been proven successful on steel and aluminium alloys. [3][4][5] This project aims to evaluate the anti-corrosive properties of a hybrid coating obtained via the sol-gel route and deposited on a cast Elektron 21 magnesium alloy (El21) and, secondly, to identify its mechanisms.

3-(trimethoxysilyl)propyl-methacrylate (MAP), deionized water and ethanol with a molar ratio of 11 : 1 : 60 : 80, under constant stirring and at room temperature. In order to adjust the pH of the sol to 4, nitric acid (HNO 3 ) was drop-added to the mixture when required. The production of cerium-doped sols was performed by adding cerium nitrate (Ce(NO 3 ) 3 • 6H 2 O) at four different concentrations: (0.005, 0.01, 0.05 and 0.1) mol L -1 . The addition of the cerium salt was carried out by previously dissolving this compound in the corresponding water volume of the formulation. After maturing for 24 h, the sols were deposited on the magnesium El21 substrates using the dip-coating technique, at a controlled withdrawal speed of 200 mm min -1 . They were then dried at 60 °C for 20 min.

Characterization techniques

The microstructures of the El21 substrate and the sol-gel coatings were analyzed with scanning electron microscopy (SEM) using a JEOL JSM-6510LV microscope, at an operating voltage of 20 kV. Electrochemical tests of the open circuit potential (E ocp ) and electrochemical impedance spectroscopy (EIS) were performed in 0.05 mol L -1 of a NaCl corrosive solution at room temperature, using a Bio-Logic SP-150 potentiostat. The electrochemical cell consisted of a one-chamber threeelectrode cell, the working electrode having an exposed area of 2 cm², delimited with an insulating tape. The reference and auxiliary electrodes included a saturated calomel electrode (SCE) and a platinum-foil electrode, respectively. The EIS spectra were drawn using the potentiostatic mode and a frequency ranging from 100 mHz and 10 mHz, with an applied voltage oscillation of 10 mV vs. OCP. For each test, three samples were analyzed in order to check the reproducibility of the tests.

RESULTS AND DISCUSSION

Morphology of the hybrid coatings

The microstructure of a hybrid sol-gel coating was first observed with SEM (Figure 1). The surface of the coating (a) shows a homogeneous surface, with the presence of some cracks and defects spotted in the neodymium-rich zones. The origin of these defects may be attributed to the formation of a galvanic couple between this intermetallic phase and the alpha phase of magnesium 2 , and to the internal stresses of the hybrid coating that lead to fracture. A cross-sectional observation of the substrate in the BSE (back-scattered electrons) mode (b) allows the thickness of the hybrid coating to be measured at around 1 μm.

Electrochemical characteristics of the hybrid coatings

Firstly, the open-circuit potential of the samples was recorded during the immersion of the substrates in a corrosive solution containing 0.05 mol L -1 of NaCl (Figure 2). All the samples exhibit a similar behavior, except the hybrid coatings that were non-doped or doped with 0.005 mol L -1 of Ce. These show a stabilized potential from the beginning of the immersion, at around -1.63 V and -1.65 V, respectively, attributed to the insulating effect of the protective coating. In contrast, the sol-gel coatings containing 0.01 mol L -1 or higher concentrations, present a behavior similar to that of the bare El21 alloy. This is related to the growth of the passive layer and corrosion products at the surface of the substrates, due to the reaction with the electrolyte. 6 Secondly, the hybrid coatings were tested with EIS after 1 h of the immersion in the corrosive solution, right after the OCP recording. The Bode plots of the EIS spectra obtained for different protective systems are shown in Figure 3. It is worth noting that only the hybrid coating doped with 0.1 mol L -1 presents a behavior similar to that of the bare El21 substrate. The Bode phase-angle diagram (a) shows that the last group presents two time constants, at the high and low frequencies (10 kHz and 1 Hz, respectively). The first is normally attributed to the capacitive response of a hybrid sol-gel coating, which indicates that this film has a physical barrier effect. 7 The second time constant is attributed to the presence of a porous layer in the corrosion products. 8 It is important to observe that the phase angle of a hybrid coating is the highest when it is doped with 0.01 mol L -1 of cerium. On the other hand, the impedance modulus (b) obtained at a low frequency (10 mHz) is typically assigned to the resistance of the electrochemical system, and so to its corrosion resistance. 9 The cerium concentrations lower than 0.1 mol L -1 , especially 0.01 mol L -1 , exhibit higher impedance values. Figure 4 presents the evolution of a hybrid coating doped with 0.01 mol L -1 during its immersion in the corrosive solution. The time constant attributed to the coating gradually disappears with the immersion time (Figure 4a, the time constant at 10 kHz), simultaneously with the shift of the time constant at a low frequency from 1 Hz to 20 Hz. After 48 h of immersion, both curves, representing the bare El21 substrate and the hybrid coating, are superimposed, meaning that the coating lost its protective properties. However, the corrosion resistance of the hybrid film, depicted by the impedance modulus at a low frequency (Figure 4b) shows a progressive decrease with the time. Here, it is proven that the hybrid coating doped with 0.01 mol L -1 of cerium has the best corrosion resistance. Higher or lower concentrations of cerium decrease the coating resistance and capacitance, which leads to a rapid loss of the protective properties. 10 The inhibiting properties of this element are strongly determined by the ion concentration inside the hybrid film, 11 showing the existence of the optimum cerium concentration. This is due to the formation of insoluble compounds such as CeO 2 and Ce(OH) 3 12 that temporally block the passage of the corrosive species through the hybrid coating to the metallic substrate.

CONCLUSION

A hybrid sol-gel coating offers a slight protection to a magnesium substrate during the first hours of an immersion. Moreover, an addition of cerium inside the coating with the optimum concentration of 0.01 mol L -1 leads to the increase of its anti-corrosion properties. This is due to the corrosion-inhibiting effect of cerium ions that allows a formation of insoluble compounds, enhancing the resistance of the hybrid coating to corrosive species.

Figure 2 :

 2 Figure 2: Evolution of the open-circuit potential (E ocp ) of the hybrid coatings doped with cerium of different concentrations, during the immersion in a corrosive solution of 0.05 mol L -1 Slika 2: Razvoj potenciala odprtega kroga (E ocp ) hibridnega nanosa, dopiranega s cerijem v razli~nih koncentracijah, med namakanjem v korozijski raztopini 0,05 mol L -1

Figure 1 :

 1 Figure 1: SEM images of the hybrid sol-gel coating on the El21 magnesium alloy: a) surface of the coating, b) cross-section of the substrate in the BSE mode Slika 1: SEM-posnetka sol-gel hibridnega nanosa na magnezijevi zlitini El21: a) povr{ina nanosa, b) pre~ni prerez podlage v BSE-na~inu

Figure 4 :

 4 Figure 4: Bode plots of the results obtained with EIS for the hybrid coating doped with 0.01 mol L -1 of cerium, during immersion in the corrosive solution (0.05 mol L -1 ): a) phase-angle diagram, b) impedance modulus Slika 4: Bodejev diagram rezultatov, dobljenih z EIS hibridnega nanosa, dopiranega z 0,01 mol L -1 cerija, med namakanjem v korozijski raztopini (0,05 mol L -1 ): a) fazni kotni diagram, b) impedan~ni modul

Figure 3 :

 3 Figure 3: EIS spectra of the El21 magnesium alloy covered with the hybrid coating, doped with different concentrations of cerium. Results obtained after 1 hour of immersion in 0.05 mol L -1 of NaCl: a) phase-angle diagram, b) impedance modulus. Slika 3: EIS-spektri magnezijeve zlitine El21 s hibridnim nanosom, dopiranim z razli~nimi koncentracijami cerija. Rezultati, dobljeni po 1 h namakanja v 0,05 mol L -1 NaCl: a) fazni kotni diagram, b) impe-dan~ni modul.

Table 1 :

 1 Chemical composition of the Elektron 21 cast alloy in mass fractions, w/%

	Tabela 1: Kemijska sestava livne zlitine Elektron 21 v masnih dele`ih,
	w/%
	Element	Nd Gd Zr Zn	Other rare earths	Mg
	w/%	3.1 1.7	1	0.5	< 0.4	Balance

Acknowledgments

The FDA and the OSEO are gratefully acknowledged for the funding provided for this project. The authors would like to thank the partners in the FUI CARAIBE project and its coordinator: the SAFRAN group.