**2.3 Modeling and linearization of the WWTP**

The use of the ASM1 model [8] is proposed to obtain the mathematical models (mainly for the biological model and the dissolved oxygen model). Among the main variables to be controlled in a WWTP are the control of the organic content in the effluent (BOD) and the control of nutrients (nitrogen and phosphorus). However, given the low concentrations found of these last two compounds in the wastewater to be treated, this chapter will focus on DO control in the bioreactor. This leads us to conclude that the proposed model can be simplified, eliminating for example the equations that have to do with the rate of change of nitrogen concentration and some other differential equations and variables that would not be necessary to take into account in this study.

**Figure 11** shows the configuration and ratio of the activated sludge WWTP inlet and outlet flows, where part of the effluent biomass (XR) is recycled to the bioreactor. The effluent from the bioreactor (QBS) feeds the clarifier (settler), used to separate substrate and biomass. On the other hand, part of the biomass in the clarifier is fed back to the bioreactor (QR), while the excess biomass (QW) is removed from the process.

**Figure 11.** *Configuration and flow ratio at the WWTP.*

*Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment… DOI: http://dx.doi.org/10.5772/intechopen.108816*

where, Qf = Inflow (300 l/min = 18 m3/h). QBS = Flow from bioreactor to settler (QBS = 1.5 Qf). QR = Biomass recirculation flow (QBS – Qf), approximately 50% del Qf. Qe = Output flow (Qe = 0.6QBS). QW = Residual biomass flow (QW = QBS – QR – Qe), approximately 20% of QR. QA = Airflow. V = Bioreactor volume (1.2 m3).

The mass balances for substrate, dissolved oxygen, and heterotrophic biomass concentrations in the bioreactor are represented by the differential equations Eqs. (20)-(22).

• Substrate balance:

$$\frac{d\mathbf{S}s}{dt} = \frac{\mathbf{Q}\_f}{V} \left( \mathbf{S} \mathbf{s}\_f - \mathbf{S}s \right) - \frac{\mu\_H}{\mathbf{Y}\_H} \left( \frac{\mathbf{S}s}{\mathbf{K}\_S + \mathbf{S}s} \right) \left( \frac{\mathbf{S}o}{\mathbf{K}\_{OH} + \mathbf{S}o} \right) \mathbf{X}\_H \tag{20}$$

• Oxygen balance:

$$\frac{d\text{So}}{dt} = \frac{\text{Q}\_f}{V}\text{So}\_f - \frac{\text{Q}\_f + \text{Q}\_R}{V}\text{So} + \frac{\text{Y}\_H - 1}{\text{Y}\_H}\mu\_H \left(\frac{\text{Sa}}{\text{K}\_S + \text{Sa}}\right) \left(\frac{\text{So}}{\text{K}\_{OH} + \text{Sa}}\right) \text{X}\_H + a \left(1 - e^{-\frac{\text{Q}\_d}{\text{b}}}\right) (\text{So}\_{\text{sat}} - \text{So}\_{\text{sat}}) \tag{21}$$

• Heterotrophic biomass balance:

$$\frac{dX\_H}{dt} = \frac{Q\_f}{V}X\_{H\sharp} - \frac{Q\_w}{V} \left(\frac{Q\_f + Q\_R}{Q\_w + Q\_R}\right)X\_H + \mu\_H \left(\frac{\text{S}\text{s}}{K\_\text{S} + \text{S}\text{s}}\right) \left(\frac{\text{So}}{K\_{OH} + \text{So}}\right)X\_H - b\_HX\_H \tag{22a}$$

These equations represent a nonlinear system, where *KLa* ¼ *a* 1 � *e* �*QA b* in accor-

dance with the foregoing, being QR, Qw y QA the manipulated variables, where,

Ss = Substrate concentration in the bioreactor.

Ssf = Substrate concentration in inflow.

So = Oxygen concentration.

Sof = DO concentration in the inflow.

XH = Heterotrophic biomass concentration.

XHf = Heterotrophic inflow biomass.

YH = Yield Heterotrophic biomass.

*μH*= Maximum specific biomass growth rate.

KS = Average Monod saturation constant for the substrate.

bH = Biomass decay rate.

