Predictive Control, a Strategy for Dissolved Oxygen Control in a Wastewater Treatment Plant

Abstract

This chapter presents a strategy for controlling the concentration of dissolved oxygen (DO) in the bioreactor of a pilot wastewater plant (WWTP). The control strategy being developed is model-based predictive control (MPC). To apply the control algorithm, the estimation of the oxygen transfer function (KLa) is first performed, then the model and linearization technique are determined and finally the MBPC controller is designed. The results are simulated in MATLAB® and executed in the plant control and supervision system Supervisory Control and Data Acquisition (SCADA) in LabVIEWTM. This chapter is organized as follows: Section 1 presents a brief introduction, and Section 2 determines and describes the model to be used and its respective linearization, as well as the results obtained for the KLa parameter. Finally, Section 3 describes the design methodology of the generalized predictive control (GPC) controller proposed by Clarke, using the Model Predictive Control Toolbox and the EPSAC strategy developed by De Keyser. It should be noted that the simulations in each of the sections were performed in MATLAB® and executed in the control and supervision system with the MATLAB® script interface in LabVIEWTM.

1. Introduction

The WWTP treats water from a pipeline of a sector in a community, and the activated sludge used comes from the wastewater plant (WWTP) of a nearby residential complex. Its main objective is to serve as a didactic and experimental means for learning, simulation, and research in the area of control of this biological process. Oxygen transfer is an important factor for aerobic biological water purification processes; hence, the need to obtain a good estimation of this parameter. Consequently, this chapter discusses a strategy for the control of dissolved oxygen (DO). The proposed control strategy is the Model-Based Predictive Control (MPC). In this strategy, a model of the process is used to predict how the system will behave under a proposed sequence of control actions u(t). Traditionally, a linear model is used where an optimal control action can be calculated; however, when nonlinear models are used in the algorithm, a numerical search algorithm is used to calculate an optimal sequence uopt(t).

2. Process and modeling of wastewater systems

In this part, the generalities of the process are described, and some techniques for obtaining important parameters in the modeling process of a wastewater treatment plant are proposed and a linear version of the ASM1 model is proposed.

2.1 General concepts

In the application of predictive control techniques based on models, it is necessary to have quality models; otherwise, you will have erroneous predictions and the control system will not be good. The predictions are carried out using a model of the process, with which the specification of a reference trajectory is obtained. De Keyser proposed some process models that can be used in the implementation of predictive control strategies (EPSAC Model) [1]. The output of the model can be seen in Figure 1.

Figure 1.

Process model representation.

Figure 1 shows schematically the representation of the system taking into account the effect of noise n (colored noise), where the response to an input u can be obtained with the Eq. (1):

$$y\left(t\right)=x\left(t\right)+n\left(t\right)$$E1

The disturbance is represented by n(t) and corresponds to colored noise, a random signal that removes the offset. x is the result of the process model, while y corresponds to the measured value. The colored noise signal n(t) can be represented by Eq. (2):

$$n\left(t\right)=\frac{C\left(q^{-1}\right)}{D\left(q^{-1}\right)}e\left(t\right)$$E2

With, e(t): White noise (uncorrelated noise with zero mean value).

The transfer function (TF) between e(t) and n(t) describes the disturbance class and corresponds to the noise filter and is given by Eq. (3):

$$\frac{C\left(q^{-1}\right)}{D\left(q^{-1}\right)}=\frac{\left(1+{\mathit{cq}}^{-1}\right)}{\left(1+{\mathit{dq}}^{-1}\right)\left(1-q^{-1}\right)}$$E3

The output of model x can be written as Eq. (4):

$$x\left(t\right)=f\left[x\left(t-1\right),x\left(t-2\right),\dots ,u\left(t-1\right),u\left(t-2\right),\dots \right]$$E4

A special kind of EPSAC Model is the CARIMA Model [2]. This model links the process model and the disturbance model, as shown in Eq. (5):

