**Condition 1.**


Therefore, given the nature of phenomenology, the independent variable throughout this work will be the position *L*. We have the following Streeter–Phelps mathematical model respect position [35]:

$$\begin{aligned} \frac{d\mathbf{x}\_1}{dL} &= -\frac{k\_1(L)}{lI}\mathbf{x}\_2 + \frac{k\_2(L)}{lI}(D\_s - \mathbf{x}\_1) \\ \frac{d\mathbf{x}\_2}{dL} &= -\frac{k\_1(L)}{lI}\mathbf{x}\_2. \end{aligned} \tag{1}$$

where the parameters *k*1(*L*) > 0 are the BOD removal rate coefficient, *k*2(*L*) > 0 is the reaeration coefficient, and *Ds* is the oxygen saturation level. This type of system tends to be extremely sensitive to parameter changes, in this case, parameters, such as linear velocity *U* and the variation of constants *k*1,2 [33]. Therefore, it is necessary to perform a parametric sensitivity analysis, which will be substantial to design a state observer. In this paper, we propose such analysis by using stability concepts in the sense of Lyapunov. For the above, it is necessary to make some conceptual assumptions to study stability analysis.

**Condition 2.** *It is proposed that the model parameters (1) fulfilled following condition:*

$$\begin{aligned} k\_1(L) &= k\_1 + \beta\_1(L) \\ k\_2(L) &= k\_2 + \beta\_2(L), \end{aligned} \tag{2}$$

*where k*1,2 *are nominal known parameters; β*1(*L*) *and β*2(*L*) *are continuously bounded functions* |*β*1(*L*)| < *<sup>α</sup>*1*l and* |*β*2(*L*)| < *<sup>α</sup>*2*L; and the <sup>α</sup>*1,2 > 0 *are Lipschitz constants respect to longitude L. In this paper, we propose that β*1(*L*), *β*2(*L*) *functions are unknown, unexpected but bounded, which cause unwanted parametric distortion <sup>δ</sup>*(*L*)*.*

Substituting (2) on the mathematical model (1) and reducing:

$$\begin{split} \frac{d\mathbf{x}\_1}{dL} &= -\frac{k\_1}{\underline{U}}\mathbf{x}\_2 - \frac{k\_2}{\underline{U}}\mathbf{x}\_1 + \delta\_1 + \frac{k\_2}{\underline{U}}D\_s\\ \frac{d\mathbf{x}\_2}{dL} &= -\frac{k\_1}{\underline{U}}\mathbf{x}\_2 + \delta\_2\\ \delta\_1 &= -\frac{1}{\underline{U}}(-\beta\_1(L)\mathbf{x}\_1 - \beta\_2(l)(D\_s - \mathbf{x}\_2))\\ \delta\_2 &= -\frac{1}{\underline{U}}\beta\_1(L)\mathbf{x}\_1. \end{split} \tag{3}$$

The previous system (3) can be represented in its state form as a autonomous linear system in bounded disturbance presence:

$$\frac{d\mathbf{x}}{dL} = A\mathbf{x} + \delta(L) \tag{4}$$

for *x<sup>T</sup>* = [*<sup>x</sup>*1, *<sup>x</sup>*2]*<sup>T</sup>* , *δ<sup>T</sup>*(*L*) = 1*U* [*δ*1 + *k*2*Ds*, *<sup>δ</sup>*2]*<sup>T</sup>*, and *A* = −*k*2 1*U* −*k*<sup>1</sup> 1*U* 0 −*k*<sup>2</sup> 1*U* . Hence, a general condition to the global disturbance *δ*(*L*) is proposed.

**Condition 3.** *Let δ*(*L*) *be a bounded unknown nonlinear function such that the following is fulfilled:*

$$\|\|\delta(L)\|\| \le D\_{\text{max}} > 0. \tag{5}$$

**Note 1**. *In this paper, we will assume that the linear velocity remains constant, as we focus on the possible change in kinetic, chemical parameters (k*1,2*). The parametric distortions due to linear velocity (Uriver) changes will not be treated in this work due to an extension of it. But, it is easy to see that this kind of distortions (linear velocity changes) would cause limited delta distortions; therefore, for Uriver* = *Unom* + *δ*(*U*)*U* :

