*2.1. Procedure 1—WNP*

According to Perelmann et al. (2015) [41], WNP is carried out in two main phases: (a) *clustering*, in which the optimal shape and size of the clusters are defined by minimizing the number of edge cuts (boundary pipes) and by simultaneously balancing the number of nodes of each cluster, and (b) *dividing*, in which clusters are separated from each other by closing isolation valves at some boundary pipes and installing flow meters at the remaining boundary pipes.

In this work, the clustering layout is obtained exploiting the properties of the normalized Laplacian matrix **L** = **D** − **A**, in which **D** is the diagonal matrix containing the node degree *ki* of each node, and **A** is the adjacency matrix. In this matrix, the elements *aij* = *aji* = 1 if nodes *ni* and *nj* are connected by a pipe; otherwise, *aij* = *aji* = 0. Shi and Malik (2000) [42] demonstrated that through the first *C* smallest eigenvector of the normalized Laplacian matrix, the relaxed version of the min-cut problem can be solved. In fact, it corresponds to the minimization of the Rayleigh quotient. If *C* is the number of clusters in which the network must be divided, the first *C* smallest eigenvectors of the Laplacian matrix are considered and used to create a new matrix **UnxC**. A *k*-means algorithm is applied to the rows of **UnxC** for grouping the nodes of the network in *C* clusters. The main trick is to change the representation of the nodes in the eigenspace of the first *C* eigenvectors, which enhances the cluster-properties of the nodes in such a way that they can be trivially detected in the new representation. The spectral clustering algorithm proved to show a superior performance to other clustering procedures, in that the provided clustering layout features both a well-balanced cluster size and a minimum number of edge cuts [43]. The main spectral clustering steps in the case of a WDN are described by Di Nardo et al. (2018b) [44]. The graph of the WDN can be considered un-weighted (every connection between the nodes has the same importance, *aij* = *aji* = 1) or weighted (the value *aij* = *aji* can be related to pipe features, such as diameter *d* and length *l*). In the applications of this work, *aij* and *aji* were set at 1. The optimal number of clusters *C* (from a topological point of view) in which to subdivide the network is chosen as a function of the number *n* of nodes, according to the relationship *Copt* = *n*0.28 [45]. The clustering phase provides the optimal cluster layout and, as a result, the edge-cut set, consisting of a group of *Nec* boundary pipes between clusters. In correspondence to each boundary pipe, the flow transfer between the adjacent clusters must always be known, if it is larger than zero, in order to make the dividing effective. Therefore, the choice must be made whether either a gate valve must be closed, or a flow meter must be installed in the generic boundary pipe. Following this choice, the sum of closed gate valves (as numerous as *Ngv*) and installed flow meters (as numerous as *Nfm*) must be equal to *Nec*. Closing gate valves has the effect of reducing the service pressure and, therefore, leakage in the WDN. However, if service pressure falls below the desired threshold value *hdes*, this negatively impacts on WDN reliability. In this work, the trade-off between leakage and WDN reliability was explored through the bi-objective optimization, performed through the NSGAII genetic algorithm [46]. In this optimization, several decisional variables equal to *Nec* was considered, to encode, inside individual genes, gate valve closure (gene value equal to 1) or flow meter installation (gene value equal to 0) at boundary pipes. The first objective function *f1* to minimize was the daily leakage volume *Vl* (m3):

$$f\_1 = V\_1 \tag{1}$$

where *Vl* is calculated as the sum of the temporal integral of the nodal leakage outflows, evaluated as a function of nodal pressure heads through the Tucciarelli et al. (1999) [47] formula.

The second objective function *f2* relates to the global resilience failure index *GRF* index proposed by Creaco et al. (2016) [48] to represent the instantaneous power surplus/deficit conditions of the WDN. In fact, *GRF* is dimensionless and is the sum of the resilience (*Ir*) and failure (*If*) indices evaluated at the generic instant of WDN operation:

$$GRF = I\_r + I\_f = \frac{\max\left(q\_{user}^T H - d^T H\_{des}, 0\right)}{Q\_0^T H\_0 - d^T H\_{des}} + \frac{\min\left(q\_{user}^T H - d^T H\_{des}, 0\right)}{d^T H\_{des}}\tag{2}$$

where *d* and *quser* are the vectors of nodal demands *d* (m3/s) and water discharges *quser* (m3/s) delivered to users, respectively, at WDN demanding nodes. In this work, *quser* was evaluated as a function of *d* and pressure head *h* (m) at each node though the pressure-driven formula of Tanyimboh and Templeman (2010) [49], with calibration proposed by Ciaponi et al. (2014) [50]. *H* and *H0* are the vectors of nodal heads (m) at demanding nodes and sources, respectively. *Hdes* is the vector of desired nodal heads, which are the sum of nodal elevations and desired pressure heads *hdes* (m). Finally, *Q0* is the vector of the water discharges leaving the sources. The *GRF* index has the advantage of being within range [−1, 1]. Higher values of *GRF* indicate higher power delivered to WDN users and, therefore, higher service pressure. With reference to WDN daily operation, the second objective function *f2* to maximize was calculated with the following relationship, as suggested by Creaco et al. (2016) [48]:

$$f\_2 = \text{median(CRF)}\tag{3}$$

The choice of the median value of *GRF* is because Creaco et al. (2016) [48] proved it to give a suitable and concise representation of a sequence of operation scenarios in the extended period simulation of the WDN. Both *f1* and *f2* can be calculated by applying a pressure-driven WDN solver (e.g., that of Creaco et al. 2016 [48]). They are mutually contrasting objectives: in fact, as the number of closed gates grows, *f1*, which has to be minimized, decreases. At the same time, *f2*, which has to be maximized, decreases as well due to the decreasing service pressure. This creates a trade-off between the two objectives, which takes the form of a Pareto front of optimal solutions, that is a group of solutions from which to select the final solution for the partitioning. To this end, additional criteria, such as the partitioning cost or demand satisfaction, can be adopted. In fact, the Pareto front of optimal solutions can be re-evaluated in terms of other functions, such as number *Nfm* of flow meters and demand satisfaction rate *Ids*. In fact, *Nfm* is a surrogate for the partitioning cost [34], whereas *Ids* represents the effectiveness of the service to WDN users. The latter index can be calculated as:

$$I\_{ds} = \frac{w\_d}{w\_{tot}}\tag{4}$$

where *wd* (m3) and *wtot* (m3) are the delivered water volume and the WDN demand, respectively. Variable *wd* can be calculated starting from the temporal integral of the water discharge *quser* delivered to the users at each node.