$$\begin{array}{c}A\left(q^{-1}\right)y\left(t\right)=B\left(q^{-1}\right)u\left(t\right)+\frac{C\left(q^{-1}\right)}{1-q^{-1}}\\ y\left(t\right)=\frac{B\left(q^{-1}\right)}{A\left(q^{-1}\right)}u\left(t\right)+\frac{C\left(q^{-1}\right)}{A\left(q^{-1}\right)\left(1-q^{-1}\right)}\\ y\left(t\right)=x\left(t\right)u\left(t\right)+n\left(t\right)e\left(t\right)\end{array}$$E5

x(t) represents a pulse transfer function and B(q-1) and A(q-1) are given by Eq. (6):

$$\begin{array}{c}B\left(q^{-1}\right)=b_1{q}^{-1}+\dots +b_{\mathit{nb}}{q}^{-\mathit{nb}}\\ A\left(q^{-1}\right)=1+a_1{q}^{-1}+\dots +b_{\mathit{na}}{q}^{-\mathit{na}}\end{array}$$E6

The model can also be obtained using neural networks or Fuzzy modeling techniques. The model obtained for the wastewater treatment plant is presented below.

2.2 Oxygen transfer function estimation, KLa

In a wastewater treatment system, oxygen must be available at a rate equivalent to the oxygen demand load exerted by the wastewater entering the plant. The process consists of bringing the wastewater in contact with oxygen, transferring it across the gas-to-liquid interface to dissolve it in the liquid, and then transferring the dissolved oxygen through the liquid to the microorganisms. The dynamics of DO concentration change can be represented by Eq. (7):

$$\frac{\mathit{dSo}}{\mathit{dt}}=K_{L}a\left({\mathit{So}}_{\mathit{sat}}-\mathit{So}\right)$$E7

where So is the dissolved oxygen concentration DO, So_sat is the oxygen saturation concentration, K_L is the oxygen transfer coefficient, and a is the total interfacial contact area per unit volume of liquid. Since it is admittedly impossible to measure the interfacial area a, the total term KLa is estimated [3]. To obtain the KLa parameter, the integration, differentiation, and with oxygen consumption methods were initially used in transient regime [4], assuming negligible biomass concentrations, proceeding first to deoxygenate the wastewater, bringing the DO to a value close to zero. Then aeration is restarted, measuring the increase in DO concentration over time using the DO sensor. In Figure 2, it can be seen how the DO concentration increases at different airflow rates until reaching a steady state.

Figure 2.

DO concentration. Left: Function on time. Right: Function on time and airflow.

The determination of KLa with the integration method proposes to separate variables and integrate Eq. (7). Assuming that KLa does not depend on the sampling time, it is obtained Eq. (8):

$$ln\left({\mathit{So}}_{\mathit{sat}}-\mathit{So}\right)=-K_L\mathit{at}$$E8

With this equation, a straight line can be obtained from a semi-logarithmic graph (Sosat - So) as a function of time, where KLa is its slope, as shown in Eq. (9):

$$K_{L}a=-{\frac{ln\left(C_f/C_i\right)}{t_f-t_i}}^{\ast }60$$E9

where ti and tf are the initial and final time selected for the slope, Ci = (Sosat - So) in ti and Cf = (Sosat - So) in tf. It is multiplied by 60 to convert from minutes to hours, when the samples have been taken in minutes. The results obtained can be seen in Figure 3. Here, it is observed that the greater the airflow, the greater the slope (KLa).

Figure 3.

Slope (-KLa). Left: Function on time. Right: Function on time and airflow.

With these results, Figure 4 shows the curve obtained from KLa, whose behavior resembles a first-order system.

Figure 4.

KLa integration method.

The differentiation method is based on the difference of (Sof - So), where Sof is the last sampled OD data. KLa is obtained from Eq. 10:

$$K_{L}a={\frac{d\left(ln\left({\mathit{So}}_f-\mathit{So}\right)\right)}{\mathit{dt}}}^{\ast }60$$E10

The results are shown in Figure 5, where it is observed that the slope for the different levels of airflow can be obtained in the first 300 seconds, during which period the slope remains constant.

Figure 5.

Slope ln(Sof - So). Left: Function on time. Right: Function on time and airflow.

In a similar way, the KLa curve is obtained, as a function of the airflow, observing a behavior similar to that of a first-order system, as shown in Figure 6.

Figure 6.

KLa differentiation method.

In activated sludge processes, when the wastewater under treatment has a significant concentration of microorganisms, it is necessary to take into account the oxygen consumption (respiration) [4], since the mass balance presented in Eq. (7) is affected as show in Eq. (11):

$$\frac{\mathit{dSo}}{\mathit{dt}}=K_{L}a\left({\mathit{So}}_{\mathit{sat}}-\mathit{So}\right)-R$$E11

where R is the rate of oxygen utilization (respiration rate). Then Eq. (11) can be written as Eq. (12):

$$\frac{\mathit{dSo}}{\mathit{dt}}=\left(K_L{\mathit{aSo}}_{\mathit{sat}}-R\right)-K_L\mathit{aSo}$$E12

Eq. (12) indicates that the derivative term dSo/dt as a function of So provides a straight line whose slope is equal to KLa and its cutoff point with the ordinate is (KLa.Sosat - R), value from which R can also be calculated. Figure 7 shows the derivative dSo/dt with respect to So, the value of KLa (slope), and the value of R (cutoff point on the ordinate).

Figure 7.

Slope (dSo/dt). Left: Function on time. Right: Function on time and airflow.

Figure 8 shows that the behavior of KLa as a function of the airflow obtained by this method also follows a dynamic similar to that of a first-order system.

Figure 8.

KLa oxygen consumption method.

The oxygen transfer function KLa describes the rate at which oxygen is transferred to the activated sludge by the aeration system. This function is nonlinear and depends on several factors, the main one being the airflow rate. Eq. (7) can be written as Eq. (13):

$$\frac{\mathit{dSo}}{\mathit{dt}}=K_{L}a\left(q_A\left(t\right)\right)\left({\mathit{So}}_{\mathit{sat}}-\mathit{So}\right)$$E13

where qA(t) is the airflow rate entering the bioreactor. Here, it is assumed that KLa depends on the nonlinearity of the airflow rate. A typical function of KLa is shown in: Figure 4, Figure 6, and Figure 8, where it is observed that the slope of KLa changes when the airflow rate changes [3]. Based on the behavior of KLa observed in these figures, it can be assumed that it corresponds to a first-order system that can be represented by the differential equation, Eq. (14):

$$b\frac{{\mathit{dK}}_{L}a\left(q_A\left(t\right)\right)}{d\left(q_A\left(t\right)\right)}+K_{L}a\left(q_A\left(t\right)\right)=\mathit{au}\left(q_A\left(t\right)\right)$$E14

where b is the time constant of the system, a is the gain of the system, and u(qA(t)) is the input of the system. Applying the Laplace transform to Eq. (14), the transfer function of the system is modeled with Eq. (15):

$$\frac{K_{L}a\left(s\right)}{u\left(s\right)}=\frac{a}{\mathit{bs}+1}$$E15

The step response of the system is obtained with of the inverse Laplace transform, as shown in Eq. (16):

$$K_{L}a\left(q_A\left(t\right)\right)=a\left(1-e^{\frac{-q_A\left(t\right)}{b}}\right)$$E16

The values of the parameters a and b are determined by identification using the Smith model [5]. In this model, it is proposed that the steady state gain a can be easily obtained from the graph of the step response of the system, while the system constant b corresponds approximately to the time in which 63.2% of the final steady state response is reached [6].

The gains (a) at steady state are determined to be 30.05, 61.38, and 61.66 for the integration, differentiation, and oxygen consumption methods, respectively. On the other hand, the constants of qA(t) are 1.386, 0.214, and 0.2042. These constants are calculated by linear interpolation or by a cubic spline. The identified models for integration, differentiation, and with oxygen consumption correspond to the following transfer functions, as shown in Eq. (17), Eq. (18), and Eq. (19):

$$\frac{K_{L}a\left(s\right)}{u\left(s\right)}=\frac{30.05}{1.386s+1}$$E17

$$\frac{K_{L}a\left(s\right)}{u\left(s\right)}=\frac{61.38}{0.214s+1}$$E18

$$\frac{K_{L}a\left(s\right)}{u\left(s\right)}=\frac{61.66}{0.2042s+1}$$E19

The verification of the models is performed by comparing the process signal (red) with the model signal (blue), as shown in Figure 9, in which it can be seen that there is a good tracking and consequently the error is low. The value of KLa is found by applying Eq. (16), for a given value of airflow qA(t).

Figure 9.

Transfer function KLa. Left: integration. Center: differentiation. Right: with oxygen consumption.

For data acquisition, supervision, and control of the WWTP, a Supervisory Control and Data Acquisition (SCADA) system is designed in LabVIEWTM [7]. The block diagram and the front panel of the virtual instrument (VI) in the control and monitoring system are shown in Figure 10.

Figure 10.

Left: block diagram integration. Right: front panel of the virtual instrument (VI).

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.

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{\mathit{dSs}}{\mathit{dt}}=\frac{Q_f}{V}\left({\mathit{Ss}}_f-\mathit{Ss}\right)-\frac{{\mu }_H}{Y_H}\left(\frac{\mathit{Ss}}{K_S+\mathit{Ss}}\right)\left(\frac{\mathit{So}}{K_{\mathit{OH}}+\mathit{So}}\right)X_H$$E20

• Oxygen balance:

  $$\frac{\mathit{dSo}}{\mathit{dt}}=\frac{Q_f}{V}{\mathit{So}}_f-\frac{Q_f+Q_R}{V}\mathit{So}+\frac{Y_H-1}{Y_H}{\mu }_H\left(\frac{\mathit{Ss}}{K_S+\mathit{Ss}}\right)\left(\frac{\mathit{So}}{K_{\mathit{OH}}+\mathit{So}}\right)X_H+a\left(1-e^{\frac{-Q_A}{b}}\right)\left({\mathit{So}}_{\mathit{sat}}-\mathit{So}\right)$$E21

• Heterotrophic biomass balance:

$$\frac{{\mathit{dX}}_H}{\mathit{dt}}=\frac{Q_f}{V}{X}_{\mathit{Hf}}-\frac{Q_w}{V}\left(\frac{Q_f+Q_R}{Q_w+Q_R}\right)X_H+{\mu }_H\left(\frac{\mathit{Ss}}{K_S+\mathit{Ss}}\right)\left(\frac{\mathit{So}}{K_{\mathit{OH}}+\mathit{So}}\right)X_H-b_H{X}_H$$E22a

These equations represent a nonlinear system, where $K_{L}a=a\left(1-e^{\frac{-Q_A}{b}}\right)$ in accordance 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.

${\mu }_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{array}{l}D=\frac{Q_f}{V},D_1=\frac{{\mu }_H}{Y_H},D_3=\frac{Q_f+Q_R}{V},D_4=\frac{Y_H-1}{Y_H}{\mu }_H,D_5=\frac{Q_w}{V}\left(\frac{Q_f+Q_R}{Q_w+Q_R}\right)\\ K_1=\frac{\mathit{Ss}}{K_S+\mathit{Ss}},K_2=\frac{\mathit{So}}{K_{\mathit{OH}}+\mathit{So}}\end{array}$$E22b

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

$$\frac{\mathit{dSs}}{\mathit{dt}}=D\left({\mathit{Ss}}_f-\mathit{Ss}\right)-D_1{K}_1{K}_2{X}_H$$E23

$$\frac{\mathit{dSo}}{\mathit{dt}}={\mathit{DSo}}_f-D_3\mathit{So}+D_4{K}_1{K}_2{X}_H+a\left(1-e^{\frac{-Q_A}{b}}\right)\left({\mathit{So}}_{\mathit{sat}}-\mathit{So}\right)$$E24

$$\frac{{\mathit{dX}}_H}{\mathit{dt}}={\mathit{DX}}_{Hf}-D_5{X}_H+{\mu }_H{K}_1{K}_2{X}_H-b_H{X}_H$$E25

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).

Figure 12.

Concentrations nonlinear system Top: OD. Center: substrate. Bottom: biomass.

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

$$\begin{array}{l}\mathrm{\Delta }\stackrel{·}{\mathit{Ss}}=-\left(D+D_1{K}_{2o}{K}_3{X}_{H_o}\right)\mathrm{\Delta }\mathit{Ss}-\left(D_1{K}_{1o}{K}_4{X}_{H_o}\right)\mathrm{\Delta }\mathit{So}\\ -\left(D_1{K}_{1o}{K}_{2o}\right)\mathrm{\Delta }{X}_H\end{array}$$E26

$$\begin{array}{l}\mathrm{\Delta }\stackrel{·}{\mathit{So}}=\left(D_4{K}_{2o}{K}_3{X}_{H_o}\right)\mathrm{\Delta }\mathit{Ss}+\left[D_4{K}_{1o}{K}_4{X}_{H_o}-D_3-a\left(1-e^{\frac{-Q_{\mathit{Ao}}}{b}}\right)\right]\mathrm{\Delta }\mathit{So}+\dots \\ \left(D_4{K}_{1o}{K}_{2_O}\right)\mathrm{\Delta }{X}_H+\left[\frac{a}{b}\left({\mathit{So}}_{\mathit{sat}}-{\mathit{So}}_o\right)e^{\frac{-Q_{\mathit{Ao}}}{b}}\right]\mathrm{\Delta }{Q}_A\end{array}$$E27

$$\begin{array}{l}\mathrm{\Delta }\stackrel{·}{X_H}=\left({\mu }_H{K}_{2o}{K}_3{X}_{H_o}\right)\mathrm{\Delta }\mathit{Ss}+\left({\mu }_H{K}_{1o}{K}_4{X}_{H_o}\right)\mathrm{\Delta }\mathit{So}\\ +\left({\mu }_H{K}_{1o}{K}_{2o}-D_5-b_H\right)\mathrm{\Delta }{X}_H\end{array}$$E28

Here ${\mathit{Ss}}_o,{\mathit{So}}_o,X_{\mathit{Ho}}$, Q_Ao, $K_{1o}$, and $K_{2o}$ 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: $\mathrm{\Delta }\mathit{Ss}=\mathit{Ss}-{\mathit{Ss}}_o$, $\mathrm{\Delta }\mathit{So}=\mathit{So}-{\mathit{So}}_o$ y $\mathrm{\Delta }{X}_H=X_H-X_{\mathit{Ho}}$.

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

$$\left[\begin{array}{c}\mathrm{\Delta }\stackrel{·}{\mathit{Ss}}\\ \mathrm{\Delta }\stackrel{·}{\mathit{So}}\\ \mathrm{\Delta }\stackrel{·}{X_H}\end{array}\right]=\left[\begin{array}{ccc}{\mathit{Ss}}_1& {\mathit{Ss}}_2& {\mathit{Ss}}_3\\ {\mathit{So}}_1& {\mathit{So}}_2& {\mathit{So}}_3\\ X_{H1}& X_{H2}& X_{H3}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }\mathit{Ss}\\ \mathrm{\Delta }\mathit{So}\\ \mathrm{\Delta }{X}_H\end{array}\right]+\left[\begin{array}{cccc}{D}_1& 0& 0& 0\\ 0& D_1& 0& {\mathit{So}}_4\\ 0& 0& D_1& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{\mathit{Ss}}_f\\ \mathrm{\Delta }{\mathit{So}}_f\\ \mathrm{\Delta }{X}_{\mathit{Hf}}\\ \mathrm{\Delta }{Q}_A\end{array}\right]$$E29

where,

$$\begin{array}{l}{\mathit{Ss}}_1=-\left(D+D_1{K}_{2o}{K}_3{X}_{H_o}\right)\\ {\mathit{Ss}}_2=-\left(D_1{K}_{1o}{K}_4{X}_{H_o}\right)\\ {\mathit{Ss}}_3=-\left(D_1{K}_{1o}{K}_{2o}\right)\\ {\mathit{So}}_1=\left(D_4{K}_{2o}{K}_3{X}_{H_o}\right)\\ {\mathit{So}}_2=\left[D_4{K}_{1o}{K}_4{X}_{H_o}-D_3-a\left(1-e^{\frac{-Q_{\mathit{Ao}}}{b}}\right)\right]\\ {\mathit{So}}_3=\left(D_4{K}_{1o}{K}_{2_O}\right)\\ {\mathit{So}}_4=\left[\frac{a}{b}\left({\mathit{So}}_{\mathit{sat}}-{\mathit{So}}_o\right)e^{\frac{-Q_{\mathit{Ao}}}{b}}\right]\\ X_{H1}=\left({\mu }_H{K}_{2o}{K}_3{X}_{H_o}\right)\\ X_{H2}=\left({\mu }_H{K}_{1o}{K}_4{X}_{H_o}\right)\\ X_{H3}=\left({\mu }_H{K}_{1o}{K}_{2o}-D_5-b_H\right)\end{array}$$E30a

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

$$\left[\begin{array}{c}\mathrm{\Delta }\stackrel{·}{\mathit{Ss}}\\ \mathrm{\Delta }\stackrel{·}{\mathit{So}}\\ \mathrm{\Delta }\stackrel{·}{X_H}\end{array}\right]=\left[\begin{array}{ccc}‐14.8622& ‐0.1117& ‐0.0028\\ ‐4.8220& ‐0.5977& ‐0.0009\\ 9.7902& 0.0748& ‐0.0856\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }\mathit{Ss}\\ \mathrm{\Delta }\mathit{So}\\ \mathrm{\Delta }{X}_H\end{array}\right]+\left[\begin{array}{cccc}0.25& 0& 0& 0\\ 0& 0.25& 0& 143.9339\\ 0& 0& 0.25& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{\mathit{Ss}}_f\\ \mathrm{\Delta }{\mathit{So}}_f\\ \mathrm{\Delta }{X}_{\mathit{Hf}}\\ \mathrm{\Delta }{Q}_A\end{array}\right]$$E30b

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):

$$x_{t+1}=f\left(x_t,u_t\right)$$E31

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\left(x_{1:T},u_{1:T}\right)=\sum _{t\in \left[T\right]}{g}_t\left(x_t,u_t\right)$$E32a

where,

$$x_{1:T}≔\left(x_1,\dots ,x_T\right);u_{1:T}≔\left(u_1,\dots ,u_T\right)$$

Restrictions codify the allowed domains for states and entries:

$$\begin{array}{l}{x}_{t+1}=f\left(x_t,u_t\right),\forall t\in \left[T-1\right]\\ x_t\in {\mathcal{X}}_t,u_t\in {\mathcal{U}}_t;\forall t\in \left[T\right]\\ x_1=x_{\mathit{\text{init}}}\end{array}$$E32b

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:

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 13.

General block diagram of a MBPC strategy.

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.

General scheme of a MBPC strategy.

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 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 ${\widehat{u}}_{1:T};$

  1. $u_{1:T}\leftarrow {\widehat{u}}_{1:T}$;

  2. While task not completed do

     $x_{\mathit{\text{init}}}\leftarrow$ Get current state ();

     $u_{1:T\leftarrow }$ Solve optimization problem

     $u\leftarrow \mathit{\text{first}}\left(u_{1:T}\right)$

     $\mathit{\text{ApplyInput}}\left(u\right)$ to the system

  3. $\mathit{\text{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:

$$\left\{\mathrm{y}\left(\mathrm{t}+\mathrm{k}/\mathrm{t}\right),\mathrm{k}=1,2,\dots ,\mathrm{N}2\right\}$$E33a

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

$$\left\{\mathrm{y}\left(\mathrm{t}\right),\mathrm{y}\left(\mathrm{t}-1\right)\dots ,\mathrm{u}\left(\mathrm{t}-1\right),\mathrm{u}\left(\mathrm{t}-2\right),\dots \right\}$$E33b

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:

$$\left\{\mathrm{u}\left(\mathrm{t}/\mathrm{t}\right),\mathrm{u}\left(\mathrm{t}+1\right)/\mathrm{t},\mathrm{u}\left(\mathrm{t}+2/\mathrm{t}\right)\dots \right\}$$E33c

• 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\left(t+k/t\right)=x\left(t+k/t\right)+n\left(t+k/t\right)$$

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:

Figure 15.

Parallel model prediction.

$$\begin{array}{l}{n}_f=\frac{D\left(q^{-1}\right)}{C\left(q^{-1}\right)}n\left(t\right)\\ =-c_1{n}_f\left(t-1\right)-c_2{n}_f\left(t-2\right)-\dots \\ \dots +n\left(t\right)+d_{1}n\left(t-1\right)+d_{2}n\left(t-2\right)+\dots \end{array}$$E34

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):

$$n\left(t\right)=y\left(t\right)-x\left(t\right)$$E35

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

$$n_f\left(t+k/t\right)\equiv 0,k=1\dots N_2$$E36

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

$$n_f\left(t\right)=\frac{D\left(q^{-1}\right)}{C\left(q^{-1}\right)}\frac{C\left(q^{-1}\right)}{D\left(q^{-1}\right)}e\left(t\right)=e\left(t\right)$$E37

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 ≠ 0). The forward prediction of n(t) is given by Eq. (38):

$$\begin{array}{l}\mathit{\text{for}}k=1:N_2\\ n\left(t+k/t\right)=\frac{C\left(q^{-1}\right)}{D\left(q^{-1}\right)}{n}_f\left(t+k/t\right)\end{array}$$E38

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

$$\begin{array}{c}n\left(t+k/t\right)=d_{1}n\left(t+k-1/t\right)-d_{2}n\left(t+k-2/t\right)-\dots \\ +n_f\left(t+k/t\right)+c_1{n}_f\left(t+k-1/t\right)\end{array}$$E39

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

$$y\left(t+k/t\right)=y_{\mathit{\text{base}}}\left(t+k/t\right)+y_{\mathit{\text{optimize}}}\left(t+k/t\right)$$E40

• ${\mathit{y}}_{\mathit{\text{base}}}\left(\mathit{t}+\mathit{k}/\mathit{t}\right):$ it takes into account the effects of past control actions, {u(t-1), u(t-2) …}, the effects of a basic future control scenario (u_base (t + k/t)), which is chosen for linear systems, and in nonlinear systems it is different, and finally the effect of future perturbations n(t + k/t).

• ${\mathit{y}}_{\mathit{\text{optimize}}}\left(\mathit{t}+\mathit{k}/\mathit{t}\right)$_: it takes into account the effect of optimization of future control actions, {$\mathit{\delta u}\left(t/t\right),\mathit{\delta u}\left(t+1/t\right),\dots \mathit{\delta u}\left(t+N_u-1/t\right)$}.

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.

Figure 16.

u, past, base, and optimum control actions.

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{array}{c}{y}_{\mathit{\text{optimize}}}\left(t+k/t\right)=h_k\mathit{\delta u}\left(t/t\right)+h_{k-1}\mathit{\delta u}\left(t+1/t\right)+\dots \\ \dots +g_{k-N_u+1}\mathit{\delta u}\left(t+N_u-1/t\right)\end{array}$$E41

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

$$\mathit{Y}=\overline{\mathit{Y}}+\mathit{G}.\mathit{U}$$E42

