# Analytical Fuzzy Predictive Control Applied to Wastewater Treatment Biological Processes

Authors: Pedro M. Vallejo LLamas; Pastora Vega
Journal: Complexity
 (2019)

---

## Abstract

A novel controlfuzzy predictive control law is proposed and successfully applied to a wastewater treatment process in this paper. The proposed control law allows us to evaluate the control signal in an analytical way, each sampling time being a nonlinear and fuzzy alternative to other classic predictive controllers. The control law is based on the formalization of the internal fuzzy predictive model of the process as linear time-varying state space equations that are updated every discrete time instant to take into account the nonlinearity effects due to disturbance action and changes in the operating point with time. The model is then used to evaluate the predictions, and, taking them as a starting point and considering them as a paradigm of the predictive functional control strategy, a control law, it is derived in an analytical and explicit way by imposing on the outputs of the follow-up of certain reference trajectories previously established. The work presented here addresses the application of this particular strategy of intelligent predictive control to the case of an activated sludge wastewater treatment process successfully in a simulation environment of a real plant taking into account real data for the disturbance records. Such a process is multivariable, nonlinear, time varying, and difficult to control due to its biological nature. The proposed control law can be straightforwardly used within a dual-mode MPC scheme to handle constraints, as a nonlinear and fuzzy alternative to the classic state feedback control law.


---

## Body

## 1. Introduction

The traditional strategy ofmodel-based predictive control (MBPC or MPC) [1–3] consists basically of the use of a prediction model to determine the necessary control actions at each instant, imposing the minimization of a cost function (which will generally include a term dependent on the error and another dependent on the control efforts, among other possible scenarios). This strategy includes however multiple variants or approaches depending on various factors, mainly on the type of model that will be used to calculate the predictions and the mathematical algorithm used to determine the control law. There are various fundamental methods based on linear models of the process which determine the control variable by means of optimization (the minimization of the cost function chosen). Other methods consider nonlinear mathematical models (fuzzy models [4], models formalized by means of artificial neuronal networks [5], or other alternatives), but they also determine the control signal by means of optimization in the same way as the first ones. In [6], an orderly bibliographic review of the evolution of the linear and nonlinear MPC is carried out, detailing numerous contributions and approaches in this field of process control. Our study falls within the category of predictive control strategies based on nonlinear models and more specifically on fuzzy models: fuzzy model-based predictive control (FMBPC). For the internal model of the process, a fuzzy model of the Takagi-Sugeno (TS) [7] type was chosen in which the premises of the rules are diffuse logical expressions while the conclusions are linear numerical combinations of the consequents. TS-type fuzzy models are suitable for the identification and description of complex nonlinear processes from numerical input-output data and, if available, expert knowledge of the process. The parameters of our TS fuzzy model were deduced by identification as from input-output numerical data previously obtained in simulation. In [6, 8–15], several process control strategies can be seen, based on Takagi-Sugeno fuzzy models, some of them framed in the FMBPC field, with the identification made from available input-output data.Concerning the calculation of the control variable, in our study, such a calculation is carried out following a methodology that could be considered as an extension of the so-calledpredictive functional control (PFC) [2, 16–18] (which was initially designed for linear systems) to nonlinear systems. The approach that we have chosen uses a fuzzy model of the process (for the calculation of the predictions) developed in the state space form and follows a PFC strategy for the derivation of a control law in an analytical explicit form, which can be an advantage in comparison with optimization-based control schemes [10, 18]. In relation to this issue, it is possible to say that some strategies or approaches of MPC could be considered as analytical MPC algorithms, such as dynamic matrix control, generic model control, or predictor-corrector control (in case of unconstrained linear MPC problems with quadratic cost functions, the control rules can be expressed in an analytical form and even for nonlinear MPC problems, if the model is linearized). But in the case of these algorithms, the analytical obtaining of a control law implies the statement and solution of a problem of optimization of a cost function. However, the PFC strategy offers us a method for the obtaining of the control law in which it is not necessary to solve an optimization problem, which is mathematically less complex. In addition, our approach deals with nonlinear models in a direct and practical way, by formalizing the fuzzy model in the form of time-varying state space equations (with coefficients that must be updated at each discrete time instant) and without this supposing a great increment in the complexity of the mathematical derivation of the control law.The main contribution of our work to theFMBPC field is the application of this approach to a multivariable, with disturbances, strongly nonlinear, with a complex dynamic and of a biological nature case study (carried out by simulation), with the aim of a future generalization.No parameters have been included to handle constraints in the proposed control law, that is, it is an unconstrained control law. But our control law can be used, as base law, within aMBPC scheme that is designed to satisfy previously fixed constraints (on the control action, on the control action increment, on the plant output, etc.). In particular, following the approach of the MPC schemes proposed by Rossiter [3] could be used (OLP or CLP dual-mode MPC schemes).To obtain the control law by means of the procedure mentioned above, some relevant mathematical tasks must be carried out. It will be necessary to deduce the mathematical expression that relates the model outputs, for a certain prediction horizon, to the control variable, that is, the appropriate expression for calculating the predictions, making use of the fuzzy model of the process. In addition, it will also be necessary to specify the control objective for the outputs of the plant, which in our case will consist in the imposition of the tracking of the previously established reference trajectories. And, of course, we must also formalize mathematically themodel-based predictive control characteristic or nature, by means of some kind of relationship between the model and the process (the plant). We will establish such a relationship using the incremental equivalence principle between the model and the plant (equivalence between the plant output-desired increment and the model output increment) [2, 16, 17]. Finally, combining appropriately all the relationships and equations mentioned, it will be necessary to find the control variable that guarantees the control objective established at each sampling instant. Due to the multivariable character of the process considered and the intrinsic mathematical complexity of the fuzzy models, in order to obtain the control variable in an analytical and explicit manner, we need to start from a (fuzzy) mathematical model expressed compactly and clearly. In addition to working with matrix expressions, therefore, in this study, we have formalized the expressions of the fuzzy model with a format similar to that of state equations, with the peculiarity that the coefficients of the various terms are not constant but depend on the instant of sampling. This is due to the fact that they depend on the instantaneous premise vector (more specifically on the levels of compliance with the various rules by the premise vector at each instant); it is therefore necessary for these coefficients to be updated (recalculated) during each period. This particular strategy of predictive control considered in this paper has previously been approached by other authors; to be precise, it is developed in [10] for a case study with a manipulated input and a single controlled output and without considering disturbances (previously, something similar was also developed in [18], for another case study). This work however approaches the case of a multivariable system and disturbances in the input are considered. To be precise, our system has three inputs: a manipulated input and two disturbances. It also has two outputs, both of which are controlled. Its multivariable character involves working with matrices, which leads to an increase in mathematical complexity. Moreover, due to the existence of two disturbances and a single manipulated input, it is assumed that it will not be easy to control the process, especially taking into account that the process considered is of a biological nature (significantly more unpredictable than many industrial physical-chemical processes). Tackling all these complexities is another of the contributions of this paper.The present work has focused on the application of our fuzzy predictive control strategy to wastewater treatment biological processes. In our case study, the plant to be controlled was awastewater treatment plant (WWTP) that originally had a relatively simple architecture (something that is clearly advantageous for the purposes of the study) and whose depuration method was the well-known activated sludge process. The entire study was implemented through simulation.The control of the biological processes present in wastewater treatment plants is a rather complex problem, due to the nonlinear character of the corresponding physical-chemical reactions and also to the abrupt changes that can occur (sometimes unpredictably) in the influent input flow rate and in its degree of contamination (disturbances). The proposals that have been made to improve the control and operation of biological processes (generic or related toWWTPs) have been numerous and different, highlighting in our case those that include some variant of nonlinear model-based predictive control (NLMPC) and, especially, the proposals framed in the field of the fuzzy model-based predictive control (FMBPC). In [19–34], some of these contributions can be seen. Our work tries to be, precisely, a contribution more in this last area, proposing a strategy of predictive control based on a fuzzy model (obtained from numerical input-output data), but formalizing the law of control in an analytical and explicit way.The paper is organized as follows. Section2 describes the case study chosen: a municipal WWTP with a relatively simple architecture and with a biological purification method, the so-called activated sludge. The subsection related to the mathematical model is complemented by the Appendix D. In Section 3, the fuzzy modeling and identification procedure of the treatment plant is developed, starting only from input-output numerical data and also the subsequent formalization of the fuzzy model in an equivalent model in the state space (with time-dependent coefficients). In Section 4, we address the problem of control, in the FMBPC field, adopting a nonlinear PFC strategy based on a fuzzy model. In this section, all the mathematical tasks necessary to deduce an analytical and explicit control law for our case study are carried out (using Appendices A, B, and C). In Section 5, we explain the FMBPC scheme that will be implemented and we present and explain some experiments developed by simulation and their corresponding results. And, finally, Section 6 contains the conclusions of the work.


## 2. Case Study: Wastewater Treatment Processes

Our case study consists of applying the previously introduced theFMBPC strategy to a multivariable biological process, highly nonlinear and with complex dynamics, but with a simple architecture that can facilitate the extraction of useful conclusions. For this reason, a municipal wastewater treatment plant (WWTP) located in Manresa (in the province of Barcelona, Spain) [35], equipped with only the elements and processes necessary for the purification of basic organic water pollution, that is, for the purification of the substrate, was chosen. The study was carried out in a simulation environment, for obvious reasons of feasibility and availability.


### 2.1. Plant Description and Input-Output Configuration

The wastewater treatment plant chosen was originally designed to achieve a relatively simple purification process, without including a nitrification-denitrification process. So, the central objective of the overall process of purification in this plant was mainly the reduction of the concentration of the substrate, carried out through the so-calledactivated sludge process, of a biological nature. The plant had a single subprocess of substrate biological depuration (with six reactors working in parallel, arranged in two lines of three, which can be considered equivalent to one only for the purposes of mathematical modeling) and, following the reactors, a single decantation subprocess (with two settlers working in parallel, which can also be considered as one, for modeling purposes). The equivalent architecture of the considered plant corresponds to that shown in Figure 1 with a single aerated biological reactor followed by a secondary settler; the meaning of the main variables of the WWTP can be seen in Table 1.Figure 1
Wastewater treatment plant with activated sludge.Table 1
Main variables of the WWTP with activated sludge.
Influent
Biological reactor
Secondary settler and recirculation flows
Effluent
qi: input flow rate
q: reactor input-output flow rate
q2: total sludge recirculation flow rateqp: purge flow rateqr: sludge recirculation flow rate (to the reactor)
qout: output flow rate
si: input substrate concentrationxi: input biomass concentration
sir: substrate concentration in the input of the reactorxir: biomass concentration in the input of the reactors: output substrate concentrationx: output biomass concentrationc: output oxygen concentrationfk: aeration factor
xd: biomass concentration in the top layer of the settlerxb: biomass concentration in the middle layer of the settlerxr: biomass concentration in the lower layer of the settler
s: output substrate concentrationAs already indicated in the introductory chapter, the purification technique used in this plant is consisting of the elimination of organic pollutants by means ofactivated sludge (a mixed culture of microorganisms in suspension in the aerated biological reactor), with the recirculation of this sludge being the main control action and the sludge recirculation flow rate to the reactor (qr) being the only manipulated variable considered. The purification process that reduces the substrate concentration in the water is based on the interaction, through an aerobic reaction, between the microorganisms (the biomass) and the organic matter present in the water (the substrate). The microorganisms feed on the substrate (digest it), and consequently, its concentration decreases and therefore the contamination. This reaction takes place in the biological reactor (where the mix is supposed to be perfect) and requires a sufficient concentration of oxygen dissolved in the water. The concentration of oxygen (c) is one of the three variables involved in the basic mathematical models of the purification process and also needs to be controlled. However, in order to reduce the complexity of the case study, we will assume that the control of c will be guaranteed by the implementation of an adequate aeration system (a set of aeration turbines) and one independent control loop (with a PID algorithm for example), with an appropriate sampling period. After the biological reaction, the treated water is sent to a secondary settler, where the clean water and the activated sludge are separated. The clean water is sent to the outside of the plant and the activated sludge is recirculated and divided in two flows: one part of the sludge (qp) is purged from the bottom of the settler, and the other one (qr) is recirculated into the reactor with the objective of maintaining the population of the microorganisms.In our case study, we will control only the substrate and the biomass (s, x). The input variables are the input flow rate (qi), the input substrate concentration (si), and the sludge recirculation flow rate sent to the reactor (qr). The first two are the system disturbances and the third is the manipulated variable. Figure 2 shows the block diagram of the multivariable nonlinear system considered, with the input and output variables involved and their role in the control system (in the case of the inputs, an alternative notation has been used, between brackets, more suitable and usual in control systems).Figure 2
The multivariable nonlinear system considered.The information relating to the considered inputs and outputs of the plant is shown more detailed in Table2 including for each of them and also their biological or physical-chemical significance.Table 2
Inputs and outputs of the wastewater treatment plant.
Inputs
Disturbances
u1d1
Input flow rate (influent)
qi
u2d2
Substrate concentration in the influent
si
Manipulated variable
u3u
Sludge recirculation flow rate
qr
Outputs
Controlled variable 1
y1
Substrate concentration in the effluent
s
Controlled variable 2
y2
Biomass concentration in the reactor
xThe system chosen has two disturbances in the input, and the ultimate objective is to simultaneously control two output variables (coupled) by means of a single manipulated variable. In this study, we consider the wastewater treatment plant basically as a multivariable system case study with a complex dynamic and strongly nonlinear, the model of which (that could be unknown in a generic case) is identified from input-output numerical data.


### 2.2. Mathematical Model of the Wastewater Treatment Process

The mathematical model of the wastewater treatment process taken as reference is based on mass balances for the substrate, biomass, and oxygen, and it is founded on the classicalMonod and Maynard-Smith model (with the assumption of a perfectly mixed tank reactor). This model could be considered a simplification of the standard model denominated as activated sludge model no. 1, which is better known by its initials, ASM1 [36], but it is the corresponding model to the plant considered, a simple plant in its origin, but real.The equations that constitute the model and the values of the different parameters of the ManresaWWTP [35] can be seen in [23] and in Appendix D (in English language). We will use this model, in substitution or representation of the plant, both in the identification process and in the simulation study of the proposed control strategy. But we will only use the equations of the substrate and the biomass, taking into account the aforementioned, in the sense of assigning oxygen control to an independent loop. This simplification does not however make the study any less interesting because it approaches a general problem of interest: the possible utility of a predictive control scheme based on fuzzy models to control strongly nonlinear multivariable systems, starting only from the information implicit in input-output sets of data.


## 2.1. Plant Description and Input-Output Configuration

