www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Model-based control strategies for anaerobic digestion processes Neli Dimitrova∗, Mikhail Krastanov†∗ ∗Institute of Mathematics and Informatics Bulgarian Academy of Sciences Sofia, Bulgaria †Faculty of Mathematics and Informatics Sofia University “St Kl. Ohridski” nelid@math.bas.bg, krastanov@fmi.uni-sofia.bg Received: 28 February 2019, accepted: 12 July 2019, published: 17 August 2019 Abstract—In this paper we consider a four- dimensional bioreactor model, describing an anaer- obic wastewater treatment with methane produc- tion. Different control strategies for stabilizing the dynamics are presented and discussed. A general and practice-oriented bounded open-loop control is proposed, aimed to steer the model solutions towards an a priori given set in the phase plane. Keywords-continuous bioreactor; dynamical non- linear model; global asymptotic stabilization; feed- back control; open-loop control; model uncertainty I. INTRODUCTION Anaerobic digestion (AD) is a biological pro- cess in which organic degradable material is con- verted into biogas by microorganisms [7], [8]. Recently, AD has been evaluated as one of the most promising processes for waste recovery, en- vironmental protection and bioenergy production. The biogas is a mixture of gases composed of methane, carbon dioxide, nitrogen, oxygen, hydro- gen sulphide and traces of other gases. The biogas is classified as a renewable energy, which can be used in gas engines to produce electricity and heat energies. Storing biogas prevents greenhouse gas emissions from entering the atmosphere. Some estimates from 1997 [23] report that recovery of organic wastes and industrial effluents could reduce 20% of the global warming effect on the planet. At laboratory or industrial scales, the AD pro- cess occurs inside an anaerobic digester (reactor), where degradation of organic material holds by plenty of anaerobic microorganisms. The growth of the microorganisms (bacteria, yeasts, etc.) pro- ceeds by consumption of appropriate nutrients (substrates) involving carbon, nitrogen, oxygen, etc., under favorable conditions (temperature, pH, etc.). The mass of the living organisms (or cells) Copyright: c©2019 Dimitrova et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes, Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 1 of 17 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes is called biomass. The number, behavior and in- teraction of the included organisms pose a chal- lenge to the specialist in the field and have been extensively investigated in the literature. The re- actor configuration and environmental conditions (retention time, temperature, feedstock, stirring, etc.) influence the dynamics and composition of the different groups of bacteria responsible for the organic material degradation [7]. According to the user objectives and the nature of the biodegradable waste (solid wastes or released in wastewater), different technologies can be used. We mention below some commonly cited reactors in the liter- ature [7], [8]. Batch reactor is a reactor without inflow nor outflow. The digester is filled with the biodegrad- able materials and left until the substrate has been degraded. Then the digester is emptied and a new cycle can start again. Fed-bach reactor (also called semi-continuous or fed/sequencing batch reactor) is a reactor with- out outflow. The process is cyclic, the digester is filled gradually according to the progress of the reaction in order to ensure optimal growth conditions. At the end of the digestion phase, decantation allows to separate the liquid phase and suspended biomass. Continuous bioreactor: the tank is continuously fed at a constant rate and the digestat (the material that remains after the AD process) is evacuated by a mechanical action. Depending on the contact between the substrate and biomass or the feeding mode, the continuous bioreactors fall into several categories. Among them we mention the contin- uously stirred tank reactor, where the outflow is equal to the rate of inflow and a continuous mixing ensures the medium homogeneity. Independently of the chosen reactor type, a key parameter in biogas plants is the dilution rate. It is proportional to the speed of the input mechanisms which feeds the reactor with substrate. To avoid wash out of bacteria, the dilution rate is always constrained, i. e. u ∈ [umin,umax] with umin > 0 (cf. [8]). The AD process follows several phases: Hydrolysis is the step where polymers (macro- molecules) are hydrolysed to monomers (simple organic matter); the speed of degradation depends on the substrate type (glucide, proteins, lipids, cellulose, etc.). Acidogenesis, where monomers are degraded to Volatile Fatty Acids (VFA) and alcohol. Acetogenesis, performed by acitogenic bacte- ria which transform the VFA into acetic acid, hydrogen (H2) and carbon dioxide (CO2). The responsible bacteria for this step produce H2 and can be inhibited by an excess of H2 concentration in the digester, that’s why they live fixed to the methanogenic bacteria which consume the hydro- gen. Methanogenesis, where methanogenic bacteria reduce the specific substrate into methane. A good management and control of the AD pro- cess can be achieved via validated mathematical models—an area, which is extensively studied in recent years. The proposed models are specific to a couple of criteria such as the waste nature and its composition, the used technology, collected data and its quality, the possible experiments and changes in the operating conditions. In general, a dynamical model accounts for the time evolution of substartes and biomasses, and is based on ordinary differential equations representing mass balance within the process. The mathematical modeling of AD processes has a long history. In 1968, Andrews [5] mod- elled the methanogenic fermentation by only the final step methanogenesis. In 1973, Graef and Andrews [24] included the acidogenesis step in the description of fermentation. Later on, other re- searchers, Hill and Barth (1977, [28]), Boone and Bryant (1980, [12]), Eastman and Ferguson (1981, [20]) added a hydrolysis step to their description, and modelled a three-step process. The interested reader is referred to [32] for more details about other models. Over time the models were extended depending on the different substrates (wastewater, sludge, etc.). In 2002, a group of experts in the AD pro- cesses (IWA Task Group for Mathematical Mod- Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 2 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes elling of Anaerobic Digestion Processes) devel- oped a standard model for the AD process, called ADM1 (Anaerobic Digestion Model No 1) [9]. In order to make the ADM1 a standard platform for AD simulation, it has been decided to generalize the composition of waste, so it is measured by an unified unit Chemical Oxygen Demand (COD) and the process is supposed to occur in a continuously stirred tank reactor. ADM1 is described by 32 ordinary differential equations and its parameters are collected from different applications. Many modifications, adaptations and variations of the ADM1 have been done later, cf. [43], [44], [47], [48] and the references therein. The ADM1 and its variations are complex mod- els suitable for process knowledge and simulation, but not appropriate for process control and soft- ware sensors design, because they require a plenty of input parameters which are difficult to obtain. To overcome the ADM1 complexity, simpler models based on mass balance equations [8] have been developed, more suitable to support moni- toring or control strategies. Such a model (called AM1), including two reactions (acidogenesis and methanogenesis) is proposed by Bernard et al [11], and turns out to approximate efficiently the ADM1 for modeling anaerobic wastewater treatment. This model will be investigated in the present paper. The management and control of the AD process require a good information about the internal state of the system. Biological processes are known to be highly unstable due to the specific behavior of the system itself or to the presence of disturbances. Thus, an obvious need for an efficient control and monitoring of such systems arises. A summary and review of different sensoring approaches can be found e. g. in [31], [33], [42], [49]. The ma- jority of sensors intended to measure the process key variables often require complex equipment and careful maintenance, so that the plant costs may climb quickly, which is not desirable from industrial standpoint. Therefore, it is crucial to find a methodology which allows cost-effective and easily adopted to practice monitoring of AD plants. Such methodology is the development of efficient software sensors, also called observers. The observers are auxiliary dynamical systems that provide information on the unmeasurable state variables of the system by using its mathematical model and its input and output signals (the measur- able variables of the system), see e. g. [2], [3], [4], [8], [26], [40]. In biological processes observers are mainly used in on-line estimations for control purposes. In addition to the modeling and observer design, the AD control in biogas plants is gaining an increased importance. The main reason is the sig- nificant growth of bioenergy markets. Moreover, due to the climate-energy package of the European Commission, the produced biogas must be rich in methane to fulfil the environmental standards [51]. Controlling the AD processes is a delicate problem due to the high complexity and strong instability of the ecosystem inside the reactor [26]. Several factors are to be handled like e. g. the highly nonlinear behavior of the system itself, load disturbances, system uncertainties, constraints on manipulated and state variables and the limited on- line measurements information [39]. Moreover, the AD process involves living organisms which are very sensitive to the operating conditions and may be washed out or inhibited due to an accidental toxic feeding, leading, in the worst case, to a stop of the digester. The control design varies with the application objectives. Usually, in biogas plants, the con- troller is designed to satisfy one specific cri- terion, either economical (maximizing methane production) or ecological (minimizing COD con- centration of the effluent) or stability (VFA or dissolved hydrogen) criteria [21], [39]. The con- troller type depends on many factors such as the accuracy of the monitoring, knowledge of the system and availability and complexity of the considered model. Among the classical controllers are the proportional-integral (PI) controller, the proportional integral-differential (PID) controller, the adaptive PID and the cascade PI controls; all they have been recognized as a good alternative Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 3 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes for the regulation of AD plants (cf. [22], [40] and the references therein). Various advanced control approaches like expert systems (rule-based and fuzzy-logic-based systems, neural networks, etc.) have been recently developed for AD control. An other type of model-based control designed for the AD processes, is the linearizing control [8]. The latter is based on a nonlinear model, aimed to achieve linear closed-loop dynamics. A drawback of this method is that it relies of full knowledge of the system parameters. Later on, an adaptive linearizing control has been proposed [35], which ensures global asymptotic stability of the closed-loop system. When intervals of the model uncertainties are known a priori, robust linearizing control based on interval observers has been proposed in [1], [41] and [45]. Other recently developed approaches for controlling AD processes are based on differential geometry [29], [38]. Sliding mode approaches have been also used to control anaerobic continuous bioreactors. Further, the nonlinear adaptive control law, which is robust with respect to unknown kinetic rates has been proposed for the global stabilization of AD processes [35]. Extremum seeking control (ESC) is another technique to handle dynamic optimization problems. The goal of ESC is to find the operating set-point, a priori unknown, such that a perfor- mance function reaches its extremum value. The classical extremum seeking approach [36], [37], [50] is designed in the form of a block diagram (scheme) that is implemented on the bioreactor to tune the dilution rate of the open-loop system towards a set-point, where an optimal value of the output is achieved. The main limitation in applying this approach is that the dynamics should be open-loop stable. Otherwise, a local controller is necessary to stabilize the system around the optimal operating point. Extremum seeking numerical techniques are de- veloped in recent years to overcome the above drawbacks. Based on a mathematical model, this new approach splits the extremum seeking prob- lem into two steps: global dynamics stabiliza- tion and application of a numerical optimization method. In [46] this approach is applied to the classical two-dimensional model of methane fer- mentation. For further information about instrumentation and control of AD processes we refer the reader to the excellent review [30]. The present paper is organized as follows. Section II presents the dynamic model of the anaerobic wastewater treatment process and gives an overview of authors’ results on global stabi- lizability of the model dynamics using different control strategies. Section III reports on a new re- sult, concerning a general and practically oriented stabilization approach by means of an arbitrary measurable bounded control function. II. BASIC PROPERTIES AND GLOBAL STABILIZABILITY OF THE MODEL DYNAMICS We consider the mass balance model AM1 of anaerobic wastewater treatment in a continuous bioreactor, described by the following nonlinear system of ordinary differential equations (cf. e. g. [6], [11], [25], [27], [34]): ds1 dt = u(si1 −s1) −k1µ1(s1)x1 dx1 dt = (µ1(s1) −αu)x1 ds2 dt = u(si2−s2)+k2µ1(s1)x1−k3µ2(s2)x2 dx2 dt = (µ2(s2) −αu)x2. (1) The definition of the state variables s1, s2 and x1, x2 as well as of the model parameters is given in Table I. All coefficients are assumed to be positive. The parameter α ∈ (0, 1) represents the proportion of bacteria that are affected by the dilution; α = 0 and α = 1 correspond to a fed-batch reactor and to a continuously stirred tank reactor, respectively (cf. [2], [6], [11], [25]). The input substrate concentrations si1 and s i 2 are assumed to be constant. The dilution rate u is considered as a control (manipulated) input. The model describes two steps of the AD pro- cess: Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 4 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes TABLE I DEFINITION OF THE MODEL PHASE VARIABLES AND PARAMETERS s1 concentration of chemical oxygen demand (COD) [g/l] s2 concentration of volatile fatty acids (VFA) [mmol/l] x1 concentration of acidogenic bacteria [g/l] x2 concentration of methanogenic bacteria [g/l] u dilution rate [day−1] si1 influent concentration s1 [g/l] si2 influent concentration s2 [mmol/l] k1 yield coefficient for COD degradation [g COD/(g x1)] k2 yield coefficient for VFA production [mmol VFA/(g x1)] k3 yield coefficient for VFA consumption [mmol VFA/(g x2)] k4 coefficient [l2/g] Q methane gas flow rate m1 maximum acidogenic biomass growth rate [day−1] m2 maximum methanogenic biomass growth rate [day−1] ks1 saturation parameter associated with s1 [g COD/l] ks2 saturation parameter associated with s2 [mmol VFA/l] kI inhibition constant associated with s2 [(mmol VFA/l)1/2] (i) acidogenesis, where the organic substrate s1 is degraded into volatile fatty acids (VFA) (s2) by acidogenic bacteria (x1); (ii) methanogenesis, where VFA (s2) are degraded by methanogenic bacteria (x2) into methane CH4. The methane solubility is very low, therefore the methane produced by the second step is not stored in the liquid phase, thus the output methane flow rate Q = QCH4 is written in the form Q = k4µ2(s2)x2. (2) The model dynamics can be described schemat- ically by the following biological reaction path- ways: Acidogenesis: k1s1 r1(·)−→ x1 + k2s2 Methanogenesis: k3s2 r2(·)−→ x2 + k4Q, where rj(·), j = 1, 2, are the reaction rates, which are given by rj(·) = µj(·)xj. The functions µ1(s1) and µ2(s2) model the spe- cific growth rates of the microorganisms. Usually, the model is studied using the following particular expressions of µ1 and µ2 (cf. [2], [6], [11], [25], [27]) µ1(s1) = m1s1 ks1 + s1 (Monod law) µ2(s2) = m2s2 ks2 +s2 + ( s2 kI )2 (Haldane law), (3) where m1, ks1 , m2, ks2 and kI are positive coef- ficients (see Table I). The most crucial problem in investigating AD models is the formulation of reasonable analytic expressions for µ1(s1) and µ2(s2). In our theoret- ical studies on the model (1) we do not assume to know explicit expressions for µ1 and µ2, we only impose the following general assumption on the latter. Assumption A1. The function µj(sj) is defined for sj ∈ [0, +∞), µj(0) = 0, µj(sj) > 0 for sj > 0; µj(sj) is continuously differentiable and bounded for all sj ∈ [0, +∞), j = 1, 2. The model (1) exhibits very rich dynamics. Con- sidering u as a bifurcation parameter, the system possesses different types of bifurcations of the equilibrium points, most of them leading to wash- out of the biomass and thus to process break-down [10], [14]. Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 5 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Define s = k2 k1 s1 + s2, s i = k2 k1 si1 + s i 2. (4) The quantity s is called biological oxygen de- mand (BOD) and represents the biological equiv- alent of COD, i. e. the biological equivalent of the total amount of organic substrate in the digester. For the practical application it is worth to note that BOD is online measurable and is used as a depollution factor in wastewater treatment, see e. g. [2], [6], [11], [13] and the references therein. Define the set Ω0 = { (s1,x1,s2,x2) ∈ R4 : s1 > 0, x1 > 0, s2 > 0, x2 > 0} . (5) The basic properties of the model solutions are summarized in the next lemma, which extends assertions that are given in [25], [53] and [54]. Lemma 1. Let Assumption A1 be fulfilled. Then for each point p0 = (s01,x 0 1,s 0 2,x 0 2) ∈ Ω0 the corresponding solution (s1(t),x1(t),s2(t),x2(t)) of (1) with (s1(0),x1(0),s2(0),x2(0)) = p0 is defined for each t > 0. Moreover, for each ε > 0 there exists Tε such that for each t > Tε the following inequalities hold true: si −ε < s(t) + k3x2(t) < si α + ε, si1 −ε < s1(t) + k1x1(t) < si1 α + ε, s1(t) < s i 1 and x1(t) ≥ ε. Below we present some authors’ results related to global stabilizability of the model dynamics by means of different control strategies and satisfying different criteria—the ecological criterion in sub- sections A and B or the economical criterion in subsection C. A. Global stabilizability via input control Here we investigate the global stabilizability of the dynamics (1) using the classical approach, where the dilution rate u is considered as a con- trol variable. More precisely, we show that for any admissible value of u there exists a nontriv- ial (with positive components) equilibrium point, which is globally asymptotically stable for the system. Although the manipulated input u is the most exploited variable for control purposes, to the authors’ knowledge there is no rigorous proof so far in the literature for global stabilizability of this model. Assume that the control variable u varies in the interval u ∈ (0,ub), where ub ≤ 1 α min { µ1(s i 1),µ2(s i 2) } ≤ umax. Let for some ū ∈ (0,ub) the following assump- tion holds true: Assumption A2. There exist points s1(ū) = s̄1 ∈( 0,si1 ) and s2(ū) = s̄2 ∈ ( 0,si2 ) , such that the following equalities hold true ū = 1 α µ1(s̄1) = 1 α µ2 (s̄2) . Assumption A2 is called regulability [25] of the system: it ensures the existence of a nontrivial equilibrium of the model (1), corresponding to the chosen value of the dilution rate u. Determine the points s̄1 and s̄2 according to Assumption A2 and compute further x1(ū) = x̄1 = si1 − s̄1 αk1 , x2(ū) = x̄2 = si2 − s̄2 + αk2x̄1 αk3 . Then the point p(ū) = p̄ = (s̄1, x̄1, s̄2, x̄2) is an equilibrium point of the system (1). In practical applications, the chosen equilibrium point is also called an operating or a reference point. Let s and si be determined according to (4). The next assumption is Assumption A3. The following inequalities hold true: (i) µj(sj) < µj(s̄j) for sj ∈ (0, s̄j),j = 1, 2; (ii) µ1(s1) > µ1(s̄1) for s1 ∈ (s̄1,si1); (iii) µ2(s2) > µ2(s̄2) for s2 ∈ (s̄2,si). Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 6 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Assumption A3 is technical. It is always fulfilled when the functions µj(·), j = 1, 2, are monotone increasing (like the Monod specific growth rate). If µj(·) is not monotone increasing (like e. g. the Haldane law) then ū has to be chosen in a proper way, such that Assumption A3 will be satisfied. For biological evidence s2 satisfies the inequal- ity s2 ≤ si2. The requirement s2 < s i in As- sumption A3(iii) is motivated by the fact that si can be considered as the worst-case upper bound of the concentration s2 due to imbalance between acidogenesis and methanogenesis, leading to acid- ification (x2 = 0) in the bioreactor (cf. [27]). Let ū ∈ (0, ub) be chosen according to As- sumptions A2 and A3. Denote by Σ1 the system obtained from (1) by substituting the control vari- able u by ū. Then the following Theorem 1 reports on the global stability of the equilibrium point p̄. Theorem 1. (cf. [18]) Let the Assumptions A1, A2 and A3 be fulfilled and let p0 = (s01,x 0 1,s 0 2,x 0 2) be an arbitrary point from the set Ω0. Then the solution of Σ1 starting from the point p0 converges asymptotically towards p̄. A drawback of this control approach is that it requires exact and full knowledge of the system states and precise tuning of u to achieve the global stabilizability. This drawback will be overcome by applying an output feedback control. B. Global stabilizability via output feedback con- trol This subsection is devoted to global stabilizabil- ity of the dynamics (1) by means of a feedback control and in the presence of model uncertainties. The proposed state feedback is of the form u ≡ κ(s2,x2) = βk4µ2(s2)x2, where β is a properly chosen positive parameter. The proposed feedback law is strongly related to the output (2). The parameter β provides some degrees of freedom in choosing an admissible reference (equilibrium) point, usually determined by ecological rules. This fact is also exploited in designing an extremum seeking algorithm for maximizing the methane production in real time (cf. subsection C below). Consider the dynamical system (1) in the state space ζ = (s1,x1,s2,x2). Using the definition of si from (4) we make the following assumption: Assumption A4. Lower bounds si− and k−4 for the values of si and k4 respectively, as well as an upper bound k+3 for the value of k3 are known. Define the following feedback control law: κ(ζ) = β k4 µ2(s2) x2 (6) with β ∈ ( k+3 si− ·k−4 , +∞ ) . The feedback control law κ(·) can be written in the form κ(·) = β ·Q(·), where Q is the methane output (2). For the practical applications it is worth to note that Q is on-line measurable, so this holds true for the feedback control κ(·) as well. Denote by Σ2 the closed-loop system obtained from (1) by substituting the control variable u by the feedback κ(ζ) from (6). Choose some β ∈ ( k+3 si− ·k−4 , +∞ ) and let ξ̄ = si − k3 βk4 ; obviously, ξ̄ belongs to the interval (0,si). The next assumption is similar to the regulability As- sumption A2. Assumption A5. There exists a point s̄1 such that µ1(s̄1) = µ2 ( ξ̄ − k2 k1 s̄1 ) > 0, s̄1 ∈ ( 0,si1 ) . Find s̄1 according to Assumption A5 and define s̄2 = ξ̄ − k2 k1 s̄1, x̄1 = si1 − s̄1 αk1 , x̄2 = 1 αβk4 . It is straightforward to see that the point p̄β = (s̄1, x̄1, s̄2, x̄2) is an equilibrium point of Σ2. We shall show below that the feedback law (6) stabilizes asymptotically the closed-loop system towards p̄β (cf. [16]). Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 7 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Let Ω0 be defined according to (5). Using the definitions of s and si from (4), we define the sets Ω1 = {(s1,x1,s2,x2) : s1 + k1x1 ≤ si1/α, s + k3x2 ≤ s i/α } , Ω2 = {( s1,x1, ξ̄ − k2 k1 s1, x̄2 ) : 0 < s1 < k1 k2 ξ̄, x1 > 0 } , Ω = Ω0 ∩ Ω1. Assumption A6. Let the inequality µ′1(s̄1+θ(s1−s̄1))+ k2 k1 ·µ′2 ( s̄2−θ k2 k1 (s1−s̄1) ) >0 be satisfied for each s1, belonging to the projection of the set Ω ∩ Ω2 on the s1-axis and for each θ ∈ [0, 1]. Assumption A6 is technical. It is always fulfilled when the functions µj(·), j = 1, 2, are monotone increasing. If µj(·) is not monotone increasing then the set-point ξ̄ (or equivalently the value for β) has to be chosen in a proper way in order to satisfy Assumption A6. Theorem 2. (cf. [16]) Let the Assumptions A1, A4, A5 and A6 be satisfied. Let us fix an ar- bitrary number β ∈ ( k+3 si− ·k−4 , +∞ ) and let p̄β = (s̄1, x̄1, s̄2, x̄2) be the corresponding equi- librium point. Then the feedback control law κ(·) defined by (6) stabilizes asymptotically the control system Σ2 to the point p̄β for each starting point ζ0 = (s01,x 0 1,s 0 2,x 0 2) ∈ Ω0. C. Extremum seeking control Consider the equation (2) describing the process output, i. e. the methane production. A model- based numerical extremum seeking algorithm is proposed in authors’ publications [14]–[18] to steer and stabilize the dynamics (1) towards a steady state, where maximum methane flow rate Qmax is achieved. For that purpose the function Q is computed on the set of all equilibrium points, parameterized with respect to: (i) u in the case of the input control (Theorem 1), (ii) β in the case of the output feedback law (Theorem 2). Denote the so obtained function by Q(q), where q ∈ (q−,q+) denotes one of u or β, and q− > 0, q+ > 0 are the corresponding bounds accord- ing to Theorem 1 or Theorem 2. The function Q(q) is called input-output static characteristic of the model. Assume that Q(q), q ∈ (q−,q+), is strongly unimodal, i. e. there exists a unique point qmax ∈ (q−,q+) where Q(q) takes a maximum, Qmax = Q(qmax), the function strictly increases in the interval (q−,qmax) and strictly decreases in (qmax,q +). Denote by E(q) the equilibrium point parame- terized on q and let E(qmax) be the steady state where Qmax is achieved. The goal is to stabilize in real time the systems Σ1 and Σ2 towards this (a priori unknown) equilibrium point E(qmax) and therefore to the maximum methane flow rate Qmax. This is realized by applying a numerical model-based extremum seeking algorithm. The main idea of the algorithm is the following: a sequence of points q1,q2, . . . ,qn, . . . from the interval (q−,q+) is constructed, each qj being in the form qj = qj−1 ± hj, hj > 0, and such that {qj} tends to qmax; Theorems 1 and 2 guaran- tee that the dynamics is globally asymptotically stabilizable towards the equilibrium E(qj), j = 1, 2, . . .. Then by computing and comparing the values Q(q1),Q(q2), . . . ,Q(qn), . . ., the desired equilibrium point E(qmax) and thus Qmax are achieved. In the computer implementation the algorithm is carried out in two stages. In the first stage, “rough” intervals [q] and [Q] are found which enclose qmax and Qmax respectively; in the second stage, the interval [q] is refined using an elimination proce- dure based on the golden mean value (or Fibonacci search) strategy. The second stage produces the final intervals [qmax] = [q−max,q + max] and [Qmax] such that qmax ∈ [qmax], Qmax ∈ [Qmax] and q+max − q−max ≤ �, where the tolerance � > 0 is specified by the user. III. OPEN–LOOP CONTROL STABILIZATION The previous Section II was devoted to global stabilization of the dynamics (1) towards a previ- ously chosen equilibrium (operating) point. This Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 8 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes section proposes a different approach for stabiliz- ing the model solutions. Instead to an equilibrium point, the goal here is to steer the model solutions so that the BOD values s fall onto an interval [S−,S+], given a priori by ecological rules, and remain there for all time. This is achieved by suitably constructed bounded open-loop control. Assumption A7. Let the functions µj(sj) are strictly increasing in the intervals (0,sij), j = 1, 2. Assumption A7 is technical. It is always fulfilled when the functions µj(sj), j = 1, 2, are presented by the Monod specific growth rate. Let s and si be defined according to (4), i. e. s = k2 k1 s1 + s2, s i = k2 k1 si1 + s i 2. Let us choose an operating (reference) point s∗, s∗ ∈ (0,si). Assumption A8 (regulability). There exists a point s∗1 such that µ1(s ∗ 1) = µ2 ( s∗ − k2 k1 s∗1 ) > 0, s∗1 ∈ ( 0,si1 ) . Find s∗1 according to Assumption A8 and deter- mine further s∗2 = s ∗ − k2 k1 s∗1, x ∗ 1 = si1 −s ∗ 1 αk1 , x∗2 = si2 −s ∗ 2 + αk2x ∗ 1 αk3 = si −s∗ αk3 . It is straightforward to see that the point ζ∗ = (s∗1,x ∗ 1,s ∗ 2,x ∗ 2) is an equilibrium point of system (1). One can directly check that the equilibrium points of (1) satisfy the equalities (cf. the straight line l2 in Fig. 1) s1 + αk1x1 = s i 1, s + αk3x2 = s i. Denote u∗ = 1 α µ1(s ∗ 1) = 1 α µ2(s ∗ 2). Practically, ecological norms prescribe an ad- missible interval [S−,S+] for the BOD values s. We choose an arbitrary interval [s−,s+] con- tained in the interior of [S−,S+], i. e. [s−,s+] ⊂ (S−,S+). Assumption A9. The positive reals s−1 , s + 1 , s − 2 , s+2 , u − and u+ satisfy the relations s− = k2 k1 s−1 + s − 2 , s + = k2 k1 s+1 + s + 2 , αu− = µ1(s − 1 ) = µ2(s − 2 ), αu+ = µ1(s + 1 ) = µ2(s + 2 ), 0 < s−1 < s ∗ 1 < s + 1 < s i 1, 0 < s−2 < s ∗ 2 < s + 2 < s i 2, u− < u∗ < u+, and each point of the interval [u−,u+] is an admissible value for the dilution rate u. The imposed boundedness of u in Assumption A9 is not restrictive. Practically, the dilution rate is associated with the speed of the feeding pump of the bioreactor and so, there are always a lower bound u− and an upper bound u+ for u (see [25] for more details). It follows from (4) and Assumption A9 that 0 < s− < s∗ < s+ < si. Consider the following open-loop system Σ3 ds1 dt = χ(t)(si1 −s1) −k1 µ1(s1) x1 (7) dx1 dt = (µ1(s1) −αχ(t)) x1 (8) ds2 dt = χ(t)(si2 −s2) −k2µ1(s1)x1 −k3µ2(s2)x2 (9) dx2 dt = (µ2(s2) −αχ(t))x2, (10) obtained from system (1) after replacing u by an arbitrary bounded measurable control function χ(t) such that χ(t) ∈ [u−,u+] for each t ≥ 0. (11) Consider first the subsystem (7)–(8), which equations do not depend on s2 and x2. Let s11 and s 2 1 be arbitrary real numbers, satisfy- ing the inequalities 0 < s11 < s − 1 < s + 1 < s 2 1 < s i 1. Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 9 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Fig. 1. The parallelogram L(s1,s2) = ABCD, A = (s1,x12), B = (s 1, x̃12), C = (s 2,x22), D = (s 2, x̃22); the parallelogram L(s−,s+) = EFGH, E = (s−,x−2 ), F = (s −, x̃−2 ), G = (s +,x+2 ), H = (s +, x̃+2 ); the lines l1 : s + k3x2 = s i, l2 : s + αk3x2 = s i, l3 : s + k3x2 = si/α. Denote by L(s11,s 2 1) the parallelogram deter- mined by the points (s11,x 1 1), (s 1 1, x̃ 1 1), (s 2 1,x 2 1) and (s21, x̃ 2 1), where the reals x 1 1, x̃ 1 1, x 2 1 and x̃ 2 1 are the solutions of s11 + αk1x 1 1 −s i 1 = 0, s 2 1 + αk1x 2 1 −s i 1 = 0, s11 + k1x̃ 1 1 = s 2 1 + k1x 2 1, s 1 1 + k1x 1 1 = s 2 1 + k1x̃ 2 1. Analogously, a parallelogram L(s−1 ,s + 1 ) can be defined. For reader’s convenience we note that L(s11,s 2 1) and L(s − 1 ,s + 1 ) are similar to the paral- lelograms on Fig. 1 with s1 and x1 instead of s and x2 respectively. The following assertion holds true for the open- loop subsystem (7)–(8). Theorem 3. (cf. [19]) Let the restriction of As- sumptions A1, A7, A8 and A9, concerning s1, x1 and µ1(s1) be fulfilled. Then for each point ζ01 = (s 0 1,x 0 1) ∈ {(s1,x1) : s1 > 0,x1 > 0} and for each measurable function χ(t) ∈ [u−,u+], t ≥ 0, the corresponding solution ϕ1(t,ζ 0 1 ) = (s1(t),x1(t)) of (7)–(8) with s1(0) = s01 and x1(0) = x 0 1 is well defined for all t ∈ [0, +∞) and tends to the parallelogram L(s−1 ,s + 1 ) as t → +∞. Below we shall prove a similar result related to the whole open-loop system (7)–(10). Let us choose arbitrary real numbers s1 and s2, satisfying the inequalities 0 < s1 < s− < s+ < s2 < si. (12) Denote by L(s1,s2) the parallelogram determined by the points (s1,x12), (s 1, x̃12), (s 2,x22) and (s2, x̃22), where the reals x 1 2, x̃ 1 2, x 2 2 and x̃ 2 2 are solutions of the equations s1 + αk3x 1 2 −s i = 0, s2 + αk3x 2 2 −s i = 0, s1 + k3x̃ 1 2 = s 2 + k3x 2 2, s1 + k3x 1 2 = s 2 + k3x̃ 2 2. (13) The parallelograms L(s1,s2) and L(s−,s+) are visualized on Fig. 1. Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 10 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Let the set Ω0 be defined according to (5) and ζ0 = (s01,s 0 2,x 0 1,x 0 2) ∈ Ω0 be an arbitrary point. Denote by ϕ(t,ζ0) = (s̄1(t), s̄2(t), x̄1(t), x̄2(t)), t ≥ 0, the corresponding solution of the open-loop system Σ3 (i. e. of (7)–(10)) with (s̄1(0), s̄2(0), x̄1(0), x̄2(0)) = ζ 0. According to Lemma 1, ϕ(t,ζ0) is bounded and well defined for each t ∈ [0,∞). If we set s̄(t) := k2 k1 s̄1(t) + s̄2(t) and s0 = s(0) = k2 k1 s̄1(0) + s̄2(0), then one can verify that (s̄(t), x̄2(t)), t ≥ 0, satisfy the following open- loop system starting from the point (s0,x02): ds dt = u(si −s) −k3 µ2(s2) x2 dx2 dt = (µ2(s2) −αu) x2, (14) where χ is defined by (11). With s̃ > 0, x̃2 > 0 and δ > 0 we define B(s̃, x̃2; δ) :={(s,x2) : |s−s̃| ≤ δ, |x2−x̃2| ≤ δ}. Proposition 1. Let the Assumptions A1, A7, A8 and A9 be fulfilled, and χ : [0, +∞) → [u−,u+] be the measurable function defined by (11). Then for each point (s̃, x̃2) from the bound- ary of L(s1,s2) there exists δ > 0 such that if (s̄(τ), x̄2(τ)) ∈ B(s̃, x̃2; δ) for some sufficiently large τ ≥ 0, then there exists T > τ so that the point (s̄(T), x̄2(T)) belongs to the interior of the set L(s1,s2) \L(s−,s+). Proof: Let us fix an arbitrary point (s̃, x̃2) from the boundary of L(s1,s2) and let τ0 ≥ 0 be a sufficiently large positive number so that s̄1(τ) ∈ ( s−1 −η,s + 1 + η ) (15) for each τ ≥ τ0, where η := k1 3k2 min { s− −s1,s2 −s+ } > 0. (16) We define z := ( s x2 ) , f(z,u,s2) := ( −k3µ2(s2)x2 + u(sin −s) (µ2(s2) −αu) x2 ) . Let r > 0 be sufficiently small, so that B(s̃, x̃2; r) ∩L(s−,s+) = ∅. We set L = max{‖f ′z(z,u,s2)‖ : z ∈ B(s̃, x̃2; r), u ∈ [u−,u+], s2 ∈ [ms2,M s 2 ]}, where [ms2,M s 2 ] is an interval containing the val- ues of s̄2(t) for t ≥ 0, f ′z(z,u,s2) is the Jacobian of f with respect to z calculated at the point (z,u,s2), and ‖ ·‖ is the Euclidean norm. Let τ be an arbitrary number in (τ0, +∞). Without loss of generality we may find T0 > τ such that T0 − τ > 0 is so small that the solution z̃(·) := (s̃(·), x̃2(·)) of (14), starting from the point z̃ = (s̃, x̃2)T at the moment of time τ is well defined on the interval [τ,T0] and z̃(t) ∈ B(s̃, x̃2; r/2). Let us choose δ0 > 0 to be sufficiently small, so that for each point z ∈ B(s̃, x̃2; δ0) the solution of Σ3, starting from the point z at the moment of time τ is well defined on [τ,T0] and the following inequality holds true 3eL(T0−τ)δ0 < r. (17) Let z be an arbitrary point from B(s̃, x̃2; δ0) and z(t) be the value at the moment of time t of the solution of (14) starting from the point z at the moment of time τ. We set T1 := sup{t ∈ [τ,T0] : z(t) ∈ B(s̃, x̃2; r)}. Clearly T1 > τ. Moreover, for each t ∈ (τ,T1] we have that ‖z(t)−z̃(t)‖≤‖z+ ∫ t τ f(z(ξ),χ(ξ), s̄2(ξ))dξ−z̃ − ∫ t τ f(z̃(ξ),χ(ξ), s̄2(ξ))dξ‖ ≤ ‖z − z̃‖ + ∫ t τ ‖f(z(ξ),χ(ξ), s̄2(ξ))−f(z̃(ξ),χ(ξ), s̄2(ξ))‖dξ ≤‖z − z̃‖ + ∫ t τ L‖z(ξ) − z̃(ξ)‖dξ. Applying the Gronwall inequality we obtain ‖z(t) − z̃(t)‖≤ eL(t−τ)‖z − z̃‖≤ eL(T0−τ)δ0 < r 3 for each t ∈ (τ,T1] (18) Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 11 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes according to the choice of δ0 from (17), and hence ‖z(t)‖ ≤ ‖z̃(t)‖ + ‖z(t) − z̃(t)‖. This means that z(t) ∈ B(s̃, x̃; r/2) + B(0; r/3) = B(s̃, x̃; 5r/6) ⊂ B(s̃, x̃; r). From this inclusion, applied for t = T1, we obtain that T1 = T0. Assume first that (s̃, x̃2) = (s̃(τ), x̃2(τ)) = (s1,x12), i. e. (s̃, x̃2) coincides with the point A in Fig. 1. Without loss of generality, we may think that δ0 ∈ ( 0, 1 2 min{s− −s1,s2 −s+} ) . Let us assume that (s̄(τ), x̄2(τ)) ∈ B(s̃, x̃2; δ0). Using the relations (15) and (16), we obtain that s̄2(τ) = s̄(τ) − k2 k1 s̄1(τ) ≤ s̃ + δ0 − k2 k1 ( s−1 − k1 3k2 (s− −s1) ) = s1 − k2 k1 s−1 + s− −s1 3 + δ0 = s− − (s− −s1) − k2 k1 s−1 + s− −s1 3 + δ0 = s− − k2 k1 s−1 − (s − −s1) + s− −s1 3 + δ0 = s−2 − 2(s− −s1) 3 + δ0 < s−2 − s− −s1 6 < s−2 . Then for each T ∈ ( τ, min { τ + s− −s1 12Ms2 ,T0 }) (19) we have that s̄2(t) = s̄2(τ) + ∫ t τ ˙̄s2(σ)dσ < s−2 − s− −s1 6 + (T − τ)Ms2 < s−2 − s− −s1 12 (20) for each t ∈ (τ,T ] . We set κ := µ2(s − 2 ) −µ2 ( s−2 − s− −s1 12 ) > 0. Then, using (20) and Assumption A9, we obtain that d dt x̃2(t) = x̃2(t)(µ2(s̄2(t)) −αχ(t)) ≤ x̃2(t)(µ2(s̄2(t)) −µ2(s−2 )) ≤−κx̃2(t) < 0, and hence x̃2(t) < e −κ(t−τ)x̃2(τ) = e −κ(t−τ)x12 < x 1 2 (21) for each t ∈ (τ,T ]. The equality (s̃(τ), x̃2(τ)) = (s1,x12) implies s̃(τ) + k3x̃2(τ) − d = 0, where d := si + (1 −α)k3x12. Using the presentation d dt (s̃(t) + k3x̃2(t) −d) = −χ(t)(s̃(t) + k3x̃2(t) −d) + (1 −α)k3χ(t)(x̃2(t) −x12), we obtain that s̃(t) + k3x̃2(t) −d = − ∫ t τ e ∫ ζ t χ(ξ)dξ(1−α)k3χ(ζ)(x12−x̃2(ζ))dζ< 0 (22) for each t ∈ (τ,T ]. As it was shown above, we have that for each t ∈ (τ,T ], z̃(t) ∈ B(s̃, x̃; r/2) and B(s̃, x̃2; r/2)∩ L(s−,s+) = ∅. Moreover, the estimate (20) im- plies that µ2(s̄2(t)) < µ2(s − 2 ) −κ. Taking into account this inequality, (21) and (13), we obtain that ds̃(t) dt = −k3µ2(s̄2(t))x̃2(t) + χ(t)(si − s̃(t)) > −k3(µ2(s−2 )−κ)e −κ(t−τ)x12 +u −(si−s1) > −k3µ2(s−2 )x 1 2 + u −(si −s1) = 0 for each t ∈ (τ,T ]. Then for each T satisfying (19) the inequality s̃(t) > s1 holds true for each t ∈ (τ,T ]. So, we obtain that s̃(t) ∈ (s1,s−) for each t ∈ (τ,T ]. (23) Also, without loss of generality, we may think that L(s1,s2) ∩B(s1,x12; r) = {(s,x2) : s ≥ s1,s+k3x2 ≤ d}∩B(s1,x12; r). Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 12 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Taking into account (22) and (23), we obtain from here that the point z̃(T) = (s̃(T), x̃2(T)) belongs to the interior of the set L(s1,s2) \ L(s−,s+). Hence there exists γ > 0 so that B(s̃(T), x̃2(T); γ) is a subset of the interior of L(s1,s2) \L(s−,s+). Next we choose an arbitrary δ from the interval( 0,γe−L(T−τ) ) and assume that for some τ > τ0 the point (s̄(τ), x̄2(τ)) ∈ B(s̃, x̃2; δ) (note that the real δ > 0 does not depend on τ and on T , but on the difference T−τ). Having fixed δ and τ, we fix an arbitrary T according to (19). For this choice of the real T we apply (18) with δ instead of δ0 and obtain that the point (s̄(T), x̄2(T)) belongs to B(s̃(T), x̃2(T); γ), and hence it belongs to the interior of the set L(s1,s2) \L(s−,s+). The remainder cases concerning the location of the point (s̃, x̃2) on the boundary of the parallel- ogram L(s1,s2) can be considered in analogous way. This completes the proof of Proposition 1. Corollary 1. Let the Assumptions A1, A7, A8 and A9 be fulfilled, and χ : [0, +∞) → [u−,u+] be the measurable function defined by (11). Then for each point (s̃, x̃2) from the boundary of L(s1,s2) there exists δ > 0 such that if (s̄(τ), x̄2(τ) ∈ B(s̃, x̃2; δ) for some sufficiently large τ ≥ 0, then there exists T > τ so that (s̃(t), x̃2(t)), t ∈ (τ,T ], belongs to the interior of the set L(s1,s2) \ L(s−,s+), where (s̃(t), x̃2(t)) denotes the solution at the time moment t of (9)–(10) with initial condition s̃(τ) = s̃ and x̃2(τ) = x̃2. Theorem 4. Let the Assumptions A1, A7, A8 and A9 be fulfilled. Let ζ0 = (s01,s 0 2,x 0 1,x 0 2) ∈ Ω0 be an arbitrary point and ϕ(t,ζ0) = (s̄1(t), s̄2(t), x̄1(t), x̄2(t)), t ≥ 0, be the corre- sponding solution of the open-loop system Σ3 (i. e. of (7)–(10)) with ϕ(0,ζ0) = ζ0, where the measurable function χ is defined according to (11). Then (s̄1(t), x̄1(t)) and (s̄(t), x̄2(t)) tend to the parallelograms L(s−1 ,s + 1 ) and L(s −,s+) respectively, as t → +∞. Proof: Let ζ0 = (s01,s 0 2,x 0 1,x 0 2) ∈ Ω0 be an arbitrarily fixed point and ϕ(t,ζ0) = (s̄1(t), s̄2(t), x̄1(t), x̄2(t)), t ≥ 0, be the corre- sponding solution of the open-loop system Σ3. Theorem 3 implies that (s̄1(t), x̄1(t)) tends to the parallelogram L(s−1 ,s + 1 ) as t → +∞. One can directly check that (s̄(t), x̄2(t)) is a tra- jectory of (9)–(10) starting from the point (s0,x02), where s0 := k2 k1 s01 + s 0 2. We consider the following two cases: Case 1. s(t) + k3x2(t) ≥ si α for each t ≥ 0 (cf. Fig. 1). According to Lemma 1 we have that for each ε > 0 there exists T1 > 0 such that si α ≤ s(t) + k3x2(t) < si α + ε for each t ≥ T1. Hence s(t) + k3x2(t) → si/α as t →∞. Let (s(tk),x2(tk)) → (ŝ, x̂2) for some sequence tk → ∞ as k → ∞. Clearly, ŝ + k3x̂2 = si/α holds true. Without loss of generality we may think that s̄1(tk) → ŝ1 ≥ s−1 and s̄2(tk) → ŝ2 ≥ 0 as tk → +∞. Because s(tk) = k2 k1 s1(tk) + s2(tk), we obtain that ŝ = k2 k1 ŝ1 + ŝ2 ≥ k2 k1 s−1 > 0. Then the relations ŝ > 0 and ŝ + k3x̂2 = si α imply ŝ+αk3x̂2 = si α −(1−α)k3x̂2 > si α −(1−α) si α = si. Hence there exists δ > 0 such that each point (s,x2) satisfying the inequalities |s− ŝ| < δ and |x2 − x̂2| < δ satisfies the inequality s + αk3x2 > si, too. We set M = max { (|u(si−s)−k3µ2(s2)x2|, |x2(µ2(s2)−αu)|) : u∈[u−,u+], |s−ŝ|≤δ, |x2−x̂2|≤δ,s2∈[0,Ms2 ]} Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 13 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes and m= min { s+αk3x2−si : |s−ŝ|≤δ, |x2−x̂2|≤δ } . According to the choice of δ > 0, we have that m > 0. Next we choose the positive reals τ and ε so small that the following inequalities hold true: 2τM < δ and 2ε < τmu−. (24) Denote q(t) := s(t) + k3x2(t)− si α and choose tk so large that |s(tk)−ŝ| < δ 2 , |x2(tk)−x̂2| < δ 2 , 0 < q(tk) < ε. We set τ̄ := sup{σ̃ ∈ [0,τ] : |s(tk + σ) − ŝ| < δ, |x2(tk + σ) − x̂2| < δ for each σ ∈ [0, σ̃]}. Clearly τ̄ > 0. Then, using (24) we obtain that for each σ̃ ∈ [0, τ̄) |s(tk + σ̃) − ŝ| = ∣∣∣∣s(tk) + ∫ tk+σ̃ tk ṡ(ξ)dξ − ŝ ∣∣∣∣ ≤ |s(tk) − ŝ|+∫ tk+σ̃ tk |χ(ξ)(si−s(ξ))−k3µ2(s2(ξ))x2(ξ)|dξ ≤ δ 2 + τM < δ and |x(tk +σ̃)−x̂2| = ∣∣∣∣x2(tk)+ ∫ tk+σ̃ tk ẋ2(ξ)dξ−x̂2 ∣∣∣∣ ≤ |x2(tk) − x̂2| + ∫ tk+σ̃ tk |µ2(s2(ξ))x2(ξ) −αχ(ξ)x2(ξ)|dξ ≤ δ 2 + τM < δ. These inequalities imply that τ̄ = τ, and then for each τ̃ ∈ [0,τ] we have |s(tk + τ̃) − ŝ| < δ and |x2(tk + τ̃) − x̂2| < δ. On the other hand we have that q(tk + τ) = q(tk) + ∫ tk+τ tk q̇(ξ)dξ = s(tk) + k3x2(tk) − si α − ∫ tk+τ tk χ(ξ)(s(ξ) + αk3x2(ξ) −si)dξ < ε−mτu− < −ε. This contradicts the assumption that q(t) = s(t) + k3x2(t)−si/α ≥ 0 for each t ≥ 0 and shows that Case 1 is impossible. Case 2. s(t) + k3x2(t) ≤ si for each t ≥ 0. Similarly to the previous case 1, one can easily show that this case 2 is also impossible. Since the previous two cases 1 and 2 are im- possible, there exists t1 > 0 so that si < s(t1) + k3x2(t1) < si α and 0 < s(t1) < s i. If (s(t1),x2(t1)) ∈ L(s−,s+), we are done. If (s(t1),x2(t1)) 6∈ L(s−,s+) then there exists a parallelogram L(s1,s2) determined by (12) and (13) such that the point (s(t1),x2(t1)) belongs to the boundary of L(s1,s2), and L(s−,s+) is contained in the interior of L(s1,s2). Let L+ be the set of all ω-limit points of the tra- jectory (s̄(t), x̄2(t)) as t →∞, i. e. (s̄, x̄2) ∈ L+ iff there exists a sequence tk tending to infinity as k → ∞ and such that (s̄(tk), x̄2(tk)) tends to (s̄, x̄2) as k → ∞ (cf. [52]). According to Proposition 1, the set L+ is a nonempty compact subset of the parallelogram L(s1,s2). Let us assume the existence of a point (s̄, x̄2) ∈ L+ such that the distance between this point and the set L(s−,s+) is strictly positive. We denote by L(s̄1, s̄2) a parallelogram such that the point (s̄, x̄2) belongs to its boundary and L(s−,s+) is contained in the interior of L(s̄1, s̄2). Now, applying Proposition 1 to the point (s̄, x̄2), the parallelogram L(s̄1, s̄2) and the function χ from (11), we obtain that there exists a neighbor- hood B(s̄, x̄2; δ) of the point (s̄, x̄2) such that if (s̄(τ̄), x̄2(τ̄) ∈ B(s̄, x̄2; δ) for some sufficiently large τ̄ ≥ 0, then there exists T > τ̄ so that the point (s̄(T), x̄2(T)) belongs to the interior of the set L(s̄1, s̄2) \L(s−,s+). Because (s̄, x̄2) is a ω-limit point, there exists a moment of time τ > τ̄ so that (s̄(τ), x̄2(τ)) ∈ B(s̄, x̄2; δ). According to Proposition 1, there ex- ists T > τ such that (s̄(T), x̄2(T)) belongs to the interior of the set L(s̄1, s̄2) \L(s−,s+). Denote by L(ŝ1, ŝ2) a parallelogram containing the point (s̄(T), x̄2(T)) on its boundary, contain- Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 14 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes ing L(s−,s+) in its interior and contained in the interior of L(s̄1, s̄2). According to Corollary 1, the solution (s̄(t), x̄2(t)), t ≥ T , will remain in L(ŝ1, ŝ2). But this is impossible, because (s̄, x̄2) is a ω-limit point. This contradiction completes the proof of Theorem 4. IV. CONCLUSION We investigate a four-dimensional nonlinear dy- namic system, which models anaerobic degrada- tion of soluble organic wastes in a continuous bioreactor with methane production. The main result of the paper (Section III) is the construction of a general and practically oriented approach for stabilizing the model. The aim is to ensure in finite time the values of BOD (denoted by s) to fall onto a given interval [S−,S+], determined by known ecological norms, and to remain there for all time. The approach is based on suitably constructed open-loop control by means of an arbitrary bounded measurable function. The open- loop control approach can serve as a valuable tool for stability investigations using concrete control functions, in particular in the presence of time delays. This new technique can also be considered as a generalization of the previous authors’ results, presented in Section II, which are related to global stabilizability of the model dynamics towards an equilibrium (operating) point using different con- trol approaches. ACKNOWLEDGEMENTS The work of the first author has been accom- plished with the financial support of the Ministry of Education and Science by Grant No DO1- 161/28.08.2018 for the National Geo-Information Centre – part of the Bulgarian National Roadmap for Research Infrastructure. The work of the sec- ond author has been partially supported by the Sofia University “St. Kl. Ohridski” under Grant No 80-10-20/09.04.2019 and by the Bulgarian Science Fund under Grant No KP-06-H22/4/04.12.2018. REFERENCES [1] V. Alcaraz-González, J. Harmand, A. Rapaport, J.- P. Steyer, V. González-Alvarez and C. Pelayo-Ortiz, Robust interval-based siso regulation under maximum uncertainty conditions in an anaerobic digester, In: Intelligent Control, 2001 (ISIC’01), Proc. IEEE In- tern. Symposium, 240–245, IEEE, 2001. [2] V. Alcaraz-González, J. Harmand, A. Rapaport, J.- P. Steyer, V. González-Alvarez and C. Pelayo-Ortiz, Software sensors for highly uncertain WWTPs: a new apprach based on interval observers, Water Research, 36, 2515–2524, 2002. [3] V. Alcaraz-González, J. Harmand, A. Rapaport, J.- P. Steyer, V. González–Alvarez and C. Pelayo-Ortiz, Robust interval-based regulation for anaerobic digestion processes, Water Sci. Technology, 52, 449–456, 2005. [4] V. Alcaraz-González and V. González-Alvarez, Selected topics in dynamics and control of chemical and biolog- ical processes, Springer, Berlin, 2007. [5] J. F. Andrews, A mathematical model for the continuous culture of microorganisms utilizing inhibitory substrates, Biotechnology and Bioengineering, 10, 4, 707–723, 1968. [6] R. Antonelli, J. Harmand, J.-P. Steyer and A. Astolfi, Set-point regulation of an anaerobic digestion process with bounded output feedback, IEEE Trans. Control Systems Tech., 11, 4, 495–504, 2003. [7] L. Appels, J. Lauwers, J. Degréve, L. Helsen, B. Lievens, K. Willems, J. Van Impe and R. Dewil, Anaerobic digestion in global bio-energy production: Potential and research challenges, Renewable and Sus- tainable Energy Reviews, 15, 9, 4295–4301, 2011. [8] G. Bastin and D. Dochain, On-line estimation and adaptive control of bioreactors, Process Measurement and Control, Elsevier, 1990. [9] D. J. Batstone, J. Keller, I. Angelidaki et al, Anaerobic digestion model No. 1 (ADM1), IWA Scientific and Technical Report No. 13, IWA, 2002. [10] L. Benyahia, T. Sari, B. Cherki and J. Harmand, Bifur- cation and stability analysis of a two-step model for monitoring anaerobic digestion processes, Journal of Process Control, 22, 6, 1008–1019, 2012. [11] O. Bernard, Z. Hadj-Sadok, D. Dochain, A. Genovesi and J.-P. Steyer, Dynamical model development and parameter identification for an anaerobic wastewater treatment process, Biotechnology and Bioengineering, 75, 424–438, 2001. [12] D. R. Boone and M. P. Bryant, Propionate-degrading bacterium, syntrophobacter wolinii sp. nov. gen. nov., from methanogenic ecosystems, Applied and Environ- mental Microbiology, 40 (3):626–632, 1980. [13] I. Chang, K. Jang, G. Gil, M. Kim, H. Kim, B. Cho and B. Kim, Continuous determination of biochemical oxygen demand using microbial fuel cell type biosensor, Biosensors and Bioelectronics, 19, 607–813, 2004. [14] N. Dimitrova and M. Krastanov, Nonlinear stabi- lizing control of an uncertain bioprocess model, Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 15 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes Int. Journ. Applied Mathematics and Computer Science, 19, 3, 441–454, 2009. [15] N. S. Dimitrova and M. I. Krastanov, Adaptive asymp- totic stabilization of a bioprocess model with unknown kinetics, Int. J. Numerical Analysis and Modeling, Series B, Computing and Information, 2, 2–3, 200–214, 2011. [16] N. S. Dimitrova and M. I. Krastanov, On the asymp- totic stabilization of an uncertain bioprocess Model, In: Large-Scale Scientific Computing (LSSC’2011), I. Lirkov, S. Margenov, J. Waśniewski (eds.), Lecture Notes in Computer Sciences 7116, Springer, 115–122, 2012. [17] N. Dimitrova and M. Krastanov, Nonlinear adaptive stabilizing control of an anaerobic digestion model with unknown kinetics, Int. J. Robust Nonlinear Control, 22, 15, 1743–1752, 2012. [18] N. S. Dimitrova and M. I. Krastanov, Model-based optimization of biogas production in an anaerobic biodegradation process, Computers and Mathematics with Applications, 68, 9, 986–993, 2014. [19] N. Dimitrova and M. Krastanov, On practice–oriented stabilizability of a chemostat model via bounded open-loop control, to appear in Comptes rendus de l’Academie bulgare des sciences, 2019. [20] J. A. Eastman and J. F. Ferguson, Solubilization of particulate organic carbon during the acid phase of anaerobic digestion, Water Pollution Control Federa- tion, 352–366, 1981. [21] D. Gaida et al. Dynamic real-time substrate feed opti- mization of anaerobic codigestion plants, Leiden Insti- tute of Advanced Computer Science (LIACS), Faculty of Science, Leiden University, 2014. [22] D. Gaida, C. Wolf and M. Bongards, Feed control of anaerobic digestion processes for renewable energy pro- duction: A review, Renewable and Sustainable Energy Reviews 68, 2016, DOI: 10.1016/j.rser.2016.06.096 [23] S. Ghosh, Anaerobic Digestion for renewable energy and environmental restoration, Proceedings of the 8th IWA International Conference on Anaerobic Digestion, 1, 9–15, 1997. [24] S. P. Graef and J. F. Andrews, Mathematical modeling and control of anaerobic digestion, Water Research , 8, 261–289, 1974. [25] F. Grognard and O. Bernard, Stability analysis of a wastewater treatment plant with saturated control, Wa- ter Science Technology, 53, 149–157, 2006. [26] J. Harmand, F. Miens and J. Ph. Steyer, High gain observer for diagnosing an anaerobic fixed bed reactor, In: European Control Conference (ECC) 2001, 2829– 2834, IEEE, 2001. [27] J. Hess and O. Bernard, Design and study of a risk man- agement criterion for an unstable anaerobic wastewater treatment process, Journal of Process Control, 18, 1, 71– 79, 2008. [28] D. T. Hill and C. L. Barth, A dynamic model for simulation of animal waste digestion, Water Pollution Control Federation, 2129–2143, 1977. [29] A. Isidori, Nonlinear control systems, Springer Science & Business Media, 2013. [30] J. Jimenez, E. Latrille, J. Harmand, A. Robles, J. Ferrer, D. Gaida, C. Wolf, F. Mairet, O. Bernard, V. Alcaraz- Gonzalez, et al. Instrumentation and control of anaer- obic digestion processes: a review and some research challenges, Reviews in Environmental Science and Bio/Technology, 14, 4, 615–648, 2015. [31] C. Kravaris, V. Sotiropoulos, C. Georgiou, N. Kazantzis, M. Q. Xiao and A. J. Krener, Nonlinear observer design for state and disturbance estimation, Systems & Control Letters, 56, 11, 730–735, 2007. [32] G. Lyberatos and I. V. Skiadas, Modelling of anaerobic digestion – a review, Global Nest. Int. J., 1, 2, 63–76, 1999. [33] M. Madsen, J. B. Holm-Nielsen and K. H. Esbensen, Monitoring of anaerobic digestion processes: A review perspective, Renewable and Sustainable Energy Re- views, 15, 6, 3141–3155, 2011. [34] L. Maillert, O. Bernard and J.-P. Steyer, Robust regu- lation of anaerobic digestion processes, Water Science and Technology, 48, 6, 87–94, 2003. [35] L. Maillert, O. Bernard and J.-P. Steyer, Nonlinear adaptive control for bioreactors with unknown kinetics, Automatica, 40, 1379–1385, 2004. [36] N. I. Marcos, M. Guay and D. Dochain: Output feedback adaptive extremum seeking control of a continuous stirred tank bioreactor with Monod’s kinetics, J. Process Control, 14, 807–818, 2004 [37] N. I. Marcos, M. Guay, D. Dochain and T. Zhang, Adap- tive extremum-seeking control of a continuous stirred tank bioreactor with Haldane’s Kinetics, J. Process Control, 14, 317–328, 2004 [38] H. O. Méndez-Acosta, R. Femat, V. González-Alvarez and J.-P. Steyer, Substrate regulation in an anaerobic digester based on geometric control, Proceedings IFAC– CLCA, 2002. [39] H. O. Méndez-Acosta, B. Palacios-Ruiz, J. Ph. Steyer, V. Alcaraz-González, E. Latrille and V. González- Álvarez, Nonlinear approach for the vfa regulation in an anaerobic digester, IFAC Proceedings Volumes, 40, 4, 79–84, 2007. [40] J. Mohd Ali, N. H. Hoang, M. A. Hussain and D. Dochain, Review and classification of recent ob- servers applied in chemical process systems, Computers & Chemical Engineering, 76, 27–41, 2015. [41] E. Petre, D. Selisteanu and D. Sendrescu, Adaptive and robust-adaptive control strategies for anaerobic wastewater treatment bioprocesses, Chemical Engineer- ing Journal, 217, 363–378, 2013. [42] P. F. Pind, I. Angelidaki, B. K Ahring, K. Stamatelatou and G. Lyberatos: Monitoring and control of anaero- bic reactors, In: Biomethanation II, 135–182, Springer, 2003. [43] I. Ramirez, E. I. P. Volcke, R. Rajinikanth and J.- P. Steyer, Modeling microbial diversity in anaerobic Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 16 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Neli Dimitrova, Mikhail Krastanov, Model-based control strategies for anaerobic digestion processes digestion through an extended adm1 model, Water Re- search, 43, 11, 2787–2800, 2009. [44] I. Ramirez, A. Mottet, H. Carrére, S. Déléris, F. Ve- drenne and J.-P. Steyer, Modified adm1 disintegra- tion/hydrolysis structures for modeling batch ther- mophilic anaerobic digestion of thermally pretreated waste activated sludge, Water Research, 43, 14, 3479– 3492, 2009. [45] A. Rapaport and J. Harmand, Robust regulation of a class of partially observed nonlinear continuous biore- actors, Journal of Process Control, 12, 2, 291–302, 2002. [46] A. Rapaport, J. Sieber, S. Rodriges and M. Desroches, Extremum seeking via continuation techniques for optimizing biogas production in the chemostat, 9th IFAC Symposium onNonlinear Control Systems (NOL- COS 2013), September 2013, Touluse, France, hal- 00787510v2. [47] C. Rosén and U. Jeppsson, Aspects on adm1 implemen- tation within the bsm2 framework, TEIE, 2005. [48] A. Spagni, M. Ferraris and S. Casu, Modelling wastewater treatment in a submerged anaerobic mem- brane bioreactor, Journal of Environmental Science and Health, Part A, 50, 3, 325–331, 2015. [49] P. Vanrolleghem and D. Lee, On-line monitoring equip- ment for wastewater treatment processes: state of the art, Water Science and Technology, 47, 2, 1–34, 2003. [50] H.-H. Wang, M. Krstic and G. Bastin, Optimizing biore- actors by extremum seeking, Intern. Journal of Adaptive Control and Signal Processing, 13, 651–669, 1999. [51] P. Weiland, Biogas production: current state and per- spectives, Applied Microbiology and Biotechnology, 85, 4, 849–860, 2010. [52] S. Wiggins, Introduction to applied nonlinear dynamical systems and chaos, Springer, New York, 1990. [53] G. S. K. Wolkowicz and H. Xia, Global asymptotic behaviour of a chemostat model with discrete delays, SIAM J. Appl. Math., 57, 1281–1310, 1997. [54] H. Xia, G. S. K. Wolkowicz and L. Wang, Transient oscillations induced by delayed growth response in the chemostat, J. Math. Biol., 50, 489–530, 2005. Biomath 8 (2019), 1907127, http://dx.doi.org/10.11145/j.biomath.2019.07.127 Page 17 of 17 http://dx.doi.org/10.11145/j.biomath.2019.07.127 Introduction Basic properties and global stabilizability of the model dynamics Global stabilizability via input control Global stabilizability via output feedback control Extremum seeking control Open–loop control stabilization Conclusion References