$$\frac{d\mathbf{x}}{dL} = A\_{U\_{\text{new}}}\mathbf{x} + \delta(\mathbb{U})X \tag{6}$$

$$\left||\delta(\mathbb{U})\right|| \le \left||A\_{U\_{\text{new}}}\right||.$$

Therefore, the following theoretical result is proposed, where parametric sensitivity is analyzed via a stability analysis using Lyapunov's convergence functions.

**Theorem 1.** *Let x* = 0 *be the spatial equilibrium point of (4) and δ*(*L*) *fulfilled above Condition 3*. *The origin x* = 0 *will maintain a bonded dynamic in a ball of convergence in parametric chances presence for non-negative matrix Q and P such that: <sup>λ</sup>*min{*Q*} >> <sup>2</sup>*D*max*λ*max{*P*}.

**Proof.** The following Lyapunov function is proposed: *V* = 1*u xTPx* with independent variable *L*, if it replaces the derivatives under trajectories systems:

$$\begin{split} \frac{dV}{dL} &= \frac{1}{U} \mathbf{x}^T P \mathbf{x} + \frac{1}{U} \mathbf{x}^T P \dot{\mathbf{x}} \\ \frac{dV}{dL} &= \mathbf{x}^T (PA + A^T P) \mathbf{x} + 2\mathbf{x}^T P \delta(t). \end{split} \tag{7}$$

Since *A* of (4) is Hurwitz matrix for *U*, *k*1, *k*2 > 0, then it is Re *λi* > 0 (due to the nature non-negative parametric system). Hence, Lyapunov algebraic equation calculation *PA* + *ATP* = −*Q* is fulfilled.

Therefore, (7) is reduced, with *c* = <sup>2</sup>*D*max*λ*max{*P*}, *<sup>λ</sup>*min{*Q*} = *a* and bounded:

$$\begin{aligned} \frac{dV}{dL} &= -\mathbf{x}^T Q \mathbf{x} + 2\mathbf{x}^T P \delta(t) \\ \frac{dV}{dL} &\le -a||\mathbf{x}||^2 + c||\mathbf{x}|| \\ \frac{dV}{dL} &\le -a||\mathbf{x}|| \left( ||\mathbf{x}|| - \frac{c}{a} \right) .\end{aligned} \tag{8}$$

The higher the value of the system, the faster it will converge in the neighborhood to a ball of convergence such that

$$B\_{\mathbf{x}} = \left\{ \|\mathbf{x}\| \in \mathbb{R}^2 : \|\mathbf{x}\| < \frac{\mathcal{C}}{a} \right\} \tag{9}$$

$$B\_{\mathbf{x}} = \left\{ \|\mathbf{x}\| \in \mathbb{R}^2 : \|\mathbf{x}\| < \frac{2D\_{\text{max}}\lambda\_{\text{max}}\{P\}}{\lambda\_{\text{min}}\{Q\}} \right\}.$$

In the disturbances system presence, (4) in spatial equilibrium point will remain Ultimated Bounded [36].

In order to exemplify the previous theorem use, as well as *P*, *Q* matrix calculation, we will use real data obtained from literature and be able to make a parametric stability analysis and sensitivity in the theorem above sense. *U* = 1 is proposed, and it is assumed that the parameters of *k*1, *k*2 have a 5% maximum change from the nominal value.

First, via Condition 2, the maximum bounded is obtained |*β*1(*l*)| = 0.05*k*1,*nom* and |*β*2(*l*)| = 0.05*k*2,*nom*. Thus, by replacing this in (4) with nominal parameters, *k*1 = 0.3 day−<sup>1</sup> and *k*2 = 0.06 day−<sup>1</sup> of literature [34]. Disturbance vector *δ*(*L*) will be calculated:

$$\delta(L) = \begin{pmatrix} \beta\_1(l)\mathbf{x}\_1\\ \beta\_1(l)\mathbf{x}\_1 - \beta\_2(t)\mathbf{x}\_2 \end{pmatrix}$$

$$\delta(L) = \begin{pmatrix} -0.03\mathbf{x}\_1\\ 0.03\mathbf{x}\_1 - 0.006\mathbf{x}\_2 \end{pmatrix}.\tag{10}$$

If parameters are substituted, in matrix A (4), it is obtained:

$$A = \begin{bmatrix} -0.3 & 0\\ -0.3 & -0.06 \end{bmatrix}.$$

Therefore, there is a Λ = 0.05*A* matrix, such that *δ*(*L*) = Λ*x*, such that:

$$\delta(L) = \underbrace{\begin{bmatrix} -0.015 & 0\\ 0.015 & -0.003 \end{bmatrix}}\_{\Lambda.} \text{s.}$$

If the matrix *Q* = −*I* is chosen, Lyapunov *P* matrix calculation is obtained by matrix system solution *PA* + *ATP* = <sup>−</sup>*I*, such that:

$$P = \begin{bmatrix} 1.6667 & -1.3889 \\ -1.3889 & 15.2778 \end{bmatrix} \text{.} \tag{11}$$

where the eigenvalues of *P* are Re *λi* = [1.667 90.2778], to obtain the convergence ball of (9), in the presence of a constant perturbation of 5%, such that *<sup>λ</sup>*max{*P*} = 15.2778, *<sup>λ</sup>*max{*Q*} = 1 and *<sup>λ</sup>*max{Λ} = 0.03. In addition:

$$\begin{aligned} \|\delta(L)\| &= \|\Lambda \mathbf{x}\|\\ \|\delta(L)\| &\le 0.05 \|A\| \|\mathbf{x}\| \le 0.05 \lambda\_{\text{max}} \{\Lambda\} \left\|\mathbf{x}\right\| 0.015 \left\|\mathbf{x}\right\|\\ D\_{\text{max}} &\le 0.015 \left\|\mathbf{x}\right\|. \end{aligned} \tag{12}$$

Substituting (12) on the inequality calculation (8):

$$\begin{aligned} \frac{dV}{dL} &\leq -a\left\|\mathbf{x}\right\|^2 + c\left\|\mathbf{x}\right\| \leq -\lambda\_{\text{min}}\{Q\}\left\|\mathbf{x}\right\|^2 + 2D\_{\text{max}}\lambda\_{\text{max}}\{P\}\left\|\mathbf{x}\right\| \\\frac{dV}{dL} &\leq -\left\|\mathbf{x}\right\|^2 + 0.4581\left\|\mathbf{x}\right\|^2 \leq -0.5419\left\|\mathbf{x}\right\|^2. \end{aligned}$$

Since Lyapunov's function *V* = 1*u xTAx* and Raleigh's inequality *U <sup>λ</sup>*max{*A*}*<sup>V</sup>* ≤ *x*2, it is had:

$$\frac{dV}{dL} \le -1.80V.$$

Therefore, finally:

$$V = V(0)e^{-1.80l}$$

For a constant perturbation, the system reaches *asymptotic stability*, which is a sub-case of *practical stability*. When perturbation varies from the position, the system will converge to a ball of convergence in the sense of Reference [36].

## *Problem Statement*

The system represented in the linear system (4) has analytical properties that do not allow state estimators to be used, since the *y* output is a signal that does not permit all states reconstruction, due to general observability matrix, has not been full range, and this has already been studied in the literature. Through certain analytical assumptions about the system (4), a non-linear observable system is obtained, with which it is possible to estimate variables via the use of an observer [32]. By assuming a linear relationship between the BOD removal rate and DO, with zero DO giving no decrease in BOD, a new removal rate coefficient is defined by *k*3 = *k*1*x*1. These modifications, performed in (1), show:

$$\begin{aligned} \, \_\text{II} \frac{d\mathbf{x}\_1}{dL} &= -k\_1 \mathbf{x}\_2 + k\_2 (D\_s - \mathbf{x}\_1) + I\_D - c\_1 \mathbf{x}\_1 \\ \, \_\text{II} \frac{d\mathbf{x}\_2}{dL} &= -k\_3 \mathbf{x}\_1 \mathbf{x}\_2 + I\_B - c\_1 \mathbf{x}\_2 \end{aligned} \tag{13}$$

where *ID* is the inflow of DO, *IB* the BOD inflow, *c*1 the oxygen turnover rate due to flow through, and, on the average, *ID* = *c*1*x*1 and *IB* = *c*1*x*2. Thus, the previous system can be rewritten as follows:

$$\begin{aligned} \frac{dX}{dL} &= A\_0 X + f(X) \\ y &= x\_{1\prime} \end{aligned} \tag{14}$$

where *A*0 = 0 1 0 0, and *f*(*X*) = −*x*2 + 1*U* (−*k*1*x*<sup>2</sup> + *k*2(*Ds* − *<sup>x</sup>*1) + *ID* − *<sup>c</sup>*1*x*1) 1*U* (−*k*3*x*1*x*<sup>2</sup> + *IB* − *<sup>c</sup>*1*x*2) .

For robust BOD estimation via the above model, we propose a novel fractional observer, and its analytical design is proposed in the following section.

#### **3. Design of a Novel Fractional Observer (FOHO)**

Firstly, fractional calculus is the expansion of traditional analysis to derivation and integration operations employing non-integer orders. In this paper, a fractional-order high gain observer (FOHO) is proposed in order to estimate non-measurable variables (BOD) of a river basin environmental monitoring problem. The Caputo fractional derivative of order *β* of a function *φ*(*L*) on the positive real half axis is determined. It is important to notice that the notation *<sup>D</sup>*−*βφ*(*L*) is used, as well as to denote the fractional integral of order *β* of function *φ*(*L*), more precisely *<sup>D</sup>*−*βφ*(*L*) ≡ *Iβφ*(*L*) [37]. The definitions of fractional derivative and fractional integral, as stated above, cannot be used in practice; thus, numeric methods, such as the one based on the Grünwald–Letnikov approach, are commonly used. Here, we deal with the same note in terms of the definition of the fragmentary calculation tool as in Reference [38]. Therefore:

$$M^{\beta}\phi(L) = \frac{1}{\Gamma(\beta)} \int\_0^L \frac{\phi(\tau)}{(L-\tau)^{1-\beta}} d\tau,\tag{15}$$

where 0 < *τ* < *L* and *β*. Thus, the following fractional observer (FOHO) is proposed for system (14):

$$\begin{array}{ll}\frac{dz}{dL} = & A\_o z + f(z) + bu + \Lambda(\theta) \mathbb{C}^T \mathbb{C} e + \delta\\ \delta'(L) = & -S\Lambda\_\theta^{-1} I^\theta \text{sign}(e) \end{array} . \tag{16}$$

For *e* = *Z* − *X*, it is as follows:


where *h*1, ..., *hn* > 0 are gains. *hi* is an observer gain, *θ* is a positive high gain, and Λ(*θ*) matrix is as follows: 

$$
\Lambda(\theta) = \begin{bmatrix}
\theta h\_1 & 1 & 0 & \cdots & 0 \\
\theta^2 h\_2 & 0 & 1 & \ddots & \cdots \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
\theta^{n-1} h\_{n-1} & 0 & \cdots & 0 & 1 \\
\theta^n h\_n & 0 & \cdots & 0 & 0
\end{bmatrix}.
$$

Hence, there is a matrix *Ah* such that it is fulfilled:

$$
\Delta\_{\theta}^{-1} \Lambda(\theta) \Delta\_{\theta} = \theta A\_{\text{h}}.\tag{17}
$$

There is a symmetric Δ*θ* such that Δ*θ* = *diag*(1, *θ*, *θ*2 ... *<sup>θ</sup><sup>n</sup>*−<sup>1</sup>). Thus, *Ah* is a Hurwitz matrix by construction, such that:

$$A\_{\hbar} = \begin{pmatrix} h\_1 & 1 & 0 & \cdots & 0 \\ h\_2 & 0 & 1 & \cdots & \vdots \\ \vdots & \vdots & \ddots & \ddots & 0 \\ h\_{\hbar-1} & 0 & 0 & 0 & 1 \\ h\_{\hbar} & 0 & 0 & 0 & 0 \end{pmatrix} . \tag{18}$$

**Condition 4.** *Let e* = *Z* − *X be estimated error such that the following is fulfilled: 1.AhwillbeaHurwitzmatrixalways existsamatrixS* = *S<sup>T</sup> suchthat[39]:*

$$S A\_h + A\_h^T S = -I.\tag{19}$$