The wastewater treatment plant chosen was originally designed to achieve a relatively simple purification process, without including a nitrification-denitrification process. So, the central objective of the overall process of purification in this plant was mainly the reduction of the concentration of the substrate, carried out through the so-calledactivated sludge process, of a biological nature. The plant had a single subprocess of substrate biological depuration (with six reactors working in parallel, arranged in two lines of three, which can be considered equivalent to one only for the purposes of mathematical modeling) and, following the reactors, a single decantation subprocess (with two settlers working in parallel, which can also be considered as one, for modeling purposes). The equivalent architecture of the considered plant corresponds to that shown in Figure 1 with a single aerated biological reactor followed by a secondary settler; the meaning of the main variables of the WWTP can be seen in Table 1.Figure 1
Wastewater treatment plant with activated sludge.Table 1
Main variables of the WWTP with activated sludge.
Influent
Biological reactor
Secondary settler and recirculation flows
Effluent
qi: input flow rate
q: reactor input-output flow rate
q2: total sludge recirculation flow rateqp: purge flow rateqr: sludge recirculation flow rate (to the reactor)
qout: output flow rate
si: input substrate concentrationxi: input biomass concentration
sir: substrate concentration in the input of the reactorxir: biomass concentration in the input of the reactors: output substrate concentrationx: output biomass concentrationc: output oxygen concentrationfk: aeration factor
xd: biomass concentration in the top layer of the settlerxb: biomass concentration in the middle layer of the settlerxr: biomass concentration in the lower layer of the settler
s: output substrate concentrationAs already indicated in the introductory chapter, the purification technique used in this plant is consisting of the elimination of organic pollutants by means ofactivated sludge (a mixed culture of microorganisms in suspension in the aerated biological reactor), with the recirculation of this sludge being the main control action and the sludge recirculation flow rate to the reactor (qr) being the only manipulated variable considered. The purification process that reduces the substrate concentration in the water is based on the interaction, through an aerobic reaction, between the microorganisms (the biomass) and the organic matter present in the water (the substrate). The microorganisms feed on the substrate (digest it), and consequently, its concentration decreases and therefore the contamination. This reaction takes place in the biological reactor (where the mix is supposed to be perfect) and requires a sufficient concentration of oxygen dissolved in the water. The concentration of oxygen (c) is one of the three variables involved in the basic mathematical models of the purification process and also needs to be controlled. However, in order to reduce the complexity of the case study, we will assume that the control of c will be guaranteed by the implementation of an adequate aeration system (a set of aeration turbines) and one independent control loop (with a PID algorithm for example), with an appropriate sampling period. After the biological reaction, the treated water is sent to a secondary settler, where the clean water and the activated sludge are separated. The clean water is sent to the outside of the plant and the activated sludge is recirculated and divided in two flows: one part of the sludge (qp) is purged from the bottom of the settler, and the other one (qr) is recirculated into the reactor with the objective of maintaining the population of the microorganisms.In our case study, we will control only the substrate and the biomass (s, x). The input variables are the input flow rate (qi), the input substrate concentration (si), and the sludge recirculation flow rate sent to the reactor (qr). The first two are the system disturbances and the third is the manipulated variable. Figure 2 shows the block diagram of the multivariable nonlinear system considered, with the input and output variables involved and their role in the control system (in the case of the inputs, an alternative notation has been used, between brackets, more suitable and usual in control systems).Figure 2
The multivariable nonlinear system considered.The information relating to the considered inputs and outputs of the plant is shown more detailed in Table2 including for each of them and also their biological or physical-chemical significance.Table 2
Inputs and outputs of the wastewater treatment plant.
Inputs
Disturbances
u1d1
Input flow rate (influent)
qi
u2d2
Substrate concentration in the influent
si
Manipulated variable
u3u
Sludge recirculation flow rate
qr
Outputs
Controlled variable 1
y1
Substrate concentration in the effluent
s
Controlled variable 2
y2
Biomass concentration in the reactor
xThe system chosen has two disturbances in the input, and the ultimate objective is to simultaneously control two output variables (coupled) by means of a single manipulated variable. In this study, we consider the wastewater treatment plant basically as a multivariable system case study with a complex dynamic and strongly nonlinear, the model of which (that could be unknown in a generic case) is identified from input-output numerical data.


## 2.2. Mathematical Model of the Wastewater Treatment Process

The mathematical model of the wastewater treatment process taken as reference is based on mass balances for the substrate, biomass, and oxygen, and it is founded on the classicalMonod and Maynard-Smith model (with the assumption of a perfectly mixed tank reactor). This model could be considered a simplification of the standard model denominated as activated sludge model no. 1, which is better known by its initials, ASM1 [36], but it is the corresponding model to the plant considered, a simple plant in its origin, but real.The equations that constitute the model and the values of the different parameters of the ManresaWWTP [35] can be seen in [23] and in Appendix D (in English language). We will use this model, in substitution or representation of the plant, both in the identification process and in the simulation study of the proposed control strategy. But we will only use the equations of the substrate and the biomass, taking into account the aforementioned, in the sense of assigning oxygen control to an independent loop. This simplification does not however make the study any less interesting because it approaches a general problem of interest: the possible utility of a predictive control scheme based on fuzzy models to control strongly nonlinear multivariable systems, starting only from the information implicit in input-output sets of data.


## 3. Fuzzy Identification and State Space Modeling of a Wastewater Treatment Plant

One of the great advantages of using fuzzy models is their potential for describing nonlinear behaviour. This capacity depends on the suitability of the identification process. Fuzzy identification is therefore a data treatment procedure (with modelling objectives) of great significance both in the general sector of the modelling of nonlinear systems and in the more specific one of predictive control based on models. However, the identification of fuzzy models and their use within predictive control strategies can be studied and applied from many points of view. Our study focuses on the application to nonlinear multivariable systems of a particular fuzzy modelling methodology which consists of the following: expressing a fuzzy model (previously obtained by means of identification) in the form of state equations in such a way that it can be useful in calculating predictions and in the search for and obtaining of an explicit analytical expression for the predictive control law within the framework of afuzzy model-based predictive control (FMBPC) system [14]. It is for this reason that in this article, we will concentrate more on the analytical formalization of the fuzzy model that has already been identified than on the identification itself.In our case study, the wastewater treatment plant was identified starting from a series of input-output numerical data previously obtained by open-loop simulation and subsequently treated with the fuzzy identification software tool known asFMID (fuzzy model identification toolbox) [37]. Some functions of the toolbox were (partially) adapted and some code complements were also included (in the Matlab and Simulink environment) so as to be able to make our calculations and carry out our experiments. This tool, which is based on clustering techniques by means of the Gustafson-Kessel algorithm, was developed by Babuška et al. as software support for the theories and techniques of fuzzy modelling and identification described in the book Fuzzy Modelling for Control [9].


### 3.1. Input-Output Data

The input-output numerical data of our system were obtained by means of simulation in open loop with the wastewater treatment plant represented by the differential equations corresponding to its traditional nonlinear mathematical model ([23], Appendix D). The identification process could also have been carried out with samples from the real system, but in this case, the predictive control tests would have been carried out with the real wastewater treatment plant, which is not usually possible.The multivariable system to identify corresponding to our case study has 3 inputs and 2 outputs. Figure3 shows another simplified block diagram of the wastewater treatment plant to emphasize the multivariable nature of the said system with multiple inputs and outputs (MIMO system).Figure 3
The MIMO system to identify (3 inputs and 2 outputs).The inputs of the considered plant were the following: the incoming water flow rate (influent),qik; the substrate concentration in the influent (organic pollution of the incoming water), sik; and the activated sludge recirculation flow rate, qrk. The first two are disturbances and the third is an input variable that can be used as a control variable (manipulated variable). And the two outputs taken into account are the substrate concentration in the outgoing water (effluent), sk, and the biomass concentration (microorganisms that feed on the substrate and therefore purify the water) in the reactor, xk.The values chosen for the first two inputs (disturbances) were chosen by taking as a reference data from real experimental campaigns of an industrial wastewater treatment plant, to be precise, of the municipal wastewater treatment plant of Manresa, in the province of Barcelona, Spain (see the physical-chemical parameters of such a treatment plant in Table 12, in sub-AppendixD1). These data originated from tests and measures carried out at this plant on the occasion of the work on predictive control carried out at the time by Moreno [35].For the control variable, pseudorandom value sequences were chosen, initially taking as a reference the values of the control variable used in the aforementioned campaigns of the industrial wastewater treatment plant and subsequently incorporating modifications both in variability and in the extreme values, aiming to cover different areas of operation, all this with the aim of capturing and extracting sufficient information on the dynamics of the wastewater treatment plant. For each combination of values of the input variables (qik,sik,qrk), we determined by simulation the corresponding output values (sk,xk), thus obtaining (for each test) a matrix of input-output data with 5 columns and as many rows as samples (thousands of samples), with each line consisting of 5 values (three inputs and two outputs). The sampling period of the data set recorded in the campaign was of the order of one hour. But the data was extended by performing a mathematical interpolation between each two samples, obtaining intermediate estimated data every 0.2 hours (i.e., 12 minutes).Figures4 and 5 represent the input-output data of one of the numerous identification tests carried out, which we will refer to as case A. Figure 4 is a graphic representation of the sequences of values of the inputs, and Figure 5 is a graphic representation of the sequences of values of the outputs, the latter having been obtained by means of simulation in open loop (being the wastewater treatment plant represented by the aforementioned standard nonlinear mathematical model). The variable associated with the abscissa axis of both graphs is the sample number, not the time, but the samples were ordered temporarily, as is logical.Figure 4
Identification data of case A (inputs).Figure 5
Identification data of case A (outputs).


### 3.2. Fuzzy Model: Type, Structure, and Parameters

The internal model that we will use for to the predictions, in our model-based predictive control frame, will be aTakagi-Sugeno-type discrete time fuzzy model [7]. These models are composed of a series of if-then rules, each of which represents a linear submodel corresponding to a certain subset or partition of the universe of fuzzy values of the antecedent vector (vector premise). The antecedent or premise of each of the rules is composed of several simple propositions connected by means of and logical operators. The simple propositions compare each of the components of the antecedent vector with a certain set or associated fuzzy value (characterised by its membership function). And the consequent or conclusion of each rule allocates to the output a linear combination of the variables that form the consequent vector, plus an independent term. Algorithm 1 shows the general form of this type of discrete time fuzzy models. The mathematical expression by means of which the numeric value that will allocate to the output is calculated is a function of the consequent vector x and has been represented in the said algorithm by ϕjx.Algorithm 1: General form of the Takagi-Sugeno-type fuzzy models.
Rj: if (xa1isAj1andxa2isAj2and, …,andxap is Ajp) then:
yk=ϕjx
=αj1x1+αj2x2+⋯+αjqxq+δj
where:
j=1,2,…,mr; mr: number of rules
xa=xa1,xa2,…,xap: antecedent vector
x=x1,x2,⋯,xq: consequent vector
Aj1,Aj2,…,Ajp: fuzzy sets related to the components
of the antecedent vectorruleRj
Membership function associated to the fuzzy setAjh:
θAjh:ℝ⟶0,1 (smooth function)
xah⟶μAjhxah (membership grade of xah
with respect to the fuzzy setAjh)
h=1,…,pThe composition of both the antecedent vector and the consequent vector, as well as the number of rules, is two of the main aspects of the structure of a fuzzy TS model. The first is directly related to the dynamics of the process we want to model, that is, will depend on the input-output dynamic dependencies or relationship, and the number of rules will depend on the number of clusters observed or considered in theproduct space of the available input-output data set. Another characteristic aspect of the fuzzy model (less intuitive) is the type of membership functions associated with different sets or fuzzy values.Before carrying out the identification process, we must choose what structure we will consider, that is, what decision we will take on the characteristics mentioned above. Such a choice could be based on hypothesis or be the consequence of expert knowledge of the process to be identified (if such empirical information is available) and implies the need to set various concrete numerical parameters, such as the specific number of rules and the order of discrete time recursive models that make up each rule. Subsequently, after the appropriate identification process, based on our case in input-output numerical data, we will obtain the specific numerical values of the coefficients of the different terms of the mathematical expression of the consequent of each of the rules, as well as the centers of clusters and constants relating to the membership functions.


### 3.3. Analytical Expression of the Fuzzy Model Output

To calculate the outputs of the considered fuzzy models (given by several rules), a numerical calculation method should be applied from among those contemplated in the fuzzy logic theory, adding or combining all the rules, with the appropriate weight to the consequent of each rule; for example (in our case), by means of the centroid method:(1)y~i=∑j=1mriμAj1xa1μAj2xa2⋯μAjpixapiϕijx∑j=1mriμAj1xa1μAj2xa2⋯μAjpixapi,
(i)
i=1,2 (two outputs).
(ii)
j=1,2,…,mri.
(iii)
mri is the number of rules of the output yi.
(iv)
pi is the number of components of the antecedent vector corresponding to yi.
(v)
μAjpixapi is the membership grade of xapi with respect to the fuzzy set Ajpi.And defining for each output and rule the following membership functions (normalized) for all the antecedent vector(2)βijxa=μAj1xa1μAj2xa2⋯μAjpixapi∑j=1mriμAj1xa1μAj2xa2⋯μAjpixapi,the numerical expression of the output will remain as
(3)y~i=∑j=1mriβijxaϕijx,i=1,2,j=1,2,…,mri,where mri is the number of rules of the output yi.


### 3.4. Identification of the WWTP

The fuzzy identification process consists of searching, as from the input-output data map obtained, a fuzzy model to represent the behaviour of the system as faithfully as possible. As we said above, in order to do so, it is necessary to choose, previously, a certain type of fuzzy model with a specific structure and with determined dynamic characteristics and then to try to adjust unknown coefficients or parameters by means of an appropriate mathematical method. In our experiments, as has been mentioned in the introduction to this section, the whole of the identification process was carried out by using theFMID software tool [37], which is capable of extracting a fuzzy model of the Takagi-Sugeno type [7] for each of the outputs of the system to identify, as from the objective information implicit in the input-output data provided and according to the previous choice of certain parameters concerning the possible dynamics of the plant. Numerous identification tests were carried out with different input-output data set (corresponding to different campaigns) and with different choices for the values of various parameters (available in the tool) related to the dynamic characteristics of the recursive discrete model that will represent the plant. The result obtained in each of the tests carried out was a fuzzy model for each of the two outputs. Some of the experiments done are presented in [31, 32] and many others were subsequently carried out as the study and the research progressed. In this paper, we will refer to some of the case studied during the whole of the process, which we will refer to as case A, case B, case C, and case D.The main dynamic structural parameters that need to be chosen to carry out the identification of the WWTP with theFMID software tool can be summarized in Table 3.Table 3
Structural parameters for fuzzy identification.
Parameter
Meaning
Description
c
Number of data clusters (for each output)
Vector with two integer components
Ny
Output-output dynamic relationship
Matrix (2 × 2)
Nu
Input-output dynamic relationship
Matrix (2 × 3)
Nd
Input-output transport delays
Matrix (2 × 3)
Ts
Sample time
Real numberNext, in Table4, we will specify the choice of parameters carried out in case D and its concrete meaning, as well as the relationship with the dynamics of the model.Table 4
Structural parameters in case D.
Parameter
Value
Dynamic model
c
65
Out-y1k=sk: 6 clusters/6 rules
Out-y2k=xk: 5 clusters/5 rules
Ny
1101
[Row 1]:y1k depends on y1k−1 and y2k−1
[Row 2]:y2k depends on y2k−1
Nu
112101
[Row 1]:y1k depends on u1k−1, u2k−1, u3k−1 and u3k−2 (the parameter value of the first row and the third column, equal to 2, implies two consecutive terms of u3 (in discrete time))
[Row 2]:y2k depends on u1k−1 and u3k−1
Nd
111101
[Row 1]: they are 1 transport delay for the three inputs (u1, u2, and u3) with respect to y1k
[Row 2]: they are 1 transport delay for theu1 and u3 inputs with respect to y2k
Ts
0.2
Sample time (hours)The result of each identification test carried out consists of a fuzzy model of theTakagi-Sugeno type for each of the two outputs of the wastewater treatment plant. After the identification, all parameters of each fuzzy model are available in appropriate numerical structures and we did an adequate use of them for the automatic processing of information in our simulation study. The rules of the two fuzzy models that have been identified corresponding to case D of the WWTP identifications carried out can be seen as follows: in Algorithm 2, the fuzzy model corresponding to output 1 and in Algorithm 3, the fuzzy model corresponding to output 2. The relationship between the parameters contained in Table 4 (except Ts) and the structure of the rules of the fuzzy model of the two outputs is quite clear. Thus, the recursive dynamic models of y1k and y2k have the following structure:
(4)y1k=f1y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2,oralternatively,y1k+1=f1y1k,y2k,u1k,u2k,u3k,u3k−1,
(5)y2k=f2y2k−1,u1k−1,u3k−1,oralternatively,y2k+1=f2y2k,u1k,u3k.Algorithm 2: Takagi-Sugeno fuzzy model for substrate y1k=sk.
1. if (y1k−1isA11andy2k−1isA12andu1k−1isA13andu2k−1isA14
andu3k−1isA15andu3k−2isA16) then
y1k=8.74·10−1·y1k−1−2.72·10−3·y2k−1+3.96·10−3·u1k−1
+2.24·10−2·u2k−1+5.33·10−3·u3k−1−5.98·10−3·u3k−2−
6.47·10−1
2. if (y1k−1isA21andy2k−1isA22andu1k−1isA23andu2k−1isA24
andu3k−1isA25andu3k−2isA26) then
y1k=8.58·10−1·y1k−1−3.87·10−3·y2k−1+5.19·10−3·u1k−1
+2.31·10−2·u2k−1−3.89·10−3·u3k−1+3.69·10−3·u3k−2+
6.93·10−1
3. …
4. …
5. …
6. if (y1k−1isA61andy2k−1isA62andu1k−1isA63andu2k−1isA64
andu3k−1isA65andu3k−2isA66) then
y1k=8.70·10−1·y1k−1−6.54·10−3·y2k−1+1.33·10−2·u1k−1
+3.59·10−2·u2k−1+3.01·10−2·u3k−1−3.04·10−2·u3k−2−
1.25·101Algorithm 3: Takagi-Sugeno fuzzy model for biomass y2k=xk.
1. if (y2k−1isA∗11andu1k−1isA∗12andu3k−1isA∗13) then
y2k=9.58·10−1·y2k−1−2.22·10−2·u1k−1+1.33·10−2·u3k−1+
1.02·102
2. if (y2k−1isA∗21andu1k−1isA∗22andu3k−1isA∗23) then
y2k=1.01·100·y2k−1−4.10·10−2·u1k−1+3.54·10−2·u3k−1+
2.45·101
3. …
4. …
5. if (y2k−1isA∗51andu1k−1isA∗52andu3k−1isA∗53) then
y2k=9.35·10−1·y2k−1−5.22·10−2·u1k−1−3.96·10−2·u3k−1+
2.62·102In Tables5 and 6, the consequent parameters and the cluster centers, respectively, can be seen corresponding to output 1, and in Tables 7 and 8, the consequent parameters and the cluster centers, respectively, can be seen corresponding to output 2. All these numerical coefficients are results provided by the identification tool that uses the Gustafson-Kessel algorithm. The clustering techniques implicit in the tool are decisive when identifying the linear submodels.Table 5
Consequent parameters for substratey1k=sk.
Rule
y1k−1
y2k−1
u1k−1
u2k−1
u3k−1
u3k−2
Offset
1
8.74·10−1
−2.72·10−3
…
…
…
−5.98·10−3
…
2
8.58·10−1
−3.87·10−3
…
…
…
3.69·10−3
…
3
…
…
…
…
…
…
…
4
…
…
…
…
…
…
…
5
8.83·10−1
−3.15·10−3
…
…
…
5.72·10−3
…
6
8.70·10−1
−6.54·10−3
…
…
…
−3.04·10−2
…Table 6
Cluster centers for substratey1k=sk.
Rule
y1k−1
y2k−1
u1k−1
u2k−1
u3k−1
u3k−2
1
5.71·101
1.84·103
…
…
…
8.29·102
2
5.72·101
1.99·103
…
…
…
1.07·103
3
…
…
…
…
…
…
4
…
…
…
…
…
…
5
1.06·102
1.99·103
…
…
…
7.90·102
6
1.66·102
1.90·103
…
…
…
1.56·103Table 7
Consequent parameters for biomassy2k=xk.
Rule
y1k−1
u1k−1
u3k−1
Offset
1
9.58·10−1
−2.22·10−2
1.33·10−2
1.02·102
2
1.01·100
−4.10·10−2
3.54·10−2
2.45·101
3
…
…
…
…
4
…
…
…
…
5
9.35·10−1
−5.22·10−2
−3.96·10−2
2.62·102Table 8
Cluster centers for biomassy2k=xk.
Rule
y2k−1
u1k−1
u3k−1
1
1.80·103
…
5.68·102
2
1.82·103
…
4.81·102
3
…
…
…
4
…
…
…
5
2.12·103
…
9.96·102The identification tool also provides the necessary numerical information relating to the fuzzy sets or valuesAij, with which the components of the antecedent vectors are compared. Such information must be searched in the analytical expressions of the membership functions corresponding to such sets, which will be given in a parametric form. Using these parameters automatically and developing the suitable software, we can represent in graphic form the membership functions of each set. As an example, we show below a graph (Figure 6) with the membership functions corresponding to the six fuzzy sets (one for each rule) associated with the u3k−1=qrk−1 variable, which is one of the components of the antecedent vector of the rules of the output y1k, for the case A.Figure 6
Membership functions ofu3k−1=qrk−1.Observing the previous tables, we can see the composition of the antecedent and consequent vectors of each of the two outputs. The antecedent vector of output-1,xa∣out1, coincides with the consequent vector of output 1, x∣out1, and the same is true for the antecedent vector and consequent vector of output 2, xa∣out2 and x∣out2. Therefore, we will refer to the antecedent-consequent vector of output 1 and the antecedent-consequent vector of output 2. The composition of such vectors and the physical meaning of each component are specified below.
(6)xa=x∣out1=y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2,(7)xa=x∣out2=y2k−1,u1k−1,u3k−1,where
(8)y1k−1=sk−1;effluentsubstrateatk−1,y2k−1=xk−1;reactorbiomassatk−1,u1k−1=qik−1;inputflowrateatk−1,u2k−1=sik−1;inffluentsubstrateatk−1,u3k−1=qrk−1;sludgerecirculationatk−1,u3k−2=qrk−2;sludgerecirculationatk−2,k−N=k−N∗T,N=1,2,T=samplingperioddiscretetimesystem.As can be seen by comparing (6) and (7), initially, the antecedent-consequent vectors of both outputs do not coincide. However, in this study, we decided to treat the mathematical expressions of the fuzzy models jointly for both outputs (using matrices), defining a single mixed antecedent-consequent vector formed by the union of the two antecedent-consequent vectors (considering, naturally, null factors where appropriate in the necessary mathematical developments). Such vector, antecedent-consequent common to both outputs, is as follows:
(9)xa=x=y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2,where the involved variables are the same as those detailed in (8).


### 3.5. Validation of the Identified Fuzzy Models

An essential part of any identification process is the validation of the identified models. TheFMID software tool [37] uses a specific mathematical procedure of validation and provides the results both in graphic form, comparing the real process output with the model estimated output and, numerically, by calculating and giving the validation index known as VAF (VAF: percentile variance accounted for between two signals), which is among those habitually used for giving validity to the identified models. In our study, the real process output was obtained by means of simulation (with the wastewater treatment plant represented by the nonlinear model mentioned in Section 2.2 and detailed in Appendix D, given in the form of differential equations) and the model estimated output was obtained by applying to the identified fuzzy model the same inputs as to the wastewater treatment plant simulated.In our research process, numerous identifications were made with the same or different input-output data sets. Likewise, the use of the available input-output data was diverse, beginning with basic tests, using the same set of input-output data in the identification and validation of the model, and then carrying out more reliable tests in which a partition was made of the available input-output data, using a certain percentage of the data to identify the model and the rest to validate the model. On the other hand, and as we said at the beginning of Section3.4, identifications with different dynamic characteristics of the recursive discrete model that will represent the plant were made. The presentation, analysis, and discussion of the results of all the identifications made, together with the study of the relationship between the used data or the choice of parameters of the structure of the fuzzy model and the goodness of the identification, are beyond the scope of this article. However, it is necessary to consider at least the four identification cases mentioned at the beginning of Section 3.4, cases A, B, C, and D, commenting in this section the graphic and numerical results of the validation process of the corresponding fuzzy models identified. We summarize in Table 9 the main characteristics of the four cases.Table 9
Parameters and use of input-output data in 4 cases.
Case
Dynamical Parameters of the recursive discrete model
Identification in-out data versus validation in-out data
Ny
Nu
Nd
The same
Different (partition)
A
1101
012001
011001
●
B
1101
012001
011001
●
C
1101
112101
111101
●
D
1101
112101
111101
●For case A and case B, dependence ony1k and y2k outputs with respect to u1k−1=qik−1 input has not been considered. But it has been considered in case C and case D. We are interested in testing the effectiveness of the control algorithm, both in cases where the model is the result of more complete identification (C and D) and in cases where the identification is incomplete or less precise (A and B).


#### 3.5.1. Results of the Validation Process for Case A and Case B

The results of the validation corresponding to case A are shown in Figure7. This figure contains two graphic representations (one for each output of the wastewater treatment plant under study) and compares the values of the real process output with those of the model estimated output (for the same set of values of the inputs), including in addition the VAF index, which measures the goodness of the identification for each of the two partial identifications.Figure 7
Validation of the identified fuzzy model corresponding to case A for output 1,sk, and output 2, xk.In case B, we have the same structure of rules for the fuzzy model as in the case A and therefore the same components for the antecedent vector and the consequent vector, but with the difference that in this case, the numerical input-output data sequences used in the identification process were different from those used in the validation process (the available input-output data were partitioned). The results of the corresponding validation, for the two outputs of the wastewater treatment plant, can be seen in Figure8.Figure 8
Validation of the identified fuzzy model corresponding to case B for output 1,sk, and output 2, xk.In relation to the two identifications shown above, we can give a brief comparative analysis of both cases, taking into account the corresponding graphs and validation indexes. In case A (Figure7), the adjustment is quite precise for the two outputs, with very high VAF indexes (over 92%), while in case B (Figure 8), the adjustment is markedly lower with VAF indexes around 62% (output 1) and 42% (output 2). The difference is due to the fact that in the second case, the input-output data used to identify the model were not the same as the input-output data used to validate the model, a situation in which a more difficult adjustment can be expected. In case B, however, although the VAF indexes are lower, the response of the fuzzy model follows quite well the evolution of output 1 of the wastewater treatment plant and acceptably that of output 2, the tendency of which was at least detected by the identification. On the other hand, one of the potential complementary objectives of this study is precisely the assessing of the efficiency of fuzzy predictive control algorithms, including with imperfect models, as indeed we have checked in the control experiments (FMBPC) carried out with both fuzzy models and which we will subsequently discuss.


#### 3.5.2. Results of the Validation Process for Case C and Case D

In case C and case D, dependence ony1k and y2k outputs with respect to u1k−1=qik−1 input has been considered. Therefore, in principle, these two cases are more precise and realistic when it comes to identifying the model of the plant and one would expect a greater utility of the models within a strategy of predictive control based on models. The dynamic structure of the model considered for the identification of the plant is the same in both cases, as can be seen in Table 8, but they differ (as in the previous two cases) in that, in case D, the numerical input-output data sequences used in the identification process were different from those used in the validation process (the available input-output data were partitioned). The results of the validation, for the two outputs of the wastewater treatment plant, corresponding to case C and case D, can be seen in Figures 9 and 10, respectively.Figure 9
Validation of the identified fuzzy model corresponding to case C for output 1,sk, and output 2, xk.Figure 10
Validation of the identified fuzzy model corresponding to case D for output 1,sk, and output 2, xk.In relation to the other two identifications shown above, we can also give a brief comparative analysis of both cases, considering the corresponding graphs and validation indexes shown. In case C (Figure9), the adjustment is quite precise for the two outputs, with very high VAF indexes (around 99% for output 1 and of 87% for output 2). In case D (Figure 10), the adjustment is also quite precise for output 1 (95%) and quite acceptable for output 2 (around 66%). In addition, the responses of the two fuzzy models follow quite well the evolution of the corresponding outputs of the wastewater treatment plant (output 1 and output 2), that is, the tendency in the two outputs was detected very well by the identification. Finally, it is relevant to emphasize that the goodness of the validation of the model identified in case D is significantly better than that of case B. The difference between both cases is the consideration, in case D, of the dependence of the two outputs with respect to u1k−1=qik−1 input. In the simulation study presented in Section 5, four different cases of predictive control (FMBPC) applied to the WWTP have been considered, each of them corresponding to the use of one of the four identified models (case A, case B, case C, and case D). We are interested in analyzing the possible influence of the models on the effectiveness of the control algorithm, although a detailed study of this aspect is not possible to address it in this article.


### 3.6. State Space Modeling

By using suitable definitions and mathematical treatments, we can describe or formalize the rules of the fuzzy models of theTakagi-Sugeno type in the form of state equations. This involves the huge advantage of being able to treat the fuzzy models in an analytical manner and also to calculate the predictions analytically and ultimately being able to express the control law in an analytical and explicit manner. The application of this methodology to the case study which concerns us constitutes one of the main contributions of our study, which follows the line previously initiated by other authors [10]. The mathematical development necessary for our case study has previously been presented [31]. The discrete time fuzzy models in the state space were formalized jointly for the two outputs of the water treatment plant (with the necessary definitions), and the equations corresponding to them both were grouped together using matrices to obtain a single state equation and a single output equation. The result of the mathematical formalization was as follows:
(10)zmk+1=A¯mzmk+B¯muak+R¯m,(11)ymk=C¯mzmk,where
(i)
zmk is the extended-state vector, integrated by the outputs and disturbances of the WWTP, at the kth instant
(ii)
zmk+1 is the extended-state vector, at the k+1th instant
(iii)
ymk is the output of the model in the state space, integrated by the outputs of the WWTP, at the kth instant
(iv)
uakinput of the model in the state space, integrated by the manipulated variable of the WWTP, at the kth instant and the previous one
(v)
A¯m,B¯m,R¯m, and C¯m are the state matricesand the corresponding expressions are
(12)zmk=y1ky2kd1kd2k=y1ky2ku1ku2k=skxkqiksik,(13)uak=u3ku3k−1=ukuk−1=qrkqrk−1,(14)ymk=y1ky2k,and the state matrices are
(15)A¯m=∑j=1mrβjxaAmj,B¯m=∑j=1mrβjxaBmj,C¯m=∑j=1mrβj12xaCmj,R¯m=∑j=1mrβjxaRmj,where
(16)βjxa=β1jxa0000β2jxa0000100001,βj12xa=β1jxa00β2jxa,β26xa=0,mr=max⋅mr1,mr2,being that
(17)Amj=aj1aj2bj1bj20aj2∗bj1∗0001mr00001mr,Bmj=bj3bj4bj3∗00000,Cmj=10000100,Rmj=rjrj∗00,
(i)
aj1,aj2,bj1,bj2,bj3,bj4,rj are the coefficients of the antecedent vector and independent term in the jth rule of the output 1 fuzzy model of the WWTP.
(ii)
aj2∗,bj1∗,bj3∗,rj∗ are the coefficients of the antecedent vector and independent term in the jth rule of the output 2 fuzzy model of the WWTP (with a62∗=0, b61∗=0, b63∗=0, and r6∗=0).It is important and relevant to point out that the matrix coefficients of the state equations obtained,A¯m,B¯m,R¯m, and C¯m, depend on the antecedent vector xa (through βjxa) and therefore also depend on the kth instant, because xa depends on time. This assumes that concerning the necessary calculations to achieve the simulation, it will be necessary to update those coefficients in each iteration after having updated the antecedent vector xa and subsequently βjxa. And as a conclusion to this process of formalization of the fuzzy models identified in the state space, we can summarize by saying that the behaviour of our nonlinear multivariable system, initially identified by fuzzy models, was finally represented by a state equation system with a linear shape, but with coefficients depending on time. In [10], two theoretical references are mentioned on the association of systems with nonlinear dynamics, with variant linear systems over time (Leith, D.J. and Leithead, W.E., 1998 and 1999).


## 3.1. Input-Output Data

The input-output numerical data of our system were obtained by means of simulation in open loop with the wastewater treatment plant represented by the differential equations corresponding to its traditional nonlinear mathematical model ([23], Appendix D). The identification process could also have been carried out with samples from the real system, but in this case, the predictive control tests would have been carried out with the real wastewater treatment plant, which is not usually possible.The multivariable system to identify corresponding to our case study has 3 inputs and 2 outputs. Figure3 shows another simplified block diagram of the wastewater treatment plant to emphasize the multivariable nature of the said system with multiple inputs and outputs (MIMO system).Figure 3
The MIMO system to identify (3 inputs and 2 outputs).The inputs of the considered plant were the following: the incoming water flow rate (influent),qik; the substrate concentration in the influent (organic pollution of the incoming water), sik; and the activated sludge recirculation flow rate, qrk. The first two are disturbances and the third is an input variable that can be used as a control variable (manipulated variable). And the two outputs taken into account are the substrate concentration in the outgoing water (effluent), sk, and the biomass concentration (microorganisms that feed on the substrate and therefore purify the water) in the reactor, xk.The values chosen for the first two inputs (disturbances) were chosen by taking as a reference data from real experimental campaigns of an industrial wastewater treatment plant, to be precise, of the municipal wastewater treatment plant of Manresa, in the province of Barcelona, Spain (see the physical-chemical parameters of such a treatment plant in Table 12, in sub-AppendixD1). These data originated from tests and measures carried out at this plant on the occasion of the work on predictive control carried out at the time by Moreno [35].For the control variable, pseudorandom value sequences were chosen, initially taking as a reference the values of the control variable used in the aforementioned campaigns of the industrial wastewater treatment plant and subsequently incorporating modifications both in variability and in the extreme values, aiming to cover different areas of operation, all this with the aim of capturing and extracting sufficient information on the dynamics of the wastewater treatment plant. For each combination of values of the input variables (qik,sik,qrk), we determined by simulation the corresponding output values (sk,xk), thus obtaining (for each test) a matrix of input-output data with 5 columns and as many rows as samples (thousands of samples), with each line consisting of 5 values (three inputs and two outputs). The sampling period of the data set recorded in the campaign was of the order of one hour. But the data was extended by performing a mathematical interpolation between each two samples, obtaining intermediate estimated data every 0.2 hours (i.e., 12 minutes).Figures4 and 5 represent the input-output data of one of the numerous identification tests carried out, which we will refer to as case A. Figure 4 is a graphic representation of the sequences of values of the inputs, and Figure 5 is a graphic representation of the sequences of values of the outputs, the latter having been obtained by means of simulation in open loop (being the wastewater treatment plant represented by the aforementioned standard nonlinear mathematical model). The variable associated with the abscissa axis of both graphs is the sample number, not the time, but the samples were ordered temporarily, as is logical.Figure 4
Identification data of case A (inputs).Figure 5
Identification data of case A (outputs).


## 3.2. Fuzzy Model: Type, Structure, and Parameters

The internal model that we will use for to the predictions, in our model-based predictive control frame, will be aTakagi-Sugeno-type discrete time fuzzy model [7]. These models are composed of a series of if-then rules, each of which represents a linear submodel corresponding to a certain subset or partition of the universe of fuzzy values of the antecedent vector (vector premise). The antecedent or premise of each of the rules is composed of several simple propositions connected by means of and logical operators. The simple propositions compare each of the components of the antecedent vector with a certain set or associated fuzzy value (characterised by its membership function). And the consequent or conclusion of each rule allocates to the output a linear combination of the variables that form the consequent vector, plus an independent term. Algorithm 1 shows the general form of this type of discrete time fuzzy models. The mathematical expression by means of which the numeric value that will allocate to the output is calculated is a function of the consequent vector x and has been represented in the said algorithm by ϕjx.Algorithm 1: General form of the Takagi-Sugeno-type fuzzy models.
Rj: if (xa1isAj1andxa2isAj2and, …,andxap is Ajp) then:
yk=ϕjx
=αj1x1+αj2x2+⋯+αjqxq+δj
where:
j=1,2,…,mr; mr: number of rules
xa=xa1,xa2,…,xap: antecedent vector
x=x1,x2,⋯,xq: consequent vector
Aj1,Aj2,…,Ajp: fuzzy sets related to the components
of the antecedent vectorruleRj
Membership function associated to the fuzzy setAjh:
θAjh:ℝ⟶0,1 (smooth function)
xah⟶μAjhxah (membership grade of xah
with respect to the fuzzy setAjh)
h=1,…,pThe composition of both the antecedent vector and the consequent vector, as well as the number of rules, is two of the main aspects of the structure of a fuzzy TS model. The first is directly related to the dynamics of the process we want to model, that is, will depend on the input-output dynamic dependencies or relationship, and the number of rules will depend on the number of clusters observed or considered in theproduct space of the available input-output data set. Another characteristic aspect of the fuzzy model (less intuitive) is the type of membership functions associated with different sets or fuzzy values.Before carrying out the identification process, we must choose what structure we will consider, that is, what decision we will take on the characteristics mentioned above. Such a choice could be based on hypothesis or be the consequence of expert knowledge of the process to be identified (if such empirical information is available) and implies the need to set various concrete numerical parameters, such as the specific number of rules and the order of discrete time recursive models that make up each rule. Subsequently, after the appropriate identification process, based on our case in input-output numerical data, we will obtain the specific numerical values of the coefficients of the different terms of the mathematical expression of the consequent of each of the rules, as well as the centers of clusters and constants relating to the membership functions.


## 3.3. Analytical Expression of the Fuzzy Model Output

To calculate the outputs of the considered fuzzy models (given by several rules), a numerical calculation method should be applied from among those contemplated in the fuzzy logic theory, adding or combining all the rules, with the appropriate weight to the consequent of each rule; for example (in our case), by means of the centroid method:(1)y~i=∑j=1mriμAj1xa1μAj2xa2⋯μAjpixapiϕijx∑j=1mriμAj1xa1μAj2xa2⋯μAjpixapi,
(i)
i=1,2 (two outputs).
(ii)
j=1,2,…,mri.
(iii)
mri is the number of rules of the output yi.
(iv)
pi is the number of components of the antecedent vector corresponding to yi.
(v)
μAjpixapi is the membership grade of xapi with respect to the fuzzy set Ajpi.And defining for each output and rule the following membership functions (normalized) for all the antecedent vector(2)βijxa=μAj1xa1μAj2xa2⋯μAjpixapi∑j=1mriμAj1xa1μAj2xa2⋯μAjpixapi,the numerical expression of the output will remain as
(3)y~i=∑j=1mriβijxaϕijx,i=1,2,j=1,2,…,mri,where mri is the number of rules of the output yi.


## 3.4. Identification of the WWTP

The fuzzy identification process consists of searching, as from the input-output data map obtained, a fuzzy model to represent the behaviour of the system as faithfully as possible. As we said above, in order to do so, it is necessary to choose, previously, a certain type of fuzzy model with a specific structure and with determined dynamic characteristics and then to try to adjust unknown coefficients or parameters by means of an appropriate mathematical method. In our experiments, as has been mentioned in the introduction to this section, the whole of the identification process was carried out by using theFMID software tool [37], which is capable of extracting a fuzzy model of the Takagi-Sugeno type [7] for each of the outputs of the system to identify, as from the objective information implicit in the input-output data provided and according to the previous choice of certain parameters concerning the possible dynamics of the plant. Numerous identification tests were carried out with different input-output data set (corresponding to different campaigns) and with different choices for the values of various parameters (available in the tool) related to the dynamic characteristics of the recursive discrete model that will represent the plant. The result obtained in each of the tests carried out was a fuzzy model for each of the two outputs. Some of the experiments done are presented in [31, 32] and many others were subsequently carried out as the study and the research progressed. In this paper, we will refer to some of the case studied during the whole of the process, which we will refer to as case A, case B, case C, and case D.The main dynamic structural parameters that need to be chosen to carry out the identification of the WWTP with theFMID software tool can be summarized in Table 3.Table 3
Structural parameters for fuzzy identification.
Parameter
Meaning
Description
c
Number of data clusters (for each output)
Vector with two integer components
Ny
Output-output dynamic relationship
Matrix (2 × 2)
Nu
Input-output dynamic relationship
Matrix (2 × 3)
Nd
Input-output transport delays
Matrix (2 × 3)
Ts
Sample time
Real numberNext, in Table4, we will specify the choice of parameters carried out in case D and its concrete meaning, as well as the relationship with the dynamics of the model.Table 4
Structural parameters in case D.
Parameter
Value
Dynamic model
c
65
Out-y1k=sk: 6 clusters/6 rules
Out-y2k=xk: 5 clusters/5 rules
Ny
1101
[Row 1]:y1k depends on y1k−1 and y2k−1
[Row 2]:y2k depends on y2k−1
Nu
112101
[Row 1]:y1k depends on u1k−1, u2k−1, u3k−1 and u3k−2 (the parameter value of the first row and the third column, equal to 2, implies two consecutive terms of u3 (in discrete time))
[Row 2]:y2k depends on u1k−1 and u3k−1
Nd
111101
[Row 1]: they are 1 transport delay for the three inputs (u1, u2, and u3) with respect to y1k
[Row 2]: they are 1 transport delay for theu1 and u3 inputs with respect to y2k
Ts
0.2
Sample time (hours)The result of each identification test carried out consists of a fuzzy model of theTakagi-Sugeno type for each of the two outputs of the wastewater treatment plant. After the identification, all parameters of each fuzzy model are available in appropriate numerical structures and we did an adequate use of them for the automatic processing of information in our simulation study. The rules of the two fuzzy models that have been identified corresponding to case D of the WWTP identifications carried out can be seen as follows: in Algorithm 2, the fuzzy model corresponding to output 1 and in Algorithm 3, the fuzzy model corresponding to output 2. The relationship between the parameters contained in Table 4 (except Ts) and the structure of the rules of the fuzzy model of the two outputs is quite clear. Thus, the recursive dynamic models of y1k and y2k have the following structure:
(4)y1k=f1y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2,oralternatively,y1k+1=f1y1k,y2k,u1k,u2k,u3k,u3k−1,
(5)y2k=f2y2k−1,u1k−1,u3k−1,oralternatively,y2k+1=f2y2k,u1k,u3k.Algorithm 2: Takagi-Sugeno fuzzy model for substrate y1k=sk.
1. if (y1k−1isA11andy2k−1isA12andu1k−1isA13andu2k−1isA14
andu3k−1isA15andu3k−2isA16) then
y1k=8.74·10−1·y1k−1−2.72·10−3·y2k−1+3.96·10−3·u1k−1
+2.24·10−2·u2k−1+5.33·10−3·u3k−1−5.98·10−3·u3k−2−
6.47·10−1
2. if (y1k−1isA21andy2k−1isA22andu1k−1isA23andu2k−1isA24
andu3k−1isA25andu3k−2isA26) then
y1k=8.58·10−1·y1k−1−3.87·10−3·y2k−1+5.19·10−3·u1k−1
+2.31·10−2·u2k−1−3.89·10−3·u3k−1+3.69·10−3·u3k−2+
6.93·10−1
3. …
4. …
5. …
6. if (y1k−1isA61andy2k−1isA62andu1k−1isA63andu2k−1isA64
andu3k−1isA65andu3k−2isA66) then
y1k=8.70·10−1·y1k−1−6.54·10−3·y2k−1+1.33·10−2·u1k−1
+3.59·10−2·u2k−1+3.01·10−2·u3k−1−3.04·10−2·u3k−2−
1.25·101Algorithm 3: Takagi-Sugeno fuzzy model for biomass y2k=xk.
1. if (y2k−1isA∗11andu1k−1isA∗12andu3k−1isA∗13) then
y2k=9.58·10−1·y2k−1−2.22·10−2·u1k−1+1.33·10−2·u3k−1+
1.02·102
2. if (y2k−1isA∗21andu1k−1isA∗22andu3k−1isA∗23) then
y2k=1.01·100·y2k−1−4.10·10−2·u1k−1+3.54·10−2·u3k−1+
2.45·101
3. …
4. …
5. if (y2k−1isA∗51andu1k−1isA∗52andu3k−1isA∗53) then
y2k=9.35·10−1·y2k−1−5.22·10−2·u1k−1−3.96·10−2·u3k−1+
2.62·102In Tables5 and 6, the consequent parameters and the cluster centers, respectively, can be seen corresponding to output 1, and in Tables 7 and 8, the consequent parameters and the cluster centers, respectively, can be seen corresponding to output 2. All these numerical coefficients are results provided by the identification tool that uses the Gustafson-Kessel algorithm. The clustering techniques implicit in the tool are decisive when identifying the linear submodels.Table 5
Consequent parameters for substratey1k=sk.
Rule
y1k−1
y2k−1
u1k−1
u2k−1
u3k−1
u3k−2
Offset
1
8.74·10−1
−2.72·10−3
…
…
…
−5.98·10−3
…
2
8.58·10−1
−3.87·10−3
…
…
…
3.69·10−3
…
3
…
…
…
…
…
…
…
4
…
…
…
…
…
…
…
5
8.83·10−1
−3.15·10−3
…
…
…
5.72·10−3
…
6
8.70·10−1
−6.54·10−3
…
…
…
−3.04·10−2
…Table 6
Cluster centers for substratey1k=sk.
Rule
y1k−1
y2k−1
u1k−1
u2k−1
u3k−1
u3k−2
1
5.71·101
1.84·103
…
…
…
8.29·102
2
5.72·101
1.99·103
…
…
…
1.07·103
3
…
…
…
…
…
…
4
…
…
…
…
…
…
5
1.06·102
1.99·103
…
…
…
7.90·102
6
1.66·102
1.90·103
…
…
…
1.56·103Table 7
Consequent parameters for biomassy2k=xk.
Rule
y1k−1
u1k−1
u3k−1
Offset
1
9.58·10−1
−2.22·10−2
1.33·10−2
1.02·102
2
1.01·100
−4.10·10−2
3.54·10−2
2.45·101
3
…
…
…
…
4
…
…
…
…
5
9.35·10−1
−5.22·10−2
−3.96·10−2
2.62·102Table 8
Cluster centers for biomassy2k=xk.
Rule
y2k−1
u1k−1
u3k−1
1
1.80·103
…
5.68·102
2
1.82·103
…
4.81·102
3
…
…
…
4
…
…
…
5
2.12·103
…
9.96·102The identification tool also provides the necessary numerical information relating to the fuzzy sets or valuesAij, with which the components of the antecedent vectors are compared. Such information must be searched in the analytical expressions of the membership functions corresponding to such sets, which will be given in a parametric form. Using these parameters automatically and developing the suitable software, we can represent in graphic form the membership functions of each set. As an example, we show below a graph (Figure 6) with the membership functions corresponding to the six fuzzy sets (one for each rule) associated with the u3k−1=qrk−1 variable, which is one of the components of the antecedent vector of the rules of the output y1k, for the case A.Figure 6
Membership functions ofu3k−1=qrk−1.Observing the previous tables, we can see the composition of the antecedent and consequent vectors of each of the two outputs. The antecedent vector of output-1,xa∣out1, coincides with the consequent vector of output 1, x∣out1, and the same is true for the antecedent vector and consequent vector of output 2, xa∣out2 and x∣out2. Therefore, we will refer to the antecedent-consequent vector of output 1 and the antecedent-consequent vector of output 2. The composition of such vectors and the physical meaning of each component are specified below.
(6)xa=x∣out1=y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2,(7)xa=x∣out2=y2k−1,u1k−1,u3k−1,where
(8)y1k−1=sk−1;effluentsubstrateatk−1,y2k−1=xk−1;reactorbiomassatk−1,u1k−1=qik−1;inputflowrateatk−1,u2k−1=sik−1;inffluentsubstrateatk−1,u3k−1=qrk−1;sludgerecirculationatk−1,u3k−2=qrk−2;sludgerecirculationatk−2,k−N=k−N∗T,N=1,2,T=samplingperioddiscretetimesystem.As can be seen by comparing (6) and (7), initially, the antecedent-consequent vectors of both outputs do not coincide. However, in this study, we decided to treat the mathematical expressions of the fuzzy models jointly for both outputs (using matrices), defining a single mixed antecedent-consequent vector formed by the union of the two antecedent-consequent vectors (considering, naturally, null factors where appropriate in the necessary mathematical developments). Such vector, antecedent-consequent common to both outputs, is as follows:
(9)xa=x=y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2,where the involved variables are the same as those detailed in (8).


## 3.5. Validation of the Identified Fuzzy Models

An essential part of any identification process is the validation of the identified models. TheFMID software tool [37] uses a specific mathematical procedure of validation and provides the results both in graphic form, comparing the real process output with the model estimated output and, numerically, by calculating and giving the validation index known as VAF (VAF: percentile variance accounted for between two signals), which is among those habitually used for giving validity to the identified models. In our study, the real process output was obtained by means of simulation (with the wastewater treatment plant represented by the nonlinear model mentioned in Section 2.2 and detailed in Appendix D, given in the form of differential equations) and the model estimated output was obtained by applying to the identified fuzzy model the same inputs as to the wastewater treatment plant simulated.In our research process, numerous identifications were made with the same or different input-output data sets. Likewise, the use of the available input-output data was diverse, beginning with basic tests, using the same set of input-output data in the identification and validation of the model, and then carrying out more reliable tests in which a partition was made of the available input-output data, using a certain percentage of the data to identify the model and the rest to validate the model. On the other hand, and as we said at the beginning of Section3.4, identifications with different dynamic characteristics of the recursive discrete model that will represent the plant were made. The presentation, analysis, and discussion of the results of all the identifications made, together with the study of the relationship between the used data or the choice of parameters of the structure of the fuzzy model and the goodness of the identification, are beyond the scope of this article. However, it is necessary to consider at least the four identification cases mentioned at the beginning of Section 3.4, cases A, B, C, and D, commenting in this section the graphic and numerical results of the validation process of the corresponding fuzzy models identified. We summarize in Table 9 the main characteristics of the four cases.Table 9
Parameters and use of input-output data in 4 cases.
Case
Dynamical Parameters of the recursive discrete model
Identification in-out data versus validation in-out data
Ny
Nu
Nd
The same
Different (partition)
A
1101
012001
011001
●
B
1101
012001
011001
●
C
1101
112101
111101
●
D
1101
112101
111101
●For case A and case B, dependence ony1k and y2k outputs with respect to u1k−1=qik−1 input has not been considered. But it has been considered in case C and case D. We are interested in testing the effectiveness of the control algorithm, both in cases where the model is the result of more complete identification (C and D) and in cases where the identification is incomplete or less precise (A and B).


### 3.5.1. Results of the Validation Process for Case A and Case B

The results of the validation corresponding to case A are shown in Figure7. This figure contains two graphic representations (one for each output of the wastewater treatment plant under study) and compares the values of the real process output with those of the model estimated output (for the same set of values of the inputs), including in addition the VAF index, which measures the goodness of the identification for each of the two partial identifications.Figure 7
Validation of the identified fuzzy model corresponding to case A for output 1,sk, and output 2, xk.In case B, we have the same structure of rules for the fuzzy model as in the case A and therefore the same components for the antecedent vector and the consequent vector, but with the difference that in this case, the numerical input-output data sequences used in the identification process were different from those used in the validation process (the available input-output data were partitioned). The results of the corresponding validation, for the two outputs of the wastewater treatment plant, can be seen in Figure8.Figure 8
Validation of the identified fuzzy model corresponding to case B for output 1,sk, and output 2, xk.In relation to the two identifications shown above, we can give a brief comparative analysis of both cases, taking into account the corresponding graphs and validation indexes. In case A (Figure7), the adjustment is quite precise for the two outputs, with very high VAF indexes (over 92%), while in case B (Figure 8), the adjustment is markedly lower with VAF indexes around 62% (output 1) and 42% (output 2). The difference is due to the fact that in the second case, the input-output data used to identify the model were not the same as the input-output data used to validate the model, a situation in which a more difficult adjustment can be expected. In case B, however, although the VAF indexes are lower, the response of the fuzzy model follows quite well the evolution of output 1 of the wastewater treatment plant and acceptably that of output 2, the tendency of which was at least detected by the identification. On the other hand, one of the potential complementary objectives of this study is precisely the assessing of the efficiency of fuzzy predictive control algorithms, including with imperfect models, as indeed we have checked in the control experiments (FMBPC) carried out with both fuzzy models and which we will subsequently discuss.


### 3.5.2. Results of the Validation Process for Case C and Case D

In case C and case D, dependence ony1k and y2k outputs with respect to u1k−1=qik−1 input has been considered. Therefore, in principle, these two cases are more precise and realistic when it comes to identifying the model of the plant and one would expect a greater utility of the models within a strategy of predictive control based on models. The dynamic structure of the model considered for the identification of the plant is the same in both cases, as can be seen in Table 8, but they differ (as in the previous two cases) in that, in case D, the numerical input-output data sequences used in the identification process were different from those used in the validation process (the available input-output data were partitioned). The results of the validation, for the two outputs of the wastewater treatment plant, corresponding to case C and case D, can be seen in Figures 9 and 10, respectively.Figure 9
Validation of the identified fuzzy model corresponding to case C for output 1,sk, and output 2, xk.Figure 10
Validation of the identified fuzzy model corresponding to case D for output 1,sk, and output 2, xk.In relation to the other two identifications shown above, we can also give a brief comparative analysis of both cases, considering the corresponding graphs and validation indexes shown. In case C (Figure9), the adjustment is quite precise for the two outputs, with very high VAF indexes (around 99% for output 1 and of 87% for output 2). In case D (Figure 10), the adjustment is also quite precise for output 1 (95%) and quite acceptable for output 2 (around 66%). In addition, the responses of the two fuzzy models follow quite well the evolution of the corresponding outputs of the wastewater treatment plant (output 1 and output 2), that is, the tendency in the two outputs was detected very well by the identification. Finally, it is relevant to emphasize that the goodness of the validation of the model identified in case D is significantly better than that of case B. The difference between both cases is the consideration, in case D, of the dependence of the two outputs with respect to u1k−1=qik−1 input. In the simulation study presented in Section 5, four different cases of predictive control (FMBPC) applied to the WWTP have been considered, each of them corresponding to the use of one of the four identified models (case A, case B, case C, and case D). We are interested in analyzing the possible influence of the models on the effectiveness of the control algorithm, although a detailed study of this aspect is not possible to address it in this article.


## 3.5.1. Results of the Validation Process for Case A and Case B

The results of the validation corresponding to case A are shown in Figure7. This figure contains two graphic representations (one for each output of the wastewater treatment plant under study) and compares the values of the real process output with those of the model estimated output (for the same set of values of the inputs), including in addition the VAF index, which measures the goodness of the identification for each of the two partial identifications.Figure 7
Validation of the identified fuzzy model corresponding to case A for output 1,sk, and output 2, xk.In case B, we have the same structure of rules for the fuzzy model as in the case A and therefore the same components for the antecedent vector and the consequent vector, but with the difference that in this case, the numerical input-output data sequences used in the identification process were different from those used in the validation process (the available input-output data were partitioned). The results of the corresponding validation, for the two outputs of the wastewater treatment plant, can be seen in Figure8.Figure 8
Validation of the identified fuzzy model corresponding to case B for output 1,sk, and output 2, xk.In relation to the two identifications shown above, we can give a brief comparative analysis of both cases, taking into account the corresponding graphs and validation indexes. In case A (Figure7), the adjustment is quite precise for the two outputs, with very high VAF indexes (over 92%), while in case B (Figure 8), the adjustment is markedly lower with VAF indexes around 62% (output 1) and 42% (output 2). The difference is due to the fact that in the second case, the input-output data used to identify the model were not the same as the input-output data used to validate the model, a situation in which a more difficult adjustment can be expected. In case B, however, although the VAF indexes are lower, the response of the fuzzy model follows quite well the evolution of output 1 of the wastewater treatment plant and acceptably that of output 2, the tendency of which was at least detected by the identification. On the other hand, one of the potential complementary objectives of this study is precisely the assessing of the efficiency of fuzzy predictive control algorithms, including with imperfect models, as indeed we have checked in the control experiments (FMBPC) carried out with both fuzzy models and which we will subsequently discuss.


## 3.5.2. Results of the Validation Process for Case C and Case D

In case C and case D, dependence ony1k and y2k outputs with respect to u1k−1=qik−1 input has been considered. Therefore, in principle, these two cases are more precise and realistic when it comes to identifying the model of the plant and one would expect a greater utility of the models within a strategy of predictive control based on models. The dynamic structure of the model considered for the identification of the plant is the same in both cases, as can be seen in Table 8, but they differ (as in the previous two cases) in that, in case D, the numerical input-output data sequences used in the identification process were different from those used in the validation process (the available input-output data were partitioned). The results of the validation, for the two outputs of the wastewater treatment plant, corresponding to case C and case D, can be seen in Figures 9 and 10, respectively.Figure 9
Validation of the identified fuzzy model corresponding to case C for output 1,sk, and output 2, xk.Figure 10
Validation of the identified fuzzy model corresponding to case D for output 1,sk, and output 2, xk.In relation to the other two identifications shown above, we can also give a brief comparative analysis of both cases, considering the corresponding graphs and validation indexes shown. In case C (Figure9), the adjustment is quite precise for the two outputs, with very high VAF indexes (around 99% for output 1 and of 87% for output 2). In case D (Figure 10), the adjustment is also quite precise for output 1 (95%) and quite acceptable for output 2 (around 66%). In addition, the responses of the two fuzzy models follow quite well the evolution of the corresponding outputs of the wastewater treatment plant (output 1 and output 2), that is, the tendency in the two outputs was detected very well by the identification. Finally, it is relevant to emphasize that the goodness of the validation of the model identified in case D is significantly better than that of case B. The difference between both cases is the consideration, in case D, of the dependence of the two outputs with respect to u1k−1=qik−1 input. In the simulation study presented in Section 5, four different cases of predictive control (FMBPC) applied to the WWTP have been considered, each of them corresponding to the use of one of the four identified models (case A, case B, case C, and case D). We are interested in analyzing the possible influence of the models on the effectiveness of the control algorithm, although a detailed study of this aspect is not possible to address it in this article.


## 3.6. State Space Modeling

By using suitable definitions and mathematical treatments, we can describe or formalize the rules of the fuzzy models of theTakagi-Sugeno type in the form of state equations. This involves the huge advantage of being able to treat the fuzzy models in an analytical manner and also to calculate the predictions analytically and ultimately being able to express the control law in an analytical and explicit manner. The application of this methodology to the case study which concerns us constitutes one of the main contributions of our study, which follows the line previously initiated by other authors [10]. The mathematical development necessary for our case study has previously been presented [31]. The discrete time fuzzy models in the state space were formalized jointly for the two outputs of the water treatment plant (with the necessary definitions), and the equations corresponding to them both were grouped together using matrices to obtain a single state equation and a single output equation. The result of the mathematical formalization was as follows:
(10)zmk+1=A¯mzmk+B¯muak+R¯m,(11)ymk=C¯mzmk,where
(i)
zmk is the extended-state vector, integrated by the outputs and disturbances of the WWTP, at the kth instant
(ii)
zmk+1 is the extended-state vector, at the k+1th instant
(iii)
ymk is the output of the model in the state space, integrated by the outputs of the WWTP, at the kth instant
(iv)
uakinput of the model in the state space, integrated by the manipulated variable of the WWTP, at the kth instant and the previous one
(v)
A¯m,B¯m,R¯m, and C¯m are the state matricesand the corresponding expressions are
(12)zmk=y1ky2kd1kd2k=y1ky2ku1ku2k=skxkqiksik,(13)uak=u3ku3k−1=ukuk−1=qrkqrk−1,(14)ymk=y1ky2k,and the state matrices are
(15)A¯m=∑j=1mrβjxaAmj,B¯m=∑j=1mrβjxaBmj,C¯m=∑j=1mrβj12xaCmj,R¯m=∑j=1mrβjxaRmj,where
(16)βjxa=β1jxa0000β2jxa0000100001,βj12xa=β1jxa00β2jxa,β26xa=0,mr=max⋅mr1,mr2,being that
(17)Amj=aj1aj2bj1bj20aj2∗bj1∗0001mr00001mr,Bmj=bj3bj4bj3∗00000,Cmj=10000100,Rmj=rjrj∗00,
(i)
aj1,aj2,bj1,bj2,bj3,bj4,rj are the coefficients of the antecedent vector and independent term in the jth rule of the output 1 fuzzy model of the WWTP.
(ii)
aj2∗,bj1∗,bj3∗,rj∗ are the coefficients of the antecedent vector and independent term in the jth rule of the output 2 fuzzy model of the WWTP (with a62∗=0, b61∗=0, b63∗=0, and r6∗=0).It is important and relevant to point out that the matrix coefficients of the state equations obtained,A¯m,B¯m,R¯m, and C¯m, depend on the antecedent vector xa (through βjxa) and therefore also depend on the kth instant, because xa depends on time. This assumes that concerning the necessary calculations to achieve the simulation, it will be necessary to update those coefficients in each iteration after having updated the antecedent vector xa and subsequently βjxa. And as a conclusion to this process of formalization of the fuzzy models identified in the state space, we can summarize by saying that the behaviour of our nonlinear multivariable system, initially identified by fuzzy models, was finally represented by a state equation system with a linear shape, but with coefficients depending on time. In [10], two theoretical references are mentioned on the association of systems with nonlinear dynamics, with variant linear systems over time (Leith, D.J. and Leithead, W.E., 1998 and 1999).


## 4. The Control Law in Analytical Form

The next step of ourfuzzy predictive control strategy consists of deducing an analytical and explicit control law, making use of the model in the form of state equations. As was mentioned in the introduction, the deduction method for the control law can be considered an extension of the so-called predictive functional control 1 [2, 16–18] for the multivariable case. To be able to use this method, the following will be necessary: on the one hand, to define both the desired behaviour of the closed-loop system and the control goal and, on the other hand, to formalize the relationship between the model and the plant that allows carrying out a model-based predictive control strategy. The desired behaviour will be defined by imposing on the plant outputs the follow-up of certain reference trajectories (one for each output), as faithfully as possible, and the control goal will consist on calculating the future control action so that the predicted plant output values coincide with the reference trajectory in at least one point. In our case, we will consider a single coincidence point. The difference between the current time (kth time instant) and the one corresponding to the coincidence point, measured on the number of sampling periods, is called coincidence horizon and it is usually represented by H. Therefore, starting from time instant k, the coincidence with each trajectory must occur at time instant k+H.The relationship between the model and the plant will be established assuming the so-calledequivalence principle between the plant and the model, used in PFC. This principle consists in supposing that, for each output of the plant, the plant output increment, between k and k+H, necessary to coincide with the reference trajectory, ∆p, will be equal to the predicted model output increment, also between k and k+H, ∆m. This principle is analytically expressed by the following expression:
(18)∆p=∆m,which must be satisfied for each output of the plant. If we consider the equations of both outputs simultaneously, we would have a vector equality: ∆p=∆m, where ∆p is usually the so-called objective increment vector and ∆m is the model output increment vector. The idea of the PFC strategy and equivalence principle can be shown in graphical form (Figure 11).Figure 11
Predictive functional control-equivalence principle.


### 4.1. Predictions Based on the Model

The control actions in the various existingpredictive control strategies are not calculated in the same way in all cases; there are different methods and different formulas or algorithms that are applied. However, all these strategies share the use of some type of mathematical model to represent as faithfully as possible the process to be controlled (in our case, a fuzzy model formalized in the state space) and which serves to predict the future behaviour of the process and ultimately to determine the appropriate control action for compliance with the control objectives. We can therefore say that the mathematical basis of the predictive control strategies will be the general expression corresponding to ymk+H, that is, the expression of the output at the instant k+H, predicted by the model at the kth instant, with H being the prediction horizon.The deduction of such an expression, from the state equations previously obtained, is fairly extensive in our case and has for this reason been developed in its totality in AppendixA with the objective of facilitating the follow-up of the article. The calculation of Hstep ahead prediction is done under the assumption of constant future manipulated variables, that is, uk=uk+1=⋯=uk+H−1. The final expression deduced is the following (see (A.26), at the end of Appendix A):
(19)ymk+H=C¯mA¯mHzmk+A¯mH−1B¯m+A¯mH−1−IA¯m−I−1B¯mP1010uak++A¯mH−1+A¯mH−1−IA¯m−I−1R¯m,where
(20)H∈Z+,H≥1;∃A¯mH−nifH−n≥1withn∈Z+,I=1001,theorder2identitymatrix,andP1010=1010.


