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 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.
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: 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 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
is the pollutant concentration in CSTR 1;
is the pollutant concentration in CSTR 2;
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;
is the reaction rate of CSTR 1;
is the reaction rate of CSTR 2. In addition,
is defined in
is defined in
is defined in
.
Assumptions:
Assumption 1. ,,,, andare positive and constant. This implies that,,andare constant and positive.
Assumption 2. is positive and constant,is varying but bounded,.
Assumption 3. The reaction rate rs1 only depends on, it is bounded forbounded,;only depends on, it is bounded forbounded, and.
In the case of constant
,
= 0 in Equation (4),
converges to
,
converges to
and
converges to
, where
,
and
are positive constants obtained by equating Equations (1)–(3) to zero:
The
dynamics is obtained by subtracting Equation (5) from Equation (1):
where
,
is
evaluated at
.
The
dynamics is obtained by subtracting Equation (6) from Equation (2):
The term is evaluated at .
The
dynamics is obtained by subtracting Equation (7) from Equation (3):
In addition,
is defined in [−
is defined in [−
and
is defined in [−
. 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 is defined for the first state , for the second state and for the third state , and their properties are determined. An overall Lyapunov function V is defined as a weighted sum of the quadratic forms , , and . 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 , where is the dead-zone quadratic form for . The function and its gradient are defined so as to achieve this goal.
The procedure comprises the following tasks: (i) determination of the expression for ; (ii) definition of the gradient and the dead zone quadratic form for the state ; (iii) arrangement of the expression in terms of a non-positive function of ; (iv) integration of the expression. The definition of the gradient . and requires the determination of the 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 converges asymptotically to , , where , , , and , and converges asymptotically to where , , and is provided by Equation (5). (Tii) The sets and are invariant.
Proof. Task 1. Determination of the expression.
The time derivative of
can be expressed as:
Incorporating the
dynamics (11) and arranging, yields:
where,
□
Task 2. Definition of the gradient and the dead zone quadratic form for the state .
At what follows, we examine the properties of the term
appearing in Equation (18). From Equations (12), (19) and (20) and Assumption 3 we have:
where,
To obtain non-positive
in Equation (18), we need to choose
such that:
In view of properties (26) and to fulfill conditions (27) and (28), we choose:
The main properties of
are:
A dead-zone Lyapunov function that satisfies
(Equation (17)) is:
The main properties of
are:
□
Task 3. Arrangement of the expression in terms of a non-positive function of .
From properties (31), (34), it follows that the term
appearing in Equation (18) satisfies:
From the definition of
(19) and definitions of
(24), (25), it follows that:
From Equations (38) and (39), it follows that:
From the property (21), definition of
(40), and definition of
(29), it follows that:
Therefore, Equation (41) yields:
Using this in Equation (18) yields:
Accounting for the definition of
(35), we have:
□
Task 4. Integration of the expression.
From the above expression, it follows that:
That is,
(35) converges exponentially to zero, and consequently
converges exponentially to zero. Further, accounting for the definition of
(29), it follows that
converges asymptotically to
. From this convergence result and the definition of
(8), it follows that
converges asymptotically to
; where,
is provided by Equation (5), the bounds , are defined in Equations (24) and (25). This completes the proof of Ti.
From Equation (44) and properties (31) and (32), it follows that and . Therefore, the set is invariant, and consequently, is invariant. This completes the proof of Tii. □
Remark 1. The property (42) is crucial for proving the asymptotic convergence of x1.
Remark 2. The invariance of the set stated in Theorem 1 implies that once the state is inside the convergence region , that is, , 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 and to their compact sets, is based on determining the non-positive nature of , where , is the overall Lyapunov function, is the dead zone quadratic form for the state and the dead zone quadratic form for the state , and and are positive constants. To achieve this goal, the functions and and the gradients and are defined accordingly, and the non-positive nature of the terms of the expression is determined.
The procedure comprises the following tasks: (i) determination of the expression for ; (ii) arrangement of the expression for in terms of ; (iii) definition of gradient and dead zone quadratic form ; (iv) arrangement of the expression for in terms of a non-positive function of ; (v) definition of gradient and dead zone quadratic form ; (vi) definition of the overall Lyapunov function V and determination of the non-positive nature of the expression for ; (vii) integration of the expression. The definition of the gradients and and the quadratic forms and requires the determination of the expression for and the properties of the terms involved.
Theorem 2. Consider the model (11)–(15), subject to assumptions 1 to 3. (Ti) the state converges asymptotically to , , where:converges asymptotically to
where
,
,
and
is provided by Equations (5)–(7). (Tii) the state x3 converges asymptotically to ,
, where
, and
converges asymptotically to
where
,
,
and
is provided by Equations (5)–(7). (Tiii) The sets and
are invariant. Proof. Task 1. Determination of the expression for .
The
expression (13) can be rewritten as:
where:
The time derivative of
can be expressed as:
Substituting the
expression (Equation (48)) and arranging, yields:
The time derivative of
can be expressed as:
Substituting the
expression (Equation (15)) and arranging, yields
Adding
(52) and
(56), yields
□
Task 2. Arrangement of the expression for in terms of .
The
signal appearing in Equation (57) can be expressed as the sum of
and a disturbance term. Let,
Using the definition of
(29), we have:
From Equation (58),
can be expressed as:
Substituting in expression (57) and arranging, yields:
where
The effect of the term is tackled later by considering the expression. □
Task 3. Definition of gradient and dead zone quadratic form .
At what follows, we examine the properties of and the term appearing in Equation (62), what allows us to choose and . From the definition of (63) and (59) it follows that .
The gradient of
(49) is:
Furthermore, accounting for the definition of
(14), the definition of
in Equation (49), and the properties of
stated in assumption 3, we have:
In turn, these properties lead to:
where
The term
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
, appearing in Equation (62), we need to choose
such that:
In view of properties (72) and to fulfill (73) and (74), we choose:
The main properties of
are:
A Lyapunov function that satisfies
Equation (51) is:
The main properties of
are:
□
Task 4. Arrangement of the expression for in terms of a non-positive function of .
From properties Pv (80) and Pii (77), it follows that the term
appearing in Equation (62) satisfies:
In addition, the term
satisfies the following:
where,
The main properties of
are:
From Equations (85) and (86), it follows that:
From the definitions of
(87) and
(75) it follows that:
Therefore, Equation (92) yields:
Using this, Equation (62), yields:
□
Task 5. Definition of gradient and dead zone quadratic form .
We choose
with the structure of
:
where the bounds
,
are defined as:
and
,
are defined in Equations (70) and (71). A Lyapunov function that satisfies the gradient
(Equation (54)) is:
The main properties of
are:
As a consequence of the definitions of
(75) and
(96), the property:
holds true, as proved in
Appendix A. Equation (95) can be arranged as:
where
,
are positive constants that satisfy:
□
Task 6. Definition of the overall Lyapunov function and determination of the non-positive nature of the expression for .
The term
appearing in Equation (103) can be expressed in terms of
:
Substituting into Equation (103), we have:
We tackle the effect of the
term by incorporating the
expression. Equation (44) leads to:
Combining with Equation (106) and accounting for the property (102), yields:
which can be rewritten as:
where
is the overall Lyapunov function, defined as:
□
Task 7. Integration of the expression.
Integrating Equation (108), yields:
Therefore: (BPi)
,
,
, and from the definitions of
(35),
(81),
(98) it follows that
,
,
, and consequently
; (BPii)
,
,
, what follows from the definitions of
(29),
(75),
(96); (BPiii)
. Considering the properties
and
and applying the Barbalat’s lemma (cf. [
36]), yields:
Furthermore, accounting for the definition of
(75), it follows that
converges asymptotically to
. From this convergence result and the definition of
x2 (9), it follows that
S2 converges asymptotically to
; where:
is provided by Equations (5)–(7); the bounds , 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
and from properties Bpi, BPii it follows that
. Applying the Lasalle’s theorem (cf. [
22]) yields
.
Hence, either
or
converges to zero. This result and the convergence of
to
, the convergence of
to zero, and the definition of
(96) and
(75), imply: (i) if
converges to zero, then
converges to
,
; (ii) if
converges to zero, then
converges to zero, which implies the convergence of
to
. Therefore, as is indicated by both of these possibilities,
converges asymptotically to
From this convergence result and the definition of
(10) it follows that
converges asymptotically to
, where:
is provided by Equations (5)–(7), and the bounds , are defined in Equation (97). This completes the proof of the second part of the theorem (Tii).
From Equation (109), it follows that:
The condition implies: (i) , and , as follows from definitions of V (110), V1 (35), V2 (81), and V3 (98); (ii) , , , as follows from the definition of (29), (75) and (96). Therefore, it follows from (117) that is invariant, and consequently 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 , which is arbitrarily large but bounded;
the initial values of ,
, and
can take arbitrarily large but bounded positive values, since
is defined in [− is defined in [− and
is defined in [−, as stated in Section 2.
Remark 2. In the study of the stability of , the dynamics of and (2), (3) are not accounted for, as can be noticed in the proof of Theorem 1. This is a consequence of the fact that and are not involved in the dynamics of (11). In contrast, in the study of the stability of and , the dynamics of , and are accounted for, as can be noticed in the proof of theorem 2. This is a consequence of the fact that , and are involved in the dynamics of (13).
Remark 3. The first tank influences the transient value of , which can be noticed from Equation (111). Also, the first tank influences the width of the convergence region of , that is [, ], what can be noticed from the definition of , (70), (71), the definition of δ2 (63), (58), the properties of (64), (65), and the definition of , (24), (25).
Remark 4. The presence of the third tank influences the transient value of , what can be noticed from Equation (111). However, it does not influence the width of the convergence region of the state , that is, [, ], since the bounds (, ) are not influenced by the parameters of the third tank, that is, and , which can be noticed from the definition of , (70), (71).
Remark 5. Important tasks of the stability analysis are: the definition of the overall Lyapunov-like function (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 ; the property (93) which was required for proving the convergence of . These tasks are crucial for proving the non-positive nature of . 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 stated in Theorem 2 implies that once the three states S1, S2, S3 are inside the convergence region , that is, , and 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:
The terms
and
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.
To evaluate the obtained convergence results stated in Theorems 1 and 2, model (1)–(4) is simulated with
The bounds
,
of
(47); the bounds
,
of
(114); and the bounds
,
of
(116) are computed and represented by the upper and lower horizontal dotted lines, whereas the equilibrium values
,
and
(Equations (5)–(7)) are represented by the middle horizontal dotted lines. The calculation of the aforementioned bounds requires the calculation of bounds
,
(24), (25);
,
(70), (71); and
,
(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
Figure 2a,
Figure 3a,
Figure 4a,b); the state
converges to the region within computed bounds
,
(see
Figure 2b,
Figure 3b,
Figure 4c,d); and
converges to the region within the computed bounds
,
(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 , stated in Theorem 1 and discussed in Remark 2: one can see that once S1 enters , it remains inside; (ii) the equivalence of the convergence regions and defined in Theorem 2: one can see that S2 and S3 converge to the same regions; (iii) the definitions of and , which indicate that and can be quite different, as 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 , stated in Theorem 2 and discussed in Remark 6: one can see that once the three state variables are inside the convergence region , that is, , , , simultaneously, they remain inside afterwards.