*Therefore, for any θ* > 1*, introduce the next symmetric matrix:*

$$\mathbb{S}\_{\theta} = \Delta\_{\theta}^{-1} \mathbb{S} \Delta\_{\theta}^{-1}.$$

*2. Consider the error function e* = *z* − *x of (16) and (14). Existing <sup>λ</sup>*min{*S*} = *c* > 0 *such that [40]:*

$$\left| c \left| I^{\mathbb{B}} \text{sign}(e) \right| \right| \; > \delta\_{\text{max}}.\tag{20}$$

Therefore, the main general theoretical result of this work is presented by means of the following theorem.

**Theorem 2.** *Convergence of FOHO estimator. Let an observer (16) and (14) be, it is said, that general error of state estimation e*(*t*) = *z*(*L*) − *x*(*L*) *converge asymptomatically, with high gain function θ* > 1 *in unknown nonlinear perturbation presence* |*δ*(*t*)| < *D*max *if Ah is a Hurwitz matrix with D*max *known, and matrix C* = [100...0] *if it is fulfilled Condition 4.*

**Proof.** The estimate error *de dL* = *dz dL* − *dxdL* has been substituted in the dynamics Equations (16) and (14) corresponding:

$$\begin{aligned} \frac{de}{dL} &= A\_0 z + bu + \hat{\delta} + \Lambda(\theta) \mathbf{C}^T \mathbf{C} e - A\_0 \mathbf{x} - bu - \delta + k(\cdot) \\ \frac{de}{dL} &= (A\_0 + \Lambda(\theta) \mathbf{C}^T \mathbf{C}) e + \bar{\delta} + k(\cdot) \\ \frac{de}{dL} &= (A\_0 + \Lambda(\theta) \mathbf{C}^T \mathbf{C}) e + \bar{\delta} + k(\cdot), \end{aligned} \tag{21}$$

where *k*(·) = *f*(*Z*) − *f*(*X*) via Lipschitz relationship, and the following condition can be obtained [39], such that: *k*(·) ≤ *γ*(*z* − *x*) ≤ *γ<sup>e</sup>*. In addition, it has the following remarkable matrix properties:

$$
\Lambda(\boldsymbol{\theta})\boldsymbol{\mathsf{C}}^{\mathrm{T}}\mathbb{C} = \begin{bmatrix}
\theta h\_{1} & \boldsymbol{1} & \boldsymbol{0} & \cdots & \boldsymbol{0} \\
\theta^{2}h\_{2} & \boldsymbol{0} & \boldsymbol{1} & \ddots & \cdots \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
\theta^{n-1}h\_{n-1} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{1} \\
\theta^{n}h\_{n} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{1} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0} \\
\boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \ddots & \cdots \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
\boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0} \\
\boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0}
\end{bmatrix} \\
= \begin{bmatrix}
\theta h\_{1} & \boldsymbol{0} & \boldsymbol{0} & \cdots & \boldsymbol{0} \\
\boldsymbol{\theta}^{2}h\_{2} & \boldsymbol{0} & \boldsymbol{0} & \ddots & \cdots \\
\vdots & \vdots & \ddots & \ddots & \vdots \\
\boldsymbol{\theta}^{n-1}h\_{n-1} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0} \\
\boldsymbol{\theta}^{n}h\_{n} & \boldsymbol{0} & \cdots & \boldsymbol{0} & \boldsymbol{0}
\end{bmatrix}.
$$

Therefore, *A*0 + Λ(*θ*)*CTC* = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎣ *θh*1 1 0 ··· 0 *<sup>θ</sup>*2*h*2 0 1 ... ··· . . . . . . . . . . . . . . . *<sup>θ</sup><sup>n</sup>*−<sup>1</sup>*hn*−<sup>1</sup> 0 ··· 0 1 *θnhn* 0 ··· 0 0 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎦ = <sup>Λ</sup>(*θ*).

Hence, (21) is modified:

$$\frac{de}{dL} = \Lambda(\theta)\varepsilon + \tilde{\delta} + k(\cdot). \tag{22}$$

.

It is proposed a Lyapunov function *Vθ* (*e*) = *Vθ* with its derivates along trajectories of (22):

