Functional Feasibility in Optimal Evaluation of Water Distribution Network Performances

: The traditional approach for the optimization of water distribution networks (WDNs) does not always lead to consistent solutions from an operational point of view. The latest optimization algorithms identify solutions that are “the best solutions” in mathematical terms but that can be less than robust against changes in operating conditions, resulting in the worst case in hydraulically infeasible conﬁgurations. Thus, this paper aims to provide a methodology that can synthesize the network performance capabilities under the change in operating conditions with two convergent strategies. The ﬁrst consists of the implementation of new performance indices (PIs), the demand deﬁcit and the pressure range, and the evaluation of their ability to criticality highlight in operating conditions. The second is the introduction of a new approach to weight the infeasible solutions in the ﬁnal result, which are those inconsistent with the real hydraulic network performances. The analysis shows that the use of these new indices makes it easier to understand the behavior of the network and to identify any weaknesses. This is true if these indices consider the hydraulically inconsistent solutions that may arise from the simulations of di ﬀ erent operation conditions; otherwise, results that poorly represent the real behavior of the network would be obtained.


Multi-Objective Optimization Algorithms
Early research on the optimization of water distribution network can be dated back to the end of the 1960s, where the most popular optimization methods were linear programming [1,2] and nonlinear programming [3]. A significant advancement in water distribution network (WDN) optimization was the introduction of heuristics methods in 1990s, as they are efficient in handling discrete variables and do not require derivative calculations. In addition, as their sampling is global, it reduces the tendency to become entrapped in local optima, and the dependency on a starting point is avoided. Among the most popular methods are genetic algorithms (GA) [4][5][6][7], simulated annealing (SA) [8][9][10], and Tabu search (TS) [11,12]. Most of these works solve basic design problems, i.e., single-objective least-cost design, using benchmark networks (two loop [2], Hanoi [13], and New York City Tunnel [1]). The evolution of cities and population growth have led to the modification and expansion of existing water networks. This resulted in a growing interest in the problem of the optimization of water networks for the least-cost design system without affecting the proper operation of the hydraulic system and consumer supply. This led to the concept of multi-objective optimization, where the problem has not one unique solution, but a Pareto optimal set of solutions, which represents the best trade-off

Account of Uncertainty in WDN Management
The approach used to deal with the problem of optimization of WDNs has evolved over the years, leading most recently to a widespread interest in future uncertainty. A limit of the deterministic approach is the lack of consideration of both uncertainty in the model parameters and in the future operating conditions [24]. Several approaches to deal with uncertainty in model parameters were analyzed by Savic [25], for example, the use of standard safety margins for uncertain parameters such as nodal demands and pipe roughness and the use of robustness evaluation methods. The results demonstrate that neglecting uncertainty in the design process may lead to severe under-design of water distribution networks. Regarding the future uncertainty, the solution of the optimization process, although optimal for the present scenario, may worsen if reality turns out to be different. To overcome this problem, a robust optimization method was developed [26], which consists of a scenario-based proactive approach in which uncertainty is included from the beginning of the decision process. It considers a set of scenarios, each one associated with a probability of occurrence, whose purpose is to make the network works well even in adverse conditions. Certainly, the price to be paid for a more robust and secure network in the future is an increase in costs.
The idea of considering a planning horizon consisting of several construction phases was also investigated [27]; this allows planners to actively manage the configuration of the network over time as new information becomes available. The flexible design strategy determines the value of the decision variables for the first phase and for several scenarios that arise from the possible different decision paths that can be followed in the future. Basupi and Kapelan [28] developed a flexible design using a multi-objective optimization with two objectives: the total cost and the resilience. This approach was tested on the Anytown and New York Tunnels networks, and the results demonstrate that the flexible design can lead to lower economic costs and a better hydraulic performance when compared with the corresponding deterministic staged design. Marques et al. used simulated annealing to solve a two-objective optimization model, including flexible design [29], and then considering the three-objective optimization problem, they minimized the costs, pressure deficits and carbon emissions [30].
Among the most important factors that can affect the future performances of a network are the change in water demands and the deterioration of the mechanical components of the network. Wagner et al. [31] suggest that both of these aspects should be considered during water distribution system design. Later, further researchers assessed the most appropriate way to quantify both forms of performance of the WDN, mainly through the evaluation of different performance indicators [32,33]. These studies show that the combination of different indices is essential to design a network that can best cope with future changes, both in terms of hydraulic behavior and resilience to potential breakage of mechanical components. Moreover, the choice of indices is crucial too, since their optimization highlights different aspects of the network performances.
The aim of the paper is to present a new procedure for the analysis of WDN performances in order to facilitate the reliable design of networks that may be more robust to possible future changes in operation conditions. Particularly, a new approach is proposed to better take into account the infeasible solutions (inconsistent with real hydraulic network behavior) that may emerge. Moreover, two new performance indices (PIs) are also proposed, the demand deficit and the pressure range, which are able to highlight critical scenarios in network operation. The suggested procedure is based on a weight system applied to PIs evaluated in different operating scenarios. Application to three benchmark networks is presented, where a two-objective optimization model is implemented (using the GALAXY MOEA) to obtain a Pareto front of solutions, from which a subset of network settings is selected, and PIs are evaluated.
The paper is structured as follows: in Section 2, the different steps of the methodology are presented; Section 3 contains the results of its application to three benchmark networks; conclusions and main findings are shown in Section 4.

