Next Article in Journal
Towards Productive and Ergonomic Order Picking: Multi-Objective Modeling Approach
Previous Article in Journal
Decontamination of Food Packages from SARS-CoV-2 RNA with a Cold Plasma-Assisted System
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Global Stability Analysis of the Model of Series/Parallel Connected CSTRs with Flow Exchange Subject to Persistent Perturbation on the Input Concentration

by
Alejandro Rincón
1,2,
Fredy E. Hoyos
3,* and
John E. Candelo-Becerra
4
1
Grupo de Investigación en Desarrollos Tecnológicos y Ambientales—GIDTA, Facultad de Ingeniería y Arquitectura, Universidad Católica de Manizales, Carrera 23 N. 60-63, Manizales 170002, Colombia
2
Grupo de Investigación en Microbiología y Biotecnología Agroindustrial—GIMIBAG, Instituto de Investigación en Microbiología y Biotecnología Agroindustrial, Universidad Católica de Manizales, Carrera 23 N. 60-63, Manizales 170002, Colombia
3
Facultad de Ciencias, Escuela de Física, Universidad Nacional de Colombia–Sede Medellín, Carrera 65 No. 59A, 110, Medellín 050034, Colombia
4
Departamento de Energía Eléctrica y Automática, Facultad de Minas, Campus Robledo, Universidad Nacional de Colombia–Sede Medellín, Carrera 80 No. 65–223, Medellín 050041, Colombia
*
Author to whom correspondence should be addressed.
Appl. Sci. 2021, 11(9), 4178; https://doi.org/10.3390/app11094178
Submission received: 16 March 2021 / Revised: 22 April 2021 / Accepted: 28 April 2021 / Published: 4 May 2021
(This article belongs to the Section Electrical, Electronics and Communications Engineering)

Abstract

:
In this paper, we study the convergence properties of a network model comprising three continuously stirred tank reactors (CSTRs) with the following features: (i) the first and second CSTRs are connected in series, whereas the second and third CSTRs are connected in parallel with flow exchange; (ii) the pollutant concentration in the inflow to the first CSTR is time varying but bounded; (iii) the states converge to a compact set instead of an equilibrium point, due to the time varying inflow concentration. The practical applicability of the arrangement of CSTRs is to provide a simpler model of pollution removal from wastewater treatment via constructed wetlands, generating a satisfactory description of experimental pollution values with a satisfactory transport dead time. We determine the bounds of the convergence regions, considering these features, and also: (i) we prove the asymptotic convergence of the states; (ii) we determine the effect of the presence of the side tank (third tank) on the transient value of all the system states, and we prove that it has no effect on the convergence regions; (iii) we determine the invariance of the convergence regions. The stability analysis is based on dead zone Lyapunov functions, and comprises: (i) definition of the dead zone quadratic form for each state, and determination of its properties; (ii) determination of the time derivatives of the quadratic forms and its properties. Finally, we illustrate the results obtained by simulation, showing the asymptotic convergence to the compact set.

1. Introduction

Several systems converge to a compact set instead of an equilibrium point, for instance: (i) chaotic systems [1,2,3]; (ii) closed-loop systems subject to external disturbances or model uncertainties [4,5,6,7], or input saturation [8,9,10]. For this type of system, the global stability analysis considers large but bounded initial values of the state variables, and comprises determining the convergence region of the state variables, proving asymptotic or exponential global stability, and proving the invariance of the convergence region [2,3]. To this end, the Lyapunov function can be used, which consists of a radially unbounded function, possibly a quadratic form.
There are three common approaches based on the Lyapunov function for studying the aforementioned systems converging to compact sets. In the first approach, the Lyapunov function is positive definite, and can include quadratic forms of states with some deviation. The expression of the time derivative of the Lyapunov function is arranged as a first order differential equation, and it is concluded that the Lyapunov function converges exponentially to a compact set [2,3,11]. As an example, in Zhang (2017), the stability of a Lorenz system is determined for three different ranges of a model parameter [3]. This system consists of three ordinary differential equations (ODEs). In Liao (2017), the stability of a Yang-Chen system is analyzed [2]. The system consists of three ODE’s with states x, y, z. Global stability of both compact sets and equilibrium points are studied. To study global exponential convergence to compact sets, different radially unbounded Lyapunov functions are constructed for different quadrant regions of the state plane. Also, the ranges of a model parameter that imply global stability or instability of the equilibrium points are determined.
In the second approach, the Lyapunov function is positive definite and the expression of the time derivative of the Lyapunov function is negative when the state variables are outside the convergence region. It is concluded that the system converges asymptotically [3,11].
The third approach comprises dead-zone radially unbounded positive functions. This approach has been traditionally used for robust control design, but not for open-loop systems. An early development is presented in [12,13,14], and further developments in [15,16,17]. The main advantages of this approach are: (i) it allows rigorous proof of the asymptotic convergence to compact sets, through Barbalat’s lemma; (ii) it facilitates stability analysis for the case that asymptotic but not exponential convergence occurs [16,17,18]. In Rincon (2020), this approach was used to determine the convergence of a continuously stirred tank reactor (CSTR) open-loop system that comprises three reactions in the presence of an external disturbance [19]. This system is described by three cascaded differential equations whose states converge to a compact set. In the particular case of bioreaction systems, there are several global stability studies, proving convergence to equilibrium points, but not the convergence to compact set [20,21,22].
Another result of global stability analysis for systems converging to compact sets is BIBO stability. It can be performed for the case of bounded or L2 disturbances, and it results in an input–output relationship that relates the system output or states in terms of either the input signal or the time integral of the squared input signal, where the input signal may be an external disturbance or a manipulated input. This expression is obtained by integrating the time derivative of the Lyapunov function. A BIBO stability study for an open loop system is presented in [23], where a non-linear Caputo fractional system with time varying bounded delay is considered. A quadratic Lyapunov function is defined, and the resulting fractional derivative is a function of the negative Lyapunov function. As a consequence of this expression, the upper bound of the Lyapunov function, system states and system output are a function of the upper bound of the input.
BIBO stability is also determined for closed loop systems (see [24,25,26,27,28]). In [24], a semi-active controller is formulated for a vibrating structure. The Lyapunov analysis is combined with passivity theory and uses a storage function as a Lyapunov function. The BIBO stability result relates the displacements and velocities as a function of L2 external forces. Also, the ranges of the controller parameters that lead to BIBO stability are defined. In [25], a robust observer-based controller design is formulated for a class of nonlinear systems described by Tagaki–Sugeno (TS) models, subject to persistent bounded external disturbance. A non-quadratic Lyapunov function is used. The BIBO stability result relates the closed loop states as a function of the external disturbance. In [26], an adaptive neural controller is developed for a n-link rigid robotic manipulator. Neural networks are used for approximating the unknown dynamics, and the backstepping strategy is used as a basic control framework. As a result of the controller design, the overall Lyapunov function (V2) comprises the subsystem Lyapunov functions associated to the zi states and the parameter updating error. The dV2/dt expression is a function of V2, so that it is concluded that the zi states converge to compact sets whose width depends on the neural approximation error.
In this paper, we study the stability of a network model comprising three CSTRs with the following features: (i) the pollutant concentration in the inflow to the first CSTR is time varying; (ii) the first and second CSTRs are connected in series, whereas the second and third CSTRs are connected in parallel with flow exchange; (iii) the states converge to a compact set instead of an equilibrium point. The practical applicability of the arrangement of CSTRs is to provide a simpler model of pollution removal from wastewater treatment via constructed wetlands, generating a satisfactory description of experimental pollution values with a satisfactory transport dead time. The aforementioned features make the global stability analysis complex. The convergence sets of the system states are determined as a function of the concentration of the inflow to the first tank, and the global asymptotic convergence of the states to these sets and their invariant nature are proved. Also, the effect of the presence of flow exchange on the transient value of the system states is determined, and it is concluded that it has no effect on the convergence regions. The study is based on global stability using dead zone Lyapunov functions and Barbalat’s lemma.
The contribution of this work with respect to closely related studies on global stability of continuous stirred tank reactors (CSTR) are:
  • we consider three connected CSTRs, including a side CSTR with flow exchange, whereas related studies consider a single CSTR with several states (e.g., biomass and substrate concentrations) [20,21,29,30];
  • we consider the case that the system converges to a compact set, whereas related studies consider convergence to an equilibrium point [20,21,29,30].
