Analytical Fuzzy Predictive Control Applied to Wastewater Treatment Biological Processes

A novel control fuzzy 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.


Introduction
The traditional strategy of model-based predictive control (MBPC or MPC) [1][2][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][9][10][11][12][13][14][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-called predictive functional control (PFC) [2,[16][17][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 the FMBPC 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 a MBPC 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 the model-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 a wastewater 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 physicalchemical 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 to WWTPs) have been numerous and different, highlighting in our case those that include some variant of nonlinear 2 Complexity model-based predictive control (NLMPC) and, especially, the proposals framed in the field of the fuzzy model-based predictive control (FMBPC).In [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][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.Section 2 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 timedependent 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.

Case Study: Wastewater Treatment Processes
Our case study consists of applying the previously introduced the FMBPC 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.

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-called activated 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.
As already indicated in the introductory chapter, the purification technique used in this plant is consisting of the elimination of organic pollutants by means of activated 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 (q r ) 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 (q p ) is purged from the bottom of the settler, and the other one (q r ) 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 (q i ), the input substrate concentration (s i ), and the sludge recirculation flow rate sent to the reactor (q r ).The first two q i s i , x i s ir , x ir f k q, s, x q out , s x d q 2 q r q p q 2 x b x r s, x, c q Figure 1: Wastewater treatment plant with activated sludge.
3 Complexity 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).
The information relating to the considered inputs and outputs of the plant is shown more detailed in Table 2 including for each of them and also their biological or physicalchemical significance.
The 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.

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 classical Monod 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 Manresa WWTP [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.

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 a fuzzy model-based predictive control (FMBPC) system [14].It is for this reason that in this article, we will concentrate more on the analytical u 3 = q r (u)  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 as FMID (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.Figure 3 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).
The inputs of the considered plant were the following: the incoming water flow rate (influent), q i k ; the substrate concentration in the influent (organic pollution of the incoming water), s i k ; and the activated sludge recirculation flow rate, q r k .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), s k , and the biomass concentration (microorganisms that feed on the substrate and therefore purify the water) in the reactor, x k .
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-Appendix D1).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 ( q i k , s i k , q r k ), we determined by simulation the corresponding output values ( s k , x k ), 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).
Figures 4 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.
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 a Takagi-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 Wastewater treatment plant (activated sludge)

Multiple input
Multiple output The MIMO system to identify (3 inputs and 2 outputs).Biomass concentration in the reactor x 5 Complexity 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 ϕ j x .
The 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 the product 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 Identification data (3 inputs: q i , s i , q r ) q i (input flow rate (influent)) s i (concentration of substrate in influent) q r (flow rate of sludge recirculation) R j : if (x a1 is A j1 and x a2 is A j2 and, … , and x ap is A jp ) then: , mr ; mr: number of rules x a = x a1 , x a2 , … , x ap : antecedent vector x = x 1 , x 2 , ⋯, x q : consequent vector A j1 , A j2 , … , A jp : fuzzy sets related to the components of the antecedent vector rule R j Membership function associated to the fuzzy set A jh : θ A jh ℝ⟶ 0, 1 (smooth function) x ah ⟶ μ A jh x ah (membership grade of x ah with respect to the fuzzy set A jh ) h = 1, … , p Algorithm 1: General form of the Takagi-Sugeno-type fuzzy models.6 Complexity 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.

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: (ii) j = 1, 2, … , mr i .
(iii) mr i is the number of rules of the output y i .
(iv) p i is the number of components of the antecedent vector corresponding to y i .
(v) μ A jp i x ap i is the membership grade of x ap i with respect to the fuzzy set A jp i .
And defining for each output and rule the following membership functions (normalized) for all the antecedent vector the numerical expression of the output will remain as where mr i is the number of rules of the output y i .
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 the FMID 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 the FMID software tool can be summarized in Table 3.
Next, in Table 4, 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.
The result of each identification test carried out consists of a fuzzy model of the Takagi-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 T s Sample time Real number 7 Complexity 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 T s ) and the structure of the rules of the fuzzy model of the two outputs is quite clear.Thus, Algorithm 2: Takagi-Sugeno fuzzy model for substrate Algorithm 3: Takagi-Sugeno fuzzy model for biomass y 2 k = x k .
8 Complexity the recursive dynamic models of y 1 k and y 2 k have the following structure: or alternatively In Tables 5 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.
The identification tool also provides the necessary numerical information relating to the fuzzy sets or values A ij , 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 u 3 k − 1 = q r k − 1 variable, which is one of the components of the antecedent vector of the rules of the output y 1 k , for the case A.
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, x a|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, x a|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.
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: where the involved variables are the same as those detailed in (8).

