Practical Computational Approach for the Stability Analysis of Fuzzy Model-Based Predictive Control of Substrate and Biomass in Activated Sludge Processes

: This paper presents a procedure for the closed-loop stability analysis of a certain variant of the strategy called Fuzzy Model-Based Predictive Control (FMBPC), with a model of the Takagi-Sugeno type, applied to the wastewater treatment process known as the Activated Sludge Process (ASP), with the aim of simultaneously controlling the substrate concentration in the efﬂuent (one of the main variables that should be limited according to environmental legislations) and the biomass concentration in the reactor. This case study was chosen both for its environmental relevance and for special process characteristics that are of great interest in the ﬁeld of nonlinear control, such as strong nonlinearity, multivariable nature, and its complex dynamics, a consequence of its biological nature. The stability analysis, both of fuzzy systems (FS) and the very diverse existing strategies of nonlinear predictive control (NLMPC), is in general a mathematically laborious task and difﬁcult to generalize, especially for processes with complex dynamics. To try to minimize these difﬁculties, in this article, the focus was placed on the mathematical simpliﬁcation of the problem, both with regard to the mathematical model of the process and the stability analysis procedures. Regarding the mathematical model, a state-space model of discrete linear time-varying (DLTV), equivalent to the starting fuzzy model (previously identiﬁed), was chosen as the base model. Furthermore, in a later step, the DLTV model was approximated to a local model of type discrete linear time-invariant (DLTI). As regards the stability analysis itself, a computational method was developed that greatly simpliﬁed this difﬁcult task (in a local environment of an operating point), compared to other existing methods in the literature. The use of the proposed method provides useful conclusions for the closed-loop stability analysis of the considered FMBPC strategy, applied to an ASP process; at the same time, the possibility that the method may be useful in a more general way, for similar fuzzy and predictive strategies, and for other complex processes, was observed.