The paper is organized as follows. The preliminaries and description of the model are presented in Section 2; the stability analysis is presented in Section 3, including the first tank (Section 3.1) and both the second and third tanks (Section 3.2). The numerical simulation is presented in Section 4, and the conclusions in Section 5.

2. Model Description

Modelling of constructed wetlands is complex because of the large number of processes involved in pollution removal, the hydraulics in a porous medium with the diffusive flow, and the dead time. In addition, the pollution removal is influenced considerably by the hydraulics and environmental conditions. The last generation of models (process—based models) are rigorous but overly complex, with difficult estimation of its parameters [31,32,33]. In contrast, series/parallel connection of CSTRs is simple and capable of describing transport dead time, diffusion and reaction [31].
In Marsili (2005), a simple model is proposed for describing the disperse hydraulics and pollution removal dynamics in horizontal subsurface constructed wetlands [31]. The structure of the hydraulic model comprises three series CSTRs, followed by two parallel CSTRs, and a plug flow reactor. The parallel CSTRs involve no flow exchange. The kinetics is of the Monod type, and only depends on the concentration of the modelled pollutant. The outlet and inlet flows of the system are the same, and the liquid CSTR volumes are constant. The model was fitted to both experimental tracer data and current (no tracer) data.
In Davies (2007) and Freire (2009), a simple modelling strategy is proposed for describing hydraulics and pollution removal dynamics corresponding to tracer data from vertical subsurface constructed wetlands [34,35]. The model comprises three CSTRs, the first and second CSTRs are connected in series, whereas the second and third CSTRs are connected in parallel with flow exchange. The inlet flow is considered different from the outlet flow, and no model is considered for it. The liquid volumes of the parallel CSTRs are constant, whereas the liquid volume of the first CSTR is considered varying and is described by a first-order differential equation that is fitted to experimental data. The specific reaction rate is of the first order and it only depends on the concentration of the modelled pollutant.
We use the global stability definition as follows.
Definition 1.
Consider the system dX/dt = f(X), where X = (x1, x2, …, xn) ∈ n, f: nn, X(t,to,Xo) denoted as X(t) is a solution to the system, and X0n is an initial value of X(t). Consider the compact set ΩQn. Define the distance between the solution X and the compact set ΩQ as:
ρ ( X , Ω Q ) = i n f Y Ω Q | | X Y | | .
If for all X0 satisfying X0n, X0 ∉ ΩQ the property lim t ρ ( X , Ω Q ) = 0 holds, then the compact set ΩQ is called a global convergence set of the system, and the system is called globally asymptotically stable.
We consider a series/parallel CSTRs-based model for representing the hydraulics and pollution removal dynamics in subsurface constructed wetlands, resulting from a combination of the models proposed by [31,34]. A single pollutant and a time-varying inflow pollutant concentration are considered. The CSTRs connection comprises three CSTRs, the first and second CSTRs are connected in series, whereas the second and third CSTRs are connected in parallel with flow exchange (Figure 1). The mass balances for the CSTRs are:
d S 1 d t = ( S i n S 1 ) Q v 1 r s 1 ,
d S 2 d t = ( S 1 S 2 ) Q v 2 r s 2 + Q D v 2 ( S 3 S 2 ) ,
d S 3 d t = ( S 2 S 3 ) Q D v 3 ,
S i n = S ¯ i n + δ s i n ,
where S 1 is the pollutant concentration in CSTR 1; S 2 is the pollutant concentration in CSTR 2; S 3 is the pollutant concentration in CSTR 3; S i n is the inlet pollutant concentration; Q is the wastewater flow entering CSTR 1 and CSTR 2; Q D is the wastewater flow entering and leaving CSTR 3; v i is the liquid volume of the ith CSTR; r s 1 is the reaction rate of CSTR 1; r s 2 is the reaction rate of CSTR 2. In addition, S 1 is defined in + ,   S 2 is defined in + , S 3 is defined in + .
Assumptions: 
Assumption 1.
Q , Q D , v 1 , v 2 , and v 3 are positive and constant. This implies that Q / v 1 , Q / v 2 , Q D / v 2 and Q D / v 3 are constant and positive.
Assumption 2.
S ¯ i n is positive and constant, δ s i n is varying but bounded, m a x { δ s i n } > 0 ,   m i n { δ s i n } < 0 .
Assumption 3.
The reaction rate rs1 only depends on S 1 , it is bounded for S 1 bounded, d r s 1 d S 1 0 ; r s 2 only depends on S 2 , it is bounded for S 2 bounded, and d r s 2 d S 2 0 .
In the case of constant S i n , δ s i n = 0 in Equation (4), S 1 converges to S 1 e q , S 2 converges to S 2 e q and S 3 converges to S 3 e q , where S 1 e q , S 2 e q and S 3 e q are positive constants obtained by equating Equations (1)–(3) to zero:
0 = ( S ¯ i n S 1 ) Q v 1 r s 1 ,
0 = ( S 1 S 2 ) Q v 2 r s 2 + Q D v 2 ( S 3 S 2 ) ,
0 = ( S 2 S 3 ) Q D v 3 .
Let
x 1 = S 1 S 1 e q ,
x 2 = S 2 S 2 e q ,
x 3 = S 3 S 3 e q .
The x 1 dynamics is obtained by subtracting Equation (5) from Equation (1):
d x 1 d t = D ( x 1 + r x 1 D δ s i n ) ,
r x 1 = r s 1 r s 1 e q ,
where D = Q / v 1 > 0 , r s 1 e q is r s 1 evaluated at S 1 = S 1 e q .
The x 2 dynamics is obtained by subtracting Equation (6) from Equation (2):
d x 2 d t = ( x 1 x 2 ) Q v 2 r x 2 Q D v 2 ( x 2 x 3 ) ,
r x 2 = r s 2 r s 2 e q .
The term r s 2 e q is r s 2 evaluated at S 2 = S 2 e q .
The x 3 dynamics is obtained by subtracting Equation (7) from Equation (3):
d x 3 d t = ( x 2 x 3 ) Q D v 3 .
In addition, x 1 is defined in [− S 1 e q ,   ) ,   x 2 is defined in [− S 2 e q ,   ) , and x 3 is defined in [− S 3 e q ,   ) . The model is fitted to current (no tracer) data in Section 3.

3. Global Stability Analysis

In this section, the stability analysis is performed for the network model (1), including the determination of the convergence set, the asymptotic convergence and the invariance. The stability analysis is based on Lyapunov theory and Barbalat’s lemma. Recall that early global stability studies based on dead zone Lyapunov functions are presented in [12,13,14], whereas later studies are presented in [16,17,19].
A dead zone quadratic form V 1 is defined for the first state x 1 , V 2 for the second state x 2 and V 3 for the third state x 3 , and their properties are determined. An overall Lyapunov function V is defined as a weighted sum of the quadratic forms V 1 , V 2 , and V 3 . The time derivative of each dead zone quadratic form is determined. The stability analysis of the first tank is independent of the other tanks, whereas the stability analysis of the second and third tanks is simultaneous and requires the stability results of the first tank.

3.1. Stability Analysis for the First Tank

