Bi-objective design-for-control of water distribution networks with global bounds

This manuscript investigates the design-for-control (DfC) problem of minimizing pressure induced leakage and maximizing resilience in existing water distribution networks. The problem consists in simultaneously selecting locations for the installation of new valves and/or pipes, and optimizing valve control settings. This results in a challenging optimization problem belonging to the class of non-convex bi-objective mixed-integer non-linear programs (BOMINLP). In this manuscript, we propose and investigate a method to approximate the non-dominated set of the DfC problem with guarantees of global non-dominance. The BOMINLP is first scalarized using the method of ϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document}-constraints. Feasible solutions with global optimality bounds are then computed for the resulting sequence of single-objective mixed-integer non-linear programs, using a tailored spatial branch-and-bound (sBB) method. In particular, we propose an equivalent reformulation of the non-linear resilience objective function to enable the computation of global optimality bounds. We show that our approach returns a set of potentially non-dominated solutions along with guarantees of their non-dominance in the form of a superset of the true non-dominated set of the BOMINLP. Finally, we evaluate the method on two case study networks and show that the tailored sBB method outperforms state-of-the-art global optimization solvers.


Introduction
The optimal design of water distribution networks (WDN) is a classical problem which has traditionally consisted in selecting pipe diameters, for a given network topology, with the objective of minimizing design cost. Recently, renewed interest for the problem of optimal WDN design has emerged, as a result of the need to expand, rehabilitate or sectorize networks, to cope with the deterioration of assets and increasing customer demands (Giudicianni et al. 2020). In this context, the problem formulation no longer assumes a fixed network topology, but aims to simultaneously optimize network design and control, taking into account a wider range of objectives. In the UK, the requirements of the water services regulation authority to minimize water losses and interruptions of customer supply (Ofwat 2019) have led to the implementation of pressure control and network rehabilitation measures. This work investigates the optimal design-for-control (DfC) of existing WDNs with the objectives of minimizing pressure induced background leakage, which makes up a large part of total water losses, and maximizing network resilience, which reflects the ability of a WDN to maintain continuous customer supply. To the best of the authors' knowledge, this problem has not been considered in previous literature.
The presented design-for-control problem consists in simultaneously selecting elements to be added to existing WDNs from predefined sets of candidate valves (CNV) and pipes (CNP), and optimizing the controls of new and existing pressure control valves. WDNs are modeled as directed graphs, subject to the conservation of mass and energy. Flows through network links as well as pressure heads at network nodes, and head differences across pipes and valves are represented by continuous variables, while the selection status of CNVs and CNPs is represented by binary variables. Friction losses across network links are modeled using the empirical Hazen-Williams or Darcy-Weisbach equations. Finally, we consider the objectives of minimizing pressure induced background leakage and maximizing network resilience, represented respectively by the average zone pressure (AZP) (Wright et al. 2015) and the resilience index (I r ) (Todini 2000). This results in a non-convex bi-objective mixed-integer non-linear program (BOMINLP).
As AZP and I r are both increasing functions of pressure, minimizing AZP and maximizing I r are conflicting objectives. As a result, solving the DfC problem requires identifying the best compromises between the objectives, which are characterized using the notions of efficiency and non-dominance: a solution in the decision space is said to be efficient or Pareto optimal if its image in the objective space is nondominated, i.e. no other feasible solution can improve on one of the objectives without making the others worse (Miettinen 1998). In the case of bi-objective mixed-integer linear programs, the non-dominated set (or Pareto front) consists of a subset of the convex hull of a finite set of non-dominated points (Yu and Zeleny 1975). As a result, it is possible to obtain a complete description of the Pareto front using multi-objective branch-and-bound (Mavrotas and Diakoulaki 1998). More recently, De Santis et al.
(2020) proposed a branch-and-cut method to find the finite set of non-dominated points of purely integer bi-objective programs. However, Pareto fronts of multi-objective mixed-integer non-linear programs (MOMINLP) may be composed of disconnected segments of continuous curves and isolated points (Das 2000). Hence, obtaining a com-plete description of the efficient and non-dominated sets (resp. Pareto set and Pareto front) can be impractical in the case of MOMINLPs. The heuristic algorithm proposed by Cacchiani and D'Ambrosio (2017), which returns a subset of non-dominated points, is the first general purpose method for the solution of convex MOMINLPs. In this paper, we are interested in methods to approximate the Pareto front of the considered non-convex bi-objective DfC problem.
Heuristic approaches, based on genetic algorithms (GAs) or simulated annealing, and mathematical optimization methods have been applied to the design and/or control of WDNs, including optimal pipe sizing (Xu and Goulter 1999;Prasad and Tanyimboh 2009;Cunha and Marques 2020), network rehabilitation (Walters et al. 1999;Alvisi and Franchini 2009) and valve placement and control (Nicolini and Zovatto 2009;Eck and Mevissen 2012;Dai and Li 2014;Pecci et al. 2017a, c) problems. GAs and simulated annealing require a significant number of hydraulic simulations and objective function evaluations. Moreover, they do not provide global optimality (in the singleobjective case) or efficiency/non-dominance (in the multi-objective case) bounds to guarantee the quality of the computed solutions. Global optimization methods, on the other hand, return, for single-objective problems, a feasible solution with certified global optimality bounds (Eiger et al. 1994;Gleixner et al. 2012;D'Ambrosio et al. 2015;Pecci et al. 2018). For specific cases (e.g. convex mixed-integer problems or non-convex continuous problems), global optimization methods have been extended to generate guarantees of global efficiency/non-dominance for multi-objective problems, using scalarization approaches (Fernández and Tóth 2007;Ehrgott and Gandibleux 2007) or multi-objective branch-and-bound (De Santis et al. 2019;Niebling and Eichfelder 2019). Instead of a strict subset of the Pareto front, these methods return a set of feasible solutions and a guarantee of their non-dominance, within a given tolerance, in the form of a tight superset of the efficient and/or non-dominated sets.
In this paper, we present a method to approximate the non-dominated set of the design-for-control of existing water distribution networks problem using a superset of the Pareto front. The main contributions of this work include the formulation of the bi-objective design-for-control problem for the joint minimization of AZP and maximization of I r in WDNs, the development of a method to approximate the Pareto front of the resulting BOMINLP and compute potentially non-dominated solutions with global non-dominance bounds, and, finally, an equivalent linear reformulation to facilitate the computation of global optimality bounds for the maximization of the resilience index I r .
In Sect. 2, we formulate the DfC problem as a non-convex BOMINLP, and Sect. 3 presents a novel tailored method to approximate its Pareto front. First, the BOMINLP is scalarized into a sequence of single-objective mixed-integer non-linear problems (MINLP) using an -constraint method; a tailored spatial branch-and-bound (sBB) method is then applied to compute feasible solutions of the single-objective MINLPs with certified global optimality bounds. In particular, the formulation of the resilience index includes linear-fractional terms, resulting in a non-linear function. As a consequence, previous literature has solved optimization problems involving the resilience index using meta-heuristic approaches (Todini 2000;Atkinson et al. 2014;Alvisi 2015;Creaco et al. 2016). We adapt a previous sBB algorithm from the literature ) and propose an equivalent linear reformulation of I r , which enables the lineariza-tion of the parametrized MINLPs and the computation of lower bounds throughout the sBB search. We show that the proposed approach returns a set of potentially non-dominated solutions as well as global non-dominance bounds, in the form of a superset of the Pareto front. In Sect. 4, we apply this method to two case study networks. We show that the proposed sBB method outperforms state-of-the-art global optimization solvers SCIP and BARON, and converges to good quality feasible solutions. We comment on the computational results and hydraulic characteristics of the network configurations corresponding to the obtained set of potentially non-dominated solutions, before presenting our conclusions in Sect. 5.

