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

: 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 ﬁrst and second CSTRs are connected in series, whereas the second and third CSTRs are connected in parallel with ﬂow exchange; (ii) the pollutant concentration in the inﬂow to the ﬁrst CSTR is time varying but bounded; (iii) the states converge to a compact set instead of an equilibrium point, due to the time varying inﬂow 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) deﬁnition 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.


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 , 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 openloop 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 L 2 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 L 2 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 (V 2 ) comprises the subsystem Lyapunov functions associated to the z i states and the parameter updating error. The dV 2 /dt expression is a function of V 2 , 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.

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 = (x 1 , x 2 , . . . , x n ) ∈ R n , f: R n → R n , X(t,t o ,X o ) denoted as X(t) is a solution to the system, and X 0 ∈ R n is an initial value of X(t). Consider the compact set Ω Q ⊂ R n . Define the distance between the solution X and the compact set Ω Q as: If for all X 0 satisfying X 0 ∈ R n , X 0 / ∈ Ω 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: 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 in 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 s1 is the reaction rate of CSTR 1; r s2 is the reaction rate of CSTR 2. In addition, S 1 is defined in R + , S 2 is defined in R + , S 3 is defined in R + .
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of 21 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: ℛ n →ℛ n , X(t,to,Xo) denoted as X(t) is a solution to the system, and X0 ∈ ℛ n is an initial value of X(t). Consider the compact set ΩQ ⊂ ℛ n . Define the distance between the solution X and the compact set ΩQ as: If for all X0 satisfying X0 ∈ ℛ n , X0 ∉ ΩQ the property lim →∞ ( , ) = 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: where 1 is the pollutant concentration in CSTR 1; 2 is the pollutant concentration in CSTR 2; 3 is the pollutant concentration in CSTR 3; is the inlet pollutant concentration; is the wastewater flow entering CSTR 1 and CSTR 2; is the wastewater flow entering and leaving CSTR 3; is the liquid volume of the ith CSTR; 1 is the reaction rate of CSTR 1; 2 is the reaction rate of CSTR 2. In addition, 1 is defined in ℝ + , 2 is defined in ℝ + , 3 is defined in ℝ + .  Assumption 2. S in is positive and constant,δ sin is varying but bounded, max{δ sin } > 0, min{δ sin } < 0 .

Assumption 3.
The reaction rate r s1 only depends on S 1 , it is bounded for S 1 bounded, dr s1 dS 1 ≥ 0; r s2 only depends on S 2 , it is bounded for S 2 bounded, and dr s2 dS 2 ≥ 0.
In the case of constant S in , δ sin = 0 in Equation (4), S 1 converges to S eq 1 , S 2 converges to S eq 2 and S 3 converges to S eq 3 , where S eq 1 , S eq 2 and S eq 3 are positive constants obtained by equating Equations (1)-(3) to zero: Let The x 1 dynamics is obtained by subtracting Equation (5) from Equation (1): where D = Q/v 1 > 0, r eq s1 is r s1 evaluated at S 1 = S eq 1 . The x 2 dynamics is obtained by subtracting Equation (6) from Equation (2): The term r eq s2 is r s2 evaluated at S 2 = S eq 2 . The x 3 dynamics is obtained by subtracting Equation (7) from Equation (3): In addition, . The model is fitted to current (no tracer) data in Section 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.

Stability Analysis for the First Tank
The stability analysis is based on determining the non-positive nature of dV 1 /dt, where V 1 is the dead-zone quadratic form for x 1 . The function V 1 and its gradient dV 1 /dx 1 are defined so as to achieve this goal.
The procedure comprises the following tasks: (i) determination of the expression for dV 1 /dt; (ii) definition of the gradient dV 1 /dx 1 and the dead zone quadratic form V 1 for the state x 1 ; (iii) arrangement of the dV 1 /dt expression in terms of a non-positive function of dV 1 /dx 1 ; (iv) integration of the dV 1 /dt expression. The definition of the gradient dV 1 /dx 1 . and V 1 requires the determination of the dV 1 /dt expression and the properties of the terms involved. (11) and (12), subject to assumptions 1 to 3.

Theorem 1. Consider the model in Equations
(Ti) The state x 1 converges asymptotically to Ω x1 , D , and S 1 converges asymptotically to Ω S1 where Ω S1 = S l and S eq 1 is provided by Equation (5). (Tii) The sets Ω x1 and Ω S1 are invariant.
Proof. Task 1. Determination of the dV 1 /dt expression. The time derivative of V 1 can be expressed as: Incorporating the x 1 dynamics (11) and arranging, yields: where, Task 2. Definition of the gradient dV 1 /dx 1 and the dead zone quadratic form V 1 for the state At what follows, we examine the properties of the term g x1 − δ sin appearing in Equation (18). From Equations (12), (19) and (20) and Assumption 3 we have: where, Therefore, To obtain non-positive dV 1 /dt in Equation (18), we need to choose f v1 such that: In view of properties (26) and to fulfill conditions (27) and (28), we choose: The main properties of f v1 are: A dead-zone Lyapunov function that satisfies dV 1 /dx 1 = f v1 (Equation (17)) is: The main properties of V 1 are: Task 3. Arrangement of the dV 1 /dt expression in terms of a non-positive function of From properties (31), (34), it follows that the term f v1 (g x1 − δ sin ) appearing in Equation (18) satisfies: From the definition of g x1 (19) and definitions of x l 1 , x u 1 (24), (25), it follows that: From Equations (38) and (39), it follows that: From the property (21), definition of g t1 (40), and definition of f v1 (29), it follows that: Therefore, Equation (41) yields: Using this in Equation (18) yields: Accounting for the definition of V 1 (35), we have:

Task 4.
Integration of the dV 1 /dt expression. From the above expression, it follows that: That is, V 1 (35) converges exponentially to zero, and consequently f v1 converges exponentially to zero. Further, accounting for the definition of f v1 (29), it follows that x 1 converges asymptotically to Ω x1 = x 1 : x l 1 ≤ x 1 ≤ x u 1 . From this convergence result and the definition of x 1 (8), it follows that S 1 converges asymptotically to Ω S1 ; where, S eq 1 is provided by Equation (5), the bounds x l 1 , x u 1 are defined in Equations (24) and (25). This completes the proof of Ti.

Remark 1.
The property |g t | ≥ | f v1 | (42) is crucial for proving the asymptotic convergence of x 1 .

Remark 2.
The invariance of the set Ω s1 stated in Theorem 1 implies that once the state S 1 is inside the convergence region Ω s1 , that is, S 1 ∈ Ω s1 , it remains inside afterwards.

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 dV/dt = dV 1 /dt + k 2 dV 2 /dt + k 3 dV 3 /dt, where V = V 1 + kV 2 + kV 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 dV 2 /dx 2 and dV 3 /dx 3 are defined accordingly, and the non-positive nature of the terms of the dV/dt expression is determined.
The procedure comprises the following tasks: (i) determination of the expression for dV 2 /dt + kdV 3 /dt; (ii) arrangement of the expression for dV 2 /dt + kdV 3 /dt in terms of f v1 ; (iii) definition of gradient dV 2 /dx 2 and dead zone quadratic form V 2 ; (iv) arrangement of the expression for dV 2 /dt + kdV 3 /dt in terms of a non-positive function of dV 2 /dx 2 ; (v) definition of gradient dV 3 /dx 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 dV/dt = dV 1 /dt + kdV 2 /dt + kdV 3 /dt; (vii) integration of the dV/dt expression. The definition of the gradients dV 2 /dx 2 and dV 3 /dx 3 and the quadratic forms V 2 and V 3 requires the determination of the expression for dV 2 /dt + kdV 3 /dt and the properties of the terms involved. where: and S eq 2 is provided by Equations (5)-(7). (Tii) the state x 3 converges asymptotically to Ω x3 , 3 and S eq 3 is provided by Equations (5)-(7). (Tiii) The sets Ω x = Ω x1 Ω x2 Ω x3 and Ω S = Ω S1 Ω S2 Ω S3 are invariant.
The dx 2 /dt expression (13) can be rewritten as: where: The time derivative of V 2 can be expressed as: Substituting the dx 2 /dt expression (Equation (48)) and arranging, yields: The time derivative of V 3 can be expressed as: Substituting the dx 3 /dt expression (Equation (15)) and arranging, yields Adding dV 2 /dt (52) and dV 3 /dt (56), yields Task 2. Arrangement of the expression for dV 2 /dt + kdV 3 /dt in terms of f v1 . The x 1 signal appearing in Equation (57) can be expressed as the sum of f v1 and a disturbance term. Let, Using the definition of f v1 (29), we have: From Equation (58), x 1 can be expressed as: Substituting in expression (57) and arranging, yields: where The effect of the term f v2 (Q/v 2 ) f v1 is tackled later by considering the dV 1 /dt expression.

Task 3.
Definition of gradient f v2 = dV 2 /dx 2 and dead zone quadratic form V 2 .
At what follows, we examine the properties of δ 2 and the term g x2 − δ 2 appearing in Equation (62), what allows us to choose f v2 and V 2 . From the definition of δ 2 (63) and d x1 (59) it follows that δ 2 ∈ x l 1 , x u 1 . Hence, The gradient of g x2 (49) is: Furthermore, accounting for the definition of r x2 (14), the definition of g x2 in Equation (49), and the properties of r x2 stated in assumption 3, we have: is bounded f or x 2 bounded, and sgn(g x2 ) = sgn(x 2 ).
In turn, these properties lead to: where The term δ 2 is defined in Equations (63) and (59) and satisfies Equations (64) and (65). From Equations (67)-(69) it follows that: To obtain the non-positive nature of the term (−1) f v2 (g x2 − δ 2 ), appearing in Equation (62), we need to choose f v2 such that: In view of properties (72) and to fulfill (73) and (74), we choose: The main properties of f v2 are: (Pi) f v2 is continuous with respect to x 2 , A Lyapunov function that satisfies dV 2 /dx 2 = f v2 Equation (51) is: The main properties of V 2 are: V 2 is continuous with respect to x 2 .
Task 4. Arrangement of the expression for dV 2 /dt + kdV 3 /dt in terms of a non-positive function of f v2 = dV 2 /dx 2 .
From properties Pv (80) and Pii (77), it follows that the term f v2 (g x2 − δ 2 ) appearing in Equation (62) satisfies: In addition, the term (g x2 − δ 2 ) satisfies the following: where, (87) The main properties of g 2v are: (Pi) g 2v is continuous with respect to x 2 , (88) From Equations (85) and (86), it follows that: From the definitions of g 2v (87) and f v2 (75) it follows that: Therefore, Equation (92) yields: Using this, Equation (62), yields: Task 5. Definition of gradient dV 3 /dx 3 and dead zone quadratic form V 3 . We choose f v3 with the structure of f v2 : where the bounds x u 3 , x l 3 are defined as: and x l 2 , x u 2 are defined in Equations (70) and (71). A Lyapunov function that satisfies the gradient dV 3 /dx 3 = f v3 (Equation (54)) is: The main properties of V 3 are: V 3 is continuous with respect to x 3 .
As a consequence of the definitions of f v2 (75) and f v3 (96), the property: holds true, as proved in Appendix A. Equation (95) can be arranged as: where β 1 , β 2 are positive constants that satisfy: Task 6. Definition of the overall Lyapunov function V and determination of the nonpositive nature of the expression for dV/dt = dV 1 /dt + kdV 2 /dt + dV 3 /dt.
Q v 2 f v1 appearing in Equation (103) can be expressed in terms of f 2 v1 : Substituting into Equation (103), we have: We tackle the effect of the f 2 v1 term by incorporating the dV 1 /dt expression. Equation (44) leads to: d dt Combining with Equation (106) and accounting for the property (102), yields: which can be rewritten as: where V is the overall Lyapunov function, defined as: Task 7. Integration of the dV/dt expression. Integrating Equation (108), yields: Therefore: (BPi) V 2 ∈ L ∞ , V 3 ∈ L ∞ , V 1 ∈ L ∞ , and from the definitions of V 1 (35), Considering the properties f 2 v2 ∈ L 1 and f 2 v2 ∈ L ∞ and applying the Barbalat's lemma (cf. [36]), yields: lim Furthermore, accounting for the definition of f v2 (75), it follows that x 2 converges asymptotically to Ω x2 = x 2 : x l 2 ≤ x 2 ≤ x u 2 . From this convergence result and the definition of x 2 (9), it follows that S 2 converges asymptotically to Ω S2 ; where: S eq 2 is provided by Equations (5)- (7); the bounds x l 2 , x u 2 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 ( Hence, either x 2 − x 3 or f v2 − f v3 converges to zero. This result and the convergence of x 2 to Ω x2 , the convergence of f v2 to zero, and the definition of f v3 (96) and f v2 (75), imply: (i) if x 2 − x 3 converges to zero, then x 3 converges to Ω x3 , Ω x3 = x 3 : x l 3 ≤ x 3 ≤ x u 3 ; (ii) if f v2 − f v3 converges to zero, then f v3 converges to zero, which implies the convergence of x 3 to Ω x3 . Therefore, as is indicated by both of these possibilities, x 3 converges asymptotically to Ω x3 . From this convergence result and the definition of x 3 (10) it follows that S 3 converges asymptotically to Ω S3 , where: S eq 3 is provided by Equations (5)- (7), and the bounds x l 3 , x u 3 are defined in Equation (97). This completes the proof of the second part of the theorem (Tii).

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).  | ≥ | f 2v | (93) which was required for proving the convergence of x 2 . These tasks are crucial for proving the non-positive nature of dV/dt. 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 S 1 , S 2 , S 3 are inside the convergence region Ω s , that is, S 1 ∈ Ω s1 , S 2 ∈ Ω s2 and S 3 ∈ Ω s3 simultaneously, they remain inside afterwards. (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.

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: The terms r s1 and r s2 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, BOD 5 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.
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 l 1 , S 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 convergences to the region within computed bounds , (see Figures 2a,3a,4a,4b); the state converges to the region within computed bounds , (see Figures 2b,3b,4c,4d); and converges to the region within the computed bounds , (see Figures 2b,3b,4c,4d).     Also, the three simulation cases illustrate the results stated in Theorems 1 and 2: (i) the invariant nature of Ω s1 , stated in Theorem 1 and discussed in Remark 2: one can see that once S 1 enters Ω s1 , it remains inside; (ii) the equivalence of the convergence regions Ω s2 and Ω s3 defined in Theorem 2: one can see that S 2 and S 3 converge to the same regions; (iii) the definitions of Ω s1 and Ω s2 , which indicate that Ω s1 and Ω s2 can be quite different, as Ω s2 depends on the reaction rate r x2 : one can see this remarkable difference in the regions to which S 1 and S 2 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 ∈ Ω s1 , S 2 ∈ Ω s2 , S 3 ∈ Ω s3 , simultaneously, they remain inside afterwards.

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: