Collaborative Control Applied to BSM1 for Wastewater Treatment Plants

: This paper describes a design procedure for a collaborative control structure in Plant Wide Control (PWC), taking into account the existing controllable parameters as a novelty in the procedure. The collaborative control structure includes two layers, supervisory and regulatory, which are determined according to the dynamics hierarchy obtained by means of the Hankel matrix. The supervisory layer is determined by the main dynamics of the process and the regulatory layer comprises the secondary dynamics and controllable parameters. The methodology proposed is applied to a wastewater treatment plant, particularly to the Benchmark Simulation Model No 1 (BSM1) for the activated sludge process, comparing the results with the use of a Model Predictive Controller in the supervisory layer. For determining controllable parameters in the BSM1 control, a new speciﬁc oxygen mass transfer model in the biological reactor has been developed, separating the k L a volumetric mass transfer coefﬁcient into two controllable parameters, k L and a .


Introduction
The plantwide control (PWC) approach is the most widely used solution for optimization and control of large-scale processes. This methodology provides a link between the systems and strategies necessary to control a large scale chemical plant comprised by different interacting units [1]. PWC has been a very active research field due to the significant improvements to the operation of industrial processes [2]. There are different approaches, broadly classified in centralized, decentralized, and distributed control. In Figure 1, a general arrangement for PWC methodologies is presented, and more details are given in [3].
The simplest approach is a centralized control structure, where a unique coordinator controller either selects the set points for other controllers in the lower layer or sends information directly to the final control elements (FCE) [4]. In a hierarchical structure, the coordinator controller is typically an optimization based controller, and PID controllers are in the lower level, which focuses on plant behavior and its dynamic interactions [5,6]. The optimization based controller can be an economic one including costs in the objective function, or a standard Model-based predictive controller (MPC). The lower level typically consists of several independent PID control loops but multiple MPC control loops can also be considered. Model-based predictive control (MPC) is a powerful strategy frequently used in Decentralized and Distributed control structures. The MPC algorithm performs an optimization at each sampling time to obtain the optimal values of the manipulated inputs [7]. This advanced control strategy considers an objective function that includes tracking errors and control efforts properly balanced and constraints on control and input variables in the formulation of the optimization problem. Despite of the large amount of MPC successful applications, especially in the process industries [8], they have some drawbacks, such as the computational load required to solve the optimization problems involved, in particular when the prediction model is nonlinear. Another drawback is the difficulty for tuning the controller parameters, which is often performed by a trial and error procedure. Therefore, in this work, collaborative control is proposed as a practical multivariable control strategy for implementing a hierarchical centralized control structure, and it is compared to the use of MPC as an alternative strategy. In this strategy, the key point is the determination of a dynamics hierarchy using the Hankel matrix that gives the impact of the inputs on the plant states. The main dynamics is located in the upper layer and the rest of control handles can be used as collaborative regulatory controllers aiming to achieve the global objective. Collaborative characteristics imply a direct relation among all secondary dynamics and main dynamics in order to regulate the main dynamic, unlike in common MPC strategies.
Therefore, two levels are considered in the collaborative hierarchical structure: a supervisory layer where the coordinator controller provides set points for the regulatory layer [5], and the regulatory layer comprised by the controllers receiving this information. Usually, local PID controllers are employed in the lower level, each one associated with one process unit. However, it is important to highlight that collaborative control is different from cooperative control. Collaboration occurs when individuals from different species work together for a common goal, while cooperation concept implies that individuals from the same species are working for the same objective. Thus, the conditions for the implementation of collaborative control strategies are commonly found in chemical processes, where each subsystem (e.g., reaction units, separation processes, transportation) pursues a common objective, the final product, even though dynamics of the different subsystems are different (e.g., temperature, concentration, level). On the other hand, classic examples of cooperative control are found in robotics applications, where a team of robots with similar characteristics work for a common goal.
The concept of collaboration is considered in different fields of knowledge. The measurement of collaborative intelligence in knowledge based service planning is studied in [9] who identify measures to achieve the best collaboration between different enterprises. The application of collaborative control in a robotic wheelchair is addressed in [10] introducing a collaborative control mechanism that assists users when they require some help. In [11], a collaborative control strategy is used for teleoperation systems, where the human is modeled as a collaborator in a multidisciplinary team. They say that in "collaborative control model a human operator and a robot are peers who work together, collaborating to perform tasks and to achieve common goals".
Although few applications of collaborative control in chemical processes are found, some works have employed this strategy for hierarchical control with good results. In [12], a multilayer architecture is considered to control a continuous bioethanol process. They use two kinds of PWC schemes and compare single layer with multilayer structures. In [6], a collaborative control structure with communication between the process units and a coordinator controller is presented. In this methodology, the Hankel matrix is used to classify the different dynamics. Similar approaches can be found in [12][13][14]. The works of [15][16][17] explore the possibility of introducing some kind of sensitivity terms in the objective function for the optimization of large-scale process. Although all these works demonstrate the potential of collaborative control strategy, some further improvement could be achieved including controllable parameters in the control structure.
The main contribution of the collaborative control strategy proposed in this work is the possibility of exploiting the manipulation of some process parameters to improve the overall process performance. The capabilities of the collaborative control strategy are improved introducing one or more process controllable parameters into the strategy as additional control handles. Another contribution of this work is the application of this methodology to wastewater treatment plants (WWTP), especially to the activated sludge process. The Benchmark Simulation Model No. 1 (BSM1) is selected to represent the activated sludge process, and this is a recognized benchmark plant designed for the evaluation of different control methodologies. There are many applications of centralized linear and nonlinear MPC strategies to WWTP control, some of them in a hierarchical framework [18][19][20][21][22]. However, to the authors' knowledge, no collaborative control applications to WWTP have been found in the literature.
Finally, it is important to mention that all works about WWTP control use k L a (volumetric oxygen mass transfer coefficient) in the reactor model as a unique manipulated variable, while, in this work, a more realistic model for liquid-gas mass transfer has been developed, adding some controllable parameters to the process. Particularly, the controllable parameters have been computed increasing the degrees of freedom of the process of oxygen transference by separating the volumetric mass transfer coefficient into the liquid phase mass transfer coefficient (k L ) and the specific area for mass transfer (a). This paper is organized as follows: Section 2 presents the collaborative control strategy together with the methodology. In Section 3, the MPC as alternative methodology for the supervisory layer is presented. The case study, the Benchmark Simulation Model No. 1 (BSM1) and the mass transfer model development is described in Section 4. The application of the methodology to BSM1 is described in Section 5, and Sections 6 and 7 show the results and conclusions.

Collaborative Control Strategy
In the following, the methodology to develop a collaborative control proposed is described. In summary, collaborative control is a multivariable control strategy useful for processes with very different dynamic effects. All of the individual control loops work together, in spite of each one having different effects on the process and dynamic behavior. Figure 2 illustrates the proposed methodology.
The common practice in control is to take the process states as controlled variables, considering that all states can be measured. In this work, some process parameters, called controllable parameters, are additional controlled variables, which can be directly manipulated through available control actions. They give extra degrees of freedom to the process control adding more control loops. This option of control structure is particularly interesting where the process is strongly affected by the changes on some parameters, although sometimes this fact is not evident due to some process model simplifications.
The methodology proposal is formed by three stages. The first one is the process analysis, aiming to get a full knowledge of the process. This implies obtaining a mathematical description of the process, a model. The second stage consists of designing the regulatory control layer. This stage distributes PID control loops over the whole plant, following the decentralized control approach. These PID controllers act on tracking mode on secondary dynamics, or over controllable parameters. Only one PID in regulatory mode acts over the the main dynamic of the process. Finally, the third stage designs the supervisory control layer using an optimization problem formulated for including information regarding the main and the secondary dynamic, and the controllable parameters as an additional and new part. Following the procedure of [6], but including some modifications, each stage is explained in the next subsections.

