Continuing with the algorithm, the leak magnitudes are computed with Equations (31) and (32). These results are shown in Figure 11.

**Figure 11.** Estimation of leak magnitudes (2nd experiment).

#### **5. Conclusions**

The simultaneous leak detection and isolation problem is currently an open and challenging problem. Even though extensive research in the field is currently in progress, to the best of our knowledge, only simulation results have been reported until now. One of the reasons is that the distinguishability of two (or more) simultaneous leaks depends on the input.

In this work, an *LDI* methodology for the two simultaneous leak detection and isolation has been proposed based on an algebraic observer that uses the injection, the inputs and outputs of the system, and the corresponding time derivatives, by applying an appropriate input. The time derivatives are computed using *Numerical Differentiation with Annihilators* and the approach has been successfully applied to real data. This methodology could be extended to more general cases of simultaneous leaks in real operational conditions, such as the case of the *SIAPA* aqueduct in Guadalajara, Mexico.

**Author Contributions:** A.N.-D .: software, validation, data curation, writing—original draft preparation, methodology, formal analysis, and investigation; J.A.D.-A.: conceptualization, methodology, formal analysis, writing—review and editing, and investigation; O.B.: writing—review; G.B.: writing review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Tecnológico de Monterrey through the School of Engineering and Science.

**Acknowledgments:** Authors would like to thank Tecnológico de Monterrey for the financial support and the facilities granted for fulfillment of this research article. All databases were obtained by the second author during his PhD studies at Cinvestav Guadalajara.

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

#### **References**

## *Article* **Chlorine Concentration Modelling and Supervision in Water Distribution Systems**

**Ramon Pérez 1,2, Albert Martínez-Torrents 2, Manuel Martínez 2,3, Sergi Grau 4,5, Laura Vinardell 2, Ricard Tomàs 4, Xavier Martínez-Lladó 2 and Irene Jubany 2,\***

**Abstract:** The quality of the drinking water distributed through the networks has become the main concern of most operators. This work focuses on one of the most important variables of the drinking water distribution networks (WDN) that use disinfection, chlorine. This powerful disinfectant must be dosed carefully in order to reduce disinfection byproducts (DBPs). The literature demonstrates researchers' interest in modelling chlorine decay and using several different approaches. Nevertheless, the full-scale application of these models is far from being a reality in the supervision of water distribution networks. This paper combines the use of validated chlorine prediction models with an intensive study of a large amount of data and its influence on the model's parameters. These parameters are estimated and validated using data coming from the Supervisory Control and Data Acquisition (SCADA) software, a full-scale water distribution system, and using off-line analytics. The result is a powerful methodology for calibrating a chlorine decay model on-line which coherently evolves over time along with the significant variables that influence it.

**Keywords:** chlorine; water distribution networks; modelling; supervision; decay model

#### **1. Introduction**

Disinfection is one of the most important steps in water treatment, as it must ensure the microbiological safety of the water generated, not only after treatment, but also throughout the transport process to the consumption point. Many countries use chlorine-based chemicals (sodium hypochlorite, chlorine dioxide, chloramines, etc.) to achieve this objective, as they guarantee the degree of residual disinfection potential that is required by their laws [1]. If required, booster disinfection stations are installed at different points in the network. Their need and best location can be optimized using models and tools based on estimates of the chlorine concentration. Chlorine concentration is precisely one of the most relevant parameters to consider for the water distribution network (WDN) quality management. Although chlorine ensures the absence of pathogens, it is the main cause of the formation of disinfection byproducts (DBPs) [2]. Most of these compounds are toxic or carcinogenic for human health and need to be controlled to ensure drinking water safety [3]. Thus, European legislation limits the concentration of some DBPs in drinking water [4].

Nowadays, given the lack of reliable and applicable models for predicting chlorine behavior, disinfection management is not optimal in most WDNs, since it is based on

**Citation:** Pérez, R.; Martínez-

Torrents, A.; Martínez, M.; Grau, S.; Vinardell, L.; Tomàs, R.; Martínez-Lladó, X.; Jubany, I. Chlorine Concentration Modelling and Supervision in Water Distribution Systems. *Sensors* **2022**, *22*, 5578. https://doi.org/10.3390/ s22155578

Academic Editors: Eduard Llobet and Angela Maria Stortini

Received: 9 May 2022 Accepted: 19 July 2022 Published: 26 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a point-specific control as opposed to the consideration of the whole network [5]. There are often dark points where the chlorine level may be too low, along with over-chlorination at other points (particularly in summer), with a subsequent increase in both operating costs and DBP concentrations.

The absence of robust models for predicting chlorine behavior in WDNs is fundamentally due to two aspects: (1) the complexity of modelling the hydraulics of the WDNs and (2) the need for on-line quality data. Although authors report good results in chlorine prediction in full-scale networks in some studies [6], the predictions become less accurate when the environmental conditions or the water composition change from those of the calibration. Such a situation is very common in WDNs fed with treated surface water.

Regarding the first aspect, WDNs are highly meshed and complex systems, the behavior of which is difficult to predict. The introduction of flow, level and pressure sensors, and automated metering readers (AMR) for consumption has recently increased the model accuracy [7,8]. Thus, the intense use of a large amount of hydraulic data together with hydraulic models and numerical simulators allow prediction of residence time, which is one of the main parameters needed for successful water quality prediction.

Regarding the second aspect, it is mandatory to obtain information on water quality in the effluent of the drinking water treatment plant and the relevant points in the WDNs in order to predict the behavior of chlorine. Several studies [6,9] base their decay models on parameters that are easily measured on-line, such as temperature, pH, redox potential, conductivity, turbidity, and chlorine concentration. Nevertheless, the calibration and maintenance of these models for their on-line use is seldom performed.

Common models for chlorine modelling in WDNs are first-order linear differential equations [10] such as:

$$\frac{d\mathbb{C}}{dt} = K\_b \cdot \mathbb{C}^a \tag{1}$$

where *Kb* is a constant that contains the different parameters and physicochemical phenomena that may affect chlorine decay, such as natural organic matter, inorganic compounds, or temperature [6]. More sophisticated models have also been studied, including second-order models [11].

Furthermore, knowing the effect of the parameters influencing chlorine decay is also important, as this will, in turn, allow the prediction of such decay and, in some cases, the application of corrective measures to reduce its effect.

Equations of different complexity have been used to model the effect of some of these parameters. One of the most influencing and studied parameters is temperature, which is usually based on Arrhenius's model [12] or other power models. Liu et al. [13] took pH and temperature into account in their models, thus differentiating the effect of HOCl and OCl species on pH. Similarly, Arevalo, in his doctoral dissertation [14], used a model that considered temperature and UV254 as an indication of organic matter. In this case, two decay constants were used, one related to chlorine decay on bulk water and one related to chlorine decay on water close to pipe wall.

Hassan et al. [15] studied the specific case of organic matter adsorbed onto goethite, which is the predominant iron oxide in pipe deposits, in order to see how effectively the presence of organic matter increased the decay rate. Their main conclusions were: (i) an increase in temperature causes an increase in the decay constant and therefore in the decay rate, (ii) the pH has not been seen to greatly affect decay, (iii) a higher initial chlorine concentration leads to a lower decay rate, (iv) a higher organic matter concentration, in general dissolved organic matter (DOC), causes an increase in the rate of decay, (v) an increase in the velocity of the water flow through the pipe causes an increase in the rate of decay, and (vi) the concentration of ammonium, nitrites, iron, and manganese seems to affect the rate of decay, causing its increase.

Chlorine decay first-order equations can also be used in software like EPANET, which has sufficient power to simulate and predict the concentration of chlorine in the network. EPANET is a public domain software for WDN modelling developed by the United States

Environmental Protection Agency (US EPA). This software can perform transient simulations of hydraulic behavior and water quality in pressurized pipe networks. In order to properly model water quality and its time evolution at consumer points, it is mandatory to have a reliable hydraulic model of the WDN. However, like any simulation software, EPANET depends on the availability and application of continuous data into robust models. The default built-in models have fixed calibration parameters that are not easily extrapolated in most real cases. The approach some authors take to overcome these barriers is to modify and recalibrate the default models included in the EPANET database based on the real system data to be modelled [16].

Another important aspect when modelling chemical reactions in pipes is the different behavior in wall and in bulk water. In the literature, bulk reactions are usually considered first-order and wall reactions are considered zero-order. Values for the bulk reaction coefficient are usually obtained using laboratory measurements [17]. Nevertheless, changing water characteristics in the network requires updating the model. There are a few approaches for the quality model calibration. This calibration requires a validated hydraulic model and water quality data. This is often carried out in a well-monitored part of the network and then generalized to the whole network [18,19].

This paper focuses on the chlorine decay process and variables that affect it, the models for concentration prediction, and their application within WDNs in a specific case-study. First, an on-line calibration procedure, with available data from the transport network, is adapted to the full-scale system and performed over a long period so that the evolution of the decay parameters can be studied. This model is used to predict the chlorine concentration in the distribution network and validated with discrete monitoring data. The dependence of chlorine decay on the relevant variables is also studied. Finally, this dependence is compared with the evolution of the parameters estimated using the on-line calibration method. The aim is to illustrate how the intense use of models and available data can provide a better understanding of the behavior of chlorine in a WDN and, thus, be used to support decision-making to improve water quality.

#### **2. Materials and Methods**

#### *2.1. Case Study Network*

The case study in this work is a WDN in Catalonia (Spain) (see location in Figure 1) managed by Aigües de Manresa, who provided the network configuration information and hydraulic and water quality data for a period of 14 months (2017–2018). The water supplied comes from the Llobregat River and goes through a prechlorination step with sodium hypochlorite (Apliclor Water Solutions S.L., Sant Martí Sesgueioles, Spain) or chlorine dioxide generated using sodium chlorite (Apliclor Water Solutions S.L., Sant Martí Sesgueioles, Spain) and hydrochloric acid (Apliclor Water Solutions S.L., Sant Martí Sesgueioles, Spain) (depending on the season), a sand filtration process, and a final disinfection step with sodium hypochlorite.

Two parts of the network were used in this study: the transport network and the district metered area (DMA). The transport network (Figure 2) consists of two water storage tanks (T1 and T2) equipped with sensors for chlorine concentration (input of T1 (Cl1) and output of T2 (Cl2)), flowmeters (outflows from the tanks, Q1 and Q2), and water level (H1 and H2). Water flows from T1 to T2 through a 6859 m main. T2 is a boosting station with known sodium hypochlorite (Apliclor Water Solutions S.L., Sant Martí Sesgueioles, Spain) addition. The geometry of the tanks (volume) and pipes (length and diameter) are known. Water from T2 is distributed to the rest of the network (through Q2) of the DMA.

The DMA corresponds to a residential area. The hydraulic model includes 572 nodes and 610 pipes with a total length of 31 km, providing water to 300 consumers. Water flows by gravity. There are two quality-sampling points where the chlorine concentration is measured weekly. Figure 3 presents the model of this DMA visualized in EPANET. The input tank (corresponding to T2 in Figure 2) and the two sampling points are highlighted (S1 and S2).

**Figure 1.** Study site location, a network in Catalonia supplied by the Llobregat River.

**Figure 2.** Outline of the network section used for the on-line calibration.

**Figure 3.** DMA network model in EPANET.

#### *2.2. On-Line Calibration*

A very well parametrized system in terms of hydraulics and chlorine concentration (at least at two points) is required to calibrate the chlorine decay constant of a WDN. As in this study, the transport network often fulfils this condition. Therefore, the network used for on-line calibration in this study was the transport network shown in Figure 2. The objective was to find the decay constants for the model that best explained the chlorine concentration measured at the output of T2.

The chosen model was a first-order model. Higher-order models could be used with no fundamental changes in the methodology. Equation (2) shows that the chlorine concentration (*Cl2*) at the outflow of T2 depends on the input chlorine concentration (*Cl1*) and the residence time in the system (*t*). The solution of Equation (1) is as follows:

$$Cl\_2 = Cl\_1 \cdot e^{-K\_b \cdot t} \tag{2}$$

where *Kb* is the decay constant and *α* is considered as 1.

This dependence is defined by the decay constant *Kb*, which was calibrated on-line using the measurements available so that it was adapted throughout the year to the different water characteristics and environmental conditions. Estimations were performed on a weekly basis, since some information was only available at this frequency (chlorine dosing in T2).

The residence time (RT) in T1 was calculated from the hydraulic information available using (3). The weekly mean residence time at tank T1 and the pipe was calculated using the flowmeter data (Q1) and the volume of this subsystem.

$$
\overline{RT}\_1 = \frac{\overline{V}\_1 + V\_{pipe}}{\overline{Q}\_1} \tag{3}
$$

In order to estimate the mean water volume in T1, the level data of the tank (H1) and the geometric information was used. The residence time in T2 was calculated using the mean values of the volume obtained from the level data (H2) and the mean values of the tank effluent (Q2), as shown in (4):

$$
\overline{RT}\_2 = \frac{\overline{V}\_2}{\overline{Q}\_2} \tag{4}
$$

The chlorine concentration increase due to rechlorination (*Cladded*) was calculated using the added volume of chlorine divided by the mean volume of water treated:

$$
\overline{Cl}\_{addcd} = \frac{\overline{\Delta V\_{Cl}} \cdot 143}{\int \, Q\_2} \tag{5}
$$

where Δ*VCl* is the volume in liters of the concentrated chlorine added weekly to the network and 143 is the concentration of the added chlorine in g/L (value obtained from the conversion of the 15% NaClO to reactive chlorine, see Section S1 in the Supplementary Material).