The stability analysis is based on determining the non-positive nature of d V 1 / d t , where V 1 is the dead-zone quadratic form for x 1 . The function V 1 and its gradient d V 1 / d x 1 are defined so as to achieve this goal.
The procedure comprises the following tasks: (i) determination of the expression for d V 1 / d t ; (ii) definition of the gradient d V 1 / d x 1 and the dead zone quadratic form V 1 for the state x 1 ; (iii) arrangement of the d V 1 / d t expression in terms of a non-positive function of d V 1 / d x 1 ; (iv) integration of the d V 1 / d t expression. The definition of the gradient d V 1 / d x 1 . and V 1 requires the determination of the d V 1 / d t expression and the properties of the terms involved.
Theorem 1.
Consider the model in Equations (11) and (12), subject to assumptions 1 to 3. (Ti) The state x 1 converges asymptotically to Ω x 1 , Ω x 1 = { x 1 : x 1 l x 1 x 1 u } , where x 1 u = { x 1 : g x 1 m a x { δ s i n } = 0   } , x 1 u > 0 , x 1 l = { x 1 : g x 1 m i n { δ s i n } = 0   } ,   x 1 l < 0 , and g x 1 = x 1 + r x 1 D , and S 1 converges asymptotically to Ω S 1 where Ω S 1 = [ S 1 l     S 1 u ] , S 1 l = x 1 l + S 1 e q , S 1 u = x 1 u + S 1 e q and S 1 e q is provided by Equation (5). (Tii) The sets Ω x 1   and Ω S 1 are invariant.
Proof. 
Task 1. Determination of the d V 1 / d t expression.
The time derivative of V 1 can be expressed as:
d V 1 d t = f v 1 d x 1 d t ,
f v 1 = d v 1 d x 1 .
Incorporating the x 1 dynamics (11) and arranging, yields:
d V 1 d t = ( 1 ) D f v 1 ( g x 1 δ s i n ) ,
where,
g x 1 = x 1 + r x 1 D ,
d g x 1 d x 1 = 1 + 1 D d r x 1 d x 1 1 .
Task 2.
Definition of the gradient d V 1 / d x 1 and the dead zone quadratic form V 1 for the state x 1 .
At what follows, we examine the properties of the term g x 1 δ s i n appearing in Equation (18). From Equations (12), (19) and (20) and Assumption 3 we have:
d g x 1 d x 1 1 ;   g x 1 = 0   f o r   x 1 = 0 ;   s g n ( g x 1 ) = s g n ( x 1 ) ,
g x 1 δ s i n > 0         f o r     x 1 > x 1 u ,
g x 1 δ s i n < 0         f o r     x 1 < x 1 l ,
where,
x 1 u = { x 1 : g x 1 m a x { δ s i n } = 0   } , x 1 u > 0 ,
x 1 l = { x 1 : g x 1 m i n { δ s i n } = 0   } ,   x 1 l < 0 .
Therefore,
s g n ( g x 1 δ s i n ) = s g n ( x 1 ) = s g n ( g x 1 )   f o r   x 1 [ x 1 l , x 1 u ] .
To obtain non-positive d V 1 / d t in Equation (18), we need to choose f v 1 such that:
f v 1 ( g x 1 δ s i n ) < 0       f o r     x 1   [ x 1 l , x 1 u ] ,
f v 1 ( g x 1 δ s i n ) = 0     f o r     x 1     [ x 1 l , x 1 u ] .
In view of properties (26) and to fulfill conditions (27) and (28), we choose:
f v 1 = { x 1 x 1 u   f o r   x 1 > x 1 u 0   f o r     x 1         [ x 1 l , x 1 u ] . x 1 x 1 l   f o r   x 1 < x 1 l
The main properties of f v 1 are:
( P i ) f v 1   i s   c o n t i n u o u s   w i t h   r e s p e c t   t o   x 1   ,
( P i i )   f v 1 = 0   f o r   x 1   [ x 1 l , x 1 u ] ,
( P i i i )   f v 1 0   f o r   x 1     [ x 1 l , x 1 u ] ,
P ( i v )   s g n ( f v 1 ) = s g n ( x 1 ) 0   f o r   x 1     [ x 1 l , x 1 u ] ,
( P v )   s g n ( f v 1 ) = s g n ( g x 1 ) = s g n ( g x 1 δ s i n ) 0   f o r   x 1     [ x 1 l , x 1 u ] .
A dead-zone Lyapunov function that satisfies d V 1 / d x 1 = f v 1 (Equation (17)) is:
V 1 = 1 2 f v 1 2 .
The main properties of V 1 are:
V 1 > 0   f o r       x 1     [ x 1 l , x 1 u ] ,
V 1 = 0   f o r   x 1   [ x 1 l , x 1 u ] .
Task 3.
Arrangement of the d V 1 / d t expression in terms of a non-positive function of f v 1 = d V 1 / d x 1 .
From properties (31), (34), it follows that the term f v 1 ( g x 1 δ s i n ) appearing in Equation (18) satisfies:
f v 1 ( g x 1 δ s i n ) = | f v 1 | | ( g x 1 δ s i n ) | 0 .
From the definition of g x 1 (19) and definitions of x 1 l ,   x 1 u (24), (25), it follows that:
g x 1 δ s i n | g t 1 | ,
g t 1 = { g x 1 m a x { δ s i n }     f o r     x 1 > x 1 u 0     f o r       x 1     [ x 1 l ,     x 1 u ] g x 1 m i n { δ s i n }     f o r     x 1 < x 1 l .
From Equations (38) and (39), it follows that:
f v 1 ( g x 1 δ s i n ) | f v 1 | | g t 1 | = f v 1 g t 1 0 .
From the property (21), definition of g t 1 (40), and definition of f v 1 (29), it follows that:
| g t 1 | | f v 1 | .
Therefore, Equation (41) yields:
f v 1 ( g x 1 δ s i n ) f v 1 2 0 .
Using this in Equation (18) yields:
d V 1 d t D f v 1 2 .
Accounting for the definition of V 1 (35), we have:
d V 1 d t 2 D V 1 0 .
Task 4.
Integration of the d V 1 / d t expression.
From the above expression, it follows that:
V 1 V 1 t o e 2 D ( t t o ) .
That is, V 1 (35) converges exponentially to zero, and consequently f v 1 converges exponentially to zero. Further, accounting for the definition of f v 1 (29), it follows that x 1 converges asymptotically to Ω x 1 = { x 1 : x 1 l x 1 x 1 u } . From this convergence result and the definition of x 1 (8), it follows that S 1 converges asymptotically to Ω S 1 ; where,
Ω S 1 = [ S 1 l S 1 u ] ,
S 1 l = x 1 l + S 1 e q ,   S 1 u = x 1 u + S 1 e q   .
S 1 e q is provided by Equation (5), the bounds x 1 l , x 1 u are defined in Equations (24) and (25). This completes the proof of Ti.
From Equation (44) and properties (31) and (32), it follows that d V 1 d t < 0   f o r   x 1     Ω x 1 and d V 1 d t 0   f o r   x 1   Ω x 1 . Therefore, the set Ω x 1 is invariant, and consequently, Ω s 1 is invariant. This completes the proof of Tii. □
Remark 1.
The property | g t | | f v 1 | (42) is crucial for proving the asymptotic convergence of x1.
Remark 2.
The invariance of the set Ω s 1 stated in Theorem 1 implies that once the state S 1 is inside the convergence region Ω s 1 , that is, S 1 Ω s 1 , it remains inside afterwards.

3.2. Stability Analysis for the Second and Third Tanks

The stability analysis and majorly the proof for asymptotic convergence of x 2 and x 3 to their compact sets, is based on determining the non-positive nature of d V / d t = d V 1 / d t + k 2 d V 2 / d t + k 3 d V 3 / d t , where V = V 1 + k V 2 + k V 3 , V is the overall Lyapunov function, V 2 is the dead zone quadratic form for the state x 2 and V 3 the dead zone quadratic form for the state x 3 , and k 2 and k 3 are positive constants. To achieve this goal, the functions V 2 and V 3 and the gradients d V 2 / d x 2 and d V 3 / d x 3 are defined accordingly, and the non-positive nature of the terms of the d V / d t expression is determined.
The procedure comprises the following tasks: (i) determination of the expression for d V 2 / d t + k d V 3 / d t ; (ii) arrangement of the expression for d V 2 / d t + k d V 3 / d t in terms of f v 1 ; (iii) definition of gradient d V 2 / d x 2 and dead zone quadratic form V 2 ; (iv) arrangement of the expression for d V 2 / d t + k d V 3 / d t in terms of a non-positive function of d V 2 / d x 2 ; (v) definition of gradient d V 3 / d x 3 and dead zone quadratic form V 3 ; (vi) definition of the overall Lyapunov function V and determination of the non-positive nature of the expression for d V / d t = d V 1 / d t + k d V 2 / d t + k d V 3 / d t ; (vii) integration of the d V / d t expression. The definition of the gradients d V 2 / d x 2 and d V 3 / d x 3 and the quadratic forms V 2 and V 3 requires the determination of the expression for d V 2 / d t + k d V 3 / d t and the properties of the terms involved.
Theorem 2.
Consider the model (11)–(15), subject to assumptions 1 to 3. (Ti) the state x 2 converges asymptotically to Ω x 2 , Ω x 2 = { x 2 : x 2 l x 2 x 2 u } , where:
x 2 u = { x 2 :   g x 2 m a x { δ 2 } = 0 } ,   x 2 u > 0 x 2 l = { x 2 :   g x 2 m i n { δ 2 } = 0 } ,   x 2 l < 0 g x 2 = x 2 + r x 2 Q / v 2 , δ 2 = d x 1 , d x 1 = { x 1 u   f o r   x 1 > x 1 u x 1   f o r   x 1     [ x 1 l ,   x 1 u ] x 1 l   f o r   x 1 < x 1 l .
S 2 converges asymptotically to Ω S 2 where Ω S 2 = [ S 2 l     S 2 u ] , S 2 l = x 2 l + S 2 e q , S 2 u = x 2 u + S 2 e q   and S 2 e q is provided by Equations (5)–(7). (Tii) the state x3 converges asymptotically to Ω x 3 , Ω x 3 = { x 3 : x 3 l x 3 x 3 u } , where x 3 u = x 2 u ;   x 3 l = x 2 l , and S 3 converges asymptotically to Ω S 3 where Ω S 3 = [ S 3 l     S 3 u ] , S 3 l = x 3 l + S 3 e q , S 3 u = x 3 u + S 3 e q   and S 3 e q is provided by Equations (5)–(7). (Tiii) The sets Ω x = Ω x 1   Ω x 2   Ω x 3 and Ω S = Ω S 1   Ω S 2   Ω S 3 are invariant.
Proof. 
Task 1. Determination of the expression for d V 2 / d t + k d V 3 / d t .
The d x 2 / d t expression (13) can be rewritten as:
d x 2 d t = Q D v 2 ( x 2 x 3 ) + ( 1 ) Q v 2 ( g x 2 x 1 ) ,
where:
g x 2 = x 2 + r x 2 Q / v 2 .
The time derivative of V 2 can be expressed as:
d V 2 d t = f v 2 d x 2 d t ,
f v 2 = d V 2 d x 2 .
Substituting the d x 2 / d t expression (Equation (48)) and arranging, yields:
d V 2 d t = f v 2 ( 1 ) Q D v 2 ( x 2 x 3 ) + ( 1 ) Q v 2 f v 2 ( g x 2 x 1 ) .
The time derivative of V 3 can be expressed as:
d V 3 d t = f v 3 d x 3 d t ,
f v 3 = d V 3 d x 3 .
Substituting the d x 3 / d t expression (Equation (15)) and arranging, yields
d V 3 d t = f v 3 ( ( x 2 x 3 ) Q D v 3 ) .
Hence,
d d t ( v 3 v 2 V 3 ) = f v 3 ( ( x 2 x 3 ) Q D v 2 ) .
Adding d V 2 / d t (52) and d V 3 / d t (56), yields
d V 2 d t + d d t ( v 3 v 2 V 3 ) = ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) + ( 1 ) Q v 2 f v 2 ( g x 2 x 1 ) .
Task 2.
Arrangement of the expression for d V 2 / d t + k d V 3 / d t in terms of f v 1 .
The x 1 signal appearing in Equation (57) can be expressed as the sum of f v 1 and a disturbance term. Let,
d x 1 = f v 1 x 1 .
Using the definition of f v 1 (29), we have:
d x 1 = { x 1 u   f o r   x 1 > x 1 u x 1   f o r     x 1         [ x 1 l , x 1 u ] x 1 l   f o r   x 1 < x 1 l .
Hence,
d x 1 [ x 1 u , x 1 l ] .
From Equation (58), x 1 can be expressed as:
x 1 = f v 1 d x 1 .
Substituting in expression (57) and arranging, yields:
d V 2 d t + d d t ( v 3 v 2 V 3 ) = ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) + ( 1 ) Q v 2 f v 2 ( g x 2 δ 2 ) + f v 2 Q v 2 f v 1 .
where
δ 2 = d x 1 .
The effect of the term f v 2 ( Q / v 2 ) f v 1 is tackled later by considering the d V 1 / d t expression. □
Task 3.
Definition of gradient f v 2 = d V 2 / d x 2 and dead zone quadratic form V 2 .
At what follows, we examine the properties of δ 2 and the term g x 2 δ 2 appearing in Equation (62), what allows us to choose f v 2 and V 2 . From the definition of δ 2 (63) and d x 1 (59) it follows that δ 2 [ x 1 l ,   x 1 u ] .
Hence,
m a x { δ 2 } = x 1 u > 0 ,
m i n { δ 2 } = x 1 l < 0 .
The gradient of g x 2 (49) is:
d g x 2 d x 2 = 1 + 1 Q / v 2 d r ¯ x 2 d x 2 .
Furthermore, accounting for the definition of r x 2 (14), the definition of g x 2 in Equation (49), and the properties of r x 2 stated in assumption 3, we have:
d g x 2 d x 2 1 ,   d g x 2 d x 2   i s   b o u n d e d   f o r   x 2   b o u n d e d ,   a n d   s g n ( g x 2 ) = s g n ( x 2 ) .  
In turn, these properties lead to:
g x 2 δ 2 > 0       f o r   x 2 > x 2 u > 0 ,
g x 2 δ 2 < 0       f o r   x 2 < x 2 l < 0 .
where
x 2 u = { x 2 :   g x 2 m a x { δ 2 } = 0 } , x 2 u > 0 ,
x 2 l = { x 2 :   g x 2 m i n { δ 2 } = 0 } , x 2 l < 0 .
The term δ 2 is defined in Equations (63) and (59) and satisfies Equations (64) and (65). From Equations (67)–(69) it follows that:
s g n ( g x 2 δ 2 ) = s g n ( x 2 ) = s g n ( g x 2 )   f o r   x 2 [ x 2 l , x 2 u ] .
To obtain the non-positive nature of the term ( 1 ) f v 2 ( g x 2 δ 2 ) , appearing in Equation (62), we need to choose f v 2 such that:
f v 2 ( g x 2 δ 2 ) < 0     f o r   x 2 [ x 2 l , x 2 u ] ,
f v 2 ( g x 2 δ 2 ) = 0   f o r   x 2 [ x 2 l , x 2 u ] .
In view of properties (72) and to fulfill (73) and (74), we choose:
f v 2 = { x 2 x 2 u     f o r     x 2 > x 2 u > 0 0     f o r     x 2       [ x 2 l ,     x 2 u ] x 2 x 2 l     f o r     x 2 < x 2 l < 0 .
The main properties of f v 2 are:
( P i )   f v 2   i s   c o n t i n u o u s   w i t h   r e s p e c t   t o   x 2 ,
( P i i )   f v 2 = 0   f o r   x 2 [ x 2 l , x 2 u ] ,
( P i i i )   f v 2 0   f o r   x 2 [ x 2 l ,         x 2 u ] ,
( P i v )   s g n ( f v 2 ) = s g n ( x 2 ) 0   f o r   x 2 [ x 2 l , x 2 u ] ,
( P v )   s g n ( f v 2 ) = s g n ( g x 2 ) = s g n ( g x 2 δ 2 ) 0   f o r   x 2 [ x 2 l ,   x 2 u ] .
A Lyapunov function that satisfies d V 2 / d x 2 = f v 2 Equation (51) is:
V 2 = 1 2 f v 2 2 .
The main properties of V 2 are:
V 2 > 0   f o r       x 2     [ x 2 l , x 2 u ] ,
V 2 = 0   f o r   x 2   [ x 2 l , x 2 u ] ,
V 2   is   continuous   with   respect   to   x 2 .
Task 4.
Arrangement of the expression for d V 2 / d t + k d V 3 / d t in terms of a non-positive function of f v 2 = d V 2 / d x 2 .
From properties Pv (80) and Pii (77), it follows that the term f v 2 ( g x 2 δ 2 ) appearing in Equation (62) satisfies:
f v 2 ( g x 2 δ 2 ) = | f v 2 | | g x 2 δ 2 | = | f v 2 ( g x 2 δ 2 ) | 0 .
In addition, the term ( g x 2 δ 2 ) satisfies the following:
| g x 2 δ 2 | | g 2 v | ,
where,
g 2 v = { g x 2 m a x { δ 2 } f o r   x 2 > x 2 u 0     f o r       x 2     [ x 2 l ,     x 2 u ] g x 2 m i n { δ 2 } f o r   x 2 < x 2 l .
The main properties of g 2 v are:
( P i )   g 2 v   i s   c o n t i n u o u s   w i t h   r e s p e c t   t o   x 2 ,
( P i i )   g 2 v = 0   f o r   x 2     [ x 2 l ,   x 2 u ] ,
( P i i i )   g 2 v 0   f o r   x 2     [ x 2 l ,   x 2 u ] ,
( P i v )   s g n ( g 2 v ) = s g n ( x 2 ) 0   f o r   x 2     [ x 2 l ,   x 2 u ] .
From Equations (85) and (86), it follows that:
f v 2 ( g x 2 δ 2 ) = | f v 2 | | g x 2 δ 2 | | f v 2 g 2 v | 0 .
From the definitions of g 2 v (87) and f v 2 (75) it follows that:
| g 2 v | | f v 2 | .
Therefore, Equation (92) yields:
f v 2 ( g x 2 δ 2 ) f v 2 2 0 .
Using this, Equation (62), yields:
d V 2 d t + d d t ( v 3 v 2 V 3 ) ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) + ( 1 ) Q v 2 f v 2 2 + f v 2 Q v 2 f v 1 .
Task 5.
Definition of gradient d V 3 / d x 3 and dead zone quadratic form V 3 .
We choose f v 3 with the structure of f v 2 :
f v 3 = { x 3 x 3 u     f o r     x 3 > x 3 u > 0 0     f o r     x 3       [ x 3 l ,     x 3 u ] x 3 x 3 l     f o r     x 3 < x 3 l < 0 .
where the bounds x 3 u , x 3 l are defined as:
x 3 u = x 2 u ;   x 3 l = x 2 l   ,
and x 2 l , x 2 u are defined in Equations (70) and (71). A Lyapunov function that satisfies the gradient d V 3 / d x 3 = f v 3 (Equation (54)) is:
V 3 = 1 2 f v 3 2
The main properties of V 3 are:
V 3 > 0   f o r       x 3     [ x 3 l ,     x 3 u ] ,
V 3 = 0   f o r   x 3   [ x 3 l ,     x 3 u ] ,
V 3   is   continuous   with   respect   to   x 3 .
As a consequence of the definitions of f v 2 (75) and f v 3 (96), the property:
( x 2 x 3 ) ( f v 2 f v 3 ) 0 ,
holds true, as proved in Appendix A. Equation (95) can be arranged as:
d V 2 d t + d d t ( v 3 v 2 V 3 ) ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) β 1 Q v 2 f v 2 2 + ( 1 ) β 2 Q v 2 f v 2 2 + f v 2 Q v 2 f v 1 ,
where β 1 , β 2 are positive constants that satisfy:
1 = β 1 + β 2 ,     β 1 ( 0 ,   1 ) ,     β 2 ( 0 ,   1 ) .
Task 6.
Definition of the overall Lyapunov function V and determination of the non-positive nature of the expression for d V / d t = d V 1 / d t + k d V 2 / d t + d V 3 / d t .
The term ( 1 ) β 2 Q v 2 f v 2 2 + f v 2 Q v 2 f v 1 appearing in Equation (103) can be expressed in terms of f v 1 2 :
( 1 ) β 2 Q v 2 f v 2 2 + f v 2 Q v 2 f v 1 = ( 1 ) β 2 Q v 2 ( f v 2 1 2 β 2 f v 1 ) 2 + 1 4 β 2 Q v 2 f v 1 2 1 4 β 2 Q v 2 f v 1 2 .
Substituting into Equation (103), we have:
d V 2 d t + d d t ( v 3 v 2 V 3 ) ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) + ( 1 ) β 1 Q v 2 f v 2 2 + 1 4 β 2 Q v 2 f v 1 2 .
We tackle the effect of the f v 1 2 term by incorporating the d V 1 / d t expression. Equation (44) leads to:
d d t ( 1 D 1 4 β 2 Q v 2 V 1 ) 1 4 β 2 Q v 2 f v 1 2 .
Combining with Equation (106) and accounting for the property (102), yields:
d V 2 d t + d d t ( v 3 v 2 V 3 ) + d d t ( 1 D 1 4 β 2 Q v 2 V 1 ) ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) + ( 1 ) β 1 Q v 2 f v 2 2 0 ,
which can be rewritten as:
d V d t ( 1 ) Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) + ( 1 ) β 1 Q v 2 f v 2 2 0 .
where V is the overall Lyapunov function, defined as:
V = V 2 + v 3 v 2 V 3 + 1 D 1 4 β 2 Q v 2 V 1 .
Task 7.
Integration of the d V / d t expression.
Integrating Equation (108), yields:
V 2 + v 3 v 2 V 3 + 1 4 β 2 1 D Q v 2 V 1 + β 1 Q v 2 t o t f v 2 2 d t + t o t Q D v 2 ( x 2 x 3 ) ( f v 2 f v 3 ) d t V 2 t o + v 3 v 2 V 3 t o + 1 4 β 2 1 D Q v 2 V 1 t o .
Therefore: (BPi) V 2 , V 3 ,   V 1 , and from the definitions of V 1 (35), V 2 (81), V 3 (98) it follows that f v 1 , f v 2 , f v 3 , and consequently f v 2 2 ; (BPii) x 1 , x 2 , x 3 , what follows from the definitions of f v 1 (29), f v 2 (75), f v 3 (96); (BPiii) f v 2 2 1 . Considering the properties f v 2 2 1   and f v 2 2 and applying the Barbalat’s lemma (cf. [36]), yields:
lim t f v 2 2 = 0 .
Furthermore, accounting for the definition of f v 2 (75), it follows that x 2 converges asymptotically to Ω x 2 = { x 2 : x 2 l x 2 x 2 u } . From this convergence result and the definition of x2 (9), it follows that S2 converges asymptotically to Ω S 2 ; where:
Ω S 2 = [ S 2 l     S 2 u ]   ,
S 2 l = x 2 l + S 2 e q ,   S 2 u = x 2 u + S 2 e q   .
S 2 e q is provided by Equations (5)–(7); the bounds x 2 l , x 2 u are defined in Equations (70) and (71). This completes the proof of the first part of the theorem (Ti).
From Equation (111), it follows that ( x 2 x 3 ) ( f v 2 f v 3 ) L 1 and from properties Bpi, BPii it follows that ( x 2 x 3 ) ( f v 2 f v 3 ) L . Applying the Lasalle’s theorem (cf. [22]) yields lim t ( x 2 x 3 ) ( f v 2 f v 3 ) = 0 .
Hence, either x 2 x 3 or f v 2 f v 3 converges to zero. This result and the convergence of x 2 to Ω x 2 , the convergence of f v 2 to zero, and the definition of f v 3 (96) and f v 2 (75), imply: (i) if x 2 x 3 converges to zero, then x 3 converges to Ω x 3 , Ω x 3 = { x 3 : x 3 l x 3 x 3 u } ; (ii) if f v 2 f v 3 converges to zero, then f v 3 converges to zero, which implies the convergence of x 3 to Ω x 3 . Therefore, as is indicated by both of these possibilities, x 3 converges asymptotically to Ω x 3 . From this convergence result and the definition of x 3 (10) it follows that S 3 converges asymptotically to Ω S 3 , where:
Ω S 3 = [ S 3 l S 3 u ] ,
S 3 l = x 3 l + S 3 e q ,   S 3 u = x 3 u + S 3 e q   .
S 3 e q is provided by Equations (5)–(7), and the bounds x 3 l , x 3 u are defined in Equation (97). This completes the proof of the second part of the theorem (Tii).
From Equation (109), it follows that:
i f   V = 0   f o r   t = t * ,       t * t o ,     t h e n   V = 0   f o r   t t * .
The condition V = 0   implies: (i) f v 1 = 0 , f v 2 = 0 and f v 3 = 0 , as follows from definitions of V (110), V1 (35), V2 (81), and V3 (98); (ii) x 1 Ω x 1 , x 2 Ω x 2 , x 3 Ω x 3 , as follows from the definition of f v 1 (29), f v 2 (75) and f v 3 (96). Therefore, it follows from (117) that Ω x = Ω x 1   Ω x 2   Ω x 3 is invariant, and consequently Ω s = Ω s 1   Ω s 2   Ω s 3   is invariant. This completes the proof of Tiii. □
Remark 1.
Important advantages of the stability result stated in theorem 2 are:
  • the widths of the convergence regions Ωx1, Ωx2, Ωx3 are arbitrarily large, as they are a function of δ s i , which is arbitrarily large but bounded;
  • the initial values of x 1 , x 2 , and x 3 can take arbitrarily large but bounded positive values, since x 1 is defined in [− S 1 e q ,   ) ,   x 2 is defined in [− S 2 e q ,   ) , and x 3 is defined in [− S 3 e q ,   ) , as stated in Section 2.