$$\begin{split} V &= \boldsymbol{\epsilon}^{T} \mathbf{S}\_{\theta} \boldsymbol{\epsilon} \\ \frac{dV}{dL} &= \boldsymbol{\epsilon}^{T} (\mathbf{S}\_{\theta} \boldsymbol{\Lambda}(\theta) + \boldsymbol{\Lambda}^{T}(\theta) \mathbf{S}\_{\theta}) \boldsymbol{\epsilon} + 2\boldsymbol{\epsilon}^{T} \mathbf{S}\_{\theta} \boldsymbol{\delta} + 2\boldsymbol{\epsilon}^{T} \mathbf{S}\_{\theta} \boldsymbol{k}(\cdot) . \end{split} \tag{23}$$

Recalling Condition 4 *Sθ* = Δ−<sup>1</sup> *θ S*Δ−<sup>1</sup> *θ* ,*eθ* = Δ−<sup>1</sup> *θ e*, and *eTθ* = (<sup>Δ</sup>−<sup>1</sup> *θ e*)*<sup>T</sup>* and reshuffling terms:

$$\frac{dV}{dL} = e\_{\theta}^{T} [S\Delta\_{\theta}^{-1}\Lambda(\theta)\Delta\_{\theta} + (\Delta\_{\theta}^{-1}\Lambda(\theta)\Delta\_{\theta})^{T}S]e\_{\theta} + 2e\_{\theta}^{T}Sk(\cdot) + 2e\_{\theta}^{T}S\Delta\_{\theta}^{-1}\delta.\tag{24}$$

Replacing (17) on (24):

$$\frac{dV}{dL} = \theta \epsilon\_\theta^T [SA\_h + A\_h^T S] \mathbf{e}\_\theta + 2 \epsilon\_\theta^T S \mathbf{k}(\cdot) - 2 \epsilon\_\theta^T S \Delta\_\theta^{-1} \tilde{\delta}.\tag{25}$$

Giving ˜ *δ* = ˆ *δ* − *δ*, if we substitute the ˆ *δ* of observer (16):

$$\frac{dV}{dL} = \theta \varepsilon\_{\theta}^{T} [\mathbb{S}A\_{h} + A\_{h}^{T}\mathbb{S}]\varepsilon\_{\theta} + 2\varepsilon\_{\theta}^{T}\mathbb{S}k(\cdot) - 2\varepsilon\_{\theta}^{T}\mathbb{S}\Delta\_{\theta}^{-1}(2\mathbb{S}\Delta\_{\theta}^{-1}I^{\overline{\theta}}\text{sign}(\varepsilon) - \delta). \tag{26}$$

Given that *Ah* is a Hurwitz matrix and *δ* ≤ *D*max and fulfilling Condition 4, it has, as follows:

$$\begin{split} \frac{dV}{dL} &= -\theta e\_{\theta}^{T} e\_{\theta} + 2e\_{\theta}^{T} \text{Sk}(\cdot) - 2e\_{\theta}^{T} \text{S} \Delta\_{\theta}^{-1} (2\text{S} \Delta\_{\theta}^{-1} I^{\theta} \text{sign}(e) - \delta) \\ \frac{dV}{dL} &= -\theta e\_{\theta}^{T} e\_{\theta} + 2e\_{\theta}^{T} \text{Sk}(\cdot) - e\_{\theta}^{T} \text{S} \Delta\_{\theta}^{-1} (\text{S} \Delta\_{\theta}^{-1} I^{\theta} \text{sign}(e) - \delta) \\ \frac{dV}{dL} &\leq -(\theta - b) \left\| \left\| e \right\|\_{\theta}^{2} + c \left\| c \right\|\_{\theta} \left( c \left\| I^{\theta} \text{sign}(e) \right\| \right) - D\_{\text{max}} \right). \end{split} \tag{27}$$

Since *<sup>λ</sup>*max{<sup>Δ</sup>−<sup>1</sup> *θ* } = 1 for one high gain *θ b*, where *b* = <sup>2</sup>*λ*min{*S*}*<sup>λ</sup>*min, *c* = *<sup>λ</sup>*max{*S*}*<sup>λ</sup>*max{<sup>Δ</sup>−<sup>1</sup> *θ* } = *<sup>λ</sup>*min{*S*}, the observer will converge asymptomatically to zero when the longitude tends to infinity, in same sense as Reference [37,40].