Introduction
In general, the control of non-linear industrial processes is not an easy task and it is further complicated in the case of processes of a biological nature, especially if the process is affected by major disturbances.One of the possible alternatives to approach the control of this type of systems could be the well-known and effective advanced control strategy called Model-Based Predictive Control (MBPC or MPC) [1][2][3][4].However, this strategy is not associated with a single control algorithm, but rather it is a methodology or procedure to search for optimal control actions based on the predictions of the behavior of the process to be controlled, provided by some mathematical model thereof (in classic MPC, the search for the optimal control actions is carried out imposing the minimization of a certain cost function).Consequently, there are many types of predictive control, depending on the type of model chosen, the method of identifying the model, the procedure or algorithm to determine the control actions, and the control law itself that is finally adopted, among other aspects.Regarding the nonlinear nature of the processes to be controlled, for systems whose nonlinearity is not very high, if there is or is possible to obtain a model or several linearized sub models of the process that are sufficiently valid under certain hypotheses and restrictions, then the well-known and proven techniques of the so-called linear MPC (highly consolidated) could be enough.However, if the system is highly nonlinear, it is necessary to consider a predictive control strategy based on nonlinear models that more accurately represent the behavior of the process.This approach is framed (along with other variants) within the so-called nonlinear model predictive control (nonlinear MPC, or NLMPC) [5][6][7][8][9].For highly nonlinear, complex, or unknown systems, one class of models that seems quite adequate for predictive control purposes are fuzzy models (FM), especially the Takagi-Sugeno (TS) type fuzzy models [10], which are very useful because the outputs of the sub models are given in the form of affine linear combinations.In addition, there are consolidated procedures for obtaining TS fuzzy models from series of numerical inputoutput data (collected by simulation or from real experiments).The fuzzy models obtained using data-based identification techniques have the enormous advantage of being able to directly capture the dynamics of such systems.In addition, this type of modeling is also able to deal with multivariable systems in a straightforward manner.In other words, within the framework of predictive control, the strategy known as Fuzzy Model-Based Predictive Control (FMBPC) and especially the variant that uses TS fuzzy models [11][12][13][14][15][16][17] can be a good alternative, in general, to approach the control of highly nonlinear and complex systems, among them, the systems of a biological nature; in particular, the socalled activated sludge process (ASP), which is a very common purification procedure in Wastewater Treatment Plants (WWTP) [18][19][20].
In this article, a certain FMBPC multivariable control strategy is analyzed [21][22][23], which is framed at the same time in the so-called Predictive Functional Control (PFC) [1,24,25].In such a strategy, the considered prediction original model is a TS fuzzy model, obtained by identification from input-output data, later represented in the state-space by means of an equivalent model of linear time-varying (LTV) type.As for the control action, it is calculated (at each sampling instant) by means of an analytical and explicit expression (with time-dependent coefficients), deduced using the state-space LTV model and applying the so-called equivalence principle (used in PFC and related to the concepts of reference trajectory and coincidence point), rather than by minimizing a cost function.This control algorithm is applied to an ASP biological process, taken as a case study, with the aim of simultaneously controlling the substrate concentration in the effluent and the biomass concentration in the reactor, with the activated sludge recirculation flow-rate being the only control variable used.The main objective of this paper is to study some issues of the closed-loop stability analysis of the specific FMBPC strategy that has just been described, by a practical computational approach.
The stability analysis of nonlinear control systems is much more complex, in general, than the linear systems stability analysis.Moreover, there are many different approaches, depending on both the type of control law and the type of system to be controlled.The objective of this work is not to analyze all these possibilities, but to focus on the field of fuzzy predictive control, especially in control strategy and in the case study presented in [23].However, it is important to mention some of the works that have addressed the stability of nonlinear systems-it was highlighted in works [26][27][28][29][30][31][32], among many others.In the case of the FMBPC strategy, the stability analysis may be even more complicated than other types of nonlinear control strategies, due, on the one hand to the particular characteristics of the nonlinear predictive control and on the other hand, to the structure of the Fuzzy Systems (FS) [11,33], initially in the form of rules, which are intrinsically very nonlinear, in general, and that logically does not facilitate the use of already consolidated classical procedures.The stability analysis of the FMBPC strategy is, moreover, more uncertain and diverse than other strategies, as a consequence of the great variety of types of fuzzy models and possible parameterizations or choices (antecedent vector and consequent vector, set of fuzzy values for each variable, form and parameters of the membership functions, number of rules, and others).Potentially, many alternatives could be explored to analyze the stability of this type of strategy, but there are two fields in which it seems more logical to focus on initially, and explore the lines of work on this topic: the field of NLMPC control systems and the field of FS systems.Each of the two fields has its particularities, but a large part of the stability analysis approaches developed in them share a theoretical framework: the so-called Lyapunov's stability theory [34,35], which could be considered as the main theoretical reference on stability.One of the most useful mathematical procedures provided by this theory is the so-called direct method, consisting of the search for certain characteristic functions, called Lyapunov Functions, which provide sufficient conditions for the internal stability, simple or asymptotic, of a system.
Regarding the first of the considered alternatives-the NLMPC control systems-it is known that stability has been an essential and recurring matter in the evolution of the MPC.From its origin in the late 1970s, one of the main problems or disadvantages of the MPC was that stability could not be guaranteed (especially in the presence of constraints) and solutions were not introduced until the end of the 1990s and early year 2000, when the conditions to guarantee stability were clearly established [36,37].Throughout this period of time, numerous works were carried out on the stability of predictive control algorithms, with the greatest efforts being directed to the search and proposal of procedures that could guarantee closed-loop stability (predictive controllers with guaranteed stability).The use of restrictive terminal regions for states [38] and the consideration of infinite prediction horizons and inclusion of terminal penalty terms in the cost function [39][40][41][42] constitute some of the most important contributions and have been, at the same time, a reference for further lines of research in predictive controllers analysis and design, from then until now.In [3], a precise list of works on predictive controllers with guaranteed stability, carried out until 2002, is shown, including methods such as: MPC with terminal restriction of equality, MPC with terminal cost, MPC with terminal restriction of inequality, and MPC with terminal cost and terminal restriction.In addition, the author's proposals on stability and robustness of predictive controllers of nonlinear systems subject to constraints (with or without uncertainties) are also presented, based on the combination of Lyapunov's theory and invariant set theory.In [43], the advances in stability analysis of predictive control are also reviewed, within the framework of a general review or state of the art of predictive control.In [44], a brief but precise review of some methods on stability in predictive control existing at that time is included, specifically focused on the formulation of the MPC with guaranteed stability.In [45], a methodical and complete study of nonlinear model predictive control (theory and algorithms) is carried out, treating the stability aspects with an amplitude and intensity consistent with their importance.In addition, [9] reviews the state of the art of predictive control today, with special emphasis on the main objective of research works on predictive control, which is (as in the early 2000s) to ensure stability and robustness of the closed-loop system in a compatible way with the optimization approach.In all the referenced works, the use of Lyapunov's theory (direct method) tries to demonstrate that the cost function (the minimization of which allows us to calculate the optimal control signal) is a Lyapunov function.However, all the proposed closed-loop stabilization methods are not directly applicable to our FMBPC control strategy, since in it the control variable is not determined from a cost function.It is, therefore, necessary to analyze other alternatives.
Regarding the second of the alternatives considered-FS systems-it must be noted that the stability analysis of control systems in such a frame or field is not a trivial task, due, on the one hand, to the intrinsic nonlinear structure of this type of systems and, on the other, to the great diversity of modalities of use of fuzzy logic for control purposes.In relation to the latter and among other possible classifications, we could differentiate between controllers based on fuzzy plant models and controllers based on fuzzy rules (the plant models being fuzzy or not).Our strategy, FMBPC, fits the first type, but even if we focus only on that category, there can be, in principle, many possible ways to approach stability analysis, since there is no single method of stability analysis for FS systems.Many lines of research have been developed and numerous works with different approaches have been presented.Many of the works carried out have consisted of deducing sufficient (but not necessary) stability conditions.Furthermore, due to the great difficulty of generalization, many of the results have been developed, or tested, only for specific case studies.In short, there is a great diversity of proposals.Nevertheless, a large part of them share the theoretical framework of reference, which is, of course, Lyapunov's stability theory, widely used also by many other control strategies, as we have already commented previously.In [46], an exhaustive review of the state of the art of FS systems stability analysis is carried out (until 2006), highlighting its difficulty due to the non-linearity characteristic of this type of systems, and focusing the study on global and asymptotic stability.Different stability criteria and various techniques based on Lyapunov's stability theory are reviewed (considering quadratic, but also fuzzy Lyapunov functions) and a specific review of FS systems' stability using TS fuzzy models is made.Likewise, the approach of Linear Matrix Inequalities (LMIs) [47,48] is analyzed.In [49], there is also a review of several works and results related to stability analysis and control algorithms design for nonlinear systems represented by TS fuzzy models.In such works, the general procedure used to address stability and derive stabilization results is the Lyapunov's direct method.The stability conditions and the control design problem are expressed by means of constraints in the form of LMIs.Additionally, regarding the stability analysis, most of the methods studied in the mentioned review use a quadratic Lyapunov function (quadratic stability) and obtain, as a result of the demonstrations, only sufficient stability conditions (stability does not necessarily imply quadratic stability).With regard to stabilization methods, the review focuses mainly on those based on state feedback, although other alternative methods, such as those using output feedback, are also discussed.Naturally, in addition to these two great references that have just been described (where several criteria, techniques and approaches are reviewed and numerous authors are cited), there are many other works on the stability of FS systems.Thus, in the following list of works (not exhaustive, but extensive), we can see other contributions and approaches, some of them very significant: [11,.In addition, some more specific references, relative to control systems that use fuzzy models (among them, FMBPC control systems) and that include (to a greater or lesser degree) stability analysis, or stability-based design, may be, among others, the following: [11][12][13]22,[83][84][85][86][87][88][89][90][91][92][93].
Our stability analysis could follow, in principle, any of the multiple lines of work or approaches mentioned in the two fields analyzed.However, regardless of the approach, FMBPC control systems' stability analysis is not an easy task.For this reason, our choice consisted in looking for a practical computational approach, which is less complex, developed for our specific case study, but with theoretical possibilities of generalization.Our approach consists, basically, of starting from the identified model of the system, a TS fuzzy discrete model, expressing it in the state space in discrete linear time-varying (DLTV) form, and subsequently obtaining an approximate local model of type discrete linear time-invariant (DLTI) for states located in the vicinity of a certain operating point, and applying the first Lyapunov's stability method or theorem, demonstrating its (local) fulfillment computationally.As has been already previously commented, the stability analysis of control systems based on fuzzy models, starting from the original structure of these in the form of rules, is very complicated.However, the use of equivalent LTV state-space models could dramatically simplify the stability analysis task.The stability theory for state-space linear systems is well established and there are numerous works, both for continuous-time and for discrete-time.In [94][95][96][97][98] can be found the basic theoretical foundations, and in [99,100], two specific works on stability analysis of LTV systems of continuous-time and discrete-time, respectively, can be seen.The work presented in [101] can also be consulted, where a theoretical review is made regarding the so-called convex systems based on Takagi-Sugeno, linear parameter varying (LPV), or quasi-LPV models.LPV systems can be considered a generalization of certain LTV systems and therefore, the stability analysis carried out in the said paper could be useful in some way.In any case, the work developed in [100] cited above would be more suitable for our modeling approach and could be a good reference to study the stability corresponding to our case study, after having converted the TS fuzzy model of the process to an equivalent DLTV model.However, and although some necessary and sufficient stability conditions are derived in that specific work, most of the conclusions of this type of study refer only to sufficient stability conditions.Furthermore, the mathematical complexity of the proposed methods in this type of work does not facilitate their use, in a simple and methodical way, in determining the stability of other different case studies.This justifies our attempt to solve the problem by making use of an equivalent local DLTI model (valid for states close to a certain steady state) and using a computational method to analyze stability.
In Section 2, the chosen case study, the ASP wastewater treatment biological process, is presented.Furthermore, the plant identification procedure, which produces a discrete TS fuzzy model, as a result, is revised and this model is converted into another equivalent in the state space, of DLTV type.In Section 3, a local state-space DLTI model is deduced and using such a model, open-loop stability is analyzed using the appropriate criteria for DLTI systems.Section 4 is dedicated to describing the control law and the closed-loop stability criteria are derived.In Section 5, using a symbolic computational approach and an induction process, the closed-loop stability conditions are obtained as a function of the controller parameters and open-loop dynamics of the selected case study.In Section 6, several numerical results are presented: first a stability test and then a series of FMBPC control experiments performed by simulation.Finally, in Section 7, the conclusions are presented, together with some future research options.Complementarily, Appendix A contains detailed numerical information on the membership functions associated with the identification of the fuzzy models of the ASP process and Appendix B includes a complementary mathematical development of Section 4.

The Activated Sludge Process: Takagi-Sugeno Fuzzy Modeling and Representation in State Space
In this article, we chose as case study the ASP biological purification procedure, commonly used in wastewater treatment plants.This process has a complex dynamic and is highly nonlinear due to its biological nature.The identification of a fuzzy model of the process from input-output data is a good choice to represent the behavior of this type of processes (as already noted in the Introduction).Initially, the resulting models of the identification procedure are Takagi-Sugeno type fuzzy models, but through the appropriate formalization they will be expressed finally as equivalent state space models, in the form of DLTV models.This section describes the process chosen as a case study, as well as the procedure for the fuzzy identification of the said process and the representation of the resulting TS model in state space, in the form of a DLTV model.

The Activated Sludge Process: Description and Classical Mathematical Model
In our case study, the basic structure for the WWTP (an aerated biological reactor followed by a secondary settler) was considered.The corresponding purification process consists, essentially, of the elimination of organic contaminants (substrate) through a culture of bacteria (biomass) that will feed precisely on that substrate.This biological process takes place in the aerated and stirred reactor and the treated water is sent to a secondary decanter in which the clean water separates from the resulting sludge.This sludge is sedimented at the bottom of the decanter and sent back to the reactor (sludge recirculation process).The complete cycle is called the activated sludge process.The mass balance of the reactions that take place in this process can be described using the classical Monod model.Based on this model, together with the hypothesis of a perfectly mixed reactor, in [23], a certain mathematical model is specified to represent the activated sludge process.The adopted model, which can be considered as a simplification of the well-known Activated Sludge Model No. 1 (ASM1) [102], consists of a set of differential and algebraic equations that relate the temporal evolution of the masses of the substrate, biomass, and oxygen.From the analysis of the mathematical model, either the complete or the simplified model, it can be deduced that it is highly non-linear, and this is precisely one of the main reasons why it is interesting to adopt it as a case study, in the context of a control problem.Another important reason why it is an attractive system as a case study is that it is a difficult system to control (due to its biological nature).In addition, it is multivariable.
Our case study has focused on a sub model of the adopted model (mentioned above), consisting of the equations that describe the interactions of five variables: three inputs (one manipulable variable and two disturbances) and two outputs (to be controlled both, simultaneously).Figure 1 shows schematically the different inputs and outputs of the biological process taken into account, as well as its characterization.The variables to be controlled simultaneously are two: the substrate concentration in the effluent, s, and the biomass concentration in the reactor, x.The inputs variables are three: one manipulated variable, the sludge recirculation flow-rate, q r , and two disturbances, the wastewater flowrate in the input of the plant, q i , and the substrate concentration in the input of the plant, s i .This set of inputs and outputs, and its role, is more clearly detailed in Table 1.
the reactions that take place in this process can be described using the classical Monod model.Based on this model, together with the hypothesis of a perfectly mixed reactor, in [23], a certain mathematical model is specified to represent the activated sludge process.The adopted model, which can be considered as a simplification of the well-known Activated Sludge Model No. 1 (ASM1) [102], consists of a set of differential and algebraic equations that relate the temporal evolution of the masses of the substrate, biomass, and oxygen.From the analysis of the mathematical model, either the complete or the simplified model, it can be deduced that it is highly non-linear, and this is precisely one of the main reasons why it is interesting to adopt it as a case study, in the context of a control problem.Another important reason why it is an attractive system as a case study is that it is a difficult system to control (due to its biological nature).In addition, it is multivariable.
Our case study has focused on a sub model of the adopted model (mentioned above), consisting of the equations that describe the interactions of five variables: three inputs (one manipulable variable and two disturbances) and two outputs (to be controlled both, simultaneously).Figure 1 shows schematically the different inputs and outputs of the biological process taken into account, as well as its characterization.The variables to be controlled simultaneously are two: the substrate concentration in the effluent, , and the biomass concentration in the reactor, .The inputs variables are three: one manipulated variable, the sludge recirculation flow-rate,  , and two disturbances, the wastewater flowrate in the input of the plant,  , and the substrate concentration in the input of the plant,  .This set of inputs and outputs, and its role, is more clearly detailed in Table 1.

Biomass concentration in the reactor 𝑥
In our work, the municipal WWTP of Manresa (Barcelona, Spain), operating since 1985 (with an evolution of the population of the urban area of the city from about 67,000 people in 1985 to about 76,000 people today), was considered as a reference plant.The dynamics of the ASP process of this treatment plant, in its initial period of operation, is acceptably well represented by the model and sub-model considered above and, therefore, it is coherent to make use of real data from such a plant for our studies.Thus, various Substrate concentration in the influent s i

Sludge recirculation flow-rate q r
Substrate concentration in the effluent s

Biomass concentration in the reactor x
In our work, the municipal WWTP of Manresa (Barcelona, Spain), operating since 1985 (with an evolution of the population of the urban area of the city from about 67,000 people in 1985 to about 76,000 people today), was considered as a reference plant.The dynamics of the ASP process of this treatment plant, in its initial period of operation, is acceptably well represented by the model and sub-model considered above and, therefore, it is coherent to make use of real data from such a plant for our studies.Thus, various series of input and output data recorded in such a period, in different industrial campaigns of this WWTP [103], were used (or taken into account for reference) to guide the fuzzy identification procedures of the simulated ASP process and also to adjust, in the control experiments implemented by means of simulation, the ranges of values of the process inputs.Both in the identification procedures and in the control experiments, a part of the numerical values used for the inputs correspond, or are related, to the data collected in the campaigns of the reference WWTP, as follows: for the disturbances, typical values (real or weighted), and for the manipulated variable, similar values to the real values at the different operating points.

Fuzzy-TS Identification of the ASP Process
The identification of the ASP process was carried out from input-output numerical datasets, with the inputs properly selected, taking as reference data recorded in industrial campaigns of the WWTP of Manresa.The outputs were obtained in simulation, using the mathematical model chosen for the treatment plant (specified in the previous subsection).The result of the identification carried out is a multivariable TS fuzzy model.
The tasks and procedures of the fuzzy identification were performed using the software tool known as FMID (Fuzzy Model Identification Toolbox) [104], complemented with some adaptations made within the framework of this work, programmed with Matlab ® & Simulink ® .The FMID tool was developed to support the fuzzy modeling and identification techniques described in [11].The fuzzy identification functions of this tool allow the implementation of clustering techniques based on the Gustafson-Kessel algorithm.The numerical data of input-output provided to the FMID tool for the identification and validation of the process were composed as follows: for the input data, as we said in the previous subsection, typical values or similar to the real values of the WWTP of Manresa; as for the output data, they were obtained by simulation in open loop, with the mentioned inputs and with the plant represented by its continuous model in the form of differential equations.In the numerous identification tests carried out, the various available data series were generally divided into two subsets: identification data and validation data (although, for comparison purposes, in part of the tests, the same data set was used for both identification and validation).Our case study is a multivariable system (MIMO system), with three inputs and two outputs.This multivariable architecture can be seen schematically in Figure 2: control experiments implemented by means of simulation, the ranges of values of the process inputs.Both in the identification procedures and in the control experiments, a part of the numerical values used for the inputs correspond, or are related, to the data collected in the campaigns of the reference WWTP, as follows: for the disturbances, typical values (real or weighted), and for the manipulated variable, similar values to the real values at the different operating points.

Fuzzy-TS Identification of the ASP Process
The identification of the ASP process was carried out from input-output numerical datasets, with the inputs properly selected, taking as reference data recorded in industrial campaigns of the WWTP of Manresa.The outputs were obtained in simulation, using the mathematical model chosen for the treatment plant (specified in the previous subsection).The result of the identification carried out is a multivariable TS fuzzy model.
The tasks and procedures of the fuzzy identification were performed using the software tool known as FMID (Fuzzy Model Identification Toolbox) [104], complemented with some adaptations made within the framework of this work, programmed with Matlab ® & Simulink ® .The FMID tool was developed to support the fuzzy modeling and identification techniques described in [11].The fuzzy identification functions of this tool allow the implementation of clustering techniques based on the Gustafson-Kessel algorithm.The numerical data of input-output provided to the FMID tool for the identification and validation of the process were composed as follows: for the input data, as we said in the previous subsection, typical values or similar to the real values of the WWTP of Manresa; as for the output data, they were obtained by simulation in open loop, with the mentioned inputs and with the plant represented by its continuous model in the form of differential equations.In the numerous identification tests carried out, the various available data series were generally divided into two subsets: identification data and validation data (although, for comparison purposes, in part of the tests, the same data set was used for both identification and validation).Our case study is a multivariable system (MIMO system), with three inputs and two outputs.This multivariable architecture can be seen schematically in Figure 2: Consistent with such architecture, the input-output datasets required for the FMID tool must be arranged as arrays with as many rows as samples and with five columns: three columns for the inputs and two columns for the outputs.
The identification tool provides TS fuzzy models (as well as mechanisms for their validation), expressed in the form of logical rules of the if-then type, that is: if (antecedent) then (consequent), where the consequent is a linear combination of the components of the consequent vector, plus a constant term (affine function).In Equation (1), we show the structure of a generic rule,   (j-th rule), corresponding to a certain output, (), of the identified process: Consistent with such architecture, the input-output datasets required for the FMID tool must be arranged as arrays with as many rows as samples and with five columns: three columns for the inputs and two columns for the outputs.
The identification tool provides TS fuzzy models (as well as mechanisms for their validation), expressed in the form of logical rules of the if-then type, that is: if (antecedent) then (consequent), where the consequent is a linear combination of the components of the consequent vector, plus a constant term (affine function).In Equation (1), we show the structure of a generic rule, R j (j-th rule), corresponding to a certain output, y(k), of the identified process: where: is the consequent vector • A jp is the fuzzy set (or fuzzy value) associated with the p-th component of the antecedent vector (for j-th rule); this fuzzy set will associate a certain membership function µ A jp (x) : R → [0, 1] , which, applied to any value of the variable x ap , will determine its membership grade with respect to the fuzzy set A jp , that is: µ A jp x ap (for each component of the antecedent vector, there will be as many membership functions as rules are) • k is the k-th instant (discrete-time system) • φ j (x) is an affine function of x • α jq is the coefficient, in the affine function, of the q-th component of the consequent vector (for j-th rule) • δ j is the independent term or offset (in the affine function) Each rule of the TS fuzzy model is equivalent to a linear sub-model associated with a certain region of the input space of the model, in a greater degree than the sub models of the other rules.The output of each sub model is calculated by the corresponding affine function.In addition, the global output of the TS fuzzy model will be calculated by combining or aggregating the outputs of all the rules in a weighted way, according to the degree to which the antecedent vector fulfills or satisfies the premise or antecedent of each of the rules.This calculation will be called global weighted output and will be represented by y(k).In the fuzzy logic literature, there are several methods for combining or aggregating rules; among them is the well-known centroid method.Applying such a method to the fuzzy model represented in (1), for a general case with multiple outputs, the weighted global output can be expressed as shown in (2) (see [23]), where y i (k) (i = 1, 2, . . ., l) is the expression corresponding to the global i-th output (having omitted, for the antecedent and consequent vectors, x a and x, respectively, and also for their respective components, the i index, i.e., indication of the dependence with respect to the i-th output, except in the dimension of the vectors (p i and q i , respectively), in order to simplify the notation; for the same reason, the express indication of the dependency of the two vectors with respect to k, has also been discarded): where: • y i (k) is the i-th weighted global output • k is the k-th instant (discrete-time system) • i = 1, 2, . . ., p (number of outputs) • j = 1, 2, . . ., mr i (number of rules of fuzzy model for output y i ) • β ij (x a ) are the normalized membership functions of the antecedent vector, x a (for j-th rule of i-th output) • φ ij (x) is the affine function of the consequent vector, x (for j-th rule of i-th output) In the specific case of our ASP process, the output variables considered are two: y 1 (k), the substrate in the effluent, and y 2 (k), the biomass in the reactor.Therefore, each of the fuzzy identification tests carried out led to a fuzzy model constituted by two TS fuzzy submodels (one TS model for each output).In the present work, the results obtained in two of the many identification tests carried out (see [23]) were selected, denoting the two obtained fuzzy models as FM 1 and FM 2 , respectively.The criteria for their selection were two: on the one hand, the dynamic structure of the models and on the other, the identification validation index known as VAF (percentile Variance Accounted For between two signals).Thus, two different dynamic structures were chosen (see Dynamical Parameters in Table 2) to test the possible influence of this structural aspect on the performance of the control strategy, and from among the identified models, two models were chosen, one for each structure, with sufficiently high VAF indices.Specifying more, each of these fuzzy models is characterized by the choice (in the design phase) of different values of several identification parameters, which are shown below, in Table 2.Among them, N y , N u and N d (Dynamical Parameters of the recursive discrete model) have special relevance in the identification process and its results, since they contain implicit information related to the dependencies of the process outputs at any k-th instant with respect to the different inputs and outputs in previous time instants (the legend of Table 2 explains the particular meaning for the model FM 1 ; the meaning is analogous for FM 2 ): • where the N y , N u , and N d parameters and its concrete values, for FM 1 (for FM 2 it would be analogous), mean the following: • N y : output-output dynamic relationship matrix (2 × 2) N y o f FM 1 (row 1) mean that y 1 (k) depends on y 1 (k − 1) and ) mean that they is 1 transport delay for the three inputs (u 1 , u 2 and u 3 ) with respect to y 1 (k) N d o f FM 1 (row 2) mean that they is 1 transport delay for the u 1 and u 3 inputs with respect to y 2 (k) From the parametric specifications contained in Table 2 and with a suitable selection of the available input and output data, the chosen identification tool (FMID) will determine the two fuzzy models of the ASP process (FM 1 and FM 2 ).Each of these models will be characterized by the rules corresponding to the two TS sub-models, making it necessary to know, for each rule, the fuzzy sets associated with each of the antecedent vector components (that is, the identifiers of those sets and their corresponding membership functions), and the affine function of the consequent, φ ij (x) (or equivalently the numerical coefficients that multiply the components of the consequent vector, and the independent term or offset).The structure of the models identified by using the FMID tool, FM 1 and FM 2 , their rules, and the corresponding fuzzy and numerical elements that have just been mentioned, are shown in detail in different tables.For the FM 1 model, Table 3 show the fuzzy sets of the antecedent and Table 4 the affine function of the consequent; and for the FM 2 model, Table 5 presents the fuzzy sets of the antecedent, and Table 6 shows the affine function of the consequent.In all cases, the input and output variables of the ASP process are the ones shown in Table 1.It can be observed that all these tables contain information corresponding to the two outputs of the ASP process (the two TS submodels), but a different notation has been used for each one.Thus, in the case of the FM 1 model, the rules are represented by por R j , for y 1 (k), and R * j , for y 2 (k); as for the fuzzy sets of the antecedent vector, they are represented by A jp , for y 1 (k), and A * jp , for y 2 (k).In the case of the FM 2 model, the established notation differences are analogous (R j and R * j , for the rules of y 1 (k) and y 2 (k), respectively, and A jp and A * jp , for the fuzzy sets of the antecedent vector of y 1 (k) e y 2 (k), respectively).Finally, precise and detailed information related to the membership functions µ A jp (x) defining each fuzzy set A jp (using a general notation), has been specified, given its extension, in Appendix A (Tables A1 and A2), in parametric form.In our case, membership functions are piece-wise exponential type functions.Table 3. Fuzzy sets of the antecedent in the rules of the identified FM 1 fuzzy model.