Materials and Methods
The first step is the application of a multi-objective optimization algorithm to obtain a Pareto set of solutions. From the Pareto front, a subset is selected to be tested under different scenarios of demands and pipe breakage. Network performances are evaluated through a system of PIs by also weighing the infeasible solutions that may derive from the scenario simulations. The methodology adopted is set out in the flowchart in Figure 1.

Optimization Model and Solution Algorithm
The two-objective formulation considers the minimization of the cost and the maximization of the network resilience [34]. The resilience is defined as the capability of the system to adapt to internal and external changes by maintaining its operative conditions. The resilience index was first introduced by Todini [35] as a measure of network benefits. Then, an improved version of the index, the network resilience index, was proposed by Prasad and Park [36]. It considers the uniformity of nodes, because if the pipes connected to a node are not widely varied in diameter, reliable loops are ensured. The total cost of the investment for the construction of a water distribution network includes several components. In addition to the cost related to the purchase and installation of pipe components, the life-cycle cost, i.e., the total money that has to be spent on the components and processes of a WDN during its existence, has to be considered. Given the complexity of such analysis and the little significance in this type of study, only the expenditure of pipe components will be considered for the total cost. To take into account the redundancy of a pipe network, the network resilience index is used as the second objective. The corresponding formulation can be written as where C is the total cost; np is the number of pipes; a and b are constants depending on the specific problem; Dj and Lj are, respectively, the diameter and length of pipe j.
In Equation (2), In is the network resilience; nn is the number of demand nodes; Ci, qi * , hi and hi * are the uniformity coefficient, demand, actual head and minimum required head at node i, respectively; nk is the number of reservoirs; Qk and Hk are the discharge and the head relevant to each reservoir k; npu is the number of pumps in the network; Pj is the power supplied by pump j; γ is the specific weight of water. In Equation (3), npi is the number of pipes connected to node i, and Dj is the diameter of pipe j connected to demand node i.
The two objective functions (1) and (2) are subjected to a set of constraints.

Optimization Model and Solution Algorithm
The two-objective formulation considers the minimization of the cost and the maximization of the network resilience [34]. The resilience is defined as the capability of the system to adapt to internal and external changes by maintaining its operative conditions. The resilience index was first introduced by Todini [35] as a measure of network benefits. Then, an improved version of the index, the network resilience index, was proposed by Prasad and Park [36]. It considers the uniformity of nodes, because if the pipes connected to a node are not widely varied in diameter, reliable loops are ensured. The total cost of the investment for the construction of a water distribution network includes several components. In addition to the cost related to the purchase and installation of pipe components, the life-cycle cost, i.e., the total money that has to be spent on the components and processes of a WDN during its existence, has to be considered. Given the complexity of such analysis and the little significance in this type of study, only the expenditure of pipe components will be considered for the total cost. To take into account the redundancy of a pipe network, the network resilience index is used as the second objective. The corresponding formulation can be written as where C is the total cost; np is the number of pipes; a and b are constants depending on the specific problem; D j and L j are, respectively, the diameter and length of pipe j.
In Equation (2), I n is the network resilience; n n is the number of demand nodes; Ci, q i * , h i and h i * are the uniformity coefficient, demand, actual head and minimum required head at node i, respectively; n k is the number of reservoirs; Q k and H k are the discharge and the head relevant to each reservoir k; n pu is the number of pumps in the network; P j is the power supplied by pump j; γ is the specific weight of water. In Equation (3), np i is the number of pipes connected to node i, and D j is the diameter of pipe j connected to demand node i. The two objective Functions (1) and (2) are subjected to a set of constraints.
where Q i is the inflows or outflows of node i, q* i is the demand at node i, np n is the set of pipes entering or leaving node i, np loop is the pipes in a loop, ∆hj is the head loss in pipe j, h i and h imin are the actual head and the minimum requested head at node i, and D j and {D} are the diameter of pipe j and the set of commercial diameters.
Equations (4) and (5) represent the conservation laws of mass and energy, expressed in terms of algebraic summation. Equation (6) represents the set of constraints for the minimum pressure at each node. Equation (7) expresses that the value of the decision variables, i.e., diameter of each pipe must be selected from the available set of commercial diameters.
The GALAXY algorithm (genetically adaptive leaping algorithm for approximation and diversity)-version available in the database of the Centre for Water Systems at University of Exeter during the period March to July 2019-is used to solve the two-objective optimization problem. GALAXY is a hybrid MOEA developed to deal with a discrete, multi-objective design of water distribution systems. In the generation step, six search operators are used in a genetically adaptive way. Then, in the replacement process, the Pareto dominance and the -dominance concept are combined. Given an intermediate population whose individuals have been ranked in accordance to the non-dominated sorting, if the number of the top rank individuals is less than the population size, the Pareto dominance replacement is performed. If the number of the top rank individuals is higher than the population size, the -replacement is performed, where the -non-dominated solutions are transferred in the new population and any empty spaces are filled with some -dominated solutions, those nearest to the global optimum. Thus, this allows some -dominated solutions to not be completely removed, and despite the crowding distance procedure, it avoids the elimination of solutions belonging to the Pareto optimal front.
The strength of the GALAXY parametrization issue lies in the fact that some simplifications have been made to allow for as few parameters as possible to be implemented in GALAXY as few. Normally, the search operators require several parameters to be set before use. Here, the effective components of each operator are identified and tailored to suite the discrete nature of the WDN design problem. As a consequence, GALAXY only needs to specify the two general parameters required by any genetic algorithm: the population size (N) and the number of function evaluations (NFEs). These values are chosen based on the complexity of different WDNs.
Algorithm solutions are spread on a Pareto front, from which a subset of suboptimal elements is selected. This selection is aimed to reduce the number of network simulations and is performed along the front to cover different cases.