Problem formulation
In this paper, we consider a water distribution network with n n demand nodes or junctions, n 0 water sources and n p links, among which n V existing pressure control valves, n CNV candidate links for valve installation and n CNP candidate pipes for installation. Pumps are not considered and we assume known the location of existing valves, CNVs and CNPs, as well as the length, diameter and material of CNPs. The problem consists in simultaneously selecting, from the predefined sets of CNVs and CNPs, the candidate valves and pipes to add to the network, and optimizing the control of existing and new valves.
We model the network as a directed graph with n n + n 0 nodes (vertices) and n p links (edges). A positive flow through a link is defined as directed from start to end node. We refer to V as the set of network nodes (vertices) and E as the set of network links (edges). Index subsets E CNV and E CNP correspond to CNV and CNP locations, respectively. We consider n t different time steps of network operation and, for k = 1, . . . , n t , the vector of known nodal demands is given by d k ∈ R n n and known hydraulic heads at water sources by h k 0 ∈ R n 0 . Let ζ and h min ∈ R n n be the vectors of known elevations and minimum allowed hydraulic heads at network nodes, respectively. For the case study networks considered here, we set, for i = 1, . . . , n n , where A 12 ∈ R n p ×n n and A 10 ∈ R n p ×n 0 are the link node incidence matrices for demand nodes and water sources respectively. Friction head losses associated with pipe flows q k = [q k 1 . . . q k n p ] T , represented by the vector φ(q k ) = [φ 1 (q k 1 ) . . . φ n p (q k n p )] T , are commonly described either by the Darcy-Weisbach (D-W) or Hazen-Williams (H-W) equations. As neither equation is smooth, it is convenient to use a quadratic approximation of the formulae over the operational range of flows in each link. This produces a pair of positive coefficients (a j , b j ) ∈ (R + ) 2 for each pipe j in the network and the vector of friction head loss at time step k is then given by . This is a common approach for the formulation of head loss constraints in WDN optimization literature (Dai and Li 2014;Eck and Mevissen 2015;Pecci et al. 2017b;Zamzam et al. 2019). Moreover, the results presented in "Appendix F" suggest that the use of a quadratic approximation does not affect the quality of the solutions computed for the case study networks considered in this work. The error in the objective values, in particular, does not exceed 2%.
To isolate the non-convex friction head loss term, we introduce the vector of auxiliary variables θ k ∈ R n p and constraint (1a) becomes: where (2b) is the non-convex friction loss equality constraint, while constraints (1b) and (2a) are linear. Constraint (3) ensures water is flowing out from the network sources. We also define vectors q k U , q k L , h k U , h k L , θ k U , θ k L , η k U and η k L ∈ R n p of upper and lower bounds on variables q k , h k , θ k and η k , respectively, such that and η k L are defined for pipes, valves and CNVs, and CNPs by the expressions (A.1), (A.2) and (A.3) in "Appendix A".
To model the status of the candidate valves and links, we introduce the binary variable z ∈ R n CNV +n CNP : for j ∈ E CNV ∪ E CNP , the status of j is represented by z l j , l j ∈ {1 . . . n CNV + n CNP }. For j ∈ E CNV , either z l j = 1 and the candidate link j is selected for valve installation and modeled as a pressure control valve, or z l j = 0 and link j is unselected for valve installation, i.e. modeled as an open pipe. For j ∈ E CNP , either z l j = 1 and the candidate location j is selected for pipe installation and modeled as an open pipe, or z l j = 0 and j is modeled as a closed pipe. Additional head losses across unselected CNVs and flows across closed pipes are forced to zero, while hydraulic heads at the start and end nodes of closed pipes are decoupled. This is enforced by the following big-M constraints on variables q k , θ k and η k : . Further constraints on binary variables z, to prevent specific combinations of CNVs and CNPs from being selected for instance, can be enforced with where A ∈ R n A ×(n CNV +n CNP ) is a constant matrix of binary inequalities representing combinatorial constraints on the feasible sets of candidate valves and pipes and e ∈ R n A is a n A by 1 vector with fixed entries. Finally, we define the vector of decision variables x = (z, q, h, η, θ), where q = (q k ) k=1,...,n t , h = (h k ) k=1,...,n t , η = (η k ) k=1,...,n t , θ = (θ k ) k=1,...,n t , and the set F(Q) ⊂ {0, 1} n CNV +n CNP × R n p × R n n × R n p × R n p given by (1b), (2a), (3), (4), (5) and (6), which depends on the bounds Q:={q k L , q k U } k=1,...,n t .

Objective functions
In this work, we consider the bi-objective problem of the optimal design-for-control of water distribution networks, with the objectives of minimizing pressure induced background leakage and maximizing resilience. Pressure induced background leakage Water utilities in the UK are incentivized, by means of financial penalties, to reduce leakage throughout their systems. Background leakage, in particular, which represents a large part of total water losses, depends on pressure levels in the distribution network (Van Zyl and Clayton 2007). Here, we use average zone pressure (Wright et al. 2015) as a surrogate measure for pressure induced leak flow: where w ∈ R n n and W ∈ R are the vector of AZP weights and the AZP normalization factor, respectively. For i ∈ V , we define J i , the set of all pipes incident to node i. Then w i = j∈J i L j /2 and W = n n i=1 w i , where L ∈ R n p is the vector of pipe lengths.
Resilience Another key performance commitment concerns the reduction of interruptions of supply. The resilience of a WDN, which refers to its ability to maintain continuous customer supply under failure conditions, depends on the availability of alternative supply paths (topological redundancy) and pressure surplus (energy redundancy). The resilience of a WDN can be measured directly as its performance under failure conditions. This approach however suffers from combinatorial limitations, as it requires the simulation of a wide range of failure events. In comparison, surrogate measures of resilience depend only on the structure of the network, or its performance under normal operating conditions. A commonly adopted surrogate measures of resilience is the resilience index I r (Todini 2000), which is defined as the ratio of surplus energy to total energy input measured under normal operating conditions: Both the numerator and denominator in (8) are affine functions of the vectors of hydraulic heads h and flows q respectively: in the general case, I r is a linear fractional function of the problem variables. Moreover, for all x = (z, q, h, η, θ) satisfying (1b), (2), (3), (4), (5), (6) and z ∈ {0; 1} n CNV +n CNP , the denominator of I r (x) is strictly positive and 0 ≤ I r (x) ≤ 1-see "Appendix B". We observe, from (7) and (8), that both AZP and I r are increasing functions of the hydraulic heads at the network nodes. As a result, minimizing pressure induced background leakage and maximizing WDN resilience are conflicting objectives. The optimal DfC problem for the joint minimization of AZP and maximization of I r (i.e. minimization of −I r ) can be formulated as a bi-objective (non-convex) mixed-integer non-linear program: minimize x= (z,q,h,η,θ) (z,q,h,η,θ) where f 1 and f 2 represent either objective functions AZP or −I r , (2b) is the non-convex potential-flow coupling constraint and F(Q) is the set defined by the linear constraints (1b), (2a), (3), (4), (5) and (6). Solving BOMINLP(Q) requires identifying the best compromises between the objective functions f 1 and f 2 , or non-dominated set: if x * is an efficient solution of BOMINLP(Q). The set of non-dominated solutions of BOMINLP(Q) is called the non-dominated set, or Pareto front of BOMINLP(Q).
In the following section, we present a method to approximate the Pareto front of BOMINLP(Q).