Process Analysis
This obtains a phenomenological based semi-physical model for the treated process. Manipulated variables, states, and parameters are determined and a set of controllable parameters is established in accordance with a process dynamic response. Finally, the dynamic hierarchy is stated for the given process operation point.
Step 1: Obtaining a phenomenological based semi-physical model (PBSM) for the treated process. Developing a PBSM for the process is the first step for designing the proposed collaborative control structure. It should be noted that, if a PBSM is not available, the model can be obtained following the proposal presented in [23,24].
Step 2: Determinate available manipulated variables, states and parameters. For the treated process, the available manipulated variables u and the state variables x are stated from the mathematical model. During model simulation, it is possible to consider any variable or parameter as measured. However, in a real process, the available information depends on installed instruments or designed observers for these variables or parameters.
Step 3: Determinate the controllable parameters. If controllable parameters can be used in the process, this step determines the parameters to be considered as controllable ones. These parameters give the possibility of being controlled and fulfill the characteristics previously mentioned. That set of parameters are used to implement the collaborative control approach proposed here.
Step 4: Define the operating point and span of states and input variables. The real operation of the process is analyzed in this step. The stationary states of the process and its operation conditions are chosen. Such values will be used to determine the main dynamic reference value, the only fixed set point. In this step, the span of process states and input variables is defined in accordance with process operative conditions.
Step 5: Determine the dynamic hierarchy. This is the final step of process analysis, which results in being crucial for the collaborative control structure due to the result obtained defining the design of the supervisory layer. For this, the procedure proposed by [6,25] is used here. This procedure is based on singular value decomposition (SVD) of the Hankel matrix of the process, using the state impactability index (SI I x k ).
The Hankel matrix, H, is the input-output behavior of the process, establishing a relation between a sequence of past inputs and a sequence of future outputs of the process [6]. H is defined by the following Equation: The SVD of Hankel matrix, showed in Equation (2), allows for calculating a measure for the dynamic effect of each process variable over process response and to determine the input-output pairings.
The state impactability index, presented in (3), is the impact of the manipulated inputs u of the process (taken as a whole) over the k-th process output variable y k : where σ 2 ii and U are found using a SVD of the Hankel matrix H of the process model. p is defines as the number of non-negative singular values, n the number of state variables, σ ii the i-th singular value and U and V the orthonormalized eigenvectors HH T and H T H, respectively. In this way, SI I x k is used to determine the main and secondary dynamics. Then, the main dynamic will be the state with higher SI I x k . The rest of states will be the secondary dynamics.
In summary, the Hankel matrix is calculated in order to determine the dynamics hierarchy in the system and to select the dynamics for the supervisory layer. This procedure is performed offline, and no further results are obtained using calculations based on the Hankel matrix, which prevents obtaining unreliable results due to a large matrix condition number. In fact, the condition number is usually large, allowing for the division of dynamics between regulatory and supervisory layers.

Regulatory Layer Design
This task is dedicated to analyzing each one of the individual closed loops to be implemented on the plant. Such control loops include those controllers designed at the commissioning of the process to operate in regulatory mode, following the common practice in decentralized control strategy. Collaborative control strategy changes the secondary dynamic controllers from their regulatory mode to operate in tracking mode. These loops should track the set point given by the coordinator controller. In this level, the current proposal uses PID controllers due to this control strategy being easy to implement and understand.
Step 6: Identification of control objective. It consists of defining the control objective for the process as a whole and establishing the control objective of the pieces of equipment forming the process. The control objective is determined in accordance with hierarchical dynamics defined in Step 5 using Equation (3). In this sense, main dynamic is related directly with the main control objective of the process. It should be noted that it is possible to require additional control loops, which will be determined as it is declared in the following step.
Step 7: Selection of input-output pairing for control loop at each process unit, based on steps 3 to 5. In congruence with the control objectives of the process, this step states the control loops included in the regulatory layer. It is possible to use diverse strategies to determinate an adequate input-output pairings. For example, the RGA approach [26] or the SVD of the Hankel matrix as in [6].
Step 8: Preliminary control tuning. To finish the design of regulatory layer, the tuning of each proposed individual controller is executed here. For this task, any available tuning method for single-input/single-output (SISO) loops can be used, considering that any interaction among control loops should be considered on the supervisory layer. Nevertheless, these controllers should be individually evaluated regarding their performance, previously to design the supervisory layer.