Remark 2.
In the study of the stability of x 1 , the dynamics of x 2 and x 3 (2), (3) are not accounted for, as can be noticed in the proof of Theorem 1. This is a consequence of the fact that x 2 and x 3 are not involved in the dynamics of x 1 (11). In contrast, in the study of the stability of x 2 and x 3 , the dynamics of x 1 , x 2 and x 3 are accounted for, as can be noticed in the proof of theorem 2. This is a consequence of the fact that x 1 , x 2 and x 3 are involved in the dynamics of x 2 (13).
Remark 3.
The first tank influences the transient value of x 2 , which can be noticed from Equation (111). Also, the first tank influences the width of the convergence region of x 2 , that is [ x 2 l , x 2 u ], what can be noticed from the definition of x 2 l , x 2 u (70), (71), the definition of δ2 (63), (58), the properties of δ 2 (64), (65), and the definition of x 1 l , x 1 u (24), (25).
Remark 4.
The presence of the third tank influences the transient value of x 2 , what can be noticed from Equation (111). However, it does not influence the width of the convergence region of the state x 2 , that is, [ x 2 l , x 2 u ], since the bounds ( x 2 l , x 2 u ) are not influenced by the parameters of the third tank, that is, Q D and v 3 , which can be noticed from the definition of x 2 l , x 2 u (70), (71).
Remark 5.
Important tasks of the stability analysis are: the definition of the overall Lyapunov-like function V (110) consisting on the weighted sum of three dead-zoned quadratic forms; the determination of the property (102), which expresses the effect of the flow exchange on the non-positive nature of d V / d t ; the property | g 2 v | | f 2 v | (93) which was required for proving the convergence of x 2 . These tasks are crucial for proving the non-positive nature of d V / d t . Also, the whole stability analysis and especially these tasks are major contributions to the study of convergence of biological processes to compact sets.
Remark 6.
The invariance of the set Ω s stated in Theorem 2 implies that once the three states S1, S2, S3 are inside the convergence region Ω s , that is, S 1 Ω s 1 , S 2 Ω s 2 and S 3 Ω s 3 simultaneously, they remain inside afterwards.
Remark 7.
The developed stability analysis can be extended to other types of nonlinear connected systems, featuring a higher number of state variables, and a vector field with different non-linear terms, by using the following general procedure: (i) determination of the equilibrium points corresponding to the case of no time varying external disturbance; (ii) definition of the new states as the differences between the current state variables and their equilibria; (iii) rewriting of the system dynamics in terms of the new states; (iv) definition of the subsystem Lyapunov functions, corresponding to each state variable, and determination of its properties; (v) determination of the time derivative of each subsystem Lyapunov function, arrangement in terms of non-positive functions and definition of the convergence regions; (vi) definition of the overall Lyapunov function V as a weighted sum of the subsystem Lyapunov functions; (vii) determination of the dV/dt expression and arrangement in terms of non-positive functions; (viii) integration of the dV/dt expression, determination of the boundedness properties of the state variables and their functions; (ix) application of the Barbalat’s lemma. As part of this procedure: the subsystem Lyapunov functions can be defined as dead-zone quadratic forms, or they can be defined according to the non-linear terms of the vector field; the subsystem Lyapunov functions corresponding to the state variables featuring connection must be chosen so as to obtain a non-positive effect in the dV/dt expression; if there are one or more state variables whose vector field is independent of the other states, an independent stability analysis can be performed for each of them.

4. Simulation and Analysis

In this section, some simulations of the system (1)–(3) are performed to illustrate the results stated in Theorems 1 and 2.
We consider the following reaction rate expressions:
r s 1 = μ m x 1 S 1 K 1 + S 1 ,   h e n c e   r s 1 = μ m x 1 x 1 + S 1 e q K 1 + x 1 + S 1 e q , r s 2 = μ m x 2 S 2 K 2 + S 2 ,   h e n c e   r s 2 = μ m x 2 x 2 + S 2 e q K 2 + x 2 + S 2 e q .
The terms r s 1 and r s 2 satisfy assumption 3. The model (1)–(3) with these reaction rate expressions was calibrated with data of suspended solids concentration from vertical subsurface flow constructed wetlands (CWs) reported in [37]. The CWs are located in Cuenca (Ecuador) and receive domestic wastewater coming out from primary treatment. The influent wastewater comprises a suspended solids concentration of 88.83 mg/L, BOD5 concentration of 95.75 mg/L, and temperature of 24.6 °C. The hydraulic loading rate (HLR) is 0.2 md−1.
Model (1)–(3) was fitted by minimizing the sum of the square residuals between experimental and model data of suspended solids concentration. The parameter values obtained were used in the first simulation case and are shown in Table 1.
To evaluate the obtained convergence results stated in Theorems 1 and 2, model (1)–(4) is simulated with
S ¯ i n = 100   mg / L ;   δ s i n = A 1 s i n ( 2 π t T s i n ) ,   A 1 = 10 ;   T s i n = 20   days .
The bounds S 1 l , S 1 u of Ω S 1 (47); the bounds S 2 l , S 2 u of Ω S 2 (114); and the bounds S 3 l , S 3 u of Ω S 3 (116) are computed and represented by the upper and lower horizontal dotted lines, whereas the equilibrium values S 1 e q , S 2 e q and S 3 e q (Equations (5)–(7)) are represented by the middle horizontal dotted lines. The calculation of the aforementioned bounds requires the calculation of bounds x 1 l , x 1 u (24), (25); x 2 l , x 2 u (70), (71); and x 3 l , x 3 u (97). Three simulations are performed, using definitions (118) and parameter values shown in Table 1.
In the three simulations, convergence to the regions within the calculated bounds is observed: the state S 1 convergences to the region within computed bounds S 1 l , S 1 u (see Figure 2a, Figure 3a, Figure 4a,b); the state S 2 converges to the region within computed bounds S 2 l , S 2 u (see Figure 2b, Figure 3b, Figure 4c,d); and S 3 converges to the region within the computed bounds S 3 l , S 3 u (see Figure 2b, Figure 3b, Figure 4c,d).
Also, the three simulation cases illustrate the results stated in Theorems 1 and 2: (i) the invariant nature of Ω s 1 , stated in Theorem 1 and discussed in Remark 2: one can see that once S1 enters Ω s 1 , it remains inside; (ii) the equivalence of the convergence regions Ω s 2 and Ω s 3 defined in Theorem 2: one can see that S2 and S3 converge to the same regions; (iii) the definitions of Ω s 1 and Ω s 2 , which indicate that Ω s 1 and Ω s 2 can be quite different, as Ω s 2 depends on the reaction rate rx2: one can see this remarkable difference in the regions to which S1 and S2 converge; (iv) the invariant nature of Ω s , stated in Theorem 2 and discussed in Remark 6: one can see that once the three state variables are inside the convergence region Ω s , that is, S 1 Ω s 1 , S 2 Ω s 2 , S 3 Ω s 3 , simultaneously, they remain inside afterwards.