$$\left\{\begin{array}{c}\begin{array}{c}Y={\left[y\left(t+N_1/t\right)\dots y\left(t+N_2/t\right)\right]}^T\\ \overline{Y}={\left[y_{\mathit{\text{base}}}\left(t+N_1/t\right),\dots y_{\mathit{\text{base}}}\left(t+N_2/t\right)\right]}^T\end{array}\\ U={\left[\mathit{\delta u}\left(\frac{t}{t}\right)\dots \mathit{\delta u}\left(t+N_u-1/t\right)\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\left(t+k/t\right)-y\left(t+k/t\right)\right]}^2$$E43

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

Figure 17.

Objective of the predictive control strategy.

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

$$\begin{array}{c}\sum _{k=N_1}^{N_2}{\left[r\left(t+\frac{k}{t}\right)-y\left(t+\frac{k}{t}\right)\right]}^2={\left[\mathit{R}-\overline{\mathit{Y}}-\mathit{G}.\mathit{U}\right]}^T\left[\mathit{R}-\overline{\mathit{Y}}-\mathit{G}.\mathit{U}\right]\\ \end{array}$$E44

$$U^{\ast }={\left[G^{T}G\right]}^{-1}\left[G^T\left(\mathit{R}-\overline{\mathit{Y}}\right)\right]$$E45

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

$$u\left(t\right)=u_{\mathit{\text{base}}}\left(t/t\right)+\mathit{\delta u}\left(t/t\right)=u_{\mathit{\text{base}}}\left(t/t\right)+U^{\ast }\left(1\right)$$E46

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\left(t\right)=x\left(t\right)+n\left(t\right)=\frac{B\left(z^{-1}\right)}{A\left(z^{-1}\right)}u\left(t-1\right)+\frac{C\left(z^{-1}\right)}{\mathrm{\Delta }D\left(z^{-1}\right)}e\left(t\right)$$E47

where,

x(t) = process transfer function, n(t) = disturbance transfer function, Δ = integrator in the disturbance, $e\left(t\right)$ = white noise, and u(t-1) = control action (intrinsic delay in the process). With A(z⁻¹), B(z⁻¹), C(z⁻¹), and D(z⁻¹), polynomials are obtained by identification. Following the procedure described in (Aguado A), Eq. (47) is transformed into Eq. (48):

$$\widehat{y}\left(t+i/t\right)=G_i\mathrm{\Delta }u\left(t+i-1\right)+{\mathrm{\Gamma }}_i\mathrm{\Delta }{u}^f\left(t-1\right)+F_i{y}^f\left(t\right)$$E48

This expression allows us to know at instant t, the value of the predicted output at instant t + i. Here, G_i y${\mathrm{\Gamma }}_i$ y F_i are polynomials FIR in z⁻¹. Δu^f and y^f are filtered inputs and outputs, respectively. The vector expression of the prediction model is shown in Eq. (49):

$$Y=\mathit{Gu}+\mathrm{\Gamma \Delta }{U}^f+{\mathit{FY}}^f$$E49a

where,

$$\mathrm{\Delta }{U}^f={\left[\mathrm{\Delta }{u}^f\left(t-1\right)\dots \mathrm{\Delta }{u}^f\left(t-2\right)\dots \mathrm{\Delta }{u}^f\left(t-\mathit{nt}\right)\right]}^T$$E49b

$$Y^f={\left[y^f\left(t\right)\dots y^f\left(t-1\right)\dots y^f\left(t-\mathit{na}\right)\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\left(u\right)=\sum _{i=1}^N{\alpha }_i{\left[\widehat{y}\left(t+i/t\right)-w\left(t+i\right)\right]}^2+\sum _{j=1}^N{\lambda }_j{\left[\mathrm{\Delta }u\left(t+j-1\right)\right]}^2$$E50

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

$$J\left(u\right)={\left(Y-W\right)}^T\alpha \left(Y-W\right)+u^T\mathit{\lambda u}$$E51

$\alpha$ and $\lambda$ are diagonal matrices of NxN, with:

$$Y={\left[\widehat{y}\left(t+1/t\right)\dots \widehat{y}\left(t+2/t\right)\dots \widehat{y}\left(t+N/t\right)\right]}^T$$E52

$$u={\left[\mathrm{\Delta }u\left(t\right)\dots \mathrm{\Delta }u\left(t+1\right)\dots \mathrm{\Delta }u\left(t+N-1\right)\right]}^T$$E53

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

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):

$$u={\left(G^T\mathit{\alpha G}+\lambda \right)}^{-1}{G}^T\alpha \left(W-\mathrm{\Gamma \Delta }{U}^f-{\mathit{FY}}^f\right)$$E55

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

$$\mathrm{\Delta }u\left(t\right)=\mathit{hW}-h\mathrm{\Gamma \Delta }{U}^f-{\mathit{hFY}}^f$$E56

where h is the first row of ${\left(G^T\mathit{\alpha G}+\lambda \right)}^{-1}{G}^T\alpha$.

The Z-transform of $\mathrm{\Delta }u\left(t\right)$ allows us to obtain the configuration of the GPC controller, as shown in Figure 18, represented by Eq. (57):

Figure 18.

GPC controller (green).

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

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

$$\mathit{\text{KMPC}}=\left[0.00581561.9218\mathrm{e}‐006‐6.0503\mathrm{e}‐006‐1.209\mathrm{e}‐0071.8134\mathrm{e}‐0062.4445\mathrm{e}‐006\right]$$

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.

Figure 19.

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

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

Figure 20.

GPC controller, with disturbances.

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.

Figure 21.

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

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 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.

Conflict of interest

The authors declare no conflict of interest.