KOH = Mean Monod saturation constant for oxygen for heterotrophs.

To simplify the system of equations, the following constants are formed:

$$\begin{aligned} D &= \frac{\mathbf{Q}\_f}{V}, D\_1 = \frac{\mu\_H}{Y\_H}, D\_3 = \frac{\mathbf{Q}\_f + \mathbf{Q}\_R}{V}, D\_4 = \frac{Y\_H - \mathbf{1}}{Y\_H} \mu\_H, D\_5 = \frac{\mathbf{Q}\_w}{V} \left(\frac{\mathbf{Q}\_f + \mathbf{Q}\_R}{\mathbf{Q}\_w + \mathbf{Q}\_R}\right) \\ K\_1 &= \frac{\mathbf{S}\mathbf{s}}{K\_S + \mathbf{S}\mathbf{s}}, K\_2 = \frac{\mathbf{S}\mathbf{s}}{K\_{OH} + \mathbf{S}\mathbf{s}} \end{aligned} \tag{22b}$$

When replacing, you have Eqs. (23)–(25):

$$\frac{d\mathbf{Ss}}{dt} = D\left(\mathbf{Ss}\_f - \mathbf{Ss}\right) - D\_1 K\_1 K\_2 X\_H \tag{23}$$

$$\frac{d\text{So}}{dt} = D\text{So}\_f - D\_3\text{So} + D\_4K\_1K\_2X\_H + a\left(\mathbf{1} - e^{\frac{-Q\_4}{b}}\right)(\text{So}\_{\text{sat}} - \text{So})\tag{24}$$

$$\frac{d\mathbf{X}\_H}{dt} = D\mathbf{X}\_{Hf} - D\_5\mathbf{X}\_H + \mu\_H \mathbf{K}\_1 \mathbf{K}\_2 \mathbf{X}\_H - b\_H \mathbf{X}\_H \tag{25}$$

The results obtained with the nonlinear model are shown in **Figure 12**, which shows the accumulation of biomass until the adaptation of the microorganisms in the bioreactor, once they begin to grow the substrate decreases. On the left without airflow and on the right supplying 3.3985 m3/h, for integration (blue), differentiation (red), and with oxygen consumption (green).

Eqs. (23)–(25), linearized by Taylor series, result in Eqs. (26)–(28):

$$\begin{split} \Delta \dot{\mathbf{S}} &= -(D + D\_1 K\_{2o} K\_3 X\_{H\_o}) \Delta \mathbf{S} - (D\_1 K\_{1o} K\_4 X\_{H\_o}) \Delta \mathbf{S} o \\ -(D\_1 K\_{1o} K\_{2o}) \Delta X\_H &\end{split} \tag{26}$$

$$
\Delta \dot{\mathbf{S}} = (D\_4 K\_{2o} K\_3 X\_{H\_o}) \Delta \mathbf{S} + \left[ D\_4 K\_{1o} K\_4 X\_{H\_o} - D\_3 - a \left( \mathbf{1} - e^{\frac{-Q\_{do}}{\theta}} \right) \right] \Delta \mathbf{S} \mathbf{o} + \dots \tag{27}
$$

$$
(D\_4 K\_{1o} K\_{2o}) \Delta X\_H + \left[ \frac{d}{b} (\text{So}\_{\text{sat}} - \text{So}\_o) e^{\frac{-Q\_{do}}{b}} \right] \Delta Q\_A
$$

$$
\begin{split}
\Delta \dot{X}\_H &= (\mu\_H K\_{2o} K\_3 X\_{H\_o}) \Delta \mathbf{S} + (\mu\_H K\_{1o} K\_4 X\_{H\_o}) \Delta \mathbf{S} \\
&+ (\mu\_H K\_{1o} K\_{2o} - D\_5 - b\_H) \Delta X\_H
\end{split} \tag{28}
$$

Here *Sso*, *Soo*, *XHo*, *QAo*, *K*1*<sup>o</sup>*, and *K*2*<sup>o</sup>* are the parameters and constants evaluated at the point of operation. The constants K3 y K4 correspond to the partial derivatives of K1 and K2, respectively. The incremental variables are: Δ*Ss* ¼ *Ss* � *Sso*, Δ*So* ¼ *So* � *Soo* y Δ*XH* ¼ *XH* � *XHo*.

