Model-based Control Strategies for Anaerobic Digestion Processes

In this paper we consider a four-dimensional bioreactor model, describing an anaerobic wastewater treatment with methane production. 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 thephase plane.


I. INTRODUCTION
Anaerobic digestion (AD) is a biological process in which organic degradable material is converted into biogas by microorganisms [7], [8].
Recently, AD has been evaluated as one of the most promising processes for waste recovery, environmental protection and bioenergy production.
The biogas is a mixture of gases composed of methane, carbon dioxide, nitrogen, oxygen, hydrogen 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 process 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.) proceeds 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) is called biomass. The number, behavior and interaction of the included organisms pose a challenge to the specialist in the field and have been extensively investigated in the literature. The reactor 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 literature [7], [8].
Batch reactor is a reactor without inflow nor outflow. The digester is filled with the biodegradable 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 without 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 continuously 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 ∈ [u min , u max ] with u min > 0 (cf. [8]).
The AD process follows several phases: Hydrolysis is the step where polymers (macromolecules) 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 bacteria which transform the VFA into acetic acid, hydrogen (H 2 ) and carbon dioxide (CO 2 ). The responsible bacteria for this step produce H 2 and can be inhibited by an excess of H 2 concentration in the digester, that's why they live fixed to the methanogenic bacteria which consume the hydrogen.
Methanogenesis, where methanogenic bacteria reduce the specific substrate into methane.
A good management and control of the AD process 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] modelled 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 researchers, 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 processes (IWA Task Group for Mathematical Mod-elling of Anaerobic Digestion Processes) developed 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 models suitable for process knowledge and simulation, but not appropriate for process control and software 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 monitoring 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 majority 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 measurable 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 significant 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 online 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 controller is designed to satisfy one specific criterion, either economical (maximizing methane production) or ecological (minimizing COD concentration of the effluent) or stability (VFA or dissolved hydrogen) criteria [21], [39]. The controller 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 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 performance 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 developed in recent years to overcome the above drawbacks. Based on a mathematical model, this new approach splits the extremum seeking problem into two steps: global dynamics stabilization and application of a numerical optimization method. In [46] this approach is applied to the classical two-dimensional model of methane fermentation.
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 stabilizability of the model dynamics using different control strategies. Section III reports on a new result, concerning a general and practically oriented stabilization approach by means of an arbitrary measurable bounded control function.

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]): The definition of the state variables s 1 , s 2 and x 1 , x 2 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 s i 1 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 process:  (i) acidogenesis, where the organic substrate s 1 is degraded into volatile fatty acids (VFA) (s 2 ) by acidogenic bacteria (x 1 ); (ii) methanogenesis, where VFA (s 2 ) are degraded by methanogenic bacteria (x 2 ) into methane CH 4 .
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 = Q CH4 is written in the form (2) The model dynamics can be described schematically by the following biological reaction pathways: where r j (·), j = 1, 2, are the reaction rates, which are given by r j (·) = µ j (·)x j . The functions µ 1 (s 1 ) and µ 2 (s 2 ) model the specific 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]) where m 1 , k s1 , m 2 , k s2 and k I are positive coefficients (see Table I).
The most crucial problem in investigating AD models is the formulation of reasonable analytic expressions for µ 1 (s 1 ) and µ 2 (s 2 ). In our theoretical 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.
The model (1) exhibits very rich dynamics. Considering u as a bifurcation parameter, the system possesses different types of bifurcations of the equilibrium points, most of them leading to washout of the biomass and thus to process break-down [10], [14].
The quantity s is called biological oxygen demand (BOD) and represents the biological equivalent 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 The basic properties of the model solutions are summarized in the next lemma, which extends assertions that are given in [25], [53] and [54].
Moreover, for each ε > 0 there exists T ε such that for each t > T ε the following inequalities hold true: 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 subsections 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 control variable. More precisely, we show that for any admissible value of u there exists a nontrivial (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.
Let for someū ∈ (0, u b ) the following assumption holds true: Assumption A2. There exist points s 1 (ū) =s 1 ∈ 0, s i 1 and s 2 (ū) =s 2 ∈ 0, s i 2 , such that the following equalities hold truē 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 pointss 1 ands 2 according to Assumption A2 and compute further Then the point is an equilibrium point of the system (1). In practical applications, the chosen equilibrium point is also called an operating or a reference point.
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 s 2 satisfies the inequality s 2 ≤ s i 2 . The requirement s 2 < s i in Assumption A3(iii) is motivated by the fact that s i can be considered as the worst-case upper bound of the concentration s 2 due to imbalance between acidogenesis and methanogenesis, leading to acidification (x 2 = 0) in the bioreactor (cf. [27]).
Letū ∈ (0, u b ) be chosen according to Assumptions A2 and A3. Denote by Σ 1 the system obtained from (1) by substituting the control variable u byū. Then the following Theorem 1 reports on the global stability of the equilibrium pointp. 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 control
This subsection is devoted to global stabilizability 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 ≡ κ(s 2 , x 2 ) = βk 4 µ 2 (s 2 )x 2 , 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 ζ = (s 1 , x 1 , s 2 , x 2 ). Using the definition of s i from (4) we make the following assumption: Assumption A4. Lower bounds s i− and k − 4 for the values of s i and k 4 respectively, as well as an upper bound k + 3 for the value of k 3 are known. Define the following feedback control law: 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).
, +∞ and let obviously,ξ belongs to the interval (0, s i ). The next assumption is similar to the regulability Assumption A2.
Assumption A5. There exists a points 1 such that Finds 1 according to Assumption A5 and definē It is straightforward to see that the point is an equilibrium point of Σ 2 . We shall show below that the feedback law (6) stabilizes asymptotically the closed-loop system towardsp β (cf. [16]).
Let Ω 0 be defined according to (5). Using the definitions of s and s i from (4), we define the sets Assumption A6. Let the inequality be satisfied for each s 1 , belonging to the projection of the set Ω ∩ Ω 2 on the s 1 -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. , +∞ and let p β = (s 1 ,x 1 ,s 2 ,x 2 ) be the corresponding equilibrium point. Then the feedback control law κ(·) defined by (6) stabilizes asymptotically the control system Σ 2 to the pointp β for each starting point ζ 0 = (s 0 1 , 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 modelbased 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 Q max 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 according 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 q max ∈ (q − , q + ) where Q(q) takes a maximum, Q max = Q(q max ), the function strictly increases in the interval (q − , q max ) and strictly decreases in (q max , q + ).
Denote by E(q) the equilibrium point parameterized on q and let E(q max ) be the steady state where Q max is achieved. The goal is to stabilize in real time the systems Σ 1 and Σ 2 towards this (a priori unknown) equilibrium point E(q max ) and therefore to the maximum methane flow rate Q max . 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 q 1 , q 2 , . . . , q n , . . . from the interval (q − , q + ) is constructed, each q j being in the form q j = q j−1 ± h j , h j > 0, and such that {q j } tends to q max ; Theorems 1 and 2 guarantee that the dynamics is globally asymptotically stabilizable towards the equilibrium E(q j ), j = 1, 2, . . .. Then by computing and comparing the values Q(q 1 ), Q(q 2 ), . . . , Q(q n ), . . ., the desired equilibrium point E(q max ) and thus Q max 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 q max and Q max respectively; in the second stage, the interval [q] is refined using an elimination procedure based on the golden mean value (or Fibonacci search) strategy. The second stage produces the final intervals 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 previously chosen equilibrium (operating) point. This section proposes a different approach for stabilizing 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 (s j ) are strictly increasing in the intervals (0, s i j ), j = 1, 2. Assumption A7 is technical. It is always fulfilled when the functions µ j (s j ), j = 1, 2, are presented by the Monod specific growth rate.
Let s and s i be defined according to (4), i. e.
Let us choose an operating (reference) point s * , Assumption A8 (regulability). There exists a point s * 1 such that Find s * 1 according to Assumption A8 and determine further It is straightforward to see that the point is an equilibrium point of system (1). One can directly check that the equilibrium points of (1) satisfy the equalities (cf. the straight line l 2 in Fig. 1) Assumption A9. The positive reals s − 1 , s + 1 , s − 2 , s + 2 , u − and u + satisfy the relations 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).
Consider the following open-loop system Σ 3 obtained from system (1) after replacing u by an arbitrary bounded measurable control function χ(t) such that Consider first the subsystem (7)- (8), which equations do not depend on s 2 and x 2 .
Analogously, a parallelogram L(s − 1 , s + 1 ) can be defined. For reader's convenience we note that L(s 1 1 , s 2 1 ) and L(s − 1 , s + 1 ) are similar to the parallelograms on Fig. 1 with s 1 and x 1 instead of s and x 2 respectively.
Below we shall prove a similar result related to the whole open-loop system (7)-(10).
Proof: Let us fix an arbitrary point (s,x 2 ) from the boundary of L(s 1 , s 2 ) and let τ 0 ≥ 0 be a sufficiently large positive number so that for each τ ≥ τ 0 , where We define Let r > 0 be sufficiently small, so that B(s,x 2 ; r) ∩ L(s − , s + ) = ∅. We set is an interval containing the values ofs 2 (t) for t ≥ 0, f z (z, u, s 2 ) is the Jacobian of f with respect to z calculated at the point (z, u, s 2 ), and · is the Euclidean norm.
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 [τ, T 0 ] and the following inequality holds true 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 Clearly T 1 > τ . Moreover, for each t ∈ (τ, T 1 ] we have that Applying the Gronwall inequality we obtain 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 = T 1 , we obtain that T 1 = T 0 . Assume first that (s,x 2 ) = (s(τ ),x 2 (τ )) = (s 1 , x 1 2 ), i. e. (s,x 2 ) coincides with the point A in Fig. 1.
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(s 1 , s 2 ) \ L(s − , s + ).
The remainder cases concerning the location of the point (s,x 2 ) on the boundary of the parallelogram L(s 1 , s 2 ) can be considered in analogous way. This completes the proof of Proposition 1.
On the other hand we have that q(t k + τ ) = q(t k ) + This contradicts the assumption that q(t) = s(t) + k 3 x 2 (t) − s i /α ≥ 0 for each t ≥ 0 and shows that Case 1 is impossible.
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 impossible, there exists t 1 > 0 so that s i < s(t 1 ) + k 3 x 2 (t 1 ) < s i α and 0 < s(t 1 ) < s i .
Let L + be the set of all ω-limit points of the trajectory (s(t),x 2 (t)) as t → ∞, i. e. (s,x 2 ) ∈ L + iff there exists a sequence t k tending to infinity as k → ∞ and such that (s(t k ),x 2 (t k )) 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(s 1 , s 2 ).
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 neighborhood 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 + ).
Denote by L(ŝ 1 ,ŝ 2 ) a parallelogram containing the point (s(T ),x 2 (T )) on its boundary, contain-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 dynamic system, which models anaerobic degradation 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 openloop 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 control approaches.