Components of the Antecedent Vector, x a , and Corresponding Fuzzy Sets ASP Output
Rule where: • A uv is the fuzzy set (or fuzzy value) corresponding to the v-th component of the antecedent vector, for the u-th rule, R u , of the output y 1 (k) (substrate); the corresponding membership functions are of piece-wise exponential type.

•
A * wz is the fuzzy set (or fuzzy value) corresponding to the z-th component of the antecedent vector, for the w-th rule, R * w , of the output y 2 (k) (biomass); the corresponding membership functions are of piece-wise exponential type.

•
The membership functions corresponding to the fuzzy sets A uv and A * wz are provided by FMID tool in parametric form; these parameters, as well as the mathematical expression of this type of membership functions, are detailed in Appendix A (Table A1).
Table 5. Fuzzy sets of the antecedent as per the rules of the identified FM 2 fuzzy model.

Components of the Antecedent Vector, x a , and Corresponding Fuzzy Sets ASP Output
Rule where: • A uv is the fuzzy set (or fuzzy value) corresponding to the v-th component of the antecedent vector, for the u-th rule, R u , of the output y 1 (k) (substrate); the corresponding membership functions are of piece-wise exponential type.

•
A * wz is the fuzzy set (or fuzzy value) corresponding to the z-th component of the antecedent vector, for the w-th rule, R * w , of the output y 2 (k) (biomass); the corresponding membership functions are of piece-wise exponential type.

•
The membership functions corresponding to the fuzzy sets A uv and A * wz are provided by FMID tool in parametric form; these parameters, as well as the mathematical expression of this type of membership functions, are detailed in Appendix A (Table A2).Table 6.Coefficients of the consequent in the rules of the identified FM 2 fuzzy model, and offset terms.