**Figure 12.** *Concentrations nonlinear system Top: OD. Center: substrate. Bottom: biomass.*

*Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment… DOI: http://dx.doi.org/10.5772/intechopen.108816*

As a result, the linear model is obtained with Eq. (29):

$$
\begin{bmatrix}
\Delta \dot{\mathbf{S}} \\
\Delta \dot{\mathbf{S}} \\
\Delta \dot{\mathbf{X}}\_{H}
\end{bmatrix} = \begin{bmatrix}
\text{S} \mathbf{s}\_{1} & \text{S} \mathbf{s}\_{2} & \text{S} \mathbf{s}\_{3} \\
\text{S} \mathbf{o}\_{1} & \text{S} \mathbf{o}\_{2} & \text{S} \mathbf{o}\_{3} \\
\text{X}\_{H1} & \text{X}\_{H2} & \text{X}\_{H3}
\end{bmatrix} \begin{bmatrix}
\Delta \mathbf{S} \\
\Delta \mathbf{S} \\
\Delta \mathbf{X}\_{H}
\end{bmatrix} \\
+ \begin{bmatrix}
D\_{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & D\_{1} & \mathbf{0} & \mathbf{S} \mathbf{o}\_{4} \\
\mathbf{0} & \mathbf{0} & D\_{1} & \mathbf{0}
\end{bmatrix} \begin{bmatrix}
\Delta \mathbf{S} \mathbf{s}\_{f} \\
\Delta \mathbf{S} \mathbf{o}\_{f} \\
\Delta \mathbf{X}\_{Hf} \\
\Delta \mathbf{Q}\_{A}
\end{bmatrix} \\
\tag{29}
$$

where,

$$\begin{aligned} \text{Sr}\_{1} &= -(D + D\_{1}K\_{2o}K\_{3}X\_{H\_{s}}) \\ \text{Sr}\_{2} &= -(D\_{1}K\_{1o}K\_{4}X\_{H\_{s}}) \\ \text{Sr}\_{3} &= -(D\_{1}K\_{1o}K\_{2o}) \\ \text{Sr}\_{1} &= (D\_{4}K\_{2o}K\_{3}X\_{H\_{s}}) \\ \text{Sr}\_{2} &= \left[D\_{4}K\_{1o}K\_{4}X\_{H\_{s}} - D\_{3} - a\left(1 - e^{\frac{-Q\_{ds}}{b}}\right)\right] \\ \text{Sr}\_{3} &= (D\_{4}K\_{1o}K\_{2o}) \\ \text{Sr}\_{4} &= \left[\frac{a}{b}(\text{So}\_{\text{at}} - \text{So}\_{o})e^{\frac{-Q\_{ds}}{b}}\right] \\ X\_{H1} &= (\mu\_{H}K\_{2o}K\_{3}X\_{H\_{s}}) \\ X\_{H2} &= (\mu\_{H}K\_{1o}K\_{4}X\_{H\_{s}}) \\ X\_{H3} &= (\mu\_{H}K\_{1o}K\_{2o} - D\_{5} - b\_{H}) \end{aligned} \tag{30a}$$

The linear model obtained in incremental variables is Eq. (30).

$$
\begin{bmatrix}
\Delta \dot{\mathbf{s}} \\
\Delta \dot{\mathbf{o}} \\
\Delta \dot{\mathbf{x}}\_{H}
\end{bmatrix} = \begin{bmatrix}
9.7902 & 0.0748 & -0.0856
\end{bmatrix} \begin{bmatrix}
\Delta \mathbf{s} \\
\Delta \mathbf{s} \\
\Delta \mathbf{X}\_{H}
\end{bmatrix} \\
+ \begin{bmatrix}
0.25 & 0 & 0 & 0 \\
0 & 0.25 & 0 & 143.9339 \\
0 & 0 & 0.25 & 0
\end{bmatrix} \begin{bmatrix}
\Delta \mathbf{s}\_{f} \\
\Delta \mathbf{s}\_{f} \\
\Delta \mathbf{X}\_{H}
\end{bmatrix} \tag{30b}
$$

#### **3. Predictive control strategy**

#### **3.1 Concepts**

Model-based predictive control, MBPC, corresponds to a control strategy that uses a model of the system dynamics to predict the future behavior of the system over a finite time window called horizon. The predictive control strategy corresponds to an optimal control algorithm, which asks how best to control the system and calculates the future control action "u(t)" as a function of a penalty or cost function. The optimization of predictive control is limited to a moving time interval and is carried out continuously online.

Based on the model predictions and the actual measurement or estimation of the system state, the optimal control inputs are calculated, based on the defined control objective, and subject to the imposed constraints. After a time interval, the measurement, estimation, and calculation processes are repeated using a shifted horizon.

Among the advantages is that the controller can anticipate future disturbances and can be applied to systems with high nonlinearities without the need to perform system liberalization and finally can explicitly consider operational, physical, or safety constraints of the system.

The following are the three ingredients of an optimal control problem.

Model: the model will be used to predict the evolution of the state for a given sequence of inputs, as shown in Eq. (31):

$$\mathbf{x}\_{t+1} = f(\mathbf{x}\_t, \mathbf{u}\_t) \tag{31}$$

Given a sequence of inputs, if you have a model, you can predict the evolution of the trajectory, that is, the evolution of the state in the future.

**Objetive** is used to assign a cost to a given trajectory and qualifies how good a given trajectory is for the task to be performed. The objective is written as a cost function, which depends on the sequence of states x and inputs u, mapping it to a scalar, as shown in Eq. (32):

$$J(\mathbf{x}\_{1:T}, \boldsymbol{u}\_{1:T}) = \sum\_{t \in [T]} \mathbf{g}\_t(\mathbf{x}\_t, \boldsymbol{u}\_t) \tag{32a}$$

where,

$$\begin{array}{c} \pi\_{1:T} := (\pi\_1, \dots, \pi\_T); \qquad \mu\_{1:T} := (\mu\_1, \dots, \mu\_T) \end{array}$$

**Restrictions** codify the allowed domains for states and entries:

$$\begin{aligned} \boldsymbol{\chi}\_{t+1} &= \boldsymbol{f}(\boldsymbol{\chi}\_t, \boldsymbol{u}\_t), \forall t \in [T - \mathbf{1}] \\ \boldsymbol{\chi}\_t &\in \mathcal{X}\_t, \boldsymbol{u}\_t \in \mathcal{U}\_t; \forall t \in [T] \\ \boldsymbol{\chi}\_1 &= \boldsymbol{\chi}\_{init} \end{aligned} \tag{32b}$$

#### **3.2 Solution of the optimization problem**

In general, there is no closed solution to this type of problem, and in general numerical methods must be used to solve MATLAB and Python.

Using an optimal control strategy does not guarantee success due to several reasons. Firstly, the prediction model will always have errors, since it is an abstraction of the real system. Errors accumulate over time resulting in divergent predictions. If the control sequence u is applied to the open-loop prediction model, it will follow a different trajectory than the actual one and the task will not be carried out accurately. Secondly, even if you have a perfect model, knowing how the system behaves to the inputs that are applied, but you have a very long task horizon, if you want to apply an optimal control strategy to the entire sequence that is too long, you have the problem of having many time steps to consider in the optimization. Thus, you have a problem that is too long and too difficult to solve in a limited amount of time.

#### **3.3 Model-based predictive control**

Model-based predictive control adds an idea to the optimization problem explained in the previous section. Instead of optimizing for the complete trajectory in the future, the following receding horizon control is done:

*Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment… DOI: http://dx.doi.org/10.5772/intechopen.108816*

You start at the current time step, considering only the current state.

From there, you only look ahead and only consider a limited preview, e.g. the next 50 steps.

Only the first entry is applied, then plan again.

**Figure 13** shows a general scheme of the predictive control strategy, where it is observed that first the state is measured, then the task is assigned to the optimizer, to optimize for example the next 50 steps in time in the future, and from these 50 steps in time, the optimal solution is taken, taking only the first input to be applied to the system, to then measure the new state, that is, re-planning from there. The advantage is that errors can be taken into account because replanning introduces feedback and reduces the size of the problem, because you only need to plan for the next 50 time steps, for example.

**Figure 14** shows a general scheme corresponding to the predictive control strategy, MBPC. The method consists of a computerized algorithmic control software based on optimization, where the best control policy is chosen so that it meets a defined criterion, for example, after 10 samples is desired that the output y is on the set point.

**Figure 14** shows that past information is available, past control policy data (u(t-1), u(t-2) … ) and past process output information (y(t-1), y(t-2) … ). Two cases are shown on how the process output will evolve in the future for two future inputs u to the process and the reference output r is shown. N2 corresponds to the prediction

**Figure 13.**

*General block diagram of a MBPC strategy.*

**Figure 14.** *General scheme of a MBPC strategy.*

horizon and corresponds to one of the MBPC design parameters. In this strategy, u is selected such that the shaded area corresponding to the evolution of the process output in the future is minimal.

The MBPC algorithm can be described as follows:

• **Input:** Objective (cost function) J, Dynamics model f, horizon T, initial guess *u*^1:*T*;

1.*u*1:*<sup>T</sup> u*^1:*T*;

2.**While** task not completed do

*xinit* Get current state (); *u*<sup>1</sup>:*<sup>T</sup>* Solve optimization problem *u first u*ð Þ <sup>1</sup>:*<sup>T</sup> ApplyInput u*ð Þ to the system

3.*Start again*

The Diophantine equations propose the generalized predictive control (GPC) strategy, while the EPSAC strategy [1] proposes filtering techniques, where the prediction of the process output is given by:

$$\{\mathbf{y}(\mathbf{t} + \mathbf{k}/\mathbf{t}), \mathbf{k} = \mathbf{1}, \mathbf{2}, \dots, \mathbf{N2}\}\tag{33a}$$

and is based on available measurements of the control input u and the output y at instant t:

$$\left\{ \mathbf{y}(\mathbf{t}), \mathbf{y}(\mathbf{t}-\mathbf{1}) \dots, \mathbf{u}(\mathbf{t}-\mathbf{1}), \mathbf{u}(\mathbf{t}-\mathbf{2}), \dots \right\} \tag{33b}$$

What the control algorithm must calculate is the last value of the input u(t). Similarly, future values of the control input u must be postulated:

$$\{\mathbf{u}(\mathbf{t}/\mathbf{t}), \mathbf{u}(\mathbf{t}+\mathbf{1})/\mathbf{t}, \mathbf{u}(\mathbf{t}+\mathbf{2}/\mathbf{t})\dots\}\tag{33c}$$

• Output prediction: The process output prediction is given by two terms, namely the model output prediction and the disturbance prediction, as shown in Eq. (33):

$$y(t + k/t) = \mathbf{x}(t + k/t) + n(t + k/t)$$

Model prediction can be obtained using the parallel method as it is simpler and is used in stable processes, as shown in **Figure 15**.

• **Disturbance prediction:** In Eq. (34), the disturbance prediction can be obtained by defining a hypothetical filter nf (the signal does not exist), and this value must be stored in the computer memory, as it will be needed at time instant t + 1:

$$\begin{split} n\_f &= \frac{D(q^{-1})}{C(q^{-1})} n(t) \\ &= -c\_1 n\_f(t-1) - c\_2 n\_f(t-2) - \dots \\ &\dots + n(t) + d\_1 n(t-1) + d\_2 n(t-2) + \dots \end{split} \tag{34}$$

*Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment… DOI: http://dx.doi.org/10.5772/intechopen.108816*

**Figure 15.** *Parallel model prediction.*

All the information in Eq. (34) is known, it is from the past and is in the computer's memory, and the value of n(t) can be found from (1) for k = 0, as shown in Eq. (35):

$$m(t) = \mathcal{y}(t) - \mathcal{x}(t) \tag{35}$$

The future prediction of nf is given by Eq. (36):

$$n\_{\!\!\!/}(t+k/t) \equiv 0, k = 1 \dots N\_2 \tag{36}$$

In Eq. (37), the average of this signal is zero (0), because,

$$m\_f(t) = \frac{D(q^{-1})}{C(q^{-1})} \frac{C(q^{-1})}{D(q^{-1})} e(t) = e(t) \tag{37}$$

White noise is an uncorrelated signal, so it cannot be predicted (μe = 0), while n(t) corresponds to a correlated colored noise signal, where there is trend, then its value can be predicted (the prediction corresponds to the mean, μe 6¼ 0). The forward prediction of n(t) is given by Eq. (38):

$$\begin{aligned} \text{for } k &= \mathbf{1}:N\_2\\ n(t+k/t) &= \frac{\mathbf{C}(q^{-1})}{D(q^{-1})} n\_\uparrow(t+k/t) \end{aligned} \tag{38}$$

The best prediction of the future random n signal is obtained using Eq. (39):

$$\begin{split} n(t+k/t) &= d\_1 n(t+k-1/t) - d\_2 n(t+k-2/t) - \dots \\ &+ n\_f(t+k/t) + c\_1 n\_f(t+k-1/t) \end{split} \tag{39}$$

• **Base/Optimizing response:** according to the superposition principle, the forward response of the system can be written as Eq. (40):

$$y(t + k/t) = y\_{base}(t + k/t) + y\_{optimize}(t + k/t) \tag{40}$$


**Figure 16** shows the past, base, and sought control actions u(t + k/t). In order to reduce the number of computations and the size of the matrices, another design parameter is introduced, the control horizon Nu.

The control algorithm finds by optimization (minimizing the prediction errors), the future values of δu(t + k/t)(vector with Nu elements), adds them to the control signal ubase, to obtain the values of u(t + k/t).

yoptimize(t + k/t) can be calculated using the impulse response coefficients up to the value of Nu, with the last coefficient given by the value of the step response coefficient given by Eq. (41):

$$\begin{split} \mathbf{y}\_{\text{optimize}}(t+k/t) &= h\_k \delta u(t/t) + h\_{k-1} \delta u(t+1/t) + \dots \\ &\dots + \mathbf{g}\_{k-N\_u+1} \delta u(t+N\_u-\mathbf{1}/t) \end{split} \tag{41}$$

According to [1], the basic equation for the EPSAC algorithm can be written in Eq. (42):

$$\mathbf{Y} = \overline{\mathbf{Y}} + \mathbf{G}.\mathbf{U} \tag{42}$$

$$\left\{ \begin{array}{l} Y = \left[ y(t + N\_1/t) \dots y(t + N\_2/t) \right]^T \\\\ \overline{Y} = \left[ y\_{\text{bare}}(t + N\_1/t), \dots y\_{\text{bare}}(t + N\_2/t) \right]^T \\\\ U = \left[ \delta u \left( \frac{t}{t} \right) \dots \delta u (t + N\_u - \mathbf{1}/t) \dots \right]^T \end{array} \right\}$$

The objective of the predictive controller consists of finding the control vector, {u (t + k/t), k = 0 … N2–1}, that minimizes the following cost function, as shown in Eq. (43):

$$\sum\_{k=N\_1}^{N\_2} \left[ r(t+k/t) - y(t+k/t) \right]^2 \tag{43}$$

**Figure 17** shows the error between the forward trajectory of the process and the postulated reference.

EPSAC proposes the following solution for an unconstrained linear system, as shown in Eqs. (44) and (45):

$$\sum\_{k=N\_1}^{N\_2} \left[ r \left( t + \frac{k}{t} \right) - \chi \left( t + \frac{k}{t} \right) \right]^2 = \left[ \mathbf{R} - \overline{\mathbf{Y}} - \mathbf{G} \, \mathbf{U} \right]^T \left[ \mathbf{R} - \overline{\mathbf{Y}} - \mathbf{G} \, \mathbf{U} \right] \tag{44}$$

**Figure 16.** *u, past, base, and optimum control actions.*

*Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment… DOI: http://dx.doi.org/10.5772/intechopen.108816*

**Figure 17.** *Objective of the predictive control strategy.*

$$\boldsymbol{U}^\* = \left[\boldsymbol{G}^T \boldsymbol{G}\right]^{-1} \left[\boldsymbol{G}^T (\mathbf{R} - \overline{\mathbf{Y}})\right] \tag{45}$$

The actual control action applied to the process is shown in Eq. (46):

$$u(\mathbf{t}) = u\_{base}(\mathbf{t}/\mathbf{t}) + \delta u(\mathbf{t}/\mathbf{t}) = u\_{base}(\mathbf{t}/\mathbf{t}) + U^\*(\mathbf{1})\tag{46}$$

#### **3.4 Predictive control of DO concentration in the WWTP**

The basic MPC strategy consists of estimating the error of the future to determine the value of the control signal, predicting the future outputs, by means of the past outputs and the control signals. GPC control will be used, starting from the CARIMA prediction model based on a transfer function (TF) model [9], Eq. (44):

$$y(t) = \varkappa(t) + n(t) = \frac{B(z^{-1})}{A(z^{-1})} u(t-1) + \frac{C(z^{-1})}{\Delta D(z^{-1})} e(t) \tag{47}$$

where,

x(t) = process transfer function, n(t) = disturbance transfer function, Δ = integrator in the disturbance, *e t*ð Þ = white noise, and *u*(*t-1*) = control action (intrinsic delay in the process). With A(*z*�*<sup>1</sup>* ), *B*(*z*�*<sup>1</sup>* ), *C*(*z*�*<sup>1</sup>* ), and *D*(*z*�*<sup>1</sup>* ), polynomials are obtained by identification. Following the procedure described in (Aguado A), Eq. (47) is transformed into Eq. (48):

$$
\hat{y}(t+i/t) = G\_i \Delta u(t+i-1) + \Gamma\_i \Delta u^f(t-1) + F\_i y^f(t) \tag{48}
$$

This expression allows us to know at instant t, the value of the predicted output at instant *t+i*. Here, *Gi* yΓ*<sup>i</sup>* y *Fi* are polynomials FIR in *z*�*<sup>1</sup>* . Δ*u<sup>f</sup>* and *y f* are filtered inputs and outputs, respectively. The vector expression of the prediction model is shown in Eq. (49):

$$Y = Gu + \Gamma \Delta U^f + FY^f \tag{49a}$$

where,

$$
\Delta \mathbf{U}^{\zeta} = \left[ \Delta \boldsymbol{\mu}^{\zeta} (t - 1) \dots \Delta \boldsymbol{\mu}^{\zeta} (t - 2) \dots \Delta \boldsymbol{\mu}^{\zeta} (t - nt) \right]^T \tag{49b}
$$

$$
\mathbf{Y}^{\zeta} = \left[ \mathbf{y}^{\zeta}(t) \dots \mathbf{y}^{\zeta}(t - 1) \dots \mathbf{y}^{\zeta}(t - nu) \right]^T
$$

The N present and future control actions are calculated from the minimization of the quadratic cost index, as shown in Eq. (50):

$$J(u) = \sum\_{i=1}^{N} a\_i [\dot{\mathbf{y}}(t + i/t) - w(t + i)]^2 + \sum\_{j=1}^{N} \lambda\_j [\Delta u(t + j - 1)]^2 \tag{50}$$

The index in vector form is shown in Eq. (51):

$$J(u) = \left(Y - \mathcal{W}\right)^T a (Y - \mathcal{W}) + u^T \lambda u \tag{51}$$

*α* and *λ* are diagonal matrices of NxN, with:

$$Y = \begin{bmatrix} \hat{y}(t+1/t) \dots \hat{y}(t+2/t) \dots \hat{y}(t+N/t) \end{bmatrix}^T \tag{52}$$

$$\boldsymbol{u} = \begin{bmatrix} \Delta \boldsymbol{u}(t) \dots \Delta \boldsymbol{u}(t+1) \dots \Delta \boldsymbol{u}(t+N-1) \end{bmatrix}^T \tag{53}$$

$$\mathcal{W} = \left[\mathcal{W}(t+1) \dots \mathcal{W}(t+2) \dots \mathcal{W}(t+N)\right]^T \tag{54}$$

Substituting Eq. (49) into equation Eq. (51), we obtain the expression that calculates the N future changes of the control action that minimizes the quadratic cost index, as shown in Eq. (55):

$$\boldsymbol{u} = \left(\boldsymbol{G}^T \boldsymbol{a} \mathbf{G} + \boldsymbol{\lambda}\right)^{-1} \boldsymbol{G}^T \boldsymbol{a} \left(\boldsymbol{W} - \boldsymbol{\Gamma} \boldsymbol{\Delta} \boldsymbol{U}^f - \boldsymbol{F} \mathbf{Y}^f\right) \tag{55}$$

Although the N control actions are computed, the linear controller only implements the first Δu(t), converting u into:

$$
\Delta u(t) = h\mathcal{W} - h\Gamma \Delta \mathcal{U}^f - hF\mathcal{Y}^f \tag{56}
$$

where h is the first row of *<sup>G</sup><sup>T</sup>α<sup>G</sup>* <sup>þ</sup> *<sup>λ</sup>* � ��<sup>1</sup> *G<sup>T</sup>α*.

The Z-transform of Δ*u t*ð Þ allows us to obtain the configuration of the GPC controller, as shown in **Figure 18**, represented by Eq. (57):

$$u(\mathbf{z}) = \frac{T(\mathbf{z}^{-1})}{(T(\mathbf{z}^{-1}) + R(\mathbf{z}^{-1})\mathbf{z}^{-1})\Delta} \left[ (H(\mathbf{z})W(\mathbf{z}) - \frac{\mathbf{S}(\mathbf{z}^{-1})}{T(\mathbf{z}^{-1})}y(\mathbf{z}) \right] \tag{57}$$

From the linear system found, the gain matrix of the KMPC controller was found using the Model Predictive Control Toolbox:

*KMPC* <sup>¼</sup> <sup>½</sup> <sup>0</sup>*:*0058156 1*:*9218e‐<sup>006</sup> ‐6*:*0503e‐<sup>006</sup> ‐1*:*209e‐007 1*:*8134e‐006 2*:*4445e‐<sup>006</sup> �

**Figure 18.** *GPC controller (green).*

*Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment… DOI: http://dx.doi.org/10.5772/intechopen.108816*

#### **Figure 19.**

*Left: no valve restrictions. Right: with valve restrictions.*

#### **Figure 20.**

*GPC controller, with disturbances.*

To determine the robustness of the GPC, simulations were performed with and without perturbation for different reference values of the DO. In all of them, the performance achieved was adequate. By way of explanation, the case is reproduced when there is no disturbance in **Figure 19**, where the set point (red) of 2, 1, 4, and 1 mg/l of DO is reached in the upper part, while the lower part shows the control actions sent to the valve that regulate the airflow.

If disturbances to the controller and to the process output are involved, the block diagram shown in **Figure 20** [10] is obtained.

When the plant is subjected to disturbances, as shown in **Figure 21**, a faster control action is observed, causing the response to have a slight overshoot. However, the DO set point is reached.

The MPC controller design methodology has a great acceptance in the control of the different variables that are present in the wastewater treatment processes, thanks to its great robustness and the availability of software that allows the development and execution of the algorithms necessary for its implementation.

This chapter shows the possibility of controlling the concentration of DO in a WWTP from a GPC, whose algorithm has been programmed in MATLAB and executed in a LabVIEW supervisory system.

#### **4. Conclusions**

The MPC controller design methodology is widely accepted in the control of the different variables that occur in wastewater treatment processes, thanks to its great

**Figure 21.** *Left: no valve restrictions. Right: with valve restrictions.*

robustness and the availability of software that allows the development and execution of algorithms necessary for its implementation. One of its advantages is to guarantee lower energy consumption due to the fact that the steps of the control action are smaller and of an anticipated type, resulting in lower opening values of the airflow valve.

This work shows how to control the concentration of DO in a WWTP from a GPC, whose algorithm has been programmed in MATLAB and executed in a LabVIEW SCADA system. It was based on the step response of the system and the description in state space, with or without restrictions in the airflow (manipulated variable) and/or in the concentration of DO (controlled variable), either in the presence or not of disturbances in the system (concentrations of dissolved oxygen, substrate, and biomass in the inlet flow to the bioreactor). In the different simulations and tests in the implemented controller, it was found that when there are no restrictions neither in the airflow valve nor in the DO sensor, even if disturbances are present, the set point is adequately reached in all cases. On the contrary, when physical restrictions are present in the airflow and/or the oxygen concentration, there are moments in which the desired value of the controlled variable is not adequately reached.
