Optimal Number of Pressure Sensors for Real-Time Monitoring of Distribution Networks by Using the Hypervolume Indicator

: This article proposes a novel methodology to determine the optimal number of pressure sensors for the real-time monitoring of water distribution networks based on a quality hypervolume indicator. The proposed methodology solves the optimization problem for different numbers of pressure sensors, assesses the gain of installing each set of sensors by means of the hypervolume indicator and determines the optimal number of sensors by the variation of the hypervolume indicator. The methodology was applied to a real case study. Several robustness analyses were carried out. The results demonstrate that the methodology is hardly inﬂuenced by the method parameters and that a reasonable estimation of the optimal number of sensors can be easily achieved. optimization method based on the maximization of nodal sensitivities to both pipe burst events and variations of the pipe roughness coefﬁcient. The obtained results are discussed and the most relevant conclusions are drawn.


Introduction
Water distribution networks (WDN) constitute vital infrastructure to ensure drinking water in sufficient quantity and quality for users. Real-time monitoring of these networks is essential to ensure the quality of the service provided. This is achieved by employing on-site sensors to measure certain water quality (e.g., pH, temperature, contaminants) and hydraulic (i.e., pressure and flowrate) parameters [1]. These sensors produce different types of data (e.g., flowrate, pressure, contaminants concentration) that can be used in different activities. Only a limited number of sensors can be installed in a given WDN due to budget constraints. Moreover, these WDN may rapidly increase in complexity (as a function of the area served), resulting in a challenging task of real-time monitoring of large networks when few sensors are installed [2].
Numerous methods have been developed over the recent decades to find the optimal location for sensors to detect contamination events in real time and mitigate the associated effects. In the Battle of the Water Sensor Networks [3], fifteen different approaches of optimal sensor location for contaminant detection were compared according to four different objectives. Aral et al. [4] blended four distinct criteria into one single objective which was solved using genetic algorithms. According to Weickgenannt et al. [5], the risk of contamination is explicitly evaluated as the product of likelihood of not detecting the contaminant intrusion and the corresponding consequence (water consumed). In a different approach, Zhao et al. [6] attempted to minimize the consumption of contaminated water prior to contamination detection. More recently, Naserizade et al. [7] used multi-objective optimization and a multicriteria decision-making technique whilst Ponti et al. [8] proposed a new multi-objective evolutionary algorithm which showed improvements over NSGA-II. However, sensor location techniques aiming at monitoring water quality parameters may not provide the best solutions for hydraulic parameters since they have different goals in the utility activities.
In particular, pressure data are widely used by water utilities in the real-time monitoring of WDN for creating alerts and use in hydraulic modelling calibration and pipe burst detection. Alarm systems can be designed to respond to the occurrence of abnormal pressure values (low or high pressure) at specific locations of the network [9]. The calibration process of hydraulic models requires pressure data from specific locations of the network to calibrate high uncertainty variables such as nodal demands and pipe roughness coefficient [2,10]. Pipe burst detection and location techniques use pressure data in a variety of ways, for instance, to carry out inverse analysis of the pipe burst location (using hydraulic simulation and optimization) [11,12] or by using a data-driven classifier approach [13].
Distinct methods have been proposed in the literature to optimally locate pressure sensors in WDN, most of them focusing on the total gain of sensor placement considering different aims. Cao et al. [14] positioned pressure sensors in such a way as to represent the pressure patterns of the homogeneous area of the WDN. Zhao et al. [15] aimed at maximizing the detection coverage rate of pipe burst events whilst Raei et al. [16] aimed at reducing the detection time of these events. Both Casillas et al. [17] and Steffelbauer et al. [18] focused on maximizing the percentage of leak scenarios correctly identified (according to the introduced criteria). Sarrate et al. [19] focused on maximizing the leak detectability performance. Furthermore, the maximization of nodal pressure sensitivities was considered both by de Schaetzen et al. [20] and Francés-Chust et al. [21]. The maximization of the accuracy of the hydraulic model was also considered by Kapelan et al. [22] and Behzadian et al. [23].
Ideally, the optimal location and number of pressure sensors should be determined simultaneously. Nonetheless, most of the previously presented methods are not able to directly provide the optimal number of pressure sensors to be installed. As a result, the number of pressure sensors to be installed in a WDN is usually defined using different criteria. Sarrate et al. [19] and Sanz et al. [24] defined the number of pressure sensors for pipe burst detection by budget limitations whilst both Sophocleous et al. [11] and Wu et al. [25] used a metric of the number of households per pipe unit length. For the same objective of pipe burst detection, Soldevila et al. [26] used a metric based on the network size, and Quintiliani et al. [27] used engineering good sense. Regarding pressure sensors for model calibration, de Schaetzen et al. [20] used the number of households per unit length to decide the number of sensors whilst Wéber et al. [2] emphasized the minimum percentage of calibration errors to be achieved.
This article presents a novel methodology to determine the optimal number of pressure sensors in a WDN. The major advantage is the possibility to be coupled with the existing optimization methods of sensor location (such as those previously presented), thus providing the optimal number according to the considered objective functions. Furthermore, the methodology does not require previous knowledge on sensor installation costs nor the consideration of costs in the optimization problem itself. The methodology uses the hypervolume indicator [28] to measure the total gain of installing different numbers of pressure sensors (according to the objective functions being considered) with a single value. In order to reduce the number of solved optimization problems, a discrete set of numbers of sensors was established in a distributed and representative way for which the optimization problem would be solved. A trade-off function was then fitted to the hypervolume data, which can be interpreted as a compromise between the total gain (using the hypervolume) and the number of sensors. Finally, the optimal number of pressure sensors was determined by assessing the evolution of the hypervolume indicator as a function of the number of sensors, specifically by identifying the point of maximum curvature of the hypervolume data. From that point onwards, the inclusion of another sensor was rather not worth the increment on the total gain. This methodology is demonstrated through application to a real WDN. The pressure sensors were located using an optimization method based on the maximization of nodal sensitivities to both pipe burst events and variations of the pipe roughness coefficient. The obtained results are discussed and the most relevant conclusions are drawn.
The main novel contributions are as follows: (i) the proposal and application of the hypervolume indicator to quantify the total gain associated with a given number of pressure sensors and (ii) the application of different techniques to detect the point of maximum curvature aiming to determine the optimal number of pressure sensors.

General Approach
A flowchart of the proposed methodology is present in Figure 1 and is composed of six main steps: (1) definition of the maximum number of sensors; (2) establishment of a discrete set of numbers of sensors; (3) optimization of the pressure sensor locations; (4) assessment of the total gain of sensors; (5) trade-off function calculation; (6) determination of the optimal number of pressure sensors. The following sections further detail each step of the proposed methodology.
Water 2021, 13, x FOR PEER REVIEW This methodology is demonstrated through application to a real WDN. The pre sensors were located using an optimization method based on the maximization of sensitivities to both pipe burst events and variations of the pipe roughness coefficien obtained results are discussed and the most relevant conclusions are drawn.
The main novel contributions are as follows: (i) the proposal and application hypervolume indicator to quantify the total gain associated with a given number of sure sensors and (ii) the application of different techniques to detect the point of maxi curvature aiming to determine the optimal number of pressure sensors.

General Approach
A flowchart of the proposed methodology is present in Figure 1 and is compos six main steps: (1) definition of the maximum number of sensors; (2) establishmen discrete set of numbers of sensors; (3) optimization of the pressure sensor location assessment of the total gain of sensors; (5) trade-off function calculation; (6) determin of the optimal number of pressure sensors. The following sections further detail each of the proposed methodology.

Definition of the Maximum Number of Sensors
The optimal number of pressure sensors is determined by analyzing the variat the total gain (using the results of optimization problems) with the number of sen However, solving optimization problems for an incremental number of sensors (e.g tween 1 and 150) can be computationally expensive due to the exponential increase search space. For example, considering 20, 21 and 22 sensors in a total set of 150 po locations results in 3.63 × 10 24 , 2.24 × 10 25 and 1.31 × 10 26 possible combinations of se respectively. Therefore, the maximum number of installed sensors, Nmax, is firstly de aiming at minimizing the computational burden by limiting the number of optimiz problems to be solved. Different approaches can be used for the definition of the maximum number of sure sensors in a WDN, for instance, by using economic criteria based on the size WDN [18] or by determining when the gain of installation of an additional sensor h significant improvement [15], taking into account the total number of possible loca [16] or the available budget [23,27].

Definition of the Maximum Number of Sensors
The optimal number of pressure sensors is determined by analyzing the variation of the total gain (using the results of optimization problems) with the number of sensors. However, solving optimization problems for an incremental number of sensors (e.g., between 1 and 150) can be computationally expensive due to the exponential increase in the search space. For example, considering 20, 21 and 22 sensors in a total set of 150 possible locations results in 3.63 × 10 24 , 2.24 × 10 25 and 1.31 × 10 26 possible combinations of sensors, respectively. Therefore, the maximum number of installed sensors, N max , is firstly defined, aiming at minimizing the computational burden by limiting the number of optimization problems to be solved. Different approaches can be used for the definition of the maximum number of pressure sensors in a WDN, for instance, by using economic criteria based on the size of the WDN [18] or by determining when the gain of installation of an additional sensor has no significant improvement [15], taking into account the total number of possible locations [16] or the available budget [23,27].

Establishment of the Discrete Set of Numbers of Sensors
A discrete set of numbers of sensors is established between 1 and the maximum number of sensors N max , for which the optimization problem will be solved. The established set does not have to be a continuous sequence of natural numbers, for instance, it could be 1, 5, 10, 15, . . . , N max . The aim is to reduce the number of optimization problems to be solved whilst leaving enough results to characterize the total gain as a function of the number of sensors between 1 and N max . As such, different sets of the number of sensors can be considered, for instance, with an evenly distributed number of sensors between 1 and N max or with a higher density of observations in the lower number of sensors. The effect of different discrete sets of the numbers of sensors is assessed in this article.

Optimization of the Pressure Sensor Locations
The total gain with the installation of pressure sensors depends on the method used for locating those sensors. Different pressure sensor location methods can be coupled to the proposed methodology, based on different optimization problem formulations [2,[15][16][17][18][19][20][21][22][23]29].
A method to optimally locate pressure sensors based on the maximization of nodal sensitivities is considered to demonstrate the performance of the proposed methodology. This method is based on the formulation of an unconstrained multi-objective optimization problem, in which the decision variables are the pressure sensor locations (using nodes as possible locations). It requires the prior computation of two pressure sensitivity matrices S1 and S2 (one for the pipe roughness coefficient and the other for the pipe burst size, respectively). The size of the obtained matrices is the number of pipes × the number of nodes and the number of nodes × the number of nodes, respectively. As such, S1 i,j refers to the variation of the pressure in node j with respect to the variation of the pipe roughness coefficient of pipe i. Similarly, S2 i,j refers to the variation of the pressure in node j with respect to the variation in the pipe burst size in node i. Further details regarding pressure sensitivity analysis can be found in the works by de Schaetzen et al. [20] and Lansey et al. [30].
Accordingly, two objective functions are formulated for a given sampling design X of sensor locations. The first function f 1 is defined according to de Schaetzen et al. [20] by using a compromise programming formulation. It aims at maximizing the sensitivity to the pipe roughness coefficient that the sensors cover (f A ) whilst ensuring the evenly spread geographical distribution of sensor locations by maximizing the entropy according to Shannon's definition (f B ). Functions f A and f B are defined as follows: where Np is the number of pipes. Finally, the two functions f A and f B are combined together into f 1 by taking the weighted sum of the two functions.
The second function f 2 aims at directly maximizing the sensitivity of sensors to pipe burst events and can be computed as follows: where Nn is the number of nodes. Further details regarding the formulation of the objective functions can be found in the work by de Schaetzen et al. [20].
The optimization problem should be solved for each number of sensors in the discrete set. The results of each optimization problem are a single set of pressure sensor locations (for single objective optimization) or a set of optimal combinations of sensors (for multiple objective functions). In the latter case, a Pareto front of solutions is obtained; each solution of a Pareto front refers to a specific combination of sensors for which no other combination exists that presents better results for both objectives. Each optimal combination is characterized by the values associated with each objective function. These values are used to Water 2021, 13, 2235 5 of 16 assess the gain of installing the different number of sensors in the WDN. Note that the same objective functions and optimization method should be used for all numbers of sensors in the discrete set.

Assessment of the Total Gain of Sensors
The result of the previous step is a Pareto front for each number of sensors in the discrete set. Each Pareto front embraces all the possible combinations of sensors for that specific number of sensors.
The assessment of the total gain for a given number of sensors is not a trivial task as the gain can be described by multiple objective functions and, therefore, can have multiple different values; consider the example presented in Figure 2a, where a single optimal configuration of five sensors is depicted with a red dot. The gain of that specific configuration of sensors is characterized by three different values in accordance with some hypothetical objective functions 1, 2 and 3 to maximize (these are deliberately different from f 1 and f 2 since the methodology can be adopted to different optimization problem formulations with different dimensions). Furthermore, multiple optimal configurations might exist for a given number of sensors due to a trade-off between the objectives. This is represented in Figure 2b with a Pareto front of four optimal combinations of five sensors depicted with red dots. Note that the increase of objective function 3 will lead to the decrease of objective functions 1 and 2 as a trade-off between these contradictory objectives.

Trade-Off Function Fitting
The result of the previous step is a set of hypervolume values, with one value associated with each number of sensors in the discrete set previously defined in Section 2.3. As the number of installed pressure sensors increases, the hypervolume also increases since more sensors lead to a higher total gain. This can be interpreted as a trade-off between the total gain of installing pressure sensors (using the hypervolume) and the number of sensors.
Therefore, observations of the hypervolume are used to derive a trade-off function that describes the behavior of the hypervolume (interpreted as the total gain) as a function of the number of sensors. By doing so, it is possible to estimate the total gain associated with the number of sensors for which the optimization problem is not solved.
Different algorithms can be used to fit a function to the hypervolume values by solving the nonlinear least-squares problem. A review on numerical methods for nonlinear least-squares optimization problems can be found in the work by Yuan [43]. The use of the Levenberg-Marquardt method is proposed to optimally fit a function by finding its optimal set of parameters. It is a well-established method with fast convergence, although it may be sensitive to the initial parameters guess.
Distinct trade-off functions (with different mathematical formulations and numbers of parameters) can be fitted to the same hypervolume data. In this study, the well-known root-mean-square error (RMSE) is used to assist in deciding which function (after the fitting) best describes the hypervolume data. The smaller the value of the RMSE is, the better the fitted function represents the behavior of the hypervolume. In order to assess the total gain of installing a different number of pressure sensors, the quality of each obtained Pareto front is assessed using a quality measure [28,[31][32][33][34][35][36]. This calculation leads to a single value for each Pareto front, characterizing the total gain associated with each specific number of sensors.
The use of the hypervolume indicator [28] is proposed herein to measure the quality of each Pareto front. It is one of the most widely used quality measures, the main advantages whereof being its easy interpretation and properties, such as guaranteeing strict monotonicity regarding Pareto dominance [37]. Nonetheless, this indicator also has disadvantages, such as sensitivity to the presence or absence of extreme points in a Pareto front [38]. In sum, this indicator considers the "volume" (i.e., in as many dimensions as the number of the objective functions) of the region of the objective space dominated by the set of optimal solutions in relation to the global worst point. Note that the hypervolume has no defined units, and its magnitude is related to the objective functions used. The hypervolume is graphically depicted in Figure 2 as a grey volume and considering the global worst point as (0,0,0). The formal mathematical definition of the hypervolume can be found in the work by Guerreiro et al. [37].
Different methods of the hypervolume calculation have been developed [37][38][39][40][41][42], with different efficiency levels against high dimensionality problems [37]. The efficiency of the hypervolume calculation in this particular methodology is not a major concern because no high dimensionality is expected and a reduced number of hypervolume calculations is required (specifically, one calculation for each Pareto front).

Trade-Off Function Fitting
The result of the previous step is a set of hypervolume values, with one value associated with each number of sensors in the discrete set previously defined in Section 2.3. As the number of installed pressure sensors increases, the hypervolume also increases since more sensors lead to a higher total gain. This can be interpreted as a trade-off between the total gain of installing pressure sensors (using the hypervolume) and the number of sensors.
Therefore, observations of the hypervolume are used to derive a trade-off function that describes the behavior of the hypervolume (interpreted as the total gain) as a function of the number of sensors. By doing so, it is possible to estimate the total gain associated with the number of sensors for which the optimization problem is not solved.
Different algorithms can be used to fit a function to the hypervolume values by solving the nonlinear least-squares problem. A review on numerical methods for nonlinear leastsquares optimization problems can be found in the work by Yuan [43]. The use of the Levenberg-Marquardt method is proposed to optimally fit a function by finding its optimal set of parameters. It is a well-established method with fast convergence, although it may be sensitive to the initial parameters guess.
Distinct trade-off functions (with different mathematical formulations and numbers of parameters) can be fitted to the same hypervolume data. In this study, the well-known root-mean-square error (RMSE) is used to assist in deciding which function (after the fitting) best describes the hypervolume data. The smaller the value of the RMSE is, the better the fitted function represents the behavior of the hypervolume.
The selected fitted function is used to derive estimations of the hypervolume for the complete set of numbers of sensors between 1 and N max .

Determination of the Optimal Number of Pressure Sensors
The hypervolume estimations are analyzed in order to determine the optimal number of sensors. This optimal number corresponds to the point where the inclusion of another sensor does not significantly improve the increment on the total gain. Such a point is defined as the point of maximum curvature in the hypervolume curve. The mathematical formulation for the point of maximum curvature can be found in the work by Satopaa et al. [44]. This point is defined herein as the "knee/elbow" of the trade-off function and can be identified using different automatic techniques [44][45][46][47]. The two techniques used herein are based on the works of Satopaa et al. [44] and Salvador and Chan [46], respectively. The performance of both techniques is assessed later in the article.
The first technique, the Kneedle method, starts by rescaling the estimated hypervolume data to [0,1], leading to the normalized hypervolume values, HV norm . The number of sensors between 1 and N max should also be rescaled to [0,1], leading to the normalized number of sensors values, N norm . These values are represented in Figure 3 as black circles and squares, respectively, considering N max = 15 in this example. Finally, the differences ∆ between HV norm and N norm are calculated, as exemplified with grey arrows and triangular markers in Figure 3. The optimal number of sensors, N opt , is selected as the number of sensors that presents the maximum difference. and squares, respectively, considering Nmax = 15 in this example. Finally, the differences Δ between HVnorm and Nnorm are calculated, as exemplified with grey arrows and triangular markers in Figure 3. The optimal number of sensors, Nopt, is selected as the number of sensors that presents the maximum difference. The second technique, the L-method, determines the "knee/elbow" as the point that best divides the hypervolume data into two straight lines. Pairs of lines are fitted to the hypervolume data using a linear regression procedure, and their junction point is the looked-for "knee/elbow". Each line must contain at least two points and start at either end of the hypervolume data. Figure 4 shows all the possible pairs of straight lines for an example of the hypervolume data considering Nmax = 7. The seven hypervolume values are presented both as blue or orange circles (either they are on the left or right side of the junction point). The corresponding fitted straight lines are represented with blue or red dashed lines, respectively. The number of sensors tested for the optimal number is in a blue diamond. The optimal number of sensors is the junction point that minimizes the weighted RMSE for the two linear parts of the hypervolume data. The final optimal locations can be found in the Pareto front related to the optimal number of sensors Nopt. Nonetheless, the Pareto front might not have yet been obtained (i.e., when Nopt was not considered in the discrete set). In these cases, the optimization problem should be solved for the number of sensors equal to Nopt. The second technique, the L-method, determines the "knee/elbow" as the point that best divides the hypervolume data into two straight lines. Pairs of lines are fitted to the hypervolume data using a linear regression procedure, and their junction point is the looked-for "knee/elbow". Each line must contain at least two points and start at either end of the hypervolume data. Figure 4 shows all the possible pairs of straight lines for an example of the hypervolume data considering N max = 7. The seven hypervolume values are presented both as blue or orange circles (either they are on the left or right side of the junction point). The corresponding fitted straight lines are represented with blue or red dashed lines, respectively. The number of sensors tested for the optimal number is in a blue diamond. The optimal number of sensors is the junction point that minimizes the weighted RMSE for the two linear parts of the hypervolume data. and squares, respectively, considering Nmax = 15 in this example. Finally, the differences Δ between HVnorm and Nnorm are calculated, as exemplified with grey arrows and triangular markers in Figure 3. The optimal number of sensors, Nopt, is selected as the number of sensors that presents the maximum difference. The second technique, the L-method, determines the "knee/elbow" as the point that best divides the hypervolume data into two straight lines. Pairs of lines are fitted to the hypervolume data using a linear regression procedure, and their junction point is the looked-for "knee/elbow". Each line must contain at least two points and start at either end of the hypervolume data. Figure 4 shows all the possible pairs of straight lines for an example of the hypervolume data considering Nmax = 7. The seven hypervolume values are presented both as blue or orange circles (either they are on the left or right side of the junction point). The corresponding fitted straight lines are represented with blue or red dashed lines, respectively. The number of sensors tested for the optimal number is in a blue diamond. The optimal number of sensors is the junction point that minimizes the weighted RMSE for the two linear parts of the hypervolume data. The final optimal locations can be found in the Pareto front related to the optimal number of sensors Nopt. Nonetheless, the Pareto front might not have yet been obtained (i.e., when Nopt was not considered in the discrete set). In these cases, the optimization problem should be solved for the number of sensors equal to Nopt. The final optimal locations can be found in the Pareto front related to the optimal number of sensors N opt . Nonetheless, the Pareto front might not have yet been obtained (i.e., when N opt was not considered in the discrete set). In these cases, the optimization problem should be solved for the number of sensors equal to N opt .

Case Study
This methodology is demonstrated in a real WDN with a total network extension of ca. 70 km, 2200 service connections and an average inlet flowrate of about 60 m 3 /h. The supplied area consists mainly of single-family houses with irrigated gardens, most of them with a swimming pool. The WDN also supplies a few golf courses. The hydraulic simulation model was developed in EPANET [48] and included four storage tanks, 4474 pipes and 4429 nodes, of which 2200 were consumers with an individual hourly demand pattern. Figure 5 presents the layout of the hydraulic network.
Water 2021, 13, x FOR PEER REVIEW 8 of 16

Case Study
This methodology is demonstrated in a real WDN with a total network extension of ca. 70 km, 2200 service connections and an average inlet flowrate of about 60 m 3 /h. The supplied area consists mainly of single-family houses with irrigated gardens, most of them with a swimming pool. The WDN also supplies a few golf courses. The hydraulic simulation model was developed in EPANET [48] and included four storage tanks, 4474 pipes and 4429 nodes, of which 2200 were consumers with an individual hourly demand pattern. Figure 5 presents the layout of the hydraulic network. The maximum number of installed sensors, Nmax, was determined based on the intrinsic network characteristics, namely, on network length. A maximum of one sensor per kilometer of pipe network was adopted according to several Portuguese expert opinions, leading to Nmax = 70 sensors (i.e., one sensor per km).
Five discrete sets of numbers of sensors were defined. The objective was to assess if major differences in the optimal number of sensors were found when considering distinct sets of numbers of sensors. These five sets presented distinct characteristics and were defined as follows: Note that both Set 1 and Set 2 attempted to describe evenly distributed observations of the hypervolume between 1 and Nmax whilst Set 3, Set 4 and Set 5 incorporated a higher density in the lower number of sensors. Furthermore, the number of observations highly varied between the sets: Set 2 and Set 3 contained only five numbers of sensors (leading to five optimization problems to be solved) whilst Set 5 contained 26 numbers of sensors (leading to a total of 26 optimization problems).
The method described in Section 2.4 was used to optimally locate the pressure sensors. Note that other methods can be used with distinct objective functions (both in number and formulation). Two pressure sensitivity matrices were computed prior to the optimization problem itself. A variation of Hazen-Williams's pipe roughness coefficient of 10 was considered for the computation of the first matrix. The second matrix was obtained The maximum number of installed sensors, N max , was determined based on the intrinsic network characteristics, namely, on network length. A maximum of one sensor per kilometer of pipe network was adopted according to several Portuguese expert opinions, leading to N max = 70 sensors (i.e., one sensor per km).
Five discrete sets of numbers of sensors were defined. The objective was to assess if major differences in the optimal number of sensors were found when considering distinct sets of numbers of sensors. These five sets presented distinct characteristics and were defined as follows: Note that both Set 1 and Set 2 attempted to describe evenly distributed observations of the hypervolume between 1 and N max whilst Set 3, Set 4 and Set 5 incorporated a higher density in the lower number of sensors. Furthermore, the number of observations highly varied between the sets: Set 2 and Set 3 contained only five numbers of sensors (leading to five optimization problems to be solved) whilst Set 5 contained 26 numbers of sensors (leading to a total of 26 optimization problems).
The method described in Section 2.4 was used to optimally locate the pressure sensors. Note that other methods can be used with distinct objective functions (both in number and formulation). Two pressure sensitivity matrices were computed prior to the optimization problem itself. A variation of Hazen-Williams's pipe roughness coefficient of 10 was considered for the computation of the first matrix. The second matrix was obtained by generating a burst of fixed size for every node of the hydraulic model with a single emitter coefficient of 0.25 (leading to an average burst of 5 L/s). The hydraulic simulations were carried in EPANET [48].
The optimal sensor location for a given number of sensors was formulated as an unconstrained multi-objective optimization problem. The decision variables were the nodes where pressure sensors could be potentially installed, and all the nodes were considered as possible locations. The two objective functions presented in Section 2.4 were considered, Water 2021, 13, 2235 9 of 16 aiming at the maximization of nodal pressure sensitivities to both pipe roughness coefficient variations (f 1 ) and pipe burst events (f 2 ).
The problem was solved once for each number of sensors in each of the five discrete sets. When the number of sensors appeared in more than one set (e.g., 70), the problem was only solved once. The NSGA-II algorithm was applied in the Python environment using the Pymoo package [49] to solve several optimization problems. Discrete variables and integer encoding were used for problem formulation; each discrete location (i.e., node) was translated by an integer value. A population member is a set of locations for the pressure sensors and each variable of a population member represents a possible pressure sensor location.
The following NSGA-II parameters were used: random integer sampling and selection operators; integer polynomial mutation with the probability p m = 0.05 and index parameter n m = 20; simulated binary crossover with probability p c = 0.95 and index parameter n c = 20.
A population size of 100 was considered and all NSGA-II runs were carried out for 500 generations. As such, each NSGA-II run (i.e., solving an optimization problem for a given number of sensors in a given set) resulted in 500 × 100 = 50,000 objective function evaluations. The optimization problems were solved using an Intel Core i5-8250U processor of 1.80 GHz and 8 GB of memory, with a total running time of around 1 h per optimization problem. Table 1 presents the number of objective function evaluations associated with each of the five discrete sets; compare the computation effort (i.e., the number of objective function evaluations) in Set 2 and Set 3 (both with five optimization problems to be solved) with that of Set 5 (with 26 optimization problems to be solved).  The result of each optimization problem was a Pareto front of optimal pressure sensor locations, according to the objective functions f 1 and f 2 . For the sake of simplicity, the obtained Pareto fronts for 10, 20, 30 and 40 sensors are depicted in Figure 6. Each circle represents an optimal configuration of either 10, 20, 30 or 40 sensors. The hypervolume indicator was calculated for each Pareto front of each set of numbers of sensors, leading to a total of five sets of Hypervolume values (i.e., one for each set of numbers of sensors). The global worst point was considered equal to (0,0) since both functions aimed at the maximization of the objective function. The calculation method for the hypervolume was based on variant 3 of the algorithm proposed by Fonseca et al. [42]. The five sets of hypervolume values are depicted in Figure 7 as black triangles.

Objective Function Evaluations
Five different functions, Fi, were considered to describe the hypervolume data as a function of the number of sensors, N: The hypervolume indicator was calculated for each Pareto front of each set of numbers of sensors, leading to a total of five sets of Hypervolume values (i.e., one for each set of numbers of sensors). The global worst point was considered equal to (0,0) since both functions aimed at the maximization of the objective function. The calculation method for the hypervolume was based on variant 3 of the algorithm proposed by Fonseca et al. [42]. The five sets of hypervolume values are depicted in Figure 7 as black triangles.  Table 2 presents the five functions' specific parameters found using the Levenberg-Marquardt method for Set 1 of hypervolume values.  The RMSE was calculated for each function Fi of each set of numbers of sensors in order to assist in deciding which function Fi best describes the hypervolume data. The obtained results are presented in Table 3, with the best fitted function for each set of numbers of sensors presented in bold and shaded areas. Finally, the optimal number of sensors Nopt (considered as the point of maximum curvature) was determined using an automatic detection technique. The two distinct techniques presented in Section 2.7 were used for each function Fi of each set of numbers of sensors, with the obtained results presented in Table 3. Figure 8 presents the application of both Kneedle method and L-method techniques to the estimated hypervolume for the F5 function of Set 1.  Five different functions, F i , were considered to describe the hypervolume data as a function of the number of sensors, N: where a, b, c and d are the function parameters. These parameters, once optimized, can represent the curve of the hypervolume as a function of the number of sensors. Note that different functions can be considered to describe the increasing behavior one can see in Figure 7 in black triangles. The Levenberg-Marquardt method was used to fit each function F i to each set of hypervolume values by finding the optimal set of parameters. This was done by using MATLAB's Curve Fitting Toolbox. A reasonable initial parameter guess was found by coarsely gridding the parameter space. This was carried out to cope method sensitivity to the initial parameters. Figure 7 presents in colored lines the estimated hypervolume values obtained by using each model F i . The original hypervolume values are represented as triangular black markers. Table 2 presents the five functions' specific parameters found using the Levenberg-Marquardt method for Set 1 of hypervolume values.  The RMSE was calculated for each function F i of each set of numbers of sensors in order to assist in deciding which function F i best describes the hypervolume data. The obtained results are presented in Table 3, with the best fitted function for each set of numbers of sensors presented in bold and shaded areas. Finally, the optimal number of sensors N opt (considered as the point of maximum curvature) was determined using an automatic detection technique. The two distinct techniques presented in Section 2.7 were used for each function F i of each set of numbers of sensors, with the obtained results presented in Table 3. Figure 8 presents the application of both Kneedle method and L-method techniques to the estimated hypervolume for the F 5 function of Set 1. Table 3. RMSE and the optimal number of sensors N opt for the different fitted functions and discrete sets of sensors. The results for the best-fitted function of each discrete set are presented in bold and shaded area. 10, 20, 30, 40, Table 2 presents the five functions' specific parameters found using the Levenberg-Marquardt method for Set 1 of hypervolume values. The RMSE was calculated for each function Fi of each set of numbers of sensors in order to assist in deciding which function Fi best describes the hypervolume data. The obtained results are presented in Table 3, with the best fitted function for each set of numbers of sensors presented in bold and shaded areas. Finally, the optimal number of sensors Nopt (considered as the point of maximum curvature) was determined using an automatic detection technique. The two distinct techniques presented in Section 2.7 were used for each function Fi of each set of numbers of sensors, with the obtained results presented in Table 3. Figure 8 presents the application of both Kneedle method and L-method techniques to the estimated hypervolume for the F5 function of Set 1.

Discussion
The increase in the number of installed sensors led to a higher total gain. Nonetheless, the marginal gain decreased when considering additional sensors. This is visible in Figure 7 in black triangles. This figure shows that regardless of the considered number of sensors in the set, the function F 1 (grey line) clearly fails to mathematically describe the trend of the hypervolume as a function of the number of sensors, whereas the remaining functions can approximately describe the hypervolume variation. Table 3 presents better insights regarding the fitting capabilities of the different functions F i for each discrete set of numbers of sensors. The best fit corresponds to the smallest values of the RMSE. Complex models F 4 and F 5 (i.e., both with a larger number of parame-ters) presented the best fitting results (in bold) for most sets of numbers of sensors, except for Set 3, for which F 3 led to the lowest RMSE values.
The first robustness analysis was carried out to determine the effect of the different functions F on the optimal number of sensors. For each technique (Kneedle method and Lmethod) of each set of numbers of sensors, it was possible to verify that the optimal number of sensors did not significantly vary with the fitted functions F i to estimate the hypervolume value (note that results from F 1 were excluded from analysis). Technique 1 presented higher variability in the optimal number of sensors, ranging from 16 to 18. Technique 2 presented smaller variability, between 16 and 17. Therefore, and as long as selected functions F i (and respective parameters) can properly represent the hypervolume as a function of the number of sensors, a good estimation of N opt can be obtained (considering either the Kneedle method or the L-method).
The second robustness analysis was carried out to analyze the effect of different discrete sets of numbers of sensors on the optimal number of sensors. For each technique and function F i , it was possible to verify that the optimal number was relatively stable for the different sets, and no differences greater than one were found. Thus, and as long as the selected discrete set of numbers of sensors can properly represent the hypervolume as a function of the number of sensors, a good estimation of N opt can be obtained. This allows for a significant reduction in the number of optimization problems to be solved as similar results are obtained when both five (in Set 2 and Set 3) or 26 (in Set 5) optimization problems are solved.
The third robustness analysis was carried out to assess the effect of the two different "knee/elbow" detection techniques (Kneedle method and L-method) on the optimal number of sensors. No differences greater than two in the optimal number of sensors were found by comparing the pairs of results for both techniques in each discrete set of numbers of sensors and function F i as the optimal number varied between 16 and 18 sensors. Furthermore, when the best function F i was considered for each set (results in bold in Table  3), the optimal number of sensors was equal to 16 or 17 (depending on the "knee/elbow" detection technique) regardless of the discrete set of numbers of sensors. Based on these results, the Kneedle method was recommended given the straightforward implementation and easier concept understanding.
Overall, the proposed methodology was robust regarding different discrete sets of numbers of sensors, fitting functions and "knee/elbow" automatic detection techniques. The optimal number of sensors for the best fitting function was equal to 16 or 17 (depending on the "knee/elbow" detection method) regardless of the chosen discrete set of numbers of sensors. This allowed for a great reduction in the computational burden associated with the overall analysis as similar results were obtained when considering both Set 2 and Set 3 (with five optimization problems to be solved) and Set 5 (26 optimization problems to be solved). Furthermore, a good approximation between 16 and 18 sensors was found considering any combination of function F (with the exception of F 1 ), any discrete set of numbers of sensors and any "knee/elbow" detection technique.
Based on these results, the optimal number of sensors was selected as 16. The optimization problem was already solved for 16 sensors (as part of Set 5), and the obtained Pareto front is presented in Figure 9a. Each grey point in Figure 9a represents a combination of 16 sensors. Figure 9b depicts a possible optimal location for the 16 sensors in green triangles, associated with the highlighted solution from the Pareto front.
considering any combination of function F (with the exception of F1), any discrete set of numbers of sensors and any "knee/elbow" detection technique.
Based on these results, the optimal number of sensors was selected as 16. The optimization problem was already solved for 16 sensors (as part of Set 5), and the obtained Pareto front is presented in Figure 9a. Each grey point in Figure 9a represents a combination of 16 sensors. Figure 9b depicts a possible optimal location for the 16 sensors in green triangles, associated with the highlighted solution from the Pareto front.

Conclusions
This article proposes a six-step methodology to determine the optimal number of pressure sensors. The major advantage of the proposed methodology is that it does not require neither prior knowledge about sensor installation costs nor the consideration of costs in the optimization problem itself to decide on the optimal number of sensors. Furthermore, the methodology can be used with different optimization methods of the pressure sensor location.
In the proposed methodology, the optimal number of pressure sensors is determined by analyzing the relationship between the results of optimization problems and the number of sensors. The total gain of installing a given number of sensors is characterized by

Conclusions
This article proposes a six-step methodology to determine the optimal number of pressure sensors. The major advantage of the proposed methodology is that it does not require neither prior knowledge about sensor installation costs nor the consideration of costs in the optimization problem itself to decide on the optimal number of sensors. Furthermore, the methodology can be used with different optimization methods of the pressure sensor location.
In the proposed methodology, the optimal number of pressure sensors is determined by analyzing the relationship between the results of optimization problems and the number of sensors. The total gain of installing a given number of sensors is characterized by using the hypervolume indicator. A trade-off function is derived to describe the gain of installing sensors as a function of the number of sensors. Finally, a technique is used to determine the optimal number of sensors as the point where the inclusion of another sensor is not worth the increment on the total gain.
A real case study was used to demonstrate, assess and discuss the performance of the proposed methodology. Three distinct robustness analyses were carried out and, according to the obtained results, the following key conclusions related to the proposed methodology can be drawn:

•
The hypervolume indicator can be used to characterize the total gain of installing pressure sensors whose locations are obtained by solving multi-objective optimization problems. • A trade-off function can be derived, allowing the characterization of the total gain of installing sensors (i.e., hypervolume) as a function of the number of sensors. This trade-off function fitting process allows for a great reduction in the computational effort associated with the overall analysis, as similar results are obtained when solving both five and 26 multi-objective optimization problems.

•
The obtained results are not affected by considering different trade-off functions or different sets of numbers of sensors as long as the selected function and the selected discrete set can properly represent the hypervolume as a function of the number of sensors. • Both Techniques 1 and 2 lead to similar results. The use of Technique 1 (Kneedle method) is further recommended given the straightforward implementation and easier concept understanding.
In future work, the robustness of the proposed methodology could be further assessed by comparing the obtained results with those obtained by using two distinct methods of locating those sensors. Furthermore, different quality measures to compare consecutive Pareto fronts and additional techniques to automatically detect the "knee/elbow" could be compared.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.