Supervisory Layer Design
This step states an optimization problem considering the main and secondary dynamics, and the selected controllable parameters. The decision variables are the reference points for each one of the regulatory individual controllers. The values of such reference points look for reaching the global control objective through the minimization of main dynamic error, which is the only one variable with a control loop operating in regulatory mode.
Step 9: Formulation of the optimization problem. It is considered that the performance of regulatory level is good using fixed set points, but it is not enough regarding desired objectives of the plant wide control (PWC). The supervisory layer has a direct relation with the main dynamic and the stated optimization, with a direct collaboration of secondary dynamics and controllable parameters in order to fulfill the objectives of the plant as a whole. It should be noted that the effect of every secondary dynamic and controllable parameter over the main dynamic is exploited using the dynamics relations described in the process model. In this way, the set points of the control loops of secondary dynamics and controllable parameters are degrees of freedom in this optimization problem [6]. The optimization problem is formulated as: where x j,sp are the controlled variables and θ k are the controlled parameters. α, β i , and γ k correspond to weights, which operates as tuning factors, x * sp represents the set point of the main dynamic (x * ) of the process, and x * ss is the value of x * at the steady-state that would be achieved if the process reached its steady-state over all the x j,sp values determined by the optimization and the rest of states finished at the current value. S i corresponds to a sensitivity function of main dynamic regarding changes on i-th secondary dynamic at the current values of the process state, Ss k corresponds to the sensitivity function of main dynamic to changes on k-th controllable parameter at the current values of the process state. x i,sp min , x i,sp max , θ k,sp min and θ k,sp max represent limit values stated for the set points of secondary dynamics and controllable parameters, looking for guaranteeing the stability of the process. ∆x j,sp min , ∆x j,sp max , ∆θ k,sp min and ∆θ k,sp max indicate the maximum and minimum possible changes to apply in the references of secondary dynamics and controllable parameters.
The term S i represents the effect of every secondary dynamic over the main dynamic and it is calculated using the Equation (5), with x i for the i-th secondary dynamic. Note that j = i because one of the process dynamics is the main one: In addition, the new sensitivity added to the cost function proposed in [27], which is related to the main dynamic and the controllable parameters, is formulated using the Equation (6).
It is noted that the stability of the plant as a whole is assumed as guaranteed provided all the local controllers are stable. Therefore, the design of these controllers should guarantee their stability, considering that the variations on their set point are limited by operation conditions, transport parameters, capabilities of the final control elements, and physical and chemistry properties, looking for guaranteeing that the control actions of local controller are into the stable region of the process plant considered as a whole. The superiority of the collaborative proposal consist on using PID as local controllers. In this way, a wide list of PID design techniques guaranteeing closed loop stability are available in published works.
Additionally, some characteristics of the objective function are: (i) The effect of disturbances is considered because the sensitivities S i and Ss k are computed at each sampling time; (ii) The sensitivities are algebraic expressions; for this reason, they are calculated with current measurements. In this sense, the measurements include the dynamic time evolution of the system, and a multi-step prediction as in Model Predictive Control is not necessary; and (iii) the process model is transformed considering explicit differentiation, while the process information is maintained [6].
Step 10: Establishment of triggers. This step determines some triggers for the optimization problems. When the main dynamic is near enough to its set point, the collaborative action is not required, being the main dynamic control loop which provides all control effort. Therefore, either the maximum set point deviation for the main dynamic, or the maximum cost value deviation in some number of past time steps, are useful as particular triggers [6].

Model Predictive Control (MPC)
In order to assess the results provided for the methodology proposed in this paper, an MPC has been considered in the supervisory layer, in particular, a centralized nonlinear MPC that provides the set points for the secondary layer, with the following objective function: Here, y and u are the outputs and control variables, r are the reference signals, P , Q , R are tuning matrices to guarantee stability and to weight the outputs and the control movements, respectively, H c is the control horizon, H w and H p are the initial and final prediction horizons, u min , u max , y min , y max , ∆u min , and ∆u max are bounds to be imposed on the system variables [28].
It should be noted that the MPC uses a receding horizon strategy with an internal prediction model. The control actions are calculated to minimise an index involving set point and predictions deviations in addition with the control movements over two horizons. In this work, the MPC scheme is based on the use of the following nonlinear discrete time prediction model: where k is the current sampling point, and y(k + i | k) are the future control moves at time k + i, with measurements up to time k.

BSM1 Process Description Including the Mass Transfer Model
The simulation benchmark used in this work is presented below together with some specific modeling of some process belonging to the aerobic section. In particular, phenomenological based semi-physical models, based on conservation principles, with constitutive and assessment equations for model parameters are used [23].
A plant for water treatment (WWTP) is formed by three main sections beginning with the primary treatment, which is included to reduce suspended solids and large pollutants, preparing wastewater for a secondary treatment. This last one comprises anoxic and aerobic zones. The last step is the tertiary treatment, for additional improvement of the effluent quality [29].

BSM1 Model
The BSM1 is a well known simulation benchmark for the biological process taking place in secondary treatment of WWTPs. It includes some influent sets for different weather conditions (dry, storm, and rain weather) and defines some performance indexes. The plant layout comprises two anoxic reactors where no dissolved oxygen is present (only bonded oxygen is possible), and three aerobic reactors (Figure 3). Finally, a secondary settler separate sludge from the clean water [30,31]. The aeration in the aerobic tanks are required for nitrification process and it is also relevant to determine the microorganisms metabolic route to degrade pollutants [31]. Therefore, the oxygen concentration is a key variable that must be maintained at a proper level. Note that that oxygen transfer rate depends mainly on reactor geometry, aerators arrangement, sludge characteristics, stir rate and air flow characteristics as a dispersed phase.
The BSM1 model has thirteen states: soluble and particulate inert organic matter (S I , X I ), readily and slowly biodegradable substrate (S S , X S ), active heterotrophic and autotrophic biomass (X BH , X BA ), particulate products from biomass decay (X P ), dissolved oxygen (S O ), nitrite and nitrate nitrogen (S NO ), NH + 4 plus NH 3 nitrogen (S NH ), soluble and particulate biodegradable organic nitrogen (S ND , X ND ). Then, mass balances are obtained for each state, forming the mathematical model presented below, where the subindexes represent the BSM1 state and the reactor number, respectively, w represents mass fractions,ṁ k is mass rate in the flow k, V is the tank volume, S O and S * O are the oxygen and the oxygen saturation concentration, respectively, r i is the reaction rate, and M i is total mass for component i:

•
Mass balance for reactor number one: • From the second anoxic reactor to the fifth reactor (second reactor as an example): • Special case for oxygen (fifth reactor as an example): The main constitutive equations for this mathematical model are all the reaction rates and the mass transfer equation. The reaction rates r i,k are conformed by biological processes and other parameters, and the details are given in [31]. The accuracy of the BSM1 model considered as a case study to apply the methodology proposed is based on the wide use of this benchmark in the literature in the last decade, and its recognition by practitioners in the wastewater treatment field. There are a large amount of research work supporting its validity, and there are also some extensions, such as BSM2, enlarging the activated sludge process with the primary and sludge treatment in a WWTP. Moreover, this benchmark is based on the well known ASM1 model (Activated Sludge Model No 1), developed by [30].
It is important to describe the operational objectives of BSM1, which are considered to evaluate the performance of the methodology proposed [31]. Particularly, in this work, pumping energy (PE), aeration energy (AE), and effluent quality (EQ) are evaluated. Additionally, ISE (Integral of Squared Error) and IAE (Integral of Absolute Error) are considered: (β p,SS · SS e (t) + β p,COD · COD e (t) + β p,Nkj · S Nkj,e (t) + β p,NO · S NO,e (t) + β p,BOD5 · BOD e (t)) · Qe(t)dt where Θ is the period of observation, t 0 is the initial time, Q a and Q w are the internal and purge recycling flows, S * O is the oxygen saturation concentration, V k is the volume of reactor k, k L a is the volumetric oxygen mass transfer coefficient, β p,i are weighting factors for the different types of pollutants, SS e represents the suspended solids, COD e is the effluent chemical oxygen demand, S Nkj,e is the Kjeldahl nitrogen concentration, S NO,e is nitrite and nitrate concentration, BOD e is the biological oxygen demand, and Q e is the effluent flow.