Validation of the Identified Fuzzy
Models.An essential part of any identification process is the validation of the identified models.The FMID 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 Section 3.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.
For case A and case B, dependence on y 1 k and y 2 k outputs with respect to  including in addition the VAF index, which measures the goodness of the identification for each of the two partial identifications.
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 Figure 8.The same Different (partition) A

Complexity
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 (Figure 7), 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.

Results of the Validation Process for Case C and Case D.
In case C and case D, dependence on y 1 k and y 2 k outputs with respect to u 1 k − 1 = q i k − 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 inputoutput 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.
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 (Figure 9), 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 u 1 k − 1 = q i k − 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    12 Complexity 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 the Takagi-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: and the state matrices are where 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 x a (through β j x a ) and therefore also depend on the kth instant, because x a 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 x a and subsequently β j x a .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).

The Control Law in Analytical Form
The next step of our fuzzy 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][17][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-called equivalence 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: 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).

Predictions
Based on the Model.The control actions in the various existing predictive 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 y m k + 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.14 Complexity 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 Appendix A with the objective of facilitating the followup of the article.The calculation of Hstep ahead prediction is done under the assumption of constant future manipulated variables, that is, The final expression deduced is the following (see (A.26), at the end of Appendix A): where , the order 2 identity matrix, and P 1010 = 1 0 1 0 20 4.2.Deduction of the Control Law in an Analytical and Explicit Way.As from the expression of y m k + 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 u a k in (19), obtaining an expression that will be a function of y m k + H and we will subsequently extract, of the vector variable u a k obtained, the scalar variable u k , which will also be a function of y m k + 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 y m k + H .And as from the relations obtained, we will deduct the expression of the control action u k 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.

Manipulated Variable u k as a Function of y m k + H .
Leaving at (19) to one side of equality, the term in which u a k is multiplying and grouping together terms for simplification, defining matrix M a as follows: we will have ( 21) simplified, obtaining and multiplying on the left, on both sides of equality, by the matrix M a −1 , we will have the following matrix expression: and replacing u a k by its generic expression specified in ( 12), we will have where u k − 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, u k − 1 will be a data, a specific numerical value saved, while u k will be the unknown quantity to find or determine.
If matrix (25), with u k − 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 u k by means of the following expression: being and replacing u a k by its expression given by (24), we would finally have 15 Complexity which is the expression that will allow the determining of the variable u k which will be necessary at the kth sampling instant so that the model can attain, H periods later (i.e., at the k + H th instant), a certain vector valuey m k + H .

Control Action u k
Necessary for the Follow-Up of the Output Reference Trajectories.The control action u k 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: where y set point i k is the reference or desired value of the ith output y i k at the kth instant, y r i k and y r i k + 1 , the values in k and in k + 1, respectively, of the reference trajectory, and a r i and b r i 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 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 single coincidence point.Such a point will correspond to the instant k + H, where H is called the coincidence horizon.
The control action at each kth instant, u k , should therefore be such that the plant outputs match the corresponding reference trajectories, H sampling periods later, such that y i k + H equals y r i k + H . Therefore, the desired plant output increment should be equal to y r i k + H − y i k , for each ith output (i = 1, 2), that is,

number of outputs 31
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-called equivalence 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: and matching ( 31) and (32) (i.e., matching the desired plant output increments with the predicted model output increments: ), we will have

number of outputs 33
The control action at each kth instant, u k , 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: Finding now y m k + H in (34), we will have that the predicted model output in the k + 1 th time instant should satisfy the following expression: Replacing (28) y m k + H by the second member of equality (35), which we have just obtained, we will finally have the expression that we were looking for u k , 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 u k is the following, where M a is given by ( 22) and P 10 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
where the term y r k + 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 16 Complexity reference on the prediction horizon in Appendix C. The result is as follows (see (C.10) at the end of said appendix): 36), complemented with (37), constitutes the analytical and explicit expression of the control law u k 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.

