{'fraction_non_alphanumeric': 0.05910664527569028, 'fraction_numerical': 0.009710377438149117, 'mean_word_length': 4.13261105092091, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hybrid Control of a Bioreactor with Quantized Measurements: Extended Version

We consider the problem of global stabilization of an unstable bioreactor model (e.g. for anaerobic digestion), when the measurements are discrete and in finite number ("quantized"), with control of the dilution rate. The model is a differential system with two variables, and the output is the biomass growth. The measurements define regions in the state space, and they can be perfect or uncertain (i.e. without or with overlaps). We show that, under appropriate assumptions, a quantized control may lead to global stabilization: trajectories have to follow some transitions between the regions, until the final region where they converge toward the reference equilibrium. On the boundary between regions, the solutions are defined as a Filippov differential inclusion. If the assumptions are not fulfilled, sliding modes may appear, and the transition graphs are not deterministic.


I. INTRODUCTION
Classical control methods are often based on the complete knowledge of some outputs y(t) of the system [22]. By complete, we mean that any output y i is a real number, possibly measured with some noise δ i . The control is then built with this (noisy) measurement. These tools have been successfully applied in many domains of science and engineering, e.g. in the domains of biosystems and bioreactors [10]. However, in these domains, detailed quantitative measurements are often difficult, too expensive, or even impossible. A striking example is the measurements of gene expression by DNA-chips, giving only a Boolean measure equal to on (gene is expressed) or off (not expressed). In the domain of bioprocesses, it frequently happens that only a limited number or level of measurements are available (e.g. low, high, very high ...) because the devices only give a discretized semi-quantitative or qualitative measurement [6]. The measure may also be quantized by some physical device (as the time given by an analogical clock), and give as a result some number among a finite collection.
For this case of quantized outputs, the problem of control has to be considered in a non-classical way: the control cannot be a function of the full continuous state variables anymore, and most likely will change only when the quantized measurement changes. Moreover, the control itself could be quantized, due to physical device limitations. The above framework has been considered by numerous works, having their own specificity: quantized output and control with adjustable "zoom" and different time protocols, [20], hybrid systems abstracting continuous ones (cf. [17] for many examples of theories and applications).
In this paper, we consider a classical problem in the field of bioprocesses: the stabilization of an unstable bioreactor model, representing for example anaerobic digestion, towards a working set point. Anaerobic digestion is one of the most employed process for (liquid) waste treatment [10]. Considering a simplified model with two state variables (substrate and biomass), the system has two stable equilibria (and an unstable one, with a separatrix between the two basins of attractions of the respective stable equilibria), one being the (undesirable) washout of the culture ( [15]). The goal is to globally stabilize the process toward the other locally stable reference equilibrium. The (classical) output is the biomass growth (through gaseous production), the control is the dilution rate (see [23] for a review of control strategies). There exists many approaches based on well-accepted models ( [7]), for scalar continuous output [18], [3] or even for delayed and piecewise constant measurements [19]. Some original approaches make use of a supplementary competitor biomass [21].
In this paper, we suppose that the outputs are discrete or quantized: there are available in the form of finite discrete measurements. The precise measurements models are described later: roughly, the simplest one is "perfect", without noise, meaning that the true value is supposed to be in one of the discrete measurements, and that the transitions between two contiguous discrete measures are perfectly known. The next model is an uncertain model where the discrete measurements may overlap, and the true value is at the intersection between two quantized outputs. Remark that the model of uncertainty is different from the interval observers approaches ( [2], [13]) for the estimation or regulation [1], where some outputs or kinetics are not well known, but upper or lower bounds are known. Moreover, in the interval observer case, the variables are classical continuous variables.
For this problem, the general approaches described above do not apply, and we have to turn to more tailored methods, often coming from the theory of hybrid systems, or quantized feedbacks (see above) ... We here develop our adapted "hybrid" approach. It has also some relations with the fuzzy modeling and control approach: see e.g. in a similar bioreactor process the paper [11]. We provide here a more analytic approach, and prove our results of stability with techniques coming from differential inclusions and hybrid systems theory [12]. Our work has some relations with theoretical qualitative control techniques used for piecewise linear systems in the field of genetic regulatory networks ( [9]). The approach is also similar to the domain approaches used in hybrid systems theory, where there are some (controlled) transitions between regions, forming a transition graph [5], [14].
The paper is organized as follows: Section II describes the bioreactor model and the measurements models. The next section is devoted to the model analysis in open-loop (with a constant dilution rate). In Section IV, we propose a control law and show its global stability through the analysis of the transitions between regions (with the help of the Filippov definition of solutions of differential inclusion). Finally, in Section V, we explain how to choose the dilution rates by a graphical approach, and we end by giving some simulations, with or without uncertainty.

A. Model presentation
In a perfectly mixed continuous reactor, the growth of biomass x limited by a substrate s can be described by the following system (see [4], [10]): where s in is the input substrate concentration, u(t) the dilution rate, k the pseudo yield coefficient, and µ(s) the specific growth rate. Given ξ = (s, x), let us rewrite System (1) aṡ ξ = f (ξ, u(t)), where the dilution rate u(t) is the manipulated input.
The specific growth rate µ(s) is assumed to be a Haldane function (i.e. with substrate inhibition) [7]: whereμ, k S , k I are positive parameters. This function admits a maximum for a substrate concentration s = √ k S k I :=s, and we will assumes < s in . Proof: It is easy to check that the solutions stay positive. Now consider z = s + kx whose derivative writeṡ It follows that z is upper bounded by max(z(0), s in ), and so is kx. Finally, if s(t) > s in , thenṡ(t) < 0, therefore s is upper bounded by max(s(0), s in ).
In the following, we will assume initial conditions within the interior of the positive orthant.

B. Quantized measurements
We consider that a proxy of biomass growth y(ξ) = αµ(s)x is monitored (e.g. through gas production), but in a quantized way, in the form of a more or less qualitative measure: it can be levels (high, medium, low...) or discrete measures. Finally, we only know that y(ξ) is in a given range, or equivalently that ξ is in a given region (parameter α is a positive yield coefficient): where 0 = y 1 < y 2 < . . . < y n and y 1 < y 2 < . . . < y n−1 . We will consider two cases: • (A1). Perfect quantized measurements: This corresponds to the case where there is no overlap between regions. The boundaries are perfectly defined and measured.
. Uncertain quantized measurements: In this case, we have overlaps between the regions. In these overlaps, the measure is not deterministic, and may be any of the two values.
and (closed) switching domains: For perfect measurements (A1), we haveỸ i = intY i , and the switching domains Y i|i+1 correspond to the lines y(ξ) = y i = y i+1 . For uncertain measurements (A2), the switching domains Y i|i+1 become the regions {ξ ∈ R 2 + : Unless otherwise specified, we consider in the following uncertain measurements.

C. Quantized control
Given the risk of washout, our objective is to design a feedback controller that globally stabilizes System (1) towards a set-point. Given that measurements are quantized, the controller should be defined with respect to each region: Here D i is the positive dilution rate in region i. This control scheme leads to discontinuities in the vector fields. Moreover, in the switching domains, the control is undetermined. Thus, solutions of System (1) under Control law (3) are defined in the sense of Filippov, as the solutions of the differential inclusion [12]:ξ ∈ H(ξ) where H(ξ) is defined on regular domainsỸ i as the ordinary function H(ξ) = f (ξ, D i ), and on switching domains Y i|i+1 as the closed convex hull of the two vector fields in the two domains i and i + 1: Following [8], [9], a solution of System (1) under Control We first consider System (1) when a constant dilution rate D is applied in the whole space (i.e. u(t) = D, ∀t ≥ 0). In this case, the system is a classical ordinary differential equation which can present bistability, with a risk of washout. Let us denote s a (D) and s b (D) the two solutions, for D ∈ (0, µ(s)), of the equation µ(s) = D, with 0 < s a (D) <s < s b (D). For the Haldane growth rate defined by (2), we have:

III. MODEL ANALYSIS WITH A CONSTANT DILUTION
The asymptotic behavior of the system can be summarized as follows: k , see Figure 1.
Proof: See [15]. For D ∈ (0, µ(s)), let us define: y a (D) and y b (D) are the growth proxy obtained respectively at the equilibria ξ a (D) and ξ b (D) (if it exists) 1 .
In order to design our control law, we need to provide some further properties of the system dynamics. In particular, we need to characterizeẏ(ξ), the time derivative of y(ξ) along a trajectory of System (1) with a constant dilution rate u(t) = D: Let consider the following functions: defined respectively on (0,s) ∪ (s, +∞) and (0, +∞).
In the (s, x) plane, g D (s), h a D (s) and h b D (s) represent respectively the nullclineẏ(ξ) = 0 and the isolines y(ξ) = y a (D) and y(ξ) = y b (D) (i.e. passing through the equilibria ξ a (D) and ξ b (D)), see Figure 1. Knowing that the nullclinė y(ξ) = 0 is tangent to the isoline y(ξ) = y a (D) (resp. y(ξ) = y b (D)) at the equilibrium point ξ a (D) (resp. ξ b (D)), we will determine in the next lemma the relative positions of these curves, see Fig. 1.
Proof: See Appendix. This allows us to determine the monotonicity of y(ξ) in a region of interest (for the design of the control law).
Proof: See Appendix.

A. Control design
Our goal is to globally stabilize the system towards an equilibrium with a high productivity (the productivity is the output αµ(s)x), corresponding to a high dilution rate D n where there is bistability in open loop (case (ii) of Proposition 1). Number n is the number of measurements (see section II-B) and the control is such that D 1 < D 2 < . . . < D n . We consider the following control law, based on the quantized measurements y(ξ), and constant within a given region: given that the following conditions are fulfilled: These conditions make the equilibrium ξ a (D n ) globally stable, as we will see below. In Section V-A, we will precise how to choose the D i such that these conditions hold. In order to prove the asymptotic behavior of System (1) under Control law (4-7), the study will be divided into three steps: • the dynamics in one region, • the transition between two regions, • the global dynamics. This approach is similar to those deducing the global dynamics from a "transition graph" of possible transitions between regions [5].

B. Dynamics in one region with a given dilution: exit of domain
We first focus on a region Y i , i < n. A constant dilution D i -such that Conditions (5-6) for i hold -is applied. These conditions guarantee that the stable operating equilibrium for this dilution (see Proposition 1) is located in an upper region Y j , j > i, while the saddle point is located in a lower region Y k , k < i. This allows us to establish the following lemma: Proof: Let us consider the function V (ξ) = y a (D n ) − y(ξ) on Y i . Given Conditions (5-6), we get ∀ξ ∈ Y i : Since a constant dilution rate D i is applied, we can apply Lemma 3 to conclude thatẏ(ξ) > 0. Thus, V (ξ) is decreasing on Y i . Recalling that the trajectories are also bounded (Lemma 1), we can apply LaSalle invariance theorem [16] on the domain Ω 1 := {ξ ∈ Y i | x ≤ max(z(0), s in )/k, s ≤ max(s(0), s in )}. Given that the set of all the points in Ω 1 whereV (ξ) = 0 is empty, any trajectory starting in Ω 1 will leave this region. The boundaries x = max(z(0), s in )/k and s = max(s(0), s in ) are repulsive (see Proof of Lemma 1). Finally, the boundary y(ξ) = y i corresponds to the maximum of V (ξ) on Ω 1 , so every trajectory will leave this domain, crossing the boundary y(ξ) = y i .

C. Transition between two regions
Now we will characterize the transition between regions (as we have seen above, the intersection can be either a simple curve in the case of perfect measurements, or a region with non empty interior in the uncertain case): Proof: First, we consider i = n − 1. We will follow the same reasoning as for the previous lemma, applying LaSalle theorem on a domain Similarly, V (ξ) is also decreasing on Y i+1 whenever u = D i+1 . Now under Control law (4), we have shown that V (ξ) is decreasing on the regular domainsỸ i andỸ i+1 . V (ξ) is a regular C 1 function, and can be differentiated along the differential inclusion. On the switching domains Y i|i+1 , we have:V Thus, V (ξ) is decreasing on Ω 2 . Following the proof of Lemma 4 concerning the boundaries, we can deduce that every trajectory will reach the boundary y(ξ) = y † , i.e. it will enter Y i+1 . For i = n − 1, taking any y † ∈ (y n−1 , y a (D n )), we can show similarly that V (ξ) is decreasing on Ω 2 so every trajectory will enter the regionỸ n . Following the same proof, we can show that the reverse path is not possible, in particular for the last region: Lemma 6. Consider System (1) under Control law (4) with Conditions (5) for n − 1, n, Condition (6) for n − 1, and Condition (7). The regular domainỸ n is positively invariant.

D. Global dynamics
Now we are in a position to present the main result of the paper:

Proposition 2. Control law (4) under Conditions (5-7) with perfect or uncertain measurements (A1 or A2) globally stabilizes System (1) towards the point ξ a (D n ).
Proof: From Lemmas 5 and 6, we can deduce that every trajectory will enter the regular domainỸ n , and that this domain is positively invariant.
System (1) under a constant control u(t) = D n has two non-trivial equilibria (see Proposition 1): ξ a (D n ), and ξ b (D n ). The growth proxy at these two points satisfy y a (D n ) > y n−1 and y b (D n ) < y n (Conditions (5,7)), so there is only one equilibrium inỸ n : ξ a (D n ). Moreover, it is easy to check thatỸ n is in the basin of attraction of ξ a (D n ), therefore all trajectories will converge toward this equilibrium.

A. How to fulfill Conditions (5-7)
The global stability of the control law is based on Conditions (5-7). We now wonder how to easily check if these conditions hold, or how to choose the dilution rates and/or to define the regions in order to fulfill these conditions. In this purpose, a graphical approach can be used.
As an example, we will consider the case where the regions are imposed (by technical constraints) and we want to find the different dilution rates D i such that Conditions (5-7) hold. Our objective is to globally stabilize the equilibrium point ξ a (D * ) ∈Ỹ n , with µ(s in ) < D * < µ(s ⋄ ) (s ⋄ is defined just after). Note that we impose D * < µ(s ⋄ ) given that for any D * > µ(s ⋄ ), the same productivity can be achieved with a smaller dilution rate, leading to a reduced risk of instability.
Let us denote s c (y) and s d (y) the two solutions, for y ∈ (0, φ(s ⋄ )), of the equation φ(s) = y, with 0 < s c (y) < s ⋄ < s d (y) < s in . For the Haldane growth rate, we have: This can be done analytically or graphically, as shown on Figure 2A. Whenever we choose D i ∈ (D(y i ), D(y i )), we have  (6)). so Conditions (5-6) for i hold.
If D(y i ) ≥ D(y i ) (see Figure 2B), it is not possible to fulfill the conditions and thus to implement the control law with this measurement range.
How to choose all the D i : The procedure proposed in the previous subsection should be repeated for all the regions. We can depict two particular cases: • For the first region Y 1 , given y 1 = 0, we actually impose D 1 < µ(s in ). • For the last region Y n , one should check that D(y n−1 ) < D n = D * < D(y n ). This approach is illustrated on Figure 3.

Increasing measurement resolution:
We have seen that for a given region, it is not always possible to fulfill Conditions (5)(6)(7). This gives rise to a question: is there any constraint on the measurement that guarantees the possibility to implement the control law?
For perfect measurements (A1) with equidistribution, we will show that it will always be possible to implement the control law increasing the number of regions.
First, we can arbitrarily define the lower bound of the last region: y n ∈ (y b (D * ), y a (D * )), recalling nonetheless that n is unknown. Now, we will determine the limit range of a measurement region. Letŷ such that D(ŷ) = µ(s in ). For y i < y i <ŷ, we have D(y i ) < µ(s in ) < D(y i ) so it is always possible to choose a D i in order to fulfill Conditions (5)(6).
For y ∈ [ŷ, y n ], we define (see Figure 2C): It gives a lower bound on y i : whenever one choose y i > ψ(y i ), then we have: so one can choose any D i ∈ (D(y i ), D(y i ) and Conditions (5-6) for i will hold. We can now define the mapping ∆ : y → y − ψ(y) on [ŷ, y n ]. ∆(y) defines the maximal range of the region with upper bound y (in order to be able to find a dilution rate such that Conditions (5-6) hold). ∆ admits a minimum: Thus, for n > y n ∆m + 1, the regions Y i defined by: y n , i = 1, ..., n − 1 allow the implementation of the control law. In conclusions, whenever the measurement resolution is good enough (i.e. the number of regions is high enough), it is always possible to find a set of dilution rates D i such that Conditions (5-7) hold in the case of perfect measurements (A1) with equidistribution. For uncertain measurements (A2), the same result holds, but the proof is omitted for sake of brevity.

B. Simulations
As an example, we consider the anaerobic digestion process, where the methane production rate is measured. Parameters, given in Table I, are inspired from [7] (considering only the methanogenesis step). For uncertain measurements, we use discrete time simulation. At each time step t k (with ∆t = 0.05d), when ξ(t k ) is in a switching region Y i|i+1 , we choose randomly the control u(t k ) between D i and D i+1 . In this case, we perform various simulations for a same initial condition. Our objective is to stabilize the equilibrium ξ a (D * ), with D * = 0.47 d −1 (which corresponds to a productivity of 92% of the maximal productivity). We first consider the following perfect measurement set (with equidistant region): (6) with y n = 4 L CH 4 .L −1 .d −1 . With four regions (n = 4), we can define a set of dilution rates such that Conditions (5-7) are fulfilled (see Fig. 3): For uncertain measurements, we increased each upper bound and decreased each lower bound by 10%. It appears that the same dilution rates can be chosen.
Trajectories for various initial conditions are represented in the phase portrait for perfect and uncertain measurements, see Fig. 4. In accordance with Proposition 2, all the trajectories converge towards the set-point. Thus, the transition graph is deterministic (there is only one transition from a region to the upper one), see Fig. 5a.
If the number of regions is reduced (three regions only), it is not possible in this example to choose dilution rates such that Conditions (5-7) hold for all i. In this case, some trajectories do not converge towards the set-point. Some regions have transitions towards the upper region, but also towards the lower one. There are sliding modes. This aspect will be further discussed in the next subsection.

C. When conditions are not verified: risk of failure
We here detail what happens if Conditions (5-7) are not fulfilled, and in particular if there is a risk of washout. This point is illustrated by Fig. 6 and Fig. 7.
First, given the previous analysis of the system, one can easily see that only the condition y b (D 1 ) < y 1 , i.e. D 1 < µ(s in ) is necessary to prevent a washout, so D 1 can be chosen with a safety margin in order to avoid such situation. Now, if Condition (5) is not fulfilled for some i > 1, the unstable equilibrium ξ b (D i ) will be located in the region Y i . Thus, the regionỸ i have transitions towards the lower region, and a trajectory can stay in the switching domain Y i−1|i . Given that z = s + kx converges towards s in (cf. proof of Lemma 1), such trajectory will converge towards the intersection between the switching domain and the invariant manifold z = s in (see Fig. 5b and Fig. 6): • For perfect measurements (A1), this gives rise to a sliding mode and the convergence towards a singular equilibrium point. • For uncertain measurements (A2), all the trajectories converge towards a line segment. In our simulation, they actually also converge towards a singular equilibrium point. This situation can be detected by the incessant switches between two regions. In such case, the dilution rate D i should be slightly decreased.
On the other hand, if Condition (6) does not hold for some i, the stable equilibrium ξ a (D i ) will be located in the region Y i , so some trajectories can converge towards this point instead of going to the next region, see Fig. 5c and Fig. 7. To detect such situation, a mean escape time for each region can be estimated using model simulations. A trajectory which stay much more than the escape time in one region should have reached an undesirable equilibrium. In this case, the dilution rate D i should be slightly increased.
In all the cases, the trajectories converge towards a point or a line segment. Although it is not desired, this behavior is particularly safe (given that there is theoretically no risk of washout). Moreover, as explained above, these situations can easily be detected and the dilution rates can be changed accordingly (manually or through a supervision algorithm).

VI. CONCLUSION
Given the quantized measurements, we were able to design (under some conditions) a control based on regions and transition between regions. These tools are similar to the ones of piecewise linear systems, and it is possible to draw a transition graph showing all the possible transitions. Moreover, we have seen that for some undesirable cases, singular behaviors (sliding modes) are possible on the boundaries between regions. We think that this kind of control on domains, and the design of the resulting transition graph, is a promising approach, that we want to deepen in future works. This approach could be generalized to other classical systems, e.g. in mathematical ecology. Given that ϕ b D (s) is continuous on (s, +∞), we finally conclude that on this interval, ϕ b D (s) ≤ 0 , i.e. g D (s) ≤ h b D (s).