Mass Transfer Model
The mass transfer coefficient represents a constitutive equation of the model, although this coefficient is not modeled in detail in BSM1. There it is only represented by a volumetric oxygen mass transfer coefficient (k L a), and it is calculated only with dependence on the temperature and some other empirical parameters. Particularly, it is represented by Equation (17): where T is temperature in the system and α, F and θ are some empirical parameters. Obviously, this expression for k L a does not represent faithfully the real process because it does not describe other physical phenomena with a paramount influence on the mass transference, such as the size of bubbles and the gas hold up [32,33]. Moreover, in real practice, it is not possible to manipulate the k L a because it does not represent a physical handle. In fact, k L a is the product between the liquid phase mass transfer coefficient (k L ) and the specific interfacial area (a). In order to understand the gas-liquid mass transfer procedure, it is necessary to separate these parameters because they are differently affected by manipulated variables. In this work, an improved model for k L a has been developed and it is presented below. The bubble formation is an equilibrium between gas pressure (expanding the bubble) and superficial tension (aiming to burst the bubble). This formation takes place in two steps [34], the expansion step and the detachment step. Note that the equilibrium not only depends on the physico-chemical properties of the materials, but also on the energy applied to stir the fluid [35]. The bubble formation is relevant because the size and quantity of bubbles define the oxygen transfer in the aerobic reactors. Bubble diameter, d b is obtained by (18) [36], where d 0 is the orifice diameter of the diffuser membrane (m), σ is the surface tension (N/m), ∆ρ is the difference of the liquid and gas densities (kg/m 3 ), and g is the gravity: ∆ρ is calculated according to: where ρ L is the liquid density, considered here as pure water density, and ρ G the gas density, obtained using the ideal gas equation with the internal pressure of the bubble P int : with Wm the molecular weight of the gas, R the Ritcher constant, and T G is the gas temperature. The bubble internal pressure (Pa) is calculated in (21) with the gas injection pressure and the gas drop pressure in the diffuser orifice: where P inj is the injection pressure and ∆P or the drop pressure in the diffuser orifice [37]: with Q air the gas flow rate through the orifice and k 0 the orifice constant (m 3 /Pa 0.5 s). This constant is a function of gas flow rate, density (ρ g ), and viscosity (µ g ), orifice plate thickness (L), and orifice diameter (d 0 ). It is evaluated with the following expression: where C g is a dimensionless constant: with f F the Fanning friction factor. For laminar flow, this factor is obtained as: where Re 0 is the Reynold number. Then, the mass transfer coefficient k L a is obtained considering the liquid-phase mass transfer coefficient k L and the gas-liquid interfacial area a separately. As this coefficient quantifies the velocity of species movement across an interphase [38], it is useful to estimate the mass transfer rate. The k L a dependence on the presence of pollutants is very high [38,39]. Moreover, k L considers the effects of the liquid movement around the gas bubbles and the specific area a reflects the behavior of the bubbles [40]. In [32,41], some empirical correlations to obtain k L are given. In this case, Kawase correlation has been selected because it has been applied in similar conditions to BSM1, and it complies with the required operation range for BSM1 [32]: where D L is the gas diffusivity in liquid (m 2 /s), µ L the liquid viscosity (Pa /cdot s), U G the superficial gas velocity (m/s) and g the gravity constant (m/s 2 ). The superficial gas velocity is obtained as volumetric flow rate (Q air ) passing through the active cross sectional area of the diffuser. In turn, this area represents the empty spaces of diffuser, calculated from the membrane porosity. The diffusers chosen in this work are membrane type (wastewater treatment plant El Marin in Salamanca, Spain), due to its high available active area (A active ). Then, the superficial gas velocity is: As for the specific area (a), it can be obtained as the interfacial area (gas-fluid) per unit of volume [42]. Then, a is obtained according to (28) given by [43], where the bubbles are considered as spheres [44]: where G is the gas hold-up and d b is the bubble diameter. The gas hold-up can be seen as the gas retention in the fluid. More precisely, G is affected by the diameter of the diffuser orifice, liquid density, superficial gas velocity, and surface tension [45]. There are many empirical expressions to obtain G because it affects the mass transfer directly [46]. One of them is [32], where different systems and fluids have been considered. Here, the physical conditions are similar because a column was considered to represent a vertical slide of an activated sludge reactor. However, according to the operating conditions and the of mass transfer taken place, the following equation is chosen: The model of mass transfer given by Equations (18) to (29) is used to design and implement separate control loops manipulating the specific area and the liquid phase mass transfer coefficient. This submodel is essential for the development of the methodology presented in the next section.

Application of the Methodology to BSM1
The collaborative control strategy previously described is applied to BSM1 oxygen control in this section. Note that, for the oxygen mass transfer, separated expressions for k L and a are considered.
Step 1. Development of a PBSM. The wastewater treatment plant is simulated in accordance with BSM1 model equations where a high variability of the influent dynamic characteristics (flow and concentration variables) is implemented [31]. It should be highlighted that proposed collaborative control approach is focused on the aerobic section, for this reason, the mass transfer model presented previously is required.
Step 2. Manipulated variables, states and parameters. According to BSM1, the states considered are thirteen, corresponding to the compounds studied in BSM1. Moreover, there exist many parameters in this process: kinetic and stoichiometric parameters, physicochemical properties as density and viscosity, temperature, design parameters as volume, diffuser specifications, and finally k L and a. The manipulated variables are related with recirculation flows (Q r and Q a ), inlet flow (Q 0 ), air flow (Q air ) and air injection pressure (P inj ).
Step 3. Selection of Controllable parameters. In BSM1, among other existing parameters, k L and a in the aerobic reactors are controllable ones because they are process parameters that can be directly manipulated using a control action. In this sense, k L and a can be controlled by air flow (Q air ) and air injection pressure (P inj ) in each reactor, respectively.
Step 4. Operation conditions and range of variation of states and inputs. The operation conditions are fixed by BSM1 conditions. In this way, the set point for oxygen concentration and total nitrogen concentration are taken form [31]. In this sense, in the second anoxic reactor, the reference for the nitrite+nitrate concentration is stated in 1 g/m 3 . The oxygen concentration in aerobic section is 1.63, 2.47 and 2 g/m 3 in the third, fourth and fifth reactor, respectively.
Step 5. Dynamic hierarchy. The design of the collaborative control structure requires the definition of the dynamic hierarchy. This is determined in accordance with the state impactability index SI I x k . The state (x k ) presenting the highest SI I xk defines the main dynamic while the rest of the states are the secondary ones. The SI I xk values are shown in Table 1. According to the values presented in the table, the process main dynamic is the oxygen concentration being the secondary ones defined by the other variables considered in the table.
Step 6. Control objective. The oxygen concentration is the main dynamic in each reactor; consequently, the collaborative control objective is to regulate the oxygen concentration on a prescribed value defined by the set point. For this reason, the oxygen concentration is regulated using the collaborative control strategy while nitrite + nitrate concentration is controlled by means of an independent PID controller. In addition, other process goal is to fulfill the environmental regulations by keeping the total nitrogen, COD, TSS, and oxygen concentrations below some well established legal limits. Step 7. Input-output pairing. The control methodology is implemented in the three aerated reactors. In this sense, the system has two controllable parameters in each reactor, k L and a. The main goal of the coordinator control implemented in the supervisory level is to control the oxygen concentration, this being a fixed objective. Moreover, a pair of simple controllers, in the lower level, cooperate to fulfill this goal. These two simple control loops are related to the parameters k L and a. The first ones changes the air flow and the second the pressure of the injected air to control k L and a, respectively (see Figure 4). As mentioned before, to keep nitrite and nitrate concentration at a prescribed value given by the set point, an independent PID controller is used, where the manipulated variable is the internal recycling flow. Step 8. Preliminary control tuning. The preliminary control tuning is focused in a regulatory layer. PIDs controllers are tuning using the values given by the Ziegler-Nichols method as seeds, but with a further fine manual tuning. The PIDs controllers are implemented using: u(t k ) being the manipulated variable at time t k , K p the proportional gain, e(t k ) the difference between set point and process output variable at time t k . T i is the integral time, T D the derivative time, and t s is the sampling time. Table 2 presents the corresponding values for the current application. Step 9. Optimization problem. Supervisory layer is based on one optimization problem. In this sense, Equation (4) is implemented taking into account the dynamic hierarchy established in step 5. The optimization problem is reduced to Equation (31). In this case, only controllable parameters are in the objective function and secondary dynamics are not taken into account because they do not have available manipulate variables in order to achieve the control objective: The tuning parameters used are α, β k L and β a , with values 1, 1 and 0.1, respectively. It is important to highlight that S O * ss is the stationary value of oxygen concentration variable when the lower layer control loops maintain the k L and a at the values evaluated by solving the optimization problem previously stated, keeping the rest of process states at their current values. By solving the optimization problem defined above, the set points for k L and a are evaluated and the PIDs manipulate the corresponding control variables to respond to set point changes. Note that S k L and S a are calculated by means of Equations (32) and (33). These expressions were obtained from the process model (BSM1) as follows: Step 10. Triggers of optimization. In this case, deviations greater than 5% on the main dynamic (O 2 concentration) were defined as optimization triggers.

Results and Discussion
In this section, some results of collaborative control applied to BSM1 are presented and compared with the use of an MPC strategy in the supervisory layer. The controller structure for collaborative control is shown in Figure 4. In the case of the MPC controller, the controlled variables are the oxygen concentration in the fifth reactor and the nitrate+nitrite concentration in the second reactor, while the set points of k L 5 , a 5 and the internal recycling flow Q a are the manipulated variables. Then, the MPC controller structure is analogous to the collaborative control one, except for the case of dissolved oxygen in reactors 3 and 4, which have been removed from the MPC structure by simplicity. Note that, in both strategies, there are PID local controllers for regulation in the lower layer, with P inj , Q air and Q a as manipulated variables. For all simulations, the dry weather influent of BSM1 has been considered, having a daily peak of flow and pollutants that is somehow reflected in all concentrations throughout the plant. The simulation time is eight days.
Regarding the MPC tuning, although there are many tuning methodologies available, in this paper, the focus is on the new collaborative control strategy developed, and therefore the MPC tuning has been performed systematically only by trial and error. As for PID tuning, the Ziegler-Nichols method has been applied, with a further fine-tuning to prevent aggressive behavior. This is the current industry practice, taking Ziegler-Nichols settings as seed parameters but applying a trial and error procedure to obtain the final tuning.
In Figure 5, the dynamic evolution of the controlled variables for both strategies is presented, specifically the oxygen concentration in the fifth reactor and the nitrite+nitrate concentration in the second reactor, showing good set point tracking with some periodic oscillations due to the daily high variability of the plant influent. Although MPC presents better tracking for nitrite+nitrate concentration, collaborative control presents better behavior for oxygen concentration because this is the main dynamics in this strategy. Moreover, in collaborative control, the PID for nitrite+nitrate concentration is independent of the collaborative loops. In Figure 6, the oxygen concentration in the third and fourth reactors using collaborative control is presented, also with good tracking.

Oxygen concentration (reactor 3)
Set point Oxygen concentration in third reactor (b) Oxygen concentration in reactor 3. The set point tracking of controllable parameters (Figure 7) is also good, following the references given by the supervisory layer. Additionally, in Figure 8, the control actions P inj , Q air are presented. Note that the set points provided by the MPC are more aggressive than those provided by the collaborative control, due to a more frequent update performed by the MPC optimization problem, while in collaborative control only the optimization triggers command the set point change by solving Equation (31). Finally, the internal recirculation flow Q a is similar for both strategies, due to the relative independence of control loop in the second reactor for collaborative control strategy.         In Figure 9, the ammonium in the effluent is presented, where, in the case of collaborative control, some violations of the legal limit of 4 mg/L can be seen, but only in very short periods of time. On the other hand, for the MPC, there are no violations of the limit, but at the expense of larger variations in the oxygen concentration ( Figure 5b) required in the biological processes to decrease the effluent ammonium concentration.  For a better comparison of the proposed methodologies, Table 3 presents the following performance indexes: ISE and IAE for oxygen set point tracking, energy for pumping (PE), energy for aeration (AE), and the quality of the effluent (EQ). Moreover, these indexes are compared with [18], where an economic oriented MPC controller and a default PI strategy without supervisory layer [31] is implemented, in order to evaluate the general performance of the methodology proposed in our work. From these indices, it can be seen that collaborative control strategy presents better ISE and IAE than MPC controller for oxygen control, as can also be seen in figures with the dynamic evolution of this concentration. On the other hand, PE and AE are better in MPC due to smaller control actions required by the supervisory layer in this case. Finally, EQ is shown here as an additional index for performance, though it has not been included directly in the methodology. Although the EQ value for collaborative control is a bit higher than for the MPC structure, it is still acceptable, particularly if it is compared to [18].

Conclusions
This paper proposes a novel collaborative control structure taking into account in its design the plant dynamical characteristics, establishing a dynamic hierarchy by using the Hankel matrix of the process model and the control of specific parameters by means of simple PID control loops, which is especially attractive for industrial applications. The methodology proposed has been tested in the control of a specific process of wastewater treatment plant, particularly using the BSM1 benchmark.
The good results obtained for set point tracking show the efficiency of the proposed methodology, when it is compared to other structures reported in the literature such as the ones based on MPC. The method complexity and the computational load reductions are some of the most remarkable characteristics, in comparison with MPC, which makes the presented proposal very suitable to be applied to large-scale processes. This computational reduction is mainly due to the fact that the calculation of inputs and states predictions over some time horizon is not further needed. In the proposed methodology, the decision variables are evaluated without solving the complete model at each sampling time.
Although it might be possible to achieve better performance with the MPC due to higher complexity of the optimization problem, the results obtained should be evaluated considering the difficulties with MPC tuning. In fact, one of the main advantages of the collaborative control strategy is the simplicity, including an easier tuning because only weights α and β must be selected. Note also the particular difficulties in considering BSM1 process as a case study due to remarkable different dynamics (nitrogen, oxygen) and disturbances. Moreover, the inclusion of the nitrate+nitrite concentration as a controlled variable within the multivariable MPC framework further increases its complexity and therefore its performance is worsened.
The collaborative control outperforms the conventional MPC/PID control structure depending obviously on the process considered because both methodologies have their advantages and drawbacks. Anyway, a general criterion to apply the collaborative control is that it may be valuable for industrial complex processes with difficult dynamics, where the use of other prediction-based controllers (such as MPC) need a time-consuming tuning and an optimization routine executed each sampling time.
Finally, it is important to mention that the use of a detailed model for the mass transference is also a novelty of the work. The advantage compared to other simpler models is the improved description of the real phenomenon making the choice of controllable parameters easier. In addition to this, it represents the relationship between controllable parameters and manipulated variables in a straightforward way, increasing the degrees of freedom to control the process. The validation of this mass transfer has been performed in comparison with the model proposed in [31], with adequate results.