Stabilization and Robustness Analysis for a Chemostat Model with ...
Proceedings of the 46th IEEE Conference on Decision and Control New Orleans, LA, USA, Dec. 12-14, 2007
ThC01.6
Stabilization and Robustness Analysis for a Chemostat Model with Two Species and Monod Growth Rates via a Lyapunov Approach Fr´ed´eric Mazenc
Michael Malisoff
Abstract— We study a two species chemostat model with one limiting substrate. We design feedback controllers so that an equilibrium with arbitrary prescribed species concentrations becomes globally asymptotically stable. The novelty of our work is that we assume that only a linear combination of the species concentrations is available for measurement, combined with our use of Lyapunov function methods to generate stable coexistence and quantify the effects of disturbances. Our closed-loop error dynamics are integral input-to-state stable with respect to small perturbations of the controllers. Hence, no matter what initial positive levels for the species concentrations and substrate are selected, the long term species and substrate levels remain close to the equilibrium, even when there are small unexpected changes in the dilution rate and input nutrient concentration. This is a highly desirable robustness property, because the dilution rate is prone to actuator errors. We illustrate our approach using a numerical example. Index Terms—Bioreactors, chemostats, global stabilization, uncertain states, integral input-to-state stability
I. I NTRODUCTION This note continues our search (begun in [14], [15]) for ways to stabilize prescribed equilibrium behaviors in chemostats with one limiting substrate. The chemostat is known in engineering as the continuously stirred tank reactor. See for example [6], [18] for the fundamental role of chemostats in bioengineering, ecology, and population biology; and see Sections II-IV below for the relevant definitions and discussion of the known results. The novelty of the present work is in our design of controllers that stabilize prescribed equilibria (with arbitrary species concentrations) when only an arbitrary linear combination of the species concentrations is available for measurement, combined with our use of a strict Lyapunov function to quantify the effects of actuator errors. Hence, our work leads to species coexistence and persistence (via new strict Lyapunov function and feedback designs) in important cases where there is uncertainty regarding both the current state and the feedback controller. The basic model for a (well mixed) chemostat with two competing species and one limiting substrate is ( 1 (s) X1 − µ2Y(s) X2 , s˙ = D(sin − s) − µY 1 2 (1) X˙ i = [µi (s) − D]Xi , i = 1, 2, Supported by NSF Grant 0424011. Mazenc is with Projet MERE INRIA-INRA, UMR Analyse des Syst`emes et Biom´etrie INRA, 2, pl. Viala, 34060 Montpellier, France (email: [email protected]). Malisoff is with the Department of Mathematics, Louisiana State University, Baton Rouge, LA 70803-4918 USA (email: [email protected]). J. Harmand is with Institut National de la Recherche Agronomique (INRA), UR050, Laboratoire de Biotechnologie de l’Environnement, Narbonne, F-11100, France (email: [email protected]).
1-4244-1498-9/07/$25.00 ©2007 IEEE.
J´erˆome Harmand
evolving on X := (0, ∞)3 ; see [18, Chapter 1] for the derivation of the model. Here and in the sequel, dots indicate time derivatives; s(t) is the concentration of the substrate; X1 (t) and X2 (t) are the concentrations of the two species of organism; the dilution rate D(·) and the input nutrient concentration sin (·) are nonnegative controllers we will specify; µ1 and µ2 are given uptake functions; and Y1 and Y2 are positive constant yield coefficients (but see Section III below for more general models with uncertainties). We assume throughout that the µi s have the Monod form µi (s) =
Ki s Li + s
(2)
where Ki and Li are positive constants we indicate below. Since D is the ratio of the volumetric flow rate (with units of volume over time) to the constant reactor volume, it is proportional to the speed of the pump that supplies the reactor with fresh medium containing the nutrient. The well known “competitive exclusion principle” implies that in classical chemostats with one limiting substrate and constant D and sin , at most one species can survive [11]. This is at odds with the observation that in real ecological systems, it is common for many species to coexist in equilibrium, even if there is only one limiting nutrient. This paradox has motivated a great deal of ongoing research. Much of the research is devoted to choosing D and sin to be time-varying or state-dependent controllers that force a stable coexistence in which multiple species survive [6], [8]. Another important issue is that the species concentrations may not be available for measurement, and there may also be actuator errors caused e.g. by unexpected variability in the speed of the pump supplying the reactor with fresh nutrient. This note addresses all of these issues by designing sin and a controller D(·), depending only on a linear combination Y = X1 + AX2
(3)
of the species concentrations (where A is a given positive constant), that globally stabilize a prescribed equilibrium (s? , x1? , x2? ) ∈ X (with arbitrary positive prescribed species concentrations x1? and x2? , and with s? determined by the µi s using Assumption 1 below). Therefore, no matter what initial values are selected, the chemostat trajectories in closed loop with our feedbacks all converge to (s? , x1? , x2? ). Another important feature of our work is that our use of an explicit strict Lyapunov function construction makes it possible to quantify the effects of disturbances using the input-to-state stability (ISS) paradigm that was introduced by Sontag in his seminal paper [19]. See Section III for the
3933
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 relevant definitions; and see [3], [21] for the fundamental role of ISS in nonlinear control theory and engineering applications. In fact, to our knowledge, our work provides the first ISS analysis for multiple species chemostats whose input nutrient concentrations, outputs, dilution rates, and uptake functions are all perturbed by small noise. In Section II, we state our assumptions and stability result for (1). Section III explains how the stability provided by our feedback is robust to small actuator errors. We compare our work with the known results in Section IV. In Section V, we prove our result for (1). We prove our robustness result in Section VI and we illustrate our methods in Section VII. We close in Section VIII with suggestions for further research. II. A SSUMPTIONS
AND
S TABILITY T HEOREM
Let SAT denote the set of all continuously differentiable functions σ : R → R satisfying (I) r 7→ rσ(r) is positive definite (i.e., zero at zero and strictly positive at all other points) and (II) |σ(r)| ≤ min{|r|, 1} for all r ∈ R (which we refer to as saturations). For example, σ ∈ SAT if r . (4) σ(r) = √ 1 + r2 Given µi s as in (2) (where the Ki s and Li s are given positive constants), set χ(s) = µ2 (s) − µ1 (s). We always assume: Assumption 1: There is a constant s? > 0 such that (i) χ(s? ) = 0, (ii) χ(s) < 0 when 0 < s < s? , and (iii) χ(s) > 0 when s > s? . In particular, µ1 (s? ) = µ2 (s? ). Simple calculations and our Assumption 1 readily yield Li µi (s) (s − r) ∀s, r > 0 Li + r s for i = 1, 2 and (using the formula for χ) µi (s) − µi (r) =
s? =
L1 K2 −L2 K1 K1 −K2
, µ1 (s? ) =
L1 K2 −L2 K1 L1 −L2
,
(5)
again evolving on X = (0, ∞)3 (where sin and D are controllers to be chosen) and the output y = Y /Y1 satisfying
(8)
We always assume that a 6= 1. Given arbitrary prescribed x1? > 0 and x2? > 0, we set sin? = s? + x1? + x2? ; x˜i = xi − xi? for i = 1, 2; s˜ = s − s? ; ξ˜i = ln(˜ xi + xi? ) − ln(xi? ) = ln(xi /xi? ) for ˜ = ln(˜ i = 1, 2; and Σ s + s? ) − ln(s? ) = ln(s/s? ).
β : [0, ∞) × [0, ∞) → [0, ∞) for which (1) β(·, t) ∈ K∞ for each t ≥ 0 and (2) β(r, ·) is non-increasing and β(r, t) → 0 as t → +∞ for each r ≥ 0. Let δB2 denote the closed ball of radius δ > 0 centered at 0 ∈ R2 . We use the constants ω1 = ω3 = and
, ω2 =
µ1 (s? ) 4s? La − L1 +1 1 2
,
µ1 (s? )s? , 20 (a−1)sin? +s2? La − L1 +1 1 2 0.72µ1 (s? )s2? L1 − L1 1 2 h i ω4 = 2 (a−1)sin? +s2? La − L1 +1
(10) .
2
The following yields (s(t), x1 (t), x2 (t)) → (s? , x1? , x2? ) as t → +∞ for all trajectories (s(t), x1 (t), x2 (t)) of (7) (i.e., for all initial values in X ) when the indicated controllers are used, hence globally asymptotically stable coexistence: Theorem 1: Let ε ∈ (0, min{ω1 , ω2 , ω3 , ω4 }], σ ∈ SAT , and x1? , x2? > 0 be given. Let s? > 0 satisfy the requirements of Assumption 1. Consider (7) with output (8) with a 6= 1. Then (s? , x1? , x2? ) is a globally asymptotically stable equilibrium of (7) when sin = sin? := s? + x1? + x2? and D(y) = µ1 (s? ) − ε(a − 1)σ (y − x1? − ax2? ). More precisely, there exists β ∈ KL such that for all trajectories (s, x1 , x2 )(t) of (7) with sin and D (for all initial values ˜ ξ˜1 , ξ˜2 )(t) from (9) in X ), the transformed error vector (Σ, ˜ ˜ ˜ ˜ ˜ ˜ satisfies |(Σ, ξ1 , ξ2 )(t)| ≤ β(|(Σ, ξ1 , ξ2 )(0)|, t) for all t ≥ 0. Remark 1: The controller D from Theorem 1 consists of µ1 (s? ) plus a term that can be made arbitrarily small by picking a sufficiently small ε > 0, and it yields coexistence because x1? > 0 and x2? > 0. However, we cannot pick ε = 0 since that would make D and sin constants, and so yield competitive exclusion [6]. III. ROBUSTNESS
TO
D ISTURBANCES
Consider the perturbed chemostat error dynamics s˜˙ = [D(y)+d2 ](sin +d1 −s) − µ1 (s)x1 − µ2 (s)x2 x ˜˙ i = [µi (s) − D(y) − d2 ]xi , i = 1, 2
(11)
with output (8) for a 6= 1, where s˜ = s−s?; x ˜i = xi −xi? for i = 1, 2; sin and D are from Theorem 1; σ is from (4); the µi s are as before (so Assumption 1 holds); and the unknown measurable disturbance d(t) = (d1 (t), d2 (t)) is bounded in sup norm |·|∞ by a constant that we specify below. For each ρ > 0, let D2 (ρ) denote the set of all measurable essentially bounded d : [0, ∞) → ρB2 . Fix arbitrary x1? > 0 and x2? > 0, and let ε > 0 be as in the statement of Theorem 1. To simplify the notation, we set 2? L2 +s? , β¯ = xx1? , L2 ˜ ˜ , A0 (˜ s) = s˜ − s? ln 1 + ss˜? = s? eΣ − 1 − Σ ˜ and Ai (ξ˜i ) = eξi − 1 − ξ˜i = xx˜i?i − ln 1 + xx˜i?i
α ¯=
(9)
For each trajectory (s, x1 , x2 )(t) of (7), we call the vector ˜ ξ˜1 , ξ˜2 )(t) as defined by (9) the (corresponding trans(Σ, formed) error vector. Let K∞ denote the set of all continuous ρ : [0, ∞) → [0, ∞) for which (i) ρ(0) = 0 and (ii) ρ is strictly increasing and unbounded. If ρ1 , ρ2 ∈ K∞ , then ρ−1 1 , ρ1 ◦ ρ2 ∈ K∞ . Let KL denote the set of all continuous
µ1 (s? ) 20|a−1|
1
(6)
and therefore K2 > K1 , K1 L2 > K2 L1 , and L2 > L1 . We transform (1) with the output Y from (3), using x1 = X1 /Y1 and x2 = X2 /Y2 . Doing so produces the system s˙ = D(sin − s) − µ1 (s)x1 − µ2 (s)x2 , (7) x˙ i = [µi (s) − D]xi , i = 1, 2
y = x1 + ax2 , where a = AY2 /Y1 .
ThC01.6
L1 +s? L1
(12) (13)
for i = 1 and 2, where we used (9). We show that (11) is integral input-to-state stable (iISS) ˜ ξ˜1 , ξ˜2 ) (from (9) above) for distur[20] in the variable (Σ, ¯ ¯ is defined in (15), in the sense bances d ∈ D2 (∆) when ∆
3934
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 that there exist α ∈ K∞ and β ∈ KL and a constant K > 0 ¯ and all corresponding closed-loop so that for all d ∈ D2 (∆) trajectories (˜ s(t), x ˜1 (t), x ˜2 (t)) of (11) for our controllers D and sin (for all initial values), the corresponding transformed ˜ error vector (Σ(t), ξ˜1 (t), ξ˜2 (t)) satisfies ˜ ξ˜1 , ξ˜2 )(t)|) α(|(Σ, R ˜ ξ˜1 , ξ˜2 )(0)|, t) + K t |d(r)|dr ∀t ≥ 0 . ≤ β(|(Σ, 0
(iISS)
(This is a stronger condition than the usual iISS condition from [20] since it involves the integral of |d(r)| instead of the usual integral of a K∞ function of |d(r)|.) We also show that (11) is input-to-state stable (ISS) in ˜ ξ˜1 , ξ˜2 ) when d ∈ D2 (∆ ¯ [ ) and ∆ ¯ [ is in (41), meaning (Σ, ¯ [) there are β ∈ KL and γ ∈ K∞ so that for all d ∈ D2 (∆ and all closed-loop trajectories (˜ s, x ˜1 , x ˜2 )(t) of (11) for our ˜ ξ˜1 , ξ˜2 )(t) satisfies D and sin (for all initial values), (Σ, ˜ ξ˜1 , ξ˜2 )(t)| |(Σ, ˜ ξ˜1 , ξ˜2 )(0)|, t) + γ(|d|∞ ) ∀t ≥ 0. ≤ β(|(Σ,
(ISS)
See [3], [21] for the fundamental roles of iISS and ISS in nonlinear control and engineering. The ISS and iISS estimates give (s, x1 , x2 )(t) → (s? , x1? , x2? ) as t → +∞ for all trajectories (s, x1 , x2 )(t) of (11) when the indicated controllers D and sin are used and d = 0, and they also quantify the effects of actuator errors d. Also, (ISS) guarantees persistence, since if for example we had x2 (t) → ˜ 0 as t → t− ? for some t? ∈ (0, +∞], then ξ2 (t) = ln(x2 (t))− ln(x2? ) → −∞, which would violate (ISS). Similarly, (iISS) implies persistence when d is integrable. To specify the bounds on d, we introduce the constants C1 = s10? L11 − L12 ε 2 L1 +s? L2 +s? ε2 a − + 1 C2 = 19µ20 4C1 L1 L2 1 (s? ) (14) n o ?) κ1 = min 1, µ1 (s +x +x +µ (s )+ε|a − 1| 1? 2? 1 ? 5 n o κ1 1 κ2 = 2s2 + 2C2 + sin? κ1 C2 +max α, ¯ xx1? β¯ . 2? in?
See Section VII below for examples where the bounds on d are computed explicitly. Our iISS result is as follows (but see Remark 3 for ISS results under appropriate conditions): Theorem 2: Under the preceding assumptions with ¯ = min 1, 2µ1 (s? )s? , 1 ∆ (15) 5(3κ1 + 2s? ) 2κ2
˜ ξ˜1 , ξ˜2 ) for disturbances d ∈ D2 (∆). ¯ (11) is iISS in (Σ, Remark 2: A much more general scenario is where there are measurement errors in D(·) and uncertain uptake functions, as well as actuator errors. This gives the error dynamics s˜˙ x ˜˙ 1 x ˜˙ 2
= [D(y + r3 ) + u2 ](sin + u1 − s)
−(1 + r1 )µ1 (s)x1 − (1 + r2 )µ2 (s)x2 ,
= [(1 + r1 )µ1 (s) − D(y + r3 ) − u2 ]x1 , = [(1 + r2 )µ2 (s) − D(y + r3 ) − u2 ]x2 .
(16)
ThC01.6 d = (r, u) for suitable small values of the disturbance. Due to space constraints, we leave this extension to the reader. IV. C OMPARISON
WITH
K NOWN R ESULTS
The importance of coexistence and stability in chemostats has motivated a great deal of significant research at the interface of biology and mathematical control theory [2], [4], [5], [6], [23]. Much of this work uses monotonicity or reduction methods to establish stability. By contrast, Lyapunov functions have not frequently been used to study stability in chemostats with multiple species, and where they are used they are often weak (i.e., nonstrict) Lyapunov functions which do not lend themselves to robustness analysis. An exception is the famous one from [12]; see also Theorem 4.1 on [18, p.35] and [10]. Strict Lyapunov functions for two species chemostats were explicitly constructed in [15], but no robustness to disturbances was established there. However, [15] provides simple linear feedback stabilizers that yield local asymptotic stability of a prescribed periodic trajectory for the closed-loop dynamics (i.e. tracking). See [13], [23] where weak Lyapunov functions are used in conjunction with variants of the LaSalle Invariance Principle. This raises the important question of whether strict Lyapunov functions can be explicitly constructed for multispecies chemostats that are globally stabilized through static output feedbacks, and if so, whether the Lyapunov functions can be used to quantify the effects of uncertainty. For chemostats with only one species that are made oscillating through a suitable choice of D(·), this type of problem was solved in [14]. There, a Lyapunov function is constructed and yields ISS relative to small actuator uncertainty. Here we extend these results to two-species models with disturbance perturbations, or where only a linear combination of the species concentrations is known. Our work owes a great deal to [6], [7], which stabilize chemostats in which only the sum of the species concentrations can be measured, using an appropriate dilution rate. However, [6] and [7] do not include our work since [6] does not rely on a Lyapunov approach and the Lyapunov functions from [7] are nonstrict and so do not lend themselves to ISS. See [9] for results where only the substrate level (multiplied by an a priori bounded error) is available for measurement. V. P ROOF
OF
T HEOREM 1
In all of what follows, all (in)equalities should be understood to hold globally unless otherwise indicated, we omit the arguments of our functions whenever this would not lead to confusion, sin = s? + x1? + x2? , and D? = µ1 (s? ). Since s? is constant and µ1 (s? ) = µ2 (s? ), the change of feedback D = D? + v, (5), and (9) transform (7) into ( s˜˙ = (D? + v)(sin − s) − µ1 (s)x1 − µ2 (s)x2 , (17) ˙˜ µi (s) i ξi = LiL+s ˜ − v, i = 1, 2 . s s ? Step 1: Construction of a weak Lyapunov function
A variant of our proof of Theorem 2 shows that (16) is ISS in ˜ ξ˜1 , ξ˜2 ) with respect to the disturbance vector the variable (Σ,
3935
We first show that in terms of (12) and (13), 1 ¯ 2 (ξ˜2 ) V (˜ s, ξ˜1 , ξ˜2 ) = A0 (˜ s) + α ¯ A1 (ξ˜1 ) + βA x1?
(18)
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 is a weak Lyapunov function for (17) in the sense that 1 s˜ x˜1 x ˜2 V˙ = −H(˜ s)˜ s2 + (sin −s)v − α ¯ + β¯ v (19) x1? s x1? x2?
holds along all trajectories of (17), where H as defined by h 1 (s? ) H(˜ s) = sx11? µ1 (s? ) + µ1 (s)−µ x1? s˜ i (20) µ2 (s)−µ2 (s? ) x , s ˜ = 6 0 + 2? s˜
is everywhere positive. To this end, first notice that by (9), s˜ x ˜1 ˜˙ x˜2 ˜˙ A˙ 0 = s˜˙ , A˙ 1 = ξ1 , and A˙ 2 = ξ2 . s x1? x2?
(21)
Hence, along the trajectories of (17), we easily get h i 1 2 V˙ = x11? ss˜ (D? + v)(sin − s) − α ¯ xx˜1? + β¯ xx˜2? v h i µ (s) L1 ˜1 x11? 1s s˜ + −x1 + α ¯ L1 +s? x h i L2 1? + −x2 + β¯ xx2? ˜2 x11? µ2s(s) s˜. L2 +s? x =
1 x1?
[D? (sin − s) − xh1? µ1 (s) − x2?iµ2 (s)] ss˜ 1 2 ¯ xx˜1? + β¯ xx˜2? v. + x11? ss˜ (sin − s)v − α
Step 2: Construction of a strict Lyapunov function Set y˜ = x ˜1 + a˜ x2 , Z˜ = s˜ + x ˜1 + x ˜2 , Ψ(r) = 21 r2 , and
2s?
L1
−L
(23)
where V is from (18). We show that V2 is a strict Lyapunov function for (17) in the sense that its time derivative along the trajectories of (17) with v = −ε(a − 1)σ(˜ y ) satisfies 4s
10
1 L1
i − L12 ε˜ yσ (˜ y)+ Z˜ 2
(24)
and x˜i and ξ˜i are related by (9). Then W2 is positive definite (since a 6= 1, L2 > L1 , and r 7→ rσ(r) is positive definite). To this end, first notice that by our choices of ω1 , ε, and D, we get D(y) ≥ 19 20 D? everywhere. Therefore, (7) gives Z˜˙ = D(y)[sin − s − x1 − x2 ] = −D(y)Z˜ d ˜ = −D(y)Z˜ 2 ≤ − 19 D? Z˜ 2 . hence dt Ψ(Z) 20
(25)
Moreover, our assumption that a 6= 1 gives
a(Z˜ − s˜) − y˜ y˜ − Z˜ + s˜ x ˜1 = and x ˜2 = . a−1 a−1
−L
2
Applying the general relation ab ≤ 41 a2 + b2 , we get n o (a−1)sin +s2 a − 1 s˜ p ? L1 L2 √ µ1 (s? ) √ εσ(˜ y ) s µ1 (s? )s µ1 (s? ) 2 ˜ 4s s
+
h
(a−1)sin +s2?
a L1
µ1 (s? )s
− L1
2
i2
ε2
+
h
(a−1)sin +s2?
a L1
− L1
2
µ1 (s? )s
i2
ε2
(28)
2
σ (˜ y)
for all s > 0. Combining (27) and (28), it follows that h i (s? ) 2 1 1 ? V˙ 2 ≤ − µ12s s˜ − 9s y σ (˜ y ) − Z˜ 2 10 L1 − L2 ε˜
(29)
σ 2 (˜ y) .
We next distinguish between two cases. Case A: If s ≥ then (29) and our choices of ε and ω4 give h i (s? ) 2 1 1 ? V˙ 2 ≤ − µ12s s˜ − 9s − y σ (˜ y ) − Z˜ 2 10 L1 L2 ε˜
9 10 s? ,
+
h i2 10 (a−1)sin +s2? La − L1 ε2 1
2
9µ1 (s? )s?
(s? ) 2 ≤ − µ12s s˜ −
s? 10
h
1 L1
−
1 L2
i
σ 2 (˜ y)
(30)
ε˜ y σ (˜ y ) − Z˜ 2
9 because σ 2 (˜ y ) ≤ y˜σ(˜ y ) everywhere. Case B: If s ≤ 10 s? , s? then |˜ s| ≥ 10 , so (29) and our choices of ε and ω3 give h i (s? ) 2 1 (s? ) 2 V˙ 2 ≤ − µ14s s˜ − µ400s s? − s2? L11 − L12 ε˜ y σ (˜ y)
−Z˜ 2 +
h
(a−1)sin +s2?
(s? ) 2 ≤ − µ14s s˜ −
(26)
Substituting (12) and (26) into (19), we easily get h i s? s? 1 V˙ = −H(˜ s)˜ s2 + x11? 1s sin + a−1 a − s˜v + L L2 h i h 1 i s? L1 +s? L2 +s? ˜ 1 1 1 − y ˜ v − a − Zv. x1? (a−1) L1 L2 x1? (a−1) L1 L2
L1
It follows from (25) and our formula (23) for V2 that h i 1 1 ? − y σ (˜ y ) − Z˜ 2 V˙ 2 ≤ − 3µ14s(s? ) s˜2 − 9s 10 L L2 ε˜ h 1 i (27) s˜ − (a − 1)sin + s2? La1 − L12 εσ (˜ y ) . s
≤
2
V˙ 2 ≤ −W2 , where h W2 (˜ s, ξ˜1 , ξ˜2 ) = µ1 (s? ) s˜2 + s?
By our choices of ε and ω2 , the fact that σ 2 (˜ y) ≤ y˜σ(˜ y) δ 2 5 2 everywhere, and the general relation ab ≤ a + b for 2δ 10 √ b = εσ(˜ y ) and δ = s? ( L11 − L12 ) > 0, we then easily get h i y σ (˜ y) x1? V˙ ≤ − 3µ14s(s? ) s˜2 − s? L11 − L12 ε˜ h i s˜ − (a − 1)sin + s2? La1 − L12 y) s εσ (˜ n√ h i o √ ? ? ε L1L+s a − L2L+s Z˜ { εσ (˜ y)} + 1 h 2 i 1 1 ? ≤ − 3µ14s(s? ) s˜2 − 9s y σ (˜ y) 10 L − L2 ε˜ h 1 i s˜ − (a − 1)sin + s2? La1 − L12 y) s εσ (˜ h i2 ? ? + 15ε 1 L1L+s a − L2L+s Z˜ . 1 2
(22)
The formulas sin = s? + x1? + x2? and D? = µ1 (s? ) and the fact that µ1 (s? ) = µ2 (s? ) therefore easily give (19).
V2 (˜ s, ξ˜1 , ξ˜2 ) = x1? V (˜ s, ξ˜1 , ξ˜2 ) + 2 20 ˜ 5ε L1 +s? a − L2 +s? + 1 Ψ(Z) 19D? 1 L1 L2 1
Taking v = −ε(a − 1)σ(˜ y ) and recalling (20) gives h i (a−1)sin s˜+s? s? s? ?) 2 s ˜ − + a− s˜εσ (˜ y) x1? V˙ ≤ − µ1 (s s s s L1 L h i h i 2 ? ? ˜ (˜ −s? L11 − L12 ε˜ yσ (˜ y ) + L1L+s a − L2L+s Zεσ y) . 1 2
2s?
Since −xi + x ˜i = −xi? for i = 1, 2, our choices (12) give V˙
ThC01.6
a L1
− L1
2
i2
ε2
h µ1 (s? )s i s? 1 1 y σ (˜ y) 2 L1 − L2 ε˜
(31)
− Z˜ 2
where we used |σ(˜ y )| ≤ 1. In either case, we get (24). Since W2 is positive definite, V2 is the desired strict Lyapunov function. However, V2 becomes unbounded as s˜ → −s? so extra care is needed to obtain β ∈ KL from the statement of the theorem; see [16] for details. This proves Theorem 1.
3936
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007
ThC01.6
VI. P ROOF OF T HEOREM 2 ¯ and let ∆ denote its sup norm, so We fix d ∈ D2 (∆) ¯ We use the notation from the proof of Theorem 1. ∆ ≤ ∆. We take sin = sin? := s? + x1? + x2? throughout and ¯ 1? /x2? sin? C3 = 2κ1 + 2 2C2 κ1 + 2 max α, ¯ βx +8C2 s2in? and, in terms of the relations (9), (32) (s? ) 2 W4 (˜ s, ξ˜1 , ξ˜2 ) = µ120s s˜ + C1 y˜σ (˜ y ) + 21 Z˜ 2 ,
¯ 1? /x2? |˜ ¯ 1? /x2? }Z˜ 2 /sin? . ∆{¯ α|˜ x1 |+ βx x2 |} ≤ ∆ max{¯ α, βx Substituting into (38) and recalling (14),
where σ is from (4). Rearranging terms in (11) gives s˜˙ = D(y)[sin − s] − µ1 (s)x1 − µ2 (s)x2 +g1 (˜ s, y, d), ˜˙ ξi = µi (s) − D(y) − d2 , i = 1, 2
¯ ≤ 1/(2κ2 ), we get V˙ 2 ≤ −W4 (˜ Since ∆ s, ξ˜1 , ξ˜2 ); see (32). 2nd case: |˜ s + x˜1 + x˜2 | ≤ 2sin? . One easily verifies [16]: Lemma 3: If |˜ s+x ˜1 + x ˜2 | ≤ 2sin? , then |˜ s|+|˜ x1 |+|˜ x2 | ≤ 4sin? . Using Lemma 3 and (38), we deduce in the second case that V˙ 2 ≤ −W3 (˜ s, ξ˜1 , ξ˜2 ) + C3 ∆ ≤ −W4 (˜ s, ξ˜1 , ξ˜2 ) + C3 ∆, so
(33)
with g1 (˜ s, y, d) = d2 (−˜ s + x1? + x2? ) + [D(y) + d2 ]d1 . The time derivative of V2 along the trajectories of (33) satisfies: i h ∂V2 ∂V2 2 V˙ 2 ≤ −W2 + ∂V g (˜ s , y, d) − d + , (34) 1 2 ˜ ˜ ∂ s˜ ∂ξ ∂ξ 1
2
by (24). In terms of the constants from (14) and the relations (9), we have |g1 (˜ s, y, d)| ≤ ∆(κ1 + |˜ s|), (s? ) 2 W2 (˜ s, ξ˜1 , ξ˜2 ) = µ14s s˜ + C1 y˜σ (˜ y) + Z˜ 2 , ∂V2 ˜ ∂V2 = α˜ ˜ 1, = s˜ + C2 Z, ¯ x1 + C2 Zx ∂ s˜
and
s ∂V2 ˜ ∂ ξ2
=
x1? ¯ ˜2 x2? β x
∂ ξ˜1
(35)
˜ 2, + C2 Zx
¯ ≤ min{1, µ1 (s? )/5} (which holds rewhere we used ∆ gardless of the choice of κ1 > 0, by (15)) to get κ1 . Set v˜ = |˜ s| + |˜ x1 | + |˜ x2 |. From (34), we get V˙ 2
2 1? ¯ ≤ −W2 + ∆ s˜s + κ1 ∆ |˜ss| + ∆ xx2? β|˜ x2 | ˜ (|x1 | + |x2 | + κ1 + |˜ + ∆C2 |Z| s|) + ∆¯ α|˜ x1 | |˜ s| x1? ¯ s˜2 ≤ −W2 + ∆ s + κ1 ∆ s + ∆ x2? β|˜ x2 | ˜ (2κ1 + v˜) + ∆¯ + ∆C2 |Z| α|˜ x1 |,
(36)
where we used |x1 | + |x2 | + |˜ s| ≤ v˜ + x1? + x2? ≤ v˜ + κ1 . One easily shows that [16]: Lemma 1: For all s > 0, we have 3˜ s2 |˜ s| ≤2+ . (37) s 2s? s Combined with (36), Lemma 1 gives 2 s˜ 1 V˙ 2 ≤ −W2 (˜ s, ξ˜1 , ξ˜2 ) + ∆ 1 + 3κ 2s? s + 2κ1 ∆ ¯ x2 | . ˜ +∆C2 |Z| (2κ1 + v˜) + ∆¯ α|˜ x1 | + ∆ x1? β|˜
V˙ 2
∆κ1 ≤ −W3 (˜ s, ξ˜1 , ξ˜2 ) + 2s (˜ s+x ˜1 + x˜2 )2 2 h i in? 1 +∆C2 sκin? + 2 (˜ s + x˜1 + x ˜ 2 )2 ¯ 1? /x2? }(˜ + s∆ max{¯ α, βx s + x˜1 + x ˜ 2 )2 in? = −W3 (˜ s, ξ˜1 , ξ˜2 ) + ∆κ2 (˜ s + x˜1 + x ˜ 2 )2 .
V˙ 2 ≤ −W4 (˜ s, ξ˜1 , ξ˜2 ) + C3 ∆
holds in both cases. Define V and the positive definite W by ˜ ξ˜1 , ξ˜2 ) = V2 ([eΣ˜ − 1]s? , ξ˜1 , ξ˜2 ) and V(Σ, ˜ ξ˜1 , ξ˜2 ) = W4 ([eΣ˜ − 1]s? , ξ˜1 , ξ˜2 ). W(Σ,
1? ¯ V˙ 2 ≤ −W3 (˜ s, ξ˜1 , ξ˜2 )+2κ1∆+∆¯ α|˜ x1 |+∆ xx2? β|˜ x2 |
+∆C2 |˜ s+x ˜1 + x˜2 | [2κ1 + v˜]
(38)
(s? ) 2 where W3 (˜ s, ξ˜1 , ξ˜2 ) = µ120s s˜ + C1 y˜σ (˜ y) + Z˜ 2 . We consider two cases, based on the magnitude of Z˜ = s˜ + x ˜1 + x ˜2 : ˜ 1st case: |Z| ≥ 2sin? . One easily checks that [16]: ˜ ≥ 2sin? , then |˜ ˜ Lemma 2: If |Z| s| + |˜ x1 | + |˜ x2 | ≤ 2|Z|. 2 2 We deduce that ∆κ1 ≤ ∆κ1 (˜ s + x˜1 + x˜2 ) /{4sin? }, h i 1 ˜ and ∆C2 [2κ1 + v˜] ≤ ∆C2 sκin? + 2 |Z|,
(40)
By (39), V˙ ≤ −W +C3 ∆. We can also find α1 , α2 ∈ K∞ so ˜ ξ˜1 , ξ˜2 )|) ≤ V(Σ, ˜ ξ˜1 , ξ˜2 ) ≤ α2 (|(Σ, ˜ ξ˜1 , ξ˜2 )|) [16]. that α1 (|(Σ, Hence, a standard argument [1, p.1088] gives our estimate (iISS) with α = α1 and K = 2C3 , proving Theorem 2. Remark 3: One can prove that the system (11) is also ISS ˜ ξ˜1 , ξ˜2 ) when d ∈ D2 (∆[ ), where, in terms of (15), in (Σ, 1 [ ¯ ¯ ∆ = min ∆, mini α ˜ 1 (xi? /2) ; (41) 2C3 α ˜ 1 ∈ K∞ is chosen so that α ˜ 1 (r) ≤ r2 and W4 (˜ s, ξ˜1 , ξ˜2 ) ≥ ˜ α ˜ 1 (|(˜ s, x ˜1 , x˜2 )|) hold everywhere; x ˜i and ξi are again related by (9); and C3 is from (32). For the proof, see [16]. Remark 4: An important special case of Theorem 2 is where only the dilution rate is affected by noise; i.e., d1 = 0 and d2 6= 0. This often occurs when the speed of the pump that supplies the chemostat with fresh nutrient is subjected to small fluctuations [14]. When d1 = 0, our proof of Theorem 2 remains valid if κ1 is replaced by x1? + x2? . Another important case is where the noise term is only on the input nutrient concentration sin , i.e., d2 = 0 and d1 6= 0, in which ¯ is replaced by case Theorem 2 remains true if ∆
x2?
From the formula for W2 in (35), Z˜ = s˜ + x ˜1 + x˜2 , our ¯ and our choice (15) of ∆, ¯ we get assumption that ∆ ≤ ∆,
(39)
¯ = ∆
0.16D?s? ; D? + ε|a − 1|
(42)
see [16] for the proof. We illustrate these special cases below. Remark 5: Our bound (15) is conservative. It can be enlarged by reducing κ2 . For example, arguing as in the proof of Theorem 2 except separately considering the cases |˜ s + x˜1 + x ˜2 | ≥ M sin? and |˜ s + x˜1 + x ˜2 | ≤ M sin? (where M ≥ 2 is any constant), one checks [16] that Theorem 2 still holds if κ2 is replaced by 2κ1 1 2x1? ¯ κ2 = 2 2 + 2C2 + 2κ1 C2 + 2¯ α+ β M sin? M sin? x2? so κ2 ≈ 2C2 for large M , if C3 is appropriately enlarged.
3937
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 VII. E XAMPLE We use the following special case of (7)-(8): 0.05sx1 .052sx2 2 )(sin +d1 −s)− 20+s − 25+s s˙ = (D(y)+d h i .05s x˙ 1 = 20+s − D(y) − d2 x1 h i x˙ 2 = .052s − D(y) − d2 x2 25+s
(43)
with y = x1 + 0.8x2 . We take x1? = 0.05 and x2? = 0.02. Using the notation from Theorem 1, a = 0.8, s? = 105, and ε ∈ (0, .00753]. We take sin = 105.07 and D(y) = .042 + 0.001506σ(y − 0.066) where σ is from (4). Theorem 1 implies that regardless of the initial values, all trajectories of (43) with our controllers D and sin converge to (105, 0.05, 0.02) when d = 0. We next consider the important special case where d1 = 0 but d2 6= 0. In that case, we can take κ1 = x1? + x2? ; see ¯ as in (15) and κ2 as in Remark 5 with Remark 4. Taking ∆ ¯ ≈ 0.20D? . By Theorem 2, the stability M = 25 gives ∆ has the highly desirable property of robustness (in the iISS sense) to the noise d = (0, d2 ), provided |d2 |∞ is no more than (about) 20% of D; in particular, if d2 is also integrable, then we have persistence of both species. If instead the disturbance is only added to sin , then we can use (42) to get robustness to disturbances d1 that are bounded ¯ ≈ 16, or about 15% of sin = 105.07. We simulated by ∆ the corresponding dynamics (43) with d(t) ≡ (1, 0) and the initial value (s, x1 , x2 )(0) = (103, 2, 1) and obtained: Substrate 105.5 105 104.5 104 103.5 1000
2000
Species 1
3000
4000
5000
Species 2
2
1 0.8
1.5
0.6 1 0.4 0.5
0.2
1000
2000
3000
4000
5000
1000
2000
3000
4000
5000
Our simulation illustrates the persistence of (x1 , x2 )(t) and the convergence (s(t), x1 (t), x2 (t)) → (105, 0.05, 0.02), with an overshoot determined by the iISS estimate and the magnitude of d1 , and so validates our theorems. We emphasize that in general, our iISS estimate (iISS) will not guarantee persistence unless we choose d to be integrable. VIII. C ONCLUSIONS We provided an output feedback design for stabilizing equilibria with arbitrary prescribed species concentrations in two species chemostats. The novelty of our work is our construction of a Lyapunov function, which made it possible
ThC01.6 to quantify the effects of uncertainties such as actuator errors, measurement noise, and uncertain growth rates using ISS. It would be of interest to extend our work to chemostats with delays between the instant the micro-organism uptakes the substrate and the “instant” the biomass integrates it, as well as to tracking. See [15] for tracking in two species chemostats where the species concentrations are known. Another desirable extension would be to models with three or more competing species, or more than one limiting substrate. R EFERENCES [1] D. Angeli, E.D. Sontag, and Y. Wang, “A characterization of integral input to state stability,” IEEE Trans. Automat. Control, vol. 45, pp. 1082-1097, 2000. [2] O. Bernard and J.-L. Gouz´e, “Nonlinear qualitative signal processing for biological systems: Application to the algal growth in bioreactors,” Mathematical Biosciences, vol. 157, pp. 357–372, 1999. [3] M. Chaves, “Input-to-state stability of rate-controlled biochemical networks,” SIAM J. Control Optim., vol. 44, pp. 704–727, 2005. [4] P. De Leenheer, B. Li, and H.L. Smith, “Competition in the chemostat: some remarks,” Canadian Applied Mathematics Quarterly, vol. 11, pp. 229–248, 2003. [5] P. De Leenheer and S.S. Pilyugin, “Feedback-mediated oscillatory coexistence in the chemostat.” In Positive Systems, C. Commault and N. Marchand, Eds. Berlin: Springer-Verlag, 2006, pp. 97–104. [6] P. De Leenheer and H.L. Smith, “Feedback control for chemostat models,” J. Math. Biol., vol. 46, pp. 48-70, 2003. [7] J.-L. Gouz´e and G. Robledo, “Feedback control for competition models with different removal rates in the chemostat,” INRIA rapport de Recherche 5555, 2005. [8] J.-L. Gouz´e and G. Robledo, “Feedback control for nonmonotone competition models in the chemostat,” Nonlinear Analysis: Real World Applications, vol. 6, pp. 671–690, 2005. [9] J.-L. Gouz´e and G. Robledo, “Robust control for an uncertain chemostat model,” International J. Robust Nonlinear Control, vol. 16, pp. 133-155, 2006. [10] F. Grognard, F. Mazenc, and A. Rapaport, “Polytopic Lyapunov functions for the stability analysis of persistence of competing species.” In Proc. 44th IEEE Conf. Decision Control and European Control Conference ECC 2005, Seville, Spain, 2005, pp. 3699-3704. [11] G. Hardin, “The competitive exclusion principle,” Science, vol. 131, pp. 1292-1298, 1960. [12] S.B. Hsu, “Limiting behavior for competing species,” SIAM J. Applied Mathematics, vol. 34, pp. 760-763, 1978. [13] B. Li, “Asymptotic behavior of the chemostat: General responses and different removal rates,” SIAM J. Applied Mathematics, vol. 59, pp. 411-422, 1998. [14] F. Mazenc, M. Malisoff, and P. De Leenheer, “On the stability of periodic solutions in the perturbed chemostat,” Mathematical Biosciences and Engineering, vol. 4, no. 2, pp. 319-338, 2007. [15] F. Mazenc, M. Malisoff, and J. Harmand, “Stabilization of a periodic trajectory for a chemostat with two species.” In Proc. 2007 American Control Conf., New York, NY, 2007, pp. 6128-6132. [16] F. Mazenc, M. Malisoff, and J. Harmand, “Stabilization in a twospecies chemostat with Monod growth functions,” preprint. [17] N. Rao and E. Roxin, “Controlled growth of competing species,” SIAM J. Applied Mathematcs, vol. 50, pp. 853–864, 1990. [18] H.L. Smith and P. Waltman, The Theory of the Chemostat. Cambridge: Cambridge University Press, 1995. [19] E.D. Sontag, “Smooth stabilization implies coprime factorization,” IEEE Trans. Automat. Control, vol. 34, pp. 435-443, 1989. [20] E.D. Sontag, “Comments on integral variants of ISS,” Systems Control Lett., vol. 34, pp. 93-100, 1998. [21] E.D. Sontag, “Input-to-state stability: Basic concepts and results.” In Nonlinear and Optimal Control Theory, P. Nistri and G. Stefani, Eds. Berlin: Springer-Verlag, 2006, pp. 163-220. [22] E.D. Sontag and Y. Wang, “On characterizations of the input-to-state stability property,” Systems Control Lett., vol. 24, pp. 351-359, 1995. [23] G. Wolkowicz and Z. Lu, “Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates,” SIAM J. Applied Mathematics, vol. 52, pp. 222-233, 1992.
3938
ThC01.6
Stabilization and Robustness Analysis for a Chemostat Model with Two Species and Monod Growth Rates via a Lyapunov Approach Fr´ed´eric Mazenc
Michael Malisoff
Abstract— We study a two species chemostat model with one limiting substrate. We design feedback controllers so that an equilibrium with arbitrary prescribed species concentrations becomes globally asymptotically stable. The novelty of our work is that we assume that only a linear combination of the species concentrations is available for measurement, combined with our use of Lyapunov function methods to generate stable coexistence and quantify the effects of disturbances. Our closed-loop error dynamics are integral input-to-state stable with respect to small perturbations of the controllers. Hence, no matter what initial positive levels for the species concentrations and substrate are selected, the long term species and substrate levels remain close to the equilibrium, even when there are small unexpected changes in the dilution rate and input nutrient concentration. This is a highly desirable robustness property, because the dilution rate is prone to actuator errors. We illustrate our approach using a numerical example. Index Terms—Bioreactors, chemostats, global stabilization, uncertain states, integral input-to-state stability
I. I NTRODUCTION This note continues our search (begun in [14], [15]) for ways to stabilize prescribed equilibrium behaviors in chemostats with one limiting substrate. The chemostat is known in engineering as the continuously stirred tank reactor. See for example [6], [18] for the fundamental role of chemostats in bioengineering, ecology, and population biology; and see Sections II-IV below for the relevant definitions and discussion of the known results. The novelty of the present work is in our design of controllers that stabilize prescribed equilibria (with arbitrary species concentrations) when only an arbitrary linear combination of the species concentrations is available for measurement, combined with our use of a strict Lyapunov function to quantify the effects of actuator errors. Hence, our work leads to species coexistence and persistence (via new strict Lyapunov function and feedback designs) in important cases where there is uncertainty regarding both the current state and the feedback controller. The basic model for a (well mixed) chemostat with two competing species and one limiting substrate is ( 1 (s) X1 − µ2Y(s) X2 , s˙ = D(sin − s) − µY 1 2 (1) X˙ i = [µi (s) − D]Xi , i = 1, 2, Supported by NSF Grant 0424011. Mazenc is with Projet MERE INRIA-INRA, UMR Analyse des Syst`emes et Biom´etrie INRA, 2, pl. Viala, 34060 Montpellier, France (email: [email protected]). Malisoff is with the Department of Mathematics, Louisiana State University, Baton Rouge, LA 70803-4918 USA (email: [email protected]). J. Harmand is with Institut National de la Recherche Agronomique (INRA), UR050, Laboratoire de Biotechnologie de l’Environnement, Narbonne, F-11100, France (email: [email protected]).
1-4244-1498-9/07/$25.00 ©2007 IEEE.
J´erˆome Harmand
evolving on X := (0, ∞)3 ; see [18, Chapter 1] for the derivation of the model. Here and in the sequel, dots indicate time derivatives; s(t) is the concentration of the substrate; X1 (t) and X2 (t) are the concentrations of the two species of organism; the dilution rate D(·) and the input nutrient concentration sin (·) are nonnegative controllers we will specify; µ1 and µ2 are given uptake functions; and Y1 and Y2 are positive constant yield coefficients (but see Section III below for more general models with uncertainties). We assume throughout that the µi s have the Monod form µi (s) =
Ki s Li + s
(2)
where Ki and Li are positive constants we indicate below. Since D is the ratio of the volumetric flow rate (with units of volume over time) to the constant reactor volume, it is proportional to the speed of the pump that supplies the reactor with fresh medium containing the nutrient. The well known “competitive exclusion principle” implies that in classical chemostats with one limiting substrate and constant D and sin , at most one species can survive [11]. This is at odds with the observation that in real ecological systems, it is common for many species to coexist in equilibrium, even if there is only one limiting nutrient. This paradox has motivated a great deal of ongoing research. Much of the research is devoted to choosing D and sin to be time-varying or state-dependent controllers that force a stable coexistence in which multiple species survive [6], [8]. Another important issue is that the species concentrations may not be available for measurement, and there may also be actuator errors caused e.g. by unexpected variability in the speed of the pump supplying the reactor with fresh nutrient. This note addresses all of these issues by designing sin and a controller D(·), depending only on a linear combination Y = X1 + AX2
(3)
of the species concentrations (where A is a given positive constant), that globally stabilize a prescribed equilibrium (s? , x1? , x2? ) ∈ X (with arbitrary positive prescribed species concentrations x1? and x2? , and with s? determined by the µi s using Assumption 1 below). Therefore, no matter what initial values are selected, the chemostat trajectories in closed loop with our feedbacks all converge to (s? , x1? , x2? ). Another important feature of our work is that our use of an explicit strict Lyapunov function construction makes it possible to quantify the effects of disturbances using the input-to-state stability (ISS) paradigm that was introduced by Sontag in his seminal paper [19]. See Section III for the
3933
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 relevant definitions; and see [3], [21] for the fundamental role of ISS in nonlinear control theory and engineering applications. In fact, to our knowledge, our work provides the first ISS analysis for multiple species chemostats whose input nutrient concentrations, outputs, dilution rates, and uptake functions are all perturbed by small noise. In Section II, we state our assumptions and stability result for (1). Section III explains how the stability provided by our feedback is robust to small actuator errors. We compare our work with the known results in Section IV. In Section V, we prove our result for (1). We prove our robustness result in Section VI and we illustrate our methods in Section VII. We close in Section VIII with suggestions for further research. II. A SSUMPTIONS
AND
S TABILITY T HEOREM
Let SAT denote the set of all continuously differentiable functions σ : R → R satisfying (I) r 7→ rσ(r) is positive definite (i.e., zero at zero and strictly positive at all other points) and (II) |σ(r)| ≤ min{|r|, 1} for all r ∈ R (which we refer to as saturations). For example, σ ∈ SAT if r . (4) σ(r) = √ 1 + r2 Given µi s as in (2) (where the Ki s and Li s are given positive constants), set χ(s) = µ2 (s) − µ1 (s). We always assume: Assumption 1: There is a constant s? > 0 such that (i) χ(s? ) = 0, (ii) χ(s) < 0 when 0 < s < s? , and (iii) χ(s) > 0 when s > s? . In particular, µ1 (s? ) = µ2 (s? ). Simple calculations and our Assumption 1 readily yield Li µi (s) (s − r) ∀s, r > 0 Li + r s for i = 1, 2 and (using the formula for χ) µi (s) − µi (r) =
s? =
L1 K2 −L2 K1 K1 −K2
, µ1 (s? ) =
L1 K2 −L2 K1 L1 −L2
,
(5)
again evolving on X = (0, ∞)3 (where sin and D are controllers to be chosen) and the output y = Y /Y1 satisfying
(8)
We always assume that a 6= 1. Given arbitrary prescribed x1? > 0 and x2? > 0, we set sin? = s? + x1? + x2? ; x˜i = xi − xi? for i = 1, 2; s˜ = s − s? ; ξ˜i = ln(˜ xi + xi? ) − ln(xi? ) = ln(xi /xi? ) for ˜ = ln(˜ i = 1, 2; and Σ s + s? ) − ln(s? ) = ln(s/s? ).
β : [0, ∞) × [0, ∞) → [0, ∞) for which (1) β(·, t) ∈ K∞ for each t ≥ 0 and (2) β(r, ·) is non-increasing and β(r, t) → 0 as t → +∞ for each r ≥ 0. Let δB2 denote the closed ball of radius δ > 0 centered at 0 ∈ R2 . We use the constants ω1 = ω3 = and
, ω2 =
µ1 (s? ) 4s? La − L1 +1 1 2
,
µ1 (s? )s? , 20 (a−1)sin? +s2? La − L1 +1 1 2 0.72µ1 (s? )s2? L1 − L1 1 2 h i ω4 = 2 (a−1)sin? +s2? La − L1 +1
(10) .
2
The following yields (s(t), x1 (t), x2 (t)) → (s? , x1? , x2? ) as t → +∞ for all trajectories (s(t), x1 (t), x2 (t)) of (7) (i.e., for all initial values in X ) when the indicated controllers are used, hence globally asymptotically stable coexistence: Theorem 1: Let ε ∈ (0, min{ω1 , ω2 , ω3 , ω4 }], σ ∈ SAT , and x1? , x2? > 0 be given. Let s? > 0 satisfy the requirements of Assumption 1. Consider (7) with output (8) with a 6= 1. Then (s? , x1? , x2? ) is a globally asymptotically stable equilibrium of (7) when sin = sin? := s? + x1? + x2? and D(y) = µ1 (s? ) − ε(a − 1)σ (y − x1? − ax2? ). More precisely, there exists β ∈ KL such that for all trajectories (s, x1 , x2 )(t) of (7) with sin and D (for all initial values ˜ ξ˜1 , ξ˜2 )(t) from (9) in X ), the transformed error vector (Σ, ˜ ˜ ˜ ˜ ˜ ˜ satisfies |(Σ, ξ1 , ξ2 )(t)| ≤ β(|(Σ, ξ1 , ξ2 )(0)|, t) for all t ≥ 0. Remark 1: The controller D from Theorem 1 consists of µ1 (s? ) plus a term that can be made arbitrarily small by picking a sufficiently small ε > 0, and it yields coexistence because x1? > 0 and x2? > 0. However, we cannot pick ε = 0 since that would make D and sin constants, and so yield competitive exclusion [6]. III. ROBUSTNESS
TO
D ISTURBANCES
Consider the perturbed chemostat error dynamics s˜˙ = [D(y)+d2 ](sin +d1 −s) − µ1 (s)x1 − µ2 (s)x2 x ˜˙ i = [µi (s) − D(y) − d2 ]xi , i = 1, 2
(11)
with output (8) for a 6= 1, where s˜ = s−s?; x ˜i = xi −xi? for i = 1, 2; sin and D are from Theorem 1; σ is from (4); the µi s are as before (so Assumption 1 holds); and the unknown measurable disturbance d(t) = (d1 (t), d2 (t)) is bounded in sup norm |·|∞ by a constant that we specify below. For each ρ > 0, let D2 (ρ) denote the set of all measurable essentially bounded d : [0, ∞) → ρB2 . Fix arbitrary x1? > 0 and x2? > 0, and let ε > 0 be as in the statement of Theorem 1. To simplify the notation, we set 2? L2 +s? , β¯ = xx1? , L2 ˜ ˜ , A0 (˜ s) = s˜ − s? ln 1 + ss˜? = s? eΣ − 1 − Σ ˜ and Ai (ξ˜i ) = eξi − 1 − ξ˜i = xx˜i?i − ln 1 + xx˜i?i
α ¯=
(9)
For each trajectory (s, x1 , x2 )(t) of (7), we call the vector ˜ ξ˜1 , ξ˜2 )(t) as defined by (9) the (corresponding trans(Σ, formed) error vector. Let K∞ denote the set of all continuous ρ : [0, ∞) → [0, ∞) for which (i) ρ(0) = 0 and (ii) ρ is strictly increasing and unbounded. If ρ1 , ρ2 ∈ K∞ , then ρ−1 1 , ρ1 ◦ ρ2 ∈ K∞ . Let KL denote the set of all continuous
µ1 (s? ) 20|a−1|
1
(6)
and therefore K2 > K1 , K1 L2 > K2 L1 , and L2 > L1 . We transform (1) with the output Y from (3), using x1 = X1 /Y1 and x2 = X2 /Y2 . Doing so produces the system s˙ = D(sin − s) − µ1 (s)x1 − µ2 (s)x2 , (7) x˙ i = [µi (s) − D]xi , i = 1, 2
y = x1 + ax2 , where a = AY2 /Y1 .
ThC01.6
L1 +s? L1
(12) (13)
for i = 1 and 2, where we used (9). We show that (11) is integral input-to-state stable (iISS) ˜ ξ˜1 , ξ˜2 ) (from (9) above) for distur[20] in the variable (Σ, ¯ ¯ is defined in (15), in the sense bances d ∈ D2 (∆) when ∆
3934
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 that there exist α ∈ K∞ and β ∈ KL and a constant K > 0 ¯ and all corresponding closed-loop so that for all d ∈ D2 (∆) trajectories (˜ s(t), x ˜1 (t), x ˜2 (t)) of (11) for our controllers D and sin (for all initial values), the corresponding transformed ˜ error vector (Σ(t), ξ˜1 (t), ξ˜2 (t)) satisfies ˜ ξ˜1 , ξ˜2 )(t)|) α(|(Σ, R ˜ ξ˜1 , ξ˜2 )(0)|, t) + K t |d(r)|dr ∀t ≥ 0 . ≤ β(|(Σ, 0
(iISS)
(This is a stronger condition than the usual iISS condition from [20] since it involves the integral of |d(r)| instead of the usual integral of a K∞ function of |d(r)|.) We also show that (11) is input-to-state stable (ISS) in ˜ ξ˜1 , ξ˜2 ) when d ∈ D2 (∆ ¯ [ ) and ∆ ¯ [ is in (41), meaning (Σ, ¯ [) there are β ∈ KL and γ ∈ K∞ so that for all d ∈ D2 (∆ and all closed-loop trajectories (˜ s, x ˜1 , x ˜2 )(t) of (11) for our ˜ ξ˜1 , ξ˜2 )(t) satisfies D and sin (for all initial values), (Σ, ˜ ξ˜1 , ξ˜2 )(t)| |(Σ, ˜ ξ˜1 , ξ˜2 )(0)|, t) + γ(|d|∞ ) ∀t ≥ 0. ≤ β(|(Σ,
(ISS)
See [3], [21] for the fundamental roles of iISS and ISS in nonlinear control and engineering. The ISS and iISS estimates give (s, x1 , x2 )(t) → (s? , x1? , x2? ) as t → +∞ for all trajectories (s, x1 , x2 )(t) of (11) when the indicated controllers D and sin are used and d = 0, and they also quantify the effects of actuator errors d. Also, (ISS) guarantees persistence, since if for example we had x2 (t) → ˜ 0 as t → t− ? for some t? ∈ (0, +∞], then ξ2 (t) = ln(x2 (t))− ln(x2? ) → −∞, which would violate (ISS). Similarly, (iISS) implies persistence when d is integrable. To specify the bounds on d, we introduce the constants C1 = s10? L11 − L12 ε 2 L1 +s? L2 +s? ε2 a − + 1 C2 = 19µ20 4C1 L1 L2 1 (s? ) (14) n o ?) κ1 = min 1, µ1 (s +x +x +µ (s )+ε|a − 1| 1? 2? 1 ? 5 n o κ1 1 κ2 = 2s2 + 2C2 + sin? κ1 C2 +max α, ¯ xx1? β¯ . 2? in?
See Section VII below for examples where the bounds on d are computed explicitly. Our iISS result is as follows (but see Remark 3 for ISS results under appropriate conditions): Theorem 2: Under the preceding assumptions with ¯ = min 1, 2µ1 (s? )s? , 1 ∆ (15) 5(3κ1 + 2s? ) 2κ2
˜ ξ˜1 , ξ˜2 ) for disturbances d ∈ D2 (∆). ¯ (11) is iISS in (Σ, Remark 2: A much more general scenario is where there are measurement errors in D(·) and uncertain uptake functions, as well as actuator errors. This gives the error dynamics s˜˙ x ˜˙ 1 x ˜˙ 2
= [D(y + r3 ) + u2 ](sin + u1 − s)
−(1 + r1 )µ1 (s)x1 − (1 + r2 )µ2 (s)x2 ,
= [(1 + r1 )µ1 (s) − D(y + r3 ) − u2 ]x1 , = [(1 + r2 )µ2 (s) − D(y + r3 ) − u2 ]x2 .
(16)
ThC01.6 d = (r, u) for suitable small values of the disturbance. Due to space constraints, we leave this extension to the reader. IV. C OMPARISON
WITH
K NOWN R ESULTS
The importance of coexistence and stability in chemostats has motivated a great deal of significant research at the interface of biology and mathematical control theory [2], [4], [5], [6], [23]. Much of this work uses monotonicity or reduction methods to establish stability. By contrast, Lyapunov functions have not frequently been used to study stability in chemostats with multiple species, and where they are used they are often weak (i.e., nonstrict) Lyapunov functions which do not lend themselves to robustness analysis. An exception is the famous one from [12]; see also Theorem 4.1 on [18, p.35] and [10]. Strict Lyapunov functions for two species chemostats were explicitly constructed in [15], but no robustness to disturbances was established there. However, [15] provides simple linear feedback stabilizers that yield local asymptotic stability of a prescribed periodic trajectory for the closed-loop dynamics (i.e. tracking). See [13], [23] where weak Lyapunov functions are used in conjunction with variants of the LaSalle Invariance Principle. This raises the important question of whether strict Lyapunov functions can be explicitly constructed for multispecies chemostats that are globally stabilized through static output feedbacks, and if so, whether the Lyapunov functions can be used to quantify the effects of uncertainty. For chemostats with only one species that are made oscillating through a suitable choice of D(·), this type of problem was solved in [14]. There, a Lyapunov function is constructed and yields ISS relative to small actuator uncertainty. Here we extend these results to two-species models with disturbance perturbations, or where only a linear combination of the species concentrations is known. Our work owes a great deal to [6], [7], which stabilize chemostats in which only the sum of the species concentrations can be measured, using an appropriate dilution rate. However, [6] and [7] do not include our work since [6] does not rely on a Lyapunov approach and the Lyapunov functions from [7] are nonstrict and so do not lend themselves to ISS. See [9] for results where only the substrate level (multiplied by an a priori bounded error) is available for measurement. V. P ROOF
OF
T HEOREM 1
In all of what follows, all (in)equalities should be understood to hold globally unless otherwise indicated, we omit the arguments of our functions whenever this would not lead to confusion, sin = s? + x1? + x2? , and D? = µ1 (s? ). Since s? is constant and µ1 (s? ) = µ2 (s? ), the change of feedback D = D? + v, (5), and (9) transform (7) into ( s˜˙ = (D? + v)(sin − s) − µ1 (s)x1 − µ2 (s)x2 , (17) ˙˜ µi (s) i ξi = LiL+s ˜ − v, i = 1, 2 . s s ? Step 1: Construction of a weak Lyapunov function
A variant of our proof of Theorem 2 shows that (16) is ISS in ˜ ξ˜1 , ξ˜2 ) with respect to the disturbance vector the variable (Σ,
3935
We first show that in terms of (12) and (13), 1 ¯ 2 (ξ˜2 ) V (˜ s, ξ˜1 , ξ˜2 ) = A0 (˜ s) + α ¯ A1 (ξ˜1 ) + βA x1?
(18)
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 is a weak Lyapunov function for (17) in the sense that 1 s˜ x˜1 x ˜2 V˙ = −H(˜ s)˜ s2 + (sin −s)v − α ¯ + β¯ v (19) x1? s x1? x2?
holds along all trajectories of (17), where H as defined by h 1 (s? ) H(˜ s) = sx11? µ1 (s? ) + µ1 (s)−µ x1? s˜ i (20) µ2 (s)−µ2 (s? ) x , s ˜ = 6 0 + 2? s˜
is everywhere positive. To this end, first notice that by (9), s˜ x ˜1 ˜˙ x˜2 ˜˙ A˙ 0 = s˜˙ , A˙ 1 = ξ1 , and A˙ 2 = ξ2 . s x1? x2?
(21)
Hence, along the trajectories of (17), we easily get h i 1 2 V˙ = x11? ss˜ (D? + v)(sin − s) − α ¯ xx˜1? + β¯ xx˜2? v h i µ (s) L1 ˜1 x11? 1s s˜ + −x1 + α ¯ L1 +s? x h i L2 1? + −x2 + β¯ xx2? ˜2 x11? µ2s(s) s˜. L2 +s? x =
1 x1?
[D? (sin − s) − xh1? µ1 (s) − x2?iµ2 (s)] ss˜ 1 2 ¯ xx˜1? + β¯ xx˜2? v. + x11? ss˜ (sin − s)v − α
Step 2: Construction of a strict Lyapunov function Set y˜ = x ˜1 + a˜ x2 , Z˜ = s˜ + x ˜1 + x ˜2 , Ψ(r) = 21 r2 , and
2s?
L1
−L
(23)
where V is from (18). We show that V2 is a strict Lyapunov function for (17) in the sense that its time derivative along the trajectories of (17) with v = −ε(a − 1)σ(˜ y ) satisfies 4s
10
1 L1
i − L12 ε˜ yσ (˜ y)+ Z˜ 2
(24)
and x˜i and ξ˜i are related by (9). Then W2 is positive definite (since a 6= 1, L2 > L1 , and r 7→ rσ(r) is positive definite). To this end, first notice that by our choices of ω1 , ε, and D, we get D(y) ≥ 19 20 D? everywhere. Therefore, (7) gives Z˜˙ = D(y)[sin − s − x1 − x2 ] = −D(y)Z˜ d ˜ = −D(y)Z˜ 2 ≤ − 19 D? Z˜ 2 . hence dt Ψ(Z) 20
(25)
Moreover, our assumption that a 6= 1 gives
a(Z˜ − s˜) − y˜ y˜ − Z˜ + s˜ x ˜1 = and x ˜2 = . a−1 a−1
−L
2
Applying the general relation ab ≤ 41 a2 + b2 , we get n o (a−1)sin +s2 a − 1 s˜ p ? L1 L2 √ µ1 (s? ) √ εσ(˜ y ) s µ1 (s? )s µ1 (s? ) 2 ˜ 4s s
+
h
(a−1)sin +s2?
a L1
µ1 (s? )s
− L1
2
i2
ε2
+
h
(a−1)sin +s2?
a L1
− L1
2
µ1 (s? )s
i2
ε2
(28)
2
σ (˜ y)
for all s > 0. Combining (27) and (28), it follows that h i (s? ) 2 1 1 ? V˙ 2 ≤ − µ12s s˜ − 9s y σ (˜ y ) − Z˜ 2 10 L1 − L2 ε˜
(29)
σ 2 (˜ y) .
We next distinguish between two cases. Case A: If s ≥ then (29) and our choices of ε and ω4 give h i (s? ) 2 1 1 ? V˙ 2 ≤ − µ12s s˜ − 9s − y σ (˜ y ) − Z˜ 2 10 L1 L2 ε˜
9 10 s? ,
+
h i2 10 (a−1)sin +s2? La − L1 ε2 1
2
9µ1 (s? )s?
(s? ) 2 ≤ − µ12s s˜ −
s? 10
h
1 L1
−
1 L2
i
σ 2 (˜ y)
(30)
ε˜ y σ (˜ y ) − Z˜ 2
9 because σ 2 (˜ y ) ≤ y˜σ(˜ y ) everywhere. Case B: If s ≤ 10 s? , s? then |˜ s| ≥ 10 , so (29) and our choices of ε and ω3 give h i (s? ) 2 1 (s? ) 2 V˙ 2 ≤ − µ14s s˜ − µ400s s? − s2? L11 − L12 ε˜ y σ (˜ y)
−Z˜ 2 +
h
(a−1)sin +s2?
(s? ) 2 ≤ − µ14s s˜ −
(26)
Substituting (12) and (26) into (19), we easily get h i s? s? 1 V˙ = −H(˜ s)˜ s2 + x11? 1s sin + a−1 a − s˜v + L L2 h i h 1 i s? L1 +s? L2 +s? ˜ 1 1 1 − y ˜ v − a − Zv. x1? (a−1) L1 L2 x1? (a−1) L1 L2
L1
It follows from (25) and our formula (23) for V2 that h i 1 1 ? − y σ (˜ y ) − Z˜ 2 V˙ 2 ≤ − 3µ14s(s? ) s˜2 − 9s 10 L L2 ε˜ h 1 i (27) s˜ − (a − 1)sin + s2? La1 − L12 εσ (˜ y ) . s
≤
2
V˙ 2 ≤ −W2 , where h W2 (˜ s, ξ˜1 , ξ˜2 ) = µ1 (s? ) s˜2 + s?
By our choices of ε and ω2 , the fact that σ 2 (˜ y) ≤ y˜σ(˜ y) δ 2 5 2 everywhere, and the general relation ab ≤ a + b for 2δ 10 √ b = εσ(˜ y ) and δ = s? ( L11 − L12 ) > 0, we then easily get h i y σ (˜ y) x1? V˙ ≤ − 3µ14s(s? ) s˜2 − s? L11 − L12 ε˜ h i s˜ − (a − 1)sin + s2? La1 − L12 y) s εσ (˜ n√ h i o √ ? ? ε L1L+s a − L2L+s Z˜ { εσ (˜ y)} + 1 h 2 i 1 1 ? ≤ − 3µ14s(s? ) s˜2 − 9s y σ (˜ y) 10 L − L2 ε˜ h 1 i s˜ − (a − 1)sin + s2? La1 − L12 y) s εσ (˜ h i2 ? ? + 15ε 1 L1L+s a − L2L+s Z˜ . 1 2
(22)
The formulas sin = s? + x1? + x2? and D? = µ1 (s? ) and the fact that µ1 (s? ) = µ2 (s? ) therefore easily give (19).
V2 (˜ s, ξ˜1 , ξ˜2 ) = x1? V (˜ s, ξ˜1 , ξ˜2 ) + 2 20 ˜ 5ε L1 +s? a − L2 +s? + 1 Ψ(Z) 19D? 1 L1 L2 1
Taking v = −ε(a − 1)σ(˜ y ) and recalling (20) gives h i (a−1)sin s˜+s? s? s? ?) 2 s ˜ − + a− s˜εσ (˜ y) x1? V˙ ≤ − µ1 (s s s s L1 L h i h i 2 ? ? ˜ (˜ −s? L11 − L12 ε˜ yσ (˜ y ) + L1L+s a − L2L+s Zεσ y) . 1 2
2s?
Since −xi + x ˜i = −xi? for i = 1, 2, our choices (12) give V˙
ThC01.6
a L1
− L1
2
i2
ε2
h µ1 (s? )s i s? 1 1 y σ (˜ y) 2 L1 − L2 ε˜
(31)
− Z˜ 2
where we used |σ(˜ y )| ≤ 1. In either case, we get (24). Since W2 is positive definite, V2 is the desired strict Lyapunov function. However, V2 becomes unbounded as s˜ → −s? so extra care is needed to obtain β ∈ KL from the statement of the theorem; see [16] for details. This proves Theorem 1.
3936
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007
ThC01.6
VI. P ROOF OF T HEOREM 2 ¯ and let ∆ denote its sup norm, so We fix d ∈ D2 (∆) ¯ We use the notation from the proof of Theorem 1. ∆ ≤ ∆. We take sin = sin? := s? + x1? + x2? throughout and ¯ 1? /x2? sin? C3 = 2κ1 + 2 2C2 κ1 + 2 max α, ¯ βx +8C2 s2in? and, in terms of the relations (9), (32) (s? ) 2 W4 (˜ s, ξ˜1 , ξ˜2 ) = µ120s s˜ + C1 y˜σ (˜ y ) + 21 Z˜ 2 ,
¯ 1? /x2? |˜ ¯ 1? /x2? }Z˜ 2 /sin? . ∆{¯ α|˜ x1 |+ βx x2 |} ≤ ∆ max{¯ α, βx Substituting into (38) and recalling (14),
where σ is from (4). Rearranging terms in (11) gives s˜˙ = D(y)[sin − s] − µ1 (s)x1 − µ2 (s)x2 +g1 (˜ s, y, d), ˜˙ ξi = µi (s) − D(y) − d2 , i = 1, 2
¯ ≤ 1/(2κ2 ), we get V˙ 2 ≤ −W4 (˜ Since ∆ s, ξ˜1 , ξ˜2 ); see (32). 2nd case: |˜ s + x˜1 + x˜2 | ≤ 2sin? . One easily verifies [16]: Lemma 3: If |˜ s+x ˜1 + x ˜2 | ≤ 2sin? , then |˜ s|+|˜ x1 |+|˜ x2 | ≤ 4sin? . Using Lemma 3 and (38), we deduce in the second case that V˙ 2 ≤ −W3 (˜ s, ξ˜1 , ξ˜2 ) + C3 ∆ ≤ −W4 (˜ s, ξ˜1 , ξ˜2 ) + C3 ∆, so
(33)
with g1 (˜ s, y, d) = d2 (−˜ s + x1? + x2? ) + [D(y) + d2 ]d1 . The time derivative of V2 along the trajectories of (33) satisfies: i h ∂V2 ∂V2 2 V˙ 2 ≤ −W2 + ∂V g (˜ s , y, d) − d + , (34) 1 2 ˜ ˜ ∂ s˜ ∂ξ ∂ξ 1
2
by (24). In terms of the constants from (14) and the relations (9), we have |g1 (˜ s, y, d)| ≤ ∆(κ1 + |˜ s|), (s? ) 2 W2 (˜ s, ξ˜1 , ξ˜2 ) = µ14s s˜ + C1 y˜σ (˜ y) + Z˜ 2 , ∂V2 ˜ ∂V2 = α˜ ˜ 1, = s˜ + C2 Z, ¯ x1 + C2 Zx ∂ s˜
and
s ∂V2 ˜ ∂ ξ2
=
x1? ¯ ˜2 x2? β x
∂ ξ˜1
(35)
˜ 2, + C2 Zx
¯ ≤ min{1, µ1 (s? )/5} (which holds rewhere we used ∆ gardless of the choice of κ1 > 0, by (15)) to get κ1 . Set v˜ = |˜ s| + |˜ x1 | + |˜ x2 |. From (34), we get V˙ 2
2 1? ¯ ≤ −W2 + ∆ s˜s + κ1 ∆ |˜ss| + ∆ xx2? β|˜ x2 | ˜ (|x1 | + |x2 | + κ1 + |˜ + ∆C2 |Z| s|) + ∆¯ α|˜ x1 | |˜ s| x1? ¯ s˜2 ≤ −W2 + ∆ s + κ1 ∆ s + ∆ x2? β|˜ x2 | ˜ (2κ1 + v˜) + ∆¯ + ∆C2 |Z| α|˜ x1 |,
(36)
where we used |x1 | + |x2 | + |˜ s| ≤ v˜ + x1? + x2? ≤ v˜ + κ1 . One easily shows that [16]: Lemma 1: For all s > 0, we have 3˜ s2 |˜ s| ≤2+ . (37) s 2s? s Combined with (36), Lemma 1 gives 2 s˜ 1 V˙ 2 ≤ −W2 (˜ s, ξ˜1 , ξ˜2 ) + ∆ 1 + 3κ 2s? s + 2κ1 ∆ ¯ x2 | . ˜ +∆C2 |Z| (2κ1 + v˜) + ∆¯ α|˜ x1 | + ∆ x1? β|˜
V˙ 2
∆κ1 ≤ −W3 (˜ s, ξ˜1 , ξ˜2 ) + 2s (˜ s+x ˜1 + x˜2 )2 2 h i in? 1 +∆C2 sκin? + 2 (˜ s + x˜1 + x ˜ 2 )2 ¯ 1? /x2? }(˜ + s∆ max{¯ α, βx s + x˜1 + x ˜ 2 )2 in? = −W3 (˜ s, ξ˜1 , ξ˜2 ) + ∆κ2 (˜ s + x˜1 + x ˜ 2 )2 .
V˙ 2 ≤ −W4 (˜ s, ξ˜1 , ξ˜2 ) + C3 ∆
holds in both cases. Define V and the positive definite W by ˜ ξ˜1 , ξ˜2 ) = V2 ([eΣ˜ − 1]s? , ξ˜1 , ξ˜2 ) and V(Σ, ˜ ξ˜1 , ξ˜2 ) = W4 ([eΣ˜ − 1]s? , ξ˜1 , ξ˜2 ). W(Σ,
1? ¯ V˙ 2 ≤ −W3 (˜ s, ξ˜1 , ξ˜2 )+2κ1∆+∆¯ α|˜ x1 |+∆ xx2? β|˜ x2 |
+∆C2 |˜ s+x ˜1 + x˜2 | [2κ1 + v˜]
(38)
(s? ) 2 where W3 (˜ s, ξ˜1 , ξ˜2 ) = µ120s s˜ + C1 y˜σ (˜ y) + Z˜ 2 . We consider two cases, based on the magnitude of Z˜ = s˜ + x ˜1 + x ˜2 : ˜ 1st case: |Z| ≥ 2sin? . One easily checks that [16]: ˜ ≥ 2sin? , then |˜ ˜ Lemma 2: If |Z| s| + |˜ x1 | + |˜ x2 | ≤ 2|Z|. 2 2 We deduce that ∆κ1 ≤ ∆κ1 (˜ s + x˜1 + x˜2 ) /{4sin? }, h i 1 ˜ and ∆C2 [2κ1 + v˜] ≤ ∆C2 sκin? + 2 |Z|,
(40)
By (39), V˙ ≤ −W +C3 ∆. We can also find α1 , α2 ∈ K∞ so ˜ ξ˜1 , ξ˜2 )|) ≤ V(Σ, ˜ ξ˜1 , ξ˜2 ) ≤ α2 (|(Σ, ˜ ξ˜1 , ξ˜2 )|) [16]. that α1 (|(Σ, Hence, a standard argument [1, p.1088] gives our estimate (iISS) with α = α1 and K = 2C3 , proving Theorem 2. Remark 3: One can prove that the system (11) is also ISS ˜ ξ˜1 , ξ˜2 ) when d ∈ D2 (∆[ ), where, in terms of (15), in (Σ, 1 [ ¯ ¯ ∆ = min ∆, mini α ˜ 1 (xi? /2) ; (41) 2C3 α ˜ 1 ∈ K∞ is chosen so that α ˜ 1 (r) ≤ r2 and W4 (˜ s, ξ˜1 , ξ˜2 ) ≥ ˜ α ˜ 1 (|(˜ s, x ˜1 , x˜2 )|) hold everywhere; x ˜i and ξi are again related by (9); and C3 is from (32). For the proof, see [16]. Remark 4: An important special case of Theorem 2 is where only the dilution rate is affected by noise; i.e., d1 = 0 and d2 6= 0. This often occurs when the speed of the pump that supplies the chemostat with fresh nutrient is subjected to small fluctuations [14]. When d1 = 0, our proof of Theorem 2 remains valid if κ1 is replaced by x1? + x2? . Another important case is where the noise term is only on the input nutrient concentration sin , i.e., d2 = 0 and d1 6= 0, in which ¯ is replaced by case Theorem 2 remains true if ∆
x2?
From the formula for W2 in (35), Z˜ = s˜ + x ˜1 + x˜2 , our ¯ and our choice (15) of ∆, ¯ we get assumption that ∆ ≤ ∆,
(39)
¯ = ∆
0.16D?s? ; D? + ε|a − 1|
(42)
see [16] for the proof. We illustrate these special cases below. Remark 5: Our bound (15) is conservative. It can be enlarged by reducing κ2 . For example, arguing as in the proof of Theorem 2 except separately considering the cases |˜ s + x˜1 + x ˜2 | ≥ M sin? and |˜ s + x˜1 + x ˜2 | ≤ M sin? (where M ≥ 2 is any constant), one checks [16] that Theorem 2 still holds if κ2 is replaced by 2κ1 1 2x1? ¯ κ2 = 2 2 + 2C2 + 2κ1 C2 + 2¯ α+ β M sin? M sin? x2? so κ2 ≈ 2C2 for large M , if C3 is appropriately enlarged.
3937
46th IEEE CDC, New Orleans, USA, Dec. 12-14, 2007 VII. E XAMPLE We use the following special case of (7)-(8): 0.05sx1 .052sx2 2 )(sin +d1 −s)− 20+s − 25+s s˙ = (D(y)+d h i .05s x˙ 1 = 20+s − D(y) − d2 x1 h i x˙ 2 = .052s − D(y) − d2 x2 25+s
(43)
with y = x1 + 0.8x2 . We take x1? = 0.05 and x2? = 0.02. Using the notation from Theorem 1, a = 0.8, s? = 105, and ε ∈ (0, .00753]. We take sin = 105.07 and D(y) = .042 + 0.001506σ(y − 0.066) where σ is from (4). Theorem 1 implies that regardless of the initial values, all trajectories of (43) with our controllers D and sin converge to (105, 0.05, 0.02) when d = 0. We next consider the important special case where d1 = 0 but d2 6= 0. In that case, we can take κ1 = x1? + x2? ; see ¯ as in (15) and κ2 as in Remark 5 with Remark 4. Taking ∆ ¯ ≈ 0.20D? . By Theorem 2, the stability M = 25 gives ∆ has the highly desirable property of robustness (in the iISS sense) to the noise d = (0, d2 ), provided |d2 |∞ is no more than (about) 20% of D; in particular, if d2 is also integrable, then we have persistence of both species. If instead the disturbance is only added to sin , then we can use (42) to get robustness to disturbances d1 that are bounded ¯ ≈ 16, or about 15% of sin = 105.07. We simulated by ∆ the corresponding dynamics (43) with d(t) ≡ (1, 0) and the initial value (s, x1 , x2 )(0) = (103, 2, 1) and obtained: Substrate 105.5 105 104.5 104 103.5 1000
2000
Species 1
3000
4000
5000
Species 2
2
1 0.8
1.5
0.6 1 0.4 0.5
0.2
1000
2000
3000
4000
5000
1000
2000
3000
4000
5000
Our simulation illustrates the persistence of (x1 , x2 )(t) and the convergence (s(t), x1 (t), x2 (t)) → (105, 0.05, 0.02), with an overshoot determined by the iISS estimate and the magnitude of d1 , and so validates our theorems. We emphasize that in general, our iISS estimate (iISS) will not guarantee persistence unless we choose d to be integrable. VIII. C ONCLUSIONS We provided an output feedback design for stabilizing equilibria with arbitrary prescribed species concentrations in two species chemostats. The novelty of our work is our construction of a Lyapunov function, which made it possible
ThC01.6 to quantify the effects of uncertainties such as actuator errors, measurement noise, and uncertain growth rates using ISS. It would be of interest to extend our work to chemostats with delays between the instant the micro-organism uptakes the substrate and the “instant” the biomass integrates it, as well as to tracking. See [15] for tracking in two species chemostats where the species concentrations are known. Another desirable extension would be to models with three or more competing species, or more than one limiting substrate. R EFERENCES [1] D. Angeli, E.D. Sontag, and Y. Wang, “A characterization of integral input to state stability,” IEEE Trans. Automat. Control, vol. 45, pp. 1082-1097, 2000. [2] O. Bernard and J.-L. Gouz´e, “Nonlinear qualitative signal processing for biological systems: Application to the algal growth in bioreactors,” Mathematical Biosciences, vol. 157, pp. 357–372, 1999. [3] M. Chaves, “Input-to-state stability of rate-controlled biochemical networks,” SIAM J. Control Optim., vol. 44, pp. 704–727, 2005. [4] P. De Leenheer, B. Li, and H.L. Smith, “Competition in the chemostat: some remarks,” Canadian Applied Mathematics Quarterly, vol. 11, pp. 229–248, 2003. [5] P. De Leenheer and S.S. Pilyugin, “Feedback-mediated oscillatory coexistence in the chemostat.” In Positive Systems, C. Commault and N. Marchand, Eds. Berlin: Springer-Verlag, 2006, pp. 97–104. [6] P. De Leenheer and H.L. Smith, “Feedback control for chemostat models,” J. Math. Biol., vol. 46, pp. 48-70, 2003. [7] J.-L. Gouz´e and G. Robledo, “Feedback control for competition models with different removal rates in the chemostat,” INRIA rapport de Recherche 5555, 2005. [8] J.-L. Gouz´e and G. Robledo, “Feedback control for nonmonotone competition models in the chemostat,” Nonlinear Analysis: Real World Applications, vol. 6, pp. 671–690, 2005. [9] J.-L. Gouz´e and G. Robledo, “Robust control for an uncertain chemostat model,” International J. Robust Nonlinear Control, vol. 16, pp. 133-155, 2006. [10] F. Grognard, F. Mazenc, and A. Rapaport, “Polytopic Lyapunov functions for the stability analysis of persistence of competing species.” In Proc. 44th IEEE Conf. Decision Control and European Control Conference ECC 2005, Seville, Spain, 2005, pp. 3699-3704. [11] G. Hardin, “The competitive exclusion principle,” Science, vol. 131, pp. 1292-1298, 1960. [12] S.B. Hsu, “Limiting behavior for competing species,” SIAM J. Applied Mathematics, vol. 34, pp. 760-763, 1978. [13] B. Li, “Asymptotic behavior of the chemostat: General responses and different removal rates,” SIAM J. Applied Mathematics, vol. 59, pp. 411-422, 1998. [14] F. Mazenc, M. Malisoff, and P. De Leenheer, “On the stability of periodic solutions in the perturbed chemostat,” Mathematical Biosciences and Engineering, vol. 4, no. 2, pp. 319-338, 2007. [15] F. Mazenc, M. Malisoff, and J. Harmand, “Stabilization of a periodic trajectory for a chemostat with two species.” In Proc. 2007 American Control Conf., New York, NY, 2007, pp. 6128-6132. [16] F. Mazenc, M. Malisoff, and J. Harmand, “Stabilization in a twospecies chemostat with Monod growth functions,” preprint. [17] N. Rao and E. Roxin, “Controlled growth of competing species,” SIAM J. Applied Mathematcs, vol. 50, pp. 853–864, 1990. [18] H.L. Smith and P. Waltman, The Theory of the Chemostat. Cambridge: Cambridge University Press, 1995. [19] E.D. Sontag, “Smooth stabilization implies coprime factorization,” IEEE Trans. Automat. Control, vol. 34, pp. 435-443, 1989. [20] E.D. Sontag, “Comments on integral variants of ISS,” Systems Control Lett., vol. 34, pp. 93-100, 1998. [21] E.D. Sontag, “Input-to-state stability: Basic concepts and results.” In Nonlinear and Optimal Control Theory, P. Nistri and G. Stefani, Eds. Berlin: Springer-Verlag, 2006, pp. 163-220. [22] E.D. Sontag and Y. Wang, “On characterizations of the input-to-state stability property,” Systems Control Lett., vol. 24, pp. 351-359, 1995. [23] G. Wolkowicz and Z. Lu, “Global dynamics of a mathematical model of competition in the chemostat: general response functions and differential death rates,” SIAM J. Applied Mathematics, vol. 52, pp. 222-233, 1992.
3938