Components of the Consequent Vector, x, and Corresponding Coefficients ASP Output
Rule The above tables specifically show the dynamic relationship between the outputs of the ASP process and the consequent vector x, consistent with the configuration defined by the values of N y , N u , and N d parameters shown in Table 2. Furthermore, it can also be observed in the above tables that, in the fuzzy TS models obtained by using the used identification tool, the antecedent vector is coincident with the original consequent vector, i.e., x a = x.All these relationships or dependencies can also be functionally expressed as described below.
Thus, in the case of the FM 1 model, the observed dependencies for the output y 1 (k), and for the global output y 1 (k), which is a weighted linear combination of the outputs of all rules (Equation ( 2)), as well as the components of the corresponding antecedent and consequent vectors, x a = x| (FM 1 , out 1 ) , are the following (Equation ( 3)): Regarding the output y 2 (k) and the global output y 2 (k), as well as the components of the corresponding antecedent and consequent vectors, x a = x| (FM 1 , out 2 ) , the observed relationships are the following (Equation ( 4)): Note that from Equations ( 3) and ( 4) can be observed that the consequent vector corresponding to the y 1 (k) output is not exactly the same as the consequent vector corresponding to the y 2 (k) output, but the components of the first vector include the components of the second vector.We will consider then the set formed by the union of the components of both vectors and we will form a consequent vector common to both outputs, x| FM 1 , to be able to mathematically deal together the fuzzy models of both outputs, within a matrix mathematical framework (by considering, logically, null factors where necessary in the numerical expressions related to output y 2 (k)).Equation ( 5) summarizes the functional relationships of the two weighted global outputs, with respect to the common consequent vector, x| FM 1 : where: • f 1c (•) and f 2c (•) are the trivial adaptations of f 1 (•) and f 2 (•), respectively, after the definition of the common consequent vector, x| FM 1 • x| FM 1 is the consequent vector, common to both outputs For the FM 2 model, after a process of observation and deduction analogous to that of the FM 1 model, we observe the following dependency relationships for the two outputs of the ASP process, as well as the corresponding consequent vectors (Equations ( 6) and ( 7)): Again, as in the case of the previous model, we observe that the consequent vectors of both outputs are not coincident, but given that the one corresponding to the output y 1 (k) includes among its components those of the one corresponding to the output y 2 (k), we can define a common consequent vector that is the union of both, x| FM 2 (considering null factors in the numerical expressions related to output y 2 (k) where necessary).Equation (8)  summarizes the functional relationships of the two weighted global outputs, with respect to the common consequent vector, x| FM 2 : (8) where: • f 1c (•) and f 2c (•) are the trivial adaptations of f 1 (•) and f 2 (•), respectively, after the definition of the common consequent vector, x| FM 2 • x| FM 2 is the consequent vector, common to both outputs Generalization: Finally, and given that the components of the consequent vector x| FM 1 include those of the consequent vector x| FM 2 , or equivalently, those of the latter are a subset of those of the first, we can consider the union of both as a consequent vector common to both models (in our case, it would be coincident with the first), representing it as x.In this way, we can express the functional dependency relationships, for either of the two models, with respect to a single common consequent vector, x, naturally considering null factors, where necessary, in the numerical expressions corresponding to the FM 2 model.The summarized general expressions, formally valid for both models, would be the following (Equation ( 9)): or, simpli f ying the notation (denoting to the f inal weighted outputs in the same way as to the individual rules outputs : where: • f a f 1 (•) and f a f 2 (•) are the affine functions corresponding to the outputs y 1 (k) and y 2 (k), respectively (the coefficients and offset term of these affine functions depend on the fuzzy model (FM 1 or FM 2 )) • x is the consequent vector, common to both fuzzy models (and common to two outputs) Remark: Both the outputs of the process, as well as the different antecedent and consequent vectors, have been expressed in the original form used by the FMID software tool, concerning temporal dependencies.Thus, the expressions corresponding to the outputs at the k-th time instant have been explicitly expressed.In order to specify the expressions of the outputs at the (k + 1)-th time instant, as well as those of the corresponding antecedent and consequent vectors, we will simply have to substitute k by (k + 1) in the expressions we have shown previously.Thus, for example, in the case of the FM 1 model, the antecedent and consequent vectors corresponding to the outputs y 1 (k + 1) and y 2 (k + 1), as well as the fuzzy rules, in their generic form, would be as is detailed below (before adopting a common consequent vector): FM 1 model (outputs at the (k + 1)-th time instant).

State-Space DLTV Equivalent Model
Making use of the necessary mathematical basis (mainly, basic matrix calculation) and by an adequate treatment, the TS fuzzy models obtained can be transformed or expressed in the state space, in the form of discrete, linear, and time-varying models, that is, as DLTV type models.
We will start from Equation (2), the mathematical expression that allows us to calculate the outputs of the process (at every k-th time instant) from the knowledge of the membership functions of the antecedent and of the numerical coefficients of the consequent, obtained as a result of the identification of fuzzy models (FM 1 or FM 2 , in our case study).If we develop Equation (2) for a generic TS fuzzy model (of the type identified), for each of the two global outputs of the ASP process at the k-th time instant (that is, returning to the form provided by the software tool FMID), explicitly showing the affine function (as a linear combination of the components of the consequent vector, plus the independent term), we would obtain the following two equations (where, so as not to overload the notation, as in Equation (2), it has been omitted again the formal indication of the dependency of x a with respect to the i index, that is, with respect to the i-th output, and also the dependency of x a and x with respect to k): By introducing now suitable mathematical definitions (mr y β j (x a )), we can group Equations ( 11) and ( 12) into a single matrix equation: , where : is the normalized membership f unction o f the antecedent vector corresponding to the f irst output (y 1 ), f or the j-th rule β 2j (x a ) is the normalized membership f unction o f the antecedent vector corresponding to the second output (y 2 ), f or the j-th rule • β 26 (x a ) = 0; α * 6i ≡ 0 (i = 1, . . ., 6), δ * 6 ≡ 0(there are only 5 rules f or y 2 ) and by substituting k for (k + 1) in Equation (13), that is, considering the outputs at the (k + 1)-th time instant (see (10)), we will obtain the equivalent Equation ( 14) (having mr y β j (x a ) the same or analogous meaning, respectively, as they have in Equation (13), and where x a represents, in β 1j (x a ), the corresponding antecedent vector to the output y 1 (k + 1), and in β 2j (x a ), x a represents the corresponding antecedent vector to the output y 2 (k + 1); all these meanings and the numerical peculiarities for j = 6 shown in (13) will be kept from now on): Next, to improve the characterization of the different terms of Equation ( 14), we will introduce the following notation changes: Now, rearranging the terms of Equation ( 14) (exchanging the positions of the second and third terms) and applying the notation changes introduced in (15), the following equation will be obtained: ), r * 6 ≡ 0 (there are only 5 rules f or y 2 ) As a further step, in the process of development and transformation of the original Equation ( 2), with integration of the two outputs, we will make the following vector and matrix definitions and formalizations (where the particularities corresponding to j = 6, indicated in (16), are maintained): where: • k represents (k•T) and T is the sampling period • z mN (k) is the definition of state vector at the k-esimo time instant (which groups the two outputs of the process at the k-th time instant) • u a (k) is the input vector or the extended input vector (with two components: the manipulated process variable at the k-th time instant and the manipulated variable at the(k − 1)-th time instant, respectively) • d(k) is the disturbances vector (which groups the two unmanipulated process variables at the k-th time instant) • y mN (k) is the output vector (that matches the state vector definition) and where the different input and output process variables involved (whose physical meaning can be seen in Table 1) are the following: Finally, by substituting in Equation ( 16) the definitions and formalizations made in (17), we will have the following state equation, where the time dependence of the antecedent vector, x a , (x a = x a (k)) has been shown explicitly due to its importance in the model characterization: On the other hand, taking into account that ∑ mr i j=1 β ij (x a ) = 1, for any i (see Equation ( 2)), it is then true that ∑ mr j=1 β 1j (x a ) = 1 and ∑ mr j=1 β 2j (x a ) = 1 (since mr = max.(mr 1 , mr 2 ) and β 26 (x a ) = 0).Taking into account these last two equalities and grouping the outputs y 1 (k) and y 2 (k) in a single vector, it follows: and doing now the following vector and matrix definitions and formalizations: where: • k represents (k•T) and T is the sampling period • y mN (k) is the output vector (which groups the two outputs of the process at the k-th time instant and matches the state vector definition) and by substituting Equation (20) in Equation ( 19), we will have the following output equation, where the time dependence of the antecedent vector, x a , (x a = x a (k)) has been shown explicitly (as in the state equation) due to its importance in the model characterization: Next, we will represent by means of a single matrix each of the matrix sums of the different terms of Equations ( 18) and ( 21), according to the following definitions: The matrices defined in (22) can also be expressed in the following symbolic and generic way (which will be useful later to simplify the stability analysis), representing the elements of the matrices in a simplified way, taking into account the dynamic characteristics of each one of the models (see the identification parameters N y , N u and N d , in Table 2), and by specifying its relationship with the original coefficients of the starting fuzzy models, FM 1 and FM 2 : (1) f or both f uzzy models, a * j1 = 0, ∀j; then : ( f or FM 1 ; f or FM 2 : c = 0 and g = 0 (4)) (3) f or both f uzzy models, d * j2 = 0, ∀j; then : Finally, and independently of the possible formal matrix representations shown in ( 23)), if we substitute the matrices defined in (22) in Equations ( 18) and ( 21), and we consider both equations together, we will obtain the final equivalent state space model (Equations ( 24) and ( 25)), which is of DLTV type, as we will reason below: where: • k represents (k•T) and T is the sampling period • z mN (k) is the state vector • u a (k) is the input vector (or the extended input vector) • d(k) is the disturbances vector • y mN (k) is the output vector • A mN , B mN , D mN , R mN , and C mN are the system matrices As can be seen in the equations of the obtained model, the state Equation ( 24) and the output Equation ( 25), the final equivalent model is a state-space discrete model, with input disturbances, and formally linear.However, observing Equation (22), it is immediate to deduce that the system matrices (A mN , B mN , D mN , R mN , and C mN ) depend on time, since they depend on β j (x a (k)), and x a (k) depends directly of k.That is: Therefore, the equivalent state-space model just described is time-varying and, more specifically, it is of DLTV type (which is what we want to demonstrate).Furthermore, as a global conclusion of the transformation process, we can highlight that the highly non-linear system under study (the ASP process), initially identified through fuzzy models, has finally been represented by a model in state space, with a linear formal appearance, although with varying coefficients over time.
Remark: The dependence of the system matrices, A mN , B mN , D mN , R mN , and C mN , with respect to the antecedent vector x a (k), and therefore also with respect to time, implies that it is not possible to characterize the fuzzy models of the ASP process using a set of unique constant numerical matrices (as is the case of time-invariant linear systems), except for small environments around certain operating points or for small time intervals, in the case of sufficiently slow systems.On the other hand, as a consequence of this dependence of the system matrices with respect to x a (k), in each iteration of the simulation experiments with the FMBPC control strategy, it is necessary to carry out the corresponding update sequence.That is, at the k-th time instant, both x a (k) and β j (x a (k)) must be updated, as well as, finally, the system matrices, A mN (k), B mN (k), D mN (k), R mN (k), and C mN (k).
State-space DLTV model of ASP process, in matrix form and with the elements of all matrices expressed in a generic symbolic way: Replacing in ( 24) and ( 25) the state, input, disturbances and output vectors by their corresponding expressions (Equations ( 17) and ( 20)) and the system matrices by their corresponding generic symbolic expressions (defined in ( 23)), we would have, for the specific case of our identified ASP process, the state-space DLTV model in matrix form and with the elements of the system matrices expressed in a generic way, as shown below (which will be useful, as we already said in the paragraph before Equation (23), to simplify the stability analysis of the system and its interpretation): y mN (k) this formalization being valid for both FM 1 and FM 2 (with c = 0 and g = 0, for FM 2 ).

Open-Loop Local Stability
We will now assume that for the model represented by Equations ( 24) and ( 25), there will be some equilibrium point, that is, we will assume that there will be some steady state, z mNss , for a certain input, u ass , for certain values of the disturbances, d ss , for a certain antecedent vector x ass and consequently, for certain matrices A mNss , B mNss , D mNss , R mNss and C mNss and with a certain output, y mNss .Such a state must satisfy the steady state condition for a discrete-time system, i.e.: (28) and on the other hand, it must also satisfy the equations of the model (( 24) and ( 25)), that is: y mNss = C mNss z mNss (30) Let us now consider a generic state, z mN (k), belonging to the model represented by ( 24) and ( 25) and suppose that the disturbances, d(k), are the same or very similar to those corresponding to the steady state, that is: and let us also suppose that the antecedent vector, x a (k), is close enough to x ass (in its corresponding space), to be able to consider that the system matrices are approximately equal to those of the steady state, that is: The assumption made, x a (k) close to x ass , would imply the following, considering the more general case of the antecedent corresponding to the output y 1 (k + 1) of the FM 1 model (for the others cases, it would be analogous): for this case, (17)); consequently, the as- sumption x a (k) ∼ = x ass would imply that all variables (state variables, disturbances and the extended input) should be close to the corresponding values associated to the considered stationary state, i.e.,: z mN (k) ∼ = z mNss , d(k) ∼ = d ss (hypothesis already considered in (31)), and u a (k) ∼ = u ass .We will assume the fulfillment of these conditions in our case study (for all cases) and will refer to them, in abbreviated form, simply as: proximity condition of state z mN (k) with respect to the steady state z mNss .
Subtracting now the Equalities ( 24) and ( 29), on the one hand, and ( 25) and (30), on the other, we will obtain the equality relations that are shown in the following equations, by also considering some groupings of terms and some numerical considerations derived from the hypothesis mentioned above (d(k) − d ss ∼ = 0 y R mN − R mNss ∼ = 0): Processes 2021, 9, 531 20 of 50 (y mN (k) − y mNss ) The Equations ( 33) and ( 34) describe an incremental model, valid for small deviations with respect to a certain operating point, certain steady state taken as a reference (as just explained above), that can be simplified and expressed more compactly as follows (using the hypotheses and the complementary notation included in the same expressions (33) and ( 34)): and making the following change of notation: the incremental model will be finally expressed in the following standard form: which constitutes a state-space local model, valid for states close to zero state, of DLTI type (considering the matrices A mN , B mN and C mN as constants, approximately equal to the corresponding matrices associated with the steady state, z mNss ).
Open-loop local stability: The form of the approximate (local) model that we have just deduced allows us to directly apply the more known and accepted stability criterion (necessary and sufficient condition) for DLTI systems described in the state-space: a DLTI system is asymptotically stable in the sense of Lyapunov (internal stability) if and only the eigenvalues of the state matrix are strictly within the unit circle.Therefore, taking into account that the state matrix corresponding to our case, A mN , is triangular (see Equation ( 23)), the system described by Equations ( 38) and (39), which we will denote by S OL below, will be stable if and only the absolute value of all the elements of the main diagonal of A mN is strictly less than unity, that is: S OL is asymptotically stable the eigenvalues o f A mN are all strictly within the unit circle → f or particular cases with triangular state matrix, A mN = a b 0 f , we would have : and in our case study (both f or FM 1 and FM 2 ), replacing the elements o f the main diagonal o f the A mN matrix by their respective original expressions, based on the original coe f f icients o f the identi f ied Fuzzy Models (see ( 23)), we would have :

Closed-Loop Local Stability of ASP Process Controlled by FMBPC Law
The objective of this section is to analyze the closed-loop local stability of our ASP process, making use of a certain FMBPC control strategy.We will start from the equations of DLTI local model developed in Section 3 (Equations ( 38) and ( 39)) and will replace in them the FMBPC control law deduced in [23], thus obtaining the corresponding closed-loop model.

FMBPC Control Law
In [23], an analytical and explicit FMBPC control law was deduced following a method that could be considered an extension to the multivariable case of the well-known PFC strategy (already mentioned in the introductory section of this article), one of whose main bases is the so-called equivalence principle.In Figures 3 and 4, we summarize the FMBPC strategy implemented in such work and the mentioned principle, respectively: 021, 9, x FOR PEER REVIEW 21 of 50

FMBPC Control Law
In [23], an analytical and explicit FMBPC control law was deduced following a method that could be considered an extension to the multivariable case of the well-known PFC strategy (already mentioned in the introductory section of this article), one of whose main bases is the so-called equivalence principle.In Figures 3 and 4, we summarize the FMBPC strategy implemented in such work and the mentioned principle, respectively:  The control law was derived by imposing on the outputs of the prediction model the following certain reference trajectories, along the so-called coincidence horizon (), and applying the equivalence principle to the end of that horizon.The result of the deductive process carried out in [23] was an analytical and explicit expression for the control action () (which is the sludge recirculation flow rate of the ASP process, i.e., () =  ()) and also, implicitly, the mathematical expression corresponding to the vector variable   () (reviewing the mathematical development and making some trivial considerations will suffice).This control law is the one that has been chosen to be implemented and obtain the closed-loop control system considered in this work.The original form of the law obtained in [23] is a matrix algebraic expression, deduced from a certain global model of predictions

FMBPC Control Law
In [23], an analytical and explicit FMBPC control law was deduced following a method that could be considered an extension to the multivariable case of the well-known PFC strategy (already mentioned in the introductory section of this article), one of whose main bases is the so-called equivalence principle.In Figures 3 and 4, we summarize the FMBPC strategy implemented in such work and the mentioned principle, respectively:  The control law was derived by imposing on the outputs of the prediction model the following certain reference trajectories, along the so-called coincidence horizon (), and applying the equivalence principle to the end of that horizon.The result of the deductive process carried out in [23] was an analytical and explicit expression for the control action () (which is the sludge recirculation flow rate of the ASP process, i.e., () =  ()) and also, implicitly, the mathematical expression corresponding to the vector variable   () (reviewing the mathematical development and making some trivial considerations will suffice).This control law is the one that has been chosen to be implemented and obtain the closed-loop control system considered in this work.The original form of the law obtained in [23] is a matrix algebraic expression, deduced from a certain global model of predictions The control law was derived by imposing on the outputs of the prediction model the following certain reference trajectories, along the so-called coincidence horizon (H), and applying the equivalence principle to the end of that horizon.The result of the deductive process carried out in [23] was an analytical and explicit expression for the control action u(k) (which is the sludge recirculation flow rate of the ASP process, i.e., u(k) = q r (k)) and also, implicitly, the mathematical expression corresponding to the vector variable u a (k) (reviewing the mathematical development and making some trivial considerations will suffice).
This control law is the one that has been chosen to be implemented and obtain the closed-loop control system considered in this work.The original form of the law obtained in [23] is a matrix algebraic expression, deduced from a certain global model of predictions in the state space, characterized by a certain extended state vector and certain system matrices.Consequently, the deduced control law is a function of those system matrices and of the aforementioned extended state vector.Therefore, in order to use such a law in the local model developed in Section 3 of this article (Equations ( 38) and ( 39)), must be first adapted.This local model derives from the global model in the state space obtained in Section 2.3 (Equations ( 24) and ( 25)), which is characterized by a state vector and system matrices different from those of the model used in [23].However, the two state space global models considered, the one deduced in [23] and the one deduced in Section 2.3 of this article, both of type DLTV, are actually equivalent, since both are transformations of TS fuzzy models of the ASP process (previously identified).Due to this equivalence and also taking into account that the extended state vector of the first model includes the state vector of the second, it will be possible to transform the expression of the control law into an equivalent one, which is a function of the matrices of the system and of the state vector corresponding to the global model deduced in this article.On the other hand, and prior to this adaptation, the expression corresponding to the extended manipulated input variable, u a (k), which is really implicit in the development done in [23], must be specified.The original expression of u a (k) and its subsequent adaptation have been detailed in Appendix B of this article.The resulting final expression for this variable, already adapted and therefore suitable to be substituted in our local open-loop model, is the following (see details of the matrices involved in Appendix B): , λ mN : matrix f unctions o f system matrices; γ : submatrix o f order 2;

Closed-Loop Local Stability Analysis
In the closed-loop stability study, we will also assume the existence of some equilibrium point for the open-loop plant model, given for Equations ( 24) and ( 25), and we will represent the corresponding steady state also by z mNss .We will consider that the objective of our control system is precisely to drive the system towards that state, whose output, y mNss , must therefore be the reference for the output of the closed-loop system.Furthermore, we will assume that the control strategy will achieve the system reaching this state and, consequently, the entry of the steady state input will be compatible with the control law.Finally, we will suppose that the model described in Equations ( 24) and ( 25) is a perfect model of the process that it represents, or what is the same, that for each k (k-th instant), the following will be verified: y(k) = y mN (k) and, therefore, −y(k) + y mN (k) = 0.
Taking into account now the previous assumptions regarding the perfection of the model and the compatibility of the steady state input with the control law, we can first simplify the expression of the control law, u aN (k) (Equation ( 41)), and secondly substitute in the resulting expression the values of the different variables and coefficients corresponding to the equilibrium point, obtaining the two following relationships: The control law expressed in Equation ( 42) is the one that corresponds to the k-th instant and is a function of the generic state, z mN (k), but also depends on the following variables and parameters: the coincidence horizon, H (which will have been previously chosen); the reference trajectory at (k + H), y r (k + H); the disturbances, d(k); and the system matrices, which depend on the antecedent vector, x a (k).We will assume, as in the study of open-loop stability, that the disturbances will be the same or similar to those corresponding to the steady state (d(k) − d ss ∼ = 0 or d(k) ∼ = d ss ; see Equation (31) and that the antecedent vector, x a (k), will be close enough (in its corresponding space) to x ass , such that the system matrices (and also other matrices that appear in Equation (42) and that are a function of them) are approximately equal to those of the steady state (see Equation (32)).Furthermore, we will assume that what we called in Section 3 as proximity condition of state z mN (k) with respect to the steady state z mNss , will be fulfilled (with its corresponding implications).Finally, we will assume that both the curvature of the reference trajectories, y r (k) (function of the parameters a r i and b r i , i = 1, 2), and the value of H, will be adequate (H, large enough) to be able to consider that y r (k + H) is sufficiently close to the output reference value, which is, as we said, the process output corresponding to the steady state, y mNss ; that is, we will suppose that it is verified: y r (k + H) ∼ = y mNss .Considering all these assumptions, if we subtract expressions (42) and ( 43), we will have the following expression for the increase in the control action (with respect to that of the steady state) or incremental control law: ∆z mN (k) = (z mN (k) − z mNss ) and : and using the change of notation introduced in Equation (37) (i.e., ∆z mN (k) ≡ x inc (k) and ∆u a (k) ≡ u inc (k)), it will be as follows: Remark: It should be noted that, formally, the deduced incremental control law coincides with the known state-feedback control law (u = −Kx), although in our case it was obtained from a previously identified fuzzy model and under the PFC approach.
Once the expression of the incremental control law has been obtained (Equation ( 45)), the next step will be to substitute it in the state equation of the open-loop incremental DLTI model (Equation (38)) and adequately group terms, obtaining the following: Now we will group in a single matrix the matrix expression that multiplies to x inc (k) in Equation ( 46), which we will denote with A mNCL and we will name as a closed-loop state matrix of the incremental model: and, by analogy, we will also define the following matrix: being the relationship between the matrices defined in Equations ( 47) and (48): A mNCL ∼ = A mNCLss , according to the proximity specifications between matrices summarized in Equation ( 46)).Substituting Equation (47) in Equation ( 46), we obtain Equation (49) (below) and repeating the output equation of the open-loop incremental DLTI model (Equation (39), in which u inc (k) does not intervene), we have Equation (50) (below).Ex- pressing both together (Equations ( 49) and ( 50)), we finally have the following closed-loop incremental DLTI model: which constitutes a state-space local model, valid for states close to zero state, of DLTI type (considering the matrices A mNCL and C mN as constants, approximately equal to the corresponding matrices associated with the steady state, z mNss , i.e.,:A mNCL ∼ = A mNCLss and C mN ∼ = C mNss ).Closed-loop local stability: The form of the obtained local DLTI model, given by Equations ( 49) and ( 50) and which we will denote by S CL below, will allow us to apply the same stability criterion in the sense of Lyapunov (internal stability) that was applied in open-loop analysis.Such criterion will now be expressed as follows: S CL is asymptotically stable the eigenvalues o f A mNCL are all strictly within the unit circle (51) that is, the asymptotic stability of S CL will depend on the eigenvalues of A mNCL = A mN − B mN M aN −1 C mN A mN H .However, the analysis of the eigenvalues of this matrix expression is not as simple as in the case of the open-loop, both because of the presence of multiple operations between matrices and the dependence on H, and not in a simple way (H is the exponent of a power of matrix base).To try to solve this problem, a computational approach is proposed in the next section.

Practical Determination of Closed-Loop Local Stability: A Computational Approach
The purpose of this section is to analyze by computation the closed-loop local stability of our case study, carrying out the study through a direct numerical analysis, which can be considered as an alternative approach to the classical analytical analysis (which is not applicable in a direct way in the case of multivariable fuzzy systems).
The stability analysis will be carried out with the help of the so-called symbolic calculation considering generic expressions for the matrices involved in the previously deduced mathematical models of our system.We will use the software tool Symbolic Math Toolbox ™ belonging to the Matlab ® programming environment (The MathWorks Inc., Natick, MA, USA).
The procedure consists of determining the position in the plane of the eigenvalues of the closed-loop state matrix, A mNCL Equation (47), which depends on the coincidence horizon, H, for a succession of increasing values of this parameter, and afterwards the tendency when H tends to infinity is also determined by means of an inductive process.The position of these eigenvalues will be compared with that of the eigenvalues of the open-loop state matrix, A mN (see Equations ( 22) and ( 23)), and conclusions will be drawn about the closed-loop stability, taking into account the stability criterion expressed in Equation (51).
Closed-loop stability analysis by means of symbolic computation: We will start from the following generic specifications for the matrices involved in the considered mathematical models, adding the symbol suffix to the subscript of the matrices names (consistent with the generic name for the system matrices introduced in Equation ( 23)): Of the above matrices, the most significant initially is the open-loop state matrix, A mN_symb , whose eigenvalues (with Matlab code: eig A mN_symb ) are the following: and therefore, the open-loop system will be asymptotically stable ⇔ |a| < 1 ∧ | f | < 1, as already established in Equation ( 40).
Next, we will assign increasing values to H and determine, for each value of H, the closed-loop state matrix and its eigenvalues, by means of symbolic calculation.The use of this type of calculation will allow us to carry out the necessary matrix operations (see the original expression of A mNCL , in Equation ( 47)) with matrices expressed in generic form, including sums, products, and powers with exponent H (H ∈ Z + , H ≥ 1), even for very large H values, and later determine the eigenvalues of the matrix obtained.
The results obtained using the above cited software tool are (for some selected H values) the following (where A mNCL_symb_Hn is the, generically expressed, closed-loop state matrix corresponding to case H = n and f is the second of the eigenvalues shown in Equation ( 53), corresponding to the open-loop state matrix, A mN_symb ): From observation of the results obtained and reasoning by induction, we can deduce, with respect to the eigenvalues of A mNCL_symb_Hn (for all H > 1, H increasing), that the first is always 0 and that the second fits a sequence whose general term depends on H analytically, according to the following expression: That is, the eigenvalues of A mNCL_symb_Hn are the following: where we can observe, for the eig 2 , that the summations in the numerator and denominator satisfy (both) the shape of a geometric series (GS m ) of the following type: k being a nonzero constant coefficient, r the ratio of the series and m the upper level of the summation (in our case: k = 1, r = f and m = H − 2 or m = H − 1).The geometric series of the type shown in Equation ( 60), with k ∈ R, k = 0 and ratio r ∈ R, are convergent (when m tends to infinity) if and only |r| < 1 and, if this is verified, its sum will be: k/(1 − r).That is: and therefore, in our case: if | f | < 1, both the numerator and the denominator of the expression corresponding to the eig 2 eigenvalue, in Equation (59), will tend to 1/(1 − f ) when H tends to infinity.That is, we will have the following: Thus, for the limit case corresponding to H tending to infinity, we will finally have: Next, in Table 7, we will summarize the results obtained in the calculation of the eigenvalues corresponding to the closed-loop state matrix, together with the eigenvalues of the open-loop state matrix: The main characteristic of the trend that we can observe in Table 7 is that, for H → ∞ , the eigenvalues of the closed-loop state matrix, A mNCL , are: 0 and a second eigenvalue, f , which coincides with the second eigenvalue of the open-loop state matrix, A mN .Consequently, we can state the following conclusions: Conclusion 1.When the coincidence horizon tends towards infinity, H → ∞ , one of the two eigenvalues of the closed-loop state matrix is zero and the other coincides with one of the eigenvalues of the open-loop state matrix (the second eigenvalue of ( 53)) if the absolute value of this eigenvalue is less than 1.

Conclusion 2.
If the open-loop plant is locally asymptotically stable (around some equilibrium point), then for a large coincidence horizon ( H → ∞ ), the corresponding closed-loop plant will also be locally asymptotically stable (around from that equilibrium point).

Proof of the Conclusion 2.
If the open-loop plant is locally asymptotically stable, then all the eigenvalues of the open-loop state matrix will be strictly within the unit circle Equation (40).On the other hand, if H tends towards infinity, the eigenvalues of the closed-loop state matrix will be (Conclusion 1): the zero (which is clearly within the unit circle) and a second eigenvalue, which will match some eigenvalue of the open-loop state matrix with absolute value less than 1, if it exists, and this happens, since its eigenvalues are all within the unit circle, as we just established in the starting hypothesis.Therefore, all the eigenvalues of the closed-loop state matrix will also be within the unit circle and, consequently, taking into account the criterion expressed in Equation ( 51), we can affirm that the closed-loop plant will also be locally asymptotically stable.

Stability Test and FMBPC Experiments
This section is dedicated to the presentation of some practical and experimental tests carried out, showing both their description and their results.First, we carry out a stability test, consisting of the verification of the inductive process developed in the previous section, for a specific numerical example, that is, considering numerical values for the elements of the system matrices.Second, as an experimental test of the behavior of the closed loop control system, we analyze the evolution of the controlled variables in two of the multiple FMBPC simulation experiments carried out in the framework of the research corresponding to the present work.

Stability Test
To carry out this test, a specific case was chosen, determined by simulation in the vicinity of an equilibrium point, for a certain antecedent vector, and characterized by certain open-loop system matrices, with specific numerical values.From these matrices, we will determine, for each value of the coincidence horizon, H, the corresponding closed-loop state matrix, A mNCL (according to expression shown in Equation ( 47)), expressed with specific numerical values as well.The numerical open-loop system matrices corresponding to our example are the following: A mN = 0.8000 0.7700 0 0.8800 ; B mN = 0.9000 0.5700 0.8100 0 ; D mN = 0.2900 0.4800 0.2900 0 ; R mN = 0.7900 0.8800 ; C mN = 1 0 0 1 (64) where the most significant matrix for our study is the open-loop state matrix, A mN , whose eigenvalues (eig A mN in Matlab code) are: 0.8000 0.8800 (65) and therefore, since the two eigenvalues are strictly within the unit circle, in this case, the open-loop system would be asymptotically stable (see (40)).
Next, we will repeat the systematic calculation process performed in Section 5 (now with numbers), giving increasing values to H and determining, for each value of H (H = n), the corresponding closed-loop state matrix, A mNCL_Hn Equation ( 47), as well as its eigenvalues.The results obtained in the Matlab ® environment are the following (where we have approximated the numbers using only four decimal places):  From observation of previous calculations, summarized in Table 8, we can conclude that as H increases and tends towards large values, one of the eigenvalues of the closedloop state matrix, A mNCL , is 0.0000, in all cases, and the other increases slowly, always with a value less than 1 (and greater than 0), until it reaches the value 0.8800, which coincides with the second eigenvalue (the greater) of the open-loop state matrix, A mN .For our particular numerical example, Conclusion 1 established in Section 5 has been verified.On the other hand, in our case, we started from a plant that is locally asymptotically stable in open-loop and it has been found that, for sufficiently large values of H, the eigenvalues of the corresponding closed-loop state matrix are all strictly within the unit circle and therefore, according to criterion expressed in Equation ( 51), the corresponding closed-loop plant will also be locally asymptotically stable.That is, Conclusion 2 established in Section 5 has also been verified for our particular numerical example.

FMBPC Control Experiments
Within the framework of this work, many control experiments were carried out applying the considered FMBPC control-law to a simulated ASP process (represented by a continuous model given by differential equations).The objective of these experiments was to test the behavior of our FMBPC control strategy, both for constant reference values and for variable reference values, in the neighborhood of certain operating points in terms of the control performance and the stability of the closed loop system.The simulations were performed using the MATLAB & Simulink software environment (Matlab ® : programming environment of The MathWorks, Inc.).The controlled variables are the substrate concentration in the effluent and the biomass concentration in the reactor (s and x, respectively), the only manipulated variable being the sludge recirculation flow-rate (q r ).
The general control objective is to drive the system to an appropriate operating point (s re f , x re f ) despite the strong variations in disturbance signals (the input flow-rate of contaminated water, q i , and the organic substrate concentration in it, s i ).More precisely, the substrate concentration in the effluent must be kept at a certain reference value fulfilling legal regulations, while the biomass concentration approaches some convenient reference values according to biochemical criteria to guarantee a proper purification ability of the plant, avoiding certain inappropriate dynamic behaviors of the purification process (dependent on the work area) among other operational aspects.
In this section, a set of control experiments considering different changes in the biomass concentration reference while keeping the substrate concentration reference at a constant value, are presented.The action of strong disturbances is also considered.The experiments with the FMBPC strategy are structured in three cases: Case 1, Case 2, and Case 3. The performance of the FMBPC control algorithm will be qualitatively studied in terms of disturbances rejection and references tracking from the closed-loop simulation model.The stability of the closed-loop system, with the FMBPC control strategy, can be qualitatively analyzed if the design of the experiment is appropriate and, both the openloop stability conditions of the plant and the conditions of the controller parameter H, are fulfilled.Consequently, tests to evaluate the performance and stability of the closed-loop control system were developed within an operating range in which the plant is locally open-loop stable and high-enough values of H are chosen.
Classical-PID control Tests: In the reference industrial plant [103], the control algorithm used to control the substrate concentration was the classical-PID algorithm.For this reason, it has been considered that it could be useful to simultaneously carry out tests with a monovariable classical PID control algorithm, for some of the FMBPC cases considered.Specifically, for Cases 1 and 2, two tests with the classical PID control algorithm were carried out: one, focusing the PID to substrate concentration control and the other, focusing the PID to biomass concentration control.The results of these PID tests, as well as a small qualitative study of the stability and performance of this control strategy applied to the control of the ASP process, are also included in this section.
Configuration of the FMBPC experiments and classical-PID control tests: For implementation of the FMBPC strategy, the identified fuzzy models FM 1 and FM 2 , described in Section 2.2 and whose identification parameters, as well as their meaning, are shown in Table 2, were used.Three cases of multivariable FMBPC control are presented, with the substrate and biomass being the variables simultaneously controlled.In the first two cases, the same control experiment is also carried out with the classic PID methodology, in duplicate, considering that the controlled variable is either the substrate or the biomass.The information about the configuration, several characteristics, and parameters considered for these FMBPC and PID control experiments is shown in Table 9.In this table, we will refer to the three FMBPC control experiments as Case 1, Case 2, and Case 3 and to the PID control experiments as PID tests (1a and 1b, corresponding to Case 1, and 2a and 2b, corresponding to Case 2).This table contains (from left to right): an identifier (FM 1 or FM 2 ) for the original fuzzy model of the ASP process, chosen from several models, previously identified; the FMBPC control algorithm design parameters, a r j j=1, 2 and H, which are the outputs reference trajectories parameters and the coincidence horizon, respectively; the information regarding which are the controlled variables in the case of the FMBPC strategy (the substrate, s, and the biomass, x); an identifier (D A or D B ) for the chosen input disturbances, which are certain sets of typical values, or logical variations thereof, of q i and s i variables that were measured in the real WWTP of reference; the simulation interval (0 to 166 h, in all cases); the reference values, s re f and x re f , of the controlled variables, s and x, respectively; the time instants scheduled for the changes in the biomass reference (step time); and, finally, the choice of the variable that will be the control target, that is, the variable controlled, for the implemented classical monovariable PID controller, considering two options: either the substrate (s) or the biomass (x), depending on the PID test.The reference value s re f is kept constant over time with the same value for all the tests.However, the reference value x re f is changed throughout the entire simulation interval, going from a constant value to another in two different time instants for all the tests (two sequences were planned: one for the two sub-cases of Case 1 and another for the two subcases of Case 2 and for the Case 3).For the PID controller, the reference for the controlled variable (s or x) is the same (s re f or x re f , respectively) that considered for the FMBPC in the corresponding case.
The PID controller tuning was carried out in two phases: a first approximation was obtained by applying the Ziegler-Nichols method and later it was adjusted more finely by a trial-and-error procedure to improve the performance of the control system.The PID tuning parameters used were: k p = 1 and k i = 0.1 (in practice, a PI controller) and a bias value for the control signal equal to 775 m 3 /h.Results of the simulation: The results of the different experiments carried out are shown graphically.For all cases corresponding to FMBPC strategy (Cases 1, 2 and 3), the time evolution of each of the two controlled variables (substrate and biomass) is represented, separately, together with their corresponding references and reference trajectories.The two input disturbances are also included (in the graphs of the two controlled variables) to show the magnitude of the values that the controller must compensate.The time evolution of the control variable (sludge recirculation flow-rate) is also represented in a separate graphic.The objective of all the representations is to show the evolution of the controlled variables in terms of disturbance rejection and reference tracking.
Furthermore, for cases 1 and 2 and for each controlled variable, the temporal evolution corresponding to the FMBPC strategy and that corresponding to the classical monovariable PID control algorithm (in its two modalities or tests: substrate control and biomass control) have been represented jointly, although avoiding showing the input disturbances in these graphs, so as not to impair the observation of the controlled variables.The time evolutions of the control variables corresponding to both strategies have also been shown jointly, but in different graphs than those of the controlled variables.
The graphical results obtained are shown below, organized according to the different cases and tests specified in Table 9: 6.2.1.FMBPC-Case 1 and PID-Tests (1a and 1b) Case 1 is characterized by the use of the fuzzy model FM 1 , the disturbance pattern D A , and a certain sequence of changes in the biomass concentration reference signal, within an operation range between 1800 y 2200 (mg/L).The remainder parameters of this experiment are shown in Table 9.The output responses of the system controlled by the FMBPC strategy, as well as corresponding to the two tests performed with the PID (PID-test 1a and PID-test 1b), are presented below in Figures 5-7, respectively: 6.2.1.FMBPC-Case 1 and PID-Tests (1a and 1b).
Case 1 is characterized by the use of the fuzzy model  , the disturbance pattern  , and a certain sequence of changes in the biomass concentration reference signal, within an operation range between 1800 y 2200 (mg L ⁄ ).The remainder parameters of this experiment are shown in Table 9.The output responses of the system controlled by the FMBPC strategy, as well as corresponding to the two tests performed with the PID (PID-test 1a and PID-test 1b), are presented below in Figures 5-7, respectively:  In Figure 5 shown above, the behavior of the ASP process controlled by using the FMBPC strategy is represented, for Case 1. From the observation of the graphic representations presented in this figure, we can conclude that the values of the substrate concentration along the time remain relatively close to their set points and, at the same time, the values of the biomass concentration follow quite closely the jump changes in the corresponding reference signal, following very precisely the predetermined reference trajectory.Furthermore, it can be noted that the control variable varies adequately to reject the disturbances, with moderate control efforts, except in certain cases, in which the two disturbances change simultaneously and abruptly.Moreover, stable behavior of the closed-loop system is observed, despite the strong changes in the biomass reference.In Figure 5 shown above, the behavior of the ASP process controlled by using the FMBPC strategy is represented, for Case 1. From the observation of the graphic representations presented in this figure, we can conclude that the values of the substrate concentration along the time remain relatively close to their set points and, at the same time, the values of the biomass concentration follow quite closely the jump changes in the corresponding reference signal, following very precisely the predetermined reference trajectory.Furthermore, it can be noted that the control variable varies adequately to reject the disturbances, with moderate control efforts, except in certain cases, in which the two disturbances change simultaneously and abruptly.Moreover, stable behavior of the closed-loop system is observed, despite the strong changes in the biomass reference.(c) Control variables calculated using both the FMBPC algorithm and the PID algorithm.Figure 6 shows the time evolution of the two outputs of the ASP process, when it is controlled with two strategies: the FMBPC algorithm (Case 1) and the classic PID control designed, in this particular case, to control the substrate concentration.We can observe that the values of the substrate concentration produced when the PID algorithm is implemented are a little higher than the values of the substrate concentration produced by the FMBPC strategy, but in both cases, the dynamic response is similar.Note that, the tracking of the biomass reference is not guaranteed with this particular PID control while the FMBPC strategy (Figure 6b) allows us to do so.Specifically, it is observed that, during the last time-interval, the biomass values were significantly lower than the corresponding reference signal, what is not recommendable from an operational point of view in this type of plant.This can be avoided when implementing FMBPC because this strategy uses a model that simultaneously considers the dynamic behavior of both the substrate and the biomass with respect to the manipulated variable, as well as the interaction between both variables, and consequently is capable of adequately controlling the substrate, while maintaining the biomass at values prescribed by your reference signal (that is, it allows both variables to be controlled simultaneously).(c) Control variables calculated using both the FMBPC algorithm and the PID algorithm.Figure 6 shows the time evolution of the two outputs of the ASP process, when it is controlled with two strategies: the FMBPC algorithm (Case 1) and the classic PID control designed, in this particular case, to control the substrate concentration.We can observe that the values of the substrate concentration produced when the PID algorithm is implemented are a little higher than the values of the substrate concentration produced by the FMBPC strategy, but in both cases, the dynamic response is similar.Note that, the tracking of the biomass reference is not guaranteed with this particular PID control while the FMBPC strategy (Figure 6b) allows us to do so.Specifically, it is observed that, during the last time-interval, the biomass values were significantly lower than the corresponding reference signal, what is not recommendable from an operational point of view in this type of plant.This can be avoided when implementing FMBPC because this strategy uses a model that simultaneously considers the dynamic behavior of both the substrate and the biomass with respect to the manipulated variable, as well as the interaction between both variables, and consequently is capable of adequately controlling the substrate, while maintaining the biomass at values prescribed by your reference signal (that is, it allows both variables to be controlled simultaneously).In Figure 7, the responses of the ASP process controlled using two strategies are shown together in the same graphic representation.Particularly, the FMBPC algorithm (Case 1) is considered again, and a classic PID oriented, in this case, to the biomass concentration control (the PID changes its target variable from control).In this test, we can observe better PID behavior compared to the previous test (Figure 6).The values of the substrate concentration with the PID are very similar to those of the FMBPC, both with regard to the maximum values and the evolution over time.Regarding the values of the biomass concentration, the PID achieves much more acceptable levels over time than in In Figure 5 shown above, the behavior of the ASP process controlled by using the FMBPC strategy is represented, for Case 1. From the observation of the graphic representations presented in this figure, we can conclude that the values of the substrate concentration along the time remain relatively close to their set points and, at the same time, the values of the biomass concentration follow quite closely the jump changes in the corresponding reference signal, following very precisely the predetermined reference trajectory.Furthermore, it can be noted that the control variable varies adequately to reject the disturbances, with moderate control efforts, except in certain cases, in which the two disturbances change simultaneously and abruptly.Moreover, stable behavior of the closed-loop system is observed, despite the strong changes in the biomass reference.
Figure 6 shows the time evolution of the two outputs of the ASP process, when it is controlled with two strategies: the FMBPC algorithm (Case 1) and the classic PID control designed, in this particular case, to control the substrate concentration.We can observe that the values of the substrate concentration produced when the PID algorithm is implemented are a little higher than the values of the substrate concentration produced by the FMBPC strategy, but in both cases, the dynamic response is similar.Note that, the tracking of the biomass reference is not guaranteed with this particular PID control while the FMBPC strategy (Figure 6b) allows us to do so.Specifically, it is observed that, during the last time-interval, the biomass values were significantly lower than the corresponding reference signal, what is not recommendable from an operational point of view in this type of plant.This can be avoided when implementing FMBPC because this strategy uses a model that simultaneously considers the dynamic behavior of both the substrate and the biomass with respect to the manipulated variable, as well as the interaction between both variables, and consequently is capable of adequately controlling the substrate, while maintaining the biomass at values prescribed by your reference signal (that is, it allows both variables to be controlled simultaneously).
In Figure 7, the responses of the ASP process controlled using two strategies are shown together in the same graphic representation.Particularly, the FMBPC algorithm (Case 1) is considered again, and a classic PID oriented, in this case, to the biomass concentration control (the PID changes its target variable from control).In this test, we can observe better PID behavior compared to the previous test (Figure 6).The values of the substrate concentration with the PID are very similar to those of the FMBPC, both with regard to the maximum values and the evolution over time.Regarding the values of the biomass concentration, the PID achieves much more acceptable levels over time than in the previous test, following quite well the selected changes in the reference signal, although the tracking capability is less than that obtained with the FMBPC strategy.By comparing the control actions of both strategies, we see again that the PID control actions are not able to make the necessary efforts to achieve the biomass concentration reference values properly in the face of high disturbances action (while maintaining the substrate concentration values around the reference values, acceptably).In the case of FMBPC, the biomass concentration accurately follows the reference trajectory.The multivariable control objective is successfully achieved with our strategy, even using a single manipulated variable, the sludge recirculation flow-rate.
Note, in the case of the PID, that a better tuning (more aggressive) would allow better tracking of the biomass reference but causing greater variability in the values of the substrate concentration, with the risk of even exceeding values not allowed by environmental regulations.This is to be expected because the PID does not incorporate information on the interaction between substrate and biomass.Consequently, the tuning implemented for the PID was the result of the search for an equilibrium in the variations of both variables, placing the accent on the fulfillment of the substrate restrictions.

FMBPC-Case 2 and PID-Tests (2a and 2b)
Case 2 is characterized by the use of the fuzzy model FM 2 , the disturbance pattern D B , and a particular sequence of changes in the reference of the biomass concentration, within an operation range between 500 × 1300 (mg/L).The remainder parameters of this experiment are shown in Table 9.The responses of the system controlled by the FMBPC strategy, as well as the corresponding to the two tests performed with the PID (PID-test 2a and PID-test 2b), are presented below in Figures 8-10, respectively:  In Figure 8, shown above, the behavior of the ASP process controlled by using the FMBPC strategy is represented, for Case 2. From the observation of the time evolution of the two outputs represented in this figure, we can conclude that, as well as for Case 1, the values of the substrate concentration remain relatively close to its reference signal despite the disturbance action and, at the same time, the biomass concentration values follow quite fast the abrupt jumps in the corresponding reference variable, following very accurately the predefined reference trajectory.Moreover, it can be observed that the changes in the control variable and the control efforts to reject disturbances are softer than for Case 1 all the time.Note that, although the input flow variations are fairly large and, besides, are subjected to strong oscillations, the control system presents a very good disturbance rejection capability.This fact can also be observed by comparing the profile of the control signal time evolution and the input flow profile, particularly at times when the changes in the input flow are the largest ones.It must be pointed out that the control actions respond very quickly after each oscillation.It can also be observed that the closed-loop system keeps a stable behavior in the presence of strong changes in the biomass reference.A conclusion from this second case study, with important differences with respect to the first (in the fuzzy model, in the disturbances and in the operating zone), can be drawn: the FMBPC strategy also allows to obtain very good performance, in terms of disturbance rejection and reference tracking, while ensuring stable behavior of the closed-loop system.In Figure 9, the responses of the ASP process provided by the two control strategies are shown together: the FMBPC (Case 2) and the classic PID oriented to control the substrate concentration, respectively.In this case, we observe that the values of the substrate concentration provided by the PID control algorithm are lower than the values of the substrate concentration corresponding to the FMBPC strategy.Now, if we look at the graphic   In Figure 10, the responses of the ASP process controlled using the two strategies are shown in the same graphic.Particularly, the FMBPC algorithm corresponding to Case 2 is implemented again but, this time, a classic PID oriented to biomass concentration control is considered (the PID changes its target variable from control).Regarding the time evolution of biomass concentration (Figure 10b), we can observe a better PID behavior compared to the previous test (Figure 9).Since the biomass is now the controlled variable for the PID control, we can observe that the closed-loop system has quite a good capability to follow the changes in the biomass reference.On the other hand, we can observe a similar reduction in substrate concentration for both strategies (Figure 10a), although a little more favorable to FMBPC.Finally, it is also interesting to observe the control actions of both strategies (Figure 10c).Roughly speaking, they are quite similar; nevertheless, some of the differences are interesting.Firstly, the PID controller presents a greater number of situations in which the control efforts are larger, which could be bad for the actuator.Secondly, the operation under the PID control presents saturation in the control signal, with values equal to zero during a period of several hours (in the time interval between 20 h and 40 h of simulation time), a circumstance that does not occur in the case of the FMBPC strategy.

FMBPC-Case 3
Case 3 is basically characterized by the use of the fuzzy model  , the disturbance pattern  , and a particular sequence of changes in the reference of the biomass concen- In Figure 8, shown above, the behavior of the ASP process controlled by using the FMBPC strategy is represented, for Case 2. From the observation of the time evolution of the two outputs represented in this figure, we can conclude that, as well as for Case 1, the values of the substrate concentration remain relatively close to its reference signal despite the disturbance action and, at the same time, the biomass concentration values follow quite fast the abrupt jumps in the corresponding reference variable, following very accurately the predefined reference trajectory.Moreover, it can be observed that the changes in the control variable and the control efforts to reject disturbances are softer than for Case 1 all the time.Note that, although the input flow variations are fairly large and, besides, are subjected to strong oscillations, the control system presents a very good disturbance rejection capability.This fact can also be observed by comparing the profile of the control signal time evolution and the input flow profile, particularly at times when the changes in the input flow are the largest ones.It must be pointed out that the control actions respond very quickly after each oscillation.It can also be observed that the closed-loop system keeps a stable behavior in the presence of strong changes in the biomass reference.A conclusion from this second case study, with important differences with respect to the first (in the fuzzy model, in the disturbances and in the operating zone), can be drawn: the FMBPC strategy also allows to obtain very good performance, in terms of disturbance rejection and reference tracking, while ensuring stable behavior of the closed-loop system.
In Figure 9, the responses of the ASP process provided by the two control strategies are shown together: the FMBPC (Case 2) and the classic PID oriented to control the substrate concentration, respectively.In this case, we observe that the values of the substrate concentration provided by the PID control algorithm are lower than the values of the substrate concentration corresponding to the FMBPC strategy.Now, if we look at the graphic showing the time evolution of the biomass concentration (Figure 9b), a tracking capability of its reference trajectory fairly optimal is observed when the FMBPC strategy is used, while with the PID algorithm, the biomass reference tracking cannot be ensured, as shown in the figure, and was expected.It must be noted at this point that the PID is a monovariable control law and, particularly for this test, it was designed with the aim of controlling just the substrate concentration.Nevertheless, it must be highlighted that the FMBPC strategy accomplishes the multivariable control objectives very well, keeping the two controlled variables, the substrate and the biomass concentrations, close to their reference values using just one manipulated variable.The reason why in this test the substrate values slightly higher using the FMBPC than using PID control is clear and a logical consequence of having forced the biomass to follow too low values, so that the substrate cannot be reduced as much as it was expected.This fact can be sort out by setting higher values of the biomass reference as in Case 1.An important conclusion is drawn from this test: a compromised selection of the two reference variables is needed to keep the biomass at high-enough levels to ensure a good purification process of the water, while keeping the substrate at low-enough values according to its reference signal and legal specifications.
In Figure 10, the responses of the ASP process controlled using the two strategies are shown in the same graphic.Particularly, the FMBPC algorithm corresponding to Case 2 is implemented again but, this time, a classic PID oriented to biomass concentration control is considered (the PID changes its target variable from control).Regarding the time evolution of biomass concentration (Figure 10b), we can observe a better PID behavior compared to the previous test (Figure 9).Since the biomass is now the controlled variable for the PID control, we can observe that the closed-loop system has quite a good capability to follow the changes in the biomass reference.On the other hand, we can observe a similar reduction in substrate concentration for both strategies (Figure 10a), although a little more favorable to FMBPC.Finally, it is also interesting to observe the control actions of both strategies (Figure 10c).Roughly speaking, they are quite similar; nevertheless, some of the differences are interesting.Firstly, the PID controller presents a greater number of situations in which the control efforts are larger, which could be bad for the actuator.Secondly, the operation under the PID control presents saturation in the control signal, with values equal to zero during a period of several hours (in the time interval between 20 h and 40 h of simulation time), a circumstance that does not occur in the case of the FMBPC strategy.