Performance Indices
The evaluation of WDN performances under different operating scenarios can be achieved by the assessment of the performance indices. The PIs quantify the reliability of a structure, seen as the capability of the network to satisfy the requested demand, with sufficient pressure. Together with the cost and the resilience index implemented in the optimization algorithm, two new PIs are proposed in this study:

•
Demand deficit (DD) evaluated as the percentage of undelivered demand. • Pressure range (PR) evaluated as the variation of the nodal pressure out of the optimal operating range.
The two indices are developed to highlight the network capability of achieving satisfactory performances. The demand deficit index is a measure of the vulnerability of the network, which represents the unacceptability of the operation conditions of the system. It is evaluated for the node with the lowest pressure value as the difference between the requested nodal discharge (Q req ) and the real delivered discharge (Q del ).
The pressure range index is a measure of the network capability to meet the design constraints for an optimal hydraulic operating condition. This range refers to the minimum pressure that allows the delivery of the requested nodal discharge and to the maximum pressure that is linked both to the maximum pipe resistance and to the working standards of the hydraulic devices connected to the network. The pressure range index is calculated for each scenario through Equation (9): where nn is the number of nodes, h i is the actual pressure at node i, and β and h ref are evaluated as follows: where h min and h max are the minimum and the maximum pressure constraints.
Although the DD highlights the local worst operating condition, i.e., the node with the lowest pressure, the PR considers all the nodes of the network, thus giving an overview of the feasibility of the operating conditions.
For the correct assessment of these PIs, a new way to consider the infeasible solutions that may derive from scenario simulations is also proposed. For each solution, the total number of infeasible configurations among all the scenarios should be considered. The infeasible configurations are those that exhibit negative pressures, thus making it impossible to provide any outflow. Therefore, the PR and DD values are multiplied for an IC index (infeasible configurations), resulting in two new indices, weighted pressure range (PR*) and weighted demand deficit (DD*): With the number of scenarios n tot and the number of infeasible configurations n inf .

Network Simulation: Pressure-Driven Analysis
For the correct evaluation of the indices, it is necessary to consider the real delivered flow rate based on the node pressure, i.e., to simulate the network through a pressure-driven analysis [37]. The optimized network configurations are obtained by running the optimization algorithm simultaneously with the hydraulic simulator EPANET [38], which works as demand-driven software (DDA). This type of analysis produces correct results since the GALAXY algorithm imposes Water 2020, 12, 3404 7 of 20 pressure restrictions to be greater than the minimum required at each node. When increasing the demand scenario, the head losses will also increase, and this may result in the pressure in some nodes being less than the allowed minimum. For this reason, a pressure-driven analysis (PDA) is needed, imposing the dependency of the flow rate delivered at the node i with the pressure (p i ). The idea is to consider the original system of equations consisting of two equations, i.e., the equation of continuity and the equation of energy conservation, and to add a third equation which relates the delivered flow to the nodal pressure. This is possible by including in the DD version an element which models the PD behavior, i.e., the emitter. The implementation in EPANET consists of carrying out the DDA and checking the nodal pressure by highlighting all the nodes in which the pressure is lower than the minimum one. For these nodes, the demand is set to zero and an emitter is located. To model an emitter to be introduced in EPANET, it is necessary to define the emitter coefficient (C i ) related to the node i, obtained from the following equation: where Q dem is the delivered flow rate at the node i, H i and z i are the total head and the elevation at node i, γ is an exponent which considers the layout of the secondary network fed by the node-for the sake of simplicity here, it is imposed equal to 0.5. Three different coefficients can be evaluated based on the total head value. Table 1 shows these different situations. Table 1. Operating range of emitter coefficients.

Benchmark Set-Up
The algorithm has been applied on three benchmark networks of increasing level of complexity; the two-loop network [2], the Hanoi network [13], and the Balerma irrigation network [39], which are classified as small, medium, and large problems, respectively [34] (Figure 2). We chose these networks because they are benchmarks widely used in the literature to compare different problems and methods in WDN studies. Table 2 gives a summary of the networks' characteristics.
Water 2020, 12, x FOR PEER REVIEW 7 of 20 PD behavior, i.e., the emitter. The implementation in EPANET consists of carrying out the DDA and checking the nodal pressure by highlighting all the nodes in which the pressure is lower than the minimum one. For these nodes, the demand is set to zero and an emitter is located. To model an emitter to be introduced in EPANET, it is necessary to define the emitter coefficient (Ci) related to the node i, obtained from the following equation: where Qdem is the delivered flow rate at the node i, Hi and zi are the total head and the elevation at node i, γ is an exponent which considers the layout of the secondary network fed by the node-for the sake of simplicity here, it is imposed equal to 0.5. Three different coefficients can be evaluated based on the total head value. Table 1 shows these different situations.

Benchmark Set-Up
The algorithm has been applied on three benchmark networks of increasing level of complexity; the two-loop network [2], the Hanoi network [13], and the Balerma irrigation network [39], which are classified as small, medium, and large problems, respectively [34] (Figure 2). We chose these networks because they are benchmarks widely used in the literature to compare different problems and methods in WDN studies. Table 2 gives a summary of the networks' characteristics.    The GALAXY algorithm needs two parameters to be set, N and NFEs, whose values are taken from a previous study [34]. To cover a wide range of the best-known Pareto fronts, three different population sizes are implemented, different for each of the three benchmark problems. For each group, 10 runs are performed, maintaining a fix number of function evaluations, which increases with the network complexity. For each network, a total of 30 runs are performed (Table 3). The input parameters for the optimizer are (i) the set of commercial diameters, different for the three networks; (ii) diameter costs, expressed in USD per meter; and (iii) minimum pressure constraints. Table 4 shows the detailed values.

Pareto Fronts and Solution Subset Selection
For the analysis of the scenario simulations, Pareto fronts are needed, thus, the GALAXY algorithm is implemented to identify the Pareto optimal solutions for the three problems.
The whole analysis is developed by using the EPANET-Matlab toolkit [40]. For each network, solutions obtained by GALAXY algorithm are compared with the BKPF (best-known Pareto front) [34] (numerical results can be found at Design/Resiliance | Engineering | University of Exeter), and only the solutions which are better or equal in one or both of the objectives with respect to the BKPF are retained. This operation results in three new GALAXY fronts that are made up of 111 solutions for the two-loop network (TLN) (Figure 3), 618 solutions for the Hanoi network (HAN) (Figure 4) and 529 solutions for the Balerma irrigation network (BIN) ( Figure 5).
An analysis of the results shows that GALAXY can assess all of the solutions for the two loop and even more solutions for the Hanoi when compared to the BKPF. For the Balerma irrigation, it is unlikely that the algorithm will fail to produce better solutions than those already known in the area of lower resilience value. In fact, it produces better solutions when compared to the BKPF only for a cost higher than 5.72 × 10 6 USD. Hence, we retain the GALAXY solutions only in the upper part of the front. In contrast, in the part of the front below 5.72 × 10 6 USD, since we need network configurations for the subsequent evaluations, we chose to consider the solutions produced by the latest version of the simulated annealing algorithm (MOSA-GR solutions), which gave rise to better fronts with respect to the BKPF [23] (Figure 4). Following the completion of the fronts calculation, the solutions to be used in the following evaluations can be selected. Well-spread solutions along the fronts are chosen with the aim of providing a range of cases as general as possible. The chosen solutions are represented in black in the following graphs, and their objective values are reported in Table 5. In total, seven solutions for the TLN, six solutions for the HAN and seven solutions for the BIN are selected.
Water 2020, 12, x FOR PEER REVIEW 9 of 20

Operation Scenarios and Performance Assessment
To assess the reliability of the system, two different aspects have to be considered: the hydraulic reliability and the mechanical reliability. The former reflects how well the system can cope with changes in water demands, and the latter express the suitability of the system to continue to work under breakage in mechanical components, e.g., pipe failure. Therefore, the networks are tested under different demand scenarios and pipe failures.
For each chosen solution, three different demand scenarios are proposed:  A 10% increase in the demand in all nodes;  A 30% increase in the demand in 1/3 nodes of the network with the maximum demand;  A 30% increase in the demand in 1/3 nodes of the network with the minimum demand.
For the BIN, which has an equal demand set to all nodes, the 30% increase in demand is carried out in one-third of the nodes, with the maximum and minimum pressure recorded in the optimal configurations.
To examine the effect of individual pipe failure against the available level of supply, several pipes belonging to different areas of the networks (generally close to the reservoirs), are selected in a heuristic way. In Figure 6 and Table 6, the selected pipes for the three networks and their characteristics are reported. Those pipes are closed one at a time and the network is simulated. First, scenarios of increase demand and pipe closures are considered individually, then combined. Table 7 shows the details of scenario implementation. Scenarios 4, 5, 6, 7 are divided in turn into other scenarios for the three networks considering one pipe closure at a time.
It has to be pointed out that the simulations are performed in a single time step for the different scenarios. The single-phase simulation is used to simplify the analysis, as the aim is to assess through PIs the adaptability of existing networks to satisfy the water demand under predefined scenarios. Thus, some illustrative scenarios are considered, and if the network is able to front the increasing demands and pipe failures in all the situations, it will be considered adequate; otherwise, the network will be classified as inadequate for some aspects.

Operation Scenarios and Performance Assessment
To assess the reliability of the system, two different aspects have to be considered: the hydraulic reliability and the mechanical reliability. The former reflects how well the system can cope with changes in water demands, and the latter express the suitability of the system to continue to work under breakage in mechanical components, e.g., pipe failure. Therefore, the networks are tested under different demand scenarios and pipe failures.
For each chosen solution, three different demand scenarios are proposed: • A 10% increase in the demand in all nodes; • A 30% increase in the demand in 1/3 nodes of the network with the maximum demand; • A 30% increase in the demand in 1/3 nodes of the network with the minimum demand.
For the BIN, which has an equal demand set to all nodes, the 30% increase in demand is carried out in one-third of the nodes, with the maximum and minimum pressure recorded in the optimal configurations.
To examine the effect of individual pipe failure against the available level of supply, several pipes belonging to different areas of the networks (generally close to the reservoirs), are selected in a heuristic way. In Figure 6 and Table 6, the selected pipes for the three networks and their characteristics are reported. Those pipes are closed one at a time and the network is simulated. First, scenarios of increase demand and pipe closures are considered individually, then combined. Table 7 shows the details of scenario implementation. Scenarios 4, 5, 6, 7 are divided in turn into other scenarios for the three networks considering one pipe closure at a time.
It has to be pointed out that the simulations are performed in a single time step for the different scenarios. The single-phase simulation is used to simplify the analysis, as the aim is to assess through PIs the adaptability of existing networks to satisfy the water demand under predefined scenarios. Thus, some illustrative scenarios are considered, and if the network is able to front the increasing demands and pipe failures in all the situations, it will be considered adequate; otherwise, the network will be classified as inadequate for some aspects.    30% increase in the demands in 1/3 nodes with the maximum demand + pipe failure; 7 30% increase in the demands in 1/3 nodes with the minimum demand + pipe failure.

Results and Discussion
In this section, the results of the analysis of networks behavior under different scenarios of water demands and pipe failures are presented. To carry out the performance analysis, first, the Pareto fronts of the three benchmark networks are evaluated, then subsets of the optimal design solutions are selected to be used as a new domain for scenario simulation.
For all simulations of each network, the PIs are evaluated and compared with the optimal configuration values ( Table 8). In particular, the network resilience index (I n ) is already inside the optimization algorithm, whereas the pressure range (PR), the minimum nodal pressure (Pmin) and the corresponding demand deficit (DD) are evaluated ex post for each solution (Sol). For the simulated configurations, the percentage of infeasible solutions (IC) among all the scenarios is also considered, and it is used to weight the pressure range and the demand deficit indices, resulting in weighted pressure range (PR*) and weighted demand deficit (DD*), respectively. In order to provide a global view on the ability of the three networks to cope with variability in their operation, in the following paragraphs, we discuss for each network the results of the simulations, providing the average values of the PIs considering the scenarios altogether. Table 9 provides the average results of the scenario simulations. The configuration associated with the lower cost shows an overall average reduction of 12.9% of the value of resilience with respect to the value of the optimize configuration (Table 8), 0.216 compared to 0.248, respectively. Despite the criticalities characterizing this configuration and the fact that the resilience index is the lowest among all the solutions, its reduction in regard to the initial value is lower compared to the two subsequent solutions, i.e., 28.4% for solution 2 and 34.87% for the solution 3. By increasing the cost, the resilience index also increases, which means higher reliability in general. For example, the average I n of configuration 2 is 0.251; even if this value is higher than the first solution, this configuration shows greater criticalities related to the failure of pipe 4 in terms of minimum average pressure and average demand deficit of 24.73 m and 11.49%, respectively. This is because the cost increase in solution 2 in comparison to the solution 1 is caused by the diameter of pipe 4 passing from 25.4 mm to 254 mm ( Figure 7). Thus, when this pipe is closed, the higher diameter of pipe 6 in solution 1 (152.4 mm compared to 101.6 mm of solution 2) causes high pressure at node 6 for the first solution. Even for the third solution, which is characterized by a resilience of 0.550, the situation remains worse than the first configuration with an average Pmin of 27.96 m, and an average demand deficit of 7.34%. It has to be pointed out that since this configuration does not show infeasible solutions, the DD* and the PR* are both zero. Thus, by weighting the PIs with the IC, the third solution looks more reliable than the solution 1. It should be stressed that the higher demand deficit of solutions 2 and 3 with respect to the solution 1 is due to the fact that, despite solution 1, the two solutions do not exhibit infeasibility deriving from the failure of pipe 2. Therefore, the low pressures recorded in both configurations, with corresponding high demand deficit, result in lower average pressures and higher demand deficit compared to the first configuration. The difference is that solution 2 still presents an infeasibility of 27.27% among the simulated scenarios, whereas this is not the case in the third solution, resulting in it possessing overall more reliability than the first ones. Solution 4 is the first solution that does not exhibit demand deficit even in the worst case of pipe failure, in addition to the increase in water demand. This is due to the fact that, despite the fact that solution 3 shows most of the pipes' diameters of the network to be equal to the previous solution (solution 2), solution 4 is characterized by an increase between 10% and 28% in all of the diameters with respect to solution 3, which ensures a better arrangement of the network to face variability in its working conditions. Therefore, from the fourth solution on, the analysis does not reveal criticalities, providing that such configurations can be considered reliable. for the third solution, which is characterized by a resilience of 0.550, the situation remains worse than the first configuration with an average Pmin of 27.96 m, and an average demand deficit of 7.34%. It has to be pointed out that since this configuration does not show infeasible solutions, the DD* and the PR* are both zero. Thus, by weighting the PIs with the IC, the third solution looks more reliable than the solution 1. It should be stressed that the higher demand deficit of solutions 2 and 3 with respect to the solution 1 is due to the fact that, despite solution 1, the two solutions do not exhibit infeasibility deriving from the failure of pipe 2. Therefore, the low pressures recorded in both configurations, with corresponding high demand deficit, result in lower average pressures and higher demand deficit compared to the first configuration. The difference is that solution 2 still presents an infeasibility of 27.27% among the simulated scenarios, whereas this is not the case in the third solution, resulting in it possessing overall more reliability than the first ones. Solution 4 is the first solution that does not exhibit demand deficit even in the worst case of pipe failure, in addition to the increase in water demand. This is due to the fact that, despite the fact that solution 3 shows most of the pipes' diameters of the network to be equal to the previous solution (solution 2), solution 4 is characterized by an increase between 10% and 28% in all of the diameters with respect to solution 3, which ensures a better arrangement of the network to face variability in its working conditions. Therefore, from the fourth solution on, the analysis does not reveal criticalities, providing that such configurations can be considered reliable.

Hanoi Network
The low value of resilience among all of the scenarios is initially notable when observing the results (Table 10). It reaches at best a value of 0.326, far below the unit. This means that the network, even in the most expensive configuration, will still be characterized by low resilience; moreover, to increase the resilience, a substantial cost increase is needed. The reason for this result may be the

Hanoi Network
The low value of resilience among all of the scenarios is initially notable when observing the results (Table 10). It reaches at best a value of 0.326, far below the unit. This means that the network, even in the most expensive configuration, will still be characterized by low resilience; moreover, to increase the resilience, a substantial cost increase is needed. The reason for this result may be the configuration of the network, which has a mesh with consecutive nodes on the same branch, without sufficient ramifications. This means that any change in diameter of a single pipe affects the whole mesh without any possibility for the flow to find alternative routes, resulting in a low resilience. Starting from the solution 1, e.g., with a cost of 6.258 × 10 6 USD and a network resilience of 0.219 (Table 8), the simulation of the scenarios shows an infeasibility of 60% and six solutions with demand deficit. The average DD* among the solutions is 18.54%, and the infeasibility derives mainly from the combined scenarios of increase in demands and pipe failures. In the other tail of the front, the most expensive solution still shows a hydraulic inconsistency only in the scenario of pipe 20 failure combined with the 30% increase in water demand in one-third of the nodes with the highest demand. Despite the previous configurations, in the other demand scenarios combined with the failure of pipe 20, the operating conditions are ensured; however, with an average demand deficit of 20.74%, which is almost the same as that of solution 1, but weighted for the IC, it decreases to 1.99%. In these cases, the pressures below the low limit are almost recorded in nodes 20-33; this is because the failure of pipe 20 leads to inhomogeneity in water distribution mainly in the right area of the network (Figure 6b has an extra deficit configuration-2 against 3, respectively. Nevertheless, a significant difference between the two solutions is the average percentage of demand deficit, 4.21% for solution 4 and 1.48% for solution 5, due to a combination of higher diameters close to the reservoir which result in a lower demand deficit. It is interesting to note that the failure of pipe 20 alone (scenario 4) involves infeasible configurations, with negative pressures mainly in the nodes in the left part of the network (nodes 20-32), for all the solutions except the sixth one. Here, the importance of weighting the indices with respect to the infeasible solutions is evident. In fact, without considering the weight, the average PR and DD of the simulations would be 0.332 and 20.74%, respectively, basically the worst among all of the configurations. In contrast, by multiplying these PIs for the IC, their values decrease.
It may be useful to provide a more detailed analysis of the behavior of a single configuration. Therefore, we chose to analyze the worst configuration, which, in this case, is solution 3 subjected to the simultaneous increment of 30% in water demand in one-third of the most demanding nodes and the failure of pipe 13 (Figure 8). In this condition, an absolute minimum pressure of 16.45 m at node 13 is observed, resulting in a demand deficit of 35.97%. Despite the significant local failure, it is important to evaluate the pressure trend all over the network to identify how much the network is affected by low pressures. For this purpose, we use a histogram of absolute frequencies of the nodal pressures of the network, representing the number of nodes with pressure values falling within a given class (Figure 9). By observing the graph, we can argue that, globally, the network presents pressures with values inside the optimal operating range (30-80 m)-27 of 31 nodes with respect to the pressure limitations. Only one node exhibits a pressure higher than 90 m (node 2, near the reservoir), and four nodes are characterized by demand deficit. To conclude, the worst pressure value of 16.45 m is isolated and concerns only one node, thus, excluding a few local failures, the configuration can be considered quite reliable.

Balerma Irrigation Network
From the results of the simulations, the one-by-one breakage of the four pipes linked to the sources (pipes 188, 194, 51, 338, Figure 6c) causes negative pressures in some nodes for all the seven solutions. This is true both considering the pipe failures alone and combining them with the increase in water demands. We report in Table 11, only the outcomes for the most expensive solution as a representative, because all of the solutions show 12 infeasibilities with respect to the same scenarios. Therefore, since the IC is the same for all of the solutions, its weight on the PIs is less relevant with respect to the previous networks. On the other hand, the breakage of pipes far from the sources (pipes 480, 119 and 75; Figure 6c) does not cause infeasible solutions but only solutions with pressures below the minimum required. From Table 12, we can see that in addition to inconsistent solutions, solution

Balerma Irrigation Network
From the results of the simulations, the one-by-one breakage of the four pipes linked to the sources (pipes 188, 194, 51, 338, Figure 6c) causes negative pressures in some nodes for all the seven solutions. This is true both considering the pipe failures alone and combining them with the increase in water demands. We report in Table 11, only the outcomes for the most expensive solution as a representative, because all of the solutions show 12 infeasibilities with respect to the same scenarios. Therefore, since the IC is the same for all of the solutions, its weight on the PIs is less relevant with respect to the previous networks. On the other hand, the breakage of pipes far from the sources (pipes 480, 119 and 75; Figure 6c) does not cause infeasible solutions but only solutions with pressures below the minimum required. From Table 12, we can see that in addition to inconsistent solutions, solution

Balerma Irrigation Network
From the results of the simulations, the one-by-one breakage of the four pipes linked to the sources (pipes 188, 194, 51, 338, Figure 6c) causes negative pressures in some nodes for all the seven solutions. This is true both considering the pipe failures alone and combining them with the increase in water demands. We report in Table 11, only the outcomes for the most expensive solution as a representative, because all of the solutions show 12 infeasibilities with respect to the same scenarios. Therefore, since the IC is the same for all of the solutions, its weight on the PIs is less relevant with respect to the previous networks. On the other hand, the breakage of pipes far from the sources (pipes 480, 119 and 75; Figure 6c) does not cause infeasible solutions but only solutions with pressures below the minimum required. From Table 12, we can see that in addition to inconsistent solutions, solution 1 and solution 2 feature the mostly simulations in which some nodes show pressures below 20 m with an average demand deficit of 39.4% and 9.52%, respectively. From the third solution onwards, the pressure recorded in each node is higher than the minimum; therefore, there are no solutions with demand deficit, but the pressure at some nodes is even greater than the maximum pressure of 80 m usually imposed by the regulation (Table 13). This is reflected in a PR* that increases from solution 3 onwards. This is a criticism of the PR index since it considers likewise the upper and the lower pressure exceeded out from the optimal range, even if the cases of pressure below the lower limit are more critical than the others; the first ones result in demand deficit, increasing the vulnerability of the system, whereas in the second, the inclusion of valves can be dealt with. Solution 3, instead, due to the increase in the diameter of 209 pipes (generally near the reservoirs), with respect to solution 2, is the first configuration that is able to ensure the lack of demand deficit under every scenario, while ensuring lower operating pressure (104.5 m) when compared to the subsequent solutions and fairly lower costs.  By proposing the same analysis carried out for the HAN network, we consider the most critical solution to analyze the trend of the nodal pressures in the network. The most demanding condition is configuration 1 with the increase of 30% in demand in one-third of the highest pressure nodes in addition to the failure of pipe 480. It exhibits the highest demand deficit of 88.65 % due to a minimum pressure of 9.52 m. From the histogram of the absolute frequencies of nodal pressures (Figure 10), it can be observed that only one node shows a pressure below 10 m and that 29 nodes have pressures between 10-20 m and thus are still below the minimum required. Therefore, a total of 30 nodes are characterized by demand deficit. However, considering that the BIN network has 443 nodes, we can state that most of the nodes work at acceptable pressure values of between 20 and 70 m. Although the deficit cases are outnumbered, they may be localized and possibly considered if they affect sensitive targets, for example. Once again, despite the heavy deficit of 88.65%, it represents an isolated case, and the network performs well overall in facing the changes in working conditions.
Water 2020, 12, x FOR PEER REVIEW 17 of 20 is configuration 1 with the increase of 30% in demand in one-third of the highest pressure nodes in addition to the failure of pipe 480. It exhibits the highest demand deficit of 88.65 % due to a minimum pressure of 9.52 m. From the histogram of the absolute frequencies of nodal pressures (Figure 10), it can be observed that only one node shows a pressure below 10 m and that 29 nodes have pressures between 10-20 m and thus are still below the minimum required. Therefore, a total of 30 nodes are characterized by demand deficit. However, considering that the BIN network has 443 nodes, we can state that most of the nodes work at acceptable pressure values of between 20 and 70 m. Although the deficit cases are outnumbered, they may be localized and possibly considered if they affect sensitive targets, for example. Once again, despite the heavy deficit of 88.65%, it represents an isolated case, and the network performs well overall in facing the changes in working conditions.

Conclusions
The optimal design of water distribution networks still represents a challenge issue. This is due, on the one hand, to the continuous growth of urban centers and the corresponding evolution of WDNs and, on the other hand, to the high rate of uncertainty that characterized the WDNs' operation during their life cycle. Among the main factors that may be considered, only the breakage of pipes and the change in water demands are analyzed. Other factors, such as the operating conditions of pumps and valves and the networks extension, are considered of minor importance, and therefore are not considered.
In this complex framework, the challenge is to find a design solution which is the best in terms of performance in the present and, at the same time, also in the future. The traditional approach based on the evaluation of the Pareto optimal solutions and the successive use of classic PIs to rank them seems not always satisfactory. Particularly, the analysis of resilience and cost indices is not sufficient when choosing among optimal or sub-optimal solutions. Moreover, there is the risk that "optimal" solutions reveal themselves to be unreliable in terms of hydraulic modelling or unacceptable in terms of service compliance, for example, due to violation of pressure constrains or non-fulfilment of service goals. Thus, this study aims to outline a methodological approach to better understand the behavior of the WDNs with respect to changes in working conditions. This approach focuses on two important aspects, the evaluation of the demand fulfilment and the range of working pressures, through two new indices, the demand deficit (DD), and the pressure range (PR). The two PIs represent a measure of the vulnerability of the network, and the more they tend to zero, the more robust the network is from an operational point of view. Moreover, in the proposed approach, the importance of also taking into account the hydraulic and service feasibility of solutions in all the testing scenarios is stressed, and the use of a weighting factor, i.e., the percentage of infeasible solutions (IC), is suggested.

Conclusions
The optimal design of water distribution networks still represents a challenge issue. This is due, on the one hand, to the continuous growth of urban centers and the corresponding evolution of WDNs and, on the other hand, to the high rate of uncertainty that characterized the WDNs' operation during their life cycle. Among the main factors that may be considered, only the breakage of pipes and the change in water demands are analyzed. Other factors, such as the operating conditions of pumps and valves and the networks extension, are considered of minor importance, and therefore are not considered.
In this complex framework, the challenge is to find a design solution which is the best in terms of performance in the present and, at the same time, also in the future. The traditional approach based on the evaluation of the Pareto optimal solutions and the successive use of classic PIs to rank them seems not always satisfactory. Particularly, the analysis of resilience and cost indices is not sufficient when choosing among optimal or sub-optimal solutions. Moreover, there is the risk that "optimal" solutions reveal themselves to be unreliable in terms of hydraulic modelling or unacceptable in terms of service compliance, for example, due to violation of pressure constrains or non-fulfilment of service goals. Thus, this study aims to outline a methodological approach to better understand the behavior of the WDNs with respect to changes in working conditions. This approach focuses on two important aspects, the evaluation of the demand fulfilment and the range of working pressures, through two new indices, the demand deficit (DD), and the pressure range (PR). The two PIs represent a measure of the vulnerability of the network, and the more they tend to zero, the more robust the network is from an operational point of view. Moreover, in the proposed approach, the importance of also taking into account the hydraulic and service feasibility of solutions in all the testing scenarios is stressed, and the use of a weighting factor, i.e., the percentage of infeasible solutions (IC), is suggested.
The methodology is tested on three case studies, the two loop, Hanoi, and Balerma irrigation networks, which are standard references in the international literature.
The results show that the two proposed indices are able to highlight the low reliability in terms of minimum pressure and demand fulfilment of some solutions, in addition to the evaluation of the resilience index alone. However, even for the most "critical scenarios", solutions characterized by higher diameters, and thus by higher values of the resilience index, are expected to be associated with more satisfactory ranges of nodal pressures, but this is not always true due to their higher degree of spatial variation. This is particularly clear for the more complex network, the Balerma irrigation.
Moreover, a significant percentage of infeasible solutions is also observed for the scenarios that drive to more expensive and thus more resilient networks. It is interesting to note that for the networks of Hanoi and Balerma, the values of index IC are high and quite constant for all the scenarios, highlighting the limits of pure mathematical optimization procedures. In fact, the results show that carrying out a performance analysis only on average values of classical PIs, without also considering the number of infeasible configurations, may not lead to fully reliable outcomes. The joint use of the new PIs, DD, PR, and IC may be helpful to improve the analysis of functionality of WDNs, widening the range of factors to be considered in optimization procedures. This advantage may be confirmed in application to other case studies that will be performed in future research.