Solution method: -constraints and sBB
The bi-objective design-for-control problem aims to identify trade-off solutions representing efficient compromises between conflicting objectives, as no single solution jointly optimizes all objectives. In particular, we are interested in deterministic methods for computing potentially non-dominated solutions of BOMINLP(Q) with guarantees of global non-dominance. In Fernández and Tóth (2007); Ehrgott and Gandibleux (2007), the authors present scalarization-based approaches to compute outer approximations of the Pareto set of bi-objective non-linear problems (NLP) and the Pareto front of bi-objective mixedinteger linear problems (MILP). Previous studies (Cacchiani and D'Ambrosio 2017;De Santis et al. 2019;Niebling and Eichfelder 2019) obtain similar guarantees of global efficiency and/or non-dominance by extending the early work of Mavrotas and Diakoulaki (2005); Fernández and Tóth (2007) on multi-objective branch-andbound to compute supersets of the Pareto sets and/or Pareto fronts of a wider range of multi-objective optimization problems, including convex MINLPs and continuous NLPs with non-convex objective functions. None of these approaches has, however, considered BOMINLPs with non-convex constraints.
This work investigates the application of scalarization, together with global optimization methods, to approximate the Pareto front of BOMINLP(Q). We implement the method of -constraints (Ehrgott 2006) to scalarize the bi-objective problem (see Sect. 3.1) and use spatial branch-and-bound to compute upper and lower bounds on the optimal values of the resulting sequence of parametrized single-objective problems (see Sect. 3.2). We show that, using this approach, we can compute a superset of the Pareto front of BOMINLP(Q), i.e. potentially non-dominated solutions as well as guarantees of their global non-dominance (see Sect. 3.3).

Scalarization using the method of -constraints
Scalarization approaches convert a multi-objective problem into a sequence of parametrized single-objective problems. By varying parameter values, different solutions from the efficient and/or non-dominated sets are produced. For a review of scalarization methods, see Ehrgott (2006). In this work, we focus in particular on the method of -constraints. In -constraints, one objective function is selected to be minimized while upper bounds are enforced on other objective function. We scalarize BOMINLP(Q) by introducing an -constraint on the second objective function f 2 and reformulate it into a sequence of n (non-convex) parametrized MINLPs: minimize x= (z,q,h,η,θ) where j ∈ [0, 1], for all j = 1, . . . , n , and the anchor point x * f i , i = 1, 2, is an optimal solution of the single-objective non-convex problem MINLP f i (Q): minimize x= (z,q,h,η,θ) As a result, the quality of the approximation of the Pareto front depends on the order of the objective functions in BOMINLP(Q) and choice of scalarization parameters j , j = 1, . . . , n . We discuss this aspect, as well as the scalability challenges associated with the extension of the proposed method to multi-objective problems with more than two objectives, in Sect. 4.3.
The use of the -constraints scalarization method has two main benefits. First, we show (see Sect. 3.2.2) that relaxing the potential-flow coupling constraints in MINLP f i (Q) or MINLP f 1 , f 2 , j (Q) produces a problem which can be equivalently reformulated into an mixed-integer linear program (MILP). This provides a convenient mean to compute tight lower bounds throughout the sBB procedure. In addition, we show that the sets of lower and upper bounds on the optimal value of the single-objective problems MINLP f i (Q), i = 1, 2, and MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , computed by sBB define the boundary of a superset of the Pareto front of BOMINLP(Q)-see Sect. 3.3.

Spatial branch-and-bound
In this section, we present a tailored spatial branch-and-bound (sBB) algorithm to compute good quality feasible solutions, with certified global optimality bounds, for a non-convex mixed-integer non-linear program MINLP(Q), of the form MINLP f i (Q) or MINLP f 1 , f 2 , j (Q), defined on the initial bound set Q = Q init . The sBB method described in Algorithm 1 recursively divides the initial bound set Q init into increasingly smaller sets Q stored in the working list Q, creating a tree of descendant subproblems MINLP(Q). Lower bounds on the optimal value of the subproblems are obtained by solving mixed-integer linear relaxations. In particular, we propose a tailored lower bounding procedure, which relies on an equivalent linear reformulation of the resilience index −I r . In comparison, upper bounds correspond to locally optimal solutions of the NLPs obtained by fixing the values of the integer variables in the non-convex MINLP subproblems. Global optimality is then certifiably reached when the relative optimality gap between upper bound and lower bound is smaller than a specified tolerance, δ (here we set δ = 1e −2 -see Algorithm 1).
For i = 1, 2 (resp. j = 1, . . . , n ), the sBB algorithm returns a feasible solution x f i (resp. x f 1 , f 2 , j ), an upper bound UB f i (resp. UB f 1 , f 2 , j ) and a lower bound LB f i (resp. LB f 1 , f 2 , j ) on the optimal value of MINLP f i (Q init ) (resp. MINLP f 1 , f 2 , j (Q init )), and the set of nodes Q f i (resp. Q f 1 , f 2 , j ) explored during the search. Finally, to compute the right-hand side of the -constraint in MINLP

Selection strategy, branching strategy and bound tightening
After initializing the sBB search, and as long as the relative optimality gap is greater than the set tolerance δ and the time limit T max has not been exceeded, the set of bounds Q s corresponding to the subproblem MINLP(Q s ) with the smallest lower bound L(Q s ) is removed from the working list Q and divided. The resulting descendant bound sets are then reduced using an Optimality-Based Bound Tightening (OBBT) methodsee Puranik and Sahinidis (2017) for a review of domain reduction techniquesand stored in the working list Q. Although OBBT is associated with a significant computational overhead, the implementation of domain reduction methods throughout the sBB search was shown to improve final bounds on the optimal value of non-convex problems more efficiently than branching alone (Castro and Grossmann 2014;Ulusoy et al. 2019). In this work, we implement the domain reduction procedure described in Algorithm 1 in Pecci et al. (2018)

Lower bounding
To obtain lower bounds L(Q) on the optimal value of the subproblems MINLP(Q) generated throughout the sBB search, we relax the subproblems by computing linear outer approximations of the non-convex potential-flow coupling constraint (2b) and reformulate the resilience index −I r to produce relaxed MILPs. The solution of the resulting linear problems provides lower bounds on the optimal value of subproblems MINLP(Q).
Relaxation of the potential-flow coupling constraints To relax MINLP(Q), we compute a polyhedral envelope of the potential-flow coupling constraint (2b) implementing a method for monomials of odd degree developed by Liberti and Pantelides (2003). The original potential-flow coupling constraint (2b) is replaced by the linear constraint: where matrices R k and E k ∈ R m×n p and vector r k ∈ R m depend on the bounds Q and number of linearizations used to compute the linear approximation of the convex envelope of constraint (2b). For the sake of brevity, we omit a detailed description of these relaxations, which can be found in Pecci et al. (2018); Ulusoy et al. (2019). The relaxation of constraint (2b) in MINLP(Q) results either directly in an MILP or, for In this case, we propose a linear reformulation of the objective function −I r for (10). Linear reformulation of −I r In Charnes and Cooper (1962), the authors propose a linear reformulation for the continuous non-linear problem (NLP) of maximizing a linear fractional function subject to linear constraints. Similarly, we reformulate the linear fractional objective function −I r in (10), and introduce new variablest ∈ R * + andw ∈ R n CNV +n CNP to produce the following MILP: ), (C.6), (C.8), (C.9), (C.10) and (C.11) in "Appendix C". We can then show that is a feasible solution of (10), where z =z, q =q/t, h =h/t, η =η/t and θ =θ/t. Proof See Proposition 1 in "Appendix C".
In addition, it follows from the above proposition that Corollary 1 If x is a global optimal solution of (10), thenx, as given by Proposition 1, is a global optimal solution of (11), and vice versa.
Proof See Corollary 1 in "Appendix C".
For the sake of brevity, we refer the reader to the "Appendix" for the full description of the MILP reformulation and technical proofs of Proposition 1 and Corollary 1. Hence, to solve the linear fractional problem (10), it is necessary and sufficient to solve the corresponding linear reformulation (11). Note that a similar reformulation can be applied to problems MINLP −I r ,AZ P, j , j = 1, . . . , n .
As a result, we are able to compute lower bounds L(Q) on the optimal value of subproblems MINLP(Q) throughout the sBB search by relaxing potential-flow flow constraints (2b), reformulating the resilience index −I r if necessary, and applying an MILP solver to the resulting linear problem.

