Distributed Model Predictive Control Applied to a Sewer System

In this work, a Distributed Model Predictive Control (MPC) methodology with fuzzy negotiation among subsystems has been developed and applied to a simulated sewer network. The wastewater treatment plant (WWTP) receiving this wastewater has also been considered in the methodology by means of an additional objective for the problem. In order to decompose the system into interconnected local subsystems, sectorization techniques have been applied based on structural analysis. In addition, a dynamic setpoint generation method has been added to improve system performance. The results obtained with the proposed methodology are compared to those obtained with standard centralized and decentralized model predictive controllers.


Introduction
Currently, it is a fact that urban wastewater treatment plants (WWTPs) have a high degree of interconnection with other elements, such as sewers, conduction systems, and storage tanks, which are part of the urban water system. Urban drainage systems (UDS) generally catch and carry both urban wastewater and that which comes from rainfall to WWTPs for treatment before being discharged into the environment, constituting a combined urban drainage system. During periods of heavy rain, the wastewater resulting from the mix can overload the urban system and produce overflows that can be harmful to the environment. Until now, UDS and WWTPs in general have been considered separately to solve their control problems.
In the case of UDS, the research has essentially consisted of the development of detailed process models that facilitate the design and simulation of different control strategies. The operational objectives are basically avoiding wastewater leaks from the collection network before being treated by the treatment plants, due to overflows in the storage tanks and at the entrance of the WWTP, optimizing the degree of use of the station by maximizing the inlet flow, and reducing the economic cost of operation [1,2]. Other problems that appear are the "first flush" problems, which occur after a dry season when abundant precipitation occurs, causing a large increase in the concentration of pollutants in the network water.
Some simulation models constitute benchmarks for the control systems of the sewer, such as [3]. On the one hand, the models try to describe with great precision the processes that take place in a drainage network (hydraulic and chemical processes of degradation of discharged pollutants) and, on the other, they have to be simple enough to be the target of simulations in the long term. In the development of these models, the characteristics of the wastewater collected according to the area (urban and industrial waters with different concentrations of pollutants), the rainwater reaching the network, and the time of year are taken into account, based on a large volume of collected data. As a result, typical patterns that attempt to cover a large number of situations have been developed [4]. From these models, simpler ones can be obtained to be used in the elaboration of certain types of control algorithms.
Regarding the control systems for UDS [5], they are usually classified into off-line controllers, which apply static rules, and improved on-line controllers that apply real-time control actions [3,[6][7][8]. Table 1 shows a summary of features for the real-time control systems: In particular, simple RBC (Rule-Based Control) strategies have been applied in [9,10], offering as main disadvantage the increase of rules when the complexity of the system increases. Alternatively, FLC (Fuzzy Logic Control) based control systems that combine simple rules with an expert system and a flexible specification of the output parameters can be used [9,[11][12][13][14]. An FLC is also used to reduce the volume and number of overflows in the Wilhelmshaven (Germany) UDS achieving both objectives [15]. The principles of FLC have been used to control the operation of the pumps in the collection network od the city of Taipei, achieving a more efficient drainage of rainwater, and thus preventing flooding [16]. Another type of algorithm, the LQR regulator, has been applied in [17] to flow control in a collection network in Bavaria (Germany). In this case, the control objectives were to minimize the overflows in the system by optimally using all the storage capacity of the UDS and emptying the network as soon as possible. Evolutionary Strategies (EA) mimic evolutionary principles to find optimal solutions [18]. For example, Genetic Algorithms have been applied to water quality management systems in [19] and to the optimal multiobjective control of UDS [20]. In [10], EA combined with self-adaptive algorithms has achieved an improvement in the quality of the water discharged into rivers, with reduced costs. In the control techniques based on Population Dynamics (PD), the designed controller poses a problem of resource location, and its interpretation can be related to biological evolutionary processes, without considering a system model, using game theory concepts. These strategies have been applied in [21,22], achieving a better use of network volume and avoiding overflows.
Regarding the architecture of the control systems used in UDS, one of the most common global control configurations is centralized control. However, the local control scheme may be a suitable solution if there are few actuators in the system [23]. In complex large-scale systems, it is common to have both levels of control (global and local). In this case, there can be up to three levels of control configuring a hierarchical control structure [24,25].
In the urban water systems, the interaction between the UDS and the WWTP is clear, since the outflow from the UDS is the inlet of the treatment plants. Consequently, to improve performance, the development and evaluation of control strategies for WWTPs should be addressed from a broader perspective [26]. However, there are only few works in the literature that take both systems into account in an integrated way and focus on the performance of the WWTP, using simplified models for the sewer networks [27,28]. In particular, the consideration of the integrated system allows for a bidirectional coupling of the conceptual models of the subsystems including WWTPs [29]. Regarding the applied control strategies, in [30], a heuristic-type control system based on rules that covers both situations of scarcity and abundance of water to be treated is described, also taking into account the operating costs of the entire system. Others, like [31], present strategies that combine RBC algorithms with Model Predictive Control (MPC).
MPC is a control methodology that uses a process prediction model to calculate the manipulated variables on a future horizon in order to optimize a certain cost function [32][33][34][35][36]. An MPC controller has four basic elements: a control-oriented mathematical model of the system, a cost function that expresses the control objective to be achieved, a set of system constraints, and a finite horizon optimization problem that is solved in a receding-horizon way [33]. In a hierarchical control structure, the MPC is normally used to generate optimal references for local controllers. The characteristics of MPC controllers have certain advantages for their application to UDS, particularly the ability to anticipate the system's response to future rainfall events and the ability to consider delays and disturbances [32]. This controller has been applied to the UDS in [25], where a centralized predictive control has been developed according to a proposed model in which the state variables considered are the volumes of the network tanks, the control or manipulated variables are flow rates in the network pipes, and the measurable disturbances are the flows of rainwater. The restrictions in this case are given by the available volume of tanks and pipes, among others.
A prioritization of multi-objective cost functions is applied by means of the lexicographic technique in the control of a UDS in [37]. In [38], an MPC controller has been simulated in a part of Barcelona's drainage network, in Spain, achieving significant reductions in floods and overflows [25,[38][39][40]. Other cases in which MPC control algorithms have been studied and/or applied to UDS are [23,41,42], where both operation costs and overflows have been reduced using non-linear prediction models.
In large-scale systems, it may be advisable to divide the global process into simpler subsystems to facilitate the application of MPC algorithms [7]. In this case, prediction models and local cost functions are used, obtaining local solutions to the global control problem. If there is no exchange of information between subsystems, Decentralized Model Predictive Control is a possible solution, but a more efficient alternative is the Distributed Model Predictive Control (DMPC), where local controllers exchange information to calculate their own local solution, as in the cooperative (or Coordinated) Distributed Model Predictive Control [8,[43][44][45]. Although there is a large amount of DMPC strategies applied to different processes [46][47][48], only a few deal with water management, such as the level control in tanks [49] and the coordination of drinking water supply networks [50][51][52]. To the knowledge of the authors, no applications to sewer systems have been developed.
The main contribution of this work consists in the development and application of a practical cooperative distributed model predictive control to a UDS, based on local linearized models of the system and fuzzy negotiation among subsystems [53]. Another contribution is the inclusion of the WWTP in the control strategy as an additional objective-more concretely, the optimization of the WWTP inlet flow. The benchmark described in [3] has been considered, with different disturbances, and the results have been compared with a centralized and decentralized MPC, to assess the suitability of the proposed methodology and its application to this system.
The article is structured as follows. After an introduction, the description of the system under study begins, along with the mathematical model used for the simulation and the model used by the control algorithm. In the following section, its sectorization is detailed. The article continues with the description of the centralized and distributed MPC control algorithms, showing the corresponding results in each case in two representative scenarios. The conclusions are presented at the end.

System Description
The sewer system used as a benchmark [3] (Figure 1) is made up of six rainwater and wastewater catchment areas (numbered from 1 to 6 in the figure), six wastewater storage tanks (ST1, . . . , ST6; one of them, ST5, is off-line), wastewater pipes, five valves and a pump for flow control, and a wastewater treatment plant (WWTP). The objective of the sewer is to collect all the residual water and lead it to the WWTP, maintaining a supply flow with the least possible variability and as close to its nominal value as possible. This is achieved by retaining the volume collected in the tanks and releasing that volume when the collected flow is small. In addition, efforts are made to reduce as much as possible the overflows in the reservoirs in the face of extraordinarily intense precipitation, which can cause flooding and the dispersion of polluting substances potentially harmful to the environment.
The simulator has algorithms to generate different scenarios that take into account both the characteristics and the volume of urban wastewater discharged into the network and the runoff and filtrations collected by the network, produced by rainfall of different intensity and duration, in relation to the day of the week or season of the year considered.
To allow for the design and validation of the control algorithms, a simplified model has been developed. Firstly, only the hydraulic part of the benchmark has been considered. The ST5 tank has been removed, since, in the current configuration, it is not possible to exercise control over it. Thus, the ST6 tank has been renamed as tank 5 in Figure 2. The ST4 tank is considered to be equal to the rest, that is to say, the inlet valve is fixed so that all the water that enters the tank and the outlet flow is produced by opening another valve, rather than activating the pump. The simplified mathematical model of the process is made up of the following elements [54]: -Water catchment areas: all the water collected in the area constitutes an inlet flow to the system that is treated as a disturbance. The flows collected in each zone are denoted by q ri , i = 1, . . . , 6. -Link elements: gravity residual water pipes in open channels, which connect the flows collected in each area with the tanks and the tanks with each other and with the treatment plant. They are modeled as first-order systems with very slow dynamics, according to their length. Their discrete mathematical model is the following: where: q i (k) is the output flow of the element i τ i is the time constant of the element i T is the sampling period q u,i (k) is the sum of inflows to the link element i -Storage tanks: reservoirs for storing wastewater. They are located in specific places throughout the network, with the possibility of overflow if a maximum level is reached. Their discrete model for the water level is: where all parameters are referred to tank i and instant k: is the volume in use with water V max,i is the maximum volume of the tank u in, i (k) is the input flow rate q ov, i (k) is the overflow flow rate d 0i is an experimental parameter related to overflow u i (k) is the output flow rate h i (k) is the water level h max,i is the tank height A i is the tank area v i (k) is the opening of the tank outlet valve (control variable: v ∈ [0, 1]) c 0i is an empirically calculated discharge coefficient for each tank that depends on the section of the outlet pipe -Nodes: they represent places of confluence of various residual water pipes. The resulting flow is considered the sum of the incoming flows.
In Figure 2, the complete system is shown by means of a block diagram built from the basic elements, where the catchments are represented by ovals, the link elements by rectangles, and the storage tanks by triangles. This system is analogous to the benchmark presented in Figure 1.
Then, the mathematical model corresponding to the diagram in Figure 2 was obtained. The states considered are the levels of the deposits 1, 2, 3, 4, and 5 (x 1 , . . . ,x 5 ) and the flow rates of the link elements 3, 7, 8, and 9, which correspond to the states (x 6 , . . . ,x 9 ). The output flow of link element 9 (state x 9 ) is in fact the inlet flow to the WWTP, provided that there is no overflow, that is, as long as it does not exceed the nominal value of the inlet flow to the WWTP. The output flows of the link elements of the catchments 1, 2, 4, 5, and 6, and the flow collected in zone 3, will be considered as measurable disturbances on the process: (d 1 , . . . ,d 6 ). The overflowed volume in tank 1 is collected again by the network and led to tank 4 thanks to the arrangement of the tank itself and that of the sewage network. The manipulated variables are the desired flows at the tank outputs, that in turn would constitute set points of local PID flow regulators at a lower level: (u 1 , . . . ,u 5 ). All states will be considered as outputs of the system, with a special interest in the input flow to the treatment plant. Hence, the states vector, the manipulated variables vector, and the disturbances vector are defined as follows in terms of the benchmark variables: The set of discrete differential equations of the non-linear mathematical model is: From this model, a linear model of the system is obtained removing the overflow flow rate terms, to be used as a prediction model in the linear MPC algorithm. The state space model equations are the following: And the matrices of the model are:

Evaluation Criteria
The following evaluation criteria are used to analyze the behavior of the system and the effects of the applied control strategies. The evaluation considers different points of overflow in the sewer network and at the entrance of the treatment plant [3]. Some of them only consider hydraulic magnitudes, while others also consider variables related to water quality. Among the indices corresponding to the first group, the most used are shown below, where the subscript i represents the specific point in the system:

1.
Number of overflows (N ov,i ): this criterion represents the total number of overflow events that have taken place at a certain point in the network. Two overflow events that are less than an hour apart are considered a single event.

2.
Duration of overflow (T ov,i ): this criterion represents the cumulative sum of the duration of all overflow events at a certain point in the system at simulation time T sim (d). If each event j lasts t j days and n events take place, the overflow duration is given by the expression: 3.
Overflow volume (V ov,i ) (m 3 ): it is the total volume of wastewater discharged into the receiving waters receivers from a certain overflow point i. Considering simulation time T sim and overflow q ov,i (t): and the total overflow volume is: 4. Degree of usage of WWTP (G u,i ) (%): the ratio of the average flow rate in the simulation time and the nominal flow rate q WWTPnom,i of the plant i, in percent:

5.
Smoothness in the application of control signals (flows) on the system (m 6 /d 2 ): it can be measured from variations in control actions u(k) between the instant k and k − 1 in the simulation time, by means of the following expression:

System Sectorization
Applying control techniques to large-scale systems usually involves dividing the global system into simpler subsystems. To this end, a structural analysis of the system has been carried out to determine the best way to sectorize the system, so that the subsystems obtained have a minimum degree of coupling and are controllable [55]. With this end, the reachability from the entrance of the system has been checked, and the direct graph of the system has been obtained. Definition: A system S with a direct graph D = (U∪X∪Y, E) is input-reachable if X⊆ n is a reachable set of U⊆ m , that is, for each node x∈X there is at least one node u∈U with a direct path to it. Definition: A system S with a direct graph D = (U∪X∪Y, E) is output-reachable if X⊆ n is an antecedent set of Y⊆ r , that is, for each node x∈X there is at least one node y∈Y with a direct path to it. Then, a system S is reachable if it is input-and output-reachable. The study of the reachability of the system allows us to determine the possibilities of sectoring while maintaining the reachability of the subsystems obtained. It begins with obtaining the interconnection matrix E, which is defined from the matrices that determine the dynamic model of the system, A, B p , and C, as: , where A, B p , and Care given by : Reachability matrix R can be computed as R = E ∨ E 2 ∨ . . . ∨ E s , with s = dimE, and it has the following form: where F, G, H, and θ are binary matrices whose dimensions are consistent with E.
The system will be input-and output-reachable if and only if the binary matrix θ has neither zero rows nor zero columns [55]. For the sewer system considered here, the following matrix θ has been obtained, resulting in a system that is reachable from the input and the output: Moreover, obtaining the direct graph of the system ( Figure 3) will allow us to know the interactions between the different variables that are involved in the process and determine the best way to divide the global system into subsystems with minimal interaction [36], maintaining the reachability of the subsystems obtained.
The graph shows a system coupled by the inputs. The sectoring has been carried out avoiding interactions between states of different subsystems and at the same time ensuring their reachability ( Figure 4). The controllability of each subsystem can be verified. Accordingly, this work considers two subsystems, the first one englobing the tanks 1 to 4 together with link element 3, and the second one tank 5 and link elements 7, 8, and 9.
To apply the distributed MPC developed in Section 7, the state space local models of each subsystem are expressed as follows, where the coupling input u 4 is considered to belong to subsystem 1.
where : where : Matrices A 1 , B p11 , B p12 , A 2 , B p21 , and B p22 are formed by suitably selecting the rows and columns of A and B from the matrices (6). Note that B p12 and B d12 are null; B p21 only has non-zero in the last column, and B d21 is null, due to the plant's configuration.

Control Objectives
The following evaluation criteria are used to analyze the behavior of the system, and the effect of the main control objective is to ensure an input flow to the WWTP as close as possible to its nominal value, making the most of its capacity, avoiding as much as possible overflows in the tanks and at the wastewater plant itself, and also minimizing operating costs. To achieve the proposed control objective, the outlet flow rates of the tanks (manipulated variables) are calculated to minimize the difference between the flow rate of the inlet to the treatment plant and its nominal value. Moreover, the volume of water at each instant is evenly distributed among all the tanks in the network, which is achieved by minimizing the difference between the water level of each tank and a dynamically calculated reference level to achieve that goal [54]. The uniform distribution of the water stored in the tanks will reduce the effects of the disturbances (collected flows at the catchment areas), minimizing the overflows. Mathematically, the control objective can be formulated by means of a cost function that includes the partial objectives previously exposed [54,56]: where N is the prediction horizon, ϕ j is a partial objective, and w j is the weight corresponding to each partial objective ϕ j , with j = 1, ..., M.
In the problem considered, the following partial objectives are taken into account: • Minimization of overflows and uniform distribution of the stored volume of water: where V G is the total occupied volume in the sewer at instant k and v i is a factor that represents the weight of the tank capacity i in the total available storage volume in the sewer. The variable weight q i (k) allows for the overflows' penalization, increasing its value proportionally to the corresponding overflow. Parameters f i and α i are used to prioritize overflows in some tanks.

•
Maximum usage of the treatment plant and minimum overflow at the WWTP influent: where Q WWTP is the inlet flow to the WWTP at instant k and Q WWTPmax its nominal value. • Control efforts minimization: where u iref are the output flows needed to keep the desired level in a tank, obtained by Bernoulli's law. Weights r i are tuning parameters.

Hierarchical Control System
The control system presents a hierarchical structure as shown in Figure 5 [23,25]. At the upper level, the level set point signals are generated for each tank to achieve the indicated control objectives. This will be done following the strategy of uniform distribution of the total volume occupied among all the deposits, called EFD (Equal Filling Degree) [57], which has been shown to be among the most effective dealing with overflow problems. Subsequently, the MPC controller calculates the set points to be applied to local controllers by solving a constrained optimization problem, maximizing the use of WWTP and reducing operating costs. Finally, the local controllers are in charge of the effective control of the process variables [58].

Optimization Problem
The Centralized Predictive Control considered in this work includes a linear state space model for prediction, considering the effect of perturbations. The cost function (objective function) of the MPC is a quadratic form that takes into account both the tracking errors in the states and the deviations in the control sequence from the flow reference (therefore penalizing control efforts), assuming that the prediction and control horizons coincide and their value is N: where x re f (k) = (x 1re f , x 2re f , . . . , x 9re f ) and u re f (k) = (u 1re f , u 2re f , . . . , u 5re f ) are the states and inputs references, respectively, and U(k) = u(k) u(k + 1) . . . u(k + N − 1) T . The optimization problem for the centralized MPC is: where q maxj and u max j are the upper bounds for flow in the linking elements and the tank outputs, respectively. Matrices Q(k) and R and the horizon N are tuning parameters for the controller, and matrix P is a terminal penalty for MPC stability obtained by means of the Riccati equation [59]. Note that a variable diagonal matrix Q(k) has been considered to somehow prevent overflows in the tanks, excepting tank 1, because the overflow returns to the sewer. If an overflow is produced in a tank, the corresponding weight is increased to further penalize the deviation from the reference for that tank. The non-zero elements of Q(k) are q 1 , . . . , q 5 , and q 9 , and are obtained by (17), where f i and α i are tuning parameters. The rest of the elements are zero because they correspond to state variables of link elements that are not necessary to follow any set point. Matrix R is a diagonal matrix that penalizes the variations in the flow set points with respect to their reference value, which is related to the cost of operation.
This optimization problem corresponds to a quadratic optimization problem (QP) with constraints [32]. The implemented algorithm keeps the last applied solution, so that if the optimization problem is not feasible, the last feasible solution is applied.

Set Point Calculations
To achieve optimal operation, we designed a hierarchical controller that modifies, in a higher layer, the tanks' level set points according to the strategy of distributing the current volume of water among all the tanks in the most uniform way possible considering their capacity [57]. This is achieved by calculating for each one its reference level, based on the total capacity of the network and the capacity of that deposit at each sampling instant: where x ire f (k) is the reference level for tank i at instant k. The references for the states x 6re f ,x 7re f , and x 8re f would be zero because they are outflows of link elements that are not considered. Moreover, x 9re f = 60.000m 3 /d, which is the nominal influent flow to the WWTP. The flow setpoints are calculated taking into account the desired level for each tank, according to expression: The references are updated every two sampling periods, which is sufficient considering the dynamics of the system and the disturbances.

Distributed Model Predictive Control (DMPC)
The DMPC algorithm considered in this work is based on [60], consisting of the minimization of local cost functions that depend on the future evolution of the inputs of each subsystem and the neighbor, using prediction models and local constraints. At each sampling period, agents get an optimal control sequence for their subsystem, keeping the neighbor control sequence constant, and then each agent solves another optimization problem that provides the control sequence for the neighbor, keeping theirs constant.
Due to the particular characteristics of this system, since there is no interaction of the inputs of subsystem 2 in subsystem 1, the procedure is simplified, and some of the optimization problems of the general method are not solved. Figure 6 gives an overview of the algorithm.
And it consists of several steps: STEP 1: The state measurements x 1 (k) and x 2 (k) corresponding to each subsystem are obtained at instant k. Two local predictive control problems (MPC1 and MPC2) are solved, obtaining local solutions U * 1 and U * 2 for each subsystem: MPC1 objective function: The optimization problem is: T subject to : where U S 2 is the previous optimal control sequence extended to the current instant, keeping constant the last value in the prediction horizon.
MPC2 objective function: The optimization problem is: T subject to : where U S 1 is the previous optimal control sequence extended to the current instant, as in (25). STEP 2: • Agent 1 (MPC 1) does not solve any optimization problem because U 2 does not affect this subsystem (there is no coupling). • Agent 2 (MPC 2) minimizes J 2 with respect to U 1 , keeping constant U * 2 , obtained in step 1. Note that u 4 is the only variable optimized, because the rest of the inputs do not affect the state in the prediction model, due to the zero columns in B 21 .
The optimization problem is: T subject to : STEP 3: Information exchange: agent 2 sends U w 1 to agent 1. STEP 4: Feasibility test: It is necessary to check if the states of subsystem 1 satisfy the constraints of this subsystem, provided that U w 1 is obtained by the other agent. If the constraints are not fulfilled, this control action will not be considered for negotiation in the following step. STEP 5: Agent 1 applies a fuzzy negotiation between the two available control signals, which takes into account the average input flow rate to the wastewater treatment plant, q WWTP , and the average overflow, q ov .

Fuzzy Negotiation
The fuzzy negotiation generates the final control actions for the first subsystem U f 1 , because U f 2 = U * 2 has already mentioned. The starting point for the negotiation are two solutions obtained in Steps 1 and 2 of the algorithm, U * 1 and U w 1 . In order to apply the fuzzy criteria, the average input flow to the WWTP, q WWTP,j and the volume of water that has overflowed on average in the prediction horizon q ov,j are calculated for each of the negotiating control actions. For this, local prediction models are used and local solutions are taken into account, as well as disturbances:

Fuzzification
Fuzzification consists of the transformation of imprecise knowledge into fuzzy sets. Specific numerical values of a property are transformed into fuzzy values for each linguistic variable (degrees of membership) to make decisions as humans do [61,62]. In this case, the fuzzy sets proposed for each variable have been obtained in a heuristic way, in order to get the best performance of the system. Therefore, the values for the flow rates a, b, and c that determine the precise shape of the fuzzy sets are the following (in m 3 /d): The graphic representation of these sets is shown in Figures 7 and 8.  The numerical values considered in each case appear in Table 2.

Fuzzy Rules
The evaluation of the rules provides the fitness of a possible control action based on the imprecise knowledge defined by the values of the specific membership functions defined before. The set of rules represents all the possibilities of combinations of linguistic labels in the two antecedents, to provide results for all the universe of discourse:

1.
IF q WWTP is ideal AND q ov is negligible, THEN the solution is excellent.

2.
IF q WWTP is ideal AND q ov is noticeable, THEN the solution is good.

3.
IF q WWTP is low AND q ov is negligible, THEN the solution is good.

4.
IF q WWTP is low AND q ov is noticeable, THEN the solution is poor.

5.
IF q WWTP is high AND q ov is negligible, THEN the solution is good.

6.
IF q WWTP is high AND q ov is noticeable, THEN the solution is fair.
These rules are applied to the two possible solutions available in agent 1 for negotiation. The aggregated antecedent (logical AND) is obtained by using function min. Then, an example of an equation is: where i = 1, . . . 6 is the number of rule, and j = 1, 2 is the solution considered for negotiation (U * 1 or U w 1 ).

Defuzzification
Defuzzification consists of obtaining a characterization of each control action by means of a precise numerical value from the fuzzy values obtained in the evaluation of the rules. There are numerous procedures for this, and in this work it was done by applying the centroid method to the membership sets of Figure 9 where the values of U p and U m are selected taking into account which solution, U * 1 or U w 1 , provides a smaller value of the global cost function: and vice versa if J k, U * 1 , U * 2 > J k, U w 1 , U * 2 . Then, the degree of membership for the consequents of the rules is calculated using the RSS method (Root Sum Squared [61]), considering that the rules are applied to both control actions to negotiate. The degree of membership for the excellent, good, fair, and poor membership sets, are denoted as µ e , µ b , µ r , and µ m , respectively: Finally, U F 1 is obtained as follows (centroid method [61]):

Results
In this section, some simulation results considering two different scenarios obtained from the sewer benchmark are presented. In order to reduce the computational time, two periods of 10 days have been taken out from a series that extends to two years, choosing periods where the flows variations are more relevant.
Both scenarios represent specific periods of a humid season with rainfall ( Figure 10), the second one with more rainfall collected ( Figure 11).  For each scenario, to validate the results of the methodology proposed in the paper, four cases have been developed. The first case (without control, CASE 1) is equivalent to keeping all the valves at their maximum opening. The second case (Decentralized MPC, CASE 2) considers two local MPCs according to the models obtained in (15). The third case (DMPC with fuzzy negotiation, CASE 3) is the methodology proposed in the paper, and CASE 4 is a centralized MPC. The received water comes from the normal discharge of urban wastewater and more or less intense rains, depending on the season, which are the ones that mainly cause the overflow in the tanks and at the entrance of the WWTP. For all simulations, the prediction horizon selected is N = 5, and the weights in the cost function are shown in Table 3. Table 3. Tuning parameters.

Weights for (12)
Weights for (13) Weights for (15) w i = 1, i = 1, 2, 3 f i = 10, i = 1, 2, 3, 4, 5 r i = 10 −8 , i = 1, . . . , 11 α i = 10, i = 2, 3, 4, 5 q 9 = 10 −5 The system parameters have been taken from [3] and are shown in Table 4. Some simulation graphic results are shown below for scenario 1. As can be seen in Figure 12, the inflow to WWTP is very similar (b), as well as the overflow flow rate (a), except in the case without control. The differences between the cases considered are better revealed in the numerical results of Table 5. For brevity, only tank 4 level (a) and outflow set point (b) in each case are shown in Figure 13.  The overflow volume in tank 1 returns to the network due to the sewer configuration. Therefore, it is not added to V ov . The evolution of the water level in tank 4 can be seen together with the reference proposed by the upper layer of the hierarchical controller. The reference tracking is considerably better for cases 2, 3, and 4 than for case 1 (without control). Therefore, if no control is applied, more overflow is produced at the inlet of the WWTP, because not enough water is maintained in the sewer tanks. The controlled cases also improve the total number of overflows at the inlet of the WWTP and the total volume of overflow water (Table 4). Regarding the inlet flow to the WWTP, it is difficult to keep the nominal value q WWTPnom , because the water received in the catchments has large variations and most of the time is not enough to reach that value. However, when control is applied, more water is stored in tank 5, allowing for a better regulation of the nominal value of the WWTP inlet flow (Figure 12a).
In Figure 14, level (a) and output flow rate set point (b) can be observed for every tank in case 3, corresponding to a Distributed MPC with fuzzy negotiation. Levels are represented together with the reference level calculated at each instant according to (22). Some of them are closer than others to the reference value to minimize the volume of overflowed water and optimize the inflow to the WWTP.   The results for scenario 2 present a larger number of overflow events (N ov ) and a larger overflow (V ov ) due to the increase of the water received (disturbances). On the other hand, the usage of the WWTP (G u ) and the Q WWTP average are due to the larger availability of water to keep the WWTP working close to the nominal operation point.
In Table 5 (for scenario 1) and Table 6 (for scenario 2), the performance indices presented in Section 2.2 are presented for the four cases. Case 1 (without control) is the worst for all indices, as expected, and the centralized MPC presents the best performance for indices such as the average flow entering the WWTP (Q WWTP ), the degree of usage of WWTP (G u ), and the total volume overflow (V ov ). However, the DMPC with fuzzy negotiation does not show an important decrease of those indices, with the advantage of using local simpler models and optimization problems. Moreover, the use of DMPC provides the smallest S, due to the smoothing effect of the fuzzy sets considered for negotiation. The case of Decentralized MPC provides the worst performance indices because the local controllers do not exchange information and the existing coupling in the network is not considered.  Finally, to analyze the influence of the fuzzy set design on the performance of the distributed control algorithm, the location of the fuzzy sets has been changed according to Table 7, but their shape was preserved. Results are presented in Table 8. The DMPC1 and DMPC2 cases present practically identical results. Therefore, it follows that the influence of the position variation of the fuzzy sets is negligible. By comparing DMPC2 and DMPC3, it can be seen that the degree of usage of the WWTP is slightly smaller in the case of DMPC3 because the center of the fuzzy set has been moved to the WWTP constraint of 60,000 m 3 /d. By comparing DMPC3 and DMPC4, the overflow (V ov ) is larger for DMPC4 because the fuzzy set considers as negligible overflows of larger amount. Anyway, the results are very similar because the limits of the fuzzy sets are very next to each other. 6.0324 × 10 10 6.0144 × 10 10 6.0583 × 10 10 6.0322 × 10 10 6.0325 × 10 10

Conclusions
In this document, a DMPC with fuzzy negotiation has been developed and applied to a sewer network, with good results in comparison with centralized and decentralized MPC. As expected, centralized MPC provides the best performance because the controller makes use of the full linearized model of the process. However, the DMPC performance is similar, dealing with simpler local optimization problems and other common advantages of distributed strategies, such as fault tolerance. As for the comparison with a decentralized MPC, where there is no communication among agents, the negotiation improves the results significantly. The methodology developed is based on linearized local models, which further simplifies the MPC optimization problems solved by each subsystem. Moreover, the fuzzy negotiation allows for the inclusion of expert process knowledge by means of fuzzy sets whose precise shape is defined by flows and intuitive fuzzy rules that can be adjusted depending on the specific needs of the network. In this way, no other cumbersome tuning parameters have been added to the negotiation procedure.
The methodology proposed here is currently valid only for two subsystems, but it can easily be extended to a larger number of local systems due to the use of standard sectorization techniques and introducing, for example, the concept of neighborhood, where the coupling among subsystems is high. Due to the use of the benchmark [3], it would also be possible to include the concentration of the most relevant pollutants in the procedure. Although UDS control is mainly commanded by flow and level dynamics, due to the irrelevant biological degradation of pollutants in the sewer, future works will consider some of them, e.g., the Total Suspended Solids concentration, with the aim of improving the industrial applicability of the controller.
Author Contributions: P.V. and M.F. conceived and designed the general approach; A.C. and M.F. performed the simulations; M.F. and P.V. contributed to the simulation interpretations; writing-review and editing, M.F. and A.C.; supervision, P.V.; A.C. and M.F. wrote the paper. All authors have read and agreed to the published version of the manuscript.
Funding: This work has received financial support from the Spanish Government through the MICINN project PID2019-105434RB-C31 and the Samuel Solórzano Foundation through project FS/20-2019.

Conflicts of Interest:
The authors declare no conflict of interest.