### 4.2. Deduction of the Control Law in an Analytical and Explicit Way

As from the expression ofymk+H, which allows the prediction at the kth instant of the evolution of the outputs of the model, H periods further on, we can search an analytical and explicit expression for the control law. This is precisely the main objective of our paper. We detail below the necessary mathematical process. In the first place, we will find uak in (19), obtaining an expression that will be a function of ymk+H and we will subsequently extract, of the vector variable uak obtained, the scalar variable uk, which will also be a function of ymk+H. This first stage is developed in Section 4.2.1. Secondly, we will establish the desired behaviour of our system in a closed loop, imposing the follow-up by the outputs of the WWTP of certain reference trajectories, which will also condition to ymk+H. And as from the relations obtained, we will deduct the expression of the control action uk which is necessary to comply with the control objective, which will be a function of the reference (set point) of the future outputs (at k+H). The second stage is developed in Section 4.2.2.


#### 4.2.1. Manipulated Variableuk as a Function of ymk+H

Leaving at (19) to one side of equality, the term in which uak is multiplying
(21)ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m=C¯mA¯mH−1B¯m+A¯mH−1−IAm−I−1B¯mP1010uak,and grouping together terms for simplification, defining matrix Ma as follows:
(22)Ma=C¯mA¯mH−1Bm+A¯mH−1−IA¯m−I−1B¯mP1010,we will have (21) simplified, obtaining
(23)ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m=Mauak,and multiplying on the left, on both sides of equality, by the matrix Ma−1, we will have the following matrix expression:
(24)uak=Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,and replacing uak by its generic expression specified in (12), we will have
(25)ukuk−1=Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,where uk−1 will be the value of the variable u at the previous sampling instant to the kth and should therefore have been conveniently memorized (or for the first sampling instant have a predetermined initial value). In short, uk−1 will be a data, a specific numerical value saved, while uk will be the unknown quantity to find or determine.If matrix (25), with uk−1 given, has a single numerical solution (or if it has none, assuming that it is possible to determine an approximate solution by computer), then we will be able to determine the numerical value of uk by means of the following expression:
(26)uk=P10uak,being
(27)P10=10,and replacing uak by its expression given by (24), we would finally have
(28)uk=P10Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,which is the expression that will allow the determining of the variable uk which will be necessary at the kth sampling instant so that the model can attain, H periods later (i.e., at the k+Hth instant), a certain vector valueymk+H.


#### 4.2.2. Control Actionuk Necessary for the Follow-Up of the Output Reference Trajectories

The control actionuk which we search will be that guaranteeing compliance with the control objectives for the process output. Thus, following the same approach as in [10, 18], we will establish the desired behaviour of the closed-loop system imposing for the outputs the follow-up of discrete reference trajectories, which must gradually approach to the corresponding references. These trajectories constitute the reference model and will be given formally by means of the following recursive equations:
(29)yrik+1=ariyrik+briyset_pointik,i=1,2numberofoutputs,where yset_pointik is the reference or desired value of the ith output yik at the kth instant, yrik and yrik+1, the values in k and in k+1, respectively, of the reference trajectory, and ari and bri the parameters of the reference model, the values of which must be such so as to ensure that the gain of the reference model is the unit (for both outputs). For this to occur, as can be seen in Appendix B expression (B.5), the following mathematical relationship should be complied with
(30)1−ari−1bri=1,i=1,2numberofoutputs.And as we have introduced at the beginning of the section, to achieve the follow-up of the reference trajectories, the control goal fixed will consist on to calculate the future control action so that the predicted plant output values coincide with the reference trajectory (for each output of the plant) in a singlecoincidence point. Such a point will correspond to the instant k+H, where H is called the coincidence horizon.The control action at eachkth instant, uk, should therefore be such that the plant outputs match the corresponding reference trajectories, H sampling periods later, such that yik+H equals yrik+H. Therefore, the desired plant output increment should be equal to yrik+H−yik, for each ith output (i=1,2), that is,
(31)Δpi=yrik+H−yik,i=1,2numberofoutputs.On the other hand, we need now to do use of the relationship between the model and the plant established in the introduction of this section, consisting of adopting the so-calledequivalence principle (following again the same approach as in [10, 18]), which is a way to introduce the main idea or mechanism implicit in the model-based predictive control. We must therefore consider equality (18), ∆p=∆m, for each ith output (i=1,2), where ∆m is the predicted model output increment, whose expression is as follows:
(32)Δmi=ymik+H−ymik,i=1,2numberofoutputs,and matching (31) and (32) (i.e., matching the desired plant output increments with the predicted model output increments: ∆p∣i=∆m∣i), we will have
(33)yrik+H−yik=ymik+H−ymik,i=1,2numberofoutputs.The control action at eachkth instant, uk, should be such that equality (33) is satisfied for each of the two outputs. And considering together the two equalities implicit in (33), we will have the following single vector expression:
(34)yrk+H−yk=ymk+H−ymk.Finding nowymk+H in (34), we will have that the predicted model output in the k+1th time instant should satisfy the following expression:
(35)ymk+H=yrk+H−yk+ymk.Replacing (28) ymk+H by the second member of equality (35), which we have just obtained, we will finally have the expression that we were looking for uk, dependent on terms that can be determined using the plant model or the reference model or that can be measured at the kth instant. The final expression for uk is the following, where Ma is given by (22) and P10 by (27) and where the matrix coefficients that intervene must be updated in each iteration of the implementation of the control algorithm (due to their dependence on the antecedent vector and therefore on time)
(36)uk=P10Ma−1yrk+H−yk+ymk−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯mH∈Z+,H≥1,where the term yrk+H would remain to be specified, which would be developed by induction. This development is shown in an abbreviated form together with the result obtained under the hypothesis of maintenance of the reference on the prediction horizon in Appendix C. The result is as follows (see (C.10) at the end of said appendix):
(37)yrk+H=ArHyrk+I−ArHyset_pointk,ArH=ar1H00ar2H,I=1001,H∈Z+,H≥1.Formula (36), complemented with (37), constitutes the analytical and explicit expression of the control law uk which would be needed, at each kth instant, to guarantee the desired behaviour of our multivariable closed-loop system, which is precisely what we were searching in this section.


## 4.1. Predictions Based on the Model

The control actions in the various existingpredictive control strategies are not calculated in the same way in all cases; there are different methods and different formulas or algorithms that are applied. However, all these strategies share the use of some type of mathematical model to represent as faithfully as possible the process to be controlled (in our case, a fuzzy model formalized in the state space) and which serves to predict the future behaviour of the process and ultimately to determine the appropriate control action for compliance with the control objectives. We can therefore say that the mathematical basis of the predictive control strategies will be the general expression corresponding to ymk+H, that is, the expression of the output at the instant k+H, predicted by the model at the kth instant, with H being the prediction horizon.The deduction of such an expression, from the state equations previously obtained, is fairly extensive in our case and has for this reason been developed in its totality in AppendixA with the objective of facilitating the follow-up of the article. The calculation of Hstep ahead prediction is done under the assumption of constant future manipulated variables, that is, uk=uk+1=⋯=uk+H−1. The final expression deduced is the following (see (A.26), at the end of Appendix A):
(19)ymk+H=C¯mA¯mHzmk+A¯mH−1B¯m+A¯mH−1−IA¯m−I−1B¯mP1010uak++A¯mH−1+A¯mH−1−IA¯m−I−1R¯m,where
(20)H∈Z+,H≥1;∃A¯mH−nifH−n≥1withn∈Z+,I=1001,theorder2identitymatrix,andP1010=1010.


## 4.2. Deduction of the Control Law in an Analytical and Explicit Way

As from the expression ofymk+H, which allows the prediction at the kth instant of the evolution of the outputs of the model, H periods further on, we can search an analytical and explicit expression for the control law. This is precisely the main objective of our paper. We detail below the necessary mathematical process. In the first place, we will find uak in (19), obtaining an expression that will be a function of ymk+H and we will subsequently extract, of the vector variable uak obtained, the scalar variable uk, which will also be a function of ymk+H. This first stage is developed in Section 4.2.1. Secondly, we will establish the desired behaviour of our system in a closed loop, imposing the follow-up by the outputs of the WWTP of certain reference trajectories, which will also condition to ymk+H. And as from the relations obtained, we will deduct the expression of the control action uk which is necessary to comply with the control objective, which will be a function of the reference (set point) of the future outputs (at k+H). The second stage is developed in Section 4.2.2.


### 4.2.1. Manipulated Variableuk as a Function of ymk+H

Leaving at (19) to one side of equality, the term in which uak is multiplying
(21)ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m=C¯mA¯mH−1B¯m+A¯mH−1−IAm−I−1B¯mP1010uak,and grouping together terms for simplification, defining matrix Ma as follows:
(22)Ma=C¯mA¯mH−1Bm+A¯mH−1−IA¯m−I−1B¯mP1010,we will have (21) simplified, obtaining
(23)ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m=Mauak,and multiplying on the left, on both sides of equality, by the matrix Ma−1, we will have the following matrix expression:
(24)uak=Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,and replacing uak by its generic expression specified in (12), we will have
(25)ukuk−1=Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,where uk−1 will be the value of the variable u at the previous sampling instant to the kth and should therefore have been conveniently memorized (or for the first sampling instant have a predetermined initial value). In short, uk−1 will be a data, a specific numerical value saved, while uk will be the unknown quantity to find or determine.If matrix (25), with uk−1 given, has a single numerical solution (or if it has none, assuming that it is possible to determine an approximate solution by computer), then we will be able to determine the numerical value of uk by means of the following expression:
(26)uk=P10uak,being
(27)P10=10,and replacing uak by its expression given by (24), we would finally have
(28)uk=P10Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,which is the expression that will allow the determining of the variable uk which will be necessary at the kth sampling instant so that the model can attain, H periods later (i.e., at the k+Hth instant), a certain vector valueymk+H.


### 4.2.2. Control Actionuk Necessary for the Follow-Up of the Output Reference Trajectories

The control actionuk which we search will be that guaranteeing compliance with the control objectives for the process output. Thus, following the same approach as in [10, 18], we will establish the desired behaviour of the closed-loop system imposing for the outputs the follow-up of discrete reference trajectories, which must gradually approach to the corresponding references. These trajectories constitute the reference model and will be given formally by means of the following recursive equations:
(29)yrik+1=ariyrik+briyset_pointik,i=1,2numberofoutputs,where yset_pointik is the reference or desired value of the ith output yik at the kth instant, yrik and yrik+1, the values in k and in k+1, respectively, of the reference trajectory, and ari and bri the parameters of the reference model, the values of which must be such so as to ensure that the gain of the reference model is the unit (for both outputs). For this to occur, as can be seen in Appendix B expression (B.5), the following mathematical relationship should be complied with
(30)1−ari−1bri=1,i=1,2numberofoutputs.And as we have introduced at the beginning of the section, to achieve the follow-up of the reference trajectories, the control goal fixed will consist on to calculate the future control action so that the predicted plant output values coincide with the reference trajectory (for each output of the plant) in a singlecoincidence point. Such a point will correspond to the instant k+H, where H is called the coincidence horizon.The control action at eachkth instant, uk, should therefore be such that the plant outputs match the corresponding reference trajectories, H sampling periods later, such that yik+H equals yrik+H. Therefore, the desired plant output increment should be equal to yrik+H−yik, for each ith output (i=1,2), that is,
(31)Δpi=yrik+H−yik,i=1,2numberofoutputs.On the other hand, we need now to do use of the relationship between the model and the plant established in the introduction of this section, consisting of adopting the so-calledequivalence principle (following again the same approach as in [10, 18]), which is a way to introduce the main idea or mechanism implicit in the model-based predictive control. We must therefore consider equality (18), ∆p=∆m, for each ith output (i=1,2), where ∆m is the predicted model output increment, whose expression is as follows:
(32)Δmi=ymik+H−ymik,i=1,2numberofoutputs,and matching (31) and (32) (i.e., matching the desired plant output increments with the predicted model output increments: ∆p∣i=∆m∣i), we will have
(33)yrik+H−yik=ymik+H−ymik,i=1,2numberofoutputs.The control action at eachkth instant, uk, should be such that equality (33) is satisfied for each of the two outputs. And considering together the two equalities implicit in (33), we will have the following single vector expression:
(34)yrk+H−yk=ymk+H−ymk.Finding nowymk+H in (34), we will have that the predicted model output in the k+1th time instant should satisfy the following expression:
(35)ymk+H=yrk+H−yk+ymk.Replacing (28) ymk+H by the second member of equality (35), which we have just obtained, we will finally have the expression that we were looking for uk, dependent on terms that can be determined using the plant model or the reference model or that can be measured at the kth instant. The final expression for uk is the following, where Ma is given by (22) and P10 by (27) and where the matrix coefficients that intervene must be updated in each iteration of the implementation of the control algorithm (due to their dependence on the antecedent vector and therefore on time)
(36)uk=P10Ma−1yrk+H−yk+ymk−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯mH∈Z+,H≥1,where the term yrk+H would remain to be specified, which would be developed by induction. This development is shown in an abbreviated form together with the result obtained under the hypothesis of maintenance of the reference on the prediction horizon in Appendix C. The result is as follows (see (C.10) at the end of said appendix):
(37)yrk+H=ArHyrk+I−ArHyset_pointk,ArH=ar1H00ar2H,I=1001,H∈Z+,H≥1.Formula (36), complemented with (37), constitutes the analytical and explicit expression of the control law uk which would be needed, at each kth instant, to guarantee the desired behaviour of our multivariable closed-loop system, which is precisely what we were searching in this section.


## 4.2.1. Manipulated Variableuk as a Function of ymk+H

Leaving at (19) to one side of equality, the term in which uak is multiplying
(21)ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m=C¯mA¯mH−1B¯m+A¯mH−1−IAm−I−1B¯mP1010uak,and grouping together terms for simplification, defining matrix Ma as follows:
(22)Ma=C¯mA¯mH−1Bm+A¯mH−1−IA¯m−I−1B¯mP1010,we will have (21) simplified, obtaining
(23)ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m=Mauak,and multiplying on the left, on both sides of equality, by the matrix Ma−1, we will have the following matrix expression:
(24)uak=Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,and replacing uak by its generic expression specified in (12), we will have
(25)ukuk−1=Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,where uk−1 will be the value of the variable u at the previous sampling instant to the kth and should therefore have been conveniently memorized (or for the first sampling instant have a predetermined initial value). In short, uk−1 will be a data, a specific numerical value saved, while uk will be the unknown quantity to find or determine.If matrix (25), with uk−1 given, has a single numerical solution (or if it has none, assuming that it is possible to determine an approximate solution by computer), then we will be able to determine the numerical value of uk by means of the following expression:
(26)uk=P10uak,being
(27)P10=10,and replacing uak by its expression given by (24), we would finally have
(28)uk=P10Ma−1ymk+H−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯m,which is the expression that will allow the determining of the variable uk which will be necessary at the kth sampling instant so that the model can attain, H periods later (i.e., at the k+Hth instant), a certain vector valueymk+H.


## 4.2.2. Control Actionuk Necessary for the Follow-Up of the Output Reference Trajectories

The control actionuk which we search will be that guaranteeing compliance with the control objectives for the process output. Thus, following the same approach as in [10, 18], we will establish the desired behaviour of the closed-loop system imposing for the outputs the follow-up of discrete reference trajectories, which must gradually approach to the corresponding references. These trajectories constitute the reference model and will be given formally by means of the following recursive equations:
(29)yrik+1=ariyrik+briyset_pointik,i=1,2numberofoutputs,where yset_pointik is the reference or desired value of the ith output yik at the kth instant, yrik and yrik+1, the values in k and in k+1, respectively, of the reference trajectory, and ari and bri the parameters of the reference model, the values of which must be such so as to ensure that the gain of the reference model is the unit (for both outputs). For this to occur, as can be seen in Appendix B expression (B.5), the following mathematical relationship should be complied with
(30)1−ari−1bri=1,i=1,2numberofoutputs.And as we have introduced at the beginning of the section, to achieve the follow-up of the reference trajectories, the control goal fixed will consist on to calculate the future control action so that the predicted plant output values coincide with the reference trajectory (for each output of the plant) in a singlecoincidence point. Such a point will correspond to the instant k+H, where H is called the coincidence horizon.The control action at eachkth instant, uk, should therefore be such that the plant outputs match the corresponding reference trajectories, H sampling periods later, such that yik+H equals yrik+H. Therefore, the desired plant output increment should be equal to yrik+H−yik, for each ith output (i=1,2), that is,
(31)Δpi=yrik+H−yik,i=1,2numberofoutputs.On the other hand, we need now to do use of the relationship between the model and the plant established in the introduction of this section, consisting of adopting the so-calledequivalence principle (following again the same approach as in [10, 18]), which is a way to introduce the main idea or mechanism implicit in the model-based predictive control. We must therefore consider equality (18), ∆p=∆m, for each ith output (i=1,2), where ∆m is the predicted model output increment, whose expression is as follows:
(32)Δmi=ymik+H−ymik,i=1,2numberofoutputs,and matching (31) and (32) (i.e., matching the desired plant output increments with the predicted model output increments: ∆p∣i=∆m∣i), we will have
(33)yrik+H−yik=ymik+H−ymik,i=1,2numberofoutputs.The control action at eachkth instant, uk, should be such that equality (33) is satisfied for each of the two outputs. And considering together the two equalities implicit in (33), we will have the following single vector expression:
(34)yrk+H−yk=ymk+H−ymk.Finding nowymk+H in (34), we will have that the predicted model output in the k+1th time instant should satisfy the following expression:
(35)ymk+H=yrk+H−yk+ymk.Replacing (28) ymk+H by the second member of equality (35), which we have just obtained, we will finally have the expression that we were looking for uk, dependent on terms that can be determined using the plant model or the reference model or that can be measured at the kth instant. The final expression for uk is the following, where Ma is given by (22) and P10 by (27) and where the matrix coefficients that intervene must be updated in each iteration of the implementation of the control algorithm (due to their dependence on the antecedent vector and therefore on time)
(36)uk=P10Ma−1yrk+H−yk+ymk−C¯mA¯mHzmk−C¯mA¯mH−1+A¯mH−1−IA¯m−I−1R¯mH∈Z+,H≥1,where the term yrk+H would remain to be specified, which would be developed by induction. This development is shown in an abbreviated form together with the result obtained under the hypothesis of maintenance of the reference on the prediction horizon in Appendix C. The result is as follows (see (C.10) at the end of said appendix):
(37)yrk+H=ArHyrk+I−ArHyset_pointk,ArH=ar1H00ar2H,I=1001,H∈Z+,H≥1.Formula (36), complemented with (37), constitutes the analytical and explicit expression of the control law uk which would be needed, at each kth instant, to guarantee the desired behaviour of our multivariable closed-loop system, which is precisely what we were searching in this section.


## 5. FMBPC Applied to WWTP

The main objective of this work is the deduction of afuzzy predictive control (FMBPC) law, expressed in an analytical and explicit way and to apply it to the control of wastewater treatment biological processes that are multivariable, highly nonlinear, time-varying, and complex processes. These processes are difficult to control due to their biological nature, and, in addition, our control system will use a single manipulated variable to control two output variables. Therefore, the possible obtaining of positive results would be quite interesting.Next, we show the configuration of the implemented control system and the experiments carried out by means of simulation.


### 5.1. FMBPC Scheme

The implemented control strategy is in the field ofnonlinear predictive control (NLMPC) and uses principles and methodology of functional predictive control (PFC) in the deduction of the control law. The basic mechanisms of operation of this strategy have been reflected in Figure 12.Figure 12
FMBPC scheme with an internal fuzzy model.As can be seen, ourFMBPC scheme uses an internal fuzzy model to carry out the predictions and uses PFC principles for the calculation of the control law. The controller needs the information related to the set point of the outputs and the corresponding reference trajectories. In addition, in each iteration (sampling period), state matrices must be updated, because they depend on the premise or antecedent vector (and ultimately also of time). Therefore, the controller must also have the instantaneous information of all the variables that constitute the antecedent vector.The task sequence that must be executed for the implementation of ourFMBPC strategy and that constitute the control algorithm is shown in Algorithm 4.Algorithm 4: Control algorithm for FMBPC implementation.
Repeat at eachkth instant (sampling period, Ts):
Step 1. Update the current premise vector. Measurement of
all the necessary variables for the determination of
the current premise or antecedent vector,xa:
y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2
Step 2. Update the state matrices. The TS fuzzy model was
formalized in the form of linear time-varying state space
equations, and the state matrices depend on the current
premise vector. Therefore, at each instant, the following
parameters must be updated, using (14):
A¯m,B¯m,R¯m, and C¯m
Step 3. Calculate the predicted model outputs. The use of a
model to predict future behaviour is essential in MPC.
By means of the appropriate mathematical treatment,
the outputs predicted by the model can be calculated. In
our case, using (18), we will get:
ymk+H
Step 4. Compute the current control action. By establishing
the desired behaviour of our system through the previous
definition of certain reference trajectories, which must
be followed by the outputs of the plant, and applying
principles of PFC, we can determine the expression of
the necessary control action to achieve such behaviour.
Thus, using (36), we will obtain:
uk
Step 5. Apply the calculated control action to WWTP.


### 5.2. Simulation Study

The predictive control experiments were carried out by means of simulation in theMatlab and Simulink environment. Of all the developed software, the most important part was that corresponding to the implementation of the control algorithm, that is, the necessary code for the calculation in each iteration of the control variable, using the analytical and explicit expression obtained ((36) and (37)), including the necessary updating of the state matrices, which depend on time.The study included numerous experiments carried out at different periods of time, some of which are presented in [31, 32]. Numerous tests were carried out with different fuzzy models, different input disturbances, and different output references (mainly between 45 (mg/l) and 65 (mg/l) for the substrate and between 700 (mg/l) and 2000 (mg/l) for the biomass). For each case, different tests with different values of H, the coincidence horizon, were carried out. The simulation time interval chosen was from 0 to 166 hours. For this study, four predictive control experiments were selected, based on fuzzy models that were the result of four identifications considered in Section 3 of this paper: case A, case B, case C, and case D. The dynamic characteristics of these models were summarized in Table 9. In addition to implementing for all cases our FMBPC scheme applied to the WWTP, in the last two cases (C and D), a closed-loop control system for the substrate was also implemented (in parallel), based on a previously tuned PID controller.We will show the results of the selected experiments by means of several graphic representations. For each experiment, the temporal evolutions of the different variables involved are included: the two controlled variables, that is, the substrate concentration in effluent (sk) and the biomass concentration in the reactor (xk), and the control variable, that is, the activated sludge recirculation flow rate (qrk), whose numerical values are precisely the result of the application of the control law obtained in Section 4 of this article. On the other hand, in some of the representations, the simultaneous temporal evolutions of the two perturbations considered have also been included: the input flow rate (qik) and the substrate concentration in the influent (sik). The reason is that it is relevant to show the high level of the disturbances that are being tried to compensate with the control action and, at the same time, to see the evolution of the controlled variables. That is, the objective is to relate the response of the controlled variables to the level and variability of the disturbances. In the case, the graph of the temporal evolution of the substrate is even more relevant, because the control system manages to reduce the value of the substrate concentration to significantly lower levels. The disturbances considered in case A were different from those considered in case B and also different (both) to the perturbations of C and D cases (being the same in case C and in case D).The set points ofs and x, as well as the coincidence horizons considered in each of the four experiments are detailed in Table 10.Table 10
Set points andH of experimental cases.
Case
Set point value
Coincidence horizon
Substrate,ssp (mg/l)
Biomass,xsp (mg/l)
H∈ℤ
A
55
1837
6
B
55
750
3
C
55
1900
6
D
55
1837
180


#### 5.2.1. FMBPC Using Fuzzy Models of A and B Cases

The three graphic representations for the experiment corresponding to case A are shown in Figure13 (sk), Figure 14 (xk), and Figure 15 (qrk).Figure 13
Substrate in the effluent and disturbances (case A).Figure 14
Biomass in the reactor and disturbances (case A).Figure 15
Activated sludge recirculation flow rate: the calculated predictive control variable (case A).And the three graphic representations for the experiment corresponding to case B are shown in Figure16 (sk), Figure 17 (xk), and Figure 18 (qrk).Figure 16
Substrate in the effluent and disturbances (case B).Figure 17
Biomass in the reactor and disturbances (case B).Figure 18
Activated sludge recirculation flow rate: the calculated predictive control variable (case B).The graphic representations corresponding to case A show that the control law obtained is able to maintain the levels of the substrate in effluent around the reference value (with acceptable deviations), despite the strong disturbances of the input flow rate and the pollution of the incoming water (the substrate concentration in the influent). At the same time, the concentration of biomass in the reactor is led to its reference value by following the preestablished reference trajectory. All of this is with a single manipulated variable. On the other hand, from the observation of the evolution of the control actions and the comparison with the disturbances and the outputs, we can say that the variations in the activated sludge recirculation flow rate, determined by the control law, timely counteract the output deviations due to the disturbances and by means of acceptable control efforts.The graphic representations corresponding to case B show, as in case A, that the control law obtained is also able to maintain the levels of the substrate in effluent around the reference value, despite the disturbances of the input flow rate and the pollution of the incoming water (which are different from those of case A). At the same time, the concentration of biomass in the reactor is also led to its reference value (which is also different from that of case A), following the preestablished reference trajectory fairly faithfully (better even than in case A). All of this is also with a single manipulated variable. From the observation of the control actions, we can say that the activated sludge recirculation flow rate, calculated by means of the control law deduced, also compensates in this case, quite well, the output deviations due to the disturbances and also by means of acceptable control efforts. The maximum values of the control variable (as well as its range of variation) are lower than in case A, but also the pollution of the incoming water is much lower in this case than in case A.In none of the cases shown, it was necessary to impose limits on the increase of the control action. However, in some cases, it is necessary to do so to avoid instability. The analysis of this problem requires further study. In relation to the prediction horizon, many tests were carried out with different prediction horizons, but herein, only the results corresponding to anH value, for each of the cases studied, have been shown. A specific study of the influence of the prediction horizon, which is a very important design parameter in predictive control, is needed.


#### 5.2.2. FMBPC Using Fuzzy Models of C and D Cases

The dynamic characteristics of C and D cases are different from those of A and B cases (see Table9). In addition, a previously tuned PID controller has been implemented with the objective of evaluating the FMBPC performance, not only individually but also comparatively. The PID tune parameters were kp=1 and ki=0.1 (with an offset or threshold for the control signal value equal to 775 m3/h).The results of the experiments corresponding to case C can be seen in the graphic representations shown in Figures19–23. And the results of the experiments corresponding to case D can be seen in the graphic representations shown in Figures 24–28. From the observation of such figures, we can summarize the results obtained in the experiments corresponding to C and D cases. In both cases (for the two models used), the substrate is controlled in an acceptable manner considering the strong and permanent disturbances in the input flow rate and in the substrate present in the effluent. In addition, the responses of the controlled substrate with our FMBPC strategy and with the PID strategy are also similar. The same does not occur with the biomass present in the reactor. On the one hand, the answer in case C is a little different from that in case D: tracking the reference trajectory, it is more precise in case C than in case D, probably due to the best VAF indexes in the first case, especially for biomass. However, the case D model is more realistic, due to the validation procedure and, therefore, it would be expected that in other areas of operation, further away from those of the identification, the case D model would respond better than the case C model. In any case, as we said at the beginning, this type of interpretation would require a broader study. The verification of the usefulness of the proposed FMBPC strategy has more weight in this work. Thus, with respect to the comparison of our FMBPC strategy and the PID strategy in the control of biomass, the differences are clear: by our FMBPC strategy, the biomass is acceptably controlled and simultaneously also the substrate, acting with a single manipulated variable, while the PID, acting with the same manipulated variable, manages to control the substrate, but at the cost of losing control of the biomass.Figure 19
Effluent substrate byFMBPC and disturbances (case C).Figure 20
Effluent substrate byFMBPC and PID (case C).Figure 21
Reactor biomass byFMBPC and disturbances (case C).Figure 22
Reactor biomass byFMBPC and PID (case C).Figure 23
Control variables calculated byFMBPC and PID strategies (case C).Figure 24
Effluent substrate byFMBPC and disturbances (case D).Figure 25
Effluent substrate byFMBPC and PID (case D).Figure 26
Reactor biomass byFMBPC and disturbances (case D).Figure 27
Reactor biomass byFMBPC and PID (case D).Figure 28
Control variables calculated byFMBPC and PID strategies (case D).


#### 5.2.3. Performance Evaluation

The goal of the proposed control strategy is to follow the output references as closely as possible. To evaluate the degree of performance in the experiments carried out, theintegral square (or quadratic) error index, known as the ISE index, was taken as a criterion. In the discrete case, the ISE corresponding to each output is the sum of the differences between the reference of the output (set point) and the values of the outputs, squared. The results obtained in the experiments corresponding to C and D cases can be seen in Table 11.Table 11
Performance evaluation by theISE index.
Case
Out 1 (substrate)
Out 2 (biomass)
Both outputs (sum)
FMBPC
PID
FMBPC
PID
FMBPC
PID
C
3.9824e + 05
4.0415e + 05
1.1959e + 05
4.3241e + 06
5.1783e + 05
4.7282e + 06
D
3.9619e + 05
4.0415e + 05
1.2063e + 06
6.7338e + 06
1.6025e + 06
7.1380e + 06Analyzing the numerical results of Table11, we can conclude that, for the two cases analyzed, the performance of the FMBPC strategy is a little better than the PID strategy, for the substrate, and rather better for biomass, something that was also evident after the qualitative analysis of the graphics shown. For biomass, our FMBPC strategy reduces 97.2% of the ISE of the PID in case C and 82% in case D.


## 5.1. FMBPC Scheme

The implemented control strategy is in the field ofnonlinear predictive control (NLMPC) and uses principles and methodology of functional predictive control (PFC) in the deduction of the control law. The basic mechanisms of operation of this strategy have been reflected in Figure 12.Figure 12
FMBPC scheme with an internal fuzzy model.As can be seen, ourFMBPC scheme uses an internal fuzzy model to carry out the predictions and uses PFC principles for the calculation of the control law. The controller needs the information related to the set point of the outputs and the corresponding reference trajectories. In addition, in each iteration (sampling period), state matrices must be updated, because they depend on the premise or antecedent vector (and ultimately also of time). Therefore, the controller must also have the instantaneous information of all the variables that constitute the antecedent vector.The task sequence that must be executed for the implementation of ourFMBPC strategy and that constitute the control algorithm is shown in Algorithm 4.Algorithm 4: Control algorithm for FMBPC implementation.
Repeat at eachkth instant (sampling period, Ts):
Step 1. Update the current premise vector. Measurement of
all the necessary variables for the determination of
the current premise or antecedent vector,xa:
y1k−1,y2k−1,u1k−1,u2k−1,u3k−1,u3k−2
Step 2. Update the state matrices. The TS fuzzy model was
formalized in the form of linear time-varying state space
equations, and the state matrices depend on the current
premise vector. Therefore, at each instant, the following
parameters must be updated, using (14):
A¯m,B¯m,R¯m, and C¯m
Step 3. Calculate the predicted model outputs. The use of a
model to predict future behaviour is essential in MPC.
By means of the appropriate mathematical treatment,
the outputs predicted by the model can be calculated. In
our case, using (18), we will get:
ymk+H
Step 4. Compute the current control action. By establishing
the desired behaviour of our system through the previous
definition of certain reference trajectories, which must
be followed by the outputs of the plant, and applying
principles of PFC, we can determine the expression of
the necessary control action to achieve such behaviour.
Thus, using (36), we will obtain:
uk
Step 5. Apply the calculated control action to WWTP.


## 5.2. Simulation Study

The predictive control experiments were carried out by means of simulation in theMatlab and Simulink environment. Of all the developed software, the most important part was that corresponding to the implementation of the control algorithm, that is, the necessary code for the calculation in each iteration of the control variable, using the analytical and explicit expression obtained ((36) and (37)), including the necessary updating of the state matrices, which depend on time.The study included numerous experiments carried out at different periods of time, some of which are presented in [31, 32]. Numerous tests were carried out with different fuzzy models, different input disturbances, and different output references (mainly between 45 (mg/l) and 65 (mg/l) for the substrate and between 700 (mg/l) and 2000 (mg/l) for the biomass). For each case, different tests with different values of H, the coincidence horizon, were carried out. The simulation time interval chosen was from 0 to 166 hours. For this study, four predictive control experiments were selected, based on fuzzy models that were the result of four identifications considered in Section 3 of this paper: case A, case B, case C, and case D. The dynamic characteristics of these models were summarized in Table 9. In addition to implementing for all cases our FMBPC scheme applied to the WWTP, in the last two cases (C and D), a closed-loop control system for the substrate was also implemented (in parallel), based on a previously tuned PID controller.We will show the results of the selected experiments by means of several graphic representations. For each experiment, the temporal evolutions of the different variables involved are included: the two controlled variables, that is, the substrate concentration in effluent (sk) and the biomass concentration in the reactor (xk), and the control variable, that is, the activated sludge recirculation flow rate (qrk), whose numerical values are precisely the result of the application of the control law obtained in Section 4 of this article. On the other hand, in some of the representations, the simultaneous temporal evolutions of the two perturbations considered have also been included: the input flow rate (qik) and the substrate concentration in the influent (sik). The reason is that it is relevant to show the high level of the disturbances that are being tried to compensate with the control action and, at the same time, to see the evolution of the controlled variables. That is, the objective is to relate the response of the controlled variables to the level and variability of the disturbances. In the case, the graph of the temporal evolution of the substrate is even more relevant, because the control system manages to reduce the value of the substrate concentration to significantly lower levels. The disturbances considered in case A were different from those considered in case B and also different (both) to the perturbations of C and D cases (being the same in case C and in case D).The set points ofs and x, as well as the coincidence horizons considered in each of the four experiments are detailed in Table 10.Table 10
Set points andH of experimental cases.
Case
Set point value
Coincidence horizon
Substrate,ssp (mg/l)
Biomass,xsp (mg/l)
H∈ℤ
A
55
1837
6
B
55
750
3
C
55
1900
6
D
55
1837
180


### 5.2.1. FMBPC Using Fuzzy Models of A and B Cases

The three graphic representations for the experiment corresponding to case A are shown in Figure13 (sk), Figure 14 (xk), and Figure 15 (qrk).Figure 13
Substrate in the effluent and disturbances (case A).Figure 14
Biomass in the reactor and disturbances (case A).Figure 15
Activated sludge recirculation flow rate: the calculated predictive control variable (case A).And the three graphic representations for the experiment corresponding to case B are shown in Figure16 (sk), Figure 17 (xk), and Figure 18 (qrk).Figure 16
Substrate in the effluent and disturbances (case B).Figure 17
Biomass in the reactor and disturbances (case B).Figure 18
Activated sludge recirculation flow rate: the calculated predictive control variable (case B).The graphic representations corresponding to case A show that the control law obtained is able to maintain the levels of the substrate in effluent around the reference value (with acceptable deviations), despite the strong disturbances of the input flow rate and the pollution of the incoming water (the substrate concentration in the influent). At the same time, the concentration of biomass in the reactor is led to its reference value by following the preestablished reference trajectory. All of this is with a single manipulated variable. On the other hand, from the observation of the evolution of the control actions and the comparison with the disturbances and the outputs, we can say that the variations in the activated sludge recirculation flow rate, determined by the control law, timely counteract the output deviations due to the disturbances and by means of acceptable control efforts.The graphic representations corresponding to case B show, as in case A, that the control law obtained is also able to maintain the levels of the substrate in effluent around the reference value, despite the disturbances of the input flow rate and the pollution of the incoming water (which are different from those of case A). At the same time, the concentration of biomass in the reactor is also led to its reference value (which is also different from that of case A), following the preestablished reference trajectory fairly faithfully (better even than in case A). All of this is also with a single manipulated variable. From the observation of the control actions, we can say that the activated sludge recirculation flow rate, calculated by means of the control law deduced, also compensates in this case, quite well, the output deviations due to the disturbances and also by means of acceptable control efforts. The maximum values of the control variable (as well as its range of variation) are lower than in case A, but also the pollution of the incoming water is much lower in this case than in case A.In none of the cases shown, it was necessary to impose limits on the increase of the control action. However, in some cases, it is necessary to do so to avoid instability. The analysis of this problem requires further study. In relation to the prediction horizon, many tests were carried out with different prediction horizons, but herein, only the results corresponding to anH value, for each of the cases studied, have been shown. A specific study of the influence of the prediction horizon, which is a very important design parameter in predictive control, is needed.


### 5.2.2. FMBPC Using Fuzzy Models of C and D Cases

The dynamic characteristics of C and D cases are different from those of A and B cases (see Table9). In addition, a previously tuned PID controller has been implemented with the objective of evaluating the FMBPC performance, not only individually but also comparatively. The PID tune parameters were kp=1 and ki=0.1 (with an offset or threshold for the control signal value equal to 775 m3/h).The results of the experiments corresponding to case C can be seen in the graphic representations shown in Figures19–23. And the results of the experiments corresponding to case D can be seen in the graphic representations shown in Figures 24–28. From the observation of such figures, we can summarize the results obtained in the experiments corresponding to C and D cases. In both cases (for the two models used), the substrate is controlled in an acceptable manner considering the strong and permanent disturbances in the input flow rate and in the substrate present in the effluent. In addition, the responses of the controlled substrate with our FMBPC strategy and with the PID strategy are also similar. The same does not occur with the biomass present in the reactor. On the one hand, the answer in case C is a little different from that in case D: tracking the reference trajectory, it is more precise in case C than in case D, probably due to the best VAF indexes in the first case, especially for biomass. However, the case D model is more realistic, due to the validation procedure and, therefore, it would be expected that in other areas of operation, further away from those of the identification, the case D model would respond better than the case C model. In any case, as we said at the beginning, this type of interpretation would require a broader study. The verification of the usefulness of the proposed FMBPC strategy has more weight in this work. Thus, with respect to the comparison of our FMBPC strategy and the PID strategy in the control of biomass, the differences are clear: by our FMBPC strategy, the biomass is acceptably controlled and simultaneously also the substrate, acting with a single manipulated variable, while the PID, acting with the same manipulated variable, manages to control the substrate, but at the cost of losing control of the biomass.Figure 19
Effluent substrate byFMBPC and disturbances (case C).Figure 20
Effluent substrate byFMBPC and PID (case C).Figure 21
Reactor biomass byFMBPC and disturbances (case C).Figure 22
Reactor biomass byFMBPC and PID (case C).Figure 23
Control variables calculated byFMBPC and PID strategies (case C).Figure 24
Effluent substrate byFMBPC and disturbances (case D).Figure 25
Effluent substrate byFMBPC and PID (case D).Figure 26
Reactor biomass byFMBPC and disturbances (case D).Figure 27
Reactor biomass byFMBPC and PID (case D).Figure 28
Control variables calculated byFMBPC and PID strategies (case D).


### 5.2.3. Performance Evaluation

The goal of the proposed control strategy is to follow the output references as closely as possible. To evaluate the degree of performance in the experiments carried out, theintegral square (or quadratic) error index, known as the ISE index, was taken as a criterion. In the discrete case, the ISE corresponding to each output is the sum of the differences between the reference of the output (set point) and the values of the outputs, squared. The results obtained in the experiments corresponding to C and D cases can be seen in Table 11.Table 11
Performance evaluation by theISE index.
Case
Out 1 (substrate)
Out 2 (biomass)
Both outputs (sum)
FMBPC
PID
FMBPC
PID
FMBPC
PID
C
3.9824e + 05
4.0415e + 05
1.1959e + 05
4.3241e + 06
5.1783e + 05
4.7282e + 06
D
3.9619e + 05
4.0415e + 05
1.2063e + 06
6.7338e + 06
1.6025e + 06
7.1380e + 06Analyzing the numerical results of Table11, we can conclude that, for the two cases analyzed, the performance of the FMBPC strategy is a little better than the PID strategy, for the substrate, and rather better for biomass, something that was also evident after the qualitative analysis of the graphics shown. For biomass, our FMBPC strategy reduces 97.2% of the ISE of the PID in case C and 82% in case D.


## 5.2.1. FMBPC Using Fuzzy Models of A and B Cases

The three graphic representations for the experiment corresponding to case A are shown in Figure13 (sk), Figure 14 (xk), and Figure 15 (qrk).Figure 13
Substrate in the effluent and disturbances (case A).Figure 14
Biomass in the reactor and disturbances (case A).Figure 15
Activated sludge recirculation flow rate: the calculated predictive control variable (case A).And the three graphic representations for the experiment corresponding to case B are shown in Figure16 (sk), Figure 17 (xk), and Figure 18 (qrk).Figure 16
Substrate in the effluent and disturbances (case B).Figure 17
Biomass in the reactor and disturbances (case B).Figure 18
Activated sludge recirculation flow rate: the calculated predictive control variable (case B).The graphic representations corresponding to case A show that the control law obtained is able to maintain the levels of the substrate in effluent around the reference value (with acceptable deviations), despite the strong disturbances of the input flow rate and the pollution of the incoming water (the substrate concentration in the influent). At the same time, the concentration of biomass in the reactor is led to its reference value by following the preestablished reference trajectory. All of this is with a single manipulated variable. On the other hand, from the observation of the evolution of the control actions and the comparison with the disturbances and the outputs, we can say that the variations in the activated sludge recirculation flow rate, determined by the control law, timely counteract the output deviations due to the disturbances and by means of acceptable control efforts.The graphic representations corresponding to case B show, as in case A, that the control law obtained is also able to maintain the levels of the substrate in effluent around the reference value, despite the disturbances of the input flow rate and the pollution of the incoming water (which are different from those of case A). At the same time, the concentration of biomass in the reactor is also led to its reference value (which is also different from that of case A), following the preestablished reference trajectory fairly faithfully (better even than in case A). All of this is also with a single manipulated variable. From the observation of the control actions, we can say that the activated sludge recirculation flow rate, calculated by means of the control law deduced, also compensates in this case, quite well, the output deviations due to the disturbances and also by means of acceptable control efforts. The maximum values of the control variable (as well as its range of variation) are lower than in case A, but also the pollution of the incoming water is much lower in this case than in case A.In none of the cases shown, it was necessary to impose limits on the increase of the control action. However, in some cases, it is necessary to do so to avoid instability. The analysis of this problem requires further study. In relation to the prediction horizon, many tests were carried out with different prediction horizons, but herein, only the results corresponding to anH value, for each of the cases studied, have been shown. A specific study of the influence of the prediction horizon, which is a very important design parameter in predictive control, is needed.


## 5.2.2. FMBPC Using Fuzzy Models of C and D Cases

The dynamic characteristics of C and D cases are different from those of A and B cases (see Table9). In addition, a previously tuned PID controller has been implemented with the objective of evaluating the FMBPC performance, not only individually but also comparatively. The PID tune parameters were kp=1 and ki=0.1 (with an offset or threshold for the control signal value equal to 775 m3/h).The results of the experiments corresponding to case C can be seen in the graphic representations shown in Figures19–23. And the results of the experiments corresponding to case D can be seen in the graphic representations shown in Figures 24–28. From the observation of such figures, we can summarize the results obtained in the experiments corresponding to C and D cases. In both cases (for the two models used), the substrate is controlled in an acceptable manner considering the strong and permanent disturbances in the input flow rate and in the substrate present in the effluent. In addition, the responses of the controlled substrate with our FMBPC strategy and with the PID strategy are also similar. The same does not occur with the biomass present in the reactor. On the one hand, the answer in case C is a little different from that in case D: tracking the reference trajectory, it is more precise in case C than in case D, probably due to the best VAF indexes in the first case, especially for biomass. However, the case D model is more realistic, due to the validation procedure and, therefore, it would be expected that in other areas of operation, further away from those of the identification, the case D model would respond better than the case C model. In any case, as we said at the beginning, this type of interpretation would require a broader study. The verification of the usefulness of the proposed FMBPC strategy has more weight in this work. Thus, with respect to the comparison of our FMBPC strategy and the PID strategy in the control of biomass, the differences are clear: by our FMBPC strategy, the biomass is acceptably controlled and simultaneously also the substrate, acting with a single manipulated variable, while the PID, acting with the same manipulated variable, manages to control the substrate, but at the cost of losing control of the biomass.Figure 19
Effluent substrate byFMBPC and disturbances (case C).Figure 20
Effluent substrate byFMBPC and PID (case C).Figure 21
Reactor biomass byFMBPC and disturbances (case C).Figure 22
Reactor biomass byFMBPC and PID (case C).Figure 23
Control variables calculated byFMBPC and PID strategies (case C).Figure 24
Effluent substrate byFMBPC and disturbances (case D).Figure 25
Effluent substrate byFMBPC and PID (case D).Figure 26
Reactor biomass byFMBPC and disturbances (case D).Figure 27
Reactor biomass byFMBPC and PID (case D).Figure 28
Control variables calculated byFMBPC and PID strategies (case D).


## 5.2.3. Performance Evaluation

The goal of the proposed control strategy is to follow the output references as closely as possible. To evaluate the degree of performance in the experiments carried out, theintegral square (or quadratic) error index, known as the ISE index, was taken as a criterion. In the discrete case, the ISE corresponding to each output is the sum of the differences between the reference of the output (set point) and the values of the outputs, squared. The results obtained in the experiments corresponding to C and D cases can be seen in Table 11.Table 11
Performance evaluation by theISE index.
Case
Out 1 (substrate)
Out 2 (biomass)
Both outputs (sum)
FMBPC
PID
FMBPC
PID
FMBPC
PID
C
3.9824e + 05
4.0415e + 05
1.1959e + 05
4.3241e + 06
5.1783e + 05
4.7282e + 06
D
3.9619e + 05
4.0415e + 05
1.2063e + 06
6.7338e + 06
1.6025e + 06
7.1380e + 06Analyzing the numerical results of Table11, we can conclude that, for the two cases analyzed, the performance of the FMBPC strategy is a little better than the PID strategy, for the substrate, and rather better for biomass, something that was also evident after the qualitative analysis of the graphics shown. For biomass, our FMBPC strategy reduces 97.2% of the ISE of the PID in case C and 82% in case D.


## 6. Conclusions

In this article, a fuzzy predictive control law in an analytical and explicit way has been developed and has been applied to a multivariable, with disturbances, strongly nonlinear, with a complex dynamic and of a biological nature process. The process was an urban wastewater treatment plant (WWTP) with purification by activated sludge. The case study was developed in a simulation environment. The main conclusion of this study is the capacity of the proposed nonlinear MPC strategy to control a strongly nonlinear and multivariable system, in the presence of strong disturbances, even with adaptation to changes in the operating point with time. OurFMBPC strategy has been able to control two variables of a WWTP (substrate and biomass) simultaneously, making use of a single manipulated variable (sludge recirculation). The second important conclusion is that the performance of the proposed strategy improves that of a PID controller, in a very appreciable way for the case of biomass and with a similar performance for the substrate. The result of the evaluation made on the performance of each strategy (using the ISE index) is, in summary, that the proposed FMBPC approach reduces between 82% and 97.2% the ISE of the PID for the biomass variable (82% for one of the studied cases and 97.2% for the other).Another line of future work could be the search for a variant of the proposed control algorithm that incorporates parameters (degrees of freedom), to be determined by means of optimization, with the aim of avoiding instabilities and improving the operation of the controlled plant. Likewise, the possibility, already mentioned in the introductory section of this article, that the proposed fuzzy control law can be straightforwardly used within adual-mode MPC scheme to handle constraints if needed opens another interesting line of future work both in the field of nonlinear predictive control as well as in the field of intelligent control.



---

*Source: 5720185-2019-01-03.xml*