5. Conclusions

This paper presented the analysis of the stability of a network model comprising three CSTRs with the following features: (i) the pollutant concentration in the inflow of the first CSTR is time varying but bounded; (ii) the first and second CSTRs are connected in series, whereas the second and third CSTRs are connected in parallel with flow exchange; (iii) the states converge to a compact set instead of an equilibrium point. The practical applicability of the arrangement of CSTRs is the creation of a simpler model of pollution removal from wastewater treatment via constructed wetlands, generating a satisfactory description of experimental pollution values with a satisfactory transport dead time.
The convergence sets of the states of the CSTR model were determined, and the global asymptotic convergence to these compact sets and their invariance were proved. Also, the effect of the side tank (third tank) on the transient value of the system states was determined, and it was concluded that it had no effect on the convergence regions. To this end, the proposed stability analysis comprises:
(i)
the definition of a dead-zone Lyapunov-like function for each state, the determination of its properties, and the definition of the overall Lyapunov function as the sum of the dead zone quadratic forms;
(ii)
the determination of the time derivative of the quadratic forms and their properties;
(iii)
the use of these properties in the time derivative of the overall Lyapunov-like function.

Author Contributions

Conceptualization, A.R.; methodology, A.R.; validation, F.E.H. and J.E.C.-B.; formal analysis, A.R., F.E.H. and J.E.C.-B.; investigation, A.R.; writing—original draft preparation, A.R.; writing–review and editing, A.R., F.E.H. and J.E.C.-B.; visualization, F.E.H. and J.E.C.-B. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The work of A. Rincón was supported by Universidad Católica de Manizales. The work of F.E. Hoyos and J.E. Candelo-Becerra was supported by Universidad Nacional de Colombia–Sede Medellín.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The expression ( x 2 x 3 ) ( f v 2 f v 3 ) can be expressed as:
( x 2 x 3 ) ( f v 2 f v 3 ) = | ( x 2 x 3 ) ( f v 2 f v 3 ) |   s g n ( x 2 x 3 )   s g n ( f v 2 f v 3 ) ,
or
( x 2 x 3 ) ( f v 2 f v 3 ) = = | ( x 2 x 3 ) ( f v 2 f v 3 ) |   s g n ( x 3 x 2 )   s g n ( f v 3 f v 2 ) .
Thus, the signum of ( x 2 x 3 ) ( f v 2 f v 3 ) depends on the value of s g n ( x 2 x 3 )   s g n ( f v 2 f v 3 ) , or, equivalently, the signum of s g n ( x 3 x 2 )   s g n ( f v 3 f v 2 ) . To determine it, three different cases are analyzed as follows.
Case 1. | f v 2 | 0 ,   | f v 3 | 0 . In this case,
f v 2 0 , f v 3 0 .
From the definitions of f v 2 (75) and f v 3 (96) it follows that:
x 2     [ x 2 l , x 2 u ] , x 3     [ x 3 l , x 3 u ] .
Hence,
s g n ( f v 2 ) = s g n ( x 2 ) 0 ,
s g n ( f v 3 ) = s g n ( x 3 ) 0 .
Therefore, accounting for the definition of f v 2 and f v 3 , we have:
s g n ( f v 2 f v 3 ) = s g n ( x 2 x 3 ) .
Therefore, from Equation (A1), it follows that:
( x 2 x 3 ) ( f v 2 f v 3 ) 0 .
Case 2. f v 2 = 0 ,   f v 3 0 . In this case,
x 2     [ x 2 l ,   x 2 u ] , x 3     [ x 3 l ,   x 3 u ] .
Moreover,
s g n ( f v 2 ) = 0 ,
s g n ( f v 3 ) = s g n ( x 3 ) 0 ,
s g n ( x 3 x 2 ) = s g n ( x 3 ) 0 .
Furthermore, accounting for the definition of f v 2 (75) and f v 3 (96), we have:
s g n ( f v 3 f v 2 ) = s g n ( f v 3 ) = s g n ( x 3 ) 0 .
Therefore,
s g n ( f v 3 f v 2 ) s g n ( x 3 x 2 ) = + 1 .
Therefore, from Equation (A2), it follows that:
( x 2 x 3 ) ( f v 2 f v 3 ) 0 .
Case 3. f v 2 0 ,   f v 3 = 0 . In this case, with a procedure similar to that used for case 2, it is found that:
( x 2 x 3 ) ( f v 2 f v 3 ) 0 .