FMBPC-Case 3
Case 3 is basically characterized by the use of the fuzzy model FM 2 , the disturbance pattern D A , and a particular sequence of changes in the reference of the biomass concentration, within an operation zone between 500 y 1300 (mg/L).The remainder parameters can be seen in Table 9.The configuration of this case is the same as Case 2, except for the disturbance pattern, which is coincident with that of Case 1. Precisely, the objective of the experiments in this case is to study the robustness of the FMBPC strategy by using a pattern of disturbances and operating conditions different from the ones used for the fuzzy model identification considered in the experiments.It must be taken into account that fuzzy models are identified with certain patterns of disturbances and in certain areas of operation and therefore, although the ranges of the values for different variables in the identification are wide, the validity of the models will be higher for such specific areas and will decrease when we move to remote ones.The responses of the system controlled by the FMBPC strategy are shown below, in Figure 11: pattern of disturbances and operating conditions different from the ones used for the fuzzy model identification considered in the experiments.It must be taken into account that fuzzy models are identified with certain patterns of disturbances and in certain areas of operation and therefore, although the ranges of the values for different variables in the identification are wide, the validity of the models will be higher for such specific areas and will decrease when we move to remote ones.The responses of the system controlled by the FMBPC strategy are shown below, in Figure 11:  Figure 11, shown above, shows the behavior of the ASP process controlled with the FMBPC strategy for Case 3. From observation of the figure, we can conclude that, as for Case 1 and Case 2, the substrate concentration values remain relatively close to the reference and, at the same time, the biomass concentration values also change accordingly the jumps in reference signal presents, following, especially in this case, an accurately predetermined reference trajectory.From a more detailed analysis of Figure 11, it can be noted that the substrate levels in Case 3 are higher than the corresponding ones in Case 2, due to the fact that the variations of the substrate in the inflow water for Case 3 are much greater than the corresponding ones for Case 2. Furthermore, although the control actions should compensate for the excess substrate by providing an adequate sludge recirculation Figure 11, shown above, shows the behavior of the ASP process controlled with the FMBPC strategy for Case 3. From observation of the figure, we can conclude that, as for Case 1 and Case 2, the substrate concentration values remain relatively close to the reference and, at the same time, the biomass concentration values also change accordingly the jumps in reference signal presents, following, especially in this case, an accurately predetermined reference trajectory.From a more detailed analysis of Figure 11, it can be noted that the substrate levels in Case 3 are higher than the corresponding ones in Case 2, due to the fact that the variations of the substrate in the inflow water for Case 3 are much greater than the corresponding ones for Case 2. Furthermore, although the control actions should compensate for the excess substrate by providing an adequate sludge recirculation flow-rate, as previously reasoned, this action is compromised by the requirement to follow certain low reference values for biomass (multivariable control).Regarding the control variable, it is observed that it acts adequately to reject the disturbances, as in Case 1 and Case 2, with moderate control efforts in general, although it requires greater efforts in some specific moments (initially and in the reference changes of biomass).Finally, we can observe that the sludge recirculation flows used to control the two variables are different (lower) in Case 3 than in Case 2. This is due to the differences between the patterns of disturbances in both cases.
Overall analysis of the results: From the observation of Figures 5, 8 and 11, we can conclude that our FMBPC closed-loop multivariable control system is stable in the face of reference changes, presenting a satisfactory reference tracking capability.In all cases studied, the FMBPC-strategy reacts adequately to steps in the biomass reference, in such a way that this variable follows its reference trajectory fairly closely and, at the same time, the substrate remains around its reference value, with acceptable deviations despite the strong disturbances present at the WWTP input.This behavior occurs for different patterns of disturbances (all with many oscillations), different steps sequences of the biomass reference and different operating points of the ASP process (the variation ranges of the references of output variables were: s ref = 55 (mg/L), for the substrate, and x ref ∈ [500, 1300] (mg/L) or x ref ∈ [1800, 2200] (mg/L), for the biomass).
The other graphic representations (Figures 6, 7, 9 and 10) show the evolution of the different variables both for the FMBPC strategy and for the PID, simultaneously.Looking at these figures, we can deduce that both strategies are stable within the considered operating range.Regarding the performance of both control algorithms, the performance of the FMBPC algorithm is, in general, better than that of the monovariable PID algorithm, this difference being much greater in the cases in which the PID is oriented to substrate concentration control and not to the biomass concentration control.In these cases, the control of the substrate concentration with both strategies is similar, while the control of the biomass concentration with the FMBPC algorithm is much better than with the PID, as is logical, since the FMBPC strategy is multivariable and is oriented to simultaneously control both substrate and biomass.For the cases in which the PID is oriented to biomass concentration control, the differences are smaller, but the performance of the FMBPC algorithm continues to be better, mainly concerning biomass concentration control (in relation to substrate concentration control, the performances of both algorithms are closer).As mentioned above, the performance of the PID in relation to the tracking of the biomass reference could be improved, but at the cost of increasing variations in substrate concentration, with the consequent risk of not complying with environmental restrictions.
As a summary, we can highlight that the FMBPC multivariable advanced control strategy considered in this work would be competent to address the control of the substrate concentration in the effluent (reduction of contamination) in WWTP with ASP processes, in a stable way and with a performance similar or better to that of a classical PID, and simultaneously, also control the biomass concentration in the reactor, with high performance.
Remark: It is necessary to remember again here that the classical PID control algorithm was the technique implemented in the reference industrial plant [103], whose simulating model together with real data of the disturbance records are used in this work.For this reason, it has also been studied here.Being aware that the comparison with the FMBPC is not fair according to the different degrees of complexity of both approaches, the aim is just to show how the performance of the plant can be improved by using advanced control techniques and to encourage their use in a real environment.
Comparison with other more advanced PID techniques is also possible and recommendable.For instance, PID methods that consider a feedforward action to take into account the effect of measurable disturbances (see [105]) could be appropriate.In addition, comparison could be made with multivariable PID control or fuzzy PID control, and with other controllers using explicit control laws (as in the case of our FMBPC strategy).However, a controller benchmarking is out of the scope of this paper, although it could be very interesting for future research.

Conclusions
In this paper, a computational approach of closed-loop stability analysis of a specific FMBPC control strategy, applied to an activated sludge biological process, was presented.The original model used for the predictions is a TS type fuzzy discrete model, previously identified and later formalized in an equivalent DLTV state-space model.Based on this model but considering a certain generic steady-state as the operating point, in this work we have deduced, for analysis purposes, a local incremental state-space model of the DLTI type, valid for states close enough to steady-state (that is, for small variations or increases with respect to the steady-state).The great advantage of working with models of that type is the possibility of applying the existing stability criteria for DLTI models, which are well defined and widely used, and demonstrating their compliance for the chosen case study.However, the demonstration procedures of compliance with such stability criteria for systems with complex dynamics, as is our case study, are often mathematically quite laborious and also difficult to generalize.For this reason, in the present work, it was decided to carry out stability analysis using a computational approach, alternative to the usual procedures existing in the literature.
The solution adopted consists of determining, first, the generic closed-loop state matrix (with algebraic variables) and its eigenvalues, as a function of the coincidence horizon H, by means of symbolic numerical calculation and an induction process (considering H, increasing).Second, deduce under what conditions the eigenvalues will be within the unit circle and consequently (stability criteria) when the closed-loop plant will be asymptotically stable and its relationship with the stability of the open-loop plant.The result is satisfactory, in the sense that by means of this procedure it has been possible to establish an important conclusion regarding closed-loop (local) stability, which is the following: if the open-loop plant is (locally) asymptotically stable, then for values large enough of H, the closed-loop plant will also be (locally) asymptotically stable.
As a final assessment, we consider that the work carried out constitutes an appreciable contribution in the field of stability analysis of FMBPC control systems, for two reasons.On the one hand, due to the simplicity of the method developed, in relation to other methods, most of them are more complicated.On the other hand, it has been considered a more complex case study and more difficult to approach than those chosen in previous similar works (SISO systems with not very complicated dynamics).In the present work, a useful conclusion regarding the stability of FMBPC control systems has been deduced, for a case study of a multivariable nature (MIMO system), highly nonlinear, with complex dynamics (biological system), and subject to strong disturbances.
As possible future research works, we propose the following: to generalize the proposed FMBPC methodology for any MIMO system and to integrate it within a Closed-loop Predictive Control scheme to manage input and output constraints, considering at the same time the stability analysis of the resulting system.This will allow the use of the proposed methodology for the control of more complex systems and for complete plant control (Plantwide Control).• mr = max.(mr 1 , mr 2 ): common number of rules • a ji , b ji , and r j , the coefficients of the consequent vector, and the independent term, respectively, in the j-th rule of the fuzzy model corresponding to output y 1 (k) • a * ji , b * ji , and r * j , the coefficients of the consequent vector, and the independent term, respectively, in the j-th rule of the fuzzy model corresponding to output y 2 (k) (with the following particularity: a State-space global model of the ASP process, corresponding to [23], with the elements of the vectors and system matrices expressed in a generic way: the two state space models, the following relations can be easily deduced (see all the equations previously shown in this Appendix and especially Equations (A9), (A10), and (A13), (A14)): In a similar way, we could transform the expressions of two other matrices present in Equation (A16), M a (specified in Equation (A15)) and λ m (specified in Equation (A16)), in terms of the matrices of the model of Section 2.3 of this article, making use of the relations included in Equation (A17).However, it is not necessary to detail these transformations and it is preferable to keep each of these expressions (matrix operations between system matrices) as a single compact matrix, adding a suffix N referring to the corresponding transformation.We will denote those adapted matrices as M aN and λ mN , respectively (and M aN −1 for the inverse of M aN ).On the other hand, we will adapt the corresponding notation to the output of the model: y m (k) ≡ y mN (k).Finally, considering all these notation adaptations and the development made in (A18), and replacing all of it in Equation (A16), the expression of u a (k), as a function of the vector and matrix variables corresponding to the model obtained in Section 2.3 of this article, will be the next:

Figure 1 .
Figure 1.The activated sludge process (ASP) and the inputs and outputs taken into account (multivariable nonlinear system).

Figure 1 .
Figure 1.The activated sludge process (ASP) and the inputs and outputs taken into account (multivariable nonlinear system).

Figure 2 .
Figure 2. The multivariable architecture of the ASP process (3 inputs and 2 outputs).

Figure 2 .
Figure 2. The multivariable architecture of the ASP process (3 inputs and 2 outputs).

j being : mr 2 :
number o f rules o f f uzzy model, f or output y 2(12)
(a) Substrate concentration in the effluent and disturbances.(b) Biomass concentration in the reactor and disturbances.Processes 2021, 9, x FOR PEER REVIEW 32 of 50 (c) Control variable calculated using the FMBPC algorithm.
(a) Substrate concentration in the effluent.(b)Biomass concentration in the reactor.
(a) Substrate concentration in the effluent.(b) Biomass concentration in the reactor.
(a) Substrate concentration in the effluent.(b) Biomass concentration in the reactor.
(a) Substrate concentration in the effluent.(b) Biomass concentration in the reactor.Processes 2021, 9, x FOR PEER REVIEW 34 of 50 (c) Control variables calculated using both the FMBPC algorithm and the PID algorithm.
(a) Substrate concentration in the effluent and disturbances.(b) Biomass concentration in the reactor and disturbances.(c)Control variable calculated using the FMBPC algorithm.
(a) Substrate concentration in the effluent.(b) Biomass concentration in the reactor.(c)Control variables calculated using both the FMBPC algorithm and the PID algorithm.
(a) Substrate concentration in the effluent and disturbances.(b) Biomass concentration in the reactor and disturbances.(c) Control variable calculated using the FMBPC algorithm.

Table 1 .
Inputs and outputs of the activated sludge process (ASP) and its role.

Table 1 .
Inputs and outputs of the activated sludge process (ASP) and its role.

Table 2 .
Identification parameters of the FM 1 and FM 2 fuzzy models.

Table 4 .
Coefficients of the consequent as per the rules of the identified FM 1 fuzzy model, and offset terms.

of the Consequent Vector, x, and Corresponding Coefficients ASP Output
Rule y 1 (k−1)

Table 7 .
Eigenvalues of the state matrices obtained by means of symbolic computation.

Table 8 .
Eigenvalues of the state matrices obtained by means of direct computation.

Table 9 .
Configuration of the FMBPC control experiments and classical PID tests.

Table A1 .
Antecedent fuzzy sets and membership functions of the identified FM 1 fuzzy model.

Table A2 .
Antecedent fuzzy sets and membership functions of the identified FM 2 fuzzy model.