Finally, a decay *Kb* constant was calculated which explained the chlorine concentration *Cl*2 at the outflow of T2 given the residence time calculated using (6)

$$
\overline{\mathbb{C}I}\_2 = \overline{\mathbb{C}I}\_1 \cdot e^{-\mathbb{K}\_b(\overline{\mathbb{R}T}\_1 + \mathbb{R}\overline{\mathbb{R}T}\_2)} + \overline{\mathbb{C}I}\_{addcd} \cdot e^{-\mathbb{K}\_b\overline{\mathbb{R}T}\_1} \tag{6}
$$

where *Cl*2 was considered equal to 0.6 ppm, which is the set point of the chlorine control system in the boosting station. The algorithm for *Kb* calibration is shown in the Supplementary Material (Section S2).

The first order decay model is the most used in the water industry. Its decay constant includes all the dependencies related to environmental and water characteristics. Thus, the continuous updating of this constant is the guarantee of its reliability. The limitation of this methodology is the information required. Hydraulic information, that allows the determination of the residence time, must be available. Multiple chlorine concentration measurements and the exact volume of added chlorine between these measurements are also mandatory data.

#### *2.3. Chlorine Decay Calibration and Validation*

The chlorine decay first-order model validation for the transport network was carried out using available the on-line data of the chlorine concentration in the output of T2 considering the residence time in this tank. The data used covered the period from February 2017 to April 2018. The chlorine decay model was also validated for its use in the distribution network in the section where chlorine concentration is monitored.

Applying the calibrated decay model directly to the distribution network produced poor results. This was expected, due to the difference between the transport network and the distribution network regarding pipe size, materials, age, etc. To adjust the model, the available data period was divided into two sets: one for training the new distribution quality model and the remaining data for validation. There were 35 samples available in S1 and 11 samples in S2. Thus, the first 21 samples in S1 were used for the training and the remaining ones for the validation. The algorithm used for this calibration is shown in the Supplementary Material (Section S3).

#### *2.4. Parametrised Chlorine Decay Model*

The decay constant, determined from the available on-line data, evolved clearly throughout 2017. The question arose if this could be due to the effect of the available variables such as temperature, the initial chlorine concentration, the cumulative precipitation, and the turbidity at the drinking water treatment plant or not.

This suggested the idea of analyzing the variables that influence chlorine decay in order to generate an empirical model based on the available independent variables. The availability of considerable data (temporal and spatial) implies dealing with large amounts of data, multiple variables, and experimental noise, which hinders the direct extraction of valuable information.

Principal component analysis (PCA) is a multivariate statistical technique that allows the description of the data according to the variance [20]. This method transforms data in noncorrelated new variables by linear transformation, decreasing the data dimensions. New data description is more condensed and can describe patterns that are hard to identify in multivariable datasets. PCA has already been used to determine the physical and chemical parameters influencing chlorine decay [21]. Therefore, for being a powerful, reliable, and globally accepted tool when dealing with big data, PCA was selected to extract the main trends, patterns, and correlations among the variables (dimensions) [11].

Based on the PCA results, the chlorine decay constant was modelled using the available variables and a potential multiparametric model (7).

$$K\_b = K \cdot parameter1^a \cdot parameter2^b \cdot parameter3^c \cdot \dots \tag{7}$$

Specifically, a power model (8) and an Arrhenius model (9) [12] where calibrated using experimental data (temperature in 2017) and *Kb* obtained from the on-line calibration using the least square error fitting method implemented in the "Solver" function in Excel.

$$K\_b = K\_{power} \cdot T^t \tag{8}$$

$$K\_{\Phi} = A \cdot e \exp\left(-E\_{\mathfrak{a}}/RT\right) \tag{9}$$

where, *Kpower*, *a*, and *A* are constants, *Ea* is the activation energy (Jmol−1), *R* is the universal gas constant, and *T* is the temperature. Finally, the parametrized chlorine decay model was compared with that obtained in the on-line calibration to assess its coherence throughout the year.

#### **3. Results and Discussion**

#### *3.1. On-Line Calibration*

The decay constants *Kb* obtained are presented in Figure 4. In the upper graphic, the weekly evolution between February 2017 and April 2018, can be observed. A different icon was used for the data of each trimester to clearly identify the season of the year. In the lower graphic obtained, *Kb* are grouped by month to observe how this parameter evolves throughout the year (some months include estimations of both years). It seems clear that there may be a seasonal variation related to temperature.

**Figure 4.** Up: Mean Kb calculated weekly indicating the season of the year (by trimester). Down: Data of mean Kb grouped by month.

#### *3.2. Chlorine Decay Validation*

The calibrated model was applied to the peak episodes observed at the output of T2 due to the rechlorination and mixing effect. The dataset used in this validation was not used for the estimation. The dataset for calibration consisted of the mean values corresponding to the stationary state. Figure 5 shows the chlorine concentration data and the model prediction. It can be observed how this high-frequency dynamic is adjusted with the model obtained with the mean values. For this prediction, Kb evolves weekly.

**Figure 5.** Chlorine decay model validation at the output of T2. Up: data from January 2017. Down: data from March 2017. Due to high frequency of sampling measurement, data appears like a thick line.

For the distribution network simulation, a relation between the transport *Kb*, estimated by on-line calibration, with the distribution *Kb* \* was obtained by adjusting the concentration in the training set of chlorine sampling. This relation was applied to the entire period and the predicted concentration was compared with the measurements for the validation set of samples. A total of 40 days were simulated. The decay constant for both the bulk and wall were fitted using the first 21 samples of the chlorine concentration in S1. These are the first samples of upper graphic in Figure 6.

**Figure 6.** Simulation results for the two sampling points (**S1** and **S2**) compared with experimental samplings.

The result was that both decay constants minimize the error when the original *Kb* obtained in the transport system was divided by 2, as if the calibrated effect was distributed in the two phenomena (K\* b,bulk = K\* b,wall = Kb/2). The results obtained are compared with the available experimental data in Figure 6. The mean absolute percentage error was 16% for S1 (including calibration and validation samples) while it was 17% using only validation samples. Therefore, not significantly different deviations were obtained for the calibration and the validation steps. Graphically, the fit may seem poor; however, the concentration is lower in S2 than in S1 both in prediction and measurements, and the mean values in both sampling points are coherent between the prediction and measurements. One aspect that may justify part of the mismatching is that the exact hour of the day of the manual measurements was not available and, therefore, each experimental data may not be in its exact position. This difficulty could be overcome with on-line chlorine sensors instead of manual analysis. Figure 7 presents the measured chlorine concentration at the source (T2) and the chlorine prediction in the two sampling points (S1 and S2). The chlorine concentration decreases with the residence time, since the concentration in S2 is lower than in S1, and both are lower than in T2.

Finally, Figure 8 shows the network nodes colored by their chlorine concentration: green, black, or red, depending on whether their concentration is too low, acceptable, or too high, respectively. In fact, no red points exist in this area and period. Such a representation is very useful in order for the network operator to make decisions. Nodes in this figure correspond to those in Figure 3 and are presented with north at the top.

**Figure 7.** Measured chlorine concentration at the source (T2) and prediction at the two sampling points (S1 and S2).

**Figure 8.** Distribution of network nodes with low concentration (<0.4 ppm, green), acceptable concentration (0.4 ppm > chlorine < 1 ppm, black), and excessive concentration (>1 ppm, red).

Results from the PCA applied to the decay constant (K\_Cl\_decay) determined in the distribution network and other data available (temperature as Tavg\_C, initial chlorine concentration as Initial\_Cl, cumulative precipitation as Cumulative\_Prec, and turbidity at the drinking water treatment plant) are shown in Figure 9.

**Figure 9.** Variables for the two main dimensions of the PCA.

In this case, dimension one was related to the temperature, and dimension two was related to the initial chlorine concentration. Therefore, the decay constant was closely related to the temperature. Figure 10 shows the variables that had the most influence on dimension five, which were again the temperature and the chlorine decay, demonstrating the clear strong relation of the temperature on the decay constant.

**Figure 10.** Contribution of variables to dimension 5. Variables with contributions below the dotted line are considered not significant for that dimension.

Thus, it was concluded that the variable that had a higher effect on the decay constant was temperature. It was observed that the variables turbidity, precipitation, and initial chlorine did not excessively improve the fit between the decay constants of the model and the decay constants obtained. Therefore, only the temperature was used, since the effort required to obtain values for the rest of the parameters did not compensate the improvement of the model adjustment.

Experimental data from 2017 and *Kb* obtained with the on-line calibration were used to calibrate the parameters of the two equations, an Arrhenius model (8) and a power model (9), to predict the temperature effect on the decay constant. The following Equations (10) and (11) show the results obtained.

$$K\_b = 5.477 \cdot 10^{-8} \cdot T^{1.524} \left(\text{s}^{-1}\right) \tag{10}$$

$$K\_b = 3950 \cdot \exp\left(-49873/RT\right) \left(\text{s}^{-1}\right) \tag{11}$$

The fit obtained using the Arrhenius model and the power model were similar, although the Arrhenius one was slightly better. The Arrhenius model also determines the activation energy (J/mol), which is the minimum energy that the system needs for the reaction to take place. The ratio Ea/R obtained in this study was 5999 K, which is in accordance with other authors. For example, Courtis et al. [22] estimated 5388 K and 6701 K for two different water distribution systems, Powell et al. [23] obtained a range between 7500 and 9600 K, and Hua et al. [12] obtained a range between 8203 and 8727 K, depending on the type of water. The variability of the values for this ratio suggests that this is a water-specific parameter that might depend significantly on the natural organic matter composition [24].

Figure 11 shows the fit of the two models to data from 2017, and Figure 12 their forecast for the first days of 2018. In these figures, K Chlorine is the chlorine decay constant determined previously, K power the constant determined following the potential model, and K Arrhenius the constant obtained from the Arrhenius model. As it can be seen, the Arrhenius model provides good predictions while using only the temperature as an input parameter.

**Figure 11.** Chlorine decay constant fitting to temperature dependent models using data from 2017.

**Figure 12.** Chlorine decay constant forecast for the first days of 2018.

#### **4. Conclusions**

The literature review shows that only the simplest models of chlorine decay are applied to water distribution networks where the hydraulic behavior is complex enough. Even so, these models are seldom used due to the lack of proper calibration. In this paper, the performance of a decay model was evaluated when the parameters were calibrated using state-of-the-art techniques.

The calibration was first carried out in the transport network, where the on-line data allowed an on-line calibration. The relation between the decay constant in the transport and distribution networks used analytical data and, thus, it could not be done on-line.

The prediction error in the validation data was 17% and quite similar to the error obtained for the training set (16%), which meant that there was no overfitting. The decay constant obtained changed during the year following the assumed dependence of the chlorine decay on the temperature. This result suggested the possibility of using available data for predicting this decay constant.

The principal component analysis determined that the temperature was the parameter with higher effect on the decay of chlorine. The chlorine decay constant was obtained using temperature as an independent variable. The obtained constants were compared with the data-driven model obtained in the on-line calibration, showing a high correlation. While the dominant dependence on the temperature is not a novelty, it is important to ensure this unique dependency, as it guarantees that other characteristics of the water source will not be relevant. This has been studied throughout one year, and the results obtained by both models are coherent.

This procedure could also be applied to other quality parameters, such as disinfection byproduct concentrations, which are currently under investigation by the authors. The variables analyzed for chlorine decay estimation are being studied for the trihalomethanes formation prediction.

The final aim of this study is to increase knowledge within the network in order to enable decision-making processes regarding chlorine dosing (quantity and frequency) both in the disinfection process and the boosting stations, in addition to identifying whether other control systems are required to ensure the continuous good quality of supplied water to the final user at a minimum cost.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/s22155578/s1, Table S1: Pseudocode to estimate the decay chlorine constant; Table S2: Pseudocode to estimate the chlorine decay constant in the distribution network.

**Author Contributions:** Conceptualization, R.P., R.T. and X.M.-L.; Data curation, R.P., A.M.-T., M.M. and S.G.; Formal analysis, R.P., A.M.-T. and M.M.; Funding acquisition, R.T., X.M.-L. and I.J.; Investigation, A.M.-T., M.M. and L.V.; Methodology, R.P., A.M.-T. and M.M.; Project administration, X.M.-L. and I.J.; Resources, S.G. and R.T.; Software, R.P., A.M.-T. and M.M.; Supervision, R.P., X.M.-L. and I.J.; Validation, R.P. and X.M.-L.; Visualization, R.P., A.M.-T. and L.V.; Writing—original draft, R.P., A.M.-T. and L.V.; Writing—review and editing, R.P., L.V. and I.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data is not publicly available because the owners want to keep track of its use, as it is sensitive information from a full-scale system (hydraulic models). Nevertheless, software for data treatment and chlorine decay model estimation and validation including nonsensitive data (model\_estimation.m and model\_validation.m) are available at https://cs2ac.upc.edu/en (18 July 2022) in Matlab format.

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

#### **References**

## *Article* **Economic Linear Parameter Varying Model Predictive Control of the Aeration System of a Wastewater Treatment Plant †**

**Fatiha Nejjari 1, Boutrous Khoury 1, Vicenç Puig 1,2,\*, Joseba Quevedo 1, Josep Pascual 1 and Sergi de Campos 3**

**Abstract:** This work proposes an economic model predictive control (EMPC) strategy in the linear parameter varying (LPV) framework for the control of dissolved oxygen concentrations in the aerated reactors of a wastewater treatment plant (WWTP). A reduced model of the complex nonlinear plant is represented in a quasi-linear parameter varying (qLPV) form to reduce computational burden, enabling the real-time operation. To facilitate the formulation of the time-varying parameters which are functions of system states, as well as for feedback control purposes, a moving horizon estimator (MHE) that uses the qLPV WWTP model is proposed. The control strategy is investigated and evaluated based on the ASM1 simulation benchmark for performance assessment. The obtained results applying the EMPC strategy for the control of the aeration system in the WWTP of Girona (Spain) show its effectiveness.

**Keywords:** economic model predictive control; linear parameter varying modelling; wastewater treatment process

#### **1. Introduction**

Biological wastewater treatment plants (WWTPs) are complex nonlinear systems with large variations in their flow rates and feed concentrations. These plants have to be operated continuously taking care of strict environmental regulations. Thus, the use of advanced control strategies becomes necessary to make them more efficient.

The most widely used biological wastewater treatment is the activated sludge process (ASP). In the ASP, microorganisms are mixed with wastewater. The pollutants of the wastewater constitute the nutrient of the microorganisms. As the organisms feed on the organic pollutants in the wastewater, the pollutants are converted to more organisms, biomass, and some by-products. Following an adequate amount of treatment time, the mixture of microorganisms and wastewater, the mixed liquor flows from the aeration tank to a clarifier or settler where the sludge is separated from the treated water. Some of the settled sludge is continuously recirculated from the clarifier to the aeration tank to ensure the maintenance of adequate amounts of microorganisms in this tank. The microorganisms are again mixed with incoming wastewater where they are reactivated to consume organic nutrients. There are five major groups of microorganisms generally found in the aeration basin of the activated sludge process: (i) aerobic bacteria responsible for removing the organic nutrients, (ii) protozoa to remove and digest dispersed bacteria and suspended particles, (iii) metazoa to dominate longer age systems and clarify effluent, (iv) filamentous bacteria or bulking sludge, which are present when operating conditions change, (v) algae and fungi, which are photosynthetic organisms that are present with pH changes and older sludge.

**Citation:** Nejjari, F.; Khoury, B.; Puig, V.; Quevedo, J.; Pascual, J.; de Campos, S. Economic Linear Parameter Varying Model Predictive Control of the Aeration System of a Wastewater Treatment Plant. *Sensors* **2022**, *22*, 6008. https://doi.org/ 10.3390/s22166008

Academic Editor: Assefa M. Melesse

Received: 10 May 2022 Accepted: 9 August 2022 Published: 11 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The majority of the culture is mixed and reused with inlet wastewater to keep high reaction rates and sludge age characteristics. In particular, nitrogen is eliminated as follows: First, ammonium is oxidized producing nitrate under nitrification in the aerobic step. The nitrate that is produced is then converted into nitrogen gas by means of denitrification in the anoxic step. Thus, the control of aeration is very important because a low amount of dissolved oxygen can cause the biomass death. On the other hand, an excess of dissolved oxygen could cause the sludge to settle insufficiently. Moreover, because 60% to 80% of the global energy consumption is due to aeration and the operating costs of a WWTP [1], an excessive aeration is not desirable regarding economic efficiency.

The models that are usually considered for characterizing the WWTP processes are the ones developed by the International Association on Water Quality (IAWQ) known as Activated Sludge Models (ASMs) [2].

In this paper, optimal economic operation of the aeration system is considered to improve the efficiency and reliability of an ASP with intermittent aeration, which is used for the removal of nitrogen from domestic wastewater. The objective of the control is to design an aeration strategy (air-on and air-off periods) which minimizes the energy dissipated by the aeration system, with adherence to the limits of the effluent requirements and the operating constraints. The implementation of optimal operation strategies is therefore interesting because WWTPs face the challenge of treating water properly albeit ensuring the minimization of operational costs. This has been the driving force for the active research in the development of advanced control techniques and hierarchical control schemes to improve the operation of the WWTPs, see for example [3,4].

Model predictive control (MPC) has been the most successful advanced control approach applied to control WWTPs. This is due to the fact that MPC controllers allow in a straightforward manner the different operational requirements, the multivariate nature of the control problem (that could even include delay) and directly handling constraints on the control inputs, system outputs and/or internal states [5]. It can also include disturbance prediction, allowing to anticipate the appropriate control actions (feedforward) to achieve optimal performance according to defined criteria in the cost function, which can include different quality criteria and operational costs. Adjusting the MPC control strategy is carried out by suitable manipulating prioritization of different objectives of the performance index that could also include the use of soft constraints. In this way, MPC has become an attractive control strategy for a considerable number of WWTP applications in the last few years. Some examples of MPC control of WWTP can be found in [6,7]. In [6], a benchmarking of different hierarchical control structures for WWTPs that combines static and dynamic real-time optimization (RTO) and nonlinear model predictive control (NMPC) is presented. In [8], a procedure to find the best controlled variables in an economic sense for the activated sludge process in a wastewater treatment plant, despite the large load disturbances, is introduced.

Classical MPC formulation considers pre-established set points, and the objective functions related with error and energy effort have quadratic forms [5]. However, the determination of optimal and reachable reference set points in real time is not an easy task because of the existence of disturbances, set-point changes, time-varying parameters and model uncertainties, among others. This constitutes one of the main limitations of classical MPC. To remedy this issue, real-time optimizers (RTO) or steady-state target optimizers (SSTO) are used to pre-compute the reference set- points at a supervisory layer in the control hierarchy. Then, these pre-computed set points are sent to the lower layer, where a classical MPC behaves as a regulatory controller, forcing the process to follow the desired set points. However, in spite of the use of an RTO, not reachable trajectories might be generated because of the appearance of unexpected disturbances or set points variations, among others. Moreover, there is a delay between the different layers, because the lower layer receives the reference set points determined from the upper layer before its execution. These problems can be avoided by using economic MPC (EMPC) that optimizes process performance directly (e.g., by means of economic objective functions), eliminating the need

of generating reachable reference set points [9]. The first results of the application of EMPC to DO concentration control in WWTPs have been presented in a previous work from the authors [10] considering the nonlinear model of the plant. However, this leads to a nonlinear optimization problem.

Alternatively, this paper proposes an EMPC strategy using the linear parametervarying (LPV) framework to optimize the effluent quality and minimize the operational cost of a WWTP under operating and physical constraints. The objective is to minimize the energy used by the aeration system with the control of the dissolved oxygen (DO) concentrations in the aerated reactors and maintain the effluent concentration under the required limits. The proposed approach is based on real-time dynamic optimization methods. Optimization in MPC with nonlinear models presents a non-convex problem which is computationally demanding, especially when dealing with large-scale plants with complex dynamics such as the WWTP. Thus, the LPV framework allows the embedding of these nonlinearities in scheduling variables, which are functions of system states (i.e., qLPV). This allows obtaining a pseudo-linear model which is linear in state space but nonlinear in the parameter space and deriving a less demanding convex MPC optimization problem, since convex quadratic optimization tools can be applied. The stability and recursive feasibility of MPC with LPV models has been studied (see [11] for a review of the recent results). The application of dynamic optimization methods requires a sufficiently accurate mathematical model describing the wastewater treatment process. The present work uses the Activated Sludge Model No. 2 (ASM2) [12]. To illustrate the proposed approach a WWTP located in Girona (Spain) is considered as a case study.

In Section 2, the WWTP is described and modeled using a reduced ASM2 model, which is then represented in a qLPV form. The proposed EMPC strategy is introduced and described in Section 3, while the proposed MHE approach is presented in Section 4. The results are presented in Section 5, with simulation scenarios obtained from the application of the EMPC strategy on the Girona WWTP. Finally, some conclusions are given in Section 6.

#### **2. WWTP Description and Modeling**

#### *2.1. WWTP Description*

The Girona WWTP is a biological treatment plant designed to treat the wastewater generated by 200,000 inhabitant equivalents with a medium daily inflow of 35,000 m3/d. The processes of the plant can be divided into two main treatment lines: water and sludge (see Figure 1). The water line is separated into three phases: pre-treatment, primary treatment and secondary treatment. The secondary treatment is designed to convert biodegradable, organic wastewater constituents and certain inorganic fractions into new cell mass and by-products. The plant uses an activated sludge system and has three lines composed of three main reactors that are divided into various compartments and three clarifiers. Each line is made of two anoxic reactors located at the beginning, three aerated tanks and an anoxic tank followed by an aerated one. With this configuration, the plant can nitrify and denitrify with great efficiency. The anoxic and aerobic tanks have volumes of 1335, 4554, 1929, and 1929 m3 for anoxic and 1929, 1276, and 1409 m3 for aerobic, respectively. Oxygen is supplied to aerated tanks by the aeration system, which delivers air to each of the aeration tanks. The wastewater and activated sludge are separated into three parallel secondary settlers. The volume of each secondary settler is approximately 5024 m3. The activated sludge is internally recirculated from the last aerobic zone to the anoxic tank (210% of influent waste). Additionally, the wastewater is recirculated from the secondary settlers to the anoxic tank (45 to 100% of influent waste).

Figure 2 shows a standard WWTP technological layout. The wastewater flow enters into the biological part after the mechanical treatment. The nutrient removal takes place in the activated sludge reactor through the biological treatment. The first zone in this treatment is anaerobic, where phosphorus is released. The mixed liquor internal recirculation originates from the anoxic zone. The denitrification occurs in the second zone. The activated sludge returned from the clarifiers bottom, and the internal recirculation from the aerobic zones end is directed toward the anoxic zone.

#### *2.2. WWTP Modeling*

The Benchmark Simulation Model (BSM1), developed within the framework of COST Actions 624 and 682 [2], has been adapted to represent the Girona WWTP (see Figure 1).

**Figure 1.** Girona Wastewater Treatment Plant.

The Activated Sludge Model No. 1 (ASM1) describes the biological phenomena that takes place in the biological reactors, and it is supposed that no biological reactions take place in the settlers. Due to the complexity of the nonlinear model describing the different complex processes in the plant, various reduced models have been proposed in the literature [13–15] to aid in the online implementation of certain modern control schemes (e.g., MPC), which would have otherwise presented ill-conditioned or stiff numerical problems due to slow and fast dynamic interactions. In [15,16], one can see some successful implementations using reduced WWTP models in various areas of control applied to WWTPs. The reduced model as suggested in [14], which primarily involves certain simplification criteria for a reduced order of the rigorous high dimensional WWTP model, has been adapted to conditions representing the Girona WWTP. This basically involves the derivation of the reactor model based on mass balances of the wastewater species, which are generally expressed as follows:

$$Accuracy = Inflow - Outflow + Reaction$$

Validation of the reduced model considering data from the ASM1 and the reduced model has been undertaken in [14]. In simplifying the complex model, a systematic reduction process of the high-dimensional model considers some assumptions, with the principal conditions given as follows:

Under these conditions, the resultant state variables of the reduced model are therefore the chemical oxygen demand (*XCOD*), the dissolved oxygen concentration, (*SO*), heterotrophic biomass, *XBH*, ammonia concentration (*SNH*), nitrate concentration (*SNO*) and autotrophic biomass (*XBA*). The control of oxygen concentration (*S*0) in the aerobic tanks is via the manipulation of the control input, the oxygen transfer coefficient *KLa*(*t*).

The states and input vectors are thus given as:

$$\mathbf{x}(t) = \left[ \mathbf{X}\_{COD}(t), \mathbf{S}\_O(t), \mathbf{X}\_{BH}(t), \mathbf{S}\_{NH}(t), \mathbf{S}\_{NO}(t), \mathbf{X}\_{BA}(t) \right]^T$$

$$\boldsymbol{\mu}(t) = \mathbf{K}\_{La}(t)$$

The WWTP process is therefore described by the following dynamic equations of the reduced model:

$$\mathcal{X}\_{\rm COD}(t) = \frac{1}{Y\_{\rm h}} \left[ \theta\_1(t) + \theta\_2(t) \right] + \left( 1 - f\_p \right) \left( \theta\_4(t) + \theta\_5(t) \right) + \theta\_1(t), \tag{1}$$

$$
\dot{S}\_{\rm O}(t) = \frac{Y\_h - 1}{Y\_h} \theta\_1(t) + \frac{Y\_d - 4.57}{Y\_d} \theta\_3(t) + \theta\_2(t),
\tag{2}
$$

$$\dot{S}\_{NH}(t) = -i\_{\rm xb} \left[ \theta\_1(t) + \theta\_2(t) \right] - \left[ i\_{\rm xb} + \frac{1}{Y\_a} \right] \theta\_3(t) + \left( i\_{\rm xb} - f\_p i\_{\rm xp} \right) \left[ \theta\_4(t) + \theta\_5(t) \right] + \theta\_3(t), \tag{3}$$

$$\mathcal{S}\_{\rm NO}(t) = \frac{Y\_{\hbar} - 1}{2.86Y\_{\hbar}} \theta\_2(t) + \frac{1}{Y\_{\hbar}} \theta\_3(t) + \theta\_4(t), \tag{4}$$

$$
\mathcal{X}\_{BH}(t) = \theta\_1(t) + \theta\_2(t) - \theta\_4(t) + \theta\_5(t),
\tag{5}
$$

$$
\dot{X}\_{BA}(t) = \theta\_3(t) - \theta\_5(t) + \theta\_6(t). \tag{6}
$$

where

$$\begin{aligned} \theta\_1(t) &= \mu\_h \frac{\mathcal{X}\_{\text{COD}}(t)}{\mathcal{K}\_{\text{COD}} + \mathcal{X}\_{\text{COD}}(t)} \frac{\mathcal{S}\_{\text{O}}(t)}{\mathcal{K}\_{\text{OH}} + \mathcal{S}\_{\text{O}}(t)} \mathcal{X}\_{\text{BH}}(t) \\ \theta\_2(t) &= \mu\_h \eta\_{NO} \frac{\mathcal{X}\_{\text{COD}}(t)}{\mathcal{K}\_{\text{COD}} + \mathcal{X}\_{\text{COD}}(t)} \frac{\mathcal{S}\_{\text{NO}}(t)}{\mathcal{K}\_{\text{NO}} + \mathcal{S}\_{\text{NO}}(t)} \frac{\mathcal{K}\_{\text{OH}}}{\mathcal{K}\_{\text{OH}} + \mathcal{S}\_{\text{O}}(t)} \mathcal{X}\_{\text{BH}}(t) \\ \theta\_3(t) &= \mu\_d \frac{\mathcal{S}\_{\text{NH}}(t)}{\mathcal{K}\_{\text{NH},A} + \mathcal{S}\_{\text{NH}}(t)} \frac{\mathcal{S}\_{\text{O}}(t)}{\mathcal{K}\_{\text{O},A} + \mathcal{S}\_{\text{O}}(t)} \mathcal{X}\_{\text{BA}}(t) \\ \theta\_4(t) &= b\_H \mathcal{X}\_{\text{BH}}(t) \\ \theta\_5(t) &= b\_A \mathcal{X}\_{\text{BA}}(t) \end{aligned}$$

With the flow rate given as *Qin*(*t*), *Vo* as the volume of the aerobic tank and considering that *S*0*in* (*t*), *SNOin* (*t*), *XBAin* (*t*) are equal to zero. *ϑ*1(*t*), *ϑ*2(*t*), ··· , *ϑ*6(*t*) are given as follows:

$$
\begin{split}
\theta\_{1}(t) &= \frac{\underline{Q\_{in}}(t)}{V\_{o}} \left[ \mathbf{X}\_{\text{COD}\_{\text{in}}}(t) - \mathbf{X}\_{\text{COD}}(t) \right] \\
\theta\_{2}(t) &= \frac{\underline{Q\_{in}}(t)}{V\_{o}} \left[ -S\_{O}(t) \right] + K\_{Lx}(t) \left[ S\_{O\_{\text{sdf}}} - S\_{O}(t) \right] \\
\theta\_{3}(t) &= \frac{\underline{Q\_{in}}(t)}{V\_{o}} \left[ S\_{N\underline{H}\_{\text{in}}}(t) - S\_{N\underline{H}}(t) \right] \\
\theta\_{4}(t) &= \frac{\underline{Q\_{in}}(t)}{V\_{o}} \left[ -S\_{NO}(t) \right] \\
\theta\_{5}(t) &= \frac{\underline{Q\_{in}}(t)}{V\_{o}} \left[ X\_{BH\_{\text{in}}}(t) - \frac{f\_{\text{w}}(1+f\_{r})}{f\_{r}+f\_{w}} X\_{BH}(t) \right] \\
\theta\_{6}(t) &= \frac{\underline{Q\_{in}}(t)}{V\_{o}} \left[ X\_{BA}(t) - \frac{f\_{\text{w}}(1+f\_{r})}{f\_{r}+f\_{w}} X\_{BA}(t) \right]
\end{split}
$$

where *YH*, *YA*, *fr*, *fw*, *bh*, *bA*, *ixb*, and *fp* are the stoichiometric parameters and *μh*, *KCOD*, *KOH*, *μa*, *KNH*,*A*, and *KO*,*A* are the kinetic parameters.

**Figure 2.** Layout of Girona WWTP.

#### *2.3. LPV Representation of the WWTP*

For ease of computational burden, the nonlinear reduced model is represented in a LPV form which involves the embedding of nonlinearities in varying parameters, resulting in a linear representation in state space. This procedure offers benefits when applied to MPC over its nonlinear MPC [15] and linear MPC [17] counterparts as applied on the WWTP by providing a faster run time and the avoidance of numerical problems with respect to the former and the ability to operate in a wide range of operating points with regard to the latter. The nonlinear model in this case is defined by linear systems at each time instance based on some time-varying parameters *σ*(*t*) ∈ R*nσ* , with an assumption that the parameters *σ*(*t*) are not known a priori but can be measured or estimated online [18]. The dynamic behavior of the LPV model is therefore described as:

$$
\dot{\mathbf{x}}(t) = A(\sigma(t))\mathbf{x}(t) + B(\sigma(t))\boldsymbol{u}(t) \tag{7}
$$

$$y(t) = \mathbb{C}(\sigma(t))x(t) + D(\sigma(t))u(t) \tag{8}$$

where *x*(*t*) ∈ R*nx* and *u*(*t*) ∈ R*nu* are the states and inputs, respectively, with *y*(*t*) ∈ R*ny* as the measured signals. *A*(*σ*(*t*)), *B*(*σ*(*t*)), *C*(*σ*(*t*)) and *D*(*σ*(*t*)) are time-varying matrices of appropriate dimensions that are affine in *σ*(*t*) ∈ R*nσ* . In the quasi LPV case, the scheduling parameters are dependent on measured signals, *ys*(*t*) ∈ R*k* ⊂ *y*(*t*) ∈ R*ny* , such that

$$
\sigma(t) = f(y\_s(t))
$$

where *f* : R*k* → R*nσ* is a continuous mapping [19]. With observed states and exogenous inputs (*w*(*t*)), nonlinearities involving the system states can be "hidden" in the varying parameters, *σ*(*t*, *ys*(*t*), *w*(*t*)).

Therefore, from the generic nonlinear form

$$\begin{aligned} \dot{x}(t) &= f(x(t), u(t), w(t)) \\ y(t) &= g(x(t), u(t)) \end{aligned} \tag{9}$$

a linear quadruple (*A*(*σ*(*t*, *ys*(*t*), *w*(*t*))), *B*(*σ*(*t*, *ys*(*t*), *w*(*t*))), *C*(*σ*(*t*, *ys*(*t*), *w*(*t*))), *and D*(*σ*(*t*, *ys*(*t*), *w*(*t*))) estimate is formulated and incorporated into the EMPC for a convex optimization problem. In the following, the function *σ*(*t*, *ys*(*t*), *w*(*t*))) will simply be represented as *σ*(*t*). The choice of scheduling parameters considering the origin of nonlinearities in the reduced model (1) are

$$\begin{aligned} \sigma\_1(t) &= Q\_{\text{in}}(t), \qquad \sigma\_2(t) = \frac{X\_{\text{COD}}(t)}{K\_{\text{COD}} + X\_{\text{COD}}(t)} \frac{X\_{BH}(t)}{K\_{OH} + S\_O(t)},\\ \sigma\_3(t) &= \frac{X\_{\text{COD}}(t)}{K\_{\text{COD}} + X\_{\text{COD}}(t)} \frac{S\_{\text{NO}}(t)}{K\_{\text{NO}} + S\_{\text{NO}}(t)} \frac{K\_{\text{OH}}}{K\_{\text{OH}} + S\_O(t)},\\ \sigma\_4(t) &= \frac{1}{K\_{\text{OA}} + S\_O(t)} \frac{S\_{NH}(t)}{K\_{NH} + S\_{NH}(t)} X\_{BH}(t), \qquad \sigma\_5(t) = S\_O(t). \end{aligned}$$

The dynamic LPV model is thus given as:

$$
\dot{\mathbf{x}} = A(\sigma(t))\mathbf{x}(t) + B(\sigma(t))u(t) + Ew(t). \tag{10}
$$

with the time-varying matrices, *A*(*σ*(*t*)), *B*(*σ*(*t*)), and time-invariant disturbance matrix *E* as:

$$A(\sigma(t)) = \begin{bmatrix} a\_{11}(t) & 0 & 0 & 0 & a\_{15}(t) & a\_{16} \\ 0 & a\_{22}(t) & 0 & 0 & a\_{25}(t) & 0 \\ 0 & a\_{32}(t) & 0 & 0 & a\_{35}(t) & a\_{36} \\ 0 & a\_{42}(t) & 0 & 0 & a\_{45}(t) & 0 \\ 0 & 0 & 0 & 0 & a\_{55}(t) & 0 \\ 0 & a\_{62}(t) & 0 & 0 & 0 & a\_{66}(t) \end{bmatrix},$$

$$B(\sigma(t)) = \begin{bmatrix} 0 \\ b\_{12}(t) \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} \text{and } E = \frac{1}{V\_{\odot}} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}.$$

where

$$\begin{aligned} a\_{11}(t) &= \frac{\sigma\_1(t)}{V\_o}, & a\_{15}(t) &= -\frac{\mu\_h}{Y\_h}\sigma\_2(t) + (1-f\_f)b\_h - \frac{\mu\_h \eta\_{NO\_3}}{Y\_h}\sigma\_3(t), & a\_{16} &= (1-f\_f)b\_{a\_1}, \\ a\_{22}(t) &= -\frac{\sigma\_1(t)}{V\_O} - \frac{4.57 - Y\_a}{Y\_a}\mu\_4\sigma\_4(t), & a\_{25}(t) &= \frac{Y\_h - 1}{Y\_h}\mu\_h\sigma\_2(t) \\ a\_{32}(t) &= -\left(i\_{xb} + \frac{1}{Y\_a}\right)\mu\_{47}\sigma\_4(t), & a\_{35}(t) &= (i\_{xb} - f\_{p1x})b\_{h} - i\_{xb}\mu\_h\sigma\_2(t) - (i\_{xb}\mu\_h\sigma\_3(t)), \\ a\_{36} &= \left(i\_{xb} - f\_{p1x}b\right)b\_{a}, & a\_{42}(t) &= \frac{1}{Y\_a}\mu\_{47}(t), & a\_{45}(t) &= \frac{Y\_h - 1}{2.86Y\_h}\mu\_h\eta\_{NO\_3}\sigma\_3(t), \\ a\_{55}(t) &= \mu\_h\sigma\_2(t) - b\_{h} - \left(\frac{\sigma\_1(t)}{V\_O} - \frac{f\_{w1}(1+f\_r)}{f\_r+f\_w}\right) - b\_{a}, & a\_{62}(t) &= \mu\_a\sigma\_4(t), \\ a\_{66}(t) &= \frac{\sigma\_1(t)}{V\_O}\left(-\frac{f\_{w1}(1+f\_r)}{f\_r+f\_w} - 1 - b\_{a}\right), & b\_{12}(t) &= S\_{\rm sat} - \sigma\_5(t). \end{aligned}$$

The input concentrations are

$$w(t) = \begin{bmatrix} Q\_{in}(t)X\_{COD\_{in}}(t) & Q\_{in}(t)S\_{NH\_{in}}(t) & Q\_{in}(t)X\_{BH\_{in}}(t) \end{bmatrix}^T$$

.

**Remark 1.** *In this work, it is assumed that all the concentrations are measured online, but it must be noted that in practice, not all the concentrations, such as, e.g., XCODin can be measured.*

#### **3. EMPC of a WWTP**

#### *3.1. Operational Goals*

The immediate control goal of a WWTP is to meet water quality levels established by regulators while operating efficiently by reducing operational cost. As discussed in the introduction, predictive control techniques may be used to compute strategies which achieve this goal while at the same time optimizing the system performance in terms of different operational indices. To achieve this objective, the control of dissolved oxygen concentration as well as nitrates within certain limits is necessary. The MPC presents the advantage of being a non-conservative control strategy, such that in periods of low influents, with a minimal level of pollutants, the effluent quality can be achieved by regulating the levels of *S*0 and *SNO* below the stipulated reference point to avoid waste of energy. Subsequently, during periods of high influents levels, it is then important to meet the predefined set points to reduce pollutants, avoiding the violation of the standard effluent quality set by authorities [15]. In this work, a PI-EMPC control strategy is employed: PI designed by authors of the BSM1 for the regulation of *SNO* and a designed EMPC for the control of *SO* in the aeration tank. In the proposed LPV EMPC, the following objectives are then considered:

• **Economic costs.** The main economic costs associated with WWTP are primarily due to treatment and electricity costs. Water through the WWTP involves important electricity costs in pumping stations in charge of internal and external water recirculations as well as aeration in the aerobic tanks. In our case, only the aeration energy is considered with an objective of minimizing the cost associated with supply of oxygen for controlled culture growth. The performance index is described as follows

$$J\_{\rm esc}(k) = \frac{S\_{\rm ost}}{1800} V\_o K\_{La}(k) \left[\frac{\rm kwh}{\rm day}\right]. \tag{11}$$

• **DO concentration control.** In order to control the *So* within some bounds in the EMPC during the aeration process, slack variables are introduced in the optimization problem, which seek to penalize the dissolved oxygen states, such that they are maintained in a range to maintain effluent quality. Selecting slack variables, (*λ*+ > 0 and *λ*− > 0), additional terms of soft constraints (see (16c) and (16d)) and a quadratic objective index are introduced with *xsp* as the selected DO concentration value. The introduction of the slack variables ensures that the DO concentration varies within a boundary around *xsp* aided by the appropriate selection of weights in the objective function. The performance index is thus given as

$$J\_{\lambda}(k) = \|\lambda(k)\|\_{2'}^{2} \tag{12}$$

where *λ*(*k*) = *λ*−(*k*), *λ*+(*k*))*T* .

• **Smooth set points for equipment conservation.** The operation of WWTP and main valves and pumps usually requires smooth flow set-point variations. To obtain such a smoothing effect, the proposed MPC controller includes a third term in the objective function to penalize the control signal variation between consecutive time intervals. This term is expressed as

$$J\_{smo}(k) = \Delta u(k)^T \mathcal{W}\_u \Delta u(k). \tag{13}$$

Therefore, the performance function *J* considering the aforementioned control objectives has the form

$$J = w\_1 \sum\_{k=0}^{H\_p - 1} J\_{\text{coo}}(k) + w\_2 \sum\_{k=0}^{H\_p - 1} J\_{\text{smo}}(k) + w\_3 \sum\_{k=1}^{H\_p} J\_\lambda(k). \tag{14}$$

#### *3.2. Control Strategy Computation*

The control strategy is determined by the computation of an optimal sequence of control actions for a prediction horizon, *Hp*.

$$\mathfrak{u}\_k = \left(\mathfrak{u}(k|j)\right)\_{j=0}^{H\_{p-1}} = \left(\mathfrak{u}(k|0), \mathfrak{u}(k|1), \dots, \mathfrak{u}(k|H\_{p-1})\right). \tag{15}$$

We solve at each time instance *k*, the following optimal control problem with initial state obtained from measurements (or state estimation) of the dynamics WWTP model and prediction in the MPC loop with the qLPV plant model (10),

$$\min\_{\vec{u}\_k} J(\vec{u}\_k, k) \tag{16a}$$

subject to

$$\mathbf{x}(i+1|k) = A(\sigma(k))\mathbf{x}(i|k) + B(\sigma(k))\mathbf{u}(i|k) + E\mathbf{w}(i|k) \quad i = 0, \dots, H\_p - 1,\tag{16b}$$

$$\mathbf{x}\_{\mathcal{S}\_0}(i|k) > = \mathbf{x}\_{\mathcal{sp}} + \lambda^+(i|k), \qquad i = 1, \cdots, H\_{\mathcal{P}'} \tag{16c}$$

$$\mathbf{x}\_{\mathfrak{s}\_0}(i|k) <= \mathbf{x}\_{sp} - \lambda^-(i|k), \qquad i = 1, \cdots, \mathbf{l}, \mathbf{H}\_{\mathfrak{p}}.\tag{16d}$$

$$u(i|k) \in \mathcal{U} \quad i = 0, \cdots, H\_p - 1,\tag{16e}$$

$$\mathbf{x}(i|k) \in \mathcal{X} \qquad i = 1, \cdots, H\_{\mathbf{p}} \tag{16f}$$

$$\mathcal{Y}(i|k) \in \mathcal{Y} \qquad j = 0, \cdots, H\_{\mathbb{P}}, \tag{16g}$$

$$
\lambda^+(i|k), \lambda^-(i|k) \;> = 0 \tag{16h}
$$

where *xso* is the dynamic state representing the soluble oxygen. (16c–f) are described by the box constraints:

$$\begin{array}{l} \mathcal{M} = \left\{ \boldsymbol{\mu} \in \mathbb{R}^{n\_{\boldsymbol{\mu}}} | \boldsymbol{\mu}^{\text{min}} \le \boldsymbol{\mu} \le \boldsymbol{\mu}^{\text{max}} \right\}, \\ \mathcal{X} = \left\{ \boldsymbol{x} \in \mathbb{R}^{n\_{\boldsymbol{\mu}}} | \boldsymbol{x}^{\text{min}} \le \boldsymbol{x} \le \boldsymbol{x}^{\text{max}} \right\}, \\ \mathcal{Y} = \left\{ \boldsymbol{y} \in \mathbb{R}^{n\_{\boldsymbol{\mu}}} | \boldsymbol{y}^{\text{min}} \le \boldsymbol{y} \le \boldsymbol{y}^{\text{max}} \right\}. \end{array} \tag{17}$$

which are determined from the maximum residual concentrations imposed in order to cope with the European Union effluent standards on chemical oxygen demand COD, suspended solids *SS* and total nitrogen *TN*:

$$COD \lessapprox COD\_{\max} = 125 \,\text{gm}^{-3}/$$

$$S\_S \lessapprox S\_{S\max} = 35 \,\text{gm}^{-3}/$$

$$T\_N \lessapprox T\_{N\max} = 10 \,\text{gm}^{-3}.$$

The first control action of the sequence *u*( *k*|0) is applied to the WWTP plant to obtain the system measurements and/or MHE estimated states, which are then used in the succeeding optimization problem, resulting in a recursive procedure. Not all the state variables are measured as stated earlier; the moving horizon estimator (MHE), which is the dual of the MPC controller, estimates the unmeasurable states.

#### **4. Moving Horizon Estimation**

Since some states cannot be measured online in the operation of the WWTP, a design of an estimator, in our case the MHE, is necessary for the prediction of system outputs, bearing in mind that apart from purposes of feedback control, the quasi-LPV formulation relies on information of the system states for the model construction. By solving a constrained optimization problem, the MHE utilizes a limited *N*-prediction horizon of past measurements through an error minimization scheme aided by information of the system model in a prediction window to estimate the system states. The optimization problem is therefore set up with the discretized plant model as:

$$\begin{aligned} \min\_{\{\hat{\mathbf{x}}(i|k)\}\_{i=-N}^{0}} & \quad \left(\hat{\mathbf{x}}(-N|k)-\mathbf{x}\_{o}\right)^{T} \boldsymbol{P}\_{o} \Big(\hat{\mathbf{x}}(-N|k)-\mathbf{x}\_{o}\Big) + \sum\_{i=-N}^{k} \left(\boldsymbol{\varepsilon}(i|l)^{T}\mathbf{Q}\boldsymbol{\varepsilon}(i|k) + \boldsymbol{s}(i|k)^{T}\mathbf{R}\boldsymbol{s}(i|k)\right) \\ \text{s.t.} & \quad \hat{\mathbf{x}}(i+1|k) = A(\boldsymbol{\sigma}(i|k))\hat{\mathbf{x}}(i|k) + B(\boldsymbol{\sigma}(i|k))u(i|k) + E\boldsymbol{w}(i|k) + \boldsymbol{\epsilon}(i|k) \quad i = -N\_{r} \cdots -1, \\ & \quad y(i|k) = \mathbf{C}\boldsymbol{x}(i|k) + s(i|k), \\ & \quad \pounds\_{k} \in \mathcal{X}. \end{aligned} \tag{18}$$

where *R* = *RT* ∈ R*ny*×*ny* > 0, *Q* = *QT* ∈ R*nx*×*nx* ≥ 0 and *Po* = *PT o* ∈ R*nx*×*nx* ≥ 0 are the weighting matrices that are defined according to uncertainty levels induced respectively by the noise, disturbance and unknown initial conditions (*xo*). X bounds the estimated states. At every iteration, *N* sets of control inputs, {*u*(*i*|*k*)}−1 *i*=−*N* ∈ R*nu*×*N*, measurements {*y*(*i*|*k*)}−1 *i*=−*N* ∈ R*ny*×*N* and *N* sets of LPV matrices {*Ai*}−1 *i*=−*N* ∈ R(*nx*×*nx*)*N*, {*Bi*}−1 *i*=−*N* ∈ R(*nx*×*nu*)*N* are taken as inputs into the optimization problem to predict the state sequence {*x*ˆ(*i*|*k*)}0 *i*=−*N* ∈ R*nx*×(*N*+1) by solving the dynamical optimization problem (18). The last element of the sequence {*x*ˆ(*i*|*k*)}0 *i*=−*N* is subsequently chosen as the estimated states, the measurements and inputs are then discarded, and the procedure is repeated. The ammonia concentration (*SNH*), nitrate concentration (*SNO*) and the soluble oxygen (*So*) are supposedly measurable; therefore, the MHE is designed for the estimation of [*XCOD*, *XBH*, *XBA*] *T* as shown in Figures 3–5 .

**Figure 3.** MHE estimate of oxygen demand concentration (*XCOD*) for 7 days.

**Figure 4.** MHE estimate of heterotrophic biomass (*XBH*) for 7 days.

**Figure 5.** MHE estimate of autotrophic biomass (*XBA*) for 7 days.

#### **5. Simulation Results**

#### *5.1. LPV EMPC Implementation Details*

To illustrate LPV EMPC approach presented in this paper, the Girona WWTP case study presented in Section 2 is used. The constituents of the influent wastewater of Girona WWTP varies during the day between the following bounds :

The inflow of Girona WWTP is shown in Figure 6.

**Figure 6.** WWTP inflow.

With a quasi-linear approximation of the nonlinear WWTP via the LPV representation, the constrained optimization problem (16) is solved using quadratic programming formulation using the CPLEX® solver in MATLAB® on an Intel Core i7, 8 GB of RAM PC. A sampling time of 15 min and a prediction horizon of 6 h is chosen for simulation. The process is simulated for 7 days in a Simulink environment representing the dynamics of the Girona WWTP, as shown in Figure 1.

Using the weights *wi* associated with the multiobjective EMPC cost function, (14) is tuned using the procedure as performed in [20,21] with the aim of maintaining the quality of the exit water at some levels within the current regulations regardless of the entry at a minimum cost.

Some control scenarios are selected to show different behaviors of the proposed scheme by altering *Xsp* and manipulating weights *wi*, ideally to illustrate the different actions of aeration corresponding to different dissolved oxygen requirements for a quality effluent.

#### *5.2. First Scenario*

The first scenario consists of controlling the dissolved oxygen concentration in the exit of the biological treatment plant between the bounds (1.5, 2.5). Figure 7 shows the dynamics of the DO concentration (above) and its corresponding aeration energy (below). The operation of the aeration, as stated in the preceding section, corresponds to the variation of the influents during the day; therefore, the DO concentration varies between the defined bounds in relation to the amount of pollutants at each time instance in the influents, which can be inferred from Figure 6.

**Figure 7.** (**Above**): DO concentraton variation. (**Below**): Aeration flow for Scenario 1.

#### *5.3. Second Scenario*

The second scenario also consists of controlling the DO concentration between the ranges of 0.5 to 1.2 mg/L with minimum aeration energy consumption.

From Figure 8, a similar behavior of oxygen in the tanks as in the first scenario is realized with an expected less aeration energy, as less DO is required for treatment. The nitrates in the exit of the WWTP (Figure 9) range approximately between 5 and 7 mg/L.

**Figure 8.** (**Above**): DO concentration variation and (**Below**): Aeration flow for Scenario 2.

**Figure 9.** Nitrate concentration variation.

#### **6. Conclusions**

In this paper, an LPV EMPC strategy for the control of dissolved oxygen concentration in the aerated reactors of a WWTP is proposed and applied to the Girona (Spain) case study. The proposed approach combines two improvements with respect to the existing approaches in the literature: First, differently from standard tracking MPC, the proposed EMPC strategy optimizes the economic performance of the plant instead of following some pre-established set points. Second, a reduced model of the WWTP is represented in a quasi-LPV form allowing the real-time implementation of the controller thanks to the use of quadratic programming optimization tools. If otherwise, the nonlinear model plant was used, nonlinear programming algorithms are required that usually prevent the realtime implementation because of the large computational time. Moreover, an LPV moving horizon state estimation scheme has also been proposed that allows the implementation of the LPV EMPC with the available sensors in the WWTP. The effectiveness of the proposed scheme has been illustrated in the considered case study with two scenarios aiming at keeping the DO within some bounds.

As future work, real testing in the WWTP plant will be conducted to further validate the performance of the proposed solution. Another issue to take into consideration is the application of the proposed methodology for aerobic conditions maintenance in sewer networks [22].

**Author Contributions:** Conceptualization, F.N., B.K. and V.P.; methodology, F.N., B.K. and V.P.; software, B.K. and J.P.; writing—original draft preparation, F.N.; writing—review and editing, F.N., B.K, V.P., J.P., J.Q. and S.d.C.; supervision, F.N., B.K and V.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been co-financed by the Spanish State Research Agency (AEI) and the European Regional Development Fund (ERFD) through the project SaCoAV (ref. MINECO PID2020- 114244RB-I00), by the European Regional Development Fund of the European Union in the framework of the ERDF Operational Program of Catalonia 2014–2020 (ref. 001-P-001643 Looming Factory), and by the DGR of Generalitat de Catalunya (SAC group ref. 2017/SGR/482).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the required data are included in the paper.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

## *Article* **Comparison of Optimisation Algorithms for Centralised Anaerobic Co-Digestion in a Real River Basin Case Study in Catalonia**

**David Palma-Heredia 1,\*, Marta Verdaguer 1, Vicenç Puig 2,3, Manuel Poch 1 and Miquel Àngel Cugueró-Escofet 2**

**Abstract:** Anaerobic digestion (AnD) is a process that allows the conversion of organic waste into a source of energy such as biogas, introducing sustainability and circular economy in waste treatment. AnD is an intricate process because of multiple parameters involved, and its complexity increases when the wastes are from different types of generators. In this case, a key point to achieve good performance is optimisation methods. Currently, many tools have been developed to optimise a single AnD plant. However, the study of a network of AnD plants and multiple waste generators, all in different locations, remains unexplored. This novel approach requires the use of optimisation methodologies with the capacity to deal with a highly complex combinatorial problem. This paper proposes and compares the use of three evolutionary algorithms: ant colony optimisation (ACO), genetic algorithm (GA) and particle swarm optimisation (PSO), which are especially suited for this type of application. The algorithms successfully solve the problem, using an objective function that includes terms related to quality and logistics. Their application to a real case study in Catalonia (Spain) shows their usefulness (ACO and GA to achieve maximum biogas production and PSO for safer operation conditions) for AnD facilities.

**Keywords:** anaerobic co-digestion; ant colony optimisation; particle swarm optimisation; genetic algorithms; waste management; circular economy

#### **1. Introduction**

In the context of global climate change with rising and more extreme events—such as droughts and floods—which will likely provide growing uncertainty to water demand and jeopardise the availability of specific resources, there is a growing interest in the adaptation and use of technologies related to the circular economy that promote environmental sustainability. In this framework, resource recovery is a key issue for industrial and environmental processes and shows a wide spectrum of study possibilities. In water sanitation, wastewater treatment plants (WWTPs) offer a wide range of possibilities for resource recovery, mainly related to sludge treatment processes [1–7] as biogas generation via the substrate co-digestion process, which can be an alternative source for thermal and electrical energy production [8–14]. This potential for biogas generation could be translated as well to a source of renewable natural gas, which has specific composition requirements that demand high-tech sensors to assure its quality no matter its origin, as those developed in [15,16]. Due to their potential for resource recovery and the further implications in the water–food–energy nexus, WWTPs have been a research focus from different areas of

**Citation:** Palma-Heredia, D.; Verdaguer, M.; Puig, V.; Poch, M.; Cugueró-Escofet, M.À. Comparison of Optimisation Algorithms for Centralised Anaerobic Co-Digestion in a Real River Basin Case Study in Catalonia. *Sensors* **2022**, *22*, 1857. https://doi.org/10.3390/s22051857

Academic Editor: Ernest W. Tollner

Received: 27 January 2022 Accepted: 24 February 2022 Published: 26 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

expertise: from modelling and engineering design [17–24] to process dynamics, simulation and integration [25–28].

Anaerobic digestion (AnD), a complex process involved in biogas production, has a delicate balance of substrate composition. The optimal performance requires avoiding process inhibition and maximising biogas generation. The optimal balance may be achieved with the correct mixture of available substrates, but this task is challenging and difficult to achieve manually due to the high combinatorial possibilities and changing availability of substrates of heterogeneous nature. The complexity increases when the process is codigestion, with the addition of residual substrates produced by agro-food and similar industries, each with its own dynamics of substrate generation and composition [29–32]. Additionally, not all WWTPs have an anaerobic digester. Therefore, optimisation also requires logistical challenges to process a maximum volume of the available substrates in a certain geographical area and its travel logistics (of sludge and co-substrates) from its origin to the destination digester.

Hence, dealing with such complexity is a former step to tackle optimal co-digestion in a complex network composed of many substrate sources—including WWTPs without AnD processes and industrial producers—and several co-substrate receptors. These will be located in different geographical places. As a result, the logistics of substrates will be affected by the geographical distance between actors involved and the restrictions related to the receptors.

Optimisation of the individual digester feed requires optimal blending of different co-substrates in order to fulfil the volumetric and compositional requirements of the anaerobic procedure. This problem can be understood as a multidimensional knapsack problem (MKP) [33–35]. The MKP is an NP-hard problem [36] and has been widely studied in the literature. To solve this type of problem, the use of combinatorial optimisation metaheuristics is proposed in [37,38], mainly when a high number of restrictions are presented [39].

Many tools have been developed to this end, either focused on modelling, control of the optimum co-substrate blending, or system operation, as shown in [37,39–43]. In [37,40], identification and modelling of critical parameters are performed; in [39] control schemes based on the composition qualities are developed; and in [42,43], optimised control strategies are implemented according to blend composition. In [41], logistics are also included to optimise the performance of a single anaerobic digester with co-digestion strategies.

However, in real-world installations, most of these systems are managed and supervised not in a single fashion but in a network fashion. Thus, proper system management requires simultaneous consideration of the entire AnD network to select which combination is the best for each digester to maximise the potential of the overall infrastructure. Besides, literature on this matter is relatively scarce due to its ad-hoc nature. There is literature related to optimised placement of new AnD plants, such as in [44], but it lacks optimisation of the operational part involved in the feeding of the anaerobic digesters. Very specific works can be found about optimisation of supply chain networks in the field of waste valorisation, such as in [45], where an integrated geographical information system (GIS)-based optimisation is performed, but it requires highly detailed and tailored data, so its implementation becomes time-consuming and highly dependent on data availability; furthermore, it does not tackle process optimisation regarding waste processing facilities. Regarding logistics, other works can be found for path planning optimisation such as in [46], where truck routes are traced based on GIS-oriented algorithms, or in [47], where a smart waste bin prototype is developed for sensor-based waste classification. As it can be seen, there is a gap in the literature regarding network optimisation of existing waste management facilities (such as AnD plants) that would include both logistics (i.e., minimising route impact and length) and quality (i.e., improving process performance) optimisation. This is a necessary gap to fulfil, since, as stated before, AnD networks are currently managed in an ad-hoc, manual fashion by the practitioners, which is dramatically time-consuming and needs highly qualified personnel. Although there are currently different approaches

that could help overcoming specific parts of this challenge (i.e., those observed in [44–47]), none of them can currently successfully accomplish the overall task.

The approach presented in this study continues the work introduced in [48] where the optimisation problem of blending in anaerobic co-digestion (AnCD) is handled by an ant colony optimisation (ACO) algorithm in a synthetic case study with realistic conditions, simulating a centralised AnD single-stage reactor that received feed once a day. In [41], the work is extended while considering the quality, social and travel logistics of the cosubstrates, analysing its importance to the overall optimisation. In addition, ACO has been implemented in real-world waste sector case studies, e.g., [48,49]. Here, the work is extended to a similar real-world case study considering multiple receptors in the geographical area of the Besòs River basin (Catalonia, Spain). The data used correspond to the real operation conditions of this area. For these conditions, the authors have also evaluated the results obtained using different optimisation approaches such as ACO, genetic algorithms (GA) and particle swarm optimisation (PSO).

ACO, GA and PSO algorithms were selected as convenient approaches to tackle a problem of the nature stated here after reviewing applications of similar nature in the literature. In [50], a review of nature-inspired algorithms is performed, including GA, ACO and PSO—among others—for AnD modelling and optimisation, showing how in the field of AnD, PSO obtains better performance in substrate feed optimisation for agricultural biogas plants than other evolutionary methods. For example, in [51] genetic algorithms are used to minimize the environmental impact caused by mine water. The main drawback for PSO pointed out in [50] is premature convergence, since particles may become trapped in local optima or suffer stagnation, but this may be solved by a partial restart of the process introducing new particles in the search space. In [52], ACO and GA are applied to optimise the route of waste collection vehicles for municipal waste collection and transportation—the highest cost of the entire waste management system—with similar performance attained by both algorithms; however, only the problem of travel logistics is considered, not the blending of municipal waste. Ref. [53] proposes a nonlinear model predictive control strategy using the MATLAB BioOptim toolbox, developed by the same authors, for optimal control of substrate feed to AnD operation of an agricultural biogas plant, with a graphical user interface (GUI) integrating a fitness function including different operating constraints and parameters such as pH, solids or methane concentration, and using evolutionary optimisers such as PSO, covariance matrix adaptation evolution strategy (CMAES). Alternatively, they propose differential evolution (DE), which achieved better performance with PSO but without considering the substrate travel logistics in the optimisation. Ref. [54] presents a prediction and optimisation method using a multi-layer perceptron artificial neural network (ANN) and PSO for the maximisation of biogas generation in a real wastewater treatment facility. A similar approach is presented in [55], where modelling and optimisation of biogas production with mixed substrates are obtained with a combination of ANN and GA methods.

Additionally, ref. [55] and references therein point out how stochastic global optimisation algorithms (SGOAs), such as PSO, the ACO, and GA, among others, are considered efficient alternatives in the design of optimal production media and optimal process operating conditions in fermentation research and can significantly reduce the process development time. Regarding the comparison of the optimisation algorithms selected here solving combinatorial optimisation (CO) problems, in [56], ACO and GA are compared, both achieving good performance but with GA exhibiting slightly better performance than ACO. In the latter reference, it is also mentioned how trimming of specific parameters for both optimisers—e.g., number of iterations, evaporation coefficient and number of ants for ACO, or chromosome population, crossover and mutation probabilities for GA—is required to achieve good performance in both cases. In [57], relationships between GA and ACO-type algorithms are detailed, presenting their similitudes and showing how they use similar principles to succeed in CO problems with globally convex structure of its solution space. Overall, SGOAs such as ACO, GA or PSO have shown good performance

in the type of applications presented and hence demonstrated suitability for non-convex nonlinear multidimensional optimisation problems, as presented here. Optimal blending for AnCD is considered, e.g., in [40,55], but to the knowledge of the authors, the optimisation of such blending combined with the travel logistics of co-substrates in a centralised multi-receptor co-digestion strategy has not yet been studied. Additionally, a comparison between different suitable optimisation algorithms for such applications, i.e., ACO, GA and PSO, is presented.

The application of optimisation strategies in AnD allows a significant enhancement of co-digestion strategies [30] maximising biogas production and minimising associated risks to each AnD operation (e.g., overdosing or acidification). In this work, the performance of each optimisation approach considered is evaluated on a real case study in the area of the Besòs River basin in Catalonia, including a network of substrate generators and three anaerobic digesters. Hence, the objective of this study is to develop a tool that is able to optimise the centralised digestion process of an AnD network with multiple waste sources and waste receptors by means of three evolutionary optimisation algorithms—namely, ACO, GA and PSO. Such a tool is tested in a real case study to further analyse and compare the performance of each algorithm in the overall AnD network optimisation.

#### **2. Material and Methods**

#### *2.1. Optimisation Algorithms Considered*

The optimisation algorithms presented here fall within the set of SGOA, concretely in the subset of evolutionary algorithms (EA) for GA—which use mechanisms inspired by biological evolution, e.g., mutation or recombination to achieve the goal of optimisation and in the subset of swarm intelligence methods for ACO and PSO based on the collective behaviour of self-organised decentralised systems, respectively. SGOA algorithms have been widely used to solve NP-hard combinatorial optimisation problems, such as that presented in this study, which deterministic optimisation methods fail to handle due to their complexity.

Regarding each proposal, ACO is a metaheuristic approach that has been shown to be effective in solving a variety of NP-hard problems [58]. The algorithm is based on simulation of the behaviour of real ants in their search for food. When ants find food, they leave a pheromone trail on their path. Then, new ants follow that trail. In this way, an increasing number of ants are concentrated in places where there is food. In a similar way, the virtual ants construct a solution moving through the graph that represents the search space of solutions. Their paths are guided by a probabilistic state transition rule, which is based on pheromone trails and specific heuristic information. The algorithmic procedure is iterative. At each iteration, the pheromone trails are updated by applying an evaporation coefficient (when the value selected is not part of a feasible solution). To avoid rapid stagnation of the solution, the ACO algorithms can use several strategies [59], such as that related to the limitation of the pheromone trails between maximum and minimum values. Max–Min Ant System [60] uses this procedure.

The GA is a metaheuristic approach also used in combinatorial optimisation problems. It is based on the mechanics of natural selection and natural genetics. GA applications cover a range of combinatorial optimisations, e.g., hydraulic model calibration [61], performance of photovoltaic systems under variable atmospheric conditions [62], and sensor placement for leak detection in water distribution networks [39,63]. The GA is based on three main parameters: selection, crossover and mutation. The population matrix is randomly generated and consists of the design variables, and the best variables are selected according to their fitness value. From these solutions, new solutions are produced via the crossover operator [64]. The mutation operator is finally employed to avoid the algorithm converging to local optima (i.e., to maintain the genetic diversity).

The GA cycle is repeated through a number of generations until a stopping criterion is met. It is worth noting that elitism is not generally considered an operator in the canonical GA. However, it is deemed a robust and effective operator because it leads the optimisation procedure towards the optimal solution. Accordingly, this operator stops the best solutions from being mutated. In this way, the best solutions of each generation would pass to the next, unaltered. Over the course of the algorithm and through a sufficient number of generations, the traits of these solutions would transfer to their offspring, increasing the chance of producing new solutions whose fitness function values might be better than their parents [64]. Some drawbacks of GAs are noted in [61], e.g., achieving a global optimum for large and complex systems is not guaranteed, which is also a drawback for ACO. In [57], the relation between GAs and ACO is noted.

PSO is a recently developed EA that includes features such as easy implementation for solving practical problems, high accuracy and fast convergence of the solution as some of its main advantages [65,66]. While similarities exist in the iterative nature of PSO and GAs, conversely, in PSO, there are no, e.g., "crossover" or "mutation" operations. Instead, PSO is based on a population of candidate solutions, defined as particles. The set of particles composes a swarm, where each individual flows through the parameter space. The flow of such particles is defined by trajectories, which are driven by the best performance of the particle and the neighbouring particles in the parameter space. The initialization of particle swarm is random. The initial solution of each particle represents an alternative solution; that is, each particle has its own initial position and speed and is randomly distributed in each position of the feasible solution space to be searched. Therefore, the initialization of the particle swarm represents the preparation of the particle swarm search. Its size is determined by its speed and position, and the particle update is based on the comparison of the fitness values between each search particle and its neighbouring search particles to determine the necessity of updating a particle. The updated particle adjusts its speed and position according to the particle's new flight path model, which is based on the best results achieved by its neighbouring particles. These conditions yield different optimal experiences for different particle subgroups, which dynamically evolve according to the current position of particles, the particle current velocities, the distance between each particle of the subgroup and its best position and the distance between each particle of the subgroup and the best position of the whole subgroup.

The PSO algorithm does not need cross-mutation or other genetically inspired operations, so the algorithm has fewer parameters and is still high efficiency [65]. These properties are suitable for both engineering applications and scientific research, and a significant number of research results have been produced in recent years [67]. For example, in [68], PSO is applied for function optimisation regarding eco-economics modelling and assessment; in [69], it is used as part of fuzzy systems developed to optimise fuel consumption of hybrid vehicles; and in [70], PSO is used to train neural network models and perform real-time optimisation.

#### *2.2. Centralised Co-Digestion as an Optimisation Problem*

Mathematical optimisation involves the selection of one solution amongst a set, according to some criterion and constraints (that is, the optimisation problem). This optimisation problem can be stated as:

$$\min\_{\mathbf{x}\in\mathcal{X}} f(\mathbf{x}) \text{ subject to } \operatorname{g}(\mathbf{x}), \tag{1}$$

where *f*(*x*) is the objective or cost function, *X* is a feasible region and *g*(*x*) are the constraints that have to hold to find a minimiser *x\** of *f*(*x*) such that *f*(*x*∗) = min*x*∈*X f*(*x*). ACO, GA and PSO introduced in Section 2.1 are algorithms aimed at finding the optimal solution of the optimisation problem posed in (1)—i.e., minimise the objective function *f*(·) subject to the set of constraints *g*(*x*) that apply—which here consists of the selection of the best substrates and volumes according to a set of restrictions related to the operation of the anaerobic digester. In addition, the cost function allows quantifying each alternative potential solution according to (1), involving the calculation of a value or "cost" associated with each alternative considered to find the optimal solution.

The problem statement is similar to that presented in [41], although the number of waste receptors increases from one to three, thus making necessary a reformulation of the optimisation problem involved. Specifically, it is required to increase the dimensions of all data vectors that define each generator-receptor interaction and their subsequent calculations. In addition, matrix operations are repeated per new dimension (i.e., waste receptor) added.

This optimisation problem, which can be understood as a MKP and is of combinatorial nature, can be represented as a matching problem. It is defined with a graph *G =* (*N, E*) that summarises all the possible combinations. The graph consists of *N* vertices (or nodes) and *E* edges or pairs of vertices. Specifically, for the case of AnD optimisation, a bipartite graph can be used to represent the posed optimisation problem to differentiate between the set of waste generators (*N*1), containing *W* nodes, and the set of waste receptors (*N*2), containing *R* nodes. For the specific optimisation problem, all *W* nodes are connected to each *R* node, thus resulting in a total of *W*·*R=E* edges. Figure 1 shows a generic representation of the defined matching problem applied to AnD co-digestion optimisation.

**Figure 1.** Generic representation of the posed matching problem.

A set of substrate generators *w* ∈ {1, . . . , *N*} is considered. The volume of each substrate *Vw* can be selected as a contribution to any of the AnD systems. The binary decision variable *ys w* allows generating array volumetric possibilities (*Vs w*, with *s* ∈ {0..., *lw*} that are determined as a multiple of a number (e.g., 1000 by default) such that 1000*lw* = *Vw*. The selection of each volumetric possibility is determined by the corresponding value of the binary decision variable, *ys w*, where *y* ∈ {0, 1}, with *ys w* = 0 when the corresponding volumetric configuration is not selected, and *ys w* = 1 when it is selected. Note that for each waste generator *w*, there are *lw* different volumetric configurations in *ys w*, but only one is selected at a time, i.e., ∑*lw s*=1 *ys w* = 1 ∀*w* ∈ {1..., *N*}.

To normalise the objective function, selected volumes *Vs w* are divided by the maximum volume from their corresponding waste generator (*Vw*). This approach provides values of the cost function between 0 and 1, where the closest to 1 the better the solution. However, note that the ACO algorithm looks for a maximum of the objective function, while GA and PSO look for a minimum. This behaviour is considered using the constant *K* ∈ {−1, 1}, which depends on the algorithm considered: for the ACO algorithm *K* = 1 and for the GA and PSO algorithms *K* = −1. Hence, cost index *B* would take positive values between 0 and 1 for ACO and negative values between 0 and −1 for GA and PSO. Instead, the absolute value is taken for all three algorithms.

*Fc w* (*c* = 1, ... , 3) and *Tw* are the set of dimensionless coefficients corresponding to the substrate characterisation and the quality term ∑3 *c*=1 *Fc w ρq*, already used and explained previously in [41,46] and defined as shown in Figure 2.

**Figure 2.** (**A**) *F*1 *w*, (**B**) *F*2 *w*, (**C**) *F*3 *w* and (**D**) *Tw* equations used for dimensionless coefficient calculation.

*F*1 *w* is a coefficient related to the potential biogas production, measured as a function of the Chemical Oxygen Demand (*COD*) content. *F*2 *w* indicates the ratio of *COD*/*TN* (where *TN* refers to Total Nitrogen), a useful measure to prevent situations of acidification and other undesired reactions of the AnD process, as long as it is maintained around the range of 20–60. *F*3 *w* is linked to the alkalinity (Alk) concentration, and it is associated with a restriction ranging from 2500 to 6000 mg CaCO3/L integrated within all optimisation algorithms. *Tw* is a coefficient of the utmost importance since it describes the toxicity level of all waste fluxes, which should be kept at the lowest level possible (specifically below 2.1 mg Pb/L).

The *N* different substrate generators are located at different distances (*dw*) from each anaerobic digester. The conveyance of the selected volumes implies a travel distance *dw* (in km) with an economic cost *xw* (in €/km) and a social impact *Iw* = 1, ... , 3 (dimensionless). The higher the value of *Iw*, the higher the social impact of the related route (e.g., proximity to sensitive areas due to pollution, traffic density, or pedestrian presence). Since each route is different for each generator, different values are assigned to approach the logistic impact of the corresponding waste generator, so a value for *Iw* is assigned for each sludge/substrate generator depending on its route to the ST.

The coefficient weight *ρq* (dimensionless) is related to the quality term ∑3 *c*=1 *Fc w ρq*, and the coefficient *ρx* (dimensionless) is the coefficient that weights the logistics term *ρx Xwdw Iw* . Each weight is given a value of 0.5 to provide a balance between the quality and logistics terms in the optimisation. Selected volumes of each substrate to each receptor contribute to the input to the AnD network, and the aforementioned parameters constitute the objective or cost function *f*(*x*).

Additionally, the optimisation problem presented considers a set of restrictions *g*(*x*) related to each of the total inputs to each of the receptor systems, based on those presented in [41,48]. The first restriction is the sum of accepted substrates ∑*N w*=1 ∑*lw s*=0 *ys wVs w* must not exceed the maximum acceptable volume *V* for each AnD system. Moreover, the *COD*/*TN* ratio, related to the dimensionless coefficient *F*2 *w*, must be kept within the range *C*2 *min*, *C*2 *max* . The alkalinity concentration, related to the dimensionless coefficient *F*3 *w*, must also be kept within the range *C*3 *min*, *C*3 *max* . The toxicity level, related to *Tw*, does not require restriction since the corresponding coefficient *Tw* is considered to be restrictive enough, as shown in Figure 2.

In addition, an estimation of the produced biogas is made assuming a conversion factor of 0.268 m3 biogas/kg *COD*. Finally, note that the cost function presented in this work is adapted from [41,48], where the ACO algorithm was used for waste management optimisation in a similar fashion but limited to one AnD receptor.

The objective function *f*(*x*) for the presented optimisation problem is as follows in (2). However, note that the performance comparison of the ACO, GA and PSO algorithms is not conducted directly on the value of the optimised objective function, *B* , but on its absolute value, *B*, as shown in (3).

$$B' = K \left\{ \sum\_{w=1}^{N} \sum\_{s=0}^{l\_w} y\_w^s \frac{V\_w^s}{V\_w} T\_w \left[ \left( \sum\_{c=1}^3 F\_w^c \right) \rho\_q + \frac{\rho\_x}{X\_w d\_w I\_w} \right] \right\} \tag{2}$$

$$B = \begin{vmatrix} B' \end{vmatrix} \tag{3}$$

#### **3. Results**

*3.1. Case Study*

The case study includes a network of 19 organic waste generators and three organic waste receptors. These 22 locations (i.e., 19 generators and 3 receptors) are part of the wastewater treatment system managed by Consorci Besòs Tordera (CBT), a public local water administration composed of 64 municipalities in four different regions of Catalonia (Spain) with a population of approximately 470,000 inhabitants. This case study and its anaerobic network system were also considered in [41]. Figure 3 shows the corresponding bipartite graph of the case study.

**Figure 3.** Bipartite graph of the case study.

The three organic waste receptors (R1–R3, or nodes 1–3 of Figure 3) refer to three separate WWTPs that produce their own sewage sludge, but that also have available AnD technology. Due to oversized design, which is a usual practice in WWTP design [71], these AnD systems in R1–R3 have available capacity. This free excess capacity can be used to accept wastes from external sources, such as the undigested sewage sludge of W1–W12 or the industrial substrates from C1–C7.

The 19 waste generators consist of 12 WWTPs that produce undigested sewage sludge (W1–W12, or nodes 4–15 of Figure 3) and seven industrial substrate generators (C1–C7, or nodes 16–22 of Figure 3), which were considered suitable sources of organic waste for the AnD network under study. Each of these locations is a separate and independent system that must manage its own waste produced as best as possible. Additionally, seven industrial substrate generators have been previously verified as feasible substrates for AnD by CBT technical services.

#### *3.2. Simulation Methodology*

The algorithms used in this work have been implemented in the MATLAB environment. Simulations were performed with a Lenovo ThinkPad (Lenovo Group, Ltd., Girona, Spain) L14 Gen1-20U10016SP ×64 using the OS Microsoft Windows 10 Pro and an Intel(R) Core(TM) i7-10510U CPU processor (1.80 GHz, 2304 MHz) consisting of four main processors and eight logic processors.

The main optimisation parameters of both GA and PSO algorithms were trimmed in an attempt to select the most suitable array to provide reliable results. Accordingly, the same procedure was already performed for the ACO algorithm to determine the values of its corresponding optimisation parameters in [48], where the same optimisation parameters are used in this work.

For the ACO algorithm, an initial population of 100 individuals (or ants) and 500 iterations per repetition is set, and the values used for the algorithm optimisation parameters are *α* = 1, *β* = 2 and *ρ* = 0.98, each corresponding to the importance assigned to the pheromone trail, the importance assigned to the heuristic information and the persistence degree or pheromone evaporation, as explained in [38,58,60]. For GA, the initial population is set to 100, the total number of iterations (or generations) is set to 500, the crossover fraction is set to 0.8 and the fraction of elite children is set to 5% of the corresponding total children. For PSO, the initial population was set to 100, the total number of iterations was set to 500, cognitive attraction was set to 0.8 and the social attraction factor is set to 1.25. Tables 1 and 2 summarise trimming tests for GA and PSO, respectively, where the best results are obtained for higher values of objective index *B*.

**Table 1.** Summary of trimming tests for the GA.

**Table 2.** Summary of trimming tests for PSO.

For the sake of performance comparison, some parameters were fixed for the three algorithms. The fixed parameters are the number of independent simulations (set to 10, the best result is selected), the population (set to 100 individuals), and the maximum number of iterations (set to 500). With these constraints on algorithm trimming, a performance comparison of ACO, GA and PSO was conducted.

The comparison of ACO, GA and PSO performances is based on the value of the fitness function, execution time and an array of technical variables related to the total expected performance of the optimised AnD network: total daily biogas production (in Nm3), average organic load (in kg of COD per m3 of volume of the digestion system and day), average carbon to nitrogen ratio (C/N), and average alkalinity (in mg of CaCO3). All algorithms are tested with data from a real case study as the main simulation scenario. However, other synthetic scenarios are tested to further compare the performance of each algorithm under different scenario conditions.

In the approach presented here, simulated scenarios are based on the waste generator data in Table 3, alongside route distance and receptor system characterisation. For all the 19 waste generators (i.e., the 12 WWTPs without AnD and the seven substrate generators), the addition to the AnD network is optimised. For each of the three AnD systems (i.e., receptors R1, R2 and R3), different volume constraints have been determined, according to operational data and assuming a limit to the hydraulic retention time of 20 days (below that retention time, AnD efficiency is expected to greatly decrease).

**Table 3.** Waste generator dataset, including distance between waste generators and receptors and characterisation of each receptor of the case study (Baseline Scenario or Scenario 0).

Each simulation for ACO, GA and PSO is repeated 10 times since these algorithms have probabilistic, iterative-based search methods. The best solution among these runs is selected for further analysis, although the average fitness function is also registered for discussion.

The data obtained for ACO, GA and PSO (each comprising 10 repetitions of the corresponding algorithm) are compared for every simulated scenario. The baseline scenario (i.e., Scenario 0) corresponds to the real case study, as described in Table 3. Additional synthetic Scenarios 1–4 are simulated, and their corresponding data are created from alterations of the baseline scenario, as described in Table 4.

**Table 4.** Synthetic scenarios created from original Scenario 0 in Table 3. Description of data alteration procedure.

While for Scenarios 1 and 2 any modification is viable, for Scenario 3, an increase in volume involves a significant increase in execution time. This is because the optimisation problem works around combinations of fixed volumes, and an increase in volume would involve a higher number of possible combinations for the algorithms to consider (i.e., an increase in the search space), hence the expected increase in execution time. Thus, volume modification for Scenario 3 was limited to a triple increase in the baseline scenario volume. Alternatively, for Scenario 4, the C/N ratios were modified while being kept below 60 to facilitate the algorithms in finding a viable solution. This measure was adopted because the C/N ratio was the most limiting optimisation parameter in previous applications of a similar optimisation problem in [41]. Scenario 2 was designed with both linear and nonlinear distance modifications (Scenarios 2a and 2b, respectively) to discuss the effect of distance distribution, as pointed out in [41]. Note that trimming tests were carried out only for GA and PSO using the baseline scenario, assuming that trimmed parameters would suffice for simulation of other synthetic scenarios similar to the baseline scenario.

The optimisation results are presented as a sequence of contributions from all the generators to each anaerobic digester. This optimised contribution sequence can be considered a suggested logistic plan for the co-substrate distribution as follows: once enough substrate has been produced and stocked on a waste generator, a truck of 20 metric tonnes capacity would be fully loaded with substrate from the corresponding waste generator, disregarding the truck waiting time before starting each route; once fully loaded, the truck is assumed to travel to the waste receptor without further stops (assuming it always follows the same route). As long as the cycle of supply routes of all involved waste generators is completed within the AnD retention time of 20 days, the properties of the resulting blending should not vary significantly, especially considering that every waste receptor would have a receiving system for these external organic substrates, where they would be stored and blended before being added to the AnD system. The specific start and finish time for each route along the day have not been considered; this does not affect the optimisation, although it has been noted that it has considerable impact on real-world implementation.

#### *3.3. Algorithm Performance Comparison and Scenario Analysis*

The simulation results for each scenario are shown in Table 5. For every scenario and for ACO, GA and PSO, this table shows the best cost index (B) achieved, elapsed optimisation time, and additional parameters related to the performance of the AnD systems: total biogas production, average organic load, average carbon/nitrogen ratio and average alkalinity.

In the baseline scenario, ACO and GA show higher biogas production than PSO (23% and 30% higher, respectively). However, they show a slightly lower *B* index achieved (4% and 11% lower). This result indicates that although one of the main goals of AnD optimisation involves maximising biogas production, it is not all that matters because there are other parameters also subjected to optimisation. PSO appears to find a solution with lower biogas production but better optimises other quality parameters, such as the C/N ratio and alkalinity. However, it is remarkable that out of the three algorithms, ACO and GA find a "similar solution" (prioritising high biogas production), and PSO fins a significantly different solution (prioritising other quality-related parameters).

**Table 5.** Summary of algorithm performance. The best value *B* is highlighted. Scenario 2a feasible results (\*) are associated with a poor solution, so no direct comparison is conducted.

In Scenario 1, the COD concentration was increased tenfold. This was done to compare the efficiency of algorithms to optimise substrates with high organic loads, which is especially meaningful for maximisation of biogas production. For this scenario, ACO and GA have better performance, according to the best index *B* achieved. As a natural consequence of substantial COD increases, biogas production also dramatically increases. However, PSO is unable to achieve a competitive solution in relation to both ACO and GA within this scenario and the baseline scenario, respectively.

In Scenario 2a, a tenfold lineal increase in the geographical distances between facilities was conducted. This scenario allows comparing how well each algorithm can handle situations where most substrates have long distances. For this scenario, ACO is unable to find a feasible solution. On the other hand, GA and PSO find a solution, but the corresponding biogas production is far lower than that obtained in the baseline scenario (48% and 41% lower biogas production for GA and PSO, respectively).

Alternatively, Scenario 2b shows the optimisation results when nonlinearly modifying geographical distances between facilities by the square root of the original distances. This modification allows understanding which algorithm would be more favoured by a more equally distributed geographic location of plants. All the algorithms tested are able to find a feasible solution, showing that GA has the best performance (both in terms of best *B* and biogas production). On the other hand, ACO shows the worst best index *B* achieved.

Scenario 3 was modified by a threefold increase in available volume from all sources. The presented modification allows studying the performance of each algorithm when the total number of possible solutions is much greater. ACO shows noticeably poor performance, below the best index *B* achieved by ACO in former scenarios. Although GA shows better performance than PSO in terms of the best index *B* achieved, biogas production appears similar to that in other scenarios.

In Scenario 4, an increase in the C/N ratio for waste generators W1–W12 was conducted. This modification would test the ability of each algorithm when one of the restrictions (i.e., C/N ratio) requires more adjustments. In this case, PSO shows the best index *B* achieved, although it presents the lowest biogas production. Similar to the baseline scenario, ACO shows slightly better performance than GA in terms of the best *B* achieved, but GA still has slightly better biogas production.

Geographical distance modification was performed with two alternative scenarios. Scenario 2a includes a lineal modification of the distance matrix (tenfold), and Scenario 2b considers a distance modified by the square root of the original distance. The relative locations of all involved waste generators and receptors in the case study are shown in Figure 4. For the baseline scenario, waste generators are homogeneously geographically distributed, but receptors are located in a relatively small area—i.e., the geographical distance difference of each receptor from the emitters might be negligible by the optimisation—which may be interpreted as a single receptor with higher volume capacity, caused by the geographical overlapping of waste receptors, or the "big dot" effect. The linear modification of geographical distances in Scenario 2a does not alter this relative distribution, but the nonlinear modification in Scenario 2b does so, avoiding this "big dot" effect by dispersing Receptors A, B, and C in the geographical space. It is important to note that for Scenario 2a, ACO was unable to find a viable solution, and both GA and PSO achieved a relatively poor solution compared to the corresponding solutions for the baseline scenario.

**Figure 4.** Map of waste generators and waste receptors R1–R3 for Baseline Scenario (**A**), Scenario 2a (**B**) and Scenario 2b (**C**). Distance is expressed as longitudinal distance (X-axis) and latitudinal distance (*Y*-axis) with respect to the R1 plant.

Figure 5 shows the resulting blending profile for the baseline scenario and Scenario 2b. For the baseline scenario, PSO tends to balance the blending of substrates with a low organic load content—i.e., from W1–W12—and selects noticeably lower amounts of high organic load substrates—i.e., from C1–C7— than ACO or GA. This observed behaviour is similar between the three waste receptors A, B and C. On the other hand, ACO and GA tend towards selective blending, showing similar preferences for both receptors B and C. For receptor A, the GA algorithm tends towards slightly more homogeneous blending. In any case, both ACO and GA include more substrates of high organic load—i.e., from C1–C7—except for receptor C.

For Scenario 2b, the ACO blending profiles are similar to those obtained in the baseline scenario—showing a certain tendency to include particular waste generators, —although varying the substrates that the algorithm selects. The GA blending profiles obtained in Scenario 2b are the most affected by geographical distance distortion. GA appears to balance the blending of substrates from all waste generators, much like PSO for both

scenarios 2a and 2b. Additionally, the GA blending profile for Scenario 2b accounts for more industrial, high organic load wastes—i.e., from C1–C7—than the PSO blending profile, which remains relatively similar between the baseline scenario and Scenario 2b.

**Figure 5.** Blending profiles for every waste receptor and ACO, GA and PSO algorithms for the baseline scenario (**left**) and Scenario 2b (**right**).

#### **4. Discussion**

Simulations with the optimisation algorithms ACO, GA and PSO were performed, showing successful optimisation results in almost all scenarios. Data from a real case study were used to carry simulations of centralised anaerobic co-digestion blending. As detailed in Section 3, these datasets are composed of 19 organic waste generators and three organic waste receptors within the context of a sanitation network in an area of high industrial activity in Catalonia. This case study composes the baseline scenario. In that previous work, the potential impacts of optimising AnD with wastes from external sources were already demonstrated, bearing up to 77% cost savings regarding waste management. Different modifications were made to this dataset to compare the performance of the ACO, GA, and PSO algorithms under different conditions to assess the performance of each optimisation algorithm in relevant situations. Regarding the optimisation problem, the C/N ratio is the dominant restriction, as was previously seen in [41]. This is the reason why this parameter is included in the discussion of the results, together with biogas production.

For the baseline scenario as seen in Table 5, PSO shows the best index *B* achieved but also the lowest biogas production. However, PSO also shows the lowest C/N ratio, which might play a role in achieving the best solution, compensating for the lack of biogas produced. If biogas production is increased by the design of a particular setup, this could be trimmed by the corresponding weight in the objective function *B* as a trade-off among the different parameters involved. The results obtained have been considered convenient for the installation under study and improved dramatically performance obtained in the baseline scenario [41]. However, both ACO and GA generally show higher amounts of biogas production, but their best index *B* values achieved are below that of PSO, and their C/N ratios are above 50.

As shown in Figure 5, ACO and GA show similar behaviours for the baseline scenario, prioritising specific substrates. A first hypothesis suggests that prioritised substrates would be those with higher COD since they would allow higher biogas production. On the other hand, PSO shows a different strategy blending more available substrates and tends to exclude industrial substrates. This trend may point to PSO performing a conservative strategy where it is avoided in all cost situations where the operation of the AnD would be put at risk. Therefore, the general trend is that ACO and GA solve the presented optimisation problem by maximising biogas production and pushing restrictions to the limit, while PSO tends to balance biogas maximisation and the C/N ratio trade-off. In addition, note that PSO has the shortest execution times and ACO the largest, which is observed for all scenarios, indicating PSO to be more computationally efficient, where even here, the execution time is not a drawback for real implementation with the values obtained.

The similarities between ACO and GA and the differences between those and PSO could be partially explained by the nature of these algorithms. Both ACO and GA tend to explore the search space of solutions around the borders, thus increasing the number of nonfeasible solutions but also increasing the chances of finding a "rare" solution with a higher best index [72,73]. Thus, these algorithms appear to be based on relatively independent behaviour between particles so that each one can explore separate areas of the border search space and be able to find different non-redundant solutions. On the other hand, PSO algorithm exploration of the search space is based on dependent behaviour between neighbouring particles, which does not encourage particles to explore the limits of the search space. Instead, it promotes the exploration of other mid-term areas between the centre and the borders of the search space. This could help explain why the PSO algorithm attains solutions within shorter execution times but also with generally lower biogas production. Hence, PSO would tend to be a conservative strategy where instead of selecting the most promising solution, single ant or particle, it would prioritise a consensus between the best neighbourhoods.

As detailed in Section 3, Scenario 1 is modified by increasing the organic load of all substrates tenfold. Thus, the dominant condition, in this case, is that organic waste valorisation is fostered, leading to higher biogas production. As observed in Table 5, ACO and GA show better performance than PSO in this scenario, but GA is more efficient since its attained biogas production is noticeably higher than that achieved with ACO.

Additionally, as detailed in Section 3, Scenarios 2a and 2b include a geographical location modification of the involved facilities. As observed in Figure 4, the relative distances between waste generators and receptors (R1, R2, R3) are not modified by lineal modification of the distance when the map plot of the baseline scenario is compared to that of Scenario 2a. However, nonlinear modification of geographical distances in Scenario 2b leads to a different map plot, where waste receptors are more dispersed between them in relation to waste generators (Figure 4). The effect of this distortion of distances is that waste receptors are more separated, thus avoiding geographical overlapping of waste receptors, or the "big dot" effect, i.e., assimilating closer plants as a single centralised plant from a geographical perspective.

Hence, Scenario 2a is modified by tenfold increasing the geographical distances between waste generators and receptors, making geographical distance a dominant condition for optimisation. In this case, ACO is unable to find a solution, and both GA and PSO show extremely poor performance when compared with the baseline scenario, as shown in

Table 5. The linear modification of distances of Scenario 2a shows the performance of each algorithm under the pressure of cases with high geographical distances. This pressure case of Scenario 2a was especially relevant to test because it can significantly impact the logistics processes. On the other hand, Scenario 2b presents a different trend due to the nonlinear modification of distances. As detailed in Table 5, GA shows the best performance, and PSO shows better performance than ACO even in terms of biogas production.

As observed in Figure 5, ACO and PSO maintain similar blending profiles, while GA and PSO also exhibit similar blending profiles, but including GA results in a greater volume of industrial substrates. This shared behaviour between GA and PSO is exclusive to Scenario 2b, but it might indicate that GA behaves similarly to PSO in this case. However, from the operational point of view, GA solutions involve major risks since they tend to include more industrial substrate than PSO solutions.

In Scenario 3, a threefold increase in the volume of all waste generated was performed. First, this result implies that the search space—i.e., the total number of combinations and possible solutions—drastically increases. In this case, Table 5 confirms a similar trend observed for previous scenarios, where GA obtained the best performance and ACO the worst. Again, this finding is consistent with the observation that the ACO algorithm attains weaker performance than GA and PSO for this particular case and that GA and PSO attain similar performance in this study, although GA appears to generally provide better performance than PSO.

Finally, Scenario 4 was composed of increasing the C/N ratio of W1–W12 substrates. These substrates originally conformed to sewage sludge with a low nitrogen load, but in Scenario 4, the drastic increase in the C/N content of sewage sludge was the dominant condition to be tested. Table 5 also presents a summary of the results for Scenario 4, where PSO shows the best performance and GA the worst. The main observation is that the PSO algorithm is more able than the GA to manage situations with high nitrogen loads or major restrictions, while the GA has more potential to maximise biogas production. However, it is more sensitive to high nitrogen loads because it reduces available space to acquire industrial wastes with both high organic loads and high nitrogen loads.

The developed algorithms have successfully optimised the AnD network of the case study, and their performances have been tested under different conditions (i.e., Scenarios 1–4). Simultaneous logistics and quality optimisation of a network of existing waste management facilities is a gap in the current state of the art due to its ad-hoc nature and its interdisciplinarity: there are specialized works for logistics optimisation such as in [45,46], but they do not include process optimisation. The present study implements this logistic optimisation by minimising a cost function designed to this end. The reason for choosing this approach is also based on the need for professionals who manage the AnD network considered here to have a decision support tool capable of integrating logistics and process performance optimisation. And, in addition, to have the ability to handle changing operation conditions and scenarios, as it actually happens in real facilities.

It is also worth noting that it exists a variety of sensors for the determination of physical-chemical parameters that could complement the sensor network considered in this case study, such as a variation of the ones presented in [47]. These sensors could provide additional insight, especially if combined with GIS and process optimisation, and also facilitate the real-time implementation of the presented approach. Additionally, they could also be used as control mechanisms for those cases where ACO and GA optimisation is applied, since attained optimised outcomes pushed quality restrictions of the AnD process close to their thresholds. However, there is a trade-off between information (i.e., data gathered from new sensors) and resources (e.g., implementation, maintenance) which has to be taken into account when considering new sensors.