References

  1. Zhou, G.; Liao, X.; Xu, B.; Yu, P.; Chen, G. Simple algebraic necessary and sufficient conditions for Lyapunov stability of a Chen system and their applications. Trans. Inst. Meas. Control 2018, 40, 2200–2210. [Google Scholar] [CrossRef]
  2. Liao, X.; Zhou, G.; Yang, Q.; Fu, Y.; Chen, G. Constructive proof of Lagrange stability and sufficient—Necessary conditions of Lyapunov stability for Yang–Chen chaotic system. Appl. Math. Comput. 2017, 309, 205–221. [Google Scholar] [CrossRef]
  3. Zhang, F.; Liao, X.; Zhang, G.; Mu, C. Dynamical Analysis of the Generalized Lorenz Systems. J. Dyn. Control Syst. 2017, 23, 349–362. [Google Scholar] [CrossRef]
  4. Meng, R.; Chen, S.; Hua, C.; Qian, J.; Sun, J. Disturbance observer-based output feedback control for uncertain QUAVs with input saturation. Neurocomputing 2020, 413, 96–106. [Google Scholar] [CrossRef]
  5. Boulkroune, A.; M’saad, M.; Farza, M. Adaptive fuzzy system-based variable-structure controller for multivariable nonaffine nonlinear uncertain systems subject to actuator nonlinearities. Neural Comput. Appl. 2017, 28, 3371–3384. [Google Scholar] [CrossRef]
  6. Zhang, L.; Wang, Z.; Li, S.; Ding, S.; Du, H. Universal finite-time observer based second-order sliding mode control for DC-DC buck converters with only output voltage measurement. J. Franklin Inst. 2020, 357, 11863–11879. [Google Scholar] [CrossRef]
  7. Shao, X.; Wang, L.; Li, J.; Liu, J. High-order ESO based output feedback dynamic surface control for quadrotors under position constraints and uncertainties. Aerosp. Sci. Technol. 2019, 89, 288–298. [Google Scholar] [CrossRef]
  8. Askari, M.R.; Shahrokhi, M.; Khajeh Talkhoncheh, M. Observer-based adaptive fuzzy controller for nonlinear systems with unknown control directions and input saturation. Fuzzy Sets Syst. 2017, 314, 24–45. [Google Scholar] [CrossRef]
  9. Gao, S.; Ning, B.; Dong, H. Fuzzy dynamic surface control for uncertain nonlinear systems under input saturation via truncated adaptation approach. Fuzzy Sets Syst. 2016, 290, 100–117. [Google Scholar] [CrossRef]
  10. Min, H.; Xu, S.; Ma, Q.; Zhang, B.; Zhang, Z. Composite-Observer-Based Output-Feedback Control for Nonlinear Time-Delay Systems With Input Saturation and Its Application. IEEE Trans. Ind. Electron. 2018, 65, 5856–5863. [Google Scholar] [CrossRef]
  11. Zhang, F.; Liao, X.; Zhang, G. On the global boundedness of the Lü system. Appl. Math. Comput. 2016, 284, 332–339. [Google Scholar] [CrossRef]
  12. Slotine, J.-J.E.; Li, W. Applied Nonlinear Control; Prentice Hall: Englewood Cliffs, NJ, USA, 1991; ISBN 0130408905. [Google Scholar]
  13. Koo, K.-M. Stable adaptive fuzzy controller with time-varying dead-zone. Fuzzy Sets Syst. 2001, 121, 161–168. [Google Scholar] [CrossRef]
  14. Wang, X.-S.; Su, C.-Y.; Hong, H. Robust adaptive control of a class of nonlinear systems with unknown dead-zone. Automatica 2004, 40, 407–413. [Google Scholar] [CrossRef]
  15. Su, C.-Y.; Feng, Y.; Hong, H.; Chen, X. Adaptive control of system involving complex hysteretic nonlinearities: A generalised Prandtl–Ishlinskii modelling approach. Int. J. Control 2009, 82, 1786–1793. [Google Scholar] [CrossRef]
  16. Ranjbar, E.; Yaghubi, M.; Abolfazl Suratgar, A. Robust adaptive sliding mode control of a MEMS tunable capacitor based on dead-zone method. Automatika 2020, 61, 587–601. [Google Scholar] [CrossRef]
  17. Rincón, A.; Hoyos, F.E.; Candelo-Becerra, J.E. Adaptive Control for a Biological Process under Input Saturation and Unknown Control Gain via Dead Zone Lyapunov Functions. Appl. Sci. 2021, 11, 251. [Google Scholar] [CrossRef]
  18. Rincón, A.; Piarpuzán, D.; Angulo, F. A new adaptive controller for bio-reactors with unknown kinetics and biomass concentration: Guarantees for the boundedness and convergence properties. Math. Comput. Simul. 2015, 112, 1–13. [Google Scholar] [CrossRef]
  19. Rincón, A.; Florez, G.Y.; Olivar, G. Convergence Assessment of the Trajectories of a Bioreaction System by Using Asymmetric Truncated Vertex Functions. Symmetry 2020, 12, 513. [Google Scholar] [CrossRef] [Green Version]
  20. Meadows, T.; Weedermann, M.; Wolkowicz, G.S.K. Global Analysis of a Simplified Model of Anaerobic Digestion and a New Result for the Chemostat. SIAM J. Appl. Math. 2019, 79, 668–689. [Google Scholar] [CrossRef]
  21. Mairet, F.; Ramírez, C.H.; Rojas-Palma, A. Modeling and stability analysis of a microalgal pond with nitrification. Appl. Math. Model. 2017, 51, 448–468. [Google Scholar] [CrossRef]
  22. Sari, T.; Mazenc, F. Global dynamics of the chemostat with different removal rates and variable yields. Math. Biosci. Eng. 2011, 8, 827–840. [Google Scholar] [CrossRef] [PubMed]
  23. Almeida, R.; Hristova, S.; Dashkovskiy, S. Uniform bounded input bounded output stability of fractional-order delay nonlinear systems with input. Int. J. Robust Nonlinear Control 2021, 31, 225–249. [Google Scholar] [CrossRef]
  24. Garrido, H.; Curadelli, O.; Ambrosini, D. On the assumed inherent stability of semi-active control systems. Eng. Struct. 2018, 159, 286–298. [Google Scholar] [CrossRef]
  25. Vafamand, N.; Asemani, M.H.; Khayatian, A. Robust L1 Observer-Based Non-PDC Controller Design for Persistent Bounded Disturbed TS Fuzzy Systems. IEEE Trans. Fuzzy Syst. 2018, 26, 1401–1413. [Google Scholar] [CrossRef]
  26. Kong, L.; He, W.; Dong, Y.; Cheng, L.; Yang, C.; Li, Z. Asymmetric Bounded Neural Control for an Uncertain Robot by State Feedback and Output Feedback. IEEE Trans. Syst. Man Cybern. Syst. 2019, 1–12. [Google Scholar] [CrossRef] [Green Version]
  27. Yu, J.; Zhao, L.; Yu, H.; Lin, C. Barrier Lyapunov functions-based command filtered output feedback control for full-state constrained nonlinear systems. Automatica 2019, 105, 71–79. [Google Scholar] [CrossRef]
  28. Niu, B.; Liu, Y.; Zhou, W.; Li, H.; Duan, P.; Li, J. Multiple Lyapunov Functions for Adaptive Neural Tracking Control of Switched Nonlinear Nonlower-Triangular Systems. IEEE Trans. Cybern. 2020, 50, 1877–1886. [Google Scholar] [CrossRef] [PubMed]
  29. De Battista, H.; Jamilis, M.; Garelli, F.; Picó, J. Global stabilisation of continuous bioreactors: Tools for analysis and design of feeding laws. Automatica 2018, 89, 340–348. [Google Scholar] [CrossRef]
  30. El Hajji, M.; Chorfi, N.; Jleli, M. Mathematical modeling and analysis for a three-tiered microbial food web in a chemostat. Electron. J. Differ. Equ. 2017, 2017, 1–13. [Google Scholar]
  31. Marsili-Libelli, S.; Checchi, N. Identification of dynamic models for horizontal subsurface constructed wetlands. Ecol. Modell. 2005, 187, 201–218. [Google Scholar] [CrossRef]
  32. Langergraber, G. Applying Process-Based Models for Subsurface Flow Treatment Wetlands: Recent Developments and Challenges. Water 2016, 9, 5. [Google Scholar] [CrossRef] [Green Version]
  33. Defo, C.; Kaur, R.; Bharadwaj, A.; Lal, K.; Kumar, P. Modelling approaches for simulating wetland pollutant dynamics. Crit. Rev. Environ. Sci. Technol. 2017, 47, 1371–1408. [Google Scholar] [CrossRef]
  34. Davies, L.C.; Vacas, A.; Novais, J.M.; Freire, F.G.; Martins-Dias, S. Vertical flow constructed wetland for textile effluent treatment. Water Sci. Technol. 2007, 55, 127–134. [Google Scholar] [CrossRef] [PubMed]
  35. Freire, F.G.; Davies, L.C.; Vacas, A.M.; Novais, J.M.; Martins-Dias, S. Influence of operating conditions on the degradation kinetics of an azo-dye in a vertical flow constructed wetland using a simple mechanistic model. Ecol. Eng. 2009, 35, 1379–1386. [Google Scholar] [CrossRef]
  36. Ioannou, P.A.; Sun, J. Robust Adaptive Control; Prentice Hall: Upper Saddle River, NJ, USA, 1996; ISBN 9780486498171. [Google Scholar]
  37. García-Ávila, F.; Patiño-Chávez, J.; Zhinín-Chimbo, F.; Donoso-Moscoso, S.; Flores del Pino, L.; Avilés-Añazco, A. Performance of Phragmites Australis and Cyperus Papyrus in the treatment of municipal wastewater by vertical flow subsurface constructed wetlands. Int. Soil Water Conserv. Res. 2019, 7, 286–296. [Google Scholar] [CrossRef]