#### **4. Materials and Methods**

In order to carry out the numerical simulation study, it is necessary to mention the conditions and system parameters to be studied, in this case, an urban river in Mexico. All these data were obtained from the literature; thus, it is possible to carry out a theoretical numerical simulation to be able to appreciate the *β* parameter observer effect on the BOD estimate in the river simulation. The Culiacán river, located in the State of Sinaloa, is of grea<sup>t</sup> regional importance since this country region is known as the granary of Mexico, and this river feeds the agricultural industry. Its course is 87.5 km long, and its basin covers 17,200 km2, with an annual flow of 3280 million hm3. It has an average depth of 1335 m, with an average width of 50 m and a flow of 0.5 m3/s [41,42].

The river is formed at the confluence Tamazula and Humaya rivers. It runs along the Pacific coastal plain, initially flowing through a large part of the urban area, which is why it has been indicated that the Culiacán River has a certain level of pollution resulting from the discharge of contaminated water from industrial processes [43]. Therefore, it is necessary to be able to monitor water quality and its DO and BOD levels. In this section, we sugges<sup>t</sup> a numerical simulation study where we suggested estimation BOD and monitoring of DO in presence to Streeter-Phelps model parameters changes. The modelling parameters are taken from those published in Reference [32], and we will have 100% *k*1 and *k*2 nominal parameters along the river position *L*, simulating the course of an aquatic boat robot at constant velocity in the environmental monitoring task by measuring DO (*x*1) and estimating BOD (*x*2) using the fractional observer shown (16).

In this paper, we performed a simulation study on the use of the Streeter–Phelps model for DO measurement and BOD estimation of the Culiacán river. The simulation was carried out in MATLAB-Simulink language, via a numerical method ODE45 with variable step. Based on the physical data model river parameters (14), the study was carried out using the simulation parameters summarised in the following Table 1.


**Table 1.** Model (14) and observer (16) simulation parameters.

To test FOHO robustness (16) and Luenberger classical [32] observer algorithms, we noted a parameter change throughout the simulation, recreating changes along river length (see to Figure 1). We sugges<sup>t</sup> changes on the Culiacán nominal Streeter–Phelps parameters from the position or length of the river. It is proposed that the DO signal (*x*1) recreates more realistic conditions; therefore, we add a white Gaussian noise signal. Since a theoretical numerical study is presented in this work, the *L* = 0 is arbitrary, and the initial conditions of BOD and DO are based on the state-of-the-art. Similarly, it is assumed that there are no additional pollution sources within the 2 km proposed in this numerical simulation.

**Figure 1.** Proposed changes in parameters.

## **5. Results**

In this brief chapter, we show the numerical simulation results carried out. We study the effect of the fractional observer order in presence of parameters changed. Simulating the spatial change dynamics in DO and BOD concentrations (see Figures 2 and 3), through that of the FOHO and the classical Luenberger observer [32], we can appreciate *β* parameter change of FOHO; this is unperceivable at a glance. Still, we can significantly understand the FOHO obtains better performance over its classical counterpart. The FOHO (16) rejects Streeter–Phelps parameter changes with significant performance.

**Figure 2.** Comparison of Dissolved oxygen (DO) monitoring estimate: the fractional observer (16) and the extended Luenberger classical observer.

**Figure 3.** Biochemical Oxygen Demand (BOD) monitoring using observer (16) and its changes depending on *β* fractional-order value and classic observer.

To quantify the performance and accuracy of the state observers, criteria, such as the Integral Distance Absolute error (IDAE), can be used to evaluate the estimation of states in processes of chemical or biochemical character usually used in the evaluation of observers [44–46]. To evaluate the effect of the *β* performance and accuracy of the state observers, we use two different types of error criteria, in same sense as Reference [44].

1. Steady-state performance factor (SSPF) ranges from 0 to 1, where zero indicates perfect process state reconstruction. Perfect reconstruction in this context means that the process variables in the estimated state vector coincide with the true process. The steady-state performance factor is defined as:

$$SSPF = \frac{|\mathbf{x}(t) - \mathbf{\hat{x}}(t)|}{\mathbf{x}(t)},\tag{28}$$

where *x* and *x*ˆ are evaluated at the steady-state point. Therefore, the SSPF determines the accuracy of a state observer at the steady state [44].

2. Integral Distance Absolute error (IDAE), which is inspired by the criteria shown in Reference [44]. It measures the observer dynamic performance and penalises the errors that persist for long position. In this work, the independent variable is not the time, but it is the length; therefore, the definition is the following:

$$IDAE = \int\_{0}^{L} L \left| \mathfrak{x}(L) - \mathfrak{x}(L) \right| dL. \tag{29}$$

To evaluate the effect of the *β* of (16) parameter on state estimation (see to Figure 4), we propose using IDAE and SSPF (see to Figure 5). We can appreciate that, in the presence of noise, a minor fractional-order redeems the value of the IDAE. In the case of classic mode absence, the IDAE is duplicated, proving the superiority of the observer proposed in this work over the classic one.

**Figure 5.** Observers Steady-state performance factor.

In the case of the SSPF, it appreciated a slight effect that has the fractional-order value (*β*), and the value of 0.1 is the one that maintains the best performance, while the classic Luenberger extended has the lowest performance. IDAE criterion importance using specific case is that the independent variable is the *L* length, and this criterion reflects the net yield from the river length to measure. Before reaching the steady-state, even the extended Luenberger observer [32] maintains a better performance. After *L* > 1, when the steady-state is reached, it loses effectiveness due to the parameter changes *k*1and *k*2 affecting the final value of the steady-state SSPF results, as in Figure 5.

The latest generation of fractional gain observers are superior to the Luenberger Extended versions, in the presence of perturbations at the kinetic constants *k*1 and *k*2. However, high-gain observers

depend mainly on the structure of the mathematical model. If there is a change in structure, it will result in a type of undefined perturbation, which could not be compensated by the estimator proposed in this paper. Thus, if maximum distortion is greater than the state's maximum, our algorithm would have problems because it is not a fulfilled primary theorem condition, which would be a problem or weakness of our algorithm. Therefore, it could be said that any disturbance that makes the system not comply with Conditions 2 and 3 will make the Fractional observer provided in this paper not perform correctly. When perturbations in parameters, such as *U* (volumetric flow rate), occur, as long as the condition shown is met, the fractional observer will maintain an excellent performance, since the condition of asymptotic stability of the main theorem is already fulfilled. Fortunately, the parametric and model distortions in most of the real practical cases meet this primary condition.

## **6. Conclusions**

In this work, a stability study for Streeter–Phelps model water quality parameters was proposed, where it could be obtained by using Lyapunov functions. The main proposal of this work was to design a new fractional observer to estimate BOD and compare it with the classic form of an extended Luenberger observer. It has been shown that the use of FOHO improves state estimation performance and that the fractional-order plays a critical role in evaluation by analysing the error by the IDAE criteria. This type of algorithm can be used and programmed in waterborne autonomous vehicles to estimate real-time BOD and monitor DO in case of distortions along the water source position. Our laboratory is currently working on this task.

#### *Future Work Perspectives*

The main contribution presented in this work is to provide an algorithm that estimates variables, such as BOD, since there are no on-line sensors for this type of variable [26,47]. This estimating algorithm (dynamic observers) in a robotic water vehicle could be programmed, and, then, the boat robot must include a position sensing unit; some works have shown the ability to have autonomy in the position control via GPS signal [26,48]

The vehicle maintains an automatic GPS position control (*L* longitude) that allows it to keep a constant speed ( *U* volume/time), similar to those of the river ( *Uriver* volume/time), as well as an intelligent system of obstacle avoidance that allows the vehicle to stay in the average radius of the river width; hence, *U* = *Uriver*, see Figure 6. In addition, we propose that these results be presented in a comparative study with renowned modern, robust techniques [29,30] where the performance at parametric and river spatial velocity presence is studied.

**Figure 6.** Outline of the application of environmental monitoring.

**Author Contributions:** Conceptualization, A.E.R.-M.; methodology, V.G.-H., P.A.L.-P. and O.H.-G.; writing—original draft preparation, Y.B.-T. and L.E.A.-S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