Upper bounding
Finally, upper bounds U(Q) on the optimal value of subproblems MINLP(Q) are computed throughout the sBB search by fixing the values of the binary variables z in MINLP(Q) to the solution of the corresponding lower bounding MILP relaxations and computing a locally optimal solution to the resulting non-convex NLPs. This is done directly using an NLP solver.

Note on the initialization of the sBB search for MINLP
returned by the sBB procedure. Then, if j+1 is greater than or equal to j , since x f 2 is the best feasible solution found for MINLP f 2 (Q), we have that Since all other constraints are the same, As a result, we define the sequence of parameters ( j ) j=1,...,n such that j ≤ j+1 , for all j = 1, . . . , n − 1, and we apply sBB to the sequence of parametrized problems MINLP f 1 , f 2 , j (Q) for j = 1, . . . , n , in increasing order of j , initializing the upper bound UB of every problem MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , at the best feasible solution of the previous problem in the sequence, i.e. UB ←

Approximation of the Pareto front
In this section, we show that, combining the methods described in Sects. 3.1 and 3.2, we obtain a superset of the Pareto front of BOMINLP(Q). We define L FS f 1 , f 2 as the image in the objective space of the set of feasible solutions of MINLP f i (Q), i = 1, 2, and MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , returned by Algorithm 1: Next, we perform a pairwise comparison of the elements of L FS f 1 , f 2 , called a Pareto filter, to eliminate dominated solutions. The resulting subset of potentially non-dominated solutions, L NDS f 1 , f 2 , is then given by: For a detailed description of the Pareto filter algorithm, see Section 4.1 in Messac et al. (2003). We also define the sets L LBS f 1 , f 2 , T f 1 , f 2 and U f 1 , f 2 : We refer the reader to Figs. 1 and 3 for an illustration of the sets L LBS f 1 , f 2 , L NDS f 1 , f 2 , T f 1 , f 2 and U f 1 , f 2 . Then, the following Proposition holds.
We distinguish three cases.
Asp is a feasible solution of BOMINLP(Q), this contradicts LB f 2 lower bound on the optimal value of f 2 .
Proposition 2 implies that the set R 2 \ T f 1 , f 2 ∪ U f 1 , f 2 , defined by the best feasible solutions and lower bounds returned for MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , by Algorithm 1, is a superset of the non-dominated set of BOMINLP(Q). In particular, as Proposition 2 holds for both f 1 = AZP and f 2 = −I r , and f 1 = −I r and f 2 = AZP, the sets T , U , L NDS and L LBS , given by define a tighter superset of the Pareto front of BOMINLP(Q), corresponding to the intersection of the supersets Similarly, the sets of explored nodes returned by Algorithm 1 for problems MINLP f i (Q), i = 1, 2, and MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , define a superset of the efficient set of BOMINLP(Q)-see "Appendix D". The quality of the global efficiency bounds on solutions x i , i = 1, 2, and x j , j = 1, . . . , n then depend on the size of the sets F(Q) in the tree of explored nodes returned by Algorithm 1.

Application and results
We investigate the optimal design-for-control of two WDNs, namely Net25 and pescara. For each network, we consider the bi-objective optimization problem BOMINLP(Q), with the objectives of minimizing average zone pressure and maximizing network resilience, and apply the method described in Sect. 3. As the scalarization method of -constraints is sensitive to the order of objectives functions, we propose to solve both sequences of parametrized problems MINLP AZP,−I r , j and MINLP −I r ,AZP, j , for j = 1, . . . , n . We consider different values of n , and for j = 1, . . . , n , we define j = j/n .
For both Net25 and pescara, we first solve the anchor point problems MINLP AZP (Q) and MINLP −I r (Q) and compare the performance of Algorithm 1 on these instances to the off-the-shelf solvers SCIP and BARON. For each solution method, the time limit is fixed to 5 min (T max = 5 min), including pre-processing steps. In particular, all bound tightening steps in Algorithm 1 are performed within the fixed time limit. The results of the comparison, presented in Tables 3, 4, 5 and 6 in "Appendix E", show that the spatial branch-and-bound algorithm outperforms state-of-the-art global optimization solvers on all instances of MINLP AZP (Q) and MINLP −I r (Q). Moreover, a posteriori analysis of the convergence of Algorithm 1 for the solution of the anchor point problems MINLP AZP (Q) and MINLP −I r (Q) shows that both upper and lower bounds improve rapidly at first. Past 5 min, however, the improvement in the optimality gap slows down significantly-see "Appendix E". As a result, in the rest of this section, we implement Algorithm 1, with T max = 5 min, for the solution of the sequence of scalarized problems MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , and define the sets T f 1 , f 2 and U f 1 , f 2 based on the best feasible solutions and lower bounds returned. Within Algorithm 1, GUROBI (v8.1) (Gurobi Optimization 2017) is used to solve the lower bounding and bound tightening MILPs, while upper bounds are computed with IPOPT (v3.12.9) (Wachter and Biegler 2006) using the Matlab interface (Currie et al. 2012). Taking advantage of the sparsity of the NLPs, we directly provide pre-computed gradients and Jacobians to the solver. The tolerance and time limit in Algorithm 1 are set to 10 −2 and 5 min. All computations are performed using MATLAB 2018b for Windows, on a 2.26GHz Intel ® Xeon ® CNPU E5520 with 8 cores.

Net25
Sectorization is a common pressure and leakage management practice which consists in closing off links in operational networks using boundary valves (BV) to create isolated supply areas, called district metered areas (DMA). Although sectorization reduces the topological redundancy of WDNs, improvements in resilience can be achieved at minimal cost in sectorized networks by re-opening closed BVs and aggregating DMAs. We illustrate this on the modified network Net25 from Ulusoy et al. (2019), where nodes with positive demand represent DMAs. We formulate the optimal aggregation of isolated DMAs (represented in Net25 by nodes with positive demand) as a DfC problem, where the set of CNPs correspond to closed BVs and DMA aggre-gation is limited to pairing, to preserve the structure of the network. The selection of CNPs for installation can then lead to the aggregation of a DMA with at most one other DMA, and the pairing of two DMAs excludes pairing opportunities with other neighboring DMAs.
The network model, provided in the supplementary material, consists of 3 reservoirs, 17 nodes and 26 links (15 pipes, 3 existing valves and 8 closed BVs, representing CNPs for installation). Friction head losses across network links are modelled using the H-W equation, and we consider 24 time steps. The maximum allowed velocity is fixed to 2 m/s for all network pipes (based on preliminary hydraulic simulations of the network) and a quadratic approximation φ(·) of the H-W formula is adopted to model head loss across the pipes (see Sect. 2). We set the pressure requirement at nodes with positive demand to 10 m and model the pairing limit on DMA aggregation with constraint (6). For the details of the computation of matrix A corresponding to the constraint of DMA pairing, see Appendix A in Ulusoy et al. (2019). Note that the combinatorial constraint limits the number of CNPs which can be simultaneously selected for installation.
Anchor points The computation of the anchor points of BOMINLP(Q) for Net25 Solution of MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), j = 1, . . . , n We apply Algorithm 1 to the sequences of parametrized problems MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), for j = 1, . . . , n and n = 10, 20, 100. After discarding the dominated solutions, we obtain the sets L NDS AZP,−I r and L NDS −I r ,AZP of potentially non-dominated solutions, represented in Fig. 1a and b for n = 20. On Fig. 1a and b, we also represent the sets L LBS AZP,−I r , T AZP,−I r and U AZP,−I r , and L LBS −I r ,AZP , T −I r ,AZP and U −I r ,AZP , respectively (see Sect. 3.3). The staircase-shaped areas R 2 \ T AZP,−I r ∪ U AZP,−I r and R 2 \ T −I r ,AZP ∪ U −I r ,AZP then represent the supersets of the Pareto front of BOMINLP(Q) obtained with each scalarization approach. The differences between the two scalarization approaches are discussed in greater detail in Sect. 4.3. As observed in Sect. 3.3, a tighter superset of the Pareto front is obtained by taking the intersection of R 2 \ T AZP,−I r ∪ U AZP,−I r and suggests that the potentially non-dominated solutions obtained for MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), j = 1, . . . , 20, provide a good approximation of the true nondominated set of BOMINLP(Q).
Next, we discuss the trade-off between the objectives of AZP and I r , their evolution along the Pareto front and the hydraulic characteristics of the configurations corresponding to the potentially non-dominated solutions represented in Fig. 1c. Potentially non-dominated solutions p ∈ L NDS with lower AZP values, represented in the top left corner in Fig. 1c, correspond to more sectorized network designs, resulting in greater friction head loss, and the implementation of aggressive pressure reduction. Such configurations minimize surplus pressure at customer nodes, resulting in lower AZP and I r values. Moving down along the Pareto front in Fig. 1c towards more resilient network configurations in the bottom right corner, pressure reduction is relaxed and more CNPs are selected for installation. This results in better flow distribution, reduced friction head loss across the network, excess pressure at customer nodes and, eventually, greater AZP and I r values.
Finally, we focus on the geometry of the Pareto front of BOMINLP(Q), which is expected to be disconnected and/or composed of different branches. Figure 1d represents the potentially non-dominated solutions p ∈ L NDS from Fig. 1c, where solutions corresponding to different network designs are represented by different colors. Note that several points in L NDS correspond to the same network design: L NDS represents the union of 6 fronts, corresponding to 6 different network designs, i.e. sets of CNPs candidate valve locations candidate pipe locations  Bragalli et al. (2008). The original network model is available at http://www.or.deis.unibo.it/research_pages/ORinstances/ORinstances.htm selected for addition. As a result, when selecting the most suitable network design, the decision-maker is not limited to a single trade-off point, but rather a set of trade-offs corresponding to different PCV control settings. For example, if the decision maker's priority is minimizing AZP, then design 5, which corresponds to opening a single BV in the network, represents a better trade-off than design 1, which results in 8% larger minimum AZP values. On the other hand, design 1, which corresponds to the creation of 3 pairs of DMAs, results in more resilient network configurations.

pescara
The pescara model is a reduced version of the WDN of an Italian medium-size city, studied by Bragalli et al. (2008). The original network model counts 71 junctions, 3 reservoirs, 99 pipes and is defined on a single time step. Friction losses are modeled according to the H-W equation and maximum allowed velocities are fixed to 3 m/s for all network links. We investigate the problem of expanding the network by installing new pressure control valves and pipes in predefined candidate locations (see Fig. 2 and supplementary material), determined based on engineering judgement. For each candidate link location, we consider 3 CNPs with diameters 100, 150 and 200 mm. The material of the CNPs is cast iron, as for existing pipes, and we include a combinatorial constraint in the optimization problem in the form (6) to limit the selection of CNPs to 1 per location.
Solution of MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), j = 1, . . . , n We apply Algorithm 1 to the sequences of parametrized problems MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), for j = 1, . . . , n and n = 10, 20, 100. For each approach, dominated solutions p ∈ L FS AZP,−I r or L FS −I r ,AZP are discarded. The resulting subsets of potentially non-dominated solutions L NDS AZP,−I r and L NDS −I r ,AZP obtained for n = 20 are represented in Fig. 3a and b, along with the sets L LBS AZP,−I r , T AZP,−I r and U AZP,−I r , and L LBS −I r ,AZP , T −I r ,AZP , and U −I r ,AZP . We comment on the differences between the two approaches in Sect. 4.3. Moreover, we represent, on Fig. 3c, the sets T = T AZP,−I r ∪ T −I r ,AZP and U = U AZP,−I r ∪ U −I r ,AZP obtained for n = 20. The resulting staircase-shaped set R 2 \ (T ∪ U ), which corresponds to the intersection of the sets R 2 \ T AZP,−I r ∪ U AZP,−I r and R 2 \ T −I r ,AZP ∪ U −I r ,AZP represented in Fig. 3a and b, provides a tighter superset of the set of non-dominated solutions of BOMINLP(Q).
Potentially non-dominated solutions with lower AZP values in Fig. 3c correspond to the installation of fewer or smaller CNPs and more CNVs, and the implementation of aggressive pressure reduction, which result in greater head loss across the network, reduced surplus pressure at customer nodes and lower I r values. In particular, solutions in the top left corner of Fig. 3c require the installation of all CNVs. Moving down along the Pareto front in Fig. 3c towards the bottom right corner, more resilient potentially non-dominated solutions p ∈ L NDS correspond to the installation of fewer CNVs and CNPs with larger diameters, which increase the number and capacity of alternative customer supply paths. For network configurations corresponding to solutions p ∈ L NDS in the intermediate segment of the Pareto front, only the CNV located at the outlet of the source of greatest fixed hydraulic head h 0 is selected for installation (and the implementation of pressure reduction). The effect is to increase the relative importance of other supply sources, improving flow distribution across the network. As the energy efficiency of the network grows, so does the pressure surplus at network nodes, increasing AZP and the resilience index, I r .
Design variables in the optimal design-for-control problem of pescara include the placement of PCVs, and 3 different diameters to choose from, for each candidate pipe. As this results in a higher number of possible design combinations than for Net25, it is expected that fewer potentially non-dominated solutions p ∈ L NDS correspond to the same network design. From the perspective of the decision-maker, it might be interesting to take more objectives (e.g. cost) into account when deciding between two design solutions offering similar trade-offs between AZP and −I r .

Scalarization approach
Scalarization methods, including the method of -constraints, are known to be sensitive to the order of the objective functions (Logist and Van Impe 2012) and the choice of scalarization parameters (Kim and De Weck 2005;Burachik et al. 2017). In the method we propose in Sect. 3, the first objective function f 1 is chosen to be minimized while parametrized upper bound constraints are enforced on the value of f 2 . As a result, the quality of the Pareto front approximation depends on the order of the objective functions-see Figs. 1 and 3. In this section, we discuss the differences between the results presented in Tables 1 and 2, obtained for Net25 and pescara with both sequences of parametrized problems MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q),  for j = 1, . . . , n and n = 10, 20, 100. We also comment on the choice of the uniform sequence of scalarization parameters j = j/n , for j = 1, . . . , n .

Order of objective functions
To Tables 1 and 2 show that the quality of the Pareto front approximation obtained by solving the sequence of scalarized problems MINLP f 1 , f 2 , j (Q), j = 1, . . . , n , depends on the order of the objective functions, AZP and −I r . In particular, the solutions of the sequence of parametrized problems MINLP AZP,−I r , j (Q), j = 1, . . . , n , define a tighter superset of the non-dominated set of BOMINLP(Q) over the range of values taken by AZP and −I r , as n(L NDS AZP,−I r ) ≥ n(L NDS −I r ,AZP ) and n(L LBS AZP,−I r ) ≥ n(L LBS −I r ,AZP ), for both networks and for all values of n . This is illustrated, for pescara, by the difference in the width of the supersets R 2 \ T AZP,−I r ∪ U AZP,−I r and R 2 \ T −I r ,AZP ∪ U −I r ,AZP represented in Fig. 3a and b.
The difference between the obtained Pareto front approximations can be explained by the number of sBB iterations completed by Algorithm 1 for the scalarized problems MINLP −I r ,AZP, j (Q) and MINLP AZP,−I r , j (Q). For instance, for pescara and n = 20, a single sBB iteration takes on average 5 times longer to complete for problems MINLP −I r ,AZP, j (Q) than for MINLP AZP,−I r , j (Q). Therefore, half as many sBB iterations, on average, are completed within the fixed time limit (T max = 5 min) for each parametrized problem in the former sequence, resulting in larger gaps. By fixing a limit on the maximum number of sBB iterations completed (instead of computational time) in Algorithm 1, the results obtained with the different scalarization approaches are more balanced-see "Appendix G". However, the results presented in Tables 1 and 2 show that it is more efficient to solve the sequence of parametrized problems (MINLP AZP,−I r , j (Q)) j=1,...,n to produce a tight approximation of the Pareto front of BOMINLP(Q), within a prescribed tolerance or time limit.

Note on the choice of scalarization parameters j
The results presented in this section for the optimal design-for-control of Net25 and pescara were obtained by solving the sequences of parametrized problems MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), for j = 1, . . . , n , using a predefined, uniformly distributed set of scalarization parameters j = j/n .
Depending on the geometry of the Pareto front, traditional scalarization approaches relying on predefined distributions of parameters might fail to explore regions of the non-dominated set, resulting in redundant points and poorly distributed solutions on the Pareto front. Traditional weighted sum or -constraint methods for instance struggle to obtain a representative approximation of Pareto fronts with steep segments (i.e. knees), while the normalized boundary intersection (NBI) method (Das and Dennis 1998) requires the non-dominated set to be simply connected (i.e. not composed of disconnected branches) to guarantee finding a Pareto optimal solution for all parametrized problems (Burachik et al. 2017). Previous works have shown that these limitations could be successfully addressed by using adaptive choices of scalarization parameters (Kim and De Weck 2005) or grid generation techniques (Burachik et al. 2017(Burachik et al. , 2019, to focus on unexplored regions of the Pareto front and obtain representative approximations of non-dominated sets with irregular geometries. In this case, however, we observe (see Figs. 1 and 3) that the choice of scalarization parameters j results in an even spread of potentially non-dominated solutions L NDS and lower bounds L LBS , providing a good approximation of the whole non-dominated set of BOMINLP(Q). As a result, in this work, we do not further investigate adaptive methods for the generation of scalarization parameters j . The implementation of grid generation techniques could however be investigated in future work, in particular for multi-objective problems with more than two objectives.

Conclusion
To the best of our knowledge, this study is the first to consider the bi-objective optimization for the minimization of pressure induced background leakage and maximization of resilience in water distribution networks. In particular, we investigate the bi-objective optimal design-for-control of existing WDNs, which consists in simultaneously identifying new valves and/or pipes to add from a predefined set of candidates and optimizing the control of new and existing valves in the network. The mathematical formulation of the problem results in a non-convex BOMINLP.
The exact identification or approximation of the Pareto front of non-convex BOMINLPs with global bounds are challenging tasks, which have not been previously addressed in the literature. In this work, we present a method to compute an approximation of the non-dominated set of the considered bi-objective DfC problem with global non-dominance bounds. The BOMINLP is first scalarized, using the method of -constraints, into a sequence of parametrized single-objective non-convex MINLPs. Spatial branch-and-bound is then applied to compute feasible solutions of the individual non-convex MINLPs with certified global optimality bounds. In particular, we propose a novel linear reformulation of the resilience index, I r , which allows us to derive MILP relaxations of the parametrized MINLPs and, as a result, lower bounds on the optimal value of the parametrized problems.
We show that the implementation of sBB in combination with -constraints returns a set of potentially non-dominated solutions as well as a superset of the Pareto front of the original BOMINLP based on the feasible solutions and guarantees of global optimality returned by sBB for the individual parametrized MINLPs. In particular, the superset of the Pareto front is a guarantee of global non-dominance of the potentially non-dominated solutions found. The proposed method and results can also be extended, in future work, to study multi-objective design-for-control problems with additional objective functions (e.g. cost, water quality). The implementation of grid generation techniques might be required to address the scalability challenges associated with the method of -constraints for multi-objective problems with more than two objectives.
We apply this method to solve the bi-objective DfC problem for Net25 and pescara, and compute a representative set of trade-offs between the objectives of pressure induced background leakage minimization and resilience maximization. Network configurations corresponding to the installation of more candidate PCVs and the implementation of aggressive pressure reduction minimize pressure surplus at network nodes and result in lower AZP and I r values. The installation of more (and larger) candidate pipes, on the other hand, improves network efficiency by increasing the number (and capacity) of alternative supply paths, and results in greater AZP and resilience index values. The results provide insights into the complex structure of the Pareto fronts of non-convex BOMINLPs and show that the quality of Pareto front approximations obtained by solving a sequence of -parametrized problems depends on the order of the objective functions.
Acknowledgements The authors acknowledge the financial support of EPSRC (EP/P004229/1, Dynamically Adaptive and Resilient Water Supply Networks for a Sustainable Future), EPSRC CDT in Sustainable Civil Engineering and Suez for conducting this research. The authors would also like to thank the anonymous reviewers for their valuable comments and suggestions, resulting in the inclusion of additional proofs and numerical analysis, as well as the correction of an error in the numerical implementation of our sBB algorithm.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Problem formulation (Sect. 2.1)
The formulation of the multi-objective optimal DfC problem MOMINLP f 1 , f 2 (Q) relies on the set F(Q) ⊂ {0; 1} n CNV +n CNP × R n p × R n n × R n p × R n p , defined by constraints (1b), (2a), (3), (4), (5) and (6)-see Sect. 2 of the manuscript. While constraints (1b), (2a), (3), (6) are described in the manuscript, constraints (4) and (5) rely on the definition of (constant) vectors of upper and lower bounds q k U , q k L , θ k U , θ k L , η k U and η k L ∈ R n p on variables q k , θ k and η k and matrices and vectors of big-M constraints Q k U , Q k L , k U , k L , N k U and N k L ∈ R n p ×(n CNV +n CNP ) , andη k L andη k U ∈ R n p which are defined below.
Let v max ∈ R n p be the vector of maximum allowed velocities across the network links, and q max = diag([π D 2 1 /4 . . . π D 2 n p /4] T )v max the corresponding vector of maximum flows across network links, where D j is the diameter of link j. Recall where h k 0 is the vector of known hydraulic heads at network sources, and h k L = h min , the vector of minimum allowed hydraulic heads at network nodes. Index subsets E P , E V , E CNV and E CNP , correspond to network pipes, existing valves, CNV and CNP locations, respectively. Observe that E P , E V , E CNV and E CNP represent a partition of E. Moreover, the notation i 1 j − → i 2 represents that a link j positively connects node i 1 (start node) to i 2 (end node).

A.1 Upper and lower bound vectors
We first define the (constant) vectors of upper and lower bounds q k U , q k L , θ k U , θ k L , η k U and η k L ∈ R n p on variables q k , θ k and η k . Pipes Both positive and negative flows (and, consequently, friction head losses) are allowed across network links j ∈ E P . Flows and friction losses are only limited in absolute value by the maximum allowed velocities. We assume only friction head loss occurs across network pipes, so η k j = 0, ∀ j ∈ E P , k = 1, . . . , n t . This is reflected in the upper and lower bounds on the values of q, θ and η associated with the network pipes. For all k = 1, . . . , n t : Valves and candidate valves In both networks studied in Sect. 4 of the main manuscript, we consider the direction of operation of CNVs, which are oriented in the positive flow direction, to be known. Therefore, we constrain flows across valves and CNVs to be positive (directed from start to end node). This constraint is also reflected by the upper and lower bounds of the friction loss and additional head loss across the valves. For all k = 1, . . . , n t : Candidate pipes Finally, the values of individual elements of constant vectors q k U , q k L , θ k U , θ k L , η k U and η k L for all k = 1, . . . , n t are defined as follows for CNPs:

A.2 Big-M constraints
We also define the n p by n CNV + n CNP diagonal matrices Q k U , Q k L , k U , k L , N k U and N k L of big-M constraints (5), for all k = 1, . . . , n t : And vectorsη k L andη k U by:

Appendix B: Resilience objective function (Sect. 2.2)
In Sect. 2.2 of the main manuscript, we introduce I r , the surrogate measure of WDN resilience defined by (8): Here, we make observations about the values taken by I r (x) and its denominator. In particular, we show that the denominator of I r (x) is strictly positive and that possible values of I r (x) range from 0 to 1 for all vectors x = (z, q, h, η, θ) satisfying (1b), (2), (3), (4), (5), (6) and z ∈ {0; 1} n CNV +n CNP , where q = (q k ) k=1,...,n t , h = (h k ) k=1,...,n t , η = (η k ) k=1,...,n t , θ = (θ k ) k=1,...,n t . By multiplying Eq. (2a) by q k and rearranging, we obtain: Then, Eq. (1b) yields and, by rearranging, For all k = 1, . . . , n t , we observe that (q k ) T η k ≥ 0, because η k j = 0 for j ∈ E P , q k j ≥ 0 and η k j ≥ 0 for j ∈ E V , and q k j = 0 or η k j = 0 for j ∈ E CNP . Moreover, unless d k = 0 or a j = b j = 0 for all j ∈ E P , there exist j ∈ E P such that As we consider all customer nodes to have positive demands, ∀i ∈ V , k = 1, . . . , n t , d k i ≥ 0, and, for all feasible solutions, h k ≥ h min . Therefore, we have that (d k ) T h k ≥ (d k ) T h min , for all k = 1, . . . , n t , and As a result, the following holds: (B.6) or, equivalently, As a result, for all x = (z, q, h, η, θ) satisfying constraints (1b), (2), (3), (4), (5), (6) (see Sect. 2.1 in the manuscript) and z ∈ {0; 1} n CNV +n CNP , the denominator of I r (x) is strictly positive and 0 ≤ I r (x) ≤ 1. A value of I r (x) = 1 would correspond to there being no friction losses in the system, i.e. n t k=1 (d k ) T h k = n t k=1 (A 10 h k 0 ) T q k . Values of I r less than 0 on the other hand would mean that the pressure at the customer nodes is lower than the minimum requirement. If I r is exactly 0, the system can just about meet the minimum pressure requirements in the network.

Appendix C: Linear reformulation of −I r (Sect. 3.2.3)
The resilience index I r (x) given by (8) is a linear fractional function: both its numerator and denominator are linear (affine) functions of the problem variables. As a result, the relaxation of the potential-flow coupling constraint (2b) in MINLP −I r (Q) into (9) produces the mixed-integer linear fractional problem (10): minimize x= (z,q,h,η,θ) where F(Q) ⊂ {0, 1} n CNV +n CNP × R n p × R n n × R n p × R n p is the set of variables x defined by constraints (1b), (2a), (3), (4), (5) and (6), which depends on Q = {q k L , q k U } k=1,...,n t . In Charnes and Cooper (1962), the authors propose a linear reformulation for the continuous non-linear problem (NLP) of maximizing a linear fractional function subject to linear constraints. They show that, to solve the original problem, it is sufficient to solve two LPs (one, if the sign of the denominator of the linear fractional function is known). We provide an illustration of Charnes and Cooper's reformulation for the single variable function 3x−1 2x−3 , for x ∈ [2; 4], in Fig. 4. Unlike the problem studied in Charnes and Cooper (1962), the mixed-integer linear fractional problem (10) includes integer variables. Nonetheless, the particular form of the constraints in which the binary variables appear in the original problem (i.e. big-M constraints) allows us to apply a similar reformulation. We introduce a new variablet ∈ R * + , and multiply both the numerator and denominator of −I r (x) byt: By introducing the variablesq k =tq k ,h k =th k ,η k =tη k ,θ k =tθ k ,z = z, we can then reformulate (C.1) into: If we introduce the following constraint on the denominator of (C.2): then (C.2) becomes a linear function of variablesq k andh k : By multiplying byt, constraints (1b), (2a), (3) and (9) become: The combinatorial constraint (6) is not affected by the change of variables, while new upper and lower bounds on variablesq,h,θ ,η andt are given by: where t min = 1/M and t max = 1/m, with (M, m) ∈ (R * + ) 2 defined by Finally, we introduce the variablew ∈ R n CNV +n CNP and reformulate the big-M constraints on variablesq k ,θ k andη k , for all k = 1, . . . , n t (Belotti et al. 2013): and introduce the following constraint on variablesh k , for k = 1, . . . , n t : For solutions of the linear reformulation to be feasible solutions of (10), we must also ensure that and Az ≤ e. (C.11) Observe that (C.10) implies thatw =tz, whenz ∈ {0, 1} n CNV +n CNP . This produces the following MILP reformulation of (10): (z,q,h,θ,η,t,w) whereq = (q k ) k=1,...,n t ,h = (h k ) k=1,...,n t ,η = (η k ) k=1,...,n t ,θ = (θ k ) k=1,...,n t , F(Q) ⊂ {0; 1} (n CNV +n CNP ) × R n p × R n n × R n p × R n p × R × R n CNV +n CNP is the set defined by (C.3), (C.5a)-(C.5c), (C.6), (C.8), (C.10) and (C.11).
Before we show that a feasible solution of (10) can be produced from any feasible solution of (11) (Proposition 1), we establish the following result: for all x = (z, q, h, η, θ) feasible solutions of (10).

is a superset of the efficient set of flow variables of BOMINLP(Q).
Proof Let x = (z, q, h, η, θ) be an efficient solution of BOMINLP(Q init ). Then ) is a non-dominated solution of BOMINLP(Q init ) and, from Proposition 3.3., we have that f (x) must be in We then distinguish 3 cases: -there exists j = 1, . . . , n such that By definition of x f 1 , f 2 , k , we then have that , which contradicts x efficient. As a result, we must have that f 1 (x) ≤ UB f 1 , f 2 , k , for all k = 0, . . . , j − 1, and, in particular, . Moreover, as neither the branching nor the bound tightening steps in Algorithm 1 exclude feasible solutions during the search, there exists a node Q belonging to Q f 1 , f 2 , j , the tree of nodes explored during the solution of MINLP f 1 , f 2 , j (Q init ) such that q ∈ Q and x ∈ F(Q). As we have that L(Q) ≤ f 1 (x), then L(Q) ≤ min k=0,..., j−1 As a result, F(Q) ∈ F(Q f 2 ). Finally, we can conclude that, if x is an efficient solution of BOMINLP(Q init ), it belongs to a feasible set F(Q) defined by at least one of the explored nodes Q returned by Algorithm 1 when applied to the solution of MINLP f 1 , f 2 , j (Q init ), j = 1, . . . , n , or MINLP f i (Q init ), i = 1, 2, such that The notation UB * indicates that the solver terminates without finding a feasible solution

Appendix E: Solution of the anchor point problems and comparison with off-the-shelf optimization solvers (Sect. 4)
In this section, we evaluate the performance of the tailored sBB Algorithm 1 for solving the single-objective anchor point problems MINLP AZP (Q) and MINLP −I r (Q), compared to state-of-the-art global optimization solvers. We solve MINLP AZP (Q) and MINLP −I r (Q) for networks Net25 and pescara using the off-the-shelf solvers SCIP and BARON (with their default settings). Similarly to Algorithm 1, SCIP and BARON implement a spatial branch-and-bound approach and perform additional bound tightening on subproblems within the branch-and-bound procedure (Vigerske and Gleixner 2018;Tawarmalani and Sahinidis 2005). Unlike Algorithm 1, however, the off-the-shelf solvers produce lower bounds on the solution of the original problem by solving a convex continuous relaxation of the MINLP subproblems in the branch-andbound tree: non-convex constraints are relaxed using a convex or polyhedral envelope and integrality constraints are ignored. In addition, BARON implements advanced range reduction techniques, where inequalities based on the primal and dual optimal solutions of the convex lower bounding problems are added during the search to exclude feasible solutions with suboptimal objective value. We compare the performance of the different solution methods for a running time of 5 min (T max = 300s). In 3 out of the 4 experiments conducted, BARON was unable to produce a feasible solution within the time limit. As a result, we also present the results of the solution of MINLP −I r (Q) and MINLP AZP (Q) for networks Net25 and pescara using BARON, for a running time of 12 h (T max = 43200s). Finally, for all solvers, the relative optimality gap is defined as UB−LB min(|LB|,|UB|) . The results presented in Tables 3, 4, 5 and 6 show that, for all considered networks and problem instances, the spatial branch-and-bound Algorithm 1 consistently outperforms SCIP and BARON: in every experiment, Algorithm 1 returns a better feasible solution and optimality gap than both off-the-shelf solvers. This is consistent with previous experiments from the literature, where tailored algorithms are reported to outperform general purpose MINLP solvers for non-linear network optimization problems (Sahinidis 2019). In particular, previous literature on the optimal design and/or control of WDNs has shown that both SCIP and BARON might fail to pro- Net25 The notation UB * indicates that the solver terminates without finding a feasible  (Misener and Floudas 2013) and are outperformed by tailored spatial branch-and-bound algorithms ). Thanks to the linear reformulation of the resilience index −I r implemented in the lower bounding procedure, Algorithm 1 is able to produce a good quality feasible solution and tight lower bounds (compared to SCIP and BARON) for both instances of MINLP −I r (Q). In comparison, BARON, which was able to produce a lower bound and feasible solution for MINLP AZP (Q), both for Net25 and pescara, is unable to return a feasible solution for MINLP −I r (Q) or compute any useful lower and upper bounds in under 12 h (recall that, by definition, −1 ≤ −I r ≤ 0). As the problems MINLP AZP (Q) and MINLP −I r (Q) share the same constraints, these results show that the minimization of the resilience index −I r proved particularly challenging for BARON. SCIP, on the other hand, returns upper bounds close to the ones obtained with Algorithm 1 in 3 out of 4 instances. However, the final optimality gaps remain large. This might be because the lower bounding problems solved by SCIP drop the integrality constraints and the resulting relaxations are not as tight as the MILP relaxations derived in Algorithm 1.
Although the performance of the state-of-the-art solvers SCIP and BARON might improve significantly if non-default options are selected or different starting points are  Tables 3, 4, 5 and 6 support the implementation of Algorithm 1 for the solution of the scalarized problems MINLP f 1 , f 2 , j (Q) in Sect. 4.
Finally, we discuss the convergence of the spatial branch-and-bound Algorithm 1, and, in particular, the effect of the maximum computation time T max on the quality of the upper and lower bounds returned by Algorithm 1. We solve MINLP AZP (Q) and MINLP −I r (Q) for networks Net25 and pescara using Algorithm 1, this time with T max = 30 min. Figs. 5 and 6 represent the evolution of the upper and lower bounds for each problem as a function of CPU time/number of sBB iterations completed. In each figure, the red dotted cut-off line represents the upper and lower bounds returned by Algorithm 1 after 5 min, which correspond to the solutions presented in Tables 3  to 6. Figures 5 and 6 show that Algorithm 1 converges to a good feasible solution, or upper bound, after only a few sBB iterations and that the remaining iterations mainly consist in improving the lower bound. Similarly, lower bounds increase quickly at first, before the improvement in the optimality gap slows down. The results show that opti-mality gaps between the best solution found and the lower bound in MINLP AZP (Q) and MINLP −I r (Q) are more than halved after only 5 min of sBB, from 33% to 16% (i.e. 17% reduction in optimality gap) and 5.6% to 2.3% (i.e. 3.3% reduction in optimality gap) for Net25 (see Fig. 5). In comparison, running Algorithm 1 for 30 min only brings an additional 4.3% and 0.6% reduction respectively in the optimality gaps of MINLP AZP (Q) and MINLP −I r (Q). The same is observed for pescara in Fig. 6. Finally, the results presented in Figs. 5 and 6 are consistent with previous sBB performance reports from the literature on the optimal design-for-control of WDNs Ulusoy et al. 2019).
Although the quality of the upper and lower bounds returned by Algorithm 1 for the anchor point problems MINLP AZP (Q) and MINLP −I r (Q) depends on the computational time, the reduction in the optimality gap resulting from each additional sBB iteration decreases rapidly. In Figs. 5 and 6, in particular, there is little improvement in the optimality gaps past 5 min. Moreover, as problems MINLP AZP (Q) and MINLP −I r (Q) (but also MINLP f 1 , f 2 , j (Q)) are NP-hard, there is no guarantee that the optimality gap can be closed in a finite amount of time. As a result, we limit the computation time of Algorithm 1 to T max = 5 min for the solution of the scalarized problems MINLP f 1 , f 2 , j (Q), for networks Net25 and pescara, in Sect. 4. For a different case study network, however, it may be necessary to fix a different time limit in order to obtain good quality feasible solutions and optimality gaps. f i , η * f i , θ sim f i ) and θ sim f i = φ sim (q sim f i ). We perform a hydraulic simulation of the optimal network configurations obtained, for Net25 and pescara, by solving the anchor point problems MINLP AZP (Q) and MINLP −I r (Q) with Algorithm 1. In Figs. 7, 8, 9 and 10, we present the cumulative distribution functions F P and F Q corresponding to network configurations (z * AZP , η * AZP ) and (z * −I r , η * −I r ), defined for Net25 and pescara, as in Pecci (2018), by: The results in Figs. 7,8,9 and 10 show that the errors in network pressures and flows between the optimization and hydraulic simulation are small compared to the average pressure at network nodes and the total flow across network pipes in Net25 and pescara, respectively. For pescara, for instance, for both configurations (z * AZP , η * AZP ) and (z * −I r , η * −I r ), pressures are off by less than 0.4 m for more than 90% of network nodes, and flow errors do not exceed 0.3 l s −1 across over 90% of network pipes. In comparison, the total customer demand or cumulative flow across pipes in pescara is 563 ls −1 , whereas the minimum average zone pressure is 16.6 m-see Network Configuration  Table 8. Moreover, for both case studies networks Net25 and pescara, all simulated pressures for the resilience maximizing configurations (z * −I r , η * −I r ) satisfy minimum pressure requirements. For the AZP minimizing configurations (z * AZP , η * AZP ), more than 90% of all simulated pressures violate minimum pressure requirements by less than 0.04 m in Net25, while all simulated pressure violations are below 0.02 m in pescara. Finally, all simulated flows satisfy maximum flow bounds for each pipe, for both case studies networks Net25 and pescara.
Based on the small errors in the pressure and flow distributions between x * f i and x sim f i , we expect the quadratic approximation of the H-W friction head loss constraint to result in limited deviations in objective values f i (x * i ) compared to f i (x sim i ) for Net25 and pescara. Tables 7 and 8 present the optimal values returned by Algorithm 1 for MINLP AZP (Q) and MINLP −I r (Q), and the simulated AZP and −I r values corresponding to the optimal configurations (z * AZP , η * AZP ) and (z * −I r , η * −I r ). In particular, the results in Tables 7 and 8 show that, for both configurations (z * AZP , η * AZP ) and (z * −I r , η * −I r ), the deviations between the optimized and simulated objective values do not exceed 1.6%.
As shown by the results presented in Figs. 7, 8, 9 and 10, and Tables 7 and 8, the errors in network pressures and flows resulting from the quadratic approximation of the H-W friction head loss constraints in Net25 and pescara do not affect the quality of the solutions returned by Algorithm 1 for the optimal DfC problems  Table 9 Results of the solution of the sequence of parametrized problems MINLP AZP,−I r , j (Q) and MINLP −I r ,AZP, j (Q), j = 1, . . . , n , using Algorithm 1 and 2, for the network pescara and n = 20 In conclusion, the quality of the Pareto front approximation obtained for the biobjective optimal DfC problem BOMINLP(Q) of pescara by solving the sequence of scalarized problems MINLP f 1 , f 2 , j (Q) depends on the number of sBB iterations completed. As a single sBB iteration takes longer to complete for parametrized problems MINLP −I r ,AZP, j (Q) than for MINLP AZP,−I r , j (Q), solving the sequence of problems (MINLP AZP,−I r , j (Q)) j=1,...,n produces, in a given time, a tighter Pareto front approximation for BOMINLP(Q).