Figure 1. Schematic structure of the series/parallel continuously stirred tank reactor (CSTR)-based model.
Figure 1. Schematic structure of the series/parallel continuously stirred tank reactor (CSTR)-based model.
Applsci 11 04178 g001
Figure 2. Simulation of the CSTR network model (1)–(4), first case. (a) Time course of S 1 ; the upper and lower horizontal dotted lines indicate the bounds of Ω S 1 . (b) Time course of S 2 and S 3 , the upper and lower horizontal dotted lines indicate the bounds of Ω S 2 , and Ω S 3 .
Figure 2. Simulation of the CSTR network model (1)–(4), first case. (a) Time course of S 1 ; the upper and lower horizontal dotted lines indicate the bounds of Ω S 1 . (b) Time course of S 2 and S 3 , the upper and lower horizontal dotted lines indicate the bounds of Ω S 2 , and Ω S 3 .
Applsci 11 04178 g002
Figure 3. Simulation of the CSTR network model (1)–(4), second case. (a) Time course of S 1 ; the upper and lower horizontal dotted lines indicate the bounds of Ω S 1 . (b) Time course of S 2 and S 3 , the upper and lower horizontal dotted lines indicate the bounds of Ω S 2 and Ω S 3 .
Figure 3. Simulation of the CSTR network model (1)–(4), second case. (a) Time course of S 1 ; the upper and lower horizontal dotted lines indicate the bounds of Ω S 1 . (b) Time course of S 2 and S 3 , the upper and lower horizontal dotted lines indicate the bounds of Ω S 2 and Ω S 3 .
Applsci 11 04178 g003
Figure 4. Simulation of the CSTR network model (1)–(4), third case. (a) Time course of S 1 ; the upper and lower horizontal dotted lines indicate the bounds of Ω S 1 . (b) Detail of the time course of S 1 . (c) Time course of S 2 and S 3 , the upper and lower horizontal dotted lines indicate the bounds of Ω S 2 and Ω S 3 . (d) Detail of the time course of S 2 and S 3 .
Figure 4. Simulation of the CSTR network model (1)–(4), third case. (a) Time course of S 1 ; the upper and lower horizontal dotted lines indicate the bounds of Ω S 1 . (b) Detail of the time course of S 1 . (c) Time course of S 2 and S 3 , the upper and lower horizontal dotted lines indicate the bounds of Ω S 2 and Ω S 3 . (d) Detail of the time course of S 2 and S 3 .
Applsci 11 04178 g004
Table 1. Parameter values for the CSTR network model (1)–(3).
Table 1. Parameter values for the CSTR network model (1)–(3).
ParameterValue for the First Simulation CaseValue for the Second Simulation CaseValue for the Third Simulation Case
Q1/v10.702 day−10.702 day−10.702 day−1
Q2/v20.702 day−10.702 day−10.702 day−1
QD/v20.0001 day−10.0001 day−10.0001 day−1
QD/v31 day−11 day−11 day−1
μmax1101 mg /(L day)101 mg /(L day)101 mg /(L day)
K1303.1 mg/L303.1 mg/L303.1 mg/L
μmax2101 mg/(L day)202 mg/(L day)404 mg/(L day)
K2303.1 mg/L151.55 mg/L30.31 mg/L
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rincón, A.; Hoyos, F.E.; Candelo-Becerra, J.E. Global Stability Analysis of the Model of Series/Parallel Connected CSTRs with Flow Exchange Subject to Persistent Perturbation on the Input Concentration. Appl. Sci. 2021, 11, 4178. https://doi.org/10.3390/app11094178

AMA Style

Rincón A, Hoyos FE, Candelo-Becerra JE. Global Stability Analysis of the Model of Series/Parallel Connected CSTRs with Flow Exchange Subject to Persistent Perturbation on the Input Concentration. Applied Sciences. 2021; 11(9):4178. https://doi.org/10.3390/app11094178

Chicago/Turabian Style

Rincón, Alejandro, Fredy E. Hoyos, and John E. Candelo-Becerra. 2021. "Global Stability Analysis of the Model of Series/Parallel Connected CSTRs with Flow Exchange Subject to Persistent Perturbation on the Input Concentration" Applied Sciences 11, no. 9: 4178. https://doi.org/10.3390/app11094178

APA Style

Rincón, A., Hoyos, F. E., & Candelo-Becerra, J. E. (2021). Global Stability Analysis of the Model of Series/Parallel Connected CSTRs with Flow Exchange Subject to Persistent Perturbation on the Input Concentration. Applied Sciences, 11(9), 4178. https://doi.org/10.3390/app11094178

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop