Reducing the Environmental Impact of Sewer Network Overflows Using Model Predictive Control Strategy

This paper proposes a method for reducing the environmental impact of sewer network (SN) overflows. The main objective of the paper is to minimize the wastewater quantity and the pollutant loads that overflow from the SN. The proposed algorithm to achieve this goal is Model Predictive Control using Particle Swarm Optimization as optimization method. It was tested in simulation using a simplified model of the network based on Benchmark Simulation Modelsewer as prediction model, and a forecasted influent. Three cases have been considered: (a) the fitness function is defined as the global yearly overflow volume calculated using equal weights for each tank; (b) the fitness function uses different weights for each tank depending on the medium loads and (c) integrating a penalty term related to the system state at the end of the prediction horizon in the previous fitness function. The simplified model determined a significant reduction of the integration time minimizing the optimization time.

are integrated multiple times, the advantage of these simplified models consists in the significant reduction of the optimization time.In (Balla et al., 2022;Troutman et al., 2017) the authors present a data-driven identification/ learning toolchain to be used for combined SNs.
Usually, in order to increase the SN efficiency, research has been oriented toward the use of optimization algorithms, by accounting for various performance indicators.The main objective regarding the SN operation is the minimization of its environmental impact by minimizing the quantity of water discharged directly into natural receptors during overflow events, as well as minimizing the pollutant concentrations of the discharged water.The aim is to avoid the situations in which environmentally damaging events (such as floods or soil/water pollution) occur.At the same time, the SN effluent, which is influent of the WWTP, must be "calibrated" so that it operates with a high efficiency.Thus, the papers (Joseph-Duran et al., 2014, 2015;Ocampo-Martinez et al., 2007;Sun et al., 2021;Troutman et al., 2020) propose predictive algorithms in order to minimize in real time the amount of polluted water.A series of papers (Minzu & Serbencu, 2016;Minzu et al., 2015;Ochoa et al., 2019;Rathnayake & Faisal Anwar, 2019;Safavi & Geranmehr, 2017;Vasiliev et al., 2022bVasiliev et al., , 2022c) ) presents different optimization algorithms such as Particle Swarm Optimization (PSO), genetic algorithms, FminSearch optimization etc. aiming to simultaneously minimize two objectives: the pollutant loads in the wastewater discharged from the SN directly into the environment and the cost of the WWT and of the energy consumption.In (Barbu et al., 2016;Mounce et al., 2020) the authors use fuzzy-based algorithms for SN control.The previously mentioned algorithms were tested by numerical simulation or/and applied in practical situations.
Summing up, the literature presents various control methods for SNs that have at their base a range of models.The use of complex models leads to an increased computational effort, limiting the capabilities of their use for realtime control.On the other hand, simplified models are successfully used for real-time control, with the drawback that the ability to assess pollutant concentrations in the overflowed wastewater is affected.Thus, the wastewater quantity could be better minimized, but without having any control on its quality.The control of the discharged wastewater quality could be regained by measuring or estimating pollutant concentrations in various key points of the SN.Advancements in research facilitated the use of more affordable sensors for online monitoring of such pollutant concentrations (Capodaglio et al., 2016) and provided methods for estimating some of them based on other compound concentrations (e.g., (Viviano et al., 2017) used caffeine and turbidity to estimate the phosphorus load in the SN overflowed wastewater).
The current paper aims to increase the operational efficiency of the systems that collect wastewater from the served collecting areas, such that the resulted wastewater flux (the effluent of the SN) to be efficiently directed to the WWT plant.The considered control algorithm is Model Predictive Control (MPC), combined with PSO.The main contribution of the current paper is the development and use of a simplified flow model in a MPC strategy for minimizing the quantity of overflowed wastewater, while at the same time, trying to keep a good quality of it by weighting each tank's overflow quantity with a quality index that is determined offline, as presented in Figure 1.In this way, there is no need to invest in equipping the SN with sensors to measure pollutant concentrations, these measurements being made offline by human operators on regular time intervals.
The paper contains four sections: the second section presents the structure, the influent, the mathematical model of SN and a brief description of MPC and PSO algorithms.The third section shows the results obtained with 10.1029/2023WR035448 3 of 14 SN control algorithm, alongside their analysis.The last section is dedicated to the conclusions and consequent research.

The Sewer Network Structure
A mathematical model of a SN for wastewater collection from a medium-sized city with a population of approximately 250,000 inhabitants is used.Figure 2 shows the structure of the SN.The network is sized for a city like Galați, Romania, considering that it combines three types of wastewater: domestic, pluviometric and industrial.The results were obtained by numerical simulation.
The network is composed of 7 storage tanks, noted as TKi, with  = 1..7 .Each one of the storage tanks has an influent input (Q in,TKi ), an overflow output (Q ovf,TKi ) and an effluent output whose flow can be controlled by applying a normalized command U TKi ∈ [0, 1] to the output valve or pump.As it can be seen in Figure 2 the outflows of five of the tanks are controlled by valves, while the outflows of the other two tanks (TK4 and TK6) are controlled by pumps.
The SN serves five collecting areas described by a surface and a population equivalent noted as S i , respectively pe i in Figure 2. The fifth collecting area serves the industrial part of the city.The influents of the tanks TK1, TK2, TK3, TK4 and TK5 are represented by the wastewater collected from each of the collecting areas, while the influent of the tank TK6 cumulates the effluents of TK3, TK4 and TK5.The influent of the tank TK7 is composed of the sum of the effluents of TK1, TK2 and TK6.The effluent of tank TK7 represents the effluent of the whole SN, which is in fact the influent of the WWTP.
The storage tanks characteristics described in Figure 2 (volume and maximum tank outflow) have been computed using the population equivalent of each collecting area as follows: 1.For the volume of each tank: a) If the tank directly collects water from a collecting area, a value of 0.05 m 3 for each population equivalent has been considered.For the tank that serves the fifth collecting area (the one containing the industrial part of the city) an additional capacity equal to the average industrial wastewater production per day (Q ind = 2500 m 3 /day) has been added.b) If the tank collects water from other storage tanks, a value of 0.05 m 3 for each average population equivalent of the upstream tanks has been considered.2. For the maximum tank outflow: a) If the outflow of the tank is controlled by a pump, the maximum flow of the pump is the maximum tank outflow.b) If the outflow of the tank is controlled by a valve, the maximum flow of it is computed to be 3 times higher than the average domestic wastewater production (Q pe = 0.15 m 3 /(day•inhabitant)) of the population living in the area served by the tank.

The Influent of the Sewer Network
The influent of the SN is composed of an influent defined for each one of the collecting areas.Each collecting area's influent has two components: a domestic component and a pluviometric component.Besides this, the fifth collecting area has an additional industrial component.

The Forecasted Influent
As the predictive control of the SN requires a prediction of the influent for a certain time span in the future, a forecasted influent has been generated, following the next procedures for each one of the components: 1.The domestic component has been generated as described in (Saagi et al., 2016): for each collecting area the domestic component flow/loads are obtained by multiplying the population equivalent of that area with the average flow/loads per inhabitant.Further, these values are adjusted with: a) A daily profile-flow/loads are adjusted to daily human activities (increased values during morning and evening, decreased values during afternoon and night); b) A weekly profile-flow/loads are adjusted to weekly human activities (decreased values during weekend); c) A yearly profile-flow/loads are adjusted to acquire the reduced values during the annual holidays.2. The pluviometric component is described only by flow and is generated based on the weather forecast, while accounting for each of the collecting area's surface.For a given value Q event measured in m 3 /day representing the forecasted pluvial flow over the entire city, the flow of the pluviometric component of the ith collecting area which has a surface equal to S i is described by the following Equation 1: 3. The industrial component is considered to have an average daily flow Q ind = 2500 m 3 /day and loads that are specific to a medium size brewery (Jern, 2006).Further, the flow and loads of this component are adjusted according to (Saagi et al., 2016): a) A weekly profile-an increase of the flow and loads is considered on Friday afternoon (weekly cleaning), followed by a decrease during the weekend; b) A yearly profile-flow/loads are adjusted to acquire the reduced values during the holidays.
Following these steps, a forecasted influent scenario has been generated considering an horizon of 28 days to be used for testing the proposed control strategy.The forecasted influent scenario accounts two pluviometric events: a rain and a storm event.

The Disturbed Influent
The forecasted influent is not 100% accurate, deviations being possible on the forecast of each component.For this reason, a disturbed influent scenario has been generated starting from the forecasted influent by altering: 1. the domestic influent flow of each hour with a Gaussian distributed random number of zero mean and a variance equal to  ( 0.2 ⋅ pe  ⋅ pe ) 2 where pe i is the population equivalent of the ith collecting area; 2. the industrial influent flow of each hour with a Gaussian distributed random number of zero mean and a variance equal to  (0.2 ⋅ ind) 2 .

The Mathematical Model of the Sewer Network
The mathematical model was previously developed in (Vasiliev et al., 2022d), being implemented in BSMsewer Software Package (Saagi et al., 2016).The same model has been used in (Vasiliev et al., 2022c) for developing a series of control strategies aiming to optimize the SN operation.As noticed in (Vasiliev et al., 2022c), the only tanks that are overflowing when no control is considered (all controls are set at 100%) are TK 2 ,TK 3 and TK 6 .
Analyzing the network structure, it can be concluded that changing the controls of TK 1 ,TK 2 ,TK 6 and TK 7 will not lead to any improvements of the environmental impact and therefore keeping the control action values of these tanks at 100% is the best solution.Consequently, only the output controls of tanks TK 3 ,TK 4 and TK 5 will be considered when designing a controller that would be able to minimize the quantity of wastewater overflowing during heavy pluviometric events.
As also observed in (Vasiliev et al., 2022c), numerical integration of this model over a time horizon of 28 days using a mid-range computer (Intel® Core‫‬ ™ i5-6200U CPU @ 2.30 GHz with 8 GB of RAM) takes approximatively 30 s, which can be a problem when implementing the MPC controller because the optimization over all prediction horizon must be completed within one step of the MPC algorithm.For this reason, a simplified model is necessary to decrease the time required for its numerical integration by keeping only the flow part of the model.
As the model that needs to be simplified is implemented in BSMsewer (Saagi et al., 2016), the imposed changes to be made will consider its structure.In BSMsewer, the mathematical model of SN contains three components: the transport component, the flushing component and the storage component, each one considering the wastewater flow and its loads (soluble and particulate COD, ammonia and phosphate concentrations).Thus, the simplified model is obtained by removing from BSMsewer model all states and dynamics related to the pollutant storage and transport.The resulted simplified model decreased the integration time from 30 s to approximately 3 s (for the same time horizon of 28 days).
The implementation of the simplified model in Simulink considered the possibility to run the model multiple times in parallel.Therefore, all configuration variables, all input variables, all control inputs and all output variables have been packed in data structures that could be separately passed to each of the simulations.All the initial states can be set to any desired value and all the input vectors containing data varying in time have been augmented with specific timestamps making possible to integrate the model for any specific time interval [t x ,t y ].
The possibility to run the model multiple times in parallel has a series of advantages.One of them is that the metaheuristic optimization could be solved in a reduced time, as it requires integrating the model many times, which could be made in parallel.However, the decision to use parallel computation or not should be taken after analyzing some factors regarding the used computational machine and the MATLAB Parallel Computing Toolbox.This is because the MATLAB Parallel Computing Toolbox needs some extra time to load and clean the parallel pool that must be compensated by the parallel computing time reduction.Therefore, the computer running the optimization procedure should have enough CPU cores for parallel computing, so the reduction of the optimization time is higher than the delays added by MATLAB Parallel Computing Toolbox.
Such an analysis will be presented, after implementing the MPC algorithm, in Section 3.2 devoted to the analysis of the computational effort of the algorithm.

Performance Indicators
BSMsewer defines a series of 10 performance indicators (Saagi et al., 2016) to be used for the performance assessment of the control strategies applied to SNs.These indicators can be computed individually for each tank, but also globally for the whole SN.Table 1 presents a short description of them.
As the simplified model accounts only the wastewater flow through the network, only three of the 10 performance indicators defined for BSMsewer can be calculated for each storage tank of the SN, but also globally for the whole network: • N ovf [events/year]-yearly overflow frequency; • T ovf [days/year]-yearly overflow duration; • V ovf [m 3 /year]-yearly overflow volume.
Consequently, using the simplified model, only the overflow quantity can be evaluated.For this reason, to also assess the overflow quality when implementing optimal control strategies that make use of this simplified model, the following methodology is proposed: using a global performance indicator obtained by summing the overflow volume of each of the tanks weighted by a factor determined by medium loads expected for the evaluated time interval in each of the tanks.The medium loads will be measured offline by the operators of the SN on a regular schedule and, when required, the weights will be updated consequently.

Model Predictive Control Formulation
The objective of controlling the SN is to minimize its environmental impact, mostly caused by overflow events during heavy pluviometric events, by finding the control strategy that minimizes the volume of overflow.
MPC (Camacho & Bordons, 2007) is an advanced process control method used for controlling complex systems, while, at the same time, satisfying a set of constraints.MPC requires an iterative optimization procedure of a system model over a finite horizon (prediction horizon-T H ) which keeps being shifted forward with a certain step (T S ) such as .Only the first step of the computed optimal controls vector is implemented within the real system.After this, the new state of the system is sampled and the procedure repeated.The MPC algorithm is illustrated in Figure 3a.
It needs three essential elements: (a) a dynamic model of the system, (b) a definition of the fitness function over the chosen prediction horizon, and (c) an optimization algorithm to compute the optimal control actions in order to minimize the fitness function without violating any of the defined constraints (such as low and high limits).
For this purpose, the previously presented simplified model of the SN will be used as a dynamic model for the MPC algorithm, while the fitness function is defined as the global yearly overflow volume (Equation 2): (2) Performance Indicators Defined by Benchmark Simulation Model (BSM)sewer (Saagi et al., 2016)

Optimization Using Particle Swarm Optimization Method
For the computation of the optimal controls that minimize the fitness function (   * = argmin   ) PSO algorithm (Bonyadi & Michalewicz, 2017) is used.It is a metaheuristic optimization algorithm that defines a solution of the optimization problem as a particle that is moving through the solution space.At every optimization step, each of the particle's velocities in the solution space is updated with a cognitive term and a social one (Equation 3).These velocities are then used to compute the new particle's positions (Equation 4).
Equations 3 and 4 x represents the position of the particle,  is the velocity of the particle,  is the best-known position of the current particle, while  is the best solution in the particle neighborhood (the particle neighborhood is randomly chosen at each iteration and consists of 25% of the swarm's population).The dynamics of the particles are adjusted using the following parameters: -the particle adaptive inertia (a value dynamically increasing from 0.1 to 1.1), c 1 -the self-adjustment weight (chosen to be c 1 = 1.49), c 2 -the social adjustment weight (chosen to be c 2 = 1.49).Besides this, two uniformly distributed random values between 0 and 1 are involved in the computation of the particle's velocity: r 1 and r 2 .
As stated before, only the output controls of the tanks TK 3 ,TK 4 and TK 5 will be considered when designing the controller, so, in consequence, the particle position in the solution space will have 3•n T dimensions, corresponding to the control actions of the 3 tanks over the n T steps in the prediction horizon.Thus, a particle position (a solution of the optimization problem) will be described by: The PSO algorithm has the following characteristics: • swarm size (number of particles in the population): 20; • initial population: • one particle encoding the case when all the controls are set to 100%; • one particle defined as the solution found at the previous MPC step; • 18 randomly generated particles; • stop conditions: • reaching the maximum number of iterations (100); • reaching 30 consecutive iterations with no improvement; • fitness value reaching 0.

Model Predictive Control Parameters
The two parameters of the algorithm, the prediction horizon time T H and the step size T S , should be chosen in such a way as to ensure the best performance of the controller.When choosing these parameters, a series of factors need to be considered: 1. Optimization time should be lower than T S .Considering the PSO algorithm characteristics, the maximum number of fitness evaluations at each step of MPC is the swarm size multiplied by the maximum number of iterations which means 2,000 evaluations.Also, because each fitness evaluation requires numerical integration of the model, at each step of MPC, the model will be integrated 2,000 times.Consequently, T S should be chosen such as the optimization time falls within it.As the simplified model integration takes about 3 s for the entire time horizon of 28 days and during optimization the model will be integrated only over the prediction horizon [t,t + T H ] which is smaller than 28 days, a minimum step size of 1 hr should be enough to fulfill the optimization time constraint.2. T S should be low enough to have a good resolution of the solution space.3. T H must be multiple of step size T S .4. T H must be low enough to represent a reliable weather prediction interval.5. T H must be high enough to have a larger solution space.6.There must be a compromise between T S and T H such as the optimization algorithm keeps a reduced number of dimensions of the solution ensuring the algorithm convergence to the global minimum.

Simulation Tests
The simulation tests were performed using the simplified model as the model required by MPC and the model implemented in BSMsewer as the process to be controlled, as presented in Figure 3b.For the influent of the simplified model, the forecasted one (see Section 2.2.1) has been considered, while, for the process, the disturbed one (see Section 2.2.2) has been used.
Three cases have been considered:

Case No. 1
The first set of tests considers T S = 1 hr, while T H takes the following values 3, 4, 6, 8, 12 and 24 hr.Each of the 6 tests has been run 20 times.For each run, the reduction of the overflow volume, V ovf , and of the overflow quality index, OQI, comparing to the "no-control" performance has been computed and analyzed (minimum value-MIN, maximum-MAX, average reduction value-AVG and the standard deviation value-STD), resulting the data presented in the first 6 lines of Table 2.
As it can be observed different values of the T H do not affect the average reduction of V ovf .However, the standard deviation gets lower when the prediction horizon is smaller.This could be explained by the fact that increasing the value of T H leads to an increase in the number of dimensions of the solution space and, therefore, the chance to get stuck in local minima gets higher.Analyzing the overflow quality index, OQI, an average reduction of about 40% can be observed, but with a higher standard deviation, the reduction ranging from 12% to 51%.The obtained results show that the time T H can be kept as low as possible in order not to increase the integration time of the simplified model.
For the next set of tests, the step size was increased to T S = 2 hr, while keeping the number of steps in the prediction horizon the same as in the case of the previous tests, thus, using the following values for T H : 6, 8, 12, 16, 24 and 48 hr.Each of the 6 tests has been run 20 times.The statistics of the results are presented in the next 6 lines of Table 2.
Comparing the results obtained with both test sets, it can be noticed that doubling the step size did not provide any significant improvement of neither the indicators.Besides this, a higher value of T S leads to an increase in the standard deviations of the average V ovf reduction, which means that the obtained solutions are not so well grouped around the average.
The highest average percentual reduction of V ovf has been obtained for T S = 1 hr and T H = 6 hr (50.63%), however, the differences between this and the previous test are not significant.On the other hand, the highest average percentual reduction of OQI has been obtained for T S = 1 hr and T H = 12 hr (44.12%).
Figure 4 shows an example of optimal commands evolution during a pluviometric event and presents the TK 3 , TK 4 , TK 5 and TK 6 filling degrees when the optimal commands are used compared with the case when no control is used.The green line shows the start of the pluviometric event, while the red markings show the overflow event that occurs when no control is considered.As noticed, using the optimal commands prevents this event, by using the free space in the other three tanks.

Case No. 2
When analyzing results obtained in the previous case, it could be observed that the standard deviation is high for all cases, OQI reduction ranging from 12% to 52%, meaning that the solution space has many local minimums with similar reduction of the volume of overflow but different reduction of the quality index.
For this reason, it is necessary to discriminate between the local minima with similar reduction of V ovf but different reductions of the OQI.Integrating OQI in the simplified model will increase the complexity of the model and consequently, the computational effort of MPC, but also it will require sensors for measuring various loads in the SN which will require additional costs.To avoid this, a second solution is proposed: a preliminary analysis of the medium loads of each of the storage tanks, and based on this analysis, setting up weights for each of the tank's overflows (w TKi ).This way, the fitness function will become Equation 6: For setting up the weights w TKi , the medium loads of the inflow for each of the storage tanks when the forecasted influent is used (without the weather component) have been computed and converted to Activated Sludge Model (ASM)2d variables (Gernaey et al., 2011).The ASM2d variables were used to calculate a quality index for each of the tanks in the same way the effluent quality index is calculated for BSM (Alex et al., 2008) WWTPs.The obtained indexes measured in kg-polluting units/day were further divided by the medium flow resulting quality indexes measured in kg-polluting units/m 3 : • QI TK1 = 5.2884 kg-polluting units/m 3 • QI TK3 = 5.2884 kg-polluting units/m 3 • QI TK4 = 5.2884 kg-polluting units/m 3 • QI TK5 = 6.1748 kg-polluting units/m 3 • QI TK6 = 6.1074 kg-polluting units/m 3 • QI TK7 = 6.3009 kg-polluting units/m 3 The quality indexes for the first four tanks have the same value as only the forecasted domestic influents of the corresponding collecting areas have been used in their computation.As it was presented in Section 2.2.1, the forecasted domestic influent corresponding to an area considers the flow/loads being obtained by multiplying the population equivalent of that area with the average flow/loads per inhabitant.Thus, the pollutant concentrations profile remains the same no matter the area's population equivalent, resulting in the same quality index.
The computed quality indexes were scaled to  [1 max] range resulting in the weights  TKi , according to Equation 7: Considering T S = 1 hr and T H = 6 hr (n T = 6), 3 tests were carried out for 3 values of w max .Each of the 3 tests has been run 20 times, obtaining the statistical results presented in the rows related to case no. 2 in Table 2.
As observed, weighting the overflow volume of each tank in the fitness function leads to a small improvement of the average reduction of V ovf and an improvement of 10% of OQI reduction compared to case no. 1, proving First column shows the tanks filling evolution when no control is considered, the middle column shows the tanks filling evolution when the proposed control strategy is used, while the third column shows the control actions computed by the optimization algorithm.
that the new fitness function provides better results in terms of the environmental impact.The average reduction of the volume of overflow decreases as w max increases (thus, the best V ovf reduction occurs for w max = 2), while the best average reduction of OQI (50.39%) occurs for w max = 3.In results it can be observed that the differences between w max = 2 and w max = 3 are not significant.Consequently, w max = 2 can be chosen if the main objective is V ovf reduction (e.g., to avoid floods).If the main objective is reducing the environmental impact through pollutant loads, then choosing w max = 3 is recommended.

Case No. 3
Further improvements could be made by leaving, at the end of the current prediction horizon, the SN storage tanks in a state that could be more advantageous for the performance of the next MPC step.This can be achieved by trying to keep the wastewater level in the storage tanks as low as possible at the end of the prediction horizon.
For this reason, a penalty term (P) has been added in the fitness function, weighted by a value w p , according to Equation 8: As only the outputs of three tanks are subject to optimization, two versions for computing the penalty have been implemented: one that uses the levels for all the tanks Equation 9 whereas the other one uses only the levels in the three tanks being subject to optimization Equation 10.
For each of the penalty versions, two tests were carried out for two different values of w p resulting 4 tests.Each of them has been run 20 times.The statistics of the results can be seen in the last 4 rows of Table 2.
As observed, integrating a penalty in the fitness function does not lead to any further improvements of neither the volume of overflow, nor the overflow quality index, but, on the contrary, the computational effort of the optimization algorithm increases (as the fitness function will never be able to reach 0, therefore one of the PSO stop conditions will never get reached).

Computational Effort
One of the most important aspects when implementing MPC algorithms is the integration time of the prediction model.It is imperative to solve the optimization problem within one step of MPC.In addition, the shorter the integration time, the greater the flexibility of the parameterization of MPC algorithm, which could lead to improved performance.
Consequently, three directions could be followed for reducing the computational effort: 1. Using a simplified model with a shorter integration time.
2. Using parallel computing for integrating the model.3. Using a more powerful computational machine for integrating the model.
The first direction was discussed in Section 2.3 where the SN model was simplified in order to decrease the computational effort of its integration.By simplifying the model, the integration time decreased from 30 s to approximately 3 s for a time horizon of 28 days.Of course, acknowledging the slow process dynamics, using the full model in the proposed MPC algorithm on a powerful machine could still allow for solving the optimization problem within one step of MPC.However, when using the full model, at each step of MPC, the full state of the SN, including pollutant loads, should be sampled, which will require sensors in various points of the SN, leading to additional costs.Moreover, considering that the simplified model implements exactly the same hydrological equations as the full model, and that the fitness function relies on the volume of overflow and some weights computed offline, using the full model will lead to the same control system performances.
Thus, for the proposed strategy, using either the full model or the simplified one will provide the same performances.Both of them will allow further improvements of MPC performances: • using the simplified model will allow to achieve a better resolution of the solution space (by decreasing the step size T_S) and to improve PSO algorithm performances by increasing the population size or the maximum number of iterations.• using the full model will allow to adapt the fitness function described by Equation 6 such as the weights for each of the tank's overflows will be calculated based on predicted loads.
The other two directions for reducing the computational effort will be discussed in what follows.
The implementation of the simplified model in Simulink allows the parallel running of the model aiming at reducing the time of the metaheuristic optimization algorithm as it requires integrating the model multiple times, which could be made in parallel.However, the MATLAB Parallel Computing Toolbox takes some extra time to load and clean the parallel pool, time that must be compensated by time reduction due to the parallel computing.Consequently, the computational machine performing the optimization should have enough CPU cores, so that the reduction of the optimization time is higher than the delays added by the MATLAB Parallel Computing Toolbox.
For this reason and to analyze how much the integration time decreases when a more powerful machine is used, tests were carried away on 4 machines with different specifications.Before starting the optimization algorithm, the seed of the pseudorandom number generator has been set to a constant value, the same for all the tests, so they converge to the same solution.Tests were performed by running the optimization algorithm for only one prediction horizon (T S = 1 hr, T H = 6 hr) on both cases: with and without parallel computing on each machine and by measuring the optimization times.
To make sure that all eight cases have the same number of fitness evaluations, for the PSO algorithm described in Section 2.4.2,only one stop condition has been kept: maximum number of iterations (100) reached.Consequently, each test requires an exact number of 2,000 fitness evaluations, that means 2,000 integrations of the simplified model.
Table 3 presents the results of the computational effort tests: As observed, neither of the computational machines achieved a reduced optimization time when using parallel computing because of the times required for loading and cleaning the parallel pool.However, when no parallel computing was used, the optimization time depends on the machine's performance, the obtained times ranging from 32 to 6 min which are more than enough for an MPC step.Using a more performant computer leads to a significant decrease in the optimization time (that is, the most performant Personal Computer obtained a total optimization time of 6 min, which means that 3.479242 s were required for each iteration of the PSO algorithm and 173 milliseconds were required for a fitness function evaluation.
Obviously, if a performant computer is used for MPC, because of the significant decrease in the optimization time, further improvements of MPC performances can be achieved in future work by: • decreasing the step size T S -achieving a better resolution of the solution space.
• improving PSO algorithm performances by increasing the population size or the maximum number of iterations.

Conclusions
A method for minimizing the environmental impact of SN overflows is proposed, based on MPC algorithm, and using PSO as the optimization method.A simplified model of the SN has been developed for reducing the computational effort of the MPC algorithm, allowing a smaller optimization step.Concluding the results, the case that minimized the globally V ovf achieved good reduction (≈50%) in terms of overflow quantity, but pollutant quantity minimization is not ensured.This has been improved in the second case, were offline computed tank weights depending on the medium loads were integrated in the fitness function.This case achieved ≈50% reduction of both V ovf and OQI.Notably, the methodology to compute tank weights based on offline measurements of pollutant loads could readily be adapted for online measurements, enhancing, this way, the efficiency of the proposed optimal control algorithm.A third case that integrates the state of the system at the end of the prediction horizon in the fitness function as a penalty did not lead to any further improvements.

Figure 1 .
Figure 1.The proposed control strategy.

Figure 2 .
Figure 2. Schematic representation of the sewer network structure.

Figure 3 .
Figure 3.The control structure: (a) the proposed model predictive control (MPC) structure; (b) the MPC structure used in simulation.

Figure 4 .
Figure 4. Example of optimal command evolution during a pluviometric event.Green line marks the beginning of the rain event.First column shows the tanks filling evolution when no control is considered, the middle column shows the tanks filling evolution when the proposed control strategy is used, while the third column shows the control actions computed by the optimization algorithm.

Table 1
(Alex et al., 2008)11)to Activated Sludge Model (ASM) state variables using the method proposed by(Gernaey et al., 2011).The ASM state variables are further used to calculate the index in a similar way as the effluent quality index of BSM WWTPs gets calculated(Alex et al., 2008)

Table 2
Results Obtained in Simulation in the 3 Cases Considered • QI TK2 = 5.2884 kg-polluting units/m 3 Case no.J T S [h] T H[h]