FMBPC Applied to WWTP
The main objective of this work is the deduction of a fuzzy 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.

FMBPC Scheme.
The implemented control strategy is in the field of nonlinear 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.
As can be seen, our FMBPC 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 our FMBPC strategy and that constitute the control algorithm is shown in Algorithm 4.

Simulation Study. The predictive control experiments
were carried out by means of simulation in the Matlab 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 (s k ) and the biomass concentration in the reactor (x k ), and the control variable, that is, the activated sludge recirculation flow rate (q r k ), 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 (q i k ) and the substrate concentration in the influent (s i k ).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 17 Complexity 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 of s and x, as well as the coincidence horizons considered in each of the four experiments are detailed in Table 10.

FMBPC Using Fuzzy
Models of A and B Cases.The three graphic representations for the experiment corresponding to case A are shown in Figure 13 (s k ), Figure 14 (x k ), and Figure 15 (q r k ).
And the three graphic representations for the experiment corresponding to case B are shown in Figure 16 (s k ), Figure 17 (x k ), and Figure 18 (q r k ).
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 Repeat at each kth instant (sampling period, T s ): Step 1. Update the current premise vector.Measurement of all the necessary variables for the determination of the current premise or antecedent vector, x a : 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): By means of the appropriate mathematical treatment, the outputs predicted by the model can be calculated.In our case, using ( 18), we will get: 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: u k Step 5. Apply the calculated control action to WWTP.18 Complexity 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  19 Complexity 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 an H 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.

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 Table 9).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 k p = 1 and k i = 0 1 (with an offset or threshold for the control signal value equal to 775 m 3 /h).
The results of the experiments corresponding to case C can be seen in the graphic representations shown in Figures 19-23.And the results of the experiments  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 Predictive control variable FMBPC: q r (m 3 /h) Control variable PID: q r (m 3 /h)  21 Complexity 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.

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, the integral 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.
Analyzing the numerical results of Table 11, 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.

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.Our FMBPC 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 a dualmode 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.Predictive control variable FMBPC: q r (m 3 /h) Control variable PID: q r (m 3 /h)  10) and ( 11), we will develop the expression for a finite number of cases, giving to H consecutive numerical values, integers (beginning with H = 1).Then we will reason by means of induction.We detail the deductive-inductive process below, beginning with the case H = 1 .
If we replace k with k + 1 in state equation ( 11), we will have and then replacing z m k + 1 with the expression specified in state equation (10), this will give If we replace k with k + 2 in state equation ( 10) and then we perform the change of variable (intermediate) k ′ = k + 1 and then use twice consecutively state equation ( 10), we will have the following, under the hypothesis of the invariability of the coefficients of state equations for the predictions (from the kth instant onwards): Imposing now as a hypothesis that u k is maintained constant during the prediction horizon, that is, it can easily be proved that where and now, making use of equality (A.4), the corresponding development to y m k + 2 initiated in (A.3) will finally be By replacing k with k + 3 in state equation ( 11) and using successively state equation (10) and performing the necessary changes of the variable, we will have the following, also under the hypothesis of the invariability of the coefficients of state equations for the predictions (from the kth instant onwards): Considering now again the hypothesis that u k is maintained constant during the prediction horizon, expressed in abbreviated form as (A.4) and using P 1010 (the matrix specified in (A.6)), it can be proved that Taking into account once again the hypothesis of the constant maintenance of u k in the prediction horizon, (A.4) and using again the P 1010 matrix, (A.6), it can be proved that Carrying out the necessary development, totally analogous to that of the previous cases and which we will omit in the interests of simplification, and using equality (A.11), the following final expression would be reached for the output predicted by the model for the instant k + 4 : Induction Hypothesis.Once the development corresponding to y m k + H has been made for the first four values of H, we will assume that we can generalize the expressions obtained and apply them to the pth case.For the case H = p, therefore, the final expression would be where P 1010 is the matrix specified in (A.6) and I is the order 2 identity matrix, that is, Next, we will reason by means of induction.Assuming that the expression corresponding to case H = p is correct, it will need to be shown that it is also correct for H = p + 1, that is, that it would also be complied with in it if we replace p with p + 1.In order to do so, we will use the state equations (in a similar way to what was done in the previous cases and with the same hypotheses) and (A.13), which is corresponding to case H = p.
As from state equation (11), making a trivial change of variable, we will have and comparing (A.13) and (A.15) and equalling the righthand side of both equalities, we will have the following: Taking into account again the state equations and also the previous result, we can now approach the development corresponding to case H = p + 1.First, we will use state equation (11) and then we will make the intermediate change of variable k p = k + p and then we will use state equation (10) to develop z m k p + 1 On the other hand, considering again the hypothesis that u k is maintained constant during the prediction horizon, expressed in (A.4), and with P 1010 being the matrix specified in (A.6), the compliance with the following equality sequence can easily be proved: and taking into account (A.16) and (A.18), replacing the terms which correspond with (A.17), and operating in 24 Complexity an appropriate manner, we will obtain the development corresponding to y m k + p + 1 , initiated in (A.17), as follows: This last result compared with (A.13) shows that the predictions in case H = p + 1 satisfy the same formula as those of case H = p.We can therefore conclude that the validity of (A.13) can be extended to p ∈ ℤ + , p ≥ 1, and understanding that the factors A m p−n which will appear (with n ∈ ℤ + ), one must comply with p − n ≥ 1.This expression may therefore be considered to be the general formula of the calculation of predictions at the kth instant, based on our statespace fuzzy model.
And now, once the demonstration has been concluded, we will formalize the mathematical expression of the calculation of predictions in a more general manner, representing the prediction horizon with H.In short, the general expression of y m k + H obtained will be as follows (A.20): Expression (A.20) may be presented in another way with the objective being the appearance as main factors of z m k , u a k , and R m .Regrouping terms, therefore, we will have and making use of the following development of matrix algebra (for square matrixes) A m 0 = I with I being the order 2 identity matrix ,

A 22 or what is the same
then we will have (considering that p = H − 1) and replacing (A.24) in (A.21) in the two terms where it appears to leave finally being that , the order 2 identity matrix and P 1010 = 1 0 1 0

A 25
The final expression obtained (A.26) is the general analytical expression that must be used for calculating the predictions, that is, for the calculation at instant k of the outputs predicted by the fuzzy model for the instant k + H , with H being the prediction horizon.It can easily be proved that the final formula obtained for y m k + H is also satisfied for each of the four particular cases developed in the induction process (H = 1, 2, 3, 4).

B. Relationship between the Parameters of the Reference Model
The reference model of each of the two outputs of our system is given by reference trajectories that should gradually approach the corresponding references (set points).The 25 Complexity discrete equations of these trajectories, already introduced in (29), are as follows: y r i k + 1 = a r i y r i k + b r i y set point i k , i = 1, 2 number of outputs B 1 As is logical, for the reference trajectories to follow their corresponding references, the reference model gain must be the unit.The imposition of this condition will cause a certain relationship between the parameters a r i and b r i , which is precisely what we wish to determine.In order to do so, firstly, we will calculate the discrete transfer function inz of the reference models and then we will impose the unit gain condition.
By applying the ztransform to equality (B.1), we will have and grouping terms together and operating properly And now, making use of the expression of the discrete gain that derives from the application of the final value theorem and imposing the unit gain condition, we will have and in this way, we will finally obtain the relationship between the parameters of the reference model that we were searching (for each of the two outputs)

C. Reference Trajectory on the Prediction Horizon
The mathematical model of the reference trajectories, introduced in (B.1) and in (29), is a discrete time recursive model that allows the calculation of the value of the output variable (of that model) at the instant k + 1 .The objective of this appendix is to show in an abbreviated manner the procedure for obtaining the expression corresponding to the output of the reference model at the instant k + H .In order to determine the expression corresponding to y r i k + H , for i = 1, 2, we will reason by means of induction (omitting the indication of the values of i i = 1, 2 , except at the beginning and at the end so as not to overload the development).
For this first case, we will use the expression defining the reference trajectories, (B.1) and also (B.5), that is, the relationship between the parameters of the trajectories y r i k + 2 = ⋯ = a 2 r i y r i k + 1 + a r i 1 − a r i y set point i k = a 2 r i y r i k + 1 − a 2 r i y set point i k C 6 Induction Hypothesis.The expression corresponding to the following cases would be obtained in a similar way, deducing (with the consideration of the maintaining of the constant reference with constant values on the prediction horizon) the following generic expression, which will be taken as the induction hypothesis: H = p ; p ∈ ℤ + , p ≥ 1 y r i k + p = a p r i y r i k + 1 − a p r i y set point i k C 7 If we suppose the above formula to be correct, it can be proved easily that it is also complied with for p + 1 and reasoning by induction, the general expression will be the following (considering the reference on the prediction horizon to be constant), for each of the two outputs:

D. Mathematical Nonlinear Model of the Wastewater Treatment Process with Activated Sludge
The mathematical model of the wastewater treatment process taken as reference in this article is founded on the classical Monod and Maynard-Smith model (with the assumption of a perfectly mixed tank reactor) and it is obtained taking into account the corresponding mass balances of the substrate, biomass, and oxygen.The equations of the model, as well as its variables and parameters, are the following: (i) Aerated biological reactor x ir = x i q i + x r q r q s ir = s i q i + sq r q D 1 (ii) Secondary settler (three layers with increasing concentrations of biomass) (iii) Equations of equilibrium of the flow rates q = q i + q r q out = q i − q p q 2 = q r + q p D 3 (iv) Variables and parameters (with the areas expressed in m 2 , the volumes in m 3 , the concentrations in mg/l, and the flow rates in m 3 /h): A: Settler area V: Reactor volume 27 Complexity D.1.Parameters of the Manresa WWTP.The numerical values of the Manresa wastewater treatment plant (taken as reference) [35] are specified in Table 12.Height of the lower layer of the settler q i : Input flow rate (contaminated water/ influent) q: Reactor input-output flow rate q 2 : Total sludge recirculation flow rate q p : Purge flow rate (excess sludge) q r : Sludge recirculation flow rate (to the reactor) q out : Output flow rate (purified water/ effluent) Equivalence coefficient between cell growth and oxygen consumption rate.

[
the parameter value of the first row and the third column, equal to 2, implies two consecutive terms of u 3 (in discrete time)) [Row 2]: y 2 k depends on u 1 k − 1 and u 3 k − 1 Row 1]: they are 1 transport delay for the three inputs (u 1 , u 2 , and u 3 ) with respect to y 1 k [Row 2]: they are 1 transport delay for the u 1 and u 3 inputs with respect to y 2 9.35•10 −1 −5.22•10 −2 −3.96•10 −2 2.62•10 2

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 Figure 7.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), Validation of identified fuzzy model (wastewater treatment plant) wastewater treatment plant (classic model equations) Output 1 of the fuzzy model identified Output 2 of the fuzzy model identified

Figure 7 :
Figure 7: Validation of the identified fuzzy model corresponding to case A for output 1, s k , and output 2, x k .

Figure 8 :
Figure 8: Validation of the identified fuzzy model corresponding to case B for output 1, s k , and output 2, x k .
the fuzzy model identified Output 2 of the wastewater treatment plant (classic model equations) Output 2 of the fuzzy model identified

Figure 9 :
Figure 9: Validation of the identified fuzzy model corresponding to case C for output 1, s k , and output 2, x k .
the wastewater treatment plant (classic model equations) Output 1 of the fuzzy model identified Output 2 of the wastewater treatment plant (classic model equations) Output 2 of the fuzzy model identified

Figure 10 :
Figure 10: Validation of the identified fuzzy model corresponding to case D for output 1, s k , and output 2, x k .
) z m k is the extended-state vector, integrated by the outputs and disturbances of the WWTP, at the kth instant (ii) z m k + 1 is the extended-state vector, at the k + 1 th instant (iii) y m k is the output of the model in the state space, integrated by the outputs of the WWTP, at the kth instant (iv) u a k input 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 matrices and the corresponding expressions are

β 26 x
a = 0, mr = max ⋅ mr 1 , mr 2 a j1 , a j2 , b j1 , b j2 , b j3 , b j4 , r j are the coefficients of the antecedent vector and independent term in the jth rule of the output 1 fuzzy model of the WWTP.(ii) a * j2 , b * j1 , b * j3 , r * j are the coefficients of the antecedent vector and independent term in the jth rule of the output 2 fuzzy model of the WWTP (with a * 62 = 0, b * 61 = 0, b * 63 = 0, and r * 6 = 0).

Figure 12 :
Figure 12: FMBPC scheme with an internal fuzzy model.

and C m Step 3 .
Calculate the predicted model outputs.The use of a model to predict future behaviour is essential in MPC.

Figure 13 :
Figure 13: Substrate in the effluent and disturbances (case A).

2 H+ 1 − a r i y set point i k + 1 = a 2 + 1 − a r i y set point i k + 1 C 4
y r i k + 1 = a r i y r i k + b r i y set point i k , b r i = 1 − a r i , i = 1, 2 number of outputs , C1and replacing b r i in the expression of the trajectoryy r i k + 1 = a r i y r i k + 1 − a r i y set point i k C = 2 k + H = k + 2 We develop k + H for H = 2, considering k + 2 = k + 1 + 1 andconsidering the trivial change of variable k′ = k + 1, which would leave k + 2 = k ′ + 1.Using again (B.1) that define the reference trajectories, but with k ′ instead of k, and making b r i = 1 − a r i , we will havey r i k + 2 = y r i k′ + 1 = a r i y r i k′ + 1 − a r i y set point i k′ , C 3and reversing the change of variable (k′ = k + 1) and again using (B.1) to develop y r i k + 1 , we will havey r i k + 2 = a r i y r i k + 1 + 1 − a r i y set point i k + 1 = a r i a r i y r i k + 1 − a r i y set point i k r i y r i k + a r i 1 − a r i y set point i kNow considering the reference on the prediction horizon to be constanty set point i k = y set point i k + 1 = ⋯ = y set point i k + H − 1 , C 526Complexity the development of y r i k + 2 , initiated in (C.4), will continue as follows:

H
2 number of outputs , C 9 expressions that can be formalized jointly by means of the following (single) matrix expression:y r k + H = A rH y r k + I − A rH y set point k ∈ ℤ + , H ≥ 1, C 10with y r k + H , y r k e y set point k being size-two columnvectors that group together the components corresponding to both outputs.
= q out x b − q out x d − Aν s x d Al b dx b dt = qx − q out x b − q 2 x b + Aν s x d − Aν s x b Al r dx r dt = q 2 x b − q 2 x r + Aν s x b ν s x d = nnr x d exp aar x d ν s x b = nnr x b exp aar x b D 2 Biomass concentration in the reactor s:Substrate concentration in the reactor c:Oxygen concentration in the reactor x d :Biomass concentration in the top layer of the settler x b :Biomass concentration in the middle layer of the settler x r :Biomass concentration in the lower layer of the settler l d :Height of the top layer of the settler l b :Height of the middle layer of the settler l r :

Table 1 :
Main variables of the WWTP with activated sludge.

Table 2 :
Inputs and outputs of the wastewater treatment plant.

Table 3 :
Structural parameters for fuzzy identification.

Table 4 :
Structural parameters in case D.

Table 5 :
Consequent parameters for substrate y 1

Table 7 :
Consequent parameters for biomass y 2

Table 8 :
Cluster centers for biomass y 2

Table 9 :
Parameters and use of input-output data in 4 cases.

Table 10 :
Set points and H of experimental cases.
Figure 28: Control variables calculated by FMBPC and PID strategies (case D).The objective of this appendix is the obtaining of the general expression of the predictions of the wastewater treatment plant outputs, as from the fuzzy model identified and subsequently formalized in the state space.It is a case of deducing the general mathematical expression corresponding to y m k + 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.Making use of the model obtained, that is, of the state equations (

Table 11 :
Performance evaluation by the ISE index.

Table 12 :
Parameters of the Manresa WWTP.